custom_reactions.py
(Source)
"""
An example demonstrating how to use custom reaction objects.
For benchmark purposes, an ignition test is run to compare simulation times.
Requires: cantera >= 2.6.0
Keywords: kinetics, benchmarking, user-defined model
"""
from timeit import default_timer
import numpy as np
from math import exp
import warnings
import cantera as ct
mech = 'gri30.yaml'
fuel = 'CH4'
gas0 = ct.Solution(mech)
species = gas0.species()
reactions = gas0.reactions()
# construct custom reactions: replace 2nd reaction with equivalent custom reaction
custom_reactions = [r for r in reactions]
custom_reactions[2] = ct.CustomReaction(
equation='H2 + O <=> H + OH',
rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T),
kinetics=gas0)
gas1 = ct.Solution(thermo='ideal-gas', kinetics='gas',
species=species, reactions=custom_reactions)
# old framework - use xml input (to be removed after Cantera 2.6)
with warnings.catch_warnings():
# suppress warning: XML input is deprecated
warnings.simplefilter("ignore")
gas2 = ct.Solution(mech.replace(".yaml", ".xml"))
# construct test case - simulate ignition
def ignition(gas):
# set up reactor
gas.TP = 1000., 5 * ct.one_atm
gas.set_equivalence_ratio(0.8, fuel, 'O2:1.0, N2:3.773')
r = ct.IdealGasReactor(gas)
net = ct.ReactorNet([r])
net.rtol_sensitivity = 2.e-5
# time reactor integration
t1 = default_timer()
net.advance(.5)
t2 = default_timer()
return 1000 * (t2 - t1)
# output results
repeat = 100
print("Average time of {} simulation runs for '{}' "
"({})".format(repeat, mech, fuel))
sim0 = 0
for i in range(repeat):
sim0 += ignition(gas0)
sim0 /= repeat
print('- New framework (YAML): '
'{0:.2f} ms (T_final={1:.2f})'.format(sim0, gas0.T))
sim1 = 0
for i in range(repeat):
sim1 += ignition(gas1)
sim1 /= repeat
print('- One Python reaction: '
'{0:.2f} ms (T_final={1:.2f}) ... '
'{2:+.2f}%'.format(sim1, gas1.T, 100 * sim1 / sim0 - 100))
sim2 = 0
for i in range(repeat):
sim2 += ignition(gas2)
sim2 /= repeat
print('- Old framework (XML): '
'{0:.2f} ms (T_final={1:.2f}) ... '
'{2:+.2f}%'.format(sim2, gas2.T, 100 * sim2 / sim0 - 100))