III - 0D computations¶
The following link may help you to find the answers or the documentation you need :
https://cantera.org/documentation/docs-3.0/sphinx/html/cython/zerodim.html
https://cantera.org/science/reactors/reactors.html
The first one is about the function that can be used for computing a reactor and the second one explains the equations that are computed.
1. Knowledge of 0D computations¶
ReactorBase class¶
0D Reactors are often used in combustion to model auto-ignition time and the chemistry associated to it. For reactors and reservoirs, the class used is the following :
ReactorBase(ThermoPhase contents, name=None, arguments)
Here are the different types of reactors you can have for a 0D computation :
Reservoir
Reactors with a constant state. The temperature, pressure, and chemical composition in a reservoir never change from their initial values.Reactor
- IdealGasReactor
Constant volume, zero-dimensional reactor for ideal gas mixtures. - ConstPressureReactor
Homogeneous, constant pressure, zero-dimensional reactor. The volume of the reactor changes as a function of time in order to keep the pressure constant. - IdealGasConstPressureReactor
Homogeneous, constant pressure, zero-dimensional reactor for ideal gas mixtures. The volume of the reactor changes as a function of time in order to keep the pressure constant. - FlowReactor
A steady-state plug flow reactor with constant cross sectional area. Time integration follows a fluid element along the length of the reactor. The reactor is assumed to be frictionless and adiabatic.
FlowDevice class¶
Reactors can be added inlets and outlets, to be able to loop or chain several reactors, control mass flow rates, volume or pressure. Here are the properties of the common baseclass FlowDevice for flow controllers between reactors and reservoirs :
FlowDevice(upstream, downstream, name=None, arguments)
The FlowDevice controls the fluids passage between an upstream and a downstream object, which should be specified. The arguments depend upon the different types of flow controllers available in Cantera :
- MassFlowController
Mass flow controllers maintain a specified mass flow rate independent of upstream and downstream conditions. Unlike a real mass flow controller, a MassFlowController object will maintain the flow even if the downstream pressure is greater than the upstream pressure. - Valve
Valves are flow devices with mass flow rate that is a function of the pressure drop across them. 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 Kv to a sufficiently large value, very small pressure differences will result in flow between the reactors that counteracts the pressure difference. - PressureController
A pressure controller is designed to be used in conjunction, typically, with a MassFlowController. The master flow controller is installed on the inlet of the reactor, and the corresponding PressureController is installed on the 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.
Walls¶
Reactors can also add heat transfer and heterogeneous reactions at the walls, through a special object "wall". Walls separate 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. They are stateless objects in Cantera, meaning that no differential equation is integrated to determine any wall property. A heterogeneous reaction mechanism may be specified for one or both of the wall surfaces.
ReactorNet¶
Now, the evolution of a reactor, or a network of reactors -which are 0-D objects- with time, is performed
through an integrator object ReactorNet :
... (define gases) ...
r1 = Reactor(gas1)
r2 = Reactor(gas2)
... (install walls, inlets, outlets, etc)...
reactor_network = ReactorNet([r1, r2])
time = 1 #s
reactor_network.advance(time)
There are two possibilities to advance a reactor simulation in time :
- the advance(self, double t) option (shown in the example above)
This will advance the state of the reactor network in time from the current time to the specified 't' time, taking as many integrator timesteps as necessary. - the step(self, double t) option.
This will take a single internal time step toward the specified time 't' [s]... The use of 'advance' is recommended, unless you need to elucidate a bug.
2. Simple closed vessel¶
Import statements¶
import sys
import numpy as np
import cantera as ct
import matplotlib.pyplot as plt
from matplotlib import *
Set the mechanism properties¶
gas = ct.Solution('gri30.yaml')
gas.TPX = 1000.0, ct.one_atm, 'CH4:0.5,O2:1,N2:3.76'
Create the reactor and the ReactorNet associated to it¶
# Create Reactor and fill with gas
r = ct.Reactor(gas)
# Prepare the simulation with a ReactorNet object
sim = ct.ReactorNet([r])
time = 4.e-1
Compute the reactor¶
# Arrays to hold the datas
n_points = 400
times = np.zeros(n_points)
data = np.zeros((n_points, 4))
# Advance the simulation in time
# and print the internal evolution of temperature, volume and internal energy
print(('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'vol [m3]', 'u [J/kg]')))
for n in range(n_points):
time += 5.e-3
sim.advance(time)
times[n] = time # time in s
data[n, 0] = r.T # set the temperature in the first column
data[n, 1:] = r.thermo['O2', 'CO2', 'CH4'].X # set the molar fractions in the other column
print(('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T,
r.thermo.v, r.thermo.u)))
Plot results¶
rcParams['figure.figsize'] = (14, 10)
plt.clf()
plt.subplot(2, 2, 1)
plt.plot(times, data[:, 0])
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2, 2, 2)
plt.plot(times, data[:, 1])
plt.xlabel('Time (ms)')
plt.ylabel('O2 Mole Fraction')
plt.subplot(2, 2, 3)
plt.plot(times, data[:, 2])
plt.xlabel('Time (ms)')
plt.ylabel('CO2 Mole Fraction')
plt.subplot(2, 2, 4)
plt.plot(times, data[:, 3])
plt.xlabel('Time (ms)')
plt.ylabel('CH4 Mole Fraction')
plt.show()
There are the characteristic evolutions that one can observe when simulating 0D cases. You should be aware that your case must auto-ignite (for some temperature, it is not sufficiently hot to auto-ignite, therefore nothing happens) and that the simulated time should be sufficient to capture the time where it ignites.
- at 960 K
- at 1040 K
What do you observe ?
As you can see, slightly shifting the temperature up or down moves the auto-ignition time. As you never know a priori the order of magnitude of the ignition time, it is good to use the step version until you reach burnt gases.
3. A simple constant pressure reactor¶
Here, we want to create a simple constant pressure reactor. To do so, it is necessary to create a reactor and its environment (which will be a Reservoir). The interface between the two objects created is handled by a wall, of which the expansion rate can be defined by the user.
Set the mechanism properties¶
# Mechanisms used for the process
gri3 = ct.Solution('gri30.yaml')
air = ct.Solution('air.yaml')
# Gas state
gri3.TPX = 1000.0, ct.one_atm, 'CH4:0.5,O2:1,N2:3.76'
Set the reservoir¶
- an IdealGasReactor object with the gri3 object created above
- a Reservoir representing the environment with the air object created above
# Reactor and environment
r = ct.IdealGasReactor(gri3)
env = ct.Reservoir(air)
# Wall
w = ct.Wall(r, env)
w.expansion_rate_coeff = 1.0e6 # set expansion parameter. dV/dt = KA(P_1 - P_2)
w.area = 1.0
# Prepare the simulation with a ReactorNet object
sim = ct.ReactorNet([r])
time = 4.e-1
As explained in the lecture, note that the environment defines the air in which the reactor is set in. We also define a wall between the reactor and the environment, and make it flexible, so that the pressure in the reactor is held at the environment pressure.
Simulate the reactor¶
# Arrays to hold the datas
times = np.zeros(200)
data = np.zeros((200, 4))
# Advance the simulation in time
print(('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [Pa]', 'h [J/kg]')))
for n in range(200):
time += 5.e-3
sim.advance(time)
times[n] = time # time in s
data[n, 0] = r.T
data[n, 1:] = r.thermo['O2', 'CO2', 'CH4'].X
print(('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T,
r.thermo.P, r.thermo.h)))
rcParams['figure.figsize'] = (14, 10)
plt.clf()
plt.subplot(2, 2, 1)
plt.plot(times, data[:, 0])
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2, 2, 2)
plt.plot(times, data[:, 1])
plt.xlabel('Time (ms)')
plt.ylabel('O2 Mole Fraction')
plt.subplot(2, 2, 3)
plt.plot(times, data[:, 2])
plt.xlabel('Time (ms)')
plt.ylabel('CO2 Mole Fraction')
plt.subplot(2, 2, 4)
plt.plot(times, data[:, 3])
plt.xlabel('Time (ms)')
plt.ylabel('CH4 Mole Fraction')
plt.show()
Reset gas state and set the new IdealGasConstPressureReactor¶
# Mechanisms used for the process
gri3 = ct.Solution('gri30.yaml')
# Gas state
gri3.TPX = 1000.0, ct.one_atm, 'CH4:0.5,O2:1,N2:3.76'
r = ct.IdealGasConstPressureReactor(gri3)
# Prepare the simulation with a ReactorNet object
sim = ct.ReactorNet([r])
time = 4.e-1
Simulate the new case¶
# Arrays to hold the datas
times = np.zeros(200)
data = np.zeros((200, 4))
# Advance the simulation in time
print(('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [Pa]', 'h [J/kg]')))
for n in range(200):
time += 5.e-3
sim.advance(time)
times[n] = time # time in s
data[n, 0] = r.T
data[n, 1:] = r.thermo['O2', 'CO2', 'CH4'].X
print(('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T,
r.thermo.P, r.thermo.h)))
Plot the results and compare with the previous one¶
rcParams['figure.figsize'] = (14, 10)
plt.clf()
plt.subplot(2, 2, 1)
plt.plot(times, data[:, 0])
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2, 2, 2)
plt.plot(times, data[:, 1])
plt.xlabel('Time (ms)')
plt.ylabel('O2 Mole Fraction')
plt.subplot(2, 2, 3)
plt.plot(times, data[:, 2])
plt.xlabel('Time (ms)')
plt.ylabel('CO2 Mole Fraction')
plt.subplot(2, 2, 4)
plt.plot(times, data[:, 3])
plt.xlabel('Time (ms)')
plt.ylabel('CH4 Mole Fraction')
plt.show()
4. An example of application : mixing two streams¶
Up until now, we have seen how to generate simple closed vessels. However, it is sometimes interesting to mix different streams, or to feed the exhaust of a reservoir with constant parameters to a reacting reactor. In this exercise, we will design a script that simulates the stoichiometric constant volume mixing of a stream of air with a stream of pure gaseous methane.
Set the two gases and the molar mass associated¶
- set a first gas a with an air mixture at 0.21 of O2, 0.01 of AR and 0.78 of N2
- set a second gas b with a gri30 mixture at 1 of CH4
gas_a = ct.Solution('air.yaml')
gas_a.TPX = 300, ct.one_atm, 'O2:0.21, N2:0.78, AR:0.01'
mw_a = gas_a.mean_molecular_weight/1000 #kg/mol
gas_b = ct.Solution('gri30.yaml')
gas_b.TPX = 300.0, ct.one_atm, 'CH4:1'
mw_b = gas_b.mean_molecular_weight/1000 #kg/mol
Creation of two reservoirs and an exhaust one (ReactorBase)¶
- set a first reservoir from gas a
- set a second reservoir from gas b
- set a downstream reservoir filled with gas b
res_a = ct.Reservoir(gas_a)
res_b = ct.Reservoir(gas_b)
downstream = ct.Reservoir(gas_b)
Generate the reactor that will receive the mixed stream (ReactorBase)¶
- change the condition of the gas b to respect an air mix.
- create and IdealGasReactor with the energy equation set to 'on'.
gas_b.TPX = 300, ct.one_atm, 'O2:1., N2:3.78, CH4:0.5'
mixer = ct.IdealGasReactor(gas_b, energy='on')
Watch out here, the mixture should be the one of gas_b, as the reactions will occur with the mechanism from gri30 (and not from the air).
Creation of the mass flow controller that will control the mixing (FlowDevice)¶
Now that you have specified all the elements of the reactor, you need to detail the link between them. This is what is done in the following lines :
mfca = ct.MassFlowController(res_a, mixer)
mdota = mfca.mass_flow_rate = mw_a*2./0.21
mfcb = ct.MassFlowController(res_b, mixer)
mdotb = mfcb.mass_flow_rate = mw_b*1.0
Creation of the valve between the mixer and the downstream (FlowDevice)¶
Finally, the valve is created.
outlet = ct.Valve(mixer, downstream, K=10.0)
This valve will enable the system to keep a constant pressure at the inlet of the mixer and at the outlet.
Creation of the net reactor to compute simulation¶
sim = ct.ReactorNet([mixer])
# Since the mixer is a reactor, we need to integrate in time to reach steady
# state. A few residence times should be enough.
print('{0:>14s} {1:>14s} {2:>14s} {3:>14s} {4:>14s}'.format('t [s]', 'T [K]', 'h [J/kg]', 'P [Pa]', 'X_CH4'))
t = 0.0
for n in range(100):
tres = mixer.mass/(mdota + mdotb)
t += 2*tres
sim.advance(t)
print('{0:14.5g} {1:14.5g} {2:14.5g} {3:14.5g} {4:14.5g}'.format(t, mixer.T, mixer.thermo.h, mixer.thermo.P, mixer.thermo['CH4'].X[0]))
print(mixer.thermo.report())
Notice that the gas did not ignite, and that the composition is indeed stoichiometric. We could trigger
the ignition, by assuming that the reactor is 'already lit'. To do so, you need to initialize the content
of the reactor with hot gases. The simplest way to do so is by filling it with stoichiometric mixture of
air and methane at equilibrium. This is easily done by replacing the lines :
gas_b.TPX = 300.0, ct.one_atm, 'O2:0.21, N2:0.78, AR:0.01'
mixer = ct.IdealGasReactor(gas_b)
with :
gas_b.TPX = 300.0,ct.one_atm,'O2:1., N2:3.78, CH4:0.5'
gas_b.equilibrate("HP")
mixer = ct.IdealGasReactor(gas_b)
5. Autoignition timing¶
An interesting feature of 0D simulations is that it allows to compute the temporal evolution of a
mixture under specific conditions towards its equilibrium state. That evolution can be viewed as
a sort of autoignition of the mixture.
To go further, that last exercise will guide you through
such a computation of the autoignition timing of a methane/air mixture. Several definitions of the
autoignition timing can be found in the literature. The most commonly accepted definition relies
on the temperature gradient : the time of the sharpest temperature increase is spotted as being the
autoignition point. In order to catch this time with precision, the time step of the simulation should
be as small as possible.
Create the gas object¶
gas = ct.Solution('gri30.yaml')
gas.TPX = 1250, ct.one_atm, 'CH4:0.5, O2:1, N2:3.76'
Initialisation of the different values¶
# Initial temperatures
Temperature_range = list(range(800, 1700, 100))
# Specify the number of time steps and the time step size
nt = 100000
dt = 1.e-4 # s
# Storing auto ignitions
auto_ignitions = []
for index, Temperature in enumerate(Temperature_range):
#################################################################
# Initial temperature, Pressure and stoichiometry
gas.TPX = Temperature, ct.one_atm, 'CH4:0.5, O2:1, N2:3.76'
# Create the batch reactor
r = ct.IdealGasReactor(gas)
# Now create a reactor network consisting of the single batch reactor
sim = ct.ReactorNet([r])
# Storage space
mfrac = []
# ...
time = []
temperature = []
HR = []
# Run the simulation
# Initial simulation time
current_time = 0.0
# Loop for nt time steps of dt seconds.
for n in range(nt):
current_time += dt
sim.advance(current_time)
time.append(current_time)
temperature.append(r.T)
mfrac.append(r.thermo.Y)
HR.append(- np.dot(gas.net_production_rates, gas.partial_molar_enthalpies))
#################################################################
# Catch the autoignition timing
#################################################################
# Get the ignition delay time by the maximum value of the Heat Release rate
auto_ignition = time[HR.index(max(HR))]
print('For T = ' + str(Temperature) + ', Autoignition time = ' + str(auto_ignition) + ' s')
# Posterity
FinalTemp = temperature[nt - 1]
auto_ignitions.append(auto_ignition)
# #################################################################
# # Save results
# #################################################################
# # write output CSV file for importing into Excel
# csv_file = '3-Output/Phi-1_P-1_T-' + str(Temperature) + '_UV.csv'
# with open(csv_file, 'w') as outfile:
# writer = csv.writer(outfile)
# writer.writerow(['Auto ignition time [s]', 'Final Temperature [K]'] + gas.species_names)
# writer.writerow([auto_ignition, FinalTemp] + list(mfrac[:]))
# print('output written to ' + csv_file)
T_invert = [1000 / Temperature for Temperature in Temperature_range]
#################################################################
# Plot results
#################################################################
# create plot
plt.plot(Temperature_range, auto_ignitions, 'b-o')
plt.xlabel(r'Temperature [K]', fontsize=20)
plt.ylabel("Auto ignition [s]", fontsize=20)
plt.yscale('log')
plt.title(r'Autoignition of $CH_{4}$ + Air mixture at $\Phi$ = 1, and P = 1 bar',
fontsize=22, horizontalalignment='center')
plt.axis()
plt.grid()
plt.show()
This is a typical curve of the auto-ignition time as a function of the initial temperature, meaning that when the temperature increases, this time is reduced.