轴承转子齿轮系统非线性动力学MATLAB仿真与故障特征分析

做旋转机械动力学仿真的人，迟早会撞上一整套连环问题：轴承转子系统怎么建模、齿轮传动的时变刚度怎么处理、裂纹故障怎么引入、非线性振动算出来之后怎么判断它是周期解还是混沌。这几个问题单独拎出来每一个都有大量文献，但真正落到MATLAB里跑通一条完整链路，从建模到数值求解，再到庞加莱截面和分岔图，能少走弯路的人并不多。这篇内容就是把我自己在这条线上踩过的坑、验证过的做法和可以直接抄的代码框架整理出来，给同样在做轴承转子动力学、齿轮动力学和非线性振动分析的朋友一个可复现的起点。

先说清楚这套东西能干什么。轴承-转子-齿轮构成的旋转机械是汽轮机、航空发动机、风电齿轮箱、压缩机这些设备的核心，而裂纹、磨损、齿面剥落这类故障本质上都会改变系统的刚度和阻尼特性，让系统表现出明显的非线性行为。传统的线性振动分析只能给出固有频率和振型，一旦系统里存在齿侧间隙、滚动轴承游隙、裂纹呼吸效应这些强非线性因素，线性方法就失效了，必须上非线性动力学的工具。MATLAB在这里的优势非常直接：ODE求解器成熟、矩阵运算快、画图方便，从建立运动方程到输出庞加莱截面，全程可以在一个脚本里完成，不需要在多个软件之间倒腾数据。这套方法适合机械工程、力学专业的研究生，也适合在企业里做故障诊断算法预研的工程师。

1. 项目定位与整体设计思路

1.1 研究背景：为什么必须上非线性分析

轴承转子系统和齿轮传动系统在理想工况下可以近似看成线性系统，但实际运行中处处是非线性源。滚动轴承的游隙导致滚珠与滚道之间的接触力呈现分段特性，间隙大的时候转子处于自由状态，间隙小的时候又是强接触，这种“时有时无”的约束本身就是典型的非光滑非线性。齿轮传动更明显：齿侧间隙是制造和安装必然存在的，啮合过程中轮齿交替承载导致啮合刚度随时间周期性变化，再加上齿面摩擦、润滑膜动态行为，齿轮副的运动方程里天然包含强非线性项。

如果系统里再出现裂纹故障，情况就更复杂。裂纹在旋转过程中会周期性张开和闭合，也就是所谓的“呼吸效应”，导致转轴的等效刚度在每个旋转周期内变化两次甚至更多次。裂纹深度不同、位置不同，对刚度的削弱程度也不同。这种刚度周期性变化会让系统在某些转速下产生次谐波共振、超谐波共振，甚至直接进入混沌运动。这些现象在频域里的表现就是边频带丰富、分数谐波成分突出，跟正常轴承故障的特征完全不一样。

所以说这个课题的本质不是“用MATLAB解一个微分方程”，而是“建立一套能够反映真实故障特征的动力学模型，然后用合适的数值方法把非线性行为完整地算出来，最后用庞加莱截面这类工具把混沌运动识别出来”。模型是地基，数值方法是工具，庞加莱截面是判据，三者缺一不可。

1.2 方案选型逻辑：代码方案为什么优于Simulink

很多初学者拿到这类问题第一反应是开Simulink搭积木，拖几个积分模块、加几个函数模块就把框图搭起来了。我的建议是：做科研分析、要扫参数、要画分岔图，老老实实写m脚本，不要用Simulink。原因很实在。

Simulink适合做控制系统仿真，因为它天然处理的是信号流和反馈回路。但动力学方程里那些分段函数、时变刚度、积分变量的相互耦合，在Simulink里表达起来特别别扭。比如齿侧间隙的分段函数，在Simulink里需要用Switch模块配合逻辑判断搭一整套东西，改一个参数要在好几个对话框里跳来跳去。而写脚本的话，分段函数就是几行if语句的事，改刚度表达式也就改一行代码。

更重要的是扫参和分岔图计算。分岔图需要让某个参数（通常是转速或激励频率）在很大范围内连续变化，每一个参数点都要跑完完整的时域积分、去掉瞬态段、提取稳态响应。这个流程在Simulink里几乎没法高效实现，得反复启动停止仿真模型，非常浪费时间。而脚本里一个for循环就能完成全部分岔计算，再配合parfor做并行加速，效率差距是数量级的。

我自己的习惯是：所有理论研究、参数扫描、特征提取全部用m脚本完成；只有当最终需要做实时仿真或者硬件在环测试时，才会把已验证过的模型转成Simulink。这条经验在几个项目里都验证过，效率最高。

1.3 项目整体分析流程

整个分析链路可以分成五步。第一步是做系统的物理建模，把轴承转子系统和齿轮传动系统的运动微分方程写出来，明确系统的自由度、激励形式和非线性项；第二步是方程的量纲一化处理，把物理参数转换成无量纲参数，这一步非常关键，它决定数值求解的稳定性和结果的可比性；第三步是编写MATLAB数值积分脚本，求解系统的时域响应；第四步是在时域响应的基础上计算庞加莱截面、分岔图、相轨迹、FFT频谱等特征量；第五步是改变关键参数，观察系统从周期运动到拟周期再到混沌的演化路径，提炼故障特征。

这五步里面，第一步和第二步是最耗精力的，也是决定结果对不对的核心。后面具体展开讲每一个环节的实现细节。

2. 动力学建模与非线性因素解析

2.1 轴承转子系统的基础模型

轴承转子系统里面有大量的非线性因素，而建立模型的核心是决定如何简化和取舍。如果做的是单盘转子，最简单的模型是Jeffcott转子模型：把转轴简化成无质量弹性轴，把轴上的圆盘简化成一个集中质量点，用两个互相垂直的横向位移$ x$ 和 $ y$ 表示转子的运动。运动方程为：

$$m\ddot{x} + c\dot{x} + k_{xx}x + k_{xy}y = F_x$$$$m\ddot{y} + c\dot{y} + k_{yx}x + k_{yy}y = F_y$$

其中 $k_{xy}$ 和 $k_{yx}$ 是交叉刚度项，来源于油膜轴承或滚动轴承的各向异性。到了这一步，线性问题已经可以求解了，但真正要做非线性分析的话，$F_x$ 和 $F_y$ 里必须包含轴承非线性的贡献。

滚动轴承的非线性主要是Hertz接触力和游隙。根据Hertz接触理论，滚珠与滚道之间的接触力与变形量的3/2次方成正比。如果轴承存在径向游隙 $\delta$，那么只有当径向变形超过游隙的时候才产生接触力，否则接触力为零。所以轴承反力可以写成：

$$F_b = k_b (r - \delta)^{3/2} \cdot \mathrm{sgn}(r - \delta)$$

其中 $r = \sqrt{x^2 + y^2}$ 是轴心位移幅值。这个式子天然是分段非线性的，而且指数是3/2，不是整数次幂，给数值求解带来了一定的挑战。更重要的是，这个非线性项对轴心位移非常敏感：转子小幅振动的时候，轴承力时有时无，系统就体现出松紧交替的特性，大幅振动的时候，接触力急剧增大，相当于刚度突然变硬。这种“软-硬交替”是非线性转子系统最典型的动力学行为来源之一。

我建议建模时给滚动轴承加一个变刚度系数，因为滚珠在滚道内的位置不断变化导致轴承刚度随保持架转角周期性波动，这个波动频率叫变柔度频率（VC频率）。把变柔度项写进刚度表达式里，模型的非线性程度会更接近真实情况，庞加莱截面上会出现更丰富的频率成分，也更有利于后面的故障特征分析。

2.2 齿轮副的时变啮合刚度与齿侧间隙

齿轮动力学里最核心的量是时变啮合刚度。一对齿轮啮合时，同时参与啮合的轮齿对数在单齿啮合区和双齿啮合区之间交替变化，导致啮合刚度随着啮合位置周期性波动。啮合频率是

$$f_m = \frac{z_1 n_1}{60}$$

其中 $z_1$ 是小齿轮齿数，$n_1$ 是小齿轮转速（r/min）。时变啮合刚度 $k_m(t)$ 可以近似写成一个平均刚度加一个谐波波动项的形式：

$$k_m(t) = k_{m0} + \sum_{i=1}^{n} k_i \cos(i \omega_m t + \phi_i)$$

工程中通常取前几阶谐波就够用了。$k_{m0}$ 是平均啮合刚度，$k_1$ 是一阶波动幅值，通常占平均刚度的10%到30%。如果精度要求不高，也可以用方波函数近似，因为本质上啮合刚度变化的周期性远大于具体波形的细节。

齿侧间隙引入的非线性是强非线性的分段函数。设齿轮副的相对位移为 $\delta = x_1 - x_2$，齿侧间隙为 $2b$，那么非线性间隙函数为

$$f(\delta) = \begin{cases} \delta - b, & \delta > b \ 0, & -b \leq \delta \leq b \ \delta + b, & \delta < -b \end{cases}$$

这个分段函数反映的是齿轮传动中的“脱啮-啮合-冲击”过程。当振动幅值小于间隙时，主从动齿轮根本不接触，力无法传递；当超过间隙后突然撞击，产生冲击载荷。这个非光滑特性会让系统在某些转速下非常容易进入混沌。把时变刚度、齿侧间隙和齿轮阻尼放在一起，齿轮副的动力学方程就是典型的参数激励非线性系统，其运动行为要比固定刚度的线性系统复杂得多。

2.3 裂纹故障的等效建模方法

齿轮裂纹和转轴裂纹在建模上有本质区别。转轴裂纹的建模通常采用呼吸裂纹模型，其核心思想是裂纹在旋转过程中随着应力状态周期性张开和闭合，导致转轴的截面惯性矩在一个旋转周期内发生变化，进而引起刚度变化。对于圆截面转轴，假设存在深度为 $a$ 的横向裂纹，裂纹处的截面惯性矩降低，等效刚度可表示为

$$k_{crack}(\theta) = k_0 \left[ 1 - \Delta \bar{k} \cdot g(\theta) \right]$$

其中 $\Delta \bar{k}$ 是平均刚度降低率，$g(\theta)$ 是描述裂纹呼吸效应的周期函数，通常用余弦函数或者更复杂的多谐波函数来近似。裂纹每旋转一圈，刚度经历“降低-恢复”的变化，呼吸函数可以写成

$$g(\theta) = \frac{1 + \cos\theta}{2}$$

这个形式表示裂纹在受压侧时闭合（刚度恢复），在受拉侧时张开（刚度降低）。更精细的模型会用分段函数模拟“完全闭合-部分张开-完全张开”的过程。

齿轮裂纹则主要表现为齿根裂纹导致的啮合刚度变化。齿根裂纹会使该齿在啮合过程中的有效截面减少，对应的啮合刚度在啮合到裂纹齿时显著下降。在数值实现中，最直接的方法是把时变啮合刚度波形在裂纹齿对应的相位处叠加上一个刚度凹陷。设正常时变啮合刚度为 $k_m(t)$，裂纹影响区对应的刚度为 $k_{crack}$，则实际啮合刚度为

$$k_m^*(t) = k_m(t) - \Delta k \cdot h(t - t_0)$$

其中 $h(\cdot)$ 是定义在裂纹啮合相位附近的窗函数，$\Delta k$ 是裂纹导致的刚度损失量。裂纹扩展程度越大，$\Delta k$ 越大，齿面啮合时的冲击也越剧烈。

做建模时有几个很实用的参数选择。裂纹深度比 $\mu = a / h$（$a$ 为裂纹深度，$h$ 为齿高或轴半径）一般取0.1到0.5之间，太浅了响应特征不明显，太深了模型已经接近断裂，没有工程意义。刚度损失率 $\Delta k / k_0$ 在裂纹扩展初期一般是5%到15%，中期可以达到20%到40%。我建议仿真时从5%起步逐步增大，观察庞加莱截面和频谱图的变化，这样能看到故障从轻微到严重的完整演化过程。

2.4 方程的量纲一化处理

写完运动方程后不要急着求解，先做量纲一化。量纲一化的意义在于把物理参数从具体的单位系统中抽象出来，使方程组在数学上只依赖于无量纲参数组合，这样既避免了数值计算中可能出现的巨大量级差异，也让不同工况、不同尺寸的系统之间的结果具有可比性。

以齿轮转子系统为例，设参考长度为齿侧间隙的半间隙 $b$，参考时间为系统的固有周期 $T_n = 1 / \omega_n$，那么位移和时间可以分别无量纲化为

$$\bar{x} = \frac{x}{b}, \quad \tau = \omega_n t$$

代入原方程以后，原来的质量、阻尼、刚度参数就变成了无量纲阻尼比 $\zeta$、无量纲激励频率 $\Omega = \omega / \omega_n$、无量纲刚度比 $\beta = k_{mesh} / k_{shaft}$ 等参数组合。这样做还有一个额外的好处：在扫参的时候，直接扫无量纲频率比 $\Omega$ 就可以覆盖不同转速下的共振、次谐波、混沌等行为，不需要每次换算回物理量纲。

我踩过的坑是量纲一化之后经常忘了检查参数是否落在合理的物理范围。比如无量纲阻尼比 $\zeta$ 如果取得太小（小于0.01），系统会非常容易发散，数值积分步长必须压得很小，计算时间成倍增加；如果取得太大（大于0.1），混沌区域会被明显抑制，什么都看不出来。一般取0.02到0.05之间比较合适，计算效率和混沌特征能兼顾。

3. 仿真实现与庞加莱截面计算

3.1 微分方程数值求解的MATLAB实现

把运动方程写成状态空间形式是标准的做法。设状态向量为 $\mathbf{X} = [x_1, \dot{x}_1, x_2, \dot{x}_2, \ldots]^T$，那么系统可以写成一阶常微分方程组

$$\dot{\mathbf{X}} = f(\mathbf{X}, t)$$

直接写成一个函数文件。下面给一个典型的轴承-齿轮耦合系统的运动方程参考写法。

function dX = rotor_gear_dynamics(t, X, params) % 状态变量: X = [x1, dx1, y1, dy1, x2, dx2, y2, dy2] % x1,y1: 转子轴心位移; x2: 齿轮扭转位移; 加载 % 参数解包 m1 = params.m1; % 转子等效质量 c1 = params.c1; % 转子阻尼 kb = params.kb; % 轴承刚度系数 delta = params.delta; % 轴承游隙 m2 = params.m2; % 齿轮等效质量 c2 = params.c2; % 齿轮啮合阻尼 km0 = params.km0; % 齿轮平均啮合刚度 km1 = params.km1; % 时变啮合刚度一阶波动幅值 omega_m = params.omega_m; % 啮合频率 b = params.b; % 齿侧间隙半宽 F0 = params.F0; % 外载荷力幅值 omega = params.omega; % 激励频率 x1 = X(1); dx1 = X(2); y1 = X(3); dy1 = X(4); x2 = X(5); dx2 = X(6); % 轴承非线性接触力（Hertz接触 + 游隙） r = sqrt(x1^2 + y1^2); if r > delta Fb = kb * (r - delta)^1.5; else Fb = 0; end % 齿轮时变啮合刚度 km = km0 + km1 * cos(omega_m * t); % 齿侧间隙非线性函数 s = x1 - x2; if s > b g = s - b; elseif s < -b g = s + b; else g = 0; end % 运动方程 dX = zeros(6, 1); dX(1) = dx1; dX(2) = (F0 * cos(omega * t) - c1 * dx1 - Fb * x1 / r) / m1; dX(3) = dy1; dX(4) = (F0 * sin(omega * t) - c1 * dy1 - Fb * y1 / r) / m1; dX(5) = dx2; dX(6) = (km * g - c2 * dx2) / m2; end

求解的时候用ode45是最常规的选择，因为它对大多数非刚性问题都适用且精度适中。但要注意一个关键细节：默认的ode45输出时间点是不均匀的，而后面做庞加莱截面和FFT分析时通常需要等间隔采样，所以求解时最好指定固定的输出时间步。

% 参数设置 params.m1 = 10; % kg params.c1 = 20; % N.s/m params.kb = 2e6; % N/m^(3/2) params.delta = 1e-4; % m params.m2 = 5; % kg params.c2 = 30; % N.s/m params.km0 = 8e6; % N/m params.km1 = 1.2e6; % N/m params.b = 1e-4; % m params.F0 = 50; % N params.omega = 150; % rad/s % 用啮合频率作为庞加莱截面的激励周期 params.omega_m = 2 * params.omega; % 时间设置：至少要足够多个周期，保证瞬态衰减完毕 T = 2 * pi / params.omega; tspan = 0 : T/200 : 600 * T; % 初始条件（零初始即可） X0 = zeros(6, 1); % 求解 [t, X] = ode45(@(t, X) rotor_gear_dynamics(t, X, params), tspan, X0);

这里有几个坑必须提醒。第一，积分时间一定要够长，不能只看前面几十个周期，因为非线性系统在参数接近临界值时的瞬态过程非常漫长。第二，固定时间步长用的是tspan = 0 : T/200 : 600 * T的写法，200个点覆盖一个激励周期，画庞加莱截面时精度基本足够，FFT的采样率也够了。第三，ode45对带有1.5次方这种非光滑项的方程有时会报计算慢或者失败，如果遇到这种情况，优先检查是不是参数导致系统发散，再用ode15s试试，后者对付刚性方程更稳定。

3.2 庞加莱截面的代码实现

庞加莱截面的本质是：把连续时间系统的高维流形降维到离散映射上，通过观察映射点的分布结构来识别系统的运动类型。实现上不需要复杂的几何算法，最简单可靠的方法就是在时间序列上按激励周期等间隔采样。

对于周期激励系统，庞加莱截面就是在每个激励周期的整数倍时刻取样。采样点的集合 $\Sigma = { (x_1(nT), \dot{x}_1(nT), x_2(nT)) }$ 就构成了庞加莱映射。如果系统做严格的周期运动，每周期都回到同一个点，那么庞加莱截面上只有一个点；做二周期运动就有两个点；做拟周期运动时点在截面上连续分布成一条闭合曲线；做混沌运动时点呈现奇怪吸引子结构，在截面上形成分形图案。

代码实现的核心是根据激励周期确定采样索引。因为ode45输出时间点是自己指定的tspan，所以直接按下标取点即可：

% 提取稳态响应（丢弃前半段瞬态） n_skip = floor(length(t) / 2) + 1; X_steady = X(n_skip:end, :); t_steady = t(n_skip:end); % 计算庞加莱截面点：按激励周期 T 等间隔采样 N_poincare = 200; % 采样点数 poincare_x1 = zeros(N_poincare, 1); poincare_x2 = zeros(N_poincare, 1); for i = 1 : N_poincare % 找到第 i 个激励周期结束时刻对应的索引 target_t = t_steady(1) + i * T; [~, idx] = min(abs(t_steady - target_t)); poincare_x1(i) = X_steady(idx, 1); poincare_x2(i) = X_steady(idx, 5); end % 画庞加莱截面图 figure; plot(poincare_x1, poincare_x2, '.', 'MarkerSize', 6); xlabel('x_1 位移 (m)'); ylabel('x_2 位移 (m)'); title('庞加莱截面'); grid on;

取点时的注意事项非常关键。采样点必须严格对应激励周期时刻，不能随便选。如果激励是转频 $\omega$，那就每 $2\pi/\omega$ 取一次；如果激励是啮合频率 $\omega_m$，那就每 $2\pi/\omega_m$ 取一次。取错周期会导致庞加莱图完全失真。另外要丢弃足够长的瞬态段数据，我一般丢一半，如果系统瞬态特别长就丢三分之二。庞加莱截面点的数量建议选200到500个，太少看不清结构，太多计算量上去了但信息量不会增加多少。

实际使用中我发现一个非常有用的技巧：不要只画一个截面的二维图，可以把三个正交截面都画出来，分别取 $x_1$-$y_1$ 平面、$x_1$-$x_2$ 平面、$y_1$-$x_2$ 平面。这样虽然庞加莱点数量一样，但不同投影面上混沌吸引子的形态差异很大，有些藏在某一投影内的特征会在另一投影中非常明显。这个技巧在判断周期解的倍数时特别有用。

3.3 分岔图的计算与绘制

分岔图是研究非线性动力学最重要的工具之一，它展示的是当某个控制参数连续变化时系统稳态响应从周期到混沌的演化全过程。以无量纲转速比 $\Omega$ 作为控制参数时，分岔图的横轴是 $\Omega$，纵轴是庞加莱截面上的位移值。

实现上是一个嵌套循环：外层循环遍历参数，内层循环对每个参数做完整积分和庞加莱采样。

% 分岔图计算 Omega_list = 0.5 : 0.005 : 3.0; % 无量纲转速比范围 bifur_x = zeros(length(Omega_list), 50); % 每参数取50个庞加莱点 for k = 1 : length(Omega_list) Omega = Omega_list(k); params.omega = Omega * sqrt(params.km0 / params.m1); params.omega_m = 2 * params.omega; T = 2 * pi / params.omega; tspan = 0 : T/200 : 600 * T; [t, X] = ode45(@(t, X) rotor_gear_dynamics(t, X, params), tspan, X0); % 丢弃瞬态 n_skip = floor(length(t) / 3); X_s = X(n_skip:end, :); t_s = t(n_skip:end); % 每周期取一个点，共取50点 for i = 1 : 50 target_t = t_s(1) + i * T; [~, idx] = min(abs(t_s - target_t)); bifur_x(k, i) = X_s(idx, 1); end % 简单的进度提示 if mod(k, 20) == 0 fprintf('已完成 Omega = %.3f / %.3f\n', Omega, Omega_list(end)); end end % 画分岔图 figure; plot(Omega_list, bifur_x, '.', 'MarkerSize', 2); xlabel('无量纲转速比 \Omega'); ylabel('庞加莱截面位移 x_1'); title('系统分岔图'); grid on;

分岔图的绘制有几个经验性的设置。横轴参数步长不能太大，0.005算比较合适，再小计算量成倍增加，再大一些细小的分岔结构会被漏掉。每参数取50个庞加莱点足够看出分布结构了。画图时点的大小要调小，不然密集区域会糊成一片黑。

这个代码跑起来会有点慢，尤其是Omega_list有500多个参数点时。可以加一个parfor并行循环（需要Parallel Computing Toolbox），把外层循环改成parfor即可，注意循环体内的随机数或写入顺序问题要处理干净。实测下来四核机器上加速比能到三倍左右，效果还是明显的。

4. 结果解读与故障特征分析

4.1 庞加莱截面上的运动类型判别

拿到庞加莱图之后，第一步是判断系统处于什么运动状态。我把判别的口诀总结成一句话：一个点是周期一，N个点是N周期，闭合曲线是拟周期，杂乱无章是混沌。

具体来说，当庞加莱截面上只有一个孤立点时，系统做周期一运动，也就是频率成分里只有激励频率及其整数倍谐波；当有N个孤立点时，系统做N倍周期运动，说明出现了分岔，频率成分里有激励频率的1/N次谐波分量；当点在截面上连续分布、连成一条平滑闭合曲线时，系统做拟周期运动，通常意味着系统里存在两个不可公约的激励频率；当点的分布形成具有一定自相似结构的奇怪吸引子时，系统进入混沌状态，这时频谱上是连续谱叠加离散峰值，时域波形没有重复性。

庞加莱图上的点形状也有讲究。如果是周期运动但振动幅值非常大，庞加莱点可能出现一条横向的短线，这是因为ode45的插值误差在陡峭响应处被放大了。这种情况不用紧张，把采样点数量减少或者用RelTol、AbsTol收紧绝对误差就能改善。

我自己的经验是，庞加莱截面必须跟相轨迹图、FFT频谱图配合来看，单看庞加莱截面容易被误导。比如拟周期运动在庞加莱截面上是闭合曲线，但有些参数下闭合曲线退化得很窄，看起来像一簇点，这时看频谱如果存在两个不可公约的频率峰值就能确认是拟周期。再比如混沌运动在某些投影面上可能看起来像一堆点，但放到另一个投影面就能看出奇怪吸引子的折叠结构。三张图配合起来判断，准确率要高得多。

4.2 裂纹故障的动力学特征与识别

引入裂纹故障之后，系统响应的特征变化非常明显，这也是裂纹故障诊断的理论基础。我把裂轴和裂齿两种情况分别说。

转轴裂纹导致的刚度周期性变化会在频谱中产生两个典型特征：一是在转频的1/2倍频处出现次谐波峰，这是因为裂纹呼吸函数的主要频率成分是转频的两倍，与系统某阶固有频率发生参数共振时会出现亚临界共振；二是转频的高次谐波幅值会明显增大，尤其以2倍频和3倍频的幅值增长最快，这与裂纹引起的刚度不对称直接相关。

齿轮齿根裂纹的动力学特征跟转轴裂纹不太一样，主要表现是啮合频率及其边频带的改变。正常齿轮的啮合频率附近边带幅值较小，而裂纹齿投入啮合时会产生一个幅值明显的冲击，这个冲击在频谱上表现为啮合频率两侧的边频带幅值增大，且边频带间隔等于转轴转频。随着裂纹扩展，啮合频率的高次谐波特别是2倍啮合频率、3倍啮合频率处的边带幅值急剧增大，同时系统进入混沌运动的参数区间变宽，庞加莱截面上从周期点变成奇怪吸引子的转速范围明显扩大。

从工程诊断的角度，我总结了一个相对实用的判断方法：看庞加莱截面上点的数量随转速的变化趋势。正常齿轮系统在较宽的转速范围内保持周期一或周期二运动，庞加莱点数量少且固定；裂纹故障系统则会在更多转速段内出现高倍周期的点群和混沌吸引子，点的数量随转速变化更加频繁。这种“周期结构随转速剧烈变化”本身就是故障的指纹特征。

4.3 从理论仿真到工程诊断的映射

仿真分析的最终目的是服务诊断，所以做完理论分析后，一定要把仿真特征和实际测试特征对应起来。这里要特别强调一点：仿真中的庞加莱截面是理想化的，实际振动信号里的“庞加莱截面”对应的是同步采样下的轴心轨迹或齿轮振动信号的等间隔采样。

具体操作中，如果现场有键相传感器（每转一个脉冲），就以外触发信号为基准采集振动信号，这样采集到的每个数据块都正好对应转轴一周，把同一相位位置的振动幅值提取出来，按顺序排列，就相当于实现了工程版庞加莱截面。用这个思路，即使不做复杂的信号处理，也能在现场识别出周期运动和混沌运动：混沌工况下同一相位处的振动值每次都不一样，而且没有规律性，这跟仿真中庞加莱点的分布特征是一一对应的。

还有一个非常实用的映射关系：仿真中裂纹故障导致的分岔提前现象，在实验台架上可以验证为振动频谱中边频带幅值的提前增大。举个例子，某型齿轮箱仿真中健康状态在转速比2.3附近才进入拟周期运动，而裂纹深度比0.3时，在转速比1.8附近就出现了明显的次谐波成分，这个转速提前量可以作为裂纹早期预警的特征指标。做诊断算法时，不需要完全复现理论混沌，只需要提取提前出现的次谐波幅值变化率就可以做出趋势预警。

5. 常见问题排查与实操心得

5.1 数值求解中的典型问题与对策

这里我要专门强调一个新手最容易犯的错误：把分岔图横轴的突然跳变当成程序bug，其实那恰恰是分岔点。分岔点处系统结构发生突变，比如从周期一突然跳到混沌，这是非线性动力学里非常有价值的信息，千万不要“修”掉它。

另一个常见问题是参数没有量纲一化导致的数值病态。如果直接使用实际物理参数，比如刚度是 $10^7$ 量级，位移是 $10^{-4}$ 量级，两者的乘积和阻尼项、惯性项在数值上会跨越十几个数量级，ode45的误差控制会被极端量级搅乱，导致步长自动调节异常。这个问题我在第一次做齿轮动力学仿真时遇到过，花了整整一天排查参数，最后发现只是没做量纲一化，做完量纲一化后同样的方程跑起来飞快。

5.2 参数选择的实用经验

非线性动力学仿真的参数选择要遵循“先简单后复杂，先稳定后混沌”的原则。具体来说有几种实际经验值得分享。

第一，初始参数选在系统稳定周期运动区间，先跑通流程再改参数。比如无量纲转速比从0.8开始，通常周期一运动比较稳定，庞加莱图就是一个点，所有代码全部跑通确认无误后，再逐步增大转速比扫描分岔区间。直接设置一个可能处于混沌的参数，结果图出来很漂亮但很难校对对不对。

第二，激励幅值从零开始逐渐增大。非线性系统的幅频响应有所谓跳跃现象，激励幅值不同，系统的分岔结构完全不同。从零幅值开始等于先算线性系统，容易验证线性部分的正确性，再慢慢增加非线性强度观察变化。

第三，阻尼比保持在0.01以上。非线性系统在低阻尼下数值计算极其吃力，庞加莱截面上的点会出现大范围的伪散布，因为数值误差沿混沌轨道的指数发散被放大。阻尼比0.02到0.05的区间内，系统既能展示丰富的非线性现象又不会把数值计算逼到死胡同。

5.3 从零开始搭建仿真的避坑指南

给刚开始做这类课题的朋友一个可以照抄的路径规划，这条路我自己走过，也看很多学生走过，踩坑率最低。

第一步，先做纯线性系统的验证。把轴承刚度当常数，齿轮啮合刚度当常数，齿侧间隙设为零，跑通ode45，和解析解对比，确认程序基本框架没有bug。第二步，把齿轮啮合刚度改为时变项，观察是否出现参数激励共振，此时系统应该出现周期二或周期三的子谐波共振，这和文献结果可以对上。第三步，引入齿侧间隙分段函数，系统会出现脱啮-冲击现象，庞加莱截面上开始出现非光滑映射的特征。第四步，加入裂纹刚度变化模型，对比健康系统和故障系统的分岔图差异。第五步，切换控制参数做完整的参数扫描，得到分岔图后确定混沌区域，最后在混沌区域里画庞加莱截面、相图、频谱图和Lyapunov指数，完成整套非线性特征分析。

每一步都确认无误再进入下一步，不要跳步。我见过太多同学一上来就搞完整模型，结果庞加莱图上出现了某种结构，但根本说不清是间隙引起的还是裂纹引起的还是数值误差造成的。分步验证的好处在于，每一个非线性因素从无到有引入时，你能直接看到它对应的影响是什么，这种“因果对应”的经验积累起来之后，后面的故障特征分析会非常顺手。

还有一个细节想提醒大家：保存仿真数据时，把时间序列、庞加莱点、参数值、程序版本全部存在同一个文件名里。非线性系统的结果对参数极其敏感，你三个月后回来看一个庞加莱图，如果不知道当初用的具体参数值，这张图等于白画。我自己吃过这个亏，后来写了个简单的数据管理脚本，每次仿真自动将参数保存为同名mat文件，再也不会出现“这张图是哪个参数跑出来的”这种尴尬问题。

5.4 代码性能优化的几个实操手段

非线性动力学扫描参数次数多，性能优化非常值得投入时间。第一个手段是适当地降低输出数据量。前面例子中tspan用T/200的步长，如果参数扫描范围很大，可以改成T/50先粗扫一遍，找到分岔结构的大致区域后，再对感兴趣的区域用T/200细扫。粗扫的庞加莱图可能有点毛糙，但分岔趋势完全能看清，计算量能省四倍。

第二个手段是用solver的OutputFcn来提前终止判定。当系统已经明显进入混沌时，继续积分只是浪费时间，可以设置一个判断条件：如果相邻庞加莱点之间的最大距离在连续几百个周期内都超过某个阈值，就认为已经进入稳态混沌，可以提前跳出循环。这个优化在扫描整个大参数范围时效果显著，有些参数点本来需要跑完600个周期，实际上150个周期就可以判断出来了。

第三个手段是在画分岔图时先用稀疏矩阵存储，最后一次性画图。直接把一个500乘2000的矩阵塞进plot不会出问题，但如果中间还有后续数据处理，建议只在最后保存分岔图的横坐标数组和对应的庞加莱点矩阵，减少内存占用。我扫过最多的一次是Omega_list有2000个点，每点取100个庞加莱点，用上述手段在普通笔记本上跑20分钟完成，内存占用不到10MB，这个量级对大多数研究场景都够用了。

做完这套流程之后，我对非线性动力学仿真最大的体会是：关键不在工具，而在模型和分析逻辑。MATLAB的代码几天就能写完，真正花时间的是理解每个非线性项对系统的贡献。用最简单的模型验证最核心的逻辑，再用逐步增加复杂度的方法逼近真实系统。分岔图、庞加莱截面和频谱图是三个互相印证的工具，一个结果有疑问，永远用另外两个去交叉验证。