郭 婭，唐文勇，鄭紹文，李德聰，蘇怡然

（1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.中國艦船研究設計中心，武漢 430064）

水下近場爆炸載荷主要包括爆炸前期沖擊波載荷與爆炸后期氣泡載荷。局部結構響應主要由二者所致。為此，沖擊波載荷、氣泡載荷及與結構相互作用已成研究熱點。數(shù)值模擬方法有基于雙漸近近似DAA方法［1］及水下沖擊波分析程序［2］，后者對水下爆炸領域影響深遠。文獻［3－4］對水下爆炸沖擊波載荷及對結構毀傷進行研究。而對氣泡的研究，自Reyleigh建立不可壓縮流場中球形氣泡運動方程以來氣泡潰滅理論不斷完善與發(fā)展，張阿漫等［5－7］基于邊界積分法對氣泡脈動過程進行過深入研究。

理想的水下近場爆炸數(shù)值模擬需同時考慮結構受沖擊波作用的流固耦合效應及受射流沖擊的流固耦合效應；但因軟件功能限制尚不能模擬沖擊波與氣泡射流兩階段完整的流固耦合作用。通常為單獨考慮沖擊波載荷及射流載荷分別進行流固耦合模擬，再將兩計算結果疊加［8］。而實際結構在沖擊波后階段已發(fā)生變形或破損，射流載荷實際作用于已變形或破損結構上，因而對結構造成進一步?jīng)_擊。因此，簡單疊加并不能真實反映結構響應結果。

為真實模擬整個水下近場爆炸過程，本文用MSC.Dytran軟件據(jù)給定爆炸初始參數(shù)，利用基于歐拉法的氣泡動力學程序模擬氣泡近壁脈動過程［9］，獲得射流基本參數(shù)；基于能量等效原則將射流對結構沖擊簡化為等效射流載荷。通過編制MSC.Dytran軟件子程序，在沖擊波數(shù)值模擬階段后自動添加等效射流載荷，從而實現(xiàn)結構受沖擊波－射流聯(lián)合作用的動態(tài)響應數(shù)值模擬。在整個計算過程中沖擊波、射流載荷對結構影響連貫統(tǒng)一，計算結果會更符合實際情況。

1 水下近場爆炸沖擊波及射流數(shù)值模擬

1.1 水下近場爆炸載荷作用過程

水下近場爆炸過程中按時間順序結構先后會受三種載荷，即沖擊波載荷、氣泡脈動壓力及氣泡射流沖擊載荷。沖擊波載荷壓力峰值高、作用時間短，結構受沖擊波作用后會在短時間內產(chǎn)生巨大速度、加速度，導致受損。沖擊波作用后氣泡周圍水體受氣泡內壓力影響產(chǎn)生壓力場，形成氣泡脈動壓力。氣泡脈動后期受壁面影響，在氣泡潰滅時會形成高速射流，對已損傷結構造成進一步破壞。據(jù)研究［10］，氣泡脈動載荷的最大壓力小于沖擊波壓力的10%～20%；據(jù)試驗、數(shù)值仿真研究知，氣泡脈動載荷主要引起中遠場水下爆炸時結構整體鞭狀運動，對近場爆炸的局部損傷影響較小，故本文數(shù)值模擬水下近場爆炸時僅考慮沖擊波載荷及射流載荷。

1.2 沖擊波階段數(shù)值模擬

MSC.Dytran求解流固耦合問題時需定義固體或流體壓力及密度與比內能間之關系，即狀態(tài)方程。常規(guī)水下爆炸情況用TAIT狀態(tài)方程描述水體為

水下爆炸，水與爆炸產(chǎn)物氣體將在歐拉網(wǎng)格中共存。為描述該現(xiàn)象，應使用多物質歐拉單元。炸藥及流體的初始分布、條件參數(shù)可通過多物質歐拉單元初始狀態(tài)定義，見表1。流場邊界均定義為水介質可流出的無反射邊界。

表1 沖擊波數(shù)值模擬模型歐拉網(wǎng)格初始狀態(tài)Tab.1 Initial conditions for the Euler mesh in blast wave simulation model

1.3 射流階段數(shù)值模擬

求解水下氣泡脈非球狀脈動過程的數(shù)值算法主要有基于勢流理論的邊界積分法、基于NS方程的歐拉法。邊界積分法計算效率高，但適用范圍有限；歐拉法適用范圍廣，但計算效率低。本文采用改進的歐拉法模擬氣泡近壁脈動過程［9］。該方法在歐拉法基礎上提出，可有效提高計算速度的近場壓力邊界條件。將近場壓力邊界外部流體設為理想流體，通過勢流理論求解。外部流體為內部流體提供動態(tài)壓力邊界條件，實現(xiàn)歐拉域減小、有效提高歐拉法模擬水下氣泡動態(tài)特性的計算效率。

圖1 射流時刻氣泡形狀簡化描述Fig.1 The simplification of the bubble shape at jet time

據(jù)給定的爆炸初始參數(shù)，采用改進的歐拉法編制氣泡動力學程序模擬氣泡近壁脈動過程，獲得氣泡脈動過程中流場壓力分布及氣泡射流基本參數(shù)。與射流沖擊載荷計算相關參數(shù)為射流形狀、射流速度及射流沖擊時刻氣泡離結構距離?？紤]射流形狀、尺寸參數(shù)時可將射流簡化為凸臺［8］，用3個尺寸描述射流形狀，見圖1。圖中H為凸臺高度；d1，d2分別為凸臺上、下表面直徑；Dist為射流距壁面距離；Vel為射流產(chǎn)生沖擊時刻速度。

2 結構受沖擊波－射流聯(lián)合作用數(shù)值模擬方法

考慮沖擊波載荷與射流載荷聯(lián)合作用時進行沖擊波作用于結構的流固耦合模擬。據(jù)1.3節(jié)所得射流基本參數(shù)進行單獨射流沖擊流固耦合數(shù)值模擬，獲得結構不同位置受沖擊載荷的時歷曲線及結構變形能W0。由于沖擊壓力時歷曲線波動性較大，故需適當簡化方可用于Msc.Dytran二次開發(fā)加載。壓力曲線簡化方法為

設在射流沖擊時對特定時間步及板格其等效射流沖擊壓力為恒定值，式（4）中壓力簡化主要思想為在時間、面積上對壓力取平均，對給定時間區(qū)間Ti及面積區(qū)域（環(huán)形）Sj，等效射流壓力記Pi，j。式（4）積分計算為

獲得等效射流載荷后施加于結構進行流固耦合分析，獲得結構變形能W1。將W1與射流水柱流固耦合模型所得變形能W0進行對比，兩次所得變形能應保持一致。若兩者相差較大，則修正等效射流載荷大小，直至W1與W0相等，此時等效射流載荷對結構作用效果與射流水柱直接作用效果應一致。

計算模型以沖擊波流固耦合模型為基礎，通過調用Msc.Dytran用戶子程序完成結構等效射流沖擊載荷施加。計算過程中可通過用戶子程序標記失效單元標號，統(tǒng)計失效區(qū)域面積。計算流程見圖2。

圖2 結構受沖擊波－射流聯(lián)合作用的數(shù)值模擬流程Fig.2 Process for the numerical simulation of structure subjected to both explosive wave and jet impact

Msc.Dytran用戶子程序用FORTRAN語言編寫，其中EXFUNC函數(shù)允許用戶定義隨時間變化的外部函數(shù)，可用于描述隨時間變化載荷。Msc.Dytran通過TABLEEX選項卡在每個時步內多次調用EXFUNC函數(shù)對結構加載。同一模型可對結構不同受載區(qū)域定義不同TABLEEX選項卡，該選項卡通過不同標識名區(qū)別彼此。在用戶子程序中通過判斷結構單元應變是否超過許用應變篩選破損單元，并統(tǒng)計破損區(qū)域及塑性區(qū)域面積。

3 沖擊波－射流聯(lián)合作用數(shù)值模擬實例

3.1 沖擊波－射流聯(lián)合作用結構響應數(shù)值模擬

圖3 板架結構有限元模型Fig.3 The frame structure finite element model

據(jù)上節(jié)計算流程計算典型舷側板架結構遭受水下接觸爆炸作用后的動態(tài)響應。在水深5 m、距板架5 m處裝500 kg當量TNT炸藥，距離參數(shù)（爆距與氣泡最大半徑比值）為0.5。板架單層結構長寬均6 m，板厚12 mm，長寬方向各設用5根T材⊥加固，見圖3。結構單元網(wǎng)格尺寸100 mm×100 mm。板結構彈性模量210 GPa，泊松比0.3，密度7.8 g／cm3，屈服強度440 MPa；材料破壞準則采用最大等效應變失效準則，即等效塑性應變超過0.3時結構發(fā)生破壞［3］，或取 0.25［12］，本文用 0.25作為臨界失效應變。應變率強化模型用CowperSymond模型［4］，即

據(jù)爆炸初始參數(shù)用氣泡動力學程序對氣泡脈動過程進行數(shù)值模擬見圖4，獲得射流基本參數(shù)見表2。由表2知，射流時刻氣泡貼在結構壁面上直接對結構造成沖擊，射流速度120 m／s。經(jīng)比較，射流參數(shù)與經(jīng)驗、試驗數(shù)據(jù)規(guī)律相符［13－16］。

表2 水下爆炸近壁射流數(shù)值模擬計算結果Tab.2 Numerical results of the nearwall jet caused by underwater explosion

圖4 結構受射流沖擊作用的數(shù)值模型Fig.4 Numerical model for the simulation of the structure under jet impact load

射流流固耦合模型耦合面由結構單元與啞元構成。用多物質歐拉單元，歐拉區(qū)域初始狀態(tài)見表3。

區(qū)域 狀態(tài)方程密度／（kg·m－3）比內能／（J·kg－1）優(yōu)先級全歐拉域 TAIT 1000 8.4×1041射流水柱 TAIT 1000 8.4×1042

在板架結構距中心3 m內由中心向外圍均勻取30個監(jiān)測點進行壓力監(jiān)測，記為P1～P30。部分監(jiān)測點壓力時歷曲線見圖5，可見壓力變化規(guī)律與經(jīng)驗相符［8，15］。

圖5 板架結構監(jiān)測點壓力時歷曲線Fig.5 Time history graph of the pressure in several monitoring points

圖6 等效射流載荷加載區(qū)域示意圖Fig.6 Equivalent jet load loading area

圖7 水下接觸爆炸作用下板架結構吸能曲線Fig.7 Time evolution of structural energy for the shell structure

射流沖擊時間取0.2 ms，以0.5 m為半徑，最終得板架結構由中心向外5層圓環(huán)區(qū)域等效射流載荷分別為：第一層81 MPa，第二層73 MPa，第三層65 MPa，第四層4.8 MPa，第五層 0.07 MPa，其它區(qū)域載荷較小，可忽略，見圖6。

獲得等效射流沖擊壓力隨時間、空間分布后，用上節(jié)用戶子程序方法將該載荷數(shù)據(jù)添加到?jīng)_擊波流固耦合模型中，實現(xiàn)結構受沖擊波－射流聯(lián)合作用的數(shù)值模擬。結構受沖擊波作用的流固耦合數(shù)值模型與射流作用的流固耦合數(shù)值模型類似，耦合面亦為由外殼拉格朗日單元及啞元組合的閉合面，啞元不提供剛度。

3.2 數(shù)值模擬計算結果

用MSC.Dytran軟件計算3.1節(jié)數(shù)值模型，獲得板架結構在水下近場爆炸載荷作用下的動力響應結果。板架結構在沖擊波－射流聯(lián)合作用下的結構變形能及動能隨時間變化曲線見圖7，爆炸過程中結構等效塑性應變見圖8。由圖7看出，結構變形能有兩次大幅度提升，分別為沖擊波作用階段與射流沖擊階段，且結構在沖擊波作用階段吸收的變形能累積到射流階段；結構在射流階段吸收的變形能小于沖擊波作用階段；結構動能經(jīng)歷兩次先增大后減小過程：第一次結構動能增加因沖擊波作用所致，且在阻尼作用下逐漸恢復靜止；第二次動能增加因射流作用所致，且射流載荷對結構沖擊小于沖擊波作用。由圖8可知，沖擊波作用的板架結構等效塑性應變最大為0.157，射流沖擊作用的板架等效塑性應變?yōu)?.175，因此若不考慮射流載荷影響，對結構評估結果偏危險。

由計算結果知，板架結構受沖擊波瞬間作用而獲得較大動能。雖沖擊波作用時間極短，但沖擊作用后10 ms內板架動能持續(xù)轉化為結構變形能。圖8（a）為結構受沖擊波作用后板架的等效塑性應變云圖，板架結構以小板格為單位發(fā)生局部彎曲，塑性變形在邊界附近略高，此由結構邊界約束偏剛性導致。射流沖擊載荷的實際作用時間約在爆炸后1 s。由計算知沖擊波作用后約15 ms時結構塑性變形能已趨于穩(wěn)定，動能趨于零，可認為結構破壞與變形已充分發(fā)展。為節(jié)約計算時間可將射流沖擊作用提前至沖擊波作用后15 ms。板架結構受射流沖擊作用后響應過程與沖擊波階段相似，即板架獲得較大動能后逐漸轉化為變形能。板架受射流沖擊中心變形、應變進一步增大。圖8（b）～（d）為板架受射流沖擊作用后的變形與破壞形式。

圖8 水下近場爆炸作用板架結構等效塑性應變云圖Fig.8 Fringe figure of the effective plastic strain for the shell structure of underwater nearfield explosion

4 結 論

本文利用基于歐拉法的氣泡動力學程序與Msc.Dytran軟件對結構受沖擊波－射流聯(lián)合作用的動態(tài)過程進行數(shù)值模擬，通過算例對計算結果進行分析，結論如下：

（1）近場水下爆炸時射流沖擊載荷不可忽略，射流沖擊會使結構變形、應變進一步加劇，進而導致受沖擊波作用未損壞結構在射流作用下發(fā)生破壞。

（2）本文綜合考慮沖擊波載荷、射流載荷作用，簡化處理射流載荷，使其在沖擊波流固耦合模擬后期自動添加等效射流沖擊載荷，能較真實模擬水下近場爆炸過程，適合工程應用。

［1］Geers T L.Doubly asymptotic approximations for transient motions of submerged structures［J］.The Journal of the Acoustical Society of America，1978，64（5）：1500－1508.

［2］DeRuntz Jr J A.The underwater shock analysis code and its applications［C］.／／60th Shock and Vibration Symposium，SAVIAC，1989：89－107.

［3］陳永念，譚家華.固支方形板在水下爆炸沖擊波載荷作用下動態(tài)響應數(shù)值仿真研究［J］.艦船科學技術，2008，30（1）：131－135.CHEN Yongnian，TAN Jiahua.Dynamical response of square plates under shock load［J］.Journal of Ship Science and Technology，2008，30（1）：131－135.

［4］袁建紅，朱錫，張振華.水下爆炸載荷數(shù)值模擬方法［J］.艦船科學技術，2011，33（9）：18－24.YUAN Jianhong，ZHU Xi，ZHANG Zhenhua.Numerical simulation method study of underwater explosion load［J］.Journal of Ship Science and Technology，2011，33（9）：18－24.

［5］張阿漫.水下爆炸氣泡三維動態(tài)特性研究［D］.哈爾濱：哈爾濱工程大學，2006.

［6］姚熊亮，汪玉，張阿漫.水下爆炸氣泡動力學［M］.哈爾濱：哈爾濱工程大學出版社，2012.

［7］宗智，何亮，孫龍泉.水下爆炸氣泡對水面艦船載荷的數(shù)值研究［J］.船舶力學，2008，12（5）：733－739.ZONG Zhi，HE Liang，SUN Longquan.Numerical study of loading on the ship near an underwater explosion bubble［J］.Journal of Ship Mechanics，2008，12（5）：733－739.

［8］崔杰，李世銘，黃超，等.射流沖擊模型在水下爆炸實驗室中的應用研究［J］.高壓物理學報，2012，26（5）：523－530.CUI Jie，LI Shiming，HUANG Chao，et al.Applied research of jet impact model in underwater explosion experiment［J］.Chinese Journal of High Pressure Physics，2012，26（5）：523－530.

［9］蘇怡然.水下接觸爆炸作用下雙層防護結構瞬態(tài)動力分析方法研究［D］.上海：上海交通大學，2013.

［10］姚熊亮，陳建平，任慧龍.水下爆炸氣泡脈動壓力下艦船動態(tài)響應分析［J］.哈爾濱工程大學學報，2000，21（1）：1－5.YAO Xiongliang，CHEN Jianping，REN Huilong.The analysis of dynamic response of ship hull subjected to gas bubble impulsive pressure of underwater explosions［J］.Journal of Harbin Engineering University，2000，21（1）：1－5.

［11］Popinet S，Zaleski S.Bubble collapse near a solid boundary：a numerical study of the influence of viscosity［J］.Journal of Fluid Mechanics，2002，464：137－163.

［12］李海濤，朱錫，張振華，等.彈塑性邊界下水下爆炸氣泡動態(tài)特性的數(shù)值仿真研究［J］.系統(tǒng)仿真學報，2009，21（14）：4532－4535.LI Haitao，ZHU Xi，ZHANG Zhenhua，et al.Numerical simulation study of pulsation dynamics of underwater explosion bubble located under an elastoplastic plane［J］.Journal of System Simulation，2009，21（14）：4532－4535.

［13］張阿漫，王詩平，白兆宏，等.不同環(huán)境下氣泡脈動特性實驗研究［J］.力學學報，2011，43（1）：71－82.ZHANG Aman，WANG Shiping，BAI Zhaohong，et al.Experimental study on bubble pulse features under different circumstances［J］.Chinese journal of Theoretical and Applied Mechanics，2011，43（1）：71－82.

［14］汪斌，張遠平，王彥平.水中爆炸氣泡與水底邊界相互作用的水射流現(xiàn)象［J］.爆炸與沖擊，2011，31（3）：251－256.WANG Bin，ZHANG Yuanping，WANG Yanping.Water jet phenomena induced by the interaction between underwater explosion bubbles and water bottom boundaries［J］.Explosion and Shock Waves，2011，31（3）：251－256.

［15］張偉，岳永威，張阿漫，等.基于AUTODYN的氣泡與固定壁面相互作用數(shù)值模擬［J］.中國艦船研究，2012，7（6）：24－32.ZHANG Wei，YUE Yongwei， ZHANG Aman， et al.Numerical simulation of underwater explosion bubble interactions with a solid boundary based on Autodyn［J］.Chinese Journal of Ship Research，2012，7（6）：24－32.

［16］張阿漫，姚熊亮.水深和藥量的變化對水下爆炸氣泡射流的影響研究［J］.工程力學，2008，25（3）：222－230.ZHANG Aman，YAO Xiongliang.The effect of charge and water depth on the underwater explosion bubble［J］.Engineering Mechanics，2008，25（3）：222－230.

［17］梅志遠，朱錫，侯海量.雙層防爆艙壁結構動態(tài)吸能特性分析［A］.無錫：船舶力學學術會議暨《船舶力學》創(chuàng)刊十周年紀念學術會議論文集，2007.