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 'element potential' 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 - Is currently deprecated.
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.cti') # create an object representing the gas phase
gas.TPX = 300, 100000,{'CH4':1, 'O2':2, 'N2':7.52} # 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.xml')
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 whether 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¶
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', maxsteps=1500)
# gas.equilibrate('TP', solver = 'gibbs') # the gibbs solver works also
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 an ethylene/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.
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 learn how to properly store the results in a csv file.
gas = ct.Solution('Mechanisms/mech.cti')
T = 300.0 # Temperature
P = 101325.0 # Pressure
phi_min = 0.3 # Minimal equivalence ratio
phi_max = 3.5 # Maximal equivalence ratio
npoints = 50 # Point in-between the two preceeding values
fuel_species = 'C2H4' # fuel species
air_N2_O2_molar_ratio = 3.76 # ratio representing the air
Loop of equilibrium¶
phi = np.zeros(npoints) # 1D array
tad = np.zeros(npoints) # 1D array
xeq = np.zeros((gas.n_species, npoints)) # 2D array
for i in range(npoints):
gas.TP = T, P
phi[i] = phi_min + (phi_max - phi_min) * i / (npoints - 1)
gas.set_equivalence_ratio(phi[i], {fuel_species: 1}, {'O2': 1, 'N2': air_N2_O2_molar_ratio})
gas.equilibrate('HP') # Equilibrate the mixture adiabatically at constant P with the solver vcs
xeq[:, i] = gas.X
tad[i] = gas.T
print("At phi = ","%10.4f"% (phi[i])+ " Tad = ","%10.4f"% (tad[i]))
Save files¶
csv_file = '2-Output/yourfile.csv'
with open(csv_file, 'w') as outfile:
writer = csv.writer(outfile)
writer.writerow(['Phi', 'T (K)'] + gas.species_names)
for i in range(npoints):
writer.writerow([phi[i], tad[i]] + list(xeq[:, i]))
print("Output written to", "%s" % csv_file)
Plot the results¶
rcParams['figure.figsize'] = (14, 10)
plt.plot(phi, tad, '-')
plt.title(r'Tad vs. $\Phi$ for $C_{2}H_{4}/Air$ flames')
plt.xlabel(r'$\Phi$', fontsize=20)
plt.ylabel("Adiabatic flame temperature (K)")
plt.grid()
plt.show()
plt.savefig('2-Output/plot_tad.png', bbox_inches='tight')
for i, cas in enumerate(gas.species_names):
if cas in ['O2','CO2','CO']:
plt.plot(phi,xeq[i,:], label = cas)
plt.xlabel('Equivalence ratio')
plt.ylabel('Mole fractions')
plt.legend(loc='best')
plt.show()
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.