基于流体-磁流体-粒子混合方法的高空核爆炸碎片云模拟

彭国良 张俊杰

1)(北京理工大学机电学院,北京 100081)

2)(西北核技术研究所,西安 710024)

提出了描述高空核爆炸碎片云运动的流体-磁流体-粒子(particle-in-cell,PIC)混合模型,相较目前的主流模型,该模型能够计算更加广泛的空间尺度.根据碎片云运动涉及的高温离子、低温离子和中性大气的不同性质,采用三种模型进行联合求解:高温离子用PIC 粒子模型计算,低温离子用磁流体模型计算,中性大气用流体模型计算,并将三者之间的相互作用作为源项加入相应的控制方程.最后,计算了美国Starfish 试验中碎片云的扩展情况,与试验结果进行了比对,并验证了求解方案的可靠性.此外,还给出了不同投影角度下碎片云形状随时间的变化,并分析了影响碎片云运动的主要因素,包括大气阻力、磁压、槽型不稳定性和霍尔电流等.

1 引 言

碎片云是高空核爆炸重要的环境现象之一.高空核爆后,弹体材料和裂变产物在高温下形成高速运动的等离子体碎片云.碎片离子在扩展时沿垂直磁场和平行磁场方向产生性质不同的激波[1]和不稳定性[2,3],最终损失能量在大气中沉积形成槽型结构[4−6].高空核爆产生的效应与地面核爆完全不同[7−9].向高空运动的裂变碎片会持续衰变并产生大量电子,这些电子可以运动到同步卫星的轨道高度[10,11];与此同时,向地面运动的碎片等离子体可以电离并加热低海拔处的大气,并产生缓慢变化的电磁信号[12−14].高空核爆碎片云的运动和衰变形成的附加电离层、人工辐射带等地球物理现象可以对在轨航天器和地面电站构成威胁[7].因此,研究高空核爆炸碎片云的运动在物理机理和国防安全两个方面均有重要意义.

对爆高在100 km 以下的碎片云运动,由于大气密度较高,碎片粒子自由程较短,且磁压远小于大气压,通常可采用雪犁模型等流体力学方法计算[15,16].对爆高在100 km 以上的碎片云运动,杨斌等[17]和John 等[18]采用磁流体力学模型,能够解释部分试验观测结果.核爆炸早期的碎片云为热等离子体,初始运动速度一般在1000 km/s 量级[18].磁流体模型的主要问题在于当大气密度较低时,碎片自由程极大,不满足流体的局域热平衡条件.为克服这个问题,混合模型[19]被引入到碎片云的计算中.混合模型用PIC(particle in cell)方法描述离子的运动,将电子假设为无质量流体,再结合电中性条件来建立运动方程.Brecht 等[20,21]应用混合模型研究了扩散碎片与大气离子的碰撞,对观测到的核爆现象进行了模拟,研究了碎片速度谱、离子化水平及辐射带中的动力学问题,并进一步模拟了海拔400 km高度1 Mt 当量核爆炸碎片云的扩展情况,得到多种核素离子的较精细分布图像.此外,Morrow[22]也利用2.5 维混合模型ZMR(流体在XOY平面、磁场沿Z轴)计算了碎片云的无碰撞激波和电荷交换.不过,混合模型在电子密度较低、磁场较高的近真空区域计算电磁场时容易发生数值发散,常见的处理方法(如强制设置高密度背景电子)[23]会导致非物理的结果.目前,混合模型依然无法处理电荷高度分离时产生的电磁效应,因而无法描述部分高频等离子体不稳定过程[24],这些过程对末期碎片离子的分布、激波形成及高轨道区域的电子分布都有直接影响.

综合上面三种模型的优缺点,本文参考Bai等[25]计算宇宙射线的方法,提出了一种流体-磁流体-PIC 混合方法.高温的碎片等离子体用混合模型进行描述,低温的大气等离子体用磁流体进行描述,中性大气成分用流体力学进行描述.最后我们模拟了Starfish 试验的过程,并给出了与试验一致的碎片云空间分布.

2 计算模型和方法

在模拟碎片云运动时,离子分为两类:高温的碎片离子和低温的环境离子.高温的碎片离子可以使用混合模型进行描述.混合模型兼顾了计算资源和离子碎片的粒子属性,能够给出低频电荷分离效应.对于低温的环境等离子体,为避免混合粒子模拟中近真空区域计算导致的困难,我们用磁流体模型进行描述.虽然在大气极端稀疏时粒子的动理学(kinetic)效应远大于流体效应(Knudsen’s 数大于1),如果我们关注的主要是碎片离子的大体分布,依然可以使用流体的连续性方程来描述该过程,以此保证计算时不会出现真空区域.另外,为考虑中性大气对离子运动的影响,还需计算中性大气的运动.可以看到流体-磁流体-PIC 混合模型可以保证计算能够自洽地覆盖更广泛的空间尺度.

碎片离子用PIC 模型描述,运动方程为

式中,m为质量,v为速度,q为碎片离子携带的电荷(本文中假设所有离子均为一价离子),E和B表示电场和磁场,f为碰撞阻力,n为电荷密度,下标e 表示电子,下标d 表示碎片离子,下标i 表示环境离子,下标a 表示中性大气,Seb,Sn,Sef分别表示束缚电子阻止本领、原子核阻止本领、自由电子阻止本领.阻止本领计算公式可以参考文献[26].粒子的初速度假设服从麦克斯韦分布,则给出初始碎片温度就可得到粒子的速度分布.为保证稳定性,运动方程的时间推进采用半隐式格式:

式中,上标n代表第n时间步的值.

环境离子用磁流体模型描述,由于环境磁流体温度较低,大部分空间位置的磁流体密度极低,其压强相比于磁压是小量,很容易因数值误差造成压强为负值导致计算无法进行,本文中对磁流体不直接计算能量方程.当大气密度较大时,可设离子温度与大气温度相等;而大气密度较小时,离子密度也较小,设离子温度与大气温度相等导致的压强误差对磁流体的运动也不会产生大的影响.环境离子控制方程为

式中,ρ为密度;k为玻尔兹曼常数;µ0为真空磁导率;T为环境温度;υia为离子和大气的碰撞频率,计算公式为

为便于与常用磁流体算法格式统一,对(3)式进行化简.假设等离子体满足电中性条件:

等离子电流可表示为

式中,jd为碎片离子电流.由广义欧姆定律

可得到电场表达式

式中,η为电阻.令

E0即为普通磁流体方程中的电场,Fd为PIC 方程导致的附加项,则有

将(11)式和(12)式代入环境离子磁流体方程(3),则运动方程和磁场方程可改写为

与普通磁流体方程相比,仅增加了源项,故可用常见的磁流体求解格式来计算,本文中利用HLL 有限体积格式来求解低温磁流体方程[27].初始的环境离子包括空间环境离子和X射线电离离子,由于X射线电离的初级光电子在海拔300 km 以上的自由程很长,会很快逃逸,对环境离子密度贡献很小,故忽略300 km 以上X射线电离对环境离子密度的贡献.大气层的空间离子密度可由国际参考电离层模型得到.X射线电离离子密度的计算公式为

式中,P(ν) 为X射线能谱,本文取为1keV黑体谱;µ(ν)为大气对X射线的谱吸收系数,取值可参考文献[28];r为距爆心的距离;ρair为大气密度;E为X射线总能量,一般取总爆炸能量的70%;Iion为平均电离能,一般取33 eV/ion.

中性大气满足流体方程:

式中,R为气体常数,g为重力加速度,I为中性大气的能量密度.初始的温度按标准大气模型选取,初始密度则需从标准大气中减去离子密度.计算方法选常见的HLL 有限体积格式.

值得注意的是,实际的物理过程中还包含离子复合、电荷交换、中性大气的电离等复杂过程,由于高空相关参数的缺乏,本文暂不予考虑.

3 分析与讨论

为验证计算模型,以高空核爆炸试验数据最完整的Starfish 试验为例进行计算.Starfish 试验爆炸地点为北纬17°,西经170°,高度为400 km,当量1.4 Mt.初始条件的选取参考文献[17,18],即初始碎片质量取1000 kg,碎片动能占总爆炸能量的20%,碎片原子质量数取铁原子的质量数56,电荷数取1,地磁场取爆点的偶极子磁场.计算区域取西东方向为X方向,计算范围为[–800 km,800 km],网格数240;南北方向为Y方向,计算范围为[–800 km,800 km],网格数240;沿高度方向为Z方向,计算范围为[100 km,800 km],网格数120;爆心坐标为[0,0,400 km].图1 为试验拍摄到的爆炸场景照片[18],计算得到的碎片密度在子午面上的投影面密度云图如图2 所示,沿磁力线方向的投影二维面密度云图如图3 所示,沿高度方向的投影二维面密度云图如图4 所示.文献[18]中给出的试验观测数据在爆后0.1 s 左右时碎片云水平扩展半径约200 km,碎片云底部约在220 km高度附近.由图2(a)和图4(a)可知,本文的计算结果在0.1 s 时碎片云在南北方向的半径略大于200 km,碎片云底部略大于250 km,与试验结果基本一致.X射线火球在80 km高度附近,会遮挡试验光学观测结果,但碎片缓发β射线在海拔70 km 左右的沉积区形状与爆炸区沿磁力线方向的投影较为相似,利用这个特点可以通过观测沉积区形状反推碎片云形状.观测数据表明0.2 s 时碎片云南北方向水平半径小于500 km,下边界在180 km 附近.由图2(b)和图4(b)可知,在0.2 s时本文计算的碎片云在南北方向的半径大约为400 km,碎片云底部的半径大约为150 km,与试验结果基本一致.

图1 Starfish 试验照片:爆后4 s β射线沉积区图像Fig.1.Photos of Starfish:β patch photo after burst time 4 s.

图4 不同时刻Starfish 试验模拟计算的三维碎片云密度沿高度方向的投影面密度云图 (a) 0.1 s;(b) 0.2 s;(c) 0.5 s;(d) 1 sFig.4.Contour plot of three-dimensional debris density in the XOY plane(perpendicular to altitude direction) in Starfish experiment at different time:(a) 0.1 s;(b) 0.2 s;(c) 0.5 s;(d) 1 s.

由图2 可知,碎片离子在垂直于磁场的方向受到磁场的约束和大气阻力的影响,而在沿磁场方向,离子基本不受磁场影响,因此碎片云的形状整体为两头尖、中间厚的“帽子形”,与试验观测结果一致.由于大气密度随高度增加而增加,故碎片云向下运动会受到更大的阻力,向上运动的碎片离子运动速度较快.文献[4,5]指出碎片运动时存在槽型不稳定性,且这种不稳定随离子回旋半径Rd与磁泡半径R之比变小而变得更强,其表达式为

图2 不同时刻Starfish 试验模拟计算的三维碎片云密度沿子午面方向的投影面密度云图 (a) 0.1 s;(b) 0.2 s;(c) 0.5 s;(d) 1 sFig.2.Contour plot of three-dimensional debris density in the geomagnetic meridian plane in Starfish experiment at different time:(a) 0.1 s;(b) 0.2 s;(c) 0.5 s;(d) 1 s.

式中,N为碎片离子数密度.显然,离子速度越大、密度越小时,更易发生槽型不稳定性.0.5 s 以后,等离子体受槽型不稳定性的影响,向上运动的碎片离子已散开.

设磁倾角为θ,则沿磁力线方向投影时坐标变换关系为

式中,x0为西东方向,y0为南北方向,z0为海拔高度,下标1 表示投影坐标系下的坐标.由图3 可知,碎片云沿磁场方向的投影基本为椭圆,与图1 所示沉积区观测的形状一致.由于大气密度的差异,碎片云向左上方向运动阻力最小,向右下方向运动阻力越大,因此碎片云椭圆的偏心率会越来越大.

图3 不同时刻Starfish 试验模拟计算的三维碎片云密度沿磁力线方向的投影面密度云图 (a) 0.1 s;(b) 0.2 s;(c) 0.5 s;(d) 1 sFig.3.Contour plot of three-dimensional debris density in the XOY plane(perpendicular to the geomagnetic fields) in Starfish experiment at different time:(a) 0.1 s;(b) 0.2 s;(c) 0.5 s;(d) 1 s.

由图4 可知,碎片云在西东方向受磁场约束较大,扩展速度明显慢于南北方向.碎片云在西东方向的扩展并不是对称的,向西方向速度大于向东方向的速度,出现不稳定性的时间也较早.考虑霍尔电流的作用时,电场表达式为

式中,v为离子的平均速度,vH为考虑霍尔电流时的等效速度,J/ne也是速度的量纲,本文中称其为霍尔速度(Hall velocity).初始时刻,垂直于X轴(西东方向)的磁子午面是对称面,即

则电流项为

等效速度

可见,当霍尔速度相对离子速度不可忽略时,等效速度不再具有对称性,进而导致电场、磁场、速度场等不再具备对称性.

图5 给出了爆后0.2 s 时子午面上x0方向离子速度与霍尔速度的分布,可以看到离子速度较大的地方霍尔速度也较大,且霍尔速度略小于离子速度,但其影响不可忽略.这表明在高空核爆炸中必须考虑霍尔效应,在霍尔效应作用下各场量空间对称性被破坏.

图5 爆后0.2 s 子午面上的x0 方向速度分布: (a)离子速度;(b)霍尔速度Fig.5.Velocity distributions in the x0 direction in meridian plane at t=0.2 s:(a) Ion velocity;(b) Hall velocity.

4 结 论

将高空核爆炸碎片云运动涉及的流体分类描述,用适用于高温低密度流体的PIC 粒子模型描述高温碎片云,用适用于低温流体的磁流体模型描述环境离子和电磁场,用流体方程描述中性大气,进而得到了碎片云运动的计算模型.利用该模型计算了Starfish 试验的碎片云运动,给出了不同观测方向的碎片云运动规律,并与试验观测结果进行比对,验证了模型的合理性.进一步的分析表明,碎片云向下运动主要受地磁场和大气阻力作用;向上运动主要受地磁场作用,且由于运动速度较快,易发生槽型不稳定性;由于霍尔电流的影响,东西方向的运动是不对称的,西侧的运动速度会大于东侧.

感谢磷虾科技束庆邦博士在编程方面的帮助及中国科技大学陆全明教授关于混合模拟的讨论.