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 | |
Contenido
1 Introducción
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.
2 Explicar PVI y definir por qué tiene solución única
3 ¿Qué sucedería si el proceso fuese reversible?
4 Resolución Numérica por el Método de Euler
Ahora, nos disponemos a resolver numéricamente el problema de valor inicial descrito en el apartado anterior, suponiendo a0 = 3 mol/l, b0 = 1 mol/l y k1 = 1 mol/s. 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:
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: a0= 3 mol/l , b0= 1 mol/l , c0= 0 mol/l.
Las funciones exponenciales tienden a valores diferentes para un tiempo t infinito según 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 para un tiempo t que tiende a infinito, 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.
5 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:
6 Resolucion por trapecio y runge kutta
6.1 Trapecio
%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
6.2 Runge Kutta
%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
7 Apartado 6
%Establecemos los valores iniciales
t0=0;
tN=10;
y0=0;
k1=1;
a=3;
b=1;
k2=1/5;
%Discretizamos
h=0.1;
N=(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-yy2));
yyd=yyd+h*(k2*(y1(n)-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
8 Apartado 7
%C como reactivo y D como producto
t0=0;
tn=10;
h=0.01;
N=(tn-t0)/h;
t=t0:h:tn;
y1=zeros(1,N+1);
yd=zeros(1,N+1);
k1=1;
k2=5;
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)=yk(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
%C como reactivo y D como producto
t0=0;
tn=10;
h=0.01;
N=(tn-t0)/h;
t=t0:h:tn;
y1=zeros(1,N+1);
yd=zeros(1,N+1);
k1=1;
k2=1/5;
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)=yk(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
