可压缩气体与不可压缩液体多介质界面追踪数值方法

摘要

本发明公开一种可压缩气体与不可压缩液体多介质界面追踪数值方法，单一的气体和液体介质区域均采用间断伽辽金方法进行计算，同时由于多介质流场的时间步长一般是非均匀变化的，本发明在采用间断伽辽金方法求解不可压缩Navier-Stokes方程时，推导出基于非均匀时间步长的时间离散格式。多介质界面边界条件采用修正虚拟流体方法进行定义。在物质界面法向构造一种新的可压缩与不可压缩黎曼问题来预测界面状态，预测的界面状态不仅用于推进多介质界面，而且可以直接更新界面附近可压缩气体一侧的密度并外推得到对应的虚拟流体状态。通过多个经典数值算例的验证，发现对于气体中出现激波的多介质问题本发明能得到更为精确的数值模拟结果。

说明书

技术领域：

本发明涉及一种可压缩气体与不可压缩液体多介质界面追踪数值方法，其属于计 算流体力学，多相流数值模拟领域。

背景技术：

由于可压缩气体与不可压缩液体的物质属性存在较大差异，数值模拟相关气液多 介质问题存在较大困难。在早期的研究中，为降低计算难度气体和液体均作为可压缩介质 并且液体的状态方程为Tait方程。然而，由于界面两侧的物质属性差异较大，在物质界面处 易产生非物理振荡，而且由于液体中的声速较大会导致计算时间步长过小从而降低计算效 率。当气体和液体均作为不可压缩介质时，对于含有激波的高速流动问题难以求解，因而这 一方法只适应于求解低速多介质流动问题。而当气体作为可压缩介质而液体作为不可压缩 介质时，可以较好地考虑到物质属性的差异以及能够处理气体中含有激波的高速流动问 题。当然这一方法的前提是正确地处理多介质界面。

近二十年来，间断伽辽金(DG)方法在求解单介质流场上逐渐成为研究热点。在对 DG方法进行的研究中，较著名的是由Cockburn和Shu等人提出的龙格库塔间断伽辽金 (RKDG)方法：在时间离散上采用总变差减小的(TVD)高阶龙格库塔(RK)方法，在空间离散上 通过求解精确或近似黎曼问题获得网格界面处的数值通量，应用总变差有界(TVB)限制器 控制数值解的虚假振荡。相比之下，DG方法在不可压缩流场的求解上起步较晚。Liu等人最 初以涡量函数的形式研究高阶不可压缩DG方法。后来Cockburn等人提出了求解Stokes， Oseen，Navier-Stokes方程的局部间断伽辽金(LDG)方法。Ferrer等人提出了一种更简单高 效的求解非定常不可压Navier-Stokes方程的DG方法，该方法采用InteriorPenalty(IP) 空间离散方法求解椭圆形方程。在多介质界面处理上，早期基于γ模型、质量分数模型、 levelset函数的方法在多介质界面处往往会产生较大的数值耗散。然而，对于界面追踪方 法，界面是显式追踪的而且界面位置在计算过程中能够始终保持非常清晰。原始的虚拟流 体方法(GFM)主要用于定义界面边界条件从而使得每种介质可以单独求解。但是当压力或 者速度在界面附近存在较大的梯度时，该方法的表现效果较差。后来，在处理水气多介质界 面问题时，Fedkiw等提出新的虚拟流体方法(newGFM)。在该方法中，虚拟流体速度取决于 液体而虚拟流体压力由气体决定。NewGFM的优点是简单且易于推广到高维情况。然而界面 边界条件应当非线性地取决于界面两侧的物质属性，所以直接把真实流体状态作为虚拟流 体状态并不是一种严格的方法。而且，对于强激波打击物质界面问题，应当首先预测界面状 态来考虑界面处波系的相互作用。其后又有许多学者提出了一些改进版本的虚拟流体方 法，例如修正虚拟流体方法(MGFM)。该方法通过求解界面黎曼问题来定义虚拟流体状态。 MGFM在处理可压缩多介质问题上是非常有效的。

当可压缩气体与不可压缩液体相互作用时，由于物质界面两侧控制方程和物质属 性存在较大的差异，需要给出更为严格的物质界面边界条件。NewGFM给出的界面边界条件 过于简单，没有充分考虑到界面两侧物质属性的非线性影响。MGFM通过在物质界面处构造 黎曼问题来预测界面状态，在求解可压缩多介质问题上取得了较大的成功。然而当液体作 为不可压缩介质时，需要在界面处构造可压缩与不可压缩黎曼问题。所以本发明在界面的 法向构造出新的可压缩与不可压缩黎曼问题来预测界面状态，并且预测的界面状态可以直 接用于推进界面以及定义虚拟流体状态。在定义了物质界面边界条件后，可压缩气体和不 可压缩液体分别采用DG方法进行求解。同时由于现有的不可压缩DG方法的时间步长为均匀 不变的，这往往难以适应多介质问题变时间步长的要求，所以本发明推导出基于非均匀时 间步长的不可压缩时间离散格式。采用可压缩与不可压缩DG方法的目的是由于其具有较好 的空间紧致性，这使得所需定义的虚拟流体状态数目大为减少，在处理复杂多介质界面问 题上具有较大优势。

发明内容：

为了克服现有技术中的不足，本发明提供了一种可压缩气体与不可压缩液体多介 质界面追踪数值方法。

本发明采用如下技术方案：一种可压缩气体与不可压缩液体多介质界面追踪数值 方法，其包括如下步骤：

步骤1:利用一系列标志点离散多介质界面，在标志点法向构造可压缩与不可压缩 黎曼问题预测界面状态；

步骤2:直接利用标志点上预测的界面状态来更新界面附近可压缩气体一侧的密 度，并分别外推得到对应的虚拟流体状态；

步骤3:采用龙格库塔间断伽辽金方法求解可压缩气体，对于不可压缩Navier- Stokes方程，推导基于非均匀时间步长的时间离散格式，采用不可压缩间断伽辽金方法对 液体介质进行计算；

步骤4:由标志点上预测的界面状态推进多介质界面，重构界面，得到新的界面位 置和标志点并根据新的界面确定整个多介质流场的数值解；

步骤5:判断是否满足计算终止条件，如不满足此条件，返回步骤1继续进行下一时 刻的数值计算。

进一步地，在可压缩与不可压缩黎曼问题的构造上，把不可压缩液体当作声速趋 于无穷大的可压缩液体，当所有介质为可压缩时，状态方程可以统一表示为：

p＝(γ-1)ρe-γB

其中B是物质参数，对于理想气体，B为零，考虑黎曼问题的初始条件：

$Q_0=\left(\begin{array}{c}{Q}_L\\ Q_R\end{array}\right)$

其中Q_L＝[ρ_L,u_L,p_L,γ_L,B_L]^T,Q_R＝[ρ_R,u_R,p_R,γ_R,B_R]^T，不失一般性，假设左侧为气 体，右侧为液体，首先采用双激波近似方法求解该黎曼问题：

${\rho }_L^I=\frac{{\gamma }_L(p_L+B_L)+{\gamma }_L(p_I+B_L)+p_I-p_L}{{\gamma }_L(p_L+B_L)+{\gamma }_L(p_I+B_L)+p_L-p_I}{\rho }_L$

${\rho }_R^I=\frac{{\gamma }_R(p_R+B_R)+{\gamma }_R(p_I+B_R)+p_I-p_R}{{\gamma }_R(p_R+B_R)+{\gamma }_R(p_I+B_R)+p_R-p_I}{\rho }_R$

$\frac{p^I-p_L}{W_L}+\frac{p^I-p_R}{W_R}+u_R-u_L=0$

其中：

$W_L=\sqrt{\frac{{\rho }_L{\rho }_L^I(p^I-p_L)}{{\rho }_L^I-{\rho }_L}},W_R=\sqrt{\frac{{\rho }_R{\rho }_R^I(p^I-p_R)}{{\rho }_R^I-{\rho }_R}}$

而预测的界面速度定义为：

$u^I=\frac{1}{2}(u_L+u_R)+\frac{1}{2}\left(f_R\right(p^I,Q_R)-f_L(p^I,Q_L\left)\right)$

其中：

$f_L(p^I,Q_L)=\frac{\sqrt{2}(p^I-p_L)}{\sqrt{{\rho }_L[{\gamma }_L(p_L+B_L)+{\gamma }_L(p^I+B_L)+p^I-p_L]}}$

$f_R(p^I,Q_R)=\frac{\sqrt{2}(p^I-p_R)}{\sqrt{{\rho }_R[{\gamma }_R(p_R+B_R)+{\gamma }_R(p^I+B_R)+p^I-p_R]}}$

由于液体的声速：

$a_R=\sqrt{k_R/{\rho }_R}=\sqrt{{\gamma }_R(p_R+B_R)/{\rho }_R}$

假设为无穷大，但是ρ_R为有限值，所以：

k_R＝γ_R(p_R+B_R)→+∞

将其代入界面速度的定义式以及基于双激波近似方法得到的黎曼解中，可以得 到：

$\underset{k_R\to +\infty }{\lim }{\rho }_R^I={\rho }_R$

$\underset{k_R\to +\infty }{\lim }{W}_R=+\infty$

f_R(p^I,Q_R)＝0

$\frac{p^I-p_L}{W_L}+u_R-u_L=0$

$u^I=\frac{1}{2}(u_L+u_R)-\frac{1}{2}{f}_L(p^I,Q_L)$

对于左侧气体有：

u^I＝u_L-f_L(p^I,Q_L)

可以得到：

u^I＝u_R

采用迭代方法求解式可以得到预测的界面压力p^I，代入式 ${\rho }_L^I=\frac{{\gamma }_L(p_L+B_L)+{\gamma }_L(p_I+B_L)+p_I-p_L}{{\gamma }_L(p_L+B_L)+{\gamma }_L(p_I+B_L)+p_L-p_I}{\rho }_L$得到同时有${\rho }_R^I={\rho }_R,$u^I＝u_R，这样即求得可压 缩与不可压缩黎曼问题的解。

进一步地，采用不可压缩间断伽辽金方法求解液体介质时，推导出一种基于非均 匀时间步长的时间离散格式，对Navier-Stokes方程的时间导数项采用二阶后退差分离散：

${\left(\frac{\partial V}{\partial t}\right)}^{n+1}\approx {\gamma }_0{V}^{n+1}-{\alpha }_0{V}^n-{\alpha }_1{V}^{n-1}$

其中γ₀＝α₀+α₁，下标“n-1”,“n”和“n+1”分别表示t^n-1，tⁿ和tⁿ⁺¹时刻，参数γ₀，α₀， α₁为：

${\gamma }_0=\frac{2{\mathrm{\Delta t}}^n+{\mathrm{\Delta t}}^{n-1}}{{\mathrm{\Delta t}}^n({\mathrm{\Delta t}}^{n-1}+{\mathrm{\Delta t}}^n)}$

${\alpha }_0=\frac{{\mathrm{\Delta t}}^n+{\mathrm{\Delta t}}^{n-1}}{{\mathrm{\Delta t}}^n{\mathrm{\Delta t}}^{n-1}}$

${\alpha }_1=\frac{-{\mathrm{\Delta t}}^n}{{\mathrm{\Delta t}}^{n-1}({\mathrm{\Delta t}}^{n-1}+{\mathrm{\Delta t}}^n)}$

其中Δt^n-1＝tⁿ-t^n-1,Δtⁿ＝tⁿ⁺¹-tⁿ，时间导数项可以进一步表示为：

$\left(\begin{array}{c}{\left(\frac{\partial V}{\partial t}\right)}^{n+1}={\alpha }_0(V^{n+1}-V^n)+{\alpha }_1(v^{n+1}-V^{n-1})\\ ={\alpha }_0{\int }_{t^n}^{t^{n+1}}\frac{\partial V}{\partial t}dt+{\alpha }_1{\int }_{t^{n-1}}^{t^{n+1}}\frac{\partial V}{\partial t}dt\\ ={\alpha }_0{\int }_{t^n}^{t^{n+1}}fdt+{\alpha }_1{\int }_{t^{n-1}}^{t^{n+1}}fdt\end{array}\right)$

其中f是Navier-Stokes方程的非线性项、粘性项和压力项，非线性项采用显式处 理，粘性项和压力项采用隐式处理，可以得到不可压缩Navier-Stokes方程的时间离散格式 为：

${\gamma }_0{V}^{n+1}-{\alpha }_0{V}^n-{\alpha }_1{V}^{n-1}=-{\beta }_0{(V·▿V)}^n-{\beta }_1{(V·▿V)}^{n-1}-▿{\bar{p}}^{n+1}+v{▿}^2{V}^{n+1}$

其中：

${\beta }_0={\alpha }_0{\mathrm{\Delta t}}^n=\frac{{\mathrm{\Delta t}}^n+{\mathrm{\Delta t}}^{n-1}}{{\mathrm{\Delta t}}^{n-1}}$

${\beta }_1={\alpha }_1({\mathrm{\Delta t}}^{n-1}+{\mathrm{\Delta t}}^n)=-\frac{{\mathrm{\Delta t}}^n}{{\mathrm{\Delta t}}^{n-1}}$

通过设置Δtⁿ≠Δt^n-1，即可得到基于非均匀时间步长的不可压缩Navier-Stokes 方程的时间离散格式，满足多介质数值模拟中非均匀时间步长的要求。

本发明具有如下有益效果：

(1)通过构造新的可压缩与不可压缩黎曼问题来预测物质界面状态，可以在液体 为不可压缩时应用MGFM定义更为精确的界面边界条件；

(2)推导出基于非均匀时间步长的不可压缩时间离散格式，使之能够应用于多介 质中不可压缩流体的数值计算；

(3)采用空间紧致性较好的可压缩与不可压缩DG方法进行计算使得在物质界面附 近所需虚拟流体状态信息大为减少，提高了流场计算的精确性；

(4)预测的物质界面状态不仅用于推进界面，而且可直接用于定义虚拟流体状态， 编程简单，计算效率较高，且易于向高维推广。

附图说明：

图1：标志点处黎曼问题的构造。

图2：可压缩流体密度的更新。

图3：数值结果(算例8.1)：压力(左)和速度(右)。

图4：数值结果(算例8.2)：压力(左)和速度(右)。

图5：数值结果(算例8.3)：压力(左)和速度(右)。

图6：数值结果(算例8.4)：压力(左)和速度(右)。

图7：数值结果(算例8.5)：压力(左)和速度(右)。

图8：不可压密度ρ＝1000kg/m³的数值结果(算例8.6)：压力(左)和x方向速度 (右)。

图9：不可压密度ρ＝1000kg/m³的压力云图(算例8.6)。

图10：不可压密度ρ＝10kg/m³的数值结果(算例8.6)：压力(左)和x方向速度(右)。

图11：不可压密度ρ＝10kg/m³的界面形状(算例8.6)。

图12：不可压密度ρ＝1000kg/m³的数值结果(算例8.7)：压力(左)和x方向速度 (右)。

图13：不可压密度ρ＝1000kg/m³的压力云图(算例8.7)。

图14：不可压密度ρ＝10kg/m³的数值结果(算例8.7)：压力(左)和x方向速度(右)。

图15：不可压密度ρ＝10kg/m³的界面形状(算例8.7)。

具体实施方式：

下面结合附图和算例对发明内容作进一步说明。

(1)控制方程

可压缩气体的控制方程为Euler方程：

$\frac{\partial U}{\partial t}+▿·F\left(U\right)=0---\left(1\right)$

其中：U＝[ρ,ρu,ρv,E]^T,F(U)＝[F₁(U),F₂(U)],F₁(U)＝[ρu,ρu²+p,ρuv,(E+p)u]^T, F₂(U)＝[ρv,ρuv,ρv²+p,(E+p)v]^T。ρ是密度，u是x方向速度，v是y方向速度，p是压力，E是单 位体积的总能量，其定义为：

E＝ρe+ρ(u²+v²)/2(2)

其中e是单位质量的内能。状态方程为：

p＝(γ-1)ρe(3)

对于理想气体γ是比热比。

不可压缩液体的控制方程为Navier-Stokes方程：

$\frac{\partial V}{\partial t}=-V·▿V-▿\bar{p}+\frac{\mu }{\rho }{▿}^{2}V---\left(4\right)$

$▿·V=0---\left(5\right)$

其中，V＝(u,v)是速度矢量，μ是粘性系数。ρ和μ在不可压缩液体中假设 为常数。Navier-Stokes方程的解记为：

(2)自适应时间步长

多介质全流场的时间步长取为可压缩气体与不可压缩液体时间步长的最小值：

Δt＝min(Δt_c,Δt_i)(6)

对于可压缩气体，时间步长定义为：

${\mathrm{\Delta t}}_c=\frac{C_f}{\left(\right(\left|u\right|+a)/\Delta x+(\left|v\right|+a)/\Delta y)}---\left(7\right)$

其中C_f是CFL系数，a是局部声速，Δx和Δy是网格步长。对于不可压缩液体，时间

步长定义为：

${\mathrm{\Delta t}}_i\approx O\left(\frac{\bar{L}}{\bar{U}{k}_i^2}\right)---\left(8\right)$

其中是特征长度(网格尺寸),是特征速度，k_i是不可压缩DG方法的多项式阶 数。

(3)求解Euler方程的RKDG方法

对于矩形网格I_i,j＝[x_i-1/2,x_i+1/2]×[y_j-1/2,y_j+1/2]，网格中心记为(x_i,y_j)，网格步 长为Δx＝x_i+1/2-x_i-1/2，Δy＝y_j+1/2-y_j-1/2。数值解U^h(x,y,t)和基函数由 给定(其中是网格I_i,j上多项式次数小于或等于k_c的多项 式集合)。网格I_i,j上的当地正交基函数定义为：

$\begin{array}{c}{\phi }_{i,j}^0(x,y)=1\\ {\phi }_{i,j}^1(x,y)=\frac{x-x_i}{\Delta x/2}\\ {\phi }_{i,j}^2(x,y)=\frac{y-y_j}{\Delta y/2}\\ {\phi }_{i,j}^3(x,y)={\phi }_{i,j}^1(x,y){\phi }_{i,j}^2(x,y)\\ {\phi }_{i,j}^4(x,y)={\phi }_{i,j}^1(x,y){\phi }_{i,j}^1(x,y)-\frac{1}{3}\\ {\phi }_{i,j}^5(x,y)={\phi }_{i,j}^2(x,y){\phi }_{i,j}^2(x,y)-\frac{1}{3}\\ \mathrm{...}\end{array}---\left(9\right)$

数值解U^h(x,y,t)可以表示为：

$U^h(x,y,t)=\underset{l=0}{\overset{L}{\Sigma }}{U}_{i,j}^{\left(l\right)}\left(t\right){\phi }_{i,j}^l(x,y),\forall (x,y)\in I_{i,j}---\left(10\right)$

其中定义为：

$U_{i,j}^{\left(l\right)}\left(t\right)=\frac{1}{a_l}{\int }_{I_{i,j}}{U}^h(x,y,t){\phi }_{i,j}^l(x,y)dxdy,l=0,1,...,L---\left(11\right)$

其中把式(10)代入式(1),对式(1)乘以基函数并在网格 I_i,j上分部积分，可以得到：

$\begin{array}{c}\frac{{\mathrm{dU}}_{i,j}^{\left(l\right)}\left(t\right)}{dt}=\frac{1}{a_l}(-\underset{e\in \partial I_{i,j}}{\Sigma }{\int }_{e}F(U^h\left(x,t,t\right))·n_{e,I_{i,j}}{\phi }_{i,j}^l(x,y)d\Gamma \\ +{\int }_{I_{i,j}}F\left(U^h\right(x,y,t\left)\right)·{\mathrm{grad\phi }}_{i,j}^l(x,y)dxdy),l=0,1,\mathrm{...},L\end{array}---\left(12\right)$

其中是网格I_i,j上边e的单位法向量。采用Lax- Friedrichs通量近似。右端积分项采用高斯求积法求解。式(12)的半离散形式表示为：

U_t＝L(U)(13)

采用三阶TVDRK方法进行时间离散：

$\begin{array}{c}{U}^{\left(1\right)}=u^n+\Delta tL\left(U^n\right)\\ U^{\left(2\right)}=\frac{3}{4}{U}^n+\frac{1}{4}{U}^{\left(1\right)}+\frac{1}{4}\Delta tL\left(U^{\left(1\right)}\right)\\ U^{n+1}=\frac{1}{3}{U}^n+\frac{2}{3}{U}^{\left(2\right)}+\frac{2}{3}\Delta tL\left(U^{\left(2\right)}\right)\end{array}---\left(14\right)$

当流场出现间断时应用斜率限制器克服可能出现的虚假振荡问题。对于标量方 程，记：

$U_{i+1/2,j}^{-}=U_{i,j}^{\left(0\right)}\left(t\right)+{\bar{U}}_{i,j}^x,U_{i-1/2,j}^{+}=U_{i,j}^{\left(0\right)}\left(t\right)-{\overline{\stackrel{‾}{U}}}_{i,j}^x$

$U_{i,j+1/2}^{-}=U_{i,j}^{\left(0\right)}\left(t\right)+{\bar{U}}_{i,j}^y,U_{i,j-1/2}^{+}=U_{i,j}^{\left(0\right)}\left(t\right)-{\overline{\stackrel{‾}{U}}}_{i,j}^y$

其中：

$\begin{array}{c}{\bar{U}}_{i,j}^x=\underset{l=1}{\overset{L}{\Sigma }}{U}_{i,j}^{\left(l\right)}\left(t\right){\phi }_{i,j}^l(x_{i+1/2},y_j),{\overline{\stackrel{‾}{U}}}_{i,j}^x=-\underset{l=1}{\overset{L}{\Sigma }}{U}_{i,j}^{\left(l\right)}\left(t\right){\iota }_{i,j}^l(x_{i-1/2},y_j)\\ {\bar{U}}_{i,j}^y=\underset{l=1}{\overset{L}{\Sigma }}{U}_{i,j}^{\left(l\right)}\left(t\right){\phi }_{i,j}^l(x_i,y_{j+1/2}),{\overline{\stackrel{‾}{U}}}_{i,j}^y=-\underset{l=1}{\overset{L}{\Sigma }}{U}_{i,j}^{\left(l\right)}\left(t\right){\phi }_{i,j}^l(x_i,y_{j-1/2})\end{array}---\left(15\right)$

对进行修正：

$\begin{array}{c}{\bar{U}}_{i,j}^{x\left(\mathrm{mod}\right)}=\bar{m}({\bar{U}}_{i,j}^x,U_{i+1,j}^{\left(0\right)}(t)-U_{i,j}^{\left(0\right)}(t),U_{i,j}^{\left(0\right)}(t)-U_{i-1,j}^{\left(0\right)}(t\left)\right)\\ {\overline{\stackrel{‾}{U}}}_{i,j}^{x\left(\mathrm{mod}\right)}=\bar{m}({\overline{\stackrel{‾}{U}}}_{i,j}^x,U_{i+1,j}^{\left(0\right)}(t)-U_{i,j}^{\left(0\right)}(t),U_{i,j}^{\left(0\right)}(t)-U_{i-1,j}^{\left(0\right)}(t\left)\right)\\ {\bar{U}}_{i,j}^{y\left(\mathrm{mod}\right)}=\bar{m}({\bar{U}}_{i,j}^y,U_{i,j+1}^{\left(0\right)}(t)-U_{i,j}^{\left(0\right)}(t),U_{i,j}^{\left(0\right)}(t)-U_{i,j-1}^{\left(0\right)}(t\left)\right)\\ {\overline{\stackrel{‾}{U}}}_{i,j}^{y\left(\mathrm{mod}\right)}=\bar{m}({\overline{\stackrel{‾}{U}}}_{i,j}^x,U_{i,j+1}^{\left(0\right)}(t)-U_{i,j}^{\left(0\right)}(t),U_{i,j}^{\left(0\right)}(t)-U_{i,j-1}^{\left(0\right)}(t\left)\right)\end{array}---\left(16\right)$

其中是修正minmod函数：

其中是TVB限制器常量。对于和上式中Δx应替换为Δy。对于方 程组，需要对限制器进行当地特征分解。

(4)求解Navier-Stokes方程的DG方法

不可压缩区域记为Ω，其边界可以为Dirichlet类型(Γ_D)或者Neumann类型 (Γ_N)等。对不可压缩区域进行矩形划分Ω_h＝{K}，外边界记为内边界记为Γ_h。数值解 W^h(x,y,t)和基函数由给定(其中是网格K上多项式次数 小于或等于k_i的多项式集合)。记网格K＝I_i,j，数值解W^h(x,y,t)可以表示为：

$W^h(x,y,t)=\underset{m=0}{\overset{M}{\Sigma }}{W}_{i,j}^{\left(l\right)}\left(t\right){\phi }_{i,j}^l(x,y),\forall (x,y)\in I_{i,j}---\left(18\right)$

其中是多项式系数，基函数的定义与式(9)一致。在现有的不可压 缩DG方法中，求解单介质流动问题的时间步长是均匀不变的。然而对于多介质问题，总的时 间步长(式(6))应由可压缩介质与不可压缩介质共同决定，所以有必要推导出基于非均匀 时间步长的不可压缩时间离散格式。对Navier-Stokes方程的时间导数项采用二阶后退差 分离散：

${\left(\frac{\partial V}{\partial t}\right)}^{n+1}\approx {\gamma }_0{V}^{n+1}-{\alpha }_0{V}^n-{\alpha }_1{V}^{n-1}---\left(19\right)$

其中γ₀＝α₀+α₁。下标“n-1”,“n”和“n+1”分别表示t^n-1，tⁿ和tⁿ⁺¹时刻。参数γ₀，α₀， α₁为：

${\gamma }_0=\frac{2{\mathrm{\Delta t}}^n+{\mathrm{\Delta t}}^{n-1}}{{\mathrm{\Delta t}}^n({\mathrm{\Delta t}}^{n-1}+{\mathrm{\Delta t}}^n)}---\left(20\right)$

${\alpha }_0=\frac{{\mathrm{\Delta t}}^n+{\mathrm{\Delta t}}^{n-1}}{{\mathrm{\Delta t}}^n{\mathrm{\Delta t}}^{n-1}}---\left(21\right)$

${\alpha }_1=\frac{-{\mathrm{\Delta t}}^n}{{\mathrm{\Delta t}}^{n-1}({\mathrm{\Delta t}}^{n-1}+{\mathrm{\Delta t}}^n)}---\left(22\right)$

其中Δt^n-1＝tⁿ-t^n-1,Δtⁿ＝tⁿ⁺¹-tⁿ。式(19)进一步表示为：

$\begin{array}{c}{\left(\frac{\partial V}{\partial t}\right)}^{n+1}={\alpha }_0(V^{n+1}-V^n)+{\alpha }_1(v^{n+1}-V^{n-1})\\ ={\alpha }_0{\int }_{t^n}^{t^{n+1}}\frac{\partial V}{\partial t}dt+{\alpha }_1{\int }_{t^{n-1}}^{t^{n+1}}\frac{\partial V}{\partial t}dt\\ ={\alpha }_0{\int }_{t^n}^{t^{n+1}}fdt+{\alpha }_1{\int }_{t^{n-1}}^{t^{n+1}}fdt\end{array}---\left(23\right)$

其中f是式(4)的右端项。非线性项采用显式处理，粘性项和压力项采用隐式处理。 这样可以得到不可压缩Navier-Stokes方程的时间离散格式：

${\gamma }_0{V}^{n+1}-{\alpha }_0{V}^n-{\alpha }_1{V}^{n-1}=-{\beta }_0{(V·▿V)}^n-{\beta }_1{(V·▿V)}^{n-1}-▿{\bar{p}}^{n+1}+v{▿}^2{V}^{n+1}---\left(24\right)$

其中：

${\beta }_0={\alpha }_0{\mathrm{\Delta t}}^n=\frac{{\mathrm{\Delta t}}^n+{\mathrm{\Delta t}}^{n-1}}{{\mathrm{\Delta t}}^{n-1}}---\left(25\right)$

${\beta }_1={\alpha }_1({\mathrm{\Delta t}}^{n-1}+{\mathrm{\Delta t}}^n)=-\frac{{\mathrm{\Delta t}}^n}{{\mathrm{\Delta t}}^{n-1}}---\left(26\right)$

可以发现，如果令Δtⁿ≠Δt^n-1，式(24)即为非均匀时间步长的时间离散格式。引 入中间变量和式(24)被分裂成：

${\gamma }_0\bar{V}-{\alpha }_0{V}^n-{\alpha }_1{V}^{n-1}=-{\beta }_0{(V·▿V)}^n-{\beta }_1{(V·▿V)}^{n-1}---\left(27\right)$

${\gamma }_0(\overline{\stackrel{‾}{V}}-\bar{V})=-▿{\bar{p}}^{n+1}---\left(28\right)$

${\gamma }_0(V^{n+1}-\overline{\stackrel{‾}{V}})=v{▿}^2{V}^{n+1}---\left(29\right)$

对式(28)求散度并考虑得到：

$-{▿}^2{\bar{p}}^{n+1}=-{\gamma }_0▿·\bar{V}---\left(30\right)$

式(29)重新表示为：

$-{▿}^2{V}^{n+1}+\frac{{\gamma }_0}{\nu }{V}^{n+1}=\frac{{\gamma }_0}{v}\overline{\stackrel{‾}{V}}---\left(31\right)$

由式(30)求解更新中间变量并根据式(31)求解Vⁿ⁺¹。采用对称IP伽辽金 (SIPG)方法求解式(30)和(31)：

其中，q∈H¹(Ω)为标量函数(可以直接推广到矢量函数)，α是波数，n是上外法 向矢量，g∈L²(Ω)，L_D∈H^1/2(Γ_D)和L_N∈L²(Γ_N)表示边界条件。式(32)的弱解可以通过下式 求得：

$a(q^h,\phi )=l\left(\phi \right),\forall \phi \in D_h^{k_i}\left({\Omega }_h\right)---\left(33\right)$

$\begin{array}{c}a(q^h,\phi )=\underset{K\in {\Omega }_h}{\Sigma }{\int }_K(▿q^h·▿\phi +{\mathrm{\alpha q}}^h\phi )dxdy\\ -\underset{e\in {\Gamma }_h\cup {\Gamma }_D}{\Sigma }{\int }_e\left(\right\{\{▿q^h\}\}·n_e[\left[\phi \right]]+\{\{▿\phi \}\}·n_e[\left[q^h\right]]-\lambda [\left[q^h\right]\left]\right[\left[\phi \right]\left]\right)ds\end{array}---\left(34\right)$

$l\left(\phi \right)=\underset{K\in {\Omega }_h}{\Sigma }{\int }_K\left(g\phi \right)dxdy-\underset{e\in {\Gamma }_D}{\Sigma }{\int }_e+(▿\phi ·n_e-\lambda \phi )L_{D}ds\underset{e\in {\Gamma }_N}{\Sigma }{\int }_e{\mathrm{\phi L}}_{N}ds---\left(35\right)$

其中，n_e是边e的外法向矢量，[[·]]和{{·}}表示相邻网格K₁和K₂之间的跳跃值 和平均值：

$\left[\right[·\left]\right]=(·{|}_{K_1})-(·{|}_{K_2}),\left\{\right\{·\left\}\right\}=\left(\right(·{|}_{K_1})+\left(\right(·{|}_{K_2}))/2,\forall e=\partial K_1\cap \partial K_2---\left(36\right)$

在边界上有：

$\left[\right[·\left]\right]=\left\{\right\{·\left\}\right\}=(·{|}_{K_1}),\forall e=\partial K_1\cap \partial {\Omega }_h---\left(37\right)$

参数λ定义为：

$\lambda =\frac{(k_i+1)(k_i+2)}{2}{\max }_K\left({\eta }_K\right)---\left(38\right)$

其中η_K为网格K的面积与周长之比。

(5)界面追踪方法

如图1所示，在tⁿ时刻可压缩气体(介质1)和不可压缩液体(介质2)被界面Π分开。 标记点是界面和网格线的交点。以标记点E为例，点A(x_A,y_A)和点B(x_B,y_B)在不同介质中,与 标记点E相距为Δn。这里：

$\Delta n={[{\left(\frac{N_{Ex}}{\Delta x}\right)}^2+{\left(\frac{N_{Ey}}{\Delta y}\right)}^2]}^{\frac{1}{2}}---\left(39\right)$

x_A＝x_E-Δn·N_Ex,y_A＝y_E-Δn·N_Ey

(40)

x_B＝x_E+Δn·N_Ex,y_B＝y_E+Δn·N_Ey

其中是标记点E的单位法向矢量。点A和点B的状态由式(10)和式 (18)得到：

$\begin{array}{c}{U}_A=\underset{l=0}{\overset{L}{\Sigma }}{U}_{K_A}^{\left(l\right)}\left(t\right){\phi }_{i,j}^l(x_A,y_A)\\ W_A=\underset{l=0}{\overset{M}{\Sigma }}{W}_{K_B}^{\left(l\right)}\left(t\right){\phi }_{i,j}^l(x_B,y_B)\end{array}---\left(41\right)$

其中K_A和K_B分别是包含点A和点B的网格。计算点A和点B的密度、法向速度和压力， 分别表示为和沿着标志点E的法线方向构造可压缩与 不可压缩黎曼问题，其初始条件为：

$Q_E=\left(\begin{array}{c}{Q}_A\\ Q_B\end{array}\right)---\left(42\right)$

求解该黎曼问题并记其解为其中上标I表示界面，下标L和R 表示界面左右两侧。由于不可压缩介质的刚性较强，界面的切向速度应由液体决定。此处标 志点E的切向速度定义为：其中是点B的切向速度。这样标志点E以速度v_f＝(u^I, v^I)被推进到下一时刻。其他标志点的推进方式是类似的。

(6)可压缩与不可压缩黎曼问题

在修正虚拟流体方法中，所有的介质为可压缩的，通过构造可压缩黎曼问题来预 测界面状态。然而当液体作为不可压缩介质时，应当在界面处构造可压缩与不可压缩黎曼 问题。考虑到当马赫数趋于零时，可压缩流体的解趋于不可压缩流体的解。本发明把不可压 缩液体当作声速趋于无穷大的可压缩液体来构造可压缩与不可压缩黎曼问题。当所有介质 为可压缩时，状态方程可以统一表示为：

p＝(γ-1)ρe-γB(43)

其中B是物质参数。对于理想气体，B为零，上式退化为式(3)。考虑黎曼问题的初始 条件：

$Q_0=\left(\begin{array}{c}{Q}_L\\ Q_R\end{array}\right)---\left(44\right)$

其中Q_L＝[ρ_L,u_L,p_L,γ_L,B_L]^T,Q_R＝[ρ_R,u_R,p_R,γ_R,B_R]^T。不失一般性，假设左侧为气 体，右侧为液体。首先采用双激波近似方法求解该黎曼问题：

${\rho }_L^I=\frac{{\gamma }_L(p_L+B_L)+{\gamma }_L(p_I+B_L)+p_I-p_L}{{\gamma }_L(p_L+B_L)+{\gamma }_L(p_I+B_L)+p_L-p_I}{\rho }_L---\left(45\right)$

${\rho }_R^I=\frac{{\gamma }_R(p_R+B_R)+{\gamma }_R(p_I+B_R)+p_I-p_R}{{\gamma }_R(p_R+B_R)+{\gamma }_R(p_I+B_R)+p_R-p_I}{\rho }_R---\left(46\right)$

$\frac{p^I-p_L}{W_L}+\frac{p^I-p_R}{W_R}+u_R-u_L=0---\left(47\right)$

其中：

$W_L=\sqrt{\frac{{\rho }_L{\rho }_L^I(p^I-p_L)}{{\rho }_L^I-{\rho }_L}},W_R=\sqrt{\frac{{\rho }_R{\rho }_R^I(p^I-p_R)}{{\rho }_R^I-{\rho }_R}}---\left(48\right)$

而预测的界面速度定义为：

$u^I=\frac{1}{2}(u_L+u_R)+\frac{1}{2}\left(f_R\right(p^I,Q_R)-f_L(p^I,Q_L\left)\right)---\left(49\right)$

其中：

$f_L(p^I,Q_L)=\frac{\sqrt{2}(p^I-p_L)}{\sqrt{{\rho }_L[{\gamma }_L(p_L+B_L)+{\gamma }_L(p^I+B_L)+p^I-p_L]}}---\left(50\right)$

$f_R(p^I,Q_R)=\frac{\sqrt{2}(p^I-p_R)}{\sqrt{{\rho }_R[{\gamma }_R(p_R+B_R)+{\gamma }_R(p^I+B_R)+p^I-p_R]}}---\left(51\right)$

由于液体的声速：

$a_R=\sqrt{k_R/{\rho }_R}=\sqrt{{\gamma }_R(p_R+B_R)/{\rho }_R}---\left(52\right)$

假设为无穷大，但是ρ_R为有限值，所以：

k_R＝γ_R(p_R+B_R)→+∞(53)

将其代入式(46)，(48)，(51)，(47)和(49)得到:

$\underset{k_R\to +\infty }{\lim }{\rho }_R^I={\rho }_R---\left(54\right)$

$\underset{k_R\to +\infty }{\lim }{W}_R=+\infty ---\left(55\right)$

f_R(p^I,Q_R)＝0(56)

$\frac{p^I-p_L}{W_L}+u_R-u_L=0---\left(57\right)$

$u^I=\frac{1}{2}(u_L+u_R)-\frac{1}{2}{f}_L(p^I,Q_L)---\left(58\right)$

对于左侧气体有：

u^I＝u_L-f_L(p^I,Q_L)(59)

可以得到：

u^I＝u_R(60)

对于式(57)通过迭代方法求解p^I。可由式(45)得到。在上述可压缩与不可压缩 黎曼问题的解中，有这一结果满足不可压缩介质密度保持不变的特点。u^I＝u_R表明 了多介质界面的速度应当由不可压缩介质决定。

(7)修正虚拟流体方法

在修正虚拟流体方法中，在界面附近通过构造黎曼问题来预测界面状态。由于在界 面追踪部分已经构造标志点处的黎曼问题，预测的界面状态可以直接用于定义虚拟流体状态。 如图2所示，可压缩气体(介质1)和不可压缩液体(介质2)被界面Π分开。是标志点的法向矢量，是界面附近可压缩气体中网格P的法向矢量。网格P的密度可由与 其法向夹角最小的标志点上的界面状态进行更新。如果标志点E是与网格P法向夹角最小的 标志点，那么对于标记点E通过等熵修正更新网格P密度的平均值。界面附近其它可压缩气 体的密度更新方法类似。虚拟流体状态通过求解如下方程得到：

其中表示密度，法向速度，切向速度和压力。如图2所示，在求解介质1的边界条 件时采用“+”号，否则采用“-”号。采用迭代方法求解该方程，空间上采用一阶迎风格式离 散，时间上采用三阶TVDRK方法进行离散。由于DG方法空间紧致性很好，此处仅采用2～3个 网格作为虚拟流体区域。通过把虚拟流体状态投影到基函数空间，得到虚拟网格状态的平 均值，这样即可确定界面边界条件。

(8)新方法验证

数值算例采用均匀结构网格。可压缩与不可压缩DG方法的多项式阶数均取为2。式 (7)中的CFL数C_f取为0.18。RKDG方法中的TVB限制器常量取为0.1。比热比γ取为1.4，粘 性系数μ＝0.001137kg/(ms)。数值算例中同时给出了基于newGFM的计算结果作为比较。 为了进行更为细致的比较，在压力和速度结果下方给出了计算结果的局部放大图。

一维算例

所有一维算例的计算域为[0m,1m]，网格数为200。如果气体和液体均作为可压缩 介质，算例8.1～8.3实际上为可压缩黎曼问题，可以得到其相应的解析解以作进一步比较。

算例8.1.气液界面位于x＝0.4m，初始条件为：气体：ρ＝1kg/m³，u＝1000m/s,p＝1 ×10⁵Pa，液体：ρ＝1009kg/m³，u＝0m/s，p＝2×10⁷Pa。图3给出了t＝2.8×10^-4s时压力和速 度的计算结果。观察可压缩黎曼问题的解析解，可以发现在液体中形成了稀疏波(但在速度 结果中不是很明显)。然而，当液体为不可压缩介质时，稀疏波立刻传出计算域。基于MGFM得 到的激波位置与可压缩解析解基本一致，而基于newGFM得到的激波在可压缩气体中传播 较慢。这是由于在MGFM中，通过求解可压缩与不可压缩黎曼问题预测界面状态。对于在可压 缩气体中产生激波的问题，界面状态事先被预测出来。这导致在计算的最初几个时间步中 产生强度较大的激波，所以基于MGFM得到的激波传播较快。

算例8.2.气液界面位于x＝0.4m，初始条件为：气体：ρ＝1kg/m³，u＝800m/s,p＝1 ×10⁵Pa，液体：ρ＝1000kg/m³，u＝0m/s，p＝1×10⁵Pa。图4给出了t＝2.8×10^-4s时压力和速 度的计算结果。可以发现在可压缩黎曼问题的解析解中，液体中形成了激波(但在速度结果 中不是很明显)。但是当液体为不可压缩介质时，该激波立刻以无穷大的速度传出计算域。 基于MGFM得到的激波位置与可压缩解析解基本一致，而基于newGFM得到的激波在可压缩 气体中传播较慢。

算例8.3.气液界面位于x＝0.5m，初始条件为：气体：ρ＝10kg/m³，u＝-1000m/s,p ＝1×10⁷Pa，液体：ρ＝1002kg/m³，u＝0m/s，p＝5×10⁶Pa。图5给出了t＝2×10^-4s时压力和速 度的计算结果。在可压缩黎曼问题的解析解中，稀疏波仍然在液体中传播(但在速度结果中 不是很明显)。可以发现，当在可压缩气体中出现稀疏波时，这些方法得到的稀疏波位置基 本一致。

算例8.4.液滴位于0.4<x<0.6的区域中。初始条件为：气体：ρ＝1.226kg/m³，u＝ 0m/s，p＝1×10⁵Pa，液滴：ρ＝1000kg/m³，u＝100m/s，p＝1×10⁵Pa。图6给出了t＝7.5×10^-4s 时压力和速度的计算结果。由于液滴向右运动，在右侧的气体中会形成激波，而在左侧的气 体中会形成稀疏波。基于MGFM和newGFM得到的稀疏波位置基本一致，而基于MGFM得到的激 波位置几乎比newGFM领先一个网格步长。

算例8.5.气体中位于x＝0.1m的右行激波打击位于0.4<x<0.6的液滴。初始条件 为：波前气体：ρ＝1.58317kg/m³，u＝0m/s,p＝98066.5Pa，波后气体：ρ＝2.124kg/m³，u＝ 89.981m/s,p＝148407.3Pa，液滴：ρ＝1000kg/m³，u＝0m/s，p＝98066.5Pa。图7给出了t＝ 1.75×10^-3s时压力和速度的计算结果。由于液滴的刚性较强，在左侧的气体中会反射激波， 而在液滴中形成速度很快、强度较弱的传递激波。可以发现基于MGFM得到的激波位置几乎 仍然比newGFM领先一个网格步长。

二维算例

二维多介质问题的计算域为[0m,1m]×[0m,1m]，网格数为100×100。半径为0.2m 的圆形液滴位于计算域的中央。由于流场上下对称，在二维算例中仅仅模拟下面半个流场 的流动问题。

算例8.6.初始条件为：气体：ρ＝1.226kg/m³，u＝0m/s，v＝0m/s，p＝1×10⁵Pa，液 滴：ρ＝1000kg/m³，u＝100m/s，v＝0m/s，p＝1×10⁵Pa。图8给出了t＝5×10^-4s时沿着y＝ 0.5m压力和x方向速度的计算结果。由于液滴向右运动，在右侧的气体中会形成激波，而在 左侧的气体中会形成稀疏波。基于MGFM得到的激波位置几乎比newGFM领先一个网格步长。 图9给出了相同时刻的压力云图(图中虚线表示界面位置)。可以看出由于液滴的刚性较大， 界面的形状基本保持不变。为了进一步考察液滴的密度对于计算结果的影响，这里将液滴 密度替换为ρ＝10kg/m³，压力和x方向速度的计算结果如图10所示。基于newGFM得到的激 波位置几乎落后一个网格步长。图11给出在t＝2.5×10^-3s、网格尺寸为0.02m，0.1m，0.005m 时的界面形状(图中虚线表示初始时刻界面形状)。可以看出较轻的液滴更易发生变形且速 度降低较快。密网格上已经出现了Kelvin-Helmholtz不稳定现象。这主要是因为在界面处 没有考虑粘性力的影响从而不能保证界面两侧切向速度连续。然而粗网格的数值粘性在一 定程度上弱化了这一现象。对上述不同网格尺寸下的液滴损失率进行研究，计算结果为： 0.0065，0.0032，0.0011。可以看出随着网格的不断细化，水滴的损失率逐渐变小。

算例8.7.气体中位于x＝0.1m的右行激波打击液滴。初始条件为：波前气体：ρ＝ 1.58317kg/m³，u＝0m/s，v＝0m/s，p＝98066.5Pa，波后气体：ρ＝2.124kg/m³，u＝89.981m/s， v＝0m/s，p＝148407.3Pa，液滴：ρ＝1000kg/m³，u＝0m/s，v＝0m/s，p＝98066.5Pa。图12给出 了t＝1.25×10^-3s时沿着y＝0.5m压力和x方向速度的计算结果。基于newGFM得到的激波位 置几乎落后一个网格步长。相同时刻的压力云图如图13所示(图中虚线表示界面位置)，可 以较清晰地看出入射激波和反射激波。由于液滴的刚性较强，液滴几乎保持静止。和算例 8.6相似，这里把液滴的密度替换为ρ＝10kg/m³，压力和速度的结果如图14所示。可以发现 基于MGFM得到的激波位置领先一个网格步长。图15给出了在t＝2.5×10^-3s、网格尺寸为 0.02m，0.1m，0.005m时的界面形状(图中虚线表示初始时刻界面形状)。通过观察发现， Kelvin-Helmholtz不稳定现象再次出现在细网格上。计算上述不同网格尺寸下的液滴损失 率：0.0175，0.0041，0.003，可以看出损失率随着网格的细化而变小。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人 员来说，在不脱离本发明原理的前提下还可以作出若干改进，这些改进也应视为本发明的 保护范围。