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

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

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

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

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

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

1.1 幾何模型

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

圖1 幾何模型Fig.1 Geometry model

1.2 數(shù)值方法

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

圖2 翼型周?chē)W(wǎng)格分布Fig.2 Mesh distribution around the airfoil

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

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

表1 計(jì)算方案Tab.1 Calculation plan

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

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

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

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

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

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

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

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

表2 網(wǎng)格對(duì)預(yù)測(cè)轉(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方法計(jì)算結(jié)果的對(duì)比

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

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

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

圖8為DU00-W2-401翼型中間對(duì)稱(chēng)面上渦量云圖.從圖8 可以看出，DES方法可以模擬出更細(xì)致的漩渦結(jié)構(gòu).從DES計(jì)算結(jié)果可以看出在分離區(qū)的上游有組織、有規(guī)律地脫落出分離渦，而在翼型的尾緣處，翼型上表面分離區(qū)內(nèi)是低壓區(qū)，翼型下表面是高壓區(qū)，翼型下表面的流體會(huì)繞過(guò)葉片尾緣卷入分離區(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可以看出，隨著攻角的增大，翼型上表面的湍流分離渦尺度越來(lái)越大，分離區(qū)尺度越來(lái)越大.由于大厚度翼型極易發(fā)生流動(dòng)分離，分離之后RANS方法不能模擬出細(xì)致的漩渦結(jié)構(gòu)，導(dǎo)致其預(yù)測(cè)氣動(dòng)性能不準(zhǔn)確，而DES方法能較好地模擬出細(xì)致的湍流脈動(dòng)漩渦，所以DES 方法較RANS 方法更適合大厚度翼型的氣動(dòng)計(jì)算.

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

3 結(jié) 論

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

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

［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方法模擬不同類(lèi)型翼型的失速特性［J］.航空學(xué)報(bào)，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ù)對(duì)風(fēng)力機(jī)翼型轉(zhuǎn)捩特性的影響［J］.動(dòng)力工程學(xué)報(bào)，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)力機(jī)定常與非定常氣動(dòng)問(wèn)題的數(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.