程杰 来永斌 王龙 武浩 李亮

摘要：针对风力机叶片的结冰会严重影响风电机组的功率输出、危害风电机组的安全运行的情况，采用计算流体力学方法研究了典型的风力机翼型的结冰特性。首先，对攻角分别为0°，4°，8°的NACA0012翼型水滴撞击特性进行了计算，计算结果与文献实验结果吻合良好，表明了计算方法的正确性。在此基础上，对不同风速和不同温度下的典型风力机叶片表面的覆冰进行了模拟研究，得到了水收集系数和结冰厚度随温度和速度的变化趋势。结果表明：风速增加导致叶片表面结冰厚度增加，结冰附着面积增加;温度降低导致叶片表面结冰厚度增加;S801翼型的抗冰效果明显大于S802翼型的抗冰效果。研究结果有利于在低温环境中合理选择翼型，为降低叶片结冰危害程度提供参考。

关键词：发电工程;风力机;翼型结冰;水滴收集系数;计算流体力学

中图分类号：V211.3文献标志码：A

CHENG Jie，LAI Yongbin，WANG Long，et al.Study on the icing characteristics of typical wind turbine[J].Journal of Hebei University of Science and Technology，2019，40（1）：9-14.Study on the icing characteristics of typical wind turbine

CHENG Jie1， LAI Yongbin1， WANG Long1，  WU Hao1， LI Liang2

（1 School of Mechanical Engineering， Anhui University of Science and Technology， Huainan， Anhui 232001， China;2 School of Mechanics and Optoelectronic Physics， Anhui University of Science and Technology，  Huainan，  Anhui 232001， China）

Abstract：The icing of wind turbine blade will seriously affect the power output of wind turbine and endanger the safe operation of wind turbine. Computational fluid dynamics （CFD） is used to study the icing characteristics of typical wind wing. Firstly， the impact characteristics of water droplets on NACA0012 type turbine under different attack angle of 0°，4°，8°are calculated. The calculated results are in good agreement with the experimental results， which shows that the calculation method is correct， and， the ice coating on the blade surface of typical wind turbine with different wind speed and temperature was simulated， and the variation trend of water collection coefficient and ice thickness with temperature and velocity was obtained. The results show that： 1） with the increase of wind speed， the ice thickness of blade surface increases and the icing adhesion area increases; 2） with the decrease of temperature， the ice thickness of blade surface increases; and 3） the anti-ice effect of S801 airfoil is obviously greater than that of S802 airfoil. The research result helps logically select turbine under low temperature conditions， providing reference for lowering the blade icing hazard.

Keywords：project for power generation; wind turbine; turbine blade icing; water drop collection coefficient; computational fluid dynamics

風力发电是一种清洁无污染、可持续的新型能源，可用来缓解对常规能源的依赖。中国风电行业起步较晚，却取得了惊人的成绩，风电场遍布全国各地。风力机通常工作在野外，工作环境恶劣，尤其在中国北方地区，风电机组经常会在低温环境下工作，并极有可能会遭遇结冰现象，增加了部件的不平衡性和疲劳程度。因此研究风力机翼型结冰特性，可为风力机防除冰提供数据支持，降低风力机在结冰条件下发生故障的可能性。

数值模拟方法是结冰研究最常用的方法，国外首先采用了数值模拟方法对结冰翼型进行了预测研究，国内数值模拟研究也在不断发展。SHIN等[1]提供了一种翼型前缘结冰表面粗糙度的测量方法，证实了光滑和粗糙冰区的存在，并指出结冰表面粗糙度受温度和空气的液态含水量的影响较大，受风速的影响不大;SON等[2]对翼型结冰增长进行了定量分析，并开发出相关的翼型结冰数值模拟软件;BOSE等[3-4]研究了小型水平轴风力机覆冰的过程， 对叶片表面的覆冰做了分析，给出了覆冰过程中风力机输出功率的变化;FORTIN 等[5]利用CFD技术对风力机叶片旋转过程中的覆冰过程进行了模拟，对覆冰后风力机叶片的性能进行了相应的分析。FU等[6]通过Fluent计算空气和水的两相流动，进行了霜冰的外形和质量积累的仿真计算。HOMOLA等[7]利用结冰软件模拟翼型结冰，并通过Fluent软件对比了结冰前后翼型气动性能的变化。KRAJ等[8]通过实验对涡轮机叶片结冰进行研究，指出叶片上冰的黏附力与附着面积与环境的温度成比例关系。马强等[9]采用Fluent中的欧拉多相流模型，通过添加质量传递以及对流换热系数计算模型等，实现了对风力机叶片表面覆冰的二维数值模拟，得到了叶片表面覆冰状况。孔满昭等[10]应用fensapice结冰计算软件，对NACA0012翼型进行了结冰过程的模拟研究，同时对翼型结冰前后的气动特性作了比较和分析，发现翼型表面在结冰后，会导致翼型前缘气流提前分离，并严重破坏其气动性能。桑为民等[11]对于霜状结冰， 在给定的结冰气候条件下， 采用数值方法， 预测整个结冰过程，并分析不同形状的冰型对翼型绕流及气动特性的影响。邵晓海等[12]采用欧拉两相流法研究了冷水滴撞击翼型表面的收集特性，结果表明，平均水滴直径越大，局部水滴收集系数和撞击区域越大。李岩等[13]利用风洞试验和数值模拟相结合的手段，研究了小型垂直轴风力机叶片在旋转状态下的结冰特性以及结冰后翼型气动特性变化，指出叶片结冰导致不同旋转角下叶片翼型周围的压力场和速度场发生了不同程度的变化，从而气动特性发生变化，影响了风力机性能。

河北科技大学学报2019年第1期程杰，等：典型风力机翼型结冰特性分析笔者采用计算流体力学方法对风力机翼型表面进行结冰情况的模拟研究，分析温度和来流速度对翼型表面冰层增长的影响，对比不同翼型在相同條件下的结冰情况。

1结冰计算的原理及过程

1.1空气流场的计算

通过求解空气流体的控制方程得到翼型的外绕流流场，采用无量纲化积分形式的N-S方程[14-15]：t∫∫ΩQdV+∫αΩF（Q）·ndS=∫αΩG（Q）·ndS，（1）式中：Ω是控制体;Ω是对应的边界;Q是守恒变量;F是无黏通量;G是黏性通量;n是控制体外法线向量，其中：Q=ρ

1.2水滴运动方程与冰形计算

分别跟踪每一个水滴的运动并做出如下假设：1）水滴在运动过程中，体积保持不变;2）水滴在撞击过程中不破裂，不溅射;3）水滴初始速度与自由来流的速度相等，水滴在运动过程中不影响流场的性质。在计算水滴运动轨迹时，考虑作用在水滴上的重力、浮力和阻力的影响，则根据牛顿第二定律，水滴运动方程可写作式（5）。Mdd2xddt2=12ρaAdCd|u-ud|（u-ud）+（ρd-ρa）Vdg。（5）上述各项含义可参见文献＼[12＼]。

结冰过程采用Messinger 控制容积思想，结冰厚度计算可参见文献＼[16＼]。

2计算结果和分析

2.1流场网格划分

使用Ansys软件划分计算区域网格，对弦长为1 m的S801和S802翼型分别进行了网格的划分，计算区域尺寸为30 m×18 m×18 m，其中入口边界距离叶片10 m，出口边界距离叶片20 m，上下边界距离叶片9 m。网格总数约为50万。网格的划分结果如图1所示。

2.2算例验证

首先对NACA012翼型进行了水滴特性的计算，翼型弦长c=0533 4 m，来流马赫数Ma=0302，攻角AOA分别选取0°，4°，8°，水滴的当量直径MVD=20 μm，空气中液态水含量为1 g/m3。图2给出了不同攻角下的水滴表观密度ρd的分布云图。从图2中可见，随着攻角的增加，上翼型表面的阴影区区域逐渐增加，下翼型表面的阴影区域范围逐渐减小。图3给出了不同攻角下翼型的水收集系数的分布，横坐标中的S为该单元的中心点到驻点的距离，本文将计算结果与文献＼[17＼]结果进行比对。图3中可见，水滴局部收集系数的计算结果与文献结果吻合良好，验证了本文计算方法的正确性。通过图中不同攻角下的水收集系数的结果对比，可以发现最大水收集系数不受攻角变化的影响，均保持在0.7附近。攻角增大，最大水收集系数点将向右偏离驻点，实际情况则是沿着叶片上表面向上移动。

2.3S801和S802翼型计算结果及分析

在0°攻角时，对温度为263.15 K，来流速度为10，15和20 m/s，以及来流速度为20 m/s、不同温度两种情况分别对S801和S802翼型进行模拟计算。模拟计算过程中其他计算参数如表1所示。

结冰时间/s弦长/m来流压强/Pa空气含水量/（g·m-3）水滴当量直径/μm4201101 325120

2.3.1不同风速下的结冰对比

局部水收集系数对结冰计算至关重要，分别对不同风速下叶片局部的水收集系数进行了对比，图4给出了S801和S802翼型在温度为263.15 K，来流速度分别为10，15和20 m/s时，其他计算参数和表1相同条件下叶片的水收集系数变化趋势。为方便比较，横坐标采用翼型表面单元的中心点到驻点的距离与翼型弦长的比值，用“S/c”表示。从图4中可以看出，随着速度的增加水收集系数也在不断增加，这与水滴计算的假设条件相关。在假设条件中，水滴不破裂，水滴不会溅射，随着速度增大水滴的惯性也会增大，水滴撞击翼型表面时的趋势也会增加，最终导致收集系数增大[18]。因为只考虑0°攻角的情况，故水收集系数最大位置点均集中在驻点附近。另外，对比S801和S802翼型的水收集系数，可以看出S802翼型的水收集系数最大值在相同的条件下，均大于S801翼型的水收集系数的最大值，两者水收集系数的最大值相差1.5倍左右。

图4水收集系数

Fig.4Coefficient of local water collection

風速是影响叶片表面结冰的重要因素，图5给出了在温度为263.15 K时，来流速度分别为10，15和20 m/s，其他参数与表1相同的条件下的S801和S802翼型前缘表面的结冰厚度变化图。从图中可以看出，随着速度不断增加翼型的最大结冰厚度也在不断增加，结冰最大位置点都在集中驻点附近。随着风速的增大，单位时间内叶片表面收集到的水滴数就越大，从而导致结冰量增加。随着风速增加，水滴在空气剪切应力和压力作用下沿翼型上下表面运动距离更远，此外低温导致水滴凝结成冰，从而冰形覆盖范围更广[19]。对比不同翼型的最大结冰厚度可以明显看出，在相同条件下S802翼型最大结冰厚度均比S801翼型的最大结冰厚度大， S801翼型前缘结冰附着面较大，S802翼型前缘结冰较集中，两种翼型的最大结冰厚度相差1.5倍左右。

2.3.2不同温度下的结冰对比

温度是影响结冰厚度的重要因素，温度的不同会对结冰速度产生影响。在来流风速为20 m/s，其他条件和表1相同，分别对S801和S802翼型在温度不同的情况下进行了模拟研究，温度条件设置为243.15，248.15，253.15，258.15，263.15和268.15 K。

图6给出了S801和S802翼型叶片表面在不同温度下结冰厚度的变化趋势，从图中可以看出，随着温度的降低翼型的表面最大结冰厚度在不断增加，在263.15～273.15 K之间冰厚会变化比较明显，计算条件温度的变化使得结冰的类型发生了改变[20]，温度的不断降低使得结冰类型从明冰向霜冰发生转变，冰面的粗糙程度不断加大，最终导致了温度在263.15和273.15 K时结冰厚度的剧烈变化。

笔者首先在算例部分给出了软件计算与试验对比结果，验证了计算过程的正确性，再经过数值模拟方法对S801翼型叶片和S802翼型叶片进行模拟，计算了不同风速下两种翼型水收集系数的变化趋势，分析了在不同温度和不同风速下的叶片前缘结冰厚度的变化，得到如下结论。

1）风速对风力机前缘叶片表面的水收集系数影响较大，随着速度的增加，叶片表面水收集系数在不断增大。对比两种翼型的计算结果可以看出，S802翼型的水收集系数要大于S801翼型。

2）风速对风力机叶片前缘表面结冰厚度有较大的影响，风速的增加会导致叶片前缘结冰厚度增加，结冰附着面积加大，且前缘表面结冰厚度变化的趋势和表面水收集系数变化趋势一致，反映出水收集系数对叶片表面结冰的影响。

3）温度是影响叶片表面结冰的另一个重要因素，温度的降低会导致叶片结冰厚度的增加。

4）不同种翼型对结冰的抗性不同，在相同条件下S801翼型的抗冰性能要优于S802翼型的抗冰性能。

本文为简化计算过程，采用二维截断模型得到了温度和来流风速对两种翼型的影响结果，为解决工程实际问题提供了参考依据。但在计算过程中，水滴假设为球形，与实际情况可能存在一定偏差，对结果有一定的影响，后续工作尚需改进。

参考文献/References：

[1]SHIN J. Characteristics of surface roughness associated with leading-edge ice accretion[J]. Journal of Aircraft， 2015， 33（2）：316-321.

[2]SON C， OH S， YEE K. Quantitative analysis of a two-dimensional ice accretion on airfoils[J]. Journal of Mechanical Science and Technology， 2012， 26（4）：1059-1071.

[3]BOSE N. Icing on a small horizontal-axis wind turbine part Ⅰ： Glaze-ice profiles[J]. Journal of Wind Engineering and Industrial Aerodynamics， 1992， 45（1）： 75-85.

[4]BOSE N. Icing on a small horizontal-axis wind tubine part Ⅱ： Three dimensional ice and wet snow formation[J].Journal of Wind Engineering and Industrial Aerodynamics，1992， 45（1）：87-96.

[5]FORTIN G， PERRON J. Wind Turbine Icing and De-Icing[R]. Florida： AIAA 2009-274，2009.

[6]FU P， FARZANEH M. A CFD approach for modeling the rime-ice accretion process on a horizontal-axis wind turbine[J]. Journal of Wind Engineering & Industrial Aerodynamics， 2010， 98（4）：181-188.

[7]HOMOLA M C， VIRK M S， WALLENIUS T， et al. Effect of atmospheric temperature and droplet size variation on ice accretion of wind turbine blades[J]. Journal of Wind Engineering & Industrial Aerodynamics， 2010， 98（12）：724-729.

[8]KRAJ A G， BIBEAU E L. Measurement method and results of ice adhesion force on the curved surface of a wind turbine blade[J]. Renewable Energy， 2010， 35（4）：741-746.

[9]马强， 吴晓敏， 陈彦. 风力机叶片覆冰的数值模拟[J]. 工程热物理学报， 2014， 35（4）：770-773.

MA Qiang， WU Xiaomin， CHEN Yan. Numerical simulation of ice accretion on wind turbine blades[J]. Journal of Engineering Thermophysics， 2014， 35（4）：770-773.

[10]孔满昭， 蒲宏斌. 翼型结冰过程数值模拟验算与分析[J]. 航空工程进展， 2011， 2（3）：267-272.

KONG Manzhao， PU Hongbin. Numerical simulation check and analysis of ice accretion on airfoil[J]. Advances in Aeronautical Science and Engineering， 2011， 2（3）：267-272.

[11]桑为民， 蒋胜矩， 李凤蔚. 翼型冰增长和结冰影响的数值模拟研究[J]. 应用力学学报， 2008， 25（3）：371-374.

SANG Weimin， JIANG Shengju， LI Fengwei. Numerical simulation for ice accretion and icing effects on airfoils[J]. Chinese Journal of Applied Mechanics， 2008， 25（3）：371-374.

[12]邵晓海， 刘星. 数值模拟过冷水滴撞击翼型表面的收集特性[J]. 应用力学学报， 2015，32（2）：266-270.

SHAO Xiaohai， LIU Xing. Numerical method to simulate collection efficiency of supercooled droplets impinging on an airfoil[J]. Chinese Journal of Applied Mechanics， 2015，32（2）：266-270.

[13]李岩， 刘钦东， 王绍龙，等. 小型垂直轴风力机叶片结冰风洞试验与数值计算[J]. 空气动力学学报， 2016， 34（5）：568-572.

LI Yan， LIU Qiandong， WANG Shaolong， et al. Wind tunnel test and numerical simulation on blade icing of small-scaled vertical axis wind turbine[J].Acta Aerodynamica Sinica， 2016， 34（5）：568-572.

[14]李鑫， 白俊强， 王昆，等. 欧拉两相流法在翼型积冰中的应用[J]. 中国科学：物理学 力学 天文学， 2014， 44（3）：258-266.

LI Xin， BAI Junqiang， WANG Kun， et al. The application of eulerian two-phase flow method in airfoil ice accretion[J]. Scientia Sinica Physica， Mechanica & Astronomica， 2014， 44（3）：258-266.

[15]朱程香， 王珑， 孙志国，等. 风力机叶片翼型的结冰数值模拟研究[J]. 空气动力学学报， 2011， 29（4）：522-528.

ZHU Chengxiang， WANG Long， SUN Zhiguo， et al. Numerical study of wind turbine blade airfoil ice accretion[J]. Acta Aerodynamica Sinica， 2011， 29（4）：522-528.

[16]GENT R W， DART N P， CANSDALE J T. Aircraft icing[J]. Philosophical Transactions Mathematical Physical & Engineering Sciences， 2000， 358：2873-2911.

[17]KIM J W， DENNIS P G， SANKAR L N， et al. Ice Accretion Modeling using an Eulerian Approach for Droplet Impingement[R]. Reston： AIAA 2013-0246，2013.

[18]鲁天， 桑为民， 刘晓宇，等. 三维翼面结冰过程及其影响数值研究[J]. 航空动力学报， 2014， 29（6）：1339-1345.

LU Tian， SANG Weiming， LIU Xiaoyu， et al. Numerical simulation of icing accretion and influence on three-dimensional wing configuration[J]. Journal of Aerospace Power， 2014， 29（6）：1339-1345.

[19]陳晓圆. 飞机结冰多参数影响研究[D]. 南京：南京航空航天大学， 2015.

CHEN Xiaoyuan. Study of Multi Parameter Influence of Aircraft Icing[D]. Nanjing： Nanjing University of Aeronautics and Astronautics， 2015.

[20]任鹏飞. 结冰风力机叶片的空气动力学特性数值研究[D]. 北京：中国科学院研究生院（工程热物理研究所）， 2014.

REN Pengfei. Numerical Study on Aerodynamics of Wind Turbine Blades with Icing[D]. Beijing： Institute of Engineering Thermophysics Chinese Academy of Sciences， 2014.第40卷第1期河北科技大学学报Vol.40，No.1

2019年2月Journal of Hebei University of Science and TechnologyFeb. 2019