FPSO液艙晃蕩與船舶時域耦合運動數(shù)值模擬

操戈，李旭，張詠鷗，王康

1 海軍裝備部 駐武漢地區(qū)軍事代表局，湖北武漢430064

2 華中科技大學(xué)船舶與海洋工程學(xué)院，湖北武漢430074

3 華中科技大學(xué)船舶與海洋工程水動力學(xué)湖北省重點實驗室，湖北武漢430074

0 引 言

隨 著LNG（Liquefied Natural Gas），F(xiàn)PSO（Floating Production Storage and Offloading）和FLNG（Floating Liquefied Natural Gas）等液貨船的發(fā)展，液艙晃蕩對船舶運動的影響成為關(guān)注的熱點。此外，為滿足市場需求，船舶尺寸越來越大，船舶運動與液艙晃蕩的耦合影響也愈發(fā)重要。

與船舶耦合的液艙晃蕩在過去的研究中被簡化為自由液面或者質(zhì)點彈簧，這能精確求解小幅晃蕩。幾十年來，人們針對耦合運動進行了許多研究，主要分為2 類：基于線性勢流理論的頻域研究和基于非線性粘性流的時域研究。Molin 等［1］基于線性勢流理論的頻域方法研究了船舶與液艙晃蕩的耦合運動，研究假設(shè)船舶液艙的運動是線性的，但當(dāng)液艙晃蕩很激烈時，線性理論就有一定的局限性。船舶與液艙晃蕩的時域耦合運動是最近研究的熱點問題，脈沖響應(yīng)法（Impulse Response Function，IRF）被用于求解時域耦合問題。Cummins［2］利 用IRF法將船舶頻域結(jié)果轉(zhuǎn)化為時域結(jié)果進行了分析。Kim 等［3］將IRF 法擴展到求解船舶與液艙晃蕩耦合運動，并驗證得到仿真計算結(jié)果與試驗結(jié)果具有較好的一致性。李裕龍等［4］基于OpenFOAM 研究了船舶與液艙晃蕩在波浪中的時域耦合運動，數(shù)值模擬計算能夠清晰地顯現(xiàn)出液艙晃蕩對船體全局運動的影響，且船體運動計算結(jié)果與模型試驗結(jié)果吻合良好。洪亮等［5］對S175 加載方形液艙在迎浪、橫浪等不同工況下液艙流體晃蕩及其與船體運動時域耦合分別進行了計算模擬與驗證研究，耦合運動模擬結(jié)果能清晰地反映液艙晃蕩對船體運動的影響，且數(shù)值結(jié)果與試驗結(jié)果吻合良好，并具有較高的計算效率。黃碩等［6］對LNG 船和FPSO 液艙晃蕩與船舶耦合運動問題的研究思路及進展進行了綜述，分別對求解該類問題的模型試驗和數(shù)值計算方法的最新研究成果進行了總結(jié)，分析了模型試驗的特點和數(shù)值計算方法的應(yīng)用范圍，探討了試驗研究與數(shù)值模擬的關(guān)鍵技術(shù)及優(yōu)缺點，并提出了加強相關(guān)研究的建議。

試驗是研究船舶與液艙晃蕩耦合運動最有用的方法之一。模型試驗結(jié)果常用于證實理論和仿真計算結(jié)果的準(zhǔn)確性。為研究液艙晃蕩對船舶運動的影響，Molin，Clauss 以及李裕龍等［1，4，7］進行了一系列的相關(guān)試驗。

近年來，人們對由液艙晃蕩引起的載荷也做了很多研究，提出了許多數(shù)值和理論的方法。Loots 等［8］提出了改進的VOF（Volume of Fluid）法以計算液艙晃蕩的沖擊載荷，并與試驗數(shù)據(jù)進行了比較。Lee 等［9］利用CFD 軟件計算了LNG 的矩形液艙在外部激勵力作用下的響應(yīng)，并和試驗結(jié)果做了對比，研究表明液體粘性、液體與氣體的壓縮性，以及密度比對沖擊載荷的影響很大。張書誼等［10］采用CFD 軟件Fluent 模擬了二維矩形液艙在不同艙內(nèi)水深、不同激振頻率時的橫蕩，發(fā)現(xiàn)在距自由面較深點處，流體載荷的計算結(jié)果與試驗值相符合。金晶等［11］根據(jù)改進的VOF 法建立了晃蕩載荷計算方法，可用于對補給艦船的液艙晃蕩設(shè)計載荷進行快速、準(zhǔn)確的預(yù)報。但有關(guān)FPSO液艙晃蕩與船舶耦合運動時域及頻域計算方法的系統(tǒng)性分析還很少，液艙晃蕩對船舶運動的影響在過去的研究中經(jīng)常被忽略或線性化。因此，為準(zhǔn)確預(yù)報FPSO 運動，本文擬計算不同液艙液量下液體粘性對液艙晃蕩與船舶時域耦合運動的影響，并對時域及頻域計算方法的適用性進行系統(tǒng)性分析。

本文將通過求解時域方程，采用能考慮液艙內(nèi)流體粘性的CFD 軟件Fluent 及IRF 方法分析船舶與液艙晃蕩耦合運動，并將時域、頻域方程求解結(jié)果與試驗數(shù)據(jù)進行對比。一方面，通過CFD 計算液艙晃蕩所產(chǎn)生的力及力矩，并以外力形式作用于船體；另一方面，將船舶運動仿真計算結(jié)果代入CFD 以求解液艙晃蕩運動。通過兩方面步驟的迭代計算，即可求解出船舶與液艙晃蕩的耦合運動。隨后，將試驗結(jié)果與仿真結(jié)果進行對比以得到阻尼比的范圍。最后，通過時域和頻域方法求出不同液艙液量下船舶幅值響應(yīng)算子（RAO）。

1 船舶與液艙晃蕩耦合運動

1.1 船體運動與液艙晃蕩固有頻率

勢流理論是計算船舶在誘導(dǎo)波中運動的有效方法。勢流理論假設(shè)流體無粘無旋，那么速度勢可用來表示整個流域，如式（1）所示：

式中，Φ 為流體速度勢。由于船體表面的復(fù)雜性，直接求解整個流域的速度勢較困難。速度勢可以分解為入射勢、繞射勢和輻射勢，可根據(jù)繞射理論和輻射理論分別予以求解。已知總的速度勢就可求得到波作用于船舶的力。船體可視為剛性體，船體運動方程可寫為

式中：M 為船體質(zhì)量矩陣；C 為船體阻尼系數(shù)矩陣；K 為回復(fù)力剛度矩陣；Fw(t)為船舶受到的力；ζ 為船舶位移。

具有一定水深液艙的固有頻率可以從波的色散關(guān)系中得出

式中：h 為水深；k 為波數(shù)；g 為重力加速度。液艙寬度和長度分別表征橫向和縱向振型自由液面寬度，推導(dǎo)可得液艙每階模態(tài)的固有頻率

式中：ωn為第n 階模態(tài)的固有頻率；B 為液艙自由液面寬度；n 為頻率階次。

1.2 頻域耦合運動

液艙晃蕩會影響到船舶的運動，晃蕩引起的力作為外部激勵將作用于船體。受晃蕩載荷作用的船舶6 個自由度的運動方程可寫為

式中：Ma(ω)為船舶附加質(zhì)量；Fs(ω)為液艙對船舶的作用力。質(zhì)量矩陣M 包括液體的質(zhì)量。忽略液艙內(nèi)液體的慣性力，晃蕩引起的載荷可寫為

式中：Mas(ω)為液艙的附加質(zhì)量；Cs(ω)為液艙的阻尼系數(shù)；Ks為液艙修正的回復(fù)力剛度，可以寫為

式中：Is為液艙自由面對x 軸的二階慣性矩；ρs為液艙內(nèi)液體的密度。

橫搖運動時粘性的作用不能忽略，可采取添加附加線性阻尼的辦法來填補粘性的影響。

將式（6）和式（8）代入式（5），頻域的耦合運動方程可改寫為如下形式：

船舶運動位移和波浪力可以寫成如下形式：

將式（10）代入式（9），可以得到頻域耦合運動方程的另一種形式：

1.3 IRF 法

和附加質(zhì)量、阻尼系數(shù)一樣，所有的水動力系數(shù)均是在頻域中求得，相應(yīng)的力應(yīng)轉(zhuǎn)化到時域。利用IRF 法，經(jīng)過推導(dǎo)，可得一階波浪輻射力的表達式

在實際的數(shù)值計算中，式（14）積分上限不可能是無窮大，存在數(shù)值計算截斷誤差，為減小這種截斷誤差，采用Kim 等［3］提出的方法，將改寫為式（15）的形式，并對截斷誤差做相應(yīng)的變換。

圖1 給出了截止頻率和截斷誤差對垂蕩方向（Heave-Heave）遲滯方程數(shù)值的影響。圖中：橫坐標(biāo)t 表示時間；縱坐標(biāo)R33表示垂蕩方向遲滯方程數(shù)值；ωmax為截止頻率；L 為船舶長度。在時間較早時，遲滯方程的值較大，截斷誤差的影響較小，而到后期，截斷誤差的影響較大，可見有必要考慮截斷誤差的影響。

1.4 時域耦合運動

當(dāng)考慮液艙對船舶運動影響時，船舶運動方程可寫為

式中：Fext(t)表示波浪對船體表面的作用力；Fsa(t)表示液艙內(nèi)液體晃蕩引起的對液艙壁的作用力。將液艙內(nèi)液體的慣性力從Fsa(t)分離出來，則Fsa(t)可改寫為如下形式：

式中：Ms為液艙內(nèi)液體的質(zhì)量矩陣；為液艙內(nèi)液體對液艙壁的靜水力、水動力以及粘性力之和。

綜合以上推導(dǎo)，船舶與液艙晃蕩時域耦合運動方程可寫為

圖1 截止頻率和截斷誤差對垂蕩方向遲滯函數(shù)的影響Fig.1 Effects of cut off frequency and truncation error on heave-heave retardation function

其中，

式中：Vliquid為液艙內(nèi)液體體積；LGg為Gm與gs間的垂向距離，其中Gm為原始船舶的重力，gs為液艙內(nèi)液體的重力。

2 計算方法及模型

2.1 數(shù)值計算方法

本文求解了船舶與液艙晃蕩運動的頻域方程和時域方程，求解這2 種方程涉及的程序以及相應(yīng)的流程如圖2 和圖3 所示。

圖2 船舶與液艙晃蕩頻域耦合運動計算流程圖Fig.2 Flow chart of solving ship motion and tank sloshing coupling in frequency domain

圖3 船舶與液艙晃蕩時域耦合運動計算流程圖Fig.3 Flow chart of solving ship motion and tank sloshing coupling in time domain

在頻域中，利用法國船級社（BV）的三維水動力學(xué)軟件Hydrostar 計算船舶以及液艙的水動力學(xué)系數(shù)（附加質(zhì)量、阻尼系數(shù)等），通過求解船舶與液艙的頻域耦合運動方程，得到船舶運動RAO。在時域中，利用Hydrostar 計算船舶的水動力系數(shù)以及波浪載荷，由IRF 法得到時域波浪輻射力?；赩OF 法，采用CFD 軟件Fluent 模擬液艙晃蕩，在Fluent 的UDF 中迭代求解船舶與液艙晃蕩耦合運動的時域方程。數(shù)值求解得到的液艙晃蕩力以及力矩作為外部激勵作用于船體，船舶運動仿真計算結(jié)果逐步代入求解液艙晃蕩運動。

2.2 FPSO 模型

引用MARINE 實驗室完成的FPSO 試驗作為算例求取船舶與液艙晃蕩的耦合運動。圖4 為船舶及液艙的幾何外形。FPSO 的長、寬、吃水分別為285.0，63.0 和13.0 m，排水量為220 017.6 m3。圖5 所示為FPSO 在Hydrostar 中的網(wǎng)格劃分圖，有5 024 個單元，圖中液艙液量為56%。試驗探討了2 個液艙液量均為18%，37%和56%這3 種工況，后部液艙和前部液艙的長、寬、高分別為49.68，46.92，32.23 m 以及56.62，46.92，32.23 m，液艙底部至龍骨的距離為3.3 m。

圖4 FPSO 外形與液艙分布圖Fig.4 Sketch of FPSO with internal tanks

圖5 FPSO 與液量為56%液艙網(wǎng)格劃分圖Fig.5 Mesh of FPSO and tanks（56% filling level）

3 結(jié)果及分析

3.1 固有頻率及附加質(zhì)量和阻尼系數(shù)

表1 給出了液艙液量為18%，37%和56%時橫向與縱向振型的固有頻率。除長度以外，2 個液艙的幾何尺寸一樣，因此，2 個液艙的橫向固有頻率相等，而結(jié)合式（4）可知，縱向固有頻率不同。圖6 和圖7 分別為2 個液艙的橫搖、縱搖附加質(zhì)量和阻尼系數(shù)。圖中橫坐標(biāo)ω表示頻率，縱坐標(biāo)分別表示橫搖、縱搖附加質(zhì)量和阻尼系數(shù)。

表1 液艙固有頻率Tab.1 Natural frequencies of tanks

圖6 前部液艙的橫搖、縱搖附加質(zhì)量和阻尼系數(shù)Fig.6 Roll and pitch added mass，damping coefficients of the front tank

圖7 后部液艙的橫搖、縱搖附加質(zhì)量和阻尼系數(shù)Fig.7 Roll and pitch added mass，damping coefficients of the rear tank

3.2 阻尼比分析

橫搖運動是最受關(guān)注的問題，當(dāng)考慮液艙影響時，橫搖運動變得更加復(fù)雜，求解式（11）可以得出船舶與液艙晃蕩耦合運動幅值響應(yīng)算子（RAO）。圖8 給出了3 種液艙液量下船舶橫搖運動RAO 試驗值與頻域計算結(jié)果的對比。圖中：橫坐標(biāo)ω 表示頻率；縱坐標(biāo)表示橫搖RAO；帶圓圈的曲線表示MARINE 實驗室完成的FPSO 試驗結(jié)果。

圖8 不同液量下橫搖運動RAO 試驗值與計算值對比（橫浪）Fig.8 Comparison of roll RAO of experimental results and computational results at different filling levels（beam sea）

從圖中可以看出，船舶液艙液量為0%時，一階固有頻率在0.5 rad/s 附近，但是一階固有頻率附近的計算結(jié)果和試驗結(jié)果相比偏大，為了能使計算結(jié)果更加準(zhǔn)確，加入了線性阻尼項，如式（8）所示。忽略液艙液量，系統(tǒng)阻尼比γ =0.1～0.15 時能夠得到較好的頻域結(jié)果。當(dāng)液艙液量為18%時，液艙橫搖一階固有頻率為0.49 rad/s，與船舶的一階固有頻率很接近，使得在0.42 rad/s 處RAO 波峰值變?yōu)椴ü戎担ǚ逯灯频搅?.6 rad/s 處。當(dāng)液艙液量為37%和56%時，液艙的一階固有頻率分別為0.66 rad/s 和0.74 rad/s，遠(yuǎn)高于船舶的一階固有頻率，其第2 個峰值分別在0.76 rad/s 和0.8 rad/s 附近。隨著液量的增加，第1 個峰值增加，第2 個峰值不斷減小。在固有頻率附近，液艙可起到減搖水艙的作用。由圖中可以看到，當(dāng)液艙液量較低時，頻域計算結(jié)果與試驗計算結(jié)果相差較大，這是由于粘性的影響很大，盡管加入了液艙的阻尼系數(shù)來減小忽略粘性帶來的影響，但還是存在一定的誤差；當(dāng)液艙液量較高時，慣性力占主導(dǎo)作用。

3.3 耦合結(jié)果

圖9 給出了3 種液艙液量下船舶縱搖運動的RAO。從圖中可以看出，計算結(jié)果與試驗結(jié)果吻合較好，不同液艙液量之間差異很小，這是由于縱搖的縱向慣性矩很大，遠(yuǎn)大于液艙晃蕩引起的波動影響，液艙的加入對船舶縱搖的影響較小。

圖9 不同液量下縱搖運動RAO 實驗值與計算值對比（迎浪）Fig.9 Comparison of pitch RAO of experimental results and computational results at different filling levels（head sea）

液艙內(nèi)液體粘性對船舶橫搖的影響顯著，因此有必要計算船舶與液艙晃蕩時域耦合運動。船舶運動的RAO 可以通過傅立葉變換來與前面的計算結(jié)果相比較。

圖10 給出了3 種液艙液量下船舶橫搖運動RAO 時域、頻域計算結(jié)果和試驗數(shù)據(jù)三者的對比?？梢钥闯?，頻域計算對與某些頻域點對應(yīng)的RAO 值預(yù)測過高，而時域計算則考慮了液艙內(nèi)液體的非線性影響，起到了很好的削峰作用；當(dāng)液艙液量為18%時，液艙的固有頻率為0.49 rad/s，與船舶固有頻率很相近, 其峰值的位置從0.5 rad/s 移動到0.6 rad/s，但峰值幾乎沒變。當(dāng)液艙液量為37%和56%時，結(jié)合圖8 能夠觀察到明顯的2 個峰值的分離。當(dāng)液艙的固有頻率遠(yuǎn)離船舶固有頻率時，液艙能起到減搖水艙的作用，并且隨著兩種頻率相差越大，減搖效果更明顯。當(dāng)液艙液量較低時，由于考慮了液艙內(nèi)流體粘性以及自由液面的非線性影響，RAO 的時域耦合計算結(jié)果和頻域耦合計算結(jié)果相比與試驗結(jié)果更接近。當(dāng)液艙液量較高時，艙內(nèi)流體慣性的影響相對粘性影響占主導(dǎo)作用，2 種方法的求解準(zhǔn)確度差異不大。

圖10 不同液量下橫搖運動RAO 實驗結(jié)果、時域和頻域計算結(jié)果對比（橫浪）Fig.10 Comparison of roll RAOs of experimental and computational results（frequency domain and time domain）at different filling levels（beam sea）

4 結(jié) 語

船舶在波浪中運動會激勵液艙內(nèi)液體產(chǎn)生晃蕩，而液艙晃蕩力又會對船舶運動姿態(tài)產(chǎn)生影響，相互間產(chǎn)生耦合作用。本文建立了在波浪載荷下考慮液艙內(nèi)流體粘性阻尼時船體與液艙晃蕩的時域耦合運動數(shù)值計算方法。首先，基于頻域勢流理論得到水動力系數(shù)以及波浪載荷，基于IRF 法在時域積分得到波浪輻射力。隨后，一方面通過CFD 計算獲取液艙晃蕩所產(chǎn)生的力及力矩加載于船體，另一方面，在船舶運動仿真計算結(jié)果的基礎(chǔ)上計算液艙晃蕩運動。最后，對FPSO 船型在不同液量下的時域耦合運動進行計算，并與頻域和試驗結(jié)果進行了對比分析。

對采用IRF 法求解時域波浪輻射力過程中截斷誤差影響的研究表明，在時間較早時，遲滯方程數(shù)值較大，截斷誤差的影響較小，到后期，截斷誤差的影響逐漸增大。此外，與試驗數(shù)據(jù)的對比研究表明，在不同液艙液量條件下，系統(tǒng)阻尼比γ =0.1～0.15 時能夠得到較好的計算結(jié)果。在耦合運動方面，液艙內(nèi)液體的粘性對耦合作用有較大影響：低液量時，粘性的影響相對慣性影響占主導(dǎo)作用，橫搖運動RAO 的時域耦合計算結(jié)果和頻域耦合結(jié)果相比與試驗值更接近；高液量時，慣性影響較粘性影響更加重要，時域和頻域這2 種計算方法的求解精度差異不大。本文中計算方法可作為FPSO 或相關(guān)船型液艙晃蕩和船舶耦合運動模擬的參考。

［1］MOLIN B，REMY F，LEDOUX A，et al. Effect of roof impacts on coupling between wave response and slosh?ing in tanks of LNG-carriers［C］//ASME 2008 27th In?ternational Conference on Offshore Mechanics and Arc?tic Engineering.Estoril，Portugal，2008：15-24.

［2］CUMMINS W E. The impulse response function and ship motions［R］. The Symposium on Ship Theory，Hamburg，Germany，1962.

［3］KIM Y，NAM B W，KIM D W，et al. Study on cou?pling effects of ship motion and sloshing［J］. Ocean En?gineering，2007，34（16）：2176-2187.

［4］李裕龍，朱仁傳，繆國平，等.基于OpenFOAM 的船舶與液艙流體晃蕩在波浪中時域耦合運動的數(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 OpenFOAM［J］. Journal of Ship Mechanics，2012，16（7）：750-758.

［5］洪亮，朱仁傳，繆國平，等.波浪中船體與液艙晃蕩耦合運動的時域數(shù)值計算［J］. 哈爾濱工程大學(xué)學(xué)報，2012，33（5）：635-647.HONG Liang，ZHU Renchuan，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-647.

［6］黃碩，段文洋，馬慶位，等.液艙晃蕩及其與船舶耦合運動問題的研究進展［J］. 船舶力學(xué)，2013，17（10）：1209-1220.HUANG Shuo，DUAN Wenyang，MA Qingwei，et al.Research progress on sloshing and coupling motion of sloshing and ship［J］. Journal of Ship Mechanics，2013，17（10）：1209-1220.

［7］CLAUSS G F，TESTA D，SPRENGER F. Coupling ef?fects between tank sloshing and motions of a LNG carri?er［C］// Proceedings of the ASME 2010 29th Interna?tional Conference on Ocean，Offshore and Arctic Engi?neering.Shanghai，China，2010.

［8］LOOTS E，PASTOOR W，BUCHNER B，et al. The numerical simulation of LNG sloshing with an im?proved volume of fluid method［C］// Proceedings of OMAE04 23rd International Conference on Offshore Mechanics and Arctic Engineering. Vancouver，British Columbia，Canada，2004.

［9］LEE D H，KIM M H，KWON S H，et al. A parametric sensitivity study on LNG tank sloshing loads by numeri?cal simulations［J］.Ocean Engineering，2007，34（1）：3-9.

［10］張書誼，段文洋. 矩形液艙橫蕩流體載荷的Fluent 數(shù)值模擬［J］. 中國艦船研究，2011，6（5）：73-77.ZHANG Shuyi，DUAN Wenyang. Numerical simula?tion of sloshing loads on rectangular tank based on Fluent［J］. Chinese Journal of Ship Research，2011，6（5）：73-77.

［11］金晶，薛鴻祥，唐文勇，等.大型補給艦船液艙晃蕩載荷計算方法［J］.中國艦船研究，2012，7（6）：50-56.JIN Jing，XUE Hongxiang，TANG Wenyong，et al.Calculation procedure of sloshing loads for large-scale depot ships［J］. Chinese Journal of Ship Research，2012，7（6）：50-56.