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

%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__}")
Using Cantera version: 2.6.0a4

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
    dot_file = "reaction_paths.dot"
    png_file = "reaction_paths.png"
    diagram.write_dot(dot_file)
    subprocess.run(f"dot {dot_file} -Tpng -o{png_file} -Gdpi=100".split())
    img = Image(filename=png_file)
    display(img)
interactive(children=(IntSlider(value=100, description='plot_step', max=1510, step=10), FloatSlider(value=0.1,…
  • Find reactions containing the specie of interest
    • C2H6 in 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()
interactive(children=(FloatSlider(value=0.01, description='annotation_cutoff', max=4.0, min=0.01), Output()), …

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()
interactive(children=(FloatLogSlider(value=1e-05, description='annotation_cutoff', max=-4.0, min=-5.0), Output…