values spiral out of control - Printable Version +- Python Forum (https://python-forum.io) +-- Forum: Python Coding (https://python-forum.io/Forum-Python-Coding) +--- Forum: Data Science (https://python-forum.io/Forum-Data-Science) +--- Thread: values spiral out of control (/Thread-values-spiral-out-of-control) values spiral out of control - Tinfoil - Jan-16-2020 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]