平面和圆柱面构形的磁流体力学计算

赵继波，孙承纬，谷卓伟，王桂吉，蔡进涛

(中国工程物理研究院流体物理研究所冲击波物理与爆轰物理重点实验室，四川 绵阳 621999)

平面和圆柱面构形的磁流体力学计算

赵继波，孙承纬，谷卓伟，王桂吉，蔡进涛

(中国工程物理研究院流体物理研究所冲击波物理与爆轰物理重点实验室，四川 绵阳 621999)

在一维弹塑性反应流体力学Lagrange 编码SSS的基础上，研制了多介质、多空腔、与集中参数外电路耦合的实验构形弹塑性反应磁流体力学一维计算编码SSS-MHD。在编码中，通过对构形边界与外电路的耦合以及不同构形磁场边界条件的统一处理等措施，解决了其中的关键问题。本文中编码的优势体现在两个方面：能够统一处理不同的磁流体力学实验构形，包括平面及圆柱面几何构形；可以同时进行磁流体力学和含能材料反应流动的计算。对平面准等熵压缩炸药和炸药爆轰驱动套筒压缩磁通量等两类典型实验进行了数值计算，结果与实验结果一致。

流体力学；计算编码；多场耦合；准等熵压缩

磁驱动实验构形的数值模拟已进行了较多研究，主要用于实验结果预估和结构参数优化设计等工作。对磁驱动准等熵压缩实验进行较成熟数值模拟的二维磁流体力学(MHD)计算程序是ALEGRA编码[1]，这是一种采用有限元和任意拉格朗日-欧拉方法(ALE)的MHD代码。其主要特点为：耦合了加载装置精细的电路模型，以便不依靠实验测量结果，在计算样品力学运动的同时准确地计算该实验的实际负载电流。R.W.Lemke等[2]、M.D.Knudson等[3-4]就利用该程序开展了很多理论计算工作，在他们建立的磁驱动准等熵加载一维模型中，样品表面速度的计算值和测量值的相对偏差已在1%以内，这表明由模拟计算得到的样品动力学行为及其状态是符合实际物理过程的。RAVEN[5]适用于一维Lagrange辐射磁流体力学计算的MHD编码。其磁流体力学方程组的能量方程，根据电子离子之间是否达到热平衡，又可分为单温方程和双温方程，如问题涉及到高温非平衡辐射时，则使用三温方程。该程序考虑了磁扩散模型和变化的电导率模型，耦合了外电路，具备对包括磁通量聚积发生器(MC-1)等复杂构形的模拟能力[6]。在对MC-1装置进行实验的磁流体力学数值计算方面有了比较多的工作，其中SMOG-DISK程序[7]是典型的可用于MC-1发生器的一维模拟程序，使用解析形式的状态方程和电导率函数，引入的混合物的状态方程和电导率模型能处理由特殊混合物制成的套筒。除此之外，还有MACH2和MACH3编码[8-10]以及LS-DYNA软件中建立的EM模块[11]能进行MHD计算。以上程序或软件有力地配合了磁驱动实验的进展，但由于此类编码保密甚严，一直未见开源编码。

在我国，研究人员在各自的工作基础上针对特定研究目标发展了不同的磁流体力学编码[12-14]，如MHD-ICE和MC11D等。这些编码有各自的特点和侧重点，能得到与实验符合得较好的计算结果，但是目前仅能对各自构形的模型进行研究，还不能进行多构形建模，同时也不具备外电路耦合功能。本文中，在SSS编码[15]的基础上，拓展了磁流体力学计算功能，不仅可解决磁腔与力学构形的一体化计算、与外电路的耦合计算等关键问题，还具备统一处理不同的磁流体力学实验构形、以及同时计算磁流体力学和反应流体动力学的能力。

1 SSS-MHD的基本框架

磁流体力学考察磁场对流体运动的影响，反之流体运动又影响电磁场，这就要求必须将控制导电流体介质运动的流体力学方程与控制电磁学现象的方程耦合求解。磁流体力学方程分为两部分：一组是控制力学变量的流体力学方程；另外一组是控制电磁学变量和电路参数的方程。

1.1 动力学方程

连续性方程、动量方程和能量方程分别为：

(1)

(2)

(3)

在动力学方程中，已包含了电磁学量的影响项，如动量方程中含有洛仑兹力项，能量方程中含有焦耳热项，表示磁场对流场的影响。另外，物质的本构和物态方程等表述内在变量、热力学量的关系在原SSS编码中已存在，在此不再单独列出。

1.2 磁场和电路方程

安培定律、磁扩散方程和电路方程分别为：

(4)

(5)

(6)

(7)

式中：Q、C0、L0、LL、I、R、Vs、Vg和U0分别为电容器剩余电荷量、储能电容量、传输线路电感量、负载(动态)电感、总电流(即负载电流)、传输线路电阻量、开关电压降、负载电压降和初始充电电压。

SSS-MHD编码中，主要从磁通量变化的角度来探讨构形电感对回路电流的影响，即每一时刻构形的动态电感，是由其总磁通量Φ对负载电流的导数计算的：

(8)

每一时刻磁流体力学负载构形的反电动势和电感量由磁扩散方程给出，每一时刻电路方程给出的电流值将被提供给磁扩散方程计算磁场边界值。这里需要特别说明的是，在一些磁流体力学构形计算中没有显示的电流值，即没有外电路的负载电流参与上述计算，此时可通过无源电路的涡电流来进行耦合。

1.3 爆轰反应流动计算

SSS-MHD能进行磁流体力学和反应流体动力学的一体化计算。反应流体动力学的一大特色就是反应率方程，即不基于详细的化学反应动力学机理，依靠变量宏观唯象表述。如果固态质量分数W代表反应进程的内变量，则W=0表示反应完成，处理反应流动的方法就是确定W的唯象方法。在爆轰数值模拟中，常用Forest Fire律来描述非均质炸药的冲击起爆过程，它以压力作为流体动力学热点的统计表征，表达式如下：

(9)

式中：C0、…、Cn称为Forest Fire系数，通常为14个。当ppCJ，W=0。

爆轰产物HOM物态方程是以爆轰产物BKW方程确定的CJ等熵线为参考线的变β物态方程，计算爆轰产物沿CJ等熵线膨胀释放的相对能量与圆筒实验结果相当一致。共有17个参数，除了定容比热cV和能量基准调整常数e0，其余均为拟合参数。

2 关键问题处理

2.1 磁流体力学构形描述

SSS-MHD编码可以进行平面和柱面一维的磁流体力学计算，相应的坐标系和坐标方向如图1所示。按电流的走向，可以将样品的MHD构形计算划分为θ-MHD计算、z-MHD计算，以及两者混合的Screw-MHD计算。

图1 MHD计算的构形和坐标系Fig.1 MHD configurations and coordinates

2.2 构形边界与电路耦合方式

在磁流体力学计算中，与外电路的耦合是解决问题的关键。然而，耦合的方式多种多样，对应的磁场分布也不相同，需要建立统一的描述方法。大致说来，如果磁场和电流的相对位置已知，界面或边界处的磁场可按总电流作理论计算。处理外电路的耦合主要应事先确定MHD构形中的有源区，可以是一个(只考虑主构形，不计回流导体)，也可以是两个(同时考虑主构形和回流导体)，分别称为单边和双边耦合。其中，单边耦合又分为导体区左边耦合电路和导体区右边耦合电路等2种情形。

边界与电路耦合的几何说明，如图2所示。图中，iMHD定义为计算分区的电磁物性指数，其值为0、1、2、3时，分别表示该计算分区不计电磁作用、存在磁场的空腔(或电介质)，理想导体和非理想导体等4种情况。MHD片区是指具有相同iMHD指数的彼此邻接的几个连片分区(包括空腔)。为了便于计算磁通量和程序的编制和运行，标示了左、右导体区，以及左、右电路点。

图2 边界与电路的耦合方式Fig.2 Modes of coupling with circuit

可以看出，电路区自身由相连接的导体区构成，其前后必须是非导体区。在双边耦合情形，磁场全部被封闭在两个电路区所围的范围内，他们的外面无论任何物质或空腔都不能存在磁场；在单边耦合情形，在与导体片区耦合边的邻区存在磁场或为磁腔，非耦合边的各分区中不进行磁流体力学计算。

2.3 磁场边界条件及其处理

根据安培定律，一个MHD构形载流导体周围的磁场只与该构形内部通过的总电流有关，与其内部电流分布细节无关，据此可得导体外面磁场的解。然而，由上述的磁场边界条件推导却表明，在任何导体/电介质(或空腔)的界面两侧：磁场切向间断(导体表面不存在线电流时)、法向连续，电场切向连续、法向间断(介质表面不存在面电荷时)。一维MHD计算无论是z构形或θ构形，界面处的磁场总是切向的(轴向或环向)，物质粒子的运动总是径向的，电场总是在电流的方向。除了等离子体双流体模型的场合，一般不考虑电场的计算。磁扩散计算中，导体表面只能有面电流密度，磁场切向分量也应是连续的，适合于磁场空间导数的差分运算。

图3表明，平面构形在一对电极板之外的区域、圆柱形构形在回流罩(蓝色柱筒)之外的区域，磁场都是零，磁扩散要算的区域只在图中红、蓝部件自身以及他们之间的区域。

图3 平面、圆柱面MHD构形的磁场围线Fig.3 Magnetic surrounding lines on planar and cylindrical MHD configurations

下面，用安培定律给出z构形和θ构形情形下界面处磁场的值。图3标出了有关导体周边的围线，一半处于成对导体之间的空腔(或电介质，还可以有感应导体，称为磁腔)中，另一半处于导体外侧磁场为零的空腔(或电介质)中，所包围的导体截面流过的电流是外电路提供的总负载电流。除了圆柱形z构形，其余情形中围线都可视为矩形的，主边长度等于导体的“高度”l，侧边与磁感应强度矢量垂直，因此线积分的值就是磁腔中主边所在处的磁感应强度B乘以该边长度l。从安培定律得出，这个值等于-μ0I，I是总电流Iz或Iθ，也就是：

(10)

上式表明，成对导体构成的磁腔中磁场是均匀的，只取决于负载电流的线密度i。容易看出，平面情形的z构形和θ构形实际上一样，只不过位相差90°而已，所以上式中把磁场和电流都加了下标说明。

圆柱形z构形中的围线是与构形同轴心的圆周线，磁场显然是轴对称的，线积分值等于圆周线长乘以当地的磁场值，安培定律给出：

(11)

编码中，用变量Cθ代替Bθ， 则Cθ在磁腔中为常数：

(12)

外空间磁场为零以及内磁腔磁场为常数，都不是导体自身中的磁场值。理想导体中磁场也为零，但流过的电流产生外磁场。非理想导体中的磁场和电流，要依靠磁扩散计算确定。此时，由总电流给出的上述磁场值就起到了边界值的作用。

2.4 磁通量计算

由于磁通量是一种面积分，或者是一种求和式，对于各个区域具有可加性。为了避免在计算中因为理想导体、非理想导体和磁腔区的差别(他们的电阻率相差较大)而出现“刚性问题”，对磁通量进行分区计算。即根据各区的性质和耦合情况，分别用总电流和理论公式，或者求解磁扩散方程组，得出各区的磁场分布和区磁通量，再相加得出总的磁通量。

3 算 例

3.1 平面磁驱动准等熵加载炸药实验

3.1.1 实验装置和计算模型

在CQ-1.5装置上进行了塑料黏结炸药JO9159的台阶靶样品ICE实验，炸药JO9159的质量比为w(HMX)∶w(黏接剂)∶w(钝感剂)=95∶4.3∶0.7。利用两个相近平行导体平面上的强电流与其感应磁场相互作用产生的洛仑兹力，对电极板施加压力脉冲，得到平滑上升的压力波，实现对炸药样品的准等熵加载。采用线圈测量流经负载的电流曲线；DPS测量炸药/窗口界面的粒子速度历程。图4(a)为在电极板上已安装好的炸药样品、窗口和DPS光测探头实验照片。

在磁驱动台阶靶构形实验中，由一发实验中获取多种厚度样品的自由面粒子速度历程，在建立模型时左右两边同时计算。模型结构示意图见图4(b)所示，不同物质为不同的计算区域，空腔单独作为一个区。其中，空腔距离设定为0.5 mm，LiF窗口为3 mm，两个炸药样品的厚度为0.599和0.789 mm，对应的铝电极板厚度为0.353 mm。铝和炸药均采用固态HOM状态方程，初始电压设定为60 kV。

图4 平面磁驱动准等熵加载实验Fig.4 The planar magnetic driven isentropic loading experiment

3.1.2 实验结果验证

图5(a)给出了同一发实验中两个厚度的炸药样品/LiF窗口界面粒子速度的计算结果与实验结果；图5(b)给出了回路放电电流的计算和测试波形。

从图5(a)中可看出，总体上实验值和计算值符合得较好，速度曲线的峰值和上升沿斜率基本一致。在实验和计算的结果中，均显示未产生冲击波，0.599 mm样品粒子速度峰值的实验值和计算值分别为304.6和303.1 m/s，0.789 mm样品粒子速度峰值的实验值和计算值分别约为332.8和335.4m/s，其最大偏差约为0.8%。厚样品/窗口界面的粒子速度峰值高于薄样品，这是因为薄样品加载还未完成，界面反射波就到达加载面改变了加载历史，而厚样品还处于加载状态，造成其速度峰值更高。粒子速度曲线的前沿上升时间偏差较大，最大约为8.4%，这是因为实验的粒子速度曲线起跳点附近不如计算结果平滑。这可能是多种因素造成的，一种可能的原因是实验系统的安装误差造成的，但更主要的原因可能是加载电流波形的变化造成的。从图5(b)中可看出，虽然电流周期和峰值的实验与计算结果基本一致，但实验波形受到回路各种耦合因素的影响，其电流1/4周期内的前半部分上升斜率明显高于计算得到的电流波形，而电流波形是决定加载压力的重要因素，因此实验中速度曲线上升初期的斜率较大。

图5 计算结果与实验结果Fig.5 Calculated and experimental results

3.2 炸药爆轰驱动圆柱套筒压缩磁通量实验

3.2.1 实验装置和计算模型

MC-1实验装置照片如图6(a)所示，它把炸药爆轰释放出来的能量转换为内爆导体套筒的动能，然后压缩其空腔中的预置磁通量，再转变为磁能。主要参数如下：炸药装药为RHT-901，其质量比为w(TNT)∶w(RDX)=4∶6，爆速约为7.79 km/s，厚度为55 mm；驱动套筒为不锈钢，外径为100 mm，厚度1.5 mm，其内壁镀40 μm的银；样品管为铜，外径为8 mm，厚度为1 mm；中心磁腔轴向初始磁场6 T。采用磁探针测量轴向磁场的变化和分布情况；采用小型激光干涉测速探头测量磁压力加载下样品管内壁径向速度的历史。

模型结构示意图见图6(b)，共分为5个计算区块：第1区为铜材料的样品管，厚度为0.1 cm；第2区为磁腔区，根据具体实验构型设置长度为4.75 cm，设置6 T的初始磁场；第3区为银镀层；第4区为不锈钢，厚度设置为0.15 cm；第5区为炸药装药区，起爆方式设定为右端5点起爆(模拟聚心爆轰)，反应速率采用Forest Fire律，爆轰产物状态方程采用气态HOM状态方程。炸药反应速率参数分别为：

C0=-1.035×101，C1=4.734×102，C2=-1.675×104，C3=4.476×105，C4=-8.493×106，

C5=1.156×108，C6=-1.140×109，C7=8.207×109，C8=-4.299×1010，C9=1.618×1011，

C10=-4.261×1011，C11=7.438×1011，C12=-7.279×1011，C13=3.617×1011；

状态方程参数分别为：

A1=-3.526，A2=-2.334，A3=0.597，A4=0.003，A5=-0.175，

K1=-1.561，K2=0.533，K3=0.081，K4=0.003，K5=-6.844×10-4，

Q1=7.503，Q2=-0.441，Q3=0.151，Q4=0.068，Q5=-0.024，

cV=2.09 J/(g·K)，e0=10 kJ/g。

图6 内爆圆柱套筒磁通量压缩实验Fig.6 Magnetic flux compression experiment of imploding cylindrical liners

3.2.2 实验结果验证

计算和实验得到的空腔磁场随时间的变化曲线，样品管内表面速度历程曲线，见图7。

由于在小管道中的测试比较困难，同时内爆收缩后期重复性较差，因此即使在相同状态下，两发实验值的磁感应强度峰值仍相差约254 T(见图 7(a))。计算得到的磁感应强度峰值约为1 489 T，如果以两发实验的平均值作为基准，计算的磁感应强度峰值与实验结果偏差约为13.1%；如果取两发实验中的最大值1 443 T作为基准(偏低的实验结果可能是由于磁探针损坏而没有记录到更多的信号所致)，计算值与实验结果非常接近，偏差约为3.1%。图7(b)显示计算的样品管内表面速度历史曲线，与实验测试曲线基本一致，速度从零开始平滑上升，计算得到的峰值约为6.87 km/s，两个通道的测量结果峰值分别为6.79和5.98 km/s，测试结果之间偏差11.9%，这也许与安装精度和测试误差有关。如果以两个测量结果的平均值作为标准，计算值比实验值高7.5 %左右。总体说来，计算结果与实验测量结果符合较好。

图7 计算结果与实验结果Fig.7 Calculated and experimental results

4 结 论

改造后的SSS-MHD编码将力学样品和电磁学负载合二为一，能够统一处理平面和圆柱面等多种磁流体力学实验构形，并能同时进行磁流体力学和含能材料的反应流动计算。通过SSS-MHD编码，对平面准等熵压缩炸药和炸药爆轰驱动圆柱套筒压缩磁通量这两类典型磁流体力学构形实验进行的数值计算表明，电流和磁场变化、速度历史等计算结果与直接测量得到的实验结果符合。

[1] Robinson A C, Brunner T A, Carrol S, et al. ALEGRA: An arbitrary Lagrangian-Eulerian multimaterial, multiphysics code[R]. AIAA 2008-1235, 2008.

[2] Lemke R W, Knudson M D, Robinson A C, et al. Self-consistent, two-dimensional, magnetohydrodynamic simulations of magnetically driven flyer plates[J]. Physics of Plasmas, 2003,10(5):1867-1874.

[3] Knudson M D, Lemke R W, Hayes D B, et al. Near-absolute Hugoniot measurements in aluminum to 500 GPa using a magnetically accerlerated flyer plate technique[J]. Journal of Applied Physics, 2003,94(7):4420-4431.

[4] Knudson M D, Hanson D L, Bailey J E, et al. Principal Hugoniot, reverberating wave, and mechanical reshock measurement of liquid deuterium to 400 GPa using plate impact techniques[J]. Physical Review B, 2004,69(14):144209.

[5] Oliphant T A, Witte K H. RAVEN[R]. NM: Los Alamos National Laboratory, 1987.

[6] Sheppard M G, Fowler G V, Freeman B L. Generation of ultra-high magnetic fields for AGEX[R]. NM: Los Alamos National Laboratory, 1994.

[7] Kozlov M B, Aseeva V V, Boriskov G V, et al. Computational investigational of operating conditions of the cascade generator MC-1 with large HE charge[C]∥Proceedings of 28th IEEE International Conference on Plasma Science and 13th IEEE International Pulsed Power Conference. Las Vegas, NV, 2001:1182-1184.

[8] Frese M H. MACH2: A two-dimensional magnetohydrodynamic simulation code for complex experimental configurations[R]. AMRC-R-874, 1987.

[9] MacGillivray J, Peterkin R E, Burke D A, et al. MACH3: A CHSSI code for computational magnetohydrodynamics of general materials in generalized coordinates[C]∥IEEE Computer Society. Proceedings of the users group 2003 conference. Bellevue, WA, 2003:9-13.

[10] Frese S D, Frese M H. Recent improvements to MACH2 and MACH3 for fast Z-pinch modeling[C]∥Proceedings of 5th International Conference on dense Z-pinch. Albuquerque, NM, 2002:380-383.

[11] L’Eplattenier P, Cook G, Ashcraft C, et al. Introduction of an electromagnetism module in LS-DYNA for coupled mechanical-thermal-electromagnetic simulations[J]. Steel Research International, 2009,80(5):351-358.

[12] 王刚华.磁驱动准等熵压缩和高速飞片实验、计算和反积分数据处理技术[D].绵阳：中国工程物理研究院，2008.

[13] 张恒第.MC-1型爆磁压缩发生器一维内爆:磁流体力学模拟[D].绵阳：中国工程物理研究院，2012.

[14] 张旭平.磁驱动超高速飞片发射技术研究[D].绵阳：中国工程物理研究院，2012.

[15] 孙承纬.一维冲击波和爆轰波计算程序SSS[J].计算物理，1986,3(2):142-154. Sun Chengwei. SSS: A code computing one dimensional shock and detonation wave propagation[J]. Chinese Journal of Computational Physics, 1986,3(2):143-154.

(责任编辑 丁 峰)

Magnetic hydrodynamic calculation on planar and cylindrical configurations

Zhao Jibo, Sun Chengwei, Gu Zhuowei, Wang Guiji, Cai Jintao

(NationalKeyLaboratoryofShockWaveandDetonationPhysics,InstituteofFluidPhysics,ChinaAcademyofEngineeringPhysics,Mianyang621999,Sichuan,China)

In this work we developed a one-dimensional elastic-plastic reactive magnetic hydrodynamic code SSS/MHD with multi-medium, multi-cavities and coupling with external circuits on experimental configurations, based on the Lagrange code SSS. The key problems in coding were resolved through defining the configuration boundary coupling with circuits and unifying the treatment of the magnetic field boundary conditions of various configurations. The advantages of our coding are mainly shown in that, on the one hand, the various magnetic hydrodynamic experiment configurations-including both planar and cylindrical ones-can be handled in a unified fashion and, on the other, that the magnetic and reactive hydrodynamic calculation can be simultaneously conducted. Thus, following our way of coding, two typical experiments were calculated, i.e. the magnetic driven isentropic compression planar explosive and the explosive detonation driven liner for flux compression, the results of which accord with the experimental results.

fluid mechanics; calculation code; multi-field couple; isentropic compression

10.11883/1001-1455(2016)01-0009-08

2014-07-16； < class="emphasis_bold">修回日期：2014-12-21

2014-12-21

国家自然科学基金NSAF联合基金重点项目( 11176002)

赵继波(1977— )，男，博士，副研究员，zhaojibo77@hotmail.com。

O361.3 <国标学科代码：1302547 class="emphasis_bold"> 国标学科代码：1302547 文献标志码：A国标学科代码：1302547

A