This is a quick note on how to do single-model simulations with solid-state kinetic models using Python. Isothermal simulation is done rearranging the isothermal rate law, the non-isothermal case is solved using an optimization method with scipy's optimize submodule (Levenberg-Mardquart). The Avrami-Erofeev A2 kinetic model will be used for both simulations, setting $E_{a} = 100 kJ/mol$ and $A = 10^{13} min^{-1}$.
We first do all required imports and define our kinetic triplet (model, A, E):
In [31]:
from __future__ import division
import time
import functools
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
from scipy.constants import R
E = 100 * 1000 # kJ/mol
A = 10**13 # min ** -1
def iA2(a):
"""Integral form of Avrami-Erofeev model (n=2)"""
return (-np.log(1-a))**(1/2)
The integral isothermal rate law is expressed as:
$$ \begin{align*} g(\alpha) = A e ^ {-\frac{E_{a}}{RT}} t \\ \end{align*} $$which can be rearranged as:
$$ \begin{align*} t = \frac{g(\alpha)}{A e ^ {-\frac{E_{a}}{RT}}} \end{align*} $$Accordingly, we can get the time associated with each $\alpha$ value by setting a model $g(\alpha)$ and kinetic parameters $E_{a}$ and $A$ as constants.
In [32]:
alphas = np.arange(0.01, 1, 0.01) # Transformed fraction, numbers between 0 and 1
temps = [340, 345, 350, 355, 360] # Simulated temperatures, in Kelvins
times = []
for T in temps:
curve = iA2(alphas) / A * np.exp(E/(R*T))
times.append(curve)
plots = [plt.plot(t, alpha)[0] for t in times]
plt.title('Simulated Isothermal Reaction')
plt.xlabel('Time (min)')
plt.ylabel('Transformed Fraction')
plt.legend(plots, map(str, temps), loc='lower right')
Out[32]:
As presented on this paper by Kwawam & Flanagan, non-isothermal data can be simulated with an optimization method by calculating the temperature (T) for transformed fraction values ($\alpha$) according to:
$$ \begin{align*} \psi = \mid {g(\alpha) \frac{\beta R}{A E_{a}} - p(x)} \mid \end{align*} $$where A, E, $g(\alpha)$ are fixed constants, the heating rate $\beta$ changes with every curve, and the exponential integral p(x)
$$ \begin{align*} p(x) = \int_{x}^{\infty} \frac{e^{-x}}{x^2} dx \end{align*} $$where $x = \frac{E}{RT}$, can be approximated using a 3rd degree Senum-Yang polynomial. In the following code $\psi$ is expressed as an objetive function where the keyword argument $g$ has to be set up explicitly (avoiding the need to make multiple psi functions). This can be done multiple times, or usign functools.partial as depicted here.
In [33]:
def senumyang(x, d=3):
"""Senum-Yang 3rd degree temperature integral approximation. x = Ea / (R*T)"""
return np.exp(-x) / x * (x**2 + 10*x + 18) / (x**3 + 12*x**2 + 368*x + 24)
def psi(T, alpha, beta, A, Ea, g=None):
"""Objetive function for nonisothermal simulation."""
x = Ea/(R*T)
return np.abs((g(alpha) * (beta*R)/(A*Ea)) - senumyang(x))
# we are using the same alpha values from the isothermal simulation
rates = [1, 2, 4, 8, 16] # Heating Rate in Kelvin Degrees
guess = 300 # Initial guess for local optimization. Your mileage may vary.
ofunc = functools.partial(psi, g=iA2)
start = time.time()
temperatures = []
for b in rates:
curve = [leastsq(ofunc, guess, args=(a, b, A, E))[0][0] for a in alphas]
temperatures.append(curve)
print "Done all curves in", time.time()-start, "seconds"
plots = [plt.plot(T, alpha)[0] for T in temperatures]
plt.title('Simulated Non-isothermal Reaction')
plt.xlabel('Temperature (K)')
plt.ylabel('Transformed Fraction')
plt.legend(plots, map(str, rates), loc='upper left')
Out[33]:
In [33]: