interactive_path_diagram.ipynb (Source)
Interactive Path Diagrams¶
This example uses ipywidgets to create interactive displays of path diagrams from Cantera simulations.
In [1]:
from IPython.display import Image, display
from ipywidgets import widgets, interact
import cantera as ct
import numpy as np
import graphviz
%matplotlib inline
%config InlineBackend.figure_formats = ["svg"]
from matplotlib import pyplot as plt
from collections import defaultdict
import subprocess
plt.rcParams["figure.dpi"] = 120
print(f"Using Cantera version: {ct.__version__}")
When using Cantera, the first thing you usually need is an object representing some phase of matter. Here, we'll create a gas mixture using GRI-Mech:
In [2]:
gas = ct.Solution("gri30.yaml")
Shock tube ignition delay measurement conditions (Spadaccini, L.J., and Colket, M.B., (1994) Prog. Energy Combust. Sci. 20, 431. )¶
- $CH_4-C_2H_6-O_2-Ar$ (3.29%-0.21%-7%-89.5%)
- phi = 1.045
- P = 6.1 - 7.6 atm
- T = 1356 - 1688 K
- Set temperature,pressure and composition
In [3]:
gas.TPX = 1550.0, 6.5 * ct.one_atm, "CH4:3.29, C2H6:0.21, O2:7 , Ar:89.5"
- Residence time is close to ignition delay reported by (Spadaccini, L.J., and Colket, M.B., (1994) Prog. Energy Combust. Sci. 20, 431. )
In [4]:
residence_time = 1e-3
- Create a batch reactor object and set solver tolerances
In [5]:
reactor = ct.IdealGasConstPressureReactor(gas, energy="on")
reactor_network = ct.ReactorNet([reactor])
reactor_network.atol = 1e-12
reactor_network.rtol = 1e-12
- Store time, pressure, temperature and mole fractions
In [6]:
profiles = defaultdict(list)
time = 0
steps = 0
while time < residence_time:
profiles["time"].append(time)
profiles["pressure"].append(gas.P)
profiles["temperature"].append(gas.T)
profiles["mole_fractions"].append(gas.X)
time = reactor_network.step()
steps += 1
Interactive reaction path diagram¶
- Plot steps, threshold and element can be changed using the slider provided by IPyWidgets
In [7]:
@interact(
plot_step=widgets.IntSlider(value=100, min=0, max=steps, step=10),
threshold=widgets.FloatSlider(value=0.1, min=0, max=1, step=0.01),
details=widgets.ToggleButton(),
species=widgets.Dropdown(
options=gas.element_names,
value="C",
description="Element",
disabled=False,
),
)
def plot_reaction_path_diagrams(plot_step, threshold, details, species):
P = profiles["pressure"][plot_step]
T = profiles["temperature"][plot_step]
X = profiles["mole_fractions"][plot_step]
time = profiles["time"][plot_step]
gas.TPX = T, P, X
print("time = {:.2g} s".format(time))
diagram = ct.ReactionPathDiagram(gas, species)
diagram.threshold = threshold
diagram.show_details = details
graph = graphviz.Source(diagram.get_dot())
display(graph)
- Find reactions containing the specie of interest
-
C2H6in this case
-
In [8]:
C2H6_stoichiometry = np.zeros_like(gas.reactions())
for i, r in enumerate(gas.reactions()):
C2H6_moles = r.products.get("C2H6", 0) - r.reactants.get("C2H6", 0)
C2H6_stoichiometry[i] = C2H6_moles
C2H6_reaction_indices = C2H6_stoichiometry.nonzero()[0]
The following cell calculates net rates of progress of reactions containing the species of interest and stores it¶
In [9]:
profiles["C2H6_production_rates"] = []
for i in range(len(profiles["time"])):
X = profiles["mole_fractions"][i]
t = profiles["time"][i]
T = profiles["temperature"][i]
P = profiles["pressure"][i]
gas.TPX = (T, P, X)
C2H6_production_rates = (
gas.net_rates_of_progress
* C2H6_stoichiometry # [kmol/m^3/s]
* gas.volume_mass # Specific volume [m^3/kg].
) # overall, mol/s/g (g total in reactor, same basis as N_atoms_in_fuel)
profiles["C2H6_production_rates"].append(
C2H6_production_rates[C2H6_reaction_indices]
)
Interactive plot of instantaneous fluxes¶
- Threshold for annotating of reaction strings can be changed using the slider provided by IPyWidgets
In [10]:
@interact(
annotation_cutoff=widgets.FloatSlider(value=1e-2, min=1e-2, max=4, steps=10),
profiles=widgets.fixed(profiles),
)
def plot_instantaneous_fluxes(profiles, annotation_cutoff):
profiles = profiles
fig = plt.figure(figsize=(6, 6))
plt.plot(profiles["time"], np.array(profiles["C2H6_production_rates"]))
for i, C2H6_production_rate in enumerate(
np.array(profiles["C2H6_production_rates"]).T
):
peak_index = abs(C2H6_production_rate).argmax()
peak_time = profiles["time"][peak_index]
peak_C2H6_production = C2H6_production_rate[peak_index]
reaction_string = gas.reaction_equations(C2H6_reaction_indices)[i]
if abs(peak_C2H6_production) > annotation_cutoff:
plt.annotate(
(reaction_string).replace("2", "$_2$").replace("<=>", "="),
xy=(peak_time, peak_C2H6_production),
xytext=(
peak_time * 2,
(
peak_C2H6_production
+ 0.003
* (peak_C2H6_production / abs(peak_C2H6_production))
* (abs(peak_C2H6_production) > 0.005)
* (peak_C2H6_production < 0.06)
),
),
arrowprops=dict(
arrowstyle="->",
color="black",
relpos=(0, 0.6),
linewidth=2,
),
horizontalalignment="left",
)
plt.xlabel("Time (s)", fontsize=16)
plt.ylabel("Net rates of C2H6 production", fontsize=16)
plt.tight_layout()
plt.show()
Integrating over time using scipy.integrate.cumptraz¶
In [11]:
from scipy import integrate
integrated_fluxes = integrate.cumtrapz(
np.array(profiles["C2H6_production_rates"]),
np.array(profiles["time"]),
axis=0,
initial=0,
)
Interactive plot of integrated fluxes¶
- Threshold for annotating of reaction strings can be changed using the slider provided by iPyWidgets
In [12]:
@interact(
annotation_cutoff=widgets.FloatLogSlider(
value=1e-5, min=-5, max=-4, base=10, step=0.1
),
profiles=widgets.fixed(profiles),
integrated_fluxes=widgets.fixed(integrated_fluxes),
)
def plot_integrated_fluxes(profiles, integrated_fluxes, annotation_cutoff):
profiles = profiles
integrated_fluxes = integrated_fluxes
fig = plt.figure(figsize=(6, 6))
plt.plot(profiles["time"], integrated_fluxes)
final_time = profiles["time"][-1]
for i, C2H6_production in enumerate(integrated_fluxes.T):
total_C2H6_production = C2H6_production[-1]
reaction_string = gas.reaction_equations(C2H6_reaction_indices)[i]
if abs(total_C2H6_production) > annotation_cutoff:
plt.text(final_time * 1.06, total_C2H6_production, reaction_string)
plt.xlabel("Time (s)", fontsize=16)
plt.ylabel("Integrated net rate of progress", fontsize=16)
plt.title("Cumulative C2H6 formation", fontsize=16)
plt.tight_layout()
plt.show()