水平轴风力机叶片气弹建模与响应分析

陆 洋， 周桂林

（南京航空航天大学 旋翼动力学国防科技重点实验室，江苏 南京210016）

0 引 言

大型水平轴风力机系统是强非线性流－刚柔耦合的周期时变多体系统，结构和运动关系非常复杂，对风力机叶片进行动力学建模时必须考虑叶片的几何非线性以及非定常气动载荷等因素，理论推导和数值计算都比较困难。另一方面，随着风力机额定功率不断增加，叶片尺寸越来越大，细长叶片在交变载荷下的挥舞变形更大。为防止叶尖与塔架碰撞，出现了预弯叶片设计概念，增加了风力机叶片气弹建模的难度。

风力机叶片气弹建模包括叶片结构动力学建模和空气动力学建模。文献［1］采用模态分析法进行叶片结构动力学建模，采用Greenberg's二维准定常气动力模型计算气动力，研究了风力机气弹稳定性；文献［2］采用5节点18自由度梁单元对叶片进行有限元离散，采用二维准定常气动模型，利用拟线性法计算了水平轴风机工作状态下的气动弹性响应；文献［3］以Marc有限元软件为基础平台，采用动量叶素理论二维气动力模型，计算叶片的气弹响应；文献［4］采用二节点梁单元对叶片进行有限元离散，采用Leishman－Beddoes非定常模型计算气动力，分析了叶片的动态响应。

可见，为了保证叶片响应的计算精度，有限元方法已逐渐替代模态法成为风力机叶片结构动力分析的主流方法，气动力建模则仍以动量叶素理论和二维准定常气动力模型为主。

然而，现有的叶片气弹建模方法仍存在以下局限性：1）基于有限元方法的风力机动力学方程通常建立在非惯性系下，即处理动能部分时将其分为旋转坐标系的动能和牵连运动引起的动能两部分，破坏了推导和编程的一致性；2）建模时，没有考虑叶片预弯对叶片气动弹性动力响应的影响；3）气动力模型多采用二维准定常空气动力计算模型，没有考虑叶片动态入流以及动态失速影响。这些局限性，使得传统建模方法的计算精度和适用性受到了制约。

针对前述叶片建模方法的局限性，本文开展了如下工作：1）引入了Hartenberg－Denavit增广转换矩阵［5］，通过递归计算的方式得到惯性坐标系下的位置矢量，从而在惯性系下建立了风力机叶片动力学方程，将动能部分合为一体，避免了传统建模中的繁琐；2）引入上反坐标系方法描述风力机叶片的预弯变形，考虑了预弯变形影响；3）入流模型采用Suzuki提出的针对风机叶片的广义动态尾迹入流模型（GDW）［6］，并采用Pierce修正的风机Leishman－Beddoes模型［7］计算翼型非定常气动力，综合考虑了叶片动态入流以及动态失速影响。利用Newmark数值积分方法求解叶片动力学方程获得稳态周期解。最后，分别以美国可再生能源实验室（NREL）Phase VI非定常空气动力学实验风力机［8－9］及其公开的1.5MW风力机叶片［10］为算例计算了有／无预弯叶片的气弹响应，验证了本文所建模型的正确性和有效性。

1 系统坐标系

1.1 增广转换矩阵

为描述风力机叶片复杂的结构形式及运动特点，引入多个坐标系，表1给出了建立风力机系统的坐标系定义和描述。图1为风力机系统坐标系示意图。对于开环拓扑结构，如果要描述梁剖面上的任意一点P，在变形后的坐标系下，可以将位置矢量表示为：

将位置坐标行向量S增广为：

在未变形坐标下，增广后P点的位置增广向量及增广转换矩阵为：

根据式（3）类推，在n个坐标系中第i个坐标系下，P点的位置增广向量及增广转换矩阵为：

依次类推，在惯性坐标系下，叶片上任意点P的位置增广向量表示为：

由上面的递归表达式（4）可以看出，只要确定了各坐标系原点的位置和转换关系，就可以通过递归计算的方式得到惯性坐标系下的位置矢量，从而简化了刚柔耦合计算，使推导和编程也得到简化，提高了算法对风力机结构和运动形式的适用性。

图1 风力机系统坐标系示意图Fig.1 Schematic of coordinate system of wind turbine

表1 风力机坐标系统Table 1 Coordinate system of wind turbine

1.2 上反坐标系

叶片设计时预先向偏离塔架方向弯曲，可以增大叶尖与塔架之间的净空间隙，减少二者发生碰撞的可能性。大型风力机叶片长度通常超过30m，叶片一般从1／3长度位置处开始预弯。对于长度为35m左右的1.5MW风机叶片，预弯高度一般取值为600～750mm。由于预弯值远小于叶片长度，为了描述预弯部分叶片的位置矢量，本文提出一种简化处理办法，如图2所示，用直线叶片BC′来替代连续预弯叶片部分BB′C′，直线叶片BC′的上反角度为ζ0。因此，为了描述整个叶片的位置矢量，需要将叶片划分为AB和BC′两段进行处理，在BC′段引入上反坐标系Rs：osxsyszs。

图2 预弯叶片示意图Fig.2 Sketch of pre－curved rotor blade

上反坐标系与变距坐标系之间的转换关系Ts为：

对于预弯量较大的叶片，可以采用类似方法，引入多个上反坐标系。引入的上反坐标系越多，则预弯变形情况描述越精确，综合考虑计算效率，可选择适当数量的上反坐标系进行计算。

2 叶片动力学建模

2.1 叶片单元弹性势能贡献项

采用中等变形梁理论得到Green应变表达式［11］，带入单元弹性势能表达式并对其求变分，得到弹性广义力及弹性势能贡献项：

2.2 叶片单元动能贡献项

单片叶片上的动能变分［12］可以表示为：

经推导得动能产生的第i个广义力为：

动能所产生的切线质量、阻尼和刚度阵为：

为了计算质量、阻尼和刚度阵，还需求出增广位置矢量的相关导数和偏导数，利用递推算法后，可使推导通用性大大提高，避免了传统建模中的繁琐。

2.3 重力功贡献项

在计算风力机风轮叶片气弹响应时，必须考虑叶片的重力影响，叶片重力引起的外力功和重力广义力分别为：

2.4 气动力贡献项

诱导速度计算采用Suzuki广义动态尾迹模型［6］。用任意谐波次数和任意阶次径向型函数的级数形式来描述风轮平面垂向诱导速度：

其中，Δ为环量项，Δ为非环量项。利用Kirchhoff理论处理翼型后缘分离的非线性问题，翼型法向力系数和弦向力系数与分离点的关系为：

式中、ψ分别是无量纲半径、方位角和无量纲时间。是径向函数，是基于时间的状态量。

翼型剖面非定常气动力计算采用Pierce修正后应用于风力机计算的Leishman－Beddoes模型［7］。该模型通过阶跃响应的叠加模拟附着流条件下的非定常效应：

2.5 叶片动力学方程

应用Hamilton原理建立系统运动方程：

其中，CN是法向力系数，CC是弦向力系数，f是分离点位置。根据动态失速的特性，模拟动态失速涡沿翼型上表面的移动及因而引起的压力中心的移动，从而计算动态失速涡诱导产生的升力和俯仰力矩［9］。空气动力FA及力矩MA所做虚功的变分可以表示为：

对式（16）进行变分运算和有限元离散，组集动能、应变能及空气动力等产生的虚功，得到基于广义力形式的叶片非线性动力学隐式微分方程：

式中右上标T、E、A、G分别表示动能、应变能、气动力及重力引起的广义力，q、、为广义位移、速度、加速度向量。

由于建模过程中考虑了叶片整体的刚性运动与叶片自身弹性运动间的刚柔耦合特点及非定常和动态失速等非线性因素，最终采用适应性较强的Newmark数值积分方法进行运动方程求解。

3 算例验证

算例1：对NREL Phase VI非定常空气动力学实验中叶片长度为10.058m的小型水平轴风力机进行计算并与Bladed软件以及实验结果比较。本文分别计算了来流速度分别为 7m／s、10m／s、15m／s、20m／s、25m／s，偏航角为0°时，下风向风机径向位置为63%截面的法向力系数Cn，切向力系数Ct，俯仰力矩系数Cm以及叶片根部的挥舞力矩。图3～图6为本文计算结果与实验数据及Bladed计算结果的对比曲线。

从图3～图5可以看出，相比较Bladed计算结果，本文计算得到的截面气动力系数与实验结果吻合程度更佳，表明了本文叶片气弹建模中空气动力学建模方法的正确性和精确性。图6中计算结果与实验值的叶片根部挥舞弯矩吻合程度较高，而叶片根部挥舞力矩主要由气动力贡献，再次证明了本文气动力模型的正确性。

算例2：以NREL 1.5MW下风向风机叶片为算例进行气弹响应计算。该风机基本参数如下：风轮仰角η＝0°，风轮锥角χ＝0°，叶片长度33.25m，轮毂半径1.75m，来流速度10m／s，风轮转速20rpm。叶片无预弯，剖面特性具体参数见文献［8］。为了分析叶片预弯对叶片气弹响应的影响，假设该叶片在1／3位置处开始连续预弯，至叶尖处预弯高度为3.5m，即上反角10°，计算该预弯叶片的气弹响应。将上述计算结果与NREL公开文献的计算结果一并进行对比分析，以验证本文算法的正确性。图7～图12给出了叶片旋转系下单片叶片作用于轮毂三个方向的载荷及力矩沿方位角变化对比曲线。

图7 叶片根部径向载荷Fig.7 Blade axial force at blade root

对比无预弯情况下叶片在旋转系下的载荷和力矩，本文计算结果与NREL公开文献的计算结果的变化趋势以及幅值符合较好，其中叶片径向载荷的相位存在约50°的差异。由于叶片径向载荷主要由风轮旋转产生的离心力和周期变化的重力载荷分量构成，考虑重力载荷在叶片径向周期变化，可知本文计算结果更加合理。通过无预弯情况的算例对比，验证了本文建模算法的正确性。

图12 叶片根部摆振力矩Fig12 Blade in－plane moment at blade root

对比有／无预弯两种情况下基于本文模型计算所得的叶片载荷和力矩，可以看出对于下风向风机预弯对叶片的径向载荷以及摆振方向的载荷影响不大（图7、图8），这是因为叶片摆振方向的载荷主要是周期变化的重力载荷，预弯对重力载荷基本没有影响。同理，由于叶片摆振力矩主要为摆振载荷贡献，预弯对摆振力矩影响也很小（图12）。由图9和图11可知，预弯使得叶片挥舞方向载荷、力矩的稳态值变小。挥舞方向载荷的稳态值主要为气动力，表明预弯降低了叶片气动力载荷。由图10叶片变距力矩对比曲线可以看出，预弯使得叶片变距力矩的绝对值增大。影响叶片变距力矩的因素很多，其中最重要的是由叶片自身质量引起的离心力作用而产生的惯性力矩和变距结构各种摩擦副产生的摩擦阻力距。下风向风机叶片由于预弯增加了叶片离心力作用产生的惯性力矩，从而使叶片变距力矩显著增加。通过上述分析可知，本文采用上反坐标系处理叶片预弯变形的方法所得到的预弯影响与理论分析相符，表明该处理方法是合理的。

4 结 论

本文基于Hamilton原理，建立了水平轴风力机叶片的气弹模型。通过数学建模和算例分析，可以得到以下结论：1）引入增广转换矩阵，在惯性系下建立叶片动力学方程，可以简化刚柔耦合计算，避免传统建模中的繁琐；2）气动力建模时，采用Suzuki广义动态尾迹方法计算诱导速度，采用Pierce修正的风机Leishman－Beddoes模型计算翼型气动力，可提高了气动力计算精度，比Bladed计算结果更符合实验值；3）本文所建模型预测的无预弯叶片的载荷与NREL计算结果符合较好；4）采用上反坐标系处理叶片预弯变形，计算所得的预弯影响与理论分析相符，验证了本文所建立的可考虑叶片预弯的气弹动力学模型的合理性。

［1］李本立，安玉华.风力机气动弹性稳定性研究［J］.太阳能学报，1996，17（4）：314－320.

［2］王介龙，陈彦，薛克宗.风力发电机耦合转子／机舱／塔架的气弹响应［J］.清华大学学报（自然科学版），2002，42（2）：211－215.

［3］傅程，王延荣.风力发电机叶片气动弹性响应分析［J］.机械设计与研究，2009，25（1）：68－71.

［4］LIU X，ZHANG X M.Dynamic response analysis of the rotating blade of horizontal axis wind turbine［J］.WindEngineering，2010，34（5）：543－559.

［5］GERADIN M，CARDONA A.Kinematics and dynamics of rigid and flexible mechanisms using finite elements and quaternion algebra［J］.ComputationalMechanics，1989，18（4）：651－659.

［6］SUZUKI.Application of dynamic inflow theory to wind turbine rotors［D］.［Doctoral Dissertation］.Salt Lake City：Department of Mechanical Engineering，University of Utah，2000.

［7］KIRK GEE PIERCE.Wind turbine load prediction using the Beddoes－Leishman model for unsteady aerodynamics and dynamic stall［D］.［Doctoral Dissertation］.Salt Lake City：Department of Mechanical Engineering，University of Utah，1996.

［8］HAND M M，SIMMS D A.Unsteady aerodynamics experiment phase VI：wind tunnel test configurations and available data campaigns［R］.NREL／TP－500－29955，2001.

［9］HAND M M，SIMMS D A.NREL unsteady aerodynamics experiment in the NASA－Ames Wind Tunnel：A comparison of predictions to measurements［R］.NREL／TP－500－29494，2001.

［10］JASON M.FAST user's guide［R］.Technical Report NREL／EL－500－38230，2005.

［11］王浩文，高正，郑兆昌.前飞状态下直升机旋翼系统气弹响应及稳定性分析［J］.振动工程学报，1999，12（4）：521－528.

［12］王益锋.直升机旋翼桨叶动力学建模与弹性碰撞动响应分析［D］.［博士学位论文］.南京：南京航空航天大学，2009.