II - Equilibrium calculations¶
Equilibrium explanations¶
Equilibrium calculations are usually performed to obtain the adiabatic flame temperature, the equilibrium composition, and the thermodynamic state of a specific mixture under given conditions. These are virtually performed in every simulation.
For example, Cantera will call its equilibrium solver to initialize the gas state before trying to obtain a solution to the equations for a free flame. As such, it is interesting to understand how Cantera proceeds.
There are 2 different types of solver currently implemented for equilibrium calculation in Cantera that deserves our attention :
- The element potential 'ChemEquil' solver
The class ChemEquil implements a chemical equilibrium solver for single-phase solutions. It is a "non-stoichiometric" solver in the terminology of Smith and Missen, meaning that every intermediate state is a valid chemical equilibrium state, but does not necessarily satisfy the element constraints.
Non-stoichiometric methods are faster when they converge, but stoichiometric ones tend to be more robust. - The 'VCS' chemical equilibrium solver
The other type of solver is designed to be used to set a mixture containing one or more phases to a state of chemical equilibrium. It uses a "stoichiometric" algorithm, in which each intermediate state satisfies the element constraints but is not a state of chemical equilibrium.
More specifically, it implements the VCS algorithm, described in Smith and Missen, "Chemical Reaction Equilibrium". It finds a set of component species and a complete set of formation reactions for the non-components in terms of the components.
As expected, the ChemEquil solver is the fastest of the Cantera equilibrium solvers for many single-
phase equilibrium problems (particularly if there are only a few elements but very many species), but
can be less stable.
Problem situations include low temperatures where only a few species have non-zero
mole fractions, precisely stoichiometric compositions (we will see an example shortly). In general, if
speed is important, this solver should always be tried first before falling back to another one in case of
failure.
The function equilibrate¶
The default setting in Cantera, when launching an equilibrium calculation without specifying
the solver, is to try the 'vcs' before falling back to another vcs solver labelled 'gibbs' :
gas.equilibrate('TP')
The equilibrate function can be applied on a single phase or on a mixture. Here, we recall its definition:
equilibrate(self, XY, solver, double rtol, int maxsteps, int maxiter, int loglevel)
Parameters:
XY - A two-letter string, which must be one of the set: ['TP','TV','HP','SP','SV','UV'].
solver - Specifies the equilibrium solver to use. May be one of the following :
element_potential = A fast solver using the element potential method.
gibbs = A slower but more robust Gibbs minimization solver.
vcs = The VCS non-ideal equilibrium solver.
auto = The element potential solver will be tried first, then if it fails the gibbs solver will be
tried.
rtol - The relative error tolerance.
maxsteps - Maximum number of steps in composition to take to find a converged solution.
maxiter - This specifies the number of outer iterations on T or P when some property pair other than TP is
specified (only for 'gibbs').
loglevel - Determines the amount of output displayed during the solution process. 0 indicates no output, while larger numbers produce successively more verbose information.
Useful python imports¶
import cantera as ct
import numpy as np
import csv
from matplotlib import *
import matplotlib.pyplot as plt
import sys
1. Simple homogeneous equilibrium of a gas¶
- the mechanism used is the gri30
- this is an air/CH4 mix
- the mix is at stoichiometry
- the temperature is 300K and the pressure 1 bar.
gas = ct.Solution('gri30.yaml') # create an object representing the gas phase
# set the composition of the gas
gas.TP = 300,100000
gas.set_equivalence_ratio(1.0,'CH4:1','O2:1, N2:3.79') # set initial state
gas.equilibrate('TP') # equilibrate using Temperature (T) and Pressure (P)
print(gas())
As you see, the gas has been equilibrated since it now shows only quantities for the product of the reaction (H2O and CO2). You can try to set yourself out of the perfect mixing (for example set CH4 to 0.4 and to 0.6) and see the impact on the species in the mix at the end.
2. Failure of a solver¶
Cantera has 3 different equilibrium solvers, 2 of them are worth mentionning:
- The 'ChemEquil' solver uses the element potential method for homogeneous equilibrium in gas mixtures. It is fast, but sometimes doesn't converge.
- The 'VCS' solver uses the VCS algorithm (Gibbs minimization), which is slower but more robust. It can also handle multiple phases.
Here we'll solve a problem for which the ChemEquil solver fails, but the VCS solver has no problem.
Properties of the gas¶
pressure = 1.0e5 # pressure
temperature = 400.0 # unburned gas temperature
comp = 'CH4:0.5, O2:1, N2:3.76' # premixed gas composition
gas = ct.Solution('gri30.yaml')
gas.TPX = temperature, pressure, comp
Initial state of the gas¶
print("******************************************************** ")
print(" Initial state :")
print("******************************************************** ")
print("P = ", "%10.4e " % (gas.P) + " Pa")
print("T = ", "%10.4e " % (gas.T) + " K")
print("V = ", "%10.4e " % (gas.volume_mass) + " m3/kg")
print("U = ", "%10.4e " % (gas.int_energy_mass) + " J/kg")
print("H = ", "%10.4e " % (gas.enthalpy_mass) + " J/kg")
print("S = ", "%10.4e " % (gas.entropy_mass) + " J/kg/K")
print("")
print("")
Comparing chemical potentials and element potentials¶
Here, the chemical potentials (noted mu_xx) are compared to the corresponding calculated values with the element potentials (noted lambda_xx). For instance, mu_H2 = lambda_H x 2. This is a good way to check wether the solver has managed to compute the results correctly. The chemical potentials are the one of the vapor phase.
chemeq = np.zeros(gas.n_species)
chemeq = gas.chemical_potentials
mu_H2 = chemeq[gas.species_index("H2")]
mu_OH = chemeq[gas.species_index("OH")]
mu_H2O = chemeq[gas.species_index("H2O")]
mu_O2 = chemeq[gas.species_index("O2")]
lambda_H = chemeq[gas.species_index("H")]
lambda_O = chemeq[gas.species_index("O")]
print()
print("Comparison between Chem potentials and element potentials:")
print()
s_mu_H2 = "%11.4e" % mu_H2
s_lam_mu_H2 = "%11.4e" % (2.0 * lambda_H)
print("mu_H2 : ", s_mu_H2, ", ", s_lam_mu_H2)
s_mu_O2 = "%11.4e" % mu_O2
s_lam_mu_O2 = "%11.4e" % (2.0 * lambda_O)
print("mu_O2 : ", s_mu_O2, ", ", s_lam_mu_O2)
s_mu_OH = "%11.4e" % mu_OH
s_lam_mu_OH = "%11.4e" % (lambda_H + lambda_O)
print("mu_OH : ", s_mu_OH, ", ", s_lam_mu_OH)
s_mu_H2O = "%11.4e" % mu_H2O
s_lam_mu_H2O = "%11.4e" % (2.0 * lambda_H + lambda_O)
print("mu_H2O : ", s_mu_H2O, ", ", s_lam_mu_H2O)
Program equilibrate¶
By default, the method equilibrate can be used without giving the solver. It will try first the element_potential solver than the vcs solver. Just below, the default behavior of Cantera is illustrated. Don't use that in your everyday calculation, just call gas.equilibrate('TP'), Cantera will handle the rest.
try:
# print("0")
gas.equilibrate('TP', solver='element_potential') # use the ChemEquil solver
except:
print("")
print("ChemEquil solver failed! Trying the vcs solver...")
gas.equilibrate('TP', solver='vcs', max_steps=1500)
Compare the results with the initial values¶
print("")
print("******************************************************** ")
print(" Final state :")
print(" Tadiabatique = " + str(gas.T) + " K")
print("******************************************************** ")
print("P = ", "%10.4e " % (gas.P) + " Pa")
print("T = ", "%10.4e " % (gas.T) + " K")
print("V = ", "%10.4e " % (gas.volume_mass) + " m3/kg")
print("U = ", "%10.4e " % (gas.int_energy_mass) + " J/kg")
print("H = ", "%10.4e " % (gas.enthalpy_mass) + " J/kg")
print("S = ", "%10.4e " % (gas.entropy_mass) + " J/kg/K")
print("")
print("")
Comparing chemical and element equilibrium for the equilibrate mixture¶
chemeq = gas.chemical_potentials
mu_H2 = chemeq[gas.species_index("H2")]
mu_OH = chemeq[gas.species_index("OH")]
mu_H2O = chemeq[gas.species_index("H2O")]
mu_O2 = chemeq[gas.species_index("O2")]
lambda_H = chemeq[gas.species_index("H")]
lambda_O = chemeq[gas.species_index("O")]
print()
print("Comparison between Chem potentials and element potentials:")
print()
s_mu_H2 = "%11.4e" % mu_H2
s_lam_mu_H2 = "%11.4e" % (2.0 * lambda_H)
print("mu_H2 : ", s_mu_H2, ", ", s_lam_mu_H2)
s_mu_O2 = "%11.4e" % mu_O2
s_lam_mu_O2 = "%11.4e" % (2.0 * lambda_O)
print("mu_O2 : ", s_mu_O2, ", ", s_lam_mu_O2)
s_mu_OH = "%11.4e" % mu_OH
s_lam_mu_OH = "%11.4e" % (lambda_H + lambda_O)
print("mu_OH : ", s_mu_OH, ", ", s_lam_mu_OH)
s_mu_H2O = "%11.4e" % mu_H2O
s_lam_mu_H2O = "%11.4e" % (2.0 * lambda_H + lambda_O)
print("mu_H2O : ", s_mu_H2O, ", ", s_lam_mu_H2O)
Saving the results¶
csv_file = '2-Output/all_mole_fractions.csv'
with open(csv_file, 'w') as outfile:
writer = csv.writer(outfile)
writer.writerow(['phi', 'T (K)'] + gas.species_names)
writer.writerow(['1', gas.T] + list(gas.X))
print(('Output written to {0}'.format(csv_file)))
3. Perform adiabatic flame calculations as a function of equivalence ratio¶
Set gas composition and interesting parameters¶
We will now perform several constant pressure equilibrium calculations of a methane/air mixture, at
an initial temperature of 300K and under atmospheric pressure, in order to obtain the adiabatic flame temperature and burnt gas state for several equivalence ratios (from 0.3 to 10.0).
The goal of this exercise is to find a way to loop through several gaseous composition, in order to perform several computations in a single script and to plot the result.
Time to code a real cantera script !¶
Then it is time for you to play ! Open another terminal folder and create your own Cantera python script. The goal is to compare the adiabatic flame temperature between gri 2.11 and gri3.0 (go search for them on https://chemistry.cerfacs.fr), for the following operating conditions :
- a temperature of 400K
- a pressure of 2 bars
- a range of equivalence ratio between 0.3 to 3.5 You should obtain the following results:
import numpy as np
import matplotlib.pyplot as plt
nb_pts = 100
temp = 400
pres = 200000
phi = np.linspace(0.3,10,nb_pts)
Tad_gas1 = np.zeros(nb_pts)
Tad_gas2 = np.zeros(nb_pts)
gas1 = ct.Solution('Mechanisms/GRI-Mech2.11.yaml')
gas2 = ct.Solution('gri30.yaml')
for i in range(0,len(phi)):
gas1.TP = temp,pres
gas2.TP = temp,pres
gas1.set_equivalence_ratio(phi[i],'CH4:1',{'O2':1,'N2':3.76})
gas2.set_equivalence_ratio(phi[i],'CH4:1',{'O2':1,'N2':3.76})
gas1.equilibrate('HP')
gas2.equilibrate('HP')
Tad_gas1[i] = gas1.T
Tad_gas2[i] = gas2.T
plt.figure
plt.plot(phi,Tad_gas1,label='GRI-Mech2.11')
plt.plot(phi,Tad_gas2,label='gri30')
plt.grid()
plt.legend()
# L'image présentée n'est pas pour de méthane mais de l'éthylène (C2H4)
4. Compute the low heating value and high heating value of a fuel¶
In combustion, the low heating value (LHV or PCI in French) corresponds to the energy released when the fuel and $O_2$ are converted into $CO_2$ and $H_2O$. It is directly linked to the formation enthalpy of the species involved and given by :
$$ LHV = \sum_i^c n_i \bar{h_i}(T_ref) $$Here, we want to compute the LHV for methane. The complete reaction can be written as:
$$ CH_4 + 2 O_2 \to CO_2 + 2 H_2O $$To compute the LHV, we will use the thermodynamics data from the gri30 mechanism. The mixture will be fixed at the reference state (298K, 1 atm).
First, we need to compute the enthalpy of the reactants :
import cantera as ct
gas = ct.Solution('gri30.yaml')
# Set the reactants state
gas.TP = 298, 101325
gas.set_equivalence_ratio(1.0, 'CH4:1.0', 'O2:1.0')
gas()
# Get the enthalpy and the mass fractions of the fuel
h_reactants = gas.enthalpy_mass
Y_CH4 = gas["CH4"].Y[0]
Now that we have the enthalpy of the reactants, we need to compute the enthalpy of the products :
# Set the products state
X_products = {'CO2': gas.elemental_mole_fraction('C'), 'H2O': 0.5*gas.elemental_mole_fraction('H')}
gas.TPX = 298, 101325, X_products
#or
#gas.equilibrate('TP') (sometimes you will get minor species that can alter a bit the LHV calculation)
gas()
# Get the products enthalpy
h_products = gas.enthalpy_mass
Now, we can compute the LHV for methane :
LHV = -(h_products - h_reactants) /Y_CH4
print("LHV [MJ/kg] = ", LHV/1e6)
Conclusion¶
We have just generated the skeleton of a script to perform a series of common equilibrium calculations
to obtain the constant pressure equilibrium composition of a fuel/air mixture. Starting from there,
you could modify your initial conditions, plot the mole/mass fractions of other species, change the
solver or even try another fuel (methane, acetylene) without changing your mechanism.
Technically, adiabatic flame calculations could also be performed at constant volume: simply invoke
the good equilibrate option of your equilibrate function, 'UV' (see 3.1.2), in your script.
The computation of the low heating value for a fuel has also been investigated in this tutorial