熊耀宇 段 宏 姜治芳 陳 立

中國艦船研究設(shè)計中心，湖北武漢 430064

數(shù)值模擬現(xiàn)代水面艦船帶自由面的湍流繞流場

熊耀宇 段 宏 姜治芳 陳 立

中國艦船研究設(shè)計中心，湖北武漢 430064

為了研究帶自由面的船舶湍流繞流場，選擇了一艘ITTC推薦的公開船模作為模擬對象，數(shù)值求解RANS方程。此模型為帶聲吶導(dǎo)流罩和方形尾封板的復(fù)雜水面艦船。采用Ogrid塊拓撲算法生成質(zhì)量較高的純六面體多塊結(jié)構(gòu)化網(wǎng)格，計算中采用RNG k－ε湍流模式和標(biāo)準壁面函數(shù)，并使用VOF算法來捕捉自由面?？臻g離散采用QUICK格式，壓力－速度解耦采用PISO算法。將阻力和波形的數(shù)值結(jié)果與實驗數(shù)據(jù)相比較，總阻力系數(shù)誤差約2%，船側(cè)波形、船首自由面吻合良好，顯示了CFD方法在船舶水動力學(xué)中預(yù)測帶自由面湍流繞流場的有效性。

船舶阻力；自由面；波形；CFD；Ogird；湍流；水面艦船

1 引言

由于存在著不斷增加的設(shè)計目標(biāo)，現(xiàn)代水面艦船在設(shè)計過程中遇到了越來越多的困難［1］，比如艦首型線不僅要考慮減小阻力，還需降低噪聲，而且要有用于聲吶設(shè)備安裝的空間構(gòu)型。

當(dāng)前，在船舶水動力學(xué)中，CFD已經(jīng)成為了不可或缺的工具，不僅作為學(xué)術(shù)研究的手段，還可用做常規(guī)的設(shè)計工具［2］。其重要性不斷增加，而帶自由面的船舶湍流繞流場的數(shù)值模擬是其中的一個重要領(lǐng)域。

在船舶水動力學(xué)的數(shù)值計算中，CFD主要以如下幾種方式得到應(yīng)用：忽略粘性的流場計算［3］，在船模［4］或者實船［5］尺度下的粘性流場模擬以及無粘／有粘混合的分區(qū)計算方法［6］。無粘假設(shè)至今還被經(jīng)常使用，不過RANS、DES和LES正在船舶粘性流場計算中發(fā)揮越來越大的作用。

本文選擇了公開的DTMB 5415艦船模型用于CFD方法驗證，此模型是ITTC推薦的艦船類型中的唯一公共數(shù)值模擬平臺［7］。模型帶有方形尾封板和用作聲吶導(dǎo)流罩的球首。三家研究機構(gòu)的船模實驗數(shù)據(jù)得到了ITTC的認可和采用，分別是 ：DTMB 5415，IIHR 5512 和 INSEAN 2340。DTMB 5415 是縮尺比為 1／24.8 的模型，INSEAN 2340的幾何構(gòu)型與其完全相同，而IIHR 5512的縮尺比為 1／46.6。ISEAN 2340 模型的示意見圖 1，幾何尺寸見表1。

表1 艦船模型的主尺度

INSEAN與IIHR合作出具了一份公開的官方技術(shù)報告［8］，本文選其實驗數(shù)據(jù)用作CFD結(jié)果驗證。

本文中生成的多塊結(jié)構(gòu)化純六面體網(wǎng)格具有良好的正交性，用其求解瞬態(tài)RANS方程來模擬船模繞流場，并與實驗數(shù)據(jù)進行比較，主要關(guān)注點為阻力、船側(cè)波形和首部自由面。

2 理論模型

2.1 控制方程

普適的連續(xù)性方程和N－S方程的張量形式［9］如下：

在船舶水動力學(xué)中，為了簡化上述方程，采用不可壓縮流體近似，并假設(shè)流體物性參數(shù)均勻且不變（密度、粘性系數(shù)等）。再假定船舶繞流場滿足各態(tài)遍歷假設(shè)，使用雷諾分解。在笛卡爾坐標(biāo)系下，式（1）、式（2）簡化為：

式中，ui，ui′分別為時均速度和脈動速度的i方向分量（i＝ 1，2，3）；為雷諾應(yīng)力張量；下標(biāo)滿足愛因斯坦求和約定。

2.2 湍流模式

式（4）引入的雷諾應(yīng)力張量有6個新的未知量，數(shù)值求解中必須加以封閉。最常用的封閉模式是Boussinesq的湍流渦粘假設(shè)［10］，忽略雷諾應(yīng)力張量的歷史效應(yīng)，將其與速度導(dǎo)數(shù)張量的關(guān)系近似描述為：

方程中的經(jīng)驗常數(shù)取值見表2。

表2 RNG k－ε湍流模式的經(jīng)驗常數(shù)取值

方程中各項的詳細物理意義和表達式參看文獻［11-12］。

2.3VOF 方法

為解決流場計算中的自由邊界問題，Hirt和Nichols提出了VOF方法［13］。船舶繞流場的自由面是典型的兩相流態(tài)自由邊界。當(dāng)流體增多一相時，網(wǎng)格的每一個控制體積就會增加一個變量用于描述此相流體的體積分數(shù)。用αn表示第n相的體積分數(shù)，對于網(wǎng)格內(nèi)的任意控制體積，可能存在以下3種情況：

αn＝0：控制體積內(nèi)沒有第n相的流體。

0＜αn＜1：控制體積內(nèi)包含第n相流體與其他相流體的分界面。

αn＝1：控制體積內(nèi)充滿第n相的流體。

當(dāng)流域中各網(wǎng)格內(nèi)的相體積分數(shù)確定后，對位于兩相交界面的控制體采用插值算法來確定自由面。

3 數(shù)值模擬

3.1 流域

建立笛卡爾坐標(biāo)系如下：以基平面、中站面、中線面的交點作為原點O，X軸指向船尾為正，Y軸指向右舷為正，Z軸指向上方為正。

忽略尾流中可能存在的不對稱渦街，假設(shè)船舶繞流場以中縱剖面為鏡面對稱。將中縱剖面作為流域的一個邊界，以此減少流域體積和網(wǎng)格數(shù)量。

流域尺寸的設(shè)置參考文獻［14］，并適當(dāng)延長出流面和增高頂面以擴大計算流域。長方體流域的邊界位置參看圖2和表3。

表3 流域邊界的位置

3.2 網(wǎng)格

ITTC阻力委員會的最新報告中推薦使用純六面體的結(jié)構(gòu)化網(wǎng)絡(luò)［2］，對上述流域使用專業(yè)網(wǎng)格軟件ICEM CFD［12］生成多塊結(jié)構(gòu)化網(wǎng)格，得到的純六面體網(wǎng)格包含了大約64.8萬個有限控制體和67.8萬個節(jié)點。

使用自上而下的剖分方法，將一個單塊剖分為大約100個分塊。在船表曲面形狀復(fù)雜的地方使用Ogrid算法以提高網(wǎng)格質(zhì)量，Ogrid拓撲分為內(nèi)Ogrid和外Ogrid兩種，其主要思想都是在一個塊的四周創(chuàng)造數(shù)個環(huán)繞的塊拓撲結(jié)構(gòu)，詳細的拓撲理論參看文獻［15］。

在船表面所有地方，沿法線創(chuàng)造外Ogrid塊作為邊界層網(wǎng)格。對于首部導(dǎo)流罩區(qū)域，由于曲率過大，需要反復(fù)嘗試調(diào)整節(jié)點位置和外Ogrid的形狀以獲得較好的網(wǎng)格質(zhì)量，見圖3。

在尾封板區(qū)域和上甲板前部區(qū)域，使用內(nèi)Ogrid算法處理曲率過大的曲線，網(wǎng)格見圖4和圖5。

使用正交性算法 Determinant 2×2 ×2［14］檢查網(wǎng)格質(zhì)量，以此作為網(wǎng)格調(diào)整優(yōu)劣的參考。參看圖6，全部網(wǎng)格控制體的正交因子中，最小為0.262，大于 0.2，適合求解器使用［15］。 正交因子小于 0.5的網(wǎng)格控制體總數(shù)大約為100個，約占總數(shù)（64.8萬）的 0.02%，其余 99.98%的網(wǎng)格控制體的正交因子大于0.5，顯示了良好的網(wǎng)格質(zhì)量。

由于在邊界層內(nèi)使用壁面函數(shù)法，邊界層網(wǎng)格內(nèi)首層節(jié)點需要布置在對數(shù)律層。本文先采用平板湍流邊界層近似，求出首層網(wǎng)格節(jié)點的量級范圍，約為mm級。在流場計算結(jié)束后，再檢查船表面各處y＋的值，若不符合壁面函數(shù)法要求，就重新調(diào)整網(wǎng)格。邊界層網(wǎng)格首層節(jié)點的法向高度取1.5 mm，并以1.1的等比因子遞增15層。邊界層網(wǎng)格的示意可以參看圖3和圖4，其合理性需要以流場的計算結(jié)果作為判據(jù)，參看下一節(jié)。

3.3 數(shù)值格式

使用PISO算法來處理壓力—速度解耦；使用基于節(jié)點的格林高斯方法來處理梯度項；使用QUICK格式進行空間離散；使用二階隱式格式進行時間離散。詳細的數(shù)學(xué)機理參看文獻［11，16］。

3.4 邊界條件

在入流面，使用速度入口邊界條件給定來流速度和湍流參數(shù)，布置兩相流體積分數(shù)的分布。在出流面，使用壓力出口邊界條件編寫UDF來描述靜壓貫穿兩相并隨高度變化的性質(zhì)。在頂部面、底部面和右側(cè)面，使用運動壁面邊界條件給定其平移速度為來流速度。在中縱面，給定對稱邊界條件假定所有的物理量在此面為零通量。

3.5 初始條件

按照來流速度初始化整個流域，以水線面為界限給定兩相流體積分數(shù)的分布并分區(qū)賦值湍流參數(shù)。

由于對時間項采用隱式算法，時間步長的給定并沒有嚴格的判據(jù)。本文以柯朗數(shù)等于1為參照并經(jīng)過多次試算，給出合適的時間步長。在迭代計算的初始階段，Δt＝ 0.005 s， 待流場穩(wěn)定后，再適當(dāng)增加時間步長以節(jié)約計算的硬件資源。

4 計算結(jié)果及分析

4.1 計算的實驗工況

選取的用于CFD驗證的實驗工況［8］見表3，此工況下，INSEAN重復(fù)做過10次測試，實驗值為10次重復(fù)所取的平均值。

表3 選取的實驗工況

4.2 阻力系數(shù)的數(shù)值結(jié)果

使用 Fluent軟件［16］，數(shù)值求解非穩(wěn)態(tài) RANS方程，在數(shù)值時間約600 s后，總阻力系數(shù)平穩(wěn)收斂，適當(dāng)增加時間步長，繼續(xù)計算至1 100 s。采用八核并行計算，單核主頻2.6 GHz，耗費的計算資源約為500 CPU hours。得到的總阻力系數(shù)收斂曲線見圖7。

求解瞬態(tài)RANSE時，由于自由面的存在，總阻力系數(shù)呈現(xiàn)周期性的振蕩，并逐漸收斂到穩(wěn)定值。

計算得到的阻力系數(shù)與實驗數(shù)據(jù)比較見表4。

表4 阻力系數(shù)的數(shù)值模擬結(jié)果與實驗值的比較

計算得到的總阻力系數(shù)與實驗值的偏差約為2.0%，具有很高的精度，顯示了CFD方法的有效性和應(yīng)用前景。計算得到的摩擦阻力系數(shù)與ITTC 57公式的偏差約為4.3%，在工程的合理范圍內(nèi)。計算得到的壓差阻力系數(shù)與實驗的剩余阻力系數(shù)偏差約為－3.2%。

4.3 波形

自由面采用體積分數(shù)等值面作為分界面判據(jù)，取水的體積分數(shù)為50%，比較船側(cè)波形曲線和首部波形自由面。

船側(cè)波形曲線的計算值與實驗值的對比如圖8所示。首部自由面波形的計算、實驗照片分別如圖9、圖10所示。

在船側(cè)，波形曲線的計算值與實驗值符合較好。首部的興波形態(tài)能準確地模擬，首部橫波波峰值偏差5%，其波形趨勢良好。尾部橫波的前半段（波高上升區(qū)域）能在曲線上定性地展現(xiàn)。船模表面其余各處的波形與實驗值近似相同，數(shù)值結(jié)果清晰地捕捉到了兩個波谷，與實驗數(shù)據(jù)吻合。

4.4 y＋云圖

壁面函數(shù)法忽略粘性底層的求解，通過工程公式來聯(lián)系固壁處和邊界層核心區(qū)的物理參量，因此壁面法線的第一個網(wǎng)格節(jié)點的高度有尺寸要求。

詳細理論參看文獻［10，16］。

利用云圖檢查船表曲面各處的y＋，其取值如圖11所示。

在自由面以下的水相區(qū)域，30≤y＋≤179，參看公式（8），y＋分布滿足壁面函數(shù)法的要求。

對于空氣相，有部分網(wǎng)格 y＋≤11.5。高速行駛的船舶，其總阻力中的空氣阻力比例很小，可忽略空氣相網(wǎng)格的偏差所其帶來的總阻力計算誤差。

y＋云圖表明了本文邊界層網(wǎng)格的合理性。

5 結(jié) 論

為了對現(xiàn)代水面艦船的繞流場進行數(shù)值模擬和分析，本文選取了具有復(fù)雜曲面的ITTC DTMB 5415模型，利用RANS方程，求解高雷諾數(shù)下帶自由面的湍流繞流場，計算結(jié)果與實驗數(shù)據(jù)符合良好，顯示了CFD方法在水動力學(xué)研究和船型設(shè)計領(lǐng)域的準確性和實用性，得到的結(jié)論如下：

1）數(shù)值模擬的結(jié)果與實驗數(shù)據(jù)吻合良好，表明CFD可以為艦船設(shè)計人員提供可靠的設(shè)計參考，降低模型實驗的成本，以加快新型艦船的研發(fā)進程。

2）求解RANS方程所得到的總阻力系數(shù)可以精確到2%，也能較準確地得到摩擦阻力（4.3%）分量和壓差阻力（－3.2%）分量，表明 RANS方法在阻力計算方面的有效性。

3）采用VOF算法，計算得到的波形和自由面與實驗吻合良好，表明此方法能比較精確地構(gòu)造出自由界面的位置和形狀。

4）多塊結(jié)構(gòu)化網(wǎng)格的質(zhì)量需要通過反復(fù)調(diào)整節(jié)點位置來得到優(yōu)化，Ogrid塊拓撲算法具有強大的優(yōu)勢，能得到正交性很高的網(wǎng)格。

［1］CAMPANA E F， PERI D， TAHARAK Y，et al.Shape optimization in ship hydrodynamics using computational fluid dynamics［J］.Computer methods in applied mechanics and engineering，2006，196：634－651.

［2］International Towing Tank Conference ［R］.Proceedings of 25th ITTC.Fukuoka，2008.

［3］NOBLESSE F， YANG C.Alternative representations of steady potential flow about a ship ［C］//Proceeding of the 7th International Conference on Hydrodynamics （ICHD）.Italy，2006.

［4］WILSON W，F(xiàn)U T C，F(xiàn)ULLERTON A，et al.The measured and predicted wave field of model 5365：an evaluation of current CFD capability［C］//Proceeding of 26th Symposium on Naval Hydrodynamics.Italy，2006.

［5］BHUSHAN S，XING T，CARRICA P，et al.Model and full－scale URANS ／DES simulations for athena R ／V resistance， powering and motions［C］//Proceeding of 9th International Conference on Numerical Ship Hydrodyanmics.Ann Arbor， MI， USA，2007.

［6］HUAN J C， HUANG T T.Surface ship total resistance prediction based on a nonlinear free surface potential flow solver and a Reynolds－Averaged Navier－Stokes viscous correction［J］.Journal of Ship Research，2007，51：47－64.

［7］International Towing Tank Conference.CFD，Resistance and flow benchmark database for CFD validation for resistance and propulsion［R］，1999.

［8］OLIVIERI A，PISTANI F，AVANZINI A，et al.Towing tank experiments of resistance， sinkage and trim， boundary layer，wake and free surface flow around a naval combatant INSEAN 2340 Model ［R］.IIHR Technical Report No.421，IIHR，2001.

［9］BATCHELOR G K.An introduction to fluid dynamics［M］.Cambridge： Cambridge University Press， 2000.

［10］STEPHEN B P.Turbulent flows ［M］.Cambridge：Cambridge University Press， 2000.

［11］VERSTEEG H K，MALALASEKERA W.An introduction to computationalfluid dynamics， the finite volume method［M］.Person Education Ltd，2007.

［12］ORSZAG S A， YAKHOT V， FLANNERY W S， et al.Renormalization group modeling and turbulence simulations［C］//International Conference on Near－Wall Turbulent Flows， Tempe， Arizona，1993.

［13］HIRT C W，NICHOLS B D.Volume of fluid （VOF）method for the dynamics of free boundaries［J］.Journal of computational Physics， 1981， 39： 201－225.

［14］ZHANG Z R，ZHAO F，LI B Q.Numerical calculation of viscous free－surface flow about ship hull ［J］.Journal of Ship Mechanics，2002（6）：10－17.

［15］ANSYS Inc.ANSYS ICEM CFD Manual［S］， 2009.

［16］ANSYS Inc.ANSYS FLUENT Guide［S］， 2009.

Numerical Simulation of Modern Surface Ship in Free Surface Turbulent Flow

Xiong Yao－yu Duan Hong Jiang Zhi－fang Chen Li
China Ship Development and Design Center， Wuhan 430064， China

A public ITTC recommended complex modern surface ship model with sonar dome and transom stern was chosen for numerical simulation of free surface turbulent flow.Pure hexahedral multiblock structured grids with fine quality were generated by using Ogrid topology.Transient RANSE were numerically calculated and the RNG k－ε turbulence model with standard wall functions was employed.The VOF algorithm was applied to track free surface.QUICK scheme was used for spatial discretization and PISO scheme for pressure－velocity coupling.Computational results including resistance coefficients and wave profile were compared to experimental data with good agreement.The total resistance coefficient has a bias of two percentage，wave profile and bow free surface accord the data with a good precision，showing that CFD method can effectively predict the free surface turbulent flow in ship hydrodynamics.

ship resistance；free surface；wave profile；turbulent flow；CFD；Ogrid；surface ship

U661.1

A

1673－3185（2010）06－16－05

10.3969/j.issn.1673-3185.2010.06.004

2010-02-19

熊耀宇（1985－），男，碩士研究生。研究方向：船舶水動力學(xué)。E-mail：xiong－yaoyu＠ hotmail.com

段 宏（1954－），男，研究員。研究方向：艦船總體研究設(shè)計