向 玲， 高 楠， 唐 亮， 郭鹏飞(华北电力大学 机械工程系, 河北 保定 071003)

风能作为一种资源丰富、分布广、无污染的绿色可再生资源，越来越受到世界各地的重视。风力发电机齿轮箱是风机的重要组成部分，但由于外部恶劣的工作环境，加上风力发电机齿轮传动系统本身传动比大、扭矩高等特点，导致风机齿轮箱频繁出现故障问题，降低了风力发电的效率。因此，风电齿轮传动系统的动力学研究，对降低齿轮传动系统故障有重要意义。常见的兆瓦级风力机齿轮箱是由行星齿轮传动和直齿轮传动组成的，行星齿轮传动具有结构紧凑、传动效率高、运动平稳和抗冲击能力强等优点，但行星齿轮传动系统含有内外同时啮合的齿轮，在内外激励下存在复杂的非线性因素，所以其非线性振动问题一直是研究的重点。卜忠红等[1]从模型、固有特性和动态相应计算等多个方面综述了近年来行星齿轮传动动力学的研究进展；邱星辉等[2]也对风力机组行星齿轮传动系统动力学的研究进行了综述；Kahrarman[3]建立了行星齿轮传动系统的线性纯扭转动力学模型，求解获得该系统的固有特性和振动响应；孙智民等[4]建立了2K-H行星齿轮传动系统的非线性动力学模型，研究了系统在不同参数条件下的简谐、非简谐单周期、次谐波、准周期和混沌稳态强迫响应，结果表明齿侧间隙参数影响系统的非线性动力学行为；张慧博等[5]研究了考虑多间隙耦合关系影响的齿轮副系统非线性动力学特性；孙涛等[6-7]采用谐波平衡法对行星齿轮传动系统的弯扭耦合非线性动力学进行了研究；李晟等[8]分析了激励频率和啮合阻尼比对两级行星齿轮传动系统分岔与混沌特性的影响；李同杰等[9]建立了考虑时变啮合刚度、综合啮合误差和齿侧间隙非线性因素下的行星齿轮传动系统纯扭转动力学模型，分析了转速、齿侧间隙和阻尼比对该系统动力学特性的影响；Zhao等[10]研究了静态传递误差、平均交变力比和时变啮合刚度变化对风电齿轮传动系统动力学特性的影响；Wang等[11]以风电齿轮传动系统的纯扭转非线性动力学模型为研究对象，对比并分析了不同阻尼比对该系统动力学特性的影响。众多学者对风电齿轮传动系统的动力学特性进行了研究，本论文在此基础上进一步研究风电齿轮传动系统扭转振动的分岔及混沌特性。

以1.5 MW双馈式风力发电机组的齿轮箱为研究对象，建立此齿轮传动系统的纯扭转振动非线性动力学模型，分析在内外激励变化下系统的非线性动力学特性。以时间历程图、相图、Poincaré截面图、FFT频谱图、全局分岔图及最大Lyapunov指数图作为分析手段，详细说明系统随着激励频率和综合啮合误差变化下的动力学特性。结果为风电齿轮传动系统的设计、减振和降噪提供理论参考。

1 非线性动力学模型

风电齿轮传动系统扭转动力学模型如图1所示，它是由一级行星齿轮传动加两级平行轴齿轮传动组成，所有齿轮为直齿轮。行星齿轮传动系统有太阳轮s，行星架c，行星轮pi(i=1,2,…,N)，内齿圈r；g1、g2为低速级齿轮；g3、g4为高速级齿轮，整体采用集中参数模型对系统建模。在行星齿轮传动系统中内齿圈固定，行星架c作为输入端，输入扭矩为Tin；高速级g4为输出端，输出扭矩为Tout。θm(m=c,s,pi,g1,g2,g3,g4)为行星架(各齿轮)的扭转角度；kj、cj、ej(j=spi,rpi,g1g2,g3g4)分别表示为各齿轮副间的时变啮合刚度、啮合阻尼及综合啮合误差。假设系统全部构件为刚体，不考虑各轴系扭转刚度和齿面摩擦。

在行星齿轮传动中，行星轮同时与太阳轮和内齿圈啮合传动，故各齿轮副啮合频率相同。在风电齿轮传动系统中，视内齿圈固定，可根据行星齿轮传动副传动关系推导啮合频率ωspi(i=1,2,…,N)可表示为

ωspi=ωrpi=ωczr

(1)

式中:ωc为行星架转速;zr为内齿圈齿数。

图1 系统非线性动力学模型

在整个系统中，行星齿轮级、低速级齿轮系和高速级齿轮系的齿轮副间具有不同的啮合频率，根据传动比关系，平行轴齿轮副的啮合频率是ωg1g2=Λ1ωspi、ωg3g4=Λ1Λ2ωspi，式中Λ1=5，Λ2=3.9分别为行星齿轮传动系和低速级齿轮系的传动比。

由于直齿轮啮合副刚度的特点，可采用周期矩形波表示啮合刚度[12]，其表达式可利用以啮合频率为基频的Fourier级数展开表示，取一次谐波项即

多增氧：溶氧决定产量，多开增氧机增氧，曝气、调水的效果好。使用罗茨风机加纳米管或纳米盘增氧，功率可达10瓦/立方水体，产量也可达10斤/立方水体。

kj=kmj+kajsin(ωjt+φj)

(2)

式中:kmj,kaj分别为各啮合副的平均啮合刚度和刚度变化幅值;φj为各啮合副啮合刚度变化幅值的初相位;ωj为各啮合副的啮合频率。

1.1 综合啮合误差

各齿轮系的综合啮合误差包括齿轮副静态传递误差和各齿轮的偏心误差等，可表示为按啮合频率变化的正弦函数形式[13]

ej=Ejsin(ωjt+φej)

(3)

式中:Ej为各齿轮副综合啮合误差幅值；φej为各齿轮副综合误差初相位。

1.2 齿侧间隙函数

齿轮副间隙非线性函数[14]可表示为

(4)

式中:2b为齿侧间隙;xj为各齿轮副在啮合线上的相对位移。

2 系统动力学方程

根据拉格朗日方程建立系统的纯扭转动力学方程，但由于该方程为半正定、变参数微分方程，在计算中存在刚体位移无法直接求解。为了能够顺利求解，引入相对坐标消除刚体位移xj，公式表示为

xg1g2=rg1θg1+rg2θg2-eg1g2(t),

xg3g4=rg3θg3+rg4θg4-eg3g4(t),

xrpi=rpiθpi-rcθccosα-erpi(t) (i=1,2,…,N)

(5)

(6)

量纲一综合啮合误差为

(7)

式中：Ωj=ωj/ωn，为激励频率。

量纲一齿侧间隙函数为

(8)

在本文计算中，视太阳轮与低速级齿轮g1的转速相同，低速级齿轮g2与高速级齿轮g3的转速相同，即扭转角相等θs=θg1，θg2=θg3。

根据以上各式，可以得到本文研究系统的相对位移坐标下的量纲一化动力学方程表达式

(9)

(10)

(11)

(12)

式中：

i=1,2,…,N

本文视风电机组额定功率为系统的传递功率，其表达式为

(13)

式中：ρair为空气密度;rblade为叶片半径;νwind为平均风速;Cp为风能利用系数。系统的输入扭矩和输出扭矩可表示为

(14)

式中:Λ为风电齿轮传动系统总传动比。

3 系统求解及分析

本文利用4-5阶变步长Runge-Kutta法对系统量纲一振动方程进行求解。为了消除系统的瞬态响应，计算所得数值结果略去前1 500个周期，然后以全局分岔图、最大Lypapunov指数图(以下称LLE)、时间历程图、FFT频谱图、相图及Poincaré截面图作为分析手段，研究说明系统随参数变化的分岔及混沌特性。取风电齿轮传动系统基本参数Pin=1.5 MW，Λ=101，ξ=0.03，α=20°，N=3，ωblade=17 r/min，各齿轮基本参数如表1所示。

3.1 系统随激励频率Ω的动力学特性

取系统各齿轮副的综合啮合误差幅值Ej=150 μm，半齿侧间隙bc=300 μm，计算激励频率Ω在0.001～2之间的系统分岔图和LLE图，分别如图2和图3所示，由于xrpi的动态响应和xspi类似，故以下仅分析xspi。由两幅图可清晰看出系统具有丰富的分岔及混沌特性，且随着激励频率的增加，系统在周期、拟周期运动和混沌运动之间相互切换。激励频率Ω在0.001～0.559变化期间，由分岔图可知系统处于幅值较小的混沌运动，期间也会阵发性地出现单周期运动、2周期运动和拟周期运动。进入混沌道路的途径主要是激变。配合此期间的LLE图可知，LLE在正负间切换即系统在混沌与稳态相互切换；激励频率Ω在0.559～0.771变化期间，系统经倍分岔进入混沌运动，该范围内系统存在单周期运动、3周期运动和12周期运动，期间存在阵发性混沌，在LLE图存在LLE正负间切换的波动现象可知系统在稳态与混沌之间相互切换；激励频率Ω在0.771～0.891变化期间，系统处于混沌运动，且混沌位移增大，LLE为正；当激励频率Ω=0.891时，系统经切分岔由混沌响应进入拟4周期运动，LLE为-0.001 94，随后在Ω=0.923，系统重新进入混沌；激励频率Ω在0.923～1.081期间，系统仍然处于混沌运动；激励频率Ω在1.081～1.161期间，系统从混沌倒分岔至2周期运动，继而“锁相”长周期运动，最终进入混沌运动，由LLE图可知，LLE在0附近波动，然后激变为正值，说明系统在稳态与混沌的临界状态波动最终演变为混沌运动；激励频率Ω在1.161～1.641期间系统运动为混沌-拟2周期运动，配合LLE图可知激励频率变化下，LLE由正值随后逐渐回归负值；当激励频率Ω=1.651时，系统重新进入混沌运动，系统在此时的LLE突变至0.120 7，之后虽然LLE指数存在回落，但是其赋值为正，说明系统处于混沌运动。

表1 系统基本参数

图4～6分别为不同激励频率下系统的FFT频谱图、时间历程图、相图和Poincaré截面图。从图4可以看出，当激励频率Ω=1.101时系统是2周期运动，对应的Poincaré截面为两个点，其频谱图为两个尖峰，即除了其基频外还出现了分频Ω/2和倍频nΩ/2(n为整数)；而从图5可以看出，激励频率Ω=1.141时系统处于拟2周期运动，其Poincaré截面为两个环面，FFT频谱图的基频和分频出现了微小的次谐和超谐成分；在此之后，系统的Poincaré截面图会出现环面破裂形成奇怪吸引子即进入了混沌运动，但随着激励频率的增加，系统从混沌回归至拟2周期运动。最后系统重新进入混沌运动，图6所示为Ω=1.745时，系统的混沌运动，其Poincaré截面出现一个奇怪吸引子且在FFT谱图出现连续频带。

图2 激励频率Ω变化下的分岔图

图3 激励频率Ω变化下的最大Lypapunov指数图

在风电齿轮传动系统中，由于风速的随机性，导致随机输入的存在，同时会引起激励频率的变化，使系统进入混沌。混沌运动是系统不断地由某种轨道突跳到另一个轨道上去的无规则运动，表示了系统不可预测的行为和轨道的永不重复性，使系统受到影响，常会导致疲劳、产生噪声甚至发生故障，进而影响机器寿命。在设计中，应避开混沌状态的范围，采取安全、符合要求的作业范围。通过全局分岔图、时间历程图、FFT频谱图、相图及Poincaré截面图可以定性对系统进行分析，利用全局或单点最大Lyapunov指数图可以定量的对系统进行分析，从而可以更为精确的确定系统的运动状态，保证机器的正常运行。

图7 综合啮合误差变化下的分岔图

图8 综合啮合误差变化下的最大Lypapunov指数图

轮齿之间的设计误差和齿轮的偏心误差等其他的误差综合来说对于系统的作业状态具有很大的影响。在设计时要尽可能提高零件精度，在制造安装方面尽量减少误差，这对提高整个系统的可靠性及稳定性具有明显的作用，同时也可达到降低噪声和减小振动的目的。

图11 激励频率Ω变化下的分岔图

图12 激励频率Ω变化下的分岔图

4 结 论

(1) 本文建立了风电齿轮传动系统纯扭转非线性动力学模型，以风力发电机额定功率1.5 MW为传递功率，采用4-5阶变步长Runge-Kutta法进行求解得到系统非线性动态响应结果。

(2) 利用多种工具对系统响应结果进行分析，在时变啮合刚度、综合啮合误差和齿侧间隙强非线性因素的作用下，系统在激励频率的变化下表现出丰富的分岔特性及混沌现象，如周期运动、拟周期运动和混沌运动，进入混沌运动的途径以倍分岔或激变途径为主。

(3) 综合啮合误差对系统运动分岔特性的影响显著。随着综合啮合误差的增加，系统由周期响应进入混沌响应最后进入拟周期响应；取不同的综合啮合误差进行对比，结果表明综合啮合误差的减小能够减少混沌现象并减弱系统的混沌的区域，使系统更倾向于稳态。因此，在设计中，尽量减小整个系统误差以提高设备的可靠性和寿命。

[1] 卜忠红, 刘更, 吴立言. 行星齿轮传动动力学研究进展[J]. 振动与冲击, 2010, 29(9)：161-166.

BU Zhonghong, LIU Geng, WU Liyan. Research advances in planetary gear trains dynamics[J]. Journal of Vibration and Shock, 2010, 29(9)：161-166.

[2] 邱星辉, 韩勤锴, 褚福磊. 风力机行星齿轮传动系统动力学研究综述[J]. 机械工程学报,2014, 50(11)：23-36.

QIU Xinghui, HAN Qinkai, CHU Fulei. Review on dynamic analysis of wind turbine geared transmission system[J]. Journal of Mechanical Engineering, 2014, 50(11)：23-36.

[3] KAHRARMAN A. Free torsional vibration characteristics of compound planetary gear sets[J]. Mech Mach Theory, 2001, 36：953-971.

[4] 孙智民, 季林红, 沈允文. 2K-H行星齿轮传动非线性动力学[J]. 清华大学学报(自然科学版), 2003, 43(5)：636-639.

SUN Zhimin, JI Linhong, SHEN Yunwen. Nonlinear dynamics of 2K-H planetary gear train[J]. J Tsinghua Univ(Sci&Tech), 2003, 43(5)：636-639.

[5] 张慧博, 王然, 陈子坤,等. 考虑多间隙耦合关系的齿轮系统非线性动力学分析[J].振动与冲击, 2015,34(8):144-150.

ZHANG Huibo, WANG Ran, CHEN Zikun, et al. Nonlinear dynamic analysis of a gear-rotor system with coupled multi-clearance[J]. Journal of Vibration and Shock, 2015, 34(8):144-150.

[6] 孙涛, 沈允文, 孙智民，等. 行星齿轮传动非线性动力学模型方程[J]. 机械工程学报, 2002, 38(3)：6-10.

SUN Tao, SHEN Yunwen, SUN Zhimin， et al. Study on nonlinear dynamic behavior of planetary gear train dynamic model and governing equations[J]. Journal of Mechanical Engineering, 2002, 38(3)：6-10.

[7] 孙涛, 胡海岩. 基于离散傅里叶变换与谐波平衡法的行星齿轮非线性动力学分析[J]. 机械工程学报, 2002, 38(11)：58-61.

SUN Tao, HU Haiyan. Nonlinear dynamics of planetary gear transmission by harmonic balance method based on DFT[J]. Journal of Mechanical Engineering, 2002, 38(11)：58-61.

[8] 李晟, 吴庆鸣, 张志强，等. 两级行星齿轮系分岔与混沌特性研究[J]. 中国机械工程, 2014, 25(7)：931-937.

LI Sheng, WU Qingming, ZHANG Zhiqiang， et al. Bifurcation and chaos characteristics of two-stage planetary gear train sets[J]. China Mechanical Engineering, 2014, 25(7)：931-937.

[9] 李同杰, 朱如鹏, 鲍和云，等. 行星齿轮系扭转非线性振动建模与运动分岔特性研究[J]. 机械工程学报, 2011, 47(21)：76-83.

LI Tongjie, ZHU Rupeng, BAO Heyun， et al. Nonlinear torsional vibration modeling and bifurcation characteristic study of a planetary gear train[J]. Journal of Mechanical Engineering, 2011, 47(21)：76-83.

[10] ZHAO M M, JI J C. Nonlinear torsional vibrations of a wind turbine gearbox[J]. Applied Mathematical Modelling, 2015, 39：4928-4950.

[11] WANG Xiaosun, WU Shijing. Nonlinear dynamic modeling and numerical simulation of the wind turbine’s gear train[C]∥International Conference on Electrical and Control Engineering. Yichang, China: ICECE, 2011.

[12] LIN J, PARKER R G. Planetary gear parametric instability caused by mesh variation[J]. Journal of Soundand Vibration, 2002(1)： 129-145.

[13] PARKER R G. A physical explanation for the effectiveness of planet phasing to suppress planetary gear vibration[J]. Journal of Sound and Vibration, 2000, 236(4)：561-573.

[14] 李润方, 王建军. 齿轮系统动力学——振动、冲击、噪声[M]. 北京：科学出版社, 1997.