船舶非線性橫搖運(yùn)動(dòng)數(shù)值模擬研究

楊波，石愛國(guó)，王驍，吳明

（海軍大連艦艇學(xué)院航海系，遼寧大連116018）

船舶非線性橫搖運(yùn)動(dòng)數(shù)值模擬研究

楊波，石愛國(guó)，王驍，吳明

（海軍大連艦艇學(xué)院航海系，遼寧大連116018）

非線性大角橫搖運(yùn)動(dòng)是危害船舶安全的重要因素，也是船舶水動(dòng)力學(xué)關(guān)注的熱點(diǎn)及難點(diǎn)問(wèn)題。文章采用CFD方法，對(duì)某型驅(qū)逐艦船模在不同波長(zhǎng)規(guī)則波中橫浪航行時(shí)的橫搖運(yùn)動(dòng)進(jìn)行了數(shù)值模擬。通過(guò)與水池實(shí)驗(yàn)比對(duì)，數(shù)值仿真得到的穩(wěn)定橫搖幅值與實(shí)驗(yàn)值的誤差在10%左右。同時(shí)，數(shù)值模擬證實(shí)了非線性橫搖的兩種典型特征—“多值”及“跳躍”現(xiàn)象。數(shù)值模擬結(jié)果表明，CFD方法可作為研究船舶非線性橫搖運(yùn)動(dòng)的有效手段。

船舶；耐波性；非線性橫搖；CFD

0 引言

橫搖運(yùn)動(dòng)與船舶航行安全密切相關(guān)，往往是導(dǎo)致船舶傾覆的重要原因，因此一直是船舶耐波性研究所關(guān)心的問(wèn)題。但由于橫搖運(yùn)動(dòng)的復(fù)雜性，目前尚無(wú)準(zhǔn)確的理論預(yù)報(bào)方法。工程應(yīng)用中，一般采用基于實(shí)驗(yàn)獲得的經(jīng)驗(yàn)及半經(jīng)驗(yàn)公式對(duì)計(jì)算結(jié)果進(jìn)行修正[1]。

近年來(lái)，計(jì)算流體力學(xué)方法（Computational Fluid Dynamics，簡(jiǎn)稱CFD）在船舶水動(dòng)力研究中顯示出了巨大潛力，并逐步應(yīng)用于船舶的耐波性研究，而船舶的橫搖運(yùn)動(dòng)是其中的熱點(diǎn)之一，國(guó)內(nèi)外許多學(xué)者進(jìn)行了相關(guān)研究：Robert和Pablo等[2]對(duì)DTMB5512船模（帶/不帶舭龍骨）的強(qiáng)迫橫搖及橫搖衰減運(yùn)動(dòng)進(jìn)行了數(shù)值模擬，并進(jìn)行了詳細(xì)的不確定度分析，分析結(jié)果表明數(shù)值模擬達(dá)到了較高精度；國(guó)內(nèi)的張懷新[3]通過(guò)求解N-S方程，對(duì)Series 60船模不同位置的二維剖面的橫搖流場(chǎng)進(jìn)行模擬，得到了不同剖面的漩渦形狀并計(jì)算了相應(yīng)的橫搖阻尼；黃昊等[4]基于FLUENT軟件，對(duì)Series60船模的橫搖二維剖面繞流進(jìn)行了數(shù)值模擬，計(jì)算了橫搖阻尼；朱仁傳等[5]基于FLUENT軟件，對(duì)實(shí)尺度的集裝箱船S175的二維剖面橫搖繞流進(jìn)行數(shù)值模擬，通過(guò)數(shù)值模擬結(jié)果計(jì)算得了附加質(zhì)量和阻尼系數(shù)；紀(jì)東方和朱良生[6]利用FLUENT軟件對(duì)某一無(wú)航速駁船的二維剖面在靜水中的自由橫搖衰減運(yùn)動(dòng)進(jìn)行了數(shù)值模擬，計(jì)算了橫搖固有周期和阻尼系數(shù)；黃常青等[7]利用FLUENT軟件對(duì)某型起重船的強(qiáng)迫橫搖運(yùn)動(dòng)進(jìn)行了數(shù)值模擬，并計(jì)算了相應(yīng)的阻尼系數(shù)。

上述研究多限于船舶在靜水中的橫搖衰減及強(qiáng)迫橫搖運(yùn)動(dòng)模擬，且國(guó)內(nèi)的研究多限于二維或簡(jiǎn)單船型，對(duì)復(fù)雜線型船模在波浪中的橫搖運(yùn)動(dòng)，尤其是非線性大角橫搖運(yùn)動(dòng)尚無(wú)涉及。本文基于FLUENT平臺(tái)，對(duì)其進(jìn)行二次開發(fā)，首次對(duì)驅(qū)逐艦船模橫浪航行時(shí)的橫搖運(yùn)動(dòng)進(jìn)行了數(shù)值模擬。實(shí)驗(yàn)涉及多種波長(zhǎng)，最大穩(wěn)定橫搖幅值超過(guò)30度，橫搖運(yùn)動(dòng)呈現(xiàn)出明顯的非線性特征。數(shù)值模擬結(jié)果與水池實(shí)驗(yàn)結(jié)果吻合良好，證實(shí)了本文方法在船舶非線性橫搖預(yù)報(bào)方面的有效性及精度。

1 非線性橫搖運(yùn)動(dòng)方程

艦船的在規(guī)則波中的橫浪橫搖運(yùn)動(dòng)方程可表示為下述形式：

本文阻尼力矩采用如下非線性形式：

式中：B1和B2為系數(shù)，B2=2Nφφ（2Nφφ為小角度橫搖時(shí)的線性阻尼系數(shù)）。

復(fù)原力矩項(xiàng)采用如下形式：

式中：C1、C3、C5為系數(shù)，其中C1=Dh。

為得到解析解，對(duì)阻尼力矩及復(fù)原力矩進(jìn)行線性化。阻尼力矩線性化采用能量相等法。設(shè)有阻尼系數(shù)與組成的線性阻尼力矩若要用線性阻尼力矩代換非線性阻尼力矩需要保證在相同時(shí)間內(nèi)，兩者耗散的能量相同（做功相同），即：

設(shè)非線性方程的解為：

將（5）式代入（4）式，經(jīng)整理可得線性化的橫搖阻尼系數(shù)：

復(fù)原力矩線性化采取的方法是對(duì)非線性項(xiàng)展開，略去部分高階項(xiàng)。仍然采用（5）式表示的橫搖解析解，復(fù)原力矩可表示為：

將（6）式、（7）式代入（1）式，并除以（Iφφ+ΔIφφ）得：

式中：2νe為無(wú)因次化橫搖阻尼系數(shù)，ne為橫搖固有頻率，nφ為線性橫搖時(shí)的固有頻率。（8）式即為阻尼及復(fù)原力矩非線性時(shí)的橫搖運(yùn)動(dòng)方程。

將（5）式代入（8）式并進(jìn)行迭代求解，可得方程的解析解：

由（9）式可得穩(wěn)定橫搖幅值φa隨波頻ω的變化規(guī)律，即橫搖幅頻響應(yīng)曲線。

2 控制方程

對(duì)船舶運(yùn)動(dòng)來(lái)說(shuō)，可將水視作不可壓縮粘性流體，其控制方程主要包括連續(xù)性方程、動(dòng)量方程（NS方程）、湍流方程，各方程張量表示如下：

式中：xi、xj為笛卡爾坐標(biāo)方向，i、j=1，2，3；ui、ui′為xi方向速度時(shí)均值、脈動(dòng)值；P為壓力；t為時(shí)間；μ為運(yùn)動(dòng)粘性系數(shù)；k為湍動(dòng)能；ε為湍流耗散率

采用流體體積法（Volume of Fluid，簡(jiǎn)稱VOF）模擬自由面，體積分?jǐn)?shù)Cq的控制方程為：

式中：mpq為從q相傳輸?shù)絧相的質(zhì)量；ρq為第q相密度；t為時(shí)間。

3 數(shù)值模擬方案

3.1 船模選擇及實(shí)驗(yàn)參數(shù)

文獻(xiàn)[8]對(duì)某驅(qū)逐艦船模在規(guī)則波中橫浪航行時(shí)的系列橫搖運(yùn)動(dòng)進(jìn)行了水池實(shí)驗(yàn)研究，并采用ITTC的推薦方法對(duì)水池實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了不確定度分析，具有較高的可信度。為驗(yàn)證數(shù)值模擬方法精度，本文采用與文獻(xiàn)中相同的船模及實(shí)驗(yàn)參數(shù)。該型驅(qū)逐艦型船模與實(shí)船比例為1：50，其參數(shù)如表1所示，型線圖如圖1所示。共進(jìn)行20組實(shí)驗(yàn)，波頻范圍3.394≤ω≤5.459 rad/s，波高波長(zhǎng)比H/A=1/30，實(shí)驗(yàn)參數(shù)如表2所示。

表1 船模參數(shù)Tab.1 Parameters of ship model

表2 實(shí)驗(yàn)參數(shù)Tab.2 Parameters of simulations

圖1 船模型線圖Fig.1 Body plan of the ship model

圖2 計(jì)算域尺寸（頂視圖）Fig.2 Scale of computational domain(top view)

3.2 計(jì)算域及網(wǎng)格劃分

計(jì)算域設(shè)置為長(zhǎng)方體，在船體后部及右側(cè)（假設(shè)左舷來(lái)波）設(shè)有消波區(qū)。計(jì)算域尺寸為3L×5L×1.2L（長(zhǎng)×寬×深，L為船長(zhǎng)），其中水深為1L，具體尺寸如圖2所示。對(duì)計(jì)算域進(jìn)行分區(qū)處理，如圖3所示，各分區(qū)網(wǎng)格劃分情況為：

（1）船體近流場(chǎng)區(qū)域（Z1）：為船體表面到一個(gè)包裹船體的長(zhǎng)方體之間的區(qū)域，長(zhǎng)方體尺寸為1.5L×1.5B×1.5B（L為船長(zhǎng)，B為船寬）；該區(qū)域生成非結(jié)構(gòu)性網(wǎng)格，在船體表面布設(shè)20層邊界層，船體面網(wǎng)格尺寸為0.01 m，該區(qū)域網(wǎng)格最大尺寸為0.02 m。

（2）與近流場(chǎng)區(qū)域等高的圓柱形區(qū)域（Z2），將Z1“包裹”在其內(nèi)部，圓柱底面半徑為2B，軸線為船舶橫搖軸；該區(qū)域生成結(jié)構(gòu)性網(wǎng)格，內(nèi)邊界網(wǎng)格尺寸及類型與Z1相同，外邊界（底面及側(cè)面）網(wǎng)格最大尺寸約為0.035 m。

（3）外流場(chǎng)區(qū)域（Z3）為除Z1、Z2之外的計(jì)算域。該區(qū)域生成結(jié)構(gòu)性網(wǎng)格，在自由面處加密Z方向網(wǎng)格，以生成高質(zhì)量波浪環(huán)境，Z網(wǎng)格尺寸為A/2（A為波幅），其余網(wǎng)格由自由面向上下漸疏。

設(shè)置船體僅作橫搖運(yùn)動(dòng)，通過(guò)動(dòng)網(wǎng)格技術(shù)實(shí)現(xiàn)。將Z1及Z2設(shè)置為剛體，隨船體轉(zhuǎn)動(dòng)，Z3相對(duì)于坐標(biāo)系保持靜止，Z1、Z2與Z3的內(nèi)部邊界設(shè)置為interface，船模運(yùn)動(dòng)通過(guò)船體所受水動(dòng)力驅(qū)動(dòng)。

圖3 計(jì)算域分區(qū)處理Fig.3 Multi-blocks of computational domain

圖4 邊界名稱Fig.4 Boundary conditions

3.3 邊界條件

各邊界條件名稱如圖4所示，具體設(shè)置如下（假設(shè)左舷來(lái)波）：

（1）入口邊界（in，port）—速度入口，給定三個(gè)方向的波速及水的體積分?jǐn)?shù)；

（2）出口邊界（out，stab）—壓力出口，設(shè)置靜水壓力；

（3）外邊界（top，bot）—速度入口，給定三個(gè)方向流速（u=v=w=0）及水的體積分?jǐn)?shù)（0）；

（4）船體（hull）—壁面，有剪切力無(wú)滑移。

4 數(shù)值模擬結(jié)果分析

本文數(shù)值模擬得到的各波頻的穩(wěn)定橫搖幅值與水池實(shí)驗(yàn)值的對(duì)比如表3所示（CFD為仿真值，EXP為水池實(shí)驗(yàn)值）。從表中可以看出，仿真值與水池實(shí)驗(yàn)值的誤差均在10%以內(nèi)。

表3 橫搖幅值Tab.3 Roll amplitude

續(xù)表3

4.1 “偏頂”現(xiàn)象

圖5中虛線為船模的固有頻率曲線，CFD（三角）為本文數(shù)值仿真值，THE（實(shí)線）為（9）式的理論解。從圖中可以看出：

（1）本文數(shù)值模擬值隨頻率變化規(guī)律與理論解較為接近，且所捕捉到的橫搖最大值及多值區(qū)所對(duì)應(yīng)的頻率范圍，亦與理論解相似。這間接證明了本文數(shù)值模擬的可信度。

（2）幅頻曲線分布于固有頻率曲線兩側(cè)，且產(chǎn)生了非常明顯的“偏頂”現(xiàn)象，固有頻率曲線基本通過(guò)幅頻曲線最大值點(diǎn)。

圖5 幅頻曲線Fig.5 Roll amplitude-frequency curve

4.2 多值及跳躍現(xiàn)象

從表3可以看出，當(dāng)波浪頻率ω=3.642和3.672時(shí)，對(duì)應(yīng)2個(gè)穩(wěn)定橫搖幅值。圖6為圖5的局部放大圖，從圖中可以清楚地觀察到多值區(qū)域。

圖6 多值區(qū)域Fig.6 Frequency domain of multi-amplitude

雖然存在多值區(qū)域，但是其頻率區(qū)域非常小，大約在3.60～3.70 rad/s范圍內(nèi)；同時(shí)，在實(shí)驗(yàn)中發(fā)現(xiàn)，雖然在兩個(gè)頻率點(diǎn)產(chǎn)生了多值現(xiàn)象，但是其實(shí)現(xiàn)過(guò)程并不相同。下面通過(guò)仿真時(shí)歷進(jìn)行說(shuō)明。

圖7 橫搖角時(shí)歷圖（ω=3.642 rad/s）Fig.7 Time history of roll angle(ω=3.642 rad/s)

（1）圖7為ω=3.642 rad/s的橫搖仿真時(shí)歷圖。在仿真開始時(shí)刻，通過(guò)已經(jīng)生成好的橫浪流場(chǎng)初始化橫搖運(yùn)動(dòng)流場(chǎng)，船模很快進(jìn)入穩(wěn)定橫搖狀態(tài)，橫搖幅值穩(wěn)定在17.29°左右（見表3）。在t=42 s左右，對(duì)船加一瞬時(shí)激勵(lì)，此激勵(lì)幅值為穩(wěn)定橫搖力矩的10倍，作用時(shí)間為T/20（T為波浪周期）。這樣船模會(huì)產(chǎn)生一個(gè)非常大的激勵(lì)幅值，約為50°左右，大約經(jīng)3～4個(gè)周期后，船模即進(jìn)入另一穩(wěn)定橫搖狀態(tài)，橫搖幅值約為36°（見表3）。但是無(wú)論再對(duì)船模再次進(jìn)行何種激勵(lì)，船模也無(wú)法回到初始時(shí)的穩(wěn)定橫搖狀態(tài)。

（2）圖8為ω=3.672 rad/s的橫搖仿真時(shí)歷圖，與圖7采用相同的初始化方法。經(jīng)過(guò)小幅振蕩后，橫搖角幅值穩(wěn)定在37°左右。采用與圖7相同的激勵(lì)方法，橫搖角振蕩后仍然穩(wěn)定在37°左右，也就是說(shuō)，在此種工況下難以出現(xiàn)多值現(xiàn)象。

圖8 橫搖角時(shí)歷圖（ω=3.672 rad/s）Fig.8 Time history of roll angle(ω=3.672 rad/s)

（3）圖9同為ω=3.672 rad/s的橫搖仿真時(shí)歷圖，采用與圖7相同的初始化方法。但是在計(jì)算開始階段，約束船模不發(fā)生橫搖運(yùn)動(dòng)；當(dāng)t=5 s左右，去掉橫搖約束，船模開始運(yùn)動(dòng)，在經(jīng)過(guò)4～5個(gè)周期后，橫搖幅值慢慢趨于穩(wěn)定，此時(shí)橫搖幅值為20°左右；穩(wěn)定一段時(shí)間后加激勵(lì)，此時(shí)船模很快穩(wěn)定于另一橫搖幅值，約為37°左右，出現(xiàn)多值現(xiàn)象。

圖9 橫搖角時(shí)歷圖（ω=3.672 rad/s）Fig.9 Time history of roll angle(ω=3.672 rad/s)

從上述分析可以看出：

（1）雖然此型船模存在橫搖多值區(qū)域，但是由于靜穩(wěn)性曲線及阻尼力矩的特點(diǎn)，產(chǎn)生多值的頻段很窄，較難發(fā)生“多值”及“跳躍”現(xiàn)象。

（2）如果船舶在某一波頻上會(huì)產(chǎn)生多值現(xiàn)象，初始較小橫搖角的穩(wěn)定狀態(tài)很容易被外部激勵(lì)改變，而穩(wěn)定于較大的橫搖角幅值；反之，較大橫搖角的穩(wěn)定狀態(tài)則要穩(wěn)定得多，不易變化。

（3）多值現(xiàn)象的出現(xiàn)及其狀態(tài)特點(diǎn)，與船模的初始條件和擾動(dòng)力矩特點(diǎn)密切相關(guān)，即使處于理論多值頻段內(nèi)，多值現(xiàn)象也不一定會(huì)出現(xiàn)。

5 結(jié)論

本文采用CFD方法，對(duì)某型驅(qū)逐艦船模在不同波長(zhǎng)規(guī)則波中橫浪航行時(shí)的橫搖運(yùn)動(dòng)進(jìn)行了數(shù)值模擬。數(shù)值模擬結(jié)果與水池實(shí)驗(yàn)結(jié)果吻合良好，穩(wěn)定橫搖幅值的誤差在10%以內(nèi)，證實(shí)了橫搖幅頻曲線的“偏頂”、“多值”及“跳躍”現(xiàn)象。上述模擬結(jié)果說(shuō)明，采用CFD方法研究船舶的非線性橫搖運(yùn)動(dòng)是可行的，且具有較高精度，可以作為理論研究及水池實(shí)驗(yàn)研究的有力補(bǔ)充。

[1]Chakrabarti S.Empirical calculation of roll damping for ships and barges[J].Ocean Engineering,2001,28:915-932.

[2]Wilson R V,Carrica P M,Stern F.Unsteady RANS method for ship motions with application to roll for a surface combatant[J].Computers&Fluids,2006,35:501-524.

[3]張懷新,劉應(yīng)中,繆國(guó)平.船體各種剖面的橫搖阻尼與漩渦的形狀[J].水動(dòng)力學(xué)研究與進(jìn)展,2001,16(3):382-389. Zhang Huaixin,Liu Yingzhong,Miao Guo-ping.Vortex patterns and roll damping at various cross sections of ship[J]. Journal of Hydrodynamics,2001,16(3):382-389.

[4]黃昊,郭海強(qiáng),朱仁傳,繆國(guó)平.粘性流中船舶橫搖阻尼計(jì)算[J].船舶力學(xué),2008,12(4):578-573. Huang Hao,Guo Haiqiang,Zhu Renchuan,Miao Guoping.Computations of ship roll damping in viscous flow[J].Journal of Ship Mechanics,2008,12(4):578-57.

[5]朱仁傳,郭海強(qiáng),繆國(guó)平,余建偉.一種基于CFD理論船舶附加質(zhì)量與阻尼的計(jì)算方法[J].上海交通大學(xué)學(xué)報(bào),2009, 43(2):198-203. Zhu Renchuan,Guo Haiqiang,Miao Guoping,Yu Jianwei.A computational method for evaluation of added mass and damping of ship based on CFD theory[J].Journal of Shanghai Jiaotong University,2009,43(2):198-203.

[6]紀(jì)東方,朱良生.粘性流中船舶自由橫搖衰減運(yùn)動(dòng)數(shù)值模擬[J].科學(xué)技術(shù)與工程,2009,9(23):7061-7065. Ji Dongfang,Zhu Liangsheng.Numerical simulation of ship free roll decay motion in viscous flow[J].Science Technology and Engineering,2009,9(23):7061-7065.

[7]黃常青,王學(xué)林,胡于進(jìn).基于CFD的起重船水動(dòng)力系數(shù)數(shù)值模擬[J].中國(guó)機(jī)械工程,2011,22(17):2076-2079. Huang Changqing,Wang Xuelin,Hu Yujin.Numerical simulation of hydronamic coefficients of crane ship based on CFD [J].China Mechanical Engineering,2011,22(17):2076-2079.

[8]Alberto Francescutto,Giorgio Contento.Bifurcations in ship rolling:Experimental results and parameter identification technique[J].Ocean Engineering,1999,26:1095-1123.

Simulation of ship’s non-linear roll motion

YANG Bo,SHI Ai-guo,WANG Xiao,WU Ming
(Dept.of Navigation,Dalian Naval Academy,Dalian 116018,China)

Non-linear roll motion is a crucial factor in endangering ship’s security and a hot issue in ship’s hydrodynamics study.In this paper,ship’s roll motion in different regular beam sea waves is simulated. The simulation results show good agreement with the tank test results,and the error is less than 10%.The typical phenomena of nonlinear roll-‘multi-amplitude’and‘jumping’are also realized in the simulation. The simulation results show that CFD method can be an effective tool in studying ship’s non-linear motion.

ship;seakeeping;non-linear roll;CFD(Computational Fluid Dynamics)

U661.74

A

10.3969/j.issn.1007-7294.2015.05.005

1007-7294（2015）05-0509-09

2014-12-15

遼寧省博士科研啟動(dòng)基金項(xiàng)目（20111037）

楊波（1983-），男，博士研究生，E-mail：bisecond@163.com；

石愛國(guó)（1956-），男，教授，博士生導(dǎo)師。