2018, Vol. 50
引用本文 | doi: 10.16356/j.1005-2615.2018.03.013 |
2. 中国航天空气动力技术研究院, 北京, 100074;
3. 山东交通学院航空学院, 济南, 250357
2. China Academy of Aerospace Aerodynamics, Beijing, 100074, China;
3. School of Aeronautics, Shandong Jiaotong University, Jinan, 250357, China
风力机的气动性能研究对现代风能利用有重要意义。对于大多数未解决的问题，都有相应发展并经过验证的工程规律。所有的这些规律适用范围都相当有限，越来越需要用物理机制和物理模型来取代这些工程规律，这也是世界风能气动研究机构转向更基础研究的一个重要原因[1]。
不管是在研究过程还是在设计过程，准确地预估风力机风轮叶片的动力输出具有特别重要的意义[2]。要做到对动力输出能够准确地预估，不管是使用带有三维修正的改进动量叶素理论(Blade-element theory, BEMT)方法，还是基于涡系理论的升力线、面法[3, 4]，或者基于NS方程的计算流体力学(Computational fluid dynamics, CFD)方法[5]，都需要对叶片周围的流场结构细节有所了解。同时叶片周围流场结构决定了风力机风轮的气动性能，直接影响着风力机的风能利用率、结构动态特性、叶片载荷分布和噪声，因此有必要采用多种手段对叶片周围流场的分布规律进行研究。
文献[6]对美国国家能源部可再生能源实验室(National renewable energy laboratory, NREL)第Ⅵ期风轮叶片的流场进行了数值模拟研究。Hirahara等[7]利用PIV技术在风洞闭口段对一个微型风力机进行了来流、叶尖和尾迹区的流场测量。胡丹梅等[8-11]使用PIV技术和三维热线风速仪研究了旋转叶片的尾迹。高志鹰等[12]在直流式低速风洞中采用PIV手段对水平轴风力机叶尖涡流动进行了研究。
本文对某三维水平轴风力机风轮模型进行风洞试验和数值模拟研究，试验在来流风速为8 m/s，尖速比，λ=4, λ=9.8工况下，主要采用PIV系统对风轮叶片5个径向截面的瞬时流场进行了测量，通过多次采样取平均获得了截面上速度场的分布，并与数值计算结果进行了比较。
1 试验方案 1.1 试验件试验风轮包含3个标准STORK 5.0 WPX叶片组成的定桨距风轮，直径为0.75mm。该风轮为上风向水平轴风力机风轮，采用8个翼型截面确定叶片。由于PIV试验需要，在叶片表面喷涂亚光黑漆，如图 1所示。
|
图 1 风轮UG模型(左)与实物(右) Figure 1 Model of turbine |
1.2 试验装置
试验系统基本结构简图如图 2所示。
|
图 2 试验系统 Figure 2 Test system |
(1) 试验是在南京航空航天大学能源与动力学院叶轮机实验室的直流式低速风洞的开口实验段进行。该实验段直径为900 mm，最大风速10 m/s。
(2) 风轮的转速由一台测功机控制。通过转速传感器、扭矩传感器测得当前的转速和扭矩，将转速和扭矩信号通过测功机处理后来控制磁粉制动器的电流从而对风轮转速进行控制，如图 3所示。
|
图 3 转速测量与控制系统 Figure 3 Rotating speed measure and control system |
(3) 本文使用的PIV系统是美国TSI公司的产品(图 4)，该系统能够用于一个空间的二维或三维速度场的瞬态信息，测量对象为气体或水流流动。本系统由以下主要部件组成：双钇铝石榴石(Yttrium aluminum garnet, YAG)脉冲激光器，片光源透镜组，专用于PIV测量的200万像素CCD相机，高速图像数据接口板，同步控制器，立体图像分析系统以及计算机系统。
|
图 4 PIV系统的激光器(左)及光臂(右) Figure 4 Laser generator (left) and waveguide arm (right) in PIV |
(4) 由于风力机叶片是转动部件，为了实现叶片周围流场的PIV测量，必须使得所测叶片通过测量平面的频率与PIV系统激光脉冲发射频率相同，并控制两者的相位差。为此，使用了本实验室设计的一套简易的锁相触发装置(图 5)。该装置可根据安装在风轮轴上的光电编码器输出的脉冲数进行计数及延时触发。在试验时设定当风轮叶片每转一圈时，该装置输出一个5 V的TTL正跳变脉冲给PIV系统的同步控制器，从而进行图像采集。同时该装置还有延迟触发功能，并且延迟时间可以微调，可以确保将叶片锁定在PIV测量范围的中心。
|
图 5 锁相触发装置 Figure 5 Phase locked trigger device |
(5) 示踪粒子的选择。示踪粒子的光色散性和跟随性的好坏决定了最终成像品质和速度矢量计算的成败。这两个特性是互相制约的：色散性随尺寸增大而增大，但是尺寸越大粒子随流体流动的跟随性就越差。同时，粒子浓度也非常关键，它与激光强度的衰减、流体的流动、信噪比、测量区域大小和光学系统的放大率等诸多因素有关[13]。权衡上述因素，本试验中采用舞台烟机来产生示踪粒子，示踪粒子为烟油在高温下雾化而成的颗粒微团，粒子直径在1~2 μm之间。在风洞壁面开口将导流管深入风洞中进行粒子释放，其位置及PIV测量现场如图 6所示。
|
图 6 PIV系统 Figure 6 PIV system |
1.3 试验方案
来流风速为8 m/s。首先通过测功机控制系统将风轮转速从最大转速调至最小转速获得1条该风速下的特性线。在该特性线上选取λ=4和λ=9.8(对应转速分别为800 r/min和2000 r/min)两个点来进行PIV测量。为获得叶片周围流场的情况，沿叶片叶高方向布置5个PIV测量截面S1~S5(图 7)，其相对位置见表 1。
|
图 7 PIV系统组成和结构 Figure 7 Sections for PIV measurement |
|
表 1 PIV测量径向截面分布 Table 1 Section for measurement in PIV |
2 数值模拟计算结果及分析 2.1 计算网格和边界条件
风轮叶片流场数值模拟网格在ICEM CFD中生成，由UG生成的叶片造型导入到ICEM中，计算域如图 8所示。整个计算域轴向长度为4倍的风轮直径。固壁边界直径为5倍的风轮直径，经计算验证该尺寸下洞壁效应可以忽略。计算采用非结构化四面体网格，在叶片前缘、尾缘和叶尖部分进行加密，并在叶片表面生成3层棱柱体网格来控制第1层网格厚度，如图 9所示。整个区域计算网格单元数为173万。采用稳态计算，湍流模型采用k-ω SST模型。流场区域采用旋转参考坐标系。
|
图 8 计算网格及边界设置图 Figure 8 Computational domain and boundary |
|
图 9 叶片区域网格细节图 Figure 9 Details on blade |
2.2 计算结果
图 10给出了风轮叶片在4个来流风速下的风能利用系数随尖速比的变化曲线(Cp-λ曲线)。从图 10可以看出，各型风轮在不同风速下的Cp-λ曲线比较集中。因为Cp-λ曲线反映叶轮的几何特性，对于定桨距风力机，风轮几何构型一定，Cp-λ曲线基本确定，不会随转速或风速变化产生明显的变化，并且只要在风轮的正常工作转速下，风轮都会在同一尖速比达到其最佳工作状态。因此在试验过程中仅对某一个风速进行测量即可。
|
图 10 数值模拟的Cp-λ曲线 Figure 10 CP-λ curves by numerical simulation |
3 数值计算与PIV试验结果对比分析 3.1 风轮叶片气动特性对比分析
图 11，12给出了测得的风轮在8 m/s风速下的性能曲线及其与CFD数值计算结果的比较。图 11为风轮的扭矩随转速的变化曲线，从图 11, 12中可以看出：风轮的输出扭矩随着转速呈先增大后减小的变化趋势；在小转速n＜1 400 r/min时，风轮的计算结果与试验结果存在较大差别，而在大转速n＞1 400 r/min时, 两者比较相符。
|
图 11 8m/s风速风轮下的T-n曲线 Figure 11 T-n at wind speed with 8 m/s |
|
图 12 8 m/s风速风轮下的Cp-λ曲线 Figure 12 Cp-λ at wind speed with 8 m/s |
由测得的扭矩计算出功率，进而计算出风轮的风能利用效率为Cp。图 12为风轮在8 m/s风速下的Cp-λ曲线与CFD数值计算结果的对比。在Cp＞7.4范围，试验测得的Cp与计算值符合性很好，在λ＜7.4范围由于测得的扭矩与计算值有偏差，那么相应的Cp也比计算值低。
3.2 风轮叶片周围流场分析在来流风速为8 m/s下，λ=4和λ=9.8时叶片周围的流场如图 13，14所示。
|
图 13 风速为8 m/s，λ=4时S1~S5截面上的速度矢量图 Figure 13 Velocity vector distribution on S1—S5 at V=8 m/s and λ=4 |
|
图 14 风速为8 m/s, λ=9.8时S1~S5截面上的速度矢量图 Figure 14 Velocity vector distribution on S1—S5 at V=8 m/s and λ=9.8 |
对比PIV获得的截面上的速度矢量分布和CFD计算的相对速度矢量分布，两者在趋势和数值上都比较一致。对比尖速比为λ=4和λ=9.8两个工况点下的叶片绕流流场可以发现，当λ=4即风轮转速小、风速相对较大时，试验测得的S1, S2截面叶背区域有明显的流动分离, 如图 13(a)所示，由此产生了较大的尾迹区二次流损失，从而影响了风轮的输出扭矩。CFD计算仅在S1截面有叶背有分离，S2截面叶背速度较小但是没有分离，如图 13(b)所示。风轮工作在λ=9.8时，如图 14所示，从试验和计算得出的速度矢量分布上看，5个测量截面上都没有出现流动分离，风轮的风能利用率也较高，这与图 12风轮的特性曲线是相一致的。
图 15是λ=4时，5个截面相对速度矢量和云图分布。在该组截面上，获得风轮流场中轴向速度和切向速度分量。从图 15速度云图可以看出S1~S5截面，叶片尾迹区均比较明显，其中S4，S5低速区的范围较大。从速度矢量分布来看S1和S2截面的尾迹区终有明显的流动分离。从风轮运行的工况来看，此时转速较低，气流流经叶片的相对迎角较大，因此在牵连速度较低的叶根部更容易发生流动分离。
|
图 15 风速为8 m/s, λ=4时S1~S5的平均流场相对速度矢量和云图分布 Figure 15 Velocity vector distribution on S1—S5 at V=8 m/s and λ=4 |
图 16是λ=9.8时，沿径向分布的5个截面的平均流场相对速度分布。与λ=4不同的是，在高转速下，由于从叶根到叶尖的牵连速度都比较高，在叶片尾迹区仅能看到有稍微偏低的速度分布，但是没有流动分离出现。
|
图 16 风速为8 m/s, λ=9.8时S1~S5截面上的速度矢量图 Figure 16 Velocity vector distribution on S1—S5 at V=8 m/s and λ=9.8 |
初步分析，造成在小转速(或者小λ)范围两者差别较大的原因有:
(1) 计算分离流动时计算方法的影响。
在小转速区对应的叶片当地迎角较大，在叶背气流有较大的分离。对于这种情况，在数值模拟的时候要选择非定常方法并结合更加适用于分离流的湍流模型，如大涡模拟(Large eddy simulation, LES)模型。而本文计算方法在全工况下采用的都是定常的计算方法，对湍流的模拟采用的也全都是k-ω SST模型。
(2) 低雷诺数的影响。
在小转速区，转速较低导致线速度较低，因而叶片的当地雷诺数较低。当雷诺数低于临界雷诺数值时会造成附面层向湍流的转捩推迟，而层流附面层比湍流附面层更容易分离，导致气流损失增大，风轮实际输出扭矩偏小。而在数值模拟计算的时候，使用的湍流模型无法模拟近壁附近的层流，更没有模拟出附面层的转捩。
综合上述两个原因，发现数值模拟过程对风轮叶片上的流动损失预估偏小，而使得数值模拟得到的风轮输出扭矩偏大、风能利用率偏高。
3.3 考虑转捩过程的数值模拟在这种情况下开展数值模拟计算，如果想更加准确地反应流场的流动现象，不能忽略层流流动，要引入转捩模型。Sezer采用LES的方法对风轮进行了模拟[14]，LES具有模拟转捩过程的能力，同时可以得到风轮的非定常流场，但是该方法耗费的资源大，由于计算条件限制，本文没有做此方面的尝试。早期的研究通过经验或半经验方法来确定转捩，比如经验关联方法，eN法以及基于间歇因子的预测方法。但这些方法都不是基于当地变量的，其模拟精度依赖经验参数，也难以与现有的CFD计算方法相结合。Menter等[15]提出完全基于当地变量关联的γ-Reθ转捩模型，该方法可以与现代CFD计算，方法比如非结构化网格相结合，也适用于大规模并行计算。该模型基于两个输运方程：一个是间歇因子的，另一个是动量厚度雷诺数的。该模型在风力机翼型转捩过程的模拟中得到了应用，本文针对小尺寸模型风轮的低雷诺数问题，在湍流模型中引入γ-Reθ转捩模型来改进计算方法。
图 17给出了采用带有转捩模型的计算结果和不带转捩模型的扭矩的对比。在考虑转捩后，在小于1 800 r/min转速范围内，计算所得的扭矩值相对全湍流状态偏低，在大于1 800 r/min范围，两者几乎一样。图 18中结果显示，在λ=2~7的范围内，考虑转捩后Cp与试验值更加接近。在小转速范围，以800 r/min为例进行分析，如图 19所示，为50%叶高截面的湍动能分布。可以看出，在吸力面附近，采用转捩模型计算获得的流场湍动能要低于全湍流的。
|
图 17 8 m/s风速风轮下的T-n曲线 Figure 17 T-n at wind speed with 8 m/s |
|
图 18 8 m/s风速风轮下的Cp-λ曲线 Figure 18 Cp-λ at wind speed with 8 m/s |
|
图 19 50%叶高截面的湍动能分布 Figure 19 Turbulent kinetic energy at 50% span |
图 20为50%叶高截面翼型前缘吸力面细节，从流线分布来看，带有转捩模型的结果，因为考虑了层流流动，其抗分离能力更弱，导致在更靠近前缘的位置发生分离，其分离区较大，因此使得风轮的输出功率降低，更加接近实际的流动状态。
|
图 20 50%叶高截面的吸力面前缘流线 Figure 20 Suction side streamline at 50% span |
4 结论
本文使用CFD软件对某三维水平轴风力机模型进行了数值模拟，得到了风轮叶片在不同工况下的气动特性参数和速度场分布，并且通过试验的手段，采用PIV系统测量翼型在相同雷诺数不同迎角下的流场，并且通过分析对比，得到以下结论：
(1) 在8 m/s来流风速下，尖速比大于7.4时，试验测得的风轮扭矩和风能利用率与数值模拟结果趋于一致；尖速比小于7.4时，试验测得的扭矩值低于计算值，其风能利用效率也较低。
(2) 将PIV试验测量得到的速度矢量分布和CFD数值计算结果进行比较，在λ=4时，PIV测得靠近叶根的两个截面S1，S2在叶背有明显的流动分离，计算结果中仅在S1截面叶背存在流动分离，S2截面叶背存在低速区。
(3) 在λ=9.8时，5个测量截面上都没有出现流动分离。
(4) 初步尝试采用γ-Reθ转捩模型后，流场的湍动能要低于全湍流的，尤其是在小转速情况下，如800 r/min下，50%叶高处其分离位置比全湍流计算结果提前。考虑层流向湍流的转捩过程，使得数值模拟结果更加接近试验值。
| [1] |
SNEL H.
Review of the present status of rotor aerodynamics[J]. Wind Energy, 1998, 1(S1): 46–69.
|
| [2] |
杜朝辉, 刘富斌.
旋转对风力涡轮叶片边界层的影响[J]. 上海交通大学学报, 1999, 33(8): 962–973.
DU Chaohui, LIU Fubin. Effect of rotation on blade boundary layer of wind turbine[J]. Journal of Shanghai Jiaotong University, 1999, 33(8): 962–973. |
| [3] |
吕超, 王同光, 许波峰.
三维动态失速模型在风力机气动特性计算中的应用[J]. 南京航空航天大学学报, 2011, 43(5): 707–712.
LV Chao, WANG Tongguang, XU Bofeng. Application of 3-D dynamic stall model in wind turbine aerodynamic performance prediction[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(5): 707–712. |
| [4] |
许波峰, 王同光.
基于自由涡尾迹法和面元法全耦合风力机气动特性计算[J]. 南京航空航天大学学报, 2011, 43(5): 592–597.
XU Bofeng, WANG Tongguang. Wind turbine aerodynamic performance prediction basedon free-wake/panel model coupled method[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(5): 592–597. |
| [5] |
钟伟, 王同光.
风力机叶尖涡的数值模拟[J]. 南京航空航天大学学报, 2011, 43(5): 640–644.
ZHONG Wei, WANG Tongguang. Numerical analysis of the wind turbine blade-tip vortex[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(5): 640–644. |
| [6] |
SØRENSEN N N, MICHELSEN J A, SCHRECK S.
Navier-Stokes predictions of the NREL phase Ⅵ rotor in the NASA Ames 80 ft×120 ft wind tunnel[J]. Wind Energy, 2002, 5(2/3): 151–169.
|
| [7] |
HIRAHARA H, ZAKIR HOSSAIN M, et al.
Testing basic performance of a very small wind turbine designed for multi-purposes[J]. Renewable Energy, 2005, 30(99): 1279–1297.
|
| [8] |
胡丹梅, 田杰, 杜朝辉.
水平轴风力机尾迹流场PIV实验研究[J]. 太阳能学报, 2007, 28(2): 200–206.
HU Danmei, TIAN Jie, DU Chaohui. PIV experimental study on the wake flow of horizontal axis wind turbine model[J]. Acta Energiae Solaris Sinica, 2007, 28(2): 200–206. |
| [9] |
胡丹梅, 竺晓程, 杜朝辉.
水平轴风力机尾迹三维流场的热线测量[J]. 太阳能学报, 2006, 27(1): 7–13.
HU Danmei, ZHU Xiaocheng, DU Chaohui. Wake measurements of a model horizontal-axis wind turbine using hot-wire technique[J]. Acta Energiae Solaris Sinica, 2006, 27(1): 7–13. |
| [10] |
胡丹梅, 竺晓程, 杜朝辉.
水平轴风力机尾迹的测量与分析[J]. 动力工程, 2006, 26(5): 751–760.
HU Danmei, Zhu Xiaochen, DU Chaohui. Measurement and analysis of wake behind horizontally orientated air turbines[J]. Journal of Power Engineering, 2006, 26(5): 751–760. |
| [11] |
胡丹梅, 欧阳华, 杜朝辉.
水平轴风力机尾迹流场试验[J]. 太阳能学报, 2006, 27(6): 606–612.
HU Danmei, OUYANG Hua, DU Chaohui. An experimental study of the wake structure of a model horizontal-axis wind turbine[J]. Acta Energiae Solaris Sinica, 2006, 27(6): 606–612. |
| [12] |
高志鹰, 汪建文, 东雪青, 等.
水平轴风力机叶尖涡流动的PIV测试[J]. 工程热物理学报, 2010, 31(3): 414–418.
GAO Zhiying, WANG Jianwen, DONG Xueqing, et al. PIV experiment on tip vortex flow of horizontal axis wind turbine[J]. Journal of Engineering Thermophysics, 2010, 31(3): 414–418. |
| [13] |
李军.
PIV技术在涡轮静叶测量中的应用[J]. 激光与红外, 2008, 38(2): 105–108.
LI Jun. Application of PIV measurement technology on turbomachinery[J]. Laser & Infrared, 2008, 38(2): 105–108. |
| [14] |
SEZERU. 3-D time-accurate CFD simulations of wind turbine rotor flow fields[C]//44th AIAA Aerospace Sciences Meeting. Reno, NV, United States: [s. n.], 2016: 394-417. |
| [15] |
MENTER F R, LANGTRY R B.
A correlation-based transition model using local variables-Part Ⅰ:Model formulation[J]. Journal of Turbomachinery, 2006, 128(3): 413–422.
DOI:10.1115/1.2184352
|