基于重疊網格方法的船舶迎浪增阻與運動數值預報

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

(1. 江蘇科技大學 船舶與海洋工程學院，江蘇 鎮江 212003；2. 江蘇現代造船技術有限公司，江蘇 鎮江 212003；3. 廣州船舶及海洋工程設計研究院，廣東 廣州 510250)

0 引 言

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

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

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

1 數值波浪水池

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

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

1.1 控制方程及離散格式

數值波浪水池的數學模型以連續性方程和N-S方程為控制方程：

式中：ui為流體質點在i方向的速度分量；p為流體壓力；fi為質量力；為流體密度；為動力粘性系數。

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

式中：a1和a2分別為空氣相和水相的體積分數，并定義aq=0.5處為自由波面。

1.2 造波與消波

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

通過在離散的N-S方程中添加源項，可以實現強迫消波[13]。

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

圖2 強迫消波技術Fig.2 Wave-forcing technology

1.3 波浪模擬結果與分析

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

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

2 規則波中船模迎浪航行數值模擬

2.1 帶舵的 KCS 船模算例描述

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

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

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

2.2 重疊網格方法

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

圖3 各方案監測點波高的時間歷程Fig.3 Time history of wave elevation for various wavelengths

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

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

圖5 重疊網格方法Fig.5 Overset mesh method

2.3 船體運動方程

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

2.4 網格生成

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

2.5 數值計算結果與分析

為了方便結果分析，對試驗結果進行無量綱化處理，船舶在規則波中迎浪航行的垂蕩與縱搖方程[15]:

無量綱化垂蕩與縱搖幅值表達式：

圖6 區域與船體網格劃分Fig.6 Meshing of overlap and hull

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

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

如圖7所示，由于計算工況較多，本文僅顯示了當船體運動穩定后Case 3三個周期的阻力系數、垂蕩與縱搖時歷曲線。

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

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

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

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

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

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

3 結 語

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

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

2）研究表明重疊網格技術可以準確處理波浪中船模的運動問題，能夠獲得較為精確的預報結果；

3）所采用的強迫消波技術獲得了較好的效果，能夠有效消除邊界的波浪反射，縮短了模擬時間，同時可以實現計算域入口處的消波。

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

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

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