石留斌 李廷秋 蘇 焱

(武漢理工大學(xué)船海與能源動(dòng)力工程學(xué)院 武漢 430063)

0 引 言

隨著液化石油氣(LPG)船舶，液化天然氣(LNG)船舶及超大型液化天然氣(VLGC)船舶的需求不斷增加，船舶液艙晃蕩現(xiàn)象逐漸成為研究熱點(diǎn).

Faltinsen[1]采用解析、數(shù)值,以及實(shí)驗(yàn)方法對(duì)液艙晃蕩運(yùn)動(dòng)問題進(jìn)行了系統(tǒng)研究.Ebrahimian等[2]采用邊界元方法確定帶有擋板的軸對(duì)稱容器中對(duì)稱和非對(duì)稱晃蕩固有頻率和振型.Cho等[3]運(yùn)用有限元方法(FEM)，計(jì)算分析了在不同充液率和激勵(lì)振幅時(shí)二維液艙晃蕩運(yùn)動(dòng)，并將數(shù)值計(jì)算結(jié)果與線性理論計(jì)算結(jié)果進(jìn)行了比較.Wu等[4]基于二維Navier-Stokes方程運(yùn)用有限差分法(FDM)，研究在縱蕩和橫蕩耦合激勵(lì)條件下，內(nèi)部擋板結(jié)構(gòu)對(duì)晃蕩運(yùn)動(dòng)的影響并通過實(shí)驗(yàn)驗(yàn)證數(shù)值結(jié)果的可靠性.衛(wèi)志軍等[5]采用光滑粒子水動(dòng)力學(xué)法(SPH)對(duì)二維矩形液艙晃蕩運(yùn)動(dòng)進(jìn)行研究，結(jié)果表明：SPH方法可以很好地模擬液體晃蕩時(shí)水躍、破碎等自由液面的大變形運(yùn)動(dòng).陳翔等[6]將移動(dòng)粒子半隱式法(MPS)與圖形處理器(GPU)并行加速技術(shù)相結(jié)合，對(duì)LNG船舶液艙晃蕩進(jìn)行了數(shù)值模擬，并將LNG船舶液艙與方型液艙的晃蕩運(yùn)動(dòng)進(jìn)行對(duì)比.結(jié)果表明：在高充液率下LNG型液艙可以有效地減小晃蕩運(yùn)動(dòng)幅值和壁面砰擊壓力，但在中低充液率下，LNG型液艙則會(huì)加劇晃蕩運(yùn)動(dòng).王慶豐等[7]采用去奇異邊界元方法，在時(shí)域內(nèi)建立液艙內(nèi)不規(guī)則晃蕩運(yùn)動(dòng)數(shù)值模型，模擬了規(guī)則激勵(lì)下液艙晃蕩運(yùn)動(dòng)并與解析解進(jìn)行對(duì)比驗(yàn)證其準(zhǔn)確性，在此基礎(chǔ)上完成了非規(guī)則激勵(lì)下液艙晃蕩運(yùn)動(dòng)的數(shù)值模擬.寧德志等[8]采用完全非線性邊界條件的時(shí)域數(shù)學(xué)模型對(duì)縱搖容器中的液體晃蕩運(yùn)動(dòng)進(jìn)行研究，結(jié)果表明：相對(duì)于偶數(shù)階固有頻率來說，液體晃蕩運(yùn)動(dòng)對(duì)奇數(shù)階固有頻率更為敏感.薛米安等[9]基于振動(dòng)臺(tái)對(duì)液艙晃蕩運(yùn)動(dòng)進(jìn)行實(shí)驗(yàn)研究，發(fā)現(xiàn)在波浪破碎等強(qiáng)非線性作用下，實(shí)際一階共振頻率大于理論推導(dǎo)的一階固有頻率.楊志勛等[10]采用縮尺比為1∶1、1∶2、1∶3，載液率為20%的系列GTT液艙，進(jìn)行長(zhǎng)時(shí)間激勵(lì)下二維液艙晃蕩運(yùn)動(dòng)尺度效應(yīng)研究，結(jié)果表明：氣液密度比是導(dǎo)致模型和原型結(jié)果之間存在較大差異的重要原因.

文中采用高精度全非線性Boussinesq型方程建立淺水液艙晃蕩運(yùn)動(dòng)數(shù)值模型，并分析其在非規(guī)則外部激勵(lì)作用下液艙內(nèi)部晃蕩運(yùn)動(dòng)特征，對(duì)深入理解非線性淺水液艙晃蕩運(yùn)動(dòng)機(jī)理具有重要意義.

1 數(shù)值計(jì)算模型

選取長(zhǎng)度為l、靜水水深為h的二維矩形液艙，見圖1.慣性坐標(biāo)系oyz原點(diǎn)位于靜水面中點(diǎn)，固定于大地上，z軸豎直向上.假定矩形水箱內(nèi)為無旋無黏不可壓縮流體，流場(chǎng)中速度勢(shì)Φ滿足下列條件:

圖1 二維液艙晃蕩示意圖

2Φ=0

(1)

(2)

ηt-Φz(mì)+ηyΦy=0z=η

(3)

(4)

隨體坐標(biāo)系固定于液艙上，隨液艙一起運(yùn)動(dòng).對(duì)于自由表面條件，將慣性坐標(biāo)系下定義的時(shí)間導(dǎo)數(shù)通過以下表達(dá)式轉(zhuǎn)化為隨體坐標(biāo)系下的時(shí)間導(dǎo)數(shù).

(5)

(6)

故隨體坐標(biāo)系下的自由表面運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)邊界條件為

ηt-v·η-Φz(mì)+ηyΦy=0z=η

(7)

(8)

本節(jié)后續(xù)方程均定義在隨體坐標(biāo)系下.

總速度勢(shì)Φ分為兩個(gè)部分:Φ=φ+φ，其中速度勢(shì)特解φ滿足拉普拉斯方程和物面不可穿透條件.若僅考慮橫蕩運(yùn)動(dòng)和升沉運(yùn)動(dòng),那么

φ=yvb+zvc

(9)

另一部分速度勢(shì)φ采用Boussinesq方程進(jìn)行求解.

(10)

(11)

代入隨體坐標(biāo)系下自由表面邊界條件可得：

ηt=v·η+Φz(mì)-ηy·Φy=

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

將式(19)和式(20)聯(lián)立可得矩陣如下.

(21)

(22)

2 非規(guī)則激勵(lì)計(jì)算方法

采用北海聯(lián)合波浪譜(JONSWAP)，其波能譜密度函數(shù)S(ω)表達(dá)式為

(23)

式中：

(24)

其中：ωp為JONSWAP譜圖像峰值對(duì)應(yīng)頻率，當(dāng)頻率ω<ωp時(shí)，參數(shù)σ=0.07；當(dāng)頻率ω>ωp時(shí)，參數(shù)σ=0.09.參數(shù)γ取標(biāo)準(zhǔn)值3.3,參數(shù)Hs為有義波高,系數(shù)α的取值要保證等式(25)成立.橫坐標(biāo)為頻率ω，縱坐標(biāo)為譜密度S.

(25)

不規(guī)則激勵(lì)由規(guī)則激勵(lì)疊加獲得，其位移和速度表達(dá)式為

(26)

(27)

式中：A為幅值；ψ為隨機(jī)相位.

A與譜密度函數(shù)S(ω)存在以下關(guān)系：

(28)

式中：Δω為頻率間隔.

考慮到開始階段速度變化過大，添加漸進(jìn)函數(shù)f(t)，對(duì)前五個(gè)譜峰周期進(jìn)行修正，以保證計(jì)算結(jié)果的穩(wěn)定性.

(29)

t0=10π/ωp

(30)

式中：t0為前五個(gè)譜峰周期.

則有最終結(jié)果如下.

(31)

(32)

艙內(nèi)液體一階固有頻率(ω1=2.04 rad/s，根據(jù)色散關(guān)系式求得)作為譜峰頻率，則有譜峰頻率ωp=ω1=2.04 rad/s.頻率ω取值范圍為0.01～5 rad/s，頻率間隔Δω=0.01 rad/s.參數(shù)γ取標(biāo)準(zhǔn)值3.3，有義波高Hs=0.012 m.波能譜密度函數(shù)S(ω)隨頻率變化(見圖2)，不規(guī)則速度激勵(lì)圖像見圖3，并選取其中t=0～30 s表明漸進(jìn)函數(shù)f(t)在前五個(gè)周期的緩沖效果(見圖4).

圖2 波能譜密度函數(shù)S(ω)隨頻率變化圖

圖3 不規(guī)則速度激勵(lì)隨時(shí)間變化圖

圖4 不規(guī)則速度激勵(lì)中添加與不添加漸進(jìn)函數(shù)f(t)對(duì)比圖

3 數(shù)值模擬結(jié)果

3.1 規(guī)則激勵(lì)驗(yàn)證

為驗(yàn)證該數(shù)值模型的有效性，在相同工況下與文獻(xiàn)結(jié)果(采用固有頻率模塊疊加求解方法)進(jìn)行對(duì)比.二維矩形液艙長(zhǎng)度L=1.175 m，水深h=0.06 m，受到長(zhǎng)度方向簡(jiǎn)諧激勵(lì)，激勵(lì)幅值A(chǔ)=0.003 9 m.基于線性色散關(guān)系式獲得液艙晃蕩運(yùn)動(dòng)一階固有頻率ω1=2.04 rad/s，激勵(lì)頻率按文獻(xiàn)選取ωp1=0.98ω1，ωp2=1.02ω1，ωp3=1.08ω1.

圖5為激勵(lì)頻率分別為ωp1，ωp2，ωp3時(shí)二維液艙壁面處自由表面高度變化對(duì)比圖，橫、縱坐標(biāo)分別為時(shí)間與激勵(lì)周期之比t/T和自由表面高度與水深之比η/h.在激勵(lì)頻率ωp1=0.98ω1時(shí)，液艙中可以觀察到三個(gè)行進(jìn)波波峰在艙壁之間來回運(yùn)動(dòng)，艙壁處自由面波高時(shí)歷曲線見圖5a).在ωp2=1.02ω1和ωp3=1.08ω1時(shí)，液艙分別觀察到兩個(gè)行進(jìn)波波峰和單個(gè)行進(jìn)波波峰，自由面波高時(shí)歷曲線見圖5b)～c).由圖5可知，文中所采用數(shù)值模型能準(zhǔn)確地捕捉固有頻率附近自由液面高度變化規(guī)律.

圖5 艙壁處自由液面高度對(duì)比圖

3.2 不規(guī)則激勵(lì)數(shù)值模擬結(jié)果

基于色散關(guān)系式計(jì)算獲得液艙晃蕩運(yùn)動(dòng)一階至五階固有頻率分別為ω1=2.04 rad/s，ω2=4.03 rad/s，ω3=5.93 rad/s，ω4=7.71 rad/s，ω5=9.34 rad/s.下面對(duì)海浪譜參數(shù)對(duì)液艙晃蕩運(yùn)動(dòng)的影響進(jìn)行研究，采用matlab進(jìn)行數(shù)值模擬，程序運(yùn)行時(shí)間t=600 s，網(wǎng)格點(diǎn)數(shù)取N=39(長(zhǎng)度L=1.175 m)，時(shí)間步長(zhǎng)Δt=0.02 s，matlab自帶能量譜密度函數(shù)pwelch函數(shù)對(duì)液艙艙壁處的自由液面高程時(shí)歷曲線進(jìn)行譜分析.

選取海浪譜譜峰頻率ωp=ω1=2.04 rad/s，參數(shù)γ=3.3，有義波高Hs=0.008 m進(jìn)行研究.圖6為艙壁處自由液面高程時(shí)歷曲線圖，圖6中自由液面高程變化整體呈現(xiàn)出不規(guī)則趨勢(shì).

圖6 艙壁處自由液面高程時(shí)歷曲線圖(ωp=ω1)

改變海浪譜參數(shù)有義波高Hs的取值研究其對(duì)液艙晃蕩運(yùn)動(dòng)的影響，圖7為ωp=ω1=2.04 rad/s，參數(shù)γ=3.3，三種有義波高Hs艙壁處自由液面高程時(shí)歷曲線譜分析圖.由圖7可知:三種工況下呈現(xiàn)出相似規(guī)律.波能譜密度圖像譜峰頻率位于一階固有頻率附近，其他階固有頻率附近存在部分波能，且表現(xiàn)出波能逐階遞減的趨勢(shì)，在五階固有頻率之后波能可忽略不計(jì).三種工況對(duì)比發(fā)現(xiàn)波能大小與有義波高成正相關(guān)，有義波高越大波能越大.

圖7 不同Hs取值艙壁處自由液面高程時(shí)歷曲線譜分析圖

研究海浪譜譜峰頻率ωp與液艙晃蕩運(yùn)動(dòng)的關(guān)系，圖8為海浪譜譜峰頻率分別選取一階至五階固有頻率ωp=ω1，ωp=ω2，ωp=ω3，ωp=ω4，ωp=ω5，參數(shù)γ=3.3，有義波高Hs=0.008 m時(shí)艙壁處自由液面高程時(shí)歷曲線譜分析圖.相對(duì)于偶數(shù)階固有頻率而言，液艙晃蕩運(yùn)動(dòng)對(duì)奇數(shù)階固有頻率更加敏感，液艙晃蕩運(yùn)動(dòng)波能主要集中于一、三、五階.

圖8 不同ωp取值艙壁處自由液面高程時(shí)歷曲線譜分析圖

對(duì)譜峰頻率不等于共振頻率進(jìn)行研究，以ωp=1.8 rad/s和ωp=2.3 rad/s為例(見圖9)，與一階共振情況ωp=ω1=2.04 rad/s進(jìn)行對(duì)比，參數(shù)γ=3.3，有義波高Hs=0.008 m.結(jié)果表明， 三種工況下譜分析圖像呈現(xiàn)出類似規(guī)律，譜峰頻率位于一階固有頻率附近，其他階固有頻率附近存在部分波能，且表現(xiàn)出波能逐階遞減的趨勢(shì).

圖9 艙壁處自由液面高程時(shí)歷曲線譜分析圖

4 結(jié) 論

1) 當(dāng)譜峰頻率等于一階固有頻率時(shí)，波能大小與有義波高大小成正相關(guān)，有義波高越大波能越大.當(dāng)譜峰頻率等于或偏離一階固有頻率時(shí)，波能主要集中于各階固有頻率附近，波能在一階固有頻率附近最大且從一階至五階波能逐漸降低，五階以后可忽略不計(jì).

2) 當(dāng)譜峰頻率等于二階及以上固有頻率時(shí)，淺水液艙晃蕩波能集中于三階或五階固有頻率附近.當(dāng)譜峰頻率為二階、三階固有頻率時(shí)，淺水液艙晃蕩運(yùn)動(dòng)波能均集中于三階固有頻率附近.當(dāng)譜峰頻率為四階、五階固有頻率時(shí)，淺水液艙晃蕩運(yùn)動(dòng)波能均集中于五階固有頻率附近.