预变形Jeffcott裂纹转子超谐波共振的动力学响应

Dynamics Response of Super-harmonic Resonance of a Pre-deformed Jeffcott Cracked Rotor

摘要

以转子动力学和非线性动力学理论为基础，针对刚性支承的水平Jeffcott 裂纹转子，研究1∶1内共振条件下系统的2/3超谐波共振的动力学响应。考虑了重力引起的预变形效应和裂纹呼吸行为，建立了由转子预变形引入的平方非线性和立方非线性的系统动力学方程。采用多尺度法推导出系统的演化方程，详细讨论了系统参数变化对系统动力学行为的影响规律，表明随着裂纹深度和横向阻尼系数的增大，系统频响曲线会出现“频率岛”现象；阐明了平方非线性对系统有软特性，立方非线性对系统有硬特性。对系统动力学方程采用Runge-Kutta法进行了数值积分来观察其跳跃现象，并对多尺度的近似解进行了验证。研究结果为Jeffcott转子系统裂纹故障的非线性动力学分析提供了理论指导。

Abstract

Based on the theories of rotor dynamics and nonlinear dynamics， the dynamics response of 2/3 super-harmonic resonance of the system was studied for the rigidly supported horizontal Jeffcott cracked rotors under the condition of 1∶1 internal resonance. Considering the gravity-induced pre-deformation effect and crack respiration behavior， a system dynamical equation with quadratic nonlinearity introduced by pre-deformation and cubic nonlinearity was established. The multi-scale method was used to derive the evolution equation of the system， and the influences of the changes of system parameters on the dynamics behaviors of the system were discussed in detail， and it is shown that with the increasing of crack depths and lateral damping coefficients， the frequency response curves of the system exhibit “frequency island” phenomenon. It is clarified that the square nonlinearity exhibits a softening effect on the system， while the cubic nonlinearity exhibits a hardening effect. Numerical integration of the system's dynamics equation was conducted using the Runge-Kutta method to observe the jumping phenomenon， and the approximate solutions obtained through the multi-scale method were validated. The findings provide theoretical guidance for the nonlinear dynamics analysis of crack faults in Jeffcott rotor systems.

Graphical abstract

Author summay

张博^*（通信作者），男，1989年生，副教授。研究方向为非线性振动与控制、多场环境下复合材料壁板动力学行为、转子动力学。E-mail： [e-mail]。

0 引言

在我国的电力和化学工业中，涡轮和压缩机是最重要和最普遍使用的机械设备。转子作为这类旋转机械的核心部件，常常由于工况复杂，容易引起疲劳裂纹^［1］。为了避免发生灾难性事故而造成经济损失和意外风险，必须及时检测和识别疲劳裂纹，因此，建立准确的转子动力学模型，并深入研究其不同共振形式下的动力学响应机理具有重要意义。

Jeffcott转子最早由Foppl于1895年提出^［2］，其结构组成为一根两端简支无质量的轴与安装在轴中央的圆盘。裂纹转子^［3］的非线性动力学行为的复杂性主要源于裂纹引起的时变刚度。裂纹模型目前已有的研究主要分为常开型和呼吸型。对于裂纹开闭规律的演绎模型，目前常见的主要有方波模型^［4］、余弦模型^［5］和高朱模型^［6］。SINOU等^［7］研究的呼吸裂纹模型是通过改进得到的一种三角函数模型，更细致地描述了裂纹开合的过程。相比于简单的方波模型，三角函数模型能得到呼吸裂纹更多的动力学行为。

在裂纹转子动力学响应的研究方面，国内外学者通过仿真与实验等方法做了大量的研究。DIMAROGONAS等^［8］对具有横向裂纹转子的动力学行为进行了全面的总结，并对其动力学特性进行了综述。路振勇等^［9］对具有裂纹的航空发动机转子进行了非线性振动分析。胡亮等^［10］对带有横向裂纹的拉杆转子进行了非线性动力分析。ZHANG等^［11］研究了由横向呼吸裂纹转子引发的一系列谐波共振的六自由度耦合动力学响应。基于呼吸裂纹模型，JUN等^［12］建立了一个简单转子系统的运动方程，研究发现二次水平谐波分量可以识别出裂纹引起的振动行为。XIE等^［13］在研究工作中发现了裂纹转子中的超谐波分量。李常有等^［14］的研究表明超谐波响应成分可作为转子裂纹故障的诊断依据。CHEN等^［15］讨论了旋转时裂纹方向和不平衡偏心率之间的取向角对Jeffcott裂纹转子水平振动的影响。HOU等^［16］通过多尺度法研究了系统在1∶2和1∶3的超谐波共振下裂纹转子的非线性振动，发现了由呼吸裂纹和惯性激励之间的相互作用产生的谐波共振。

Jeffcott转子模型由于具有对称性使得转子竖直方向和水平方向的固有振动频率接近，对于考虑非线性轴承刚度的转子模型，1∶1内共振可公度条件满足。内共振条件的存在会使得转子呈现出丰富且具有研究价值的复杂动力学现象。王宗勇等^［17］研究了质量、激励频率、阻尼以及刚度突变对转子动力学特性的影响。SINOU等^［18］分析了具有呼吸裂纹转子系统的动态响应，发现在内共振条件下激发的次谐波分量可以用来检测系统中的裂纹。YABUNO等^［19］研究了水平支撑式Jeffcott转子在内共振条件下一阶主共振的动力学行为和局部分岔情况。SAEED等^［20］讨论了具有横向裂纹的垂直支撑式Jeffcott转子系统的非线性动力学行为和水平支撑式Jeffcott转子的一阶主共振的动力学行为。

在过往的研究中，也有学者报道了自重引起的预变形对裂纹转子动力学响应的影响。何成兵等^［21］研究了以重力占优的呼吸裂纹转子，对裂纹转子的升速瞬态响应进行了仿真，结果表明自重与质量偏心等因素有着必然的联系。杨丹等^［22］证明了自重会直接或间接地影响Jeffcott转子的动力学行为及稳定性。

通过调研不难发现，裂纹转子领域中不考虑内共振或不依赖重力占优的研究较多，但有关同时考虑预变形及内共振条件下谐波共振的动力学行为的研究报道较少，因此本文研究了含横向裂纹的预变形Jeffcott转子的2/3超谐波共振的动力学行为，揭示了自重引起的预变形效应对裂纹转子超谐波共振的影响规律。

1 模型描述和系统动力学方程

本文选取理想的Jeffcott模型由一根两端刚支的无质量轴与位于轴中部的圆盘组成，转轴长度为l，圆盘质量为m。在转轴正中间假设具有一条横向裂纹，如图1a所示。以水平轴线和圆盘中心交点为原点O，建立惯性坐标系OXY，其中X、Y轴分别表示水平方向（横向）和竖直方向（竖向）。如图1b和图1c所示，建立OX₁Y₁为圆盘的静止坐标系，OX₂Y₂为旋转坐标系。几何中心为G （x₁， y₁），转子角速度为Ω，轴承的弹性恢复力F是其弯曲时圆盘几何中心挠度r的三次非线性函数，可表示为^［23］

F(r)=kr+\beta r^{\mathrm{3}}（1）

式中：k为线性刚度系数；β为非线性刚度系数。

在静平衡状态下，由重力引起轴的静挠度y_st需满足

k{y}_{\mathrm{s}\mathrm{t}}+\beta y_{\mathrm{s}\mathrm{t}}^{\mathrm{3}}=-mg（2）

式中：g为重力加速度。

当圆盘旋转时，几何中心G在X₁和Y₁方向的运动因静平衡状态下产生的位移而发生变化，分别用x₁、y₁表示，则有x=0+x₁，y=y_st+y₁，结合式（1）、式（2），仅考虑转轴上的横向裂纹，忽略裂纹对系统质量和阻尼的影响，利用Lagrange方程建立如下系统的运动方程^［23］：

\left.\begin{array}{l}m{\ddot{x}}_{\mathrm{1}}+c_{\mathrm{1}}{\dot{x}}_{\mathrm{1}}+(k+\beta y_{\mathrm{s}\mathrm{t}}^{\mathrm{2}})x_{\mathrm{1}}+\mathrm{2}\beta y_{\mathrm{s}\mathrm{t}}{x}_{\mathrm{1}}{y}_{\mathrm{1}}+\\ \beta x_{\mathrm{1}}{y}_{\mathrm{1}}^{\mathrm{2}}+\beta x_{\mathrm{1}}^{\mathrm{3}}-F_{\mathrm{1}}=m{e}_{\mathrm{d}}{\Omega }^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}(\Omega t+\gamma )\\ m{\ddot{y}}_{\mathrm{1}}+c_{\mathrm{2}}{\dot{y}}_{\mathrm{1}}+(k+\mathrm{3}\beta y_{\mathrm{s}\mathrm{t}}^{\mathrm{2}})y_{\mathrm{1}}+\beta y_{\mathrm{s}\mathrm{t}}{x}_{\mathrm{1}}^{\mathrm{2}}+\\ \beta x_{\mathrm{1}}^{\mathrm{2}}{y}_{\mathrm{1}}+\mathrm{3}\beta y_{\mathrm{s}\mathrm{t}}{y}_{\mathrm{1}}^{\mathrm{2}}+\beta y_{\mathrm{1}}^{\mathrm{3}}-F_{\mathrm{2}}=m{e}_{\mathrm{d}}{\Omega }^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}(\Omega t+\gamma )\end{array}\right\}（3）

其中，c₁、c₂为阻尼系数，e_d为圆盘的偏心率，F₁、F₂为转子旋转过程中呼吸裂纹对转轴施加的参数力，γ为转轴OX₂与不平衡方向Ge之间的转角，如图1b所示。不难看出，式（3）等号左端第三项中两个方向的刚度系数不同，这是由非线性刚度系数β和静挠度y_st引起的。运动方程同时包含平方和立方非线性项，其中平方非线性是由静挠度引起的。在重力占优的情况下裂纹转子旋转时，裂纹面随转角变化而周期性地开闭呈现呼吸行为，即用f（Ωt）将裂纹呼吸描述为圆盘转角的周期函数，其表达式如下：

f(\Omega t)=(\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}(\Omega t))/\mathrm{2}（4）

将式（4）代入式（3）可得

\begin{array}{c}m{\ddot{x}}_{\mathrm{1}}+c_{\mathrm{1}}{\dot{x}}_{\mathrm{1}}+(k+\beta y_{\mathrm{s}\mathrm{t}}^{\mathrm{2}})x_{\mathrm{1}}+\mathrm{2}\beta y_{\mathrm{s}\mathrm{t}}{x}_{\mathrm{1}}{y}_{\mathrm{1}}+\\ \beta y_{\mathrm{1}}^{\mathrm{2}}{x}_{\mathrm{1}}+\beta x_{\mathrm{1}}^{\mathrm{3}}=m{e}_{\mathrm{d}}{\Omega }^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}(\Omega t+\gamma )+\\ \frac{\mathrm{\Delta }k}{\mathrm{4}}(-\frac{\mathrm{1}}{\mathrm{2}}{x}_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}(\Omega t)-\frac{\mathrm{1}}{\mathrm{2}}{y}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}(\Omega t))+\\ \frac{\mathrm{\Delta }k}{\mathrm{4}}(y_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\Omega t)+x_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\Omega t))+\end{array}

\frac{\mathrm{\Delta }k}{\mathrm{4}}(\frac{\mathrm{1}}{\mathrm{2}}{y}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{3}\Omega t)-\frac{\mathrm{1}}{\mathrm{2}}{x}_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{3}\Omega t))（5）

\begin{array}{c}m{\ddot{y}}_{\mathrm{1}}+c_{\mathrm{2}}{\dot{y}}_{\mathrm{1}}+(k+\mathrm{3}\beta y_{\mathrm{s}\mathrm{t}}^{\mathrm{2}})y_{\mathrm{1}}+\beta y_{\mathrm{s}\mathrm{t}}{x}_{\mathrm{1}}^{\mathrm{2}}+\beta x_{\mathrm{1}}^{\mathrm{2}}{y}_{\mathrm{1}}+\\ \mathrm{3}\beta y_{\mathrm{s}\mathrm{t}}{y}_{\mathrm{1}}^{\mathrm{2}}+\beta y_{\mathrm{1}}^{\mathrm{3}}=m{e}_{\mathrm{d}}{\Omega }^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}(\Omega t+\gamma )+\\ \frac{\mathrm{\Delta }k}{\mathrm{4}}(-\frac{\mathrm{1}}{\mathrm{2}}{x}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}(\Omega t)-\frac{\mathrm{3}}{\mathrm{2}}{y}_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}(\Omega t))+\\ \frac{\mathrm{\Delta }k}{\mathrm{4}}(x_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\Omega t)-y_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\Omega t))+\end{array}

\frac{\mathrm{\Delta }k}{\mathrm{4}}(\frac{\mathrm{1}}{\mathrm{2}}{y}_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{3}\Omega t)+\frac{\mathrm{1}}{\mathrm{2}}{x}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{3}\Omega t))（6）

其中，Δk为裂纹引起的线性刚度系数k的变化，令

\begin{array}{c}{c}_x=c_{\mathrm{1}}/\sqrt[]{km}{c}_y=c_{\mathrm{2}}/\sqrt[]{km}E=e_{\mathrm{d}}/y_{\mathrm{s}\mathrm{t}}{\beta }_{\mathrm{2}}=\beta y_{\mathrm{s}\mathrm{t}}^{\mathrm{2}}/k\\ \lambda =\mathrm{\Delta }k/k\omega =\Omega /\sqrt[]{k/m}{\omega }_x^{\mathrm{2}}=\mathrm{1}+{\beta }_{\mathrm{2}}-\lambda /\mathrm{4}\\ {\omega }_y^{\mathrm{2}}=\mathrm{1}+\mathrm{3}{\beta }_{\mathrm{2}}-\lambda /\mathrm{4}\end{array}

其中，λ为裂纹深度因子。则可得如下量纲一化的动力学方程：

\begin{array}{c}{\ddot{x}}_{\mathrm{1}}+c_x{\dot{x}}_{\mathrm{1}}+{\omega }_x^{\mathrm{2}}{x}_{\mathrm{1}}+\mathrm{2}{\beta }_{\mathrm{2}}{x}_{\mathrm{1}}{y}_{\mathrm{1}}+{\beta }_{\mathrm{2}}(x_{\mathrm{1}}^{\mathrm{3}}+x_{\mathrm{1}}{y}_{\mathrm{1}}^{\mathrm{2}})=\\ E{\omega }^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}(\omega t+\gamma )+\\ \frac{\mathrm{1}}{\mathrm{8}}\lambda x_{\mathrm{1}}(-\mathrm{s}\mathrm{i}\mathrm{n}(\omega t)+\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\omega t)-\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{3}\omega t))-\end{array}

\frac{\mathrm{1}}{\mathrm{8}}\lambda y_{\mathrm{1}}(\mathrm{c}\mathrm{o}\mathrm{s}(\omega t)-\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\omega t)-\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{3}\omega t))（7）

\begin{array}{c}{\ddot{y}}_{\mathrm{1}}+c_y{\dot{y}}_{\mathrm{1}}+{\omega }_y^{\mathrm{2}}{x}_{\mathrm{1}}+{\beta }_{\mathrm{2}}(x_{\mathrm{1}}^{\mathrm{2}}+\mathrm{3}{y}_{\mathrm{1}}^{\mathrm{2}})+{\beta }_{\mathrm{2}}(y_{\mathrm{1}}^{\mathrm{3}}+y_{\mathrm{1}}{x}_{\mathrm{1}}^{\mathrm{2}})=\\ E{\omega }^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}(\omega t+\gamma )-\\ \frac{\mathrm{1}}{\mathrm{8}}\lambda x_{\mathrm{1}}(\mathrm{c}\mathrm{o}\mathrm{s}(\omega t)-\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\omega t)-\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{3}\omega t))+\end{array}

\frac{\mathrm{1}}{\mathrm{8}}\lambda y_{\mathrm{1}}(-\mathrm{3}\mathrm{s}\mathrm{i}\mathrm{n}(\omega t)-\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\omega t)+\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{3}\omega t))（8）

2 多尺度法的摄动分析

对式（7）、式（8）进行重刻度，并引入三个时间尺度：

\left.\begin{array}{l}{c}_x\leftrightarrow {\epsilon }^{\mathrm{2}}{c}_x{c}_y\leftrightarrow {\epsilon }^{\mathrm{2}}{c}_y\lambda \leftrightarrow {\epsilon }^{\mathrm{2}}\lambda E\leftrightarrow {\epsilon }^{\mathrm{3}}E\\ T_{\mathrm{0}}=t{T}_{\mathrm{1}}=\epsilon t{T}_{\mathrm{2}}={\epsilon }^{\mathrm{2}}t\end{array}\right\}（9）

其中，ε为量纲一小参数。假设式（7）、式（8）的解为

\left.\begin{array}{l}{x}_{\mathrm{1}}(t)=\epsilon x_{\mathrm{11}}(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}})+{\epsilon }^{\mathrm{2}}{x}_{\mathrm{12}}(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}})+\\ {\epsilon }^{\mathrm{3}}{x}_{\mathrm{13}}(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}})\\ y_{\mathrm{1}}(t)=\epsilon y_{\mathrm{11}}(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}})+{\epsilon }^{\mathrm{2}}{y}_{\mathrm{12}}(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}})+\\ {\epsilon }^{\mathrm{3}}{y}_{\mathrm{13}}(T_{\mathrm{0}},T_{\mathrm{1}},T_{\mathrm{2}})\end{array}\right\}（10）

将式（9）、式（10）代入式（7）、式（8），并对比同幂次系数，包括ε项：

\left.\begin{array}{l}{D}_{\mathrm{0}}^{\mathrm{2}}{x}_{\mathrm{11}}+{\omega }_x^{\mathrm{2}}{x}_{\mathrm{11}}=\mathrm{0}\\ D_{\mathrm{0}}^{\mathrm{2}}{y}_{\mathrm{11}}+{\omega }_y^{\mathrm{2}}{y}_{\mathrm{11}}=\mathrm{0}\end{array}\right\}（11）

和ε²项：

\left.\begin{array}{l}{D}_{\mathrm{0}}^{\mathrm{2}}{x}_{\mathrm{12}}+{\omega }_x^{\mathrm{2}}{x}_{\mathrm{12}}=-\mathrm{2}{D}_{\mathrm{0}}{D}_{\mathrm{1}}{x}_{\mathrm{11}}-\mathrm{2}{\beta }_{\mathrm{2}}{x}_{\mathrm{11}}{y}_{\mathrm{11}}\\ D_{\mathrm{0}}^{\mathrm{2}}{y}_{\mathrm{12}}+{\omega }_y^{\mathrm{2}}{y}_{\mathrm{12}}=-\mathrm{2}{D}_{\mathrm{0}}{D}_{\mathrm{1}}{y}_{\mathrm{11}}-{\beta }_{\mathrm{2}}{x}_{\mathrm{11}}^{\mathrm{2}}-\mathrm{3}{\beta }_{\mathrm{2}}{y}_{\mathrm{11}}^{\mathrm{2}}\end{array}\right\}（12）

以及ε³项：

\left.\begin{array}{l}{D}_{\mathrm{0}}^{\mathrm{2}}{x}_{\mathrm{13}}+{\omega }_x^{\mathrm{2}}{x}_{\mathrm{13}}=-\mathrm{2}{D}_{\mathrm{0}}{D}_{\mathrm{1}}{x}_{\mathrm{12}}-(D_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{D}_{\mathrm{2}}{D}_{\mathrm{0}})x_{\mathrm{11}}-\\ c_{\mathrm{1}}{D}_{\mathrm{0}}{x}_{\mathrm{11}}-\mathrm{2}{\beta }_{\mathrm{2}}{x}_{\mathrm{12}}{y}_{\mathrm{11}}-\mathrm{2}{\beta }_{\mathrm{2}}{y}_{\mathrm{12}}{x}_{\mathrm{11}}-{\beta }_{\mathrm{2}}{y}_{\mathrm{11}}^{\mathrm{2}}{x}_{\mathrm{11}}-{\beta }_{\mathrm{2}}{x}_{\mathrm{11}}^{\mathrm{2}}+\\ \frac{\mathrm{1}}{\mathrm{8}}\lambda x_{\mathrm{11}}(-\mathrm{s}\mathrm{i}\mathrm{n}(\omega t)+\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\omega t)-\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{3}\omega t))-\\ \frac{\mathrm{1}}{\mathrm{8}}\lambda y_{\mathrm{11}}(\mathrm{c}\mathrm{o}\mathrm{s}(\omega t)-\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\omega t)-\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{3}\omega t))+\\ E{\omega }^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}(\omega t+\gamma )\\ D_{\mathrm{0}}^{\mathrm{2}}{y}_{\mathrm{13}}+{\omega }_y^{\mathrm{2}}{y}_{\mathrm{13}}=-\mathrm{2}{D}_{\mathrm{0}}{D}_{\mathrm{1}}{y}_{\mathrm{12}}-(D_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{D}_{\mathrm{2}}{D}_{\mathrm{0}})y_{\mathrm{11}}-\\ c_{\mathrm{2}}{D}_{\mathrm{0}}{y}_{\mathrm{11}}-\mathrm{2}{\beta }_{\mathrm{2}}{x}_{\mathrm{12}}{x}_{\mathrm{11}}-\mathrm{6}{\beta }_{\mathrm{2}}{y}_{\mathrm{12}}{y}_{\mathrm{11}}-{\beta }_{\mathrm{2}}{x}_{\mathrm{11}}^{\mathrm{2}}{y}_{\mathrm{11}}-{\beta }_{\mathrm{2}}{y}_{\mathrm{11}}^{\mathrm{2}}-\\ \frac{\mathrm{1}}{\mathrm{8}}\lambda x_{\mathrm{11}}(\mathrm{c}\mathrm{o}\mathrm{s}(\omega t)-\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\omega t)-\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{3}\omega t))+\\ \frac{\mathrm{1}}{\mathrm{8}}\lambda y_{\mathrm{11}}(-\mathrm{3}\mathrm{c}\mathrm{o}\mathrm{s}(\omega t)-\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\omega t)+\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{3}\omega t))+\\ E{\omega }^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}(\omega t+\gamma )\end{array}\right\}（13）

其中，D₀表示对T₀求导，D₁表示对T₁求导，D₂表示对T₂求导，则式（11）的通解为

\left.\begin{array}{l}{x}_{\mathrm{11}}=A_x(T_{\mathrm{1}},T_{\mathrm{2}})\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_x{T}_{\mathrm{0}})+V_{\mathrm{c}\mathrm{c}}\\ y_{\mathrm{11}}=A_y(T_{\mathrm{1}},T_{\mathrm{2}})\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_y{T}_{\mathrm{0}})+V_{\mathrm{c}\mathrm{c}}\end{array}\right\}（14）

其中，A_x （T₁，T₂）、A_y （T₁，T₂）分别表示T₁、T₂的复函数，V_cc表示前面各项的共轭复数（complex conjugate）。将式（14）代入式（12）可得

\left.\begin{array}{l}{D}_{\mathrm{0}}^{\mathrm{2}}{x}_{\mathrm{12}}+{\omega }_x^{\mathrm{2}}{x}_{\mathrm{12}}=-\mathrm{2}\mathrm{i}{\omega }_x{D}_{\mathrm{1}}{A}_x\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_x{T}_{\mathrm{0}})-\\ \mathrm{2}{\beta }_{\mathrm{2}}(A_x{A}_y\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}({\omega }_x+{\omega }_y)T_{\mathrm{0}}))-\\ \mathrm{2}{\beta }_{\mathrm{2}}({\overline{A}}_x{A}_y\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(-{\omega }_x+{\omega }_y)T_{\mathrm{0}}))+V_{\mathrm{c}\mathrm{c}}\\ D_{\mathrm{0}}^{\mathrm{2}}{y}_{\mathrm{12}}+{\omega }_y^{\mathrm{2}}{y}_{\mathrm{12}}=-\mathrm{2}\mathrm{i}{\omega }_y{D}_{\mathrm{1}}{A}_y\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_y{T}_{\mathrm{0}})-\\ {\beta }_{\mathrm{2}}(A_x^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{2}\mathrm{i}{\omega }_x{T}_{\mathrm{0}})+A_x{\overline{A}}_x)-\\ \mathrm{3}{\beta }_{\mathrm{2}}(A_y^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{2}\mathrm{i}{\omega }_y{T}_{\mathrm{0}})+A_y{\overline{A}}_y)+V_{\mathrm{c}\mathrm{c}}\end{array}\right\}（15）

其中，{\overline{A}}_x、{\overline{A}}_y分别为A_x和A_y对应的共轭项，并得到ε²中的长期项为

D_{\mathrm{1}}{A}_x=\mathrm{0}{D}_{\mathrm{1}}{A}_y=\mathrm{0}（16）

这意味着复函数A_x 、A_y 是仅依赖于T₂的函数，可用极坐标形式重写A_x 、A_y，其表达式如下：

\left.\begin{array}{l}{A}_x=\frac{\mathrm{1}}{\mathrm{2}}{a}_x{T}_{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\phi }_x{T}_{\mathrm{2}})\\ A_y=\frac{\mathrm{1}}{\mathrm{2}}{a}_y{T}_{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\phi }_y{T}_{\mathrm{2}})\end{array}\right\}（17）

式中：a_x 、a_y 分别为横向和竖向的振幅；φ_x 、φ_y 分别为横向和竖向的角频率。

将式（14）、式（16）代入式（12）得到解为

\left.\begin{array}{l}{x}_{\mathrm{12}}=\frac{\mathrm{2}{\beta }_{\mathrm{2}}}{{\omega }_y}(\frac{{\overline{A}}_x{A}_y}{{\omega }_y-\mathrm{2}{\omega }_x}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{T}_{\mathrm{0}}({\omega }_y-{\omega }_x)))+\\ \frac{\mathrm{2}{\beta }_{\mathrm{2}}}{{\omega }_y}(\frac{A_x{A}_y}{{\omega }_y+\mathrm{2}{\omega }_x}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{T}_{\mathrm{0}}({\omega }_y+{\omega }_x)))+V_{\mathrm{c}\mathrm{c}}\\ y_{\mathrm{12}}=\frac{{\beta }_{\mathrm{2}}{A}_x^{\mathrm{2}}}{\mathrm{4}{\omega }_x^{\mathrm{2}}-{\omega }_y^{\mathrm{2}}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{2}\mathrm{i}{T}_{\mathrm{0}}{\omega }_x)+\frac{{\beta }_{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}}{A}_y^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{2}\mathrm{i}{T}_{\mathrm{0}}{\omega }_y)-\\ \frac{{\beta }_{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}}(A_x{\overline{A}}_x+\mathrm{3}{A}_y{\overline{A}}_y)+V_{\mathrm{c}\mathrm{c}}\end{array}\right\}（18）

将式（14）、式（18）代入式（13）可得

\begin{array}{c}{D}_{\mathrm{0}}^{\mathrm{2}}{x}_{\mathrm{13}}+{\omega }_x^{\mathrm{2}}{x}_{\mathrm{13}}=(-\mathrm{2}\mathrm{i}{\omega }_x{D}_{\mathrm{2}}{A}_x-\mathrm{i}{c}_x{\omega }_x{A}_x)\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{T}_{\mathrm{0}}{\omega }_x)-\\ (\mathrm{3}{\beta }_{\mathrm{2}}+\frac{\mathrm{2}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{4}{\omega }_x^{\mathrm{2}}-{\omega }_y^{\mathrm{2}}}-\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}})A_x^{\mathrm{2}}{\overline{A}}_x\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{T}_{\mathrm{0}}{\omega }_x)-\\ [{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y-\mathrm{2}{\omega }_x)}+\frac{\mathrm{2}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}}]A_x^{\mathrm{2}}{\overline{A}}_x\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{T}_{\mathrm{0}}(\mathrm{2}{\omega }_y-{\omega }_x))-\\ [{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y+\mathrm{2}{\omega }_x)}+\frac{\mathrm{2}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}}]A_x^{\mathrm{2}}{A}_x\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{T}_{\mathrm{0}}(\mathrm{2}{\omega }_y+{\omega }_x))-\\ ({\beta }_{\mathrm{2}}+\frac{\mathrm{2}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{4}{\omega }_x^{\mathrm{2}}-{\omega }_y^{\mathrm{2}}})A_x^{\mathrm{3}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{3}\mathrm{i}{T}_{\mathrm{0}}{\omega }_y)+\\ \frac{\mathrm{1}}{\mathrm{2}}E{\omega }^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega T_{\mathrm{0}}+\gamma ))-(\mathrm{2}{\beta }_{\mathrm{2}}+\frac{\mathrm{8}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}-\mathrm{4}{\omega }_x^{\mathrm{2}}}-\frac{\mathrm{12}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}})·\\ A_x{A}_y{\overline{A}}_y\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_x{T}_{\mathrm{0}})+\frac{\mathrm{1}}{\mathrm{16}}\lambda A_x(\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega +{\omega }_x)T_{\mathrm{0}})+\\ \mathrm{2}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega +{\omega }_x)T_{\mathrm{0}})+\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega +{\omega }_x)T_{\mathrm{0}}))+\\ \frac{\mathrm{1}}{\mathrm{16}}\lambda {\overline{A}}_x(\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega -{\omega }_x)T_{\mathrm{0}})+\mathrm{2}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega -{\omega }_x)T_{\mathrm{0}})+\\ \mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega -{\omega }_x)T_{\mathrm{0}}))-\frac{\mathrm{1}}{\mathrm{16}}\lambda A_y(\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega +{\omega }_y)T_{\mathrm{0}})+\\ \mathrm{2}\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega +{\omega }_y)T_{\mathrm{0}})-\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega +{\omega }_y)T_{\mathrm{0}}))-\\ \frac{\mathrm{1}}{\mathrm{16}}\lambda {\overline{A}}_y(\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega -{\omega }_y)T_{\mathrm{0}})+\mathrm{2}\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega -{\omega }_y)T_{\mathrm{0}})-\end{array}

\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega -{\omega }_y)T_{\mathrm{0}}))+V_{\mathrm{c}\mathrm{c}}（19）

\begin{array}{c}{D}_{\mathrm{0}}^{\mathrm{2}}{y}_{\mathrm{13}}+{\omega }_y^{\mathrm{2}}{y}_{\mathrm{13}}=(-\mathrm{2}\mathrm{i}{\omega }_y{D}_{\mathrm{2}}{A}_y-\mathrm{i}{c}_y{\omega }_y{A}_y)\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_y{T}_{\mathrm{0}})-\\ (\mathrm{3}{\beta }_{\mathrm{2}}-\frac{\mathrm{30}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}})A_y^{\mathrm{2}}{\overline{A}}_y\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_y{T}_{\mathrm{0}})-({\beta }_{\mathrm{2}}+\frac{\mathrm{6}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}})A_y^{\mathrm{3}}·\\ \mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\mathrm{3}{\omega }_y{T}_{\mathrm{0}})-[{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y-\mathrm{2}{\omega }_x)}+\frac{\mathrm{6}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{4}{\omega }_x^{\mathrm{2}}-{\omega }_y^{\mathrm{2}}}]·\\ A_x^{\mathrm{2}}{\overline{A}}_y\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}{\omega }_x-{\omega }_y)T_{\mathrm{0}})-[{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y-\mathrm{2}{\omega }_x)}+\\ \frac{\mathrm{6}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{4}{\omega }_x^{\mathrm{2}}-{\omega }_y^{\mathrm{2}}}]A_x^{\mathrm{2}}{A}_y\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}{\omega }_x+{\omega }_y)T_{\mathrm{0}})-\\ (\mathrm{2}{\beta }_{\mathrm{2}}+\frac{\mathrm{8}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}-\mathrm{4}{\omega }_x^{\mathrm{2}}}-\frac{\mathrm{12}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}})A_y{A}_x{\overline{A}}_x\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}{\omega }_y{T}_{\mathrm{0}})-\end{array}

\begin{array}{c}\frac{\mathrm{1}}{\mathrm{16}}\lambda A_x(\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega +{\omega }_x)T_{\mathrm{0}})+\mathrm{2}\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega +{\omega }_x)T_{\mathrm{0}})-\\ \mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega +{\omega }_x)T_{\mathrm{0}}))-\frac{\mathrm{1}}{\mathrm{16}}\lambda A_x(\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega -{\omega }_x)T_{\mathrm{0}})+\\ \mathrm{2}\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega -{\omega }_x)T_{\mathrm{0}})-\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega -{\omega }_x)T_{\mathrm{0}}))+\\ \frac{\mathrm{1}}{\mathrm{16}}\lambda A_y(\mathrm{3}\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega +{\omega }_y)T_{\mathrm{0}})-\mathrm{2}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega +{\omega }_y)T_{\mathrm{0}})-\\ \mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega +{\omega }_y)T_{\mathrm{0}}))+\frac{\mathrm{1}}{\mathrm{16}}\lambda {\overline{A}}_y(\mathrm{3}\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega -{\omega }_y)T_{\mathrm{0}})-\\ \mathrm{2}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{2}\omega -{\omega }_y)T_{\mathrm{0}})-\mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}\omega -{\omega }_y)T_{\mathrm{0}}))-\end{array}

\frac{\mathrm{1}}{\mathrm{2}}\mathrm{i}E{\omega }^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\omega T_{\mathrm{0}}+\gamma ))+V_{\mathrm{c}\mathrm{c}}（20）

本文讨论1:1内共振条件下的2/3超谐波共振，引入解谐参数σ₁和σ₂，用σ₁表示ω_y 接近ω_x 的程度，用σ₂表示ω接近2ω_x /3的程度，则有

\left.\begin{array}{l}{\omega }_y={\omega }_x+{\epsilon }^{\mathrm{2}}{\sigma }_{\mathrm{1}}\\ \omega =\frac{\mathrm{2}}{\mathrm{3}}{\omega }_x+{\epsilon }^{\mathrm{2}}{\sigma }_{\mathrm{2}}\end{array}\right\}（21）

将式（21）代入式（20）中，为避免长期项的出现，需要满足

\left.\begin{array}{l}\mathrm{2}\mathrm{i}{\omega }_x{D}_{\mathrm{2}}{A}_x=\mathrm{i}{c}_x{A}_x+{\Gamma }_{\mathrm{11}}{A}_x^{\mathrm{2}}{\overline{A}}_x+\\ {\Gamma }_{\mathrm{12}}{A}_y{\overline{A}}_y{A}_x+{\Gamma }_{\mathrm{13}}{\overline{A}}_x{A}_y^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{2}\mathrm{i}{\sigma }_{\mathrm{1}}{T}_{\mathrm{2}})+\\ \frac{\mathrm{1}}{\mathrm{16}}\lambda \mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{3}\mathrm{i}{\sigma }_{\mathrm{2}}{T}_{\mathrm{2}})+\frac{\mathrm{1}}{\mathrm{16}}\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}{\sigma }_{\mathrm{2}}-{\sigma }_{\mathrm{1}})T_{\mathrm{2}})\\ \mathrm{2}\mathrm{i}{\omega }_y{D}_{\mathrm{2}}{A}_y=\mathrm{i}{c}_y{A}_y+{\Gamma }_{\mathrm{21}}{A}_y^{\mathrm{2}}{\overline{A}}_y+\\ {\Gamma }_{\mathrm{22}}{A}_x{\overline{A}}_x{A}_y+{\Gamma }_{\mathrm{23}}{\overline{A}}_y{A}_x^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{2}\mathrm{i}{\sigma }_{\mathrm{1}}{T}_{\mathrm{2}})-\\ \frac{\mathrm{1}}{\mathrm{16}}\lambda \mathrm{i}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}{\sigma }_{\mathrm{2}}-\mathrm{2}{\sigma }_{\mathrm{1}})T_{\mathrm{2}})+\frac{\mathrm{1}}{\mathrm{16}}\lambda \mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}(\mathrm{3}{\sigma }_{\mathrm{2}}-{\sigma }_{\mathrm{1}})T_{\mathrm{2}})\end{array}\right\}（22）

\begin{array}{c}{\Gamma }_{\mathrm{11}}=-(\mathrm{3}{\beta }_{\mathrm{2}}+\frac{\mathrm{2}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{4}{\omega }_x^{\mathrm{2}}-{\omega }_y^{\mathrm{2}}}-\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}})\\ {\Gamma }_{\mathrm{12}}=-[\mathrm{2}{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y-\mathrm{2}{\omega }_x)}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y+\mathrm{2}{\omega }_x)}-\frac{\mathrm{12}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}}]\\ {\Gamma }_{\mathrm{13}}=-[{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y-\mathrm{2}{\omega }_x)}+\frac{\mathrm{2}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}}]\\ {\Gamma }_{\mathrm{21}}=-(\mathrm{3}{\beta }_{\mathrm{2}}-\frac{\mathrm{30}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}})\\ {\Gamma }_{\mathrm{22}}=-[\mathrm{2}{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y-\mathrm{2}{\omega }_x)}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y+\mathrm{2}{\omega }_x)}-\frac{\mathrm{12}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y^{\mathrm{2}}}]\\ {\Gamma }_{\mathrm{23}}=-[{\beta }_{\mathrm{2}}+\frac{\mathrm{4}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{{\omega }_y({\omega }_y-\mathrm{2}{\omega }_x)}+\frac{\mathrm{6}{\beta }_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{4}{\omega }_x^{\mathrm{2}}-{\omega }_y^{\mathrm{2}}}]\end{array}

将式（17）代入式（22）中分离实部和虚部，得到系统自治演化方程：

\frac{\mathrm{d}{a}_x}{\mathrm{d}{T}_{\mathrm{2}}}=-\frac{\mathrm{1}}{\mathrm{2}}{c}_x{a}_x-\frac{{\Gamma }_{\mathrm{13}}}{\mathrm{8}{\omega }_x}{a}_x{a}_y^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}({\psi }_x+{\psi }_y)+

\frac{\lambda a_y}{\mathrm{32}{\omega }_x}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{\psi }_y-{\psi }_x}{\mathrm{2}}+\frac{\lambda a_x}{\mathrm{32}{\omega }_x}\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_x（23）

\frac{\mathrm{d}{\psi }_x}{\mathrm{d}{T}_{\mathrm{2}}}=\mathrm{3}{\sigma }_{\mathrm{2}}-\frac{{\Gamma }_{\mathrm{11}}{a}_x^{\mathrm{2}}}{\mathrm{4}{\omega }_x}-\frac{{\Gamma }_{\mathrm{12}}{a}_y^{\mathrm{2}}}{\mathrm{4}{\omega }_x}-\frac{\lambda }{\mathrm{16}{\omega }_x}\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_x+

\frac{\lambda a_y}{\mathrm{16}{\omega }_x{a}_x}\mathrm{c}\mathrm{o}\mathrm{s}\frac{{\psi }_y-{\psi }_x}{\mathrm{2}}+\frac{{\Gamma }_{\mathrm{13}}{a}_y^{\mathrm{2}}}{\mathrm{8}{\omega }_x}\mathrm{c}\mathrm{o}\mathrm{s}({\psi }_x+{\psi }_y)（24）

\frac{\mathrm{d}{a}_y}{\mathrm{d}{T}_{\mathrm{2}}}=-\frac{\mathrm{1}}{\mathrm{2}}{c}_y{a}_y-\frac{{\Gamma }_{\mathrm{23}}}{\mathrm{8}{\omega }_y}{a}_y{a}_x^{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}({\psi }_x+{\psi }_y)+

\frac{\lambda a_x}{\mathrm{32}{\omega }_y}\mathrm{s}\mathrm{i}\mathrm{n}\frac{{\psi }_y-{\psi }_x}{\mathrm{2}}-\frac{\lambda a_y}{\mathrm{32}{\omega }_y}\mathrm{c}\mathrm{o}\mathrm{s}{\psi }_y（25）

\frac{\mathrm{d}{\psi }_y}{\mathrm{d}{T}_{\mathrm{2}}}=\mathrm{2}{\sigma }_{\mathrm{1}}-\mathrm{3}{\sigma }_{\mathrm{2}}+\frac{{\Gamma }_{\mathrm{21}}{a}_y^{\mathrm{2}}}{\mathrm{4}{\omega }_y}+\frac{{\Gamma }_{\mathrm{22}}{a}_x^{\mathrm{2}}}{\mathrm{4}{\omega }_y}+\frac{\lambda }{\mathrm{16}{\omega }_y}\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_y-

\frac{\lambda a_x}{\mathrm{16}{\omega }_y{a}_y}\mathrm{c}\mathrm{o}\mathrm{s}\frac{{\psi }_y-{\psi }_x}{\mathrm{2}}+\frac{{\Gamma }_{\mathrm{23}}{a}_x^{\mathrm{2}}}{\mathrm{8}{\omega }_y}\mathrm{c}\mathrm{o}\mathrm{s}({\psi }_x+{\psi }_y)（26）

{\psi }_x=\mathrm{3}{T}_{\mathrm{2}}{\sigma }_{\mathrm{2}}-\mathrm{2}{\phi }_x{\psi }_y=\mathrm{2}{T}_{\mathrm{2}}{\sigma }_{\mathrm{1}}-\mathrm{3}{T}_{\mathrm{2}}{\sigma }_{\mathrm{2}}+\mathrm{2}{\phi }_y

为研究该系统的稳态响应，可以将系统演化方程等号左侧设为0，进而求解所得到的非线性代数方程组。由此系统的稳态响应就可以通过求解该方程组得到，其解可分为两种情况：①当a_x =0、a_y =0时称为零解，用笛卡儿坐标系下的演化方程来判断零解的稳定性；②当a_x ≠0、a_y ≠0时称为非零解。由于本文讨论的共振形式对应的演化方程较为复杂，解析解难以推导，因此本文运用 Mathematic中的“Nsolve”命令寻求方程组的近似数值解，并利用李雅普诺夫稳定性理论研究演化方程的稳定性。

3 参数变化对系统频响曲线的影响

采用文献［19］由实验测定的转子系统各参数值进行分析：c_x =0.015，c_y =0.025，β₂=0.05，E=0.02，σ₂=0，λ=0.25。

图2为不同裂纹深度下系统的频率响应曲线。由图2a可知，当裂纹深度较小时系统只有零解，裂纹深度增大到一定水平后，系统超谐波响应会被激发，如图2b~图2e所示。如图2b所示，当裂纹深度因子λ达到0.2左右时，频率响应曲线由相互孤立的一条零解分支和一条凸起非零解分支组成，本文称为“频率岛”现象。随着裂纹深度的增大，“频率岛”不断扩大，当裂纹深度增大至0.25左右时，凸起的非零解的分支逐渐与零解连接在一起，“频率岛”现象消失。如图2b~图2d所示，裂纹深度的增大导致零解的共振带宽扩大，非零解的曲线分为两支，其中一支稳定，另一支不稳定。两支非零解之间的零解不稳定，其他的零解为稳定解。如图2e~图2f所示，随着裂纹深度的增大，非零解不稳定分支出现缠绕、绕圈现象，且绕圈逐渐扩大并与非零解交汇，系统的多值现象逐渐复杂。 此外由图2可见，本文讨论的裂纹转子的2/3超谐波共振在完全调谐附近有较大的响应带宽，可作为裂纹识别的选项。

图3给出了不同横向阻尼系数c_x 下系统振幅的频响曲线。当c_x 较小时，系统有零解和非零解，且非零解的曲线分为两支，不稳定分支出现了缠绕现象，多值区域较大。c_x 增大至0.02左右时，频响曲线出现了“频率岛”现象，c_x 增大至0.04左右时，“频率岛”现象消失，系统只有零解。随着阻尼系数的增大，系统响应的不稳定区域减小，双跳跃的峰值逐渐减小直到消失。可以得到c_x 的增大对系统的跳跃现象会有抑制作用。在共振区域附近，响应振幅随着c_x 增大而减小，但在较远的共振区域，响应振幅的变化并不明显。

图4给出了不同竖向阻尼系数c_y 下系统振幅的频响曲线。由图可见，c_y 对系统稳定性以及响应幅值的影响与c_x 影响大致类似。与c_x 影响不同的是，当c_y 增大到足够大时，系统总存在非零解即系统的超谐波共振总能被激发出来，受篇幅所限，这里仅展示c_{y≤0.1的结果。这种物理不对称性是由于本文考虑了重力引起的预变形效应，给系统动力学方程引入了平方非线性。}

图5显示了不同解谐参数下系统稳态响应振幅值随着裂纹深度的演变曲线。如图5a所示，完全调谐（σ₂=0）时，裂纹深度因子λ≥0.22时出现了非零解，同时零解转为不稳定。较大的解谐参数导致稳态响应的多值性变得更为复杂；而发生完全调谐共振时，稳态响应幅值没有跳跃区间，非零解的多值现象消失。

4 平方非线性和立方非线性对系统的影响

图6比较了平方非线性和立方非线性对系统超谐波共振响应的影响。由图6可见，立方非线性对裂纹转子的频率响应有硬化效应，而平方非线性有软化效应。并且立方非线性引起的硬化效应强于平方非线性引起的软化效应。平方非线性缩小了系统动力响应的不稳定区间。

5 数值验证

为验证通过多尺度法获得的系统稳态解的准确性，本节将利用Runge-Kutta法对原始系统进行足够时长的数值积分，以确保系统能够达到稳态。在此基础上，对原系统进行了正向和反向的扫频分析。如图7所示，圆圈表示正向扫频过程中的数值积分点，加号表示反向扫频时的数值积分点，同时，多尺度法得出的近似解以曲线的形式呈现。观察发现，在正向和反向扫频过程中，系统发生跳跃现象的跳跃点不同，原系统数值积分的结果与多尺度的近似解吻合性良好，验证了本文采用的多尺度法的正确性。

图8为裂纹转子在不同裂纹深度和解谐参数条件下的时域图、相图、庞加莱图和横向位移的频谱图，频谱图上包含了自由振动频率、转速频率以及两者的线性组合，其中时间T=ωt/（2π）。当固定σ₂ =0.02时（图8a~图8d），随着λ的增大，系统经历了单周期→多周期→拟周期→多周期的变化历程。由图8a中的频谱图可知，1倍频主要是由质量偏心引起的，而其他部分则主要由裂纹形成。对比图8a和图8b可知，考虑裂纹参数力的作用使得系统由单周期演变为倍周期运动，由此可见裂纹参数力开始对系统的动力学响应有影响。随着λ的增大，系统响应的高阶谐波分量被激发出来。当固定λ=0.3时（图8c、图8e、图8f），随着σ₂的增大，系统经历了多周期→拟周期→多周期的变化历程。由图8e可得，当系统完全调谐时，与裂纹引起的频率相比，质量偏心引起的1倍频比较弱，因此频谱图上包含比较丰富的裂纹引起的频率成分。当解谐参数增大至一定水平时，由质量偏心引起的1倍频起主导作用。由此可以得到系统的动力学行为显著依赖于σ₂、λ，随着这些参数的变化，系统动力学响应会展现出分岔等动力学现象。

6 结论

本文基于呼吸裂纹模型，考虑了重力引起的预变形效应，研究了Jeffcott裂纹转子在1∶1内共振条件下的2/3超谐波共振的动力学行为。采用多尺度法推导出了系统的演化方程，探讨了裂纹深度、横向及竖向阻尼系数对系统稳态振幅的影响规律，揭示了平方和立方非线性对系统稳态振幅的影响。采用数值积分对近似解进行了验证。绘制了系统时程图、相图、庞加莱图及频率能谱图随解谐系数和裂纹深度的变化。得出以下结论：

1）裂纹深度增大时频率响应曲线硬特性增强，同时响应幅值与不稳定区域增大。由于本文考虑了重力引起的预变形效应，使得阻尼系数的影响有了不对称性。

2）立方非线性对转子系统的频率响应有硬化效应，而平方非线性有软化效应，前者的效应更强。

3）本文讨论的转子模型所对应的2/3超谐波共振在完全调谐附近有较大的响应带宽，且稳态响应与圆盘偏心率和夹角无关，可作为裂纹识别的选项。

参考文献

[1] WANG C H， BROWN M W. Life Prediction Techniques for Variable Amplitude Multiaxial Fatigue［J］. Journal of Engineering Materials and Techno-logy， 1996， 118（3）：367-370.

[2] JEFFCOTT H H. The Lateral Vibration of Loaded Shafts in the Neighborhood of a Whirling Speed［J］. Philosophical Magazine Series， 1919， 37：304-314.

[3] AL-SHUDEIFAT M A. On the Finite Element Modeling of the Asymmetric Cracked Rotor［J］. Journal of Sound and Vibration， 2013， 332（11）：2795-2807.

[4] GASCH R. A Survey of the Dynamic Behavior of a Simple Rotating Shaft with a Transverse Crack ［J］. Journal of Sound and Vibration， 1993， 160（2）：313-332.

[5] MAYES I W， DAVIES W G R . et al. Analysis of the Response of a Multi-rotor-bearing System Containing a Transverse Crack in a Rotor ［J］. Journal of Sound and Vibration， 1984， 106（7）：139-145.

[6] 高建民，朱晓梅.转轴上裂纹开闭模型的研究［J］．应用力学学报，1992，9（1）：108-112.

[7] GAO Jianmin， ZHU Xiaomei. Study on the Model of the Shaft Crack Opening and Closing［J］. Chinese Journal of Applied Mechanics， 1992，9（1）：108-112.

[8] SINOU J J， LEES A W. The Influence of Cracks in Rotating Shafts［J］. Journal of Sound and Vibration， 2005，285（45）： 1015-1037.

[9] DIMAROGONAS D A， PAIPETIS S A. Analytical Methods in Rotor Dynamics［M］. London：Applied Science Publisher， 1983.

[10] 路振勇，侯磊，侯升亮，等.含裂纹故障的航空发动机转子系统动力学特性分析［J］. 振动与冲击， 2018， 37（3）：40-46.

[11] LU Zhenyong， HOU Lei， HOU Shengliang， et al. Dynamic Characteristics of an Aero-engine Rotor System with Crack Faults［J］. Journal of Vibration and Shock，2018，37（3）：40-46.

[12] 胡亮，柳亦兵，徐晓星，等. 带横向裂纹的拉杆转子非线性动力学特性研究［J］．噪声与振动控制，2016，36（5）：11-14.

[13] HU Liang， LIU Yibin， XU Xiaoxing， et al. Nonlinear Dynamic Response of a Rod Fastening Rotor with a Transverse Crack［J］. Noise and Vibration Control， 2016， 36（5）：11-14.

[14] ZHANG B， LI Y M. Nonlinear Vibration of Rotating Pre-deformed Blade with Thermal Gradient［J］. Nonlinear Dynamics， 2016， 86 （1）：459-478.

[15] JUN O S， EUN H J， EARMME Y Y， et al. Modelling and Vibration Analysis of a Simple Rotor with Breathing Crack［J］.Journal of Sound and Vibration， 1992， 155 （2）：273-290.

[16] XIE J， ZI Y， CHENG W， et al. Mechanism Explanation and Experimental Verification of a New Modulation Frequency Characteristic in a Disturbed Crack Rotor［J］. Nonlinear Dynamics，2019，95（1）：597-616.

[17] 李常有，徐敏强，郭耸， 基于有限元的横向裂纹转子系 统 的 动 力 学 分 析［J］. 振 动 工 程 学 报，2009，22（5）：486-491.

[18] LI Changyou， XU Minqiang， GUO Song， et al. Dynamic Analysis of Transversely Cracked Rotor System Based on Finite Element Method［J］. Journal of Vibration Engineering， 2009，22（5）：486-491.

[19] CHEN C， DAI L， FU Y. Nonlinear Response and Dynamic Stability of a Cracked Rotor［J］. Communications in Nonlinear Science and Numerical Simulation，2007，12（6）： 1023-1037.

[20] HOU L， CHEN Y， LU Z， et al. Bifurcation Analysis for 2∶1 and 3∶1 Super-harmonic Resonances of an Aircraft Cracked Rotor System due to Maneuver Load［J］. Nonlinear Dynamics， 2015， 81（9）：531-547.

[21] 王宗勇，林伟，闻邦椿.质量突变转子系统的动力学研究［J］.中国机械工程，2009，20（13）：1586-1589.

[22] WANG Zongyong， GONG Bin， WEN Bangchun. Dynamic Analysis of Rotor System with Abrupt Quality［J］. China Aircraft Engineering， 2009， 20（13）：1586-1589.

[23] SINOU J J， LEES A W. A Non-linear Study of a Cracked Rotor［J］. European Journal of Mechanics—A/Solids， 2007， 26（1）： 152-170.

[24] YABUNO H， KASHIMURA T， INOUE T， et al. Nonlinear Normal Modes and Primary Resonance of Horizontally Supported Jeffcott Rotor［J］. Nonlinear Dynamics， 2011， 66 （3）：377-387.

[25] SAEED N A， EISSA M. Bifurcation Analysis of a Transversely Cracked Nonlinear Jeffcott Rotor System at Different Resonance Cases［J］. International Journal of Nonlinear Mechanics， 2017，58（7）：1924-1934.

[26] 何成兵， 顾煜炯， 宋光雄. 裂纹转子弯扭耦合振动非线性特性分析［J］. 振动与冲击， 2012， 31（9）：33-38.

[27] HE Chenbing， GU Yutong， SONG Guangxiong. Analysis of Nonlinear Characteristics of Flexural and Torsional Coupled Vibration of Cracked Rotor［J］. Journal of Vibration and Shock，2012，31（9）：33-38.

[28] 杨丹，甘春标，杨世锡，等.含横向裂纹Jeffcott转子刚度及动力学特性研究［J］.振动与冲击，2012，31（15）：121-126.

[29] YANG Dan， GAN Chunbiao， YANG Shixi， et al. Stiffness and Dynamical Behavior of Jeffcott Rotor with Cross Crack［J］.Journal of Vibration and Shock，2012，31（15）：121-126.

[30] YAMAMOTO T， ISHIDA Y. Linear and Nonlinear Rotor Dynamics［M］. New York： Wiley，2001.