电爆炸金属桥箔早期过程中电磁-热-力多物理场耦合建模与分析

2024-02-24 09:54王刚华唐久棚欧海彬阚明先段书超

高压物理学报 2024年1期

王刚华，谢龙，肖波，王强，唐久棚，欧海彬，阚明先，段书超

（中国工程物理研究院流体物理研究所, 四川绵阳 621999）

利用大电流流经金属箔/金属丝时的欧姆加热作用，可以将金属箔/丝汽化电离，产生急速膨胀的高能量密度等离子体，进而可以作为高压驱动源。该类高压驱动源应用广泛，例如：基于该原理发展的电炮可以用于开展超高速飞片发射和材料动力学特性研究[1-4]；利用金属箔电爆炸发射小尺寸飞片来制造高安全性、低感度的冲击起爆雷管[5-10]等。

电爆炸过程中，导体被快速加热，短时间内经历固-液-气/等离子体多个状态，欧姆电阻可能增长几个数量级。因此，采用数值模拟方法研究这一过程并非易事。早期的研究主要基于电路方程和电导率模型以实现对实验电流的模拟和温度的估算[11-12]。1981 年，董玉斌等[12]将金属箔电爆炸过程分为沸点前和沸点后两个阶段，分阶段计算了放电电流和桥箔电阻的变化以及金属箔电爆炸驱动飞片的加速过程，实现了一维流体动力学方程、放电回路方程以及沙哈方程的耦合求解。1989 年，Osher 等[13]报道了用于计算金属箔电爆炸全过程的一维磁流体动力学程序，该程序与2010 年Saxena 等[14]报道的基于Burgess 电导率模型的一维磁流体程序SEEFA 一样，均对金属箔区域做零维处理。后来，罗斌强[15]、贺佳等[16]和Luo 等[17]结合大量实验完善了电箔爆炸一维磁流体程序，该程序能较好地模拟爆炸箔的放电波形，同时给出与实验基本一致的飞片速度波形。

然而，对于冲击片雷管而言，提升雷管工作的可靠性需要对桥箔电爆炸驱动飞片的过程有更为深入的认识，一维简化的磁流体模型虽然能够预估飞片的速度，但不能给出飞片的形状、温度、速度的空间分布这些可靠起爆炸药的关键信息，也无法给出爆炸箔上电流密度、磁场的空间分布，因而无法为具体的改进设计提供参考。因此，发展更全面的物理模型和更高维度的模拟程序是深入理解冲击片雷管爆炸物理过程的基础。2017 年，英国原子武器研究中心（Atomic Weapons Establishment，AWE）和美国圣地亚实验室（Sandia National Laboratories，SNL）报道了他们利用三维磁流体模拟程序ALEGRA MHD 开展的爆炸箔起爆的数值仿真工作[18]，首次获得了与实验照相结果基本相似的图像，是目前报道的最先进的冲击片雷管模拟能力。

近年来，我国的爆炸箔数值仿真工作取得了不少进展，在工程应用中发挥了较大作用[19-23]。然而，目前国内相关数值模拟工作还存在物理模型不全、空间维度降阶等问题，对于最为关键的桥箔三维结构和电磁场分布造成的等离子体空间分布等问题缺乏研究模型和研究工具。例如，北京理工大学与中国工程物理研究院化工材料研究所合作开展了三维冲击片雷管的数值模拟研究，对桥箔进行了简化处理，使用的是等效电阻模型，无法分析电磁场和电流的三维空间分布。虽然基于磁流体的数值模型在物理上较为全面，但是现有的报道还是一维的，中国工程物理研究院流体物理研究所在二维磁流体方面的研究正在进行，发展二维、三维磁流体数值模拟能力尚需时日。

尽管桥箔电爆炸的全过程模拟极其困难，但是其后续复杂的物理过程是在早期能量沉积和分布的基础上演化而来的，因此，在大量等离子体形成和爆炸膨胀之前，弄清楚桥箔上的电流分布和温度分布演化过程极其重要，有利于开展桥箔性能改进设计。本研究综合考虑电磁场在三维结构上的扩散、欧姆加热效应、材料受热膨胀效应，拟建立三维电爆炸桥箔的电磁-热-固体力学耦合模型，开展桥箔早期电磁能量沉积过程模拟。这一模型充分考虑爆炸桥箔发生剧烈汽化膨胀前的各种物理效应，是一种宏观尺度上的全三维、全物理模拟方法，对清晰理解爆炸桥箔早期过程和指导工程设计具有重要意义。

1 电磁-热-固体力学耦合模型

桥箔的早期能量沉积过程取决于电流密度和电导率，而电导率又依赖于温度场分布。由于早期温度相对较低，桥箔膨胀较小，可以使用固体力学来描述桥箔的膨胀过程。为此，建立了电磁-热-固体力学耦合的物理模型以模拟桥箔早期的能量沉积过程。采用电磁模块进行桥箔的电磁和电流演化分析，采用固体传热模块考虑桥箔的温度变化，采用固体力学模块模拟桥箔的早期受热膨胀情况。

电磁模块求解如下控制方程。

对于桥箔区域，使用A形式的矢量势公式（安培定律）描述

对于空气区域，使用标量势公式（磁通量守恒）描述

式中：H为磁场强度，A/m；B为磁通密度，T；E为电场，V/m；A为磁矢势，Wb/m； σ为电导率，S/m； µ0为真空磁导率，H/m； µr为相对磁导率；Vm为磁标势，A。

固体传热模块求解如下控制方程

式中： ρ为密度，kg/m3；cp为比定压热容，J/(kg·K)；T为温度，K；u为速度，m/s；k为热导率，W/(m·K)；Qe为热源，W/m3，耦合计算中热源为欧姆加热，即Qe=J·E。

固体力学模块求解如下控制方程

式中：s为位移，m；S为应力张量，Pa；Fv为外力，N。

耦合计算中调用热膨胀模块，即考虑由热效应引起的热应变

式中： α(T) 为随温度变化的热膨胀系数，K-1； εth为弹性应变；C为弹性模量张量，Pa；Tref为参考温度，K。

计算初始时刻，电场、磁场强度、电流密度等均为零，温度为常温。电磁边界条件主要包括：对称面上指定磁标势V0，区域边界上法相磁场为零，空气与桥箔交界面上切向磁场连续。固体传热模块边界条件：将所有面均设定为绝热条件。力学边界条件主要包括：对称面上位移为零，其余面为自由面。计算中将图1 所示的实验电流作为激励电流，对应的回路电阻为280 mΩ，回路电感为310 nH，电容为0.25 μF，充电电压为5.1 kV。桥箔材料为铜，考虑了铜材料的热导率（k）随温度（T）和密度（ ρ）的变化以及电导率随温度和密度的变化情况（如图2 和图3 所示）。

图1 实验获得的电流历史曲线Fig.1 Experimental current history curve

图2 铜的热导率Fig.2 Thermal conductivity of copper

图3 铜的电导率Fig.3 Conductivity of copper

具体求解过程中，3 个模块的耦合关系如下：电磁求解模块通过电流密度J将电流分布传递给热传导模块，传热模块通过温度T与固体力学模块进行耦合；同时，材料温度的改变通过电导率模型传递给电磁模块，用于计算实时的电磁场分布；固体力学计算产生的几何变形，也会耦合传递到电磁模块和传热模块。

为了减少计算量，采用如图4 所示的1/2 几何模型进行分析，并通过镜像还原整体情况（图4 中未显示外围空气区）。由于物理模型只适用于小变形，因此计算时使用拉氏网格，网格划分采用自适应四面体网格，在所关心的桥箔区采用密网格，空气区使用较粗网格，总网格数约为2.084×107，如图5 所示。图4 中，由于桥箔的厚度（z方向）很薄，仅为5 μm（为了便于显示，图4 中将厚度放大了10 倍），远小于x、y方向上的尺寸（毫米量级），为了减小计算量，划分网格时，铜箔区在z方向只划分了2 层网格。为便于后续沿标记线进行物理量分布分析，在如图4 所示的桥面上标记两条线L1、L2，其中：L1 为桥面上对称轴线，L2 为边缘轮廓线，弧长为l。

图4 爆炸箔几何1/2 模型示意图（单位：mm）Fig.4 Schematic diagram of the 1/2 model of explosion foil (Unit: mm)

图5 计算使用的网格Fig.5 Grid diagram used for calculation

2 计算分析

2.1 基于常数电导率的简要分析

采用常温电导率模型分析金属桥箔的磁场和电流密度分布，不考虑材料电导率、热导率随温度和密度的变化。由于常温电导率模型只能用于分析放电早期的电磁场，为了更接近实际情况，需要进一步考虑欧姆加热的影响，因此，只对常数电导率模型的计算结果作简要介绍。

鉴于材料性质不发生变化，可采用归一化电流进行分析。图6 给出了桥箔区附近的磁场强度分布。受几何构型的影响，磁场强度在桥区和过渡区的交界处表现出明显的非均匀性，特别是边缘附近非均匀性更为显著。因此，边缘附近的电流密度也很大，如图7 所示。

图6 表面磁场强度分布Fig.6 Magnetic field intensity distribution of upper and lower surfaces

图7 电流密度与方向分布Fig.7 Current density and directional distribution

依据图7 的电流密度分布，绘制电流密度的高度图，如图8 所示。图9 给出了沿L1、L2 的电流密度随弧长l的分布，弧长起点为桥箔左端。从图8 和图9 可以非常直观地看出：随着过渡区变窄，电流密度逐步增加，在桥区和过渡区交界处的4 个角点上电流最为集中；在桥区主体部分虽然还存在一定的边界效应，但桥面主体部分的电流分布较为均匀。

图8 电流密度高度Fig.8 Current density height

图9 线上电流密度分布Fig.9 Current density distribution on lines

爆炸箔上的电流分布决定了各处产生的电功率密度Pe（图10），也随即决定了各处的温升和早期等离子体的分布，因为在电功率高的地方，材料被率先加热，进而汽化形成等离子体。

图10 电功率密度高度Fig.10 Electric power density height

本节的所有分析是基于常数电导率模型开展的，即爆炸箔的电导率在空间上均匀分布，因此，电流和磁场只与爆炸箔的几何尺寸相关，各个时刻获得的磁场强度、电流密度的空间分布类似，不同时刻电流的大小仅影响其绝对值。然而，实际上，金属的电导率与其自身的温度和密度密切相关，需要考虑这种因素的影响。

2.2 基于Burgess 电导率的计算分析

常数电导率模型只能分析放电早期的电磁场，在中后期，为了更接近实际情况，需要进一步考虑欧姆加热的影响。事实上，随着电流加热形成等离子体，材料的电参数会发生复杂的变化，为此，在电导率模型中引入随温度变化的Burgess 电导率模型，并引入固体传热模块和固体力学模块，以分析桥箔放电一段时间内的电磁场分布情况。

图11 给出了引入Burgess 电导率模型、固体传热模块和固体力学模块后计算得到的不同时刻爆炸箔上的磁场强度分布。可以看出，与常数电导率模型类似，改进后的模型计算结果显示，磁场先集中分布在爆炸箔4 个角点和爆炸箔边缘，这些地方的磁场梯度较大，磁场分布极不均匀。在桥区，由于宽度变窄，磁场相对较大，边缘效应更加明显，特别是桥区与过渡区衔接的4 个角点的磁场强度很大。与常数电导率模型计算结果不同的是，随着时间的增长，磁场向爆炸箔中间扩散，中间区域的磁场逐步增大。磁场扩散呈现两种扩散特征：一是从4 个角点出发的点扩散，二是从边缘向里的线扩散。

图11 不同时刻的磁场强度分布Fig.11 Magnetic field intensity distribution at different instants

电流分布主要由两个因素决定：构型和电导率。前述的磁场角点和边缘效应是构型引起的。由于电流密度对应磁场强度的旋度，因此，磁场不均匀的地方电流密度更大。图12 给出了由改进后的模型计算获得的不同时刻的电流密度分布。可以看出：电流一开始集中在边缘，然后向中间扩散，此时构型是电流分布的主导因素，在70 ns 时，桥区主体的电流分布基本均匀；此后，电导率因素占据主导，边缘由于温升快，其电导率低于中心处的电导率，使得电流向中心集中；后续电流分布在两种因素的竞争中不断演化。图13 给出了90 ns 时刻沿线L1、L2 的电流密度分布，可以看出，其大致特征与常数电导模型率给出的结果一致，只是桥区的电流分布更均匀。

图12 不同时刻的电流密度分布Fig.12 Current density distribution at different instants

图13 90 ns 时刻沿线的电流密度分布Fig.13 Current density distribution on lines at 90 ns

爆炸箔的主要温升机制是欧姆加热，温度分布取决于电流密度、电导率和时间。图14 给出了温度分布的演化。可以看出，早期边缘加热效应明显，4 个角点的温升相比于其他地方更快，温度也更高。但是，随着磁场和电流的扩散作用，电流向温度较低的中间区域集中，桥区主体部分的温度分布逐渐变得匀滑。图15 给出了沿线L1、L2 的温度分布。两条曲线的平台区对应桥区，两条曲线的平台有错位，这是L1、L2 的弧长不一致导致的。两条曲线的平台有高度差，反映了桥区边缘与桥区中心有温度差。92 ns 时，桥区大部分区域的温度均在1 245～1 345 K之间，温度差在100 K 以内，这意味着在该结构设计下爆炸箔的温度分布已十分均匀。通常情况下，在设计冲击片雷管时，将桥箔爆炸时刻设计在电流峰值附近，虽然计算时刻与电流峰值时刻相距时间较长，尚不能确定后续温度分布是否更加均匀，但是从目前的计算结果来看，温度均匀性已经较好。此外，爆炸箔紧邻桥区的过渡区中有相当多的区域温升较为明显，意味着该处沉积了较多能量。在冲击片雷管设计中，可以充分利用该处产生的等离子体能量来驱动飞片，使飞片获得更多动能。因此，从提高能量利用率的角度来讲，设计冲击片雷管时可考虑使炮筒口径略大于桥区。

图14 不同时刻的温度分布Fig.14 Temperature distribution at different times

图15 92 ns 时刻线上的温度分布Fig.15 Temperature distribution on lines at 92 ns

图16 给出了90 ns 时爆炸箔沿z轴的膨胀速度分布，4 个角点附近的膨胀速度相对较快，桥区边缘的膨胀速度比桥区中心慢，其原因是边缘处金属箔除了z向膨胀以外，还向宽度方向膨胀。

图16 90 ns 时刻沿z 轴的速度分布Fig.16 Velocity distribution alone z-axis at 90 ns

图17 显示了采用改进模型计算的温度分布与文献[16]报道结果的对比，其中：右图为3 个时刻的温度分布灰度图，左图为文献[16]中ALEGRA MHD3D 的模拟结果，中间图为文献[16]给出的实验照片。可以看出，文献[16]的模拟结果没有抓拍到实验观测到的边缘加热效应，而本研究的模

图17 模拟获得的温度分布与文献[16]的对比Fig.17 Comparison between simulated temperature distribution and Ref.[16]

拟结果较为清晰地显示了前期的边缘加热效应和持续时间较长的角加热效应，4 个角点附近的扩散现象也与实验结果定性相符。

以上分析只适应于电爆炸箔的初步设计及其工作过程的早期情况。在角点出现明显的受热膨胀以后，建立的模型受制于固体力学模型所能描述的范围，无法进行后续过程的模拟。事实上，随着电流加热形成等离子体，材料的电参数会发生复杂的变化，其动力学行为将是复杂的磁流体力学过程，需要开展更为深入的磁流体力学建模研究。

3 结论

建立了三维电磁-热-固体力学耦合求解模型，模拟了爆炸箔在大电流加载下的早期受热膨胀过程。基于常数电导率模型，给出了爆炸箔的磁场、电流密度、电功率密度的空间分布。基于Burgess 电导率模型，分析了早期膨胀过程中桥箔的磁场、电流密度、温度的演化情况，指出影响其分布的决定因素为构型和桥箔电导率分布。模拟了爆炸箔的膨胀速度，计算结果较好地捕捉到了实验中观测到的边缘加热效应和角加热效应。