一种中心型间断有限元MMALE方法

卿 芳, 蔚喜军, 赵晓龙, 邹世俊, 贾祖朋,*

(1. 北京应用物理与计算数学研究所, 计算物理实验室, 北京 100088;2. 中国工程物理研究院研究生院, 北京 100088)

0 引 言

多介质可压缩流体动力学问题的数值计算方法主要分为两类：欧拉方法(Eulerian method)和拉氏方法(Lagrangian method)。欧拉方法采用固定的空间网格进行计算。拉氏方法的网格跟随流体一起运动。这两种方法各有利弊。欧拉方法易于处理流场大变形，但是难以精确地捕捉物质界面。拉氏方法能清晰地刻画、描述物质界面，但是当流场中发生大变形时，网格会严重扭曲，导致计算精度下降甚至计算终止。为了结合这两种方法的优点，Hirt等[1]提出了ALE(Arbitrary Lagrangian-Eulerian)方法，其中网格可以以任意指定的方式移动。大多数ALE算法[2-7]包括三步：拉氏步、重分步和重映步。传统的ALE方法用单元边描述物质界面，并且要求每个单元只有一种介质。当物质界面变形严重时，很难生成质量良好的新网格。为了解决这个问题，Peer等[8]将传统的ALE方法与界面捕捉方法相结合，提出了MMALE(Multi-material ALE)方法。这种方法引入了混合单元(一个单元中存在多种介质)，在混合单元中重构物质界面，重分网格时可以跨过界面生成新的网格，因此比传统的ALE方法提供了额外的灵活性。MMALE方法的计算流程与传统的ALE方法类似。首先，在拉氏步更新每种介质的物理量以及网格坐标，其中对混合单元需要用热力学封闭模型来更新各种介质的物理量。当网格变形严重时，需要重构混合单元的物质界面，为后面重映步提供界面信息。然后，在重分步生成一个几何品质较好的新网格。最后,在重映步将旧网格的物理量映射到新网格上。

在拉氏步，常用的离散方法有两种。一种是交错型格式[9-12]，指的是速度定义在节点上，而其他变量(密度、压力和内能等)定义在单元中心。另一种是中心型格式[13-17]，所有变量都定义在单元中心。中心型格式比交错型格式有一些优点。例如，交错型格式对变量使用不同的控制体,很难对所有变量构建一致的高阶精度。本文采用中心型格式。中心型格式可以分为两类: 中心型有限体积方法[16-17]和中心型间断有限元方法[18-20]。有限体积方法需要较大的单元模板来构造高阶格式；而间断有限元是利用单元内的高阶插值多项式来构造高阶格式。因此，间断有限元方法比有限体积方法更具有发展高阶格式的前景。重映步通常采用两种方法。一种是对流重映[21-25]，它要求新旧网格离得比较近，计算效率比较高.另一种是积分重映[26-31],它对新旧网格的关联要求不高，但是需要计算新旧网格的交点，计算结果比较精确。近年来，在MMALE方法的研究上取得了一些进展[32-35]。Galera和Maire等[33-34]提出了一种单元中心型二维MMAIE方法，界面重构采用VOF(Volume of Fluid)或MOF方法。贾祖朋等[29]提出了一种基于MOF界面重构的交错型非结构MMALE方法，模拟涉及强剪切变形的多介质可压缩流动。后来，贾祖朋等[36]发展了一种基于改进的MOF方法的中心型MMALE方法。上面提到的MMALE方法在拉氏步都是用有限体积方法。

本文提出了一种二维中心型间断有限元MMALE方法。在拉氏步，使用间断有限元方法对流体力学方程组进行离散，选取文献[19]中的基函数。为了构造二阶积分重映算法，对拉氏步获得的物理量单元均值做分片线性重构。经典的限制方法是对重构多项式的梯度进行限制[26,37]。近些年,Blanchard和Loubere[38]提出了基于MOOD (Multi-dimensional Optimal Order Detection)限制策略[39-40]的后验校正。本文采用后验校正，并做了一些改动以适应多介质计算。对新旧网格相交的计算，借鉴文献[31]的三维“剪裁投影”算法的基本思想，设计了二维“剪裁投影”算法，显著降低了多边形相交算法的复杂度。

1 MMALE方法的计算步骤

本文的MMALE方法的计算步骤如图1所示。首先是初始化步：定义初始网格上所有物理量的分布，使用文献[41]中2.3.1节的方法定义初始网格中物质体积分数和中心坐标。在拉氏步，采用单元中心型间断有限元方法求解流体动力学方程，并将物理量和网格从时间tn更新到时间tn+1=tn+Δt上。混合单元的封闭模型采用Tipton压力松弛模型[42-43]。使用文献[44]中等参坐标方法更新物质的中心坐标。界面重构采用卿芳等提出的一种健壮的MOF算法[45]。在网格重分步，利用Knupp算法[46]对拉氏网格进行光滑处理；或者采用初始网格作为重分网格。最后，在重映步，使用二阶积分守恒重映方法，将时间tn+1的拉氏网格中的所有物理量映射到重分后的新网格上。

图1 MMALE 方法流程图

2 间断有限元拉氏方法

欧拉框架下的二维无粘可压缩Euler方程组的向量形式为：

(1)

其中U=(ρ,ρu,ρv,ρE)T，F=(0,p,0,pu)T，G=(0,0,p,pv)T，∇·是散度算子，即∇·=∂/∂x+∂/∂y。ρ是流体密度,u=(u,v)是流体速度，p是压力，E是总能量，相应的比内能为e=E-u2/2，状态方程为p=p(ρ,e)。

下面给出式(1)的弱形式。式(1)两边同时乘以函数φ，并经过微分运算，得：

+∇·(φ(F,G))-(F,G)·∇φ=0

(2)

对上式第一项用Reynolds输运定理以及对第二项用散度定理，最终得到：

其中n=(nx,ny) 表示单元边界S的单位外法向量。

∀Wk,k=1,…,N}

(5)

其中Pn(Wk)为单元Wk上不超过n次多项式的全体集合。那么，对任意的单元Wk，下面的方程成立：

定义标准单元Ωst=[-1,1]2，它的坐标记为(ξ,η)，ξ,η∈[-1,1]。对于任意的时间t，存在一个从Ωst到Wk的双线性映射Fk:(ξ,η)∈Ωst→(x,y)∈Wk:

标准单元Ωst的四个基函数为σ1=1,σ2=ξ,σ3=η,σ4=ξη。由于Fk是光滑的可逆映射，因此：

i=1,…,4

(8)

数值解在单元的交界面上可以是间断的。因此在单元边界上定义一个数值通量函数F*和G*来代替原来的通量函数F和G[20]。边[Mr,Mr+1]上的数值通量定义为：

图2 单元的记号

将(13)-(16)代入(12)：

m=1,2,3,4

(18)

m=1,2,3,4

(19)

最后，如果令：

uk1=[ρk1,ρk2,ρk3,ρk4]

(22)

uk2=[ρuk1,ρuk2,ρuk3,ρuk4]

(23)

uk3=[ρvk1,ρvk2,ρvk3,ρvk4]

(24)

uk4=[ρEk1,ρEk2,ρEk3,ρEk4]

(25)

其中Lkj(j=1,2,3,4)是(17)～(20)式空间离散算子得到的向量。

半离散格式(26)的时间离散采用二阶显式Runge-Kutta方法[48]。时间步长的估计遵循标准的CFL准则和每个时间步上限制单元面积的变化，详情请参考文献[16-17]。空间二阶重构采用最小二乘算法[17]，并采用Barth-Jesperson限制器[49]抑制非物理的数值振荡。

3 界面重构

图3 MOF方法示意图

其中xh(θ)表示直线截取的多边形的中心点。文献[50]用最优化方法求解上述极小化问题，有可能收敛到局部极小值。因此，近年有学者对经典的MOF方法做了一些改进。

注意到f(θ)的一阶导数[50-51]，

文献[51]中提到f′(θ)是关于θ的连续函数。因此，f(θ)的最小值点一定是f′(θ)的零点或者区间端点。文献[45]提出了一种健壮的算法计算f′(θ)的零点。本文采用这种健壮的MOF方法[45]。

4 重映

4.1 多项式重构

4.2 新旧网格相交

图4 新单元和拉氏网格{W}的相交部分。拉氏网格和新网格分别用蓝线和黑线表示。所要求的相交部分分别用不同颜色(a-f)标记

图5 单元内物质用物质界面(红线)分隔开。单介质子多边形用Pgreen，Pcyan表示

下面介绍求相交子多边形的一个算法。这个算法遵循陈享等[31]的三维“剪裁投影”算法基本思想，并改为退化的二维“剪裁投影”算法。由于新单元是四边形，可以先把它分解成两个三角形。然后用这两个三角形分别和旧网格求相交子多边形，因而，新单元的面积等于这两个三角形分别得到的相交子多边形的面积之和。这个新旧网格相交的算法归结为求三角形与多边形的交点。算法分为以下三步：

1) 仿射变换

2) “剪裁投影”算法

图6 单位三角形和多边形的交点

3) 相交多边形的面积和中心坐标

最后，通过坐标变换，将参考坐标系中相交多边形的面积和中心坐标转换为物理坐标系中面积和中心坐标：

4.3 积分

4.4 后验校正

一般而言，重构多项式(28)需要一个限制器来抑制间断附近的振荡。本文采用文献[38]中提出的基于MOOD限制策略的后验校正。后验校正的原理是检测出“坏”单元，然后降低“坏”单元上重构多项式的次数，重复之前的步骤，直到这个单元是“好”单元为止。值得注意的是，文献[38]中提出的方法只适用于单介质流。为了使其适应于多介质流，我们对它做了一些小改动。对于纯(单介质)新单元，对新单元进行检测；对于混合(多介质)新单元，对所有单介质子多边形进行检测，而不是整个新单元。

检测的目的是过滤“好”单元，这些单元已经使用无限制的重构多项式进行精确积分守恒重映；而把存在一些问题(振荡，非物理等)的“坏”单元标识出来。检测准则有两个：

1) 物理相容性检验准则(PAD)

该准则的依据是密度、内能和压强的必须为正：

2) 数值相容性检验准则(NAD)

该准则是为了检验数值解是否局部光滑：

DNAD=

图7 基于后验校正的重映算法示意图

5 数值算例

5.1 二维周期涡问题[36]

表1 误差及收敛阶(拉氏方法)

表2 误差及收敛阶(MMALE方法)

5.2 圆柱内爆中的Rayleigh-Taylor不稳定性问题[53]

这是一个类似ICF的测试问题。初始状态如图8所示，一团轻气体(R∈[0,1])被一团重气体包围(R∈[0,1.2]))。两种气体均为理想气体(绝热指数γl=γh=1.5)。轻气体的初始状态为(ρl,pl,Ul)=(0.05,0.1,0),重气体为(ρh,ph,Uh)=(1,0.1,0)。在重气体外表面给定随时间变化的压力边界条件，

图8 初始区域的几何数据

图9 初始网格

对于MMALE方法，在开始时初始网格没有混合单元，因此仅需执行拉氏步。当网格由于界面不稳定而发生严重变形时，需要进行重分和重映。如果网格中存在混合单元，采用Tipton压力松弛模型进行计算，并在重映步之前对这些混合单元执行MOF界面重构算法。

当取a=10-3时，图10给出t=0.3时用拉氏方法计算的网格和密度图。注意到最终的网格质量保持很好，没有翻转单元。图11给出的是t=0.3时用MMALE方法计算网格与界面和密度等值线图。从图10和图11可见，拉氏方法和MMALE方法的计算结果是非常相似的。

图10 a=2×10-4，t=0.3时用间断有限元拉氏方法计算的网格(左)和密度图(右)

图11 a=2×10-4，t=0.3时用MMALE方法计算的网格(左)和密度(右)

当取a=10-3时，图12给出了t=0.25时用拉氏方法计算的网格和局部放大图，此时网格是没有翻转单元的。当t=0.3时，用拉氏方法计算的结果是不准确的，此时网格扭曲严重而且存在翻转单元，如图13所示。然而，用MMALE方法运行到最终时刻没有任何问题。图14给出的是t=0.3时用MMALE方法计算的网格与界面和密度等值线图。该算例的计算结果表明了MMALE方法具有较好的灵活性和鲁棒性。

图12 a=10-3，t=0.25时用间断有限元拉氏方法计算的网格(左)和局部放大图(右)

图13 a=10-3，t=0.3时用间断有限元拉氏方法计算的网格(左)和局部放大图(右)

图14 a=10-3，t=0.3时用MMALE方法计算的网格(左)和密度(右)

5.3 Sedov 问题[54]

这是一个点爆炸问题。在原点处给定一个能量。随着时间的推进，一个圆柱型激波会向外传播。初始密度为ρ=1，速度为(u,v)=(0,0)。在原点处，存在单位内能，并且在其它点处初始压力为p=0。这个问题的精确解由Sedov在文献[54]中给出：在时间t=1时，激波运动到离原点距离为1处，密度峰值达到6。使用本文的中心型MMALE方法进行计算。计算区域为Ω=[0,1.125]×[0,1.125]，初始剖分为45×45的一致网格。在靠近原点的第一个单元上给定内能e0=400。在初始时刻，界面位于以原点为圆心，半径为0.43的圆上。注意到在混合单元中，这两种介质都是理想气体，具有相同的绝热指数γ=1.4。我们将它们视为混合单元，以便将数值解与精确解进行比较。图15 给出了t=0.5,0,75,1时的网格和界面。图16是t=1时刻的密度等值线图和以半径为变量的密度图。从图16可知。计算得到的激波位置是准确的，数值解与精确解很接近。

图15 t=0.5,0,75,1时的网格和界面(从左到右)

图16 t=1时的密度等值线图(左)和以半径为变量的密度图(右)

5.4 二维Rayleigh-Taylor 不稳定性问题[6，36]

计算区域为矩形[0,1/3]×[0,1].初始时刻将它剖分为35×100的均匀网格。初始时刻，密度为ρh=2的重气体和密度为ρl=1的轻气体被扰动界面分隔开，界面方程为yi(x)=0.5+0.01cos(6πx)。重气体在轻气体之上。这两种气体都有相同的绝热指数γ=1.4。重力场垂直向下，重力加速度为g=(gx,gy)T=(0,-0.1)T。初始时刻两种气体都处于静止状态。初始压力分布通过静压平衡条件推导出来:

众所周知，这种状态是不稳定的。由于正弦界面的存在，在界面附近会产生旋涡。随着时间的推移，较重的气体会落入较轻的气体中形成一个尖钉；较轻的气体会上升到较重的气体中形成一个气泡。图17给出了t=5,7,9,11,13时的网格和界面。图18给出了t=5,7,9,11,13时的密度等值线图。从图17和图18可以看出，我们的结果与文献[6,29]中的结果非常接近，这里的计算结果对称性好，界面清晰。

图17 t=5,7,9,11,13时的网格和界面(从左到右)

图18 t=5,7,9,11,13时的密度等值线图(从左到右)

5.5 激波与氦气泡相互作用问题[6，36]

计算区域为[0,0.65]×[-0.089,0.089]，初始时刻将它剖分成260×72个均匀单元。如图19所示，气泡是由圆心(xc,yc)=(0.32,0)和半径r=0.025定义的圆盘。在右边界(初始时刻位于x=0.65)上给出了一个类似活塞的边界条件，它以速度(u*,0)向左移动。在其他边界上给固壁边界条件。入射激波的马赫数为Ms=1.22。氦气泡和空气在初始时刻处于静止状态。氦气泡和空气的初始数据分别为(ρ1,p1)=(0.182,105),(ρ2,p2)=(1,105)；摩尔质量和绝热指数分别为(M1,γ1)=(5.269×10-3,1.648)，(M2,γ2)=(28.963×10-3,1.4)。由Rankine-Hugoniot关系可以得到活塞在x方向的速度u*=-124.824。入射激波的x方向速度是Dc=-456.482。入射激波在ti=668.153×10-6时刻撞击氦气泡。计算终止时间tfinal=ti+674×10-6=1342.153×10-6。图20给出了t1=ti+245×10-6=913.153×10-6，t2=ti+4247×10-6= 1095.153×10-6和tfinal三个时刻的网格和界面以及这三个时刻的试验纹影图[55]。图21 给出了tfinal时刻的密度等值线图以及局部放大图。从图20和图21可以看出，这里的计算结果与文献[55]中的试验纹影图符合的很好，与文献[6，36]的数值结果非常接近。

图19 初始区域的几何数据

图20 t1，t2, tfinal时刻(从左到右)用MMALE方法计算的密度(上)和试验纹影图(下)

图21 tfinal时刻的密度等值线图(左) 和局部放大图(右)

5.6 三重点问题[6]

考虑如图22所示矩形域中的三种状态的二维黎曼问题。初始区域Ω=[0,7]×[0,3]分为三个子区域Ω1=[0,1]×[0,3]；Ω2=[1,7]×[0,1.5]和Ω3=[1,7]×[1.5,3]。Ω1中是初始状态为(ρ1,p1,U1)=(1,1,0)的高压高密度气体。Ω2中是初始状态为(ρ2,p2,U2)=(1,0.1,0)的低压高密度气体。Ω3中是初始状态为(ρ3,p3,U3)=(0.125,0.1,0)的低压低密度气体。Ω1和Ω3中物质的绝热指数均为γ1=γ3=1.5；Ω2中物质的绝热指数为γ2=1.4。边界条件均为固壁边界条件,时间计算到tfinal=10。初始网格数为210×90。由于声阻抗的不同，Ω2和Ω3中的两个激波以不同的速度传播。这在Ω2和Ω3的界面处产生强剪切力。这种剪切作用产生Kelvin-Helmholtz不稳定性，并形成涡旋。用拉氏方法计算该问题时，在涡旋发展前由于网格扭曲缠结而导致计算终止。当使用经典的ALE方法时，捕获涡流是困难的。我们使用新的中心型MMALE方法计算此问题。图23 给出了t=5,6时的网格和界面。

图22 初始区域的几何数据

图23 t=5,6的网格和界面(从上到下)

6 结 论

本文提出了一种具有二阶精度的中心型间断有限元MMALE方法。在该方法中，拉氏步采用间断有限元方法求解可压缩欧拉方程。对于混合网格，采用Tipton压力松弛模型更新物理量。采用具有鲁棒性的MOF方法进行界面重构。针对重映步的多边形相交问题，提出了一种“裁剪投影”算法；在此基础上构建了一个基于后验校正的二阶积分守恒重映方法。数值试验的结果表明，该方法具有较好的鲁棒性和灵活性。