郭小锋^1,2, 陈进¹, 李松林¹
1. 重庆大学 机械传动国家重点实验室， 重庆 400044
2. 中原工学院 机电学院， 河南 郑州 450007
收稿日期: 2014-11-11
基金项目: 国家高技术研究发展计划项目(2012AA051301)； 国家自然科学基金资助项目(51175526).
作者简介: 郭小锋(1983-)，男，河南南阳人，重庆大学博士研究生，中原工学院讲师。

摘要: 基于修正的叶素动量理论提出了一种适应于水平轴风力机叶片的静气弹性分析方法；运用该方法对不同稳态风速下的叶片变形曲线进行研究，揭示了叶片预弯设计的机理，并给出了一种预弯型线的设计方法.通过对预弯和非预弯两种叶片在相同风速下的静气弹性变形和气动功率的对比和分析表明，预弯叶片的气动性能受轴向变形、扭转变形和叶片截面形状变化三个方面的影响，叶片弹性扭转变形使叶片实际气动功率偏离设计值.该研究为预弯叶片的设计提供了理论依据和方法.
关键词：风力机预弯叶片静气弹性预弯机理预弯型线设计
Mechanism and Method of Pre-bending Design for the Wind Turbine Blades
GUO Xiao-feng^1,2, CHEN Jin¹, LI Song-lin¹
1. State Key Laboratory of Mechanic Transmission, Chongqing University, Chongqing 400044, China
2. School of Mechanical Science & Engineering, Zhongyuan University of Technology, Zhengzhou 450007, China
Corresponding author: GUO Xiao-feng, E-mail: mexf@foxmail.com
Abstract: A new static aeroelastic analysis method for the HAWT was proposed based on the modified blade element momentum (BEM) theory. The deformation curves of the blade under different steady wind speeds were analyzed. The prebending design principles and a profile design method were developed. The static aeroelastic analyzing result of the prebending and unprebending blade shows that, the aerodynamic performance is influenced by the cross section shape change, blade axial and torsional deformation. The prebending blade’s actual power has certain decrease because of the blade torsional deformation. This study provides theoretical basis and design method for the blade.
Key Words: wind turbinepre-bending bladestatic aeroelasticpre-bending mechanismpre-bending profile design
为了减小大型水平轴风力机的发电成本，避免叶片因弹性变形后与塔筒发生干涉和碰撞，MW级叶片多进行了预弯设计.在预弯叶片研究方面，LM公司最早提出了预弯设计的概念，并申请了发明专利^[1]；Riziotis基于Euler-Bernoulli 梁理论，对预弯叶片的振型和气弹变形等问题进行了研究^[2]；Kooijman等给出了6 MW风力机及其配套预弯叶片的设计原型和相关参数^[3]；郭婷婷等对不同预弯尺寸的叶片进行CFD气动分析^[4]，研究了预弯尺寸对叶片发电功率输出的影响.然而，在预弯机理、预弯型线表达及预弯叶片的参数化设计方法方面，相关研究文献很少.风力机叶片的气动外形和铺层结构之间存在着气弹耦合关系，对叶片静气弹性平衡点下的气弹性能进行分析可预测叶片实际工作中的气动性能.
1 参数化静气弹性分析方法1.1 叶片外形的参数化表达叶片的外形是若干个翼型截面扫掠而成的，对各个截面的空间坐标进行参数化转换就能实现对整个叶片的参数化表达.叶片各截面的三维坐标变换过程可用式(1)表达，其转换过程如图 1所示，翼型1为标准翼型，将翼型1的参考坐标原点变换至翼型的桨矩轴中心得到翼型2，将翼型2按当地弦长进行等比变换得到翼型3，将翼型3绕参考坐标原点按当地扭角进行旋转变换得到翼型4.
图 1(Fig. 1)

图 1 叶片截面坐标变换Fig.1 Coordinate transformation of the blade section

在翼型4到翼型5的预弯变换过程中，为避免叶片上表面产生扭曲变形和屈曲失效，翼型的预弯位移变换不是将翼型4由桨心平移到对应的预弯型线对齐点处，而是将翼型4由平行于弦长线的1/2厚度线与过翼型桨矩轴中心的弦长垂直线的交点(等厚度点)平移到预弯对齐点处，即将图 2中的等厚度点预弯曲线上的预弯对齐点重合.变换后的翼型5仍与z轴垂直，为保证最终的变换翼型面与预弯型线垂直，需要将翼型5进行扭转变换得到最终的截面形状6.
图 2(Fig. 2)

图 2 叶片截面预弯对齐图Fig.2 Aligned figure of the blade section

(1)

式中：[X₀,Y₀,Z₀] 为基本翼型坐标；[X_i,Y_i,Z_i] 为变换后的叶片截面坐标；T₁,T₂,T₃,T₄,T₅为变换矩阵；X_R为相对于翼型前缘点的桨矩轴中心坐标；c为当地弦长；β为当地扭角；Y_B为预弯对齐点与z轴的垂直距离；α为翼型预弯扭转角.
1.2 铺层结构模型风力机叶片主要由5种材料铺设而成，分别为：三轴向玻璃纤维布(Triax)、双轴向玻璃纤维布(Biax)、单轴向玻璃纤维布(UD)、腹板泡沫(Foam-W) 和翼面泡沫(Foam-A).叶片的截面结构主要有箱型梁、I型梁、II型梁等形式，本文研究的叶片采用I型梁结构，其截面结构如图 3所示.
图 3(Fig. 3)

图 3 叶片截面结构铺层Fig.3 Structure of the blade section

叶片的上下翼面具有相似的铺层结构，从结构功能上，可将其划分为6个区域：前缘加强(LER)、前缘面板(LEP)、梁冒(CAP)、后缘面板(TEP)、后缘单向布加强(TEUD)和后缘加强(TER).在叶片的上下翼面之间安装有腹板，用以支撑上下翼面和传递剪切载荷，叶片的主腹板安装于上下翼面的梁帽之间，在叶根部位只有主腹板，在其他部位为了增加尾缘的强度，增设了副腹板.
1.3 分布载荷施加叶片在实际工作中，气动外形和铺层结构存在着复杂的流固耦合关系，为了准确分析气动外形与铺层结构的相互影响，一些学者运用流体计算软件如FLUENT对叶片进行气动和气弹性能分析^[5-6]，但这种方法需要较多的时间开销，中小型计算机很难胜任.为了实现对叶片气动结构耦合性能的快速分析，本文采用MATLAB和ANSYS APDL混合编程的方法对叶片进行载荷施加^[7].对于叶片上的每个翼型截面，当初始条件确定以后，其表面压力系数可以通过风力机翼型气动分析专业软件RFOIL计算得到.将本文研究的1.5 MW风力机叶片划分为30个截段，调用RFOIL软件从叶尖至叶根依次计算每个截段翼型的压力系数分布C_p，并用式(2)求得每个叶片截面上的压力分布p.为了在ANSYS APDL程序中实现对叶片每个单元的压力施加，首先对压力分布进行多项式拟合，然后读取叶片截段上每个单元的弦向坐标，插值得到叶片有限单元上的压力值.

(2)

其中：ρ为空气密度,取1.205 kg/m³；p _∞ 为标准大气压强；v_rel为作用在叶片上的相对速度.
对于1.5 MW叶片，距叶根20 ~ 21 m处截段的离散压力分布和拟合压力分布如图 4所示，施加拟合压力的叶片截面受力情况如图 5所示.
图 4(Fig. 4)

图 4 叶片展向位置r=20 m处截面的压力分布Fig.4 Pressure distribution at r=20 m along the blade span-wise

图 5(Fig. 5)

图 5 叶片展向位置r=20 m处截面的气动载荷加载Fig.5 Aerodynamic load applied at r=20 m along the blade span-wise

此外，叶片在转动中还要承受重力和离心力.重力可以通过指定重力加速度实现，对于惯性力，本文将叶片各个截段的离心力平均加载在各截段的节点上.
1.4 静气弹性分析叶片的气动力与结构弹性力之间存在着耦合关系，从静态的角度分析，气动载荷使叶片产生弹性变形，弹性变形使叶片的气动载荷减小，两者相互影响，最终使得叶片的弹性变形维持在静气弹性平衡点状态；从结构动力学的角度分析，叶片实际工作中的弹性变形不是静止的，而是围绕着静态平衡点作动态变化，对叶片在静态平衡点下的气弹变形及气动性能进行分析(即静气弹性分析)，可近似等价于对叶片实际工作中平均气动性能的分析，实现对叶片在实际工况下的气动性能评估.
基于叶片的外形和铺层结构参数，通过MATLAB和ANSYS APDL语言编程参数化生成叶片的有限元模型，在生成的叶片模型下，采用迭代算法求解叶片的静气动弹性平衡状态.叶片的静气弹性分析流程可概况如下：第1步，调用RFOIL软件计算叶片各个截面的气动压力分布，对压力分布进行多项式拟合，按叶片截段各单元的弦向坐标进行面压力加载和静态求解分析，从计算结果文件中提取各个截面的扭转变形角度β₀；第2步，考虑叶片各截面翼型形状变化、扭角变化及叶片长度变化三个因素对叶片气动载荷的影响，计算得到叶片各截面上的气动载荷分布，按步骤1的方法在新的载荷分布下对叶片进行加载和分析，得到叶片各个截面的扭转变形角度β₁；重复步骤2计算得到叶片各个截面的扭转变形角度β₂，β₃…，直至前后两次计算中叶片各个截面处扭转变形角度的差值不大于给定容许偏差为止.
2 预弯原理与预弯型线设计本文研究的1.5 MW风力机叶片工作于IEC/GL IIA风区^[8]，切入风速为3 m/s，切出风速为 25 m/s，风场年平均风速为6 ~ 8.5 m/s；叶片总长为37.5 m，风轮额定转速为17 r/min.叶片在几种典型稳态风况下的变形情况如图 6所示，从图中可以看出，当风速为6.5 m/s时，叶片变形后的桨矩轴线基本与水平轴重合；当风速达到额定风速(11 m/s)时，叶尖的变形量达到4.307 m.分析叶片的变形情况，根据Euler-Bernoulli梁理论可以做出如下推理：将预弯叶片在6.5 m/s风速下的载荷施加在具有相同气动布局和铺层结构的非预弯叶片上，非预弯叶片的变形形状应和原预弯叶片的形状关于水平轴近似呈对称形态，即具有相同气动外形和铺层结构的非预弯叶片在特征风速下的变形曲线关于水平轴的镜像线可作为叶片的预弯型线.
图 6(Fig. 6)

图 6 不同风速下的叶片变形Fig.6 Blade deflection under various wind speeds

为了验证该推理，构建与原预弯叶片具有相同气动布局和铺层结构的非预弯曲叶片(后文中简称之为非预弯叶片)，对非预弯叶片在6.5 m/s风速下进行静气弹性分析，得到了非预弯叶片的变形曲线如图 7所示.分析该曲线可见，该曲线的叶尖挠度为1.745 m，与原叶片的叶尖预弯量1.72 m基本相同，叶片中后部(R/2 ~ R，R为叶片总长)的变形曲线与原叶片预弯型线关于水平轴的镜像线基本吻合；在叶片根部，非预弯叶片变形线与原预弯型线镜像线有较大的差别，原预弯叶片在前R/3段预弯量为0，而非预弯叶片的变形线却有较大的相对变形.而事实上，预弯叶片在叶片根部前段一般不进行预弯，分析其原因，一是由于叶片在这些部位的变形较小，二是为了便于生产制造和储存运输.因此，在求得非预弯叶片在特征风速下的变形后，须对该曲线修正以满足实际生产制造工艺需要.
图 7(Fig. 7)

图 7 叶片预弯型线的生成Fig.7 Generation of the pre-bending line

在对非预弯叶片特征风速下的弯曲型线型进行修正时，为了实现对特征风速下叶片变形曲线的拟合，同时又能满足预弯叶片的设计要求，本文提出了一种新的预弯型线表达方法，其表达式为

(3)

式中：A为幂函数预弯型线的指数；R为叶片的水平长度；B为叶尖的预弯量；b(x) 为叶片各个截面外形在图 1中由翼型4到翼型5垂直于z轴的偏移量.预弯叶片的弯曲变形起始于距离叶片根部约R/3处，预弯段的型线函数为幂函数，由幂函数的性质可知，式(3)表示的预弯型线在预弯起始点处与水平轴(即桨矩轴)相切，因此能保证叶片在预弯过渡段处实现平滑过渡.通过改变式(3)中A和B的值就可以得到不同的预弯型线.
对6.5 m/s风速下非预弯叶片的变形曲线运用式(3)进行拟合，得到最优参数A=2.18，B=1.792，拟合型线关于x轴的镜像即为设计的预弯型线.比较分析原叶片的预弯型线和新设计的预弯型线，沿叶片展向各点处最大相对偏差为4.1%，最大相对偏差在很小的范围内，验证了该预弯型线设计方法的可行性和正确性.
3 预弯叶片的气弹性能分析预弯叶片在气弹变形影响下，气动性能受轴向变形、扭转变形和叶片截面形状变化的影响会与初始设计值产生偏离，考虑三种因素的影响，通过对预弯叶片和非预弯叶片的静气弹性变形及叶片静气弹性平衡点下发电功率的分析和讨论，研究预弯设计对叶片气动性能的影响.
表 1为叶片在6.5和11 m/s的稳态风速下对原预弯叶片和与原预弯叶片具有相同气动布局和铺层结构的非预弯叶片进行气弹分析的结果，从中可以看出在相同风速下，两只叶片在y方向的变形几乎相同，在z方向上，预弯叶片在6.5和11 m/s风速下的变形量分别为0.207和 0.482 m，叶片变形后的长度较原叶片增长了0.552%和1.283%，而非预弯叶片仅增长了0.181%和0.450%.分析预弯叶片在z方向上有较大变形的原因是在与非预弯叶片具有相同叶片长度(叶片在z轴上的投影长度)的情况下，预弯叶片的曲线长度较非预弯叶片长，在特征风速6.5和11 m/s下，叶片的曲线形状得以伸展，因此预弯叶片较非预弯叶片具有较大的z向变形.
表 1(Table 1)

表 1 气弹变形和气动功率对比Table 1 Comparison of the aeroelastic deflection and power 风速 叶片 类型 y方向 变形量/m z方向 变形量/m 气动功率 /kW (m·s-1) 6.5 预弯 1.609 0.207 359.7 非预弯 1.604 -0.007 357.1 11 预弯 5.649 -0.482 1 516.1 非预弯 5.582 -0.017 1 493.2

表 1 气弹变形和气动功率对比 Table 1 Comparison of the aeroelastic deflection and power

图 8为两只叶片在6.5和11 m/s的稳态风速下的扭转变形图，预弯叶片与非预弯叶片在相同风速下叶片各截面的扭转变形基本相同，叶片扭转变形后各翼型截面绕气动中心顺时针变形，从而减小了实际攻角.
图 8(Fig. 8)

图 8 叶片扭转变形图Fig.8 Torsional deflection of the blade

在叶片的气动功率计算中，叶片长度的增加有助于气动功率的提高，而叶片截面翼型形状变化和扭角变化则会使叶片的气动功率减小.预弯叶片在6.5和11 m/s风速下的设计气动功率分别为363.4和1 500 kW，基于普朗特叶尖修正模型和Shen修正模型^[9]对经典的叶素动量理论进行改进，依据式(4)计算得到预弯叶片和非预弯叶片的实际气动功率如表 1所示.

(4)

式中：R为叶片直径；B为叶片个数；ρ为空气密度；c为当地弦长；v_rel为当地风速；C_t为切向力系数；F₁为Shen修正因子.
从表 1中可以看出，在考虑扭转变形时，因为叶片z向变形、叶片截面翼型形状和扭转变形三个因素的相互影响，非预弯叶片的气动功率在两种风速下均比设计值小，但相差不大；而预弯叶片的气动功率在额定风速时由于叶片z方向的变形，叶片长度增加，叶片的实际气动功率较设计功率有一定的增加.
4 结论对风力机叶片的预弯机理进行了研究，提出了一种实用的风力机叶片预弯曲型线设计方法.通过对预弯叶片的静气弹性分析可知，预弯叶片的轴向弹性变形有利于实际功率的提高，而叶片截面的扭转变形则使风力机的实际气动功率减小，在叶片设计中对叶片外形进行预扭修形设计可以提高叶片的年发电量.
参考文献

[1]	Olesen E E.Method of manufacturing pre-bent wind turbine blades:US,8216500B2[P].2012-07-10.(0)
[2]	Riziotis V A,Voutsinas S G,Manolas D I,et al.Aeroelastic analysis of pre-curved rotor blades[C]//Proceedings of EWEC 2010.Warsaw:EWEC,2010:2-10.(0)
[3]	Kooijman H J T,Lindenburg C,Winkelaar D,et al.DOWEC 6 MW pre-design:aero-elastic modelling of the DOWEC 6 MW pre-design in PHATAS[R].Petten:Energy Research Center of the Netherlands,2003.(0)
[4]	郭婷婷, 吴殿文, 张强, 等. 风力机叶片预弯技术的数值模拟[J].太阳能学报, 2011, 32(7): 1011–1025. ( Guo Ting-ting, Wu Dian-wen, Zhang Qiang, et al. Numerical researches of the wind turbine blade’s preflex technology[J].Acta Energiae Solaris Sinica, 2011, 32(7): 1011–1025.)(0)
[5]	Hansen M O L, S?rensen J N, Voutsinas S, et al. State of the art in wind turbine aerodynamics and aeroelasticity[J].Progress in Aerospace Sciences, 2006, 42: 285–330.(0)
[6]	Rethore P E,Niels N,S?rensen J N,et al.Mexico wind tunnel and wind turbine modelled in CFD[C]//29th AIAA Applied Aerodynamics Conference.Hawaii,2011:2011-3373.(0)
[7]	Lee Y J, Jhan Y T. Fluid-structure interaction of FRP wind turbine blades under aerodynamic effect[J].Composites Part B, 2012, 43(5): 2180–2191.(0)
[8]	Lloyd G.Guideline for the certification of wind turbines[S].Hamburg:Germanischer Lloyd,2010.(0)
[9]	Shen W Z, Mikkelsen R, S?rensen J N, et al. Tip loss correction for wind turbines computations[J].Wind Energy, 2005, 8(4): 457–475.(0)