CFD在水平轴风力机研究中的运用文献调研

风力机叶片气动弹性分析的气动理论模型常用的有3种：1，叶素动量理论

（BBM）；2，动力失速模型；3，3D－CFD【１】。

最近关于采用求解N-S方程的CFD方法研究风力机翼型气动性能的主要工作

集中在翼型的静态失速和动态失速方面。从物理意义上说,数值求解

Navier2Stokes方程的3D CFD 方法是最能全面准确地计算风力机气动特性的方

法。美国科罗拉多博耳德的国家可再生能源实验室于2000年12月首次将3D CFD

方法应用于NREL Phase2V I风力机气动性能研究,通过与风洞实验数据比较,结果

表明即使在高度的三维效果和深度失速状态下, 3D CFD数值模拟方法与风洞实

验数据也高度吻合【２，３】。但由于其计算工作量大、紊流模型的选取复杂等

原因,目前3D CFD方法还远不能作为风力机气动设计和研究的日常工具，所以２

D CFD方法是目前风力机研究的主流手段。

文献【４】根据当前计算流体动力学(CFD)数值模拟的研究，对风力机专用

翼型Ｓ816模型进行数值模拟，主要针对翼型的升阻特性、失速、流动分离等方

面内容，研究了翼型升阻力系数随迎角的变化、翼型在非定常、定常数值模拟结

果对比等因素对翼型气动特性的影响，总结出一些规律。

具体实现：

1 几何模型

采用FLUENT软件的前处理软件GAMBIT进行几何建模。取弦长为1m的几何

模型，网格为结构化C型网格，网格数50100。

图1 S816局部网格

2 控制方程及湍流模型

二维不可压；RNG k-8两方程湍流模型，与标准K-E模型相似，但比其精度更高。

3 边界条件

进口：设定无穷远来流风速做为进口的边界条件，压力为大气压，进口气流的湍

流度可根据具体风场状况确定，这里选用常规的3％来计算。

出口：采用自由出流作为出口边界条件，但必须保证足够的计算域。

壁面：将翼型表面设为壁面，在近表面处要使用无滑移条件和无渗透条件。

4 二维非定常计算

Re=4 000 000

5 结果与分析

图2 S816升、阻力系数随攻角的变化

从图中可以看出，该翼型具有较高的升阻比，失速点出现在14。迎角附近，最大升阻比大约为46．5，出现在6度左右；而6度之后开始出了流动分离，升力系数与实验开始出现偏差，随着分离的增大，差别也增大。对存在较大分离的流动模拟不够精准也是二维模拟的缺陷所在。

图3 S816定常与非定常计算结果对比

定常模拟结果要比非定常模拟有偏差，主要表现在定常模拟结果升力系数偏大而阻力系数偏小。

图4 S816(一16度迎角)时的速度矢量图

图5 S816(8迎角)时的速度矢量图

图6 S816(16度迎角)时的速度矢量图

图7 S816(24。迎角)时的速度矢量图

主要结论：

１S816翼型是300 kW以上大中型风力机所采用的翼型，翼型的气动参数比较缺乏；

２随着迎角的增大分离点前移，涡的范围和强度增大；迎角为负角度并进一步减小时，翼型的下表面发生流动分离，随着迎角的进一步减小，分离点前移，涡的范围增大，强度增大。

文献【5】探讨了弯度、最大厚度及其位置、前缘粗糙度和雷诺数等因素对翼型气动性能的影响机理，并开发出了在Re数为3 000 000时CAS-W1-XXX系列水平轴风力机专用翼型族，期望获得最大升阻比、最大升力系数和降低前缘粗糙度的敏感性。

几点结论：

前缘半径的影响：前缘半径对翼型的气动特性影响显著,尤其是对最大升力系数和最大升阻比。"当翼型其它的几何参数保持不变时,前缘半径较大,翼型的最

大升力系数也较大；前缘半径较小时,气流会在较靠近前缘的地方发生从层流到湍流的转挨。

.最大厚度及其位置的影响：最大厚度及其位置对翼型的最大升力系数影响显著。最大厚度的位置距离前缘较近时,最大升力系数较大,但这会造成边界层层流段的长度较小,湍流段的长度较大,使得翼型的阻力系数较大"适当加大最大厚度的位置和前缘之间的距离,可以有效加大层流段的长度,降低翼型的阻力系数"

弯度的影响：当其它几何参数保持不变时,随着弯度的增加,翼型吸力面的负压峰值逐渐增大,升力系数逐渐增加,但增加的程度逐渐减小;阻力随弯度的增加而增加;最大升阻比的变化趋势依据不同的翼型而表现出不同的规律。

表面粗糙度的影响:叶片表面粗糙度的增加,这会导致翼型边界层提前发生转抿,边界层动量厚度增加,尾缘逆压梯度加大,引发气流较早分离,继而导致翼型提前发生失速,升力线的斜率下降,最大升力系数下降.

雷诺数的影响：雷诺数较小时，对雷诺数的小小改变就会非常敏感，而对雷诺数较大时则不明显。

目前广泛使用的风力机叶片翼型族有两大类：1，航空翼型广泛使用在水平轴风力机上的翼型族包括两大类,一类是传统的航空翼型族,主要是NACA翼型族,如NACA44XX、NACA23XXX、NACA63一XXX、NACA64-XXX、NACA65一XXX等翼型族;

2，上世纪八、九十年代开始出现的风力机专用翼型族,如NRELS系列、DU 系列、RsiΦ形50系列、FFA系列等翼型族"

用XFOIL数值模拟的可靠性的结论：

低雷诺数(Re

XFOIL的计算结果和实验结果存在的偏差较大"升力特性的预测结果偏大,表现为升力线的斜率较大,失速延迟,最大升力系数较大;阻力特性的预测结果比实验结果偏小,偏差随雷诺数的增大而增大,但平均偏差一般保持在18%之内;

.低雷诺数(Re<1,000,000)情况下,对于预测层流翼型(以NAcA64一418为例)的气动特性,XFOIL预测升力特性的结果良好,和实验结果较为吻合"阻力特性的预测结果的准确度和雷诺数有一定关系,雷诺数较小时(如Re=300,000),和实验结果较为一致;

雷诺数较大时(如Re=700,000)和实验结果有一定偏差,偏差最大达到了25%;

.雷诺数较大时(Re>l,000,000),重点讨论了XFOIL对风力机专用翼型的预测结果准确度问题。"较之低雷诺数的情况,预测结果的准确度大为提高,除失速区大攻角下的预测结果和实验结果有出入外,升力特性的预测结果和实验结果基本保持一致,阻力特性的预测结果仍比实验结果偏小,但平均偏差一般保持在18%之内: .固定转抿时,XFOIL基本上可以准确预测升力特性,但阻力特性的预测准确性较低,一般比实验结果偏大,最高偏差可达40%,因而计算结果只能作为定性的参考"

NUMECA软件模拟结果

根据边界层理论,边界层分为粘性底层、缓冲层和对数率层,粘性底层紧挨着固壁。低雷诺数湍流模型可较好地描述壁面粘性底层的剪切应力和换热系数,对网格密度的要求较高,Y+值一般需在1一10之间;高雷诺数湍流模型可较好地描述对数率层,Y+值一般需在20一50之间"。

S809翼型

计算模型

计算域网格在网格生成器IGG中生成,网格类型为C型拓扑结构化网格,网格节点数为297x49个,翼型表面的网格节点为201个,在翼型的前缘和尾缘处进行了网格的局部加密,翼型弦长为1,计算域为10倍弦长,壁面Y+值取为1左右,壁面最小网格间距为2x10- 5"网格正交性最小角度56度,延展比为2.6,长宽比为1600

网格局部图

数值模拟时采用了四种湍流模型,分别是0方程的B一L模型，一方程的S一A 模型，二方程的k一ε(extendedwallfunetion)模型和k一ε(Yang一Shih)模型。使用这些湍流模型时分别计算了来流攻角为时的翼型特性

结果分析及结论

攻角较小时,四种湍流模型对升力特性以及俯仰力矩特性的预测良好,基本和实验结果保持一致;对阻力特性预测不佳,两种k一ε模型的计算结果和实验值之间的偏差最大。对阻力特性预测不佳的原因在于四种模型都是全湍流模型,认为整个翼型表面都处于湍流状态,而实验结果显示吸力面的流动有相当一部分仍为层流流动,按照全湍流模型进行计算显然会高估阻力系数"

攻角较大、未发生失速时,翼型表面的流动基本都是湍流流动,四种全湍流模型得到的计算结果对升力，阻力以及俯仰力矩特性的预测良好,基本和实验结果保持一致；

发生失速以后,翼型表面的气流出现较大分离,四种湍流模型预测的气流分离点的位置比试验值靠后,计算值和实验值有较大偏差,升力系数较实验值偏大,阻力系数偏小；

影响计算结果精度的因素很多,对于相对厚度较大!前缘曲率半径较小的翼型,层流和湍流混合模型!转挨点和分离点位置的准确预估是提高数值预测结果准确度的关键。

两种数值方法的比较：

各有好坏，没有统一的规律。

低雷诺数时，B-L和S-A湍流模型的精度高于k-E湍流模型；

大攻角下，使用k-w（sst)模型得到的结果精度比SA湍流模型要高，对气流

分离点的位置预测较准。

雷诺数小于1,000,000时, XFOIL对NACA湍流翼型的预测结果和实验结果偏

差较大,对层流翼型的预测结果较为理想;雷诺数大于100,000时,xFOIL的预测结

果良好,但对吸力面上气流分离点的位置预测不准确;

文献【6】用CFD软件Fluent 6．0，S-A湍流模型对NRELS809翼型的二维动态流场进行了数值模拟，得到了翼型攻角在9度～31度范围内按正弦周期变化时的绕流流场。计算结果显示：动态失速下翼型的绕流流场与相同工况下的静态绕流流场有着十分明显的差别，同时也引起翼型升力、阻力系数的显著变化。

结论：

1)绕动态失速翼型的流场结构比较复杂，翼型动态失速与静态失速在流动特性上有很大区别。翼型在上仰过程出现的失速主要由翼型后缘分离引起，分离涡的影响范围较小，主要是在后缘附近，失速后翼型尾部压差变化不大可以充分说明这一点。当失速进一步深化，翼型前缘出现分离涡，逐步后移，形成大尺度分离流动，并在后缘诱发二次涡，从而引起翼型升力、阻力系数显著变化，呈现出较为明显的动态特性。

2)翼型在动态失速工况下，其升力系数比静态失速有很大提高，由此可以推出，当水平轴风力机在动态失速工况下运行时，叶片得到的升力将大大高于额定工况，使风力机输出功率远远超过额定功率。如果在风力机的设计过程中，缺乏对动态失速工况的考虑，很可能造成对风力机的性能预估不足。．

3)由于动态失速通常在来流状态变化的情况下发生，包括紊流流动、塔架阴影和偏航来流等，这些决定了动态失速的复杂性

文献【7】本文以 1.5MW 风力机叶片的六种翼型为研究对象，通过Gambit 提取叶片几何结构，运用CFD 软件建立三维模型，并模拟计算了五种工况下叶片的气动特性。计算结果以压力云图、流线图、速度矢量图、三维流场图等直观形象的形式描述了风力机流场的规律。风力发电机叶片气动性能主要评价指标有：功率系数、力矩系数、推力系数、叶尖速比等

FINE TM/ Turbo 中的Autogrid5进行几何和网格的处理

S-A湍流模型

文献【8】用Fluent研究了风轮附近的流场，以等值线图、云图、流动矢量、流线、迹线、二维曲线图等方式较为形象直观地描述了流场运动规律。主要控制参数：升阻力系数和部分流场图。

对NACA4415翼型进行了2D研究。用的三角形非结构化网格。近壁面第一层网格距离10-4. 。

文献【9】用FINE TM/ Turbo中的Auto grid5建立模型并针对叶片表面灰尘厚度对风力机气动特性的影响进行了数值研究。结果表明，叶片表面积灰厚度对风力机叶片气动特性有较大影响：随着积灰厚度的增加，叶片两侧大部分区域风速减小，气流运动无规律性增强，压力面与吸力面压差减小，风力机效率降低。

S-A湍流模型和k-w湍流模型。所有网格总数大约为200万如图1所示，其中翼型弦向网格数为135，叶片展向网格数为67；叶片表面大部分区域Y 小于2。

边界：进口边界给定速度分量和大气温度，出口边界给定大气压力，叶片表面是无滑移边界，轮毂表面是滑移边界。残差下降3个数量级及以上作为收敛条件。

比较了压力分布和诱导速度分布。

文献【10】采用带转捩修正的k－ωSST湍流模型对失速控制型NREL Phase VI 风轮在5~10 m/s的范围内，转速为固定的72 r/min的无偏航工况进行了数值模拟。计算得到了各个风速下的风轮功率、翼展方向5个截面的压力系数和叶片吸力面的极限流线的分布，并与NASA Ames风洞试验结果以及无失速延迟修正的升力面方法得到的结果进行对比分析。分析表明：CFD方法能够比较准确地预测风力机的功率以及展向空气动力载荷，在风力机设计和性能预测领域将得到更广泛的应用。

用的FLUENT软件计算。计算区域的整场结构化网格采用专业网格生成器ANSYS ICEM生成，如图3所示。整场网格总数约为2×106个网格单元。叶片表面第一层网格在保证壁面y+<1.

结论：在风速较低的情况下，流动基本为附着流时，利用升力面方法和带有过渡流修正的k－ωSST湍流模型能够得到和试验数据吻合较好的三维旋转叶轮的气动特性，即使在有较弱的失速存在情况下，三维CFD 方法也能获得比较准确的结果。但是在较高风速的情况下，流动出现复杂的三维分离流动的时候，CFD 预测值和试验值之间也存在一定的差距

李栋等【11】（2005）应用DES(Detached-Eddy Simulation)方法数值模拟了3种不同失速类型的翼型的升力特性，对此模型使用了LU-SGS隐式格式求解，和实验结果对比，显示了这种方法可以有效地预测翼型的失速特性。2006年李栋

【12】等又对比了RANS和DES方法数值模拟翼型失速特性的能力，发现对大分离流动的数值模拟，DES方法体现出更强的能力。

动态失速方面的CFD研究主要有：

钱炜祺等【13】（2001）分别采用k-ω模型和k-ωSST湍流模型对翼型的动态失速进行了数值模拟。振荡翼型轻失速和深失速算例的计算结果表明：(1)绕动态失速翼型的流场结构十分复杂，轻失速和深失速在流动特性上有很大区别。计算结果显示轻失速主要是由于后缘分离引起，分离涡的影响范围主要是在后缘附近。而深失速则首先形成很大的前缘分离涡，该分离涡在翼型表面上运动，并诱发出二次分离涡，引起翼型升、阻力系数的显著变化。(2)对于动态失速的翼型绕流，

k-ωSST湍流模式是较为有效的，计算出的气动力系数曲线变化趋势与实验结果符合得比较好。

陈旭等

【14】（2003）应用CFD软件Fluent6.0 求解N-S方程和S-A湍流模型对NRELS809翼型的二维动态流场进行了数值模拟，得到了翼型攻角在9°～31°范围内按正弦周期变化时的绕流流场。计算结果显示：动态失速下翼型的绕流流场与相同工况下的静态绕流流场有着十分明显的差别，同时也引起翼型升力、阻力系数的显著变化。

【15】F. Bertagnolio等（2005）分别采用二维和三维的Navier-Stocks方程并结合SSTk-ω和DES湍流模型以及半经验动态失速Beddoes–Leishman模型对风力机翼型RIS?-B1-18在静止和做俯仰震动两种条件下的绕流进行了计算，二维RANS 和k-ωSST湍流已经足够精确来模拟附着的流动，流动分离时，差距巨大，三维DES模型能够给出较好的结果，在叶片静止条件下比实验预测的失速攻角要小，三维RANS模型同样能得到较好的结果，流场即使稳定时也会表现出三维特性。

风力机叶轮气动方面，全湍流计算的主要有EllipSys3D和混合算法的Georgia Tech 两种软件。

文献【16】用Gambit画叶片的几何，用FLUent6.2计算的，对湍流分别采用Spalart-Allmaras湍流模型和k-ωSST湍流模型进行处理，计算区域离散采用四边形或六面体单元进行结构化网格划分，壁面处湍流采用增强壁面函数法进行处理。对翼型进行数值模拟时，计算模型的攻角的变化范围为1°至20°，来流的雷诺数和马赫数分别为1×106和0.07333。数值计算结果与来自DTU的实验数据进行了比较，分析了翼型边界层分离流动现象，统计了升阻力系数、流场中的压力截面、速度矢量图，并得出了三维模拟在有大分离时效果比二维的要好，且在三维模型下采用Spalart-Allmaras湍流模型时翼型失速性能的数值计算结果较为准确。

三维翼型网格

文献【17】对水平轴风力机CFD计算的湍流模型进行了讨论。对NREL 第6 水平轴风力机在风速13 m/ s 情况下,采用定常k-ωSST 和非定常k-ωSST 以及分离流DES 3 种湍流模型进行了三维旋转流场数值模拟,计算结果与NASA 非稳定气动力学风洞试验结果进行了对比分析,计算结果表明非定常k-ωSST 湍流模型计算结果比，分离流DES 湍流模型计算结果与试验结果在大部分位置更接近. 可能原因是DES对小分离计算可能并不是太适合，也有可能作者没有对DES进行验证，采用非定常k-ωSST 湍流模型可以满足水平轴风力机流场计算要求.

用的ICEM化分的网格，FLuent计算的。

文献【18】作者基于最大能量系数为目标函数，研发了适用于水轮机的新型叶片，并用数值模拟方法进行了研究，结果表示：最大能量系数提高了5.2%。

A CAD image of GEM hydro-turbine.

用SCoPE进行的数值模拟。

Mesh around the arifoil

文献【19】对水平轴风力机边界层的旋转效应进行了数值研究，计算模型为S806

翼型和NREL Phase VI风力机，对3D叶轮和2D翼型用Ansys ICEM 12.0划分网格，

Fluent进行数值计算。计算结果显示：相对于2D计算，3D叶片上的失速推迟了，原因是随着旋转速度的增加或者随着展向位置的减少，分离点推迟了。

对于2D翼型，采用C型网格，共55200个网格点。用的K-w SST湍流模型

2D计算网格

3D计算拓扑和计算域

3D结构化网格

3D非结构化网格

文献【20】提供了可用于风力机上的翼型的目录，将用EllipSys2D计算出的结果与xfoil计算出的2维结果以及实验结果对比，为风力机叶片的设计提供参考。

k-ωSST湍流模型，第一层网格为翼型弦长的10-5，用随攻角变化的升力、阻力和俯仰力矩系数，以及压力与表面粗糙度分布来考察气动特性。

文献【21】作者用叶素理论（BEM）来找到最佳的理论值。给出了不同叶片角度、升阻比、升力系数、叶片实度和叶尖速比对能量系数的影响，并获得了最佳设计值。最佳设计值为叶片实度0.15，攻角5度，叶尖速比8，阻力系数与升力系数的比为0.025。

文献【22】作者针对NREL Phase VI系列风力机，用SST-k-w湍流模型和Langtry-Menter transitional models对8个工况进行模拟，将推力、流线和压力系数与实验值进行对比，分析了近壁面网格（弦长和径向网格数）对计算结果的影响和水平轴风力机的气动特性。结果表明sst-k-w模型过早了预测了流动分离，Gamma–Theta transitional 模型更贴近实验值。用的Ansys ICEM 11划分的网格，Ansys CFX11计算的。

文献【23】用OpenFOAM对一个三叶片小型风力机进行了数值模拟，通过与实验值对比，验证SST湍流模型。结果表明，能量利用与实验所测拟合得比较好。

参考文献

【１】杨树莲, 侯志强, 任勇生, 等. 风力机叶片气动弹性和颤振主动控制研究进展[J]. 机械设计, 2009 (9): 1-3.

【２】S? rensen N N, Michelsen J A, Schreck S. Navier2stokes p redictions of the NREL phase V I rotor in the NASA Ames 802by2120 wind tun2nel[C ] / /A

IAA2200220032, 2002: 94 - 105.

【３】Xu G, Sankar L N. Application of a viscous flow methodology to the NREL phase V I ROTOR[C] / /A IAA - 2002 - 0030, 2002: 83 - 93.

【４】杨从新, 李春辉, 巫发明. 水平轴风力机专用翼型的数值模拟研究[J]. 科学技术与工程, 2008, 8(17): 4994-4998.

【５】李宏利. 水平轴风力机专用翼型族的设计及其气动性能研究[D]. 北京: 中国科学院工程热物理研究所, 2009.

【6】陈旭, 郝辉, 田杰, 等. 水平轴风力机翼型动态失速特性的数值研究[J]. 太阳能学报, 2004, 24(6): 735-740.

【7】张耀华. 水平轴风力机叶片气动性能数值模拟[D]. 重庆大学, 2013.

【8】朱德臣. 水平轴风力机叶片附近区域流场的数值研究[J]. 内蒙古工业大学学报, 2007.

【9】吴小飞, 唐胜利, 周文平. 水平轴风力机叶片表面积灰厚度对三维气动特性影响的数值模拟[J]. 流体机械, 2012, 39(11): 32-37.

【10】俞国华, 杜朝辉. 水平轴风力机气动性能的三维数值计算[J]. 中国科技论文在线, 2009, 4(10): 752-757.

【11】李栋，焦予秦，Igor Men’shov，中村佳朗．Detached-Eddy Simulation方法模拟不同类型翼型的失速特性．航空学报，2005，26（4）：406-410 【12】李栋，Igor Men’shov，中村佳朗．翼型失速特性的RANS方法与DES方法数值模拟对比分析．西北工业大学学报，2006，24（2）：228-231 【13】钱炜祺，符松，蔡金狮．翼型动态失速的数值研究．空气动力学学报2001，19（4）：427-433

【14】陈旭，郝辉，田杰，杜朝辉．水平轴风力机翼型动态失速特性的数值研究．太阳能学报，2003，24（6）：735-740

【15】 F. Bertagnolio, N. N. S?rensen, F. Rasmussen. New Insight Into the Flow Around a WindTurbine Airfoil Section, Journal of Solar Energy Engineering,

Transactions of the ASME,2005, 127(): 214-222

【16】张义华, 李隆键. 水平轴风力机空气动力学数值模拟[D]. 重庆: 重庆大学, 2007.

【17】杨瑞, 李仁年, 张士昂, 等. 水平轴风力机CFD 计算湍流模型研究[J].

甘肃科学学报, 2009 (4).

【18】Daniele E, Coiro D P. Optimization of diffuser geometry for an horizontal axis shrouded hydroturbine[C]//Clean Electrical Power (ICCEP), 2013

International Conference on. IEEE, 2013: 240-247.

【19】Gao X, Hu J. Numerical simulation to the effect of rotation on blade boundary layer of horizontal axial wind turbine[C]//World Non-Grid-Connected Wind

Power and Energy Conference (WNWEC), 2010. IEEE, 2010: 1-4.

【20】Bertagnolio F, S?rensen N N, Johansen J, et al. Wind turbine airfoil catalogue[M]. 2001.

【21】Bavanish B, Thyagarajan K. Optimization of power coefficient on a horizontal axis wind turbine using bem theory[J]. Renewable and Sustainable Energy

Reviews, 2013, 26: 169-182.

【22】Moshfeghi M, Song Y J, Xie Y H. Effects of near-wall grid spacing on

SST-K-ω model using NREL Phase VI horizontal axis wind turbine[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2012, 107: 94-105.