CN102663251A - 基于计算流体力学模型的风电场功率物理预测方法 - Google Patents

Abstract

本发明公开了风电场预测预报技术领域中的一种基于计算流体力学模型的风电场功率物理预测方法。包括：建立计算流体力学模型；对风电场的风况进行离散化处理，将各离散风况作为边界条件进行计算流体力学模型数值模拟，获得各离散风况的风电场空间流场分布；建立各离散风况下各台风电机组轮毂高度风速、风向及发电功率数据库；以数值天气预报参数为输入数据，利用数据库求得各台风电机组的风速和风向，进而计算出各台风电机组的发电功率，从而得到风电场功率的预测值。本发明适用各种风电场，功率预测阶段计算量小且计算时间短。

Description

基于计算流体力学模型的风电场功率物理预测方法

技术领域

本发明属于风电场功率预测预报技术领域，尤其涉及一种基于计算流体力学模型的风电场功率物理预测方法。

背景技术

我国风力发电经历了跨越式发展，未来依然将保持较快发展势头，十二五期间我国风电装机预计每年新增1000万千瓦以上。但是，风力发电固有的间歇性、随机波动性严重威胁电力系统运行的经济性、稳定性和供电质量。风电场功率预测是减轻风电对电网造成的不利影响、提高电网中风电装机比例的一种有效途径。国家对风电场功率预测已提出强制性要求，《风电场预测预报管理暂行办法》规定：新建风电场要同步建设风电预测预报体系，各风电场预测预报系统从2012年7月1日起正式开始运行。

根据是否需要风电场历史数据的支持，风电场功率预测可以分为两大类方法：一类是基于历史数据的统计模型方法，一类是不需要历史数据的物理模型方法。统计模型是在若干个历史数据和未来风电场功率之间建立一种映射关系，由历史数据或其他输入参数(如数值天气预报数据)外推得到未来一段时间的风电场功率。近年来，国内外统计模型方法的相关研究较多，常用方法有时间序列法、持续性算法、神经网络及支持向量基等。国内现有的预测系统均属于统计模型方法，但是，这种方法需要大量历史数据的支持，且对历史数据变化规律的一致性有很高的要求，否则将使功率预测的误差增大，甚至导致样本训练失败。因此，统计模型无法用于缺少历史记录的大量新建风电场。

基于物理模型的风电场功率预测方法不需要风电场历史功率数据的支持，它采用气象部门提供的数值天气预报(numerical weatherprediction--NWP)数据，通过考虑地形条件等因素，建立空气动力学模型描述流场，求取风电机组轮毂高度处的风速和风向，然后根据风电机组功率曲线进行风电场功率预测。因此，它既适用于新建风电场、也适用于已有风电场的功率预测。风电场功率物理预测方法最关键的环节在于：将NWP给出的风速、风向等参数准确地计算为风电机组轮毂高度水平的风速和风向，即NWP数据的精细化处理。它是风电功率预测的基础，也直接决定了风电功率预测的精度。对NWP数据的精细化处理主要采用以下两种方法：(1)采用诊断模型与解析算法分析风电场局地效应对流场的影响，即在假设的基础上解析求解简化的大气运动方程，这种方法计算量小，但是结果精度不高。(2)中尺度气象模式结合计算流体力学(Computational Fluid Dynamics-CFD)模型，动态模拟流场在风电场内的发展变化过程。这种方法可以获得比解析法更准确的流场分布，从而提高风电场功率预测精度。

计算流体力学模型(CFD模型)已经成功地应用于换热器内流动、汽车及飞行器的绕流等多种流体工程中，在风资源的评估中也有人利用CFD模型做了初步的尝试。但是，当它用于风电场功率预测时，每次预测都需要求解Navier-Stokes方程(N-S方程)进行流场计算，计算量巨大，难以满足预测时效的要求。因此，尽管CFD模型可以提高预测的精度，到目前为止还没有成功用于风电场功率预测领域的例子。

为此，本发明提出一种基于CFD模拟数据库的风电场功率预测新方法：它将风电场可能出现的风况离散化，以离散化的风况为依据建立“风况-各台风电机组轮毂高度风速风向及发电功率数据库”；然后以中尺度NWP数据为输入数据，利用该数据库完成风电场功率预测工作。这种方法有如下创新之处：

(1)本预测方法采用微尺度CFD模型求解N-S方程计算风速和风向分布，风速预测精度高；同时将耗时的CFD数值计算工作放到预测之前完成，预测过程计算量小，有效地解决了CFD方法的预测时效问题。

(2)采用单台风电机组功率预测、运行风电机组功率叠加获得风电场发电功率的预测原理，有效地解决了风电机组可利用率问题。风电机组可用率是影响风电场输出功率的重要因素之一，目前国内外的预测系统都未考虑风电机组可用率的问题。所谓风电机组可利用率，是指由于限电、故障等原因，每个风场内运行的风电机组数量是不断变化的，往往有部分风电机组处于停机状态，不能同时全部运行。随着风电机组的装机容量不断增加，在某些地方的限电现象时有发生，尤其是在东北和内蒙等地区的采暖季，为保证供热，不得不采取限制风电机组出力的调度措施，严重时一个风电场会有一半甚至更多的机组被限制出力。因此，在进行风电场功率预测时，限电或故障造成的风电机组可利用率已经成为一个必须考虑的问题，否则会严重影响预测精度。

发明内容

本发明的目的在于，针对当前预测风电场功率的方法存在的问题，提供一种基于计算流体力学模型的风电场功率物理预测方法。

为实现上述目的，本发明提供的技术方案是，一种基于计算流体力学模型的风电场功率物理预测方法，其特征是所述方法包括：

步骤1：建立计算流体力学模型；

步骤2：对风电场的风况进行离散化处理，将各离散风况作为边界条件进行计算流体力学模型数值模拟，获得各离散风况的风电场空间流场分布；

步骤3：建立各离散风况下各台风电机组轮毂高度风速、风向及发电功率数据库；

步骤4：求得各台风电机组的风速和风向，进而计算出各台风电机组的发电功率，从而得到风电场功率的预测值。

所述建立计算流体力学模型具体是：

步骤101：收集和测绘风电场及其周边设定范围内的地形等高线数据及粗糙度数据；

步骤102：根据步骤101获得的数据，建立预测风电场及其周边设定范围内的地形等高线数据及粗糙度的物理模型，然后确定风电场上方的空气流场区域；

步骤103：根据风电场及其周边设定范围内的地形，对步骤102确定的空气流场区域进行网格划分；

步骤104：确定风速轮廓线形状，所述风速轮廓线形状为计算流体力学模型的入口速度边界条件；

步骤105：初始化空气流场区域，进行计算流体力学模型数值计算，得到空气流场区域的风速；

步骤106：将步骤105计算的空气流场区域的风速与实际测量的空气流场区域的风速进行比较，根据比较结果优化计算流体力学模型。

所述步骤2具体是：

步骤201：对风电场的风况进行离散化处理，得到离散风况的风速和风向；

步骤202：分别将各个离散风况的风速和风向代入风速轮廓线，作为空气流动区域入口速度边界条件，利用优化后的计算流体力学模型进行模拟计算，得到各个离散风况对应的空气流场区域的风速和风向；

所述步骤3具体是：

步骤301：确定各台风电机组在风电场中的空间坐标；

步骤302：在各个离散风况对应的空气流场区域中，读出各台风电机组轮毂高度周边网格的风速和风向，然后利用线性插值法求得各个离散风况对应的风电场中各台风电机组轮毂高度的风速和风向；

步骤303：针对每个离散风况对应的风电场中各台风电机组轮毂高度的风速，利用风电场空气密度调节风电场中各台风电机组轮毂高度的风速值，得到各台风电机组轮毂高度的风速调节值；

步骤304：根据各台风电机组之间的相对位置，利用尾流模型求解尾流效应引起的各台风电机组轮毂高度的风速衰减值；

步骤305：根据各台风电机组轮毂高度的风速调节值和各台风电机组轮毂高度的风速衰减值，计算各台风电机组轮毂高度的风速修正值；

步骤306：根据各台风电机组轮毂高度的风速修正值，按照风电场各台风电机组的功率曲线，计算出各台风电机组的发电功率，从而获得各个离散风况对应的风电场中各台风电机组发电功率；

步骤307：将各个离散风况对应的风电场中各台风电机组的风速、风向和发电功率存入数据库。

所述步骤4具体是：

步骤401：计算待预测时间段内的风电机组轮毂高度的平均空气密度，在所述数据库中查询所述平均空气密度对应的数据表；

步骤402：以中尺度数值天气预报的风速和风向时间序列参数作为输入数据，在步骤401查询到的数据表中，查询与设定时刻所述输入数据输入的风速和风向值最相近的四个风况，读出四个风况的各风电机组轮毂高度风速、风向及发电功率；

步骤403：利用分段线性插值法求出所述设定时刻各风电机组轮毂高度的风速及发电功率预测值；

步骤404：确定所述设定时刻正在运行的风电机组；

步骤405：将正在运行的风电机组的功率叠加，获得所述设定时刻整个风电场的功率预测值；

步骤406：改变设定时刻，重复上述步骤401-405，得到整个待测时间段的风电场的功率预测值。

所述利用风电场空气密度调节风电场中各台风电机组轮毂高度的风速值具体采用公式其中为第n台风电机组轮毂高度风速调节值，u_n为第n台风电机组轮毂高度的风速，ρ₀为标准空气密度且ρ₀＝1.225kg/m³，ρ_m为设定时间段内风电场轮毂高度的平均空气密度。

所述利用尾流模型求解尾流效应引起的各台风电机组轮毂高度的风速衰减值具体利用公式 $\mathrm{\Δu}=-\frac{u_{\mathrm{WTn}}^m}{9}{\left(C_T{\mathrm{Ax}}^{-2}\right)}^{\frac{1}{3}}{[R_w^{\frac{3}{2}}{\left({3c}_1^2{C}_T\mathrm{Ax}\right)}^{-\frac{1}{2}}-\left(\frac{{35}^{\frac{3}{10}}}{2\mathrm{\π}}{\left({3c}_1^2\right)}^{-\frac{1}{5}}\right)]}^2,$ 其中，Δu为风电机组轮毂高度的风速衰减值，

为第n台风电机组轮毂高度风速调节值，A为风力机扫风面积，C_T为风电机组推力系数，c₁为无量纲混合长， $R_w={\left[\frac{35}{2\mathrm{\π}}\right]}^{\frac{1}{5}}{\left[{{3c}_1}^2\right]}^{\frac{1}{5}}{\left[C_T\mathrm{Ax}\right]}^{\frac{1}{3}}.$

所述计算各台风电机组轮毂高度的风速修正值具体利用公式

其中，

为各台风电机组轮毂高度的风速修正值，

为各台风电机组轮毂高度的风速调节值，Δu为各台风电机组轮毂高度的风速衰减值。

所述计算出各台风电机组的发电功率具体利用公式

其中，

为第n台风电机组的发电功率，C_p为风能利用系数，ρ₀为标准空气密度且ρ₀＝1.225kg/m³，A为风力机扫风面积，

为第n台风电机组轮毂高度风速调节值。

本发明无需历史数据的支持，既适用新建风电场，也适用于已运行风电场；同时，本发明将耗时的计算流体力学模型数值计算工作放到预测之前完成，功率预测阶段计算量小且计算时间短。

附图说明

图1是基于计算流体力学模型的风电场功率物理预测方法流程图；

图2是基于计算流体力学模型的风电场功率物理预测方法结构图；

图3是根据各台风电机组轮毂高度风速、风向及发电功率数据库得到风电场功率的预测值的流程图；

图4是使用本预测方法得到的预测误差表；

图5是预测功率与实测功率对比图。

具体实施方式

下面结合附图，对优选实施例作详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。

实施例1

图1是基于计算流体力学模型的风电场功率物理预测方法流程图，图2是基于计算流体力学模型的风电场功率物理预测方法结构图。根据图1和图2所示，本发明提供的基于计算流体力学CFD模型的风电场功率预测方法包括：

步骤1：建立计算流体力学模型。

其具体实现过程如下：

步骤101：借助GIS地理信息系统等工具，收集和测绘风电场及其周边设定范围内的地形等高线数据及粗糙度数据。

步骤102：根据步骤101获得的数据，建立预测风电场及其周边设定范围内的地形等高线数据及粗糙度的物理模型，然后确定风电场上方的空气流场区域。

步骤103：根据风电场及其周边设定范围内的地形，对步骤102确定的空气流场区域进行网格划分。

首先，根据风电场区域的地形情况确定合理的网格划分方案，网格划分方案指根据风电场区域的地形，确定网格的形状、尺寸及数量。不同地形适用的方案不同，根据CFD工程计算经验确定网格方案。

之后，对步骤102建立的空气流场区域划分计算网格。为减小网格数量，以六面体网格为主，风电场区域水平方向网格在十米的数量级上，周边区域网格更大。高度方向保证200米以内有十层以上网格。

步骤104：确定风速轮廓线形状，所述风速轮廓线形状为计算流体力学模型的入口速度边界条件。

近地层中，风速随高度显著变化，风速随高度的变化曲线称为风速轮廓线。假定大气层结为中性，风速轮廓线可按照幂定律公式计算：

式中，u_n为Z_n高度处风速，u₁为Z₁高度处风速，α为风切变指数。

风速轮廓线的形状取决于α值，α值与地形和地面粗糙度有关，通常可由两种方法确定：(1)根据风电场当地测风塔不同高度的测风数据拟合出α的数值；(2)根据地面粗糙度等级及粗糙度长度确定α经验值。

步骤105：初始化空气流场区域，进行计算流体力学CFD模型数值计算，得到空气流场区域的风速和风向分布。

进行计算流体力学CFD模型数值计算时，首先设定边界条件及湍流模型。边界条件设置为计算域顶部及侧边界为压力边界条件，设定为大气压力，入口处为速度边界条件，设定为入口风速轮廓线。

其次，根据风电场地形特点的不同，湍流模型通常采用标准k-ε模型或RNG k-ε双方程模型。迭代收敛后，可获得该计算工况风电场空间区域空气流场分布。通常，CFD数值模拟工作可以利用商用流体力学计算软件完成。

对划分的所有网格点赋给N-S方程待求解参数的迭代初始值，即是流场的初始化，理论上初始值只会影响收敛的速度，不会影响最终的计算结果。

由于风电场区域的空气流动为典型的湍流流动，由Navier-Stokes方程描述。将其时均化后，采用标准k-ε模型或RNG k-ε双方程模型完成方程的闭环，具体为：

$\frac{\∂\mathrm{\ρ}}{\∂t}+\frac{\∂\left(\mathrm{\ρ}{\stackrel{\‾}{u}}_i\right)}{{\∂x}_i}=0$

$\frac{\∂\left(\mathrm{\ρ}{\stackrel{\‾}{u}}_i\right)}{\∂t}+\frac{\∂\left(\mathrm{\ρ}{\stackrel{\‾}{u}}_i{\stackrel{\‾}{u}}_j\right)}{{\∂x}_j}=-\frac{\∂\stackrel{\‾}{p}}{{\∂x}_i}+\frac{\∂}{{\∂x}_j}({\stackrel{\‾}{\mathrm{\τ}}}_{\mathrm{ij}}-\mathrm{\ρ}\stackrel{\‾}{{u_i}^{\′}{u_j}^{\′}})+F_i$

式中：ρ为空气密度，t为时间，P为流体的静压力，u_i为速度在某一方向上的分量，x_i为坐标分量，F_i为i方向的体积力，τ_ij为粘性应力张量，可以用下式表示：

${\stackrel{\‾}{\mathrm{\τ}}}_{\mathrm{ij}}={2\mathrm{\μS}}_{\mathrm{ij}}-\frac{2}{3}\mathrm{\μ}\frac{\∂{\stackrel{\‾}{u}}_k}{{\∂x}_k}{\mathrm{\δ}}_{\mathrm{ij}}$

上式中μ为层流粘性系数，变形率张量

${\mathrm{\δ}}_{\mathrm{ij}}=\left\{\begin{array}{c}0,i\≠j\\ 1,i=j\end{array}\right.,$ 为kerodiler函数。

雷诺应力可用下式表示：

$-\mathrm{\ρ}\stackrel{\‾}{{u_i}^{\′}{u_j}^{\′}}={\mathrm{\μ}}_t(\frac{\∂{\stackrel{\‾}{u}}_i}{{\∂x}_j}+\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i})-\frac{2}{3}(\mathrm{\ρk}+{\mathrm{\μ}}_t\frac{\∂{\stackrel{\‾}{u}}_i}{{\∂x}_i})$

上式中，k为单位质量流体湍流脉动动能：

μ_t为湍流粘性，由以下求解：

(1)标准k-ε模型

$\mathrm{\ϵ}=\frac{\mathrm{\μ}}{\mathrm{\ρ}}\stackrel{\‾}{\left(\frac{{{\∂\mathrm{\μ}}_i}^{\′}}{{\∂x}_k}\right)\left(\frac{{{\∂\mathrm{\μ}}_i}^{\′}}{{\∂x}_k}\right)}$

引入ε后，k方程可表示为：

$\frac{\∂\left(\mathrm{\ρk}\right)}{\∂t}+{\mathrm{\ρu}}_j\frac{\∂k}{{\∂x}_j}=\frac{\∂}{{\∂x}_j}\left[(\mathrm{\μ}+\frac{{\mathrm{\μ}}_t}{{\mathrm{\σ}}_k})\frac{\∂k}{{\∂x}_j}\right]+{\mathrm{\μ}}_t\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i}(\frac{\∂{\stackrel{\‾}{u}}_i}{{\∂x}_j}+\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i})-\mathrm{\ρ\ϵ}$

确定耗散率ε的微分方程是：

$\frac{\∂\left(\mathrm{\ρ\ϵ}\right)}{\∂t}+{\mathrm{\ρu}}_j\frac{\∂\mathrm{\ϵ}}{{\∂x}_j}=\frac{\∂}{{\∂x}_j}\left[(\mathrm{\μ}+\frac{{\mathrm{\μ}}_t}{{\mathrm{\σ}}_{\mathrm{\ϵ}}})\frac{\∂\mathrm{\ϵ}}{{\∂x}_j}\right]+\frac{c_1\mathrm{\ϵ}}{k}{\mathrm{\μ}}_t\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i}(\frac{\∂{\stackrel{\‾}{u}}_i}{{\∂x}_j}+\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i})-c_2\mathrm{\ρ}\frac{{\mathrm{\ϵ}}^2}{k}$

相应地，μ_t的表达式写成：

${\mathrm{\μ}}_t=c_{\mathrm{\μ}}\mathrm{\ρ}\frac{k^2}{\mathrm{\ϵ}}$

上式中，σ_ε为ε的湍流Prandlt数，一般取1.3，c₁、c₂和c_u为经验常数，分别取1.44、1.92和0.09。

(2)RNG k-ε模型

RNG k-ε模型中，k方程依然采用前述公式，在ε方程中添加了一个附加项，ε方程的形式是：

$\frac{\∂\left(\mathrm{\ρ\ϵ}\right)}{\∂t}+{\mathrm{\ρu}}_j\frac{\∂\mathrm{\ϵ}}{{\∂x}_j}=\frac{\∂}{{\∂x}_j}\left[(\mathrm{\μ}+\frac{{\mathrm{\μ}}_t}{{\mathrm{\σ}}_{\mathrm{\ϵ}}})\frac{\∂\mathrm{\ϵ}}{{\∂x}_j}\right]+\frac{c_1\mathrm{\ϵ}}{k}{\mathrm{\μ}}_t\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i}(\frac{\∂{\stackrel{\‾}{u}}_i}{{\∂x}_j}+\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i})-c_2\mathrm{\ρ}\frac{{\mathrm{\ϵ}}^2}{k}-R$

其中，

R为变形率引起的附加源项，无量纲参数η＝Sk/ε，η₀＝4.38，β＝0.012，

步骤106：将步骤105计算的空气流场区域的风速与实际测量的空气流场区域的风速进行比较，根据比较结果优化计算流体力学模型。

模拟风速与实测风速进行对比，是为了验证模拟计算方法是否正确。湍流模型有成熟的多种，k-ε模型、RNG k-ε双方程模型等等，根据经验以及模拟风速与实测风速进行对比的结果试算确定。

步骤2：对风电场的风况进行离散化处理，将各离散风况作为边界条件进行计算流体力学模型数值模拟，获得各离散风况的风电场空间流场分布。

步骤201：对风电场的风况进行离散化处理，得到离散风况的风速和风向。

根据风电机组的有效风速范围及风电场的风资源特性，对风电场可能出现的风速、风向等风况信息进行统计并分析，以对风电场功率影响较大的参数(如风速和风向)为依据对整个风况范围进行离散化，从中筛选出能够覆盖整个风况范围的若干典型风况，作为下一步骤GFD数值计算的边界条件。通常，选择的离散工况越多，预测精度越高，但是，相应的计算工作量也越大。

典型风况的较好的选取方法是：风速的选择为：切入风速和额定风速之间风速值间隔较小，大于额定风速后选择的风速值间隔可以较大；风向的选择为：16个均分的标准风向，分别为0°、22.5°、45°、67.5°、90°、112.5°、135°、157.5°、180°、202.5°、225°、247.5°、270°、292.5°、315°、337.5°。选取的每个风速值分别给定16个风向，每个风速和风向的组合为一个离散风况。例如：风电机组的切入风速4m/s，额定风速13m/s，切出风速25m/s。若选定6个离散风况的风速值，分别为4、6、8、10、12、14m/s；每个风速和风向的组合为一个离散风况，共计96个离散风况。

步骤202：分别将各个离散风况的风速和风向代入风速轮廓线，作为空气流动区域入口速度边界条件，利用优化后的计算流体力学模型进行模拟计算，得到各个离散风况对应的空气流场区域的风速和风向分布。

将步骤201确定的一个离散风况的风速值和风向，赋给步骤104中确定的风速轮廓线，作为风电场空气流动区域入口速度边界条件；利用步骤106中确定的GFD方案进行数值模拟计算，获得该风况相对应的风电场区域空气流场分布。针对每一个离散风况的风速值和风向，重复步骤202，分别以其他各个离散工况的风速和风向作为边界条件，进行GFD数值模拟，获得所有离散工况下风电场区域空气流场的风速和风向分布。

步骤3：建立各离散风况下各台风电机组轮毂高度风速、风向及发电功率数据库。

步骤301：确定各台风电机组在风电场中的空间坐标。

在风电场中对风电机组进行数字化定位，确定各台风电机组在风电场中的具体空间坐标及其相对于受尾流影响的相邻风机的相对位置。

步骤302：在各个离散风况对应的空气流场区域中，读出各台风电机组轮毂高度周边网格的风速和风向，然后利用线性插值法求得各个离散风况对应的风电场中各台风电机组轮毂高度的风速和风向。

CFD计算结束后，会获得各个网格点的风速值(三个方向)，不同的CFD商业软件有不同的数据存储格式，此时风速和风向都是已知量。

步骤303：针对每个离散风况对应的风电场中各台风电机组轮毂高度的风速，利用风电场空气密度调节风电场中各台风电机组轮毂高度的风速值，得到各台风电机组轮毂高度的风速调节值。

根据风电场每年空气密度的变化范围，选择若干密度值ρ_m，m＝1，2，3…，密度修正以修正风速的方式体现，进行密度修正后第n台风电机组轮毂高度风速度值为

式中，u_n为CFD模拟获得的第n台风电机组轮毂高度的风速值，ρ₀为标准空气密度ρ₀＝1.225kg/m³，ρ_m为某一时期风电场轮毂高度水平的平均空气密度。

例：某风电场风机轮毂高度全年的空气密度ρ在0.944-1.202kg/m³的范围内，可取ρ_m分别为0.95、1.0、1.05、1.1、1.15、1.2(m＝1，2，3，4，5，6)，针对每个密度值，调节风机轮毂高度的风速。

步骤304：根据各台风电机组之间的相对位置，利用尾流模型求解尾流效应引起的各台风电机组轮毂高度的风速衰减值。

可采用拉森解析模型(Larsen wake model)、杰森模型(Jensen wake model)等计算尾流效应。Larsen模型是基于普朗特边界层方程的渐近表达式，假定下风向不同位置的风速衰减具有相似性，并且风速只会发生中等程度的衰减，则下风向L＝x处的尾流影响区域半径R_w及Larsen尾流模型的风速衰减Δu表达式分别为：

$R_w={\left[\frac{35}{2\mathrm{\π}}\right]}^{\frac{1}{5}}{\left[{3c}_1^2\right]}^{\frac{1}{5}}{\left[C_T\mathrm{Ax}\right]}^{\frac{1}{3}}$

$\mathrm{\Δu}=-\frac{u_{\mathrm{WTn}}^m}{9}{\left(C_T{\mathrm{Ax}}^{-2}\right)}^{\frac{1}{3}}{[R_w^{\frac{3}{2}}{\left({3c}_1^2{C}_T\mathrm{Ax}\right)}^{-\frac{1}{2}}-\left(\frac{{35}^{\frac{3}{10}}}{2\mathrm{\π}}{\left({3c}_1^2\right)}^{-\frac{1}{5}}\right)]}^2$

上式中，

为风电机组轮毂高度的平均风速，A为风力机扫风面积，C_T为风电机组推力系数，c₁为无量纲混合长。

步骤305：根据各台风电机组轮毂高度的风速调节值和各台风电机组轮毂高度的风速衰减值，计算各台风电机组轮毂高度的风速修正值。

将风速衰减值Δu叠加到各风电机组位置的CFD模拟风速

上，即可获得该离散风况下计及密度修正和尾流效应后各台风电机组轮毂高度的修正风速值 $u_n^m=u_{\mathrm{WTn}}^m+\mathrm{\Δu}.$

步骤306：根据各台风电机组轮毂高度的风速修正值，按照风电场各台风电机组的功率曲线，计算出各台风电机组的发电功率，从而获得各个离散风况对应的风电场中各台风电机组发电功率。

根据步骤305获得的调节风速值

按照风电场各台风电机组的功率曲线，计算出各台风电机组的发电功率，从而获得该离散风况条件下各台风电机组发电功率数据库。空气密度ρ_m条件下第n台风电机组的发电功率为

则有

$P_n^m=\frac{1}{2}{C}_p{\mathrm{\ρ}}_0{\left(u_n^m\right)}^{3}A,$

上式中，C_p为风能利用系数；ρ₀为标准空气密度ρ₀＝1.225kg/m³；A为风轮扫风面积，单位为m²。

步骤307：将各个离散风况对应的风电场中各台风电机组的风速、风向和发电功率存入数据库。

数据库包括风电场设备信息表、风机运行信息表、风机CFD模拟数据表、历史数据记录表等。其中最重要的是风机CFD模拟数据表，字段有：空气密度、边界风速、边界风向、测风塔风向、测风塔风速、风电机组1风速、风向、发电功率、风电机组2风速、风向、发电功率、风电机组3风速、风向、发电功率等等。

步骤4：求得各台风电机组的风速和风向，进而计算出各台风电机组的发电功率，从而得到风电场功率的预测值。

图3是根据各台风电机组轮毂高度风速、风向及发电功率数据库得到风电场功率的预测值的流程图。图3中，计算出各台风电机组的发电功率，从而得到风电场功率的预测值的具体过程是：

步骤401：计算待预测时间段内的风电机组轮毂高度的平均空气密度，在所述数据库中查询所述平均空气密度对应的数据表。

首先根据NWP温度、压力和湿度等参数，计算得到待预测时间区间风电机组轮毂高度的平均空气密度，在步骤307建立的数据库中查询代表该空气密度的ρ_m所对应的数据表。

τ_i时刻空气密度： $\mathrm{\ρ}\left({\mathrm{\τ}}_i\right)=\frac{1.276}{1+0.00366t}\×\frac{(p-0.378e)}{1000}$

平均空气密度： $\stackrel{\‾}{\mathrm{\ρ}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\ρ}\left({\mathrm{\τ}}_i\right)$

式中：p-大气压力(hPa)；

t-温度(℃)；

e-水汽分压力(hPa)。

例：假定数据表中空气密度ρ_m＝0.95、1.0、1.05、1.1、1.15、1.2(m＝1，2，3，4，5，6)分别代表0.925-0.975、0.975-1.025、1.025-1.075、1.075-1.125、1.125-1.175、1.175-1.225的密度区间。假设计算得到待预测期间平均空气密度

则查询m＝4即ρ_m＝1.1kg/m³的数据表。

步骤402：以中尺度数值天气预报的风速和风向时间序列参数作为输入数据，在步骤401查询到的数据表中，查询与设定时刻所述输入数据输入的风速和风向值最相近的四个风况，读出四个风况的各风电机组轮毂高度风速、风向及发电功率。

步骤403：利用分段线性插值法求出所述设定时刻各风电机组轮毂高度的风速及发电功率预测值。

在密度ρ_m的条件下，第n台风电机组轮毂高度处风速由下式计算：

$u_n^m(V,d)=\frac{V-V_i}{V_{i+1}-V_i}(u_{i+1}-u_i)+u_i$

$u_i=\frac{d-d_j}{d_{j+1}-d_j}(u_{i,j+1}-u_{i,j})+u_{i,j}$

$u_{i+1}=\frac{d-d_j}{d_{j+1}-d_j}(u_{i+1,j+1}-u_{i+1,j})+u_{i+1,j}$

式中：

V为NWP设定时刻的风速；d为NWP某时刻的风向；

V_i，V_i+1为数据库中的两个相邻风速值(V_i≤V≤V_i+1)；

d_j，d_j+1为数据库中的两个相邻风向值(d_j≤d≤d_j+1)；

u_i，j为数据库中V_i，d_j对应的第n台风机轮毂高度水平风速；

u_i，j+1为数据库中V_i，d_j+1对应的第n台风机轮毂高度水平风速；

u_i+1，j为数据库中V_i+1，d_j对应的第n台风机轮毂高度水平风速；

u_i+1，j+1为数据库中V_i+1，d_j+1对应的第n台风机轮毂高度水平风速。

发电功率预测值可由两种方法获得：(1)可以在数据表中，利用线性插值法计算发电功率预测值，与上述风速的插值计算同理；(2)可以将上述计算得到的风速

代入步骤306中所述的发电功率计算公式获得。

步骤404：确定所述设定时刻正在运行的风电机组。

步骤405：将正在运行的风电机组的功率叠加，获得所述设定时刻整个风电场的功率预测值。

步骤406：改变设定时刻，重复上述步骤401-405，得到整个待测时间段的风电场的功率预测值。

实施例2

以中国北方某风电场为例进行验证分析，预测风电场所在区域海拔高度1200～1900m，风电场占地面积约45平方平公里，其中坐落着若干村庄，区域内存在大量农田，种植玉米等高秆作物。风电场装机容量183MW，风电机组额定容量1.5MW。

数据包括：风电场地形、粗糙度数据以及风电机组信息数据，用于建立CFD模型；从SCADA系统中提取的每15min的平均风电场功率实测值、平均风速，用于验证预测结果的精度；以及同期数值中尺度天气预报结果，用作预测的输入数据。预测时间为2010年。

预测后，利用下述公式进行精度评估。

$\mathrm{RMSE}=\frac{\sqrt{{\mathrm{\Σ}}_{i=1}^n{(P_{\mathrm{Mi}}-P_{\mathrm{Pi}})}^2}}{\mathrm{Cap}\×\sqrt{n}}$

$\mathrm{MAE}=\frac{{\mathrm{\Σ}}_{i=1}^n|P_{\mathrm{Mi}}-P_{\mathrm{Pi}}|}{\mathrm{Cap}\×n}$

式中：

RMSE为均方根误差；

MAE为平均绝对误差；

P_Mi为i时刻的实际功率；

P_Pi为i时刻的预测功率；

Cap为风电场开机(或风力发电机组)总容量；

n为预测点数。

采用本预测方法获得的预测结果年均方根误差为15.086％，年平均绝对误差10.693％，图4所示为本方法预测结果的各月均方根误差和平均绝对误差。由图4可见，整体精度较高，满足国家《风电场功率预测管理办法》的要求。图5所示为7月21日-27日一周时间中该风电场预测功率与实测功率的对比。由图可见，本预测方法在预测功率变化规律方面表现出色。

本发明采用中尺度数值天气预报结合CFD模型、建立了一种基于数据库的物理预测新方法，它具有以下重要优势：

1、本预测方法将风电场可能出现的风况条件离散化，以风速和风向两个风况参数为依据，建立流场及风电机组功率数据库，以数据库为基础完成风电场功率预测工作。它无需历史数据的支持，既适用新建风电场，也适用于已运行风电场。

2、对任何一个特定地形风的电场，这种风电功率预测方法将耗时的CFD数值计算工作放到预测之前完成，功率预测阶段计算量小、计算时间短、预测时效性强，是一种适用于工程应用的鲁棒性强的预测方法。

3、本预测方法利用CFD数值求解N-S方程获得流场风速分布，能够充分体现地形和粗糙度造成的湍流对风电机组功率的影响，预测精度高，尤其适用于复杂地形。

4、采用单台风电机组功率预测、运行风电机组功率叠加获得风电场发电功率的预测原理，有效地解决了由于限电、维修维护等原因造成的风电机组可利用率问题。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (9)

1.一种基于计算流体力学模型的风电场功率物理预测方法，其特征是所述方法包括：

步骤1：建立计算流体力学模型；

步骤2：对风电场的风况进行离散化处理，将各离散风况作为边界条件进行计算流体力学模型数值模拟，获得各离散风况的风电场空间流场分布；

步骤3：建立各离散风况下各台风电机组轮毂高度风速、风向及发电功率数据库；

步骤4：求得各台风电机组的风速和风向，进而计算出各台风电机组的发电功率，从而得到风电场功率的预测值。

2.根据权利要求1所述的物理预测方法，其特征是所述建立计算流体力学模型具体是：

步骤101：收集和测绘风电场及其周边设定范围内的地形等高线数据及粗糙度数据；

步骤102：根据步骤101获得的数据，建立预测风电场及其周边设定范围内的地形等高线数据及粗糙度的物理模型，然后确定风电场上方的空气流场区域；

步骤103：根据风电场及其周边设定范围内的地形，对步骤102确定的空气流场区域进行网格划分；

步骤104：确定风速轮廓线形状，所述风速轮廓线形状为计算流体力学模型的入口速度边界条件；

步骤105：初始化空气流场区域，进行计算流体力学模型数值计算，得到空气流场区域的风速；

步骤106：将步骤105计算的空气流场区域的风速与实际测量的空气流场区 域的风速进行比较，根据比较结果优化计算流体力学模型。

3.根据权利要求2所述的物理预测方法，其特征是所述步骤2具体是：

步骤201：对风电场的风况进行离散化处理，得到离散风况的风速和风向；

步骤202：分别将各个离散风况的风速和风向代入风速轮廓线，作为空气流动区域入口速度边界条件，利用优化后的计算流体力学模型进行模拟计算，得到各个离散风况对应的空气流场区域的风速和风向分布。

4.根据权利要求3所述的物理预测方法，其特征是所述步骤3具体是：

步骤301：确定各台风电机组在风电场中的空间坐标；

步骤302：在各个离散风况对应的空气流场区域中，读出各台风电机组轮毂高度周边网格的风速和风向，然后利用线性插值法求得各个离散风况对应的风电场中各台风电机组轮毂高度的风速和风向；

步骤303：针对每个离散风况对应的风电场中各台风电机组轮毂高度的风速，利用风电场空气密度调节风电场中各台风电机组轮毂高度的风速值，得到各台风电机组轮毂高度的风速调节值；

步骤304：根据各台风电机组之间的相对位置，利用尾流模型求解尾流效应引起的各台风电机组轮毂高度的风速衰减值；

步骤305：根据各台风电机组轮毂高度的风速调节值和各台风电机组轮毂高度的风速衰减值，计算各台风电机组轮毂高度的风速修正值；

步骤305：根据各台风电机组轮毂高度的风速修正值，按照风电场各台风电机组的功率曲线，计算出各台风电机组的发电功率，从而获得各个离散风况对应的风电场中各台风电机组发电功率；

步骤307：将各个离散风况对应的风电场中各台风电机组的风速、风向和发电功率存入数据库。

5.根据权利要求4所述的物理预测方法，其特征是所述步骤4具体是：

步骤401：计算待预测时间段内的风电机组轮毂高度的平均空气密度，在所 述数据库中查询所述平均空气密度对应的数据表；

步骤402：以中尺度数值天气预报的风速和风向时间序列参数作为输入数据，在步骤401查询到的数据表中，查询与设定时刻所述输入数据输入的风速和风向值最相近的四个风况，读出四个风况的各风电机组轮毂高度风速、风向及发电功率；

步骤403：利用分段线性插值法求出所述设定时刻各风电机组轮毂高度的风速及发电功率预测值；

步骤404：确定所述设定时刻正在运行的风电机组；

步骤405：将正在运行的风电机组的功率叠加，获得所述设定时刻整个风电场的功率预测值；

步骤406：改变设定时刻，重复上述步骤401-405，得到整个待测时间段的风电场的功率预测值。

6.根据权利要求4所述的物理预测方法，其特征是所述利用风电场空气密度调节风电场中各台风电机组轮毂高度的风速值具体采用公式 

其中 

为第n台风电机组轮毂高度风速调节值，u_n为第n台风电机组轮毂高度的风速，ρ₀为标准空气密度且ρ₀＝1.225kg/m³，ρ_m为设定时间段内风电场轮毂高度的平均空气密度。

7.根据权利要求6所述的物理预测方法，其特征是所述利用尾流模型求解尾流效应引起的各台风电机组轮毂高度的风速衰减值具体利用公式 

其中，Δu为风电机组轮毂高度的风速衰减值， 为第n台风电机组轮毂高度风速调节值，A为风力机扫风面积， C_T为风电机组推力系数，c₁为无量纲混合长， 

8.根据权利要求7所述的物理预测方法，其特征是所述计算各台风电机组轮毂高度的风速修正值具体利用公式 

其中， 

为各台风电机组轮毂高度的风速修正值， 

为各台风电机组轮毂高度的风速调节值，Δu为各台风电机组轮毂高度的风速衰减值。

9.根据权利要求8所述的物理预测方法，其特征是所述计算出各台风电机组的发电功率具体利用公式 

其中， 

为第n台风电机组的发电功率，C_p为风能利用系数，ρ₀为标准空气密度且ρ₀＝1.225kg/m³，A为风力机扫风面积， 

为第n台风电机组轮毂高度风速调节值。 

Images