Methane combustion¶

The main reference for this model is: Andersen, J. et al. (2009) ‘Global Combustion Mechanisms for Use in CFD Modeling under Oxy-Fuel Conditions’, Energy & fuels, 23(3), pp. 1379–1389. doi: 10.1021/ef8003619.

In [ ]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

The variables $x_0,\dotsc,x_4$ are the concentrations of the molecules $CH_4$ (methane), $O_2$ (oxygen), $CO$ (carbon monoxide), $H_2O$ (water) and $CO_2$ (carbon dioxide) respectively. The first row of the matrix $E$ gives the numbers of carbon atoms in each of these molecules, the second row gives the numbers of oxygen atoms, and the third row gives the numbers of hydrogen atoms.

In [ ]:
# This is the modified Westbrook-Dryer model as in Table 4 of Andersen et al.
# Their Table 3 contains Chemkin code and indicates that the parameters are to
# be interpreted according to Chemkin conventions.  The Chemkin manual at 
# https://www3.nd.edu/~powers/ame.60636/chemkin2000.pdf indicates (in Equation 52)
# that the reaction rate should be given by k = A * T^beta * ey(-Ea / (RT)).
# Equation 52 actually has R_c instead of R, but the footnote eylains that R_c
# is morally the same as R except that it is in cal/mol-K instead of J/mol-K.

R   = 1.987                 # Universal gas constant, cal/mol-K
P0  = 101325. / 4.184 / 1e6 # Atmospheric pressure, cal / cm^3

def molarity_from_fraction(p, T, P=P0):
    """Convert molar fractions to molarities at temperature T and pressure P."""
    return np.array(p) * P / (R * T)

def fraction_from_molarity(x):
    """Convert molarities to molar fractions."""
    return np.array(x) / np.sum(x, axis=0)

A    = np.array([1.59e13, 3.98e8, 6.16e13])
beta = np.array([0, 0, -0.97])
E    = np.array([47.8, 10.0, 78.4]) * 1e3

def rate_constants(T):
    """Calculate the rate constants for the three reactions at temperature T."""
    return A * T ** beta * np.exp(-E / (R * T))

N = np.array([[-1, 0, 0],[-1.5,-0.5, 0.5],[ 1,-1, 1],[ 2, 0, 0],[ 0, 1,-1]])
U = np.array([[1, 0, 1, 0, 1],
              [0, 2, 1, 1, 2],
              [4, 0, 0, 2, 0]])

names  = [r"$CH_4$", r"$O_2$", r"$CO$", r"$H_2O$", r"$CO_2$"]
masses = np.array([16.04, 32.00, 28.01, 18.02, 44.01]) # g/mol

def x_dot(t, x, k):
    """
    The rate of change of the molar concentrations of the species in the
    Westbrook-Dryer model.  The input x is a vector of molar concentrations
    of the species, and the output is the time derivative of x.
    The argument t is ignored, but is required by solve_ivp.
    The argument k is a vector of rate constants for the three reactions,
    as calculated by rate_constants().
    """
    y = np.maximum(x, 1e-9)
    r = np.array([k[0] * (y[0] **  0.70) * (y[1] ** 0.80),         # [CH4] ** 0.7 * [O2] ** 0.8
                  k[1] * (y[1] **  0.25) * y[2] * (y[3] ** 0.50),  # [O2] ** 0.25 * [CO] * [H2O] ** 0.5
                  k[2] * (y[1] ** -0.25) * (y[3] ** 0.50) * y[4]]) # [O2] ** -0.25 * [H2O] ** 0.5 * [CO2]
    return N @ r

def solve_combustion(T, p0, t_max=0.1, n=1000):
    k = rate_constants(T)
    x0 = molarity_from_fraction(p0, T)
    t = np.linspace(0, t_max, n)
    sol = solve_ivp(x_dot, [0, t_max], x0, args=(k,), t_eval=t, atol=1e-9, method='BDF')
    x1 = sol.y
    p1 = fraction_from_molarity(x1)
    return sol.t, p1

def show_combustion(T, p0, t_max=0.1, n=1000):
    t, p1 = solve_combustion(T, p0, t_max, n)
    for i in range(5):
        plt.plot(t, p1[i], label=names[i])
    plt.xlabel("Time (s)")
    plt.xscale("log")
    plt.ylabel("Molar fraction")
    plt.legend()
    plt.show()

The matrix product $UN$ is zero, as shown below. It follows that $d(Ux)/dt=U\dot{x}=UNr=0$, so the vector $Ux$ is constant. The reason is that $U_{ij}$ is the number of atoms of element $i$ in molecule $j$ (where the elements are $C$, $O$ and $H$, and the molecules are $CH_4$, $O_2$, $CO$, $H_2O$ and $CO_2$). Thus $(Ux)_i$ is the molar concentration of atoms of element $i$. We are doing chemistry rather than nuclear physics, so atoms cannot be created, destroyed or changed to other elements, so $(Ux)_i$ is constant.

In [ ]:
U @ N
In [ ]:
show_combustion(1400., [0.1, 0.25, 0.01, 0.01, 0.63], n=10000)
In [ ]: