李德源， 池志强， 汪显能， 倪晨峰
(广东工业大学 机电工程学院， 广州 510006)
摘 要：为了研究大型风力机叶片在静止和转动状态下的振动模态及其变化特点，通过叶素动量理论和复合材料的叶片设计方法完成了10 MW风力机叶片的设计.基于多体系统动力学理论和超级单元模型，结合动力学分析软件ADAMS对静止状态下叶片的线性特征值进行了分析，考虑叶片的弹性变形和旋转，应用刚性积分方法对叶片的非线性控制方程进行数值求解，通过傅里叶谱分析方法，实现风轮旋转条件下的运转模态识别.结果表明，在动力刚化效应作用下，叶片的固有频率会随着转速的增加而增大.
关 键 词：10 MW风力机； 柔性叶片； 超级单元； 多体动力学； 固有频率； 振型； 频谱分析； 动力刚化
随着风电技术的发展，现代风电机组的单机功率不断增大，投入商业运行的主流产品以1.5～6.0 MW的机型为主，而10 MW的巨型风电机组也已成为世界各大风电公司和研究机构竞相研发的对象 [1].而大型风力机柔性叶片在运行过程中所产生的较大变形对其动力特性和气弹稳定性的影响日益突出 [2]，准确模拟柔性叶片在静止与风轮转动过程中的各阶模态并理清其变化规律成为备受关注的研究课题.
叶片气动外形设计主要是选择合适的翼型并确定弦长、安装角以及厚度等气动参数的分布 [3].本文在参考文献[4]所设计的5 MW风力机模型和参考文献[1]中的10 MW概念型风力机叶片基础上，基于叶素动量理论，考虑叶尖损失等因素的影响，结合SolidWorks三维建模软件设计了机组功率为10 MW的风力机叶片的气动外形和内部结构.
目前风力机叶片的模态分析一般采用两大类方法 [5]：一类是直接的特征值分析，这需要在某个平衡状态对非线性方程进行线性化，建立特征方程，如ADAMS等动力学分析软件的模态分析模块；另一类是用一组已知的力函数激励叶片，基于时域响应信号，用系统模态识别的方法提取出叶片的模态信息，这对于叶片转动情况下的模态识别是可行的.
考虑大型风力机柔性叶片变形的不可忽略性，本文采用超级单元法 [6]，将柔性叶片离散为由运动副、弹簧和阻尼器连接的多体系统，在多体动力学分析软件ADAMS中建立了叶片的多体动力学模型，应用ADAMS的模态分析模块和动力学分析与后处理模块分别对所设计的10 MW风力机叶片在静止和各转速下的固有特性进行了分析.
设计的10 MW风力机叶片是对美国可再生能源实验室(NREL)的5 MW风力机模型 [4]进行相应的放大，参考了DTU(丹麦技术大学风能实验室)所设计的10 MW风力机模型 [1]，并对部分参数进行调整和修正，应用SolidWorks三维建模软件进行设计.其中，机组额定功率10 MW，切入风速4 m/s，切出风速25 m/s，额定风速11.4 m/s，额定转速9.6 r/min，叶片数量为3.
1.1 10 MW风力机叶片气动外形设计
叶片长度根据输出功率、风轮直径和来流风速之间的关系式确定，即
(1)
式中： P为机组输出功率，取为10 MW； ρ为空气密度，取值1.225 kg/m 3； V r为额定风速，取值11.4 m/s； D为风轮直径； C p 为风能利用系数，取值0.47； η 1为传动系统效率，取值0.97； η 2为发电机效率，取值0.96.由式(1)可得风轮直径约为179 m，取轮毂直径5.6 m，取叶片长度为86.3 m.
叶片气动设计的方法通常以叶素动量理论为基础，同时引入Prandt叶尖和轮毂损失修正，利用叶栅理论对攻角的修正以及Glauert对轴向诱导因子的修正，建立风电机组气动设计和计算的数学模型.在选择沿展向所采用的翼型时，考虑到叶根部分与轮毂连接的结构要求以及叶中到叶尖部分的功率输出要求，从最大弦长附近开始选择了具有一定尾缘厚度的FFA-W3系列5种厚度的翼型，分别是FFA-W3-241、FFA-W3-301、FFA-W3-360、FFA-W3-480和FFA-W3-600，靠近叶尖部分选择了最大相对厚度为24.1%.由于大型风力机叶片一般具有较大的雷诺数，因此采用了参考文献[1]中修正过的翼型数据.为了保证叶片表面的光顺性，对弦长、扭角和相对厚度的设计结果进行了多项式拟合，最终得到的弦长、扭角和相对厚度的分布如图1～3所示( r为叶片沿展向位置与叶片总长比).
图1 弦长分布
Fig.1 Distribution of chord
图2 扭角分布
Fig.2 Distribution of twist angle
图3 相对厚度分布
Fig.3 Distribution of relative thickness
1.2 叶片的铺层结构设计
叶片铺层设计的任务就是根据铺层的目标力学性能，确定铺层中纤维布的类型、铺设方向、铺设次序、铺设层数等.铺设方向一般为0°、90°和±45°方向，通常0°铺层承受轴向荷载，±45°铺层承受剪切荷载，90°铺层承受横向荷载和控制泊松比.本文设计的叶片铺层中，主梁的结构形式采用了矩形梁，如图4所示.采用玻璃钢材料，最终得到的质量、截面刚度和DTU模型基本一致，分布如图5～7所示，叶片总质量为41 724.8 kg.
图4 叶片内部结构
Fig.4 Internal structure of blade
图5 质量分布
Fig.5 Distribution of mass
图6 挥舞刚度分布
Fig.6 Distribution of flapwise stiffness
图7 摆振刚度分布
Fig.7 Distribution of edgewise stiffness
为了描述叶片较大的弹性变形，在多体系统模型中引入“超级单元”对叶片离散，使其成为由刚体、铰、力元和外力组成的多体系统 [6].借助机械系统动力学分析软件ADAMS来建立柔性叶片的动力学模型.
水平轴风力机的叶片固定安装在轮毂上，力学特性类似于悬臂梁，将其离散成为若干个超级单元，每个超级单元由4个刚体组成，每两个刚体间采用万向节或转动铰连接，其结构如图8所示.两个相互垂直平面的弯曲由万向节定义，弹簧和阻尼器约束刚体的相对运动.
图8 超级单元
Fig.8 Super element
根据solidworks中建立的叶片三维模型，将叶片离散为6个超级单元，前一个单元最后一个刚体与下一个单元第一个刚体刚性连接，合并为一个刚体，因此叶片共分割为19个刚体，共31个自由度，模型如图9所示.
图9 叶片的超级单元模型
Fig.9 Super element model for blade
根据梁的弯曲与扭转理论，计算出各弹簧的刚度系数 [6]，计算公式为
(2)
(3)
(4)
(5)
(6)
式中： L s为超级单元长度； E为综合弹性模量； G为综合剪切模量； I为截面惯性矩； C x1 为扭转弹簧系数，其余为弯曲刚度系数.根据Rauh的结论 [7]，当1/5< k<1/4时，能够用为数不多的超级单元得到固有频率的近似解，本文 k取0.217 6.
据式(2)～(6)及上节叶片截面数据计算叶片上各个刚体之间的弹簧刚度系数如表1所示.
表1 10 MW风力机叶片弹簧刚度系数
Tab.1 Spring stiffness coefficients of blade of 10 MW wind turbine
3.1 叶片动力学方程与求解
在ADAMS中，用某刚体的质心笛卡尔坐标和反映刚体方位的欧拉角作为广义坐标，由于采用了不独立的广义坐标，系统动力学方程虽然数量较大，但却是高度稀疏耦合的微分代数方程，可表示为
q=[ x， y， z， ψ， θ， φ] T
(7)
考虑约束方程，ADAMS利用带乘子的第一类拉格朗日方程的能量形式得到计算式为
(8)
式中： T为系统广义坐标所表达的系统动能； q j 为广义坐标； Q j 为与广义坐标 q j 对应的广义力；最右端项涉及约束方程和拉格朗日乘子 λ i 表达式在广义坐标 q j 方向的约束反力.
式(8)所建立的系统运动微分方程一般是一组关于
N个广义坐标
q
j
的二阶非线性微分
代数方程组(DAE)，结合相应的初始条件和数值积分法，求出方程组的数值解，这些数值解就代表了系统的运动规律.
3.2 采用ADAMS/Vibration模块的模态计算
在ADAMS中建立10 MW风力机叶片的多体动力学模型后，将叶片根部第一个刚体与惯性参考系之间采用固定副进行约束.使用ADAMS/Vibration振动分析模块对静止状态下的叶片进行模态计算，即根据线性化动力学方程进行特征值分析，得出的叶片振型和前若干阶频率如图10、11和表2所示.该风力机风轮额定转速为9.6 r/min，即转频为0.16 Hz，可知叶片各阶频率均不会引起共振.
图10 叶片第1阶振型
Fig.10 First order vibration mode of blade
图11 叶片第2阶振型
Fig.11 Second order vibration mode of blade
表2 本文模型与DTU模型对比
Tab.2 Comparison between proposed
model and DTU model
3.3 时
频域分析方法确定固有频率
确定叶片的固有频率还可以对叶片多体系统施加特定的外激励，通过求解控制方程，并对响应进行频谱分析，由系统的传递函数特性来获得.
对于一个机械系统，其传递函数定义为输出 x( t)与输入 f( t)的傅里叶变换之比，即
(9)
式中：
；
t.
一般外激励(输入函数 f( t))采用三角波脉冲载荷，当脉冲激励时间足够短时完全可以激励出叶片前若干阶频率，对其进行傅里叶变换，得到频域上的输入函数 F i .通过求解式(8)，可得到叶片系统在该脉冲激励下的动力学响应，如各刚体速度、加速度响应及位移响应.对响应进行傅里叶变换，得到频域上的输出函数 X i ，按式(9)可求得各刚体的传递函数 H i .
在叶片刚体的挥舞和摆振方向上同时施加三角波脉冲载荷，数值求解出叶片各刚体的动力响应，图12为叶尖挥舞和摆振方向的位移响应.
图12 叶尖挥舞位移和摆振位移
Fig.12 Flapwise and edgewise displacements
of blade tip
通过ADAMS处理模块对响应进行快速傅里叶变换，并结合式(9)可得出系统各低阶频率，如图13所示.其中前两阶挥舞为0.62 Hz和1.77 Hz，前两阶摆振为0.96 Hz和2.83 Hz，与表2中结果符合.
图13 叶片固有频率
Fig.13 Natural frequencies of blade
风力机运转时，由于柔性叶片存在弹性变形及振动，各广义坐标和速度间存在耦合，叶片的转动和其弹性变形是耦合的，这种耦合会产生所谓的动力刚化或弱化效应 [8].同时，从系统的非线性微分方程要导出相应的特征方程有一定困难，各叶片在不同时刻处于不同位置，即系统的质量与刚度分布是时变的，此时求系统的模态一般比较困难.
本文通过在旋转叶片上施加激励，应用前述时
频域分析方法来确定叶片在各转速下的频率.通过将叶片的旋转速度引入叶片多体系统约束方程，利用固定在叶根处的坐标系来观察叶片的响应，由频域分析方法获得叶片的固有频率随转速的变化情况.
分析时，对叶尖处的刚体施加角位移驱动函数，使系统从静止状态下稳步旋转到相应的转速，函数表达式 [9-10]为
(10)
式中： φ为角位移； ω 1为输入角速度； t为时间(0～50 s)，求解步长控制为0.001 s.
通过求解非线性动力方程分别得到转速为3、6、9.6、12、15 r/min的时域响应，选取叶尖的响应信号进行频谱分析，得到各阶模态所对应的频率如表3所示，叶片频率随转速变化的根迹图如图14所示.
表3 不同转速下的固有频率
Tab.3 Natural frequencies under
different rotational speeds
图14 叶片在旋转状态下的固有频率
Fig.14 Natural frequencies of blade
under rotational condition
从表3和图14中可以看出，随着转速的不断增加，叶片的动力刚化现象会越来越明显.其中一阶、二阶挥舞频率分别增加了0.09 Hz和0.08 Hz，增加幅度为14.5%和4.5%；摆振方向分别增加了0.02 Hz和0.07 Hz，增加幅度为2%和2.5%.摆振方向的刚化效应没有挥舞方向明显，原因在于叶片转动过程中，旋转平面内沿叶片轴向的哥氏惯性力弱化了摆振方向的刚度，即部分抵消了摆振方向的动力刚化效应.
应用风力机空气动力学理论和设计理论，基于SolidWorks三维建模软件设计了10 MW风力机叶片的气动外形和内部铺层结构；通过超级单元法将柔性叶片离散为由运动副、力元弹簧和阻尼器连接的有限个刚体，利用ADAMS建立叶片的多体动力学模型.分别通过ADAMS/Vibration模块和频域分析完成了叶片在静止和转动状态下固有频率的求解.考察了基于叶片旋转的动力刚化效应，模拟了叶片频率随转速的变化趋势，挥舞方向的变化会大于摆振方向的变化，分析模型为下一阶段研究叶片的气弹稳定性 [11]打下了基础.
参考文献(References)：
[1]Christian B，Frederik Z，Robert B，et al.Description of the DTU 10 MW reference wind turbine [R].Denmark：DTU Wind Energy Laboratory，2013.
[2]胡燕平，戴巨川，刘德顺.大型风力机叶片研究现状与发展趋势 [J].机械工程学报，2013，49(20)：140-146.
(HU Yan-ping，DAI Ju-chuan，LIU De-shun.Research status and development trend on large scale wind turbine blades [J].Journal of Mechanical Engineering，2013，49(20)：140-146.)
[3]李春，叶舟，高伟.现代陆海风力机计算与仿真 [M].上海：上海科学技术出版社，2012：108-133.
(LI Chun，YE Zhou，GAO Wei.Calculation and simulation of modern offshore and onshore wind turbine [M].Shanghai：Shanghai Science and Technology Press，2012：108-133.)
[4]Jonkman J，Butterfield S，Muaial W，et al.Definition of a 5-MW reference wind turbine for offshore system development [R].Springfield：US National Renewable Energy Laboratory，2009.
[5]Bir G，Jonkman J.Aeroelastic instabilities of large offshore and onshore wind turbines：preprint [J].Journal of Physics：Conference Series，2007，75(2)：12069-12087.
[6]Meng F Z，Ozbek，M，Rixen D J，et al.Aero-elastic stability analysis for large-scale wind turbines [D].Delft：Delft University of Technology，2011：339-349.
[7]李德源，莫文威，严修红，等.基于多体模型的水平轴风力机气弹耦合分析 [J].机械工程学报，2014，50(12)：140-143.
(LI De-yuan，MO Wen-wei，YAN Xiu-hong，et al.Aeroelastic analysis of horizontal axis wind turbine based on multi-body model [J].Journal of Mechanical Engineering，2014，50(12)：140-143.)
[8]潘柏松，陈栋栋，韩斌，等.基于不同厚度兆瓦级风力机叶片动态特性研究 [J].浙江工业大学学报，2013，41(4)：364-366.
(PAN Bo-song，CHEN Dong-dong，HAN Bin，et al.Study of dynamic characteristics based on different thickness for MW wind turbine blade [J].Journal of Zhejiang University of Technology，2013，41(4)：364-366.)
[9]Gebhardt C G，Roccia B A.Non-linear aeroelasticity：an pproach to compute the response of three-blade large-scale horizontal-axis wind turbines [J].Renewable Energy，2014(66)：495-514.
[10]李德源，严小辉，莫文威.基于脉冲激励响应的风力机叶片结构阻尼计算 [J].沈阳工业大学学报，2014，36(6)：619-624.
(LI De-yuan，YAN Xiao-hui，MO Wen-wei.Structural damping calculation of wind turbine blade based on pulse excitation response [J].Journal of Shenyang University of Technology，2014，36(6)：619-624.)
[11]Gebhardt C G，Preidikman S，Jgensen M H，et al.Non-linear aeroelastic behavior of large horizontal-axis wind turbines：a multibody system approach [J].Journal of Hydrogen Energy，2012，37：14719-14724.
(责任编辑：景 勇 英文审校：尹淑英)
LI De-yuan, CHI Zhi-qiang, WANG Xian-neng, NI Chen-feng
(Faculty of Electromechanical Engineering, Guangdong University of Technology, Guangzhou 510006, China)
Abstract：In order to study the vibration modes and change characteristics of blade of large wind turbine under stationary and rotational conditions, the blade design of 10 MW wind turbine was completed through blade element momentum theory and blade design method of composites. Based on the multi-body system dynamics theory and super element model, the linear eigenvalue of blade under stationary condition was analyzed in combination with the dynamics analysis software ADAMS. With the consideration of elastic deformation and rotation of blade, the nonlinear control equations of blade were numerically solved with the rigid integral method. Through Fourier spectrum analysis method, the operation mode recognition was accomplished under the rotational condition of wind wheel. The results indicate that under the effect of dynamic stiffening, the natural frequency of blade will increase with increasing the rotational speed.
Key words：10 MW wind turbine; flexible blade; super element; multi-body dynamics; natural frequency; vibration mode; spectral analysis; dynamic stiffening
收稿日期：2015-12-24.
基金项目：国家自然科学基金资助项目(51276043).
作者简介：李德源(1965-)，男，重庆人，教授，博士，主要从事风力机气动与结构分析等方面的研究.
doi：10.7688/j.issn.1000-1646.2016.03.01
中图分类号：TK 83
文献标志码：A
文章编号：1000-1646(2016)03-0241-06
*本文已于2016-03-02 16∶48在中国知网优先数字出版. 网络出版地址： http：∥www.cnki.net/kcms/detail/21.1189.T.20160302.1648.046.html
风力发电技术