Reacciones Complejas Grupo11C

Trabajo realizado por estudiantes
Título	Reacciones complejas (Grupo 11C)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Martín Salmerón, Carlos Molina de Francia, Víctor Pérez Abellán, Ignacio David Rueda, Juan Manuel Vallino, Carlos
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Consideramos una reacción química irreversible en una solución bien mezclada, suponiendo que la reacción ocurre a volumen y temperatura constantes y que cumple la ley de acción de masas. Ésta última nos indica la relación de proporcionalidad entre la velocidad de reacción y las concentraciones de los reactivos.

                                                              A + B → C

Como se puede observar, en nuestra reacción inicialmente intervienen dos reactivos A y B que forman un producto C mediante una reacción bimolecular, es decir, para formar una molécula de C hará falta una de A y una de B.

2 Ecuación diferencial: validación e interpretación

Supuesta la reacción bimolecular irreversible anteriormente mencionada y dado que se cumple la acción de masas, podemos establecer la ecuación diferencial:

                                            y′(t)=k1(a0−y(t))(b0−y(t)) para t>0

Cuyos elementos podemos interpretar como:

.y′(t) velocidad de reacción a lo largo del tiempo.

.y(t) concentración del producto C a lo largo del tiempo.

.a0 es la concentración inicial del reactivo A.

.b0 será la concentración inicial de B.

.k1 es la constante específica de la velocidad de reacción, cuyas unidades se adaptan en base al número de reactivos [s^(-1)*mol^(1-n)*L^(n-1)], siendo n el número de reactivos.En nuestro caso, la condición de que la temperatura sea constante, implica que la constante específica de la velocidad de reacción es constante también.

Problema de valor inicial: Atendiendo al teorema de la existencia y unicidad de Cauchy, podemos asegurar que nuestra ecuación diferencial lo cumple, puesto que la función es continua en nuestro dominio (t>0), lo que implica que admite solución, y su primera derivada al ser también continua en su dominio, admitirá una única solución.

3 ¿Cómo cambiaría la ecuación si el proceso fuera reversible?

Dado que se sigue cumpliendo la ley de acción de masas, podemos relacionar las velocidades de reacción del proceso directo e inverso mediante sus respectivas constantes específicas de velocidad de reacción k1 y k2, obteniendo finalmente como ecuación diferencial para trabajar:

                                             y'(t)=k1*(a0-y(t))*(b0-y(t))-k2*y(t)

Como comentario, podemos añadir que las constantes específicas de velocidad de reacción se relaciónan entre sí mediante la constante de equilibrio de la reacción K (K=k1/k2).

4 PVI por Euler

Resolvemos el PVI de la ecuación diferencial planteada al inicio por el método de Euler:

a0=3;%concentración inicial del reactivo a
b0=1;%concentración inicial del reactivo b
k1=1;%valor de la constante específica de la velocidad de reacción
t0=0;
tN=2;
y0=0;
h=0.1;%tamaño del paso
N=round((tN-t0)/h);%número de subintervalos
t=t0:h:tN;
y=zeros(1,N+1);%creamos un vector vacío de tamaño N
y(1)=y0;%valor inicial
for i=1:N
    y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i))
    end
%sacamos tabla de resultados
[t',y']
%gráfico
hold on
plot(t,y)
plot(t,(a0-y),'r')
plot(t,(b0-y),'g')
legend('productos','reactivo a','reactivo b','location','best');
xlabel('tiempo');
ylabel('concentracion');
hold off

Finalmente, dibujamos los resultados:

Evolución en el tiempo (0,2) de las concentraciones de los reactivos A y B y del producto C

5 Límite de las concentraciones

Hallamos el límite de la Ecuación Diferencial, y nos da 1, ya que la concentración final del elemento C debe ser el resultado de la concentración del reactivo limitante.

%DATOS DEL PROBLEMA
t0=0;
y0=0;
h=0.1;%tamaño del paso
t(1)=t0;
y(1)=y0;
i=1;
%el límite es 1, el objetivo es alcanzar 1
lim=1;
while abs(lim-y(i))>1*lim
 y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));%aplicamos euler
 t(i+1)=t(i)+h;
 i=i+1;
end
[t',y']
disp('Tiempo final:')
disp(t(end))
plot(t,y)%gráfico de resultados

%Vemos que cuando el tiempo es igual a 2, se puede aproximar al valor que obtendríamos para C con t=infinito. %Sí, es lógico lo que se observa, dado que como hemos visto en la gráfica del apartado anterior, la curva tiende a estabilizarse en el valor 1.

Evolución en el tiempo (0,2) de la concentración de C. Observamos que tiende a 1.

6 PVI con Trapecio y Runge-Kutta

Resolvemos el anterior PVI por el método del Trapecio y Runge-Kutta de orden 4, manteniendo también la resolución por el Método de Euler:

%DATOS DEL PROBLEMA
a0=3;
b0=1;
k1=1;
t0=0;
tN=2;
y0=0;
h=0.1;%número de subintervalos
N=round((tN-t0)/h);%número de pasos
t=t0:h:tN;
y=zeros(1,N+1);%inicializamos el vector y de tamaño N para el código de Euler
y(1)=y0;
z=zeros(1,N+1);%inicializamos el vector z de tamaño N para el código del trapecio
z(1)=y0;
r=zeros(1,N+1);%inicializamos evector r de tamaño N para el código de RK4
r(1)=y0;
for i=1:N
    y(i+1)=y(i)+h*(k1*(a0-y(i))*(b0-y(i))) %Euler
    z(i+1)=12-sqrt(12^2-(z(i).^2+16*y(i)+3))%Trapecio
    K1=k1*(a0-r(i))*(b0-r(i))%RK4
    K2=k1*(a0-(r(i)+K1*h*0.5))*(b0-(r(i)+K1*h*0.5))
    K3=k1*(a0-(r(i)+K2*h*0.5))*(b0-(r(i)+K2*h*0.5))
    K4=k1*(a0-(r(i)+K3*h))*(b0-(r(i)+K3*h))
    r(i+1)=r(i)+(h/6)*(K1+2*K2+2*K3+K4)
    
end
%sacamos tabla de resultados
[t',y']
[t',z']
[t',r']
%gráfico
hold on
plot(t,y,'k')%Euler
plot(t,z,'y')%Trapecio
plot(t,r)%RK4
plot(t,(a0-y),'r')
plot(t,(b0-y),'g')

legend('euler','trapecio','rk','reactivo a','reactivo b');
xlabel('tiempo');
ylabel('concentracion');
hold off

En la gráfica resultado, dibujamos las soluciones de la concentración de producto por los métodos de Euler, Trapecio y Runge-Kutta4 (además de las concentraciones de los reactivos A y B) para comparar así visualmente las soluciones por los diferentes procedimientos:

Gráfica comparativa de las soluciones de la ecuación diferencial por los métodos de Euler, Trapecio y Runge-Kutta (de orden 4).

7 Reacción consecutiva: obtención de la ecuación diferencial y PVI

Basándonos en la ley de acción de masas, y sabiendo que se produce la reacción consecutiva que se nos indica en el enunciado, podemos plantear un sistema con estas dos ecuaciones para obtener el PVI asociado.

y'1=k1*(a0-y1(i))*(b0-y1(i))

y'2=k2*(y1(i)-y2(i))

8 Resolución e interpretación del sistema por Euler y Runge-Kutta con k2=5

Haciendo uso del sistema de ecuaciones del apartado anterior; con los valores indicados en este apartado (k2=5, en los primeros 10 segundos y utilizando un paso de 0.1), y a través de los métodos de Euler y Runge-Kutta, podemos hacer una primera interpretación de la gráfica resultante. En el caso de realizar este ejercicio con un paso de 0.3, se puede observar como la gráfica se distorsiona, por lo cual llegamos a la conclusión que en los métodos utilizados, es preferible emplear pasos pequeños, dado que en los métodos explícitos con pasos grandes provocan que la gráfica se empiece a separar de la solución real.

centro

centro

9 Resolución e interpretación del sistema por Euler y Runge-Kutta con k2=1/5

Haciendo uso del sistema de ecuaciones del apartado 6; con los valores indicados en este apartado (k2=1/5, en los primeros 10 segundos y utilizando un paso de 0.1) y a través de los métodos de Euler y Runge-Kutta, podemos hacer una primera interpretación de la gráfica resultante.

centro