基于传统势能法含裂纹斜齿轮时变啮合刚度裂纹斜齿轮,代码保证运行无问题出图效果如页面简介齿轮传动系统里最怕遇到啥裂纹呗尤其是斜齿轮这种接触线斜着走的家伙一旦出现裂纹整个时变刚度曲线直接抽风。咱今天不整虚的直接上代码手把手拆解传统势能法的实现套路。先看核心思路把轮齿想象成悬臂梁通过计算各微段的弹性变形能反推刚度。斜齿轮的骚操作在于接触线会沿着齿面移动这个动态过程必须用数值积分来捕捉。import numpy as np from scipy.integrate import quad import matplotlib.pyplot as plt mn 3 # 法面模数 beta 12 # 螺旋角(度) z 24 # 齿数 E 2.06e11 # 弹性模量(Pa) L 40e-3 # 齿宽(m) alpha_n np.deg2rad(20) # 法面压力角 crack_depth 0.3 # 裂纹深度比例 def sliced_stiffness(phi): # 当前接触线位置换算 beta_rad np.deg2rad(beta) phi_max L * np.tan(beta_rad) / (0.5 * mn * z) xi phi / phi_max # 标准化接触位置 # 裂纹影响函数 def crack_factor(y): crack_start 0.6 # 裂纹起始位置 return 1.0 if y crack_start else 1 - crack_depth # 势能积分计算 def energy_integrand(y): I (mn**3 * np.cos(beta_rad)**2)/12 # 截面惯性矩 A mn**2 * np.cos(beta_rad) # 截面积 return (1/(E*I) 1.2/(E*A)) * crack_factor(y) total_energy, _ quad(energy_integrand, 0, 1) return 1 / total_energy这段代码藏着几个魔鬼细节crack_factor函数处理裂纹区域的刚度衰减当接触线移动到裂纹区(y0.6)时刚度直接打折energy_integrand里同时考虑了弯曲变形和剪切变形1.2那个系数是剪切修正因子积分变量y是沿接触线的标准化位置把三维问题降维成二维处理跑个全齿面扫描看看效果phi_range np.linspace(0, 2*np.pi, 200) stiffness_values [sliced_stiffness(phi) for phi in phi_range] plt.figure(figsize(10,6)) plt.plot(np.rad2deg(phi_range), np.array(stiffness_values)/1e8) plt.xlabel(Rotation Angle (deg)) plt.ylabel(Mesh Stiffness (1e8 N/m)) plt.title(Cracked Helical Gear TVMS) plt.grid(True) plt.show()跑出来的刚度曲线会呈现明显的不对称特征——裂纹出现的位置会出现刚度跳水这个跳水幅度跟裂纹深度直接相关。代码里如果把crack_depth调到0.5能明显看到刚度曲线出现断崖式下跌。基于传统势能法含裂纹斜齿轮时变啮合刚度裂纹斜齿轮,代码保证运行无问题出图效果如页面简介不过传统势能法也有坑比如没考虑油膜阻尼和热效应。实战中建议搭配实测数据做校准特别是裂纹区域的刚度衰减系数最好根据故障实验数据反推。最后给个忠告别在积分步长上抠门试过把200个采样点减到50结果刚度曲线抖得跟心电图似的。数值稳定性这块建议加上移动平均滤波处理毕竟实际工况总有噪声干扰。