R

#### Renato

##### Guest

Code:

```
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
" Dependant variables : Fet, Fac, Fo, Fp, Fa, Fc, Ft, Pet, Pac, Po, Pp, Pa, Pc, P"
def f(W, y):
Fet, Fac, Fo, Fp, Fa, Fc, Ft, Pet, Pac, Po, Pp, Pa, Pc, P = y
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 # [ lb / min ]
G = m / A_c # [ lb / (min * ft ^2) ]
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)
dPet_dW = (Ft*P*dFet_dW-Fet*P*dFt_dW+Fet*Ft*dP_dW)/Ft**2
dPac_dW = (Ft*P*dFac_dW-Fac*P*dFt_dW+Fac*Ft*dP_dW)/Ft**2
dPo_dW = (Ft*P*dFo_dW -Fo*P*dFt_dW +Fo*Ft*dP_dW)/Ft**2
dPp_dW = (Fp*P*dFp_dW -Fp*P*dFt_dW +Fp*Ft*dP_dW)/Ft**2
dPa_dW = (Fa*P*dFa_dW -Fa*P*dFt_dW +Fa*Ft*dP_dW)/Ft**2
dPc_dW = (Fc*P*dFc_dW -Fc*P*dFc_dW +Fc*Ft*dP_dW)/Ft**2
return np.array([dFet_dW, dFac_dW, dFo_dW, dFp_dW, dFa_dW, dFc_dW, dFt_dW, dPet_dW, dPac_dW, dPo_dW, dPp_dW, dPa_dW, dPc_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
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
Pp_0, Pa_0, Pc_0 = 0, 0, 0
Pet_0 = (Fet_0/Ft_0)*P_0
Pac_0 = (Fac_0/Ft_0)*P_0
Po_0 = (Fo_0/Ft_0)*P_0
P_0 = Pet_0 + Pac_0 + Po_0
y_0 = np.array([Fet_0, Fac_0, Fo_0, Fp_0, Fa_0, Fc_0, Ft_0, Pet_0, Pac_0, Po_0, Pp_0, Pa_0, Pc_0, P_0])
# Constants
phi = 0.5 # Void fraction
A_c = 0.0122718463 # straight section area of the catalyst [ ft ^2 ]
rho_c = 60 # dry catalyst density [ lb / ft ^3 ]
g_c = 115826.4 # gravity conversion factor [ lb * ft / (min ^2 *lbf) ]
D_p = 0.020833333 # average particle diameter [ ft ]
M_0 = (28.05 * Fet_0 + 60.052 * Fac_0 + 31.998 * Fo_0) / (Ft_0) # molar mass of the initial mixture [ lb / lbmol]
rho_0 = (128 * M_0) / (10.73 + 731.07) # Initial density [ lb / ft ^3 ]
# Solution IVP
sol = solve_ivp(f, W_span, y_0, t_eval=points)
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]
Pet = sol.y[7]
Pac = sol.y[8]
Po = sol.y[9]
Pp = sol.y[10]
Pa = sol.y[11]
Pc = sol.y[12]
P = sol.y[13]
```

I tried to write some print() along the code while trying to understand what could be happening and realized that some values like Fet, Fac, Fo, Fp, Fa, Fc, Ft, m, G and beta didn't change. On the other hand, r1, r2, Pet, Pac, Po, Pp, Pa, Pc and P changed cyclically.

<p>I'm trying to solve a problem involving a fixed-bed reactor that uses a system of ordinary differential equations, but for some reason it gets stuck in an infinite loop. Here's the code:</p>

<pre><code>import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

" Dependant variables : Fet, Fac, Fo, Fp, Fa, Fc, Ft, Pet, Pac, Po, Pp, Pa, Pc, P"

def f(W, y):

Fet, Fac, Fo, Fp, Fa, Fc, Ft, Pet, Pac, Po, Pp, Pa, Pc, P = y

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 # [ lb / min ]

G = m / A_c # [ lb / (min * ft ^2) ]

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)

dPet_dW = (Ft*P*dFet_dW-Fet*P*dFt_dW+Fet*Ft*dP_dW)/Ft**2

dPac_dW = (Ft*P*dFac_dW-Fac*P*dFt_dW+Fac*Ft*dP_dW)/Ft**2

dPo_dW = (Ft*P*dFo_dW -Fo*P*dFt_dW +Fo*Ft*dP_dW)/Ft**2

dPp_dW = (Fp*P*dFp_dW -Fp*P*dFt_dW +Fp*Ft*dP_dW)/Ft**2

dPa_dW = (Fa*P*dFa_dW -Fa*P*dFt_dW +Fa*Ft*dP_dW)/Ft**2

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

return np.array([dFet_dW, dFac_dW, dFo_dW, dFp_dW, dFa_dW, dFc_dW, dFt_dW, dPet_dW, dPac_dW, dPo_dW, dPp_dW, dPa_dW, dPc_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

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

Pp_0, Pa_0, Pc_0 = 0, 0, 0

Pet_0 = (Fet_0/Ft_0)*P_0

Pac_0 = (Fac_0/Ft_0)*P_0

Po_0 = (Fo_0/Ft_0)*P_0

P_0 = Pet_0 + Pac_0 + Po_0

y_0 = np.array([Fet_0, Fac_0, Fo_0, Fp_0, Fa_0, Fc_0, Ft_0, Pet_0, Pac_0, Po_0, Pp_0, Pa_0, Pc_0, P_0])

# Constants

phi = 0.5 # Void fraction

A_c = 0.0122718463 # straight section area of the catalyst [ ft ^2 ]

rho_c = 60 # dry catalyst density [ lb / ft ^3 ]

g_c = 115826.4 # gravity conversion factor [ lb * ft / (min ^2 *lbf) ]

D_p = 0.020833333 # average particle diameter [ ft ]

M_0 = (28.05 * Fet_0 + 60.052 * Fac_0 + 31.998 * Fo_0) / (Ft_0) # molar mass of the initial mixture [ lb / lbmol]

rho_0 = (128 * M_0) / (10.73 + 731.07) # Initial density [ lb / ft ^3 ]

# Solution IVP

sol = solve_ivp(f, W_span, y_0, t_eval=points)

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]

Pet = sol.y[7]

Pac = sol.y[8]

Po = sol.y[9]

Pp = sol.y[10]

Pa = sol.y[11]

Pc = sol.y[12]

P = sol.y[13]

</code></pre>

<p>I tried to write some print() along the code while trying to understand what could be happening and realized that some values like Fet, Fac, Fo, Fp, Fa, Fc, Ft, m, G and beta didn't change. On the other hand, r1, r2, Pet, Pac, Po, Pp, Pa, Pc and P changed cyclically.</p>