CN104331621A - 一种风资源计算方法 - Google Patents

Abstract

本发明提供了一种风资源计算方法，通过对原始测风数据进行分析和基于全高度大气层边界建立CFD计算所需模型；使用CFD求解器计算得到各个风向条件下的风速和湍流分布；并根据各个风向条件下的风速和湍流分布、测风塔坐标、机位坐标、风机功率数据以及风况数据计算出风资源分布情况，确定风电场中适合安装风电机组的位置；也可以使用本发明的方法对已经设置好风电机组的风电场的风资源进行分析。

Description

一种风资源计算方法

技术领域

本发明涉及风电场技术领域，尤其涉及一种用于解决复杂地形风电场微观选址的风资源计算方法。

背景技术

微观选址是风电场建设的最重要阶段，主要任务是根据宏观选址的风场进行测风工作，使用测风点数据进行风场内风资源评估，计算出整个风场的风功率，根据地形条件和风分布进行风机选型和风机布点工作。风电场微观选址是在宏观选址的基础上布置风力发电机组，使整个风电场具有较好的经济效益。国内外的经验教训表明，风电场微观选址的失误造成的发电量损失和增加的维修费用将远远大于对场址进行详细调查的费用。因此，风电场的微观选址对风电场的建设至关重要。目前国内外的微观选址通常采用CFD软件，通过输入风能、气象、地形、地貌等各种数据，经过计算机的复杂计算来完成。

其中，风能是由地球表面大量空气流动所产生的动能，而人类可以利用的风能主要集中在大气边界层，大气边界层是大气层底部流动状态受到地表影响的部分，其中最下面的称为表面层，目前多数模型都是基于表面层设置的，但是，大气边界层表面层通常只有一百米左右的高度，而目前大型风电机组的叶轮所在空域通常超过这个高度，所以基于风电场风资源评估技术，设计一种专门用于复杂地形风电场微观选址的方法是有必要的。

发明内容

本发明提供一种用于优化风电场设计、风电机组选型和排布的风资源计算方法，通过模拟已知空间点的测风数据来推算其他空间点的风况，进而根据所有空间点的风况确定适合安装风电机组的位置，也可以对已经设置好风电机组的风电场的风资源进行分析。

本发明解决其技术问题所采用的技术方案是：一种风资源计算方法，包括以下步骤：

S1，对原始地形数据进行处理，生成三维立体的计算域网格；

S2，对原始测风数据进行统计分析得到全年的风况统计数据，所述风况统计数据包括：风速分布、风速频率、风向频率；

S3，根据所述风况数据生成多个风向边界条件，通过所述计算域网格、风向边界条件和纬度、地表粗糙度得到CFD计算所需模型；所述CFD计算模型为全高度大气边界层数值计算模型，包括：壁面函数、湍流模型和风廓线计算模型；

S4，根据CFD计算所需模型，使用CFD求解器计算得到各个风向条件下的风速和湍流分布；

S5，根据各个风向条件下的风速和湍流分布、测风塔坐标、机位坐标、风机功率数据以及风况数据计算出风资源分布情况；

1、壁面函数

壁面函数用来定义第一层网格在CFD计算所需模型；所述壁面函数包括湍流区域模型；所述湍流区域模型由传统的壁面函数：

$u^{+}=\frac{1}{\mathrm{\κ}}\ln \left(\frac{z_p}{z_0}\right)---\left(1\right),$

在考虑修正系a后得到等式： $u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{z_p}{z_0}\right)---\left(2\right);$

其中，u⁺是无量纲壁面切线速度；z₀为地表粗糙度；κ为冯卡门常数；z_p是从近壁单元中心到壁面的距离；a为修正系数由试验测量所得；

优选的，

a1)将所述等式(2)，考虑空气的运动粘度v、摩擦速度u_τ0后得到等式：

$u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{z_p}{z_0}\frac{u_{\mathrm{\τ}0}}{v}\frac{v}{u_{\mathrm{\τ}0}}\right)---\left(3\right);$

a2)将所述等式(3)经变换得到等式：

$u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{u_{\mathrm{\τ}0}{z}_p}{v}\frac{v}{z_0{u}_{\mathrm{\τ}0}}\right)---\left(4\right);$

a3)将无量纲壁面切线距离和地面无量纲粗糙长度带入所述等式(4)中，得到等式：

$u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{z^{+}}{z_0^{+}}\right)---\left(5\right);$

其中，z⁺为无量纲壁面切线距离；

优选的，所述第一层网格在水平面上的投影为矩形，网格的边长为30m-50m；所述第一层网格的垂直高度一般在10m以内。

2、湍流模型

湍流模型常数的标定为 $C_{\mathrm{\μ}}={\left(\frac{u^{*2}}{k}\right)}^2---\left(7\right);$

其中，u^*为地表摩擦速度；k为湍流脉动动能；所述u^*和k均为测量得到。

3、风廓线计算模型

所述风廓线计算模型为所述计算域网格边界部的计算模型，取所述计算域网格除第一层以外的区域，并且在每层计算域网格上都进行相同的计算；

所述风廓线计算模型包括：湍动能耗散率ε；湍动能耗散率由传统的大气边界层表面层的湍动能耗散率在考虑系数后得到

$\mathrm{\ϵ}=\frac{u^{*3}}{\mathrm{\κz}}{(1-\frac{z}{z_i})}^3---\left(8\right);$

其中，u^*为地表摩擦速度，由测量得到；κ为冯·卡门常数；κ＝0.4；z为所述计算域网格中心点的高度；z_i为在所述计算域大气边界层的厚度；

优选的，所述风廓线计算模型还包括：湍流脉动动能k的计算模型；所述湍流脉动动能在大气边界层表面层的模型是由常用等式

$k=\frac{u^{*2}}{\sqrt{C_{\mathrm{\μ}}}}---\left(9\right),$

在考虑系数： $(-20.31{(1-\frac{z}{z_i})}^{4}35.03{(1-\frac{z}{z_i})}^3-17.88{(1-\frac{z}{z_i})}^2-3.67(1-\frac{z}{z_i})-0.136)$ 后得到全高度大气边界层的模型：

$k=\frac{u^{*2}}{\sqrt{C_{\mathrm{\μ}}}}(-20.31{(1-\frac{z}{z_i})}^{4}35.03{(1-\frac{z}{z_i})}^3-17.88{(1-\frac{z}{z_i})}^2-3.67(1-\frac{z}{z_i})-0.136)---\left(10\right);$

其中，u^*为地表摩擦速度，C_μ为湍流模型常数；模型系数 $(-20.31{(1-\frac{z}{z_i})}^{4}35.03{(1-\frac{z}{z_i})}^3-17.88{(1-\frac{z}{z_i})}^2-3.67(1-\frac{z}{z_i})-0.136)$ 由实测数据分析和理论研究获得；

优选的，所述C_μ为常数；优选的C_μ＝0.036；

优选的，所述风廓线计算模型还包括：在全高度大气边界上的风廓线模型；所述在全高度大气边界上的风廓线模型由大气边界层表面层风廓线模型：

$u=\frac{u^{*}}{\mathrm{\κ}}\ln \left(\frac{z}{z_0}\right)---\left(11\right)$

在考虑全高度大气边界层修正后，得到：

$u=\frac{u^{*}}{\mathrm{\κ}}(\ln \left(\frac{z}{z_0}\right)+\frac{z}{L_M}-\frac{z}{z_i}\·\left(\frac{z}{2L_M}\right))---\left(12\right)$

其中，u为来流风速，u^*为地表摩擦速度由测量得到，κ为冯·卡门常数，z为所述计算域网格中心点的高度，z_i为在所述计算域大气边界层的厚度,L_M为大气边界层中部的长度尺度；z₀为地表粗超度长度；

对原始测风数据进行统计分析还包括对缺失数据检测，所述缺失数据检测包括：依据测风时间起始、终了点逐月建立测风记录完整时间表，通过完整时间表与实际测风记录时间对照查询缺失测风记录的数量和位置，依照所述原始测风数据中与缺失数据相似条件推算恢复出缺失数据；所述相似条件包括：时间、温度、风力、风向、大气压强；

根据所述风向边界条件、风廓线、湍流模型、壁面函数使用C语言编译软件编译成可被所述CFD求解器执行的动态运行库；其中，风向边界条件为16个，所述每个风向间相隔22.5度；

优选的，所述CFD求解器根据每个所述风向，依次运行；优选的，所述求解器运行16次并生成16个算例文件和对应的数据文件，得到与所述来风最为接近一个风向。

所述测风塔坐标的计算方法：根据所述测风塔的位置、高度经过坐标平移后得到所述测风塔在所述计算域网格中的三维坐标；

所述机位坐标为所述机位在所述计算域网格中的三维坐标，根据所述机位的位置、高度经过坐标平移后得到所述机位在所述计算域网格中的三维坐标；

使用所述CFD求解器计算得到所述风机功率数据方法为：

1)根据测风仪和风电机组轮毂在计算域网格上的三维坐标，运行求解器导出空间风电机组的参数，包括风电机组轮毂三维坐标值、风速、湍流动能、速度分量、耗散率、速度角；

2)根据测风塔测风仪位置计算数据和实测数据、指定风电机组轮毂高度位置计算数据以及风电机组功率曲线，生成机组功率数据；

优选的，所述步骤1调用CFD求解器16次，获取16个风向所在扇区的数值计算结果，同时确定所述风电机组的风向所在的风向扇区；

优选的，根据所述测风统计数据、风电机组功率曲线以及CFD数值计算结果，计算风电机组的相关数据，包括全年的平均风速、风电机组年出力和风电机组容量系数；

所述风电机组出力时间序列和风资源分布计算方法为：

1)根据测风数据和风电机组的功率曲线，得到风电机组轮毂位置的平均风速表、风电机组年出力表、风电机组容量系数表和所述风电机组出力时间序列，所述风电机组出力时间序列为与测风数据时间相对应的10分钟或1小时平均出力的时间序列；

2)根据测风数据和风电机组的功率曲线，计算等高面数据，包括曲面平均风速、曲面年出力、曲面容量系数；

优选的，所述机组为多个；

优选的，对缺失的测风数据可按用户要求进行线性插补。

本发明的有益效果是，

1)提供一种用于优化风电场设计、风电机组选型和排布的风资源计算方法，通过模拟已知空间点的测风数据来推算其他空间点的风况，进而根据所有空间点的风况确定适合安装风电机组的位置，也可以对已经设置好风电机组的风电场的风资源进行分析；

2)改进后的基于RANS方程和k－ε湍流模型的CFD计算方法更加适应大气边界层流动的特点；

3)经过实际测风数据修正的湍流模型和改进后的壁面函数、风廓线模型可以应用于风电场的流场计算，基于流场计算的风资源评估结果与实际状况相符；

4)只需操作几个简单的界面，调用相应的程序和求解器即可很方便的实现风电场风资源计算前期的数据录入、生成网格、建模，中期的计算和结果生成，后期对生成结果数据进行组织和诠释，并以直观可视的图形形式显示出来；

5)生成的矩形计算域网格质量高，设计研发了先进的大气边界层微尺度数值计算模型，使CFD计算结果更精确、效率更高，风向边界条件设置合理，运用CFD模型对风况进行计算流体力学(CFD)模拟,并通过CFD计算得出风电场出力的理论值，完成风电机组的最优布局和评估风场的运行情况、预测风场出力，最终达到充分利用风资源、获得最大发电量的目的；

6)本发明的计算模型采用全高度大气层边界，这样不仅可以满足表面层模型同时也延伸到更高空间区域，适合于50-300m高度区间的风速计算。

附图说明

下面结合附图对本发明所述的风资源计算方法进行具体说明。

图1是本发明计算方法的流程图；

图2是本发明风资源计算方法的计算域的坐标平移图；

图3是本发明风电场流场与风资源计算软件各个模块之间的关系图；

图4是本发明16个风向扇形区示意图。

具体实施方式

下面结合附图对本发明风资源计算方法作进一步详细说明，但本发明的实施方式不限于此。

如图1为本发明风资源计算方法的流程图，该方法包括以下步骤：

(一)对原始地形数据进行平滑、坐标平移等处理生成计算域网格，并根据计算域网格获得网格点间距和计算域侧面面积。

该计算域网格为三维立体网格，按水平高度分为多层，其中，第一层网格区域为矩形区域，该矩形区域形成的面与地形相配合，最高一层的网格形成的面为水平面。

在矩形区域界定的范围内，距离实际地形表面一定高度的位置，取一定数量的点，使用绘图工具将所述一定数量的点与所述的矩形框拟合成一个面；其中按照矩形边界生成网格，这些网格为规则的正方形，所有正方形的边长可以一致，也可以不一致，当正方形边长不一致时，在测风场中心区域内，正方形的边长较小(例如，边长可以取30m左右)，随着对测风场中心区域距离的增加，正方形的边长逐渐增大(例如，边长可以取60m左右)。第一层网格内所有正方形沿与地心垂直的高度向上延伸至计算域的最高层形成多个柱状体，沿每个柱状体的高度，从低到高等比设置一定数量的点，并将相应的点连接形成多层矩形网格。

对第一层网格的处理，包括对不规则的地形数据的四边进行平滑处理、坐标平移，并生成计算域网格。

1)准备好风电场地形数据文件，(通常文件名为combine*.xyz，即文件为xyz坐标点格式)，如果风电场原始测绘数据文件为其他格式，要先转换为xyz格式文件(xyz格式文件为多行三列文本文件格式，每行三个数据表示一个空间坐标点的x、y、z坐标，数据点按x-y坐标从南向北、从西向东排列)。CFD计算域为矩形区域，依据风电场周边地形状况向四周外延数公里，若测绘数据不包括风电场周边地形数据，可采用公开的卫星测绘数据，例如：采用ASTER GDEM30卫星测量数据，或SRTM90卫星测量数据。若采用卫星测量数据，需要将风电场测绘数据和卫星测量数据进行合并处理，生成矩形域的地形数据文件以备本系统使用。

2)输入边缘平滑参数；

CFD计算采用的计算域地形四边是水平的直线，而实际地形是高低起伏的，因此实际地形的边缘需要进行平滑处理，即第一层矩形网格的四个端点在同一水平高度，同时矩形的四边也在与端点相同的水平高度。被平滑处理的区域距离风电场边界要有数公里的距离，以便在风电场边界形成接近实际风况的流场。

计算域边缘平滑处理的参数包括边缘高度(Brim elevation)和地形分辨率(Resolution)，边缘高度也就是风电场地形的基准高度，该参数可以设定在边缘线的最低高度和平均高度之间，为了提高计算速度，可以假定边缘高度为0，其他点的高度参照边缘高度取相对高度地形分辨率在生成原始地形数据文件时可以设定，地形分辨率和网格分辨率一致，因为地形底边高低不同，需要把四边进行平滑处理成直线边，这里平滑步长都设为40，即矩形区域内的正方形边长为40m(当然，根据具体的地形，也可设为其他值)，分辨率设为40。

3)坐标平移处理

同样的，为了提高计算速度，可以将计算域网格的地形数据由GPS得到的绝对坐标调整为相对坐标，这里，我们可以假定第一层网格矩形西南顶点为坐标原点(也可以选取其他点为坐标原点)，对计算域网格内所有网格点进行坐标平移，并将经过平滑和坐标平移后的计算域数据保存。同时将三个坐标方向平移的距离也进行保存，以便在CFD计算后处理环节中还原地形坐标位置而使用。

4)输入网格划分参数；

当计算域网格进行完平滑处理、坐标平移后，需要重新生成计算域网格，生成计算域网格可以利用本系统进行，也可以利用网格生成软件来进行；所述利用本系统进行首先输入网格划分参数，网格划分参数包括：计算域从水平面起算的高度(Depth of theDomain)、x、y、z三个方向的网格尺度(也就是正方形的边长Grid Spacing)；高度方向的网格步长递增比例(正方形沿与地心垂直的高度向上延伸，至计算域的最高层等比设定的点z successive ratio)；还包括输出的地形数据文件包括的坐标点的列数(columns)、总共点数(total number of points)，以及地形数据的延展范围，即东西距离lx，南北距离ly和高度范围lz。

平滑和坐标平移后的地形数据文件包括的坐标点数不宜过多，如果坐标点数过多，则网格生成速度可能会超过30分钟甚至因内存不足而崩溃，一旦出现上述情况时，需要调整网格划分参数。

5)计算域网格生成

包括地表第一层网格的高度(First Cell Highness，第一层网格的垂直高度一般在10m以内)，网格总数(Total Cells)，计算域的水平投影面积(Total Horizontal Area)，南、北侧面面积(N-S Vertical Area)和东、西侧面面积(W-E Vertical Area)，如果第一层网格高度或网格总数不合适，可以高度参数或者网格间距等参数重新生成计算域网格；另外，还可以对风电场区域进行网格加密处理；加密处理可以使用现有软件，例如Gambit或ICEM。

第一层网格高度的控制算法如下：

$\frac{1}{Z_P}=\frac{1}{2z_p}\ln \left(\frac{{2z}_p}{z_0}\right)\⇒z_p=3.69z_0$

其中，z_p为第一层网格中心点至地表的距离，z₀为地表粗糙度，若第一层网格尺度需要满足上式要求，则标准壁面函数中的湍流动能生成率是合理的，如果z_p＞3.69z₀，则湍流动能生成率偏小。

第一层网格在水平面上的投影为矩形，网格的边长为30m-50m。

所生成的计算域网格保存到文件terrain.msh中，该文件可以被CFD求解器读入进行数值计算。

计算域网格侧面面积根据所述计算域网格中每个投影为正方形的矩形的四个端点在所述计算域网格上的相对坐标计算出每个投影为正方形的矩形的面积。

(二)对原始测风数据进行统计分析得到全年的风况统计数据(风况统计数据包括：风速分布、风速频率、风向频率)

1)导入测风塔水平坐标文件，或者手动输入坐标数据；

测风塔水平坐标文件中每一行数据代表一个测风塔的(x、y)坐标，如果只有一个测风塔，则测风塔文本文件只有一行数据；

然后输入测风塔的z坐标，测风塔z坐标为测风塔基座的地表高度，测风仪的坐标可以取测风塔的坐标加上测风仪距地高度。

2)导入原始测风数据文件，在导入原始测风数据前，先要对缺失数据进行检测，即执行缺失记录检查和补齐的工作，缺失数据检测可以，依照所述原始测风数据中与缺失数据相似条件(如：时间、温度、风力、风向、大气压强等)推算恢复出缺失数据，然后按月测风记录标识，数据统计分析

如果测风数据文件包括多个高度的测风数据，则可以进行对数风廓线拟合计算。

3)对测风数据的分析，包括全部数据和各个扇区的数据统计，即平均风速、威布尔分布参数、概率密度曲线、风向频率、bins数据统计、风能玫瑰图等，相应的参数值在文本显示区显示。

原始的测风数据文件，文件扩展名是.csv，如文件mast21.csv，其形式如下：

Mast No.：0021	V_No_40	D_No_40	H_40
				2004/10/6 15：20	5.7	79	40
2004/10/6 15：30	5.8	83	40
				2004/10/6 15：40	5.2	79	40
2004/10/6 15：50	5.1	83	40
				2004/10/6 16：00	5.2	84	40
2004/10/6 16：10	5.3	90	40
				2004/10/6 16：20	5.4	87	40
2004/10/6 16：30	6.2	90	40
				…	…	…	…

该文件的第一行为说明文字，第一列为文本格式的时间，第二列为该时段的平均风速，第三列为该时段的平均风向(以正北为0度按顺时针方向)，第四列为测风仪的安装高度。如果包括同一测风塔不同安装高度的测风仪数据，数据的时间间隔应相同，较低高度的数据放在前面。测风时间一般应持续一年，测风时间间隔可以为10分钟或60分钟。其中，测风数据可以以如下方式显示：

a.只有风速和风向数据；

b.时间(时间数据为日期)、风速和风向；

c.时间向量(将时间按照年、月、日、时、分秒分别独立显示)、风速和风向。

本实施例使用第一种输出数据文件进行CFD结果处理，注意输出文件不包含测风仪高度和测风时间间隔，这两个数据在后处理时需要手动输入。

4)将测风分析结果保存成风频统计数据文件、测风数据统计文件，

其中本发明分的风频统计数据文件为16个，按16个风向统计(16个风向如图4所示，图4中仅标出几个风向角度，其他省略，可以参考后面的数据)，每个风向所在扇形区域的角度为22.5度，依次为：11.25、33.75、56.25、78.75、101.25、123.75、146.25、168.75、191.25、213.75、236.25、258.75、281.25、303.75、326.25、348.75，这里的风频就是对某一个风向的次数占总的观测统计次数的百分比。风频也可以用玫瑰图表示，风频玫瑰图是一个地区多年的风向频率的统计图。

本发明的测风数据统计文件：包括各个扇区的平均风速、风向频率、威布尔概率密度函数的形状参数和尺度参数等，同时生成一个全年平均风速的文本文件meanspeed.txt，该文件可用于CFD初次计算的边界条件设置。

(三)根据测风数据生成多个风向边界条件，通过所述计算域网格、风向边界条件和纬度、地表粗糙度得到CFD计算所需的模型；CFD计算模型文件包括计算湍流模型、壁面函数和风廓线；

大气边界层是大气层底部流动状态受到地表影响的部分，其中最下面的称为表面层，目前多数模型都是基于表面层设置的，但是其局限性为，大气边界层表面层通常只有一百米左右的高度，目前大型风电机组的叶轮所在空域通常超过这个高度，所以本发明的计算模型采用全高度大气层边界，这样不仅可以满足表面层模型同时也延伸到更高空间区域，适合于50-300m高度区间的风速计算。

1、描述地表影响的壁面函数

这里提到的地表是指在计算域网格上的第一层，这里又可以分为：层流和湍流，主要根据距离壁面的远近进行区分的，在近壁处的流动与远离边壁区域的流动存在很大的差异：在靠近边壁的流层内，由于边壁的约束，流体质点基本不能垂直于边壁方向运动，而且流速梯度较大，粘滞切应力起主导作用，该薄层称为层流底层；而在远离边壁区的某一范围时，当空气流动的雷诺数大于临界雷诺数即为湍流，自然界的流动多为紊流。

因此层流底层的函数与湍流区域的函数表示也不相同，其中：

(a)在层流底层的函数为：

u⁺＝z⁺

其中，u⁺是无量纲壁面切线速度；z⁺为无量纲壁面切线距离；

(b)在充分湍流区域也称为对数区域：

$u^{+}=\frac{1}{\mathrm{\κ}}\ln \left(\frac{z_p}{z_0}\right)---\left(1\right)$

其中，式(1)中的u⁺是无量纲壁面切线速度；z₀为地表粗糙度；κ为冯卡门常数；z_p是从近壁单元中心到壁面的距离，在计算域网格上表示为第一层网格中心点至地表的距离，但该模型的实测数据与数值计算的风廓线有一定的偏差，因此采用系数a做为修正，得到：

$u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{z_p}{z_0}\right)---\left(2\right)$

(这里 $a=0.257z_p^{0.799}+0.985,$ 由实验测量得到)；

再将式(2)进行等式变换得到：

$u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{z_p}{z_0}\frac{u_{\mathrm{\τ}0}}{v}\frac{v}{u_{\mathrm{\τ}0}}\right)---\left(3\right)$

其中，v是空气的运动粘度，u_τ0是摩擦速度

然后将式(3)调整为：

$u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{u_{\mathrm{\τ}0}{z}_p}{v}\frac{v}{z_0{u}_{\mathrm{\τ}0}}\right)---\left(4\right)$

其中，将无量纲壁面切线距离和地面无量纲粗糙长度带入式(4)中，得到：

$u^{+}=\frac{a}{\mathrm{\κ}}\ln \left(\frac{z^{+}}{Z_0^{+}}\right)---\left(5\right)$

2、湍流模型

湍流模型通用CFD求解器中采用以粗糙度高度为参数的壁面函数，适合工业流动计算，不适合大气流动计算。本发明壁面函数直接使用地表粗超度长度尺度z0作为计算参数，和风工程中对地表的描述方法一致；

本发明使用的湍流模型为k-ε双方程湍流模型，湍流模型测风数据分析的目的是针对中性大气边界层流动确定湍流模型的常数C_μ，在测风数据对应的时期，取测风塔所在地的日出时间、日落时间和湍流强度最大时间的数据，有大约44％的数据可以满足“近中性”大气的条件，其中，计算式如下：

$k=\frac{1}{2}({\mathrm{\σ}}_u^2+{\mathrm{\σ}}_v^2+{\mathrm{\σ}}_w^2)=\frac{3}{2}{\mathrm{\σ}}_u^2---\left(6\right)$

式(6)中，k为湍流脉动动能，σ_u、σ_v、σ_w分别为纵向、横向和垂向瞬时速度的均方差，测风数据处理中假定测风高度均在大气边界层的表面层以内，不考虑湍流动能k沿高度的变化，k的取值为x、y、z3个方向分量(这里x、y、z3个方向是指：纵向、横向和垂向)的平均值，同时假定湍流流动是各向同性的，即纵向、横向和垂向的脉动速度(这里为标准差)都相等。

湍流模型常数的标定为：

$C_{\mathrm{\μ}}={\left(\frac{u^{*2}}{k}\right)}^2---\left(7\right)$

式(7)中，u^*为地表摩擦速度，其计算方法为：根据一年内测到测风仪所在点的平均数据，如果测风仪为多个，则取多个测风仪的平均值；

相比之下，国外现有测量数据计算的C_μ的数值在0.013～0.09之间，较常采用的数值为0.03，而现有数据与国内实际测量值有出入，本发明为了减少模型与实测值之间差异，因此是采用大量实测数据对模型进行修订，并且经过CFD验证计算，采用本发明湍流模型常数的数值计算模型与中性大气边界层的平均风速对数分布模型吻合较好。

3、风廓线计算模型

与壁面函数不同，风廓线计算模型为边界外部，可以取计算域网格除第一层以外的区域，并且在每层计算域网格上都进行相同的计算，这里介绍任意一层的网格计算。

(a)湍动能耗散率ε

对于大气边界层表面层，通常采用表示湍动能耗散率的变化，而对于全高度大气边界层，考虑到表面层以上湍动能的耗散率小于该式表示的值，设计了系数因此湍动能耗散率的变化在全高度大气边界层的模型为：

$\mathrm{\ϵ}=\frac{u^{*3}}{\mathrm{\κz}}{(1-\frac{z}{z_i})}^3---\left(8\right)$

式(8)中，u^*为地表摩擦速度，由测量得到，κ为冯卡门常数这里取0.4，z为计算域网格的高度，z_i为在计算域大气边界层的厚度；

(b)湍流脉动动能

通常大气边界层表面层的湍动能采用常数假设为：

$k=\frac{u^{*2}}{\sqrt{C_{\mathrm{\μ}}}}---\left(9\right),$

但对于全高度大气边界层，在近地区域湍动能趋近于0，在大气边界层上界湍动能也应趋近于0，因此建立了系数：

$(-20.31{(1-\frac{z}{z_i})}^{4}35.03{(1-\frac{z}{z_i})}^3-17.88{(1-\frac{z}{z_i})}^2-3.67(1-\frac{z}{z_i})-0.136),$ 将系数代入式(9)中，得到全高度大气边界层的湍动能采用假设，

$k=\frac{u^{*2}}{\sqrt{C_{\mathrm{\μ}}}}(-20.31{(1-\frac{z}{z_i})}^{4}35.03{(1-\frac{z}{z_i})}^3-17.88{(1-\frac{z}{z_i})}^2-3.67(1-\frac{z}{z_i})-0.136)---\left(10\right);$

式(10)中，C_μ为湍流模型常数，模型系数由实测数据分析和理论研究获得，u^*为地表摩擦速度；

(c)风廓线模型：

通常采用的大气边界层表面层风廓线模型为：

$u=\frac{u^{*}}{\mathrm{\κ}}\ln \left(\frac{z}{z_0}\right)---\left(11\right)$

式(11)中，u为来流风速，u^*为地表摩擦速度，κ为冯卡门常数，z为所在层低一层网格的高度，z₀为地表粗超度长度；

大气边界层表面层通常只有一百米左右的高度，目前大型风电机组的叶轮所在空域通常超过这个高度，对于全高度大气边界层，建立了修正后的风速廓线模型：

$u=\frac{u^{*}}{\mathrm{\κ}}(\ln \left(\frac{z}{z_0}\right)+\frac{z}{L_M}-\frac{z}{z_i}\·\left(\frac{z}{2L_M}\right))---\left(12\right)$

式(12)中，u为来流风速，u^*为地表摩擦速度，κ为冯卡门常数，z为所述计算域网格中心点的高度，z_i为在所述计算域大气边界层的厚度，L_M为大气边界层中部的长度尺度，z₀为地表粗超度长度。

4、平均风速数据包括：x坐标、y坐标、全年平均风速、测风仪距离地面高度，在本实施例中，测量到的数据如下表：

	x坐标	y坐标	全年平均风速	测风仪距离地面高度
					1	3.95068e+007	2.59506e+006	7.6045	40

5、在本实施例中，测量到的风电场的纬度和地表粗糙度如下表：

Latitude(纬度)：	40.0	Deg
			Roughness Length(z0)(地表粗糙	0.03	m

度)：

(四)根据CFD计算模型条件，调用第三方求解器计算得到CFD计算结果，将上述风向边界条件、风廓线、湍流模型、壁面函数使用C语言编译软件编译成动态运行库，该动态运行库与CFD求解器相链接，可以在CFD求解中调用。

1)导入前面获得的计算域网格文件、16个风向的边界条件、风廓线、湍流模型、壁面函数等CFD模型；

2)根据CFD模型，依次运行16次CFD计算程序，该程序模块外部调用CFD求解器，计算完成后，生成16个指令文件和对应的数据文件，所述指令文件和数据文件通过CFD软件或CFD后处理软件打开和提取数据。

(五)根据各个风向条件下的风速和湍流分布、测风塔坐标、机位坐标、机组功率曲线以及测风数据计算出风资源分布情况和机组出力时间序列

1、根据以上步骤获得的CFD模型，调用CFD求解器进行CFD数值计算，所述CFD数值计算的结果包括计算域内空间点在每个时间间隔内的风速、风向、湍流动能等数据；

2、测风塔坐标

根据所述测风塔的位置、高度经过坐标平移后得到所述测风塔在所述计算域网格中的三维坐标

3、机位坐标

根据所述机位的位置、高度经过坐标平移后得到所述机位在所述计算域网格中的三维坐标

4、机组功率曲线

1)根据测风仪和风电机组轮毂在计算域网格上的三维坐标，运行求解器导出空间风电机组的参数，包括风电机组轮毂三维坐标值、风速、湍流动能、速度分量、耗散率、速度角；

2)根据指定风电机组轮毂高度，生成机组功率曲线数据；

3)根据曲线数据再次运行求解器获取等高曲面上的风速等数据，由于有16个风方向，因此程序将调用CFD求解器16次，获取16个风向所在扇区的数值计算结果，在16个计算结果中确定当风电机组的风向属于哪个风向扇区；

4)指定风功率曲线文本文件之后，再输入风电机组额定功率值，输入功率值(比如1000)后，会显示该风电机组的出力曲线，如图

5、风资源分布情况

根据测风统计数据、风电机组功率曲线以及CFD数值计算结果，计算风电机组的相关数据，包括全年的平均风速、风电机组年出力、风频玫瑰图和风电机组容量系数；

再统计风电场所有风电机组的年出力情况，计算出当前风电场的风资源情况；

6、机组出力时间序列

1)由风电机组的风资源数据，可以得到风电机组平均风速表、风电机组年出力表、风电机组容量系数表；

2)同时根据风电机组的风资源数据可以计算等高面的相关数据，包括曲面平均风速、曲面年出力、曲面容量系数、风电机组切出时间分布；

3)根据等高面风资源情况可以得到：等高面平均风速分布图、等高面年出力分布图、等高面容量系数分布图(所述登高面容量系数为同一高度面上风电机组的发电量分布)；

4)最后，输出文件包括三个：曲面坐标数据文件(x，y，z)、曲面平均风速数据文件(x，y，Vm)，曲面容量系数分布文件(x，y，Cf)。

以上所述仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专利的技术人员在不脱离本发明技术方案范围内，当可利用上述提示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明方案的范围内。

Claims (9)

1.一种风资源计算方法，包括以下步骤： 

S1，对原始地形数据进行处理，生成三维立体的计算域网格； 

S2，对原始测风数据进行统计分析得到全年的风况统计数据，所述风况统计数据包括：风速分布、风速频率、风向频率； 

S3，根据所述风况数据生成多个风向边界条件，通过所述计算域网格、风向边界条件和纬度、地表粗糙度得到CFD计算所需模型； 

S4，根据CFD计算所需模型，使用CFD求解器计算得到各个风向条件下的风速和湍流分布； 

S5，根据各个风向条件下的风速和湍流分布、测风塔坐标、机位坐标、风机功率数据以及风况数据计算出风资源分布情况； 

其特征在于，所述CFD计算模型为全高度大气边界层数值计算模型。 

2.如权利要求1所述的风资源计算方法，其特征在于，所述全高度大气边界层数值计算模型包括：壁面函数；所述壁面函数用来定义第一层网格在CFD计算所需模型；所述壁面函数包括湍流区域模型；所述湍流区域模型由传统的壁面函数(1) 在考虑修正系a后得到等式(2)：

其中，u⁺是无量纲壁面切线速度；z₀为地表粗糙度；κ为冯卡门常数；z_p是从近壁单元中心到壁面的距离；a为修正系数由试验测量所得； 

优选的， 

a1)将所述等式(2)，考虑空气的运动粘度v、摩擦速度u_τ0后得到等式(3)： 

a2)将所述等式(3)经变换得到等式(4)：

a3)将无量纲壁面切线距离和地面无量纲粗糙长度带入所述等式(4)中，得到等式(5)：

其中，z⁺为无量纲壁面切线距离； 

优选的，所述第一层网格在水平面上的投影为矩形，网格的边长为30m-50m；所述第一层网格的垂直高度一般在10m以内。 

3.如权利要求1所述的风资源计算方法，其特征在于，所述全高度大气边界层模型还包括：湍流模型常数的标定，其中，u^*为地表摩擦速度；k为湍流脉动动能；所述u^*和k均为测量得到。 

4.如权利要求1所述的风资源计算方法，其特征在于，所述全高度大气边界层模型还包括：风廓线计算模型；所述风廓线计算模型为所述计算域网格边界部的计算模型，取所述计算域网格除第一层以外的区域，并且在每层计算域网格上都进行相同的计算； 

所述风廓线计算模型包括：湍动能耗散率ε；湍动能耗散率由传统的大气边界层表面层的湍动能耗散率在考虑系数后得到

其中，u^*为地表摩擦速度，由测量得到；κ为冯·卡门常数；κ＝0.4；z为所述计算域网格中心点的高度；z_i为在所述计算域大气边界层的厚度； 

优选的，所述风廓线计算模型还包括：湍流脉动动能k的计算模型；所述湍流脉动动能在大气边界层表面层的模型是由常用等式(9)在考虑系数： 后得到全高度大气边界层的模型： 

其中，u^*为地表摩擦速度，C_μ为湍流模型常数；模型系数 由实测数据分析和理论研究获得； 

优选的，所述C_μ为常数；优选的C_μ＝0.036； 

优选的，所述风廓线计算模型还包括：在全高度大气边界上的风廓线模型；所述在 全高度大气边界上的风廓线模型由大气边界层表面层风廓线模型(11)：在考虑全高度大气边界层修正后，得到： 

其中，u为来流风速，u^*为地表摩擦速度由测量得到，κ为冯·卡门常数，z为所述计算域网格中心点的高度，z_i为在所述计算域大气边界层的厚度,L_M为大气边界层中部的长度尺度；z₀为地表粗超度长度。 

5.如权利要求1-4任一项所述的风资源计算方法，其特征在于，对原始测风数据进行统计分析还包括对缺失数据检测，所述缺失数据检测包括：依据测风时间起始、终了点逐月建立测风记录完整时间表，通过完整时间表与实际测风记录时间对照查询缺失测风记录的数量和位置，依照所述原始测风数据中与缺失数据相似条件推算恢复出缺失数据；所述相似条件包括：时间、温度、风力、风向、大气压强。 

6.如权利要求1所述的风资源计算方法，其特征在于，根据所述风向边界条件、风廓线、湍流模型、壁面函数使用C语言编译软件编译成可被所述CFD求解器执行的动态运行库；其中，风向边界条件为16个，所述每个风向间相隔22.5度； 

优选的，所述CFD求解器根据每个所述风向，依次运行；优选的，所述求解器运行16次并生成16个算例文件和对应的数据文件，得到与所述来风最为接近一个风向。 

7.如权利要求1-6任一项所述的风资源计算方法，其特征在于，所述测风塔坐标的计算方法：根据所述测风塔的位置、高度经过坐标平移后得到所述测风塔在所述计算域网格中的三维坐标； 

所述机位坐标为所述机位在所述计算域网格中的三维坐标，根据所述机位的位置、高度经过坐标平移后得到所述机位在所述计算域网格中的三维坐标。 

8.如权利要求1-7任一项所述的风资源计算方法，其特征在于，使用所述CFD求解器计算得到所述风机功率数据方法为： 

1)根据测风仪和风电机组轮毂在计算域网格上的三维坐标，运行求解器导出空间风电机组的参数，包括风电机组轮毂三维坐标值、风速、湍流动能、速度分量、耗散率、速度角； 

2)根据测风塔测风仪位置计算数据和实测数据、指定风电机组轮毂高度位置计算数据以及风电机组功率曲线，生成机组功率数据； 

优选的，所述步骤1调用CFD求解器16次，获取16个风向所在扇区的数值计算结果，同时确定所述风电机组的风向所在的风向扇区； 

优选的，根据所述测风统计数据、风电机组功率曲线以及CFD数值计算结果，计算风电机组的相关数据，包括全年的平均风速、风电机组年出力和风电机组容量系数。 

9.如权利要求1-8任一项所述的风资源计算方法，其特征在于，所述风电机组出力时间序列和风资源分布计算方法为： 

1)根据测风数据和风电机组的功率曲线，得到风电机组轮毂位置的平均风速表、风电机组年出力表、风电机组容量系数表和所述风电机组出力时间序列，所述风电机组出力时间序列为与测风数据时间相对应的10分钟或1小时平均出力的时间序列； 

2)根据测风数据和风电机组的功率曲线，计算等高面数据，包括曲面平均风速、曲面年出力、曲面容量系数； 

优选的，所述机组为多个； 

优选的，对缺失的测风数据可按用户要求进行线性插补。