不同叶尖速比下VAWT的叶片冰型及性能影响

Blade Ice Shape and Performance Influence of VAWT at Different Tip Speed Ratios

摘要: 风力机叶片结冰后会造成大量的能量损失，本文在Fensap-ice中应用多时间步的方法，模拟静叶片表面的结冰，验证结果表明，本文冰型计算的方法具有准确性，模拟得到不同相对风速下的冰型。然后，构建垂直轴风力机整机模型，应用Ansys-cfx计算不同叶尖速比下垂直轴风力机结冰前后的功率系数，发现不同叶尖速比下结冰对功率系数的影响不同，随着叶尖速比的增大，功率系数下降程度增加。

Abstract: Wind turbine blades will cause a lot of energy loss after icing. In this paper, the multi-time step method is applied in Fensap-ice to simulate the icing on the surface of static blades. The verification results show that the method of ice shape calculation in this paper is accurate, and the ice shapes under different relative wind speeds are simulated. Then, the whole model of vertical axis wind turbine is built, and the power coefficient of vertical axis wind turbine before and after icing under different tip speed ratios is calculated by Ansys-cfx. It is found that icing has different effects on the power coefficient under different tip speed ratios, and the power coefficient decreases with the in-crease of tip speed ratios.

1. 引言

随着大量化石能源的使用，全球变暖、化石燃料枯竭等问题日益严重，2021年10月24号，中共中央国务院印发《关于完整准确全面贯彻新发展理念做好碳达峰碳中和工作的意见》。文件提出：积极发展非化石能源，实施可再生能源替代行动 [1]。在所有可再生能源中，风能被认为是最具成本效益的能源，增长迅速 [2]。

根据旋转轴，风力机可以分为水平轴风力机和垂直轴风力机。相比于水平轴风力机，垂直轴风力机外形简单、制造及维修成本低、系统稳定、无需偏航、气动噪声低等优点 [3]，目前是一个研究热点。

然而，对于安装在寒冷地区的风力机，在寒冷潮湿的气候条件下，叶片结冰是一个非常严重的问题，影响风力机的输出性能，并导致安全事故 [4]。对于叶片表面的结冰类型，受外界环境及自身工况的影响，主要包括温度、压强、风速、液滴水含量、水滴大小、叶片曲率、弦长、水滴收集系数和结冰时间等 [5]：当水滴撞击到叶片表面，在驻点处收集并形成回流水，在对流冷却区域回流水结冰并形成冰角，此种冰型称为明冰 [6]；水滴撞击后，立即在叶片表面结冰，没有形成不规则的冰角，此种冰型称为霜冰 [7]。

国内外对于风力机积冰的模型研究多是在飞机积冰的基础上开展的，早期研究主要采用风洞试验，随着CFD技术的迅速发展，关于积冰的数值模拟计算也开始兴起，Lewice、Turbice、Fensapice等结冰数值模拟软件都相继发展，本文研究采用的是Fensap-ice软件 [8] 进行流场计算、水滴撞击、积冰的综合计算，在计算过程中，对水滴做如下假设：水滴在空气中分布均匀，运动过程中尺寸保持不变，水滴的重力、表观质量力等于粘性阻力相比可忽略 [9]。

叶片积冰后，叶片表面气流边界层分离，会改变风机的气动性能 [10]，引起叶片较大的振动甚至叶片之间的共振，造成机组的发电效率降低 [11]。目前，国内外通过数值模拟的方法对于风力机结冰的影响研究主要借鉴航空领域，集中于单叶片的研究 [12] [13]；对于叶片结冰之后在垂直轴风力机整机中的研究鲜少。Rustem Manatbayev等人应用静态叶片结冰的方法，研究了不同冰型对VAWT的性能影响，发现霜冰对气动性能没有显著影响，而明冰会损失VAWT 60%的功率性能 [14]。

研究冬季湿冷地区不同叶尖速比下垂直轴风力机整机的结冰及气动载荷分析，会对防冰、除冰提供有效指导。为此，本文的研究应用静态叶片结冰方法 [14]，内容分为两部分：一、叶尖速比不同，即对叶片而言，不同相对风速下的结冰冰型研究；二、不同叶尖速比下，垂直轴风力机叶片结冰后的性能影响。

2. 计算模型与方法

数值模型验证

叶片结冰模拟采用Fensap-ice软件，其计算域如图1所示，整体为C型，上游设定为R = 20c的半圆，下游为L = 40C的正方形流场域。湍流模型采用Spalart-Allmaras一方程模型，采用欧拉法计算液滴撞击特性 [15]；采用多时间步长的计算方法，将总结冰过程分为多次循环叠加，每次循环中包括流场计算、水滴撞击特性、结冰计算，且进行粗糙度及网格的更新。

Figure 1. Computational domain

图1. 计算域示意图

根据文献 [16] NASA报告中的实验结果进行数值模拟验证。实验是以NACA0012翼型为研究对象，翼型弦长(C) 533.4 mm，攻角(α)为4˚，温度(T)为262.04 K，来流风速(V)为102.8 m/s，液滴水含量(LWC)为0.55 g/m³，水滴直径(MVD)为20 μm，结冰时间(t) 420 s。

采用 Pointwise进行网格绘制，整体网格为C型，采用结构化网格，如图2所示。为了确保网格无关性，采用了四套不同加密网格，分别为3万，7万，13万，20万，对应的y+分别为5.9、2.5、1.1、0.94。以NACA0012翼型升力系数为监测参数，对应值分别为1.20、1.09、1.00、0.99，对比来看，综合精确度与计算量，本文选用第3套网格。

(a) 整体计算域网格 (b) 翼型周围网格

Figure 2. Mesh for numerical simulation

图2. 数值计算网格

依据以上条件，图3将本文模拟NACA0012翼型结冰形状与实验形状 [16] 以及用LEWICE程序模拟结冰的形状 [17] 进行对比，从中可知，本文的模拟结果与实验值吻合基本吻合，证实了本文数值模拟计算方法的正确性。

风力机整机计算域俯视图如图4所示，由一个20 D × 30 D的矩形静止域，一个环形转动域(其中包含两根叶片，位于0˚和18˚相位角)和一个环形静止域组成。环形旋转域，外径为1.5 D，内径为0.5 D；环形静止域，外径为0.5 D，内径为0.04 D (中心轴直径)；速度入口和压力出口分别距旋转中心5 D，25 D。D为风轮直径，D = 1 m。

Figure 3. Validation for NACA0012 airfoil icing

图3. 冰型数值验证

Figure 4. Computational domain of VAWT

图4. 垂直轴风力机计算域

3. 计算结果与分析

3.1. 相对风速对叶片结冰冰型的影响

从垂直轴风力机的实际出发，风机尺寸普遍较小，安装环境中的风速相对较低。模拟研究中选用弦长为60 mm的叶片，风机来流速度保持9.3 m/s，静态叶片结冰攻角保持0˚，其余结冰参数与文献 [15] 保持一致，温度(T)为262.04 K，液滴水含量(LWC)为0.55 g/m³，水滴直径(MVD)为20 μm，结冰时间(t) 420 s。

根据风机叶尖速比及来流风速，根据公式(1)计算气流经过静叶片的相对风速，算例具体参数如表1所示。

$W=V\sqrt{{\sin }^2\theta +{\left(\lambda +\cos \theta \right)}^2}$ (1)

式(1)中，V——来流风速，m·s⁻¹；θ——相位角，˚； $\lambda$ ——叶尖速比。

图5显示了静叶片在不同相对速度下的积冰形状，可以观察到，在0˚攻角下，叶片结冰后，冰型对称分布在翼型两侧，从驻点到两侧的结冰厚度逐渐减小。图5(a)~(c)中，冰型曲线比较光滑，符合霜冰的特征，图5(d)中可以观察到驻点两侧已经开始出现角冰，冰型类似明冰。4个算例的工况保持一致，归因此种现象是由于相对风速不同，导致水滴收集系数不同，进而影响冰型的类型。在较高的相对风速下，水滴收集系数较高，水滴撞击到叶片表面存在溢流，驻点两侧的局部对流系数较大，较多的溢流水就会在驻点两侧的位置冻结，形成双角冰 [18]。

从图中可以观察到，相对风速不同，翼型表面结冰存在很大的差异。图中横纵坐标代表翼型前缘X、Y相对于弦长的无量纲坐标，h为驻点处的结冰厚度，s为结冰宽度。对比图5(a)~(d)，发现随着相对风速的增加，结冰厚度越厚，结冰宽度越宽。厚度依次增加2.6 mm、2.2 mm、1.8 mm，厚度依次增加2.9 mm、2.7 mm、2.4 mm，增加幅度均稳定减小。

图6表示了叶片表面的积冰质量和相对风速的关系，发现在本文研究风速范围内，积冰质量随相对风速的增加而增加，两者为正相关的关系。

Table 1. Example parameters

表1. 算例参数

(a) (b) (c) (d)

Figure 5. Ice types under different relative wind speeds

图5. 不同相对风速下的冰型

Figure 6. Variation of leaf icing quality with relative wind speed

图6. 叶片积冰质量随相对风速的变化

3.2. 叶尖速比对结冰后垂直轴风力机的性能影响

为了分析积冰对垂直轴风力机性能的影响，使用Ansys-cfx进行了数值分析。上节中从Fensap-ice中获得的静态叶片的冰型被用作Ansys-cfx的输入，研究在不同叶尖速比下叶片积冰对垂直轴风力机性能输出的影响。

在Ansys-cfx中，进行非定常数值模拟，以定常流场作为初始流场，采用二阶格式对所有方程进行空间和时间的离散，湍流模型选用SST_γ-Re_θ模型 [19]。采用双时间步进行计算，每个物理时间步为转动域转动2.5˚，每个时间步迭代步数为20，残差收敛标准为1.0 ´ 10⁻⁵且获得周期性扭矩值。经过20圈后得到收敛解，取后5圈进行平均，比较垂直轴风力机的功率系数。

为了评估结冰对垂直轴风力机性能的影响，比较结冰前后风机的功率系数。叶片结冰后，垂直轴风力机的功率系数相应降低，随着叶尖速比的增加，叶片相对风速增加，叶片结冰面积增大，功率系数降低越明显。这主要是由于叶片形状发生严重变化导致叶片表面气流分离严重，改变了风机的气动性能 [10]。在图7中，当风机系数为0.12时，结冰前，风机的功率系数随叶尖速比先增大后减小，存在最优叶尖速比，为5；结冰后，风机功率系数随叶尖速增大而减小，如图中红色曲线，功率系数的减小幅度随叶尖速比的增大而增大。

比较不同算例运行20圈后的瞬时流场的涡量图，如图8所示，叶片结冰后，随着叶尖速比的增大，叶片周围气流的分离程度加大，流场中涡量的损失增大。

Figure 7. Variation of power coefficient of VAWT with tip speed ratio

图7. VAWT的功率系数随叶尖速比的变化

Figure 8. Vorticity diagram of wind turbine flow field before and after icing under different tip speed ratios

图8. 不同叶尖速比下结冰前后风轮流场涡量图

图9中比较不同叶尖速比下结冰前后的垂直轴风力机的功率系数随相位角的变化。发现在叶尖速比

(a) (b) (c) (d)

Figure 9. Variation of power coefficient with phase angle before and after VAWT icing under different tip speed ratios

图9. 不同叶尖速比下VAWT结冰前后功率系数随相位角的变化

为1时，叶片结冰对风机整体性能影响不大，甚至在50˚~130˚相位角内还有一定功率的提升。在叶尖速比为3时，风机转动一圈，从0˚到360˚相位角内，叶片结冰后功率系数均存在一定程度的降低，但功率系数随相位角变化的趋势并未收到结冰的影响。随着叶尖速比的增大，当叶尖速比为5和7时，叶片结冰使垂直轴风力机功率系数产生显著下降，除此之外，功率系数随相位角的变化趋势发生翻天覆地的变化，结冰前在90˚和180˚相位角时，为功率系数的峰值，在0˚、180˚和360˚相位角时为功率系数的谷值；而结冰后，正好相反，峰值与谷值所对应的相位角发生了互换。

4. 结论

本文运用数值模拟的方法，研究相对风速对叶片结冰的影响及垂直轴风力机功率系数降低随叶尖速比的变化，得到以下结论：

1) 本文应用公式，计算不同叶尖速比下对应的叶片相对风速，发现随着相对风速的增加，叶片表面的积冰质量随之等比例增加。

2) 在相同的工况条件下，随着叶尖速比的增加，垂直轴风力机结冰后的性能下降更为明显。

3) 在高叶尖速比下，结冰除了会降低相应相位角下的功率系数，还会改变功率系数随相位角变化的趋势。