Cantera Lecture

As a TA for CHEN 6153, I gave a lecture on the installation and use of Cantera. I made a video screencast of the lecture, which is available here: http://files.charlesmartinreid.com/CanteraScreencast.mov

You can download the example Matlab scripts here as well:

• Constant volume batch reactor (Matlab) http://files.charlesmartinreid.com/CanteraLecture_Batch.m
• Constant pressure batch reactor (Matlab) http://files.charlesmartinreid.com/CanteraLecture_BatchConstP.m
• Continuously stirred combustor (Matlab) http://files.charlesmartinreid.com/CanteraLecture_CSTR.m
• Adiabatic flame temperature (Python) http://files.charlesmartinreid.com/CanteraLecture_Adiabatic.py

Cantera Overview

Cantera is a code written in C++ that can be used to perform chemical kinetic and thermodynamic calculations.

The way it does this is by providing the user with a set of building blocks, and allowing the user to piece these building blocks together to form complex kinetic/thermodynamic systems or networks

Users can interface with Cantera through C++, or through the other interfaces:

• Matlab toolbox
• Python
• Fortran

Visit the Cantera (software) page for a Cantera installation guide for either Mac, Linux, or Windows.

Cantera Objects

Reservoir - fixed state

• the state of the reservoir never changes
• used to impose fixed upstream/inlet conditions
• chemistry is frozen - so no reactions will take place

Reactor - canonical chemical reactor

• contains some mixture of species
• has a set of characteristics
• e.g. T, V, P, composition, enthalpy, density, MW, etc...

Fluids - gases or liquids

Chemical Mechanisms - external files "fed" to Cantera

• chemical mechanisms specify two things:
  • species
  • reactions
• all thermodynamic data for species is specified (e.g. heat capacity as a polynomial function - "NASA coefficients")
• all thermodynamic data for chemical reactions (pre-exponential factors, activation energies, etc.)

Walls - separate reactors

• can have heat transfer via the usual mechanisms - convection, conduction, radiation
• can have a set heating rate (a constant, or a function of time)
• can move - due to a pressure difference, a displacement rate expression, etc.)

Flow controllers - control flow... 3 types

Mass flow controller - maintains a fixed mass flowrate

Valves - the mass flowrate is a function of $ \Delta P $

Pressure flow controller - maintains a fixed $ P $ in the reactor (think of this as a pressure relief valve)

Matlab Examples

Batch Reactor Example

How might we combine the above Cantera objects to make a batch reactor?

e.g. where would the reactor go? Where would the flow controllers go?

Any reservoirs needed? No

Any flow controllers needed? No

What will we need? A reactor

What are we going to fill the reactor with? We also need a fluid phase, and a mechanism to go along with it (if we want reactions), or thermodynamic data to go along with it (if we want thermodynamic state changes)

We want to know what happens when the reactor is changing from state 1 to state 2

For a constant volume reactor:

What is going into our reactor?

% Import a phase to fill the reactor with
% (The H2 combustion mechanism is contained in the file h2o2.cti, which comes with Cantera)
gas = importPhase('h2o2.cti');

What is the thermodynamic state of the material going into the reactor?

% Set the gas temperature at 1500 K
% Set the pressure at 1 atm
% Set the composition to be 2 parts hydrogen, 1 part oxygen, and 20 parts dilute, inert Argon
set(gas, 'T',1500, 'P',oneatm, 'X','H2:2, O2:1, AR:20')

How to create a batch reactor? Create a Cantera Reactor object

% We have already set the thermodynamic state of the "gas" object
% But if we wanted to, we could set it when we create the Reactor
batch = Reactor(gas);

We didn't specify the volume of the Reactor object - so how do we find it?

How do we get information about the Reactor object?

What information can we get about the Reactor object?

% We can use the volume() function to find the reactor's volume
volume(batch)

% We can find the Matlab object type like this:
class(batch)
% This returns "Reactor" type

% Show the functions available for Reactor objects
methods Reactor

The method "setInitialVolume()" looks good, so we can use it to set the initial volume of the reactor:

setInitialVolume( batch, 10.0 );

Another interesting method is "massFractions()":

% Display the mass fractions in the reactor
massFractions(batch)

This returns a 9-element vector, but it isn't clear what mass fractions are for which species.

How can we figure out what species are in the reactor?

Is this a reactor property?

Or is it a gas/fluid property?

So use the same methodology:

% Find the class type of the gas phase object
class(gas)
% This returns type "Solution"

% This shows the methods available to Solution types
methods Solution 

% This list looks a little short - but we can show inherited methods like this:
methods Solution -full

The "speciesNames()" function looks like what we want:

speciesNames(gas)

This returns a 9-element vector of strings - names of species in the mechanism file 'h2o2.cti'.

This plot of Cantera batch reactor results is a stub. A team of expertly-trained librarian ninjas has been dispatched to research this topic, finish this plot of Cantera batch reactor results, and assassinate all eyewitnesses who contradict their account. They should be done pretty soon. Template:Stub Category:Stubs

Constant Pressure Reactor Example

Now for a more interesting (relevant) combustion problem:

How might we use reactors, walls, reservoirs, etc. to create a constant-pressure batch reactor?

We can create a piston-cylinder system - this would be a constant-pressure system, with combustion occurring inside the piston

The constant-pressure, constant-temperature, constant-property atmosphere is a reservoir

The piston lid is a wall

The reactants and products inside the piston are gases

So first, we'll create a gas phase to fill the piston with

% Create the H2-O2 gas phase
gas = importPhase('h2o2.cti');

% Set the thermodynamic state of the gas
set(gas, 'T', 1500, 'P', oneatm, 'X', 'H2:2,O2:1,AR:20');

Next, create the reservoir that represents the atmosphere, and fill it with air:

% Create an ideal gas mix of air
a = IdealGasMix('air.cti');
atmosphere = reservoir(a);

A batch reactor must be created to represent the piston reactor:

batch = Reactor(gas);

The next step is to create a wall representing the piston between the reservoir and the gas. Walls in Cantera can be thought of as pistons with a certain amount of friction, or weight.

The wall expansion parameter is defined as:

$ \frac{dV}{dt} = K_{\mbox{expansion}} A_{\mbox{wall}} ( P_1 - P_2 ) $

What happens when $ K=0 $?

What happens when $ K=\infty $?

The larger the $ K $ value, the faster the response of the wall (piston).

% Create a wall between the reactor and the atmosphere
w = Wall;
install(w, batch, atmosphere);
setArea(w,1.0);

% Set the expansion coefficient K
setExpansionRateCoeff(w, 1.0e10);

Remember, we can also have walls with heat transfer, so we could have heat flowing in, or flowing out, through the wall

By default, the wall is adiabatic

Next, a reactor network can be created, in order to advance the reactor in time.

network = ReactorNet( {batch} );

CSTR Combustor

This example is for a continuously-fired combustor.

The reactor is a combustor with a fuel and oxidizer reservoir, as well as a reservoir of igniter to initiate combustion.

There is also a reservoir for the exhaust.

So we will need some reservoirs, a reactor...

What else? Flow controllers

We want to run the CSTR with various fuel/oxidizer ratios, and see how it affects the concentrations, temperatures, etc.

We might want to vary residence time, too - see the extent of combustion as a function of residence time

Create a reservoir of fuel, a reservoir of air, and a reservoir of hydrogen radicals (radicals are sure to initiate combustion reactions).

% Create a fuel from the GRI combustion mechanism
fuel = GRI30;
set(fuel, 'T', 300, 'P', oneatm, 'X', 'C3H8:1.0');
fuel_mw = meanMolecularWeight(fuel);

% Create the fuel reservoir
fuel_in = Reservoir(fuel);

% Create the oxidizer
oxidizer = Air();
set(oxidizer, 'T', 300, 'P', oneatm);
oxidizer_mw = meanMolecularWeight(oxidizer);

% Create the oxidizer reservoir
oxidizer_in = Reservoir(oxidizer);

% Create the igniter reservoir
igniter = GRI30;
set(igniter, 'T', 300, 'P', oneatm, 'X', 'H:1.0');
igniter_in = Reservoir(igniter);

The combustor must be initially filled with something - an inert gas is a good idea

inert = GRI30;
set(inert, 'T', 300, 'P', oneatm, 'X', 'N2:1.0');
combustor = Reactor(inert);

The next step is to set the equivalence ratio, and set the mass flowrates accordingly.

ER = 0.7;

% Fuel mass flowrate
fuel_mdot = fuel_mw*ER;

% Oxidizer mass flowrate
oxidizer_mdot = (1/ER)*(1/.21)*oxidizer_mw;

% Igniter mass flowrate
igniter_mdot = 0.05;

The mass flowrates could also be a Gaussian function in time, or a polynomial, e.g.

igniter_mdot = Func('gaussian', 0, [1.0, 0.0, 0.1] );
igniter_mdot = Func('polynomial', 2.0, [0.0, 0.0, 0.1] );

You can type "help Func" for more information.

(However, this does not seem to be working...)

The next step is to create the mass flow controllers entering the reactor, and the pressure exhaust valve, and finally to set up the reactor network.

% Create the mass flowrate controllers feeding into the combustor
m1 = MassFlowController(oxidizer_in, combustor);
setMassFlowRate(m1, oxidizer_mdot);

m2 = MassFlowController(fuel_in, combustor);
setMassFlowRate(m2, fuel_mdot);

m3 = MassFlowController(igniter_in, combustor);
setMassFlowrate(m3, igniter_mdot);

% And create an exhaust valve regulating the pressure in the combustor
exhaust_valve = Valve(combustor,exhaust);
setValveCoeff(exhaust_valve,1.0);

% Now the reactor network can be created
network = ReactorNet( {combustor} );

Python Examples

Before You Begin

Before you start to use Cantera from Python, I highly recommend installing Py4Sci, a collection of 4 very handy packages that greatly extend Python and give a more integrated development environment. Essentially, it makes Python behave more like Matlab by giving it numerical capabilities, a nicer user shell, and some great plotting functions/tools.

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: