圓柱裝藥近地空爆沖擊波峰值超壓空間轉換模型

郗洪柱，孔德仁，樂貴高，史 青，彭泳卿

（1.南京理工大學機械工程學院，江蘇 南京 210094；2.北京遙測技術研究所，北京 100076）

高能戰(zhàn)斗部爆炸毀傷威力大，產生的高能量沖擊波超壓是其殺傷目標的主要毀傷元之一，受到相關研究人員的重視。學者們已提出一系列爆炸沖擊波峰值超壓的經驗公式，如葉曉華公式[1]、Mills 公式[2]、Henrych 公式[3]、Sadovskyi公式[4]、Brode公式[5]及國防工程設計規(guī)范[6]等。上述公式描述了自由場峰值超壓的一維分布規(guī)律。在近地空爆（本研究中均指比例爆高不小于0.35且不大于3的空中自由場爆炸）中，存在入射波、反射波、馬赫波及三波疊加現(xiàn)象，使波陣面峰值超壓的空間分布不均勻，無法直接應用上述公式。

根據峰值超壓數(shù)值大小關系，最大者為三波疊加即三波點超壓，其次為馬赫波超壓。三波點軌跡隨沖擊波傳播的時間演化而逐漸升高。通常有3種確定三波點軌跡、超壓及馬赫超壓的方法。第1種方法是試驗法，又可分為兩類：一類是通過設置足夠多的超壓傳感器測點，根據測點上的超壓變化關系，利用多次試驗得到三波點近似位置，同時得到馬赫超壓[7]；另一類是利用高速紋影技術得到爆炸波傳播過程圖像，從中解讀出三波點位置，并推算三波點跡線[8]。該方法的試驗成本較大，不易推廣。第2種方法是基于Whitham 的幾何激波動力學理論[9]，利用邊界條件，計算得到近似解[10–11]。第2種方法及其改進方法已有計算機程序，可在較短的時間內獲取三波點跡線及馬赫波超壓[10]。第3種方法是鏡像法[12–14]，即將近地爆炸視為真實空間中爆炸與相對于地面對稱的虛擬爆炸的疊加，可求得沖擊波反射流場參數(shù)。與前兩種方法相比，第3種方法不僅試驗成本較低，而且方便建模分析，尤其適合中小當量裝藥和爆高較大的工況[15]。

應用上述方法得到馬赫反射超壓，結合自由場超壓經驗公式，可在一定程度上描述爆炸沖擊波空間分布規(guī)律。然而在近地空爆中，目前少有直接研究馬赫反射超壓與入射波超壓空間分布關系的工作。揭示超壓空間分布關系，對沖擊波超壓毀傷元空間威力評估具有重要意義。

為了建立與空間位置相關的近地自由場沖擊波超壓轉換模型，本研究采用鏡像法確定與爆高平齊時的三波點位置，將其作為建模起始邊界，將馬赫反射終點作為建模結束邊界。將測點與爆心連接起來，該連線與過爆心的水平面的銳角夾角定義為測點角度，簡稱測角。在起始和結束邊界之間的區(qū)域，應用測角等分和超壓歸一化值分段線性假設得到超壓空間歸一化關系，從而建立超壓空間轉換模型，并將模型推廣至不同爆高、當量、長徑比的圓柱裝藥近地空爆中。利用符合經驗公式和實爆結果的AUTODYN-2D數(shù)值模型獲取試驗數(shù)據，確定模型形式并求解系數(shù)；通過將模型轉換結果與數(shù)值結果和實爆數(shù)據對比，驗證該模型的可靠性。

1 理論分析

近地空中爆炸產生的沖擊波在到達地面前符合自由場傳播規(guī)律，到達地面后依次經歷正反射、正規(guī)反射，入射角超過一定角度后發(fā)生馬赫反射。以長徑比為1的圓柱裝藥TNT 為研究對象，采用鏡像法分析近地自由場爆炸時空中各點超壓關系。如圖1所示，爆高為h，爆心為O，爆心在剛性地面投影點為P，虛擬爆源OM與O相對于地面鏡像對稱。剛好發(fā)生馬赫反射時，入射波Ⅰ與反射波Ⅱ在地面的交點為S。Th為入射波Ⅰ、反射波Ⅱ和馬赫波Ⅲ的交點，位于過爆心的水平線OO′上；TM為虛擬三波點，位于虛擬水平線OMO′M上。弧線STd為真實三波點軌跡，STM為虛擬三波點軌跡。直線ThOM與OO′的銳角夾角記為αⅡ，該角與反射波Ⅱ的反射角互為余角。

圖1 近地爆炸沖擊波近距離傳播流場示意圖Fig.1 Near-earth blast shockwave propagation flow field at close range

水平線OO′上有測點2，與爆心O的直線距離為R2；測點1與測點2在同一波陣面上，與爆心O的直線距離為R1，兩測點高度差為H1，夾角為θ；測點3與測點1在過爆心O的同一直線上，與爆心O的直線距離為R3；過水平線OO′上有測點4，與測點3在同一波陣面上，高度差為H2，與爆心O的直線距離為R4。顯然，測點4在三波點以下，其超壓大于同一波陣面上的測點3；而測點1、測點2和測點3均在三波點軌跡以上，前兩者超壓接近。測點2處峰值超壓可用Sadovskyi 公式[4]計算

式中：rˉ 為比例距離，是爆心距R與炸藥當量m的三次根之比，m·kg?1/3。用RⅠ表 示OTh長度，RⅡ表示OMTh長度，則根據圖1的限定條件，點Th處反射角余角 αⅡ滿足

式中：DI和DII分別為入射和反射波陣面?zhèn)鞑ニ俣?，MI和MII為對應的馬赫數(shù)。馬赫數(shù)M與超壓Δp滿足

故可建立 αⅡ與 入射波陣面超壓（?pⅠ） 及反射波陣面超壓（?pⅡ）的關系

因此式(4)可表示為

將一系列比例距離代入經驗公式（如Sadovskyi 公式）中可得到一系列超壓值，利用 αⅡ的約束方程

通過迭代法可求得Th處爆心距RⅠ。 按照圖1規(guī)定坐標軸方向，以P為坐標原點時Th點坐標為(RⅠ,h)。

隨著時間推移，Th點之后三波點的縱坐標大于Th處的縱坐標。當馬赫反射超壓與兩倍當量自由場爆炸沖擊波峰值超壓相等時達到馬赫反射終點[16]，將馬赫反射終點記為Td，如圖2所示。

圖2 近地爆炸沖擊波遠距離傳播流場示意圖Fig.2 Near-earth blast shockwave propagation flow field over long distance

理想空氣中馬赫波地面超壓公式[16]為

式中：A1、A2和A3為無量綱系數(shù)，基于試驗數(shù)據確定（或直接應用經驗公式）。

聯(lián)立式（8）、式（9）、式（10）可求得Td處R值。為簡便起見，用Td在水平線OO′上的垂足Td′表示馬赫反射終點界限，并將Td′坐標約定為(R,h)。

針對區(qū)域二超壓分布的復雜性，以O為中心，對90°平面進行每15°等分。建立7種測角下峰值超壓與比例距離的關系式，然后以90°超壓為參考，其他方向峰值超壓與參考值求比值，測角θ相同時的比值稱為峰值超壓歸一化值，記為yθ。建立該值與測角 θ的線性關系式

根據歸一化值的定義，在已知某測角處超壓值和式（11）的基礎上，可得到任意測角α處相同比例距離的超壓，即

對式（12）進行擴展，不考慮地面材質、不平度、環(huán)境等因素，區(qū)域1和區(qū)域2中不同空間位置的峰值超壓與圓柱裝藥長徑比k、爆高h、當量m及由測角θ和比例距離rˉ表示的空間位置有關。長徑比增加引起入射波形狀改變，即等壓線曲率變化，應對空間不同測角處峰值超壓進行修正。以kˉ表示與長徑比有關的修正因子，以比例爆高hˉ表 示爆高和當量的影響，以?pθ表示已知測角θ處峰值超壓，將待求測角α處峰值超壓?pα表示為上述影響因素的函數(shù)，即

2 數(shù)據來源

AUTODYN 在爆炸問題的求解上有廣泛應用[17–20]，利用該軟件建立典型環(huán)境圓柱裝藥空中爆炸數(shù)值模型并獲取數(shù)據，利用經驗公式和實爆試驗數(shù)據驗證數(shù)值模型，利用數(shù)值結果完成式（12）的構建，利用數(shù)值結果和實爆數(shù)據評估所構建的超壓轉換模型性能。

2.1 AUTODYN數(shù)值模型

利用沖擊波傳播的對稱性特點，在AUTODYN-2D中構建二維1/2平面模型。研究對象為不同當量、長徑比、中心起爆的圓柱裝藥TNT 近地自由場爆炸（比例爆高不小于0.35）。x軸沿TNT 模型圓柱軸線方向，如圖3(a)所示；地面模型沿y軸放置，并位于空氣模型左側，如圖3(b)所示。

圖3 數(shù)值模型Fig.3 Numerical model

空氣模型網格粗細與爆炸沖擊波求解的精確性正相關，與求解速度負相關。當網格尺寸減小為10和5 mm 時，兩者的超壓時程曲線基本一致[21]。故網格不大于10 mm 時可較準確計算超壓。為確保沖擊波快速精確計算，建立兩種尺度和網格粗細的空氣模型。第1個空氣模型小，僅包裹TNT，網格尺寸為1 mm；第2個空氣模型較大，包裹第1個空氣模型，網格尺寸為10 mm?？諝夂偷孛婺Ｐ途捎枚S多物質歐拉方法求解。利用歐拉和拉格朗日方法實現(xiàn)各模型間的網格自動連接。對稱面不設置邊界條件，自由邊界采用壓力流出（Flow-out）類型。計算時間根據TNT當量及所需傳播的距離確定，設置為20～150 ms不等。

TNT采用材料庫中的“TNT”，密度為1.630 g/cm3；空氣模型采用材料庫中的“AIR”，內能設置為2.068×105J，密度為1.225 kg/m3；地面模型采用材料庫中的“SAND”；其余參數(shù)均按照AUTODYN默認設置。采用JWL 狀態(tài)方程描述TNT 的爆轟過程。JWL 狀態(tài)方程為

式中：p為爆轟壓力，V為相對體積，E為單位體積內能，A、B、R1、R2和 ω為材料常數(shù)。JWL 方程參數(shù)列于表1[22]。

表1 TNT爆炸的JWL狀態(tài)方程參數(shù)[22]Table 1 Parameters of JWL equation of TNT detonation[22]

2.2 數(shù)值方法驗證

以爆心為起點，測點呈輻射狀分布，將90°范圍的二維仿真平面進行分隔。從0°開始，每15°設置一系列角度相同但比例距離逐漸增大的測點。

經驗公式與TNT實爆試驗的入射波峰值超壓相符[2]，可用經驗公式評估數(shù)值模型的正確性。將1 kg TNT 爆高4 m 的數(shù)值仿真結果同經驗公式求解結果對比，如圖4所示，圖中橫坐標為比例距離，縱坐標為峰值超壓；數(shù)值結果用五角星表示，經驗公式計算結果用其他形狀表示。其中，數(shù)值仿真結果與Henrych 公式的一致性好，平均相對偏差僅為8.9%；與葉曉華公式、Sadovskyi 公式及國防工程設計規(guī)范的一致性均較好，其中相對于Sadovskyi 公式的偏差最小，平均相對偏差為9.8%；相對于國防工程設計規(guī)范和葉曉華公式的平均相對偏差分別為15.0%和16.8%。

圖4 數(shù)值模型的驗證Fig.4 Validation of numerical model

2020年10月中下旬，團隊在中部某靶場進行了10 kg 長徑比為1.0的圓柱裝藥TNT空中實爆（爆高1.5 m）試驗，地面為平整沙土地，現(xiàn)場微風，天氣晴朗，如圖5所示。中間懸吊TNT，下方為平整鋼板，自由場傳感器（PCB137A）與爆高等高并指向爆心。實爆、數(shù)值仿真結果及仿真誤差如表2所示。表中實測超壓值為有效重復試驗測試結果的均值，其中工況1和工況2為3次重復試驗測試結果的均值，工況3和工況4為兩次重復試驗結果均值。傳感器標定后線性度不大于滿量程的1%；以線纜將超壓信號傳輸至存儲測試儀中，采樣率為1 MHz。相同比例距離時仿真結果偏??；仿真較實爆值的最大偏差約為?15%；其余測點偏差在?10%以內。

表2 實爆及數(shù)值模擬結果對比Table 2 Real blast and numerical simulation results

綜上，數(shù)值模型與經驗公式和實爆結果均較吻合，表明構建的數(shù)值模型及相關設置合理，可代替實爆試驗作進一步研究。

圖 5 試驗現(xiàn)場Fig.5 Testing site

3 模型參數(shù)求解

3.1 峰值超壓與比例爆高及測點位置的關系

利用100 kg長徑比為1的圓柱裝藥TNT在爆高3 m 下的試驗數(shù)據，建立7種測角的峰值超壓與比例距離的關系式。7 個方程系數(shù)項如表3所示。

表3 系數(shù)匯總表Table 3 Summary of coefficients

以90°測點結果為參考，其他方向峰值超壓與參考值求比值。依此方法分別對比例爆高為0.538、0.646、2.480和2.977的4種工況進行計算，所得結果見圖6。橫坐標表示測角，縱坐標為相同比例距離處峰值超壓與90°測點超壓之比，即峰值超壓歸一化值。圖6中各曲線對應的比例距離為1.5～8.0。

圖6 不同比例爆高時峰值超壓的歸一化值Fig.6 Normalized value of peak overpressure under different scaled height of burst (HOB)

假設在比例爆高不變的條件下，歸一化值與測角呈線性關系；在測角不變時，歸一化值與比例爆高呈線性關系；此外，由于45°前后歸一化值與比例爆高的關系明顯不同，以45°為界分為兩段。

固定比例爆高，計算45°內的平均歸一化值，將該值對應的測角記為22.5°，然后計算45°處平均歸一化值，由此確定45°內歸一化值與測角的關系。通過不同比例爆高與22.5°和45°處平均歸一化值的線性擬合關系，將比例爆高引入。最終建立的歸一化值方程為

基于線性最小二乘法將圖6所示典型工況的比例爆高、測角和對應歸一化值數(shù)據代入式(15)，得到對應系數(shù)。將45°～90°內平均歸一化值對應測角設定為67.5°，同理得到對應方程和系數(shù)。完整的歸一化值方程為

3.2 峰值超壓與長徑比及測點位置的關系

建立1 kg 當量圓柱裝藥TNT 爆高0.5 m，長徑比分別為0.5、1.0、1.5和2.0時的數(shù)值模型。以長徑比1.0時各測點的峰值超壓為基準，計算其他長徑比時測點峰值超壓相對基準的偏差，以該偏差為縱坐標，以比例距離為橫坐標，得到圖7(a)所示結果。近場最大正偏差約15%，最小負偏差約?7%；中遠場最大正偏差約7%，最小負偏差約?3%；遠場最大正偏差約10%，最小負偏差在?10%以內。

顯然，近場峰值超壓受長徑比影響最明顯；中遠場受影響稍??；遠場由于其峰值超壓接近零，即使差值不大，經比值運算后也會表現(xiàn)出一定誤差。本研究主要修正45°內比例距離小于4時的峰值超壓轉換誤差。比例距離小于4時，長徑比為1.5、2.0、3.0和5.0與長徑比為1.0時對應測角和比例距離處測點的峰值超壓求比值，結果記為Qk′，下標k′表示參與求比值的兩個長徑比，如長徑比為5.0與長徑比為1.0時對應0°處測點峰值超壓之比記為Q5.0/1.0。為便于比較，取長徑比0.5與1.0的峰值超壓比值之倒數(shù)為Qk′，Qk′與比例距離的關系如圖7(b)所示。圖中展現(xiàn)出的變化趨勢較一致，長徑比相差小的Qk′較長徑比相差大者隨比例距離的變化更小。進一步，用線性模型構建Qk′與比例距離的關系。

圖7 相對于長徑比（k）1.0的變量Fig.7 Variables that relative to the length diameter ratio(k)of 1.0

基于多項式擬合法建立0°處Q1.0/0.5與比例距離的多項式關系

3.3 轉換模型驗證

綜合3.1節(jié)和3.2節(jié)，同時考慮爆高和長徑比的峰值超壓轉換關系式可表示為

在已知某長徑比的圓柱裝藥空爆時測角 α某些比例距離處的峰值超壓，求約束范圍內某工況下各測角處的峰值超壓。首先利用數(shù)值擬合得到三次多項式形式的超壓預測公式，然后利用式（21）實現(xiàn)由已知峰值超壓轉換為所約束范圍內任意長徑比和爆高工況下 β測角處的峰值超壓。利用仿真數(shù)據和實爆結果驗證上述超壓空間轉換模型的可靠性。

基于仿真數(shù)據，在已知長徑比1.0、比例爆高0.5的圓柱裝藥在90°處峰值超壓，求該工況其他測角峰值超壓，并求長徑比1.5、比例爆高1.5的圓柱裝藥各測角處峰值超壓，將所得結果同原值比較，見圖8。

圖8 峰值超壓轉換誤差Fig.8 Conversion errors of peak overpressure

對同一工況其他空間位置峰值超壓的轉換精度較不同工況時稍高，前者平均相對誤差為8.0%，后者為9.6%。兩者在水平方向誤差大于15%的點數(shù)多于其他方向，這是由于水平方向受馬赫反射的影響大，超壓同其他方向的差異明顯。

基于在中部某靶場的同一批次實爆試驗結果，求0°處相同比例距離的峰值超壓，并與試驗數(shù)據對比，如表4所示?？梢?，誤差均在20%以內，轉換效果良好。

表4 模型計算結果對比Table 4 Comparison of model calculation results

需要注意的是，轉換模型構建方法采用了諸多假設，對于近地爆炸沖擊波流場的非線性效應明顯的區(qū)域，模型轉換結果與真實值會存在一定差異。

4 結 論

基于鏡像法，利用與爆高平齊的三波點處特殊幾何關系及馬赫反射終點條件，確定了沖擊波流場分布界限；基于角等分和超壓歸一化值分段線性假設，建立了超壓空間轉換模型，并推廣至圓柱裝藥不同長徑比、當量和爆高的工況。

通過將數(shù)值模型與經驗公式和實爆結果對比，驗證了所用數(shù)值模型的可靠性，同時也驗證了超壓空間轉換模型的合理性和準確性，可為沖擊波超壓毀傷元空間威力評估提供參考。