一种耦合电路分析的磁驱动飞片数值计算方法＊

张旭平，赵剑衡，谭福利，王桂吉，罗斌强，莫建军，种 涛，孙承纬，刘仓理

（1.中国工程物理研究院流体物理研究所，四川 绵阳621999；2.中国工程物理研究院，四川 绵阳621999）

磁驱动加载技术是近十年发展起来的一种新的脉冲载荷加载技术，在高能量密度物理、冲击动力学和航空航天等领域具有重要应用［1－3］。对磁驱动飞片的数值模拟研究主要是在对飞片加载历史计算［1，4］、飞片击靶前状态确定［5］、负载优化设计［6］等方面。磁驱动过程中极板受洛伦兹力作用会产生严重的变形，并且极板加载面受焦耳热而有很大的温升，导致极板加载面有烧蚀、相变［7－8］。这些负载变化都会反馈到回路中影响磁驱动装置对负载的放电电流。如果磁流体计算软件不能耦合电路计算，程序的计算范围只能是利用实验数据进行计算和物理分析，则一般的磁流体动力学模拟计算只有在给定电流是实验电流的情况下，得到的结果才接近于真实情况，这样对一些预测性的计算分析将非常困难。并且研究磁驱动装置负载参数对电路的影响和对实验放电物理过程的认识都需要采用耦合电路的磁流体动力学计算［9－10］。相对于国外磁流体动力学软件的快速研发，国内数值模拟方面工作还处在起步阶段。国内的模拟研究都是在磁流体动力学软件中直接导入电流曲线进行计算的［11－12］。

鉴于动态电路分析的重要性，本文中利用LS－DYNA980 MHD计算软件和自编的电路分析程序，通过2个计算程序计算结果的迭代，建立一种等效耦合电路的磁流体动力学计算方法，并且利用建立好的计算方法对实验放电电流曲线进行模拟，并分析磁驱动相同尺寸飞片时极板长度的选择。

1 耦合电路计算方法和流程

1.1 电路计算模块建立和参数确定

对磁驱动实验平台电路分析时可以简化为一个R－L－C电路［13］。脉冲功率加载装置CQ－4的等效电路如图1所示，第1部分为主电容、开关及第1段传输线的总等效电容C1、电感L1和电阻R1，第2部分为峰值电容部分的等效电容C2、电感L2和电阻R2，第3部分为负载区的等效电阻R3和电感L3。

根据电容、开关等的额定参数和磁驱动实验，可以确定R1、L1、C1、R2、L2、C2的值。R1中包含气体开关的火花电阻r1，r1在放电过程中不断变化，在放电初期r1较高，在放电回路电流很大时则很低。

根据等效的电路我们编了电路求解程序，程序主要采用四阶龙格－库塔方法求解微分方程组：

图1 脉冲加载装置电路图Fig.1 Circuit of pulsed power generator

这个程序的优点是电路中的各参数可以赋固定的值，也可以赋值为时变的表达式。

利用程序和拟定的R、L、C参数对CQ－4短路放电实验 CQ－4 shot－1和 CQ－4 shot－2做了模拟计算，结果如图2所示，分别是在CQ－4充电电压73、80 k V下的短路电流曲线。程序计算的电流曲线和实验值符合较好。

图2 模拟和实验电流曲线对比Fig.2 Comparison of simulated and measured load currents

1.2 计算流程和方法

磁驱动飞片的计算是一个典型的耦合电路的磁流体动力学计算问题。计算中以磁驱动装置的充电电压和负载区参数（极板构型、飞片大小）作为初始条件，得到装置的放电电流和飞片速度等结果。一般耦合电路的磁流体动力学计算都是按时间步计算的（见图3（a））。以装置的充电电压和静态电路参数为初始条件，通过程序的电路模块计算出加载到负载的电流值，利用程序的磁流体动力学模块计算得到该时间步负载的电磁－力－热响应结果和新的负载电感、电阻值。接着将得到的新的负载参数反馈到下一个时间步的计算，作为下一个时间步的输入条件，依次循环直到计算结束。图中V是充电电压，T为计算时间，R、L是负载电阻和电感，Lt、Rt、it分别表示t时刻的电感、电阻和电流值；MHD即电磁－力－热耦合的磁流体动力学计算，加入材料的本构关系、状态方程、热状态方程、电导率模型等，通过求解流体动力学方程组、麦克斯韦方程组和热传导方程，得到负载的结构变形、受力、温度、电磁场分布等结果；CIRCUIT指电路分析计算，已知t时刻的参数Lt、Rt等，求解电路方程组（式（1）～（4）），得到电流值it。

由于建立以上耦合计算方法的磁流体动力学程序很困难，本文中在磁流体动力学程序LS－DYNA980 MHD和自建的电路计算程序的基础上，利用实验时各参数（电流、电感L（t）、电阻R（t）等）对时间t的单值性，通过2个程序计算结果的迭代实现了耦合电路的磁流体动力学计算。单值指实验中任一时刻的所有参数（电流、电感、电阻等参数）都只有唯一值。用于数值计算即用t＋Δt时刻计算输入的电流值I（t＋Δt）和用t时刻的整个电路的参数（电感L（t）、电阻R（t）等）求解的电流值是相等的，具体计算流程如图3（b）所示。图中t为计算时间，T为总计算时间，Δt为时间步长，n为循环次数，N 为最大循环次数，MHD是磁流体动力学程序LS－DYNA980 MHD，CIRCUIT为自编的外部电路计算程序，Iin是输入磁流体计算软件的电流曲线，Iout是CIRCUIT程序计算出的电流，In（t＋Δt）是第n次迭代电流曲线的（t＋Δt）时刻的电流值。

图3 计算流程示意图Fig.3 Schematic diagram of calculation flow

本文的计算流程（见图3（b））与以上计算（见图3（a））的主要区别是：第1，初始输入条件是一条电流曲线；第2，不能在计算中将磁流体动力学计算的负载电路参数变化及时反馈到下一时间步的计算，而是磁流体和电路模块分别计算后迭代。计算中以｜In（t＋Δt）out－In（t）in｜≤η作为单值原理的判据，表示两曲线上对应时刻的电流值相差小于η为止。否则下一循环输入的电流曲线为

这种计算的收敛速度很快，一般循环次数N≤5时，η就已经很小，能得到很好的结果。

2 负载动态参数计算

2.1 负载动态电感计算

耦合电路的磁流体动力学计算的难点是对负载时变电感和电阻的计算，磁驱动是一个电磁－力－热耦合问题，计算中涉及到极板变形、磁扩散和材料相变等问题。电流、电磁场的分布很难准确计算，对电感和电阻计算也有一定的影响。

磁驱动飞片实验负载的原理示意图如图4所示，电流沿极板的内表面传播，在极板之间形成磁场，

流过极板的电流与极板加载 面磁场相互作用后在两极 板上产生相互排斥的磁压 力。图4中h为极板到短路端的长度，d为极板 之间的距离，a为从极板 电流入口处到飞片处的长 度，b为从飞片到极板短路处距离。在飞片发射过程 中两极板之间间隙增大和 极板变形都会影响磁场的 分布和大小，导致负载电感变化。

图4 磁驱动飞片负载物理模型Fig.4 Physics model of magnetically driven flyer

F.G.Steven等［9］在计算电感时采用了2种方法。第1种是加载的时变电感与加载时变电流有如下的关系：

式中：μ为磁导率，I（t）为电流，B（t）为磁场。利用2D－MHD对计算区域内磁场能量积分的方法计算时变电感L（t）。第2种是等效计算一个完全电导体（PEC）条片模型的电场，PEC条片模型中条片之间中点的磁场强度是MHD计算的不同时间点的磁场强度，通过计算模型的电场能量计算电感：

式中：E/V是归一化电场。以上2种方法电流为零时计算都存在奇点。

采用假定电流在极板表面均匀分布，对一端短路的平行条片的电感计算公式［14］：

式中：W 为极板宽度，h为极板长度，d为极板之间的距离。而其中公式极板之间的时变距离则通过三维的LS－DYNA980 MHD计算软件计算。电流均匀分布于极板表面的假定，忽略了磁扩散的影响和极板宽度方向电流的不均匀性，计算中需要对极板间隙d（t）做修正。

2.2 负载动态电阻计算

电阻计算采用下式：

式中：电流流过的横截面积高度δ的计算采用磁扩散深度代替，假设在这一磁扩散深度内阻值均匀。

从欧姆定律和麦克斯韦方程导出导体中的磁通密度的演化方程［14］：

假设只有宽度方向的磁场且电导率是固定值，宽度方向的磁场表示为

向电极板厚度方向磁扩散深度和扩散时间的关系也可表示为：

式中：σ为电导率，μ为磁导率。

将式（12）代入电阻计算公式（9）推出

实验中极板加载面受焦耳热影响温度升高，甚至极板受焦耳热烧蚀发生相变，导致材料电导率在加载中发生变化，推导前面的公式是假定电导率不变。模拟中首先取2个电导率的区间节点计算了时变电阻（见图5），一个值是电导率初值（曲线1），另一个是MHD计算时间终点时温度和密度变化后的加载面电导率值（曲线2）。这2个电导率的值相当于实验中的最大值和最小值，实验中负载的电阻值应该在图中计算的2条时变电阻的区间内。而对于加载初期电导率接近于初值，对于加载后期电导率接近于下降后的值。计算的极板中电流流过的等效厚度内，能量与温度Tm成正比，能量与时间也成正比，温度与时间成正比，电导率随温度的变化曲线可以拟合为一个衰减函数式

即在前面假设的条件下可以推出电导率随时间的变化曲线，将其代入电阻计算公式（13），在前面的假定下磁驱动飞片过程电阻按图中曲线3变化。电阻计算假定宽度方向电流均匀和在计算的磁扩散厚度内各处电导率保持均匀一致。

在前面电感和电阻的计算假定下，在磁扩散较小和极板电流分布均匀的情况下能给出较好的计算结果。我们采用建立的耦合电路的磁驱动飞片的计算方法，对CQ－4的驱动飞片实验shot－13和shot－32的电流结果进行了数值计算（见图6），实验CQ－4 shot－13中飞片尺寸为10 mm×5 mm×0.7 mm，实验CQ－4 shot－32中飞片尺寸为Ø5 mm×0.35 mm。从图6可以看出，计算结果和实验结果在飞片能测到速度的时间范围内符合较好。后期飞片脱离极板在镗孔中飞行和极板加载面受焦耳热烧蚀后形成等离子体层等，这些原因均可导致实验中电流不再沿飞片加载面传播，而是直接击穿在等离子体中传播，这时计算的电感和电阻值都出现偏差。

图5 负载电阻变化曲线Fig.5 Resistance varied with time

图6 驱动飞片实验电流和计算电流比较Fig.6 Measured and calculated currents of flyer plate driven experiments

3 负载极板长度设计

磁驱动飞片极板长度的设计主要受2个因素影响，第1是飞片沿极板长度方向的平面性，第2是回路电感。实验中电流在极板流入端有一个汇流过程，电流密度由小变大；在负载的短路端，电流流过垂直极板的短路板，形成的磁场会叠加到极板之间的磁场；这2种原因都导致极板加载面从电流入口端到短路端磁场逐渐增强，在发射飞片时会造成飞片加载面压力不均匀，从而使飞片平面性变差。通过增加极板长度可以减弱这种不均匀性，但是极板长度增加会造成回路电感增大，导致加载电流减小，影响飞片速度。对只给定加载电流曲线计算极板长度影响加载面磁压力均匀性的模拟计算，可以给出均匀性结果，但对飞片速度或加载电流的大小不能比较，不能起到优化计算的效果。耦合电路的磁流体动力学计算不仅可以给出飞片平面性信息，也可以给出飞片速度历史。

下面分别对相同飞片长度（l＝10 mm）、极板宽度（w＝6 mm），不同极板长度情况下飞片的平面性和速度进行数值计算。第1种情况是极板长度h＝20 mm、a＝5 mm、b＝5 mm；第2种情况h＝25 mm、a＝10 mm、b＝5 mm；第3种情况h＝30 mm、a＝15 mm、b＝5 mm。与第1种情况相比后2种情况只是增加了入口端到飞片处的极板长度a（见图4）。

图7是3种极板长度在飞片发射实验充电电压73 k V时发射飞片的自由面速度比较，图8是飞片长度方向的平面性比较。由图7～8可知，在加载时间t＝0.72μs时，极板长度h＝25 mm的飞片比h＝20 mm的飞片速度低约1.2 km/s，但是飞片的平面性有明显的提高。极板长度h增加到30 mm时，相比于h＝25 mm时对飞片的平面性提高较小。所以认为实验中极板长度取h＝25 mm比较合适。

图7 不同极板长度h时飞片自由面速度比较Fig.7 Comparison of free surface velocities under different strip line lengths

图8 极板长度h不同时沿极板长度方向飞片的平面性比较Fig.8 Comparison of the flatness under different strip line lengths

4 结 论

分析了磁驱动飞片实验和数值计算的研究现状。根据实验装置等效电路建立了可带入参数为表达式的电路计算程序。利用建立的电路计算程序结合LS－DYNA980 MHD计算软件，建立了一种可以耦合电路分析的磁驱动飞片计算方法。采用这种方法模拟了CQ－4的发射飞片实验，计算结果与实验基本符合。最后计算了极板长度h分别为20、25、30 mm在充电73 k V时磁驱动飞片的速度和平面性，认为实验中极板长度取25 mm比较合适。

磁驱动飞片涉及到高压、高温、等离子体烧蚀等问题，数值计算中普遍存在的问题是材料本构关系、状态方程、电导率模型的选择对计算结果的影响较大。文中对电感和电阻的计算由于涉及到磁扩散和焦耳热烧蚀，耦合电路的磁驱动计算结果的误差也主要来自于电阻和电感的计算的误差。计算流程为耦合电路的磁流体动力学计算模拟提供了一种新的耦合计算思路。

本文工作得到陈学秒、吴刚、蔡进涛、税荣杰、胥超、马骁、邓顺益的帮助，在此特别感谢！

［1］Lemke R E，Knudson M D，Davis J P.Magnetically driven hyper－velocity launch capability at the Sandia Z accelerator［J］.International Journal of Impact Engineering，2011，38（6）：480－485.

［2］孙承纬，赵剑衡，王桂吉，等.磁驱动准等熵平面压缩和超高速飞片发射实验技术原理、装置及应用［J］.力学进展，2012，42（2）：206－218.Sun Cheng－wei，Zhao Jian－heng，Wang Gui－ji，et al.Progress in magnetic loading techniques for isentropic compression experiments and ultra－high velocity flyer launching［J］.Advances in Mechanics，2012，42（2）：206－218.

［3］赵剑衡，孙承纬，谭福利，等.一维平面磁驱动等熵加载发射飞片技术［J］.爆炸与冲击，2005，25（4）：303－308.Zhao Jian－heng，Sun Cheng－wei，Tan Fu－li，et al.Launch technique for isentropic compression flyer plates magnetically driven by using fast pulsed power［J］.Explosion and Shock Waves，2005，25（4）：303－308.

［4］Lemke R W，Knudson M D，Robinson A C，et al.Self－consistent，two－dimensional，magneto－hydrodynamic simulations of magnetically driven flyer plates［J］.Physics of Plasmas，2003，10（5）：1867－1874.

［5］Lemke R W，Knudson M D，Bliss D E，et al.Magnetically accelerated，ultrahigh velocity flyer plates for shock wave experiments［J］.Journal of Applied Physics，2005，98（7）：073730.

［6］Ao T，Asay J R，Chantrenne S，et al.A compact strip－line pulsed power generator for isentropic compression experiments［J］.Review of Scientific Instruments，2008，79（1）：013903.

［7］Gael L B，Petit J，Chanal P Y，et al.Modelling the dynamic magneto－thermomechanical behaviour of materials using a multi－phase EOS［C］∥7th European LS－DYNA Conference.Salzburg，Austrian，2009.

［8］Gael L B，Chanal P Y，Herei P L.Ramp wave compresion in a copper strip line：Comprarison between MHD numerical simulations（LS－DYNA）and experimental results（GEPI device）［C］∥10th International LS－DYNA users conference.Dearborn，MI，USA，2008.

［9］Glover S F，Schneider L X，Reed K W，et al.Genesis：A 5 MA programmable pulsed power driver for isentropic compression experiments［C］∥IEEE Pulsed Power Conference（PPC）.Washington DC，2009：763－767.

［10］Glover S F，Davis J P，Schneider L X，et al.Impact of time－varying loads on the programmable pulsed power driver called genesis［J］.IEEE Transacts on Plasma Science，2012，40（10）：2588－2596.

［11］王刚华，孙承纬，赵剑衡，等.磁驱动平面飞片的一维磁流体力学计算［J］.爆炸与冲击，2008，28（3）：261－264.Wang Gang－hua，Sun Cheng－wei，Zhao Jian－heng，et al.One－demensional，magnetohydrodynamic simulations of magnetically driven flyer plates［J］.Explosion and Shock Waves，2008，28（3）：261－264.

［12］王桂吉，蒋吉昊，孙承纬，等.磁驱动飞片的一维磁流体动力学数值研究［J］.计算力学学报，2008，25（6）：776－781.Wang Gui－ji，Jiang Ji－hao，Sun Cheng－wei，et al.One－demensional，magneto－hydrodynamic simulation of magnetically driven flyer plates［J］.Chinese Journal of Computational Mechanics，2008，25（6）：776－781.

［13］王桂吉，谭福利，王刚华，等.小型脉冲功率发生器放电电流波形的调节模式［J］.高电压技术，2007，33（11）：223－226.Wang Gui－ji，Tan Fu－li，Wang Gang－hua，et al.Shaping models of discharging current pulse for the compact pulsed power generator［J］.High Voltage Engineering，2007，33（11）：223－226.

［14］孙承纬.磁驱动等熵压缩和高速飞片的实验技术［J］.爆轰波与冲击波，2005（3）：125－138.Sun Cheng－wei.The technique of magnetically driven isentropic compression and high－velocity flyer plates［J］.Detonation and Shock Waves，2005（3）：125－138.