近年来,“能源危机”越来越引起人们的重视,能源短缺使得可再生能源得到空前发展。风能作为取之不尽用之不竭的可再生能源在近几年得到了迅速发展,上不少都把开发利用风能作为一项能源政策。风机叶片是风力机的关键部件之一,目前大型风机叶片的材料主要是轻质高强、耐腐蚀性好、具有可设计性的复合材料,由于叶片采用复合材料铺层设计，结构异常复杂,单纯的经典理论解析计算已难以精确计算出叶片的强度和刚度,为此需要进行有限元的仿真模拟。本文采用ANSYS的参数化语言APDL直接建模,然后赋材料属性、划分网格,进行叶片模态分析,较好地模拟了叶片的真实结构,计算了叶片在自由状态下的固有频率和在20RPM下的预应力频率。后加载极限载荷校核了叶片各截面稳定性。
    1　叶片模型的建立
    1.1　建立几何模型
    叶片截面的翼型数据通过CATIA导出,结合弦长和扭转角计算出实际叶片截面的坐标。在AN-SYS的程序中形成如下格式:

    k,, 838. 309405 , -83. 92648 , 0
    k,, 771. 63901 , -338. 19972 , 0

    根据1MW叶片翼型的特点,将叶片分为45个截面,每个截面上有86个关键点。通过Bspline命令将每个翼型上的关键点连成18条曲线,然后将叶片翼型上的曲线通过纵向直线连接。后通过Askin命令建立曲面,每两个翼型截面之间就有18个曲面,建完所有曲面就生成了叶片蒙皮的几何模型。如图1(a)所示,再布置上主梁、腹板,形成整体叶片的几何模型。

    1.2　建立有限元模型
    在单元类型的选择上,根据叶片特点,主要采用shell91和shell99单元,其中shell91单元用于模拟夹芯结构。在定义材料性能参数时,主要采用定义实常数的方式来模拟材料的性能、铺层角和铺层厚度。之后选择合适的单元尺寸进行网格划分,终形成叶片的有限元模型。该模型共有21295个节点,共划分了7414个高精度的壳单元。如图1(b)所示。

    2　模态分析
    2.1　理论依据
    通用运动方程为:
    [M]{u..}+[C]{u.}+[K]{u}={F(t)} (1)
    假定自由振动并忽略阻尼,则:
    [M]{u..}+[K]{u}=0(2)
    其中, [M ]、[K]分别为叶片的质量和刚度矩阵; {u}为节点位移向量。
    式(2)的解为如下简谐运动:
    {u}={U}sinωt (3)
    式中, {U}为模态形状(无量纲位移);ω为圆频率。
    式(3)代入式(2)得:
    ([K]-ω2[M]){U}={0} (4)
    方程(4)中{U}要有非0解,系数行列式为0,得:
    det([K]-λ[M])=0(5)
    其中,λ=ω2,上述行列式为λ的多项式,有根λ1,λ2,…,λn。
    代入λi得方程:
    ([K]-λi[M]){Ui}={0}, i=1,2,…,n (6)
    可以求得{Ui},即模态。fi=ωi/2π为系统固有频率
    2.2　计算结果
    2.2.1　无预应力自然模态分析

 [-page-]

    2.2.2　预应力模态分析
   假设叶片在20RPM的转速下工作,计算此时叶片的前10阶模态,见表1。

    2.2.3　小　结
    通过对1MW叶片模态计算,结论如下:
    (1)阶振动形式主要为挥舞,第二阶振动形式主要为摆振,这说明弯曲振动(包括挥舞和摆振)频率较低,是风力机叶片的主要振动形式。叶片低阶固有频率f1=0・987Hz=59・22rpm,远大于叶片启动和正常工作转速(12~21・5rpm),因此启动过程与正常工作时叶片不会出现共振现象;
    (2)从前10阶模态计算结果可见,前5阶模态都为弯曲振动,仅第6阶和第10阶出现扭转振动,说明叶片抗扭转振动的能力较强;
    (3)风力机叶片在额定转速时的模态频率比静止时的固有频率要高,这是因为离心力引起动力刚化,从而导致叶片的模态频率增加。
    3　稳定性分析
    3.1　理论依据
    屈曲稳定性分析是在结构的线性刚度矩阵中引入微分刚度的影响。微分刚度是从应变-位移关系式中的高阶项导出的。设结构线性刚度矩阵为[Ka],考虑应变-位移的高阶非线性项的微分刚度矩阵为[Kd],一般[Kd]与所施加载荷Pa成比例,即:
    [Kd]=Pa[-Kd] (7)
    则结构的总刚度矩阵为:
    [K]=[Ka]+[Kd] (8)
    总应变能为:
    U=12{X}T[Ka]{X}+12{X}T[Kd]{X} (9)
    其中, {X}为各节点的位移向量。为使系统达到静力平衡,总应变能必须有一个驻值,即:
    UX=[Ka]{X}+[Kd]{X}={0} (10)
    将方程(7)代入方程(10)可得:
    [Ka]+Pa[-Kd] ){X}={0} (11)
    为使方程(11)有非0解,则方程(11)的系数行列式为0,因此:
    det([Ka]+Pa[-Kd])=0 (12)
    方程(12)只有对特定的Pa才成立,这样的Pa称为临界屈曲载荷Pcr,记:λi=PcriP(13)
    则方程(12)可以表示为:
    det([Ka]+λi[-Kd])=0 (14)
    可见,求解屈曲临界载荷Pcri转化为求解特征值问题即式(14),所求屈曲临界载荷为:
    Pcr=min(λi)Pa(15)
    min(λi)为失稳临界特征值,又称为失稳屈曲因子,为失稳临界载荷与设计载荷之比。

    3.2　计算结果

    3.3　小　结
    分析表明:
    (1)载荷为Mx(min)极限载荷的1.01倍时,在4m截面处先出现失稳
    (2)载荷为Mx(max)极限载荷的1.34倍时,在4m截面处先出现失稳;
    (3)载荷为My(min)极限载荷的1.25倍时,在4.5m截面处先出现失稳;
    (4)载荷为My(max)极限载荷的1.11倍时,在4.5m截面处先出现失稳;
    (5)在My(max),My(min),Mx(max),Mx(min)四种工况下屈曲因子都大于1,说明这时在给定的安全系数和设计载荷下,结构不会失稳。 [-page-] 
    4　结　论
    本文主要介绍通过ANSYS建模技术对叶片进行有限元分析。通过采用ANSYS参数化语言APDL建立风力机叶片的几何模型,通过实常数赋值来实现对叶片铺层的模拟,较为真实地模拟了叶片的实际结构。在剖分网格、建立叶片有限元模型之后,对叶片进行无预应力模态分析和预应力模态分析,计算得出了叶片前10阶的频率,分析得出叶片在启动和正常工作时不会出现共振。后分析了叶片在极限载荷下的整体稳定性,计算表明,在My(max),My(min),Mx(max),Mx(min)四种工况下屈曲因子都大于1,说明结构稳定,不会发生失稳。(赵　娜,李军向,李成良
)