垂直轴风力机在不同运动状态下的电加热融冰策略研究

张汉文，罗日成，刘志勇

（长沙理工大学电气与信息工程学院，湖南长沙 410114）

0引言

近年来我国中南、西南等地新建了多个风电基地，但这些地区春冬季节的湿冷气候使得风力机面临着严峻的覆冰问题[1]，[2]。覆冰后，叶片气动性能下降，功率输出降低，荷载增大造成主轴磨损加剧，以及其它安全事故[3]。目前，风力机防除冰技术主要有涂料、热空气、电加热、电磁脉冲和超声波振动等技术，其中电加热被认为是最有效和最经济的技术[4]。国内外相关学者在风力机电加热防除冰技术的研究上已取得了一定成果。文献[5]通过ANSYS软件对风力机叶片表面电热元件的布局形式、几何形状进行了研究，结果表明，圆形加热片除冰速度更快，交错布局可以降低叶片表面的热应力。文献[6]对风力机叶片表面的电阻丝布置方式进行了试验研究。文献[7]对风力机电加热系统的控制理论进行了研究，通过仿真分析了周期性加热对防除冰的影响。文献[8]提出了一种风力机叶片循环控制加热的除冰方法，并利用数值仿真进行了验证。

垂直轴风力机具有低成本、占地空间小、无需偏航装置、易于维护和安全性高等优点[9]。随着垂直轴风力发电技术的不断进步，已有一些技术较为完善的兆瓦级垂直轴风力机推向风电市场。现有的风力机电加热防除冰技术主要是基于水平轴风力机的结构特点进行研发的，而垂直轴风力机的结构存在较大不同，在气流中进行电加热融冰时，两者的对流散热存在较大不同。因此，本文通过分析垂直轴风力机的空气流场和叶片表面冰层的融冰过程，依据叶片在不同攻角下的融冰速度差异和叶片两侧的融冰速度差异，提出了垂直轴风力机在不同运动状态下的电加热融冰策略。该策略可以加快风力机的融冰速度，减少电能损失，具有一定的工程应用价值。

1 数理模型及控制方程

在风力机电加热融冰过程中，电热单元产生的焦耳热逐渐向四周传递，附近冰层的温度缓慢升高。一段时间后，冰层表面与周围气流间产生了温度差，冰层表面出现对流散热和辐射散热。此外，当靠近叶片表面的冰层的温度上升到0℃时，接触面的冰层将逐渐融化，并产生空气间隙[10]。当一定量的内部冰层被融化后，剩余冰层与叶片表面间的附着力极大降低，剩余冰层将极易在重力、离心力等作用下发生脱落。

基于上述物理过程，本文建立了风力机的电加热融冰数值模型，分别对空气流场、热量传递过程以及相变过程进行计算。

1.1 空气流场分布

不考虑外界气流的波动以及散热项对叶片周围气流流动的影响，本文假定风力机周围的气流场为稳态分布。空气流场的控制方程选用低速粘流的雷诺平均N-S方程，湍流模型选用对边界层湍流和自由剪切湍流均有较好模拟效果的SST k-ω湍流模型，该模型能够较好地还原风力机周围的空气流场分布[11]。输运变量为湍流动能k和比耗散率ω，其输送方程为

式中：u为气流速度场；ρ为空气密度；μ为空气动力粘度；t为时间；Pk为k的有效生成率；Pω为ω的有效生成率；β，β*，σk，σω和σω2均为经验系数；F1为混合函数。

1.2 传热计算

在融冰过程中，电热单元产生的焦耳热主要有3个去向：①叶片、冰层吸收热量温度升高；②冰层表面与周围环境间的对流散热和辐射散热；③部分冰层吸收热量相变融化。

1.2.1温度场分布

叶片、冰层的传热形式属于固体域传热，周围气流的传热形式属于流体域传热，这两种传热形式存在一定差异。

周围气流的温度场可由流体域内的时变传热微分方程描述。

式中：Tf为流体域的温度场；Vf为流体域的速度场；Cp，λ分别为空气流的比热容和热导系数。

在不考虑部分冰层相变融化的情形下，叶片、冰层的温度场可由固体域内的时变传热微分方程描述。

式中：TS为固体域的温度场；qin为固体域内的净热源。

当接触面的冰层融化成液态水时，该区域的温度场表达式为

式中：▽mice为融化冰层的质量；Lm为冰的融化潜热，Lm=330 kJ/kg；qnet为该融化区域内的净热源。

1.2.2对流散热项与辐射散热项

当冰层表面的温度为TW，环境的温度为Ta时，冰层表面的对流散热项Δqf和辐射散热项Δqr分别为

式中：TW，Ta分别为冰层表面和环境的温度；hc为冰层表面的对流换热系数，其值与局部空气流场有关；εi为冰层表面辐射率，取0.95；σ为斯特藩常量，取5.67×10-8W/（m2·K4）；S1为冰层表面与周围环境的接触面积。

1.2.3冰水相变过程

实际的冰水相变过程中，冰水混合物的温度会一直保持在0℃，而数值计算模型会假设冰水相变过程发生在0℃附近一个较小的温度区间，否则无法求解。本文使用显热容法来描述冰水相变过程（图1）。

图1 冰水相变数值模型Fig.1 Ice/water phase change numerical model

图1中：体积分数θ1，θ2分别为冰和液态水的含量；冰相和水相之间的相变温度Tpc，1→2设为0℃；转变间隔ΔT1→2设为0.001℃。相变过程中潜热L1→2设为330 kJ/kg，相变过程中涉及的一些等效物理参数计算如下：

式中：ρeq，Cp_eq，λeq分别为等效密度、等效热容和等效导热系数。

2 风力机电加热融冰模型

2.1 几何模型

文献[12]对NACA0018翼型垂直轴风力机进行了覆冰研究，得到了多种覆冰工况下风力机叶片表面冰层的形貌和增长规律。这些不同覆冰工况下风力机叶片表面冰层的形貌均较为相似，覆冰均主要集中在叶片前缘区域。其中一种覆冰工况下的叶片表面冰层如图2所示。

图2 叶片表面冰层的形貌和增长规律Fig.2 Morphology and growth law of ice on blade surface

依据上述垂直轴风力机叶片表面冰层的形貌特点，通过多物理场耦合仿真软件COMSOL建模，本文的仿真模型如图3所示。

图3 风力机电加热融冰几何模型Fig.3 VAWT's electric heating ice-melting geometric model

覆冰叶片由叶片自身、电热元件和冰层（分为外冰层和内冰层，内冰层又分为内冰层_outside和内冰层_inside）构成。整个计算域尺寸为10 m×8 m，风力机直径为3.6 m，叶片弦长为0.5 m，叶片翼型为NACA0018。模型存在如下假设条件：传热过程不影响叶片附近的空气流场；忽略冰水相变过程中两者密度的变化。

2.2 参数设置

融冰模型为二维模型，面外厚度为0.2 m，所有叶片的热功率均设为60 W，对应的热功率密度为580 W/m2。环境参数：来流风速为5 m/s，气流温度为-5℃，各物体的初始温度均为-5℃，液态水含量（LWC）为0 g/m3，压力出口为0 MPa，内冰层设为相变材料。模型中各材料的物理参数见表1。其中外冰层最厚处为14 mm，最薄处为2 mm。模型网格划分如图4所示。

表1 各材料的物理参数Table 1 Physical parameters of each material

图4 网格划分Fig.4 Mesh division of model

3 空气流场分布

通过“冻结转子”求解器得到风力机在λ为0和0.2两种叶尖速比下的空气流场（图5）。由图5可知，当风力机处在不同方位时，叶片附近区域的空气流场存在差异，风力机处于静止状态（λ=0）和处于低转速状态（λ=0.2）的空气流场存在一定程度的相似。

图5 不同方位下风力机周围的气流流场分布Fig.5 In different directions，distribution of air flow field around the VAWT

以V外，V内分别计算叶片表面两侧的平均风速，在λ为0和0.2两种叶尖速比下，叶片两侧的V外，V内变化曲线如图6所示。

图6 叶片两侧风速的变化规律Fig.6 Variation law of wind speed on both sides of blade

从图6可知：当叶片位于迎风侧时，叶片外表面的空气流场流速高于叶片内表面；当叶片位于背风侧时，叶片外表面的空气流场流速多数情况下低于叶片内表面。

4 风力机处于静止状态下的融冰策略研究

4.1 温升融化过程分析

加热640 s后，不同方位和不同α下风力机周围的温度场分别如图7，8所示。为了方便比较，所有温度场云图均采用相同标尺。

图7 不同方位下风力机周围的温度场分布Fig.7 In different directions，distribution of temperature field around the VAWT

图8 不同α下叶片周围的温度场分布Fig.8 Temperature distribution around blades at different angles of attack

从图7，8可知：不同α下叶片周围的温度场存在较大差异，叶片前缘和叶片尾缘存在一条明显的温升带；当α分别为240°和300°时，叶片的背风侧均存在明显的升温区域，这是由于在加热过程中，通过对流散热和辐射散热从叶片表面所释放出来的热量在叶片附近这些流速较低的区域积聚较多。

依据文献[10]，本文假设当内冰层融化90%时，认为叶片表面剩余冰层可在重力、离心力等作用下脱落。从开始加热至内冰层即将开始融化所经历时长为内冰层的升温期，从内冰层开始融化至融化90%所经历时长为内冰层的90%融化期。不同方位下风力机叶片内冰层的受热过程如图9所示。不同α下风力机叶片内冰层的受热过程如图10所示。

图9 不同方位下风力机的融冰情况Fig.9 Ice-melting of the VAWT in different directions

图10 在不同α下叶片的融冰情况Fig.10 Ice-melting of blade at different angles of attack

从图9可知，由于风力机的融冰时间取决于该风力机中所需加热时间最长的那个叶片，则C方位下风力机的融冰速度最快，各叶片间的融冰速度差异也小，其次为B方位和D方位，A方位下风力机的融冰速度最慢。从图10可知：α靠近0°的叶片的融冰速度较慢，当叶片位于迎风侧时，叶片内表面的融冰速度快于叶片外表面；而当叶片位于背风侧时，叶片外表面的融冰速度多数情况下快于叶片内表面，这是由于叶片两侧的对流散热情况不同所导致的。

4.2 融冰策略分析

当对在气流中处于静止状态下的垂直轴风力机进行电加热融冰时，可采用下述几种措施：①安装特定机械装置，利用该装置将风力机调整至最佳融冰方位，C方位下风力机融冰速度最快，A方位下风力机融冰速度最慢；②α靠近0°的叶片的融冰速度较慢，通过控制系统来适当增大这些α下的叶片的电热功率；③当叶片处于迎风侧或背风侧时，叶片表面两侧的融冰速度存在差异，并且这种差异会随着风速的增大、环境温度的下降而变大，通过控制系统来适当调节叶片表面两侧的电热功率来缩短叶片两侧的融冰差异。

5 风力机处于低速旋转状态下的融冰策略研究

5.1 融冰策略分析

当风力机处于旋转状态时，叶片两侧会周期性交替出现高速、低速气流场，将持续影响冰层表面的散热项。本文提出一种电热功率动态调节方式，通过跟踪风力机中各叶片α的变化来动态调节叶片表面两侧的电热功率：当叶片处于迎风侧时，减少叶片外表面的对流散热损耗；当叶片处于背风侧时，减少叶片内表面的对流散热损耗。

5.2 建立等效模型

考虑到在旋转状态下叶片两侧的空气流场变化情况，图11为等效简化模型。利用两个风速函数分别控制入口A和入口B的气流场，以5 m/s的风速下风力机叶尖速比为0.2为例（图6），高速气流设为5 m/s，低速气流设为0.2 m/s，周期为T，前T/2周期内叶片表面A侧为高速气流，B侧为低速气流，后T/2周期内A侧为低速气流，B侧为高速气流。

图11 等效模型Fig.11 The equivalent model

5.3 仿真分析

叶片两侧的电热功率分别为PA和PB，存在两种功率调节方式：①当风速高时，电热功率减弱，当风速低时，电热功率加强；②当风速高时，电热功率加强，当风速低时，电热功率减弱。以第一种功率调节方式融冰时，在前T/2周期内，A侧为低电热功率，B侧为高电热功率，在后T/2周期内，A侧为高电热功率，B侧为低电热功率；以第二种功率调节方式融冰时，A，B两侧电热功率与上述情况相反。加入恒功率组（PA，PB均为75 W）作为对照组，第一种功率调节方式和第二种功率调节方式的PA，PB分别在75 W的基础上进行±25%的调整，且任意时刻PA+PB=150 W。

研究风力机在λ分别为0.2，0.1，0.05，0.025，0.012 5和0.006 25下的融冰过程，环境参数：高速气流为5 m/s，低速气流为0.2 m/s，气流温度为-8℃，各物体的初始温度均设为-8℃，LWC为0 g/m3，压力出口为0 MPa。

图12为3种电热功率调节方式下叶片的融冰情况。从图12可知：当λ分别为0.2和0.1时，第二种功率调节方式的融冰速度均最快，第一种功率调节方式的融冰速度均明显慢于恒功率方式；当λ为0.05时，第二种功率调节方式的融冰速度依旧最快，第一种功率调节方式快于恒功率方式；当λ分别为0.025，0.012 5和0.006 25时，3种功率调节方式的融冰速度均较为相近，恒功率方式具有微弱优势。

图12 3种电热功率调节方式下叶片的融冰情况Fig.12 Ice-melting of blade under three electrothermal power regulation modes

6 结论

①风力机处于静止状态（λ=0）和处于低转速状态（λ=0.2）的空气流场存在一定程度的相似。

②对在气流中处于静止状态的风力机进行电加热融冰时，叶片在不同α下的融冰速度存在明显差异，其中α靠近0°时的叶片融冰速度明显较慢；当叶片位于迎风侧或背风侧时，叶片两侧的融冰速度存在差异，当叶片位于迎风侧时，叶片内表面的融冰速度快于叶片外表面，当叶片位于背风侧时，叶片外表面快于叶片内表面。因此，安装机械装置对风力机进行方位调整，通过控制系统来适当调节特定攻角的叶片的电热功率和叶片两侧的电热功率分配可以提升融冰速度。

③基于垂直轴风力机在旋转状态下的叶片表面两侧的空气流场差异提出了一种电热功率动态调节方式，在该策略中，当叶片位于迎风侧时，减少叶片外表面的对流散热损耗；当叶片位于背风侧时，减少叶片内表面的对流散热损耗。利用所建立的等效模型仿真发现，当λ为0.05～0.2时，所提出的电热功率动态调节方式的融冰速度快于恒功率方式。