单一水平轴风电机组尾迹的模拟方法与流动机理研究综述1)

杨晓雷

(中国科学院力学研究所非线性力学国家重点实验室,北京 100190)

(中国科学院大学工程科学学院,北京 100049)

引言

化石能源的使用导致地球气温升高,给环境带来了巨大甚至不可逆的影响,比如极端天气、冰川融化等[1].根据国家发改委能源研究所发布的《2020 年中国可再生能源展望报告》,到2050 年,中国风电将占到能源消费的38.5%[2].在风电占如此大的比重下,其成本需进一步降低以保持经济上的竞争力[3].从19 世纪末到现在,基于空气动力学理论的叶片设计技术不断进步,促进风力机发电效率不断提高.同时,通过采用更高的塔筒和更长的叶片,使得单个机组的发电量不断提升,从千瓦到十兆瓦甚至更高[3].另一方面,风电场中的风力机通常按照一定间距排布成阵列,其性能受风力机尾迹影响[3].风力机尾迹影响发电量、风电质量、及风电机组维护费用.风机尾迹风速低,带来的平均发电量损失约为20%[4],高可达80%[5],湍流强度高,影响下游风力机所受疲劳载荷,从而增加维护费用[6-7],是影响风电成本的关键因素,尾迹湍流脉造成下游机组电量产出脉动,影响风电并网性能.

当前风力机控制和设计以最优自身性能为目标,然而这通常导致整个风电场性能不佳.将风电场作为一个整体进行设计和控制优化有潜力显著提高风电场电量产出,降低运营与维护成本.对风力机尾迹机理进行深入研究,发展尾迹快速预测模型,是实现风电场整体设计和控制的流体力学基础.随着计算机计算能力的不断提升,数值模拟在尾迹机理研究和模型发展中发挥越来越重要的作用.

本文介绍风力机尾迹的数值模拟方法和机理研究进展,分为3 部分,第1 节尾迹模拟方法,第2 节尾迹机理,第3 节结语与展望.

1 尾迹模拟方法

风力机尾迹模拟需考虑多个尺度,从叶片表面边界层(～ 10-2m)到风轮直径(D～ 102m)再到大气边界层厚度(～ 103m)及局部气象尺度(105m)[8-9],是一个典型的多尺度问题.风力机尾迹模型大致可分有两类: (1)行业广泛采用的工程模型,主要用于风力机布置和控制;(2)精度高、但计算量大的高可信度模型,主要用于机理研究,为发展工程模型提供理论依据.

1.1 解析模型与低阶模型

Jensen 模型是常见的尾迹解析模型[10],其形式如下

其中a为轴向诱导因子,U0为来流风速,k为夹带常数k=0.5/ln(zh/z0) (其中,zh为轮毂高度,z0为地面粗糙度长度),r1为初始尾迹宽度,通过一维动量理论给定,需要说明的是,式(1)来自Frandsen 文章[11],其夹带常数和初始尾迹宽度的计算方式与Jensen[10]有所不同.后续工作对文献[11]的模型作了进一步改进[12-13].考虑速度分布二维效应的Jensen 模型可以参考文献[14].

类似Jensen 模型的尾迹解析模型可计算尾迹速度亏损,但无法考虑尾迹与大气边界层的相互作用.另一方面,将风电场视为等效粗糙度长度的模型,可模化风力机阵列对大气边界层的影响及风电场内水平方向的平均速度.常见风电场等效粗糙度长度模型包括双层对数模型[11]、三层对数模型[15-16]及考虑流向和展向风力机间距不同作用的模型等[17].结合等效粗糙度长度模型和尾迹解析模型有望更好预测尾迹在大气边界层中的演化特征.为在统一框架下模拟风电场的发展区域和充分发展区域,Frandsen 等[18]将风电场分为尾迹相互独立的区域、尾迹发生相互影响的区域以及风电场与大气边界层相平衡的区域,联合尾迹解析模型、尾迹相互作用模型和等效粗糙度长度模型模化.该模型被Rathmann 等[19]进一步完善.Yang 和Sotiropoulos[20]利用内边界层概念,结合适合风电场的尾迹解析模型[17]和Frandsen[11]等效粗糙度长度模型发展了可预测任意分布和大小风电场功率的耦合模型.Stevens 等[21]通过迭代方式确定模型常数,发展并验证了耦合尾迹解析模型和等效粗糙度长度模型的耦合模型.Zhang 等[22]基于大涡模拟结果发展新的耦合模型,可预测不同风力机排布的等效粗糙度长度.

上述解析模型可快速预测不同位置的速度亏损,但预测精度依赖于参数选取,无法解析更多的物理.另一类方法通过求解简化Navier-Stokes 方程[23],计算尾迹速度分布.薄边界层方程是风力机尾迹模拟常用的简化方程,其形式如下

其中U和V为流向和径向的平均速度,u′和v′为速度脉动,雷诺应力可通过涡黏模型模化[24],μt∂U/∂r(μt为涡黏系数).

注意到上述模型只能预测风力机尾迹的时间平均特性,无法预测尾迹的时间和空间脉动特性,比如,远尾迹蜿蜒.丹麦技术大学的动态尾迹蜿蜒模型(dynamic wake meandering model)[25-26]假定来流大尺度涡是尾迹蜿蜒的主要成因,进一步通过泰勒流动冻结假设[27],将蜿蜒模拟成随来流大尺度运动的被动标量.近年来,机器学习在流体力学领域被广泛应用[28],人工神经网络也被用于预测尾迹蜿蜒,并取得较好结果[29].

1.2 大涡模拟和风力机参数化模型

相比于解析模型和基于简化Navier-Stokes 方程的模型,更高可信度的风力机尾迹计算方法包括离散涡方法[30]、雷诺平均方法 (RANS)[31-32]和 大涡模拟方法 (LES)[33].由于风力机尾迹的雷诺数很高,直接数值模拟(DNS)所需计算量极大,很难开展风力机尾迹的直接数值模拟.离散涡方法和雷诺平均方法具有计算效率高的优点,但无法准确计算尾迹湍流脉动.大涡模拟方法直接模拟湍流含能尺度,模化未解析小尺度,可较好捕捉风力机尾迹湍流脉动,目前广泛用于风力机尾迹模拟[34].由于所涉及流动尺度跨度大,难以直接解析所有尺度,风力机尾迹的大涡模拟通常采用参数化模型模化风力机与来流的相互作用,常用控制方程为不可压Navier-Stokes 方程

其中,xi(i=1,2,3) 为笛卡尔坐标,ui(i=1,2,3) 为滤波后的速度,ρ 为密度,p为压力,ν 为运动黏性系数,vt为湍流黏性系数,可采用经典Smagorinsky 模型[35-36]和动态Smagorinsky 模型[37]等大涡模拟亚格子模型模化,fi(i=1,2,3) 为风力机参数化模型引入的体积力.1.2.1 节将简要介绍常用风力机参数化模型.来流湍流对风力机尾迹特性有关键影响,1.2.2节将介绍风力机尾迹模拟中常采用的来流湍流生成方法.

1.2.1 风力机参数化模型

常用的风力机参数化模型大致可分为3 类,致动盘模型、致动线模型和致动面模型(如图1 所示).其中,致动盘和致动面的表面采用三角形网格离散;致动线采用分布点离散

图1 风力机参数化模型示意图Fig.1 Schematic for parameterization models for wind turbines.The surface of actuator disks and actuator surfaces are discretized using triangle cells;actuator lines are discretized using distributed points

致动盘模型,顾名思义,将整个风轮模化成可穿透的圆盘,其对来流的作用通过分布体积力表征[15].体积力分为轴向的推力和沿着转动方向的切向力.大部分致动盘模拟只考虑沿轴向的推力,可采用如下方式确定

其中,CT为推力系数,可通过一维动量理论给定CT=4a(1-a),A为风轮面积A=πR2(R为风轮半径),U0为该风力机的来流风速.对于风电场中的风力机尾迹模拟,来流风速U0不方便直接定义,可通过一维动量理论中来流和风轮面上的风速U1的关系确定U0=U1/(1-a) .后续工作中,人们在致动盘模型中考虑了转动的作用,发现可以提高尾迹预测精度[38].针对单一风力机和风力机阵列,Yang 和Sotiropoulos[39]将致动盘模型的预测结果和实验结果进行了对比,发现致动盘模型结果和实验结果的主要差别出现在近尾迹,而在远尾迹较好吻合.Li 和Yang[40]针对均匀来流和湍流来流两种情形,将致动盘模型预测结果与致动面结果进行了对比,发现在湍流来流时,两者在远尾迹差别较小,在均匀来流时,速度亏损、湍动能和雷诺应力在近尾迹和远尾迹都有明显差别.

致动线模型将风力机叶片模化为转动的线,不同位置叶片对来流的作用,通过致动线上的分布力表征[33,41].相比致动盘模型,致动线模型可以预测叶尖涡和中心涡等特征,但捕捉这些流动结构需要更高的空间解析度.在致动线模型中,参数化对象包括叶片不同径向位置的弦长c、扭角、翼型类型等几何特征,以及各个翼型在不同雷诺数、不同攻角下的升阻力系数(CL,CD)等空气动力学特征.根据来流风速、风轮转速及该径向位置的扭角,可以确定相对来流速度Urel和攻角,进而通过以下公式计算该径向位置的升阻力

致动面模型进一步将风力机叶片模化为转动的面,从而直接解析弦长方向几何的作用,减少需要模化的参数,提高对叶片几何的刻画精度.致动面模型首先由Shen 等[42]提出.该致动面模型通过预先得到的表面压力分布确定致动面上的体积力,可以较为准确地模化弦长方向力的分布特征.但这通常需要弦长方向被一定数量的网格点解析,同时,需要不同攻角、不同雷诺数的压力分布(需要通过实验或模拟确定),这使得模化弦长方向压力分布的致动面方法难以用于实尺度风电场的风力机尾迹模拟(主要有两方面困难: 解析弦长方向力分布特征的网格量巨大;叶片径向的翼型类型分布难以获得).Yang 和Sotiropoulos[9]提出基于叶素理论的致动面方法,和致动线方法采用同样的方式计算不同径向位置的力,进一步再将得到的力均匀分布到弦长方向.相比于致动线方法,该致动面方法可在一定程度上反映弦长方向的几何特征.同时,Yang 和Sotiropoulos 提出了机舱的致动面方法.类似于浸没边界方法,该方法采用法向无穿透条件计算机舱致动面上的法向力,通过来流速度和指定的摩擦力系数计算切向力.该致动面方法已在不同算例中得以验证,相比于致动线方法,可更准确预测远尾迹蜿蜒.Liao 等[43]进一步发展了适合螺旋桨的致动面模型,采用RANS计算叶片表面力,并将得到的力系数用于螺旋桨尾迹的大涡模拟.

在风力机参数化模型中,需要将得到的力分布到流场求解的背景网格.力的分布通常采用高斯函数实现,但这需要分布到周围较多(＞10 个网格宽度)网格才能保证在力在分布过程中力和力矩的守恒.Yang 和Sotiropoulos[39]采用浸没边界方法中采用的离散Delta 函数进行力的分布,可以在很少数目(2,3,4 或5 个网格宽度)的网格上满足力和力矩的守恒.为满足力和力矩守恒,离散Delta 函数需满足一定矩条件.对于动边界问题,离散Delta 函数导数也需满足一定矩条件.Yang 等[44]发展的离散Delta 函数,其本身和导数都满足相应矩条件,其中,4 个网格宽度的离散Delta 函数的形式如下

其中r=(x-X)/h,x和X分别为流场求解的网格坐标和参数化模型的网格坐标.

需要注意的是,该文着重论述风轮尾迹及机舱影响,未考虑塔筒尾迹及其模型.关于塔筒及其尾迹的相关研究,读者可以参考文献[45-47].

1.2.2 来流湍流生成

来流湍流影响风力机尾迹.生成接近真实大气边界层环境的来流湍流对风力机尾迹的大涡模拟至关重要.风力机尾迹模拟通常采用两种方式生成来流湍流: 人工湍流和前置模拟.可以考虑垂直方向均匀剪切的Mann 方法[48-49]在风力机模拟中广泛应用.在此类方法中,来流的速度脉动可以通过如下方式确定

其中,k为波数,nj(k) 为高斯随机复变量,其方差为1.上述公式的参数Cij通过如下方法确定

人工湍流方法只能在一定程度上反映大气边界层湍流特性(即满足相近的能量谱),前置模拟方法[50]生成的湍流更接近真实的大气边界层湍流,但需要和风力机尾迹湍流模拟相近的计算量.前置模拟方法通过以下步骤生成来流湍流.

(1)在前置模拟中,将湍流充分发展(统计定常情形)或发展到可以生成满足一定条件来流湍流的流动状态;

(2)在一定长度的时间段内,保存某一流向位置截面上的瞬时速度场;

(3)将得到的瞬时速度场进行适当的时间-空间插值,生成与风力机尾迹模拟的网格和时间步长一致的来流湍流.

采用前置模拟方法生成来流湍流多用于来流方向不变的统计定常情形,考虑风向改变的情形可以参照该工作[51].另一方面,对于真实的各向异性复杂地形,采用前置模拟方法生成与实际情形相近的来流湍流有较大困难.在这些工作[52-53]中,作者在入口处将复杂地形逐渐过渡到平坦地形,再将通过平坦地形前置模拟生成的来流湍流施加在入口处,得到的风力机电量产出结果和实测结果较好吻合.这些算例中的地形较为平缓,海拔高度变化不大.对于地形更为复杂的算例,可通过扩大计算区域的方式尽可能降低入口条件的影响.

2 尾迹机理

风力机尾迹可分为近尾迹和远尾迹.来流接近风力机时,压力升高,速度降低,风力机将风能转化为电能,压力发生台阶式下降,在其后方形成低压区,随着离风力机距离的增加,压力逐渐恢复,速度亏损继续增加,在某一位置达到最大后,速度开始恢复.由于尾迹和环境风速不同,在风力机尾迹边缘会形成剪切层,随着离风力机距离的增加,剪切厚度增加并在尾迹中心相遇,这一位置通常认为是近尾迹结束、远尾迹开始的位置.由叶片旋转导致尾缘涡卷起、脱落形成的螺旋状叶尖涡和中心涡是近尾迹的主要特征.假设叶尖涡的环量为 Γ,对于有个叶片的风轮,由于环量守恒,在中心会形成环量为的中心涡.在风力机远尾迹,尾迹中心通常不会停留在风轮轴线位置,而是以较低频率(远低于风轮转动频率)左右、上下摆动,形成尾迹蜿蜒.

图2 显示了数值模拟得到的叶尖涡、中心涡和蜿蜒示意图.以下章节讲分别介绍尾迹的时均特性,叶尖涡、中心涡和尾迹蜿蜒的主要特征和机理.本文将局限于风向不发生改变的情形,对于风向发生改变时的尾迹特性,读者可参考文献 [54-55].

图2 风力机叶尖涡、中心涡和蜿蜒示意图.云图显示瞬时流向速度.叶尖涡通过λ-2 准则可视化.中心涡通过流向涡量显示Fig.2 Schematic for tip vortices,hub vortex and meandering of turbine wakes.The contours show the instantaneous streamwise velocity.Tip vortices are visualized using the λ-2 criterion.The hub vortex is shown using the streamwise vorticity

2.1 尾迹的时均特性

尾迹具有风速低、湍流强度高两大主要特性.近尾迹速度亏损受风力机设计影响,比如叶片径向的升阻力分布和机舱,靠近叶片根部,由于较低升阻力和较低实度,会形成局部相对速度较高的区域,同时受机舱尾迹影响,具有复杂的径向分布特征.在远尾迹,速度亏损受风力机设计影响小,可以采用高斯分布近似.风力机导致的尾迹湍流强度主要分布在远尾迹.在近尾迹,周期性叶尖涡会带来名义湍流强度,若采用相位平均可以避免这部分湍流强度.在风力机下游,高湍流强度首先出现在与翼尖平行的位置,随着与风力机距离增加,高湍流区域在径向扩张,最终分布在整个尾迹及周围区域.在地面垂直方向存在流动剪切时,高湍流强度常出现在尾迹上边界.

尾迹时均特性受风力机运转工况和来流湍流等因素影响[56-59].风力机运转工况决定了风力机近尾迹速度亏损,比如,根据一维动量理论,近尾迹速度亏损为 ΔU1=2aU∞,其中a为轴向诱导因子.受大气稳定度和地面粗糙度等因素影响的来流湍流强度影响远尾迹速度亏损恢复.Xie 和Archer[56]研究了不同风速、不同风力机工况下的速度亏损自相似特性和尾迹湍流特性.钱耀如等[57]采用致动线模型模拟了低来流湍流度下的风力机尾迹.胡丹梅等[58]研究了动态来流对风力机尾迹的影响.韩玉霞等[59]研究了来流湍流强度对尾迹速度恢复的影响.实验和数值模拟结果[56,58]显示,来流湍流强度越高,尾迹速度亏损恢复越快.

2.2 叶尖涡和中心涡

叶尖涡位于区分尾迹区和自由流区的剪切层,影响尾迹与自由流的相互作用,相邻叶尖涡相互诱导,并受剪切层影响.Lignarolo 等[60]研究了叶尖涡涡对的不稳定机制,分析了平均动能通量,发现了具有不同动量掺混特性的两个区域,即,掺混受到叶尖涡抑制的近尾迹区域,及叶尖涡失稳后的高效掺混区域.叶尖涡稳定性依赖于来流湍流.在均匀来流或低湍流度来流时,在远尾迹仍可观察到叶尖涡[61],但在高湍流强度来流时,叶尖涡局限于近尾迹区,在风力机下游2 到5 个风轮直径处叶尖涡特征已不明显[60,62-64].叶尖涡通常认为具有较规则的几何形状.Yang 等[65]通过2.5 MW EOLOS 风力机的大尺度PIV[66]和大涡模拟研究发现具有尾巴状结构和二次涡的复杂叶尖涡,并分析表明该复杂叶尖涡结构与离心不稳定性相关.

中心涡强度与叶尖涡相关,并受机舱尾迹影响.对风洞实验结果进行稳定性分析发现中心涡特征在远尾迹仍然存在[67-68].Kang 等[69]通过对水动力机的几何解析大涡模拟和致动盘/致动线模拟发现,中心涡在向下游移动的过程中,不断向外扩张,与剪切层作用,激发或增强远尾迹蜿蜒,而没有机舱的风力机参数化模型不能准确预测中心涡作用.1.2.1 节介绍的机舱和叶片的致动面模型在较粗网格时可较准确模化机舱作用,预测尾迹蜿蜒,计算得到的湍动能和蜿蜒特征频率与几何解析大涡模拟和实验测量结果很好吻合[9].

2.3 尾迹蜿蜒

蜿蜒是远尾迹的低频大尺度运动.一方面,蜿蜒增加尾迹与周围高速流掺混加速尾迹恢复,从而提高下游风力机电量产出,另一方面,蜿蜒导致下游风力机在自由来流和上游风力机尾迹间切换,从而增加下游风力机所受疲劳载荷.蜿蜒受多种因素影响,如风向改变、来流湍流强度、风力机尾迹等.对风向不变的情形,尾迹蜿蜒有两种产生机制: 来流大尺度涡机制和剪切层不稳定机制.大尺度涡机制认为,大气边界层湍流中的大尺度涡结构对风力机尾迹的对流输运导致了尾迹蜿蜒.该机制得到了实地观测和数值模拟的验证[26,70-72],并被用于发展动态尾迹蜿蜒模型[26].剪切层不稳定机制认为,与钝体绕流相似,尾迹剪切层失稳导致了尾迹蜿蜒.该机制得到了不稳定性分析、实验和模拟结果的证实[73-78].观测和模拟结果[79-80]显示两种机制同时存在于风力机尾迹.

尾迹蜿蜒特性研究主要关注其特征频率、幅值和波长等特征.研究显示,由于剪切层失稳导致的蜿蜒特征频率所对应的斯特劳哈尔数(S t=f D/U)在0.1～ 0.5 的范围[73-79]内,与钝体绕流S t数相近,而大尺度涡导致的尾迹蜿蜒频率远低于剪切层失稳导致的尾迹蜿蜒频率[79-80].尾迹蜿蜒幅值决定了风力机会受尾迹影响的横向范围.Foti 等[74]采用浸没边界方法直接解析叶片几何,模拟了直径为0.128 m 的模型风力机尾迹.模拟结果显示在远尾迹(＞6D),蜿蜒幅值为风轮直径10%～ 30%,最大的尾迹蜿蜒幅值约为风轮直径的50%,蜿蜒波长为风轮直径的1～ 2倍,最大的尾迹蜿蜒波长约为风轮直径的3.5 倍.进一步,Foti 等[81]研究了不同尺寸、不同设计风力机尾迹蜿蜒,考虑的风机包括风轮直径0.128 m 的模型风力机[74]、1.1 m 的G1 模型风力机[72]、风轮直径27 m 的SWiFT 风力机[73]、风轮直径96 m 的EOLOS风力机[65-66].研究结果显示,在以风轮直径作为无量纲参数时,不同尺寸风力机的尾迹蜿蜒幅值有一定差异,尤其是SWiFT 风力机(远尾迹平均蜿蜒幅值约为0.15D)与0.128 m 模型和EOLOS风力机(远尾迹平均蜿蜒幅值约为0.1D)的对比,而在以推力和尾迹蜿蜒频率定义的特征长度其中T为推力,fm为尾迹蜿蜒频率)进行无量纲化时,不同设计风力机的尾迹蜿蜒幅值的差异显著降低.不同设计风力机的尾迹蜿蜒波长均在1D到2D之间.在文献[81]中,Foti 等研究了G1 风力机在不同工况下的尾迹蜿蜒,发现不同运转工况下的风力机尾迹蜿蜒的幅值和波长具有相似性.Kang 等[69]的研究显示机舱对于远尾迹蜿蜒的发生、发展有重要作用.通过数值实验,Foti 等[75]指出机舱对于大规模风电场中的尾迹蜿蜒也有重要作用,结果显示有机舱数值模拟得到的尾迹蜿蜒幅值和波长在远尾迹高于没有机舱的模拟结果.Yang 等[82-83]采用大涡模拟和致动面方法研究了不同偏航角的尾迹蜿蜒特性,基于偏航风力机尾迹分解,推导了不同尾迹特征的特征速度和特征长度与偏航角的关系式,发现在以这些特征速度和特征长度进行无量纲化时不同偏航角的风力机的尾迹速度亏损、湍动能、尾迹蜿蜒存在相似性.

3 结语与展望

本文综述了风力机尾迹的模拟方法,主要流动结构和流动机理.采用风力机参数化模型的大涡模拟是目前模拟风力机尾迹湍流的主要方法.对于单个风力机尾迹,致动面模型是合适的选择.对于大规模风电场,如果主要关注点在于风力机远尾迹,对时间/空间解析度要求低的致动盘模型是更为经济的选择.尾迹的主要流动结构包括叶尖涡、中心涡和蜿蜒.蜿蜒主要发生在远尾迹,对风电场整体性能有重要影响.本文介绍了尾迹蜿蜒机理和主要特征.对于风力机尾迹的平均特性,人们已有较好的理解.在选择合适的参数后,工程模型可以较为准确地预测尾迹平均特性.然而,对于风力机尾迹湍流,人们还缺乏深入理解,更没有具备预测能力的低阶模型.

随着风电在整个能源系统中所占比重的不断增加,对于风电场电量产出预测与风力机布置和控制有了更为精细的时间和空间解析度要求[3],以便进一步降低风电度电成本.在这个时间和空间尺度上,大气边界层和风力机尾迹湍流起着至关重要的作用.深入研究相关机理,发展可以预测大气边界层和风力机尾迹湍流特性的快速模型变得尤为迫切.这主要面临两方面的难题: 一方面,对于不同工况的风力机尾迹机理仍缺乏深入理解,比如,复杂地形的风力机尾迹[52]、浮式风力机尾迹[84]、风电场尾迹等;另一方面,如何将不同理想工况的结果推广到实际工况,并发展工程可用的快速低阶模型.