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

彭國良 張俊杰

1)(北京理工大學(xué)機(jī)電學(xué)院,北京 100081)

2)(西北核技術(shù)研究所,西安 710024)

提出了描述高空核爆炸碎片云運(yùn)動的流體-磁流體-粒子(particle-in-cell,PIC)混合模型,相較目前的主流模型,該模型能夠計算更加廣泛的空間尺度.根據(jù)碎片云運(yùn)動涉及的高溫離子、低溫離子和中性大氣的不同性質(zhì),采用三種模型進(jìn)行聯(lián)合求解:高溫離子用PIC 粒子模型計算,低溫離子用磁流體模型計算,中性大氣用流體模型計算,并將三者之間的相互作用作為源項(xiàng)加入相應(yīng)的控制方程.最后,計算了美國Starfish 試驗(yàn)中碎片云的擴(kuò)展情況,與試驗(yàn)結(jié)果進(jìn)行了比對,并驗(yàn)證了求解方案的可靠性.此外,還給出了不同投影角度下碎片云形狀隨時間的變化,并分析了影響碎片云運(yùn)動的主要因素,包括大氣阻力、磁壓、槽型不穩(wěn)定性和霍爾電流等.

1 引 言

碎片云是高空核爆炸重要的環(huán)境現(xiàn)象之一.高空核爆后,彈體材料和裂變產(chǎn)物在高溫下形成高速運(yùn)動的等離子體碎片云.碎片離子在擴(kuò)展時沿垂直磁場和平行磁場方向產(chǎn)生性質(zhì)不同的激波[1]和不穩(wěn)定性[2,3],最終損失能量在大氣中沉積形成槽型結(jié)構(gòu)[4?6].高空核爆產(chǎn)生的效應(yīng)與地面核爆完全不同[7?9].向高空運(yùn)動的裂變碎片會持續(xù)衰變并產(chǎn)生大量電子,這些電子可以運(yùn)動到同步衛(wèi)星的軌道高度[10,11];與此同時,向地面運(yùn)動的碎片等離子體可以電離并加熱低海拔處的大氣,并產(chǎn)生緩慢變化的電磁信號[12?14].高空核爆碎片云的運(yùn)動和衰變形成的附加電離層、人工輻射帶等地球物理現(xiàn)象可以對在軌航天器和地面電站構(gòu)成威脅[7].因此,研究高空核爆炸碎片云的運(yùn)動在物理機(jī)理和國防安全兩個方面均有重要意義.

對爆高在100 km 以下的碎片云運(yùn)動,由于大氣密度較高,碎片粒子自由程較短,且磁壓遠(yuǎn)小于大氣壓,通?？刹捎醚├缒Ｐ偷攘黧w力學(xué)方法計算[15,16].對爆高在100 km 以上的碎片云運(yùn)動,楊斌等[17]和John 等[18]采用磁流體力學(xué)模型,能夠解釋部分試驗(yàn)觀測結(jié)果.核爆炸早期的碎片云為熱等離子體,初始運(yùn)動速度一般在1000 km/s 量級[18].磁流體模型的主要問題在于當(dāng)大氣密度較低時,碎片自由程極大,不滿足流體的局域熱平衡條件.為克服這個問題,混合模型[19]被引入到碎片云的計算中.混合模型用PIC(particle in cell)方法描述離子的運(yùn)動,將電子假設(shè)為無質(zhì)量流體,再結(jié)合電中性條件來建立運(yùn)動方程.Brecht 等[20,21]應(yīng)用混合模型研究了擴(kuò)散碎片與大氣離子的碰撞,對觀測到的核爆現(xiàn)象進(jìn)行了模擬,研究了碎片速度譜、離子化水平及輻射帶中的動力學(xué)問題,并進(jìn)一步模擬了海拔400 km高度1 Mt 當(dāng)量核爆炸碎片云的擴(kuò)展情況,得到多種核素離子的較精細(xì)分布圖像.此外,Morrow[22]也利用2.5 維混合模型ZMR(流體在XOY平面、磁場沿Z軸)計算了碎片云的無碰撞激波和電荷交換.不過,混合模型在電子密度較低、磁場較高的近真空區(qū)域計算電磁場時容易發(fā)生數(shù)值發(fā)散,常見的處理方法(如強(qiáng)制設(shè)置高密度背景電子)[23]會導(dǎo)致非物理的結(jié)果.目前,混合模型依然無法處理電荷高度分離時產(chǎn)生的電磁效應(yīng),因而無法描述部分高頻等離子體不穩(wěn)定過程[24],這些過程對末期碎片離子的分布、激波形成及高軌道區(qū)域的電子分布都有直接影響.

綜合上面三種模型的優(yōu)缺點(diǎn),本文參考Bai等[25]計算宇宙射線的方法,提出了一種流體-磁流體-PIC 混合方法.高溫的碎片等離子體用混合模型進(jìn)行描述,低溫的大氣等離子體用磁流體進(jìn)行描述,中性大氣成分用流體力學(xué)進(jìn)行描述.最后我們模擬了Starfish 試驗(yàn)的過程,并給出了與試驗(yàn)一致的碎片云空間分布.

2 計算模型和方法

在模擬碎片云運(yùn)動時,離子分為兩類:高溫的碎片離子和低溫的環(huán)境離子.高溫的碎片離子可以使用混合模型進(jìn)行描述.混合模型兼顧了計算資源和離子碎片的粒子屬性,能夠給出低頻電荷分離效應(yīng).對于低溫的環(huán)境等離子體,為避免混合粒子模擬中近真空區(qū)域計算導(dǎo)致的困難,我們用磁流體模型進(jìn)行描述.雖然在大氣極端稀疏時粒子的動理學(xué)(kinetic)效應(yīng)遠(yuǎn)大于流體效應(yīng)(Knudsen’s 數(shù)大于1),如果我們關(guān)注的主要是碎片離子的大體分布,依然可以使用流體的連續(xù)性方程來描述該過程,以此保證計算時不會出現(xiàn)真空區(qū)域.另外,為考慮中性大氣對離子運(yùn)動的影響,還需計算中性大氣的運(yùn)動.可以看到流體-磁流體-PIC 混合模型可以保證計算能夠自洽地覆蓋更廣泛的空間尺度.

碎片離子用PIC 模型描述,運(yùn)動方程為

式中,m為質(zhì)量,v為速度,q為碎片離子攜帶的電荷(本文中假設(shè)所有離子均為一價離子),E和B表示電場和磁場,f為碰撞阻力,n為電荷密度,下標(biāo)e 表示電子,下標(biāo)d 表示碎片離子,下標(biāo)i 表示環(huán)境離子,下標(biāo)a 表示中性大氣,Seb,Sn,Sef分別表示束縛電子阻止本領(lǐng)、原子核阻止本領(lǐng)、自由電子阻止本領(lǐng).阻止本領(lǐng)計算公式可以參考文獻(xiàn)[26].粒子的初速度假設(shè)服從麥克斯韋分布,則給出初始碎片溫度就可得到粒子的速度分布.為保證穩(wěn)定性,運(yùn)動方程的時間推進(jìn)采用半隱式格式:

式中,上標(biāo)n代表第n時間步的值.

環(huán)境離子用磁流體模型描述,由于環(huán)境磁流體溫度較低,大部分空間位置的磁流體密度極低,其壓強(qiáng)相比于磁壓是小量,很容易因數(shù)值誤差造成壓強(qiáng)為負(fù)值導(dǎo)致計算無法進(jìn)行,本文中對磁流體不直接計算能量方程.當(dāng)大氣密度較大時,可設(shè)離子溫度與大氣溫度相等;而大氣密度較小時,離子密度也較小,設(shè)離子溫度與大氣溫度相等導(dǎo)致的壓強(qiáng)誤差對磁流體的運(yùn)動也不會產(chǎn)生大的影響.環(huán)境離子控制方程為

式中,ρ為密度;k為玻爾茲曼常數(shù);μ0為真空磁導(dǎo)率;T為環(huán)境溫度;υia為離子和大氣的碰撞頻率,計算公式為

為便于與常用磁流體算法格式統(tǒng)一,對(3)式進(jìn)行化簡.假設(shè)等離子體滿足電中性條件:

等離子電流可表示為

式中,jd為碎片離子電流.由廣義歐姆定律

可得到電場表達(dá)式

式中,η為電阻.令

E0即為普通磁流體方程中的電場,Fd為PIC 方程導(dǎo)致的附加項(xiàng),則有

將(11)式和(12)式代入環(huán)境離子磁流體方程(3),則運(yùn)動方程和磁場方程可改寫為

與普通磁流體方程相比,僅增加了源項(xiàng),故可用常見的磁流體求解格式來計算,本文中利用HLL 有限體積格式來求解低溫磁流體方程[27].初始的環(huán)境離子包括空間環(huán)境離子和X射線電離離子,由于X射線電離的初級光電子在海拔300 km 以上的自由程很長,會很快逃逸,對環(huán)境離子密度貢獻(xiàn)很小,故忽略300 km 以上X射線電離對環(huán)境離子密度的貢獻(xiàn).大氣層的空間離子密度可由國際參考電離層模型得到.X射線電離離子密度的計算公式為

式中,P(ν) 為X射線能譜,本文取為1keV黑體譜;μ(ν)為大氣對X射線的譜吸收系數(shù),取值可參考文獻(xiàn)[28];r為距爆心的距離;ρa(bǔ)ir為大氣密度;E為X射線總能量,一般取總爆炸能量的70%;Iion為平均電離能,一般取33 eV/ion.

中性大氣滿足流體方程:

式中,R為氣體常數(shù),g為重力加速度,I為中性大氣的能量密度.初始的溫度按標(biāo)準(zhǔn)大氣模型選取,初始密度則需從標(biāo)準(zhǔn)大氣中減去離子密度.計算方法選常見的HLL 有限體積格式.

值得注意的是,實(shí)際的物理過程中還包含離子復(fù)合、電荷交換、中性大氣的電離等復(fù)雜過程,由于高空相關(guān)參數(shù)的缺乏,本文暫不予考慮.

3 分析與討論

為驗(yàn)證計算模型,以高空核爆炸試驗(yàn)數(shù)據(jù)最完整的Starfish 試驗(yàn)為例進(jìn)行計算.Starfish 試驗(yàn)爆炸地點(diǎn)為北緯17°,西經(jīng)170°,高度為400 km,當(dāng)量1.4 Mt.初始條件的選取參考文獻(xiàn)[17,18],即初始碎片質(zhì)量取1000 kg,碎片動能占總爆炸能量的20%,碎片原子質(zhì)量數(shù)取鐵原子的質(zhì)量數(shù)56,電荷數(shù)取1,地磁場取爆點(diǎn)的偶極子磁場.計算區(qū)域取西東方向?yàn)閄方向,計算范圍為[–800 km,800 km],網(wǎng)格數(shù)240;南北方向?yàn)閅方向,計算范圍為[–800 km,800 km],網(wǎng)格數(shù)240;沿高度方向?yàn)閆方向,計算范圍為[100 km,800 km],網(wǎng)格數(shù)120;爆心坐標(biāo)為[0,0,400 km].圖1 為試驗(yàn)拍攝到的爆炸場景照片[18],計算得到的碎片密度在子午面上的投影面密度云圖如圖2 所示,沿磁力線方向的投影二維面密度云圖如圖3 所示,沿高度方向的投影二維面密度云圖如圖4 所示.文獻(xiàn)[18]中給出的試驗(yàn)觀測數(shù)據(jù)在爆后0.1 s 左右時碎片云水平擴(kuò)展半徑約200 km,碎片云底部約在220 km高度附近.由圖2(a)和圖4(a)可知,本文的計算結(jié)果在0.1 s 時碎片云在南北方向的半徑略大于200 km,碎片云底部略大于250 km,與試驗(yàn)結(jié)果基本一致.X射線火球在80 km高度附近,會遮擋試驗(yàn)光學(xué)觀測結(jié)果,但碎片緩發(fā)β射線在海拔70 km 左右的沉積區(qū)形狀與爆炸區(qū)沿磁力線方向的投影較為相似,利用這個特點(diǎn)可以通過觀測沉積區(qū)形狀反推碎片云形狀.觀測數(shù)據(jù)表明0.2 s 時碎片云南北方向水平半徑小于500 km,下邊界在180 km 附近.由圖2(b)和圖4(b)可知,在0.2 s時本文計算的碎片云在南北方向的半徑大約為400 km,碎片云底部的半徑大約為150 km,與試驗(yàn)結(jié)果基本一致.

圖1 Starfish 試驗(yàn)照片:爆后4 s β射線沉積區(qū)圖像Fig.1.Photos of Starfish:β patch photo after burst time 4 s.

圖4 不同時刻Starfish 試驗(yàn)?zāi)M計算的三維碎片云密度沿高度方向的投影面密度云圖 (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 可知,碎片離子在垂直于磁場的方向受到磁場的約束和大氣阻力的影響,而在沿磁場方向,離子基本不受磁場影響,因此碎片云的形狀整體為兩頭尖、中間厚的“帽子形”,與試驗(yàn)觀測結(jié)果一致.由于大氣密度隨高度增加而增加,故碎片云向下運(yùn)動會受到更大的阻力,向上運(yùn)動的碎片離子運(yùn)動速度較快.文獻(xiàn)[4,5]指出碎片運(yùn)動時存在槽型不穩(wěn)定性,且這種不穩(wěn)定隨離子回旋半徑Rd與磁泡半徑R之比變小而變得更強(qiáng),其表達(dá)式為

圖2 不同時刻Starfish 試驗(yàn)?zāi)M計算的三維碎片云密度沿子午面方向的投影面密度云圖 (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為碎片離子數(shù)密度.顯然,離子速度越大、密度越小時,更易發(fā)生槽型不穩(wěn)定性.0.5 s 以后,等離子體受槽型不穩(wěn)定性的影響,向上運(yùn)動的碎片離子已散開.

設(shè)磁傾角為θ,則沿磁力線方向投影時坐標(biāo)變換關(guān)系為

式中,x0為西東方向,y0為南北方向,z0為海拔高度,下標(biāo)1 表示投影坐標(biāo)系下的坐標(biāo).由圖3 可知,碎片云沿磁場方向的投影基本為橢圓,與圖1 所示沉積區(qū)觀測的形狀一致.由于大氣密度的差異,碎片云向左上方向運(yùn)動阻力最小,向右下方向運(yùn)動阻力越大,因此碎片云橢圓的偏心率會越來越大.

圖3 不同時刻Starfish 試驗(yàn)?zāi)M計算的三維碎片云密度沿磁力線方向的投影面密度云圖 (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 可知,碎片云在西東方向受磁場約束較大,擴(kuò)展速度明顯慢于南北方向.碎片云在西東方向的擴(kuò)展并不是對稱的,向西方向速度大于向東方向的速度,出現(xiàn)不穩(wěn)定性的時間也較早.考慮霍爾電流的作用時,電場表達(dá)式為

式中,v為離子的平均速度,vH為考慮霍爾電流時的等效速度,J/ne也是速度的量綱,本文中稱其為霍爾速度(Hall velocity).初始時刻,垂直于X軸(西東方向)的磁子午面是對稱面,即

則電流項(xiàng)為

等效速度

可見,當(dāng)霍爾速度相對離子速度不可忽略時,等效速度不再具有對稱性,進(jìn)而導(dǎo)致電場、磁場、速度場等不再具備對稱性.

圖5 給出了爆后0.2 s 時子午面上x0方向離子速度與霍爾速度的分布,可以看到離子速度較大的地方霍爾速度也較大,且霍爾速度略小于離子速度,但其影響不可忽略.這表明在高空核爆炸中必須考慮霍爾效應(yīng),在霍爾效應(yīng)作用下各場量空間對稱性被破壞.

圖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 結(jié) 論

將高空核爆炸碎片云運(yùn)動涉及的流體分類描述,用適用于高溫低密度流體的PIC 粒子模型描述高溫碎片云,用適用于低溫流體的磁流體模型描述環(huán)境離子和電磁場,用流體方程描述中性大氣,進(jìn)而得到了碎片云運(yùn)動的計算模型.利用該模型計算了Starfish 試驗(yàn)的碎片云運(yùn)動,給出了不同觀測方向的碎片云運(yùn)動規(guī)律,并與試驗(yàn)觀測結(jié)果進(jìn)行比對,驗(yàn)證了模型的合理性.進(jìn)一步的分析表明,碎片云向下運(yùn)動主要受地磁場和大氣阻力作用;向上運(yùn)動主要受地磁場作用,且由于運(yùn)動速度較快,易發(fā)生槽型不穩(wěn)定性;由于霍爾電流的影響,東西方向的運(yùn)動是不對稱的,西側(cè)的運(yùn)動速度會大于東側(cè).

感謝磷蝦科技束慶邦博士在編程方面的幫助及中國科技大學(xué)陸全明教授關(guān)于混合模擬的討論.