Using Cantera

Main page: Cantera

This page covers some examples of Cantera usage for several of the languages for which Cantera provides an interface.

Matlab

The example Cantera Matlab scripts are available here:

• Constant volume batch reactor (Matlab): http://files.charlesmartinreid.com/CanteraLecture_Batch.m (file is gone)
• Constant pressure batch reactor (Matlab): http://files.charlesmartinreid.com/CanteraLecture_BatchConstP.m (file is gone)
• Continuously stirred combustor (Matlab): http://files.charlesmartinreid.com/CanteraLecture_CSTR.m (file is gone)

A walkthrough of these example scripts is given on the Cantera Lecture#Matlab Examples page.

Python

The example Python script that I will be walking through is an adiabatic flame temperature calculation example.

Before You Begin

You will need Python, as well as Numpy, Scipy, Jupyter notebook, and some other essential scientific computing libraries for Python. These are handy packages that will simplify the development and scripting of Cantera work.

The below example makes use of these packages.

Cantera Python Scripts

To make use of the above packages, and to make use of Cantera, I put the following in all Cantera Python scripts:

from Cantera import *
from pylab import *
import sys

Adiabatic Flame Temperature

Problem Description: given a fuel-oxidizer mixture at a given temperature and pressure, determine the adiabatic flame temeprature. Specifically, determine the adiabatic flame temperature as a function of equivalence ratio.

Setting Up

First you'll want to define the conditions and state of your fuel:

# Edit these parameters to change the initial temperature, the
# pressure, and the phases in the mixture
temp = 298.0
pres = 101325.0

# phases
gas  = importPhase('gri30.cti')
carbon = importPhase('graphite.cti')

# the phases that will be included in the calculation, and their
# initial moles
mix_phases = [ (gas, 1.0), (carbon, 0.0) ]

# gaseous fuel species
fuel_species = 'C3H8'

# air composition
air_N2_O2_ratio = 3.76

# create a mixture object out of the fuel and oxidizer
mix = Mixture(mix_phases)
nsp = mix.nSpecies()

Next, define the range of equivalence ratios over which to calculate the adiabatic flame temperature:

# equivalence ratio range
phi_min = 0.8
phi_max = 1.2
npoints = 10

It is useful to know the index of the fuel and oxidizer species, as well as nitrogen. These can be obtained by using the speciesIndex() method of the GRI Mech phase created above:

# find fuel, nitrogen, and oxygen indices
ifuel, io2, in2 = gas.speciesIndex([fuel_species, 'O2', 'N2'])

if ifuel < 0:
    raise "ERROR: fuel species "+fuel_species+" not present!"

if gas.nAtoms(fuel_species,'O') > 0 or  gas.nAtoms(fuel_species,'N') > 0:
    raise "ERROR: only hydrocarbon fuels are supported."

Finally, it will be important to have the stoichiometric coefficient for O2 for the complete combustion reaction (i.e. assuming reactants go completely to products, and products consist only of CO2 and H2O):

stoich_O2 = gas.nAtoms(fuel_species,'C') + 0.25*gas.nAtoms(fuel_species,'H')
#           ^^^^^^^^^^                      ^^^^^^^^^^^^^^^
#         accounts for formation of CO2     accounts for formation of H2O (4 H's per O2)

Oh, and to speed things up, you can pre-allocate some arrays to store values of $\Phi ,T_{\mbox{adiabatic}},x_{i}$ (where the concentrations are the equilibrium concentrations)

# create some arrays to hold the data
phi = zeros(npoints,'d')
tad = zeros(npoints,'d')
xeq = zeros([nsp,npoints],'d')

Calculation of Adiabatic Flame T

The next section takes place inside a loop over all the values of equivalence ratios $\Phi$ specified:

for i in range(npoints):
   phi[i] = phi_min + (phi_max - phi_min)*i/(npoints - 1)

Above, the mix_phases object was created, which was a mixture of 1 mole of the GRI Mech gas phase, and 0 moles of the solid carbon phase:

mix_phases = [ (gas, 1.0), (carbon, 0.0) ]

But this doesn't set the composition of the 1 mole of gas. This can be done by creating a vector of compositions x, then specifying the thermodynamic state of the gas phase by specifying the temperature, pressure, and composition:

   x = zeros(nsp,'d')
   x[ifuel] = phi[i]
   x[io2] = stoich_O2
   x[in2] = stoich_O2*air_N2_O2_ratio

   # set the gas state
   gas.set(T = temp, P = pres, X = x)

Next, a Cantera Mixture object can be created, and its temperature and pressure set:

   # create a mixture of 1 mole of gas, and 0 moles of solid carbon.
   mix = Mixture(mix_phases)
   mix.setTemperature(temp)
   mix.setPressure(pres)

The mixture can then be equilibrated at constant pressure and enthalpy

   mix.equilibrate('HP', maxsteps = 1000,
                   err = 1.0e-6, maxiter = 200, loglevel=0)

Once the mixture is equilibrated, the temperature (the adiabatic flame temperature) can be obtained using the temperature() method:

   tad[i] = mix.temperature();
   print 'At phi = %12.4g, Tad = %12.4g'  % (phi[i],tad[i])

and the composition of the equilibrated phase can be obtained using the speciesMoles() method:

   xeq[:,i] = mix.speciesMoles()

Outputting Data File

Data generated can be output to a .csv file using the writeCSV() method:

# create array with species names
neq = mix.speciesNames()

# write output CSV file for importing into Excel
csvfile = 'adiabatic.csv'
f = open(csvfile,'w')
writeCSV(f,['phi','T (K)']+mix.speciesNames())
for n in range(npoints):
    writeCSV(f,[phi[n], tad[n]]+list(xeq[:,n]))
f.close()
print 'output written to '+csvfile

Create Plots

There are three different plots that would be useful for the adiabatic flame temperature case:

• Adiabatic flame temperature vs. equivalence ratio $\Phi$
• Concentration of major species vs. equivalence ratio
• Concentration of minor species vs. equivalence ratio

Start by creating a list of species that are interesting:

species = list()
species.append('CO2')
species.append('CO')
species.append('H2O')
species.append('H2')
species.append('OH')
species.append('H')
species.append('O2')
species.append('O')
species.append('NO')
species.append('N2')
species.append('N')

and define a mole fraction limit to distinguish major and minor species:

xlim = 0.04

Create a figure for each of the plots:

fig=figure(1)

You can use Matplotlib to create the plot of temperature vs. equivalence ratio. The syntax is very similar to Matlab.

a=fig.add_subplot(131)
a.plot(phi,tad)
title(r'$T_{adiabatic}$ vs. $\Phi$')
xlabel(r'$\Phi$', fontsize=20)
ylabel("Adiabatic Flame Temperature [K]")
a.xaxis.set_major_locator(MaxNLocator(5)) # this controls the number of tick marks on the axis

The next plot is concentration of major species as a function of equivalence ratio:

b=fig.add_subplot(132)
print('\n\nSpecies | Conc. at Phi = 0.8 | Conc. at Phi = 1.2')
print('\nMajor species:')
for j in range(size(neq)):
    #for each species...
    if (neq[j] in species):
        #if it is a species of interest for our problem...
        if (False in (xeq[j] < xlim) ):
            # if j is a MAJOR species:
            print(neq[j],round(xeq[j,0],6),round(xeq[j,npoints-1],6))
            b.plot(phi,xeq[j],label=neq[j])
            #legend(loc='best')
            hold(True)

hold(False)
legend(loc='best')
title(r'Major species')
xlabel(r'$\Phi$', fontsize=20)
ylabel("Mole fraction")
b.xaxis.set_major_locator(MaxNLocator(5))

The last plot is concentration of minor species:

c=fig.add_subplot(133)
print('\nMinor species:')
for j in range(size(neq)):
    #for each species...
    if (neq[j] in species):
        #if it is a species of interest for our problem...
        if (False not in (xeq[j] < xlim) ):
            # if j is a MINOR species:
            # xeq[row,col]
            print(neq[j],round(xeq[j,0],6),round(xeq[j,npoints-1],6))
            c.plot(phi,xeq[j],label=neq[j])
            hold(True)

hold(False)
legend(loc='best')
title(r'Minor Species')
xlabel(r'$\Phi$', fontsize=20)
c.xaxis.set_major_locator(MaxNLocator(5))

Finally, a title can be added to the plot:

fig.text(0.5,0.95,r'Adiabatic Flame Temperature for $C_{3}H_{8}$ + Air', ...
 fontsize=22,horizontalalignment='center')

show()

This results in the following plot:

C++

(Examples forthcoming)

There are some C++ examples included with Cantera, located in /path/to/cantera/1.8.0/demos/cxx

I have also detailed the implementation of the Jacobian matrix in Cantera: Jacobian in Cantera