In [12]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from IPython.html.widgets import interact
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
plt.style.use('bmh')
$ E-S+G = A $
$ -dn_A+r_A'dM_{cat} = 0 $
$ dn_A = -n_{A0}dX_A $
La velocidad de reacción es de orden 2:
$ r_A' = -aKCa^2 $
$ Ca = \frac{n_A}{Qv} = \frac{n_{A0}(1-X_A)}{Qv_0\frac{P_0T}{PT_0}(1+\epsilon Xa)} = \frac{C_{A0}(1-X_A)}{\frac{P_0T}{PT_0}(1+\epsilon Xa)} $
Se comporta como un gas ideal, por tanto:
$ C_{A0} = \frac{P_0}{RT_0} $
La constante de reacción varía con la ecuación de Arrhenius:
$ K = k_0 e^\frac{-E}{RT} $
La desactivación varía con el tiempo:
$ \frac{da}{dt} = -k_da $
$ a = e^{-k_dt} $
$ M_{cat} = V_{cat}\rho _{cat} = V_{lecho}(1-\epsilon _{lecho})\rho _{cat} = SL(1-\epsilon _{lecho})\rho _{cat} $
$ dM_{cat} = S(1-\epsilon _{lecho})\rho _{cat}dL $
$\frac{dM_{cat}}{dL} = S(1-\epsilon _{lecho})\rho _{cat} $
Teniendo en cuenta la difusión interna:
$ \eta = \frac{tagh(mL)}{mL} $
$ L = \frac{R}{3} $
$ m = ({\frac{KaC_a3}{2De}})^{0.5} $
Sustituyendo terminos en:
$ -dn_A+r_A'dM_{cat} = 0 $
$ n_{A0} = Q_{v0}C_{A0} = v_0SC_{A0} $
Obtenemos:
$ v_0SC_{A0}dX_A-\eta e^{-k_dt}k_0 e^\frac{-E}{RT} (\frac{n_{A0}(1-X_A)}{Qv_0\frac{P_0T}{PT_0}(1+\epsilon Xa)})^2 S(1-\epsilon _{poro})\rho _{cat}dL = 0 $
Reorganizarndo:
$ \frac{dX_A}{dL} =\frac{\eta a K Ca^2 S(1-\epsilon _{poro})\rho _{cat}}{v_0SC_{A0}} = \frac{\eta e^{-k_dt}k_0 e^\frac{-E}{RT} (n_{A0}(1-X_A))^2 S(1-\epsilon _{poro})\rho _{cat}}{(Qv_0\frac{P_0T}{PT_0}(1+\epsilon Xa))^2v_0SC_{A0}} = \frac{\eta e^{-k_dt}k_0 e^\frac{-E}{RT} (C_{A0}(1-X_A))^2 (1-\epsilon _{poro})\rho _{cat}}{(\frac{P_0T}{PT_0}(1+\epsilon Xa))^2v_0C_{A0}} $
$ d(\sum n_jh_j) -dQ = 0 $
Al ser adiabatico:
$ d(\sum n_jh_j)= 0 $
Y $ \Delta C_p = 0 $:
$ dT(\sum n_{j0}C_p)+\Delta H_rdX_A = 0 $
$ \frac{dT}{dX_A} = \frac{-\Delta H_r}{C_{pA}} $
Usando la regla de la cadena:
$ \frac{dT}{dL} = \frac{dT}{dX_A}\frac{dX_A}{dL} $
E = 7.20e4 J/mol
$k_0$ = 7.7e5
$T_0$ = 250ºC
$\Delta H_r$ = -800 J/mol
$C_{pA}$ = 15 J/(mol·K)
$v_0$ = 3 m/s
P = 5.0e5 Pa
R = 8.31 J/(mol·K)
$\epsilon$ = 2
$\epsilon_{lecho}$ = 0.4
$\rho_{cat} = 2000$ $kg/(m^3)$
$R_{cat} = \frac{6·10^{-3}}{2}m$
$D_e$ = $2.66·10^{-8} m^2/s$
In [13]:
class reactor():
def __init__(self,Xa0,T0,kd,t,adiabatico):
self.adi = adiabatico
self.Xa0 = Xa0
self.T0 = T0
self.kd = kd
self.t = t
self.E = 7.20e4
self.k0 = 7.7e5
self.Hr = -800
self.Cp = 15
self.v0 = 3
self.P = 5e5
self.R = 8.31
self.poro_lecho =0.4
self.rho_cat = 2000.
self.R_cat = 6e-3/2
self.De =2.66e-8
self.expansion = 2.
def actividad(self):
return np.exp(-self.kd*self.t)
def concentracion_ini_A(self):
return self.P/(self.R*self.T0)
def concentracion_A(self, Xa, T):
a = self.concentracion_ini_A()*(1-Xa)
b = (1+self.expansion*Xa)*(T/self.T0)
return a/b
def constante_reaccion(self, T):
return self.k0*np.exp(-self.E/(self.R*T))
def mL(self, Xa, T):
L = self.R_cat/3
a = (self.constante_reaccion(T)*self.actividad()*self.concentracion_A(Xa,T)*3)
b = 2*self.De
m = (a/b)**0.5
return m*L
def efectividad(self, Xa, T):
return np.tanh(self.mL(Xa,T))/self.mL(Xa,T)
def ode(self, y, l):
if self.adi:
Xa = y[0]
T = y[1]
#######BM#######
a = self.efectividad(Xa,T)*self.constante_reaccion(T)*self.actividad()
b = self.concentracion_A(Xa, T)**2*(1-self.poro_lecho)*self.rho_cat
c = self.v0*self.concentracion_ini_A()
dXadL= (a*b)/c
#######BE#######
dTdXa = -self.Hr/self.Cp
dTdL = dTdXa*dXadL
return [dXadL, dTdL]
else:
Xa = y[0]
#######BM#######
a = self.efectividad(Xa,self.T0)*self.constante_reaccion(self.T0)*self.actividad()
b = self.concentracion_A(Xa, self.T0)**2*(1-self.poro_lecho)*self.rho_cat
c = self.v0*self.concentracion_ini_A()
dXadL= (a*b)/c
return [dXadL,0]
In [14]:
def RFP_reactor(kd, tmax, Lmax, N, Adiabatico):
nt = 50
nl = 50
t = np.linspace(0,tmax,nt)
l = np.linspace(0,Lmax,nl)
Xa0 = 0
T0 = 273+250
Xa = np.zeros((nt,nl))
T = np.zeros((nt,nl))
Ca = np.zeros((nt,nl))
Kr = np.zeros((nt,nl))
r = np.zeros((nt,nl))
actividad = np.zeros(nt)
efecti = np.zeros((nt,nl))
for j in range(nt):
reac = reactor(Xa0,T0,kd,t[j],Adiabatico)
y = odeint(reac.ode,[Xa0, T0], l)
Xa[j,:] = y[:,0]
T[j,:] = y[:,1]
actividad[j] = reac.actividad()
Ca[j,:] = reac.concentracion_A(y[:,0],y[:,1])
Kr[j,:] = reac.constante_reaccion(y[:,1])
efecti[j,:] = reac.efectividad(y[:,0],y[:,1])
r[j,:] = efecti[j,:]*actividad[j]*Kr[j,:]*(Ca[j,:]*Ca[j,:])
if N==0:
####################
#Representar Xaf vs t
plt.plot(t,Xa[:,-1])
plt.ylabel('Xa')
plt.xlabel('tiempo / s')
plt.title('Conversión de A en la posición final del reactor \n'+
'frente al tiempo')
elif N==1:
####################################
#Representar Tf frente al tiempo
plt.plot(t,T[:,-1])
plt.ylabel('Tf / K')
plt.xlabel('tiempo / s')
plt.title('Temperatura en la posición final del reactor \n'+
'frente al tiempo')
elif N==2:
####################################
#Representar velocidad de reaccion frente al tiempo y longitud
ll, tt = np.meshgrid(l,t)
ax = plt.axes(projection='3d')
ax.plot_wireframe(tt, ll, r)
plt.ylabel('L / m')
plt.xlabel('tiempo / s')
plt.title('velocidad de reacción de A \n'+
'frente al tiempo y longitud')
plt.figure()
plt.contourf(tt, ll, r)
plt.ylabel('L / m')
plt.xlabel('tiempo / s')
plt.title('velocidad de reacción de A \n'+
'frente al tiempo y longitud')
plt.colorbar()
elif N==3:
####################################
#Representar Ca frente al tiempo y longitud
ll, tt = np.meshgrid(l,t)
ax = plt.axes(projection='3d')
ax.plot_wireframe(tt, ll, Ca)
plt.ylabel('L / m')
plt.xlabel('tiempo / s')
plt.title('Concentración de A \n'+
'frente al tiempo y longitud')
plt.figure()
plt.contourf(tt, ll, Ca)
plt.ylabel('L / m')
plt.xlabel('tiempo / s')
plt.title('Concentración de A \n'+
'frente al tiempo y longitud')
plt.colorbar()
elif N==4:
####################################
#efectiviad frente al tiempo y longitud
ll, tt = np.meshgrid(l,t)
ax = plt.axes(projection='3d')
ax.plot_wireframe(tt, ll, efecti)
plt.ylabel('L / m')
plt.xlabel('tiempo / s')
plt.title('Factor deeficacia frente al tiempo y longitud')
plt.figure()
plt.contourf(tt, ll, efecti)
plt.ylabel('L / m')
plt.xlabel('tiempo / s')
plt.title('Factor deeficacia frente al tiempo y longitud')
plt.colorbar()
In [15]:
interact(RFP_reactor, kd=[0.01,1,0.01],
tmax = [10,1000,10], Lmax = [0.001,0.1,0.001],
N = {"Xaf vs t":0, "Tf vs t":1, "Velocidad de reaccion":2,
"Concentracion de A":3, "Factor eficacia":4},
Adiabatico = True)
Out[15]:
In [ ]: