Reacciones complejas - Grupo 1 A

Trabajo realizado por estudiantes
Título	Reacciones complejas. Grupo 1A
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	María Ramírez Ignacio Posada Antonio López-Mateos Pablo Bueno
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción y explicación del PVI

El problema que vamos a tratar es un problema que aborda el tema de Reacciones Complejas. Para ello consideramos, según el enunciado dado, una reacción química irreversible en una solución bien mezclada. Tal y como se nos indica inicialmente en los datos del problema, suponemos también que dicha reacción dada ocurre a un volumen y temperatura constantes. La reacción química con la que trabajaremos responde a la estructura que a continuación se presenta:

                                                      A + B → C

Ésta nos indica que tenemos inicialmente 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, tal y como expresamos en el esquema de la reacción.

Otro dato también a suponer es que dicha reacción química satisface la Ley de Acción de Masas, ley que establece que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos.

Una vez conocidos todos estos datos sobre nuestro problema, nos disponemos a ir resolviéndolo respondiendo a las preguntas que se nos han planteado y que a continuación en los siguientes apartados desarrollaremos.

1.1 ¿Qué sucedería si el proceso fuese reversible?

Hasta ahora hemos planteado el PVI suponiendo el proceso irreversible. En cambio, si lo suponemos reversible, únicamente tendríamos que añadir un componente negativo en la velocidad de reacción del compuesto C, proporcional a la concentración de C, destinada a la creación de A y B.

y'(t) = k1 * (A0 - y(t)) * (B0 - y(t)) - k2*y(t)

A(t) = A0 - y(t)

B(t) = B0 - y(t)

Concentración reactivos k1=1 k2=2

Concentración reactivos k1=1 k2=3

En las gráficas vemos cómo las reactivos acaban estabilizándose.

Si k2=∞ ó k1=0 No habría reacción. En el primer caso se crearía A y B más rápido que su destrucción. En el segundo nunca se llegaría a destruir nada de A y B.

Si k2=0 ó k1=∞ La reacción no sería reversible. En el primer caso no se crearía nada de A y B por lo que sólo se crearía C. En el segundo caso la creación de C es mucho mayor que la de A y B.

2 Resolución Numérica del problema asociado

Ahora, nos disponemos a resolver numéricamente el problema de valor inicial descrito en el apartado anterior, suponiendo:

\left. \begin{array}{rcl}

    a_0 & = & 3 mol/l  \\  
    b_0 & = & 1 mol/l  \\
    k_1 & = & 1 mol/s 

\end{array} \right\}

Utilizaremos, en este caso, el Método de Euler como método numérico para la resolución del problema, eligiendo un paso h = 0.1 en los primeros 2 segundos. El código Matlab resultante tras la aplicación numérica de este método a nuestro problema de valor inicial es el siguiente:

%PVI por Euler con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    y(i+1)=y(i)+ h*(k*(3 - y(i))*(1 - y(i)));
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Una vez establecido dicho código Matlab, dibujamos en una misma gráfica tanto las concentraciones de los dos reactivos como la del producto, en función de su evolución en el tiempo, quedando las mismas de la siguiente manera:

Gráfica de las dos concentraciones y el producto. Utilizando el método de Euler y con una h = 0.1

Observamos ahora en la gráfica las curvas que describen la variación de las concentraciones de los reactivos A y B y del producto C en función del tiempo. Las tres concentraciones varían exponencialmente según el paso del tiempo y, como era de esperar, las concentraciones de A y B disminuirán a medida que vaya avanzando el mismo mientras que la concentración de C, por el contrario, irá aumentando de manera exponencial. Esto es debido a que los reactivos necesitan irse consumiendo para poder dar lugar a C, que va creciendo ya que, a medida que el tiempo transcurre, se va formando más cantidad del mismo. El punto de inicial de las tres concentraciones para un tiempo t=0 representado en la gráfica nos indica la concentración de cada elemento que tenemos en el instante inicial, antes de que comience a darse lugar la reacción, que coincide con los datos que se nos planteaban en el problema de valor inicial:

\left. \begin{array}{rcl}

    a_0 & = & 3 mol/l  \\  
    b_0 & = & 1 mol/l  \\
    c_0 & = & 0 mol/l 

\end{array} \right\}

Las gráficas de las funciones exponenciales de las concentraciones tienden a valores diferentes a medida que se avanza en el tiempo, valores que dependen de la concentración inicial de cada elemento. Observamos que la gráfica que representa la variación de la concentración de A con el tiempo (color verde) tiende a 2, la de B (color rojo) tiende a 0 y la de C (color azul) tiende a 1. Esto es así debido a tratarse de una reacción bimolecular, que, como ya se ha explicado anteriormente, significa que una molécula de A y una de B producen una de C. Por tanto, esa es la razón por la que la concentración de A pasa de a(0)= 3 mol/l a a(t)= 2 mol/l en función del tiempo t transcurrido, y así ocurre también con B y C, disminuyendo B de b(0)= 1 mol/l a b(t)= 0 mol/l y, por el contrario, aumentando C de c(0)= 0 mol/l a c(t)= 1 mol/l.

Cuando t tiende a infinito

Una vez analizada cómo es la variación de las concentraciones en función del transcurso del tiempo, resolveremos numéricamente el mismo sistema para un tiempo lo suficientemente grande como para ser capaces de establecer una aproximación de los límites de las concentraciones cuando t → ∞. Para ello, lo único a modificar en el programa anterior es el dato del tiempo, que lo aumentaremos razonablemente para lograr ver hacia dónde tienden las gráficas con el paso del tiempo, mientras que el resto de datos se mantienen.

Ahora, la gráfica del Método de Euler para un tiempo t → ∞ queda de la siguiente manera:

Gráfica de las dos concentraciones y el producto. Utilizando el método de Euler y con una h = 0.1

Como ya se podía intuir en el anterior gráfico a este, las gráficas de las funciones exponenciales de las concentraciones tienen sus límites para un tiempo t → ∞ determinados según la reacción bimolecular en la que están participando. Por tanto, tal y como ya se ha comentado, la concentración de A tiende a 2 mol/l, la de B tiende a 0 mol/l y la de C tiende a 1 mol/l, cuando la reacción tiende a un tiempo t → ∞.

2.1 Método del Trapecio

Al ser el Método de Euler un método numérico de aproximación, resolvemos ahora nuestro mismo problema de valor inicial por otros dos métodos numéricos también de aproximación, que son el Método del Trapecio y el Método de Runge Kutta. De esta manera, podemos comparar la solución al problema obtenida por los tres métodos numéricos empleados, debido a que, al ser métodos aproximados, las soluciones serán más precisas en función del método empleado. Mantenemos todos los datos iguales, utilizando también el mismo paso de tiempo.

Tendremos, por un lado, el problema de valor inicial resuelto por el Método del Trapecio, quedando de la siguiente manera:

%PVI por Trapecio con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
a = 3;
b = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    %Definimos tres variables para reducir la ecuación
    Q = 1 + (h*k*(a+b))/2;
    W = -(h*k)/2;
    E = -y(i)- h*k*(1/2)*((a - y(i))*(b - y(i)) + a*b);
    y(i+1)=(-Q + sqrt(Q^2- (4*W*E)))/(2*W);
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Acompañado de su gráfica:

Gráfica de las dos concentraciones y el producto. Utilizando el método del Trapecio y con una h = 0.1

2.2 Método de Runge-Kutta de Orden 4

Y, por otro lado, el problema de valor inicial resuelto empleando el Método de Runge Kutta, quedando así:

%PVI por Runge Kutta orden 4 con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
a = 3;
b = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    K1 = k*(a - y(i))*(b - y(i));
    K2 = k*(a-(y(i)+K1*h/2))*(b - (y(i)+K1*h/2));
    K3 = k*(a-(y(i)+K2*h/2))*(b - (y(i)+K2*h/2));
    K4 = k*(a-(y(i)+(K3*h)))*(b - (y(i)+K3*h));
    y(i+1)= y(i) + h/6*(K1+(2*K2)+(2*K3)+K4);
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Acompañado también de su gráfica:

Gráfica de las dos concentraciones y el producto. Utilizando el método de Runge Kutta de orden 4

            HASTA AQUÍ ESTÁ TODO REVISADO Y COMPROBADO Y ESTÁ BIEN

3 Segunda reacción. A partir de C obtenemos otro compuesto D

                                                    A + B → C → D

La y(t) Calculada anteriormente ahora pasaría a ser la suma de las concentraciones y(t)= yreal(t)+z(t) es decir: y(t)= [C] + [D]

  yreal(t)=y(t)- z(t)

Teniendo ahora este PVI

 yreal`(t)=y`(t)- z`(t)
 z`(t)=k*yreal(t)
 yreal`(t)=k*(a0-y(t))*(b0-y(t))-k*(y(t)- z(t))
 a0=3
 b0=1
 y(0)=0
 z(0)=0

3.1 Resolución numérica

3.1.1 Euler

%Establecemos los valores iniciales
t0=0;
tN=10;
y0=0;
k1=1;
a=3;
b=1;
k2=input('introduce k2: ');
%Discretizamos
h=input('introduce longitud de paso: ');
N=round((tN-t0)/h);
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
y1=zeros(1,N+1);
yd=zeros(1,N+1);

y1(1)=y0;
yd(1)=y0;

yy1=y0;
yyd=y0;

%Utilizamos el 'for' para crear los vectores solución
for n=1:N
  yy1=yy1+h*(k1*(a-yy1)*(b-yy1)-k2*(yy1-yyd));
  yyd=yyd+h*(k2*(yy1-yyd));
  y1(n+1)=yy1;
  yd(n+1)=yyd;
end

%Definimos las variables de las concentraciones para dibujarlas
ya=3.-y1;
yb=1.-y1;
yc=y1-yd;

hold on

plot(t,yc,'b','linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
plot(t,yd,'k','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Concentracion de D','Location','best');
hold off

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=o.1

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=o.3

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=o.1

3.1.2 Runge-Kutta de Orden 4

%C como reactivo y D como producto
t0=0;
tn=10;
h=input('introduce longitud de paso: ');
N=round((tn-t0)/h);
t=t0:h:tn;
y1=zeros(1,N+1);
yd=zeros(1,N+1);
k1=1;
k2=input('introduce k2: ');
a=3;
b=1;
for n=1:N
    
    %runge kutta 
    kc1=k1*(a-y1(n))*(b-y1(n));
    kc2=k1*(a-y1(n)-h/2*kc1)*(b-y1(n)-h/2*kc1);
    kc3=k1*(a-y1(n)-h/2*kc2)*(b-y1(n)-h/2*kc2);
    kc4=k1*(a-y1(n)-h*kc3)*(b-y1(n)-h*kc3);
    y1(n+1)=y1(n)+(h/6)*(kc1+2*kc2+2*kc3+kc4);
    kd1=k2*(y1(n)-yd(n));
    kd2=k2*(y1(n)-yd(n)-h/2*kd1);
    kd3=k2*(y1(n)-yd(n)-h/2*kd2);
    kd4=k2*(y1(n)-yd(n)-h*kd3);
    yd(n+1)=yd(n)+(h/6)*(kd1+2*kd2+2*kd3+kd4);
 
end
%concentraciones
ya=3.-y1;
yb=1.-y1;
yc=y1-yd;

hold on

plot(t,yc,'b','linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
plot(t,yd,'k','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Concentracion de D','Location','best');
hold off

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.1 k2=5

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.3 k2=5

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.1 k2=1/5