田鑫 ，萬德成*

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

2上海交通大學(xué)海洋工程國家重點實驗室，上海200240

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

0 引 言

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

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

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

1 數(shù)值方法

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

1.1 控制方程

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

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

1.2 核函數(shù)

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

式中：r為兩個粒子間的距離；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é)果與擴散方程的解析解相一致。

1.5 壓力泊松方程

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

1.6 自由面判斷方法

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

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

1.7 GPU加速技術(shù)

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

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

2 模型驗證

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

圖2 幾何模型Fig.2 Geometric model

液艙橫蕩激勵的運動方程為

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

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

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

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

圖4 液艙所受水平半徑方向合力的時歷曲線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的橫蕩激勵進行了數(shù)值模擬。圖5給出了不同激勵頻率下晃蕩流體的運動情況。

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

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

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

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

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

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

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

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

當(dāng)f＜1.1 Hz時，受力峰值與激勵頻率呈正相關(guān)。當(dāng)f=1.0 Hz時，流體的運動幅度較小，晃蕩較為平緩，液艙的受力幅值在100 N左右。隨著激勵頻率的增大，受力峰值增大，在1.1～1.2 Hz時達到峰值，約為220 N，此時晃蕩最劇烈，說明此時處于共振頻率。當(dāng)f＞1.2 Hz時，受力峰值與激勵頻率呈負相關(guān)。隨著激勵頻率的增大，液艙所受水平半徑方向合力的峰值減小，晃蕩幅度減小。當(dāng)f=1.3 Hz時，受力的幅值降至120 N左右。當(dāng)f=1.4 Hz時，受力幅值降至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求解器，研究了圓柱形液艙在橫蕩運動中激勵頻率對液艙受力的影響，主要得到以下結(jié)論：

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

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

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

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

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