偏航风力机尾流特性和功率性能的大涡模拟研究

Study on Large-eddy Simulation of Wake Characteristics and Power Performance of Yawed Wind Turbines

摘要

对偏航风力机尾流特性更好地了解有助于风电场偏航策略的使用，本文开展了偏航风力机尾流大涡模拟研究. 首先，建立基于伪谱法的大涡模拟耦合新发展的滤波致动盘模型的模拟方法，使用EPFL偏航风力机尾流风洞实验数据验证了该方法的可行性和准确性. 然后，模拟了不同偏航角下的风力机尾流，发现滤波致动盘模型对远尾流区的速度亏损、湍流度和尾流中心偏移有较高的预测精度. 与非偏航情形不同，偏航情形大涡模拟重现了实验中观察到的“卷曲尾流”现象，这种现象造成了速度亏损和湍流度在垂直面出现不对称的分布.最后，模拟了两台风力机的尾流状况，并分析上游风力机偏航角对功率性能的影响.模拟结果显示风力机间距不宜过短，一方面，上游风力机尾流干扰造成下游风力机的功率下降，另一方面，在较短间距情况下上游风力机尾流中心偏移值在下游风力机处较小，总功率增加量减少.研究内容对在风电场使用偏航技术具有指导意义.

Abstract

A better understanding of yawed wind turbine wake characteristics could facilitate the use of a yaw-control strategy in wind farms. In this research， large-eddy simulations of the yawed wind turbine wakes were carried out. First， a simulation method based on the pseudo-spectral-based large-eddy simulation， as well as the recently developed filtered actuator disk model， was established. The feasibility and accuracy of the method were validated by the EPFL wind tunnel measurements of the yawed wind turbine wakes. Then， the turbine wakes under different yaw angles were simulated. It was found that the filtered actuator disk model could predict the velocity deficit， turbulence intensity， and wake center deflection in the far-wake region with good accuracy. Different from the non-yaw case， a “curling wake” was observed in the simulated yaw cases， resulting in an asymmetric distribution of the velocity deficit and turbulence intensity in the vertical plane. Finally， the wake flows of two wind turbines were simulated， and the yaw angles of the upstream turbine on the power performance was analyzed. The results suggest that the wind turbine spacing should not be too short. On one hand， the power of the downstream wind turbines decreases due to the wake interference induced by the upstream wind turbine. On the other hand， in the case of shorter turbine spacing， the wake center deflection of the upstream wind turbines is smaller at the location of the downstream wind turbine， and therefore， the total power increase becomes reduced. The present analysis is of guiding significance for the use of yaw-control technology in the wind farms.

Graphical abstract

关键词

风力机 / 大涡模拟 / 大气边界层 / 偏航尾流 / 功率性能

Key words

wind turbine / large-eddy simulation / atmospheric boundary layer / yawed wakes / power performance

风能作为实现碳中和目标不可或缺的一部分，近年来发展迅速^［1］. 在风电场中，下游风力机会受上游风力机的尾流影响，导致显著的功率损失和疲劳荷载的增加，因此风力机尾流是风能领域的主流研究之一^［2-3］. 为了解决风力机尾流干扰导致风电场产能下降这一问题，近年来风力机偏航控制技术得到了研究人员的广泛关注^［4］.

针对偏航风力机尾流的研究手段主要有解析建模、实验（风洞和现场实测）以及数值模拟^［5-7］. Jiménez等^［8］首先提出尾流模型，其模型假设速度亏损服从“顶帽”分布. 随后Bastankhah等^［9］提出了一种基于高斯分布速度亏损的模型，发现该模型比Jiménez等^［8］提出的精度更高. Howland等^［10］、Bastankhah等^［9］以及Zong等^［11］开展了偏航尾流的风洞实验，发现了尾流形成两个偏转涡对. 美国国家可再生能源实验室的Fleming等^［12-13］在较小型风电场上开展了风力机偏航的实测实验，证实了风力机偏航对改善风场发电性能的有效性.

随着计算机技术的日益发展，大涡模拟已经成为风能领域研究风力机尾流效应和尾流干扰的常用手段^［14］. 由于大气边界层的尺度和风轮叶片的尺度差异较大，现今大涡模拟大气边界层下的风轮叶片尺度的计算量巨大^［15］. 因此，大多数关于风力机尾迹的数值研究，是用不同的参数化方法来表示风力机在大气边界层流动中引起的力.致动理论模型主要有致动盘^［16］、致动线^［17］和致动面模型^［18］. 曹九发等^［19］使用了致动线模型模拟了阵风工况下多台风力机尾流效应的非定常特性. Yang等^［18］使用致动面模型模拟了MEXICO风力机的尾流，与实验结果吻合较好. 后两者模型虽然被认为拥有更高的精度，但由于它们需要极大的网格量，对大型风电场的尾流模拟仍然是一个挑战^［20］. 在尾流信息捕捉方面，制动盘模型与其他模型相比略显不足，但其具备计算量少、精度高的特点，能够满足科研工程的需求，在风能尾流模拟中得到广泛运用^［20］. 关于致动盘模型最新的研究来自Shapiro等^［21］提出的滤波致动盘模型，该模型拥有较准确的功率预测性能而无需精细的数值网格或经验校正因子，被成功运用到大型风电场、城市风能和山脉风电场^［22-25］中，但目前尚缺乏在偏航风力机尾流情形下对该数值模型的有效性验证.

鉴于以上分析，本文使用大涡模拟耦合最新发展的滤波致动盘模型对偏航风力机尾流开展数值模拟，利用EPFL的偏航风力机尾流风洞实验数据验证模拟框架的准确性，并分析不同偏航角下的风力机尾流特性. 同时，模拟了两台风力机的尾流情形，并评估偏航角对风力机风电效率的影响.

1 数值模拟方法

1.1 控制方程

本文研究不考虑热分层和风偏转的影响，因此求解的三维不可压缩滤波N-S方程具体形式为：

\frac{\partial {\tilde{u}}_i}{\partial x_i}=\mathrm{0}（1）

\frac{\partial {\tilde{u}}_i}{\partial t}+{\tilde{u}}_j\frac{\partial {\tilde{u}}_i}{\partial x_j}=-\frac{\partial {\tilde{p}}^{\mathrm{*}}}{\partial x_i}-\frac{\partial {\tau }_{ij}^{\mathrm{d}}}{\partial x_j}+f_i+{\delta }_{i\mathrm{1}}{F}_{\mathrm{p}}（2）

式中：“˜”表示在尺度\tilde{\Delta }上进行空间滤波；t为时间；ũ_i 为瞬时速度，其中i=1、2、3分别为对应的顺流方向x、展向y和垂直方向z；{\tilde{p}}^{\mathrm{*}}为滤波压力；f_i 代表风力机对流体的作用力；F_{\mathrm{p}}=u_{\mathrm{*}}^{\mathrm{2}}/L_z是用来驱动预前模拟流动的压力梯度驱动力； {\tau }_{ij}代表亚格子应力，{\tau }_{ij}^{\mathrm{d}}为其偏应力. 因为求解的大气边界层流动具有高雷诺数的性质，所以黏性项从方程中省去. 本文使用的求解器程序来源于约翰·霍普金斯大学的湍流研究组，该开源求解器被广泛用于大气边界层^［26］、大涡模拟湍流模型^［27］和风电场流动^［20］研究. 为了封闭方程，亚格子应力采用Smagorinsky模型：

{\tau }_{ij}^{\mathrm{d}}=-\mathrm{2}{\upsilon }_{\mathrm{t}}{\tilde{\mathit{S}}}_{ij}（3）

{\upsilon }_{\mathrm{t}}={\left(C_{\mathrm{s}}\tilde{\Delta }\right)}^{\mathrm{2}}|\tilde{S}|（4）

式中：{\tilde{\mathit{S}}}_{ij}=\mathrm{0.5}(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i})为应变速率张量，其强度为|\tilde{S}|=\sqrt[]{\mathrm{2}}{\tilde{\mathit{S}}}_{ij}；C_{\mathrm{s}}为Smagorinsky系数. 在靠近底部壁面的区域采用阻尼函数：

\frac{\mathrm{1}}{{\lambda }^n}=\frac{\mathrm{1}}{{\lambda }_{\mathrm{0}}^n}+\frac{\mathrm{1}}{[k(z+z_{\mathrm{0}}{)]}^n}（5）

{\lambda }_{\mathrm{0}}=C_{\mathrm{0}}\tilde{\Delta }（6）

式中：k为冯·卡门常数；z₀为粗糙度长度；λ₀为Smagorinsky长度尺度.

1.2 风力机模拟方法

关于风力机的模拟，使用Shapiro等^［21］发展的滤波致动盘模型模拟风力机对流体的推力：

F_{\mathrm{t}}=-\frac{\mathrm{1}}{\mathrm{2}}\rho C_{\mathrm{T}}{A}_{\mathrm{D}}{\overline{u}}_{\mathrm{H}}^{\mathrm{2}}（7）

式中：F_t为风轮对流体的推力；C_T为推力系数；A_D为风轮扫掠面积；{\overline{u}}_{\mathrm{H}}^{}为轮毂高度处的平均风速. 该模型理论中，推力风轮的有效风速为：

u_{\mathrm{d}}=(\mathrm{1}-a){\overline{u}}_{\mathrm{H}}^{}\mathrm{c}\mathrm{o}\mathrm{s}\gamma（8）

a=\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}-\sqrt[]{\mathrm{1}-C_{\mathrm{T}}}\right)（9）

式中：a为诱导因子.

推力可被修改为：

F_{\mathrm{t}}=-\frac{\mathrm{1}}{\mathrm{2}}\rho {C_{\mathrm{T}}}^{\mathrm{\text{'}}}{A}_{\mathrm{D}}{u}_{\mathrm{d}}^{\mathrm{2}}（10）

式中：{C_{\mathrm{T}}}^{\mathrm{\text{'}}}=\frac{C_{\mathrm{T}}}{{(\mathrm{1}-a)}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\gamma }为修改推力系数.

为了防止数值求解出现不稳定状况，该模型使用指标函数R将推力分布到风轮周边网格点，其推导过程和具体公式都较为复杂，可参考文献［21］. 风力机的功率可由式（11）求得：

P=-F_{\mathrm{t}}{u}_{\mathrm{d}}（11）

1.3 离散方法

该求解器在各个方向上采用均匀的网格，因此网格尺寸为：

\mathrm{\Delta }x=\frac{L_x}{N_x}，\mathrm{\Delta }y=\frac{L_y}{N_y}，\mathrm{\Delta }z=\frac{L_z}{N_z}（12）

式中：L_x、L_y、L_z分别为计算域各个方向的长度；N_x、N_y、N_z分别为各个方向的网格点数量.垂直速度（w）的网格存储与其他变量（w、v、p）的网格在垂直方向上交错. 求解器在水平方向（x-y）采用拟谱法离散，在垂直方向上采用二阶中心差分法离散：

{\partial }_{x}A(x,y,z)=\sum _{k}i{k}_x\widehat{A}\left(k_x,k_y,z\right){\mathrm{e}}^{i\left(k_{x}x+k_{y}y\right)}（13）

{\partial }_{y}A(x,y,z)=\sum _{k}i{k}_y\widehat{A}\left(k_x,k_y,z\right){\mathrm{e}}^{i\left(k_{x}x+k_{y}y\right)}（14）

{\partial }_z{A}_{i,j,k}=\frac{A_{i,j,k}-A_{i,j,k-\mathrm{1}}}{\mathrm{\Delta }z}（15）

式中：A为任意变量；\widehat{A}为其傅里叶级数系数. 速度场采用二阶精度Adams-Bashforth格式对时间项进行离散，使用投影法来满足速度散度为零的条件.具体步骤为先计算中间速度：

\frac{{\tilde{u}}_i^{\mathrm{*}}-{\tilde{u}}_i^t}{\mathrm{\Delta }t}=\frac{\mathrm{3}}{\mathrm{2}}\mathrm{R}\mathrm{H}{{{\mathrm{S}}_i}^{\mathrm{\text{'}}}}^t-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{R}\mathrm{H}{\mathrm{S}}_i^{t-\mathrm{\Delta }t}（16）

式中：\mathrm{\Delta }t为时间步长；\mathrm{R}\mathrm{H}{{\mathrm{S}}_i}^{\mathrm{\text{'}}}和\mathrm{R}\mathrm{H}{\mathrm{S}}_i分别如下式所示：

\mathrm{R}\mathrm{H}{{\mathrm{S}}_i}^{\mathrm{\text{'}}}=-{\tilde{u}}_j\frac{\partial {\tilde{u}}_i}{\partial x_j}-\frac{\partial {\tau }_{ij}^{\mathrm{d}}}{\partial x_j}+f_i+{\delta }_{i\mathrm{1}}{F}_{\mathrm{p}}（17）

\mathrm{R}\mathrm{H}{\mathrm{S}}_i=\mathrm{R}\mathrm{H}{\mathrm{S}}_i^{\text{'}}-\frac{\partial \tilde{p}}{\partial x_i}（18）

然后，下一时间步的速度可以表示为：

\frac{{\tilde{u}}_i^{t+\mathrm{\Delta }t}-{\tilde{u}}_i^{\mathrm{*}}}{\mathrm{\Delta }t}=-\frac{\mathrm{3}}{\mathrm{2}}\cdot \frac{\partial {\tilde{p}}^t}{\partial x_i}（19）

由无散条件，可以得到压力泊松方程：

\frac{{\partial }^{\mathrm{2}}{\tilde{p}}^t}{\partial x_i\partial x_i}=\frac{\mathrm{2}}{\mathrm{3}\mathrm{\Delta }t}\cdot \frac{\partial {\tilde{u}}_i^{\mathrm{*}}}{\partial x_i}（20）

该方程可以写成以下形式：

\left[{\partial }_z^{\mathrm{2}}-\left(k_x^{\mathrm{2}}+k_y^{\mathrm{2}}\right)\right]\widehat{{\tilde{p}}^t}=\frac{\mathrm{2}}{\mathrm{3}\mathrm{\Delta }t}\left(i{k}_x\widehat{{\tilde{u}}_i^{\mathrm{*}}}+i{k}_y\widehat{{\tilde{v}}^{\mathrm{*}}}+{\partial }_z\widehat{{\tilde{w}}^{\mathrm{*}}}\right)（21）

该方程可以采用三对角线矩阵算法（TDMA）求解，相应的边界条件为：

{\partial }_z\widehat{p}=-\widehat{{\partial }_j{\tau }_{jz}}；z=0 或L_z（22）

1.4 边界条件

由于大气边界层具有高雷诺数的特性，在底面采用壁面函数方法，将最靠近底面的速度和底面切应力{\tau }_{i\mathrm{3}|\mathrm{w}}联系起来：

{\tau }_{i\mathrm{3}|\mathrm{w}}=-{\left(\frac{\kappa }{\mathrm{l}\mathrm{n}(z/z_{\mathrm{0}})}\right)}^{\mathrm{2}}\sqrt[]{{\widehat{\tilde{u}}}_{\mathrm{1}}^{\mathrm{2}}+{\widehat{\tilde{u}}}_{\mathrm{2}}^{\mathrm{2}}}{\widehat{\tilde{u}}}_i（23）

式中：“^”表示在尺度\mathrm{2}\tilde{\Delta }进行空间滤波.对于计算域顶部，采用零速度梯度和零应力边界条件：

\frac{\partial \tilde{u}}{\partial z}=\frac{\partial \tilde{v}}{\partial z}=\mathrm{0}，{\tau }_{xz}={\tau }_{yz}=\mathrm{0}（24）

求解器使用了同步预前模拟法^［28］生成对应实验的入口湍流. 为了消除在水平方向使用谱方法离散引入的周期边界条件的影响，求解器设置了一个缓冲区域，用一个三角函数将主域出口速度平滑过渡到预前模拟生成的入口速度，每一步的入口湍流由预前模拟域（不含风力机）通过MPI发送到主域（包含风力机）.为了方便起见，下文中滤波符号\tilde{\cdot }省略，\overline{\cdot }和〈\cdot 〉分别表示时间和水平方向取平均.

1.5 算例设置

Bastankhah等^［9］在洛桑联邦理工学院使用立体粒子图像测速（S-PIV）系统地开展了大气边界层偏航风力机尾流风洞实验.实验段高2.0 m，宽2.6 m，长28 m.实验中，轮毂高度处的平均风速{\overline{u}}_{\mathrm{H}}约为 \mathrm{4.88} m/s，粗糙度和摩擦速度分别为z₀=0.000 022 m和u_*=0.19 m/s（基于将对数速度剖面拟合到实验中边界层下部15%部分），大气边界层高度约为0.4 m. 使用的风力机模型是三叶片的水平轴风力机，风轮直径为D=0.15 m，轮毂高度为H=0.125 m［图1（a）］. 计算域尺寸设置如图1（b）所示，L_x =4.2 m（28D），L_y =0.825 m（5.5D），L_z =0.4 m（2.67D），采用的网格数量为N_x =232，N_y =56，N_z =36. 风力机布置在距离大气湍流入口0.6 m（\mathrm{4}D）位置. 采用的网格尺寸在先前的研究工作中已被证实可以达到网格无关性的要求^［25］.

2 风力机尾流模拟分析与讨论

2.1 大气边界层入口模拟

图2展示了预前模拟生成的大气边界层入口特性. 由于模拟的大气边界层流动是由压力梯度驱动力F_{\mathrm{p}}=u_{\mathrm{*}}^{\mathrm{2}}/L_z驱动的，所以入口应符合理论公式：

-\frac{〈\overline{u^{\text{'}}{w}^{\text{'}}}〉}{u_{\mathrm{*}}^{\mathrm{2}}}-\frac{〈\overline{{\tau }_{\mathrm{13}}^{\mathrm{d}}}〉}{u_{\mathrm{*}}^{\mathrm{2}}}=\mathrm{1}-\frac{z}{L_z}（25）

图2（a）展示了模拟值与理论值吻合得很好，其可以作为模拟的大气边界层达到准平衡状态的依据.从图2（b）可以发现，靠近大气边界层顶部的速度模拟值与实验值相比有所降低，但在轮毂高度附近整体吻合良好.图2（c）展示了流向湍流度的结果，其和实验值吻合较好.为了检验预前模拟生成的速度入口的自保持性，图3展示了流向速度和流向湍流度在不同位置处的垂直剖面.平均速度剖面［图3（a）］在下游各处保持不变，但在靠近地面的湍流度［图3（b）］峰值随着下游距离的增加而小幅衰减，而在其他高度变化不明显.因此，可以认为模拟的大气边界层能满足自保持性要求.正确模拟和生成大气边界层入口是数值模拟成功的前提，风工程领域也常用人工合成法生成大气边界层入口湍流（如CDRFG），可参考文献［29］.

2.2 风力机尾流速度亏损和湍流度验证

本文使用大涡模拟了实验中\gamma =-\mathrm{10}°、-\mathrm{20}°、-\mathrm{30}°的情况，实验中测得风轮推力系数为C_{\mathrm{T}}^{\mathrm{E}\mathrm{x}\mathrm{p}}=0.788、0.738、0.666. 但由于参数设置须指定局部推力系数，根据公式{C_{\mathrm{T}}}^{\mathrm{\text{'}}}=\frac{C_{\mathrm{T}}}{\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\gamma {(\mathrm{1}-a)}^{\mathrm{2}}}可分别设置为1.524、1.463、1.427.图4展示了流向平均速度亏损在轮毂高度处不同下游位置的展向剖面.可以发现，在不同偏航角情况下近尾迹区域（x/D=4）速度亏损的模拟值较实验值偏低.

图5为流向湍流度在轮毂高度处不同下游位置的展向剖面.同样地，湍流度的模拟值在近尾迹区域明显低于实验值.这是因为实际状况该区域会受轮毂和叶尖涡的影响，尾流结构较复杂，而滤波致动盘数值模型不能反映它们的影响.但在更下游的位置，它们的模拟值与实验值吻合较好.由于风力机的间距在实际的风电场布局一般保持在x/D≥5，远尾迹区的尾流特性是研究关注的重点，因此该大涡模拟框架具有较高的精度.需要注意的是，对于偏航风力机尾流中的湍流度特性，其在轮毂高度处的展向剖面不会呈现对称性分布，最大值出现在尾流偏移的一侧.未来针对湍流度尾流模型的建立需要考虑这一特性.

图6为不同偏航角下的尾流最大速度亏损和尾流中心偏移. 随着偏航角的增大，速度亏损强度变小，尾流中心偏移增大，这对下游风力机的功率提高是有益的. 另外，尾流中心偏移量并不会随着x/D的增加而一直增加，而是存在一个渐近极限.尾流中心偏移量会直接影响下游风力机的功率输出. 从图6可以观察到，滤波致动盘模型在近尾迹区会较低预测最大速度亏损，其在远尾迹区的尾流中心偏移预测值与实验值吻合较好.本文采用Bastankhah模型预测最大速度亏损和尾流中心偏移，发现其在远尾迹区的最大速度亏损和尾迹中心偏移的预测较好，在一定程度上也反映了数值模拟的准确性.尾流解析模型属于流体降阶模型，其优点是比数值模拟的计算速度快，能被运用于实际工程项目，其缺点是该模型基于尾流速度亏损呈轴对称性的假设，但实际情况偏航尾流速度亏损不会呈轴对称性分布（见图7）.

2.3 风力机尾流云图

图7为不同偏航角下不同下游位置的尾流速度亏损云图（白色线为展向平面流线）.非偏航情况下{C_{\mathrm{T}}}^{\mathrm{\text{'}}}=1.631，非偏航和偏航情况下的风力机的速度亏损剖面明显不同，偏航情况下表现为左右不对称的肾形，称为“卷曲尾流”（curled wakes）.该现象首次被Howland等^［10］的风洞实验观察到，该现象能被大涡模拟技术重现.观察尾流的流线情况，可以看到在偏航状态风力机的上下两端附近会形成一个对转涡对.该对转涡对在尾流中心处诱导了一个非均匀的横向速度，将高速度亏损区域输运到一侧，形成了向内“卷曲”的尾流剖面.现有的解析模型能准确预测最大尾流亏损和尾流中心偏移，但如何采用数学模型去描述这种尾流剖面是目前的难题.

图8为不同偏航角下不同下游位置的流向湍流度云图. 在非偏航情况下，可以观察到风轮边缘附近湍流度较高，风轮上部的湍流度比下部高，这是因为受尾流速度亏损影响，风轮上部的流向速度剪切较大，促使湍动能生成，湍流度剖面呈左右对称而上下不对称. 在偏航情况下，由于尾流发生了卷曲，湍流度展现出左右非对称性，在尾流偏转一侧更高. 在相同的下游位置，观察到尾流湍流度随着风力机偏航角的增大而降低. 这是因为较高偏航角的风力机推力较小，对气流干扰也变小.

3 偏航尾流对下游风力机功率性能的影响

为了衡量上游风力机的偏航角对下游风力机功率性能的影响，本文考虑两台间距为x/D=4、6、8的风力机.两台风力机的局部推力系数设置为{C_{\mathrm{T}}}^{\mathrm{\text{'}}}=\mathrm{1.463}.图9为风力机间距为x/D=6的尾流速度云图. 由于上游风力机尾流发生偏转，下游风力机有效风速增加.下游风力机虽然没有偏航，但其尾流中心仍会偏移，且随着上游风力机偏航角的增大而增大.原因是上游风力机尾流形成的旋转涡在下游流场诱导了横向速度.这对大型风电场的偏航控制是有益的，能降低下游风力机受上游风力机尾流的影响，且上游风力机不必因为调整过大的偏航角而降低自身的发电性能.

图10比较了不同风力机间距下的风力机总平均功率、下游风力机平均功率和下游风力机平均推力.图中{\overline{P}}_{\mathrm{0}}、{\overline{P}}_{\mathrm{1}}、{\overline{P}}_{\mathrm{2}}分别表示上游风力机无偏航平均功率、上游风力机平均功率和下游风力机平均功率，{\overline{T}}_{\mathrm{0}}、{\overline{T}}_{\mathrm{2}}分别表示上游风力机无偏航的平均推力和下游风力机平均推力.如图10（a）所示，上游风力机采取了偏航策略后总功率提高，这意味着下游风力机的功率提高补偿了上游风力机的功率损失.x/D=4时，总功率提高最少，这是因为上游风力机的尾流偏转量在下游风力机处较小，上游风力机的尾流仍较大程度地影响下游风力机的功率. x/D=4、6、8时，最大总功率分别在γ=20°、25°、20°出现（以5°为分辨率）.该模拟结果提示若使用静态尾流模型对风电场进行偏航优化，风力机的偏航角范围可设置在30°以内，一方面，若偏航角范围设置得过大，利用梯度算法搜寻最优偏航角会耗费更多的计算时间；另一方面，偏航尾流解析模型是基于小偏航角的假设，风力机偏航角过大会造成尾流模型预测出现偏差. 如图10（b）所示，随着间距的增大，上游风力机尾流在下游风力机处速度亏损减小，下游风力机的平均功率增加，当x/D=4时，下游风力机的功率增量比其他情况少，这是因为在风机间距较小情况下，上游风力机的尾流中心偏移值在下游风力机处较小. 图10（c）所示的下游风力机平均推力有着和其平均功率相似的变化规律，这是因为滤波致动盘数值模型中风力机的功率是用公式P=-F_{\mathrm{t}}{u}_{\mathrm{d}}计算得到. 虽然本文关注的重点是偏航风力机的尾流特性和功率性能，但偏航风力机和受偏航尾流干扰的风力机的结构动力响应值得引起关注.

4 结 论

本研究建立基于伪谱法的大涡模拟耦合滤波致动盘模型的模拟方法，使用EPFL风洞实验数据验证了数值模拟的可行性和准确性，并研究了不同偏航角下风力机尾流特性和功率性能. 结论如下：

1）基于滤波致动盘模型模拟风轮作用力的大涡模拟框架，预测的近尾迹区的速度亏损和湍流度偏小，但对预测远尾迹区速度亏损、湍流度和尾流中心偏移具有较高精度. 现实中由于下游风力机布置在远尾迹区域，因此该大涡模拟框架具有较高的实用价值.

2）随着偏航角的增大，风力机尾流速度亏损变小，而尾流中心偏移变大，这有利于下游风力机捕获风能.

3）该大涡模拟框架很好地复现偏航情况下的“卷曲尾流”现象. 由于尾流发生卷曲，平均速度亏损和湍流度在垂直面分布不对称，湍流度在尾流偏转一侧更高.

4）采取偏航策略能提高总功率，下游风力机功率的提高能补偿上游风力机功率的损失. 风力机间距不宜过短，一方面是因为上游风力机尾流干扰造成下游风力机的功率下降；另一方面是因为上游风力机尾流中心偏移值在下游风力机处较小，导致总功率增加量较小.

参考文献

[1] PORTÉ-AGEL F， BASTANKHAH M， SHAMSODDIN S ．Wind-turbine and wind-farm flows：a review［J］．Boundary-layer Meteorology， 2020， 174（1）： 1-59．

[2] 王同光，田琳琳，钟伟，等 ．风能利用中的空气动力学研究进展Ⅰ：风力机气动特性［J］．空气动力学学报，2022，40（4）：1-21．

[3] WANG T G， TIAN L L， ZHONG W，et al. Aerodynamic research progress in wind energyⅠ：wind turbine aerodynamic characteristics［J］. Acta Aerodynamica Sinica， 2022， 40（4）： 1-21．（in Chinese）

[4] 王同光， 田琳琳， 钟伟， 等 ．风能利用中的空气动力学研究进展Ⅱ：入流和尾流特性［J］．空气动力学学报，2022， 40（4）：22-50．

[5] WANG T G， TIAN L L， ZHONG W，et al ．Aerodynamic research progress in wind energy Ⅱ：inflow and wake characteristics［J］．Acta Aerodynamica Sinica，2022，40（4）：22-50．（in Chinese）

[6] 宗豪华，孙恩博 ．水平轴风力机主动尾流控制综述［J］．空气动力学学报，2022，40（4）：51-68．

[7] ZONG H H， SUN E B ．Reivew of active wake control for horizontal-axis wind turbines［J］．Acta Aerodynamica Sinica，2022，40（4）：51-68．（in Chinese）

[8] 张绍海， 高晓霞， 朱霄珣， 等. 考虑风切变影响的三维尾流模型风场实验［J］. 空气动力学学报， 2023， 41（11）： 71-79.

[9] ZHANG S H， GAO X X， ZHU X X， et al. Experimental study on wind field of three-dimensional wake model considering the influence of wind shear［J］. Acta Aerodynamica Sinica， 2023， 41（11）： 71-79. （in Chinese）

[10] ZONG H H， PORTÉ-AGEL F ．Experimental investigation and analytical modelling of active yaw control for wind farm power optimization［J］．Renewable Energy，2021，170：1228-1244．

[11] ISHIHARA T， QIAN G W ．A new Gaussian-based analytical wake model for wind turbines considering ambient turbulence intensities and thrust coefficient effects［J］．Journal of Wind Engineering and Industrial Aerodynamics，2018，177：275-292．

[12] JIMÉNEZ Á， CRESPO A， MIGOYA E ．Application of a LES technique to characterize the wake deflection of a wind turbine in yaw［J］．Wind Energy，2010，13（6）：559-572．

[13] BASTANKHAH M， PORTÉ-AGEL F. Experimental and theoretical study of wind turbine wakes in yawed conditions［J］. Journal of Fluid Mechanics， 2016， 806： 506-541．

[14] HOWLAND M F， BOSSUYT J， MARTÍNEZ-TOSSAS L A，et al ．Wake structure in actuator disk models of wind turbines in yaw under uniform inflow conditions［J］. Journal of Renewable and Sustainable Energy， 2016， 8（4）：043301．

[15] ZONG H H， PORTÉ-AGEL F ．A point vortex transportation model for yawed wind turbine wakes［J］. Journal of Fluid Mechanics， 2020， 890： A8．

[16] FLEMING P， KING J， DYKES K， et al. Initial results from a field campaign of wake steering applied at a commercial wind farm-part 1［J］. Wind Energy Science，2019，4（2）：273-285．

[17] FLEMING P， KING J， SIMLEY E， et al. Continued results from a field campaign of wake steering applied at a commercial wind farm-part 2［J］．Wind Energy Science， 2020， 5（3）： 945-958．

[18] STEVENS R J A M， MENEVEAU C. Flow structure and turbulence in wind farms［J］．Annual Review of Fluid Mechanics，2017，49：311-339．

[19] CHURCHFIELD M J， LEE S， MICHALAKES J，et al. A numerical study of the effects of atmospheric and wake turbulence on wind turbine dynamics［J］. Journal of Turbulence， 2012， 13：N14．

[20] WU Y T， PORTÉ-AGEL F ．Large-eddy simulation of wind-turbine wakes：evaluation of turbine parametrisations［J］．Boundary-Layer Meteorology，2011，138（3）：345-366．

[21] SO̸RENSEN J N， SHEN W Z ．Numerical modeling of wind turbine wakes［J］．Journal of Fluids Engineering，2002，124（2）：393-399．

[22] YANG X L， SOTIROPOULOS F ．A new class of actuator surface models for wind turbines［J］．Wind Energy，2018，21（5）：285-302．

[23] 曹九发，宋佺珉，王超群，等 ．阵风工况下多台风力机尾流效应的非定常特性［J］．空气动力学学报，2022，40（4）：247-255．

[24] CAO J F， SONG Q M， WANG C Q，et al ．Unsteady characteristics of wake effect for multiple wind turbines under gust wind condition［J］．Acta Aerodynamica Sinica，2022，40（4）：247-255．（in Chinese）

[25] STEVENS R J A M， MARTÍNEZ-TOSSAS L A， MENEVEAU C ．Comparison of wind farm large eddy simulations using actuator disk and actuator line models with wind tunnel experiments［J］．Renewable Energy，2018，116：470-478．

[26] SHAPIRO C R， GAYME D F， MENEVEAU C ．Filtered actuator disks：theory and application to wind turbine models in large eddy simulation［J］．Wind Energy，2019，22（10）：1414-1420．

[27] GE M W， GAYME D F， MENEVEAU C ．Large-eddy simulation of wind turbines immersed in the wake of a cube-shaped building［J］．Renewable Energy，2021，163：1063-1077．

[28] ZHANG S B， DU B W， GE M W，et al ．Study on the operation of small rooftop wind turbines and its effect on the wind environment in blocks［J］．Renewable Energy，2022，183：708-718．

[29] LIU L Q， STEVENS R J A M ．Effects of two-dimensional steep hills on the performance of wind turbines and wind farms［J］．Boundary-Layer Meteorology，2020，176（2）：251-269．

[30] ZHANG Z Y， HUANG P， BITSUAMLAK G，et al ．Large-eddy simulation of wind-turbine wakes over two-dimensional hills［J］．Physics of Fluids，2022，34（6）：065123．

[31] ALBERTSON J D， PARLANGE M B ．Surface length scales and shear stress：implications for land-atmosphere interaction over complex terrain［J］．Water Resources Research，1999，35（7）：2121-2132．

[32] BOU-ZEID E， MENEVEAU C， PARLANGE M ．A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows［J］．Physics of Fluids，2005，17（2）：025105．

[33] STEVENS R J A M， GRAHAM J， MENEVEAU C ．A concurrent precursor inflow method for large eddy simulations and applications to finite length wind farms［J］．Renewable Energy，2014，68：46-50．

[34] 周桐，闫渤文，杨庆山，等 ．大气边界层大涡模拟入口湍流生成方法研究［J］．工程力学，2020，37（7）：68-76．

[35] ZHOU T， YAN B W， YANG Q S， et al ．Study of inflow turbulence generation methods with large eddy simulation for atmospheric boundary layer［J］. Engineering Mechanics，2020， 37（7）：68-76．（in Chinese）

基金资助

国家自然科学基金资助项目(52178500)

National Natural ScienceFoundation of China(52178500)

国家留学基金委资助项目(202106260167)

Chinese Scholarship Council Foundation(202106260167)