Diferencia entre revisiones de «Reaccion Autocatalisis Grupo 2B»

Revisión del 17:22 26 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 si 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, es decir, 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, nos ayudaremos del principio de conservación de la materia, por el cual la suma de la masa de los reactivos y los productos se mantendrá constante a lo largo del tiempo. Así pues, darán como producto un tercer compuesto:

                                                      A+B → k * 2B

La resolución de la ecuación diferencial para y se apoya en las dos leyes mencionadas con anterioridad. De esta forma, determinamos las siguientes ecuaciones:

                                                        x=A+B          
                                                        y=2B

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

A partir del principio de conservación de la masa obtenemos que

                                                      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=kAB, por lo tanto, v=y’=-x’. Por lo que concluimos que y’=kAB en la que sustituimos 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) concluimos que: x=cte-y.

Y sustiyendo en la fórmula anteriormente demostrada, obtenemos que:

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

A partir de la cual se forma el PVI correspondiente: [math] \left\{\begin{matrix} y'=-k*y^2+k*y*cte\\ y(0)=y0 \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

Según los datos teóricos, A=2 y B=0.05 mol/L, 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

2 Resolución del PVI

La resolución del PVI se ha llevado acabo por diferentes métodos. Para ellos 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:

\[ 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étodo 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:

% 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),y(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 reactivos(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 aumentara 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,0.88) . 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.

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. Las ecuaciones son las siguientes:

                                               x'=k1*A*x-k2*x*y
                                                y'=k2*x*y-k3*y
                                                   B'=k3*y
                          A'+x'+y'+B'=0 → A'=-x'-y'-B'=-k1*A*x+k2*x*y-K2*x*y+k3*y-k3*y=-k1*A*x

A partir de este modelo se puede crear un PVI asociado al sistema de ecuaciones anteriormente descrito, que a continuación queda resuelto numéricamente mediante el método de Euler. Para ello, las condiciones iniciales elegidas son las cantidades iniciales de los compuestos, que son: A0=4.7, B0=2, X0=2 e Y0=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

Apartado 6

Apartado 7

% 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