非規(guī)則激勵下淺水液艙晃蕩數(shù)值模擬

石留斌 李廷秋 蘇 焱

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

0 引 言

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

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

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

1 數(shù)值計算模型

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

圖1 二維液艙晃蕩示意圖

2Φ=0

(1)

(2)

ηt-Φz+ηyΦy=0z=η

(3)

(4)

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

(5)

(6)

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

ηt-v·η-Φz+ηyΦy=0z=η

(7)

(8)

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

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

φ=yvb+zvc

(9)

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

(10)

(11)

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

ηt=v·η+Φz-ηy·Φy=

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

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

(21)

(22)

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

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

(23)

式中：

(24)

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

(25)

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

(26)

(27)

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

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

(28)

式中：Δω為頻率間隔.

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

(29)

t0=10π/ωp

(30)

式中：t0為前五個譜峰周期.

則有最終結(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ī)則速度激勵圖像見圖3，并選取其中t=0～30 s表明漸進(jìn)函數(shù)f(t)在前五個周期的緩沖效果(見圖4).

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

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

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

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

3.1 規(guī)則激勵驗證

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

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

圖5 艙壁處自由液面高度對比圖

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

基于色散關(guān)系式計算獲得液艙晃蕩運(yùn)動一階至五階固有頻率分別為ω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.下面對海浪譜參數(shù)對液艙晃蕩運(yùn)動的影響進(jìn)行研究，采用matlab進(jìn)行數(shù)值模擬，程序運(yùn)行時間t=600 s，網(wǎng)格點數(shù)取N=39(長度L=1.175 m)，時間步長Δt=0.02 s，matlab自帶能量譜密度函數(shù)pwelch函數(shù)對液艙艙壁處的自由液面高程時歷曲線進(jìn)行譜分析.

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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