行星齿轮传动具有体积小、重量轻、速比大、效率高等特点，在航空、船舶传动、汽车、起重机械及其他机械传动中获得了越来越广泛的应用，吸引了众多专家学者的研究目光^{[1, 2, 3, 4, 5]}。齿轮副间存在的时变啮合刚度、齿侧间隙等因素，决定了行星齿轮传动系统强非线性的本质^[6]，文献[7, 8, 9, 10, 11, 12]的研究结论也证实了行星齿轮传动系统存在诸如运动状态的稳定性分岔、幅值跳跃等非线性动力学特征。

在工程实际中，机器的启动、停车以及工作状态的转换都会使得行星齿轮系统所传递的功率频繁变化，因而非常有必要研究行星齿轮传动系统随传递功率的分岔特性，以揭示传递功率的变化对其运动状态稳定性的影响规律，获得各种稳定的周期轨道与相应传递功率区间的映射关系。目前，有关行星齿轮传动系统分岔特性方面的研究还比较少见。文献[6]采用直接数值积分法，定性研究了行星齿轮传动系统随齿侧间隙等参数的全局分岔规律。文献[8]采用改进的PNF方法结合Floquet原理研究了行星齿轮传动系统的共存周期轨道的稳定性随转速的局部精细分岔规律。而关于行星齿轮传动系统随传递功率的分岔特性研究尚未见报道。

考虑到PNF(Poincaré-Newton-Floquet)方法需要在每次迭代计算前单独预估迭代初始值，计算效率较低，笔者将基于虞烈等^[13]提出的CPNF(continuous Poincaré-Newton-Floquet)方法研究传递功率对行星齿轮传动系统周期轨道稳定性的局部精细分岔规律，并采用直接数值积分的方法绘制系统随功率的全局分岔图，以验证仿真结果的可靠性。

本文力学模型采用文献[6]中的2K-H型行星齿轮传动系纯扭转非线性动力学模型，以下是该模型的简单介绍。2K-H型行星齿轮系统有1个太阳轮(以S表示)，N个行星齿轮(第i个齿轮以p_i表示)，1个行星架(以C表示)和1个内齿圈(以R表示，固定于基础之上)组成，各齿轮为直齿圆柱齿轮，系统的扭转振动模型如图 1所示。

图1(Fig. 1)

图1 系统扭转非线性振动模型 Fig. 1 Nonlinear vibration model of the system

图中，太阳轮、行星架、第i个行星轮的角位移分别以θ_s、θ_c、θ_pi表示；太阳轮、第i个行星轮、内齿圈的基圆半径分别以r_bs、r_bpi、r_br表示；行星架半径以r_c表示，其值为太阳轮与行星轮的节圆半径之和，标准安装下亦即太阳轮与行星轮的分度圆半径之和；太阳轮与第i路行星轮组成的外啮合副的啮合刚度、啮合阻尼系数、半齿侧间隙、综合啮合误差分别以k_spi、c_spi、b_spi、e_spi表示；内齿圈与第i路行星轮组成的内啮合副的啮合刚度、啮合阻尼系数、半齿侧间隙、静传递误差分别以k_rpi、c_rpi、b_rpi、e_rpi表示，各齿轮的压力角均为α。

引入如下量纲一的参数$\Omega =\omega /{{\omega }_{n}}$, $\tau ={{\omega }_{n}}t$, ${{\omega }_{n}}=\sqrt{\frac{{{k}_{\text{msp}i}}}{{{m}_{\text{eq}1}}}}$, ${{m}_{\text{eq1}}}=\frac{{{m}_{\text{c}}}{{m}_{\text{s}}}}{{{m}_{\text{c}}}+{{m}_{\text{s}}}},$${{m}_{\text{eq2}}}=\frac{{{m}_{\text{c}}}{{m}_{\text{s}}}}{{{m}_{\text{c}}}+2{{m}_{\text{s}}}}$, ${{m}_{\text{s}}}=\frac{{{J}_{\text{s}}}}{r_{\text{bs}}^{2}}$, ${{m}_{\text{p}i}}=\frac{{{J}_{\text{p}i}}}{r_{\text{bp}i}^{2}}$, ${{m}_{\text{c}}}=\frac{{{J}_{\text{c}}}+\sum\limits_{i=1}^{3}{{{m}_{\text{p}i}}r_{\text{c}}^{2}}}{r_{\text{bc}}^{2}}$, $X=\bar{X}{{b}_{\text{c}}}$, $\dot{X}=\dot{\bar{X}}{{b}_{\text{c}}}{{\omega }_{n}}$, $\ddot{X}=\ddot{\bar{X}}{{b}_{\text{c}}}\omega _{n}^{2}$, $b={{b}_{\text{c}}}\bar{b}$, $e={{b}_{\text{c}}}\bar{e}$。其中，J_s、J_pi、J_c分别为太阳轮转动惯量、第i个行星轮转动惯量、行星架转动惯量；b_c为位移标称尺度，X表示相应齿轮副的啮合点相对位移；e表示相应齿轮副的静传递误差；ω表示齿轮的啮合频率，ω_n表示固有频率；k_mspi表示系统外啮合副的平均啮合刚度，字母上加“-”表示相应物理量量纲一的结果。

根据刚体平面运动微分方程并结合上述量纲一的参数，容易列写行星轮系的纯扭转量纲一的振动微分方程为

方程(1)中，f表示相应齿轮副上的间隙非线性函数；k_mrpi表示系统内啮合副的平均啮合刚度；k_aspi表示系统外啮合副刚度波动幅值；k_arpi表示系统内啮合副刚度波动幅值；φ_spi表示系统第i路外啮合副刚度波动初相位；φ_rpi表示系统第i路内啮合副刚度波动初相位；F_D为作用在太阳轮上的等效扭矩，F_L为作用在行星架上的等效扭矩。

CPNF^{[13, 14]}法是将确定参数下用来判断系统周期轨道稳定性的PNF方法与参数延续算法相结合形成的一种计算方法。这种方法利用预测校正原理，根据上一个参数值下求得的周期轨道对下一个参数值下的迭代初始值进行预估，很好地解决了迭代的初始值问题。该方法的另一优点是在迭代计算周期轨道的过程当中可以同步获得周期轨道的状态转移矩阵，使周期轨道的求解与稳定性判断工作可以同步进行，因而大大提高了周期轨道稳定性局部分岔的计算效率。

对于如下含有参数ω的n维非线性动力系统

应用庞加莱映射，可以把求解该非线性动力系统周期解曲线的问题转变为求解点映射系统

在庞加莱截面上的不动点

的问题，即求解非线性代数方程

的解曲线问题。CPNF方法将该问题的求解过程分成了预测和校正两步。

假设已经求得ω=ω^k时周期解在庞加莱截面上的不动点u^k，那么采用欧拉型积分公式，可以得到如下预测公式

式中:P′_u(ω_k,u_k)=DP(ω_k,u_k)为当参数取ω=ω_k时Poincaré映射(3)在u_k处的Jacobi矩阵;P′_ω(ω_k,u_k)为Poincaré映射(3)在参数ω=ω^k处的梯度向量;P′_u(ω_k,u_k)及P′_ω(ω_k,u_k)可以在第k步校正过程中获得，u_(k+1)⁰为第(k+1)步校正计算的迭代初始值。

第k步校正过程实际就是利用PNF方法求解参数ω=ω_k时系统周期解的牛顿迭代过程，不同的是，牛顿迭代的初始值已经由前一步的预测过程给出，而不需要由计算者人为给定，从而大大提高了求解的成功率，节省了计算时间。其系统第k步周期解的迭代公式为

其中，DP(ω_k ,u_k⁽ⁱ⁾)以及P(ω_k,u_k⁽ⁱ⁾)由联立微分方程组

以(E,u_kⁱ)为初始值积分一个庞加莱映射周期T得到。式(8)中Df_u(u,ω,t)为非线性动力系统(2)的Jacobi矩阵。当利用公式(7)迭代求出满足精度要求的第k步周期解u_k时，就同时得到了相应的状态传递矩阵DP(ω_k,u_k)，可以用于周期解u_k的稳定性分析以及第(k+1)步的预测。

为了实现第(k+1)步的预测，还需要计算梯度向量P′_ω(ω_k,u_k)，可以通过如下微分方程组

由以(0,u_k)为初始值积分一个庞加莱映射周期T，即可得到(P′_ω(ω_k,u_k),P(u_k))。式(9)中，Df_ω(u,ω,t)为非线性动力系统(2)的关于参数ω的梯度向量。

以下运用CPNF方法研究行星齿轮非线性系统的传递功率P在196~220 kW变化时，系统周期运动稳定性的分岔特性。行星齿轮系统的基本参数如下：模数m=1.75 mm，压力角α=20°，太阳轮齿数z_s=24，行星轮齿数z_pi=39，内齿圈齿数z_r=106，标称尺度b_c=10 μm，行星轮个数N=3，内外啮合副齿侧间隙b_spi=b_rpi=50 μm，量纲一的转速Ω=0.958 6。

CPNF方法在进行校正时，涉及到计算非线性动力系统(1)对应状态方程的Jacobi矩阵Df_u(u,ω,t)。考虑到纯扭转非线性动力学模型的非光滑性，Jacobi矩阵Df_u(u,ω,t)采用有限差分形式近似代替^[14]。

当传递功率P=220 kW时，经PNF方法计算，只获得了一个稳定的周期一不动点P₁(0.993 218 5,-2.647 681 1,5.401 011 1,-1.933 860 6)，其轨道相图如图 2所示。

图2(Fig. 2)

图2 P=220 kW时稳定的周期一轨道 Fig. 2 Stable phase thajectory of P1 when P=220 kW

由此考虑在小于P=220 kW的某段功率区间内，系统的运动可能以周期一轨道的形态稳定运行，为了确定周期一轨道保持稳定的功率范围以及具体的失稳分岔点功率值，以下通过CPNF方法计算了周期一轨道的Floquet乘子的最大模随功率P的变化曲线，如图 3所示。

图3(Fig. 3)

图3 周期一轨道最大Floquet乘子模随功率的变化曲线 Fig. 3 The changing curve of P1 Floquet multiplier with power

从图 3中可以看到，在P为211~220 kW范围内，周期一轨道的Floquet乘子的最大模小于1，即运动状态能够保持稳定；在P为196~211 kW范围内，Floquet乘子的最大模大于1，即周期一运动依然存在，但已经发生了Lyapunov意义下的失稳^[15]，变成不稳定的运动形态。由此可确定，周期一轨道稳定性的分岔功率为P=211 kW。

为了确定周期一轨道的稳定性在P=211 kW处具体的分岔途径，笔者模拟了在功率P=211 kW前后，周期一轨道具有最大模的Floquet乘子穿出复数域上单位圆的途径，如图 4所示，图中箭头方向指示了当传递功率由大变小时，周期一轨道Floquet乘子的变化走向。显然，Floquet乘子是由实轴负半轴穿出单位圆的，根据周期轨道分岔理论^[15]，当传递功率由大变小经过分岔点功率P=211 kW时，周期一轨道会发生倍周期分岔而失稳裂变成周期二轨道。

图4(Fig. 4)

图4 周期一轨道Floquet乘子随功率穿出单位圆的途径 Fig. 4 The way across the unit circle of P1 Floquet multiplier

由以上周期一运动的稳定性分岔分析可知，在小于211 kW的某个功率区间内，系统的运动将会以稳定的周期二轨道的形态存在。

当传递功率P=211 kW时，经PNF方法计算，获得了一个不稳定的周期一不动点P₁(1.122 263 389 05,-2.681 586 539 89,5.462 945 802 425,-1.903 188 660 61)和一个稳定的周期二不动点P₂(0.956 373 419 33,-2.724 380 076 71,4.664 356 134 64,-1.953 595 730 27)，其周期轨道相图分别如图 5和图 6所示。

图5(Fig. 5)

图5 P=211 kW时不稳定的周期一轨道 Fig. 5 Unstable phase thajectory of P1 when P=211 kW

图6(Fig. 6)

图6 P=211 kW时稳定的周期二轨道 Fig. 6 Stable phase thajectory of P2 when P=211 kW

以上计算结果验证了CPNF方法关于周期一轨道在P=211 kW处发生倍周期分岔的预测。为了确定周期二轨道保持稳定状态的功率范围以及具体的失稳分岔点功率值，通过CPNF方法计算的周期二运动的Floquet乘子的最大模随功率P的变化曲线如图 7所示。

图7(Fig. 7)

图7 周期二轨道最大Floquet乘子模随功率的变化曲线 Fig. 7 The changing curve of P2 Floquet multiplier with power

从图 7中可以看到，在P为205.8~211 kW范围内，周期二轨道的最大Floquet乘子模小于1，即运动状态能够保持稳定；在P为196~205.8 kW范围内，最大Floquet乘子模大于1，即周期二运动依然存在，但已经发生了Lyapunov意义下的失稳^[9]，变成了不稳定的运动形态。由此可以确定，周期二轨道稳定性的分岔功率为P=205.8 kW。

为了确定周期二轨道稳定性在P=205.8 kW时具体的分岔途径，笔者模拟了功率由大变小经过P=205.8 kW时，周期二轨道具有最大模的Floquet乘子穿出复数域上单位圆的途径，如图 8所示。显然，Floquet乘子是由实轴负半轴穿出单位圆的，根据周期轨道分岔理论^[15]，当传递功率由大变小经过分岔点功率P=205.8 kW时，周期二轨道会发生倍周期分岔而失稳裂变成周期四轨道。

图8(Fig. 8)

图8 周期二轨道Floquet乘子随功率穿出单位圆的途径 Fig. 8 The way across the unit circle of P2 Floquet multiplier

由以上周期二运动的稳定性分岔分析可知，在小于205.8 kW的某个功率范围内，系统的运动将会以稳定的周期四轨道的形态存在。

当传递功率为205.8 kW时，经PNF方法计算，获得了一个不稳定的周期一不动点P₁ (1.246 607 377 329,-2.707 429 571 563,5.523 882 489 132,-1.867 422 471 809)、一个不稳定的周期二不动点P₂ (2.098 566 305 065,-2.734 138 541 796,8.311 191 163 223,-1.323 076 184 707)和一个稳定的周期四不动点P₄(2.797 053 954 0,-2.606 852 968 25,8.572 761 368 432 6,-1.587 082 095 98)，其周期轨道相图分别如图 9~11所示。

图9(Fig. 9)

图9 P=205.8 kW时不稳定的周期一轨道相图 Fig. 9 Unstable phase thajectory of P1 when P=205.8 kW

图10(Fig. 10)

图10 P=205.8kW时稳定的周期二轨道相图 Fig. 10 Unstable phase thajectory of P2 when P=205.8 kW

图11(Fig. 11)

图11 P=205.8 kW时不稳定的周期四轨道相图 Fig. 11 Stable phase thajectory of P4 when P=205.8 kW

以上计算结果验证了CPNF方法关于周期二轨道在P=205.8 kW处发生倍周期分岔的预测。为了确定周期四轨道保持稳定的功率范围以及具体的失稳分岔点功率，以下通过CPNF方法计算了周期四轨道的Floquet乘子的最大模随功率P的变化曲线，如图 12所示。

图12(Fig. 12)

图12 周期四轨道最大Floquet乘子模随功率的变化曲线 Fig. 12 The changing curve of P4 Floquet multiplier with power

从图 12中可以看到，在196.8~205.8 kW范围内，周期四轨道Floquet乘子的最大模小于1，即运动状态能够保持稳定；在196.0~196.8 kW范围内，Floquet乘子的最大模大于1，即周期四轨道依然存在，但已经发生了Lyapunov意义下的失稳^[15]，变成不稳定的运动形态。由此可确定，周期四轨道稳定性的分岔功率为P=196.8 kW。

为了确定周期四轨道稳定性在P=196.8 kW时具体的分岔途径，笔者模拟了功率由大变小经过P=196.8 kW时，周期四轨道具有最大模的Floquet乘子穿出复数域上单位圆的途径，如图 13所示。显而易见，Floquet乘子是由实轴负半轴穿出单位圆的，也就是说随着功率的由大变小，周期四轨道在P=196.8 kW处会发生倍周期分岔而失稳裂变成周期八轨道^[15]。

图13(Fig. 13)

图13 周期四轨道Floquet乘子随功率穿出单位圆的途径 Fig. 13 The way across the unit circle of P4 Floquet multiplier

当传递功率为196.0 kW时，采用PNF方法计算获得了如下4个共存的周期不动点：

P₁(1.326 969 56,-2.720 465 9,5.563 583 4,-1.840 685 9)；

P₂(2.478 602 02,-2.779 455 9,8.898 100 5,-0.994 215 5)；

P₄(-0.007 560 10,-2.361 369 9,2.470 765 6,-1.816 708 5)；

P₈(3.313 389 72,-2.629 371 0,8.819 356 8,-1.698 036 1)。

同时得到了以上各周期不动点的状态转移矩阵P₁、P₂、P₄、P₈。

进一步分析以上4矩阵的特征值模的最大值，亦即Floquet乘子模的最大值，分别为λ(P₁)=1.278 438 866,λ(P₂)=1.605 355 509,λ(P₄)=1.320 196 013 62,λ(P₈)=0.691 605 07，显然在该组参数下只有周期八运动稳定。

总结以上分析，可以大致勾勒出传递功率在196.0~220.0 kW范围内系统周期运动状态稳定性的分岔情况，如图 14所示。图中实线表示稳定的周期运动，点划线表示不稳定周期运动。作为对比，笔者采用直接数值积分的方法绘制了行星齿轮传动系统随功率的全局分岔图，如图 15所示。可以发现在功率区间196.0~220.0 kW内，两图所显示的分岔规律以及分岔点功率值非常吻合。需要指出的是，基于CPNF方法预测的局部分岔图可以获得所有共存周期轨道(包括不稳定的周期轨道)的分岔规律，而基于直接数值积分的全局分岔图只能显示渐进稳定的周期轨道之间的分岔关系，这是CPNF方法预测的局部分岔图的优点之一。此外，由于局部分岔图中14所显示的分岔点参数是由Floquet乘子定量判断而得，所以比全局分岔图 15定性显示的分岔点参数值更加精确。

图14(Fig. 14)

图14 周期运动稳定性随传递功率的局部分岔图 Fig. 14 The local bifurcation diagram with power

图15(Fig. 15)

图15 行星齿轮系统随传递功率P的全局分岔图 Fig. 15 The global bifurcation diagram with power

从分岔图 14中可以看到，在功率196.0~196.8 kW范围内，行星齿轮非线性系统共存一个不稳定的周期一运动、一个不稳定的周期二运动、一个不稳定的周期四运动和一个稳定的周期八运动；在功率196.8~205.8 kW范围内，行星齿轮非线性系统共存一个不稳定的周期一运动、一个不稳定的周期二运动和一个稳定的周期四运动；在功率205.8~211.0 kW范围内，行星齿轮非线性系统共存一个不稳定的周期一运动和一个稳定的周期二运动；在功率211.0~220.0 kW范围内，行星齿轮非线性系统呈现单一而稳定的周期一运动。

总而言之，在传递功率196.0~220.0 kW范围内，随着功率值的逐渐增大，行星齿轮传动系统的各周期轨道均是通过倍周期倒分岔的途径在相应功率分岔点处发生稳定性突变的。全局分岔图 15显示，经过充分的倍周期分岔后，在小于196.0 kW的某个功率值上，行星齿轮非线性系统会呈现混沌运动状态。

1) 采用CPNF方法研究获得的行星齿轮传动系统的周期运动轨道稳定性随传递功率的局部分岔规律，与通过直接数值积分法得到的的全局分岔规律基本吻合，证明了CPNF方法在研究行星齿轮传动系统分岔特性方面的有效性。

2) 在某些确定的参数组合下，行星齿轮传动系统会共存多个周期轨道，这些轨道有的稳定，有的不稳定；随着功率值的逐渐增大，行星齿轮传动系统的渐进稳定的周期轨道会通过倍周期倒分岔的途径发生失稳；在轻载工况下，行星齿轮非线性系统容易呈现混沌运动状态。