海浪流場作用下的郵輪實船運動仿真

周新聰，肖仲歧，張 聰，歐陽武，黃 健

(1. 武漢理工大學 能源與動力工程學院，湖北 武漢 430063；2. 船舶動力工程技術交通行業重點實驗室，湖北 武漢 430063；3. 國家水運安全工程技術研究中心可靠性與新能源研究所，湖北 武漢 430063)

0 引 言

大型豪華郵輪作為一種擁有高附加值、高技術要求以及高難度的船舶，一直被譽為“海上流動五星級賓館”，對于中國船舶行業來說這也是一塊未經開發的空白領域。大型豪華郵輪的建造意味著我國在相關領域上的重大突破，而隨著國際海事組織對于能效日益嚴格的要求以及郵輪本身對于經濟性以及舒適性的追求，有必要追求最佳的能效表現。因此，有必要研究船舶在真實海況下航行的運動性能。

Zhao等[1]建立了評估船舶在真實海況下航行能效的綜合化船舶推進模型。Yang等[2]對LNG船進行了考慮液艙晃蕩與船體耦合的迎浪和斜浪中的數值模擬。廖煒昊[3]利用簡化模型對波浪增阻與推進主機功率變化關系進行了研究。本文針對某大型郵輪基于粘性流的CFD理論建立了數值波浪水池，對郵輪在基于JONSWAP海浪譜的不規則波浪中航行進行流體仿真，模擬了一定的海況下船體所受阻力變化以及運動響應，通過對計算結果的對比為郵輪能效的提高提供一定依據。

1 數值模擬

1.1 控制方程

本文采用商用仿真軟件STAR-CCM+對數值波浪水池中流場及船舶運動進行仿真，其流體運動依托于RANS的k-ε湍流模型，船體在三維數值波浪水池中的運動響應屬于不可壓縮多相流流動同時不涉及熱量的傳遞，只需考慮質量以及動量的守恒，故其基本控制方程組如下：

質量守恒

動量方程

其中：ρ為密度，即單位體積的質量；v為連續體速度；?表示克羅內克乘積；fb為作用于連續體的單位體積的各種體積力（例如重力和離心力）的合力；p為壓力；T為粘性應力張量，I為單位張量，－?·(pI)+?·T=?σ ，σ=－pI+T意為應力張量為法向應力與剪切應力之和。

自由液面追蹤采用VOF法，該方法假設網格大小足以求解各相之間的交界面位置及形狀。通過對每個網格單元中各相的體積分數αi來描述交界面相的分布與位置，當 αi=0 時代表該網格單元中不存在相i；αi=1 時代表該網格單元完全由相i填充；而0＜αi＜1時代表該網格處于相的交界面上。

1.2 邊界造波與消波

為了盡可能模擬真實海況，本文選用基于JONSWAP海浪譜生成的不規則波，同時隨機數發生器將使用固定種子以提供相同的波形圖案從而方便對比。

其頻譜表示為[4]：

其中： γ=3.3， 屬于無量綱峰形參數；Aγ=1－0.287 ln(γ) ，為歸一化因子；σ為頻譜寬度參數，小于或等于譜峰頻率ωp時 ，取 σa=0.07，大于譜峰頻率時，取σb=0.09。

為了避免VOF波在邊界發生反射對計算結果產生干擾，于出口及計算域側面設置具有波阻尼的消波區，其阻力方程為[5]：

其中，xsd為波阻尼區域在x方向上的起點；xed為其在x方向上的終點，即邊界；f1，f2和nd為阻尼模型的控制參數；w為垂直方向上的速度分量。

1.3 六自由度剛體運動

在動態流體相互作用(DFBI)分析中，流體區域與六自由度剛體耦合，根據作用力計算六自由度剛體的運動。然后，流體網格嚴格按照計算得出的剛體運動進行移動。在求解中具有2個坐標系，即基于全局的慣性坐標系以及以船體質心為原點的局部坐標系。

由于求解對象為迎浪航行的船舶，其運動包括沿z軸方向的平移及以y軸為軸線的旋轉，即升沉與縱傾，故其運動方程包括平移及旋轉。其平移方程依據全局慣性坐標系建立：

其中：m為體質量；v為 平移速度；f為船體所受合力。

其體旋轉方程依托于以船體質心為原點的局部坐標系：

其中：M為慣性矩張量；ω為剛體角速度；n為船體所受合力矩。

2 仿真模型

本文以某大型郵輪為計算對象，其相關參數見表1。

表1 目標郵輪船體相關參數Tab. 1 Relevant parameters of target cruise ship

建立實尺度模型后，其具體外形如圖1所示。

圖1 船體模型Fig. 1 Hull model

根據模型尺度劃分計算域，由于實尺度模型較大，考慮船體的對稱性，計算域取船舶左側一半，不規則波的最大波長 λmax=60·H1/3，H1/3為有義波高，其寬度及水深約為2倍最大波長，水面以上部分高度約為1倍最大波長，水池前段既其波浪水流入口處位于船首部上游約2倍最大波長處，而尾端也即水池出口位于船尾部4倍最大波長處，見圖2。

圖2 計算域劃分Fig. 2 Division of computing domain

網格劃分時針對船體周圍水域及自由水面進行局部加密，根據文獻[6]中對于網格收斂性的論證，對以自由液面為中心，上下1.5倍H1/3范圍內的網格加密，其中X和Y方向上的網格尺寸不得超出 λ/n1，Z方向上尺寸約為H1/3/n2，時間步長小于Tp/2.4·n1，其中H1/3為有義波高，λ 為波長，n1=80～100,n2=20，Tp為譜峰周期。為了在精度與計算時間之間取得平衡，對于遠離自由液面及船體表面的區域過渡至較粗的網格，最終網格劃分結果如圖3所示。

圖3 局部加密網格劃分Fig. 3 Local mesh refinement

不同于船舶縮比模型尺度下的仿真，為確保模擬結果的準確性，船體壁面第一層網格厚度應使壁面Y+值大體保持在30～300之間，實船尺度下由于雷諾數的增加，壁面Y+值往往會大于300，此時船體Y+值除少數位置外應該保持不大于10 000[7]，同時為確保模擬精度，設定實船平均沙粒粗糙度hs為0.03 mm。

在此基礎上，對網格尺寸對精度的影響進行研究，劃分了3套網格，其基礎尺寸分別為 Sbasic， 1.225Sbasic以及1 .5Sbasic。以工況6為目標，計算結果與《港口工程荷載規范》提供的公式計算[8]所得理論結果對比，如表2所示。

表2 網格收斂性分析Tab. 2 Analysis of grid convergence

可以看出，隨著網格的加密，其計算精度逐漸上升，而當其基礎尺寸為 Sbasic時，誤差已足夠小，同時若再加密網格，網格數大于1 100萬，其計算時間將進一步拉長，故最終選定的網格尺寸其劃分的網格總數為1 150萬左右，在這個基礎上時間步長度向下取整，最終長度為0.04 s。此時，船體表面Y+值分布如圖4所示。

圖4 船體表面Y+值分布Fig. 4 Distribution of Y + value on hull surface

3 計算結果

數值模擬時通過動態流體相互作用（DFBI）分析其運動響應，模型釋放時間為1 s，緩沖時間為10 s，入口處持續生成波浪，同時計算時間步長取0.04 s，工況見表3。以郵輪在海流流速為1.3 m/s，有義波高為2.8 m且譜峰周期為10 s的海況下航行的仿真結果為例，計算結果如圖5和圖6所示。

表3 工況列表Tab. 3 List of working conditions

圖5 運動響應時歷變化曲線Fig. 5 Time series curve of motion response

圖6 剪切及壓差阻力時歷變化曲線Fig. 6 Time series curve of shear and differential pressure resistance

由圖5和圖6可知，當該船在海浪水面頂浪航行時，阻力隨運動響應發生變化，其中剪切阻力基本維持穩定，壓差阻力根據船舶運動響應狀態變化。最終通過模擬得到郵輪船身在不同海況下頂浪逆流時的流場示意如圖7所示，阻力平均值R、縱傾均值θ以及升沉均值h如表4所示。

圖7 船舶頂浪航行興波Fig. 7 Making waves by sailing on top of waves

表4 仿真計算結果Tab. 4 simulation results

4 結果對比分析

4.1 海流流速影響

對比海流流速不同的3個工況，其阻力以及根據靜水工況修正后的運動響應時歷曲線如圖8～圖10所示。

圖8 工況1～工況3阻力時歷對比Fig. 8 Comparison of resistance time histories of working conditions 1, 2 and 3

圖9 工況1～工況3縱傾時歷對比Fig. 9 Comparison of trim time histories of working conditions 1, 2 and 3

圖10 工況1～工況3升沉時歷對比Fig. 10 Comparison of heave time histories of working conditions 1, 2 and 3

通過分析可得，隨著海流流速的減小，其最開始阻力及運動響應的變化趨勢基本不變，然而在經過一段時間之后，流速較慢的工況在時歷變化上開始出現明顯的滯后，且滯后隨著時間的推移越來越大。具體原因可能為海流流速直接影響到了船體遇浪的相對速度，使其在遭遇海浪的頻率上出現了變化。對船體穩定之后的運動響應曲線進行分布估計，如圖11所示。

圖11 運動響應幅值分布估計Fig. 11 Estimation of amplitude distribution of motion response

由圖11可知，當海流流速發生變化時，郵輪的縱傾運動只在平衡位置附近發生一定的變化，而其升沉運動隨著海流變慢，越來越傾向于減小吃水，同時也更易偏離平衡位置。

4.2 海浪周期影響

對于海浪譜峰周期變化的工況1、工況4、工況5，其阻力及修正后的運動響應曲線對比以及分布估計如圖12所示。在其他條件不變的情況下，隨著譜峰周期的變短，其波長隨之變短，此時其阻力及運動響應總體趨勢不發生變化，但其上下波動的范圍隨著海浪周期的縮短而變小，趨近于短時間內的平衡位置。其分布估計也表明，隨著周期的縮短，其縱傾及橫搖越來越多地集中于平衡位置附近。

圖12 海浪譜峰周期變化下的對比Fig. 12 Comparison of changes caused by wave spectrum peak period

其剪切壓差阻力隨工況變化如圖13所示。

圖13 阻力構成占比隨工況的變化Fig. 13 Comparison of resistance composition ratio with working conditions

隨著海流流速以及譜峰周期的減小，其剪切阻力占比都會增加。海流流速越小，其剪切阻力占比增加的越多，而隨著譜峰周期的減小，其剪切壓力占比增加的量越少。

5 結 語

1）船舶迎浪航行時波浪起伏造成的船體運動響應其縱傾與升沉變化的周期一致且正負相反，運動響應與阻力變化之間有很大的相關性，但其影響集中于壓差阻力，屬于船體運動興波帶來的阻力變化，對于剪切阻力影響較小。

2）對于海浪流的迎浪航行在一定程度上會加深吃水且加劇其升沉，對于縱傾影響較小，而波浪周期足夠短時更有利于船舶航行的穩定。

3）基于真實海況的實船頂浪航行數值模擬為郵輪耐波性及吊艙的選型奠定了基礎，后續可開展基于真實海況的吊艙船體耦合分析。