波浪中液艙流體漂蕩耦合運(yùn)動(dòng)分析

當(dāng)液體船航行時(shí)，液艙隨著船的波浪運(yùn)動(dòng)而旋轉(zhuǎn)。同時(shí)，液艙的旋轉(zhuǎn)力用于船體，改變了船舶的運(yùn)動(dòng)模式。近年來，隨著海洋油氣資源的開發(fā)和應(yīng)用，lkg和lg等特殊類型溶液貨運(yùn)船的開發(fā)和應(yīng)用，液體船和船舶之間的作用越來越受到重視。嚴(yán)重的室內(nèi)液體搖晃可能會(huì)給船體帶來巨大的破壞。尤其是在裝載液體船的過程中，由于巨大的搖晃，巨大的擾動(dòng)壓力可能會(huì)嚴(yán)重影響船體的結(jié)構(gòu)。因此，在船舶載載設(shè)計(jì)中，液體船沉陷和船舶運(yùn)動(dòng)之間的耦合作用也是值得考慮的重要因素之一。許多國內(nèi)外科學(xué)家對船體運(yùn)動(dòng)與船沉運(yùn)動(dòng)的結(jié)合效果進(jìn)行了分析和研究。例如，kim、roninrougebacke、faltin和新man等。動(dòng)態(tài)船沉壓法（例如新man）和非線性時(shí)區(qū)法（例如kim）等。從以往的研究結(jié)果來看，線性船體運(yùn)動(dòng)的理論方法基本上可以滿足船舶運(yùn)動(dòng)和空腔運(yùn)動(dòng)之間的耦合效應(yīng)分析和研究的需要。本文對船舶在波浪上運(yùn)動(dòng)和液艙液體晃蕩,即對船體內(nèi)外流場問題均采用了勢流理論方法求解,其中波浪中船體水動(dòng)力和時(shí)延函數(shù)采用切片法和脈沖響應(yīng)函數(shù)法計(jì)算獲得,艙內(nèi)液體非線性晃蕩采用時(shí)域邊界元法計(jì)算,并最終建立了在波浪中船體與液艙流體晃蕩耦合的運(yùn)動(dòng)方程.文中就S175加載方形液艙在迎浪、橫浪等不同工況下液艙流體晃蕩及其與船體運(yùn)動(dòng)耦合分別進(jìn)行了計(jì)算模擬與驗(yàn)證研究.研究中發(fā)現(xiàn),當(dāng)不出現(xiàn)液面破碎等強(qiáng)非線性現(xiàn)象時(shí),非線性時(shí)域邊界元法能夠給出較好的液艙流體晃蕩波形和壓力;S175的液艙加載50%的液體時(shí),船體與液艙晃蕩耦合運(yùn)動(dòng)時(shí)歷結(jié)果能清晰地反映液艙晃蕩對船體運(yùn)動(dòng)的影響,運(yùn)動(dòng)RAO能反映出不同頻率液艙對船體運(yùn)動(dòng)的影響程度;載液S175船橫搖運(yùn)動(dòng)RAO能準(zhǔn)確給出船體有無加載液艙時(shí)共振頻率的偏移現(xiàn)象.1基于液體的船舶耦合運(yùn)動(dòng)的基本理論1.1船舶運(yùn)動(dòng)方程航行船舶的運(yùn)動(dòng)計(jì)算所采用的參考坐標(biāo)系為oxyz,如圖1所示,oxy平面與靜水面重合,oz軸垂直向上,船舶以定速U0航行,在波浪激勵(lì)下作六自由度運(yùn)動(dòng).直接求解滿足定解條件的速度勢不是一件很容易的事,基于切片理論的STF法可以不必求解三維有航速相應(yīng)的邊值問題,只要求解船舶各個(gè)橫剖面的二維零航速水動(dòng)力問題即可,STF法的詳細(xì)介紹可見文獻(xiàn).若將船體橫剖面的速度勢?沿船長方向積分即可獲得整船的附加質(zhì)量和阻尼系數(shù):μij(ω)+iωλij(ω)=ρ∫Ldx∫C(x)?jnidl.(1)μij(ω)+iωλij(ω)=ρ∫Ldx∫C(x)?jnidl.(1)那么船體五自由度運(yùn)動(dòng)方程可表示如下:Σj=26[(mij+μij(ω)ξ¨j+λij(ω)ξ˙j+Cijξj]=Fwi(ω),(i,j=2,3,?,6).(2)Σj=26[(mij+μij(ω)ξ¨j+λij(ω)ξ˙j+Cijξj]=Fiw(ω),(i,j=2,3,?,6).(2)式中:ω為頻率;ξj為船體j模態(tài)的運(yùn)動(dòng)位移;mij為包含了液艙內(nèi)液體質(zhì)量的船體廣義質(zhì)量矩陣;μij(ω)、λij(ω)為頻域附加質(zhì)量和阻尼系數(shù);Cij為船體回復(fù)力系數(shù);FiW(ω)為作用在船體上的廣義波浪力;j=2,3,…,6分別對應(yīng)的運(yùn)動(dòng)模態(tài)為橫蕩、垂蕩、橫搖、縱搖和艏搖.切片法不能計(jì)及船體縱蕩運(yùn)動(dòng).1.2laplace方程的基本思想采用三維勢流理論對液艙內(nèi)非定常流體晃蕩問題進(jìn)行求解,液艙流體晃蕩是采用時(shí)域邊界元法進(jìn)行模擬計(jì)算的.液艙內(nèi)流體運(yùn)動(dòng)計(jì)算采用的坐標(biāo)系如圖2所示,是固定在液艙內(nèi)的局部坐標(biāo)系,其坐標(biāo)原點(diǎn)固定在靜止時(shí)自由液面的幾何中心,x軸正方向與液艙長度方向平行指向右端,y軸正方向與液艙寬度方向平行指向前端,z軸正方向?yàn)樨Q直向上,在計(jì)算過程中該坐標(biāo)系始終與液艙長寬高保持平行,其原點(diǎn)保持不變.假定流體不可壓縮,無粘,流動(dòng)無旋,則流體運(yùn)動(dòng)速度勢滿足拉普拉斯方程,即控制方程滿足:?2??x2+?2??y2+?2??z2=0?2??x2+?2??y2+?2??z2=0,在流場內(nèi).(3)三維流體晃蕩運(yùn)動(dòng)的流場的自由表面上動(dòng)力學(xué)邊界條件為???t+12|??|2?{(Ω×r)??}?+dVdt?r+(Ω×V)?r?gζ=0.(4)???t+12|??|2-{(Ω×r)??}?+dVdt?r+(Ω×V)?r-gζ=0.(4)其運(yùn)動(dòng)學(xué)邊界條件:nx?ζ?t=???n,(5)nx?ζ?t=???n,(5)壁面條件:???n=0.(6)???n=0.(6)式中:?為艙內(nèi)流體運(yùn)動(dòng)速度勢,Ω為液艙角速度矢量,V為液艙線速度矢量,r為位置矢量,n為邊界外法向,ζ為自由面升高.本文采用簡單格林函數(shù)法來求解Laplace方程,取1/r為格林函數(shù),則流場中任意一點(diǎn)的速度勢可以表述如下:α?(x0,y0,z0;t)=∫∫S[???n1r????n(1r)]ds.(7)t)=∫∫S[???n1r-???n(1r)]ds.(7)式中:α為源點(diǎn)處的歐拉角,通過對空間的離散求得每個(gè)時(shí)間步中的速度勢,再通過對自由表面動(dòng)力學(xué)條件進(jìn)行中心差分獲得下一個(gè)時(shí)間步的自由面條件,從而使得液艙晃蕩模擬可以在時(shí)域中反復(fù)進(jìn)行下去.本文的空間離散采用了線性面元,積分計(jì)算中源點(diǎn)與場點(diǎn)重合時(shí)積分為假性奇點(diǎn),采用解析法進(jìn)行了計(jì)算,源點(diǎn)與場點(diǎn)不重合時(shí)采用了七點(diǎn)高斯積分公式計(jì)算.計(jì)算中根據(jù)每個(gè)時(shí)間步獲得的速度勢,通過拉格朗日積分計(jì)算獲得靜壓力與動(dòng)壓力.1.3船舶回復(fù)力和附加質(zhì)量在線性系統(tǒng)中,任意激勵(lì)可以寫為脈沖響應(yīng)函數(shù)和激勵(lì)的卷積積分形式:x(t)=∫0tF(t?τ)h(τ)dτ+x(t=0).(8)x(t)=∫0tF(t-τ)h(τ)dτ+x(t=0).(8)式中:x(t)為在輸入h(t)下的系統(tǒng)響應(yīng),F(t)為單位脈沖輸入下的脈沖響應(yīng)函數(shù).將以上概念推廣到船舶運(yùn)動(dòng)問題,時(shí)域船體運(yùn)動(dòng)方程有著如下形式:Σj=26[Mij+μij(∞)]x¨j+∫?∞tKij(τ)x˙j(t?τ)dτ+b41x˙4+b42x˙24+[Cij]xj=FWi(t).(9)Σj=26[Μij+μij(∞)]x¨j+∫-∞tΚij(τ)x˙j(t-τ)dτ+b41x˙4+b42x˙42+[Cij]xj=FiW(t).(9)式中:Mij、μij(∞)和Cij分別代表載液船體質(zhì)量,無窮大遭遇頻率的附加質(zhì)量以及船舶回復(fù)力系數(shù);Kij(τ)為時(shí)延函數(shù);FWiiW(t)為作用在船體上的外部時(shí)域波浪力;b41和b42為考慮粘性作用的非線性船舶橫搖阻尼系數(shù),可以用經(jīng)驗(yàn)公式、物理實(shí)驗(yàn)等方法獲得.時(shí)延函數(shù)與無窮大遭遇頻率附加質(zhì)量可由以下頻域轉(zhuǎn)時(shí)域方法獲得:Kij(t)=2π∫∞0λij(ω)cos(ωt)dω,(10)μij(∞)=μij(ω)+1ω∫t0Kij(τ)sin(ωτ)dτ.(11)Κij(t)=2π∫0∞λij(ω)cos(ωt)dω,(10)μij(∞)=μij(ω)+1ω∫0tΚij(τ)sin(ωτ)dτ.(11)式中:μij(∞)指頻率無窮大時(shí)船體的附加質(zhì)量.求解運(yùn)動(dòng)方程時(shí),橫蕩與首搖模態(tài)的運(yùn)動(dòng)可采用數(shù)值彈簧技術(shù)來抑制其慢漂現(xiàn)象.所加入的橫蕩和首搖兩種運(yùn)動(dòng)模態(tài)彈簧剛度表達(dá)如下:Cij=[Mij+μij(∞)](2πTi2),i=2,6Cij=[Μij+μij(∞)](2πΤi2),i=2,6.(12)式中,周期Ti遠(yuǎn)大于波浪激勵(lì)周期.2船舶全局運(yùn)動(dòng)時(shí)域耦合船舶與液艙晃蕩耦合運(yùn)動(dòng)方程可表示為Σj=26[Mij+μij(∞)]x¨j+∫?∞tKij(τ)x˙j(t?τ)dτ+b41x˙4+b42x˙24+[Cij]xj=FWi(t)+FSloshi(t).(13)Σj=26[Μij+μij(∞)]x¨j+∫-∞tΚij(τ)x˙j(t-τ)dτ+b41x˙4+b42x˙42+[Cij]xj=FiW(t)+FiSlosh(t).(13)在計(jì)算模擬過程中,船體當(dāng)前時(shí)刻的運(yùn)動(dòng)規(guī)律傳遞到液艙,然后液艙晃蕩模擬程序計(jì)算出該時(shí)刻液艙內(nèi)部流場的速度勢,再通過對艙壁上的壓力積分得到液艙的晃蕩力和力矩FSlosh,并添加至式(13)的右端,采用四階龍格庫塔法進(jìn)行耦合求解,此時(shí)獲得的整船運(yùn)動(dòng)位移又作為下一時(shí)刻液艙流體晃蕩的條件進(jìn)行計(jì)算.由于船體運(yùn)動(dòng)預(yù)報(bào)與液艙晃蕩是在不同坐標(biāo)系中進(jìn)行的,所以二者之間力與運(yùn)動(dòng)需要在這2個(gè)坐標(biāo)系中進(jìn)行轉(zhuǎn)換,轉(zhuǎn)換方法可參考文獻(xiàn).這樣反復(fù)循環(huán)便可以得到船體與液艙耦合運(yùn)動(dòng)的時(shí)歷,這就是時(shí)域下耦合液艙晃蕩的船舶全局運(yùn)動(dòng)的求解方法.具體的運(yùn)算步驟如下:1)讀取船型液艙參數(shù)及頻域水動(dòng)力數(shù)據(jù);2)根據(jù)來波工況計(jì)算波浪力時(shí)歷;3)計(jì)算時(shí)延函數(shù),確定計(jì)算時(shí)長;4)在時(shí)域下求解船體液艙耦合運(yùn)動(dòng)方程;5)轉(zhuǎn)換整船運(yùn)動(dòng)規(guī)律到液艙,實(shí)時(shí)給出液艙晃蕩的邊值條件;6)進(jìn)行液艙流體晃蕩數(shù)值計(jì)算;7)傳遞液艙激勵(lì)力(矩)至船體液艙耦合運(yùn)動(dòng)方程,轉(zhuǎn)到步驟4)或結(jié)束.3通過計(jì)算和分析液體艙的運(yùn)動(dòng)及其與船體的運(yùn)動(dòng)連接3.1產(chǎn)品模擬結(jié)果液艙流體晃蕩力對船體運(yùn)動(dòng)直接產(chǎn)生影響,因此非線性液艙流體晃蕩問題的準(zhǔn)確計(jì)算是船舶與液艙晃蕩耦合運(yùn)動(dòng)系統(tǒng)模擬的前提.基于2.2節(jié)中所述時(shí)域邊界元法,對2種不同尺寸的方形液艙進(jìn)行了編程計(jì)算模擬,并將計(jì)算結(jié)果與實(shí)驗(yàn)、計(jì)入粘性的CFD結(jié)果進(jìn)行了比較.3.1.1自由土地自由面波面抬高時(shí)歷及與試驗(yàn)結(jié)果的比較方形液艙模型尺寸為1m×1m×1m,所裝的液體深度為0.5m,本文對3種工況分別進(jìn)行了模擬,其中模型液艙受到的激勵(lì)分別為2個(gè)單自由度振蕩和一個(gè)組合振蕩,具體見表1.圖3(a)中所示的是工況1中液艙自由面最右端中心點(diǎn)的波面升高時(shí)歷及其與Faltinsen的試驗(yàn)結(jié)果的對比,該工況所選擇的振蕩頻率接近液艙的固有頻率.可以看到,在自由液面出現(xiàn)破碎之前數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好.圖3(b)、(c)分別是工況2、3中液艙自由液面最右端中心點(diǎn)處的波面升高時(shí)歷,以及考慮了粘性影響的CFD計(jì)算結(jié)果.可以看出2工況中勢流計(jì)算結(jié)果和考慮了粘性的真實(shí)流體的模擬結(jié)果吻合良好,兩工況模擬計(jì)算過程中皆未出現(xiàn)自由液面破碎現(xiàn)象.3.1.2模擬計(jì)算實(shí)際用來進(jìn)行壓力驗(yàn)證的方形液艙模型尺寸為1.2m×1.2m×0.6m,其液體裝載深度為0.36m,受到的激勵(lì)為橫蕩運(yùn)動(dòng):x=0.015sin(2.475t).記液艙右端壁面中線上2點(diǎn)為P1、P2,其中P1距離底面0.3m,P2距離底面0.426m.模擬計(jì)算所得P1、P2壓力時(shí)歷及其與實(shí)驗(yàn)值的比較,見圖4所示.由圖可以看出,液艙晃蕩的數(shù)值程序計(jì)算結(jié)果與實(shí)驗(yàn)值吻合良好.上述計(jì)算結(jié)果說明基于時(shí)域勢流理論的邊界元法在不出現(xiàn)自由液面破碎等強(qiáng)非線性現(xiàn)象時(shí)可以較準(zhǔn)確的進(jìn)行液艙流體晃蕩時(shí)域模擬,從而為船體與液艙流體晃蕩耦合運(yùn)動(dòng)的準(zhǔn)確計(jì)算奠定了基礎(chǔ).3.2根據(jù)船體運(yùn)動(dòng)與液體艙傾斜的耦合時(shí)間間隔值3.2.1模型的建立和求解模型試驗(yàn)是在中國船舶科學(xué)研究中心耐波性水池中進(jìn)行的,試驗(yàn)?zāi)Ｐ褪悄筒ㄐ詷?biāo)準(zhǔn)船模S175,模型與實(shí)船的縮尺比為1∶55,主尺度參數(shù)見表2,型線圖船模加載的方形液艙尺寸和安裝位置見圖5.液艙長600mm,寬300mm,高250mm,艙內(nèi)液體的深度為125mm液艙位于第9站到13站,重心位置與船模重心重合.如上文所述,求解船體與液艙流體晃蕩耦合時(shí)域運(yùn)動(dòng),需預(yù)先得到加載液艙的S175船模的頻域計(jì)算結(jié)果,進(jìn)而利用脈沖響應(yīng)函數(shù)方法獲得時(shí)延函數(shù).鑒于篇幅,本文只給出部分阻尼系數(shù)和時(shí)延函數(shù)計(jì)算結(jié)果,圖6為橫搖頻率阻尼系數(shù)以及橫搖和縱搖的時(shí)延函數(shù).耦合運(yùn)動(dòng)方程中的非線性橫搖阻尼系數(shù)b41和b42是根據(jù)實(shí)驗(yàn)所得到的模型自搖衰減曲線通過最小二乘法得到的.3.2.2高頻入射波激勵(lì)下液艙流體搖蕩特性圖7給出了迎浪工況下載液S175船?？v搖和垂蕩的運(yùn)動(dòng)時(shí)歷.可以看出無論是高頻還是低頻的入射波激勵(lì),液艙流體晃蕩對船舶縱向運(yùn)動(dòng)的影響不大.圖8給出了迎浪工況下,數(shù)值和實(shí)驗(yàn)所得到的S175船?？v搖RAO,從圖中可以看出數(shù)值計(jì)算與模型實(shí)驗(yàn)的結(jié)果吻合較好,同時(shí)也可以看出迎浪工況下液艙晃蕩對船舶縱搖運(yùn)動(dòng)的影響不大.圖9分別給出了橫浪工況下載液S175船模橫搖和垂蕩的運(yùn)動(dòng)時(shí)歷.后者處于S175船模的橫搖共振頻率附近.可以看出對于遠(yuǎn)離自振頻率的高頻入射波激勵(lì),液艙流體晃蕩對船體橫向運(yùn)動(dòng)的影響不大,而對于自振頻率附近的低頻入射波激勵(lì),液艙流體晃蕩明顯減小了船體橫搖運(yùn)動(dòng)的幅值.出現(xiàn)這種現(xiàn)象的原因在于不同入射波激勵(lì)下,液艙流體晃蕩力和波浪誘導(dǎo)力相比有數(shù)量級(jí)上的差別和相位差.在高頻入射波激勵(lì)時(shí),晃蕩力遠(yuǎn)小于波浪誘導(dǎo)力,對船體運(yùn)動(dòng)的影響有限;而在低頻入射波激勵(lì)時(shí),液晃蕩力與波浪誘導(dǎo)力處于同一數(shù)量級(jí),而且存在180°左右的相位差,在耦合運(yùn)動(dòng)的計(jì)算中二者相互抵消,明顯減小了S175船模的橫搖幅值.圖10給出了橫浪工況下,數(shù)值和實(shí)驗(yàn)所得到的S175船模橫搖RAO,從圖中可以看出數(shù)值計(jì)算與模型實(shí)驗(yàn)的結(jié)果吻合較好,同時(shí)液艙流體晃蕩使得S175船模的橫浪共振頻率區(qū)間發(fā)生了偏移.由于對液艙流體晃蕩問題采用了基于時(shí)域勢流理論的邊界元法求解,無法很好處理液艙大幅度晃蕩時(shí)所出現(xiàn)的液面波隨等強(qiáng)非線性現(xiàn)象,耦合運(yùn)動(dòng)的計(jì)算無法達(dá)到穩(wěn)態(tài),圖11給出了加載液艙后S175船模共振頻率附近的某個(gè)頻率下橫搖運(yùn)動(dòng)的時(shí)歷,可以看出S175船模橫搖運(yùn)動(dòng)幅值在不斷增加,當(dāng)幅值達(dá)到一定程度之后計(jì)算便無法進(jìn)行下去.此外較低的共振頻率對應(yīng)的波浪周期較長,S175船模達(dá)到穩(wěn)態(tài)運(yùn)動(dòng)所需時(shí)間也隨之增加,計(jì)算誤差的積累也是導(dǎo)致計(jì)算發(fā)散的原因之一.如果能夠計(jì)算到穩(wěn)態(tài),其實(shí)際運(yùn)動(dòng)幅值預(yù)計(jì)會(huì)比圖107中所示的幅值要大.所以在圖10中,加載液艙的S175船模的橫搖RAO的數(shù)值結(jié)果在共振頻率附近的值比實(shí)驗(yàn)值要小.4非線性動(dòng)力系統(tǒng)的穩(wěn)定性本文對船舶在波浪上運(yùn)動(dòng)和液艙液體晃蕩,即對船體內(nèi)外流場問題均采用了勢流理論方法求解,其中波浪中船體水動(dòng)力和時(shí)延函數(shù)采用切片法和脈沖響應(yīng)函數(shù)法計(jì)算獲得,艙內(nèi)液體非線性晃蕩采用時(shí)域邊界元法計(jì)算,并最終建立了在波浪中船體與液艙流體晃蕩耦合的運(yùn)動(dòng)方程.文中就S175加載方形液艙在迎浪、橫浪等不