ABSTRACT

To analyze the vibrational frequencies and dynamic response of wind turbine blade with wind loads, the zeroth-order approximation model was adopted and the centrifugal stiffening effect was included to derive the dynamic equation of the rotating blade. The formula of stiffness matrix and mass matrix was given in the paper. Finally, the vibrational frequency and wind-induced response of a blade were calculated. The results show that the rigid-flexible coupling effect has few effects on the out-plane stiffness of blade. The in-plane stiffness of blade, however, was affected by this effect. The centrifugal stiffening effect both strengthened the in-plane and out-plane stiffness of blade. The results of the paper have some reference value to the dynamic design of wind turbine blades.

Keywords:Dynamics of Blade, Rigid-Flexible Coupling Effect, Zeroth-Order Approximation Model, Centrifugal Stiffening Effect

考虑刚柔耦合效应的风力机叶片动力学分析

高凯强，张志强

东南大学，土木工程学院，江苏 南京

收稿日期：2018年3月2日；录用日期：2018年3月16日；发布日期：2018年3月23日

摘 要

为了研究风力机叶片的固有频率以及随风动力响应，本文采用了刚–柔耦合中的零次近似模型，并考虑了离心刚化效应，推导了叶片在旋转时的动力学方程。并给出了叶片的刚度和质量矩阵的表达式。研究了某风力机叶片在不同模型下的固有频率以及随风振动响应。研究表明，刚柔耦合效应对叶片的面内刚度具有弱化作用，而对面外刚度无影响。而离心刚化效应则增强了叶片的面内外刚度。研究结果对风力叶片的动力学设计具有一定参考价值。

关键词 :叶片动力学，刚–柔耦合，零次近似模型，离心刚化

Copyright © 2018 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

叶片作为大型水平轴风力机的重要吸能部件，其动力学研究是风力机叶片设计的重要内容。叶片动力学研究主要包括叶片的旋转固有频率研究以及叶片的随风振动响应分析。

在涉及风力机叶片固有频率的研究当中，Putter等 [1] 采用了模态近似法研究了旋转梁的自振频率，Yoo等 [2] 采用瑞利–里兹法研究了预弯叶片的旋转频率，Choi等 [3] 在计算约束多体系统旋转的自振频率时，根据多体系统理论，推导了运动方程，并将之线性化。李德源 [4] 与孙宝苍 [5] 等在研究风力机叶片的模态时，采用了有限元方法并且考虑了离心刚化效应。信伟平 [6] 采用零次近似模型分析了叶片的旋转固有频率时，并考虑了刚–柔耦合效应。而在涉及风力机叶片的随风振动分析中，柯世堂等 [7] [8] 利用有限元软件建立了“叶片–塔架”耦合模型，进行了风力机的风振时域分析，并对风机关键部位的风振疲劳寿命进行了预测研究。

本文将叶片等效为变截面悬臂梁来进行研究，采用柔性体刚–柔耦合理论中的零次近似理论来对叶片进行分析，并考虑了叶片的离心刚化效应，研究了叶片的固有频率以及在风载下的动力响应。

2. 柔性梁的动力有限元方程

2.1. 刚柔耦合建模理论

首先采用有限元方法对柔性梁进行离散，将柔性梁离散成n个单元，每个单元由两个节点i和j，每个节点具有3个平动自由度和3个转动自由度，单元f节点的位移向量可表示为：

$u_f={\left\{u^i,v^i,w^i,{\theta }_1^i,{\theta }_2^i,{\theta }_3^i,u^j,v^j,w^j,{\theta }_1^j,{\theta }_2^j,{\theta }_3^j\right\}}^{\text{T}}$ (1)

在右手螺旋坐标系下，不考虑剪切变形的影响，可以利用以下形函数矩阵来表示梁中线上的任一点k的位移和转角。

$\begin{array}{l}{u}_k=N_1\left({\xi }_k\right)u_f,v_k=N_2\left({\xi }_k\right)u_f,w_k=N_3\left({\xi }_k\right)u_f,\\ {\theta }_{1k}=N_4\left({\xi }_k\right)u_f,{\theta }_{2k}=-\frac{\partial w_k}{\partial \xi }=-\frac{\partial N_3\left({\xi }_k\right)}{\partial \xi }{u}_f,{\theta }_{3k}=\frac{\partial w_k}{\partial \xi }=\frac{\partial N_2\left({\xi }_k\right)}{\partial \xi }{u}_f\end{array}$ (2)

其中， ${\xi }_k$ 为单元坐标系下，点k的纵坐标。设单元长度为l，则 ${\xi }_k=\frac{x_k-x_i}{x_j-x_i}=\frac{x_k-x_i}{l}\in \left(0,1\right)$ ，形函数矩阵为：

$\begin{array}{l}{N}_1=\left[a_1,0,0,0,0,0,a_2,0,0,0,0,0\right],\\ N_2=\left[0,b_1,0,0,0,b_2,0,b_3,0,0,0,b_4\right],\\ N_3=\left[0,0,b_5,0,b_6,0,0,0,b_7,0,b_8,0\right],\\ N_4=\left[0,0,0,a_1,0,0,0,0,a_2,0,0,0\right],\end{array}$ (3)

其中：

$\begin{array}{l}{a}_1=1-\xi ,a_2=\xi ,b_1=1-3{\xi }^2+2{\xi }^3,\\ b_2=l\left(\xi -2{\xi }^2+{\xi }^3\right),b_3=3{\xi }^2-2{\xi }^3,\\ b_4=l\left(-{\xi }^2+{\xi }^3\right),b_5=1-3{\xi }^2+2{\xi }^3,\\ b_6=l\left(-\xi +2{\xi }^2-{\xi }^3\right),b_7=3{\xi }^2-2{\xi }^3,\\ b_8=l\left(2{\xi }^2-{\xi }^3\right)\end{array}$ (4)

下面建立叶片的动力学方程，如图1所示，xyz为惯性坐标系，x’y’z’为连体坐标系，z轴和z’轴平行，连体坐标系原点固结在未变形的梁的一端上，连体坐标系绕z轴作定轴转动，角速度为ω。r₀为连体坐标系原点关于惯性坐标系的矢径，模为R。ρ₀为未变形时柔性梁非中线上任意一点k关于连体坐标系原点的矢径，p_k为连体坐标系下的变形位移矢量，则在惯性坐标系下，矢径r可表示为：

$r=r_0+A\left({\rho }_0+p_k\right)$ (5)

其中：

$\begin{array}{l}{r}_0={\left[R\cos \omega t,R\sin \omega t,0\right]}^{\text{T}},{\rho }_0={\left[x,y,z\right]}^{\text{T}}\\ A=\left[\begin{array}{ccc}\cos \omega t& -\sin \omega t& 0\\ \cos \omega t& \cos \omega t& 0\\ 0& 0& 1\end{array}\right]\\ p_k=\left[\begin{array}{c}{N}_1-y\frac{\partial N_2}{\partial \xi }-z\frac{\partial N_3}{\partial \xi }\\ \begin{array}{l}{N}_2-z{N}_4\\ N_3+y{N}_4\end{array}\end{array}\right]u_k\end{array}$ (6)

式中，A为连体坐标系关于惯性坐标系的方向余弦矩阵；u_k为p_k所在单元的节点位移，求导可得k点速

图1. 柔性梁作定轴转动

度以及加速度在惯性坐标系下的表达式为 [9] ：

$\begin{array}{l}{v}_k=\dot{r}={\dot{r}}_0+A\tilde{\alpha }\left({\rho }_0+p_k\right)+A{\dot{p}}_k\\ a_k={\ddot{r}}_0+\left(\tilde{\dot{\alpha }}A+\tilde{\alpha }\tilde{\alpha }A\right)\left({\rho }_0+p_k\right)+2\tilde{\alpha }A{\dot{p}}_k+A{\ddot{p}}_k\end{array}$ (7)

其中： $(.)$ 表示对时间求导； $\tilde{\alpha }$ 为反对称矩阵。

$\tilde{\alpha }=\left[\begin{array}{ccc}0& -\omega & 0\\ \omega & 0& 0\\ 0& 0& 0\end{array}\right]$ (8)

对于k点处的微元体dv，利用Kane方程 [10] ：

$f^{*}+f=0$ (9)

式中：f^*为微元体的广义惯性力：

$f^{*}={\left\{\frac{\partial v_k}{\partial \dot{u}}\right\}}^{\text{T}}{F}^{*},F^{*}=-\rho \cdot a_k\text{d}v$ (10)

其中，ρ为微元体的密度。

f为作用在微元体上的其他广义惯性力，应该包括该点的外力，弹性恢复力，阻尼力以及结构的内力。

最后推导得到梁的动力学方程为：

$M\ddot{u}+C\dot{u}+Ku=Q$ (11)

具体表达式见文献 [10] 。

2.2. 离心刚化效应

如图2所示，考虑叶片在xy平面内的弯曲，则微元在x轴方向上的长度改变可由下面式子给出 [11] ：

$\text{d}s-\text{d}x=\sqrt{1+{\left(\frac{\partial V}{\partial x}\right)}^2}\text{d}x-\text{d}x\approx \frac{1}{2}{\left(\frac{\partial V}{\partial x}\right)}^2\text{d}x$ (12)

则由于xy平面内弯曲而产生的离心力势能为：

$U_{exy}=\frac{1}{2}{\displaystyle \underset{0}{\overset{L}{\int }}F\left(x\right)}{\left(\frac{\partial V}{\partial x}\right)}^2\text{d}x$ (13)

其中：F(x)表示为距离梁端x处截面的离心力 [12] ：

图2. 柔性梁在xy平面内弯曲

$F\left(x\right)={\displaystyle \underset{x}{\overset{L}{\int }}\rho {\omega }^2\left(x+R\right)\text{d}V}={\displaystyle \underset{x}{\overset{L}{\int }}\rho A\left(x\right){\omega }^2\left(x+R\right)\text{d}x}$ (14)

同理由于xz平面内弯曲而产生的离心力势能为：

$U_{exz}=\frac{1}{2}{\displaystyle \underset{0}{\overset{L}{\int }}F\left(x\right)}{\left(\frac{\partial W}{\partial x}\right)}^2\text{d}x$ (15)

所以总的离心力势能为：

$U_e=\frac{1}{2}{\displaystyle \underset{0}{\overset{L}{\int }}F\left(x\right)}\left[{\left(\frac{\partial V}{\partial x}\right)}^2+{\left(\frac{\partial W}{\partial x}\right)}^2\right]\text{d}x$ (16)

已知离心惯性力势能U_e，可以采用最小势能原理来求出离心刚度矩阵K_e [13] 。

离心刚度矩阵K_e的表达式为：

$\begin{array}{c}{K}_e={\displaystyle \underset{0}{\overset{l_1}{\int }}F\left(\xi \right)\left(\frac{\partial N_2^{\text{T}}}{\partial \xi }\frac{\partial N_2}{\partial \xi }+\frac{\partial N_3^{\text{T}}}{\partial \xi }\frac{\partial N_3}{\partial \xi }\right)\text{d}\xi }\\ +{\displaystyle \sum _{i=2}^n{\displaystyle \underset{0}{\overset{l_i}{\int }}F\left(\xi +{\displaystyle \sum _{j=1}^{i-1}{l}_j}\right)\left(\frac{\partial N_2^{\text{T}}}{\partial \xi }\frac{\partial N_2}{\partial \xi }+\frac{\partial N_3^{\text{T}}}{\partial \xi }\frac{\partial N_3}{\partial \xi }\right)\text{d}\xi }}\end{array}$ (17)

其中，l_j为柔性梁第j个单元的长度。

3. 叶片的旋转固有频率分析

由式(11)可得叶片的特征方程为：

$M\ddot{u}={\beta }^{2}Ku$ (18)

将叶片等效成变截面矩形梁，截面沿长度线性变化。且截面具有两个对称轴，则可以求得：

$\begin{array}{l}M={\displaystyle \sum _{i=1}^n\rho {\displaystyle \underset{vi}{\int }{N}_1^{\text{T}}{N}_1+N_2^{\text{T}}{N}_2+N_3^{\text{T}}{N}_3+\left(y^2+z^2\right)N_4^{\text{T}}{N}_4}\text{d}V}\\ +{\displaystyle \sum _{i=1}^n\rho {\displaystyle \underset{vi}{\int }{y}^2\frac{\partial N_2^{\text{T}}}{\partial \overline{x}}\frac{\partial N_2}{\partial \overline{x}}+z^2\frac{\partial N_3^{\text{T}}}{\partial \overline{x}}\frac{\partial N_3}{\partial \overline{x}}\text{d}V}}\\ K={\displaystyle \sum _{i=1}^n{K}_{fi}}+K_d,\\ K_d=-{\displaystyle \sum _{i=1}^n{\omega }^2{\displaystyle \underset{vi}{\int }{N}_1^{\text{T}}{N}_1+N_2^{\text{T}}{N}_2+z^2{N}_4^{\text{T}}{N}_4}\text{d}V}\\ -{\displaystyle \sum _{i=1}^n{\omega }^2{\displaystyle \underset{vi}{\int }{y}^2\frac{\partial N_2^{\text{T}}}{\partial \overline{x}}\frac{\partial N_2}{\partial \overline{x}}+z^2\frac{\partial N_3^{\text{T}}}{\partial \overline{x}}\frac{\partial N_3}{\partial \overline{x}}\text{d}V}}\end{array}$ (19)

式中，ρ为叶片密度，K_fi为第i个梁单元的弹性刚度矩阵；K_d为零次近似模型中的附加刚度矩阵，它代表了叶片的刚柔耦合效应。

按照式(18)便可以计算叶片的旋转固有频率。此外，当考虑叶片的离心刚化效应时，还应将K_e计入总的刚度阵K当中。

某风力机叶片长度30 m，叶片根部宽度3 m，根部厚度0.8 m，叶尖宽度0.3 m，叶尖厚度0.08 m，叶片密度为1900 kg/m³，展向弹性模量为4.26 × 10¹⁰ N/m²，剪切弹性模量为5.5 × 10⁹ N/m²。将叶片等效为变截面矩形梁并忽略叶片的阻尼，计算叶片在不同情况下的固有频率，其结果主要如下面的表和图所示。

由表1以及图3和图4可以得出如下主要结论：

表1. 在转速为2.856 rad/s，轮毂半径为3 m时，叶片的固有频率比较(单位：Hz)

图3. 第1阶固有频率随角速度的变化

图4. 第1阶固有频率随轮毂半径的变化

1) 在标准工况下，零次近似模型的1、2、4、5阶频率与静力模型几乎完全相同。而它的3阶频率比静力模型要小，此时叶片的1、2、4、5阶频率是z向(旋转面外)弯曲，3阶频率是y向(旋转面内)弯曲。则说明K_d刚度阵对叶片旋转面内的刚度具有弱化作用，而对旋转面外的刚度影响微乎其微。

2) 在标准工况下，零次近似 + 离心刚化模型的前5阶频率都要比静力和零次近似模型大，说明K_e刚度阵对叶片旋转面内和面外都产生刚化作用，而且在面内盖过了K_d矩阵弱化的影响。

3) 观察叶片频率随轮毂半径的变化可知，零次近似 + 离心刚化模型随着轮毂半径的增大而增大，这是因为K_e矩阵的表达式与轮毂半径R有关，且轮毂半径越大，它们对叶片所起的刚化作用越强。而K_d矩阵与轮毂半径无关，所以零次近似模型和静力模型相同。

4) 观察叶片频率随旋转角速度的变化可知，当角速度的增大至30 rad/s左右时，零次近似模型中1阶频率开始减小，这是因为叶片的1阶振型已由原来的z向弯曲变成y向弯曲，即由面外弯曲变为面内弯曲，其振动形态发生了改变。即随着角速度的增大，受到K_d的影响，叶片的面内弯曲频率开始减小，当其值小于面外弯曲时，就变成了1阶频率。而且随着角速度的进一步增大，最后叶片的旋转固有频率减小为零，此时叶片在面内的弯曲变形已经发散。

4. 叶片随风动力响应分析

4.1. 风荷载计算

作用于叶片的风荷载可分为顺风向荷载和垂直于风向上的荷载，本文仅考虑叶片顺风向的荷载。风荷载由平均风速和脉动风速所产生，根据风压理论，风荷载可以表达为 [14] ：

$f_w=\frac{1}{2}{C}_d{\rho }_{\text{air}}B{\left(u_m+u_f\right)}^2$ (20)

其中：C_d为拖曳力系数，ρ_air为空气密度，u_m和u_f中分别为叶片上的平均风速和脉动风速。B为风荷载作用面积。

$u_m=V_h{\left(\frac{r\cdot \sin \left({\omega }_{0}t\right)+h}{h}\right)}^{\alpha }$ (21)

其中：V_h为轮毂中心处的风速；ω₀为叶片角速度；r为叶片上一点到轮毂中心的距离；h为轮毂中心的高度。

叶片的脉动风速采用谐波叠加法进行计算，请参考文献 [15] ，此处不再赘述。

4.2. 叶片随风振动分析

计算得到了叶片的风荷载之后，便可以进行叶片的风振响应研究。叶片的基本参数同上一节，其中，叶片的拖曳力系数C_d [14] 为2，空气密度1.25 kg/m³，叶片匀速转动的角速度为ω₀ = 1.988 rad/s，阻尼比取0.008。与风有关的参数为：轮毂中心高65 m，该中心处的风速为10.76 m/s，风廓线指数α为0.2，风谱为Von Karman谱。

图5是叶片叶尖在不同模型下计算得到的随风响应时程结果。

由图5可以看出：

1) 零次近似模型和静力模型计算的结果基本一致，这是因为K_d刚度阵只对旋转面内的刚度有影响。

2) 零次近似 + 离心刚化模型的z向位移相比静力模型要小，这说明K_e矩阵对叶片旋转面外的刚度具有增强效果，能够降低叶片在随机风载作用下的位移。

图5. 不同模型下，叶尖的随风响应时程曲线

5. 结论

1) 叶片频率的研究表明，K_d对叶片的旋转面内刚度产生软化作用，降低了叶片的振动频率，且易使得叶片在旋转面内的弯曲变形增大。当转速达到一定程度时，结构便遭到了破坏。而K_e对叶片旋转面内和面外的刚度均具有增强效果，提高了结构的振动频率。

2) 叶片随风振动的分析表明，零次近似模型的计算结果和静力模型基本相同，而零次近似 + 离心刚化模型的风振响应较前两者要小，说明了K_d刚度矩阵对叶片的面外刚度不产生影响作用，而离心刚化效应使叶片的面内外刚度都得到了增强，所以能够降低叶片在随机风载作用下的位移。因此，今后的叶片风振响应分析当中(即垂直于平面的运动)，可以不考虑K_d的影响，而叶片的离心刚化效应则需要加以考虑。

基金项目

国家自然科学基金(51278106, 51438002)，江苏省产学研前瞻性项目(BY2016076-11)，论文编号是2750499。

文章引用

高凯强,张志强. 考虑刚柔耦合效应的风力机叶片动力学分析
Dynamic Analysis of Wind Turbine Blade with the Rigid-Flexible Coupling Effect[J]. 土木工程, 2018, 07(02): 205-213. https://doi.org/10.12677/HJCE.2018.72025

参考文献