CN106121913A - 一种海上风电场的风机排布方法 - Google Patents

Abstract

本发明涉及一种海上风电场的风机排布方法，包括步骤：确定风机的排布阵列；获取风电场的约束条件数据；根据风机的排布阵列和风电场的约束条件筛选候选排布方案；计算各方案的尾流损失；检验最优方案。相对于现有技术，本发明的海上风电场的风机排布方法，充分考虑风机排布在形态、方位、间距等三个方面的变化，能够灵活适应多主导风向的情况与风电场的复杂几何边界，阵列筛选范围更广，效率更高。

Description

一种海上风电场的风机排布方法

技术领域

本发明属于海上风力发电场设计技术领域，具体涉及一种风电机组排布方案的快速生成与优化方法。

背景技术

风电机组的排布是海上风电场的关键参数。在设计与建造阶段，风机选型、交通运输、电气设计、土建工程、施工组织、环境保护等多个环节离不开风机排布作为输入条件；在生产运行阶段，风机偏航与其他控制策略也需要与风机排布方案相适应。

海上风电场由于地形平坦、下垫面类型单一，风能资源分布的空间差异性小，因此通常采取保持风机间距相对固定的规则化排布方案。形成的风机阵列在充分利用风能资源的基础上，控制用海面积，并减少随机排布带来的视觉冲击。国内外已建成的海上风电场所采用的主流风机阵列有矩形、平行四边形和扇形等三大类。

风机排布的优劣直接反映在风机间相互干扰产生的尾流强弱，并最终影响整个风电场的发电效率与风机安全。目前在大型海上风电场设计过程中，风机排布主要遵从如下原则：

(1)阵列方向与风电场主导风向垂直，多列排布时前后行错开呈“梅花”型；

(2)行间距不低于7倍风机叶轮直径，列间距不低于5倍风机叶轮直径；

(3)单台风机因尾流导致的发电量损失不超过15％。

上述做法存在若干缺陷，如无法考虑两个甚至多个主导风向，无法指导扇形阵列，风机间距的确定带有经验性，无法保证尾流损失最低等。因此，海上风电场风机排布的生成与优化面临代表性有限和效率低等问题。

发明内容

本发明是通过以下技术方案实现的：

一种海上风电场的风机排布方法，包括步骤：

S1：确定风机的排布阵列；

S2：获取风电场的约束条件数据；

S3：根据风机的排布阵列和风电场的约束条件筛选候选排布方案；

S4：计算各方案的尾流损失；

S5：检验最优方案。

相对于现有技术，本发明的海上风电场的风机排布方法，充分考虑风机排布在形态、方位、间距等三个方面的变化，通过定义风机排布基本单元，生成平行四边形和扇形阵列风机坐标，利用枚举法筛选出满足风电场约束条件的所有候选排布方案，引入尾流数学模型，对候选方案逐一进行计算，选出平均尾流损失最低的推荐排布方案，能够灵活适应多主导风向的情况与风电场的复杂几何边界，阵列筛选范围更广，效率更高。

进一步，所述风机的排布阵列包括平行四边形和扇形两类。

进一步，风电场的约束条件数据包括风机排布区域边界数据和风机的数目。

进一步，所述步骤S3进一步包括步骤：

S31：确定西南、西北、东北、东南共4个固定点S1～S4；

S32：确定阵列排布的起始机位P0，坐标记为(X0，Y0)；

S32：确定阵列的行数n和列数m；

S33：计算阵列机位P(i,j)的坐标(Xi，Yj)，其中，i＝1,2,…n，j＝1,2,…m；

S34：剔除位于内、外边界之外的无效机位P″，有效机位P′(i,j)总数记为N′；

S35：循环控制，判断N′是否等于N，如果是，将P′列为候选排布方案，如果否，转到步骤S31。

进一步，所述步骤S4具体包括步骤：

S41：计算单个扇区风速损失；

S42：计算单个扇区功率损失；

S43：尾流计算循环控制，判断是否全部扇区的尾流计算都已经完成，如果是，转到步骤S44，如果否，转到步骤S41；

S44：计算平均尾流损失；

S45：排布方案循环控制，判断是否全部排布方案都已经计算完平均尾流损失，如果是，转到步骤S46，如果否，转到步骤S41；

S46：获得平均尾流损失最低的N个方案，N>1。

进一步，步骤S5进一步包括步骤：

S51：建立风向频率变动模型；

S52：计算平均尾流损失；

S53：评价变动度指标；

S54:推荐风机排布方案。

进一步，所述步骤S4具体为：利用成熟的N.O.Jensen模型，在0°至360°平均划分144个扇区对候选排布方案的逐台风机进行风速损失和功率损失计算，并用平均功率损失表示整个阵列的平均尾流损失。

进一步，所述步骤S5具体为：选出平均尾流损失最低的5个候选排布方案，利用风向频率变动模型重新计算各自的尾流损失，确定平均尾流损失的变化范围，选择变动度最小的候选排布方案作为推荐排布方案。

进一步，所述风机的排布阵列具体为：风机阵列的基本单元定义为风机排布方案中每4台风机的相对位置所形成的平面形态，

若基本单元为平行四边形，长对角线AC称为主对角线，短对角线BD称为次对角线，控制参数包括主对角线的方位角α，主、次对角线间的夹角β，以及主、次对角线的长度d1、d2；

若基本单元为扇形，短弧AD与长弧BC共圆，圆心为O；半径OA称为短径，OB称为长径，直线AB称为始边，直线CD称为终边，控制参数包括扇形始边的方位角α，始、终边间的夹角β，以及长、短径的长度d1、d2。

进一步，所述步骤S3具体的计算方法为：

若风机阵列为平行四边形，P0坐标由两条虚拟线L1、L2相交确定，L1、L2的斜率记为k1、k2，

当0°≤α<90°，L1、L2分别经过S1、S4；当90°≤α<180°，L1、L2分别经过S2、S1；当180°≤α<270°，L1、L2分别经过S3、S2；当270°≤α<360°，L1、L2分别经过S4、S3，P0(X0，Y0)则根据解析几何中的两直线交点公式确定；

n、m由两组虚拟平行线L1、L′1和L2、L′2共同确定，当0°≤α<90°，L′1、L′2分别经过S3、S2；当90°≤α<180°，L′1、L′2分别经过S4、S3；当180°≤α<270°，L′1、L′2分别经过S1、S4；当270°≤α<360°，L′1、L′2分别经过S2、S1；

记L1、L′1在东西方向上投影之间的距离为dX，记L2、L′2在南北方向上投影之间的距离为dY，平行四边形的边AB、边AD长度分别为x、y，可由α、β、d1、d2根据三角形余弦定理计算得到；

n、m表达式如下：

n＝dX/x

m＝dY/y

机位P_(i,j)的坐标(X_i，Y_j)表达式如下：

X_i＝X₀+(i–1)×x i＝1,2,…n，

Y_j＝Y₀+(j–1)×y j＝1,2,…m

若风机阵列为扇形，P₀坐标由S₁～S₄经过坐标旋转后确定，坐标旋转角为α+β/2。记旋

转后的S₁～S₄为S′₁～S′₄，重新计算X′_min、X′_max、Y′_min、Y′_max。旋转后的(X₀，Y₀)记为(X′₀，

Y′₀)，(X′₀，Y′₀)和(X₀，Y₀)表达式如下：

X′₀＝(X′_min+X′_max)/2

Y′₀＝Y′_min–d₂

X₀＝X′₀cos(α+β/2)+Y′₀sin(α+β/2)

Y₀＝X′₀cos(α+β/2)-Y′₀sin(α+β/2)

n、m表达式如下：

$n=\frac{\sqrt{{(X_{\max }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\max }^{\&prime;}-Y_0^{\&prime;})}^2}-\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\min }^{\&prime;}-Y_0^{\&prime;})}^2}}{d_1-d_2}$

$m=\&lsqb;a\tan \left(\frac{X_{\max }^{\&prime;}-X^{\&prime;}}{{(X_{\max }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\max }^{\&prime;}-Y_0^{\&prime;})}^2}\right)-a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)\&rsqb;\frac{180}{\mathrm{\&pi;}\mathrm{\&beta;}}$

机位P_(i,j)的坐标(X_i，Y_j)表达式如下：

$\begin{array}{c}{X}_i=X_0+\&lsqb;d_2+(i,1)(d_1-d_2)\&rsqb;\sin \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\&rsqb;\cos (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)+\\ \&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\sin (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)\end{array}$

$\begin{array}{c}{Y}_j=Y_0+\&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\&rsqb;\cos (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)-\\ \&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\sin (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)\end{array}$

进一步，所述步骤S4具体计算方法为：单台风机在单个扇区所受实际风速损失dv′为：

${\mathrm{dv}}^{\&prime;}=\sqrt{\stackrel{n}{\mathrm{\&Sigma;}}{\left(dv\right)}^2}$

$dv=\mathrm{\&phi;}\&lsqb;1-\sqrt{1-C_T}\&rsqb;{\left(\frac{R}{R+kL}\right)}^2$

其中，n为上风向尾流个数；CT为风机推力系数；R为风机叶轮半径；L为风机间距在

风向上的投影距离；k为尾流衰减系数，Φ为尾流遮蔽系数，表明风机被上风向尾流遮

蔽的程度，

单台风机在单个扇区所受实际功率损失dP′为：

${\mathrm{dP}}^{\&prime;}={\&Integral;}_{v_s}^{v_e}g\left(v\right)f\left(v\right)dv-{\&Integral;}_{v_s}^{v_e}g\left(v\right)f\&lsqb;v(1-{\mathrm{dv}}^{\&prime;})\&rsqb;dv$

其中，v_s、v_e为风机切入和切出风速；f(v)为风机功率曲线，g(v)为风速在该扇区的概率

密度分布，平均尾流损失dW′为：

${\mathrm{dW}}^{\&prime;}=\frac{\underset{j=1}{\overset{144}{\&Sigma;}}{\mathrm{\&omega;}}_j\&CenterDot;\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{dP}}_i^{\&prime;}}{N}$

其中，N为风机数目；ω为各个风向扇区的权重，根据海上风电场测风塔在12、16或48个方位的风向频率分布结果平均分配。

进一步，步骤S5具体为：

记测风塔的风向频率分布为J(θi)，i＝1,2,…h，h为风向扇区数目；θ为扇区风向，单位为度；主导风向记为θp，相应频率为J(θp)。风向频率变动模型设置8种变动情景，模拟风向在年际、年代际尺度上的变化：

(1)J(θp)增加30％，其余扇区的频率相应平均减少；

(2)J(θp)减少30％，其余扇区的频率相应平均增多；

(3)θp左移一个扇区，J(θp)保持不变；

(4)θp右移一个扇区，J(θp)保持不变；

(5)θp左移一个扇区，J(θp)增加30％，其余扇区的频率相应平均减少；

(6)θp左移一个扇区，J(θp)减少30％，其余扇区的频率相应平均增多；

(7)θp右移一个扇区，J(θp)增加30％，其余扇区的频率相应平均减少；

(8)θp右移一个扇区，J(θp)减少30％，其余扇区的频率相应平均增多。

对以上8种情景重新计算平均尾流损失dWi′,i＝1,2,…8，形成尾流损失区间dWmin′～dWmax′；变动度指标定义为dWmax′/dWmin′。

为了能更清晰的理解本发明，以下将结合附图说明阐述本发明的具体实施方式。

附图说明

图1是本发明的海上风电场的风机排布方法的步骤流程图。

图2是本发明的风机阵列的基本单元的示意图。

图3是本发明的风机阵列潜在基本单元的示意图。

图4是本发明的风电场约束条件的示意图。

图5是本筛选候选排布方案的步骤流程图。

图6是本发明的风机阵列生成的示意图。

图7是计算各方案的尾流损失的步骤流程图。

图8是检验最优方案的步骤流程图。

具体实施方式

本发明的海上风电场的风机排布方法整体技术方案主要包含以下几个核心过程：

(1)风机阵列基本单元的定义；

(2)风电场约束条件的定义；

(3)候选排布方案的筛选；

(4)风机尾流损失的计算；

(5)最优排布方案的检验。

如图1所示，本发明的海上风电场的风机排布方法，包括步骤：

S1：确定风机的排布阵列；

S2：获取风电场的约束条件数据；

S3：根据风机的排布阵列和风电场的约束条件筛选候选排布方案；

S4：计算各方案的尾流损失；

S5：检验最优方案。

以下对各步骤作详细介绍：

S1：确定风机的排布阵列。

风机阵列的基本单元定义为指风机排布方案中每4台风机的相对位置所形成的平面形态，包括平行四边形和扇形两类，均由两个角度和两个长度控制，如图2所示。

若基本单元为平行四边形，如图2所示，长对角线AC称为主对角线，短对角线BD称为次对角线。控制参数包括主对角线的方位角α，主、次对角线间的夹角β，以及主、次对角线的长度d1、d2。

若基本单元为扇形，如图2所示，短弧AD与长弧BC共圆，圆心为O；半径OA称为短径，OB称为长径，直线AB称为始边，直线CD称为终边。控制参数包括扇形始边的方位角α，始、终边间的夹角β，以及长、短径的长度d1、d2。

α取值范围为0°～359°，步长为1°。β取值范围为1°～179°(平行四边形，90°为矩形)或1°～359°(扇形)，步长为1°。d1、d2的取值范围为风机叶轮直径的任意倍数，步长为风机叶轮直径的十分之一；d2≤d1(平行四边形)或<d1(扇形)。

利用计算机编程，对α、β、d1、d2的不同取值进行组合，形成所有潜在的基本单元，如图3所示。

S2：获取风电场的约束条件数据。

风电场约束条件的定义包括风机排布区域边界定义和风机数目定义，如图4所示。

风电场资料包括：风电场外边界、限制性区域、预留区域、规划容量，由这些数据可以计算出风机排布区域；风机资料包括单机额定功率，结合规划容量，可以计算出风机数目。

风机排布区域由外边界和内边界组成，外边界为整个风电场的边界，内边界为风电场内部的限制性区域(不良地质区、急流区等)或预留区域(如海上升压站、测风塔等)的边界，例如采用北京54、西安80或UTM坐标系进行描述。风机排布限定在内、外边界之间。

风机排布的风机数目N决定于风电场规划容量E和所选机型的单机额定容量e：

N＝E/e

S3：根据风机的排布阵列和风电场的约束条件筛选候选排布方案。

如图5所示，候选排布方案的筛选包括以下步骤：

S31：确定西南、西北、东北、东南共4个固定点S1～S4；

S1～S4，坐标依次为(Xmin，Ymin),(Xmin，Ymax)，(Xmax，Ymax)，(Xmax，Ymin)；其中Xmin、Xmax为风电场边界在东西方向坐标的最小、最大值；Ymin、Ymax为南北方向坐标的最小、最大值。

S32：确定阵列排布的起始机位P0，坐标记为(X0，Y0)；

S32：确定阵列的行数n和列数m；

S33：计算阵列机位P(i,j)的坐标(Xi，Yj)，其中，i＝1,2,…n，j＝1,2,…m；

S34：剔除位于内、外边界之外的无效机位P″，有效机位P′(i,j)总数记为N′；

S35：循环控制，判断N′是否等于N，如果是，将P′列为候选排布方案，如果否，转到步骤S31。

对n,m还有P(i,j)的具体如下：

若阵列为平行四边形，如图6所示，P0坐标由两条虚拟线L1、L2相交确定，L1、L2的斜率记为k1、k2(与图2中平行四边形斜边AB、AD平行)，可由α、β、d1、d2根据三角形余弦定理计算得到。

当0°≤α<90°，L1、L2分别经过S1、S4；当90°≤α<180°，L1、L2分别经过S2、S1；当180°≤α<270°，L1、L2分别经过S3、S2；当270°≤α<360°，L1、L2分别经过S4、S3。(X0，Y0)则根据解析几何中的两直线交点公式确定。

n、m由两组虚拟平行线L1、L′1和L2、L′2共同确定，当0°≤α<90°，L′1、L′2分别经过S3、S2；当90°≤α<180°，L′1、L′2分别经过S4、S3；当180°≤α<270°，L′1、L′2分别经过S1、S4；当270°≤α<360°，L′1、L′2分别经过S2、S1。

记L1、L′1在东西方向上投影之间的距离为dX，记L2、L′2在南北方向上投影之间的距离为dY，图2中AB、AD长度分别为x、y，可由α、β、d1、d2根据三角形余弦定理计算得到，则n、m表达式如下：

n＝dX/x

m＝dY/y

机位P(i,j)的坐标(Xi，Yj)表达式如下：

Xi＝X0+(i–1)×x i＝1,2,…n，

Yj＝Y0+(j–1)×y j＝1,2,…m

若阵列为扇形，如图6所示，P0坐标由S1～S4经过坐标旋转后确定，坐标旋转角为α+β/2。记旋转后的S1～S4为S′1～S′4，重新计算X′min、X′max、Y′min、Y′max。旋转后的(X0，Y0)记为(X′0，Y′0)，(X′0，Y′0)和(X0，Y0)表达式如下：

X′0＝(X′min+X′max)/2

Y′0＝Y′min–d2

X0＝X′0 cos(α+β/2)+Y′0 sin(α+β/2)

Y0＝X′0 cos(α+β/2)-Y′0 sin(α+β/2)

n、m表达式如下：

$n=\frac{\sqrt{{(X_{\max }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\max }^{\&prime;}-Y_0^{\&prime;})}^2}-\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\min }^{\&prime;}-Y_0^{\&prime;})}^2}}{d_1-d_2}$

$m=\&lsqb;atan\left(\frac{X_{\max }^{\&prime;}-X^{\&prime;}}{\sqrt{{(X_{\max }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\max }^{\&prime;}-Y_0^{\&prime;})}^2}}\right)-atan\left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)\&rsqb;\frac{180}{\mathrm{\&pi;}\mathrm{\&beta;}}$

机位P(i,j)的坐标(Xi，Yj)表达式如下：

$\begin{array}{c}{X}_i=X_0+\&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\sin \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\&rsqb;\cos (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)+\\ \&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\sin (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)\end{array}$

$\begin{array}{c}{Y}_j=Y_0+\&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\&rsqb;\cos (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)-\\ \&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\sin (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)\end{array}$

S4：计算各方案的尾流损失。

如图7所示，步骤S4具体包括以下步骤：

S41：计算单个扇区风速损失；

S42：计算单个扇区功率损失；

S43：尾流计算循环控制，判断是否全部扇区的尾流计算都已经完成，如果是，转到步骤

S44，如果否，转到步骤S41；

S44：计算平均尾流损失；

S45：排布方案循环控制，判断是否全部排布方案都已经计算完平均尾流损失，如果是，转到步骤S46，如果否，转到步骤S41；

S46：获得平均尾流损失最低的N个方案，N>1。

具体地，利用经典的N.O.Jensen模型，在0°～360°平均划分144个扇区对候选排布方案的逐台风机进行风速损失和功率损失计算，并用平均功率损失表示整个阵列的平均尾流损失。

单台风机在单个扇区所受实际风速损失dv′为：

${\mathrm{dv}}^{\&prime;}=\sqrt{\stackrel{n}{\mathrm{\&Sigma;}}{\left(dv\right)}^2}$

$dv=\mathrm{\&phi;}\&lsqb;1-\sqrt{1-C_T}\&rsqb;{\left(\frac{R}{R+kL}\right)}^2$

式中，n为上风向尾流个数；CT为风机推力系数；R为风机叶轮半径；L为风机间距在风向上的投影距离；k为尾流衰减系数，根据海面条件取为0.04；Φ为尾流遮蔽系数，表明风机被上风向尾流遮蔽的程度，取值范围为0～1。

单台风机在单个扇区所受实际功率损失dP′为：

${\mathrm{dP}}^{\&prime;}={\&Integral;}_{v_s}^{v_e}g\left(v\right)f\left(v\right)dv-{\&Integral;}_{v_s}^{v_e}g\left(v\right)f\&lsqb;v(1-{\mathrm{dv}}^{\&prime;})\&rsqb;dv$

其中，v_s、v_e为风机切入和切出风速；f(v)为风机功率曲线；g(v)为风速在该扇区的概率密度分布，采用Weibull分布形式。

平均尾流损失dW′为：

${\mathrm{dW}}^{\&prime;}=\frac{\underset{j=1}{\overset{144}{\&Sigma;}}{\mathrm{\&omega;}}_j\&CenterDot;\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{dP}}_i^{\&prime;}}{N}$

其中，N为风机数目；ω为各个风向扇区的权重，根据海上风电场测风塔在12、16或48个方位的风向频率分布结果平均分配。这里N可以根据实际情况而定，例如可以取5，获得平均尾流损失最低的5个方案。

S5：检验最优方案。

选出平均尾流损失最低的5个候选排布方案，利用风向频率变动模型重新计算各自的尾流损失，确定平均尾流损失的变化范围，选择变动度最小的候选排布方案作为推荐排布方案。如图8所示，步骤S5具体包括以下步骤：

S51：建立风向频率变动模型；

S52：计算平均尾流损失；

S53：评价变动度指标；

S54:推荐风机排布方案。

具体地，记测风塔的风向频率分布为J(θi)，i＝1,2,…h，h为风向扇区数目；θ为扇区风向，单位为度；主导风向记为θp，相应频率为J(θp)。风向频率变动模型设置8种变动情景，模拟风向在年际、年代际尺度上的变化：

(1)J(θp)增加30％，其余扇区的频率相应平均减少；

(2)J(θp)减少30％，其余扇区的频率相应平均增多；

(3)θp左移一个扇区，J(θp)保持不变；

(4)θp右移一个扇区，J(θp)保持不变；

(5)θp左移一个扇区，J(θp)增加30％，其余扇区的频率相应平均减少；

(6)θp左移一个扇区，J(θp)减少30％，其余扇区的频率相应平均增多；

(7)θp右移一个扇区，J(θp)增加30％，其余扇区的频率相应平均减少；

(8)θp右移一个扇区，J(θp)减少30％，其余扇区的频率相应平均增多。

对8种情景重新计算平均尾流损失dWi′,i＝1,2,…8，形成尾流损失区间dWmin′～dWmax′；变动度指标定义为dWmax′/dWmin′。

相对于现有技术，本发明的海上风电场的风机排布方法，充分考虑风机排布在形态、方位、间距等三个方面的变化，通过定义风机排布基本单元，生成平行四边形和扇形阵列风机坐标，利用枚举法筛选出满足风电场约束条件的所有候选排布方案，引入尾流数学模型，对候选方案逐一进行计算，选出平均尾流损失最低的推荐排布方案，能够灵活适应多主导风向的情况与风电场的复杂几何边界，阵列筛选范围更广，效率更高。

本发明并不局限于上述实施方式，如果对本发明的各种改动或变形不脱离本发明的精神和范围，倘若这些改动和变形属于本发明的权利要求和等同技术范围之内，则本发明也意图包含这些改动和变形。

Claims (12)

1.一种海上风电场的风机排布方法，其特征在于：包括步骤：

S1：确定风机的排布阵列；

S2：获取风电场的约束条件数据；

S3：根据风机的排布阵列和风电场的约束条件筛选候选排布方案；

S4：计算各方案的尾流损失；

S5：检验最优方案。

2.根据权利要求1所述的海上风电场的风机排布方法，其特征在于：所述风机的排布阵列包括平行四边形和扇形两类。

3.根据权利要求2所述的海上风电场的风机排布方法，其特征在于：风电场的约束条件数据包括风机排布区域边界数据和风机的数目。

4.根据权利要求3所述的海上风电场的风机排布方法，其特征在于：所述步骤S3进一步包括步骤：

S31：确定西南、西北、东北、东南共4个固定点S1～S4；

S32：确定阵列排布的起始机位P0，坐标记为(X0，Y0)；

S32：确定阵列的行数n和列数m；

S33：计算阵列机位P(i,j)的坐标(Xi，Yj)，其中，i＝1,2,…n，j＝1,2,…m；

S34：剔除位于内、外边界之外的无效机位P″，有效机位P′(i,j)总数记为N′；

S35：循环控制，判断N′是否等于N，如果是，将P′列为候选排布方案，如果否，转到步骤S31。

5.根据权利要求4所述的海上风电场的风机排布方法，其特征在于：所述步骤S4具体包括步骤：

S41：计算单个扇区风速损失；

S42：计算单个扇区功率损失；

S43：尾流计算循环控制，判断是否全部扇区的尾流计算都已经完成，如果是，转到步骤S44，如果否，转到步骤S41；

S44：计算平均尾流损失；

S45：排布方案循环控制，判断是否全部排布方案都已经计算完平均尾流损失，如果是，转到步骤S46，如果否，转到步骤S41；

S46：获得平均尾流损失最低的N个方案，N>1。

6.根据权利要求1至5任一权利要求所述的海上风电场的风机排布方法，其特征在于：步骤S5进一步包括步骤：

S51：建立风向频率变动模型；

S52：计算平均尾流损失；

S53：评价变动度指标；

S54:推荐风机排布方案。

7.根据权利要求6所述的海上风电场的风机排布方法，其特征在于：所述步骤S4具体为：在0°至360°平均划分144个扇区对候选排布方案的逐台风机进行风速损失和功率损失计算，并用平均功率损失表示整个阵列的平均尾流损失。

8.根据权利要求7任一权利要求所述的海上风电场的风机排布方法，其特征在于：所述步骤S5具体为：选出平均尾流损失最低的5个候选排布方案，利用风向频率变动模型重新计算各自的尾流损失，确定平均尾流损失的变化范围，选择变动度最小的候选排布方案作为推荐排布方案。

9.根据权利要求8所述的海上风电场的风机排布方法，其特征在于：所述风机的排布阵列具体为：风机阵列的基本单元定义为风机排布方案中每4台风机的相对位置所形成的平面形态，

若基本单元为平行四边形，长对角线AC称为主对角线，短对角线BD称为次对角线，控制参数包括主对角线的方位角α，主、次对角线间的夹角β，以及主、次对角线的长度d1、d2；

若基本单元为扇形，短弧AD与长弧BC共圆，圆心为O；半径OA称为短径，OB称为长径，直线AB称为始边，直线CD称为终边，控制参数包括扇形始边的方位角α，始、终边间的夹角β，以及长、短径的长度d1、d2。

10.根据权利要求9所述的海上风电场的风机排布方法，其特征在于：所述步骤S3具体的计算方法为：

若风机阵列为平行四边形，P0坐标由两条虚拟线L1、L2相交确定，L1、L2的斜率记为k1、k2，

当0°≤α<90°，L1、L2分别经过S1、S4；当90°≤α<180°，L1、L2分别经过S2、S1；当180°≤α<270°，L1、L2分别经过S3、S2；当270°≤α<360°，L1、L2分别经过S4、S3，P0(X0，Y0)则根据解析几何中的两直线交点公式确定；

n、m由两组虚拟平行线L1、L′1和L2、L′2共同确定，当0°≤α<90°，L′1、L′2分别经过S3、S2；当90°≤α<180°，L′1、L′2分别经过S4、S3；当180°≤α<270°，L′1、L′2分别经过S1、S4；当270°≤α<360°，L′1、L′2分别经过S2、S1；记L1、L′1在东西方向上投影之间的距离为dX，记L2、L′2在南北方向上投影之间的距离为dY，平行四边形的边AB、边AD长度分别为x、y，可由α、β、d1、d2根据三角形余弦定理计算得到；

n、m表达式如下：

n＝dX/x

m＝dY/y

机位P_(i,j)的坐标(X_i，Y_j)表达式如下：

X_i＝X₀+(i–1)×x i＝1,2,…n，

Y_j＝Y₀+(j–1)×y j＝1,2,…m

若风机阵列为扇形，P₀坐标由S₁～S₄经过坐标旋转后确定，坐标旋转角为α+β/2。记旋转后的S₁～S₄为S′₁～S′₄，重新计算X′_min、X′_max、Y′_min、Y′_max。旋转后的(X₀，Y₀)记为(X′₀，Y′₀)，(X′₀，Y′₀)和(X₀，Y₀)表达式如下：

X′₀＝(X′_min+X′_max)/2

Y′₀＝Y′_min–d₂

X₀＝X′₀cos(α+β/2)+Y′₀sin(α+β/2)

Y₀＝X′₀cos(α+β/2)-Y′₀sin(α+β/2)

n、m表达式如下：

$n=\frac{\sqrt{{(X_{\max }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\max }^{\&prime;}-Y_0^{\&prime;})}^2}-\sqrt{{(X_{min}^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{min}^{\&prime;}-Y_0^{\&prime;})}^2}}{d_1-d_2}$

$m=\&lsqb;a\tan \left(\frac{X_{\max }^{\&prime;}-X^{\&prime;}}{\sqrt{{(X_{\max }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{\max }^{\&prime;}-Y_0^{\&prime;})}^2}}\right)-a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)\&rsqb;\frac{180}{\mathrm{\&pi;}\mathrm{\&beta;}}$

机位P_(i,j)的坐标(X_i，Y_j)表达式如下：

$\begin{array}{c}{X}_j=X_0+\&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\sin \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\&rsqb;\cos (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)-\\ \&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\sin (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)\end{array}$

$\begin{array}{c}{Y}_j=Y_0+\&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\&rsqb;\cos (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)-\\ \&lsqb;d_2+(i-1)(d_1-d_2)\&rsqb;\cos \&lsqb;a\tan \left(\frac{X_0^{\&prime;}-X_{\min }^{\&prime;}}{\sqrt{{(X_{\min }^{\&prime;}-X_0^{\&prime;})}^2+{(Y_{inx}^{\&prime;}-Y_0^{\&prime;})}^2}}\right)+(j-1)\mathrm{\&beta;}\sin (\mathrm{\&alpha;}+\mathrm{\&beta;}/2)\end{array}$

11.根据权利要求10所述的海上风电场的风机排布方法，其特征在于：所述步骤S4具体计算方法为：单台风机在单个扇区所受实际风速损失dv′为：

${\mathrm{dv}}^{\&prime;}=\sqrt{\stackrel{n}{\mathrm{\&Sigma;}}{\left(dv\right)}^2}$

$dv=\mathrm{\&phi;}\&lsqb;1-\sqrt{1-C_T}\&rsqb;{\left(\frac{R}{R+kL}\right)}^2$

其中，n为上风向尾流个数；CT为风机推力系数；R为风机叶轮半径；L为风机间距在风向上的投影距离；k为尾流衰减系数，Φ为尾流遮蔽系数，表明风机被上风向尾流遮蔽的程度，

单台风机在单个扇区所受实际功率损失dP′为：

${\mathrm{dP}}^{\&prime;}={\&Integral;}_{v_s}^{v_e}g\left(v\right)f\left(v\right)dv-{\&Integral;}_{v_s}^{v_e}g\left(v\right)f\&lsqb;v(1-{\mathrm{dv}}^{\&prime;})\&rsqb;dv$

其中，v_s、v_e为风机切入和切出风速；f(v)为风机功率曲线，g(v)为风速在该扇区的概率密度分布，

平均尾流损失dW′为：

${\mathrm{dW}}^{\&prime;}=\frac{\underset{j=1}{\overset{144}{\&Sigma;}}{\mathrm{\&omega;}}_j\&CenterDot;\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{dP}}_i^{\&prime;}}{N}$

其中，N为风机数目；ω为各个风向扇区的权重，根据海上风电场测风塔在12、16或48个方位的风向频率分布结果平均分配。

12.根据权利要求11所述的海上风电场的风机排布方法，其特征在于：步骤S5具体为：记测风塔的风向频率分布为J(θi)，i＝1,2,…h，h为风向扇区数目；θ为扇区风向，单位为度；主导风向记为θp，相应频率为J(θp)。风向频率变动模型设置8种变动情景，模拟风向在年际、年代际尺度上的变化：

(1)J(θp)增加30％，其余扇区的频率相应平均减少；

(2)J(θp)减少30％，其余扇区的频率相应平均增多；

(3)θp左移一个扇区，J(θp)保持不变；

(4)θp右移一个扇区，J(θp)保持不变；

(5)θp左移一个扇区，J(θp)增加30％，其余扇区的频率相应平均减少；

(6)θp左移一个扇区，J(θp)减少30％，其余扇区的频率相应平均增多；

(7)θp右移一个扇区，J(θp)增加30％，其余扇区的频率相应平均减少；

(8)θp右移一个扇区，J(θp)减少30％，其余扇区的频率相应平均增多。

对以上8种情景重新计算平均尾流损失dWi′,i＝1,2,…8，形成尾流损失区间dWmin′～dWmax′；变动度指标定义为dWmax′/dWmin′。