波浪中船體與液艙晃蕩耦合運動的時域數值計算

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

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

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

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

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

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

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

圖1 船體運動坐標系Fig.1 The coordinate system for ship motion

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

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

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

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

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

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

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

圖2 液艙坐標系示意Fig.2 The coordinate system of the tank

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

其運動學邊界條件:

壁面條件:

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

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

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

1.3 脈沖響應函數方法的船舶運動時域計算

在線性系統中，任意激勵可以寫為脈沖響應函數和激勵的卷積積分形式［10］:

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

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

時延函數與無窮大遭遇頻率附加質量可由以下頻域轉時域方法獲得:

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

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

式中，周期Ti遠大于波浪激勵周期.

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

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

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

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

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

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

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

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

6)進行液艙流體晃蕩數值計算;

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

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

3.1 液艙晃蕩數值模擬與驗證

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

3.1.1 液艙晃蕩數值模擬

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

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

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

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

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

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

用來進行壓力驗證的方形液艙模型尺寸為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所示.由圖可以看出，液艙晃蕩的數值程序計算結果與實驗值吻合良好.

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

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

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

3.2.1 模型試驗與時延函數

模型試驗［11］是在中國船舶科學研究中心耐波性水池中進行的，試驗模型是耐波性標準船模S175，模型與實船的縮尺比為1∶55，主尺度參數見表2，型線圖船模加載的方形液艙尺寸和安裝位置見圖5.液艙長600 mm，寬300 mm，高250 mm，艙內液體的深度為125 mm液艙位于第9站到13站，重心位置與船模重心重合.

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

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

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

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

3.2.2 耦合計算結果與分析

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

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

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

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

4 結論

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

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

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

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

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

［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，FALTINSEN 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，FALTINSEN O M.Ship motions and sea loads［J］.Trans.of Society of Naval Architects and Marine Engineers，1970，78:250-287.

［9］劉應中，繆國平.船舶在波浪上的運動理論［M］.上海:上海交通大學出版社，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］.江蘇:江蘇科技大學，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.