风力发电机叶片的雾凇覆冰数值模拟

黎芷毓，蒋兴良，韩兴波，任晓东，王洋洋

(1. 输配电装备及系统安全与新技术国家重点实验室(重庆大学)， 重庆 400044；2. 重庆交通大学 机电与车辆工程学院， 重庆 400074)

风能作为可再生能源，越来越受到人们的重视，近年来，中国风电行业发展迅速，风力发电机装机容量逐年攀升，目前已处于世界第一[1-3]。然而，中国近年来极端气象灾害频发，在重庆、云南、湖南等华中和西南地区内陆风电场，冬季低温且潮湿，风力发电机易遭到覆冰灾害的影响[4-5]。结冰的风力机叶片不仅会使得输出功率降低，严重的会导致停机[6-7]。文献[8]显示1996年至2002年期间因为叶片覆冰导致芬兰某风电场风力发电机出现问题的时长达到8 000小时以上，约占低温故障总时长的70%。研究发现，风力发电机结冰导致年发电量(AEP)降低17%，空气动力学特性降低20%～50%[9]。

由于对风力发电机叶片进行人工覆冰试验耗费大量物力及财力，近年来，利用计算机技术对风力发电机覆冰进行预测研究成为了研究热点。利用计算机信息技术对叶片覆冰进行模拟计算，不仅可以为风力发电机翼型的防冰设计提供技术参考，对后续防冰除冰工作及设备维护也有重要的指导意义。早期芬兰技术研究中心Makkonen等[10-11]借鉴飞机覆冰数值模型开发了一套风力机覆冰代码TURBICE，该模型可计算得到不同温度下风力机的覆冰形状及质量。文献[12]对大型水平轴“NREL 5MW”风力发电机叶片的覆冰增长进行了数值计算发现，随着叶片几何参数的变化，覆冰的增长速率及冰形都有一定的改变。文献[13]对大型风力发电机叶片覆冰进行预测，并分析了覆冰后转矩及功率的损失，但其预测仍为二维模型，无法获得叶片整体覆冰增长趋势。文献[14]运用商用软件Fluent对水平轴风力发电机三维雨凇覆冰进行了数值计算及人工覆冰试验，然而该项研究仅针对小型风力发电机。文献[15]对常用于大型风力发电机的翼型覆冰进行了试验及仿真计算，但计算仍然停留在单翼型，无法了解整体叶片覆冰趋势。

风力发电机叶片的覆冰预测模型不仅要考虑其准确性，同时也应考虑时间和开销，目前研究其气流场的方法有：有限元法、有限体积法、有限差分法及边界元法。本文采用的边界元法计算风力发电机叶片覆冰流场，通过降维的方法，可将三/二维问题化为二/一维问题进行分析，具有计算量小，计算时间短的优点[16]。由于风力发电机叶片外部气流场可以近似看作二维流动，而大型风力发电机叶片的尺寸较大，进行三维叶片空气流场计算不仅耗费时间长，且需要较大的计算机存储空间，导致其三维覆冰模拟较难实现[17]。本文提出了一种切片-重构的方法，对风力发电机叶片进行切片得到数个二维截面，基于边界元法对叶片外部气流场进行计算，使用拉格朗日法跟踪水滴轨迹，通过优化选择时间步长，仿真模拟各二维风力发电机叶片切片的覆冰增长，最终通过截面覆冰重构，得到三维叶片雾凇覆冰图像。

1 风力发电机叶片气液两相流数值模型

当风力发电机叶片在运转时，其外部气流场合速度为风速与该截面所拥有的线速度的相反值之和，称之为相对风速，也即来流合速度。同时，风力发电机在运行的过程中，其外部空间存在气液(空气-水滴)两相流，随着温度的降低，空气中的过冷却水滴撞击到风力机叶片表面冻结形成覆冰。研究叶片外部空气流场是研究过冷却水滴运动轨迹的基础，因此，首先要对外部流场的速度分布进行计算，再将水滴视为离散相，根据空气外流场速度分布进一步计算每个水滴的运动轨迹，得到水滴在叶片表面的局部碰撞率。

边界层理论表明，结构物大气绕流场可分为边界层内的层流或湍流区域，以及边界层以外的非黏性势流区域[18-19]。

1.1 外部空气流场求解

如图1所示，假设P为外部流场中的任一点(a、b、c、d、e、f)，Q为边界上的点，在边界层以外的势流区，流体可以看作是理想流体，根据拉普拉斯方程：

(1)

式中：φ为速度势；vx，vy分别为沿x轴方向，y轴方向的来流速度。

图1 风力发电机叶片截面边界元划分

边界条件：风力发电机叶片边界为壁面，流体不可穿透，因此流体在边界的法向速度为0，即

(2)

(3)

根据格林公式得到如下边界积分方程，即流场势函数求解控制方程：

(4)

(5)

式中：φ*(P,Q)为拉普拉斯算子基本解，r(P,Q)为P点与Q点之间的距离。

需要确定边界上的势函数值，对于边界点P，为避免出现P与Q重合而产生奇异性，式(5)化为

(6)

式中C(P)是与P点处的边界几何形状有关的常数，对于光滑边界，C(P)=1/2。

本文采用的边界元法，将风力发电机叶片的边界分割成n个边界单元，整个边界上的积分以n个边界单元上的积分的和来表示，且将节点取为微元中点。通过常量单元离散，假定每个边界单元上的势函数值与其法向导数值在边界单元上为常数，且等于节点的值。式(6)可离散为

(7)

其中：

(8)

(9)

(10)

式中：i为P点所在微元序号，j为Q点所在微元序号，δij为克罗内克δ。

对于n个节点，得到联立的一次方程组，用矩阵形式可表示为

H·φ=G·Q

(11)

式中：H和G是n×n阶的系数矩阵，φ和Q分别是边界单元节点的扰动势函数值和扰动势函数的法向导数值的列向量。

对于非边界单元，式(6)可离散为

(12)

通过式(2)～(12)，采用高斯积分公式求解边界积分方程，即可求出流域内任意位置的气流速度。

1.2 边界层内流场求解

临界边界层动量厚度雷诺数可作为判断流体从层流发展为湍流的准则数[17]：

(13)

式中：θ为风力发电机叶片表面边界层动量厚度；θtr为层流发展为湍流的分界点的边界层动量厚度；μ为空气黏度系数；ρa为空气的密度；Ue为边界层速度，此处定义为99%来流速度。

卡门边界层动量积分关系式[20]为

(14)

式中：p为压强，δ为边界层厚度。

边界层厚度δ可由边界层动量厚度θl估算[21]：

δ=8.5θl

(15)

边界层动量厚度θl[22]的表达式为

(16)

式中：str表示层流与湍流分界点位置，θi(str)为分界点处的边界层动量厚度。

摩擦切应力τ0表达式[18-19]为

(17)

边界层以内的气流速度[23-24]为

(18)

式中：h为距离风力发电机叶片边界的垂直距离，ds为边界微元单位长度。

求解式(13)～(18)，即可求出边界层内气流速度。

1.3 过冷却水滴轨迹求解

本文计算水滴轨迹时利用拉格朗日法，为了得到过冷却水滴在风力发电机叶片表面的局部碰撞率，需要对过冷却水滴运动轨迹进行计算。假设水滴在距离障碍物较远时(设水滴从距离叶片前缘10倍弦长处发射)与空气拥有相同的速度，靠近障碍物时，气流沿障碍物表面绕流，过冷却水滴轨迹由于惯性的作用，与气流轨迹产生偏差，此时对过冷却水滴进行受力分析后可知，过冷却水滴在运动过程中主要受到重力FG、气流浮力Fb、水滴前后的压差阻力Fg、表观质量力FM、空气黏性阻力Fd的作用，其中占主要作用的是空气黏性阻力Fd，此处仅考虑空气黏性阻力的作用[25]，根据牛顿第二定律，得

(19)

相对雷诺数表示为

(20)

将式(20)代入式(19)，得

(21)

(22)

k为运动水滴惯性的斯托克斯数，计算公式为

(23)

使用差分法对式(19)求解，根据速度和加速度的一阶和二阶差分格式可得到水滴在tn时刻沿x轴方向的速度及加速度的表达式为

(24)

式中：Δt为时间步长，xn+1、xn和xn-1分别为水滴在tn+1时刻、tn时刻和tn-1时刻的位置的横坐标。同理可得水滴沿y轴方向的速度及加速度。

水滴轨迹追踪方程：

(25)

式中：va_x,n和va_y,n分别为水滴在tn时刻所在位置的气流沿x、y轴方向的速度，vd_x,n和vd_y,n分别为水滴在tn时刻沿x、y轴方向的速度，yn+1、yn和yn-1分别为水滴在tn+1时刻、tn时刻和tn-1时刻的位置的纵坐标。

通过式(25)可根据水滴初始位置和初始速度求出水滴轨迹。

1.4 过冷却水滴碰撞率求解

风力发电机叶片水滴局部碰撞率α1计算公式：

(26)

式中：dL为两相邻水滴在叶片上碰撞位置在叶片表面的距离，dY为相邻水滴纵坐标之差。

基于以上算法计算得出的叶片表面气流场和水滴运动轨迹如图2所示。

2 风力发电机叶片雾凇覆冰模型计算

2.1 叶片表面水滴冻结特性

在过冷却水滴碰撞到风力机叶片表面的过程中，忽略水滴飞溅的影响，因此捕获率α2≈1。覆冰是一种热交换的过程，在叶片表面水滴冻结过程中，可以通过传质传热平衡方程计算叶片表面的冻结率，传质传热方程表示为

(27)

其中：mim为撞击到叶片表面水滴的质量；min和mout为流入和流出的水质量；meva为蒸发水质量，mice为冻结水质量；Qf为空气摩擦叶片表面产生的热量；Qim为碰撞到叶片的过冷却水滴带来的动能转化的热量；Qi为水滴冻结成冰后释放的热量；Qin和Qout为流入水带来的能量和流出水带走的能量，设风力发电机叶片表面水膜温度Tw=0 ℃，则水膜的流动在微元之间不会产生能量交换，Qin、Qout均为0；Qeva为蒸发升华热损失；Qc为空气对流换热损失的热量；Qw为将叶片表面捕获的过冷却水滴加热到水膜温度所需要的热量；Qk为冰面传导热损失；Qr为长波辐射热损失。

忽略式(27)中的较小的热量项摩擦热Qf和传导热损失Qk，将各热量控制项表达式代入式(27)后可得到冻结率的计算方法：

(28)

2.2 叶片冰形计算

风力发电机叶片表面微元控制体在覆冰速率dm/dt计算方法如下：

dm/dt=V·α1·α2·α3·w·ds

(29)

式中：V为水滴碰撞在叶片上的速度；w为空气液态水含量，kg/m3；ds为控制体微元的长度。

在覆冰不脱落的情况下，覆冰沿叶片法向生长，每个控制体微元的覆冰厚度增长速率为

(30)

式中：dt为覆冰时间步长；ρice为冰密度，其经验公式[26]为

(31)

其中R为Macklin冰密度参数，R=-URd/TS，TS为叶片表面温度。

本文基于某型300 kW风力发电机叶片(叶片全长14.5 m，额定风速13 m/s)进行建模，将叶片等分为数个翼型截面进行分析，分别计算获得气-液两相流轨迹，并根据不同的截面选取合适的覆冰时间步长dt，计算覆冰厚度d，进而得出覆冰冰形，迭代直至达到总覆冰时长Tice=30 min，计算流程见图3。

图3 三维翼型覆冰仿真流程图

3 模型计算结果及分析

3.1 算法验证

水滴局部碰撞系数是覆冰计算中较为关键的参数，为了验证本文提出的计算空气流场和水滴轨迹最终获得碰撞率的方法的有效性，采用文献[27]中的相关参数对二维NACA0012翼型进行仿真。其中，来流合速度为44.39 m/s，液态水含量为0.78 g/m3，温度为-7.65 ℃，水滴中值直径为20 μm。并将试验结果及本文计算出的碰撞率计算结果进行比较。

试验结果及计算结果如图4所示，s表示碰撞点距离驻点的位置，翼型的水滴局部碰撞率呈现中间大两头小的现象。试验中s>0即叶片下表面碰撞范围略大于上表面而碰撞率略小于上表面，这是因为水滴在运动的过程中受重力的作用会稍向下偏移，但计算时由于忽略了重力的影响导致碰撞范围和碰撞率是上下对称的。通过对比可知，本文使用边界元法数值计算出的碰撞率与实验结果吻合良好。计算值的最大局部碰撞率略大于试验值，其他部分基本吻合。但计算的碰撞范围小于试验值，产生误差的可能原因有二：一是因为本文采用的离散方法求解流场积分的方程和求解水滴轨迹所取的时间步长导致的计算误差；二是因为文献中的覆冰试验不易控制，测量精度有一定的分散性，此外，在实际的覆冰试验中，水滴直径的取值分布在一定范围内，其中大尺寸水滴会导致试验中水滴碰撞范围增大，而本文采用单一水滴直径进行计算，水滴碰撞范围偏小，但从整体上看，本文数值计算的水滴局部碰撞率与试验值非常接近，表明使用该方法计算风力发电机叶片外部空气流场、水滴运动轨迹和水滴碰撞是合理的。

图4 NACA0012翼型水滴局部碰撞率对比

3.2 风力发电机叶片的覆冰增长特性

为了分析沿叶展方向的整体叶片覆冰增长特性和风速变化对覆冰形态及覆冰量的影响，本文选取了3个风速进行了冰形计算，其中，风速U分别为3、6、9 m/s，对应的风轮转速分别为27、32、38 r/min，液态水含量为0.3 g/m3，温度为-15 ℃，水滴中值直径为20 μm。图5～7为计算得到的三维覆冰叶片。图5显示，覆冰区域基本主要集中在迎风面，P—P截面(叶尖部)处背风面尚存在一定的覆冰，而J—J截面(叶中部)和E—E截面(叶根部)处，冰基本只存在于叶片的迎风面，在图5(b)中几乎看不到冰的存在，这是因为沿着叶展方向翼型的不断变化和攻角的逐渐减小。同时，仿真结果显示，覆冰厚度沿着叶展的方向逐渐增大，距离叶根4 m处(E—E截面)的覆冰冰厚仅为0.15 mm，而距离叶根14.44 m处(P—P截面)的叶尖部位，冰厚达到22 mm。过冷却水滴在空气中受到惯性和阻力的影响向前运动，如果气流阻力占主要作用，则过冷却水滴将沿着气流的方向运动而绕过叶片，沿着叶展方向的翼型弦长和厚度逐渐减小，使阻力的作用逐渐减小而惯性的影响逐渐增大，同时由于沿着叶展方相对速度的增加，更加大了惯性对水滴运动的影响，最终导致了越靠近于叶尖部位截面的覆冰越厚。

图5 覆冰后三维叶片成像(U=9 m/s)

(a)叶尖冰形 (b)叶中冰形 (c)叶根冰形

图7 不同风速下叶尖成像

3.3 风速对风力发电机覆冰的影响分析

图8显示了在不同风速下沿叶展方向不同部位的覆冰情况。图8表明，风速的增大会使得水滴局部碰撞率增加，且单位时间内将有更多的水滴碰撞在截面上，导致覆冰厚度的增加，同时也会改变覆冰的面积和区域，在风速为3 m/s时，叶尖部位的覆冰区域明显整体小于风速为6 m/s和9 m/s时，这是由于随着风速的增大，叶片转速也将逐渐增大，气流合速度随之增大，过冷却水滴的惯性对水滴运动的影响逐渐增大，在风速较小时本会绕过叶片运动的水滴在风速增大后可能会碰撞在叶片上，导致碰撞面积的增大。但随着风速的增加，背风面覆冰区域逐渐减小，而迎风面覆冰区域逐渐增加，这是由风力发电机叶片的空气动力学特性决定的，随着风速和转速的增加，气流的攻角也会逐渐增加，使水滴碰撞区域从背风面向迎风面移动。同时，在叶中部和叶根部也出现了这种趋势，随着覆冰面积增大的同时，覆冰区域从背风面向迎风面移动，且随着截面位置向轮毂靠近，覆冰区域将只存在于叶片的迎风面，但与叶尖相比，叶根附近受风速的影响明显降低，在距离轮毂8 m(I—I截面)处，尚能看到明显的厚冰，而更靠近叶根的两个截面(F—F、C—C截面)几乎已经看不到冰。另外，风速的变化对冰形几乎没有影响，同一截面下，3种风速形成的覆冰冰形基本相同。为了分析风速对冰厚及覆冰增长速率的影响，设置沿截面法向的最大冰厚度为该处截面的冰厚。

(a) O—O截面(叶尖) (b)L—L截面(叶尖-叶中) (c)I—I截面(叶中)

(d)F—F截面(叶中-叶根) (e)C—C截面(叶根)

图9显示了沿叶展方向的不同风速的冰厚变化，在3种风速下，覆冰厚度沿着叶展方向逐渐增大，在叶尖处达到最大厚度，从叶根部到叶片中部(距离轮毂0～6 m)，3种风速下的冰厚度相差不大，即冰厚受风速的影响较小，且冰厚增长较缓，此时的过冷却水滴运动主要受阻力作用的影响，同时相对气流合速度差别不大，且该部分并没有明显的覆冰现象，从叶片中部到叶尖部分(距离轮毂6～14.6 m)，水滴惯性逐渐占据主要作用，冰厚逐渐增大。随着风速的增大，冰厚沿叶展的增长也随之加快，冰厚差增大，一直到叶尖位置，当风速从3 m/s增加到9 m/s时，冰厚从16 mm增长到22 mm，冰厚增长速率从0.53 mm/min增长到0.73 mm/min。

图9 不同风速下冰厚随叶展的变化

综合图8、9，受覆冰影响最大的区域是从叶片中部到叶片尖部，此趋势与文献[12]对兆瓦级大型风机的覆冰仿真中计算出的冰厚增长趋势基本一致。

4 结 论

基于切片-重构的思想，将风力发电机叶片沿叶展切为数个截面，对于每个截面，采用边界元法对风力发电机叶片周围大气流场进行计算，并用拉格朗日法对过冷却水滴的轨迹进行追踪。仿真得到每个截面的冰形，再将所有截面及冰形进行复原重构，最终建立了风力发电机叶片三维雾凇计算模型，同时对边界元算法用于计算大气流场的有效性进行了检验，分析了风速对叶片雾凇覆冰的影响规律。通过对比和分析可得出如下结论：

1)在风力发电机叶片雾凇覆冰模拟中，碰撞率的计算至关重要，本文方法所得计算值和试验值十分接近，证明所提出的计算方法是合理且有效的。

2)对某型300 kW风电设备叶片雾凇覆冰进行了仿真模拟，结果表明，雾凇覆冰冰形呈现流线型。覆冰主要聚积在叶片的前缘和迎风面靠近前缘位置，沿着叶展的方向，覆冰逐渐增厚，覆冰区域占比增大。

3)风速对雾凇覆冰的影响主要在冰厚和覆冰区域，对冰形的变化无明显影响，且风速的变化对覆冰的影响在叶尖部位更加显著，在靠近叶根的位置，不同风速产生的冰厚都较薄，从叶中部到叶根，3种风速下冰厚的差别越来越大，直到叶尖冰厚度相差6 mm，且随着风速的增加，覆冰区域变大的同时，覆冰区域从背风面向迎风面移动。

4)叶片的几何参数(弦长，厚度和扭转角)的变化，都会对覆冰的增长有一定的影响，这意味着通过优化叶片的几何参数，可以一定程度上降低其覆冰速率，从而减少寒冷潮湿地区安装在风电设备上的防冰和除冰装置的数量。