基于重疊網(wǎng)格方法的船舶迎浪增阻與運(yùn)動數(shù)值預(yù)報(bào)

李冬琴，李國煥，戴晶晶，章易立，李 鵬

(1. 江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003；2. 江蘇現(xiàn)代造船技術(shù)有限公司，江蘇 鎮(zhèn)江 212003；3. 廣州船舶及海洋工程設(shè)計(jì)研究院，廣東 廣州 510250)

0 引 言

船舶在波浪中的快速性和耐波性是衡量船舶性能的重要指標(biāo)，一般而言，與靜水中的船舶阻力相比，波浪中阻力增加了約15%～30%[1]，導(dǎo)致船舶產(chǎn)生額外的燃料消耗和失速現(xiàn)象，而船舶在波浪中所產(chǎn)生的縱搖和垂蕩運(yùn)動對船舶安全性和舒適性等造成了不利影響，鑒于目前國際海事組織已經(jīng)推出了船舶能效設(shè)計(jì)指數(shù)（EEDI）和船舶能效營運(yùn)指數(shù)（EEOI），準(zhǔn)確預(yù)報(bào)船舶在波浪中航行的阻力增值和運(yùn)動響應(yīng)十分必要。

目前研究波浪增阻的方法主要有船模試驗(yàn)、勢流理論計(jì)算和全粘流CFD軟件數(shù)值模擬等。船模試驗(yàn)因其復(fù)雜性與成本較高等原因，往往僅作為對最終計(jì)算結(jié)果的驗(yàn)證。而勢流理論計(jì)算相對較為簡單高效，但沒有考慮粘性的影響。勢流理論主要包括切片法、三維面元法等，最早開始研究波浪增阻的是Havelock[2]，其通過研究零航速的圓柱體在波浪中的受力，提出一個(gè)簡單的公式計(jì)算船舶在波浪中的阻力增值。Maruo[3]提出了利用勢流方法進(jìn)行求解，并指出主要是由船舶在波浪中的垂蕩和縱搖運(yùn)動引起的阻力增加。而隨著計(jì)算機(jī)技術(shù)與計(jì)算流體動力學(xué)的高速發(fā)展，CFD方法因其充分考慮粘性作用，且對非線性船舶運(yùn)動和力的響應(yīng)能夠作出較為準(zhǔn)確的預(yù)報(bào)，已越來越多地應(yīng)用到船舶與海洋工程水動力學(xué)領(lǐng)域。國外的Orihara和Miyata[4]采用了 WISDAM-X 軟件和 overset mesh 方法分析了SR-108船舶的阻力增值，結(jié)果表明該船首能夠有效減小波浪中的阻力。Simonsen等[5]利用URANS代碼計(jì)算了KCS船模（縮尺比為52.667）在規(guī)則波中的增阻與運(yùn)動響應(yīng)，并進(jìn)行了相應(yīng)的模型試驗(yàn)研究。Seonguk Seo等[6]利用開源代碼Open FOAM平臺對KCS船模（縮尺比為37.9）在規(guī)則波中迎浪航行進(jìn)行模擬，同時(shí)評估了網(wǎng)格的不確定性影響，并和模型試驗(yàn)數(shù)據(jù)進(jìn)行比較。而國內(nèi)的吳乘勝等[7]建立了數(shù)值波浪水池，并對高速三體船波浪中航行進(jìn)行數(shù)值模擬。沈志榮等[8 - 10]利用Open FOAM工具箱計(jì)算分析了Wigley III型、DTMB5415、S175等船型的阻力與耐波性能。查若思[11]使用naoe-FOAM-SJTU求解器分別對單體船和雙體船在典型規(guī)則波海況下航行進(jìn)行數(shù)值模擬。曹陽等[12]利用重疊網(wǎng)格方法研究了KVLCC2船模在規(guī)則波中的增阻與運(yùn)動響應(yīng)，分析了波浪增阻成分，并同勢流理論計(jì)算結(jié)果及試驗(yàn)結(jié)果進(jìn)行對比。

本文基于CFD軟件對KSC船模在規(guī)則波中迎浪航行進(jìn)行數(shù)值模擬，首先建立了數(shù)值波浪水池，利用入口處設(shè)置波速函數(shù)的方法和強(qiáng)迫消波技術(shù)，成功模擬了5個(gè)不同波長規(guī)則波的生成與傳播；然后結(jié)合重疊網(wǎng)格方法，對航速為2.017 m/s、縮尺比為37.9的（帶舵的）KCS船模在靜水和波浪中的阻力與運(yùn)動響應(yīng)進(jìn)行數(shù)值計(jì)算，并將計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)進(jìn)行對比分析，最后對艦船在規(guī)則波中航行模擬的CFD方法進(jìn)行了初步探討。

1 數(shù)值波浪水池

本文的數(shù)值模擬是在長方形的三維數(shù)值波浪水池中進(jìn)行，其入口距離船首約1倍船長，出口距離船尾約2.5倍船長（距離隨著入射波的波長變化），側(cè)邊界距離舷側(cè)約2倍船長，上邊界距離波面約1倍船長，下邊界距離波面約2倍船長。如圖1所示，上部為空氣，下部為水，兩者間的交界面即為自由波面，波浪沿著x軸的負(fù)方向傳播。

圖1 三維數(shù)值波浪水池Fig.1 3D numerical wave tank

1.1 控制方程及離散格式

數(shù)值波浪水池的數(shù)學(xué)模型以連續(xù)性方程和N-S方程為控制方程：

式中：ui為流體質(zhì)點(diǎn)在i方向的速度分量；p為流體壓力；fi為質(zhì)量力；為流體密度；為動力粘性系數(shù)。

本文使用Realizable k-ε湍流模型，數(shù)值計(jì)算中采用有限體積法離散控制方程，選擇三維非定常分離隱式求解器和歐拉多項(xiàng)流模型，對流項(xiàng)采用二階迎風(fēng)格式，擴(kuò)散項(xiàng)為中心差分格式，采用半隱式方法（SIMPLE法）求解壓力耦合方程組；自由波面的追蹤使用VOF（Volume of Fluid，流體體積）方法處理，VOF 方法主要采用網(wǎng)格單元中流體的體積與網(wǎng)格總體積的比值函數(shù)來確定自由波面的位置和形狀，其方程為

式中：a1和a2分別為空氣相和水相的體積分?jǐn)?shù)，并定義aq=0.5處為自由波面。

1.2 造波與消波

目前常見的造波方法主要包括源造波方法和邊界模擬方法。本文采用在入口處給定波速函數(shù)的方法進(jìn)行造波，根據(jù)線性理論，考慮無限水深的情況，規(guī)則波的波面方程和速度場可表示為：

通過在離散的N-S方程中添加源項(xiàng)，可以實(shí)現(xiàn)強(qiáng)迫消波[13]。

由圖2（b）可見，阻尼消波方法是逐漸消除靠近出口的波浪，而強(qiáng)迫消波方法是將強(qiáng)迫區(qū)域內(nèi)的波浪強(qiáng)迫保持原來的輪廓；強(qiáng)迫消波方法的邊界條件設(shè)置與阻尼消波有所不同：一般入口、出口以及側(cè)面的邊界條件都設(shè)置為速度入口，而頂部邊界條件設(shè)置為壓力出口。

圖2 強(qiáng)迫消波技術(shù)Fig.2 Wave-forcing technology

1.3 波浪模擬結(jié)果與分析

進(jìn)行波浪中船模航行模擬前，需要在數(shù)值水池中實(shí)現(xiàn)波浪的生成與傳播。本文共模擬了5個(gè)不同波長的規(guī)則波，包括短波長（Case 1，Case 2）、中波長（Case 3）、長波長（Case 4，Case5），波高比波長（H/λ）固定為1/60，水深認(rèn)為是無限的，各方案波浪要素如表1所示。

本文是在三維數(shù)值水池中對規(guī)則波進(jìn)行模擬，為避免數(shù)值耗散引起波浪幅值的沿程衰減較大，在波高范圍內(nèi)需要設(shè)置足夠的網(wǎng)格數(shù)，一般設(shè)置10～20個(gè)網(wǎng)格單元，而在波長范圍內(nèi)網(wǎng)格單元一般不少于80個(gè)，網(wǎng)格單元尺寸的比例取為：；通過在數(shù)值波浪水池中設(shè)置波高監(jiān)測點(diǎn)來獲得5個(gè)不同波長入射波的波面變化時(shí)間歷程，從圖3可以看出，所模擬的規(guī)則波波幅與一階Stokes波理論波幅基本吻合。

2 規(guī)則波中船模迎浪航行數(shù)值模擬

2.1 帶舵的 KCS 船模算例描述

本文采用CFD方法計(jì)算了不同波長的規(guī)則波中帶舵的KCS船模迎浪增阻與運(yùn)動響應(yīng)（縱搖和垂蕩），并將計(jì)算結(jié)果與試驗(yàn)結(jié)果[14]進(jìn)行對比分析。由于整個(gè)流場和船體幾何關(guān)于船舶中縱剖面對稱，因此本文計(jì)算域僅采用流場的一半進(jìn)行模擬。

KCS船是由韓國船舶與海洋工程研究所（KRISO）設(shè)計(jì)的3 600 TEU集裝箱船，該KCS船模（包括舵）主尺度參數(shù)如表2所示，船體幾何如圖4所示。

表1 各方案波浪要素Tab.1 Wave conditions of each case

2.2 重疊網(wǎng)格方法

本文采用重疊網(wǎng)格方法求解船模在波浪中的運(yùn)動響應(yīng)。重疊網(wǎng)格（overset mesh，又稱為嵌套網(wǎng)格），是一種區(qū)域分割與網(wǎng)格組合的策略，其包含2種區(qū)域：背景區(qū)域和重疊區(qū)域，如圖5（a）所示，整個(gè)計(jì)算域?yàn)楸尘皡^(qū)域，嵌套在其中的小區(qū)域?yàn)橹丿B區(qū)域，重疊區(qū)域包裹著船體并隨船體同步運(yùn)動；背景區(qū)域內(nèi)產(chǎn)生的網(wǎng)格稱為背景網(wǎng)格，同樣，重疊區(qū)域內(nèi)產(chǎn)生的網(wǎng)格稱為重疊網(wǎng)格，如圖5（b）所示；該套網(wǎng)格包含2種類型的單元：有效單元和無效單元。有效單元是指參與離散求解控制方程的網(wǎng)格單元，而不參與計(jì)算的即為無效單元。沿著有效單元的第1層無效單元被稱為受者單元，而貢獻(xiàn)單元和受者單元相鄰卻位于不同區(qū)域。2套網(wǎng)格之間通過形成交界面（Interface）進(jìn)行信息交換，通過受者單元和貢獻(xiàn)單元之間插值來完成，一般采用拉格朗日插值的思想，進(jìn)行線性插值，具有以下形式[13]：

圖3 各方案監(jiān)測點(diǎn)波高的時(shí)間歷程Fig.3 Time history of wave elevation for various wavelengths

表2 KCS 船模的主尺度參數(shù)Tab.2 Main dimensions of KCS model

圖4 帶舵的 KCS 船型計(jì)算模型Fig.4 KCS ship computational model with rudder

圖5 重疊網(wǎng)格方法Fig.5 Overset mesh method

2.3 船體運(yùn)動方程

船舶在規(guī)則波中迎浪航行時(shí)，受到波浪的擾動，船舶將產(chǎn)生搖蕩運(yùn)動，而其中縱搖和垂蕩運(yùn)動對波浪增阻的影響較大，因此本文僅考慮縱搖和垂蕩運(yùn)動，其運(yùn)動方程為：

2.4 網(wǎng)格生成

本文采用軟件中切割體網(wǎng)格技術(shù)和棱柱層網(wǎng)格技術(shù)進(jìn)行網(wǎng)格生成，在船體近壁面處設(shè)置有5層棱柱層網(wǎng)格，平均y+取60左右，而舵附近僅設(shè)置2層，為了減少計(jì)算量，甲板可以不設(shè)置棱柱層網(wǎng)格；對于船體曲率變化大、流場變化比較劇烈的局部區(qū)域，例如球首、球尾和舵等，需要進(jìn)行局部的網(wǎng)格加密；在船身周圍和開爾文興波區(qū)域內(nèi)物理量隨著空間位置的變化梯度較大，因此在這一范圍內(nèi)網(wǎng)格要精細(xì)一些，而遠(yuǎn)場的網(wǎng)格可以稀疏一些；而尾部消波區(qū)的網(wǎng)格設(shè)置為沿縱向逐漸變稀，這樣可以起到數(shù)值消波的作用；由于采用重疊網(wǎng)格技術(shù)實(shí)現(xiàn)船模在波浪中的縱搖和垂蕩運(yùn)動，背景區(qū)域與重疊區(qū)域網(wǎng)格密度要相互匹配，以保證兩區(qū)域網(wǎng)格之間的信息交換更加精確。圖6為各區(qū)域網(wǎng)格和船體表面網(wǎng)格。

2.5 數(shù)值計(jì)算結(jié)果與分析

為了方便結(jié)果分析，對試驗(yàn)結(jié)果進(jìn)行無量綱化處理，船舶在規(guī)則波中迎浪航行的垂蕩與縱搖方程[15]:

無量綱化垂蕩與縱搖幅值表達(dá)式：

圖6 區(qū)域與船體網(wǎng)格劃分Fig.6 Meshing of overlap and hull

式中：R為船舶阻力；S為濕表面積；U為船速。

船舶在波浪中增加的阻力系數(shù)表示為:

如圖7所示，由于計(jì)算工況較多，本文僅顯示了當(dāng)船體運(yùn)動穩(wěn)定后Case 3三個(gè)周期的阻力系數(shù)、垂蕩與縱搖時(shí)歷曲線。

圖7 Case 3 三個(gè)周期的阻力系數(shù)、垂蕩與縱搖時(shí)間歷程Fig.7 Time history for resistance coefficient, heave and pitch motions over 3 wave periods of Case 3

波浪中船舶阻力的平均值通過對阻力時(shí)間歷程取平均得到，而靜水中船舶阻力可以通過靜水中自由船模阻力與運(yùn)動響應(yīng)數(shù)值計(jì)算獲得，Case 0（靜水中工況）的計(jì)算結(jié)果如表3所示。

圖8給出了不同波長的規(guī)則波中KCS船模的阻力系數(shù)、垂蕩和縱搖的時(shí)歷曲線，并將CFD計(jì)算結(jié)果與船模試驗(yàn)（EFD）數(shù)據(jù)進(jìn)行對比。由圖8可見，在短波長（Case 1，Case 2）中，垂蕩運(yùn)動響應(yīng)與實(shí)驗(yàn)數(shù)據(jù)略有不同，特別是在Case 2中，阻力計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)差異較大；中、長波長（Case 3，Case 4，Case 5）中的運(yùn)動響應(yīng)與實(shí)驗(yàn)數(shù)據(jù)相對一致。由于Case 3模型試驗(yàn)的阻力值存在較大變化，因此其阻力系數(shù)取為平均值。數(shù)值計(jì)算與試驗(yàn)結(jié)果存在差異的潛在原因可能是模型試驗(yàn)存在不確定性，模型試驗(yàn)與數(shù)值計(jì)算存在不同初始條件，以及反射波的影響等。

圖9給出了KCS船模波浪增阻系數(shù)以及垂蕩、縱搖傳遞函數(shù)曲線，并與試驗(yàn)數(shù)據(jù)進(jìn)行對比。由于波高比波長（H/λ）固定為1/60，波幅很小，因此船模阻力增加相對較小。從圖9可以看出，在短波范圍內(nèi)波浪增阻系數(shù)較小，隨著波長增大到稍大于船長（λ/L=1.15附近），波浪增阻出現(xiàn)峰值，而對于長波長的規(guī)則波，其垂蕩運(yùn)動的幅值接近于相應(yīng)規(guī)則波的波幅，此時(shí)船舶處于“隨波逐流”狀態(tài)，盡管運(yùn)動響應(yīng)較大，而增加的阻力卻較小?？偟膩碇v，靜水中阻力計(jì)算的差異和生成的波幅差異增加了波浪增阻的計(jì)算誤差，而垂蕩與縱搖傳遞函數(shù)隨波長與船長比（λ/L）的變化與試驗(yàn)結(jié)果具有相同的變化趨勢。

圖10給出了計(jì)算穩(wěn)定后KCS船模在靜水和5個(gè)不同波長的規(guī)則波中的瞬時(shí)自由波面波形圖，從圖中可以看出波浪與船行波的相互作用隨波長與船長比（λ/L）的變化。

表3 靜水中船模阻力系數(shù)與垂蕩、縱搖幅值Tab.3 Resistance coefficient and heave and pitch motions in calm water

3 結(jié) 語

本文基于CFD軟件建立了數(shù)值波浪水池，生成的規(guī)則波波幅與理論波幅吻合良好，并結(jié)合重疊網(wǎng)格技術(shù)與強(qiáng)迫消波技術(shù)對帶舵的KCS船模在靜水和5個(gè)不同波長的規(guī)則波中迎浪航行進(jìn)行了數(shù)值模擬，計(jì)算了其在靜水和波浪中的阻力、垂蕩與縱搖運(yùn)動響應(yīng)。本文數(shù)值計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)總體上吻合良好，說明了CFD方法預(yù)報(bào)船舶在靜水或波浪中的阻力和運(yùn)動的可靠性，并得出如下結(jié)論：

1）為避免數(shù)值耗散引起波浪幅值的沿程衰減較大，在波高范圍內(nèi)可以設(shè)置10～20個(gè)網(wǎng)格單元，而在波長范圍內(nèi)網(wǎng)格單元一般不少于80個(gè)，網(wǎng)格單元尺寸的比例（）可以取為1/2或1/4；

2）研究表明重疊網(wǎng)格技術(shù)可以準(zhǔn)確處理波浪中船模的運(yùn)動問題，能夠獲得較為精確的預(yù)報(bào)結(jié)果；

3）所采用的強(qiáng)迫消波技術(shù)獲得了較好的效果，能夠有效消除邊界的波浪反射，縮短了模擬時(shí)間，同時(shí)可以實(shí)現(xiàn)計(jì)算域入口處的消波。

圖8 各方案阻力系數(shù)、垂蕩和縱搖的時(shí)間歷程Fig.8 Time history for resistance coefficient, heave and pitch motions of each case

圖9 波浪增阻系數(shù)、垂蕩與縱搖傳遞函數(shù)Fig.9 Added resistance coefficient, heave and pitch motions in head waves

圖10 各方案瞬時(shí)自由波面波形圖Fig.10 Instantaneous wave figure of each case