简介本资源是一份面向机械工程、动力学建模与传动系统优化领域的专业MATLAB脚本聚焦行星齿轮传动系统的扭振特性分析适用于高校研究生、传动系统设计工程师及振动控制方向研究人员。它通过建立行星齿轮机构含太阳轮、行星轮、行星架与齿圈的集中参数动力学模型求解其扭转振动模态辅助识别共振频率、评估系统稳定性并支撑减振结构优化。压缩包仅含1个核心文件NewMyModel.m3KB为可直接运行的MATLAB函数脚本内含动力学方程线性化、质量/刚度矩阵构建、特征值求解及模态频率输出等完整流程代码结构清晰、注释充分便于理解行星齿轮扭振建模逻辑与数值实现方法。目前已有367人学习下载适合需要快速掌握行星齿轮传动动态建模与模态分析实操能力的进阶学习者。 如果你也在做齿轮传动系统的动力学分析或者正在被行星齿轮的扭振问题折磨那这篇东西应该对你有用。前几天我把手头代号为 NewMyModel 的模型又重新梳理了一遍这个模型的核心就是行星齿轮传动系统的扭振动力学建模与模态分析。说实话这类项目说难不难说简单也不简单——难的是把物理过程转成正确的动力学方程更难的是一堆参数放在一起后模态结果能不能解释得通。这篇博文就围绕这个模型展开从思路、方程、参数到实操记录和踩坑经验一次性讲清楚。1. 齿轮传动扭振分析的思路拆解1.1 传动扭振问题的工程背景齿轮传动系统在旋转机械里太常见了风电齿轮箱、直升机主减速器、汽车自动变速箱、工业机器人关节核心都是齿轮。正常情况下系统平稳转动但只要转速、负载稍有波动整个轴系就会发生扭转振动。所谓扭振就是各旋转构件沿着旋转方向的角向振动表现出来的是转速波动、扭矩波动、齿轮敲击严重时导致齿面疲劳、断齿、轴承损伤。做传动系统设计的时候静态强度校核只是第一步动态特性如果不过关照样出大问题。比如某型齿轮箱在某个转速区间噪声特别大拆开一看齿面有微动磨损这就是扭振模态被激励起来了系统发生了共振。所以我们需要在设计阶段就建立动力学模型算出系统的固有频率和振型也就是扭振模态再判断工作转速范围是否避开了共振区间。NewMyModel 这个项目就是做这个事的。它面向的是一个包含行星齿轮级的传动系统建模目标是获得系统的动力学方程求解扭振模态再进一步分析不同参数对固有频率的影响。1.2 为什么聚焦行星齿轮行星齿轮和普通定轴齿轮最大的区别在于多个行星轮同时参与啮合功率分流行星架作为输出或输入构件又有公转运动。这个构型的好处是体积小、传动比大、承载能力强但也带来了一个麻烦——动力学行为比定轴齿轮复杂得多。定轴齿轮每一对啮合的路径是固定的齿轮轴的位置固定分析相对简单。而行星齿轮里太阳轮、齿圈、行星架、行星轮四类构件全都耦合在一起行星轮既要自转又要公转啮合副有两类太阳轮-行星轮外啮合行星轮-齿圈内啮合。如果有 N 个行星轮就有 2N 条啮合路径。这些啮合路径之间还互相影响行星轮的位置差异、载荷不均都会反映到动力学方程里。做扭振分析时如果我们只是想知道“这个减速器轴系的扭转固有频率是多少”最核心的建模对象是各构件的扭转自由度而行星轮的额外数量会显著增加系统的自由度和方程规模也需要额外小心处理啮合相位关系。这也是 NewMyModel 选择行星齿轮作为对象的原因它比定轴齿轮更有代表性也更贴近实际工程。1.3 建模选型集中质量法还是有限元法很多新手拿到问题第一反应是直接上有限元软件画个全三维模型算个模态。这个思路没错但代价太大。齿轮系统尤其是行星轮系接触非线性、时变啮合刚度、多体运动关系有限元算一轮得几小时甚至几天参数优化的时候根本跑不动。NewMyModel 采用的是集中质量法。把每个齿轮构件简化成具有转动惯量的刚体只保留旋转自由度啮合关系用弹簧-阻尼单元表示。相当于把整个传动系统抽象成一个“圆盘-轴-弹簧”的扭振系统。这种方法物理概念清晰方程规模小计算速度快特别适合做参数扫描和方案对比。集中质量法的缺点是忽略构件的弹性变形、齿向载荷分布不均等细节精度不如有限元。但它的价值在于快速判断固有频率区间和振型趋势。在工程上我们通常先用集中质量模型扫一圈发现问题后再用有限元或试验验证关键工况。这是性价比最高的路子。2. 核心细节解析与实操要点2.1 行星齿轮系统的自由度与坐标定义建立动力学方程之前第一步是确定自由度。对纯扭转模型来说每个旋转构件只取一个扭转角度作为广义坐标。一个典型的单级行星排包括太阳轮转角记为 θs齿圈转角记为 θr行星架转角记为 θc第 i 个行星轮转角记为 θpii 1, 2, ..., N如果齿圈固定那 θr 0不是自由度。但如果齿圈浮动或者作为输出构件就要保留。NewMyModel 里我将齿圈作为固定件所以自由度构成是太阳轮、行星架、N 个行星轮。对于 3 个行星轮的系统一共是 5 个自由度如果齿圈不固定就是 6 个自由度。坐标不一定要直接用角度。为了和啮合刚度单位统一建议把角度量纲转成线位移。啮合线方向上的等价位移定义为x_s r_bs * θsx_ci r_bc * θcx_pi r_bp * θpi其中 r_bs 是太阳轮基圆半径r_bp 是行星轮基圆半径r_bc 是行星架中心到行星轮中心的距离。这样处理后所有广义坐标的量纲都是米刚度项的单位统一为 N/m装配质量矩阵时不会出现量纲混乱。2.2 转动惯量与刚度参数怎么定集中质量模型里转动惯量是基础参数。太阳轮、齿圈、行星轮、行星架的转动惯量可以通过三维模型测量或者用解析公式估算。工程上如果用 CAD 模型直接查质量属性即可如果没有模型可以按齿轮几何尺寸估算。啮合刚度是另一个关键参数。齿轮啮合时齿对的啮合刚度会随啮合位置周期性变化也就是时变啮合刚度。这个变化是齿轮振动的主要激励来源之一。NewMyModel 里我用了两种刚度取值方式解析公式估算和有限元接触计算。解析估算可以用 ISO 6336 标准里的啮合刚度公式或者简化计算。比如单齿对啮合刚度的范围大概在 1×10^7 ~ 5×10^8 N/m具体与模数、齿宽、材料有关。我们做模态分析时可以先取平均啮合刚度忽略时变波动这样得到一个“静态”固有频率再考虑刚度上下限看固有频率带有多宽。阻尼参数的确定最难。实测阻尼很麻烦设计阶段常用模态阻尼比来近似。钢制齿轮箱的扭转模态阻尼比一般在 0.01~0.05 之间。阻尼对固有频率影响很小主要通过复特征值法计算衰减率。NewMyModel 里我最初直接忽略阻尼只做无阻尼模态分析后来才根据试验数据修正了阻尼矩阵。2.3 时变啮合刚度的处理策略上面提到啮合刚度是时变的但模态分析通常是线性时不变理论。如果要严格考虑时变效应就不能简单地用常数刚度矩阵做特征值求解因为刚度矩阵每个时刻都在变。这会引出周期时变系统的稳定性问题比如参数共振。不过在工程设计的初始阶段我们一般做两个层次的近似第一层把啮合刚度取为平均值直接算系统固有频率用于判断是否存在共振风险。这个方法简单快速适合初步筛选。如果工作频率离固有频率足够远基本可以放心。第二层把啮合刚度考虑为方波或简谐波动的周期函数用多体动力学仿真或者 Floquet 理论分析稳定性边界。这能判断系统是否会发生参数激振尤其是行星轮个数较多时啮合相位关系会让某些谐波激励变得很强。NewMyModel 在模态分析这步用的是平均刚度但在后续瞬态分析里保留了完整时变刚度曲线。两部分配合既保证模态结果稳定可解释又不会丢失关键激励特征。2.4 扭振模态的物理意义我们求出的模态每一个都对应一个固有频率和一组振型向量。振型向量表示各自由度在这个频率下的相对振幅和相位。用来判断哪个构件振动大哪个构件是节点。比如某个扭振模态下太阳轮的振幅很大而行星架几乎不动说明这个模态主要发生在太阳轮-行星轮啮合副附近。如果工作转速的啮合频率刚好和这个固有频率重合那太阳轮侧的振动会显著放大齿面载荷波动也大。对于行星齿轮还会出现一种“纯扭转模态”所有构件只做扭转振动中心构件没有横向位移此时所有行星轮的运动是等相位的。还有一种“行星轮模态”不同行星轮之间的振动相位不一致振型更加复杂。NewMyModel 的计算结果里我特意把这两种模态区分开不然单看固有频率不知道到底哪个构件危险。3. 实操过程与核心环节实现3.1 建立单级行星排的动力学方程直接上方程。对于齿圈固定、N 个行星轮的单级行星排广义坐标向量取为q [x_s, x_c, x_p1, x_p2, ..., x_pN]^T其中 x_pi 是第 i 个行星轮在啮合线方向的等价位移。太阳轮与行星轮之间的啮合刚度为 k_spi行星轮与齿圈之间的啮合刚度为 k_rpi。对于第 i 个行星轮它与太阳轮的相对位移沿啮合线可以写为δ_spi x_s - x_pi - x_c * cos(α)这里的 α 是啮合角。注意行星架的位移也会带入到啮合线方向的相对位移里因为行星架的公转运动影响了太阳轮和行星轮的中心距投影。行星轮与齿圈的相对啮合位移写为δ_rpi x_pi - x_c * cos(α)有了相对位移啮合力就是刚度乘以相对位移然后按虚功原理组装到各自由度的运动方程。写出行星轮 i 的力矩平衡方程后再综合所有构件得到整体矩阵形式的动力学方程M q C q K q F其中 M 是质量矩阵主要由各构件转动惯量换算后的等效质量组成K 是刚度矩阵由啮合刚度项组合而成C 是阻尼矩阵F 是外部扭矩、负载波动等激励向量。需要注意由于啮合刚度的时变性K 实际上是随时间变化的矩阵 K(t)但在模态分析中我把它取为时均值所以可以写成定常矩阵。3.2 组装质量矩阵、刚度矩阵以 N3、齿圈固定的行星排为例我给出实际组装过程中的矩阵片段。质量矩阵 M 是对角阵M diag(m_s, m_c, m_p1, m_p2, m_p3)其中 m_s J_s / r_bs^2m_c J_c / r_bc^2m_pi J_p / r_bp^2。由于我们用了等价线位移这里换算后的等效质量单位是 kg和刚度单位 N/m 配套。刚度矩阵 K 是一个对称矩阵。对于第 i 个行星轮太阳轮-行星轮啮合副的刚度 k_spi 会同时贡献到太阳轮和行星轮对应的行列行星轮-齿圈副的刚度 k_rpi 会贡献到行星轮对应的行列两个啮合副还通过行星架位移耦合起来。整理后的 K 可以写成三部分叠加K K_sp K_rp K_cK_sp 包含太阳轮-行星轮啮合刚度K_rp 包含行星轮-齿圈啮合刚度K_c 是行星架引起的耦合刚度项。具体的元素不在这里全部展开但组装时有一个原则每个啮合副的刚度按“位移差”二次型展开后填入对应的位置。比如某个相对位移 δ x_a - x_b那么刚度 k 在矩阵中会在 (a,a) 位置加 k(b,b) 位置加 k(a,b) 和 (b,a) 位置加 -k。这个规则非常简单但非常容易犯错尤其是行星架耦合项很多需要细心整理。3.3 求解扭振模态特征值问题忽略阻尼求无阻尼自由振动模态。将方程设为M q K q 0设解为 q φ e^{jωt}代入后得到广义特征值问题K φ λ M φ其中 λ ω^2ω 是圆频率f ω/(2π) 是对应固有频率φ 是振型向量。求解这一特征值问题我用的是 Python 的 scipy.linalg.eigh 函数可以直接处理广义特征值问题并且返回按特征值升序排列的结果。如果矩阵规模不大自由度几十以内这个方式非常快。早期版本我用过 MATLAB 的 eig(K, M)结果也是一样的。求解得到的固有频率中通常会有零频率对应刚体模态。对于纯粹的扭转模型如果系统没有任何约束至少会有一个刚体模态——整体旋转不产生弹性变形。齿圈固定后自由度被约束刚体模态会消失但要注意检查坐标定义是否真的消除了刚体位移。下面是我在 NewMyModel 里的一段伪代码用于组装矩阵后求解模态import numpy as np from scipy.linalg import eigh # 自由度个数 n 5 # 质量矩阵对角 M np.diag([m_s, m_c, m_p1, m_p2, m_p3]) # 刚度矩阵由啮合刚度组装这里为示意 K np.zeros((n, n)) # ... 填充 K ... # 求解广义特征值问题 eigenvalues, eigenvectors eigh(K, M) # 固有频率Hz freqs np.sqrt(np.maximum(eigenvalues, 0)) / (2 * np.pi) # 输出前几阶模态 for i in range(len(freqs)): print(fModal {i1}: {freqs[i]:.2f} Hz) print(eigenvectors[:, i])注意eigh返回的是列向量为特征向量的矩阵使用时要留意每一列的归一化方式。振型向量可以乘以任意常数所以归一化只是方便比较相对振幅。3.4 NewMyModel 案例参数与结果示例为了直观我给出一组示例参数。单级行星排齿圈固定3 个行星轮均匀分布。基本参数模数 m 2 mm太阳轮齿数 zs 24行星轮齿数 zp 36齿圈齿数 zr 96zr zs 2zp压力角 α 20°齿宽 b 20 mm转动惯量经三维模型测量后近似取值太阳轮 Js 0.002 kg·m²行星轮 Jp 0.005 kg·m²行星架 Jc 0.05 kg·m²啮合刚度取平均时变刚度太阳轮-行星轮啮合刚度 k_sp 8.5 × 10^7 N/m行星轮-齿圈啮合刚度 k_rp 1.2 × 10^8 N/m换算成等效质量m_s Js / r_bs^2其中 r_bs m * zs / 2 * cos(α) 2 * 24 / 2 * 0.94 ≈ 22.56 mm 0.02256 mm_s 0.002 / 0.02256^2 ≈ 3.93 kg行星轮基圆半径 r_bp m * zp / 2 * cos(α) 236/20.94 ≈ 33.84 mm等效质量 m_p 0.005 / 0.03384^2 ≈ 4.36 kg。行星架等效半径取中心距 r_c m*(zszp)/2 2*(2436)/2 60 mm 0.06 mm_c 0.05 / 0.06^2 ≈ 13.89 kg。求解得到的前几阶固有频率大约在阶次固有频率 (Hz)振型特征1约 620 Hz太阳轮与行星轮主振动2约 890 Hz行星架参与明显3约 1530 Hz行星轮之间反相振动这个数值只是示例实际结果会和啮合刚度、转动惯量偏差有关系。但趋势可以参考固有频率随刚度增大而增大随惯量增大而减小。系统中如果存在多个相同的行星轮会因为对称性出现重根模态也就是两个频率相同但振型不同的模态。4. 常见问题与排查技巧实录4.1 刚度矩阵奇异出现多余零频第一次组装 NewMyModel 时我求出来的固有频率里多了两个零频明显不对。后来排查发现行星架与行星轮之间的耦合项少加了一部分。虽然齿圈固定但如果行星架的位移和行星轮位移之间存在未被约束的组合系统仍可能存在刚体自由度或机构自由度。排查技巧检查刚度矩阵的秩。如果自由度数为 n排除刚体模态后刚度矩阵的秩应该接近 n 或者 n-1根据约束情况。用 Python 的np.linalg.matrix_rank(K)可以直接看到如果秩比理论值小说明缺少约束或组装有误。另外可以给每个模态画一下振型图观察是否存在所有构件同向转动的刚体模式。如果存在即使频率不为零也要检查是不是某个刚度元素漏填了。4.2 单位不一致导致频率数量级离谱早前时候我直接用 Jkg·m²和 kN/m混着填矩阵结果固有频率算出来是几千赫兹甚至上万赫兹完全不符合常理。问题出在如果广坐标直接取角度转动方程里应该是扭矩 J * θ而啮合力产生的力矩需要乘以力臂刚度也要换算成扭转刚度 N·m/rad。如果我把线位移形式的刚度和角位移坐标混用单位就乱了。解决方法有两类统一用角坐标刚度用扭转刚度 k_t k * r_b^2质量直接用 J。统一用线位移坐标刚度用线刚度 k质量用等效质量 J / r_b^2。我建议首选后者因为啮合刚度通常以 N/m 给出线位移坐标更不容易出错。组装完以后检查一下量纲M 的单位是 kgK 的单位是 N/mω² 的单位应该是 1/s²这样频率才是 Hz。4.3 模态频率密集振型难以区分行星齿轮系统因为对称性会出现频率非常接近的模态对比如两个行星轮反相模态的频率可能只差几赫兹。如果模态分析时只提取频率不仔细看振型容易把两个不同物理意义的模态搞混。一种实用方法把振型按自由度分类分别查看中心构件和行星轮分量的幅值比例。比如某个模态如果所有行星轮幅值几乎相等且同相就是纯扭转模态如果行星轮之间有 180° 相位差就是行星轮模态。NewMyModel 里我加了一个简单的后处理按相位关系自动分组省了不少事。如果后续要做响应分析还要注意这些密集模态对在激励下会同时被激起响应可能是两个模态叠加而不是单一模态。这个时候光看固有频率不够要算振型参与因子或模态力。4.4 时变刚度导致模态分析结果不稳定如果直接把某个瞬时时刻的时变刚度矩阵拿去算模态会发现固有频率随着时间在一个区间内波动。这本身是正常现象不是数值 bug。但如果波动范围太大说明系统参数设计存在潜在共振风险需要对啮合刚度波动幅值做控制。实际工程中齿轮修形、齿廓设计都会影响时变刚度波动。NewMyModel 里我做了一个参数扫描把啮合刚度在最低、平均、最高三个值下分别算模态观察固有频率带。如果某个激励频率落在带内就得进一步做瞬态分析而不是只靠平均刚度结论。4.5 模型验证的常用手段模型算完了总要回答一个问题算得准不准。如果手头有试验台最直接的方法是用锤击模态试验或者工作变形分析ODS测出实际固有频率。如果只有仿真也可以用有限元模型交叉验证。我在 NewMyModel 项目中做过的验证有三个层次第一简化模型验证。把行星轮个数设为 1自由度减少用解析公式或手算校核部分特征值。第二与有限元模型对比。在 ANSYS 里建立简化的实体齿轮系做模态分析对比前几阶频率。只要误差在 10% 以内基本可以接受。第三与试验数据对比。实测时需要注意齿轮箱边界条件比如是否固定底座、联轴器刚度等如果试验约束和模型不一致频率偏差会很大。这个验证思路不仅适用于扭振模型也适用于其他传动系统的建模项目。模型的价值在于预测和优化如果连基准工况都对不上后面所有的参数影响分析都不可信。5. 一些值得记住的经验写到这里我想起 NewMyModel 项目初期犯过的最大错误一开始太想追求“全维度”的动力学模型把太阳轮、齿圈、行星架的横向自由度都加进去自由度一下多了几十个结果调试周期被拖得很长最后砍掉横向自由度、只保留扭转自由度模型反而更清晰了。做工程动力学分析不是模型越复杂越好而是要用最简模型回答最关键问题。扭振模态分析的目标是避开共振、识别危险模态纯扭转模型足够完成这个任务。另外参数准备工作一定不能省。我吃过亏以为转动惯量随便估个量级就行结果算出来的频率和后来实测差了一倍。后来老老实实从 CAD 模型提取参数细节到轴承位置的附加惯量都算进去误差才降下来。齿轮系统的动力学分析参数质量直接决定结果质量这一步偷懒后面要花十倍的精力去排查。最后想分享一个小技巧在计算完成后把固有频率和对应振型整理成一张表再叠加到 Campbell 图或者频率-转速图上一眼就能看出工作转速下哪些激励阶次会命中固有频率。NewMyModel 里最后输出就是一张这样的图比单纯列一串数字直观得多。做齿轮传动的人最终看的不是模型文件有多漂亮而是能不能在设计阶段提前发现问题。这个模型的价值就是把那些看不见的扭振风险变成屏幕上可量化的数字。本文还有配套的精品资源点击获取