—
Interactive Reaction Path Diagrams for CH₄ plasma chemistry.#
This example uses ipywidgets to create interactive displays of reaction path
diagrams from Cantera simulations.
Requires: cantera >= 3.0.0, matplotlib >= 2.0, ipywidgets, graphviz, scipy
Import the required libraries.#
import cantera as ct
import graphviz
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from scipy import integrate
from rizer.misc.utils import get_path_to_data
sns.set_theme("poster")
print(f"Using Cantera version: {ct.__version__}")
# Determine if we're running in a Jupyter Notebook. If so, we can enable the interactive
# diagrams. Otherwise, just draw output for a single set of inputs.
try:
from IPython.core.getipython import get_ipython
config = get_ipython()
if config is None:
raise ImportError("console")
except (ImportError, AttributeError):
is_interactive = False
else:
is_interactive = True
if is_interactive:
from IPython.display import display
from matplotlib_inline.backend_inline import set_matplotlib_formats
set_matplotlib_formats("pdf", "svg")
from ipywidgets import interact, widgets
print(
"Running in interactive mode"
if is_interactive
else "Running in non-interactive mode"
)
Using Cantera version: 3.2.0
Running in non-interactive mode
Create the gas object and set the initial conditions.#
# Load the mechanism for the Fincke GRC model of CH₄ plasma chemistry.
gas = ct.Solution(get_path_to_data("mechanisms", "Fincke_GRC.yaml"))
# Set temperature, pressure, and composition to 2000 K, 1 atm, and 1 mole of CH₄.
gas.TPX = 2000.0, 1.0 * ct.one_atm, "CH4:1.0"
# Set the residence time for the reactor network.
residence_time = 1.0 # s
# Create a batch reactor object and set solver tolerances
reactor = ct.IdealGasConstPressureReactor(gas, energy="on")
reactor_network = ct.ReactorNet([reactor])
reactor_network.atol = 1e-12
reactor_network.rtol = 1e-12
/home/runner/work/rizer/rizer/examples/cantera/plot_reaction_path_diagram.py:70: DeprecationWarning: ReactorBase.__init__: After Cantera 3.2, the default value of the `clone` argument will be `True`, resulting in an independent copy of the `phase` being created for use by this reactor. Add the `clone=False` argument to retain the old behavior of sharing `Solution` objects.
reactor = ct.IdealGasConstPressureReactor(gas, energy="on")
Store time, pressure, temperature and mole fractions.#
times: list[float] = []
pressures: list[float] = []
temperatures: list[float] = []
mole_fractions: list[np.ndarray] = []
time = 0.0
steps = 0
while time < residence_time:
times.append(time)
pressures.append(gas.P)
temperatures.append(gas.T)
mole_fractions.append(gas.X)
time = reactor_network.step()
steps += 1
Plot temperature evolution with time.#
fig, ax = plt.subplots(1, 1, figsize=(20, 10))
ax.plot(times, temperatures)
ax.set_ylabel("Temperature [K]")
ax.set_title("Temperature evolution")
ax.set_xlabel("Time [s]")
ax.set_xscale("log")
ax.set_xlim(1e-5, 1.0)
(1e-05, 1.0)
Plot species evolution with time.#
fig, ax = plt.subplots(constrained_layout=True)
latex_species = {
"CH4": r"$\mathrm{CH_4}$",
"H2": r"$\mathrm{H_2}$",
"C2H2": r"$\mathrm{C_2H_2}$",
"C2H4": r"$\mathrm{C_2H_4}$",
"C2H6": r"$\mathrm{C_2H_6}$",
"C(s)": r"$\mathrm{C_{(s)}}$",
}
for species in ["CH4", "H2", "C2H2", "C2H4", "C2H6", "C(s)"]:
mole_fraction_index = gas.species_index(species)
x = [x[mole_fraction_index] * 100 for x in mole_fractions]
ax.plot(times, x)
idx_max = np.argmax(x)
time_idx_max = times[idx_max] if idx_max > 0 else 1e-5
ax.text(
x=time_idx_max,
y=x[idx_max] + 5,
s=latex_species[species],
color=ax.lines[-1].get_color(),
horizontalalignment="center",
verticalalignment="center",
)
# ax.legend()
ax.set_xlabel("Time [s]")
ax.set_ylabel("Mole fraction [%]")
ax.set_title("Mole fraction evolution with time.")
ax.set_xscale("log")
ax.set_xlim(1e-5, 1.0)
(1e-05, 1.0)
Interactive reaction path diagram.#
When executed as a Jupyter Notebook, the plotted time step, threshold and element can be changed using the slider provided by IPyWidgets.
def plot_reaction_path_diagrams(plot_step, threshold, details, element):
P = pressures[plot_step]
T = temperatures[plot_step]
X = mole_fractions[plot_step]
time = times[plot_step]
gas.TPX = T, P, X
diagram = ct.ReactionPathDiagram(gas, element)
diagram.threshold = threshold
diagram.title = f"time = {time:.2g} s"
diagram.show_details = details
if is_interactive:
graph = graphviz.Source(diagram.get_dot())
display(graph)
else:
graph = graphviz.Source(diagram.get_dot(), format="svg")
return graph
if is_interactive:
interact(
plot_reaction_path_diagrams,
plot_step=widgets.IntSlider(value=100, min=0, max=steps - 1, step=10),
threshold=widgets.FloatSlider(value=0.1, min=0.001, max=0.4, step=0.01),
details=widgets.ToggleButton(),
element=widgets.Dropdown(
options=gas.element_names,
value="C",
description="Element",
disabled=False,
),
)
diagram = ""
else:
# For non-interactive use, just draw the diagram for a specified time step
diagram = plot_reaction_path_diagrams(
plot_step=520, threshold=0.1, details=False, element="C"
)
class PlotGraphviz:
# See https://stackoverflow.com/questions/65008861/capturing-graphviz-figures-in-sphinx-gallery.
def __init__(self, dot_string):
self.dot_string = dot_string
def _repr_html_(self):
return graphviz.Source(self.dot_string)._repr_image_svg_xml()
PlotGraphviz(str(diagram))
<__main__.PlotGraphviz object at 0x7f54d9373230>
Interactive plot of instantaneous fluxes.#
Find reactions containing the species of interest, C₂H₆ in this case.
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: list[int] = C2H6_stoichiometry.nonzero()[0].tolist()
Net rates of progress of reactions containing interested species.#
The following cell calculates net rates of progress of reactions containing the species of interest and stores them.
C2H6_production_rates_list: list[np.ndarray] = []
for i in range(len(times)):
X = mole_fractions[i]
t = times[i]
T = temperatures[i]
P = pressures[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)
C2H6_production_rates_list.append(C2H6_production_rates[C2H6_reaction_indices])
# Create the instantaneous flux plot. When executed as a Jupyter Notebook, the threshold
# for annotating of reaction strings can be changed using the slider provided by IPyWidgets.
plt.rcParams["figure.constrained_layout.use"] = True
def plot_instantaneous_fluxes(annotation_cutoff):
_ = plt.figure(figsize=(6, 6))
plt.plot(times, np.array(C2H6_production_rates_list))
for i, C2H6_production_rate in enumerate(np.array(C2H6_production_rates_list).T):
peak_index = abs(C2H6_production_rate).argmax()
peak_time = times[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("<=>", "="),
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)")
plt.ylabel("Net rates of C2H6 production")
plt.xscale("log")
plt.xlim(1e-7, 1.0)
plt.show()
if is_interactive:
interact(
plot_instantaneous_fluxes,
annotation_cutoff=widgets.FloatSlider(value=0.1, min=1e-2, max=4, steps=10),
)
else:
plot_instantaneous_fluxes(annotation_cutoff=0.1)
Interactive plot of integrated fluxes.#
When executed as a Jupyter Notebook, the threshold for annotating of reaction strings can be changed using the slider provided by iPyWidgets
# Integrate fluxes over time
integrated_fluxes = integrate.cumulative_trapezoid(
np.array(C2H6_production_rates_list),
np.array(times),
axis=0,
initial=0,
)
def plot_integrated_fluxes(integrated_fluxes, annotation_cutoff):
_ = plt.figure(figsize=(6, 6))
plt.plot(times, integrated_fluxes)
final_time = times[-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, fontsize=8
)
plt.xlabel("Time (s)")
plt.ylabel("Integrated net rate of progress")
plt.title("Cumulative C2H6 formation")
plt.xscale("log")
plt.xlim(1e-5, 1.0)
plt.show()
if is_interactive:
interact(
plot_integrated_fluxes,
annotation_cutoff=widgets.FloatLogSlider(
value=1e-5, min=-5, max=-4, base=10, step=0.1
),
integrated_fluxes=widgets.fixed(integrated_fluxes),
)
else:
plot_integrated_fluxes(integrated_fluxes=integrated_fluxes, annotation_cutoff=1e-5)
Total running time of the script: (0 minutes 0.933 seconds)