低空強爆炸中火球的一維數(shù)值模擬研究

閆 凱，劉 鈺，田 宙（西北核技術研究所，陜西西安 710024）

低空強爆炸中火球的一維數(shù)值模擬研究

閆 凱，劉 鈺，田 宙
（西北核技術研究所，陜西西安 710024）

摘要：對海平面環(huán)境下當量為1Mt TNT的強爆炸進行一維輻射流體力學數(shù)值計算，得到了較為完整的火球擴張過程數(shù)值模擬結果。數(shù)值給出了火球不同發(fā)展階段物質(zhì)參數(shù)和輻射參數(shù)的空間分布，該結果與實際強爆炸中一些經(jīng)驗公式的計算結果吻合較好。給出了火球陣面和輻射功率隨時間的變化關系，較合理地解釋了低空強爆炸火球二次極大現(xiàn)象。

關鍵詞：火球；輻射流體力學；數(shù)值模擬；二次極大

Numerical Simulation of Fireball in Strong Explosion at Low Altitude

YAN Kai，LIU Yu，TIAN Zhou
（Northwest Institute of Nuclear Technology，Xi’an710024，China）

火球現(xiàn)象是強爆炸早期最顯著的特征之一。根據(jù)維恩位移定律，強爆炸瞬間輻射波以高頻波為主，約75%～85%的爆炸能量為X射線［1］。在低空環(huán)境下，這種X射線的平均自由程只有m量級。X射線在傳播過程中加熱周圍空氣，形成一團高溫高壓的氣體，即火球。火球迅速向外膨脹，壓縮周圍空氣形成沖擊波，同時向外發(fā)出光輻射。

低空強爆炸火球的重要特征是二次極大現(xiàn)象，即火球光輻射會出現(xiàn)兩個脈沖。該現(xiàn)象的形成和空氣在被加熱過程中引起的光學性質(zhì)變化有關［2－4］，文獻［2］指出，隨著沖擊波卷入的空氣層增厚，空氣的不透明度增加，火球的亮度逐漸下降，達到極小值。但當火球陣面溫度降至2 000K以下，波后空氣對輻射的吸收越來越弱，逐漸變?yōu)橥该鳎锩娴幕鹎蛞环矫嬷饾u冷卻，另一方面卻越來越看得清楚。目前國內(nèi)尚未見到對該現(xiàn)象的數(shù)值模擬結果。

數(shù)值模擬低空強爆炸火球現(xiàn)象的主要困難有兩點。首先是輻射輸運模型問題，由于低空環(huán)境下溫度對空氣不透明度有重要的影響，爆炸初期火球邊緣附近的空氣不透明度變化極為劇烈，因此輻射輸運模型需能同時描述光學厚和光學薄介質(zhì)中的光子輸運過程。其次是剛性問題，輻射流體方程組涉及3種重要的動力學尺度：物質(zhì)波波速（聲速）、輻射波波速（光速）及源項［5］。當輻射與物質(zhì)相互作用較強時，源項代表最快的尺度；當輻射與物質(zhì)相互作用較弱時，源項代表最慢的尺度。因此如何處理源項問題是提高程序計算效率的關鍵。針對上述兩點困難，本文采用一種最大熵變Eddington因子的M1近似來描述輻射場［6－7］，并設計一種分裂格式來處理源項問題［8－9］。通過對當量為1Mt TNT的強爆炸火球在海平面大氣中的數(shù)值模擬，計算火球參數(shù)在不同發(fā)展階段的時空分布，給出火球二次極大現(xiàn)象的數(shù)值解釋。

1　輻射流體方程

本文采用文獻［10］給出的一維球?qū)ΨQEuler型輻射流體動力學方程組。其中，流體力學方程忽略重力和黏性，輻射輸運方程采用灰體近似和局部熱力學平衡近似。由于爆炸初期火球的擴張速度很快，因此必須考慮相對論效應的影響。常用的做法是方程組采用互動坐標系或混合坐標系。本文考慮采用混合坐標系的形式，將相對論修正項當成源項的一部分，這樣輻射流體方程組可寫成雙曲守恒律形式，即：

其中

此時輻射輸運方程是不封閉的，本文引入

Levermore等［6－7］構造的變Eddington因子χ，令：

其中：r為空間坐標；t為時間坐標；p、ρ、u、T分別為氣體的壓力、密度、速度和溫度；e、eI分別為單位體積氣體的總能量和內(nèi)能為單位體積氣體的輻射能密度；為輻射能流；為輻射壓力；c為光速；φ為平衡輻射能密度

2　計算方法

輻射流體方程組中源項的尺度對方程組的計算效率影響很大，本文應用以下算子分裂方法將耦合項從方程組中分裂開來：

其中

式（4）為雙曲守恒律方程組，式（5）為常微分方程組。對于雙曲守恒律方程組，本文采用有限體積5階WENO格式進行數(shù)值求解；由于常微分方程組具有強剛性的特點，一般的顯式迭代法時間步長必須非常苛刻才能得到收斂解。本文設計一種高效的溫度迭代方法，數(shù)值計算的具體過程可參考文獻［9－10］。

3　數(shù)值計算結果及分析

計算選取當量為1Mt TNT的強爆炸，火球半徑取為0．75m。采用文獻［11－12］的做法，數(shù)值計算中將源區(qū)物質(zhì)等價為空氣，源區(qū)總質(zhì)量為1 000kg。強爆炸源區(qū)設為等壓球，強爆炸能量均勻分布在等壓球內(nèi)，等壓球外為實際空氣。

初始時刻，等壓球處于熱力學平衡狀態(tài)，等壓球內(nèi)輻射能與氣體能量之和等于爆炸總能量。等壓球內(nèi)外的密度為大氣密度，速度為0，溫度、壓力、內(nèi)能利用狀態(tài)方程給出。輻射能密度利用溫度求出，輻射能流為0，吸收系數(shù)的計算采用灰體近似法。

本文計算采用一維球坐標系，其中內(nèi)邊界采用對稱邊界條件。對流體部分，壓力、密度、能量對稱相同，速度對稱相反；對輻射輸運部分，輻射能、輻射壓力對稱相同，輻射能流對稱相反。外邊界則采用透射邊界條件。

從時間上來說，一個完整的火球發(fā)展過程應包括X射線火球階段、輻射擴張階段、過渡階段及沖擊波擴張階段（也有文獻將過渡階段劃分到輻射擴張階段）。這4個階段是以火球擴張的主要作用力的不同而劃分的。在火球發(fā)展的早期階段，約在前1～2μs，稱為X射線火球階段。X射線火球階段尚難以用灰體近似法來描述，必須采用多群方法。由于本文僅采用灰體近似法，基于這種考慮，本文未計算X射線火球階段，而是將X射線加熱空氣形成的火球作為初始條件。本文數(shù)值模擬的初始條件可看作輻射擴張階段的開始。由于強爆炸火球參數(shù)的實際測量較困難，國內(nèi)、外公開發(fā)表的文獻中給出的結果大多數(shù)是定性的。本文將火球發(fā)展不同階段的計算結果與文獻［2］進行定性對比。

3．1 輻射擴張階段

圖1　 輻射擴張階段火球的物質(zhì)參數(shù)Fig．1 Fireball’s material parameters in radiation transport period

圖1為輻射擴張階段火球的物質(zhì)參數(shù)。在該階段，火球擴張主要是由輻射輸運造成的，其特點是火球以輻射熱波的形式向外傳播。文獻［2］指出該階段波陣面上溫度、壓力和速度均有躍變，且波后溫度大體是常數(shù)，即等溫球。圖1中各物理參數(shù)在5μs以前的空間分布基本反映了這一特點。圖2為火球輻射參數(shù)的時空分布。可看出，火球陣面附近輻射能密度急劇降低，這是由于爆炸初期源區(qū)的高頻射線不斷被周圍空氣吸收并加熱空氣的結果。輻射能流在火球陣面附近有突變，這是由火球吸收系數(shù)的變化引起的。雖然火球內(nèi)部的溫度很高，但是由于源區(qū)內(nèi)的密度較大，火球內(nèi)部的吸收系數(shù)（灰體近似條件下）仍較大，因此火球內(nèi)部并不透明；在火球陣面附近，氣體密度較低，火球在該區(qū)域吸收系數(shù)較小；但當火球溫度處于20 000～100 000K之間時，火球的吸收系數(shù)又會突然增大。隨著火球的擴張，火球溫度和輻射能密度迅速降低，其中輻射能密度下降的更快。該物理量實際反映火球的輻射擴張能力，它的減小預示著流體力學過程的作用開始增加，火球開始進入過渡階段。

3．2 過渡階段

隨著時間的推移，流體力學過程的作用不能再忽略。表現(xiàn)在火球陣面的溫度、壓力不斷下降，火球陣面逐漸向強沖擊波過渡，即所謂過渡階段。該階段火球的增長是輻射擴展和流體力學過程的共同作用。文獻［2］將火球發(fā)展過程分為輻射擴張階段和沖擊波擴張階段，本文采用文獻［12］的觀點，認為在輻射擴張階段和沖擊波擴張階段之間存在一流體和輻射共同作用的階段，即過渡階段，該階段對應于文獻［2］中輻射擴張階段后期。

過渡階段流體力學過程不能忽略，表現(xiàn)在火球陣面附近的密度逐漸有躍變，且密度壓縮比逐漸趨向經(jīng)典激波的極限，同時彈殼沖擊波（又稱等溫激波）開始形成。圖3為過渡階段火球的物質(zhì)參數(shù)。從圖3中的溫度和壓力曲線上可明顯看出彈殼沖擊波的形成過程。在彈殼沖擊波的陣面上各物理量均有突變，波后的溫度、壓力基本是均勻的，而且在火球內(nèi)部會產(chǎn)生多個二次激波，圖3中30、100、220μs時火球中心溫度和壓力均有突變，二次激波對火球中心物質(zhì)參數(shù)的影響較大，但對輻射參數(shù)的影響并不大，如圖4所示。隨著時間的推移，二次激波對火球陣面的影響越來越小，在220μs時對火球陣面的影響已非常微弱。根據(jù)文獻［2］，形成沖擊波時相應的波陣面速度約為62km／s，陣面溫度約為180 000K，大致可認為220μs左右時火球開始進入沖擊波擴張階段。

圖2　 輻射擴張階段火球的輻射參數(shù)Fig．2 Fireball’s radiation parameters in radiation transport period

圖3　 過渡階段火球的物質(zhì)參數(shù)Fig．3 Fireball’s material parameters in transition period

圖4　 過渡階段火球的輻射參數(shù)Fig．4 Fireball’s radiation parameters in transition period

3．3 沖擊波擴張階段

圖5為沖擊波擴張階段火球的物質(zhì)參數(shù)。當?shù)葴丶げㄚs上輻射波后，即進入火球發(fā)展的強沖擊波擴張階段。當沖擊波超過輻射波后，在溫度空間剖面上出現(xiàn)由沖擊波陣面溫度Ts到輻射波陣面溫度TI的階梯，形成沖擊波擴張階段火球特有的結構，即火球內(nèi)部在輻射作用下為等溫球，火球邊緣基本對應沖擊波波峰位置。隨著火球的進一步擴張，沖擊波陣面的溫度壓力不斷下降，隨著沖擊波波速不斷下降，約在10ms時沖擊波陣面密度壓縮比達到頂峰，隨后開始下降。

圖6為沖擊波擴張階段火球的輻射參數(shù)。對比圖5、6可發(fā)現(xiàn)，此階段輻射能密度的陣面已落后于沖擊波陣面，其陣面大致與等溫球陣面一致。雖然此時輻射對沖擊波傳播的影響已減弱（表現(xiàn)在輻射能流迅速減小），但對火球現(xiàn)象仍要用輻射流體力學來描述。這是因為火球內(nèi)部溫度仍較高，吸收系數(shù)仍很大，因此輻射仍是傳熱的主要方式，火球內(nèi)部基本處于輻射平衡狀態(tài)。在等溫球內(nèi)部，密度、速度均是均勻分布的，且隨著時間單調(diào)下降。火球表面與等溫球表面殼層間，壓力和密度明顯上升，溫度則明顯下降。

圖5　 沖擊波擴張階段火球的物質(zhì)參數(shù)Fig．5 Fireball’s material parameters in shock transport period

圖6　 沖擊波擴張階段火球的輻射參數(shù)Fig．6 Fireball’s radiation parameters in shock transport period

3．4 火球沖擊波陣面與文獻［2］結果的對比

根據(jù)實測沖擊波參數(shù)的分析，文獻［2］提出應用以下經(jīng)驗公式描述強爆炸沖擊波的傳播規(guī)律。

式中：Q為爆炸總當量，kt；Rd為沖擊波波峰距爆心的距離，m；p0為壓力的峰值，105Pa。

圖7為沖擊波擴張階段后（此時沖擊波波峰距爆心約50m）的壓力峰值與經(jīng)驗公式的對比，可看出，兩者差異較小，說明本文的算法能較好地描述強爆炸沖擊波。

圖7　 沖擊波壓力峰值隨沖擊波半徑的變化Fig．7 Peak pressure of shock wave versus radius

同時文獻［2］提到：在沖擊波擴張階段，如果對模擬解不同時刻用不同γ的理想氣體來逼近真實大氣，那么沖擊波陣面R（t）可寫為：

圖8為本文計算得到的沖擊波陣面與經(jīng)驗公式計算結果的對比。可看出，在爆炸早期，本文計算結果較經(jīng)驗公式計算的結果略高，這是由于本文采用的空氣吸收系數(shù)與實際有差異引起的；在沖擊波陣面擴張至400m后，本文計算結果較經(jīng)驗公式計算的結果略低，這是由于隨著計算區(qū)域的擴大網(wǎng)格尺寸變大（導致沖擊波波峰不夠鋒銳，沖擊波速度變慢）引起的。總體而言，兩者計算結果差異不大，說明本文采用的物理模型可描述火球沖擊波的擴張過程。

圖8　 沖擊波陣面隨時間的變化Fig．8 Shock wave front versus time

3．5 火球的光輻射參數(shù)

光輻射的源是火球，因此火球的輻射功率與遠方各點接收到的光輻射的能量有關。輻射功率可表示為：

式中：4πR2f為火球的面積；B為全波亮度，在灰體近似的條件下，可將輻射能流近似認為全波亮度，即：

圖9　 輻射功率隨時間的變化Fig．9 Radiation power versus time

圖9為本文計算的輻射功率與文獻［11］計算結果的對比。文獻［13］主要介紹了RADFLO程序?qū)姳ɑ鹎虻囊恍?shù)值模擬結果，該程序采用多群方法，輻射方面只考慮輻射能流對物質(zhì)內(nèi)能的影響。本文計算得到的第一極大，第一極小和第二極大時間分別為0．002、0．095和0．84s。文獻［2］給出λ≈0．65μm的紅光亮度的第一極小和第二極大到來時間的經(jīng)驗公式為：

對于當量為1Mt的強爆炸，計算得到紅光亮度的第一極小和第二極大時間約為0．08s 和0．62s。此結果較本文的計算結果略小，這是由于本文采用灰體近似法，導致吸收系數(shù)的計算不準確引起的。總體來看，本文的計算結果與文獻［2］中的經(jīng)驗公式及RADFLO的計算結果差異不大，驗證了本文方法的合理性。

由于物理模型、空氣狀態(tài)方程、吸收系數(shù)等處理的不同，本文與RADFLO程序得到的輻射功率有差異。此外，在火球第一極大前，兩者的結果均有震蕩，這是由于此時火球溫度的變化極為劇烈，而溫度又影響網(wǎng)格的光子平均自由程造成的，原則上網(wǎng)格尺寸要小于光子平均自由程才能得到合理的數(shù)值解，但由于計算效率的原因網(wǎng)格尺寸不可能滿足該條件。在火球第一極大后，隨著火球陣面溫度的降低，輻射功率的變化較為光滑。

圖10為火球、沖擊波及等溫球半徑隨時間的變化。可看出，沖擊波形成后，火球表面以沖擊波速度向外擴張，而等溫球界面的擴張速度明顯減慢，約在2s等溫球消失。在0．13s左右，沖擊波脫離火球，此時對應的沖擊波半徑為500m左右，文獻［2］給出的沖擊波脫離火球的半徑和時間的估算公式為：

圖10　 沖擊波陣面、火球陣面及等溫球陣面隨時間的變化Fig．10 Shock wave front，fireball front and isothermal sphere front versus time

對于當量為1Mt的強爆炸，可計算出脫離半徑約為524m，脫離時間約為0．12s，此結果與本文的計算結果基本符合。該時間基本對應火球輻射功率極小值時間。

4　結論

采用Euler坐標系下的能描述光學厚和光學薄介質(zhì)中光子輸運的輻射輸運方程，連同流體力學方程構成輻射流體力學方程組。數(shù)值計算給出了強爆炸火球現(xiàn)象中火球陣面、彈殼沖擊波的形成發(fā)展過程及其物理參數(shù)的分布規(guī)律。計算結果與文獻結果符合較好，證明了本文所采用物理模型、數(shù)值方法及計算程序的合理性。最后給出了低空強爆炸火球光輻射參數(shù)的長時間計算結果，計算結果表明：強爆炸初期，由于源區(qū)物質(zhì)不透明度較大，X射線基本被吸收，輻射功率較小；隨著強沖擊波的形成，火球陣面附近的不透明度逐漸減小，輻射功率達到第一極大點；隨著沖擊波卷入空氣的增厚，由于空氣的不透明度（灰體近似）在20 000～100 000K時達到最大，因此空氣又變得不透明，輻射功率達到了第一極小點；沖擊波脫離火球后，火球陣面本身發(fā)光已微弱，且對波后輻射的吸收也逐漸微弱，火球外面有透明層出現(xiàn)，此時里面的火球逐漸露出，輻射功率達到第二極大點。

參考文獻：

［1］ 歐陽建明，馬燕云，邵福球，等．高空核爆炸X射線電離的時空分布數(shù)值模擬［J］．物理學報，2012，61（24）：242801．OUYANG Jianming，MA Yanyun，SHAO Fuqiu，et al．Numerical simulation of X－ray ionization and atmospheric temporal evolutions with high－altitude nuclear explosions［J］．Acta Phys Sin，2012，61（24）：242801（in Chinese）．

［2］ 喬登江．核爆炸物理概論［M］．北京：國防工業(yè)出版，2003：169－262．

［3］ 吳健輝．核爆炸光輻射特性及探測技術的理論與實驗研究［D］．武漢：華中科技大學，2009．

［4］ GOLDBLAT J，COX D．Nuclear weapon tests：Prohibition or limitation［M］．ENG：Oxford University Press，1988：10－18．

［5］ SEKORA M，STONE J．A higher order Go－dunov method for radiation hydrodynamics：Ra－Computational Physics，2013，30（3）：379－386 diation subsystem［J］．Communications in Ap－（in Chinese）．plied Mathematics and Computational Science，

［10］SEKORA M．Algorithms for hyperbolic balance 2009，4（1）：135－152．laws with multiscale behavior：Applications in

［6］ LEVERMORE C D．Relating Eddington factorsradiation hydrodynamics［D］．United States：to flux limiters［J］．Journal of Quantitative Spec－Princeton University，2010．troscopy and Radiative Transfer，1984，31（2）：

［11］田宙，喬登江，郭永輝．不同高度強爆炸早期火149－160．球數(shù)值研究［J］．兵工學報，2009，30（8）：1 078－

［7］ ANILE A M，PENNISI S，SAMMARTINO M．1 084．A thermodynamical approach to Eddington fac－TIAN Zhou，QIAO Dengjiang，GUO Yonghui．tors［J］．Journal of Mathematical Physics，1991，Numerical investigation of early fireball of strong 32（2）：544－555．explosion for different altitudes［J］．Acta Arma－

［8］ LOWRIE R B，MOREL J E．Issues with high－mentarii，2009，30（8）：1 078－1 084（in Chiresolution Godunov methods for radiation hydro－nese）．dynamics［J］．Journal of Quantitative Spectrosco－

［12］BRODE H L．Thermal radiation phenomena，py and Radiative Transfer，2001，69（4）：475－VolumeⅤ：Radiation hydrodynamics of high 489．temperature air，AD672837［R］．United States：

［9］ 閆凱，李若，田宙，等．強爆炸火球數(shù)值模擬中的Lockheed Missiles and Space Company Sunny－算子分裂方法［J］．計算物理，2013，30（3）：379－vale，1967．386．

［13］EUGENE M D，JOHN Z，RODNEY W．RADYAN Kai，LI Ruo，TIAN Zhou，et al．An oper－FLO physics and algorithms，LA－12988－MS［R］．a(chǎn)tor splitting method for numerical simulationofUnited States：Los Alamos National Laboratory，strong explosion fireball［J］．Chinese Journal of1995．

作者簡介：閆 凱（1984—），男，河南民權人，助理研究員，碩士，從事輻射流體力學及磁流體力學數(shù)值模擬研究Abstract：A one－dimensional radiation hydrodynamics numerical computation was presented for a strong explosion at sea level，of which the equivalent is 1 Mt TNT．The spatial distributions of material and radiation parameters in strong explosion，which reflect a integrate process of fireball evolution，were given in terms of the simulation．The calculation results agree well with those of some empirical formulas．The results of optical observables－radiation power versus time and fireball radius versus time were described，and the second maximum phenomenon of the fireball at low altitude was given．Key words：fireball；radiation hydrodynamics；numerical simulation；second maximum

基金項目：國家自然科學基金資助項目（91330205）

收稿日期：2014－04－18；修回日期：2014－07－24

doi：10．7538／yzk．2015．49．08．1345

文章編號：1000－6931（2015）08－1345－09

文獻標志碼：A

中圖分類號：O38；O242