新型双风轮风力机气动特性的三维流场数值模拟

周云龙， 杨承志， 李律万

（1.东北电力大学 能源与动力工程学院，吉林132012；2.中科院广州能源研究所，广州510640）

风力发电在洁净可再生能源的开发利用中占有 重要的地位，风轮是风力发电机捕捉风能的关键部件，它的设计直接决定了风力机的出力和风能转换效率，因此风轮的研究是风力机整机设计的重点.传统的单风轮风力机叶片少则迎风面积小，叶片多则转速低，都不能获得较高的输出功率，而新型的双风轮风力机（见图1）则弥补了单风轮风力机的缺陷，在发电机的两侧装设前后2个风轮，后风轮捕捉前风轮的漏风，增大了风力机的捕风面积，并且各自保持了较高的转速，同时应用了新型、高效的对转式异步发电机［1］，前、后2个风轮分别带动发电机的内、外转子进行相对旋转运动，增大了线圈切割磁感线的速度，从而获得了较大的发电功率.

图1 新型双风轮风力机结构示意图Fig.1 Structural diagram of a new double－rotor wind turbine

然而，目前双风轮风力机作为一种新型的风力机，各国对其的研究还处于起步阶段.日本的Kubo［2］以及国内的安源等［3］对2个风轮同侧布置的双风轮风力机已经有了初步的研究成果，但是这种布置方式在实际运行中存在缺陷，例如转动部件摩擦严重、整个风力机维持受力平衡困难以及2个风轮由于间距小、可调性差而造成碰撞风险等，影响了风力机的运行安全和使用寿命.2个风轮异侧布置的新型双风轮风力机可以有效地解决同侧布置存在的问题，并且可以同样实现增大发电功率的目的.不过，这方面的研究刚开始，其流动特性尚不被人们熟知.近些年，随着空气动力学和CFD数值模拟理论的发展，Fluent软件被广泛应用于风力机流场的流动特性分析等领域，并得到了充分的验证［4－7］.笔者利用Fluent软件对这种2个风轮异侧布置的新型双风轮风力机的流场进行探索性的研究，并结合实验测试结果加以分析.

1 风力机模型的建立

所研究的风力机为小型风力机，叶径R为1.5 m，以吉林某地10年的风速数据为依据，确定风力机的额定风速为11.26m／s，选用NACA4412翼型数据，根据各种风轮的叶片数确定相应叶片的叶尖速比，利用Wilson叶片优化设计方法，基于Matlab自主编程完成风力机叶片的设计.风力机的塔架高度为3.3m，双风轮风力机的前、后风轮间距为0.4 m，利用 Pro／E建立各类型风力机模型［8－9］，模型如图2所示.

图2 各类型风力机模型Fig.2 Models of different wind turbines

根据模拟的需要，构建了5种类型的风力机：3叶片单风轮风力机，记为N3；前轮3叶片、后轮3叶片双风轮风力机，记为N前3，后3；前轮3叶片、后轮4叶片双风轮风力机，记为N前3，后4；前轮3叶片、后轮5叶片双风轮风力机，记为N前3，后5；前轮3叶片、后轮6叶片双风轮风力机，记为N前3，后6.

2 数值计算方法

2.1 计算域和计算网格划分

将风力机模型导入前处理软件Gambit，并对流场进行网格划分，为了能完整地查看风力机流场的气动性能，计算域的横向上风向长为10 R，下风向长为20 R，纵向关于风力机中心面对称，宽为10 R.由于研究对象是扭曲叶片数目较多的风力机，物理模型比较复杂，尽量减少计算网格数有利于数值计算.在保证网格质量的前提下，在对计算域进行网格划分时采取如下处理：一方面在不同的区域采用不同的网格生成方法，首先利用Size Function合理分布风力机各边和面的网格，然后利用TGrid方法在风轮的旋转区域以及风力机附近的小区域生成加密网格，并利用Cooper方法在风力机前后两边的大区域生成简单网格；另一方面设置较大的网格节距，网格节距的增大在一定程度上增加了数值模拟结果曲线的拐点数目，大范围的使用可能会降低曲线的平滑性，但保持合理的限度并不会影响数值计算的收敛性和结果的正确性.网格划分结果如图3所示，各类型风力机流场的总网格数均为350万左右，Equiangle skew小于0.89，满足了网格数量和质量的要求.

图3 各类型风力机模型的网格划分Fig.3 Grid division for models of different wind turbines

2.2 边界条件

对各类型风力机在额定工况（风速v取额定风速11.26m／s）和5个工况（v分别取2m／s、5m／s、8m／s、10m／s和12m／s）下的运行情况进行数值模拟.模型入口为VELOCITY＿INLET，出口为OUTFLOW.风轮旋转流体选用MRF模型，单风轮风力机的风轮附近为逆时针方向的旋转域，双风轮风力机的前、后风轮附近分别为逆时针和顺时针方向的旋转域.假设后风轮旋转受前风轮的影响较小，根据式（1），可求得不同类型的风轮在不同风速下的转速［10］，在Fluent软件中进行旋转域旋转速度设定

式中：n为风轮转速，r／min；λ为叶尖速比，3、4、5、6叶片风轮的λ一般分别取6、4.5、4、3；v为风轮上风向大气流速，m／s；R 为风轮半径，取1.5m.

2.3 求解方法

假定无传热现象，模拟采用Segregated（分离式）求解器隐式求解三维稳态不可压缩流动，紊流模型使用SST k－ω 模型［11］，压力－速度耦合采用Simplic算法，对流项差分格式采用二阶迎风格式，参考压力为1.01×105Pa.求解器的时间步长设为“自动调整时间步长”，最大迭代次数为2 000步，残差类型为均方根，前后计算的残差余量设为1×10－5.

2.4 模拟结果处理方法

风力机性能分析的主要参数包括风力机风轮输出的机械功率P、风力机发电功率P0和风能利用系数CP，其计算公式［10－12］（对于单风轮风力机，可以认为是只有前风轮的双风轮风力机进行计算）分别为

式中：M前、M后分别为前、后风轮的转矩，N·m；n前、n后分别为前、后风轮的转速，r／min；η1为发电机效率，%；η2为传动效率，%，小型实验性风力机一般取η1η2＝72%；ρ为通常状况下的大气密度，取值为1.205kg／m3；A 为风轮的扫掠面积，m2.

由式（4）可推导出双风轮风力机较单风轮风力机的发电功率增加百分比αPo与风能利用系数增加百分比αCp的表达式为

式中：Po，单、Po，双分别为单风轮、双风轮风力机的发电功率，W.

3 模拟结果及分析

3.1 静压与速度云图分析

在额定工况下，风力机风轮最有代表性的位置0.85 R处，5种类型风力机环形截面上的静压（相对压力）云图与速度（绝对速度）云图见图4.

对比图4左侧各类型风力机的静压云图可知，大气绕流各个叶片翼型时，在翼型的上表面附近大气流动受阻，动能转化为静压能，形成表压为700 Pa左右的高压区；在翼型的下表面边界层发生分离，附近形成了表压为－50Pa左右的低压区，翼型这2个表面间的压力差提供了叶片持续快速旋转的升力.这样，前风轮叶片后的低压区与后风轮叶片前的高压区综合作用形成了2个风轮间相互影响的流场.当双风轮风力机后风轮的叶片数目少时，这种影响作用较小，随着后风轮叶片数目的增加，这种影响作用越明显.

比较图4右侧各类型风力机的速度云图可知，各风力机风轮旋转区域入口处的大气流速接近额定风速11.26m／s，而出口出现了风速为3m／s左右的低风速区域.与单风轮风力机相比，随着双风轮风力机后风轮叶片数目的增加，低风速区域增大，风速降低，并且无传热现象，这表明大气通过对各个叶片的做功推动了风轮的旋转，动能发生了形式上的转化，双风轮风力机转化的能量较为明显，从而使输出较多的电能成为可能.

3.2 等速线图分析

在额定工况下，5种类型风力机在沿大气来流方向（x轴方向）与前风轮某叶片展向（叶片的叶根轴线方向）所形成的平面截面、沿大气来流方向与前风轮该叶片展向垂直截面上的等速线图见图5.

图4 各类型风力机的静压云图（左）与速度云图（右）Fig.4 Static pressure（left）and speed（right）contour of different wind turbines

图5 各类型风力机的等速线图Fig.5 Equal speed lines of different wind turbines

受风轮旋转的影响，各类型风力机都在风轮前部出现了空气压缩流场（阻塞效应），但相差不大，可见双风轮风力机虽然增加了后风轮，但只要叶片数目不是很多，其对前风轮流场的影响就不会很大，这对保持前风轮较大的输出功率是非常重要的.各类型风力机的旋转风轮都在叶片径向外侧的后部形成了漩涡，双风轮风力机受到后风轮反向旋转的影响，这部分涡流的大小比单风轮风力机的涡流有所减小，同时位置向下游延迟，这一区域为高速等值线区，风能较大，由于风力机叶片的出力主要是由这部分叶片贡献，涡流的减小和延迟非常有利于后风轮捕捉前风轮的漏风.最后，大气经过各类型双风轮风力机的前风轮叶片时都会产生涡旋回流，并对后风轮造成一定程度的影响，在某些区域这种影响对后风轮的旋转是有利的，在另外一些区域则不利.只有通过进一步的性能参数分析，确定后风轮叶片的合理数目，才能尽可能地增加对有利区域流场的利用，从而合理布置后风轮，使得后风轮可以较好地捕捉到前风轮的漏风，因而在增大迎风面积的同时，保持2个风轮较高的转速，高效地实现风能的两级利用.

3.3 湍动能云图分析

在额定工况下，5种类型风力机流场的湍动能云图如图6所示.图6中单风轮风力机的下风向流场的湍动能很快趋近于来流风速流场的湍动能E＝2.0m2／s2，而各类型双风轮风力机的下游流场出现了湍动能为5.0m2／s2左右的高湍动能区域，且随后风轮叶片数目的增加，该区域不断扩大，N前3，后6类型双风轮风力机下游流场部分区域的湍动能甚至达到了8.0m2／s2左右.由此可见，双风轮风力机较单风轮风力机的湍动能有明显地增大，且后风轮叶片数越多，流场的湍动能越大，相应的湍流强度也越大，N前3，后6类型双风轮风力机的湍动能已经非常大了.所以，为了保证风力机的安全稳定运行，双风轮风力机在设计加工时，后风轮叶片数目不能过多，而且在设计加工过程中应提高设备的机械强度.

图6 各类型风力机的湍动能云图Fig.6 Turbulent kinetic energy contours of different wind turbines

3.4 性能参数分析

通过数值计算，得到了各类型风力机风轮在不同工况下的输出转矩M，如表1所示.由表1可知，各类型的双风轮风力机由于后风轮的影响，前风轮的转矩比N3类型单风轮风力机风轮的转矩小.但是，当后风轮叶片布置合理时，在大气经过前风轮后的涡旋回流作用下，后风轮对前风轮漏风的利用能力有了较大的提高，其输出转矩甚至超过前风轮.前、后风轮的共同作用使双风轮风力机总的输出转矩有所增加，进而在很大程度上提高了发电功率和增大了风能利用系数.

表1 各类型风力机风轮在不同工况下的输出转矩Tab.1 Output torque Mof wind rotors for different wind turbines under varying working conditions

将相应工况下各类型风力机风轮的转速n分别代入式（2）～式（5），可得到各工况下评价各类型风力机的性能分析参数，即发电功率Po、风能利用系数CP和双风轮较单风轮风力机的发电功率增加百分比及风能利用系数增加百分比，然后在Matlab中进行数据拟合得到各个参数的曲线（见图7）.

由图7可知，数值模拟得出的 N前3，后3、N前3，后4、N前3，后5与 N前3，后6类型双风轮风力机的性能分析参数较N3类型单风轮风力机的性能分析参数均有一定的提高，其风能利用系数分别为29.4%、32.0%、34.3%和27.7%左右，相对于单风轮风力机分别增大了43.7%、56.3%、67.4%和35.3%，这与文献［13］中的实验测试结果一致，说明了模拟结果的正确性.

由数值模拟结果可知，由于后风轮叶片数目较多，N前3，后6类型双风轮风力机流场的湍动能很大，使得运行稳定性较低，另外，受到前风轮的不利影响增加，其发电功率与风能利用系数相对于其他类型的双风轮风力机反而有减小的趋势，不适合实际应用，所以后风轮叶片数目应以5叶片以下为宜.其中，N前3，后5类型双风轮风力机的流场相对平稳，发电功率及风能利用系数较大，可以作为主要研究方向.

图7 模拟得出的各类型风力机的性能分析参数Fig.7 Performance parameters of various simulated wind turbines

4 结 论

（1）对比单风轮风力机的流场，当新型双风轮风力机的后风轮叶片数目合理时，后风轮对前风轮的影响较小，而且能够充分地利用大气经过前风轮所产生的涡旋回流的有利影响，有效地捕捉到前风轮的漏风，这样可以在增大迎风面积的同时，保持2个风轮较高的转速，进而达到提高发电功率和增大风能利用系数的目的.当双风轮风力机的后风轮在5叶片以内时，叶片数目越多，发电功率和风能利用系数的增加越明显.

（2）双风轮风力机比单风轮风力机的流场更加复杂，随着后风轮叶片数目的增加，湍流强度变大，风力机运行稳定性变差.基于运行稳定性和高效性，N前3，后5类型双风轮风力机将成为主要研究方向.

（3）数值模拟的性能参数变化趋势与实验测试的结果一致，说明与单风轮风力机相比，新型的双风轮风力机在一定程度上提高了发电功率和增大了风能利用系数.

［1］BOOKER J D，MELLOR P H，WROBEL R，et al.A compact，high efficiency contra－rotating generator suitable for wind turbines in the urban environment［J］.Renewable Energy，2010，35（9）：2027－2033.

［2］KUBO K，KANEMOTO T.Development of intelligent wind turbine unite with tandem wind rotors and double rotational armatures［J］.Journal of Fluid Science and Technology，2008，3（3）：370－378.

［3］安源，韩凤琴，久保田乔，等.反转叶轮风力机非定常尾流对下游环境的影响［J］.水电能源科学，2010，28（10）：158－160.AN Yuan，HAN Fengqin，TAKASHI Kubota，et al.Effect of unsteady wake from counter－rotating wind rotor on downstream environment［J］.Water Resources and Power，2010，28（10）：158－160.

［4］BARTHELEMIE R J，RATHMANN O，FRANDSEN S T，et al.Modelling and measurements of wakes in large wind farms［J］.Journal of Physics：Conference Series，2007，75（3）：12－19.

［5］DUQUETTE M M，VISSER K D.Numerical implications of solidity and blade number on rotor performance of horizontal－axis wind turbines［J］.Journal of Solar Energy Engineering，2003，125（4）：425－432.

［6］李宇红，张庆麟.风力机风轮三维流动特性与气动性能的数值分析［J］.太阳能学报，2008，29（9）：1172－1174.LI Yuhong，ZHANG Qinglin.Numerical simulation of flow field and aerodynamic performance of a wind turbine blade［J］.Acta Energiae Solaris Sinica，2008，29（9）：1172－1174.

［7］祝贺，徐建源，腾云，等.风力机风轮气动性能三维流场数值模拟［J］.中国电机工程学报，2010，30（17）：85－90.ZHU He，XU Jianyuan，TENG Yun，et al.3Dflow field numerical aerodynamic performance test of wind turbine rotor［J］.Proceedings of the CSEE，2010，30（17）：85－90.

［8］张果宇，蒋劲，刘长陆.风力发电机整机气动性能数值模拟计算与仿真研究［J］.华东电力，2009，37（3）：449－451.ZHANG Guoyu，JIANG Jin，LIU Changlu.Numerical simulation of aerodynamic performance for wind turines［J］.East China Electric Power，2009，37（3）：449－451.

［9］李国宁，杨福赠，杜白石，等.基于 MATLAB与Pro／E的风力机风轮设计及造型［J］.机械设计，2009，26（6）：3－7.LI Guoning，YANG Fuzeng，DU Baishi，et al.Design and modeling of wind wheel of wind mill based on MATLAB and Pro／E［J］.Journal of Machine Design，2009，26（6）：3－7.

［10］何显富，卢霞，杨跃进，等.风力机设计、制造与运行［M］.北京：化学工业出版社，2009：73－81.

［11］郭婷婷，吴殿文，王成荫，等.风力发电机叶片预弯设计及其数值研究［J］.动力工程学报，2010，30（6）：450－455.GUO Tingting，WU Dianwen，WANG Chengyin，et al.Pre－bend design and numerical simulation of wind turbo－generator［J］.Journal of Chinese Society of Power Engineering，2010，30（6）：450－455.

［12］李少华，岳巍澎，匡青峰，等.双机组风力机尾流互扰及阵列的数值模拟［J］.中国电机工程学报，2011，31（5）：101－107.LI Shaohua，YUE Weipeng，KANG Qingfeng，et al.Numerical simulation of wake interaction and array of double wind turbine［J］.Proceedings of the CSEE，2011，31（5）：101－107.

［13］周云龙，杨承志，岳巍澎.一种新型双风轮风力发电装置的特性分析［J］.中国电力，2011，44（4）：75－78.ZHOU Yunlong，YANG Chengzhi，YUE Weipeng.Characteristic analysis for a new type of double－rotor wind turbine［J］.Electric Power，2011，44（4）：75－78.