Diferencia entre revisiones de «Reacciones con autocatálisis 4-C»

Revisión del 00:08 13 mar 2015

1 Introdución: Reacción bimolecular: A+B[math]\longrightarrow[/math] C

Tenemos una reacción autocatalítica bimolecular irreversible que ocurre a una presión y temperatura concretas. Esto implica que la constante k que rige la reacción permanecerá constante. Autocatalítica implica que A se va transformando en B para dar lugar a B. Bimolecular significa que una molecula de A con una de B producirán dos de B. Irreversible que solo ocurre en una dirección.

                                                                             [math]A+B \longrightarrow C [/math].

Esta reacción obviamente cumplirá el principio de conservación de la masa, es decir, la suma de las masas de los reactivos deberá ser igual a la de los productos, o dicho de otra forma, la velocidad de destrucción es igual a la de creación. Además sabemos que la velocidad de reacción es igual al producto de los reactivos multiplicado por la constante cinética.

[math] X_{0}=1 mol/l[/math] e [math]Y_{0} =0.01 mol/l[/math] ,siendo [math]Y=B[/math] , [math]X=A [/math] y [math] K=1 [/math] La ley de acción de masas nos dará la siguiente ecuación:

[math] y'(t)=kx(t)y(t)[/math] 

De la ley de la conservación de la masa obtenemos la siguiente ecuación:

[math] y'(t)+x'(t)=0[/math]

Integrando esta última ecuación, sustituyendo el valor de x(t) en la ecuación de la acción de masas y con la concentración inicial del producto B obtenemos el siguiente PVI:

[math] \begin{cases} y'(t)=f(t,y)=y(t)(1.01-y(t))\\y_{0}(t=0) =0.01  mol/l\end{cases} [/math]                                                                                                                                  

Para estudiar la existencia e unicidad de nuestro PVI habrá que estudiar la intersección de una bola centrada en el punto B((0,0.01),r) con el dominio, que al ser un polinomio es todo R², y existirá un r>0 tal que f es continua en dicha intersección, por lo que admite por lo menos una solución. Su derivada también es continua en el intervalo nombrado por lo que tendrá una única solución.

2 Resolución del P.V.I mediante los métodos numéricos: Euler, trapecio y Rounge-Kutta

En este apartado resolveremos el PVI mediante diferentes métodos numéricos y compararemos en una gráfica los errores cometidos entre ellos y con la solución exacta que hemos obtenido mediante bernoulli y con un cambio de variable.

figure 1

clear all
clf
clc
%Trabajo 2: Reacciones con autocatálisis

%Datos
h=0.1;
t0=0;
y0=0.01;
x0=1;
tN=10;

%Calculamos los vectores(Apartados 2 y 3)
t=t0:h:tN;          % vector de tiempos
N=(tN-t0)/h;
y=zeros(1,length(t)); %vector de concentración de B según euler
x=zeros(1,length(t)); %vector de concentración de A según euler
z=zeros(1,length(t)); %vector de concentración de B según trapecio
r=zeros(1,length(t)); %vector de concentración de B según rounge-kutta
e=zeros(1,length(t)); %vector de concentración de B exacto
y(1)=y0;
x(1)=1;
z(1)=y0;
r(1)=y0;

for i= 1:N;
    %Euler
    y(i+1)=y(i)+h*(y(i)*(1.01-y(i)));
    x(i)=1.01-y(i);
    %Trapecio
    z(i+1)=((1.01*(h/2)-1)+sqrt((1-1.01*(h/2))^2-2*h*((h/2)*(z(i)^2)-(1+1.01*(h/2))*z(i))))/h;
    %Rounge-kutta
    K1=r(i)*(1.01-r(i));
    K2=(r(i)+1/2*K1*h)*(1.01-(r(i)+1/2*K1*h));
    K3=(r(i)+1/2*K2*h)*(1.01-(r(i)+1/2*K2*h));
    K4=(r(i)+(K3*h))*(1.01-r(i)+(K3*h));
    r(i+1)=r(i)+(h/6)*(K1+2*K2+2*K3+K4);
    %Exacta
    e(i+1)=(1.01/(100*exp(-1.01*t(i+1))+1)); %obtenida mediante bernoulli mediante el cambio u=1/y
end

figure(1)
hold on
plot(t,y);
plot(t,x,'r');
xlabel('tiempo(s)'); ylabel('concentración(mol/l)');
legend('Elemento B','Elemento A','Location','best');
title('Gráfico de las concentraciones de A y B');
hold off


figure(2)
hold on
plot(t,y,'b') %euler
plot(t,z,'g') %trapecio
plot(t,r,'k') %rounge-kutta
plot(t,e,'r') %exacta
xlabel('tiempo(s)'); ylabel('concentración de B(mol/l)');
legend('Euler','Trapecio','Rounge-kutta','Exacta','Location','best');
title('Gráfico de la concentración de B según los diferentes métodos')
hold off

figure 2

Zoom figure 2

figure(3)
%gráfico de errores
ey=zeros(1,length(t));
et=zeros(1,length(t));
er=zeros(1,length(t));
ey(1)=0; et(1)=0; er(1)=0;
for i=1:N
ey(i+1)=abs(y(i+1)-e(i+1)); %errores euler
et(i+1)=abs(z(i+1)-e(i+1)); %errores trapecio
er(i+1)=abs(r(i+1)-e(i+1)); %errores rounge-kutta
end

hold on
plot(t,ey,'b')
plot(t,et,'g')
plot(t,er,'k')
legend('Errores euler','Errores trapecio','Errores rounge-kutta','Location','best');
xlabel('tiempo(s)'); ylabel('error cometido(mol/s)');
title('Grafico de errores')
hold off

figure 3

Interpretación: Observamos que el método más exacto es el del trapecio porque es un método ímplicito.Esto lo podemos observar en el gráfico figure 3.

También hay que resaltar que nuestro valor de k es 1, que es un valor alto, como debe de ser al ser una reacción autocatalítica, ya que estas reacciones tienen la característica de que ocurren muy rápidamente. En nuestro caso en 10 segundos se produce por completo la reacción. Sin embargo no es lo único que influirá en la velocidad de la reacción. También hay que tener en cuenta la concentración de los productos en cada instante. Esto viene muy bien representado en la gráfica figure 2, al comienzo hay muy poca concentración de B, lo que impedirá una gran transformación al principio, sin embargo, en cuanto empieza a haber más cantidad de B, al ser autocatalítica, una mayor concentración provoca una gran rapidez de reacción, una mayor pendiente de las rectas, y en muy poco tiempo desaparece gran cantidad de A.

Cuando los valores de A ya son pequeños, la baja concentración de este producto produce una disminución en la velocidad de reacción, que sería la última parte del gráfico. En nuestra gráfica podemos observar que la cantidad que se genera de B y la que se consume de A son iguales, ya que los gráficos son simétricos, por lo que la suma de la cantidad de los dos componentes en cualquier momento daría 1.01 que es también la cantidad inicial suma de las dos. En el momento en el que tenemos la misma cantidad de A y de B será justamente cuando los dos tengan (1+0.01)/2, como podemos observar en el gráfico (punto de corte de ambas graficas). Las conclusiones que podemos sacar de esto son dos 1.Se conserva la masa en todo instante, lo que se crea de uno se consume de otro. 2.Como lo que se crea de uno es lo que se consume de otro significa que no hay ningún otro factor exterior que afecte a nuestra reacción y que pudiera aumentar o disminuir la cantidad de nuestros componentes(conservación de la masa).

3 Resolución del P.V.I mediante un sistema de ecuaciones por los métodos numéricos: Euler y Rounge-Kutta

En este apartado, queremos obtener los mismos resultados pero mediante un sistema de ecuaciones, que resolveremos tanto por Euler como por Rounge-kutta.

%Por sistemas de ecuaciones (Apartado 4)

xx=zeros(1,length(t));  %vector de concentración de A según euler
yy=zeros(1,length(t));  %vector de concentración de B según euler
ra=zeros(1,length(t));  %vector de concentración de A según rounge-kutta
rb=zeros(1,length(t));  %vector de concentración de B según rounge-kutta
xx(1)=x0;
yy(1)=y0;
ra(1)=x0;
rb(1)=y0;

for i=1:N
    %Sistema con Euler
    xx(i+1)=xx(i)-h.*(xx(i).*yy(i));
    yy(i+1)=yy(i)+h.*(xx(i).*yy(i));
    %Sistema con Rounge-kutta de orden 4
    K1a=-1*ra(i)*rb(i);
    K1b=1*ra(i)*rb(i);
    K2a=-1*(ra(i)+1/2*K1a*h)*(rb(i));
    K2b=1*(rb(i)+1/2*K1b*h)*(ra(i));
    K3a=-1*(ra(i)+1/2*K2a*h)*(rb(i));
    K3b=1*(rb(i)+1/2*K2b*h)*(ra(i));
    K4a=-1*(ra(i)+K3a*h)*(rb(i));
    K4b=1*(rb(i)+K3b*h)*(ra(i));
    ra(i+1)=ra(i)+(h/6)*(K1a+2*K2a+2*K3a+K4a);
    rb(i+1)=rb(i)+(h/6)*(K1b+2*K2b+2*K3b+K4b);
end

figure (4)
hold on
plot(t,xx,'b')
plot(t,yy,'r')
plot(t,ra,'g')
plot(t,rb,'k')
xlabel('tiempo(s)'); ylabel('concentración(mol/l)');
legend('Concentracion de A con Euler','Concentracion de B con Euler','Concentración de A con Rounge-kutta','Concentración de B con Rounge-kutta','Location','best');
hold off

figure4

4 Reacción consecutiva de Lokta

Ahora tenemos una nueva situación con tres reacciones. Las dos primeras son autocatalíticas y por lo tanto actúan como la estudiada anteriormente. La reacción total consume A para producir B, mientras que X e Y intervendrán en la velocidad y la cantidad de cada producto en los procesos intermedios.

[math] A + X \rightarrow _{k1} 2X [/math]

[math] X+ Y \rightarrow _{k2} 2Y [/math]

[math] Y \rightarrow _{k3} B [/math]

Las siguientes ecuaciones provienen también de la ley de conservación de masas y de la ley de acción de masas, pero en este caso influyen unas sobre otras. De la ley de conservación de la masa sabemos que la masa conjunta de todos nuestros compuestos será constante, por lo que derivando obtendremos:

[math] A'(t)+y'(t)+x'(t)+B'(t)=0[/math]

De la ley de acción de masas obtenemos las otras tres. Pensando de atrás alante la cantidad que se forme de B sera proporcional a la de Y, con la constante k3. Sin embargo, la cantidad de Y vendrá dada por la suma de lo que se forme, es decir, proporcional a las concentraciones de X*Y pero habrá que restarle lo que se consume para formar B. De igual manera la cantidad de X nos viene dada por una parte por lo que se crea en la primera ecuación menos lo que se consume en la segunda, todo ello con sus correspondientes proporciones dadas por los valores de las k_{i}.

[math] x'(t)=k_{1}x(t)A(t)-k_{2}x(t)y(t)[/math]

[math] y'(t)=k_{2}x(t)y(t)-k_{3}y(t)[/math]

[math] B'(t)=k_{3}y(t)[/math]

Para nuestro PVI utilizaremos un intervalo de tiempo de entre 0 y 200 segundos y con [math]k_{1}=k_{2}=2k_{3}=0.1[/math] y las concentraciones de A,x,y,B serán: [math]A_{0}=5, x_{0}=5*10^{-4}, y_{0}=10^{-5} B_{0}=0[/math]

Utilizaremos el método de Euler con h= 0.01 y h=0.001; y posteriormente el método de Heun con h=0.01.

GRAFICOS Y PROGRAMAS

A0=5;  %Primero introducimos valores iniciales y constantes
X0=5*10^(-4);
Y0=10^(-5);
B0=0;
k1=0.1;
k2=0.1;
k3=0.05;
t0=0;    %Intervalo de de tiempo y matrices
tN=200;
h=0.01
N=(tN-t0)/h;
t=t0:h:tN;
A=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1);
A(1)=A0;
X(1)=X0;
Y(1)=Y0;
B(1)=B0;  %Bucle para rellenar las matrices
for i=1:N
X(i+1)=X(i)+h.*(k1.*A(i).*X(i)-k2.*X(i).*Y(i));
Y(i+1)=Y(i)+h.*(k2.*X(i).*Y(i)-k3.*Y(i));
B(i+1)=B(i)+h.*(k3*Y(i));
A(i+1)=A(i)+h.*(-k1.*X(i).*A(i));
end

figure(5)
hold on

plot(t,X)   %Graficas
plot(t,Y,'r')
plot(t,B,'y')
plot(t,A,'g')
xlabel('tiempo(s)'); ylabel('concentración(mol/l)');
legend('Elemento X','Elemento Y','ElementoB','Elemento A','Location','best');
title('Gráfico de las concentraciones');
 
hold off

A continuación hemos hecho el mismo programa pero usando un tamaño de paso h=0.001, el código de matlab usado es el mismo pero con esa pequeña diferencia por lo que con el fin de no eternizar el trabajo, solo mostraremos la gráfica.

figure 5: Euler con paso h=0.01

figure 6: Euler con paso h=0.001

%Datos
A0=5;
X0=5*10^(-4);
Y0=10^(-5);
B0=0;
k1=0.1;
k2=0.1;
k3=0.05;
h=0.01;
t0=0;
tN=200;
%Calculamos los vectores
N=(tN-t0)/h;
t=t0:h:tN; %vector de tiempos
A=zeros(1,N+1);
B=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
%valores iniciales
A(1)=A0;
B(1)=B0;
X(1)=X0;
Y(1)=Y0;
%Bucle para rellenar los vectores
for i=1:N
    %Método de Heun
    K1A=(-k1.*X(i).*A(i));
    K1B=(k3*Y(i));
    K1X=(k1.*A(i).*X(i)-k2.*X(i).*Y(i));
    K1Y=(k2.*X(i).*Y(i)-k3.*Y(i));
    K2A=(-k1.*(X(i)+K1X.*h).*(A(i)+K1A.*h));
    K2B=(k3*(Y(i)+K1Y.*h));
    K2X=(k1.*(A(i)+K1A.*h).*(X(i)+K1X.*h)-k2.*(X(i)+K1X.*h).*(Y(i)+K1Y.*h));
    K2Y=(k2.*(X(i)+K1X.*h).*(Y(i)+K1Y.*h)-k3.*(Y(i)+K1Y.*h));
    
    X(i+1)=X(i)+(h/2).*(K1X+K2X);
    Y(i+1)=Y(i)+(h/2).*(K1Y+K2Y);
    B(i+1)=B(i)+(h/2).*(K1B+K2B);
    A(i+1)=A(i)+(h/2).*(K1A+K2A);
end

figure(7)
hold on
plot(t,X)
plot(t,Y,'r')
plot(t,B,'y')
plot(t,A,'g')
xlabel('tiempo(s)'); ylabel('concentración(mol/l)');
legend('Elemento X','Elemento Y','ElementoB','Elemento A','Location','best');
title('Gráfico de las concentraciones según Heun (h=0.01)');
hold off

figure 7

Interpretación: En primer lugar habrá que decir que la reacción desde el punto de vista químico sí que es estable, ya que como podemos observar terminaremos con que todo se transforma en B. Hay otros tipo de reacciones químicas, como las llamadas reacciones reloj que no se estabilizan en un valor, sino que pueden volver a estados anteriores varias veces hasta llegar finalmente a estabilizarse. En nuestro caso observamos que en unos 200 segundos podemos decir que es completamente estable. También es importante fijarnos en el tamaño de paso, esto sería decir si es estable desde el punto de vista de las ecuaciones diferenciales.El tamaño de paso puede influir mucho haciendo que una gráfica no represente en absoluto a la ecuación exacta, produciendo que esta no sea estable, sin embargo, en nuestro caso particular el cambio entre un tamañano de paso y otro es inapreciable y los dos son considerados estables.

Por otro lado lo que observamos en el gráfico es completamente concordante con el tipo de reacciones que tenemos. Partimos de una alta concentración de A, un poco de X e Y; y nada de B. En t=0 empieza la primera reacción autocatalítica: al principio a una pequeña velocidad por falta de concentración de x y poteriormente a una gran velocidad para dar lugar en su totalidad a X, que a su vez empezará otra reacción autocatalítica con Y, pero esto no se producirá con rapidez hasta que la concentración de los dos sea suficientemente alta. En último lugar se producirá la ultima reacción, la que transforma Y en B. Sin embargo algo distinto ocurre en esta reacción, ya que podemos observar que la concentración de Y no llega al máximo posible al igual que X, sino que su máximo se queda en un valor inferior. Esto se debe a que esta reacción no es autocatalítica,sino que simplemente Y se convierte en B, por lo que no es necesario ninguna concentración inicial de B que empiece la reacción ni una concentración muy alta de Y para empezar a producirse como ocurría en las reacciones anteriores. Esto se ve reflejado en el gráfico, dando lugar finalmente a que todo reacciona y solo obtendremos componente B. Sin embargo, la velocidad de reacción para producirse B es menor, esto se debe a que no es una reacción autocatalítica.