Diferencia entre revisiones de «Reacciones Complejas Grupo11C»

Revisión actual del 14:38 6 mar 2015

1 Introducción

El fin de esta exposición es el planteamiento de los resultados obtenidos en el estudio de la ecuación diferencial de la velocidad de reacción entre dos sustancias químicas mediante el uso de distintos métodos numéricos para aproximar y representar dichas soluciones. 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 para formar un producto C mediante una reacción bimolecular, es decir, para formar una molécula de C reaccionan una molécula de A con una de B.

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

La ley de acción de masas establece que: "La velocidad de una reacción, es directamente proporcional al producto de los moles por litro (concentración molar) de cada uno de los reactantes, elevadas a una potencia igual a su coeficiente estequiométricos y multiplicadas por una constante (k) de proporcionalidad, cuyo valor, depende de la naturaleza química de los reactantes y de la temperatura". La forma general de representar a una reacción química es la siguiente:

aA + bB → cC

En la expresión anterior, las letras minúsculas representan los coeficientes estequiométricos (número de moles) de cada una de las sustancias que participan en la reacción mientras que, las letras mayúsculas representan las fórmulas de los reactantes (A, B) y productos (C). En este caso la reacción es bimolecular, por lo que las letras minúsculas desaparecen.

Con base en lo anterior podemos entonces establecer la expresión matemática para la ley de acción de masas en una reacción química como una ecuacion 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 al menos una solución, y su primera derivada al ser también continua en su dominio, admitirá una única solución. Suponiendo la reacción dada en la que dos reactivos A y B reaccionan para formar C, podremos describir ésta situación mediante el siguiente problema de valor inicial (para un dominio t>0):

[math]
   \left . 
      \begin{matrix} 
         y'(t)={k_{1}}({a_{0}}-y(t))({b_{0}}-y(t))=0\\
         y(0)=0
      \end{matrix}
   \right \}
[/math]

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

En una reacción reversible existen dos reacciones químicas diferentes, la primera es la que se lleva a cabo entre los reactivos para formar a los productos, la segunda en la cual los productos reaccionan entre sí, una vez que se han generado, para formar de nuevo a los reactivos. A la primera reacción se le llama reacción directa y a la segunda, reacción inversa. 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)

La velocidad de creación de C será la velocidad con que se genera el compuesto C por parte de A y B restándole la velocidad de destrucción de C. Las constantes son diferentes dado que su valor, tal y como se establece en la ley de acción de masas, depende de la temperatura y de la naturaleza o propiedades químicas de las sustancias que reaccionan. Las propiedades químicas de los reactantes son diferentes a la de los productos. Es utilizada la relación entre ambas mediante la constante de equilibrio de la reacción K (K=k1/k2).

4 PVI por Euler

Suponiendo ahora que la reacción parte con cantidades iniciales de a0= 3 mol/l,b0= 1 mol/l y la constante especifica es k1= 1 mol/s se procede a realizar una primera aproximación del PVI mediante métodos numéricos, esta vez por el método de Euler, eligiendo un paso h= 0:1, en los primeros 2 segundos. El código de Matlab resultante tras la aplicación numérica de este método a nuestro problema de valor inicial es el siguiente:

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

Cuyo resultado queda representado en la siguiente gráfica mostrando las concentraciones en cada momento de los elementos de la reacción con los productos de azul, el reactivo A de rojo y el reactivo B en verde:

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

Por ejemplo, cuando una reacción inicial, en el momento en que se conjugan todas las condiciones (concentración de reactantes, temperatura, presión, catalizadores, etcétera), para que la reacción se efectúe, en ese preciso momento, la concentración o cantidad de productos es igual a cero y la de los reactivos igual a la cantidad que se tenga para la reacción. Conforme pasa el tiempo, la concentración de reactantes empieza a disminuir mientras que, la de los productos aumenta. Se diferencia perfectamente que el reactivo limitante de esta reacción química es el B, que disponía de una cantidad inicial menor que el A y a los dos segundos se ha acabado. También se observa a simple vista que la reacción es bimolecular, ya que cada reactivo pierde una molécula para producir una molécula del producto C. Éste compuesto parte de concentración inicial nula, produciéndose a partir de las concentraciones iniciales de A y B hasta el momento se acaba alguno.

5 Límite de las concentraciones

Para asegurarnos del comportamiento de la reacción tras largos periodos de tiempo, procedemos a calcular el límite de la ecuación de la concentración de C cuando el tiempo tiende a infinito. 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

El resultado es 1, siendo esto bastante lógico dado que al ser una reacción bimolecular la concentración final del elemento C debe ser el resultado de la concentración del reactivo limitante. Esto significa que pasados dos segundos, con la desaparición del compuesto C comienza el equilibrio químico en el recipiente al quedar solo compuesto de A y C que no reaccionan entre sí. Por tanto, 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 → ∞.

Gráfico 2. 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

Para resolver nuestro problema de valor inicial utilizaremos otros métodos para realizar aproximaciones distintas. Dado que elmétodo del Trapecio es un método implícito, debemos despejar la componente [math]z_{n+1}[/math]: [math]z_{n+1}=12-sqrt(12^2-z_{n}^2+16∙z{n}+3)[/math] A continuación se presenta el programa informático utilizado para aproximar mediante el método del trapecio y el método de Runge-Kutta de orden 4, en la gráfica mantendremos también la resolución por el método de Euler para poder compararlas de manera sencilla:

%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*z(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

El Método de Runge-Kutta es una forma numérica de resolver problemas de valor inicial que trabaja de la misma manera que lo hace el Método de Euler, dando soluciones aproximadas al problema de valor inicial propuesto. Por ello, las gráficas de ambos métodos serán muy similares, pero sabiendo de antemano que el Método de Runge-Kutta consigue hacer una mejor aproximación de las soluciones de la que puede hacer Euler. En la gráfica adjunta, quedan plasmadas 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 poder comprobar visualmente las diferencias de las soluciones por los diferentes procedimientos:

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

Se puede apreciar que el método de Runge-Kutta es el mas preciso de todos, dando un resultado muy similar al de Euler, distanciándose ambos notablemente del resultado obtenido por el método del trapecio.

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:

                                                          A + B → C→ D

podemos plantear un sistema con estas dos ecuaciones para obtener el PVI asociado.

[math] \left . \begin{matrix} y'_1(i)={k_{1}}({a_{0}}-y_1(i))({b_{0}}-y_1(i))\\ y'_2(i)={k_{2}}(y_1(i)-y_2(i))\\ y_1(0)=0\\y_2(0)=0 \end{matrix} \right \} [/math]

Observamos que la primera es la misma que obtuvimos anteriormente, ya que los compuestos A y B se siguen descomponiendo independientemente de lo que haya al otro lado de la igualdad, solo depende de la concentración de ambos compuestos. Y la segunda ecuación del sistema es de carácter parecido solo que, la concentración de C mas D es la antigua concentración de C cuando no era consecutiva la reacción. k2 hace referencia a una nueva constante específica de velocidad de reacción, independiente de la primera ya que trata otra reacción química totalmente distinta, la de disgregación de C, que no tiene nada que ver con A y B.

Por otro lado, un posible PVI con unos valores tN,k1 y k2 genéricos (ponemos los datos del siguiente ejercicio) sería:

t0=0;%tiempo inicial
tN=10;%tiempo final
h=0.1;%tamaño del paso
N=(tN-t0)/h; %número de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
y1=zeros(1,N+1); %primera incógnita
y2=zeros(1,N+1); % segunda incógnita
y1(1)=0;
y2(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
end
hold on 
plot(t,y1-y2,'r')
plot(t,y2)
hold off
%sería un tN,k1 y k2 genéricos. ponemos los valores del siguiente ejercicio.

A continuación tantearemos que ocurre cuando modificamos las constantes específicas del sistema.

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, y utilizando los siguientes valores: k2=5, lo que supone una segunda reacción mucho más rápida que la primera, t= 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 representación grafica de los compuestos C y D mediante el siguiente método numérico:

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

Gráfico 4. 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).

Como se puede apreciar, el compuesto C comienza con gran velocidad de creación. En el momento en que empieza a haber C se comienza a generar D, que como tiene una velocidad de creación aún mayor, en cuanto pasa un poco de tiempo el compuesto D termina por hacer desaparecer al C. 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.

Gráfico 5. 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 de la reacción consecutiva pero esta vez con los valores k2=1/5 es decir una segunda reacción mucho más lenta que la primera, en los primeros 10 segundos y utilizando un paso de 0.1. Utilizando para ello 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

Gráfico 6. 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 Gráfico 4, la diferencia es la reducción de la constante específica de la velocidad de la segunda reacción (k2), y observamos que: .La velocidad de la reacción de formación de D también se reduce, 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. Las concentraciones de C y D varían más lentamente, por lo que la reacción tiende a estabilizarse más tarde.