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

2 间断有限元拉氏方法

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

$ \frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \nabla \cdot (\mathit{\boldsymbol{Uu}}) + \nabla \cdot (\mathit{\boldsymbol{F}},\mathit{\boldsymbol{G}}) = 0 $

(1)

其中$\mathit{\boldsymbol{U}} = {(\rho , \rho u, \rho v, \rho E)^{\rm{T}}}$，$\mathit{\boldsymbol{F}} = {(0, p, 0, pu)^{\rm{T}}}$，$\mathit{\boldsymbol{G}} = {(0, 0, p, pv)^{\rm{T}}}$，$\nabla \cdot $是散度算子，即$\nabla \cdot = \partial /\partial x + \partial /\partial y$。ρ是流体密度, u=(u, v)是流体速度，p是压力，E是总能量，相应的比内能为e=E-u²/2，状态方程为p=p(ρ, e)。

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

$ \begin{array}{l} \frac{{\partial (\varphi \mathit{\boldsymbol{U}})}}{{\partial t}} + \nabla \cdot (\varphi \mathit{\boldsymbol{Uu}}) - U\left( {\frac{{\partial \varphi }}{{\partial t}} + \mathit{\boldsymbol{u}} \cdot \nabla \varphi } \right)\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} + \nabla \cdot (\varphi (\mathit{\boldsymbol{F}},\mathit{\boldsymbol{G}})) - (\mathit{\boldsymbol{F}},\mathit{\boldsymbol{G}}) \cdot \nabla \varphi = 0 \end{array} $

(2)

由于$\frac{{{\rm{d}}\varphi }}{{{\rm{d}}t}} = \frac{{\partial \varphi }}{{\partial t}} + \mathit{\boldsymbol{u}} \cdot \nabla \varphi $，并在Ω(t)上积分，得：

$ \begin{array}{*{20}{l}} {\int_{\mathit{\Omega }(t)} {\frac{{\partial (\varphi \mathit{\boldsymbol{U}})}}{{\partial t}}} {\rm{d}}\mathit{\Omega } + \int_{\mathit{\Omega }(t)} \nabla \cdot (\varphi \mathit{\boldsymbol{Uu}} + \varphi (\mathit{\boldsymbol{F}},\mathit{\boldsymbol{G}})){\rm{d}}\mathit{\Omega }}\\ {\quad - \int_{\mathit{\Omega }(t)} \mathit{\boldsymbol{U}} \frac{{{\rm{d}}\varphi }}{{{\rm{d}}t}}{\rm{d}}\mathit{\Omega } - \int_{\mathit{\Omega }(t)} {(\mathit{\boldsymbol{F}},\mathit{\boldsymbol{G}})} \nabla \varphi {\rm{d}}\mathit{\Omega } = 0} \end{array} $

(3)

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

$ \begin{array}{c} \frac{{\rm{d}}}{{{\rm{d}}t}}\int_{\mathit{\Omega }(t)} \varphi \mathit{\boldsymbol{U}}{\rm{d}}\mathit{\Omega } = \int_{\mathit{\Omega }(t)} \mathit{\boldsymbol{U}} \frac{{{\rm{d}}\varphi }}{{{\rm{d}}t}}{\rm{d}}\mathit{\Omega } + \int_{\mathit{\Omega }(t)} {(\mathit{\boldsymbol{F}},\mathit{\boldsymbol{G}})} \cdot \nabla \varphi {\rm{d}}\mathit{\Omega }\\ - \int_{\partial \mathit{\Omega }(t)} {\left( {\varphi \mathit{\boldsymbol{F}}{n^x} + \varphi \mathit{\boldsymbol{G}}{n^y}} \right)} {\rm{d}}S \end{array} $

(4)

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

接下来对空间进行离散。假设初始时刻的计算区域为W_t=0=Ω₀，Ω₀被剖分成N个不重叠的四边形单元{W_k⁰，k=1, …, N}。在任意时刻t，这些四边形单元记为{W_k，k=1, …, N}。定义一个分片多项式有限元空间：

$ \begin{array}{c} V_h^n = \left\{ {{\varphi _h} \in {L_2}(\mathit{\Omega }(t)):{\varphi _h}\mid {W_k} \in {P_n}\left( {{W_k}} \right),} \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {\forall {W_k},k = 1, \cdots ,N} \right\} \end{array} $

(5)

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

$ \begin{array}{*{20}{c}} {\frac{{\rm{d}}}{{{\rm{d}}t}}\int_{{W_k}} {{\varphi _h}} {\mathit{\boldsymbol{U}}_h}{\rm{d}}{W_k} = \int_{{W_k}} {{\mathit{\boldsymbol{U}}_h}} \frac{{{\rm{d}}{\varphi _h}}}{{{\rm{d}}t}}{\rm{d}}{W_k} + }\\ {\int_{{W_k}} {\left( {{\mathit{\boldsymbol{F}}_h},{\mathit{\boldsymbol{G}}_h}} \right)} \cdot \nabla {\varphi _h}{\rm{d}}{W_k}}\\ { - \int_{\partial {W_k}} {\left( {{\varphi _h}{\mathit{\boldsymbol{F}}_h}{n^x} + {\varphi _h}{\mathit{\boldsymbol{G}}_h}{n^y}} \right)} {\rm{d}}l} \end{array} $

(6)

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

$ x = \sum\limits_{i = 1}^4 {{x_i}} {\psi _i}(\xi ,\eta ),\quad y = \sum\limits_{i = 1}^4 {{y_i}} {\psi _i}(\xi ,\eta ) $

(7)

其中(x_i, y_i)表示单元W_k的顶点M_i的坐标，${\psi _i} = \frac{1}{4}(1 \pm \xi )(1 \pm \eta ), i = 1, \cdots , 4$是双线性形函数。

标准单元Ω_st的四个基函数为σ₁=1, σ₂=ξ, σ₃=η, σ₄=ξη。由于F_k是光滑的可逆映射，因此：

$ \begin{array}{*{20}{c}} {{\sigma _i}(\xi ,\eta ) = {\sigma _i}\left( {F_k^{ - 1}(x,y)} \right) \buildrel \Delta \over = {B_i}(x,y,t),}\\ {i = 1, \cdots ,4} \end{array} $

(8)

选取B_i, i=1, …, 4为有限元空间V_hⁿ的基函数。这个基函数有个非常重要的性质，它的物质导数为零，即$\frac{{{\rm{d}}{B_i}}}{{{\rm{d}}t}} = 0$，参见文献[47]。U_h和φ_h可表示为：

$ {{\mathit{\boldsymbol{U}}_h} = \sum\limits_{i = 1}^4 {{\mathit{\boldsymbol{U}}_{ki}}} {B_{ki}},\quad (x,y) \in {W_k}} $

(9)

$ {{\varphi _h} = \sum\limits_{i = 1}^4 {{\varphi _{ki}}} {B_{ki}},\quad (x,y) \in {W_k}} $

(10)

用间断有限元方法对(1)进行空间离散，就是寻找数值解U_h=(ρ_h, (ρu)_h, (ρv)_h, (ρE)_h)^T∈V_hⁿ，使得对所有的试探函数φ_h∈V_hⁿ和所有单元W_k, k=1, …N，满足：

$ \begin{array}{*{20}{l}} {\frac{{\rm{d}}}{{{\rm{d}}t}}\int_{{W_k}} {{\varphi _h}} {\mathit{\boldsymbol{U}}_h}{\rm{d}}{W_k} = \int_{{W_k}} {\left( {{\mathit{\boldsymbol{F}}_h},{\mathit{\boldsymbol{G}}_h}} \right)} \cdot \nabla {\varphi _h}{\rm{d}}{W_k}}\\ { - \int_{\partial {W_k}} {\left( {{\varphi _h}{\mathit{\boldsymbol{F}}_h}{n^x} + {\varphi _h}{\mathit{\boldsymbol{G}}_h}{n^y}} \right)} {\rm{d}}l} \end{array} $

(11)

由于对所有试探函数φ_h∈V_hⁿ，上式都成立。因此上式等价于：

$ \begin{array}{*{20}{c}} {\frac{{\rm{d}}}{{{\rm{d}}t}}\sum\limits_{i = 1}^4 {\int_{{W_k}} {{\mathit{\boldsymbol{U}}_{ki}}} } {B_{ki}}{B_{km}}{\rm{d}}{W_k} = \int_{{W_k}} {{\mathit{\boldsymbol{F}}_h}} \cdot {\partial _x}{B_{km}}{\rm{d}}{W_k} + }\\ {\int_{{W_k}} {{\mathit{\boldsymbol{G}}_h}} \cdot {\partial _y}{B_{km}}{\rm{d}}{W_k}}\\ { - \int_{\partial {W_k}} {\left( {{\mathit{\boldsymbol{F}}_h}{n^x} + {\mathit{\boldsymbol{G}}_h}{n^y}} \right)} {B_{km}}{\rm{d}}l,m = 1,2,3,4} \end{array} $

(12)

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

$ \left\{ {\begin{array}{*{20}{l}} {F_{r,r + 1}^{1,*} = 0}\\ {F_{r,r + 1}^{2,*} = p_{r,r + 1}^*}\\ {F_{r,r + 1}^{3,*} = 0}\\ {F_{r,r + 1}^{4,*} = (pu)_{r,r + 1}^*} \end{array}\left\{ {\begin{array}{*{20}{l}} {G_{r,r + 1}^{1,*} = 0}\\ {G_{r,r + 1}^{2,*} = 0}\\ {G_{r,r + 1}^{3,*} = p_{r,r + 1}^*}\\ {G_{r,r + 1}^{4,*} = (pv)_{r,r + 1}^*} \end{array}} \right.} \right. $

(13)

引入与文献[16]相同的记号。如图 2所示，在每条边上引入四个压力，其中每侧引入两个压力。$p_{r, r + 1/2}^{*, i}\left( {p_{r + 1/2, r + 1}^{*, i}} \right)$表示“半边”[M_r, M_r+1/2]([M_r+1/2, M_r+1])在单元W_k的压力。V_r^*=(u_r^*, v_r^*)表示在顶点M_r的速度。令：

$ {\mathit{\boldsymbol{V}}_{r,r + 1}^* = \frac{1}{2}\left( {\mathit{\boldsymbol{V}}_r^* + \mathit{\boldsymbol{V}}_{r + 1}^*} \right)} $

(14)

$ {p_{r,r + 1}^{*,k} = \frac{1}{2}\left( {p_{r,r + \frac{1}{2}}^{*,k} + p_{r + \frac{1}{2},r + 1}^{*,k}} \right)} $

(15)

$ {(p\mathit{\boldsymbol{V}})_{r,r + 1}^{*,k} = \frac{1}{2}\left( {p_{r,r + \frac{1}{2}}^{*,k}\mathit{\boldsymbol{V}}_r^* + p_{r + \frac{1}{2},r + 1}^{*,k}\mathit{\boldsymbol{V}}_{r + 1}^*} \right)} $

(16)

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

$ {\frac{{\rm{d}}}{{{\rm{d}}t}}\sum\limits_{i = 1}^4 {\int_{W_{k}} {{\rho _{ki}}} } {B_{ki}}{B_{km}}{\rm{d}}{W_k} = 0,m = 1,2,3,4} $

(17)

$ \begin{array}{*{20}{c}} {\frac{{\rm{d}}}{{{\rm{d}}t}}\sum\limits_{i = 1}^4 {\int_{W_{k}} {{{(\rho u)}_{ki}}} } {B_{ki}}{B_{km}}{\rm{d}}{W_k} = \int_{W_{k}} {{p_k}} {\partial _x}{B_{km}}{\rm{d}}{W_k}}\\ { - \sum\limits_{r = 1}^4 {\frac{1}{2}} \left( {p_{r,r + \frac{1}{2}}^{*,k} + p_{r + \frac{1}{2},r + 1}^{*,k}} \right)n_{r,r + 1}^x\int_{r,r + 1} {{B_{km}}} {\rm{d}}l,}\\ {m = 1,2,3,4} \end{array} $

(18)

$ \begin{array}{*{20}{c}} {\frac{{\rm{d}}}{{{\rm{d}}t}}\sum\limits_{i = 1}^4 {\int_{W_{k}} {{{(\rho v)}_{ki}}} } {B_{ki}}{B_{km}}{\rm{d}}{W_k} = \int_{W_{k}} {{p_k}} {\partial _y}{B_{km}}{\rm{d}}{W_k}}\\ { - \sum\limits_{r = 1}^4 {\frac{1}{2}} \left( {p_{r,r + \frac{1}{2}}^{*,k} + p_{r + \frac{1}{2},r + 1}^{*,k}} \right)n_{r,r + 1}^y\int_{r,r + 1} {{B_{km}}} {\rm{d}}l,}\\ {m = 1,2,3,4} \end{array} $

(19)

$ \begin{array}{c} \frac{\mathrm{d}}{\mathrm{d} t} \sum\limits_{i=1}^{4} \int_{W_{k}}(\rho E)_{k i} B_{k i} B_{k m} \mathrm{~d} W_{k}= \\ \int_{W_{k}}(p u)_{k} \partial_{x} B_{k m}+(p v)_{k} \partial_{y} B_{k m} \mathrm{~d} W_{k} \\ -\sum\limits_{r=1}^{4} \frac{1}{2}\left[\left(p_{r, r+\frac{1}{2}}^{*, k} u_{r}^{*}+p_{r+\frac{1}{2}, r+1}^{*, k} u_{r+1}^{*}\right) n_{r, r+1}^{x}+\right. \\ \left.\left(p_{r, r+\frac{1}{2}}^{*, k_{1}} v_{r}^{*}+p_{r+\frac{1}{2}, r+1}^{*, k} v_{r+1}^{*}\right) n_{r, r+1}^{y}\right] \\ \int_{r, r+1} B_{k m} \mathrm{~d} l, \quad m=1,2,3,4 \end{array} $

(20)

关于单元W_k的顶点速度V_q^*和压力p_{q, k}^{*, k}(p_{q, k}^{*, k-1})的计算，参考Maire的节点求解器^[16-17]。

最后，如果令：

$ {{\mathit{\boldsymbol{M}}_k} = \int_{{W_k}} {{B_{ki}}} {B_{km}}{\rm{d}}{W_k}(i = 1,2,3,4;m = 1,2,3,4)} $

(21)

$ {{\mathit{\boldsymbol{u}}_{k1}} = \left[ {{\rho _{k1}},{\rho _{k2}},{\rho _{k3}},{\rho _{k4}}} \right]} $

(22)

$ {{\mathit{\boldsymbol{u}}_{k2}} = \left[ {\rho {u_{k1}},\rho {u_{k2}},\rho {u_{k3}},\rho {u_{k4}}} \right]} $

(23)

$ {{\mathit{\boldsymbol{u}}_{k3}} = \left[ {\rho {v_{k1}},\rho {v_{k2}},\rho {v_{k3}},\rho {v_{k4}}} \right]} $

(24)

$ {{\mathit{\boldsymbol{u}}_{k4}} = \left[ {\rho {E_{k1}},\rho {E_{k2}},\rho {E_{k3}},\rho {E_{k4}}} \right]} $

(25)

将顶点速度(u_r^*, v_r^*)和压力p_{q, k}^{*, k}(p_{q, k}^{*, k-1})代入(17)~(20)，可以得到：

$ \frac{{\rm{d}}}{{{\rm{d}}t}}{\mathit{\boldsymbol{M}}_k} \cdot {\mathit{\boldsymbol{u}}_{kj}} = {L_{kj}},\quad j = 1,2,3,4 $

(26)

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

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

3 界面重构

在拉氏步，需要对混合网格中的每种介质更新物质中心点坐标以及体积分数，以便在重映步之前用MOF方法重构物质界面。先来回顾一下经典的MOF方法^[50]。如图 3所示，给定混合网格内的某种介质物质中心点坐标x_c^*和体积分数f^*。n表示待确定的直线界面的单位法向量, θ表示向量n与x轴的夹角，即n=(cosθ, sinθ)。需要重构一条直线，使得该直线所截取的多边形的体积分数为f^*，而中心点与x_c^*尽可能的接近。求重构直线的问题可以转化为极小化目标函数:

$ f(\theta ) = \left\| {{x_h}(\theta ) - x_c^*} \right\|_2^2 $

(27)

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

注意到f(θ)的一阶导数^[50-51]，f′(θ)=2(x_h(θ)-x_c^*)·x′_h(θ), 0≤θ≤2π文献[51]中提到f′(θ)是关于θ 的连续函数。因此，f(θ)的最小值点一定是f′(θ)的零点或者区间端点。文献[45]提出了一种健壮的算法计算f′(θ)的零点。本文采用这种健壮的MOF方法^[45]。

4 重映

为了便于描述，引入一些记号。假设计算区域划分成两套不重叠的网格。拉氏(旧)网格{W}由单元W组成；重分(新)网格$\{ \tilde W\} $由单元${\tilde W}$组成。C′(W)表示单元W的相邻单元的集合。C(W)表示包含单元W的相邻单元的集合，即C(W)=C′(W)∪W。本文采用积分守恒重映方法，计算分为四步：多项式重构，新旧网格相交，积分和后验校正。

4.1 多项式重构

为了构造二阶重映算法，需要在旧单元上对每种介质的物理量g=ρ, ρu, ρe重构一个分片线性多项式。如果旧单元W_k是一个只含有介质i的纯网格单元，那么把它记为P_{i, k}；如果旧单元W_k是一个混合网格，需要用MOF方法重构物质界面，那么包含介质i的多边形记为P_{i, k}。设(x_{c, i, k}, y_{c, i, k})是多边形P_{i, k}的中心坐标。物理量g在P_{i, k}的均值用g_{i, k}表示.在P_{i, k}重构一个多项式g_{i, k}(x, y)：

$ \begin{array}{l} {g_{i,k}}(x,y) = {{\bar g}_{i,k}} + {\left( {\frac{{\partial g}}{{\partial x}}} \right)_{i,k}}\left( {x - {x_{c,i,k}}} \right) + \\ {\left( {\frac{{\partial g}}{{\partial y}}} \right)_{i,k}}\left( {y - {y_{c,i,k}}} \right) \end{array} $

(28)

其中梯度${\left( {\frac{{\partial g}}{{\partial x}}} \right)_{i, k}}$、${\left( {\frac{{\partial g}}{{\partial y}}} \right)_{i, k}}$用最小二乘法计算^[52]。

4.2 新旧网格相交

由于本文采用积分守恒重映方法，因此需要计算新单元${\tilde W}$和拉氏网格{W}相交的交点。如图 4所示，新单元${\tilde W}$由拉氏网格{W}与之相交的若干子多边形组成，即：

$ \tilde{W}=\bigcup\limits_{i, k} \tilde{W} \cap P_{i, k}, \tilde{W} \cap P_{i, k} \neq \varnothing $

(29)

因此，新单元${\tilde W}$的面积等于新单元${\tilde W}$与拉氏网格{W}相交子多边形的面积之和：

$ \tilde{V}_{\widetilde{W}}=\sum\limits_{i, k} V_{\widetilde{W} \cap P_{i, k}} $

(30)

值得注意的是，每一个相交子多边形只包含一种介质。因此，如果旧单元是混合网格，那么需要计算新单元与混合网格中单介质子多边形的交点，而不是新单元与整个混合网格的交点。如图 5，3×3的网格由两种介质组成，物质界面用中间的粗实线来表示。界面的左边是第一种物质，界面的右边是第二种物质。正中间的新单元${\tilde W}$由标记为a-h的子多边形组成。

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

1) 仿射变换

设参考坐标系为X，物理坐标系为x，参考坐标系中单位三角形${\mathbb{T}} $的顶点分别为O(0, 0)，E(1, 0)，F(0, 1)。通过仿射变换X=Ax+b，将物理空间的三角形T转换成参考坐标系中的单位三角形${\mathbb{T}} $。在此仿射变换下，物理空间的多边形P转换成参考坐标系中的多边形$ {\mathbb{P}}$。

2) “剪裁投影”算法

我们提出了一种“剪裁投影”算法，高效、鲁棒地计算参考坐标系中单位三角形$ {\mathbb{T}}$和多边形$ {\mathbb{P}}$的交点。在该算法中，首先多边形$ {\mathbb{P}}$的每条边都被直线X=0, Y=0和X=1依次切割。然后，将不在单位三角形$ {\mathbb{T}}$内部或者上面的部分丢弃。最后，将单位三角形$ {\mathbb{T}}$上方的外部部分投影到直线X轴上，得到最终的相交多边形$ {\mathbb{T}}$∩$ {\mathbb{P}}$。

以图 6(a)为例, 单位三角形$ {\mathbb{T}}$是OEF，多边形$ {\mathbb{P}}$是ABCD。图 6(b)-6(e)分别分析了多边形ABCD的每一条边与单位三角形相交。在图 6(b)中，标记了两个关键线段AP和PQ。线段AP是AB在单位三角形内部的部分，线段PQ是AB从+Y方向到单位三角形边EF上的投影。线段AP和PQ到X轴投影所组成的多边形记为S₁。同理，多边形ABCD另外三条边的关键线段如图 6(c)~6(e)所示，到X轴的投影所组成的多边形分别为S₂，S₃，S₄。因此, 单位三角形$ {\mathbb{T}}$和多边形$ {\mathbb{P}}$的相交多边形面积S=$ {\mathbb{T}}$∩$ {\mathbb{P}}$为$\int_s {\rm{d}} X{\rm{d}}Y = - \int_{s1} {\rm{d}} X{\rm{d}}Y + \int_{s2} {\rm{d}} X{\rm{d}}Y + \int_{s3} {\rm{d}} X{\rm{d}}Y - \int_{s4} {\rm{d}} X{\rm{d}}Y$使用多边形每条边的单位外法向矢量的Y分量来表示上述等式右边积分号前的正负号，它可以写成：

$ \begin{array}{*{20}{c}} {\int_S {{\rm{d}}} X{\rm{d}}Y = {\mathop{\rm sign}\nolimits} \left( {{n_{1Y}}} \right)\int_{{S_1}} {{\rm{d}}} X{\rm{d}}Y + {\mathop{\rm sign}\nolimits} \left( {{n_{2Y}}} \right)\int_{{S_2}} {{\rm{d}}} X{\rm{d}}Y}\\ {\quad + {\mathop{\rm sign}\nolimits} \left( {{n_{3Y}}} \right)\int_{{S_3}} {{\rm{d}}} X{\rm{d}}Y + {\mathop{\rm sign}\nolimits} \left( {{n_{4Y}}} \right)\int_{{S_4}} {{\rm{d}}} X{\rm{d}}Y = }\\ {\sum\limits_i {{\mathop{\rm sign}\nolimits} } \left( {{n_{iY}}} \right)\int_{{S_i}} {{\rm{d}}} X{\rm{d}}Y} \end{array} $

(31)

相交多边形${\mathbb{T}}$∩${\mathbb{P}}$的面积为：

$ V(\mathbb{T} \cap \mathbb{P})=\sum\limits_{i} \operatorname{sign}\left(n_{i Y}\right) \int_{S_{i}} \mathrm{~d} X \mathrm{~d} Y $

(32)

同理，可以得到${\mathbb{T}}$∩${\mathbb{P}}$的中心坐标：

$ X(\mathbb{T} \cap \mathbb{P})=\frac{\sum\limits_{i} \operatorname{sign}\left(n_{i Y}\right) \int_{S_{i}} \mathrm{~d} X \mathrm{~d} Y}{V(\mathbb{T} \cap \mathbb{P})} $

(33)

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

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

$ V(\boldsymbol{T} \cap \boldsymbol{P})=\operatorname{det}(\boldsymbol{A}) \times V(\mathbb{T} \cap \mathbb{P}) $

(34)

$ x(\boldsymbol{T} \cap \boldsymbol{P})=\boldsymbol{A}^{-1}(\boldsymbol{X}(\mathbb{T} \cap \mathbb{P})-b) $

(35)

4.3 积分

由(29)可知，新单元${\tilde W}$中第i种介质(如果只含有一种介质，则为纯网格)的面积为：

$ \tilde{V}_{\widetilde{W}, i}=\sum\limits_{k} \int_{\widetilde{W} \cap P_{i, k}} 1 \mathrm{~d} x \mathrm{~d} y $

(36)

因而新单元${\tilde W}$的面积为：

$ \tilde{V}_{\widetilde{W}}=\sum\limits_{i} \tilde{V}_{\widetilde{W}, i} $

(37)

新单元${\tilde W}$中第i种介质的体积分数为：

$ \tilde{V}_{\widetilde{W}}=\sum\limits_{i} \tilde{V}_{\widetilde{W}, i} $

(38)

同理, 新单元${\tilde W}$中第i种介质的中心坐标为：

$ \begin{array}{l} \tilde{x}_{c, \widetilde{W}, i}=\frac{\sum_{i} \int_{\widetilde{W} \cap P i, k} x \mathrm{~d} x \mathrm{~d} y}{\sum_{i} \int_{\widetilde{w} \cap P i, k} 1 \mathrm{~d} x \mathrm{~d} y}, \\ \tilde{y}_{c, \tilde{w}, i}=\frac{\sum_{i} \int_{\widetilde{w} \cap P_{i, k}} y \mathrm{~d} x \mathrm{~d} y}{\sum_{i} \int_{\tilde{w} \cap P_{i, k}} 1 \mathrm{~d} x \mathrm{~d} y} \end{array} $

(39)

对重构多项式(28)的密度进行积分，新单元${\tilde W}$中第i种介质的质量为：

$ \tilde{m}_{\widetilde{W}, i}=\sum\limits_{k} \int\limits_{\widetilde{W} \cap P_{i, k}} \rho_{i, k}(x, y) \mathrm{d} x \mathrm{~d} y $

(40)

因此，新单元${\tilde W}$中第i种介质的密度为：

$ {\tilde \rho _{\tilde W,i}} = \frac{{{{\tilde m}_{\tilde W,i}}}}{{{{\tilde V}_{\hat W,i}}}} $

(41)

同理, 分别对重构多项式(28)的动量和内能进行积分, 新单元${\tilde W}$中第i种介质的动量和内能分别为：

$ {{{(\tilde m\tilde u)}_{\tilde W,i}} = \sum\limits_k {\int\limits_{\tilde W \cap {P_{i,k}}} {{\rho _{i,k}}} } (x,y){\mathit{\boldsymbol{u}}_{i,k}}(x,y){\rm{d}}x{\rm{d}}y} $

(42)

$ {{{(\tilde m\tilde e)}_{\tilde W,i}} = \sum\limits_k {\int\limits_{\tilde W \cap {P_{i,k}}} {{\rho _{i,k}}} } (x,y){e_{i,k}}(x,y){\rm{d}}x{\rm{d}}y} $

(43)

因此, 新单元${\tilde W}$中第i种介质的速度和内能分别为：

$ \tilde{\boldsymbol{u}}_{\widetilde{W}, i}=\frac{(\tilde{m} \tilde{\boldsymbol{u}})_{\widetilde{W}, i}}{\tilde{m}_{\widetilde{W}, i}}, \quad \tilde{e}_{\tilde{W}, i}=\frac{(\tilde{m} \tilde{e})_{\widetilde{W}, i}}{\tilde{m}_{\tilde{W}, i}} $

(44)

新单元${\tilde W}$的平均速度为：

$ \tilde{\boldsymbol{u}}_{\widetilde{W}}=\frac{\sum_{i}(\tilde{m} \tilde{\boldsymbol{u}})_{\widetilde{W}, i}}{\sum_{i} \tilde{m}_{\tilde{W}, i}} $

(45)

4.4 后验校正

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

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

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

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

$ C^{P A D}=\left\{\begin{aligned} 1, & \tilde{\rho}_{\widetilde{W}, i}>0 \ \text { and } \ \tilde{e}_{ \widetilde{w}, i} >0 \\ & \ \ \ \ \ \ \text { and } \ \ \tilde{p}_{\widetilde{w}, i} >0 \\ 0, & \ \ \ \ \ \ \ \ \ \ \ \ \ \text { else } \end{aligned}\right. $

(46)

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

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

$ \begin{array}{c} D^{N A D}= \\ \left\{\begin{array}{l} 0, \widetilde{G}_{\widetilde{W}, i}<\min \limits_{k \in \widetilde{W} \cap P_{i, k} \neq \varnothing} \bar{g}_{i, k} \text { and } \max\limits _{k \in \widetilde{W} \cap P_{i, k} \neq \varnothing} \bar{g}_{i, k}<\widetilde{G}_{\widetilde{W}, i} \\ 1, \text { else } \end{array}\right. \end{array}. $

(47)

其中g_{i, k}的定义在(28)中给出，${{\tilde G}_{\tilde W, i}} = {{\tilde \rho }_{\tilde W, i}}, {{\tilde \rho }_{\tilde W, i}}{{\tilde u}_{\tilde W, i}}{\rm{ }}, {{\tilde \rho }_{_{\tilde W, i}}}{{\tilde e}_{_{\tilde W, i}}}$。

因此，如果C^PAD和D^NAD都等于1，那么单介质子多边形被认为是“好的”。如果一个新单元中所有单介质子多边形都是“好的”，那么它就被标记为“好”单元；否则，新单元是“坏”单元。当一个新单元${\tilde W}$被标记为“坏”单元时，它必须降阶来计算。降阶是指用拉氏网格的单元均值作为重构多项式来进行重映。请参见图 7中的简单说明。

5 数值算例 5.1 二维周期涡问题^[36]

有一个平均对角流ρ=1，p=1和(u, v)= (1, 1)。在这个平均流的速度(u, v)和温度T=ρ/p上加一个等熵涡扰动，对熵S=p/ρ^r不加扰动，即$，$(\delta u, \delta v) = \frac{\varepsilon }{{2{\rm{ \mathsf{ π} }}}}{{\rm{e}}^{0.5\left( {1 - {r^2}} \right)}}( - \bar y, \bar x), \delta T = - \frac{{(r - 1){\varepsilon ^2}}}{{8r{{\rm{ \mathsf{ π} }}^2}}}{e^{\left( {1 - {r^2}} \right)}}$, δS=0其中$(\bar x, \bar y) = (x - 5, y - 5), {r^2} = {{\bar x}^2} + {{\bar y}^2}$，涡的强度ε=5, 该问题存在解析解^[36]。计算区域取为[0, 10]×[0, 10]。在两个方向上都使用周期边界条件。时间计算到t=1。用本文的间断有限元拉氏方法和MMALE方法分别在四套网格上进行计算。表 1和表 2分别给出了密度，动量和总能的L₁和L₂误差以及相应的收敛阶。从表 1和表 2中可见，间断有限元拉氏方法和MMALE方法均能达到2阶精度。

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

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

$ p^{*}(t)=\left\{\begin{array}{ll} 13, & t \in[0,0.04] \\ 13-12.5 \frac{t-0.04}{0.125-0.04}, & t \in[0.04,0.12] \\ 0.5, & t \in[0.125,0.3] \end{array}\right. $

为了研究圆柱内爆中Rayleigh-Taylor不稳定性的时间演化，我们给出了轻、重气体界面的初始扰动，r_i^pert=r_i[1+acos(nθ)]。这里，r_i和r_i^pert分别表示无扰动和受扰动的半径。a表示扰动的振幅，整数n对应于扰动的模态。我们取模态n=8。两个不同的初始振幅a=2×10^-4和a=10^-3。初始时刻，采用射线网格把计算区域剖分成90×40个单元。如图 9所示，在径向上，轻气体均匀剖分成50份，重气体则剖分成不均匀的40份。将重气体中的单元按几何级数沿径向上进行分级，以便在界面处获得质量匹配网格。利用本文的间断有限元拉氏方法和MMALE方法分别计算这个问题。

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

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

当取a=10^-3时，图 12给出了t=0.25时用拉氏方法计算的网格和局部放大图，此时网格是没有翻转单元的。当t=0.3时，用拉氏方法计算的结果是不准确的，此时网格扭曲严重而且存在翻转单元，如图 13所示。然而，用MMALE方法运行到最终时刻没有任何问题。图 14给出的是t=0.3时用MMALE方法计算的网格与界面和密度等值线图。该算例的计算结果表明了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的一致网格。在靠近原点的第一个单元上给定内能e₀=400。在初始时刻，界面位于以原点为圆心，半径为0.43的圆上。注意到在混合单元中，这两种介质都是理想气体，具有相同的绝热指数γ=1.4。我们将它们视为混合单元，以便将数值解与精确解进行比较。图 15给出了t=0.5, 0, 75, 1时的网格和界面。图 16是t=1时刻的密度等值线图和以半径为变量的密度图。从图 16可知。计算得到的激波位置是准确的，数值解与精确解很接近。

5.4 二维Rayleigh-Taylor不稳定性问题^{[6, 36]}

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

$ \left\{\begin{aligned} p_{h}(x, y)=& 1+\rho_{h} g(y-1), \quad y>y_{i}(x) \\ p_{l}(x, y)=& 1+\rho_{l} g_{y}\left(y_{i}(x)-1\right)+\\ & \rho_{l} g_{y}\left(y-y_{i}(x)\right), \quad y \leqslant y_{i}(x) \end{aligned}\right. $

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

5.5 激波与氦气泡相互作用问题^{[6, 36]}

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

5.6 三重点问题^[6]

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

6 结论

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