吕计男，刘子强，赵 玲，冉景洪

（中国航天空气动力技术研究院，北京 100074）

0 引 言

随着风力发电事业的蓬勃发展，风力机的功率已经由初期的千瓦级小型风力机上升为现在的兆瓦级大型风力机，风力机叶片长度和重量都呈现快速增加的趋势。风力机叶片的气动弹性分析需要气动力载荷及结构数据。目前，气动力载荷计算主要以动量叶素理论（BEM）及计算空气动力学（CFD）方法为主。BEM方法简单快捷，应用方便，理论成熟度高，一直是国内外风力机气动载荷及气动弹性问题研究的重要手段。CFD方法可以提供精确的流场描述，但对计算机硬件的要求很高，计算效率较低，特别是对气动弹性分析所需的非定常气动力计算尚存在很大难度。

本文选用南京航空航天大学设计的NH1500风轮叶片作为研究对象。该叶片长度40.5m，三维构型由南航优化后NH02系列翼型叠加组成。针对大型风力机叶片气动弹性响应计算，本课题组开发了一套基于BEM方法的气动弹性分析程序，该程序可以考虑偏航工况、计及叶尖及轮毂损失、考虑翼型厚度和宽度、风剪切的影响，并可以进行失速修正。本文针对NH1500叶片，比较了本程序计算的气动载荷结果与商用软件GH Bladed的结果，验证了程序可靠性，并给出风力机整机结构在气动载荷作用下的响应结果。

1 叶片动量叶素理论气动力建模

1.1 理论基础

动量理论［1］主要用来估算风力机的理想输出功率，该理论描述的是风力机叶轮整体与空气气流的能量交换情况。叶素理论［1］将风力机叶片简化为有限个叶素叠加的形式，叶片三维气动特性由叶素的气动性能沿径向积分得到，该理论以叶素表面气体流动的状态分析叶片的受力和能量交换。考虑到采用动量理论与叶素理论描述的叶片气动力与力矩在物理本质上的同一性，可以建立平衡方程得出动量叶素理论中轴向及切向速度诱导因子的显式表达式［2］。

在理论的基础上，可以通过叶尖损失系数、轮毂损失系数考虑由于叶尖及根部气流沿叶片产生的二次流动而引起的力矩的减小［1，3］。

叶栅理论（Cascade Theory）研究叶片厚度和宽度的影响导致的叶素攻角的改变［1］。当周向速度诱导因子大于0.5时，风轮工作在湍流尾流状态，这时需要对动量理论进行修正。本文的计算模型采用的是 Glauert提出的修正表达式［4－5］：

1.2 风力机坐标系统建模

根据动量叶素理论得到的风力机叶片所受的气动力与力矩是相当于叶素局部坐标系统的。由于风力机气动性能最优设计的要求，一般情况下，风轮转轴倾角，桨叶锥角不为零，这样使得气动力／力矩需要完成从叶素局部坐标到全局坐标的坐标变换。此外，叶片旋转过程中，工作方位角对坐标系统转换关系的影响也需要考虑。文献［2］中对该部分内容进行了详尽的讨论。

大型风力机的叶片很长，桨叶不同位置的风速相差很大，一般情况下需要考虑风剪切的影响。本文采用的是指数模型［4，6］：

其中，H代表叶素高度；HR为参考高度；UR为参考高度处的风速；ε为经验风剪切指数，本文计算模型取值为0.1667。

1.3 计算流程

利用动量叶素理论获得叶片气动力主要分为速度诱导因子迭代求解及叶片气动性能计算两部分内容，文献［7］中详尽了叙述了速度诱导因子迭代的求解过程，这里由于篇幅限制，不再累述。程序的总流程如图1所示。

图1 风力机结构响应求解流程图Fig.1 Structure response solution scheme

1.4 风力机叶片气动力模型结果分析

南京航空航天大学针对NH1500叶片采用GH Bladed软件进行了叶片气动性能的校核。NH1500叶片为南航自主设计，长度40.5m，截面翼型沿叶片展向分别为 NH02＿40、NH02＿35、NH02＿30、NH02＿25、NH02＿21、NH02＿18、NH02＿15，翼型形状及气动性能数据为南航提供。针对该叶片，通过编写程序，完成了叶片气动力的建模工作，对比结果显示自编程序的计算结果与南航提供的采用GH Bladed软件计算的结果吻合的很好，验证了自编程序的正确性，计算结果如图2所示。

图2 功率系数计算结果对比图Fig.2 Power coefficient contrast

图2中纵坐标表示功率系数，横坐标表示来流风速。CpCAAA为本课题组自编程序得到的功率系数数据，CpNH为南航采用GH Bladed软件计算得到的功率系数数据。

2 风力机全机有限元建模

2.1 叶片结构有限元建模

风力机整体结构在气动力作用下的响应问题是风力机气动弹性研究的一个重要方向。结构设计需要根据响应结果获得弯矩载荷等重要数据。

风力机叶片外形沿叶片展向变化且内部结构复杂，进行结构三维有限元建模难度很大，使得分布气动力作用下的结构响应计算的代价极大。对响应问题，局部变形对整体响应结果的影响通常不予考虑。借鉴航空气动弹性问题的做法，将大展弦比机翼结构简化为梁模型处理，这里将大型风力机叶片也做简化处理。利用南航提供的结构刚度及质量分布数据通过模态分析可以获得叶片的结构动力学特性。叶片简单梁模型与三维有限元模型的结构动力学特性对比如表1所示。

对比结果显示，除了第四阶固有振动的频率误差略大之外，其它各阶频率误差都在有限元模型建模允许误差范围之内。

2.2 塔架结构有限元建模

进行风力机全机的响应计算，柔性塔架的影响需要考虑。目前，南航设计的大型风力机只提供了叶片的结构属性，无塔架属性可供参考。本部分内容主要根据塔架刚度设计准则［8］，对塔架构成部件及边界条件进行简化，建立了塔架的有限元模型。

表1 叶片简单梁模型与三维有限元模型的结构动力学特性对比图Table1 Structural dynamic characteristic of beam and 3DFE model

塔架结构主要包括塔体及塔顶集中质量。根据刚度设计准则，塔架结构的一阶固有频率以及其上10%，和三阶固有频率以及其下10%，需要避开风力机工作频率的1倍频和3倍频。南航设计的风力机额定转速17.2RPM／min，即工作频率为0.287Hz，三倍频为0.861Hz。根据刚度设计准则，塔架固有振动频率需要避开的频率范围为：0.287Hz～0.316Hz，0.775Hz～0.861Hz。

本文关于塔架的有限元建模没有考虑法兰上的螺栓、塔架内部的附属结构等构件；在塔架顶端创建一个质量单元节点，模拟塔顶质量，质量点的位置设置在机舱、轮毂和叶片的合重心位置处；此外，塔架基础按固支边界条件处理。

按照以上条件，塔架的结构动力学特性如表2所示。

表2 塔架的结构动力学特性Table2 Structural dynamic characteristic of tower

2.3 风力机全机建模

通过商用有限元软件建立风力机全机模型，叶片结构采用梁单元模拟，塔顶质量采用点质量单元模拟。叶片工作方位角每转过10°动量叶素气动力求解模块求解一次气动力并保存结果到数据文件。结构分析软件读入气动力数据，气动载荷以数据表形式施加到结构单元上，含边界条件的有限元模型如图3所示。

图3 风力机全机有限元模型示意图Fig.3 FE model of total wind turbine

3 计算结果与讨论

采用软件的瞬态动力学分析模块，结构动力学方程求解采用直接积分法中的Newmark方法，结构阻尼采用Rayleigh阻尼形式。Rayleigh阻尼中的α，β为不依赖于频率的常数，且阻尼对稳态振动的振幅没有影响。

考虑柔性塔架作用时，参考高度处来流风速为12m／s时，叶片叶尖节点振动位移随时间的变化如图4－图5所示。由于初始条件设置为零位移、零速度，振动的初始阶段会产生跳跃，是强迫振动的过渡阶段，随着时间的进行，振动逐渐进入简谐激励下强迫振动的稳态阶段。

3.1 挥舞方向振动

叶片X方向及Z（挥舞）方向振动的时程如图4所示。X方向位移最大值27.27cm，最小值23.91cm，平衡位置为25.59cm，振幅1.68cm。挥舞方向位移最大值20.07cm，最小值17.77cm，平衡位置为18.92cm，振幅1.15cm。

图4 叶片叶尖X及Z方向振动位移时程图Fig.4 Response of blade in Xand Zdirection

3.2 摆振方向振动

Y（摆振）方向位移最大值4.51m，最小值4.14m，平衡位置为4.33m，振幅18.85cm，时程如图5所示。

图5 叶片叶尖Y方向振动位移时程图Fig.5 Response of blade in Y－direction

3.3 塔架振动

塔架的振动产生结构内部载荷，塔架的振动情况将影响结构的弯矩载荷及疲劳设计。图6给出塔架顶端节点的时域振动情况。不考虑塔架对叶片的影响，如叶片经过塔架时加速逆风效应时，塔架顶端迎风方向位移最大值7.37cm，最小值6.97cm，平衡位置为7.17cm，振幅0.2cm，时程如图6所示。

图6 塔架顶端振动位移时程图Fig.6 Response of tower in time domain

3.4 几何非线性

大型风力机叶片长度非常大，受力时会发生大位移小应变的几何非线性效应，这时平衡方程和几何关系都是非线性的。早期的非线性有限元分析基本上是线性分析的扩展，只能针对个别具体问题进行分析。近年来基于非线性连续介质力学原理的有限元分析有了很大的发展，商用软件也有相应成熟的模块求解该类问题。本部分利用商用有限元软件分别进行了叶片线性响应分析与几何非线性响应分析，对比结果如图7和表3所示。

结果显示，对叶片叶尖节点，考虑几何非线性时叶片的变形较小，呈现出比线性分析结果刚硬的性质，这是叶片初位移效应及初应力效应共同作用的结果。

图7 几何非线性对叶片振动的影响Fig.7 Influence of geometry nonlinearity

表3 几何非线性对叶片振动的影响Table3 Influence of geometry nonlinearity

4 结 论

本文建立了一套适用于大型风力机气动弹性响应的快速计算方法。该方法可以对叶片所受气动力进行快速预测，并快速的给出结构在气动力作用下的响应情况，提供了一种可供风力机设计人员借鉴的技术手段。

［1］ROBERT E WILSON，PETER B SLISSAMAN，STEL N WALKER.Aerodynamic performance of wind turbines［M］.Department of Mechanical Engineering，Oregon State University，1976.

［2］刘雄，陈严，叶枝全.水平轴风力机气动性能计算模型［J］.太阳能学报，2005，26（6）：792－800.

［3］PRANDTL L，TIETJENS O G.Applied hydro and aeromechanics［M］.Dover Publications，1957.

［4］DET NORSKE VERITAS and RIS NATIONAL LABORATORY.Guidelines for design of wind turbines［M］.Det Norske Veritas and Ris National Laboratory，2002.

［5］GLAUERT H.The analysis of experimental results in the windmill brake and vortex ring states of an airscrew［R］.Reports and Memoranda，1926.

［6］SPERA D A.Wind turbine technology［M］.New York ASME Press，1994.

［7］ROBERTS E.WILSON，PETER B.S.LISSAMAN.Applied aerodynamics of wind power machines［M］.Or－egon State University Corvallis，1974.

［8］赵立新.风力发电机塔架的有限元分析与优化设计［D］.吉林：吉林大学，2008.