約束空間內爆超壓載荷影響因素二維數值模擬

徐維錚，吳衛國

1武漢理工大學高性能艦船技術教育部重點實驗室，湖北武漢430063

2武漢理工大學交通學院，湖北武漢430063

0 引 言

海洋是支撐未來發展的戰略制高點，已成為世界各國競相發展的領域，導致各國間的海洋權益摩擦、領海和島嶼主權爭端不斷。大型水面艦船是海上國防力量建設最為關鍵的部分之一，然而其在執行作戰任務時，不可避免地面臨著日益先進的反艦武器的打擊，包括水下威脅（魚雷、水雷等）和水上威脅（反艦導彈、自殺式飛機等）。

目前，在反艦導彈的研發和制造中，半穿甲爆破型戰斗部被各國廣泛采用，如美國的“魚叉”、法國的“飛魚”、德國的“鸕鶿”以及俄羅斯的“白蛉”等反艦導彈。從攻擊模式和毀傷機理來看，采用半穿甲爆破型戰斗部的反艦導彈距海面3～8 m掠海飛行，憑借彈體的動能侵徹艦船舷側外板，同時觸發延遲引信，穿入艦船艙室內部爆炸［1］。艙內爆炸由于局限在約束空間內，戰斗部爆炸形成的沖擊波經多次壁面反射及疊加，會形成持續時間較長的準靜態壓力，將對艙室內人員、設備及結構產生較大的毀傷。因此，深入研究約束空間內爆炸的超壓載荷及其影響因素具有重要意義。

針對約束空間內爆炸的準靜態超壓載荷問題，國內外學者均有所研究。例如，Anderson等［2］采用量綱分析理論，在大量實驗數據的基礎上，提出了約束空間內爆炸的準靜態超壓載荷計算公式；Feldgun等［3］基于伯努利方程，推導出了泄壓過程中的準靜態超壓公式；徐維錚等［4-6］基于自主開發的三維內爆炸波高精度數值計算程序，探討了泄壓口、WENO格式（Weighted Essentially Non-oscillation Scheme）的精度、炸藥量等因素影響約束空間內爆炸準靜態超壓載荷的規律；Xu等［7］在文獻［4，6］研究的基礎上，進一步通過量綱分析理論給出了計算準靜態超壓載荷的簡化公式。然而，上述研究尚存在一定的缺陷，研究結果與真實的爆炸過程也存在某些差異。

從真實的爆炸過程來看，在載體搭載下戰斗部處于運動狀態，即使其通過艦船舷側穿甲進入艙內，仍保留了一定的剩余速度。從近期國內相關文獻來看，有學者針對空中動態爆炸（動爆）超壓場進行了研究。例如，聶源等［8］研究指出，由于炸藥隨載體運動，故對爆炸后產生的沖擊波超壓進行測試較為復雜，目前主要以地面靜態爆炸（靜爆）實驗測量為主，在動爆與靜爆沖擊波場有著明顯差異的情況下，靜爆實驗無法精確測量出戰斗部的實際威力，所以研究炸藥動爆沖擊波超壓場及其計算模型非常有必要。鑒于此，為了獲得球形炸藥動爆沖擊波超壓場的計算模型，在計算靜爆沖擊波超壓的Baker公式中加入修正因子進行修正，建立了包含炸藥運動速度、對比距離和方位角的修正因子函數的方法；蔣海燕等［9］利用AUTODYN商用軟件，對炸藥在不同運動速度下的空中爆炸（空爆）沖擊波場進行數值計算，以定量研究運動狀態下炸藥空爆沖擊波場的特性，分析了動爆與靜爆沖擊波場的關聯特性，建立了動爆沖擊波超壓工程模型，其計算結果與動爆實驗及仿真結果吻合較好。然而，對于炸藥在運動狀態下艙內爆炸的準靜態超壓載荷的影響規律鮮有研究。

值得注意的是，導彈戰斗部殼體破裂后，會對爆轟產物陣面產生一定的初始擾動，因此在進行等效裸藥數值模擬過程中，需要給予初始界面一定的隨機擾動，文獻［10-11］在球形炸藥爆炸數值模擬過程中，均在初始界面處添加了一定的隨機擾動，使數值計算更符合實際，但有關添加的初始擾動對內爆炸超壓載荷的影響規律仍需進一步探討。此外，在數值模擬過程中，有關網格尺寸對內爆炸超壓載荷的影響規律也需考察。

綜上所述，在當前的數值模擬研究中，對于炸藥運動、爆炸初期界面的初始擾動及網格尺寸這些因素對內爆炸超壓載荷的影響規律關注尚少，因此，本文擬通過自主開發的高精度二維數值計算程序對上述3種因素進行初步探討，為真實爆炸過程中的抗爆設計和毀傷評估提供一定的指導及依據。

1 控制方程及數值求解

基于自主開發的約束空間內爆炸波高精度數值計算程序，通過求解二維可壓縮歐拉方程［12］，分別采用五階WENO有限差分格式和三階TVD龍格—庫塔（TVD-RK）法，對其空間項及時間項進行數值離散［13］。

2 內爆炸載荷數值計算

本節主要采用自主開發的計算程序開展不同內爆炸工況下的數值模擬研究，探討炸藥運動、炸藥界面初始擾動及網格尺寸對艙內爆炸超壓載荷的影響規律。

2.1 炸藥運動的影響

為了探討炸藥運動對艙內爆炸超壓載荷的影響，假設炸藥自左向右沿艙室長度方向運動，并在艙內中部瞬時起爆。在數值模擬中，炸藥運動速度v0分別取為0，600和1 200 m/s。

2.1.1 艙室尺寸及爆炸初場

圖1所示為數值模擬計算用的長方形內爆炸艙室示意圖。本文在瞬時爆轟假定［14］的基礎上，將初始的方形炸藥等效為均勻的高壓爆轟產物，其邊長為160 mm，密度為1 630 kg/m3，壓力為3.057 9×109Pa。周圍的空氣密度為 1.0 kg/m3，壓力為1.0×105Pa。艙室壁面邊界條件設置為剛性邊界（不考慮艙室結構變形），泄壓口處的邊界設置為透射邊界［13］。計算網格數為180×80。在艙室右壁面中心處設置一個典型的測點1，以輸出內爆炸超壓載荷的時程曲線。

2.1.2 爆炸壓力場分布

圖2所示為典型時刻下3種運動速度的炸藥在密閉艙內爆炸的壓力場分布云圖。由圖可以看出，靜爆壓力場分布沿艙室長度方向呈對稱分布，炸藥運動狀態改變了艙內壓力場的分布形態，破壞了對稱性；而對于動爆壓力場，炸藥運動狀態改變了艙內壓力分布形態。

2.1.3 爆炸載荷分析

圖3所示為不同運動速度的炸藥在艙室（密閉艙和泄壓艙）內爆炸后壁面測點1處的超壓載荷時程曲線。由圖可以看出，炸藥運動速度對沖擊波峰值及到達時刻的影響較大，隨著炸藥運動速度的增大，沖擊波峰值明顯增大，初始沖擊波峰值到達時刻提前。

為定量分析圖3中靜爆和動爆沖擊波的差異性，并驗證所開發數值計算程序對動爆沖擊波計算的可靠性，參考文獻［9］，利用AUTODYN商用程序對不同運動速度下炸藥空中爆炸的沖擊波場進行數值計算，在綜合分析動爆與靜爆沖擊波場的關聯特性的基礎上，建立了動爆沖擊波超壓的工程計算模型。動爆沖擊波超壓峰值pp，d與炸藥運動速度v0、靜爆沖擊波超壓峰值pp，s的關系表達式為

式中：R為沖擊波陣面至爆心的距離；c0為波陣面前的空氣聲速；α為沖擊波陣面與炸藥爆心的夾角。

由式（1）可知，增大炸藥運動速度，將增大向右方向運動的沖擊波強度。針對本節算例，在3種炸藥運動速度下，分別測得內爆炸壁面測點1的初始反射沖擊波峰值為184.4，327.8，485.1 MPa。根據強沖擊波壁面正反射關系pr/pi=(3γ-1) /(γ-1)（其中pi為入射沖擊波峰值，pr為反射沖擊波峰值，γ為氣體絕熱指數），分別得到入射沖擊波峰值為 23.05，40.98，60.64 MPa。由上式，分別得到600，1 200 m/s運動速度下炸藥的動爆初始入射沖擊波峰值為37.68，55.90 MPa，其數值解與擬合公式的相對誤差分別為-8.05%，-7.82%。這初步驗證了使用本文數值計算程序計算動爆沖擊波的可靠性。

根據能量守恒定律，推導出考慮炸藥運動初速的密閉空間內爆準靜態超壓峰值Δps的計算公式為［6］

式中：m為炸藥質量；V為艙室體積；，為爆炸過程中釋放的總能量（含運動狀態下炸藥動能），其中QTNT為炸藥爆熱，QTNT=4.69×106J/kg；p0為初始大氣壓力；ρE為炸藥密度，ρE=1 630 kg/m3。

由式（2）可知，增大炸藥運動速度將使艙內最終形成的內爆炸準靜態壓力增加，炸藥初始動能最終轉化為艙內氣體內能。由式（2），得3種運動速度下對應的內爆炸準靜態超壓峰值分別為60，62，69 MPa，最大增加幅值約15%。然而，由于圖3中爆炸初期前幾個沖擊波峰值很大，使得準靜態超壓峰值間的差異不能很好地體現。為此，本文單獨給出了3種運動速度下炸藥在密閉艙內爆炸的準靜態超壓峰值與理論解的對比，如圖4所示，從圖中可以明顯看出兩者的差異。

2.2 界面初始擾動的影響

為了探討界面初始擾動對艙內爆炸超壓載荷的影響，本文選取了圓形炸藥進行數值計算研究。在保證炸藥質量一定的前提下，圓形炸藥界面初始擾動采用如下方式添加［15］：

式中：r為圓形炸藥的半徑；r0為無擾動時炸藥的半徑；ξ為擾動幅值；N為擾動波數；θ為界面上一點與圓心的夾角。選取2種擾動波數（8和16）進行數值計算。本算例中，ξ的取值為10 mm。

2.2.1 爆炸初場

本文在瞬時爆轟假定［14］的基礎上，將初始圓形炸藥等效為均勻高壓爆轟產物，其半徑為80 mm，密度為 1 630 kg/m3，壓力為 3.057 9×109Pa。計算的艙室尺寸、炸藥位置、周圍空氣參數、邊界條件及測點設置與3.1.1節相同。對于爆炸初場計算，網格數取為360×160。在該算例中，針對泄壓艙，泄壓口設置在艙室上壁面中心處，尺寸取為200 mm。圖5所示為不同界面擾動下的炸藥初始形態。

2.2.2 爆炸壓力場分布

圖6所示為典型時刻下不同初始界面擾動時炸藥在泄壓艙內爆炸的壓力場分布云圖。從圖可以看出，界面初始擾動主要改變了艙內壓力分布的局部形態，對艙內整體壓力分布形態的影響較小。

2.2.3 爆炸載荷分析

圖7和圖8分別給出了不同初始界面擾動下密閉、泄壓艙壁面典型測點1處的超壓載荷時程曲線。從圖中可以看出，不同初始界面擾動下超壓時程曲線具有相似的變化規律，主要差異在于沖擊波峰值的大小及相位。由于求解的方程為雙曲型強非線性方程，即使是微小的初始界面擾動，均可使得爆炸沖擊波的峰值與相位產生差別。這說明界面初始擾動主要影響瞬態沖擊波載荷，而對于持續時間較長的準靜態超壓載荷的影響很小。

從圖中還可以看出，當初始擾動波數N=8時，超壓峰值數值較大；當初始擾動波數N=16時，超壓峰值與無界面擾動結果的接近程度增加。這說明，隨著擾動波數的增大，界面上的總體擾動趨于均勻化，使得計算結果與無擾動時的更接近。若增加初始擾動波數至N=32，計算得到的局部放大圖，如圖7（b）和圖 8（b）所示。從這兩個圖來看，同樣可說明上述影響規律。

由上述分析可知，對于某爆炸工況，存在某一特定的擾動波數N可給出較高沖擊波峰值的計算結果。

2.3 網格尺寸的影響

這里需要說明的是，在數值計算程序中，網格對炸藥填充質量的影響較大，在分析網格對爆炸載荷的影響規律時，需要保證炸藥質量一定。鑒此，本文在保證炸藥質量一定的前提下，選取了3種網格數（80×80，160×160，320×320）進行數值計算。

2.3.1 艙室尺寸及爆炸初場

由于網格數增大會使計算時間迅速增加，所以在計算艙室的爆炸初場時，將艙室選為正方形，如圖9所示。在瞬時爆轟假定［14］的基礎上，本文將正方形炸藥等效為均勻高壓氣團，其邊長為20 mm，密度為1 630 kg/m3，壓力為3.057 9×109Pa。炸藥位置、周圍空氣參數、邊界條件及測點設置與3.1.1節相同。

2.3.2 爆炸載荷分析

圖10所示為不同網格尺寸下密閉艙和泄壓艙壁面典型測點1處的超壓載荷時程曲線。

從圖10可以看出，隨著網格尺寸的減小，沖擊波峰值增大（尤其是初始沖擊波峰值和第2次沖擊波峰值），到達時間提前，不過網格尺寸對形成的準靜態壓力載荷影響很小。

為了定量分析沖擊波峰值隨網格數變化的規律，圖11給出了壁面測點1處的超壓載荷時程局部放大圖和初始及第2次沖擊波峰值隨網格數的變化曲線。

這里需要說明的是，當網格數取為640×640時，計算非常緩慢，若計算到6 ms，需消耗大量時間。因此，圖11中只在對比初始和第2次沖擊波峰值時給出了網格數為640×640的計算結果，且僅計算到t=0.20 ms。

分析圖11可知，艙壁面反射沖擊波載荷的計算與網格尺寸相關性很大，隨著網格加密，反射沖擊波載荷迅速增加，且收斂速度較慢，不僅如此，計算時間迅速增加。因此，在數值計算時不能將網格劃分得過細。由于沖擊波陣面非常陡峭，在數值計算中會跨越2～3個網格節點，至于采用多大的網格尺寸既能捕捉沖擊波陣面，又能保證峰值不失真，目前還沒有合適的準則來判斷，作者也將在后續研究中對此進一步探討。

3 結論

通過本文的研究，得到如下主要結論：

1）炸藥運動改變了艙內爆炸壓力場的分布形態，增大炸藥運動速度能較明顯地提高沖擊波峰值，但對內爆炸準靜態超壓載荷的影響較小。

2）界面初始擾動主要影響瞬態沖擊波峰值，而對持續時間較長的準靜態超壓載荷的影響很??；對于內爆炸工況，存在某一特定的擾動波數可給出較高沖擊波峰值的計算結果。

3）網格尺寸對沖擊波峰值及到達時間的影響較大，但對形成準靜態壓力載荷的影響很小。隨著網格尺寸的減小，沖擊波峰值增大，到達時間提前。尤其是在進行結構與沖擊波耦合響應計算時，需要專門對網格進行對比分析，以選取合適的網格尺寸。