0 引 言

通常ALE(Arbitrary Lagrangian-Eulerian)方法求解流体力学方程组时包含三个步骤：1) 拉氏计算步，即求解Lagrangian形式的流体力学方程组，将时间推进一步。在此过程中，随着流体的运动，网格会发生变形，甚至可能出现网格扭曲或翻转的情况；2) 网格重分步，即生成一个质量较好的新网格；3) 物理量重映步，即将旧网格上的物理量映射到新网格上。Li 等人^{[1, 2, 3]}研究了拉氏计算步，主要是针对可压缩欧拉方程的半拉格朗日微分形式，推导出它的积分弱形式，并用间断有限元方法对其进行空间离散，最后得到了中心型拉格朗日格式，避免了求解拉格朗日坐标系映射到欧拉坐标系下的雅可比矩阵，并且计算网格随着流体运动。当网格变形较小时，可以只进行拉氏计算步，来追踪物质界面；但当网格变形较大时，需要将网格重分，使重分后的网格具有更好的几何品质，并将拉氏计算步得到的物理量映射到新网格上。可见物理量重映步是ALE^{[4, 5]}方法的重要组成部分。

重映算法可以分为两类：积分重映和对流重映。积分重映对新旧网格的关联性要求不高，可以出现跨网格的情形。所谓跨网格是指下面的情形之一：一是新旧网格有相同的拓扑结构，单元数目相等，但新网格移动的距离较大，出现了跨网格情形；二是新旧网格有相同的拓扑结构，但单元数目不相等，比如说粗网格重分为细网格；三是新旧网格的拓扑结构不一样，比如说结构四边形网格重分为任意多边形网格。积分重映比较自然，但需要计算新旧网格的交点，计算量较大。对流重映的计算效率比较高，但它要求新旧网格上物理量的交换只在相邻单元中进行，因此要求新旧网格离得比较近，不能出现跨网格的情形。在实际计算中，我们还希望重映算法满足以下几个条件：一是守恒性，对于气动力学方程来讲，要求满足质量、动量和能量守恒；二是保界性，是指重映后物理量的值不能超过重映前物理量的最小和最大值；三是相容性，当新旧网格重合时，重映前后物理量的值不变；四是保证二阶或二阶以上收敛。

Dukowicz^[6]和 Miller^[7]分别给出了二维空间中的二阶重映方法；Grandy^[8]给出了三维空间中的一阶重映方法；Dukowicz^[9]和 Mosso^[10]分别在各自的文章中给出了三维空间中的二阶方法；Leslie^[11]和 Leonard^[12]则分别利用一维重映算法给出了多维空间中重映的分离算法。在这些重映方法中，具有代表性的一类局部ALE重映方法是通过适当的对流方程来构造重映过程。研究发现这类方法可以克服许多实际困难，所以在处理一些特殊应用问题上十分有效。此类方法的文献有很多^{[10, 13, 14, 15, 16, 17, 18, 19, 20, 21]}，其中Dukowicz^[20]表明当时间步充分小时，重映过程可以写成一个通量形式的对流算法，这使得重映方法变为一个简单、有效的标准对流格式，并且这种格式保留了一般重映方法的很多优势。

虽然上述用对流方程描述的重映算法有很多优势，但是对流方程和守恒重映之间的联系不易理解。Margolin和 Shashkov^[22]构造了与对流方程无关的重映算法，首先将新单元表示为新旧单元相交部分的组合，并以此为基础给出了一种相应的守恒积分重映算法。接下来Margolin 和 Shashkov又将新单元质量表示为相应旧单元质量加上其与周围单元的质量交换，再通过对密度函数的分片常数重构，得到了一阶格式的对流重映算法。为了得到二阶精度的重映算法，Margolin 和 Shashkov 分析了上述一阶对流重映算法的误差，通过对误差做补偿，得到了二阶、保正的重映算法。Kucharik等人^[23] 是首先利用已知的旧网格上单元平均值，重构出旧网格上密度的分片一次多项式逼近，再利用该分片一次多项式重映出新网格上密度的单元平均值，最后通过修补算法得到了一种保线性且保界的重映格式。

上述重映算法^{[22, 23]}可总结为：已知旧网格上单元平均值，重映后得到新网格上单元平均值。当ALE方法的拉氏步采用有限体积法离散时(例如文献^[24])，只需知道重分网格上物理量的单元平均值，就可以更新tⁿ⁺¹时刻物理量的值，因此可用文献^[22]和文献^[23]的重映算法。但当ALE方法的拉氏步采用间断有限元法离散时(例如文献^{[1, 2, 3]})，必须知道重分网格上物理量的分片多项式，才能更新出tⁿ⁺¹时刻物理量的值，此时上述重映算法不再适用。不过我们可以以上述重映算法为基础，先重映出物理量的单元平均值，再通过重构或者其他方法得到重分网格上的分片多项式。

因此当ALE方法的拉氏步采用间断有限元法离散时，以文献^{[1, 2, 3]}中的二阶拉氏格式为例，本文给出了相应重映步的对流守恒重映算法。首先在已知旧网格上分片一次多项式的基础上，采用文献^[22]和文献^[23]的重映和修补算法，得到重分网格上物理量的单元平均值，再通过极小化问题和Van Leer限制器，得到重分网格上物理量的梯度，由此得到重分网格上的分片一次多项式。

1 问题描述

考虑二维计算区域Ω，假设它是一个任意的多边形区域，给定Ω上一个多边形剖分{C_i,i=1,...,i_max}，其中单元C_i可以是非凸的。与单元C_i共点或共边的单元称为它的邻居，单元C_i所有邻居的集合称为它的邻域，记为C(C_i)。我们称{C_i,i=1,...,i_max}为旧网格，旧网格通过节点的小位移移动得到新网格{$\tilde{C}$_i,i=1,...,i_max}。在接下来的讨论中，假设新单元$\tilde{C}$_i包含在旧单元C_i和它的邻域中。下面以密度函数为例，叙述重映问题。

假设有一个定义在区域Ω上的正函数ρ(x,y)>0，称为密度，但是ρ(x,y)的具体表达式未知，仅知道ρ(x,y)在旧网格上的分片一次多项式逼近：

其中，(x_i,y_i)是单元C_i的中心，。V(C_i)>0是单元C_i的面积，(ρ_ix,ρ_iy)是密度梯度。显然单元质量为： 这表明近似函数ρ_i(x,y)保证ρ(x,y)单元平均值不变，ρ_i=$\frac{{{m}_{i}}}{V\left( {{C}_{i}} \right)}$是单元密度平均值。区域Ω上总质量为： 本文最终目的是在保证总质量M不变的情况下，找到新网格上密度的分片一次多项式逼近： 其中$\overline{{\tilde{\rho }}}$_i是新单元$\tilde{C}$_i上密度的近似平均值，($\overline{{\tilde{x}}}$_i,$\overline{{\tilde{y}}}$_i)是新单元$\tilde{C}$_i的中心，(${\tilde{\rho }}$_ix,${\tilde{\rho }}$_iy)是新单元上密度函数的梯度。

2 重映新网格上单元平均值 2.1 重映单元平均值

本节采用文献^[22]中对流重映算法，计算新单元上密度的平均值。下面我们简要回顾一下该重映算法，具体细节参见文献^[22]。

首先如文献^[22]给定一组结构四边形网格，因此需适当调整相应记号。网格顶点记为P_i,j=(x_ij,y_ij),i=1,...,m;j=1,...,n。单元标号记为(i+$\frac{1}{2}$,j+$\frac{1}{2}$),i=1,...,m－1;j=1,...,n－1，单元C_{i+$\frac{1}{2}$,j+$\frac{1}{2}$}的四个顶点是：{P_i,j,P_i+1,j,P_i+1,j+1,P_i,j+1}。那么旧单元C_{i+$\frac{1}{2}$,j+$\frac{1}{2}$}上密度函数的一次多项式逼近相应表示为：

假设定向边F_{i,j+$\frac{1}{2}$}={P_i,j+1,P_i,j}是单元C_{i+$\frac{1}{2}$,j+$\frac{1}{2}$}与C_{i-$\frac{1}{2}$,j+$\frac{1}{2}$}的公共边，位移向量场使顶点P_i,j移动到新位置$\tilde{P}$_i,j,P_i,j+1移动到新位置$\tilde{P}$_i,j+1,因此形成一个定向四边形δF_{i,j+$\frac{1}{2}$}={P_i,j,$\tilde{P}$_i,j,$\tilde{P}$_i,j+1,P_i,j+1}称为边移动扫过的区域。由文献^[22]可知边F_{i,j+$\frac{1}{2}$}移动扫描区域上带符号的面积|δF_{i,j+$\frac{1}{2}$}|可通过如下线积分来表示：

式(6)的符号由定向四边形δF_{i,j+$\frac{1}{2}$}={P_i,j,$\tilde{P}$_i,j,$\tilde{P}$_i,j+1,P_i,j+1}的顶点顺序决定。由定向的选择可知图 1中的情形为|δF_{i,j+$\frac{1}{2}$}|<0。

同样可以定义任意多项式在多边形上的符号积分。在本节下面描述中，线性函数在多边形上积分都是指在这个多边形上的符号积分。

定义新单元上近似质量$\tilde{m}$_{i+$\frac{1}{2}$,j+$\frac{1}{2}$}为相应旧单元上质量m_{i+$\frac{1}{2}$,j+$\frac{1}{2}$}加上与其邻近单元的质量交换：

这里f是边数值通量。与文献^[22]不同，本文已知的是旧网格上分片一次多项式逼近，所以以f_{i,j+$\frac{1}{2}$}为例，定义为： 其中， 这意味着在边移动扫过的区域中，面积符号决定了被积函数取为边哪一侧单元上的函数。

最后，由式(7)计算出近似质量$\tilde{m}$_{i+$\frac{1}{2}$,j+$\frac{1}{2}$}，那么新单元上密度平均值为：

2.2 修补算法

具体计算时，一般要求重映算法具有保界性，也就是说：如果一个新单元包含在相应旧单元的邻域中，则有:

其中I(C_i)是C(C_i)中所有单元标号的集合。文献^[23]给出了一种基于质量守恒的、局部单元质量重新分配的修补方法——局部调整不在范围内的值到给定界限内。将这种方法用在本文重映过程(重映出单元平均值之后)，下面给出该方法的简洁描述。

首先，对每个单元C_i选择一个决定它单元平均值范围的邻近单元集合R(C_i)，R(C_i)包含单元C_i本身，那么密度平均的局部最小、最大值分别为：

其中r_i是R(C_i)中所有单元标号的集合。

修补过程中首先要检查新单元平均值$\overline{{\tilde{\rho }}}$_i是否在最小、最大值范围内。如果在范围内就不进行下面的修补过程。如果不在局部极值范围内，则对它进行修补，使修补后的值在最小、最大值范围内。以$\overline{{\tilde{\rho }}}$_i<ρ_i^min为例：

如果$\overline{{\tilde{\rho }}}$_i<ρ_i^min，则需要向这个单元添加质量：

因为要求总质量守恒，所以这个需要添加的质量必须从其它单元获得。首先检查在不超过邻近单元平均值最下界的情况下，可以从周围取出的总质量值。也就是说对任意k∈r_i，可以安全从这个单元取走的质量为： 那么总的可用质量是 如果总的可用质量足够提供单元$\tilde{C}$_i的需求质量，即： 那么单元$\tilde{C}$_i的质量和密度平均值可以修补到最下界， 邻近单元质量适当的减少为： 并且单元$\tilde{C}$_i和它邻近单元的总质量不变。对于密度平均值$\overline{{\tilde{\rho }}}$′_k： 这就完成了单元$\tilde{C}$_i的修补工作。

更多修补细节参见^[23]，并且修补算法中相邻单元间的质量转移，与共有边上的质量通量无关，只是将临近单元中可以安全取走的质量直接加入当前单元。

3 重构分片一次多项式 3.1 极小化问题求解梯度

通过第2节的算法，可以得到新网格上密度的单元平均值$\overline{{\tilde{\rho }}}$′_i，仍记为$\overline{{\tilde{\rho }}}$_i。但是根据文献^{[1, 2, 3]}中描述的二阶拉氏方法可知，为了下一步拉氏运算，必须由新网格上密度的单元平均值，重构出新单元上近似线性函数${\tilde{\rho }}$_i(x,y)。下面通过极小化问题得到该函数的梯度(${\tilde{\rho }}$_ix,${\tilde{\rho }}$_iy)。

以新单元$\tilde{C}$_i为例，已知C($\tilde{C}$_i)是单元$\tilde{C}$_i的邻域，令

将H看成是${\tilde{\rho }}$_ix、${\tilde{\rho }}$_iy的函数，通过求函数H的极小值，可求出${\tilde{\rho }}$_ix、${\tilde{\rho }}$_iy。这仅须满足下面两个条件： 对C($\tilde{C}$_i)中任何一个单元，比如$\tilde{C}$_k，记$\overline{{\tilde{x}}}$_ki=$\overline{{\tilde{x}}}$_k－$\overline{{\tilde{x}}}$_i，$\overline{{\tilde{y}}}$_ki=$\overline{{\tilde{y}}}$_k－$\overline{{\tilde{y}}}$_i。

那么将式(20)代入式(21)可得一元二次方程组：

解上述一元二次方程组便可得到新单元上梯度(${\tilde{\rho }}$_ix,${\tilde{\rho }}$_iy)。至此，就得到了新单元$\tilde{C}$_i上的一次多项式逼近式(4)。

3.2 Van Leer 限制器

在整个重映算法中，我们还希望新单元上分片一次多项式不出现新的极值。因此使用文献^[20]中 Van Leer 限制器对式(4)中梯度(${\tilde{\rho }}$_ix,${\tilde{\rho }}$_iy)进行限制。下面以单元$\tilde{C}$_i为例，给出限制器的描述。

设$\overline{{\tilde{\rho }}}$_i^min、$\overline{{\tilde{\rho }}}$_i^max是密度平均值在邻域单元中的最小值和最大值，${\tilde{\rho }}$_i^min、${\tilde{\rho }}$_i^max是由式(4)计算${\tilde{\rho }}$_i(x,y)在单元$\tilde{C}$_i上的最小、最大值。由于${\tilde{\rho }}$_i(x,y)是一次多项式，所以它的最小、最大值应该在单元顶点上取得。令

那么最终得到新单元$\tilde{C}$_i上一次多项式为：

至此整个重映过程完成，它把旧网格上分片线性数值解重映为新网格上分片线性数值解，并且满足保界性和守恒性。

4 数值算例

将上述加修补和限制器的重映方法记为“SRRL”方法，将不加修补和限制器的重映方法记为“SRUN”方法。为了检验SRRL方法的有效性，对给定函数在一系列网格上用SRRL方法做重映，并与SRUN方法的结果作比较。在下面误差计算中 L₁模和L_∞模定义为：

其中ρ(x,y)是真解，${\tilde{\rho }}$(x,y)是重映方法得到的分片一次多项式逼近。

4.1 交替系列网格

初始网格是[0,1]×[0,1]上的一致均匀网格。下一套网格是通过对该一致均匀网格作一个小的随机位移得到的随机网格。交替系列网格就是由一致网格和相应随机网格交替产生的系列网格，即一致网格——随机网格——再恢复到一致网格，这样交替进行的一系列网格。在数值算例中作偶数次网格交替和重映，使最终网格恢复为初始网格，计算出整个重映过程的累积误差。当细分网格时，增加重映的步数，从而使总的位移保持近似相同。

算例 1。分别用SRRL、SRUN方法在交替网格上重映函数(x－0.5)⁴+(y－0.5)⁴。选取三套不同水平的网格：i_max=j_max=50,100,200和相应时间步数：n_max=10,20,40来考查SRRL方法的收敛性。

表 1和表 2分别给出了SRUN和SRRL方法的数值结果。从表中可以看出，在误差允许范围内，这两种方法对于光滑函数的误差基本相等。文献^[23] L₁模误差在剖分为100、200时的收敛阶分别是1.862、1.927； L_∞模收敛阶分别是1.784、1.741。显然本文SRRL方法收敛更快一些，更趋于2。

算例 2。在50×50剖分的交替系列网格上用SRRL和SRUN方法重映一个三层塔式函数，函数取值范围是[0,1]。令d(x,y)=|x－0.5|+|y－0.5|，那么这个函数表达式为：

表 3给出了在最终网格上用SRRL和SRUN方法重映得到的函数最大值和最小值。从表 3中可以看出，该结果与文献^[23]中数值结论一致：加过修补和Van Leer限制器的SRRL方法是保界的，并且SRRL方法L₁模误差比SRUN方法的误差小。

4.2 随机网格

本节使用一系列的随机网格，其中每一套网格都是通过对一致均匀网格做一个随机扰动得到：

这里－0.5≤r_iⁿ,r_jⁿ≤0.5是随机数，γ=0.7。在这种随机情况下，网格不能回到初始网格状态，所以严格上来说，不能用来分析方法的收敛性。

算例 3。在此用SRRL方法和SRUN方法在随机系列网格上重映激波函数:

图 2给出了SRUN和SRRL方法计算得到的激波函数曲面图和等值线图，可以看出SRUN方法的数值解在间断处出现了数值振荡，并且数值解超出了给定范围，出现负值。而SRRL方法的数值解在间断处没有振荡，并且数值解在给定界限内。表 4给出了在最终网格上两种重映方法得到的函数最大和最小值，结果显示SRRL方法是保界的，并且SRRL方法L₁模误差比SRUN方法的误差小。文献^[22]中mc方法在γ=0.7时出现负值，da和mb方法的数值解仍为正值。通过表 4对比SRRL、da以及mb方法的结果可以发现，本文SRRL方法比da和mb方法的结果更精确。

4.3 初始随机网格的连续光滑化

在本次测试中，选取与ALE数值模拟方法相接近的一系列网格，其中每一组网格都是由前一组网格光滑化得到：

y方向的节点选取类似。初始网格是通过对一致均匀网格作随机扰动得到的：

这里－0.25≤r_i,r_j≤0.25是随机数，h=1/(i_max－1)。边界上节点只有一个坐标被扰动。

算例 4。用SRRL方法和SRUN方法在连续光滑化网格上重映激波函数式(31)。图 3给出了SRUN和SRRL方法计算的激波函数侧面图，可以看出SRUN方法的数值解在间断处出现了小振荡，并且数值解超出给定范围，而SRRL方法的数值解在间断处没有振荡。表 5给出了两种方法的数值结果，结果显示SRRL方法是保界的，并且SRRL方法L₁模误差比SRUN方法误差小。比较算例3和算例4的结果，可以发现同样的函数在连续光滑化网格上重映结果更好，这表明本文重映算法适合用于ALE方法。

5 结 论

本文针对ALE方法中由间断有限元法求解流体方程的拉氏步，给出了一种新的重映算法。算法由两步组成：第一步在已有对流重映算法基础上，计算出新网格上物理量的单元平均值，并通过修补算法，使之满足保界性；第二步是通过重构得到新单元上函数的近似梯度，并通过Van Leer限制器对梯度进行限制，使之不出现新极值。这样就得到了新网格上的分片一次多项式。最后通过一系列数值算例表明，加上修补和限制器的重映算法和不加修补和限制器的重映算法，对于光滑函数L₁和L_∞模误差都是二阶收敛，且误差相差不大。但不加修补和限制器时，会出现数值结果越界的情形，加上修补和限制器之后，重映算法是保界的，并且误差更小。本文重映算法可以直接用于Euler方程或NS方程，在重映物理量ρV和ρE时(其中V是速度，E是单位质量的总能量)，只是将文中的密度ρ替换为ρV和ρE，并且可以保证动量和能量守恒，以及整个方程组的相容性。