考虑时变压力角和时变齿侧间隙的直齿轮六自由度平移-扭转耦合非线性动力学程序，包括时域图，相图，FFT图，庞加莱图，分岔图。 要想学好齿轮动力学，需要有扎实的齿轮动力学理论和非线性动振动理论。 齿轮啮合刚度建模是齿轮动力学求解的第一步。
齿轮动力学这玩意儿，真不是背几个公式就能上手的。前段时间为了搞定一个六自由度平移-扭转耦合模型，硬是折腾了两个月才摸清门道。今天咱们就聊点实在的——如何从零开始构建这个带时变压力角和齿侧间隙的非线性动力学程序。
先说说齿轮啮合刚度的建模吧，这玩意儿简直就是整个系统的地基。传统方法用抛物线函数近似刚度变化，但在实际工况下，接触线长度的动态变化会让刚度曲线出现明显的"台阶"特征。我一般习惯用分段函数来处理这种非线性：
def time_varying_stiffness(t, omega): phase = omega * t % (2*np.pi) if phase < np.pi/3: return 1.2e8 + 5e7*np.sin(3*phase) elif np.pi/3 <= phase < 2*np.pi/3: return 8e7 * (1 + 0.5*np.cos(3*phase)) else: return 9e7 * np.exp(-0.1*(phase - 2*np.pi/3))这个模型考虑了接触线长度突变导致的刚度跳跃，比单一的正弦函数更贴近实测数据。注意指数项用来模拟冲击衰减，这对后续的振动谐波分析很关键。
动力学方程部分要特别注意惯性耦合项。六自由度模型包含两个齿轮的三个平移和三个扭转自由度，微分方程组的刚性特征非常明显。用SymPy生成符号方程能省不少事：
from sympy import symbols, Matrix # 定义自由度向量 x1, y1, z1, theta1, x2, y2, z2, theta2 = symbols('x1 y1 z1 theta1 x2 y2 z2 theta2') q = Matrix([x1, y1, z1, theta1, x2, y2, z2, theta2]) # 构造质量矩阵（以简化版为例） M = Matrix.diag([m1, m1, m1, I1, m2, m2, m2, I2]) # 刚度矩阵要考虑时变项 K = Matrix([ [kxx, kxy, 0, kxtheta, -kxx, -kxy, 0, -kxtheta], [...], # 完整矩阵需补充其他元素 ])实际求解时发现，当齿侧间隙超过50μm时，常规的Runge-Kutta法会发散。后来改用变步长的Radau IIA算法，配合事件检测处理间隙碰撞，稳定性提升明显：
options = odeset('RelTol',1e-8,'AbsTol',1e-10,'Events',@backlash_events); [t,y,te,ye,ie] = ode15s(@gear_odefun, [0 0.5], y0, options);分岔图的绘制要注意参数扫描策略。固定转速下逐渐改变阻尼比，每次仿真需要跑够200个啮合周期后再采集数据。这个循环在Python里可以这样实现：
bifurcation = [] for c_ratio in np.linspace(0.02, 0.2, 100): sol = solve_ivp(fun, [0, 200*T], y0, args=(c_ratio,)) # 取最后10个周期的峰值 peaks, _ = find_peaks(sol.y[0][-10000:]) bifurcation.extend([(c_ratio, sol.y[0][i]) for i in peaks])庞加莱截面最能暴露系统的混沌特性。采样点选在驱动轴每转整周期时刻，能清晰看到吸引子结构。有个坑要注意：必须用精确的事件检测，直接按时间间隔采样会导致截面失真。
FFT分析建议用Welch方法平滑频谱。某次测试中在3倍啮合频率处发现了明显的边频带，排查发现是时变刚度中的三次谐波与齿频调制产生交互作用。这个现象用常规的线性模型根本捕捉不到。
调试过程中最抓狂的是参数敏感性问题。啮合阻尼系数差0.001，分岔图就能从周期3直接跳变到混沌状态。后来引入Sobol序列进行全局参数扫描，才摸清了各参数的敏感区域。现在回头看，那些熬夜调参的日子，都是理解非线性耦合本质的必经之路啊。
