周廣禮，姚朝幫

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

0 引 言

小展弦比三維舵翼的水動力預報是船舶操縱性預報的前提和基礎，在以往進行舵初步設計時主要采用相關翼型已有的試驗資料進行展弦比換算［1］或基于勢流理論的方法來得到翼的水動力性能［2］，然而采用展弦比換算方法所得計算結果誤差較大，基于勢流理論的方法并不能準確預估大攻角翼的水動力性能而且無法對幾何形狀較為復雜的舵(如半懸掛舵、襟翼舵)的水動力進行預估。

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

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

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

1 控制方程及適用條件

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

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

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

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

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

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

2.1 計算對象

文中以NACA0018 小展弦比翼作為計算對象，翼平面形狀為矩形，模型主尺度見表1。

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

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

由于計算對象關于過展長中點剖面對稱，為減小計算量，設置過展長中點剖面為對稱面取1/2 翼進行計算，計算域以及邊界條件(見圖1)設置如下:

1)入口:入口距離翼導邊4 倍弦長，邊界條件設置為速度入口;

2)出口:出口距離翼隨邊7 倍弦長，邊界條件設置為壓力出口，壓力為未擾動時邊界壓力;

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

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

圖1 計算域Fig.1 Computational domain

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

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

3 數(shù)值計算分析

3.1 計算模型適用性

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

3.1.1 Re ＜106時計算模型適用性

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

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

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

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

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

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

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

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

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

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

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

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

3.3 翼失速后的流動分離模擬

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

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

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

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

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

4 結 語

通過對比分析不同計算條件下基于不同湍流模型的RANS、URANS 與DES 方法在大攻角翼型水動力預報中的宏觀力計算結果與流場信息可以得出以下結論:

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

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

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

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

［2］陳材侃，劉華.解三維水翼繞流的下潛渦環(huán)柵格法［J］.船舶力學，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］黃少鋒，劉登成，張志榮.基于精細流場分析的大攻角翼型的數(shù)值計算與試驗研究［C］//2013 年船舶水動力學學術會議，2013:365 -370.

［6］李鋒，汪翼云，崔爾杰.翼型大攻角繞流的數(shù)值模擬［J］.航空學報，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方法應用［J］. 北京航空航天大學學報，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］.北京:清華大學出版社有限公司，2005.

［12］陸惠生，朱文蔚，費乃振，等.敞水舵的試驗研究［J］.上海交通大學學報，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.