基于雙漸進方法的水下爆炸氣泡載荷作用下艦船的動態響應分析

張 弩，宗 智，張文鵬，孫 雷，李章銳

(大連理工大學 工業裝備結構分析國家重點實驗室 船舶工程學院，大連 116024)

艦船在水下爆炸載荷作用下的動態響應是一個非常復雜的問題，需要涉及水下爆炸載荷的模擬，爆炸氣泡的動力學特性，流固耦合分析及水中結構的非線性動態響應等多項研究內容。目前，對于這類三維結構動態問題的主要研究方法是邊界元/有限元結合的方法。主要采用流固耦合方程來表征耦合界面上的載荷，而不需要對水域進行建模，因此減少了大量的計算量。這種方法主要經歷了平面波理論[1]，曲面波理論[2]，虛質量理論[3]到雙漸近(DAA)理論[4－6]的發展過程。特別是由Geers所提出的DAA理論，綜合并改進了前面幾種方法，在高、中、低頻段都有較好的精度。DAA方法經過眾多學者們的改進與發展[7－10]，目前是一種解決流固耦合問題的非常有效的計算方法。

水下爆炸中對艦船主要產生兩種載荷:一是瞬態的沖擊波，其壓力峰值很高但持續時間很短，通常造成水中結構的局部破壞;二是氣泡脈動載荷，氣泡中大約含有47%的爆炸能量，具有非常強的破壞力。這一階段的明顯特征是達到的壓力峰值較沖擊波的壓力峰值低，而且持續時間很長，約為幾百毫秒，會造成水中結構的總體響應和局部響應[11－12]。

對于氣泡的動力學，眾多學者已經進行了大量的研究工作[13－19]。對于氣泡載荷作用下的艦船的動態響應，目前的主要研究方法都是將艦船簡化為船體梁，研究船體梁的總體鞭狀響應[20－26]。然而，對于氣泡作用下三維全船結構的動態響應的細節分析，開展的研究還很少。

本文致力于建立一套有限元方法與DAA方法相結合的計算程序，用以研究水下爆炸氣泡作用下水面艦船的動態響應。首先，建立一個考慮遷移效應，自由面效應和氣泡阻力的氣泡模型和一個水面艦船結構模型，并闡述了流固耦合與結構響應相關的理論分析;然后以該船模作為算例，研究了氣泡作用下船體模型的總體響應和局部響應。比較了不同位置的加速度，速度和位移時程曲線。最后詳細討論了氣泡作用下船體模型動態響應的特征與機理。

1 理論背景

1.1 控制方程

一個彈性結構受到外部激勵的運動方程可以表示為:

其中:Ms，Cs，Ks分別為N×N階的質量矩陣，阻尼矩陣和剛度矩陣。N×1階向量分別為結構的位移，速度和加速度。列向量F代表外力。N為結構的自由度總數量。

對于浸沒在流體的結構所受的外部激勵可以表示為:

式中:pi和ps分別為濕表面上入射流和輻射流作用下產生的節點壓力向量，Af為濕表面單元的面積對角矩陣，G為結構節點力和濕表面節點力的坐標轉換矩陣。

當一個水中結構如水面艦船或潛艇浸沒在無限聲學流體介質中，結構濕表面的殼單元的控制方程可以由DAA方法來表示。對于流體中的彈性結構，DAA方法表示的結構表面流體的運動即為各正交流體邊界模態的線性組合。一階雙漸近DAA方程的矩陣形式表達為:

式中:Mf為N×N階的流體質量矩陣;ρ和c分別為流體密度和流體中聲速;us是輻射流的流體法向速度向量。

根據在流固耦合濕表面上，結構和流體的法向速度相等的條件，得到:

式中:ui為入射流的流體法向速度向量。

將式(2)代入式(1)，式(4)代入式(3)，得到了下面的流固耦合方程組:

將式(8)左右兩邊同乘以GAfM－1f，得到:

這樣式(7)和式(9)就可以聯立進行數值耦合求解。

1.2 氣泡載荷

沖擊波之后，大約47%的爆炸能量仍然存留在氣泡中。假設流體為無旋，不可壓縮的，氣泡中心在自由水面下深度d處。由于本文考慮的是中遠場爆炸，假設船體對氣泡沒有影響，并且氣泡在運動過程中保持球形。為考慮氣泡的遷移效應，自由面效應和氣泡阻力的影響，采用Vernon[27]的無量綱方程組:

式中:x是無量綱氣泡半徑，ζ0為無量綱初始壓頭。ζ為無量綱壓頭，τ為無量綱時間;δ為無量綱深度，k是無量綱能量參數(對于 TNT炸藥，k≈0.074 3(z0)1/4，z0為初始壓頭)，γ 為絕熱氣體參數，取為 1.25[11]，Cd為阻尼系數，取為2.5。α為氣泡遷移控制系數(α只取值0或1;α=0表示不考慮氣泡的遷移效應)。在本文中取α=1。β為自由面效應控制系數，同樣只取值為0或1;β=0表示不考慮自由面效應。在本文中取β=1(即考慮自由面效應)。

接著定義L為長度尺度因子，T為時間尺度因子。

式中:E0為爆炸的總能量，ρ為流體的密度，g為重力加速度，W為炸藥量。可以得到以下無量綱參數:

其中:z為壓力水頭。

只要給定初始條件，式(10a)～(10d)可以由四階龍格庫塔法進行求解。選定初始條件為x=x0，ζ=ζ0，σ=0，λ =0。而其中 x0可以由下面能量守恒方程[25]得到:

方程(10a)～(10d)的解即為x(t)和ζ(t)。速度勢函數可以表示為:

式中:r為所取源點到氣泡中心的徑向距離，r1為該源點對應的偶極到氣泡中心的徑向距離。θ為r方向與垂直方向的夾角，θ1為r1方向與垂直方向的夾角。e1為源強系數，e2為偶極強度系數，定義如下:

對于勢流，速度為速度勢的負梯度，即為:

式中:X和Y分別為以氣泡中心作為原點的直角坐標系的橫坐標和縱坐標。

流體加速度可以表示為:

其中:ν為氣泡的垂向平均速度。

式(20)中的兩項可以分別表達為下面兩式:

2 求解方法

本文的程序中，首先利用有限元方法提取結構質量矩陣和剛度矩陣，然后利用一階DAA方法來求解流體方程。一階DAA方法實質上是在高頻頻域段和低頻頻域段分別采用平面波近似理論和虛質量近似理論進行逼近近似，中頻段采用線性的過度。這樣就可以在從高頻頻域到低頻頻域的區域都具有較高精度。文獻[6]通過對浸沒在水中的圓形球殼的數值計算結果與精確解的比較，驗證了一階DAA方法在求解一般問題上的精度可以滿足要求。而程序中的三維流體質量矩陣利用邊界積分方法進行求解[28]。在進行耦合計算時，采用分部交錯迭代(Staggered solution)的方法，在每一個時間步內，采用預報－代入－判斷－同步的方式進行迭代。結構的計算采用Wilson－θ法，考慮到計算的穩定性，流體的計算采用隱式的單步求解方法。計算程序流體圖如圖1所示。

圖1 計算程序流程圖Fig.1 Flow chart of the calculation procedure

3 算例分析

取一艘船模作為算例來分析其在氣泡載荷作用下的動態響應。示意圖如圖2所示。此船模總長4.5 m，型寬0.6 m，型深0.45 m，吃水T=0.265 m，排水量390 kg。船模共設7道橫艙壁、兩層平臺，有中縱艙壁，平臺2與船底外底板組成雙層底。甲板、平臺1、平臺2的板厚4 mm，其余結構的板厚3 mm。采用普通鋼建造，密度 ρ=7 850 kg/m3，彈性模量 E=210 GPa，泊松比 μ=0.3，屈服極限約為235 MPa。

圖2 船模的結構示意圖Fig.2 Schematic diagram of the ship model

船模的有限元模型如圖3所示。船模共有2 373個節點，2 648個單元，其中包括661個濕表面單元。對于非接觸水下爆炸，破壞力最大的情況是炸藥在艦船正下方發生爆炸。所以本文中，假設藥包在船中正下方爆炸，如圖3所示，爆炸的位置可以由爆炸深度h來表示。

圖3 船模的有限元模型Fig.3 Finite element model of the ship-model

取6 kg TNT炸藥在自由表面下10 m處爆炸。為忽略沖擊波效應的影響，取爆炸后0.05 s作為初始時間。氣泡半徑，氣泡壓力和氣泡中心深度的時間歷程曲線如圖4所示。從圖中可以看到，氣泡半徑先增大到最大值，然后在t=0.26 s時減小到最小半徑。在這個過程中，隨著氣泡半徑的減小，氣氣泡壓力隨著氣泡半徑的減小而增加。當氣泡半徑減小到最小值之后，氣泡開始回彈，因為此時氣泡內部壓力已經變得非常大。隨后氣泡半徑隨著時間而增大。在回彈過程中，氣泡壓力迅速衰減到零以下。當氣泡半徑在最大值附近時，氣泡中心位置上升的很緩慢，而當氣泡半徑要達到最小值時，氣泡開始快速地向上遷移。

圖4 氣泡半徑，深度和壓力的時間歷程曲線Fig.4 Time histories of bubble radius，depth and pressure

圖5 給出了一系列不同時刻的船模的垂線位移響應。其中云圖給出的是全船的響應分布情況，云圖下方的曲線表示船底龍骨的位移響應。從圖5中可以清楚地觀察到船模在氣泡在作用下的位移響應主要為船體的總體響應。船模按梁的一階垂線振型進行總體鞭狀運動。這主要因為氣泡載荷的頻率較低，與船體的低階垂向固有頻率相接近，因此激起了船體的低階振型，使船體發生鞭狀運動。船中處位移最大，圖中可以看到最大位移可達到2.4 cm。這樣的鞭狀振動使船中反復遭受作用力，嚴重時會對船體造成較大的總體損傷。

為進一步研究船模在不同位置處的響應特性，在船模的中線面上取幾個典型的位置作為測點。其中包括船底龍骨上的測點(B1，B2，B3)和主甲板上的測點(D1，D2，D3)。測點的具體位置如6所示。

圖5 不同時刻船模的垂向位移響應Fig.5 Vertical displacement responses of the ship model at different times

圖6 不同測點的分布位置Fig.6 The locations of different test points

圖7 給出了不同測點的加速度，速度和位移響應的時程曲線。從圖中可以觀察到所有測點的加速度，速度和位移響應均開始于大約0.25 s時，這是氣泡載荷開始迅速增長至峰值的時刻。然后響應的幅值開始緩慢平滑的衰減。由于炸藥位于船中的正下方，所以船中的測點B2和D2的加速度，速度和位移的峰值要比其他測點大很多。如B2的最大加速度值為1 869.6 m/s2，分別約為 B1(726.8 m/s2)和 B3(1 474.8.8 m/s2)的2.57倍和1.27倍。速度和位移響應也遵循同樣的規律。在同一縱向位置，主甲板和船底龍骨處的響應相差不大。如B1和D1，B2和D2，B3和D3的加速度，速度和位移的峰值都比較接近。

圖7 船模不同測點的加速度、速度和位移響應Fig.7 Acceleration，velocity and displacement responses of the ship model at different test points

除了船體的總體響應，船模上的一些區域也出現了明顯的局部響應，如圖8和圖9所示。圖8給出了在t=0.260 s時刻的垂向位移響應。可以觀察到在船中處的舷側外板上有明顯的局部大位移響應區域。在此區域取兩個測點T1和T2，如圖所示，此時T1的垂向位移為3.19 cm，T2的垂向位移為2.73 cm。兩者相差17%。結合圖5(a)～(c)中的位移云圖，可以看到沿著船長方向，在舷側外板上還有幾個位置均出現了局部大位移響應區域，分布在船中兩側，比船中處小一些。

圖8 船模的局部垂向位移響應Fig.8 Local vertical displacement response of the ship model

圖9 給出的是t=0.271 s時刻船模的垂向加速度響應。可以觀察到在主甲板的艙口處，出現了明顯的局部加速度集中區域。同樣取測點T3和T4，如圖所示，此時T3的垂向加速度為3 166.12 m/s2，T2的垂向加速度為1 520.66 m/s2。T3為T4的2.08倍。這是因為船舶甲板艙口角隅，由于形狀不連續，在船體發生鞭狀運動受到較大的面內載荷時，使局部的應力梯度升高，產生應力集中，嚴重時可能會造成塑性變形與屈服。因此在結構設計時，為了降低這種應力集中程度，應該采取加厚板、復板或形狀優化設計等使艙口角隅的最大加速度極小化。

圖9 船模的局部垂向加速度響應Fig.9 Local vertical acceleration response of the ship model

通常認為，水下爆炸主要引起船體的垂直方向的響應。其實在某些位置，船體的橫向響應也較為明顯。圖10給出的是t=0.258 s時刻船模的橫向位移響應。可以觀察到沿著船長方向，在舷側外板上有多處較大的局部大位移區域，尤其在艏尖艙壁位置尤為明顯，在艏尖艙壁處出現了面積非常大的局部大位移響應區域，而且左右兩側的位移是相反方向的。在左右兩側分別取兩個對稱的測點T5和T6，它們的加速度時間歷程曲線如圖11所示。可見兩個測點一直在做相位相反的振動，即按殼體的呼吸模態振動。而兩個測點T5和T6的橫向加速度峰值分別為 2 899.9 m/s2和2 784.7 m/s2。它們已經和上面分析的垂向加速度在同一量級，嚴重時會造成船體結構的局部破壞，必須引起重視。

圖10 船模的局部橫向位移響應Fig.10 Local transverse displacement response of the ship model

圖11 測點T5和T6處的橫向加速度響應比較Fig.11 Comparison of transverse direction accelerations at T5 and T6

4 結論

本文通過建立一套有限元方法與DAA方法相結合的計算程序，研究了水下爆炸氣泡作用下水面艦船的動態響應。主要結論為:

(1)水下爆炸氣泡脈動載荷作用下，船體結構的響應主要以總體響應為主。船體發生了顯著的鞭狀運動。且主要表現為低階模態響應，即一階彈性振型。

(2)氣泡作用下，在船體的某些位置，局部垂向響應也非常顯著。如船中處的舷側外板位置有明顯的局部大位移響應;在主甲板的艙口處，出現了明顯的局部加速度集中區域。

(3)在船模的一些典型區域，如舷側外板，局部橫向響應也比較明顯，主要表現為兩側外板按殼體的呼吸模態振動。尤其艏尖艙壁區域，橫向響應很大，加速度已經達到和垂向響應相同的量級。

[1]Mindlin R D，Bleich H H.Response of an elastic cylindrical shell to a transverse step shock wave[J].Journal of Applied Mechanics，1953，20:189 －195.

[2]Haywood J H.Response of an elastic cylindrical shell to a pressure pulse[J].The Quarterly Journal of Mechanics and Applied Mathematics，1958，11:129 －141.

[3]Chertock G.The transient flexural vibrations of ship-like structures exposed to underwater explosions[J].Journal of the Acoustical Society of America，1970，48(l):170 －180.

[4]Geers T L.Residual potential and approximate methods for three-dimensional fluid-structure interaction problems[J].The Journal of the Acoustical Society of America，1971，49:1505－1510.

[5]Geers T L.Doubly asymptotic approximations for transient motions of submerged structures[J].The Journal of the Acoustical Society of America，1978，64:1500－1508.

[6]Geers T L，Felippa C A.Doubly asymptotic approximations for vibration analysis of submerged structures[J].The Journal of the Acoustical Society of America，1983，73:1152 －1159.

[7]Ergin A.The response behaviour of a submerged cylindrical shell using the doubly asymptotic approximation method(DAA)[J].Computers and Structures，1997，62(6):1025－1034.

[8]Liang C C，Hsu C Y，Lai W H.A study of transient responses of a submerged spherical shell under shock waves[J].Ocean Engineering，2000，28:71 －94.

[9]劉建湖.船舶非接觸水下爆炸動力學的理論和應用[D].無錫:中國船舶科學研究中心，2002.LIU Jian-hu.Theory and its applications of ship dynamic responses to non-contact underwater explosions[D].Wuxi:China Ship Scientific Research Center，2002.

[10]姚熊亮，孫士麗，陳 玉，等.非線性雙漸進法應用于水中結構瞬態運動的研究[J].振動與沖擊，2010，10(29):9－15.YAO Xiong-liang，SUN Shi-li，CHEN Yu，et al.Transient motions of submerged structures with nonlinear fluid-structure interaction method [J].Journal of Vibration and Shock，2010，10(29):9 －15.

[11]Cole R H.Underwater explosion[M].New Jersey:Princeton University Press，1948.

[12]Keil A H.The response of ships to underwater explosions[J].Transactions of the Society of Naval Architects and Marine Engineers，1961，69:366 －410.

[13]Rayleigh L.On the pressure developed in a liquid during the collapse of a spherical void [J].Philosophical Magazine，1917，34:94 －98.

[14]Benjamin T B，Ellis A T.Cavitation:The collapse of cavitation bubbles and the pressures thereby produced against solid boundaries[J].Philosophical Transactions of the Royal Society of London，1966，260:221－240.

[15]Plesset M S，Chapman R B.Collapse of an initially spherical vapour cavity in the neighbourhood of a solid boundary [J].Journal of Fluid Mechanics，1971，47:283 －290.

[16]Gibson D C，Blake J R.The growth and collapse of bubbles near deformable surfaces[J].Applied Scientific Research，1982，38:215－224.

[17]Best J P，Kucera A A.Numerical investigation of nonsphericalrebounding bubbles [J]. JournalofFluid Mechanics，1992，245:137 －154.

[18]Zhang S，Duncan J H，Chahine G L.The final stage of the collapse of a cavitation bubble near a rigid wall[J].Journal of Fluid Mechanics，1993，257:147－181.

[19]Klaseboer E，Khoo B C，Hung K C.Dynamics of an oscillating bubble near a floating structure[J].Journal of Fluids and Structure，2005，21:395－412.

[20]Hicks A N.Explosion induced hull whipping.In:Advances in marine structures(Smith，C.S.，Clarke，J.D.，eds.)[C].Int.Conf.on Advances in Marine Structures.London:Elsevier 1986.

[21]Smiljanic B，Bobanac N，Senjanovic I.Bending moment of ship hull girder caused by pulsating bubble of underwater explosion[J]. In:Faltinsen， O.(Ed.)， Proceedings Hydroelasticity in Marine Technology. A.A.Balkema，Trondheim，1994.pp.149－156.

[22]Zong Z.Dynamic plastic response of a submerged free-free beam to an underwater gas bubble[J].Acta Mechanics，2003，161:179－214.

[23]Zong Z.A hydroplastic analysis of a free-free beam floating on water subjected to an underwater bubble[J].Journal of Fluids and Structures，2005，20:359 －372.

[24]Zhang N，Zong Z.The effect of rigid-body motions on the whipping response of a ship hull subjected to an underwater bubble[J].Journal of Fluids and Structures，2011，27:1326－1336.

[25]姚熊亮，陳建平，任慧龍.水下爆炸氣泡脈動壓力下艦船動態響應分析[J].哈爾濱工程大學學報，2001，42(2):48－55.YAO Xiong-liang， CHEN Jian-ping， REN Hui-long. The analysis of dynamic response of ship hull subjected to gas bubble impulsive pressure of underwater explosions[J].Journal of Harbin Engineering University，2001，42(2):48－55.

[26]李玉節，張孝慈，吳有生，等.水下爆炸氣泡激起的船體鞭狀運動[J].中國造船，2001，42(3):1－7.LIYu-jie， ZHANG Xiao-ci， WU You-sheng， etal.Whipping response of ship hull induced by underwater explosion bubble[J].Shipbuilding of China，2001，42(3):1－7.

[27]Vermon T A. Whippingresponse ofship hullsfrom underwaterexplosion bubble loading[J]. Technical Memorandum 86/255，Defence Research Establishment Atlantic，1986:1 －41.

[28]Deruntz J A，Geers T L.Added mass computation by the boundary integral method[J].International Journal for Numerical Methods in Engineering，1978，12:531－550.