# Python Forum

Full Version: values spiral out of control
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Hi users,
for a project at my university i have to optimize a reactor. To find out the concentrations of the 3 phases in the reactor tube we have written this python script. It is an euler integration of 4 reactions and the last one is based on an energy balance. I would love it if one of you could point out a mistake i am making that leads to the results i have added below the code.

```
import numpy as np
import matplotlib.pyplot as plt

#defining the variables
z = 2 #length of tube
t0 = 0
dz = 0.01 #size of intergration nodes
uc = 0.2 #coolant lineair velocity
us = 0.5 #superficial velocity of reacting gas
dt = dz/us #amount of time per integration node
k01 = 3.62E9
k02 = 7.99E5
Ca_0 = 0.375
rho = 1791 #coolant density

t1 = z/uc #using time to calculate with
dH1 = 1.54E6 #heat effect van eerste reactie J/mol
dH2 = 1.06E6 #heat effect van de tweede reactie
d = 0.04 #inside diameter of the tube
Tc = 450
Ea1 = 1.34E5#activation energy of the first reaction
R = 8.314 #J / mol·K. ideal gas constant
Ea2 = 8.80E4#activation energy of the second reaction
Ti = 450
U_m = 13.333E3
dout= 0.003 + d

t = np.linspace(t0, t1, int(1 + (t1-t0)/dt))

Ca = np.zeros(len(t))
Cb = np.zeros(len(t))
Cc = np.zeros(len(t))
Ct = np.zeros(len(t))
T = np.zeros(len(t))
k1 = np.zeros(len(t))
k2 = np.zeros(len(t))

Ca = 0.375
Cb = 0.00001
Cc = 0
Ct = 0.375
T = Ti

#calculating U_out for triangular and squared
Re_c = 1.8853E5 * dout
U_out_t = 0.2626* ((Re_c**0.6)/dout)
U_out_sq = 0.1968*((Re_c**0.6)/dout)

#calculating U_in
ef = 352.98/Ti #density
Re_p = (96.774*ef)*us
U_in = (0.1020*(Re_p**0.8))/(d*np.exp(0.018/d))
U = U_in + U_m + U_out_t

for n in range(len(t)-1):

#post-exponential factors reaction rate constants
k1[n] = k01*np.exp(-Ea1/(R*T[n]))
k2[n] = k02*np.exp(-Ea2/(R*T[n]))

#concentrations
Ca[n+1] = Ca[n]*dt*-k1[n]+ Ca[n]
Cb[n+1] = (Ca[n]*k1[n]- k2[n]*Cb[n])*dt + Cb[n]
Cc[n+1] = Cb[n]*dt*k2[n] + Cc[n]

#calculating the temperature in the tube with energy balance
T[n+1] = dt*(k1[n]*Ca[n]*dH1 + (k1[n]*Ca[n]-k2[n]*Cb[n])*dH2 - (4/d)*U*(T[n]-Tc))/(rho*Cb[n])+T[n]

#the total concentration to check it stays constant
Ct[n] = Ca[n]+Cb[n]+ Cc[n]

print("The inlet temperature used:", Ti)
print("The U that was used:", U)

plt.figure(1)
plt.plot(t, Ca, "g--")
plt.plot(t, Cb, "y--")
plt.plot(t, Cc, 'r--')

plt.figure(2)
plt.plot(t, Ct)

plt.figure(3)
plt.plot(t, T)
```
Results:

The inlet temperature used: 450
The U that was used: 14715.315745904882