空中強(qiáng)爆炸沖擊波地面反射規(guī)律數(shù)值模擬研究*

姚成寶，王宏亮，浦錫鋒，壽列楓，王智環(huán)

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

炸藥在距地面有限高度內(nèi)爆炸時(shí)，沖擊波將會(huì)在地面、建筑等障礙物的表面發(fā)生反射和繞射現(xiàn)象，傳播規(guī)律變得更加復(fù)雜。相比于化學(xué)爆炸，強(qiáng)爆炸產(chǎn)生的初始?jí)毫蛇_(dá)上萬(wàn)吉帕，初始溫度高達(dá)數(shù)十萬(wàn)攝氏度[1-2]，空氣在沖擊波作用下將發(fā)生劇烈的壓縮和溫升，熱力學(xué)性質(zhì)發(fā)生顯著變化，需要建立高溫高壓下真實(shí)氣體的狀態(tài)方程。此外，由于強(qiáng)爆炸沖擊波的傳播距離遠(yuǎn)，其在傳播過(guò)程中會(huì)受到空氣不均勻性的影響，空氣的初始?jí)簭?qiáng)和密度越低，沖擊波的衰減越明顯[2-3]?？罩袕?qiáng)爆炸問(wèn)題具有的上述特點(diǎn)，對(duì)數(shù)值計(jì)算的精度和效率提出了很高的要求。

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

在爆炸沖擊波地面反射規(guī)律研究方面，Crowl[10]給出了TNT 等化學(xué)炸藥在不同爆高下的地面反射沖擊波載荷分布規(guī)律，喬登江[1]和Glasstone 等[2]依據(jù)實(shí)驗(yàn)結(jié)果分別給出了標(biāo)準(zhǔn)大氣狀態(tài)下1 kt TNT 當(dāng)量的空中強(qiáng)爆炸沖擊波在大尺度范圍內(nèi)的等超壓曲線。目前公開發(fā)表的文獻(xiàn)大多針對(duì)TNT 等化學(xué)爆炸產(chǎn)生的沖擊波在局部范圍內(nèi)的傳播[11-15]，較少見到能夠完整給出空中強(qiáng)爆炸沖擊波在地面反射過(guò)程的大尺度范圍數(shù)值模擬。

針對(duì)空中強(qiáng)爆炸問(wèn)題的強(qiáng)對(duì)流性，本文在前期工作[16]的基礎(chǔ)上進(jìn)一步改進(jìn)基于歐拉坐標(biāo)系的多介質(zhì)流體數(shù)值方法，并考慮高溫、高壓下的真實(shí)氣體狀態(tài)方程和空氣隨高度不均勻分布的影響[17]，結(jié)合網(wǎng)格自適應(yīng)技術(shù)，計(jì)算1 kt TNT 當(dāng)量空中強(qiáng)爆炸產(chǎn)生的沖擊波在不同爆高下的地面反射完整過(guò)程，得到地面上距爆心投影點(diǎn)5 km 范圍內(nèi)的沖擊波載荷分布。通過(guò)和實(shí)測(cè)結(jié)果進(jìn)行比對(duì)，驗(yàn)證計(jì)算結(jié)果的正確性，并在此基礎(chǔ)上研究地面反射沖擊波載荷隨爆高的變化規(guī)律。

1 物理模型和數(shù)值方法

1.1 物理模型

根據(jù)空中強(qiáng)爆炸問(wèn)題的對(duì)稱性，采用二維柱對(duì)稱計(jì)算模型，控制方程組為：

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

由于強(qiáng)爆炸伴隨著強(qiáng)烈的放熱、電離等效應(yīng)，爆炸產(chǎn)物迅速膨脹至氣體狀態(tài)，且周圍空氣在受到?jīng)_擊波、熱輻射等因素的綜合作用下發(fā)生劇烈的加熱和壓縮，物理狀態(tài)極其復(fù)雜。本文中采用基于Saha 電離平衡理論給出的、以數(shù)值表形式給出的狀態(tài)方程[1]來(lái)描述強(qiáng)爆炸產(chǎn)物的熱力學(xué)行為，并采用真實(shí)氣體狀態(tài)方程來(lái)描述空氣：

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

圖 1 不同參考密度下真實(shí)氣體狀態(tài)方程的熱力學(xué)參數(shù)關(guān)系Fig. 1 Plots of equations of state for real gas at different reference densities

實(shí)際空氣的初始?jí)毫兔芏入S著高度的增加逐漸降低，且滿足下列關(guān)系[1,17]：

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

1.2 數(shù)值方法

采用基于歐拉坐標(biāo)系、具有互不相溶特征的二維柱對(duì)稱多介質(zhì)流體計(jì)算模型來(lái)進(jìn)行數(shù)值離散。在任意時(shí)刻，整個(gè)計(jì)算區(qū)域 Ω 可以分成2 個(gè)子區(qū)域(t) 和(t) ，且滿足：

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

在前期工作[16]的基礎(chǔ)上，采用Level Set 方法捕捉爆炸產(chǎn)物與空氣的物質(zhì)界面，根據(jù)Level Set 函數(shù)分片線性的特征清晰重構(gòu)出具有分片線性的物質(zhì)界面，并通過(guò)在界面的法向上求解多介質(zhì)Riemann 問(wèn)題的精確解來(lái)構(gòu)造爆炸產(chǎn)物與空氣之間的數(shù)值通量。整個(gè)計(jì)算流程如下。

1.2.1 物質(zhì)界面的追蹤和重構(gòu)

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

式中：u? 為由物質(zhì)界面的運(yùn)動(dòng)速度延拓得到的速度場(chǎng)[21]。

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

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

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

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

當(dāng)?shù)玫絫n時(shí)刻的Level Set 函數(shù) φn后，根據(jù)其分片線性的特征重構(gòu)出界面單元內(nèi)分片線性流形的物質(zhì)界面，如圖2 所示，圖中nKi,n為物質(zhì)界面的單位法向。假設(shè)單元Ki,n與單元Kj,n具有共邊Sij,n，且 包含物質(zhì)界面ΓKi,n。ΓKi,n將Ki,n和Sij,n依次分為2 個(gè)部分Ki±,n和Δtn，此時(shí)該單元包含2 種流體，且每種流體的守恒量分別定義為：

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

1.2.2 單元邊界數(shù)值通量

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

本文中采用local Lax-Friedrich 數(shù)值通量函數(shù)：

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

1.2.3 物質(zhì)界面數(shù)值通量

物質(zhì)界面數(shù)值通量是物質(zhì)界面兩側(cè)不同流體之間的數(shù)值通量，且滿足：

式中： | ΓKi,n| 為 物質(zhì)界面的長(zhǎng)度；和分別為物質(zhì)界面上的壓力和法向速度，且通過(guò)在物質(zhì)界面法向上求解局部一維多介質(zhì)Riemann 問(wèn)題的精確解來(lái)獲得，具體求解過(guò)程見文獻(xiàn)[22]。

1.2.4 守恒量的更新

當(dāng)計(jì)算得到Ki,n的單元邊界數(shù)值通量和物質(zhì)界面數(shù)值通量后，可對(duì)該單元中每種流體的守恒量進(jìn)行更新：

2 空中強(qiáng)爆炸沖擊波地面反射

2.1 空中強(qiáng)爆炸沖擊波反射數(shù)值模擬

利用數(shù)值模擬程序?qū)Ξ?dāng)量為1 kt TNT、爆高為100 m 的空中強(qiáng)爆炸沖擊波傳播過(guò)程進(jìn)行數(shù)值模擬。其中，強(qiáng)爆炸產(chǎn)物采用等溫等壓球模型，取爆炸總當(dāng)量的85%作為爆炸產(chǎn)物的力學(xué)初始能量[1]，產(chǎn)物質(zhì)量設(shè)為50 kg，初始半徑設(shè)為30 cm。空氣采用真實(shí)氣體狀態(tài)方程，取水平面處的初始?jí)毫?01.3 kPa，初始密度為1.29 kg/m3，并按照式(4)和(5)對(duì)不同高度處空氣的初始?jí)毫兔芏冗M(jìn)行設(shè)置。

計(jì)算區(qū)域在r=0 和z=0 設(shè)為固壁反射邊界條件，其余邊界設(shè)為無(wú)反射邊界條件。數(shù)值計(jì)算中采用h-網(wǎng)格自適應(yīng)技術(shù)[23]，根據(jù)相鄰網(wǎng)格的壓力梯度來(lái)進(jìn)行網(wǎng)格的局部加密或稀疏，得到典型時(shí)刻的沖擊波壓力等值線以及相應(yīng)的網(wǎng)格自適應(yīng)圖，如圖3 所示。

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

圖 3 典型時(shí)刻的沖擊波壓力等值線圖和網(wǎng)格自適應(yīng)圖Fig. 3 Pressure contours and adaptive meshes at typical times

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

2.2 不同爆炸高度的沖擊波反射規(guī)律

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

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

圖5 的計(jì)算結(jié)果表明，對(duì)1 kt TNT 當(dāng)量的空中強(qiáng)爆炸，在距爆心投影點(diǎn)300 m 范圍以內(nèi)，地面上的峰值超壓和沖量隨著爆高的增大而逐漸減小。在距爆心投影點(diǎn)300 m 范圍以外，隨著距爆心投影點(diǎn)的距離的增大，不同爆高條件下地面反射沖擊波載荷的差距逐漸變小。其中，爆高約為200 m 的空中強(qiáng)爆炸產(chǎn)生的地面峰值超壓和沖量處于最大值，其他爆高條件下的沖擊波載荷隨著爆高的進(jìn)一步增大或減小均呈現(xiàn)逐漸減弱的趨勢(shì)，其原因主要是由于空中強(qiáng)爆炸地面反射而引發(fā)的復(fù)雜沖擊波傳播引起的。當(dāng)空中強(qiáng)爆炸形成的沖擊波達(dá)到地面時(shí)，隨著入射角的增大，會(huì)依次在地面形成規(guī)則反射和馬赫反射等復(fù)雜的波系結(jié)構(gòu)。其中，在規(guī)則反射區(qū)，反射沖擊波載荷由入射波強(qiáng)度和入射角決定。對(duì)于距爆心投影點(diǎn)相同距離的觀察點(diǎn)，爆高越大，沖擊波載荷越弱。

當(dāng)入射角增大到臨界角附近時(shí)，沖擊波將會(huì)發(fā)生馬赫反射。馬赫反射的特點(diǎn)是，反射波陣面和入射波陣面的交點(diǎn)離開地面，形成垂直于地面的馬赫桿。馬赫反射對(duì)地面沖擊波載荷分布起到了雙重效果，一方面馬赫反射的出現(xiàn)使得沖擊波載荷的強(qiáng)度局部增大，另一方面由于馬赫桿高度的不斷增大，加快了沖擊波波陣面能量的擴(kuò)散速度，使得馬赫反射區(qū)地面的沖擊波載荷衰減速度加快。因此，空中強(qiáng)爆炸沖擊波在地面反射后產(chǎn)生的大尺度空間載荷分布與爆炸當(dāng)量以及爆炸高度密切相關(guān)，在實(shí)際應(yīng)用中需要根據(jù)預(yù)定的目的，正確地選擇爆炸方式，可以有效地?cái)U(kuò)大沖擊波破壞區(qū)域。

3 結(jié) 論

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

利用計(jì)算程序首先對(duì)空中強(qiáng)爆炸沖擊波遇到剛性地面時(shí)的反射過(guò)程進(jìn)行了數(shù)值模擬，得到了與實(shí)測(cè)數(shù)據(jù)一致的數(shù)值結(jié)果，驗(yàn)證了數(shù)值模擬程序的可靠性。在此基礎(chǔ)上進(jìn)一步對(duì)不同爆高條件下強(qiáng)爆炸沖擊波的反射過(guò)程進(jìn)行了數(shù)值模擬，給出了不同爆高條件下地面反射沖擊波峰值超壓和沖量在不同反射區(qū)域的分布情況，得到了沖擊波載荷隨爆高的變化規(guī)律，并基于數(shù)值模擬結(jié)果給出了強(qiáng)爆炸沖擊波對(duì)典型目標(biāo)結(jié)構(gòu)毀傷破壞的最佳爆炸高度，為空中強(qiáng)爆炸沖擊波的殺傷破壞效果預(yù)測(cè)以及爆炸方式的選取提供了研究基礎(chǔ)。