1 前言
    振动是工程结构中经常遇到的现象，风力机也不例外。当风力机在自然条件下运营时，作用在风力机叶片上的空气动力、惯性力和弹性力等交变荷载会使叶片结构产生振动，有可能影响风力机的正常运行，为防止灾难性振动的发生，在设计风力机时就必须进行动力学分析，包括对叶片结构自振特性的分析。
    风力机叶轮叶片的振动主要有三种振动形式：①挥舞，是叶片在垂直于旋转平面方向上的弯曲振动；②摆振，是叶片在旋转平面内的弯曲振动；③扭振，是绕叶片变距轴的扭转振动。此外，在空气动力、惯性力和弹性力的耦合作用下，这三种形式的振动还会发生耦合，引起气动弹性问题。一般风力机叶片气动弹性动态不稳定包括挥舞-扭转不稳定、扭转-摆振不稳定、失速颤振和失速诱导的摆振等^[1]。
    实际上作用在叶片上的荷载为动荷载，且具有交变性和随机性，当叶片的固有频率与激振力的频率相同时，就会产生共振，从而影响叶片的安全。由于动荷载很难用数学式来描述，所以工程上一般并不对叶片在动荷载作用下的动态响应进行动态分析，而是着重于分析计算叶片的自振频率，以防发生共振。一般叶片自振频率不与叶轮的转速频率的整数倍重合，就可以避免发生共振，实践中重要的是避免低阶自振频率与叶轮的转速频率的整数倍重合。因此，精确分析计算叶片结构的自振频率，对大型风力机叶片的设计具有重要的意义。
    关于风力机叶片自振频率的计算方法，常用的方法有瑞利法、迭代法、矩阵分析法和有限元法。对于迭代法和矩阵分析法^[2]，主要求解弯曲振动频率及扭转振动频率，不能够求解耦合振动频率。瑞利法求解一阶挥舞和摆振方向自振频率具有很高的精度，但计算叶片弯扭耦合振动频率方法不可靠，误差很大^[3]。有限元法求解叶片结构的自振特性是一种可靠而有效的方法，但是风力机叶片的结构设计过程往往是一个逐步调整结构布局和材料铺层方案反复计算的过程，再加上叶片结构的复杂性，三维有限元模型建模相当困难，数值分析的效率较低，因此在叶片的初步设计中很少直接应用三维壳单元或实体单元进行有限元分析，只是到详细设计阶段才会选用该法进行分析。
    由于叶片的剖面尺寸远较叶片的长度为小，故可以把叶片结构的计算模型简化为一根端固定的悬臂梁。本文拟把复杂的叶片结构简化为一端固定一端自由的悬臂梁进行结构动力特性分析，从而使创建叶片结构的有限元模型大大简化。同时，应用商业有限元软件对风力机叶片进行三维有限元分析，以证明此方法的可靠性。 [-page-]
2 计算结构固有频率和振型的有限元法
    由结构动力有限元可知，叶片总体结构有限元模型的运动方程为：

          
    其中，{u}为单元节点系位移向量；[C]为阻尼矩阵；[K]为总体刚度矩阵[M]为质量矩阵；｛P(t)｝为外荷载向量。
     对于无阻尼自由振动，可得叶片结构有限元运动方程为：

          
    由于结构的自由振动总可分解为一系列的简谐振动，把这些简谐振动相叠加即构成式(2)的解，因此可设式(2)的一系列简谐振动解中的第i个解为

        
    其中，ω_i和φ_i分别为第i个振型的固有频率和初相位角;A_i为对应于第i个振型的诸位移的振幅大值组成的振幅向量,A_i={A₁,A₂,…An}^T。
    将式(3)代入式(2)得：

        
    对于n个自由度的振动系统，振幅向量不会全部为0，即A_i≠{0}，因此式(4)中必有行列式：

        
    式(5)即为结构振动系统的特征方程，展开后得：

      
    根据上式即可求得各阶固有频率ω_i及对应的振幅比（即振型）。
3 采用梁单元分析叶片自振特性
    根据国内200kW风力机叶片技术资料^[3]，应用大型通用性有限元软件，按梁单元进行建模分析，200kW风力机叶片技术资料如下。
    图1为国产200kW风力机叶片示意图^[9]，由定叶和动叶两部分组成。定叶采用O型主梁和上、下薄壳组合式结构，在O型梁中间装有叶尖扰流器主轴，其结构形式见图2。动叶采用上、下薄壳组合式结构，中间装有叶尖扰流器轴套，其结构形式见图3。

              
    叶片技术参数为额定功率200kW;叶轮直径23m；起动风速4.5m/s；额定风速14m/s；关机风速28m/s；大抗风60m/s；重量800kg；叶片长10.8m；根部弦长1.445m；尖部弦长0.565m。翼型为NA-CA632XX^[3]。 [-page-] 
    按等代设计方法计算时等效材料参数选用弹性模量E＝20GPa；泊松比μ=0.14;剪切模量G=3.6GPa。
    叶片各剖面几何特性见表1。

    求得叶片挥舞方向及摆振方向前三阶自振频率，见表2。

             
4 采用壳单元分析叶片的自振特性
    根据上述200kW风力机叶片技术资料，采用壳单元分析计算上述该风力机叶片的自振特性。图4为有限元网格图，单元数为13886个，建模中根据各个剖面的翼型数据，采用三次样条曲线拟合，并根据各剖面O型梁的坐标确定出主梁位置，后分别由各边界曲线生成空间曲面。对叶片结构进行模态分析图，求得各振型及对应的频率，图5(a~f)分别为叶片的前六阶振型，由振型图可以看出，阶振型为挥舞方向的振动，第二阶振型为摆振方向的振动，第三阶为挥舞方向的二阶振型，第四阶第五阶振型为挥舞-扭转耦合振动，第六阶振型为摆振-扭转耦合振动。表3为计算所得到的各方向振动所对应的频率。

              
5 分析与结论
    表4为三种方法计算结果的比较。
    从表4可以看出，采用壳单元分析得出的叶片前三阶自振频率与采用梁单元分析得到的结果符合得很好，究其原因从叶片振型看，发现叶片的前三阶振动表现为纯弯曲振动。采用瑞利法得到的阶自振频率也具有足够的精度，第二阶自振频率误差稍大。可见在计算叶片的前三阶弯曲自振频率时，完全可以把复杂叶片结构简化为一悬臂梁并采用梁单元进行分析计算，甚至可以运用瑞利法计算叶片挥舞和摆振方向的一阶频率而不失良好的精度，没有必要在叶片的初步结构设计阶段即采用壳体单元分析叶片的自振特性。
    风力机叶片结构的初步设计是一个反复试算的过程，鉴于采用梁单元分析计算叶片的前三阶弯曲自振频率及振型具有足够的精度，且易于实现计算的程序化，从而使计算风力机叶片自振频率的过程大大简化，因此该方法对于大型水平轴风力机叶片的结构设计具有一定意义。
    诚然,进一步分析叶片的高阶频率及振型,按三维壳单元或实体单元进行建模分析计算的优势还是明显的。
                    参考文献
[1] 李本立，宋宪耕，贺德馨等.风力机结构动力学[M].北京:北京航空航天大学出版社,1999.
[2] 刘晶波,杜修力.结构动力学[M].北京:机械工业出版社,2005.
[3] 袁国青.200kW风力机叶片结构设计计算报告[R].上海:同济大学,1993．
[4] 陈云程，陈孝耀，朱成名.风力机设计与应用[M].上海:上海科学技术出版社,1990.
[5] 张锦南.200kW风力机叶片[J].太阳能,1997,3:23-24.
[6] Tony Burton，David Sharpe，Nick Jenkins，Ervin Bossanyi. Wind Energy Handbook[M].Tohn Wiley＆Sons Ltd,2001．
[7] H.Snel. Review of the Present Status of Rotor Aerodynamics[J]．Wind Energy,1998,1:46-69.