非理想情况下磁化球形重质气体物理爆炸的数值模拟

林震亚; 陈志华; 刘迎; 洪延姬

doi:10.11883/1001-1455(2017)03-0422-09

非理想情况下磁化球形重质气体物理爆炸的数值模拟

doi: 10.11883/1001-1455(2017)03-0422-09

1.
南京理工大学瞬态物理国家重点实验室, 江苏南京 210094
2.
解放军装备学院激光推进及应用国家重点实验室, 北京 101416

基金项目:

国家自然科学基金项目 11272156

详细信息

作者简介:
林震亚(1990—)，男，博士

通讯作者:
陈志华, chenzh@mail.njust.edu.cn

中图分类号: O381
计量
- 文章访问数: 4357
- HTML全文浏览量: 1301
- PDF下载量: 474
- 被引次数: 0
出版历程
- 收稿日期: 2015-10-26
- 修回日期: 2016-01-05
- 刊出日期: 2017-05-25

Influence of nonideal magnetic field on physical explosion of spherical heavy gas

1.
State Key Laboratory of Transient Physics, Nanjing University of Science and Technology, Nanjing 210094, Jiangsu, China
2.
State Key Laboratory of Laser Propulsion and Application, Academy of Equipment, Beijing 101416, China

摘要

摘要: 采用磁流体力学(magnetohydrodynamics, MHD)方程组对非理想情况下磁场对球形重质气体物理爆炸的作用过程进行数值模拟，同时，为了保证每一步中磁场散度为零，引入由12-solve CTU算法衍生出的CTU+CT算法。结果清晰地显示了重质气体团在磁场影响下的整个爆炸过程。在非理想情况下，爆炸过程中气体团界面产生液滴状结构，随着气体团被压缩，不稳定性现象最终得到抑制。从结果中看出，电阻和双极扩散效应会阻碍磁场对气体团的作用，同时，双极扩散效应还会增大磁压力的作用范围。
- 不稳定性 /
- 爆炸 /
- MHD /
- 双极扩散 /
- 激波
Abstract: Magnetohydrodynamics equations are adopted to simulate the process of spherical heavy gas explosion. Meanwhile, in order to ensure the divergence of magnetic field is zero in each step, we use the CTU+CT algorithm which is derived by 12-solve CTU algorithm. The results clearly show the process of spherical heavy gas physical explosion with the influence of magnetic field. In the non-ideal case, the droplet-like structure appears on the interface of the gas cluster. With the gas cluster being compressed, the instabilities are being restrained in the end. As can be seen from the results, resistance and ambipolar diffusion effect will hinder the influence of magnetic field on gas cluster, at the same time, the ambipolar diffusion effect will increase the scope of the magnetic pressure.
- instability /
- explosion /
- magnetohydrodynamics /
- ambipolar diffusion /
- shock wave

HTML全文

当2类不同的流体被激波加速，界面上的小扰动会线性增长，接着会形成以“尖钉”和“气泡”为特征的非线性结构。这类现象被称为Richtmyer-Meshkov(RM)不稳定性。这种不稳定性先由R.D.Richtmyer^[1]在1960年理论预测，10年后，E.E.Meshkov^[2]用实验证明了R.D.Richtmyer的预言。RM不稳定性在惯性约束聚变(inertial confinement fusion, ICF)中扮演着重要的角色，由于界面不稳定性产生的流体混合是ICF中限制能量获取的重要因素。RM不稳定性同样存在于其他很多自然以及人造现象中，如超新星爆发、爆燃及爆炸中等。此外，RM不稳定性可以诱导流场的湍流混合^[3]。研究不稳定性对理解湍流机理很有帮助。

RM不稳定性和Rayleigh-Taylor(RT)不稳定性关系密切，它们有相似的特征，如“气泡”和“尖钉”的形成与增长。其中，“气泡”是轻质流体渗入重质流体的部分，“尖钉”是重质流体渗入轻质流体的部分。此外，由于不同的外部作用力，这2种不稳定性表现出不同的特征。RT不稳定性的驱动机制是均匀重力场对流体的作用，它使界面的小扰动指数增长。相对的，RM不稳定性由激波驱动，无论激波是由重质流体流入轻质流体时的界面碰撞还是轻质流体流入重质流体时的界面碰撞产生^[4-6]，它都不稳定。初始界面沉积的涡度控制不稳定性的增长率，接下来，该涡度由追赶上界面的反射波进行修正。

不稳定性的演变与系统的压缩能力有密切联系，实际上，固有的压缩性物理是其中重要的因素^[7]。压缩性现象，如激波的相互作用、激波的折射以及界面上的不稳定性，包括线性及非线性的增长和随后过渡到湍流的过程，对于理论学家和实验学家一直是巨大的挑战^[7-8]。C.A.Zoldi^[9]利用激波管开展了激波加载SF₆气柱RM界面不稳定性实验，并采用RAGE程序模拟了SF₆气柱在空气冲击波作用下界面的演化、发展过程以及后期空气和SF₆气体的混合。G.Layes等^[10]采用高速摄影方法研究了不同密度的气泡在激波作用后的变形过程。邹立勇等^[11]利用高速摄影测试技术，对弱激波冲击下空气中的SF₆气柱和气帘界面的演变过程进行了实验研究，对演变过程中涡的产生和发展进行了初步解释。

然而，在大多数实际应用中，如ICF或超新星爆发，RM不稳定性呈现弯曲的几何结构，显著增加了系统的复杂度。首先，非微扰系统没有弯曲几何区域中的解析解，然而该系统在平面几何区域中存在解析解; 再者，由于弯曲的几何结构，初始扰动即使在RM不稳定性的初始阶段也呈现非线性增长; 更进一步的，在反射激波的作用下，物质界面产生二次加速，从而使整个系统更复杂。目前，没有在圆形或球形几何下RM不稳定性的解析理论，数值模拟结果也很少。Zhang Qiang等^[12]开展了关于SF₆气体内爆及爆炸产生的圆形激波的数值模拟研究。同时，他们^[13]采用界面追踪的方法研究了圆形几何下由强激波驱动的不稳定界面的比例定律并发现了临界马赫数。S.Dutta等^[14]、J.Glimm等^[15]采用界面追踪法对球面几何下坐标对称流动的RM不稳定性进行了数值模拟，证明了关于激波马赫数对于流体混合统计的标度不变性，如增长率和体积分数。

RM不稳定性在很多领域^[16]的应用中，如惯性约束融合^[17]、天体物理现象^[18]等，流体可能是离子化的，因此会受磁场的影响。R.Samtaney^[19]通过数值模拟证明了RM不稳定性的增长会受到磁场的抑制, 他研究的问题关于二维理想磁流体力学(magnetohydrodynamics, MHD)中不同密度下由斜平面接触间断(contact discontinuity, CD)分离的导电流体中的激波相互作用; 结果表明了加入磁场时不稳定性的抑制作用主要由接触间断处激波折射过程的变化引起^[20]，这些变化阻止了接触间断处沉积物的流通。本文中，同样采用磁场来控制重质气体爆炸过程中的RM不稳定性; 同时，为了保证每一步中磁场散度为零，引入CTU+CT算法对其进行数值求解; 由于现实情况下不可能出现电导率为无穷的情形，拟随后加入电阻率及双极扩散效应并和之前的结果进行对比，从而分析每个因子对RM不稳定性的影响。

1. 方程与数值算法

除满足流体力学的条件，带电荷的流体团有形成团块更有利的条件，主要原因如下：(1)带电粒子之间的碰撞频率远大于中性粒子之间的碰撞频率，因此碰撞平均自由程较短。(2)当有磁场时，磁场对带电粒子有侧向约束。洛伦兹力使带电粒子在垂直于磁场方向只能保持在回转半径范围内运动。等离子体通常比中性粒子有更低的密度与更高的温度，由于爆炸可在短时间内完成, 爆炸的速度非常快，本文中只考虑高磁雷诺数的情况，采用MHD方程^[21]来研究爆炸中出现的不稳定性:

$\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \mathit{\boldsymbol{\nu }}} \right) = 0$

(1)

$\frac{{\partial \left( {\rho \mathit{\boldsymbol{\nu }}} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho \mathit{\boldsymbol{\nu \nu }} - \frac{{\mathit{\boldsymbol{BB}}}}{\mu } + {p^*}} \right) = - \nabla \cdot \mathit{\boldsymbol{ \boldsymbol{\varPi} }}$

(2)

$\frac{{\partial E}}{{\partial t}} + \nabla \cdot \left[ {\left( {E + {p^*}} \right)\mathit{\boldsymbol{\nu }} - \frac{{\mathit{\boldsymbol{B}}\left( {\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{\nu }}} \right)}}{\mu }} \right] = - \nabla \cdot \mathit{\boldsymbol{Q}} - \nabla \cdot \left( {\mathit{\boldsymbol{ \boldsymbol{\varPi} }} \cdot \mathit{\boldsymbol{\nu }}} \right)$

(3)

$\frac{{\partial \mathit{\boldsymbol{B}}}}{{\partial t}} - \nabla \times \left[ {\mathit{\boldsymbol{\nu }} \times \mathit{\boldsymbol{B}} - \eta \mathit{\boldsymbol{J}} - {\eta _{\rm{H}}}\left( {\mathit{\boldsymbol{J}} \times \mathit{\boldsymbol{B}}} \right) - {\eta _{{\rm{AD}}}}\left( {\mathit{\boldsymbol{J}} \times \mathit{\boldsymbol{B}}} \right) \times \mathit{\boldsymbol{B}}} \right] = 0$

(4)

${p^*} = P + \frac{{\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{B}}}}{{2\mu }}$

(5)

$\mathit{E = }\frac{p}{{\gamma - 1}} + \frac{{\rho \left( {\mathit{\boldsymbol{\nu }} \cdot \mathit{\boldsymbol{\nu }}} \right)}}{2} + \frac{{\mathit{\boldsymbol{B}} \cdot \mathit{\boldsymbol{B}}}}{{2\mu }}$

(6)

$\mathit{\boldsymbol{J}} = \frac{1}{\mu }\nabla \times \mathit{\boldsymbol{B}}$

(7)

$\mathit{\boldsymbol{Q}} = - {\kappa _0}\nabla T - {\kappa _\parallel }\mathit{\boldsymbol{bb}} \cdot \nabla T$

(8)

$\mathit{\boldsymbol{ \boldsymbol{\varPi} }} = - {\nu _0}\nabla \mathit{\boldsymbol{\nu }} - 3{\nu _\parallel }\left( {\mathit{\boldsymbol{bb}} - \frac{1}{3}\mathit{\boldsymbol{I}}} \right)\left( {\mathit{\boldsymbol{bb}} - \frac{1}{3}\mathit{\boldsymbol{I}}} \right):\nabla \mathit{\boldsymbol{\nu }}$

(9)

式中:ρ为质量密度，v为速度矢量，p为压力，E为能量，B为磁场强度，μ为磁导率，Π为包含各向同性与各向异性部分的黏性应力张量，分别由系数ν₀和ν_‖表示。在各向异性的情况下，黏性通量限制在平行于磁场线方向。Q为热通量，同样包含各向同性与各向异性的情况，分别由系数κ₀和κ_‖控制。在各向异性的情况下，热通量同样限制在平行于磁场线方向。同时，感应方程包含多种非理想MHD效应，包括欧姆损耗(由电阻率η控制)、霍尔效应(由系数η_H控制)和双极性扩散(由系数η_AD控制)。

为了求解该方程，采用CTU+CT算法。CTU算法是为了求解守恒双曲系统的一种不分离的2D有限体积算法，它的3D变形随后由J.Saltzman^[22]提出。不幸的是，传统的12-solve CTU算法不能成为理想MHD方程的有效算法。这是可以理解的，而一般的理由指出算法的中间步骤使用基于守恒形式方程的维度分解的部分更新，这反过来又忽略了不同方向通量梯度的潜在平衡，即忽略了在MHD方程中始终存在约束 $\nabla$ ·B=0。因此，构建一个Godunov通量中网格单元边缘的平均电场的算法是必要的，这种算法通常称为CT算法。该算法可以描述为以Godunov电场为预测值、以CT电场为校正值的预测过程。T.A.Gardiner等^[23]提出了一个构建CT算法的简单框架，并构建和测试了一些CT算法。鉴于ε^c CT算法显示的最佳性能^[23]，本文中采用ε^c CT算法，ε^c CT算法包含逆风偏差，并为网格对齐的平行平面流动降低到正确的Godunov电动势。

为了使这一点更具体，注意到平行通量梯度条件(x-界面的x-通量梯度等)都包含在使用尺寸分割的PPM(piecewise parabolic method)界面状态算法中，这也是MHD方程的原始形式。与此同时，横向通量梯度条件包含在使用保守形式的方程中。由于在尺寸上分解的原始和守恒形式的MHD方程不相称，这相当于忽略了某些MHD源项导致的一阶准确的集成算法。此外，这种算法也显示了Gedanken实验中垂直于磁场回路的磁场分量的长期演化。这种相对简单的3D综合算法有二阶精度，且没有源项需要包含在垂直于磁场的分量的界面演化中，这也是和12-solve CTU算法相比的优势所在。同时，该算法对于网格对齐流动的相关对称问题可以恰好减少到二维CTU或是一维PPM算法。该算法的缺点是，该算法在Courant-Friedrichs-Lewy(CFL)因数 < 1/2时才是稳定的。因此，在很多情况下，6-solve和12-solve算法表现出类似的计算成本: 6-solve算法在CFL因数为1/2时的2个时间步几乎相当于12-solve算法在CFL因数为1时的一个时间步。

2. 计算域与初始条件

本文中采用加入等离子体的SF6气体团，爆炸初始气体团密度为150 kg/m³，比热比为1.09，高密度气体团气压为10.13 MPa; 外界空气密度为1.205 kg/m³，比热比为1.4，外界气压为101.3 kPa。同时，初始SF₆气体团的半径为0.1 m，黏性系数约为1.57×10^-5 Pa·s，整个计算域为边长为1.5 m的正方形。这里，采用1 200×1 200的方形网格，并且通过检验，满足网格收敛性条件。由于加入了等离子体，电导率、热导率以及双极扩散系数由下面的方法确定。

首先计算电导率，这里，电导率采用玻尔兹曼方程的Spitzer公式计算。它从碰撞积分的玻尔兹曼方程出发，在方程中忽略了短程碰撞的影响，考虑剩余的碰撞项，得到了在弱电场的情况下，离子和电子所组成的等离子体的电导率：

$\sigma = \frac{{2m{\gamma _0}}}{{3{{\rm{ \mathsf{ π} }}^{3/2}}{e^2}{j^3}\ln \wedge }}$

(10)

式中：m为等离子体质量，γ₀=0.1Z+0.49，e为电子电荷，Z为离子电荷; ${{j}^{2}}=\frac{m}{2kT}$ ，k为热力学常数，T为温度; ∧为库仑对数项。

如果碰撞积分到德拜长度为止，则库仑对数项为：

$\wedge = \frac{3}{{{{\rm{ \mathsf{ π} }}^{1/2}}n_{\rm{e}}^{1/2}{e^3}}}{\left( {\frac{{kT}}{2}} \right)^{3/2}}$

(11)

式中：n_e为电子数密度。根据式(10)计算，本文中初始电导率约为10⁶ S/m。

接下来计算热导率。当气体中相当部分的原子被电离，静电力将对粒子间碰撞起主要作用。若考虑远程作用，则离子碰撞将不能被视为双体方式。通过求解玻尔兹曼方程，可推得电离气体热导率的一阶近似：

${\mathit{k}_{\rm{c}}} = \frac{5}{2}\nu {c_\mathit{V}}$

(12)

式中：ν为一阶近似得到的黏度。若采用动力论有关比热的关系式， ${{c}_{V}}=\frac{{{N}_{0}}}{2}\frac{k}{{\bar{m}}}$ ，其中N₀为常数，m为等离子体平均质量。根据S.Chapman等^[24]的理论，可以导出k_c、ν的二阶近似。在作用力为平方反比作用律的时候，气体各个系数一阶近似值的误差就显著大于作用力指数较高时的误差，还有一种情况就是m₂/m₁为一个可以忽略的小量，它对应于离子-电子混合气，m₁为离子质量，m₂为电子质量。Landshoff研究过此种情况更高阶的近似，从他得到的数据表明：尽管二阶近似值与一阶近似值有显著的差异，但以后各阶近似所增加的修正量相对来说都是很微小的。在此假定黏度都是由重粒子引起的，电子黏度很小。由于电子的热运动速度远大于离子的热运动速度，因此热扩散效应主要是由电子引起的。经过估算，本文中，热导率约为1.2 W/(m·K)。

由于准电中性条件，电子和离子在等离子体内独立运动是不可能的，电子很快离开等离子体某一单元，不可避免地引起电场的出现，这个电场将进一步阻碍它们离开所在的体积单元，并促使离子更快地离去。这种扩散状态被称为双极型状态。由于电子的扩散系数远大于离子的扩散系数，因此可以得到^[25]：

${\eta _{{\rm{AD}}}} = {\eta _{\rm{i}}}\left( {1 + \frac{{{T_{\rm{e}}}}}{{{T_{\rm{i}}}}}} \right)$

(13)

式中：ηi为正离子热扩散系数，T_e、T_i分别为电子温度和离子温度。由此可见，双极扩散系数恒小于电子自由扩散系数，而大于离子扩散系数，因此双极性电场极大地影响了电子的定向速度。本文中，双极扩散系数取 $\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s。

3. 结果与讨论

在数值模拟过程中，为了更清楚地看到不稳定性，在爆炸界面处加入正弦微扰。对于没有加磁场的情况，当激波从重质气体传向轻质气体时，会产生向内传播的稀疏波和向外传播的透射激波，从而引起失稳，加速界面混合。激波作用不同密度气体界面后，界面演化一般经历3个阶段^[26]：(1)振幅线性增长阶段; (2)振幅非线性增长阶段; (3)湍流混合阶段。在膨胀初期，随机扰动在气体界面发生，扰动振幅在界面线性增长，随着膨胀的加剧，振幅经历非线性增长，气泡状和尖钉状结构出现，随着膨胀的结束以及球心处低压区的形成，圆球界面开始向球心回流，尖钉状结构破碎，不同流体开始形成湍流混合状态。下面给出外加磁场为0.05 T时，理想磁流体的演化情况，首先是密度分布图。

这里，加入方向为右上方45°的初始均匀强磁场并忽略电阻。从可以发现，与不加磁场的情况相比^[27]，由于磁场的抑制作用，垂直于磁场方向爆炸界面的尖钉状及气泡状结构明显减少，而平行于磁场方向的爆炸界面并没有受到明显影响。从可以看出磁场的作用好比垂直于磁场方向的压力，并阻碍界面处的轻质流体流入重质流体，反之亦然。为了更清楚地说明这一现象，将动量方程中与磁场相关的项拆开，分为磁压力 $\nabla$ (B·B/(2μ))和磁张力 $\nabla$ ·(BB/μ)，它们都是二阶张量。由于结果中磁压力的值比磁张力的值大很多，下面只讨论磁压力的变化情况。

图 1 理想情况下激波结构及流动界面的演化

Figure 1. The shock wave structure and the evolution of the interface under the ideal conditions

下载: 全尺寸图片幻灯片

首先，从公式角度进行分析。根据张量公式:

$\nabla \cdot \frac{{\mathit{\boldsymbol{BB}}}}{\mu } = \frac{1}{{2\mu }}\mathit{\boldsymbol{bb}} \cdot \nabla {B^2} + \frac{1}{\mu }{B^2}\mathit{\boldsymbol{\kappa }}$

(14)

式中：等号后第1项表示磁场能量在磁感线方向上的变化的投影，κ=b· $\nabla$ b=-b×( $\nabla$ ×b)为磁感线上某点的法向曲率，其大小为 $\left| \frac{\partial \mathit{\boldsymbol{b}}}{\partial {{l}_{\parallel }}} \right|$ ，方向(与切线斜交)指向法向曲率中心。

初始时刻，磁场强度处处相等，由磁感应方程(4)，开始时磁场在速度垂直于磁场的分量的梯度较大处变化较快，如垂直于磁场方向的爆炸界面处。因此，随着时间的推移，爆炸界面处的磁场强度越来越大，并妨碍轻、重质气体的相互混合。渐渐的，SF₆气体团受到y=-x方向磁压力的作用被压缩，空气和SF₆开始混合并沿着初始磁场方向形成尖钉状和气泡状结构，同时，磁压力作用于尖钉状结构及气泡状结构的横向部分，从而减缓其进一步演化。由公式(14)曲线坐标系下的形式可以看出，磁压力只作用于垂直于磁场的方向上，这也可以解释为什么气体团最终被压扁，而平行于初始磁场的方向并没有受到明显的影响。

由于初始磁场沿直线y=x方向，因此沿直线y=-x处受磁场影响最大。图 2显示了沿着y=-x的压力分布，与不加磁场的情况相比^[26]，发现在界面处出现小的突起，这一现象预示着磁场的作用。同时，与不加磁场的情况相比二次激波衰减速度更快。

图 2 在B₀=0.05 T, η=0, η_AD=0的情况下，不同时刻沿着y=-x的压力分布

Figure 2. Pressure distribution along y=-x at different times when B₀=0.05 T, η=0, η_AD=0

下载: 全尺寸图片幻灯片

图 3展示了不同时刻沿着y=-x的磁压力分布情况。初始时刻，磁压力很小且对称，随着透射激波的传播，界面在内部高压的作用下向外传播，因此，磁场最大值出现在边界上，同时，可以获得界面在y=-x方向的移动轨迹。可以看到界面从t=3.6 ms起开始向内移动，二次激波在t=6.9 ms于爆炸中心产生，并与爆炸界面在t=9.6 ms左右相遇。同样，从图中可以清楚地看到，二次激波传播速度非常快，同时在与界面交汇后，产生了巨大的磁压力。

图 3 在B₀=0.05 T, η=0, η_AD=

$\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s的情况下，不同时刻沿y=-x的磁压力分布

Figure 3. Magnetic pressure distribution along y=-x at different times when B₀=0.05 T, η=0, η_AD=

$\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s

下载: 全尺寸图片幻灯片

接下来研究加入电阻和双极扩散效应后的爆炸情况，首先是密度分布图。从图 4可以看出，与理想情况相比，磁场对气体团的压缩作用明显减弱了，在垂直于初始磁场的方向上出现了RM不稳定性现象。与之前的情况不同的是，此时RM不稳定性并没有出现尖钉状结构，不稳定性末端出现了球状突触，同时，在磁场的作用下，该滴状结构没有破裂，而是逐渐沿着垂直于初始磁场的方向被压缩，最终，不稳定性也得到了一定程度的控制。

图 4 B₀=0.05 T时非理想情况下激波结构及流动界面的演化

Figure 4. The shock wave structure and the evolution of the interface under the nonideal conditions when B₀=0.05 T

下载: 全尺寸图片幻灯片

下面是垂直于磁场方向的磁压力分布情况。从图 5可以看出，磁压力与之前相比明显减小，同时，在二次爆炸的过程中，磁压力的作用范围增大了，可以理解为电阻和双极扩散效应阻碍了磁场的作用，其张力减缓了界面突变，从而使边界处过渡更平缓。二次爆炸后，二次激波与界面交汇时受到磁压力的作用明显减小。由此也可以看出，在考虑电阻以及双极扩散效应时，磁场对不稳定性的抑制作用明显减弱了。

图 5 在B₀=0.05 T, η=10^-6 Ω·m, η_AD=

$\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s的情况下，不同时刻沿着y=-x的磁压力分布情况

Figure 5. Magnetic pressure distribution along y=-x at different times when B₀ = 0.05 T, η=10^-6 Ω·m, η_AD=

$\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s

下载: 全尺寸图片幻灯片

为了更清楚地看到2种作用分别对磁压力的影响，下面分别研究只考虑电阻和只考虑双极扩散效应的情况。由于除压缩程度不同外整个流场结构类似，因此接下来重点讨论垂直于初始磁场方向的磁压力分布图像。首先是加入电阻的情形，结果如图 6所示。从图 6可以看出，在二次爆炸前和二次爆炸后，界面处的磁压力减小了一半以上。因此，可以认为电阻对磁场也是起阻碍作用的。在二次爆炸前(图 6(a))，电阻不仅显著减小了磁压力，同时延缓了磁压力增大的速度以及界面向内移动的时间; 在二次爆炸时(图 6(b))，由于之前磁压力的减小产生的不稳定性现象，磁压力增大，外侧界面向内移动速度减缓; 在二次爆炸后(图 6(c))，磁压力再次得到了抑制，尤其是在激波与界面交汇处，并在距离爆炸中心0.5 m附近取得最大值。从图中也可以看出，电阻并没有增大磁压力的作用区域，仅从数值上减小了磁场对流体的作用强度。

图 6 在B₀=0.05 T, η=10^-6 Ω·m, η_AD=0的情况下，不同时刻沿着y=-x的磁压力分布

Figure 6. Magnetic pressure distribution along y=-x at different times when B₀=0.05 T, η=10^-6 Ω·m, η_AD=0

下载: 全尺寸图片幻灯片

下面是忽略电阻只考虑双极扩散效应的情况。从图 7可以看出，整体趋势和同时考虑电阻和双极扩散效应的情况类似，只是磁压力的大小比考虑电阻的情形要大。在二次爆炸前(图 7(a))，磁压力的值在t=1.2 ms后迅速减小，磁压力的作用区域在界面开始向内移动时逐渐增大(t=3.6 ms后); 在二次爆炸时(图 7(b))，磁压力进一步减小，磁压力作用区域继续增大; 在二次爆炸后(图 7(c))，磁压力同样在距离爆炸中心0.5 m附近取得最大值，由于忽略了电阻作用，该值和同时考虑电阻和双极扩散效应的结果相比有所增大，与只考虑电阻的情形相比还是大大减小了。综上，可以判断双极扩散效应不仅可以显著减小磁压力，还能加大磁压力的作用区域。

图 7 在B₀=0.05 T, η=0, η_AD=

$\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s的情况下，不同时刻沿着y=-x的磁压力分布情况

Figure 7. Magnetic pressure distribution along y=-x at different times when B₀=0.05 T, η=0, η_AD=

$\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s

下载: 全尺寸图片幻灯片

4. 结论

基于CTU+CT算法，对非理想磁场环境下球形重质气体的物理爆炸进行数值模拟，并与不考虑电阻以及双极扩散效应的情况进行对比, 得到以下结论：

(1) RM不稳定性会在重质气体爆炸的过程中产生，并会在剪切速度的影响下产生尖钉状及气泡状结构。同时，二次激波加速流场的湍流过度。在磁场作用下，RM不稳定性得到了有效的控制，尤其是沿着垂直于初始磁场的方向，整个高密度气体团会随着时间的推移逐渐被压缩。

(2) 在考虑电阻和双极扩散效应的情况下，磁场对气体团的作用明显减弱，不稳定性会逐渐产生。在爆炸后期，由于速度梯度的增加引起的界面磁场强度增大，不稳定性的发展再次受到磁压力控制，尖钉状与气泡状结构没有出现，取而代之的是液滴状结构。最终，液滴状结构在磁场的作用下沿垂直于初始磁场的方向被压缩，不稳定性得到控制。

(3) 在只考虑电阻的情况下，磁压力与理想情况相比有所减小。同时，它还延缓了磁压力增大的速度以及界面向内移动的时间。在这种情况下，气体团爆炸的整个过程与理想情况类似，与之前不同的是磁场对不稳定性的作用相对减弱，气体团被压缩的过程变慢。

(4) 在只考虑双极扩散效应的情况下，不仅磁压力的大小显著减弱，磁压力的作用区域还明显增大了。从中可以看出，双极扩散效应有明显的耗散作用，并且能减缓气体爆炸过程中的界面突变。在这种情况下，二次激波对气体团的作用大大减弱，气体团产生液滴状不稳定性结构。最终，随着液滴状结构的压缩，不稳定性得到控制。

图 1 理想情况下激波结构及流动界面的演化

Figure 1. The shock wave structure and the evolution of the interface under the ideal conditions

下载: 全尺寸图片幻灯片

图 2 在B₀=0.05 T, η=0, η_AD=0的情况下，不同时刻沿着y=-x的压力分布

Figure 2. Pressure distribution along y=-x at different times when B₀=0.05 T, η=0, η_AD=0

下载: 全尺寸图片幻灯片

图 3 在B₀=0.05 T, η=0, η_AD= $\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s的情况下，不同时刻沿y=-x的磁压力分布

Figure 3. Magnetic pressure distribution along y=-x at different times when B₀=0.05 T, η=0, η_AD= $\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s

下载: 全尺寸图片幻灯片

图 4 B₀=0.05 T时非理想情况下激波结构及流动界面的演化

Figure 4. The shock wave structure and the evolution of the interface under the nonideal conditions when B₀=0.05 T

下载: 全尺寸图片幻灯片

图 5 在B₀=0.05 T, η=10^-6 Ω·m, η_AD= $\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s的情况下，不同时刻沿着y=-x的磁压力分布情况

Figure 5. Magnetic pressure distribution along y=-x at different times when B₀ = 0.05 T, η=10^-6 Ω·m, η_AD= $\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s

下载: 全尺寸图片幻灯片

图 6 在B₀=0.05 T, η=10^-6 Ω·m, η_AD=0的情况下，不同时刻沿着y=-x的磁压力分布

Figure 6. Magnetic pressure distribution along y=-x at different times when B₀=0.05 T, η=10^-6 Ω·m, η_AD=0

下载: 全尺寸图片幻灯片

图 7 在B₀=0.05 T, η=0, η_AD= $\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s的情况下，不同时刻沿着y=-x的磁压力分布情况

Figure 7. Magnetic pressure distribution along y=-x at different times when B₀=0.05 T, η=0, η_AD= $\frac{1}{{{10}^{4}}\rho }$ m²·Pa/s

下载: 全尺寸图片幻灯片

参考文献(27)

[1]	Richtmyer R D. Taylor instability in shock acceleration of compressible fluids[J]. Communications on Pure and Applied Mathematics, 1954, 13(2):297-319. doi: 10.1002-cpa.3160130207/
[2]	Meshkov E E. Instability of a shock wave accelerated interface between two gases[J]. NASA Technical Translation F-13, 1970:74. http://d.old.wanfangdata.com.cn/Periodical/wlxb201723027
[3]	Cohen R H, Dannevik W P, Dimits A M, et al. Three-dimensional simulation of a Richtmyer-Meshkov instability with a two-scale initial perturbation[J]. Physics of Fluids, 2002, 14(10):3692-3709. doi: 10.1063/1.1504452
[4]	Yang Yumin, Zhang Qiang, Sharp D H. Small amplitude theory of Richtmyer-Meshkov instability[J]. Physics of Fluids, 1994, 6(5):1856-1873. doi: 10.1063/1.868245
[5]	Velikovich A L. Analytic theory of Richtmyer-Meshkov instability for the case of reflected rarefaction wave[J]. Physics of Fluids, 1996, 8(6):1666-1679. doi: 10.1063/1.868938
[6]	Velikovich A L. Richtmyer-Meshkov-like instabilities and early-time perturbation growth in laser targets and Z-pinch loads[C]//41st Annual Meeting of the Division of Plasma Physicsics Meeting. American Physical Society, 1999.
[7]	Brouillette M. The Richtmyer-Meshkov instability[J]. Annual Review of Fluid Mechanics, 2001, 34(34):445-468. http://d.old.wanfangdata.com.cn/Periodical/wlxb201723027
[8]	Sturtevant B. Rayleigh-Taylor instability in compressible fluids[J]. California Institute of Technology Report, 1989:92. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_cond-mat%2f0402318
[9]	Zoldi C A. A numerical and experimental study of a shock-accelerated heavy gas cylinder[D]. New York: State University of New York, 2002.
[10]	Layes G, Métayer O L. Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[J]. Physics of Fluids, 2007, 19(4):501-100. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=debe921b528f2bf99453d0abeabfbad9
[11]	邹立勇, 刘金宏, 谭多望, 等.弱激波冲击无膜重气柱和气帘界面的实验研究[J].高压物理学报, 2010, 24(4):241-247. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201001828991 Zou Liyong, Liu Jinhong, Tan Duowang, et al. Experimental study on the membraneless heavy gas cylinder and gas curtain interfaces impacted by a weak shock wave[J]. Chinese Journal of High Pressure Physics, 2010, 24(4):241-247. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201001828991
[12]	Zhang Qiang, Graham M J. A numerical study of Richtmyer-Meshkov instability driven by cylindrical shocks[J]. Physics of Fluids, 1998, 10(4):974-992. doi: 10.1063/1.869624
[13]	Zhang Qiang, Graham M J. Scaling laws for unstable interfaces driven by strong shocks in cylindrical geometry[J]. Physical Review Letters, 1997, 79(14):2674-2677. doi: 10.1103/PhysRevLett.79.2674
[14]	Dutta S, Glimm J, Grove J W, et al. Spherical Richtmyer-Meshkov instability for axisymmetric flow[J]. Mathematics and Computers in Simulation, 2004, 65(4/5):417-430. http://www.sciencedirect.com/science/article/pii/S0378475404000266
[15]	Glimm J, Grove J, Zhang Y, et al. Numerical study of axisymmetric richtmyer-meshkov instability and azimuthal effect on spherical mixing[J]. Journal of Statistical Physics, 2002, 107(1):241-260. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=cac9e80fe96ff73075d3d6ce99a8e6c3
[16]	Holmes R L, Dimonte G, Fryxell B, et al. Richtmyer-Meshkov instability growth: Experiment, simulation and theory[J]. Journal of Fluid Mechanics, 1999, 389(4):55-79. http://d.old.wanfangdata.com.cn/NSTLQK/NSTL_QKJJ0211419717/
[17]	Lindl J D, Mccrory R L, Campbell E M. Progress toward ignition and burn propagation in inertial confinement fusion[J]. Physics Today, 1992, 45(9):32-40. doi: 10.1063/1.881318
[18]	Arnett D. The role of mixing in astrophysics[J]. Physics, 2000, 127(2):213. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_1012.1925
[19]	Samtaney R. Suppression of the Richtmyer-Meshkov instability in the presence of a magnetic field[J]. Physics of Fluids, 2003, 15(8):53-56. doi: 10.1063/1.1591188
[20]	Wheatley V, Samtaney R, Pullin D I. The Richtmyer-Meshkov instability in magnetohydrodynamics[J]. Physics of Fluids, 2009, 21(8):445. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_1206.2087
[21]	胡希伟.等离子体理论基础[M].北京:北京大学出版社, 2006:6-9.
[22]	Saltzman J. An unsplit 3D upwind method for hyperbolic conservation laws[J]. Journal of Computational Physics, 1994, 115(1):153-168. doi: 10.1006/jcph.1994.1184
[23]	Gardiner T A, Stone J M. An unsplit Godunov method for ideal MHD via constrained transport in three dimensions[J]. Journal of Computational Physics, 2007, 205(2):509-539. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_0712.2634
[24]	Chapman S, Cowling T G. Non-uniform gases: Scientific books: The mathematical theory of non-uniform gases[J]. Science, 1940, 91(2372):575. doi: 10.1126/science.91.2372.575
[25]	Brown S C, Holt E H. Introduction to electrical discharges in gases[J]. American Journal of Physics, 1968, 36(9):854-854. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_0711.1672
[26]	薛大文, 陈志华, 韩珺礼.球形重质气体物理爆炸特性[J].爆炸与冲击, 2014, 34(6):759-763. http://d.old.wanfangdata.com.cn/Periodical/bzycj201406017 Xue Dawen, Chen Zhihua, Han Junli. Physical characteristics of circular heavy gas cloud explosion[J]. Explosion and Shock Waves, 2014, 34(6):759-763. http://d.old.wanfangdata.com.cn/Periodical/bzycj201406017
[27]	Zhuang Zhihong, Xue Dawen, Chen Zhihua, et al. The eruption of a high-pressure cylindrical heavy gas cloud[J]. Canadian Journal of Physics, 2013, 91(10):850-854. doi: 10.1139/cjp-2013-0014