非均匀风影响下风力机三维气动粘性流场的数值模拟

周 胡，赵文超，万德成

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

在可再生能源领域中，风能的开发和利用技术是最为成熟也是最具大规模商业开发前景的，受到了世界各国的广泛关注。特别是最近几年海上风能发电的兴起，再一次掀起了研究和利用风能的热潮。风力机将风能转化为机械能，风资源的品质对风力机的功率、寿命和运行等起着重要影响。风力机的来流风由于受地球表面的影响，最显著的特征是在时间或空间上的分布不均。这种非均匀的来流风可以用风速廓线的形式予以量化表示。由于风速廓线引起的风切变效应的存在，叶片在旋转过程中将经历风速的周期性变化，导致风轮受到风力载荷的不均，这种不均会对风力机的使用寿命和运行安全产生深远影响。同时，在未来海上风电发展过程中，为了降低成本，叶轮直径会随着额定功率的增加而增加，风切变对风力机功率的影响将更为重要，必须在设计时引起注意。

随着计算机性能和计算方法的飞速提高和发展，使用计算流体力学理论数值模拟风力机空气动力学性能的方法日益受到了重视，数值的方法不仅能够准确预测风力机性能，还能够清楚地观察流场细节，深化对绕流场流动问题本质的理解。使用基于雷诺时均Navier－Stokes 方程求解风力机绕流场的研究已有很多［1－4］，但是使用数值方法模拟非均匀风或剪切风作用下风力机空气动力性能的工作还不是很多［5－6］，这也正是本文的立足点。

Nilay Sezer－Uzol 等［5］使用基于自由涡流理论的势流方法研究了稳态和瞬态风剪切对水平轴风力机的性能特别是涡结构的影响，文中着重研究了三种典型工况，分别是均匀风、指数率的风速廓线及瞬态极限剪切风，研究表明风速廓线的存在对风力机的尾涡结构和叶片表面的压力将产生很大非对称影响。廖明夫等［7］基于对数律分布的风切变模型从理论上对非均匀风对风力机功率影响进行了详细的研究和分析，得出了非均匀风引起的风力机设计功率的损失与风力机叶轮直径的大小相关等有用结论。刘磊等［6］依据指数率的风速分布曲线，使用了CFD 方法对非均匀风作用下风力机三维非定常气动特性进行了研究，同时还研究了不同风切变指数对风力机载荷波动振幅的影响。文中基于OpenFOAM 开源平台，选用任意网格界面元法瞬态求解器进行非均匀风作用下风力机空气动力性能预报方法的研究。

1 OpenFOAM 简介

OpenFOAM 全称为Open Field Operation and Manipulation，即场的操作和处理的开源计算平台。它是一个基于有限体积法，可用于对连续介质力学问题进行数值计算的面向对象的C+ +库，同时还提供了许多预编译好的求解器、辅助工具和模型库等。由于使用了许多C+ +语言的高级工具，例如模板类、操作符重载、多态等［8］，OpenFOAM 具有强大的可定制性和可拓展性，使用者和开发者可以根据所求问题的特殊性编写自己的求解器，重点关注求解的流程，而不需要关注离散和求解的最底层知识，这是许多商业软件无法比拟的。将OpenFOAM 作为底层库类来构建自己的求解器是许多研究者选择OpenFOAM 作为求解器开发平台的重要原因。

上海交通大学万德成课题组已经在OpenFOAM 的开发和使用上做了大量丰富工作。查晶晶、曹洪建等［9－10］基于OpenFOAM 求解器interDyMFoam，开发实现了数值粘性水池造波和阻尼消波，对圆柱波浪爬高等进行了广泛研究。沈志荣等［11］将六自由度运动模块植入OpenFOAM，开发了naoe－FOAM－SJTU 求解器，实现了对船舶在波浪上运动的数值模拟。王强、周胡等［12－14］基于MRFSimpleFoam 和pimpleDyMFoam求解器对风力机风轮的气动力性能、风力机叶片与支撑塔架相互作用的耦合流场进行了数值模拟。

2 计算模型与网格生成

Phase VI 实验风力机两叶片模型是美国国家新能源实验室于2000年在美国国家航空航天局Ames 120×80 ft 风洞中进行系列实验的模型。由于该模型具有详细的实验数据［15］，所以模型常被用来进行风力机叶片空气动力特性数值模拟结果的验证。通过求解RANS 方程结合k－ω SST 模型对实验风力机模型外流场非定常流动进行数值计算。计算模型为上风向风力机，从x 轴正向看，风轮逆时针转动，转速为72 转/分钟。风力机直径10.058 m，轮毂高度为12.192 m，叶片翼型为NREL S809，详细模型数据可以参考文献［16］。

网格的生成流程大致可以分为两过程，首先借助ICEM CFD 网格软件划分背景网格，再使用OpenFOAM自带的SnappyHexMesh 工具通过调整system 文件夹下的snappyHexMeshDict 字典参数自动生成最终网格。最终划分的网格情况如图1 所示，风力机距入口5 m，距出口20 m。在来流风向，靠近风力机前后0.8 m 范围内进行了加密，同时为了捕捉尾涡，风力机后方5 m 内也进行了加密。最终网格量大概为80 万左右。滑移界面(AMI)网格如图2 所示，交接面的具体生成过程可以参考下文任意网格界面元法的介绍。

计算工况中的空气密度、运动粘度、转速和桨矩角参考Sequence S 系列实验，选取轮毂处5、10、15 和25 m/s四个风速进行数值计算，四个风速下的初始参数设置参如表1 所示。

图1 整体网格情况Fig.1 Global mesh

图2 AMI 交界面网格Fig.2 Mesh of AMI Interface

表1 实验各工况参数Tab.1 Experiment parameter of different computation cases

3 数值模拟

3.1 控制方程

在惯性参考系中，风力机风轮以设定的角速度绕固定轴转动，流场为非定常。由于风力机外流场的风速较低，绕流场可以看成不可压。其时均控制方程可以表述为:

3.2 湍流模型

为了计算雷诺应力项，需要引入k－ω SST 模型［17］使控制方程封闭。其中k 和ω 分别代表湍动能和涡量脉动强度。两者的输运方程为:

3.3 离散方法

OpenFOAM 软件包基于有限体积法，采用空间网格的形式将计算区域划分为若干的控制体，在每个控制体上分别求解连续性方程、动量方程和能量方程。

在控制方程各项离散中，时间项采用隐式欧拉方法，对流项采用一阶高斯迎风格式，扩散项采用修正的高斯线性方法。控制方程中的速度和压力的解耦使用PIMPLE 算法(PISO(pressure－implicit split－operator)和SIMPLE(semi－implicit method for pressure－linked equations)混合算法)，PIMPLE 算法的主要结构是从PISO 中继承的，主要的区别表现在对PISO 算法的每个时间步再用SIMPLE 算法求解以得到更稳定的解。速度与压力解耦算法的深入解析可以参考文献［18］。

3.4 边界条件

3.4.1 入口边界条件

风速廓线能确定风速沿高度的变化规律，是确定给定高度处风力机输出功率的前提。风速沿高度的变化是风的重要特性之一。风速廓线一般有两种描述方式:一种是按边界层理论得到的对数风速廓线;一种是按实测结果得到的指数风速廓线。文中采用指数风速廓线来描述入口的速度分布。

风速廓线的指数分布规律可以近似表示为［19］:

图3 轮毂处两种风速下的入口风速廓线Fig.3 Wind profile of 5 m/s and 15 m/s at hub location

在OpenFOAM 中，没有现成的指数分布速度边界条件，需要自己编写。一般有两种边界条件编写的方法:其一是将边界条件加到某一个特定求解器，在使用求解器的时候就可以使用编译好的新的边界条件;其二是将需要添加的边界条件编译成一个动态链接库，需要使用该边界条件的时候，只需要在OpenFOAM的system 文件夹里的controlDict 中添加相关语句即可。文中将使用一个名为windprofileVelocity 类型的边界条件，通过编写动态链接库的方法添加。

在windprofileVelocity 边界条件编译好后，边界条件的使用就像使用其它已经编译好的边界条件一样。0文件夹下的U 文件里新的边界设置如下:

其中，inlet 是入口边界的名字，n 是来流风速的方向，meanValue 是轮毂处的风速，value 为入口速度初始值。最终计算时入口的速度分布如图4 所示，与预定义的速度廓线一致，说明入口边界条件编写正确。

图4 轮毂处5 m/s 和15 m/s 速度计算模型入口速度分布Fig.4 Inlet velocity distribution at speeds of 5 m/s and 15 m/s

入口处的压力设为零梯度，湍动能k 和涡量脉动强度ω 设为固定值，固定值的取法可根据经验公式求得，涡粘度也通过计算求得。

3.4.2 出口边界条件出口处的速度边界条件为零梯度，压力为固定值0，湍动能k 和涡量脉动强度ω 设为固定值，湍动能k 和涡量脉动强度ω 设为零梯度，涡粘度通过计算求得。

3.4.3 壁面条件

计算域外围的速度设为固定值0;压力设为零梯度;湍动能k、涡量脉动强度ω 和涡粘度采用壁面函数。

3.4.4 风力机模型边界条件

风力机的塔架和不动的轮毂的速度设为固定值0，叶片和转动的轮毂设为movingWallVelocity;整个风力机的压力边界条件设为零梯度;湍动能k、涡量脉动强度ω 和涡粘度采用壁面函数。

3.5 任意网格界面元法

在OpenFOAM 新的版本中，采用任意网格界面元法来处理旋转的交界面。任意网格界面元法的本质是一种滑移网格技术。任意网格界面元法和之前版本中的通用网格界面法原理基本类似，但是也存在一些细微区别。任意网格界面元法主要是通过插值实现动静区域流场参数和信息的交换。为描述方便，界面两边运动和静止的面分别定义为主面(master)和从面(slave)。其求解过程主要步骤如下:

图5 主从面的关系示意Fig.5 Sketch map of relation between rotor and stator

3.5.1 创建AMI 交界面

在OpenFOAM 中，创建AMI 交界面可以通过网格的拓扑运算完成，例如使用topoSet 工具来选取相应的面域。之后需要把该交界面复制一份，使得一个面属于主面，一个属于从面，如图5 所示，图中rotor side 代表主面，stator side 代表从面。这一步操作可以使用OpenFOAM 中的createBaffles 和splitOrMergeBaffles两个工具完成。

3.5.2 寻找相邻面

当AMI 交界面生成后，就可以开始计算了。计算的第一步便是寻找相邻面。寻找相邻面的过程是已知主面寻找相邻从面的过程。

如图6 所示，运动部分的主面O 与静止部分的从面a、b、c 三个面相邻，故可得这三个面为主面O 的相邻面。具体搜寻相邻面的算法在此不作介绍。

3.5.3 计算权重

权重的计算是相互的，即需要计算从主面到从面的权重，也需要计算从从面到主面的权重。权重由相重叠区域占该面的面积来表示:

主面到从面的权重:

图6 相互滑移网格交界面的处理［21］Fig.6 Interaction of different sliding meshes

从面到主面的权重:

式中:i 为某一个主面的第i 个相邻从面;j 为某一个从面的第j 个相邻主面;为第j 个主面和从面重叠区域面积大小，第i 个从面和主面重叠面积大小;为某一从面第j 个主面的面积和某一主面第i 个从面的面积。

例如对于图6，想要计算从面b 到主面O 的权重，那么可以通过以上公式计算得到权重。由此可以得到与主面O 相邻的所有从面到主面的权重。

3.5.4 差值

有了权重，那么就可以通过带权重差值的方式把速度场、压力场等在滑移面两侧实现数据交互。

由主面中变量差值得出从面中变量值:

ΦS= ∑WMj_to_SΦMj

由从面中变量差值得出主面中变量值:

ΦM= ∑WSi_to_MΦSi

式中:ΦMj，ΦSi为某一从面第j 个主面中Φ 值和某一主面第i 个从面中的Φ 值。

通过这种带权重的差值，就实现了网格相对运动时主面和从面之间数据的交互。

4 计算结果与分析

4.1 推力的时历曲线

风力机在非均匀风的作用下，叶片在旋转过程处于不同的位置将经历风速的周期性变化，这种变化会导致叶片在整个的受风面上风载荷的不均匀，会对风力机的使用年限和运行安全产生很大影响。图7(a)～(d)分别表示一圈内不同风速均匀风和剪切风作用下推力随方位角的变化对比图。为方便比较，此处的推力为单个叶片所受的推力，非均匀风下的叶片经历了从高风速区到低风速区再回到高风速区的过程。从图中可以明显观察到推力也经历了一个从高到低再到高的过程，并且推力在风速最低处(对应方位角180°)出现了最低值，这与理论分析完全一致。经过计算发现一圈内的叶片所受推力的平均值变化不超过5%，但是由于非均匀风引起的脉动值却差距很大。例如10 m/s 风速下，均匀风下推力的脉动量占平均推力的7%左右，但是非均匀风作用下推力的脉动量却占了平均推力的14.3%，几乎是均匀风的两倍。对于25 m/s 风速，均匀风和非均匀风对推力脉动量的影响差异将更大，均匀风影响下推力的脉动量约是平均推力的4.2%，非均匀风影响下推力的脉动量约是平均推力的11.1%，约是均匀风的三倍。

从上面的分析可以看出，非均匀风影响下风力机叶片所受推力的脉动值会比均匀风影响下的脉动值大很多，这对于结构的疲劳将产生非常不利的影响。所以在风力机的疲劳设计阶段，需要根据具体风场的风速廓线情况，对非均匀风所引起的疲劳载荷予以考虑。

图7 四种不同风速下剪切风与均匀风推力时历曲线Fig.7 Comparison of time histories of thrust for different uniform and nonuniform winds

4.2 涡的分布

尾涡区的紊流结构以及叶尖涡的捕获是风力机气动力学研究所感兴趣的，使用速度梯度张量的二阶不变量Q 对尾流场的涡结构进行可视化处理，Q 的定义为:

其中，Ωij代表漩涡的强度，Sij表示剪切应变率，分别是速度梯度反对称和对称分离，可以理解为剪切应变率和涡量间的局部平衡量［22］。

2.6 s 轮毂处风速为5、10、15 和25 m/s 时，非均匀风影响下风力机尾涡情况如图8 所示。图中Q 使用上面定义的速度梯度张量的二阶不变量的等值线表达，并且使用速度染色。风轮上部区域涡的颜色要比下部涡的颜色深，这与非均匀风入口设置是一致的，也验证了文中计算模型的正确性。从图中可以清晰看到叶尖涡和叶片根部过渡区的涡的脱落，同时在高风速下还能看到叶片处涡的分离。此外，塔架和叶尖涡相互作用也可以明显看出，而塔架下半部分没有明显的涡生成，可能是下半部分网格变稀疏的缘故造成的。

图8 不同风速非均匀风影响下风力机尾涡情况Fig.8 Vortex structure for uniform and nonuniform winds at different speeds

4.3 尾流场分析

尾迹区分析是风力机空气动力学分析的重要部分。图9(a)给出了轮毂处5 m/s 风速非均匀风作用下不同顺风向位置中剖面处的风速曲线，分别截取了距叶片4、8、12 和16 m 处与塔架平行的直线上风速随高度的变化关系。从图中可以清晰的观察到，距离叶片越近尾流场越紊乱，在距离叶片4 m 纵坐标为0 时，风力机轮毂的地方风速有急剧的下降，这主要是由于轮毂的遮蔽作用，同时也注意到在距离叶片16 m 处风速基本恢复且接近于入口处的非均匀分布。图9(b)给出了轮毂处10 m/s 风速非均匀风下不同顺风向位置中剖面处的风速曲线。与图9(a)一样，在距离轮毂较近的尾流区风速骤降，随着据叶片距离的增大，风速逐渐恢复但是10 m/s 的恢复比5 m/s 风速下要慢很多，可以推出高风速下尾流区的影响更长，需要更长的距离进行恢复。在进行风场风力机布置时，一定要着重考虑风力机之间的间距以减小上一风力机尾流场对下一风力机功率的影响。

图9 轮毂处5 m/s 和10 m/s 横风向尾迹区风速曲线Fig.9 Cross－wake speed distribution at 5 m/s and 10 m/s

为了进一步探讨非均匀风影响下塔架与叶片的相互作用，以10 m/s 风速为例，取2.5 s 时风力机下方叶片r/R=0.15 和r/R=0.95 位置的水平截面进行研究，如图10 所示。图中背景以速度值染色，涡由Q 等值面表达并由速度染色，白色圆圈代表塔架，白色翼型代表的是叶片。首先观察速度场，可以看到塔筒前后的速度值有明显下降，而塔筒后的速度场有明显的波浪形变化，这是典型的圆柱绕流的速度场。从图中还可以看出叶片尖部(z= －4.78 m)的攻角要比叶片根部(z= －1.51 m)的攻角小，叶片尖部泄涡更为明显。叶片尖部的泄涡与塔架的尾涡相互作用可以从图10(b)中看出，而叶根处塔架与叶片尾涡相互作用的过程却不是特别明显。在轮毂10 m/s 非均匀风速下叶片分离涡中的叶尖涡与塔架尾涡相互作用的过程更明显，影响更大。

图10 10 m/s 非均匀风速下风力机下叶片不同截面塔架与叶片相互作用尾涡情况Fig.10 Blade tip vortices and tower vortices interaction for uniform and nonuniform winds at two different sections(cut plane at z= －1.51 m and －4.78 m)

4.4 叶片截面的压力系数

通过将均匀风和非均匀风作用下叶片各个截面压力分布情况进行对比可以进一步探讨非均匀风对风力机空气动力性能的影响的细节。选取了2.5 s 时(即风力机正好旋转3 圈时)对四个不同截面(r/R =0.466，0.633，0.8，0.95)的压力情况做对比研究。文中用到的压力系数定义如下:

其中，P0为叶片表面附近的压力值;P∞为无穷远处的压力值，文中取0;U 代表风速;ω 代表风轮角速度;r 代表截面距转动中心的距离。

从图11 中能够清晰的看出在5 m/s 低风速时，非均匀风的影响很小，叶片截面的压力系数只有细微差别，同时随着截面从根部向叶尖处发展，这种影响会越来越小，这从10 m/s 风速下的不同截面的压力系数分布能更清晰的看出。

图11 5 m/s 和10 m/s 均匀风和非均匀风叶片各截面压力系数分布曲线Fig.11 Pressure coefficient distribution for uniform and nonuniform winds at four sections at 5，10 m/s

5 结 语

在考虑非均匀风影响的情况下，对风力机整体结构的三维非定常气动问题进行了数值模拟，得出以下结论:

1)在非均匀风影响下，风力机所受推力的脉动值会比均匀风影响下的脉动值大很多，这对结构的疲劳将产生非常不利的影响。所以在风力机的疲劳设计阶段，需要根据具体风场的风速廓线情况，对非均匀风所引起的疲劳载荷予以考虑。

2)非均匀风的存在会引起尾涡结构的非对称，同时对尾流场的分析中可以发现尾流区风速分布非常紊乱，且高风速下风速的恢复更慢。

3)非均匀风影响下的叶片截面压力系数在低风速时与均匀风相比变化不大，但是高风速时有很大差距，且随着截面从叶片根部到叶片尖部发展，这种影响会逐渐变小。

［1］SøRENSEN N N，MICHELSEN J，SCHRECK S.Navier－Stokes predictions of the NREL phase VI rotor in the NASA Ames 80 ft× 120 ft wind tunnel［J］.Wind Energy，2002，5(2－3):151－169.

［2］李宇红，张庆麟.风力机叶片三维流动特性与气动性能的数值分析［J］.太阳能学报，2008，29(9):1172－1176.(LI Yuhong，ZHANG Qinglin.Numerical simulation of flow field and aerodynamic performance of a wind turbine blade［J］.Acta Energiae Solaris Sinica，2008，29(9):1172－1176.(in Chinese))

［3］张果宇，蒋劲，刘长陆.风力发电机整机气动性能数值模拟计算与仿真研究［J］.华东电力，2009，37(3):449－452.(ZHANG Guoyu，JIANG Jin，LIU Changlu.Numerical simulation of aerodynamic performance for wind turbines［J］.East China Electric Power，2009，37(3):449－452.(in Chinese))

［4］任年鑫，欧进萍.大型海上风力机尾迹区域风场分析［J］.计算力学学报，2012，29(3):327－331.(REN Nianxin，OU Jinping.Numerical analysis for the wake zone of large offshore wind turbine［J］.Chinese Journal of Computational Mechanics，2012，29(3):327－331.(in Chinese))

［5］SEZER－UZOL N，UZOL O.Effect of steady and transient wind shear on the wake structure and performance of a horizontal axis wind turbine rotor［J］.Wind Energy，2013，16(1):1－17.

［6］刘磊，石可重，杨科，等.风切变对风力机气动载荷的影响［J］.工程热物理学报，2010，31(10):1667－1670.(LIU Lei，SHI Kezhong，YANG Ke，et al.Effect of wind shear on the aerodynamic load of wind turbine［J］.Journal of Engineering Thermophysics，2010，31(10):1667－1670.(in Chinese))

［7］廖明夫，徐可，吴斌，等.风切变对风力机功率的影响［J］.沈阳工业大学学报，2008，30(2):163－167.(LIAO Mingfu，XU Ke，WU Bin，et al.Effect of wind shear on wind turbine power［J］.Journal of Shengyang University of Technology，2008，30(2):163－167.(in Chinese))

［8］JASAK H，JEMCOV A，TUKOVIC Z.Openfoam:A c + + library for complex physics simulations［C］//Proceedings of International Workshop on Coupled Methods in Numerical Dynamics.2007.

［9］查晶晶，万德成.用OpenFOAM 实现数值水池造波和消波［J］.海洋工程，2011，29(3):1－12.(CHA Jingjing ，WAN Decheng.Numerical wave generation and absorption based on OpenFOAM［J］.The Ocean Engineering，2011，29(3):1－12.(in Chinese))

［10］CAO H J，WAN D C.Development of multidirectional nonlinear numerical wave tank by naoe－FOAM－SJTU solver［J］.International Journal of Ocean System Engineering，2014，4(1):52－59.

［11］SHEN Z R，WAN D C.RANS computations of added resistance and motions of ship in head waves［J］.International Journal of Offshore and Polar Engineering，2013，23(4):263－271.

［12］周胡，王强，万德成.风机叶片三维绕流场数值模拟［C］//第十一届全国水动力学学术会议暨第二十四届全国水动力学研讨会并周培源教授诞辰110 周年纪念大会.2012:627－636.(ZHOU Hu，WANG Qiang，WAN Decheng.Numerical Simulation of the 3D viscous flow field over wind turbine Blade［C］//Proceedings of 11th National Conference of Hydrodynamics.2012:627－636.(in Chinese))

［13］WANG Q，ZHOU H，WAN D C.Numerical simulation of wind turbine blade－tower interaction［J］.Journal of Marine Science and Application，2012，11(3):321－327.

［14］周胡，万德成.下风向风力机塔影效应的非定常数值模拟［C］//第二十五届全国水动力学研讨会暨第十二届全国水动力学学术会议.2013:277－283.(ZHOU Hu，WAN Decheng.Unsteady numerical simulation of tower shadow of downwind wind turbine［C］//Proceedings of 12th National Conference of Hydrodynamics.2013:277－283.(in Chinese))

［15］FINGERSH L J，SIMMS D，HAND M，et al.Wind tunnel testing of NREL’s unsteady aerodynamics experiment［J］.AIAA paper，2001:35－46.

［16］HAND M M，SIMMS D，FINGERSH L，et al.Unsteady aerodynamics experiment phase v:test configuration and available data campaigns［R］.National Renewable Energy Laboratory，2001.

［17］MENTER F R.Two －equation eddy －viscosity turbulence models for engineering applications［J］.AIAA Journal，2012，32(8):76－82.

［18］JASAK H.Error analysis and estimation for the finite volume method with applications to fluid flows［D］.London:University of London，1996.

［19］贺德馨.风工程与工业空气动力学［M］.北京:国防工业出版社，2006.(HE Dexin.Wind engineering and industrial aerodynamics［M］.Beijing:National Defense Industry Press，2006.(in Chinese))

［20］GB.建筑结构荷载规范［S］.2002.(GB.Load code for the design of building structures［S］.2002.(in Chinese))

［21］JASAK H，BEAUDOIN M.OpenFOAM TURBO TOOLS:from general purpose CFD to turbomachinery simulations［C］//Proceedings of ASME－JSME－KSME Joint Fluids Engineer－ing Conference (AJK2011－FED).2011.

［22］JEONG J，HUSSAIN F.On the identification of a vortex［J］.Journal of Fluid Mechanics，1995，285(69):69－94.