0 引　言

随着科学技术的快速变更，资源的合理利用受到更多的关注，不可再生能源的缺点，使可再生能源优势更突出，风能作为一种可再生能源，受到很多学者的青睐。风力机使得风能转换为电能，大功率的风力机转化效率突出，经济效益高，大功率风力机比小功率风力机优势更明显。风电场中当来流风速经过旋转的叶片，受到旋转叶片的粘性阻碍，空气在叶片附近流动，进而被吸收发生能量损耗，造成速度损失，湍流度变强，涡能量损耗，造成尾流效应。根据贝兹理论极限值，风能转化为59.3%，实际工程中，这一数值会更低，此外尾流效应带来的湍流度增强也造成下游风力机叶片疲劳，影响风力机使用寿命，因此研究风力机尾流有其重要的意义，可为风电场地资源的合理利用和提高风能利用率提供参考。风力机尾流分布存在近尾流区和远尾流区，近尾流区流动情况复杂，远尾流区是除了近尾流区外的区域，主要受到风力机尾流的相互作用以及尾流效应带来的湍动能变化等。

以往很多学者研究风力机的方法主要分为实验研究法和数值模拟法两种，实验研究是通过风洞实验，给定流场条件，研究风力机在风洞中的各种湍流特性。在早期，李万润等^[1] 就通过实地测算风场中尾流的实际速度数值，对风力机的尾流进行分析；随着计算机技术的快速发展，数值模拟方法越来越受到欢迎。在数值模拟方面，杨立国等^[2] 结合风洞试验与数值模拟，两者数据结果对比，得出数值模拟的正确性；孙芳锦等^[3] 用CFD方法对单台风力机不同变桨角度的尾流特性进行定性分析；刘智益等^[4] 通过对串列与错列风力机模型，修正尾流区域膨胀系数来验证模型准确性。王兵等^[5] 用大涡模拟比较了六种亚格子模型，分析了六种亚格子模型在湍流运动中瞬时场的压力和粘性大小的分布情况，定性给出了六种亚格子模型选取的导向；张亚光^[6] 用大涡模拟结合致动线风力机几何模型，对单台风力机的尾流特性进行分析。致动线是一种将风力机等效为沿叶片径向分布的体积力，以线表达的一种技术，致动线几何模型极大减少网格数量，节省计算资源，一般结合大涡模拟运用，但致动线几何模型会忽略叶片翼型细节，是一种简化模型，不符合实际情况，故本文采用未简化风力机几何模型，结合大涡模拟进行尾流特性研究；在尾流研究方面，杨瑞等^[7] 对水平轴风力机尾流研究揭示尾流产生平面机理及在叶片末端膨胀发展和向下游风力机扩散的机理；钟宏民^[8] 耦合大涡模拟与致动线的技术研究了风力机尾流湍流运动的变化情况；温文等^[9] 利用数值模拟对2 MW风力机串列下不同间距尾流进行分析，得到2 MW风力机不同间距下尾流的影响，间距越大，影响越小，但是目前主要趋势是针对大功率风力机进行研究，提高经济效益。

本文采用数值模拟的方法，构建风力机模型及流场空间，结合陈荣盛等^[10] 提出的相似准则，对大功率NREL-5 MW风力机模型按1∶10比例缩小，解决因流域大导致计算机资源浪费的问题。对单个风力机模型进行不同风速下风力机功率和水平推力验证，同美国可再生实验室风洞实验做对比，分析了不同风速下模型的功率和水平推力数据，验证模型的可行性。并以此为基础，在风速入口采用UDF函数输入剪切风模型，对风力机采用正对串列的形式，控制了上下游风力机的间距，采用大涡模拟湍流模型，分析不同间距下5 MW大功率风力机的尾流特性变化。

1 CFD数值模拟 1.1 风力机模型

模型采用美国可再生实验所给出的翼型数据，在Solidworks建模软件中，调整翼型坐标，对各个翼型建立一个翼型平面，组合成风力机单叶片，最终结合轮毂生成三叶片模型。基本参数见表1，叶片及风力机模型见图1和图2。5 MW风力机直径为126 m，额定风速为11.4 m/s，切入风速3 m/s，切出风速25 m/s，额定转速12.1 r/min，轮毂高度为90 m。

1.2 计算域与网格划分

计算域为长方体流域，流域入口处两侧距风力机轮毂中心为1.5D，风力机轮毂处距底部9 m（未缩放实际高度为90 m），即轮毂高度，轮毂距上部对称壁面为2D，上游风力机距入口为1.5D，上下游风力机间距设为${{xD}}$，${{x}}$取5、8、11，风力机下游距出口处为10D间距，其中D=126 m为风力机直径，具体如图3所示。

网格划分对计算结果有着重要的影响，采用 ICEM软件对流域进行网格处理，大涡模拟对网格要求很高，考虑到结构化网格精度更高更好计算，但由于风力机表面比较复杂，较难生成结构化网格，基于外网格跟内网格分块思想，分成两个块，内网格块采用非结构网格，做法为风力机表面过渡到外空间部分逐层加密，边界层网格Y+小于1，边界层网格划分5层，保证捕抓到精细的流场变化，由于外流场分布均匀，故外网格采用更高质量的结构化网格。综上所述，用结构化外网格混合非结构内网格，形成混合网格的结构形式，如图4所示。

1.3 大涡模拟湍流模型

目前，湍流数值模拟方法分为直接数值模拟法(DNS)、雷诺平均模拟法(RANS)和大涡模拟法三种方法，直接模拟法用不可压缩流动的连续性N-S方程求解湍流运动中湍流细节，雷诺平均法是用时均N-S方程求解时均运动，小的脉动忽略，得到接近平均的结果，但很难准确预测尾流内的详细湍流特性^[11-13] ，且对比LES，RANS的计算精度较低，综合考虑，本文采用大涡模拟，大涡模拟是一种介于DNS和RANS方法的一种方法。大涡模拟的整体思想是认为大涡被流场影响较大，小涡认为是各向同性，因此将湍流运动紊流涡的空间平均，通过滤波函数将大尺度的涡与小尺度的涡过滤，大涡直接用方程模拟，小涡通过亚格子直接模拟。大涡模拟的控制方程如下式所示：

式中：$ {{\boldsymbol{\tilde u}}_{i,}}{{\boldsymbol{\tilde u}}_j} $——被过滤后的速度；

$ {{{\boldsymbol{\tau}}}_{ij}} $——亚格子应力；

$ \tilde p $ ——被过滤后的压力；

$ \mu $ ——动力粘度；

$ {F_i} $ ——源项，代表风力机引起的体积外力对动量方程的影响。

其中亚格子应力$ {{{\boldsymbol{\tau}}}_{ij}} $表示为

式中：${{\boldsymbol{\mu }}_{\rm t}}$——亚格子湍流度；

$ {{\boldsymbol{\delta }}_{ij}} $——克罗内则符号；

$ {{\boldsymbol{\tilde S}}_{ij}} $——求解亚格子的应变率张量，定义如下式所示。

为实现大涡模拟，还应该具有亚格子应力模型，本文大涡模拟研究采用的是动态Smagorinsky-Lilly模型，利用Smagorinsky-Lilly模型来求解亚格子湍流粘度$ {{\boldsymbol{\mu }}_{\text{t}}} $。

其中${L_{\rm s}}$定义为：

式中：$ {\text{κ} } $——冯卡门常数；

$ \delta $——近壁面距离；

$ V $——计算单元的体积；

$ {{\text{C}}_{\text{S}}} $−Smagorinsky常数，根据Oka和Ishihara的研究，${{{\rm{C}}}_{{\rm{S}}}}$=0.032^[13] 。

1.4 模拟边界条件设置

本文流场的入口条件设置为速度入口，速度采用UDF程序输入剪来流风，产生的风速随着高度而变化，即x方向的速度随着z方向高度变化而变化，出口边界为压力出口边界，风力机计算域上部与左右部壁面采用对称边界条件，为现实叶轮旋转，叶轮旋转区域建立旋转域，旋转域表面跟外流域交界面采用interface面静动网格交界面，设定旋转域旋转，叶轮相对于旋转域静止以实现叶轮转动，故风力机叶轮采用无滑移壁面，底部也设定无滑移壁面。具体如图3所示，采用大涡模拟湍流模型，亚格子模型采用动态Smagorinsky-Lilly模型，算法采用simple算法对速度跟压力方程进行耦合计算，扩散项用中心差分格式，时间步长采用风力机旋转3°所用的时间即0.04 s，收敛精度为0.001，先进行稳态计算，为瞬态流场提供良好的初始流场信息，在进行瞬态大涡模拟计算，每种工况计算20 s。

在地表附近，由于空气会形成一系列不规则的涡流，形成湍流，对空气流动造成很大的阻碍作用，由于受到近地层的阻碍作用，风速随着高度会呈现出指数函数形式变化，即高度越高，风速越大，一定高度后风速会趋于常数，这种现象被称为风剪切，本文采用FLUENT中UDF功能，编写风速函数程序，形成垂直风剪切来流，具体如图5所示，入口处的风速随着高度变大，风速也变大，在轮毂处形成额定风速，具有明显的速度分层，相对于均匀来流，符合实际的入流条件，具体公式为：

式中：U₀——轮毂处的风速，取11.4 m/s；

Z₀——轮毂处的高度为90 m，模型缩小为1/10，取9 m；

Z——速度入口边界任意高度；

$ \alpha $——地面粗糙度系数。

A类地貌取0.12，具体为沙漠、海岛等；B类地貌取0.16，具体为田野、乡村、丛林、丘陵以及房屋比较稀疏的乡镇和城市郊区。本文针对B类粗糙度系数进行模拟，取$ \alpha $=0.16，分析尾流。

1.5 相似准则

陈荣盛等^[10] 以流体机械理论为基础，提出了风力机械相似准则，具体为风力机模型与风力机原型在空气能量传递过程与风流经过风力机过程中流动过程相似，包括几何相似，运动相似，力相似。风力机械相似准则将风力机特性用相关参数将模型和原型风力机联系起来，其中功率，水平推力之间的联系表现在比例因子$\alpha $，具体如下式：

式中：P——功率；

T——水平推力；

$ \alpha $——比例因子，$\alpha = {{{R}}_{{\rm{S}}}}/{{{R}}_{{\rm{M}}}}$为风力机原型大小跟风力机模型大小的比例。

本文采用原型风力机大小与模型大小比例为10的比例因子，先求出模型风力机的力矩，用比例因子求出原型风力机的力矩，进一步通过力矩与功率的关系求出风力机的功率。

1.6 网格无关性验证

采用瞬态计算不同数量的网格在额定风速下功率，提取不同数量网格与不同时间点的功率数据，如图6所示，结果表明：同一时间点上，当437万网格至766万之间，功率相比较437万之前越来越趋于平缓。同一数量网格上，437万网格相比较530万跟766万网格三条曲线的点较拟合，说明在10~20 s时间段功率较为平稳，综合考虑可认为437万网格之后，网格数量对功率影响变小，且在时间段上平稳，故本文采用437万网格进行模拟，保证网格无关性的同时提高计算效率。

1.7 功率和水平推力验证

本文用功率和水平推力参数进行模型验证。以单台风力机创建流场域，布置流场，采用稳态计算，对模型导入Fluent求解不同风速下5 m/s，8 m/s，11.4 m/s（额定风速），20 m/s风轮对旋转轴的转矩，计算风力机模型的功率和水平推力，依据美国可再生实验给出的5 MW风力机功率和水平推力，验证模型的可靠性。当风速变化时，特别是超过额定风速后，调整叶片的角度，实现变桨，控制风力机功率，在本文中，对超过额定风速的功率验证工况做改变风力机变桨角度的处理，具体变桨过程如图7所示。

依据文献[10]，得出了风力机对旋转轴力矩跟功率之间的公式关系。如下式：

式中：$ P $——输出功率，W；

$ M $——转矩，N·m；

${{\varOmega }}$——叶片转速，${\rm{r}}/\min $;

$ n $——叶片数，3。

经Fluent模拟，整个过程对旋转轴的力矩趋于稳定，模拟结果收敛后，求得风轮对x轴的力矩。依据式(8)、式(9)相似准则，力矩为缩小模型后的风轮对x旋转轴的力矩，原型风力机力矩需依据比例因子$\alpha $=10求出。具体如图8所示。

根据上图数据，结果表明，风力机功率数据与实验数据具有很小的误差，水平推力跟实验数据均偏大一些，五种具有代表性的风速工况下，均能满足风力机经过Fluent求解后的要求，五种风速下，整体误差都控制在10%内，产生误差的可能原因为风力机叶片表面翼型弯曲复杂，生成的非结构网格造成网格质量存在一定缺陷，未考虑风力机叶片截面的约束，都可能会造成误差，但误差控制在可接受的20%内，保证模型的有效性。

2 结果分析 2.1 工况设置

为研究不同间距对大型风力机尾流效应的影响，设置如下工况：大功率风力机采用串列方式，分别置于5D、8D、11D三种间距下的工况，其中D为风力机直径，尾流速度衰减在尾流发展4D～6D，设置最短间距5D开始，每一个工况间隔3D^[14] 。图9~图11分别表示5D、8D、11D下游风力机旋转区域的压力云图。图12~图14分别表示两台风力机在串列下前后距离为5D、8D、11D的速度云图，工况布置见表2。

2.2 下游旋转域压力分析

风场流入剪切来流风后，形成不同高度的风速层。从图9~图11可知来流风在经过有旋转速度的风力机后，形成类似于叶片形状的压力轮廓，有着明显的压力分层，形成这种现象的原因为风力机对风气流有一定的阻碍和扰动作用，扰动作用会使风力机叶片两面形成压力差，这种压力差和粘性力作用会造成风力机表面形成阻力，由于这种阻力一旦脱离风力机叶片，会使得粘性力在风力机表面损失，形成流动损失，使得下游风力机的尾流区速度亏损且持续发生剪切作用，进一步造成速度损失；相比于11D的间距，5D和8D受到上游风力机的气流扰动和压力影响更大，说明间距越小，这种压力影响作用越大；在11D的工况下，下游风力机入流扰动相比5D和8D更加趋于平缓；此外，三种工况下，共同特点是风力机轮毂处的压力大于叶表面，由于压差阻力之和大于其他部分，造成压力较大，导致轮毂处的能量转换比较强烈。

2.3 侧面速度云图分析

真实的风场中，存在差异性较大的速度梯度，进而形成湍流较强的剪切层；为了研究在剪切来流风情况下串列风力机不同间距对尾流特性的影响，因此截取两台风力机间距为5D、8D、11D三种工况下xz面，变量随着高度z变化的速度云图。由图12和图13对比可知：5D对比8D下，下游风力机入流处风速还在受到上游风速损失的影响，恢复不明显，来流速度的损失会影响下游风力机入流风速大小，同时下游风力机处于自身与上游风力机产生的湍流叠加区域，导致下游尾流形成膨胀流，虽然风力机间距大小不同，但是轮毂处的压强差很大，能量交换剧烈，导致中心处尾流损失较为明显；下游风力机6D之后由于不同梯度的速度的剪切层，尾流恢复区与未恢复区域的上层与下层区能量交换比较明显，会形成速度分离，与不同层速度梯度混合，最终速度趋于平缓，尾流恢复处（下游风力机速度损失恢复到来流环境的位置）发生在10D之后；从图13可得知：8D下游风力机的入流速度明显高于5D，下游风力机由于受到流动损失小于5D，导致了8D下游风力机尾流恢复处发生在9D左右，对比图12~图14得知：随着间距的增大，下游风力机受到的影响减小；综上所述：随着风力机上下游间距的增大，下游风力机受到上游风力机速度影响会越来越小，导致下游风力机在剪切层环境下，尾流恢复处发生的位置会越来越近，趋于正常的剪切层会越来越快。

3 结束语

1）使用剪切风入流，风速会随着高度变化而变化，风经过风力机扰动，会在风力机叶片形成压力分层，且这种压力扰动随着上下游风力机间距增大会越来越小。

2）风力机间距不同，但是共同特点是轮毂处压力大，能量转换强烈，使得轮毂后速度损失比叶轮其他地方严重，这一处速度损失应进一步加强研究。

3）基于大涡模拟，在剪切来流风下两台风力机采用串列方式布置，随着间距增大，尾流恢复处越近，说明间距越大，下游风力机受到上游风力机尾流速度影响越来越小。因此在风电场布置中，可建议参考尾流恢复处的距离来合理布置风力机间的距离。