磁驱动加载技术是近十年发展起来的一种新的脉冲载荷加载技术, 在高能量密度物理、冲击动力学和航空航天等领域具有重要应用^[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段传输线的总等效电容C₁、电感L₁和电阻R₁, 第2部分为峰值电容部分的等效电容C₂、电感L₂和电阻R₂, 第3部分为负载区的等效电阻R₃和电感L₃。

图  1  脉冲加载装置电路图

Figure  1.  Circuit of pulsed power generator

下载: 全尺寸图片 幻灯片

根据电容、开关等的额定参数和磁驱动实验, 可以确定R₁、L₁、C₁、R₂、L₂、C₂的值。R₁中包含气体开关的火花电阻r₁, r₁在放电过程中不断变化, 在放电初期r₁较高, 在放电回路电流很大时则很低。

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

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

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

图  2  模拟和实验电流曲线对比

Figure  2.  Comparison of simulated and measured load currents

下载: 全尺寸图片 幻灯片

1.2   计算流程和方法

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

图  3  计算流程示意图

Figure  3.  Schematic diagram of calculation flow

下载: 全尺寸图片 幻灯片

由于建立以上耦合计算方法的磁流体动力学程序很困难, 本文中在磁流体动力学程序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为自编的外部电路计算程序, I_in是输入磁流体计算软件的电流曲线, I_out是CIRCUIT程序计算出的电流, I_n(t+Δt)是第n次迭代电流曲线的(t+Δt)时刻的电流值。

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

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

2.   负载动态参数计算

2.1   负载动态电感计算

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

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

图  4  磁驱动飞片负载物理模型

Figure  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条时变电阻的区间内。而对于加载初期电导率接近于初值, 对于加载后期电导率接近于下降后的值。计算的极板中电流流过的等效厚度内, 能量与温度T_m成正比, 能量与时间也成正比, 温度与时间成正比, 电导率随温度的变化曲线可以拟合为一个衰减函数式

图  5  负载电阻变化曲线

Figure  5.  Resistance varied with time

下载: 全尺寸图片 幻灯片

即在前面假设的条件下可以推出电导率随时间的变化曲线, 将其代入电阻计算公式(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可以看出, 计算结果和实验结果在飞片能测到速度的时间范围内符合较好。后期飞片脱离极板在镗孔中飞行和极板加载面受焦耳热烧蚀后形成等离子体层等, 这些原因均可导致实验中电流不再沿飞片加载面传播, 而是直接击穿在等离子体中传播, 这时计算的电感和电阻值都出现偏差。

图  6  驱动飞片实验电流和计算电流比较

Figure  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时飞片自由面速度比较

Figure  7.  Comparison of free surface velocities under different strip line lengths

下载: 全尺寸图片 幻灯片

图  8  极板长度h不同时沿极板长度方向飞片的平面性比较

Figure  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比较合适。

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