曹陽， 朱仁傳， 蔣銀， 洪亮

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,海洋工程國家重點實驗室，高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

船舶在波浪上的運動與阻力增加一直是船舶與海洋工程領(lǐng)域研究的熱點，準確預(yù)報船舶在波浪上的運動是船舶在水上航行、作業(yè)安全性和舒適性的基礎(chǔ)，船舶波浪增阻在船型優(yōu)化乃至航線優(yōu)化問題中是經(jīng)濟性的重要參數(shù)指標。

目前，關(guān)于波浪中航行船舶運動與增阻研究主要有船模試驗(experimental fluid mechanics, EFD)、基于勢流理論的理論計算[1]和考慮黏性的計算流體力學(xué)(CFD)數(shù)值模擬等方法。基于勢流理論的方法主要有切片法、三維面元法等，相對簡單且計算效率高，但勢流理論沒有考慮粘性的影響，且對諸如船體大幅度運動、波浪破碎等強非線性因素難以適用。而船模試驗往往需要進行費時費力的系列模型試驗，成本較高，試驗頻率也有限制，由于儀器設(shè)備和測試技術(shù)等方面的原因，目前還難以準確地監(jiān)測出船體周圍復(fù)雜流場的信息。

CFD方法由于考慮了粘性能對非線性問題作較為準確的數(shù)值模擬，隨著計算機計算能力的提高日漸成為求解船舶與海洋工程水動力問題的一種重要手段，近年來已經(jīng)取得了一系列重要的進展。Carrica等[2-3]利用重疊網(wǎng)格采用單相 level-set方法捕捉自由面，考慮了縱搖和垂蕩兩個自由度，對船模DTMB 5512在波浪中的自由運動的進行了數(shù)值模擬。GUO等[4-5]研究了KVLCC2在短波中的增阻，驗證了Faltinsen漸進公式預(yù)報船舶在短波中增阻的準確性，同時還對KVLCC2的耐波性做了研究，并分析了波浪對螺旋槳盤面處流場的影響。Sadat-Hosseini等[6]考慮縱蕩的影響，數(shù)值預(yù)報超大型油輪KVLCC2在不同波長迎浪工況下的波浪增阻和運動響應(yīng)，結(jié)果表明縱蕩對縱搖和波浪增阻幾乎沒有影響。國內(nèi)基于CFD方法關(guān)于船舶在波浪上運動與增阻的研究也取得了很多成果。方昭昭等[7-9]基于CFD方法建立了數(shù)值波浪水池，使用Fluent軟件結(jié)合動網(wǎng)格技術(shù)對不同航速頂浪航行的Wigley-III型船模的水動力系數(shù)、垂蕩、縱搖、波浪增阻等進行了數(shù)值計算，結(jié)果與實驗結(jié)果[10]吻合較好。沈志榮等[11-12]利用基于開源代碼OpenFOAM平臺開發(fā)計算分析了Wigley-III 型、DTMB5415、S175等船型在迎浪中的運動響應(yīng)以及波浪增阻。

CFD模擬船體運動需要使用動網(wǎng)格技術(shù)，涉及到網(wǎng)格的再生與變形，若船體運動響應(yīng)較大，網(wǎng)格畸變?nèi)菀滓l(fā)數(shù)值計算發(fā)散或是計算結(jié)果誤差大。重疊網(wǎng)格方法對于復(fù)雜曲面離散以及物體大幅運動具有無可比擬的優(yōu)勢，能方便處理復(fù)雜三維船體曲面，容易保證局部網(wǎng)格質(zhì)量，動態(tài)重疊網(wǎng)格在處理結(jié)構(gòu)物具有大幅運動的繞流問題時不需要網(wǎng)格再生，提高了動態(tài)網(wǎng)格處理效率。

本文基于粘性流理論使用CFD軟件建立了數(shù)值波浪水池并結(jié)合重疊網(wǎng)格方法，對航速為1.047 m/s(Fn=0.142 )縮尺比為58的KVLCC2模型在波浪上的運動進行了數(shù)值模擬，同時研究了KVLCC2在波浪中的阻力增加，分析了波阻成分。

1 數(shù)值波浪水池與重疊網(wǎng)格技術(shù)

本文的數(shù)值模擬是在一個三維數(shù)值波浪水池中進行，水池尺度隨入射波長的不同而不同，在出口處和兩側(cè)邊界附近設(shè)有人工阻尼消波區(qū)，消波區(qū)段約為2倍波長，船側(cè)消波區(qū)距船舷側(cè)約1倍船長，水池前端距船艏約1倍波長，尾端消波區(qū)距船艉約1倍波長。數(shù)值波浪水池及坐標系如圖1所示，固定坐標系O0x0y0z0原點位于船體中縱剖面、中橫剖面以及靜水面交點，x0軸指向船艏，y0軸指向船舷左側(cè)，z0軸垂直向上；參考坐標系Oxyz原點位于重心處，各坐標軸方向與固定坐標系對應(yīng)一致。

圖1 數(shù)值波浪水池及坐標系Fig.1 Numerical wave tank and the coordinate system

1.1 流體控制方程

整個流場以連續(xù)性方程和N-S方程為控制方程：

(1)

(2)

湍流模式為SSTk-ω兩方程模型，其控制方程為

(i=1,2,3)

(3)

式中：Γk、Γω為湍動能k和比耗散率ω的有效擴散系數(shù)，Yk、Yω為k和ω湍流耗散，Gk為平均速度梯度引起k的產(chǎn)生項，Gω為ω的產(chǎn)生項，Sω為交叉擴散項。

流體體積輸運方程為

(4)

本文在后面的數(shù)值計算中采用有限體積法離散上述控制方程，對流項使用二階迎風插值格式，擴散項的離散采用中心差分格式，應(yīng)用SIMPLE算法分離式求解；考慮重力影響，使用多重網(wǎng)格方法迭代求解離散代數(shù)方程組。

1.2 數(shù)值造波與消波

船舶在波浪中航行時的運動模擬，需要提供滿足計算要求的波浪環(huán)境，即首先要構(gòu)建數(shù)值波浪水池，本文采用速度邊界法造波。假定波浪沿x軸負方向傳播，根據(jù)波浪理論，參考坐標系下入射勢、遭遇頻率以及波面表達式為

(5)

參考坐標系下流體質(zhì)點的速度表達式為

(6)

在數(shù)值波浪水池的尾部設(shè)置有人工阻尼消波區(qū)，其長度約為入射波長的2倍。在消波區(qū)，對流體質(zhì)點垂向速度做強迫衰減，為保證流動的連續(xù)性，水平速度不做衰減處理[13]。強迫衰減的公式為

式中：wr(x,y,z;t)為衰減后的垂向速度；μ(x,z)為衰減函數(shù)；xs≤x≤xe，zb≤z≤zf；下標s和e分別代表阻尼區(qū)沿x方向的起點和終點；b和f分別代表沿z方向的底部和自由面；α為阻尼控制參數(shù)。

當流場中有船體存在時，還要在數(shù)值波浪水池的側(cè)面設(shè)置人工阻尼消波區(qū)，同樣只對流體質(zhì)點垂向速度做強迫衰減，強迫衰減公式為

式中：η(y,z)為衰減函數(shù)；yw≤y≤yn，zb≤z≤zf；下標w和n分別代表阻尼區(qū)沿y方向的起點和終點；b和f分別代表沿z方向的底部和自由面；α為阻尼控制參數(shù)。

1.3 重疊網(wǎng)格技術(shù)

為更好地模擬船體在波浪上的運動，本文采用重疊網(wǎng)格技術(shù)。重疊網(wǎng)格，也叫嵌合體網(wǎng)格，是一種區(qū)域分割與網(wǎng)格組合的策略，涉及到兩種區(qū)域：背景區(qū)域和重疊區(qū)域，兩種區(qū)域以任意方式疊加。如圖2所示，整個流場稱之為背景區(qū)域，該區(qū)域生成的網(wǎng)格為背景網(wǎng)格；包裹著船體隨船體做剛性運動的小區(qū)域為重疊區(qū)域，這個區(qū)域產(chǎn)生的網(wǎng)格稱為重疊網(wǎng)格，兩套網(wǎng)格之間通過交界面來識別。

圖2 重疊網(wǎng)格圖解Fig.2 The sketch of multi-region

重疊網(wǎng)格技術(shù)將網(wǎng)格分為有效單元和無效單元。只有有效單元參與離散控制方程的求解，無效單元指的是重疊區(qū)域內(nèi)被挖去的背景網(wǎng)格，不參與控制方程的求解，當重疊區(qū)域隨船體做剛性運動時背景網(wǎng)格中的部分有效單元和無效單元會相互轉(zhuǎn)化。重疊網(wǎng)格與背景網(wǎng)格的信息交換通過受者單元和貢獻單元之間插值來完成的，一般采用拉格朗日插值的思想，進行線性插值，具有以下形式[14]：

φrec=∑ξiφi

(11)

式中：φrec為受者單元的流動物理量；φi為其相鄰貢獻單元的流動物理量；對于二維問題，ξi為以相鄰貢獻單元中心為頂點組成三角形對應(yīng)的形函數(shù)；如果是三維問題，ξi是四面體所對應(yīng)的形函數(shù)，該四面體的頂點由相鄰貢獻單元的中心組成。

1.4 船體運動方程

本文數(shù)值模擬了KVLCC2在長波中迎浪航行工況下的運動，縱蕩運動較小且對其他運動以及波浪增阻幾乎沒有影響，此處只考慮垂蕩和縱搖兩個自由度，其運動方程為

(12)

為避免作用于船體的力或力矩過大造成船體劇烈的運動，在初始時給船體受到的力和力矩乘上一個阻尼函數(shù)fr，定義如下

其中，是時刻t時相關(guān)矩陣C的估計。我們可以使用指數(shù)加權(quán)或者滑動加窗的來估計。其中，關(guān)于的最簡單的選擇就是在最小均方算法中的。因此，信號子空間的表達式可更新為：

(13)

式中：ts為船體可以開始運動的時間，tr為船體開始運動后作用于船體的力和力矩受到人工抑制的時間。

2 波浪中航行船舶的數(shù)值模擬

本文數(shù)值模擬的是迎浪工況，流場關(guān)于船體中縱剖面對稱，計算域取流場的一半即可。數(shù)值模擬計算主要考慮船體垂蕩和縱搖兩個自由度，部分短波工況，也采用固定模進行了數(shù)值模擬。

2.1 船體幾何參數(shù)

KVLCC2是II型MOERI(maritime and ocean engineering research institute)油輪，相較于廣泛研究的Wigley、S60以及S175船模具有更復(fù)雜的船體外形，其主要船型參數(shù)如表1所示，各站橫剖面如圖3所示。

表1 KVLCC2及模型主尺度

2.2 邊界條件與網(wǎng)格劃分

波浪中航行船舶水動力計算域的邊界條件設(shè)定如下：前端及上下面均為速度入口；左右為對稱邊界條件；長波時尾端采用壓力出口，短波時尾端可以視為遠場，用速度入口作為邊界條件；船體設(shè)定為壁面邊界條件。

圖3 KVLCC2各站橫剖面圖Fig.3 The cross section of KVLCC2

計算區(qū)域網(wǎng)格的劃分主要從如下幾個方面考慮：沿波浪傳播方向保證足夠數(shù)量的網(wǎng)格，以避免數(shù)值耗散引起的波浪幅值的沿程衰減；垂向沿自由面附近要保證足夠數(shù)量的網(wǎng)格，以精準捕捉自由液面；為準確捕捉尾跡，在開爾文興波區(qū)范圍內(nèi)也進行了網(wǎng)格加密；在流場變化比較劇烈的地方如球鼻艏以及船艉附近要進行網(wǎng)格加密；為了獲得準確的船體周圍流場信息，船體附近網(wǎng)格劃分較遠場網(wǎng)格劃分應(yīng)精細一些；在消波區(qū)，保證網(wǎng)格平緩過渡的情況下適當拉大x和y方向網(wǎng)格的尺度，既減少計算量又可以起到數(shù)值消波的作用。船體壁面設(shè)置5層邊界層網(wǎng)格，y+取50左右，為減少網(wǎng)格量，甲板不設(shè)邊界層。

主要工況計算所用網(wǎng)格參數(shù)及時間步如表2所示，計算所用網(wǎng)格劃分如圖4所示。

2.3 數(shù)值模擬結(jié)果的處理

本文采用傅里葉級數(shù)展開法，使用傅里葉級數(shù)對數(shù)值模擬得到的垂蕩、縱搖以及阻力時歷曲線進行擬合。

表2 主要工況及計算所用網(wǎng)格參數(shù)(Fn=0.142)Table 2 The main conditions simulated and the parameters used for mesh generation(Fn=0.142)

圖4 計算域網(wǎng)格劃分Fig.4 Mesh used for numerical simulation

Rwave=Rt+R1stsin(wet+γ1)+…

(14)

式中：Rwave為波浪力；ωe為遭遇頻率；Rt為擬合得到的傅里葉級數(shù)的定常值，即平均波浪力；R1st為一階波浪力幅值；γ1為一階波浪力的相位。

圖5 傅里葉級數(shù)展開法示例(Fn=0.142、L=4.965 m、ξa=0.075 m)Fig.5 Example of Fourier series expansion method (Fn=0.142，L=4.965 m，ζa=0.075 m)

3 KVLCC2在波浪中的運動與增阻

勢流理論使用的是由上海交通大學(xué)水動力學(xué)實驗室基于二維勢流理論開發(fā)的切片程序計算的結(jié)果，實驗結(jié)果取自大阪大學(xué)(OU)、意大利船舶研究所(INSEAN)以及挪威科技大學(xué)(NTNU)三家單位發(fā)表在相關(guān)文獻上的實驗數(shù)據(jù)[15]。

3.1 垂蕩和縱搖

采用傅里葉級數(shù)擬合運動時歷曲線分析可得不同頻率遭遇波下船體運動響應(yīng)的幅值和相位。KVLCC2垂蕩和縱搖幅值分別按式(15)進行無因次化，其幅頻響應(yīng)曲線如圖6、7所示。

(15)

式中：ξa為遭遇波波幅，K為波數(shù)，X3和X5分別為垂蕩和縱搖的幅值。

圖6 Fn=0.142時KVLCC2的垂蕩幅值Fig.6 RAO heave of KVLCC2 in head waves, Fn=0.142

圖7 Fn=0.142時KVLCC2的縱搖幅值Fig.7 RAO pitch of KVLCC2 in head waves, Fn=0.142

從圖6、7中可以看出，數(shù)值模擬的結(jié)果與實驗結(jié)果整體吻合較好。在0.7

圖8、9給出了KVLCC2垂蕩和縱搖運動與遭遇波的相對相位隨Lpp/λ變化情況，數(shù)值模擬的結(jié)果同實驗結(jié)果整體符合較好。1.01兩個區(qū)間內(nèi)分別隨Lpp/λ增大而單調(diào)遞減；0.5

圖8 Fn=0.142時KVLCC2的垂蕩相位Fig.8 Heave phase of KVLCC2 in head waves, Fn=0.142

圖9 Fn=0.142時KVLCC2的縱搖相位Fig.9 Pitch phase of KVLCC2 in head waves, Fn=0.142

3.2 KVLCC2在波浪中的增阻

(16)

圖10為Fn=0.142 時KVLCC2波浪增阻的CFD模擬結(jié)果與切片法計算結(jié)果以及實驗結(jié)果的對比，其中切片法采用輻射能量法計算波浪增阻，可以看到，CFD模擬的結(jié)果同實驗結(jié)果吻合較好。切片法計算的增阻值整體上小于實驗值，增阻幅值同實驗值差異較大，主要原因是：勢流方法沒有考慮粘性的影響，流體的粘性也貢獻了一部分增阻；本文勢流切片方法采用輻射能量法計算波浪增阻，沒有考慮繞射對增阻的貢獻。

圖10 Fn=0.142時KVLCC2在波浪上的增阻Fig.10 Added resistance coefficient of KVLCC2 in head waves, Fn=0.142

從作用力方向分析，總的波浪增阻由壓力增阻、摩擦增阻兩部分組成；從物理問題角度出發(fā)，波浪增阻又可分解為輻射增阻和繞射增阻。圖11給出了遭遇頻率對增阻成分的影響。由圖11(a)可以看到，在全波長范圍內(nèi)波浪增阻基本上由壓力增阻貢獻，由摩擦力引起的增阻很小。從圖11(b)可以看出，波長對繞射增阻影響較小，輻射增阻在Lpp/λ>1.3時趨于0，在長波范圍，輻射增阻隨著波長的增大而減小。

3.3 KVLCC2在短波中的增阻

由上述分析可知，船舶在短波中的增阻主要由繞射增阻貢獻，Lpp/λ>1.3時輻射增阻可以忽略不計，現(xiàn)在對KVLCC2在短波中的增阻成分作進一步驗證。

圖12中，CFD方法計算的是Fn=0.142 時KVLCC2固定模在短波中的增阻值，實驗測量的是Fn=0.142 時KVLCC2自由模在短波中的增阻值。其中“Faltinsen”是用Faltinsen的漸近公式[14]計算得到的波浪增阻，頂浪工況下漸進公式為

(17)

式中:C_l表示非遮蔽水線，ω0是波浪的頻率，θ是水線的切線與x軸的夾角，ζa表示波幅。

可以看到，KVLCC2在短波范圍內(nèi)的增阻值隨波長基本變化不大，不計入輻射增阻的CFD方法計算值以及Faltinsen漸進公式的計算值都同實驗值吻合的非常好，也說明了本文數(shù)值模擬方法預(yù)報船舶在短波中的波浪增阻相當有效。

圖11 Fn=0.142時KVLCC2在波浪上增阻的構(gòu)成Fig.11 Added resistance components of KVLCC2 in head waves, Fn=0.142

圖12 Fn=0.142時KVLCC2在短波中的增阻Fig.12 Added resistance coefficient of KVLCC2 in head short waves, Fn=0.142

4 結(jié)論

1)結(jié)合重疊網(wǎng)格技術(shù)的CFD方法能夠準確處理船舶在波浪上的運動響應(yīng)，相較于勢流理論， CFD方法能夠在全波長范圍內(nèi)更準確地預(yù)報船舶在波浪上增阻；

2)船體遭遇短波時，船體輻射運動非常小，可以忽略，隨著波長增加，垂蕩幅值和縱搖幅值都迅速增大，船體在長波中航行時無因次的運動幅值趨于1；遭遇足夠長的入射波，船體垂蕩與波浪同步，相位差0°，縱搖運動相位則落后波浪四分之一周期，船體處于“隨浪逐流”狀態(tài)；

3)Lpp/λ=0.9左右時波浪增阻達到峰值；在短波中船體在波浪中的運動響應(yīng)比較小，波浪增阻主要由繞射增阻貢獻，輻射增阻可以忽略不計；長波范圍，船體輻射運動引起的增阻隨波浪增大逐漸減?。淮霸诓ɡ酥械脑鲎柚饕蓧毫龅母淖円?，摩擦對增阻的貢獻非常小。

[1] 繆國平, 劉應(yīng)中. 船舶在波浪上的運動理論[M]. 上海: 上海交通大學(xué)出版社, 1989.

[2] CARRICA P M, WILSON R V, NOACK R W, et al. Ship motions using single-phase level set with dynamic overset grids[J]. Computers & fluids, 2007, 36(9): 1415-1433.

[3] CARRICA P M, WILSON R V, STERN F. Unsteady RANS simulation of the ship forward speed diffraction problem [J]. Computers & fluids, 2006, 35(6): 545-570.

[4] GUO B, STEEN S. Evaluation of added resistance of kvlcc2 in short waves [J]. Journal of hydrodynamics, Ser B, 2011, 23(6): 709-722.

[5] GUO B J, STEEN S, DENG G B. Seakeeping prediction of KVLCC2 in head waves with RANS [J]. Applied ocean research, 2012, 35: 56-67.

[6] SADAT-HOSSEINI H, WU P, CARRICA P M, et al. CFD verification and validation of added resistance and motions of KVLCC2 with fixed and free surge in short and long head waves[J]. Ocean engineering, 2013, 59: 240-273.

[7] 方昭昭, 趙丙乾, 朱仁傳. 頂浪中船舶運動的數(shù)值模擬與波浪增阻計算[J]. 中國造船, 2014, 55(2): 8-17.

FANG Zhaozhao, ZHAO Bingqian, ZHU Renchuan. Numerical simulation of ship motion and calculation of added resistance in heading waves [J]. Shipbuilding of China, 2014, 55(2): 8-17.

[8] 方昭昭, 朱仁傳, 繆國平, 等. 基于數(shù)值波浪水池的波浪中船舶水動力計算[J]. 水動力學(xué)研究與進展A輯, 2012, 27(5): 515-524.

FANG Zhaozhao, ZHU Renchuan, MOU Guoping, et al. Numerical calculation of hydrodynamic forces for a ship in regular waves based on numerical wave tank[J]. Journal of hydrodynamics, 2012, 27(5): 515-524.

[9] 方昭昭, 朱仁傳, 繆國平, 等. 數(shù)值波浪水池中航行船舶繞射問題的數(shù)值模擬[J]. 上海交通大學(xué)學(xué)報, 2012, 46(8): 1203-1209.

FANG Zhaozhao, ZHU Renchuan, MOU Guoping, et al. Numerical simulation of diffraction problems of moving vessels in numerical wave tank[J]. Journey of Shanghai JiaoTong University, 2012, 46(8): 1203-1209.

[10] OURNéE J M J. Experiments and calculations on 4 Wigley hull forms in head waves[J]. Delft University of Technology, 1992： 909.

[11] YE H, SHEN Z, WAN D. Numerical prediction of added resistance and vertical ship motions in regular head waves [J]. Journal of marine science and application, 2012, 11(4): 410-416.

[12] 沈志榮, 葉海軒, 萬德成. 船舶在迎浪中運動響應(yīng)和波浪增阻的RANS數(shù)值模擬[J]. 水動力學(xué)研究與進展A輯, 2012, 27(6): 621-633.

SHEN Zhirong, YE Haixuan, WAN Decheng. Motion response and added resistance of ship in head waves based on RANS simulations[J]. Journal of hydrodynamics, 2012, 27(6): 621-633.

[13] 齊鵬, 王永學(xué). 三維數(shù)值波浪水池技術(shù)與應(yīng)用[J]. 大連理工大學(xué)學(xué)報, 2003,43(6): 825-830.

QI Peng, WANG Yongxue. 3-D numerical-wave-tank technology and its application [J]. Journey of Dalian University of Technology, 2003, 43(6): 825-830.

[14] FALTINSEN O M, MINSAAS K J, LIAPIS N, et al. Prediction of resistance and propulsion of a ship in a seaway[C]∥Proceeding of 13th Symposium on Naval Hydrodynamics, Tokyo, Japan, 1980.

[15] LARSSON L, STERN F, VISONNEAU M. Numerical ship hydrodynamics [C]//A Workshop on Numerical Ship Hydrodynamics. Gothenburg, Sweden, 2010.

本文引用格式：

曹陽， 朱仁傳， 蔣銀， 等. 波浪中KVLCC2運動與阻力增加的CFD計算及分析[J]. 哈爾濱工程大學(xué)學(xué)報， 2017， 38(12): 1828-1835.

CAO Yang, ZHU Renchuan, JIANG Yin, et al. CFD verification and analysis of added resistance and motions of KVLCC2 in head waves[J]. Journal of Harbin Engineering University， 2017， 38(12): 1828-1835.