PRACTICAS DE SIMULACION DE PROCESOS
PRACTICA N° 3
SIMULACIÓN DE REACTOR TUBULAR
1. MARCO TEORICO DE DIFUSION
El objetivo aclarar algunos aspectos relacionados con el modelamiento y simulación
de un reactor tubular para un fluido incompresible y de fase homogénea, en el que
se tomen en cuenta los fenómenos que pueden influenciar notablemente su
rendimiento.
En lugar de que fluya fluido por una tubería, suponga que la tubería es un reactor
tubular en el que tiene lugar la reacción A B. A medida que una rodaja de
material se mueve a lo largo del reactor, la concentración de CA reaccionante
disminuye a medida que se consume A. La densidad , la velocidad V y la
concentración CA pueden variar con el tiempo y la posición axial z. Seguimos
asumiendo condiciones de flujo homogéneo de modo que no haya gradientes
radiales en velocidad, densidad o concentración.
La concentración de A alimentada a la entrada del reactor en z = 0 se define como
𝐶𝐴(𝑡,0) = 𝐶𝐴0(𝑡)
La concentración de A en el fluido del reactor en z = L se define como
𝐶𝐴(𝑡,𝐿) = 𝐶𝐴𝐿(𝑡)
Ahora queremos aplicar la ecuación de continuidad del componente para el reactivo
A a un pequeño segmento diferencial de ancho dz, como se muestra en la figura.
Los términos de entrada se pueden dividir en dos tipos: flujo masivo y difusión. La
difusión puede ocurrir debido al gradiente de concentración en la dirección axial.
Por lo general, es mucho menos importante que el flujo masivo en la mayoría de
los sistemas prácticos, pero lo incluimos aquí para ver qué aporta al modelo.
Diremos que el flujo difusivo de A, NA, (moles de A por unidad de tiempo por unidad
de área), está dado por el tipo de relación de la ley de Fick.
Que después de aplicar un balance de materia nos dará la siguiente ecuación o
modelo matemático que gobierna el sistema de Reactor Tubular, que es lel
siguiente:
𝜕𝐶𝐴 𝜕(𝑣𝐶𝐴 ) 𝜕 𝜕𝐶𝐴
+ + 𝑘𝐶𝐴 = (𝐷𝐴 )
𝜕𝑡 𝜕𝑧 𝜕𝑧 𝜕𝑧
Que despejado queda como sigue:
1
DR. ROLANDO S. BASURCO CARPIO
𝜕𝐶𝐴 𝜕 2 𝐶𝐴 𝜕(𝑣𝐶𝐴 )
= (𝐷𝐴 ) − − 𝑘𝐶𝐴
𝜕𝑡 𝜕𝑧 2 𝜕𝑧
Suponiendo que la densidad como la velocidad son constantes.
𝜕𝐶𝐴 𝜕 2 𝐶𝐴 𝜕(𝐶𝐴 )
= (𝐷𝐴 ) − 𝑣 − 𝑘𝐶𝐴
𝜕𝑡 𝜕𝑧 2 𝜕𝑧
2. MÉTODO DE RESOLUCIÓN
Para la solución de las EDP pueden utilizarse una gran variedad de técnicas. Entre
otras técnicas pueden mencionarse, los métodos de diferencias finitas, los métodos
de elementos finitos, los métodos de volúmenes finitos, los métodos espectrales,
etc. Por ser el más sencillo, y no requerir transformar la EDP, se analizará el método
de las diferencias finitas (MDF). Consiste en aproximar la derivada parcial por
cocientes incrementales.
Para este caso usaremos el método de Crank Nicholson Modificado
El método de diferencias finitas se basa en las fórmulas para aproximar las
derivadas primera y segunda de una función.
ui , j +1 − ui , j
ut
h
ux
(ui −1, j +1 − ui , j +1 ) θ + (ui −1, j − ui , j ) (1 − θ )
u xx
(ui −1, j +1 − 2ui , j +1 + ui +1, j +1 ) θ + (ui −1, j − 2ui , j + ui +1, j ) (1 − θ )
2
k
1. Discretización del sistema
2. Aproximación de las ecuaciones diferenciales por diferencias finitas
Para i = 2 a i =N-1
c ci , j +1 − ci , j
=
t h
2c
= D 2 (ci −1, j − 2ci , j + ci +1, j )(1 − ) + (ci −1, j +1 − 2ci , j +1 + ci +1, j +1 )
1
D
x 2
x
2
PRACTICAS DE SIMULACION DE PROCESOS
𝜕𝑐 1
𝑢 =𝑢 [(𝑐 − 𝑐𝑖,𝑗 )(1 − 𝜃) + (𝑐𝑖+1,𝑗+1 − 𝑐𝑖,𝑗+1 )𝜃] − 𝑘𝑐𝑖,𝑗
𝜕𝑥 𝛥𝑥 𝑖+1,𝑗
lado derecho de la igualdad
𝜌(1 + 𝜃)𝑐𝑖−1,𝑗 + [1 − (1 − 𝜃)(2𝜌 + 𝑃𝑒) − 𝑘] ⋅ 𝑐𝑖,𝑗 + [(1 − 𝜃)(𝜌 − 𝑃𝑒) ⋅ 𝑐𝑖+1,𝑗
lado izquierdo de la igualdad
−𝜌𝜃 ⋅ 𝑐𝑖−1,𝑗+1 + [1 + 𝜃(2𝜌 − 𝑃𝑒)] ⋅ 𝑐𝑖,𝑗+1 − 𝜃(𝜌 − 𝑃𝑒) ⋅ 𝑐𝑖+1,𝑗+1
Para i = 1
c c1, j +1 − c1, j
=
t h
2c
D 2 = D 2 (2c0 − 3c1, j + c2, j )(1 − ) + (2c0 − 3c1, j +1 + c2, j +1 )
1
x x
𝜕𝑐 2
𝑢 =𝑢 [(𝑐 − 𝑐0 )(1 − 𝜃) + (𝑐1,𝑗+1 − 𝑐0 )𝜃] − 𝑘𝑐0,𝑗
𝜕𝑥 𝛥𝑥 1,𝑗
lado derecho de la igualdad
2𝜌(1 + 𝑃𝑒) ⋅ 𝑐0 + [1 − 𝜌(1 − 𝜃)(3 + 2𝑃𝑒)] ⋅ 𝑐1,𝑗 + 𝜌(1 − 𝜃) ⋅ 𝑐2,𝑗 − 𝑘𝑐0,𝑗
lado izquierdo de la igualdad
[1 + 𝜌𝜃(3 + 2𝑃𝑒)] ⋅ 𝑐1,𝑗+1 − 𝜌𝜃 ⋅ 𝑐2,𝑗+1
Para i = N
c cN , j +1 − cN , j
=
t h
𝜕 2𝑐 1
𝐷 2 = 𝐷 2 [(𝑐𝑁−1,𝑗 − 3𝑐𝑁,𝑗 + 2𝑐𝑁+1,𝑗 )(1 − 𝜃) + (𝑐𝑁−1,𝑗+1 − 3𝑐𝑁,𝑗+1 + 2𝑐𝑁+1,𝑗+1 )𝜃]
𝜕𝑥 𝛥𝑥
𝜕𝑐 2
𝑢 =𝑢 [(𝑐 − 𝑐𝑁+1,𝑗 )(1 − 𝜃) + (𝑐𝑁,𝑗+1 − 𝑐𝑁+1,𝑗+1 )𝜃] − 𝑘𝑐𝑁,𝑗
𝜕𝑥 𝛥𝑥 𝑁,𝑗
lado derecho de la igualdad
𝑢 ⋅ 𝛥𝑥
𝜌(1 − 𝜃) ⋅ 𝑐𝑁−1,𝑗 + [1 − 𝜌(1 − 𝜃)(3 + 2𝑃𝑒) − 𝑘] ⋅ 𝑐𝑁,𝑗 + 2𝜌(1 − 𝜃) (1 + ) ⋅ 𝑐𝑁+1,𝑗
𝐷
lado izquierdo de la igualdad
Dependiendo de la condición de contorno en x = finito o infinito se tiene:
3
DR. ROLANDO S. BASURCO CARPIO
Caso 1: Para la condición de borde c =0 en x = L se tiene que cN+1,j = 0. Luego:
lado derecho de la igualdad
𝑢 ⋅ 𝛥𝑥
𝜌(1 − 𝜃) ⋅ 𝑐𝑁−1,𝑗 + [1 − 𝜌(1 − 𝜃) (3 + 2 ) − 𝑘] ⋅ 𝑐𝑁,𝑗
𝐷
lado izquierdo de la igualdad
𝑢 ⋅ 𝛥𝑥
−𝜌𝜃 ⋅ 𝑐𝑁−1,𝑗+1 + [1 + 𝜌𝜃 (3 + 2 )] ⋅ 𝑐𝑁,𝑗+1
𝐷
Caso 2: Para la condición de borde c =0 en x = ∞ esta es reemplazada por La
condición c = 0 en x = L , luego cN+1,j = cN,j . Por tanto:
x
lado derecho de la igualdad
𝜌(1 − 𝜃) ⋅ 𝑐𝑁−1,𝑗 + [1 − 𝜌(1 − 𝜃) − 𝑘] ⋅ 𝑐𝑁,𝑗
lado izquierdo de la igualdad
−𝜌𝜃 ⋅ 𝑐𝑁−1,𝑗+1 + (1 + 𝜌𝜃) ⋅ 𝑐𝑁,𝑗+1
A= =
.
{1 + 𝜃(3𝜌 + 2𝑃𝑒)} −𝜌𝜃
−𝜌𝜃 {1 + 𝜃(2𝜌 − 𝑃𝑒)} −𝜃(𝜌 − 𝑃𝑒)
−𝜌𝜃 {1 + 𝜃(2𝜌 − 𝑃𝑒)} −𝜃(𝜌 − 𝑃𝑒)
.... .... ....
−𝜌𝜃 {1 + 𝜃(2𝜌 − 𝑃𝑒)} −𝜃(𝜌 − 𝑃𝑒)
[ −𝜌𝜃 {1 + 𝜃(3𝜌 + 2𝑃𝑒)}]
B= =
{1 − (1 − 𝜃)(3𝜌 + 2𝑃𝑒)} 𝜌(1 − 𝜃)
𝜌(1 − 𝜃) {1 − (1 − 𝜃)(2𝜌 − 𝑃𝑒) − 𝑘} (1 − 𝜃)(𝜌 − 𝑃𝑒)
𝜌(1 − 𝜃) {1 − (1 − 𝜃)(2𝜌 − 𝑃𝑒) − 𝑘} (1 − 𝜃)(𝜌 − 𝑃𝑒)
.... .... ....
𝜌(1 − 𝜃) {1 − (1 − 𝜃)(2𝜌 − 𝑃𝑒) − 𝑘} (1 − 𝜃)(𝜌 − 𝑃𝑒)
[ 𝜌(1 − 𝜃) {1 − (1 − 𝜃)(3𝜌 + 2𝑃𝑒) − 𝑘)}]
(2(𝜌 + 𝑃𝑒) − 𝑘)𝑐0
0
D= 0
:
[ 0 ]
Luego el sistema de ecuaciones a resolver:
A c j+1 = B c j + D
3. EJEMPLO A RESOLVER
La reacción de primer orden
A→B
4
PRACTICAS DE SIMULACION DE PROCESOS
se lleva a cabo en un PFR de 10 cm de diámetro y 6.0 m de largo. La velocidad de
reacción específica es de 0.25 min-1 = 0.0042 s-1.
Datos adicionales:
CA0 = 0.5 mol/dm3 = 500 mol/m3, U0 = 1.24 m/min= 2.1*10-2 m/s, Da = 1.05 m2/min.
DAB = 7.6*10-5 m2/min =1.3*10-7 m2/s.
Realizar la simulación
Programa: ReacTubular
clear all
% Transporte de Soluto Horizontal de un Reactor Tubular C=Co a C=0
% Definicion de parametros
De=1.3*10^-7;
u=2.1*10^-2;
Czero=500;
Long=6.0;
Tfinal=300;
delt=3;
Ncel=20;
delx=Long/Ncel;
Nstep=Tfinal/delt;
rho=De*delt/delx^2;
Adu=u*delt/delx;
theta=0.5;
kreaccion= 0.0042;
% formar matrices y vectores
%Matriz A
A(1,1)=1 + theta*(3*rho+2*Adu);
A(1,2)=-rho*theta;
A(Ncel,Ncel-1)=-rho*theta;
A(Ncel,Ncel)=1+ theta*(3*rho+2*Adu);
for i=2:Ncel-1
A(i,i-1)=-rho*theta;
A(i,i)=1+ theta*(2*rho-Adu);
A(i,i+1)=-theta*(rho-Adu);
end
% Matriz B
B(1,1)=1 - (1-theta)*(3*rho+2*Adu);
B(1,2)=rho*(1-theta);
B(Ncel,Ncel-1)=rho*(1-theta);
B(Ncel,Ncel)=1+ (1-theta)*(3 *rho +2*Adu)-kreaccion;
for i=2:Ncel-1
B(i,i-1)=rho*(1-theta);
B(i,i)=1-(1-theta)*(2*rho-Adu)-kreaccion;
B(i,i+1)=(1-theta)*(rho-Adu);
dist(i)=delx*(i-0.5);
end
% Vector D
D(1)=(2*(rho+ Adu)-kreaccion)*Czero;
D(2:Ncel)=0;
C(1:Ncel)=0;
5
DR. ROLANDO S. BASURCO CARPIO
dist(1)=delx/2;
dist(Ncel)=delx*(Ncel-0.5);
% Calculo de Concentraciones
D=D';
C=C';
dist=dist'
for k=1:Nstep
Z=B*C+D;
C=A\Z;
tiempo(k)=k*delt;
Conc(:,k)=C/Czero;
end
plot(dist,Conc(:,1:1000:Nstep))
title('Conc v/s. distancia')
xlabel('Distancia,m')
ylabel('Concentracion, C/Co')
Resultados
Conc v/s. distancia
1
0.9
0.8
0.7
Concentracion, C/Co
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Distancia,m
4. TAREA DE LA PRÁCTICA
1.- Resolver el algebra de los polinomios de la celda 1, celda 2 a N-1 y de la celda
N.
2.- Investigar cómo se determina el tiempo de residencia en un reactor tubular y
generar la simulación en el programa que desee (Excel, Matlab o Visual Basic)