基于高精度Boussinesq方程的三維淺水晃蕩數(shù)值研究

袁心怡， 蘇 焱， 劉祖源

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

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

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

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

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

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

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

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

圖1 晃蕩運(yùn)動(dòng)模型Fig.1 Sloshing motion model

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

(1)

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

(2)

(3)

(4)

?Φ/?n=v·n

(5)

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

則可以給出隨艙運(yùn)動(dòng)坐標(biāo)系下，自由面滿足運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)邊界條件：

(8)

(9)

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

φ=xvx+yvy+zvz

(10)

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

(11)

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

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

(12)

(13)

2 數(shù)值計(jì)算方法

(14)

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

(15)

(16)

(17)

式中：參數(shù)δ?1，表示展開面是緩慢變化的.由上述遞推關(guān)系式可推導(dǎo)速度勢和垂向速度的表達(dá)式如下：

(18)

(19)

將表達(dá)垂向速度的無窮級(jí)數(shù)截?cái)嘣?階時(shí)，以上兩式中：

由此可得到速度勢和垂向速度的表達(dá)式，將其代入到底部運(yùn)動(dòng)學(xué)邊界條件a1、a3的表達(dá)式中得：

(25)

整理后可得底部方程的表達(dá)式為

(26)

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

3 二維工況驗(yàn)證

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

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

圖2 液面高度時(shí)歷對(duì)比圖Fig.2 Comparisons of time history of surface height

可以看出，三維數(shù)值計(jì)算程序能很好地捕捉自由面高度的變化和晃蕩運(yùn)動(dòng)形式對(duì)于外部激勵(lì)頻率變化的敏感性.

4 三維工況

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

數(shù)值計(jì)算過程中，將靜止液面劃分成了正方形網(wǎng)格，x和y方向都有25個(gè)節(jié)點(diǎn).在二維液艙中，波高觀測點(diǎn)為左壁面處第1號(hào)節(jié)點(diǎn)，而在三維液艙中，波高觀測點(diǎn)為混合運(yùn)動(dòng)合成矢量方向上的第1號(hào)節(jié)點(diǎn).

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

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

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

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

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

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

5 結(jié)論

本研究基于高精度Boussinesq型速度勢方程對(duì)二維和三維液艙中在小幅激勵(lì)下的淺水晃蕩現(xiàn)象進(jìn)行了數(shù)值模擬.主要結(jié)論有：

(1) 首次采用Boussinesq方程實(shí)現(xiàn)了三維淺水晃蕩模擬，建立的數(shù)值模型能很好地捕捉三維液艙內(nèi)的強(qiáng)非線性晃蕩運(yùn)動(dòng)以及晃蕩運(yùn)動(dòng)形式對(duì)于外部激勵(lì)頻率變化的敏感性，且具有較高的計(jì)算效率.

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

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

(4) 目前數(shù)值模型僅考慮單方向的平動(dòng)或兩方向耦合運(yùn)動(dòng)，未來將對(duì)其他運(yùn)動(dòng)形式(單方向轉(zhuǎn)動(dòng)或耦合)作進(jìn)一步的研究.