周新聰，肖仲歧，張 聰，歐陽(yáng)武，黃 健

(1. 武漢理工大學(xué) 能源與動(dòng)力工程學(xué)院，湖北 武漢 430063；2. 船舶動(dòng)力工程技術(shù)交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430063；3. 國(guó)家水運(yùn)安全工程技術(shù)研究中心可靠性與新能源研究所，湖北 武漢 430063)

0 引 言

大型豪華郵輪作為一種擁有高附加值、高技術(shù)要求以及高難度的船舶，一直被譽(yù)為“海上流動(dòng)五星級(jí)賓館”，對(duì)于中國(guó)船舶行業(yè)來(lái)說(shuō)這也是一塊未經(jīng)開(kāi)發(fā)的空白領(lǐng)域。大型豪華郵輪的建造意味著我國(guó)在相關(guān)領(lǐng)域上的重大突破，而隨著國(guó)際海事組織對(duì)于能效日益嚴(yán)格的要求以及郵輪本身對(duì)于經(jīng)濟(jì)性以及舒適性的追求，有必要追求最佳的能效表現(xiàn)。因此，有必要研究船舶在真實(shí)海況下航行的運(yùn)動(dòng)性能。

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

1 數(shù)值模擬

1.1 控制方程

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

質(zhì)量守恒

動(dòng)量方程

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

自由液面追蹤采用VOF法，該方法假設(shè)網(wǎng)格大小足以求解各相之間的交界面位置及形狀。通過(guò)對(duì)每個(gè)網(wǎng)格單元中各相的體積分?jǐn)?shù)αi來(lái)描述交界面相的分布與位置，當(dāng) αi=0 時(shí)代表該網(wǎng)格單元中不存在相i；αi=1 時(shí)代表該網(wǎng)格單元完全由相i填充；而0＜αi＜1時(shí)代表該網(wǎng)格處于相的交界面上。

1.2 邊界造波與消波

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

其頻譜表示為[4]：

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

為了避免VOF波在邊界發(fā)生反射對(duì)計(jì)算結(jié)果產(chǎn)生干擾，于出口及計(jì)算域側(cè)面設(shè)置具有波阻尼的消波區(qū)，其阻力方程為[5]：

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

1.3 六自由度剛體運(yùn)動(dòng)

在動(dòng)態(tài)流體相互作用(DFBI)分析中，流體區(qū)域與六自由度剛體耦合，根據(jù)作用力計(jì)算六自由度剛體的運(yùn)動(dòng)。然后，流體網(wǎng)格嚴(yán)格按照計(jì)算得出的剛體運(yùn)動(dòng)進(jìn)行移動(dòng)。在求解中具有2個(gè)坐標(biāo)系，即基于全局的慣性坐標(biāo)系以及以船體質(zhì)心為原點(diǎn)的局部坐標(biāo)系。

由于求解對(duì)象為迎浪航行的船舶，其運(yùn)動(dòng)包括沿z軸方向的平移及以y軸為軸線的旋轉(zhuǎn)，即升沉與縱傾，故其運(yùn)動(dòng)方程包括平移及旋轉(zhuǎn)。其平移方程依據(jù)全局慣性坐標(biāo)系建立：

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

其體旋轉(zhuǎn)方程依托于以船體質(zhì)心為原點(diǎn)的局部坐標(biāo)系：

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

2 仿真模型

本文以某大型郵輪為計(jì)算對(duì)象，其相關(guān)參數(shù)見(jiàn)表1。

表1 目標(biāo)郵輪船體相關(guān)參數(shù)Tab. 1 Relevant parameters of target cruise ship

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

圖1 船體模型Fig. 1 Hull model

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

圖2 計(jì)算域劃分Fig. 2 Division of computing domain

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

圖3 局部加密網(wǎng)格劃分Fig. 3 Local mesh refinement

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

在此基礎(chǔ)上，對(duì)網(wǎng)格尺寸對(duì)精度的影響進(jìn)行研究，劃分了3套網(wǎng)格，其基礎(chǔ)尺寸分別為 Sbasic， 1.225Sbasic以及1 .5Sbasic。以工況6為目標(biāo)，計(jì)算結(jié)果與《港口工程荷載規(guī)范》提供的公式計(jì)算[8]所得理論結(jié)果對(duì)比，如表2所示。

表2 網(wǎng)格收斂性分析Tab. 2 Analysis of grid convergence

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

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

3 計(jì)算結(jié)果

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

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

圖5 運(yùn)動(dòng)響應(yīng)時(shí)歷變化曲線Fig. 5 Time series curve of motion response

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

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

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

表4 仿真計(jì)算結(jié)果Tab. 4 simulation results

4 結(jié)果對(duì)比分析

4.1 海流流速影響

對(duì)比海流流速不同的3個(gè)工況，其阻力以及根據(jù)靜水工況修正后的運(yùn)動(dòng)響應(yīng)時(shí)歷曲線如圖8～圖10所示。

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

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

圖10 工況1～工況3升沉?xí)r歷對(duì)比Fig. 10 Comparison of heave time histories of working conditions 1, 2 and 3

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

圖11 運(yùn)動(dòng)響應(yīng)幅值分布估計(jì)Fig. 11 Estimation of amplitude distribution of motion response

由圖11可知，當(dāng)海流流速發(fā)生變化時(shí)，郵輪的縱傾運(yùn)動(dòng)只在平衡位置附近發(fā)生一定的變化，而其升沉運(yùn)動(dòng)隨著海流變慢，越來(lái)越傾向于減小吃水，同時(shí)也更易偏離平衡位置。

4.2 海浪周期影響

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

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

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

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

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

5 結(jié) 語(yǔ)

1）船舶迎浪航行時(shí)波浪起伏造成的船體運(yùn)動(dòng)響應(yīng)其縱傾與升沉變化的周期一致且正負(fù)相反，運(yùn)動(dòng)響應(yīng)與阻力變化之間有很大的相關(guān)性，但其影響集中于壓差阻力，屬于船體運(yùn)動(dòng)興波帶來(lái)的阻力變化，對(duì)于剪切阻力影響較小。

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

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