First, you’ll need to install Cantera on your computer. We have instructions for many platforms in our Installation section.
When using Cantera, the first thing you usually need is an object representing some phase of matter. Here, we’ll create a gas mixture. Start MATLAB, and at the prompt type:
If you have successfully installed the Cantera toolbox, you should see something like this:
gri30: temperature 300 K pressure 101325 Pa density 0.0818891 kg/m^3 mean mol. weight 2.01588 amu 1 kg 1 kmol ----------- ------------ enthalpy 26470.1 5.336e+04 J internal energy -1.21087e+06 -2.441e+06 J entropy 64913.9 1.309e+05 J/K Gibbs function -1.94477e+07 -3.92e+07 J heat capacity c_p 14311.8 2.885e+04 J/K heat capacity c_v 10187.3 2.054e+04 J/K X Y Chem. Pot. / RT ------------- ------------ ------------ H2 1 1 -15.7173 [ +52 minor] 0 0
What you have just done is to create an object (
implements GRI-Mech 3.0, the 53-species, 325-reaction natural gas
combustion mechanism developed by Gregory P. Smith, David M. Golden,
Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer, Mikhail
Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, William
C. Gardiner, Jr., Vitali V. Lissianski, and Zhiwei Qin. (See
http://combustion.berkeley.edu/gri-mech/ for more information about
gas1 object has properties you would expect for a gas mixture—it has a
temperature, a pressure, species mole and mass fractions, etc. As we’ll soon
see, it has many more properties.
The summary of the state of
gas1 printed above shows that new objects
created from the
gri30.cti input file start out with a temperature of 300 K,
a pressure of 1 atm, and have a composition that consists of only one species,
in this case hydrogen. There is nothing special about H2, it just happens to
be the first species listed in the input file defining GRI-Mech 3.0. In
general, whichever species is listed first will initially have a mole fraction
of 1.0, and all of the others will be zero.
Setting the State¶
The state of the object can easily be changed. For example:
sets the temperature to 1200 K (Cantera always uses SI units). After this
gas1() results in:
gri30: temperature 1200 K pressure 405300 Pa density 0.0818891 kg/m^3 mean mol. weight 2.01588 amu 1 kg 1 kmol ----------- ------------ enthalpy 1.32956e+07 2.68e+07 J internal energy 8.34619e+06 1.682e+07 J entropy 85227.6 1.718e+05 J/K Gibbs function -8.89775e+07 -1.794e+08 J heat capacity c_p 15377.9 3.1e+04 J/K heat capacity c_v 11253.4 2.269e+04 J/K X Y Chem. Pot. / RT ------------- ------------ ------------ H2 1 1 -17.9775 [ +52 minor] 0 0
Notice that the temperature has been changed as requested, but the pressure has changed too. The density and composition have not.
When setting properties individually, some convention needs to be adopted to specify which other properties are held constant. This is because thermodynamics requires that two properties (not one) in addition to composition information be specified to fix the intensive state of a substance (or mixture).
Cantera adopts the following convention: only one of the set (temperature, density, mass fractions) is altered by setting any single property. In particular:
Setting the temperature is done holding density and composition fixed. (The pressure changes.)
Setting the pressure is done holding temperature and composition fixed. (The density changes.)
Setting the composition is done holding temperature and density fixed. (The pressure changes).
If you want to set multiple properties at once, use the
set() method. (Note: a
method is just the term for a function that acts on an object. In MATLAB,
methods take the object as the first argument.):
This statement sets both temperature and pressure at the same
time. Any number of property/value pairs can be specified in a
set(). For example, the following sets the mole fractions too:
set() method also accepts abbreviated property names:
Either version results in:
gri30: temperature 900 K pressure 100000 Pa density 0.369279 kg/m^3 mean mol. weight 27.6332 amu 1 kg 1 kmol ----------- ------------ enthalpy 455660 1.259e+07 J internal energy 184862 5.108e+06 J entropy 8529.31 2.357e+05 J/K Gibbs function -7.22072e+06 -1.995e+08 J heat capacity c_p 1304.4 3.604e+04 J/K heat capacity c_v 1003.52 2.773e+04 J/K X Y Chem. Pot. / RT ------------- ------------ ------------ O2 0.190114 0.220149 -27.9596 CH4 0.095057 0.0551863 -37.0813 N2 0.714829 0.724665 -24.935 [ +50 minor] 0 0
Other properties may also be set using
set(), including some that
can’t be set individually. The following property pairs may be
set: (Enthalpy, Pressure), (IntEnergy, Volume), (Entropy,
Volume), (Entropy, Pressure). In each case, the values of the
extensive properties must be entered per unit mass.
Setting the enthalpy and pressure:
The composition above was specified using a string. The format is a
comma-separated list of <species name>:<relative mole numbers>
pairs. The mole numbers will be normalized to produce the mole
fractions, and therefore they are ‘relative’ mole numbers. Mass
fractions can be set in this way too by changing
'Y' in the
The composition can also be set using an array, which can be either a column vector or a row vector but must have the same size as the number of species. For example, to set all 53 mole fractions to the same value, do this:
To set the mass fractions to equal values:
Importing multiple phases or interfaces¶
A Cantera input file may contain more than one phase specification,
or may contain specifications of interfaces (surfaces). Here we
import definitions of two bulk phases and the interface between them
>> gas2 = Solution('diamond.cti', 'gas'); % a gas >> diamond = Solution('diamond.cti','diamond'); % bulk diamond >> diamonnd_surf = importInterface('diamond.cti','diamond_100',... gas2, diamond);
Note that the bulk (i.e., 3D) phases that participate in the surface
reactions must also be passed as arguments to
The following command clears all Matlab objects created:
and this clears all Cantera objects created:
Working with input files¶
Previously, we used the function
GRI30() to create an object that models an ideal
gas mixture with the species and reactions of GRI-Mech 3.0. Another way to do
this is shown here:
Solution() constructs an object representing a phase of
matter by reading in attributes of the phase from a file, which in
this case is
gri30.cti. This file contains several phase
specifications; the one we want here is
gri30, which is specified
by the second argument. This file contains a complete specification
of the GRI-Mech 3.0 reaction mechanism, including element data
(name, atomic weight), species data (name, elemental composition,
coefficients to compute thermodynamic and transport properties), and
reaction data (stoichiometry, rate coefficient parameters).
CTI files distributed with Cantera¶
Several reaction mechanism files in this format are included in the
Cantera distribution, including ones that model high-temperature
air, a hydrogen/oxygen reaction mechanism, and a few surface
reaction mechanisms. Under Windows, these files may be located in
depending on how you installed Cantera and the options you
specified. On a Unix/linux/macOS machine, they are usually kept
data subdirectory within the Cantera installation directory.
If for some reason Cantera has difficulty finding where these files
are on your system, set the environment variable
CANTERA_DATA to the
directory where they are located. Alternatively, you can call function
adddir() to add a directory to the Cantera search path:
Note that Cantera has two input formats for data, the human-readable and -writable CTI format, and the lower-level XML format. All of the CTI files distributed with Cantera are also available as XML files; using the XML files may be somewhat faster and does not invoke the Python interpreter to read the CTI file. More information on why Cantera uses two file formats is available in the input files tutorial.
Interfaces can be imported from XML files too:
Let’s clear out all our Matlab and Cantera objects, before we move on:
To learn more about the cti files already available with Cantera and how to create new cti files, see Working With Input Files
Suppose you have created a Cantera object and want to know what methods are available for it, and get help on using the methods.
The first thing you need to know is the MATLAB class object
belongs to. Type:
This tells you that
g belongs to a class called
Solution. To find
the methods for this class, type
This command returns only a few method names. These are the ones
directly defined in this class. But
Solution inherits many other
methods from base classes. To see all of its methods, type
Now a long list is printed, along with a specification of the class
the method is inherited from. For example,
inherited from a class
ThermoPhase. Don’t be concerned at this
point about what these base classes are—we’ll come back to them later.
Now that you see what methods are available, you can type
<method_name> to print help text for any method. For example,
For help on how to construct objects of a given class, type
Now that you know how to get help when you need it, you can explore using the Cantera Toolbox on your own. But there are a few more useful things to know, which are described in the next few sections.
To set a gas mixture to a state of chemical equilibrium, use the
The statement above sets the state of object
g to the state of
chemical equilibrium holding temperature and pressure
fixed. Alternatively, the specific enthalpy and pressure can be held fixed:
>> disp('fixed H and P:'); >> set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2.0,N2:7.52'); >> equilibrate(g,'HP')
Other options are:
UVfixed specific internal energy and specific volume
SVfixed specific entropy and specific volume
SPfixed specific entropy and pressure
>> disp('fixed U and V:'); >> set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52'); >> equilibrate(g,'UV') >> disp('fixed S and V:'); >> set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52'); >> equilibrate(g,'SV') >> disp('fixed S and P:'); >> set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52'); >> equilibrate(g,'SP')
How can you tell if
equilibrate() has correctly found the
chemical equilibrium state? One way is verify that the net rates of
progress of all reversible reactions are zero.
Here is the code to do this:
>> set(g,'T',2000.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52'); >> equilibrate(g,'TP') >> rf = rop_f(g); >> rr = rop_r(g); >> format short e; >> for i = 1:nReactions(g) >> if isReversible(g,i) >> disp([i, rf(i), rr(i), (rf(i) - rr(i))/rf(i)]); >> end >> end
You might be wondering how
equilibrate() works. (Then again, you might not.)
equilibrate() invokes Cantera’s chemical equilibrium solver, which uses
an element potential method. The element potential method is one of a class of equivalent
‘nonstoichiometric’ methods that all have the characteristic that the problem reduces to solving a
set of \(M\) nonlinear algebraic equations, where \(M\) is the number of elements (not
species). The so-called ‘stoichiometric’ methods, on the other hand, (including Gibbs minimization),
require solving K nonlinear equations, where \(K\) is the number of species (usually \(K >>
M\)). See Smith and Missen, “Chemical Reaction Equilibrium Analysis” for more information on the
various algorithms and their characteristics.
Cantera uses a damped Newton method to solve these equations, and does a few other things to generate a good starting guess and to produce a reasonably robust algorithm. If you want to know more about the details, look at the C++ code in ChemEquil.h.
Reaction information and rates¶
Methods are provided that compute many quantities of interest for kinetics. Some of these are:
>> set(g,'T',1500,'P',oneatm,'X',ones(nSpecies(g),1)); >> nu_r = stoich_r(g) % reactant stoichiometric coefficient mstix >> nu_p = stoich_p(g) % product stoichiometric coefficient mstix >> nu_net = stoich_net(g) % net (product - reactant) stoichiometric % coefficient mstix
For any of these, the
(k,i) matrix element is the stoichiometric
coefficient of species \(k\) in reaction \(i\). Since these coefficient
matrices are very sparse, they are implemented as MATLAB sparse matrices.
Reaction rates of progress¶
This plots the rates of progress
Species production rates¶
column vectors containing the creation, destruction, and net
production (creation - destruction) rates, respectively, for all species.
>> cdot = creationRates(g); >> ddot = destructionRates(g); >> wdot = netProdRates(g); >> rates = [cdot, ddot, wdot]
This plots the production rates:
For comparison, the production rates may also be computed directly from the rates of progress and stoichiometric coefficients.
The equilibrium constants are computed in concentration units, with concentrations in kmol/m^3.
For each reaction, a multiplier may be specified that is applied to the forward rate coefficient. By default, the multiplier is 1.0 for all reactions.
Let’s clear out the Matlab and Cantera objects, before moving on:
Methods are provided to compute transport properties. By default, calculation of transport properties is not enabled. If transport properties are required, the transport model must be specified when the gas mixture object is constructed.
Currently, two models are implemented. Both are based on kinetic theory expressions, and follow the approach described in Dixon-Lewis (1968) and Kee, Coltrin, and Glarborg (2003). The first is a full multicomponent formulation, and the second is a simplification that uses expressions derived for mixtures with a small number of species (1 to 3), using approximate mixture rules to average over composition.
To use the multicomponent model with GRI-Mech 3.0, call function GRI30 as follows:
To use the mixture-averaged model:
Both models use a mixture-averaged formulation for the viscosity.
The thermal conductivity differs, however.
Binary diffusion coefficients
Mixture-averaged diffusion coefficients. For convenience, the multicomponent model implements mixture-averaged diffusion coefficients too.
Multicomponent diffusion coefficients. These are only implemented if the multicomponent model is used.
Thermal diffusion coefficients. These are only implemented with the multicomponent model. These will be very close to zero, since the composition is pure H2.
Now change the composition and re-evaluate
Note that there are no singularities for pure gases. This is because a very small positive value is added to all mole fractions for the purpose of computing transport properties.
Let’s clear out the Matlab and Cantera objects, before moving on:
A variety of thermodynamic property methods are provided.
Temperature, pressure, density:
Species non-dimensional properties:
Mixture properties per mole:
>> hmole = enthalpy_mole(gas) >> umole = intEnergy_mole(gas) >> smole = entropy_mole(gas) >> gmole = gibbs_mole(gas)
Mixture properties per unit mass:
>> hmass = enthalpy_mass(gas) >> umass = intEnergy_mass(gas) >> smass = entropy_mass(gas) >> gmass = gibbs_mass(gas)
Le’ts do one final clearing of the workspace: