Estas son las simulaciones presentes en el artículo Isoconversional Kinetic Analysis of Complex Solid State Processes: Parallel and Succesive Reactions (URL). El artículo presenta dos simulaciones diferenciales de ejemplo, que son resueltas mediante el método Runge-Kutta. Para realizar estas simulaciones realizamos nuestros imports:
In [1]:
from __future__ import division
import numpy as np
from scipy.constants import R
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
En primer lugar definimos el modelo de primer orden y las ecuaciones diferenciales que lo emplean:
$$ \begin{align*} \frac{d\alpha_{1}}{dT} &= \frac{A_{1}}{\beta} e ^ \frac{-E_{1}}{RT} f_{1}(1- \alpha_{1}) \\ \frac{d\alpha_{2}}{dT} &= \frac{A_{2}}{\beta} e ^ \frac{-E_{2}}{RT} f_{2}(\alpha_{1} - \alpha_{2}) \\ \frac{d\alpha}{dT} &= \frac{1}{2} (\frac{d\alpha_{1}}{dT} + \frac{d\alpha_{2}}{dT}) \\ \end{align*} $$Que se programan como:
In [2]:
def F1(a):
return 1 - a
def diffeq(y, T, b,
E1 = 58.5 * 1000,
A1 = 900 * 60,
E2 = 125.4 * 1000,
A2 = 5 * 10**8 * 60):
a1, a2, a3 = y
k1 = A1/b * math.exp(-E1/(R*T))
k2 = A2/b * math.exp(-E2/(R*T))
da1 = k1 * F1(a1)
da2 = k2 * (1 - F1(a1) - a2)
return da1, da2, (da1+da2)/2
Definimos el rango de temperatura y ejecutamos la simulación junto con la gráfica.
In [3]:
rates = [2, 4, 8, 16]
temps = np.arange(400, 700, 1)
y0 = [1e-3, 1e-3, 1e-3]
curves = np.zeros((len(rates), len(temps)))
for n, b in enumerate(rates):
out = odeint(diffeq, y0, temps, args=(b,))
curves[n] = out[:, 2]
plt.plot(temps, out[:, 2])
Esto se puede hacer de igual manera con las tasas [10, 20, 30, 40], así que se omitirá. De los datos anteriores ahora obtenemos las derivadas:
In [4]:
for c in curves:
plt.plot(temps, np.gradient(c))
Aquí tengo algunas dudas, pero vamos a poner esto tal cual.
$$ \begin{align*} \frac{d\alpha_{1}}{dT} &= \frac{A_{1}}{\beta} e ^ \frac{-E_{1}}{RT} f_{1}(\alpha_{1}) \\ \frac{d\alpha_{2}}{dT} &= \frac{A_{2}}{\beta} e ^ \frac{-E_{2}}{RT} f_{2}(\alpha_{2}) \\ \end{align*} $$donde se cumple la igualdad $l_{1} + l_{2} = 1$ (la fracción de contribución del primer proceso y el segundo) y $\alpha = l_{1} \alpha_{1} + l_{2} \alpha_{2}$ para el resultado final de la reacción. Ambos valores de $l$ son asumidos como 0.5.
In [5]:
def f1(a):
out = (1-a) * (-np.log(1-a))**(1/2)
if out < 1:
return out
else:
return 0
def f2(a):
out = 1.5 * (1-a)**(2/3) / (1-(1-a)**(1/3))
return out
def diffeq2(y, T, b,
E1 = 111 * 1000,
A1 = 2.2 * 10**10,
E2 = 202 * 1000,
A2 = 7.2 * 10**16):
#print T, y
a1, a2, a3 = y
k1 = A1/b * np.exp(-E1/(R*T))
k2 = A2/b * np.exp(-E2/(R*T))
da1 = f1(a1) * k1
da2 = f2(a2) * k2
return da1, da2, da1/2 + da2/2
En caso de que se alcance el valor 1 en la fracción transformada antes de llegar al final de las temperaturas aparecerán valores NaN antes de números sin sentido. Por eso, la función find_nan encuentra la posición del primer NaN, para así saber hasta donde tenemos que llegar en la gráfica. De este modo nos aseguramos igualmente que el límite superior de la fracción transformada sea 1.
In [6]:
def find_nan(array):
for n, v in enumerate(array):
if math.isnan(v):
return n
return False
y0 = [1e-3, 1e-3, 1e-3]
rates = [2, 20]
temps = np.arange(450, 650, 1)
plots = []
for b in rates:
out = odeint(diffeq2, y0, temps, args=(b,))
first_nan = find_nan(out[:, 1])
if first_nan:
plots.append(plt.plot(temps[:first_nan], out[:, 2][:first_nan])[0])
else:
plots.append(plt.plot(temps, out[:, 2])[0])
plt.legend(plots, map(str, rates), loc='upper left')
Out[6]: