洪亮，朱仁傳，繆國平，范菊，李裕龍

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國家重點實驗室，上海 200240)

載液船舶航行時，液艙隨船體在波浪中運動使其內(nèi)部的液體產(chǎn)生晃蕩，與此同時液艙晃蕩力作用于船體，使得船舶的運動姿態(tài)發(fā)生變化.近年來，隨著海洋油氣資源開發(fā)力度的加大，LNG和LPG等特種液貨船型的研制和廣泛應(yīng)用，液艙與船舶之間的作用越來越受人們重視，劇烈的艙內(nèi)液體晃蕩可能會給船體帶來巨大的危害，特別是液艙部分裝載時，大幅晃蕩引起的砰擊壓力對船體結(jié)構(gòu)可能造成嚴(yán)重的破壞.因此，在載液船舶設(shè)計階段液艙晃蕩與船舶運動之間的耦合作用也是需要考慮的重要因素之一.對于船體運動與液艙晃蕩耦合作用的問題，國內(nèi)外已有許多學(xué)者進(jìn)行了分析和研究，如 Kim［1］，Rognebakke和Faltinsen［2-3］，Newman［4］等，對液艙流體晃蕩的模擬方法大體分為線性化的頻域法(如Newman［5］)和非線性時域法(如Kim［6-7］).從現(xiàn)有的研究結(jié)果來看，線性船體運動理論方法基本能夠滿足船舶運動與液艙晃蕩之間耦合效應(yīng)分析研究的需要.

本文對船舶在波浪上運動和液艙液體晃蕩，即對船體內(nèi)外流場問題均采用了勢流理論方法求解，其中波浪中船體水動力和時延函數(shù)采用切片法和脈沖響應(yīng)函數(shù)法計算獲得，艙內(nèi)液體非線性晃蕩采用時域邊界元法計算，并最終建立了在波浪中船體與液艙流體晃蕩耦合的運動方程.文中就S175加載方形液艙在迎浪、橫浪等不同工況下液艙流體晃蕩及其與船體運動耦合分別進(jìn)行了計算模擬與驗證研究.研究中發(fā)現(xiàn)，當(dāng)不出現(xiàn)液面破碎等強非線性現(xiàn)象時，非線性時域邊界元法能夠給出較好的液艙流體晃蕩波形和壓力;S175的液艙加載50%的液體時，船體與液艙晃蕩耦合運動時歷結(jié)果能清晰地反映液艙晃蕩對船體運動的影響，運動RAO能反映出不同頻率液艙對船體運動的影響程度;載液S175船橫搖運動RAO能準(zhǔn)確給出船體有無加載液艙時共振頻率的偏移現(xiàn)象.

1 載液船舶耦合運動的基本理論

1.1 載液船舶在波浪上的運動

航行船舶的運動計算所采用的參考坐標(biāo)系為oxyz，如圖1所示，oxy平面與靜水面重合，oz軸垂直向上，船舶以定速U0航行，在波浪激勵下作六自由度運動.

圖1 船體運動坐標(biāo)系Fig.1 The coordinate system for ship motion

直接求解滿足定解條件的速度勢不是一件很容易的事，基于切片理論的STF法可以不必求解三維有航速相應(yīng)的邊值問題，只要求解船舶各個橫剖面的二維零航速水動力問題即可，STF法的詳細(xì)介紹可見文獻(xiàn)［8-9］.

若將船體橫剖面的速度勢φ沿船長方向積分即可獲得整船的附加質(zhì)量和阻尼系數(shù):

那么船體五自由度運動方程可表示如下:

式中:ω為頻率;ξj為船體j模態(tài)的運動位移;mij為包含了液艙內(nèi)液體質(zhì)量的船體廣義質(zhì)量矩陣;μij(ω)、λij(ω)為頻域附加質(zhì)量和阻尼系數(shù);Cij為船體回復(fù)力系數(shù)(ω)為作用在船體上的廣義波浪力;j= 2，3，…，6分別對應(yīng)的運動模態(tài)為橫蕩、垂蕩、橫搖、縱搖和艏搖.切片法不能計及船體縱蕩運動.

1.2 基于勢流理論的液艙晃蕩問題

采用三維勢流理論對液艙內(nèi)非定常流體晃蕩問題進(jìn)行求解，液艙流體晃蕩是采用時域邊界元法進(jìn)行模擬計算的.

液艙內(nèi)流體運動計算采用的坐標(biāo)系如圖2所示，是固定在液艙內(nèi)的局部坐標(biāo)系，其坐標(biāo)原點固定在靜止時自由液面的幾何中心，x軸正方向與液艙長度方向平行指向右端，y軸正方向與液艙寬度方向平行指向前端，z軸正方向為豎直向上，在計算過程中該坐標(biāo)系始終與液艙長寬高保持平行，其原點保持不變.假定流體不可壓縮，無粘，流動無旋，則流體運動速度勢滿足拉普拉斯方程，即控制方程滿足:

圖2 液艙坐標(biāo)系示意Fig.2 The coordinate system of the tank

三維流體晃蕩運動的流場的自由表面上動力學(xué)邊界條件為

其運動學(xué)邊界條件:

壁面條件:

式中:φ為艙內(nèi)流體運動速度勢，Ω為液艙角速度矢量，V為液艙線速度矢量，r為位置矢量，n為邊界外法向，ζ為自由面升高.

本文采用簡單格林函數(shù)法來求解Laplace方程，取1/r為格林函數(shù)，則流場中任意一點的速度勢可以表述如下:

式中:α為源點處的歐拉角，通過對空間的離散求得每個時間步中的速度勢，再通過對自由表面動力學(xué)條件進(jìn)行中心差分獲得下一個時間步的自由面條件，從而使得液艙晃蕩模擬可以在時域中反復(fù)進(jìn)行下去.本文的空間離散采用了線性面元，積分計算中源點與場點重合時積分為假性奇點，采用解析法進(jìn)行了計算，源點與場點不重合時采用了七點高斯積分公式計算.計算中根據(jù)每個時間步獲得的速度勢，通過拉格朗日積分計算獲得靜壓力與動壓力.

1.3 脈沖響應(yīng)函數(shù)方法的船舶運動時域計算

在線性系統(tǒng)中，任意激勵可以寫為脈沖響應(yīng)函數(shù)和激勵的卷積積分形式［10］:

式中:x(t)為在輸入h(t)下的系統(tǒng)響應(yīng)，F(xiàn)(t)為單位脈沖輸入下的脈沖響應(yīng)函數(shù).將以上概念推廣到船舶運動問題，時域船體運動方程有著如下形式:

式中:Mij、μij(∞)和Cij分別代表載液船體質(zhì)量，無窮大遭遇頻率的附加質(zhì)量以及船舶回復(fù)力系數(shù);Kij(τ)為時延函數(shù);(t)為作用在船體上的外部時域波浪力;b41和b42為考慮粘性作用的非線性船舶橫搖阻尼系數(shù)，可以用經(jīng)驗公式、物理實驗等方法獲得.

時延函數(shù)與無窮大遭遇頻率附加質(zhì)量可由以下頻域轉(zhuǎn)時域方法獲得:

式中:μij(∞)指頻率無窮大時船體的附加質(zhì)量.

求解運動方程時，橫蕩與首搖模態(tài)的運動可采用數(shù)值彈簧技術(shù)來抑制其慢漂現(xiàn)象.所加入的橫蕩和首搖兩種運動模態(tài)彈簧剛度表達(dá)如下:

式中，周期Ti遠(yuǎn)大于波浪激勵周期.

2 載液船舶在波浪上的時域耦合算法及計算步驟

船舶與液艙晃蕩耦合運動方程可表示為

在計算模擬過程中，船體當(dāng)前時刻的運動規(guī)律傳遞到液艙，然后液艙晃蕩模擬程序計算出該時刻液艙內(nèi)部流場的速度勢，再通過對艙壁上的壓力積分得到液艙的晃蕩力和力矩FSlosh，并添加至式(13)的右端，采用四階龍格庫塔法進(jìn)行耦合求解，此時獲得的整船運動位移又作為下一時刻液艙流體晃蕩的條件進(jìn)行計算.由于船體運動預(yù)報與液艙晃蕩是在不同坐標(biāo)系中進(jìn)行的，所以二者之間力與運動需要在這2個坐標(biāo)系中進(jìn)行轉(zhuǎn)換，轉(zhuǎn)換方法可參考文獻(xiàn)［9］.這樣反復(fù)循環(huán)便可以得到船體與液艙耦合運動的時歷，這就是時域下耦合液艙晃蕩的船舶全局運動的求解方法.具體的運算步驟如下:

1)讀取船型液艙參數(shù)及頻域水動力數(shù)據(jù);

2)根據(jù)來波工況計算波浪力時歷;

3)計算時延函數(shù)，確定計算時長;

4)在時域下求解船體液艙耦合運動方程;

5)轉(zhuǎn)換整船運動規(guī)律到液艙，實時給出液艙晃蕩的邊值條件;

6)進(jìn)行液艙流體晃蕩數(shù)值計算;

7)傳遞液艙激勵力(矩)至船體液艙耦合運動方程，轉(zhuǎn)到步驟4)或結(jié)束.

3 液艙流體晃蕩及其與船體運動耦合的計算與分析

3.1 液艙晃蕩數(shù)值模擬與驗證

液艙流體晃蕩力對船體運動直接產(chǎn)生影響，因此非線性液艙流體晃蕩問題的準(zhǔn)確計算是船舶與液艙晃蕩耦合運動系統(tǒng)模擬的前提.基于2.2節(jié)中所述時域邊界元法，對2種不同尺寸的方形液艙進(jìn)行了編程計算模擬，并將計算結(jié)果與實驗、計入粘性的CFD結(jié)果進(jìn)行了比較.

3.1.1 液艙晃蕩數(shù)值模擬

方形液艙模型尺寸為1 m×1 m×1 m，所裝的液體深度為0.5 m，本文對3種工況分別進(jìn)行了模擬，其中模型液艙受到的激勵分別為2個單自由度振蕩和一個組合振蕩，具體見表1.

表1 模型液艙的受到激勵Table 1 Working condition of tank

圖3(a)中所示的是工況1中液艙自由面最右端中心點的波面升高時歷及其與Faltinsen的試驗結(jié)果［10］的對比，該工況所選擇的振蕩頻率接近液艙的固有頻率.可以看到，在自由液面出現(xiàn)破碎之前數(shù)值計算結(jié)果與試驗結(jié)果吻合良好.

圖3 液艙右端中心點自由面波高時歷Fig.3 Time history of the free surface elevation of tank

圖3(b)、(c)分別是工況2、3中液艙自由液面最右端中心點處的波面升高時歷，以及考慮了粘性影響的CFD計算結(jié)果.可以看出2工況中勢流計算結(jié)果和考慮了粘性的真實流體的模擬結(jié)果吻合良好，兩工況模擬計算過程中皆未出現(xiàn)自由液面破碎現(xiàn)象.

3.1.2 晃蕩液艙艙壁壓力計算驗證

用來進(jìn)行壓力驗證的方形液艙模型尺寸為1.2 m×1.2 m×0.6 m，其液體裝載深度為0.36 m，受到的激勵為橫蕩運動:x=0.015sin(2.475t).

記液艙右端壁面中線上2點為P1、P2，其中P1距離底面0.3 m，P2距離底面0.426 m.模擬計算所得P1、P2壓力時歷及其與實驗值的比較，見圖4所示.由圖可以看出，液艙晃蕩的數(shù)值程序計算結(jié)果與實驗值吻合良好.

圖4 P1、P2點處壓強時歷Fig.4 Time history of the pressure at P1and P2

上述計算結(jié)果說明基于時域勢流理論的邊界元法在不出現(xiàn)自由液面破碎等強非線性現(xiàn)象時可以較準(zhǔn)確的進(jìn)行液艙流體晃蕩時域模擬，從而為船體與液艙流體晃蕩耦合運動的準(zhǔn)確計算奠定了基礎(chǔ).

3.2 船體運動與液艙晃蕩的耦合時域數(shù)值計算

3.2.1 模型試驗與時延函數(shù)

模型試驗［11］是在中國船舶科學(xué)研究中心耐波性水池中進(jìn)行的，試驗?zāi)Ｐ褪悄筒ㄐ詷?biāo)準(zhǔn)船模S175，模型與實船的縮尺比為1∶55，主尺度參數(shù)見表2，型線圖船模加載的方形液艙尺寸和安裝位置見圖5.液艙長600 mm，寬300 mm，高250 mm，艙內(nèi)液體的深度為125 mm液艙位于第9站到13站，重心位置與船模重心重合.

表2 S175實船和船模主尺度Table 2 Principal dimensions of S175 and its modelm

如上文所述，求解船體與液艙流體晃蕩耦合時域運動，需預(yù)先得到加載液艙的S175船模的頻域計算結(jié)果，進(jìn)而利用脈沖響應(yīng)函數(shù)方法獲得時延函數(shù).鑒于篇幅，本文只給出部分阻尼系數(shù)和時延函數(shù)計算結(jié)果，圖6為橫搖頻率阻尼系數(shù)以及橫搖和縱搖的時延函數(shù).耦合運動方程中的非線性橫搖阻尼系數(shù)b41和b42是根據(jù)實驗［11］所得到的模型自搖衰減曲線通過最小二乘法得到的.

圖5 S175船型線及所加載液艙的尺寸與位置Fig.5 Body plan of S175 and size and location of tank

圖6 線性橫搖阻尼系數(shù)以及橫搖，縱搖模態(tài)的時延函數(shù)Fig.6 Linear damping coefficient of roll and retardation functions of roll and pitch

3.2.2 耦合計算結(jié)果與分析

圖7給出了迎浪工況下載液S175船?？v搖和垂蕩的運動時歷.可以看出無論是高頻還是低頻的入射波激勵，液艙流體晃蕩對船舶縱向運動的影響不大.圖8給出了迎浪工況下，數(shù)值和實驗所得到的S175船?？v搖RAO，從圖中可以看出數(shù)值計算與模型實驗的結(jié)果吻合較好，同時也可以看出迎浪工況下液艙晃蕩對船舶縱搖運動的影響不大.

圖7 迎浪工況下船體縱搖與垂蕩運動時歷Fig.7 Time history of pitch and heave motion of ship in head sea

圖9分別給出了橫浪工況下載液S175船模橫搖和垂蕩的運動時歷.后者處于S175船模的橫搖共振頻率附近.可以看出對于遠(yuǎn)離自振頻率的高頻入射波激勵，液艙流體晃蕩對船體橫向運動的影響不大，而對于自振頻率附近的低頻入射波激勵，液艙流體晃蕩明顯減小了船體橫搖運動的幅值.出現(xiàn)這種現(xiàn)象的原因在于不同入射波激勵下，液艙流體晃蕩力和波浪誘導(dǎo)力相比有數(shù)量級上的差別和相位差.在高頻入射波激勵時，晃蕩力遠(yuǎn)小于波浪誘導(dǎo)力，對船體運動的影響有限;而在低頻入射波激勵時，液晃蕩力與波浪誘導(dǎo)力處于同一數(shù)量級，而且存在180°左右的相位差，在耦合運動的計算中二者相互抵消，明顯減小了S175船模的橫搖幅值.

圖8 迎浪工況下數(shù)值與實驗所得的S175船?？v搖RAOFig.8 Comparison of pitch RAO of ship in head sea by experiment and calculation

圖9 橫浪工況下船體橫搖與垂蕩運動時歷Fig.9 Time history of roll and heave motion of ship in beam sea

圖10 橫浪時S175船模橫搖RAOFig.10 Roll RAO of S175 model in beam sea

圖11 橫浪時船舶橫搖時歷(周期T=4.735 s，波幅=0.1 m)Fig.11 Time history of roll of S175 model in beam sea (T=4.735s，ζa=0.1m)

圖10給出了橫浪工況下，數(shù)值和實驗所得到的S175船模橫搖RAO，從圖中可以看出數(shù)值計算與模型實驗的結(jié)果吻合較好，同時液艙流體晃蕩使得S175船模的橫浪共振頻率區(qū)間發(fā)生了偏移.由于對液艙流體晃蕩問題采用了基于時域勢流理論的邊界元法求解，無法很好處理液艙大幅度晃蕩時所出現(xiàn)的液面波隨等強非線性現(xiàn)象，耦合運動的計算無法達(dá)到穩(wěn)態(tài)，圖11給出了加載液艙后S175船模共振頻率附近的某個頻率下橫搖運動的時歷，可以看出S175船模橫搖運動幅值在不斷增加，當(dāng)幅值達(dá)到一定程度之后計算便無法進(jìn)行下去.此外較低的共振頻率對應(yīng)的波浪周期較長，S175船模達(dá)到穩(wěn)態(tài)運動所需時間也隨之增加，計算誤差的積累也是導(dǎo)致計算發(fā)散的原因之一.如果能夠計算到穩(wěn)態(tài)，其實際運動幅值預(yù)計會比圖107中所示的幅值要大.所以在圖10中，加載液艙的S175船模的橫搖RAO的數(shù)值結(jié)果在共振頻率附近的值比實驗值要小.

4 結(jié)論

本文對船舶在波浪上運動和液艙液體晃蕩，即對船體內(nèi)外流場問題均采用了勢流理論方法求解，其中波浪中船體水動力和時延函數(shù)采用切片法和脈沖響應(yīng)函數(shù)法計算獲得，艙內(nèi)液體非線性晃蕩采用時域邊界元法計算，并最終建立了在波浪中船體與液艙流體晃蕩耦合的運動方程.文中就S175加載方形液艙在迎浪、橫浪等不同工況下液艙流體晃蕩及其與船體運動耦合分別進(jìn)行了計算模擬與驗證研究.得出如下結(jié)論:

1)當(dāng)不出現(xiàn)液面破碎等強非線性現(xiàn)象時，非線性時域邊界元法能夠給出較好的液艙流體晃蕩波形和壓力.

2)S175的液艙加載50%的液體時，船體與液艙晃蕩耦合運動時歷結(jié)果能清晰地反映液艙晃蕩對船體運動的影響.在迎浪工況下，船體運動以縱搖和垂蕩為主，由于縱向波浪誘導(dǎo)力矩較大，液艙晃蕩力矩與之相比很小，此工況下液艙流體晃蕩對于船舶縱向運動的影響不大;在橫浪工況下，船舶運動以橫搖和垂蕩為主，當(dāng)液艙晃蕩力矩與波浪誘導(dǎo)力矩數(shù)量級相同時，液艙流體晃蕩船體橫搖運動產(chǎn)生較大的影響，二者的相位差則決定了液艙流體晃蕩的影響效果(增大或減小).從運動RAO中能反映出不同頻率液艙對船體運動的影響程度.其中載液S175船橫搖運動RAO能準(zhǔn)確給出船體有無加載液艙時共振頻率的偏移現(xiàn)象.

3)本文只是針對單個液艙且其重心與船舶重心重合的情況進(jìn)行了計算分析，就已經(jīng)可以看出液艙流體晃蕩對于船體橫搖運動有著不小的影響.可以預(yù)計，當(dāng)液艙重心與船舶重心不重合甚至加載多個液艙時，液艙流體晃蕩力的變化更大，液艙流體晃蕩對船舶橫搖運動也會有著更大的影響.

4)本文方法具有較高的計算效率，數(shù)值計算與試驗結(jié)果吻合良好，為設(shè)計載液船舶或減搖水艙等的前期設(shè)計提供了快速有效的分析方法和技術(shù)手段.這種方法的局限性在于無法處理自由液面破碎等強非線性現(xiàn)象，雖然可以較好的預(yù)報船體與液艙晃蕩耦合運動的共振頻率區(qū)間，但是無法準(zhǔn)確預(yù)報共振頻率附近的運動幅值.

［1］KIM Y.A numerical study on sloshing flows coupled with ship motion——the anti-rolling tank problem［J］.Journal of Ship Research，2002，46(1):52-62.

［2］FALTINSEN O M，ROGNEBAKKE O F.Sloshing［C］// International Conference on Ship and Shipping Research，Venice，Italy，2000.

［3］ROGNEBAKKE O F，F(xiàn)ALTINSEN O M.Coupling of sloshing and ship motions［J］.Journal of Ship Research，2003，47(3):208-221.

［4］NEWMAN J N.Wave effects on vessels with internal tanks［C］//The 20th Workshop on Water Waves and Floating Bodies.Longyearbyen，Norway，2005.

［5］LEE C H，NEWMAN J N.Computation of wave effects using the panel method［M］.Cambridge:MIT Press，2005: 211-251.

［6］KIM Y，SHIN Y S，LIN W M，et al.Study on sloshing problem coupled with ship motion in waves［C］//The 8th International Conference on Numerical Ship Hydrodynamics.Busan，Korea，2003.

［7］KIM Y，SHIN Y S，LEE K H.Numerical study on slosh-induced impact pressures on three-dimensional prismatic tanks［J］.Applied Ocean Research，2004，26(5):213-226.

［8］SALVESEN N，TUCK E O，F(xiàn)ALTINSEN O M.Ship motions and sea loads［J］.Trans.of Society of Naval Architects and Marine Engineers，1970，78:250-287.

［9］劉應(yīng)中，繆國平.船舶在波浪上的運動理論［M］.上海:上海交通大學(xué)出版社，1986:151-171，28-31.

［10］FALTINSEN O M.A numerical nonlinear method of sloshing in tanks with two-dimensional flow［J］.Journal of Ship Research，1978，22(3):193-202.

［11］鄒康，繆泉明，朱仁慶.耦合液艙晃蕩的船舶運動性能研究［D］.江蘇:江蘇科技大學(xué)，2009:37-55.

ZOU Kang，MIAO Quanming，ZHU Renqing.Motion Performance Research of Ship with Sloshing Tanks［D］.Jiangsu:Jiangsu University of Science And Technology，2009: 37-55.