不同涡脱落模式下垂直轴风力机叶片的气动响应

郑小波，王红亮，徐文浩，郜志腾，3，冷 峻，李 晔，*

(1.兰州理工大学 能源与动力工程学院，兰州 730050；2.上海交通大学 船舶海洋与建筑工程学院，上海 200240；3.上海交通大学 多功能拖曳水池实验室，上海 200240)

0 引言

风能是可再生能源的重要组成部分，对节能减排具有重要意义，因此风能发电已经在全世界范围内快速发展[1]。相比陆上风电，海上风电的风能资源丰富，能量密度高，而且对环境影响较小，因此受到越来越多的关注[2]。我国海上风能资源开发潜力巨大，高效环保地推进我国海上风电的发展，是实现“碳达峰、碳中和”目标的重要举措[3]。

为了捕获更多风能资源，海上风力机叶轮直径已经突破200 m，大型化已经成为海上风力机发展的一个重要趋势[4]。虽然经过几十年的技术发展，三叶片逆风水平轴风力机（horizontal-axis wind turbine,HAWT）已经成为陆上多兆瓦风力发电的最普遍范式，但却不适用于海上风电的开发利用[5]。大型HAWT 重约400 吨的机舱布置在约100 m 的海拔高度，会使处于波浪、洋流和强风复杂环境中的浮式HAWT 产生极大的不稳定性，大大增加了安装和维护成本[6]。相比而言，浮式垂直轴风力机（vertical-axis wind turbine,VAWT，本文专指升力型Darrieus 风力机）避免了这一结构问题，将发电机、齿轮箱和制动器等传动系统布置在风力机塔筒根部，提升了风力机寿命时长。此外，VAWT 还具有结构简单、对风向不敏感、适应湍流和阵风环境等优点[7]，可以降低风力机的设计和建造成本，并提升工作时长。

辽阔的海面允许布置更大范围的风电场，但是为了提高风能密度（即单位面积的风能输出），如何布置多台风力机吸引了众多研究者。Sørensen[8]指出，HAWT 风电场要求风力机流向间距为5～10 倍的叶轮直径，因为在沿流向的紧凑分布中，各台HAWT 的功率系数（机械功率输出与通过叶轮横截面积的自由流功率之比）会随着间距减小而剧烈下降[9]。但Rajagopalan 等[10]也发现，在紧凑的VAWT 风电场中，各台风力机的功率系数降低有限，甚至会有所提升。Nes[11]提出了浮式风墙的概念，通过增大扫掠面积，在有限的海域提升了风能发电能力，在降低安装维护成本的同时获得高输出风能密度。但在这种风墙概念中，阵列单元仅采用了小型HAWT。此外，Thomas[12]发现一对反向旋转的VAWT 之间的气动干涉作用，能够提升各台VAWT 的功率系数。Whittlesey 等[13]发现，基于鱼群尾流反卡门涡街空间特征布置的VAWT 风电场，输出风能密度比HAWT风电场提升了一个量级。

综上所述，在海上风能的开发利用中，VAWT 比HAWT 更具优势。但是，相比单台HAWT 的功率系数，单台VAWT 的功率系数较低[14]。虽然已经有大量针对单台VAWT 和VAWT 风电场的研究[15]，但是提升输出功率和风能密度的最优方法还没有达成共识[7]，还有必要从叶片、单台风力机和风电场等多个尺度研究VAWT 对流动的气动响应机理。

本文针对直叶片VAWT，首先从叶片尺度出发，采用正弦俯仰振动对VAWT 转动叶片的攻角变化进行了比拟。通过对正弦俯仰振动叶片受力进行直接测量和对周围流场进行锁相粒子图像测速（particle image velocimetry,PIV），揭示了不同VAWT 操作参数下，叶片受力在不同涡脱落模式的动态响应。然后基于叶片尺度的比拟结果，从风电场阵列尺度出发，提出了一种采用双叶片VAWT 作为阵列单元的风墙构型，并简述了其主要优点。

1 VAWT 叶片攻角变化比拟

1.1 VAWT 叶片攻角变化特征

直叶片VAWT 由若干矩形截面叶片和竖直转动轴构成。图1 为一个VAWT 叶片逆时针旋转一周轨迹的俯视图，叶片方位角为θ，角速度为 ω，叶轮半径为R，叶片弦长c垂直于径向。自由来流从左至右，绝对速度为VA。叶片的尖速比定义为：

图1 不同方位角VAWT 叶片速度矢量和受力矢量Fig.1 Velocity and force vectors of the VAWT foil at different azimuth

自由来流相对叶片的速度为：

其中 ω×R为叶片转动的线速度。叶片所受的升力FL垂直于VR指向吸力侧，阻力FD沿着VR方向。

叶片攻角α为叶片弦长方向（从前缘至尾缘，与ω×R反方向）和自由来流相对速度VR的夹角。由于叶片的弦长方向和线速度都随方位角变化，所以叶片攻角α可以表达为方位角 θ的函数：

图2 展示了VAWT 叶片攻角随叶尖速比和方位角变化的曲线。如图2（a）所示，转动叶片攻角幅值可以写成叶尖速比的函数：

图2 攻角随尖速比和方位角的变化Fig.2 Variation of angle of attack with tip-speed ratio and azimuth

当λ=2时，叶片攻角幅值αm约为30°。随着 λ增大，αm逐渐减小。当 λ=10时，叶片攻角幅值αm减小至约6°。使用攻角幅值αm对叶片攻角进行无量纲化，可以看到随着 λ 增大，α/αm曲线逐渐趋近于正弦曲线：

1.2 正弦俯仰振动比拟

在VAWT 正常工作的尖速比范围内，随着尖速比增大，α/αm～ θ曲线与正弦曲线公式（7）的差别越来越小。因此，可以使用均匀来流中攻角0°附近正弦俯仰振动的叶片运动：

进行比拟。该比拟基于攻角变化的一致性，忽略了VAWT 中的离心力效应。如图3 所示，叶片选用NACA0012 翼型截面，围绕c/4转轴振动，ϕ为叶片振动相位角，其与VAWT 方位角 θ的关系为：

图3 正弦俯仰振动叶片Fig.3 Sinusoidally pitching foil

在进行比拟的叶片振动和VAWT 叶片转动两种流动中，自由来流速度、叶片攻角、攻角幅值、叶片弦长、升力和阻力的定义相同。此外，定义叶片振动的缩减频率k=πfc/VA，其中f为叶片振动频率。VAWT 方位角速度 ω与叶片振动频率f的关系为：

VAWT 叶尖速比 λ与振动叶片缩减频率k的关系为：

如果知道VAWT 中叶轮半径与叶片弦长之比R/c的值，则可以知道叶尖速比 λ与振动叶片缩减频率k的比例关系。

2 比拟实验

为了获知由叶片攻角变化引起的VAWT 叶片在不同涡脱落模式的动力响应，进行了正弦俯仰振动叶片的比拟实验。实验在上海交通大学大尺度循环水槽中进行，自由来流名义湍流度不大于2%。刚性叶片采用NACA0012 翼型截面，且叶片形状沿展向均匀分布，弦长c=0.2 m，展长S=0.8 m，完全浸没水中，见图4（a）。叶片进行正弦俯仰振动α(t)=αmsin(2πft)，振动幅值αm和频率f可独立调节。叶片、测量传感器和运动机构的布置如图4（b）所示。实验坐标系原点位于攻角0°叶片尾缘，x轴正向指向下游流向，y轴沿着横向，z轴正向沿叶片展向向上。

图4 实验设置Fig.4 Experimental setup

本实验使用测力天平测量叶片所受的动态升力和阻力，使用静态扭矩传感器独立测量叶片转矩。测量得到整个叶片的升力FL、阻力FD和转矩M之后，使用来流动压/2、弦长c和展长S计算得到对应的无量纲系数：

其中 ρ为密度。PIV 测量使用Nd-YAG 脉冲激光器，以直径55 μm 的空心玻璃球为示踪粒子，使用两个沿流向并排布置的CCD 相机，搭配Nikon AF NIKKOR镜头。PIV 图像的采集和分析，使用LaVision 公司的DaVis 软件包实现。使用伺服电机控制器的位置比较信号，实现锁相PIV 测量。不确定度分析等更多测力和PIV 测量细节请参考Zheng 等[16]的研究。

实验在VA=0.25、0.5、0.8、1.2 m/s 的自由流速下进行，振动攻角幅值设置为αm=6°、8°、12°、14°、16°、24°、34°，振动频率设为f=0.25、0.5、1.0 Hz，弦长雷诺数为Re=5.9×104、1 .2×105、1 .9×105、2.4×105。通过考察攻角0°静态叶片阻力系数，可以确定叶片表面边界层在本实验雷诺数范围内已经转捩为湍流状态。缩减频率范围为 0.16 ≤k≤2.51，对应R/c=2的VAWT 叶尖速比范围为 0.64 ≤λ ≤10.04。此外，还定义了基于叶片振动幅值的Strouhal 数：

其中A为叶片尾缘总行程。St是由k和αm合成的单一变量，用以刻画尾流涡型[17]，其范围为0.010 ≤St≤0.488 。图4（c）在k～αm相平面内标出了上述测点和作为参考的St等值线云图。

3 结果与讨论

3.1 涡脱落模式

对不同缩减频率k和攻角幅值αm的锁相PIV 涡量场进行相位平均，并使用自由来流速度VA和叶片弦长c进行无量纲化，可以识别尾流涡型结构，从而揭示不同涡脱落模式。在本实验参数空间识别到了3 种涡型结构：前缘离散涡（leading-edge vortex,LEV）、蜿蜒尾流（undulating wake,UW）和反卡门涡街（reverse von Kármán vortex street,RvKVS），对应的无量纲相位平均涡量场 〈ω〉c/VA如图5 所示。UW 涡型和RvKVS 涡型自叶片尾缘产生并脱落，而LEV 涡型从叶片前缘附近产生并在输运到尾缘之前从叶片表面脱落。在UW 涡型和RvKVS 涡型中，叶片表面边界层的涡层在叶片振动过程中一直附着在叶片表面。不同于LEV 涡型和RvKVS 涡型，UW 涡型没有离散涡，主要流动结构是具有一定曲率的涡层。两种离散涡涡型中，LEV 涡型负涡量离散涡位于横向中线以上而正涡量离散涡位于横向中线以下；RvKVS涡型中正负涡量的离散涡位置与之相反。

图5 尾流不同涡型的无量纲相位平均涡量场〈ω〉c/VAFig.5 Normalized phase-averaged vorticity fields 〈 ω〉c/VA of different vortex patterns in the wake

图6 展示了无量纲尾流特征频率fw/f随幅值Strouhal数变化的规律。尾流特征频率fw的定义请参考Zheng 等[16]的研究。UW 涡型（圆形符号）和RvKVS 涡型（方形符号）的尾流特征频率fw表现出两种相反的趋势：UW 涡型的fw从 2f单调递减至f；RvKVS 涡型的fw在f的邻域内缓慢增大，意味着RvKVS 涡型的特征频率近似锁定在叶片的俯仰振动频率。当雷诺数确定，幅值Strouhal 数可以作为一个区分UW 涡型和RvKVS 涡型的单值参数：当Re=5.9×104（深色），尾流从UW 涡型转变为RvKVS 涡型的临界St在 0.122＜St＜0.167范围内；当Re=1.2×105（浅色），尾流从UW 涡型转变为RvKVS 涡型的临界St在 0.041＜St＜0.084范围内；当Re=1 320、2 640，Schnipper 等[17]在肥皂膜流动显示实验中，得到尾流转变为RvKVS 涡型的临界St为0.18。本实验较高Re结果与Schnipper 等较低Re结果的趋势一致，临界St随Re增大而减小。此外，在Re=1.2×105的UW-RvKVS 涡型转变区域上界（浅色条带区域右边界），存在相同St的LEV 涡型参数点（浅色五角星符号）和RvKVS 涡型参数点（浅色方形符号），意味着单独使用St不足以区分LEV 涡型。

图6 无量纲尾流特征频率随Strouhal 数的变化Fig.6 Variation of nondimensional characteristic frequency in the wake withSt

图7 展示了不同雷诺数下，LEV（五角星符号）、UW（圆形符号）和RvKVS（方形符号）3 种涡型在k～αm平面内的参数点分布，符号颜色随Re增大而变浅。LEV 涡型控制低k大αm区域，UW 涡型控制低k小αm区域，RvKVS 涡型控制高k区域。由式（5）、式（6）和式（11）可推导出符合VAWT 叶片转动特点的缩减频率和攻角幅值的关系曲线。图7 的红色实线、划线和点线分别为VAWT 叶轮半径与叶片弦长之比R/c=2、3、4 的k～αm关系曲线。3 条曲线缩减频率范围均对应VAWT 叶尖速比范围 2 ≤λ ≤10，且均通过3 种涡型控制的参数区域，表明在常规操作参数下，3 种涡脱落模式均可能出现在VAWT 转动叶片附近。图7 的蓝色实线、划线和点线分别为VAWT叶尖速比 λ=2、4、6的k～αm关系曲线。在VAWT 叶片数量N一定的情况下，依据符合转动叶片运动特征的k～αm曲线分布情况，可以知道：R/c越小（红色实线），叶轮实度Nc/2R越大，在低 λ操作条件下产生LEV 涡型的可能性较小，在高 λ操作条件下产生RvKVS 涡型的可能性较大；R/c越大（红色点线），叶轮实度越小，在低 λ操作条件下产生LEV 涡型的可能性较大，在高 λ操作条件下产生RvKVS 涡型可能性较小。

图7 三种尾流涡型在k～ αm平面的控制区域Fig.7 Parameter regions dominated by three wake patterns on k ～ αm plane

3.2 动力相位响应

时变受力系数反映了叶片受力响应的动态特征，可以分解为两个分量：由叶片振动引起的确定性分量和由小尺度湍流引起的非确定性分量。为研究叶片受力对叶片运动相位的响应，对时变系数进行了相位平均。

图8 展现了不同缩减频率k和攻角幅值αm下，脉动升力和转矩系数的相位平均曲线为更好地比较，使用αm对系数进行无量纲化，使用叶片振动周期T对时间变量进行无量纲化。攻角幅值αm=6°、8°、12°、16°、24°的结果分别对应浅色实线、深色实线、黑色实线、黑色划线和黑色点线。从相位响应曲线可以看到k和αm对受力相位响应的影响：的幅值随k的增大而增大；在低缩减频率k=0.31时，受力响应曲线为正弦曲线，并随着k的增大逐渐偏离正弦曲线，其本质是叶片运动和受力响应的解耦[18]。

图8 相位平均脉动升力和转矩系数Fig.8 Phase-averaged fluctuating lift and moment coefficients

4 基于小型VAWT 单元的紧凑型风墙

相比UW 涡型和LEV 涡型，RvKVS 涡型对VAWT 叶轮旋转转矩和输出功率的提升有利，这点除了从不同涡脱落模式的受力相位响应角度（见第3 节）可以知道，从振动叶片平均受力角度也可以知道。由于俯仰振动叶片截面和运动的对称性，平均升力和转矩系数近似为0。图9 展示了平均阻力系数随缩减频率k的变化，其趋势与基于线性无粘理论的Garrick 模型[20]和Fernandez 模型[21]相符。值得注意的是，当k=2.51、αm≥12◦（St=1.5ksinαm/π ≥0.25）时，出现了<0 的情况。随着St的增大，振动叶片的涡脱落模式会从UW 涡型转变为RvKVS 涡型，叶片所受平均阻力降低，直到黏性阻力CD0被克服，叶片开始受到推力；当St进一步增大，推力也会增大[16]。对于升力型VAWT，除了叶片升力对叶轮转矩的贡献，RvKVS 涡型条件下叶片平均推力的产生也会增大叶轮转矩，从而提升VAWT 性能。

图9 平均阻力系数随缩减频率的变化Fig.9 Variation of mean drag coefficient with reduced frequency

结合上述讨论和图7 不同涡脱落模式的参数分布可知，在构型确定（R/c确定）的VAWT 中，随着尖速比 λ增大，尾流涡型向RvKVS 涡型转变，旋转叶片发生动态失速的可能性降低；同时叶片所受阻力会转变为推力，对叶轮转矩产生正向贡献，利于提升VAWT功率。而欲使VAWT 旋转叶片在尽量宽的叶尖速比λ范围产生RvKVS 涡型，要求较小的叶轮叶片几何尺度比R/c。

本文采用Nes[11]的浮式风墙概念，但将其中的HAWT 阵列单元替换成小型VAWT 单元，提出了一种基于小型VAWT 单元的紧凑型风墙构型，如图10所示。该构型在继承Nes 浮式风墙构型优点的基础上，还有以下特点：

图10 基于小型VAWT 单元的浮式风墙Fig.10 Floating wind-catching wall based on small-scale VAWT units

1）因为紧凑布置的各VAWT 功率系数下降有限甚至会有所上升[10]，基于VAWT 的浮式风墙方便实现更加紧凑的布置，能够进一步增大扫风面积，从而实现较Nes 浮式风墙更大的功率系数和输出密度。

2）因为一对反向旋转的VAWT 能够提升平均功率系数[12]，所以紧凑型浮式风墙相邻VAWT 单元设置为反向旋转。相邻单元的反向旋转设置也能够进一步提升风墙结构整体的稳定性。

3）相同材质的叶轮转动惯量是长度尺度的5 次幂量级，而在相同来流条件下所受到的转矩是长度尺度的3 次幂量级，所以自启动能力是长度尺度的−2 次幂量级。小型VAWT 叶轮转动惯量较小，在相同来流条件下，具有更好的自启动能力。

4）VAWT 旋转叶片RvKVS 涡型尾流对应较小的叶轮叶片几何尺度比R/c。在紧凑型风墙叶轮实度Nc/2R有限的条件下，要求降低叶片数量N。为兼顾风墙结构的平衡性，紧凑型VAWT 浮式风墙的风轮采用两片直叶片。

5 结论

基于攻角变化的相似性，本文通过正弦俯仰振动叶片的水动力实验，比拟研究了VAWT 旋转叶片受力在不同尾流涡脱落模式的动态响应，并提出了一种基于小型VAWT 的紧凑型风墙构型。主要结论如下：

1）在Re～O(105)范围内，正弦俯仰振动叶片尾流存在3 种涡型结构：LEV、UW 和RvKVS，分别控制k～αm平面内低k大αm区域、低k小αm区域和高k区域。升力和转矩响应幅值会随k增大。LEV 涡型导致的轻失速会在转矩响应中产生高频脉动。RvKVS 涡型的离散涡核从叶片尾缘脱落，造成升力和转矩响应在峰谷值附近发生与αm有关的偏离。此外，在尾流转变为RvKVS 涡型后，随着St进一步增大，叶片平均流向力从阻力转变为推力。

2）随着 λ增大，VAWT 叶片的k增大，而αm减小。在常规 λ条件下，3 种涡脱落模式均可能出现在VAWT 转动叶片附近。较小的叶轮叶片几何尺度比R/c，在低 λ操作条件下产生LEV 涡型的可能性较小，在高 λ操作条件下产生RvKVS 涡型的可能性较大；较大的R/c，在低 λ操作条件下产生LEV 涡型的可能性较大，在高 λ操作条件下产生RvKVS 涡型的可能性较小。LEV 涡型导致的轻动态失速，会造成VAWT叶片发生高频俯仰振动，但对叶轮转矩和VAWT 功率影响不大。RvKVS 涡型的出现，伴随VAWT 叶片升力和转矩幅值增大以及平均推力的产生，会使叶片扭矩载荷增大，但同时会使叶轮转矩和输出功率提升。

3）通过研究VAWT 叶片在不同涡脱落模式下的气动响应，提出了一种基于VAWT 的新式风墙构型。在继承Nes 浮式HAWT 风墙构型优点的基础上，每个单元较小的R/c使RvKVS 涡型较易发生，有利于提升单元输出功率；小尺度单元使风墙具有较好的自启动能力；相邻单元反向旋转进一步提升平均功率系数和结构稳定性；更紧凑的阵列布置大大提升了构型整体的功率系数和风能输出密度。

本文的比拟研究，建立在两种流动叶片攻角变化的相似性基础上，而没有考虑VAWT 在旋转过程中的离心力效应，叶片的受力响应也仅考虑了自身涡脱落模式的影响，没有考虑叶片之间的相互作用。这些问题都有待进一步研究。本文提出的新式风墙构型目前处于概念设计阶段，其工作性能和安全性的定量评估还需要进行一系列原型机测试。

致谢：感谢上海交通大学海洋工程国家重点实验室的贡建国、徐昊、魏方以和代燚工程师等在实验过程中的帮助和支持。