近場水下爆炸氣泡射流載荷沖擊船體外板的動響應分析

姜忠濤， 李 燁， 龐學佳， 王詩平

(1. 哈爾濱工程大學 船舶工程學院， 哈爾濱 150001； 2. 湖南磁浮技術研究中心有限公司， 長沙 410000；3. 中國船舶重工集團公司第七0三研究所, 哈爾濱 150078)

船體結構在遭受近場水下爆炸氣泡射流的沖擊作用時，首當其沖的應當是其船體外板，船體板架根據不同的水下工況自然會產生不同的彈性或者塑性響應。關于氣泡載荷[1]及氣泡射流[2]的研究通常從實驗和理論兩個方面進行， Brunton等[3-6]針對不同類型的材料進行了相應的研究，認為射流載荷下的結構變形是由平板上的沖擊流體產生的“水錘”壓力及高速橫向射流對材料表面進行腐蝕和沖刷作用而產生，引起易變形材料結構表面的塑性凹陷及脆性材料的剪切破壞。

關于船體外板在水下爆炸作用下的動響應多從理論、試驗及數值仿真等方面進行研究。Taylor等[7-9]從理論上研究了板承受水下爆炸載荷的結構響；吳成等[10](2007)從理論出發對固支方板剛塑性動態響應進行了分析，并給出相應的運動控制方程和結構變形解析解；張振華等[11]對水下爆炸沖擊波作用下艦局部板架結構的動響應的相似規律進行研究；牟金磊等[12]提出了水下爆炸載荷作用下加筋板結構的計算方法。在試驗研究方面，Ming等[13]對水下接觸爆炸作用下加筋平板的毀傷特性開展了試驗研究，并將試驗結果與SPH (Smoothed Particle Hydrodynamics)計算結果進行對比；崔杰等[14]針對實尺度艙段模型進行了水下爆炸實驗，并對氣泡射流載荷的影響范圍進行分析；Gifford等[15](1988)結合試驗對存在裂紋的焊接厚板進行了近場爆炸分析，并分析了其動態響應。在數值模擬部分，近場爆炸數值模擬方法大多借助于任意歐拉-拉格朗日算法(Arbitrary Lagrangian-Eulerian， ALE) 和耦合歐拉-拉格朗日算法[16]( Coupling Eulerian-Lagrangian， CEL)。王龍侃等[17]采用LS-DGF-DG方法對船體板架結構在近場水下爆炸載荷作用下的毀傷特性進行研究。在具體到內部氣泡載荷的分析過程而言，在計算量允許的前提下邊界附近氣泡的運動形態更多的是采用邊界元法[18-19]進行模擬。將氣泡射流進行幾何等效后借助邊界元法計算氣泡射流的相關幾何數值特征，最后結合CEL或ALE對水射流沖擊結構的響應進行模擬計算[20]。

本文采用ABAQUS軟件中的CEL算法對氣泡射流沖擊簡單平板的過程進行數值模擬以規避材料的高應變率問題，單元大變形等很多常規問題。

1 基本理論和方法

1.1 耦合歐拉-拉格朗日算法

ABAQUS/Explicit可以實現將劃分網格后的歐拉域和拉格朗日域裝配到同一個計算模型中。其中射流對應的水域采用歐拉網格進行離散，而平板或者板架結構則遵循傳統的拉格朗日網格離散方式。兩種網格之間可以交叉，并且可以約束材料防止其進入歐拉域。在歐拉域和拉格朗日域的交界處，后者的網格為前者的網格提供了幾何流動邊界，而前者的網格為后者網格提供了壓力邊界。

在ABAQUS中，耦合歐拉拉格朗日算法(CEL)的原理基于中心差分法。求解過程中用位移的形式來表示相應的加速度和速度，其各自的表達式以及運動方程分別如式(1)和式(2)所示

(1)

(2)

將t時刻的速度和加速度結合式(1)與式(2)即可得到中心差分表達式

(3)

式中:M為質量陣;C為阻尼陣;K為阻尼陣；Q為載荷。

(4)

1.2 射流載荷針對CEL算法的等效原則

由于氣泡射流速度非常快，且氣泡凹陷形成的射流近似為圓臺形狀，為了便于工程應用，本文在采用CEL算法施加射流載荷時將氣泡射流圓臺形狀近似為圓柱，如圖1所示。當氣泡以橫剖圖顯示時，S表示射流形狀橫剖面面積，L表示氣泡射流輪廓的高度，D表示射流形狀近似后圓柱的直徑。表達式為

(5)

圖1 氣泡射流形狀等效方法

在給定初始條件的前提下，采用程序結合邊界元算法可直接確定特定工況下圓柱形氣泡射流的速度、寬度以及高度，為后續CEL計算射流沖擊提供基本參數。

2 數值結果及討論

2.1 有效性驗證

本文采用ABAQUS/Explicit中的CEL算法對射流載荷的相關機理和載荷特性進行分析，勢必需要考慮其計算的精度問題。本小節在此以Obara等[21]的射流試驗數據以及射流載荷的經驗公式為基礎，通過在同樣的工況條件下與CEL算法的計算結果進行對比驗證數值算法的準確性。

Obara等的實驗采用高速水射流沖擊一定距離下的玻璃鋼材料板，通過在特定位置設置壓力感應裝置獲取壓力載荷，同時結合高速攝影裝置捕捉實驗過程中射流沖擊產生的沖擊波、稀疏波現象。高速水射流通過裝置可將射流速度控制在90～1 000 m/s內。試驗中射流沖擊速度為570 m/s。靶板板材規格為50 mm×50 mm×5.9 mm。

針對CEL算法的數值模擬而言，模型如圖2所示。為了與上述試驗條件保持一致，平板結構采用實體單元以拉格朗日網格進行離散，拉格朗日域大小為50 mm×50 mm×5.9 mm，水域以歐拉域代替，所建域底部與拉格朗日域重合，高度為20 mm。射流載荷用圓柱體代替，其高度為12 mm，射流區域的直徑為3 mm，以570 m/s的初始速度沖向平板，并得到結構和射流流體的動響應。

圖2 計算模型示意圖

圖3 射流域在耦合模型中的位置

圖4 板架有限元模型

經驗公式、CEL耦合算法以及試驗結果的壓力時歷曲線對比如圖5所示。其中試驗結果峰值大小為850 MPa，CEL算法結果的峰值大小為817 MPa，三條曲線的峰值相差較小。同時，實驗值在駐壓段基本維持在180 MPa附近， CEL算法與經驗值和試驗值相差較小，能夠滿足相應工程使用精度。

圖5 三種算法下射流載荷中心點壓力曲線對比

2.2 基于CEL的射流載荷沖擊鋼板的動響應分析

2.2.1 數值模型

板架材料采用彈性鋼板，屈服極限為235 MPa，密度為7 850 kg/m3，彈性模量為2.1×1011Pa，泊松比為0.3，采用拉格朗日網格進行離散。射流部分的密度為1 000 kg/m3，采用圓柱形水柱對射流進行替代，圓柱體直徑10 cm，高度50 cm，水柱沖擊速度為200 m/s，水的狀態方程為線性Us-Up方程，水中聲速為1 500 m/s。歐拉域尺寸為156 cm×156 cm×70 cm，整個水域采用154 560個歐拉單元進行離散。耦合模型及平板模型分別如圖6和圖7所示。

圖6 歐拉域與平板耦合模型

圖7 加筋平板模型

2.2.2 射流的載荷壓力特性

當設定板的厚度為10 mm時，在200 m/s的速度沖擊下板架僅發生彈性變形，板架大部分區域應力最大值不超過235 MPa，沒有出現塑性變形，計算結果云圖如圖8所示，且加強筋未出現明顯的應力分布，僅僅是加筋結構與面板交接中心區域出現應力集中。而當速度設定為1 000 m/s時在射流中心區域出現部分塑性變形，如圖9所示，其局部應力已經超過235 MPa，但因射流半徑較小并未達到塑性破壞所需應變，因此沒有出現破口。

圖8 未出現塑性應變時板架應力云圖

圖9 出現塑性應變時板架應力云圖

射流載荷在水下對艦船結構進行沖擊時射流區域的沖擊壓力的分布一直以來都是業界十分關心的問題。因此在射流半徑范圍內流體與板架接觸的區域選定部分參考點提取其壓力進行直觀分析，測點位置示意圖如圖10所示，圖中V0表示射流初始速度，Vs表示水流受到平板約束后橫向流動后產生的橫向速度。

圖10 射流中心附近測點示意圖

根據圖10所示沿徑向在射流半徑范圍內均勻布點，提取其壓力時歷曲線如圖11所示，其中R表示射流半徑。根據文獻[22]計算射流載荷的經驗公式可知，在該速度下“水錘”壓力和水動壓力經驗公式值分別為1.5 GPa和0.5 GPa，從圖11可知數值結果的“水錘”壓力峰值在1.65 GPa左右，t=0.15 ms后出現壓力衰減，隨后進入水動壓力階段，圖中壓力值在0.5 GPa附近振動，與經驗公式結果基本吻合。

同時從圖11可以看出，當射流以圓柱形沖擊平板時，射流半徑范圍內的外圍流體出現壓力爬升的時間比中心附近流體較早，即射流中心產生壓力的時間較晚。同時，對比水動壓力段可以發現，在射流作用的水平壓力輸出端，射流中心區域附近的壓力明顯高于外圍接觸部分，射流中心處的水動壓力峰值近似為半徑處(r=R)附近區域的2倍。另外，且當水柱對平板作用結束后整個半徑區域內的壓力迅速衰減到零。

圖11 射流半徑范圍內典型徑向接觸面節點壓力時歷曲線

2.2.3 射流的流體速度分布特性

以上述水柱沖擊平板為例，選定V0=200 m/s，在流體和結構接觸的表面上的部分歐拉單元中設置速度輸出參考點，分別以射流中心點為中心沿著徑向分散，測點位置與圖10相同，并以射流半徑為無量綱標準量進行無量綱化，各個參考點的速度時間曲線如圖12所示。

圖12 射流中心附近歐拉域各徑向測點速度時歷曲線

從圖12中曲線可以看出，在射流與結構接觸后，測點速度在極短時間內達到最大值后開始衰減，速度衰減過程中會出現少量的振蕩，這種現象可能與平板出現彈性振動有關，但振蕩幅值不大。對比射流半徑范圍內的三個參考點可知，隨著遠離射流中心點，流體最終的穩定速度也會逐漸變大。而對于射流半徑以外的流體，其能達到的速度峰值明顯大于射流沖擊速度。對于r=1.36R的歐拉節點，其最大速度達到663.78 m/s，甚至超過射流速度的3倍，這一現象與Brunton的研究結果相吻合。隨著進一步遠離射流中心點，流體速度峰值有所下降，但依然大于原始沖擊速度，且隨著遠離射流中心點距離較大，流體最終趨于穩態流動，速度保持在180 m/s附近，因此距離射流中心較遠的節點的最終穩定速度保持一致。射流半徑以內的流體最終速度明顯小于半徑以外區域的速度。

為了直觀表示不同射流沖擊時刻射流區域及其附近的流體速度在徑向的分布情況，在流體徑向擴展的方向上設置一系列參考點，將同一時刻下的流體在徑向的分布情況如圖13所示給出。其中橫坐標表示距離射流中心點的無量綱距離，以射流半徑表示參考量，由上圖可知，在t=0.3 ms時射流中心速度有所增加但中心附近半徑范圍內的流體速度出現衰減，這一現象與圖12的結果向吻合，此時最大速度出現在r=1.5R的環形區域附近，但在最大速度以外區域速度則全部為零。隨著時間進一步推進，當t=0.6～0.9 ms時最大速度區域已經位于r=3R附近，但隨著遠離射流中心，速度峰值出現明顯的衰減，且這一階段射流半徑范圍內的速度基本保持穩定狀態。當t=1.2 ms以后雖然速度進一步向邊緣傳播，但在r≤3R范圍內速度基本保持不變，這一結果也與圖12所示結果相同。此后，因為載荷的逐步衰減以及流體的逐漸穩定，速度響應也基本不再出現變化。

圖13 不同時刻射流中心附近徑向流體速度分布

2.2.4 射流引起的結構剪應力分布特性

當高速射流流體在受到結構板架的阻礙作用而產生橫向流動時會產生比射流原始沖擊速度更大的速度，甚至一度能達到沖擊速度的3～4倍，在結構上出現如此高速的流動勢必會在結構表面上產生極大的剪切應力。圖14所示為在1 000 m/s射流沖擊下平板沿Y方向的剪應力分布云圖，由圖可知剪應力梯度分別沿從y軸旋轉45°方向各自以加強筋為界分成四個區域擴展。其中y軸正向45°為負值，負向45°為正值。由于板架邊界設定為簡支，當載荷壓力結束后，剪應力最大值及應力集中現象最終位于平板四個角落處。

(a) t=0.173 ms

(b) t=0.36 ms

(c) t=0.76 ms

(d) t=1.5 ms

為了更加細致地分析射流沖擊平板結構時結構上的剪應力分布，如圖15所示在Y軸正向和XY之間的45°方向以射流半徑為單位選取一系列參考點，并對比參考點上的剪應力大小。由圖16可知45°方向上的剪應力明顯大于Y軸正向上應力，且在射流作用的初始階段會出現一個比較大的瞬時剪應力峰值，且該峰值的出現時間緊隨載荷的峰值壓力而出現。在峰值過后剪應力則在一個較小的應力值附近振蕩，而Y軸正方向剪應力則沒有瞬時峰值的出現，即整個沖擊過程剪應力都很小。

通過對比不同徑向位置剪應力分布可知，r=2R和r=3R位置處45°方向剪應力峰值明顯大于r=R處的峰值，且平穩振蕩階段的剪應力也大于半徑范圍處節點。由此可見，在射流沖擊結構后不僅產生速度極高的橫向射流，同時在射流徑向的某些方位上會對結構產生極高的剪切應力，不同方位的應力之間會產生2～3倍的差別，且射流半徑范圍外的流體由于其速度更高，產生的剪應力也更高，其對結構表面的材料可能造成十分嚴重的剝蝕現象。因此，針對水下爆炸氣泡射流產生的毀傷效應，不僅需要關注其垂向變形和破口，同時也需要重點關注橫向剪切剝蝕毀傷。

圖15 方位示意圖

(a) r=R

(b) r=2R

(b) r=3R

3 結 論

本文以射流載荷沖擊船體外板結構為背景，采用CEL算法對水下爆炸氣泡射流載荷沖擊簡單平板的過程進行了數值模擬。將計算結果與試驗進行了對比并分析了相關規律，本文結論具有一般性，具體如下：

(1) 采用CEL算法計算了射流載荷沖擊PMMA平板的過程，結合試驗和經驗公式對比了其載荷曲線，結果表明CEL算法可有效地模擬水射流對平板的沖擊過程。

(2) 采用CEL算法模擬了射流載荷沖擊彈性鋼板的過程，對射流半徑范圍內結構上的接觸載荷壓力進行了分析，發現載荷壓力在0.6R(R為射流半徑)范圍內基本保持一定，在0.6R以外區域出現衰減，在r=0.8R附近壓力出現較小值。

(3) 通過分析射流半徑范圍內速度分布發現，流體速度最大峰值出現在半徑以外區域(約r=1.36R)，通常能達到射流沖擊速度的3～4倍，且超過一定范圍后速度峰值隨遠離射流中心而逐漸變小，但最終速度趨于一致。

(4) 對平板結構上的剪應力分布進行分析，發現在與平板骨材交叉方向存在較大剪應力，約為與骨材平行方向大小的2～3倍，即射流的存在對平板平面形成巨大的剪切應力和剝蝕效應。

參 考 文 獻

[1] ZHANG A M, CUI P, CUI J, et al. Experimental study on bubble dynamics subject to buoyancy[J]. Journal of Fluid Mechanics, 2015, 776: 137-160.

[2] FIELD J E, LESSER M B. On the mechanics of high speed liquid jets[C]∥Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. London:The Royal Society, 1977: 143-162.

[3] BRUNTON J H. High speed liquid impact[J]. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1966, 260(1110): 79-85.

[4] BOWDEN F P, BRUNTON J H. The deformation of solids by liquid impact at supersonic speeds[C]∥Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. London: The Royal Society, 1961: 433-450.

[5] BOWDEN F P, FIELD J E. The brittle fracture of solids by liquid impact, by solid impact, and by shock[C]∥Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. London: The Royal Society, 1964: 331-352.

[6] SMITH D G, KINSLOW R. Pressure due to high-velocity impact of a water jet[J]. Experimental Mechanics, 1976, 16(1): 21-25.

[7] TAYLOR G I. The pressure and impulse of submarine explosion waves on plates[J]. The Scientific Papers of GI Taylor, 1963, 3: 287-303.

[8] COX A D, MORLAND L W. Dynamic plastic deformations of simply-supported square plates[J]. Journal of the Mechanics and Physics of Solids, 1959, 7(4): 229-241.

[9] HUANG H. Transient bending of a large elastic plate by an incident spherical pressure wave[J]. Journal of Applied Mechanics, 1974, 41(3): 772-776.

[10] 吳成, 倪艷光, 郭磊, 等. 水下爆炸載荷作用下氣背固支方板的動態響應分析[J]. 北京理工大學學報, 2007, 27(3): 205-209.

WU Cheng, NI Yanguang, GUO Lei, et al. Dynamic response of rectangular air-back plate to underwater explosion loading[J]. Transaction of Beijing Institute of Technology, 2007, 27(3): 205-209.

[11] 張振華，陳平毅，漆萬鵬, 等.艦船局部板架結構在水下爆炸沖擊波下動態響應的相似率研究[J]. 振動與沖擊, 2008, 27(6): 81-86.

ZHNAG Zhenhua, CHEN Pingyi, QI Wanpeng, et al. Scaling law of dynamic response of stiffened plates for a ship subjected to underwater shock[J]. Journal of Vibration and Shock, 2008, 27(6): 81-86.

[12] 牟金磊，朱錫，張振華，等.爆炸沖擊作用下加筋板結構變形研究[J]. 海軍工程大學學報, 2007, 19(6): 12-19.

MU Jinlei, ZHU Xi, ZHANG Zhenhua, et al. A study on deformation of blast-loaded stiffened plates[J].Journal of Naval University of Engineering, 2007, 19(6): 12-19.

[13] MING F R, ZHANG A M, XUE Y Z, et al. Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions[J]. Ocean Engineering, 2016, 117: 359-382.

[14] 崔杰，張阿漫，郭君, 等.艙段結構在氣泡射流作用下的毀傷效果[J]. 爆炸與沖擊, 2012, 32(4): 355-361.

CUI Jie, ZHANG Aman, GUO Jun, et al. Damage effect of cabin structure subjected to bubble jet[J]. Explosive and Shock Waves, 2012, 32(4): 355-361.

[15] GIFFORD L N, CARLBERG J R, WIGGS A J, et al. Explosive testing of full thickness precracked weldments[C]∥Fracture Mechanics Twenty First Symposium, ASTM-STP-1074. Philadelphia: ASTM， 1990: 157-77.

[16] 王耀輝, 陳海龍, 岳永威, 等. 水下接觸爆炸作用下的船體板架結構毀傷研究[J]. 中國艦船研究, 2012, 7(4): 76-82.

WANG Yaohui, CHEN Hailong, YUE Yongwei, et al. Damage study on ship plate frame subjected to the underwater contact explosion[J]. Chinese Journal of Ship Research, 2012, 7(4): 76-82.

[17] 王龍侃，張之凡，郎濟才.基于LS-DGF-DG方法的船體板架結構近場水下爆炸毀傷特性研究[J].振動與沖擊,2016,35(16):64-71.

WANG Longkan,ZHANG Zhifan,LANG Jicai,et al.Damage characteristics of hull grillage subjected to near-field underwater explosion based on an LS-DGF-DG method[J].Journal of Vibration and Shock,2016,35(16):64-71.

[18] ZHANG A M, LIU Y L. Improved three-dimensional bubble dynamics model based on boundary element method[J]. Journal of Computational Physics, 2015, 294: 208-223.

[19] ZHANG A M, LI S. CUI J. Study on splitting of a toroidal bubble near a rigid boundary[J]. Physics of Fluids, 2015, 27(6): 809-822.

[20] KLASEBOER E, KHOO B C, HUNG K C. Dynamics of an oscillating bubble near a floating structure[J]. Journal of Fluids and Structures, 2005, 21(4): 395-412.

[21] OBARA T, BOURNE N K, FIELD J E. Liquid-jet impact on liquid and solid surfaces[J]. Wear, 1995, 186(95): 388-394.

[22] COOK S S. Erosion by water-hammer[C]∥Proceedings of the Royal Society of London A: Containing Papers of a Mathematical and Physical Character. London: The Royal Society, 1928: 481-488.