基于SPH-FEM方法的水下近場爆炸數值模擬研究

姜忠濤， 王　雷， 孫鵬楠， 黃　瀟

(1.哈爾濱工程大學 船舶工程學院，哈爾濱　150001； 2.哈爾濱船舶鍋爐渦輪機研究所，哈爾濱　150001)

?

基于SPH-FEM方法的水下近場爆炸數值模擬研究

姜忠濤1， 王雷2， 孫鵬楠1， 黃瀟1

(1.哈爾濱工程大學 船舶工程學院，哈爾濱150001； 2.哈爾濱船舶鍋爐渦輪機研究所，哈爾濱150001)

摘要：針對結構水下近場爆炸載荷作用響應求解難點，通過改進的三維軸對稱光滑粒子流體動力學方法(Smoothed Particle Hydrodynamics, SPH)計算獲得近場爆炸載荷后傳輸給非線性有限元軟件ABAQUS，利用聲固耦合模型對結構響應進行時域非線性計算，形成預報水下近場爆炸載荷對結構毀傷的SPH-FEM模型，實現從藥包起爆、結構大幅變形、局部撕裂直至完全剪切破壞的全過程模擬，對載荷時歷曲線進行試驗驗證。計算背空矩形鋼板在近場爆炸載荷的響應表明，數值結果與試驗值吻合良好。SPH-FEM模型計算效率高、可操作性強，易推廣至大型復雜結構受水下近場爆炸毀傷的分析與評估。

關鍵詞：艦船工程；矩形鋼板；近場爆炸；SPH-FEM；結構毀傷

水下近場爆炸嚴重威脅海戰的艦船[1]。大當量爆炸會導致艦船局部破壞甚至失去動力，因此對水下近場爆炸研究具有重大意義[2]。

研究近場爆炸作用的結構響應有兩個難點，即爆炸載荷強度確定及結構瞬態強非線性響應計算。較水下爆炸試驗，理論計算、數值模擬因成本低、效率高、周期短頗受關注[3]。傳統方法中對近場爆炸毀傷效果用半經驗公式估算，但只對簡單板架有效[4]，不適用如水面艦船、潛艇等復雜結構形式。而功能強大的流固耦合分析軟件可用于水下近場爆炸結構毀傷預報，如LS-DYNA、AUTODYN等。計算區域較大時兩軟件計算量極大，效率較低。ABAQUS在工程中應用程度較高，主要因其計算效率高，適用于水面艦船、潛艇等復雜結構[5]。而ABAQUS目前只局限于中遠場爆炸計算，原因為傳統方法中采用半經驗公式確定水中爆炸載荷[6]，而半經驗公式只對爆距大于6倍裝藥半徑有效，因此需通過其它有效手段如試驗、數值模擬等獲取近場爆炸載荷。

考慮工程中水下爆炸裝藥形式一般為球形或柱形，符合軸對稱特點[7]。本文用新的無網格軸對稱SPH法[8-9,16-17]計算水中近場爆炸載荷。較純三維SPH法，軸對稱SPH法計算效率極大提高。由于離散區域小，將粒子劃分得很密，可提高計算精度。用SPH方法獲得爆炸載荷后在ABAQUS中將載荷施加到結構表面進行響應計算，形成適用于近場爆炸計算的SPH-FEM方法。

矩形板作為海洋結構物中最普遍的結構形式，研究其在近場爆炸載荷作用下的響應具有重要意義[10-11]。由于矩形板結構形式簡單，試驗數據豐富，可為數值驗證提供基礎。基于SPH-FEM模型，本文計算不同厚度矩形板受近場爆炸載荷響應，并與試驗數據進行對比，驗證SPH-FEM模型的精度及有效性。數值結果分析表明，矩形板在近場爆炸作用下的失效模式主要為塑性大變形、邊緣局部張力撕裂及完全剪切破壞。

1理論模型

1.1三維軸對稱SPH方法

傳統SPH中場函數f(r)的核近似〈f(r)〉[12-13]為

f(r)≈〈f(r)〉=∫Ωf(r′)W(r-r′,h)dr′

(1)

式中：r為所求點矢量坐標；r′為r緊支域Ω內插值點；W(r-r′,h)為光滑核函數。

本文所用三次樣條核函數為

(2)

所求三維問題具有軸對稱特性時可將式(1)在柱坐標中表示，先沿周向θ積分，再沿徑向r及高度z向積分，即

(3)

式中：

(4)

據式(2)光滑核函數定義0≤R<2時，W有意義，因此可推導θ取值范圍為

(5)

(6)

將式(6)代入式(3)，得柱坐標下核近似積分式為

〈f(r,z)〉=∫SWCf(r′,z′)r′dr′dz′

(7)

由式(7)知，三維積分方程已簡化為坐標平面r-z內的二維積分，已極大降低SPH方法計算量，提高效率。由于可在r-z平面內布置較密集粒子，計算精度亦極大提高。因所求問題具有軸對稱性，故f(r,z)在θ方向無變化，只需求解f(r,z)在r軸、z軸方向偏導數，即

(8)

核函數的偏導數為

(9)

本文計算程序中，式(9)用數值積分計算。

1.2控制方程

利用式(8)將N-S方程進行離散，即

(10)

式中：ρ，v，e，m為粒子密度、速度、內能及質量；∏ij為人工粘性。

據質量守恒，mjrj=ρjΔrjΔzjrj在整個計算過程中保持不變。此時仍需引入狀態方程將連續方程、動量方程及能量方程解耦。水的Mie-Grüneisen狀態方程為：

(1)μ>0時水處于壓縮狀態，有

(γ0+aμ)e

(11)

(2)μ<0時水處于膨脹狀態，有

p=ρ0C2μ+(γ0+aμ)e

(12)

式中：ρ0為流體質點初始密度；C為聲速；γ0為Grüneisen系數；S1,S2,S3為材料Hugonoit常數；a為體積修正系數。

Mie-Grüneisen狀態方程中具體參數見表1。

表1　Mie-Grüneisen狀態方程中相關參數

TNT爆轟物采用JWL(Jones-Wilkins-Lee)狀態方程，即

(13)

式中：ei為單位質量內能；η=ρ/ρ0；A，B，R1，R2，ω為由實驗數據擬合的常數，具體見表2，其中E0為初始單位質量內能。

表2　JWL狀態方程中相關參數

1.3變光滑長度

模擬TNT爆轟過程中粒子間距逐漸增大，此時需更新光滑長度，否則會出現緊支域內粒子數量不足而產生巨大離散誤差。本文采用光滑長度隨時間變化公式[14]，即

(14)

式中：d為空間維度。

盡管該軸對稱SPH方法用于模擬三維問題，但核近似只在r-z平面內進行，因此d=2，求出dρ/dt后可利用上一時間步h及ρ求出dh/dt，每個時間步對光滑長度進行更新。

1.4顯示積分方法和對稱軸處理

1.4.1預測校正積分法

本文用預測校正法對離散控制方程進行顯示積分計算。令任意場函數f(r)，預測校正積分法為

(15)

式(15)表示在第n個時間步計算場函數f(r)n變化率Df(r)n/Dt，時間前進Δt/2獲得f(r)n+1/2，完成預測步計算。在此基礎上計算預測步的場函數變化率Df(r)n+1/2/Dt，在f(r)n上前進一整時間步長獲得f(r)n+1，完成校正步計算。本文程序中用恒定時間步長，即Δt=0.06h/C。采用OpenMP并行技術在Inter(R) Core(TM) i7-3770 處理器上實現8線程并行計算，具體并行程序設計見文獻[15]。

1.4.2對稱軸處理

由于對稱軸附近粒子數量較少，離散誤差較大，粒子穿透對稱軸現象難以避免。采用動態邊界法在每個時間步關于對稱軸鏡像分布相應的虛粒子，其參數設置為

(16)

虛粒子參與式(10)的控制方程求解，為防止虛粒子對稱軸附近真實粒子密度、加速度計算產生非物理性影響，設其速度、壓力為0。為防止個別粒子穿透對稱軸導致計算中斷，本文提出反彈算法。即在每個時間步判斷粒子是否穿透對稱軸，將穿透的粒子位置、速度關于對稱軸鏡像到坐標分量r的正向。

1.5SPH-FEM聯合算法

水下爆炸尤其近場爆炸具有瞬態強非線性特點，完成精確實時耦合計算非常困難。文獻[2-3,5]成功將ABAQUS聲固耦合算法用于結構受中遠場水下爆炸載荷響應計算，其中載荷確定用文獻[6]的半經驗公式。本文SPH-FEM聯合算法計算過程見圖1。利用已建三維軸對稱SPH方法模擬裝藥爆轟，其過程通過時域內逐步改變裝藥粒子種類實現；爆轟波傳到水汽界面時高壓氣團進一步膨脹，在水中輻射出沖擊波；數值模型中記錄結構所在位置的壓力載荷時歷曲線。壓力載荷得出后程序轉到ABAQUS，將載荷施加在結構與水域的耦合面，計算結構非線性動態響應。ABAQUS中選散波(scattered wave)公式進行計算，結構響應由入、散射波響應疊加，能充分考慮結構存在對流場中載荷分布影響。ABAQUS中具體操作過程見文獻[3]。

圖1　SPH-FEM聯合算法求解近場爆炸示意圖Fig.1 Sketch of the solving of near-field underwater explosion with SPH-FEM model

結構的彈塑性模型在近場爆炸計算中十分重要，因在此類工況下結構往往產生巨大塑形變形甚至斷裂失效。本文采用線性強化彈塑性模型，即

(17)

式中：E為彈性模量；E1為切線模量；σs為屈服應力；εs為應力等于σs時的應變；σf為抗拉強度。

對鋼材而言，材料的應變率效應不能忽略。本文采用Cowper-Symonds本構方程計算動態屈服應力，即

(18)

式中：D=40.4 s-1；q=5。

2數值結果與討論

2.1變光滑長度對精度影響

軸對稱SPH模型中對稱軸處壓力不穩定性問題最難處理。而用變光滑長度技術能較好改善對稱軸處氣泡膨脹不均勻及沖擊波壓力不穩定問題。SPH測試模型中初始粒子間距Δx=1.458 mm，球形TNT裝藥半徑r0=29.2 mm，半徑方向布置20個粒子，TNT粒子數585。計算水域為球形半徑R=291.6 mm，水域粒子數約6.3萬個。0.15 kg球形TNT藥包在球心起爆0.09 ms時爆炸氣泡形狀及水域沖擊波壓力分布見圖2。由圖2看出，用變光滑長度后球狀藥包從中心引爆后均勻膨脹，對稱軸處沖擊波壓力分布較均勻，因此模擬裝藥水下爆炸可壓縮問題時，變光滑長度能大幅提高模擬結果精度。

圖2　爆炸氣泡形狀對比Fig.2 The comparison of bubble shapes

2.2自由場水下爆炸載荷驗證

為驗證本文三維軸對稱SPH方法的水中爆炸載荷精度，將軸對稱SPH方法模擬自由場水下爆炸沖擊波壓力時歷曲線與試驗數據進行對比。試驗在中國工程物理研究院進行，工況為8 kg球形TNT藥包在自由場爆炸，測量爆距4 m處壓力時歷曲線。SPH模型中初始粒子間距Δx=5.32 mm，球形TNT裝藥半徑r0=106.4 mm，半徑方向布置20個粒子，TNT粒子數618。計算水域為球形，半徑R=6.384 mm，水域粒子約227萬個。裝藥在球心起爆，起爆后不同時刻水域沖擊波壓力云圖見圖3。由圖3看出，TNT爆轟后形成的高壓氣泡沿各方向均勻膨脹，氣泡呈球形，沖擊波壓力云圖在r-z平面內沿各角度分布均勻，對稱軸處未見壓力不穩定現象。

圖3　不同時刻高壓氣泡膨脹及沖擊波壓力云圖Fig.3 The shock wave pressure contour and the expansion of the high-pressure bubble

爆距4 m處沖擊波壓力時歷曲線見圖4。由圖4看出，SPH模擬結果與試驗結果吻合良好，未出現傳統數值算法中壓力波峰非物理性振蕩現象，由此驗證三維軸對稱SPH方法的有效性。

圖4　沖擊波壓力時歷曲線試驗與SPH結果對比Fig.4 The comparison of the time evolution of the shock pressure between experimental data and numerical results

2.3矩形板近場爆炸驗證

表3　矩形板受近場爆炸試驗工況

表4　試驗所用矩形鋼板材料參數

數值模型中用三維軸對稱SPH方法計算工況1～4的沖擊波載荷。忽略氣泡脈動壓力及氣泡水射流影響。獲得沖擊波載荷后用ABAQUS建立背空矩形板模型，用四邊形網格離散，網格邊長5 mm×5 mm，見圖5(a)。流場建立半徑1 m的半球形水域，利用四面體網格離散。矩形板位于球形水域中心。水域與矩形板重疊部分網格種子間距與矩形板相同；水域外圍網格種子間距為0.05 m，見圖5(b)。矩形板與水域施加TIE條件，其余表面施加無反射邊界條件，矩形板四周施加夾支邊界條件。

圖5　ABAQUS中建立的背空矩形板近場爆炸模型Fig.5 The numerical model created in ABAQUS

各工況矩形板失效情況見圖6。左側為試驗結果，右側為數值結果。

圖6　矩形鋼板受近場爆炸載荷作用毀傷情況Fig.6 The damaging of the rectangular plate subjected to near-field underwater explosion

工況1中矩形板厚度4 mm，在40 g裝藥爆轟作用下發生大幅度永久塑形凹陷。矩形板中心凹陷位移U隨時間變化曲線見圖7(a)，可見不考慮應變率效應時矩形板凹陷結果較試驗值偏大70%，考慮矩形板應變率效應的數值結果與試驗值較吻合，稍偏小10%。工況2中矩形板厚度2 mm，在20 g裝藥爆轟作用下亦發生大幅度塑形永久變形。由圖7(b)可知，此時不考慮應變率效應結果仍偏大較多，考慮應變率效應結果偏大5.17%。工況3將裝藥量提高到50 g，由圖6可知，試驗中矩形板長邊中部出現局部張力撕裂，數值結果較好再現出試驗現象。對比圖7(c)矩形板中心位移時間曲線知，考慮應變率效應的數值、試驗結果誤差僅1.87%。由工況1～3數據知，不考慮材料應變率效應時數值結果始終偏大，考慮時相當于增強結構的抗沖擊能力，結構響應變小，與試驗值更接近，因此應變率效應在本構模型中不可忽略。

工況4中將裝藥量提高到80 g，此時矩形板發生嚴重剪切破壞，僅在4個直角處有小塊殘留破片與夾支邊界相連。實驗、數值結果中均可發現，脫落鋼板形狀呈突出狀，即鋼板在巨大沖擊載荷作用下先發生大幅度塑形凹陷后，邊緣發生張力局部撕裂，并在張力、剪切力綜合作用下發生整體脫落。

圖7　矩形鋼板中心凹陷變形隨時間變化曲線Fig.7 Time evolution of the bulge deformation on the center of the plate

表5為試驗值、計算值誤差統計，考慮材料應變率效應的數值結果誤差均保持在10%以內，可見本文SPH-FEM方法能較準確預報矩形板在近場爆炸載荷作用下的響應。而對不同工況，計算值在試驗值上下浮動，并非全部偏大或偏小，可見考慮應變率效應的本構模型是有效的。

由數值結果顯示的鋼板變形隨時間演化過程分析知，鋼板先發生凹陷變形，隨后邊緣達到屈服極限，邊緣張力達到強度極限時出現局部張力撕裂。長邊中點最先達到強度極限斷裂(圖6(c))。隨局部撕裂擴展及剪應力共同作用下矩形板完全脫落，即在近場爆炸載荷作用下的失效模式為隨沖擊因子增加鋼板失效程度逐漸遞增，依次表現為大幅度永久凹陷，邊緣局部張力撕裂及剪切破壞。

表5　試驗值與計算值誤差分析

3結論

(1) 通過采用變光滑長度技術，改進對稱軸處的數值處理，可有效提高三維軸對稱SPH方法預報水下近場爆炸載荷精度，載荷結果與試驗值吻合良好。 將載荷曲線施加到ABAQUS軟件中，計算結構非線性動態響應，形成預報水下近場爆炸載荷對結構毀傷效應的SPH-FEM模型，并基于該模型計算矩形板在近場爆炸載荷作用響應。

(2) 在近場爆炸載荷作用下，矩形板邊緣為薄弱區域。在矩形板凹陷變形過程中四邊受張力最大，其中長邊中點最先達到強度極限。沖擊因子較高時矩形板邊緣在張力、剪切力共同作用下發生剪切破壞，易形成巨大破口，威脅結構物生命力。本文數值結果可為工程結構物抗爆抗沖擊設計、評估提供參考。

參 考 文 獻

[1] Ramajeyathilagam K, Vendhan C P. Deformation and rupture of thin rectangular plates subjected to underwater shock[J]. International Journal of Impact Engineering, 2004,30: 699-719.

[2] 姚熊亮,張阿漫,許維軍. 聲固耦合方法在艦船水下爆炸中的應[J]. 哈爾濱工程大學學報, 2005,26(6): 707-712.

YAO Xiong-liang, ZHANG A-man, XU Wei-jun. Application of coupled acoustic-structure analysisto warship underwater explosion[J].Journal of Harbin Engineering University, 2005,26(6): 707-712.

[3] 宗智,趙延杰，鄒麗. 水下爆炸結構毀傷的數值計算[M]. 北京:科學出版社,2014.

[4] 朱錫,白雪飛,黃若波,等. 船體板架在水下接觸爆炸作用下的破口試驗[J]. 中國造船, 2003,44(1):46-52.

ZHU Xi,BAI Xue-fei,HUANG Ruo-bo,et al. Crevasse experiment research of plate membrance in vessels subjected to underwater contact explosion[J]. Shipbuilding of China, 2003,44(1): 46-52.

[5] 姚熊亮,張阿漫,許維軍,等.基于ABAQUS軟件的艦船水下爆炸研究[J]. 哈爾濱工程大學學報,2006,27(1):37-41.

YAO Xiong-liang, ZHANG A-man, XU Wei-jun, et al. Research on warship underwater explosion with ABAQUS software [J].Journal of Harbin Engineering University, 2006,27(1): 37-41.

[6] Geers T L,Hunter K S. An integrated wave-effects model for an underwater explosion bubble[J]. J. Acoust. Soc. Am., 2002,111(4): 1584-1601.

[7] 李金河,趙繼波,譚多望,等. 炸藥水中爆炸的沖擊波性能[J]. 爆炸與沖擊, 2009,29(2): 172-176.

LI Jin-he,ZHAO Ji-bo,TAN Duo-wang, et al. Underwater shockwave performances of explosive [J]. Explosion and Shock Waves, 2009,29(2): 172-176.

[8] Ming F R, Sun P N, Zhang A M.Investigation on charge parameters of underwater contact explosion based on axisymmetric SPH method [J]. Applied Mathematics and Mechanics, 2014,35(4): 453-468.

[9] Liu M B, Liu G R. Smoothed particle hydrodynamics (SPH): an overview and recent developments [J]. Archives of Computational Methods in Engineering, 2010,17(1): 25-76.

[10] Rajendran R, Paik J K, Lee J M. Of underwater explosion experiments on plane plates[J]. Experimental Techniques, 2007,31(1):18-24.

[11] Rajendran R, Narasimhan K.Performance evaluation of hsla steel subjected to underwater explosion [J]. Journal of Materials Engineering and Performance, 2001,10(1): 66-74.

[12] Liu G R, Liu M B. Smoothed particle hydrodynamics: a meshfree particle method [M]. World Scientific,2003.

[13] Colagrossi A,Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics [J]. Journal of Computational Physics, 2003,191(2): 448-475.

[14] Benz W. Smoothed particle hydrodynamics: a review [J]. NATO Workshop, Les, Arcs, France, 1990,302:269-288.

[15] Zhang A M, Cao X Y, Ming F R, et al. Investigation on a damaged ship model sinking into water based on three dimensional SPH method [J]. Applied Ocean Research, 2013, 42: 24-31.

[16] 明付仁，張阿漫，楊文山，等. 艦船水下接觸爆炸的SPH算法研究[J].振動與沖擊, 2012, 31 (10): 147-151.

MING Fu-ren, ZHANG A-man, YANG Wen-shan, et al. SPH algorithm to deal with the problem of underwater contact explosion of warship [J]. Journal of Vibration and Shock, 2012, 31 (10): 147-151.

[17] 初文華,張阿漫,明付仁,等.SPH-FEM耦合算法在爆炸螺栓解鎖分離過程中的應用[J].振動與沖擊, 2012,31(23): 197-202.

CHU Wen-hua, ZHANG A-man, MING Fu-ren, et al. Application of three-dimensional SPH-FEMcoupling method in unlocking process of an explosion bolt [J]. Journal of Vibration and Shock, 2012, 31 (23): 197-202.

基金項目：國家自然科學基金資助(51279038；51479041)

收稿日期：2014-11-04修改稿收到日期：2014-12-26

中圖分類號:U663；O383

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.02.021

Numerical investigation on near-field underwater explosion using SPH-FEM method

JIANG Zhong-tao1, WANG Lei2, SUN Peng-nan1, HUANG Xiao1

(1. College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China;2. Harbin Marine Boiler and Turbine Research Institute, Harbin 150001, China)

Abstract:The response of structures subjected to underwater shock wave is a highly nonlinear problem, the main difficulties of which are in the determination of the magnitude of shock wave load and the response of structures under transient load. Based on an improved axisymmetric Smoothed Particle Hydrodynamics (SPH) method, the near-field load of the underwater explosion was calculated and then it was used as the input to the FEM package ABAQUS. By using a coupled acoustic-structure model in ABAQUS, the nonlinear response of the structure was obtained. Based on the SPH-FEM model, the whole process of the near-field underwater explosion was simulated, starting from the detonation of the explosive charge to the large deformation, local tearing and complete damaging of the structure. The theoretical explanation of the axisymmetric SPH method was presented in detail and the calculated shock wave pressure load was validated by the experimental data. Then, the responses of rectangular plates subjected to near-field explosion were numerically investigated. The numerical results agree well with the experimental observations. Some useful conclusions regarding near-field explosion were drawn. The SPH-FEM model introduced here is robust and efficient, which is suitable for engineering applications.

Key words:ship engineering; rectangular plate; near-field explosion; SPH-FEM; structural damage

第一作者 姜忠濤 男，博士，高級工程師，1977年生