对转式垂直轴风力机气动性能研究

曹 宇， 韩兆龙,b， 周 岱,b,c， 雷 航

(上海交通大学 a. 船舶海洋与建筑工程学院； b. 海洋工程国家重点实验室上海； c. 高新船舶与深海开发装备协同创新中心， 上海 200240)

垂直轴风力机因控制系统简单、整体受力合理等优势在海上风电产业中应用广泛[1].但是，复杂的海洋风环境，如风速过大、风向多变会导致垂直轴风力机的转子产生极大的转矩，造成底部浮式平台失稳.对此，学者们对海上垂直轴风力机的稳定性进行了研究.徐应瑜等[2]在单柱式平台的底部附加3根约束锚链，研究了10 MW浮式风机在极端风浪条件下的运动特性.刘中柏等[3]通过控制垂荡板的边长和面积，研究了海上风电浮式基础的运动特性.但是，以上方法均会增加海上风电的工程成本.同轴对转式垂直轴风力机是在同一支撑轴上设置上、下共2个转动方向完全相反的风力机转子，通过结构优化，在保证浮式平台稳定性的基础上，控制工程的成本.

目前，针对同轴对转式垂直轴风力机的研究较少，本文的部分研究借鉴Lam等[4]关于同平面对转式垂直轴风力机的成果：该类型风力机的结构的稳定性更高，尾端脱落涡消散更快，风能利用率更高.基于该成果，采用计算流体力学(CFD)数值模拟方法，对比在不同的叶尖风速比(TSR)下，对转式与独立式垂直轴风力机的空气动力响应情况，包括转子的转矩和风力机的风能利用效率以及涡量场分布情况，进而研究对转式垂直轴风力机的气动和稳定性能.

1 模型和数值方法

1.1 计算模型

基于计算精确度和计算成本的要求，选择雷诺时均剪切应力传输RANS SSTk-ω模型作为湍流模型模拟对转式垂直轴风力机的气动性能.RANS SSTk-ω模型可以通过时均化处理将求解复杂的Navier-Stokes方程转化为求解简单的平均流方程.此外，在求解涡旋黏性时，用比耗散率ω代替耗散率ε，使得模型在低雷诺数、过渡段和高雷诺数时均有较好的结果[5]，保证了结果的精确度.RANS SSTk-ω湍流模型由湍动能k和ω两项输运方程组成.

(1)

(2)

(3)

(4)

(5)

(6)

式中：vx、vy、vz分别为x、y、z方向的速度分量；t为来流时间；μ为运动黏度；μT为涡旋黏度；σk为k的湍流普朗特数；Pk为k方程的分项；β*、α和β为参数项，且β*=0.09，α=1，β=1.5；S为广义源项；σω为ω的湍流普朗特数；F1为ω方程的协调方程.方程的具体求解过程见文献[6].

图1 对转式垂直轴风力机三维模型图Fig.1 A 3D model sketch of counter-rotating vertical axis wind turbine

图2 对转式垂直轴风力机转子的平面俯视图Fig.2 Top view of rotor of counter-rotating vertical axis wind turbine

1.2 几何模型

图1所示为对转式垂直轴风力机的三维几何模型，风力机包括同轴的上、下转子，上转子顺时针转动，下转子逆时针转动，并利用中心支撑将叶片与底部浮式平台连接.图2为上、下转子的平面俯视图，叶片的安装角(弦线与切线的夹角)β′=6°.转子直径D=2 m，高度为1.32 m；H型叶片翼型采用NACA0021，叶片高度H=1.2 m，弦长为0.265 m.由于连接杆和中心支撑对风力机气动和力学性能的影响有限，所以本文模型不包括风力机连接杆和中心支撑[7].

1.3 计算域

利用CFD数值模拟方法和SSTk-ω湍流模型，研究在风场环境下的对转式垂直轴风力机的空气动力响应，并建立计算域，如图3所示.图3(a)中，计算域前端为速度进口，来流速度v∞=8 m/s；后端为压力出口，压力大小为一个标准大气压.由于滑移壁面对流场的摩擦阻力为0，流速和涡量分布固定不变，考虑数值模拟中计算域的尺寸远小于场地试验的流场环境，且计算域外的风速和涡量分布与壁面的情况不同，所以将计算域的4个壁面均设置为滑移壁面，可以更真实地反映实际场地试验的流场情况，保证模拟的准确性.图3(b)和(c)中，计算域包括旋转域和外流场域.旋转域的直径为3.2 m，高度为 1.32 m，上、下转子的间距为0.265 m，外流场域的尺寸为 31 m×12 m×10 m.计算模型中用重叠网格模拟转子的旋转特性，外流场域为背景区域，旋转域为重叠区域，两者之间的边界为重叠界面.每1个时间步内，背景区域和重叠区域通过界面交换数据，

图3 计算域和边界的参数说明Fig.3 Parameters of computational domain and boundary

因此，重叠界面上2个区域的网格始终一致，确保了数据交换的准确性.旋转中心轴在x方向距离速度进口6 m，距离压力出口25 m，以保证模拟的准确性和尾涡扩散的完整度[8].

1.4 网格划分

计算网格的拓扑结构如图4所示.为了加快界面数据的交换速度，提高数值模拟的效率，模型的网格划分采用切割体网格形式：图4(a)中，重叠网格背景区域的外流场域基础网格尺寸设置为 0.5 m，最小网格尺寸设置为0.25 m；重叠网格重叠区域的旋转域基础网格尺寸设置为0.05 m，最小网格尺寸设置为0.025 m；重叠边界的网格尺寸设置为 0.05 m.外流场域与旋转域的网格尺寸相差较大，为了保证2个区域网格的平稳过渡，在外流场域和旋转域之间外加1个网格尺寸为0.2 m的圆柱.同时，为了研究旋转域后侧远端的尾涡发展过程，对外流场的尾涡区域进行加密处理，加密区的网格尺寸为 0.2 m.图4(b)为叶片的网格划分情况，旋转域网格的基础尺寸为0.05 m，叶片边缘的网格尺寸为 0.01 m.旋转域到叶片边缘的网格呈逐渐加密的趋势，且由于叶片前缘的曲率较大，为保证计算结果的精确性，前缘网格尺寸更小，为0.005 m.为了精确分析风力机的气动性能和叶片上涡的情况，在模型的叶片表面附加棱柱层网格，如图4(c)和(d)所示.棱柱层的总厚度为0.005 m，增长率为120%，层数为30.基于模型的雷诺数Re=1.293×105，模型首层棱柱层厚度为3×10-5m，首层网格与壁面间距的无量纲参数

y+=u*y/ν<1

式中：u*为近壁面摩擦速度；y为首层网格与壁面的间距；ν为运动黏度.以保证叶片边缘附着涡的稳定性和模拟的准确性[9].

图4 计算网格的拓扑结构图Fig.4 Topology structure of computational grids

1.5 求解器设置

采用隐式非定常流动方法求解连续性方程，采用二阶时间差分法求解动量方程.在非定常计算中，将求解器的时间步长设置为0.002 s，保证数值模拟在转子每转动2° 时进行1个时间步计算；1次时间步包含10次迭代，保证模拟的湍动能残差始终在10-5左右波动[10].整个模拟过程设置为10个周期，总时长为3.74 s，以保证模拟结果收敛完全.

2 模型验证和网格独立性分析

2.1 验证方法

通过对比数值模拟与试验得到的风能利用效率(Cp)，验证风力机模拟的精确性.Cp是衡量风力机气动性能的关键因素，计算公式为

其中：Q为风力机转矩；Ω为转子角速度；ρ为来流密度.

图5为不同TSR时，CFD数值模拟与风洞试验[11]的Cp对比曲线和相对误差曲线.其中：

TSR=ΩD/2v∞

图5(a)中，当TSR<1.80时，风洞试验与数值模拟的Cp值的相似度较高；当TSR>1.80时，Cp值的相似度较低，在曲线的最高点处，数值模拟曲线偏左的误差率达6.7%.这是由于本模型的CFD数值模拟不考虑包含在风洞试验中的塔架、支撑等连接件，同时模拟的结果忽略了机械摩擦耗能的影响，所以数值模拟的Cp值较大.即在相同的Cp下，CFD数值模拟对应的TSR更小，对应的曲线偏左.当 TSR=2.10时，数值模拟和风洞试验的Cp值均达到最大，约为0.18.图5(b)中，当TSR>1.30时，Cp的相对误差均小于12%， 模拟结果与试验结果吻合较好；当TSR=1.05时，Cp相对误差最大，约为28.5%.这是由于当TSR=1.05时，风力机处于动态失速状态，叶片边缘出现了大量脱落涡，而风力机叶片、连接杆和中心支撑等构件上的脱落涡会导致叶片升力和风力机转矩的降低，从而引起Cp值的降低.CFD数值模拟忽略了连接件、中心支撑的摩擦耗能的影响，未考虑忽略构件上的脱落涡影响，导致Cp值偏大，精度出现较大偏差.综上可知，当TSR=2.10时，Cp值最大，且Cp相对误差约为0.因此，后文数值模拟的分析均基于TSR=2.10的情况.

图5 CFD数值模拟与风洞试验的Cp值对比Fig.5 Comparison of Cp values of CFD numerical simulation and wind tunnel test

2.2 网格独立性分析

表1为TSR=2.10时， 基于Cp相对误差的网格独立性验证.随着网格由粗糙变精细(网格尺寸由大变小)，Cp值增大，Cp相对误差减小，边界层的y+值减小.在精细网格中，风洞试验的Cp相对误差和y+值均最小.网格划分越细，Cp相对误差越小，y+越小，模拟结果的收敛性和精确性越好[12].但是，随着网格由中等变精细，总网格数增加了37%，Cp相对误差仅降低了2%.因此，综合考虑模型的计算成本和精确度要求，后文研究均基于中等网格的情况.

表1 网格独立性验证Tab.1 Grid independence verification

3 分析与讨论

3.1 转矩对比

图6为TSR=2.10时，对转式与独立式垂直轴风力机的总转矩(Qt)随来流时间的对比关系.绿色框中为一个周期内总转矩的变化规律[13]，当风力机的叶片位于上风区且方位角θ=60°～120° 时，叶片升力和风力机转矩取得最大值.而本文的研究对象为两叶片垂直轴风力机，在同一周期内，2个叶片会分别转动到上风区θ=60°～120° 的位置.因此，1个周期内风力机转矩会取得2次最大值，且两值的大小相近.由图可知，独立式风力机转矩的最大值为30 N·m，最小值为-10 N·m，两值相差40 N·m.而对转式风力机的总转矩曲线在Qt=0附近小幅度波动，最大值与最小值相差小于3 N·m.可知，同等风场条件下，对转式风力机传递到浮式平台的的总转矩更小，结构的稳定性更强.这是由于对转式风力机存在完全相同的上、下转子，在对转过程中，转子传递到浮台的转矩近似为大小相等、方向相反，即转矩几乎抵消，所以总转矩在0附近波动且波动幅度很小.

图6 风力机的转矩对比Fig.6 Comparison of torques of different wind turbines

3.2 气动性能分析

图7为不同TSR的对转式与独立式垂直轴风力机的Cp对比关系.其中，Cp1、Cp2分别为对转式风力机上、下转子的风能利用效率，Cps为独立式风力机的风能利用效率.图中，Cp1和Cp2曲线基本重合，相对误差小于1%.这是由于对转式风力机的转子为反对称结构，上、下转子的气动性能相似，流场分布和脱涡情况也基本一致，所以Cp值基本相同.而Cp1、Cp2与Cps的差异明显，相对误差大于10%.当TSR=1.30时，对转式与独立式风力机的Cp曲线出现1个交叉点，交叉点两侧的Cp特征完全相反.当TSR<1.30时，Cps大于Cp1和Cp2，最大差值点在TSR=1.05处，相对误差约为55.8%.当TSR>1.30时，Cps小于Cp1和Cp2，最大差值点在TSR=2.10处，相对差值约为10.85%.

图7 转子的Cp值对比Fig.7 Comparison of Cp values of different rotors

图8为不同TSR的叶片攻角α与θ的关系曲线，图9为独立式和对转式垂直轴风力机的三维涡量图.基于图8和9，对图7中TSR<1.30时，Cps大于Cp1和Cp2的情况进行分析(取TSR=1.05时的风况).

图9 三维涡量图Fig.9 3D vorticity diagrams

图8中α与θ为正弦关系，随着TSR的增大，α曲线的幅值减小.风力机的失速情况在α>15° 时发生[14]，图中水平虚线为α=15° 的情况，点A、A′为TSR=1.05时水平虚线与关系曲线的交点，点B、B′为TSR=2.40时水平虚线与关系曲线的交点.可以看出，线段AA′明显长于线段BB′，即当TSR=1.05、α>15° 时，风力机对应的方位角的范围更大，失速时间更长.因此，TSR较小时，风力机会提前进入失速状态，且失速时间更长，涡脱落现象更严重.图9中，黑色虚线框内为风力机转子内部和叶片边缘的涡流情况.当TSR<1.30时，相较于独立式风力机，对转式风力机的叶片边缘和转子内部的脱涡现象更严重.由于叶片边缘过多的脱涡会导致叶片升力降低、阻力上升，使得风力机的升阻比减小[15]，所以当TSR<1.30时，对转式风力机的Cp值更低.

图10为TSR=2.10时，对转式与独立式风力机的转矩系数(CM)随θ变化的关系.其中，

图11为独立式风力机顺时针转动的俯视图.基于图10和11，对图7中TSR>1.30时，Cps小于Cp1、Cp2的情况进行分析(取TSR=2.10时的风况).

图10中，在θ=0°～360° 的1个周期内，独立式与对转式风力机的CM随θ的变化过程相似.CM在θ=60°，270°处存在极大值，且在θ=270° 处取得最大值.当θ=60° 时，对转式风力机的CM大于独立式风力机的CM，相对误差为54%，且对转式的CM曲线整体高于独立式风力机的CM曲线.这是由于对转式风力机转子内脱落涡的扰动更强烈，所以转矩系数更高[16].图11中，当θ=60° 时，叶片位于下风区；当θ=270° 时，叶片位于上风区.因此，当θ=270° 时，风力机的转矩系数更大.

图10 风力机的转矩系数对比Fig.10 Comparison of torque coefficients of different wind turbines

图12 转子中平面涡量图Fig.12 Vorticity in middle plane of rotors

图12为TSR=2.10、θ=60° 时，对转式风力机上转子与独立式风力机转子的中平面涡量对比关系.其中，中平面指位于风力机转子高度1/2处的平面.基于图12，对图10中θ=60° 时CM的情况进行分析.图12(a)中，对转式风力机的远端涡流长度较短且脱涡的强度较低，由外流场域前端进入的风能更多地被转子叶片吸收，转子后端的损失风能较少.这是由于风力机远端脱涡的比率和强度对转子的转矩系数有一定影响.风力机来流方向的远端涡流长度较短且脱涡强度较低，说明风力机尾部涡脱落较少，能量散失较低，从而引起风力机转矩系数增加[17].而对转式风力机尾部脱涡较少，这可能是由于在风力机对转过程中，上、下转子间的脱涡交错会产生涡流互扰，一定程度上阻碍流体通过风力机转子向流场域尾部流动，导致更多的风能被风力机转子吸收，使得尾部脱涡较少.此外，由于风能利用效率与转矩系数存在正比关系，即

Cp=(DΩ/2v∞)CM

所以，当TSR>1.30时，对转式风力机的风能利用效率更高.

4 结论

针对海上漂浮独立式垂直轴风力机的失稳问题，提出同轴对转式垂直轴风力机的设计理念.基于计算流体力学理论，借助RANS SSTk-ω湍流模型对对转式风力机进行数值模拟，并对比对转式与独立式垂直轴风力机的转矩和风能利用效率.结合涡流理论，进一步分析不同TSR时，对转式与独立式风力机的气动性能差异性的原因.主要结论如下：

(1) 同等外流场条件下，相较于独立式风力机，对转式风力机的转子传递到浮式平台的的总转矩在0附近波动且波动幅度很小，其浮式平台的稳定性更强.

(2) 当TSR<1.30时，风力机叶片进入失速状态的时间长、程度深.在长时间的失速状态下，对转式风力机的叶片边缘和转子内部脱涡更严重，风能利用效率更低.当TSR>1.30时，相较于独立式风力机，对转式风力机的远端涡流长度较短且脱涡强度较低，风能利用效率更高.

(3) 本文提出的新型结构理念和分析研究方法可以为海上漂浮垂直轴风力机的气动和稳定性能的优化提供一定的参考价值.