1. 引言
叶片作为大型水平轴风力机的重要吸能部件,其动力学研究是风力机叶片设计的重要内容。叶片动力学研究主要包括叶片的旋转固有频率研究以及叶片的随风振动响应分析。
在涉及风力机叶片固有频率的研究当中,Putter等 [1] 采用了模态近似法研究了旋转梁的自振频率,Yoo等 [2] 采用瑞利–里兹法研究了预弯叶片的旋转频率,Choi等 [3] 在计算约束多体系统旋转的自振频率时,根据多体系统理论,推导了运动方程,并将之线性化。李德源 [4] 与孙宝苍 [5] 等在研究风力机叶片的模态时,采用了有限元方法并且考虑了离心刚化效应。信伟平 [6] 采用零次近似模型分析了叶片的旋转固有频率时,并考虑了刚–柔耦合效应。而在涉及风力机叶片的随风振动分析中,柯世堂等 [7] [8] 利用有限元软件建立了“叶片–塔架”耦合模型,进行了风力机的风振时域分析,并对风机关键部位的风振疲劳寿命进行了预测研究。
本文将叶片等效为变截面悬臂梁来进行研究,采用柔性体刚–柔耦合理论中的零次近似理论来对叶片进行分析,并考虑了叶片的离心刚化效应,研究了叶片的固有频率以及在风载下的动力响应。
2. 柔性梁的动力有限元方程
2.1. 刚柔耦合建模理论
首先采用有限元方法对柔性梁进行离散,将柔性梁离散成n个单元,每个单元由两个节点i和j,每个节点具有3个平动自由度和3个转动自由度,单元f节点的位移向量可表示为:
(1)
在右手螺旋坐标系下,不考虑剪切变形的影响,可以利用以下形函数矩阵来表示梁中线上的任一点k的位移和转角。
(2)
其中,
为单元坐标系下,点k的纵坐标。设单元长度为l,则
,形函数矩阵为:
(3)
其中:
(4)
下面建立叶片的动力学方程,如图1所示,xyz为惯性坐标系,x’y’z’为连体坐标系,z轴和z’轴平行,连体坐标系原点固结在未变形的梁的一端上,连体坐标系绕z轴作定轴转动,角速度为ω。r0为连体坐标系原点关于惯性坐标系的矢径,模为R。ρ0为未变形时柔性梁非中线上任意一点k关于连体坐标系原点的矢径,pk为连体坐标系下的变形位移矢量,则在惯性坐标系下,矢径r可表示为:
(5)
其中:
(6)
式中,A为连体坐标系关于惯性坐标系的方向余弦矩阵;uk为pk所在单元的节点位移,求导可得k点速
![](//html.hanspub.org/file/11-2750499x17_hanspub.png)
Figure 1. The beam rotates around a fixed axis
图1. 柔性梁作定轴转动
度以及加速度在惯性坐标系下的表达式为 [9] :
(7)
其中:
表示对时间求导;
为反对称矩阵。
(8)
对于k点处的微元体dv,利用Kane方程 [10] :
(9)
式中:f*为微元体的广义惯性力:
(10)
其中,ρ为微元体的密度。
f为作用在微元体上的其他广义惯性力,应该包括该点的外力,弹性恢复力,阻尼力以及结构的内力。
最后推导得到梁的动力学方程为:
(11)
具体表达式见文献 [10] 。
2.2. 离心刚化效应
如图2所示,考虑叶片在xy平面内的弯曲,则微元在x轴方向上的长度改变可由下面式子给出 [11] :
(12)
则由于xy平面内弯曲而产生的离心力势能为:
(13)
其中:F(x)表示为距离梁端x处截面的离心力 [12] :
![](//html.hanspub.org/file/11-2750499x18_hanspub.png)
Figure 2. The bending of flexible beam in the xy plane
图2. 柔性梁在xy平面内弯曲
(14)
同理由于xz平面内弯曲而产生的离心力势能为:
(15)
所以总的离心力势能为:
(16)
已知离心惯性力势能Ue,可以采用最小势能原理来求出离心刚度矩阵Ke [13] 。
离心刚度矩阵Ke的表达式为:
(17)
其中,lj为柔性梁第j个单元的长度。
3. 叶片的旋转固有频率分析
由式(11)可得叶片的特征方程为:
(18)
将叶片等效成变截面矩形梁,截面沿长度线性变化。且截面具有两个对称轴,则可以求得:
(19)
式中,ρ为叶片密度,Kfi为第i个梁单元的弹性刚度矩阵;Kd为零次近似模型中的附加刚度矩阵,它代表了叶片的刚柔耦合效应。
按照式(18)便可以计算叶片的旋转固有频率。此外,当考虑叶片的离心刚化效应时,还应将Ke计入总的刚度阵K当中。
某风力机叶片长度30 m,叶片根部宽度3 m,根部厚度0.8 m,叶尖宽度0.3 m,叶尖厚度0.08 m,叶片密度为1900 kg/m3,展向弹性模量为4.26 × 1010 N/m2,剪切弹性模量为5.5 × 109 N/m2。将叶片等效为变截面矩形梁并忽略叶片的阻尼,计算叶片在不同情况下的固有频率,其结果主要如下面的表和图所示。
由表1以及图3和图4可以得出如下主要结论:
![](//html.hanspub.org/file/11-2750499x36_hanspub.png)
Figure 3. The first natural frequency varied with the change of rotation speed
图3. 第1阶固有频率随角速度的变化
![](//html.hanspub.org/file/11-2750499x35_hanspub.png)
Figure 4. The first natural frequency varied with the change of radius of hub
图4. 第1阶固有频率随轮毂半径的变化
1) 在标准工况下,零次近似模型的1、2、4、5阶频率与静力模型几乎完全相同。而它的3阶频率比静力模型要小,此时叶片的1、2、4、5阶频率是z向(旋转面外)弯曲,3阶频率是y向(旋转面内)弯曲。则说明Kd刚度阵对叶片旋转面内的刚度具有弱化作用,而对旋转面外的刚度影响微乎其微。
2) 在标准工况下,零次近似 + 离心刚化模型的前5阶频率都要比静力和零次近似模型大,说明Ke刚度阵对叶片旋转面内和面外都产生刚化作用,而且在面内盖过了Kd矩阵弱化的影响。
3) 观察叶片频率随轮毂半径的变化可知,零次近似 + 离心刚化模型随着轮毂半径的增大而增大,这是因为Ke矩阵的表达式与轮毂半径R有关,且轮毂半径越大,它们对叶片所起的刚化作用越强。而Kd矩阵与轮毂半径无关,所以零次近似模型和静力模型相同。
4) 观察叶片频率随旋转角速度的变化可知,当角速度的增大至30 rad/s左右时,零次近似模型中1阶频率开始减小,这是因为叶片的1阶振型已由原来的z向弯曲变成y向弯曲,即由面外弯曲变为面内弯曲,其振动形态发生了改变。即随着角速度的增大,受到Kd的影响,叶片的面内弯曲频率开始减小,当其值小于面外弯曲时,就变成了1阶频率。而且随着角速度的进一步增大,最后叶片的旋转固有频率减小为零,此时叶片在面内的弯曲变形已经发散。
4. 叶片随风动力响应分析
4.1. 风荷载计算
作用于叶片的风荷载可分为顺风向荷载和垂直于风向上的荷载,本文仅考虑叶片顺风向的荷载。风荷载由平均风速和脉动风速所产生,根据风压理论,风荷载可以表达为 [14] :
(20)
其中:Cd为拖曳力系数,ρair为空气密度,um和uf中分别为叶片上的平均风速和脉动风速。B为风荷载作用面积。
(21)
其中:Vh为轮毂中心处的风速;ω0为叶片角速度;r为叶片上一点到轮毂中心的距离;h为轮毂中心的高度。
叶片的脉动风速采用谐波叠加法进行计算,请参考文献 [15] ,此处不再赘述。
4.2. 叶片随风振动分析
计算得到了叶片的风荷载之后,便可以进行叶片的风振响应研究。叶片的基本参数同上一节,其中,叶片的拖曳力系数Cd [14] 为2,空气密度1.25 kg/m3,叶片匀速转动的角速度为ω0 = 1.988 rad/s,阻尼比取0.008。与风有关的参数为:轮毂中心高65 m,该中心处的风速为10.76 m/s,风廓线指数α为0.2,风谱为Von Karman谱。
图5是叶片叶尖在不同模型下计算得到的随风响应时程结果。
由图5可以看出:
1) 零次近似模型和静力模型计算的结果基本一致,这是因为Kd刚度阵只对旋转面内的刚度有影响。
2) 零次近似 + 离心刚化模型的z向位移相比静力模型要小,这说明Ke矩阵对叶片旋转面外的刚度具有增强效果,能够降低叶片在随机风载作用下的位移。
![](//html.hanspub.org/file/11-2750499x38_hanspub.png)
Figure 5. The time history curves of wind-induced response of blade tip under different models
图5. 不同模型下,叶尖的随风响应时程曲线
5. 结论
1) 叶片频率的研究表明,Kd对叶片的旋转面内刚度产生软化作用,降低了叶片的振动频率,且易使得叶片在旋转面内的弯曲变形增大。当转速达到一定程度时,结构便遭到了破坏。而Ke对叶片旋转面内和面外的刚度均具有增强效果,提高了结构的振动频率。
2) 叶片随风振动的分析表明,零次近似模型的计算结果和静力模型基本相同,而零次近似 + 离心刚化模型的风振响应较前两者要小,说明了Kd刚度矩阵对叶片的面外刚度不产生影响作用,而离心刚化效应使叶片的面内外刚度都得到了增强,所以能够降低叶片在随机风载作用下的位移。因此,今后的叶片风振响应分析当中(即垂直于平面的运动),可以不考虑Kd的影响,而叶片的离心刚化效应则需要加以考虑。
基金项目
国家自然科学基金(51278106, 51438002),江苏省产学研前瞻性项目(BY2016076-11),论文编号是2750499。