风力机结构动力学新模型与动态响应求解计算

宁 伟，孙文磊

（新疆大学机械工程学院，新疆 乌鲁木齐 830049）

1 引言

为了在大规模利用风能的同时降低生产成本，大型化是风力发电机组目前的发展趋势。于此同时，风力机高度以及体积的持续增大，使其制造难度与成本不断升高，其承受的载荷也更加难以控制并对风力发电机组的安全运行造成影响，甚至产生破坏。因此，风力发电机固有频率、动态响应等研究对于风力发电机组的设计与优化具有重要意义。

目前，在风力机结构力学研究方面，文献[1-2]将叶片与塔架看成柔性悬臂梁，使用二节点梁单元进行离散，分别考虑了气动阻尼、离心刚化的影响。在此基础上，文献[3]采用kane方法建立风力机动力学模型，通过假设模态法离散进行柔性化。文献[4]考虑材料非线性以及截面偏心的影响，建立一种叶片非线性单元，并应用于整机分析。针对上述方法建模复杂、重复性操作过多的问题，文献[5-6]采用刚体有限元法研究结构阻尼对振动特性的影响。文献[7]采用超级单元方法将柔性叶片通过力元弹簧及阻尼器进行离散推导出叶片多体动力学方程组。

根据以上调研，在柔性多体动力学理论基础上，选取新颖、简单、计算速度快的刚体有限元方法（Rigid Finite Element Method，RFEM），建立了系统动力学模型，编写了相关程序。计算结果与GH-Bladed数据结果对比，验证了刚体有限元方法和模型建立的准确性，适用于大型风力发电机的动力学分析计算。

2 数学模型

2.1 刚体有限元方法

将塔架和叶片简化为质量与刚度连续分布的柔性悬臂梁，采用刚体有限元二次划分方法，通过刚体单元与弹性阻尼单元分别对塔架和叶片进行离散。刚体有限元法，如图1所示。

图1 刚体有限元法Fig.1 Rigid Finite Element Method

刚体有限元法中，刚体单元存储位移与速度，弹簧阻尼节点则存储变形能。考虑每个刚体单元有6个自由度，则任意刚体单元的自由度表示为：

在局部坐标系｛｝E中任意刚体单元i上某点A的位置矢量r′与该点在全局坐标系｛｝0下的位置矢量r有以下关系：

由此推导出刚体单元的动能Ei：

将系统所有刚体单元与弹性阻尼节点相关量整合，得出Lagrange方程建立的动力学模型[8]：

式中：Q—系统所受的广义力；假设刚体单元的3个旋转自由度方向的转角较小，则各阻尼节点的势能与阻尼耗散能可推导为：

式中：ki=diag（），ci=diag（）；弹簧刚度系数及阻尼系数基于Kelvin-Voigt模型推导。

各刚体单元动能的Lagrange方程变换形式为：

式中：（lj）=1…6—对该刚体段变换矩阵偏导的各自由度方向，Mi=diag（mi，mi，mi）。

将式（10）～式（11）写为矩阵的形式：

同理，由式（8）可得系统中各弹性阻尼节点的阻尼矩阵为：

最终，推导出式（6）线性化后的统一形式的结构动力学方程组[9]：

2.2 固有动力特性

结构作自由振动时，系统的激励为0，式（14）可写为：

令节点位移 q（t）=φsin（ωt+θ），得特征矩阵方程：

对系统固有频率及振型的求解就归结到式（16）的特征值与特征向量的求解，上式有非零解的条件为系数行列式为0，即：

得出其特征值与特征向量：

3 动态响应计算

由于模态叠加法相比于其它动力学分析方法具有高效、计算速度快的优点，动态响应计算采用了模态叠加法，通过在模态计算中所获得的固有频率与振型来描述结构的动态响应。

3.1 模态叠加原理

将式（20）代入式（14），则系统的强迫振动方程：

由于振型具有正交性，假设阻尼矩阵也符合正交性。则可将式（22）离散化为非耦合的n个单自由度方程：

式中：Mi—第 i阶模态质量；Ci—第 i阶模态阻尼；Ki—第 i阶模态刚度，在式（23）两边同时除以模态质量Mi：

经坐标变换，物理坐标q转换为模态坐标x：

式中：ξi—第 i阶阻尼比；ωi=—第i阶振动频率。

在求解每个振型下的位移分量后带入式（20）最终通过模态叠加法得到总位移：

3.2运动方程求解

由Newmark法可知各刚体段在每个时间步长下的位移、速度表达式：

式中：h—步长，α=0.25，β=0.5，将上式带入式（23），得到加速度表达式：

求解式（26）～式（28），取时间步长 0.01s，得出位移、速度、加速度的动态响应。

总结动态响应的求解过程：

（1）计算固有频率及振型。

（2）由振型的正交特性，分别计算各阶模态质量、模态阻尼与模态刚度。

（3）由每个时间步长下所受的载荷，计算得出各阶模态载荷。

（4）离散方程为非耦合的n个单自由度方程。

（5）分别给定位移、速度和加速度的初始值，通过Newmark积分法求解各个单自由度方程。

（6）由式（25）得到动态响应。

4 风轮与塔架耦合

由于受风轮载荷的影响，塔架将产生前后、左右方向的摆动，同时系统集合形状也在不断变化，若用标准的模态分析来分析风力机的整体动态特性会变得越来越复杂。因此，需要将模态振型与各阶频率分开考虑，塔架模态中将风轮视为刚性的，而风轮模态中风轮类似于悬挂在刚性主轴上，如图2所示。

图2 风轮和塔架的基本形状Fig.2 The Basic Shape of Wind Wheel and Tower

由于风轮与塔架模态不成直角，因此叶片与塔架的轨迹方程彼此关联，则由图2可知在某时刻叶片J的偏转可写为[10]：

式中：φTJ—相对于刚性塔架的叶片模态振型；

φTJ—叶片J的位置偏差；

φTJ=1cosψ，ψ—叶片 J的方位角；

φT—塔架模态振型。

将式（29）带入式（14），则叶片的运动方程为：

其中，塔架模态质量与塔架、机舱、风轮的质量及风轮的转动惯量相关，Mi=MTφT+MN+MR+IR/L2。

式（30）和式（31）的动态分析过程为：

（1）带入初始时刻的位移、速度、加速度和空气动力载荷。

（2）假设方程右边耦合项为不随时间改变的常数，以便在解耦后的方程中消去同类项。

（3）解方程算出位移与速度的增量，从而得出叶片J的叶尖位移与速度增量。

（4）解运式（30）、式（31），得出最后一个时间步长下的加速度。

（5）改变方程右边的耦合项，重新求解运动增量方程，得到位移和速度增量的修正值。

5 算例分析

以美国国家新能源实验室的某NREL1.5WM风力发电机[11]为研究对象，对其进行模态与动力学分析，主要参数，如表1所示。

表1 风力机结构参数Tab.1 Structural Parameters of Wind Turbine

5.1 模态分析

由于在风力发电机模态分析中，低频振动要比高频振动危险，高阶模态对响应的影响较小，同时结构阻尼的作用使响应中高阶部分衰减很快，因此忽略了高阶模态的影响[12]。首先计算叶片平面外与平面内前三阶固有频率，叶片固有频率和振型的计算结果同GH-Bladed仿真结果对比，如表2、图3所示。

表2 叶片模态Tab.2 Natural Frequencies of Blade

图3 叶片固有振型Fig.3 Mode Shapesof Blade

计算塔架前后与左右方向的前两阶弯曲固有频率，其弯曲固有频率和振型结果的对比，如表3、图4所示。

表3 塔架模态分析Tab.3 Natural Frequencies of Tower

图4 塔架固有振型Fig.4 Mode Shapesof Tower

由于塔架的对称性，一阶前后向与左右向、二阶前后向与左右向的振型相同。从塔架振型图可以看出：一阶振型塔架顶端的振幅最大，二阶振型塔架顶端的振幅较小，中间的振幅较大。

通过模态分析对比，初步验证了模型建立的正确性，同时为动态响应的计算提供了计算参数。

5.2 动态响应分析

风力发电机组运行时，对于风轮所产生的时变动态载荷，通过GH-Bladed来获取，轮毂高度处的纵向风速，如图5所示。湍流风采用von-karman模型，平均风速为11.8m/s，仿真时间70s，步长0.05s。

图5 纵向风速Fig.5 Longitudinal Wind Speed

由于风力机在来流方向上的变形远远大于其左右方向变形，因此先考虑了叶尖平面外和塔架前后向的变形情况，计算对叶片的前三阶模态与塔架的前两阶模态进行了叠加。

叶片与塔架的动态响应计算结果分别，如图6、图7所示。从分析结果可以看出实际工作状况下，叶尖在风轮平面外的偏移量较大，最大值为2.59 m，由于数值计算时初始位移、速度、加速度的值均从零开始，因此在起始时刻叶片与塔架的位移与速度均出现剧烈变化7s后趋于稳定，计算同GH-Bladed仿真结果保持一致。

图6 叶片尖端动态响应Fig.6 The Dynamic Analysis of Bladed Tip

图7 塔架顶部动态响应Fig.7 The Dynamic Analysis of Tower Top

叶片与塔架在来流方向较大的变形量不仅影响气动载荷的变化，更对风力发电机组的运行安全造成较大影响，甚至可能会产生破坏。因此，设计过程中对其来流方向的动态响应需要充分考虑。

仿真结果表明：刚体有限元法在兼顾较少的建模计算时间的情况下，能够较为准确地进行风力发电机动力学分析。

6 结论

（1）针对风力发电机组动力学问题，利用刚体有限元法进行了离散化建模，建立了叶片与塔架耦合的结构动力学模型，给出了刚体有限元法的模态分析与动态响应求解方法与步骤。

（2）利用模型对1.5MW风力机的模态与动态响应进行计算。计算结果同GH-Bladed进行了对比验证，保持了良好的一致性，证明建立的模型简化合理、计算准确，能较好预测风力机组对于时变载荷的动态响应，具有良好的实用价值，适用于大型风力机的动力学分析计算。