Diferencia entre revisiones de «Reaccion Autocatalisis Grupo 2B»

Revisión actual del 10:50 28 abr 2017

1 Introducción

Una reacción de autocatálisis es un proceso mediante el cual un compuesto químico induce y controla una reacción química sobre sí mismo. En este caso, se estudiará una reacción regida por la constante k, bimolecular e irreversible, en la cual se suponen unos reactivos A y B que presentan un volumen y una temperatura constantes:

                           				A+B → C

Dichas condiciones satisfacen la Ley de acción de masas, justifica que "La velocidad de la reacción en cuestión será proporcional al producto de las concentraciones de A y B". Para desarrollar el análisis, es necesario apoyarse en el Principio de conservación de la materia, por el cual "la suma de la masa de los reactivos y la de los productos se mantendrá constante a lo largo del tiempo". Así pues, se obtendrá un tercer compuesto que actuará como producto:

                                                      A+B → k * 2B

La resolución de la ecuación diferencial se apoya en las dos leyes mencionadas con anterioridad. De esta forma, es posible determinar las siguientes ecuaciones: [math] \left\{\begin{matrix} x=A+B\\ y=2*B \end{matrix}\right. [/math]

Por el método de sustitución se deducen las variables dependientes A y B: [math] \left\{\begin{matrix} A=x-\frac{y}{2}\\ B=\frac{y}{2} \end{matrix}\right. [/math]

En base al Principio de conservación de la masa, se conoce la igualdad:

                                                      x+y=cte                                                                       (1)

A su vez, también se sabe que la velocidad es proporcional al producto de las concentraciones de los reactivos: v=k*A*B, por lo tanto, v=y’=-x’. De esta forma se concluye que y’=k*A*B en la que se sustituyen los valores de los reactivos previamente determinados:

    y’=k*[x-y/2]*[y/2] → y'=k*[x-y/2]*y → y'=k*[cte-y-y/2]*y → y'=k*[cte-3/2*y]*y → y'=3/2*k*[2/3*cte-y]*y → y'=k*[cte-y]*y → y'=k*x*y            

De (1) se puede afirmar que:

                                                     x=cte-y

Y sustituyendo en la fórmula anteriormente demostrada, se obtiene que:

                                                 y'=-k*y^2+k*y*cte

A partir de la anterior ecuación diferencial se puede formar el PVI asociado: [math] \left\{\begin{matrix} y'=-k*y^2+k*y*cte\\ y(0)=y_0 \end{matrix}\right. [/math]

Para poder resolver el PVI se necesita calcular el valor de la cte:

                                                 x=A+B → A=x-y/2
                                                   y=2B → B=y/2

A partir de los datos numéricos, A=2 y B=0.05 mol/L, y despejando en las ecuaciones anteriores, se obtienen los valores de x e y, siendo 1.95 y 0.1 respectivamente. Por lo tanto:

                                              cte=x+y=1.95+0.1=2.05

Para comprobar que el PVI tiene solución única se aplica el Teorema de Picard, siendo:

                                              f=-k*y^2+k*y*cte

Al ser f una función polinómica está definida y es continua en RxR. A continuación, se calcula la derivada parcial de f con respecto a y, dando:

                                              f'=-2*k*y+k*cte

Análogamente, al ser la derivada una función polinómica, está definida y es continua en RxR, concluyendo, así, que el PVI tiene solución y es única.

2 Resolución del PVI

La resolución del PVI se ha llevado acabo por diferentes métodos. Para ello, se ha escogido un paso h=0.1 y una k=1.4 mol/s. El estudio se realiza durante los primeros 10 segundos en los que se produce la reacción.

Método de Euler, cuyo esquema numérico se indica a continuación:

\[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + h*f(t_n,y_n) \\ \end{cases}\quad \]

% Datos
A=2;
B=0.05;
k=1.4;
h=0.1;
t0=0;
tN=10;
c=2.05;
f=@(t,y)(-k*y^2+k*c*y);
% Discretizamos
N=round((tN-t0)/h);
% Creamos el vector de tiempo t 
t=linspace(t0,tN,N+1);
% Vector con la solucion por aproximacion
y=zeros(size(t));
y(1)=B;
% Utilizamos el metodo de Euler 
for i=1:N
    y(i+1)=y(i)+h*f(t(i),y(i));
end
x=c-y;
% Dibujamos la solucion aproximada
hold on
plot(t,y,'b')
plot(t,x,'g')
hold off
legend('Cantidad de productos','Cantidad de reactivos')
xlabel('Tiempo (s)');
ylabel('Concentración (mol/L)');
% Calculamos el momento en el que coinciden las dos concentraciones:
for i=1:N
    if abs(y(i)-x(i))<0.01
        fprintf('Corte en punto (%.2f,%.2f)\n',t(i),y(i))
    end
end

Métodos del Trapecio y Runge-Kutta:

El método del Trapecio sigue siguiente método numérico: \[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + \frac{h}{2}*(f(t_n,y_n)+f(t_{n+1},y_{n+1}))\\ \end{cases}\quad \]

Mientras que el de Runge-Kutta: \[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + \frac{h}{6}*(K_1+2*K_2+2*K_3+k_4) \\ K_1=f(t_n,y_n) \\ K_2=f(t_n+\frac{1}{2}*h,y_n+\frac{1}{2}*K_1*h) \\ K_3=f(t_n+\frac{1}{2}*h,y_n+\frac{1}{2}*K_2*h) \\ K_4=f(t_n+h,y_n+K_3*h)\\ \end{cases}\quad \]

% Datos
A=2;
B=0.05;
k=1.4;
h=0.1;
t0=0;
tN=10;
c=2.05;

f=@(t,y)(-k*y^2+k*c*y);

% Discretizamos
N=round((tN-t0)/h);

% Creamos el vector de tiempo t 
t=linspace(t0,tN,N+1);

% Metodo del trapecio
yt=zeros(1,N+1);
yt(1)=B;
for i=1:N
    yt(i+1)=(-1+0.5*h*k*c+sqrt((1-0.5*h*k*c)^2+2*h*k*(yt(i)+0.5*h*f(t(i),yt(i)))))/(h*k);
end

xt=c-yt;

hold on
plot(t,yt,'k')
plot(t,xt,'b')
hold off
legend('Cantidad de productos (trapecio)','Cantidad de reactivos (trapecio)')

% Metodo de Runge-Kutta
yr=zeros(1,N+1);
yr(1)=B;
for j=1:N
    K1=f(t(j),yr(j));
    K2=f(t(j)+0.5*h,yr(j)+0.5*h*K1);
    K3=f(t(j)+0.5*h,yr(j)+0.5*h*K2);
    K4=f(t(j)+h,yr(j)+K3*h);
    yr(j+1)=yr(j)+h*(K1+2*K2+2*K3+K4)/6;
end
xr=c-yr;

% Dibujamos la solucion aproximada
figure
hold on
plot(t,yr,'r')
plot(t,xr,'g')
hold off
legend('Cantidad de productos (Runge-kutta)','Cantidad de reactivos (Runge-kutta)')

Una vez realizados los programas y obtenidas las gráficas correspondientes,se puede concluir que la velocidad de la reacción aumentará hasta llegar a su máximo valor en la intersección de ambas curvas, disminuyendo dicha velocidad hasta agotar la concentración de los reactivos.

Tras realizar el sistema se procede al cálculo correspondiente, representando gráficamente ambos métodos, donde se puede observar una variación de ambos. El punto de cruce de las curvas es diferente en respectivos métodos, siendo en Euler (1.40,1.02) y en Runge Kutta 4 (1.30,1.02), como se puede observar en la gráfica. Se argumenta que el método Euler aporta un desarrollo del sistema más preciso debido a que este método es de orden 1, dando menos error que el método de Runge Kutta 4, de orden 4. Sin embargo, esta diferencia entre ambas gráficas queda reducida a medida que se aumenta el paso h.

3 Reacciones consecutivas de Lotka

Lotka propone un modelo de reacciones consecutivas en las que se consume A, y tras determinadas reacciones se produce B. En las reacciones intermedias se crean dos nuevos compuestos, X e Y que marcarán la velocidad de reacción.

La reacción global queda definida por tres reacciones intermedias que suceden de forma consecutiva. Cabe destacar, por un lado, el hecho de que dos de ellas son autocatalíticas y, por otro lado, que todas ellas están determinadas por diferentes constantes k_i que determinarán la velocidad de cada reacción. Las tres reacciones son: [math] \left\{\begin{matrix} A+X→k_1*2*X\\ X+Y→k_2*2*Y\\ Y→k_3*B \end{matrix}\right. [/math]

Para determinar las ecuaciones diferenciales correspondientes a la velocidad de formación (o descomposición) de cada sustancia, se hace un breve estudio de cada una de ellas. En nuestro caso, se ha seguido el orden contrario, y se empieza por la última de ellas para acabar con la primera. Además, es necesario utilizar el Principio de conservación de masas para poder determinar las cuatro ecuaciones correspondientes a los cuatro compuestos presentes. El análisis seguido es el siguiente:

• Y→k₃*B. A partir de Y se formará B, regido por una constante k₃. Mediante dicha reacción se podrá definir una velocidad de transformación de B: B'=k₃*Y.
• X+Y→k₂*2*Y. En la formación de Y hay que tener en cuenta la parte proporcional a la descomposición de X (k₂*X*Y) y descontar la cantidad de Y que se descompone para la formación de B (k₃*B). Por tanto, la velocidad de transformación de Y queda determinada por: Y'=k₂*Y*X-k₃*Y.
• A+X→k₁*2*X. La formación de X tiene una parte proporcional a la descomposición de A (k₁*A*X), al que se descuenta aquella parte de X que reacciona para formar el compuesto Y (k₂*X*Y). Así pues, se obtendrá una velocidad de transformación regida por X'=k₁*A*X-k₂*X*Y.
• Derivado la ecuación correspondiente al Principio de conservación de masas (A+X+Y+B=cte) se obtiene la siguiente ecuación: A'+X'+Y'+B'=0. Sustituyendo los valores obtenidos anteriormente de X', Y' y B' se obtiene la ecuación diferencial que rige la velocidad de A: A'+x'+y'+B'=0 → A'=-x'-y'-B'=-k₁*A*x+k₂*x*y-k₂*x*y+k₃*y-k₃*y=-k₁*A*x → A'=-k₁*A*X

Por lo tanto, a partir de este sistema de ecuaciones se puede formar el siguiente PVI asociado, al que hay que añadir las condiciones iniciales, que para este caso son las concentraciones iniciales de cada compuesto: \[ PVI = \begin{cases} X'=k_1*A*X-k_2*X*Y\\Y'=k_2*X*Y-k_3*Y\\B'=k_3*Y\\A'=-k_1*A*X\\X(0)=X_0\\Y(0)=Y_0\\B(0)=B_0\\A(0)=A_0\\ \end{cases}\quad \]

Este PVI queda a continuación resuelto numéricamente mediante el método de Euler. Para ello, las condiciones iniciales elegidas son: A₀=4.7, B₀=2, X₀=2 e Y₀=1. El estudio se realiza en los primeros 200 segundos en los que se producen las reacciones. Además, está resuelto usando dos pasos, primero h=0.01 y después h=0.005.

% Datos
A0=4.7;
B0=2;
X0=2;
Y0=1;
k=1.4;
t0=0;
tN=200;
c=2.05;
k1=0.1;
k2=0.15;
k3=0.06;

y0=[A0;B0;X0;Y0];
ne=4;

% Funcion:
f=@(t,y)[-k1*y(1)*y(3);k3*y(4);k1*y(1)*y(3)-k2*y(3)*y(4);k2*y(3)*y(4)-k3*y(4)];

% Usamos h=0.01:
h=0.01;
% Discretizamos
N=round((tN-t0)/h);

% Creamos el vector de tiempo t 
t=linspace(t0,tN,N+1);

% Matriz de soluciones
y=zeros(ne,N+1);
y(:,1)=y0;

% Aplicamos el método de Euler:
for k=1:N
    y(:,k+1)=y(:,k)+h*f(t(k),y(:,k));
end

% Dibujamos la solucion aproximada
subplot(2,1,1)
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'b')
plot(t,y(3,:),'g')
plot(t,y(4,:),'k')
legend('A','B','X','Y');
hold off

% Usamos h=0.005:
h=0.005;
% Discretizamos
N=round((tN-t0)/h);

% Creamos el vector de tiempo t 
t=linspace(t0,tN,N+1);

% Matriz de soluciones
y=zeros(ne,N+1);
y(:,1)=y0;

% Aplicamos el método de Euler:
for k=1:N
    y(:,k+1)=y(:,k)+h*f(t(k),y(:,k));
end

% Dibujamos la solucion aproximada
subplot(2,1,2)
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'b')
plot(t,y(3,:),'g')
plot(t,y(4,:),'k')
legend('A','B','X','Y');
hold off

Como se puede observar en los gráficos de detalle, existe una pequeña perturbación entre ambas gráficas, las cuales han sido tomadas arbitrariamente como ejemplo. Esta perturbación es producto de las condiciones iniciales que se introducen en el programa realizado. A pesar de la existencia de dichas discordancias y a la inexistencia de un incremento del error, se podría decir, incluso concluir, que las reacciones estudiadas son estables. Esto se debe a la necesidad de aplicar un gran zoom para poder llegar a apreciar la diferencia entre las curvas de las gráficas una vez superpuestas.

A continuación se procede a hacer un breve estudio comparativo de la reacción que tiene lugar cuando, manteniendo las ecuaciones que rigen las reacciones, se cambia el valor que toma la constante k₁ y se mantienen iguales el resto de valores. Para este estudio se han tomado cuatro valores de k₁ en el intervalo (0.02,2), y se ha resuelto el PVI numéricamente mediante el método de Heun.

% Datos
A0=4.7;
B0=2;
X0=2;
Y0=1;
k=1.4;
t0=0;
tN=200;
c=2.05;
k1=0.1;
k2=0.15;
k3=0.06;

y0=[A0;B0;X0;Y0];
ne=4;

% Usamos h=0.01:
h=0.01;
% Discretizamos
N=round((tN-t0)/h);

% Creamos el vector de tiempo t 
t=linspace(t0,tN,N+1);

% Matriz de soluciones
y=zeros(ne,N+1);
y(:,1)=y0;

% Aplicamos el método de Heun
a=1;
b=linspace(0.02,2,4);
for j=1:4
    k1=b(j);
    f=@(t,y)[-k1*y(1)*y(3);k3*y(4);k1*y(1)*y(3)-k2*y(3)*y(4);k2*y(3)*y(4)-k3*y(4)];
    for i=1:N
        K1=f(t(i),y(:,i));
        K2=f(t(i)+h,y(:,i)+K1*h);
        y(:,i+1)=y(:,i)+0.5*h*(K1+K2);
    end
    subplot(2,2,a)
    hold on
    plot(t,y(1,:),'r')
    plot(t,y(2,:),'b')
    plot(t,y(3,:),'g')
    plot(t,y(4,:),'k')
    legend('A','B','X','Y');
    title(k1);
    hold off
    a=a+1;
end

La resolución de este problema se caracteriza por los distintos valores que toma k₁ en el intervalo (0.002-2). A raíz de esto se aprecian en las gráficas, también adjuntas, los cambios que se producen en las reacciones. Además, en una de las ecuaciones definidas por Lotka se puede observar la presencia de k₁: A+X=k₁*2*X. Debido a esta relación, y a partir del estudio de las gráficas, se puede concluir que a medida que el valor de la constante k₁ aumenta:

• El valor de X experimenta un crecimiento instantáneo, que es reducido inmediatamente y cuya curva va transformándose en un asíntota.
• El valor de Y experimenta una reacción muy parecida a la de X pero presenta una forma menos brusca.
• El valor de A experimenta una desaparición prácticamente instantánea al inicio la reacción .
• El valor de B experimenta un crecimiento bastante rápido.

Estas diferencias se pueden encontrar claramente entre cada una de las cuatro gráficas, de forma que se podría afirmar que a partir de k₁= 0.68 la reacción mantiene tanto las velocidades como los cambios prácticamente iguales. Para valores pequeños de k₁, como en el caso en el que su valor es 0.02 o incluso 0.1 (estudiado anteriormente), se puede apreciar una inestabilidad que comienza en torno a los 120 segundos de los dos compuestos intermedios, X e Y. Esta inestabilidad se reduce conforme el valor de k₁ aumenta, hasta hacerse inapreciable cuando la constante alcanza valores de 0.68 y superiores. Además, se observa que el tiempo en el que ocurre la totalidad de la reacción se ve reducidoa medida que el valor de k₁ aumenta. Por lo tanto, cuanto mayor sea dicho valor, menos tiempo necesita la reacción global para llevarse a cabo en su totalidad, es decir, que todo el compuesto A se transforme en compuesto B, y los compuestos intermedios X e Y desaparezcan por completo.