強爆炸早期火球光輻射能譜的數值計算*

高銀軍,閆 凱,田 宙,劉 峰

(西北核技術研究所,陜西 西安710024)

強爆炸早期火球光輻射能譜的數值計算*

高銀軍,閆 凱,田 宙,劉 峰

(西北核技術研究所,陜西 西安710024)

基于強爆炸火球光輻射的多群輻射流體力學方法,采用算子分裂方法將方程組分裂為對流項和剛性源項,其中源項部分根據方程形式,進一步分裂為各群內的單獨求解。數值計算表明:該方法克服了直接求解過程中輻射與流體耦合所帶來的強不穩(wěn)定性,時間步長大幅提高,給出的火球光輻射能譜特征與已有規(guī)律一致?？蔀槎糠治龉廨椛淠茏V特征提供有效手段。

爆炸力學;算子分裂算法;光輻射能譜;強爆炸火球

在強爆炸過程中,由于強爆炸所釋放的巨大能量迅速加熱周圍冷空氣,從而形成高溫高壓的火球?；鹎蛟诎l(fā)展過程中,不斷向外輻射光和熱,稱為光輻射(也稱為熱輻射)。光輻射是強爆炸的重要毀傷效應之一,研究其能譜特征對于深入分析光輻射與物質相互作用、評估相應的毀傷效應具有重要意義[1]。

強爆炸火球光輻射是一個涵蓋近紫外到可見光直至遠紅外的寬譜帶,同時伴隨火球的發(fā)展變化過程呈現不同的特征。早期的研究是基于實驗測量,結合相應的理論分析,給出了光輻射各譜段較為接近的經驗公式[1-2]。隨著 相關理 論及數 值模擬 研究的 發(fā) 展,特 別 是 對 火 球 發(fā) 展 過 程 的 深 入 分 析,較 普 遍 認為采用灰體近 似下的 輻射流 體力學 方法可 以較好 的描述 火球發(fā) 展過程[3-7],并 給出了 光輻射 總能 量 的 時空分布規(guī)律[8-11]。但是灰體近似在處 理 光 輻 射 輸 運 過 程 中 ,假 定 輻 射 系 數 (空 氣 吸 收 系 數 )與 光 子 能 量無關,采用全譜范圍內的平均計算,因而無法給出光輻射能譜特征。

為能夠給出光輻射能譜,需采用多群方法,即在不同能譜范圍內求解光輻射輸運過程,分群越多,計算得到光輻射能譜特征越精細[12]。顯 然,這 種 處 理 方 法 的 代 價 使 得 輻 射 方 程 的 計 算 更 為 復 雜 ,需 要 輻射參數(分群吸收系數)更多,與流體耦合后求解,更容易由于方程的強非線性、強剛性[2](尤其是輻流耦合項)而 出 現 計 算 不 穩(wěn) 定(除 非 將 時 間 步 長 限 制 非 常 小)。 國 外 有 關 文 獻[4,12]僅 給 出 了 部 分 結 果 ,但 并 未涉及方程求解 中的具 體 處 理 過 程。 國 內 光 輻 射 能 譜 方 面 的 計 算 較 少,早 期 計 算[13]中 利 用 小 步 長 (約10-15s)直接求解的方法,給出了光輻射較短時間內的能量分布情況。

本文中基于輻射流體力學方法,采用多群方法求解光輻射輸運過程。方程求解過程中,為克服方程的強非線性、輻流耦合項的強剛性可能產生的計算不穩(wěn)定問題,通過算子分裂方法,將多群方程逐步分裂為對流項和源項單獨求解[14-15],其中源項部分 的求解需要在各群內進一步分裂,使之具 有單群 方程的形式。在這樣的處理方法下,多群輻流方程的求解比直接求解更穩(wěn)定,即使在大時間步長(約10-11s)下也能得到較好的結果,有效提高了計算效率,為大規(guī)模問題求解提供了條件。利用該計算方法,計算得到了強爆炸早期光輻射能譜特征分布,與經驗關系給出的規(guī)律較為一致,為定量分析光輻射能譜特征提供了有效手段。

1 計算方程及參數

在局域熱動力平衡(LTE)假定下,描述強爆炸火球光輻射發(fā)展的多群輻射流體力學方程如下:

式(1)~(3)分 別 為 質 量 守 恒 、動 量 守 恒 和 能 量 守 恒 方 程 。 式 中 :ρ為 空 氣 密 度 ,v為 空 氣 速 度 ,p為 空 氣壓 力 ,e為 空 氣 內 能 。為 各 群 內 光 輻 射 輸 運 的 總 動 量為各 群 內 光 輻 射 輸 運 的 總 能 量,上 標 “0”和 “1”表 示 輻 射 輸 運 方 程 的0階矩 和1階 矩 方 程,下 標g=1,2,3… 即 分 群 數 目 ,分別為第g群輻射能密度、輻射能流、輻射壓力張量及光輻射輸運能量和光輻射輸運動量。具體計算如下:

式 中 :c為 光 速 ,κg為 第g 群 吸 收 系 數 ,γ為 光 子 頻 率 。 由 于 各 群 光 子 能 量 具 有 上 限 和 下 限 ,因 此 采 用γg、γg-1分 別 標 記 第g群 光 子 頻 率 的 上 限 和 下 限 。E(γ),F(γ)及P(γ)分 別 為 輻 射 能 密 度 、輻 射 能 流 及輻 射 壓 強 張 量 函 數,與關 系 為為 黑 體 輻 射 譜 分 布 ,h為 普 朗 克 常 數 ,k為 波 爾 茲 曼 常 數 ,Bg為 黑 體 輻 射 下 第g 群 輻 射 能 密 度 :

圖1 吸收系數與光子能量(分群)和溫度關系Fig.1 Absorption coefficient of air varied with the temperature and photon energy

吸收系數采用文獻[17-18]給出的21群 Rosseland分群吸收系數κg。多群方法下,吸收系數為空氣溫度、密度以及光子分群能量的函數。數值求解中,通過雙線性插值計算相應的溫度、密度下,空氣的分群吸收系數。圖1給出空氣密度為,分群吸收系數與光子能量(分群)和溫度的關系。

對 于 各 群 內 分 群 輻 射 壓 強 張 量Pg,采 用 最 大 熵 變 Eddington 因 子 近進 行 計 算 ,空 氣 的 狀 態(tài)方 程 采 用 實 際 空 氣 狀 態(tài) 方。 光 輻 射 按 照 對 應 的 光 子 能 量pe分 為 21,如 表 1 所 示 。

表1 光子分群能量Table 1 Photon energy of each group

2 數值求解

2.1 分裂方法過程

一維球對稱下多群輻流方程組展開后也可寫成如下矩陣形式:

各項代表的物理量如下:

式中:Ψ為對流項,Φ為輻射與流體耦合項,即剛性源項部分。利用算子分裂方法,將方程分2步求解:

第1步求解對流項,采用有限體積法,構造五階 WENO 格式,數值通量的計算采用局部 Lax-Friedrichs方法,由于不 含剛性 源項,因而時間步長可以提高到10-11s,相比于直接法求解[13],極 大提高了計算效率。第2步求解源項,其初始 時刻(t=t0)的值 為第1 步對流 項方程 的解f(1)*。 源 項 求 解 較為復雜,首先根據方程特征,寫為:

分裂成為這樣的步驟進行求解,方程形式大為簡化,其過程本身也具有自身的物理意義:每步求解認為流體僅與該群光子進行能量和動量交換,二者組成的體系中,總能量和總動量守恒。這樣的假設,在方程分裂的數學處理過程中是嚴格滿足的:因為其他群內輻射能密度與輻射能流隨時間的微分都等于0。進一步的求解,則轉化為常微分方程組的求解,可以通過多種方法實現快速、高精度求解。

2.2 初始條件和邊界條件

求解初始條件假定爆炸總能量集中于等壓火球內,輻射能與空氣內能總合等于爆炸總能量,計算邊界條件采用對稱邊界。光輻射分群能量的初始分布通過Bg給出,分群輻射能流為0。

在數值求解中,計算區(qū)域的邊界處,假定物理量都處于未擾動狀態(tài),即認為:流體速度以及輻射能流均為0,而其他狀態(tài)參量取初始值。

3 計算結果與分析

取當量為1 kt,高度在海平面(h=0 km),求解相應的強爆炸火球光輻射輸運過程,空氣初始狀態(tài)參量如表2所示。

表2 不同高度空氣初始狀態(tài)參數Table 2 Air state parameters at different altitude

強爆炸火球發(fā)展過程中,沖擊波陣面和輻射波陣面是一個描述火球發(fā)展的重要參量。利用上述算子分裂方法計算得到的強爆炸火球陣面走時,如圖2所示,符合火球發(fā)展的物理過程,與 H.L.Brode[4]計算的結果符合較好,說明該方法在處理過程中是穩(wěn)定可靠的。

圖2 算子分裂法求解多群輻流方程組得到的早期火球陣面走時Fig.2 Calculational result of 1 kt nuclear fireball front by splitting method

多群方法的運用,能夠給出火球光輻射在特定時刻向外輻射的光輻射能譜特征。圖3所示為t=0.023 s時 刻 火 球 光 輻 射 分 群 (21 群)能 量 分布,I為光輻射強度。從圖3中可以看出,火球光輻射能量集中于第2~8群,對應光子波長在0.2 ~2μm。

根據火球發(fā)展過程,在光輻射強度第1個極大 值 后,火 球 有 效 溫 度 從 約 20 000 K 降 低 至3 000 K,在 光 輻 射 強 度 第2個 極 大 值 時,火 球 有效 溫 度 略 低 于10 000 K。 在 這 個 溫 度 范 圍 內,火球光輻射能譜大部分能量都集中與紫外(0.22~ 0.36μm)、可見(0.36~0.64μm)和紅外(0.64~ 4.5μm)部分,與圖3所示基本一致。

為具體分析火球光輻射能譜特征,利用文獻[1-2]給出的火球有效溫度走時關系:

圖3 光輻射分群能譜Fig.3 Fireball radiation energy of each group

式 中 :Te為 火 球 有 效 溫 度 ,σ為 斯 特 潘-玻 爾 茲 曼 常 數 ,Φe為 與 火 球 輻 射 功 率 有 關 的 函 數 表 達 式 ,r為 火球半徑,各個參數可以通過已有規(guī)律計算給出。以t=0.01、0.02 s時刻為例,利用上述理論方法計算得到該爆炸條件下火球有效溫度約為3 500K 和7 453 K,據此給出0.01~2.5μm 波長范圍內光輻射強度隨波長關系與對應有效溫度下的黑體譜對比分析,如圖4所示,圖中實線為本文中方法的計算結果,虛線為由文獻理論計算的黑體譜分布。

已 有 研 究 結 果 表 明[1-2],在 火 球 光 輻 射 強 度 第 1 個 極 大 值 后 的 整 個 發(fā) 光 階 段,波 長 在 (0.4~ 0.6μm)范圍內,可以近似看做黑體,與圖4給出的結果在變化規(guī)律上基本一致,數值大小上的差異,與所采用的空氣吸收系數等有關,需要進一步工作進行細致改進。

圖4 不同時刻光輻射強度與波長變化關系Fig.4 Relation between intensity of fireball radiation and wavelength

4 結 論

(1)基于強爆炸火球的多群輻射流體力學模型,采用算子分裂方法對方程組進行數值求解。利用空氣分群(21群)吸收系數,計算給出了1 k T 當量下火球光輻射能譜特征。數值計算結果驗證了光輻射能譜主要集中在0.2~2μm(紫外、可見到紅外波段),與已有結果和經驗規(guī)律符合的較為一致;

(2)分裂求解的處理方法,一方面克服了直接求解過程中輻射與流體耦合可能帶來的強不穩(wěn)定性,另一方面擴大了時間步長,提高了計算效率,為類似方程的數值求解提供了一定借鑒。

[1]喬 登 江 .核 爆 炸 物 理 概 論 [M].北 京:國 防 工 業(yè) 出 版 社,2003:225-262.

[2]屠 琴 芬 .核 爆 炸 火 球 的 輻 射 流 體 力 學 計 算 中 的 幾 個 問 題[R].西 安:西 北 核 技 術 研 究 所 ,1986.

[3]Pomraning G C.The equations of radiation hydrodynamics[J].Astrophysical Radiation Hydrodynamics,1986 (188):45-69.

[4]Brode H L.Fireball phenomenology[R].AD0612197,1964.

[5]Crowley B K,Glenn H D,Marks R E.An analysis of marvel:A nuclear shock-tube experiment[J].Journal of Geophysical Research,1971,76(14):3356-3374.

[6]Marrs R E,Moss W C,Whitlock B.Thermal radiation from nuclear detonations in urban environments[R]. UCRL-TR-231593,2007.

[7]Lowrie R B,Edwards J D.Radiative shock solutions with grey nonequilibrium diffusion[J].Shock Waves,2008, 18(2):129-143.

[8]陳 健 華 ,王 心 正,謝 龍 生 ,等.均 勻 大 氣 中 的 強 爆 炸 一 維 輻 射 流 體 力 學 數 值 解[J].爆 炸 與 沖 擊 ,1981,1(2):37-49. Chen Jian-hua,Wang Xin-zheng,Xie Long-sheng,et al.A one-dimensional radiation hydrodynamic numerical solution for a strong explosion in uniform atmosphere[J].Explosion and Shock Waves,1981,1(2):37-49.

[9]田 宙 , 喬 登 江 ,郭 永 輝 . 不 同 當 量 強 爆 炸 早 期 火 球 現 象 的 數 值 模 擬[J].爆 炸 與 沖 擊 ,2009,29(4):408-412.Tian Zhou,Qiao Deng-jiang,Guo Yong-hui.Numerical simulation on early fireball phenomenology of strong explosions for different yields[J].Explosion and Shock Waves,2009,29(4):408-412.

[10]田 宙,喬 登江 ,郭永 輝.不同高 度強爆炸早 期火球數值 研究[J].兵 工學 報,2009,30(8):1078-1083. Tian Zhou,Qiao Deng-jiang,Guo Yong-hui.Numerical investigation of early fireball of strong explosion for different altitudes[J].Acta Armamentarii,2009,30(8):1078-1083.

[11]田 宙,喬 登江 ,郭永 輝.強爆炸 早期火球現 象的一維數 值研究[J].計算 物理,2010,27(1):9-14. Tian Zhou,Qiao Deng-jiang,Guo Yong-hui.A one dimensional numerical study on early fireball in strong explosion[J].Chinese Journal of Computational Physics,2010,27(1):9-14.

[12]Symbalisty E M D,Zinn J,Whitaker R W.Radflo physics and algorithms[R].LA-12988-MS,1995.

[13]高 銀軍,田宙 ,劉峰 ,等 .強爆 炸早期火球 光輻射能譜 的分群計算[J].四川兵 工學 報,2011,32(3):21-24. Gao Yin-jun,Tian Zhou,Liu Feng,et al.Calculation of energy spectrum of early fireball radiation in strong explosion with multi-group method[J].Journal of Sichuan Ordnance,2011,32(3):21-24.

[14]閆 凱.二 維輻 射流體 動力學方程 組的數值求 解[D].西安 :西北核技 術研究所,2011.

[15]閆 凱.二 維輻 射流體 動力學方程 組的數值求 解[C]∥第六屆 全國青年計 算物理學術 會議.太原 ,2011.

[16]Minerbo G N.Maximum entropy eddington factors[J].Journal of Quantitative Spectroscopy and Radiative Transfer,1977,20(6):541-545.

[17]王 文高,張建 泉.物 質輻射不透 明性Ⅰ[R].西 安:西北核 技術研究所 ,1978.

[18]王 文高,張建 泉.物 質輻射不透 明性Ⅱ[R].西 安:西北核 技術研究所 ,1978.

Numerical calculation of early fireball radiation spectrum in strong explosion

Gao Yin-jun,Yan Kai,Tian Zhou,Liu Feng
(Northwest Institute of Nuclear Technology,Xi’an 710024,Shaanxi,China)

On the basis of multi-group radiation hydrodynamics method of fireball radiation in strong explosion,operator splitting method is used to split the equations into convection items and source items,which are split into radiation groups due to the equation formation and solved individually.Numerical calculations show that the method used here overcomes strong instability when solving the equations directly because of the coupling items between radiation and fluid.In the meantime,the time step in the calculation is increased obviously.Fireball radiation spectrum is obtained in fine accordance with the result in the literature.

mechanics of explosion;splitting method;radiation spectrum;strong explosion fireball

O381國標學科代碼:13035

:A

10.11883/1001-1455-(2015)03-0289-07

(責任編輯 王易難)

2013-03-13;

2013-05-31

高銀 軍(1983— ),男,碩士研 究生,助理 研究員,gyj@mail.ustc.edu.cn。