1.5MW风力机整机定常及非定常数值模拟*

刘一帆，钟淋涓，杨　涛，黄树红，（.华中科技大学 中欧清洁与可再生能源学院，武汉 40074；2.河海大学 能源与电气学院，南京 200；.华中科技大学 能源与动力工程学院，武汉 40074）

1.5MW风力机整机定常及非定常数值模拟*

刘一帆1，钟淋涓2†，杨涛3，黄树红1，3
（1.华中科技大学 中欧清洁与可再生能源学院，武汉 430074；2.河海大学 能源与电气学院，南京 211100；3.华中科技大学 能源与动力工程学院，武汉 430074）

本文运用计算流体力学方法，以1.5 MW风力机为例，对风力机整机三维模型的空气动力学特性开展研究。针对三翼型风力机叶片，利用改进的Wilson方法进行气动设计，并通过寻找各截面最佳雷诺数的方法进行优化修正。建立了整机三维模型，设计流域并划分网格，定义边界及区域。最后对上述模型进行额定工况下定常与非定常数值模拟，利用模拟结果开展有关压力、失速特性等空气动力学特性的分析。结果表明：非定常模拟在风轮背面上的平均压力比定常小，使风轮前后压差变大，输出功率加大，其主要原因是叶尖出力的增加；旋转使得风力机叶片发生流体分离延迟，且产生更高的升力系数。

风力机；气动设计及优化；整机三维模型；数值模拟；空气动力学特性

0　引　言

能源是人类社会经济发展及文明进步的重要动力。进入21世纪以来，随着经济的迅速发展，能源消费也在逐年增加，然而常规能源资源正面临着日益枯竭的窘境，人类迫切需要清洁、无污染、可再生的新能源。风力发电作为新能源中最为成熟且最具有规模开发条件的能源利用方式之一，得到了广泛的重视。风力机是把风能转化为电能的关键设备，其气动特性直接决定了输出功率。

目前国内外对风力机气动特性的研究主要集中在叶片。巴西学者VAZ等［1］在低叶尖速比的情况下，结合叶素动量理论与涡流理论进行了风力机的叶片设计。德国的DÍAZ-CASÁS等［2］基于遗传算法对叶片进行了优化。此外，还有考虑实验结果［3］、不同风况［4］等影响因素的叶片优化设计。而已有的一些对风力机整机气动特性的研究也只考虑了定常情况下的数值模拟：韩中和等［5］考虑了剪切风的情况，对1.3 MW风力机整机在8 m/s和13 m/s来流风速下的绕流流场进行了全三维定常数值模拟；李少华等［6］建立了1.2 MW风力机三维模型，并通过数值模拟分析了其压力、速度及湍动能对风力机功率的影响。但在风力机实际运行中，空气的流动状态会随着流动的发展而发生改变，这种变化过程就需要非定常的数值模拟。因此，本文运用计算流体力学的研究方法，以1.5 MW风力机为例进行定常与非定常数值模拟，对MW级风力机整机三维模型的空气动力学特性开展研究。

1　风力机叶片气动设计

根据给定的来流风速和设计功率，本文应用叶素动量理论进行风力机叶片气动设计，确定风轮风能利用系数、叶尖速比、叶片数目、风轮直径、转速等整体参数。基本参数见表1。

表1　风轮基本参数Table 1　Basic parameters of wind turbine

风力机叶片使用三种翼型组合建立简化模型。其中叶片主体选择气动性能较优的NACA 4415，该翼型气动性能较好，升阻比较高；叶片的尖端部位则选择气动性能更好的NACA 4412；叶片根部则考虑到结构的强度问题，采用最大厚度更大的NACA 4418。

对于风力机整机来说，随着轮毂半径的增加，轮毂对叶片性能的影响也逐渐显现。因此本文根据PRANDTL等的研究结果，在Wilson方法上提出了修正，添加了轮毂损失，提出了考虑总的损失因子的设计方法［7］。研究表明，对于同一截面处，当轮毂半径变大时，同一倾角对应的轮毂损失系数就会变小，而且倾角越大，这种变化越大。因此，轮毂损失系数也可以用Prandtl叶尖损失修正模型表示。

式中：Ftip为叶尖损失系数；N为叶片数目；R为风轮半径；r为叶素在叶片展向的位置；φ为叶素倾角。

则轮毂损失系数：

式中，Rhub为轮毂半径。

因此，总的损失因子应为：

把Wilson设计法中的叶尖损失系数F换成用叶尖损失系数和轮毂损失系数乘积得到的总的损失系数F，其他理论同Wilson设计法。

在叶片的优化过程中，深入考虑了来流条件对于风力机气动性能的影响，寻找各截面的最佳雷诺数，进而得到最佳攻角以优化弦长与安装角。此外，多翼型的叶片在初始优化完成后往往不满足几何连续性的要求，这就需要进行进一步的气动外形修正，以达到提高线性度、便于加工和制造的目的。目前比较简单的修正方法是直接调节设计点来满足其几何连续性，但是这一方法不可避免地带来了系统误差。本文利用拟合的方法，可减小误差并能保证多翼型叶片的几何连续性。

经过校核验证，在额定工况，即风速为13 m/s时，理想设计条件下的风轮模型功率为1.467 9 MW，与设计值1.5 MW相比，误差仅为2.14％，可见风轮设计合理。

2　风力机整机几何建模与网格划分

2.1风力机整机建模

本文利用翼型设计软件Profili，将设计风轮所采用的三个翼型对应的形状参数和叶素在叶片展向的位置坐标导入Gambit中，生成三维立体，并通过旋转得到其余两个叶片。采用三叶片风轮常用的圆锥体型式固定式轮毂［8］，并通过半圆与圆柱拼接的方式简化轮毂模型。绘制直径为2.8 m的球体和长度为3.5 m的等直径圆柱，作为轮毂必不可少的两部分构件。并简化机舱为长8 m、宽4 m、高4 m的长方体模型，选择高度为65 m的塔筒，确定塔筒上下底面圆直径分别为2.4 m和4 m。最后得到风力机整机三维模型如图1所示。

图1　整机三维模型示意图Fig.1　Sketch map of the whole turbine’s 3D model

2.2风力机流域设计

为了更好地研究边界层对风轮性能的影响，选取一大一小两个流域，分别称为外流域和内流域。内流域主要为了研究风轮叶片三维边界层的影响，所以只需比叶片稍长，呈圆饼状包含整个风轮。因此设定内流域为半径35 m、厚5 m与风轮同心的圆饼。外流域主要研究三维尾涡系的影响，由于风力机安装在地面上，存在限制流体和固体的壁面区域，而模型在地面以上的区域则为具有镜像特征的外流域。因此，应采用沿Y轴方向拉伸的长方体型的外流域。同时，为了在模型网格划分的过程中加密靠近风轮模型处的外流域网格，把外流域划分为三个小长方体，分别为风轮前、风轮区、风轮后长方体。由于风力机尾流影响的范围在8倍的转轮直径后基本可以忽略［9］，故外流域具体尺寸选取如表2所示。

表2　外流域基本参数Table 2　Basic parameters of outer flow field

最终得到流域模型如图2所示。

图2　流域模型示意图Fig.2　Sketch map of the watershed model

2.3风力机网格划分

划分网格的目的是要对空间上连续的计算区域进行剖分，把它划分成许多个子区域，并确定每个区域中的节点，以便数值计算利用这些离散的点建立控制方程，求解得到离散方程组。本文利用Gambit中的网格划分算法，对所设计的模型和流域划分网格，如表3所示，最终得到3 550 432个计算网格。

表3　网格的划分Table 3　Grid division

3　计算方法与边界条件设定

对额定工况下的整机模型进行定常与非定常数值模拟。模拟均采用基于压力的求解器，这种求解器采用常规意义上的投影方法，首先通过动量方程求解速度场，进而通过压力方程的修正使速度场满足连续性条件。整个流场的模拟结果同时满足质量守恒和动量守恒，而动量方程与压力方程所组成的控制方程存在非线性和相互耦合作用，因此需要重复迭代多次直至收敛。

模拟过程是基于时均N-S 方程（RANS）进行的，紊流模型使用k-ω SST模型。这一模型考虑了低雷诺数、可压缩性和剪切流传播等情况，且不包含k-ε模型中所需要的复杂的非线性阻尼函数，因而具有更广泛的应用和更高的精度。

假设整机为刚体，即数值模拟过程不考虑整机表面的变形。入口为额定工况下的Y轴向自由来流，给定来流速度，出口设置为出流边界条件，设定压力速度耦合算法为SIMPLE算法，方程采用隐式格式，对流项差分格式全部采用二阶迎风差分格式。SIMPLE算法属于压力修正法的一种，是目前工程上应用最为广泛的流场计算方法，主要用于求解不可压缩流场。

设置外流域入口面边界类型为速度入口边界，出口面为压力出口边界，内、外流域接触面为接触面边界。叶片、轮毂、机舱、塔筒均为壁面边界。对外流域交界面设置为内部边界，流域地面为壁面，地面以上的外壁为对称边界。内、外流域指定类型为FLUID。

4　结果与分析

4.1压力分析

由叶片展向翼型表面的压力云图可知，非定常数值模拟在风轮背面的平均压力比定常模拟小，因此，风轮前后压差变大。而叶片翼型截面上的压力云图表明：在非定常模拟中，叶片沿展向上的压力变化在叶根和叶中处与定常模拟相同，而在叶尖位置则开始显现出较大的差别，如图3所示。

图3　叶片沿展向上的压力变化图Fig.3　Pressure variation of the blade along the development direction

由图3可知，非定常计算在叶尖翼型下翼面的压力大于定常情况，导致了压力差增加，进而升力增加，输出功率增大，即考虑流体时间特性的非定常计算功率加大的主要原因是叶尖出力的增加。

4.2失速特性分析

分析叶片失速的重点是翼型上表面的气流分离现象，图4给出了叶片吸力面的流线图。

图4　叶片吸力面流线变化图Fig.4　Streamline change diagram of blade suction surface

图4表明，在额定工况下，非定常数值模拟的叶片吸力面流体分离区更小，叶片中部与尖部的绝大部分翼型截面的前缘和后缘流动仍充分附着。由此可知：旋转使得风力机叶片吸力面相比于定常情况发生分离延迟。这种现象的出现有两个原因：一是由旋转产生的离心力使得边界层中的低能气体向叶尖流动；二是旋转时叶片上的科氏力使边界层中的气流得到一个向后缘方向流动的加速度，这个加速度会使气流边界层变薄［10］。这两个力均使得风力机叶片失速点向后移动，发生失速延迟。

4.3叶片升阻特性分析

为进一步研究叶片翼型的升阻特性，本文利用ANSYS CFD-Post软件自定义函数，对定常以及非定常模拟结果进行求解，得到叶片各翼型截面所受的轴向气动力dFx和切向气动力dFy，如表4所示。利用线性方程组的高斯法计算升力与阻力，进而得到升力系数与阻力系数沿翼型攻角的变化，即数值模拟情况下的升力特性曲线与阻力特性曲线，与Profili软件提供的二维翼型（NACA4418、NACA4415、NACA4412）升阻特性曲线进行对比，如图5与图6所示。

表4　额定工况下整机模型沿叶片展向气动力变化Table 4　Aerodynamic change of the whole turbine’s model along the direction of development under rated condition

图5　翼型升力特性曲线图Fig.5　Airfoil lift characteristic curve

图6　翼型阻力特性曲线图Fig.6　Airfoil drag characteristic curve

对比非定常与定常数值模拟下的升力特性曲线与阻力特性曲线，可知设计叶片攻角α在4°～12°的变化过程中，非定常数值模拟始终表现出略高的升力系数，这是由于旋转而造成的合理结果。但是在攻角α＜8°时，非定常情况下原本较低的阻力系数开始出现升高，甚至大于定常情况。这一现象的原因是在叶尖处旋转效应造成叶尖涡旋度增加，范围扩大，阻力增加。表明当叶片旋转时，翼型的气动性能与未旋转时相比发生了很大的变化，证明了非定常数值模拟的必要性。

对比二维翼型和三维翼型（由于整机模型的三维特性，数值模拟的叶片翼型可视作三维翼型）的升力特性曲线可知，攻角α在5°～9.5°时，二维翼型和三维翼型可以较好地吻合；但当攻角大于9.5°时，二维翼型的升力系数较三维翼型高，攻角越大，这一差别也越大。因此，从三维翼型实验模拟的结果来看，二维模型并不适用于翼型在较大攻角时的性能计算。此外，三维翼型在所有的攻角范围内都显示出了更大的阻力系数。

5　结　论

（1）非定常数值模拟在风轮背面上的平均压力比定常小，使风轮前后压差变大，输出功率加大，其主要原因是叶尖出力的增加。

（2）额定工况下非定常数值模拟的叶片吸力面流体分离区更小，叶片中部与尖部的绝大部分翼型截面的前缘和后缘流动仍充分附着，证明旋转使得风力机叶片相比于定常情况发生分离延迟。

（3）旋转效应导致较高的升力系数，但在叶尖处，也会造成叶尖涡旋度增加，范围扩大，使原本较低的阻力系数开始出现升高，阻力增加。

［1］VAZ J R P， PINHO J T， MESQUITA A L A.An extension of BEM method applied to horizontal-axis wind turbine design［J］.Renewable energy， 2011， 36（6）： 1734-1740.DOI： 10.1016/j.renene.2010.11.018.

［2］DÍAZ-CASÁS V， BECERRA J A， LOPEZ-PEÑA F， et al.Wind turbine design through evolutionary algorithms based on surrogate CFD methods［J］.Optimization and engineering，2013， 14（2）： 305-329.DOI： 10.1007/s11081-012-9187-1.

［3］LANZAFAME R， MESSINA M.Fluid dynamics wind turbine design： critical analysis， optimization and application of BEM theory［J］.Renewable energy， 2007，32（14）： 2291-2305.DOI： 10.1016/j.renene.2006.12.010.

［4］STORK C H J， BUTTERFIELD C P， HOLLEY W， et al.Wind conditions for wind turbine design proposals for revision of the IEC 1400-1 standard［J］.Journal of wind engineering and industrial aerodynamics， 1998， 74-76：443-454.DOI： 10.1016/S0167-6105（98）00040-3.

［5］韩中和.考虑风剪切的1.3 MW风力机整机三维定常流动数值研究［J］.动力工程学报， 2011， 31（10）：779-783， 808.

［6］李少华， 1.2 MW风力机整机流场的数值模拟［J］.动力工程学报， 2011， 31（7）： 551-556.

［7］廖明夫， 宋文萍， 王四季， 等.风力机设计理论与结构动力学［M］.西安： 西北工业大学出版社， 2014.

［8］丹麦RISØ国家实验室，挪威船级社.风力发电机组设计导则［M］.北京： 机械工业出版社， 2011.

［9］邓力， 李龙， 许昌， 等.基于CFD的水平轴风力机数值模拟［J］.电子测试， 2014（1）： 42-44.

Steady and Unsteady Numerical Simulation on Flow Field of a 1.5MW Wind Turbine

LIU Yi-fan1， ZHONG Lin-juan2， YANG Tao3， HUANG Shu-hong1，3
（1.China-EU Institute for Clean and Renewable Energy， Huazhong University of Science and Technology， Wuhan 430074， China；2.College of Energy and Electrical Engineering， Hohai University， Nanjing 211100， China；3.School of Energy and Power Engineering， Huazhong University of Science and Technology， Wuhan 430074， China）

This paper conducted a research on the aerodynamic characteristics of megawatt wind turbine’s 3D model using the Computational Fluid Dynamics（CFD）method.A 1.5 MW wind turbine’s 3D model was proposed as an example.For the blades including three airfoils， the improved Wilson’s method was used for aerodynamic design.The blades were optimized by finding the best Reynolds number of various sections.By using FLUENT software， a 3D model of the whole turbine was established.Flow field was designed and divided into grids， and boundary conditions were set.By using steady and unsteady CFD methods， the model under rated wind speed condition was calculated to study its pressure， stalling characteristics and other aerodynamic characteristics.The results show that： average pressure on the back of the rotor for unsteady simulation is lower than that for steady simulation， which is caused by tip output increase and increases the pressure difference and the output power of the rotor； rotation leads to fluid separation delay on the wind blades， and produces a higher lift coefficient.

wind turbine； aerodynamic design and optimization； turbine’s 3D model； numerical simulation； aerodynamic characteristic

TK83

A

10.3969/j.issn.2095-560X.2016.04.002

2095-560X（2016）04-0266-06

2016-06-11

2016-07-11

国家科技支撑计划（2015BAA06B02）

钟淋涓，E-mail：zhonglinjuan@hhu.edu.cn

［10］张义华.水平轴风力机空气动力学数值模拟［D］.重庆： 重庆大学， 2007. 10.7666/d.y1139007.

刘一帆（1993-），女，硕士研究生，主要从事基于CFD的流场数值模拟、风力机状态监测与故障诊断研究。

钟淋涓（1981-），女，博士，讲师，主要从事新能源规划、优化及技术经济研究。

杨涛（1978-），男，博士，副教授，主要从事发电设备控制与仿真、状态监测与故障诊断研究。

黄树红（1958-），男，博士，教授，主要从事火电厂设备强度、振动、状态预测与故障诊断及新能源与节能技术研究。