水平轴风力机多涡格升力面涡尾迹法的研究

王渊博，李春,2，栾忠骏

(1.上海理工大学 能源与动力工程学院，上海 200093； 2．上海市动力工程多相流与传热重点实验室，上海 200093)



水平轴风力机多涡格升力面涡尾迹法的研究

王渊博1，李春1,2，栾忠骏1

(1.上海理工大学 能源与动力工程学院，上海 200093； 2．上海市动力工程多相流与传热重点实验室，上海 200093)

针对NREL Phase VI实验双叶片水平轴风力机，采用基于多涡格升力面的自由尾迹法模拟其低风速及高风速气动性能。由于高风速失速延迟导致尾缘分离滞后，建立Kirchhoff-Helmholz尾缘分离预估模型与Du-Selig失速延迟模型耦合的三维尾缘分离预估模型。计算低风速及高风速不同的偏航角工况，对比分析不同涡格数对模拟结果的影响。研究结果表明：涡格数对低风速工况影响甚小，对高风速影响很大，且采用两涡格的三维尾缘分离预估模型对法向力系数和弦向力系数的模拟最为精确。

水平轴风力机；多涡格升力面；自由涡尾迹；涡格数；偏航；失速延迟；尾缘分离

网络出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20160624.1127.014.html

为准确预估水平轴风力机在不同风速工况下所受气动力及气动力矩，以提高风力机设计水平，建立准确的气动模型是其根本前提[1]。常用的风力机气动研究方法有：动量叶素法[2]、涡尾迹法[3-4]和基于求解N-S方程的CFD方法。涡尾迹法计算精度高于动量叶素法，计算效率优于CFD方法，可准确用于计算低速不可压流动[5]。其中，自由尾迹法根据尾迹各点合速度更新尾迹位置，能模拟叶尖涡卷起和尾迹畸变现象[6]。文献[7]基于自由涡尾迹方法对四种动态来流风下的NREL Phase VI风力机性能和结构进行计算，发现非稳态来流风会导致风力机性能下降、疲劳载荷增加。文献[8]采用非线性升力面方法，建立轴向和偏航风条件下的水平轴风力机气动载荷及尾迹的数值模拟方法，证明了该方法可实现较好的工程数值模拟。文献[9]利用涡尾迹方法进行风力机气动载荷计算，结果表明大型风力机叶片几何非线性可较为明显地减小静气动弹性位移，同时降低动气动弹性的响应幅值。以上研究均采用单涡格升力面模型，当叶片弦长较大时误差较大，而本文时间步进自由尾迹法采用多涡格升力面模型，可精确模拟大弦长叶片的气动性能。此外，由于涡尾迹法基于不可压缩势流理论[10]，物面边界方程认为叶片表面流动不发生分离，适用于低风速小攻角工况，而对于高风速大攻角工况，叶片尾缘流动将发生分离，且分离点受三维旋转效应影响发生滞后现象，文献[7-9]并没有考虑这一问题。针对高风速大攻角工况，笔者建立了Kirchhoff-Helmholz尾缘分离预估模型[11]与Du-Selig三维失速延迟模型[12]耦合的三维尾缘分离预估模型，考虑高风速下各径向位置尾缘流动分离情况，并通过引入分离角偏转矩阵对原物面边界方程做出相应修正。在此基础上分别计算不同涡格数对各工况气动性能的影响，得到叶片最佳涡格数布置。

1　多涡格升力面模型

采用离散的涡分布升力面模型，即所谓涡格法(vortex lattice method, VLM)，将叶片分为径向和弦向若干个涡格，环量沿径向和弦向都存在变化，更好的描述了叶片表面的三维流动。

升力面涡格法网格划分采用在叶片弦向布置M个涡格，径向布置N个涡格[5-6]，如图1所示。沿径向的涡线为附着涡，弦向的涡线为自由涡，控制点布置在涡格中心。每个涡格四条涡线具有相同环量值，涡格环量方向根据右手定则与涡格法向量同向。

图1　叶片升力面涡格划分Fig.1　Vortex lattice distribution of blade lifting surface

涡尾迹法通过对每个涡格控制点满足物面边界条件求解叶片涡线环量。对于低风速工况，满足原物面边界方程式(1)即可:

(1)

式中：Vk为相对来流速度，Vind为尾迹及叶片涡线对控制点的诱导速度，n为涡格的法向量。对于高风速工况，需引入分离因子f进行修正，其定义为分离点距尾缘的弦向距离与弦长的比值，如图2所示。数值模拟时，对于高风速10m/s，取分离角为5°。

图2　高风速叶片尾缘流动分离示意图Fig.2　Trailing edge separation of high wind speed

在图2尾缘分离的情况中，分离点之前的涡格运用式(1)即可，而对于分离点之后的涡格，原物面边界方程已不适用。假设流动分离角η，将原法向力旋转分离角η，得到n′，令合速度与之垂直，通过在原物面边界方程中引入一分离角η偏转矩阵，以满足尾缘流动分离求解环量。修正后的物面边界方程：

(2)

为考虑粘性对诱导速度的影响，且避免计算诱导速度时出现奇值，采用Lamb-Ossen涡核模型，并引入时间补偿参数改变涡核起始半径，式(3)为随时间增大的涡核半径：

(3)

式中：α为经验系数，取常数1.256 43；ν为空气运动粘度；t为时间。根据文献[13]，时间补偿参数Sc取值0.1，扰动粘性系数δν取100。

2　三维尾缘分离预估模型

对于高风速工况，为考虑叶片尾缘分离，准确计算出叶片各径向位置分离点的位置，建立了Kirchhoff-Helmholz尾缘分离预估模型与Du-Selig失速延迟模型耦合的三维尾缘分离预估模型。

2.1Du-Selig三维失速延迟模型

文献[14]对比了6种三维失速延迟模型，其中Du-Selig模型的预测结果与实验值吻合最好，因此，本文采用该模型考虑风力机三维旋转效应。

Du和Selig通过求解旋转坐标系下风力机叶片表面三维积分边界层方程，在二维翼型升、阻力系数上添加一个增量进行三维修正：

(4)

(5)

式中：Cl,2D和Cd,2D分别为二维升、阻力系数；Cl,P=2π(α-α0)，α0为零升力攻角，α为实际攻角；Cd,0为二维零攻角阻力系数。经过与实验数据的关联，得到上式中fl、fd经验修正系数：

(6)

(7)

(8)

采用NRELPhaseVI实验叶片[15]为算例，以OSU大学在雷诺数Re=106下所测S809翼型的二维升、阻力系数实验数据为基础，在V=10 m/s，Ω=72 r/min工况下，将叶片各径向位置处弦长、当地尖速比代入式(4)和式(5)，计算经过Du-Selig三维失速延迟模型修正后的叶片各径向位置升、阻力系数，得到图3。分析三个径向位置，叶片无量纲半径(r/R)分别为0.3、0.47和0.8。

图3　Du-Selig模型修正后的三维升、阻力系数Fig.3　3D lift and drag coefficient corrected with Du-Selig model

从图3可知，Du-Selig三维失速延迟模型对小攻角线性增长区几乎无影响，而对于大攻角失速区，从叶尖到叶根，三维失速延迟模型的修正效果增强。

对于r/R=0.3，在大攻角失速区，升力系数几乎随攻角增大一直保持增长趋势，阻力系数相对原二维翼型阻力系数也降低很多，说明叶片径向流动对靠近叶根处的作用最为明显，尾缘分离点得到很大程度的推移；在r/R=0.45中间叶高处，失速攻角相应推迟，升力系数明显增大，阻力系数也相应减小；在r/R=0.8接近叶尖处，三维修正后的升、阻力系数与二维实验数据几乎一致，说明叶尖处周围流场影响较大，叶片径向流动效果明显减弱，几乎不造成尾缘分离点的推移。

2.2Kirchhoff-Helmholz尾缘分离预估模型

Kirchhoff-Helmholz模型基于二维翼型气动力系数，尾缘分离因子f可表示为翼型法向力系数Cn,2D的函数：

(9)

式中：α为当地攻角，α0为零法向力系数对应攻角，Cn,0对应法向力系数为零处的梯度。

将Kirchhoff-Helmholz尾缘分离预估模型与Du-Selig三维失速延迟模型耦合，考虑三维旋转情况下叶片径向流动对尾缘分离点位置的影响。在上节对升、阻力系数进行三维失速延迟修正后，通过下式得到相应修正后的三维法向力系数：

(10)

将Du-Selig模型得到的三维法向力系数Cn,3D代替式(9)Kirchhoff-Helmholz模型中的翼型二维法向力系数Cn,2D，得到耦合后的三维失速延迟特性的尾缘分离预估模型，用Kirchhoff-DuSelig模型代称。计算二维Kirchhoff-Helmholz模型及耦合后的三维Kirchhoff-DuSelig模型尾缘分离因子随攻角变化规律，并与二维翼型风洞实验数据[16]进行对比，如图4所示。

图4　预估模型分离因子Fig.4　Separation factor of prediction models

从图4可知，二维Kirchhoff-Helmholz模型分离因子与实验值非常接近，说明不考虑三维旋转效应时，Kirchhoff-Helmholz模型可准确模拟大攻角时二维翼型尾缘分离点位置。而耦合后的三维尾缘分离预估模型在不同径向处的分离因子与二维预估模型各不相同，体现出实际三维旋转效应使得叶片各径向位置的尾缘分离有不同程度的推移。攻角为-5°～8°，二维和三维耦合模型的分离因子相同，因为攻角在8°之前，升、阻力系数处于线性增长区，Du-Selig模型对升、阻力系数无影响，如图3所示。且在攻角小于6°时，分离因子都为0，尾缘流动不发生分离。随着攻角增大，分离因子相应增大，则尾缘流动分离发生的越快。

在r/R=0.3处，分离因子增大的幅度很小，因为叶片径向流动引起失速延迟现象，分离点推移，在靠近叶根处尤为明显。在r/R=0.8处，三维耦合模型的分离点预测几乎与原始二维模型一样，说明在靠近叶尖处，尾缘分离情况与二维翼型相差不大，没有发生明显的失速延迟现象。

3　低风速气动性能模拟

以NRELPhaseVI实验双叶片风力机为研究对象，叶片长5.03m，桨距角3°，风轮转速72r/min，叶片弦长及扭角数据参考文献[15]。计算风速5m/s，偏航角θ分别为10°、30°、45°和60°工况，叶片弦向涡格数M分布采用1、2和3，得到不同涡格数所模拟的叶片径向法向力系数Cn和弦向力系数Ct，图5为模拟值与实验数据对比图。

(a)V=5 m/s ,θ=10°

(b)V=5 m/s, θ=30°

(c)V=5 m/s ,θ=45°

(d)V=5 m/s ,θ=60°图5　风速5 m/s时不同偏航角径向Cn和Ct分布Fig.5　Normal force and tangent force coefficient of 5 m/s wind speed under different yaw angle

由图5可知，多涡格升力面法对低风速气动性能的模拟值与实验值吻合度很高。随着偏航角增大，径向法向力系数和弦向力系数减小，尤其当偏航角为60°时，能够模拟出靠近叶根处出现负攻角导致的负法向力系数，说明涡尾迹法对低风速气动性能的模拟精度很高。升力面采用不同的涡格数，对低风速模拟结果影响不大，因为低风速时，有效攻角很小，叶片尾缘未发生分离，弦向采用一个涡格已可准确表示叶片附着涡环量，亦无需采用三维尾缘分离预估模型求解分离点位置。

4　高风速气动性能模拟

4.1叶片气动力系数

模拟高风速10m/s，偏航角为10°和30°工况，采用基于Kirchhoff-DuSelig三维尾缘分离预估模型多涡格升力面法，涡格数M分别为1、2和3。图6对应叶片径向法向力系数和弦向力系数与实验数据对比图。

(a)V=10 m/s，θ=10°，Cn

(b)V=10 m/s，θ=10°，Ct

(c)V=10 m/s，θ=30°，Cn

(d)V=10 m/s，θ=30°，Ct图6　风速10 m/s时不同偏航角径向Cn和Ct分布Fig.6　Normal force and tangent force coefficient of 10 m/s wind speed under different yaw angle

图6说明，采用弦向涡格数M=2的Kirchhoff-DuSelig三维尾缘预估分离模型得到的法向力和弦向力模拟值与实验值吻合度最好。

当弦向涡格数M=1时，只在叶片弦向布置一个涡格略显粗糙，即认为叶片从1/4弦向处就发生分离，这样叶片表面分离过早，计算得到的法向力系数相对于实验值偏小。而Kirchhoff-DuSelig、M=2模型适当的增加了涡格数，更准确的描述尾缘分离点位置，所预测的法向力系数和弦向力系数与实验值最接近。Kirchhoff-DuSelig、M=3模型在叶片上布置了3个涡格，虽然能更准确的描述尾缘分离点位置，但是叶片上的涡格增多会引起对叶片附着涡诱导速度的增大，尤其是轴向诱导速度，从而导致靠近叶根模拟值偏高，靠近叶尖处模拟值偏低。

当叶片弦长增大时，涡格数需相应增加，具体弦长所对应的M值可通过与实验或精确CFD模拟结果对比求得。

4.2叶片分离因子

为了清晰直观地了解风速10m/s时，在不同的偏航角下，叶片表面各径向位置尾缘分离情况。图7给出了4个旋转周期内，在各个方位角叶片径向尾缘分离因子分布云图。

图7表明，各径向位置尾缘分离因子随当地攻角呈周期性变化，在靠近叶根处，由于三维旋转效应引起的失速延迟现象，尾缘分离点后移，分离因子较小，在靠近叶尖处，由于攻角很小，分离因子也较小。而在中间叶高处，失速延迟效应没有近叶根处明显，且攻角较大，因此，分离因子最大。而在每个旋转周期内，分离因子总在90°方位角取得最小值，在270°方位角取最大值，因为在90°方位角，来流速度在x轴的偏航分速度与相对旋转速度同向，造成攻角大幅减小，尾缘分离减弱。而在270°方位角，来流速度在x轴的偏航分速度与相对旋转速度反向，造成攻角大幅增大，尾缘分离加强。

(a)V=10 m/s，θ=10°

(b)V=10 m/s，θ=30°图7　风速10 m/s时不同偏航角周向尾缘分离因子分布Fig.7　Azimuthal trailing edge separation factor of 10 m/s wind speed under different yaw angle

4.3叶片气动力矩

根据图7中尾缘分离因子随方位角变化关系，分析尾缘分离因子对叶片主轴扭矩和挥舞力矩的影响。图8为风速10m/s，偏航角分别为10°和30°，主轴扭矩和挥舞力矩在4个周期内随方位角变化的关系。

(a)主轴扭矩

(b)挥舞力矩图8　风速10 m/s不同偏航角主轴扭矩和挥舞力矩Fig.8　Shaft torque and flap-wise moment of 10 m/s under different yaw angle

计算叶片气动力矩需查询升阻力系数表，其中升力系数呈正弦分布，阻力系数不呈正弦分布，由于攻角较小，挥舞力矩受升力系数影响大，主轴力矩受阻力系数影响大，因此挥舞力矩呈正弦波动，而主轴扭矩呈锯齿状周期波动。此外，随着偏航角的增大，主轴扭矩和挥舞力矩相应减小，受图7中尾缘分离因子的影响，主轴扭矩在方位角90°左右取最大，在方位角270°左右尾缘分离较剧烈，扭矩值最小，因此，Kirchhoff-DuSelig三维耦合尾缘分离预估模型较准确地预估了不同方位角各径向位置尾缘分离点，从而表现出尾缘分离情况对主轴扭矩的影响，主轴扭矩直接影响风力机的输出功率。

5　结论

1)多涡格升力面模型对低风速工况的气动性能模拟精度很高，且涡格数对结果影响不大；

2)涡格数对高风速气动性能模拟结果影响较大，采用两涡格的Kirchhoff-DuSelig三维尾缘分离预估模型，考虑了尾缘分离所引起的气动力系数降低，且适当增加涡格能准确地表示尾缘分离点位置，模拟结果与实验值最为接近。

3)得出风速10m/s时不同偏航角下，叶片径向分离因子随方位角分布规律。中间叶高处分离因子较大，且在每个旋转周期内，分离因子在90°方位角最小，在270°方位角最大，因此，叶片主轴扭矩在270°方位角输出最小。

[1]李春, 叶舟, 高伟, 等. 现代大型风力机设计原理[M]. 上海: 上海科学技术出版社, 2013: 110-147.

LI Chun, YE Zhou, GAO Wei, et al. Modern large-scale wind turbine design principle[M]. Shanghai: Shanghai Scientific and Technical Publishers, 2013: 110-147.

[2]OKULOV V L, VAN KUIK G A M. The Betz-Joukowsky limit: on the contribution to rotor aerodynamics by the British, German and Russian scientific schools[J]. Wind energy, 2012, 15(2): 335-344.

[3]KOMNINOS K. Modeling considerations of the optimum rotor using vortex method[D]. Denmark: Technical University of Denmark, 2008: 7-119.

[4]SHENKAR R. Design and optimization of planar and nonplanar wind turbine blades using vortex methods[D]. Denmark: Technical University of Denmark, 2010: 1-86.

[5]BOUATEM A, ALMERS A, BOUTAMMACHTE N. Load evaluation of horizontal-axis wind turbine rotor using coupled Beddoes near-wake model and free-wake method[J]. International journal of energy and environmental engineering, 2013, 4: 35.

[6]CHKIR S. Unsteady loads evaluation for a wind turbine rotor using free wake method[J]. Energy procedia, 2011, 6: 777-785.

[7]周文平, 唐胜利, 吕红. 风剪切和动态来流对水平轴风力机尾迹和气动性能的影响[J]. 中国电机工程学报, 2012, 32(14): 122-127.

ZHOU Weiping, TANG Shengli, LV Hong. Effect of transient wind shear and dynamic inflow on the wake structure and performance of horizontal axis wind turbine[J]. Proceedings of the CSEE, 2012, 32(14): 122-127.

[8]仇永兴, 康顺. 基于非线性升力面方法的风力机尾迹数值模拟[J]. 工程热物理学报, 2012, 33(12): 2072-2075.

QIU Yongxing, KANG Shun. Numerical simulation of wind turbine wake based on nonlinear lifting surface method[J]. Journal of engineering thermophysics, 2012, 33(12): 2072-2075.

[9]唐迪, 陆志良, 郭同庆. 大型水平轴风力机叶片气动弹性计算[J]. 应用数学和力学, 2013, 34(10): 1091-1097.

TANG Di, LU Zhiliang, GUO Tongqing. Aeroelastic analysis of large horizontal axis wind turbine blades[J]. Applied mathematics and mechanics, 2013, 34(10): 1091-1097.

[10]陈浩, 张亮. 漂浮式潮流电站载体水动力性能分析[J]. 应用科技, 2014, 41(1): 65-68.

CHEN Hao, ZHANG Liang. Hydrodynamic performance analysis of floating current power station[J]. Applied science and technology, 2014, 41(1): 65-68.

[11]LEISHMAN J G, BEDDOES T S. A generalised model for airfoil unsteady aerodynamic behaviour and dynamic stall using the indicial method[C]//Proceedings of the 42nd Annual forum of the American Helicopter Society. Washington DC, 1986: 243-265.

[12]杜朝辉. 水平轴风力涡轮设计与性能预估方法的三维失速延迟模型——Ⅲ. 模型改进[J]. 太阳能学报, 2000, 21(2): 145-150.

DU Zhaohui. A 3-D stall-delay model for HAWT performance prediction: Ⅲ. model improvement[J]. Acta energiae solaris sinica, 2000, 21(2): 145-150.

[13]SANT T, VAN KUIK G, VAN BUSSEL G J W. Estimating the angle of attack from blade pressure measurements on the NREL phase VI rotor using a free wake vortex model: axial conditions[J]. Wind energy, 2006, 9(6): 549-577.

[14]许波峰, 王同光, 张震宇. 风力机三维旋转效应模型研究[J]. 太阳能学报, 2014, 35(4): 562-567.

XU Bofeng, WANG Tongguang, ZHANG Zhenyu. Investigation on three dimensional rotational effect model for wind turbine[J]. Acta energiae solaris sinica, 2014, 35(4): 562-567.

[15]Hand M M, Simms D A, Fingersh L J, et al. Unsteady aerodynamics Experiment Phase VI: Wind Tunnel Test Configurations and Available Data Campaigns[M]. Golden, Colorado: National Renewable Energy Laboratory, 2001: 5-45.

[16]SHENG Wanan, GALBRAITH R A M, COTON F N. On the S809 airfoil′s unsteady aerodynamic characteristics[J]. Wind energy, 2009, 12(8): 752-767.

本文引用格式：

王渊博，李春，栾忠骏.水平轴风力机多涡格升力面涡尾迹法的研究[J]. 哈尔滨工程大学学报, 2016, 37(8): 1070-1075.

WANG Yuanbo, LI Chun, LUAN Zhongjun.Research on multi-lattice lifting surface vortex wake method in horizontal-axis wind turbine[J]. Journal of Harbin Engineering University, 2016, 37(8): 1070-1075.

Research on multi-lattice lifting surface vortex wake method in a horizontal-axis wind turbine

WANG Yuanbo1, LI Chun1,2, LUAN Zhongjun1

(1.School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China；2. Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, Shanghai 200093, China)

In this study, we used a free vortex wake method based on a multi-lattice lifting surface to simulate aerodynamic performance at low and high wind speeds of the NREL Phase VI experimental two-blade horizontal axis wind turbine. To address the trailing edge separation lag caused by stall delay at high wind speeds, we coupled the Kirchhoff-Helmholz trailing edge separation prediction model and the Du-Selig 3D static stall delay model. We then calculated the aerodynamic performance at both low and high wind speeds for different yaw angle conditions. We also discuss the influence of different vortex lattice numbers on calculation accuracy. Our simulation results demonstrate that, at low wind speeds, the vortex lattice number has little influence, while it has great impact at high wind speeds. We obtained the best simulation results of the normal force and tangent force coefficients for the set of two vortex lattices on the lifting surface.

horizontal-axis wind turbine; multi vortex lattice lifting surface; free wake method; number of vortex lattices; yaw; stall delay; trailing edge separation

2015-06-12.网络出版日期：2016-06-24.

国家自然科学基金项目(E51176129)；上海市科学技术委员会项目资助(13DZ2260900)；教育部高等学校博士学科点专项科研基金(博导类)项目(20123120110008).

王渊博(1991-)，男，硕士研究生;

李春(1963-)，男，教授，博士生导师.

王渊博，E-mail:shlgwyb@163.com.

10.11990/jheu.201506038

TK83

A

1006-7043(2016)08-1070-06