实尺度大型水平轴风机气动特性数值模拟

王新茹， 陈 刚， 肖龙飞， 温斌荣， 田新亮

(上海交通大学 a.海洋工程国家重点实验室； b.高新船舶与深海开发装备协同创新中心， 上海 200240)

0 引 言

近年来，随着海上风电的快速发展，海上风能的利用已成为可再生能源发展的重要方向[1]。风机的气动力性能直接影响其发电效率和使用寿命。因此，对风机周围流场的精细化模拟和分析作为海上风机组设计和研究的基础，逐渐成为人们关注的关键技术之一[2-3]。

随着计算机性能的提高和数值计算方法的发展，人们对风机周围流场的数值分析方法进行了大量的研究和探索。

INGRAM[4]在一维动量理论的基础上，提出叶素动量理论，该理论是计算风机叶片气动力性能的经典方法。但是，叶素动量理论在不同雷诺数和攻角下对二维翼型升力系数和阻力系数的依赖使其难以精确模拟复杂工况下风机的三维旋转运动。为减小误差，人们引入叶尖损失修正模型、动态失速修正模型、轮毂损失修正模型和偏航尾迹修正模型等诸多修正方法[5]。这些方法虽基本实现了快速高效地计算风机载荷情况，但在风机周围流场分布方面却显得无能为力。为进一步精确计算风机气动力性能，并提供风机周围流场信息，KEITH等[6]和LANDAHI等[7]分别提出各种气动力学势流模型，如面源法、升力面理论等，对风机的三维旋转进行较为准确的近似模拟。这类势流方法在大型风机非定常气动分析方面取得了较好的结果[8-10]。然而，这类方法通过工程模型的方式近似模拟叶片周围的流体黏性和流动分离等关键现象，不能提供叶片周边精密细致的流场分布。

近几年，CFD方法应运而生，基于雷诺平均的Navier-Stokes方程，选取适当的湍流模型对风机叶片周围的三维流场进行数值模拟，通过计算机求解代数方程得到结果，最终还可将计算结果通过图形或图像进行可视化显示，以便更好地理解和分析问题。但是，由于风机的翼型结构及其绕流场形状十分复杂，风机的实型建模和高质量的网格划分变得困难，因此前期的CFD计算大多采用SORENSEN等[11]提出的制动线和制动面模型，将风机的实型简化成原本并不存在的线或圆盘并离散成多个单元，求解所得的单位叶片所受升力和阻力便通过这些虚拟载体代入叶素理论从而获得整个风机所受的升力和阻力。这种近似模型计算方法虽在计算精度上有所改善，但忽略了风机的结构对其气动力性能和周围流场影响，使结果与现实情况相比仍存在一定的差距。

ZHOU等[12]对风机模型进行原型建模，采用滑移网格方法考虑风机旋转，探究风机结构细节及叶片旋转对其周围流场产生的影响。TRAN等[13]、LIU 等[14]和CAI 等[15]对完全考虑结构特点的NREL 5 MW风机进行实尺度建模，通过重叠网格的方法对旋转风机进行数值模拟，探究风机的流场特征及其在风作用下的运动响应。但是，模型所采用网格数量较少，对流场细节的捕捉还不够精细，使其在揭示物理现象和物理机制时略粗糙。

本文在对NREL 5 MW参考模型进行实尺度精细化建模的基础上，采用极大规模的网格处理模型特征，还原叶片扭转角、翼型过渡和叶尖等部分结构细节，采用OpenFOAM计算软件，结合自主开发的旋转动网格处理方法，利用pimpleDyMFoam求解器对风机叶片的气动力性能和周围流场进行全方位、多角度的精细化模拟和分析。

图1 计算模型示例(非正确比例)

1 计算模型与方法

1.1 计算模型

所采用的风机模型为NREL 5 MW参考风机[16]，这是一款上风向三叶片水平轴风机。出于计算简化考虑，未考虑风机塔柱的影响。风机盘面直径D=126 m，其所处圆柱型流域如图1所示。风机旋转平面到入口的距离约2D，距出口为10D。来流风速取10 m/s，坐标轴x的正方向与风向一致，风机转速为顺时针(+x)12 r/min。大气参数分别取：密度ρ=1.29 kg/m3，运动学黏性系数ν=14.8×10-6m2/s。

1.2 网格划分

CFD的本质是将连续介质进行区域化离散，通过插值运算将物理场中的连续量近似到各离散点上，再求解各离散点上的物理量从而获得所求流场中各物理量的连续性分布。也就是说，计算是在生成的各离散网格点上进行的。因此，将一个给定形状的物体在所需研究的流场内进行合适的网格划分十分重要。

本文首先通过ANSYS中的网格划分软件ICEM将圆柱型流域内部嵌套一个O型网格用于模拟加密的叶片网格，对风机叶片尾流部分的网格进行适当加密，以期捕捉到更真实、精细的尾涡结构。而后，利用snappyHexMesh工具对网格进行捕捉、分割、移动，以及添加边界层网格。由于风机为实体，因此在流场计算中需将其内部的网格掏空，并增加相应的边界层。整个计算域的网格划分结果如图2所示，风机叶片及边界层划分如图3所示，网格总数量为1 520万个。

图2 全计算域网格划分结果

图3 风机叶片附近网格及边界层

同时，基于OpenFOAM 1.6版本，在求解器中引入课题组研发的整体旋转网格的动网格方法[17]，避免滑移网格法信息交互误差，使整个计算网格随风机叶片一起旋转，保证长时间高精度高效率的计算。

1.3 控制方程

风机外流场的马赫数一般小于0.3，故可将流体视为牛顿不可压缩流体。非定常流场的控制方程(N-S方程)[18]为

(1)

(2)

式(1)和式(2)中：p为压力；ρ为流场密度；xi、xj为笛卡尔坐标系下点的坐标(i、j=1、2、3，分别表示x、y、z方向)；ui为速度在xi方向上的分量；t为时间。

计算流体动力学开源软件OpenFOAM采用有限体积法对N-S方程进行离散化，基于张量法和面向对象技术[19]，用于解决连续介质中的流体力学问题。采用PIMPLE(PISO和SIMPLE混合算法)的离散格式求解速度与压力的耦合，用于插值、梯度、拉普拉斯和散度计算的空间格式分别为线性、高斯线性、高斯线性修正和高斯线性。对流项采用二阶迎风格式，时间项采用二阶Crank-Nicholson格式。以0.000 7 s步长在上海交通大学超级计算中心平台中以64核计算两周得到计算结果。

2 计算结果分析

2.1 风机气动力特性

风机的基本工作原理是利用空气流经风机叶片时叶片所产生的升力或阻力推动叶片转动，通过传动系统将转矩传到发电机，将风能转化为电能。因此，风机作为一种流体机械存在效率问题。通常，主要通过推力系数和风能利用系数评价风机的气动力性能[20]。

风机叶片的叶尖线速度与对应风速之比称为翼尖速比，所研究风机模型的翼尖速比为

(3)

式中：ω为风机转速；R为风机半径；U为来流风速。

推力是指风沿来流方向对风机表面的作用力，推力系数CT为

(4)

式中：T为风向方向的风载荷分量；ρ0为空气密度。

风能利用系数Cp用于表征风轮对风能的捕获能力，其表达式为

(5)

式中：M为风力扭矩。

通过对OpenFOAM中翼尖速比为8时计算得到的风机所受推力和风能利用系数与其他计算方法在相同翼尖速比时所得结果[8， 21-22]进行对比，分析风机气动力性能细节以及所采用模拟方法的实用性，如表1所示。

观察可知，风机所受推力系数CT和风能利用系数CP的数值计算结果与其他计算方法所得结果基本吻合，间接证明所采用数值模拟方法的实用性。

表1 风机翼尖速比为8时不同计算方法推力系数和风能利用系数对比

2.2 尾涡分布

当水平轴风机的叶片被限制在一定长度时，叶片旋转会导致在风机下游部分产生螺旋形涡。产生的涡主要存在于两个区域：风机的叶尖处，也称叶尖涡；风机下游旋转轴的水平方向处。根据距离风机叶片的远近，一般将风机尾涡区划分为近尾涡区(x/D≤ 1)和远尾涡区(x/D>1)，近尾涡区的各气动力学参数会对风机自身的运行成本、发电功率和使用寿命产生较大影响，而远尾涡区则会影响风电场中其他风机的气动力性能。通过涡的可视化体现，可清楚地了解风机的尾涡结构及其下游流场形态。采用Q准则对尾涡区的涡结构进行表示，Q的定义为

(6)

式中：Ω为旋涡强度；S为剪切应变率。图4中Q= 0.01，并用风速方向速度分量值染色。

图4 风机尾涡结构示例

图5给出了流场中顺风向距离风机叶片x/D=-1.0、0.5、1.0、1.5、2.0时，风速在z轴方向上的对比曲线。z/R用于表征z方向距离与风机半径的比值，同时引入α=u/U因子表征风机叶片不同截面高度处流场的速度损失情况。可以看出，流场的风速分布关于x轴基本对称。在远尾涡区，流场速度受风机轮毂和叶片阻隔的影响而有所下降，且距离越近，速度下降越明显；在近尾涡区，由于风机轮毂和叶片(尤其是叶尖处)的旋转带动作用，流场风速在z= 0 m和z>40 m处已基本恢复。

图5 不同截面流场速度轴向分布图

2.3 瞬时流线图

2.3.1 叶片表面正负压分界线

风流过旋转的风机后，其周围流场在风机的前后表面会形成正负压区域，图6所示的即为风机表面流场正负压分界线情况。图中黑色表示正压，灰色表示负压。由观察可知，风机叶片的压力面成正负压交替分布状，而在其吸力面，流场压力分布以负压为主。

图6 叶片表面正负压分界线示例

在叶片表面用三维流线表征流场的流动特性，清晰可见在分界线周围流线分支、汇聚，且在风机轮毂后方形成涡，如图7所示。

图7 叶片表面流线图

2.3.2 来流方向流线

为更形象细致地研究风机对于风机流场的干扰作用，沿来流方向分别取x/D= -1.0、0.5、1.0、1.5、2.0处的流场截面，背景以压力染色，分别观察不同截面处流线与压力特征，圆形表征风机盘面大小，y、z轴为各截面处真实坐标值，如图8所示。

图8 来流方向上不同截面的瞬时流线

由观察可得：在叶片截面处，叶片周围流场的压力增大，瞬时流线在叶片吸力面处形成旋涡；距离风机x/D= 0.5处的流场最为明显的特征是截面流线在旋转轴中心处形成旋涡，流场压力降为负值，且越靠近旋转轴中心，三叶片投影区域的反向压力最大，说明该处截面受旋转叶片的影响显著；在x/D= 1处，流场压力仍为负值，但流线与之前不同，呈类海螺状，旋涡形状与之前相比变大；近尾涡区(见图8(d))之后，压力负值逐渐减小，流线螺旋态减小，风机叶片对流场的影响减弱。

2.3.3 翼型截面流线

选取5个典型截面的压强系数，截面分布(r/R=0.30、0.47、0.63、0.80、0.95)如图9所示，分析翼型截面的流场流线分布。

图9 5 MW风机叶片5个典型截面

如图10所示，流场压力在翼型处存在明显分界，压力面压力增大，吸力面压力减少，过渡区与翼型处于同一截面，压力成阶梯状由大至小分布，流线在压力面和引力面分别形成涡，随着翼型截面的减小，产生的涡也相应减小直至消失。

图10 不同翼型截面的瞬时流线

3 结 论

基于OpenFOAM，采用整体旋转动网格处理方法对大型风机周围三维流场进行了模拟与分析。从数值模拟的流场来看，由于翼型的存在，风机周围流场结构复杂，叶尖涡与尾流区分界明显，来流方向上不同截面的流场特性随着其与风机表面距离远近的变化而变化，涡的产生对风机自身部件和风电场其他风机都会产生一定的作用。

数值计算结果的对比说明，该数值计算方法对于计算风机气动力特性和流场真实情形的模拟较为准确，计算结果吻合较好，可作为今后风机研究的方法和手段。