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

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

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

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

1.2 頻域耦合運動

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

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

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

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

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

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

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

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

1.3 IRF 法

和附加質量、阻尼系數一樣，所有的水動力系數均是在頻域中求得，相應的力應轉化到時域。利用IRF 法，經過推導，可得一階波浪輻射力的表達式

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

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

1.4 時域耦合運動

當考慮液艙對船舶運動影響時，船舶運動方程可寫為

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

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

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

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

其中，

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

2 計算方法及模型

2.1 數值計算方法

本文求解了船舶與液艙晃蕩運動的頻域方程和時域方程，求解這2 種方程涉及的程序以及相應的流程如圖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）的三維水動力學軟件Hydrostar 計算船舶以及液艙的水動力學系數（附加質量、阻尼系數等），通過求解船舶與液艙的頻域耦合運動方程，得到船舶運動RAO。在時域中，利用Hydrostar 計算船舶的水動力系數以及波浪載荷，由IRF 法得到時域波浪輻射力。基于VOF 法，采用CFD 軟件Fluent 模擬液艙晃蕩，在Fluent 的UDF 中迭代求解船舶與液艙晃蕩耦合運動的時域方程。數值求解得到的液艙晃蕩力以及力矩作為外部激勵作用于船體，船舶運動仿真計算結果逐步代入求解液艙晃蕩運動。

2.2 FPSO 模型

引用MARINE 實驗室完成的FPSO 試驗作為算例求取船舶與液艙晃蕩的耦合運動。圖4 為船舶及液艙的幾何外形。FPSO 的長、寬、吃水分別為285.0，63.0 和13.0 m，排水量為220 017.6 m3。圖5 所示為FPSO 在Hydrostar 中的網格劃分圖，有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%液艙網格劃分圖Fig.5 Mesh of FPSO and tanks（56% filling level）

3 結果及分析

3.1 固有頻率及附加質量和阻尼系數

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

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

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

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

3.2 阻尼比分析

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

圖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 附近，但是一階固有頻率附近的計算結果和試驗結果相比偏大，為了能使計算結果更加準確，加入了線性阻尼項，如式（8）所示。忽略液艙液量，系統阻尼比γ =0.1～0.15 時能夠得到較好的頻域結果。當液艙液量為18%時，液艙橫搖一階固有頻率為0.49 rad/s，與船舶的一階固有頻率很接近，使得在0.42 rad/s 處RAO 波峰值變為波谷值，波峰值偏移到了0.6 rad/s 處。當液艙液量為37%和56%時，液艙的一階固有頻率分別為0.66 rad/s 和0.74 rad/s，遠高于船舶的一階固有頻率，其第2 個峰值分別在0.76 rad/s 和0.8 rad/s 附近。隨著液量的增加，第1 個峰值增加，第2 個峰值不斷減小。在固有頻率附近，液艙可起到減搖水艙的作用。由圖中可以看到，當液艙液量較低時，頻域計算結果與試驗計算結果相差較大，這是由于粘性的影響很大，盡管加入了液艙的阻尼系數來減小忽略粘性帶來的影響，但還是存在一定的誤差；當液艙液量較高時，慣性力占主導作用。

3.3 耦合結果

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

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

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

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

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

4 結 語

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

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

［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 的船舶與液艙流體晃蕩在波浪中時域耦合運動的數值模擬［J］.船舶力學，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］洪亮，朱仁傳，繆國平，等.波浪中船體與液艙晃蕩耦合運動的時域數值計算［J］. 哈爾濱工程大學學報，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］. 船舶力學，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 數值模擬［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.