上海交通大学学报

上海交通大学学报, 2022, 56(9): 1176-1187 doi: 10.16183/j.cnki.jsjtu.2021.258

船舶海洋与建筑工程

风机叶片结冰对其一体化结构动态响应影响的数值分析

闯振菊,, 李春郑, 刘社文

大连海事大学 船舶与海洋工程学院,辽宁 大连 116026

Numerical Analysis of Influence of Blade Icing on Dynamic Response of Integrated Wind Turbine Structure

CHUANG Zhenju,, LI Chunzheng, LIU Shewen

College of Naval Architecture and Ocean Engineering, Dalian Maritime University, Dalian 116026, Liaoning, China

责任编辑: 石易文

收稿日期: 2021-07-14  

基金资助: 大连市科技创新基金重点学科重大项目(2020JJ25CY016)
辽宁省航运联合基金(2020-HYLH-31)
中央高校基本科研业务费专项资金(3132019306)
中央高校基本科研业务费专项资金(3132021121)

Received: 2021-07-14  

作者简介 About authors

闯振菊(1983-),女,辽宁省本溪市人,副教授,研究方向为海上新能源及海洋工程结构;E-mail:zhenjuchuang@dlmu.edu.cn.

摘要

基于美国可再生能源实验室的导管架式一体化海上风机模型,将计算流体力学(CFD)方法与风机一体化分析方法进行耦合,研究叶片结冰过程及其结冰后对风机整体动态性能的影响.首先将一体化分析得到的叶片运动姿态输入到CFD中,采用离散多相模型和融化-凝固模型进行三维风机叶片覆冰增长模拟,然后应用k-ε湍流模型计算结冰前后的气动性能,最后将叶片覆冰后的气动载荷结果返回到一体化分析方法中,分析叶片结冰对风机整体响应产生的影响.计算结果表明,叶片沿叶展方向结冰呈线性增加,结冰主要集中在叶片前缘,叶尖处积冰最厚;覆冰后叶片各剖面翼型升力系数降低、阻力系数升高.叶片结冰会降低整机功率、转子的转矩和转速,叶尖和塔顶产生额外振动响应,风机达到额定功率所需风速增大.

关键词： 冰区风电; 叶片结冰; 风机一体化分析; 功率损失; 耦合分析

Abstract

Based on the integrated jacket-support offshore wind turbine model of the National Renewable Energy Laboratory (NREL), the computational fluid dynamics (CFD) method is coupled with the wind turbine integrated analysis method to study the blade icing process and its influence on the overall dynamic performance of the wind turbine. First, the blade motion attitude calculated by the integrated analysis method is input into CFD. The discrete multiphase model and melting solidification model are used to simulate the icing growth of three-dimensional wind turbine blades. The k-ε turbulence model is used to calculate the aerodynamic performance before and after icing. Finally, the aerodynamic results after blade icing are returned to the integrated analysis method to analyze the influence of blade icing on the overall response of the wind turbine. The results show that the blade icing increases linearly along the blade span. The icing is mainly concentrated on the leading edge of the blade with the thickest ice accumulation at the tip. The lift coefficient decreases and the drag coefficient increases after icing. Blade icing will reduce the power of the whole machine, the torque, and the rotor speed. At the same time, it will lead to additional vibration response at the blade tip and tower top, and increase the wind speed required by the wind turbine to reach the rated power.

Keywords： wind power in ice area; blade icing; analysis of wind turbine integration; power loss; coupling analysis

PDF (10181KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

闯振菊, 李春郑, 刘社文. 风机叶片结冰对其一体化结构动态响应影响的数值分析[J]. 上海交通大学学报, 2022, 56(9): 1176-1187 doi:10.16183/j.cnki.jsjtu.2021.258

CHUANG Zhenju, LI Chunzheng, LIU Shewen. Numerical Analysis of Influence of Blade Icing on Dynamic Response of Integrated Wind Turbine Structure[J]. Journal of Shanghai Jiaotong University, 2022, 56(9): 1176-1187 doi:10.16183/j.cnki.jsjtu.2021.258

风能作为可再生清洁能源,近年来在全球范围内得到了快速发展.而寒区蕴含着大量的优质风能,到2020年,寒冷气候地区风电装机量预计达到 186 GW[1].但因海上低温潮湿的气候特征而导致的叶片结冰问题,是在冰区开发风能源要解决的关键问题之一.叶片结冰会改变叶片表面形状,影响叶片气动性能,降低发电效率[2].同时,导致叶片产生额外振动,增加风机运行中的疲劳负载,甚至降低风机的使用寿命[3].另外,粗糙的结冰表面会增加气动噪声,冰块脱落也会对环境和机械设备造成巨大危害.渤海冬季冰面上的大气相对湿度约为50%,而水面上的相对湿度约为90%[4],易造成结构结冰.因此,若在渤海海域建设海上风力发电机,则风机结构的结冰问题不容忽视.

为减小风机叶片结冰带来的负面影响,有必要对其结冰机理、结冰过程及特性进行研究.已有针对风机叶片结冰的研究是在飞机机翼结冰的基础上逐渐开展起来的,前期主要通过大量的风洞试验来研究结冰的形成及影响因素.文献[5]在美国国家航空航天局(NASA)的结冰研究隧道中对NACA0012翼型进行了结冰实验,研究了液滴含量、温度等不同结冰条件下形成的不同冰型.文献[6]设计旋转实验台研究了不同条件下的美国可再生能源实验室(NREL)试验叶片的结冰情况,发现攻角、温度、水滴含量和结冰时间影响叶片结冰程度,结冰程度越严重造成的风机转矩越大.

随着计算机和计算流体力学(CFD)的发展,数值方法也成为了研究叶片结冰的主要手段之一.文献[7]采用Fluent软件和气动弹性程序FAST软件,研究了不对称结冰和对称结冰对风力机载荷的影响.文献[8]利用FENSAP-ICE软件对复杂的风机结冰时间进行数值模拟.文献[9]模拟了风机叶片的结冰,发现结冰量随风速增大而增多,结冰集中在叶尖区域,结冰降低了风机叶片的气动性能.文献[10]采用拉格朗日法计算叶片覆冰,模拟了风机的雨淞和雾凇覆冰过程,并研究结冰后风机转矩及功率损失情况.文献[11]进行风力机叶片雨淞覆冰的三维数值模拟,并采用试验方法研究温度对覆冰的影响以及覆冰对风力机空气动力特性的影响.文献[12]对霜冰条件下风力机翼型结冰的数值计算进行预测,研究攻角、水滴直径、来流速度和液态水含量对翼型结冰的影响.文献[13]基于自由尾流提升线模型和有限面积法开发风机结冰计算模型,将三维结冰问题转化为二维条件进行计算,通过对NACA0012翼型和NREL风机进行计算,验证了模型的有效性;并且为模拟偏航结冰过程,提出一种改进的多点结冰计算模型,并进行NACA0012机翼风洞结冰实验验证该模型,发现偏航结冰特征导致不同的空气动力性能,这与流分离密切相关[14].

虽然国内外以飞机翼型结冰为基础对风机叶片结冰问题进行了大量研究,但风机叶片结冰对整机的影响仍处于技术空白区.本文基于美国可再生能源实验室“Offshore Code Comparison Collaboration Continuation(OC4)”项目的导管架式一体化海上风机模型,将计算流体力学与风机多体动力学计算方法进行耦合分析,研究叶片结冰过程及其结冰后对风机整体动态性能的影响.所建立的CFD与风机一体化分析耦合方法可为覆冰状态下整机动态响应特征的研究以及冰区一体化海上风机的安全性能评估提供参考依据.

1 数值模拟方法

CFD方法采用Star CCM+软件模拟计算风机叶片结冰情况.利用离散多相模型模拟空气中水滴的运动情况,以及液膜的融化-凝固模型模拟叶片积冰增长,采用k-ε湍流模型计算气动性能.风机多体动力学主要通过风机计算软件FAST[15]来完成后续计算.

1.1 CFD结冰计算理论模型

1.1.1 水滴运动模型

采用离散多相模型将拉格朗日多相和欧拉多相模型的多个特性相结合,以欧拉法模拟空气中的水滴运动,以水滴占空气的体积分数来反映水滴在求解域中的分布.而空气为连续相,使用典型的单相模型进行求解.空气中水滴体积分数较低,水滴相对空气相的影响可以忽略,故离散多相模型采用单项耦合模式,此时离散多相不会计算连续相的空气, 此方法与拉格朗日多相模型一致.而离散多相模型建立的水滴运动计算的连续性方程和动量方程[16]可分别表示为

∂∂t∫Vφ(αd)ρddV+∮Aφ(αd)ρdvd·dA=0
∂∂t∫Vφ(αd)ρdvddV+∮Aφ(αd)ρd(vd·vd)·dA=Fad+∫Vφ(αd)ρdgdV

式中:φαd为水滴的体积分数;ρd为水滴密度;vd为水滴速度;t为时间;V为体积;A为面积;g为重力加速度;Fad为空气作用于水滴的阻力,可以表示为

Fad=12CD6φ(αd)4Dρdvrvr

式中: D为水滴直径;vr=vd-va为水滴和空气之间的相对速度,va为空气速度;CD为阻力系数.根据Schiller-Naumann方法[17],阻力系数计算如下:

CD=24Re1+16Re2/3,Re≤10000.424,Re>1000

式中: Re为雷诺数,可以表示为

Re=ρavrDμa

式中:ρa为空气密度;μa为空气动力黏度.在流场运动过程中,空气中水滴会与风机叶片表面发生碰撞,用β表征叶片表面处撞击水滴的收集能力,即水滴收集系数:

β=φ(αd)φ(α∞)vd·nv∞

式中:n为叶片表面法矢量;v∞为来流速度;φ(α∞)为来流水滴体积分数.

1.1.2 结冰增长模型

空气中水滴撞击到叶片后不完全冻结,部分水滴形成液膜发生流动现象;温度较低时,水滴撞击到翼型表面立即冻结.本文选取前一种结冰形式,水滴不完全冻结生成液膜流动,结冰情况考虑温度项的影响,结冰过程中存在热力学变化.基于Messinger方法[18]的液膜质量和能量守恒方程如下:

∂∂t∫VρfdV+∫Aρfvf·dA=mimp-mice-∂∂t∫VΔmahfdV
mimp=βφ(α∞)v∞=ρdφ(αd)vd·n
∂∂t∫VρfEfdV+∫AρfHfvf·dA=∫Aqf·dA+Qimp+Qice+∂∂t∫VSEhfdV

式中:ρf为叶片上液膜密度;vf为液膜速度;hf为液膜厚度;mimp为撞击到叶片表面的水滴质量;mice为结冰量;Δma为单位面积内由液膜蒸发、剥离等引起的质量变化;Ef为液膜总能量;Hf为液膜总焓;qf为液膜热通量;Qimp为撞击叶片的水滴能量;Qice为水结冰释放的能量;SE为单位面积内由液膜蒸发、剥离、摩擦损失等引起的能量变化.

由液膜或直接撞击生成的结冰增长厚度的表达式如下:

Δhice=miceΔtρiceΔS

式中:ρice为叶片上结成的冰密度;Δt为单位时间;ΔS为单位面积.

随着结冰厚度的增加,水滴部分不再和翼型表面接触,而是与冰层表面接触.结冰增长模型主要通过液膜的质量和能量守恒来控制结冰形成,结冰后将冰等效成结构表面,此时忽略水滴与冰层接触对结果的影响,水滴的冻结条件不发生改变.

1.2 风机一体化分析方法

风机一体化分析方法[19]采用Kane’s dynamics方法[20],通过牛顿运动定律直接推导得到,Kane运动方程由广义主动力和广义惯性力两部分组成,可以表示为

Fr+Fr*=0(r=1,2,…,P)

式中:Fr为广义主动力;Fr*为广义惯性力;r为广义坐标系; P为广义坐标系的总个数.Fr和Fr*都是作用于刚体质心点Xi上的力,表达式如下:

Fr=∑i=1wvNi·FNi+ωNi·MNi

Fr*=∑i=1wvNi·-miaNi+ωNi·-H·Nir=1,2,…,P(13)

式中: w为刚体组数量;Ni为一个刚体组中刚体的个数;FNi和MNi分别为施加在质心点Xi上的主动力和力矩;mi为刚体质量;aNi为质心点Xi上加速度;vNi为线速度;ωNi为角速度;H·Ni为刚体Ni关于质心点Xi的角动量时间导数.

2 数值模型描述

2.1 一体化导管架式海上风机模型

本文主要选取OC4项目[21]中的固定导管架式海上风机,如图1所示,对其进行风机叶片结冰研究.图1(a)为导管架风机整体结构,整个一体化风机系统由导管架基础、塔筒、机舱和3个叶片组成,导管架基础的4个支腿由插入海床固定的桩腿支撑,4层X型导管用于加固垂直支腿.图1(b)为NREL 5 MW风机叶片,叶片沿叶展方向分为19个截面,由8种类型的翼型组成,其中:Cylinder为圆形翼型;DU为代尔夫特大学(Delft University)翼型;NACA为美国国家航空咨询委员会(National Advisory Committee for Aeronautics)翼型;A17为翼型展弦比为17[22].风机的主要参数如表1所示,叶片具体参数如表2所示.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   导管架式海上风力发电机

Fig.1   Jacket-structure offshore wind turbine


表1   风力发电机的主要参数

Tab.1  Main parameters of wind turbinem

名称取值
水深50
导管架基础高度70
叶轮直径126
塔高68
轮毂直径3

新窗口打开| 下载CSV


表2   叶片结构分布和翼型特性

Tab.2  Blade structure distribution and airfoil characteristics

叶展位置/m弦长/m扭转角/(°)翼型类型
03.54213.308Cylinder1
1.36673.54213.308Cylinder1
4.13.85413.308Cylinder1
6.83334.16713.308Cylinder2
10.254.55713.308DU40_A17
14.354.65211.48DU35_A17
18.454.45810.162DU35_A17
22.554.2499.011DU30_A17
26.654.0077.795DU25_A17
30.753.7486.544DU25_A17
34.853.5025.361DU21_A17
38.953.2564.188DU21_A17
43.053.013.125NACA64_A17
47.152.7642.319NACA64_A17
51.252.5181.526NACA64_A17
54.66672.3130.863NACA64_A17
57.42.0860.37NACA64_A17
60.13331.4190.106NACA64_A17
61.51.4190.106NACA64_A17

新窗口打开| 下载CSV


2.2 CFD与风机多体动力学计算方法耦合系统

本文在叶片覆冰问题计算过程中将CFD方法与风机多体动力学方法进行双向耦合.

首先,采用风机多体动力学方法进行风机整体计算,将风机叶片的运动计算结果传递到CFD方法中,根据叶片运动情况,采用离散多相(DMP)和融化-凝固模型进行叶片覆冰的仿真模拟,并利用k-ε湍流模型计算叶片气动性能.然后,将CFD方法计算的叶片覆冰结果耦合到多体动力学方法中,包括各剖面翼型结冰后的弦长(c)、表面形状、积冰厚度、升力阻力系数和改变后的气动中心位置.最后,进行结冰后的风机一体化分析,研究叶片覆冰对风机整机效率、发电功率、叶片受力和振动、转子的转矩和转速、风机塔筒振动等性能的影响,具体流程如图2所示.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   结冰仿真计算流程示意图

Fig.2   Flow chart of icing simulation calculation


3 计算结果与分析

3.1 叶片结冰的计算结果

3.1.1 计算域网格及收敛性分析

选取NACA0012翼型进行网格收敛性分析,在Star CCM+中采用多面体网格划分3种网格进行结冰计算分析.网格数量及计算工况[23]的选择如表3所示.

表3   网格数量及工况[23]

Tab.3  Number of grids and working conditions[23]

网格
类型		网格
数量		攻角
α'/(°)		来流速度/
(m·s-1)		温度/
℃		水滴含量/
(g·m-3)		水滴直径/
μm		时间/s
1586433.567.1-6.421.020360
21749933.567.1-6.421.020360
34937783.567.1-6.421.020360

新窗口打开| 下载CSV


图3为3种网格划分情况和3种网格结冰前后轮廓对比图,其中: x/c为翼型在x方向上的分布与弦长之比;y/c为翼型在y方向上的分布与弦长之比.对3种类型的网格进行结冰计算后,结冰轮廓图基本一致.考虑三维风机结冰计算采用1∶1模型,为保证计算效率和计算精度,加快计算速度,现选取网格2进行三维风机结冰计算的网格划分,包括速度进口、压力出口、远场、旋转域和风机叶片表面,网格总数量为 3186650,如图4所示.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   网格划分及结冰轮廓对比

Fig.3   Meshing and comparison of icing profiles


图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   计算域网格

Fig.4   Grid of computational domain


3.1.2 结冰计算工况选择

风机在运行过程中,沿叶展方向从叶根到叶梢,风机叶片局部线速度逐渐递增,速度越快造成的结冰情况越严重[24].而二维的翼型无法真实反映风机旋转造成的结冰效果,因此本文对NREL 5 MW风机叶片进行三维结冰数值仿真.将一体化风机计算的叶片运动响应传递到CFD方法中进行计算,在此采用风机达到刚额定功率时的风速(v)和叶片转速.以冬季渤海海域海上气象环境[4]为计算工况,确定水滴平均直径和空气中的水滴含量,根据现实风机叶片结冰情况设置结冰时间[3],具体参数如表4所示.

表4   叶片结冰计算工况

Tab.4  Blade icing calculation conditions

计算参数数值
环境温度/K258
风速/(m·s-1)11.4
水滴平均直径/μm30
空气中水滴含量/(g·m-3)8
叶片转速/(r·min-1)12.1
结冰时间/min300

新窗口打开| 下载CSV


3.1.3 结冰计算结果

CFD计算过程中忽略塔筒和风机基础的大气结冰影响,重点考虑叶片在旋转运行过程中的结冰现象.根据离散多相模型,采用欧拉法计算水滴的运动过程.将叶片表面结冰后形成的冰层固化到叶片上,叶片表面网格根据固化冰层向外变形,以此来采用动网格方法模拟叶片的结冰过程,风机尾流涡量云图和水滴撞击量计算结果分别如图5和6所示.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   三维旋转风机覆冰涡量云图

Fig.5   Three dimensional vorticity cloud diagram of icing on rotating wind turbine


图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   空气水滴在叶片上的撞击量

Fig.6   Impact of air water drop on blade


水滴撞击量表示单位时间内撞击在叶片表面上的水滴质量.由图6可知,水滴主要撞击在叶片的前缘位置,叶梢的水滴撞击量最大,从叶梢到叶根水滴撞击量逐渐减少.这主要是因为叶片上水滴的撞击量受叶片与空气的相对速度影响,叶片上叶展位置越靠近叶梢的部位,因叶片旋转其线速度越大;在风速不变情况下,较大的线速度引起结构与水滴的撞击概率增大.撞击量较大部位产生的结冰量较多,同时导致叶片表面粗糙度增加,如图7所示.

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   叶片结冰结果图

Fig.7   Diagram of results of blade icing


为了将结冰厚度和气动性能数据耦合到一体化计算方法中,需要计算叶片19个截面结冰后的几何形状.在此选取图1(b)标注的前6个截面绘制的结冰表面与干净表面的表面轮廓线对比图,如图8所示.根据每个截面的结冰情况,对每个表面结冰的厚度进行计算,结果如图9所示.

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   风机截面翼型覆冰前后轮廓图

Fig.8   Profiles of wind turbine airfoil before and after icing


图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   结冰厚度沿展长的变化图

Fig.9   Variation of ice thickness along blade span


从图8和图9可以看出,叶根到叶梢方向,风机叶片结冰厚度呈线性增长.这主要是因为叶片结冰情况受叶片与空气的相对速度影响,叶片上越靠近叶梢的部位,因叶片旋转其线速度越大;在风速不变情况下,较大的线速度引起结构与水滴的撞击量增大,导致结冰严重.NACA64-A17翼型属于叶梢部位翼型,其结冰轮廓区域明显增大.叶梢处结冰厚度最大,结冰使弦长增加了15.67%,结冰主要集中在叶片前缘位置,尾缘下表面会产生轻微的结冰现象.叶片结冰情况与水滴撞击量分布规律相符,撞击量较大的部位造成的结冰量较多.由于叶片旋转过程中每个截面与空气来流方向形成非0° 攻角,同时翼型尾缘下表面处于迎风侧,在叶片旋转过程中会在叶片尾缘下表面处产生轻微结冰现象.

图10~12为叶片各翼型结冰前后升、阻力系数对比图.其中:CL为升力系数.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   NACA64翼型截面升力系数和阻力系数

Fig.10   Lift coefficient and drag coefficient of NACA64 airfoil


图11

新窗口打开| 下载原图ZIP| 生成PPT

图11   DU25翼型截面升力系数和阻力系数

Fig.11   Lift coefficient and drag coefficient of DU25 airfoil


图12

新窗口打开| 下载原图ZIP| 生成PPT

图12   DU40翼型截面升力系数和阻力系数

Fig.12   Lift coefficient and drag coefficient of DU40 airfoil


除了各剖面翼型结冰后的弦长、表面形状、积冰厚度外,还需将结冰后翼型的升力系数CL和阻力系数CD进行计算,同时耦合到一体化风机分析方法中.选取图1(b)所示的3个截面NACA64_A17、DU25_A17、DU40_A17绘制升力系数和阻力系数曲线图,并与文献[22]给出的干净翼型升、阻力系数进行对比.一体化风机分析方法考虑到风机叶片处于动态过程,风速过大时风机控制系统会通过对叶片进行变螺距控制改变叶片翼型攻角,从而保证风机安全且能达到最高发电效率.为此需要将各剖面翼型在-180° 到180° 攻角范围内的升力系数和阻力系数全部计算出来,升、阻力系数采用k-ε湍流模型进行计算,并根据Selig and Eggars方法和Viterna方法[25]对升、阻力系数进行扩展和修正.

在攻角90° 时翼型阻力系数达到峰值,升力系数减小至0.在叶梢部位,NACA64_A17截面结冰翼型,相比于干净翼型,0° 攻角时升力系数降低了82.02%,阻力系数增加了214%,最大升力系数降低了49.02%,最大阻力系数增加了5.35%.DU25和DU40翼型结冰后的最大升力系数也相应降低,但仅有15.1%和10.5%,阻力系数增加了4.87%和4.58%.可以看出,结冰厚度对翼型的升阻力系数有较大影响,在叶梢部位结冰最大,其造成的升力系数降低也最为明显.

3.2 一体化风机计算结果

通过CFD方法对叶片结冰进行数值计算,并将各剖面翼型结冰后的弦长、表面形状、积冰厚度以及升、阻力系数进行统计计算传递到一体化风机分析发方法中.在一体化计算过程中,为研究结冰对海上风机整体性能产生的影响,从3~24 m/s风速区间中,以步长为1 m/s选取22个风速作为计算工况,将叶片结冰前后风机运行中的多个性能参数进行对比分析.

图13和图14为风机叶尖和塔顶的运动时程曲线,包含8 m/s和21 m/s风速工况下的结构动态响应.其中,t'为时域分析时间.风机叶片受结冰影响导致叶片升力系数降低,图13(a)和图14(a)显示风机在达到额定功率前的8 m/s 风速工况下,升力系数降低导致叶尖和塔顶纵荡位移减小较为明显,干净叶片的纵荡位移约为结冰叶片的2倍;叶尖和塔顶的横荡位移变化较少,结冰后横荡位移减少.图13(b)和图14(b)显示,在 21 m/s 风速工况下,叶片结冰后叶尖和塔顶纵荡位移明显增加,叶尖横荡位移也有所增加,与8 m/s风速工况下的结构响应情况相反.

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13   叶尖偏移距离

Fig.13   Offset distance of blade tip


图14

新窗口打开| 下载原图ZIP| 生成PPT

图14   塔顶偏移距离

Fig.14   Offset distance of tower top


图15为风力机结冰前后各项参数达到稳态时的对比,将计算的3~24 m/s风速工况按叶片结冰前的风机工作状态分为3个阶段,分别为图15中区域1、2、3.区域1为风电机组启动阶段,3 m/s≤v<7.8 m/s,该区域用于设置发电机速度下限,以限制风机运行的速度范围;区域2为优化风机功率阶段,7.8 m/s≤v≤11.4 m/s(额定风速),在该区域时叶片保持最佳的叶尖速比,风机最大限度获取风能,逐步达到额定功率;区域3为变螺距控制器工作阶段,11.4 m/s<v≤24 m/s,该区域保持风机在额定功率下工作.图15(a)~15(c)叶片结冰前,转子速度随着区域2的风速线性增加,以保持恒定的叶尖速比和最佳风力转换效率.发电功率和低速轴转矩随着区域2中的风速显著增加,分别呈3次和2次增加.在风速11.4 m/s时达到额定值,高于额定功率时,发电功率通过变螺距控制调节到固定速度来保持恒定.而叶片结冰后,转子速度、发电功率和低速轴转矩随风速增加而增长的趋势变缓,额定值由于叶片结冰影响而后移,在风速17 m/s时达到额定值,之后风机才通过变距控制进行调节.叶片结冰明显降低风机的发电效率,在8 m/s风速时,风机发电效率降低了51.11%.

图15

新窗口打开| 下载原图ZIP| 生成PPT

图15   结冰前后风机稳态响应对比

Fig.15   Comparison of wind turbine responses before and after icing


图15(d)叶片结冰前叶片叶尖纵荡位移在风机达到额定功率(风速为11.4 m/s)时达到最大值,随后随风速增加,纵荡减少,这种响应特性是风机在达到额定功率时,转子所受推力达到峰值,之后受到变螺距控制器影响,转子所受推力随风速增加而降低所产生的;而叶尖横荡、塔顶纵荡和横荡位移中,该响应峰值也是可见的.塔顶横荡较小曲线不太明显.对于结冰后的叶尖和塔顶响应来说,在风速17 m/s时这些运动响应达到峰值,而受到叶片结冰的影响,叶片受到的升力降低、阻力增大,导致叶尖纵荡运动和横荡运动分别减少了32.98%和增加了25.59%,塔顶纵荡运动减少了3.98%.

结合图13~15可以看出,在风机达到额定功率前,结冰导致风机功率大幅度降低;额定功率后,风机保持5 MW功率工作,此时叶片结冰会导致风机叶片和塔筒的运动增大,导致风机产生额外振动,对风机安全产生不利影响.

4 结论

本文深入研究三维叶片旋转结冰机理以及结冰对风机整体动态响应的影响.基于NREL 5 MW导管架式海上风机一体化三维数值模型,结合计算流体力学和风机多体动力学计算方法,建立5 MW海上风机叶片结冰耦合计算系统.得到各剖面翼型结冰后的弦长、表面形状、积冰厚度、升力、阻力系数,以及发电功率、叶片受力和振动、转子的转矩和转速、风机塔筒振动等性能,研究结冰过程及结冰后对风机整体性能的影响.研究成果可为冰区一体化海上风机安全性能以及覆冰状态下整机功率损失的研究提供方法依据,具体结论如下:

(1) 风机运行过程中,叶片转动使得叶尖处线速度最大,导致该处水滴撞击量最多;三维叶片结冰厚度沿叶展方向呈现线性增长,与水滴撞击量分布情况相吻合,叶尖处结冰最为严重.结冰主要集中在叶片前缘,尾缘处有轻微结冰.

(2) 叶片结冰导致叶片前缘外凸,叶片表面粗糙度增加.结冰后翼型的升力系数明显降低,阻力系数增大,叶尖处积冰情况最为严重,导致叶尖的翼型升、阻力系数变化最大.NACA64-A17为叶梢翼型,其在 -180° 到180° 攻角范围内最大升力系数降低了49.02%,最大阻力系数增加了5.35%;0°攻角时升力系数降低了82.02%,阻力系数增加了214%.

(3) 对NREL 5 MW风机结冰前后的结构动态响应分析得出,叶片结冰对风机结构及性能的影响在风机额定功率前后表现不同.达到额定功率前,风机响应主要表现为叶片受力减少,转子转速降低,风机发电功率严重降低;达到额定功率后,虽然风机仍保持 5 MW 额定功率工作,但此时叶片结冰主要影响体现在叶尖和塔筒位移增大,风机结构产生额外振动,风机安全生产受到严重威胁.

(4) 对本文工况下叶片覆冰后的OC4风机进行一体化分析的过程中发现,结冰可以导致风机发电功率、转子转速和转矩显著降低.8 m/s风速时风机发电效率降低了51.11%,而该风机达到额定功率所需风速由原来的11.4 m/s变为17 m/s,增大了44.7%.升力降低导致叶尖和塔顶纵荡位移减少,阻力增加导致叶尖横荡位移额外增大,而塔顶横荡位移变化较小.

参考文献

LEHTOMÄKI V.

Emerging from the cold

[J]. Windpower Monthly, 2016, 32(8): 32-34.

[本文引用: 1]

HOCHART C, FORTIN G, PERRON J, et al.

Wind turbine performance under icing conditions

[J]. Wind Energy, 2008, 11(4): 319-333.

DOI:10.1002/we.258      URL     [本文引用: 1]

LEHTOMÄKI V, RISSANEN S, WADHAM-GAGNON M, et al.

Fatigue loads of iced turbines: Two case studies

[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2016, 158: 37-50.

DOI:10.1016/j.jweia.2016.09.002      URL     [本文引用: 2]

王可光, 吴辉碇, 王彩欣, 等.

渤海冰期的基本水文气象参量研究

[J]. 海洋通报, 1999, 18(2): 17-28.

[本文引用: 2]

WANG Keguang, WU Huiding, WANG Caixin, et al.

A study of basic hydrologic and meteorological parametersin the ice-covered Bohai Sea

[J]. Marine Science Bulletin, 1999, 18(2): 17-28.

[本文引用: 2]

SHIN J, BOND T.Results of an icing test on a NACA 0012

airfoil in the NASA Lewis Icing Research Tunnel

[C]// 30th Aerospace Sciences Meeting and Exhibit. Reston, Virginia: AIAA, 1992: 647.

[本文引用: 1]

HAN Y Q, PALACIOS J, SCHMITZ S.

Scaled ice accretion experiments on a rotating wind turbine blade

[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2012, 109: 55-67.

DOI:10.1016/j.jweia.2012.06.001      URL     [本文引用: 1]

胡良权, 陈进格, 沈昕, 等.

结冰对风力机载荷的影响

[J]. 上海交通大学学报, 2018, 52(8): 904-909.

[本文引用: 1]

HU Liangquan, CHEN Jinge, SHEN Xin, et al.

Load of wind turbine affected by icing

[J]. Journal of Shanghai Jiao Tong University, 2018, 52(8): 904-909.

[本文引用: 1]

SWITCHENKO D, HABASHI W, REID T, et al.

FENSAP-ICE simulation of complex wind turbine icing events, and comparison to observed performance data

[C]// 32nd ASME Wind Energy Symposium. Reston, Virginia: AIAA, 2014: 1399.

[本文引用: 1]

蒋维, 李亚冬, 李海波, 等.

水平轴风力机桨叶覆冰数值模拟

[J]. 太阳能学报, 2014, 35(1): 83-88.

[本文引用: 1]

JIANG Wei, LI Yadong, LI Haibo, et al.

Simulation of icing on horizontal-axis wind turbine blade

[J]. Acta Energiae Solaris Sinica, 2014, 35(1): 83-88.

[本文引用: 1]

郝艳捧, 刘国特, 阳林, 等.

风力机组叶片覆冰数值模拟及其气动载荷特性研究

[J]. 电工技术学报, 2015, 30(10): 292-300.

[本文引用: 1]

HAO Yanpeng, LIU Guote, YANG Lin, et al.

Study on ice numerical simulation and its power loss characteristics for the blades of wind turbine

[J]. Transactions of China Electrotechnical Society, 2015, 30(10): 292-300.

[本文引用: 1]

梁健, 舒立春, 胡琴, 等.

风力机叶片雨淞覆冰的三维数值模拟及试验研究

[J]. 中国电机工程学报, 2017, 37(15): 4430-4436.

[本文引用: 1]

LIANG Jian, SHU Lichun, HU Qin, et al.

3-D numerical simulations and experiments on glaze ice accretion of wind turbine blades

[J]. Proceedings of the CSEE, 2017, 37(15): 4430-4436.

[本文引用: 1]

刘杰, 杨倩, 吴涛, 等.

霜冰条件下风力机翼型结冰的数值计算预测

[J]. 机电一体化, 2020, 26(5): 3-11.

[本文引用: 1]

LIU Jie, YANG Qian, WU Tao, et al.

Numerical simulation prediction of icing on airfoil of wind turbine blade under rime ice conditions

[J]. Mechatronics, 2020, 26(5): 3-11.

[本文引用: 1]

WANG Q, XIAO J P, ZHANG T T, et al.

A new wind turbine icing computational model based on Free Wake Lifting Line Model and Finite Area Method

[J]. Renewable Energy, 2020, 146: 342-358.

DOI:10.1016/j.renene.2019.06.109      URL     [本文引用: 1]

WANG Q, YI X, LIU Y, et al.

Simulation and analysis of wind turbine ice accretion under yaw condition via an Improved Multi-Shot Icing Computational Model

[J]. Renewable Energy, 2020, 162: 1854-1873.

DOI:10.1016/j.renene.2020.09.107      URL     [本文引用: 1]

National Renewable Energy Laboratory. FAST v8[CP/OL]. (2016-07-21) [2021-07-14]. https://www.nrel.gov/wind/nwtc/fastv8.html.

URL     [本文引用: 1]

WIROGO S, SRIRAMBHATLA S.

An eulerian method to calculate the collection efficiency on two and three dimensional bodies

[C]// 41st Aerospace Sciences Meeting and Exhibit. Reston, Virginia: AIAA, 2003: 1073.

[本文引用: 1]

SCHILLER L V, NAUMANN Z Z.

Über die grundlegenden berechnungen bei der schwerkraftaufbereitung

[J]. Zeitschrift Des Vereines Deutscher Ingenieure, 1933, 77: 318-321.

[本文引用: 1]

MESSINGER B L.

Equilibrium temperature of an unheated icing surface as a function of air speed

[J]. Journal of the Aeronautical Sciences, 1953, 20(1): 29-42.

DOI:10.2514/8.2520      URL     [本文引用: 1]

JONKMAN J, BUHL M.

New developments for the NWTC’s FAST aeroelastic HAWT simulator

[C]// 42nd AIAA Aerospace Sciences Meeting and Exhibit. Reston, Virginia: AIAA, 2004: 504.

[本文引用: 1]

KANE T R, LEVINSON D A. Dynamics, theory and applications[M]. New York, USA: McGraw Hill, 1985.

[本文引用: 1]

POPKO W, VORPAHL F, ZUGA A, et al.

Offshore code comparison collaboration continuation (OC4), Phase 1—Results of coupled simulations of an offshore wind turbine with jacket support structure

[C]// The twenty-second international offshore and polar engineering conference. Rhodes, Greece: OnePetro, 2012: 337-346.

[本文引用: 1]

JONKMAN J, BUTTERFIELD S, MUSIAL W, et al.

Definition of a 5-MW reference wind turbine for offshore system development

[R]. Colorado: Office of Scientific and Technical Information, 2009.

[本文引用: 2]

FORTIN G, LULIANO E, MINGIONE G, et al.

CIRAAMIL ice accretion code improvements

[C]// 1st AIAA Atmospheric and Space Environments Conference. Reston, Virginia: AIAA, 2009: 3968.

[本文引用: 3]

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 and Industrial Aerodynamics, 2010, 98(4/5): 181-188.

DOI:10.1016/j.jweia.2009.10.014      URL     [本文引用: 1]

HANSEN C. AirfoilPrep: An excel workbook for generating airfoil tables for AeroDyn[EB/OL].(2004-11-1) [2021. 07. 14]. https://www.nrel.gov/wind/nwtc/airfoil-prep.html.

URL     [本文引用: 1]

/

〈 〉