Reacciones Complejas. Grupo 25C

Trabajo realizado por estudiantes
Título	REACCIONES COMPLEJAS. GRUPO 25C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Gálvez Aparici, Antonio Megino León, Guillermo Popa, Silviu Sistac Ara, Alejandro Veiga López, Roberto
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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:

centro

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:

centro

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:

centro

4 Resolución del PVI por el método 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)) aplicaremos la siguiente formula:

centro

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 Código 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 Resolución del PVI por el método de Euler cuando t→∞

Una vez realizado el método de Euler para un intervalo de 2 segundos nos damos cuenta de que las concentraciones convergen y podemos presuponer que se mantendrán constantes, pero vamos a comprobarlo para el límite t→∞. Para ello basta con realizar el mismo programa que el anterior pero con la salvedad de modificar el tiempo hasta un valor que se pueda considerar comparable a ∞ para comprobar que se mantienen constantes las concentraciones, tomaremos 10 segundos ya que la reacción tarda poco tiempo:

%Euler
clear all
%Datos del problema
a0=3;
b0=1;
t0=0;
tN=10;
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

centro

Como habíamos previsto las concentraciones tienden a A=2 mol/l; B=0 mol/l; C=1 mol/l. Este es debido a que el compuesto B es el limitante de la reacción, por lo que al agotarse, se mantienen estables las demás concentraciones.

6 Resolución del PVI 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

Una vez realizados los tres métodos numéricos de resolución de PVI podemos observar la gran similitud entre ellos, tendiendo todos a la misma concentración a los dos segundos manteniéndose esta estable. El método del Trapecio y el de Runge-Kutta son más similares entre ellos porque tienen mayor orden que el de Euler (segundo orden y cuarto orden respectivamente), a su vez siendo más precisos como se demuestra en la gráfica teniendo ambos una curvatura mas similar (más suave que la de Euler).

7 Reacción consecutiva

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

centro

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.

centro

• k₂: Constante de la velocidad de reacción D, en función de la temperatura, que al ser constante también lo sera k2
• y₁'(t): Es la derivada de la concentración en el tiempo, es decir, la velocidad de la reacción química de C
• y₂'(t): Velocidad de la reacción química de D

7.1 Resolucion por Euler

%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

Con paso h=0.1 y k₂=5

centro

Con paso h=0.3 y k₂=5

centro

Con paso h=0.1 y k₂=1/5

centro

Partiendo de los compuestos A y B se van a formar los productos C y D (dependiendo este último de C). Poco antes del segundo la concentración de C alcanza su máximo y a partir de aquí comienza a descender (cada vez a menor velocidad) hasta desaparecer convirtiéndose en todo en D. Como ocurría anteriormente el compuesto B desaparece a los dos segundos y el compuesto A se estabiliza en la concentración de 2 mol/l. Finalmente el compuesto D se mantiene estable con una concentración de 1 mol/l.

Podemos observar que con h=0.1 y K₂=5 surge todo bastante deprisa, ajustándose B y C a los dos segundos (debido a que tenemos una alta velocidad de reacción hacia D).

En la segunda gráfica al aumentar el paso a h=0.3 se distorsiona la aproximación dándonos valores poco precisos que no se relacionan con la realidad a tiempos bajos, pero si cuando sobrepasamos los tiempos de estabilidad (al llegar a los límites de la reacción).

Finalmente volviendo al paso de h=0.1 (precisión suficiente) y con una velocidad de reacción k₂=1/5 (mucho menor que las anteriores) obtenemos una gráfica mucho más tendida y suave en el tiempo debido a que tarda más en transformarse los compuestos en el producto D. En esta ocasión los compuestos A y B se estabilizan a los dos segundos porque k₁ no varía pero en cambio el compuesto C tarda bastante en transformarse en D (quedando todavía 0.2 mol/l de C a los 10 segundos) pero con una clara tendencia a completarse la segunda reacción.

7.2 Resolucion por Runge Kutta

%Runge Kutta
clear all
%Datos del problema
a0=3;
b0=1;
t0=0;
tN=10;
y0=0;
h=0.1;
k1=1;
k2=5;
%Calculamos el número de subintervalos
N=round((tN-t0)/h);
t=t0:h:tN;
y=zeros(1,N+1);
z=zeros(1,N+1);
y(1)=y0;
z(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);
    z1=k2*(y(i)-z(i));
    z2=k2*(y(i)-(z(i)+(h/2)*z1));
    z3=k2*(y(i)-(z(i)+(h/2)*z2));
    z4=k2*(y(i)-(z(i)+h*z3));
    z(i+1)=z(i)+(h/6)*(z1+2*z2+2*z3+z4);
end
%Soluciones
yc=y-z;
plot(t,y,'c',t,yc,'b','linewidth',1.5);

a=3.-y;
b=1.-y;
hold on
plot(t,a,'k','linewidth',1.5);
hold on
plot(t,b,'r','linewidth',1.5);
legend('D','C','A','B');
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off

Con h=0.1 y k₂=5

centro

Con h=0.3 y k₂=5

centro

Con h=0.1 y k₂=1/5

centro

Por el método de Runge-Kutta obtenemos unas gráficas bastante similares a las obtenidas por ek método de Euler, salvando la exactitud y la precisión entre ambos métodos. Donde realmente se nota la diferencia es en el segundo caso que como ya dijimos anteriormente el paso de h=0.3 distorsiona la solución, aunque en menor medida que con Euler. Los resultados son los mismos estabilizándose A en 2mol/l y D en 1 mol/l (B y C desaparecen con el tiempo, a los 2 segundos si k₂=5 y a los 20 si k₂=1/5).