朱倩云，梁玉珂，沈林維

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

0 引 言

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

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

1 研究流程及船舶模型

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

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

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

船舶在海上航行時，受到各種作用力和力矩。基于牛頓第二定律建立橫搖運(yùn)動方程式[4]：

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

表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+對船舶運(yùn)動進(jìn)行數(shù)值模擬。下面對船舶橫搖數(shù)值模擬的物理建模過程進(jìn)行說明。

（1）物理模型

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

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

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

（2）計(jì)算域

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

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

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

（4）時間步長

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

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

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

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

船舶勢能變換為

圖2 橫搖衰減曲線對比圖Fig.2 Comparison of roll attenuation curves obtained by test and simulation

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

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

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

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

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

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

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

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

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

2.2 橫搖模擬

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

2.2.1 模擬工況

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

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

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

2.2.2 結(jié)果分析

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

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

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

（1）頻率

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

（2）航速

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

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

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

（3）浪向

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

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

如圖7（a）所示，在同樣遭遇頻率下，艏斜浪中橫搖幅值較橫浪情況下大幅減小。穩(wěn)定橫搖時，不同浪向下船舶所受總力矩、壓力力矩和剪切力力矩幅值對比分別如圖7（b）、圖7（c）和圖7（d）所示。同一遭遇頻率下，斜浪中船體所受橫搖壓力力矩較橫浪中?。ㄒ妶D7（c））；而斜浪中由于航速和波浪傳播速度的疊加，自由表面水相對船體流速增大，因此橫搖中船體所受摩擦力矩減?。ㄈ鐖D7（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ù)船模遭遇波浪頻率對應(yīng)的實(shí)船遭遇波浪頻率為0.6 rad/s、0.62 rad/s、0.64 rad/s和0.66 rad/s。繪制實(shí)船橫搖RAO2曲線如圖8所示。

3 海浪譜

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

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

數(shù)學(xué)建模技術(shù)的發(fā)展為海浪自相關(guān)函數(shù)的研究提供了便利，Yang 等[10]借助Maple 軟件推導(dǎo)出PM波譜對應(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ù)是一個可以表示大部分特殊函數(shù)的通用函數(shù)，其表達(dá)式[11]定義如下：

式中，0≤m≤q，0≤n≤p，且對于任意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的波浪時歷，時間步長為0.04 s，數(shù)據(jù)長度為20 000。

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

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

圖9 海浪時歷曲線Fig.9 Time history of wave height

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

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

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

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

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

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

設(shè)置γ1、γ2和γ3三個修正系數(shù)對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ù)曲線與目標(biāo)海域?qū)嶋H自相關(guān)函數(shù)曲線對比如圖12所示，三參數(shù)修正得到的自相關(guān)函數(shù)曲線與目標(biāo)海域海浪時歷得到的實(shí)際自相關(guān)函數(shù)曲線擬合良好。

圖11 兩參數(shù)修正自相關(guān)函數(shù)與實(shí)際自相關(guān)函數(shù)對比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 中的修正波譜曲線。

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

圖13 海浪譜Fig.13 Ocean wave spectrum

4 橫搖譜密度

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

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

5 結(jié) 論

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

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

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

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