CN104239622B - 风力机尾流计算方法 - Google Patents

Abstract

本发明公开了风力机尾流计算方法，属于风力发电的技术领域。计算方法，利用制动盘模拟风力机，建立包含制动盘动量源项、机舱附加轴向动量源项、湍流动能源项、湍流动能耗散率源项的尾流求解模型，设置中性大气边界层条件和CFD求解参数；按中性大气边界层确定常数的取值，迭代求解过程确定尾流区风速、湍流度的分布值。本发明考虑到制动盘动量源项、机舱附加轴向动量源项、湍流动能源项、湍流动能耗散率源项对风力机尾流计算的影响，修正时均雷诺N‑S方程以及方程后组成风力机尾流求解模型，在远尾流区域内自适应源项模型、固定源项模型的径向分布相同，改进型模型则进一步改进了近尾流区域的径向分布，提高了模型的计算精度。

Description

风力机尾流计算方法

技术领域

本发明公开了风力机尾流计算方法，属于风力发电的技术领域。

背景技术

当空气流经旋转的风轮时，由于受到叶片的阻碍和搅拌作用，会在风轮后形成尾流区。该区域内，风速相对于来流有所降低，湍流增强，从而影响尾流区内的后排风机。对于大型风场，这种影响可造成总出力损失10％～20％。因此在风场出力计算中，必须考虑尾流效应。

风力机尾流受到来流、风机气动特性以及下垫面等多种因素的影响。基于半经验的公式模型，如N.O.Jensen提出的Jensen模型和G.C.Larsen提出的Larsen模型，对这些因素做了很多简化，不能满足复杂地形条件下的计算需要。如果考虑全尺度的CFD模型，直接对风力机建模，操作困难，计算量大。为了简化CFD计算，研究者将风轮等效成制动盘、升力线或升力面，并结合时均雷诺方程，求解尾流场。

本文研究制动盘结合改进的k-ε湍流模型，模拟风力机尾流的数值计算方法。标准的k-ε湍流模型适用于大部分的工业计算，并不适合计算风力机尾流，导致尾流风速明显高于实验值，造成尾流恢复过快。El Kasmi和Christian Masson等人基于Chen和Kim的拓展k-ε模型，通过添加湍流耗散率源项，使湍流产生率与耗散率协调。El Kasmi模型建立在将轴流简化的基础上，即认为流动关于风轮旋转中心轴对称，没有考虑风速切边和地表粗糙度。此外，模型对动量源项的计算依赖于参考风速，在计算风场尾流时该风速难以确定。这些将阻碍模型推广到整场尾流计算。丹麦实验室的Pierre-Elouan Réthoré将森林冠层模型^[13]用于尾流模拟。与El Kasmi模型不同是，冠层模型除了耗散率源项外，还附加了湍流动能源项，并且两种源项都与冠层拖拽系数成正比。D.Cabezón等人^[9]考虑了湍流在空间上的各向异性，使用雷诺应力模型求解风机尾流，得到了与测量较匹配的结果。但是应用雷诺应力模型计算量要比二方程湍流模型增加很多，实用性差。

发明内容

本发明所要解决的技术问题是针对上述背景技术的不足，提供了风力机尾流计算方法。

本发明为实现上述发明目的采用如下技术方案：

风力机尾流计算方法，

步骤1，利用制动盘模拟风力机，建立包含制动盘动量源项S_u、机舱附加轴向动量源项S_d、湍流动能源项S_k、湍流动能耗散率源项S_ε的尾流求解模型：

其中，ρ为流体密度，u_i、u_j为平均速度在各坐标轴的分量，p为压强，S_ij为平均应变张量，u_i′、u_j′为速度在各坐标轴的波动分量，μ为流体粘性系数，μ_t为湍流粘度系数，k为湍流动能，P_k为湍流动量生成率，ε为湍流动能耗散率，C_1ε、C_2ε、C_μ、C_x、C_4ε、σ_k、σ_ε为常数；

步骤2，设置中性大气边界层条件：速度入口u、湍流动能k、湍流动能耗散率ε、地面粗糙长度K_S；

步骤3，设置CFD求解参数；

步骤4，按中性大气边界层确定常数C_1ε,C_2ε,C_μ,σ_k,σ_ε的取值，开始迭代求解过程确定尾流区风速、湍流度的分布值。

作为所述风力机尾流计算方法的进一步优化方案，考虑源项的径向分布参数η_r，进一步修正制动盘动量源项S_u、湍流动能耗散率源项S_ε：

η_ε0、η_ε1为常数，η_ε0＝0.0，η_ε1＝0.9，为网格中心到制动盘径向距离r对制动盘半径D的归一量，

进一步的，所述风力机尾流计算方法步骤3中设置CFD求解参数时，湍流粘性比上限为2×10⁷。

进一步的，所述风力机尾流计算方法步骤2中的速度入口u、湍流动能k、湍流动能耗散率ε由如下表达式确定：

其中，u_*为摩擦速度，z为离地高度，z₀为大气粗糙度长度，表示近地层风速向下递减到零时的高度，κ为冯卡门常数，κ＝0.4187，λ为模型常数，对于各向同性湍流，λ取值为1.5，H为轮毂高度，U₀、I₀分别为轮毂高度H处的风速、湍流强度。

进一步的，所述风力机尾流计算方法步骤2中所述地面粗糙长度K_S的表达式为：E＝9.793，C_S＝1.0。

本发明采用上述技术方案，具有以下有益效果：考虑到制动盘动量源项、机舱附加轴向动量源项、湍流动能源项、湍流动能耗散率源项对风力机尾流计算的影响，修正时均雷诺N-S方程以及k-ε方程后组成风力机尾流求解模型，在远尾流区域内自适应源项模型、固定源项模型的径向分布相同，改进型模型则进一步改进了近尾流区域的径向分布，提高了模型的计算精度。

附图说明

图1为全尺度计算区域示意图。

图2为轴流简化计算区域示意图。

图3为参数C_4ε对中心轴速度分布的影响。

图4为轴流假设验证结果图。

图5(a)、图5(b)、图5(c)是Nibe-B风力机入流速度8.5m/s,到风轮的轴向距离相对风轮之间的比分别为2.5、6、7.5时的尾流区径向相对速度对比示意图。

图6(a)、图6(b)、图6(c)是Nibe-B风力机入流速度9.56m/s，到风轮的轴向距离相对风轮之间的比分别为2.5、6、7.5时的尾流区径向相对速度对比示意图。

图7(a)、图7(b)、图7(c)是Nibe-B风力机入流速度11.52m/s，到风轮的轴向距离相对风轮之间的比分别为2.5、6、7.5时的尾流区径向相对速度对比示意图。

图8(a)、图8(b)、图8(c)是Dawin风力机入流速度8m/s，到风轮的轴向距离相对风轮之间的比分别为1、4.15、9.4时不同轴向位置处尾流区径向相对速度对比示意图。

图9为本发明的流程图。

具体实施方式

下面结合附图对发明的技术方案进行详细说明。

1、建立风机尾流CFD计算模型

根据风力机的几何参数，如风轮直径D、轮毂高度H等信息建立风力机尾流的CFD计算区域，如图1同时考虑风速随高度切变、地形起伏、地面粗糙度等多种因素，属于三维模型，利用标准壁面函数Wall模拟地表，入口速度velocity inlet为大气边界条件，出口设置成自由出流Outflow，执行器Actuator结合标准壁面函数得到制动盘动量源项S_u、湍流动能源项S_k、湍流动能耗散率源项S_ε。

划分网格：水平方向在风轮区处加密；高度方向，从地面往上按比例尺寸递增，并控制紧邻地面的网格高度l₀满足：

l₀>2K_S (1)

式中，K_S为地面粗糙长度，由大气粗糙长度按下式计算：

其中，E＝9.793,C_S＝1.0。

2、设置中性大气边界层条件

在风电场中，风力机的风轮与大气之间相互作用，因而风力机尾流计算必须模拟大气边界层，按照中性大气边界层处理。采用速度入口：

一般入流条件只给定轮毂高度H处的风速U₀和湍流强度I₀，此时有：

k＝λ(U₀I₀)² (6)

z₀＝H/exp(κU₀/u_*) (8)

式中，λ为模型常数，对于各向同性湍流取1.5。

出口按自由出流设置；两侧面为对称边界；地表利用标准壁面函数模拟，粗糙长度按式(2)计算。

3、湍流模型改进与源项添加

3.1控制方程

本文中使用的控制方程为时均雷诺N-S方程以及k-ε方程。对于风机在某特定工况下的尾流，按定常不可压缩流场计算。因此，有：

式中，μ为流体的粘性系数；为平均应变张量；S_u,S_d为制动盘和机舱的附加的轴向动量源项，仅在相应的区域不为0；为雷诺应力，按Boussinesq假设计算：

其中，为湍流粘度系数。

P_k为湍流动能的生成率：

P_k＝μ_tS² (14)

其中，为应变率张量的模量。

模型常数C_1ε,C_2ε,C_μ,σ_k,σ_ε按中性大气边界层取值：

C_1ε＝1.176,C_2ε＝1.92,C_μ＝0.033,σ_k＝1,σ_ε＝1.3。

S_k,S_ε表示湍流动能k及其耗散率ε的方程源项，只在相应的条件和区域内有效。

3.2动量源项

风力机制动盘模型，将风轮视为空气可自由通过，具有一定渗透阻力的圆柱型薄盘。其直径与风轮直径相当，设为D；厚度记作Δx；圆盘中心与轮毂中心重合。下面应用一维制动盘模型，推导动量源项的计算式。

对于来流风速U₀，风轮上单位面积上受到的推力为：

式中，ρ为空气密度；C_T为当前来流风速下风轮的推力系数，根据推力系数曲线求出。由此得到施加在制动盘上，单位体积的应力为：

上式中，源项直接由推力系数和参考风速决定，在整个计算区域中均匀分布，且在整个迭代过程保持恒定，称为固定源项。但是对于整场尾流，处于尾流区的后排风机，其参考风速并没有标准可以确定。因此，按照(16)式计算动量源项推广到整场尾流计算时将存在很多困难。将此项改写成局部风速u的关系式。这种由当地值确定的源项，称为自适应源项。如果以a表示轴向诱导因子，则，

u＝U₀(1-a) (17)，

又，

C_T＝4a(1-a) (18)，

于是，动量源项：

其中，C_x＝4a/((1-a)·Δx)，为制动盘拖拽系数。

上面推导出的动量源项，均匀分布在制动盘内，会在制动盘边缘出现体积力突变，将导致在数值计算中速度与压强之间的脱钩，使得压强与毗邻网格无关。为此，引入径向分布函数η_r，使得：

其中，径向分布函数为：

式中，为网格中心到制动轴向的径向距离r对制动盘半径的归一量，即，

模型参数D₀，D₁由用户指定。如果已知η_ε(0)＝η₀,η_ε(1)＝η₁<1且η₀+η₁≠2，则：

机舱近似成圆盘扰流，其直径为D_nac，阻力系数为C_d，由此附加的动量源项为：

其中，C_d在0.8到1.2之间。

3.3改进k-ε湍流模型

风轮的扰动使得风轮临近区域的流场变得很复杂，加速了湍流动能的产生与耗散作用，降低湍流尺度。添加湍流耗散率源项S_ε，使湍流产生率与耗散率协调。S_ε表示从大尺度湍流向小尺度湍流传递的能量传递速率：

S_ε需要施加在风轮邻域，该临域为对称分布在风轮上下游的圆柱区域，直径与制动盘相同，且共轴，厚度为L，取L＝0.5D。

参数C_4ε默认值为0.37，与改写动量源项后的模型不匹配(参见图3，x_c表示风轮中心在x轴上的坐标，(x-x_c)/D表示到风轮的轴向距离相对风轮之间的比值)。本文在径向上对C_4ε做修正：

其中，η_ε0,η_ε1表示C_4ε在制动盘中心和边沿的取值，此处取η_ε0＝0.0,η_ε1＝0.9。

3.4源项添加

添加源项可根据动量源项的不同分成三种类型：固定源项模型、自适应源项模型和改进模型三种尾流模型(见附表1)。本专利使用改进模型进行风力机尾流计算。

4、设置其它CFD求解条件

中性大气边界层中，其湍流粘性系数非常大，湍流粘性比要大于软件的默认上限。为排出求解器影响，本文将湍流粘性比上限调到2×10⁷。数值算法采用SIMPLE，其它相关量都采用二阶迎风格式求解。

5、迭代求解与后处理

风机尾流的计算域有两种：完成以上操作后，使用CFD求解器计算风力机尾流；导出水平面上的速度分布云图等进行后处理。

6、改进尾流模型验证

6.1验证对比试验参数

使用Nibe-B为验证改进模型的准确性，本文对两种类型的风机尾流进行了数值模拟：(1)三叶片Nibe-B 630kW风力机，分别模拟三种入流风速下的尾流场；(2)Dawin 180kW风力机。两种风力机参数见附表2。

6.2轴对称流动简化

如果不考虑风速切变，认为整个入流特征(U₀,k,ε)在径向上均匀分布，取轮毂中心高度H处的入流特征(U₀(H),k(H),ε(H))作为整个计算域的入流条件。此时，模型可简化成关于制动盘中心对称的轴流进行模拟。该情况下，需要附加额外的湍流源项，从而维持风轮及其尾流场外部湍流保持均匀性：

轴对称简化的模型与全尺度模型的计算对比，结果见图4。从图中可以看出单机尾流模拟中，两种计算域对速度分布的相对比值影响较小，轴对称简化应用单机尾流模型的优化和验证是合适的。

6.3网格划分与求解设置

轴流计算域如图2所示，附加相应简化条件，将问题从三维降阶为二维，轴向长25D，径向5D直径；制动盘位于入口下游5D处。网格在径向上均匀分布，长度为1m；轴向上风轮及其临域内网格长度1m，并由此向上、下游过渡，使上、下游两部分的网格数为100和200；总网格数6.5万。

本发明涉及的风力机尾流计算方法如图7所示，

步骤1，利用制动盘模拟风力机，建立包含制动盘动量源项S_u、机舱附加轴向动量源项S_d、湍流动能源项S_k、湍流动能耗散率源项S_ε的尾流求解模型：

ρ为流体密度，u_i、u_j为平均速度在各坐标轴的分量，p为压强，S_ij为平均应变张量，u_i′、u_j′为速度在各坐标轴的波动分量，μ为流体粘性系数，μ_t为湍流粘度系数，k为湍流动能，P_k为湍流动量生成率，ε为湍流动能耗散率，C_1ε、C_2ε、C_μ、C_x、C_4ε、σ_k、σ_ε为常数；

考虑源项的径向分布参数η_r，进一步修正制动盘动量源项S_u、湍流动能耗散率源项S_ε：

步骤2，设置中性大气边界层条件：速度入口u、湍流动能k、湍流动能耗散率ε、地面粗糙长度K_S，

速度入口u、湍流动能k、湍流动能耗散率ε由如下表达式确定：

其中，u_*为摩擦速度，z为离地高度，z₀为大气粗糙度长度，表示近地层风速向下递减到零时的高度，κ为冯卡门常数，κ＝0.4187，λ为模型常数，对于各向同性湍流，λ取值为1.5，H为轮毂高度，U₀、I₀分别为轮毂高度H处的风速、湍流强度，

地面粗糙长度K_S的表达式为：E＝9.793，C_S＝1.0；

步骤3，设置CFD求解参数；

步骤4，按中性大气边界层确定常数C_1ε,C_2ε,C_μ,σ_k,σ_ε的取值，开始迭代求解过程确定尾流区风速、湍流度的分布值。

6.4模拟结果与模型验证

表1

风机类型

额定功率

额定转速

风轮直径

轮毂高度

机舱直径

邻域长度

Nibe‐B

630kW

33rpm

40m

45m

5m

20m

Dawin

180kW

40rpm

23m

31m

2m

10m

表2

对表2所示的两种风机模型，利用表1所示的三种模型进行模拟，图5(a)、图5(b)、图5(c)、图6(a)、图6(b)、图6(c)、图7(a)、图7(b)、图7(c)表示Nibe-B630kW风力机尾流在三种入流风速情况下应用不同尾流模型的计算结果对比。由于在模型调整中使用了轴对称假设，计算域简化成二维区域，因此只需要比较各轴向位置处风速沿径向的分布情况。从图中可以看出，当x-x_c＝2.5D时，自适应源项模型的计算结果与实验值存在较大差异，因此必须改进自适应源项模型。当x-x_c＝6D或x-x_c＝7.5D时，三种模型的计算结果相差很小，这说明制动盘对流场的影响随距离增大而减小，改进制动盘及其邻域内的源项设置很难影响到远尾流区。

结果比较差异最大的位置出现在x-x_c＝6D处，主要表现为尾流的扩散宽度不足，速度在径向上很快恢复到环境水平。可能的原因是，此处远离制动盘影响，尾流模型计算出的径向湍流强度高于观测值(如图7，r/D>0.8时)，导致能量交换偏大，加快速度在径向上的恢复作用。当入流风速变大时，推力系数减小，近尾流区相对速度的降低值也相应减小，如在x-x_c＝2.5D处，U₀＝8.54m/s时，轮毂中心的相对速度为0.492，而当U₀＝9.56m/s，此值变为0.530；在远尾流区内，由于距离风轮远，风轮对流场的影响变小，加上周围环境作用，导致此处相对速度变化缓慢，受入流风速量值的影响减弱。

Dawin风力机尾流数值预测结果和实验值对比如图8(a)、图8(b)、图8(c)所示。对比结果表明除x-x_c＝4.15D轴向位置外，改进模型都能很好的与实验值匹配。此结果与Nibe-B 630kW在不同入流风速下的尾流相似，由此可见改进的尾流模型不仅能适用不同工况，还对风机类型不敏感，具有普适性。

综上所述，本发明考虑到制动盘动量源项S_u、机舱附加轴向动量源项S_d，湍流动能源项S_k、湍流动能耗散率源项S_ε对风力机尾流计算的影响，修正时均雷诺N-S方程以及k-ε方程后组成风力机尾流求解模型，在远尾流区域内自适应源项模型、固定源项模型的径向分布相同，改进型模型则进一步改进了近尾流区域的径向分布，提高了模型的计算精度。

Claims (5)

1.风力机尾流计算方法，其特征在于：

步骤1，利用制动盘模拟风力机，建立包含制动盘动量源项S_u、机舱附加轴向动量源项S_d、湍流动能源项S_k、湍流动能耗散率源项S_ε的尾流求解模型：

$\left\{\begin{array}{c}\frac{\&part;u_j}{\&part;x_j}=0\\ \frac{\&part;}{\&part;x_j}\left({\mathrm{\&rho;u}}_i{u}_j\right)=-\frac{\&part;p}{\&part;x_i}+\frac{\&part;}{\&part;x_j}(2{\mathrm{\&mu;S}}_{ij}-\mathrm{\&rho;}\stackrel{\&OverBar;}{u_i^{\&prime;}{u}_j^{\&prime;}})+S_u+S_d\\ \frac{\&part;}{\&part;x_i}\left({\mathrm{\&rho;ku}}_i\right)=\frac{\&part;}{\&part;x_j}\&lsqb;(\mathrm{\&mu;}+\frac{{\mathrm{\&mu;}}_t}{{\mathrm{\&sigma;}}_k})\frac{\&part;k}{\&part;x_j}\&rsqb;+P_k-\mathrm{\&rho;}\mathrm{\&epsiv;}+S_k\\ \frac{\&part;}{\&part;x_i}\left({\mathrm{\&rho;\&epsiv;u}}_i\right)=\frac{\&part;}{\&part;x_j}\&lsqb;(\mathrm{\&mu;}+\frac{{\mathrm{\&mu;}}_t}{{\mathrm{\&sigma;}}_{\mathrm{\&epsiv;}}})\frac{\&part;\mathrm{\&epsiv;}}{\&part;x_j}\&rsqb;+C_{1\mathrm{\&epsiv;}}\frac{\mathrm{\&epsiv;}}{k}{P}_k-C_{2\mathrm{\&epsiv;}}\mathrm{\&rho;}\frac{{\mathrm{\&epsiv;}}^2}{k}+S_{\mathrm{\&epsiv;}}\\ S_u=\frac{1}{2}{\mathrm{\&rho;C}}_x{u}^2\\ S_{\mathrm{\&epsiv;}}=C_{4\mathrm{\&epsiv;}}\frac{P_k^2}{\mathrm{\&rho;}k}\end{array}\right.,$

其中，ρ为流体密度，u_i、u_j为平均速度在各坐标轴的分量，p为压强，S_ij为平均应变张量， 为速度在各坐标轴的波动分量，μ为流体粘性系数，μ_t为湍流粘度系数，k为湍流动能，P_k为湍流动量生成率，ε为湍流动能耗散率，C_1ε、C_2ε、C_μ、C_x、C_4ε、σ_k、σ_ε为常数；

步骤2，设置中性大气边界层条件：速度入口u、湍流动能k、湍流动能耗散率ε、地面粗糙长度K_S；

步骤3，设置CFD求解参数；

步骤4，按中性大气边界层确定常数C_1ε,C_2ε,C_μ,σ_k,σ_ε的取值，开始迭代求解过程确定尾流区风速、湍流度的分布值。

2.根据权利要求1所述的风力机尾流计算方法，其特征在于，考虑源项的径向分布参数η_r，进一步修正制动盘动量源项S_u、湍流动能耗散率源项S_ε：

$S_u=\frac{1}{2}{\mathrm{\&rho;C}}_x{u}^2\&CenterDot;{\mathrm{\&eta;}}_r,$

$S_{\mathrm{\&epsiv;}}=\&lsqb;\left(1-\tilde{r}\right){\mathrm{\&eta;}}_{\mathrm{\&epsiv;}0}+\tilde{r}{\mathrm{\&eta;}}_{\mathrm{\&epsiv;}1}\&rsqb;\frac{P_k^2}{\mathrm{\&rho;}k},$

η_ε0、η_ε1为常数，η_ε0＝0.0，η_ε1＝0.9，为网格中心到制动盘径向距离r对制动盘半径D的归一量，

3.根据权利要求1或2所述的风力机尾流计算方法，其特征在于，步骤3中设置CFD求解参数时，湍流粘性比上限为2×10⁷。

4.根据权利要求1或2所述的风力机尾流计算方法，其特征在于，步骤2中的速度入口u、湍流动能k、湍流动能耗散率ε由如下表达式确定：

$\left\{\begin{array}{c}u=\frac{u_{*}}{\mathrm{\&kappa;}}\ln \left(\frac{z}{z_0}\right)\\ k=\frac{u_{*}^2}{\sqrt{C_{\mathrm{\&mu;}}}}\\ \mathrm{\&epsiv;}=\frac{u_{*}^2}{\mathrm{\&kappa;}z}\\ k=\mathrm{\&lambda;}{\left(U_0{I}_0\right)}^2\\ u_{*}={\mathrm{\&lambda;}}^{0.5}{C}_{\mathrm{\&mu;}}^{0.25}{U}_0{I}_0\\ z_0=H/\exp ({\mathrm{\&kappa;U}}_0/u_{*})\end{array}\right.,$

其中，u_*为摩擦速度，z为离地高度，z₀为大气粗糙度长度，表示近地层风速向下递减到零时的高度，κ为冯卡门常数，κ＝0.4187，λ为模型常数，对于各向同性湍流，λ取值为1.5，H为轮毂高度，U₀、I₀分别为轮毂高度H处的风速、湍流强度。

5.根据权利要求4所述的风力机尾流计算方法，其特征在于，步骤2中所述地面粗糙长度K_S的表达式为：E＝9.793，C_S＝1.0。