matlab:行星齿轮动力学，集中质量参数模型，基于势能法求解齿轮时变啮合刚度，行星齿轮系统动态响应，matlab源码。

行星齿轮系统这玩意儿玩起来是真的上头，最近在搞集中质量参数模型的时候，发现用势能法计算时变啮合刚度比传统方法有意思多了。咱们今天就手撕一段MATLAB代码，看看怎么把齿轮啮合刚度算得明明白白，顺便整个行星轮系的动态响应分析。

先整几个齿轮参数热热身：

% 行星轮系基础参数
N_p = 3;           % 行星轮数量
m = 4;             % 模数(mm)
Z_s = 28;          % 太阳轮齿数
Z_p = 36;          % 行星轮齿数
pressure_angle = 20*pi/180; % 压力角

这里藏着个关键点——行星轮相位差得算准了，不然后面动态响应会抽风。每个行星轮的初始相位得按360/N_p均匀分布，这个在微分方程里会直接影响激励项的相位关系。

说到势能法求啮合刚度，核心就是计算齿轮副的势能变化。咱们把啮合线刚度拆成弯曲刚度、剪切刚度和接触刚度三部分：

function k = meshing_stiffness(E, h, L, mu, theta)
    % 弯曲刚度分量
    k_bending = E*h^3/(12*L*(1 - mu^2)*sin(theta)^2);
    
    % 剪切刚度分量
    k_shear = 5*E*h/(12*L*(1 + mu)*sin(theta)*cos(theta));
    
    % 赫兹接触刚度
    k_contact = pi*E*L/(4*(1 - mu^2)*log(4*h/(3*L)));
    
    k = 1/(1/k_bending + 1/k_shear + 1/k_contact); % 串联刚度
end

这里有个坑要注意：当齿轮进入双齿啮合区时，总刚度其实是两个啮合点刚度的并联。所以在时间维度上得做个循环判断，实时叠加刚度值。这个判断逻辑直接关系到时变刚度的计算精度。

动力学方程这块，集中质量模型得考虑太阳轮、行星轮和齿圈的平移-扭转耦合。举个微分方程组的例子：

function dx = dynamics(t, x, par)
    % 解包参数
    k_mesh = par.k_mesh(t); % 时变啮合刚度
    c = par.c;              % 阻尼系数
    
    % 系统自由度分配
    dx = zeros(12,1);
    % 太阳轮平动方程
    dx(1:2) = x(3:4);
    dx(3:4) = (k_mesh*(x(7:8) - x(1:2)) - c*(x(9:10) - x(3:4)))/par.m_s;
    
    % 行星轮转动方程
    for p = 1:N_p
        idx = 4 + 2*p;
        dx(idx+1:idx+2) = ... % 具体方程省略
    end
end

这里用了硬核的矢量写法，行星轮的位置索引得特别注意。建议用结构体数组来管理各部件参数，比直接操作索引清爽得多。

求解的时候推荐用ode45的变步长算法，但记得在options里把RelTol调到1e-6以下，不然高频振动成分会被吃掉。跑完仿真后，用短时傅里叶变换看时频谱特别带劲：

[~, sol] = ode45(@(t,x)dynamics(t,x,par), t_span, x0);
fs = 1/(t(2)-t(1));
spectrogram(sol(:,1), 256, 250, 256, fs, 'yaxis')

当看到频谱图里出现明显的边频带，八成是行星轮通过频率的调制现象。这时候回看刚度计算模块，检查是否准确捕捉了啮合相位的变化。

完整代码里其实还藏着几个骚操作：用查表法预生成刚度曲线减少实时计算量，用事件函数捕捉齿面分离现象，还有用移动矩阵法处理行星轮相位耦合。这些技巧能让计算速度提升三倍不止，不过代码量就有点大了，这里先挖个坑，下回接着唠。