基于致动线方法的风力机组风场布置优化

王渊博， 李 春,2， 缪维跑

(1.上海理工大学 能源与动力工程学院，上海 200093；2.上海理工大学 能源与动力工程学院，上海市动力工程多相流动与传热重点实验室，上海 200093)

基于致动线方法的风力机组风场布置优化

王渊博1， 李 春1,2， 缪维跑1

(1.上海理工大学 能源与动力工程学院，上海 200093；2.上海理工大学 能源与动力工程学院，上海市动力工程多相流动与传热重点实验室，上海 200093)

基于致动线方法，利用OpenFOAM对5种不同间距串列布置的具有3台5 MW风力机的风场进行了模拟.通过调整风力机组间的位置，对比不同风场布置时的总输出功率，并结合轮毂高度处流场分布及尾迹速度轮廓线，分析尾迹对风场下游风力机的影响.结果表明：串列布置时，处于尾迹中的风力机功率下降可达82.72%；改变风场中风力机组相对位置可使风场布置得到优化，虽个别风力机性能有所下降，但风场总输出功率最大可提高11.8%.

风场布置； 总输出功率； 致动线； 风力机； OpenFOAM

风能在世界各国的能源供给中发挥着十分重要的作用[1-3].风场上游风力机产生的尾迹会导致下游速度亏损及湍流强度增加[4-5]，不仅大幅降低了下游风力机的输出功率，浪费风场风和土地资源，更为重要的是还会导致叶片疲劳载荷显著提高，降低风力机寿命周期，甚至导致叶片损坏事故[6].因此，优化风场布置以减小风力机组间尾迹的影响显得尤为重要[7-8].

目前，获取风力机气动性能的主要方法有叶素动量理论(BEM)方法、涡尾迹法和求解N-S方程的计算流体力学(CFD)方法.在上述方法中，BEM方法尽管计算效率高，但不能计算风力机尾迹流场，尚无法衡量风场上游风力机尾迹对下游风力机气动性能的影响[9].与BEM方法相比，涡尾迹法有相对高的精度，且可以预测风力机尾迹，但随风力机尺寸的增大，风力机绕流场涡结构跨越多个尺度，涡尾迹法对大尺度空间内小涡发展的预测能力较低[10].因此，风力机气动性能研究逐渐转向更加精确的CFD方法[11]，即采用相关数值模拟方法求解N-S方程以获取流场变量.铁庚等[12]利用Fluent软件对FFA-W3-211、FFA-W3-301和NA-CA63-215 3种翼型进行数值模拟，并与实验数据进行对比，验证了二维黏性不可压数值模拟的可靠性.杨贺[13]利用Fluent软件，基于10 kW水平轴风力机进行了叶片外形设计和风轮气动性能数值模拟，最终得到叶片各叶素参数和风轮转矩，结果验证了设计的合理性.张昇龙[14]利用Fluent软件对某商用1.5 MW水平轴风力机在3种不同风速条件下的风轮尾迹和气动性能进行数值模拟，将计算结果与相似实验测量值进行比较，证明了三维黏性不可压模型的正确性.贾彦等[15]利用WAsP软件及风资源数据进行风场模拟计算，建立了由多台风力机组成的风场，比较不同风力机组布局下风场的年净发电量，并给出了风场最优风力机组布局.上述研究中，铁庚等[12-14]仅对风力机叶片或不超过2 MW的单台风力机进行CFD数值模拟，这是因为CFD方法需建立大量近壁面网格以求解边界层流动，模拟多台兆瓦级风力机流场需耗费巨大的计算资源，计算成本难以承受.贾彦等[15]虽对多台风力机风场进行模拟，但WAsP软件基于线性流动控制方程及半经验的尾迹模型，且假设流动不发生分离，计算结果与真实风场存在一定误差，同时无法获得风力机周围流场信息[16].

Sorensen等[17]提出致动线方法，结合了BEM方法计算速度快的优点与CFD方法计算精度高的优点，通过较少的计算资源可获得很好的流场信息，尤其适用于风力机组三维流动模拟.笔者基于致动线方法，通过开源CFD软件OpenFOAM模拟了5种不同间距串列布置的3台5 MW风力机绕流流场，研究风电机组间距对风电场整体性能的影响，通过分析风力机流场特征和尾迹速度分布，揭示上游尾迹干扰下游风力机吸收风能的流动机理，为未来大型风电场优化布局提供参考.

1 致动线方法

1.1 致动线方法控制方程

致动线方法控制方程在笛卡尔坐标系下的张量形式为：

(1)

式中：下标i、j使用了爱因斯坦求和约定；t为时间；xi、xj为笛卡尔坐标系的坐标轴；vi、vj为速度沿坐标轴方向的分量；ρ为密度；p为压强；μ为动力黏度；源项fi为体积力，表征叶片对流体的反作用力，为致动线方法控制方程特有，传统CFD方法控制方程忽略此项.

1.2 体积力

1.2.1 体积力计算

根据牛顿第三定律，可利用流场对叶片的作用力求得叶片对流场的体积力.由于致动线方法中无实体壁面，需通过其他模型获取叶片的作用力以达到与实体叶片相同的效果.以BEM方法作为附加模型，在每个时间步读取流场信息，根据翼型气动特性计算叶片的作用力，计算公式如下：

(2)

(3)

(4)

式中：W为当地翼型的相对速度；c为当地翼型弦长；Cτ为切向力系数；τ为切向单位向量；Cn为法向力系数；n为法向单位向量；Fs为当地翼展单位长度所受到的作用力.

当地翼展单位长度叶片对流场的体积力为:

(5)

1.2.2 体积力投射

叶片对流场的体积力需投射至致动线附近的计算域网格，通过致动线方法控制方程式(1)求解下一时间步的流场.体积力投射规律对模拟结果的影响极大，当投射范围较小时，数值计算将产生严重振荡，甚至无法收敛；当投射范围过大时，投射区域网格所属体积力过小，流场状态改变微弱，尤其是叶尖涡和叶根涡被过分抹平，与真实情况不相符[18].

为解决上述问题，将体积力以三维高斯分布向流场投射[17]，表达式如下：

(6)

式中：r为体积力中心至投射点的距离；ε为分布因子，控制体积力分布的集中程度.

Troldborg[19]对致动线方法中体积力的分布形式进行深入研究后指出，若采用高斯分布向流场投射体积力，分布因子ε=2Δs时的效果与真实流场最接近，其中Δs为致动线分布点间距.此外，为提高计算效率，同时兼顾精确度，需确定投射边界，边界上投射的体积力大小为最大体积力的5%，边界之外不再进行体积力投射[20].

2 构造致动线求解器

PisoFoam是OpenFOAM自带的标准求解器，可以求解湍流及层流流动，基于瞬态Piso算法编写，相比Simple算法，该算法计算效率更高[21].图1为瞬态Piso算法流程图.

图1 瞬态Piso算法流程图

基于PisoFoam求解器构造新求解器PisoFoamTurbine，专用于基于致动线方法对风力机的模拟.与标准求解器PisoFoam相比，新求解器PisoFoamTurbine建立了“turbines”类，该类包括定义体积力数组、求解体积力大小和投射体积力等一系列操作.计算过程中，通过该类求解相应范围内体积力的大小，改变动量方程，并在式(1)源项中添加体积力.

PisoFoamTurbine求解器中动量方程变为以下形式：

fvVectorMatrix UEqn

(

fvm::ddt(U)

+ fvm::div(phi, U)

+ turbulence->divDevReff(U)

- turbines.force()

);

其中,turbines.force()函数是用来实现流场中体积力计算及完成相应空间分布的操作，其余计算均与PisoFoam一致.

3 计算模型与计算域网格划分

3.1 计算模型

计算模型选用美国可再生能源实验室(NREL)设计的5 MW风力机，主要参数如表1所示.

表1 NREL 5 MW风力机主要参数

3.2 风力机组布置

以D表示风轮直径大小，在27D×6D×6D的模拟风场中串列布置3台NREL 5 MW风力机，分别为WT1、WT2和WT3，如图2所示.WT1距入口3D；WT3与WT1的距离为固定值14D；WT2与WT1的距离L为变量，取值如表2所示.由于尾迹长度小于5D时湍流强度急剧增强[22-23]，风力机疲劳载荷增大、寿命缩短，故L最小取5D以保证WT2不在强尾迹区，最大取9D以保证WT3不在强尾迹区.

图2 模拟风场中风力机组的布置

算例编号12345L5D6D7D8D9D

3.3 数值计算模型

Sarmast等[24]详细研究了分布因子及网格尺度对致动线方法计算结果的影响：当分布因子确定时，优化网格可加速计算结果收敛并提高准确性.因此，对模拟风场中风轮所在水平圆柱区域的网格进行加密，并在网格加密区与粗网格区之间添加2层网格过渡带,以避免数值计算振荡，如图3所示.网格尺度从内向外逐层加倍，其中加密区网格尺度为3.94 m×3.94 m×3.94 m，与Fleming等[25]模拟风场的网格尺度相近.

(a) 主视图

(b) 三维透视图

模拟风场的边界条件类型如图4所示.进口设为速度入口，方向垂直于进口平面，大小为11.4 m/s；出口设为压力出口，大小为标准大气压；地面和顶面设为滑移边界条件；两个侧面设为周期边界条件，即计算域一个侧面网格节点与另一个侧面网格节点一一对应，紧邻周期边界一侧计算域外的“镜像单元”信息由紧邻另一侧的计算域内的单元提供.

图4 模拟风场的边界条件类型

数值计算采用大涡模拟(LES)方法，流动变量通过网格过滤操作分为网格尺度(GS)和亚格子尺度(SGS)2部分.引入各向异性的残余应力张量模型封闭亚格子模型，利用Smagorinsky涡粘假设参数化，将模型分为2部分：线性涡粘模型和基于混合长假设的亚尺度涡粘模型[26].

4 结果与分析

通过致动线方法在每个时间步求出风力机功率并实时输出.当风从计算域入口发展到出口之后再计算200 s，风力机的输出功率即取此时间段内的平均值.图5为不同L处各风力机输出功率的数值模拟结果.

图5 风力机的输出功率

从图5可以看出，当L分别为5D、6D、7D、8D和9D时，风场中3台风力机总输出功率分别为7.719 1 MW、7.477 4 MW、7.321 6 MW、7.014 8 MW和6.808 2 MW，即随L的增大总输出功率逐渐减小.通过最小二乘法拟合总输出功率与L的关系，可得到式(7)：

P=-0.228 5l+8.867 3

(7)

式中：P为风场总输出功率，MW；l为L中的数值部分.

随着L的变化，WT2和WT3的输出功率变化明显.理论上，下游风场中的风力机要达到与上游风力机相同的输出功率，则需与上游风场风力机保持足够大的距离[27].WT1的输出功率随着L的增大稍微有所增大，证明下游风力机对WT1也有轻微影响，且随着L增大此影响逐渐减弱.WT1输出功率变化不大的原因是其处于风场上游且入流风速不变.

WT2和WT3的输出功率比见表3.输出功率比定义为WT2和WT3的输出功率分别与WT1输出功率的比值.风力机尾迹的湍流强度具有非线性特征[28]，所以WT2和WT3的输出功率比与L呈现非线性关系.由表3可见，WT2的输出功率没有超过WT1输出功率的21%，此现象在顺排风场中较为常见[29-30].当下游风力机位于上游风力机尾迹区域时，尾迹的速度亏损和湍流强度增加造成下游风力机入流风速减小且极不稳定，导致下游风力机输出功率严重下降.因为WT3所处的尾迹由WT2和WT1的尾迹叠加而成，这种叠加尾迹的不稳定性比非叠加尾迹高，可促进其与周围流体的混合和能量交换，加速尾迹恢复，使尾迹速度亏损得以快速补偿，湍流强度大幅减小，为下游风力机提供质量相对较高的入流风，所以WT3的输出功率整体上比WT2的高，这也证明了所研究内容的准确性及可靠性.

表3 WT2和WT3的输出功率比

通过最小二乘法分别拟合WT2和WT3输出功率比与L的关系，可得到如式(8)和式(9)所示的关系式.

R2=0.008 7l+0.130 1

(8)

R3=-0.055 4l+0.662 2

(9)

式中：R2为WT2的输出功率比；R3为WT3的输出功率比.

借助速度云图可更直观地看出L对风力机组的影响，图6为不同L处风力机轮毂高度处的速度云图.从图6可以看出，随着L的增大，WT2的低速尾迹区域逐渐增大.当L=5D和L=6D时，WT2与WT3间距较大，加之WT1尾迹的叠加，WT2尾迹与周围流体发生了较为充分的混合和能量交换，因此接近WT3处的尾迹速度得到一定提升，WT3入流速度较大，其输出功率较WT2的高.当L=7D时，WT2的低速尾迹发展至WT3风轮平面，导致WT3入流速度较低，是WT3输出功率降低的主要原因.随着L的进一步增大，当L=9D时，WT2与WT3间距较小，WT2尾迹尚未与周围流体发生足够的混合和能量交换，尾迹速度亏损严重，极大地降低了WT3的输出功率；与此同时WT2与WT1间距较大，WT2输出功率受WT1尾迹影响较小，WT2输出功率已高于WT3输出功率.上述分析进一步说明了文中计算的可行性与实用性.

图6 轮毂高度处速度云图

此外，从图6可以看出，WT3尾迹比WT2尾迹和WT1尾迹更为复杂，非稳定流动明显,因此与周围流体的混合和能量交换速度更充分，使得尾迹得到快速恢复.

从图5可知，风场总功率的变化主要由WT3引起，因此有必要进一步分析WT3的尾迹.图7为不同L处WT3各风轮直径轮毂高度处的速度轮廓线.从图7可以看出，L取值不同时，WT3各风轮直径轮毂高度处的速度轮廓线变化较大.L=9D时，WT3与WT2间距最小，WT2的尾迹还未充分恢复，导致WT3的入流风速非常低，此时速度轮廓线凹的幅度最大，WT3输出功率相应最小.L=5D时，速度轮廓线相对比较平缓，说明WT3所处尾迹速度亏损较小.此时,WT3距WT2较远，WT2的尾迹有足够的时间和空间与周围流体进行混合和能量交换，尾迹恢复充分.L分别为7D、8D和9D时，速度轮廓线逐步向上抬升，WT3的入流风速和输出功率相应逐渐增大.由此可以看出所得结论的统一性和正确性.此外，从图7还可以看出，尾迹附近区域风速略高于入流风速，这是风场上游风力机造成的风加速现象.

图7 速度轮廓线

5 结 论

(1) 风力机尾迹处会出现速度亏损和湍流强度增加的现象，处于尾迹中的风力机，其入流风速相应减小且极不稳定，输出功率大幅降低.

(2) 叠加尾迹的非稳定流动更加明显，与周围流体混合得更加充分，恢复更加完全.

(3) 适当调整风场中风力机的相对位置，可以优化风场布置，从而大幅提高风场总输出功率.文中模拟风场WT2与WT1的距离为L=5D时,虽然WT2的输出功率较WT1减小82.72%，但风场总输出功率提高了11.8%.

[1] 李少华, 匡青峰, 吴殿文, 等. 1.2 MW风力机整机流场的数值模拟[J]. 动力工程学报, 2011,31(7): 551-556.

LI Shaohua, KUANG Qingfeng, WU Dianwen, et al. Numerical simulation on flow field of a 1.2 MW wind turbine[J]. Journal of Chinese Society of Power Engineering, 2011, 31(7): 551-557.

[2] SAHU B K, HILOIDHARI M, BARUAH D C. Global trend in wind power with special focus on the top five wind power producing countries[J]. Renewable and Sustainable Energy Reviews, 2013, 19(1): 348-359.

[3] CHENG Shikun, LI Zifu, MANG H P, et al. A review of prefabricated biogas digesters in China[J]. Renewable and Sustainable Energy Reviews, 2013, 28: 738-748.

[4] 韩中合, 李引, 季剑. 考虑风剪切的1.3 MW风力机整机三维定常流动数值研究[J]. 动力工程学报, 2011, 31(10): 779-783.

HAN Zhonghe, LI Yin, JI Jian. Numerical study on 3D steady flow of a 1.3 MW wind turbine considering wind shear factor[J]. Journal of Chinese Society of Power Engineering, 2011, 31(10): 779-783.

[5] 黎作武, 贺德馨. 风能工程中流体力学问题的研究现状与进展[J]. 力学进展, 2013, 43(5): 472-525.

LI Zuowu, HE Dexin. Reviews of fluid dynamics researches in wind energy engineering[J]. Advances in Mechanics, 2013, 43(5): 472-525.

[6] 王强. 水平轴风力机三维空气动力学计算模型研究[D]. 北京: 中国科学院大学, 2014.

[7] OZTURK U A, NORMAN B A. Heuristic methods for wind energy conversion system positioning[J]. Electric Power Systems Research, 2004, 70(3): 179-185.

[8] INDASI V S, LYNCH M, MCGANN B, et al. WAsP model performance verification using lidar data[J]. International Journal of Energy and Environmental Engineering, 2016, 7(1): 105-113.

[9] SHENKAR R. Design and optimization of planar and nonplanar wind turbine blades using vortex methods[D]. Denmark: Technical University of Denmark, 2010.

[10] 夏巍巍. 水平轴风力机尾流结构及其对风电场布局影响的研究[D]. 兰州: 兰州理工大学, 2014.

[11] MIAO Weipao, LI Chun, JUN Yang, et al. Numerical investigation of wake control strategies for maximizing the power generation of wind farm[J]. Journal of Solar Energy Engineering, 2016, 138(3): 034501.

[12] 铁庚, 祁文军. 水平轴风力机叶片翼型的气动特性数值模拟[J]. 流体机械, 2013, 41(3): 29-33.

TIE Geng, QI Wenjun. Numerical simulation on the aerodynamic performance of aerofoil of horizontal wind turbine blade[J]. Fluid Machinery, 2013, 41(3): 29-33.

[13] 杨贺. 水平轴风力机风轮气动性能研究[D]. 呼和浩特: 内蒙古工业大学, 2014.

[14] 张昇龙. 水平轴风力机尾流场及其气动性能的研究[D]. 兰州: 兰州理工大学, 2014.

[15] 贾彦, 刘璇, 李华, 等. 考虑尾流效应对风电场机组布局的影响分析[J]. 可再生能源, 2014, 32(4): 429-435.

JIA Yan, LIU Xuan, LI Hua, et al. Analysis of wind farm units layout considering wake effect[J]. Renewable Energy Resources, 2014, 32(4): 429-435.

[16] 田琳琳, 赵宁, 武从海, 等. 复杂地形风电场的机组布局优化[J]. 南京航空航天大学学报, 2013, 45(4): 503-509.

TIAN Linlin, ZHAO Ning, WU Conghai, et al. Optimal positioning of wind turbines on complex wind farm[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2013, 45(4): 503-509.

[17] SORENSEN J N, SHEN Wenzhong. Numerical modeling of wind turbine wakes[J]. Journal of Fluids Engineering, 2002, 124(2): 393-399.

[18] 刘磊. 风力机叶片非定常气动特性的研究[D]. 北京: 中国科学院大学, 2012.

[19] TROLDBORG N. Actuator line modeling of wind turbine wakes[D]. Denmark: Technical University of Denmark, 2008.

[20] JIN Wenjie. Numerical simulation of wind turbine wakes based on actuator line method in NEK5000[D]. Stockholm, Sweden: Royal Institute of Technology, 2013.

[21] 李东岳. 浅谈Piso，Simple算法[EB/OL]. (2016-03-12) [2016-05-10]. http://www.dyfluid.Com/PISOSIMPLE.html.

[22] MO J O, CHOUDHRY A, ARJOMANDI M, et al. Large eddy simulation of the wind turbine wake characteristics in the numerical wind tunnel model[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2013, 112: 11-24.

[23] WILSON J M, DAVIS C J, VENAYAGAMOORTHY S K, et al. Comparisons of horizontal-axis wind turbine wake interaction models[J]. Journal of Solar Energy Engineering, 2015, 137(3): 031001.

[24] SARMAST S, SEGALINI A, MIKKELSEN R F, et al. Comparison of the near-wake between actuator-line simulations and a simplified vortex model of a horizontal-axis wind turbine[J]. Wind Energy, 2016, 19(3): 471-481.

[25] FLEMING P, GEBRAAD P M O, LEE S, et al. Simulation comparison of wake mitigation control strategies for a two-turbine case[J]. Wind Energy, 2014, 18(2): 2135-3143.

[26] 卞凤娇, 徐宇, 王强, 等. 基于OpenFOAM的风力机致动线模型研究[J]. 工程热物理学报, 2016, 37(1): 72-75.

BIAN Fengjiao, XU Yu, WANG Qiang, et al. Numerical study of actuator line model of wind turbine based on OpenFOAM platform[J]. Journal of Engineering Thermophysics, 2016, 37(1): 72-75.

[27] CHOWDHURY S, ZHANG Jie, MESSAC A, et al. Unrestricted wind farm layout optimization (UWFLO): investigating key factors influencing the maximum power generation[J]. Renewable Energy, 2012, 38(1): 16-30.

[28] SHAKOOR R, HASSAN M Y, RAHEEM A, et al. Wake effect modeling: a review of wind farm layout optimization using Jensen's model[J]. Renewable and Sustainable Energy Reviews, 2016, 58: 1048-1059.

[29] WU Yuting, PORTÉ-AGEL F. Simulation of turbulent flow inside and above wind farms: model validation and layout effects[J]. Boundary-Layer Meteorology, 2012, 146(2): 181-205.

[30] NEWMAN J, LEBRON J, MENEVEAU C, et al. Streamwise development of the wind turbine boundary layer over a model wind turbine array[J]. Physics of Fluids, 2013, 25(8): 231-234.

Layout Optimization of the Wind Farm in a Wind Power Station Based on Actuator Line Method

WANGYuanbo1,LIChun1,2,MIAOWeipao1

(1. School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China; 2. Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China)

Based on actuator line method, numerical simulations were conducted on the wind farm with three 5 MW wind turbines in five different aligned arrangements using OpenFOAM. By comparing the total power output, the velocity fields and velocity profiles at the hub height of wind turbines for different wind farm layouts, the effects of upstream turbine wake on the downstream wind turbine were investigated. Results show that in an aligned arrangement, the power output may be reduced by 82.72% for the wind turbine lying in the wake. Changing the relative position of wind turbines reasonably could help to optimize the wind farm layout, when the total power output may be increased by a maximum of 11.8%, even if the output of an individual wind turbine declines somewhat.

wind farm layout; total power output; actuator line; wind turbine; OpenFOAM

2016-05-26

2016-08-04

国家自然科学基金资助项目(51176129，51676131)；上海市科学技术委员会资助项目(13DZ2260900)

王渊博(1991-)，男，陕西澄城人，硕士研究生，研究方向为风力发电. 李 春(通信作者)，男，教授，博导，电话(Tel.)：15216702797；E-mail：lichun_usst@163.com．

1674-7607(2017)05-0418-07

TK83

A 学科分类号：480.60