基于BEM与CFD数值模拟的风力机叶片气动性能分析

史先路，赵 凡，赵 凯，，赵军华，3

(1.江南大学 机械工程学院，江苏 无锡 214122；2.中国空气动力研究与发展中心，四川 绵阳 621000；3.江苏省食品先进制造装备技术重点实验室，江苏 无锡 214122）

0 引言

风力机常年运行在恶劣环境下，叶片表面易表现出流动分离、失速等非定常气动特性，设计新一代高效、轻量化的风力涡轮机，需要更准确地预测气动载荷。

作为风力机叶片设计及其气动性能评估的常用方法，叶素动量理论(BEM）[1]由于其计算简单，能快速预测叶片气动性能而备受青睐。文献[2]利用BEM理论分析了风力机叶片气动载荷的分布规律以及风速大小对叶片气动载荷特性的影响。作为一种常用的数值模拟方法，计算流体动力学(CFD）可用于解决所有常见的流动问题。文献[3]，[4]利用CFD方法在均匀来流条件下对风力机叶片表面流动分离情况进行了研究，但未考虑大气边界层的影响。文献[5]基于CFD计算方法研究了一种小型风力机的气动特性，将计算结果与基于BEM计算方法得到的结果进行了对比，结果表明，两种计算方法得到的气动载荷特性相似。文献[6]基于BEM与CFD计算方法研究了风力机的功率和载荷特性，采用Spalart-Allmaras湍流模型针对海上变速变桨风力机分析了不同雷诺数下两种计算方法的结果差异，但未考虑失速状态下叶片表面的流动分离对计算结果的影响。

研究人员利用BEM及CFD数值模拟方法对风力机叶片的气动特性进行了较为深入的研究，但是，基于BEM与CFD两种计算方法对风力机叶片气动载荷的对比研究较少，且集中在两叶片风力机或者功率较小的小型风力机。本文采用修正的BEM和不同来流条件下CFD数值模拟两种方法研究大型风力机三维旋转效应导致的叶尖损失、动态失速等现象，分析对比两种风力机气动性能计算方法，为风力机的设计、功率预测和气动载荷的研究以及BEM的进一步修正提供理论支撑。

1 气动性能计算方法

1.1 修正后的叶素动量理论

BEM通过对作用在叶片展向上每个微段上的力进行积分得到作用于风轮上的推力及转矩。图1为叶片剖面和气流角受力关系图。图中：入流角φ为总速度W与风轮转动平面的夹角，由局部攻角(α）和叶素桨距角(β）组成，变量V1为入流风速，Ω为风轮转速，a为轴向诱导因子，b为周向诱导因子。根据作用在叶片截面上的升力(L）和阻力(D）的关系可以得到叶素上气动合力的轴向分量和切向分量。

图1 叶片剖面和气流角受力关系Fig.1 Force relation diagram of blade profile and flow angle

作用在叶片单元上的法向力和转矩可以通过升力系数(CL）和阻力系数(CD）的关系得到，假设叶片单元的弦长(c）已知，叶片数为N，则法向力dFA和转矩dT可以表示为

由于模型中假设每个叶片截面与相邻截面之间没有相互作用，使得该模型无法捕捉叶片上发生的三维效应，尤其是在叶尖和叶根附近。此外，气流绕风轮叶片剖面流动使得叶片表面气流分离，进而导致梢部和根部翼型表面的环量减少，相应计算的转矩也会减小。因此，利用BEM计算时需要考虑修正因子，通常使用文献[7]的叶根、叶尖损失修正因子来校正剖面数据。

式中：rn为轮毂半径；Ft为叶尖损失修正因子；Fr为叶根损失修正因子。

令叶尖叶根修正因子F=Ft•Fr，则修正后的法向力和转矩可以表示为

文献[8]在文献[7]的修正因子的基础上，引入一个新的修正项用于法向力系数Cn和切向力系数Ct，并将其代入经典的叶素理论中得到推力与转矩。

同 理，令F′=Ft′•F′，其 中g=exp[-c1(Nλ-c2）]，c1=0.125，c2=21，则修正后的法向力和转矩可以表示为

1.2 CFD数值模拟

1.2.1建立计算模型及网格无关性验证

以某1.5 MW大型水平轴风力机为研究对象，选取风机叶片常用翼型S809，设计功率为1.5 MW，风 速V为12 m/s，风 轮 转 速 为22.6 r/min，直径D为71 m。设置内域为旋转域，外域为静止域，划分非结构化四面体网格并对叶片附近以及内流域区域进行网格加密，当网格数量大于960万时，功率误差变化很小，因此最终确定数值计算的总网格数为960.4万。网格无关性验证见表1。

表1 网格无关性验证Table 1 Validation of grid independence

将模型导入求解器中进行求解，定义湍流模型为SST k-ω，设置速度入口和压力出口。压力速度耦合采用SIMPLE算法，离散方式为二阶迎风格式。将速度入口边界条件以UDF的形式加载到求解器中分别进行均匀风场和非定常流动风场的风力机三维数值模拟，在叶尖和叶根处设置压力监测点，监测叶片表面压力随时间的变化。压力监测点位置如图2所示。

图2 压力监测点位置Fig.2 Location of the pressure monitoring point

1.2.2湍流风场的模拟方法

由于风力机叶片在旋转过程中会发生流动分离、失速延迟等现象，而复杂的工作环境，例如大气边界层、湍流、风剪切、风速和风向的瞬变性以及偏航等会使翼型及叶片的气动行为呈现强烈的非定常特性，目前大多数关于水平轴风力机的研究没有考虑到这种不稳定性。

非定常流的模型复杂，模拟真实脉动风需要大量计算资源与时间，在理想情况下，可以将高频湍流风视为均匀来流、剪切来流以及湍流，其中湍流风可以被理想化为随时间变化的具有平均值、幅值和频率的正弦曲线形式的风速序列。此外，文献[9]的研究结果表明，正弦函数形式的脉动风速序列非常适合替代湍流来研究非定常流动行为，因此，本文使用式(11）表示某一局部时间点的瞬时风速。

式 中：vˉ(z）为 距 地 面 高 度z处 的 平 均 风 速；k为 湍流强度，为风速的标准差与平均风速之比；f为风的脉动频率；t为流动时间。

为了突出脉动风的波动特征，对式(11）中的相应参数进行指定，来设置较高的频率和振幅，其中f取1 Hz，k取0.1，z垂直于地面方向向上。

平均风速随高度的变化可以用指数率分布来描述。

式 中：γ为 切 变 系 数 或 粗 糙 度 指 数，取0.27；vˉ(z0）为距地面高度z0=10 m处的平均风速，取6.4 m/s。

2 结果与讨论

2.1 翼型气动性能分析

在进行三维数值模拟前，先分析S809翼型的气动性能，以验证湍流模型的准确性，俄亥俄州立大学(OSU）对翼型S809进行了不同雷诺数条件下的风洞试验[10]。本文为了改善失速状态下数值模拟的准确性，在均匀来流条件下使用SST k-ω湍流模型，选取封闭常数 β*作为调整参数。β*会改变湍动能k和湍流耗散率ω，在翼型数值模拟中校正湍流模型的封闭常数并用于叶片的数值模拟。获得不同攻角(α）下翼型的CL和CD，并将仿真结果与OSU的二维静态测试数据和BEM计算结果进行对比(图3）。

图3 不同攻角下翼型升力系数对比Fig.3 Comparison of airfoil lift coefficients at different angles of attack

由图3可知：在 α为0～8 °时，CL呈线性增长，在 α为8～16°时，CL趋于平缓，在 α超过16°后，CL急剧下降，调整封闭参数后的CFD模拟结果和实验数据几乎完全吻合；在α>8°时，BEM计算结果与CFD模拟结果和实验结果差异较大，表明BEM不能很好地预测翼型在失速状态下的气动性能，需要进一步修正。

图4为不同 α下的翼型流线图。由图4可知，随着 α的增大，失速区域逐渐增大且分离点逐渐向前缘移动，可将3个阶段称为附着流动状态(α为0～8°）、轻 失 速 状 态(α为8～16°）和 深 度失 速 状 态(α>16°）。

图4 不同攻角下翼型流线图Fig.4 Flow diagram of airfoil at different angles of attack

2.2 叶片表面极限流线图

图5为均匀来流条件下不同叶尖速比(TSR）下的极限流线。当风速分别为12，14，16.8 m/s和21 m/s时，对 应 的TSR分 别 为7，6，5和4。在 离 心力和叶片展向压力梯度的作用下，气流在叶根处发生分离并沿叶片展向流动，展向压力梯度随着攻角和叶尖速比的变化而变化。

图5 叶片表面极限流线Fig.5 Limiting streamline of blade surface

由图5可知：当风速为12 m/s时，所有流线均从前缘指向后缘，叶片压力面均为附着流动，没有发生失速，除叶根小部分区域外，叶片吸力面基本未发生分离流动；随着风速的增加，叶根部位分离点向前缘移动，叶片吸力面分离区域逐渐增大，叶片展向有明显的流动分离线；当风速增大到21 m/s时，叶片处于深度失速状态，吸力面流动完全分离，沿叶片展向有速度分量，在靠近叶根部位产生一个较大的气流漩涡。

2.3 叶片表面压力频谱特性分析

为了进一步研究非定常来流对风力机叶片表面压力波动的影响，对各叶尖和叶根部位截面监测点所测压力进行FFT变换，得到压力功率谱(图6）。压力功率谱总体上符合Kolmogorov-5/3功率定律，所有截面在频率f=0.376 Hz时均出现明显的峰值，而且该频率与风轮的旋转频率相等，风轮的旋转频率fn=n/60(n为风轮转速），说明压力脉动的主频为1倍风轮转动频率，即f=1fn。

图6 翼型压力的功率谱特性Fig.6 Power spectrum characteristics of airfoil pressure

由图6(a）可知：当f=1fn时，所有测点压力功率谱均出现明显的峰值，峰值从大到小依次是x/c等 于0.85，0.11，0.32和0.57，靠 近 后 缘 的x/c为0.85时的功率谱峰值最大，说明叶片后缘压力含有的主频能量最高；在高频段，靠近尾缘的x/c=0.85呈现明显的锯齿状分布特性，其余所有压力功率谱均呈现紊乱状态，表明叶尖部位翼型吸力面对失速旋涡的敏感程度较高，尤其是尾缘部位。

由图6(b）可知：当f=1fn时，所有测点压力功率谱均出现明显的峰值，峰值从大到小依次是x/c等于0.11，0.32，0.85和0.57，最接近前缘的x/c等于0.11的功率谱峰值最大，表明前缘对来流最敏感，后缘次之而翼型中部最弱；在频率f为2～11 Hz时，功率谱峰值点不明显，表明叶尖部位翼型压力面对失速旋涡的敏感程度较低。

由于叶根部位翼型受转频影响较小，使得叶根在2～5倍转频的功率谱峰值不明显[图6(c），(d）]，但总体上符合Kolmogorov-5/3功率定律。

总体而言，压力功率谱在高频段均呈现相对平缓的变化趋势，而且压力面压力功率谱的平缓变化趋势更明显，翼型越接近叶根，平缓变化的功率谱所占频段越宽；由于靠近叶根翼型的攻角较大，吸力面受到顺压梯度作用较强，流动分离严重，使得吸力面靠近叶根部位的功率谱的平缓变化区域较小。

2.4 BEM与CFD叶片载荷预测结果对比

本文使用QBlade进行BEM计算来预测叶片气动载荷，考虑了失速延迟效应的影响，分别利用两种修正方法对BEM进行修正并将修正后的结果与CFD计算结果进行对比。将叶片弦长、扭角等参数导入到QBlade中完成叶片模型的建立并进行气动性能分析，然后与文献[11]中的CFD模拟计算结果进行对比(图7）。

图7 叶片截面力对比Fig.7 Comparison of blade cross section force

由图7可知：两种不同风力机模型计算得到的气动载荷趋势一致，沿着展向叶片所受Fn逐渐增大，叶尖损失的存在使得在接近叶尖时Fn和Ft均有所减小，同时也从侧面证明了本文叶片设计与计算结果的合理性；Fn和Ft在叶根处修正后的BEM与原始的BEM计算结果差别不大，而在叶尖处修正过后更加接近CFD模拟结果；相较于文献[7]修正的BEM，文献[8]修正的BEM对于Fn的预测更加准确，对于Ft，二者预测结果差别不大。整体而言，修正后的BEM能够较好地预测出气动载荷随着叶片展向位置变化的趋势，但是与CFD计算结果相比仍存在一定偏差。

3 结论

本文采用修正的BEM和三维CFD数值模拟两种计算方法分析预测了1.5 MW风力机的气动性能。考虑了叶尖损失、失速延迟效应以及均匀来流和非稳态来流两种流动风场。研究正弦振荡和剪切流场下由风力机三维旋转效应产生的流动分离和动态失速效应，分析了脉动风条件下叶片表面压力频谱特性，最后利用QBlade软件建立优化设计的风力机模型，将计算结果与CFD模拟结果进行对比，得出以下结论。

①标准状况下，叶片表面基本为附着流动，在叶根和叶尖会产生分离涡。随着风速的增加，在离心力和叶片展向压力梯度的作用下，叶片吸力面开始发生流动分离，深失速区域逐渐增大，风速增大到21 m/s时，吸力面气流完全分离。

②通过叶片展向上不同位置截面压力功率谱可以看出，测点压力在1倍风轮转频时出现峰值且前缘峰值最高，由于叶尖涡的影响，叶尖吸力面后缘峰值最高；叶尖部位翼型截面吸力面对失速旋涡的敏感程度较高，压力面敏感程度较低，叶根部位翼型受转频影响较小。