王旭东 陈 进 王立存 陆群峰

1．重庆工商大学废油资源化技术与装备教育部工程研究中心，重庆，400067 2．重庆大学机械传动国家重点实验室，重庆，400030

0 引言

风力机系统作为一个多自由度系统，在对其进行动力学特性研究时就必须研究系统的稳定性。在进行风力机系统设计时，为了保证这个多自由度系统的稳定性，首先需要对每个部件的载荷和稳定性进行研究［1－2］。

风力机叶片作为风力机的关键部件，其良好的设计、可靠的质量和优越的性能是保证机组正常稳定运行的决定因素。随着风力机向大功率发展，叶片的载荷及振动偏移特性的研究就更突出其重要性，而当风力机塔架的刚度很大时，风力机系统的动力学特性主要取决于风力机叶片动力学特性，可以说，叶片作为风力机的主要部件，其振动偏移和稳定性研究直接决定了风力机系统设计的可靠性［3－5］。因此，对风力机叶片的振动偏移和稳定性进行研究就显得尤为必要。本文从结构动力学理论出发，基于风力机旋转风轮的气动弹性模型，提出了旋转叶片偏移变形特性的数值模拟计算理论和方法。

1 旋转风轮的载荷特性模型

作用在风轮上的空气动力学载荷包括各个叶片拍打方向所承受的推力和摆振方向所承受的剪力，推力和剪力的计算采用动量叶素理论，利用Prandtl叶尖损失修正因子和Shen叶尖损失修正模型，对计算过程中叶片每个叶素中的法向力和切向力系数 （Cn，Ct）进行修正［6］。那么，叶片上推力和剪力的计算公式就可以表示为

它所产生的力矩的计算公式为

式中，B为叶片个数；λ为叶尖速比；c为叶片弦长。

当风轮旋转时，将会对叶片产生离心力，并且方向沿叶片向外。同时，在风轮运转的过程中，作用在叶片上的空气动力学载荷所产生的推力使得叶片偏离旋转平面，进而使得作用在叶片上的离心力方向发生了改变，如图1所示。此时叶片所承受的离心力可以表示为

式中，R为风轮半径；Ω为风轮转速；mi为风轮叶素质量。

图1 叶片的离心力示意图

从图1可以看出，叶片离心力在拍打方向的分量为

式中，δ为叶片偏离旋转平面的角度；θ为风轮的锥角。

可以看出，由于叶片离心力在拍打方向的分量Fcn与叶片所承受的推力Tx方向相反，因此，该离心力的存在能够减小拍打弯矩带来的叶片偏离距离。

2 风轮叶片偏移变形的动力学模型

根据虚功原理，基于结构动力学理论建立叶片的动力学方程：

式中，M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；Fg为外载荷，即空气动力学载荷和离心力；x为位移向量。考虑叶片前三阶振动模态的影响，取11个自由度，即

式中，w为风轮弹性位移；φ为叶片的角度位置；xji为叶片i在第j阶模态时所产生的位移；i，j＝1，2；上标f、e分别表示挥舞方向和摆动方向。

通过对式（7）的迭代求解得到式（8）中叶片在各自由度方向的位移。对于叶片在外载荷作用下产生的位移、速度和加速度，采用叶片前三阶振动模态的线性关系表达为［7］

式中，d为叶片产生的位移；u1f为一阶挥舞方向模态单位特征值；u1e为一阶摆振方向模态单位特征值；u2f为二阶挥舞方向模态单位特征值。

3 旋转风轮叶片载荷的数值模拟

根据上述建立的风轮载荷计算模型，在MATLAB软件中采用数值迭代法对某2MW风力机风轮的载荷进行了计算分析。该风轮由3个叶片组成，每个叶片长30．56m，设计额定风速为14m／s，叶尖速比λ＝5，风轮的设计旋转速度为2．3rad／s。以叶素理论为基础，将叶片沿展向划分成20段，通过对风轮诱导因子的迭代求解，确定风轮各个叶片的空气动力学载荷［8］。图2所示为额定风速14m／s，风轮旋转12s时，即旋转了4．43r后，所求得的叶片轴向和周向诱导因子在叶片展向的分布。图3所示为该时刻叶片的空气动力学载荷法向力和切向力分布，其中切向力主要表现为风轮的输出功率，法向力主要表现为风轮所承受的推力。风轮在该推力的作用下，在垂直旋转平面的方向会发生拍打位移，该位移的产生使得叶片的离心力方向发生变化，进而使得叶片的离心力在拍打和摆振方向产生了分量，由于拍打方向的离心力和推力方向相反，因此使得叶片的推力减小，图4所示为相同时刻条件下计算得到的叶片离心力在拍打和摆振方向的分布。

图2 叶片诱导因子沿展向的分布

图3 叶片空气动力学载荷的分布

图4 叶片离心力的分布

4 风轮偏移变形特性的算例分析

为了很好地研究风力机叶片的偏移变形特性和验证该模型的可靠性，对风轮在不同风速和不同位置时刻的叶片的偏移变形特性进行了计算分析。同样针对上述2MW风力机，在已得的载荷基础上，在MATLAB软件中编写程序对叶片在该载荷下的偏移变形位移进行了数值模拟。叶片在载荷作用下的偏移变形与其所处的运转位置有关，本文以图5所示各叶片所处的实线位置为风轮的起始旋转位置，风轮以图中箭头方向即沿逆时针方向旋转。图6所示为风轮在额定风速14m／s下，运转8s时，风轮三个叶片沿展向的不同位移偏移变形。因运转时间为8s，即刚好旋转了3r，那么该时刻各叶片所处的位置和起始位置相同。图7所示为相同风速下，风轮在旋转12s以后，即旋转了4．43r，各个叶片处于图5中虚线位置时，叶片的展向偏移变形。和图6对比，尽管来流风速和风轮旋转速度相同，但是由于风轮旋转时间不同，各叶片计算时刻所处的位置不一样，因此，三个叶片的偏移变形和图6所示的偏移变形有着不同。图8所示为该风轮运转12s时，风速为10m／s下，风轮各个叶片的偏移变形分布。可以看出，图8中每个叶片的偏移变形规律和图7一致，这是由于两种工况的运转时间相同，即计算时刻每个叶片所处的位置相同，但是由于图8风速为10m／s，比图7中风速14m／s要小，即图8叶片所承受的载荷要小些，因此叶片展向所产生的偏移变形也相对较小。

图5 风轮运转时叶片所处的位置

图6 额定风速14m／s下风轮旋转3转时各个叶片的偏移变形分布

图7 额定风速14m／s下风轮旋转4．43转时各个叶片的偏移变形分布

图8 额定风速10m／s下风轮旋转4．43转时各个叶片的偏移变形分布

5 结论

（1）应用动量叶素理论对风力机风轮旋转条件下的空气动力学载荷和离心力进行了理论模型的建立和数值模拟。在充分考虑该载荷特性条件下，基于风轮的结构动力学方程，提出了叶片的偏移变形特性计算数学模型。

（2）针对建立的叶片偏移变形计算数学模型，采用迭代法对某2MW风力机在不同工况下的叶片展向不同位置的偏移变形进行了数值模拟，得到了该风轮在不同风速和不同运转时刻的各个叶片的偏移变形分布。研究结果对风力机叶片的振动特性和疲劳寿命预测具有理论指导意义。

［1］ 贺德馨．风工程与工业空气动力学［M］．北京：国防工业出版社，2006．

［2］ Maheri A，Noroozi S，Vinney J．Combined Analytical／FEA－based Coupled Aero Structure Simulation of a Wind Turbine with Bend-twist Adaptive Blades［J］．Renewable Energy，2007，32：916－930．

［3］ Larsen J W，Nielsen S R K．Nonlinear Parametric Instability of Wind Turbine Wings［J］．Journal of Sound and Vibration，2007，299：64－82．

［4］ Sarkar S，Bijl H．Nonlinear Aeroelastic Behavior of an Oscillating Airfoil During Stall－induced Vibration［J］．Journal of Fluids and Structures，2008，24：757－777．

［5］ Turhan O，Bulut G．On Nonlinear Vibrations of a Rotating Beam［J］．Journal of Sound and Vibration，2009，322：314－335．

［6］ Shen W Z，Mikkelsen R，Sorensen J N．Tip Loss Correction for Wind Turbines Computations［J］．Wind Energy，2005，8：457－475．

［7］ Hansen M O L．Aerodynamics of Wind Turbines［M］．London：James ＆ James （Science Publishers）Ltd，2000．

［8］ 陈进，王旭东，沈文忠，等．风力机叶片的形状优化设计［J］．机械工程学报，2010，46（3）：131－134．