CN111475906A - 一种风力机尾流风速的算法 - Google Patents

Abstract

本发明为一种风力机尾流风速的算法，包括：步骤一：对入口边界条件进行初始赋值，设经过风力机的入流风速为u₀、湍流动能为k、耗散率为；步骤二：采用致动盘模型下的带叶轮动量源项和机舱的动量源项的不可压N‑S方程；步骤三：采用带有湍流动能源项和耗散率源项的湍流模型；步骤四：结合上述三个步骤进行求解风速。本发明解决单独采用致动盘方法计算尾流时，尾流恢复过快的问题，以及现有的修正方法对计算精度提升幅度不大的问题，在原有湍流动能源项的基础上，添加耗散率源项，使尾流中湍流动能的生成和耗散规律更加符合实际，由此提升尾流计算精度。

Description

一种风力机尾流风速的算法

技术领域

本发明涉及风力机尾流计算方法，具体涉及一种对k-ε湍流模型修订的风力机尾流计算方法，本发明属于风力发电领域。

背景技术

随着国家环保意识的提高，清洁能源发展对促进国民经济的进步和保护环境来说至关重要，而风能作为一种重要的清洁能源，在未来将会占据重要的地位。在风力发电的过程中，风力机尾流效应的研究对于风电场机组布局，风力机功率数据的采集有着重要的指导意义，近年来引起专家学者们的密切关注。

对于风力机尾流的数值研究，一般来说，计算精度主要依赖于两部分：一是对风轮本身的模拟，也就是说采用合适的模型模拟风力机的存在及其对周围大气产生的影响；二是对大气边界层的模拟，它是指选取适当的湍流模型预测风场内大气湍流的流动状况。

对于风轮的模拟，目前的方法主要采用致动系列方法，所述致动盘模型将风轮简化为一个圆盘，与风轮同轴。在实际的建模中，该圆盘并不真实存在，实质上致动盘是一个可穿透的圆柱形薄盘。根据动量理论，致动盘与风轮同轴，厚度为Δx，直径与风轮相当设为D，假设入流风速为u₀,则致动盘处风速为：u_D＝(1-a)u₀，式中：a为诱导因子，其值可由推力系数得到，即：

于是，风轮单位体积的动量源项为：

ρ为空气密度，C_x＝4a/(1-a)为对应的致动盘的拖曳系数。

由于空气无法从机舱区域流过，因此可以认为属于固体绕流，参照致动盘假设，机舱可定义为：

式中，C_D为机舱的拖曳系数。

对大气边界层的模拟，关于大气边界层的模拟，主要是指对大气湍流的模拟。对计算条件要求不高的雷诺平均方程即N-S方程加k-ε湍流模型的方法得到了工程界的普遍认可和广泛使用。尽管N-S方法有着诸多不可比拟的优势，但在解决某些复杂流动问题时存在一些不足，如低估尾流速度的亏损，加快尾流的恢复以及缩小尾流半径等。对此，学者们提出了各种修正模型用于提高风力机尾流的模拟精度。国际上，Kasmin和Masson建议在k-ε模型方程即湍流模型方程中添加一个表示从大尺度涡到小尺度涡能量传递关系的源项。将两部分算法结合起来，计算出相应的尾流速度场，El Kasmi A较早提出通过添加耗散率源项的方法，使湍流动能生成与其耗散相互协调，修正后计算精度有所提升，但是幅度不大，需要继续改进。

发明内容

为了解决单独采用致动盘方法计算尾流时，尾流恢复过快的问题，以及现有的修正方法对计算精度提升幅度不大的问题，在原有湍流动能源项的基础上，添加耗散率源项，使尾流中湍流动能的生成和耗散规律更加符合实际，由此提升尾流计算精度，本发明提供一种风力机尾流风速的算法，包括如下步骤：

步骤一：对入口边界条件进行初始赋值，设经过风力机的入流风速为u₀、湍流动能为k、耗散率为；

步骤二：采用致动盘模型模拟风轮影响，得到带源项的不可压N-S方程，所述N-S方程带有的源项为叶轮动量源项S_u和机舱的动量源项S_d；所述带源项的不可压N-S方程为两方程：

u_i和u_j为坐标方向的速度，x_i和x_j为坐标方向，ρ为密度，p为压力，μ为粘性系数，u′_i和u′_j为湍流脉动速度，

为雷诺应力，

所述u_D为致动盘处风速，所述C_X为致动盘的拖曳系数，所述Δx为致动盘厚度，

u_D＝(1-a)u₀ (7)，

C_x＝4a/(1-a) (8)，

所述a为诱导因子，其值可由推力系数得

所述C_D为机舱的拖曳系数。

步骤三：采用k-ε湍流模型，所述k-ε湍流模型包括三个方程：

所述k-ε湍流模型方程中，k为湍流动能；ε为耗散率；μ_t为湍流粘性系数；P_k为湍流动能生成项；σ_k、σ_ε分别为对应k和ε的普朗特常数；C_1ε、C_2ε和C_μ为常数，所述S_k为湍流动能源项，

所述β_p为时均动能转化为湍流动能的系数，β_d为湍流动能损失的系数，

所述B为时均动能转化为湍流动能的系数，所述S_ε为耗散率源项，

所述C_4ε为耗散率损失系数；

步骤四：结合步骤一和步骤三对步骤二所述的带源项的不可压N-S方程进行求解，得到风力机尾流速度场。

在上述方法中，方程(10)来源于文献A note on k-εmodelling of vegetationcanopy air-flows，由Sanz在2003年发表于Boundary-Layer Meteorology杂志上，该方程应用于树冠模型。当风流过树林时，由于树林的阻力，风速会显著降低，湍流动能也会有相应的变化，描述这种树林对风产生作用的模型就被称为树冠模型。树冠模型将树林看做一个阻力源项，方程(10)不仅考虑了来流时均动能向湍流动能的转化，也考虑湍流动能的不断损失。将两种流动现象进行类比时，发现当风流经风力机的风轮时，风轮对来流产生阻力并且吸收一部分能量，因此风速也会降低，同时由于风轮的扰动，流体中的湍流动能也会出现相应变化。因此，认为可以采用方程(10)作为湍流模型中，湍流动能方程的源项。

方程(12)最早来源于美国国家航空航天局的研究报告Computation ofturbulent flows using an extended k-e turbulence closure model，由Chen和Kim于1987年发表，目的为讨论湍流模型方程的闭合方式。此方程代表了大尺度湍流向小尺度湍流过度的输运过程，此过程由湍流动能的生成和耗散的时间尺度决定。当仔细调整参数时，方程(12)可以使得耗散率能够更快地响应平均流场，从而更有效地控制湍流动能的发展。之后的2008年，El Kasmi在Journal of wind engineering杂志上发表了题为An extendedk-εmodel for turbulent flow through horizontal-axis wind turbines的文章，将方程(12)应用在了风力机的模拟中。当风流经风力机的风轮时，快速旋转的风轮对流场产生较大的扰动，因此大尺度的湍流迅速破碎成小尺度的湍流，为了描述这种变化，认为可以采用方程(12)作为湍流模型中，耗散率方程的源项。

在上述方法中，动量方程中(2)的雷诺应力项，未知数为u′_iu′_j，一共是九个变量，矩阵形式为

因为这九个变量关于主对角线对称相等，即：u′₁u′₂＝u′₂u′₁、u′₁u′₃＝u′₃u′₁、u′₂u′₃＝u′₃u′₂，所以未知变量就有六个，即：u′₁u′₁、u′₂u′₂、u′₃u′₃、u′₁u′₂、u′₁u′₃、u′₂u′₃。

式中：下角标的1、2、3分别表示的是x、y、z三个方向；

u′_i、u′_j或者是u′₁、u′₂、u′₃都被称为为湍流脉动速度。

这六个未知量的时间平均值乘以密度ρ时，就被看做是(雷诺)应力，根据1877年Boussinesq(布辛尼斯克)的假设,湍流脉动所造成的附加应力可以与应变率关联起来，公式为：

p_t为湍流脉动速度所造成的压力；

为湍流动能；δ_i,j的取值规则为：当i＝j时为0，当i≠j时为1；

根据所有的公式进行求解，最终可以得到风力机尾流速度场。

本发明通过添加了湍流动能源项和耗散率动能源项，对参数β_p进行修正，使得其根据环境参数的情况处于一个合适的取值，最终能够使得计算出的风速更接近实际数值。

优选地，所述步骤二中，0.8≤C_D≤1.2。

为了进一步地精准的计算风速，所述步骤三中，对C_4ε进行进一步的修订C_4ε＝C₁*[2(r-0.5)²+0.2]，所述C₁为分布系数，通常取为1.0，r为致动盘上任意一点到致动盘轴线的距离。

优选地，所述r∈(0,D/2)，所述D为致动盘直径。

优选地，所述步骤三中，σ_k＝1.0，σ_ε＝1.3，C_μ＝0.033，C_1ε＝1.176，C_2ε＝1.92。

上述步骤仅能计算单台风机的尾流速度，为了进一步地计算多台风机的尾流风速，在所述步骤三中，B＝[10*(k_n/k₁)]^0.4，k₁、k_n分别为第一台、第n台风力机的入流风的湍流动能。

优选地，所述步骤三中，β_d＝1.0。

附图说明

图1为本发明专利步骤示意图

图2为U＝8m/s速度扩线

图3为U＝11m/s速度扩线

图4为Danwin风力机误差分析

图5为C_T＝0.607速度扩线

图6为C_T＝0.857速度扩线

图7为GH风力机误差分析

图8为串列多排风力机模拟结果

具体实施方式

为更容易理解本发明的优点、特征以及达到技术效果的技术方法，将参照例示性实施例进行更详细地描述，且本发明可以不同形式来实现，故不应被理解为本发明仅限于此处所陈述的实施例，相反地，对本领域的技术人员，所提供的实施例将更加透彻与全面且完整地传达本发明的范畴，且本发明将以申请专利文件的权利要求确定保护范围。

以下结合附图对本申请进行进一步的说明。

实施例一：Dawin风力机尾流模拟

Dawin风力机在瑞典境内，对比工况有两种是：

(1)U＝8m/s时，C_T＝0.82，湍流强度7％；

(2)U＝11m/s时，C_T＝0.65，湍流强度6％；

按照如下步骤建模，由计算机进行计算，

步骤一：对入口边界条件进行初始赋值，设经过风机的入流风速为u₀、湍流动能为k、耗散率为；

步骤二：采用致动盘模型模拟风轮影响，得到带源项的不可压N-S方程，所述N-S方程带有的源项为叶轮动量源项S_u，S_d为机舱的动量源项；所述带源项的不可压N-S方程为两方程：

方程(1)为

方程(2)为

u_i和u_j为坐标方向的速度，x_i和x_j为坐标方向，ρ为密度，p为压力，μ为粘性系数，u′_i和u′_j为湍流脉动速度，

为雷诺应力，

所述u_D为致动盘处风速，所述C_X为致动盘的拖曳系数，所述Δx为致动盘厚度，u_D＝(1-a)u₀，C_x＝4a/(1-a)，所述a为诱导因子；

所述C_D为机舱的拖曳系数为1.0。

步骤三：采用k-ε湍流模型，所述k-ε湍流模型包括三个方程：

方程(3)为

方程(4)为

方程(5)为

所述k-ε湍流模型方程中，k为湍流动能；ε为耗散率；μ_t为湍流粘性系数；P_k为湍流动能生成项；σ_k、σ_ε分别为对应k和ε的普朗特常数；C_1ε、C_2ε和C_μ为常数，σ_k＝1.0，σ_ε＝1.3，C_μ＝0.033，C_1ε＝1.176，C_2ε＝1.92；所述S_k为湍流动能源项，

所述β_p为时均动能转化为湍流动能的系数，β_d为湍流动能损失的系数，所述S_ε为耗散率源项，

所述C_4ε为耗散率损失系数，所述

所述B为时均动能转化为湍流动能的系数，B＝[10*(k_n/k₁)]^0.4，k₁、k_n分别为第一台、第n台风力机的入流风的湍流动能。

步骤四：结合步骤一和步骤三对步骤二所述的带源项的不可压N-S方程进行求解，得到风力机尾流速度场。

从图2、图3可以看出，单独采用致动盘的标准k-epsilon(即k-ε)模型计算结果与实验值有较大偏差，说明对标准k-epsilon模型进行修正是十分必要的；El Kasmi A的修正模型提升了计算精度，但仍可继续改进。可以看出，在两种工况下，修正模型的模拟结果较为准确。特别是风轮下游1D位置，此处处于近尾流区，流场扰动强烈，从1D截面位置处较为分散的实测数据也可以看出，但模拟结果依然能够反映出正确的尾流分布。从图4的误差分析中也可以看出，修正模型的模拟结果较为准确。

实施例二：Garrad Hassan实验数据

GH公司在大气边界层风洞，对缩比风力机进行了一系列的尾流实验。原型风力机的风轮直径为43.2m，轮毂高度50m。报告中给出了风轮下游2.5D、5D、7.5D、10D位置处的风速。对比工况有两种，分别是：

(1)U＝5.3m/s时，C_T＝0.607，粗糙度z₀＝0.075；

(2)U＝5.3m/s时，C_T＝0.858，粗糙度z₀＝0.075。

按照如下步骤建模，由计算机进行计算，

步骤一：对入口边界条件进行初始赋值，设经过风力机的入流风速为u₀、湍流动能为k、耗散率为；

步骤二：采用致动盘模型模拟风轮影响，得到带源项的不可压N-S方程，所述N-S方程带有的源项为叶轮动量源项S_u，S_d为机舱的动量源项；所述带源项的不可压N-S方程为两方程：

方程(1)为

方程(2)为

u_i和u_j为坐标方向的速度，x_i和x_j为坐标方向，ρ为密度，p为压力，μ为粘性系数，u′_i和u′_j为湍流脉动速度，

为雷诺应力，

所述u_D为致动盘处风速，所述C_X为致动盘的拖曳系数，所述Δx为致动盘厚度，u_D＝(1-a)u₀，C_x＝4a/(1-a)，所述a为诱导因子；

所述C_D为机舱的拖曳系数为1.0。

步骤三：采用k-ε湍流模型，所述k-ε湍流模型包括三个方程：

方程(3)为

方程(4)为

方程(5)为

所述k-ε湍流模型方程中，k为湍流动能；ε为耗散率；μ_t为湍流粘性系数；P_k为湍流动能生成项；σ_k、σ_ε分别为对应k和ε的普朗特常数；C_1ε、C_2ε和C_μ为常数，σ_k＝1.0，σ_ε＝1.3，C_μ＝0.033，C_1ε＝1.176，C_2ε＝1.92；所述S_k为湍流动能，

所述β_p为时均动能转化为湍流动能的系数，β_d湍流动能损失的系数，所述S_ε为耗散率源项，

所述C_4ε为耗散率损失系数，所述

所述B为时均动能转化为湍流动能的系数，B＝[10*(k_n/k₁)]^0.4，k₁、k_n分别为第一台、第n台风力机的入流风的湍流动能。

步骤四：结合步骤一和步骤三对步骤二所述的带源项的不可压N-S方程进行求解，得到风力机尾流速度场。

图5、图6可以看出，单独采用致动盘的标准k-epsilon模型计算结果与实验值有较大偏差，说明对标准k-epsilon模型进行修正是十分必要的；Park模型的精度有一定的提升，但仍远低于修正模型的精度，采用本发明的模拟结果依然能够反映出正确的尾流分布。

实施例三：Horns Rev海上风电场测量

选取Horns Rev海上风电场测量数据进行对比分析，场内有80台Vestas-V80风力机，额定容量2MW，轮毂高度70m，叶轮直径80m，根据风电场的实际运行情况，统计得到表3所示三组数据。

按照如下步骤建模，由计算机进行计算，

步骤一：对入口边界条件进行初始赋值，设经过风机的入流风速为u₀、湍流动能为k、耗散率为；

步骤二：采用致动盘模型模拟风轮影响，得到带源项的不可压N-S方程，所述N-S方程带有的源项为叶轮动量源项S_u，S_d为机舱的动量源项；所述带源项的不可压N-S方程为两方程：

方程(1)为

方程(2)为

u_i和u_j为坐标方向的速度，x_i和x_j为坐标方向，ρ为密度，p为压力，μ为粘性系数，u′_i和u′_j为湍流脉动速度，

为雷诺应力，

所述u_D为致动盘处风速，所述C_X为致动盘的拖曳系数，所述Δx为致动盘厚度，u_D＝(1-a)u₀，C_x＝4a/(1-a)，所述a为诱导因子；

所述C_D为机舱的拖曳系数为1.0。

步骤三：采用k-ε湍流模型，所述k-ε湍流模型包括三个方程：

方程(3)为

方程(4)为

方程(5)为

所述k-ε湍流模型方程中，k为湍流动能；ε为耗散率；μ_t为湍流粘性系数；P_k为湍流动能生成项；σ_k、σ_ε分别为对应k和ε的普朗特常数；C_1ε、C_2ε和C_μ为常数，σ_k＝1.0，σ_ε＝1.3，C_μ＝0.033，C_1ε＝1.176，C_2ε＝1.92；所述S_k为湍流动能，

所述β_p为时均动能转化为湍流动能的系数，β_d湍流动能损失的系数，所述S_ε为耗散率源项，

所述C_4ε为耗散率损失系数，所述

所述B为时均动能转化为湍流动能的系数，B＝[10*(k_n/k₁)]^0.4，k₁、k_n分别为第一台、第n台风力机的入流风的湍流动能。

步骤四：结合步骤一和步骤三对步骤二所述的带源项的不可压N-S方程进行求解，得到风力机尾流速度场。

从图8可以看出，在270°方向上，除第二台风力机有5％的相对功率偏差外，其余结果均与测量值符合较好，并且捕捉到了第6～8台风力机相对功率的微弱回升趋势；后7台风力机平均偏差为1.79％；随着风力机间距的增大，在222°风向上，计算结果与测量值仍然较为一致，后4台风力机平均偏差为1.66％；当风向为312°时，间距进一步增大为10.5D，此时误差增加，但后4台风力机平均偏差仍然较低，为2.97％；因此可见，修正模型可以准确模拟多排风力机情况下的风电场功率。

Claims (7)

1.一种风力机尾流风速的算法，包括如下步骤：

步骤一：对入口边界条件进行初始赋值，设经过风力机的入流风速为u₀、湍流动能为k、耗散率为；

步骤二：采用致动盘模型模拟风轮影响，得到带源项的不可压N-S方程，所述N-S方程带有的源项为叶轮动量源项S_u和机舱的动量源项S_d；所述带源项的不可压N-S方程为两方程：

方程(1)为

方程(2)为

u_i和u_j为坐标方向的速度，x_i和x_j为坐标方向，ρ为密度，p为压力，μ为粘性系数，u′_i和u′_j为湍流脉动速度，

为雷诺应力，

所述u_D为致动盘处风速，所述C_X为致动盘的拖曳系数，所述Δx为致动盘厚度，u_D＝(1-a)u₀，C_x＝4a/(1-a)，所述a为诱导因子；

所述C_D为机舱的拖曳系数。

步骤三：采用k-ε湍流模型，所述k-ε湍流模型包括三个方程：

方程(3)为

方程(4)为

方程(5)为

所述k-ε湍流模型方程中，k为湍流动能；ε为耗散率；μ_t为湍流粘性系数；P_k为湍流动能生成项；σ_k、σ_ε分别为对应k和ε的普朗特常数；C_1ε、C_2ε和C_μ为常数，

所述S_k为湍流动能源项，

所述β_p为时均动能转化为湍流动能的系数，β_d湍流动能损失的系数，

所述S_ε为耗散率源项，

所述C_4ε为耗散率损失系数，所述

所述B为时均动能转化为湍流动能的系数。

步骤四：结合步骤一和步骤三对步骤二所述的带源项的不可压N-S方程进行求解，得到风力机尾流速度场。

2.根据权利要求1所述的一种风力机尾流风速的算法，其特征在于，所述步骤二中，0.8≤C_D≤1.2。

3.根据权利要求1所述的一种风力机尾流风速的算法，其特征在于，所述步骤三中，所述C_4ε＝C₁*[2(r-0.5)²+0.2]，所述C₁为分布系数，r为致动盘上任意一点到致动盘轴线的距离。

4.根据权利要求3所述的一种风力机尾流风速的算法，其特征在于：所述r∈(0,D/2)，所述D为致动盘直径。

5.根据权利要求1所述的一种风力机尾流风速的算法，其特征在于，所述步骤三中，σ_k＝1.0，σ_ε＝1.3，C_μ＝0.033，C_1ε＝1.176，C_2ε＝1.92。

6.根据权利要求1所述的一种风力机尾流风速的算法，其特征在于：所述步骤三中，B＝[10*(k_n/k₁)]^0.4，k₁、k_n分别为第一台、第n台风力机的入流风的湍流动能。

7.根据权利要求1所述的一种风力机尾流风速的算法，其特征在于，所述步骤三中，β_d＝1.0。

Images