風(fēng)力機大厚度鈍尾緣翼型數(shù)值模擬研究

李新凱，戴麗萍，康 順，2，梁思超

（1.華北電力大學(xué)電站設(shè)備狀態(tài)監(jiān)測與控制教育部重點實驗室，北京102206；2.西安現(xiàn)代控制技術(shù)研究所，西安710065）

隨著風(fēng)力機單機功率的不斷提高，風(fēng)力機葉片越來越長，為保證強度要求，葉片根部往往用大厚度鈍尾緣翼型［1-2］.葉片運行過程中根部會發(fā)生較大的流動分離，降低風(fēng)輪對風(fēng)能的捕獲功率.所以研究大厚度翼型的氣動特性，為改善葉片根部流動狀態(tài)提供理論依據(jù)非常有必要.進行大厚度翼型風(fēng)洞實驗較為困難，目前高雷諾數(shù)風(fēng)洞實驗數(shù)據(jù)較少［3-7］，而傳統(tǒng)的雷諾時均（RANS）方法對產(chǎn)生較大分離之后翼型氣動性能的預(yù)測效果較差，所以需要一種可以較準(zhǔn)確地預(yù)測大厚度翼型氣動性能的數(shù)值計算方法.

Li等［8］利用大渦模擬（LES）方法對2.5D的NACA0018翼型在全攻角下進行了數(shù)值計算，結(jié)果表明在大部分攻角下LES方法均能較好地預(yù)測翼型的升力系數(shù)與阻力系數(shù)，效果明顯優(yōu)于非定常雷諾方法（URANS）方法.李棟等［9］采用分離渦模擬（DES）計算方法對3種小厚度的NACA 翼型進行了數(shù)值計算，通過與實驗結(jié)果比較，顯示該方法可以有效地預(yù)測翼型的失速特征.Kravchenko等［10］利用LES方法分析了不同網(wǎng)格密度及展向?qū)挾鹊膱A柱繞流，通過計算發(fā)現(xiàn)圓柱展向?qū)挾燃罢瓜蚓W(wǎng)格分辨率對計算結(jié)果有很大影響.高偉等［11］采用勢流方程與邊界層耦合方法對不同厚度翼型邊界層轉(zhuǎn)捩進行了研究，結(jié)果表明轉(zhuǎn)捩位置對翼型升阻力系數(shù)有一定影響.

目前，國內(nèi)外學(xué)者針對風(fēng)力機專用大厚度鈍尾緣翼型進行的數(shù)值模擬研究還很少，筆者以具有實驗數(shù)據(jù)的DU97-W-300和DU00-W2-401翼型為研究對象，對其進行了數(shù)值模擬研究.主要研究內(nèi)容包括以下2個方面：（1）利用RANS方法對翼型進行計算，討論不同湍流模型（包括是否考慮轉(zhuǎn)捩）以及網(wǎng)格分布對RANS方法計算結(jié)果的影響；（2）分別采用RANS 和DES 方法對DU00-W2-401翼型在多個攻角下進行數(shù)值計算，對比RANS 與DES 方法的計算結(jié)果.

1 幾何模型及數(shù)值方法

1.1 幾何模型

圖1為DU97-W-300與DU00-W2-401翼型的幾何模型.其中，DU97-W-300 翼型弦長C＝0.6 mm，相對厚度為30%C，鈍尾緣，尾緣厚度為1.74%C；DU00-W2-401翼型弦長C＝0.6 mm，相對厚度為40%C，鈍尾緣，尾緣厚度為1%C.

圖1 幾何模型Fig.1 Geometry model

1.2 數(shù)值方法

用ICEM 軟件進行網(wǎng)格劃分，生成全域結(jié)構(gòu)網(wǎng)格，第一次網(wǎng)格高度0.001mm，網(wǎng)格數(shù)目分布見表1，翼型網(wǎng)格分布見圖2.

圖2 翼型周圍網(wǎng)格分布Fig.2 Mesh distribution around the airfoil

邊界條件：所有計算均采用圓形計算域，計算域半徑R＝30C，采用速度進口、壓力出口邊界條件，葉片表面設(shè)置為無滑移固體壁面.

采用商用軟件Fluent進行數(shù)值計算，定常計算采用RANS 方法，非定常計算采用DES 方法.RANS方法中湍流模型包括SA、SST 以及考慮轉(zhuǎn)捩的SST（T-SST）湍流模型，DES方法中湍流模型為SA湍流模型，具體方案設(shè)置見表1.采用有限體積方法對控制方程進行離散，壓力-速度耦合基于Simple算法.RANS方法中控制方程的各項均采用二階迎風(fēng)格式，DES方法中壓力項采用二階迎風(fēng)格式，動量項采用中心差分格式.

表1 計算方案Tab.1 Calculation plan

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

2.1 湍流模型及網(wǎng)格對RANS方法計算結(jié)果的影響

圖3 給出了不同湍流模型下DU97-W-300 翼型的計算結(jié)果.圖中，Cl為升力系數(shù)，Cd為阻力系數(shù)α為攻角.從圖3 可以看出，在翼型失速之前，TSST 湍流模型更能準(zhǔn)確地預(yù)測翼型的氣動性能，升力系數(shù)、阻力系數(shù)及升阻比更接近實驗值.但在翼型失速之后，全湍流模型（SA 和SST 模型）比T-SST湍流模型計算結(jié)果更接近實驗值.對比可以發(fā)現(xiàn)，無論是失速之前還是失速之后，相比全湍流模型，T-SST湍流模型計算得到的升力系數(shù)偏大，阻力系數(shù)偏小，升阻比偏大.這是由于考慮轉(zhuǎn)捩之后，葉片前緣部分流動變?yōu)閷恿?，從而使得葉片升力系數(shù)偏大，阻力系數(shù)偏小，這與文獻［12］和文獻［13］中得出的結(jié)論一致.全湍流模型中SA 湍流模型計算效果要優(yōu)于SST 湍流模型，計算值更加接近實驗值.

圖3 DU97-W-300翼型計算結(jié)果Fig.3 Computation result of DU97-W-300airfoil

圖4（a）為DU97-W-300翼型表面壓力系數(shù)Cp分布曲線.從圖4（a）可以看出，在9.3°、12.4°攻角下，T-SST 湍流模型的計算結(jié)果與實驗值更加吻合，SST 湍流模型計算效果最差.轉(zhuǎn)捩是翼型繞流中重要的流動現(xiàn)象，轉(zhuǎn)捩流動現(xiàn)象復(fù)雜，一直是CFD 模擬中的難點，翼型轉(zhuǎn)捩位置的預(yù)測從一定程度上反映了CFD 的計算精度.隨著對翼型繞流流動細節(jié)的研究，邊界層轉(zhuǎn)捩成為重要的研究對象.圖4（b）為CFD 預(yù)測翼型轉(zhuǎn)捩位置與實驗值的對比曲線.從圖4（b）可以看出，隨著攻角的增大，翼型吸力面轉(zhuǎn)捩位置向翼型前緣移動，壓力面轉(zhuǎn)捩位置向翼型尾緣移動.從計算轉(zhuǎn)捩位置來看，T-SST 湍流模型計算得到的翼型壓力面和吸力面轉(zhuǎn)捩位置與實驗值吻合較好，說明考慮轉(zhuǎn)捩的T-SST 湍流模型能較好地預(yù)測翼型轉(zhuǎn)捩位置.

圖4 DU97-W-300翼型表面壓力系數(shù)及轉(zhuǎn)捩位置分布Fig.4 Distribution of Cpand transition position for DU97-W-300airfoil

由上面的分析可以看出，在翼型線性段，T-SST湍流模型能更好地預(yù)測葉片升力系數(shù)、阻力系數(shù)及轉(zhuǎn)捩位置.筆者通過改變翼型一周網(wǎng)格數(shù)（方案A）及第一層網(wǎng)格高度ΔY（方案B）來研究網(wǎng)格對預(yù)測轉(zhuǎn)捩位置準(zhǔn)確性的影響.具體方案及計算結(jié)果見表2，計算攻角選為α＝12.4°.

圖5為不同網(wǎng)格數(shù)及網(wǎng)格高度下DU97-W-300翼型表面壓力系數(shù)及轉(zhuǎn)捩位置分布圖.從圖5可以看出，方案A 中1、2、3、4計算所得Cp均與實驗值吻合較好，且均能正確預(yù)測轉(zhuǎn)捩位置，所以翼型一周網(wǎng)格數(shù)對預(yù)測轉(zhuǎn)捩位置的影響較小.方案B 中5、6能較準(zhǔn)確地預(yù)測出轉(zhuǎn)捩位置，且Cp與實驗值吻合較好，但方案B 中7、8未能計算出翼型表面的轉(zhuǎn)捩位置，且Cp與實驗值相差較大，由此可見對于預(yù)測轉(zhuǎn)捩位置，第一層網(wǎng)格高度比較重要，當(dāng)ΔY／C≥1.66×10-4（即y＋≥18）時已經(jīng)不能預(yù)測出轉(zhuǎn)捩位置.

表2 網(wǎng)格對預(yù)測轉(zhuǎn)捩效果的影響Tab.2 Influence of mesh distribution on the predicted transition effects

圖5 不同網(wǎng)格數(shù)及網(wǎng)格高度下DU97-W-300翼型表面壓力系數(shù)及轉(zhuǎn)捩位置分布Fig.5 Distribution of surface pressure coefficient and transition position for DU97-W-300airfoil at different grid number and height

2.2 DES與RANS方法計算結(jié)果的對比

圖6給出了DU00-W2-401翼型DES與RANS方法計算結(jié)果的對比.其中計算攻角為5°～23°.RANS方法中湍流模型選擇SA 湍流模型，DES方法中近壁區(qū)域湍流模型也選擇SA 湍流模型.從圖6可以看出，相比于RANS方法的計算結(jié)果，DES方法能更好地預(yù)測大厚度翼型的氣動特性.相比實驗值，在翼型失速段，RANS 方法計算所得升力系數(shù)偏大，阻力系數(shù)偏小，升阻比偏大.與DES方法計算結(jié)果相比，在翼型失速之后，RANS 方法計算所得升力系數(shù)偏大，阻力系數(shù)偏小，升阻比偏大.DES方法計算得到的翼型升力系數(shù)、阻力系數(shù)及升阻比與實驗值吻合更好.從計算結(jié)果來看，在翼型失速之后，DES方法的模擬結(jié)果明顯優(yōu)于RANS方法.

圖6 DU00-W2-401翼型計算結(jié)果Fig.6 Computation result of DU00-W2-400airfoil

圖7 為DU00-W2-401 翼型表面壓力系數(shù)分布.從圖7可以看出，在15°、19°攻角時，RANS計算結(jié)果的吸力峰值均比DES方法的計算結(jié)果要高，所以計算所得氣動力較DES方法的計算結(jié)果偏大.

圖8為DU00-W2-401翼型中間對稱面上渦量云圖.從圖8 可以看出，DES方法可以模擬出更細致的漩渦結(jié)構(gòu).從DES計算結(jié)果可以看出在分離區(qū)的上游有組織、有規(guī)律地脫落出分離渦，而在翼型的尾緣處，翼型上表面分離區(qū)內(nèi)是低壓區(qū)，翼型下表面是高壓區(qū)，翼型下表面的流體會繞過葉片尾緣卷入分離區(qū)，這樣分離渦周而復(fù)始向下游脫落.

圖7 DU00-W2-401翼型表面壓力系數(shù)分布Fig.7 Distribution of surface pressure coefficient

圖8 渦量云圖Fig.8 Cloud map of vorticity

圖9 為帶速度云圖的Q等值面（Q＝1 000 s-2）.從圖9可以看出，隨著攻角的增大，翼型上表面的湍流分離渦尺度越來越大，分離區(qū)尺度越來越大.由于大厚度翼型極易發(fā)生流動分離，分離之后RANS方法不能模擬出細致的漩渦結(jié)構(gòu)，導(dǎo)致其預(yù)測氣動性能不準(zhǔn)確，而DES方法能較好地模擬出細致的湍流脈動漩渦，所以DES 方法較RANS 方法更適合大厚度翼型的氣動計算.

圖9 帶速度云圖的Q 等值面（Q＝1 000s-2）Fig.9 Isosurface of Q with velocity contour

3 結(jié) 論

（1）對于RANS方法的計算，在翼型最大升力系數(shù)攻角之前，考慮轉(zhuǎn)捩的T-SST 湍流模型的計算效果要優(yōu)于全湍流模型，在翼型最大升力系數(shù)攻角之后，全湍流模型的計算效果要優(yōu)于T-SST 湍流模型，全湍流模型中SA 湍流模型的計算效果優(yōu)于SST 湍流模型.網(wǎng)格對計算結(jié)果及轉(zhuǎn)捩預(yù)測的影響中，第一層網(wǎng)格高度更加敏感，而翼型一周網(wǎng)格數(shù)的影響不大.

（2）對比RANS與DES方法的計算結(jié)果發(fā)現(xiàn)，當(dāng)葉片存在大尺度分離時，RANS方法不能模擬出更細致的分離渦結(jié)構(gòu)，而DES方法則能在一定程度上模擬出細致規(guī)律的分離渦結(jié)構(gòu).對比氣動性能計算結(jié)果，DES 方法計算結(jié)果與實驗值更加吻合.DES方法比RANS方法更適合大厚度翼型的氣動性能計算.

［1］TPI Composites.Innovative design approaches for large wind turbine blades-final report［R］.Sandia，USA：TPI，2004.

［2］STANDISH K J，van DAM C P.Aerodynamic analysis of blunt trailing edge airfoils［J］.Journal of Solar Energy Engineering，2003，125（4）：479-487.

［3］BAKER J，MAYDA E，DAM C.Experimental and computational analysis of thick fatback wind turbine airfoils［R］.Reno，USA：AIAA，2006.

［4］LAW S P，GREGNREK G M.Wind turbine evaluation of a truncated NACA 64-621airfoil for wind turbine applications［R］.Ohio，USA：NASA，1987.

［5］TPI Composites.Flatback airfoil wind tunnel experiment［R］.Sandia，USA：TPI，2008.

［6］POST M L，JONES R B.Characterization of a flatback airfoil for use in wind power generation［R］.Reno，USA：AIAA，2008.

［7］BARONE M F，BERG D.Aerodynamic and aeroacoustic properties of a flatback airfoil：an update［R］.Oriando，USA：AIAA，2009.

［8］LI Chao，ZHU Songye，XU Youlin，etal.2.5Dlarge eddy simulation of vertical axis wind turbine in consideration of high angle of attack flow［J］.Renewable Energy，2013，51（3）：317-330.

［9］李棟，焦予秦，IGOR Men＇shov，等.Detached-Eddy Simulation方法模擬不同類型翼型的失速特性［J］.航空學(xué)報，2005，26（4）：406-410.

LI Dong，JIAO Yuqin，IGOR Men＇shov，etal.Detached-Eddy Simulation for airfoil stall［J］.Acta Aeronautica et Astronautica Sinica，2005，26（4）：406-410.

［10］KRAVCHENKO A G，MOIN P.Numerical studies of flow over a circular cylinder atReD＝3 900［J］.Physics of Fluids，2000，12（2）：403-417.

［11］高偉，李春，高月文，等.幾何參數(shù)對風(fēng)力機翼型轉(zhuǎn)捩特性的影響［J］.動力工程學(xué)報，2013，33（6）：490-496.

GAO Wei，LI Chun，GAO Yuewen，etal.Influence of geometric parameters on transition characteristics of wind turbine airfoils［J］.Journal of Chinese Society of Power Engineering，2013，33（6）：490-496.

［12］范忠瑤.風(fēng)力機定常與非定常氣動問題的數(shù)值模擬研究［D］.北京：華北電力大學(xué)，2011.

［13］MICHEL R，ARNAL D，COUSTOLS E.Stability calculations and transition criteria on two-or three-dimensional flows［J］.Laminar-Turbulent Transition，1984（7）：455-461.