基于VOF法的液艙晃蕩數值模擬及載荷計算

劉楨兵

(斯特拉斯克萊德大學，格拉斯哥G40LZ)

舶液艙晃蕩作為一種常見的流體運動現象，通常發生在部分充滿液體的液艙中［1］。當液艙與船體運動發生耦合時，將使船體運動加劇甚至傾覆，同時液體晃蕩是非常復雜的流體運動現象，呈現出強的非線性和隨機性。

1 VOF法

流體體積法最早出現在DeBar(1974)的文章中，Hirt＆Nichols［2］形成了它完整的理論體系和實現方法，現在的VOF方法都是在Hirt＆Nichols的基礎上發展起來的。VOF法的基本原理是通過研究網格單元中流體和網格體積比函數來確定自由面，追蹤流體的變化，而非追蹤自由液面上質點的運動。VOF法使得流體的質量得到守恒，并可跟蹤自由面上發生翻轉、吞并、飛濺等復雜的自由面現象。VOF法簡單而有效，得到了廣泛使用，在實際工程中得到了令人滿意的計算結果。

2 建模及網格類型

本文所有的建模和網格劃分都是在Gambit軟件上完成的。Gambit可以生成結構化網格和非結構化網格，雖然非結構化網格對復雜幾何形狀的物理模型更有吸引力，但由于其在處理邊界層還存在著不足，而且在同等精度下非結構網格的網格數是結構化網格的好幾倍。由于使用全非結構化網格來處理時對計算硬件計算能力要求較高，本文的計算模型采用結構化網格。另外網格數量與時間步長的選擇對數值計算是有影響的。網格數量太多時要求時間步長很小，這樣會大大增加計算花費的時間;網格數量太少時計算花費的時間縮短許多，但是求解所得的結果不夠精確。因此網格數量應當適中，一方面可以節省計算花費的時間，讓計算可以有效率地運行;另一方面所得的計算結果能夠達到理論研究所需要的精度。本文對離散的RANS方程都是采用分離求解器求解(segregated solve)，壓力速度耦合以PISO法計算;本文所有的計算模型都帶有自由液面，壓力項離散使用體積力加權方法(body force weighted)，用gauss-seidel方法迭代求解，以代數多重網格技術加速迭代收斂速度。

3 程序流程及數值計算模型

程序流程及數值計算模型見圖1、2。

4 計算實例

4.1 算例1:低液深大幅晃蕩下的非線性和隨機性

在滿載狀況下的晃蕩固然要引起我們的重視，但是在空載出港或到港情況下，往往液艙內液體深度很低，并且考慮到壓載艙內不同水位的壓載水的晃蕩影響，低液深的大幅晃蕩對于船舶的穩性和操縱性會帶來很大程度的影響，這種現象不容忽視。

4.1.1 計算模型及問題描述

液艙模型如圖3所示。液艙做強迫橫蕩運動，橫向激勵為x=Asin(ωt)，其中橫蕩周期為1.3 s，頻率為0.77 Hz，橫蕩幅值為0.06 m，P1為壓力監測點。表1給出了幾種計算工況。

圖3 幾何模型(單位:m)

表1 幾種工況

4.1.2 計算結果及分析

1)圖4為自由液面在不同時刻的運動情況，液面變化可以說明以下幾點:

①自由液面發生沖頂的速度極為迅速，當晃蕩比較充分時波的破碎、翻卷現象比較明顯，液體晃動比較激烈。

②晃動比較充分時，可以看出雖然液艙做周期性運動，但自由表面波的變化卻并不是完全符合周期性的變化，表現了明顯的非線性特征。

③大幅晃蕩時自由表面的運動存在太多的隨機性，因而壓力的測量可能存在很多的不確定性，影響因素較多。

④晃蕩載荷還具有離散性和飽和性［3］。離散性說明拍擊力不呈現任何規律性，需要用概率分布來描述:飽和性則表示拍擊力的大小并不與外界激勵成正比關系，存在飽和現象，在非規則的激勵下表現得更為明顯。

圖4 自由表面的變化

2)圖5～6為不同工況的數值計算結果。

可以看出晃蕩壓力曲線表現了明顯的非線性特征，此外各工況之間的壓力曲線在不同程度上存在一定的差異，分析可以得出以下結論:

①由于網格的疏密以及時間步長的不同，不同工況的計算結果均會出現一些壓力突變值，分析原因可能是時間步長相對網格間距太大或太小而產生的截斷誤差和舍入誤差。

②由于晃蕩發生沖頂時，自由表面發生的波的破碎、翻卷等現象是短暫而迅速的，進而對壁面所造成的沖擊也是快速的。因此，不同的時間步長下所捕捉的壁面壓力的大小可能會有所不同。

③由于在監測壁面壓力時采用的是“貼片法”，即在壁面監測點選取一單元作為代替該點進行監測，所得出的壓力是該單元面的平均壓力值。因此在發生壁面沖擊時，單元劃分的大小對壓力的捕捉也有一定影響。

④當晃蕩為微幅晃蕩時，數值計算的結果受網格的疏密、時間步長的變化等影響不大，而當晃蕩比較激烈(大幅晃蕩，運動周期靠近共振周期時，液體晃蕩更加劇烈［4］)，晃動液體的沖頂現象比較嚴重，進而自由表面的波的破碎、翻卷等現象也尤為明顯。

4.2 算例2:高液深微幅晃蕩下同一主尺度下的矩形液艙與菱形液艙晃蕩對比

縱觀各類液貨船，出于不同用途以工程實際應用，液艙形式從單一的矩形衍變成菱形、圓形等，尤其是隨著LNG船舶的日益盛行，菱形艙室逐漸成為晃蕩研究的一個重點。

4.2.1 計算模型及問題描述

根據文獻［5］給定液艙模型，如圖7，其中矩形和菱形的主尺度是相同的，所不同的是菱形在頂部和底部的折角有所改變。液艙受強迫橫蕩運動，橫向激勵為x=Asin(ωt)，其中橫蕩周期為1.404 s，頻率為0.71 Hz，該液深時的諧振頻率為0.69 Hz，橫蕩幅值為0.015 m。表2給出艙室結構劃分情況，其中菱形1的折角處尺寸為0.1 m×0.1 m，菱形2的折角處尺寸為0.1 m×0.2 m。

圖8 矩形液艙與菱形液艙模型(單位:m)

表2 艙室結構劃分情況

在以下的壓力圖中分別以rhombus和Rhombus命名菱形1和菱形2，以示區別。

4.2.2 計算結果及分析

1)圖9為10 s與20 s時刻自由表面的位置。在圖中可以看出，在相同的外界激勵下，菱形液艙在晃蕩過程中自由液面擺動幅度比較大，矩形液艙則相對緩和一點。

2)圖10則給出了壁面點P1、P2的壓力的時間歷程。

圖9 不同時刻自由表面的形狀

圖10 3個模型在P1和P2點壁面壓力的時間歷程

3)圖11給出了2個菱形液艙在壁面點P1和P2處的壓力對比

圖11 菱形1與菱形2在P1和P2處壁面壓力的比較

比較圖中數據可以得出:

①矩形液艙的壁面所受壓力介于2種菱形液艙所受壓力值之間，說明2種液艙形式的孰優孰劣視具體情況而定。

②菱形液艙壁面壓力的時間歷程曲線存在很多壓力突變點，而矩形液艙壁面壓力的時間歷程曲線則表現了很好的周期性。

以上分析表明，采用矩形液艙與菱形液艙的防晃性能優劣視具體情況而定，但在一定折角尺寸下，菱形艙室抑制晃蕩作用更明顯。在某個尺寸范圍外，菱形艙室的折角處面板相當于一塊擋板會對運動的來流施加一個反沖擊力，故在各點處的壓力值都有所增大;而當進入某個尺寸范圍內，折角處則起了抑制晃蕩的作用。在LNG船液艙形式上之所以選擇菱形艙室考慮的因素不僅僅是壓力，在某個壓力差額范圍內，菱形艙室可能更有其他方面的優勢，比如折角處較矩形的直角更適合船體結構，或者從避免應力集中上更勝一籌，4個折角處亦可設為底邊艙、頂邊艙或者壓載水艙。菱形艙室作為一種新型艙室，具有很大的發展前途，須不斷改善。

5 結束語

分別從小液深大幅晃蕩下的非線性和隨機性驗證以及大液深微幅晃蕩下2種液艙形式優劣比較展開討論，應用matlab強大的圖表生成功能以及tecplot的自由液面繪制功能進行數據比較討論，結果清晰準確。

［1］ZHU Ren-qing，WU You-sheng，Incecic Atilla.Numerical simulation of liquid sloshing［J］.Shipbuilding of China，2004(3):21－22.

［2］Hirt C W，Nichols B D.Volume of Fluid(VOF)method for the dynamics of free boundaries［J］.J Comput.Phys，1981，39:201－225.

［3］ 劉澤民.關于船艙中液體晃蕩問題［J］.哈爾濱船舶工程學院學報，1992，13(3):241－249.

［4］ 祁江濤.基于VOF法的液艙晃蕩數值模擬［D］.無錫:702研究所，2007:32－33.

［5］ 侯玲.液艙晃蕩與彈性防晃結構的相互耦合作用研究［D］.鎮江:江蘇科技大學，2009:47－48.