朱倩云，梁玉珂，沈林維

（浙江大學(xué)海洋學(xué)院，浙江舟山 316021）

0 引 言

對(duì)實(shí)際海浪中的船舶耐波性進(jìn)行分析，需要對(duì)不同環(huán)境因素下的船舶非線(xiàn)性橫搖運(yùn)動(dòng)響應(yīng)進(jìn)行計(jì)算。目前，雖然模型試驗(yàn)在確定船舶橫搖方面依舊處于核心地位，但隨著數(shù)值模擬技術(shù)的發(fā)展及有利的發(fā)展前景，需要更多的數(shù)值模擬結(jié)果來(lái)建立對(duì)該問(wèn)題的一般方法[1-2]。為了更好地探索船舶的橫搖運(yùn)動(dòng)，本文基于CFD 軟件Star-CCM+，利用波浪激勵(lì)船舶橫搖的方式來(lái)模擬船舶橫搖運(yùn)動(dòng)，并計(jì)算不同浪向、航速下，遭遇頻率大于固有頻率的頻段中船舶橫搖幅值響應(yīng)算子。

實(shí)際海洋環(huán)境中的波浪為不規(guī)則波，目前尚不能從理論上直接得到海浪頻譜。包括Star-CCM+在內(nèi)的CFD 計(jì)算軟件,對(duì)不規(guī)則波的模擬仍是基于傳統(tǒng)的海浪譜展式，難以反映海浪情況，因此，探索海浪譜密度的估算方法成為必要。P-M 譜是描述充分發(fā)展海洋的典型波譜[3]，本文基于P-M 譜的自相關(guān)函數(shù)解析式展開(kāi)研究，以期根據(jù)目標(biāo)海域?qū)崪y(cè)數(shù)據(jù)修正海浪譜展式，為探索海浪時(shí)歷信號(hào)的自相關(guān)函數(shù)與海洋環(huán)境參數(shù)間的關(guān)系打下基礎(chǔ)。

1 研究流程及船舶模型

本文一方面基于CFD 方法建立數(shù)值模型，分析波浪、航速及其耦合對(duì)船舶橫搖運(yùn)動(dòng)的影響，計(jì)算得到船舶橫搖幅值響應(yīng)算子；另一方面基于成熟海浪譜展式，推導(dǎo)得到海浪自相關(guān)函數(shù)表達(dá)式，結(jié)合目標(biāo)海域時(shí)歷數(shù)據(jù)對(duì)其進(jìn)行修正，反推目標(biāo)海域海浪譜；最后，結(jié)合上面兩部分內(nèi)容對(duì)船舶的橫搖譜密度進(jìn)行預(yù)報(bào)。具體流程如圖1所示。

本文選取標(biāo)模DTMB 5512 作為研究對(duì)象，模型與實(shí)船的縮尺比為1:46.6，實(shí)船與模型的主尺度如表1所示。

采用船體坐標(biāo)系分析船舶運(yùn)動(dòng)姿態(tài)，設(shè)船舶重心為原點(diǎn)（O），原點(diǎn)指向船艏為x軸正方向，原點(diǎn)指向右舷為y軸正方向，原點(diǎn)指向船底為z軸正方向。

船舶在海上航行時(shí)，受到各種作用力和力矩?；谂ｎD第二定律建立橫搖運(yùn)動(dòng)方程式[4]：

式中，Ixx+A44為橫搖質(zhì)量慣性矩系數(shù)和附加質(zhì)量慣性矩系數(shù)，N44=B44+B44v為線(xiàn)性阻尼力矩系數(shù)和非線(xiàn)性粘性阻尼力矩系數(shù)之和，S44為復(fù)原力矩系數(shù)，ME4(t)為橫搖激勵(lì)力矩。

表1 DTMB 5512船模與實(shí)船主尺度Tab.1 Principal dimensions of ship model DTMB 5512

2 橫搖幅值響應(yīng)算子計(jì)算

2.1 數(shù)值模型

在參考文獻(xiàn)[14]的基礎(chǔ)上，利用Star-CCM+對(duì)船舶運(yùn)動(dòng)進(jìn)行數(shù)值模擬。下面對(duì)船舶橫搖數(shù)值模擬的物理建模過(guò)程進(jìn)行說(shuō)明。

（1）物理模型

求解基于雷諾平均的連續(xù)方程和動(dòng)量守恒方程[5]：

式中，u=(u,v,w)代表時(shí)均速度場(chǎng)，u、v、w是速度矢量在x、y和z方向的分量，ρ是模擬中水或空氣的密度，p*代表時(shí)均壓力，f代表質(zhì)量力，μν為動(dòng)力粘度，τRe為雷諾應(yīng)力。利用SSTk-ω湍流模型來(lái)計(jì)算獲得雷諾應(yīng)力。

船舶運(yùn)動(dòng)過(guò)程中，采用二階VOF（volume of fluid）方法來(lái)捕捉自由液面。同時(shí)，采用二階迎風(fēng)格式離散對(duì)流項(xiàng)，并選取隱式方法，利用分離流模型求解非耦合條件下的控制方程。

（2）計(jì)算域

根據(jù)ITTC 船舶CFD 應(yīng)用指南[6]，建立一個(gè)長(zhǎng)方體計(jì)算域。計(jì)算域包括流入面和流出面，流入面施加一個(gè)已知的速度場(chǎng)；流出面設(shè)置為壓力出口邊界，并以阻尼消波方式抑制回流；將兩側(cè)、頂部和底部邊界設(shè)為速度入口來(lái)避免壁面和流動(dòng)間產(chǎn)生速度梯度。上述邊界條件設(shè)置能將所有橫向邊界上的流動(dòng)都指向出口邊界，保證邊界上的流動(dòng)反射速率最小。分別將流入面、流出面、兩側(cè)邊界、頂部邊界和底部邊界設(shè)置在離船重心2倍船長(zhǎng)處、4倍船長(zhǎng)處、2.5倍船長(zhǎng)處、1.5倍船長(zhǎng)處和2.5倍船長(zhǎng)處。

（3）網(wǎng)格劃分

為捕捉船舶大幅橫搖運(yùn)動(dòng)，采用重疊網(wǎng)格法，包括船體域和背景域網(wǎng)格劃分。針對(duì)船型特征和橫搖研究特點(diǎn)，首先對(duì)船體控制區(qū)域進(jìn)行整體網(wǎng)格加密，再對(duì)球鼻艏、舭龍骨等區(qū)域進(jìn)行局部加密；針對(duì)規(guī)則波問(wèn)題，自由液面應(yīng)保證每個(gè)波長(zhǎng)內(nèi)至少包含40 個(gè)單元格、每個(gè)波高范圍包含20 個(gè)單元格[7]，為準(zhǔn)確捕捉液面波動(dòng)，本文將自由液面網(wǎng)格加密區(qū)的垂向距離擴(kuò)充至兩倍波高。共計(jì)網(wǎng)格608.11 萬(wàn)個(gè)（背景網(wǎng)格165.55萬(wàn)個(gè)，重疊區(qū)域網(wǎng)格442.56萬(wàn)個(gè)）。

（4）時(shí)間步長(zhǎng)

本文考慮ITTC 推薦的時(shí)間步長(zhǎng)和庫(kù)朗數(shù)來(lái)確定時(shí)間步長(zhǎng)。船模的橫搖固有周期為1.54 s，根據(jù)ITTC推薦得時(shí)間步長(zhǎng)Δt＜0.0154 s；為保持?jǐn)?shù)值穩(wěn)定性，要求滿(mǎn)足網(wǎng)格單元的庫(kù)朗數(shù)CFL小于1，自由表面的庫(kù)朗數(shù)小于0.5[1]，波浪中的橫搖計(jì)算以自由表面最小網(wǎng)格尺寸來(lái)確定時(shí)間步長(zhǎng)，即Δt＜0.003 125 s。本文取時(shí)間步長(zhǎng)為0.003 s。

（5）數(shù)值驗(yàn)證

除文獻(xiàn)[13]提供的數(shù)值驗(yàn)證之外，本文對(duì)初始強(qiáng)迫橫搖角為10°、傅汝德數(shù)Fr=0.41的船模衰減橫搖進(jìn)行了仿真，以進(jìn)一步驗(yàn)證數(shù)值模型。模擬結(jié)果與試驗(yàn)結(jié)果對(duì)比如圖2 所示。

船舶在靜水中衰減，半個(gè)周期內(nèi)橫搖幅值由?i變?yōu)?i+1，由于橫搖在到達(dá)峰值時(shí)橫搖角速度皆為零，根據(jù)能量守恒定律可知?jiǎng)菽茏儞Q等于阻尼耗散能量。

船舶勢(shì)能變換為

圖2 橫搖衰減曲線(xiàn)對(duì)比圖Fig.2 Comparison of roll attenuation curves obtained by test and simulation

式中，D為船舶重量，-- --GM為初穩(wěn)性高

將阻尼力矩展開(kāi)為三次項(xiàng)形式，有

則阻尼做功，即消耗的能量為

用弧度制繪制消滅曲線(xiàn)并進(jìn)行三次多項(xiàng)式擬合：

得到系數(shù)a、b、c。

將式（7）代入式（4），根據(jù)ΔE=W可得

可求得等效線(xiàn)性阻尼系數(shù)：

擬合曲線(xiàn)峰值得到的阻尼系數(shù)如表2 所示。CFD 模擬結(jié)果中，一次和二次阻尼系數(shù)與試驗(yàn)結(jié)果較吻合，三次阻尼系數(shù)相差較多；橫搖幅值不大時(shí)，三次阻尼系數(shù)對(duì)阻尼力矩影響較小。穩(wěn)態(tài)周期橫搖時(shí)，線(xiàn)性阻尼系數(shù)N?1和平方阻尼系數(shù)N?2可視為常數(shù)；但大幅橫搖時(shí)，存在船舶舭龍骨出水的情況，此時(shí)N?1、N?2并非常數(shù)，Bassler[8]將其視為分段函數(shù)，傳統(tǒng)的利用橫搖衰減技術(shù)得到的橫搖衰減曲線(xiàn)計(jì)算等效線(xiàn)性阻尼的方法并不適用。

表2 阻尼系數(shù)對(duì)比Tab.2 Comparison of damping coefficients obtained by test and simulation

2.2 橫搖模擬

本節(jié)基于CFD 方法考慮不同自由度、航速、遭遇頻率、浪向等對(duì)橫搖運(yùn)動(dòng)特性的影響。為保持結(jié)果的一致性，船舶初始橫搖角設(shè)置為0°，受波浪激勵(lì)發(fā)生橫搖運(yùn)動(dòng)。

2.2.1 模擬工況

如表3 所示，選取12 種工況研究波浪頻率、浪向和航速對(duì)橫搖運(yùn)動(dòng)特性的影響，橫搖激勵(lì)采用了一階線(xiàn)性規(guī)則波激勵(lì)，波高為0.06 m，監(jiān)測(cè)了Fr=0.246（中等航速）和Fr=0.41（高航速）兩種航速狀態(tài)，以及中等航速下270°和225°浪向角（浪向角為船舶前進(jìn)方向與波浪傳播方向的夾角，且以波浪傳播方向順時(shí)針度量，如圖3所示）。在模型尺度下，選取諧搖區(qū)激勵(lì)頻率作為研究對(duì)象，研究遭遇頻率大于船模固有頻率的頻段，遭遇頻率計(jì)算式如下：

表3 橫搖模擬工況Tab.3 Simulation of rolling conditions

式中，ω為波浪頻率，U為船速，βwave為浪向角。

2.2.2 結(jié)果分析

為了便于分析，一致截取物理模擬時(shí)間前20 s進(jìn)行對(duì)比。圖4 顯示了不同航速、波浪頻率下船舶橫搖時(shí)歷變化曲線(xiàn)?？梢园l(fā)現(xiàn)，經(jīng)過(guò)幾個(gè)周期后，橫搖運(yùn)動(dòng)趨于穩(wěn)定。

圖4 不同激勵(lì)頻率、航速下船舶橫搖時(shí)歷對(duì)比圖（橫浪）Fig.4 Comparison of ship roll histories at different excitation frequencies and speeds (transverse wave)

圖3 浪向示意圖Fig.3 Schematic of wave direction

（1）頻率

不同航速下橫搖幅值隨頻率的變化情況見(jiàn)圖5?？梢园l(fā)現(xiàn)在不同頻率的橫浪條件下，船模在頻率為4.04 rad/s 下會(huì)產(chǎn)生較大的橫搖運(yùn)動(dòng)特性幅值，此時(shí)波浪頻率最接近模型的固有頻率；隨著波浪頻率遠(yuǎn)離固有頻率，穩(wěn)定周期內(nèi)橫搖幅值越來(lái)越小。

（2）航速

本文根據(jù)簡(jiǎn)諧激勵(lì)技術(shù)估算不同航速橫搖阻尼[9]。如圖6 所示的橫浪工況中，在4.19～4.49 rad/s 頻段內(nèi)，航速較高的情況下橫搖阻尼更大；高航速條件下，本文研究頻段與其固有頻率相差較大，同一頻率下橫搖幅值較中等航速中橫搖幅值??；中等航速工況中，橫搖阻尼變化也更快；Fr=0.246 時(shí)，4.04 rad/s 的遭遇頻率已接近船固有頻率，產(chǎn)生較大橫搖角，導(dǎo)致船舶舭龍骨出水，再入水時(shí)受到自由表面沖擊，同時(shí)伴有氣泡入水，橫搖阻尼進(jìn)一步增加。

圖6 無(wú)因次阻尼系數(shù)B44Fig.6 Roll damping coefficients B44

圖5 不同航速下船舶最大橫搖角隨激勵(lì)頻率變化Fig.5 Variation of maximum roll angle with wave frequency at different speeds

（3）浪向

本文以225°浪向角的艏斜浪下、Fr=0.246工況為例，對(duì)與橫浪工況中相同遭遇頻率的船舶橫搖響應(yīng)值進(jìn)行求解分析。

圖7 不同浪向下橫搖特性隨遭遇頻率變化Fig.7 Variation of rolling characteristics in different wave direction conditions with the encountered frequency

如圖7（a）所示，在同樣遭遇頻率下，艏斜浪中橫搖幅值較橫浪情況下大幅減小。穩(wěn)定橫搖時(shí)，不同浪向下船舶所受總力矩、壓力力矩和剪切力力矩幅值對(duì)比分別如圖7（b）、圖7（c）和圖7（d）所示。同一遭遇頻率下，斜浪中船體所受橫搖壓力力矩較橫浪中?。ㄒ?jiàn)圖7（c））；而斜浪中由于航速和波浪傳播速度的疊加，自由表面水相對(duì)船體流速增大，因此橫搖中船體所受摩擦力矩減小（如圖7（d）所示）。

2.3 結(jié)果統(tǒng)計(jì)

橫搖頻率響應(yīng)函數(shù)（RAO）可以表示為

式中，ζA為波幅。

模型與實(shí)船存在以下關(guān)系：

圖8 橫搖RAO2Fig.8 Roll-RAO2

式中，下標(biāo)‘s’和‘m’分別表示實(shí)船和模型。從而可得實(shí)船的RAO如下：

統(tǒng)計(jì)12種工況下橫搖響應(yīng)幅值計(jì)算結(jié)果，根據(jù)船模遭遇波浪頻率對(duì)應(yīng)的實(shí)船遭遇波浪頻率為0.6 rad/s、0.62 rad/s、0.64 rad/s和0.66 rad/s。繪制實(shí)船橫搖RAO2曲線(xiàn)如圖8所示。

3 海浪譜

通過(guò)對(duì)P-M 波譜對(duì)應(yīng)的自相關(guān)函數(shù)的解析式修正得到目標(biāo)海域海浪譜展式，有利于對(duì)特定海域的波浪進(jìn)行有效表達(dá)，可為目標(biāo)海域中的船舶橫搖譜密度預(yù)報(bào)提供新思路。

3.1 P-M譜自相關(guān)函數(shù)

數(shù)學(xué)建模技術(shù)的發(fā)展為海浪自相關(guān)函數(shù)的研究提供了便利，Yang 等[10]借助Maple 軟件推導(dǎo)出PM波譜對(duì)應(yīng)的海浪自相關(guān)函數(shù)表達(dá)式為

式中：α=8.1×10-3；β=0.74；g為重力加速度，單位m/s2；V為海面上19.5 m高度處的風(fēng)速，單位m/s；S(ω)的單位為m2·s；Meijer-G函數(shù)是一個(gè)可以表示大部分特殊函數(shù)的通用函數(shù)，其表達(dá)式[11]定義如下：

式中，0≤m≤q，0≤n≤p，且對(duì)于任意j和k，Γ(bj-s)的極點(diǎn)不得與Γ(1 -ak+s)的極點(diǎn)重合（其中，j=1,…,m;k=1,…,n）。

3.2 目標(biāo)海域自相關(guān)函數(shù)

如圖9所示，選取某實(shí)驗(yàn)室不規(guī)則波仿真程序模擬得到的有義波高為3 m、跨零周期為8.16 s的波浪時(shí)歷，時(shí)間步長(zhǎng)為0.04 s，數(shù)據(jù)長(zhǎng)度為20 000。

海浪譜估計(jì)時(shí)，常將海浪視為平穩(wěn)隨機(jī)過(guò)程，即海浪信號(hào)ηt的期望μη和標(biāo)準(zhǔn)差σ不隨時(shí)間t變化，其自相關(guān)函數(shù)可表示為時(shí)間延遲τ的函數(shù)：

利用Matlab求得其自相關(guān)函數(shù)，如圖10中紅色實(shí)線(xiàn)所示。圖10中藍(lán)色虛線(xiàn)為相應(yīng)風(fēng)速下由P-M譜自相關(guān)函數(shù)表達(dá)式所繪制的曲線(xiàn)。

圖9 海浪時(shí)歷曲線(xiàn)Fig.9 Time history of wave height

圖10 P-M波譜對(duì)應(yīng)自相關(guān)函數(shù)與實(shí)際海浪自相關(guān)函數(shù)對(duì)比Fig.10 Comparison of autocorrelation function obtained from P-M spectrum with actual wave autocorrelation function

兩條曲線(xiàn)變換趨勢(shì)基本一致，可以通過(guò)設(shè)置修正系數(shù)來(lái)修正由P-M波譜得到的自相關(guān)函數(shù)。

（1）兩參數(shù)修正

觀察由P-M 波譜得到的自相關(guān)函數(shù)曲線(xiàn)與目標(biāo)海浪實(shí)際相關(guān)函數(shù)曲線(xiàn)，設(shè)置γ1和γ2兩個(gè)修正系數(shù)對(duì)P-M波譜得到的自相關(guān)函數(shù)進(jìn)行初步修正，修正自相關(guān)函數(shù)表達(dá)式如下：

令γ1=0.708，γ2=2.05，修正后的自相關(guān)函數(shù)曲線(xiàn)與目標(biāo)海域?qū)嶋H自相關(guān)函數(shù)曲線(xiàn)的對(duì)比如圖11 所示。兩參數(shù)修正后的自相關(guān)函數(shù)曲線(xiàn)與目標(biāo)海域海浪時(shí)歷得到的自相關(guān)函數(shù)曲線(xiàn)相對(duì)吻合但還有一定差距，擬納入第三個(gè)參數(shù)γ3對(duì)由P-M波譜得到的自相關(guān)函數(shù)進(jìn)行修正。

（2）三參數(shù)修正

設(shè)置γ1、γ2和γ3三個(gè)修正系數(shù)對(duì)P-M 波譜得到的自相關(guān)函數(shù)進(jìn)行進(jìn)一步修正，修正自相關(guān)函數(shù)表達(dá)式如下：

令γ1=0.595，γ2=1.77，γ3=0.5，三參數(shù)修正后的自相關(guān)函數(shù)曲線(xiàn)與目標(biāo)海域?qū)嶋H自相關(guān)函數(shù)曲線(xiàn)對(duì)比如圖12所示，三參數(shù)修正得到的自相關(guān)函數(shù)曲線(xiàn)與目標(biāo)海域海浪時(shí)歷得到的實(shí)際自相關(guān)函數(shù)曲線(xiàn)擬合良好。

圖11 兩參數(shù)修正自相關(guān)函數(shù)與實(shí)際自相關(guān)函數(shù)對(duì)比Fig.11 Comparison between the autocorrelation function obtained by modifying the two parameters and the actual wave autocorrelation function

3.3 目標(biāo)海域海浪譜

由維納-辛欽定理可知，已知目標(biāo)海域海浪自相關(guān)函數(shù)，易得海浪譜密度函數(shù)。由式（17）推得海浪譜密度函數(shù)：

由式（18）推得海浪譜密度函數(shù)圖像如圖13 中的修正波譜曲線(xiàn)。

通過(guò)對(duì)圖9 的海浪時(shí)歷數(shù)據(jù)直接進(jìn)行譜分析來(lái)比較驗(yàn)證圖13 中的修正波譜。為減小譜分析誤差，本文使用周期圖譜估計(jì)，對(duì)數(shù)據(jù)分塊處理[12]，最后將估計(jì)得到的海浪譜密度函數(shù)進(jìn)行樣條擬合得到圖13 中的紅色虛線(xiàn)。對(duì)比圖13中的三參數(shù)修正波譜與譜估計(jì)擬合曲線(xiàn)，可以發(fā)現(xiàn)海浪譜峰值十分接近，且頻率分布接近。

圖13 海浪譜Fig.13 Ocean wave spectrum

4 橫搖譜密度

對(duì)于波浪中船舶的橫搖運(yùn)動(dòng)問(wèn)題，將船舶視為線(xiàn)性系統(tǒng)，則海浪譜為輸入，橫搖響應(yīng)幅值算子RAO為傳遞函數(shù)，橫搖譜密度為輸出且等于海浪譜和響應(yīng)幅值算子的乘積，即

結(jié)合修正海浪譜，繪制橫搖譜密度，如圖14所示。

5 結(jié) 論

圖14 橫搖譜估計(jì)Fig.14 Roll power spectral density

本文針對(duì)DTMB 5512 船模的耐波性問(wèn)題，著重研究難以獨(dú)立預(yù)報(bào)橫搖阻尼的情況中船舶的非線(xiàn)性橫搖運(yùn)動(dòng)表現(xiàn)，探究了各因素對(duì)船舶橫搖幅值的影響；同時(shí)基于P-M波譜公式推導(dǎo)出的海浪自相關(guān)函數(shù)表達(dá)式和海浪數(shù)據(jù)獲取了目標(biāo)海域海浪譜，進(jìn)而預(yù)報(bào)船舶橫搖譜密度。結(jié)論總結(jié)如下：

（1）通過(guò)對(duì)CFD 仿真結(jié)果分析發(fā)現(xiàn)，相同頻率的橫浪激勵(lì)下，增大航速會(huì)增大橫搖阻尼，并增大船體兩端壓差從而抑制橫搖，船舶橫搖幅值減??；隨著遭遇頻率遠(yuǎn)離固有頻率，船舶橫搖幅值越來(lái)越小，受渦旋的產(chǎn)生、舭龍骨出水等因素的影響，呈非線(xiàn)性變化；受波浪與航速的耦合作用影響，相同遭遇頻率下，225°浪向角的艏斜浪中船舶的橫搖幅值小于橫浪中的船舶橫搖幅值。

（2）結(jié)合目標(biāo)海域的離散海浪數(shù)據(jù)對(duì)基于P-M 波譜公式推出的海浪自相關(guān)函數(shù)進(jìn)行修正獲取目標(biāo)海域海浪譜。研究表明，基于Meijer-G 函數(shù)可以對(duì)海浪自相關(guān)函數(shù)進(jìn)行有效表達(dá)，利用海浪自相關(guān)函數(shù)推導(dǎo)海浪譜有助于解決無(wú)法用光滑的海浪譜對(duì)其海浪能量分布進(jìn)行精準(zhǔn)表達(dá)的問(wèn)題，可以靶向服務(wù)于特定目標(biāo)海域中航行船舶橫搖譜密度的預(yù)報(bào)。