Diferencia entre revisiones de «Reacciones Complejas Grupo11C»

Revisión del 13:48 5 mar 2015

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))

Observamos que la primera es la misma que obteníamos por la primera reacción, y la segunda es de carácter parecido solo que, al partir de un solo reactivo (que sería C) es más simple. k2 hace referencia a una nueva constante específica de velocidad de reacción, independiente de la primera.

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.

t0=0; 
tN=10;
h=0.1;
N=(tN-t0)/h; %número de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
%euler
y1=zeros(1,N+1); %primera incógnita
y2=zeros(1,N+1); % segunda incógnita
y1(1)=0;
y2(1)=0;
%RK4 
r1=zeros(1,N+1); %primera incógnita
r2=zeros(1,N+1); % segunda incógnita
r1(1)=0;
r2(1)=0;
k1=1;
k2=5;

for i=1:N
y1(i+1)=y1(i)+h*(k1*(a0-y1(i))*(b0-y1(i))) %euler primera ecuación del sistema
y2(i+1)=y2(i)+h*(k2*(y1(i)-y2(i))) %euler segunda ecuación del sistema
 K11=k1*(a0-r1(i))*(b0-r1(i))%RK4 para la primera ecuación
    K21=k1*(a0-(r1(i)+K11*h*0.5))*(b0-(r1(i)+K11*h*0.5))
    K31=k1*(a0-(r1(i)+K21*h*0.5))*(b0-(r1(i)+K21*h*0.5))
    K41=k1*(a0-(r1(i)+K31*h))*(b0-(r1(i)+K31*h))
    r1(i+1)=r1(i)+(h/6)*(K11+2*K21+2*K31+K41)
 K12=k2*(r1(i)-r2(i))%RK4 para la segunda ecuación
    K22=k2*(r1(i)-(r2(i)+K12*h*0.5))
    K32=k2*(r1(i)-(r2(i)+K22*h*0.5))
    K42=k2*(r1(i)-(r2(i)+K32*h))
    r2(i+1)=r2(i)+(h/6)*(K12+2*K22+2*K32+K42)
end
hold on 
plot(t,y1-y2,'r')
plot(t,y2)
plot(t,r1-r2,'k')
plot(t,r2,'r--')
legend('euler 1','euler 2','rk 1','rk 2');
hold off

Resolución del sistema de ecuaciones por los métodos de Euler y Runge-Kutta4 con paso 0.1 (el '1' hace referencia a la concentración de C y el '2' a la concentración de D).

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 los pasos grandes provocan que la gráfica se empiece a separar de la solución real.

t0=0; 
tN=10;
h=0.3;
N=(tN-t0)/h; %número de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
%euler
y1=zeros(1,N+1); %primera incógnita
y2=zeros(1,N+1); % segunda incógnita
y1(1)=0;
y2(1)=0;
%RK4 
r1=zeros(1,N+1); %primera incógnita
r2=zeros(1,N+1); % segunda incógnita
r1(1)=0;
r2(1)=0;
k1=1;
k2=5;

for i=1:N
y1(i+1)=y1(i)+h*(k1*(a0-y1(i))*(b0-y1(i))) %euler primera ecuación del sistema
y2(i+1)=y2(i)+h*(k2*(y1(i)-y2(i))) %euler segunda ecuación del sistema
 K11=k1*(a0-r1(i))*(b0-r1(i))%RK4 para la primera ecuación
    K21=k1*(a0-(r1(i)+K11*h*0.5))*(b0-(r1(i)+K11*h*0.5))
    K31=k1*(a0-(r1(i)+K21*h*0.5))*(b0-(r1(i)+K21*h*0.5))
    K41=k1*(a0-(r1(i)+K31*h))*(b0-(r1(i)+K31*h))
    r1(i+1)=r1(i)+(h/6)*(K11+2*K21+2*K31+K41)
 K12=k2*(r1(i)-r2(i))%RK4 para la segunda ecuación
    K22=k2*(r1(i)-(r2(i)+K12*h*0.5))
    K32=k2*(r1(i)-(r2(i)+K22*h*0.5))
    K42=k2*(r1(i)-(r2(i)+K32*h))
    r2(i+1)=r2(i)+(h/6)*(K12+2*K22+2*K32+K42)
end
hold on 
plot(t,y1-y2,'r')
plot(t,y2)
plot(t,r1-r2,'b')
plot(t,r2,'r--')
legend('euler 1','euler 2','rk 1','rk 2');
hold off
%para pasos muy grandes, el método de euler explícito empieza a funcionar mal.

Resolución del sistema de ecuaciones por los métodos de Euler y Runge-Kutta4 con paso 0.3(el '1' hace referencia a la concentración de C y el '2' a la concentración de D).

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.

t0=0; 
tN=10;
h=0.1;
N=(tN-t0)/h; %número de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
%euler
y1=zeros(1,N+1); %primera incógnita
y2=zeros(1,N+1); % segunda incógnita
y1(1)=0;
y2(1)=0;
%RK4 
r1=zeros(1,N+1); %primera incógnita
r2=zeros(1,N+1); % segunda incógnita
r1(1)=0;
r2(1)=0;
k1=1;
k2=1/5;

for i=1:N
y1(i+1)=y1(i)+h*(k1*(a0-y1(i))*(b0-y1(i))) %euler primera ecuación del sistema
y2(i+1)=y2(i)+h*(k2*(y1(i)-y2(i))) %euler segunda ecuación del sistema
 K11=k1*(a0-r1(i))*(b0-r1(i))%RK4 para la primera ecuación
    K21=k1*(a0-(r1(i)+K11*h*0.5))*(b0-(r1(i)+K11*h*0.5))
    K31=k1*(a0-(r1(i)+K21*h*0.5))*(b0-(r1(i)+K21*h*0.5))
    K41=k1*(a0-(r1(i)+K31*h))*(b0-(r1(i)+K31*h))
    r1(i+1)=r1(i)+(h/6)*(K11+2*K21+2*K31+K41)
 K12=k2*(r1(i)-r2(i))%RK4 para la segunda ecuación
    K22=k2*(r1(i)-(r2(i)+K12*h*0.5))
    K32=k2*(r1(i)-(r2(i)+K22*h*0.5))
    K42=k2*(r1(i)-(r2(i)+K32*h))
    r2(i+1)=r2(i)+(h/6)*(K12+2*K22+2*K32+K42)
end
hold on 
plot(t,y1-y2,'r')
plot(t,y2)
plot(t,r1-r2,'k')
plot(t,r2,'r--')
legend('euler 1','euler 2','rk 1','rk 2');
hold off

Evolución en el tiempo de las concentraciones de C (indicada con '1') y D (indicada con '2') y evaluadas por los métodos de Euler y Runge-Kutta.

Al comparar esta gráfica con la del ejercicio anterior observamos que al reducir la constante específica de la velocidad de la segunda reacción (k2) hacemos que su velocidad también se reduzca, y al ser ésta menor a su vez que la velocidad de reacción para crear C conseguimos obtener un pico de valor superior en la concentración de C y que la concentración de D crezca de manera más moderada.