1 引言
炸药爆炸伴随着电磁辐射现象的产生，这种现象不仅会对爆炸场的信号测量产生干扰，严重时还会对测量系统和设备产生破坏。戴晴等[1]用150 g铝镁工质进行实验，其爆炸产生的电磁辐射使5 m外的计算机死机。常规炸药爆炸产生电磁信号这一现象最早发现于1954年，由Kolsky发表于Nature[2]，之后不断有学者进行相关的实验研究。1958年，Cook[3]针对常规炸药爆炸产生电磁波这一现象，提出了一种解释机理，他认为是导电的爆轰产物在地球背景电磁场中产生了极化，运动至接地后放电产生电磁波。Boronin等[4-5]进行了针对性实验，来验证这个解释。在之后仍然有大量的研究，主要集中在爆炸产生的电磁信号的测量和丰富实验成果上。2001年Fine[6]提出了更为直观的理论解释，他认为爆炸产生的等离子体在激波前沿冲击波厚度内被加速引起了电磁辐射，给出了详细的计算公式，可以用来估算爆炸产生电磁场强度的量级，并给出了大量的实验数据。2014年Kuhl[7]提出了一种新的理论模型，Kuhl等认为等离子体在背景磁场中的集群运动会产生电流，而这个电流会导致周围空间电磁场的变化，他们根据这一理论模型建立了一维爆炸问题模型，模拟了球形爆轰产生的电磁辐射现象。
国内也在这方面进行了一系列研究，但相比于国外的研究，国内的研究起步较晚，研究也主要集中在实验方面。1997年，陈生玉等[8]测量了带壳装药爆炸引起的电磁辐射，用量纲分析法对实验测得的数据进行分析，得出了电磁脉冲的幅值随装药量增大而增大等结论，分析了电场强度和各个爆炸参数之间的关系。2011年，曹景阳等[9]对聚能炸药索爆炸引起的电磁干扰进行了测量，分析实验结果得到爆炸索爆炸后的数十毫秒内都会有明显的电磁脉冲出现，其频率主要集中在兆赫兹量级，在单个频点上的电场峰值可达数个V/m。王长利等[10]对B炸药和梯黑铝炸药爆炸产生的电磁效应进行测量，结果表明炸药爆炸产生的电磁辐射信号强度与炸药当量的1/3次方成正比，随距离增加而衰减。几种典型炸药的爆炸电磁辐射信号主要频率在100 MHz以内，同一种炸药信号频谱特征与当量、爆心距无关，不同种类炸药产生信号频谱不同。陈鸿等[11]对RDX基含铝炸药爆炸电磁辐射信号特性进行了实验研究。结果表明，RDX以及RDX基含铝炸药爆炸电磁辐射信号发生时刻与起爆时刻相比有明显延迟。距爆心2 m处，爆炸电磁辐射信号强度在1.87～15.20 V/m范围内，随距离的增加而衰减。目前国内外对炸药爆炸产生电磁辐射的认识仍然不够清晰，对炸药爆炸过程中电磁辐射的产生机理和影响因素缺乏全面的了解。因此本研究能加深对炸药爆炸过程中产生电磁辐射问题的认识，揭示该过程中电磁辐射产生机理与时空演化规律，并在爆炸测试电磁防护方面有一定的工程应用价值。
由于国内外在炸药爆炸产生电磁辐射领域内的数值模拟研究很少，因此对其进行数值模拟研究是十分必要的。本文将爆炸流体计算方法、热平衡电离方程、非理想磁流体动力学模型和麦克斯韦方程相结合，通过基于AUSM+-up通量分裂方法的有限体积法进行炸药爆炸产生电磁辐射数值模拟，并针对不同当量的装药进行了数值模拟研究，给出了当量对炸药爆炸产生电磁辐射的影响规律。
2 理论模型
2.1 爆炸的电离模型
从Lee[12]和栗建桥[13]等的研究可知，炸药爆炸过程中产生的等离子体处于热平衡状态，因此选择热平衡电离模型作为炸药爆炸产生等离子体的模型。在考虑了体系势能、温度和电离程度后，文献[12]给出了炸药爆炸产生等离子体的电离模型为：
(1)
式中：ne为电子的数密度；me为电子质量；k为玻尔兹曼常数；T为气体温度；h为普朗克常数；p为各个组分；Ip为组分p的一阶电离能。
可以得出电导率计算公式为:
(2)
式中：q为原电荷数；d为中性粒子的有效直径；n为电离气体所有粒子的数密度。
2.2 MHD模型
相关文献[7]分析认为，炸药爆炸产生的电磁辐射是由高温高压的导电爆轰产物在自然磁场中运动而产生的，因此采用磁流体动力学(MHD)模型来描述等离子体的运动和电磁场的演化，以实现对爆炸过程中电磁辐射的数值模拟。结合磁扩散后的MHD控制方程如式(3)所示。
(3)
式中： ρ为流体密度；v为速度；B为磁感应强度；e为比内能；σ为电导率；ε为介电常数；μ为磁导率；p为压强。从式(3)可以看出，电导率σ对最后一个方程有很大的影响。当电导率很高时，磁感应强度的时间和空间二阶导数都很小，磁场完全冻结在等离子体中，并随着等离子体的运动而运动，这对应着理想情况下的磁流体动力学模型。当等离子体的电导率和介电常数相当时，不仅要考虑到磁场随等离子体运动产生的变化，还需要考虑到磁扩散效应，电磁波效应会比较显著。
从上述分析可以看出，炸药爆炸过程中，爆轰产物电离产生的等离子体在地磁场中运动是爆炸过程中电磁辐射的根源。而等离子体的电导率是影响爆炸过程中电磁辐射特性的唯一因素，因此将2.1节得到的等离子体电导率代入全MHD模型，求解即可得到爆炸过程中的电磁辐射特性。由于爆炸过程中电导率空间分布不均匀，需要将其分为理想磁流体和非理想磁流体2种情况进行分析。
在理想磁流体情况下，忽略磁扩散和电磁波效应，由于磁场完全被冻结在等离子体中，因此磁场演化的时间尺度和流场运动的时间尺度完全一致。只需要对控制方程(3)中的前6个公式求解，即可得到理想磁流体情况下的各个物理量的变化。
在非理想磁流体情况下，由于电导率相对较低时的磁扩散和电磁波效应较为显著，需要同时考虑磁扩散、电磁波效应和求解稳定性条件，因此求解过程中的时间步长非常短，与流场求解时的时间步长存在量级差距，为了解决这个问题，求解过程被拆成2步。首先计算理想磁流体控制方程，求得并更新所有求解变量，再通过隐式时间积分法，在这个时间步内应用隐式算法求解方程(4)，之后再单独更新磁感应强度。
(4)
因为在前面的计算中已经得到了磁感应强度，仅考虑无耗无源情况，故需要考虑的Maxwell方程变为：
(5)
式中：E为电场强度。只需要求解方程(5)即可得电场强度。
3 数值模拟
随着仿真技术的发展，爆炸与高速冲击等强间断非线性问题出现了很多解决方法[14-16]。本文采用基于AUSM+-up通量分裂方法[18]的有限体积法对全MHD控制方程进行计算，给出通量描述如下：
(6)
式中：Fc为守恒通量； Fp为压力通量；FM为MHD通量； ρ为流体密度；U为速度；μe为磁导率，其余各量在式(3)中已给出计算公式。
AUSM+-up通量分裂只涉及理想磁流体方程的求解，对于电磁波和磁扩散的求解在后面用隐式方法给出。对于理想磁流体，将通量分解为守恒部分Fc、压力部分Fp和MHD部分FM，对于前两部分，采用标准的AUSM+-up通量分裂方法即可，对于FM，可按照粘性通量格式进行处理，即：
(7)
式中：下标g表示单元界面通量；i表示界面左单元；j表示界面右单元，最终界面通量为：
(8)
对于磁扩散技术过高的情况，需要在一个显式时间积分步内，采用隐式方法求解。考虑时间离散后，给出完备的离散格式为：
(9)
历遍所有单元界面，组成线性方程组，通过MKL中的PardisoLU分解稀疏矩阵，求解器直接求得最终的磁感应强度。
对比前后2个时间步内得到的磁感应强度，可得到磁感应强度的变化率，在网格上进行积分为：
(10)
式中: N为网格边界数目；ni为界面i处的单位外法线向量；
为界面i处电场强度；dsi为界面i的面积。同理可求得最终的电场强度。
为了尽量减少影响因素，开展了单纯的空中爆炸数值模拟，忽略地面反射，炸药为正方形。在空中爆炸的情况下，可以定量分析产生的电磁场参数和炸药当量之间的关系。以自然磁场作为背景磁场，进行爆炸产生的电磁辐射数值模拟研究。取自然磁场磁感应强度为0.5 Gs，其余变量取国际单位制，计算时间取0.2 ms。
炸药采用瞬时起爆模型，状态方程采用理想气体状态方程，有：
p=(γ-1)ρe
(11)
炸药和空气的参数如表1所示，其中ρ0为初始密度，e0为初始比内能，γ为比热容比，T0为初始温度，μe为磁导率。
表1 数值模拟材料参数
Mate-rialρ0/(kg·m-3)e0/(J·kg-1)γT0/Kμe/(T·m·A-1)explosive1 7178.5×1061.4300.01.257×10-6air1.29196 3661.4300.01.257×10-6
二维的计算网格模型如图1所示，计算域为边长10 m的正方形，中心位置放置炸药。不同当量(50 g、200 g、800 g和3.2 kg)的B炸药尺寸如表2所示，网格划分采用中心密集，边界稀疏，保证计算能准确描述爆炸初始阶段各个物理量变化率较大的情况。整体的网格采用非结构网格，在边界处约束磁场强度为初始值，设置初始磁场为x方向，大小为0.5 Gs，对流动过程的边界设置为流出边界。计算过程可以获得模型关键点A-F的时程曲线如图1所示。
表2 炸药质量及尺寸
Mass/gSide length/cm503.392 82005.385 88008.549 03 20013.571 4
图1 数值模拟计算模型示意图
4 数值模拟结果分析
4.1 电磁场变化情况
图2(b)给出了3.2 kg B炸药爆炸时产生在爆心处电场扰动的数值模拟结果。由于在计算中采用了瞬时爆轰模型，忽略了起爆到炸药完全爆轰的时间，因此在模拟计算一开始便产生了电场扰动信号，考虑到炸药的边长为13.57 cm，其对角线的一半长度为9.6 cm，从起爆到完全爆轰约为10 μs，因此在数值模拟结果中将-10 μs作为起爆时间，信号的出现时间与文献[7]符合较好。从整体的波形来看，数值模拟结果的波形和Kual[7]给出的理论计算结果相似，同时和任会兰[17]给出的实验结果在电场脉冲持续时间和波形上相似，由此可以证明本文中的数值模拟结果的可靠性。
图2 电场强度理论、数值模拟和实验结果曲线
图3给出了800 g B炸药爆轰完全后30 μs内x方向磁感应强度Bx的演化过程。
从图3数值模拟的结果来看，在炸药爆轰完全瞬间，由于等离子体的运动造成了磁场的剧烈变化，爆炸中心处的磁场出现了明显的下降，随着爆轰产物的运动，在y方向上形成了2个明显的磁场峰值区域，并随着时间的推移，这2个磁场峰值区域之间的距离逐渐增大。这主要是因为爆轰产物运动初期，由于高温高压爆轰产物能量密度高，形成的等离子体电导率高，处于理想磁流体情况下，磁场被冻结在等离子体内，与等离子体即爆轰产物一起运动，因此在y方向上产生了2个明显的峰值区域，并随着爆轰产物一起运动。
图3 不同时刻x方向磁感应强度Bx数值模拟云图
图4则给出了800 g B炸药爆轰完全后50 μs时x方向的流动速度、电导率、比内能及x方向磁感应强度Bx的空间分布。从上一节的理论可知，电导率的分布主要取决于流场中的能量分布，从图4(b)、图4(c)中也可以看出，电导率的峰值主要出现在流场速度和比内能最大的地方，这个结果与理论互相印证。对比图4(a)、图4(b)和图4(c)，可以发现：电导率与比内能高的地方正好是爆轰产物运动的前沿，在这些地方，由于高电导率导致磁扩散度很小，磁场被冻结在等离子体内，磁场中的扰动可以随时间积累起来。而其他地方电导率极低，磁扩散度很大，磁场很难形成明显的积累，产生的磁场扰动以很快的速度向周围扩散，从而难以形成明显的磁场扰动积累。
图4 50 μs数值模拟云图
图5给出了图1中关键点A-E上的磁感应强度Bx随时间和距离变化的情况。关键点A为爆心处，关键点E离爆心最远。
图5 各关键点处磁场扰动曲线
从图5(a)可以看出，炸药爆炸后磁感应强度先增长，在到达峰值时刻后，磁场扰动逐渐降低，到后来趋于平衡。这是因为在爆轰产物运动初期，由于其能量密度高、电导率高、磁扩散度低，磁场被冻结在爆轰产物内，产生的磁场扰动可随时间进行积累。而在爆轰产物膨胀到一定程度后，由于能量密度的下降导致电导率显著降低，从而引起磁扩散度升高，磁场积累的扰动逐渐向外扩散出去。从图5(b)可以看出，关键点离爆心越远，磁场扰动峰值越小，这符合电磁场扩散中的能量守恒原理。图6为800 g B 炸药爆炸后关键点B处的电场和磁场变化情况。
图6 关键点B处的电磁场变化曲线
从图6可以看出，磁感应强度先上升后下降并逐渐趋于平衡。电场强度先上升后下降，下降至零后有向负向的增长，增长至第2个峰值后逐渐减少并趋于平衡。
电场峰值出现在磁场上升沿斜率最大处，即磁场变化最剧烈的时刻，而当磁场峰值出现后，电场正好有正转负。这符合感生电场与磁场之间的规律。由于磁场的变化为上升沿陡，峰值之后的下降速度比上升慢，因此导致电场在负向的峰值小于正向。随着时间的变化，电场和磁场都最终趋于平衡。
4.2 炸药当量对爆炸产生的电磁场的影响
不同当量炸药爆炸产生磁场的数值模拟结果如图7所示。不同当量炸药爆炸产生电场的数值模拟结果如图8所示。
图7(a)为对不同当量(50 g、200 g、800 g、3.2 kg)炸药爆炸在1.5 m处产生磁场的数值模拟结果。从结果来看，炸药当量越大，产生的磁场扰动峰值越大，脉冲持续时间越长，磁场扰动的峰值时间也有所推迟。这主要是由于当药量增大时，爆炸气体产物维持高能量状态的时间也相应增大，电导率维持在较高水平的时间加长，使产生的磁场扰动积累时间增加，相应的磁场峰值也会增加，磁场向外扩散的时间增大。
图7(b)是1.5 m处磁感应强度峰值随炸药当量之间的关系。取拟合关系式为：
Bx=a+b×wm
(12)
从拟合结果来看，当量指数m为0.633±0.067，其相关性为0.932。在此拟合基础上，可认为磁场扰动峰值与炸药当量的2/3次方成正比。
图7 不同当量炸药爆炸产生磁场的数值模拟曲线
图8 不同当量炸药爆炸产生电场的数值模拟曲线
图8(a)为对不同当量(50 g、200 g、800 g、3.2 kg)B炸药爆炸在1.5 m处产生的电场扰动进行数值模拟的结果。从数值模拟的结果来看，电场与磁场扰动随炸药当量的变化趋势相同，即当量越大，电场扰动峰值越大，峰值到达时间略有推迟，具体原因分析与磁场相似，在此不再赘述。取拟合公式为：
Ez=a+b×wn
(13)
从拟合结果来看，当量指数n为0.341±0.031，其相关性为0.918。在此拟合基础上，可认为磁场扰动峰值与炸药当量的1/3次方成正比，这也和文献[9]中的实验结果相一致。
5 结论
1) 基于本文提出的磁流体动力学与麦克斯韦方程结合的方法，通过数值求解可以获得炸药爆炸过程中电磁辐射的结果，具有参考价值。
2) 常规炸药爆炸会产生电磁辐射，其中磁场信号只有一个峰值，上升沿比下降沿陡。电场信号有2个峰值，正向的峰值大于负向的峰值。电场信号和磁场信号持续时间相当。电场信号峰值出现时间为磁场信号上升沿中部，为整个爆轰过程中磁场变化最剧烈的时刻。
3) 对于同一种炸药而言，当炸药当量增大时，其爆炸产生磁信号峰值、电信号峰值、信号持续时间和信号峰值到达时间都随之增大，磁信号峰值与当量的2/3次方成正比，电信号峰值与当量1/3次方成正比。
[1] 戴晴，李传胪，陈国强，等.低温等离子体激励宽带电磁波信号的试验研究[J].电子信息对抗技术，2009，24(05):72-74.
[2] Kolsky H.Electromagnetic waves emitted on detonation of explosives[J].Nature，1954，173:77.
[3] Cook M A.The Science of High Explosives[M].New York:Reinhold Publishing，1958:440.
[4] Boronin A P，Velmin V A，Medvedev Y A，et al.Experimental study of the electromagnetic field in the near zone of explosions produced by solid explosives[J].Zhurnal Prikladnoi Mekhaniki I Tekhnicheskoi Fiziki，1968，9(06):99-103.
[5] Boronin A P，Kapinos V N，Krenev S A，et al.Physical mechanism of electromagnetic field generation during the explosion of condensed explosive charges[J].Combustion Explosion and Shock Waves，1990，26:597-602.
[6] Fine J E.Estimates of the electromagnetic radiation from detonation of conventional explosives[R].Army Research Laboratory Report，ARL-TR-2447，2001.
[7] Kuhl A L，White D A，Kirkendall B A.Electromagnetic waves from TNT explosions[J].Journal of Electromagnetic Analysis and Applications，2014(06):280-295.
[8] 陈生玉，孙新利，钱世平，等.化爆引起的电磁辐射[J].爆炸与冲击，1997，17(04):363-368.
[9] 曹景阳，谢树果，苏东林，等.航天火工品爆炸引起的电磁干扰测量[J].北京航空航天大学学报，2011，37(11):1384-1394.
[10] 王长利，李迅，刘晓新，等.典型炸药爆炸过程的电磁辐射实验研究[J].兵工学报，2014，35(S2):188-192.
[11] 陈鸿，何勇，潘绪超，等.RDX 基含铝炸药爆炸电磁辐射信号特性实验研究[J].含能材料，2020，28(05):475-482.
[12] Lee N，Close S，Coel A，et al.Theory and experiments characterizing hypervelocity impact plasmas on biased spacecraft materials[J].Physics of Plasmas，2013，20(03):139-311.
[13] 栗建桥，马天宝，宁建国.爆炸对自然磁场干扰机理[J].力学学报，2018，50(05):1206-1218.
[14] 马天宝，任会兰，李健，等.爆炸与冲击问题的大规模高精度计算[J].力学学报，2016，48(03):599-608.
[15] 连赟猛，吴美云，顾晓辉.密闭装置内爆炸及数值模拟精度分析[J].兵器装备工程学报，2016(08):169-172.
[16] Zeng Q L，Aydemir N U，Lien F S，et al.Extension of staggeredgrid-based AUSM-family schemes for use in nuclear safety analysis codes[J].International Journal of Multiphase Flow，2017，93:17-32.
[17] 任会兰，储著鑫，栗建桥，等.B炸药爆炸过程中电磁辐射研究[J].力学学报，2020，52(04):1199-1210.
[18] Liou M S.A sequel to AUSM，Part II:AUSM +-up for all speeds[J].Journal of Computational Physics，2006，214(01):137-170.