一种基于Liutex-Omega涡识别理论的空化流动数值模拟方法与流程

本发明属于空化流动数值模拟技术领域，具体涉及一种基于liutex-omega涡识别理论的空化流动数值模拟方法。

背景技术：

空化是一种发生在液相中特有的相变现象，在一定温度下，当水流中的压力降低至其饱和蒸汽压时将产生充满气体和蒸汽的空泡，即水流中发生空化，尤其在旋涡流动的核心区极易诱发涡空化。涡空化是一种常见却又复杂的空化形式，涡空化流动中旋涡与空化之间相互影响，涡量集中区域形成低压涡心引发空化，与此同时，空化又成为产生涡的影响机制。涡空化的出现通常是不利的，会严重影响水泵、水轮机以及螺旋桨等的工作性能及运行稳定性。涡空化流动的研究有助于认知其流动特性，对于解决涡空化引发的各种工程问题具有重要的意义。

目前，基于计算流体动力学(cfd)的数值模拟技术发展迅速，成为涡空化发生发展机理及其流动特性的重要研究手段，而有效的空化流动数值模拟方法是准确预测涡空化流动的前提与保证。空化是包含气液相间质量传输的非定常可压缩多相湍流流动现象，空化流动数值模拟的关键在于合适的湍流模型和空化模型。关于湍流模型，涡空化中旋涡会引起强烈的旋转曲率效应，而通过对传统涡粘模型进行适当修正可以反映出旋转曲率效应的作用影响。但是关于空化模型方面，普遍采用的均相流模型无法反映出旋涡运动对空化相变过程中质量传输过程的影响，从而导致涡空化流动数值精度不高，所预测的涡空化流态与真实流态偏差较大，无法反映出涡空化流动特点。因此，有必要建立一种有效反映涡空化流动特点的高精度空化流动数值模拟方法。

技术实现要素：

本发明的目的是针对现有针对涡空化流动数值模拟技术的不足，提出一种基于liutex-omega涡识别理论的空化流动数值模拟方法。

为实现本发明的目的，采用如下技术方案：

步骤一：建立计算域的三维几何模型，根据所研究的三维空化流动计算域的结构尺寸图，利用几何建模软件croe、autocad、solidworks、ug或catia等建立计算域的三维几何模型，并按照计算网格划分软件的文件格式要求输出相应的几何文件；

步骤二：计算域网格划分，将步骤一输出的几何文件导入计算网格划分软件ansysicemcfd、ansysgambit、hypermesh或truegrid中，采用高质量的六面体结构化网格方案划分计算网格，对近壁区、叶顶间隙区域以及涡空化流动范围区域的网格进行加密处理以便于数值计算获得精细化的速度场、压力及其脉动、湍动能、湍流耗散率、涡量、汽、液相的体积分数以及标准化的liutex-omega涡识别系数ωl等各种流动物理参数，采用基于理查德森外推法的gci(gridconvergenceindex)、gci-or、gci-ln、gci-r、cf(correctionfactormethod)、fs(factorandsafety)或fs1等网格不确定度估计方法对计算网格离散误差进行估计并确定计算网格的数量，生成后的计算网格按照步骤四的数值模拟计算软件所需的计算网格文件格式导出相应的网格文件；

步骤三：建立三维空化流动计算流体力学模型，三维空化流动的计算流体力学模型包括流动控制方程、湍流模型以及基于liutex-omega涡识别系数ωl分区的空化模型；

步骤三所述的流动控制方程由连续性方程、动量方程和组分输运方程构成：

连续性方程：

动量方程：

组分输运方程：

其中，ρ和μ为汽液混合相的密度和动力粘性系数且满足：

ρ＝ρvαv+ρlαl(4)

μ＝μvαv+μlαl(5)

公式中t为时间(单位s)，xi、xj为网格节点坐标且i，j＝1,2,3，δij为克罗内克符号；μt为湍流粘性系数(单位n·s/m²)，ui、uj为汽液混合相速度分量(单位m/s)，uvi为气相速度分量(单位m/s)，p为流场压力(单位pa)，ρv和ρl分别为汽相和液相的密度(单位kg/m³)，μv和μl分别为汽相的动力粘性系数(单位n·s/m²)和液相的动力粘性系数(单位n·s/m²)，αv和αl分别为汽相和液相体积分数且满足：

αv+αl＝1(6)

公式(3)中的源项和分别表示汽、液两相间的蒸发率和凝结率；

湍流粘性系数μt通过湍流模型来求解，考虑涡空化中旋涡流动具有强烈的旋转曲率效应，采用旋转曲率修正的湍流模型：

公式中，k为湍动能，ω为比耗散率；

湍动能生成项pk为：

湍流粘性系数μt求解公式为：

经验常数β^*＝0.09，a1＝0.31；对于湍流模型方程中各变量系数σk、σω、βk、βω以及γ采用wilcoxk-ω模型方程和标准k-ε模型方程的相关经验常数加权得到：

φ＝f1φ1+(1-f1)φ2，φ＝{σk,σω,βk,βω,γ}(11)

其中，wilcoxk-ω模型方程的各常数为：σk1＝0.85，σω1＝0.5，βk1＝0.09，βω1＝0.075，γ1＝5/9；标准k-ε模型方程的各常数为：σk2＝1.0，σω2＝0.856，βk2＝0.09，βω2＝0.0828，γ2＝0.44；

对于混合函数f1和f2:

公式(7)和(8)中的旋转曲率修正系数fr为：

其中cscale＝1.0，参数的表达式为：

上式中的各相关变量的表达式分别为：

r^*＝s/ω(20)

d²＝max(s²,0.09ω²)(25)

其中，sij为应变率张量，ωij为旋转率张量，εimn和εmji为勒维-契维塔符号，各相关变量参考的坐标系是惯性坐标系或旋转坐标系且旋转速度为ωm，模型常数cr1＝1.0，cr2＝2.0，cr3＝1.0；

空化汽液相变过程模拟基于简化rayleigh-plesset方程的空化模型，兼顾考虑流场局部压力脉动对空化影响，对于蒸发率和凝结率分别为：

其中，pv为液体饱和蒸汽压强，空泡半径rb＝1μm，成核区体积分数αb＝0.0005，蒸发系数fv＝50，凝结系数fc＝0.01；

采用基于liutex-omega涡识别理论方法对三维空化流场中的旋涡区进行判断，liutex-omega涡识别系数ωl为：

其中，bl＝0.001～0.002，且对于αl和βl通过下列公式求得：

式(29)、(30)中的各速度梯度是基于原始坐标系下速度梯度张量的坐标变换获得：

其中λl为的实特征值，r为旋转实特征向量；

根据三维空化流场中涡识别系数ωl分布情况对强涡区域进行辨识，(27)式中建立基于ωl分区的修正系数fωl对空化凝结率进行修正，以体现涡空化流动中旋涡对汽-液凝结过程的延迟作用机制，fωl的表达式为：

步骤四：根据所研究对象的实际运行工况，设定相应的边界条件进行三维空化流场的数值计算，计算模型的进口液相体积分数为1、湍流度采用低湍流度1％，固壁采用无滑移边界，基于步骤三建立的空化流动计算流体动力学模型，采用cfd计算软件ansyscfx、ansysfluent、openform、star-cd、star-ccm、numeca或pumplinx对三维空化流场进行计算；

步骤五：对步骤四的计算结果进行后处理，采用cfd软件自带后处理及通用cfd后处理软件，如cfx-post、tecplot、paraview、ensight、fieldview、cfview等对计算输出的结果进行分析处理，获得速度场、压力差、涡量场、湍动能、汽相和液相体积分数，以及涡空化形态及其演化发展过程等空化流场信息。

本发明的有益效果是：

本发明的一种基于liutex-omega涡识别理论的空化流动数值模拟方法，相比现有技术，兼顾考虑了涡空化流动中旋涡流动引起的旋转曲率效应以及旋涡与空化之间的相互作用影响，通过引入标准化的liutex-omega涡识别系数ωl对空化流场中的强涡区进行有效识别，并采用基于ωl分区的空化流动计算流体力学模型进行数值计算，实现了对涡空化流动区域快速高精度的数值预测，有效反映出涡空化的流动特点及其时空演化规律，显著提高了涡空化流动的预测精度，可推广用于三维水翼、水泵、水轮机、螺旋桨等具有显著涡空化现象的三维空化流动数值模拟研究，对认识涡空化流动特性及辅助解决涡空化流动所引起的诸多问题具有重要的现实意义和工程应用价值。

附图说明

图1是本发明基于liutex-omega涡识别理论的空化流动数值模拟方法的计算流程图；

图2a是本发明实施例的三维水翼模型结构示意图，2b是断面上的旋涡涡心位置对比情况；

图3是本发明实施例的三维水翼的计算网格图；

图4是本发明实施例三维水翼典型断面涡心位置对比图；其中，4a示出了三维水翼中心点o沿主流z向的0.1c断面上的旋涡涡心位置对比情况，4b示出了三维水翼中心点o沿主流z向的0.12c断面上的旋涡涡心位置对比情况；

图5是本发明实施例的三维水翼涡空化流态对比图，其中，5a为实验观测到的空化流态图，5b为现有空化数值模拟方法所预测的空化流态图，5c为本发明空化数值模拟方法所预测的空化流态图。

具体实施方式

下面结合附图与实施例对本发明作进一步详细描述。

实施例

本实施例是以带有叶顶间隙的三维水翼为对象进行相关说明，一种基于liutex-omega涡识别理论的空化流动数值模拟方法，如图1所示，主要包括以下步骤：

步骤一：带有叶顶间隙的三维水翼的计算域为矩形水槽，如图2a、2b，三维水翼截面为naca0009翼型，翼型弦长c＝0.1m，水槽长度为0.8m、宽度为0.15m，水翼的一侧与水槽边壁相接、另一侧则与其邻近边壁具有wtip＝1.98mm的间隙，利用几何建模软件croe建立计算域的三维几何模型，并输出*.igs格式文件用于步骤二的计算网格划分；

步骤二：计算域网格划分，将步骤一输出的几何文件*.igs导入网格划分软件ansysicemcfd中，采用高质量的六面体结构化网格方案划分计算网格，对水翼壁面区、叶顶间隙区域以及涡空化流动范围区域的网格进行加密处理以便于数值计算获得精细化的速度场、压力及其脉动、湍动能、湍流耗散率、涡量、汽、液相的体积分数以及标准化的liutex-omega涡识别系数ωl等各种流动物理参数，通过采用基于理查德森外推法的gci(gridconvergenceindex)网格不确定度估计方法对计算网格离散误差进行估计，确定计算网格的数量为480万，如图3所示，导出计算网格文件*.cfx5用于步骤四的cfd计算；六面体结构化网格方案为现有技术，在此不再赘述；

步骤三：建立三维空化流动计算流体力学模型，三维空化流动的计算流体力学模型包括流动控制方程、湍流模型以及基于liutex-omega涡识别系数ωl分区的空化模型；

步骤三所述的流动控制方程由连续性方程、动量方程和组分输运方程构成：

连续性方程：

动量方程：

组分输运方程：

其中，ρ和μ为汽液混合相的密度和动力粘性系数且满足：

ρ＝ρvαv+ρlαl(4)

μ＝μvαv+μlαl(5)

公式中t为时间(单位s)，xi、xj为网格节点坐标且i，j＝1,2,3，δij为克罗内克符号；μt为湍流粘性系数(单位n·s/m²)，ui、uj为汽液混合相速度分量(单位m/s)，uvi为气相速度分量(单位m/s)，p为流场压力(单位pa)，ρv和ρl分别为汽相和液相的密度(单位kg/m³)，μv和μl分别为汽相的动力粘性系数(单位n·s/m²)和液相的动力粘性系数(单位n·s/m²)，αv和αl分别为汽相和液相体积分数且满足：

αv+αl＝1(6)

公式(3)中的源项和分别表示汽、液两相间的蒸发率和凝结率；

湍流粘性系数μt通过湍流模型来求解，考虑涡空化中旋涡流动具有强烈的旋转曲率效应，采用旋转曲率修正的湍流模型：

公式中，k为湍动能，ω为比耗散率；

湍动能生成项pk为：

湍流粘性系数μt求解公式为：

经验常数β^*＝0.09，a1＝0.31；对于湍流模型方程中各变量系数σk、σω、βk、βω以及γ采用wilcoxk-ω模型方程和标准k-ε模型方程的相关经验常数加权得到：

φ＝f1φ1+(1-f1)φ2，φ＝{σk,σω,βk,βω,γ}(11)

其中，wilcoxk-ω模型方程的各常数为：σk1＝0.85，σω1＝0.5，βk1＝0.09，βω1＝0.075，γ1＝5/9；标准k-ε模型方程的各常数为：σk2＝1.0，σω2＝0.856，βk2＝0.09，βω2＝0.0828，γ2＝0.44；

对于混合函数f1和f2:

公式(7)和(8)中的旋转曲率修正系数fr为：

其中cscale＝1.0，参数的表达式为：

上式中的各相关变量的表达式分别为：

r^*＝s/ω(20)

d²＝max(s²,0.09ω²)(25)

其中，sij为应变率张量，ωij为旋转率张量，εimn和εmji为勒维-契维塔符号，各相关变量参考的坐标系是惯性坐标系或旋转坐标系且旋转速度为ωm，模型常数cr1＝1.0，cr2＝2.0，cr3＝1.0；

空化汽液相变过程模拟基于简化rayleigh-plesset方程的空化模型，兼顾考虑流场局部压力脉动对空化影响，对于蒸发率和凝结率分别为：

其中，pv为液体饱和蒸汽压强，空泡半径rb＝1μm，成核区体积分数αb＝0.0005，蒸发系数fv＝50，凝结系数fc＝0.01；

采用基于liutex-omega涡识别理论方法对三维空化流场中的旋涡区进行判断，liutex-omega涡识别系数ωl为：

其中，bl＝0.001～0.002，且对于αl和βl通过下列公式求得：

式(29)、(30)中的各速度梯度是基于原始坐标系下速度梯度张量的坐标变换获得：

其中λl为的实特征值，r为旋转实特征向量；

根据三维空化流场中涡识别系数ωl分布情况对强涡区域进行辨识，(27)式中建立基于ωl分区的修正系数fωl对空化凝结率进行修正，以体现涡空化流动中旋涡对汽-液凝结过程的延迟作用机制，fωl的表达式为：

步骤四：基于步骤三建立的空化流动计算流体动力学模型，采用cfd计算软件ansyscfx对三维水翼空化流场进行计算，将步骤二输出的计算网格文件*.cfx5导入至ansyscfx软件中，根据实验工况，对于三维水翼计算模型的进口边界进行设定，设定均匀流速为10m/s、液相体积分数为1、采用1％低湍流度，固壁为无滑移边界，出口边界设定为压力边界条件且确保进口的压力值为10⁵pa；

步骤五：对步骤四的计算结果进行后处理，采用ansyscfx后处理软件cfx-post对计算输出的结果进行分析处理，获得速度场、压力差、涡量场、湍动能、汽相和液相体积分数，以及涡空化形态及其演化发展过程等空化流场信息。

图4和图5分别给出本实施例三维水翼实验测量与采用现有数值模拟方法和本发明数值模拟方法预测获得的典型断面涡心位置与空化流态。

图4a、4b为三维水翼中心点o沿主流方向的z/c＝1.0和z/c＝1.2断面上旋涡的涡心位置对比情况，图中x、y分别为两断面上的横、纵坐标。由图4可以发现，相比现有的空化数值模拟方法(如图4中的黑色三角形)，本发明所提出的数值模拟方法(如图4中的黑色圆形)与试验测量获得的断面上泄漏涡的涡心位置(如图4中的黑色方形)更为接近，表明本发明所提出的数值模拟方法具有较好的计算精度，可以更加精确的预测旋涡流动结构和流动特性。

图5为三维水翼叶顶间隙处的空化流态图，其中数值计算获得的空化流态采用90％的汽相体积分数等值面来反映。由图5a、5b可知，相比实验观测到的空化流态，现有的空化数值模拟方法所预测的无论是叶顶间隙空化还是泄漏涡空化，空化体积(图中阴影部分)明显偏小，特别是泄漏涡空化溃灭速度较快而表现为长度偏短；通过本发明提出的方法模拟预测的叶顶间隙空化和泄漏涡空化的流态与实验观测结果吻合良好，尤其是泄漏涡空化流态表现出在强烈的旋涡作用下空泡团在向下游运动的过程中，是随着旋涡强调减弱而空化现象才逐渐消失，如图5c所示。本发明的方法能够反映出由于旋涡旋转作用影响下的相间质量传输过程，可以更好地预测涡空化流动现象。

以上结合附图对本发明的技术方案做出详细说明，但本发明不局限于所描述的技术方案。对本领域的普通技术人员而言，在本发明的原理和技术思想的范围内，对这些实施方式进行多种变化、修改、替换和改进仍落入本发明的保护范围之内。