楊 波，王 驍，吳 明

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

0 引 言

當(dāng)船舶縱浪航行時(shí)，復(fù)原力矩會(huì)隨船體與波浪相對(duì)位置的改變而發(fā)生變化，對(duì)于具有較大外飄船首和方尾的船型，這種變化尤為明顯。此時(shí)很小的初始橫向擾動(dòng)也有可能導(dǎo)致大幅橫搖運(yùn)動(dòng)，這種現(xiàn)象稱為參數(shù)橫搖。一般認(rèn)為，當(dāng)海浪波長與船長近似相等、船舶與海浪的遭遇頻率約為橫搖固有頻率的2倍時(shí)，最易發(fā)生參數(shù)橫搖。目前國際海事組織（IMO）正在推進(jìn)的“第二代完整穩(wěn)性衡準(zhǔn)”中，將參數(shù)橫搖作為船舶波浪中5種穩(wěn)性失效模式之一，是目前船舶水動(dòng)力學(xué)領(lǐng)域研究的熱點(diǎn)。

參數(shù)橫搖的研究方法主要有模型試驗(yàn)和理論研究兩種。理論研究方面，Pauling和Rosenberg（1959）[1]最早采用單自由度馬休方程對(duì)參數(shù)橫搖進(jìn)行求解，并假定波浪中初穩(wěn)性高按弦值函數(shù)變化，但該方法對(duì)參數(shù)橫搖的預(yù)報(bào)并不準(zhǔn)確。隨著研究不斷深入，研究者開始注意到橫搖阻尼、非線性復(fù)原力矩和其它自由度運(yùn)動(dòng)對(duì)參數(shù)橫搖的影響，并開展了相關(guān)研究。對(duì)于橫搖阻尼，常用方法是基于橫搖衰減實(shí)驗(yàn)數(shù)據(jù)，將其表示為橫搖角速度的高階函數(shù)[2-3]。對(duì)于非線性復(fù)原力矩，則采用橫搖角的高階多項(xiàng)式表達(dá)，然后對(duì)與波浪相關(guān)的時(shí)變項(xiàng)進(jìn)行修正[4-5]。除此之外，許多學(xué)者還利用基于勢流理論的切片法[6-7]和三維面元法[8-9]計(jì)算時(shí)變復(fù)原力矩。Spanos等（2009）[10]對(duì)14種時(shí)域勢流方法進(jìn)行了比較研究，結(jié)果表明三維面元法較切片法能更好地模擬參數(shù)橫搖，但所有方法對(duì)大幅運(yùn)動(dòng)的波浪水動(dòng)力計(jì)算均較差。對(duì)于其他自由度的影響，Bulian（2005）[11]提出了一種1.5自由度模型，以研究垂蕩和縱搖對(duì)參數(shù)橫搖的影響；更多的學(xué)者通過橫搖-縱搖-垂蕩三自由度模型[12-13]以及近年來發(fā)展起來的操縱/耐波六自由度統(tǒng)一模型[18-20]，對(duì)參數(shù)橫搖運(yùn)動(dòng)進(jìn)行研究。魯江等[14-15]對(duì)隨機(jī)波浪中的參數(shù)橫搖進(jìn)行了研究，驗(yàn)證了隨機(jī)波浪下參數(shù)橫搖的非各態(tài)歷經(jīng)特點(diǎn)。

可以看出，理論研究中，對(duì)于參數(shù)橫搖影響較大的復(fù)原力矩和非線性阻尼，均采用近似方法進(jìn)行計(jì)算，盡管做了若干修正，但在船舶大幅橫搖時(shí)仍然存在誤差。模型實(shí)驗(yàn)是研究參數(shù)橫搖的可靠方法，但較為費(fèi)時(shí)費(fèi)力，且不易開展單項(xiàng)影響因素的深入研究。目前，計(jì)算流體力學(xué)方法（CFD）已廣泛運(yùn)用于船舶水動(dòng)力計(jì)算，其基于粘性流理論，相較于勢流理論有著天然優(yōu)勢，可精確計(jì)算船舶在波浪環(huán)境中的水動(dòng)力和力矩，這為基于力學(xué)方法研究參數(shù)橫搖提供了可能。Hamid等（2010）[23]采用CFD方法模擬了一艘水面艦艇的迎浪參數(shù)橫搖運(yùn)動(dòng)，與水池實(shí)驗(yàn)結(jié)果吻合較好，但是國內(nèi)基于CFD方法的參數(shù)橫搖研究相對(duì)較少。

復(fù)原力矩和橫搖阻尼是影響參數(shù)橫搖模擬的重要因素。作者曾基于CFD方法對(duì)船模的靜水橫搖衰減運(yùn)動(dòng)進(jìn)行數(shù)值模擬研究[24]，數(shù)值模擬結(jié)果與水池實(shí)驗(yàn)吻合良好，驗(yàn)證了CFD方法用于數(shù)值計(jì)算橫搖復(fù)原力矩和阻尼的有效性。本文在前述研究的基礎(chǔ)上，基于CFD方法研究了DTMB5512船模頂浪航行時(shí)的參數(shù)橫搖運(yùn)動(dòng)。實(shí)驗(yàn)涵蓋多種航速和波高，分析了波浪遭遇周期和波高對(duì)參數(shù)橫搖的影響?；趯?shí)驗(yàn)結(jié)果，對(duì)激勵(lì)船模發(fā)生參數(shù)橫搖時(shí)復(fù)原力矩的非線性特征進(jìn)行了分析，并對(duì)產(chǎn)生非線性特征的原因進(jìn)行了探討。本文方法可對(duì)瘦削船型頂浪航行時(shí)的參數(shù)橫搖進(jìn)行預(yù)報(bào)，為參數(shù)橫搖的力學(xué)特征分析提供了新的方法，可用于船模的波浪中完整穩(wěn)性的定量評(píng)估。

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

1.1 流體控制方程

船舶水動(dòng)力研究可將水視為不可壓縮粘性流體，控制方程包括連續(xù)性方程和動(dòng)量方程（N-S方程），其張量表示為：

式中：ui、uj分別為流體速度矢量?u在xi、xj方向的分量，t為時(shí)間，P為壓力，ρ為流體密度，fi為質(zhì)量力，μ為流體動(dòng)力粘性系數(shù)。

采用RNGk-ε模型模擬湍流，使上述方程組封閉，湍動(dòng)能k及湍流耗散率ε的控制方程的張量表示為：

1.2 波浪數(shù)學(xué)模型

采用船舶耐波性水池實(shí)驗(yàn)中常用的微幅波模型模擬波浪。假設(shè)船舶靜止，建立以船舶重心為原點(diǎn)，X軸正向指向船尾，Y軸正向指向船舶右舷，Z軸正向垂直向上的船體坐標(biāo)系，基于相對(duì)運(yùn)動(dòng)原理，可得波浪波高方程為：

速度方程為：

式中：u、v、w分別為x、y、z三個(gè)方向的速度，H為波高，k為波數(shù)，ω 為波浪圓頻率，V為航速，ε0為初始相位。

2 數(shù)值模擬方案

2.1 船模選擇

選擇國際拖曳水池會(huì)議（International Towing Tank Conference，ITTC）推薦船型DTMB5512船模為研究對(duì)象，其型線圖如圖1所示，主要船型參數(shù)如表1所示。從型線圖中可以看出，該船型具有外飄船首和方形尾，是易發(fā)生參數(shù)橫搖的船型。

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

2.2 實(shí)驗(yàn)參數(shù)

參數(shù)橫搖影響因素較多，遭遇周期（Te）和波高（H）是比較重要的兩項(xiàng)，且當(dāng)波長（λ）一定時(shí)，遭遇周期只受航速（V）影響。選擇最易發(fā)生參數(shù)橫搖的波長（λ=L，L為船長，下同），波高/波長比（H/λ）為 0.04，進(jìn)行多種付汝德數(shù)（Fr）下數(shù)值模擬，具體實(shí)驗(yàn)參數(shù)如表2所示。

圖1船模型線圖Fig.1 Ship model’s body plan

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

2.3 計(jì)算域劃分及網(wǎng)格生成

（1）計(jì)算域劃分

計(jì)算域?yàn)殚L方體，長、寬、深為 4L×2L×1.7L（L為船長），其中入口距船首1L，船尾距消波區(qū)1L，另有長1L的消波區(qū)，水深1L，自由面距上邊界0.7L，船模與水池相對(duì)位置如圖2所示。

（2） 網(wǎng)格生成

采用混合網(wǎng)格，在船體周圍區(qū)域采用非結(jié)構(gòu)網(wǎng)格以較好地表現(xiàn)船型，其余區(qū)域采用結(jié)構(gòu)網(wǎng)格以減少網(wǎng)格數(shù)量并提高計(jì)算效率，在自由面附近進(jìn)行網(wǎng)格加密以滿足造波要求。船體面網(wǎng)格尺寸為0.01 m，布置5層邊界層網(wǎng)格，第一層網(wǎng)格厚度保證y+<30，總網(wǎng)格數(shù)為255萬。艦首部分面網(wǎng)格及球鼻首附近邊界層網(wǎng)格如圖3所示。

圖2計(jì)算域Fig.2 Computational domain

圖3船首面網(wǎng)格及邊界層網(wǎng)格Fig.3 Mesh of ship bow and boundary layers

2.4 邊界條件設(shè)置

采用邊界造波法模擬波浪，計(jì)算域的邊界條件具體設(shè)置如下：

入口邊界—速度入口，按照（3）-（4）式給定波高及速度值；

出口邊界—壓力出口，設(shè)置靜水壓力；

船體—壁面，有剪切力無滑移；

外邊界（包含水池底部、頂部及側(cè)壁）—壁面，剪切力為0。

2.5 其他設(shè)置

采用VOF方法追蹤自由面，采用壁面函數(shù)法模擬近壁面流動(dòng)，VOF方程離散采用改進(jìn)的HRIC格式，其余方程采用二階迎風(fēng)格式，速度壓力耦合方程求解采用SIMPLEC算法。每個(gè)時(shí)間步內(nèi)迭代求解20次，當(dāng)主要物理量殘差小于10-4時(shí)，該時(shí)間步計(jì)算收斂，進(jìn)行下一時(shí)間步計(jì)算。

數(shù)值模擬時(shí)，將船模以初始橫傾角3°置于頂浪流場中，以提供初始橫搖擾動(dòng)；待流場穩(wěn)定后，使船模做橫搖/縱搖/垂蕩三自由度耦合運(yùn)動(dòng)。

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

3.1 數(shù)值模擬結(jié)果

圖4是未發(fā)生和發(fā)生參數(shù)橫搖時(shí)典型的橫搖角變化時(shí)歷，圖中t為時(shí)間，φ為橫搖角。

可以看出，由于初始橫傾角的存在，艦船會(huì)在復(fù)原力矩作用下發(fā)生橫搖運(yùn)動(dòng)，當(dāng)不發(fā)生參數(shù)橫搖時(shí)，橫搖角幅值逐漸變小為0；當(dāng)發(fā)生參數(shù)橫搖時(shí)，橫搖角幅值則會(huì)不斷增大，直至達(dá)到最大值并保持穩(wěn)定。

通過橫搖時(shí)歷曲線可得穩(wěn)定橫搖角幅值（φ0），不同航速下的φ0如表3所示。

從表中可以看出：

（1）航速對(duì)是否發(fā)生參數(shù)橫搖有重要影響。對(duì)于該型船來說，發(fā)生參數(shù)橫搖的航速主要集中于中高速段；

（2）發(fā)生參數(shù)橫搖的四個(gè)速度對(duì)應(yīng)的波浪遭遇周期分別為0.859、0.797、0.744和0.698 s（見表2）。可見，當(dāng)航速使遭遇周期Te接近Tφ/2（Tφ為橫搖周期，見表1）時(shí)，易發(fā)生參數(shù)橫搖；且愈接近，最大橫搖角幅值愈大。這與已有研究成果相符。

圖4未發(fā)生/發(fā)生參數(shù)橫搖時(shí)的橫搖時(shí)歷曲線Fig.4 Time history of non-/parametric rolling

表3不同航速時(shí)的穩(wěn)定橫搖幅值Tab.3 Steady rolling amplitudes at different speeds

圖5是發(fā)生參數(shù)橫搖時(shí)不同航速下橫搖時(shí)歷曲線，圖中t為時(shí)間，φ為橫搖角。

從圖中可以看出：

（1）當(dāng)Fr=0.25、0.30、0.35和0.40時(shí)，達(dá)到穩(wěn)定橫搖狀態(tài)分別用了22、12、15和27個(gè)橫搖周期，可見當(dāng)遭遇周期越接近Tφ/2時(shí)，達(dá)到穩(wěn)定橫搖所需時(shí)間越少。

（2） 經(jīng)計(jì)算，穩(wěn)定橫搖時(shí)，F(xiàn)r=0.25、0.30、0.35 和 0.40 對(duì)應(yīng)的橫搖周期分別為 1.728 s、1.596 s、1.486 s和1.400 s，基本為遭遇周期的2倍。

圖5不同航速參數(shù)橫搖時(shí)歷曲線Fig.5 Time history of parametric rolling at different speeds

3.2 參數(shù)橫搖非線性力學(xué)特征分析

導(dǎo)致參數(shù)橫搖發(fā)生的直接原因在于橫搖力矩的非線性變化。圖6為發(fā)生參數(shù)橫搖時(shí)橫搖復(fù)原力矩變化時(shí)歷，圖中t為時(shí)間，Mφ為橫搖復(fù)原力矩。其中6（a）為初始至穩(wěn)定橫搖整個(gè)過程，而6（b）為一個(gè)橫搖周期內(nèi)的復(fù)原力矩與橫搖角變化對(duì)比圖。

從圖中可以看出：

（1）復(fù)原力矩幅值從由小變大，直至穩(wěn)定橫搖時(shí)保持不變；

（2）復(fù)原力矩曲線會(huì)隨橫搖角增大而出現(xiàn)2個(gè)拐點(diǎn)，呈現(xiàn)出先增大、后減小、再急劇增大的非線性特征。尤其是在參數(shù)橫搖發(fā)展初始階段，這一減小過程尤為明顯，復(fù)原力矩值會(huì)接近于0，甚至?xí)禐樨?fù)值（力矩方向發(fā)生改變）。正是這一過程，使得艦船不能及時(shí)扶正，而在慣性作用下橫搖不斷增大；而減小過程結(jié)束后，復(fù)原力矩又迅速增大，使艦船在2/4周期中加速回?fù)u，以較大的角速度通過正浮位置，搖向另一舷。在后1/2周期中，又出現(xiàn)與前1/2周期相同的過程（區(qū)別在于力矩方向改變）。如此反復(fù)，促使橫搖角越來越大，產(chǎn)生參數(shù)橫搖。

圖6復(fù)原力矩時(shí)歷曲線Fig.6 Time history of restoring moment

3.3 橫搖力矩非線性特征原因分析

基于數(shù)值模擬結(jié)果，本文發(fā)現(xiàn)，參數(shù)橫搖復(fù)原力矩曲線的非線性與縱搖、垂蕩及船/波相對(duì)位置曲線呈現(xiàn)強(qiáng)相關(guān)性，其特征點(diǎn)（零點(diǎn)、2個(gè)拐點(diǎn)及最大值點(diǎn)）可以分別用垂蕩、縱搖和船/波相對(duì)位置表征。圖7為一個(gè)橫搖周期內(nèi)復(fù)原力矩、縱搖和垂蕩時(shí)歷對(duì)比圖，圖中t為時(shí)間，Mφ、θ、Z為復(fù)原力矩、縱搖和垂蕩值（為便于比較，圖中Z值乘以100）。圖8為船模周圍瞬時(shí)波高圖，圖8（a）、（b）、（c）、（d）分別對(duì)應(yīng)圖 7 中 A、B、C、D 點(diǎn)，可見B、C為拐點(diǎn)，D為復(fù)原力矩最大點(diǎn)。

從圖中可以看出：

（1）A時(shí)刻船模橫搖角為0，復(fù)原力矩接近于0，此時(shí)波峰距船首約1/4船長，船模位置上浮，船體稍尾傾。

（2）A至B時(shí)刻，波谷逐漸向船首移動(dòng)，導(dǎo)致水面降低，同時(shí)船模向上運(yùn)動(dòng)，并處于尾傾狀態(tài)，導(dǎo)致船首部排水體積迅速減小。但由于橫搖角增大，復(fù)原力矩仍然增大。B點(diǎn)對(duì)應(yīng)垂蕩最高點(diǎn)，此時(shí)船首排水體積達(dá)到最小值，復(fù)原力矩開始變小。

（3）B至C時(shí)刻，船模開始下沉，并由尾傾過渡到首傾，但由于波谷通過船體前部，波面繼續(xù)降低，導(dǎo)致船模前部排水體積繼續(xù)減小。此時(shí)雖然橫搖角在慣性作用下繼續(xù)增大，但復(fù)原力矩不斷減小。C點(diǎn)時(shí)波谷距艦首距離約1/4船長，此時(shí)復(fù)原力矩降至最小。

圖7復(fù)原力矩、縱搖和垂蕩時(shí)歷曲線Fig.7 Time history of restoring moment,pitch and heave

圖8不同時(shí)刻瞬時(shí)波高圖Fig.8 Transient wave height at different times

（4）C至D時(shí)刻，船模繼續(xù)下沉，船體首傾繼續(xù)增大，同時(shí)波峰向船首移動(dòng)，水面上升，導(dǎo)致船首部排水體積迅速增加，同時(shí)橫搖角繼續(xù)增大，使得復(fù)原力矩迅速增大。D點(diǎn)對(duì)應(yīng)首傾最大值點(diǎn)，復(fù)原力矩達(dá)到最大值。

（5）D時(shí)刻后，船模在復(fù)原力矩作用下回?fù)u，橫搖角迅速歸零，向另一舷運(yùn)動(dòng)。同時(shí)，從圖7中可以看出，縱搖及垂蕩周期等于波浪遭遇周期，而橫搖周期為2倍遭遇周期。這使得縱搖、垂蕩及船/波相對(duì)運(yùn)動(dòng)影響周期性作用于橫搖復(fù)原力矩，從而產(chǎn)生參數(shù)橫搖。

4 結(jié) 論

本文基于CFD方法，對(duì)某瘦長型船模在規(guī)則波中頂浪航行時(shí)的參數(shù)橫搖運(yùn)動(dòng)進(jìn)行了數(shù)值模擬。模擬結(jié)果表明，當(dāng)航速使波浪遭遇周期約接近于Tφ/2時(shí)，越容易發(fā)生參數(shù)橫搖，橫搖值也越大。這與現(xiàn)有研究成果是相符的，證明了CFD方法用于參數(shù)橫搖模擬的有效性。通過對(duì)參數(shù)橫搖時(shí)復(fù)原力矩的非線性特征和產(chǎn)生原因進(jìn)行分析，可得到以下結(jié)論：

（1）船體排水體積的變化是復(fù)原力矩變化的直接原因。對(duì)于具有較大外飄角的船型，船首部的排水體積變化對(duì)復(fù)原力矩變化起主要作用，而這一變化受縱搖、垂蕩和船/波相對(duì)位置共同作用。

（2）參數(shù)橫搖力矩曲線在橫搖角變大過程中會(huì)出現(xiàn)2個(gè)拐點(diǎn)，呈現(xiàn)出先增大、后減小、再增大的非線性特征。當(dāng)波長等于船長時(shí)，曲線的4個(gè)特征點(diǎn)（零點(diǎn)、2個(gè)拐點(diǎn)和最大值點(diǎn)）可以分別由波峰距船首1/4船長、垂蕩最大值點(diǎn)、波谷距船首1/4船長、首傾最大值點(diǎn)進(jìn)行表征。

（3）頂浪運(yùn)動(dòng)使船舶縱搖和垂蕩周期等于波浪遭遇周期，而使橫搖周期基本為2倍遭遇周期。這使得縱搖、垂蕩及船/波相對(duì)位置周期性地對(duì)復(fù)原力矩產(chǎn)生影響，從而導(dǎo)致參數(shù)橫搖。

（4）基于CFD方法可以對(duì)參數(shù)橫搖時(shí)的非線性復(fù)原力矩進(jìn)行準(zhǔn)確計(jì)算和分析，這是以往方法所難以實(shí)現(xiàn)的。

本文研究僅對(duì)波長等于船長及頂浪運(yùn)動(dòng)工況進(jìn)行了數(shù)值模擬，下一步將開展不同波長及不同航向工況的參數(shù)橫搖研究。