行星齿轮组动力学建模与ode45求解实践

简介：面向机械工程、车辆工程与齿轮动力学方向的学习者和研究者，这份MATLAB代码资源聚焦行星齿轮组微分方程组的建模与数值求解。资源将太阳轮、行星轮、行星架和内齿圈构成的多体动力学问题转化为ODE初值问题，并通过ode45进行求解；代码中设置了0.001～1000s的求解区间与时间步长，同时将齿轮系统刚度矩阵以ky数组形式传入，便于考察刚度变化对动态响应的影响。压缩包内共3个M函数脚本，大小仅4KB，分别负责刚度矩阵处理、齿轮振动特性分析与微分方程组主求解，结构简洁、便于修改扩展；脚本中保留了完整的参数定义与调用关系，可作为教学演示或科研入门模板。已有898人学习借鉴。借助该资源，读者可快速掌握基于ode45的行星齿轮组动态响应仿真流程，理解步长选择、刚度参数与振动结果之间的关联，也可将其中代码框架迁移至其他齿轮传动系统，用于参数优化或故障特征分析。

1. 行星齿轮组计算的微分方程组，为什么用 ode45 即可求解

行星齿轮组的动态计算常被高估：很多人以为必须上多体软件和专用求解器。实际把广义坐标、啮合刚度、误差激励写成微分方程组之后，一个 ode45 就能覆盖绝大多数工况——成败在建模，不在求解器。

风电齿轮箱、电驱减速器、机器人关节减速器的振动校核，最常见做法是集中参数模型：构件按刚体，啮合副用变刚度与阻尼表达，得到二阶微分方程组，降阶后交给 ode45，几十个啮合周期分钟级算完。

下面按一线做法讲：行星齿轮组方程怎么写对、ode45 最小可运行代码、参数与排错。适合刚接触齿轮动力学的工程师，也适合已有方程但在容差、步长上反复踩坑的人。

2. 行星齿轮组建模：从自由度、啮合刚度到微分方程组的矩阵装配

多数人卡住的地方不是求解，而是方程组本身写不写得出来。行星齿轮组的方程不建议一条条手推，我一般按“啮合单元装配”来做：每个太阳轮–行星轮、内齿圈–行星轮的啮合副，都是一组时变刚度加阻尼加误差激励的弹簧阻尼单元；装配规则固定，改齿数只改矩阵元素，不改程序结构。

2.1 自由度与广义坐标：先建扭转模型，再考虑平移

内齿圈固定、单级行星传动，最常用的起点是扭转模型：只保留旋转自由度，广义坐标取q = [θ_s, θ_c, θ_p1, …, θ_pN]，即太阳轮转角、行星架转角、N 个行星轮自转转角，总自由度n = N + 2。先做扭转模型的原因是行星传动最核心的动态现象——齿频振动、动态啮合力、行星轮载荷分配——在这个模型里都已出现，而调试成本比带平移的模型低一个量级。

平移自由度后面再加：要算行星轮轴承受力、太阳轮浮动均载、行星架横向振动时，每个构件补 x、y 两个位移，状态数从 2n 涨到 6n 左右，对 ode45 仍是小规模问题，真正贵的是轴承刚度与侧隙参数的标定，不是求解本身。无论哪种模型，几何上必须先满足z_r = z_s + 2z_p；差一个齿，啮合关系不成立，后面所有矩阵元素全是错的。

2.2 啮合刚度、阻尼与误差激励：时变是这组方程的核心

每个啮合副的刚度随轮齿进入、退出啮合周期变化，工程上常用一阶傅里叶形式近似：

k(t) = k_m [1 + k_a cos(ω_m t + φ)]

k_m 是平均啮合刚度，k_a 是波动幅值比，φ 是初始相位。啮合角频率 ω_m 是整组方程的标尺：啮合频率f_m = z_s(n_s − n_c)/60，内齿圈固定时n_c = n_s z_s/(z_s + z_r)，因此f_m = z_s n_s z_r / [60(z_s + z_r)]。后面所有步长、FFT 横轴、稳态时长判断都要对齐到它。

相位 φ 特别容易漏：N 个行星轮均布时，太阳轮–各行星轮啮合相位差取2π z_s (i−1)/N对 2π 取模；不写这个相位，行星轮载荷分配算出来明显偏掉。误差激励同样按齿频给：e(t) = e_a sin(ω_m t + φ_e)，代表基节误差、齿形误差的累计效应，幅值按齿轮精度等级取若干微米。阻尼按啮合等效质量折算，见第 4 章。

2.3 按啮合单元装配：矩阵形式与投影向量

装配分成三步：写投影向量、算单元刚度、叠加全局矩阵。以太阳轮–行星轮 i 为例，啮合线相对变形定义为：

δ_spi = r_s(θ_s − θ_c) − r_p(θ_pi − θ_c)

这就是一个 1×(N+2) 行向量与 q 的内积，记为e_spi·q。对应代码：

% 太阳轮-第 i 个行星轮啮合的投影向量，n = Np + 2 e = zeros(1, p.Np + 2); e(1) = p.rs; % 太阳轮转角系数 e(2) = -(p.rs-p.rp); % 行星架转角系数 e(2+i) = -p.rp; % 第 i 个行星轮转角系数 p.Lsp(i,:) = e;

内齿圈–行星轮同理，只是符号习惯不同。注意：不同论文对“压缩为正”的约定差一个正负号，统一后全程序沿用，不要中途换方向。全部装配完后，运动方程写成：

M q̈ + C_b q̇ = T(t) − K_m(t) q + F_nl(q, q̇, t)

M 是惯量对角阵，C_b 是支撑阻尼，K_m(t) 由所有啮合单元叠加，F_nl 放侧隙和误差引起的非线性力。这就是 ode45 要积分的对象。扭转模型的关键参数基准如下：

齿数互质那条经验很有用：z_s 与 z_p 存在公约数时，同一批轮齿周期性重复啮合，误差激励的相位结构改变，模拟结果和实际齿轮箱的错峰特性对不上。

提示：投影向量、误差相位、啮合频率必须来自同一个 p 结构，RHS 函数、事件函数、后处理三处代码共用，否则会出现“能积分、但结果对不上”的怪异现象。

3. 把行星齿轮组微分方程组写成状态方程，用 ode45 跑通最小算例

3.1 二阶方程组降阶：状态变量取“位置 + 速度”

ode45 只吃一阶显式方程组y' = f(t, y)。行星齿轮组运动方程是二阶的，统一做法是翻倍状态：y = [q; q̇]，于是dy/dt = [q̇; q̈]。q̈ 由运动方程显式解出：质量矩阵 M 是常数对角阵（扭转模型）或对称正定阵，求逆一次缓存，不要在 RHS 里反复 inv。

3.2 求解函数与主脚本：我给的最小可运行结构

3.2.1 右侧函数 planetary_rhs 与接触力子函数

RHS 函数负责在每时刻算广义加速度：

function dydt = planetary_rhs(t, y, p) % 状态 y = [q; qd]，q 是 N+2 个转角，qd 是角速度 n = p.ndof; q = y(1:n); qd = y(n+1:end); % 所有啮合副的变形与变形速度（含误差激励） dsp = zeros(p.Np, 1); drp = zeros(p.Np, 1); dds = zeros(p.Np, 1); ddr = zeros(p.Np, 1); for i = 1:p.Np dsp(i) = p.Lsp(i,:)*q - p.es(i)*sin(p.wm*t + p.ps(i)); drp(i) = p.Lrp(i,:)*q - p.er(i)*sin(p.wm*t + p.pr(i)); dds(i) = p.Lsp(i,:)*qd - p.es(i)*p.wm*cos(p.wm*t + p.ps(i)); ddr(i) = p.Lrp(i,:)*qd - p.er(i)*p.wm*cos(p.wm*t + p.pr(i)); end % 侧隙接触力：|d| < bg 时脱离，力为 0 Fm = zeros(2*p.Np, 1); for j = 1:p.Np Fm(2*j-1) = mesh_force(dsp(j), dds(j), p.km(2*j-1), p.cm(2*j-1), p.bg); Fm(2*j) = mesh_force(drp(j), ddr(j), p.km(2*j), p.cm(2*j), p.bg); end % 装配到广义坐标并解出加速度 Fq = zeros(n, 1); for i = 1:p.Np Fq = Fq + p.Lsp(i,:)'*Fm(2*i-1) + p.Lrp(i,:)'*Fm(2*i); end qdd = p.M \ (p.T(t) - p.Cb*qd - Fq); dydt = [qd; qdd]; end

function F = mesh_force(d, dd, k, c, bg) % d 是啮合线相对变形，正方向为齿面压缩 if d > bg F = k*(d - bg) + c*dd; elseif d < -bg F = k*(d + bg) + c*dd; else F = 0; end end

两个函数构成完整右侧：planetary_rhs 里每个啮合副的变形由投影向量L·q完成，删改某一行轮齿误差只动 es/er 的对应元素；增加行星轮数量只需扩展 Lsp/Lrp 和 Jp，程序结构完全不变。

3.2.2 主脚本：ode45 调用与 deval 取数

p = define_gear_params(); % 几何、刚度、阻尼参数，见第 4 章 p.ndof = p.Np + 2; p.M = diag([p.Js, p.Jc, p.Jp]); % 惯量对角阵；Jc 含行星轮公转惯量 p.T = @(t) input_torque(t, p); % 返回 n×1 广义力：太阳轮输入+行星架负载 y0 = zeros(2*p.ndof, 1); % 从静止启动 opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'MaxStep', 1/p.fm/50); sol = ode45(@(t,y) planetary_rhs(t,y,p), [0 0.5], y0, opts); tq = linspace(0, 0.5, 5000); Y = deval(sol, tq); % 均匀网格取数，避免线性插值 omega_s = Y(p.ndof+1, :); % 太阳轮角速度

说明：M 用左除\而不是inv(M)*F，对角阵左除是 O(n)，扩展到含平移自由度的模型更稳。input_torque是带斜坡的广义力函数，避免零时刻突加扭矩把初始瞬态拉大。deval用 sol 的连续插值在任意网格取数，比[t,y]=ode45(...)之后用 interp1 重采样精度高得多，后处理做 FFT 时这里不能省。

3.3 ode45 选项：RelTol、AbsTol、MaxStep 怎么配合

ode45 默认容差在齿轮动力学里通常不够，原因有两个：转角状态幅值只有 10⁻⁴~10⁻³ rad 量级，默认 AbsTol=1e-6 会把误差预算全花在噪声上；啮合频率到千赫兹级时，默认 InitialStep 可能跨过好几个接触切换点。我常用的设置如下：

AbsTol 传向量是混合单位模型的关键：平移自由度以米计、旋转以弧度计，数量级差上万倍，单标量容差会让所有状态共用一把尺子。向量写法AbsTol = [1e-9*ones(n,1); 1e-6*ones(n,1)]，位置和速度分开给。MaxStep 与 f_m 绑定是全链路是否可信的分水岭，第 5 章会看到它如何影响结果。

4. 行星齿轮组仿真参数设定：啮合刚度、阻尼比、转速与侧隙

4.1 啮合刚度均值与波动幅值：先按平均刚度算，再叠加波动

k_m 最可靠的来源是 ISO 6336 或齿轮有限元接触分析，但前期校核经常没有这两个条件。工程做法是：钢制直齿轮、模数 2~3 mm、齿宽 15~25 mm 的经验区间，取沿啮合线的平均刚度 0.8×10⁸~2.5×10⁸ N/m；波动幅值比 k_a/k_m 由重合度 ε_α 决定，ε_α≈1.3 时约 0.3，ε_α 到 1.6 就降到 0.2 以下。教材里的齿轮刚度是分段的（单齿对/双齿对交替），一阶傅里叶近似丢掉拐点信息，但保住了齿频成分，对 5% 精度的动态因子校核足够。

调试顺序务必按“从线性到非线性”走：先把 k_a、误差激励、侧隙全置 0，确认系统是线性时不变自由振动，固有频率能对上解析解；再逐项打开。三项全开时一旦发散，你分不清是刚度、误差还是侧隙引入的问题。

4.2 阻尼比与等效质量：啮合阻尼和支撑阻尼分开取

啮合阻尼不能用 M 的倍数硬带，正确做法是按等效质量折算：m_s = J_s/r_s²，m_p = J_p/r_p²，m_eq = 1/(1/m_s + 1/m_p)，代入c_m = 2ζ√(k_m m_eq)。得到的阻尼常数单位是 N·s/m，直接乘到变形速度上。ζ 的选取：精密磨齿、油润滑取 0.02~0.05，普通滚齿或齿面粗糙度大取 0.05~0.1。支撑阻尼 C_b 用瑞利阻尼或按临界阻尼的 0.5%~2% 给，先给小的——阻尼给大了，行星轮载荷分配的差异会被抹平，模型反而失真。

4.3 转速、扭矩与侧隙：三个最容易被设错的外部参数

转速只有一个公式：f_m = z_s n_s z_r / [60(z_s + z_r)]，n_s 是太阳轮输入转速（r/min）。由它算出的 f_m 必须与第 3 章 MaxStep、第 6 章 FFT 期望峰值完全一致；f_m 对不上，后面所有频域分析都会错位。扭矩方面，太阳轮输入与行星架输出满足T_c = T_s (z_s + z_r)/z_s，两侧之比就是传动比之反比；输入给斜坡T(t) = T_0·min(1, t/0.02)，20 ms 斜坡覆盖两三个啮合周期即可。侧隙方面，沿啮合线的半间隙 b_g 取 20~80 μm，对应制造精度与润滑间隙；b_g 过大时脱啮段变长，动态因子明显上升，这就是齿轮箱“有间隙就有冲击”的数值表现。

5. ode45 求解行星齿轮组微分方程组的排错：刚性判断、容差与事件函数

5.1 用特征值判断，不做刚性猜测

“我的方程组要不要换 ode15s”是问得最多的问题。判断标准不是感觉，是特征值跨度。取平均刚度 K_mean（k_a 置 0 后的 K_m），算广义特征值：

[V, D] = eig(p.M \ p.Kmean); f_n = sqrt(abs(real(diag(D))))/2/pi; fprintf('固有频率范围: %.1f ~ %.1f Hz\n', min(f_n), max(f_n)); fprintf('最大/最小固有频率比: %.0f\n', max(f_n)/min(f_n));

特征值跨度在 10⁴ 以上才需要考虑 ode15s。行星齿轮组扭转模型里这个比值通常在 10²~10³，属于中等刚度问题，ode45 的变步长在几十个啮合周期内效率高于 ode15s。加了行星轮平移自由度后，轴承刚度软、啮合刚度硬，跨度可能升到 10⁵，那时才需要认真考虑换解算器。标题说“使用 ode45 即可求解”，适用边界就在这里：验证性计算、几十到几百个啮合周期，完全够。

5.2 解发散的排查顺序：先查装配，再查容差

发散的问题九成不在求解器。按三步排查：

1. 关掉波动、误差、侧隙，从静止启动。能量不应随时间增长；角速度单调发散时，先查投影向量 Lsp/Lrp 的符号，负刚度会把能量持续注入系统。
2. 打开刚度波动，幅值从 0.05 逐步加到位。每一步都确认位移时程仍有界。
3. 打开侧隙后出现高频毛刺是正常的，那是接触切换点附近局部步长收缩。若出现 “Unable to meet integration tolerances” 告警，优先放宽 AbsTol 而不是 RelTol。

第三步的告警很典型：理想侧隙是分段线性函数，切换点不可导，变步长算法被迫把步长压到极小去满足容差。此时把 AbsTol 从 1e-8 放到 1e-6，或在切换点附近做 5 μm 宽的三次多项式光滑过渡，步数立刻降一个量级。

5.3 用事件函数记录脱啮时刻

不需要记录时，接触切换点交给 ode45 自行处理，对精度没有影响。要统计每个行星轮的脱啮比或做啮合冲击分析时，用 Events 把切换时刻精确抓出来：

function [val, isterm, dir] = mesh_events(t, y, p) q = y(1:p.ndof); val = zeros(2*p.Np, 1); for i = 1:p.Np dsp = p.Lsp(i,:)*q - p.es(i)*sin(p.wm*t + p.ps(i)); val(2*i-1) = dsp - p.bg; % 进入接触 val(2*i) = dsp + p.bg; % 脱离接触 end isterm = zeros(2*p.Np, 1); % 只记录，不停止积分 dir = zeros(2*p.Np, 1); end

主脚本把它放进 odeset，并改用带事件输出的调用形式：

opts = odeset(opts, 'Events', @(t,y) mesh_events(t,y,p)); [t, y, te, ye, ie] = ode45(@(t,y) planetary_rhs(t,y,p), [0 0.5], y0, opts);

te 是事件时刻，ie 是事件编号，两者配合能画出“第几个行星轮在什么时刻脱啮”。isterm 全为 0 是关键语义：ode45 只记录事件、遇到事件不终止积分，这正是统计脱啮比需要的。

5.4 当 ode45 越算越慢：步数与性能对账

能跑但很慢的情况，按顺序对账：MaxStep 是否绑定了 f_m；AbsTol 是不是设了小到 1e-12；侧隙段有没有做光滑过渡。三项都查过还慢，就把仿真时长按啮合周期数定义，而不是按秒数——稳态响应通常 20 个啮合周期后不再变化，多跑全是浪费。这时看sol.stats.nsteps很直观：一个啮合周期超过 200 步，多数是容差过紧，不是 MaxStep 上限过严。

6. 从 ode45 的解里提取动态啮合力，并核对行星齿轮组传动比

6.1 动态啮合力与动态因子：一行循环恢复所有啮合副

ode45 返回的是位移和速度历程，啮合力不在状态里，需要恢复。按与 RHS 完全相同的变形定义，对每个采样时刻重算：

N = length(tq); Fsp = zeros(p.Np, N); for k = 1:N q = Y(1:p.ndof, k); qd = Y(p.ndof+1:end, k); for i = 1:p.Np dsp = p.Lsp(i,:)*q - p.es(i)*sin(p.wm*tq(k)+p.ps(i)); dds = p.Lsp(i,:)*qd - p.es(i)*p.wm*cos(p.wm*tq(k)+p.ps(i)); if dsp > p.bg Fsp(i,k) = p.km(2*i-1)*(dsp-p.bg) + p.cm(2*i-1)*dds; end end end Kv = max(Fsp, [], 2) ./ mean(Fsp, 2); % 动态因子与行星轮载荷分配

注意取均值时只用稳态段（例如后 80% 时间），避开启动瞬态。Kv 是设计里直接可用的量，对应 GB/T 3480 里的 KV 系数的工程修型。接着对比 N 个行星轮的平均力，差值超过 10% 说明均载有问题——扭转模型里这通常源于相位赋值错误，回头检查 2.2 节的相位公式即可。

6.2 传动比核验与频谱核对：判断链路是否可信

仿真是自洽的，不代表是对的。两个不依赖网格的核对手段：

% 核对 1：传动比（内齿圈固定，太阳轮输入、行星架输出） ws = Y(p.ndof+1, :); % 太阳轮角速度 wc = Y(p.ndof+2, :); % 行星架角速度 ratio_sim = mean(ws(500:end)) / mean(wc(500:end)); ratio_ana = (p.zs + p.zr) / p.zs; fprintf('仿真传动比 %.4f / 解析值 %.4f\n', ratio_sim, ratio_ana); % 核对 2：动态啮合力频谱，峰值应落在 f_m 及整数倍处 fs = 1/(tq(2)-tq(1)); S = fft(Fsp(1,:) - mean(Fsp(1,:))); f = (0:length(S)-1)/length(S)*fs; [~, pk] = max(abs(S(1:floor(length(S)/2)))); fprintf('FFT 峰值频率 %.1f Hz / 啮合频率 %.1f Hz\n', f(pk), p.fm);

传动比比对取后 80% 的数据，避开瞬态，误差应在 0.5% 以内；偏差更大时先复查投影向量和 z_r = z_s + 2z_p 的几何约束。FFT 峰值与 f_m 偏差超过 2%，优先怀疑采样网格与 f_m 的对齐，而不是模型本身。两关都过，再做一次固有频率交叉验证：关掉激励给初始位移，自由衰减响应的 FFT 峰值与 eig 结果差 2% 以内，这套“行星齿轮组微分方程组 + ode45”的链路就可以拿去改齿数、改转速、换工况了。

本文还有配套的精品资源，点击获取