寧德志，宋偉華，2，滕 斌

（1.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.中交第四航務(wù)工程勘察設(shè)計(jì)院有限公司，廣州 510230）

縱搖容器中液體晃蕩的非線性數(shù)值模擬

寧德志1，宋偉華1，2，滕 斌1

（1.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.中交第四航務(wù)工程勘察設(shè)計(jì)院有限公司，廣州 510230）

在海洋工程領(lǐng)域，液體晃蕩是一種普遍存在的物理現(xiàn)象。對(duì)于船舶而言，轉(zhuǎn)動(dòng)比平動(dòng)有著更重要的影響。該文針對(duì)縱搖容器中的液體晃蕩問題，采用高階邊界元方法建立自由水面滿足完全非線性邊界條件的時(shí)域數(shù)學(xué)模型。通過大地坐標(biāo)系和隨體坐標(biāo)系之間的坐標(biāo)變換，使得計(jì)算域僅控制在隨體坐標(biāo)系內(nèi)。求解中采用半混合歐拉—拉格朗日方法追蹤流體瞬時(shí)水面，運(yùn)用四階龍格庫(kù)塔方法更新下一時(shí)間步的波面和速度勢(shì)。通過與已發(fā)表試驗(yàn)和數(shù)值結(jié)果的對(duì)比，驗(yàn)證了建立模型的準(zhǔn)確性。進(jìn)而開展大量數(shù)值試驗(yàn)研究容器縱搖運(yùn)動(dòng)頻率、縱搖轉(zhuǎn)動(dòng)中心和容器中布置一豎直隔板對(duì)晃動(dòng)波面與荷載的影響。

液體晃蕩；縱搖；完全非線性；高階邊界元；隔板

0 引 言

在海洋工程領(lǐng)域中，液體晃蕩是一種普遍存在的物理現(xiàn)象。對(duì)于海上航行的船舶而言，縱搖運(yùn)動(dòng)對(duì)其穩(wěn)定性的影響遠(yuǎn)遠(yuǎn)大于平動(dòng)運(yùn)動(dòng)。尤其是當(dāng)運(yùn)動(dòng)幅值很大或外界運(yùn)動(dòng)頻率接近液艙的固有頻率時(shí)，液體晃蕩便會(huì)趨于劇烈，在內(nèi)壁上產(chǎn)生極大的作用力，從而對(duì)結(jié)構(gòu)造成危害，影響船舶的安全性和穩(wěn)定性。

對(duì)于縱搖液艙中的液體晃蕩，國(guó)內(nèi)外許多學(xué)者已經(jīng)進(jìn)行了研究。Faltinsen[1]早在1974年就對(duì)二維矩形液艙作縱搖運(yùn)動(dòng)時(shí)的液體晃蕩問題進(jìn)行了非線性理論分析。此非線性理論解析解為后來學(xué)者進(jìn)行數(shù)值模擬提供了很好的參照。Nakayama和Washizu[2]采用數(shù)值方法對(duì)二維矩形液艙縱搖運(yùn)動(dòng)的液體晃蕩進(jìn)行了模擬研究，他們利用有限元方法對(duì)空間進(jìn)行網(wǎng)格劃分，采用有限差分方法來計(jì)算時(shí)間遞進(jìn)。隨后他們又運(yùn)用邊界元方法對(duì)二維縱搖矩形容器中的液體晃蕩進(jìn)行了研究[3]。同時(shí)，Chen[4-8]對(duì)于液艙縱搖運(yùn)動(dòng)的液體晃蕩也做了大量工作。他針對(duì)理想流體的液體晃蕩問題展開研究，分析了運(yùn)動(dòng)頻率等因素的影響作用；同時(shí)還建立了粘性流模型，并將粘性流結(jié)果與理想流體結(jié)果進(jìn)行對(duì)比分析。此外，他還針對(duì)復(fù)合運(yùn)動(dòng)液艙內(nèi)的液體晃蕩進(jìn)行了研究，即液艙作平動(dòng)和縱搖聯(lián)合運(yùn)動(dòng)。

近幾年來，Akyildiz和Unal[9-11]又針對(duì)矩形縱搖運(yùn)動(dòng)液艙中的液體晃蕩進(jìn)行了研究。他們分別采取了物理實(shí)驗(yàn)和數(shù)值模擬兩種方法對(duì)液體晃蕩產(chǎn)生的動(dòng)水壓強(qiáng)進(jìn)行了研究對(duì)比，并且考慮了縱搖的頻率和角度幅值對(duì)液體晃蕩的影響。除此之外，Armenio[12]、Bae[13]、Celebi[14]和Xue[15]等學(xué)者也在這方面做了大量的研究工作。同時(shí)，國(guó)內(nèi)學(xué)者肖龍飛等[16]、尹立中等[17]、洪亮等[18]、徐國(guó)徽等[19]也在這個(gè)領(lǐng)域有著較為深入的研究。然而，與平動(dòng)運(yùn)動(dòng)液艙中液體晃蕩的研究相比，人們對(duì)縱搖運(yùn)動(dòng)液艙中液體晃蕩的研究還是相對(duì)較少，尤其是在減晃運(yùn)動(dòng)研究方面。

本文通過高階邊界元方法建立自由水面滿足完全非線性邊界條件的時(shí)域數(shù)學(xué)模型，針對(duì)液艙在不同轉(zhuǎn)動(dòng)中心和運(yùn)動(dòng)頻率情況下的液體晃蕩進(jìn)行了數(shù)值試驗(yàn)，并分析了各因素的影響。在此基礎(chǔ)上，考慮到對(duì)液體晃蕩的抑制，本文還針對(duì)液艙中安置垂直隔板的情況進(jìn)行了研究，給出了不同隔板位置、隔板長(zhǎng)度情況下液體晃蕩的波面幅值和液艙所受壓力的變化趨勢(shì)，可為實(shí)際工程提供一定的參考借鑒。

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

1.1 控制方程和邊界條件

如圖1所示，一個(gè)長(zhǎng)度為L(zhǎng)、水深為H的矩形容器在外界激勵(lì)作用下繞O點(diǎn)進(jìn)行縱搖運(yùn)動(dòng)。定義兩個(gè)笛卡爾坐標(biāo)系：一個(gè)是固定大地坐標(biāo)系O0X0Y0；另一個(gè)是運(yùn)動(dòng)的隨體坐標(biāo)系oxy。其中，x軸平行于容器靜止時(shí)的自由水面且向右為正方向，y軸和容器的中心線重合且向上為正方向。當(dāng)容器靜止時(shí)，兩個(gè)坐標(biāo)系是重合的。同時(shí)，定義靜止水面和x軸間的距離為e，且當(dāng)x軸在靜止水面以下時(shí)e值取為正。

假設(shè)流體是無粘、不可壓縮且運(yùn)動(dòng)無旋的理想流體，則流體運(yùn)動(dòng)滿足Laplace方程，

式中：φ表示速度勢(shì)。

運(yùn)動(dòng)角速度ω（）t定義為

式中：θ（）

t表示液艙偏離靜止位置的角度，且順時(shí)針偏離為正值。

圖1 示意圖Fig.1 Definition sketch

同時(shí)假設(shè)：一是自由表面上的壓強(qiáng)等于大氣壓強(qiáng)的自由表面邊界條件；二是自由表面上的流體質(zhì)點(diǎn)將始終保持在自由表面上。本文中將大氣壓強(qiáng)設(shè)定為0?；诘谝粭l假設(shè)，自由表面的動(dòng)力學(xué)邊界條件可寫為[2]

式中：η（x， ）

t是自由表面偏離靜止水面的位移，g是重力加速度?；诘诙l假設(shè)，自由表面的運(yùn)動(dòng)學(xué)邊界條件可寫為[2]

式中：nx和ny是自由表面上的外法向量與x、y軸夾角的余弦值。?/?n表示在外法向上的偏導(dǎo)數(shù)。

在固壁上，邊界條件可寫為如下形式[2]：

當(dāng)t＝0時(shí)，容器靜止不動(dòng)，相應(yīng)的初始條件可寫為

此外，作用在液艙上的動(dòng)水壓力可由Bernoulli方程求得：

進(jìn)而通過在容器瞬時(shí)濕表面上進(jìn)行壓強(qiáng)積分求出作用在容器上的荷載。

1.2 數(shù)值求解

為求解上節(jié)中的數(shù)學(xué)方程，在整個(gè)流域內(nèi)對(duì)速度勢(shì)應(yīng)用格林第二定理，可得到如下邊界積分方程[20]：

式中：p為源點(diǎn)，q為場(chǎng)點(diǎn)，C為固角系數(shù)，S為流域邊界，包括自由水面邊界和固體邊界。本文用三節(jié)點(diǎn)高階邊界元離散計(jì)算邊界成一些曲線單元，單元內(nèi)任一點(diǎn)的幾何坐標(biāo)和速度勢(shì)等物理量可以用形狀函數(shù)插值得到，這與曲面單元內(nèi)的做法是相同的。G是簡(jiǎn)單格林函數(shù)，可寫為

式中：r1=，H為水深。

接下來，積分方程（8）經(jīng)過三節(jié)點(diǎn)曲線高階邊界單元離散后，便可以通過求解線性方程組得到相關(guān)未知量的值。計(jì)算中認(rèn)為當(dāng)前時(shí)刻物面Sn上的速度勢(shì)法向?qū)?shù)和自由水面Sf上的速度勢(shì)是已知的，根據(jù)積分方程計(jì)算當(dāng)前時(shí)刻物面Sn上的速度勢(shì)和自由水面Sf上的速度勢(shì)法向?qū)?shù)，然后應(yīng)用四階Runga-Kutta法，根據(jù)自由水面條件式（3）、（4）計(jì)算下一時(shí)刻的水質(zhì)點(diǎn)位置和自由水面Sf上的速度勢(shì)，再對(duì)自由水面重新劃分網(wǎng)格，重新應(yīng)用積分方程計(jì)算下一時(shí)刻物面上的速度勢(shì)和自由水面上的速度勢(shì)法向?qū)?shù)。這樣計(jì)算周而復(fù)始，直到計(jì)算結(jié)束[21]。

2 模型驗(yàn)證

作為算例，首先考慮一個(gè)長(zhǎng)度為0.9m、水深為0.6m的矩形容器，運(yùn)動(dòng)方程為θ=θ0cosωt。其中，θ0是角度幅值，ω是運(yùn)動(dòng)頻率。選取θ0和ω為0.8度和5.5 rad/s。轉(zhuǎn)動(dòng)中心設(shè)定在靜止水面上，即e＝0。由運(yùn)動(dòng)方程可知，當(dāng)t＝0時(shí)，容器有一個(gè)0.8度的傾角。初始波面設(shè)為靜止不動(dòng)，即初始波面與軸也有一個(gè)0.8度的夾角。通過數(shù)值收斂性分析，選取時(shí)間步長(zhǎng)和網(wǎng)格大小為t＝0.01 s、△x＝0.045m、△z＝0.05m。圖2給出了容器右側(cè)壁處的波面時(shí)間歷程圖及本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果[2]的對(duì)比。由圖可知，本文數(shù)值結(jié)果與Nakayama和Washizu的數(shù)值結(jié)果吻合得非常好，且波面變化趨于共振。

接下來考慮另一個(gè)算例，矩形容器長(zhǎng)度L＝0.92m、水深H＝0.31m。運(yùn)動(dòng)方程為θ=θ0sinωt。其中，θ0是角度幅值，ω是運(yùn)動(dòng)頻率。此例中，選取θ0和ω為8.0度和2.0 rad/s。轉(zhuǎn)動(dòng)中心設(shè)定在靜止水面上，即e＝0。初始波面是靜止不動(dòng)的。通過數(shù)值收斂性分析，選取時(shí)間步長(zhǎng)和網(wǎng)格大小為t＝0.01 s，△x＝0.046m，△z＝0.052m。圖3給出了容器右側(cè)壁處的波面時(shí)間歷程圖及本文數(shù)值結(jié)果與Akyildiz實(shí)驗(yàn)數(shù)據(jù)[11]、Akyildiz的CFD數(shù)值結(jié)果[11]的對(duì)比。由圖可知，本文數(shù)值結(jié)果比Akyildiz的數(shù)值結(jié)果更接近實(shí)驗(yàn)數(shù)據(jù)。

圖2 余弦縱搖運(yùn)動(dòng)容器右側(cè)壁處的波面時(shí)間歷程Fig.2 Time series of surface elevation at the right wall subject to cosine pitch excitation

圖3 正弦縱搖運(yùn)動(dòng)容器右側(cè)壁處的波面時(shí)間歷程Fig.3 Time series of surface elevation at the right end subject to sine pitch excitation

圖4 右側(cè)壁底部上0.075m處的壓強(qiáng)時(shí)間歷程Fig.4 Time series of pressure ata point0.075m from the rightwall-bottom

圖5 t=5.52 s時(shí)右側(cè)壁上的壓強(qiáng)分布Fig.5 Distribution of pressure along the right wall at t=5.25 s

在上述波面驗(yàn)證的基礎(chǔ)上，進(jìn)一步開展數(shù)值模型的壓強(qiáng)驗(yàn)證。選取容器長(zhǎng)度L為1.0 m，水深H為0.35 m，運(yùn)動(dòng)方程為θ=θ0cosωt，θ0和ω為5.0度和5.19 rad/s。轉(zhuǎn)動(dòng)中心在靜止水面以下0.1 m，即e＝0.1m。圖4給出了容器右側(cè)壁上一點(diǎn)的壓強(qiáng)時(shí)間歷程圖及本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果[3]的對(duì)比，這點(diǎn)距離液艙底0.075m。由圖可知，本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果吻合良好。同時(shí)，圖5給出了t＝5.52 s時(shí)液艙右側(cè)壁上的壓強(qiáng)分布圖及本文結(jié)果與Nakayama和Washizu數(shù)值結(jié)果[3]、Higuchi實(shí)驗(yàn)數(shù)據(jù)[22]的對(duì)比。由圖可知，本文數(shù)值結(jié)果與Nakayama和Washizu數(shù)值結(jié)果吻合很好，與實(shí)驗(yàn)數(shù)據(jù)對(duì)比也吻合良好。

通過上述算例計(jì)算，驗(yàn)證了本文數(shù)值模型能夠準(zhǔn)確模擬縱搖容器中液體晃蕩的波面變化和壓強(qiáng)變化。

3 計(jì)算及分析

3.1 運(yùn)動(dòng)頻率對(duì)液體晃蕩的影響

眾所周知，運(yùn)動(dòng)頻率在很大程度上決定著液體晃蕩的劇烈程度。當(dāng)運(yùn)動(dòng)頻率接近固有頻率時(shí)，液體晃蕩將會(huì)趨于共振，在結(jié)構(gòu)上產(chǎn)生極大的壓力。因此，本文將針對(duì)運(yùn)動(dòng)頻率對(duì)液體晃蕩的影響進(jìn)行模擬研究。一般情況下，液體晃蕩的強(qiáng)度主要通過波面幅值和作用在結(jié)構(gòu)上的壓力來反映。本文中，容器左側(cè)壁處的波面幅值和作用在液艙上的水平力將被用來衡量液體晃蕩強(qiáng)度。其中，水平力方向與x軸平行，且指向x軸正方向?yàn)檎榉奖闫鹨?，本文統(tǒng)一規(guī)定：液艙尺寸L＝1.0 m，H＝0.5m，運(yùn)動(dòng)方程為θ=θ0sinωt，初始波面為靜止水面。通過數(shù)值收斂性驗(yàn)證，時(shí)間步長(zhǎng)和網(wǎng)格分別為△t=0.01 s和△x=△z= 0.05m。

對(duì)于L=1.0 m、H＝0.5 m的矩形容器，其固有頻率可由求得。其中，n表示固有頻率的階數(shù)。各階固有頻率分別為ω1＝5.32 rad/s、ω2＝7.84 rad/s、ω3＝9.61 rad/s、ω4＝11.10 rad/s、ω5＝12.41 rad/s和ω6＝13.60 rad/s。下面將針對(duì)運(yùn)動(dòng)頻率是固有頻率的情況進(jìn)行模擬分析。選取了6個(gè)固有頻率作為運(yùn)動(dòng)頻率，分別為ω1、ω2、ω3、ω4、ω5和ω6。轉(zhuǎn)角幅值統(tǒng)一為1.0度。轉(zhuǎn)動(dòng)中心固定在靜止水面上，即e＝0。

圖6給出了6個(gè)運(yùn)動(dòng)頻率下液艙左側(cè)壁處波面時(shí)間歷程圖。在圖6（a）、（c）、（e）中，波面幅值隨著時(shí)間逐漸增大；同時(shí)，運(yùn)動(dòng)頻率越大，波面幅值增長(zhǎng)的速率越小，波面幅值與時(shí)間的線性關(guān)系愈不明顯。在圖6（b）、（d）、（f）中，波面幅值隨著時(shí)間沒有增大的趨勢(shì)，而是保持在一個(gè)穩(wěn)定的數(shù)值，亦即沒有發(fā)生共振現(xiàn)象。因此，相對(duì)于偶數(shù)階固有頻率來說，液體晃蕩對(duì)奇數(shù)階固有頻率更為敏感并且當(dāng)運(yùn)動(dòng)頻率是奇數(shù)階固有頻率時(shí)，頻率越小，其對(duì)液體晃蕩的影響越大，共振現(xiàn)象越明顯。

圖6 不同運(yùn)動(dòng)頻率情況下的容器左側(cè)壁處波面時(shí)間歷程Fig.6 Time series of surface elevation at the leftwall at differentexcitation frequency

3.2 轉(zhuǎn)動(dòng)中心對(duì)液體晃蕩的影響

選定三種有代表性的轉(zhuǎn)動(dòng)中心位置e分別為-0.2 m、0、0.2 m。運(yùn)動(dòng)頻率ω和轉(zhuǎn)角幅值分別為ω1（容器的一階固有頻率為5.316 rad/s）和1.0度。

圖7 容器左側(cè)壁處的波面時(shí)間歷程 Fig.7 Time series of surface elevation at the leftwall

圖8 作用在容器上的水平力時(shí)間歷程Fig.8 Time series of horizontal force acted on the container

圖7給出了三種情況下的容器左側(cè)壁處波面時(shí)間歷程。由圖可知，當(dāng)轉(zhuǎn)動(dòng)中心在靜止水面（e＝0）時(shí)，容器左側(cè)壁處的波面幅值最小。當(dāng)轉(zhuǎn)動(dòng)中心偏離靜水面，轉(zhuǎn)動(dòng)中心在靜止水面以上（e＝-0.2m）情況下的容器左側(cè)壁處的波面幅值要遠(yuǎn)遠(yuǎn)大于轉(zhuǎn)動(dòng)中心在靜止水面下面（e＝0.2m）情況下的。同時(shí)，圖8給出了作用在容器上水平力的時(shí)間歷程圖。由圖可知，轉(zhuǎn)動(dòng)中心在靜止水面（e＝0）時(shí)的水平力最小。當(dāng)轉(zhuǎn)動(dòng)中心距離靜止水面0.2m時(shí)，轉(zhuǎn)動(dòng)中心在靜止水面以上（e＝-0.2m）情況下的水平力要明顯大于轉(zhuǎn)動(dòng)中心在靜止水面下面（e＝0.2m）情況下的，這同圖7反映了同樣的實(shí)質(zhì)，即轉(zhuǎn)動(dòng)中心在靜止水面以上更能激起容器中的液體晃蕩。

3.3 隔板對(duì)液體晃蕩的影響

對(duì)于在海上航行的液貨船，縱搖是影響其內(nèi)液體晃蕩強(qiáng)度的重要因素。通常，人們?cè)诖搩?nèi)部加上不同形式的隔板來抑制液體的劇烈晃蕩，從而保證船體的穩(wěn)定性。本文中主要考慮垂直隔板對(duì)縱搖容器中液體晃蕩的影響?？紤]液艙具有一個(gè)長(zhǎng)為0.3m、厚為0.01m的隔板。隔板是剛性的，固定在容器底部的左側(cè)。隔板中心到容器中心的距離定義為D。運(yùn)動(dòng)頻率ω為無隔板液艙的一階固有頻率ω1，轉(zhuǎn)角幅值為1.0度。隔板上的網(wǎng)格大小為△x=0.01m和△z=0.02m。圖9給出了隔板在不同位置D情況下左右側(cè)壁處的最大波高值。由圖可知，左側(cè)壁處的波高值要大于右側(cè)壁處的。兩側(cè)壁處的波高值都隨著D的增大而逐漸增大。當(dāng)隔板位置由0移動(dòng)到0.3m時(shí)，波高值只增大很小的幅度。而當(dāng)隔板位置由0.3m增大到0.4m時(shí)，波高值卻急劇增大。圖10給出了隔板在不同位置D情況下作用在液艙上的最大水平力。由圖可知，最大水平力隨著D的增大而逐漸增大。當(dāng)隔板位置由0移動(dòng)到0.3m時(shí)，最大水平力增大的幅值非常小。當(dāng)隔板位置從0.3m增大到0.4m時(shí)，最大水平力卻增幅明顯。這與圖9中波面值與D的關(guān)系是一致的。 因此，當(dāng)隔板位置在區(qū)間（0，0.3L）內(nèi)時(shí)，其位置對(duì)液體晃蕩的影響不太明顯。當(dāng)隔板位置大于0.3L時(shí)，液體晃蕩將對(duì)其位置將非常敏感。

圖9 不同D情況下的左右側(cè)壁處最大波面高度Fig.9 Maximum surface-elevation at leftand rightwalls vs. distance between the baffle and bottom center

圖10 不同D情況下的容器所受的最大水平力Fig.10 Maximum horizontal-loads on the container vs. distance between the baffle and bottom center

4 結(jié) 論

針對(duì)縱搖液艙中液體晃蕩問題，本文利用高階邊界元方法建立模擬帶有簡(jiǎn)晃隔板縱搖容器液體晃蕩的完全非線性時(shí)域數(shù)學(xué)模型，通過坐標(biāo)轉(zhuǎn)換，使得計(jì)算域僅僅控制在隨體坐標(biāo)系內(nèi)。通過縱搖容器中波面變化和作用在容器上荷載變化的模擬，并與已發(fā)表試驗(yàn)和數(shù)值結(jié)果的對(duì)比，驗(yàn)證了本文模型的準(zhǔn)確性。進(jìn)一步開展數(shù)值模擬研究發(fā)現(xiàn)，相對(duì)于容器的偶數(shù)階固有頻率來說，液體晃蕩對(duì)奇數(shù)階固有頻率更為敏感，并且當(dāng)運(yùn)動(dòng)頻率是奇數(shù)階固有頻率時(shí)，頻率越小，其對(duì)液體晃蕩的影響越大；同時(shí)，轉(zhuǎn)動(dòng)中心在靜水面以上容器中的液體晃蕩效果更加顯著，而轉(zhuǎn)動(dòng)中心在靜水面上液體晃動(dòng)效果最弱。此外，本文還研究了安置垂直隔板液艙中的液體晃蕩。對(duì)于隔板長(zhǎng)度一定的情況，隔板位置D在區(qū)間（0，0.3L）內(nèi)時(shí)，其位置對(duì)液體晃蕩的影響不太明顯；當(dāng)D大于0.3L時(shí)，液體晃蕩將對(duì)其位置將非常敏感。

參 考 文 獻(xiàn)：

[1]Faltisen OM.A nonlinear theory of sloshing in rectangular containers[J].Journal of Ship Research,1974,18:224-241.

[2]Nakayama T,Washizu K.Nonlinear analysis of liquid motion in a container subjected to forced pitching oscillation[J].International Journal for NumericalMethods in Engineering,1980,15:1207-1220.

[3]Nakayama T,Washizu K.The boundary elementmethod applied to the analysis of two-dimensional nonlinear sloshing problems[J].International Journal for NumericalMethods in Engineering,1981,17:1631-1646.

[4]Chen B F,Chiang H W.Complete 2D and fully nonlinear analysis of ideal fluid in containers[J].Journal of Engineering Mechanics,1999,125(1):70-78.

[5]Chen B F,Chiang HW.Complete two-dimensional analysis of sea-wave-induced fully non-linear sloshing fluid in a rigid floating container[J].Ocean Engineering,2000,27:953-977.

[6]Chen B F,Nokes R.Time-independent finite difference analysis of fully non-linear and viscous fluid sloshing in a rectangular container[J].Journal of Computational Physics,2005,209:47-81.

[7]Chen B F,Asce M.Viscous fluid in container under coupled surge,heave and pitchmotions[J].Journal ofWaterway, Port,Coastal and Ocean Engineering,2005,131(5):239-256.

[8]Chen B F,Wu CH.Effects of excitation angle and coupled heave-surge-swaymotion on fluid sloshing in a three-dimensional container[J].Journal ofMarine Science and Technology,2011,16:22-50.

[9]Akyildiz H,Unal E.Experimental investigation of pressure distribution on a rectangular container due to the liquid sloshing[J].Ocean Engineering,2005,32:1503-1516.

[10]Akildiz H,Unal E.Sloshing in a three-dimensional rectangular container:numerical simulation and experimental validation[J].Ocean Engineering,2006,33:2135-2149.

[11]Akildiz H.A numerical study of the effects of the vertical baffle on liquid sloshing in two-dimensional rectangular container[J].Journal of Sound and Vibration,2012,332:41-52.

[12]Armenio V,Rocca M L.On the analysis of sloshing ofwater in rectangular containers:numerical study and experimental validation[J].Ocean Engineering,1996,23(8):705-739.

[13]Bae H S,Yun YW,Park M K.Flow response ofmagnetic fluid surface by pitchingmotion[J].Journal of Mechanical Science and Technology,2010,24(2):583-592.

[14]Celebi M S,Akyildiz H.Nonlinearmodeling of liquid sloshing in amoving rectangular container[J].Ocean Engineering, 2002,29:1527-1553.

[15]Xue M A,Lin P Z.Numerical study of ring baffle effects on reducing violent liquid sloshing[J].Computers and Fluids, 2011,52:116-129.

[16]肖龍飛,楊建民,郭 彬.淺水FPSO垂蕩和縱搖運(yùn)動(dòng)的低頻響應(yīng)[J].艦船科學(xué)技術(shù),2009,31(11):120-125. Xiao L F,Yang JM,Guo B.Low frequency heave and pitch motions of FPSO in shallow water[J].Ship Science and Technology,2009,31(11):120-125.

[17]尹立中,鄒經(jīng)湘,王本利.俯仰運(yùn)動(dòng)圓柱貯箱中液體的非線性晃動(dòng)[J].力學(xué)學(xué)報(bào),2000,32(3):280-290. Yin L Z,Zou JX,Wang B L.Nonlinear sloshing of liquid in a circular cylindrical container under pitching excitation[J]. Chinese Journal of Theoretical and Applied Mechanics,2000,32(3):280-290.

[18]洪 亮,朱仁傳,繆國(guó)平,范 菊.液艙晃蕩對(duì)船舶橫搖運(yùn)動(dòng)影響的數(shù)值研究[J].中國(guó)造船,2015,1:1-10. Hong L,Zhu R Z,Miu G P,Fan J.Numerical research on effectof sloshing on ship rollmotion[J].Shipbuilding of China, 2015,1:1-10.

[19]徐國(guó)徽,胡嘉駿,顧學(xué)康.耦合運(yùn)動(dòng)下液貨艙晃蕩壓力預(yù)報(bào)研究[J].船舶力學(xué),2016,20(7):884-891. Xu G H,Hu JJ,Gu X K.Prediction of sloshing pressure in coupledmotions for Liquid Cargo Carriers[J].Journal of Ship Mechanics,2016,20(7):884-891.

[20]李玉成,滕 斌.波浪對(duì)海上建筑物的作用(第2版)[M].北京:海洋出版社,2002. Li Y C,Teng B.Wave action onmarine structures(2nd Eidtion)[M].Beijing:Ocean Press,2002.

[21]寧德志.快速多極子邊界元方法在完全非線性水波問題中的應(yīng)用[D].大連:大連理工大學(xué),2005. Ning D Z.Application of the fastmultiple boundary elementmethod to the fully nonlinear water wave problems[D].Dalian: Ph.D thesis,Dalian University of Technology,2005.

[22]HiguchiM,Tanaka T,Endo S.Study on hull vibration-induced container liquid sloshing in LPG containerers[R].Nippon Kokan Technology Report,1976,72:111-122.

Nonlinear numerical simulation of liquid sloshing in a container subjected to pitch excitation

NING DE-zhi1,SONGWei-hua1,2,TENG Bin1
(1.State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China; 2.CCCC-FHDIEngineering Co.,Ltd.,Guangzhou 510230,China)

In the field ofmarine engineering,liquid sloshing is a kind of universal physical phenomenon. For the ship,the rotation has amore important influence than translation.Therefore,a time-domain numerical model is developed by higher-order boundary elementmethod to solve the liquid sloshing in a tank subjected to pitchmotion,in which the fully nonlinear boundary conditions satisfied on the free surface.With the coordinate exchanged between the global and local coordinates,the computational domain can be governed only in the local coordinates.In the solving process,a semi-mixed Eulerian-Lagrangian technique is applied to track the transient free surface and the 4th-order Runge-Kuttamethod is used to refresh wave elevation and velocity potential on the free surface at each time-step.The proposed numerical modelwas testified by comparison with the other published experimental and numerical results.On the base of validation,lots of numerical experiments are carried out to investigate the effect of pitch motion frequency,rotational center and a vertical bafflemounted on the container bottom on the sloshing surface and loads on the container.

liquid sloshing;pitchmotion;fully nonlinear;higher order boundary element;baffle

O353.2

：Adoi:10.3969/j.issn.1007-7294.2017.01.003

2016-09-01

國(guó)家自然科學(xué)基金資助項(xiàng)目（51679036，51490672）；教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃（NCET-13-0076）

寧德志（1975-），男，教授，博士生導(dǎo)師，E-mail：dzning@dlut.edu.cn；

宋偉華（1987-），男，工程師，E-mail：songweihua9@163.com。

1007-7294（2017）01-0015-08