1.本发明涉及一种发动机气液两相流仿真方法，可以应用于工业仿真领域。
背景技术：2.两相流计算与单相流动计算的差异主要是需要计算包含体积分数的控制方程。流场中密度、黏度等参量会在相界面上产生阶跃，且相界面的位置和形状会在两相流动过程中发生变化，有时还可能发生破碎、聚合等现象。除此之外，两相间的质量、动量、能量传递也是通过相界面传递的，因此对相界面进行准确的追踪和描述是两相流数值计算的基础，也是求解其控制方程的前提。
3.vof方法是重要的界面追踪方法之一。其优点在于可以表示复杂界面的结构和变化，计算相对简单，相界面的锐利程度相对较高。但vof方法的流体体积函数在相界面处是不连续的，会导致流体参数的陡峭变化，即振荡被抹平，在计算相界面法向方向、曲率以及与曲率有关的物理量时，结果不准确。
技术实现要素：4.本发明要解决的技术问题是：vof方法的流体体积函数在相界面处是不连续的，会导致流体参数的陡峭变化，即振荡被抹平，在计算相界面法向方向、曲率以及与曲率有关的物理量时，结果不准确。
5.为了解决上述技术问题，本发明的技术方案是提供了一种航空发动机气液两相流数值模拟方法，在vof方法中引入level-set相函数其特征在于，包括以下步骤：
6.步骤1、完成vof相函数α以及level-set相函数的初始化；
7.步骤2、相界面重构：
8.采用plic方法实现相界面构造，其中，使用level-set相函数进行plic 方法中相界面法向向量n的求解；
9.网格中心到相界面的距离d则通过将两相网格中的复杂相界面构建成简单线段,进而通过点到线段距离的求解求出；
10.步骤3、level-set相函数重新距离化：
11.利用相界面重构过程中的结果，通过几何方法来实现level-set相函数的重新距离化，其过程如下：
12.对于两相网格，在完成相界面重构之后，得到网格单元中心到相界面的距离 d，距离d即为level-set相函值的绝对值，level-set相函数的正负可以由相标志函数求得：
[0013][0014]
式中，α表示体积分数，定义为液体相在所研究的局部控制容积中所占的体积份额；
[0015]
因此有：
[0016][0017]
经过上述过程可以得到各个网格单元的level-set相函数值
[0018]
步骤4、求解流输运方程，包括level-set相函数对流输运方程求解以及 vof相函数α对流输运方程求解。
[0019]
优选地，所述vof相函数α与所述level-set相函数的初始化结果相互使用，具体包括以下步骤：
[0020]
通过vof相函数α的分布将两相网格中的复杂相界面构建成简单线段，进而通过点到线段距离的求解完成level-set相函数的初始化。
[0021]
优选地，步骤2中，对于二维网格，所述相界面法向向量n求解方法如下：
[0022][0023][0024]
式中，n
x
(i,j)表示节点(i,j)的x方向的向量，ny(i,j)表示节点(i,j)的y方向的向量，δx表示x方向的单位量，δy表示y方向的单位量，表示节点(i,j)的控制体单元。
[0025]
优选地，步骤2中，对于三维空间，采用非结构网格梯度计算方法求解相界面法向向量n。
[0026]
优选地，步骤2中，本发明采用格林-高斯方法和最小二乘法进行非结构网格求解。
[0027]
优选地，采用格林-高斯方法计算网格单元i中心处的梯度其基于格心格式的公式如下：
[0028][0029]
式(1)中，是界面ij的法向量，vi是控制体体积，δs
ij
是界面面积，是网格单元i中心处的值，是与网格单元i相邻的网格单元j中心处的值；
[0030]
或者采用最小二乘法实现梯度计算，根据有限体积法，计算表达式如下：
[0031][0032]
式中，为i、j之间的偏移向量，则有：则有：表示网格单元j的偏移量，表示网格单元i的偏移量，网格单元j为网格单元i的相邻单元；
[0033]
将式(1)或式(2)推广到网格单元i的所有相邻单元，构造一个线性约束系统：
[0034]
[0035]
式(3)中，δx
ij
表示网格单元i与网格单元j的距离，表示控制体单元的梯度。
[0036]
求解式(3)所示的线性方程组可得：
[0037][0038]
式中，式中，表示控制体单元的梯度，w
ij
表示权重，由下式定义：
[0039][0040]
式中，xj表示小单元的整体j的均值；xi表示小单元的整体i的均值；xi表示网格单元i的x值，xj表示网格单元j的x值；
[0041]
优选地，步骤4中，csf模型将相界面处的表面张力等效成为作用于相界面附近一定宽度区域的体积力或称为质量力，在等效区域上，体积力是连续的，且在等效区域外的表面张力为零；表面张力作为源项添加到流输运方程中；
[0042]
跨过表面的压降依赖于通过两个半径的正交方向量度的表面曲率r1和r1和表面张力系数σ：
[0043][0044]
式中，p1表示第一节点的压强，p2表示第二节点的压强。
[0045]
表面张力表示为体积力后可以加入流输运方程的源项，形式如下：
[0046][0047]
式中，ρ表示密度，表示速度，表示受力，p表示压强，μ表示动力粘性系数，表示张量，表示取梯度，σ表示表面张力系数，κ(x)表示相界面的曲率，表示c(x)的外向量，上式中的除数，用以计算两相交界的区域宽度，c(x)表示染色函数：
[0048][0049]
式中，x表示表面上的点，c1表示第一相流体的系数，c2表示第二相流体的系数，流
体1为两相中的第一相，流体2为两相中的第二相。
[0050]
当曲率中心处于液相中时，曲率为正值；反之为负值；从相界面形状来看，曲率为正值时，相界面处于凸起；曲率为负值时，相界面处于凹陷；使用level-set 相函数对曲率进行计算，计算方法如下：
[0051][0052]
式中，表示相界面法向向量，具体实现方式如下：
[0053][0054]
本发明在vof方法中引入level-set函数，即clsvof方法，利用其距离函数追踪相界面时的连续性，改善vof方法在计算两相流动时，相界面曲率、法向向量及表面张力等参数不准确的问题，提高计算的精确度。
[0055]
本发明通过把相界面重构和距离函数φ的重新初始化使level-set方法和vof方法进行联合，提供了一种cls-vof方法。本发明提供的cls-vof方法实现气液两相界面捕捉功能，界面完整、精确性高。在本发明提供的cls-vof方法中，界面曲率、法向向量及表面张力等参数精确性高。
附图说明
[0056]
图1示意了体积分数求解流程。
具体实施方式
[0057]
本说明书中附图所显示的大小、比例等只是示意性的，用以配合说明书所描述的内容，并非用以限定本发明的实施条件，不影响本发明所产生的功效。本说明书中所述的“上”、“下”、“内”、“外”等位置关系仅是为了方便描述，而非用以限定本发明的可实施范围，其相对关系的改变，在无实质变更技术内容下亦视为本发明的可实施范畴。
[0058]
本发明提供的cls-vof两相流方法中存在vof相函数α和level-set相函数两种相函数：
[0059]
(1)vof方法
[0060]
vof相函数α又称体积分数，被定义为液体相在所研究的局部控制容积中所占的体积份额，其对流输运方程为：
[0061][0062]
对流输运方程反映了vof相函数α因流动而产生的变化，即vof相函数α随速度变化的规律。但是，vof相函数α是不连续函数，在相界面附近无法直接求导，此时vof模型的对流输运方程是无法求解的。需通过一定的方法对两相网格单元中的相界面形状进行构造。由于一般求解输运方程的数值方法容易由数值耗散引起模糊或振荡，而两相流计算是存在相界面的，为了使得界面位置和形状的精细程度达到目标，需要使用相应的计算格式求解输运方程。
[0063]
(2)level-set方法
[0064]
level-set方法采用隐式相界面捕捉方法，计算精度高，可以解决由于体积分数不连续引起的复杂相界面捕捉精度较低的问题。level-set相函数又称为距离函数或level-set函数，定义的是单元中心与相界面上任意一点连线所构成线段长度的最小值，且当单元中心位于目标相态区时level-set相函数取正值，即有：
[0065]
单元中心在目标相态内；
[0066]
单元中心在相界面上；
[0067]
单元中心在目标相态外
[0068]
为保证level-set相函数的零等值面为相界面，在任意时刻，对于相界面上的点，其level-set相函数值为零，即得到level-set相函数的对流输运方程为：
[0069][0070]
level-set相函数的对流输运方程体现的是相界面位移随时间的变化律。在使用level-set对流输运方程求解时，在计算中，level-set相函数可能不再具有距离函数的性质，若继续计算则结果将出现错误。通常情况下，在下一步计算前需要将level-set相函数重新初始化，但是在level-set相函数输运和重新初始化过程中，并不能严格遵守质量守恒定律。
[0071]
因此本发明中将vof方法与level-set方法进行结合，实现两种方法的优势互补，解决了vof方法中由于vof相函数α不连续所引起复杂相界面捕捉精度低的缺点。
[0072]
本发明求解体积因子方程的流程如图1所示：
[0073]
vof气液两相流计算模块拟采用本发明提供的cls-vof方法，主要包括以下步骤：
[0074]
步骤1、相函数初始化：
[0075]
相界面初始化是指在初始时刻，建立流体的相态分布与相函数值之间的对应关系。根据定义，vof方法和level-set方法分别是对vof相函数α及level-set 相函数进行初始化。根据本发明提供的cls-vof方法，两种相函数的初始化结果可以相互使用，具体流程如下：
[0076]
通过vof相函数α的分布将两相网格中的复杂相界面构建成简单线段，进而通过点到线段距离的求解完成level-set相函数的初始化。
[0077]
步骤2、相界面重构：
[0078]
气液相界面的构造是指在气液混合单元中将气相与液相的相界面表示出来，从而确定目标相态在流场中的分布形状。相界面构造是使用vof方法进行相界面追踪过程中不可或缺的关键环节。相界面构造的直接目的是确定两相网格中相界面的形状，为vof相函数α通量的计算及对流输运方程求解打好基础。
[0079]
本发明的数值模型中采用plic方法实现相界面构造。plic方法对应的相界面构造过程有两个关键参量：相界面法向向量n及网格中心到相界面的距离d。网格vof函数值及网格level-set函数值均要通过对流输运方程的求解给出。由于level-set相函数为连续函数，使用level-set进行相界面法向向量n求解时的计算精度将高于vof相函数α相应的求解精度。
[0080]
不失一般性，以二维网格为例，相界面法向向量n求解方法如下：
[0081][0082][0083]
式中，，n
x
(i,j)表示节点(i,j)的x方向的向量，ny(i,j)表示节点(i,j)的y方向的向量，δx表示x方向的单位量，δy表示y方向的单位量，表示节点(i,j)的控制体单元。
[0084]
对于三维空间，三元函数的梯度等同于三维等值面的法向量，因此，可以采用非结构网格梯度计算方法求解相界面法向向量n。本发明拟采用格林-高斯方法和最小二乘法进行计算。
[0085]
以三角形非结构网格为例，采用格林-高斯方法(green-gaussgradientmethod，以下简称gs法)计算网格单元i中心处的梯度其基于格心格式的公式如下：
[0086][0087]
式(1)中，是界面ij的法向量，vi是控制体体积，δs
ij
是界面面积，是网格单元i中心处的值，是与网格单元i相邻的网格单元j中心处的值。
[0088]
最小二乘法也可以实现梯度计算，根据有限体积法，计算表达式如下：
[0089][0090]
式中，为i、j之间的偏移向量，则有：表示网格单元j的偏移量，表示网格单元i的偏移量，网格单元j为网格单元i的相邻单元。
[0091]
将式(1)或式(2)推广到网格单元i的所有相邻单元，那么就可以构造一个线性约束系统：
[0092][0093]
式(3)中，δx
ij
表示网格单元i与网格单元j的距离，表示控制体单元的梯度。
[0094]
求解式(3)所示的线性方程组可得：
[0095][0096]
式中，式中，表示控制体单元的梯度，w
ij
表示权重，由下式定义：
[0097][0098]
式中，xj表示小单元的整体j的均值；xi表示小单元的整体i的均值；xi表示网格单元i的x值，xj表示网格单元j的x值；
[0099]
gs法和lsqr法是两种成熟的梯度计算方法，适用于非结构网格求解。
[0100]
网格中心到相界面的距离d则通过将两相网格中的复杂相界面构建成简单线段,进而通过点到线段距离的求解求出。
[0101]
由此可以实现结合level-set方法的高精度相界面重构。
[0102]
步骤3、level-set相函数重新距离化：
[0103]
由于在数值计算过程中无法避免的会出现数值耗散和舍入误差等情况，故 level-set相函数经过对流输运方程求解之后将不再具有距离函数的性质，这同时也是导致level-set方法计算过程中质量不能保持恒定的原因。故每经过一个时间步长计算后，需要对下一时间步的level-set相函数进行处理，保证在不改变其零等值面位置的前提下使其重新具有距离函数的性质，该处理过程称之为 level-set相函数的重新距离化。
[0104]
对于level-set方法，level-set相函数的重新距离化是利用一定的函数关系式通过迭代计算对其对流输运方程求解后的相函数值进行修正，这是一个复杂的代数求解过程。而对于本发明提出的cls-vof方法，由于距离函数的重新初始化的本质是重新计算网格单元中心与相界面上任意一点连线所构成的最短距离，可以利用相界面重构过程中的结果，通过几何方法来实现level-set相函数的重新距离化，其过程如下：
[0105]
对于两相网格，在完成相界面重构之后，网格单元中心到相界面的距离d(即 level-set相函值的绝对值)就已经得到。level-set相函数的正负可以由相标志函数求得：
[0106][0107]
式中，α表示体积分数，定义为液体相在所研究的局部控制容积中所占的体积份额；
[0108]
因此有：
[0109][0110]
经过上述过程可以得到各个网格单元的level-set相函数值
[0111]
步骤4、求解流输运方程，包括level-set相函数对流输运方程求解以及 vof相函数α对流输运方程求解。
[0112]
表面张力是气液两相流动过程中的重要参数，对气液两相的流动有着重要影响。单相流动中，分子受到各个方向上的力能相互抵消，合力为零。而在两相流动中，相界面处的分子由于界面两侧流体物理性质不同，导致其受到的各个方向上的合力不能相互抵消，合外力就以张力的形式表现出来，即两相流动过程中的表面张力，故表面张力产生的原因是两相界面处的分子引力不平衡。目前较为常用的表面张力计算模型有csf模型和pbm模型等。本发明以csf模型为例，介绍气液两相流计算中表面张力的一种计算方式，具体计算模型在详细设计和实施阶段会根据测试结果进行选取。
[0113]
csf模型将相界面处的表面张力等效成为作用于相界面附近一定宽度区域的体积力或称为质量力。在等效区域上，体积力是连续的，且在等效区域外的表面张力为零。表面张力作为源项添加到动量方程中。跨过表面的压降依赖于通过两个半径的正交方向量度的表面曲率r1和r1和表面张力系数σ：
[0114][0115]
式中，p1表示第一节点的压强，p2表示第二节点的压强，即跨过表面的第二节点。
[0116]
表面张力表示为体积力后可以加入流输运方程的源项，形式如下：
[0117][0118]
式中，ρ表示密度，表示速度，表示受力，p表示压强，μ表示动力粘性系数，表示张量，表示取梯度，σ表示表面张力系数，κ(x)表示相界面的曲率，表示c(x)的外向量，上式中的除数，用以计算两相交界的区域宽度，c(x)表示染色函数：
[0119][0120]
式中，x表示表面上的点，c1表示第一相流体的系数，c2表示第二相流体的系数，流体1为两相中的第一相，流体2为两相中的第二相。
[0121]
曲率是标量，有正负之分。当曲率中心处于液相中时，曲率为正值；反之为负值。从相界面形状来看，曲率为正值时，相界面处于凸起；曲率为负值时，相界面处于凹陷。曲率计算是表面张力计算中的关键环节，其计算的准确性直接影响到表面张力的计算精度。由于level-set相函数为连续函数，可进行直接的求导运算，用于曲率计算时无需进行特殊的转化和处理即具有较高的计算精度，拟使用该函数对曲率进行计算，计算方法如下：
[0122][0123]
式中，表示相界面法向向量，具体实现方式如下：
[0124]
[0125]
需要说明的是，曲率是标量，有正负之分。当曲率中心在液相中时，曲率取正；反之，取负。从相界面形状来看即是，相界面凸处，曲率为正；相界面凹处，曲率为负。
技术特征：1.一种航空发动机气液两相流数值模拟方法，在vof方法中引入level-set相函数其特征在于，包括以下步骤：步骤1、完成vof相函数α以及level-set相函数的初始化；步骤2、相界面重构：采用plic方法实现相界面构造，其中，使用level-set相函数进行plic方法中相界面法向向量n的求解；网格中心到相界面的距离d则通过将两相网格中的复杂相界面构建成简单线段，进而通过点到线段距离的求解求出；步骤3、level-set相函数重新距离化：利用相界面重构过程中的结果，通过几何方法来实现level-set相函数的重新距离化，其过程如下：对于两相网格，在完成相界面重构之后，得到网格单元中心到相界面的距离d，距离d即为level-set相函值的绝对值，level-set相函数的正负可以由相标志函数求得：式中，α表示体积分数，定义为液体相在所研究的局部控制容积中所占的体积份额；因此有：经过上述过程可以得到各个网格单元的level-set相函数值步骤4、求解流输运方程，包括level-set相函数对流输运方程求解以及vof相函数α对流输运方程求解。2.如权利要求1所述的一种航空发动机气液两相流数值模拟方法，其特征在于，所述vof相函数α与所述level-set相函数的初始化结果相互使用，具体包括以下步骤：通过vof相函数α的分布将两相网格中的复杂相界面构建成简单线段，进而通过点到线段距离的求解完成level-set相函数的初始化。3.如权利要求1所述的一种航空发动机气液两相流数值模拟方法，其特征在于，步骤2中，对于二维网格，所述相界面法向向量n求解方法如下：所述相界面法向向量n求解方法如下：式中，n
x
(i,j)表示节点(i,j)的y方向的向量，n
x
(i,j)表示节点(i,j)的y方向的向量，δx表示x方向的单位量，δy表示y方向的单位量，表示节点(i,j)的控制体单元。4.如权利要求1所述的一种航空发动机气液两相流数值模拟方法，其特征在于，步骤2中，对于三维空间，采用非结构网格梯度计算方法求解相界面法向向量n。5.如权利要求4所述的一种航空发动机气液两相流数值模拟方法，其特征在于，步骤2中，本发明采用格林-高斯方法和最小二乘法进行非结构网格求解。
6.如权利要求4所述的一种航空发动机气液两相流数值模拟方法，其特征在于，采用格林-高斯方法计算网格单元i中心处的梯度其基于格心格式的公式如下：式(1)中，是界面ij的法向量，v
i
是控制体体积，δs
ij
是界面面积，是网格单元i中心处的值，是与网格单元i相邻的网格单元j中心处的值；或者采用最小二乘法实现梯度计算，根据有限体积法，计算表达式如下：式中，为i、j之间的偏移向量，则有：则有：表示网格单元j的偏移量，表示网格单元i的偏移量，网格单元j为网格单元i的相邻单元；将式(1)或式(2)推广到网格单元i的所有相邻单元，构造一个线性约束系统：式(3)中，δx
ij
表示网格单元i与网格单元j的距离，表示控制体单元的梯度；求解式(3)所示的线性方程组可得：式中，式中，表示控制体单元的梯度，w
ji
表示权重，由下式定义：式中，x
j
表示小单元的整体j的均值；x
i
表示小单元的整体i的均值；x
i
表示网格单元i的x值，x
j
表示网格单元j的x值；7.如权利要求1所述的一种航空发动机气液两相流数值模拟方法，其特征在于，步骤4中，csf模型将相界面处的表面张力等效成为作用于相界面附近一定宽度区域的体积力或称为质量力，在等效区域上，体积力是连续的，且在等效区域外的表面张力为零；表面张力
作为源项添加到流输运方程中；跨过表面的压降依赖于通过两个半径的正交方向量度的表面曲率r1和r1和表面张力系数σ：式中，p1表示第一节点的压强，p2表示第二节点的压强；表面张力表示为体积力后可以加入流输运方程的源项，形式如下：式中，ρ表示密度，表示速度，表示受力，p表示压强，μ表示动力粘性系数，表示张量，表示取梯度，σ表示表面张力系数，κ(x)表示相界面的曲率，表示c(x)的外向量，上式中的除数，用以计算两相交界的区域宽度，c(x)表示染色函数：式中，x表示表面上的点，c1表示第一相流体的系数，c2表示第二相流体的系数，流体1为两相中的第一相，流体2为两相中的第二相；当曲率中心处于液相中时，曲率为正值；反之为负值；从相界面形状来看，曲率为正值时，相界面处于凸起；曲率为负值时，相界面处于凹陷；使用level-set相函数对曲率进行计算，计算方法如下：式中，表示相界面法向向量，具体实现方式如下：
技术总结本发明提供了一种航空发动机气液两相流数值模拟方法，在VOF方法中引入Level-Set相函数其特征在于，包括以下步骤：步骤1、完成VOF相函数α以及Level-Set相函数的初始化；步骤2、相界面重构：步骤3、Level-Set相函数重新距离化：步骤4、求解流输运方程。本发明通过把相界面重构和距离函数φ的重新初始化使Level-Set方法和VOF方法进行联合，提供了一种CLS-VOF方法。本发明提供的CLS-VOF方法实现气液两相界面捕捉功能，界面完整、精确性高。在本发明提供的CLS-VOF方法中，界面曲率、法向向量及表面张力等参数精确性高。及表面张力等参数精确性高。及表面张力等参数精确性高。
技术研发人员：陈灏 詹翼
受保护的技术使用者：上海索辰信息科技股份有限公司
技术研发日：2022.02.11
技术公布日：2022/7/5
