爆炸沖擊波反射流場的理論計算方法*

賈雷明，王澍霏，田 宙

（1. 清華大學工程物理系，北京 100084；

2. 西北核技術研究，陜西 西安 710024）

爆炸沖擊波是爆炸對周圍目標產生破壞效應的重要因素之一，自20 世紀40 年代以來就受到廣泛關注[1]。爆炸沖擊波在空氣中傳播，遇到地面或其他結構時會發生反射。為了準確評估破壞效應，需要確定沖擊波反射后的流場分布。本文中擬建立一種理論方法，用于對任意當量、爆高條件下沖擊波在固壁反射后的流場物理量進行快速計算。

常用的快速計算方法，大致分為2 類。一類是基于數值或實驗數據建立物理量的擬合關系式[2-5]，該方法的精度取決于數據規模、自變量數目和關系式類型等，所得到的關系式通常具有一定的適用范圍和場景。另一類方法是將固壁假想為對稱面，在固壁另一側建立與真實爆炸對稱的虛擬爆炸，然后根據一維球爆炸自由場參數確定真實和虛擬爆炸流場，再將兩者進行線性或非線性疊加，從而得到沖擊波反射后的流場分布，這種方法稱為鏡像方法[6-12]。與第一類方法相比，鏡像方法更易于建立和使用，可以給出較清晰的物理圖像，且不受限于特定的爆炸場景。

爆炸沖擊波在固壁可能發生正規或非正規反射，對應的流場分布及演化是不相同的。在鏡像方法中，目前已有的工作大都假定沖擊波始終保持球面，即僅考慮了正規反射的情形，對于非正規反射的情形則沒有涉及。另外，壓力和動壓均是評估沖擊波毀傷的重要參量，其中流場壓力可采用非線性疊加模型LAMB[12]計算，而對于密度和速度則采用線性疊加計算，所得到的結果與真實值相差較大[2]。

為了能夠將鏡像方法推廣到任意反射類型時爆炸流場的計算，本文中建立正規和非正規反射階段流場中沖擊波傳播的計算方法，并基于LAMB 模型發展新的疊加模型LAMBR (LAMB revisied)。針對任意爆炸當量、爆高條件，首先確定任一時刻流場中的沖擊波結構，然后利用鏡像方法和LAMBR 疊加模型，快速計算流場中壓力、密度和速度等物理量分布。通過與數值模擬結果和美國結構抗偶然爆炸手冊(UFC 3-340-02，后文簡稱為UFC)[13]中的數據比對，驗證該理論方法的合理性。

1 理論方法

1.1 沖擊波傳播

當爆炸沖擊波傳播到達固壁時，首先發生正規反射。在t=t0時刻，入射波陣面與固壁夾角超過某一臨界值，正規反射轉換為非正規反射。常用的轉換準則，包括脫體準則、聲速準則和力學平衡準則，考慮到爆炸沖擊波在固壁反射屬于非定常流動，本文中選取聲速準則[14]計算非正規反射的起始位置T0。在t＞t0時刻，三波點T 逐漸向著遠離固壁的方向運動，見圖1。圖1 中O 為爆炸中心，H 為爆高，固壁與y=0 重合，RR、IR 分別表示正規反射和非正規反射。

對于正規反射和非正規反射，均假定反射波是以點(0,-H)為中心的圓弧，根據不同時刻入射波傳播距離，確定反射波形狀。對于非正規反射，假定馬赫桿是以爆心在固壁投影點為中心的圓弧，通過求解三波點T 的位置，確定馬赫桿形狀。已有的計算點T 位置的方法，可以分為基于幾何相似[15]和基于流體動力學[16-19]兩類。

1.1.1 基于幾何相似的方法

圖1 中，三波點由點T0開始，隨時間發展逐漸遠離壁面，在t1時刻到達點T。此時，入射波陣面延長線、馬赫桿分別與固壁交于點B、F。點A 為直線OB 與t0時刻沖擊波陣面的交點，直線AC 與OT0平行，CE 與OB 垂直，CD、EF 與固壁垂直，則△BFE 與△BDC 是相似的。馬赫桿與固壁交點F 的坐標，可以表示為：式中：R=R(t)為 t 時刻入射波半徑，θ0、θB分別為入射波陣面在點T0、B 處的傳播方向與水平方向的夾角(以逆時針方向為正)，g 為與爆炸條件相關的常數。三波點T 的坐標為：

圖 1 爆炸沖擊波反射示意圖Fig. 1 Schematic diagram of blast wave reflection

易仰賢[15]假定常數g 僅與比例爆高H= H=W1=3有關(W 為炸藥質量)，并給出了核爆炸情況下的關系式g= g( )，但是該關系式并不能直接應用于化爆情況。因此，本文中根據UFC[13]中提供的三波點坐標數據，按照上述過程反推得到適用于化爆情況的g。定義絕對偏差：

需要說明的是，由于UFC 中沒有提供非正規反射起始點T0的坐標信息，因此本文中基于聲速準則和UFC 中爆炸自由場數據計算點T0的位置。

1.1.2 基于流體動力學的方法

在任意時刻 t (t＞t0)，如果已知三波點T 的坐標 ( xT, yT)和該時刻三波點跡線與水平方向的夾角χ，則有：

表 1 不同比例爆高條件下g 的取值Table 1 The value of g at different scaled heights of explosion

因此，給定三波點T 起始位置T0，按時間逐步推進，即求得三波點軌跡。

將三波點附近流動近似看作擬定常流動，可以采用von Neumann 的三波理論[16]求解χ。點T 的運動速度DT與χ 之間滿足：

式中：DI為入射波傳播速度。建立與點T 固連的坐標系 ，并在坐標系 中觀察流場運動，則穿過入射波后流動方向偏轉角度δI滿足：

式中：ξ 為波后與波前流場壓力之比，Ma0為波前流場馬赫數。穿過反射波和馬赫桿后流動方向偏轉角度，也按式(6)計算。根據界面兩側連續條件，聯立沖擊波關系式，建立迭代求解過程，即可得到χ。

也可基于幾何激波動力學(geometrical shock dynamics, GSD)[20]理論求解χ[17-19]。在三波點兩側，入射波與馬赫桿的傳播角度θI和θM滿足：

式中：AK= f (MaK)，下標K=I，M 分別對應入射波和馬赫桿。通過求解上式得到MaM，進而有：

在式(7)和(8)中，關系式AK=f (MaK) 可參照Itoh 等[21]文中式(26)確定。

為了記述方便，將基于幾何近似和表1 中gUFC值的方法記為Yi-UFC 方法，將基于von Neumann 三波理論的流體動力學方法記為VN 方法，將基于GSD 理論的流體動力學方法記為GSD 方法。通過與數值結果和UFC 數據比對，從中選取合理的方法。

1.2 疊加模型LAMBR

入射波后、反射波前的流場，與真實爆炸自由場無異。反射波后的流場，則利用鏡像方法計算，見圖2，將真實爆炸和虛擬爆炸自由場疊加，即可給出反射波后的流場分布。

當沖擊波在固壁發生正規反射時，見圖2(a)，對于反射波后某一點P，記真實爆炸自由場在該點處的物理量為p1、ρ1、u1，虛擬爆炸自由場在該點處的物理量為p2、ρ2、u2，流體滿足完全氣體狀態方程p = (γ-1)ρe，γ 為比熱比。則點P 處物理量按下述模型計算。

圖 2 鏡像法示意圖Fig. 2 Diagram of the method of image for blast wave reflection

1.2.1 壓力

按LAMB 模型計算流場壓力，有：

式中：p0、ρ0分別為入射波前未擾動流場的壓力、密度。

對于沖擊波在固壁發生正反射的情況，按照式（9）計算得到的壓力與精確解完全一致。當波后流場密度較小時，可能得到ρ＜0，進而給出數值極大、不符合物理的p。本文中建議選取。

1.2.2 密度

跨越沖擊波是一個非線性且熵增的過程，因此波后密度ρ 大于線性疊加得到的密度，同時又小于流體經歷等熵過程后的密度。因此，根據p 與p1的相對大小計算波后密度ρ：

1.2.3 速度

根據動量守恒計算速度預估值：

考慮到式(12)本質上是線性疊加，因此僅適用于壓力較小時的情況，需要引入一定的修正。本文中利用激波關系式得到固壁處沖擊波陣面后的速度，然后在保證該速度幅值不變的前提下，對式(12)給出的流場速度的幅值進行等比例修正，給出流場速度u。

將該模型記為LAMBR (LAMB revised),其與LAMB 模型的主要區別在于密度和速度的計算。當沖擊波在固壁發生非正規反射時，見圖2(b)，仍然采用LAMBR 模型計算點P 處的物理量。與正規反射不同，此時反射波R 超越入射波I 傳播，即反射波半徑大于入射波半徑。在確定點P 處的爆炸自由場物理量時，應先根據當前時刻反射波半徑與入射波半徑之比，對點P 到爆心距離進行等比例縮放，再根據修正后的爆心距確定自由場物理量。

2 結果與討論

利用本文中建立的理論方法，計算1 kg TNT 爆炸沖擊波在固壁反射后的流場，并與UFC 數據和動力學軟件AUTODYN 計算結果進行對比，檢驗理論方法的合理性。初始時刻空氣壓力p0=101.325 kPa，密度ρ0=1.225 kg/m3，比熱比γ=1.4，且處于靜止狀態。炸藥TNT 選用JWL 狀態方程描述：

式中： V=ρ0/ρ 為相對體積，e 為質量內能，A、B、R1、R2和w 為常數。TNT 的JWL 狀態方程參數[22]：初始密度ρ0=1 630 kg/m3, A=371.2 GPa, B=3.231 GPa, R1=4.15, R2=0.95, w=0.30, 爆速D=6.930 km/s, 初始體積能量E0=7 GPa, CJ 壓力pCJ=21 GPa。當炸藥體積膨脹到約10 倍初始體積時，為了減小計算誤差，AUTODYN 會改用理想氣體狀態方程來描述炸藥產物的力學行為，取比熱比γTNT=1.3。

2.1 爆炸自由場

圖 3 1 kg TNT 爆炸自由場沖擊波參數Fig. 3 Free field blast wave parameters for 1 kg TNT

在AUTODYN 中建立一維軸對稱楔形模型（rmin=1.00×10-3m, rmax=6.00 m），利用Euler-Multi-Material求解器，計算1 kg TNT 在空氣中的爆炸自由場，炸藥半徑r=5.27×10-2m。采用2 種網格劃分方案，網格尺寸dr 分別取為1.00×10-4m 和5.00×10-5m。圖3中給出了爆炸自由場沖擊波峰值超壓Δps、走時ta隨爆心距 r 變化的計算結果，并列出了UFC 中的數據。利用不同網格方案得到的結果基本重合，驗證了數值計算的網格收斂性。在爆心距 0.3 m ≤ r≤ 3.0 m 范圍內，計算結果與 UFC 中數據基本吻合，驗證了數值計算的合理性。在后文的理論分析和數值模擬中，均選取dr = 5.00×10-5m 的計算結果作為爆炸自由場數據。

2.2 爆炸沖擊波在固壁反射

選取4 種不同的爆高條件，見表1，利用理論方法對爆炸流場進行計算。同時，在AUTODYN中建立二維軸對稱模型，利用Euler-Multi-Material 求解器進行數值計算。計算區域為0≤x≤4H, 0≤y≤3H，以y 軸為對稱軸，在邊界x=4H, y=3H 處施加 outflow 邊界條件，在邊界y=0 處施加固壁邊界條件，爆心位于(0, H)處。采用2 種網格劃分方案，網格尺寸分別取dx=dy=H/500 和 dx=dy=H/667。初始時刻流場，由t*時刻爆炸自由場映射得到，要求沖擊波傳播距離R ＜H。4 種爆高條件對應的t*，分別為 9.60×10-5、1.88×10-4、3.16×10-4和 4.80×10-4s 。

2.2.1 三波點軌跡

圖4 給出了不同爆高條件下三波點軌跡的理論和數值結果，并列出了UFC 中數據。其中，dx=dy=H/500, dx=dy=H/667 分別為相應網格劃分方案下的數值結果，通過讀取不同時刻壓力云圖中三波點位置得到。

圖 4 不同爆高情況下三波點軌跡Fig. 4 Triple point trajectories at different heights of explosion

對比不同網格方案得到的結果，兩者基本重合，驗證了數值計算的網格收斂性。對比數值結果和UFC 數據，發現在不同爆高H 條件下，數值結果與UFC 數據在yT＜0.4～0.5 m 范圍內存在一定程度的偏差，可能是由于數值模型對爆炸產物或物理過程的描述與真實情況之間還存在差異，也可能是由于UFC 數據本身的誤差。隨著三波點遠離固壁或是爆高增加，兩者趨于吻合，表明數值結果是合理的。

將理論結果與UFC 數據、數值結果進行對比。首先，對比非正規反射起始點位置，見表2，理論與數值計算結果基本一致，表明選取聲速準則作為正規與非正規反射的轉換準則是合理的。而UFC 中缺少該數據，因此未參與比較。接下來，對比三波點軌跡。VN 與GSD 同屬流體動力學方法，僅三波點附近流場求解方法不同，因此當H=0.396 7 m 和0.595 1 m 時，兩者的結果基本一致，但與UFC 和數值結果差異較大；當H=0.793 4 m 和0.991 7 m 時，兩者的結果在靠近非正規反射起點附近基本一致，隨著三波點遠離壁面，VN 結果略大于GSD。與VN 和GSD 相比較，對于不同爆高條件，Yi-UFC 方法的結果與UFC 數據和數值結果的吻合程度均優于VN 和GSD 方法，表明Yi-UFC 方法能夠更好地描述三波點的運動。因此，選取Yi-UFC 方法用于流場沖擊波結構的理論計算。

表 2 非正規反射起始點坐標xT0Table 2 The xT0 value of the starting point of IR

2.2.2 不同時刻流場分布

以爆高H=0.793 4 m 為例，圖5、6 分別給出了不同時刻流場壓力云圖以及沿直線y=0.03 m 的物理量分布。其中，數值計算選用dx=dy=H/667 的網格劃分方案。

在t=4.76×10-4s 時刻，入射波在固壁發生正規反射。理論方法給出的流場壓力云圖與數值結果較接近，見圖5(a)，能夠反映流場中的波系結構。在真實情況中入射波后流場是非均勻的，反射波的形狀勢必不能保持圓弧；同時，反射波會在物質界面處發生折射，形成新的反射波和透射波，物質界面的運動速度也逐漸減小。理論方法近似假定反射波為圓弧，且物質界面按照原一維球對稱傳播規律運動，這與真實情況有差異，因此理論結果中反射波和物質界面的形狀與數值結果不完全一致。直線y=0.03 m 依次與反射波、入射波相交，沿該直線的物理量分布出現2 個峰值，見圖6(a)，基于同樣的原因，理論結果中反射波位置與數值結果相差約為0.03 m。除此之外，壓力、速度的理論與數值結果基本重合。

在t=9.96×10-4s 時刻，入射波在固壁發生非正規反射。對比流場壓力云圖，理論與數值結果的差異主要集中在對稱軸附近區域，而在其他區域則較接近，見圖5(b)。對比沿直線 y=0.03 m 的物理量分布，理論與數值結果中的馬赫桿位置基本重合，表明理論方法中將馬赫桿近似為圓弧是合理的，與真實情況較符合，這一點在后續時刻的比較中也可以得到說明。Hu 等[16]曾指出，當爆高小于某一臨界值時，入射波在固壁的非正規反射會依次經歷雙馬赫反射、過渡馬赫反射和單馬赫反射等3 個階段。根據數值結果判斷，此時流場中的反射為過渡馬赫反射，沿直線 y=0.03 m 的物理量分布出現了2 個峰值，見圖6(b)。基于反射波的圓弧假定，理論方法僅能構建單馬赫反射類型，所以在過渡馬赫反射的影響區域內，理論結果中并未出現第2 個峰值。

圖 5 H=0.793 4 m 時，不同時刻流場壓力云圖(黑線為物質界面)Fig. 5 Pressure contours at various time instants (black line stands for material interface) for H=0.793 4 m

圖 6 H=0.793 4 m 時，不同時刻流場物理量沿y=0.03 m 的分布Fig. 6 Parameters along y=0.03 m at various time instants for H=0.793 4 m

在t=1.92×10-3s 時刻，非正規反射類型轉換為單馬赫反射，此時流場波系結構與理論方法中的假定是一致的，因此流場壓力分布的理論與數值結果較吻合，見圖5(c)，沿直線y=0.03 m 物理量分布的理論和數值結果也趨于吻合，見圖6(c)。此時，在爆心附近出現次級沖擊波并向外傳播，形成高壓區域。受到反射波的影響，數值結果中該高壓區域的位置相較于理論結果更靠近爆心上方。考慮到流場發展本身具有極強的非線性，而本文中所用的疊加方法并不能完全精確地求解流場中的非線性作用，這就導致在非線性作用顯著的區域，諸如靠近爆心和對稱軸附近的區域，理論與數值結果存在一定的偏差。

在t=3.60×10-3s 時刻，流場中的次級沖擊波向外傳播，并在固壁發生反射，同時次級沖擊波后流場壓力逐漸減小。理論方法給出的流場壓力云圖能夠反映這一特征，但是爆心附近的壓力小于數值結果。對比沿直線y=0.03 m 物理量分布，見圖6(d)，在次級沖擊波反射點之外的區域，即x＞1.47 m 時，理論結果與數值結果基本吻合；而在x＜1.47 m 的區域中，理論結果則小于數值結果，次級沖擊波的位置相差約0.46 m。

通過以上對比可以發現，盡管理論方法對于流場發展過程進行了一定簡化，但無論流場發展是處于正規反射階段還是非正規反射階段，理論方法所得到的結果均合理地反映了流場中波系發展特征，尤其是在靠近沖擊波的區域，不同時刻的物理量分布規律與數值結果基本吻合，物理量峰值與數值結果的偏差大都小于10%，見表3，從而驗證了理論方法的合理性。表3 中同時列出了基于LAMBR 和LAMB 模型得到的理論結果，兩者所給出的壓力峰值相同，但對于密度和速度峰值，LAMBR 模型的結果明顯優于LAMB 模型。

表 3 沿直線y=0.03 m 的物理量峰值Table 3 Peak values of parameters along the line y=0.03 m

在運行時間上，以截止時刻t=3.60×10-3s 為例，采用相同的計算區域和網格劃分方案，數值計算耗時約70 h，理論方法僅需約15 s，遠遠快于數值計算。

3 結 論

爆炸沖擊波遇到固壁發生反射，使得波后流場改變。本文中基于鏡像法，在固壁另一側建立與真實爆炸對稱的虛擬爆炸，將爆炸沖擊波在固壁的反射等效為真實與虛擬爆炸流場的相互作用，建立了適用于爆炸沖擊波反射后流場的快速計算方法。

該方法首先將反射波、馬赫桿分別簡化為以虛擬爆源和地面爆心投影點為中心的圓弧，選取聲速準則作為正規與非正規反射類型的轉換準則，利用基于幾何近似的Yi-UFC 方法，計算不同時刻流場中沖擊波結構。然后，根據沖擊波結構，將流場劃分為不同的區域。對于位于反射波前、入射波后的區域，直接根據爆炸自由場數據確定流場分布。對于反射波后的流場，則通過本文新發展的流場疊加模型LAMBR 進行計算。通過與數值結果和已有數據進行對比，理論方法給出的壓力、密度和速度等物理量分布規律、峰值等與數值結果基本吻合，能夠識別流場中的次級沖擊波傳播等物理過程，驗證了理論方法的合理性。與LAMB 模型相比，LAMBR 模型得到的密度和速度等物理量的峰值與數值結果更為吻合。同時，該理論方法所需要的時間，較數值計算大幅度減少。

需要注意的是，流場發展本身是非線性極強的物理過程，鑒于理論方法所采用的諸多假定，在流場非線性作用顯著的區域中，理論結果與真實情況還存在一定差異。