Revisión actual del 22:58 3 mar 2015

1 Introducción y modelización del problema.

En este artículo vamos a examinar el comportamiento cinético de las reacciones químicas, estudiaremos mediante análisis numérico las reacciones químicas irreversibles y nos limitaremos a plantear una suposición teórica del comportamiento de las las reacciones reversibles.

Para realizar este estudio, formularemos las siguientes hipótesis sobre las condiciones y el comportamiento de las reacciones químicas:

• Los componentes de la reacción química se considerarán perfectamente mezclados.
• La reacción química se realiza en condiciones de volumen y temperatura constantes.
• Consideraremos la reacción química que planteamos perfectamente ajustada.
• Suponemos que para el instante [math]t=0[/math] la concentración de producto es nula.

Partiendo de estas hipótesis y suposiciones sobre las reacciones que vamos a estudiar, consideraremos el proceso químico:: [math]{A+B}\Rightarrow{C}[/math]

Siendo esta una reacción bimolecular, en la que los reactivos [math]A[/math] y [math]B[/math] interaccionan formando un producto [math]C[/math]. Regiremos nuestro estudio cinético de esta reacción por 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. Siguiendo esta ley podemos deducir la siguiente ecuación diferencial que nos define las cantidades de producto y reactivos a lo largo del tiempo::

[math]y'(t)=k_1(a_0-y(t))(b_0-y(t))[/math]

Donde:

• [math]y(t)[/math] nos define como evoluciona la concentración de producto [math]C[/math] a lo largo del tiempo.
• [math]y'(t)[/math] es la derivada de la concentración de [math]C[/math] con respecto al tiempo.
• [math]a_0[/math] indica la concentración inicial de producto [math]A[/math].
• [math]b_0[/math] indica la concentración inicial de producto [math]B[/math].
• [math]k_1[/math] representa la constante de proporcionalidad de la reacción.

Podemos comprobar que al estar la estequiometría de la reacción ajustada podemos decir que para una concentración [math]y(t)[/math] en cada instante de tiempo [math]t[/math], tendremos concentraciones [math]a_0-y(t)[/math] y [math]b_0-y(t)[/math] de [math]A[/math] y [math]B[/math] respectivamente. Entonces si la velocidad de reacción, [math]y'(t)[/math], es directamente proporcional a la concentración de los reactivos, suponiendo una constante de proporcionalidad [math]k_1[/math] y una concentración inicial nula de producto para [math]t=0[/math], siguiendo la ley de acción de masas nos lleva a la ecuación predicha anteriormente::

[math]y'(t)=k_1(a_0-y(t))(b_0-y(t))[/math]

Partiendo de esta ecuación podemos empezar a hacer análisis numérico del proceso.

2 Proceso reversible.

En el caso de suponer un proceso químico reversible veríamos afectada la ecuación diferencial que nos cuantifica la concentración de producto en un término reductor de esta, que nos representaría la velocidad de formación de reactivos, proporcional a la concentración actual de producto por otra constante de proporcionalidad diferente [math]k_2[/math].:

[math]y'(t)=k_1(a_0-y(t))(b_0-y(t))-k_{-1}(y(t))[/math]

3 Casos prácticos. Aplicación del análisis numérico.

Para realizar el análisis numérico del proceso, cuantificaremos las constantes que dan forma a la ecuación cinética:

• [math]k_1=1[/math]
• [math]a_0=3 mol/l[/math]
• [math]b_0=1mol/l[/math]

Entrando con estos valores en la ecuación diferencial, nos queda::

[math]y'(t)=(3-y(t))(1-y(t))[/math]

Que modelizaremos siguiendo los métodos numéricos de Euler, trapecio y Runge-Kutta.

3.1 Método de Euler

Implementamos este modelo introduciendo el siguiente código en Octave:

%Establecemos los valores iniciales
t0=0;
tN=4;
c0=0;
%Discretizamos
N=40;
h=(tN-t0)/N;
%Vector de tiempos
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
c=zeros(1,N+1);
%Damos el valor para el primer elemento
c(1)=c0;
%Creamos la variable yy para operar
cc=c0;
%Creamos el bucle que nos da soluciones de y
for n=1:N
  cc=cc+h*(1*(3-cc)*(1-cc));
  c(n+1)=cc;
end
%Dibujamos la solucione
plot(t,c)
a=3.-c;
b=1.-c;
hold on
plot(t,a,'r')
hold on
plot(t,b,'y')

Mediante este código obtenemos la gráfica de las concentraciones de los reactivos y productos, que se muestra en la siguiente imagen:

Concentraciones de los reactivos y productos calculadas por el método de Euler. A se muestra en rojo, B en amarillo, y C en azul.

Como podemos ver en la imagen la concentración de los reactivos y productos se vuelve asintótica tendiendo la de [math]A=2 mol/l[/math], la de [math]B=0 mol/l[/math] y la de [math]C=1 mol/l[/math]. Lo que por la estequiometría de la reacción es lógico, ya que el producto [math]C[/math] dejará de formarse cuando se agote uno de los dos reactivos que lo forman; en nuestro caso disminuye progresivamente la concentración de [math]B[/math], al que se denomina reactivo limitante, y a la vez se estabilizan las concentraciones, ya que se reduce la velocidad de reacción [math]y'(t)[/math].

El tiempo elegido para el primer modelo del proceso químico era [math]t\in{[0,4]}[/math], en el que vemos que ya prácticamente la velocidad de reacción es nula y las concentraciones no varían. Si aumentamos el intervalo de tiempos lo suficiente como para considerarlo infinito veremos como los valores de las concentraciones se aproximan todavía más a sus valores definitivos y estables. Introduciendo el siguiente código en lenguaje Matlab, para un intervalo de tiempos [math]t\in{[0,25]}[/math] obtenemos las gráficas de las concentraciones definitivas:

%Establecemos los valores iniciales
t0=0;
tN=25;
c0=0;
%Discretizamos
N=500;
h=(tN-t0)/N;
%Vector de tiempos
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
c=zeros(1,N+1);
%Damos el valor para el primer elemento
c(1)=c0;
%Creamos la variable yy para operar
cc=c0;
%Creamos el bucle que nos da soluciones de y
for n=1:N
  cc=cc+h*(1*(3-cc)*(1-cc));
  c(n+1)=cc;
end
%Dibujamos la solucione
plot(t,c)
a=3.-c;
b=1.-c;
hold on
plot(t,a,'r')
hold on
plot(t,b,'y')

Así las concentraciones finales de los reactivos nos quedan:

Concentraciones de los reactivos y productos calculadas por el método de Euler. A se muestra en rojo, B en amarillo, y C en azul.

Como predijimos las concentraciones estacionarias nos quedan fijas, ya que la velocidad de formación de [math]C[/math], [math]y'(t)[/math], tiende a hacerse nula.

3.2 Método del trapecio

A diferencia de los otros métodos que usaremos para aproximar las soluciones de la ecuación diferencial que rige el comportamiento de nuestro modelo matemático, el esquema numérico mediante el cual calculamos es implícito, por lo que habremos de despejar en la ecuación que nos plantea el esquema::

[math]y_{n+1}=y_{n}+h\frac{f(t_n,y_n)+f(t_{n+1},y_{n+1})}{2}[/math]

Sabiendo que:: [math]f(t,y)=k_1(a_o-y(t))(b_0-y(t))[/math]

Podemos plantear el siguiente código de matlab:

clear all
%Establecemos los valores iniciales
t0=0;
tN=4;
c0=0;
k1=1;
a0=3;
b0=1;
%Discretizamos
N=500;
h=(tN-t0)/N;
%Vector de tiempos
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
c=zeros(1,N+1);
%Damos el valor para el primer elemento
c(1)=c0;
%Creamos la variable yy para operar
cc=c0;
%Creamos el bucle que nos da soluciones de y
for n=1:N
  A=-k1*h/2;
  B=1+k1*h*(a0+b0)/2;
  C=-k1*h*(a0*b0+(a0-cc)*(b0-cc))/2-cc;
  cc=(-B+sqrt(B^2-4*A*C))/(2*A);
  c(n+1)=cc;
end
%Dibujamos la solucione
plot(t,c)
a=3.-c;
b=1.-c;
hold on
plot(t,a,'r')
hold on
plot(t,b,'y')

Que nos dibuja la siguiente gráfica de las concentraciones de los reactivos y productos frente al tiempo:

Concentraciones de los reactivos y productos calculadas por el método del trapecio. A se muestra en rojo, B en amarillo, y C en azul.

Vemos que se nos mantienen los mismos resultados que vistos con el método de Euler, aunque sabemos que este método nos da una aproximación mayor.

3.3 Método de Runge-Kutta

Al igual que el método de Euler se compone de un sistema explicito, por lo que no hay que hacer ninguna manipulación al esquema numérico. Introduciendo el siguiente código Matlab implementamos el método:

clear all
%Establecemos los valores iniciales
t0=0;
tN=4;
c0=0;
k1=1;
a0=3;
b0=1;
%Discretizamos
N=500;
h=(tN-t0)/N;
%Vector de tiempos
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
c=zeros(1,N+1);
%Damos el valor para el primer elemento
c(1)=c0;
%Creamos la variable yy para operar
cc=c0;
%Creamos el bucle que nos da soluciones de y
for n=1:N
  c1=k1*(a0-cc)*(b0-cc);
  c2=k1*(a0-(cc+c1*h/2))*(b0-(cc+c1*h/2));
  c3=k1*(a0-(cc+c2*h/2))*(b0-(cc+c2*h/2));
  c4=k1*(a0-(cc+c3*h))*(b0-(cc+c3*h));
  cc=cc+h*(c1+2*c2+2*c3+c4)/6;
  c(n+1)=cc;
end
%Dibujamos la solucione
plot(t,c)
a=3.-c;
b=1.-c;
hold on
plot(t,a,'r')
hold on
plot(t,b,'y')

Obteniendo así la gráfica de las concentraciones siguiente:

Concentraciones de los reactivos y productos calculadas por el método de Runge-Kutta. A se muestra en rojo, B en amarillo, y C en azul.

4 Reacción consecutiva.

Supongamos ahora que se produce otra reacción consecutiva a la anterior para producir otro producto, [math]D[/math], proveniente del producto de la reacción anterior, [math]C[/math].:

[math]A+B\Rightarrow{C}\Rightarrow{D}[/math]

Nos regiremos de nuevo por la ley de acción de masas para determinar el sistema de ecuaciones diferenciales que modeliza el comportamiento de esta reacción química. La velocidad de reacción de este nuevo proceso químico consecutivo al anterior también será proporcional a la concentración de su reactivo, en este caso de [math]C[/math], y la velocidad de reacción del proceso anterior se verá reducida por la velocidad de formación del nuevo compuesto [math]D[/math], siendo esta negativa cuando se vaya anulando progresivamente la cantidad de uno de los dos reactivos iniciales, [math]A[/math] o [math]B[/math]. Así prodríamos modelizar el comportamiento de la reacción anteriormente descrita con el siguiente sistema de ecuaciones diferenciales::

[math]\left. y_1'(t)=k_1(a_0-y_1(t))(b_0-y_1(t))-y_2'(t) \atop y_2'(t)=k_2(y_1(t)-y_2(t)) \right\}[/math]

Donde:

• [math]y_1(t)[/math] nos marca el producto reaccionado en la primera ecuación a lo largo del tiempo.
• [math]y_2(t)[/math] es el producto [math]D[/math], obtenido de la segunda reacción.
• [math]a_0[/math] indica la concentración inicial de producto [math]A[/math].
• [math]b_0[/math] indica la concentración inicial de producto [math]B[/math].
• [math]k_1[/math] representa la constante de proporcionalidad de la primera reacción.
• [math]k_2[/math] representa la constante de proporcionalidad de la segunda reacción.

Si nos fijamos en este caso [math]y_1(t)[/math] no nos da el valor de la concentración de [math]C[/math] directamente, sino que esta será [math]y_1(t)-y_2(t)[/math], la concentración que se forma en la primera reacción menos el producto de la segunda reacción ambas a lo largo del tiempo.

Podemos observar que se trata de un sistema de ecuaciones diferenciales no lineales de primer orden. En el siguiente apartado veremos algunos casos prácticos para resolver este sistema dando valores a las constantes que necesitamos determinar en el problema.

5 Casos prácticos. Aplicación del análisis numérico.

Para realizar el análisis numérico del proceso, determinaremos dos grupos de constantes que conformarán la deriva temporal de esta modelización.

Primeramente diremos que [math]k_1\gtk_2[/math], propondremos este grupo de valores:

• [math]k_1=1[/math]
• [math]k_2=\frac{1}{5}[/math]
• [math]a_0=3 mol/l[/math]
• [math]b_0=1mol/l[/math]

Para lo que introduciendo el siguiente código en lenguaje Matlab, conseguimos una representación gráfica de las concentraciones de las especies químicas incluidas en esta reacción:

%Establecemos los valores iniciales
t0=0;
tN=10;
y0=0;
%Discretizamos
N=1000;
h=(tN-t0)/N;
%Vector de tiempos
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
y1=zeros(1,N+1);
y2=zeros(1,N+1);
%Damos el valor para el primer elemento
y1(1)=y0;
y2(1)=y0;
%Creamos la variable yy para operar
yy1=y0;
yy2=y0;
%Creamos el bucle que nos da soluciones de y
for n=1:N
  yy1=yy1+h*(1*(3-yy1)*(1-yy1)-(1/5)*(yy1-yy2));
  yy2=yy2+h*((1/5)*(yy1-yy2));
  y1(n+1)=yy1;
  y2(n+1)=yy2;
end

%Dibujamos la solucione
a=3.-y1;
b=1.-y1;
c=y1-y2;
d=y2;
plot(t,c)
hold on
plot(t,a,'g')
hold on
plot(t,b,'y')
hold on
plot(t,d,'r')

Obteniendo el siguiente gráfico:

Concentraciones de los reactivos y productos calculadas por el método de Euler. A se muestra en verde, B en amarillo, C en azul y D en rojo

Si suponemos ahora que [math]k_1\ltk_2[/math], el ejemplo propuesto tendrá las siguientes constantes:

• [math]k_1=1[/math]
• [math]k_2=5[/math]
• [math]a_0=3 mol/l[/math]
• [math]b_0=1mol/l[/math]

Así el código nos queda:

%Establecemos los valores iniciales
t0=0;
tN=10;
y0=0;
%Discretizamos
N=1000;
h=(tN-t0)/N;
%Vector de tiempos
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
y1=zeros(1,N+1);
y2=zeros(1,N+1);
%Damos el valor para el primer elemento
y1(1)=y0;
y2(1)=y0;
%Creamos la variable yy para operar
yy1=y0;
yy2=y0;
%Creamos el bucle que nos da soluciones de y
for n=1:N
  yy1=yy1+h*(1*(3-yy1)*(1-yy1)-(5)*(yy1-yy2));
  yy2=yy2+h*((5)*(yy1-yy2));
  y1(n+1)=yy1;
  y2(n+1)=yy2;
end

%Dibujamos la solucione
a=3.-y1;
b=1.-y1;
c=y1-y2;
d=y2;
plot(t,c)
hold on
plot(t,a,'g')
hold on
plot(t,b,'y')
hold on
plot(t,d,'r')

El gráfico obtenido es el siguiente:

Concentraciones de los reactivos y productos calculadas por el método de Euler. A se muestra en verde, B en amarillo, C en azul y D en rojo

5.1 Interpretación de los resultados

Podemos ver en las gráficas anteriores, dado que los únicos cambios realizados son sobre las constantes cinéticas [math]k_1[/math] y [math]k_2[/math], cuando [math]k_1[/math] es mayor que [math]k_2[/math] la primera reacción es más rápida que la segunda, por tanto se acumula más especie [math]C[/math] antes de que la velocidad de formación de [math]D[/math], entonces la concentración de [math]C[/math] y [math]D[/math] va tendiendo más lentamente a [math]0[/math] y [math]1[/math] respectivamente. Sin embargo, cuando esta relación entre [math]k_1[/math] y [math]k_2[/math] es inversa, se produce [math]D[/math] más rápido que se disipa [math]C[/math], tendiendo así antes a sus valores estacionarios.