不同初始磁场对激波冲击R22重气柱过程影响的数值模拟*

2018-03-07 04:56林震亚郭则庆张焕好陈志华

爆炸与冲击 2018年2期

林震亚,郭则庆,张焕好,陈志华,刘迎

(南京理工大学瞬态物理国家重点实验室，江苏南京 210094)

在激波穿过物质界面过程中，由于激波区域存在较强的压力梯度，而物质界面区域存在较强的密度梯度，两者的不平衡使激波作用下的物质界面失稳形成不同尺度的涡，并最终向湍流转捩，从而增强了界面两侧流体的混合[1]，同时在界面处还会出现Richtmyer-Meshkov(RM)不稳定与Kelvin-Helmholtz(KH)不稳定现象，因而激波与界面相互作用是一个复杂的多尺度、强非线性的物理问题。此类现象广泛存在于超燃冲压发动机的燃料混合、惯性约束核聚变、水中炸药爆炸以及天体物理中的超新星爆炸等领域，因此有着重要而广泛的研究背景。

由于激波-界面相互作用问题在学术和工程领域的研究价值，对激波与不同密度的气体界面相互作用问题有了大量研究。Haas等[2]对弱激波与R22重气柱、气泡以及He轻气柱、气泡的作用过程进行了实验研究，对气柱、气泡变形过程进行了分析讨论。Tomkins等[3]对激波与SF6重气柱的作用过程进行了实验研究，研究分析了两种介质的混合机理。范美如等[4-5]对激波与矩形、椭圆、菱形以及两种三角形五种不同形状的SF6气柱作用过程进行了数值研究，对比分析了这几种形状界面的波系、涡量以及气体界面的演变。王显圣等[6]数值研究了入射激波以及反射激波与SF6重气柱作用过程。Si等[7]利用高速纹影技术对平面入射激波以及反射激波与He轻气泡以及SF6重气泡作用过程进行了实验研究，结果清晰显示了激波诱导气泡变形的过程，并对整个过程中气泡尺寸变化以及流场环量变化进行了定量分析。

非理想情况下关于RM不稳定性控制方法方面的研究尚处于起步阶段，研究表明，合适的磁场对气体界面RM不稳定性的产生和发展具有一定的抑制作用，研究仍主要集中在理论研究方面。Chandrasekhar[8]和Wheatley等[9-10]采用线性理论方法，推导了垂直或平行于界面的磁场对不可压缩流体界面不稳定性的影响，得出平行于物质界面的磁场能够抑制Rayleigh-Taylor(RT)和RM不稳定性的结论。Cao等[11]同时考虑了初始平行于界面的磁场和剪切流对界面不稳定性的影响，发现磁场能够抑制界面不稳定性，而剪切流将会促进界面不稳定性的发展，两者之间存在竞争机制。由于线性理论所得到的平面流场中的RM不稳定性与垂直于平面的磁场无关，这与实际不符，为了研究此种情况下磁场的效应，Khan等[12]在势流理论的基础上，建立了新的模型，分析发现在非线性效应下，磁场能够影响界面不稳定性的发展。Shin等[13]对不同方向磁场中平面激波与球形气泡相互作用问题进行了数值模拟，发现在界面发展初期，磁场效应很小，而到了发展后期，即使是弱磁场，对界面形状的影响也很明显，磁场越强，界面内外物质混合率越低，后期湍流发展被抑制程度越强。李源等[14]理论分析了磁场中非理想流体的RT不稳定性气泡的演化过程，在与磁场垂直的平面中，综合考虑黏性和表面张力的影响，推导了描述二维非理想磁流体RT不稳定性气泡运动的控制方程组，分析了流体黏性、表面张力和磁场对气泡界面不稳定性发展的影响。

目前，关于RM不稳定性的研究主要用于理论和机理分析，由于高维可压方程没有解析解，因此理论研究主要针对不可压缩MHD方程；对于机理分析，不论是实验还是数值模拟，目前的分析重心在激波、阿尔文波等物理现象上。由于非理想情况下无法计算阿尔文波，同时激波与磁场并无直接关联，因此本文中侧重于分析涡量、磁压力与磁能量。本文中，基于MHD方程，采用CTU+CT算法，以涡量为切入点讨论不同磁感应强度对激波冲击R22重质气柱过程中界面不稳定性的影响，揭示磁场对不稳定性控制的机理。

1 数值方法与计算模型

1.1 数值方法

采用非理想磁流体动力(MHD)方程组，模拟激波冲击R22气柱的作用过程，具体参见文献[15]，并选用6解CTU+CT算法[16-17]对MHD方程组进行求解。为了清晰地说明程序运行方式，下面对6解CTU+CT算法的主要过程进行描述。

(1)

在y和z方向上的流量也有相似的表达形式。

第3步，在每个界面通过δt/2步长的横向流量梯度生成PPM界面状态，其中流体动力学变量(质量、动量及能量密度)为：

(2)

其中，x方向的MHD动量密度源项与能量密度源项分别为：

(3)

(4)

磁场成份在预测的流量区域通过CT电场生成，磁场界面组份通过斯托克斯循环的积分形式形成：

(5)

磁场y、z方向的组分也可由类似的方法得到。这些磁场横向组份中的MHD源项起点可通过前文描述的感应方程的方向拆分处理。另外，动量与能量密度的MHD源项通过MHD方程的原始变量形式来计算PPM界面状态。同样，y、z方向的界面状态也由类似的方式得到，只要按照上面的描述改变(x,y,z)和(i,j,k)的周期排列。

第4步，对于每个通过δt/2步长更新的界面状态，计算相关的流量，给出x方向的界面流量：

(6)

y和z方向的界面流量也有相似的表达方式。

(7)

第6步，从n至n+1步更新结果，流体动力学变量(质量、动量和能量密度)通过标准流量积分关系进行迭代：

(8)

同时，磁场的界面平均组份通过斯托克斯循环积分进行迭代，x方向组份为：

(9)

y和z方向的组分也有相似的表达形式。

这相对简单的3D综合算法具有二阶精度，且在磁场界面通常组份的演化过程中不包含源项。这种算法对于均匀网格的流动，并包含相关对称性的问题，可降为2D CTU和1D PPM综合算法。从实验中观察到的下降趋势，对于该算法的CFL<1/2时稳定。

1.2 计算模型

图1为激波与R22气柱相互作用的二维计算模型，为了便于分析不同强度磁场对不稳定性的影响，选取与文献[15]相同的计算域尺寸、气体参数、初始条件、边界条件以及网格划分，并讨论了磁感应强度分别为0.01与0.05 T(即β=2 500，100，其中β是磁感应强度的量纲一量，β-1=B2/(2μ0p0))时对流场的影响。另外，为了使气体受磁场影响，先将气体进行电离。由文献[15]，其初始电导率取106S/m，热导率取1.4 W/m2，黏性系数取3.72×10-5Pa·s，霍尔系数取6×10-6m3/C，同时，双极扩散系数D随粒子数密度ρ的增加而减小，因此双极扩散系数取ρD=10-4m2·Pa/s。

2 结果与讨论

图2为无磁场时入射激波与R22气柱相互作用过程的涡量(上)与密度纹影图(下)。由图可见，计算结果中激波的反射、绕射及发展过程均与数值模拟结果[18]完全相符，且两者中气柱的变形过程相似。

当入射激波与气柱碰撞后，分别形成一道向上游传播的反射激波与一道向气柱内部传播的圆弧透射激波(见图2(a))。在激波作用下，气柱在x轴方向上被不断压缩。当入射激波绕过气柱上下两侧圆弧顶点后出现弯曲，形成绕射激波(见图2(b))。随后，入射激波在气柱尾部聚焦，形成向下游传播的圆弧二次激波，最终在气柱尾部聚焦形成局部高压区并产生射流(见图2(c)～(f))，此过程与实验结果[19]相符。另外，由涡量分布可见，当入射激波经过气柱界面后，诱导界面发生RM不稳定而出现涡量(见图2(a))，且在激波绕经界面上下两顶点时，引起该处界面的剧烈失稳而卷起形成小涡串(见图2(b))。随后，上下界面处反射形成的反射激波又与气柱的上下界面发生碰撞，加剧了界面的失稳而形成串级旋涡序列。当绕射激波在中心轴上聚焦并反射后，与气柱持续作用形成两个主涡(见图2(c)～(e))。随着激波在流场中的来回反射，使流场的波系结构复杂多变，并与气柱发生相互作用，最终使气柱界面完全失稳(见图2(f))。后期，随着界面左侧顶端扰动的发展，逐渐形成尖钉状与气泡状结构。

图3为施加磁感应强度B=0.01 T后，激波与R22气柱相互作用过程中的涡量分布图。可见，由于磁场及热传导作用，使气柱边界在与入射激波作用后会出现内外两个涡层，此现象在初始磁场垂直于来流方向(垂直磁场)时尤为明显。对于垂直磁场，气柱界面保持稳定且无涡量卷起，此时涡量均匀分布于入射激波作用过的界面处，随后内侧涡层逐渐衰减，外侧涡层的涡量较大区域出现在尾部，最终在垂直磁场作用下界面的卷起得到抑制，阻止了不稳定性结构的出现。而对于初始磁场平行于来流方向(水平磁场)时，气柱界面在上下顶端(即剪切力最大处)开始卷起，并在t=250 μs时形成4个涡。随后，气柱上下界面处的两个涡逐渐混合，最终在尾部形成两个主涡。由此可知，磁感应强度0.01 T的水平磁场虽未能抑制尾部主涡的形成，但却很好抑制了界面处次级涡的出现。

图4为B=0.05 T时的涡量分布。可知，当磁感应强度加大时，涡层分离得更快，且不再依附于气柱界面上。对于垂直磁场情况，初始时的涡量分布(见图4(a)～(c))与B=0.01 T时相似(见图3(a)～(c))，但由于两涡层出现后，外侧涡层在入射激波及透射激波先后聚焦后向右运动，同时沿界面法向方向向外扩张，形成斜向后的块状结构，最终与上下壁面碰撞并反射，如图4(d)～(f)所示。同时，内侧涡层则沿y轴向内压缩，最终变为椭圆形。此时，由于界面处涡量较小，无不稳定结构出现，因此不稳定性得到了有效的抑制。而对于水平磁场情况，物质界面上的涡量分布情况与B=0.01 T时相似。此时，气柱上下两端顶点处涡量最大，不同的是气柱前端涡层迅速分离，形成两个涡层。其中，外侧涡层右端逐渐向气柱尾端移动，逐渐包裹气柱，而内侧涡层尾部则随着射流向右运动，最终将内侧涡层拉长。在水平磁场作用下，后期气柱尾部虽有卷起，但并未形成主涡(见图4(f))，同时很好抑制了界面处的不稳定性。由此可见，当磁感应强度增大时，会加速内外两个涡层分开，并在随后的运动中脱离气柱界面，从而更好地抑制激波与界面作用过程中的不稳定性。

图5是t=945 μs时不同初始磁场情况下的密度纹影图。可见，当磁感应强度增大时，对RM不稳定性的控制效果更好。同时，垂直磁场对不稳定性的控制效果优于平行磁场，它不仅阻止了主涡的卷起以及次级涡的生成，同时对界面也有更好的压缩作用，因此对RM及KH不稳定性起到有效的抑制作用。此外，在磁场较小时，涡层附着于界面上，而随着磁场的增大，涡层逐渐与界面脱离，从而有效地控制界面处的涡量形成。

表1为不同初始条件下不同时刻的平均涡量。当B=0.05,0.01 T时，平均涡量在透射激波聚焦前逐渐增大，在聚焦时达到局部最大值后略微减小(t=200 μs)，随后开始继续增大(t=450 μs)，其增速随时间而变大。当B=0.05 T时，平均涡量始终增大，且增速随时间也逐渐增大。然而，在聚焦前，垂直磁场的平均涡量略大于水平磁场，但在聚焦后则刚好相反。此外，随着磁感应强度增大，虽然能有效抑制不稳定性，但平均涡量明显大于磁场较小的，甚至大于不加磁场的。由此可知，磁场并不能控制涡量的生成，但可以控制涡量的位置，并抑制界面的卷起。

表1 流场平均涡量Table 1 Average vorticity

为了更清楚地说明涡量变化情况，下面将采用涡量的平方进行分析，它能反映流团的旋转能量，称为渦度拟能[20]。表2为不同初始条件下不同时刻的平均涡度拟能。在不施加磁场(B=0)时，平均涡度拟能在t=200 μs内逐渐增大。当入射激波在尾部聚焦后，由于流场进入湍流转捩阶段，涡结构逐渐扩散，因此平均涡度拟能逐渐减小，并最终稳定在1.1×107s-2附近。当B=0.01 T时，平均涡度拟能明显减小，同时垂直磁场比平行磁场对涡度拟能具有更好的抑制作用。由表2可知，虽然加入磁场后平均涡度拟能变化趋势与无磁场情况类似，但随着发展后期流场的演化，加剧界面处的扰动，使平均涡度拟能逐渐增加，因而在t=945 μs时，平均涡度拟能分别增至4.533×106和9.524×106s-2。当B=0.05 T时，平均涡度拟能受到更明显的控制，在激波与气柱作用过程的前中期，平行与垂直磁场情况下的平均涡度拟能均小于B=0.01 T时的。由此可见，平均涡度拟能清晰地反映磁场对不稳定性的控制。然而，在演化后期，虽然流场进入湍流转捩阶段，但单涡由于耗散作用逐渐减弱，使平均涡度拟能小于发展前中期，因此涡度拟能并不能反映流场的无序情况。

文献[15]表明，当B=0.01 T时，磁能量与磁压力较大区域分布在界面附近，且气柱上下顶点处的磁压力及磁能量随着时间不断增大。图6～7分别为B=0.05 T时垂直与水行磁场作用下的磁压力与磁能量分布曲线。对于垂直磁场，磁压力和磁能量均增大，而当透射激波在水平轴线上发生聚焦后则稍有减小，随后，磁压力和磁能量随着界面的演化继续增大。同时，因磁场较小时涡层附着于界面，而磁场较大时涡层与界面分离并逐渐远离，因此磁压力及磁能量较大的区域位于涡层处而不是界面附近。由此可知，流场中涡量较大的位置磁压力及磁能量较大，这也与前面的结果相符。

对于水平磁场，由于初始时刻磁场曲率较小，因此来流方向上的磁压力较小。当气柱界面处由于剪切力作用开始卷起后，垂直于来流方向的磁压力对界面KH不稳定性的继续发展起到了明显的抑制作用。在这种情况下，磁压力及磁能量随时间不断增大。对于纵向磁场，磁压力在气柱界面处快速脉动，上下界面相同横坐标处磁压力并不相同，且数值相差较大，与之相对，磁能量大小基本呈上下对称分布。

3 结论

基于MHD方程及CTU+CT算法，对不同初始磁场下的激波冲击重质气柱过程进行了数值模拟，揭示了磁场抑制界面不稳定性的机理，得到以下结论。

(1) 垂直磁场及水平磁场对不稳定性均具有很好的抑制作用，其中水平磁场能很好地抑制界面的不稳定性，却不能抑制尾部主涡的形成，而垂直磁场则能同时做到以上两点。此外，随着磁感应强度的增大，对不稳定性的抑制效果会更好。

(2) 与不加磁场情况相比，加入磁场后平均涡量并未明显减小。当B=0.01 T时，平均涡量在透射激波聚焦前逐渐增大，在聚焦时达到局部最大值，随后略微减小，之后继续增大，且增速随时间呈增大趋势，可见B=0.01 T时，平均涡量略小于不加磁场情况。当B=0.05 T时，平均涡量始终增大，且增速随时间也逐渐增大。此外，当磁感应强度增大时，虽然能有效抑制界面不稳定性，但平均涡量明显大于磁场较小(B=0.01 T)的情况，甚至大于不加磁场的情形。因而磁场虽不能控制涡量的生成，但可以控制涡量的位置，并抑制界面的卷起。

(3) 平均涡度拟能可较好地反映出磁场对不稳定性的控制效果。当施加磁场后，平均涡度拟能明显小于无磁场情况，且随着磁场的增大而减小。同时，垂直磁场比平行磁场更能降低平均涡度拟能，这与不同方向初始磁场对不稳定性的控制效果相符。此外，由于涡的强度比涡的尺度对平均涡度拟能影响更大，因而平均涡量比平均涡度拟能可以更好地反映流场无序情况。

(4) 对于垂直磁场，磁场较小时的磁压力与磁能量随时间不断增大，而磁场较大时的磁压力与磁能量则在开始时较大，但在透射激波聚焦后出现减小，随后继续增大。对于水平磁场，由于初始时刻垂直于磁场方向的速度较小，此时磁压力作用不明显。随着气柱界面的卷起，磁压力迅速增大并抑制不稳定性的发展，在此情况下，磁压力及磁能量随着时间始终不断增大。此外，当磁场较小时，涡层附着于界面，涡层与界面随着磁场增大而出现分离，并随时间不断远离界面。

[1] BROUILLETTE M. The Richtmyer-Meshkov instability[J]. Annual Review of Fluid Mechanics, 2002,34:445-468.

[2] HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J]. Journal of Fluid Mechanics, 1987,181:41-76.

[3] TOMKINS C, KUMAR S, ORLICZ G, et al. An experimental investigation of mixing mechanisms in shock-accelerated flow[J]. Journal of Fluid Mechanics, 2008,611:131-150.

[4] 范美如,翟志刚,司廷,等.激波与不同形状界面作用的数值模拟[J].中国科学:物理学、力学、天文学,2011,41(7):862-869.

FAN Meiru, ZHAI Zhigang, SI Ting, et al. Numerical simulation of interaction with different shape accelerated by a planar shock[J]. Scientia Sinica: Physica, Mechanica & Astronomica, 2011,41(7):862-869.

[5] FAN M R, ZHAI Z G, SI T, et al. Numerical study on the evolution of the shock-accelerated SF 6 interface: Influence of the interface shape[J]. Science China: Physics, Mechanics & Astronomy, 2012,55(2):284-296.

[6] 王显圣,司廷,罗喜胜,等.反射激波冲击重气柱的RM不稳定性数值研究[J].力学学报,2012,44(4):664-672.

WANG Xiansheng, SI Ting, LUO Xisheng, et al. Numerical study on the RM instability of a heavy-gas cylinder interacted with reshock[J]. Chinese Journal of Theoretical and Applied Mechanics, 2012,44(4):664-672.

[7] SI T, ZHAI Z, YANG J, et al. Experimental investigation of reshocked spherical gas interfaces[J]. Physics of Fluids, 2012,24(5):054101.

[8] CHANDRASEKHAR S. Hydrodynamic and hydromagnetic stability[M]. Courier Dover Publications, 2013:441-453.

[9] WHEATLEY V, PULLIN D I, SAMTANEY R. Stability of an impulsively accelerated density interface in magnetohydrodynamics[J]. Physical Review Letters, 2005,95(12):125002.

[10] WHEATLEY V, SAMTANEY R, PULLIN D I. The magnetohydrodynamic Richtmyer-Meshkov instability: The transverse field case[C]∥The 18th Australasian Fluid Mechanics Conference. Australasian Fluid Mechanics Society, 2012.

[11] CAO J, WU Z, REN H, et al. Effects of shear flow and transverse magnetic field on Richtmyer-Meshkov instability[J]. Physics of Plasmas, 2008,15(4):042102.

[12] KHAN M, MANDAL L, BANERJEE R, et al. Development of Richtmyer-Meshkov and Rayleigh-Taylor instability in the presence of magnetic field[J]. Nuclear Instruments and Methods in Physics Research Section: Accelerators, Spectrometers, Detectors and Associated Equipment, 2011,653(1):2-6.

[13] SHIN M S, STONE J M, SNYDER G F. The magnetohydrodynamics of shock-cloud interaction in three dimensions[J]. The Astrophysical Journal, 2008,680(1):336-348.

[14] 李源,罗喜胜.黏性、表面张力和磁场对Rayleigh-Taylor不稳定性气泡演化影响的理论分析[J].物理学报,2014,63(8):277-285.

LI Yuan, LUO Xisheng. Theoretical analysis of effects of viscosity, surface tension, and magnetic field on the bubble evolution of Rayleigh-Taylor instability[J]. Acta Physica Sinica, 2014,63(8):277-285.

[15] 林震亚,张焕好,陈志华,等.磁场对激波冲击R22重气柱作用过程影响的数值模拟[J].爆炸与冲击,2017,37(4):748-758.

LIN Zhenya, ZHANG Huanhao, CHEN Zhihua, et al. Influence of magnetic field on interaction of shock wave with R22 heavy gas column[J]. Explosion and Shock Waves, 2017,37(4):748-758.

[16] SALTZMAN J. An unsplit 3D upwind method for hyperbolic conservation laws[J]. Journal of Computational Physics, 1994,115(1):153-168.

[17] GARDINER T A, STONE J M. An unsplit Godunov method for ideal MHD via constrained transport in three dimensions[J]. Journal of Computational Physics, 2008,227(8):4123-4141.

[18] 沙莎,陈志华,薛大文.激波冲击 R22 重气柱所导致的射流与混合研究[J].物理学报,2013,62(14):291-300.

SHA Sha, CHEN Zhihua, XUE Dawen. The generation of jet and mixing induced by the interaction of shock wave with R22 cylinder[J]. Acta Physica Sinica, 2013,62(14):291-300.

[19] HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J]. Journal of Fluid Mechanics, 1987,181:41-76.