电磁驱动高能量密度动力学实验的一维磁流体力学多物理场数值模拟平台：SSS-MHD*

孙承纬，陆 禹，赵继波，罗斌强，谷卓伟，王桂吉，张旭平，陈学秒，周中玉，李 牧，袁 红，张红平，王刚华，孙奇志，文尚刚，谭福利，赵剑衡，莫建军，蔡进涛，金云声，贺 佳，种 涛，赵小明，刘仓理

（1. 中国工程物理研究院流体物理研究所，四川 绵阳 621999；2. 上海激光等离子体研究所，上海 201800；3. 深圳技术大学先进材料测试技术研究中心，广东 深圳 518118；4. 深圳技术大学大数据与互联网学院，广东 深圳 518118；5. 中国工程物理研究院化工材料研究所，四川 绵阳 621999；6. 中国工程物理研究院应用电子学研究所，四川 绵阳 621999；7. 中国工程物理研究院，四川 绵阳 621999）

极端物理学顾名思义是极端条件下的物理学，但按当今科学前沿的实际理解是高能量密度状态下的物理学，即高能量密度物理[1]。高能量密度状态的公认条件是受载过程中物质内能的增量大于0.1 MJ/cm3，即其热力学（静水）压力超过100 GPa[2]，或者波长约1 μm 的辐射能束对物质的辐照度大于3 PW/cm2。数十年前高能量密度状态显然只能在天体、行星、地球数千千米深部等环境以及核武器物理过程中存在，但产生高功率短脉冲激光的啁啾技术和美国Sandia 实验室100 TW 级功率低阻抗强电流驱动器Z 机器的问世，使这种状态已可在不少实验室内实现。研究高能量密度状态下物质、辐射及其相互作用的多学科交叉的极端物理学前沿领域，对惯性聚变、天体物理、行星物理和核武器物理研究具有重要现实意义。

电磁驱动下宏观物体的高能量密度动力学实验，包括极端材料动力学、流体动力学行为和高温高密度等离子体产生等，均具有重要的技术背景和应用价值[3-4]。从驱动原理上看，这种实验可分为2 类：磁驱动和磁压缩。磁驱动系利用导体构形表面上强驱动电流脉冲与自生磁场作用的电磁力加载，实现材料准等熵（斜波）压缩或发射超高速固体飞片（统称isentropic compression experiment, ICE），可用于冲击压缩、驱动金属套筒内爆、形成高温高压等离子体等。磁压缩系利用炸药爆轰或电磁驱动金属套筒的低速内爆运动，把较大体积中的磁通量压缩到小区域中，形成很高磁场和磁压，把该小区域内的材料（尤其是体积较大的低密度材料）准等熵压缩到高密度。利用炸药爆轰为动力的磁通量压缩装置，又称为磁通量压缩（MC-1）发生器。这2 类实验均涉及极端条件下实验构形的力学与物理学耦合问题，难以从理论方程组上完善考虑，但能通过数值模拟途径用计算解决。

1980 年代，美国Los Alamos 国家实验室研制了一维辐射磁流体力学编码RAVEN[5]，具有一维任意拉格朗日-欧拉（arbitrary Lagrangian-Eulerian, ALE）功能，可以处理带真空磁腔的负载构形。RAVEN 编码利用算子分裂方式，解决了流体弹塑性材料、双流双温等离子体和高温辐射物质等强电磁加载下多物理场耦合模拟问题。目前，公认的最先进磁流体力学（magnetohydrodynamics, MHD）构形计算编码是美国Sandia 国家实验室的多介质多物理场ALE 编码ALEGRA 系列。1990 年代初，Sandia 实验室将高精度冲击动力学有限差分（欧拉）编码CTH 与改进的拉格朗日有限元编码结合，得到在无结构单元有限元网格上表述的ALE 编码，即ALEGRA 最初的版本。经过30 多年的开发，ALEGRA 已成为功能强大的计算辐射磁流体-冲击动力学-多物理场耦合的ALE 有限元编码，得到了广泛应用[6]。

本文中叙述的一维冲击爆轰和磁流体力学多物理场计算编码SSS-MHD，可对电磁驱动高能量密度动力学实验开展磁流体力学多物理场数值模拟，分析强电流或者炸药驱动的瞬态电磁场加载、受载实验构形和材料样品中发生的弹塑性磁流体力学运动和热力学变化等过程。SSS-MHD 编码的基础是1980 年代流体物理研究所研发的冲击爆轰动力学编码SSS[7]，后又扩展了激光辐照效应[8-9]、等离子体电磁驱动的计算功能[10-11]。随着磁驱动和磁压缩实验的开展[12-13]，SSS 编码也向MHD 扩展，并于2012 年提出SSS-MHD 编码的初稿[14-16]，10 年来经过许多实例计算核对和修改完善后基本定稿[17-22]。图1 所示的各类实验装置基本概括了具有重要技术应用背景的电磁驱动实验类型，均可用SSS-MHD 编码进行不同程度的模拟。

图1 适合SSS-MHD 编码模拟的各种类型电磁驱动高能量密度物理实验Fig. 1 Kinds of the magnetically-driven high-energy-density physics experiments suitably simulated with the SSS-MHD code

多物理场耦合的数值研究，主要体现在多个物理系统方程组在离散化计算中进行恰当的（数学上完备的）相互耦合计算。SSS-MHD 编码主要理论框架是弹塑性磁流体力学方程组，自身就是力学与电磁学两大系统的耦合，进一步还能扩展到与等离子体物理和辐射流体力学等学科的耦合。由于使用场景和目的不同，磁流体力学方程组形式多样。图1 的各类实验涉及低频脉冲强电流和强磁场与导电体的短时间非相对论性相互作用，在电磁力和能量表达式中磁场的份额远高于电场。这种意义下的磁流体力学近似，只需要2 个电磁学变量，如电流密度和磁感应强度，相关的方程式为磁扩散方程和欧姆定律。上述耦合体现于运动方程中的电磁力、能量方程中的焦耳热，同时磁扩散方程与介质运动速度有关，这种架构为SSS-MHD 编码结构的简化创造了有利条件。

本文组成如下。第1 节，叙述弹塑性磁流体力学多物理场耦合基本方程组；第2 节，建立作为SSSMHD 编码基础的拉格朗日一维方程组的具体形式；第3 节，叙述SSS-MHD 编码的创新结构和功能；第4 节，展示平面构形ICE 实验的SSS-MHD 模拟，包括铝、钽、PBX 炸药的准等熵斜波压缩和铝飞片高速发射实验；第5 节，展示圆柱面构形的固体套筒电磁内爆实验及炸药内爆磁通量压缩（MC-1）发生器实验的模拟计算；第4～5 节将SSS-MHD 编码计算与Sandia 实验室Z 装置的3 个高水平实验和计算实例进行比较，而且模拟流体物理研究所CQ、CJ 系列装置的4 个典型实验；第6 节的结束语中指出，上述算例表明SSS-MHD 模拟与实验数据以及ALEGRA 计算数据的相对偏差基本不超过5%。本文中的数值模拟还以丰富的物理量剖面图显示实验过程全面的物理图像。编码SSS-MHD 的创新特色以及原有的冲击动力学编码的牢固基础，使其模拟通用性强、功能灵活全面、计算速度快、精度高。编码SSS-MHD 的开发不但能为极端物理学、极端材料动力学实验提供有力的数值模拟平台，而且有助于多维MHD 多物理场编码的研发。

1 弹塑性磁流体力学基本方程组

连续介质力学运动的基本规律为质量、动量和能量守恒定律，从一般的连续介质动力学方程组可以方便地导出磁流体力学的基本方程组[26]。首先，在动量方程和能量方程中分别增加与电磁力和电阻产生热量（焦耳热）相关的项，把向量形式的连续介质动力学方程组改写如下，依次为连续性方程（质量守恒）、运动方程（动量守恒）和能量方程（能量守恒）：

将式(2)全部点乘以u，再与式(3)相加，得到以总比能E=e+u·u/2 表达的能量方程：

式(4)是SSS-MHD 编码实际使用的能量方程，其优点是简化了张量运算。后文中，能量方程中无外体力F项，外源项Q写成WL即单位质量介质吸收的激光辐射功率。

在冲击动力学或爆炸力学中，对连续介质通常采用流体弹塑性的力学模型加以描述。将应力张量Σ分解为 Σ=-pI+T，写成分量式：

式中：比容v=1/ρ 。SSS-MHD 编码计算中材料EOS 基本均采用完全物态方程数据库SESAME[27]。SESAME 库中包括热力学和低温等离子体性质，为SSS-MHD 编码多物理计算提供了重要基础。由式(5)～(6)和未写出的本构关系，温度以及应力张量各个分量都可以通过微元的比内能、比容、粒子速度分量及其导数或积分表达，此时未知的力学函数尚有e、 ρ 和u，式(1)～(3)恰好是2 个标量和1 个矢量方程，只要其中的电磁量J和B能借助于电动力学方程组耦合求解，即可实现方程组数学上的封闭。

考虑电动力学之始，先假定介质中不存在自由电荷，无需考虑位移电流和电介质的极化或磁化，并得到电场强度E和磁感应强度B都是无源的： ∇·E=0 ， ∇·B=0 。麦克斯韦方程组还有如下定律。

法拉第定律：

安培定律：

再增加广义欧姆定律：

前面方程组已含有粒子速度u为未知函数，式(7)～(9)提供了求解E、J和B等3 个电磁量的方程，其中 µ 为介质的磁导率。由于安培定律，式(1)～(3)中仅出现磁场B，在该意义下系统被称为磁流体力学方程组。进一步化简，可从式(7)～(9)中消去E和J，得到决定B的磁扩散方程：

至此，式(1)～(6)、(8)、(10)以及未写明的介质材料本构关系，组成了一般情况下多维弹塑性磁流体力学的封闭方程组，而且是用矢量、张量和微分算子表述的一般形式，提供了写出任意正交曲线坐标系中分量形式微分方程组的方便途径。

2 一维拉格朗日形式弹塑性磁流体力学方程组

根据第1 节中的基本方程组，可写出一维情形下平面、柱面、球面坐标系中欧拉形式的弹塑性磁流体力学方程组（分量式），如下。连续性方程：

运动方程：

能量方程：

在一维磁流体力学情形中，几何指数N只能取0 或1，SSS-MHD 编码的电磁学方程组（式(9)～(10)）的分量形式如下。安培定律：

磁扩散方程组：

式中：Cθ=rNBθ称为磁感应强度的广义 θ 分量，在柱对称情形中有重要意义。

在欧拉-拉格朗日形式的转换中，时间t保持不变，拉格朗日空间坐标采用质量坐标M，即计算构形自最左边界点到格点处的累计质量。欧拉空间坐标r成为格点拉格朗日坐标M的位置函数，位于r处格点的质量 ∆M则是厚度为dr、底面积分别是单位面积（平面）、r× 单位宽度（柱面）和r2（球面）的（单位弧度和立体角）体积中包含的介质质量。因此，2 个坐标系之间的转换关系是：

为了加以区分，将拉格朗日形式方程组中的r改写为R，粒子速度u改写为U，径（纵）向应力 σrr简写为 σ 。每个拉格朗日格点的质量保持不变体现了质量守恒，也省去了连续性方程。为此，需要添加一个速度方程，即U等于R的时间偏导数。

由于高速度、高压力问题参数范围的特点，冲击爆轰编码SSS 采用了“冲击波单位制”，其基本单位是质量g、长度cm 和时间μs，因而速度单位是cm/μs 即10 km/s，压力单位是g/(cm · μs2)=100 GPa。前面推导的电磁学方程组采用国际单位制（SI），因此在耦合方程组中电磁量相关项需添加修正系数，达到与力学量单位的自洽。磁导率 µ=µ′µ0，其中真空磁导率 µ0=4π×10-7H/m ， µ′为介质的相对磁导率。

经过上述推导和转换，得到SSS-MHD 编码使用的拉格朗日形式一维弹塑性磁流体力学方程组如下。

连续性方程和速度方程：

运动方程：

能量方程：

安培定律（Cθ≡rNBθ）：

磁扩散方程组：

对式(18)～(23)做离散时，时间离散格式是简单的前向差分，空间离散中格点序号自左向右排列，除了每一格的位置、速度、电流密度定义在格子中点以外，其他力学量和电磁量均定义在格子右点。因此，速度方程和安培定律是中心差分格式，其余方程都是前向差分格式。整个差分方程组的具体形式可参看文献[5]。

3 SSS-MHD 编码的结构和功能

3.1 计算编码的结构

SSS-MHD 编码的结构立足于SSS 程序原有的框架[7]，保持弹塑性流体动力学计算的各个子程序模块构架，进行适应MHD 计算的改造。首先，在导程序中增加MHD 建模的内容，在总卡（GENERL）和区卡（COMPNT）中引入与MHD 计算相关的数组、变量和指数等，然后，在专门的MHD 建模子程序（MHDZON）里再对有关区卡赋以MHD 负载构形的类型指数、构形参数以及初始条件等，完成MHD 计算建模。接着，在相关子程内部和计算循环节点处加入控制MHD 计算的语句，如在计算循环开始处设定边界条件状态及参数的同时，检查和调整空腔和磁腔的状态及参数。

随着力学计算流程的进行，在运动方程子程序（VELOC）、能量方程子程序（ENERGY）和热输运计算子程序（HETCON）中，使用上时刻电磁变量值分别做洛伦兹力和焦耳热计算。力学计算完成之后，流程进入专设的电磁计算程序包（MHDBLK，后详）。最后是数据更新（UPDATE）、后处理等，进入下一个时间循环，直至计算时间t达到设定时刻ts。编码的计算流程见图2。

图2 SSS-MHD 编码的计算流程Fig. 2 Flowchart of the SSS-MHD code

原编码SSS 是一个框架式的通用程序，具有多种物理场耦合计算的功能，可以包容或扩展各种物态方程和物性数据库，加入有关局部性物理过程的计算（例如激光与物质的热和冲量耦合、平衡电离的Saha 方程、动态断裂的NAG 模型和含能材料的燃烧爆炸反应等），适合平、柱、球对称几何的动力学问题，设有多种边界条件和初始条件可供选择，能方便地实现冲击或爆轰的加载模拟[7]，这些特点已全部为SSS-MHD 编码沿用。下面各小节将具体阐述SSS-MHD 计算编码的关键内容和特色。

3.2 MHD 负载构形的统一建模

在磁流体力学计算模型的建模中，由于力学与电磁量具有不同方向的坐标分量，首先需要设定力学构形和相应的加载磁场类型，然后赋予各个MHD 计算区域及格点处的电磁学物性参数和相应初始值，必要时在整体模型边界处补充电磁学边界条件。其中的实质性工作是负载构形及磁场的分类和表述。

图3 是SSS-MHD 编码规定的实验负载构形种类，N是一维连续性方程式(11)中的几何指数。这样的一维几何中介质速度只有径（纵）向分量，但电流和磁场可以有2 个方向，即轴向（z）和角向（θ）分量。磁扩散相应地有2 个方向的方程式，还可以用2 个驱动电路分别提供z、θ 方向的加载电流。一维磁流体力学原则上可以计算z、θ 及其组合的螺旋磁力线相关问题。负载构形是各种材料的多层平行平板或者同轴圆柱壳层组成的几何结构，其间可以包括空腔（有磁场的空腔称为磁腔）。建模依据主要的电磁驱动机制，对构形加载或被其压缩的磁场（磁力线）有轴向（z）或角向（θ）以及两者组合（z、θ 螺旋）等3 种形式。螺旋磁场中轴向与角向分量各自成体系、并行不悖，虽然图中未画出这个类别。从计算体系来说，只有图3 的4 种基本构形，其中平面z和θ 等2 种构形实际相同，构形的名称表示其加载电流方向。构形的板、壳可以是多层组合，也没有具体画出。

图3 SSS-MHD 编码中实验负载构形的种类（红线表示回流导体）[29]Fig. 3 Types of experimental configurations in the SSS-MHD code(red curves standing for return conductors)[29]

图1 中各类实验的构形和磁场均可抽象为图3 的类型。负载构形中的磁场可以是加载电流感生的（持续励磁），也可以是初始时由外界给定的（瞬态种子磁场或定常励磁）；推动构形运动的可以是加载电流自身形成的洛伦兹力（磁驱动），也可以利用炸药爆轰或冲击的力量（磁压缩）。

3.3 空腔和磁腔

SSS-MHD 编码沿用了SSS 编码的空腔功能，并扩展到带磁场的磁腔情形。计算模型通常为顺序排列的介质分区所组成的单连通区，如果队列中存在真空间隔并不妨碍欧拉形式计算，但对于拉格朗日计算则发生“零质量格点”的困难，通常需采取分区分阶段计算合成或实行一维ALE 策略，这样可能使得编码复杂、机时增多，计算精度也受影响。

SSS 编码采取物质区与空腔区一体化拉格朗日计算的方法。图4 表示把空腔（区I）处理为只有3 个零质量格点（J至J+2）的特殊区，不参加物质区运算，但其左右两岸点J和J+3 分别具有左右物质区岸点位置坐标R的数据，空腔间隙等于R(J+3)-R(J) 。若空腔间隙大于零，左右两岸点J和J+3 处均取自由边界条件，空腔中间两点J+1 和J+2 分别是J和J+3 点的鬼点（影子点）。如果空腔间隙为零或负值，意味着空腔闭合甚至两岸碰撞过冲，需要退到前一时刻，适当减小时间步长重算，使得重算的空腔间隙在规定误差下为零，并使点J和点J+3 的R值保持相等，成为同一个右岸介质内点。如果空腔原来是闭合的，但表观重合的左、右岸点相互作用的拉力已变为超过一定阈值的斥力，则空腔将被自动拉开。对于空腔开、闭变动的暂态振荡过程，程序给予致稳处理。大量算例表明，SSS 编码的空腔一体化计算功能比较成功，步骤简便，精度较高。

图4 SSS 编码中空腔区及邻区中格点的设置Fig. 4 Meshes in the cavity and nearby in the code SSS

磁流体力学计算中磁空腔较常见，其力学运动可沿用SSS 编码的处理方法，但需增加电磁学计算，并保持两岸点磁场及磁力对称相等的处理。经过改进，SSS-MHD 编码实现了一维磁流体力学多连通区的拉格朗日一体化计算，具备了对于较复杂实验及装置进行数值模拟的能力。

磁腔处理方法基于电磁学一维问题的特点并利用了SSS 编码的基础。如图3 所示，运用积分形式的安培定律和法拉第定律可以证明，实验构形中磁腔界面处的磁场总是与相应电流分量正交的，磁腔左右岸导体表面的电流只能是面电流。空腔（或电介质）中磁场的生成，无外乎馈入加载电流或初始种子磁场以及构形导体中产生感应电流等几类原因，都可由磁扩散方程组计算确定。由于不考虑静电电荷和位移电流，电介质的电磁学性质等同于磁腔。如果最右端的回流导体质量较大或可不加考虑，则可简化为最右端的单纯磁腔（力学上是自由边界，不是空腔）。此时的计算模型称为单边驱动，加载电流线只要一条，还需设定右端磁腔条件来计算该腔的磁通量。圆柱面θ 构形可能存在轴线处的中心磁腔，但只能是环形电流驱动的均匀轴向磁场，如MC-1 装置的情形。

3.4 磁腔磁场的计算

图3 表明在平面一维情形（N= 0）、以一对无限大平行电极板为界的平磁腔（即阳阴极间隙AK-gap）中，加载电流可以有2 个方向，分别产生均匀磁场Bθ或Bz，磁力线均为直线；圆柱面一维（N= 1）几何具有θ 或z等2 种构形，同轴的无限长双圆柱筒的间隙或单柱筒内部均为圆（环）磁腔，可近似推广到双多边形柱筒之间的环形磁腔（环腔）情形。因加载电流有θ 或z方向，环腔中磁力线则为轴向直线或同轴环线。按照图3 画出的磁力线围线，由安培定律给出平直腔和圆（环）腔的磁场与构形加载电流的关系，得到平直磁腔和圆（环）磁腔中的磁场分别为：

式中：I为构形的总加载电流，i为该电流经过单位宽度表面的平均电流（线密度），w为腔宽度或环腔的周线长度，下标表示电流或磁场的方向。按照一维几何的原意，平直腔是无限宽的，总电流I无限大，式(24)是用有限宽度平均电流近似无限宽的理想状况，带来的误差修正下面说明。式(25)适用于圆腔，由于非圆形环腔也是封闭的，可以近似采用圆环腔来处理，从而环腔的直线段部分相当于理想化的平直腔，虽然其电流密度因需按周线长度平均计算而大为下降。

柱筒构形按式(25)的近似很接近实际的柱面二维情况，误差不大。然而，理想平面一维的式(24)与实际构形的三维电磁场以及实验结果的差异就较显著。式(24)第一个等号右边的因子 (I/w) 即以磁腔实际宽度w计算的平均加载电流密度，相当于用电流密度 (I/w) 、宽度无限的理想磁腔做近似。而如图3所示，实际的平面一维“单条片”构形（磁驱动实验的一对加载电极板）是“折叠”后供加载电流来回的狭长金属片，其间的AK-gap 即是磁腔。从能量角度看，式(24)的近似显然是对磁腔中电流密度、磁场及磁压的一种高估。陆禹[30]利用商业电磁学软件[31]进行了二维平面电磁场计算，以频率0.35 MHz 的正弦电流加载于一对铜电极板（截面宽2 mm、厚1 mm），两板间的磁腔截面宽度（w）为2 mm、间隙（g）为0.1 mm。此无限长“单条片构形”电磁响应的二维平面高频计算结果表示，加载电流（幅值）密度分布明显不均匀，集中于有限尺度磁腔截面的4 个角点，证实磁腔中部电流密度低于边部。这种计算提供了对有限尺度平直磁腔中电流密度分布不均匀性的定量估计。定义K为等效电流因数，即利用商业电磁学软件模拟计算得到的平面磁腔中部幅值电流密度与全宽度平均幅值电流密度的比值。将式(24)（不计下标）修正为：

式中：修正后的电流密度iK=KI/w。平面一维磁腔的等效电流因数K是一个小于1 的正数，与磁腔几何参数w、g以及加载电流参数有关，如果w→∞,g→0 ，则K→1 。

严格说来，上述等效电流因数只反映实验开始时的情况，由于运动中构形动态电感的变化需要进行难以做到的至少是二维的磁扩散计算，才能对电流分布作更适当的考虑。动态电感的增大降低了加载电流的幅度，也可看作为运动过程中因数K的继续下降，然而实际上MHD 计算包括了这部分的相互作用（详见3.6 节），使得K动态变化的影响大为减弱，绝大部分计算中只需沿用初始的K值，即可得到全程与实验数据相符的良好结果。这就是SSS-MHD 编码处理平面磁腔的方法。

本节中论述的平面磁腔磁场近似计算方法对于SSS-MHD 编码的磁扩散计算具有重要意义，因为磁腔中的磁场Bz或Cθ都是常值，根据式(24)～(26)可由构形加载总电流的即时值直接确定。该总电流（时间函数）可以作为输入数据给定，也可从电路方程组计算获得。因而，这些磁腔的磁场值可以用作为磁扩散方程计算中可靠的定标值，避免了隐式计算的麻烦和误差。该问题在理论上有更严谨的方法，如Lemke 等[32]利用Poisson 方程计算了含有四边形环形磁腔的多连通区平面二维磁流体力学问题，但并不能得到有普遍意义的结果。

3.5 显式计算格式

SSS-MHD 编码主要用来模拟MHD 实验构形的力学和电磁运动，进而研究样品材料的物态方程和材料动力学问题。为了满足数值计算要求的稳定性条件（如SSS 编码中的Courant 条件），提高计算精度，并可防止邻近格点尺度差距过大引起的刚性问题，希望采取一体化的拉格朗日显式运算。力学方程组可以从给定边界条件的左端点开始做显式计算。把磁扩散方程同样离散为显式有限差分格式，但其边界条件不在左、右端点给定，需要通过迭代使得本时刻的整体磁场分布、磁通量、反电动势（构形的动态电感）与加载或励磁电流值达到高度自洽，才能获得较高的计算精度，具体做法见3.7 节的说明。

避免计算不稳定性和刚性问题的重要措施是编码中介质区可选择多种方法进行分格。计算值的振荡往往发生于几何或物理因素造成邻近格子尺度或物理量急剧变化的区域，如内爆聚心阶段的中心范围、点爆炸和散心爆轰开始阶段的核心区域、高速碰撞的接触面附近、脉冲电流的集肤层、高功率激光的吸收层和烧蚀层等。如果把这些区域做小尺度分格，并平滑过渡到周围正常尺度分格的邻区，则有可能避免振荡，得到物理上合理的计算结果。SSS-MHD 编码中除了通常的等步长、等质量分格方式外，还有在指定分区中使格子尺度（或质量）按一定比例增大或缩小的分格方式。激光吸收层处理和从冷问题开始的等离子体演化等SSS 编码计算，就是成功实例[8-9]。

3.6 SSS-MHD 计算与驱动电路方程的耦合

负载构形的动态电感变化与驱动器放电过程之间存在重要的相互作用，若要理解这些复杂过程，首先应把MHD 计算与驱动器电路计算进行耦合，建立准确预测实际加载电流的能力。与驱动电路的耦合是MHD 计算能否解决实际问题的一个关键，也是SSS-MHD 编码的重要特色。这个要求至今仍是对强电磁驱动实验计算模拟能力的重大挑战。

简单实验装置中驱动电流的传输及成形通常由RLC 电路描述，其计算是常微分方程组的初值问题。复杂的大型装置使用脉冲形成介质线技术，可由一阶偏微分方程组描述，并可离散化为线性代数方程组求解。在与磁流体力学计算耦合方面，两者没有实质性差别。下面以当前使用的集中参数RLC 电路为例，说明与MHD 计算的耦合的方法。先把RLC 电路写为：

式中：第1 个等号左边为电容器组的剩余电压，Q和C0分别为其剩余电荷量和初始电容值；U0和I分别为充电电压和负载电流；L0、LL、R和VS依次为外电路（负载构形）的固定电感、其他可变电感、固定电阻和开关电压降（伏安特性），均在数据文件中输入；ε 为负载构形对于驱动电路的反电动势，同一构形可以有z、θ 等2 个驱动电路和2 套电路参数。

驱动电路与负载构形耦合的关键是计算反电动势ε，即构形路端电压（电感性）VAK和沿其电流回路的电阻性电压降VLP之和。根据积分形式的法拉第定律，VAK为构形总磁通量 Φ 的时间导数，即VAK=dΦ/dt；根据欧姆定律，有VLP=ηlLPI（或相应的积分式），lLP为构形电流回路长度。任何时刻，电容器组剩余电压减去各项电压降之后应等于负载构形的反电动势。从物理上说，VAK就是构形磁通量变化在围绕其磁场的导体围线中感应的路端电压。因此，电路耦合不但要求解电路方程，更需要求解磁扩散方程组得到构形中的磁场分布，进行分区和全构形的磁通量计算。

3.7 SSS-MHD 计算程序包MHDBLK

力学计算之后就可计算在上一时刻应力、电磁力和各种能量作用下、在已发生力学变形并且即时负载电流作用下的构形中电磁场的变化，也就是求解本时刻的电路方程和磁扩散方程组，这就是SSSMHD 编码中体现MHD 多场耦合的程序包MHDBLK 的工作（图5）。它的功能有3 项：(1)电路方程子程序（DLRK），依据上一时刻（j）已迭代n次的反电动势ε 和负载电流I值解算电路方程（调用子程CURRNT 可设定该方程形式），得到负载电流n+1 次迭代值In+1；(2)磁扩散方程组子程序（MAGDFU），依据算出的In+1，求得磁场定标点处值，计算全部构形范围的磁场分布；(3)依据得到的磁场分布，计算负载构形n+1 次迭代的总磁通量及反电动势ε 的迭代值。若ε 迭代收敛不符合要求，则把n+1 次I和ε 的值输入DLRK 子程序，进行下一次迭代计算，直到第n次和第n+1 次的反电动势差值小于千分之一，才输出本时刻j+1 的负载电流和磁场分布的确定值。

图5 SSS-MHD 编码中计算程序包MHDBLK 的组成Fig. 5 Flow chart for the routine package MHDBLK in the SSS-MHD code

若无法获得驱动器电路方程的资料，加载电流也可选择指定电流馈入方式，如实验测量的实际负载电流或者另行估算的负载电流。此时，指定电流的数据以电流函数（CURFOM）形式馈入，直接把本时刻j+1 的指定电流值输入MAGDFU，算得磁场后直接输出，与电路无关，不进行迭代。SSS-MHD 编码整体为显式运算，但通过同一时刻的负载电流与磁腔磁场的自洽需通过迭代计算，可得到与隐式运算相同的效果，而且精度更高。按照这种方法，把上时刻j到本时刻j+1 的迭代扩大到整个力学计算范围也是可行的，但增加如此多的计算量对提高精度不一定有意义。

4 SSS-MHD 编码应用于平面构形磁驱动实验的典型算例

SSS-MHD 编码针对各类实验研究以及国外高水平实验结果，进行了大量验算和改进工作。本节列举的平面构形磁驱动实验（ICE）算例分为2 类：材料样品准等熵压缩和高速金属飞片发射。

图6 是上述2 类磁驱动实验构形的原理示意图，图6(a)～(b)都是单条片（strip-line）构形；图6(c)是方柱筒构形，其芯柱（阴极）和侧面板（阳电极）相当于4 个单条片，它们并联构成了四边筒形立体结构。样品自由面或窗口界面速度历史的测量采用激光干涉技术，如VISAR（适用于任意反射面的速度干涉仪）、PDV（光子多普勒测速仪）等。当材料样品受到很强压缩脉冲时，其自由表面会发生微喷射现象，干扰激光干涉测速过程，这种情况必须采用力学阻抗适当的光学单晶体片（LiF、蓝宝石等）作为窗口，严密覆盖样品表面以抑制喷射，测量的数据则是样品/窗口界面的速度历史。

图6 磁驱动准等熵压缩实验的典型负载构形Fig. 6 Typical loading configurations for magnetically-driven isentropic compression experiments (ICE)

4.1 磁驱动平面材料样品准等熵压缩实验（算例1～4）

算例1 是对美国Sandia 实验室Z 装置实验Z-1220[33]的验算，该实验采用类似图6(c)的6061-T6 铝合金4 层方筒构形，将一个侧面（阳电极板）加工为4 段样品，厚度h分别为0.616、0.868、1.118 和1.371 mm，相当于一个4 阶的台阶靶。每段样品外表面粘有镀增透膜的LiF 光学窗口。样品与窗口的界面速度用VISAR 测量。实验加载电流峰值约18 MA，上升沿约300 ns。美国Z-1220 实验数据以及ALEGRA-1D 编码与SSS-MHD 编码的计算模拟都示于图7。

图7 对于美国Sandia 实验室磁驱动等熵压缩实验Z-1220[33]的SSS-MHD 模拟（算例1）Fig. 7 SSS-MHD simulations for theisentropic compression experiment Z-1220[33]at the Sandia National Laboratories, USA (example 1)

SSS-MHD 计算是依据图7(a)的加载电流实验波形进行的，此例的矩形环腔周长5.8 cm，换算得知平均电流密度峰值约3.1 MA/cm，等效电流因数K为1。4 种厚度样品的2 种数值模拟速度曲线与实验数据的符合程度都较好。SSS-MHD 计算给出的材料样品最高压力为51 GPa，文献[33]给出的最高纵向应力为54 GPa。2 种模拟的差别在于SSS-MHD 计算在弹塑性转变阶段速度起跳过早，其原因可能是该计算的样品材料是纯铝而不是Z 实验的6061-T6 铝合金，此外还需要考察计算中输入电流波形初始50 ns 段落的细节以及铝样品中磁场扩散的情况，进而或可做出相应改进。对于模拟速度曲线末期与实验数据的偏差，主要是因为铝的弹塑性行为与模拟选用的材料动态本构模型之间存在差异。

算例2～4 是对流体物理研究所CQ 系列脉冲功率装置上相关实验[14,16,29]的验算，材料样品包括重金属钽、轻金属铝和以HMX 为基的塑料黏结炸药（PBX）JO-9 159，实验参数及结果列于表1。

表1 磁驱动准等熵压缩实验算例2～4 的主要参数[14,16]Table 1 Parameters of examples 2－4 for magnetically-driven isentropic compression experiments[14,16]

算例2～4 加载电流和样品表面速度的实验和计算结果示于图8。由于CQ 系列装置早期的实验样品宽度较窄，算例2（Ta）和算例4（JO-9 159）的等效电流因数均取为0.6，算例3（Al）的等效电流因数取为0.4；后续电流较高的算例中，该因数大多为0.7 以上。SSS-MHD 编码具有2 种加载电流输入方式：电路（CKT）方程组耦合计算（见图5）和实验负载电流（CUR）数据馈入，图8 将这2 种计算方式做了比较。早期工作都是使用CKT 计算的[14-16]，随着强电流测量技术的提高，本文中补充了CUR 计算，得到了更全面的认识。

图8 对于流体物理研究所CQ 装置磁驱动准等熵压缩实验的SSS-MHD 模拟（算例2～4）Fig. 8 SSS-MHD simulations for magnetically-driven quasi-isentropic compression experiments at IFP (examples 2－4)

将图8 中3 种材料样品台阶靶速度波形处理后，可获得它们的准等熵压缩线，其高端压力如表1 右端所示。算例2 中钽金属样品的密度高，同时该例的加载电流更强，使得钽的等熵线高端达80 GPa 以上。根据图8(a)可知，如果加载电流密度提高3 倍多，等熵压缩压力可提高一个量级，能得到重金属1 TPa以下范围内的等熵线，正如美国Z 装置上开展的实验[34]。当然，要把磁驱动加载电流强度提高若干倍，对驱动器硬件和实验技术都极为困难。

表1 中提及了电流计算（CUR）与电路计算（CKT）方式的比较。SSS-MHD 编码CKT 计算采用的集中参数RLC 电路可以描述很多实际驱动电路，但电路中若使用伏安特性复杂的气体开关，则RLC 电路方程就不够准确。图8(d)中CKT 计算优于CUR 方式，原因可能是该算例驱动器CQ-1.5 装置采用了伏安性能简单的爆炸开关[35]，正适合于CKT 计算的简单RLC 电路描述（早期实验使用的电流测量技术也存在较大误差）。其余2 例的驱动器CQ-4 电路采用气体开关，虽然目前简单的CKT 计算也可达到偏差低于10%的精度，但明显不如CUR 方式。本文电路耦合CKT 计算的原理是正确的，但对具体电路的描述还有待改进。

4.2 磁驱动高速金属飞片实验（算例5）

磁驱动高速金属飞片实验利用等熵压缩的磁力斜波，以极高压力（数百吉帕以至太帕）驱动作为电极板一部分的金属飞片高速射出，铝、铜飞片的速度很容易达到10 km/s 以上[21]，2011 年美国Sandia 实验室已能把亚毫米厚度的铝飞片发射到45 km/s[36]，开创了冲击压缩实验的新境界，对极端物理学研究具有重要意义。飞片运动前期由强磁压驱动，电流衰减的后期则由飞片后侧电弧烧蚀等离子体推进，这类电磁热力强耦合的复杂实验必须采取多物理场耦合计算进行设计或核算。

算例5 是Lemke 等[36]2011 年所做系列实验之“11 mm-2 s”的模拟，该系列实验研究了美国Sandia实验室Z 装置驱动高速金属飞片的能力，给定了用于直至TPa 范围的冲击压缩标准加载手段。实验“11 mm-2 s”是文献[36]中唯一报道了加载电流波形的具体实验。该实验采用单条片构形，材料为铝，电极板宽度为11 mm、长度为25 mm，构形的盲孔为矩形凹槽，槽底即矩形金属飞片，AK-gap 即磁腔间隙为1 mm。电极板凹槽顶面直接覆盖7.5 mm 厚的铜靶板，凹槽深度即飞片的飞行距离约为7.8 mm。实验测量飞片后自由面速度波形，并提供其剖面参数及高速碰撞过程的信息。计算采用的等效电流因数K为文献[36]给定值0.76（本文自算值为0.73）。图9 为SSS-MHD 编码的模拟结果及其比较。

图9 对于Sandia 实验室Z 装置铝飞片实验“11 mm-2 s”的SSS-MHD 模拟（算例5）Fig. 9 SSS-MHD simulations for the Al flyer experiment “11 mm-2 s”on the Z machine at the Sandia National Laboratories, USA (example 5)

图9(a)显示SSS-MHD 计算的“11 mm-2 s”飞片速度与原实验及ALEGRA-2D 编码模拟相当一致，但在3.2 μs 后有些偏高。图9(b)显示在磁力推动后期（3.05 μs）磁压力高达285 GPa。图9(c)表示飞行后期3.24 μs 时飞片自由面后固体区温度过高（1 000 K 左右），导致密度略低于常值，这个偏差与图9(a)后期计算速度偏高的情况相符。虽然这里没有进行飞片中熔化边界的计算，但从图9(c)中3 条剖面曲线斜率变化的拐点可以判断，飞片后自由面固态区厚度约为0.1 mm，符合Z 装置飞片实验通常的数据。图9(b)是与图9(c)时间不同的剖面，但与文献[36]中另一实验“11 mm-1 s”剖面的数据相当一致，提供了有依据的物理图像。

5 SSS-MHD 编码应用于圆筒构形电磁内爆和MC-1 发生器实验的典型算例

本节计算的金属套筒电磁内爆实验和CJ 型MC-1 发生器实验，可使材料样品中斜波压力达到太帕量级，或可使构形轴线处的压缩磁场达到700～1 000 T（相当于磁压200～400 GPa），对于低密度物质的压缩具有特殊意义。

5.1 固体金属套筒高速电磁内爆实验（算例6）

磁驱动圆筒构形的内爆实验可以高效实现样品材料的斜波压缩或套筒的高速度高对称性内爆运动，其负载构形即是作为阴、阳电极的内外两层同轴金属圆筒，其间为圆环磁腔。这种构形的主体—作内爆运动的内筒称为套筒（liner）；外筒称为回流筒（也可由若干金属柱组成），做速度不快的外向扩展运动，对内爆影响不大，然而它回流的电路功能不可或缺。由于套筒内爆运动是会聚的，其电磁驱动力和内聚压力越来越强，电流有效作用时间和效率都明显高于平面构形。金属套筒内爆速度有可能超过100 km/s，内爆压力可达到数太帕，成为高能量密度物理宏观样品实验的主要手段。

算例6 是对于美国Z 机器磁驱动Al/Cu 复合套筒超高压内爆压缩实验[34]的SSS-MHD 数值模拟。图10(a)是该实验装置的截面示意图，图中的阳极—回流筒（厚度0.45 mm、内半径13 mm）和阴极—内套筒的驱动层（厚度1 mm、外半径3.43 mm）由6061-T6 铝合金制造。外筒的膨胀速度历史由装置外围排列的VISAR 探针测量。内套筒的内层—样品层为铜材料，厚度0.53 mm、内半径1.9 mm（此系列还进行了钽和铝样品材料的内爆压缩实验）。套筒的有效高度为10 mm，内窥PDV 探针由6 根PDV 光纤组成，封装于壁厚0.15 mm、外半径0.35 mm 的铂管之中。图10(b)显示SSS-MHD 编码的计算速度与实验及ALERGRA-1D 模拟结果相符。当内爆中止即铜套筒撞上铂管时（约3.016 μs），波形优化设计的加载电流达到约16.5 MA，此时铝驱动层外半径约2 mm，电流线密度高达13 MA/cm，磁压力超过1 TPa，并且压缩过程是准等熵的。图10(c)显示铜样品中最高压力为1.25 TPa，最高密度达22.4 g/cm3，可看出SSS-MHD编码计算的剖面台阶与ALEGRA-2D 模拟结果略有差别，主要原因是介质密度接触间断的处理存在缺陷，有可能是本文计算中空间格子尺度偏大等原因造成的，需要改进。

图10 对于Sandia 实验室Z 装置上金属套筒电磁内爆实验的SSS-MHD 模拟（算例6）Fig. 10 SSS-MHD simulation of magnetically-driven liner implosions on the Z-machine at the Sandia National Laboratories, USA (example 6)

电磁内爆压缩实验的意义在于获得太帕等级的物态方程数据，关键在于如何较准确地处理实验数据。文献[34]中经过精确的物态方程计算，得到图10(c)中铜样品的最高压力、密度标定值的偏差均小于1%，是否意味着已具有状态方程参考价值，这些重要问题需要继续从实验和模拟方面进行探讨。

5.2 圆柱形炸药内爆磁通量压缩实验（算例7）

MC-1 发生器利用圆柱形炸药驱动金属套筒内爆，把套筒空腔内种子磁场的初始磁通量压缩到轴线周围小区域中，从而得到高磁场和高磁压力[37]。MC-1 发生器与5.1 节套筒电磁内爆的基本不同点在于MC-1 套筒是图3 中的圆柱形（N=1）θ 构形，而电磁内爆套筒属于z构形。MC-1 套筒的环形“励磁”电流物理上是注入套筒腔内的种子磁场Bz0在内爆套筒中感生的涡电流。从圆柱电磁构形的分析可知，z构形产生的环形磁场驱动套筒内爆具有较高效率，θ 构形压缩腔内轴向磁通量则有较高效率。

CJ-100 型MC-1 装置是中国工程物理研究院流体物理所研制的小型MC-1 发生器[38-39]，图11 是该装置的截面示意图和实物照片。CJ-100 型装置的套筒材料为304 不锈钢，外直径100 mm、厚度1.5 mm、高度210 mm。炸药为RHT-901（m(TNT)/m(RDX)=40/60），尺寸为内直径100 mm、外直径210 mm、高度65 mm，总炸药量不到3 kg。

图11 CJ-100 型MC-1 发生器的示意图和实物照片Fig. 11 Schematics and picture of the CJ-100 type MC-1 generator

算例7 是CJ-100 型MC-1 装置磁场压缩实验Shot20150630[39]的SSS-MHD 验算。该实验套筒内部的种子磁场Bz0约为5.5 T。该实验中轴线附近磁场增长历史以及套筒内、外半径变化的模拟计算示于图12(a)，在磁探针停止工作之前实测磁场与计算值十分符合。

图12 CJ-100 型MC-1 发生器磁通量压缩实验 及其性能的SSS-MHD 模拟计算（算例7）Fig. 12 SSS-MHD simulations for the magnetic flux compression experiment and the performances of the CJ-100 type MC-1 generator (example 7)

计算的套筒内层回弹半径约2.6 mm，此时计算的峰值磁场约为1 450 T。但实测磁场峰值为690 T，从图12(a)看出此时套筒腔内磁场区的半径值约3.8 mm，此值相当于实验中磁场测点的半径值。图12(b)是不同种子磁场Bz0之下SSS-MHD 计算的CJ-100 型发生器达到的峰值磁场和套筒回弹半径，可作为实验设计的参考。图12(b)中还画出了该型号发生器以前2 发实验EX1、EX2 的压缩磁场，以及对相应磁场测点半径值的估计，数据与本例相近。从图12(b)可知，设计适当的种子磁场值，可得到折中的压缩磁场值和可利用磁场的空间范围。

6 结 束 语

SSS-MHD 编码是拉格朗日形式一维平面和圆柱面几何、多介质、多组分、多连通区的弹塑性磁流体力学多物理场计算编码，是一维冲击爆轰编码SSS 向磁流体力学扩展而形成，具有25 个子程序和函数，共约6 000 条Fortran 语句。SSS-MHD 编码具有框架式结构，便于灵活增加和扩充介质类型、物性数据及物理功能。目前，除了SSS 编码原有的功能外，SSS-MHD 编码主要面向强电磁驱动的高能量密度动力学实验，特别是为极端材料动力学实验提供数值模拟平台。

本文中，第1～3 节阐述了SSS-MHD 编码的理论基础、程序结构和特色，第4～5 节给出了该编码模拟强电磁驱动实验的7 个算例，包括平面构形材料磁驱动准等熵压缩实验、磁驱动发射高速金属飞片实验、套筒电磁内爆材料压缩实验和MC-1 发生器磁通量压缩实验。总体看来，SSS-MHD 模拟与实验及其他方法数值模拟结果符合较好，相对偏差均在5%以下。这些算例的实验形式广泛多样，物理变量范围宽广，实验测量和模拟计算精度较高，为考核SSS-MHD 编码的实际能力提供了合适的试题，同时说明该编码已达到强电磁驱动实验一维数值模拟平台的要求，可以推广应用。

通过与该领域先进实验及数值模拟结果的核对，为SSS-MHD 编码的继续改进和扩展指出了方向。例如，补充流体力学模型，适应内爆动力学实验模拟的需要；在基本方程组层面引入更普遍的非线性连续介质力学本构理论，以适应大变形高应变率计算的需要；改进单温等离子体为多温多流模型，向简单辐射磁流体力学计算发展；建立更精细准确的驱动器电路和关键器件模型，实现更准确的MHD 与驱动电路的耦合计算；改进高功率激光与物质耦合计算模型，建立与SOP (streaked opticapyrometer)等先进测温技术配合的功能，进一步扩大对高能量密度物理实验的模拟能力。

本文工作得到了流体物理研究所诸多磁驱动、磁压缩实验数据的支持，谨向参加有关实验的仝延锦、唐小松、宋振飞、程诚、李建明、匡学武、吴刚、税荣杰、胥超、邓顺益、马骁等同志致以衷心感谢！