1D反应堆建模未运行（Python）

Question

我正在尝试运行此代码，但收到以下错误：

在exp中遇到 KH2_2=1.494np.exp（6025/RT） C:\Users\Ambrosios\OneDrive-uowm.gr\Desktop\draft\xxxx.py:66:RuntimeWarning:exp中遇到溢出 KH2O_4=4.73e-06np.exp（97770/RT） C:\Users\Ambrosios\OneDrive-uowm.gr\Desktop\draft\xxxx.py:69:RuntimeWarning:exp中遇到溢出 KCO_5=7.34e-06np.exp（100395/RT） C:\Users\Ambrosios\OneDrive-uowm.gr\Desktop\draft\xxxx.py:80:RuntimeWarning:在双标量中遇到无效值 r2=（k2KCO2\u2KH2pUCO2pUH2）/（（1+KCO2\u2pUCO2+KH2pUH2）2）
C:\Users\Ambrosios\OneDrive-uowm.gr\Desktop\draft\xxxx.py:83:运行时警告：在双\u标量中遇到无效值 r3=k3KCH4\u3（P\u CH4-P\uh2P\uh2/KP\u3）/
C:\Users\Ambrosios\OneDrive-uowm.gr\Desktop\draft\xxxx.py:86:RuntimeWarning:在双标量中遇到无效值 r4=（（k4/KH2O\u4）（P\uh2o/P\uh2-P\uco/KP\u4））/
C:\Users\Ambrosios\OneDrive-uowm.gr\Desktop\draft\xxxx.py:89:RuntimeWarning:在双标量中遇到无效值 r5=（k5/（KCO_5KCO 2_5）*（P_CO2/P_CO-P_CO/KP_5））/\

有人能帮我吗

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

np.seterr(divide = 'ignore')

Po = 1 # bar
To = 298.15 # K
Fin = 50.3 # mL/min
Fin_mol = Fin/(24150*60*1000) # mol/sec
Fin_m3 = 50.3/(60*1e06) # m3/sec
L = 0.22 # m
D_reactor = 0.008 # m
R_reactor = D_reactor/2 # m
Cross_section_area = np.pi*R_reactor**2 # m
U_heat = 200 # Watt/m^2
Dp = 0.00032 # m
pb = 200 # catalyst bed density (kg/m^3)
omega = 0.0002 #kg
R = 8.314 # J/mol*K
Rg = 0.083 # L*bar/K*mol
u = Fin_m3/Cross_section_area # gas velocity in m/s

dH1 = 247000 # the enthalpy of r1 J/mol 
dH2 = 41700 # the enthalpy of r2 J/mol 
dH3 = 74870 # the enthalpy of r3 J/mol 
dH4 = 131325 # the enthalpy of r4 J/mol
dH5 = 172000 # the enthalpy of r5 J/mol
Cp = 1


molar_ratio = np.array([0.0994,  # CH4
                        0.1033,  # CO2
                        0.7973]) # N2

molar_mass = np.array([16.04,  # CH4
                       44.01,  # CO2
                       28.0134]) # N2

molar_flow = molar_ratio*Fin_mol

pg = (Po/(Rg*To))*np.dot(molar_mass, molar_ratio) # gas density (kg/m^3)


def f(x, z):
    
    P_CH4 = x[0]  
    P_CO2 = x[1]
    P_CO = x[2]
    P_H2 = x[3]
    P_H2O = x[4]
    T = x[5]
    
    k1 = 1.35e-08*np.exp(-3115.22/T)
    k2 = 0.35e06*np.exp(-81030/R*T)
    k3 = 6.95e03*np.exp(-58893/R*T)
    k4 = 5.55e09*np.exp(-166397/R*T)
    k5 = 1.14e15*np.exp(-243835/R*T)
    
    KCO2_1 = 9.25e-08*np.exp(4883.32/T)
    KCH4_1 = 2.46e-07*np.exp(4606.68/T)
    KCO2_2 = 0.5771*np.exp(9262/T)
    KH2_2 = 1.494*np.exp(6025/R*T)
    KCH4_3 = 0.21*np.exp(-567/R*T)
    KH2_3 = 5.18e07*np.exp(-133210/R*T)
    KH2O_4 = 4.73e-06*np.exp(97770/R*T)
    KCH4_4 = 3.49*np.exp(0/R*T)
    KH2_4 = 1.83e13*np.exp(-216145/R*T)
    KCO_5 = 7.34e-06*np.exp(100395/R*T)
    KCO2_5 = 2.81e07*np.exp(-104085/R*T)
    KP_1 = 6.78e14*np.exp(-259660/R*T)
    KP_2 = 56.4971*np.exp(-36580/R*T)
    KP_3 = 2.98e05*np.exp(-84400/R*T)
    KP_4 = 1.3827e07*np.exp(-125916/R*T)
    KP_5 = 1.9393e09*np.exp(-168527/R*T)
    
    r1 = (k1*P_CH4*P_CO2/(KCO2_1*P_CO2 + KCH4_1*P_CH4))*\
         (1 - (P_CO*P_H2)/KP_1*(P_CH4*P_CO2))
    
    r2 = (k2*KCO2_2*KH2_2*P_CO2*P_H2)/((1 + KCO2_2*P_CO2 + KH2_2*P_H2)**2)*\
         (1 - (P_CO*P_H2O)/KP_2*(P_CO2*P_H2))
    
    r3 = k3*KCH4_3*(P_CH4-P_H2*P_H2/KP_3)/\
         (1 + KCH4_3*P_CH4 + 1/KH2_3*P_H2**1.5)**2
          
    r4 = ((k4/KH2O_4)*(P_H2O/P_H2 - P_CO/KP_4))/\
         (1 + KCH4_4*P_CH4 + (1/KH2O_4)*(P_H2O/P_H2) + (1/KH2_4)*P_H2**1.5)**2
        
    r5 = (k5/(KCO_5*KCO2_5)*(P_CO2/P_CO - P_CO/KP_5))/\
         (1 + KCO_5*P_CO + (1/KCO_5*KCO2_5)*(P_CO2/P_CO))**2
        
    dX_CH4dz = pb*omega*L*(r1 + r3)/molar_flow[0]
    dX_CO2dz = pb*omega*L*(r1 + r2 + r3)/molar_flow[1]
    dX_COdz = pb*omega*L*(2*r1 + r2 + r4 + 2*r5)/molar_flow[0]
    dX_H2dz = pb*omega*L*(2*r1 - r2 + 2*r3 + r4)/molar_flow[0]
    dX_H2Odz = pb*omega*L*(r2 - r4)/molar_flow[0]
    
    for Tw in np.linspace(723, 973):
        dTdz = (pb*L/u*pg*Cp)*(r1*dH1 +  r2*dH2 + r3*dH3 + r4*dH4 + r5*dH5) -\
                4*U_heat*L/(Dp*u*pg*Cp)*(T - Tw)
            
    return [dX_CH4dz, dX_CO2dz, dX_COdz, dX_H2dz, dX_H2Odz, dTdz]

# initial conditions
x0 = [0.7, 0.8, 0.5, 0.5, 0.5, 700]

Lspan = np.linspace(0, L, 100)

# call the function
x = odeint(f, x0, Lspan)

p_CH4 = x[:,0] 
p_CO2 = x[:,1]
p_CO = x[:,2]
p_H2 = x[:,3]
p_H2O = x[:,4]
T = x[:,5]