CN112861263A - 一种适用于可压缩两相流的计算模拟方法 - Google Patents

Abstract

本发明提出了一种适用于可压缩两相流的计算模拟方法，该方法采用五方程模型描述可压缩两相界面的流动，将五方程与刚性气体状态方程耦合使得方程组封闭。本发明采用有限体积法离散化方程组，采用三阶TVD格式对原始变量进行空间重构，采用HLLC近似黎曼求解器求解单元边界的数值通量，采用三阶显式的龙格‑库塔格式进行时间推进。本发明的计算模拟方法可以准确地捕捉到激波间断、物质接触间断和膨胀波，为两相流的数值模拟研究提供了可靠的数值模拟平台。

Description

一种适用于可压缩两相流的计算模拟方法

技术领域

本发明属于计算流体力学技术领域，具体涉及一种适用于可压缩两相流的计算模拟方法。

背景技术

可压缩两相流是一个活跃且具有挑战性的研究领域，在理论研究以及工程应用中都具有重要的意义，并且广泛存在于航空航天、化工生产、生物医学、农业生产等领域中。

在航空航天等领域，激波与液滴的相互作用是典型的两相流问题。当飞行器高速地穿过云层或者雨区时，飞机表面附近产生的激波与空气中悬浮的雨滴发生相互作用，雨滴在破碎过程中会使飞行器表面受到损伤，即雨滴侵蚀现象。深入研究激波与液滴相互作用的物理机制，可有效地减小雨滴侵蚀对飞行器带来的影响。此外，超燃冲压发动机燃烧室在超声速条件下的雾化研究是保证其在启动工况下可靠运行的关键。当气流以超声速流入燃烧室时，液体燃料喷射需要在极短的时间内完成雾化、蒸发以及复杂的燃烧过程，这涉及到激波与液体燃料的复杂相互作用。因此，研究并获得激波与液滴的物理规律对指导高速飞行器、超燃冲压发动机以及雾化装置的设计有着重要意义。

在叶轮机械与舰船结构等领域，空化是典型的两相流流动过程。空化是流体中由于压力极低而产生的气泡快速演化现象。气泡的突然产生以及突然溃灭会造成压力突变，进而导致严重的机械损伤，降低机械的使用寿命。深入研究空化现象的机理，对减缓机械磨损有着重要的指导意义。

在生物医疗等领域，可压缩两相流广泛存在于激波碎石、药物输运等方面。激波碎石是利用冲击波对体内的结石进行无创粉碎，所产生的压力脉冲与结石附近的气泡相互作用，使得结石碎裂，这是否成功取决于激波与气泡相互作用后气泡的动力学行为。因此，可压缩两相流的数值模拟具有重要的研究意义。

相比于单相流动，可压缩两相流动存在着复杂的波系结构，涉及到强激波和大密度比两相界面的时空演化过程，对数值模拟研究提出了巨大的挑战。因此，针对可压缩两相流界面问题的数值模拟方法，发展可靠的数值方法尤为重要。

发明内容

本发明要解决的问题是提供一种适用于可压缩两相流的计算模拟方法，该方法能够高效高精度地模拟两相流的流场演化过程。

为了实现上述目的，本发明提出的适用于可压缩两相流的计算模拟方法，具体步骤如下：

步骤1：根据需要计算模拟的可压缩两相流确定计算域，根据计算网格点数将计算域离散为均匀的网格点，并输入网格点信息，包括所有网格点上的速度、密度和压力；

步骤2：采用有限体积法对控制方程组进行离散，获得方程组的半离散形式；

步骤3：采用三阶TVD格式对离散后的控制方程组进行空间重构，获得单元边界两侧的原始变量；

步骤4：采用HLLC黎曼求解器求解单元边界处的通量；

步骤5：将刚性气体状态方程与流动求解器相结合以模拟气液两相；

步骤6：模拟步骤5中气液两相的混合规律；

步骤7：根据界面通量求解残差，并将半离散有限体积格式转化为时空全离散有限体积格式，全流场进行时间推进求解，得到最终的流场解。

进一步的，步骤2中，采用五方程模型描述两相流体的流动，所述五方程模型包括两相流体各自的连续方程、统一的动量方程、统一的能量方程和体积分数输运方程。

进一步的，步骤2中，不考虑粘性效应和表面张力的控制方程写为以下形式：

其中，U为控制方程的守恒变量，F为x方向的通量，G为y方向的通量，S为源项，定义如下：

ρ₁和ρ₂分别表示两相的密度，α₁和α₂分别表示两相的体积分数，下标1和2分别表示两相的物理量。ρ,u,v,p分别表示两相流体的混合密度，混合速度、混合压力，ρu和ρv表示两相流体的总动量，E表示两相流体的总动能。

进一步的，步骤2中，考虑计算的网格单元为[x_i-Δx/2,x_i+Δx/2]，不考虑体积输运方程的控制方程组的半离散形式为：

其中，

表示守恒变量的单元平均值，F_i-1/2和F_i+1/2分别表示单元左边界i-1/2和右边界i+1/2的数值通量；对于体积分数输运方程进行链式求导

得到方程的半离散形式为：

其中

为体积分数的单元平均值；公式中右端第一项为体积分数在单元边界处的数值通量，u_i+1/2和u_i-1/2为单元边界的速度；Δx为网格单元尺寸。

进一步的，步骤3中，满足TVD性质的三阶重构方法形式为

其中，χ_L、χ_R和

的计算公式如下：

进一步的，步骤4中，通过空间重构得到单元边界处的守恒变量U_l、U_r和数值通量F_l、F_r，再根据波系近似得到HLLC黎曼求解器求解的单元边界数值通量F_HLLC为：

F_l ^*和

分别为中间通量，S^*为中间波速：

U_l ^*和U_r ^*分别是中间状态向量，计算方式为：

S_l和S_r分别为左波速和右波速，根据公式

S_l＝min(S_l ^*,0)S_r＝max(S_r ^*,0)

进行计算：

而方程离散速度项u_HLLC根据公式

计算。

进一步的，步骤5中由刚性气体状态方程可知

p＝(γ-1)ρe-γπ_∞

而刚性气体状态方程的声速公式为：

得到总能E为

其中γ是比热比，π_∞为参考压力。

进一步的，步骤6中两相流体混合密度、混合体积分数和混合内能的表达式为：

ρ＝ρ₁α₁+ρ₂α₂

α₁+α₂＝1

ρe＝ρ₁α₁e₁+ρ₂α₂e₂

刚性气体状态方程中的两个状态参数Γ和Π的混合规则如下：

其中，π_∞,1和π_∞,2分别表示两相的参考压力。

有益效果

本发明提出了一种适用于可压缩两相流的计算模拟方法，该方法采用五方程模型描述可压缩两相界面的流动，将五方程与刚性气体状态方程耦合使得方程组封闭。本发明采用有限体积法离散化方程组，采用三阶TVD格式对原始变量进行空间重构，采用HLLC近似黎曼求解器求解单元边界的数值通量，采用三阶显式的龙格-库塔格式进行时间推进。本发明的计算模拟方法可以准确地捕捉到激波间断、物质接触间断和膨胀波，为两相流的数值模拟研究提供了可靠的数值模拟平台。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1为水下三氮化银***问题的示意图。

图2为t＝14μs时水下三氮化银***问题的实验结果(a)和计算模拟结果(b)。

图3为t＝28μs时水下三氮化银***问题的实验结果(a)和计算模拟结果(b)。

图4为t＝42μs时水下三氮化银***问题的实验结果(a)和计算模拟结果(b)。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

水下***是典型的可压缩两相流问题。例如本实施例中给出的水下三氮化银***问题的数值模拟计算区域如图1所示，单位为m。计算区域为{(x,y)|0≤x≤0.35,0≤y≤0.35}(单位:m)，分隔空气与水的水平面位于y＝0.175m，***核的直径为0.01m，中心位置为(0.175m，0.125m)。计算采用的网格分辨率为3000×1000，CFL数取值为0.1。该算例具体的初始条件如表1所示，其中，水和***核的参考压力π_∞均为6×10⁸，空气的参考压力π_∞为0；水和***核的比热γ取值为4.4，空气的比热比γ取值为1.4。计算区域的下边界采用反射边界条件，左右边界和上边界采用远场边界条件。

表1水下三氮化银***问题的计算参数

应用本发明的计算方法，对水下三氮化银***问题进行数值模拟，具体过程如下：

步骤1、网格信息输入以及初始化。

首先确定计算域{(x,y)|0≤x≤0.35,0≤y≤0.35}，令网格单元数为1000×1000，对计算域进行网格划分，将所有网格点上的速度、密度和压力信息输入计算方法程序中并完成初始化。

步骤2：采用有限体积法对控制方程组进行离散，获得方程组的半离散形式。

本发明采用五方程模型描述两相流体的流动，即假设两相具有相同的速度以及压力，采用各自的连续方程、统一的动量方程、统一的能量方程和体积分数输运方程。以二维为例，不考虑粘性效应和表面张力的控制方程可写为以下形式：

其中，U为控制方程的守恒变量，F为x方向的通量，G为y方向的通量，S为源项，上述4个量的具体定义如下：

上式中，ρ₁和ρ₂分别表示两相的密度，α₁和α₂分别表示两相的体积分数，下标1和2分别表示两相的物理量。ρ,u,v,p分别表示两相流体的混合密度，混合速度、混合压力，ρu和ρv表示两相流体的总动量，E表示两相流体的总动能。

以一维为例，考虑计算的网格单元为[x_i-Δx/2,x_i+Δx/2]，Δx为网格单元尺寸，则不考虑体积输运方程时，控制方程组的半离散形式如下：

其中，

表示守恒变量的单元平均值，F_i-1/2和F_i+1/2分别表示单元左边界i-1/2和右边界i+1/2的数值通量。对于体积分数输运方程，为使得其求解与欧拉方程组保持一致，对体积分数输运方程进行链式求导如下：

则上述方程的半离散形式可写为：

其中，

为体积分数的单元平均值。上式中右端第一项为体积分数在单元边界处的数值通量，u_i+1/2和u_i-1/2为单元边界的速度。

步骤3：采用三阶TVD格式对离散后的控制方程组进行空间重构，获得单元边界两侧的原始变量。

满足TVD性质的三阶重构方法形式如下：

其中，χ_L、χ_R和r_i ^L、r_i ^R的计算公式如下：

步骤4：采用HLLC黎曼求解器求解单元边界处的通量

通过使用HLLC黎曼求解器，可将网格单元附近的流动演变划分为三个波系，左右两测的波系分别为激波或膨胀波，中间的波系为接触间断。本步骤通过空间重构得到单元边界处的守恒变量U_l、U_r和数值通量F_l、F_r后，再根据波系近似可以得到HLLC黎曼求解器求解的单元边界数值通量F_HLLC为：

上式中，F_l ^*和

分别为中间通量，S^*为中间波速，计算公式如下：

F_k ^*＝F_k+S_k(U_k ^*-U_k),k＝l,r

上式中，U_l ^*和U_r ^*分别是中间状态向量，计算方式如下：

其中，S_l和S_r分别为左波速和右波速，可以采用以下公式进行计算：

S_l＝min(S_l ^*,0) S_r＝max(S_r ^*,0)

方程离散速度项u_HLLC的计算公式如下：

步骤5：将刚性气体状态方程与流动求解器相结合以模拟气液两相

刚性气体状态方程将流体的密度、压力以及内能紧密联系起来，当前被广泛应用于可压缩多组分流动模拟。本发明将刚性气体状态方程与流动求解器相结合来模拟气液两相。由刚性气体状态方程可知：

p＝(γ-1)ρe-γπ_∞

刚性气体状态方程的声速计算公式如下：

因此总能E可以表示为：

γ和π_∞均是根据经验确定的拟合参数。其中，γ是比热比，π_∞为参考压力，是与分子间引力相关的参数。考虑到所模拟的流体是不可混溶的，每一种流体均分别遵循刚性气体状态方程。对于空气，由于其参考压力为0，刚性气体状态方程即为理想气体方程。

步骤6：模拟两相的混合规律

在扩散界面区域内，流体混合物需满足一定的混合规则。对于空气、水和***核，混合密度、混合体积分数和混合内能的表达式定义如下：

ρ＝ρ₁α₁+ρ₂α₂

α₁+α₂＝1

ρe＝ρ₁α₁e₁+ρ₂α₂e₂

刚性气体状态方程中的两个状态参数Γ和Π的混合规则如下：

其中，π_∞,1和π_∞,2分别表示两相的参考压力，其为与分子之间引力相关的参数，本实施例中均为零。

步骤7：根据界面通量求解残差，并将半离散有限体积格式转化为时空全离散有限体积格式，全流场进行时间推进求解，得到最终的流场解。

图2-图4分别给出了在t＝14μs、t＝28μs和t＝42μs时刻水下三氮化银***问题的实验结果和计算模拟结果。可以看出，本发明的数值模拟方法可准确地捕捉到激波、膨胀波和两相界面等复杂的流动结构。将实验结果和数值模拟结果进行对比，可以看出数值模拟结果与实验情况吻合良好，可以很好地反应真实流动，验证了本发明所采用的数值模拟方法的有效性。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。本发明未详细描述内容为本领域技术人员公知技术。

Claims (8)

1.一种适用于可压缩两相流的计算模拟方法，其特征在于：包括以下步骤：

步骤1：根据需要计算模拟的可压缩两相流确定计算域，根据计算网格点数将计算域离散为均匀的网格点，并输入网格点信息，包括所有网格点上的速度、密度和压力；

步骤2：采用有限体积法对控制方程组进行离散，获得方程组的半离散形式；

步骤3：采用三阶TVD格式对离散后的控制方程组进行空间重构，获得单元边界两侧的原始变量；

步骤4：采用HLLC黎曼求解器求解单元边界处的通量；

步骤5：将刚性气体状态方程与流动求解器相结合以模拟气液两相；

步骤6：模拟步骤5中气液两相的混合规律；

步骤7：根据界面通量求解残差，并将半离散有限体积格式转化为时空全离散有限体积格式，全流场进行时间推进求解，得到最终的流场解。

2.根据权利要求1所述一种适用于可压缩两相流的计算模拟方法，其特征在于：步骤2中，采用五方程模型描述两相流体的流动，所述五方程模型包括两相流体各自的连续方程、统一的动量方程、统一的能量方程和体积分数输运方程。

3.根据权利要求2所述一种适用于可压缩两相流的计算模拟方法，其特征在于：步骤2中，不考虑粘性效应和表面张力的控制方程写为以下形式：

其中，U为控制方程的守恒变量，F为x方向的通量，G为y方向的通量，S为源项，定义如下：

ρ₁和ρ₂分别表示两相的密度，α₁和α₂分别表示两相的体积分数，下标1和2分别表示两相的物理量；ρ,u,v,p分别表示两相流体的混合密度，混合速度、混合压力，ρu和ρv表示两相流体的总动量，E表示两相流体的总动能。

4.根据权利要求3所述一种适用于可压缩两相流的计算模拟方法，其特征在于：步骤2中，考虑计算的网格单元为[x_i-Δx/2,x_i+Δx/2]，不考虑体积输运方程的控制方程组的半离散形式为：

其中，

表示守恒变量的单元平均值，F_i-1/2和F_i+1/2分别表示单元左边界i-1/2和右边界i+1/2的数值通量；对于体积分数输运方程进行链式求导

得到方程的半离散形式为：

其中

为体积分数的单元平均值；公式中右端第一项为体积分数在单元边界处的数值通量，u_i+1/2和u_i-1/2为单元边界的速度；Δx为网格单元尺寸。

5.根据权利要求4所述一种适用于可压缩两相流的计算模拟方法，其特征在于：步骤3中，满足TVD性质的三阶重构方法形式为

其中，χ_L、χ_R和r_i ^L、r_i ^R的计算公式如下：

6.根据权利要求5所述一种适用于可压缩两相流的计算模拟方法，其特征在于：步骤4中，通过空间重构得到单元边界处的守恒变量U_l、U_r和数值通量F_l、F_r，再根据波系近似得到HLLC黎曼求解器求解的单元边界数值通量F_HLLC为：

F_l ^*和

分别为中间通量，S^*为中间波速：

F_k ^*＝F_k+S_k(U_k ^*-U_k),k＝l,r

U_l ^*和U_r ^*分别是中间状态向量，计算方式为：

S_l和S_r分别为左波速和右波速，根据公式

S_l＝min(S_l ^*,0)S_r＝max(S_r ^*,0)

进行计算：

而方程离散速度项u_HLLC根据公式

计算。

7.根据权利要求6所述一种适用于可压缩两相流的计算模拟方法，其特征在于：步骤5中由刚性气体状态方程可知

p＝(γ-1)ρe-γπ_∞

而刚性气体状态方程的声速公式为：

得到总能E为

其中γ是比热比，π_∞为参考压力。

8.根据权利要求7所述一种适用于可压缩两相流的计算模拟方法，其特征在于：步骤6中两相流体混合密度、混合体积分数和混合内能的表达式为：

ρ＝ρ₁α₁+ρ₂α₂

α₁+α₂＝1

ρe＝ρ₁α₁e₁+ρ₂α₂e₂

刚性气体状态方程中的两个状态参数Γ和Π的混合规则如下：

其中，π_∞,1和π_∞,2分别表示两相的参考压力。

Images