大型风力机变桨角度对尾流特性影响的研究

孙芳锦， 贺伟盼， 邓富河， 卢 琛， 张大明

(1. 桂林理工大学土木与建筑工程学院，广西 桂林 541006; 2. 广西嵌入式技术与智能系统重点实验室，广西 桂林 541006;3. 桂林理工大学信息科学与工程学院，广西 桂林 541006)

0 引 言

在当今世界，石油、煤炭、天然气等不可再生能源，随着人口的增长，生产力的发展，消耗越来越多越来越快，碳中和的呼声越来越受关注。现今世界，最让人关心的问题就是能源问题，它已经成为了当今世界政治、经济、科技、军事等人类领域的重点问题。传统能源很可能会出现能源危机，即能源的供应不足、储量锐减，如果没有可代替的其他能源，那么只能有意地控制能源的开采和使用，进而加剧市场供应的不足，能源价格上涨，影响到社会运行。风力发电作为清洁能源，具有显著的环保意义，社会发展价值和巨大的经济效益，因此风电的发明有着划时代的意义。地球上的风能蕴量巨大，约为2.74×109MW，其中可利用的风能大概为2×107MW，仅仅只占总量的百分之一左右，全球可以利用的水能能源总量约为可用风能的十分之一，也大于固体燃料和液体燃料能量的总和[1]。总量虽大，但是真正被利用开发得不多，想要安全、可靠、有计划地使用好这些没有被利用的风能就需要更好地研究风力机的各种受力特性以及尾流等。任年鑫等[2]采用RAN k -ω湍流模型，三维数值模型的有效性得到了很好地验证。但是缺少对风力机变桨角度改变的研究。张亚光[3]使用AL-LES方法预测风力机尾流非定常流动特性的可行性，但是大涡模拟不是很适用于风力机的计算。刘海锋等[4]运用Fluent软件，在11.4 m/s风速下，风轮转速为12.1 r/min，变桨角为0°时，使用流固耦合理论和滑移网格技术得出尾流特性的系统分析。是在特定的情况下得到的结果，不具有一般性。杨从新等[5]采用大涡模拟与致动线模型相结合的数值方法，模拟了多台NREL 5 MW风力机的风电场情况。但是在实际的应用中，大型风力机一般是独立存在的。刘智益[6]等提出串联多风力机以及错列双风力机的三维尾流模型，并通过修正尾流区域的膨胀系数 k来验证模型正确性。前人一般都是研究在来流风和风力机之间没有变桨角的情况下的尾流等气动特性，几乎没有研究过单台大型风力机变桨情况下的尾流特性，因此本文对此问题进行了探讨。

在实际中，来流风一般不会恰好与风轮所在的截面成90°流入，本文首先用ANSYS-Fluent稳态模拟计算功率和实验数据对比，对比误差均在20%以内，说明了本文的准确性。然后采用ANSYS-Fluent的瞬态计算，得出尾流的压力和速度云图，最后得出结论，分析原因。

1 理 论

1.1 CFD技术

CFD的全称是Computational Fluid Dynamics，中文称之为计算流体力学。CFD技术是基于流体力学和数值数学的理论基础，再运用现代的计算机技术将其计算结果可视化的产物。经过前人不断地完善计算模型和多次实验数据对比分析，CFD技术可作为准确的重要数据来源，因具有成本低廉，可视化程度高，以及准确度高等优点被人所喜爱。学者们先后提出一系列的工程尾流模型，但是最适合本文的计算模型还是 k -ω SST模型。目前主流的通用求解器有 CFX、Star - CD、Fluent 、star - ccm +AUTODESK CFD flow-3D、Phoenics、Comsol，其中CFX、Star - CD、Fluent等[7]，能够解决各类流体问题。本文本着节约成本和同时兼顾准确性的原则选用了Fluent软件对大型风力机进行数值模拟。

1.2 k -ω SST湍流模型

由Menter[8]提出剪切应力输送(Shear-Stress Transport，简称SST)，该模型在近壁面自由流中具有较好的计算精度，因为其综合了近壁面区计算的优点和在远场计算的优点，在湍流粘性系数的定义中添加了湍流剪切应力的输送过程，并且增加了交叉扩散项。在风力机的计算中风力机的翼型属于低负荷的模型，仅仅在风力机的轮毂位置处属于钝体绕流，因此适用于本文的风力机数值模拟计算。该湍流模型的方程为：

2 实例验证

2.1 建立模型

本文选用位于美国科罗拉多州大丹佛地区的可再生能源实验室（NREL）5 MW大型水平轴三叶片海上风力机作为本文的计算实例，风力机风轮叶片长度为61.5 m，其轮毂直径为3 m，风力机的直径为126 m。该5 MW风力机使用变速变桨的功率控制方式，相对于功率较小风力机，是具有代表性和实际应用意义的大型海上风力机，很有研究的意义。5 MW风力机叶片各个截面几何参数[9]见表1。

表1 NREL 5MW 海上大型风力机叶片几何参数

风力机叶片的建模，首先根据表1在profili软件中获取到原始的翼型坐标，然后向左平移后坐标变换，根据安装扭角在Excel中计算得出风力机叶片的实际坐标点，为了使叶片可以更好地连成一个整体，将第一个截面为一个圆形的截面形状，最后把所有的坐标导入Solidworks中，最终建成的叶片的三维效果图如图1所示。

图1 5 MW风力机三维单叶片图

在单个叶片的基础上画出一个3 m直径的轮毂之后，再使用阵列功能得到三叶片的整体风力机如图2所示。

图2 三维叶轮图

2.2 确定计算域和划分网格

为了使得风力机的尾流流域得到充分精确计算，内流域为设置为一个圆盘形状的旋转域，尺寸为直径140 m，前后内流域壁面到风轮中心基准面距离为10 m的圆柱，将外流域设置为一个5D×5D×12.5D（D为风力机直径126 m，下文中均以D表示风力机直径）的长方体计算域，其中风力机所在基准截面到进口面的距离为，出口面距离风轮基准面，风轮中心到上下左右壁面的距离为。由于直接计算原尺寸的风力机模型，风场过大，建模困难，因此本文采用缩小为原来的模型来计算。内外流域图如图3所示。

图3 内外流域图

使用 A NSYS-ICEM划分网格，风力机计算域整体采用混合结构网格，其中内流域采用非结构化网格，外流域为结构化网格。网格的划分是决定计算准确的最重要的因素，为使得计算更加精准，需要在外流域和风力机边界层加密网格。最终的网格数目为三百多万个如图4所示。

图4 流域混合结构化网格划分

2.3 功率验证

将画好的网格导入到ANSYS-FLUENT中运算，计算的参数根据实验的设定输入，选择SST湍流模型，SIMPLE算法，稳态模拟计算。

查阅资料[9]，可以获取在不同的来流风速下，风力机风轮的转速、变桨角度和功率，选取了4个数据进行比对。风速为11.4 m/s的时候为额定风速，对应的额定功率为5 MW，具体的参考数据和本文模拟数据对比如下表2。由于尺寸的缩放为原来的1/10，计算所得的扭矩也随之变为了原来的1 / 1000[10]。

表2 功率对照表

本文需要的功率计算换算公式如下：

在Fluent中只能计算出力矩的大小，因此本文采取以下的方法[11]将力矩转化为风力机的输出功率：

大量的，而且具有一定重复性的精准良好数值计算结果均证明FLUENT计算的正确性，数值模拟计算数值与实验数据的误差均在20%以内，并且越接近额定风速误差越小，整体都是功率偏大的，可能是由于网格质量导致的。表明了数值模拟计算的可靠性和正确性。因此在接下来的案例分析中其他的参数不变的前提之下，将稳态改为瞬态。

3 大型风力机变桨角度对尾流特性影响分析

3.1 计算实例

为了使得本文的研究更具备实际性，因此研究在变桨角情况下的风力机的尾流该情况下空气流经风力机之后会产生大面积的分离，这种情况下压差阻力占主导地位，S-A模型是从翼型绕流的航空领域发展而来的，通过对涡粘性系数修正的方法来适应不同区域的流动，虽然计算的速度较快，但是不适合用的分离流动的具有变桨角的风力机计算中。同样的用大涡模拟的计算结果不如RANS精准，因此选用基于SST的k-ω模型进行计算分析。

变桨角如图5所示，定义为来流风与风力机基准面的法线夹角。

图5 变桨角示意图

3.2 计算参数设置

用缩比尺寸模型可以解决建模过程中计算域太大的问题，用缩比模型CFD计算的尾流与实际的尺寸实验尾流数据相差不大[12]。

计算使用瞬态计算，重新划分具有变桨角的风力机网格，网格数目大概四百多万，变桨角环境下的内外流域网格划分图如图6所示。计算模型选择：SST，解决方案中的方法选用SIMPLE，梯度为Least Squares Cell Based；梯度、压力、动量、湍流动能、比耗散率以及瞬态离散方案均选用二阶计算。选择混合初始化，在运行计算中，计算类型为Fixed，时间步长为：0.01 s，最大迭代步数为20。

图6 21.18°变桨角时内外流网格划分图

3.3 大型风力机变桨角对尾流特性影响分析

3.3.1 速度压力速度云图整体分析

根据NREL实验室的风力机实际工作数据，本文选取三个工况分别对应的来流风速为23 m/s,18 m/s, 13 m/s对应的变桨角度分别为21.18°，14.92°，6.60°进行尾流分析，具体如表3所示。

表3 来流风速和变桨角度对照表

如图7～图9所示可以得出：随着变桨角的增大和来流风速的增加，旋转域前后的风速分层变得明显，且风力机的尾流影响的距离越远；其中工况一的尾流可以到达距离风力机旋转域中心的10D之外，三个工况共同的特点就是尾流均在4D左右就开始明显变弱，并且在旋转域周围的风速变化均较剧烈。最后，三个工况在旋转域的前面2D和旋转域的后面3D之间的尾流较为对称。可能是进口处的空气进入后遇到具有粘性的风力机叶片之后从一开始的无旋运动转变为有旋运动导致的流体微团产生了剪切变形和转动。

图7 工况一速度云图（23 m/s, 21.18°）

图8 工况二速度云图（18 m/s, 14.92°）

图9 工况三速度云图（13 m/s, 6.60°）

3.3.2 工况一旋转域速度压力云图分析

工况一的尾流比较明显，因此对其进一步分析速度和压力云图，其来流风为23 m/s，变桨角度为21.18°。

如图10～图13所示，在速度云图中可以获知旋转域中的风力机叶片处的空气流速明显低于其他处，有着较为清晰的叶片轮廓，背面叶片轮廓比正面更加明显，可能是由于叶片背风面边界层的厚度较大，使得粘性区域也增大；由于轮毂处属于钝体绕流所以背面中心点处风速较低；变桨风使得旋转域的风速一边高一边低，可能是由于倾斜的叶轮对空气具有一定的粘性切割作用，将来流风大体分成两块风速区，最终呈现出图中效果。

图10 旋转域正面速度云图

图11 旋转域背面速度云图

图12 旋转域正面压力云图

图13 旋转域背面压力云图

在压力云图中风力机的轮廓已经几乎看不到了，在风力机周围气压比较高，没有风力机叶片的形状了，压力图较为发散，旋转域的周围存在低压漩涡，正面的压力明显比背面的要高。可能是由于迎风面阻力为摩擦阻力和压差阻力之和，从而导致的旋转域正面压力较大。在背风面可能是由于边界层是分离的，使得粘性剪切力很小，从而导致旋转域背面压力较小。

3.3.3 截面面云图分析

下列各图风力机所在基准面后的速度云图，在6D距离之后的速度均为23 m/s，并无分析的价值，因此不再这里概述。

如图14～图19所示，由上述的速度云图可以得知，随着到风力机的距离增加，空气的流速越来越接近来流风速，受到风力机的影响越来越弱。在近尾流区（1D～2D）中风速受到风力机的影响风速十分紊乱，都是大致为一个矩形的高风速包围着不规则的低风速区，在3D～4D距离处，风力机的正后方还能呈现出一定区域的高速区，在左侧方向都有一定的高速区域。在远尾流区（5D～6D）处风速的差别变得非常小，6D之后的风速已经和来流风速一样，这也进一步证明了计算域足够大，计算结果具有较高的准确性。由于受到变桨风的作用，左边的风速均比右边要高，并且越靠近风力机越明显。风力机后面的空气一开始十分紊乱，当4D距离处达到临界点，空气的流速开始越来越趋近于来流风速。

图14 1D速度云图

图15 2D速度云图

图16 3D速度云图

图17 4D速度云图

图18 5D速度云图

图19 6D速度云图

3.3.4 截面压力云图分析

如图20～图23可得，在1D～2D距离时，气压比较大，受到风力机尾流影响明显，1D处大致为一较大的圆形的低压漩涡，由于变桨风的作用，左边伴随有两个较小的高压区域；2D处气旋有一分为二的趋势，存在着高压漩涡和低压漩涡，可能是轮毂处的钝体绕流引起的高压漩涡，低压漩涡可能是由于不平衡的压差力造成的；在3D处出现转折点，此处的风压的差值开始变小，层次感开始清晰，风力机的尾流影响在此处开始明显减弱。

图20 1D压力云图

图21 2D压力云图

图22 3D压力云图

图23 9D压力云图

4 结束语

1） 随着变桨角和来流风速的增加，风力机的尾流影响的距离越远，尾流均在4D左右就开始明显变弱，并且在风力机旋转域周围的风速变化均较剧烈和对称。建议在大型风力机工作时设定的安全距离为4D以上，并且在风速较大时风力机自动停止工作。

2） 在旋转域中风力机叶片背风面的层流边界层的厚度较大，从而使得粘性区域也较大，这可能使得速度云图中有着较为清晰的叶轮轮廓；但是压力云图中压力差较小形成不出叶轮轮廓，并且正面的空气流速和空气压 力均大于背面。建议在工程中加强风力机迎风处的材料强度。

3） 受到变桨的作用，外流域尾流中左侧（-Z方向）的空气流速和空气压力较右侧（+Z方向）要高一些，并且距离风力机越远越趋近来流风风速同时空气流速差也越来越小。