水下爆炸氣泡作用下船模的動(dòng)態(tài)響應(yīng)

水下爆炸氣泡作用下船模的動(dòng)態(tài)響應(yīng)

這是一個(gè)非常復(fù)雜的問題，包括模擬水下爆炸負(fù)荷、爆炸氣泡的動(dòng)態(tài)性能、固結(jié)界面分析、水下結(jié)構(gòu)的非線性動(dòng)態(tài)響應(yīng)等。目前,對(duì)于這類三維結(jié)構(gòu)動(dòng)態(tài)問題的主要研究方法是邊界元/有限元結(jié)合的方法。主要采用流固耦合方程來(lái)表征耦合界面上的載荷,而不需要對(duì)水域進(jìn)行建模,因此減少了大量的計(jì)算量。這種方法主要經(jīng)歷了平面波理論,曲面波理論,虛質(zhì)量理論到雙漸近(DAA)理論的發(fā)展過程。特別是由Geers所提出的DAA理論,綜合并改進(jìn)了前面幾種方法,在高、中、低頻段都有較好的精度。DAA方法經(jīng)過眾多學(xué)者們的改進(jìn)與發(fā)展,目前是一種解決流固耦合問題的非常有效的計(jì)算方法。水下爆炸中對(duì)艦船主要產(chǎn)生兩種載荷:一是瞬態(tài)的沖擊波,其壓力峰值很高但持續(xù)時(shí)間很短,通常造成水中結(jié)構(gòu)的局部破壞;二是氣泡脈動(dòng)載荷,氣泡中大約含有47%的爆炸能量,具有非常強(qiáng)的破壞力。這一階段的明顯特征是達(dá)到的壓力峰值較沖擊波的壓力峰值低,而且持續(xù)時(shí)間很長(zhǎng),約為幾百毫秒,會(huì)造成水中結(jié)構(gòu)的總體響應(yīng)和局部響應(yīng)。對(duì)于氣泡的動(dòng)力學(xué),眾多學(xué)者已經(jīng)進(jìn)行了大量的研究工作[13,14,15,16,17,18,19]。對(duì)于氣泡載荷作用下的艦船的動(dòng)態(tài)響應(yīng),目前的主要研究方法都是將艦船簡(jiǎn)化為船體梁,研究船體梁的總體鞭狀響應(yīng)[20,21,22,23,24,25,26]。然而,對(duì)于氣泡作用下三維全船結(jié)構(gòu)的動(dòng)態(tài)響應(yīng)的細(xì)節(jié)分析,開展的研究還很少。本文致力于建立一套有限元方法與DAA方法相結(jié)合的計(jì)算程序,用以研究水下爆炸氣泡作用下水面艦船的動(dòng)態(tài)響應(yīng)。首先,建立一個(gè)考慮遷移效應(yīng),自由面效應(yīng)和氣泡阻力的氣泡模型和一個(gè)水面艦船結(jié)構(gòu)模型,并闡述了流固耦合與結(jié)構(gòu)響應(yīng)相關(guān)的理論分析;然后以該船模作為算例,研究了氣泡作用下船體模型的總體響應(yīng)和局部響應(yīng)。比較了不同位置的加速度,速度和位移時(shí)程曲線。最后詳細(xì)討論了氣泡作用下船體模型動(dòng)態(tài)響應(yīng)的特征與機(jī)理。1理論背景1.1流固耦合方程的建立一個(gè)彈性結(jié)構(gòu)受到外部激勵(lì)的運(yùn)動(dòng)方程可以表示為:Μs??x+Cs˙x+Κsx=F(1)Msx??+Csx˙+Ksx=F(1)其中:Ms,Cs,Ks分別為N×N階的質(zhì)量矩陣,阻尼矩陣和剛度矩陣。N×1階向量x,˙x,??xx,x˙,x??分別為結(jié)構(gòu)的位移,速度和加速度。列向量F代表外力。N為結(jié)構(gòu)的自由度總數(shù)量。對(duì)于浸沒在流體的結(jié)構(gòu)所受的外部激勵(lì)可以表示為:F=-GAf(pi+ps)(2)式中:pi和ps分別為濕表面上入射流和輻射流作用下產(chǎn)生的節(jié)點(diǎn)壓力向量,Af為濕表面單元的面積對(duì)角矩陣,G為結(jié)構(gòu)節(jié)點(diǎn)力和濕表面節(jié)點(diǎn)力的坐標(biāo)轉(zhuǎn)換矩陣。當(dāng)一個(gè)水中結(jié)構(gòu)如水面艦船或潛艇浸沒在無(wú)限聲學(xué)流體介質(zhì)中,結(jié)構(gòu)濕表面的殼單元的控制方程可以由DAA方法來(lái)表示。對(duì)于流體中的彈性結(jié)構(gòu),DAA方法表示的結(jié)構(gòu)表面流體的運(yùn)動(dòng)即為各正交流體邊界模態(tài)的線性組合。一階雙漸近DAA方程的矩陣形式表達(dá)為:Μf˙ps+ρcAfps=ρc˙us(3)Mfp˙s+ρcAfps=ρcu˙s(3)式中:Mf為N×N階的流體質(zhì)量矩陣;ρ和c分別為流體密度和流體中聲速;us是輻射流的流體法向速度向量。根據(jù)在流固耦合濕表面上,結(jié)構(gòu)和流體的法向速度相等的條件,得到:GΤ˙x=ui+us(4)GTx˙=ui+us(4)式中:ui為入射流的流體法向速度向量。將式(2)代入式(1),式(4)代入式(3),得到了下面的流固耦合方程組:Μs??x+Cs˙x+Κsx=-GAf(pi+ps)(5)Msx??+Csx˙+Ksx=?GAf(pi+ps)(5)Μf˙ps+ρcAfps=ρc(GΤ??x-˙ui)(6)Mfp˙s+ρcAfps=ρc(GTx???u˙i)(6)為了求解的方便,將式(5)和式(6)改寫成對(duì)稱形式。令ps=GΤ˙y?ps=GTy˙?這樣式(5)和式(6)可以改寫成:Μs??x+Cs˙x+Κsx=-GAf(pi+GΤ˙y)(7)ΜfGΤ??y+ρcAfGΤ˙y=ρc(GΤ??x-˙ui)(8)將式(8)左右兩邊同乘以GAfM-1f,得到:GAfGΤ??y+ρcGAfΜ-1fAfGΤ˙y=ρcGAf(GΤ??x-˙ui)(9)這樣式(7)和式(9)就可以聯(lián)立進(jìn)行數(shù)值耦合求解。1.2階龍格庫(kù)塔法沖擊波之后,大約47%的爆炸能量仍然存留在氣泡中。假設(shè)流體為無(wú)旋,不可壓縮的,氣泡中心在自由水面下深度d處。由于本文考慮的是中遠(yuǎn)場(chǎng)爆炸,假設(shè)船體對(duì)氣泡沒有影響,并且氣泡在運(yùn)動(dòng)過程中保持球形。為考慮氣泡的遷移效應(yīng),自由面效應(yīng)和氣泡阻力的影響,采用Vernon的無(wú)量綱方程組:σ=˙x?(10a)λ=˙ζ?(10b)˙σ=-3δ(2δ-βx)[σ2x(1-2βx3δ)-λ26x+ζxζ0-(γ-1)kx3γ+1+βx4δ2(Cdλ24x+σλ3-xζ0)](10c)˙λ=-3α[1ζ0+σλx-Cdλ24x+βx4δ2(3σ2+x˙σ)](10d)式中:x是無(wú)量綱氣泡半徑,ζ0為無(wú)量綱初始?jí)侯^。ζ為無(wú)量綱壓頭,τ為無(wú)量綱時(shí)間;δ為無(wú)量綱深度,k是無(wú)量綱能量參數(shù)(對(duì)于TNT炸藥,k≈0.0743(z0)1/4,z0為初始?jí)侯^),γ為絕熱氣體參數(shù),取為1.25,Cd為阻尼系數(shù),取為2.5。α為氣泡遷移控制系數(shù)(α只取值0或1;α=0表示不考慮氣泡的遷移效應(yīng))。在本文中取α=1。β為自由面效應(yīng)控制系數(shù),同樣只取值為0或1;β=0表示不考慮自由面效應(yīng)。在本文中取β=1(即考慮自由面效應(yīng))。接著定義L為長(zhǎng)度尺度因子,T為時(shí)間尺度因子。L=[3E04πρgz0]13?Τ=1.428W13z560(11)式中:E0為爆炸的總能量,ρ為流體的密度,g為重力加速度,W為炸藥量?？梢缘玫揭韵聼o(wú)量綱參數(shù):x=aL?ζ=zL?ζ0=z0L?δ=dL?τ=tΤ(12)其中:z為壓力水頭。只要給定初始條件,式(10a)～(10d)可以由四階龍格庫(kù)塔法進(jìn)行求解。選定初始條件為x=x0,ζ=ζ0,σ=0,λ=0。而其中x0可以由下面能量守恒方程得到:˙x2x3+kx3(γ+1)+x3-1=0(13)方程(10a)～(10d)的解即為x(t)和ζ(t)。速度勢(shì)函數(shù)可以表示為:φ=e1r+e2r2cosθ-e1r1+e2r21cosθ1(14)式中:r為所取源點(diǎn)到氣泡中心的徑向距離,r1為該源點(diǎn)對(duì)應(yīng)的偶極到氣泡中心的徑向距離。θ為r方向與垂直方向的夾角,θ1為r1方向與垂直方向的夾角。e1為源強(qiáng)系數(shù),e2為偶極強(qiáng)度系數(shù),定義如下:e1=L3x2Τσ(15)e2=L42Τ(x3λ-x5σ4δ2)(16)˙e1=L3Τ2(2xσ2+x2˙σ)(17)˙e2=3L4σx22Τ2[λ-x2σ4δ2]+x3L42Τ2[˙λ-(2xσ+x2˙σ4δ2)](18)對(duì)于勢(shì)流,速度為速度勢(shì)的負(fù)梯度,即為:uΖ=-?φ=1(X2+Y2+Ζ2)[e1Ζ(X2+Y2+Ζ2)12-e2(X2+Y2+Ζ2)12(1-3Ζ2(X2+Y2+Ζ2))](19)式中:X和Y分別為以氣泡中心作為原點(diǎn)的直角坐標(biāo)系的橫坐標(biāo)和縱坐標(biāo)。流體加速度可以表示為:˙uΖ=?uΖ?t-v?uΖ?Ζ(20)其中:ν為氣泡的垂向平均速度。式(20)中的兩項(xiàng)可以分別表達(dá)為下面兩式:?uΖ?t=1(X2+Y2+Ζ2)32[˙e1Ζ-˙e2(1-3Ζ2(X2+Y2+Ζ2))](21)?uΖ?Ζ=e1(X2+Y2+Ζ2)32(1-3Ζ2X2+Y2+Ζ2)+3Ζe2(X2+Y2+Ζ2)32(3-5Ζ2X2+Y2+Ζ2)(22)這樣,可以計(jì)算每一時(shí)間步的e1,e2,˙e1,˙e2并把它們代入式(21)和式(22)中,然后由式(20),就可以求解得到流體加速度˙uγ。同樣,任意點(diǎn)處的氣泡壓力p可以由伯努利方程p=-ρ?φ/?t-ρgw來(lái)計(jì)算,其中w為測(cè)點(diǎn)到自由面的垂直距離。2數(shù)值優(yōu)化算法本文的程序中,首先利用有限元方法提取結(jié)構(gòu)質(zhì)量矩陣和剛度矩陣,然后利用一階DAA方法來(lái)求解流體方程。一階DAA方法實(shí)質(zhì)上是在高頻頻域段和低頻頻域段分別采用平面波近似理論和虛質(zhì)量近似理論進(jìn)行逼近近似,中頻段采用線性的過度。這樣就可以在從高頻頻域到低頻頻域的區(qū)域都具有較高精度。文獻(xiàn)通過對(duì)浸沒在水中的圓形球殼的數(shù)值計(jì)算結(jié)果與精確解的比較,驗(yàn)證了一階DAA方法在求解一般問題上的精度可以滿足要求。而程序中的三維流體質(zhì)量矩陣?yán)眠吔绶e分方法進(jìn)行求解。在進(jìn)行耦合計(jì)算時(shí),采用分部交錯(cuò)迭代(Staggeredsolution)的方法,在每一個(gè)時(shí)間步內(nèi),采用預(yù)報(bào)-代入-判斷-同步的方式進(jìn)行迭代。結(jié)構(gòu)的計(jì)算采用Wilson-θ法,考慮到計(jì)算的穩(wěn)定性,流體的計(jì)算采用隱式的單步求解方法。計(jì)算程序流體圖如圖1所示。3船模在船底內(nèi)部的垂線位移響應(yīng)取一艘船模作為算例來(lái)分析其在氣泡載荷作用下的動(dòng)態(tài)響應(yīng)。示意圖如圖2所示。此船模總長(zhǎng)4.5m,型寬0.6m,型深0.45m,吃水T=0.265m,排水量390kg。船模共設(shè)7道橫艙壁、兩層平臺(tái),有中縱艙壁,平臺(tái)2與船底外底板組成雙層底。甲板、平臺(tái)1、平臺(tái)2的板厚4mm,其余結(jié)構(gòu)的板厚3mm。采用普通鋼建造,密度ρ=7850kg/m3,彈性模量E=210GPa,泊松比μ=0.3,屈服極限約為235MPa。船模的有限元模型如圖3所示。船模共有2373個(gè)節(jié)點(diǎn),2648個(gè)單元,其中包括661個(gè)濕表面單元。對(duì)于非接觸水下爆炸,破壞力最大的情況是炸藥在艦船正下方發(fā)生爆炸。所以本文中,假設(shè)藥包在船中正下方爆炸,如圖3所示,爆炸的位置可以由爆炸深度h來(lái)表示。取6kgTNT炸藥在自由表面下10m處爆炸。為忽略沖擊波效應(yīng)的影響,取爆炸后0.05s作為初始時(shí)間。氣泡半徑,氣泡壓力和氣泡中心深度的時(shí)間歷程曲線如圖4所示。從圖中可以看到,氣泡半徑先增大到最大值,然后在t=0.26s時(shí)減小到最小半徑。在這個(gè)過程中,隨著氣泡半徑的減小,氣氣泡壓力隨著氣泡半徑的減小而增加。當(dāng)氣泡半徑減小到最小值之后,氣泡開始回彈,因?yàn)榇藭r(shí)氣泡內(nèi)部壓力已經(jīng)變得非常大。隨后氣泡半徑隨著時(shí)間而增大。在回彈過程中,氣泡壓力迅速衰減到零以下。當(dāng)氣泡半徑在最大值附近時(shí),氣泡中心位置上升的很緩慢,而當(dāng)氣泡半徑要達(dá)到最小值時(shí),氣泡開始快速地向上遷移。圖5給出了一系列不同時(shí)刻的船模的垂線位移響應(yīng)。其中云圖給出的是全船的響應(yīng)分布情況,云圖下方的曲線表示船底龍骨的位移響應(yīng)。從圖5中可以清楚地觀察到船模在氣泡在作用下的位移響應(yīng)主要為船體的總體響應(yīng)。船模按梁的一階垂線振型進(jìn)行總體鞭狀運(yùn)動(dòng)。這主要因?yàn)闅馀葺d荷的頻率較低,與船體的低階垂向固有頻率相接近,因此激起了船體的低階振型,使船體發(fā)生鞭狀運(yùn)動(dòng)。船中處位移最大,圖中可以看到最大位移可達(dá)到2.4cm。這樣的鞭狀振動(dòng)使船中反復(fù)遭受作用力,嚴(yán)重時(shí)會(huì)對(duì)船體造成較大的總體損傷。為進(jìn)一步研究船模在不同位置處的響應(yīng)特性,在船模的中線面上取幾個(gè)典型的位置作為測(cè)點(diǎn)。其中包括船底龍骨上的測(cè)點(diǎn)(B1,B2,B3)和主甲板上的測(cè)點(diǎn)(D1,D2,D3)。測(cè)點(diǎn)的具體位置如6所示。圖7給出了不同測(cè)點(diǎn)的加速度,速度和位移響應(yīng)的時(shí)程曲線。從圖中可以觀察到所有測(cè)點(diǎn)的加速度,速度和位移響應(yīng)均開始于大約0.25s時(shí),這是氣泡載荷開始迅速增長(zhǎng)至峰值的時(shí)刻。然后響應(yīng)的幅值開始緩慢平滑的衰減。由于炸藥位于船中的正下方,所以船中的測(cè)點(diǎn)B2和D2的加速度,速度和位移的峰值要比其他測(cè)點(diǎn)大很多。如B2的最大加速度值為1869.6m/s2,分別約為B1(726.8m/s2)和B3(1474.8.8m/s2)的2.57倍和1.27倍。速度和位移響應(yīng)也遵循同樣的規(guī)律。在同一縱向位置,主甲板和船底龍骨處的響應(yīng)相差不大。如B1和D1,B2和D2,B3和D3的加速度,速度和位移的峰值都比較接近。除了船體的總體響應(yīng),船模上的一些區(qū)域也出現(xiàn)了明顯的局部響應(yīng),如圖8和圖9所示。圖8給出了在t=0.260s時(shí)刻的垂向位移響應(yīng)?？梢杂^察到在船中處的舷側(cè)外板上有明顯的局部大位移響應(yīng)區(qū)域。在此區(qū)域取兩個(gè)測(cè)點(diǎn)T1和T2,如圖所示,此時(shí)T1的垂向位移為3.19cm,T2的垂向位移為2.73cm。兩者相差17%。結(jié)合圖5(a)～(c)中的位移云圖,可以看到沿著船長(zhǎng)方向,在舷側(cè)外板上還有幾個(gè)位置均出現(xiàn)了局部大位移響應(yīng)區(qū)域,分布在船中兩側(cè),比船中處小一些。圖9給出的是t=0.271s時(shí)刻船模的垂向加速度響應(yīng)?？梢杂^察到在主甲板的艙口處,出現(xiàn)了明顯的局部加速度集中區(qū)域。同樣取測(cè)點(diǎn)T3和T4,如圖所示,此時(shí)T3的垂向加速度為3166.12m/s2,T2的垂向加速度為1520.66m/s2。T3為T4的2.08倍。這是因?yàn)榇凹装迮摽诮怯?由于形狀不連續(xù),在船體發(fā)生鞭狀運(yùn)動(dòng)受到較大的面內(nèi)載荷時(shí),使局部的應(yīng)力梯度升高,產(chǎn)生應(yīng)力集中,嚴(yán)重時(shí)可能會(huì)造成塑性變形與屈服。因此在結(jié)構(gòu)設(shè)計(jì)時(shí),為了降低這種應(yīng)力集中程度,應(yīng)該采取加厚板、復(fù)板或形狀優(yōu)化設(shè)計(jì)等使艙口角隅的最大加速度極小化。通常認(rèn)為,水下爆炸主要引起船體的垂直方向的響應(yīng)。其實(shí)在某些位置,船體的橫向響應(yīng)也較為明顯。圖10給出的是t=0.258s時(shí)刻船模的橫向位移響應(yīng)?？梢杂^察到沿著船長(zhǎng)方向,在舷側(cè)外板上有多處較大的局部大位移區(qū)域,尤其在艏尖艙壁位置尤為明顯,在艏尖艙壁處出現(xiàn)了面積非常大的局部大位移響應(yīng)區(qū)域,而且左右兩側(cè)的位移是相反方向的。在左右兩側(cè)分別取兩個(gè)對(duì)稱的測(cè)點(diǎn)T5和T6,它們的加速度時(shí)間歷程曲線如圖11所示?？梢妰蓚€(gè)測(cè)點(diǎn)一直在做相位相反的振動(dòng),即按殼體的呼吸模態(tài)振動(dòng)。而兩個(gè)測(cè)點(diǎn)T5和T6的橫向加速度峰值分別為2899.9m/s2和2784.7m/s2。它們已經(jīng)和上面分析的垂向加速度在同一量級(jí),嚴(yán)重時(shí)會(huì)造成船體結(jié)構(gòu)的局部破壞,必須引起重視。4船