For generic combustion of $C_xH_y$, products of complete combustion gives:
Lean: $$C_xH_y + (1+E)(x+y/4)(O_2 + 3.76 N_2) \rightarrow xCO_2 + \frac{y}{2}H_2O + (x+y/4)E O_2 + 3.76(x+y/4)(1+E) N_2.$$
Rich: $$\phi C_xH_y + (x+y/4)(O_2 + 3.76 N_2) \rightarrow xCO_2 + \frac{y}{2}H_2O + (\phi-1)C_xH_y + 3.76(x+y/4)(1+E) N_2.$$
When rich, include species $CO$ and $H_2$, and constrain the composition to equilibrate the water gas shift (WGS) reaction:
$$CO_2 + H_2 \rightleftharpoons CO + H_2O.$$
The rich combustion reaction is then
$$\phi C_xH_y + (x+y/4)(O_2 + 3.76 N_2) \rightarrow aCO_2 + bH_2O + cCO + dH_2 + 3.76(x+y/4)N_2.$$
For giveen $\phi$, there are four variables: $a$, $b$, $c$, $d$, and four constraints: WGS equilibrium, and balances on $C$, $H$, $O$, and $N$. (Note that the balance on $N_2$ is built in to the reaction.)
The balance equations are
These can be rearranged to solve for $c$, $b$, and $d$ in terms of $a$:
The WGS equilibrium is given as
$$K_{eq} = e^{-\Delta G^o/RT},$$
and
$$K_{eq} = \frac{n_{CO}\cdot n_{H2O}}{n_{CO2}\cdot n_{H2}} = \frac{c\cdot b}{a\cdot d}.$$
Here, $K_eq$ is the product of activities to the power of the coefficients. Activities are partial pressures divided by the reference pressure; partial pressures are pressure times mole fraction; mole fractions are moles divided by total moles. Since there is no change in moles for the WGS reaction, $K_{eq}$ simplifies as given. Inserting $b$, $c$, and $d$ gives $K_{eq}$ in terms of $a$. The result is a quadratic equation in $a$:
$$\alpha a^2 + \beta a + \gamma = 0,$$
$$\alpha = K_{eq}-1,$$ $$\beta = K_{eq}(\phi y/2 - 2(+y/4) + \phi x) + 2(x+y/4),$$ $$\gamma = -\phi x(2(x+y/4) - \phi x).$$
The solution is
$$a = \frac{-\beta \pm\sqrt{\beta^2 - 4\alpha\gamma}}{2\alpha}.$$
The first root (that is, $-\beta + \sqrt{\cdot}$) gives the physical solution. Then $b$, $c$, $d$ are recovered, from which mole fractions are found.
The results are compared below for CH4, assuming isothermal and adiabatic reaction. For the isothermal case, a constant temperature is assumed over all $\phi$ considered.
For $CH_4$ combustion, note that the WGS reaction “runs out of air” for $\phi > 4$, which is the stoichiometric point for $CO/H_2$ products:
$$C_xH_y + \frac{x}{2}(O_2 + 3.76 N_2) \rightarrow xCO + \frac{y}{2} H_2 + 3.76\frac{x}{2}N_2.$$
The corresponding mixture fraction is $f_{co}=0.1886$. Hence, WGS equilibrium products are considered for $1\le\phi\le 4$. Products of complete combustion are used for $\phi<1$.
Isothermal
import numpy as np
import matplotlib.pyplot as plt
import cantera as ct
gas = ct.Solution('gri30.yaml')
def getKeq(T):
gas.TPX = T, P, "CO2:1"
hco2 = gas.partial_molar_enthalpies[gas.species_index("CO2")]
sco2 = gas.partial_molar_entropies[gas.species_index("CO2")]
gco2 = hco2 - T*sco2
gas.TPX = T, P, "CO:1"
hco = gas.partial_molar_enthalpies[gas.species_index("CO")]
sco = gas.partial_molar_entropies[gas.species_index("CO")]
gco = hco - T*sco
gas.TPX = T, P, "H2:1"
hh2 = gas.partial_molar_enthalpies[gas.species_index("H2")]
sh2 = gas.partial_molar_entropies[gas.species_index("H2")]
gh2 = hh2 - T*sh2
gas.TPX = T, P, "H2O:1"
hh2o = gas.partial_molar_enthalpies[gas.species_index("H2O")]
sh2o = gas.partial_molar_entropies[gas.species_index("H2O")]
gh2o = hh2o - T*sh2o
dg = gco + gh2o - gco2 - gh2
return np.exp(-dg/ct.gas_constant/T)
phi = 1.5
x = 1.0 # CH4 is hardcoded below; if you change x, y, then change the fuel species below
y = 4.0 # CH4 is hardcoded below; if you change x, y, then change the fuel species below
P = 101325.0
T = 2000.0 # CHANGEME: make this 1000 to compare to adiabatic where the departure at phi>3 occurs.
gas.TPX = 300.0, P, "O2:2,N2:7.52"
h0 = gas.enthalpy_mass
gas.TPX = 300.0, P, "CH4:1"
h1 = gas.enthalpy_mass
fst = (12.011*x + 1.0079*y) / (12.011*x + 1.0079*y + (x+y/4)*4.76*29)
fco = (12.011*x + 1.0079*y) / (12.011*x + 1.0079*y + x/2*4.76*29)
phico = fco*(1-fst)/fst/(1-fco)
def getPCC(phi):
if phi <= 1.0:
E = 1.0/phi - 1.0
nco2 = x
nh2o = y/2
no2 = (x+y/4)*E
nn2 = (x+y/4)*(1+E)*3.76
nt = nco2+nh2o+no2+nn2
return nco2/nt, nh2o/nt, no2/nt, nn2/nt
else:
nco2 = x
nh2o = y/2
ncxhy = (phi-1.0)
nn2 = 3.76*(x+y/4)
nt = nco2+nh2o+ncxhy+nn2
return nco2/nt, nh2o/nt, ncxhy/nt, nn2/nt
def getWGS(phi,T):
#Keq = 4.0
Keq = getKeq(T)
alpha = Keq - 1.0
beta = Keq*(phi*y/2 - 2*(x+y/4) + phi*x) + 2*(x+y/4)
gamma = -phi*x*(2*(x+y/4)-phi*x)
a1 = (-beta + np.sqrt(beta*beta - 4*alpha*gamma))/(2*alpha)
a2 = (-beta - np.sqrt(beta*beta - 4*alpha*gamma))/(2*alpha)
#print(a1,a2)
a = np.maximum(a1, a2) # gives the positive root (moles should be positive)
a = a1
b = 2*(x+y/4)-phi*x-a
c = phi*x-a
d = phi*y/2 - 2*(x+y/4) + phi*x + a
n2 = 3.76*(x+y/4)
ntot = a+b+c+d+n2
xco2 = a/ntot
xh2o = b/ntot
xco = c/ntot
xh2 = d/ntot
xn2 = n2/ntot
return xco2, xh2o, xco, xh2, xn2
#-------------------------
nphi = 100
phi = np.linspace(0.01,phico,nphi)
xco2 = np.zeros_like(phi)
xco = np.zeros_like(phi)
xh2 = np.zeros_like(phi)
xh2o = np.zeros_like(phi)
xo2 = np.zeros_like(phi)
xn2 = np.zeros_like(phi)
xco2eq = np.zeros_like(phi)
xcoeq = np.zeros_like(phi)
xh2eq = np.zeros_like(phi)
xh2oeq = np.zeros_like(phi)
xo2eq = np.zeros_like(phi)
xn2eq = np.zeros_like(phi)
for i,p in enumerate(phi):
f = p*fst/(1.0-fst+p*fst)
h = h0*(1.0-f) + h1*f
if p <= 1.0:
xco2[i], xh2o[i], xo2[i], xn2[i] = getPCC(p)
xco[i] = xh2[i] = 0.0
gas.TPX = T,P,f'CO2:{xco2[i]}, H2O:{xh2o[i]}, O2:{xo2[i]}, N2:{xn2[i]}'
gas.equilibrate('TP')
#gas.HPX = h,P,f'CO2:{xco2[i]}, H2O:{xh2o[i]}, O2:{xo2[i]}, N2:{xn2[i]}'
#gas.equilibrate('HP')
xco2eq[i] = gas.X[gas.species_index('CO2')]
xcoeq[i] = gas.X[gas.species_index('CO')]
xh2eq[i] = gas.X[gas.species_index('H2')]
xh2oeq[i] = gas.X[gas.species_index('H2O')]
xo2eq[i] = gas.X[gas.species_index('O2')]
xn2eq[i] = gas.X[gas.species_index('N2')]
else:
co2, h2o, cxhy, n2 = getPCC(p)
#gas.HPX = h,P,f'CO2:{co2}, H2O:{h2o}, N2:{n2}, CH4:{cxhy}'
#gas.equilibrate('HP')
gas.TPX = T,P,f'CO2:{co2}, H2O:{h2o}, N2:{n2}, CH4:{cxhy}'
gas.equilibrate('TP')
xco2eq[i] = gas.X[gas.species_index('CO2')]
xcoeq[i] = gas.X[gas.species_index('CO')]
xh2eq[i] = gas.X[gas.species_index('H2')]
xh2oeq[i] = gas.X[gas.species_index('H2O')]
xo2eq[i] = gas.X[gas.species_index('O2')]
xn2eq[i] = gas.X[gas.species_index('N2')]
xco2[i], xh2o[i], xco[i], xh2[i], xn2[i] = getWGS(p,gas.T)
xo2[i] = 0.0
plt.plot(phi, xco2, label='CO2')
plt.plot(phi, xco, label='CO')
plt.plot(phi, xh2o, label='H2O')
plt.plot(phi, xh2, label='H2')
plt.plot(phi, xo2, label='O2')
plt.plot(phi, xn2, label='N2')
plt.plot(phi, xco2eq, ':', label='EQ')
plt.plot(phi, xcoeq, ':')
plt.plot(phi, xh2oeq, ':')
plt.plot(phi, xh2eq, ':')
plt.plot(phi, xo2eq, ':')
plt.plot(phi, xn2eq, ':')
plt.xlabel(r'$\phi$')
plt.ylabel(r'$X_i$')
plt.legend(frameon=False);
Adiabatic
phi = 1.5
x = 1.0
y = 4.0
P = 101325.0
T = 2000.0
gas.TPX = 300.0, P, "O2:2,N2:7.52"
h0 = gas.enthalpy_mass
gas.TPX = 300.0, P, "CH4:1"
h1 = gas.enthalpy_mass
fst = (12.011*x + 1.0079*y) / (12.011*x + 1.0079*y + (x+y/4)*4.76*29)
fco = (12.011*x + 1.0079*y) / (12.011*x + 1.0079*y + x/2*4.76*29)
phico = fco*(1-fst)/fst/(1-fco)
def getPCC(phi):
if phi <= 1.0:
E = 1.0/phi - 1.0
nco2 = x
nh2o = y/2
no2 = (x+y/4)*E
nn2 = (x+y/4)*(1+E)*3.76
nt = nco2+nh2o+no2+nn2
return nco2/nt, nh2o/nt, no2/nt, nn2/nt
else:
nco2 = x
nh2o = y/2
ncxhy = (phi-1.0)
nn2 = 3.76*(x+y/4)
nt = nco2+nh2o+ncxhy+nn2
return nco2/nt, nh2o/nt, ncxhy/nt, nn2/nt
def getWGS(phi,T):
#Keq = 4.0
Keq = getKeq(T)
alpha = Keq - 1.0
beta = Keq*(phi*y/2 - 2*(x+y/4) + phi*x) + 2*(x+y/4)
gamma = -phi*x*(2*(x+y/4)-phi*x)
a1 = (-beta + np.sqrt(beta*beta - 4*alpha*gamma))/(2*alpha)
a2 = (-beta - np.sqrt(beta*beta - 4*alpha*gamma))/(2*alpha)
#print(a1,a2)
a = np.maximum(a1, a2) # gives the positive root (moles should be positive)
a = a1
b = 2*(x+y/4)-phi*x-a
c = phi*x-a
d = phi*y/2 - 2*(x+y/4) + phi*x + a
n2 = 3.76*(x+y/4)
ntot = a+b+c+d+n2
xco2 = a/ntot
xh2o = b/ntot
xco = c/ntot
xh2 = d/ntot
xn2 = n2/ntot
return xco2, xh2o, xco, xh2, xn2
#-------------------------
nphi = 100
phi = np.linspace(0.01,phico,nphi)
xco2 = np.zeros_like(phi)
xco = np.zeros_like(phi)
xh2 = np.zeros_like(phi)
xh2o = np.zeros_like(phi)
xo2 = np.zeros_like(phi)
xn2 = np.zeros_like(phi)
xco2eq = np.zeros_like(phi)
xcoeq = np.zeros_like(phi)
xh2eq = np.zeros_like(phi)
xh2oeq = np.zeros_like(phi)
xo2eq = np.zeros_like(phi)
xn2eq = np.zeros_like(phi)
TT = np.zeros_like(phi)
for i,p in enumerate(phi):
f = p*fst/(1.0-fst+p*fst)
h = h0*(1.0-f) + h1*f
if p <= 1.0:
xco2[i], xh2o[i], xo2[i], xn2[i] = getPCC(p)
xco[i] = xh2[i] = 0.0
#gas.TPX = T,P,f'CO2:{xco2[i]}, H2O:{xh2o[i]}, O2:{xo2[i]}, N2:{xn2[i]}'
#gas.equilibrate('TP')
gas.HPX = h,P,f'CO2:{xco2[i]}, H2O:{xh2o[i]}, O2:{xo2[i]}, N2:{xn2[i]}'
gas.equilibrate('HP')
xco2eq[i] = gas.X[gas.species_index('CO2')]
xcoeq[i] = gas.X[gas.species_index('CO')]
xh2eq[i] = gas.X[gas.species_index('H2')]
xh2oeq[i] = gas.X[gas.species_index('H2O')]
xo2eq[i] = gas.X[gas.species_index('O2')]
xn2eq[i] = gas.X[gas.species_index('N2')]
TT[i] = gas.T
else:
co2, h2o, cxhy, n2 = getPCC(p)
gas.HPX = h,P,f'CO2:{co2}, H2O:{h2o}, N2:{n2}, CH4:{cxhy}'
gas.equilibrate('HP')
#gas.TPX = T,P,f'CO2:{co2}, H2O:{h2o}, N2:{n2}, CH4:{cxhy}'
#gas.equilibrate('TP')
xco2eq[i] = gas.X[gas.species_index('CO2')]
xcoeq[i] = gas.X[gas.species_index('CO')]
xh2eq[i] = gas.X[gas.species_index('H2')]
xh2oeq[i] = gas.X[gas.species_index('H2O')]
xo2eq[i] = gas.X[gas.species_index('O2')]
xn2eq[i] = gas.X[gas.species_index('N2')]
xco2[i], xh2o[i], xco[i], xh2[i], xn2[i] = getWGS(p,gas.T)
xo2[i] = 0.0
for iter in range(10): # simple fixed point iteration to converge T
gas.HPX = h,P, f"CO2:{xco2[i]}, H2O:{xh2o[i]}, CO:{xco[i]}, H2:{xh2[i]}, N2:{xn2[i]}"
xco2[i], xh2o[i], xco[i], xh2[i], xn2[i] = getWGS(p,gas.T)
xo2[i] = 0.0
TT[i] = gas.T
plt.plot(phi, xco2, label='CO2')
plt.plot(phi, xco, label='CO')
plt.plot(phi, xh2o, label='H2O')
plt.plot(phi, xh2, label='H2')
plt.plot(phi, xo2, label='O2')
plt.plot(phi, xn2, label='N2')
plt.plot(phi, xco2eq, ':', label='EQ')
plt.plot(phi, xcoeq, ':')
plt.plot(phi, xh2oeq, ':')
plt.plot(phi, xh2eq, ':')
plt.plot(phi, xo2eq, ':')
plt.plot(phi, xn2eq, ':')
plt.xlabel(r'$\phi$')
plt.ylabel(r'$X_i$')
plt.legend(frameon=False);
plt.plot(phi,TT)
plt.xlabel(r'$\phi$')
plt.ylabel(r'$T$ (K)');
Notes
The isothermal results are nearly identical for $T=1000$ K. The adiabatic results depart from equilibrium at $\phi\gtrsim 3$. As $\phi$ increases, the temperature decreases, and WGS equilibrium becomes a poorer approximation to the equilibrium solution considering all products. At $\phi=3$, $T\approx 1000$ K, and the adiabatic results are very similar to the isothermal results with that same temperature.
Plot adiabatic equilibrium over all mixture fractions
Here, we plot the mass fractions instead of mole fractions. Note that there is no simple profile, such as linear, for $\phi>4$, $f>0.188$. $CO$ and $H_2$ are linear but drop to zero much earlier than $\xi=1$; $H_2O$ and $CO_2$ are nonlinear.
gas.TPX = 300, P, 'O2:1, N2:3.76'
h0 = gas.enthalpy_mass
y0 = gas.Y
gas.TPX = 300, P, 'CH4:1'
h1 = gas.enthalpy_mass
y1 = gas.Y
nf = 1000
f = np.linspace(0,1,nf)
Y = np.zeros((nf,gas.n_species))
for i in range(nf):
hh = h0*(1-f[i]) + h1*f[i]
yy = y0*(1-f[i]) + y1*f[i]
gas.HPY = hh, P, yy
gas.equilibrate("HP")
Y[i,:] = gas.Y
plt.plot(f, Y[:,gas.species_index("CO2")], label='CO2')
plt.plot(f, Y[:,gas.species_index("CO")], label='CO')
plt.plot(f, Y[:,gas.species_index("H2O")], label='H2O')
plt.plot(f, Y[:,gas.species_index("H2")], label='H2')
plt.plot(f, Y[:,gas.species_index("O2")], label='O2')
#plt.plot(f, Y[:,gas.species_index("N2")], label='N2')
plt.plot([fco, fco], [0,1], ':', color='gray')
plt.ylim([0, 0.3])
plt.xlabel(r'$\xi$')
plt.ylabel(r'$Y_i$')
plt.legend(frameon=False);
