风力机大厚度钝尾缘翼型数值模拟研究

李新凯，戴丽萍，康 顺，2，梁思超

（1.华北电力大学电站设备状态监测与控制教育部重点实验室，北京102206；2.西安现代控制技术研究所，西安710065）

随着风力机单机功率的不断提高，风力机叶片越来越长，为保证强度要求，叶片根部往往用大厚度钝尾缘翼型［1-2］.叶片运行过程中根部会发生较大的流动分离，降低风轮对风能的捕获功率.所以研究大厚度翼型的气动特性，为改善叶片根部流动状态提供理论依据非常有必要.进行大厚度翼型风洞实验较为困难，目前高雷诺数风洞实验数据较少［3-7］，而传统的雷诺时均（RANS）方法对产生较大分离之后翼型气动性能的预测效果较差，所以需要一种可以较准确地预测大厚度翼型气动性能的数值计算方法.

Li等［8］利用大涡模拟（LES）方法对2.5D的NACA0018翼型在全攻角下进行了数值计算，结果表明在大部分攻角下LES方法均能较好地预测翼型的升力系数与阻力系数，效果明显优于非定常雷诺方法（URANS）方法.李栋等［9］采用分离涡模拟（DES）计算方法对3种小厚度的NACA 翼型进行了数值计算，通过与实验结果比较，显示该方法可以有效地预测翼型的失速特征.Kravchenko等［10］利用LES方法分析了不同网格密度及展向宽度的圆柱绕流，通过计算发现圆柱展向宽度及展向网格分辨率对计算结果有很大影响.高伟等［11］采用势流方程与边界层耦合方法对不同厚度翼型边界层转捩进行了研究，结果表明转捩位置对翼型升阻力系数有一定影响.

目前，国内外学者针对风力机专用大厚度钝尾缘翼型进行的数值模拟研究还很少，笔者以具有实验数据的DU97-W-300和DU00-W2-401翼型为研究对象，对其进行了数值模拟研究.主要研究内容包括以下2个方面：（1）利用RANS方法对翼型进行计算，讨论不同湍流模型（包括是否考虑转捩）以及网格分布对RANS方法计算结果的影响；（2）分别采用RANS 和DES 方法对DU00-W2-401翼型在多个攻角下进行数值计算，对比RANS 与DES 方法的计算结果.

1 几何模型及数值方法

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 数值方法

用ICEM 软件进行网格划分，生成全域结构网格，第一次网格高度0.001mm，网格数目分布见表1，翼型网格分布见图2.

图2 翼型周围网格分布Fig.2 Mesh distribution around the airfoil

边界条件：所有计算均采用圆形计算域，计算域半径R＝30C，采用速度进口、压力出口边界条件，叶片表面设置为无滑移固体壁面.

采用商用软件Fluent进行数值计算，定常计算采用RANS 方法，非定常计算采用DES 方法.RANS方法中湍流模型包括SA、SST 以及考虑转捩的SST（T-SST）湍流模型，DES方法中湍流模型为SA湍流模型，具体方案设置见表1.采用有限体积方法对控制方程进行离散，压力-速度耦合基于Simple算法.RANS方法中控制方程的各项均采用二阶迎风格式，DES方法中压力项采用二阶迎风格式，动量项采用中心差分格式.

表1 计算方案Tab.1 Calculation plan

2 计算结果及分析

2.1 湍流模型及网格对RANS方法计算结果的影响

图3 给出了不同湍流模型下DU97-W-300 翼型的计算结果.图中，Cl为升力系数，Cd为阻力系数α为攻角.从图3 可以看出，在翼型失速之前，TSST 湍流模型更能准确地预测翼型的气动性能，升力系数、阻力系数及升阻比更接近实验值.但在翼型失速之后，全湍流模型（SA 和SST 模型）比T-SST湍流模型计算结果更接近实验值.对比可以发现，无论是失速之前还是失速之后，相比全湍流模型，T-SST湍流模型计算得到的升力系数偏大，阻力系数偏小，升阻比偏大.这是由于考虑转捩之后，叶片前缘部分流动变为层流，从而使得叶片升力系数偏大，阻力系数偏小，这与文献［12］和文献［13］中得出的结论一致.全湍流模型中SA 湍流模型计算效果要优于SST 湍流模型，计算值更加接近实验值.

图3 DU97-W-300翼型计算结果Fig.3 Computation result of DU97-W-300airfoil

图4（a）为DU97-W-300翼型表面压力系数Cp分布曲线.从图4（a）可以看出，在9.3°、12.4°攻角下，T-SST 湍流模型的计算结果与实验值更加吻合，SST 湍流模型计算效果最差.转捩是翼型绕流中重要的流动现象，转捩流动现象复杂，一直是CFD 模拟中的难点，翼型转捩位置的预测从一定程度上反映了CFD 的计算精度.随着对翼型绕流流动细节的研究，边界层转捩成为重要的研究对象.图4（b）为CFD 预测翼型转捩位置与实验值的对比曲线.从图4（b）可以看出，随着攻角的增大，翼型吸力面转捩位置向翼型前缘移动，压力面转捩位置向翼型尾缘移动.从计算转捩位置来看，T-SST 湍流模型计算得到的翼型压力面和吸力面转捩位置与实验值吻合较好，说明考虑转捩的T-SST 湍流模型能较好地预测翼型转捩位置.

图4 DU97-W-300翼型表面压力系数及转捩位置分布Fig.4 Distribution of Cpand transition position for DU97-W-300airfoil

由上面的分析可以看出，在翼型线性段，T-SST湍流模型能更好地预测叶片升力系数、阻力系数及转捩位置.笔者通过改变翼型一周网格数（方案A）及第一层网格高度ΔY（方案B）来研究网格对预测转捩位置准确性的影响.具体方案及计算结果见表2，计算攻角选为α＝12.4°.

图5为不同网格数及网格高度下DU97-W-300翼型表面压力系数及转捩位置分布图.从图5可以看出，方案A 中1、2、3、4计算所得Cp均与实验值吻合较好，且均能正确预测转捩位置，所以翼型一周网格数对预测转捩位置的影响较小.方案B 中5、6能较准确地预测出转捩位置，且Cp与实验值吻合较好，但方案B 中7、8未能计算出翼型表面的转捩位置，且Cp与实验值相差较大，由此可见对于预测转捩位置，第一层网格高度比较重要，当ΔY／C≥1.66×10-4（即y＋≥18）时已经不能预测出转捩位置.

表2 网格对预测转捩效果的影响Tab.2 Influence of mesh distribution on the predicted transition effects

图5 不同网格数及网格高度下DU97-W-300翼型表面压力系数及转捩位置分布Fig.5 Distribution of surface pressure coefficient and transition position for DU97-W-300airfoil at different grid number and height

2.2 DES与RANS方法计算结果的对比

图6给出了DU00-W2-401翼型DES与RANS方法计算结果的对比.其中计算攻角为5°～23°.RANS方法中湍流模型选择SA 湍流模型，DES方法中近壁区域湍流模型也选择SA 湍流模型.从图6可以看出，相比于RANS方法的计算结果，DES方法能更好地预测大厚度翼型的气动特性.相比实验值，在翼型失速段，RANS 方法计算所得升力系数偏大，阻力系数偏小，升阻比偏大.与DES方法计算结果相比，在翼型失速之后，RANS 方法计算所得升力系数偏大，阻力系数偏小，升阻比偏大.DES方法计算得到的翼型升力系数、阻力系数及升阻比与实验值吻合更好.从计算结果来看，在翼型失速之后，DES方法的模拟结果明显优于RANS方法.

图6 DU00-W2-401翼型计算结果Fig.6 Computation result of DU00-W2-400airfoil

图7 为DU00-W2-401 翼型表面压力系数分布.从图7可以看出，在15°、19°攻角时，RANS计算结果的吸力峰值均比DES方法的计算结果要高，所以计算所得气动力较DES方法的计算结果偏大.

图8为DU00-W2-401翼型中间对称面上涡量云图.从图8 可以看出，DES方法可以模拟出更细致的漩涡结构.从DES计算结果可以看出在分离区的上游有组织、有规律地脱落出分离涡，而在翼型的尾缘处，翼型上表面分离区内是低压区，翼型下表面是高压区，翼型下表面的流体会绕过叶片尾缘卷入分离区，这样分离涡周而复始向下游脱落.

图7 DU00-W2-401翼型表面压力系数分布Fig.7 Distribution of surface pressure coefficient

图8 涡量云图Fig.8 Cloud map of vorticity

图9 为带速度云图的Q等值面（Q＝1 000 s-2）.从图9可以看出，随着攻角的增大，翼型上表面的湍流分离涡尺度越来越大，分离区尺度越来越大.由于大厚度翼型极易发生流动分离，分离之后RANS方法不能模拟出细致的漩涡结构，导致其预测气动性能不准确，而DES方法能较好地模拟出细致的湍流脉动漩涡，所以DES 方法较RANS 方法更适合大厚度翼型的气动计算.

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

3 结 论

（1）对于RANS方法的计算，在翼型最大升力系数攻角之前，考虑转捩的T-SST 湍流模型的计算效果要优于全湍流模型，在翼型最大升力系数攻角之后，全湍流模型的计算效果要优于T-SST 湍流模型，全湍流模型中SA 湍流模型的计算效果优于SST 湍流模型.网格对计算结果及转捩预测的影响中，第一层网格高度更加敏感，而翼型一周网格数的影响不大.

（2）对比RANS与DES方法的计算结果发现，当叶片存在大尺度分离时，RANS方法不能模拟出更细致的分离涡结构，而DES方法则能在一定程度上模拟出细致规律的分离涡结构.对比气动性能计算结果，DES 方法计算结果与实验值更加吻合.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］.航空学报，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］高伟，李春，高月文，等.几何参数对风力机翼型转捩特性的影响［J］.动力工程学报，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］范忠瑶.风力机定常与非定常气动问题的数值模拟研究［D］.北京：华北电力大学，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.