基于能量解析法的多轴系齿轮传动系统动力学模型的建立

基于能量解析法的多轴系齿轮传动系统动力学模型的建立袁勇超;王自玲;刘忠明;毛继哲;张坤
【摘要】基于能量解析法给出了转子-轴承-齿轮传动系统的动力学方程.并建立了耦合转子动力学系统的模型.建模过程中考虑了转子轴的剪切、弯曲、扭转、轴向力、陀螺效应及内阻尼的影响,考虑了不同形式的轴承即动压油膜轴承和滚动轴承的影响,同时也考虑了齿轮副的时变啮合刚度和传动误差的影响.为进一步的动力学特性分析提供了精确的模型.
【期刊名称】《机械研究与应用》
【年(卷),期】2019(032)004
【总页数】9页(P33-40,43)
【关键词】能量法;多轴系;齿轮传动;动力学模型
【作者】袁勇超;王自玲;刘忠明;毛继哲;张坤
【作者单位】郑州机械研究所有限公司,河南郑州 450052;河南省体育局,河南郑州 450044;郑州机械研究所有限公司,河南郑州 450052;郑州机械研究所有限公司,河南郑州 450052;郑州机械研究所有限公司,河南郑州 450052
【正文语种】中文
【中图分类】TH132
0 引言
齿轮转子系统的动力学模型一般包括传动轴、轴承支撑及轮盘部分，即常说的转子
-轴承-齿轮动力学系统，这是一个非线性动力系统。

如果某些参数设计不合理，在运转工况下可能会出现异常振动和噪声，这会严重影响传动系统的可靠性及整个机组的运行品质[1]。

所以，对此齿轮耦合转子系统非线性动力学模型的建立和动力学特性研究具有重要的工程实际意义。

充分考虑了转子-轴承-齿轮动力学系统的实际应用性，应用能量数值解析法建立了耦合转子系统的动力学模型，并编写了实际应用程序，为下一步对耦合转子系统非线性动力学特性的研究和分析做出了充分的准备工作。

1 柔性轴段的动力学矩阵
在高速转子系统中，轴可以设计为刚性轴，但由于空间或其他因素如轴承的要求，轴必须设计为柔性轴的几率很大，故笔者主要考虑了柔性轴的动力学分析。

将转子轴划分为若干个轴段，每个轴段单元采用两节点Euler轴单元模型。

下面对考虑了剪切、弯曲、扭转、轴向力、陀螺效应及内阻尼的轴对称Euler轴利用能量法分别求得其刚度矩阵和质量矩阵。

1.1 剪切和转动惯量的影响
在局部坐标系下x-z平面内考虑两节点四自由度细长Euler轴单元的模型如图1所示。

图1 局部坐标系下的杆单元剪切模型
杆的长度为le在杆上任意一点到局部坐标系原点的距离为ζ，如果考虑变截面的杆，需引入杆的形状函数Nei，则有：
(1)
(2)
对轴单元整体分析，其应变能为：
(3)
式中： ue杆上一节点沿x方向的位移；ψe杆上一节点沿y轴方向的角位移；βe 杆上一节点的剪切角；Nei形状函数；Ee杨氏模量；Ie截面惯性矩；Ge剪切模量；ke剪切系数；Ae截面积。

其中：
式中：ν为泊松比；r为轴的半径；Φe无量纲的剪切常量[2]。

对轴单元按自由度分别分析，其应变能又可写成：
(4)
由式(3)和式(4)可分别求得：
kij=EeIeNei″(ξ)Nej″(ξ)dξ+
(5)
从而由式(5)可求得局部坐标系下x-z平面内的轴单元模型的刚度矩阵：
Kse1=
对轴单元整体分析，其动能为：
(6)
对轴单元按自由度分别分析，其动能又可写成：
(7)
由式(6)和式(7)可分别求得：
mij=ρeAeNei(ξ)Nej(ξ)dξ+ρeIe
Nej′(ξ))dξ
(8)
从而由式(8)可求得局部坐标系下x-z平面内的轴单元模型的质量矩阵：
(9)
不难理解，式(9)中第二个矩阵是由转动惯量的影响引起。

其中：
m7=36， m8=(3-15Φe)le
同理可求得局部坐标系下y-z平面内的轴单元模型的刚度矩阵和质量矩阵。

1.2 轴向力的影响
事实上，当轴单元含有内在的拉应力时，其刚度增加；当轴单元含有内在的压应力时，其刚度减小[3]。

在局部坐标系下建立轴单元的轴向力模型如图2所示。

图2 局部坐标系下的杆单元轴向力模型
设其外部拉应力fe恒定不变。

则其应变能为：
(10)
同理按自由度分别分析，写出应变能可分别求得：
kFij=feNei′(ξ)Nej′(ξ)dξ
(11)
KsFe1=
(12)
其中：
为了便于表述计算方法以上只考虑了x-z平面内的自由度情况，下面考虑整体自由度即八个自由度情况，如图3所示。

图3中ve为杆上一节点沿y方向的位移；θe为杆上一节点沿x轴方向的角位移。

同时考虑局部坐标系下x-z和y-z两平面内的自由度，其广义位移向量为[u1 v1
θ1 ψ1 u2 v2 θ2 ψ2]T不难看出其各动力学矩阵是8×8的矩阵。

其推到过程和
1.1节类似，下面直接给出求解结果：
图3 局部坐标系下两平面内的杆单元模型
1.3 内转矩的影响
转子传递扭矩时，轴单元的内转矩对轴单元的横向位移有着显著的影响[4]，在局
部坐标系下建立杆单元的内转矩模型如图4 。

图4 局部坐标系下的杆单元内转矩模型
考虑力矩的影响可求得其应变能为：
(13)
式中：Je为扭转常数；do为轴的外径；di为轴的内径，实心轴时di等于零；
可求得其刚度矩阵为：
KsTe=
式中：te为扭矩，因为扭矩是非保守力，故此矩阵是非对称矩阵。

1.4 陀螺效应的影响
众所周知，旋转着的物体具有像陀螺一样的效应。

它有进动性和定轴性两个特点。

当高速旋转的转子遇到外力时，它的轴向是不随外力矩的方向发生改变的，而是轴
心围绕着一个定点进动。

简单来说，陀螺效应就是旋转的物体有保持其旋转方向的惯性。

尽管对于转子系统中的轮盘来说转子轴的轴颈很小，但在高速旋转的转子系统中，考虑转子轴的陀螺效应还是很有必要的。

对轴单元整体分析，其动能为：
(14)
其中：
=记为[Bij]qe
其中
Aij=-2ρeIeΩB2i(ξ)B1j(ξ)dξ
再由拉格朗日微分方程：
可得：
Gij=-2ρeIe(B2i(ξ)B1j(ξ)-B2j(ξ)B1i(ξ))dξ
(15)
式中：
1.5 内阻尼的影响
轴的内阻尼不同于轴承的阻尼，它主要来源于轴的材料阻尼和内部摩擦。

随着转子的旋转，它会产生相应的阻尼力，而阻尼力的方向往往会与旋转速度方向相反，减
小其广义位移变形量。

用应变能量法可求得其对刚度的影响矩阵结果如下：
综上所解，可得轴单元的各动力学矩阵分别书写如下：
轴单元的刚度矩阵：
Ks=Kse+KsFe+KsTe+Ksce
(16)
轴单元的质量矩阵：
Ms=Mse
(17)
轴单元的陀螺效应矩阵：
Gs=Gse
(18)
2 轴承的动力学矩阵
2.1 各向异性动压油膜轴承
动压油膜轴承工作时转子与轴承有相对运动，它们之间会形成动压油膜，油膜压力会反过来作用于转子上，利用这一工作原理，通过求解雷诺方程可对轴承进行动力学分析。

在局部坐标系下建立动压油膜轴承的工作模型如图5 。

OO1为偏心距，等于偏心率ε与半径间隙c的乘积。

根据流体动压轴承工作原理可求得轴承的径向力、切向力及合力分别为：
式中：D为轴承孔径，等于2R；r为轴半径，等于R-c；Ω为轴的旋转角速度；η为润滑油的绝对粘度；L为轴承宽度；ε为偏心率；c为半径间隙。

图5 动压油膜轴承几何尺寸示意图
其中：偏心率可由下面的一元八次方程求解：
式中：Ss为修正的索莫菲因
其刚度与阻尼为：
(19)
2.2 各向同性滚动轴承
轴承的刚度除了取决于轴承的类型和尺寸外，最重要的影响参数有滚动体的类型(球或滚子)；滚动体的数量和尺寸；接触角。

通过预加负荷能够增强轴承的刚度。

试验齿轮箱中用到的均是球轴承，故这里只讨论球轴承的动力学特性。

其结构尺寸示意图如图6所示。

图6 滚动轴承结构尺寸示意图
轴承样本手册[6]中给出了滚子与滚道间的间隙计算公式如式(20)所示，不同厂家不同精度的轴承其a的取值不同。

当滚子直径与滚子节圆直径相差较大时，可以忽略其间隙的影响，将滚动轴承转化为各向同性轴承来分析，这样其动力学矩阵可求得如式(21)所示：
δ=a×10-8d-1/3f2/3
(20)
(21)
其中：
kuu=kvv=krn2/3d1/3f1/3cos5/3α
kuv=kvu=0
kr=13×106(N2/3/m4/3)
Cbr=10-5Kbr
式中：d为滚子直径；f为轴承所受径向力；kr为刚度常数；n为滚子数量；α为轴承接触角。

3 轮盘的动力学矩阵
在转子系统中，轴上可以安装不同形式的轮盘如叶轮、齿轮、阻尼盘及联轴器等，这里先把它们理想成为短圆柱形刚性轮盘来作整体分析。

下面讨论轮盘的质量矩阵及陀螺效应矩阵。

在总体坐标系下建立轮盘的工作模型如图7 。

首先根据传递矩阵法，考虑均质刚性轮盘在总体坐标系下的动能。

在总体坐标系下轮盘沿各坐标轴的瞬时角速度如式(22)所示：
(22)
式中：ωx为绕x轴的瞬时角速度；ωy为绕y轴的瞬时角速度；ωz为绕z轴的瞬时角速度；Ω为局部坐标系下绕z轴的瞬时角速度；φ为局部坐标系下绕z轴的转动角度；θ为局部坐标系下绕x轴的转动角度；ψ为局部坐标系下绕y轴的转动角度。

对轮盘整体分析其总的动能可表示为式(23)：
(23)
考虑θ与ψ是微小量，忽略高阶无穷小后得：
(24)
式中：Td为总动能；md为轮盘质量；Id为直径转动惯量；Ip为转动惯量；u为
沿x轴的位移；v为沿y轴的位移。

不难理解上式最后一项是陀螺效应的影响。

式中广义位移向量为[u v θ ψ]T对上式利用拉格朗日微分方程可得:
其中：Md为轮盘质量矩阵；Gd为轮盘陀螺效应矩阵：
图7 总体坐标系下轮盘的广义位移
4 各动力学矩阵的叠加
对于轴对称转子，相对于总体坐标系所有自由度位置要一一对应，没有的用零补齐。

第一个轴单元的自由度是1～8，第二个轴单元的自由度是5～12，第三个轴单元
的自由度是9～16，依次类推。

轴的整体矩阵是所有单元相应矩阵的叠加。

在第
三个节点上的轮盘自由度是9～12，同理其他节点上的轮盘自由度也易得出，并
可叠加到轴模型的矩阵上去，如果轮盘是齿轮应考虑耦合叠加矩阵。

轴承自由度叠加时应考虑旋转自由度的放开位置[7]。

具体情况见示意图8。

图8 矩阵叠加原理示意图
5 齿轮副动力学模型
5.1 单对齿轮副动力学模型
图9为一对相互啮合的斜齿轮副动力学模型，端面压力角为α12，螺旋角为β12；定义齿轮1为主动齿轮，齿轮2为从动齿轮，齿轮1右旋，β12>0，反之亦然。

基圆半径分别为r1,r2；主动力矩为T逆时针为正；为了实现局部坐标系向总体坐
标系的直接转化，引入角度δ12和σ12，在总体坐标系易得：
该齿轮副共有12个自由度。

即每个齿轮都具有三个平动自由度及三个旋转自由度，
按XYZ右旋坐标系依次可记为：
q12=[u1 v1 w1 θ1 ψ1 φ1 u2 v2 w2
θ2 ψ2 φ2]T
图9 斜齿轮副动力学模型
齿轮副沿啮合线上的位移变形量，主要由齿轮副扭转引起的变形，横向振动引起的变形，轴向振动引起的变形，绕直径摆动引起的变形及误差变形五部分组成，可表示为(25)式所示：
g12=(-w1+w2)sin β12+(r1φ1-r2φ2)cos β12+
((u1-u2)sin δ12+(v1-v2)cos δ12)
cos β12+((r1θ1+r2θ2)sin δ12+(r1φ1+
r2φ2)cos δ12)·cos β12-e(t)12
(25)
式中：e(t)12为齿轮副的传动误差。

轮齿副的动态啮合力及其沿各个坐标方向的分力分别为：
(26)
式中：k(t)g1是时变啮合刚度，它是转子转频和时间的函数。

由图9可求得齿轮副的动力学方程如式(27)。

(27)
集成轴、轴承和齿轮，系统的综合动力学方程如式(28)：
(28)
式中：M12= MS+Md
C12= Cb+GS+ Gd
K12= Kb+KS+ Kd
将静态传动误差e(t)12和时变啮合刚度k(t)g1转化为傅里叶级数展开式[8]，代入上面方程中并采用四阶变步长龙格-库塔数值积分法[9]求解运动微分方程，可得各分析结果。

5.2 多轴系齿轮副动力学模型
结合上节的分析，将图9中的角码1和2改成i和j可建立多对齿轮副的动力学模型和动力学方程。

笔者试验齿轮箱中包含8根轴共有7对啮合副，其中串联啮合
副3对。

共有48个自由度如图10。

图10 多轴系齿轮副模型
其动力学方程为：
(29)
根据图10中的齿轮序号可知式中：ij={1，3；2，4；4，6；6，7；5，8；8，9；9，10}。

然后采用四阶变步长龙格-库塔数值积分法求解运动微分方程组。

6 结论
充分考虑了转子-轴承-齿轮动力学系统的实际应用性，应用能量数值解析法建立了耦合转子系统的动力学模型，并编写了实际应用程序，为进一步对耦合转子系统非线性动力学特性的研究和分析做出了充分的准备工作。

参考文献：
【相关文献】
[1] 陆春荣.齿轮-转子-轴承系统弯扭耦合非线性振动特性研究[J].振动工程学报,2018,31(2):67-72.
[2] Inman,D.J. Engineering Vibration[M].Third Edition. New Jersey: Prentice-Hall,2007:139-145.
[3] Dawe,D.J.Matrix and Finite Element Displacement Analysis of
Structures[M].Oxford:Oxford University Press,1984:134.
[4] Petyt, M. Introduction to Finite Element Vibration Analysis[M].Cambridge:Cambridge University Press,1990:148.
[5] Childs,D. Turbo machinery Rotor dynamics Phenomena Modeling and Analysis [M].New York: Wiley,1993:104.
[6] International Organization for Standardization. Rolling bearings-Dynamic load ratings and rating life:ISO281:2007(E)[S].Second edition. Switzerland: International Organization for Standardization,2007:23-25.
[7] Frederick F. Ling. Dynamics of Rotating Systems[M].Germany: Springer Group,2005:219.
[8] Chen Siyu,Tang Jinyuan.Dynamics analysis of a crowned gear transmission system with impact damping: Based on experimental transmission error[J].Mechanism and Machine Theory 2014(74):242-259.
[9] He Lin,Sanmin Wang,Earl H.Dowell,Jincheng Dong,Congma.Dynamic Characteristics of Double-Helical Planetary Gear Sets Under Time-Varying Mesh Stiffness[J].Journal of Harbin Institute of Technology,2016,23(4):44-50.。