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

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

0 引 言

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

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

1 船型特征及模型試驗

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

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

2 計算方法

2.1 控制方程

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

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

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

2.2 邊界層理論

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

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

3 計算模型建立與設(shè)置

3.1 模型及域選擇

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

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

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

3.2 網(wǎng)格劃分

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

圖1 船體模型Fig.1 Model ship

圖2 網(wǎng)格劃分圖Fig.2 Mesh information

4 數(shù)值計算結(jié)果及分析

4.1 邊界層對數(shù)值結(jié)果的影響

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

表4 邊界層對設(shè)計航速數(shù)值結(jié)果的影響Tab.4 Influence of boundary layer to numerical results

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

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

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

4.2 阻力分析

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

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

4.3 波形及自由面分析

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

5 結(jié) 論

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

（1）應(yīng)用ICEM處理網(wǎng)格繞流計算模型，實現(xiàn)內(nèi)外域交接面網(wǎng)格重構(gòu)，以提高內(nèi)外域不同網(wǎng)格類型之間的數(shù)據(jù)交換效率，同時應(yīng)用UDF控制壓強出口邊界處壓強隨水深變化以提高數(shù)值計算穩(wěn)定性，縮短計算時間。

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

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

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

[1]倪崇本,朱仁傳,繆國平,范佘明.一種基于CFD的船舶總阻力預(yù)報方法[J].水動力學(xué)研究與進展,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]劉應(yīng)中,張懷新,李誼樂,繆國平.21世紀(jì)的船舶性能計算和RANS方程[J].船舶力學(xué),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]盛振邦,劉應(yīng)中.船舶原理[M].上海:上海交通大學(xué)出版社,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的大型集裝箱船航速預(yù)報方法研究[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臨時規(guī)程的初步應(yīng)用[J].水動力學(xué)研究與進展,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.