齿轮基础参数
基于传统势能法含裂纹斜齿轮时变啮合刚度（裂纹斜齿轮）,代码保证运行无问题，出图效果如页面简介
齿轮传动系统里最怕遇到啥？裂纹呗！尤其是斜齿轮这种接触线斜着走的家伙，一旦出现裂纹整个时变刚度曲线直接抽风。咱今天不整虚的，直接上代码手把手拆解传统势能法的实现套路。
先看核心思路：把轮齿想象成悬臂梁，通过计算各微段的弹性变形能反推刚度。斜齿轮的骚操作在于接触线会沿着齿面移动，这个动态过程必须用数值积分来捕捉。
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函数处理裂纹区域的刚度衰减，当接触线移动到裂纹区(y>0.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，结果刚度曲线抖得跟心电图似的。数值稳定性这块，建议加上移动平均滤波处理，毕竟实际工况总有噪声干扰。