激波与轻质气柱作用过程的磁场抑制特性*

张升博 张焕好† 张军 毛勇建 陈志华‡ 石启陈 郑纯

1) (南京理工大学,瞬态物理全国重点实验室,南京 210094)

2) (中国工程物理研究院,总体工程研究所,绵阳 621999)

3) (南京理工大学,能源与动力学院,南京 210094)

本文采用CTU (corner transport upwind)+CT (constrained transport)算法求解理想可压缩磁流体动力学(magneto-hydro-dynamic,MHD)方程,仿真研究了不同方向磁场控制下高斯分布轻质气柱界面受平面冲击波扰动后的演化过程,揭示了磁场方向对界面不稳定性的影响机理.仿真结果探讨了有/无磁场作用下流场特性与波系结构的发展,对比分析了磁场方向对气柱的长度、高度、射流宽度和体积压缩率的影响,并结合流场上半区环量、能量分量、速度和磁场力分布,多角度分析了磁场方向对界面不稳定性的影响机理.结果表明,磁压力推动涡量远离界面,降低了涡量在密度界面上的沉积而附着在分裂后的涡层上,从而有效抑制Richtmyer-Meshkov 不稳定性对界面的影响;由于磁张力附着在被分离的涡层上,且其作用方向与界面因速度剪切而卷起涡的方向相反,因此抑制了界面因Kelvin-Helmholtz 不稳定性而形成涡串.另外,纵向磁场控制下的磁张力反作用于中轴射流方向,同样抑制了Rayleigh-Taylor 不稳定性的发展.

1 引言

当激波与不同密度流体界面相互作用时,会引起一系列复杂的流动现象,包括激波的反射、常规与非常规折射以及激波衍射等.同时,界面处还会出现Richtmyer-Meshkov 不稳定性(Richtmyer-Meshkov instability,RMI)[1,2]现象.受RMI 扰动的界面,通常会因界面两侧速度不同导致Kelvin-Helmholtz 不稳定性(Kelvin-Helmholtz instability,KHI)[3,4]的产生,并且伴随着Rayleigh-Taylor不稳定(Rayleigh-Taylor instability,RTI)[5,6]现象.界面不稳定性诱导界面扭曲变形,增强了界面两侧流体的混合[7].此现象广泛存在于爆炸物理学、超燃冲压发动机、惯性约束核聚变、以及天体物理学等领域[8].因此对其深入研究有助于理解和预测流体相互作用中的复杂现象,对相关工程应用和技术的优化具有重要学术意义和研究价值.

自RMI 被提出以来,国内外学者们对激波与不同形状、不同密度分布界面的相互过程开展了大量的实验[9-20]与仿真研究[7,8,21-24].研究表明,在某些科研领域和工程应用中,界面不稳定性具有积极的作用,如燃料与助燃剂的混合,界面不稳定性可提高燃料的燃烧效率;而在另一些领域中则是有害的甚至是危险的,如惯性约束核聚变(inertial confinement fusion,ICF)等应用中,RMI 的存在可能导致湍流混合,使不同材料之间发生混合,影响聚变反应的正常进行,这会降低聚变产额,使RMI 成为不利因素,需对其进行控制[7].因此,在工程应用领域中开展对控制RMI 的研究具有重要的价值.

在日常生活中,物质大多呈固、液、气三种状态,但在高温高压工作环境中,物质大多呈第四态,即等离子体状态.等离子体是一种很好的导电体,易受磁场的影响,而如何利用磁场来移动或加速等离子体,进而对界面不稳定性进行宏观控制,成为近年研究的热点和难点[25].现有对界面RMI 控制方法的研究发现,合适的磁场可抑制气体界面RMI 的产生,但仍处在理论研究阶段.此外,因磁流体实验研究难度大且危险系数高,国内外鲜有关于流体界面不稳定性的实验研究.因此,通过数值模拟求解磁流体动力学(magneto-hydro-dynamic,MHD)方程组[25-32]是当下对于磁控流体界面RMI研究最为可靠的方法.

2003 年,Samtaney[26]对理想磁流体的二维数值模拟表明,磁场的存在抑制了RMI 的增长,其抑制机制归因于激波引起的涡量前端Alfvén 波向远离密度界面传输.Wheatley 等[27]对不同方向磁场抑制RMI 的机理进行理论和数值研究,发现弱磁场条件下的数值结果与不可压线性理论符合较好,且密度界面与激波作用后,界面上出现的涡被分解并附着在平行和反平行于磁场传播的波上.Sano等[28]数值研究了强激波与单模界面作用过程,发现不同方向磁场控制下的RMI 会存在拉伸和增强磁场的现象,这是超新星遗迹中激波位置处形成强磁场的原因.李源和罗喜胜[29]数值模拟了磁场条件下平面激波与矩形气柱的作用过程,分析了矩形密度界面在不同磁场中的演化,结果表明磁场降低了涡量在界面上的生成.邱志勇等[30]线性分析了磁场和黏性影响下的RMI 及单模正弦扰动,发现流向磁场和黏性均起到稳定作用.最近,Tapinou等[31,32]采用五矩多流体(multi-fluid plasma,MFP)等离子体模型研究了外加磁场对碰撞等离子体RMI 的影响,发现外加磁场通过减缓关键因素(如斜压、磁场扭矩和碰撞项)在环流演化中的作用,有效降低界面扰动的生长速率.除此之外,粒子间碰撞引发的介尺度动理学行为可能对ICF 点火产生重要影响[33].在这方面,理论和模拟研究面临如下两难困境: 宏观连续模型物理功能不足、微观分子动力学模拟在尺度上无能为力;同时,相关复杂物理场分析技术严重不足,导致数据信息的严重浪费.为此,Xu 等[34-38]提出和发展了离散玻尔兹曼建模与分析方法(discrete Boltzmann modelling and analysis method,DBM),不仅关注模拟前的物理建模,而且关注模拟后的复杂物理场分析.近期,Song 等[38]发展了等离子体系统的DBM 建模和分析方法,研究了磁场对Orszag-Tang (OT) 涡问题和RMI 中的宏观流动非平衡(hydrodynamic non-equilibrium,HNE)行为和与之伴随的热力学非平衡 (thermodynamic non-equilibrium,TNE)行为.结果表明在具有外部磁场的情况下,磁场对RMI 演化具有抑制效应,而全局平均非平衡强度和热流熵产生率可用来判断抑制RMI 的临界磁场强度.

前人的工作旨在研究磁场对RMI 的抑制作用,认为磁控条件下涡量分层现象是界面致稳的原因,但先前研究仍不能完全阐明洛伦兹力对界面的作用过程,且并未说明磁场作用对流场能量分布的影响.先前,笔者团队[25]研究了不同强度流向磁场对RMI 的抑制作用,并结合磁能量分析了磁张力对界面不稳定性的影响;还分析了不同组分分布对界面不稳定性和流场特性的影响机理[24].另外,本课题组[39,40]研究了有/无磁场控制下激波冲击高斯分布重质气柱,得到圆形界面的演化过程及洛伦兹力对RMI 的抑制作用.然而,上述研究均未考虑磁场作用下激波与高斯分布轻质气柱的作用及洛伦兹力对轻质界面不稳定性的作用机理,而气柱界面密度分布规律的改变会影响激波冲击下界面不稳定性的发展及气柱形态的演化.其中,激波作用后,重质气柱内的气流速度小于外侧,而轻质气柱内的气流速度则高于外侧,使速度梯度诱导的界面形态演化过程存在差别;此次,轻质气柱受激波冲击后,气柱外侧重质气体会涌入轻质气柱内,诱导形成RTI 现象,同样导致主涡结构的演化与重质气柱存在明显的差异;因此,激波冲击下重质和轻质气柱形态结构演化的差异必将影响磁场对界面不稳定性的抑制效果.此外,在实际ICF 过程中,激波与密度界面作用过程复杂多变,同时存在激波与轻质/重质界面的相互作用,因而研究激波与轻质气柱的作用过程对ICF 实际应用具有重要的指导意义.

基于此,本文对激波冲击高斯分布轻质气柱界面进行数值研究,为保证磁场散度为零,采用CTU(corner transport upwind)+CT(constrained transport)算法[41,42]求解理想的MHD 方程组,对有/无磁场及不同方向磁场下的平面激波冲击高斯分布氦气气柱进行研究,详细分析磁场控制下流场的波系结构与界面的演化过程,定量分析了流场环量和气柱特征尺度的发展规律,并结合流场能量、磁能量、磁场梯度、磁张力和磁压力详细探讨不同方向磁场控制对界面不稳定性的影响规律.

2 计算方法和计算模型

2.1 物理模型和数值方法

基于可压缩、绝热、无黏、理想MHD 方程组[43]对磁场控制下RMI 过程进行求解,其守恒形式为

式中,ρ为气体密度,ν为速度,B为磁感应强度,μ为理想磁导率,μ=1.P*为流场总压,其由p和磁压力PB=B2/(2µ) 组成;E则为流场总能,其由内能Ep=p/(γ-1)、二维情况下的动能Ek=Ekx+Eky=ρ(ν·ν)/2和磁能量Eb=B·B/(2µ) 组成,它们表达式为

式中,p为气体压强,γ为气体比热比.

本文采用非分裂的6-solve CTU[41,42]算法求解理想的MHD 方程组.为了保证磁场散度为零(∇·B=0),结合了用于计算电场的CT[43]算法对Godunov 通量进行重构[43-46].同时采用Roe线性黎曼求解器以及带特征变量限制的三阶重构.程序采用结构化网格,其中,传统的流体力学变量,即密度、动量、能量放在网格中心位置,取周围变量的均值.磁场及其演化量放在网格面及边上.

计算流程如下: 首先计算包含MHD 源项的左右界面状态以及与之相关的界面流量;再通过CT 算法,使用上一步的数值流量和中心元胞涉及的电场计算CT 电场;随后,在每一个界面通过δt/2 步长的横向流量梯度生成PPM (piecewise parabolic method)界面状态(动量和能量密度的MHD 源项通过使用MHD 方程的原始变量形式来计算PPM 界面状态);对于每个通过δt/2 步长更新的界面状态,计算相关的流量;再通过CT 算法,使用上一步的数值流量和通过元胞平均状态(在时间步n+1/2)计算的元胞中心涉及的电场,计算CT 电场;最后,从n至n+1 步更新结果,流体动力学变量(质量,动量和能量密度)通过标准流量积分关系进行迭代,同时,磁场的界面平均组分通过斯托克斯循环积分进行迭代,更多细节可详见文献[47].

2.2 计算模型

本文的研究内容是文献[24]的延伸,因此数值模型除了磁场条件外均与之相同,包括与之相同的网格数(2500 × 890)及网格类型(笛卡尔网格)[24].图1 显示了1.22Ma的入射平面激波与半径R0=0.025 m高斯分布氦气柱相互作用的仿真模型.表1 给出了空气与氦气的气体参数.计算域为L(长) ×H(高)=0.25 m × 0.089 m,坐标系原点(0,0)位于气柱圆心上.上下边界设为壁面反射条件,右边界为流出条件.t=0 时刻,流场的压力与温度均匀分布,p0=101.325 kPa,T=300 K.图1(b)为氦气组分沿y=0 的分布,其分布方程如下[24]:

表1 气体参数Table 1.Gas parameters.

图1 (a)仿真模型;(b) 氦气组分沿y=0 的分布Fig.1.(a)Simulation model;(b) distribution of helium components along y=0.

式中,r为与气柱圆心的径向距离,a为气柱最外层氦气质量分数,此处取为a=0.8.

表2 为本文所考虑的3 种工况条件.为了探讨磁场方向对气柱不稳定性的影响,假设气体已电离成等离子态,分别选取了流向和纵向2 种方向磁场进行仿真研究,无量纲化的磁场强度β-1=(B2/μ)/(2p)=0.0005.

表2 不同工况下的初始条件Table 2.The initial conditions of different cases.

3 结果与讨论

图2 是本文仿真结果(下)与相同条件实验结果[9](上)的对比.可见,数值结果对激波结构及界面形态演变的描述均与实验结果[9]相符.

图2 平面激波与He 气柱作用过程的计算结果(下)与相关实验结果[9](上)的对比Fig.2.Comparison between simulation results of this study (down) and experimental results [9] (up) for the interaction between a shock wave and a helium gas cylinder.

图3 为不同时刻激波冲击He 气柱过程的阴影图.可见,激波与圆形气柱作用过程的流场结构呈轴对称(图2 和图3).入射激波(incident shock wave,IS)作用于气柱界面时(0.055 ms ＜t),出现非常规马赫反射,形成了向下游传播的透射激波(transmitted shock wave,T)、向上游传播的弧形反射稀疏波(reflected rarefaction wave,R)及反射激波(reflected shock wave 1,R1).由于氦气的声阻抗系数小于空气,T 的传播速度快于IS,导致其在界面上折射并诱导形成自由前体激波(freeprecursor shock wave,F).F 与IS 在界面附近的作用,与马赫杆(Mach stem,M)和反射激波R 汇聚于三波点(triple point,Tp)处,组成的马赫反射结构(图3(a)).

由于界面内侧的气流速度高于外侧,内外侧速度梯度形成了KHI.在KHI 作用下,界面变厚、失稳并逐渐卷起KH 涡串(图3(c),(d)).同时,左侧界面在激波冲击作用下向右侧压缩,逐渐过渡到近直线状(图3(b)).当T 穿过气柱右侧界面时,反射激波R1(图3(a))与计算域上下壁面发生碰撞并反射,形成向中心轴传播的壁面反射激波(wall reflected shock wave,WS).当壁反射激波在中心轴上汇聚后,在汇聚点周围形成局部高压区,迫使此处气流加速向右冲击气柱中轴界面,形成一道空气射流(jet,J),并诱导界面发生RTI.

由图3(c)可见,在空气射流高速冲击作用下,左侧气柱界面快速向右运动,气柱结构转变为“月牙形”[24].随着射流J 继续向右推进,射流头部卷起形成类似RTI 偶极子涡对的“蘑菇形”涡结构(图3(d),V1与).当射流头部到达右侧界面时,气柱中轴被压扁形成气桥结构(图3(e),bridge,B).此时,右侧界面开始被卷入涡V1()内,在偶极子涡对内形成了复杂的涡结构(图3(e)),并伴随有涡合并现象.因此,在V1()尺度不断增大的同时,逐渐融合成主涡结构(图3(f),(g),V(V')).随后,在V(V')的卷吸作用下,外部的空气和氦气团被卷入V(V')中,使气桥 B 宽度变短,见图3(f)—(h).

图4 和图5 分别为流向、纵向磁场作用下,激波冲击氦气柱过程的阴影图.与无磁控情况 (图3)相比,磁场对流场波系结构的形成与发展并无明显影响,但明显抑制了界面不稳定性,使界面厚度变薄并维持光滑.此外,磁场方向对界面后期形态的演化具有明显的影响,主要表现在影响射流结构的演化上.

图4 流向磁场作用下,激波冲击He 气柱过程的阴影图Fig.4.Shadowgraph images for interaction of a shock wave with a helium cylinder in the presence of a transverse magnetic field.

图5 纵向磁场作用下,激波冲击He 气柱过程的阴影图Fig.5.Shadowgraph images for the interaction of a shock wave with a helium cylinder in the presence of a longitudinal magnetic field.

在流向磁场控制下(图4),界面形态与无磁场情况较为类似,主要差异在于涡结构的演化上.此时,流向磁场抑制了射流头部“蘑菇形”涡结构的形成,并在射流头部前缘形成一道滑移层(slip line,S)(图4(d)).当射流头部与右侧界面碰撞后(图4(f)),滑移层受扰动影响而卷起形成涡结构V 和V'.随后,在涡结构V 和V'的卷吸作用下,气柱界面开始呈现波浪状(图4(h)),但射流头部与右侧界面碰撞形成的气桥结构始终保持光滑的弧形结构(图4(f)—(h)).

对于纵向磁场情况(图5),界面形态演化则相对较为复杂,主要表现在射流结构的演化上.在纵向磁场控制下,上、下侧界面仍会向右侧运动,但未向中心轴向翻转.在t=1.25 ms 时,射流尾部界面由圆弧转变为斜直线状,如图5(h)所示.射流头部向两侧翻转而卷起形成V 和V',但其涡量较小,对射流剪切层的卷吸作用较弱.另外,射流头部端面始终未能与右侧界面碰触,并在后期发展中,射流头部与右侧界面的宽度W不断增大.

图6 为流场涡量分布图.激波与界面作用过程中,压力梯度与密度梯度方向的不一致,形成斜压涡量.其中,上侧界面流体的逆时针旋转形成正涡量层,而下侧界面形成负涡量层.界面受KHI 影响而开始失稳,涡量层卷起系列KH 涡.随后,射流头部主涡与右侧界面涡量层的作用,使右侧界面上的涡量层在卷入主涡内的同时,与邻近小涡相互作用,使主涡整体结构十分复杂,如图6(d1)所示.

图6 不同磁场方向下的流场涡量分布Fig.6.Vorticity contour in RMI under various magnetic-field directions.

施加磁场后,界面光滑且无涡量附着在界面上,RMI 形成的涡量向气柱内外两侧传输,出现涡量分层(图6(a2)—(d3))[25].随后,WS 反向作用于气柱上并诱导形成新的涡层,最终呈现出多涡量层结构(图6(b2)—(d2),(b3)—(d3)).对于流向磁场情况,在磁场方向与界面切线平行处会形成涡量层交汇点(N点).初期,N点出现在界面的上下顶端,随着空气射流向下游的发展而被卷入到主涡内;当射流头部与右侧界面碰撞后,头部界面的正涡量层与右侧界面的负涡量层开始作用,卷起形成涡结构(图6(c2)).

在纵向磁场控制下,N点则一直位于左右侧界面的中心轴线上.当空气射流与左侧界面作用后,射流头部的上侧涡量为负,下侧则为正(图6(b3));随后,射流头部上侧负涡量层与左侧界面的正涡量层作用,涡量层失稳并卷起涡结构(图6(c3)).因此,磁场环境下,正负涡量层的耦合作用会诱导涡量层的失稳.

图7 为上半流场域的正环量、负环量及总环量随时间的变化曲线,环量为

图7 环量随时间的变化曲线(a)正环量Γ+;(b) 负环量Γ-;(c) 总环量ΓFig.7.History of the circulation: (a)Positive circulation Γ+;(b) negative circulation Γ-;(c) total circulation Γ.

由于流场的对称性,A只取上半计算域,即A=[-0.05,0.2] × [0,0.0445];ω 为涡量,Γ+和Γ-分别为正和负环量.环量可表示为流场中流体微团绕一个封闭曲线运动时所携带的旋转量,较大的环量意味着流体中存在较强的涡旋结构,因此涡量的存在可直接影响流场稳定性和能量分布.

由图7 可见,在激波与气柱作用前,无涡量产生,环量为0.在激波冲击界面瞬间,RMI 导致正负环量同时出现,且正环量大于负环量,总环量随时间呈近线性增大.当IS 扫过右侧界面后(t ＞0.12 ms),环量增速变缓.随后,WS 反向作用于气柱界面(0.18 ms ＜t),阻碍了涡量的形成,因此负环量增长速率减慢而正环量保持不变.当射流与右侧界面作用后(t ＞0.5 ms),伴随着偶极子涡对与气柱边界珠状小涡的合并,正和负环量又开始呈快速上升,并在气桥B 与主涡V 形成后趋于稳定,此过程中正负环量增速相当,总环量随时间呈现弱波动并逐渐趋于稳定.

虽然磁场的存在使涡量分层,但是初期(t＜0.18 ms)总环量与无磁场情况完全重合,这表明磁场的存在并未影响RMI 形成涡量,但磁场引起的涡量分层抑制了界面的KHI,使界面保持光滑.此外,当WS 与气柱相互作用后(t＞ 0.18 ms),纵向磁场控制下的总环量有所升高,且明显高于无磁场和流向磁场情况.这是由于纵向磁场控制下的涡量分层更为明显(图6),这增强了壁反射激波WS 与气柱界面的反向作用,减缓了负环量的增长速率,但正环量与流向磁场情况基本一致,导致总环量有所升高.另外,空气射流在磁场的作用下变慢,推迟了与右侧界面的碰触时间,使t＞ 0.5 ms 时的环量变化曲线逐渐呈现明显的差别.其中,无磁场下Γ+和Γ-的值更大,且在0.5 ms ＜t＜ 0.63 ms时,具有明显的增幅.由于磁场的作用影响了偶极子涡对的形成,使环量与无磁场情况的差距增大.其中,纵向磁场对射流的抑制最为明显,总环量降得最快,因此对抑制界面不稳定性的效果也最明显.

图8 展示了气柱演化过程中特征尺寸随时间的变化,其中气柱长度L为气柱界面左右端点的流向长度、气柱高度H为上下端点的纵向长度、中轴宽度W为气桥左右侧端点轴向长度.在IS 冲击气柱过程中(0.055 ms ＜t),L受到IS 的压缩而开始下降(图8(a)),在激波绕过右侧界面时(t=0.2 ms)达到最小值(L=0.034 m),而后空气射流J 与气柱的作用使L呈线性增大;期间,中轴宽度W则不断减小,在B 形成时接近于零(图8(c),t=0.62 ms).对于H,IS 绕过界面上/下端点时开始增加(图8(b),t=0.09 ms),而壁反射激波的反作用过程(0.09 ms ＜t＜ 0.75 ms)对气柱进行压缩,使H的增长速率减缓而呈抛物线增大.随着B的形成,加强了主涡结构对界面的卷吸作用,H开始下降(t ＞0.75 ms).

图8 激波冲击氦气柱过程中气柱特征尺寸的变化曲线(a)气柱长度L;(b)气柱高度H;(c)中轴宽度WFig.8.The variation curves of characteristic dimensions of helium cylinder during the interaction process with the shock wave:(a)Length of helium cylinder,L;(b) height of helium cylinder,H;(c) axis width of helium cylinder,W.

磁场对气柱尺寸的影响在界面失稳后(t＞0.2 ms)开始显现.其中,磁场降低了L的增长速率(图8(a)),使其与无磁场情况的偏差随时间增加.t＞ 0.55 ms 时,磁场对J 的抑制作用延迟了J 开始冲击右侧界面的时间,使L与无磁场情况的差距进一步增大.另外,纵向磁场下L的增长速率最小(0.2 ms ＜t＜ 0.8 ms),且因t＞ 0.8 ms 时无气桥B 形成,气柱右侧界面不会受到中轴涡量诱导而形成涡,因此纵向磁场对L的影响比流向磁场更强.

对于H(图8(b)),因纵向磁场控制下顶端界面并未能形成主涡,其值较小.t＞ 1.1 ms 后,由于涡V(V')的卷吸作用较弱(图6),H开始增大.至于W(图8(c)),由于磁场抑制了射流向下游的发展,B 的形成时间被延迟,因而降低了W的下降速率.其中,纵向磁场对J 的抑制作用更明显且无B 形成.因此W下降到一定值后开始升高(t＞ 0.61 ms).

为了说明磁场方向对激波驱动界面不稳定性的影响,图9 为气柱体积压缩率V(t)/V(0)随时间的变化曲线.其中气柱体积V(t)[48]为

图9 不同磁场方向下,He 气柱体积压缩率V(t)/V(0)变化Fig.9.Volume compressibility of He cylinder for various magnetic-field directions.

式中,VCi是单个网格的体积(1×10-8m2),χs(i,t)是t时刻该网格内He 所占的气体体积分数[48],V(t)/V(0)越小表示受压缩程度越高.在作用前,V(t)/V(0)=1[25].无磁场情况下,在IS 冲击界面瞬间,体积压缩率开始快速下降,并在t=0.24 ms后逐渐趋于稳定值(～0.73).由于后期大涡形成后,其卷吸作用很强,空气不断的被卷入而He 气团受其作用而逐渐减小(图6(d1)),一增一减体积压缩率在0.73 处变化.施加磁场后,V(t)/V(0)的发展趋势与无磁场情况相似,但大于无磁场情况,且流向磁场作用下的V(t)/V(0)大于纵向磁场情况.这是由于流向磁场作用下的主涡结构由S 诱导而成,其将空气自气柱下游卷入气团内,且不会与J 两侧界面相作用,因此气柱只有流入而没有流出(图6(d2)).而纵向磁场情况下的气柱形状稳定,无大涡形成,减小了卷入气柱内的空气量(图6(d3)),因而涡形成后气柱体积压缩率增长略高于纵向磁场.基于此,纵向磁场较流向磁场具有更强的形变控制效果及混合抑制效果.

为了探讨磁场方向对RMI 的影响机理及控制气柱尺度变化的原因,图10 给出上半计算域流场能量的变化曲线.随着IS 从左向右的传播过程,总能E,内能Ep及横向总动能Ekx呈线性增加,且有/无磁场条件时基本重合.当IS 传出流场后(t＞0.56 ms),无磁场及流向磁场情况下的E,Ep和Ekx快速趋于稳定值,但在纵向磁场控制下的则呈快速下降.由于Ep占E主要部分(图10(a)),图11 给出了流场内能分布.可见,磁场明显改变了流场内能分布规律,其中流向磁场提升了涡核中心的内能(图11(b2),(b3)),使低内能区由涡核中心转移到气柱界面上(图11(b3),(b4));对于纵向磁场情况(图11(c1)—(c4)),气柱界面上的内能更低且分布更广,因而阻碍了外围流体进入气柱内.此外,射流头部的高内能(高压)区(图11(c2),(c3))阻挡了空气射流向右侧界面推进.

图10 流场能量随时间的变化曲线(a)总能E、内能Ep ;(b) 磁能Eb;(c) 横向总动能Ekx;(d) 纵向总动能EkyFig.10.Time history of the flow energy: (a)Total energy E,internal energy Ep;(b) magnetic energy Eb;(c) lateral total kinetic energy Ekx;(d) vertical total kinetic energy Eky.

图11 内能分布随时间的变化Fig.11.Distribution of internal energy over time.

流场初始磁能量Eb为0.5 J (图10(b)).对于流向磁场,在IS 与气柱作用前 (0 ＜t＜ 0.055 ms),Eb维持不变;冲击气柱界面后(t＞ 0.055 ms),Eb开始上升,但上升速率逐渐减弱.对于纵向磁场,由于激波传播方向垂直于磁场方向,Eb在激波进入计算域时即开始上升且明显高于流向磁场情况.

图12 展示了磁场方向对流场磁能量分布的影响.可见,在IS 驱动阶段,流向磁场控制下的Eb主要沉积在涡量层间 (图12(c),(f)).其中,横向分量Bx(图12(a))主要在界面左侧涡层内增强,且在气柱中轴左侧被放大;纵向分量By(图12(b))出现在界面左侧涡层内,在气柱外侧也有By.对于纵向磁场情况,Eb不仅沉积在气柱界面涡层之间,还存在于气柱外激波扰动区内.其中,By(图12(d))主要存在于界面左侧涡层内,且在左侧界面涡层内下降并形成负向磁场,但在右侧界面涡层内放大.同时在激波扰动区域,By也存在磁场放大的现象,Bx(图12(e))主要存在于界面左侧涡层内.因此,纵向磁场控制下的磁场放大效果更为明显.

图12 t=0.15 ms 时,流场磁场强度与磁能量的分布云图Fig.12.The distribution of magnetic field strength and magnetic energy in the flow field at t=0.15 ms.

等离子体在磁场作用下会形成洛伦兹力,由安培环路定律和磁场无散度约束得

式(中,B)·∇B为磁张力项T;B2/2为磁压力PB,∇B2/2为磁压力梯度.研究表明[25],Eb的沉积区域会导致磁感线的扭曲和磁场的放大.磁感线的扭曲与磁场的放大在流场中形成磁张力和磁压力.

图13 为流场速度分布图,可见磁场明显改变了气柱内外侧速度分布.图14 为不同方向磁场控制下磁张力(图14(a1)—(d1))和磁压力(图14(a2)—(d2))分布矢量图.可见,无磁场情况下J 的发展速度最快,界面内外大速度梯度诱导偶极子涡对V1()快速形成(图13(a1)),而后迁移至界面外侧形成主涡(图13(b1)).施加磁场后,上半域界面涡量为正值(图6),而磁张力(图14(a1)—(d1))主要分布在涡量层上且矢量方向与速度剪切应力方向相反(图13),因而抑制了因速度剪切所引起的KHI.

图13 氦气柱附近速度矢量分布图Fig.13.Velocity vector distribution near the He cylinder.

图14 涡量层上的磁张力和磁压力分布Fig.14.Magnetic tension and magnetic pressure distributions on the vorticity layers.

磁场方向的不同改变了气柱界面受力及界面内外侧的速度梯度,使主涡结构的形成与发展出现差异(图13).对于流向磁场(图14(a1),(b1)),磁张力在涡层上的分布较少,其主要位于N点,这推动了N点自外向气柱内层运动;同时,磁张力对中轴气泡的作用(图14(b1))降低了J 的推进速度(图13(a2),(b2)),延迟了B 的形成,因而抑制了偶极子涡对的卷起(图14(b1));另外,滑移层S 以较高的速度向气团内部运动(图13(b2)),并在气团内部诱导形成主涡.对于纵向磁场,磁张力出现在气柱中轴内部,其方向与J 及气泡的运动方向相反(图14(d1)),明显抑制了气柱中轴速度(图13),使纵向磁场情况无B 形成.另外,射流运动被抑制(图13(a3))的同时,减慢了W的降低速率(图6(c)),降低了气柱内外侧速度梯度,最终仅在J 头部形成涡对(图13(b3)),因此纵向磁场对界面不稳定性的发展具有更明显的抑制作用.由此可知,磁张力通过改变速度场分布和抑制射流的运动,削弱了气柱受J 作用而拉长的同时,减缓了W的下降速度;而界面内外侧速度梯度的减小,减弱了涡的卷吸作用,影响了涡结构的形成,使不同工况下的H存在差异.

结合图6(a2),(b2)和图14(a2),(b2)可见,流向磁场情况下的双涡层在磁压力的作用下相互接触形成N点.纵向磁场情况下形成的磁压力推动涡层自界面向外迁移(图14(c2),(d2)),同时在“气泡”前端,磁压力斜向上推动射流尾端界面形成斜直线状(图14(d2)).可见,磁压力对界面涡量迁移及涡层迁移具有主要作用.

为了阐明磁场作用下速度剪切对界面不稳定性的影响,图15 为流场纵向动能Eky分布与涡量对比图,涡量图中的黑色轮廓为气柱界面边界.可见,无磁场情况下,Eky主要存在于气柱气团尾端与偶极子涡对中心(图15(a1)),Eky的存在反映出界面纵向运动的剪切应力.界面未受IS 扰动时(t＜ 0.055 ms),无纵向动能的分布(图10(d));在IS与界面作用后(0.055 ms ＜t＜ 0.18 ms),Eky主要分布于气柱受激波扰动部分的左侧,并在该时间段内快速升高.Eky在WS 接触气柱时,达到峰值(t=0.18 ms)(图10(c));在WS 与气柱相互作用过程中,Eky开始下降(0.18 ms ＜t＜ 0.29 ms)(图10(c)),Eky经过WS 在计算域上、下壁面间的多次反射(t＞ 0.29 ms),呈上下波动趋势,且Eky波动幅度随时间的推移而减小(图10(c))[24].

图15 纵向动能分布(upper)和涡量(lower)对照图Fig.15.Comparison of longitudinal kinetic energy distribution (upper) with vorticity (lower).

施加磁场后,IS 驱动段开始至WS 首次反作用时(0.055 ms ＜t＜ 0.29 ms),Eky与无磁场情况基本无差别(图10(d)).随着磁张力与磁压力对气柱形变及涡结构的作用,Eky逐渐低于无磁场情况,且WS 来回反射所导致Eky的周期性波动幅值也随之降低(图10(d)).另外,磁场对偶极子涡对形成的抑制作用,使流场Eky低于无磁场情况(图15(b2),(b3)).其中,纵向磁场作用下因无偶极子涡对的形成(图15(c1)—(c3)),使流场Eky较流向磁场情况要更低且更平缓.而施加流向磁场后,受磁张力与磁压力的作用,Eky主要分布在界面外的涡量层上,而界面位置上基本无速度梯度,因此气柱界面保持光滑.

4 结 论

本文采用CTU+CT 算法求解理想MHD 方程组,仿真研究了不同方向磁场控制下激波与高斯分布轻质气柱界面的相互作用.对于无磁场情况,激波与轻质界面作用过程中波系结构、涡以及界面不稳定性的演化与相同初始条件的实验相吻合.磁场的存在有效抑制了界面不稳定性及涡结构的形成与发展,使界面保持光滑.另外,本文通过对比分析不同方向磁场控制下的流场结构、环量、磁能量、磁张力、磁压力等变化,进一步揭示了磁场对界面不稳定性的抑制机理,得出以下结论.

1) 对于无磁场情况,激波冲击氦气柱界面过程诱导了RMI.在RMI 作用下,气柱内外侧存在明显的速度梯度,KHI 诱导界面失稳而形成珠状小涡结构.由于气柱内侧为轻质氦气,射流冲击气柱中轴过程,界面在RTI,KHI 和RTI 耦合作用下失稳卷起主涡.施加磁场后,界面出现的涡量分层现象明显抑制了KHI,但后期正负涡量层的耦合作用会诱导界面失稳.此外,RTI 的发展亦受到磁场抑制,其中纵向磁场对中轴射流的影响更为明显,使中轴射流受磁场作用而发展受阻.

2) 在不同方向磁场控制下,激波与氦气柱作用产生的涡结构出现明显的差别,改变了气柱尺寸特征.由于纵向磁场抑制了中轴射流的发展,减缓了中轴宽度的降低速率,对气柱长和高度的影响相较于流向磁场情况要更明显.此时,流向磁场作用下的气柱体积压缩率略高于纵向磁场,因此纵向磁场控制下的界面更为稳定,能更好地抑制气体的混合.

3) 磁场线受界面不稳定性的影响而发生扭曲,进而形成磁张力和磁压力.其中,磁张力主要集中在涡量层上,其作用方向与涡量卷起方向相反,有效抑制了KHI;且纵向磁场控制下磁张力对中轴射流的抑制更为明显,使其对RTI 的抑制效果优于流向磁场情况.另外,磁压力作用方向是自界面指向涡层,因此推动了涡量的迁移进而诱导出涡量分层现象,抑制了RMI 对界面的影响.

5 发展与展望

未来,对磁场控制RM 不稳定性的数值研究可朝以下几个方面开展: 1) 提高数学物理模型的复杂性和准确性,引入可模拟具有复杂边界的壁面边界处理方法,如使用离散玻尔兹法或双流体模型,以提升对其内在流动现象进行更一步揭示的能力;2) 探索多尺度和多物理场的耦合,建立更全面的模型;3) 优化MHD 控制策略,通过系统性调整控制参数找到最有效的方案.此外,由于高温高压的实验条件太过严苛,以目前的科技能力(材料学等)还无法完成该类型的中大尺度的实验研究,因此未来还需通过开展实验验证来验证仿真模型的可靠性,而多领域合作则有助于整合各个领域的专业知识,推动磁场控制RM 不稳定性仿真研究的深入发展.最后,通过探索新型材料与技术,并将研究成果引入可持续能源领域,本研究有望为未来能源的发展提供新的解决方案.