72自由度人字齿行星齿轮传动系统动力学建模。 该程序可以用于求解72自由度人字齿行星齿轮各构件间的动态响应，包括振动位移，振动速度，振动加速度，齿轮啮合力等，模型考虑了时变啮合刚度，综合误差，裂纹故障等因素。 该模型还可以用于求解系统存在不同故障时的动态响应。

行星齿轮箱这玩意儿动静大起来是真要命。前几天折腾一套72自由度的人字齿行星传动模型，差点被微分方程搞秃头。不过还真别说，这种复杂系统的动力学响应玩起来真上瘾——特别是能模拟裂纹故障的时候，那振动曲线比心电图还刺激。

先来瞅瞅模型骨架。72个自由度可不是闹着玩的，太阳轮、行星轮、齿圈各带着平移加扭转的六个自由度满场飞。核心代码里这个结构体藏着玄机：

class GearSystem:
    def __init__(self):
        self.masses = np.array([...])  # 质量矩阵
        self.k = lambda t: time_varying_stiffness(t)  # 时变刚度函数
        self.errors = [...]  # 综合误差序列
        self.crack_params = None  # 裂纹特征参数

重点是这个时变刚度函数，它直接关系到啮合冲击的计算。比如用分段谐波叠加模拟实际工况：

def time_varying_stiffness(t):
    base_k = 5e8  # 基准刚度
    fluctuation = 0.2 * np.sin(20*t) + 0.1 * np.cos(45*t)
    return base_k * (1 + 0.15 * fluctuation)

这比固定刚度模型多了个随时间波动的成分，仿真时会出现很有意思的拍振现象。不过真正要命的是组装系统方程时的稀疏矩阵操作——72x72的刚度矩阵用常规方法算会爆炸。

动力学方程用状态空间表达更带劲。下面这段向量化处理是关键：

def equations(y, t, system):
    n = len(y)//2
    displacement = y[:n]
    velocity = y[n:]
    
    K = system.k(t) * assemble_topology_matrix() 
    K += crack_effect(system.crack_params, t)  # 裂纹影响
    
    acceleration = np.linalg.solve(system.masses, 
                                  -K.dot(displacement) 
                                  - damping(velocity))
    return np.concatenate([velocity, acceleration])

这里有个骚操作：把阻尼项单独拎出来处理，避免刚度矩阵变得病态。注意np.linalg.solve的用法，比直接求逆稳定得多，特别是当行星轮数量变化时不会崩。

玩点刺激的——来点故障试试。裂纹故障的模拟其实就是在刚度矩阵里埋地雷：

def crack_effect(params, t):
    if params is None:
        return 0
    depth, pos = params
    crack_wave = np.where((t%0.1) < 0.02, depth*0.3, 0)  # 每转产生一次冲击
    return crack_wave * local_stiffness_matrix(pos)

这种脉冲式刚度衰减会在频谱上产生明显的边频带，实测数据里捕捉到这种特征就是故障的铁证。跑个带裂纹的仿真，加速度信号里立刻蹦出一串冲击峰值，跟放鞭炮似的。

最后可视化环节才是高潮。用PyQtGraph做个动态相位图，看着密密麻麻的轨迹线像烟花一样炸开：

def live_plot():
    app = pg.mkQApp()
    win = pg.GraphicsLayoutWidget()
    plot = win.addPlot()
    curve = plot.plot(pen='y')
    
    def update():
        data = get_realtime_data()  # 从求解器获取最新数据
        curve.setData(data[:,0], data[:,5])  # 太阳轮x方向位移-速度相位图
    timer = QtCore.QTimer()
    timer.timeout.connect(update)
    timer.start(50)

这比静态图直观多了，特别是观察混沌现象时，能清楚看到吸引子的形态变化。曾经在某个参数组合下看到轨迹突然从环面变成雪花状，那一刻真是汗毛倒竖。

整套模型跑下来最深的体会：行星齿轮的振动特征就像密码本，时变刚度和故障因子就是解码的关键。现在回头看看那72个自由度的方程，倒觉得它们像钢琴键——敲对了组合，故障特征曲就流淌出来了。