一种求解一维理想磁流体方程组的保正拉氏方法

邹世俊, 蔚喜军, 戴自换

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

0 引言

可压缩磁流体方程组(MHD)在天体物理、可控热核聚变、金属冶炼等众多科研领域中扮演着重要的角色。此方程组是将气体动力学方程组与Maxwell 方程组结合起来，广泛的应用于描述导电流体在磁场中的动力学现象。当忽略黏性，电阻及相对论效应时，上述方程转化为理想可压缩磁流体方程组。这组方程组十分复杂以至于很难对其进行数学分析和精确求解，因此，对于理想磁流体方程组设计精确鲁棒的数值方法有着十分重要的意义。理想磁流体方程组可以被写为如下形式

在物理上，上述系统所描述的物理现象中，密度ρ和热力学压力p总是非负的。但是，一个逼近上述系统的数值格式所得到的近似解并不能总是保持这种正性，这会引起系统失去双曲性并导致数值不稳定。为了解决这一问题，许多研究者提出了不同的数值技巧。然而，绝大多数此类方法都是在欧拉框架下提出的，即网格在空间中是固定的。在流体力学的数值模拟中，有另一种框架也被经常使用，这就是拉格朗日框架。在拉氏框架下，计算网格随流体运动。因此，拉氏方法更加适合处理多介质及自由界面问题。但是，在计算流体力学中保正拉氏格式很少被研究，更不用说是针对于MHD 系统设计的方法了。Bezard 和Desp´res[3]针对一维拉氏MHD 方程组设计了一种数值格式，这一格式满足熵不等式并在适当的CFL 条件下满足保正性质。Gallice[6]运用简单黎曼解的概念提出了一种熵保正Godunov 型格式，这一格式对拉氏和欧拉MHD 系统都可以保接触间断。近些年来，有一些保正拉氏格式[4,9]被提出。这些格式可以很好的运作于欧拉方程并能得到高阶精度。但是，所有这些格式都不能推广到MHD 系统。

本文中，我们将关注的焦点放在针对于一维理想磁流体系统的保正拉氏算法设计上。基于文献[10]中在拉氏框架下针对二维可压缩理想MHD 方程组设计的RKDG 方法，我们推导出针对一维理想MHD 方程组的拉氏格式。之后，针对拉氏格式推广了HLLD 近似黎曼解[7]，这一黎曼解可以保持密度和热力学压力的正性。这种拉氏HLLD 数值通量的表达式与其在欧拉框架下的表达式并不相同。运用这一特别的HLLD 数值通量，我们可以构造出一种保正拉氏格式。最后，我们将会给出几个一维中具有挑战性的数值算例来证明我们方法的有效性。

本文结构如下，在第1 节中，我们给出针对于一维理想MHD 方程组的拉氏格式。之后我们给出拉氏HLLD 数值通量的推导过程及其保正性质。最后，通过运用拉氏HLLD 数值通量，我们可以构造出求解一维理想MHD 方程组的保正拉氏格式。在第2 节中，通过给出数值算例可以证明我们方法的保正性质。最终，我们会在第3 节中给总结。

1 求解一维理想MHD 方程组的保正格式

1.1 一维理想MHD 方程组的拉氏格式

众所周知，黎曼问题(8)真解的求解可能既困难又需要大量的计算，因此在构造Godunov 型数值格式时经常会使用近似黎曼解。基于HLLD 黎曼解[7]的良好性质，我们在针对一维MHD 系统的方法中使用此黎曼解。

1.2 拉氏HLLD 通量及其保正性质

图1 拉氏框架下HLLD 通量的黎曼扇

引理1 拉氏HLLD 求解器具有保正性质，如果SL和SR满足

1.3 拉氏保正格式

我们将(7)中的Uh取为分片常数并运用向前欧拉时间离散，可以得到拉氏框架下积分弱形式(6)的离散格式

图2 应用HLLD 求解器的拉氏格式在Ii×[tn,tn+1]中的波系演化

2 数值算例

在这一节中，我们给出几个一维中具有挑战性的数值算例来测试之前提出方法所具有的保正性质。若不使用保正方法，则所有这些算例在计算时都会产生负压。

2.1 超快稀疏波问题

在这一小节中，考虑一个存在超快稀疏波的激波管问题。其初始条件如下

在图3 中，我们分别给出u0= 3 和u0= 3.1 时在t= 0.05 时刻密度ρ，速度分量ux和总压的图像。如文献[7]中所述，在对众多近似黎曼解进行数值测试后可以发现，虽然在Mf= 3 时，所有测试过的近似黎曼解都不会产生负的密度和压力，但是在Mf= 3.1 时有些却无法计算下去。这是因为这些无法进行计算的近似黎曼解并不保正。在图3 中可以看到我们的格式对于Mf=3 和Mf=3.1 都运算的很好，这一数值结果可以说明拉氏HLLD 黎曼解的保正性质。注意到，图3 中的图像与文献[7]中所给出的图像有所不同，这是因为本文中的结果是运用拉氏方法得到的，网格随着流体移动到了计算区域两端，在x=0 附近网格非常稀疏，x=0 附近密度趋于真空。

图3 一维超快膨胀激波管问题的结果

2.2 一维真空激波管问题

在这一小节中，我们给出文献[5]中提出的一维真空激波管问题的数值结果。问题的初始条件如下

计算区域取为[−0.5,0.5]。γ= 5/3 并应用外流边界条件。这一问题用来说明我们的拉氏保正格式可以处理拥有极低密度和压力的问题。

在图4 中，我们给出在t=0.1 时刻2 000 网格上计算得到的密度和压力的数值结果。通过此图可以观察到，在x=−0.3 附近网格非常稀疏。与之前的数值算例类似，在拉氏方法中网格随着流体移动，因此x=−0.3 附近趋于真空。通过和文献[5]中的结果进行比较，我们可以确信无论是低密度还是低压力都捕捉得很好。注意到，如果不运用保正方法来计算此问题，则计算程序将会因为非物理解的出现在几个时间步内崩溃。

图4 运用保正拉氏格式在2 000 网格上得到的一维真空激波管问题的密度和压力

2.3 旋转阿尔文波脉冲问题

本节，我们考虑文献[1]中提出的旋转阿尔文波脉冲问题。给出初始条件如下计算使用周期边界条件，γ= 5/3。算例中，初始压力比总能的万分之一还要小。除此之外，该问题的解中还存在一个很强的阿尔文波间断，大多数数值格式在计算此问题时都会产生负压。

该问题的数值模拟采用了800 个计算网格。在图5 中，我们给出拉氏保正格式在t=0.156 时刻得到的总能和压力。注意到在解中存在两个脉冲，因此成功的数值模拟应该能很好的保持这一特点。与前面的数值算例类似，如果不使用保正方法计算此问题，则计算程序会在几个时间步内崩溃。

图5 运用保正拉氏格式在800 网格上，t=0.156 时，计算得到的旋转阿尔文波脉冲问题的总能和热压

3 总结

本文中，我们给出了一维理想MHD 系统的一种拉氏形式，并提出了一种在适当的信号速度下可以保持密度和热力学压力正性的拉氏HLLD 数值通量。运用这一HLLD 通量，我们构造出了一种针对一维理想MHD 系统的保正拉氏格式。最后，我们给出了几个具有挑战性的数值算例来证明本文提出的方法所具有的保正性质。

在将来的工作中，我们将会探索如何对于多维MHD 系统构造保正拉氏格式。在对拉氏MHD 系统特别是多维情形设计保证格式时存在数个困难。这些困难使得多维MHD 系统的保正拉氏格式构造变为了一个十分复杂的工作。因此，这一工作仍然需要继续研究。