Ecuaciones diferenciales parciales: reactor químico con reacción de primer orden

Problema
Se tiene un reactor de longitud L, tal que C(x,t) es la concentración de un reactivo A cuya velocidad de reacción es de primer orden. La entrada del reactor está en x = 0, y su salida en x = L, como puede verse en la figura de arriba. Se desea calcular la concentración de A en cualquier instante y posición dentro del reactor, considerando que el reactor se encuentra operando en régimen estacionario y que el reactivo A se distribuye de manera uniforme transversalmente, para una Δx dada. Después, se calcula la concentración para un reactor en régimen no estacionario.

Planteamiento
Se realiza un balance de masa para un elemento de longitud Δx, de la forma:

ec1,
donde V = volumen (m3), Q es la tasa de flujo (m3/h), c es concentración (moles/m3), D es un coeficiente de dispersión (m2/h), Ac es el área transversal del reactor (m2) y k es el coeficiente de reacción de primer orden (h-1). Haciendo Δx → 0 y Δt → 0, la ecuación anterior puede escribirse como
ec2,
donde U = Q/Ac. Aplicando la condición de régimen estacionario, ∂c/∂t = 0, y la ecuación anterior se transforma en
ec3.
Antes de t = 0, el reactor está lleno de agua sin el reactivo A. En t = 0, se inyecta A a una concentración constante cin, de modo que se cumplen las condiciones iniciales:
ec4;
la segunda condición implica que el flujo de salida no se ve afectado por el factor de dispersión.

Solución estacionaria
Tras descomponer los términos de derivadas en componentes de diferencias finitas, la ecuación principal se escribe ahora como
ec5.
Al aplicar las condiciones iniciales para la entrada, la salida y el cuerpo del reactor, y agrupar todas las variables, se encuentran tres ecuaciones en diferencias finitas:

ec6, para el cuerpo del reactor;
ec7, para la entrada, y
ec8 para la salida.
Estas ecuaciones forman un sistema tridiagonal de n ecuaciones. Ya que en este sistema domina la diagonal principal, utilizaremos el método de Jacobi para su solución, auxiliándonos de la librería “matrix.h” (ver aquí).

======================= Código =========================
#include “matrix.h”

using namespace std;

int main(){
//declaracion de variables y parametros
float xo = 0, L = 10, D = 2, dx = 1.25, U = 1, k = 0.2, cin = 100, j = L/dx;
int n = j + 1, col = 0, i = 1;
matrix A(n,n), b(n,1), x(n,1);

//inicializacion de matrices
A = A.zeros(); b = b.zeros();

//coeficientes en la entrada
A(0,0) = (2*D/(U*dx) + k*dx/U + 2 + dx*U/D);
A(0,1) = – 2*D/(U*dx);
b(0,0) = (2 + dx*U/D)*cin;

//coeficientes a lo largo del trayecto
for (int i = 1; i < j; i++){
A(i, col) = -(D/(U*dx) + 0.5);
A(i, col + 1) = 2*D/(U*dx) + k*dx/U;
A(i, col +2) = -(D/(U*dx) – 0.5);
col++;
}

//coeficientes en la salida
A(n – 1, n – 2) = -2*D/(U*dx);
A(n – 1, n – 1) = 2*D/(U*dx) + k*dx/U;

//solucion
cout << "Matriz de coeficientes: " << endl;
A.show();
cout << endl << "Vector solucion: " << endl;
(b.t()).show();
x = A.Jacobi(b);
cout << "Solucion: " << endl;
(x.solshow());
}
Las soluciones para valores de Δx distintos pueden verse a continuación:
conc.
Puede verse que las soluciones convergen entre ellas.

Solución dinámica
Si la solución no es en régimen estacionario, entonces la ecuación en diferencias finitas es:
ec9,
donde los superíndices “l” indican variación en el tiempo y los subíndices “i” variación en el espacio. Aplicando las mismas condiciones iniciales que para el estado estacionario, se encuentran nuevamente tres ecuaciones generales a resolver:

ec10entrada, para la entrada;
ec10salida, para la salida, y
ec10trayecto para el trayecto; nótese que se hicieron los siguientes cambios de variable para simplificar los cálculos:
parametros.

Estas ecuaciones generan un sistema de ecuaciones que, en principio, ya está resuelto: cada valor de la concentración puede calcularse a partir de un valor anterior, por lo que no es necesario plantear un sistema de ecuaciones en forma matricial, basta con realizar los cálculos directamente. El código siguiente realiza esta tarea:

======================== Código ======================
#include “matrix.h”
#include “fstream”
#include “cmath” //quitar comillas al usar el código

using namespace std;

int main(){
//declaracion de variables y parametros
float xo = 0, L = 10, D = 2, dx = 1.25, dt = 0.1, U = 1, k = 0.2, cin = 100, t = 20;
int n = L/dx + 1, m = t/dt + 1;

ofstream file; file.open(“partial.txt”);

//parametros comprimidos
float B = 1-2*D*dt/pow(dx,2) – k*dt;
float E = (dt/dx)*(D/dx – U/2);
float F = (dt/dx)*(D/dx + U/2);
float G = 2*dx*U/D;
float H = 2*D*dt/pow(dx,2);

matrix C(n,m);

//inicializacion de matrices
C = C.zeros();

file << "x \t " << "t \t " << "C(x,t) \t " << endl << endl;

//desarrollo
for (int l = 0; l < m; l++){
for (int i = 0; i < n ; i++){
if (i == 0)
C(0,l+1) = C(0,l)*(B – F*G) + C(1,l)*(E + F) + F*G*cin;
else if(i == n – 1)
C(i,l+1) = C(i,l)*B + C(i-1,l)*H;
else
C(i,l+1) = C(i,l)*B + C(i+1,l)*E + C(i-1,l)*F;

file << i*dx << " \t " << l*dt << " \t " << C(i,l) << " \t " << endl;
}
file << endl;
}

C.show();
}

Nota: este código guarda los resultados en un archivo de texto para facilitar su análisis.

Los resultados pueden verse en la gráfica siguiente:
concparcial

Debe observarse que, para tiempos pequeños (menores a un segundo, aproximadamente), la concentración es cero en algunas regiones del reactor alejadas de la entrada. Conforme el tiempo pasa, el reactivo ha llegado más lejos en el reactor, por lo que la concentración es distinta de cero. Cuando el tiempo es grande (t = 20), la solución se aproxima a las soluciones obtenidas para el régimen estacionario, lo cual es de esperarse cuando se trata de un sistema con reacción química. Este comportamiento se observa abajo:
i win.png

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s