solve_ivp陷入无限循环排查：固定床反应器ODE求解异常

问题根源分析 #

从调试现象来看，流量参数（Fet、Fac等）无变化，反应速率r1/r2和分压、总压周期性振荡，核心问题集中在代码笔误、冗余状态变量、物性参数计算错误和求解器选型不当四个方面：

1. 微分方程笔误 #

dPc_dW的计算式存在逻辑错误：

# 错误代码
dPc_dW  = (Fc*P*dFc_dW -Fc*P*dFc_dW +Fc*Ft*dP_dW)/Ft**2

减号后重复使用dFc_dW，导致dPc_dW计算完全失效，直接引发分压振荡。

2. 冗余状态变量 #

将分压（Pet、Pac等）作为独立状态变量属于冗余设计——分压可通过流量和总压直接推导：Pet = (Fet/Ft)*P。冗余变量会引入代数环，导致数值求解时出现振荡或循环。

3. 物性参数计算错误 #

初始密度rho_0的计算式不符合理想气体定律，原代码分母的10.73+731.07无物理意义，导致beta计算错误，进而影响压降导数dP_dW的正确性。

4. 求解器选型不当 #

固定床反应器的ODE属于刚性方程组，默认的RK45求解器对刚性问题处理能力不足，容易陷入数值循环或无法收敛。

具体修复步骤 #

1. 修正微分方程笔误 #

将dPc_dW改为正确形式：

dPc_dW  = (Ft*P*dFc_dW - Fc*P*dFt_dW + Fc*Ft*dP_dW)/Ft**2

2. 移除冗余状态变量 #

将状态变量从14个精简为8个：[Fet, Fac, Fo, Fp, Fa, Fc, Ft, P]，分压在需要时通过流量和总压计算，无需作为状态变量求解。修改后的微分方程函数如下：

def f(W, y):
    Fet, Fac, Fo, Fp, Fa, Fc, Ft, P = y 

    # 计算分压（不再作为状态变量）
    Pet = (Fet / Ft) * P if Ft != 0 else 0
    Pac = (Fac / Ft) * P if Ft != 0 else 0
    Po = (Fo / Ft) * P if Ft != 0 else 0
    Pp = (Fp / Ft) * P if Ft != 0 else 0
    Pa = (Fa / Ft) * P if Ft != 0 else 0
    Pc = (Fc / Ft) * P if Ft != 0 else 0

    r1 = 0.1036*np.exp(-3674/T)*((Po*Pet*Pac*(1+1.7*Pa)))/((1+0.583*Po*(1+1.7*Pa))*(1+6.8*Pac))
    r2 = 1.9365*(10**5)*np.exp(-10116/T)*(Po*(1+0.68*Pa))/(1+0.76*Po*(1+0.68*Pa))

    m = Fet * 28.06 + Fac * 60.052 + Fo * 32 + Fp * 86.09 + Fa * 18.02 + Fc * 44.01
    G = m / A_c                                                                     
    beta = (1.75 * G**2 / (rho_0 * g_c * D_p)) * ((1-phi) / phi**3)

    # 流量导数
    dFet_dW = -r1 - r2
    dFac_dW = -r1
    dFo_dW  = -r1/2 - 3*r2
    dFp_dW  = r1
    dFa_dW  = r1 + 2*r2
    dFc_dW  = 2*r2
    dFt_dW  = -r1/2
    
    # 压降导数
    dP_dW   = -(beta/(A_c*(1-phi)*rho_c))*(P_0/P)*(Ft/Ft_0) if P != 0 else 0
    
    return np.array([dFet_dW, dFac_dW, dFo_dW, dFp_dW, dFa_dW, dFc_dW, dFt_dW, dP_dW])

3. 修正物性参数计算 #

按照理想气体密度公式修正rho_0，加入温度单位转换（K转°R，匹配气体常数单位）：

T_r = T * 9/5  # 转换为°R
rho_0 = (P_0 * M_0) / (10.73 * T_r)

4. 更换刚性求解器 #

调用solve_ivp时使用刚性求解器，并设置合理的容差：

sol = solve_ivp(f, W_span, y_0, t_eval=points, method='Radau', rtol=1e-6, atol=1e-8)

5. 调整初始条件 #

精简初始条件为8个变量：

y_0 = np.array([Fet_0, Fac_0, Fo_0, Fp_0, Fa_0, Fc_0, Ft_0, P_0])

完整修复后代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def f(W, y):
    Fet, Fac, Fo, Fp, Fa, Fc, Ft, P = y 

    # 计算分压
    Pet = (Fet / Ft) * P if Ft != 0 else 0
    Pac = (Fac / Ft) * P if Ft != 0 else 0
    Po = (Fo / Ft) * P if Ft != 0 else 0
    Pp = (Fp / Ft) * P if Ft != 0 else 0
    Pa = (Fa / Ft) * P if Ft != 0 else 0
    Pc = (Fc / Ft) * P if Ft != 0 else 0

    r1 = 0.1036*np.exp(-3674/T)*((Po*Pet*Pac*(1+1.7*Pa)))/((1+0.583*Po*(1+1.7*Pa))*(1+6.8*Pac))
    r2 = 1.9365*(10**5)*np.exp(-10116/T)*(Po*(1+0.68*Pa))/(1+0.76*Po*(1+0.68*Pa))

    m = Fet * 28.06 + Fac * 60.052 + Fo * 32 + Fp * 86.09 + Fa * 18.02 + Fc * 44.01
    G = m / A_c                                                                     
    beta = (1.75 * G**2 / (rho_0 * g_c * D_p)) * ((1-phi) / phi**3)

    dFet_dW = -r1 - r2
    dFac_dW = -r1
    dFo_dW  = -r1/2 - 3*r2
    dFp_dW  = r1
    dFa_dW  = r1 + 2*r2
    dFc_dW  = 2*r2
    dFt_dW  = -r1/2
    
    dP_dW   = -(beta/(A_c*(1-phi)*rho_c))*(P_0/P)*(Ft/Ft_0) if P != 0 else 0
    
    return np.array([dFet_dW, dFac_dW, dFo_dW, dFp_dW, dFa_dW, dFc_dW, dFt_dW, dP_dW])

# Integration points
W_span = np.array([0, 100])
points = np.linspace(W_span[0], W_span[1], num = 101)

# Initial conditions
P_0 = 128       # psi
T   = 423       # Kelvin
T_r = T * 9/5   # 转换为°R

Fet_0, Fac_0, Fo_0, Fp_0, Fa_0, Fc_0 = 11, 11, 1.9, 0, 0, 0
Ft_0  = Fet_0 + Fac_0 + Fo_0

y_0 = np.array([Fet_0, Fac_0, Fo_0, Fp_0, Fa_0, Fc_0, Ft_0, P_0])

# Constants
phi = 0.5                                                                       
A_c = 0.0122718463                                                              
rho_c = 60                                                                      
g_c = 115826.4                                                                  
D_p = 0.020833333                                                               
M_0 = (28.05 * Fet_0 + 60.052 * Fac_0 + 31.998 * Fo_0) / (Ft_0)                  
rho_0 = (P_0 * M_0) / (10.73 * T_r)  # 修正后的初始密度计算

# Solution IVP
sol = solve_ivp(f, W_span, y_0, t_eval=points, method='Radau', rtol=1e-6, atol=1e-8)
W   = sol.t

# 提取结果并计算分压
Fet = sol.y[0]
Fac = sol.y[1]
Fo  = sol.y[2]
Fp  = sol.y[3]
Fa  = sol.y[4]
Fc  = sol.y[5]
Ft  = sol.y[6]
P   = sol.y[7]

Pet = (Fet / Ft) * P
Pac = (Fac / Ft) * P
Po  = (Fo / Ft) * P
Pp  = (Fp / Ft) * P
Pa  = (Fa / Ft) * P
Pc  = (Fc / Ft) * P

# 可选：绘制结果
plt.figure(figsize=(12,8))
plt.subplot(2,1,1)
plt.plot(W, Fet, label='Fet')
plt.plot(W, Fac, label='Fac')
plt.plot(W, Fo, label='Fo')
plt.plot(W, Fp, label='Fp')
plt.legend()
plt.xlabel('W')
plt.ylabel('Flow Rates')

plt.subplot(2,1,2)
plt.plot(W, P, label='Total Pressure')
plt.plot(W, Pet, label='Pet')
plt.plot(W, Pac, label='Pac')
plt.plot(W, Po, label='Po')
plt.legend()
plt.xlabel('W')
plt.ylabel('Pressures')
plt.show()

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