✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
🍊个人信条：格物致知,完整Matlab代码及仿真咨询内容私信。
🔥内容介绍
一、引言
齿轮传动系统作为机械装备的核心动力传递部件，广泛应用于航空航天、汽车工业、风电装备等领域。其运行稳定性直接决定了整机的工作性能、可靠性与使用寿命。在齿轮啮合过程中，由于齿侧间隙、时变啮合刚度、啮合误差等因素的影响，系统不可避免地产生振动与噪声，严重时会导致齿轮齿面磨损、疲劳点蚀甚至断齿等故障。因此，开展齿轮动力学振动特性研究，建立精准的动力学模型并高效求解，对齿轮传动系统的优化设计、故障诊断与振动控制具有重要的理论意义与工程价值。
齿轮动力学模型的自由度划分是刻画系统振动特性的关键，四自由度模型作为经典的简化模型，能够兼顾计算效率与分析精度，主要考虑主动轮与从动轮的扭转振动和横向振动（径向振动），可有效反映齿轮啮合过程中的核心振动形态。在动力学方程求解方面，四阶龙格-库塔法（RK4）因具有精度高、稳定性好、收敛速度快等优势，被广泛应用于非线性动力学系统的数值求解。本文基于MATLAB平台，构建四自由度齿轮动力学振动模型，采用RK4法求解系统动力学方程，分析齿轮系统的振动响应特性，为齿轮传动系统的性能优化提供理论支撑。
二、四自由度齿轮动力学振动模型构建
2.1 模型假设与自由度定义
为简化模型并聚焦核心振动特性，提出以下合理假设：（1）忽略齿轮轴的弯曲变形与轴向振动，仅考虑扭转和横向（径向）两个方向的振动；（2）齿轮为刚性圆盘，质量集中于轮心；（3）齿轮轴的支撑为弹性约束，用等效刚度系数描述；（4）考虑齿侧间隙、时变啮合刚度和啮合误差等关键非线性因素。
定义四自由度如下：以齿轮啮合平面为研究平面，设主动轮的扭转角为\( \theta_1 \)（绕自身轴线旋转）、横向位移为\( x_1 \)（沿啮合线垂直方向）；从动轮的扭转角为\( \theta_2 \)、横向位移为\( x_2 \)。四个自由度相互耦合，共同构成系统的振动状态。
2.2 系统动力学方程推导
基于拉格朗日方程（\( \frac{d}{dt}(\frac{\partial T}{\partial \dot{q}}) - \frac{\partial T}{\partial q} + \frac{\partial V}{\partial q} + \frac{\partial D}{\partial \dot{q}} = Q \)）推导系统动力学方程，其中\( T \)为系统动能，\( V \)为势能，\( D \)为阻尼能，\( q \)为广义坐标（即上述四个自由度），\( Q \)为广义外力。
2.2.1 系统动能\( T \)
系统动能由主动轮和从动轮的扭转动能与横向平动动能组成：
\( T = \frac{1}{2}J_1\dot{\theta}_1^2 + \frac{1}{2}m_1\dot{x}_1^2 + \frac{1}{2}J_2\dot{\theta}_2^2 + \frac{1}{2}m_2\dot{x}_2^2 \)
其中，\( J_1、J_2 \)分别为主动轮、从动轮的转动惯量；\( m_1、m_2 \)分别为主动轮、从动轮的质量；\( \dot{\theta}_1、\dot{\theta}_2 \)分别为主动轮、从动轮的角速度；\( \dot{x}_1、\dot{x}_2 \)分别为主动轮、从动轮的横向速度。
2.2.2 系统势能\( V \)
系统势能包括齿轮轴的扭转弹性势能、横向支撑弹性势能、齿轮啮合弹性势能以及齿侧间隙引起的非线性势能：
\( V = \frac{1}{2}k_{t1}\theta_1^2 + \frac{1}{2}k_{t2}\theta_2^2 + \frac{1}{2}k_{x1}x_1^2 + \frac{1}{2}k_{x2}x_2^2 + V_m(x_1,x_2,\theta_1,\theta_2) + V_b(x_1,x_2,\theta_1,\theta_2) \)
其中，\( k_{t1}、k_{t2} \)分别为主动轮轴、从动轮轴的扭转刚度；\( k_{x1}、k_{x2} \)分别为主动轮、从动轮的横向支撑刚度；\( V_m \)为啮合弹性势能，与啮合刚度相关；\( V_b \)为齿侧间隙势能，是典型的分段非线性函数。
2.2.3 系统阻尼能\( D \)
考虑扭转阻尼和横向阻尼，采用粘性阻尼模型：
\( D = \frac{1}{2}c_{t1}\dot{\theta}_1^2 + \frac{1}{2}c_{t2}\dot{\theta}_2^2 + \frac{1}{2}c_{x1}\dot{x}_1^2 + \frac{1}{2}c_{x2}\dot{x}_2^2 \)
其中，\( c_{t1}、c_{t2} \)分别为主动轮轴、从动轮轴的扭转阻尼系数；\( c_{x1}、c_{x2} \)分别为主动轮、从动轮的横向阻尼系数。
2.2.4 广义外力\( Q \)
广义外力主要包括主动轮输入转矩\( M_1 \)和从动轮输出负载转矩\( M_2 \)，对应扭转自由度的外力，横向自由度无外力作用时\( Q \)分量为0。
2.2.5 最终动力学方程
将\( T、V、D、Q \)代入拉格朗日方程，整理得到四自由度齿轮动力学振动方程的矩阵形式：
\( M\ddot{q} + C\dot{q} + K(q)\dot{q} + F_N(q,\dot{q}) = F \)
其中，\( M = diag[J_1, m_1, J_2, m_2] \)为质量矩阵；\( C = diag[c_{t1}, c_{x1}, c_{t2}, c_{x2}] \)为阻尼矩阵；\( K(q) \)为时变刚度矩阵，与齿轮啮合位置相关；\( F_N(q,\dot{q}) \)为非线性力向量，包含齿侧间隙、干摩擦等非线性因素；\( F = [M_1, 0, -M_2, 0]^T \)为外力向量；\( q = [\theta_1, x_1, \theta_2, x_2]^T \)为广义坐标向量。
三、四阶龙格-库塔法（RK4）求解原理
3.1 RK4法基本思想
齿轮动力学方程为二阶非线性常微分方程组，无法直接求得解析解，需采用数值方法求解。RK4法通过在每个积分步长内选取4个点，利用函数在这4个点的取值加权平均来近似求解微分方程，其核心优势是在保证计算效率的同时，能够获得较高的求解精度。
3.2 一阶化转换
为适配RK4法对一阶微分方程组的求解要求，需将二阶动力学方程转化为一阶方程组。令\( \dot{q} = z \)，则\( \ddot{q} = \dot{z} \)，代入动力学方程可得：
\( \dot{z} = M^{-1}[F - C z - K(q) q - F_N(q,z)] \)
定义状态向量\( X = [q^T, z^T]^T = [\theta_1, x_1, \theta_2, x_2, \dot{\theta}_1, \dot{x}_1, \dot{\theta}_2, \dot{x}_2]^T \)，则一阶微分方程组可表示为：
\( \dot{X} = f(t,X) \)
其中，\( f(t,X) \)为状态向量的导数函数，由上述一阶化后的方程推导得到。
3.3 RK4法求解步骤
对于一阶微分方程\( \dot{X} = f(t,X) \)，给定初始条件\( X(t_0) = X_0 \)和积分步长\( h \)，RK4法的求解步骤如下：
计算第一步斜率：\( k_1 = f(t_n, X_n) \)
计算第二步斜率：\( k_2 = f(t_n + \frac{h}{2}, X_n + \frac{h}{2}k_1) \)
计算第三步斜率：\( k_3 = f(t_n + \frac{h}{2}, X_n + \frac{h}{2}k_2) \)
计算第四步斜率：\( k_4 = f(t_n + h, X_n + h k_3) \)
更新下一时刻状态向量：\( X_{n+1} = X_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \)
重复上述步骤，即可得到整个时间域内系统的状态向量\( X(t) \)，进而获得各自由度的振动响应曲线。
四、基于MATLAB的模型实现与求解
4.1 参数设置
根据某型直齿圆柱齿轮传动系统参数，设定模型关键参数如下（示例值）：
转动惯量：\( J_1 = 0.01 \, \text{kg·m}^2 \)，\( J_2 = 0.02 \, \text{kg·m}^2 \)
质量：\( m_1 = 5 \, \text{kg} \)，\( m_2 = 8 \, \text{kg} \)
刚度系数：\( k_{t1} = 1.0×10^5 \, \text{N·m/rad} \)，\( k_{t2} = 1.5×10^5 \, \text{N·m/rad} \)，\( k_{x1} = 2.0×10^6 \, \text{N/m} \)，\( k_{x2} = 2.5×10^6 \, \text{N/m} \)，啮合刚度均值\( k_{m0} = 5.0×10^7 \, \text{N/m} \)
阻尼系数：\( c_{t1} = 10 \, \text{N·m·s/rad} \)，\( c_{t2} = 15 \, \text{N·m·s/rad} \)，\( c_{x1} = 500 \, \text{N·s/m} \)，\( c_{x2} = 800 \, \text{N·s/m} \)
外部转矩：\( M_1 = 50 \, \text{N·m} \)，\( M_2 = 48 \, \text{N·m} \)
齿侧间隙：\( b = 0.0005 \, \text{m} \)，积分步长\( h = 1×10^{-4} \, \text{s} \)，仿真时间\( t = 0 \sim 0.5 \, \text{s} \)
初始条件：\( X_0 = [0, 0, 0, 0, 0, 0, 0, 0]^T \)（初始静止状态）
4.2 MATLAB程序实现
基于上述原理，编写MATLAB程序实现模型的构建与求解，核心代码包括参数定义、状态方程函数定义、RK4法求解函数定义以及结果可视化四部分。
五、仿真结果分析
运行上述MATLAB程序，得到四自由度齿轮系统在0~0.5s内的振动响应曲线，对各自由度的振动特性分析如下：
5.1 扭转振动特性
主动轮与从动轮的扭转振动响应呈现明显的衰减振荡特性。初始时刻，由于外力矩的突然施加，扭转角迅速增大并达到峰值，随后在阻尼作用下振幅逐渐减小，最终趋于稳定。主动轮的扭转角峰值小于从动轮，这是由于主动轮转动惯量较小、扭转刚度较大，系统刚度对扭转振动的抑制作用更显著。同时，振动曲线中存在微小的周期性波动，这是时变啮合刚度引起的调制效应，体现了齿轮啮合的动态特性。
5.2 横向振动特性
主动轮与从动轮的横向振动响应同样表现为衰减振荡，但其振动频率高于扭转振动频率，这是因为横向支撑刚度远大于扭转刚度。由于齿侧间隙的存在，横向振动响应在初始阶段存在明显的非线性跳跃现象，当振动幅值超过齿侧间隙后，啮合力开始作用，振动特性逐渐趋于规律。从动轮的横向位移峰值大于主动轮，主要受质量和支撑刚度的影响，质量越大、刚度越小，振动幅值越大。
5.3 耦合振动特性
四自由度模型中，扭转振动与横向振动存在明显的耦合效应。例如，主动轮的扭转振动通过啮合关系传递给从动轮，同时引发横向振动的变化；而横向振动产生的啮合力变化又会反作用于扭转振动，导致两者的振动响应相互影响。这种耦合特性是齿轮系统振动的重要特征，也是导致齿轮磨损、噪声的关键因素之一。
六、结论与展望
本文构建了四自由度齿轮动力学振动模型，考虑了时变啮合刚度、齿侧间隙等关键非线性因素，基于MATLAB平台实现了四阶龙格-库塔法的数值求解，通过仿真分析获得了系统的扭转和横向振动响应特性。结果表明，RK4法能够精准、高效地求解齿轮非线性动力学方程，所建模型可有效反映齿轮系统的动态振动特性，为齿轮传动系统的优化设计与故障诊断提供了可靠的理论工具。
未来研究方向可聚焦于以下方面：（1）完善模型细节，考虑齿轮的制造误差、装配偏差等更多实际因素，提升模型的工程适用性；（2）优化数值求解算法，结合自适应步长RK4法或其他高精度数值方法，进一步提高求解效率与精度；（3）开展实验验证，通过齿轮台架试验测试系统的振动响应，对比仿真结果与实验数据，修正并验证模型的准确性；（4）基于模型开展振动控制研究，设计合理的控制策略，降低齿轮系统的振动与噪声。
⛳️ 运行结果
🔗 参考文献
[1] 代洪华.非线性气动弹性系统求解方法及复杂响应研究[D].西北工业大学,2014.DOI:10.7666/d.D689390.
[2] 范威威.连续系统模型分布并行仿真算法研究[D].国防科学技术大学[2025-12-22].DOI:10.7666/d.D676200.
[3] 蒋珉,黄振全,王静.具有最大稳定域的实时RK4公式[J].系统仿真学报, 2006, 18(2):4.DOI:10.3969/j.issn.1004-731X.2006.02.010.
📣 部分代码
🎈 部分理论引用网络文献，若有侵权联系博主删除
👇 关注我领取海量matlab电子书和数学建模资料
🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：
🌈 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱调度、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划（2E-VRP）、充电车辆路径规划（EVRP）、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题、港口调度、港口岸桥调度、停机位分配、机场航班调度、泄漏源定位
🌈 机器学习和深度学习时序、回归、分类、聚类和降维
2.1 bp时序、回归预测和分类
2.2 ENS声神经网络时序、回归预测和分类
2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类
2.4 CNN|TCN|GCN卷积神经网络系列时序、回归预测和分类
2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类
2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类
2.7 ELMAN递归神经网络时序、回归\预测和分类
2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类
2.9 RBF径向基神经网络时序、回归预测和分类
2.10 DBN深度置信网络时序、回归预测和分类
2.11 FNN模糊神经网络时序、回归预测
2.12 RF随机森林时序、回归预测和分类
2.13 BLS宽度学习时序、回归预测和分类
2.14 PNN脉冲神经网络分类
2.15 模糊小波神经网络预测和分类
2.16 时序、回归预测和分类
2.17 时序、回归预测预测和分类
2.18 XGBOOST集成学习时序、回归预测预测和分类
2.19 Transform各类组合时序、回归预测预测和分类
方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断
🌈图像处理方面
图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知
🌈 路径规划方面
旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻
🌈 无人机应用方面
无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划
🌈 通信方面
传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化、水声通信、通信上传下载分配
🌈 信号处理方面
信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化、心电信号、DOA估计、编码译码、变分模态分解、管道泄漏、滤波器、数字信号处理+传输+分析+去噪、数字信号调制、误码率、信号估计、DTMF、信号检测
🌈电力系统方面
微电网优化、无功优化、配电网重构、储能配置、有序充电、MPPT优化、家庭用电
🌈 元胞自动机方面
交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀
🌈 雷达方面
卡尔曼滤波跟踪、航迹关联、航迹融合、SOC估计、阵列优化、NLOS识别
🌈 车间调度
零等待流水车间调度问题NWFSP、置换流水车间调度问题PFSP、混合流水车间调度问题HFSP、零空闲流水车间调度问题NIFSP、分布式置换流水车间调度问题 DPFSP、阻塞流水车间调度问题BFSP
👇
