风力机筒形塔架结构静动态特性的有限元分析

1、风力机筒形塔架结构静动态特性的有限元分析陆萍黄珊秋张俊宋宪耕(山东工业大学机械工程学院,济南250061文摘:讨论了水平轴风力机筒形塔架结构的静动态分析建模和有限单元类型的选取及计算方法,并以200k W 风力机的六棱锥筒形塔架为例,计算并给出了静动态特性。关键词:风力机,锥筒形塔架,有限元法0引言水平轴风力机的塔架要承受风轮、发电机组和传动系统等重量,它直接影响风力机的工作可靠性。为确保机组的正常运行,要求合理地设计塔架的强度和刚度,并需分析计算塔架的静动态特性。本文以200k W 风力机的六棱锥筒形塔架为例,利用AD I NA 大型有限元程序完成了自动建模和静态与动态特性的分析计算。1塔架

2、建模200k W 风力发电机是水平轴的上风式并网风力发电设备,其主要组成部分如图1所示。塔架建模时,尽可能如实反映塔架结构的主要力学特性,并力求采用较简单模型,所作的简化假设是:(1由于装有发电机、变速器及制动器等的短仓刚度远大于塔架的刚度,故将短仓简化为作用在塔架顶端的集中质量和转动惯量,并可绕塔架中心轴转动。(2风轮简化为置于塔架顶部短仓悬臂端且可绕风轮水平中心轴转动的转子。(3塔架简化为底端固定、顶端自由的空间薄壁六棱锥筒形结构。(4塔筒采用壳单元模型,塔顶采用板单元模型。取塔脚横截面中心为原点的右手坐标系,如图2所示,其中x 轴为平行于风轮水平轴的外力方向。沿塔架高度等分为40层,以每

3、层的分界线与薄壁六棱锥的交点为节点,故该模型共有247个节点(含塔顶中心节点,三维壳单元240个,板单元6个。图2为由自动建模程序生成的塔架计算模型。第18卷第4期1997年10月太阳能学报A CTA EN ER G I A E SOLA R IS S I N I CA V o l 118,N o 14O ct .,1997本文1996210228收到 图1风力发电机结构简图图2计算模型2有限单元法的应用将风力机塔架作为大型的线弹性结构处理,要直接从基本方程求解往往很困难,而有限元法是解决此类问题的非常有效的数值分析方法。211基本方程假设塔架单元为线弹性体,并处于小变形范围内。由平衡、物理、

4、几何三方程可导出结构的有限元基本方程:M u ¨+D u +K u =R (1式中:K 是整体刚度矩阵;M 是整体质量矩阵;D 是整体阻尼矩阵;R 是外载荷列阵;u 、u u ¨分别为节点的位移、速度、加速度列阵。对静力问题,方程(1可简化为K u k =R (2对动力问题,在t +t 时刻有:t +t M u ¨k +t +t D u k +t +t K u k =t +t R 1(3上式左上角表示该量所在的时间状态,右上角表示在第K 次迭代中获得的状态矢量。以上方程都是大型联立方程,一旦解得位移u k ,利用u e =N u k ;e =B u k ;e =C

5、 k ,可以求得构件内部各点的位移、应变和应力。这里N 为单元形函数矩阵;B 为应变位移关系矩阵;C 为弹性矩阵;u e 、e 、e 分别为单元的位移、应变、应力列阵。212边界条件由于风力机塔架底部为固定端,理论上底部节点的自由度被约束,即底部节点的位移、转角均为零。213计算方法空间四节点壳单元的节点均取在中曲面上,如图3所示1。V _k 为节点K 处垂直于中曲面063太阳能学报18卷的向量。壳元中任一点K 的整体坐标系为x i =2nk =1N k x k i +t 22n k =1N k h k V k ni (n =4,i =13(4每个节点最多五个自由度(三个移动,二个转动,壳内任

6、意点位移的矩阵形式为u =N u k 。其中u k =u 11u 12u 13u 41u 42u 4344 T 图3壳的局部坐标形函数矩阵N 由插值函数N k ,壳厚度坐标t 及h k 和V k ni 构造。由整体坐标变换关系式(4可导出应变位移关系矩阵B ,而后分别求出单元刚度矩阵和单元质量矩阵;K e =V eB T C B d VM e =V e N T采用离散克希霍夫理论构造线性薄板单元1,其刚度矩阵为K e =A B T C bB d A 单元的一致质量矩阵为M e =V eN TN d V 式中:A 为板单元面积,C b 为弹性矩阵。单元所受的均布横向压力q 平均分配在板单元的三个

7、节点上,外载荷节点力列阵为R e =A 3q o o q o o q o o T 21313转换矩阵由于上述矩阵是在单元局部坐标系中生成的,而整个结构的组装是在整体坐标系中进行的,因此总装前要进行一次坐标转换,即K =A T K e A ,M =A T M e A ,R =A R e 其中,K 、M 、R 为整体坐标系中的值,A 为坐标转换矩阵。2.4方程求解方法对静力方程(2,采用平衡迭代法2求角位移和应力。对方程(3采用逐步积分的数值法求解,可对系统进行动态分析并求出响应,在逐步积分法中可采用W ilson 2法2求解方程(3。3静态分析风力机塔架终年在严酷的环境下工作,长期承受风载荷和重

8、力等负载的作用,其中起主要作用的有塔架自重,发电机组的重力,来自风轮的负载和风载荷。311风轮作用于塔架的负载图4所示为风轮作用于塔架负载的模型。来自风轮的所有负载(包括气动推力,偏转力矩,1634期陆萍等:风力机筒形塔架结构静动态特性的有限元分析 陀螺力矩等转换到塔架坐标系中,可分解为沿三个坐标轴的力和力矩,见图4。计算时分别按三种工况输入负载:(a 最大设计负载;(b 最大气动负载;(c 制动力矩。自动建模程序可由这些负载分别计算出塔架顶部各节点的固端力,加到塔架上。3.2风载荷塔架受到风压作用时结构会产生偏移和振动。风可分为稳定风和脉动风两种。为便于计算,工程上常用稳定风力乘以风振系数来

9、表达脉动风振动。作用于塔架结构上的风压计算公式3为W =z s z r w o kN m 2式中:z 为高度z 处的风振系数,s 为风载体型系数,z 为高度z 处的风压变化系数,r 为重现周期系数,w o 为基本风压。图5所示是在最大设计工况下计算得到的塔架综合位移变形图,虚线为原始状态,实线为变形后的状态 。图4负载计算模型图5静态位移变形图4动态分析在进行风力机塔架的动态分析时,需要计算塔架结构的固有频率和振型。由于结构阻尼对固有频率影响很小,故在求解时忽略阻尼的影响,同时令载荷项为零,可得结构的无阻尼自由振动方程:M u ¨+K u =0(5上式为常系数线性齐次微分方程式。假定

10、塔架上各节点的振动由不同频率的简谐振动组成,则其解的形式为u =u sin t ,代入方程(5得特征矩阵方程:(K -2M u =O (6上式为齐次线性代数方程组,若有非零解,则必有系数行列式等于零,即 K -2M =0。它是方程(5的特征方程。于是结构无阻尼自由振动的解就可化为特征方程的求解问题。263太阳能学报18卷 图6塔架各阶固有振型曲线对于该塔架的求解,采用子空间迭代法2。计算中除塔架自身分布质量外,将短仓质量视为集中质量。计算得到前五阶固有频率的结果见表1。对应的前五阶二维振型曲线如图6所示。3634期陆萍等:风力机筒形塔架结构静动态特性的有限元分析364 太阳能学报 表 1固有频

11、率计算结果 阶次 1 2 3 4 5 18 卷 固有频率 (H z 1. 46346 1. 46349 6. 84966 6. 85036 7. 42214 5结论 由上述论证分析与对六棱锥筒形塔架结构静动态特性的计算结果, 可得到如下结论: ( 1 在静态分析中, 塔架在最大设计负载工况下的位移变形量最大, 制动工况下的变形量 最小。 各工况下的最大应力值在塔架的底部, 节点的最大位移值在塔架顶端。 ( 2 由于六棱锥塔架结构基本上是对称的, 使塔架的固有频率值每两阶基本相等。 塔架的 一阶频率为一阶弯曲振动, 二阶频率为另一方向的一阶弯曲振动, 三阶频率为二阶弯曲振动, 四阶频率为另一方向

12、的二阶弯曲振动, 五阶频率为扭转振动。 ( 3 塔架的固有频率与组成塔架的材料、 几何尺寸及塔顶短仓质量等有关。 对于六棱锥筒 形塔架采用壳单元模型较为合理。 转动惯量对扭转固有频率影响较大, 不容忽略。 ( 4 计算结果表明, 一、 二阶固有频率小于周期激振频率, 三阶以后各阶都远离周期激振频 率, 因此理论上不会出现共振。 总之, 上述结论与用线单元梁方法给出的结论是定性一致的, 但本方法给出的结果更加接 近结构静动态特性的真实情况。 参考文献 1郑州机械研究所. AD I NA 理论手册和使用手册. 郑州: 郑州机械研究所, 1987 2 美 丁巴斯著. 工程分析中的有限元法. 北京:

13、机械工业出版社, 1991 K 3王肇民, U. Peil 塔桅结构. 上海: 同济大学出版社, 1989 . THE F IN ITE EL EM ENT ANALY S IS O N STAT IC AND DY NAM IC CHARACTER IST IC O F THE CO N ICAL TUBE TOW ER STRUCTURE FO R THE W IND TURB INE L u P ing H uang Shanqiu Zhang J un Song X iangeng (S hand ong U n iv ersity of T echnology , J inan 25

14、0061 Abstract: T h is p ap er d iscu sses the sta t ic and dynam ic ana lysis and m odelling of the con ica l tube tow er st ructu re fo r ho rizon ta l ax is w ind tu rb ine. It a lso d iscu sses the cho ice of the typ e of fin ite elem en t and ca lcu la t ing m ethod. A fo r ca lcu la t ion hexagona l con ica l tube tow er of 200