李尚卿 王伟民 李玉同4)‡

1) (中国科学院物理研究所,北京凝聚态物理国家研究中心,北京 100190)

2) (中国科学院大学物理学院,北京 100049)

3) (中国人民大学物理学系,北京 100872)

4) (松山湖材料实验室,东莞 523808)

1 引言

激光等离子体与外加磁场相互作用是高能量密度物理的一个重要研究课题.它在磁化惯性聚变[1,2]、实验室天体物理[3-5]等领域均有应用.人们已经实施了不同构型的外加磁场与激光等离子体相互作用的实验[6-9],在这些实验中,采用的激光脉宽一般在1 ns 量级,激光强度约在1014W/cm2量级,产生等离子体的空间尺度约在1 cm 量级,时间尺度约在10 ns 量级.在这些实验参数条件下,一般用磁流体力学来描述等离子体的演化.目前可以模拟激光等离子体的磁流体力学程序有FLASH[10],GORGON[11],PERSEUS[12]等.随着研究的深入,设计了更加复杂的磁场-激光等离子体相互作用的实验构型.实验设计复杂度的提高主要体现在两个方面,一是打靶方式变复杂,如Ryutov[13]提出的环形打靶方案;二是磁场构型变复杂,如磁化惯性聚变中的会切磁场构型[14](cusp magnetic configuration).面对这些复杂的设计实验,需要开发新的模拟程序.本文发展了一款基于流体数值模拟程序OpenFOAM[15]的磁流体求解器.

OpenFOAM 是一个由C++程序编写的面向对象的开源计算流体力学(computational fluid dynamics,CFD)程序平台,全称为Open Source Field Operation and Manipulation.它拥有数十种求解器,可以实现多相流、热、电磁、化学反应等各种流动的数值模拟.由于OpenFOAM 是开源的,用户可以根据需求添加物理模型,这为开发磁流体力学求解器提供了便利.不少研究团队利用OpenFOAM发展了新的求解器来进行磁流体和等离子体模拟的研究,Singh 等[16]开发的不可压缩磁流体求解器可以模拟液态金属在管道内的流动;Xisto 等[17]用双层隐式压力分离(pressure-implicit with splitting of operators,PISO) 算法,将OpenFOAM 自带的不可压缩磁流体求解器拓展到可压缩流;Ryakhovskiy 等[18]开发了低磁雷诺数条件下的超声速磁流体求解器;Charles 等[19]开发了Kurganov-Noelle-Petrova格式[20](KNP 格式)和Kurganov-Tadmor格式[21](KT 格式)的理想磁流体黎曼求解器.此外,有研究团队将OpenFOAM 与传统的动力学模拟结合起来开展工作,譬如OpenFOAM 与粒子模拟(particle-in-cell,PIC)的结合进行模拟[22].

本文基于OpenFOAM 平台开发了一种新型磁流体力学求解器MHDFoam,用于求解二维或三维的可压缩跨音速束流.该求解器将磁场PISO算法植入OpenFOAM 自带的KT 格式的黎曼求解器rhoCentralFoam,可以控制磁场散度误差,保证模拟数值精度,避免了非物理现象的出现.MHDFoam 的收敛阶在1—2 之间,标准算例的测试结果与FLASH 等程序符合较好.然后,利用开发的求解器计算了激光等离子体分别与外加均匀轴向磁场和电容线圈产生的非均匀磁场相互作用.初步模拟结果表明,在外加均匀轴向磁场条件下,激光等离子体喷流会出现喷嘴和结节,喷嘴位置和结节之间的长度与热压比开方成线性关系;在电容线圈产生的非均匀磁场条件下,结节呈非线性分布,而喷嘴的位置受线圈电流参数调控.当线圈中心磁场相同时,小尺寸线圈产生的磁场会加快喷嘴和结节的形成,等效的均匀轴向磁场更大.此模拟结果表明,求解器MHDFoam 的特点是面向工程和实验,可以做复杂构型下的磁流体模拟.

2 程序开发

2.1 基本方程

从磁流体守恒方程出发,Xisto 等[17]推导出了以流体密度ρ、流体速度U、流体压强p、磁场B为基本变量的磁流体力学控制方程组:

其中,η为电阻率,ρe为总能量密度.对于采用理想气体模型的磁流体而言:

方程组(1)即为不考虑黏性、热传导等效应的阻抗磁流体方程组,MHDFoam 求解时可以添加这些非线性项.其他磁流体参数如电流密度J,电场E等均可以用基本变量组 (ρ,U,p,B) 结合麦克斯韦方程组、欧姆定律等导出.

2.2 数值算法

基于OpenFOAM 的磁流体力学求解器MHDFoam 的更新算法如图1 所示.它由两部分组成:第一部分是OpenFOAM 自带的中心格式密度基黎曼求解器rhoCentralFoam,可求解可压缩跨音速气体;第二部分是一个专用的PISO 算法,用来控制磁场散度误差.PISO 算法是计算流体力学中的经典算法,其可以求解不可压缩流体.PISO 算法步骤是先用上一时刻的压强计算出预测速度,再通过构建压力泊松方程更新压强,然后对预测速度进行修正.依此循环若干次,就可以保证场量组(U,p)更新时也满足∇·U ≈0 .Weller 等[15]将PISO算法迁移到不可压缩磁流体求解器中,可以让磁场B更新时满足∇·B ≈0 .Xisto 等[17]将PISO 算法拓展到了可压缩磁流体束流.参考了以上磁场求解方案,开发的MHDFoam 求解器运行过程如下.

图1 MHDFoam 求解器的更新算法示意图Fig.1.Chart flow of update algorithm in the MHDFoam solver.

步骤一求解器首先对密度ρ,动量ρU等流体场量进行KT 格式的高分辨中心差分[21],然后依次求解控制方程组(1)中的密度、动量和能量方程,得到更新后的密度ρ、速度U、压强p和总能量密度ρe.此阶段磁场B作为常量参与计算,没有更新.

步骤二求解器进入BPISO 循环.首先求解方程组(1)中的磁场演化方程,获得预测磁场B*.该磁场不满足无散条件,还需修正.根据电动力学理论,磁场可以写成B=∇×A+∇ϕ的形式,其中A为矢势,ϕ为标势.则可以构造一个虚拟“磁压”泊松方程∇2ϕ=∇·B*,解出标势ϕ.磁场需修正为B=B*-∇ϕ,再将修正后的磁场代入方程组(1)中的磁场演化方程,得到新的预测磁场B**.依此循环若干次,最终获得的更新磁场B可以满足∇·B ≈0.该算法本质是Brackbill[23]投影算法.

2.3 程序检测

为了验证本文开发的磁流体力学求解器的可靠性,采用了一些标准算例来对数值模拟结果进行了检验,检测结果表明本文发展的程序与这些标准算例结果符合得比较好.本文采用的第一个算例为奥萨格-唐磁流体涡旋(Orszag-Tang MHD vortex)问题[24].其速度和磁场的初始条件为

其中,x,y∈[0 1]2,上下两边和左右两边采用循环边界条件.γ=5/3 为绝热指数.初始速度和磁场的分布如图2 所示.初始压强和密度为均匀分布,p0=1/γ,ρ0=1.磁导率µ=1 .以上参数均为无量纲量.图3(b),(d)分别为t=0.5 时刻的密度和磁场分布.可以看到随着时间演化,模拟区域出现了许多波交叠的结构,在各个空间尺度下均呈现湍流的特征[25].图3(a),(c)分别为t=0.5 时y=0.25 处的密度和磁场分布.可以看出本文开发的求解器MHDFoam 与FLASH 程序8wave 求解器[25]和NIRVANA[26]程序的运算结果一致.

图2 奥萨格-唐磁流体涡旋的(a)初始速度场和(b)初始磁场Fig.2.Initialization of speed field (a) and magnetic field (b)in Orszag-Tang MHD vortex.

图3 t=0.5 时奥萨格-唐磁流体涡旋的模拟结果 (a),(c) y=0.25 处密度和磁场比较;(b),(d) 密度和磁场廓线Fig.3.Simulation results of Orszag-Tang MHD vortex at t=0.5:(a),(c) 1D cut comparisons of density and B-field at y=0.25;(b),(d) density and B-field contours.

采用相对误差讨论了MHDFoam 求解器的收敛阶数.对于t时刻N×N网格奥萨格-唐算例的模拟结果,变量W的相对误差定义如下[19,27]:

收敛阶数的计算公式:

其中为N0粗糙的网格数(N0＜N).模拟结束时的相对误差和收敛阶数分布如表1 所示.可以看到MHDFoam 求解器和KT-MHD 程序的求解器的收敛阶数均在1—2 之间.而无磁场情况下的KT格式黎曼求解器的收敛阶数为2,这说明收敛阶数的降低主要是磁场算法导致的.由于MHDFoam采用了磁场解耦算法,收敛效率比非磁场解耦的KT-MHD 程序的还要稍微低一些.采用这种解耦算法的目的是既能保留中心差分的黎曼求解器,又可以灵活地调整磁场的演化方程(比如考虑双流体效应时需添加的毕尔曼电势、霍尔电势等),这在传统磁流体黎曼求解器中几乎无法做到[12].需要注意的是表1 中两个程序的相对误差是分别和自身的高分辨网格的模拟结果相比,并不能直接反映和精确结果的误差.关于该误差的讨论详见文本的补充材料S1 部分,证明了MHDFoam 求解器的时间稳定性和空间收敛性.

表1 奥萨格-唐问题的相对误差和收敛阶数Table 1.Relative errors (δN) and convergence order (RN) for Orszag-Tang problem.

磁流体转子(MHD Rotor)[28]是另一个经典二维算例.其模拟空间为x,y∈[-0.5 0.5]2,四边均为零梯度边界条件(zero-gradient boundary condition)[26].初始条件为

图4 磁流体转子 (a) 初始密度廓线;(b),(c) t=0.15 时密度和磁场廓线;(d) t=0.15 时x=0 处磁场比较Fig.4.MHD rotor:(a) Initial density contours;(b),(c) density and magnetic field contours at t=0.15;(c) 1D cut comparisons of B-field at x=0 at t=0.15.

图5 是磁流体转子的磁场散度误差分布.对于600 × 600 网格的模拟结果,磁场的散度误差约在10—9量级.FLASH 程序采用了6 阶自适应网格细分技术,磁场误差可以达到10—12量级[25].由于在本算例中MHDFoam 并未采用自适应网格细分,本算例中磁场误差还达不到FLASH 的水平.对于KT-MHD 程序,磁场误差在间断面附近达到了约1 量级[19](400×400 网格,奥萨格-唐问题),该现象在图5 中并未出现.综合与参考程序FLASH 和KTMHD 的比较,可认为MHDFoam 的磁场误差控制的精度较高,可以避免非物理数值现象的出现.

图5 t=0.15 时MHDFoam 求解磁流体转子问题的磁场散度误差Fig.5.Divergence of magnetic fields using the MHDFoam solver at t=0.15 for the MHD rotor problem.

上述两个算例表明,开发的新型求解器MHDFoam 的测试结果与FLASH 和NIRVANA 等程序的一致.其他标准算例的检测结果和兼容自适应网格的算例详见本文的补充材料.

3 应 用

3.1 模拟配置

分别讨论均匀轴向磁场和电容线圈产生的非均匀磁场与激光等离子体喷流相互作用的情形.根据文献[29],当一束能量约为100 J、脉宽为1 ns、焦斑直径约为200 µm 的激光照射到CH 靶上,产生的激光等离子体的典型速度约为500 km/s,温度约为106K,压强约为108—109Pa,密度约为10—4g/cm3.此时辐射致冷效应可以忽略,可以用理想磁流体力学来模拟等离子体的演化.对于轴对称喷流,二维模拟足够描述喷流的基本动力学过程[29].因此本文选用计算量比较小的二维模拟作为算例进行分析.

如图6 所示,模拟空间为一个3 mm × 12 mm的矩形区域.它被划分为240×960 的网格,网格大小为12.5 µm.以矩形底边中心为原点,横向为x轴,纵向为z轴建立坐标系,则有—1.5 mm≤x≤1.5 mm,0≤z≤12 mm.为了更关注激光等离子体的磁流体力学行为,将激光等离子体的产生简化为一种流体注入的情况.设置入射口处(红线标注)有密度为1.0 × 10—4g/cm3、温度为1.6 × 106K、压强为1.0 ×108Pa、速度为480 km/s 的高温高密等离子体流入.绿色线处为墙边界条件,蓝线处为无反射边界条件.模拟空间内存在密度为1.0 ×10—8g/cm3,温度为300 K 的低密度背景等离子体,它和入射等离子体的质量比为1∶104.磁化背景等离子体的磁场在10—100 T之间变化,注入等离子体在进入磁场之前为非磁化的,入射口处磁场为0.

图6 二维模拟配置Fig.6.Setup for 2D simulation.

3.2 均匀轴向磁场情形

根据文献[29],当外加轴向磁场时激光等离子体喷流会出现密度、压强等参数的周期性增强结构,其中此结构的极大值处称为结节.图7 展示了MHDFoam 求解器在图6 初始配置下的模拟结果.可以看到,当t=22 ns 时,随着外加磁场的增强,结节的周期长度变短,在模拟区域内出现的次数增加.Lei 等[29]用连续马赫激波反射原理解释了结节形成的原因,并指出结节归一化长度其中L为结节长度/,D为喷流直径,βnozzle为喷嘴处的热压比为流体压强与磁场压强比值).本文模拟中也发现了类似现象.图8(a)中蓝色线、红色线、绿色线分别为图7 中磁场为40,60,80 T 时的轴线上的密度分布.本文定义原点处热压比为βo,喷流上第1 个密度极大值的位置(喷嘴位置)为S,第一个和第二个密度极大(第一个结节位置)之间的距离为L,喷流入口宽度D=1 mm.根据图8(b),(c),发现它们之间有如下近似关系:

图7 t=22 ns 时均匀磁场下的密度廓线Fig.7.Density contours at uniform magnetic fields at t=22 ns.

图 8 (a) t=22 ns 时轴线密度分布;(b),(c)参数L 和S 与关系Fig.8.(a) Density distributions at axis at t=22 ns;(b),(c) parameter L and S as a function of

该结论与文献[29]中用FLASH 计算的结果类似,这表明求解器MHDFoam 可以模拟外加均匀轴向磁场中的激光等离子体喷流.

3.3 电容线圈产生的磁场情形

2013 年,实验演示了纳秒激光作用电容线圈靶产生了寿命在纳秒量级、强度在千特量级的磁场[30],这为激光等离子体实验中的外加磁场拓展了新的参数范围.其中此种强磁场由于线圈中形成的电流产生.线圈电流磁场是非均匀的,除线圈平面法向的轴向磁场外,还具有横向分量,是一种典型的极向磁场.线圈磁场的空间分布依赖于线圈电流I和线圈半径a,由毕奥-萨伐尔定律可以推导出环形电流磁场的解析分布.为了研究线圈电流参数对等离子体喷流的影响,设计了3 种线圈磁场构型.构型(1):I=0.5 MA,a=3 mm;构型(2):I=0.25 MA,a=3 mm;构型(3):I=0.15 MA,a=1.8 mm.其中构型(2)的磁场分布和磁力线如图9(a)所示.3 种构型的磁场在x=0 处和z=0 处的分布如图9(b),(c)所示,其中橙色线代表构型(1),红色线代表构型(2),黑色线代表构型(3).构型(2)与构型(3)的线圈中心磁场相等(Bo≈52.4 T),但两者的轴线磁场梯度不同,构型(3)的磁场梯度较大.构型(1)的中心磁场(Bo≈104.7 T)为其他两个构型的2 倍,磁场大小处处为构型(2)的2 倍.

图9 线圈电流磁场 (a) xz 平面二维分布;(b) x=0 处分布;(c) z=0 处分布Fig.9.Magnetic field of coil currents:(a) 2D distributions in the xz-plane;(b) 1D cut at x=0;(c) 1D cut at z=0.

本文模拟了线圈电流产生非均匀磁场中的激光等离子体喷流,并与均匀磁场的结果进行了比对.当t=22 ns 时密度廓线的模拟结果如图10 所示,轴线处的密度分布如图11(a)所示.可见非均匀外磁场中结节位置分布也是非均匀的.对于z≤5 mm 区域,构型(1)、构型(3)和均匀外加磁场95 T 的喷嘴和结节位置近似重合,这表明在该区域内各构型磁场对喷流的作用效果近似相等.此外,构型(2)的线圈中心磁场Bo与构型(3)的相同,它的喷嘴位置和第一个结节位置却出现的较晚,作用效果相同的磁场大小约为53 T.

图11 t=22 ns 时模拟结果 (a)轴线密度分布;(b)无量纲参数S/D 在Ia 相平面上的分布Fig.11.Simulation results at t=22 ns:(a) Density distributions at axis;(b) dimensionless parameter S/D in the Iaplane.

对于非均匀线圈磁场,定义喷嘴位置值相同的均匀磁场为等效作用磁场Be,轴线上喷流密度近似相等的长度为等效作用长度为λe.则图10 的模拟结果的等效参数可总结到表2 中.

表2 图10 模拟结果的等效参数Table 2.Equivalent parameters of simulation results in Fig.10.

图10 t=22 ns 时图9 配置的非均匀磁场条件下的密度廓线及与均匀磁场的比对Fig.10.Density contours at nonuniform magnetic fields set in Fig.9 at t=22 ns compared with uniform magnetic fields.

考虑到轴线处构型(1)的初始磁场强度至少为构型(3)的2 倍以上,则小于后者的1/2,但在z＜ 5 mm 范围内两个磁场对喷流的作用效果近似相等(如图11(a)中的蓝色线和绿色线),这表明(6)式已经不再适用.可认为磁场的不均匀性起到了补偿作用.本文在非均匀磁场条件下,拟合出了(6)式的一个较为粗糙的修正关系:

初步的模拟结果表明,电容线圈产生的非均匀磁场在一定空间范围内对喷流的作用效果可能与均匀磁场的相同.在线圈中心磁场Bo相同时,线圈半径a越小,磁场不均匀性越大,等效磁场也会变大,喷流和结节出现的越早.对于约100 J 纳秒激光等离子体与线圈磁场相互作用的实验,可以设计参数a/D约为2,Bo约为50 T 的线圈磁场,它在约0.5 cm 空间范围内对激光等离子体喷流的作用效果与均匀约100 T 轴向磁场的大致相同.这是对文献[29]中均匀磁场-喷流结节结论的拓展,该初步结论还需进一步的探讨和论证.本应用算例同时还证明了MHDFoam 模拟复杂磁场构型下的磁流体问题的能力.

4 结论

为了实现激光等离子体的磁流体数值模拟,我们开发了一套基于OpenFOAM 的磁流体力学求解器MHDFoam.该求解器将一个专用的PISO 算法植入OpenFOAM 自带的求解器rhoCentralFoam中,用来保证磁场演化的同时较严格地满足磁场无散∇·B=0 的约束条件.MHDFoam 的收敛阶数在1—2 之间,奥萨格-唐涡旋和磁流体转子标准算例的运算结果与FLASH 等程序的一致.通过算法测试后,我们用求解器MHDFoam 模拟了外加均匀磁场对激光等离子体喷流的影响,给出了结节周期长度、喷嘴位置与入口热压比开方的线性关系.此外,本文还用MHDFoam 模拟了电容线圈产生的非均匀磁场中的激光等离子体喷流,讨论了磁场不均匀性对结节和喷嘴位置的影响.初步的模拟结果表明,当线圈中心磁场相同时,小尺寸线圈产生的磁场会加快喷嘴和结节的形成;在厘米空间范围内对激光等离子体喷流的作用效果与一个较强均匀轴向磁场的大致相同,从而可以替代实验中一些难以达到的强均匀磁场条件.以上两个应用算例表明,MHDFoam 可以求解在激光等离子体参数范围内的、不同磁场构型下的磁流体问题.下一步我们将在这个求解器中添加双流体和辐射模块,对激光等离子体实验进行更精细的模拟.