import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint # define CSTR model def cstr(x,t,u,Tf,Caf): # Inputs (3): # Temperature of cooling jacket (K) Tc = u # Tf = Feed Temperature (K) # Caf = Feed Concentration (mol/m^3) # States (2): # Concentration of A in CSTR (mol/m^3) Ca = x[0] # Temperature in CSTR (K) T = x[1] # Parameters: # Volumetric Flowrate (m^3/sec) q = 100 # Volume of CSTR (m^3) V = 100 # Density of A-B Mixture (kg/m^3) rho = 1000 # Heat capacity of A-B Mixture (J/kg-K) Cp = 0.239 # Heat of reaction for A->B (J/mol) mdelH = 5e4 # E - Activation energy in the Arrhenius Equation (J/mol) # R - Universal Gas Constant = 8.31451 J/mol-K EoverR = 8750 # Pre-exponential factor (1/sec) k0 = 7.2e10 # U - Overall Heat Transfer Coefficient (W/m^2-K) # A - Area - this value is specific for the U calculation (m^2) UA = 5e4 # reaction rate rA = k0*np.exp(-EoverR/T)*Ca # Calculate concentration derivative dCadt = q/V*(Caf - Ca) - rA # Calculate temperature derivative dTdt = q/V*(Tf - T) \ + mdelH/(rho*Cp)*rA \ + UA/V/rho/Cp*(Tc-T) # Return xdot: xdot = np.zeros(2) xdot[0] = dCadt xdot[1] = dTdt return xdot # Steady State Initial Conditions for the States Ca_ss = 0.87725294608097 T_ss = 324.475443431599 x0 = np.empty(2) x0[0] = Ca_ss x0[1] = T_ss # Steady State Initial Condition u_ss = 300.0 # Feed Temperature (K) Tf = 350 # Feed Concentration (mol/m^3) Caf = 1 # Time Interval (min) t = np.linspace(0,25,251) # Store results for plotting Ca = np.ones(len(t)) * Ca_ss T = np.ones(len(t)) * T_ss u = np.ones(len(t)) * u_ss # Step cooling temperature to 295 u[10:100] = 303.0 u[100:190] = 297.0 u[190:] = 300.0 # Simulate CSTR for i in range(len(t)-1): ts = [t[i],t[i+1]] y = odeint(cstr,x0,ts,args=(u[i+1],Tf,Caf)) Ca[i+1] = y[-1][0] T[i+1] = y[-1][1] x0[0] = Ca[i+1] x0[1] = T[i+1] # Construct results and save data file # Column 1 = time # Column 2 = cooling temperature # Column 3 = reactor temperature data = np.vstack((t,u,T)) # vertical stack data = data.T # transpose data np.savetxt('data_doublet.txt',data,delimiter=',',\ header='Time,Tc,T',comments='') # Plot the results plt.figure() plt.subplot(3,1,1) plt.plot(t,u,'b--',lw=3) plt.ylabel('Cooling T (K)') plt.legend(['Jacket Temperature'],loc='best') plt.subplot(3,1,2) plt.plot(t,Ca,'r-',lw=3) plt.ylabel('Ca (mol/L)') plt.legend(['Reactor Concentration'],loc='best') plt.subplot(3,1,3) plt.plot(t,T,'k.-',lw=3) plt.ylabel('T (K)') plt.xlabel('Time (min)') plt.legend(['Reactor Temperature'],loc='best') plt.show()