基于传统势能法含裂纹斜齿轮时变啮合刚度（裂纹斜齿轮）,代码保证运行无问题，出图效果如页面简介

齿轮传动系统里最怕遇到啥？裂纹呗！尤其是斜齿轮这种接触线斜着走的家伙，一旦出现裂纹整个时变刚度曲线直接抽风。咱今天不整虚的，直接上代码手把手拆解传统势能法的实现套路。

先看核心思路：把轮齿想象成悬臂梁，通过计算各微段的弹性变形能反推刚度。斜齿轮的骚操作在于接触线会沿着齿面移动，这个动态过程必须用数值积分来捕捉。

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

这段代码藏着几个魔鬼细节：

1. crack_factor函数处理裂纹区域的刚度衰减，当接触线移动到裂纹区(y>0.6)时刚度直接打折
2. energy_integrand里同时考虑了弯曲变形和剪切变形，1.2那个系数是剪切修正因子
3. 积分变量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，结果刚度曲线抖得跟心电图似的。数值稳定性这块，建议加上移动平均滤波处理，毕竟实际工况总有噪声干扰。