风力机之间的尾流效应会造成下游机组的发电量减少，从而导致海上风电场的发电量损失10%~20%[1]．交错排布能有效提升风电场的发电量，然而，现有研究主要集中在单一类型的水平轴风电场．Xie等[2]在一个大型水平轴风机(HAWT)周围安装了20个小型垂直轴风机(VAWT)，发现其发电量比传统水平轴风力发电场高出32%，且大型HAWT单独提取的功率增加了10%，原因在于小型VAWT的存在加快了尾流区域的速度恢复．Shao等[3]采用致动线模型(ALM)方法研究了HAWT和VAWT串列排布对功率的影响，结果表明：VAWT位于HAWT上游，可以增强HAWT的发电量，是由于风速在VAWT上加速所致；VAWT位于HAWT下游，会略微降低HAWT的发电效率，但总功率仍然增加．Pichandi等[4]发现：在山地地形下将VAWT放置在HAWT上游，不仅有助于提高给定空间的整体功率输出，而且还有助于提高HAWT的功率输出；功率增加与两台风力机之间的间距、叶片尖端之间的间隙及斜坡倾斜度有关．近年来，风电场中机组的交错排布引起了广泛关注，然而研究主要集中在HAWT周围安装小型VAWT及两者串列排布上，两者的尺寸差异较大，且未考虑水平轴风电场的实际排布结构．本研究以典型的大型水平轴风电场布局为基础，在其中错列排布大尺寸VAWT，使用基于OpenFOAM开发的turbinesFoam[5]对风力机进行ALM设置，探讨VAWT排布位置对整体功率及流场特性的影响规律．1 数值方法1.1　致动线基本理论ALM方法[6]是在N-S方程中添加体积力源项，以模拟叶片对流场的作用．这种方法无须建立叶片固壁边界求解叶片边界层流动，网格数量少，降低了计算时间．ALM中流体控制方程为不可压缩黏性N-S方程和连续性方程∂V/∂t+V∇V=-(1/ρ)∇p+v∇2V+fε；∇V=0，式中：V为风速；t为时间；ρ为空气密度；p为流体压力；ν为流体运动黏度；fε为叶片体积力．在ALM中，叶片被简化为由一系列体积力作用点组成的致动线，这些体积力计算由当地攻角与升阻力系数确定，叶素受到的升力L与阻力R分别为L=ClρUrel2cr/2；R=CdρUrel2cr/2，式中：r为叶素宽度；c为当地弦长；Urel为相对速度；Cl和Cd为升阻力系数，是雷诺数和攻角的函数．在实际应用中，须根据具体的雷诺数和攻角值来选择合适的翼型升阻力系数，使用的升阻力系数来自NREL 5 MW风力机[7]官方参数．图1为叶素上的气动力示意图，图中Fx和Fθ分别为合力在法向和切向上的分量．由图1可知10.13245/j.hust.250017.F001图1叶素上的气动力α=φ-β；Urel=Ux2+(ωr-Uθ)2，式中：α为当地攻角；φ为入流角；β为安放角；Ux和Uθ分别为轴向速度和切向速度；ω为旋转角速度；r为旋转半径．因此，单位长度展向体积力为f=(L,R)=ρUrel22c(CleL⃗+CdeR⃗)/2，式中eL⃗和eR⃗分别为L和R的单位方向向量．在求出每个叶素的体积力后，须考虑该力作用于流场的影响范围，且避免奇异性，对体积力进行光顺处理，采用三维高斯函数ηε进行体积力光顺，ηε=exp[-(di/ε)2]/(ε3π3/2)式中：di为网格中心到叶素点的距离；ε为体积力分布因子．因此，光顺后的体积力为fε=f⊗ηε=∑i=1Nfexp[-(di/ε)2]/(ε3π3/2)，式中N为致动点数目．1.2　湍流模型及求解器采用OpenFOAM的瞬态不可压缩流求解器pimpleFoam，该求解器为PISO算法和SIMPLE算法的结合．湍流模型采用LES，其中大尺度涡采用网格分解，小尺度涡采用亚格子模型(sub-grid scale，SGS)．滤波后的LES控制方程形式如下∂u˜i∂xi=0；∂u˜i∂t+u˜j∂u˜i∂xj=-1ρ∂p˜∂xi+v∂∂xj(∂u˜i∂xj)-∂τijSGS∂xj+fε，式中上标~表示分解量，τijSGS=uiuj̃-u˜iu˜j为SGS应力，采用Smagorinsky亚格子模型描述亚格子应力．SGS应力的偏差部分τijSGS,dev.=τijSGS-τkkSGSδij/3使用Boussinesq方法估计，τijSGS,dev.=-2vSGSSij̃，其中vSGS为亚格子涡黏系数，为可解尺度的涡流变形率，可通过下式计算vSGS=(CSΔ)2Sij̃，式中：Sij̃=12(∂u˜i∂xj+∂u˜j∂xi)，Δ=(ΔxΔyΔz)1/3为滤波尺寸；CS为Smargorinsky常数，取值0.158．2 风力机网格及模型验证分析HAWT采用NREL 5 MM风力机[7]为计算模型，参数设置：转子直径126 m；轮毂高度90 m；额定风速11.4 m/s；额定转速12.1 r/min．在ALM设置方面参考了Yu等[8]的研究，并模拟塔架及轮毂．VAWT采用额定功率为200 kW[9]的H型转子VAWT作为计算模型，参数设置：直径26 m；轮毂高度40 m；叶片长度24 m；叶片弦长0.4~0.9 m．为了简化计算，对VAWT叶片弦长做简化处理，取弦长为0.9 m，并忽略塔架及支撑杆．首先对单一HAWT进行网格无关性验证，计算域大小为3 150 m×1 260 m×630 m，而VAWT尺寸较小，对VAWT网格无关性验证时，计算域减小为2 000 m×400 m×200 m．入口边界为速度入口，出口边界为零梯度条件，顶部和两侧为滑移壁面，底部为壁面边界条件．为模拟ABL条件，速度入口采用atmBoundaryLayerInletVelocity类型，其速度剖面定义为v(z)=(v*/κ)ln(z/z0)，式中：ν*为流动摩擦速度；κ为卡曼常数；z0为地面粗糙长度，取z0=0.01 m．为同时满足HAWT及VAWT时间步长计算要求，在满足库朗数情况下选择时间步长为0.02 s．初始网格尺寸均为10 m，并采用多级加密处理方式对风轮区域网格进行加密．HAWT与VAWT风轮附近网格尺寸及不同尺寸网格下的功率系数(CP)如表1所示．根据Jonkman等[7] NREL 5MW风力机为额定风速时，CP=0.4418，对比发现中网格尺寸时CP误差最小，根据Möllerström等[10]的研究，VAWT在额定风速时，功率系数为CP=0.34，对比不同网格尺寸发现中网格尺寸时CP符合度最好．10.13245/j.hust.250017.T001表1HAWT与VAWT不同网格尺寸下的功率系数网格类型网格尺寸/m网格数量/106CPHAWT粗网格2.500 005.419 1060.439中网格1.250 006.216 0560.441细网格0.625 007.319 9280.485VAWT粗网格0.625 003.183 9720.328中网格0.312 504.428 5720.344细网格0.156 255.422 8520.373图2为三种网格尺寸下HAWT的涡量等值面，由图2可知：三种网格都捕捉到了HAWT的大尺度叶尖涡和叶根涡，图2(b)和(c)还捕捉到了塔涡．随着网格尺寸减小，HAWT叶片形状越接近一条线且对涡捕捉增多，这表明网格尺寸越小，捕获的流场信息越精准．VAWT的流场结构也证明了这一结论，如图3所示．综合功率系数、流场结构及所需的计算时间，确定HAWT与VAWT风轮附近最小网格选择为1.25 和0.312 5 m．10.13245/j.hust.250017.F002图2HAWT涡量等值面(色标单位：m/s)10.13245/j.hust.250017.F003图3VAWT涡量等值面(色标单位：m/s)3 HAWT与VAWT相互作用3.1　工况及计算域设置以大型水平轴风电场布局为背景，机组流向间距为7D (D为HAWT风轮直径)，横向间距为3D．在HAWT机组间隙排布VAWT，以研究在风电场中HAWT满发状态下，不对HAWT做变桨调节时，两种机组的相互作用，分析机组功率输出特性及尾流特性．为减小风电场尾流叠加带来的影响，突出HAWT与VAWT之间尾流的相互作用，且减少计算时间，选取一组HAWT与VAWT开展研究，计算域如图4所示．以HAWT轮毂中心为坐标原点(图中红色圆点位置)，入流方向为x轴正方向．根据风电场机组流向及横向机组间隙，使得VAWT能均匀排布在间隙中，给出六种不同位置．VAWT位置坐标如表2所示，并设置仅包含HAWT或仅包含VAWT的对照组．10.13245/j.hust.250017.F004图4计算域示意图10.13245/j.hust.250017.T002表2VAWT位置坐标工况dx/Ddy/D1-3.50.02-3.51.53-1.01.541.01.551.01.06-3.51.03.2　风力机功率输出特性分析HAWT，VAWT，HAWT+VAWT功率输出如表3所示．相较于对照风力机功率，工况2，3，6 中HAWT输出功率有所提升，其中工况6中HAWT功率提升了0.44%．然而，工况4，5中HAWT输出功率有所下降．这一现象表明：当VAWT位于HAWT前方时，由于VAWT的存在，风速经过VAWT时加速了后方风速，不进行变桨调节时，满发状态下的HAWT依然增加了功率．相反，VAWT位于后方则会降低HAWT功率输出．在工况1中，尽管VAWT位于HAWT前方，但HAWT输出功率显著下降，主要是由于HAWT的风轮底端低于VAWT的风轮底端，位于上游的VAWT产生的尾流作用于HAWT的风轮，降低了HAWT性能，导致其输出功率显著下降．10.13245/j.hust.250017.T003表3HAWT，VAWT，HAWT+VAWT功率输出工况HAWT功率VAWT功率HAWT+VAWT功率14 8711044 97625 0101055 11535 0011025 10344 9941155 10854 9911205 11165 0221055 127对照组(单一HAWT)5 000对照组(单一VAWT)106kW不同工况下VAWT功率输出存在一定的差异，且由于大气边界层(ABL)入流条件，输出功率低于额定功率．在工况1，2，3，6中VAWT位于HAWT前方，VAWT输出功率略低于对照风力机功率；而在工况4和5中VAWT位于HAWT后方，输出功率相较于对照风力机功率均有所提升，其中工况5输出功率提升了13.21%．此外，通过对比工况2与3发现：当VAWT位于HAWT前方且横向距离相同时，随着流向距离逐渐接近HAWT，VAWT的输出功率呈现降低趋势．而对比工况4与5发现：当VAWT于后方且流向距离相同时，随着横向距离与HAWT的接近，VAWT的输出功率呈现上升趋势．在水平轴风电场机组间隙中加入VAWT后，除工况1以外，总输出功率相较于单一HAWT都有提升，其中，工况5总输出功率提升了2.22%，工况6总输出功率提升了2.54%．因此，在水平轴风电场机组间隙合适位置加入VAWT，能够有效提升风电场功率输出，并实现对水平轴风电场机组间隙的充分利用．3.3　风力机尾流特性分析VAWT轮毂高度水平截面瞬时速度云图如图5所示．可以看到：VAWT尾流在经过HAWT近尾流区域时有明显的偏转，表明两者之间存在排斥作用．然而，在远尾流区域两者相互混合，这种混合作用促进了HAWT尾流蜿蜒．此外，HAWT风轮前方有低风速区域，在近尾流区域两侧存在风速加速区域，这也解释了VAWT位于HAWT前方时，功率会下降，位于后方时会有明显的提升．在工况1中，VAWT在HAWT前方串列排布，VAWT尾流完全作用于HAWT风轮，使得HAWT风轮前方气流发生紊乱，增强了湍流，破坏了HAWT的气动性能，导致HAWT风轮无法有效地捕获风能，因而HAWT尾流速度高于其他工况．工况2与3中VAWT位于HAWT前方，当VAWT尾流到达HAWT附近时，尾流发生了偏转．工况4与5位于HAWT后方，两种工况中VAWT处于HAWT影响范围内，因而VAWT尾流未发生明显的偏转，但在工况5中，VAWT尾流与HAWT尾流之间的相互作用更为强烈，VAWT尾流使得HAWT尾流蜿蜒发展更加明显，尾流速度恢复也明显加快．工况6中，VAWT尾流在经过HAWT时表现出更强的偏转效应．在两者尾流接触区域，出现了速度增大的现象，这一现象将有助于尾流的恢复．10.13245/j.hust.250017.F005图5瞬时速度云图(色标单位：m/s)HAWT尾流在6D位置的横向截面平均速度等值线如图6所示．加入VAWT后，HAWT尾流受到了一定程度的影响．工况1和6的尾流速度恢复高于对照组，而其他工况的尾流速度恢复则略低于对照组．相较于工况1，工况6不仅尾流速度恢复更快，而且HAWT功率也有所提升．这些结果表明，VAWT的存在可以加速HAWT尾流的恢复过程．10.13245/j.hust.250017.F006图6平均速度等值线图(色标单位：m/s)VAWT轮毂高度水平截面涡量云图如图7所示．在工况1中，VAWT涡流对HAWT涡流产生显著影响，当VAWT涡流到达HAWT风轮时，HAWT涡流结构被破环，叶尖涡破碎，对HAWT风轮性能产生严重影响．工况2与3中VAWT涡流都发生偏转，但涡流结构保持稳定．工况4与5中VAWT位于HAWT后方，VAWT涡流有偏向于HAWT趋势，说明HAWT对VAWT涡流存在吸引作用．工况5相较于工况4，VAWT横向距离HAWT更近，两者涡流吸引效应强烈，但依然保持各自涡结构．工况6相较于工况2横向距离HAWT更近，在VAWT涡流经过HAWT近尾流区域时表现出更强的排斥现象，在远尾流处相互作用，融于HAWT涡流之中．10.13245/j.hust.250017.F007图7涡量云图(色标单位：s-1)归一化处理HAWT轮毂中心水平和垂直截面的风速廓线如图8所示，U为HAWT轮毂高度水平速度，Zhub和Uhub分别为HAWT轮毂高度及轮毂高度风速．由图8(a)可以看出：HAWT近尾流区域平均速度廓线都出现了横向不对称性．这是由于叶片旋转导致塔架后出现斜流所致[11]．同时由速度廓线可以看出：VAWT的加入对HAWT轮毂高度处风速影响较小，在尾流6D位置处，工况3~5的速度廓线与对照风速廓线存在明显的差异．从相对位置来看，这三种工况的流向距离HAWT较近，从而导致了HAWT尾流速度较低．在远尾流区域，速度廓线与对照速度廓线相近，这表明两者尾流的相互作用，加速了尾流的恢复．图8(b)显示了HAWT轮毂中心垂直截面的风速廓线，主要差异也表现在6D位置，工况3~5降低了HAWT尾流速度，但在远尾流区域，相较于对照速度廓线，在轮毂高度以上有更高的风速．在工况6中，垂直截面的速度廓线与对照速度廓线在近尾流区域相近，在远尾流区高于对照速度廓线，这表明VAWT尾流对HAWT尾流恢复有促进作用．10.13245/j.hust.250017.F008图8HAWT风速廓线归一化处理VAWT轮毂中心水平截面风速廓线如图9所示，u为VAWT轮毂高度水平速度，d和uhub分别为VAWT风轮直径和轮毂高度风速．从图中可以看出：工况1中VAWT尾流速度廓线低于对照速度廓线，主要是HAWT风轮降低了其前方风速，进而影响了VAWT的尾流恢复．工况2中VAWT近尾流区域受HAWT影响较小，然而在远尾流区域，由于偏转效应的影响，速度廓线在靠近HAWT一侧高于对照风速廓线，而在另一侧低于对照风速廓线．工况3中VAWT尾流在2d位置时与对照风速廓线相近，在6d左右受到HAWT尾流影响，速度廓线出现明显的偏转，且尾流速度高于对照风速廓线，说明两者尾流相互作用可以加快VAWT尾流恢复．工况4和5位于HAWT后方，尾流速度廓线规律与其他工况不同，在10d左右位置发生变化，速度廓线在靠近HAWT一侧低于对照风速廓线，而在另一侧高于对照速度廓线．工况6中VAWT尾流在到达HAWT之前一直与对照风速廓线相近，在经过HAWT后，出现与工况2相似的规律．10.13245/j.hust.250017.F009图9VAWT风速廓线HAWT轮毂中心垂直截面不同位置归一化湍动能(TKE)廓线分布如图10所示．TKE被定义为k=0.5(〈u'u' 〉+〈v'v' 〉+〈w'w' 〉)，其中：u'，v'，w'分别为三个方向上的脉动速度分量；〈 〉为时间平均．在加入VAWT后，近尾流区域内TKE略微增大，这可能会增加风轮疲劳载荷，但更强的环境湍流将促进HAWT尾流的混合，从而使尾流恢复速度加快[2]．在工况1的近尾流区域，TKE显著增强，不仅轮毂高度以下的TKE得到增强，HAWT轮毂高度以上TKE也有一定的增强．值得注意的是，在远尾流区域，各工况下TKE在轮毂高度区域均有所降低，这表明在远尾流区域，VAWT与HAWT的相互作用加速了尾流的恢复过程．10.13245/j.hust.250017.F010图10HAWT湍动能廓线4 结论a．当VAWT的风轮底端高于HAWT的风轮底端位置且在HAWT前方串列排布时，VAWT的尾流会对HAWT产生不利影响，导致其功率降低．而在HAWT前方错列排布时，会对HAWT产生有利影响，提升其功率输出．b．当VAWT位于HAWT后方时，会略微降低HAWT功率，但由于HAWT影响，VAWT功率会得到显著提升，相比于单一HAWT机组，整体功率仍会增加，且最大增加2.22%．c．在水平轴风电场机组间隙加入VAWT，VAWT尾流与HAWT尾流相互作用，呈现出更快的尾流恢复速率．d．在水平轴风电场机组间隙加入VAWT，不仅可以提升风电场的发电量，而且能充分利用风电场空间资源，可以在不改变原有布局的情况下对正在运行的风电场进行改造，具有广阔的应用前景．