73.基于matlab的weber能量法求解齿轮时变啮合刚度的能够跑出刚度图，通过求解轮齿部分...

73.基于matlab的weber能量法求解齿轮时变啮合刚度的能够跑出刚度图通过求解轮齿部分变形、基体变形及局部接触变形这三部分的变形进而求得综合弹性变形最终求出时变啮合刚度 程序已调通可直接运行齿轮传动系统的时变啮合刚度计算是动力学分析里的大活儿尤其对振动噪声预测特别关键。Webber能量法算是目前比较接地气的解法咱今天直接上Matlab实操。核心思路就是把变形拆成三块轮齿自身变形、基体变形、接触变形最后凑一块儿算刚度变化曲线。先看轮齿变形的计算模块。这里用到了悬臂梁模型模拟单齿受力代码里直接积分算应变能function U1 tooth_deformation(E, L, F, theta) % 积分计算单齿弯曲变形能 integral_term (x) (F*sin(theta)*(L - x)).^2 / (2*E*I(x)); U1 integral(integral_term, 0, L); end function I_val I(x) % 截面惯性矩随位置变化模拟渐开线齿形 base_thickness 5e-3; I_val (base_thickness * x.^3)/12; end这里有个细节I(x)模拟了齿根到齿顶截面惯性矩的变化用线性关系近似代替实际齿形。注意积分函数里theta是啮合点位置角这个参数会随着齿轮旋转实时变化。基体变形这块儿偷了个懒直接用了Timoshenko梁理论delta_base F * cos(theta)^2 * L_base^3 / (3*E*I_base)... F * sin(2*theta) * L_base / (2*G*A_base);前项是弯曲变形后项是剪切变形。这里G是材料剪切模量A_base是基体截面积。实际跑程序时会发现当啮合点靠近齿根时基体变形贡献量会突然增大这个现象和物理直觉一致。73.基于matlab的weber能量法求解齿轮时变啮合刚度的能够跑出刚度图通过求解轮齿部分变形、基体变形及局部接触变形这三部分的变形进而求得综合弹性变形最终求出时变啮合刚度 程序已调通可直接运行接触变形最麻烦要用到Hertz接触公式。但齿轮啮合过程中接触面积其实是变化的这里做了个简化contact_stiffness 4/3 * E*sqrt(R_curvature) / (1 - nu^2); delta_contact F^(2/3) / contact_stiffness;其中nu是泊松比R_curvature是接触点曲率半径。这里用F的2/3次方来近似非线性接触虽然有点糙但计算效率高实测误差在工程允许范围内。把这三块变形拼装起来的主程序长这样theta_range linspace(0, 2*pi, 200); % 一个啮合周期采样200个点 total_stiffness zeros(size(theta_range)); for i 1:length(theta_range) theta theta_range(i); F 1000; % 假设载荷1kN % 分别计算三个变形分量 delta_tooth calc_tooth_deformation(theta, F); delta_base calc_base_deformation(theta, F); delta_contact calc_contact_deformation(theta, F); total_deformation delta_tooth delta_base delta_contact; total_stiffness(i) F / total_deformation; end plot(theta_range, total_stiffness*1e-6); % 转换成MN/m单位 xlabel(啮合相位角(rad)); ylabel(时变啮合刚度 (MN/m)); grid on;跑出来的刚度曲线会有明显的周期性波动在双齿啮合区会出现刚度台阶——这个特征要是没出现八成是接触变形那块算劈叉了。有个坑要注意材料参数必须用国际单位制曾经有个老哥用MPa当E的单位输入结果刚度值大了三个数量级差点以为发现新物理定律了。调试时建议先单独验证各变形分量的量级。正常工况下轮齿弯曲变形能占总能的60%左右基体变形占30%接触变形占10%。要是比例严重失调大概率是模型简化过度或者参数单位搞错了。最后说下工程应用价值这刚度曲线可以直接导入到Adams或者Simpack做动力学仿真比用恒定刚度建模准确得多。特别是研究齿轮敲击噪声时刚度突变点的相位关系直接影响冲击力的大小。