基于RANS法和邊界層理論預報三維船體阻力

胡俊明，李鐵驪，林 焰，徐雪鋒，徐利剛

（1.大連理工大學 船舶工程學院，遼寧 大連 116024；2.海軍裝備部直屬工作部，北京 100841；3.海軍駐無錫地區軍事代表室，江蘇 無錫 214061）

0 引 言

船舶阻力問題是船舶行業重要的研究課題，如何快速精確預報船體阻力成為問題研究的關鍵。模型試驗法在船舶阻力預報中一直占據著主導地位，其預報結果精確可靠，實用性強，但其經濟性較差且無法完全模擬復雜的實際海況，存在尺度效應[1]。理論計算方法由來已久，由于船體復雜的三維曲面，水的粘性和自由面的非線性，使理論求解相當困難，目前求解摩擦阻力沿用傅汝德的相當平板假定，而興波阻力理論由1898年Michell薄船理論的出現，奠定線性興波理論的基礎，之后的慢船理論Dawson法的產生，借助攝動理論使線性興波向非線性興波發展成為一種可能[2-4]。隨著技術進步，計算機性能更新和完善，計算流體力學（CFD）的興起，使數值求解復雜曲面粘性自由面繞流場成為現實，計及波浪破碎、渦分離、抨擊等非線性現象，且計算結果與船模試驗結果一致，真實可靠，經濟性好，可用范圍廣，成為研究船舶阻力性能的重要手段[5-7]。

本文應用二因次RANS法研究船舶阻力，基于粘性流理論N-S方程求解自由面數值繞流場，重點研究了船體自由面繞流模型的網格劃分方法，應用ICEM網格處理軟件處理內外域交接面，實現交界面的網格重構，以提高內外域不同網格類型之間的數據交換效率，縮短計算時間，同時應用UDF控制壓強出口邊界處壓強隨水深變化以提高數值計算穩定性，提高計算效率。數值計算中基于邊界層理論考慮邊界層第一層網格高度對數值計算結果的影響，重在研究其對船體各阻力分量的影響。基于此方法所得數值結果與實驗流體力學法（EFD）做比較分析，通過其計算精度、船體阻力各分量值和興波波形的比較得到其變化規律，結果表明，此方法預報船體阻力其數值精度滿足工程需要，具有較強的實用性。

1 船型特征及模型試驗

文中以總長為33.2 m的漁船為實船原型，進行船模試驗，實船與模型的縮尺比為10，且實船與模型在壓載狀態下主尺度如表1所示，表中LOA為總長，LWL為設計水線長，B為船寬，T為設計吃水，CB為方形系數。應用三因次EFD法對船模進行低速流試驗，根據傅汝德數Fn=0.1～0.2范圍內的試驗結果，采用普魯哈斯卡法確定形狀因子1+K的值，應用1957-ITTC平板摩擦阻力公式，計算船模平板摩擦阻力Rfm，根據測定的船模總阻力Rtm，計算船模剩余阻力Rrm=Rtm-Rfm、粘壓阻力Rpvm=KRfm、興波阻力及各阻力系數[8]。

表1 船型主尺度Tab.1 Principal dimensions of the real ship and the model

2 計算方法

2.1 控制方程

RANS二因次法基于粘性流理論N-S方程求解自由面船體阻力，考慮自由面非線性且流體不可壓縮，其流體的連續方程和雷諾平均N-S方程如下[9]：

式中：ρ是流體密度，t是時間，υ是流體運動粘性系數，p是流體壓力，fi是質量力，ui為時均速度；ui′為脈動速度為雷諾應力項。雷諾應力對于精確求解豐滿船型船體阻力至關重要，為使RANS方程封閉，應用RNG k-ε湍流模型，其輸運方程為[10]：

表2 RNG k-ε湍流模型常數Tab.2 RNG k-εturbulence model constant

2.2 邊界層理論

邊界層影響船體阻力的數值精度，特別是邊界層第一層網格高度的選擇和確定，決定數值計算結果的偏差大小。本文采用RNG k-ε模型，應用壁面函數處理粘性底層和過渡層，用一維數學模型代替求解三維N-S方程，以降低計算資源的使用，同時選擇合理邊界層第一層網格高度使其網格節點落于對數律區域內，避免網格節點落在粘性底層而導致數值計算結果的偏差。其邊界層第一層網格高度無量綱計算公式如下[11]：

式中：υ為流體的運動粘性系數，Δy為第一層網格高度，ρ為流體的密度，K和E為常數[9-10]。

3 計算模型建立與設置

3.1 模型及域選擇

因船型對稱且直航，受力對稱，文中所建模型為半船模型，應用Auto CAD可視化繪圖軟件和ICEM網格處理軟件聯合建立自由面網格繞流計算模型。因自由面興波對數值流場存在影響，為提高計算精度和效率應合理選擇計算域的大小，采用Hirt和Nichols提出的多相流模型方法（VOF）追蹤自由液面[12]。

應用水波理論中船行波波高ZA、波長λw與波速的關系選擇計算域的大小，其繞流計算模型為長方體。由公式確定水線以上計算域高度大于船行波最大波高，基于計算船行波波長，由于自由面影響，波面在船體前緣隆起，取進流段長度大于λw，尾部采用壓強出口以提高計算效率和數值結果的穩定，其長度應考慮首尾波系的影響，為實現其充分發展，文中去流段長度取為13.35 m，約大于6倍最大船行波波長。根據開爾文波系盡可能減小寬度區域邊界對船行波的反射，由計算域長度和開爾文角選取合適的區域寬度，計算域深度大于λw時認為船行波對繞流場影響甚微[13]。由模型航速范圍1.3～1.8 m/s，估算船行波波高范圍0.043 1～0.082 5 m，波長范圍1.081～2.073 m，其自由面網格繞流模型計算域尺寸劃分如表3所示。

表3 模型計算域劃分Tab.3 Size of the computational domain

3.2 網格劃分

網格計算模型考慮了邊界層和船體線型對流場的影響，文中借助網格處理軟件ICEM實現內外域交界面網格重構以合成整體域，以提高內外域不同網格類型之間的數據交換效率，縮短計算時間，在船艏艉區域進行網格加密處理，網格尺寸取為0.003L，船體表面最大網格尺寸為0.01L，體網格由船體表面按照一定的增長率向外擴展[14]，內域增長率為1.1，在外域中，采用分塊結構化網格，增長率為1.3，最大網格尺寸0.08L，考慮邊界層影響，在船體表面劃分8層邊界層，取邊界層為3 mm，增長率為1.1，在水線面附近進行網格加密，水面第一層高度為0.005L，增長率為1.2。其半船體模型如圖1所示，整個流場網格和局部邊界層網格劃分如圖2所示。

圖1 船體模型Fig.1 Model ship

圖2 網格劃分圖Fig.2 Mesh information

4 數值計算結果及分析

4.1 邊界層對數值結果的影響

本文考慮邊界層對數值計算結果的影響，對應實船設計航速，船模航速為1.5 m/s，傅汝德數Fn=0.274，邊界層船體表面第一層網格高度分別取為2-6 mm，其數值計算結果如表4。由表中可以看出邊界層第一層網格高度Δy充分影響數值結果精度，因其第一層網格高度節點未落入粘性底層，隨著Δy值的減小，其摩擦阻力值略微減小，但影響甚小，其總阻力值增加，數值結果和試驗值總體相差不大，同時基于粘性Δy值的變化改變了流線的形狀，改變了船體壓力分布，影響興波阻力，直接影響剩余阻力值，對設計航速，當Δy=3 mm時，其數值計算結果與模型試驗結果的誤差最小，達到較好的吻合狀態。此外，本文分別應用船體表面邊界層第一層網格高度3 mm和4 mm，對船體模型進行各航速數值模擬計算，其計算結果詳見表5。由表5中可以看出，第一層網格高度Δy=3 mm的數值結果的精度優于Δy=4 mm，隨值增大，總阻力數值在高傅汝德數下其誤差有增大趨勢。結合表4和5，可以得到，對應任何一個航速模態，存在較佳的第一層網格高度Δy和值，使數值結果和試驗結果達到較佳吻合狀態，表中Rf為摩擦阻力、Rr為剩余阻力、Rt為模型試驗總阻力，單位為N，error為總阻力相對誤差。

表4 邊界層對設計航速數值結果的影響Tab.4 Influence of boundary layer to numerical results

表5 邊界層對各航速數值結果的影響Tab.5 Influence of boundary layer to numerical results under different speed

圖3 摩擦和剩余阻力系數曲線圖Fig.3 Curve of friction resistance coefficient and residual resistance coefficient

圖4 粘壓和興波阻力系數曲線圖Fig.4 Curve of viscous pressure resistance coefficient and wave resistance coefficient

4.2 阻力分析

應用邊界層船體表面第一層網格高度3 mm，考慮帶自由面粘性流計算，數值模擬不同傅汝德數下的船體阻力，所得結果與EFD法做分析比較，通過基于CFD二因次法所得摩擦阻力系數Cfm0與1957-ITTC平板摩擦阻力系數Cfm隨雷諾數的變化曲線圖3（a）中可以看出基于CFD帶自由面粘性流得到的摩擦阻力較平板摩擦阻力略大，由圖3（b）中可以看出基于CFD二因次法得到剩余阻力Crm0較EFD法剩余阻力Crm略小。基于三因次EFD法得到粘壓阻力系數、興波阻力系數的變化曲線如圖4所示，Cpvm隨雷諾數增加而減小，興波阻力系數Cwm隨傅汝德數增加而增大，與實際相符。

圖5 不同航速下的自由面波形圖Fig.5 Contour of wave pattern at various speed

4.3 波形及自由面分析

自由面波形反應了興波特性，為了分析研究其變化規律，本文基于Tecplot技術處理RANS法所得自由面波形與EFD法做分析比較，取三個航速點（Fn=0.274、0.292和0.329），其自由面波形和等值線圖如5所示。由圖5中可以看出隨著傅汝德數增加，船體首部波高等值線數值增大，波形等值線覆蓋區域增加，船中及尾部區域規律變化愈發明顯，船兩側受擾動區域增加，興波區域向外擴散，波峰數值增加且范圍擴大，船首尾肩部開爾文波系明顯且波峰后移，其興波波形完全符合開爾文波系特點，與EFD法興波圖形變化趨勢吻合一致，這為船型優化提供了強有力的證據。

5 結 論

本文基于RANS法對總長為33.2 m帶球鼻艏新型漁船進行船體阻力數值模擬，考慮粘性流自由面效應及邊界層的影響，結合EFD方法對數值結果和興波波形做分析比較，得到以下結論：

（1）應用ICEM處理網格繞流計算模型，實現內外域交接面網格重構，以提高內外域不同網格類型之間的數據交換效率，同時應用UDF控制壓強出口邊界處壓強隨水深變化以提高數值計算穩定性，縮短計算時間。

（2）基于邊界層理論研究粘性流帶自由面船體繞流問題，表明邊界層對數值計算結果精度存在重要影響，邊界層第一層網格高度節點未落入粘性底層時，在相同航速下，摩擦阻力值隨邊界層第一層網格高度的增加而略微增大，剩余阻力值卻明顯減少，總阻力在高傅汝德數下其數值誤差有增大趨勢，且邊界層第一層網格高度相同時，隨傅汝德數增加其船體總阻力誤差增大。

（3）基于邊界層數值結果分析得到對應任何一個航速模態，存在較佳的第一層網格高度Δy和y+值，使數值結果和EFD試驗結果達到較佳的吻合狀態。其RANS法數值結果得到的總阻力、摩擦阻力和剩余阻力系數與EFD法吻合較好，表明此方法用于船舶阻力性能預報的合理性、可行性和可靠性。

（4）自由面波形反應興波特性，隨傅汝德數增大船體首部波高越發明顯，船體兩側受擾動區域增加，船中及尾部區域興波向外擴散，船首尾肩部開爾文波系明顯且波峰后移，其興波圖形變化趨勢與EFD法吻合一致，為船型優化提供有力參考。

[1]倪崇本,朱仁傳,繆國平,范佘明.一種基于CFD的船舶總阻力預報方法[J].水動力學研究與進展,2010,25(5):579-586.Ni Chongben,Zhu Renchuan,Miao Guoping,Fan Sheming.A method for ship resistance prediction based on CFD computation[J].Chinese Journal of Hydrodynamics,2010,25(5):579-585.

[2]Dawson C W.A practical computer method for solving ship-wave problems[C].In:Proceeding of Second International Conference on Numerical Ship Hydrodynamics,1977:30-38.

[3]Suzukik,Iokamorin.Studies on minimization of wave making resistance based on Rankine source method[J].Kansai Society of Naval Architecture in Japan,1999,185:9-19.

[4]Raven H C.A practical nonlinear method for calculating ship wave making and wave resistance[C]//Preprints 19th Symp.on Naval Hydrodynamics.Seoul,Korea,1992:60-75.

[5]劉應中,張懷新,李誼樂,繆國平.21世紀的船舶性能計算和RANS方程[J].船舶力學,2001,5(5):66-84.Liu Yingzhong,Zhang Huaixin,Li Yile,Miao Guoping.Calculation of the ship performance and solving of RANSequations in the 21st century[J].Journal of Ship Mechanics,2001,5(5):66-84.

[6]Choi J E,Min K S,Kim J H,Lee SB,Seo H W.Resistance and propulsion characteristics of various commercial ships based on CFD results[J].Ocean Engineering,2010,37:549-566.

[7]Ciortan C,Wanderley J,Guedes Soares C.Turbulent free-surface flow around a wigley hull using the slightly compressible flow formulation[J].Ocean Engineering,2007,34:1383-1392.

[8]盛振邦,劉應中.船舶原理[M].上海:上海交通大學出版社,2004.Shen Zhenbang,Liu Yingzhong.Principles of Ship.[M].Shanghai:SJTU Press,2004.

[9]Versteeg H K,Malalasekera W.An introduction to computation fluid dynamics:The Finite Volume Method[M].Wigley,New York,1995.

[10]Yakhot V,Orzag SA.Renormalization group analysis of turbulence:Basic theory[J].Scient Comput,1986(1):3-11.

[11]Schlichting H.Boundary Layer Theory[M].McGrawhill,New York,1979.

[12]Hirt C W,Nichols B D.Volume of Fluid(VOF)Method for the dynamics of free boundary[J].Journal of Computational Physics,1981(39):201-225.

[13]顧溟宇,王文華,王言英.基于FLUENT的大型集裝箱船航速預報方法研究[J].船舶工程,2012,34(2):28-31.Gu Mingyu,Wang Wenhua,Wang Yanying.Studies of FLUENT-based speed prediction method for large containers[J].Ship Engineering,2012,34(2):28-31.

[14]朱德祥,張志榮,吳乘勝,趙 峰.船舶CFD不確定度分析及ITTC臨時規程的初步應用[J].水動力學研究與進展,2007,22(3):363-370.Zhu Dexiang,Zhang Zhirong,Wu Chengsheng,Zhao Feng.Uncertainty analysis in ship CFD and the primary application of ITTCprocedures[J].Journal of Hydrodynamics,2007,22(3):363-370.