求助：Scipy odeint求解Jeffcot转子滑动轴承耦合ODE结果异常

我正在研究转子动力学，需使用Scipy的odeint函数求解带滑动轴承的Jeffcot转子的4个耦合常微分方程，目标是获取轴的位置uy、uz以及轴承内轴中心位置ym、zm。

系统控制方程 #

控制方程包含两组二阶耦合常微分方程：

• 轴的平动方程（uy、uz方向，考虑偏心激励、阻尼、轴与轴承的刚度耦合）
• 滑动轴承的运动方程（ym、zm方向，考虑轴与轴承的刚度耦合、轴承自身的刚度与阻尼）

现有求解代码 #

def dSdt(S, t):
    uy, uyt, uz, uzt, ym, ymt, zm, zmt = S
    
    duy_dt = uyt
    duyt_dt = (me * omega1**2 * np.cos(omega1 * t) - c * uyt + K * (ym - uy)) / mdisco

    duz_dt = uzt
    duzt_dt = (me * omega1**2 * np.sin(omega1 * t) - mdisco * g - c * uzt + K * (zm - uz)) / mdisco
    
    dym_dt = (K * (ym - uy) - 2 * (kyy * ym + kyz * zm + dyz * zmt)) / (2*dyy)
    dzm_dt = (K * (zm - uz) - 2 * (kzy * ym + kzz * zm + dzy * ymt)) / (2*dzz)
    
    return [uyt, duyt_dt, uzt, duzt_dt, dym_dt, ymt, dzm_dt, zmt]

问题现象 #

所有参数已在调用odeint前计算完成，但求解结果不符合预期：绘制uy-uz、ym-zm的轴中心轨迹图形态异常。当前omega1=5Hz时，已设定转子质量mdisco、偏心质量me、阻尼系数c、轴刚度K，以及滑动轴承的刚度系数kyy/kyz/kzy/kzz、阻尼系数dyy/dyz/dzy/dzz等参数。

核心排查点 #

1. ODE状态量导数映射完全错误
   代码中对轴承状态量的导数逻辑颠倒：

   • 状态量ym的导数应为其速度ymt，而非当前代码中计算的加速度表达式；
   • 状态量ymt的导数才是轴承y方向的加速度，也就是当前代码中赋值给dym_dt的表达式；
   • zm方向存在完全相同的错误。
     正确的导数逻辑应为：

   # 轴承y方向导数
   dym_dt = ymt
   dymt_dt = (K * (ym - uy) - 2 * (kyy * ym + kyz * zm + dyz * zmt)) / (2*dyy)
   # 轴承z方向导数
   dzm_dt = zmt
   dzmt_dt = (K * (zm - uz) - 2 * (kzy * ym + kzz * zm + dzy * ymt)) / (2*dzz)
   

   返回列表需对应修改为：

   return [uyt, duyt_dt, uzt, duzt_dt, dym_dt, dymt_dt, dzm_dt, dzmt_dt]
   

2. 激励频率单位错误
   omega1给定为5Hz，需转换为弧度/秒（omega1 = 5 * 2 * np.pi），否则激励项的频率量级偏差6倍以上，直接导致响应轨迹异常。

3. 控制方程符号与系数验证
   对照原始控制方程，检查：

   • 轴承阻尼项（如dyz * zmt、dzy * ymt）的符号是否与方程一致；
   • 分母2*dyy、2*dzz的系数是否符合推导逻辑，是否应为单个阻尼系数而非2倍。

4. 初始条件合理性
   确认odeint的初始状态（如初始位移、速度）是否符合物理场景，不合理的初始条件会导致非稳态响应轨迹。

5. 参数数值一致性
   核对滑动轴承刚度、阻尼与轴刚度K的数值量级，轴承刚度通常远大于轴刚度，若参数颠倒会导致系统响应不符合预期。

内容的提问来源于stack exchange，提问作者NLGenin