曹 莉， 孙文磊， 周建星

(新疆大学 机械工程学院，乌鲁木齐 830047)

由于全球性能源短缺和环境污染等问题，风能作为一种清洁的可再生能源，变得越来越重要[1]。风力发电机组常年工作在变速、变载荷的工作环境中，在随机风的作用下会受到频繁的扰动和动载激励，对其工作性能和使用寿命有很大的影响[2]。随着风电机组不断向大型化、大功率方向发展，风力发电机组的动态特性更加复杂。因此，在风力机的设计过程中，需要考虑风力机在典型工况下的动态特性。

当前，在风力机动力学分析领域，常用的方法有多体动力学方法MBS (Multibody Systems)、有限元方法FES(Finite Element System)、模态分析方法及连续系统COS(Continuous Systems)等[3]，MBS多体动力学方法将实际的机械构件视为刚体来建立系统动力学微分方程(组)，该方法模拟的精度有限[4]；COS连续系统方法建立的偏微分方程组仅在特殊、简单的几何结构及载荷下才可以求解。在风力机整机动力学方面主要采用有限元方法和模态分析方法[5]。但有限元方法具有较多的自由度，计算和分析成本较高[6]；模态叠加法通过提取结构的主模态集(主振型阵)形成坐标变换并实现降价来缩短计算时间，相对于其它方法可使用较少的自由度对风力机动态行为进行可靠分析，具有计算速度较快、效率高等优点。

近年来，国内外学者对风力发电机组的动态特性问题进行了研究： Asareh等[7]以FAST软件为平台，结合气动力分析了5 MW风力机在不同工况下的动态响应；刘雄等[8]研究了风力机在湍流风作用下的动态特性，并揭示了叶片的离心刚化和气动阻尼对分析结果的影响较大；吕计男等[9]对叶片和塔架进行简化并建立了风力机整机模型，计算了风力机在启动工况下的动态响应及叶片几何非线性对结果的影响。李明等[10]采用车载法模拟自然风场，对风力机叶片的相似模型进行测试，得到的叶片各截面的载荷数据，与BLADED仿真结果较为一致。研究中风力机所处工况多为单一工况，关于风力机在不同工况下的动态特性研究较少；研究对象多为叶片、塔架、轮毂等单个部件，鲜有文献对风力机整机的动态特性研究。

因此，笔者针对复杂工况下的大型风力机，采用了基于模态叠加法的动态特性分析方法，建立了风力机动态特性计算模型，采用该方法对2 MW大型风力机进行分析，同时总结了不同工况下的风力机的动态特性。

1 风力机动态特性计算理论

1.1 风力机叶片载荷计算理论

风力机在运行时主要通过叶片捕获风能，捕获的风能带动风轮转动，再通过传动系统将风能转化为机械能，因此叶片是主要的受力部件[11]。叶片所受的力主要包括空气动力、离心力、和重力。叶片的受力和变形坐标系如图1所示，其中 为迎风向， 为叶片展长方向， 由右手规则确定。

图1 叶片受力和变形坐标系Fig.1 Coordinate system for blade loads and deflections

对叶片的载荷进行综合可以计算出叶片单位长度所受的载荷[12]：

(1)

式中：FXg，FYa为叶片单位长度上的气动载荷；FXg，FYg，FZg为叶片单位长度上的重力；FYc，FZc为叶片单位长度上的离心力。由于笔者所考虑的风力机风轮是静止的，气动载荷和离心力可忽略不计，因此迎风向载荷FXa是主要载荷。

1.2 基于模态叠加法的风力机动态特性计算理论模型

笔者将有限元法和模态叠加法这两种方法进行了结合，其基本思想是以有限元为平台，反复利用瞬态分析求解风力机在周期性载荷作用下的稳态响应。瞬态分析包括完全法、缩减法和模态叠加法。完全法计算量大且计算时间较长；缩减法不能在时间单元上添加载荷且所有荷载必须加载用户定义的主自由度上；模态叠加法相对于其它方法计算速度较快、效率高。因此本文选用模态叠加法对风力机整机动态特性进行求解，得到风力机在时变载荷作用下的多自由度运动程[13]：

(2)

[ψ]为系统的振型矩阵，作如下坐标变换：

u=[ψ]η

(3)

式中：η=[u1u2…un]T。

则系统的强迫振动方程为：

(4)

式(4)乘以特征振型[ψ]T可得：

(5)

由于[ψ]具有正交性，则主坐标下的方程为：

(6)

由主振型的正交性可得：

(7)

将式(7)代入式(6)，得到二阶微分方程：

(8)

求出每个振型坐标上的位移分量后，采用模态叠加法可以得到各个节点自由度上的位移响应：

(9)

2 风力机动态特性有限元计算方法

2.1 风力机有限元模型

选取陆地某2 MW大型风力机为研究对象，叶轮直径为40 m，轮毂直径为2.5 m，塔架高75 m。叶片材料采用玻璃纤维，弹性模量为1.76×1010Pa，泊松比为0.17；塔架为Q345E合金钢，弹性模量为2.06×1011Pa，泊松比为0.28。叶片结构采用梁单元，塔顶质量采用点质量单元划分网格。塔架底部截面采用全约束，约束该截面节点的所有自由度。风力机有限元模型如图2所示，单元节点数为93 389，单元数为462 396。

图2 风力机有限元模型Fig.2 Finite element model of the wind turbine

2.2 风力机动态响应有限元计算方法

针对上述建立的基于模态叠加法的风力机动态特性计算模型，采用有限元法对其动态特性进行计算，进行求解的具体计算流程，如图3所示。

图3 风力机整机动态特性有限元法计算流程Fig.3 Dynamic characteristics computing process forwind turbines

图3中的外部程序1是采用叶素动量定理计算叶片上的气动力[14]，主要是叶片气动性能计算和速度诱导因子迭代，Wilson等[15]详细论述了该计算过程，本文不再赘述；主程序主要是通过二次开发语言APDL将叶片迎风向载荷读入ANSYS，并将其离散成冲击载荷，在运用模态叠加法求解风力机动态响应；外部程序2主要判断风力机在求解动响应时是否达到稳定状态，周期求解动响应所允许的误差为ε，若满足式(10)，则表明动响应求解已进入稳定状态，则可停止计算；通过主程序及外部程序2求得风力机动态响应后，根据外部程序3进行扩展可进一步得到风力机的Von Mises动应力。

(10)

2.3 风力机模态分析

模态分析是动响应分析的基础，采用Block Lanczons对风力机进行模态分析并提取前10阶固有频率，如表1所示。由于篇幅有限。仅列出风力机前五阶模态振型，如图4。

表1 风力机模态计算结果

图4 风力机模态振型(1～5阶)Fig.4 Modal vibration mode (1st to 5st order) of thewind turbine

由图5可知：一阶振型主要时风力机上面两叶片沿X轴前后弯曲振动；二阶振型主要风力机整机沿X轴前后弯曲振动；三阶、四阶、五阶振型主要是叶片一阶挥舞振动。对于三阶振型，其下叶片保持不动，其余两叶片沿X轴异向振动；四阶振型是下叶片沿X轴反向振动；五阶振型为三叶片沿Y轴同向振动。

3 典型工况下风力机动态特性分析

3.1 稳态风况下风力机动态特性分析

计算可得风力机轮毂在恒定风速(即额定风速10 m/s)下的动载激励，如图5所示。根据风力机动态特性计算流程可得到风力机轮毂的位移动响应，如图6。由图5、6可知，轮毂载荷的在266.25～274.74 kN之间呈周期性变化，载荷频率由1/3次谐波、0.9 Hz主频成分及二倍频组成；轮毂位移动响应在0.88～0.80 m之间波动，动响应主频率成分为0.30 Hz，与轮毂1/3次谐波频率(叶片主频率)一致。

图7为风力机整机的Von Mises应力云图，可以看出塔架底部应力最大，平均应力达152×106Pa，且塔架底部至塔架顶部应力逐渐减小；下叶片的叶跟至叶尖的动载激励和Von Mises应力如图8所示，可以看出叶根应力最大，最大值为2.19×107Pa；从叶根至叶片11 m处时Von Mises应力呈抛物线型减小，叶片11 m处至叶尖，Von Mise应力呈减小趋势但变化不大。

图5 10 m/s恒定风速下风力机轮毂动载激励Fig.5 Dynamic loads of hub under a constant wind speedof 10 m/s

图6 10m/s恒定风速下风力机轮毂位移动响应Fig.6 Displacement dynamic response of hub under a constantwind speed of 10m/s

图7 10 m/s恒定风速下风力机应力云图Fig.7 The stress nephogram of the wind turbine under a constantwind speed of 10 m/s

图8 叶片根部至叶尖动载激励和Von Mises应力Fig.8 The loads and stress of blade from root to tip

3.2 不同风况下风力机动态特性分析

根据风力机整机动态特性计算方法计算了风力机在不同风况下的动态特性，总结不同工况下风力机位移动响应、Von Mises应力最大值及其相应位置，如表2所示。

表2 不同工况下风力机动态特性

由表2可知，风力机在不同风况下的最大位移动响应位置在都叶尖处；额定风速下风力机的位移响应最大值在所有工况中最小，为2.59 m；正常停机工况下的位移响应最大值为2.89 m，由于停机时产生了强烈的振动使得叶尖位移响应值瞬间减小到了-1.86 m；极端运行阵风作用下的位移动响应最大值在所有工况中最大，为7.19 m，较额定风速下的位移动响应最大值增加可177.22%。风力机在启动工况下的Von Mises应力最大在在所有工况中最小，为0.44×108Pa；极端运行阵风工况下的Von Mises应力最大值在所有工况中最大，为4.05×108Pa，较额定风速下的Von Mises应力最大值增加了820.45%；且极端运行工况下的Von Mises应力最大值的位置在塔架顶部，其余工况都在塔架底部。在所有工况中，极端运行阵风的位移响应最大值和Von Mises应力最大值最大，在风力机设计中应着重考虑。

3.3 极端运行阵风下风力机动态特性分析

风力机运行环境十分复杂，由于丘陵、建筑物和森林等障碍物的阻挡空气的流动会造成许多不规则的涡旋，涡旋的流动方向与空气流动方向一致时，会产生瞬时极大风速，相反，会产生瞬时极小风速。极端运行阵风是指风速突然下降，接着陡然上升，然后再突然下降，最后又上升到初始值的过程。通过风力机动态特性计算方法可以求得风力机轮毂、叶尖的位移响应(图9)，轮毂、塔架顶部的Von Mises应力(图10)。

由图9可知在稳态风(0～15 s)作用下，轮毂位移响应在0.55～0.80 m之间波动，叶尖位移响应在0.76～2.92 m之间波动。在阵风(15～40 s)作用下，风力机发生了剧烈的振动，轮毂振幅最大值达2.08 m，较稳态振幅均值增加了208.15%，最小值达在-1.18 m，较稳态振幅均值反向增加了274.18%；叶尖振动幅值最大为7.19 m，较叶尖稳态振动均值增加了290.76%，最小为-4.78 m，较稳态均值反向增加了359.78%。阵风消失后，风力机逐渐回归到稳定状态。

图9 阵风作用下轮毂、叶尖位移动响应Fig.9 Displacement dynamic responses of hub and blade tipunder transient forcing-flurry

由图10可知在稳态风(0～15 s)作用下，风力机轮毂的Von Mises应力范围为0.30×106～2.19×106Pa，塔架顶部的Von Mises应力范围为0.08×108～1.69×108Pa；在阵风(15～40 s)作用下，轮毂的Von Mises应力最大值为5.3×106Pa，相对于其稳态应力均值增加了325.70%；塔架顶部的Von Mises应力最大值为4.05×108Pa，相对于其稳态应力均值增加了357.63%。阵风结束后，轮毂及塔架顶部的Von Mises应力减小至最小值，然后逐渐增加，最终呈周期性波动。

图10 阵风作用下风力机轮毂塔架顶部的Von Mises应力Fig.10 The stress of hub and the top of tower under transientforcing-flurry

4 结 论

针对典型工况下风力机整机动态特性问题，提出了基于模态叠加法的风力机动态响应计算模型，进而采用有限元方法对其动态响应进行计算。对2 MW风力机在不同工况下进行动态特性分析，得到以下结论：

(1)稳态风作用下风力机塔架底部应力最大，平均应力达152×106Pa，且塔架底部至塔架顶部应力逐渐减小；叶片的叶跟处应力最大，最大值达2.19×107Pa；其Von Mises应力从叶根至尖处时呈抛物线型减小。

(2)计算并总结了不同工况下的风力机动态特性，发现风力机的动响应最大的位置都在叶尖处，极端运行阵风下的Von Mises应力最大位置为塔架顶部，其余工况都在塔架底部；且极端运行阵风作用下的最大位移响应值、最大Von Mises应力值在所有工况中均为最大，在风力机设计时应着重考虑。

(3)极端运行阵风条件下，叶尖的位移动响应最大值较其稳态响应均值增加了290.76%，最小值较其稳态响应均值反向增加了359.78%；塔架顶部的Von Mises应力最大值较其稳态应力均值增加了357.63%。