沉底水雷爆炸威力的數值計算*

姚熊亮，楊文山，陳 娟，張阿漫，初文華

（哈爾濱工程大學船舶工程學院，黑龍江 哈爾濱 150001）

沉底水雷爆炸威力的數值計算*

姚熊亮，楊文山，陳 娟，張阿漫，初文華

（哈爾濱工程大學船舶工程學院，黑龍江 哈爾濱 150001）

采用沖擊波物理中的阻抗匹配鏡像線方法處理物質交界面，利用無網格光滑粒子流體動力學法（smoothed particle hydrodynamics，SPH）方法模擬沉底裝藥的水下爆炸過程，得到了不同方向不同爆距處的峰值壓力反射系數，分析了爆炸流場參量分布和變化規律，研究了不同水底砂層厚度和不同裝藥量對爆炸沖擊波傳播和流場中壓力的影響。結果表明：海底反射系數隨測點與水平線夾角的變化而變化，當測點與水平線夾角為約55°時反射系數最大，其最大值可達1.5。本文的研究可為沉底水雷爆炸載荷的確定提供參考。

爆炸力學；反射系數；SPH方法；沉底水雷；水下爆炸；阻抗匹配

炸藥沉底爆炸是一個涉及邊界效應的復雜爆炸動力學問題，爆炸沖擊波由海水介質入射到海底界面，在界面以下底質中的沖擊阻抗與海水中沖擊阻抗不同，從而產生反射和透射沖擊波或稀疏波，即在水底附近形成復雜的相互作用波系，目前這方面的理論研究還很不成熟［1］。對于沉底水雷等裝藥沉底爆炸，B.V.Zamyshlyayev［2］認為水底附近反射波與入射波相互作用非常復雜，不能用水底反射理論進行分析，一般把它簡化為無限場來處理，等效為1.70倍藥量在無限場爆炸，這樣處理對不靠近水底的測點在壓力峰值上有一定精度，但對水底附近的測點精度很低。金輝等［3］利用Autodyn軟件對裝藥沉底爆炸進行了模擬研究，得到了水底反射對藥包上方范圍起到增強作用，而對水底附近起削弱作用的相關結論，但并未對不同區域處的水底反射系數做進一步的分析。顧文彬等［4－6］利用有限元ALE算法模擬了沉底裝藥爆炸問題。已有的對沉底裝藥爆炸問題的研究均是采用基于網格的數值算法，在涉及高密度比的間斷面求解時，需要進行特殊處理，界面計算過程較復雜。無網格方法可以彌補基于網格的數值方法的一些缺陷，不受網格劃分的限制，可以解決在基于網格的數值方法中由于高壓、高能、大變形等導致網格畸變而計算崩潰的問題。為此，本文中采用阻抗匹配鏡像線方法來處理水－砂界面，研究沖擊波與海底的相互作用機理，利用SPH方法模擬沉底裝藥的水下爆炸過程，以期對認識沉底爆炸流場參量分布和變化規律有一定的幫助。

1 SPH計算方法

水中炸藥的起爆速度和沖擊波的傳播速度很快，可以假設爆炸氣體和水都是無粘性的且絕熱的。將藥包放置在海底砂層附近的沉底接觸爆炸過程中，砂層直接承受強烈的沖擊波作用，此時可認為砂層幾乎喪失傳遞剪切力的能力，呈現出流體的性質，因此可近似地將海底砂質土壤處理成流體。本文中使用Euler方程和狀態方程來模擬整個沉底水下爆炸過程

式中：m、ρ、v、p、e分別為粒子的質量、密度、速度、壓力和內能；i、j代表相互作用的2粒子；α、β表示坐標方向；Πij為Monaghan型［7］人工粘性，可防止求解結果的非物理振蕩和粒子間的非物理穿透；Wij為光滑函數，本文中采用分段五次樣條光滑核函數近似水下爆炸控制方程。

本文中，對爆炸氣體采用JWL狀態方程［8］；對水采用 Mie-Grüneisen狀態方程［9］；對水底砂層采用多項式狀態方程［10］，它是流體物態方程的經驗形式，在壓縮和膨脹狀態下的壓力分別為

式中：a1、a2、a3、b0、b1、b2為物態方程經驗系數；μ 為壓縮比，e為內能。

圖1 沖擊波的透射與反射Fig.1Transmission and reflection of shock wave

圖2 鏡像線近似原理Fig.2 Approximation principle of mirror image

2 阻抗匹配鏡像線方法

水下爆炸過程中，當沖擊波通過2種物質的界面時，會因為2種物質的沖擊阻抗不同而發生反射和透射，從而產生反射沖擊波和反射稀疏波。在SPH方法中，物質交界面處理非常重要，也是一個難點。SPH方法具有拉格朗日性質和粒子性質，在整個演變過程中，來自不同介質的相互接觸的粒子可能會隨著運動而分離，甚至可能不再成為相鄰粒子。目前，大多數利用SPH方法處理水下爆炸的氣－水界面時都采用粒子對粒子的界面算法［9］，能達到很好的效果，也能有效地防止粒子間的非物理穿透。然而將此方法運用到水砂界面的處理時卻遇到了困難，不能簡單地將氣、水粒子間的懲罰力直接施加在水砂粒子上，而關于懲罰力的形式和其中參數的取值目前還未見報道。由于水下爆炸流場壓力都處于20GPa以下，假設整個水下爆炸過程是等熵絕熱過程，所以本文中采用沖擊波物理中的阻抗匹配鏡像線方法［11］來處理水－砂界面，根據界面兩側材料的力學平衡原理和界面連續性條件，由水的沖擊絕熱線計算砂土的沖擊絕熱線。

在阻抗匹配鏡像線法中，利用介質1的沖擊絕熱線的鏡像對稱線來近似地作為等熵卸載線，進而來獲得介質2的阻抗匹配解，如圖1～2所示。

假設ρ01c01＞ρ02c02，界面處產生反射稀疏波。對于C′點，其粒子速度為u1－Δu，由于C點的壓力與C′點的壓力相等，即pC＝pC′，其中pC′＝ρ01［c01＋λ1（u1－Δu）］（u1－Δu），pC＝p2＝ρ02D2（u1＋Δu），則有

式中：ρ01和ρ02分別為介質1和2的初始密度，c01、λ1、c02、λ2為介質1和2的沖擊絕熱系數，D2為介質2中的沖擊波速度。整理式（3）得到關于Δu的二次方程

解此方程得

式中：a＝ρ01λ1，b＝－（ρ01c01＋2ρ01λ1u1＋ρ02D2），c＝ρ01c01u1＋ρ01λ1u21－ρ02D2u1。則介質2在C狀態下的粒子速度和壓力分別為

同理，當ρ01c01＜ρ02c02時，可得到與式（6）相同的阻抗匹配解。

以一維問題為例，說明阻抗匹配鏡像線法在SPH粒子系統中的實施過程，如圖3所示。入射波傳播過程中，當界面處的粒子1滿足

時，就認為入射波到達界面處，式中ε為控制系數，D1為物質1中的沖擊波速度。通過以上阻抗匹配解計算公式，可得到入射波通過界面時粒子2的運動速度u2和壓力p2，此時粒子2的狀態量變為ρ02、u2、p2。然后沖擊波以ρ02、u2、p2為初始條件繼續在介質2的粒子中傳播。這樣，入射沖擊波便完成了界面處的反射和透射。

圖3 2種物質間的運動界面Fig.3 Moving interface between two materials

3 數值計算結果分析

3.1 計算模型

為使問題簡化，假設：（1）水介質為海水；（2）海底砂層均勻平整，由各組分通過理想混合法則組成；（3）砂層下的基巖為固壁。海水計算尺寸為5m×2.5m；TNT方形裝藥質量為36kg，炸藥尺寸為0.15m×0.15m，中心坐標為（0m，0m）；海底砂層計算尺寸為5m×0.5m。整個計算區域離散成24 000個粒子。幾何構形如圖4所示，其中r為等效藥包半徑，R為測點距藥包中心的距離。

圖4 計算模型的幾何構形Fig.4 Geometrical model of example

3.2 水中沖擊波的傳播

圖5為炸藥沉底爆炸的壓力分布云圖。從圖中可清晰地觀察到，炸藥起爆后沖擊波同時向水中和砂層中傳播，由于砂土粒子的沖擊阻抗大于水粒子的沖擊阻抗，所以砂層區域的壓力高于水中的壓力，如圖5（a）所示。隨著沖擊波的不斷向前傳播，在約0.4ms時沖擊波已到達砂層以下的固壁，并向砂層中反射沖擊波，使砂層中的壓力又逐漸增高，如圖5（b）所示。由固壁反射回的反射波傳入水中，所以此后水域流場的壓力變化是一系列爆炸沖擊波和反射沖擊波共同作用的結果。

圖5 不同時刻的壓力分布Fig.5 Pressure distribution at different times

水底砂層土壤和砂底固壁反射的沖擊波和爆炸沖擊波的相互作用使水中沖擊波的峰值壓力和衰減規律與無限水域中有較大的差別，尤其是水底附近區域，如圖6所示。圖中所示為36kg的TNT炸藥沉底爆炸水域中爆距分別為0.45、0.60、0.75、0.90m，與水平線夾角分別為0°、45°、90°的點的壓力時歷曲線。在0°方向上的壓力曲線出現了2個峰值，分別是入射的爆炸沖擊波作用和從砂底固壁反射的沖擊波所致，隨著沖擊波的傳播，在0.75和0.90m處反射沖擊波的峰值壓力甚至超過了直達波；在45°和90°方向上受到反射沖擊波作用微弱得多。

圖6 不同爆距不同方向上的壓力曲線Fig.6 Pressure curves in different directions at different distances from the explosion center

3.3 爆炸流場壓力分布和變化規律

為了研究炸藥沉底爆炸時水底對流場壓力分布的影響和變化規律，提取了在不同爆距處與水平方向每隔10°的考察點的峰值壓力，即與水平方向成θ處的第1次壓力峰值，即pθ，并與自由場中經驗公式的計算峰值pCole進行比較，如圖7所示。由圖中可以看出在沉底爆炸時存在明顯的水底反射現象，且θ在0°～20°范圍內壓力削弱，30°～90°范圍內壓力增強，當θ＞50°時沖擊波壓力峰值增加到一個相對穩定的值，這一結論與文獻［3］一致。

圖7 不同方向不同爆距處的峰值壓力比較Fig.7 Comparisons of peak pressures in different directions at different distances from the explosion center

為了進一步研究砂質水底的壓力反射系數，定義峰值壓力反射系數

由于θ＝0°～10°的流場區域內受砂底固壁影響較嚴重，雙峰現象較明顯，在數值計算的時間范圍內，只考慮水底砂層反射波的作用，即pθ為第1次壓力峰值。圖8給出了不同方向不同爆距處的壓力峰值反射系數γθ的分布情況。由圖中可知，在水底0°附近區域內各測點的壓力峰值比自由場中爆炸的壓力峰值小。這是由于砂質土壤在高壓爆炸產物的作用下，迅速變形成坑，消耗了一部分爆炸能量，引起水底附近沖擊波壓力的削弱。隨著沖擊波在砂層中的傳播，砂層不斷被壓實，會在水－砂界面處反射壓縮波；當到達砂底的固壁時也會反射回沖擊波作用于水域流場中，所以水底附近區域水介質受到不同性質的波系的影響，相互作用也較復雜。隨著θ的不斷增大，γθ也不斷增大，且呈現出一定的規律。在50°～60°附近，γθ達到最大，并隨θ的繼續增大，水底砂層和砂底固壁的反射作用逐漸減弱，所以γθ也隨之減小。此外，在一定角度θ上，γθ也隨著爆距的增大呈現增大的趨勢，在R／r＝26時，即爆炸的中遠場處，各個角度θ上的γθ值基本處于1.4～1.6之間，與目前工程計算中所取的壓力峰值反射系數基本一致。

圖8 不同方向不同爆距處的反射系數Fig.8 Reflection coefficient in different directions at different distances from the explosion center

3.4 不同砂層厚度對計算結果的影響

為了研究沉底裝藥爆炸時砂層厚度對水域流場壓力的影響，選取了36kg炸藥在不同厚度的砂層上爆炸的計算模型，得到0.45m處0°、45°、90°方向上的壓力曲線（見圖9），其中δ為砂層厚度。從圖中可以看出，砂層厚度對水域介質壓力的變化的影響較大，尤其是砂層較薄（δ＝0.1，0.2m）時砂層吸收爆炸能量較少，反射沖擊波的能力也較弱，此時砂底固壁的反射作用是壓力曲線出現雙峰的主要原因，并對第2個壓力峰值作主要貢獻。隨著δ增加到1.0m后，在1.2ms內，再增加砂層厚度對此處的壓力變化幾乎沒有影響，壓力曲線幾乎完全重合。

圖9 不同砂層厚度不同方向在0.45m處的壓力曲線Fig.9 Pressure curves at at 0.45min the case of different sand thickness and directions

3.5 不同裝藥質量對計算結果的影響

本文中還研究了不同裝藥質量對沉底爆炸沖擊波壓力變化規律的影響。表1中給出了沖擊因子C2＝／R 時0°、45°、90°方向上的峰值壓力。由表1能明顯地看出各個方向上壓力峰值相差較大，且隨著藥量的增加而增大，此時自由場中藥量和爆距之間的相似定律在沉底爆炸時已不再成立，必須探索新的相似規律。

表1 不同裝藥量的峰值壓力比較Table 1 Comparisons of peak pressures for different TNT mass

4 結 論

利用無網格SPH算法對存在砂質土壤的海底沉底爆炸問題進行了數值模擬，水下爆炸時沖擊波與不同物質界面的相互作用而產生的反射和透射現象在計算中基本得到了正確的再現。得到了沉底爆炸時不同方向的峰值壓力反射系數γθ及其變化規律，在0°～5°范圍附近γθ＜1.0；在50°～60°附近γθ達到最大，并隨著θ的增大，γθ隨之減小；在中遠場處，各個角度θ上的γθ基本處于1.4～1.6之間。此外不同砂層厚度和裝藥量對計算結果的影響也較大，尤其是不同裝藥量時需要進一步探索新的相似規律。

［1］蔣國巖，金輝，李兵，等.水下爆炸研究現狀及發展方向展望［J］.科技導報，2009，27（9）：87-91.

JIANG Guo-yan，JIN Hui，LI Bing，et al.Review of studies on underwater explosion［J］.Science ＆ Technology Review，2009，27（9）：87-91.

［2］Zamyshlyayev B V.Dynamic loading in underwater explosion［R］.AD-757183，1973.

［3］金輝，張慶明，高春生，等.裝藥水下沉底爆炸壓力場特性研究［J］.科技導報，2009，27（14）：32-37.

JIN Hui，ZHANG Qing-ming，GAO Chun-sheng，et al.Characteristics of pressure field in ground explosion［J］.Science ＆ Technology Review，2009，27（14）：32-37.

［4］顧文彬，馬海洋，唐勇，等.水底對淺水中裝藥爆炸效果的影響［J］.爆破，2003，20（4）：88-92.

GU Wen-bin，MA Hai-yang，TANG Yong，et al.Influence of water bottom on the explosion effect of shallow-layer water charging［J］.Blasting，2003，20（4）：88-92.

［5］顧文彬，葉序雙，張朋祥，等.淺層水中爆炸水底影響的試驗研究［J］.解放軍理工大學學報，2001，2（2）：55-58.

GU Wen-bin，YE Xu-shuang，ZHANG Peng-xiang，et al.Experimental studies of bottom influence in shallow-layer water explosion［J］.Journal of PLA University of Science and Technology，2001，2（2）：55-58.

［6］顧文彬，陽天海，葉序雙，等.單個裝藥淺層水中沉底爆炸的數值模擬［J］.解放軍理工大學學報，2000，1（3）：41-45.

GU Wen-bin，YANG Tian-hai，YE Xu-shuang，et al.Numerical simulation of single charge explosion on seabed in shallow water［J］.Journal of PLA University of Science and Technology，2000，1（3）：41-45.

［7］Monaghan J J，Lattanzio J C.A refined particle method for astrophysical problems［J］.Astronomy and Astrophysics，1985，149（1）：135-143.

［8］Liu M B，Liu G R，Zong Z，et al.Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology［J］.Computers and Fluids，2003，32（3）：305-322

［9］Liu M B，Liu G R，Lam K Y.Investigations into water mitigations using a meshless particle method［J］.Shock Waves，2002，12（3）：181-195

［10］王繼海.多項式形式物態方程及其等熵線［J］.爆炸與沖擊，1992，12（1）：1-10

WANG Ji-hai.Polynomial form of Mie-Grüneisen equation of state and its isentropic［J］.Explosion and Shock Waves，1992，12（1）：1-10.

［11］譚華.實驗沖擊波物理導引［M］.北京：國防工業出版社，2005：1-121

Numerical calculation of explosion power of mines lying on seabed＊

YAO Xiong-liang，YANG Wen-shan，CHEN Juan，ZHANG A-man，CHU Wen-hua
（College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，Heilongjiang，China）

The underwater explosion process of mines lying on seabed was simulated by using the meshfree SPH method.And the substance interfaces were treated by adopting the impedance matching method.The reflection coefficients of the peak pressures in different directions at different distances from the explosion center were obtained.The distribution and change of the explosion field parameters were analyzed.And the effects of different sand thickness and different charge mass on shock wave propagation and pressure were discussed.Simulation results show that the seabed reflection coefficient varies with the angle between measuring points and horizontal line and it arrives at the maximum value of 1.5when the angle is about 55degree.It is helpful for determining the explosion load of mines lying on seabed.

mechanics of explosion；reflection coefficient；SPH method；mines lying on seabed；underwater explosion；impedance matching

21August 2009；Revised 4October 2010

YANG Wen-shan，wenshanyang85＠gmail.com

（責任編輯 張凌云）

O382.4；U661 國標學科代碼：130·3520

A

1001-1455（2011）02-0173-06＊

2009-08-21；

2010-10-04

國家自然科學基金項目（50939002，10976008，50809018）

姚熊亮（1963— ），男，教授，博士生導師。

Supported by the National Natural Science Foundation of China（50939002，10976008，50809018）