S809风力发电机叶片翼型雾凇覆冰程度及气动性能模拟*

2021-11-17李震万涛郝宇超于跃

内蒙古科技大学学报 2021年2期

李震，万涛，郝宇超，于跃

(内蒙古科技大学机械工程学院，内蒙古包头 014010)

李翰涛等结果表明雨凇覆冰大幅降低二维翼型的气动特性[1]；邓杰等研究了覆冰以及除冰后叶片表面气动特性，翼型前缘覆冰对气动性能影响显著，且相对厚度越小的翼型影响越大[2]；李岩等对NACA0018翼型进行风洞试验，试验分析结果表明攻角不同会导致冰形不同，来流方向的投影面积越大，结冰量与速率越大[3]；JAIWON S对NACA0012翼型实验和冰形计算表明随着冰面形状转变为釉冰，叶片的阻力系数在接近冰点时急剧增加，并且发现冰层表面的等效颗粒粗糙高度主要受液态水含量(LWC)、环境温度(TS)、水滴中值直径(MVD)等因素影响，受来流速度的影响较小；JANUSZ S显示纯雾气结冰发生在比釉面结冰更低的温度下，并导致沉积冰表面的粗糙度更高[4].

覆冰类型主要分为雨凇以及雾凇[5]，主要基于ANSYS FENSAP-ICE软件研究S809翼型硬雾凇覆冰后的覆冰程度以及气动特性，主要考虑环境温度(TS)对覆冰厚度的影响和攻角对气动特性的影响，并对覆冰模拟数据进行一定程度拟合并且进行数据验证.

1 研究方法

1.1 翼型简介

S809翼型是NERL S系列风力发电机专用翼型之一，相比NACA系列翼型可使年能源利用率得到显著提升，失速调节风力机提高23%～35%，变桨距风力机提高8%～20%，变转速风力机提高8%～10%[6].

1.2 数值计算模型

利用FENSAP-ICE软件依次进行空气绕流、水收集系数以及后续的覆冰模拟计算，风力机翼型流场控制方程采用时均Navier-Stokes方程.

选取Spalart-Allmaras湍流模型，S-A模型是1个一阶模型，于二阶模型而言，具有计算量不大、稳定性良好以及计算精度较高的特点，在航空动力学计算领域中广泛应用[7]，并且在参考文献[7]中，此模型的适用性已经得到验证.

1.3 网格划分模型以及边界条件

模拟地理环境参数是基于我国内蒙某风场进行设置的，拟定风力机功率为1.5 MW，根据风力发电机设计手册[8]计算公式，计算得叶片直径为70 m，轮毂直径设定为2 m，单叶片长度为34 m，海拔高度为2 170 m，风速设定为15 m/s.根据参考文献[9]，[10]，风力发电故障高发时间段为每年的12 m、次年1 m，于中国气象网历年数据处，查出此处每年12 m、次年1 m风场得相对湿度以及平均气温，根据饱和空气湿度对照表查得出空气液态水含量为0.8 g/m3.根据参考文献[6]，[11]，[12]，雨凇(glaze)、硬雾凇(hard rime)、软雾凇(soft rime)形成条件可用式(1)～(3)近似表达，取液滴中值直径15 mm，计算域采用(c为弦长)，通过计算得出在攻角为6°时，弦长为2.60 m[13,8]，将左侧圆弧以及上下边界设置为速度入口，右侧边界设置为压力出口，壁面设置为无滑移壁面.

网格划分四边形共计60 643个，节点数共计61 137个，经检测网格质量符合模拟需求，满足网格无关性要求，计算工况如表1所示.

(1)

(2)

(3)

式中，T为温度，℃；V为来流风速，m/s.

1.4 结冰粗糙度

结冰粗糙度主要依据美国国家宇航局结冰经验公式得到：

(4)

(5)

0.245 7LWC+1.257 1(LWC)2.

(6)

(7)

式中：基准值(ks/c)base=0.001 177，ks：表面粗糙度，mm；MVD：液态水含量，g/m3；TS：环境温度，℃；c：翼型弦长，m.

粗糙度的设定主要是为了计算表面结冰的流动解，同样会增加热通量以增加结冰速率和厚度.设置叶片表面粗糙度为0.5 mm，此数值适用范围广泛，采用此默认数值并且设为初始表面粗糙度.利用式(4)经过计算在环境温度为-15 ℃工况为表1时，表面粗糙度为0.68 mm，此数值设定为-15 ℃雾凇结冰状态下的表面粗糙度，以正交对比进行后续模拟.

表1 工况

2 结果分析与讨论

2.1 流场分析

首先对S809翼型进行流场计算，计算温度从-2 ℃开始由高到低，每2摄氏度计算一次，直至计算到-20 ℃，流场计算结果、压强数据图如图1，从翼型表面压力分布可以看出随着温度的降低，上翼型表面压力随温度一同出现梯度较少状况，下翼型表面的压力同样呈现此变化.由于来流空气视为不可压缩气体，根据查理定律：一定质量的空气在体积不变的情况下，压强与热力学温度成正比，图示结果与此结论相符，由此判定模拟结果是无误的.

图1 流场压强图

2.2 水滴场分析

在空气流场的基础上进行水滴场的计算，计算结果如图2，为水收集系数数据分布图.由水收集系数数据全局图可以看出翼型下表面水滴主要撞击在翼型前缘与弦长的2/5之间，上表面水滴主要撞击在翼型前缘与弦长的3/20之间.

由图2可以看出，上翼型下表面水滴收集量从-2 ℃开始至-20 ℃随着温度的降低呈现阶梯状增长.另外可以明显地看出水滴撞击主要集中在翼型的下表面，迎风区撞击范围远大于背风区[14]，上翼型上表面水收集系数在温度-8 ℃～-2 ℃和-20 ℃～-10 ℃区间内，均表现均匀上升，但是后者水滴收集范围更广，水滴收集量更多.

图2 水收集系数分布图

2.3 覆冰分析

最后在水滴场基础上进行覆冰模拟，由图3可知，模拟覆冰结果中冰形主要集中在翼型的前缘以及下表面.随着温度的降低，覆冰厚度呈现阶梯式均匀增长，造成此现象的原因是在结冰交界处，冰点温度保持零摄氏度不变，而外界温度的降低使得交界处的水分无法放热，因此只能通过结冰来达到温度一致，由此导致交界处的水滴结冰使冰的厚度增加.

李声茂定义了无量纲常数冰厚比：最大覆冰厚度与弦长的比值以此来表征覆冰程度[15].由于覆冰模拟是在笛卡尔坐标系下进行的，而覆冰的主要区域也是翼型的前缘以及迎风区，仅仅使用最大覆冰厚度与现场的比值虽然可以一定程度上表达覆冰程度，但是不够准确，翼型的覆冰还要考虑冰的密度以及体积对机翼造成的相关影响.

文章采用坐标转换法重新表达覆冰程度以更准确表征覆冰程度并进行量化表征.将在笛卡尔坐标系中的覆冰翼型转换到极坐标系中，由此再将覆冰翼型角度坐标与初始翼型角度坐标相减得到覆冰区域的Δθ坐标，计算覆冰的角度范围.接着再将覆冰翼型长度与初始翼型长度相减，得到各个角度上的覆冰增量Δρ，将某一温度下的覆冰增长量累加再除以此温度下的覆冰角度范围得到此温度下的覆冰程度(覆冰角度范围可从极坐标轴读出).按照此方法进行分析得到如图3结果，并经过计算得到表2数据.

图3 覆冰坐标转换图

表2 不同温度覆冰总厚度及均值

在表1工况条件下，计算得到的覆冰增长量数据表，从数据表不难看出覆冰厚度-2 ℃～-20 ℃基本呈现规律增长，每2个温度间隔冰厚增长量呈现由快之慢的现象.对数据表2用origin软件进行公式拟合，origin对多项式拟合是否准确判定，依据是相关系数R平方值(回归方程整体拟合度)是否接近1，越接近1则多项式拟合越准确，一般认定R平方大于0.99即为准确，此处计算R2>0.999.考虑到随着外界温度的降低，结冰类型改变[16]，冰的表面粗糙度不断增大，从明冰转变为霜冰[17]，最终导致结冰平均厚度骤增，对比拟合的一次、二次、三次多项式，此处选择二次多项式(8)并进行简化为最终结果式(9)：

y=(0.058±0.010 82)+(-0.107 55±0.002 26)x+(-0.001 14±1.001 16E-4)x2

y=0.058-0.107 55x-0.001 14x2.

(9)

式中：x为自变量温度，℃；y为因变量冰厚平均值，mm.

为了验证表达式的准确与否，再次计算温度-9，-13，-19和-22 ℃条件下的覆冰厚度，经过模拟以及计算，得出如表4结果，由平均绝对百分比误差MAPE(Mean Absolute Percentage Error)公式(10)，对以上数据进行预测度误差分析，计算数值越接近0模型越准确：

(10)

经过计算得出MAPE=0.457 7，因此证明以上拟合结果可行.

2.4 气动特性分析

通过模拟计算得到覆冰外形，计算在初始翼型和-15 ℃覆冰外形进行静态以及此两种状态下叶尖速比为3的气动特性，攻角范围为-10°～20°，计算结果如图4～7所示.由图4可知，在初始翼型-10°～4°攻角范围内，叶尖速比3状态下翼型气动特性优于静止状态下气动特性，在其余区间内则是静止状态气动特性优于转动状态；由图5可知，-15 ℃状态下结冰翼型3°～20°攻角范围内静态气动特性优于旋转状态，并且在超过20°攻角时静止状态气动特性降幅超过旋转状态；由图6可知，静止状态下结冰翼型在1°～20°范围内气动性能优于初始翼型，升阻比趋势表明初始翼型气动特性在超过20°攻角时将优于结冰翼型；由图7可知，-10°～5°攻角区间内雾凇状态下气动特性优于初始状态，另外5°～20°攻角范围内2种状态下气动特性基本一致，由升力系数以及阻力系数涨幅趋势可得到：在大于20°的攻角范围结冰初始翼型气动特性将优于雾凇状态.