小展弦比翼大攻角水動(dòng)力數(shù)值計(jì)算方法

周廣禮，姚朝幫

(海軍工程大學(xué) 艦船工程系，湖北 武漢430033)

0 引 言

小展弦比三維舵翼的水動(dòng)力預(yù)報(bào)是船舶操縱性預(yù)報(bào)的前提和基礎(chǔ)，在以往進(jìn)行舵初步設(shè)計(jì)時(shí)主要采用相關(guān)翼型已有的試驗(yàn)資料進(jìn)行展弦比換算［1］或基于勢(shì)流理論的方法來(lái)得到翼的水動(dòng)力性能［2］，然而采用展弦比換算方法所得計(jì)算結(jié)果誤差較大，基于勢(shì)流理論的方法并不能準(zhǔn)確預(yù)估大攻角翼的水動(dòng)力性能而且無(wú)法對(duì)幾何形狀較為復(fù)雜的舵(如半懸掛舵、襟翼舵)的水動(dòng)力進(jìn)行預(yù)估。

隨著CFD 技術(shù)的快速發(fā)展和計(jì)算機(jī)性能的不斷提高，采用商用軟件來(lái)預(yù)報(bào)小展弦比翼水動(dòng)力性能成為翼水動(dòng)力預(yù)估的重要手段。文獻(xiàn)［3-4］表明采用定常計(jì)算方式結(jié)合N -S 方程(RANS 方法)及現(xiàn)有湍流模型能夠較為準(zhǔn)確的得到小攻角下翼的升阻性能，然而對(duì)于大攻角下翼尾端出現(xiàn)非定常分離流動(dòng)后采用RANS 方法還不能合理預(yù)報(bào)翼的水動(dòng)力性能，黃少峰［5］的研究結(jié)果指出，采用考慮和不考慮物面剪切應(yīng)力的方法均不能準(zhǔn)確預(yù)報(bào)大攻角翼升力系數(shù)，而對(duì)這2 種方法的計(jì)算結(jié)果取平均可作為翼升力系數(shù)的參考值;李鋒［6］等則通過(guò)在N-S 方程中添加自定義粘性項(xiàng)得到了適用于預(yù)報(bào)特定Mach 數(shù)翼型大攻角周期性繞流的方法。在不考慮計(jì)算量的情況下，能否通過(guò)采用非定常計(jì)算方式結(jié)合N -S 方程(URANS 方法)來(lái)精確捕捉翼尾端的分離流動(dòng)及如何通過(guò)細(xì)化網(wǎng)格來(lái)提高宏觀力預(yù)報(bào)精度還有待進(jìn)一步研究。

另一方面，近年來(lái)計(jì)算機(jī)性能的快速提高使得應(yīng)用LES 和DES 等方法預(yù)報(bào)翼水動(dòng)力性能成為可能，陳江濤［7］等應(yīng)用基于S -A 湍流模型的DES 方法模擬了3.5 ×106＜Re ＜8.4 ×106時(shí)的圓柱繞流，并較為準(zhǔn)確的得到了流場(chǎng)信息;Breuer M［8］等應(yīng)用RANS、DES 與LES 方法模擬了高雷諾數(shù)下大攻角(18°)平板繞流，結(jié)果表明DES 方法與LES 方法得到的流場(chǎng)信息更為相似;Huang B［9］等應(yīng)用DES 方法較為準(zhǔn)確模擬了水翼表面的非定?？栈鲃?dòng)。然而對(duì)于水翼尾端的分離流動(dòng)能否通過(guò)DES 方法來(lái)準(zhǔn)確模擬及方法的適用條件還有待進(jìn)一步探討。

本文應(yīng)用Fluent 商用軟件以NACA0018 小展弦比翼為研究對(duì)象，重點(diǎn)探討了不同雷諾數(shù)下計(jì)算模型適用性及不同模型下網(wǎng)格離散對(duì)大攻角翼水動(dòng)力計(jì)算結(jié)果的影響規(guī)律，并對(duì)翼失速后尾端產(chǎn)生的非定常分離渦進(jìn)行模擬和探究。

1 控制方程及適用條件

本文以納維-斯托克斯方程(N -S 方程)作為求解不可壓縮牛頓流體運(yùn)動(dòng)的基礎(chǔ)方程。RANS方程是對(duì)N - S 方程進(jìn)行時(shí)均化得到，RANS 與URANS 方法分別基于RANS 方程采用定常與非定常的計(jì)算方式進(jìn)行數(shù)值求解。LES 模型是從三維非定常N-S 方程出發(fā)，通過(guò)濾波計(jì)算將小于某一格子尺度的小渦諸量過(guò)濾掉從而只針對(duì)大渦運(yùn)動(dòng)進(jìn)行模擬的數(shù)值計(jì)算方法，其中RANS 和LES 方程可以表達(dá)為同一種形式:

式中:p 為靜壓;μ 為流體粘度;ui和uj為速度分量。Re 為特征長(zhǎng)度雷諾數(shù);表征由分子運(yùn)動(dòng)引起的動(dòng)量傳輸，其具體表達(dá)見(jiàn)式(3)。在RANS 與LES 方程中(·)項(xiàng)分別代表大尺度渦的計(jì)算項(xiàng)和時(shí)均項(xiàng)，為了更好地模擬由湍流運(yùn)動(dòng)引起的動(dòng)量傳輸，在LES 與RANS 的動(dòng)量方程(2)中增加非線性對(duì)流項(xiàng)

在應(yīng)用LES 模型時(shí)對(duì)計(jì)算空間網(wǎng)格的離散有較高的要求，具體如下:1)近壁面第1 層網(wǎng)格處y +值小于1;2)各方向上的網(wǎng)格尺度擴(kuò)展比率在1.15左右;3)各方向上的網(wǎng)格尺度相差不大。

如果在整個(gè)計(jì)算域中均按此原則進(jìn)行網(wǎng)格劃分必然使得網(wǎng)格數(shù)量急劇增多，為減小網(wǎng)格數(shù)量Spalart［10］提出了DES 方法，該方法結(jié)合了RANS 和LES 方法的優(yōu)點(diǎn)，應(yīng)用DES 方法進(jìn)行計(jì)算時(shí)，物體近壁面采用基于RANS 方程的湍流模型進(jìn)行求解而針對(duì)遠(yuǎn)離物面的的分離區(qū)域則采用LES 模型進(jìn)行求解，其在分離渦計(jì)算方面的精度要高于RANS、URANS 方法而所需網(wǎng)格數(shù)量要遠(yuǎn)小于LES 方法。

本文計(jì)算中涉及RNG k －ε、SST k －ω 和SA 三種湍流模型，詳細(xì)推導(dǎo)過(guò)程和參數(shù)的取值可參考文獻(xiàn)［11］。

2 計(jì)算對(duì)象與網(wǎng)格劃分

2.1 計(jì)算對(duì)象

文中以NACA0018 小展弦比翼作為計(jì)算對(duì)象，翼平面形狀為矩形，模型主尺度見(jiàn)表1。

表1 計(jì)算模型主尺度Tab.1 Main dimensions of computational model

2.2 邊界條件及網(wǎng)格劃分

由于計(jì)算對(duì)象關(guān)于過(guò)展長(zhǎng)中點(diǎn)剖面對(duì)稱，為減小計(jì)算量，設(shè)置過(guò)展長(zhǎng)中點(diǎn)剖面為對(duì)稱面取1/2 翼進(jìn)行計(jì)算，計(jì)算域以及邊界條件(見(jiàn)圖1)設(shè)置如下:

1)入口:入口距離翼導(dǎo)邊4 倍弦長(zhǎng)，邊界條件設(shè)置為速度入口;

2)出口:出口距離翼隨邊7 倍弦長(zhǎng)，邊界條件設(shè)置為壓力出口，壓力為未擾動(dòng)時(shí)邊界壓力;

3)流域的外邊界距翼的縱向中心線四倍弦長(zhǎng)，邊界條件設(shè)置為速度入口，速度大小方向與入口相同;

4)在翼表面定義無(wú)滑移、不可穿透邊界條件。

圖1 計(jì)算域Fig.1 Computational domain

圖2 翼表面網(wǎng)格劃分Fig.2 Grid partition on hydrofoil

為更好的控制各方向上的網(wǎng)格尺度擴(kuò)展比率和提高網(wǎng)格質(zhì)量，在翼型外圍采用多塊分梯度加密的方法并在翼頂部與外圍進(jìn)行O 型拓?fù)涮幚?，翼表面網(wǎng)格劃分如圖2 所示。

3 數(shù)值計(jì)算分析

3.1 計(jì)算模型適用性

在應(yīng)用CFD 進(jìn)行預(yù)報(bào)時(shí)，計(jì)算條件尤其是雷諾數(shù)的大小是計(jì)算模型選取的重要依據(jù)，試驗(yàn)資料［12］表明當(dāng)雷諾數(shù)大于1.2 ×105時(shí)僅失速角隨雷諾數(shù)的增大而增大，在失速角范圍內(nèi)一定攻角下雷諾數(shù)的增大不會(huì)引起翼升阻系數(shù)的變化，由于試驗(yàn)測(cè)得的升力系數(shù)可靠性較高(誤差3%以內(nèi))，本節(jié)以Re=2.7 ×105時(shí)的升力試驗(yàn)結(jié)果作為不同雷諾數(shù)下宏觀力預(yù)報(bào)準(zhǔn)確性的評(píng)判依據(jù)來(lái)探討不同雷諾數(shù)下翼未失速時(shí)計(jì)算模型的適用性。

3.1.1 Re ＜106時(shí)計(jì)算模型適用性

為考察Re ＜106時(shí)計(jì)算模型適用性并減小由網(wǎng)格離散與計(jì)算方式帶來(lái)的影響，以極細(xì)網(wǎng)格(滿足DES 方法要求，網(wǎng)格數(shù)320 萬(wàn)，y + ～1)為基礎(chǔ)，采用非定常計(jì)算方式(時(shí)間步長(zhǎng)取0.0001s)，開展了Re=2.7 ×105、不同模型下10°、15°及18°攻角翼升阻系數(shù)的數(shù)值預(yù)報(bào)并與試驗(yàn)值進(jìn)行對(duì)比(見(jiàn)表2)。計(jì)算結(jié)果表明:當(dāng)Re ＜106時(shí)，采用URANS 方法的預(yù)報(bào)精度較高，其中以RNG k－ε 湍流模型的預(yù)報(bào)結(jié)果與試驗(yàn)值吻合最好，DES SST 與DES SA 模型計(jì)算所得升力系數(shù)較試驗(yàn)值偏差較大。

翼尾端流動(dòng)分離位置和分離范圍的準(zhǔn)確模擬是準(zhǔn)確預(yù)報(bào)宏觀力的關(guān)鍵，由不同湍流模型計(jì)算結(jié)果與所得的流場(chǎng)信息(見(jiàn)圖3)可知，在Re ＜106時(shí)，采用DES 方法得到的流動(dòng)分離位置明顯靠前，其模擬得到的流動(dòng)分離現(xiàn)象與實(shí)際流動(dòng)偏差較大。

表2 Re=2.7 ×105 時(shí)不同計(jì)算模型翼升阻系數(shù)計(jì)算結(jié)果Tab.2 Lift and drag coefficients calculated by different models at Re=2.7 ×105

圖3 18°攻角下不同模型計(jì)算所得展長(zhǎng)中點(diǎn)剖面處流場(chǎng)信息Fig.3 Flow fields of mid-elongation section calculated by different models at 18° incidence

3.1.2 高雷諾數(shù)下計(jì)算模型適用性

為探討高雷諾數(shù)下計(jì)算模型適用性，文中對(duì)比了1.35 ×106＜Re ＜2.7 ×107時(shí)翼升阻系數(shù)的計(jì)算結(jié)果 (見(jiàn)表3)，由表3 可知，在高雷諾數(shù)下URANS 方法的預(yù)報(bào)精度均較低，而當(dāng)Re ＞2.7 ×106時(shí)，DES 方法的計(jì)算結(jié)果與試驗(yàn)值吻合較好。

表3 高雷諾數(shù)時(shí)不同計(jì)算模型下翼升阻系數(shù)對(duì)比Tab.3 Comparison of lift and drag coefficients calculated by different models at high Re

3.2 網(wǎng)格離散的影響

在應(yīng)用Fluent 軟件進(jìn)行水動(dòng)力預(yù)報(bào)時(shí)，不同湍流模型采用的壁面處理方式存在差異而且對(duì)第1 層網(wǎng)格y+值的要求有一定的參考范圍［13］。由表3 可知，在采用加強(qiáng)型壁面函數(shù)進(jìn)行壁面處理時(shí)5 種湍流模型對(duì)近壁面y+要求較高，在網(wǎng)格不夠細(xì)(y +＞30)的情況下可以應(yīng)用RNG k － ε 湍流模型并在近壁面處施加壁面函數(shù)進(jìn)行處理，在處理復(fù)雜幾何流場(chǎng)計(jì)算問(wèn)題時(shí)此模型更適用。

表4 不同湍流模型壁面處理方式及壁面y+參考值Tab.4 Reference value of y+ based on different wall treatment methods

為了研究壁面y +值及網(wǎng)格離散密度對(duì)翼水動(dòng)力計(jì)算結(jié)果的影響，文中對(duì)比了18°攻角下5 套不同網(wǎng)格宏觀力系數(shù)的計(jì)算結(jié)果(見(jiàn)表5)，為滿足DES 模型的網(wǎng)格離散要求將網(wǎng)格進(jìn)行細(xì)化(fine-1、fine-2)，在應(yīng)用URANS 和DES 方法時(shí)分別取計(jì)算雷諾數(shù)Re=2.7 ×105和2.7 ×107。由表5 可看出，在同一模型下壁面y +值的選取對(duì)計(jì)算結(jié)果的影響較大，當(dāng)y + 在30 ～110 范圍內(nèi)時(shí)，采用URANS RNG k －ε 模型結(jié)合壁面函數(shù)的方法較其他模型所得結(jié)果相對(duì)較為合理但誤差仍較大，當(dāng)y + 在3 ～30范圍時(shí)，除URANS RNG k － ε 模型計(jì)算結(jié)果相對(duì)較為可信外，其余模型誤差均較大，因此在計(jì)算大攻角翼水動(dòng)力性能時(shí)，有必要將壁面y +增強(qiáng)至1 左右。通過(guò)對(duì)比相同y +下粗網(wǎng)格與細(xì)網(wǎng)格的計(jì)算結(jié)果可知在滿足y+ ～1 的條件下，增加網(wǎng)格密度可提高計(jì)算精度，如果模型較為復(fù)雜(如船后半懸掛舵、襟翼舵)時(shí)網(wǎng)格離散很難達(dá)到y(tǒng)+ ～1 的條件此時(shí)建議采用URANS RNG k － ε 模型并對(duì)計(jì)算結(jié)果進(jìn)行適當(dāng)修正。

表5 18°攻角不同網(wǎng)格方案下翼升阻系數(shù)計(jì)算結(jié)果Tab.5 Lift and drag coefficients based on different grid partitions at 18° incidence

3.3 翼失速后的流動(dòng)分離模擬

前面應(yīng)用非定常計(jì)算方式討論了湍流模型及網(wǎng)格離散對(duì)計(jì)算結(jié)果的影響，然而在工程應(yīng)用中采用定常計(jì)算方式(RANS 方法)可大大減少工程周期，本節(jié)探討了定常計(jì)算方式的適用性及其預(yù)報(bào)不同攻角下翼宏觀力的可行性并應(yīng)用不同方法對(duì)翼尾端分離流動(dòng)進(jìn)行模擬對(duì)比。經(jīng)計(jì)算可知，當(dāng)計(jì)算攻角小于失速角時(shí)翼尾端分離范圍較小，采用基于RNG k －ε、SST k －ω與SA 湍流模型的RANS 與URANS 方法所得計(jì)算結(jié)果相差不大且流場(chǎng)信息也基本相同，然而在實(shí)際工程應(yīng)用中RANS 方法預(yù)報(bào)的失速角存在很大偏差，有必要對(duì)翼失速前后采用RANS 與URANS 方法預(yù)報(bào)結(jié)果的偏差進(jìn)行討論。如表所示，當(dāng)計(jì)算攻角大于失速角(失速角為21°)時(shí)，采用RANS 與URANS 方法所得的升阻系數(shù)差別較大。由圖4 也可看出，翼失速后采用RANS 與URANS 方法捕捉得到的分離區(qū)大小和渦核位置存在較大差異，因此在大攻角下不宜采用RANS 方法進(jìn)行預(yù)報(bào)。從計(jì)算結(jié)果的對(duì)比可知，在未知失速角情況下可以通過(guò)比較RANS 與URANS 方法的計(jì)算結(jié)果大致估計(jì)攻角是否小于失速角。

表6 翼失速后(24°攻角)不同計(jì)算模型所得升阻系數(shù)Tab.6 Lift and drag coefficients based on different turbulence models when hydrofoil stalls

圖4 翼失速后(24°攻角)不同計(jì)算模型所得展長(zhǎng)中剖面流線圖Fig.4 Flow fields of mid-elongation section calculated by different models when hydrofoil stalls (24° incidence)

雖然翼失速后采用URANS 方法模擬得到的分離渦與RANS 有較大不同，但URANS 方法所得渦核位置及分離區(qū)范圍均不隨時(shí)間發(fā)生變化，翼的升阻系數(shù)不具有時(shí)歷波動(dòng)性。為探究DES 方法對(duì)于分離渦的捕捉能力，結(jié)合DES 方法的雷諾數(shù)適用范圍，取Re=2.7×107，攻角35°作為翼迎流條件開展計(jì)算，結(jié)果表明DES SST 與DES SA 模型計(jì)算所得宏觀力均表現(xiàn)時(shí)歷波動(dòng)性，觀察到的分離渦隨時(shí)間的生成與脫落(見(jiàn)圖5)與實(shí)際物理現(xiàn)象一致，但其模擬的準(zhǔn)確度還需試驗(yàn)驗(yàn)證，DES 模型的優(yōu)化有待進(jìn)一步研究。

圖5 35°攻角下DES SST 模型計(jì)算所得展長(zhǎng)中剖面流線圖Fig.5 Flow fields of mid-elongation section calculated by DES SST model at 35° incidence

4 結(jié) 語(yǔ)

通過(guò)對(duì)比分析不同計(jì)算條件下基于不同湍流模型的RANS、URANS 與DES 方法在大攻角翼型水動(dòng)力預(yù)報(bào)中的宏觀力計(jì)算結(jié)果與流場(chǎng)信息可以得出以下結(jié)論:

1)計(jì)算方法和湍流模型的選取對(duì)翼水動(dòng)力預(yù)報(bào)精度影響最大且與雷諾數(shù)相關(guān);當(dāng)Re ＜106時(shí)，URANS 方法的預(yù)報(bào)精度較高，其中以RNG k － ε 湍流模型的預(yù)報(bào)結(jié)果吻合最好，而DES 方法的精度較低，當(dāng)Re ＞107時(shí)，采用DES 方法能夠得到較為可靠的預(yù)報(bào)結(jié)果;

2)壁面y +對(duì)翼型水動(dòng)力預(yù)報(bào)精度有重要影響，確保y+ ～1 的同時(shí)適當(dāng)加密網(wǎng)格數(shù)量可提高計(jì)算精度，在無(wú)法滿足壁面y+要求時(shí)可采用URANS RNG k －ε 模型進(jìn)行預(yù)報(bào)并對(duì)計(jì)算結(jié)果進(jìn)行適當(dāng)修正;

3)翼型失速后，RANS 與URANS 方法計(jì)算所得升阻系數(shù)差別較大，這2 種方法均不能捕捉分離渦的時(shí)歷特性，DES 方法能夠觀察到與物理現(xiàn)象一致的漩渦結(jié)構(gòu)。905 -909.

［1］梅琴生.船用舵［M］.北京:人民交通出版社，1981.

［2］陳材侃，劉華.解三維水翼繞流的下潛渦環(huán)柵格法［J］.船舶力學(xué)，2005，9(2):41 -45.CHEN Cai-kan，LIU Hua.A submerged vortex lattice method for calculation of the flow around a three-dimension hydrofoil［J］.Journal of Ship Mechanics，2005，9(2):41-45.

［3］PARK K，LEE J. Influence of endplate on aerodynamic characteristics of low-aspect-ratio wing in ground effect［J］.Journal of Mechanical Science，2008，22(12):2578-2589.

［4］GIM O S，LEE G W. Flow characteristics and tip vortex formation around a NACA 0018 foil with an endplate［J］.Ocean Engineering，2013，60:28 -38.

［5］黃少鋒，劉登成，張志榮.基于精細(xì)流場(chǎng)分析的大攻角翼型的數(shù)值計(jì)算與試驗(yàn)研究［C］//2013 年船舶水動(dòng)力學(xué)學(xué)術(shù)會(huì)議，2013:365 -370.

［6］李鋒，汪翼云，崔爾杰.翼型大攻角繞流的數(shù)值模擬［J］.航空學(xué)報(bào)，1992.LI Feng，WANG Yi-yun，CUI Er-jie.The numerical simulation of compressible flow around an airfoil at high angle of attack［J］.Acta Aeronautica et Astronautica Sinica，1992.

［7］陳江濤，張培紅，周乃春，等. 基于SA 湍流模型的DES方法應(yīng)用［J］. 北京航空航天大學(xué)學(xué)報(bào)，2012，38(7):CHEN Jiang-tao，ZHANG Pei-hong，ZHOU Nai-chun，et al.Application of detached-eddy simulation based on spalartallmaras turbulence model［J］.Journal of Beijing University of Aeronautics and Astronautics，2012，38(7):905 -909.

［8］BREUER M，JOVIN，MAZAEV K. Comparison of DES，RANS and LES for the separated flow around a flat plate at high incidence［J］. International Journal for Numerical Methods in Fluids，2003，41(4):357 -388.

［9］HUANG B，WANG G，YU Z，et al. Detached-eddy simulation for time-dependent turbulent cavitating flows［J］.Chinese Journal of Mechanical Engineering，2012，25(3):484 -490.

［10］SPALART P R，JOU W H，STRELETS M，et al. Comments on the feasibility of LES for wings，and on a hybrid RANS/LES approach［J］.Advances in DNS/LES，1997，1:4 -8.

［11］張兆順，崔桂香，許春曉.湍流理論與模擬［M］.北京:清華大學(xué)出版社有限公司，2005.

［12］陸惠生，朱文蔚，費(fèi)乃振，等.敞水舵的試驗(yàn)研究［J］.上海交通大學(xué)學(xué)報(bào)，1981(2):001.LU Hui-sheng，ZHU Wen-wei，F(xiàn)EI Nai-zhen，et al.Experiment study on open rudders［J］.Journal of Shanghai Jiaotong University，1981(2):001.

［13］ANSYS A F.14.0 Theory Guide［J］.ANSYS inc，2011.