本发明属于流体力学技术领域，具体涉及一种多相流计算方法。
背景技术：
多相流是一种复杂的流动现象,广泛存在于石油、化工、冶金、动力以及能源等工业领域,多相流的研究对于相关行业科技的进步具有十分重要的意义。而在多相流的数值模拟中对于相界面的捕捉是其中的重点和难点。在传统的欧拉方法中为了解决相界面的捕捉问题，开发了很多的相界面追踪和重构方法如vof、levelset方法、front-tracking和voset方法等。然而这些方法重构或者构造相界面的过程中会丢失相界面的细节，并且需要复杂的网格划分操作，在处理多相流的计算过程中有很多的不便。在多相流的数值模拟中，基于拉格朗日框架的拉格朗日的方法显示出了其特有的优点。这些方法有sph方法，mps方法，以及fvp等方法。这些拉格朗日的方法被应用在计算气泡上升，溃坝，蒸汽爆炸，以及固液耦合运动等相界面变化剧烈的多相流的计算中。虽然这一类的方法能够自动的追踪相界面，不需要额外的努力来捕获或跟踪界面。但是这些方法的计算粒子数目太多，十分的消耗计算时间以及计算性能。
技术实现要素：
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种多相流计算方法，不仅能够对于多相流的相界面进行准确的追踪，而且尽可能的降低计算性能的消耗缩短计算需要的时间。
本发明采用以下技术方案：
一种多相流计算方法，首先确定计算区域，输入物性参数及相应的边界条件；使用约束插值守恒的半拉格朗日方法求解欧拉坐标系下控制方程的对流部分，求解欧拉坐标系下非对流部分的压力泊松方程，更新网格部分下个时间步的数值；将网格部分的速度传递给粒子计算部分进行修正，得到新的粒子位置，根据粒子部分计算结果更新界面附近网格的颜色函数值，达到模拟时间后完成多相流数值模拟计算。
具体的，对流部分中粘性流体的控制方程为：
其中，ρ是密度，
进一步的，对流项方程组为：
使用约束插值守恒的半拉格朗日方法求解得到下一个时间步的
具体的，求解欧拉坐标系下非对流部分的压力泊松方程具体为：
网格部分粘性流体控制方程的半离散非对流项方程为：
非对流部分计算如下：
压力的泊松方程计算如下：
其中，pn 1为下一个时刻的压力值，u□为对流部分计算得到的临时速度，ρn 1为下一个时刻的网格密度，m□是对流部分计算得到的临时动量。
进一步的，更新下一个时刻的网格的动量面积积分平均和体积积分平均具体为：
下一个时刻网格的动量面积积分平均为：
其中，
动量的下一个时间步的网格体积积分平均为：
其中，
具体的，粒子计算部分修正具体为：
粒子部分的粒子计算中，对于粘性不可压缩流体，其质量方程和动量守恒方程如下：
其中，ρ为流体密度，t为时间，
在二维有限体积法中，一个移动粒子的控制体是一个圆，其离散区域计算如下：
v＝πr2＝l2
其中，v为离散的体积，r粒子的控制体半径，l为网格的大小。
进一步的，根据高斯定理中曲面积分与体积分的关系，在二维有限体积粒子法中，梯度算子和拉普拉斯算子计算如下：
其中，φs为粒子i和粒子j之间作用的物理量，
其中，n0是初始粒子数密度，
初始的核函数积分值n0计算如下：
两粒子之间的单位向量
核函数wij计算如下：
其中，re为截断半径，当其中
梯度算子和拉普拉斯算子计算如下：
其中，φs可以表示为一个线性函数计算如下：
更进一步的，粒子部分的粒子初始速度计算如下：使用核函数从网格部分计算得到的网格的中心速度插值得到粒子的临时速度和位置计算如下：
其中，
其中，re为截断半径。
进一步的，粒子部分的粒子根据有限体积粒子法对于速度和位置的修正计算，对于不可压缩的流体，速度的散度计算如下：
根据粒子运动的动量方程确定速度的修正值与压力的梯度的关系如下：
压力泊松方程计算：
利用以上的梯度算子和拉普拉斯算子离散求解上述的压力泊松方程得到压力，修正的粒子速度与粒子压力的关系如下：
粒子的新的时刻的位置计算如下：
具体的，网格计算中的界面附近网格颜色函数值根据网格中粒子数密度计算，粒子部分计算完毕，根据粒子的位置更新界面处的网格的颜色函数的值，网格的颜色函数值θ为：
其中，n0是粒子按照4×4粒子布满网格时的粒子数密度，re为截断半径，其值为8.1倍的粒子初始距离。
与现有技术相比，本发明至少具有以下有益效果：
本发明提供一种计算多相流的混合方法，基于欧拉方法中的内部插值/多矩有限体积法和拉格朗日方法中的有限体积粒子法，结合多相流计算中对于相界面追踪的需要，把拉格朗日粒子布置在相界面附近进行计算，从而能够准确的追踪相界面，并且粒子只是布置在相界面附近，减少了计算需要的粒子数，缩短了计算的时间以及性能的消耗。
进一步的，结合欧拉方法中的多相流计算相界面捕捉的困难，以及相界面附近网格调整的不便，避免了欧拉方法中对于相界面的捕捉。
进一步的，拉格朗日的粒子只是布置在相界面附近，减少了拉格朗日方法计算需要的粒子数，从而缩短了计算的时间。
进一步的，通过把拉格朗日方法和欧拉方法结合起来，充分的结合了两者的优点，不仅充分利用了拉格朗日方法相界面捕捉准确的优点，而且利用了欧拉方法计算快速的优点。
综上所述，本发明结合了网格方法和拉格朗日方法的优点，通过计算粒子的运动来跟踪相界面的移动，无需特定的界面跟踪算法。提供了一种快速准确的多相流数值模拟计算方法。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明流程图；
图2为网格体积积分平均和面积积分的定义；
图3为有限体积粒子法的离散区域；
图4为网格与粒子速度传递的示意图；
图5为溃坝的计算区域划分图；
图6为混合方法中溃坝的模拟计算结果；
图7为距离时间计算结果；
图8为高度时间计算结果；
图9为单个气泡上升中的结算的区域；
图10为本方法对于单个气泡上升的计算结果。
具体实施方式
本发明提供了一种多相流计算方法，基于在欧拉方法中的内部插值/多矩有限体积法以及拉格朗日方法中的有限体积粒子法，结合实际数值模拟计算中欧拉方法对于多相流相界面捕捉的困难以及拉格朗日方法计算的耗时，通过把拉格朗日的粒子只是布置在相界面附近的几层网格中显著的减少了拉格朗日方法中的粒子数，并且避免了复杂的欧拉方法中的相界面的捕捉。混合方法结合了网格方法和拉格朗日方法的优点，其通过移动粒子法计算粒子的运动来跟踪相界面的移动，所以相界面由粒子法计算的粒子的分布自动确定，无需特定的前向跟踪算法，其余的部分通过网格方法计算，不会产生数值扩散，并且大大的缩短计算的时间。
请参阅图1，本发明一种多相流计算方法，使用网格法和移动粒子法混合的方法来进行多相流的计算。混合方法中网格部分使用内部约束插值/多矩的有限体积法进行求解，粒子部分使用有限体积粒子法进行求解。其中，混合方法的网格部分使用的内部约束插值/多矩有限体积法把粘性流体的控制方程分为对流项和非对流项进行分布求解，并且把网格的中心速度传递给混合方法的粒子计算部分，之后再进行粒子部分的计算对于粒子位置进行修正，修正后再根据界面附近的网格内的粒子数密度得到网格的颜色函数的值。混合方法中粒子只是布置在多相流的相界面附近，这样可以减少计算需要的拉格朗日粒子数从而缩短计算的时间，同时又可以利用拉格朗日方法追踪界面的优势而不需要任何特殊的方法来追踪界面。在数值模拟计算中大部分计算区域使用网格法计算，只是在多相流的相界面附近使用移动粒子法计算，拉格朗日部分的粒子计算完毕后根据相界面处的粒子数密度修正网格法中的颜色函数的值，从而实现快速又准确的多相流的计算，包括以下步骤：
s1、确定计算区域，以及计算区域相应的物性和边界条件
确定计算区域，以及计算区域中不同相的区域边界。在相应的网格部分的计算区域中给定物性参数，其中包括：密度、粘度、表面张力系数。根据网格部分计算区域的划分在相界面附近布置拉格朗日粒子。在拉格朗日粒子布置完毕之后根据粒子的初始布置计算得到各个网格的初始颜色函数的值。
s2、求解欧拉坐标系下控制方程的对流部分计算；
混合方法中首先使用欧拉方法中的内部插值/多矩的有限体积法计算全部的计算区域，再使用拉格朗日方法中的有限体积粒子法对于相界面附近的网格区域进行计算，根据粒子法计算得到的粒子布置更新网格法中相界面附近网格的颜色函数值，而拉格朗朗日的粒子又清楚的描述着相界面的边界。
混合方法中欧拉部分的内部插值/多矩有限体积法，把粘性流的控制方程分为两个部分，一部分是对流部分，另外一部分是非对流部分。对于对流部分的计算如下：
网格方法中粘性流体的控制方程如下：
其中，ρ是密度，
其中，对流项方程组为：
其中，ρ是密度，
对流方程使用约束插值守恒的半拉格朗日方法求解得到下一个时间步的
进一步的，在对流的计算后再计算粘滞力、表面张力以及重力对于网格的作用。
对流方程的求解使用三次约束插值函数-半拉格朗日方法。该方法是根据网格的体积积分平均和面积积分平均建立插值函数，如图2所示。
其中，网格体积积分平均定义为：
网格的面积积分平均定义：
其中，δxi和δyi分别是网格在x方向和y方向的大小，φ(x,y,t)是网格(x,y)在t时刻某个物理量的值，φ(xi 1/2,y,t)是网格(x,y)x方向的右界面的某个物理量的值，φ(x,yj 1/2,t)是网格(x,y)y方向的上界面的某个物理量的值。
该方法是根据上述网格体积积分平均和面积积分平均的定义建立迎风的插值函数，在x方向上，在u＞0情况下,基于左侧插值的构造函数为：
其中，x∈[xi-(1/2),xi (1/2)]。
其中的四个系数根据网格物理量左右界面连续、网格体积守恒以及网格中心的斜率可以确定：
根据网格的左右两边连续有：
其中
根据网格体积积分的守恒有：
其中
根据网格的中间φil的一阶导数有：
其中，斜率di采用四阶的近似表示，形式如下：
进一步的，根据这四个限定可以计算得到四个系数为：
同样的也可以得到在u＜0情况下,基于右侧的插值构造函数φir。
进一步的，采用同样的方法得到左右两侧的三次插值函数，再利用半拉格朗日方法可以得到下一个时间步的界面位置的物理量面积积分平均的值(sia)。
其中使用的半拉格朗日的方法为：
其中，
其中，
其中，ui 1是x方向上i 1网格的x方向的速度值。ui是x方向上i网格的x方向的速度值，ui (1/2)是x方向上i网格的x方向右界面的速度值。
进一步的，下一个时刻的修正值为：
其中，τ为半拉格朗日粒子在dt时间间隔中的运动轨迹。求解下一个时间步的公式中的第二项
其中，
在求得网格的面积积分平均后，下一个时间步的网格体积积分平均的值为：
其中，gi (1/2),gi-(1/2)为在相应的网格的边界位置处tn-1-tn内的通量。
其中，
这样下一步的网格体积积分平均的值就进行更新了，通过该方法可以得到下一个时刻网格体积积分平均的密度和面积积分平均的密度，以及临时的网格体积积分平均的动量和面积积分平均的动量。
s3、求解欧拉坐标系下非对流部分；
在非对流部分计算之后还需要计算粘性项和源项，在计算粘性项和源项之后可以得到半离散的非对流部分方程，所以半离散的非对流部分方程为：
其中，mn 1是下一个时刻的动量值，m□是临时动量值，pn 1是下一个时刻的压力值。
进一步，得到：
进一步，得到压力的泊松方程：
其中，pn 1为下一个时刻的压力值，u□为对流部分计算得到的临时速度，ρn 1为下一个时刻的网格密度，m□是对流部分计算得到的临时动量。
通过求解压力泊松方程可以得到下一个时刻的压力的值。再更新下一个时刻的网格的动量面积积分平均和体积积分平均。
下一个时刻网格的动量面积积分平均为：
其中，
动量的下一个时间步的网格体积积分平均为：
其中，
至此网格部分计算完毕。
s4、网格部分和粒子计算部分物理量的传递；
混合方法中粒子只是分布在相界面中密度较大的液体的一侧。如图3所示，一个网格里面分布4×4的粒子，并且粒子只是分布在两相中密度比较大的不可压缩的液体的一侧，相界面附近粒子层的厚度一般大于四层网格的厚度。
在网格部分计算完毕，需要计算界面附近的粒子，从而根据粒子的运动来追踪相界面，同时根据粒子的核函数的密度来确定相界面附近的颜色函数的值。如图4所示，根据粒子所在的网格计算部分的位置，从网格中心插值得到粒子的初始速度。
其中
其中，re为截断半径其值为8.1倍的粒子的初始距离。
s5、粒子部分的修正计算；
对于粘性不可压缩流体，其质量方程和动量守恒方程如下：
其中，ρ为流体密度，t为时间，
在二维有限体积法中，一个移动粒子的控制体是一个圆，其离散区域计算如下：
v＝πr2＝l2
其中，v为离散的体积，r粒子的控制体半径，l为网格的大小。
进一步的，根据高斯定理中曲面积分与体积分的关系，在有限体积粒子法中，梯度算子和拉普拉斯算子计算如下：
其中，φs为粒子i和粒子j之间作用的物理量，
具体的，j对粒子i的作用面积
其中，n0是初始粒子数密度，
具体的，初始的核函数积分值n0计算如下：
具体的，两粒子之间的单位向量
具体的，核函数wij计算如下：
具体的，re为截断半径，当其中
进一步的，梯度算子和拉普拉斯算子计算如下：
其中，s是粒子的面积计算区域，其值为2πr，v是粒子的离散体积，n0是初始的粒子数密度，wij是粒子i和j之间的核函数权重值，
φs表示为一个线性函数，计算如下：
粒子部分的粒子初始速度计算如下：
具体的，使用核函数从网格部分计算得到的网格的中心速度插值得到粒子的临时速度和位置计算如下：
其中，
其中，i粒子与网格(m,n)的权重计算如下：
其中，re为截断半径，其值为8.1倍的粒子的初始距离。
在插值得到粒子的初始速度后可以得到粒子运动的临时位置：
其中，
得到初始的位置后使用有限体积粒子法来调整界面位置的粒子的位置。对于不可压缩的流体，根据连续性方程可以知道速度的散度需要为0。
速度的修正值由修正的压力的梯度得出：
可以得到如下的压力泊松方程：
具体的，修正的粒子速度与粒子压力的关系如下：
通过求解上方的压力泊松方程再进行粒子位置的更新。
s6、根据粒子部分计算结果更新界面附近网格的颜色函数值，如果达到模拟时间，完成多相流数值模拟计算。
粒子部分计算完毕之后，根据粒子的位置更新界面处的网格的颜色函数的值。
其中，θ为网格的颜色函数值，其中，
其中，re为截断半径，其值为8.1倍的粒子的初始距离。
本方法具有粒子法追踪相界面的准确性，又具有网格法的快速性。本方法适用于不包含传热情况下的多相流运动，其中在相界面变化剧烈，计算区域较大的时候最能发挥其优势。
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为了验证所开发的多相流计算的混合方法正确性和有效性以及快速性，采用数值模拟的方式进行验证，选取溃坝、单个气泡上升的算例进行模拟，并且与实验或者其他方法的模拟结果进行对比。
其中，溃坝的计算中网格数为40×40，计算中的重力加速度为9.8m/s2。计算区域中网格的大小为2×10-2m,粒子的大小为5×10-3m，边界条件是无滑移的边界。表1是溃坝工况的详细参数列表，图5是溃坝的计算区域的划分。
表1溃坝的参数列表
混合方法中溃坝的模拟计算的结果如图6所示。图6(a)中呈现了溃坝计算的网格部分的计算的中间时刻的结果。图6(b)是混合方法计算的粒子部分的中间时刻的结果。可以看到在0.3s左右的时候液体撞击右壁面。
在溃坝的计算的过程，记录液体撞击到右壁面之前液面移动的距离和液面的高度，对于液面移动的距离以及液面的高度采用无量纲化的分析，无量纲化的时间取为
但是本文计算需要的网格数是40×40，本方法的对于相界面的捕捉远远的比vof方法准确，本方法可以准确的捕捉两相的界面，并且不像其他的欧拉方法需要特殊的相界面捕捉技术，另外由于本文的粒子只是布置在相界面附近，所以计算的粒子数在计算的过程中最大节省了48％，从而实现相比拉格朗日方法的快速计算。
如图9所示为单个气泡上升中的结算的区域。单个气泡上升的模拟中，气泡中心距离左右壁面至少需要2.5倍的气泡初始半径宽度，气泡中心距离下底部也需要至少1.5倍的气泡初始半径高度，以此来避免左右墙壁以及底部壁面对于气泡运动的影响。
本文的计算的区域为100×100网格的区域。网格的大小是1×10-3m，对应的气泡的直径是15mm，气相和液相的物性参数如表2所示。
表2单个气泡上升的参数列表
本方法对于单个气泡上升的计算结果如图10所示。其中图10(a)为网格部分的计算结果，图10(b)为粒子部分的计算的结果。可以看到气泡最后在相应的莫顿数和艾维数下气泡的最终形状为球冠状。这与在相应的莫顿数和艾维数下grace经验关系图符合良好，计算的结果与其他方法的结果相比也基本一致。
本方法对于单个气泡上升的模拟中，相界面捕捉准确，比传统的欧拉方法中的vof更加的真实，而且不需要特殊的相界面的捕捉的技术，本算例如果使用拉格朗日的方法需要160000个粒子，这对于粒子法是十分巨大的粒子数，计算需要耗时几天，而本方法中由于粒子只是布置在相界面的附近，所以计算需要的粒子数最多只有3208个。大大的减少了计算需要的粒子数，本算例计算也只是需要几个小时就可以完成。
综上所述，本发明方法使用网格法和移动粒子法混合的方法来进行多相流的计算。混合方法中拉格朗日的粒子只是布置在多相流的相界面附近，这样可以减少计算需要的拉格朗日粒子数从而缩短计算的时间，又可以利用拉格朗日方法追踪界面的优势而不需要任何特殊的方法来追踪界面。数值计算中大部分计算区域使用网格法计算，只是在多相流的相界面附近使用移动粒子法计算，拉格朗日部分的粒子计算完毕后根据相界面处的粒子数密度修正网格法中的颜色函数的值，从而实现快速又准确的多相流的计算。
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。
1.一种多相流计算方法，其特征在于，首先确定计算区域，输入物性参数及相应的边界条件；使用约束插值守恒的半拉格朗日方法求解欧拉坐标系下控制方程的对流部分，求解欧拉坐标系下非对流部分的压力泊松方程，更新网格部分下个时间步的数值；将网格部分的速度传递给粒子计算部分进行修正，得到新的粒子位置，根据粒子部分计算结果更新界面附近网格的颜色函数值，达到模拟时间后完成多相流数值模拟计算。
2.根据权利要求1所述的多相流计算方法，其特征在于，对流部分中粘性流体的控制方程为：
其中，ρ是密度，
3.根据权利要求2所述的多相流计算方法，其特征在于，对流项方程组为：
使用约束插值守恒的半拉格朗日方法求解得到下一个时间步的
4.根据权利要求1所述的多相流计算方法，其特征在于，求解欧拉坐标系下非对流部分的压力泊松方程具体为：
网格部分粘性流体控制方程的半离散非对流项方程为：
非对流部分计算如下：
压力的泊松方程计算如下：
其中，pn 1为下一个时刻的压力值，u□为对流部分计算得到的临时速度，ρn 1为下一个时刻的网格密度，m□是对流部分计算得到的临时动量。
5.根据权利要求4所述的多相流计算方法，其特征在于，更新下一个时刻的网格的动量面积积分平均和体积积分平均具体为：
下一个时刻网格的动量面积积分平均为：
其中，
动量的下一个时间步的网格体积积分平均为：
其中，
6.根据权利要求1所述的多相流计算方法，其特征在于，粒子计算部分修正具体为：
粒子部分的粒子计算中，对于粘性不可压缩流体，其质量方程和动量守恒方程如下：
其中，ρ为流体密度，t为时间，
在二维有限体积法中，一个移动粒子的控制体是一个圆，其离散区域计算如下：
v＝πr2＝l2
其中，v为离散的体积，r粒子的控制体半径，l为网格的大小。
7.根据权利要求6所述的多相流计算方法，其特征在于，根据高斯定理中曲面积分与体积分的关系，在二维有限体积粒子法中，梯度算子和拉普拉斯算子计算如下：
其中，φs为粒子i和粒子j之间作用的物理量，
其中，n0是初始粒子数密度，
初始的核函数积分值n0计算如下：
两粒子之间的单位向量
核函数wij计算如下：
其中，re为截断半径，当其中
梯度算子和拉普拉斯算子计算如下：
其中，φs可以表示为一个线性函数计算如下：
8.根据权利要求7所述的多相流计算方法，其特征在于，粒子部分的粒子初始速度计算如下：使用核函数从网格部分计算得到的网格的中心速度插值得到粒子的临时速度和位置计算如下：
其中，
其中，re为截断半径。
9.根据权利要求6所述的多相流计算方法，其特征在于，粒子部分的粒子根据有限体积粒子法对于速度和位置的修正计算，对于不可压缩的流体，速度的散度计算如下：
根据粒子运动的动量方程确定速度的修正值与压力的梯度的关系如下：
压力泊松方程计算：
利用以上的梯度算子和拉普拉斯算子离散求解上述的压力泊松方程得到压力，修正的粒子速度与粒子压力的关系如下：
粒子的新的时刻的位置计算如下：
10.根据权利要求1所述的多相流计算方法，其特征在于，网格计算中的界面附近网格颜色函数值根据网格中粒子数密度计算，粒子部分计算完毕，根据粒子的位置更新界面处的网格的颜色函数的值，网格的颜色函数值θ为：
其中，n0是粒子按照4×4粒子布满网格时的粒子数密度，re为截断半径，其值为8.1倍的粒子初始距离。
技术总结