王晨涛，马天宝，李坤

（北京理工大学 爆炸科学与技术国家重点实验室，北京 100081）

基于连续介质力学求解双曲守恒律方程得到的解通常具有奇异性，即解存在强间断，例如激波、接触间断，或者稀疏波之类的弱间断. 针对求解间断解发展了很多高分辨率和高精度的数值算法. 高分辨率的算法可以追溯到早期的TVD 总变差不增算法，但是算法精度不高于二阶. 随后，发展了一系列的高精度数值格式，例如ENO，RKDG，WENO[1-2]，SD，SV等. 但是高精度的数值格式数值色散较强，使得解的震荡较强. 通常需要额外的限制震荡的方法来抑制震荡. 因此，为了减少数值震荡，提高间断分辨率，从削弱方程奇异性的角度发展了弧长算法. RISK[3]通过引入弧长参数，增加了一个弧长约束方程，使得求解过程的奇异性得到减弱或消除. CHAN[4]通过引入伪弧长控制参数将其引申为伪弧长算法，并应用于求解非线性方程的的奇异性问题. 武际可[5]等人应用伪弧长算法解决了包含奇异性的Burgers 方程.TANG[6]提出了包含伪弧长参数的移动网格控制方程，并给出了在二维情况下网格最优分布的自适应函数表达式. 宁建国等发展了二阶伪弧长算法，建立了伪弧长算法稳定性[7]理论，得到了伪弧长算法的收敛性结论，并成功将伪弧长算法运用到爆炸冲击波强间断问题的求解. 由于在变形的物理空间不易构造高阶数值格式，因此上述的弧长算法限制在了二阶精度.

在多介质流体的求解过程中，围绕物质界面的追踪发展了经典的VOF[8-9]方法，Level Set[10]法，处理多介质的五方程模型[11]方法，多介质ALE[12]算法，以及基于界面形心重构物质界面的MOF[13]方法.Level Set 通过求解一个对流方程实现界面的隐式追踪，无需对界面进行繁琐的重构，因此具有计算简单，界面分辨率高等优点. 由于界面两侧流体状态方程的差异，容易在界面附近产生需要虚假振荡，因此需要合理定义界面的边界条件. RGFM[14]通过求解一个黎曼问题，修改界面附近的物理量，有效降低了界面处的非物理反射引起的数值振荡.

本文将伪弧长算法和WENO 格式以及Level-Set和RGFM 相结合，发展了多介质流的高精度伪弧长算法，为了构造高精度格式，采用坐标变换将控制方程映射到正交均匀的弧长坐标系中，计算过程在均匀正交的弧长坐标下进行计算. 为了消除变形网格几何误差，采用满足几何守恒律[15-17]的表达式计算. 针对网格移动以后距离函数的插值提出了基于泰勒插值的三阶非守恒格式. 计算结果表明，本文的算法在保持较高的精度的同时也削弱了间断处的数值振荡，大幅提升了间断的分辨率.

1 控制方程

1.1 物理空间下的控制方程

可压缩欧拉方程在物理空间的张量形式为

1.2 弧长空间下的控制方程

可以看出，所有的表达式均是守恒格式，因此方便用中心差分格式计算，从而可以消除因网格变形造成的几何误差.

2 伪弧长算法

2.1 伪弧长算法几何意义

为了削弱方程的奇异性，通过引入一个弧长约束方程，即

其中s为弧长变量. 可以发现，当弧长变量为常值时，dw越大，则 dx越小，它的物理意义就是物理空间下变量的梯度越大，网格的尺寸越小，也就是网格朝着间断点移动. 相关的几何意义参见图1，从图中可以看出，经过引入弧长约束方程，使得网格自适应朝着间断移动，原始物理空间下陡峭的物理量在弧长空间中得到了平滑的过渡，因此削弱了方程的奇异性.

图1 1D 伪弧长算法的几何意义Fig. 1 Geometric significance of one-dimensional PALM

二维空间和一维空间类似，基于维度解耦的思想，引入两个弧长约束方程，把它写成向量的形式，即

相关的几何示意图如图2 所示. 它和一维情况下的物理意义类似，都是引入约束方程使得网格朝着间断移动，在物理空间的表现形式就是网格汇聚到了间断附近，然后再弧长空间下网格依旧是正交均匀的.

图2 2D 伪弧长算法的几何意义Fig. 2 Geometric significance of two-dimensional PALM

2.2 伪弧长算法网格移动

伪弧长算法的网格满足一个能量极值函数：

然后可以给出满足弧长控制函数的网格迭代方程(8)，即

至此，网格便可以采用数值方法迭代计算，例如采用收敛较快的高斯-赛德尔迭代，格式如下.

3 数值格式

3.1 控制方程的离散

3.2 守恒变量的更新

网格移动完以后，需要更新在新网格上的物理量，本文采用守恒插值[6]来计算：

4 Level-Set 技术

4.1 控制方程的离散

4.2 距离函数的非守恒插值

网格移动完以后，需要更新新网格点上的Level-Set 值，本文基于泰勒插值，给出三阶插值的表达式，由于新旧网格间距较小，因此将新网格点在旧网格点作泰勒展开，有：

由于计算过程是在弧长坐标系下，因此将上式映射至弧长空间，可得：

4.3 多介质流高精度伪弧长算法流程

⑤根据式(10)和(18)进行控制方程以及符号距离函数在时间步上的更新. 更新完以后根据距离函数的符号来重新分配两种介质.

⑥判定计算时间是否到达. 若到达，则终止计算，否则转到步骤①.

5 数值算例

为便于计算，以下算例计算均采用无量纲计算.

5.1 精度测试

用一个包含精确解的一维的Euler 方程来验证高阶伪弧长算法PALM-5 的收敛精度. 初值条件和精确解为

从表1 的计算结果可以看出，弧长坐标系下的高阶伪弧长算法保证了较高的精度，没有因为网格的变形而损失精度. 这说明本文的算法满足了高精度的性质.

表1 固定网格下WENO 和PALM-5 的收敛阶Tab. 1 The accuracy of WENO, PALM-5 schemes

5.2 一维爆炸波

除图3(g)和3(h)算例以外其余均采用的网格数目是300，参照解是在WENO 格式下用20 000 个网格计算得到的. 从图3(a)和(3b)中可以看出，采用相同的格式伪弧长算法比固定网格的算法分辨率要高很多，对间断的捕捉能力很强. 而且5 阶伪弧长算法比2 阶伪弧长算法分辨率更高，对间断的捕捉能力要高很多，耗散很低. 对比图3(c)和3(d)可以看出采用不同的伪弧长参数得到的结果差异较大，本算例中第二套系数的分辨率要优于第一套系数，可以从图3(e)和3(f)中看出，第二套系数网格移动的更加剧烈，因此对间断捕捉的分辨率更高. 因此为了提高分辨率，可以适当的调整伪弧长控制函数. 图3(g)和3(h)均采用了220 和440 个网格数进行对比，可以看出，PALM-5 在220 个网格数目下和WENO 在440 个网格数目下的近似分辨率几乎完全相同，而从表2中可以看出采用PALM-5 算法在近似分辨率下需要的时间大大缩短.

表2 近似分辨率下WENO 和PALM-5 的时间对比Tab. 2 Time comparison of WENO and PALM-5 at approximate resolution

图3 爆炸波问题的结果Fig. 3 Results of blast wave problem

5.3 一维激波击打气-气问题

这是一个激波击打气-气界面问题的多介质算例，初值条件为

5.4 二维水下爆炸

从计算结果可以看出，两种计算结果的对称性都比较好，但是固定网格下捕捉到的界面扩散严重，而伪弧长算法下界面锐利清楚，且解的光滑性非常好，由于状态方程以及物理量的巨大差异，因此固定网格下的结果在界面附近有轻微的数值震荡，而采用伪弧长算法计算的结果非常光滑，完全消除了数值震荡,如图5 所示. 从图5(c)中可以看出网格的自适应性非常好，而且网格足够光滑，没有出现间断的网格点.

图5 水下爆炸结果Fig. 5 Results of two-dimensional underwater explosion problem

6 结 论

将伪弧长算法和高阶WENO 格式相结合发展了高精度伪弧长算法，结合RGFM 和Level-Set 将其应用到多介质的计算中. 为了构造高精度的数值格式，所有的计算均在坐标变换之后的弧长坐标系下计算，此外，在多介质的计算中，为了使的距离函数保持物理性质，采用了法向速度的表达式，针对网格移动以后距离函数的插值，给出了3 阶的非守恒插值形式.根据数值算例，得出以下结论：

①经过坐标变换以后，高精度伪弧长算法仍旧可以保持较高的精度；

②爆炸波算例说明高阶伪弧长算法比低阶伪弧长算法捕捉界面的分辨率更高，而且相对固定网格下结果提升显著，在近似分辨率下伪弧长算法需要的时间大大缩短，计算效率较高；

③多介质算例说明，高精度的伪弧长算法结合RGFM 可以成功的计算多介质问题，并且有效降低了界面以及其他间断处的数值震荡，使得解的光滑性非常好，网格自适应朝着间断移动，使得解捕捉间断的能力非常强. 解决了低阶格式耗散强，高阶格式震荡强的问题. 发展了高分辨率，低震荡的算法.