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

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗室 船舶工程學(xué)院，大連 116024)

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

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

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

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

1 理論背景

1.1 控制方程

一個彈性結(jié)構(gòu)受到外部激勵的運(yùn)動方程可以表示為:

其中:Ms，Cs，Ks分別為N×N階的質(zhì)量矩陣，阻尼矩陣和剛度矩陣。N×1階向量分別為結(jié)構(gòu)的位移，速度和加速度。列向量F代表外力。N為結(jié)構(gòu)的自由度總數(shù)量。

對于浸沒在流體的結(jié)構(gòu)所受的外部激勵可以表示為:

式中:pi和ps分別為濕表面上入射流和輻射流作用下產(chǎn)生的節(jié)點(diǎn)壓力向量，Af為濕表面單元的面積對角矩陣，G為結(jié)構(gòu)節(jié)點(diǎn)力和濕表面節(jié)點(diǎn)力的坐標(biāo)轉(zhuǎn)換矩陣。

當(dāng)一個水中結(jié)構(gòu)如水面艦船或潛艇浸沒在無限聲學(xué)流體介質(zhì)中，結(jié)構(gòu)濕表面的殼單元的控制方程可以由DAA方法來表示。對于流體中的彈性結(jié)構(gòu)，DAA方法表示的結(jié)構(gòu)表面流體的運(yùn)動即為各正交流體邊界模態(tài)的線性組合。一階雙漸近DAA方程的矩陣形式表達(dá)為:

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

根據(jù)在流固耦合濕表面上，結(jié)構(gòu)和流體的法向速度相等的條件，得到:

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

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

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

這樣式(7)和式(9)就可以聯(lián)立進(jìn)行數(shù)值耦合求解。

1.2 氣泡載荷

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

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

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

式中:E0為爆炸的總能量，ρ為流體的密度，g為重力加速度，W為炸藥量?？梢缘玫揭韵聼o量綱參數(shù):

其中:z為壓力水頭。

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

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

式中:r為所取源點(diǎn)到氣泡中心的徑向距離，r1為該源點(diǎn)對應(yīng)的偶極到氣泡中心的徑向距離。θ為r方向與垂直方向的夾角，θ1為r1方向與垂直方向的夾角。e1為源強(qiáng)系數(shù)，e2為偶極強(qiáng)度系數(shù)，定義如下:

對于勢流，速度為速度勢的負(fù)梯度，即為:

式中:X和Y分別為以氣泡中心作為原點(diǎn)的直角坐標(biāo)系的橫坐標(biāo)和縱坐標(biāo)。

流體加速度可以表示為:

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

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

2 求解方法

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

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

3 算例分析

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

圖2 船模的結(jié)構(gòu)示意圖Fig.2 Schematic diagram of the ship model

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

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

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

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

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

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

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

圖6 不同測點(diǎn)的分布位置Fig.6 The locations of different test points

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

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

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

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

圖9 給出的是t=0.271 s時刻船模的垂向加速度響應(yīng)?？梢杂^察到在主甲板的艙口處，出現(xiàn)了明顯的局部加速度集中區(qū)域。同樣取測點(diǎn)T3和T4，如圖所示，此時T3的垂向加速度為3 166.12 m/s2，T2的垂向加速度為1 520.66 m/s2。T3為T4的2.08倍。這是因為船舶甲板艙口角隅，由于形狀不連續(xù)，在船體發(fā)生鞭狀運(yùn)動受到較大的面內(nèi)載荷時，使局部的應(yīng)力梯度升高，產(chǎn)生應(yīng)力集中，嚴(yán)重時可能會造成塑性變形與屈服。因此在結(jié)構(gòu)設(shè)計時，為了降低這種應(yīng)力集中程度，應(yīng)該采取加厚板、復(fù)板或形狀優(yōu)化設(shè)計等使艙口角隅的最大加速度極小化。

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

通常認(rèn)為，水下爆炸主要引起船體的垂直方向的響應(yīng)。其實(shí)在某些位置，船體的橫向響應(yīng)也較為明顯。圖10給出的是t=0.258 s時刻船模的橫向位移響應(yīng)?？梢杂^察到沿著船長方向，在舷側(cè)外板上有多處較大的局部大位移區(qū)域，尤其在艏尖艙壁位置尤為明顯，在艏尖艙壁處出現(xiàn)了面積非常大的局部大位移響應(yīng)區(qū)域，而且左右兩側(cè)的位移是相反方向的。在左右兩側(cè)分別取兩個對稱的測點(diǎn)T5和T6，它們的加速度時間歷程曲線如圖11所示?？梢妰蓚€測點(diǎn)一直在做相位相反的振動，即按殼體的呼吸模態(tài)振動。而兩個測點(diǎn)T5和T6的橫向加速度峰值分別為 2 899.9 m/s2和2 784.7 m/s2。它們已經(jīng)和上面分析的垂向加速度在同一量級，嚴(yán)重時會造成船體結(jié)構(gòu)的局部破壞，必須引起重視。

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

圖11 測點(diǎn)T5和T6處的橫向加速度響應(yīng)比較Fig.11 Comparison of transverse direction accelerations at T5 and T6

4 結(jié)論

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

(1)水下爆炸氣泡脈動載荷作用下，船體結(jié)構(gòu)的響應(yīng)主要以總體響應(yīng)為主。船體發(fā)生了顯著的鞭狀運(yùn)動。且主要表現(xiàn)為低階模態(tài)響應(yīng)，即一階彈性振型。

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

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

[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，F(xiàn)elippa 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]劉建湖.船舶非接觸水下爆炸動力學(xué)的理論和應(yīng)用[D].無錫:中國船舶科學(xué)研究中心，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ìn)法應(yīng)用于水中結(jié)構(gòu)瞬態(tài)運(yùn)動的研究[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]姚熊亮，陳建平，任慧龍.水下爆炸氣泡脈動壓力下艦船動態(tài)響應(yīng)分析[J].哈爾濱工程大學(xué)學(xué)報，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]李玉節(jié)，張孝慈，吳有生，等.水下爆炸氣泡激起的船體鞭狀運(yùn)動[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.