Using Reactor Networks for Time Integration#
This guide explains the basics of setting up and integrating the governing equations of a transient reactor network problem.
Setting up a Reactor Network#
First, let’s take a look at a basic example to see how we might utilize Cantera’s time integration functionality. We’ll simulate an isolated reactor in Python that is filled with a homogeneous gas mixture. The gas state used in this example is arbitrary, but interesting because it’s explosive.
# import the Cantera Python module
import cantera as ct
# Create a gas mixture using the GRI-Mech 3.0 mechanism.
# This mechanism is included with Cantera for use in examples, but is not
# recommended for research use.
gas = ct.Solution("gri30.yaml")
# set gas to an interesting state
gas.TPX = 1000.0, ct.one_atm, "H2:2,O2:1,N2:4"
# create a reactor containing the gas
reac = ct.IdealGasReactor(gas)
# Add the reactor to a new ReactorNet simulator
sim = ct.ReactorNet([reac])
# View the initial state of the mixture
reac.thermo()
gri30:
temperature 1000 K
pressure 1.0132e+05 Pa
density 0.25781 kg/m^3
mean mol. weight 21.155 kg/kmol
phase of matter gas
1 kg 1 kmol
--------------- ---------------
enthalpy 1.0127e+06 2.1423e+07 J
internal energy 6.1963e+05 1.3108e+07 J
entropy 10427 2.2058e+05 J/K
Gibbs function -9.414e+06 -1.9915e+08 J
heat capacity c_p 1527.9 32322 J/K
heat capacity c_v 1134.9 24008 J/K
mass frac. Y mole frac. X chem. pot. / RT
--------------- --------------- ---------------
H2 0.027227 0.28571 -18.758
O2 0.21608 0.14286 -28.512
N2 0.7567 0.57143 -25.41
[ +50 minor] 0 0
Now, let’s advance the simulation in time to an absolute time of 1 second, and examine how the state of the gas has changed. There isn’t anything special about choosing 1 second as the end time for the integration; any end time can be chosen. We are choosing 1 second here because after that amount of time, the mixture will have almost certainly ignited for a wide range of initial conditions.
# advance the simulation to the specified absolute time, t = 1 sec
sim.advance(1)
# view the updated state of the mixture, reflecting properties at t = 1 sec
reac.thermo()
gri30:
temperature 2867.2 K
pressure 2.5924e+05 Pa
density 0.25781 kg/m^3
mean mol. weight 23.708 kg/kmol
phase of matter gas
1 kg 1 kmol
--------------- ---------------
enthalpy 1.6252e+06 3.8529e+07 J
internal energy 6.1963e+05 1.469e+07 J
entropy 11252 2.6677e+05 J/K
Gibbs function -3.0638e+07 -7.2637e+08 J
heat capacity c_p 1759.3 41709 J/K
heat capacity c_v 1408.6 33395 J/K
mass frac. Y mole frac. X chem. pot. / RT
--------------- --------------- ---------------
H2 0.0035836 0.042142 -22.913
H 0.00056302 0.013242 -11.457
O 0.0032005 0.0047426 -16.811
O2 0.015733 0.011657 -33.621
OH 0.017396 0.02425 -28.267
H2O 0.19703 0.25929 -39.724
HO2 1.0699e-05 7.6852e-06 -45.078
H2O2 9.0899e-07 6.3357e-07 -56.535
N 1.6832e-06 2.8488e-06 -13.859
NH 3.6478e-07 5.7596e-07 -25.316
NH2 1.1008e-07 1.6288e-07 -36.772
NH3 9.8884e-08 1.3765e-07 -48.229
NNH 9.2324e-08 7.5418e-08 -39.174
NO 0.010841 0.0085656 -30.67
NO2 4.0165e-06 2.0698e-06 -47.48
N2O 1.2929e-06 6.9643e-07 -44.528
HNO 1.3651e-06 1.0435e-06 -42.126
N2 0.75163 0.63609 -27.718
[ +35 minor] -3.3465e-25 -2.2393e-25
Methods for Time Integration#
The state of a ReactorNet can be advanced in time by one of the following
three methods:
step(): Thestep()method computes the state of the system after one time step. The size of the step is determined by SUNDIALS when the method is called. SUNDIALS determines the step size by estimating the local error, which must satisfy tolerance conditions. The step is redone with reduced step size whenever that error test fails. SUNDIALS also periodically checks if the maximum step size is being used. The time step must not be larger than a predefined maximum time step \(\Delta t_\t{max}\). The new time \(t_\t{new}\) at the end of the single step is returned by this function. This method produces the highest time resolution in the output data of the methods implemented in Cantera.advance(t_new): This method computes the state of the system at the user-provided time \(t_\t{new}\). \(t_\t{new}\) is the absolute time from the initial time of the system. Although the user specifies the time when integration should stop, SUNDIALS chooses the time step size as the network is integrated. Many internal SUNDIALS time steps are usually required to reach \(t_\t{new}\). As such,advance(t_new)preserves the accuracy of usingstep()but allows consistent spacing in the output data.advance_to_steady_state(max_steps, residual_threshold, atol, write_residuals)Python interface only: If the steady state solution of a reactor network is of interest, this method can be used. Internally, the steady state is approached by time stepping. The network is considered to be at steady state if the feature-scaled residual of the state vector is below a given threshold value (which by default is 10 times the time steprtol).
Limiting Timesteps#
The advance(t_new) approach is typically used when consistent, regular, time steps are
required in the output data. This usually simplifies comparisons among many simulation
results at a single time point. However, some detail, for example, a fast ignition
process, might not be resolved in the output data if the user-provided time step is not
small enough.
To avoid losing this detail, the Reactor::setAdvanceLimit() method (C++) or the
Reactor.set_advance_limit() method (Python) can be used to set the maximum
amount that a specified solution component can change between output times. For an
example of this feature’s use, see the example
reactor1.py. However, as a tradeoff, the time step
sizes in the output data are no longer guaranteed to be uniform.
Setting Tolerances#
Even though Cantera comes pre-defined with reasonable default values for tolerances and the maximum internal time step, the solver may not be correctly configured for the specific system. In this case SUNDIALS may stop with an error. To solve this problem, three parameters can be tuned:
the absolute tolerances
the relative tolerances
the maximum time step size
Reducing the third value is particularly useful when dealing with abrupt changes in the boundary conditions. One example of this is opening or closing valves; see also the IC engine example.
Adding Inlets and Outlets#
The following examples demonstrate the use of flow devices such as valves and mass flow controllers to model the inlets and outlets of reactors:
Modeling Volume Changes and Heat Transfer#
The following examples demonstrate the use of walls to model volume changes and heat transfer between reactors:
Additional examples can be found in the Python Reactor Network Examples section.
Accelerating Integration with Preconditioned Sparse Iterative Methods#
The Ideal Gas Control Volume Mole Reactor and the Ideal Gas Constant Pressure Mole Reactor are able to be solved using preconditioned sparse iterative methods, instead of the normal direct linear methods. Using these solvers can greatly speed up integration for systems with hundreds or thousands of species [Walker et al., 2023].
To enable the use of this solver, a preconditioner object needs to be created and specified for use by the reactor network. The previous example can be modified as follows:
import cantera as ct
gas = ct.Solution("gri30.yaml")
gas.TPX = 1000.0, ct.one_atm, "H2:2,O2:1,N2:4"
reac = ct.IdealGasMoleReactor(gas) # Use reactor type that supports preconditioning
sim = ct.ReactorNet([reac])
precon = ct.AdaptivePreconditioner() # Create the preconditioner
sim.preconditioner = precon # Add it to the network
The results of the time integration are the same as before:
sim.advance(1)
reac.thermo()
gri30:
temperature 2867.2 K
pressure 2.5924e+05 Pa
density 0.25781 kg/m^3
mean mol. weight 23.708 kg/kmol
phase of matter gas
1 kg 1 kmol
--------------- ---------------
enthalpy 1.6252e+06 3.853e+07 J
internal energy 6.1963e+05 1.469e+07 J
entropy 11252 2.6677e+05 J/K
Gibbs function -3.0638e+07 -7.2637e+08 J
heat capacity c_p 1759.3 41709 J/K
heat capacity c_v 1408.6 33395 J/K
mass frac. Y mole frac. X chem. pot. / RT
--------------- --------------- ---------------
H2 0.0035836 0.042142 -22.913
H 0.00056302 0.013242 -11.457
O 0.0032005 0.0047426 -16.811
O2 0.015733 0.011657 -33.621
OH 0.017396 0.02425 -28.267
H2O 0.19703 0.25929 -39.724
HO2 1.0699e-05 7.6852e-06 -45.078
H2O2 9.0899e-07 6.3357e-07 -56.535
N 1.6832e-06 2.8488e-06 -13.859
NH 3.6478e-07 5.7596e-07 -25.316
NH2 1.1008e-07 1.6288e-07 -36.772
NH3 9.8884e-08 1.3765e-07 -48.229
NNH 9.2324e-08 7.5418e-08 -39.174
NO 0.010841 0.0085656 -30.67
NO2 4.0165e-06 2.0698e-06 -47.48
N2O 1.2929e-06 6.9643e-07 -44.528
HNO 1.3651e-06 1.0435e-06 -42.126
N2 0.75163 0.63609 -27.718
[ +35 minor] 0 0
The approximate Jacobian matrices used to construct the preconditioners do not currently account for terms that link multiple reactors in the same network. For networks of this type, the iterative solver may exhibit convergence errors for systems where the default direct solver does not. If the solver converges, the solution will satisfy the error tolerances, even though the Jacobian matrix is not exact.
Common Reactor Types and their Implementation in Cantera#
Batch Reactor at Constant Volume or at Constant Pressure#
If you are interested in how a homogeneous chemical mixture changes state in time when isolated from the surroundings, a simple batch reactor can be used. Two versions are commonly considered: a rigid vessel with fixed volume but variable pressure, or a system idealized at constant pressure but varying volume.
The initial state of the solution can be specified by composition and a set of intensive
thermodynamic parameters, like temperature and pressure, as a standard Cantera
Solution object. From this base, a Reactor or a
ConstPressureReactor can be created, depending on if a constant volume or
constant pressure batch reactor should be considered, respectively. The behavior of the
solution in time can be simulated as a ReactorNet containing only the single
reactor.
An example for such a batch reactor is given in
reactor1.py.
Continuously Stirred Tank Reactor#
A Continuously Stirred Tank Reactor (CSTR), also often referred to as Well-Stirred Reactor (WSR), Perfectly Stirred Reactor (PSR), or Longwell Reactor, is essentially a single Cantera reactor with an inlet, an outlet, and constant volume.
Steady state solutions to CSTRs are often of interest. In this case, the mass flow rate \(\dot{m}\) is constant and equal at inlet and outlet. The mass contained in the control volume, \(m\), divided by \(\dot{m}\) defines the mean residence time of the fluid in the control volume.
While Cantera always solves a transient problem, if you are interested in steady-state
conditions, you can run your simulation for a long time until the states are converged.
This integration process can be managed automatically using the
ReactorNet.advance_to_steady_state() method. See
combustor.py for an example using this method.
Tip
A problem can be the ignition of a CSTR: If the reactants are not reactive enough, the
simulation can result in the trivial solution that inflow and outflow states are
identical. To solve this problem, the reactor can be initialized with a high temperature
and/or radical concentration. A good approach is to use the equilibrium composition of
the reactants (which can be computed using Cantera’s ThermoPhase.equilibrate()
function) as an initial guess. This method is demonstrated in the
combustor.py example.