Matlab Tutorial¶
Getting Started¶
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:
>> gas1 = GRI30
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 (gas1) that 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://www.me.berkeley.edu/gri_mech/ for more information about GRI-Mech 3.0.)
The 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:
>> setTemperature(gas1, 1200);
sets the temperature to 1200 K (Cantera always uses SI units). After this statement, calling 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.):
>> set(gas1, 'Temperature', 900.0, 'Pressure', 1.e5);
This statement sets both temperature and pressure at the same
time. Any number of property/value pairs can be specified in a
call to set()
. For example, the following sets the mole fractions too:
>> set(gas1, 'Temperature', 900.0, 'Pressure', 1.e5, 'MoleFractions',... 'CH4:1,O2:2,N2:7.52');
The set()
method also accepts abbreviated property names:
>> set(gas1,'T',900.0,'P',1.e5,'X','CH4:1,O2:2,N2:7.52')
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:
>> set(gas1, 'Enthalpy', 2*enthalpy_mass(gas1), 'Pressure', 2*oneatm);
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 'X' to 'Y' in the above statement.
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:
>> x = ones(53,1); % a column vector of 53 ones >> set(gas1, 'X', x)
To set the mass fractions to equal values:
>> set(gas1, 'Y', x)
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 from file diamond.cti:
>> 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 importInterface()
.
The following command clears all Matlab objects created:
>> clear all
and this clears all Cantera objects created:
>> cleanup
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:
>> gas1 = Solution('gri30.cti', 'gri30');
Function 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 ‘C:\Program Files\Cantera\data‘ depending on how you installed Cantera and the options you specified. On a Unix/linux/macOS machine, they are usually kept in the 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:
>> adddir('/usr/local/cantera/my_data_files');
XML files¶
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.
>> gas4 = Solution('gri30.xml','gri30');
Interfaces can be imported from XML files too:
>> diamonnd_surf2 = importInterface('diamond.xml','diamond_100',... gas2, diamond);
Let’s clear out all our Matlab and Cantera objects, before we move on:
>> clear all >> cleanup
To learn more about the cti files already available with Cantera and how to create new cti files, see Working With Input Files
Getting Help¶
Suppose you have created a Cantera object and want to know what methods are available for it, and get help on using the methods.
>> g = GRI30
The first thing you need to know is the MATLAB class object g belongs to. Type:
>> class(g)
This tells you that g belongs to a class called Solution. To find the methods for this class, type
>> methods Solution
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
>> methods Solution -full
Now a long list is printed, along with a specification of the class the method is inherited from. For example, setPressure is 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 help <method_name> to print help text for any method. For example,
>> help setTemperature >> help setMassFractions >> help rop_net
For help on how to construct objects of a given class, type help <classname>
>> help Solution
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.
Chemical Equilibrium¶
To set a gas mixture to a state of chemical equilibrium, use the
equilibrate()
method.
>> set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52') >> equilibrate(g,'TP')
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:
- UV fixed specific internal energy and specific volume
- SV fixed specific entropy and specific volume
- SP fixed 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.)
Method 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:
Stoichiometric coefficients¶
>> 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¶
Methods rop_f()
, rop_r()
, and rop_net()
return column vectors
containing the forward, reverse, and net (forward - reverse) rates of progress, respectively, for
all reactions.
>> qf = rop_f(g); >> qr = rop_r(g); >> qn = rop_net(g); >> rop = [qf, qr, qn]
This plots the rates of progress
>> figure(1); >> bar(rop); >> legend('forward','reverse','net');
Species production rates¶
Methods creationRates()
, destructionRates()
, and netProdRates()
return
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:
>> figure(2); >> bar(rates); >> legend('creation','destruction','net');
For comparison, the production rates may also be computed directly from the rates of progress and stoichiometric coefficients.
>> cdot2 = nu_p*qf + nu_r*qr; >> creation = [cdot, cdot2, cdot - cdot2] >> ddot2 = nu_r*qf + nu_p*qr; >> destruction = [ddot, ddot2, ddot - ddot2] >> wdot2 = nu_net * qn; >> net = [wdot, wdot2, wdot - wdot2]
Reaction equations¶
>> e8 = reactionEqn(g,8) % equation for reaction 8 >> e1_10 = reactionEqn(g,1:10) % equation for rxns 1 - 10 >> eqs = reactionEqn(g) % all equations
Equilibrium constants¶
The equilibrium constants are computed in concentration units, with concentrations in kmol/m^3.
>> kc = equil_Kc(g); >> for i = 1:nReactions(g) >> disp(sprintf('%50s %13.5g', eqs{i}, kc(i))) >> end
Multipliers¶
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.
>> for i = 1:nReactions(g) >> setMultiplier(g, i, 2*i); >> m = multiplier(g, i); >> end
Let’s clear out the Matlab and Cantera objects, before moving on:
>> clear all >> cleanup
Transport Properties¶
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:
>> g1 = GRI30('Multi')
To use the mixture-averaged model:
>> g2 = GRI30('Mix')
Both models use a mixture-averaged formulation for the viscosity.
>> visc = [viscosity(g1), viscosity(g2)]
The thermal conductivity differs, however.
>> lambda = [thermalConductivity(g1), thermalConductivity(g2)]
Binary diffusion coefficients
>> bdiff1 = binDiffCoeffs(g1) >> bdiff2 = binDiffCoeffs(g2)
Mixture-averaged diffusion coefficients. For convenience, the multicomponent model implements mixture-averaged diffusion coefficients too.
>> dmix2 = mixDiffCoeffs(g1) >> dmix1 = mixDiffCoeffs(g2)
Multicomponent diffusion coefficients. These are only implemented if the multicomponent model is used.
>> dmulti = multiDiffCoeffs(g1)
Thermal diffusion coefficients. These are only implemented with the multicomponent model. These will be very close to zero, since the composition is pure H2.
>> dt = thermalDiffCoeffs(g1)
Now change the composition and re-evaluate
>> set(g1,'X',ones(nSpecies(g1),1)); >> dt = thermalDiffCoeffs(g1)
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:
>> clear all >> cleanup
Thermodynamic Properties¶
A variety of thermodynamic property methods are provided.
>> gas = air >> set(gas,'T',800,'P',oneatm)
Temperature, pressure, density:
>> T = temperature(gas) >> P = pressure(gas) >> rho = density(gas) >> n = molarDensity(gas)
Species non-dimensional properties:
>> hrt = enthalpies_RT(gas) % vector of h_k/RT
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:
>> clear all >> cleanup
Congratulations – Next Steps¶
Congratulations—you have finished the Cantera Matlab tutorial! You should now be ready to begin using Cantera on your own. Please see the Next Steps section on the Getting Started page, for assistance with intermediate and advanced Cantera functionality. Good luck!