船舶參數(shù)橫搖非線性力學特征數(shù)值分析

楊 波，王 驍，吳 明

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

0 引 言

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

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

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

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

1 CFD數(shù)學模型

1.1 流體控制方程

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

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

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

1.2 波浪數(shù)學模型

采用船舶耐波性水池實驗中常用的微幅波模型模擬波浪。假設船舶靜止，建立以船舶重心為原點，X軸正向指向船尾，Y軸正向指向船舶右舷，Z軸正向垂直向上的船體坐標系，基于相對運動原理，可得波浪波高方程為：

速度方程為：

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

2 數(shù)值模擬方案

2.1 船模選擇

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

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

2.2 實驗參數(shù)

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

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

表2實驗參數(shù)Tab.2 Parameters of simulation

2.3 計算域劃分及網格生成

（1）計算域劃分

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

（2） 網格生成

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

圖2計算域Fig.2 Computational domain

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

2.4 邊界條件設置

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

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

出口邊界—壓力出口，設置靜水壓力；

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

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

2.5 其他設置

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

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

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

3.1 數(shù)值模擬結果

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

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

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

從表中可以看出：

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

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

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

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

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

從圖中可以看出：

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

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

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

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

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

從圖中可以看出：

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

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

圖6復原力矩時歷曲線Fig.6 Time history of restoring moment

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

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

從圖中可以看出：

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

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

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

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

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

（4）C至D時刻，船模繼續(xù)下沉，船體首傾繼續(xù)增大，同時波峰向船首移動，水面上升，導致船首部排水體積迅速增加，同時橫搖角繼續(xù)增大，使得復原力矩迅速增大。D點對應首傾最大值點，復原力矩達到最大值。

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

4 結 論

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

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

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

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

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

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