风力机复合材料柔性叶片的颤振分析

任勇生，杜向红，杨树莲

(1.山东科技大学 机械电子工程学院，山东青岛 266510;2.山东工商学院，山东烟台 264005)

为了最大限度地提高发电功率，降低发电成本，现代风力发电机叶片的尺寸显著增大。随着风力机叶片长度的增加以及复合材料在叶片制作中的广泛使用，叶片的刚度也明显降低，表现为更加柔性化的趋势。在风载荷的作用下，叶片的运动会影响周围风场从而改变作用于叶片的风载荷，存在着空气动力、惯性力和弹性力的相互作用，因此，在叶片模态之间会发生振动耦合，形成气动弹性不稳定——颤振。叶片的颤振问题属于气动弹性力学的研究范畴，是现代风力机系统设计必须解决的主要隐患之一。近年来，在气动弹性力学的理论框架内揭示叶片颤振行为发生的机理，并用于指导风力机的气动弹性设计，已经成为风力机结构动力学研究中的一个新的重要内容，日益受到风能领域的重视［1］。

在风力机叶片的颤振研究中，基于连续分布参数模型与典型截面集中参数模型相比，由于能够反映叶片的剖面、变形以及惯性沿叶片展向的变化，因此可以更好地描述叶片的颤振行为。经典的叶片分布参数模型是具有弯扭耦合变形梁的偏微分方程组，称之为Hodges-Dowell方程［2］。文献［3］基于 Hodges-Dowell方程和准定常叶素气动力模型，取单模态近似表示叶片的位移，采用特征值法分析叶片挥舞/摆振/扭转的稳定性;文献［4］引入重力、调距作用和转速变化等因素的影响，将Hodges-Dowell方程进行了扩展;文献［5］忽略叶片截面翘曲的影响，基于Euler梁理论并且采用模态叠加法，研究准定常气动力作用下的风力机叶片的受迫振动问题。但是，在上述研究中，叶片均被视为由各向同性材料构成的实心梁。事实上，典型的风力机叶片具有复合材料空心结构特点，通常需要采用薄壁闭室复合材料梁模型进行描述，它的结构动力学建模，要求恰当地表达由材料各向异性引起挥舞、摆振和扭转变形之间的弹性耦合和横截面翘曲等。文献［6］虽然涉及薄壁复合材料叶片结构的气弹稳定性问题，但没有给出具体的求解过程;文献［7］研究直升机旋翼复合材料薄壁叶片气弹稳定性问题，但所采用的是谐波平衡法(HBM)，计算过程比较复杂。

基于叶片连续分布参数模型的颤振问题研究，通常涉及高维偏微分方程组的求解，一般难以得到问题的精确解。虽然有限元方法作为一种数值近似方法，在复杂结构的力学建模已经得到广泛的应用，能够提供各向异性复合材料叶片的结构模型的详细的解答，但是由于计算成本较高，对于叶片气动弹性问题的初步设计，缺乏实用性。而振型叠加法利用叶片的低阶主要模态实现叶片颤振系统的模型降阶，能够极大地减少计算量，是风力机气动弹性力学研究中经常采用的一种有效的近似计算方法。

本文在现有研究结果的基础上，提出一个风力机叶片气动弹性力学理论模型，叶片结构采用具有预扭转的弯扭耦合闭合截面薄壁复合材料梁模型进行描述，气动载荷采用叶素动量理论和准定常气动力理论进行描述。基于Hamilton原理并结合VAM［8］，进行横截面分析，导出旋转复合材料叶片的结构力学模型，在此基础上，与叶素动量理论和准定常气动力理论进行组合，建立叶片的颤振分析模型。将位移按广义坐标进行模态展开，采用Galerkin法求解颤振分析模型。借助于特征值技术，对CAS构型叶片的颤振性能进行数值近似计算，揭示了入流比、预扭转角和纤维铺层角对风力机叶片颤振性能的影响。

1 基本理论

1.1 位移场和应变场

图1表示具有任意截面形状、长度为R的复合材料叶片，其变形前的轴线x相对旋转平面倾斜小角度βp，称之为预锥角，叶片以不变角速度Ω绕转轴旋转。为便于分析，建立三个坐标系:叶片旋转坐标系(x，y，z)，坐标原点o位于固定端;叶片局部坐标系(x，η，ζ)，η和ζ是叶片横截面的惯性主轴;叶片局部坐标系(x，s，ξ)，s沿截面中线(即，薄壁叶片的中面与横截面的交线)切向，ξ沿截面中线外法线方向。变形前的叶片截面，如图2所示，假定截面质心与弹性中心重合，弹性中心与气动中心的距离为xA，β是叶片横截面的预扭转角。

叶片坐标系和主轴坐标系之间，存在下列关系:

其中:β=β0(x/R)2;β0表示叶片尖端的预扭转角。

叶片上的任意点沿坐标系(x，y，z)的三个坐标轴方向的位移分量为［9］:

其中:截面翘曲函数g(s，x)由位移场的连续性、环向力为零以及定常剪力的条件决定。

与位移场(2)相对应的二阶近似应变场为:

其中:G、g1、g2、g3分别是与扭转、轴向拉伸、绕z轴弯曲和y轴弯曲相关的翘曲分量;()'表示对x求偏导。

1.2 运动方程

为了导出叶片的气动弹性方程，利用Hamilton原理建立叶片运动方程:

其中，U和T分别为应变能和动能，可由下式确定:

外力的虚功:

其中:Fx，Fy，Fz分别是作用于x，y，z轴方向的分布气动力，Mx是绕弹性轴的气动力矩。

其中:

叶片截面内力(矩)可以表示如下:

其中cij表示叶片截面刚度系数，它们与横截面主轴的刚度系数kij之间满足下列关系:

式中，kij的表达式见文献［9］。

将式(11)的第一式代入方程组(7)的第一个方程，假定沿叶片轴向的气动力Fx=0，从x到R做定积分，得:

如果将方程式(8)、(9)和(12)代入方程组(7)，令Ω=0和β=0，且不计外载荷的作用，则可以导出与文献［9］中的方程(27)一致的结果。

将式(11)和式(12)代入方程组(7)的后三个方程，得:

在本文的算例中，仅讨论具有CAS构型的薄壁复合材料叶片的情形，此时有:θ(y)=- θ(-y)，θ(z)=-θ(-z)，其中θ表示由正的s轴进行度量的纤维铺层角，于是得:k12=k13=k14=k24=k34=0，由此可以推出:C12=C13=C14=0，但是本文模型引入了叶片的预扭转角β，故C24=C314≠0。因此可以看出，本文模型比文献［9］的结果更具有普遍性。

1.3 气动载荷

根据文献［3］的建议，基于准定常叶素理论的叶片的分布气动力和力矩的表达式，可以导出如下:

其中:ρA为空气密度;a为升力曲线斜率;b为无量纲半弦长=xA/R;αe=β+tan-1(Up/UT)+φ为有效攻角;CD0为形状阻力系数;Up、UT分别为叶片弹性轴沿挥舞和摆振方向的相对速度分量;是叶片截面位置x和时间t的函数。不计阵风并且略去变形高阶量，则有:

其中:λ=(Vm-vi)/(RΩ)为入流比，Vm和vi分别为平均风速和诱导速度。

1.4 求解方法

其中:下标s和d分别为叶片的静力和动力分量。

将式(14)、式(15)和式(16)代入方程组(13)，导出气弹方程组如下:

为了对气弹方程组(17)进行无量纲化，引入下列变换:

叶片的弯曲变形和扭转变形可以表示如下:

根据Galerkin近似求解方法，利用假设振型(19)消除薄壁梁的空间位置变量，并结合变换关系(18)，可以将叶片气动弹性方程组(17)，转化为关于时间的常微分方程如下:

其中:［M］，［C］，［K］分别表示关于广义坐标{X}={Vj，Wj，φj}T的质量矩阵、阻尼矩阵和刚度矩阵，为了节省篇幅，有关这些矩阵元素的具体计算公式，不再列出。

叶片的气弹稳定性取决于下列6N×6N矩阵A的特征值

在叶片的其它参数给定的情况下，矩阵A的特征值可以视为λ、β0和θ的函数。显然，当特征值的实部小于零，对应的振动模态是稳定的;当特征值的实部大于零，对应的振动模态是不稳定的;而实部等于零，对应的振动模态是临界稳定的。

2 计算结果与讨论

为了验证本文建立的模型及其近似计算方法的正确性，图3给出CAS构型复合材料箱型截面薄壁悬臂梁固有频率的近似计算结果，并且与文献［10］的固有频率的精确解进行了比较，结构的几何参数和材料参数均取自文献［10］。可以看出，二者符合得很好。在计算过程中，我们发现，选取模态项数N=5，即可得到收敛的第一阶ω1、二阶挥舞为主的固有频率ω3和第一阶扭转为主的固有频率ω5。本文计算结果与实验结果的对比，如表1所示。其中的实验数据取自文献［11］中的图4，图6，图7，图10和图11，实验用悬臂复合材料箱型梁的几何尺寸、铺层方式和材料特性，详见文献［11］的表1。

在下面的叶片颤振数值计算中，选取叶片的几何尺寸为:截面宽度0.95 m，高度0.4 m，薄壁厚度h=0.047 5 m，单层厚度h/6，铺层数为6层;叶片长度R=19 m。对于CAS构型叶片，所采用铺层方式为:上壁［θ］6，下壁［-θ］6，左壁和右壁［θ/-θ］3。叶片复合材料选取石墨/环氧树脂，其性能参数为:E11=142 GPa，E22=E33=9.8 GPa，G12=G13=6.0 GPa，G23=4.83 GPa，ν12= ν13=0.42，ν23=0.50，ρ=1 601.1 kg/m3。

叶片其他参数取为:

表1 与文献［11］固有频率的实验结果的比较Tab.1 Comparison of frequencies(Hz)with experimental results from Ref.［11］

图3 CAS构型悬臂梁的第一、二阶挥舞为主和第一阶扭转为主固有频率随铺层角变化规律Fig.3 Natural frequencies vs.ply angle CAS cantilevered beam

图4表示叶片颤振性能随入流比λ的变化规律，并且还显示了叶片预扭转角的影响。其中，图4(a)、图4(c)表示叶片前两阶挥舞模态的情形(ζi，i=1，3)，图4(b)、图4(e)表示叶片前两阶摆振模态的情形(ζi，i=2，5)，图4(d)、图4(f)表示叶片前两阶扭转模态的情形(ζi，i=4，6)。由图 4 可见，当预扭转角 β0为零时，在所给定的入流比参数变化范围内，由于对应的特征值的实部位于零附近，故叶片的一阶和二阶挥舞模态和摆振模态都是临界稳定的;而前两阶扭转模态由于对应的特征值实部小于零，因此是严格稳定的。这说明，作用于叶片扭转模态的气弹阻尼相对较大，具有较强的抑制扭转模态颤振的能力。当β0增加时，叶片挥舞模态明显地由临界稳定变为严格稳定的，β0的存在对挥舞模态的稳定效果非常明显，而摆振模态则由临界稳定变为不稳定的，并且其不稳定性随着β0的增加而增加，这很可能是由于作用于摆振方向上的气动阻力非常小，并且同时也没有考虑结构阻尼的缘故;与此同时，还可以看到，β0对扭转模态的颤振性能的作用却非常微弱，不会使扭转模态的颤振性能产生本质的变化。

图5表示叶片颤振性能随纤维铺层角θ的变化规律，并且也显示了叶片预扭转角的影响。由图5可见，总体而言，纤维铺层角对挥舞模态和摆振模态的颤振的影响较大，而对扭转模态颤振的影响较小。在预扭转角不为零时，增加铺层角对第一阶挥舞模态和摆振模态的稳定性，分别产生增强和削弱的作用。除了在铺层角为10°附近，对第二阶挥舞模态的稳定性有所削弱之外，在整个铺层角变化范围内，不会带来明显的影响，而铺层角的变化对第二阶摆振模态的稳定性，几乎不会产生影响。

图4 入流比对叶片的稳定性的影响(取不同的预扭转角)Fig.4 Effect of inflow on blade stability for different pre-twist angle

图6表示叶片颤振性能随预扭转角β0的变化规律，同时也显示了叶片纤维铺层角的影响。由图6可以清楚地看到，纤维铺层角对叶片的高阶挥舞与摆振模态的颤振性能的影响很小，而预扭转角对叶片扭转模态的颤振性能的影响很小。这个结论与图5的结果是相互对应的。

3 结论

本文基于弯扭耦合闭合截面薄壁复合材料梁的分析理论，采用Hamilton原理并且借助于Galerkin法，建立了在准定常气动力作用下的风力机叶片颤振分析模型。研究结果表明，采用本文建立的模型及其近似计算方法所获得的叶片的固有振动频率数值近似解，与已有文献的精确解相当一致。针对一个长度为18 m的石墨/环氧树脂复合材料叶片，在转速 Ω=4π/3(rads-1)的情况下的颤振分析结果表明:

(1)利用所建立的模型能够对风力机叶片的颤振性能进行参数分析，以确定影响叶片颤振性能的主要因素，为进一步改善或提高叶片的气动弹性性能，提供有用的信息。

(2)叶片的预扭转角对叶片的颤振性能会产生显著的影响，但对于挥舞模态和摆振模态来说，预扭转角作用效果是相反的。在实际应用中，预扭转角是影响风力机叶片颤振性能的一个主要设计参数。以往的研究表明，叶片的预扭转角可以借助于形状记忆合金和压电材料的驱动而产生，因此，上述研究结果对于客观地评价智能材料在叶片颤振抑制中的应用效果，有一定的参考价值。

［1］ Hansen M O L ，Sorensen J N，Voutsinas S，et al.State of the art in wind turbine aerodynamics and aeroelasticity［J］.Progress in Aerospace Sciences，2006，42(4):285-330.

［2］Hodges D H，Dowell E H.Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades［R］.Technical Report TN D-7818，NASA，December 1974.

［3］Kottapalli S B R，Friedmann P P.Aeroelastic stability and response of horizontal axis wind turbine blades［J］.AIAA Journal，1979，17(12):1381-1389.

［4］ Kallesoe B S.Equations of motion for a rotor blade，including gravity，pitch action and rotor speed variations［J］.Wind Energy，2007，10:209-230.

［5］Kim K T，Lee C W.Forced vibration analysis of flexible wind turbine blades by using assumed modes method［C］.Proc.of 6thInt.Vibration Engineering Conference，2008，192-201.

［6］Wilkie W K，Belvin W K.Aeroelastic analysis of helicopter rotorblades incorporating anisotropic piezoelectric twist actuation［C］.Technical Report:NASA-96-asme-wkw Year of Publication:1996.

［7］Kim T Y.Nonlinear large amplitude structural and aeroelastic behavior of composite rotor blades at large static deflection［D］.Dorctor Dissertation， Massachusetts Institute of Technology，1992.

［8］Berdichevsky V，Armanios E A，Badir A M.Theory of anisotropic thin-walled closed-cross-section beams［J］.Composite Engineering(Special Issue)，1992，2(5-7):411-432.

［9］Armanios E A，Badir A M.Free vibration analysis of anisotropic thin-walledclosed-section beams［J］.AIAA J，1995，33(10):1905-1910.

［10］ Dancila D S，Armanios E A.The influence of coupling on the free vibration of anisotropic thin-walled closed-section beams［J］.Int J Solids Struct，1998，35(23):3105-3119.

［11］ Chandra R，Chopra I.Experimental theoretical investigation of the vibration characteristics of rotating composite box beams［J］.J Aircraft，1992，29(4):657-664.