research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Doc] Improve documentation for reactor networks and add more examples

Browse source 
Resolves #281
This commit is contained in:
Thomas Fiala 2015-06-29 17:05:54 +02:00 • committed by Ray Speth
parent ce21b08d41
commit b7f2451f06
3 changed files with 921 additions and 91 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

575
doc/sphinx/reactors.rst
View file

@ -1,74 +1,64 @@
.. default-role:: math
****************
Reactor Networks
****************
.. py:currentmodule:: cantera
Cantera's Reactor Network module is designed to simulate networks of
interconnected reactors. The contents of each reactor in the network are
assumed to be homogeneous, a model variously referred to as the Continuously
Stirred Tank Reactor (CSTR), Well-Stirred Reactor (WSR), or Perfectly Stirred
Reactor (PSR) model. Cantera solves the time-dependent governing equations
that describe the evolution of the chemical and thermodynamic state of the
reactors.
*****************************
Reactors and Reactor Networks
*****************************
The contents of each reactor can undergo chemical reactions according to a
specified kinetic mechanism, and surface reactions may occur on the reactor
walls. Each reactor in a network may be connected so that the contents of one
reactor flow into another. Reactors may be also be in contact with one another
or the environment via walls which move or conduct heat.
A Cantera Reactor represents the simplest form of a chemically reacting system.
It corresponds to an extensive thermodynamic control volume `V`, in which all
state variables are homogeneously distributed. The system is generally unsteady,
i.e. all states are functions of time. In particular, transient state changes
due to chemical reactions are possible. However, thermodynamic (but not
chemical) equilibrium is assumed to be present throughout the reactor at all
instants of time.
The purpose of this document is to describe the governing equations of reactor
models as implemented in Cantera.
Reactors can interact with the surrounding environment in multiple ways:
Wall Interactions
=================
- Expansion/compression work: By moving the walls of the reactor, its volume can
be changed and expansion or compression work can be done by or on the system,
i.e., the Reactor.
- Heat transfer: An arbitrary heat transfer rate can be defined to cross the
boundaries of the reactor.
- Mass transfer: The reactor can have multiple inlets and outlets. For the
inlets, arbitrary states can be defined. Through the outlets, fluid with the
current state of the reactor exits the reactor.
- Surface interaction: One or multiple walls can influence the chemical
reactions in the reactor. This is not just restricted to catalytic reactions,
but mass transfer between the surface and the fluid can also be modeled.
At each wall where there are surface reactions, there is a net generation (or
destruction) of homogeneous phase species. The molar rate of production for
each species `k` on wall `w` is `\dot{s}_{k,w}`. The total (mass) production
rate for species `k` on all walls is:
All of these interactions do not have to be constant, but can vary as a function
of time or state. For example, heat transfer can be described as a function of
the temperature difference between the reactor and the environment, or the wall
movement can be modeled depending on the pressure difference. Typically,
interactions of the reactor with the environment are defined on one or multiple
*walls*, *inlets*, and *outlets*.
.. math::
In addition to single reactors, Cantera is also able to interconnect reactors
into a *Reactor Network*. Each reactor in a network may be connected so that
the contents of one reactor flow into another. Reactors may also be in contact
with one another or the environment via walls which move or conduct heat.
\dot{m}_{k,wall} = W_k \sum_w A_w \dot{s}_{k,w}
where `W_k` is the molecular weight of species `k` and `A_w` is the area of
each wall. The net mass flux from all walls is then:
.. math::
\dot{m}_{wall} = \sum_k \dot{m}_{k,wall}
The total rate of heat transfer through all walls is:
.. math::
\dot{Q} = \sum_w f_w \dot{Q}_w
General Reactor
===============
Governing Equations for Single Reactors
=======================================
The state variables for Cantera's general reactor model are
- `m`, the mass of the reactor's contents
- `V`, the reactor volume
- `U`, the total internal energy of the reactors contents
- `Y_k`, the mass fractions for each species
- `m`, the mass of the reactor's contents (in kg)
- `V`, the reactor volume (in m\ :sup:`3`) (not a state variable for
*Constant Pressure Reactor* and *Ideal Gas Constant Pressure Reactor*)
- A state variable describing the energy of the system, depending on the
configuration (see `Energy Conservation`_ for further explanation):
Reactor Volume
--------------
- General *Reactor*: `U`, the total internal energy of the reactors
contents (in J)
- *Constant Pressure Reactor*: `H`, the total enthalpy of the reactors
contents (in J)
- *Ideal Gas Reactor* and *Ideal Gas Constant Pressure Reactor*: `T`, the
temperature (in K)
The reactor volume changes as a function of time due to the motion of one or
more walls:
.. math::
\frac{dV}{dt} = \sum_w f_w A_w v_w(t)
where `f_w = \pm 1` indicates the facing of the wall, `A_w` is the surface
area of the wall, and `v_w(t)` is the velocity of the wall as a function of
time.
- `Y_k`, the mass fractions for each species (dimensionless)
Mass Conservation
-----------------
@ -86,8 +76,8 @@ Species Conservation
--------------------
The rate at which species `k` is generated through homogeneous phase reactions
is `V \dot{\omega}_k W_k`, and the total rate at which species `k` is
generated is:
is `V \dot{\omega}_k W_k`, and the total rate at which species `k` is generated
is:
.. math::
@ -106,12 +96,40 @@ for `dm/dt`, the equation for each homogeneous phase species is:
.. math::
m \frac{dY}{dt} = \sum_{in} \dot{m}_{in} (Y_{k,in} - Y_k)+
m \frac{dY_k}{dt} = \sum_{in} \dot{m}_{in} (Y_{k,in} - Y_k)+
\dot{m}_{k,gen} - Y_k \dot{m}_{wall}
Reactor Volume
--------------
The reactor volume changes as a function of time due to the motion of one or
more walls:
.. math::
\frac{dV}{dt} = \sum_w f_w A_w v_w(t)
where `f_w = \pm 1` indicates the facing of the wall, `A_w` is the surface
area of the wall, and `v_w(t)` is the velocity of the wall as a function of
time.
For *Constant Pressure Reactor* and *Ideal Gas Constant Pressure Reactor*, the
volume is not a state variable, but instead takes on whatever value is
consistent with holding the pressure constant.
Energy Conservation
-------------------
The solution of the energy equation can be enabled or disabled by changing the
``energy_enabled`` flag. It is enabled by default.
The implemented formulation of the energy equation depends on which reactor
model is used.
Standard Reactor
****************
The equation for the total internal energy is found by writing the first law
for an open system:
@ -120,13 +138,36 @@ for an open system:
\frac{dU}{dt} = - p \frac{dV}{dt} - \dot{Q} +
\sum_{in} \dot{m}_{in} h_{in} - h \sum_{out} \dot{m}_{out}
Ideal Gas Reactor
=================
The Ideal Gas Reactor model is similar to the General Reactor model, with the
reactor temperature `T` replacing the total internal energy `U` as a state
variable. For an ideal gas, we can rewrite the total internal energy in terms
of the mass fractions and temperature:
Constant Pressure Reactor
*************************
For this reactor model, the pressure is held constant. The volume is not a
state variable, but instead takes on whatever value is consistent with holding
the pressure constant. The total enthalpy replaces the total internal energy
as a state variable. Using the definition of the total enthalpy:
.. math::
H = U + pV
\frac{d H}{d t} = \frac{d U}{d t} + p \frac{dV}{dt} + V \frac{dp}{dt}
Noting that `dp/dt = 0` and substituting into the energy equation yields:
.. math::
\frac{dH}{dt} = - \dot{Q} + \sum_{in} \dot{m}_{in} h_{in}
- h \sum_{out} \dot{m}_{out}
Ideal Gas Reactor
*****************
In case of the Ideal Gas Reactor Model, the reactor temperature `T` is used
instead of the total internal energy `U` as a state variable. For an ideal gas,
we can rewrite the total internal energy in terms of the mass fractions and
temperature:
.. math::
@ -150,34 +191,11 @@ significantly reduces the cost of evaluating the system Jacobian, since the
derivatives of the species equations are taken at constant temperature instead
of constant internal energy.
Constant Pressure Reactor
=========================
For this reactor model, the pressure is held constant. The volume is not a
state variable, but instead takes on whatever value is consistent with holding
the pressure constant. The total enthalpy replaces the total internal energy
as a state variable. Using the definition of the total enthalpy:
.. math::
H = U + pV
\frac{dH}{dt} = p \frac{dV}{dt} + V \frac{dp}{dt}
Noting that `dp/dt = 0` and substituting into the energy equation yields:
.. math::
\frac{dH}{dt} = - \dot{Q} + \sum_{in} \dot{m}_{in} h_{in}
- h \sum_{out} \dot{m}_{out}
The species and continuity equations are the same as for the general reactor
model.
Ideal Gas Constant Pressure Reactor
===================================
***********************************
As for the Ideal Gas Reactor, we replace the total enthalpy as a state
As for the Ideal Gas Reactors, we replace the total enthalpy as a state
variable with the temperature by writing the total enthalpy in terms of the
mass fractions and temperature:
@ -195,3 +213,378 @@ temperature:
m c_p \frac{dT}{dt} = - \dot{Q} - \sum_k h_k \dot{m}_{k,gen}
+ \sum_{in} \dot{m}_{in} \left(h_{in} - \sum_k h_k Y_{k,in} \right)
Wall Interactions
-----------------
The total rate of heat transfer through all walls is:
.. math::
\dot{Q} = \sum_w f_w \dot{Q}_w
where `f_w = \pm 1` indicates the facing of the wall.
In case of surface reactions, there is a net generation (or
destruction) of homogeneous phase species at the wall. The molar rate of
production for each species `k` on wall `w` is `\dot{s}_{k,w}` (in kmol/s/m\
:sup:`2`). The total (mass) production rate for species `k` on all walls is:
.. math::
\dot{m}_{k,wall} = W_k \sum_w A_w \dot{s}_{k,w}
where `W_k` is the molecular weight of species `k` and `A_w` is the area of
each wall. The net mass flux from all walls is then:
.. math::
\dot{m}_{wall} = \sum_k \dot{m}_{k,wall}
Reactor Networks and Devices
============================
While reactors by themselves just define the above governing equations of the
reactor, the time integration is performed in reactor networks. A reactor
network is therefore necessary even if only a single reactor is considered.
The advantage of reactor networks obviously is that multiple reactors can be
interconnected. Not only mass flow from one reactor into another can be
realized, but also heat can be transferred, or the wall between reactors can
move. To set up a network, the following components can be defined in addition
to the reactors previously mentioned:
- **Reservoir**: A reservoir can be thought of as an infinitely large volume, in
which all states are predefined and never change from their initial values.
Typically, it represents a vessel to define temperature and composition of a
stream of mass flowing into a reactor, or the ambient fluid surrounding the
reactor network. Besides, the fluid flow finally finally exiting a reactor
network has to flow into a reservoir. In the latter case, the state of the
reservoir (except pressure) is irrelevant.
- **Wall**: A wall separates two reactors, or a reactor and a reservoir. A wall
has a finite area, may conduct or radiate heat between the two reactors on
either side, and may move like a piston.
Walls are stateless objects in Cantera, meaning that no differential equation
is integrated to determine any wall property. Since it is the wall (piston)
velocity that enters the energy equation, this means that it is the velocity,
not the acceleration or displacement, that is specified. The wall velocity is
computed from
.. math:: v = K(P_{\rm left} - P_{\rm right}) + v_0(t),
where :math:`K` is a non-negative constant, and :math:`v_0(t)` is a specified
function of time. The velocity is positive if the wall is moving to the right.
The heat flux through the wall is computed from
.. math:: q = U(T_{\rm left} - T_{\rm right}) + \epsilon\sigma (T_{\rm left}^4
- T_{\rm right}^4) + q_0(t),
where :math:`U` is the overall heat transfer coefficient for
conduction/convection, and :math:`\epsilon` is the emissivity. The function
:math:`q_0(t)` is a specified function of time. The heat flux is positive when
heat flows from the reactor on the left to the reactor on the right.
A heterogeneous reaction mechanism may be specified for one or both of the
wall surfaces. The mechanism object (typically an instance of class Interface)
must be constructed so that it is properly linked to the object representing
the fluid in the reactor the surface in question faces. The surface
temperature on each side is taken to be equal to the temperature of the
reactor it faces.
Source: `Python <cython/zerodim.html#wall>`_ | :ct:`C++ <Wall>`
- **Valve**: A valve is a flow devices with mass flow rate that is a function of
the pressure drop across it. The default behavior is linear:
.. math:: \dot m = K_v (P_1 - P_2)
if :math:`P_1 > P_2.` Otherwise, :math:`\dot m = 0`. However, an arbitrary
function can also be specified, such that
.. math:: \dot m = F(P_1 - P_2)
if :math:`P_1 > P_2`, or :math:`\dot m = 0` otherwise. It is never possible
for the flow to reverse and go from the downstream to the upstream
reactor/reservoir through a line containing a Valve object.
Valve objects are often used between an upstream reactor and a downstream
reactor or reservoir to maintain them both at nearly the same pressure. By
setting the constant :math:`K_v` to a sufficiently large value, very small
pressure differences will result in flow between the reactors that counteracts
the pressure difference.
- **Mass Flow Controller**: A mass flow controller maintains a specified mass
flow rate independent of upstream and downstream conditions. The equation used
to compute the mass flow rate is
.. math:: \dot m = \max(\dot m_0, 0.0)
where :math:`\dot m_0` is either a constant value or a function of time. Note
that if :math:`\dot m_0 < 0`, the mass flow rate will be set to zero, since
reversal of the flow direction is not allowed.
Unlike a real mass flow controller, a MassFlowController object will maintain
the flow even if the downstream pressure is greater than the upstream
pressure. This allows simple implementation of loops, in which exhaust gas
from a reactor is fed back into it through an inlet. But note that this
capability should be used with caution, since no account is taken of the work
required to do this.
- **Pressure Controller**: A pressure controller is designed to be used in
conjunction with another 'master' flow controller, typically a
MassFlowController. The master flow controller is installed on the inlet of
the reactor, and the corresponding PressureController is installed on on
outlet of the reactor. The PressureController mass flow rate is equal to the
master mass flow rate, plus a small correction dependent on the pressure
difference:
.. math:: \dot m = \dot m_{\rm master} + K_v(P_1 - P_2).
Time Integration
----------------
Cantera provides an ODE solver for solving the stiff equations of reacting
systems. If installed in combination with SUNDIALS, their optimized solver is
used. Starting off the current state of the system, it can be advanced in time
by two methods:
- ``step()``: The step method computes the state of the system at the a priori
unspecified time `t_{\rm new}`. The time `t_{\rm new}` is internally computed
so that all states of the system only change within a (specifiable) band of
absolute and relative tolerances. Additionally, the time step must not be
larger than a predefined maximum time step `\Delta t_{\rm max}`. The new time
`t_{\rm new}` is returned by this function.
- ``advance``\ `(t_{\rm new})`: This method computes the state of the system at
time `t_{\rm new}`. `t_{\rm new}` describes the absolute time from the initial
time of the system. By calling this method in a for loop for pre-defined
times, the state of the system is obtained for exactly the times specified.
Internally, several ``step()`` calls are typically performed to reach the
accurate state at time `t_{\rm new}`.
The use of the ``advance`` method in a loop has the advantage that it produces
results corresponding to a predefined time series. These are associated with a
predefined memory consumption and well comparable between simulation runs with
different parameters. However, some detail (e.g. a fast ignition process) might
not be resolved in the output data due to the typically large time steps.
The ``step`` method results in much more data points because of the small
timesteps needed. Additionally, the absolute time has to be kept tracked of
manually.
Even though Cantera comes pre-defined with typical parameters for tolerances
and the maximum internal time step, the solution sometimes diverges. To solve
this problem, three parameters can be tuned: The absolute time stepping
tolerances, the relative time stepping tolerances, and the maximum time step. A
reduction of the latter value is particularly useful when dealing with abrupt
changes in the boundary conditions (e.g. opening/closing valves, see also
example :ref:`py-example-ic_engine.py`).
General Usage in Cantera
========================
In Cantera, the following steps are typically necessary to investigate a
reactor network:
1. Define ``Solution`` objects for the fluids to be flowing through your
reactor network.
2. Define the reactor type(s) and reservoir(s) that describe your system. Chose
Ideal Gas (Constant Pressure) Reactor(s) if you only consider ideal gas phases.
3. *Optional:* Set up the boundary conditions and flow devices between reactors
or reservoirs.
4. Define a reactor network which contains all the reactors previously created.
5. Advance the simulation in time, typically in a for- or while-loop. Note that
only the current state is stored in Cantera by default. If you want to observe
the transient states, you manually have to keep track of them.
6. Analyze the data.
Note that 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 (see e.g. example :ref:`py-example-surf_pfr.py`,
:ref:`py-example-combustor.py`).
Cantera comes with a broad variety of well-commented example scrips for reactor
networks. Please refer to them for further information (:ref:`Python <sec-cython-examples>`, :ref:`Matlab <sec-matlab-examples>`).
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 composition changes in time
when it is left to its own, 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.
In Cantera, such a simulation can be performed very easily. The initial state
of the solution can be specified by composition and a set of thermodynamic
parameters (like temperature and pressure) as a standard Cantera solution
object. Upon its base, a general (Ideal Gas) Reactor or an (Ideal Gas) Constant
Pressure Reactor 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 very simple Reactor Network containing
only the formerly created reactor.
An example for such a Batch Reactor is :ref:`py-example-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. Therefore, the `Governing Equations for Single Reactors`_
defined above apply accordingly.
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 confinement `m` divided by `\dot{m}` defines the mean
residence time of the fluid in the confinement.
At steady state, the time derivatives in the governing equations become zero,
and the system of ordinary differential equations can be reduced to a set of
coupled nonlinear algebraic equations. A Newton solver could be used to solve
this system of equations. However, a sophisticated implementation might be
required to account for the strong nonlinearities and the presence of multiple
solutions.
Cantera does not have such a Newton solver implemented. Instead, steady CSTRs
are simulated by considering a time-dependent constant volume reactor with
specified in- and outflow conditions. Starting off at an initial solution, the
reactor network containing this reactor is advanced in time until the state of
the solution is converged. An example for this procedure is
:ref:`py-example-combustor.py`.
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 ``equilibrate`` function) as an initial guess.
Plug-Flow Reactor
-----------------
A Plug-Flow Reactor (PFR) represents a steady-state channel with a
cross-sectional area `A`. Typically an ideal gas flows through it at a constant
mass flow rate `\dot{m}`. Perpendicular to the flow direction, the gas is
considered to be completely homogeneous. In the axial direction `z`, the states
of the gas is allowed to change. However, all diffusion processes are neglected.
Plug-Flow Reactors are often used to simulate ignition delay times, emission
formation, and catalytic processes.
The governing equations of Plug-Flow Reactors are [KCG2003]_:
- Mass conservation:
.. math:: \frac{d(\rho u A)}{dz} = P' \sum_k \dot{s}_k W_k
where `u` is the axial velocity in (m/s) and `P'` is the chemically active
channel perimeter in (m) (chemically active perimeter per unit length).
- Continuity equation of species `k`:
.. math:: \rho u \frac{d Y_k}{dz} + Y_k P' \sum_k \dot{s}_k W_k =
\dot{\omega}_k W_k + P' \dot{s}_k W_k
- Energy conservation:
.. math:: \rho u A c_p \frac{d T}{d z} =
- A \sum_k h_k \dot{\omega}_k W_k
- P' \sum_k h_k \dot{s}_k W_k
+ U P (T_w - T)
where `U` is the heat transfer coefficient in (W/m/K), `P` is the perimeter of
the duct in (m), and `T_w` is the wall temperature in (K). Kinetic and
potential energies are neglected.
- Momentum conservation in the axial direction:
.. math:: \rho u A \frac{d u}{d z} + u P' \sum_k \dot{s}_k W_k =
- \frac{d (p A)}{dz} - \tau_w P
where `\tau_w` is the wall friction coefficient (which might be computed from
Reynolds number based correlations).
Even though this problem extends geometrically in one direction, it can be
modeled via zero-dimensional reactors: Due to the neglecting of diffusion,
downstream parts of the reactor have no influence on upstream parts. Therefore,
PFRs can be modeled by marching from the beginning to the end of the reactor.
Cantera does not (yet) provide dedicated class to solve the PFR equations (The
``FlowReactor`` class is currently under development). However, there are two
ways to simulate a PFR with the reactor elements previously presented. Both
rely on the assumption that pressure is approximately constant throughout the
Plug-Flow Reactor and that there is no friction. The momentum conservation
equation is thus neglected.
PFR Modeling by Considering a Lagrangian Reactor
************************************************
A Plug-Flow Reactor can also be described from a Lagrangian point of view: An
unsteady fluid particle is considered which travels along the axial streamline
through the PFR. Since there is no information traveling upstream, the state
change of the fluid particle can be computed by a forward (upwind) integration
in time. Using the continuity equation, the speed of the particle can be
derived. By integrating the velocity in time, the temporal information can be
translated into the spatial resolution of the PFR.
An example for this procedure can be found in :ref:`py-example-pfr.py`.
PFR Modeling as a Series of CSTRs
*********************************
The Plug-Flow Reactor is spatially discretized into a large number of axially
distributed volumes. These volumes are modeled to be steady-state CSTRs.
The only reason to use this approach as opposed to the Lagrangian one is if you
need to include surface reactions, because the system of equations ends up
being a DAE system instead of an ODE system.
In Cantera, it is sufficient to consider a single reactor and march it forward
in time, because there is no information traveling upstream. The mass flow rate
`\dot{m}` through the PFR enters the reactor from an upstream reservoir. For
the first reactor, the reservoir conditions are the inflow boundary conditions
of the PFR. By performing a time integration as described in `Continuously
Stirred Tank Reactor`_ until the state of the reactor is converged, the
steady-state CSTR solution is computed. The state of the CSTR is the inlet
boundary condition for the next CSTR downstream.
An example for this procedure can be found in :ref:`py-example-pfr.py` and
:ref:`py-example-surf_pfr.py`.
Advanced Concepts
=================
In some cases, Cantera's solver is insufficient to describe a certain
configuration. In this situation, Cantera can still be used to provide chemical
and thermodynamic computations, but external ODE solvers can be applied. See
example :ref:`py-example-custom.py`.
Literature
==========
For further reading, the following books are recommended:
.. [KCG2003] Kee, Coltrin, Glarborg: *Chemically Reacting Flow*.
Wiley-Interscience, 2003
.. [Tur2000] Turns: *An Introduction to Combustion: Concepts and Applications*,
McGraw Hill, 2000

263
interfaces/cython/cantera/examples/reactors/ic_engine.py Normal file
View file

@ -0,0 +1,263 @@
# -*- coding: utf-8 -*-
"""
Simulation of a (gaseous) Diesel-type internal combustion engine.
The use of pure propane as fuel requires an unrealistically high compression
ratio.
"""
import cantera as ct
import numpy as np
#######################################################################
# Input Parameters
#######################################################################
f = 3000. / 60. # engine speed [1/s] (3000 rpm)
V_H = .5e-3 # displaced volume [m**3]
epsilon = 50. # compression ratio [-]
d_piston = 0.083 # piston diameter [m]
# turbocharger temperature, pressure, and composition
T_inlet = 300. # K
p_inlet = 1.3e5 # Pa
comp_inlet = 'O2:1, N2:3.76'
# outlet pressure
p_outlet = 1.2e5 # Pa
# fuel properties (gaseous!)
T_injector = 300. # K
p_injector = 1600e5 # Pa
comp_injector = 'C3H8:1'
# ambient properties
T_ambient = 300. # K
p_ambient = 1e5 # Pa
comp_ambient = 'O2:1, N2:3.76'
# Reaction mechanism name
reaction_mechanism = 'gri30.xml'
# Inlet valve friction coefficient, open and close timings
inlet_valve_coeff = 1.e-6
inlet_open = -18. / 180. * np.pi
inlet_close = 198. / 180. * np.pi
# Outlet valve friction coefficient, open and close timings
outlet_valve_coeff = 1.e-6
outlet_open = 522. / 180 * np.pi
outlet_close = 18. / 180. * np.pi
# Fuel mass, injector open and close timings
injector_open = 350. / 180. * np.pi
injector_close = 365. / 180. * np.pi
injector_mass = 3.2e-5 # kg
injector_t_open = (injector_close - injector_open) / 2. / np.pi / f
# Simulation time and resolution
sim_n_revolutions = 8.
sim_n_timesteps = 100000.
###################################################################
# load reaction mechanism
gas = ct.Solution(reaction_mechanism)
# define initial state
gas.TPX = T_inlet, p_inlet, comp_inlet
r = ct.IdealGasReactor(gas)
# define inlet state
gas.TPX = T_inlet, p_inlet, comp_inlet
inlet = ct.Reservoir(gas)
# define injector state (gaseous!)
gas.TPX = T_injector, p_injector, comp_injector
injector = ct.Reservoir(gas)
# define outlet pressure (temperature and composition don't matter)
gas.TPX = T_ambient, p_outlet, comp_ambient
outlet = ct.Reservoir(gas)
# define ambient pressure (temperature and composition don't matter)
gas.TPX = T_ambient, p_ambient, comp_ambient
ambient_air = ct.Reservoir(gas)
# set up connecting devices
inlet_valve = ct.Valve(inlet, r)
injector_mfc = ct.MassFlowController(injector, r)
outlet_valve = ct.Valve(r, outlet)
piston = ct.Wall(ambient_air, r)
# convert time to crank angle
def crank_angle(t):
return np.remainder(2 * np.pi * f * t, 4 * np.pi)
# set up IC engine parameters
V_oT = V_H / (epsilon - 1.)
A_piston = .25 * np.pi * d_piston ** 2
stroke = V_H / A_piston
r.volume = V_oT
piston.area = A_piston
def piston_speed(t):
return - stroke / 2 * 2 * np.pi * f * np.sin(crank_angle(t))
piston.set_velocity(piston_speed)
# create a reactor network containing the cylinder
sim = ct.ReactorNet([r])
# set up output data arrays
t_sim = sim_n_revolutions / f
t = (np.arange(sim_n_timesteps) + 1) / sim_n_timesteps * t_sim
p = np.zeros_like(t)
V = np.zeros_like(t)
T = np.zeros_like(t)
s = np.zeros_like(t)
m = np.zeros_like(t)
test = np.zeros_like(t)
mdot_in = np.zeros_like(t)
mdot_out = np.zeros_like(t)
MW = np.zeros_like(t)
d_W_v_d_t = np.zeros_like(t)
heat_release_rate = np.zeros_like(t)
species_X = np.zeros((t.size, gas.n_species))
# set parameters for the automatic time step refinement
n_last_refinement = -np.inf # for initialization only
n_wait_coarsening = 10
# do simulation
for n1, t_i in enumerate(t):
# define opening and closing of valves and injector
if (np.mod(crank_angle(t_i) - inlet_open, 4 * np.pi) <
np.mod(inlet_close - inlet_open, 4 * np.pi)):
inlet_valve.set_valve_coeff(inlet_valve_coeff)
test[n1] = 1
else:
inlet_valve.set_valve_coeff(0)
if (np.mod(crank_angle(t_i) - outlet_open, 4 * np.pi) <
np.mod(outlet_close - outlet_open, 4 * np.pi)):
outlet_valve.set_valve_coeff(outlet_valve_coeff)
else:
outlet_valve.set_valve_coeff(0)
if (np.mod(crank_angle(t_i) - injector_open, 4 * np.pi) <
np.mod(injector_close - injector_open, 4 * np.pi)):
injector_mfc.set_mass_flow_rate(injector_mass / injector_t_open)
else:
injector_mfc.set_mass_flow_rate(0)
# perform time integration, refine time step if necessary
for n2 in range(4):
if n2 is 4:
raise 'Error: Refinement limit reached'
try:
sim.advance(t_i)
except Exception:
sim.set_max_time_step(1e-6 * 10. ** -n2)
n_last_refinement = n1
# coarsen time step if too long ago
if n1 - n_last_refinement is n_wait_coarsening:
sim.set_max_time_step(1e-5)
# write output data
p[n1] = r.thermo.P
V[n1] = r.volume
T[n1] = r.T
s[n1] = r.thermo.s
m[n1] = r.mass
mdot_in[n1] = inlet_valve.mdot(0)
mdot_out[n1] = outlet_valve.mdot(0)
MW[n1] = r.thermo.mean_molecular_weight
species_X[n1] = r.thermo.X
d_W_v_d_t[n1] = - (r.thermo.P - ambient_air.thermo.P) * A_piston * \
piston_speed(t_i)
heat_release_rate[n1] = - r.volume * ct.gas_constant * r.T * \
np.sum(gas.standard_enthalpies_RT * r.thermo.net_production_rates, 0)
#####################################################################
# Plot Results in matplotlib
#####################################################################
import matplotlib.pyplot as plt
# pressure and temperature
plt.figure()
plt.clf()
plt.subplot(211)
plt.plot(t, p / 1.e5)
plt.ylabel('$p$ [bar]')
plt.xlabel('$\phi$ [deg]')
plt.xticks(plt.xticks()[0], [])
plt.subplot(212)
plt.plot(t, T)
plt.ylabel('$T$ [K]')
plt.xlabel('$\phi$ [deg]')
plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
rotation=17)
plt.show()
plt.savefig('ic_engine_t_p_T.png')
# p-V diagram
plt.figure()
plt.clf()
plt.plot(V[t > 0.04] * 1000, p[t > 0.04] / 1.e5)
plt.xlabel('$V$ [l]')
plt.ylabel('$p$ [bar]')
plt.show()
plt.savefig('ic_engine_p_V.png')
# T-S diagram
plt.figure()
plt.clf()
plt.plot(m[t > 0.04] * s[t > 0.04], T[t > 0.04])
plt.xlabel('$S$ [J/K]')
plt.ylabel('$T$ [K]')
plt.show()
plt.savefig('ic_engine_T_S.png')
# heat of reaction and expansion work
plt.figure()
plt.clf()
plt.plot(t, heat_release_rate, label='$\dot{Q}$')
plt.plot(t, d_W_v_d_t, label='$\dot{W}_v$')
plt.ylim(-1e5, 1e6)
plt.legend(loc=0)
plt.ylabel('[W]')
plt.xlabel('$\phi$ [deg]')
plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
rotation=17)
plt.show()
plt.savefig('ic_engine_Q_W.png')
# gas composition
plt.figure()
plt.clf()
plt.plot(t, species_X[:, gas.species_index('O2')], label='O2')
plt.plot(t, species_X[:, gas.species_index('CO2')], label='CO2')
plt.plot(t, species_X[:, gas.species_index('CO')], label='CO')
plt.plot(t, species_X[:, gas.species_index('C3H8')] * 10, label='C3H8 x10')
plt.legend(loc=0)
plt.ylabel('$X_i$ [-]')
plt.xlabel('$\phi$ [deg]')
plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
rotation=17)
plt.show()
plt.savefig('ic_engine_t_X.png')
#####################################################################
# Integral Results
#####################################################################
from scipy.integrate import trapz
Q = trapz(heat_release_rate, t)
W = trapz(d_W_v_d_t, t)
eta = W / Q
CO_emission = trapz(MW * mdot_out * species_X[:, gas.species_index('CO')], t) \
/ trapz(MW * mdot_out, t)
print('Heat release rate per cylinder (estimate):\t' +
format(Q / t_sim / 1000., ' 2.1f') + ' kW')
print('Expansion power per cylinder (estimate):\t' +
format(W / t_sim / 1000., ' 2.1f') + ' kW')
print('Efficiency (estimate):\t\t\t' + format(eta * 100., ' 2.1f') + ' %')
print('CO emission (estimate):\t\t' + format(CO_emission * 1.e6, ' 2.1f') +
' ppm')

174
interfaces/cython/cantera/examples/reactors/pfr.py Normal file
View file

@ -0,0 +1,174 @@
# -*- coding: utf-8 -*-
"""
This example solves a plug-flow reactor problem of hydrogen-oxygen combustion.
The PFR is computed by two approaches: The simulation of a Lagrangian fluid
particle, and the simulation of a chain of reactors.
"""
import cantera as ct
import numpy as np
#######################################################################
# Input Parameters
#######################################################################
T_0 = 1500.0 # inlet temperature [K]
pressure = ct.one_atm # constant pressure [Pa]
composition_0 = 'H2:2, O2:1, AR:0.1'
length = 1.5e-7 # *approximate* PFR length [m]
u_0 = .006 # inflow velocity [m/s]
area = 1.e-4 # cross-sectional area [m**2]
# input file containing the reaction mechanism
reaction_mechanism = 'h2o2.xml'
# Resolution: The PFR will be simulated by 'n_steps' time steps or by a chain
# of 'n_steps' stirred reactors.
n_steps = 2000
#####################################################################
#####################################################################
# Method 1: Lagrangian Particle Simulation
#####################################################################
# A Lagrangian particle is considered which travels through the PFR. Its
# state change is computed by upwind time stepping. The PFR result is produced
# by transforming the temporal resolution into spatial locations.
# The spatial discretization is therefore not provided a priori but is instead
# a result of the transformation.
# import the gas model and set the initial conditions
gas1 = ct.Solution(reaction_mechanism)
gas1.TPX = T_0, pressure, composition_0
mass_flow_rate1 = u_0 * gas1.density * area
# create a new reactor
r1 = ct.IdealGasConstPressureReactor(gas1)
# create a reactor network for performing time integration
sim1 = ct.ReactorNet([r1])
# approximate a time step to achieve a similar resolution as in the next method
t_total = length / u_0
dt = t_total / n_steps
# define time, space, and other information vectors
t1 = (np.arange(n_steps) + 1) * dt
z1 = np.zeros_like(t1)
u1 = np.zeros_like(t1)
T1 = np.zeros_like(t1)
X_H2_1 = np.zeros_like(t1)
for n1, t_i in enumerate(t1):
# perform time integration
sim1.advance(t_i)
# compute velocity and transform into space
u1[n1] = mass_flow_rate1 / area / r1.thermo.density
z1[n1] = z1[n1 - 1] + u1[n1] * dt
# write output data
T1[n1] = r1.T
X_H2_1[n1] = r1.thermo['H2'].X
#####################################################################
#####################################################################
# Method 2: Chain of Reactors
#####################################################################
# The plug flow reactor is represented by a linear chain of zero-dimensional
# reactors. The gas at the inlet to the first one has the specified inlet
# composition, and for all others the inlet composition is fixed at the
# composition of the reactor immediately upstream. Since in a PFR model there
# is no diffusion, the upstream reactors are not affected by any downstream
# reactors, and therefore the problem may be solved by simply marching from
# the first to last reactor, integrating each one to steady state.
# (This approach is anologous to the one presented in 'surf_pfr.py', which
# additionally includes surface chemistry)
# import the gas model and set the initial conditions
gas2 = ct.Solution(reaction_mechanism)
gas2.TPX = T_0, pressure, composition_0
mass_flow_rate2 = u_0 * gas2.density * area
dz = length / n_steps
r_vol = area * dz
# create a new reactor
r2 = ct.IdealGasReactor(gas2)
r2.volume = r_vol
# create a reservoir to represent the reactor immediately upstream. Note
# that the gas object is set already to the state of the upstream reactor
upstream = ct.Reservoir(gas2, name='upstream')
# create a reservoir for the reactor to exhaust into. The composition of
# this reservoir is irrelevant.
downstream = ct.Reservoir(gas2, name='downstream')
# The mass flow rate into the reactor will be fixed by using a
# MassFlowController object.
m = ct.MassFlowController(upstream, r2, mdot=mass_flow_rate2)
# We need an outlet to the downstream reservoir. This will determine the
# pressure in the reactor. The value of K will only affect the transient
# pressure difference.
v = ct.PressureController(r2, downstream, master=m, K=1e-5)
sim2 = ct.ReactorNet([r2])
# define time, space, and other information vectors
z2 = (np.arange(n_steps) + 1) * dz
t_r2 = np.zeros_like(z2) # residence time in each reactor
t2 = np.zeros_like(z2)
u2 = np.zeros_like(z2)
T2 = np.zeros_like(z2)
X_H2_2 = np.zeros_like(z2)
# iterate through the PFR cells
for n in range(n_steps):
# Set the state of the reservoir to match that of the previous reactor
gas2.TDY = r2.thermo.TDY
upstream.syncState()
# integrate the reactor forward in time until steady state is reached
sim2.set_initial_time(0) # forces reinitialization
time = 0
all_done = False
# determine steady state from H2 mole fraction
X_H2_previous = r2.thermo['H2'].X
while not all_done:
time += dt
sim2.advance(time)
if np.abs(r2.thermo['H2'].X - X_H2_previous) < 1.e-10:
# check whether surface coverages are in steady state.
all_done = True
else:
X_H2_previous = r2.thermo['H2'].X
# compute velocity and transform into time
u2[n] = mass_flow_rate2 / area / r2.thermo.density
t_r2[n] = r2.mass / mass_flow_rate2 # residence time in this reactor
t2[n] = np.sum(t_r2)
# write output data
T2[n] = r2.T
X_H2_2[n] = r2.thermo['H2'].X
#####################################################################
#####################################################################
# Compare Results in matplotlib
#####################################################################
import matplotlib.pyplot as plt
plt.figure()
plt.plot(z1, T1, label='Lagrangian Particle')
plt.plot(z2, T2, label='Reactor Chain')
plt.xlabel('$z$ [m]')
plt.ylabel('$T$ [K]')
plt.legend(loc=0)
plt.show()
plt.savefig('pfr_T_z.png')
plt.figure()
plt.plot(t1, X_H2_1, label='Lagrangian Particle')
plt.plot(t2, X_H2_2, label='Reactor Chain')
plt.xlabel('$t$ [s]')
plt.ylabel('$X_{H_2}$ [-]')
plt.legend(loc=0)
plt.show()
plt.savefig('pfr_XH2_t.png')