Diferencia entre revisiones de «Reacciones Complejas. Grupo 25C»

Revisión del 19:38 6 mar 2015

1 Introducción

Se considera una reacción química irreversible en una solución bien mezclada. Supondremos que la reacción ocurre para un volumen y temperatura constantes. Al inicio se encuentran dos reactivos A y B que van formando un producto C en lo que se conoce como una reacción bimolecular, es decir, una molécula de A y una de B producen una de C:

A + B → C

Supondremos también que se satisface la ley de acción de masas que establece que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos.

2 Interpretación del enunciado y PVI. Concentración de C

Vamos a comprobar que la concentración del producto C a lo largo del tiempo puede obtenerse a partir de esta ecuación, la de la ley de masas:

y'(t) = k₁(a₀ − y(t))(b₀ − y(t)), siendo t>0

A continuación definiremos las constantes y variables que vamos a utilizar:

• k₁: Constante de la velocidad de reacción, en función de la temperatura, que al ser constante también lo sera k1
• a₀: Concentración inicial de A
• b₀: Concentración inicial de B
• y(t): Define la evolución de la concentración de C a lo largo del tiempo
• y'(t): Es la derivada de la concentración en el tiempo, es decir, la velocidad de la reacción química

Por lo tanto, según nuestra ecuación diferencial, las concentraciones de A y B varían en el tiempo disminuyendo según, (a₀-y(t)) y (b₀-y(t)) respectivamente. Tomaremos el tiempo estrictamente mayor que 0 (t>0), de hecho para t=0 la concentración de C será nula. Estamos ante un PVI.

Problema de valor inicial: Para el intervalo t>0 la ecuación diferencial es continua y derivable en dicho dominio, por tanto tiene solución que sera única porque su primera derivada es continua en este dominio. Cumple el teorema de Cauchy (o PVI).

3 Proceso Reversible

En caso de que el proceso sea reversible nuestra ecuación necesitara un termino adicional que defina como se producen los reactivos. Para esto necesitaremos una nueva constante de velocidad de reacción (k₂) dado que suponemos que se sigue satisfaciendo la ley de masas. Este nuevo termino disminuirá la cantidad de producto, por lo tanto el termino aparecerá restando en nuestra ecuación.

La ecuación resultante es la siguiente:

y'(t) = k₁(a₀ − y(t))(b₀ − y(t))-k₂*y(t), siendo t>0

4 Resolucion del PVI por el metodo de Euler

En este apartado resolveremos el problema de valor inicial mediante el método de Euler (primer orden) que se trata de un método explícito. Sabiendo que f(t,y)=k₁(a₀-y(t))(b₀-y(t)) y aplicando: y_n+1 = y_n + h*f(t_n, y_n). Despejando obtenemos: y_n+1 = y_n + h*k₁(a₀-y(t))(b₀-y(t))

Suponiendo estos parámetros:

• a₀=3[mol/L]
• b₀=1[mol/L]
• k₁=1[mol/s]
• h=0.1(salto)
• t=[0,2] (segundos)

4.1 Codigo en Matlab

%Euler
clear all
%Datos del problema
a0=3;
b0=1;
t0=0;
tN=2;
y0=0;
h=0.1;
k1=1;
%Calculamos el número de subintervalos
N=round((tN-t0)/h);
%Definimos la variable independiente
t=t0:h:tN;
%Definimos el vector y
y=zeros(1,N+1);
y(1)=y0;
for i=1:N
    y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));
end
%Dibujamos la grafica
hold on
plot(t,y,'linewidth',2);
plot(t,(a0-y),'g','linewidth',2);
plot(t,(b0-y),'r','linewidth',2);
legend('C','A','B');
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off

4.2 Gráfica método de Euler

centro

Como se observa en la gráfica, la concentración de C aumenta en la proporción que disminuyen los reactivos A y B. Podemos ver que la velocidad disminuye en el tiempo según una ley asintótica hasta ser nula cuando se agota el reactivo limitante, en nuestro caso el B.

5 Resolucion cuando t→∞

6 Resolución numérica con los métodos del Trapecio y Runge-Kutta

Una vez resuelto el método de Euler, procedemos a resolver el problema con los métodos del Trapecio y Runge-Kutta (cuarto orden) con las siguientes condiciones iniciales:

• a₀=3[mol/L]
• b₀=1[mol/L]
• k₁=1[mol/s]
• h=0.1(salto)
• t=[0,2] (segundos)

6.1 M.Trapecio

Para resolver el problema por el método del Trapecio debemos tener en cuenta que no es un método explicito como lo era el método de Euler, sino que es implícito. Sabiendo que f(t,y)=k₁(a₀-y(t))(b₀-y(t)) aplicaremos la formula del trapecio:

centro

para despejar la componente y_n+1 que es la que utilizaremos para operar el método numéricamente.

El código de Matlab para el método del Trapecio sera el siguiente:

%trapecio
clear all
%datos del problema
t0=0;
tN=2;
h=0.1;
k1=1;
y0=0;
a0=3;
b0=1;
%creamos el vector y, que representara cada una de las soluciones
N=round((tN-t0)/h);
y=zeros(1,N+1);
% defino la variable independiente
t=t0:h:tN;
%defino el primer valor
y(1)=y0;
%resolucion de la ecuacion
for i=1:N
    R=(h*k1*a0*b0)/2+y(i)+((h*k1)/2)*(a0-y(i))*(b0-y(i));
    Q=((h*k1*(a0+b0))/2)+1;
    y(i+1)=(Q-sqrt(Q^2-2*h*k1*R))/(h*k1);
end
%grafico
hold on
plot(t,y,'c','linewidth',2);
plot(t,(a0-y),'k','linewidth',2);
plot(t,(b0-y),'r','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off

centro

6.2 M.Runge-Kutta

El método de Runge-Kutta (cuarto orden) se trata de un método explícito al igual que el de Euler, por lo que no será necesario despejar ninguna componente para operar el método numéricamente. Procedemos directamente al igual que hicimos con el método de Euler:

%Runge Kutta
clear all
%Datos del problema
a0=3;
b0=1;
t0=0;
tN=2;
y0=0;
h=0.1;
k1=1;
%Calculamos el número de subintervalos
N=round((tN-t0)/h);
t=t0:h:tN;
y=zeros(1,N+1);
y(1)=y0;
for i=1:N
    y1=k1*(a0-y(i))*(b0-y(i));
    y2=k1*(a0-(y(i)+(h/2)*y1))*(b0-(y(i)+(h/2)*y1));
    y3=k1*(a0-(y(i)+(h/2)*y2))*(b0-(y(i)+(h/2)*y2));
    y4=k1*(a0-(y(i)+h*y3))*(b0-(y(i)+h*y3));
    y(i+1)=y(i)+(h/6)*(y1+2*y2+2*y3+y4);
end
%Soluciones
plot(t,y,'c','linewidth',2);
a=3.-y;
b=1.-y;
hold on
plot(t,a,'k','linewidth',2);
plot(t,b,'r','linewidth',2);
legend('C','A','B');
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off

centro

7 Reacción consecutiva

Supongamos ahora que se produce una reacción consecutiva de la forma:

A + B →_k1 C →_k2 D

El PVI es un sistema de ecuaciones en el cual, la primera refleja la producción de C y la segunda representa la producción de D, estando esta relacionada con la cantidad de producto C que se ha producido en función del tiempo. Los valores en t=0 son las concentraciones iniciales de los productos, es decir, nulos.

• y1'(t)=k₁*(a₀-y1(t))*(b₀-y1(t))
• y2'(t)=k₂*(y1(t)-y2(t))
• y1(0)=y2(0)=0

7.1 Resolucion por Euler

Con paso h=0.1

%Euler
clear all
%Datos del problema
a0=3;
b0=1;
t0=0;
tN=10;
y10=0;
y20=0;
h=0.1;
k1=1;
k2=5;
%Calculamos el número de subintervalos
N=round((tN-t0)/h);
%Definimos la variable independiente
t=t0:h:tN;
%Definimos el vector y
y1=zeros(1,N+1);
y2=zeros(1,N+1);
y1(1)=y10;
y2(1)=y20;

for i=1:N
    y1(i+1)=y1(i)+h*(k1*(a0-y1(i))*(b0-y1(i)));
    y2(i+1)=y2(i)+h*(k2*(y1(i)-y2(i)));
end
%Dibujamos la grafica
hold on
C=y1-y2;
D=y2;
plot(t,C,'c',t,D,'b','linewidth',1.5);
plot(t,(a0-y1),'k','linewidth',1.5);
plot(t,(b0-y1),'r','linewidth',1.5);
legend('C','D','A','B');
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off