风力机旋转叶片的刚柔耦合动力学响应特性分析

赵荣珍;芦颉;苏利营

【摘 要】为获得风力机风轮在时变载荷作用下的动力学响应规律,在将轮毂假设为刚性圆盘和叶片假设为柔性悬臂梁的基础上,考虑剪切应变对叶片引起的附加位移、叶片的旋转运动与弹性变形间的耦合及离心力的作用,运用Hamilton原理建立了旋转叶片子系统的非线性动力学模型.以河西地区某风场1.2MW风电机组为例,采用有限元法和Newmark数值积分法对其叶片在时变载荷作用下的动态响应进行了计算.结果表明:考虑剪切变形影响时得到的叶片振动幅值比不考虑剪切变形影响时会平均增大约7.5％;在考虑叶片的旋转运动与弹性变形间的耦合作用时,离心力对叶片振动动态特性的影响将会被弱化. 【期刊名称】《兰州理工大学学报》 【年(卷),期】2016(042)006 【总页数】7页(P36-42)

【关键词】风力机;叶片;动态响应;刚柔耦合;剪切变形 【作 者】赵荣珍;芦颉;苏利营

【作者单位】兰州理工大学机电工程学院,甘肃兰州 730050;兰州理工大学机电工程学院,甘肃兰州 730050;兰州理工大学机电工程学院,甘肃兰州 730050 【正文语种】中 文 【中图分类】TH113

叶片是风电机组的最关键部件,它的结构动态性能优劣将直接影响到整机性能.在叶片绕风轮轴作大范围空间运动的过程中,叶片的运动与其弹性变形间[1]、弹性变形与气动力间的耦合,将使叶片的非线性振动愈加严重[2]，并且这种非线性振动会随着风电机组向大型化和柔性化发展变得更加严重和复杂.因此,研究风电机组叶片在时变载荷作用下的动态响应，具有重要的理论意义及工程应用价值[3-4]. 风力机叶片展向长、弦向短,因此目前多数风力机叶片的结构动力学分析及优化设计研究都采用欧拉梁模型[5]进行,计算方法多采用牛顿法从非线性应变位移关系式导出叶片的挥舞、摆振的全耦合非线性偏微分方程,然后进行计算求解[6-7].这种方法计算量大,而且存在着对高阶情况很难求解的缺陷.而欧拉梁模型假设变形前后垂直于中面的截面仍然保持垂直关系,忽略了横向剪切变形的影响.欧拉梁模型虽然能够在一定程度上反映叶片的整体变形性能,但是为了获得更加精确的变形,应该考虑横向剪切变形的影响.

基于此,本文在考虑剪切变形引起的附加位移、叶片的旋转运动与弹性变形间的耦合以及离心力影响作用的基础上,拟采用刚柔耦合动力学建模方法,运用Hamilton原理建立旋转叶片子系统的动力学模型.应用所建立的模型,对河西地区某风场1.2 MW风电机组工作时叶片的动态响应进行了计算,探讨剪切变形引起的附加位移、叶片的旋转运动与弹性变形的耦合及离心力对叶片动态特性的影响.

考虑刚柔耦合效应的柔性多体系统动力学称之为刚柔耦合动力学.该研究方向主要探讨柔性体的空间旋转运动与其弹性变形之间的耦合关系,以及这种耦合所导致的动力学效应[8].在刚柔耦合系统中对运动的描述多采用相对描述,即分别建立各物体的参考系,将运动分解为整体运动和变形运动两部分,然后用相对运动变量建立系统的动力学方程.方程的建立主要有两种方法:牛顿-欧拉法分别采用动量定理和动量矩定理描述柔性体的平动和转动,且都包含了柔性体的变形运动；而Hamilton原理则是用动能、势能的变分项代替弹性力和惯性力,其优点是计算过程只与纯粹的标

量有关,因此具有计算简便的优点[9].

Hamilton原理是一个积分变分原理,其数学表达式为

式中：T为体系的总动能；U为体系的位能；W为作用于体系上的非保守力所做的功；δ为指定时间区间内所取的变分.

针对叶片展弦比大的特点,研究中采用将叶片简化为悬臂梁,轮毂等效为刚性圆盘的简化方法,建立了如图1所示的运动模型.图2为叶片截面翼型图.图1中o1x1y1z1是以风轮为中心的风轮轴坐标系,假设它固定不动,o1z1为风轮转轴,oxyz为固结在叶根处绕风轮转轴转动的叶片相对坐标系.P为未变形时叶片上任意一点,柔性叶片在旋转过程中发生了变形,变形后P点到达P′点,其中u为P点在旋转平面内的横向位移,v为P点在垂直于旋转平面的横向位移.图2中o为叶片截面形心,ζ和η为形心主惯性轴,θb为截面扭转角,(°).旋转叶片子系统的参数表示如下：e为轮毂半径,m；J为轮毂的转动惯量,kg·m2；L为叶片长度,m；ρ为密度,kg/m3；E为弹性模量；G为剪切模量,Pa；I为截面惯性矩；θ为叶片旋转运动的角位移,(°). P′点关于坐标系o1x1y1z1的坐标阵可表达为

考虑剪切变形的影响,则P点在摆振和挥舞两方向的横向位移可表示为

式中:ub和vb为弯曲变形引起的横向位移；us和vs为剪切变形引起的附加横向位移.

截面的抗弯刚度按下式计算[10]：

式中：IF为截面对ζ轴的惯性矩；IFc为截面对η轴的惯性矩. 3.1 系统动能

风力机叶片固定在刚性轮毂上,并且随轮毂以一定角速度旋转,旋转叶片子系统的动能主要包括轮毂的动能和叶片的动能两部分： 式中：m为单位长度叶片的质量,kg/m. 3.2 系统势能

叶片旋转过程中系统的势能主要有弯曲应变能U1、剪切应变能U2、重力势能U3以及旋转过程中由离心力产生的离心力势能U4四部分组成: 弯曲应变能和剪切应变能的计算公式如下： 式中：A为截面积,m2；k为校正因子.

叶片所受的重力势能随着转动角θ的不同而不同.因此在叶片坐标系中,如不考虑主轴倾角和叶片锥角的影响,则系统的重力势能为 式中：g为重力加速度,m/s2.

叶片在旋转过程中产生的离心力不仅在轴向上有能量变化,而且在旋转平面内也有能量变化,其方向总是沿着叶片向外,离心力势能按下式计算： 式中：Px为叶片任意截面的轴向离心力ξ. 3.3 外力所做虚功

机组运行过程中,作用在叶片上的空气动力是整个系统的动力源.考虑失速和动态入流的影响,采用修正的叶素动量理论求解气动载荷[11].稳态工况下,单位长度叶片在弦线方向和其垂直方向上的气动力可分别表示为如下形式：

式中：CL为截面升力系数；CD为截面阻力系数；ρa为空气密度,kg/m3；va为来流速度,m/s；C为截面弦长,m；φ为来流角,(°). 叶片在摆振和挥舞两个方向的气动力如下：

以风轮旋转平面作为参照,叶片的摆振力矩和挥舞力矩分别为 考虑阻尼力的影响,非保守力所做的虚功为 式中为阻尼系数. 3.4 动力学方程的建立

分别计算动能和势能关于ub、us、vb、vs和θ的变分并代入式(1),则建立的旋转叶片系统的动力学方程组为

4.1 动力学方程的离散化处理方法

由于式(19～23)动力学方程是一组非线性、时变和强耦合的偏微分-积分方程组,通常不易得到精确的解析解.因此,采用有限元方法将叶片离散成有限自由度作为近似求解分析模型.在建模过程中,由于考虑了剪切应变的影响,因此采用修正的欧拉梁单元进行离散.设叶片某一单元长度为l,单位长度的质量为m,单元节点位移参数分别为ueb、ues、veb和ves,单元形函数矩阵为Nb和Ns[12].单元内任一点在摆振方向和挥舞方向的位移可以表示为

将式(24,25)代入式(19～23)的动力学方程组中,可得旋转叶片的离散动力学方程: 上式中的具体元素均由式(19～25)推导得出,在后续的求解过程中将会说明,此处限于篇幅不再赘述.

4.2 动力学方程的求解过程

风力机叶片在旋转过程中绕风轮轴作大范围旋转运动,其转动角速度规律是已知的不用求解,具体形式将在算例分析部分详细给出.故忽略式(26)中的第一行,则叶片在非惯性系下的动力学方程为

叶片静止时,系统的质量矩阵和刚度矩阵分别为叶片的质量矩阵和刚度矩阵,方程(27)中的各元素为

旋转过程中考虑叶片的旋转运动与弹性变形耦合引起的动力刚度矩阵和离心力引起的几何刚度矩阵对系统刚度矩阵的影响,系统的总刚度矩阵为

以上方程中其余项与叶片静止时相同.计入横向剪切变形时,考虑截面转动效应对质量矩阵的影响,以及剪切应变对刚度矩阵的影响,系统的质量矩阵和刚度矩阵分别描述为(阻尼矩阵与静止时相同) 式中

采用Newmark逐步积分法对式(27)进行数值求解.首先给定初始时刻叶片的位移、速度和加速度,即和求出在第一个时间步长结束时叶片的位移、速度和加速度和

叶片在时变载荷作用下的动态响应计算步骤设置如下：

1) 给定初始值.本研究中设叶片是从静止状态下旋转的,因此取和均为0. 2) 计算叶片的质量矩阵M,阻尼矩阵C,刚度矩阵K.

3) 选取时间步长t和参数α、β.本研究选取时间步长为0.05 s,参数α=0.25,β=0.5. 4) 对于每一个时间步长计算作用在叶片上的等效载荷Qi,其中i表示第i个时间步长.

5) 由式(39～41)计算第一个时间步长结束时叶片的位移、速度和加速度. 6) 将上述步骤循环,得到所有时刻叶片的位移.

以河西地区某风场1.2 MW水平轴风电机组叶轮子系统为例,应用所建立的动力学模型对其固有频率和动态响应进行了求解分析.该机组叶轮参数如下：轮毂半径为2 m,叶片长为29 m;其材料为玻璃钢复合材料，剪切模量G为5.5 GPa,其他几何数据及质量和刚度分布见表1.其中,挥舞刚度与摆振刚度是叶片的固有属性,其计算公式对应于式(33)中的Ky和Kz. 5.1 叶片的固有频率特性

考虑剪切应变、叶片大范围旋转运动与自身弹性变形间的耦合以及离心力对叶片固有频率的影响,在式(27)中忽略结构阻尼和外力的影响作用,计算风轮转速为12 r/min时,叶片静止、旋转及考虑剪切应变时的弯曲固有频率,计算结果见表2. 从表2可见,对于叶片的6种不同振动状态,无论摆振还是挥舞模态,叶片旋转时的固有频率都大于静止时的固有频率,说明离心力作用使叶片固有频率增大；而横向剪切应变下的叶片固有频率却比静止时的要小,由此表明,受截面转动效应的影响,叶片的质量矩阵增大将导致其固有频率减小.进一步对比可看出,挥舞模态的影响较摆振模态显著,一阶挥舞模态的频率增量为4.8%,摆振模态频率增量为1.17%,这是由于叶片大范围旋转运动与弹性变形之间的耦合削弱了离心力作用对固有频率的影响,导致摆振刚度弱化.

5.2 叶片动态响应的计算结果

假定叶片由静止状态开始启动,t1时刻开始以恒定角速度ω转动,t2时刻开始制动,制动时间为t3,叶片转动角速度规律取与文献[14]相同,即

本研究取t1=30 s,t2=90 s,t3=(t2+30)s.研究风电机组在湍流强度为16.1%、平均风速为12 m/s的湍流风场下启动、正常运转和停车等工况下叶片的动态响应.求解式(27)的动力学方程,获得叶尖挥舞方向和摆振方向的位移,结果如图3和图4所示.由图3和图4对比可知,启动过程中,无论是挥舞还是摆振叶尖振动幅值都较大,而正常运转时挥舞方向的振动明显大于摆振方向.

若不考虑剪切变形引起的附加位移,则计算结果与GH bladed软件模拟结果基本一致,说明本研究所建立的模型是准确可靠的.图5和图6显示了叶片旋转运动与弹性变形的耦合和离心力、剪切变形引起的附加位移对叶尖位移的影响,由于叶片在挥舞方向的变形比摆振方向大很多,因此仅以挥舞方向振动为例进行说明.由图3和图5对比可知,虽然大范围运动与弹性变形的耦合效应使得叶片的变形增大,但是由于离心力的作用使得叶片的变形减小,二者综合导致叶片的整体变形减小.对比图3和图6可以看出,剪切变形引起的叶片附加位移对叶尖变形影响较大,在整个仿真过程中,叶尖振动幅值平均增大了7.5%.由此可见,剪切变形引起的附加位移对旋转叶片子系统的动态响应有较大的影响,在计算过程中应予以充分考虑.

本文计入剪切变形引起的附加位移的影响,同时考虑叶片大范围旋转运动与自身弹性变形间的耦合及离心力的影响,运用Hamilton原理建立旋转叶片子系统的刚柔耦合非线性动力学方程,采用有限元法和Newmark数值积分法,对旋转叶片在时变载荷作用下的动态响应进行求解分析.分析表明,在旋转过程中,由于叶片的旋转运动与弹性变形间的耦合作用,离心力对叶片振动动态特性的影响被弱化；考虑剪切变形影响时得到的叶片振动幅值比不考虑剪切变形影响时平均增大了7.5%.由此可见,剪切变形引起的叶片附加位移、叶片大范围旋转运动与弹性变形的耦合及离心力作

用都对叶片的结构响应有较大的影响,为了更精确地预测叶片的变形和载荷,在仿真过程中应予以充分考虑.分析过程还表明,运用柔性多体系统动力学理论中的刚柔耦合动力学建模方法，以及Hamilton原理分析风力机叶片的动力学问题的方法是可行的.

【相关文献】

[1] LI D Y,YE Z Q,CHEN Y.Multi-body dynamics numerical analysis of rotating blade of horizontal axis wind turbine [J].Acta Energiae Solaris Sinica,2005,26(4):473-481.

[2] 杨从新,李 昆,季 炜.风力机叶片断面的几何参数分析 [J].兰州理工大学学报,2008,34(4):46-49. [3] HANSEN M H.Aeroelastic instability problems for wind turbines [J].Wind Energy,2007,10(6):551-577.

[4] LARSEN J W,NIELSEN S R K.Non-linear dynamics of wind turbine wings [J].International Journal of Non-Linear Mechanics,2006,41(5):629-643.

[5] 杨从新,张 强.考虑静气弹特性的风力机叶片优化设计 [J].兰州理工大学学报,2014,40(1):45-49. [6] YOUNSI R,EI-BATANONY I,TRISTCH J B,et al.Dynamic study of a wind turbine blade with horizontal axis [J].European Journal of Mechanics-A/Solids,2001,20(2):241-252. [7] LEE J W,LEE J S,HAN J H,et al.Aeroelastic analysis of wind turbine blades based on modified strip theory [J].Journal of Wind Engineering and Industrial Aerodynamics,2012,110:62-69.

[8] 洪嘉振,刘铸永.刚柔耦合动力学的建模方法 [J].上海交通大学学报,2009,42(11):1922-1926. [9] 陆佑方.柔性多体系统动力学 [M].北京:高等教育出版社,1996.

[10] 李本立,宋宪耕,贺德馨.风力机结构动力学 [M].北京:北京航空航天大学出版社,1999. [11] LANZAFAME R,MESSINA M.Fluid dynamics wind turbine design:critical analysis,optimization and application of BEM theory [J].Renewable Energy,2007,32(14):2291-2305.

[12] 王勖成.有限单元法 [M].北京; 清华大学出版社,2003.

[13] 刘 雄,李钢强,陈 严,等.水平轴风力机叶片动态响应分析 [J].机械工程学报,2010,46(12):128-134.

[14] 信伟平.风力机旋转叶片动力特性及响应分析 [D].汕头:汕头大学,2005.