Reacciones complejas - Grupo 16 A

Trabajo realizado por estudiantes
Título	Reacciones complejas. (Grupo 16-A)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Arévalo Lecanda, Javier Buitrago Peña, Marcos Chamizo Carmona, Javier La Porta, Nicoletta
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este articulo vamos estudiar un ejercicio de modelización en el cual consideramos una reacción química irreversible en una solución bien mezclada. La reacción en cuestión es en una reacción bimolecular que consiste en que una molécula de A y una de B se juntan para crear una de C.

 A + B → C

Para nuestro ejercicio, supondremos varias cosas:

• La reacción ocurre para un volumen constante.
• La reacción ocurre para una temperatura constante.
• La reacción 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 (A y B)

2 Problema de Valor Inicial (PVI) [Apartado 1]

Lo primero que debemos hacer para resolver este problema es pasar el fenómeno químico a un problema de valor inicial. Para eso usamos la ley de acción de masas.

En esta ecuación tendremos varias constantes

• k1: Esta contante es la que nos da la proporcionalidad de la reacción.
• a0: Esta constante es la que determina la concentración inicial de A.
• b0: Esta constante es la que determina la concentración inicial de B.

Ademas de estas constantes tendremos también las variables:

• y(t): Esta variable nos da la cantidad de producto C que se formaen el tiempo, como el volumen es constante coincide con la concentración de C a lo largo de t.
• y'(t): Es la derivada de y(t) respecto al tiempo ó lo que es lo mismo la velocidad con que se produce C, es decir, la velocidad de la reacción química.

Tenemos también que mencionar la concentración de A respecto al tiempo que sera (a0−y(t)) y la de B que sera (b0−y(t)). Para nuestra ecuación tomaremos el tiempo como estrictamente mayor que 0: t>0.

Podemos entonces deducir que la ecuación sera:

 y′(t)=k1(a0−y(t))(b0−y(t))

3 Ecuación diferencial si el proceso es reversible [Apartado 2]

Hasta ahora hemos supuesto que el proceso no era reversible, esto significa que A y B se juntaban para crear C. Si ahora suponemos que el proceso es reversible, ademas de A y B juntándose para crear C también tendremos C descomponiéndose para formar A y B. La ecuación diferencial se veria afectada por este cambio.

Tendríamos que añadir nuevos parámetros:

• k2: Ademas de tener k1, necesitaremos también conocer k2 que es la proporcionalidad de la reacción reversible.
• c0: Esta constante es la que determina la concentración inicial de C.

Ahora solo tendríamos que reescribir la ecuación diferencial anterior restando le la descomposición de C en A y B.

 y′(t)=k1(a0−y(t))(b0−y(t))−k2(c1+y(t))

4 PVI con el método Euler [Apartado 3]

Vamos ahora a resolver el problema de valor inicial numéricamente con el método de Euler, de aquí en adelante tomaremos la reacción como irreversible.

Primero damos valores a nuestros parametros según el enunciado del trabajo:

• a0=3 mol/l
• b0=1 mol/l
• k1=1 mol/s
• h=0.1
• t: [0,2]

Nuestro codigo sera:

a0=3; b0=1; %Concentraciones iniciales de A y B
t0=0; tf=2; %Tiempos iniciales y finales
y0=0; %Tomamos la concentracion de C inicial como 0
k1=1; 
h=0.1; %Paso
N=(tf-t0)/h;
t=t0:h:tf; 
y=zeros(1,N+1);
yy=y0;
for n=1:N
    yy=yy+h*(k1*(a0-yy)*(b0-yy));
    y(n+1)=yy;
end
A=a0-y; %Concentracion de A
B=b0-y; %Concentracion de B
hold on
plot(t,y,'r')
plot(t,A,'b')
plot(t,B,'g')
legend('C','A','B')
hold off

La gráfica obtenida es:

izquierda

En este gráfico tenemos: la concentración de A en azul, la concentración de B en verde y la concentración de C en rojo. Se puede apreciar fácilmente lo que ocurre en la reacción, las concentraciones de A y de B disminuyen en el tiempo a la misma velocidad en la que la concentración de C crece. Podemos también observar que el limitante en esta reaccion es el B.

5 Resolución numérica cuando t→∞ [Apartado 4]

Como no podemos operar con t=∞ debemos parar el programa cuando la velocidad de reacción sea cercana a cero, como nuestra reaccíon es irreversible ésta terminará en el momento en que se acabe el reactivo limitante, pero como la velocidad es directamente proporcional a esa concentración ésta nunca llegará a ser cero sino que tiende hacia él, supondremos que la reacción ha llegado a su fin cuando llegue al entorno de 1/10000 pero podrían ser menores. El siguiente programa funciona hasta que la velocidad es cero(momento en el que nos deja de interesar seguir operando) Nuestro codigo sera:

a0=3; b0=1; %Concentraciones iniciales de A y B
t0=0;  %Tiempos iniciales y finales
y0=0; %Tomamos la concentracion de C inicial como 0
k1=1; 
h=0.1; %Paso 
yy=y0;
n=1;
y(1)=y0
t(1)=t0;
v(1)=k1*a0*b0 %velocidad inicial
while v>0.0001
    yy=yy+h*k1*(a0-yy)*(b0-yy);
    t(n+1)=t(n)+h;
    v(n)=k1*(a0-yy)*(b0-yy)
    n=n+1;
end
A=a0-y; %Concentracion de A
B=b0-y; %Concentracion de B
hold on
plot(t,y,'r')
plot(t,A,'b')
plot(t,B,'g')
legend('C','B','A')
hold off

6 Resolución por los métodos del trapecio y de Runge Kutta [Apartado 5]

6.1 Método del trapecio y de Runge Kutta

En éste apartado se aproxima numéricamente a la solución real utilizando los métodos de euler, del trapecio y de Runge kutta. Dichas soluciones quedan grabadas en tres vectores distintos y, yt e yk respectivamente.

t0=0;
tn=5;
h=0.1;
N=(tn-t0)/h;
t=t0:h:tn;
y=zeros(1,N+1); %euler
yk=zeros(1,N+1);%trapecio
yt=zeros(1,N+1);%Runge kutta
y(1)=0;
yt(1)=0;
yk(1)=0;
k=1;
a=3;
b=1;
for n=1:N
    y(n+1)=y(n)+h*(k*(a-y(n))*(b-y(n))); %euler
    yk(n+1)=(-2*h-1+sqrt((2*h+1)^2-4*(-0.5*h)*(-yk(n)-0.5*h*(6-4*yk(n)+yk(n)^2))))/(-h); %trapecio
    k1=(3-yt(n))*(1-yt(n));
    k2=(3-yt(n)-0.5*h*k1)*(1-yt(n)-0.5*h*k1);
    k3=(3-yt(n)-0.5*h*k2)*(1-yt(n)-0.5*h*k2);
    k4=(3-yt(n)-h*k3)*(1-yt(n)-h*k3);
    yt(n+1)=yt(n)+(h/6)*(k1+2*k2+2*k3+k4);%runge kutta
end
A=a-yk;
B=b-yk;
A1=a-y;
B1=b-y;
A2=a-yt;
B2=b-yt;
hold on
plot(t,y,'r')
plot(t,A,'g')
plot(t,yt,'b')
plot(t,A1,'r')
plot(t,B1,'r')
plot(t,yk,'g')
plot(t,A2,'b')
plot(t,B2,'b')
plot(t,B,'g')
legend('euler','trapecio','runge kutta')

7 Sistema de dos ecuaciones diferenciales [Apartado 6]

Ahora se presenta la siguiente reacción:

                                       A + B → C → D

Son dos reacciones consecutivas con k1 y k2 como constantes de velocidad y a,b,c como concentraciones iniciales. por tanto las velocidades de reacción son respectivamente:

                                        y′(t)=k1(a−y(t))(b−y(t))

                                        z′(t)=k2(c+y(t)-z(t))

Y como y(t), z(t) designan la cantidad que se forma de producto en el instante cero

                                        y(0)=0

                                        z(0)=0

Sabemos que es un problema de valor inicial con solución y que se puede resolver en cascada puesto que la segunda función no aparece en la primera.

%en el caso de una segunda reacción con C como reactivo y D como producto
t0=0;
tn=10;
h=0.03;
N=(tn-t0)/h;
t=t0:h:tn;
y=zeros(1,N+1); 
yk=zeros(1,N+1);
z=zeros(1,N+1);
zk=zeros(1,N+1);
y(1)=0;
z(1)=0;
k1=1;
k2=5;
a=3;
b=1;
for n=1:N
    y(n+1)=y(n)+h*(k1*(a-y(n))*(b-y(n)));
    z(n+1)=z(n)+h*k2*(c+y(n)-z(n));
    %runge kutta 
    kc1=k1*(a-yk(n))*(b-yk(n));
    kc2=k1*(a-yk(n)-0.5*h*kc1)*(b-yk(n)-0.5*h*kc1);
    kc3=k1*(a-yk(n)-0.5*h*kc2)*(b-yk(n)-0.5*h*kc2);
    kc4=k1*(a-yk(n)-h*kc3)*(b-yk(n)-h*kc3);
    yk(n+1)=yk(n)+(h/6)*(kc1+2*kc2+2*kc3+kc4);
    kz1=k2*(yk(n)-zk(n));
    kz2=k2*(yk(n)-zk(n)-0.5*h*kz1);
    kz3=k2*(yk(n)-zk(n)-0.5*h*kz2);
    kz4=k2*(yk(n)-zk(n)-h*kz3);
    zk(n+1)=zk(n)+(h/6)*(kz1+2*kz2+2*kz3+kz4);
    
end
%concentraciones
C=y-z;
Ck=yk-zk;
D=z;
Dk=zk;
hold on
plot(t,C,'r')
plot(t,Ck,'g')
plot(t,D,'r')
plot(t,Dk,'g')
legend('euler','runge kutta')