0 引 言

随着经济社会的发展，能源问题越来越引起人们的关注。开发与利用风能日益受到人们的重视。垂直轴风力机因为具有无需对风装置，结构简单、容易制造、噪声小、机组整体稳定性好等优点，近年来，广泛得到人们的关注^[1]。当直线翼垂直轴风力机组运行在严寒、湿度大等恶劣气候条件下时，叶片表面经常出现覆冰现象，而覆冰后的叶片气动性能明显降低，叶片的升力减小，阻力增大，风能利用系数降低^[2-4]。因此，对风力机叶片进行防除冰技术研究是非常必要的。

目前，对于风力机桨叶防除冰主要有以下三大类：机械除冰方案、被动式除冰方案、主动式除冰方案^[5]。机械除冰方案主要有人工手动除冰、膨胀管除冰两种，其中手动除冰效率极低，且操作危险系数高，而加装膨胀管会影响叶片的气动性能，且噪声大，难以维护；被动式除冰主要有桨叶表面防结冰涂装及叶片表面喷洒化学除冰药剂等，该方案除冰效率稍优于手动除冰，但维护成本较高；主动式除冰包括电加热法、内部通热风法等，此方案优于前两种，但成本相对较高，且系统较复杂。尽管风电行业发展前景广阔，但对风力机叶片结冰缺乏有效的解决方案，一定程度上阻碍了风电的发展。

本文提出一种新的防除冰方法—气膜加热法。在风力机叶片前缘开设射流孔，孔射流热气在主流的压力和摩擦力下向下游弯曲，附着在壁面一定区域上，使叶片前端与冷气流隔离，以达到防、除冰目的。相比以往桨叶防除冰方法，该方法具有加热均匀，对叶片损伤小，安全可靠等优点。

1 基于气膜加热的风力机模型建立 1.1 气膜加热机制

本文防、除冰思想来源于在飞机机翼防冻问题上应用较广的二维槽缝热气喷射方法^[6]。利用高温射流在叶片表面形成一层薄的热气膜，使叶片前部与来流冷空气隔离，从而起到防冰作用。气膜加热有效度是度量气膜加热效果的一个重要参数，其定义如公式(1)：

式中, T_m为主流温度，设定为263K；T_c为射流温度, 设定为330K；T_w为翼型表面流体温度。计算可得η高于14.93%才能起到对叶片的保护作用(使T_w≥273K)，故在分析加热效果好坏时以14.93%为基准。气膜加热原理如图 1。

1.2 几何模型

本文叶片采用NACA四位数系列对称翼型建立模型^[7]，由于叶片前缘部位开设气膜孔，考虑到风力机对叶片强度及高升阻比的要求，为了提高叶片展向强度，故选择NACA-0020翼型。设定翼型弦长为0.25m，风轮直径D为3m，采用三叶片直翼式垂直轴达里厄风力机。

本文采用气膜加热方法在风力机叶片表面形成一层高温气体保护膜，故对翼型在叶片前缘分别开设三种不同孔型^[8-10]：圆柱孔、扩张孔、扇形孔。考虑叶片弦长大小，结合翼型前缘半径等翼型参数，参考气膜除冰孔专利，最终确定三种孔型如图 2所示。

1.3 网格划分

本文研究对象为三叶片直翼式垂直轴风力机，该翼型任意横断面截面完全相同，故在模拟不同孔型对叶片表面成膜情况及气动性能影响时可以简化为二维模型。为保证数值计算的高效性与可靠性，本文在不同计算域下进行建模并进行大量数值模拟，最终取计算域长13m和宽8m进行数值模拟研究，对比更大计算域下的大量数值模拟结果发现，其转矩系数C_m平均值偏差仅为0.0001。在对模型进行非稳态数值模拟研究时，将风力机流场划分为外部风场域、中部圆环旋转域、内部旋转域。流场主要采用分块结构化网格，中部圆环域采用非结构化网格，并对翼型表面及中部圆环域附近网格进行加密处理^[11]。整体计算域及局部网格划分如图 3所示。

2 数值方法 2.1 边界条件

计算模型的边界条件为：进口边界采用速度进口边界条件，给定来流风速为8m/s；部分滑移面采用了滑移网格技术^[12]，转动部分网格与静止部分网格交界面设置为interface，在moving type选项中设置运动区域旋转速度rpm；出口边界设为压力出口，采用充分发展的湍流条件；上下固壁区域边界采用对称边界symmetry。

2.2 控制方程及计算方法

根据对垂直轴风力机流场的模拟对比分析，本文计算湍流模型选用RNG k-ε两方程模型。RNG模型考虑了平均流动中的旋转及旋流流动。与标准k-ε模型相比可以更好地处理高应变率及流线弯曲程度较大的流体。垂直轴风力机为非定常流动，风轮旋转产生强烈扰动，选用RNG模型更加合理，而k-ε模型计算量小，稳定性高。压力速度耦合采用SIMPLE算法，基于有限容积的求解器，黏性项的离散项为中心差分，对流项为二阶迎风，收敛精度为10^-6，通用控制方程如式(2)所示：

式中，Ф为广义变量，Γ为相应于Ф的广义扩散系数，s为Ф对应的广义源项，在源项s中考虑了离心力、哥氏力和浮升力的影响。

3 数值模拟结果与分析 3.1 模型佐证

若流过风力机叶片扫掠面积的风功率为E，风力机获得的功率为P，则风能利用系数C_p，其定义如式(3)所示：

式中，M_avg为平均转矩，ω为风轮角速度，S为风轮扫掠面积。风能利用系数也称功率系数，59.3%为贝茨极限，是风力发电机组风能利用系数的最大值。

应用Fluent进行数值模拟时，当风力机运行达到稳定状态以后，风轮旋转一周转矩的平均值为M_avg，计算公式为：

式中，M(t)为任意时刻的瞬时转矩，n为转一周起始时间，m为转一周终止时间。

由于Fluent计算出的是不同时刻的瞬时转矩，故采用梯形法将公式(4)进行处理，处理后得到转矩平均值计算式为：

式中，M(t_i)为任意时刻t_i的瞬时转矩，t₀为转一周起始时间，t_n为转一周终止时间，Δt为计算所取的时间步长。

本文在来流风速固定条件下仅改变风轮圆周速度，也即使风力机在不同叶尖速比下运行，将由商业软件Fluent计算出来的C_p—λ曲线与通过双向多流管理论由MATLAB编程计算出的数据做了对比，结果如图 4所示。

由图 4以看出，对无孔翼型数值模拟的结果基本与根据多流管理论思想，利用MATLAB软件编程方法^[13]，计算出的C_p—λ分布吻合，说明本文数值模拟方案可行，并且很好地证明了采用滑移网格对直叶片垂直轴风力机进行非定常数值模拟能够很好地反映流场特性，从而证明后续数值模拟方案可靠。

从图 4中可以看出当叶尖速比为3左右时，也即风轮转速n=150r/min时，风轮的功率系数达到最佳。当叶尖速比偏高或偏低时，风轮的功率系数都将减小。由于本文主要研究翼型表面成膜情况，为保持变量的统一性，因此本文后续数值模拟均采用150r/min的风轮旋转速度。

3.2 不同孔型翼型气动性能分析

为了更好地分析翼型前缘不同射流孔及射流速度对翼型气动性能的影响，本文在来流风速固定为8m/s，风轮转速为150r/min条件下，分别对20m/s、25m/s、30m/s、35m/s、40m/s、45m/s、50m/s七种射流速度下的有孔翼型进行数值模拟^[14]。不同射流孔在不同射流速度下的转矩系数及功率系数如图 5所示。

由图 5可以看出，翼型前缘开射流孔后气动性能与无孔相比有一定的变化，开孔后的功率系数与转矩系数随孔射流速度增加呈下降趋势。当射流速度不超过35m/s时，由于开设的小孔很小，低射流速度下虽然有源源不断的热气进入主流，但是相对于整个流场来说，热气流量是微量的，孔射流在主流的作用下会随主流一起运动，且对主流的干扰较小，故对整个流场的影响非常有限，因此，随着射流速度增大，功率系数与转矩系数下降趋势较缓，且与无孔情况下的功率系数与转矩系数相比差距不大；当射流速度超过35m/s后，随着射流速度的增加，功率系数与转矩系数下降变快，气动性能恶化加快，这主要是因为当射流速度超过35m/s后，随着射流速度进一步增大，射流气体掺混进主流，对主流流场扰动强烈，进而造成较大能量损失；同时从图中对三种不同射流孔翼型气动性能的数值模拟结果中可以看出，三种孔型下的翼型功率系数及转矩系数基本相当，扇形孔的气动性能稍稍优于其它两种孔型的，这主要是扇形孔出流类似于扩张孔^[15]，孔射流气体流动较为稳定，因此本文在做流场及气膜加热有效度分析时主要针对扇形孔。

为了进一步探究开射流孔后翼型附近流场变化，在孔射流速度为35m/s下，取扇形孔翼型与无孔翼型在t=1.9s时的速度矢量及涡量分析，如图 6、图 7所示。

当翼型迎角增大到一定程度后，气流开始从翼型表面上分离。翼型气流分离有两种类型，一种是前缘分离，这种分离主要发生在薄翼型；另一种是后缘气流分离，这种分离一般在较厚的或有弯度的翼型上。本文气流分流主要是后一种，这从图 7也可以看出来。由于该位置翼型处于下风侧，计算得到该翼型迎角α=18.43°，已超出该翼型的失速迎角，翼型叶片已经失速，这从图 6中对无孔翼型流场的分析可以看出，并且气流在翼型中部偏下开始分离。从图 6中扇形孔翼型尾部流场可以看出，扇孔翼型附面层气流较无孔翼型提前离开壁面，且图 7中尾部涡范围相对较大。这可以解释为尽管射流流量不大，射流对翼型附近流场影响有限，但射流与主流掺混后扰动下游流场，改变翼型尾部气流形态，从而使翼型在原来的基础上更进一步发生轻失速状态，进而改变翼型升阻力，使转矩系数和功率系数降低，这印证了打孔后翼型气动性能的改变。通过对不同位置翼型流场的观察分析发现，这种影响主要发生在迎角较大的下风侧，对上风侧及迎角较小情况下的影响可以忽略不计，这也是当射流速度小于35m/s时有孔翼型总体气动性能减弱较小的的一个原因。

3.3 扇形孔翼型表面成膜情况

考虑射流速度过大对翼型气动性能影响较大，本文分析射流速度不高于35m/s情况下的扇形孔表面成膜情况，如图 8、图 9所示。

由图 8、图 9对比可以看出，当射流速度低于35m/s时，随着射流速度增加，翼型表面有效气膜覆盖范围变大，并且射流的中心线会稍稍远离壁面，对壁面的贴覆效果变差。考虑射流速度过大对翼型气动性能影响较大，本文对射流速度低于35m/s情况下的扇形孔表面成膜情况进行模拟分析发现，射流气体都很好的贴覆于翼型表面，因此综合考虑射流速度对翼型气动性能及成膜情况的影响，本文选择射流速度为30m/s。

气膜加热有效度η的大小直接反映射流的加热效果，已知加热有效度η为14.93%时叶片表面气膜温度为273K。从图中可以看出三翼型叶片在同一时刻表面温度高于273K的覆盖区域不同，而且翼型吸力面和压力面只有一侧有高温气膜覆盖，但是这并不影响气膜对翼型叶片的防护，因为从图中可以看出单个翼型叶片运动到不同位置时翼型吸力面和压力面交替得到高温气膜的加热。也即单个翼型叶片在旋转一周的过程中，叶片吸力面或压力面在不同时刻不同位置温度高于273K的覆盖长度不同，叶片前缘大部分区域在风轮旋转一周的过程中能交替得到有效的气膜保护。

4 结 论

本文主要采用数值模拟的方法，并且与双制动盘多流管理论模型相互佐证，针对不同孔型翼型在不同射流速度下进行详细的数值模拟，得到以下结论：

1) 通过对比双制动盘多流管模型，确定本文对无孔翼型所进行的数值模拟方案可行，并且得到翼型气动性能基本与根据MATLAB编程思想计算所得的结果吻合。

2) 通过分析翼型打孔专利及结合翼型弦长与风轮直径，确定三种孔型的具体数值参数，为本文分析翼型前缘气膜加热提供有力的保障。

3) 通过对不同孔型的数值模拟结果发现，孔型不同时，在相同设置条件下的翼型气动性能也不同，但差别很微弱，从对有孔与无孔的数值模拟结果可以看出，当扇形孔在来风速度固定为8m/s，风轮转速为150r/min，射流速度为30m/s时，有孔与无孔翼型的气动性能差别不大，但是有孔翼型前缘附近能得到很好的气膜防护。且随着旋转位置的不同，气膜覆盖长度和范围在交替改变，从而使翼型前缘附近大部分区域免于过冷气流的长时间直接冲刷，进而免于结冰的危害。

4) 气膜加热有效度η充分体现了开气膜孔后翼型表面气膜加热效果的好坏，本文以η高于14.93%为基准，从风轮三叶片在同一时刻的气膜加热有效度η的等值线图可以看出单个叶片在旋转一周的过程中翼型吸力面与压力面交替得到高温气膜的防护。

5) 本文提出的风力机叶片防除冰方法，相比于之前的方案，本文通过开设气膜孔，并根据不同情况下的风力机翼型气动特能及气膜防护情况，找到最佳工况点，可为将来风力机翼型叶片防除冰提供一定的参考。