From charlesreid1

No edit summary
(Fix math: replace \mbox with \text for MediaWiki math renderer compatibility (via update-page on MediaWiki MCP Server))
 
(9 intermediate revisions by the same user not shown)
Line 1: Line 1:
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


As a TA for CHEN 6153, I gave a lecture on the installation and use of Cantera.


= Cantera Overview =
= Cantera Overview =
Line 79: Line 74:
What is going into our reactor?
What is going into our reactor?


<syntaxhighlight lang="matlab">
<pre>
% Import a phase to fill the reactor with
% Import a phase to fill the reactor with
% (The H2 combustion mechanism is contained in the file h2o2.cti, which comes with Cantera)
% (The H2 combustion mechanism is contained in the file h2o2.cti, which comes with Cantera)
gas = importPhase('h2o2.cti');
gas = importPhase('h2o2.cti');
</syntaxhighlight>
</pre>


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


<syntaxhighlight lang="matlab">
<pre>
% Set the gas temperature at 1500 K
% Set the gas temperature at 1500 K
% Set the pressure at 1 atm
% Set the pressure at 1 atm
% Set the composition to be 2 parts hydrogen, 1 part oxygen, and 20 parts dilute, inert Argon
% 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')
set(gas, 'T',1500, 'P',oneatm, 'X','H2:2, O2:1, AR:20')
</syntaxhighlight>
</pre>


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


<syntaxhighlight lang="matlab">
<pre>
% We have already set the thermodynamic state of the "gas" 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
% But if we wanted to, we could set it when we create the Reactor
batch = Reactor(gas);
batch = Reactor(gas);
</syntaxhighlight>
</pre>


We didn't specify the volume of the Reactor object - so how do we find it?
We didn't specify the volume of the Reactor object - so how do we find it?
Line 108: Line 103:
What information ''can'' we get about the Reactor object?
What information ''can'' we get about the Reactor object?


<syntaxhighlight lang="matlab">
<pre>
% We can use the volume() function to find the reactor's volume
% We can use the volume() function to find the reactor's volume
volume(batch)
volume(batch)
Line 118: Line 113:
% Show the functions available for Reactor objects
% Show the functions available for Reactor objects
methods Reactor
methods Reactor
</syntaxhighlight>
</pre>


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


<syntaxhighlight lang="matlab">
<pre>
setInitialVolume( batch, 10.0 );
setInitialVolume( batch, 10.0 );
</syntaxhighlight>
</pre>


Another interesting method is "massFractions()":
Another interesting method is "massFractions()":


<syntaxhighlight lang="matlab">
<pre>
% Display the mass fractions in the reactor
% Display the mass fractions in the reactor
massFractions(batch)
massFractions(batch)
</syntaxhighlight>
</pre>


This returns a 9-element vector, but it isn't clear what mass fractions are for which species.
This returns a 9-element vector, but it isn't clear what mass fractions are for which species.
Line 143: Line 138:
So use the same methodology:
So use the same methodology:


<syntaxhighlight lang="matlab">
<pre>
% Find the class type of the gas phase object
% Find the class type of the gas phase object
class(gas)
class(gas)
Line 153: Line 148:
% This list looks a little short - but we can show inherited methods like this:
% This list looks a little short - but we can show inherited methods like this:
methods Solution -full
methods Solution -full
</syntaxhighlight>
</pre>


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


<syntaxhighlight lang="matlab">
<pre>
speciesNames(gas)
speciesNames(gas)
</syntaxhighlight>
</pre>


This returns a 9-element vector of strings - names of species in the mechanism file 'h2o2.cti'.
This returns a 9-element vector of strings - names of species in the mechanism file 'h2o2.cti'.
{{Stub|plot of Cantera batch reactor results}}


== Constant Pressure Reactor Example ==
== Constant Pressure Reactor Example ==
Line 185: Line 176:
So first, we'll create a gas phase to fill the piston with
So first, we'll create a gas phase to fill the piston with


<source lang="matlab">
<pre>
% Create the H2-O2 gas phase
% Create the H2-O2 gas phase
gas = importPhase('h2o2.cti');
gas = importPhase('h2o2.cti');
Line 191: Line 182:
% Set the thermodynamic state of the gas
% Set the thermodynamic state of the gas
set(gas, 'T', 1500, 'P', oneatm, 'X', 'H2:2,O2:1,AR:20');
set(gas, 'T', 1500, 'P', oneatm, 'X', 'H2:2,O2:1,AR:20');
</source>
</pre>


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


<source lang="matlab">
<pre>
% Create an ideal gas mix of air
% Create an ideal gas mix of air
a = IdealGasMix('air.cti');
a = IdealGasMix('air.cti');
atmosphere = reservoir(a);
atmosphere = reservoir(a);
</source>
</pre>


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


<source lang="matlab">
<pre>
batch = Reactor(gas);
batch = Reactor(gas);
</source>
</pre>


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 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.
Line 212: Line 203:


<math>
<math>
\frac{dV}{dt} = K_{\mbox{expansion}} A_{\mbox{wall}} ( P_1 - P_2 )
\frac{dV}{dt} = K_{\text{expansion}} A_{\text{wall}} ( P_1 - P_2 )
</math>
</math>


Line 221: Line 212:
The larger the <math>K</math> value, the faster the response of the wall (piston).
The larger the <math>K</math> value, the faster the response of the wall (piston).


<source lang="matlab">
<pre>
% Create a wall between the reactor and the atmosphere
% Create a wall between the reactor and the atmosphere
w = Wall;
w = Wall;
Line 229: Line 220:
% Set the expansion coefficient K
% Set the expansion coefficient K
setExpansionRateCoeff(w, 1.0e10);
setExpansionRateCoeff(w, 1.0e10);
</source>
</pre>


Remember, we can also have walls with heat transfer, so we could have heat flowing in, or flowing out, through the wall
Remember, we can also have walls with heat transfer, so we could have heat flowing in, or flowing out, through the wall
Line 237: Line 228:
Next, a reactor network can be created, in order to advance the reactor in time.
Next, a reactor network can be created, in order to advance the reactor in time.


<source lang="matlab">
<pre>
network = ReactorNet( {batch} );
network = ReactorNet( {batch} );
</source>
</pre>


== CSTR Combustor ==
== CSTR Combustor ==
Line 261: Line 252:
Create a reservoir of fuel, a reservoir of air, and a reservoir of hydrogen radicals (radicals are sure to initiate combustion reactions).
Create a reservoir of fuel, a reservoir of air, and a reservoir of hydrogen radicals (radicals are sure to initiate combustion reactions).


<source lang="matlab">
<pre>
% Create a fuel from the GRI combustion mechanism
% Create a fuel from the GRI combustion mechanism
fuel = GRI30;
fuel = GRI30;
Line 282: Line 273:
set(igniter, 'T', 300, 'P', oneatm, 'X', 'H:1.0');
set(igniter, 'T', 300, 'P', oneatm, 'X', 'H:1.0');
igniter_in = Reservoir(igniter);
igniter_in = Reservoir(igniter);
</source>
</pre>


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


<source lang="matlab">
<pre>
inert = GRI30;
inert = GRI30;
set(inert, 'T', 300, 'P', oneatm, 'X', 'N2:1.0');
set(inert, 'T', 300, 'P', oneatm, 'X', 'N2:1.0');
combustor = Reactor(inert);
combustor = Reactor(inert);
</source>
</pre>


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


<source lang="matlab">
<pre>
ER = 0.7;
ER = 0.7;


Line 305: Line 296:
% Igniter mass flowrate
% Igniter mass flowrate
igniter_mdot = 0.05;
igniter_mdot = 0.05;
</source>
</pre>


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


<source lang="matlab">
<pre>
igniter_mdot = Func('gaussian', 0, [1.0, 0.0, 0.1] );
igniter_mdot = Func('gaussian', 0, [1.0, 0.0, 0.1] );
igniter_mdot = Func('polynomial', 2.0, [0.0, 0.0, 0.1] );
igniter_mdot = Func('polynomial', 2.0, [0.0, 0.0, 0.1] );
</source>
</pre>


You can type "help Func" for more information.
You can type "help Func" for more information.
Line 320: Line 311:
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.
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.


<source lang="matlab">
<pre>
% Create the mass flowrate controllers feeding into the combustor
% Create the mass flowrate controllers feeding into the combustor
m1 = MassFlowController(oxidizer_in, combustor);
m1 = MassFlowController(oxidizer_in, combustor);
Line 337: Line 328:
% Now the reactor network can be created
% Now the reactor network can be created
network = ReactorNet( {combustor} );
network = ReactorNet( {combustor} );
</source>
</pre>
 
=Python Examples=
 
==Before You Begin==
 
Before you start to use Cantera from Python, I highly recommend installing [[Python#Py4Sci|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:
 
<source lang="python">
from Cantera import *
from pylab import *
import sys
</source>
 
==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:
 
<source lang="python">
# 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()
 
</source>
 
Next, define the range of equivalence ratios over which to calculate the adiabatic flame temperature:
 
<source lang="python">
# equivalence ratio range
phi_min = 0.8
phi_max = 1.2
npoints = 10
</source>
 
It is useful to know the index of the fuel and oxidizer species, as well as nitrogen.  These can be obtained by using the <code>speciesIndex()</code> method of the GRI Mech phase created above:
 
<source lang="python">
# 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."
</source>
 
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):
 
<source lang="python">
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)
</source>
 
Oh, and to speed things up, you can pre-allocate some arrays to store values of <math>\Phi, T_{\mbox{adiabatic}}, x_i</math> (where the concentrations are the equilibrium concentrations)
 
<source lang="python">
# create some arrays to hold the data
phi = zeros(npoints,'d')
tad = zeros(npoints,'d')
xeq = zeros([nsp,npoints],'d')
</source>
 
===Calculation of Adiabatic Flame T===
 
The next section takes place inside a loop over all the values of equivalence ratios <math>\Phi</math> specified:
 
<source lang="python">
for i in range(npoints):
  phi[i] = phi_min + (phi_max - phi_min)*i/(npoints - 1)
</source>
 
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:
 
<source lang="python">
mix_phases = [ (gas, 1.0), (carbon, 0.0) ]
</source>
 
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:
 
<source lang="python">
  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)
</source>
 
Next, a Cantera Mixture object can be created, and its temperature and pressure set:
 
<source lang="python">
  # create a mixture of 1 mole of gas, and 0 moles of solid carbon.
  mix = Mixture(mix_phases)
  mix.setTemperature(temp)
  mix.setPressure(pres)
</source>
 
The mixture can then be equilibrated at constant pressure and enthalpy
 
<source>
  mix.equilibrate('HP', maxsteps = 1000,
                  err = 1.0e-6, maxiter = 200, loglevel=0)
</source>
 
Once the mixture is equilibrated, the temperature (the adiabatic flame temperature) can be obtained using the <code>temperature()</code> method:
 
<source>
  tad[i] = mix.temperature();
  print 'At phi = %12.4g, Tad = %12.4g'  % (phi[i],tad[i])
</source>
 
and the composition of the equilibrated phase can be obtained using the <code>speciesMoles()</code> method:
 
<source>
  xeq[:,i] = mix.speciesMoles()
</source>
 
===Outputting Data File===
 
Data generated can be output to a .csv file using the <code>writeCSV()</code> method:
 
<source lang="python">
# 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
</source>
 
===Create Plots===
 
There are three different plots that would be useful for the adiabatic flame temperature case:
* Adiabatic flame temperature vs. equivalence ratio <math>\Phi</math>
* Concentration of major species vs. equivalence ratio
* Concentration of minor species vs. equivalence ratio
 
Start by creating a list of species that are interesting:
 
<source lang="python">
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')
</source>
 
and define a mole fraction limit to distinguish major and minor species:
 
<source lang="python">
xlim = 0.04
</source>
 
Create a figure for each of the plots:
 
<source lang="python">
fig=figure(1)
</source>
 
You can use [[Matplotlib]] to create the plot of temperature vs. equivalence ratio. The syntax is very similar to Matlab.
 
<source lang="python">
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
</source>
 
The next plot is concentration of major species as a function of equivalence ratio:
 
<source>
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))
</source>
 
The last plot is concentration of minor species:
 
<source lang="python">
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))
</source>
 
Finally, a title can be added to the plot:
 
<source lang="python">
fig.text(0.5,0.95,r'Adiabatic Flame Temperature for $C_{3}H_{8}$ + Air', ...
fontsize=22,horizontalalignment='center')
 
show()
</source>
 
This results in the following plot:
 
[[Image:adiabatic_T_plot_C3H8.png|600px]]
 
 
 
 




{{Combustion}}




[[Category:Cantera]]
{{CanteraFlag}}

Latest revision as of 23:08, 14 June 2026


As a TA for CHEN 6153, I gave a lecture on the installation and use of Cantera.

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)

Batch.jpg

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'.

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

Piston.jpg

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_{\text{expansion}} A_{\text{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.

Cstr.jpg

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} );