-
-
Notifications
You must be signed in to change notification settings - Fork 342
/
mechanism_reduction.py
97 lines (78 loc) · 3.17 KB
/
mechanism_reduction.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
"""
A simplistic approach to mechanism reduction which demonstrates Cantera's
features for dynamically manipulating chemical mechanisms.
Here, we use the full GRI 3.0 mechanism to simulate adiabatic, constant pressure
ignition of a lean methane/air mixture. We track the maximum reaction rates for
each reaction to determine which reactions are the most important, according to
a simple metric based on the relative net reaction rate.
We then create a sequence of reduced mechanisms including only the top reactions
and the associated species, and run the simulations again with these mechanisms
to see whether the reduced mechanisms with a certain number of species are able
to adequately simulate the ignition delay problem.
Requires: cantera >= 2.6.0, matplotlib >= 2.0
Keywords: kinetics, combustion, reactor network, editing mechanisms, ignition delay,
plotting
"""
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
gas = ct.Solution('gri30.yaml')
initial_state = 1200, 5 * ct.one_atm, 'CH4:0.35, O2:1.0, N2:3.76'
# Run a simulation with the full mechanism
gas.TPX = initial_state
r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r])
tt = []
TT = []
t = 0.0
# Rmax is the maximum relative reaction rate at any timestep
Rmax = np.zeros(gas.n_reactions)
while t < 0.02:
t = sim.step()
tt.append(1000 * t)
TT.append(r.T)
rnet = abs(gas.net_rates_of_progress)
rnet /= max(rnet)
Rmax = np.maximum(Rmax, rnet)
plt.plot(tt, TT, label='K=53, R=325', color='k', lw=3, zorder=100)
# Get the reaction objects, and sort them so the most active reactions are first
R = sorted(zip(Rmax, gas.reactions()), key=lambda x: -x[0])
# Test reduced mechanisms with different numbers of reactions
C = plt.cm.winter(np.linspace(0, 1, 5))
for i, N in enumerate([40, 50, 60, 70, 80]):
# Get the N most active reactions
reactions = [r[1] for r in R[:N]]
# find the species involved in these reactions. At a minimum, include all
# species in the reactant mixture
species_names = {'N2', 'CH4', 'O2'}
for reaction in reactions:
species_names.update(reaction.reactants)
species_names.update(reaction.products)
# Get the species objects
species = [gas.species(name) for name in species_names]
# create the new reduced mechanism
gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
species=species, reactions=reactions)
# save the reduced mechanism for later use
gas2.write_yaml("gri30-reduced-{}-reaction.yaml".format(N))
# Re-run the ignition problem with the reduced mechanism
gas2.TPX = initial_state
r = ct.IdealGasConstPressureReactor(gas2)
sim = ct.ReactorNet([r])
t = 0.0
tt = []
TT = []
while t < 0.02:
t = sim.step()
tt.append(1000 * t)
TT.append(r.T)
plt.plot(tt, TT, lw=2, color=C[i],
label='K={0}, R={1}'.format(gas2.n_species, N))
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.legend(loc='upper left')
plt.title('Reduced mechanism ignition delay times\n'
'K: number of species; R: number of reactions')
plt.xlim(0, 20)
plt.tight_layout()
plt.show()