袁心怡， 蘇 焱， 劉祖源

(武漢理工大學 交通學院， 武漢 430063)

晃蕩可以描述為部分填充的容器中液體的運動.液化天然氣(LNG)船、浮式液化天然氣裝置(FLNG)等采用液艙形式的液貨船在部分裝載的情況下，即使受到小幅激勵，也有可能發(fā)生劇烈的晃蕩現象，對液艙結構產生不利影響，威脅船舶航行安全.

研究晃蕩問題多采用數值方法，基于黏性流理論在流域內求解Navier-Stokes方程或基于勢流理論在流域內求解Laplace方程實現數值模擬.Wang等[1]基于比例邊界有限元法開發(fā)了一種研究擋板對艙內晃蕩影響的新方法，經過收斂性和研究，發(fā)現該方法具有計算簡單、精度高、速度快等優(yōu)點.Zhao等[2]基于勢流理論建立了LNG液艙在強制激勵下的非線性晃蕩數值模型，并加入了人工阻尼模型考慮耗能，研究了部分充液的矩形容器在平移和旋轉激勵下的晃蕩特性.Hu等[3]利用邊界元法求解控制方程，采用線性自由表面條件研究了多種形狀帶擋板的二維液艙在小幅激勵下晃蕩的固有頻率和振型，并提出了帶擋板的矩形液艙內晃蕩的固有頻率計算經驗公式.ünal等[4]研究了帶T形擋板的液艙內的晃蕩現象，考察了擋板高度、充液深度等參數對于艙壁水動力載荷和自由表面高程的影響，計算采用了兩種不同的方法，分別使用黏性的層流和紊流求解器.李金龍等[5]對一種新的幾何有限體積法isoAdvector進行修正，使之可應用于液艙晃蕩的模擬中，并基于不同的有限體積法對各種工況進行了模擬，發(fā)現修正后的方法所獲得的自由面位置和整體水動力載荷載荷精度更高，且能較好地模擬波面的翻卷和破碎.Fu等[6]提出了一種基于局部徑向基函數和二階顯式龍格庫塔法的半拉格朗日無網格模型來分析二維艙內的晃蕩現象，通過文獻對比證明了該方法的有效性和準確性.

對于強非線性的晃蕩問題，粒子法可以很好地捕捉自由面的大變形.但對于低充液率的容器，穩(wěn)定和精確地模擬容器內長時間的晃蕩具有難度.Green等[7]提出了一種有效的光滑粒子流體動力學(SPH)方法來克服這些困難，使邊界條件、耗散項和算法的誤差最小化，與試驗數據比較發(fā)現，數值模擬的結果非常精確.

布西內斯克(Boussinesq)方程是一類描述波浪傳播變形的數學模型.Bingham等[8]建立了數值離散方法，主要用于近岸變化底部處波浪場精確模擬.Antuono等[9]使用Boussinesq型平均水深方程建立了一個二維模型用于分析淺水中液體晃蕩現象.Su等[10]采用高精度Boussinesq型速度勢方程，對二維矩形艙內的晃蕩現象進行了數值模擬，結果表明基于Boussinesq型方程的數值模型在預測小激勵幅值晃蕩方面具有較好的性能.

本研究旨在探討淺水、微幅激勵條件下三維液艙內的晃蕩現象.假定液艙內為無黏無旋不可壓流體，流域內存在速度勢.將總速度勢分解為兩部分，一部分為滿足拉普拉斯方程和壁面邊界條件的特解，另一部分采用基于高精度速度勢型Boussinesq方程的數值模型求解.首先給出三維液艙中流體受強迫運動的控制方程和邊界條件，然后基于有限差分法對自由面進行空間離散，建立數值模型.通過與文獻中二維晃蕩的結果進行對比驗證了數值模型的有效性，并對三維工況的結果進行分析.

1 數學模型

考慮一個長為l、寬為b、靜止水深為h的液艙，如圖1所示，固定坐標系O-xyz位于靜止水面上，原點位于靜水面的幾何中心，z軸垂直向上.液艙底部記為z=-h(x,y)，流體自由表面記為z=η(x,y).

圖1 晃蕩運動模型Fig.1 Sloshing motion model

假定艙內為無黏無旋不可壓流體，基于勢流理論，流域內存在一個速度勢Φ.記水平速度為u，垂向速度為w，則速度勢與速度的關系為

(1)

ηt+Φxηx+Φyηy-Φz=0,z=η

(2)

(3)

(4)

?Φ/?n=v·n

(5)

其中，式(2)為自由面運動學邊界條件；式(3)為自由面動力學邊界條件；式(4)為拉普拉斯方程，Xij(i,j=x,y,z)表示物理量X對i方向求一階導數，對j方向求二階導數；式(5)為物面不可穿透條件；g為重力加速度；v=[vxvyvz]為液艙的運動速度；n為物面的法向量，由液艙內指向外.對于慣性系中定義的自由面邊界條件，利用如下關系將其轉換到隨液艙運動坐標系中：

則可以給出隨艙運動坐標系下，自由面滿足運動學和動力學邊界條件：

(8)

(9)

由于液艙的運動始終是已知量，在數值計算時，將總速度勢分為兩部分：Φ=φ+φ，其中φ是一個滿足拉普拉斯方程和物面不可穿透條件的特解.僅考慮橫蕩、縱蕩、垂蕩運動時，可以直接寫出φ的表達式如下：

φ=xvx+yvy+zvz

(10)

然后由疊加原理可導出待求解的分速度勢φ所滿足的物面邊界條件和控制方程：

(11)

將Φ=φ+φ以及式(10)引入到式(8)、(9)中，可以得到φ滿足的自由面邊界條件如下：

ηt+ηxφx+ηyφy-φz=0

(12)

(13)

2 數值計算方法

(14)

則可用新變量表示自由面條件并整理如下：

(15)

(16)

(17)

式中：參數δ?1，表示展開面是緩慢變化的.由上述遞推關系式可推導速度勢和垂向速度的表達式如下：

(18)

(19)

將表達垂向速度的無窮級數截斷在3階時，以上兩式中：

由此可得到速度勢和垂向速度的表達式，將其代入到底部運動學邊界條件a1、a3的表達式中得：

(25)

整理后可得底部方程的表達式為

(26)

數值計算過程中，空間導數計算采用有限差分法，時間步進采用四階龍格庫塔方法.計算程序基于MATLAB語言編寫，求解線性方程組使用MATLAB中的直接求解器，該求解器可根據系數矩陣的特征選擇高效的求解方法.

3 二維工況驗證

為驗證數值模型的有效性，將三維液艙的寬度設置為小量并在相同工況下與文獻中的二維結果[10]進行對比.三維液艙的寬度b=0.1 m，長度和充液高度都與二維工況相同，分別為l=1.175 m，h=0.06 m.液艙受長度方向的正弦激勵做縱蕩運動，激勵幅值A=0.39 cm，記線性理論預測的共振頻率ωr=2.042 5 s-1，激勵頻率ω′分別為0.98ωr、1.02ωr及1.08ωr.

圖2依次展示了激勵頻率為0.98ωr、1.02ωr及1.08ωr時左壁面處的液面高度變化對比圖，圖中：t′為時間與激勵周期之比，η′為自由面高度與充液高度之比.3個工況下的波高變化都呈現出強非線性現象，波峰陡峭但波谷平坦.當外部激勵頻率處于共振頻率附近時，矩形液艙內的晃蕩運動形式對外部激勵頻率的變化非常敏感，在3種工況下觀察到了完全不同的運動形式，穩(wěn)定狀態(tài)下的一個變化周期內分別包含3個波峰、2個波峰和單個波峰.相應地，在液艙內可以觀察到對應數量的行進波在艙內來回移動.

圖2 液面高度時歷對比圖Fig.2 Comparisons of time history of surface height

可以看出，三維數值計算程序能很好地捕捉自由面高度的變化和晃蕩運動形式對于外部激勵頻率變化的敏感性.

4 三維工況

計算選擇長寬相等的方型液艙，將二維液艙的長度縮小一半，三維液艙模型參數為l=b=0.587 5 m，h=0.03 m.僅考慮縱蕩、橫蕩或二者混合的正弦形式運動，A=0.1 cm，線性理論預測的艙內晃蕩運動共振頻率為ωr=2.888 1 s-1.考慮到二維情況下，外部激勵頻率在共振頻率附近變化時，晃蕩運動會出現特定的波高變化形式，三維工況計算也在共振頻率附近進行.

數值計算過程中，將靜止液面劃分成了正方形網格，x和y方向都有25個節(jié)點.在二維液艙中，波高觀測點為左壁面處第1號節(jié)點，而在三維液艙中，波高觀測點為混合運動合成矢量方向上的第1號節(jié)點.

如圖3所示，首先選取了ω′=0.98ωr，將二維工況和三維單方向運動工況進行對比.可以看出，兩種情況下自由面的高度的變化幅值和周期都是一致的，液艙的寬度并沒有對艙內流體的運動產生影響.接著將三維單方向運動工況和混合運動工況進行了對比，從圖中可以看出，兩種工況下波高變化的周期是一致的，將單方向運動的波高加倍以后，幅值的大小也有很好的一致性.

如圖4所示，計算了其他激勵頻率下，混合運動工況中波高觀測點處的波高變化.為了保證艙內的流體運動達到穩(wěn)定狀態(tài)，計算時間超過200個激勵周期.在1.02ωr激勵頻率下，觀察到了和二維工況中類似的現象：在一個變化周期內出現兩個波峰.但是單波峰的現象卻出現在1.06ωr激勵頻率下，并且在1.08ωr激勵頻率下，觀測點波高變化的非線性減弱，說明三維情況下出現強非線性行進波的激勵頻率域有所縮短.

圖3 0.98ωr下波高對比圖Fig.3 Comparisons of wave height at 0.98ωr

圖4 各頻率下波高變化圖Fig.4 Change of wave height at different frequencies

數值計算程序可以實時顯示自由面變化，以1.06ωr激勵頻率為例，圖5顯示了數值模擬過程中某一時刻整個自由面形狀.在各個運動方向上，各自形成了行進波，兩條波浪相遇的區(qū)域形成了一個波峰點.類似地，在1.02ωr激勵頻率下，會在各個方向上觀察到兩條行進波，且自由面上觀察到2×2個波峰點.而在0.98ωr激勵頻率下，會在自由面上觀察到9條行進波和3×3個波峰點.

圖5 自由面形狀圖Fig.5 Shape of free surface

5 結論

本研究基于高精度Boussinesq型速度勢方程對二維和三維液艙中在小幅激勵下的淺水晃蕩現象進行了數值模擬.主要結論有：

(1) 首次采用Boussinesq方程實現了三維淺水晃蕩模擬，建立的數值模型能很好地捕捉三維液艙內的強非線性晃蕩運動以及晃蕩運動形式對于外部激勵頻率變化的敏感性，且具有較高的計算效率.

(2) 三維數值模擬中，觀察到了4種不同的晃蕩運動形式.在各個工況下，會觀察到不同數量的行進波，行進波數量取決于外部激勵的頻率.在行進波十字交匯的區(qū)域，會疊加產生波峰.

(3) 在靠近共振頻率的工況下，即使是微幅激勵的耦合，都在艙內產生了大幅的波浪.對于遠離共振頻率的工況，相同外部激勵頻率下，與二維相比，三維情況下的非線性有所減弱.

(4) 目前數值模型僅考慮單方向的平動或兩方向耦合運動，未來將對其他運動形式(單方向轉動或耦合)作進一步的研究.