空中強爆炸沖擊波地面反射規律數值模擬研究*

姚成寶，王宏亮，浦錫鋒，壽列楓，王智環

（1. 北京大學數學科學學院，北京 100871；2. 西北核技術研究院，陜西 西安 710024）

炸藥在距地面有限高度內爆炸時，沖擊波將會在地面、建筑等障礙物的表面發生反射和繞射現象，傳播規律變得更加復雜。相比于化學爆炸，強爆炸產生的初始壓力可達上萬吉帕，初始溫度高達數十萬攝氏度[1-2]，空氣在沖擊波作用下將發生劇烈的壓縮和溫升，熱力學性質發生顯著變化，需要建立高溫高壓下真實氣體的狀態方程。此外，由于強爆炸沖擊波的傳播距離遠，其在傳播過程中會受到空氣不均勻性的影響，空氣的初始壓強和密度越低，沖擊波的衰減越明顯[2-3]。空中強爆炸問題具有的上述特點，對數值計算的精度和效率提出了很高的要求。

空中強爆炸問題屬于可壓縮多介質大變形問題的范疇，一般會存在多種介質的相互作用。由于不同介質的狀態方程和初始條件通常存在較大的差異，會給數值計算特別是介質界面附近的數值模擬帶來很大的困難。目前針對空中爆炸沖擊波問題的主流數值方法是基于Euler 坐標系下的多介質流體數值方法，其包含兩個關鍵步驟：界面追蹤和介質間相互作用的處理。其中，常見的界面追蹤方法包括VOF (volume-of-fluid)方法[4]、Level Set 方法[5]和Front Tracking 方法[6]等。當確定介質界面的位置后，需進一步準確描述不同介質間的相互作用，常見的處理方法包括虛擬流體類方法(ghost fluid method)[7-8]和切割單元法(cut cell method)[9]。

在爆炸沖擊波地面反射規律研究方面，Crowl[10]給出了TNT 等化學炸藥在不同爆高下的地面反射沖擊波載荷分布規律，喬登江[1]和Glasstone 等[2]依據實驗結果分別給出了標準大氣狀態下1 kt TNT 當量的空中強爆炸沖擊波在大尺度范圍內的等超壓曲線。目前公開發表的文獻大多針對TNT 等化學爆炸產生的沖擊波在局部范圍內的傳播[11-15]，較少見到能夠完整給出空中強爆炸沖擊波在地面反射過程的大尺度范圍數值模擬。

針對空中強爆炸問題的強對流性，本文在前期工作[16]的基礎上進一步改進基于歐拉坐標系的多介質流體數值方法，并考慮高溫、高壓下的真實氣體狀態方程和空氣隨高度不均勻分布的影響[17]，結合網格自適應技術，計算1 kt TNT 當量空中強爆炸產生的沖擊波在不同爆高下的地面反射完整過程，得到地面上距爆心投影點5 km 范圍內的沖擊波載荷分布。通過和實測結果進行比對，驗證計算結果的正確性，并在此基礎上研究地面反射沖擊波載荷隨爆高的變化規律。

1 物理模型和數值方法

1.1 物理模型

根據空中強爆炸問題的對稱性，采用二維柱對稱計算模型，控制方程組為：

式中： ρ 為密度，u、v分別為r,z方向的速度，E為體積總能量，p為壓力。

由于強爆炸伴隨著強烈的放熱、電離等效應，爆炸產物迅速膨脹至氣體狀態，且周圍空氣在受到沖擊波、熱輻射等因素的綜合作用下發生劇烈的加熱和壓縮，物理狀態極其復雜。本文中采用基于Saha 電離平衡理論給出的、以數值表形式給出的狀態方程[1]來描述強爆炸產物的熱力學行為，并采用真實氣體狀態方程來描述空氣：

式中： γ 為絕熱指數，e為比內能。Symbalisty 等[18]給出了8 組參考密度下(10-6～10 kg/m3)的溫度與內能、溫度與絕熱指數之間關系的實測結果，如圖1 所示。通過對參考曲線進行雙對數插值來獲得各熱力學變量之間的關系：首先，根據已知的比內能e，對圖1(a)進行對數插值，得到相應的溫度T；然后，對圖1(b)進行插值，獲得相應的絕熱指數 γ -1 ；最后，將 γ -1 代入式(3)進行求解，最終完成整個狀態方程的完整計算。

圖 1 不同參考密度下真實氣體狀態方程的熱力學參數關系Fig. 1 Plots of equations of state for real gas at different reference densities

實際空氣的初始壓力和密度隨著高度的增加逐漸降低，且滿足下列關系[1,17]：

式 中：Tref=288.16 K，pref=1.01325×105Pa，ρref=1.225 kg/m3分 別 為 標 準 大 氣 溫 度、壓 強 和 密 度;μ=28.966 g/mol，g=9.8 kg/m2，R=8.31 m2·s2/K；z為垂直高度，b為擬合系數。

1.2 數值方法

采用基于歐拉坐標系、具有互不相溶特征的二維柱對稱多介質流體計算模型來進行數值離散。在任意時刻，整個計算區域 Ω 可以分成2 個子區域(t) 和(t) ，且滿足：

式中： Γ (t) 為2 種互不相溶流體之間的物質界面。

在前期工作[16]的基礎上，采用Level Set 方法捕捉爆炸產物與空氣的物質界面，根據Level Set 函數分片線性的特征清晰重構出具有分片線性的物質界面，并通過在界面的法向上求解多介質Riemann 問題的精確解來構造爆炸產物與空氣之間的數值通量。整個計算流程如下。

1.2.1 物質界面的追蹤和重構

Level Set 方法[5,19-20]把隨時間運動的物質界面 Γ (t) 定義為距離函數 φ (x,t) 的零等值面。利用特征線方法求解Level Set 函數的演化方程：

式中：u? 為由物質界面的運動速度延拓得到的速度場[21]。

根據特征線 方 法 的思想，只 要 知道tn+1時 刻 網 格頂點Xi上的流體在tn時 刻的位置x(tn+1) ，就可 以 得 到新的Level Set 函數在Xi上的值 φ (x,tn+1) 。對每個頂點Xi，利用特征線方法可得：

再利用代數平均，可得新的Level Set 函數：

利用顯式正系數格式[21]求解重新初始化方程：

以保持Level Set 函數的符號距離性質，其中， φ0為重新初始化前φ (x,t) 的值，ε 為一很小的正數（例如10-6）。顯示正系數格式的求解過程[21]比較繁鎖，本文不再說細展開。

當得到tn時刻的Level Set 函數 φn后，根據其分片線性的特征重構出界面單元內分片線性流形的物質界面，如圖2 所示，圖中nKi,n為物質界面的單位法向。假設單元Ki,n與單元Kj,n具有共邊Sij,n，且 包含物質界面ΓKi,n。ΓKi,n將Ki,n和Sij,n依次分為2 個部分Ki±,n和Δtn，此時該單元包含2 種流體，且每種流體的守恒量分別定義為：

圖 2 多介質流體計算模型示意圖Fig. 2 Schematic diagram for the multi-media fluid calculation model

1.2.2 單元邊界數值通量

單元邊界的數值通量是指單元邊界上同種流體之間的數值通量，且滿足：

本文中采用local Lax-Friedrich 數值通量函數：

式中： λ =max{λKi,n,λKj,n} ， λKi,n、 λKj,n分別為和的最大信號速度。

1.2.3 物質界面數值通量

物質界面數值通量是物質界面兩側不同流體之間的數值通量，且滿足：

式中： | ΓKi,n| 為 物質界面的長度；和分別為物質界面上的壓力和法向速度，且通過在物質界面法向上求解局部一維多介質Riemann 問題的精確解來獲得，具體求解過程見文獻[22]。

1.2.4 守恒量的更新

當計算得到Ki,n的單元邊界數值通量和物質界面數值通量后，可對該單元中每種流體的守恒量進行更新：

2 空中強爆炸沖擊波地面反射

2.1 空中強爆炸沖擊波反射數值模擬

利用數值模擬程序對當量為1 kt TNT、爆高為100 m 的空中強爆炸沖擊波傳播過程進行數值模擬。其中，強爆炸產物采用等溫等壓球模型，取爆炸總當量的85%作為爆炸產物的力學初始能量[1]，產物質量設為50 kg，初始半徑設為30 cm。空氣采用真實氣體狀態方程，取水平面處的初始壓力為101.3 kPa，初始密度為1.29 kg/m3，并按照式(4)和(5)對不同高度處空氣的初始壓力和密度進行設置。

計算區域在r=0 和z=0 設為固壁反射邊界條件，其余邊界設為無反射邊界條件。數值計算中采用h-網格自適應技術[23]，根據相鄰網格的壓力梯度來進行網格的局部加密或稀疏，得到典型時刻的沖擊波壓力等值線以及相應的網格自適應圖，如圖3 所示。

由圖3 可知，當空中強爆炸產生的沖擊波傳播到地面時，最早在爆心投影點發生正反射，然后以逐漸增大的入射角發生斜反射，產生雙激波結構的規則反射。隨著沖擊波向外繼續傳播，入射角不斷增大。當入射角增大到臨界角附近時將發生馬赫反射，此時反射波陣面和入射波陣面的交點離開地面，形成垂直于地面的馬赫桿。隨著沖擊波的進一步傳播，馬赫桿迅速增長。圖4 給出了數值計算得到的地面沖擊波載荷分布(峰值超壓和沖量)與實測結果[2]的對比情況，其中在距爆心投影點1 km 以內，數值計算得到的峰值超壓與實測結果符合較好，而在1 km 以外則比實測結果略低，這是由于隨著計算區域的尺度增大，網格自適應后的單元尺寸不夠精細，使得數值計算對峰值超壓的捕捉能力減弱導致。此外，數值計算得到的沖量比實測結果略大，但最大誤差不超過50%。

圖 3 典型時刻的沖擊波壓力等值線圖和網格自適應圖Fig. 3 Pressure contours and adaptive meshes at typical times

圖 4 地面不同距離處的沖擊波峰值超壓和沖量Fig. 4 Peak overpressures and impulses at different radii

2.2 不同爆炸高度的沖擊波反射規律

利用數值模擬程序進一步對當量為1 kt TNT，爆高H依次為20、50、100、200、400 和800 m 的空中強爆炸沖擊波傳播過程進行了數值模擬，得到了地面上不同距離處的沖擊波峰值超壓和沖量，如圖5 所示。

圖 5 不同爆高下的地面沖擊波峰值超壓和沖量Fig. 5 Peak overpressures and impulses at different heights of burst

圖5 的計算結果表明，對1 kt TNT 當量的空中強爆炸，在距爆心投影點300 m 范圍以內，地面上的峰值超壓和沖量隨著爆高的增大而逐漸減小。在距爆心投影點300 m 范圍以外，隨著距爆心投影點的距離的增大，不同爆高條件下地面反射沖擊波載荷的差距逐漸變小。其中，爆高約為200 m 的空中強爆炸產生的地面峰值超壓和沖量處于最大值，其他爆高條件下的沖擊波載荷隨著爆高的進一步增大或減小均呈現逐漸減弱的趨勢，其原因主要是由于空中強爆炸地面反射而引發的復雜沖擊波傳播引起的。當空中強爆炸形成的沖擊波達到地面時，隨著入射角的增大，會依次在地面形成規則反射和馬赫反射等復雜的波系結構。其中，在規則反射區，反射沖擊波載荷由入射波強度和入射角決定。對于距爆心投影點相同距離的觀察點，爆高越大，沖擊波載荷越弱。

當入射角增大到臨界角附近時，沖擊波將會發生馬赫反射。馬赫反射的特點是，反射波陣面和入射波陣面的交點離開地面，形成垂直于地面的馬赫桿。馬赫反射對地面沖擊波載荷分布起到了雙重效果，一方面馬赫反射的出現使得沖擊波載荷的強度局部增大，另一方面由于馬赫桿高度的不斷增大，加快了沖擊波波陣面能量的擴散速度，使得馬赫反射區地面的沖擊波載荷衰減速度加快。因此，空中強爆炸沖擊波在地面反射后產生的大尺度空間載荷分布與爆炸當量以及爆炸高度密切相關，在實際應用中需要根據預定的目的，正確地選擇爆炸方式，可以有效地擴大沖擊波破壞區域。

3 結 論

本文中建立了基于歐拉坐標系的互不相溶、具有清晰銳利界面的可壓縮多介質流體數值方法，并研制了適用于模擬空中強爆炸沖擊波傳播的數值計算程序。利用Level Set 方法捕捉爆炸產物與空氣的運動界面，并清晰重構出分片線性流形的物質界面。采用有限體積方法求解每種物質的控制方程組，并通過在物質界面法向求解局部一維多介質Riemann 問題的精確解來計算物質界面的數值通量。數值計算中采用了網格自適應技術，在捕捉激波峰值的前提下，有效提高了計算效率。

利用計算程序首先對空中強爆炸沖擊波遇到剛性地面時的反射過程進行了數值模擬，得到了與實測數據一致的數值結果，驗證了數值模擬程序的可靠性。在此基礎上進一步對不同爆高條件下強爆炸沖擊波的反射過程進行了數值模擬，給出了不同爆高條件下地面反射沖擊波峰值超壓和沖量在不同反射區域的分布情況，得到了沖擊波載荷隨爆高的變化規律，并基于數值模擬結果給出了強爆炸沖擊波對典型目標結構毀傷破壞的最佳爆炸高度，為空中強爆炸沖擊波的殺傷破壞效果預測以及爆炸方式的選取提供了研究基礎。