楊云濤，朱仁傳，陸 安，高 慧

（1.江蘇科技大學船舶與建筑工程學院，江蘇張家港 215600;2.上海交通大學船舶海洋與建筑工程學院海洋工程國家重點實驗室，上海 200240）

0 引 言

船舶在海洋環(huán)境中航行時，不可避免地會遭遇來自不同方向的波浪。此時，船舶除了以某一平均速度的前進運動外，在波浪的激勵下還會產(chǎn)生六自由度的搖蕩運動。嚴重的搖蕩運動不僅會影響船舶維持其正常使用功能的能力，在某些極端情況下甚至會造成船體傾覆或損毀。因此，開展船-波作用水動力問題研究，準確可靠地預報船舶在波浪中航行時的運動響應，對船舶的設計、使用以及航行安全等具有重要意義。

研究船舶這類大尺度物體在波浪作用下的運動響應時，流體粘性和船體曲率突變處的旋渦分離影響通常是局部、次要的，勢流模型是對該問題的一個合理簡化和近似。在過去的幾十年中，基于勢流理論的切片法[1]因其計算快捷、穩(wěn)定且易于實現(xiàn)等優(yōu)勢，在船舶運動響應計算中應用廣泛。但是該方法是在諸多假定的基礎上建立的，對于高速、低頻以及存在外飄的肥大船型，其數(shù)值預報結果的精度往往差強人意[2]。為了克服這些缺陷，近年來隨著計算機技術的快速發(fā)展，出現(xiàn)了大量關于三維勢流方法的研究工作，這些三維方法通?？梢苑譃樽杂擅娓窳趾瘮?shù)法和Rankine源法兩類。自由面格林函數(shù)法采用自動滿足以輻射和均流線性化自由面條件的格林函數(shù)為邊界積分方程的核函數(shù)，數(shù)值求解時只需要在物面布源，離散量較小，是目前解決零航速船-波作用問題的一種有效方法[3]。但是當船舶以一定的航速在波浪中航行時，該方法存在一定的局限性：由于有航速自由面格林函數(shù)存在高頻振蕩特性，因而數(shù)值計算復雜且對存在外飄的非直壁船型會出現(xiàn)計算不穩(wěn)定現(xiàn)象[4]；滿足的是均流線性化自由面條件，無法考慮船體定常繞流對非定常擾動的耦合影響；在某些特殊的頻率下，計算失效、無法得到可靠的數(shù)值解，出現(xiàn)所謂的不規(guī)則頻率現(xiàn)象[5]。

采用以Laplace 方程的基本解1/r為核函數(shù)的Rankine 源法可以克服自由面格林函數(shù)法的這些局限。根據(jù)對時間因子處理方式的不同，Rankine 源法可以分為時域和頻域法。由于時域Rankine 具有使用靈活、可采用只含速度勢一階導數(shù)的自由面的運動學和動力學條件等優(yōu)勢[6]，被廣泛應用于船舶在波浪中運動問題的研究。但從公開發(fā)表的文獻[7-9]來看，這些研究大多針對的是船舶迎浪航行時的垂蕩和縱搖運動。實際上，船舶在海上多為斜浪航行，伴隨著六個自由度的耦合搖蕩運動（包括縱蕩、橫蕩、垂蕩、橫搖、縱搖和首搖），在采用上述時域Rankine源法進行數(shù)值模擬時，存在縱蕩、橫蕩和首搖模態(tài)缺少回復力的問題[10]。為了解決這一問題，學者們[11]提出了引入人工彈簧模型的方法，但該方法的模擬精度仍有待進一步提高。相比時域，頻域模型不存在沒有回復力的問題，更適合斜浪中船舶運動響應計算，且計算效率也更高。但是，由于在Brard 數(shù)小于0.25 時難以處理船前產(chǎn)生波浪反射等原因[12]，目前關于頻域Rankine源法的研究大多局限于高航速或高頻的迎浪工況。

本文針對船舶在斜浪中航行時的水動力學問題，在頻域勢流理論框架下，基于高階Rankine 源法建立考慮定常繞流和橫搖粘性阻尼影響的搖蕩運動響應計算模型。數(shù)值實現(xiàn)時，在自由面條件中引入Rayleigh 人工阻尼，并采用二階迎風差分格式進行離散處理，以保證輻射條件得到滿足；在物面條件采用計及定常擾動勢的mj項，并利用積分法進行求解，以提高數(shù)值計算的精度；采用CFD 方法模擬船舶的自由橫搖衰減運動，并根據(jù)能量法計算粘性阻尼系數(shù)，代入運動方程求解運動響應，以避免勢流理論無法考慮流體粘性所造成的橫搖運動計算誤差。在此基礎上，采用Fortran自主開發(fā)數(shù)值程序，對不同船型在規(guī)則波中以不同遭遇浪向航行時的六自由度運動響應進行數(shù)值計算，并對其結果進行分析。

1 斜浪中航行船舶運動的頻域勢流模型

1.1 坐標系統(tǒng)和搖蕩運動定義

設船舶以平均速度U和遭遇浪向β在波幅為A、圓頻率為ω0、波數(shù)為k的深水規(guī)則波中斜浪航行（存在色散關系2=gk，其中g為重力加速度），建立圖1 所示的兩組右手坐標系：坐標系o-xyz為以船速U移動、原點位于靜水面、x軸指向船舶前進方向、z軸豎直向上并通過平衡狀態(tài)船舶重心的參考坐標系；坐標系ob-xbybzb為固定于船體、在平衡狀態(tài)時與參考坐標系重合的動坐標系。在參考坐標系中，船上實際感受到的入射波頻率變?yōu)樵庥鲱l率ωe=ω0-kUcosβ，船舶在波浪激勵下產(chǎn)生的搖蕩運動可以分解為重心沿三個坐標軸的直線運動（縱蕩X1、橫蕩X2和垂蕩X3）以及動坐標繞三個坐標軸的轉動（橫搖X4、縱搖X5和首搖X6）。

圖1 坐標系統(tǒng)和搖蕩運動Fig.1 Coordinate systems and ship motions

1.2 速度勢及其邊值問題

研究船舶在波浪中的運動問題時，流體（即水）可以看作是均勻、不可壓縮且無粘的理想流體，流動可以認為是無旋的。此時，流場中任意一點的速度可以通過引入勢函數(shù)Φ(x,y,z,t)來表示。在頻域框架下，速度勢Φ可以分解為

式中：右端第一項ˉ?（x，y，z）為船舶定常運動的擾動勢；剩余三項為非定常擾動勢，其中?I和?7分別為非定常入射勢和繞射勢的空間部分，?j（j=1，2，…，6）為船舶在j方向以單位速度搖蕩產(chǎn)生的規(guī)范化輻射勢。

對于定常擾動勢ˉ?（x，y，z），在一階近似范圍內可采用疊模流勢代替[13]，它在流域內滿足Laplace方程

在線性波理論中，?I的表達式為

至于速度勢?j（j=1，2，…，7），以往為了簡便起見，通常對其進行求解時假定定常擾動勢ˉ?為小量，忽略船舶定常繞流的影響。但這一假定只適用于船體細長、前進速度是低速的情況。為了提高計算精度，本文在自由面和物面條件中計及ˉ?的影響，采用以下定解條件：

式中：r=(x,y,z)為位置矢量。

1.3 非定常流體作用力和六自由度運動方程

船舶在波浪中航行時受到的外力按是否隨時間變化可以分為定常力和非定常力。其中定常力主要包括重力、靜浮力、推進力以及定常阻力，它們之間互相平衡，不會引起船舶的搖蕩運動。因此討論船舶在波浪中的運動問題時，只需要分析作用在船體上的非定常力。

若假定流體是無旋的理想流體，非定常流體作用力可通過利用Bernoulli方程和式（1）的速度勢分解計算壓強并沿船體濕表面積分求得，它可以分為入射力FIi、繞射力FDi、輻射力FRi和靜恢復力FSi四個部分：

上述理想流體（無粘）假定在多數(shù)情況都是合理的，但是當船舶存在橫搖運動時，由于粘性對其影響較大，還需考慮粘性力的作用。通常粘性力是非線性的，但可以采用如下等效線性化的形式表示：

2 數(shù)值實現(xiàn)

2.1 基于高階Rankine源法計算繞輻射勢

由于繞輻射勢?j滿足以Laplace 方程為控制方程的邊值問題（6），根據(jù)格林第三公式，它可以采用如下的邊界積分方程進行求解：

為了克服傳統(tǒng)常值單元離散所存在的變量在邊界上不連續(xù)、離散量大（計算效率低）以及難以獲得精確的物面導數(shù)等缺陷，本文采用圖2所示的9 節(jié)點二次單元對船體表面進行離散。經(jīng)過上述的離散化處理，并將方程（6）中的自由面和物面條件代入式（14），可得如下的離散邊界積分方程組：

圖2 9節(jié)點二次單元Fig.2 9-node bi-quadratic element

式中：NB和NF分別代表離散船體表面和自由面的單元數(shù)；Nk（s，t）和J（s，t）分別為形函數(shù)和Jacobian 矩陣，它們的表達式見文獻[17]。

求解上述方程組的關鍵在于mj項的計算和?[ˉ?,?j]中?j的二階空間導數(shù)的處理。根據(jù)公式（7）可知，mj取決于定常擾動勢ˉ?的一階和二階梯度，其中一階梯度易于計算（在采用源分布法確定源強后，對定常擾動勢的邊界積分方程求梯度即可獲得ˉ?的一階梯度），難點在于二階梯度的求解。鑒于傳統(tǒng)的直接數(shù)值差分法數(shù)值實現(xiàn)困難且精度不高[18]，本文采用間接的積分法對ˉ?的二階梯度進行計算[14]。該方法將ˉ?沿三個坐標軸的偏導數(shù)看作未知變量，建立關于它們的控制方程和邊界條件，然后采用Rankine源法構建邊界積分方程：

至于?[ˉ?,?j]中?j的二階空間導數(shù)，為了反映流體相對船舶自上游流向下游（波動在上游衰減快，在下游衰減慢）的特征，保證數(shù)值計算的穩(wěn)定性和精確性，本文采用如下的二階迎風差分格式：

式中：i為自由面上離散單元節(jié)點的編號，它的數(shù)值從下游往上游逐漸增大；Ci、Ci+1、Ci+2和Ci+3為迎風差分格式的系數(shù)，它們的表達式見文獻[19]。

2.2 基于CFD法計算橫搖阻尼系數(shù)

當船舶斜浪航行時，會產(chǎn)生橫搖運動。橫搖受流體粘性影響較大，在采用頻域運動方程（13）計算它的數(shù)值時，除了需要計及興波阻尼系數(shù)外，還需要補充額外的橫搖粘性阻尼系數(shù)。而后者是無法基于前文所述的頻域勢流理論求得的。目前，模型試驗是估算橫搖阻尼最為可靠的方法，但該方法存在成本高、速度慢以及換算到實船存在尺度效應等缺陷。隨著計算機技術以及數(shù)值算法的快速發(fā)展，CFD 方法已成為研究船舶水動力問題的一種重要數(shù)值手段。CFD 方法相比模型試驗更為經(jīng)濟、高效，相比勢流方法可以考慮流體粘性。鑒于此，本文將采用CFD方法，通過仿試驗的方式獲得船舶自由橫搖衰減曲線，在此基礎上通過數(shù)值處理獲得橫搖阻尼系數(shù)。

CFD 方法模擬船舶自由橫搖衰減運動所采用的控制方程、初始和邊界條件、計算域范圍、網(wǎng)格尺寸以及時間步等參數(shù)按照本文作者公開發(fā)表的文獻[20]中的建議來設置。為了避免直接處理數(shù)值結果所產(chǎn)生的不穩(wěn)定性，采用如下的衰減函數(shù)[21]擬合橫搖衰減曲線：

式中，AD為橫搖衰減曲線的最大橫搖幅值，βi為衰減系數(shù)項，ωD為橫搖衰減運動固有頻率，ε為相位角。

根據(jù)能量守恒定律，船舶在自由橫搖衰減過程中，ti到ti+1時間段內總能量的變化應等于阻尼力矩所做的功。由此可以推得[22]

在基于CFD法確定了橫搖阻尼系數(shù)之后，本文進一步根據(jù)上一節(jié)的高階Rankine源法開發(fā)數(shù)值程序計算繞輻射勢，求解水動力系數(shù)和波浪力，并將其代入頻域運動方程（13）便可求得船舶的六自由度運動響應。

3 數(shù)值方法驗證

為了驗證本文數(shù)值方法和程序的可靠性，以ITTC 標準船模S175 集裝箱船為對象，開展船舶在首斜浪（β=150°）中垂蕩、縱搖和橫搖運動響應預報研究。S175 的船長L=175 m，船寬B=25.4 m，吃水D=9.5 m，排水量Δ=247 42 t，縱搖慣性半徑kyy=0.24L，橫搖慣性半徑kxx=0.328B。

圖3 S175自由橫搖衰減運動計算網(wǎng)格Fig.3 Meshes used for the numerical simulation of roll delay motion of S175

圖4 S175自由橫搖衰減曲線Fig.4 Free roll decay curve of S175

將采用CFD 法確定的橫搖阻尼系數(shù)輸入基于頻域高階Rankine 源法開發(fā)的數(shù)值程序，對S175 以航速Fr=0.275、遭遇浪向β=150°在規(guī)則波中航行時的運動響應進行預報。數(shù)值預報采用如圖5 所示的經(jīng)過收斂性分析的計算域和網(wǎng)格：自由面邊界距船體1.5倍遭遇波長（λe=2πg/ω2e），單位遭遇波長內的網(wǎng)格節(jié)點數(shù)為20（自由面網(wǎng)格在數(shù)值程序中自動生成），船體網(wǎng)格的節(jié)點數(shù)為682。

圖5 S175船體和自由面網(wǎng)格劃分示意圖Fig.5 Mesh of S175 and free surface

圖6給出了采用上述網(wǎng)格計算得到的垂蕩、橫搖和縱搖運動響應隨入射波長的變化曲線（紅色實線）。為了進行比較，圖中還包含了試驗值[24]以及本文作者建立的移動脈動源[4]和半解析高階移動脈動源[25]的數(shù)值計算結果。移動脈動源法和半解析高階移動脈動源法在求解勢流問題時都是以移動脈動源格林函數(shù)為邊界積分方程的核函數(shù)，其中前者在以常值元離散邊界的基礎上，采用傳統(tǒng)的高斯求積公式計算格林函數(shù)的面積分；而后者選擇了與本文相同的高階單元離散邊界，且為了避免格林函數(shù)沿水平方向的高頻振蕩所引起的數(shù)值計算的不穩(wěn)定，采用了一種半解析的公式計算格林函數(shù)的面積分。從圖中的結果可以看出，本文的頻域高階Rankine源法和另外兩種數(shù)值方法計算得到的運動響應結果與試驗值的變化趨勢基本是一致的。但由于本文方法采用的核函數(shù)1/r計算簡單且不存在振蕩性，相比移動脈動源法，計算更為穩(wěn)定（尤其是縱搖運動）；由于在自由面和物面條件中考慮了船舶定常繞流的影響，相比基于均流線性化自由面條件的半解析高階移動脈動源法，預報的垂蕩和縱搖運動響應的結果在峰值附近與試驗值更為接近。

圖6 S175以航速Fr=0.275、遭遇浪向β=150°在規(guī)則波中航行時的運動響應Fig.6 Motion responses of S175 at Fr=0.275 with heading angle β=150°

4 不同浪向下船舶六自由度搖蕩運動數(shù)值模擬與分析

上述對S175船在遭遇浪向β=150°時的垂蕩、橫搖和縱搖運動的預報證實了本文建立的頻域高階Rankine在斜浪運動響應計算中的可靠性。在此基礎上，下面進一步對不同浪向下船舶的六自由度搖蕩運動進行模擬與分析。計算船型除了細長的S175 集裝箱船，還包括一艘在KVLCC2 基礎上改進得到的方形系數(shù)CB=0.84的散貨船S-Cb84[26]。

4.1 S175船

圖7給出了S175以Fr=0.275在不同浪向的規(guī)則波中航行時運動響應的計算結果與已有試驗值[24]的比較。除了兩個首斜浪（β=120°和150°）下六自由度搖蕩運動RAO，圖中還包括了橫浪時S175的橫蕩、垂蕩和橫搖運動結果。由圖可以看出，數(shù)值計算結果與已有試驗值吻合良好，表明本文方法適用于不同浪向下船舶運動響應的預報。比較三個不同浪向下S175的運動RAO曲線可以發(fā)現(xiàn)，它們隨入射波長的變化趨勢類似，但是曲線峰值（即運動共振）對應的波長會隨著浪向角的減小逐漸減?。ù故幒涂v搖最為明顯）。這主要是因為當船舶以一定的航速在波浪中航行時，波浪是以遭遇頻率的振蕩規(guī)律作用于船舶，運動共振發(fā)生在遭遇頻率等于固有頻率的工況。根據(jù)遭遇頻率的表達式ωe=ω0-kUcosβ可知，它是隨航速U、浪向角β和入射波頻率ω0（或波長）而變化的，當船舶以一定的航速在首斜浪中航行時，若浪向角β較小，則入射波頻率ω0必須較大時（即波長較?。┎拍苁沟迷庥鲱l率ωe等于固有頻率。

圖7 S175以航速Fr=0.275在不同浪向的規(guī)則波中航行時的運動響應Fig.7 Motion responses of S175 at Fr=0.275 with different heading angles

4.2 S-Cb84散貨船

采用本文方法進一步對較為肥大的散貨船S-Cb84 的運動響應進行計算。S-Cb84 的主尺度為：船長L=320 m，船寬B=58 m，吃水D=20.8 m，排水量Δ=324 000 t，縱搖慣性半徑kyy=0.25L，橫搖慣性半徑kxx=0.35B。圖8 為運動響應計算中采用的船體網(wǎng)格（節(jié)點數(shù)為962），自由面網(wǎng)格參數(shù)的設置策略與S175相同。

圖8 S-Cb84運動響應計算采用的船體網(wǎng)格Fig.8 Meshes of S-Cb84 used in the calculation of motion responses

圖9給出了S-Cb84以航速Fr=0.099在首斜浪β=150°和尾斜浪β=30°中航行時，六自由度運動響應隨入射波長的變化。對比數(shù)值結果和試驗值可以看出，無論是首斜浪還是尾斜浪，本文方法均可以較好地預報出S-Cb84的運動響應。

圖9 S-Cb84以航速Fr=0.099在首斜浪β=150°和尾斜浪β=30°中航行時的運動響應Fig.9 Motion responses of S-Cb84 at Fr=0.275 with heading angles β=150°and 30°

圖10進一步給出了橫浪時，在數(shù)值計算中分別考慮和不考慮橫搖粘性阻尼的橫搖和垂蕩運動響應結果的對比?？梢园l(fā)現(xiàn)，對于垂蕩運動，是否考慮橫搖粘性阻尼，結果幾乎沒有差別；但是對于橫搖運動，若忽略粘性阻尼的影響，本文基于頻域高階Rankine 源法的數(shù)值計算程序將會嚴重高估共振頻率附近的運動響應。

圖10 橫浪下（β=90°）考慮和不考慮橫搖粘性阻尼的S-Cb84垂蕩和橫搖運動響應計算結果對比Fig.10 Comparison between the results of heave and roll motion responses calculated by considering and neglecting the influence of viscous roll damping of S-Cb84 in beam waves(β=90°)

5 結 論

本文針對船舶在海上斜浪航行時的水動力學問題，基于頻域高階Rankine源法并結合CFD技術獲得的橫搖粘性阻尼系數(shù)，建立了預報船舶六自由度搖蕩運動響應的數(shù)值方法和程序。通過對不同船型在不同浪向和波長規(guī)則波中運動RAO的數(shù)值計算和分析，得出以下結論：

（1）對于橫搖運動，流體粘性影響顯著，忽略橫搖粘性阻尼會嚴重高估共振頻率附近的運動響應。本文采用CFD法模擬船舶自由橫搖衰減運動，經(jīng)能量法處理獲得橫搖阻尼系數(shù)，將其與頻域高階Rankine源法相結合可以更為準確地預報橫搖運動RAO；

（2）相比基于移動脈動源格林函數(shù)的數(shù)值方法，由于本文方法在求解繞輻射勢時采用了簡單易于計算的Rankine源1/r，并且在自由面條件中計入了船舶定常擾流的耦合影響，計算更為穩(wěn)定且運動響應的預報結果在峰值附近與試驗值更為接近；

（3）無論是細長的集裝箱船S175，還是肥大的散貨船S-Cb84，數(shù)值模擬得到的不同浪向下六自由度運動響應結果均與試驗值吻合良好，表明本文建立的數(shù)值方法在斜浪航行船舶搖蕩運動預報中，對于不同船型、不同浪向均有廣泛的適用性。