在大气边界层中，风剪切使风电机组输出功率降低，机组承受载荷加剧．风剪切来流条件下，风力机尾流效应更强，对大型风电场的整体性能与安全运行有重大影响．风剪切来流使风力机下游区域尾流分布呈现不对称性，且波动性增强，在风力机设计阶段必须考虑到风剪切的影响[1]．下击暴流是一种特殊且常见的风剪切，下击暴流在垂直方向存在先增大后减小的速度梯度分布规律，对风电机组的安全运行、发电质量及尾流特性有重大影响．在雷暴天气条件下，下击暴流发生的概率高达60%～70%，下击暴流的持续时间一般在十几分钟甚至持续几个小时．文献[2-3]归纳总结了下击暴流与大气边界层风场的差异，包括径向和竖向的平均风速剖面、非平稳的风速和风向的复杂多变．下击暴流发生的频率高且持续时间长，每年都有许多的建筑物及输电塔等受到下击暴流风灾的影响[4]，风力机作为大型发电设备受到下击暴流的影响很大，因此研究下击暴流垂直风剪切来流对风力机尾流特性及转矩特性的影响十分重要．本研究基于OpenFOAM开源软件，采用大涡模拟结合致动线模型(LES-ALM)[5-6]的数值模拟方法，对均匀来流、B类风场风剪切来流及下击暴流来流垂直风剪切三种来流条件下的风力机尾流特性及转矩特性进行对比分析，讨论下击暴流来流条件下风力机尾流轴向速度特性、湍流强度特性及低速轴转矩波动特性．1 数值模拟方法1.1　大涡模拟控制方程大涡模拟[7-8]是通过滤波方法将大尺度涡通过纳维-斯托克斯方程直接求解，小尺度涡由亚格子模型建模模拟大涡，控制方程[9]为∂u¯i∂t+∂(u¯ju¯i)∂xj=-1ρ∂p¯∂xi+v∂2u¯i∂xj∂xj-∂τij∂xj+fi；(1)∂u¯i/∂xi=0, (2)式中：u¯i和u¯j为速度场的笛卡尔分量；t为时间；xi和xj为位移；ρ为当地空气密度；p¯为平均压力；τij为SGS应力；fi为对风轮效应进行等效时的体积力．为满足本研究的求解需要，选择耗散型SGS应力模型[9]，表达式为τij-δijτkk/3=-2(CsΔ)2|S¯|S¯ij=-2νtS¯ij，(3)式中：δij为狄拉克张量，当i=j时，δij=1，当i≠j时，δij=0；τkk为主对角线上雷诺应力；Cs为斯马戈林斯基系数；Δ为过滤尺度；|S¯|为应变率张量的模；S¯ij为应变率张量；νt为SGS涡黏性系数.1.2　致动线模型致动线模型(ALM)[10]是将实体风轮叶片简化成致动线，在致动线上布置若干计算点，每个计算点处叶素产生的气动力为fe=ρvref2c(Clel+Cded)/2，(4)式中：vref2为当地相对风速的平方；c为弦长；Cl和Cd分别为翼型升、阻力系数；el和ed分别为升力和阻力方向上的单位向量.高斯光顺函数表达式为ηε(d)=[1/(ε3π3/2)]exp[-(d/ε)2]; (5)体积力可以表示为fi=(fe⊗ηε(d))ei，(6)式中：fi为N-S方程中的源项；d为fi与fe之间的距离；ε为光顺函数的高斯分布因子，一般ε越大，体积力分布越平滑，体积力影响范围越大．1.3　下击暴流数学模型目前在下击暴流水平风速模型的研究中，OBV (Oseguera & Bowles/Vicroy)模型[11]是被普遍认可的，风速模型数学表达式为V(z)=1.22[exp(-0.15Z/Zmax)-exp(-3.217  5Z/Zmax)]-Vmax, (7)式中：V为来流风速；Vmax为下击暴流的最大水平风速；Z为风速对应的高度；Zmax为最大风速对应的高度．下击暴流口径在1 000～3 000 m之间[14]，风轮直径为126 m，风轮在风场中的比例在0.04～0.13之间．因此，本研究对下击暴流垂直风剪切来流风况做出以下假设：a. 下击暴流风暴中心在风轮之前，且风暴中心点在XOZ平面上，为静止型下击暴流(主流方向为X轴，Y轴为横向方向，Z方向为垂直地面方向，O为原点)；b. 下击暴流冲击地面后的径向风速在风轮前500 m处于完全发展阶段且风速保持不变；c. 下击暴流Y轴方向速度为0 m/s．参照文献[12]的比例关系，选取最大风速为11.4 m/s，其对应高度为90 m．1.4　剪切来流数学模型剪切来流根据文献[13]给定的风剖面模型U(Z)=1.77U(Z/Z(G))α, (8)式中：U=8.010 2 m/s(通过轮毂中心高度处额定风速U=11.4 m/s反推算出)；Z(G)=350 m．本研究对轮毂高度处风速都设置为11.4 m/s．1.5　边界条件计算时入口边界采用速度入口，出口边界为压力出口，壁面采用滑移边界条件．图1为流场入口处下击暴流垂直风剪切的计算值与目标值(目标函数式(5))的轴向速度对比图，由图1可知LES-ALM方法所模拟的结果与目标值对应良好．10.13245/j.hust.210613.F001图1下击暴流垂直风剪切廓线1.6　转矩验证LES-ALM方法为成熟的数值模拟方法，在文献[14]中对美国国家可再生能源实验室(NREL)实验机组Phase Ⅵ两叶片风力机功率进行了验证，通过实验数据与LES-ALM方法的计算结果进行对比分析，得出LES-ALM方法对风轮气动特性的计算精度较高，满足计算要求．2 风力机模型及网格参数确定2.1　风力机模型本研究以美国国家可再生能源实验室(NREL) 5 MW风力机作为研究对象[15]，其中叶片数为3个，额定功率为5.29 MW，风轮直径为126 m，轮毂高度为87.6 m，切入风速为3 m/s，切出风速为25 m/s，额定风速为11.4 m/s，仰角为5°．2.2　计算域和网格本研究采用长方形计算域，长3 000 m，宽600 m，高600 m，其中x轴为主流方向，y轴为横向方向，z方向为垂直地面方向，风轮位于计算域入口500 m处，风轮中心坐标为(0 m，0 m，90 m)．为了清晰地捕捉叶尖、叶根的涡特性，在风轮直径上网格节点数设置50个，因此本研究给定的网格节点数为50个[16]；另外，对风轮附近区域网格进行两次加密，加密后网格最小尺度为2.5 m，网格总数为772.5万．2.3　高斯分布因子验证高斯分布因子(ε)影响体积力分布的范围和平滑程度，ε与风轮直径上的网格尺度η (η=D/N)有关，二者关系为ε=nη，其中n为常数．本研究以给定的网格尺度η=2.5 m，变常数n来改变高斯分布因子．取n=1.5，2.0，2.5 (ε=1.5η，2.0η，2.5η)三个值进行验证，以风轮输出功率相对误差为标准确定ε的最终取值，验证得出当ε=2.0η时，其相对误差为0.96%，误差最小，取ε=2.0η．具体计算方法见文献[15]．3 计算结果与分析3.1　尾流速度特性图2分别为均匀来流、B类风场风剪切和下击暴流来流条件下瞬时速度云图，其中：Ux为x方向(轴向)的风速；X为下游区域以风轮平面为基准的轴向距离；D为风轮直径.由云图可知：均匀来流条件下，在风轮后1D～7D之间尾流边界比较稳定，尾流基本呈对称分布，在7D以后尾流边界逐渐紊乱；B类风场风剪切条件下，尾流上边界在8D位置逐渐变得比较紊乱，尾流下边界在6D位置开始变得比较紊乱，尾流分布的非对称特性明显，尾流下边界比上边界更早趋向紊乱；而下击暴流垂直风剪切条件下，尾流上边界在7D位置趋于紊乱，而下边界在9D位置尾流边界逐渐紊乱，叶尖涡耗散加剧．与均匀来流相比，B类风场风剪切和下击暴流垂直风剪切来流条件下，风力机尾流分布的非对称特性比较明显；但下击暴流垂直风剪切来流条件下，风力机尾流上、下边界的变化规律与B类风场风剪切工况下正好相反，这是因为下击暴流垂直风剪切复杂的速度梯度造成的，来流在风轮上下叶尖位置的速度差不大，但经过风轮后，尾流下边界风速由于“狭管效应”被加速，叶尖涡明显耗散的位置距风轮较远．10.13245/j.hust.210613.F002图2不同来流条件下流场瞬时速度云图(上、中、下依次为均匀来流、剪切来流和下击暴流)图3为风轮上下游不同位置处的轴向平均速度随高度变化的剖面图，其中：Z为垂直高度；Z0为风轮中心高度(Z0=90 m)；V为来流风速；V0为Z0高度处风速(V0=11.4 m/s)．由于忽略了轮毂的影响，轮毂中心处风速受挤压作用使风速偏大，出现更加明显的“双峰”分布，随着尾流向下游发展，尾流和周围流场相互掺混，尾流轴向平均速度剖面逐渐被“单峰”分布代替；在轮毂中心高度，随着尾流向下游发展，风力机尾流轴向平均速度均为先减小后增大；下击暴流垂直风剪切条件下，在轮毂中心高度位置，尾流轴向平均速度在11D位置亏损最大，其亏损率(亏损率为尾流轴向平均速度和来流风速的差值与来流风速的比值)为40.2%；均匀来流和B类风场风剪切条件下，在11D位置轴向平均速度的亏损率为35.4%和38.4%，在17D位置，下击暴流垂直风剪切、均匀来流和B类风场风剪切条件下其亏损率分别为27.2%，19.3%和22.0%，与11D相比，17D位置尾流轴向平均速度分别恢复了来流的8.2%，16.1%和16.4%，说明在下击暴流垂直风剪切条件下，风力机尾流轴向速度恢复最慢．10.13245/j.hust.210613.F003图3风轮下游不同位置处轴向平均速度剖面下击暴流垂直风剪切存在先增大后减小的速度梯度，尾流区周围的流动速度较低，风力机尾流能量补充较小，尾流亏损更大，因此尾流轴向速度恢复相对更慢．3.2　尾流湍流强度特性图4为均匀来流、B类风场风剪切和下击暴流垂直风剪切条件下的瞬时涡量云图，其中Ω为涡量值．由图4可知：在均匀来流条件下，在1D～7D之间，尾涡强度和尺度基本不变，尾流涡耗散较弱，尾流边界比较稳定，在7D附近尾涡逐渐破碎、耗散为强度较小的涡；B类风场风剪切条件下，尾流下边界较上边界相比维持稳定的距离较短，尾10.13245/j.hust.210613.F004图4不同来流条件下瞬时涡量云图(上、中、下依次为均匀来流、剪切来流和下击暴流)流下边界在6D附近叶尖涡逐渐破碎，耗散加剧，而上边界则在8D附近叶尖涡逐渐破碎，耗散加剧；下击暴流垂直风剪切条件下，风力机尾流上、下边界叶尖涡出现破碎的位置与B类风场风剪切条件下相反，尾流上边界在7D附近叶尖涡逐渐破碎，而下边界在9D附近才出现破碎现象，与图2分析的结论一致，也说明了下击暴流和剪切来流条件下风力机尾流的非对称性比较明显；与均匀来流和B类风场风剪切相比，下击暴流垂直风剪切条件下叶根涡更容易发生破碎融合现象．图5为均匀来流、B类风场风剪切和下击暴流垂直风剪切条件下风轮轴向不同位置处的湍流强度平均值随高度变化曲线，其中I为湍流强度．由图5可知：随着尾流向后发展，湍流强度分布由“双峰”分布逐渐变成“单峰”分布；在轮毂中心高度处，三种来流工况下风力机尾流的湍流强度均为先增大后减小的现象；下击暴流垂直风剪切条件下，风力机尾流湍流强度在11D位置达到最大值，为11.94%；均匀来流和B类风场风剪切来流条件下，尾流湍流强度在11D位置分别为11.24%和10.86%；在17D位置，下击暴流垂直风剪切、均匀来流和B类风场风剪切来流条件下，尾流湍流强度分别为8.34%，8.65%和9.22%，与11D相比，17D处湍流强度分别衰减(减小)30.15%，23%和15%．在湍流强度衰减过程中，下击暴流垂直风剪切条件下，风力机尾流湍流强度衰减率(17D处尾流湍流强度和11D处湍流强度的差值与11D处湍流强度的比值)是均匀来流工况下的1.3倍，是B类风场风剪切工况下的2倍，说明在该过程中，下击暴流垂直风剪切条件下风力机尾流湍流强度衰减最快．10.13245/j.hust.210613.F005图5风轮下游不同位置湍流强度时均值3.3　风轮转矩特征表1为均匀来流、B类风场风剪切和下击暴流垂直风剪切条件下风轮低速轴转矩的时均值和标准差．下击暴流垂直风剪切来流条件下，低速轴转矩时均值是均匀来流工况下的90%，是B类风场风剪切工况下的94%，但标准差是均匀来流工况下的1.96倍，是B类风场风剪切工况下的1.86倍，说明在下击暴流垂直风剪切条件下，风力机低速轴转矩的均值减小，但标准差明显增大，说明转矩波动幅值增大，导致风力机输出功率不稳定，这一现象主要是非对称来流引起的风轮平面气动力不均匀造成的．10.13245/j.hust.210613.T001表1转矩时均值和标准差工况时均值标准差均匀来流B类风场风剪切下击暴流4 137.872 34 100.393 83 854.616 521.688 122.854 342.515 3kN•m图6为转矩功率谱密度图，图中：S为转矩功率谱密度；f为频率．风轮转速为12.06 r/min，风轮旋转频率fn=12.1/60  Hz=0.200  8  Hz，风轮为3叶片，其叶片通过频率fp=3fn=0.602  4  Hz．由图6可知：不同来流条件下，风轮转矩功率谱在叶片通过频率对应位置皆出现峰值，且功率谱变化特性比较相似，说明风轮转矩波动与风轮旋转周期有关，但在下击暴流来流条件下，低频段含能最高，主要出现在叶片通过频率附近，说明下击暴流垂直风剪切条件下风轮所受到的非对称载荷最强，转矩波动幅值最大，输电质量降低．10.13245/j.hust.210613.F006图6低速轴转矩功率谱图4 结论本研究基于OpenFOAM开源软件，采取LES-ALM计算方法，在均匀来流、B类风场剪切和下击暴流条件下对NREL 5 MW机组进行数值模拟研究，分析不同来流下风力机尾流速度特性、湍流特性和转矩特性，得出结论如下．a. 与均匀来流和B类风场风剪切相比，下击暴流垂直风剪切条件下，风力机尾流轴向速度亏损率最大，且恢复最慢，尾流效应更强；下击暴流垂直风剪切条件下，风力机尾流上、下边界的变化规律与B类风场剪切相反，尾流分布的非对称特性更加明显．b. 在轮毂中心高度处，风力机尾流湍流强度在衰减过程中，与均匀来流和B类风场风剪切相比，下击暴流垂直风剪切条件下，风力机尾流湍流强度衰减最快，尾流中心涡耗散加快．c. 与均匀来流和B类风场风剪切相比，下击暴流垂直风剪切来流条件下，风力机低速轴转矩平均值降低，但标准差增大，转矩波动幅度增加，因此下击暴流垂直风剪切将会导致风力机输出功率不稳定，降低其电能质量．