MPS方法模擬三維圓柱形液艙晃蕩問題

田鑫 ，萬德成*

1上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海200240

2上海交通大學(xué)海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海200240

3高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240

0 引 言

液艙晃蕩是指在外部激勵作用下，艙內(nèi)裝載的部分液體所產(chǎn)生的波動及其與艙壁結(jié)構(gòu)相互作用的現(xiàn)象。隨著世界航運(yùn)業(yè)的發(fā)展以及對能源需求的不斷提高，相繼開發(fā)應(yīng)用了大型原油貨輪、液化石油氣船以及液化天然氣船等載液貨船。這些容積巨大的船裝載著巨量的液體貨物，其內(nèi)部產(chǎn)生晃蕩現(xiàn)象往往會對船舶產(chǎn)生巨大危害，因此晃蕩現(xiàn)象成為船舶與海洋工程研究的熱點(diǎn)問題。

船舶液艙的形狀多種多樣，以往對液艙晃蕩的研究主要集中在矩形液艙及薄膜型液艙上，例如，張書誼和段文洋［1］就曾采用CFD軟件Fluent對二維矩形液艙不同艙內(nèi)水深、不同激振頻率時(shí)的橫蕩進(jìn)行過數(shù)值計(jì)算和對比分析。對于其他形狀液艙的晃蕩，國內(nèi)外的研究則較少。Kobayashi等［2］通過實(shí)驗(yàn)方法，對水平圓柱形液艙的晃蕩進(jìn)行了研究，并對載液率、幅值、液艙長度等因素進(jìn)行了分析探討。Faltinsen等［3］對球形液艙中的液艙晃蕩現(xiàn)象予以了研究，并對其中的線性晃蕩進(jìn)行了理論分析。Gavrilyuk等［4］對直立環(huán)形圓柱的晃蕩進(jìn)行了理論及實(shí)驗(yàn)研究，分析了環(huán)形擋板的減晃作用。劉戈等［5-6］通過實(shí)驗(yàn)，研究了獨(dú)立C型液艙在縱搖激勵下的晃蕩，討論并分析了其晃蕩時(shí)自由液面的波形特征，并對該型液艙的頻域共振特性進(jìn)行了研究，發(fā)現(xiàn)0.5倍理論固有頻率為晃蕩起始頻率，終止頻率會隨著液位的升高而下降。

本文將基于課題組自主研發(fā)的無網(wǎng)格粒子法求解器MPSGPU-SJTU，對橫蕩激勵下圓柱形液艙中的晃蕩現(xiàn)象進(jìn)行研究。首先，模擬50%充水率下，激勵頻率為一階固有頻率時(shí)的晃蕩現(xiàn)象，通過對比實(shí)驗(yàn)結(jié)果，驗(yàn)證模擬的準(zhǔn)確性。然后，對不同激勵頻率下的晃蕩現(xiàn)象進(jìn)行計(jì)算，對比不同激勵頻率下的液艙受力與流場情況，以分析激勵頻率對圓柱形液艙晃蕩的影響。

1 數(shù)值方法

MPSGPU-SJTU求解器是在本課題組之前開發(fā)的CPU并行求解器MLParticle-SJTU［7-9］的基礎(chǔ)上，引入 GPU（Graphics Process Unit）加速技術(shù)開發(fā)的新一代求解器，其數(shù)值方法為經(jīng)張雨新［10］改進(jìn)的移動粒子半隱式（Moving Particle Semi-implicit，MPS）方法，GPU加速基于CUDA平臺實(shí)現(xiàn)。

1.1 控制方程

在MPS方法中，控制方程包括連續(xù)性方程和N-S方程，對于不可壓縮流體，其形式如下：

式中：ρ為流體密度；P為壓力；V為速度向量；Fm為質(zhì)量力，一般為重力；ν為流體的運(yùn)動粘性系數(shù)；t為時(shí)間。

1.2 核函數(shù)

在粒子法中，控制方程將被寫成粒子形式，兩個(gè)粒子間的相互影響是通過核函數(shù)來實(shí)現(xiàn)的，其與傳統(tǒng)的核函數(shù)不同。MPSGPU-SJTU求解器采用的核函數(shù)為

式中：r為兩個(gè)粒子間的距離；rε為粒子作用域的半徑。

1.3 梯度模型

MPSGPU-SJTU求解器采用的動量守恒的梯度模型為

式中：D為空間維數(shù)；ri，rj分別為粒子i及其鄰居粒子j的坐標(biāo)；Pi，Pj分別為粒子i和j的壓力；n0為初始粒子數(shù)密度。

1.4 拉普拉斯算子模型

在MPS方法中，拉普拉斯算子模型如下式所示：

式中，φi和φj為物理量φ在i，j處的值。

其中，

系數(shù)λ的引入是為了使數(shù)值結(jié)果與擴(kuò)散方程的解析解相一致。

1.5 壓力泊松方程

MPSGPU-SJTU求解器采用的壓力泊松方程是由Tanaka 和 Masunaga［11］提出的混合源項(xiàng)法（Mixed source term method）。該方法結(jié)合了速度散度和粒子數(shù)密度，后來被Lee等［12］寫成了更為合理的表達(dá)形式：

1.6 自由面判斷方法

在MPS方法中，對自由液面的準(zhǔn)確判斷對計(jì)算的精度和穩(wěn)定性來說十分重要。在本文采用的自由面判斷方法中，首先定義矢量

式中，F(xiàn)為粒子數(shù)密度的不確定度。計(jì)算F的模，當(dāng)粒子滿足時(shí)，即被判定為自由面粒子，其中α=0.9，為一參數(shù)。

1.7 GPU加速技術(shù)

GPU具有強(qiáng)大的計(jì)算能力，在并行計(jì)算中，應(yīng)用GPU能獲得顯著的效果。NVIDIA公司推出的CUDA技術(shù)為GPU編程提供了方便的平臺，MPSGPU-SJTU求解器即是基于此平臺而開發(fā)。本研究對MPS方法的計(jì)算流程進(jìn)行了合理的設(shè)計(jì)，實(shí)現(xiàn)了MPS方法主要計(jì)算過程的GPU加速。計(jì)算流程如圖1所示。

圖1 MPSGPU-SJTU求解器計(jì)算流程圖Fig.1 The flow chart of MPSGPU-SJTU solver

2 模型驗(yàn)證

本文構(gòu)造的數(shù)值模型是參考Kobayashi等［2］的模型而建立，其物理模型如圖2所示。液艙為圓柱體，圓柱底面直徑0.47 m，長0.94 m，載液率為50%。

圖2 幾何模型Fig.2 Geometric model

液艙橫蕩激勵的運(yùn)動方程為

式中：Ax為激勵幅值，Ax=0.015 m；f為激勵頻率，Hz。

構(gòu)造的數(shù)值模型如圖3所示。流體粒子數(shù)為79 258，總粒子數(shù)為137 380，粒子間距0.01 m。水的密度ρ=1 000 kg/m3，運(yùn)動 粘性系 數(shù)ν=10-6m2/s，重力加速度g=9.81 m/s2，時(shí)間步長dt=0.000 5 s。

圖3 數(shù)值模型Fig.3 Numerical model

本節(jié)對f=1.2 Hz、橫蕩激勵幅值A(chǔ)x=0.015 m的晃蕩進(jìn)行了模擬。圖4給出了橫蕩激勵作用下，實(shí)驗(yàn)與模擬中液艙所受水平半徑方向合力Fy的時(shí)歷曲線對比。從圖中可以看出，數(shù)值模擬給出的合力曲線與實(shí)驗(yàn)結(jié)果吻合較好，兩者的變化趨勢、周期峰值基本一致。可見，采用求解器MPSGPU-SJTU能夠很好地預(yù)測液艙所受合力。

圖4 液艙所受水平半徑方向合力的時(shí)歷曲線Fig.4 Time history curves of horizontal radial direction resultant force applied to tank

3 激勵頻率對晃蕩的影響

本節(jié)分別對激勵頻率f=1.0，1.1，1.2，1.3，1.4 Hz，幅值A(chǔ)x=0.015 m的橫蕩激勵進(jìn)行了數(shù)值模擬。圖5給出了不同激勵頻率下晃蕩流體的運(yùn)動情況。

圖5 不同激勵頻率下的流場Fig.5 Flow field at different frequencies

由圖中可以看出，激勵頻率對晃蕩流體的流動特征具有顯著影響。流場隨頻率變化的過程可分為以下4個(gè)階段：

1）如圖5（a）所示，當(dāng)激勵頻率f=1.0 Hz時(shí)，晃蕩的幅度很小，自由面非常光滑，流體運(yùn)動比較平穩(wěn)。

2）隨著激勵頻率的增加，晃蕩波的波陡迅速增大，反應(yīng)十分劇烈。當(dāng)f=1.1和1.2 Hz時(shí)（圖5（b）和圖5（c））出現(xiàn)了沖頂現(xiàn)象，晃蕩波沿壁面爬升至艙頂后下落，形成一股射流猛烈拍擊在自由面上，從而發(fā)生了波面破碎和液體飛濺現(xiàn)象。整個(gè)晃蕩過程非常劇烈，并伴隨著一些復(fù)雜的流動現(xiàn)象，具有較強(qiáng)的非線性。其中當(dāng)f=1.2 Hz時(shí)晃蕩波的運(yùn)動更加劇烈，其沿艙壁爬升的高度更高，并且流體動能更大，一部分液體甚至飛濺到了另一側(cè)艙壁。

3）隨著激勵頻率的繼續(xù)增大，晃蕩波的波高和波陡迅速減小。當(dāng)f=1.3 Hz時(shí)（圖5（d）），晃蕩雖然仍較劇烈，但晃蕩波的動能已不足以支撐其爬升至艙頂。波浪在砰擊壁面并爬升很小的高度后，波面翻卷沖擊自由面，發(fā)生了輕微的破波現(xiàn)象。在這一過程中，液體的運(yùn)動幅度是隨著激勵頻率的增加而減小的。

4）當(dāng)f=1.4 Hz時(shí)，液艙運(yùn)動較快，但晃蕩波的波幅繼續(xù)減小，流動變得平緩，不再出現(xiàn)波浪的翻卷破碎現(xiàn)象。

圖6所示為不同激勵頻率下的液艙受力時(shí)歷曲線。從整體上看，其呈現(xiàn)出與激勵頻率一致的周期性。但從細(xì)節(jié)上可以看到，液艙受力曲線的形狀和峰值也隨著不同激勵頻率而發(fā)生變化。當(dāng)激勵頻率f=1.0 Hz時(shí)，受力曲線為單峰特征（圖6（a）），受力的峰值相對較小，其大小約為100 N，此時(shí)晃蕩現(xiàn)象較為平穩(wěn)，沒有發(fā)生拍擊現(xiàn)象。當(dāng)激勵頻率增大至1.1或1.2 Hz時(shí)，晃蕩產(chǎn)生共振現(xiàn)象，受力曲線出現(xiàn)了雙峰現(xiàn)象。其中第1個(gè)波峰達(dá)到了250 N甚至是300 N以上，這是流體直接沖擊壁面形成；而第2個(gè)受力峰值大小約為100～150 N，是由爬升的流體回落，對底部流體產(chǎn)生拍擊作用所致。隨著激勵頻率繼續(xù)增大，受力峰值開始下降，當(dāng)f=1.3 Hz時(shí)，受力曲線仍然存在雙峰特征，但2個(gè)峰值已經(jīng)非常接近，此時(shí)，流場仍然存在明顯的非線性流動。當(dāng)激勵頻率增大至1.4 Hz時(shí)，只在一些周期內(nèi)出現(xiàn)了第2個(gè)受力峰值，這是由艙壁附近微弱的破波現(xiàn)象所致。

圖7給出了受力峰值隨激勵頻率的變化規(guī)律。從圖中可以看出，隨著激勵頻率的增加，受力峰值是先增大后減小。

當(dāng)f＜1.1 Hz時(shí)，受力峰值與激勵頻率呈正相關(guān)。當(dāng)f=1.0 Hz時(shí)，流體的運(yùn)動幅度較小，晃蕩較為平緩，液艙的受力幅值在100 N左右。隨著激勵頻率的增大，受力峰值增大，在1.1～1.2 Hz時(shí)達(dá)到峰值，約為220 N，此時(shí)晃蕩最劇烈，說明此時(shí)處于共振頻率。當(dāng)f＞1.2 Hz時(shí)，受力峰值與激勵頻率呈負(fù)相關(guān)。隨著激勵頻率的增大，液艙所受水平半徑方向合力的峰值減小，晃蕩幅度減小。當(dāng)f=1.3 Hz時(shí)，受力的幅值降至120 N左右。當(dāng)f=1.4 Hz時(shí)，受力幅值降至90 N。

圖6 不同激勵頻率下的液艙受力Fig.6 Forces acting on tanks at different excitation frequencies

圖7 液艙受力峰值隨激勵頻率的變化Fig.7 Variation of peak force in tank with excitation frequency

4 結(jié) 論

本文采用課題組自主開發(fā)的MPSGPU-SJTU求解器，研究了圓柱形液艙在橫蕩運(yùn)動中激勵頻率對液艙受力的影響，主要得到以下結(jié)論：

1）MPSGPU-SJTU求解器能夠較好地模擬圓柱形液艙內(nèi)的液體晃蕩現(xiàn)象，能很好地捕捉晃蕩現(xiàn)象中的液體飛濺和波浪翻卷破碎等現(xiàn)象，并且也可以較為準(zhǔn)確地計(jì)算液艙受力。

2）不同激勵頻率下，晃蕩的波形區(qū)別很大。當(dāng)f=1.1與1.2 Hz時(shí)，有沖頂現(xiàn)象出現(xiàn)；當(dāng)f=1.3 Hz時(shí)，有輕微的破波現(xiàn)象出現(xiàn)；當(dāng)f=1.1和1.4 Hz時(shí)，流動較平緩。

3）當(dāng)激勵頻率在固有頻率附近時(shí)，晃蕩最為激烈，液艙的受力幅值最大，但隨著激勵頻率遠(yuǎn)離固有頻率，晃蕩逐漸減弱。

4）在船舶運(yùn)營中，要避免在液艙固有頻率附近產(chǎn)生激勵振動。

本文研究了激勵頻率對圓柱形液艙晃蕩的影響，給出了液艙整體受力的結(jié)果，可為圓柱形液艙的設(shè)計(jì)應(yīng)用提供有效依據(jù)。在實(shí)際液艙晃蕩問題中，由于很多破壞是由于局部應(yīng)力過大所造成，所以在今后的研究中，會計(jì)算液艙的局部壓力，對液艙晃蕩中的現(xiàn)象與載荷進(jìn)行更為細(xì)致的研究。