黃 碩，段文洋，游亞戈，姜金輝，王文勝

（1.中國科學(xué)院可再生能源重點(diǎn)實(shí)驗(yàn)室，廣東廣州510000；2.中國科學(xué)院廣州能源研究所，廣東廣州510000；3.哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001；4.上海船舶運(yùn)輸科學(xué)研究所，上海200135）

液艙晃蕩與船體非線性時(shí)域耦合運(yùn)動(dòng)計(jì)算

黃 碩1，2，3，段文洋3，游亞戈1，2，姜金輝4，王文勝1，2

（1.中國科學(xué)院可再生能源重點(diǎn)實(shí)驗(yàn)室，廣東廣州510000；2.中國科學(xué)院廣州能源研究所，廣東廣州510000；3.哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001；4.上海船舶運(yùn)輸科學(xué)研究所，上海200135）

LNG船等液貨船的研制和應(yīng)用中，為了研究非線性水動(dòng)力對液艙晃蕩及船體運(yùn)動(dòng)的影響，改進(jìn)了規(guī)則波及不規(guī)則波中液艙晃蕩與船體非線性耦合運(yùn)動(dòng)的時(shí)域預(yù)報(bào)技術(shù)。通過改變液艙物面邊界條件模擬能量耗散，對液艙內(nèi)部流動(dòng)采用全非線性模型，船體外流場采用非線性瞬時(shí)物面條件的時(shí)域格林函數(shù)邊界元模型，并應(yīng)用基于組合B樣條的高階面元法求解邊值問題和迭代半隱式時(shí)域步進(jìn)法對內(nèi)外流場和船體運(yùn)動(dòng)同時(shí)求解。通過對比，驗(yàn)證了這種考慮能量耗散的方法的有效性，研究了非線性水動(dòng)力中各種因素的影響，得出能量耗散系數(shù)μ的選取規(guī)律。

液艙晃蕩；能量耗散；全非線性模型；時(shí)域耦合運(yùn)動(dòng)；LNG船；B樣條；迭代半隱式時(shí)域步進(jìn)法

液艙部分裝載時(shí)艙內(nèi)液體大幅晃蕩可能會(huì)給船體帶來巨大的危害。因此在時(shí)域內(nèi)直接求解液艙晃蕩與船體非線性耦合運(yùn)動(dòng)問題具有很大的吸引力。目前國內(nèi)外學(xué)者就該問題分別應(yīng)用線性頻域及非線性理論對液艙晃蕩進(jìn)行模擬［1-2］，對船體運(yùn)動(dòng)大多采用線性方法進(jìn)行分析［3-4］。但試驗(yàn)發(fā)現(xiàn)即使入射波是微幅線性的，某些頻率下巨大的晃蕩沖擊壓力會(huì)導(dǎo)致船體運(yùn)動(dòng)平衡位置及濕表面的大幅變化，作用在船體上的流體動(dòng)力載荷具有明顯的非線性特點(diǎn)。

為求解這種耦合問題，首先要對液艙晃蕩問題進(jìn)行精確高效的計(jì)算。部分學(xué)者應(yīng)用Navier-Stokes方程進(jìn)行求解［5］。但其計(jì)算時(shí)間較長，所以應(yīng)用加入人工能量耗散［6］的勢流理論。文獻(xiàn)［7-8］建立了一種改進(jìn)的物面邊界條件模擬能量耗散，但只對線性液藏晃蕩問題進(jìn)行了研究。

鑒于上述問題，首先將文獻(xiàn)［8］建立的能量耗散條件擴(kuò)展到全非線性液艙晃蕩數(shù)值模擬中。注意到即使入射波是線性的，船體搖蕩引起艙液晃蕩的非線性也非常顯著。因此，抓住該耦合問題的主要矛盾，對船體外部流動(dòng)采用線性自由面和非線性瞬時(shí)物面條件的時(shí)域分析方法，對液艙內(nèi)部流動(dòng)采用改進(jìn)的全非線性模型，利用B樣條高階面元法對邊界積分方程進(jìn)行離散［9］。擴(kuò)展文獻(xiàn)［10］提出的迭代半隱式時(shí)域步進(jìn)法（ISITIMFB）對內(nèi)外流場和船體運(yùn)動(dòng)同時(shí)求解。在規(guī)則波及不規(guī)則波作用下對液艙晃蕩與船體非線性耦合運(yùn)動(dòng)進(jìn)行計(jì)算，通過與文獻(xiàn)［11-12］中頻域及試驗(yàn)結(jié)果的對比，研究能量耗散條件的作用效果，得到能量耗散系數(shù)μ的選取規(guī)律。

1 數(shù)學(xué)模型

1.1 全非線性液艙晃蕩

1.1.1 速度勢定解條件

為方便求解，對全非線性液艙晃蕩的模擬在動(dòng)坐標(biāo)系oTyTzT下進(jìn)行，如圖1。其與大地坐標(biāo)系OYZ關(guān)系為

式中：α為歐拉角，T為轉(zhuǎn)化矩陣。

圖1 液艙坐標(biāo)系Fig.1 Coordinate system of the tank

假定流體為不可壓縮無粘無旋，在計(jì)算域SΩ內(nèi)速度勢φ滿足的定解條件［13］為

式中：ζb是oTyTzT下的自由表面波高，V和W是速度U在oTyT和oTzT方向的分量，單位法向量n指向流體域外部，g是重力加速度，Tij是T中的系數(shù)，SB（t）是液艙瞬時(shí)濕表面。

采用簡單格林函數(shù)法，將問題轉(zhuǎn)化為

式中：C為P點(diǎn)處流場邊界的兩條半切線之間的夾角在域內(nèi)的部分，S為除去P點(diǎn)的整個(gè)流域邊界。

值得注意的是，速度勢定解問題是在動(dòng)坐標(biāo)系下求解，因此流體對艙壁壓力可寫為

1.1.2 能量耗散的近似

為了在勢流理論模型式（2）中引入粘性效應(yīng)，文獻(xiàn)［7］在頻域分析中認(rèn)為能量耗散主要發(fā)生在壁面邊界層處，因此對液艙物面邊界條件進(jìn)行改進(jìn)。文獻(xiàn)［8］首先將該方法應(yīng)用到時(shí)域計(jì)算中，得到直接的時(shí)域表達(dá)形式為

式中：μ是阻尼系數(shù)，SB是液艙平均濕表面。Un反映能量耗散，是與式（2）中相同的物面法向速度。

數(shù)值研究發(fā)現(xiàn)式（5）中的系數(shù)μ依賴于液艙形狀及充液比例。這需要針對不同的工況選擇不同的系數(shù)μ，在實(shí)際應(yīng)用中是相對困難的。為了避免式（5）帶來的問題，文獻(xiàn)［8］將式（5）改進(jìn)為

在全非線性液艙晃蕩數(shù)值模擬中，注意到物面條件在瞬時(shí)濕表面SB（t）上滿足，并將式（1）代入式（6）得到改進(jìn)物面條件在動(dòng)作標(biāo)系oTyTzT下的表達(dá)式：

1.1.3 動(dòng)坐標(biāo)系下?φ／?t的定解條件

為解決?φ／?t的計(jì)算難題，仿照速度勢φ的求解過程，把?φ／?t看成一個(gè)調(diào)和函數(shù)，建立了?φ／?t的在大地坐標(biāo)系OYZ下的定解條件［15］：

由于求解是建立在液艙動(dòng)坐標(biāo)系oTyTzT下，因此需建立動(dòng)坐標(biāo)系下?φ／?t的定解條件，得

動(dòng)坐標(biāo)系下?φ／?t的自由面條件已由式（2）給出，至此建立了動(dòng)坐標(biāo)系下?φ／?t的定解條件。

1.2 液艙晃蕩與船體非線性耦合運(yùn)動(dòng)

考慮任意形狀帶有液艙的船體在自由水表面作大幅搖蕩的運(yùn)動(dòng)問題，定義3個(gè)坐標(biāo)系如圖2。

1）大地坐標(biāo)系OYZ，表述入射波及船體周圍的擾動(dòng)流場，建立船體的運(yùn)動(dòng)方程。

2）船體動(dòng)坐標(biāo)系oyz。表述船體表面，其隨船體一起平移和旋轉(zhuǎn)。t＝0時(shí)刻，o與O重合。搖蕩時(shí)，船體重心G′的位置矢量為S＝（η1，η2，η3）。

3）液艙局部動(dòng)坐標(biāo)系oTyTzT。t＝0時(shí)刻，oTyT軸與液艙靜水面重合，oTzT軸通過船體的重心G，且oT在液艙內(nèi)部水線面中心處。

圖2 坐標(biāo)系的建立Fig.2 Coordinate definition

未擾動(dòng)時(shí)，以上3個(gè)坐標(biāo)系互相平行，流場邊界上的單位法向量均指向流體域外部。

理想流體假設(shè)下，存在 φ ＝φ0＋φp。 其中，φ0為入射勢，φp為擾動(dòng)勢。入射波的一階速度勢為

式中：ζa是入射波波幅，ω是入射波頻率，k0是波數(shù)。

擾動(dòng)勢φp的物面非線性定解條件寫為

滿足線性自由面條件的時(shí)域Green函數(shù)為

式中：δ（t）為 delta函數(shù)，H t()為單位階躍函數(shù)，τ≤t。

求解擾動(dòng)勢φp的物面非線性分布源模型［16］為

式中：B（t）和A（t）為浮體與自由面Z＝0左右兩交點(diǎn)的橫坐標(biāo)，和為點(diǎn)B和A的橫坐標(biāo)對時(shí)間導(dǎo)數(shù)，｜A－B表示兩交點(diǎn)對應(yīng)的值相減。

液艙劇烈晃蕩與船體非線性耦合運(yùn)動(dòng)需在船體重心處建立運(yùn)動(dòng)方程：

式中：除了有外流場水動(dòng)力（矩）F（N），橫向彈簧回復(fù)力Kη2，阻尼力，船體重力及橫搖阻尼力矩Nc外，還有液艙晃蕩水動(dòng)力（矩）Fsloshing（Nsloshing）的作用；K為彈簧剛度系數(shù)；C為阻尼系數(shù)；m為裸船質(zhì)量。本文應(yīng)用文獻(xiàn)［11］的一種半解析的方法估算非線性橫搖阻尼力矩：

式中：CD為無量綱的阻尼系數(shù)，B為船寬，L為船長。

為避免每一時(shí)間步逐點(diǎn)對時(shí)間微分，再對?φ／?t積分導(dǎo)致的時(shí)域模擬失效。采用文獻(xiàn)［16］提出的先積分后微分的船體大幅運(yùn)動(dòng)水動(dòng)力計(jì)算公式：

由于在耦合運(yùn)動(dòng)計(jì)算中船體運(yùn)動(dòng)需時(shí)時(shí)傳遞至液艙，同時(shí)液艙的晃蕩力（矩）亦需時(shí)時(shí)作用于船體，故還需對這2組物理量在不同坐標(biāo)系下轉(zhuǎn)換。

2 數(shù)值模擬方法

2.1 基于組合B樣條的高階面元法

對于如式（2）描述的定解問題的求解，本文發(fā)展了一種基于B樣條的高階面元法。未知的速度勢φ和物體幾何形狀可由3次B樣條曲線定義：式中：di（t）和dj（t）分別是自由面和瞬時(shí)物面上的控制點(diǎn)，n1和 n2是相應(yīng)的控制點(diǎn)數(shù)目，φni（t）和φj（t）是相應(yīng)控制點(diǎn)上的速度勢法向?qū)?shù)和速度勢。

在自由面和物面交點(diǎn) P處采用重節(jié)點(diǎn)技術(shù)。應(yīng)用B樣條端點(diǎn)條件，交點(diǎn)處邊界及連續(xù)性條件為

將式（22）代入式（3）可得到以自由面及物面上諸控制點(diǎn)di（t）和dj（t）對應(yīng)的φj（t）和φnj（t）為未知數(shù)的線性方程，求解過程參見文獻(xiàn)［17］。

2.2 運(yùn)動(dòng)方程解耦

本文對ISITIMFB方法［18］進(jìn)行擴(kuò)展，在每一時(shí)間步，船體運(yùn)動(dòng)時(shí)時(shí)傳遞給液艙，求解該時(shí)刻液艙內(nèi)部流場，再對艙壁上的壓力積分計(jì)算出液艙的晃蕩力添加至船體運(yùn)動(dòng)方程，進(jìn)而對耦合運(yùn)動(dòng)解耦。

首先假定在t＝tn－1時(shí)刻，所有計(jì)算都已完成。tn時(shí)刻的船體及流體速度可按如下得到：

1）基于最小二乘法，由前面諸時(shí)刻已知的加速度預(yù)估tn時(shí)刻的船體加速度和液艙晃蕩能量耗散條件中的 (?φT／?t )n（0），應(yīng)用Adams-Moulton方法［19］預(yù)估船體速度Un（0）：

2）將當(dāng)前預(yù)估的Un（0）及 （?φT／?t）n（0）代入液艙和船體水動(dòng)力定解條件中求解φp和φT。

3）計(jì)算液艙內(nèi)部流體自由面表面速度勢及流體速度，速度由當(dāng)前的Un（0）轉(zhuǎn)換得到。

4）求解船體的受力及加速度：

①應(yīng)用U·n（k－1）和 Un（k－1）求解艙內(nèi)物面及自由面上的 ，n（k－1）表示t時(shí)刻的第k次n迭代；

②計(jì)算船體受力Fn（k），船體加速度及速度為

Mb是裸船慣性矩陣，

③應(yīng)用Un（k）求解船體及液艙物面速度勢；

④判斷船體加速度的相對誤差e，如誤差足夠小，繼續(xù)第5）步計(jì)算，否則重復(fù)第4）步。

5）應(yīng)用上步最后迭代得到的 Un（u）和 U·n（u）步進(jìn)船體位移，tn＋1時(shí)刻船體位移以三階泰勒展開：

應(yīng)用最后計(jì)算的液艙內(nèi)部流體自由表面速度勢計(jì)算液艙內(nèi)部流體自由表面速度。步進(jìn)下一時(shí)刻tn＋1液艙內(nèi)部自由表面上?φ／?t、自由表面、瞬時(shí)濕表面及船體瞬時(shí)濕表面位置。

3 時(shí)域數(shù)值計(jì)算及試驗(yàn)驗(yàn)證

3.1 新阻尼模型的效率

能量耗散條件提出的主要目的是在液艙晃蕩與船體耦合運(yùn)動(dòng)的求解中應(yīng)用勢流理論對晃蕩產(chǎn)生的水動(dòng)力進(jìn)行符合實(shí)際的模擬。關(guān)鍵部分就是能量耗散系數(shù)μ的選取。為研究其在長時(shí)間數(shù)值模擬中的作用，對文獻(xiàn)［1］試驗(yàn)中帶液艙駁船橫蕩運(yùn)動(dòng)（如圖3）進(jìn)行模擬。為充分反映晃蕩特性，數(shù)值計(jì)算時(shí)長都選入射波周期的60倍，對μ取0.0、0.3、0.5、0.7、0.9、1.1。典型激勵(lì)頻率下液艙晃蕩與駁船橫蕩耦合運(yùn)動(dòng)時(shí)歷如圖4所示。對時(shí)歷曲線穩(wěn)態(tài)段做快速傅里葉變換得到運(yùn)動(dòng)RAO，通過與試驗(yàn)值的比較（如圖5），得出μ的選取規(guī)律。

圖3 帶液艙駁船橫蕩運(yùn)動(dòng)模型Fig.3 Box-shaped hull section

圖4 液艙晃蕩與駁船橫蕩耦合運(yùn)動(dòng)時(shí)歷曲線Fig.4 Time history of sloshing and sway coupling motion

從圖4中可看出，在不考慮液艙晃蕩能量耗散時(shí)（μ＝0），即使在模擬60個(gè)周期后，時(shí)歷曲線仍然難以達(dá)到穩(wěn)態(tài)并且呈現(xiàn)明顯的包絡(luò)特性。一方面是由于液艙晃蕩共振效應(yīng)導(dǎo)致的。另一方面是由于數(shù)值計(jì)算中未考慮液艙晃蕩能量耗散的影響［8］。實(shí)際上，在晃蕩過程中由于流體粘性及波浪破碎等因素，能量耗散是伴隨著整個(gè)流體晃蕩過程的，在經(jīng)歷多個(gè)周期之后，液艙內(nèi)部流體運(yùn)動(dòng)趨于穩(wěn)態(tài)。增大μ的取值不但改變橫蕩運(yùn)動(dòng)幅值且可更快速的使多頻非簡諧時(shí)歷曲線衰減成單頻的穩(wěn)態(tài)簡諧曲線，這說明本文能量耗散條件的作用效率較高。

圖5 μ對帶有液艙的浮體橫蕩運(yùn)動(dòng)RAO的影響Fig.5 Sway coupling motion RAO with different μ

從圖5中可看出，當(dāng)入射波頻率ω向晃蕩一階固有頻率ω1＝8.65 rad／s靠近時(shí)，液艙晃蕩減小了駁船的橫蕩運(yùn)動(dòng)。在大于ω1附近，液艙晃蕩明顯加劇了浮體的橫蕩運(yùn)動(dòng)，可以看出ω＝9.5 rad／s附近，帶有液艙的浮體運(yùn)動(dòng)RAO出現(xiàn)峰值，增大μ可以減輕晃蕩對船體橫蕩運(yùn)動(dòng)的影響。在不考慮液艙晃蕩能量耗散影響時(shí)，數(shù)值結(jié)果與試驗(yàn)結(jié)果在峰值附近差異較大，前者要明顯大于后者。因此能量耗散系數(shù)μ的選取對計(jì)算精度的影響很大。

帶液艙船體橫蕩運(yùn)動(dòng)RAO與不帶液艙船體橫蕩運(yùn)動(dòng)RAO曲線在兩曲線交點(diǎn)處發(fā)生明顯分歧（該工況時(shí)在ω＝9.0 rad／s處），該頻率定義為交界頻率ωd。當(dāng)波頻小于ωd時(shí)，液艙晃蕩減弱了船體橫蕩運(yùn)動(dòng)；當(dāng)波頻大于ωd時(shí)，液艙晃蕩反而加劇了船體橫蕩運(yùn)動(dòng)。交界頻率 ωd處對應(yīng)的橫蕩運(yùn)動(dòng)RAO幅值幾乎不受能量耗散系數(shù) μ的影響。當(dāng)ω ＜ωd時(shí)，μ按晃蕩固有頻率處誤差最小選取μ＝0.3與試驗(yàn)值吻合較好，當(dāng)ω＞ωd時(shí)μ按RAO曲線峰值處誤差最小選取μ＝0.9與試驗(yàn)值吻合較好?？砂l(fā)現(xiàn)，規(guī)則波中對液艙晃蕩與船體耦合運(yùn)動(dòng)進(jìn)行數(shù)值模擬時(shí)，μ的取值是有規(guī)律的，可能對于某些工況來說μ需進(jìn)一步修正，但暫可認(rèn)為該規(guī)律是最合適的。

3.2 LNG船的試驗(yàn)驗(yàn)證

為了驗(yàn)證上文所述規(guī)則波中能量耗散系數(shù)μ選取的規(guī)律，選取一個(gè)典型的具有試驗(yàn)數(shù)據(jù)的LNG運(yùn)輸船如表1，3種液艙裝載工況（R0，R2，R3）進(jìn)行研究如表2所示。由于只關(guān)注橫浪規(guī)則波下LNG船的橫蕩及橫搖運(yùn)動(dòng)，船艏及船尾形狀對運(yùn)動(dòng)影響較小，因此可用二維數(shù)值模型進(jìn)行計(jì)算。計(jì)算中LNG船按排水體積及形寬不變等效為吃水相同的駁船，液艙形狀按實(shí)際LNG船液艙給定。

表1 LNG船和液艙的特征尺度Table 1 Main dimensions of LNG carrier and its tanks

表2 計(jì)算參數(shù)Table 2 Numerical parameters used in the computation

將數(shù)值結(jié)果與文獻(xiàn)［20］試驗(yàn)結(jié)果進(jìn)行對比（如圖6）可看出，橫搖運(yùn)動(dòng)RAO與試驗(yàn)值吻合較好，驗(yàn)證了能量耗散系數(shù)μ選取規(guī)律的適用性。雖然現(xiàn)有的模型不能對自由表面做精確的模擬（未考慮波浪破碎、砰擊等作用），但并不影響對液艙晃蕩與船體耦合運(yùn)動(dòng)的計(jì)算，而這才是主要目的。

圖6 LNG船橫搖運(yùn)動(dòng)RAO的比較Fig.6 Comparison of roll motion RAO of an LNG carrier

3.3 不規(guī)則波中帶液艙船體的耦合運(yùn)動(dòng)分析

為考察實(shí)際海況中的液貨船耦合運(yùn)動(dòng)效應(yīng)是否顯著，研究時(shí)域結(jié)果特性及能量耗散條件的作用效果。計(jì)算了文獻(xiàn)［11-12］中帶液艙駁船（如圖7）在不規(guī)則波中的運(yùn)動(dòng)，并與試驗(yàn)及頻域結(jié)果進(jìn)行比較。選用第17界ITTC推薦的JONSWAP海浪譜，選取中等海情的短期海況，浪向?yàn)闄M浪，譜峰周期為1.6 s，有義波高為0.06 m。液艙裝載工況如表3所示，μ分別取0.0、0.3、0.5。

圖7 帶液艙駁船試驗(yàn)?zāi)Ｐ虵ig.7 Barge model in the basin

表3 計(jì)算工況（h為艙內(nèi)部水深，L是液艙寬度）Table 3 Calculated condition（h is water depth and L is width of tank）

影響耦合運(yùn)動(dòng)的能量耗散主要有2個(gè)方面：一是晃蕩中由于液體粘性及波浪破碎等引入的能量耗散，二是由駁船外流場波浪力引起的橫搖粘性阻尼。因此橫搖阻尼系數(shù)Cd可按無液艙的駁船進(jìn)行選取，文獻(xiàn)［11］中發(fā)現(xiàn)Cd＝0.2時(shí)較適合當(dāng)前海況。由于駁船較寬，加載液艙后會(huì)明顯影響駁船橫搖運(yùn)動(dòng)，對橫蕩運(yùn)動(dòng)影響相對較小，因此重點(diǎn)研究該駁船的橫搖運(yùn)動(dòng)。圖8、9是對艙晃蕩與駁船橫搖耦合運(yùn)動(dòng)時(shí)歷曲線做交叉譜分析后得到的運(yùn)動(dòng)RAO與試驗(yàn)及頻域結(jié)果的比較。

圖8 Cd和μ對帶液艙駁船橫搖RAO的影響，h／L＝0.237 5Fig.8 Roll RAO of barge with tanks in different Cdand μ，h／L＝0.237 5

圖9 μ對帶液艙駁船橫搖運(yùn)動(dòng)RAO的影響Fig.9 Roll RAO of barge with tanks in different μ

從圖8中看出，不規(guī)則波中液艙晃蕩在ω1附近的減搖效果明顯，卻在某些頻率下加劇了駁船運(yùn)動(dòng)，這與規(guī)則波問題的分析一致。不考慮μ時(shí)，改變Cd對橫搖運(yùn)動(dòng)RAO的第1個(gè)峰值影響較大，而對第2個(gè)峰值幾乎沒有影響。這是由于雖然在耦合運(yùn)動(dòng)中增加μ與Cd都可對橫搖運(yùn)動(dòng)起到減搖的作用，但作用機(jī)理不同，產(chǎn)生的效果也不同。第1個(gè)峰值主要是駁船自身慣性引起，而第2個(gè)峰值是由液艙晃蕩引起，Cd只能對前者產(chǎn)生作用。第2個(gè)峰值的能量耗散效應(yīng)更為明顯，因此在晃蕩與船體耦合運(yùn)動(dòng)中僅僅通過加載橫搖阻尼修正運(yùn)動(dòng)是不夠的，晃蕩能量耗散效應(yīng)必須考慮。

從圖9中可以看出，μ＝0.3，Cd＝0.2時(shí)較適合當(dāng)前海況下駁船的橫搖耦合運(yùn)動(dòng)，數(shù)值與試驗(yàn)結(jié)果吻合較好。μ的選擇與充液比例無關(guān)，這在時(shí)域模擬中是十分方便及重要的。在不規(guī)則波中μ的取值是單一的，不需像在規(guī)則波中那樣根據(jù)入射波波頻與ω1的關(guān)系選取。頻域結(jié)果與本文數(shù)值及文獻(xiàn)［11-12］中的試驗(yàn)結(jié)果在第2個(gè)峰值處有較大差別。這種差別在2種充液狀態(tài)下都有發(fā)生，其原因一是頻域方法對橫搖粘性和液艙晃蕩粘性修正方法的局限性導(dǎo)致，二是艙內(nèi)液體晃蕩的多頻振蕩，使得在船體運(yùn)動(dòng)中疊加非單色入射波頻的激勵(lì)力成分，從而導(dǎo)致船體運(yùn)動(dòng)并非滿足線性單頻假設(shè)，因而對艙液流動(dòng)采用線性頻域理論處理會(huì)產(chǎn)生較大誤差。

4 結(jié)論

本文基于時(shí)域勢流理論，建立了一種新的計(jì)算液艙晃蕩與船體非線性耦合運(yùn)動(dòng)的模型。新模型給出了一種改進(jìn)的液艙壁面邊界條件以近似晃蕩能量耗散效應(yīng)。通過與試驗(yàn)及頻域結(jié)果的比較發(fā)現(xiàn)：

1）新模型成功地近似了液艙晃蕩與船體耦合運(yùn)動(dòng)中液艙晃蕩導(dǎo)致的能量耗散效應(yīng)。只需經(jīng)歷短時(shí)間的數(shù)值模擬即可達(dá)到穩(wěn)態(tài)。這種方法較之其他的基于粘性流體方程的模型效率更高。

2）分析了能量耗散系數(shù)μ修正耦合運(yùn)動(dòng)的作用及限晃效率，得到規(guī)則波及不規(guī)則波中較優(yōu)化的能量耗散系數(shù)選取規(guī)律。其與液艙的濕表面形狀及充液比例無關(guān)，且不受模型尺度的影響，這在時(shí)域數(shù)值模擬中是非常方便及高效的。

3）不規(guī)則波中非線性理論得到的帶液艙船體橫搖運(yùn)動(dòng)RAO由晃蕩引起的第2個(gè)峰值要比線性頻域理論計(jì)算的大，這與試驗(yàn)得到了一致的結(jié)果。因而對不規(guī)則波中艙液流動(dòng)采用常規(guī)船體運(yùn)動(dòng)分析的線性頻域理論處理會(huì)產(chǎn)生較大誤差。

作為基礎(chǔ)理論研究，本文模擬限于二維問題，由于耦合分析中晃蕩總體載荷是關(guān)鍵因素，且船舶橫搖受晃蕩影響最顯著，因此反應(yīng)了最主要的影響因素，結(jié)論可指導(dǎo)進(jìn)一步的三維分析。

［1］FALTINSEN O M，TIMOKHA A N.Adaptive multimodal approach to nonlinear sloshing in a rectangular tank［J］.Journal of Fluid Mechanics，2001，432：167-200.

［2］余建星，唐建飛，劉源，等.浮式液化天然氣生產(chǎn)儲(chǔ)卸裝置重氣泄漏擴(kuò)散模擬分析［J］.天津大學(xué)學(xué)報(bào)，2013，46（5）：381-386.

YU Jianxing，TANG Jianfei，LIU Yuan，et al.Simulation and analysis of dispersion of heavy gas leakage from FLNG installation［J］.Journal of Tianjin University，2013，46（5）：381-386.

［3］洪亮，朱仁傳，繆國平，等.波浪中船體與液艙晃蕩耦合運(yùn)動(dòng)的時(shí)域數(shù)值計(jì)算［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2012，33（5）：635-641.

HONG Liang，ZHU Renzhuan，MIAO Guoping，et al.Numerical calculation of ship motions coupled with tank sloshing in time domain based on potential flow theory［J］.Journal of Harbin Engineering University，2012，33（5）：635-641.

［4］李裕龍，朱仁傳，繆國平，等.基于OpenFOAM的船舶與液艙流體晃蕩在波浪中時(shí)域耦合運(yùn)動(dòng)的數(shù)值模擬［J］.船舶力學(xué)，2012，16（7）：750-758.

LI Yulong，ZHU Renchuan，MIAO Guoping，et al.Simulation of ship motions coupled with tank sloshing in time domain based on Open FOAM［J］.Journal of Ship Mechanics，2012，16（7）：750-758.

［5］FALTINSEN O M，TIMOKHA A N.A multimodal method for liquid sloshing in a two-dimensional circular tank［J］.Journal of Fluid Mechanics，2010，665：457-479.

［6］BOOKI K，YUNG S.Coupled seakeeping with liquid sloshing in ship tanks［C］／／Proceedings of the ASME 27th International Conference on Offshore Mechanics and Arctic Engineering.Estoril，Portugal，2008：15-20.

［7］MALENICA ?，ZALAR M，CHEN X B.Dynamic coupling of seakeeping and sloshing［C］／／Proceedings of the 13th International Offshore and Polar Engineering Conferece.Honolulu，Hawaii，USA，2003：486-492.

［8］HUANG S，DUAN W Y，MA Q W.An approximation to energy dissipation in time domain simulation of sloshing waves based on linear potential theory［J］.China Ocean Engineering，2011，25（2）：189-200.

［9］DATTA R，SEN D.The simulation of ship motions using a b-spline-based panel method in time domain［J］.Journal of Ship Research，2007，51（3）：267-284.

［10］YAN S，MA Q W.Numerical simulation of fully nonlinear interaction between steep waves and 2D floating bodies using the QALE-FEM method［J］.Journal of Computational Physics，2007，221：666-692.

［11］MOLIN B，REMY F，LEDOUX A，et al.Effect of roof impacts on coupling between wave response and sloshing in tanks of LNG-carriers［C］／／Proceedings of the 27th International Conference on Offshore Mechanics and Arctic Engineering.Estoril，Portugal，2008：15-24.

［12］MOLIN B，REMY F，RIGAUD S，et al.Ch.LNG-FPSO’s：frequency domain coupled analysis of support and liquid cargo motions［C］／／Proceedings of INAM Conference.Rethymnon，Greece，2002.

［13］WU G X，MA Q W，EATOCK T R.Numerical simulation of sloshing waves in a 3D tank based on a finite element method［J］.Applied Ocean Research，1998，20：337-355.［14］CHEN X B.Hydrodynamics in offshore and naval applications part 1［C］／／Proceedings of the 6th International Conference on Hydrodynamics，University of Western Australia.Perth，Australia，2004.

［15］MA Q W.Numerical simulation of nonlinear interaction between structures and steep waves［D］.London：City University，1998.

［16］段文洋.浮體大幅運(yùn)動(dòng)非線性水動(dòng)力研究［D］.哈爾濱：哈爾濱工程大學(xué)，1995：67-68.

DUAN Wenyang.Nonlinear hydrodynamic forces acting on a ship undergoing large amplitude motions［D］.Harbin：Harbin Engineering University，1995：67-68.

［17］RANADEV D，DEBABRATA S.A B-spline based method for radiation and diffraction problem［J］.Ocean engineering，2006，33：2240-2259.

［18］MA Q W，YAN S.QALE-FEM for numerical modelling of non-linear interaction between 3D moored floating bodies and steep waves［J］.International Journal for Numerical Methods in Engineering，2009，78：713-756.

［19］GEAR W.Numerical initial value problems in ordinary differential equations［M］.Englewood Cliffs：Prentice-Hall，NJ，1971.

［20］MARIN.Seakeeping tests for the JIP SALT LNG carrier［R］.Report No.18127-1＆3-SMB，Jan.1st，2003.

Nonlinear time domain simulation of sloshing and coupled ship motion

HUANG Shuo1，2，3，DUAN Wenyang3，YOU Yage1，2，JIANG Jinhui4，WANG Wensheng1，2
（1.Key Laboratory of Renewable Energy，Chinese Academy of Sciences，Guangzhou 510000，China；2.Guangzhou institute of energy conversion，Chinese Academy of Sciences，Guangzhou 510000，China；3.College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China；4.Shanghai Shipping Scientific Research Institute，Shanghai 200135，China）

With the development of LNG carriers，this paper focuses on the study of the influencing characteristics of sloshing in the coupled nonlinear ship motions.The aim is to develop the time domain prediction techniques under both regular and irregular wave conditions.The energy dissipation was modeled by changing the boundary condition on tank solid boundaries.Under the linear free surface condition，to find out the potential，the fully nonlinear model was adopted for the inner tank flow and the boundary element model of Green function satisfying transient nonlinear body surface condition was used for external flow outside the tank.The boundary value problem was solved by the B-spline higher-order panel method.The ISITIMFB（iterative semi implicit time integration method for floating bodies）was applied to solve the body’s velocity and displacements.An extended principle to determine the dissipation coefficient μ is extracted and the effectiveness of this method was demonstrated from the comparison results of the experimental and numerical situation.

sloshing；energy dissipation；fully nonlinear model；time-domain coupled motion；LNG carrier；B-spline；ISITIMFB

10.3969／j.issn.1006-7043.201307076

O352；U661.1

A

1006-7043（2014）09-1045-08

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1006-7043.201307076.html

2013-07-31. 網(wǎng)絡(luò)出版時(shí)間：2014-08-26.

國家自然科學(xué)基金資助項(xiàng)目（70271029，41106031）；國家973計(jì)劃基金資助項(xiàng)目（2012CB723804）.

黃碩（1984-），女，助理研究員，博士；段文洋（1967-），男，長江學(xué)者，教授，博士生導(dǎo)師.

黃碩，E-mail：huangshuo＠m(xù)s.giec.ac.cn.