Python Forum
values spiral out of control
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
values spiral out of control
#1
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] = 0.375
Cb[0] = 0.00001
Cc[0] = 0
Ct[0] = 0.375
T[0] = 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
C:/Users/teunv/Downloads/Multitubular reactor (1).py:74: RuntimeWarning: overflow encountered in double_scalars
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]
C:/Users/teunv/Downloads/Multitubular reactor (1).py:69: RuntimeWarning: overflow encountered in double_scalars
Ca[n+1] = Ca[n]*dt*-k1[n]+ Ca[n]
C:/Users/teunv/Downloads/Multitubular reactor (1).py:70: RuntimeWarning: overflow encountered in double_scalars
Cb[n+1] = (Ca[n]*k1[n]- k2[n]*Cb[n])*dt + Cb[n]
C:/Users/teunv/Downloads/Multitubular reactor (1).py:71: RuntimeWarning: overflow encountered in double_scalars
Cc[n+1] = Cb[n]*dt*k2[n] + Cc[n]
C:/Users/teunv/Downloads/Multitubular reactor (1).py:74: RuntimeWarning: invalid value encountered in double_scalars
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]
C:/Users/teunv/Downloads/Multitubular reactor (1).py:77: RuntimeWarning: invalid value encountered in double_scalars
Ct[n] = Ca[n]+Cb[n]+ Cc[n]
Reply


Messages In This Thread
values spiral out of control - by Tinfoil - Jan-16-2020, 03:43 PM

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020