宗 智，陳 崗，葉 帆，李海濤，趙延杰

（1大連理工大學 工業(yè)裝備結構分析國家重點實驗室 船舶與海洋工程學院，遼寧 大連 116024；2中國船舶與設計研究院，上海 200011）

1 引 言

水下爆炸主要包含兩個階段：沖擊波階段和氣泡脈動階段，但實際上除了上述兩個主要爆炸載荷外，還會有空化效應引起的二次加載對結構的影響。沖擊波在接觸到自由液面及結構濕表面時會發(fā)生反射，使附近水域壓力驟降，當水域壓力下降至水的空化極限壓力時，就會產(chǎn)生空化效應。空化效應具體可以分為區(qū)域空化和局部空化。區(qū)域空化主要是由于自由表面反射引起的反射稀疏波與入射波疊加后形成的負壓區(qū)引起的，而局部空化是由于結構物的反射波與入射波疊加后形成的負壓區(qū)引起的。

針對水下爆炸引起的空化現(xiàn)象，Aron[1]研究了區(qū)域空化的形成，理論推導了區(qū)域空化的上下邊界計算公式。Makinen[2]采用平面波近似法，對比分析了4種不同的空化模型，并對簡單例子分別利用近似方法和經(jīng)典理論求解，比較了計算結果。Rajendran[3]分析了金屬平板在遭受非接觸爆炸時受到空化效應和氣泡脈動的再加載現(xiàn)象，將試驗得到的總體損傷與理論計算得到的初次沖擊波損傷進行了對比分析。李海濤[4-5]利用平面波理論對局部空化和區(qū)域空化的形成特性進行了研究，并通過具體實驗對局部空化理論進行了驗證；理論分析了結構目標尺度的變化對空化區(qū)域形成的影響，兩者符合較好；通過試驗和數(shù)值仿真方法研究了氣泡脈動引起的錐形空化的形成特點，初步分析了錐形空化的形成原因。許斐[6]運用了ABAQUS中內嵌的聲—固耦合法對三塊不同尺寸的氣背圓板進行了水下沖擊波載荷作用下的數(shù)值模擬?？紤]了空背板材料的應變率，模型中網(wǎng)格的密度，水域中的空化和炸藥參數(shù)對數(shù)值結果的影響。葉帆[7]模擬了水下爆炸沖擊波產(chǎn)生的空化效應對船舶等水面結構的響應造成的影響，特別是片空化潰滅時會對結構造成較大的二次加載。

隨著計算機技術的發(fā)展，現(xiàn)在許多著名的商業(yè)有限元計算軟件都包含有類似的計及空化效應的聲學有限元方法。例如NASTRAN-CFA-DAA代碼[8]就采用了類似的氣穴聲學有限元方法。而ABAQUS軟件也融入了一種類似的計及空化效應的聲學有限元方法[9]，其中流體域是采用曲面波近似邊界進行截斷，以流體的壓力作為因變量來判斷是否啟動空化效應的計算。

本文運用大型商業(yè)有限元程序ABAQUS，數(shù)值模擬了水下爆炸產(chǎn)生的空化區(qū)域的膨脹擴大、收縮潰滅，得到了水中測點的壓力曲線，與實驗進行了對比。然后通過空背方板的數(shù)值模擬，捕捉到了局部空化對結構的二次加載現(xiàn)象，并與實驗值進行了對比。

2 ABAQUS水下爆炸空化效應模擬

Abaqus/Explicit計算水下爆炸時可以同時考慮空化效應的影響。Abaqus中，水域由聲學單元來模擬，水節(jié)點的絕對壓力為靜水壓力與計算所得的動壓之和，如果絕對壓力下降到了空化極限壓力，則開始計算空化效應，該部分水介質假定開始自由膨脹而動壓沒有相應的下降。一旦發(fā)生了自由膨脹，水域中的壓力分布將會重新計算?？紤]空化效應后，水介質的壓力可由下式來表述：

其中：pc為空化極限壓力，p0為靜水壓力，pv為水節(jié)點計算所得動壓，Kf為水的體積模量，x為水節(jié)點的空間坐標，uf為水節(jié)點的位移向量，θi表示i個相互獨立的場變量，例如溫度、空氣濕度、水的鹽度等等影響水介質密度和體積拖曳系數(shù)的因素。

要計及空化效應的影響，必須運用總波公式來計算水下爆炸，并且需添加兩個關鍵字。首先需要對水域設置靜水壓力初始條件，此靜水壓力只為計算空化效應提供判斷依據(jù)，而不會影響水域中動壓的計算。選擇所有水域節(jié)點創(chuàng)建一個nodeset，使用關鍵字*INITIAL CONDITIONS,TYPE=ACOUSTIC STATIC PRESSURE，對此nodeset施加靜水壓力。格式如下：

其中：Setname是水域節(jié)點nodeset的名稱，P1是第一個參考點的靜水壓力，X1、Y1、Z1分別是第一個參考點的x、y、z坐標，P2是第二個參考點的靜水壓力，X2、Y2、Z2分別是第二個參考點的x、y、z坐標，根據(jù)兩個參考點的靜水壓力，由程序自動線性計算水域節(jié)點nodeset中的所有節(jié)點的靜水壓力。兩個參考點一般設在沿著重力方向的一條直線上。其次，空化極限壓力由關鍵字*ACOUSTIC MEDIUM,CAVITATION LIMIT來設置，格式如下：

其中：PC—指定的空化極限壓力（Pa）。

3 區(qū)域空化效應

當水下爆炸產(chǎn)生的沖擊波到達水面之后會反射形成反射波，反射波與入射波疊加之后會使水域的絕對壓力驟降，當壓力下降到水的空化壓力時就會出現(xiàn)空化現(xiàn)象。

圖1 水下爆炸沖擊波載荷幾何示意圖Fig.1 Underwater shock attack geometry

圖2 沖擊波壓力的截斷Fig.2 Shock wave pressure profile with cut-off

圖3 區(qū)域空化示意圖Fig.3 Bulk cavitation zone

為了更好地反應沖擊波的反射現(xiàn)象，如圖1所示，實爆點在自由表面的另一側形成一個虛爆點，我們認為反射波壓力由虛爆點產(chǎn)生，由于入射波比反射波先前到達測點，所以入射波會由于壓力衰減而小于隨后到達的反射波壓力，引起水域絕對壓力的急劇下降，表現(xiàn)出截斷現(xiàn)象，一旦絕對壓力到達空化極限壓力，空化就會形成，見圖2。圖3是區(qū)域空化示意圖，區(qū)域空化是以爆點為中心對稱的，因此圖中只顯示了一半。其存在上下邊界，此邊界不是某時刻空化區(qū)域的邊界，而表示的是整個空化過程中會發(fā)生空化的水的區(qū)域。

如圖1所示，爆點與測點之間的距離為

虛爆點與測點之間的距離為

根據(jù)爆炸沖擊理論中經(jīng)典的經(jīng)驗公式，沖擊波壓力可以表示成

其中：t0為沖擊波到達測點處的初始時刻，t≥t0；Pmax為測點處的峰值壓力，單位MPa；W為炸藥當量，單位kg；R為爆距，單位m；K1、K2、A、B為沖擊波常數(shù)，與所使用的炸藥相關；θ為沖擊波衰減常數(shù)。

考慮靜水壓力，包括大氣壓和水壓，記為Ph，Ph=Patm+ρgy。則任意時刻任意測點處的總壓力為

可以根據(jù)Arons[1]的方法來確定區(qū)域空化的上下邊界。對于上邊界，取t=R2/c，令（8）式等于零；對于下邊界，令測點正壓力（包括入射波壓力和靜水壓力）的衰減速度等于反射波壓力的衰減速度。經(jīng)過推導，得到兩個方程：

以上兩個方程可以使用Matlab編譯程序編碼，可以計算任意當量的各種類型炸藥在任意爆炸深度下形成的區(qū)域空化的邊界。

Marcus[10]在他的報告中做了關于區(qū)域空化的實驗，如圖4所示，此實驗采用彭托利特炸藥，藥量30.8 kg，埋深21.3 m，壓力傳感器放置在水面下3.05 m，與爆點的水平距離51.8 m的地方，實驗證明水下爆炸產(chǎn)生的區(qū)域空化在閉合潰滅時會產(chǎn)生較大的沖擊壓力，并提供了相關實驗數(shù)據(jù)。

圖4 實驗示意圖Fig.4 Sketch map of the test

圖5 區(qū)域空化的范圍Fig.5 Bulk cavitation zone

首先應有上述提到的編譯代碼可以計算得到此實驗條件下區(qū)域空化的大致范圍，如圖5所示?？梢钥吹娇栈瘏^(qū)域半徑約290 m，深度約7.5 m，這個數(shù)值對數(shù)值模擬時的建模具有一定的參考價值。由于空化區(qū)域是旋轉體形的，關于爆點中心對稱的，在Abaqus中我們里采用二維模型來模擬，選取有限元模型水域寬度300 m，深10 m，幾何模型及邊界條件設置如圖6所示。模型A，最小單元尺寸0.1 m，共300 000個單元，模型B，最小單元尺寸0.05 m，共1 200 000個單元，為了得到空化區(qū)域的半徑值，取計算時間0.2 s。由于計算時間較短，氣泡脈動載荷可以忽略。本文取水的空化極限為零，即認為當水壓為負壓時就開始計算空化。兩個模型計算得到的空化區(qū)域如圖7和圖8所示，深色區(qū)域是絕對壓力為零的區(qū)域，即是空化區(qū)域，可以看到空化區(qū)域的擴大和縮小，直至潰滅。

從模型A和模型B計算得到的空化區(qū)域隨時間變化的圖示可以看出兩者的趨勢基本一致，但由于兩者的網(wǎng)格密度不同，進而影響到水域壓力分布，所以計算得到的空化區(qū)域稍微有點不同，同時可以看出模型A的空化區(qū)域在154 ms全部潰滅，此時潰滅處的半徑大概249 m，與理論計算值相比，誤差為14.1%，這是由于網(wǎng)格較大，水域中壓力衰減較快，而模型B中的空化區(qū)域在181 ms全部潰滅。此時潰滅處的半徑大概281 m，與理論計算值相比，誤差為3.1%，這是由于圖網(wǎng)格更密，模型B計算得到的空化區(qū)域顯示得更細膩，衰減也較慢。

圖6 邊界條件示意圖Fig.6 Schematic of the fluid boundary conditions

圖7 模型A計算得到的空化區(qū)域Fig.7 Bulk cavitation region calculated by Model A

圖8 模型B計算得到的空化區(qū)域Fig.8 Bulk cavitation region calculated by Model B

圖9 模型A與模型B計算得到的測點壓力時程曲線與實驗對比圖Fig.9 Comparison of the absolute pressure

表1 沖擊波壓力峰值、二次沖擊時刻及壓力峰值統(tǒng)計表Tab.1 Statistics of the shock peak value,reloading peak value and moment

模型A與模型B計算得到的測點壓力時程曲線與實驗測得的壓力時程曲線如圖9所示，模擬得到的測點處壓力時程曲線與實驗曲線相同，都表現(xiàn)出了明顯的二次沖擊壓力，二次沖擊壓力出現(xiàn)的時間與實驗基本一致，但峰值有些差別。模型A計算得到的沖擊波壓力峰值明顯小于實驗值，模型B計算得到的沖擊波壓力峰值接近于實驗值，說明網(wǎng)格越密，沖擊波在傳播過程中衰減得越少。模型A計算得到的二次壓力峰值大于模型B的計算值，但都小于實驗值。得到的持續(xù)時間比實驗值略長。模型A和模型B的二次沖擊持續(xù)時間均比沖擊波持續(xù)時間長，約為后者的2倍。

模型A中測點處進入空化發(fā)生在24.5 ms，二次沖擊發(fā)生在38.9 ms時，模型B中測點處進入空化發(fā)生在24.6 ms時，二次沖擊發(fā)生在38.2 ms時，觀察模型A和模型B在24 ms和38 ms左右測點附近空化區(qū)域變化情況，如圖10-13所示?？梢钥吹剑瑝毫y點附近空化區(qū)域24 ms開始經(jīng)過測點，38 ms發(fā)生了閉合，閉合時刻與二次沖擊時刻吻合良好。

圖10 模型A測點附近區(qū)域進入空化示意圖Fig.10 Schematic of cavitation occur near the measuring point in Model A

圖11 模型B測點附近區(qū)域進入空化示意圖Fig.11 Schematic of cavitation occur near the measuring point in Model B

圖12 模型A測點附近空化區(qū)域閉合示意圖Fig.12 Schematic of cavitation closure near the measuring point in Model A

圖13 模型B測點附近空化區(qū)域閉合示意圖Fig.13 Schematic of cavitation closure near the measuring point in Model B

4 局部空化效應

水下爆炸沖擊波接觸結構物表面后會發(fā)生反射，入射波和反射波疊加會使流固耦合面附近的水壓力下降，當達到水的空化壓力時，就會出現(xiàn)空化現(xiàn)象，這樣的空化稱為局部空化。若結構位于淺水中，則容易發(fā)生局部空化。如果結構位于足夠深度的水中，由于靜水壓力較大，局部空化一般不會發(fā)生。

Hammond[11]進行了一系列水下爆炸實驗，以研究水中空氣背襯鋼板的結構響應。其將方形鋼板密封固定在堅固的空心鋼制方箱上，方箱由繩索連接到浮體而懸浮在水中，如圖14所示。Hammond讓鋼板中心位于水面下2.3 m處，并將炸藥置于遠距離爆炸，則爆炸載荷主要由沖擊波及局部空化效應構成，而避免了自由表面、海底反射、區(qū)域空化效應及氣泡脈動的影響。

其論文中給出了某工況下鋼板中心的壁壓曲線及位移曲線，此工況采用的炸藥為0.164 kg彭托利特，置于水下2.3 m，正對鋼板中心，爆距3.63 m，如圖15所示。目標正方形鋼板邊長0.78 m，板厚 3 mm，鋼板密度 ρ=7 850 kg/m3，楊氏模量 E=2.11×1011Pa，泊松比 μ=0.3。

采用計算區(qū)域空化上下邊界程序，計算得到此工況下的區(qū)域空化范圍，如圖16所示。區(qū)域空化最深不超過2 m，且空化區(qū)域距方箱有一定距離，因此確實可以忽略區(qū)域空化效應的影響。

圖14 方箱布置圖Fig.14 Schematic of the chain arrangement used to secure the box

圖15 炸藥位置示意圖Fig.15 Schematic of the charge position

圖16 此工況下區(qū)域空化范圍示意圖Fig.16 Bulk cavitation zone

此問題具有對稱性，因此可以建立二分之一有限元模型。目標鋼板和方箱用共節(jié)點殼單元建立，單元大小為1.5 cm，共8 112個單元。根據(jù)實驗情況，近似認為目標鋼板四周是剛性固定的，因此在目標鋼板四周施加全位移約束。采用圓柱形水域，水域半徑1.56 m，鋼板附近水單元大小為1.5 cm，向外逐漸變大，共706 992個單元。以目標鋼板中心為坐標原點，x軸水平沿著鋼板表面，y軸指向鋼板的外法向，z軸豎直向上。在對稱面上施加對稱條件，水域其他表面施加無反射邊界條件。有限元模型如圖17所示。爆炸載荷加載時間0.05 s，水的空化極限設為零。

計算得到的鋼板中心位移與實驗測得曲線比較如圖18所示，可見板受到水下爆炸作用后發(fā)生了往復振動，而未產(chǎn)生塑性變形，計算得到的位移變化情況與實驗基本一致。實驗測得的板中心最大凹陷23 mm，計算最大值為15 mm；實驗測得的板中心最大外凸11 mm，計算最大值也為11 mm。

圖17 有限元模型示意圖Fig.17 Schematic of the finite element model

圖18 計算得到的鋼板中心位移與實驗比較圖Fig.18 Comparison of experimental and calculated results for the displacement of the plate centre

圖19 實驗測得的鋼板中心壁壓曲線Fig.19 Experimental pressure curve

圖20 計算得到的鋼板中心壁壓曲線Fig.20 Calculated pressure curve

圖19、20分別為鋼板中心壁壓的實驗測量曲線與計算曲線，兩者的變化趨勢一致，沖擊波到達時壁壓迅速增大，然后很快降為零，說明附近水域發(fā)生了空化現(xiàn)象，然后又受到了二次加載，二次加載壁壓峰值約為沖擊波壁壓峰值的30%，且衰減速度較慢。由表2可以看到，計算得到的壁壓值略小于實驗值，二次加載的間隔時間也略小于實驗值?？傮w來看，數(shù)值模擬的結果與實驗是比較接近的。

表2 沖擊波壁壓峰值、二次加載壁壓峰值及間隔時間的比較Tab.2 Comparison of the shock peak pressure,reloading peak pressure and their interval

圖21顯示的是計算得到的目標鋼板附近水域發(fā)生局部空化的情況，深色區(qū)域為空化區(qū)域?？梢钥吹戒摪逵娓浇虬l(fā)生了空化現(xiàn)象，空化區(qū)域迅速擴大，然后縮小、潰滅。在6.6 ms左右，板附近空化區(qū)域幾乎全部潰滅，與二次加載時刻一致。而由圖18位移曲線可以看到6.6 ms左右，無論是實驗曲線還是計算曲線都顯示鋼板中心仍往y軸反方向移動。因此可以說明二次加載是由于局部空化在流固耦合面附近閉合潰滅引起的，而不是由于鋼板回彈與水碰撞引起的。

圖21 鋼板附近水域局部空化情況Fig.21 Local cavitation region near the plate

實驗沒有測量鋼板中心的速度，計算得到的鋼板中心速度曲線如圖22所示。從圖中可以看到受到?jīng)_擊波作用后鋼板中心速度迅速增大，而后減小、振蕩，6 ms后由于受到二次加載作用而再次大幅增大，然后又減小、振蕩。這說明二次加載會對結構響應產(chǎn)生強烈的影響，二次加載引起的鋼板中心速度峰值甚至大于沖擊波引起的速度峰值。

圖22 計算得到的鋼板中心y向速度曲線Fig.22 Calculated y direction velocity curve of the plate centre

5 結 論

（1）Abaqus可以數(shù)值模擬水下爆炸產(chǎn)生的空化區(qū)域的膨脹，縮小直至潰滅的整個過程，實驗證明在空化區(qū)域閉合潰滅時附近水質點會對結構造成較大的二次沖擊，數(shù)值模擬也驗證了這一結果，并與實驗數(shù)據(jù)基本一致。

（2）水下爆炸沖擊波到達自由液面后會反射引起稀疏波，稀疏波與入射波疊加在水面下形成較大的負壓區(qū)，出現(xiàn)空化區(qū)域，根據(jù)Arons[1]的方法利用matlab編程可以預報區(qū)域空化的范圍，對流體有限元模型尺寸的確立有一定的指導意義。

（3）通過對局部空化模擬得知二次加載不是由于鋼板回彈與水碰撞引起的，而是因為局部空化區(qū)域在流固耦合面附近閉合潰滅，從而對周圍水質點及結構施加較大的二次加載。

（4）區(qū)域空化引起的二次加載相對于沖擊波雖然峰值較小，但是持續(xù)時間大約是沖擊波的兩倍，作用效果顯著。其引起的結構速度響應可以達到?jīng)_擊波相當?shù)牧考?。所以水下爆炸考慮空化效應是必須的。

（5）網(wǎng)格劃分的差異對數(shù)值模擬結果有一定的影響。較密的網(wǎng)格更能捕捉到空化區(qū)域的邊界，更精確地模擬空化效應，所以對于水下爆炸這一瞬時高度非線性問題加密水域網(wǎng)格是必要的。

[1]Arons A B,et al.Long range shock propagation in underwater explosion phenomena II[C].Underwater Explosions Compendium,1943,I.

[2]Makinen.Cavitation models for structures excited by a plane shock wave[J].Journal of Fluids and Structures,1998,12:85-101.

[3]Rajendran.Reloading effects on plane plates subjected to non-contact underwater explosion[J].Journal of Materials Processing Technology,2008,206:275-281.

[4]李海濤,朱 錫,黃曉明,牟金磊.水下爆炸沖擊波作用下空化區(qū)域形成的特性研究[J].高壓物理學報,2008,22(2):181-186.

[5]李海濤,朱 錫,牟金磊,黃曉明.水下近距爆炸作用下彈性鋼板處的空化特性研究[J].海軍工程大學學報,2008,20(1):21-24.

[6]許 斐,周 力,宗 智.鋁合金艦艇在水下沖擊波作用下動態(tài)響應的數(shù)值研究[J].艦船科學技術,2011,33(1):31-40.

[7]葉 帆,宗 智,李海濤.水下爆炸沖擊波產(chǎn)生的區(qū)域空化效應二次加載的并行數(shù)值研究[J].船科學技術,2011,33(9):11-17.

[8]Shin Y S,Santiago L D.Surface ship shock modeling and simulation:Two-dimensional analysis[J].Shock and Vibration,1998,5:129-137.

[9]Coupled acoustic-structural medium analysis[K].ABAQUS Theory Manual.

[10]Marcus M H.The response of a cylindrical shell to bulk cavitation loading[R].NSWC TR 81-295.1983.

[11]Lloyd H,Raphael G.Structural response of submerged air-backed plates by experimental and numerical analyses[J].Shock and Vibration,2000,7:333-341.