From 95f6f3ae707d62a40c1d1c69af94b1f96f5bf955 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 12 Mar 2024 18:56:00 -0400 Subject: [PATCH 01/18] Migrate nonideal equation of state example from Jupyter Co-authored-by: Gandhali Kogekar Co-authored-by: Bryan W. Weber --- doc/sphinx/conf.py | 1 + samples/python/thermo/equations_of_state.py | 280 ++++++++++++++++++++ 2 files changed, 281 insertions(+) create mode 100644 samples/python/thermo/equations_of_state.py diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py index 79f196c6ba..d595c8bdb4 100644 --- a/doc/sphinx/conf.py +++ b/doc/sphinx/conf.py @@ -52,6 +52,7 @@ "filetype_parsers": {'.h': 'C++', '.m': 'Matlab'}, 'ignore_pattern': r'(__.*__\.py|test_examples\.m)', 'image_srcset': ["2x"], + 'remove_config_comments': True, 'examples_dirs': [ '../samples/python/', '../samples/cxx/', diff --git a/samples/python/thermo/equations_of_state.py b/samples/python/thermo/equations_of_state.py new file mode 100644 index 0000000000..3790be7d34 --- /dev/null +++ b/samples/python/thermo/equations_of_state.py @@ -0,0 +1,280 @@ +""" +Non-ideal equations of state +============================ + +This example demonstrates a comparison between ideal and non-ideal equations of +state (EoS) using Cantera and CoolProp. + +The following equations of state are used to evaluate thermodynamic properties +in this example: + +1. Ideal-gas EoS from Cantera +2. Non-ideal Redlich-Kwong EoS (R-K EoS) from Cantera +3. Helmholtz energy EoS (HEOS) from CoolProp + +Requires: cantera >= 3.0.0, matplotlib >= 2.0, coolprop >= 6.0 + +.. tags:: Python, thermodynamics, non-ideal fluid +""" + +#%% +# Import required packages (Cantera and CoolProp) +# ----------------------------------------------- +# +# `CoolProp `__ [1]_ is an open-source package that contains a +# highly-accurate database for thermophysical properties. The thermodynamic properties +# are obtained using pure and pseudo-pure fluid equations of state implemented for 122 +# components. +# +# If you are using ``conda`` to manage your software environment, activate your cantera +# environment and install CoolProp by running ``conda install -c conda-forge coolprop``. +# + +import cantera as ct +import numpy as np +from CoolProp.CoolProp import PropsSI +import matplotlib.pyplot as plt + +plt.rcParams['figure.constrained_layout.use'] = True + +print(f"Running Cantera version: {ct.__version__}") + +#%% +# Helper functions +# ---------------- +# +# This example uses CO\ :math:`_2` as the only species. The function +# ``get_thermo_Cantera`` calculates thermodynamic properties based on the thermodynamic +# state (:math:`T`, :math:`p`) of the species using Cantera. Applicable phases are +# ``Ideal-gas`` and ``Redlich-Kwong``. The ideal-gas equation can be stated as +# +# .. math:: +# pv = RT, +# +# where :math:`p`, :math:`v` and :math:`T` represent thermodynamic pressure, molar +# volume, and the temperature of the gas-phase. :math:`R` is the universal gas constant. +# The Redlich-Kwong equation is a cubic, non-ideal equation of state, represented as +# +# .. math:: +# p=\frac{RT}{v-b^\ast}-\frac{a^\ast}{v\sqrt{T}(v+b^\ast)}. +# +# In this expression, :math:`R` is the universal gas constant and :math:`v` is the molar +# volume. The temperature-dependent van der Waals attraction parameter :math:`a^\ast` +# and volume correction parameter (repulsive parameter) :math:`b^\ast` represent +# molecular interactions. +# +# The function ``get_thermo_CoolProp`` utilizes the CoolProp package to evaluate +# thermodynamic properties based on the thermodynamic state (:math:`T`, :math:`p`) for a +# given fluid. The HEOS for CO\ :math:`_2` used in this example is obtained from +# http://www.coolprop.org/fluid_properties/fluids/CarbonDioxide.html. +# +# Since the standard-reference thermodynamic states are different for Cantera and +# CoolProp, it is necessary to convert these values to an appropriate scale before +# comparison. Therefore, both functions ``get_thermo_Cantera`` and +# ``get_thermo_CoolProp`` return the thermodynamic values relative to a reference state +# at 1 bar, 300 K. +# +# To plot the comparison of thermodynamic properties among the three EoS, the ``plot`` +# function is used. + +def get_thermo_Cantera(phase, T, p): + states = ct.SolutionArray(phase, len(p)) + X = "CO2:1.0" + states.TPX = T, p, X + + u = states.u / 1000 + h = states.h / 1000 + s = states.s / 1000 + cp = states.cp / 1000 + cv = states.cv / 1000 + + # Get the relative enthalpy, entropy and int. energy with reference to the first + # point + u = u - u[0] + s = s - s[0] + h = h - h[0] + + return h, u, s, cp, cv + + +def get_thermo_CoolProp(T, p): + u = np.zeros_like(p) + h = np.zeros_like(p) + s = np.zeros_like(p) + cp = np.zeros_like(p) + cv = np.zeros_like(p) + + for i in range(p.shape[0]): + u[i] = PropsSI("U", "P", p[i], "T", T, "CO2") + h[i] = PropsSI("H", "P", p[i], "T", T, "CO2") + s[i] = PropsSI("S", "P", p[i], "T", T, "CO2") + cp[i] = PropsSI("C", "P", p[i], "T", T, "CO2") + cv[i] = PropsSI("O", "P", p[i], "T", T, "CO2") + + # Get the relative enthalpy, entropy and int. energy with reference to the first + # point + u = (u - u[0]) / 1000 + s = (s - s[0]) / 1000 + h = (h - h[0]) / 1000 + cp = cp / 1000 + cv = cv / 1000 + + return h, u, s, cp, cv + + +def plot(p, thermo_Ideal, thermo_RK, thermo_CoolProp, name): + line_width = 3 + + fig, ax = plt.subplots() + ax.plot( + p / 1e5, thermo_Ideal, "-", color="b", linewidth=line_width, label="Ideal EoS" + ) + ax.plot(p / 1e5, thermo_RK, "-", color="r", linewidth=line_width, label="R-K EoS") + ax.plot( + p / 1e5, thermo_CoolProp, "-", color="k", linewidth=line_width, label="CoolProp" + ) + ax.set_xlabel("Pressure [bar]") + ax.set_ylabel(name) + ax.legend(prop={"size": 14}, frameon=False) + +#%% +# EoS Comparison based on thermodynamic properties +# ------------------------------------------------ +# +# This is the main subroutine that compares and plots the thermodynamic values obtained +# using three equations of state. + +# Input parameters +T = 300 # Temperature is constant [unit: K] +p = 1e5 * np.linspace(1, 100, 1000) # Pressure is varied from 1 to 100 bar [unit: Pa] + +# Read the ideal-gas phase +ideal_gas_phase = ct.Solution("example_data/co2-thermo.yaml", "CO2-Ideal") +h_ideal, u_ideal, s_ideal, cp_ideal, cv_ideal = get_thermo_Cantera( + ideal_gas_phase, T, p +) + +# Read the Redlich-Kwong phase +redlich_kwong_phase = ct.Solution("example_data/co2-thermo.yaml", "CO2-RK") +h_RK, u_RK, s_RK, cp_RK, cv_RK = get_thermo_Cantera(redlich_kwong_phase, T, p) + +# Read the thermo data using CoolProp +h_CoolProp, u_CoolProp, s_CoolProp, cp_CoolProp, cv_CoolProp = get_thermo_CoolProp(T, p) + +# Plot the results +plot(p, u_ideal, u_RK, u_CoolProp, "Relative Internal Energy [kJ/kg]") +#%% +plot(p, h_ideal, h_RK, h_CoolProp, "Relative Enthalpy [kJ/kg]") +#%% +plot(p, s_ideal, s_RK, s_CoolProp, "Relative Entropy [kJ/kg-K]") + +#%% +# The thermodynamic properties such as internal energy, enthalpy, and entropy are +# plotted against the operating pressure at a constant temperature :math:`T = 300` K. +# The three equations follow each other closely at low pressures (:math:`P < 10` bar). +# However, the ideal gas EoS departs significantly from the observed behavior of gases +# near the critical regime (:math:`P_{\rm {crit}} = 73.77` bar). +# +# The ideal gas EoS does not consider inter-molecular interactions and the volume +# occupied by individual gas particles. At low temperatures and high pressures, +# inter-molecular forces become particularly significant due to a reduction in +# inter-molecular distances. Additionally, at high density, the volume of individual +# molecules becomes significant. Both of these factors contribute to the deviation from +# ideal behavior at high pressures. The cubic Redlich-Kwong EoS, on the other hand, +# predicts thermodynamic properties accurately near the critical regime. + +# Specific heat at constant pressure +plot(p, cp_ideal, cp_RK, cp_CoolProp, "$C_p$ [kJ/kg-K]") +#%% + +# Specific heat at constant volume +plot(p, cv_ideal, cv_RK, cv_CoolProp, "$C_v$ [kJ/kg-K]") + +#%% +# In the case of Ideal gas EoS, the specific heats at constant pressure (:math:`C_{\rm +# p}`) and constant volume (:math:`C_{\rm v}`) are independent of the pressure. Hence, +# :math:`C_{\rm p}` and :math:`C_{\rm v}` for ideal EoS do not change as the pressure is +# varied from :math:`1` bar to :math:`100` bar in this study. +# +# :math:`C_{\rm p}` for the R-K EoS follows the trend closely with the Helmholtz EoS +# from CoolProp up to the critical regime. Although :math:`C_{\rm p}` shows reasonable +# agreement with the Helmholtz EoS in sub-critical and supercritical regimes, it +# inaccurately predicts a very high value near the critical point. However, +# :math:`C_{\rm p}` at the critical point is finite for the real fluid. The sudden rise +# in :math:`C_{\rm p}` in the case of the R-K EoS is just a numerical artifact, due to +# the EoS yielding infinite values in the limiting case, and not a real singularity. +# +# :math:`C_{\rm v}`, on the other hand, predicts smaller values in the subcritical and +# critical regime. However, it shows completely incorrect values in the super-critical +# region, making it invalid at very high pressures. It is well known that the cubic +# equations typically fail to predict accurate constant-volume heat capacity in the +# transcritical region [2]_. Certain cubic EoS models have been extended to resolve this +# discrepancy using crossover models. For further information, see the work of Span [2]_ +# and Saeed et al. [3]_. + +#%% +# Temperature-Density plots +# ------------------------- +# +# The following function plots the :math:`T`-:math:`\rho` diagram over a wide pressure +# and temperature range. The temperature is varied from 250 K to 400 K. The pressure is +# changed from 1 bar to 600 bar. + +# Input parameters +# Set up arrays for pressure and temperature +p_array = np.logspace(1, np.log10(600), 10, endpoint=True) +T_array = np.linspace(250, 401, 20) # Temperature is varied from 250K to 400K +p_array = 1e5 * np.array(p_array)[:, np.newaxis] + +# Calculate densities for Ideal gas and R-K EoS phases +states = ct.SolutionArray(ideal_gas_phase, shape=(p_array.size, T_array.size)) +states.TP = T_array, p_array +density_ideal = states.density_mass + +states = ct.SolutionArray(redlich_kwong_phase, shape=(p_array.size, T_array.size)) +states.TP = T_array, p_array +density_RK = states.density_mass + +p, T = np.meshgrid(p_array, T_array) +density_coolprop = PropsSI("D", "P", np.ravel(p), "T", np.ravel(T), "CO2") +density_coolprop = density_coolprop.reshape(p.shape) + +# Plot +import cycler +color = plt.cm.viridis(np.linspace(0, 1, p_array.size)) +plt.rcParams['axes.prop_cycle'] = cycler.cycler('color', color) + +fig, ax = plt.subplots() +ideal_line = ax.plot(density_ideal.T, T_array, "--", label="Ideal EoS") +RK_line = ax.plot(density_RK.T, T_array, "o", label="R-K EoS") +CP_line = ax.plot(density_coolprop, T_array, "-", label="CoolProp") +ax.text(-27.5, 320, "p = 1 bar", color=color[0], rotation="vertical") +ax.text(430, 318, "p = 97 bar", color=color[5], rotation=-12) +ax.text(960, 320, "p = 600 bar", color=color[9], rotation=-68) +ax.set_xlabel("Density [kg/m$^3$]") +ax.set_ylabel("Temperature [K]") +_ = ax.legend(handles=[ideal_line[0], RK_line[0], CP_line[0]]) + +#%% +# The figure compares :math:`T-\rho` plots for ideal, R-K, and Helmholtz EoS at +# different operating pressures. All three EoS yield the same plots at low pressures (0 +# bar and 10 bar). However, the Ideal gas EoS departs significantly at high pressures +# (:math:`P > 10` bar), where non-ideal effects are prominent. The R-K EoS closely +# matches the Helmholtz EoS at supercritical pressures (:math:`P \ge 70` bar). However, +# it does depart in the liquid-vapor region that exists at :math:`P < P_{\rm {crit}}` +# and low temperatures (:math:`~T_{\rm {crit}}`). +# +# .. [1] I.H. Bell, J. Wronski, S. Quoilin, V. Lemort, "Pure and Pseudo-pure Fluid +# Thermophysical Property Evaluation and the Open-Source Thermophysical Property +# Library CoolProp," Industrial & Engineering Chemistry Research 53 (2014), +# https://pubs.acs.org/doi/10.1021/ie4033999 +# +# .. [2] R. Span, "Multiparameter Equations of State - An Accurate Source of +# Thermodynamic Property Data," Springer Berlin Heidelberg (2000), +# https://dx.doi.org/10.1007/978-3-662-04092-8 +# +# .. [3] A. Saeed, S. Ghader, "Calculation of density, vapor pressure and heat capacity +# near the critical point by incorporating cubic SRK EoS and crossover translation," +# Fluid Phase Equilibria (2019) 493, https://doi.org/10.1016/j.fluid.2019.03.027 + +# sphinx_gallery_thumbnail_number = -1 From 4b84d688e35ad6e926924a39c9998f7a37056e4d Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 15 Mar 2024 14:35:34 -0400 Subject: [PATCH 02/18] [Doc] Fix error in PFR surface species equation See https://groups.google.com/g/cantera-users/c/_-pzZdv1lDM --- doc/sphinx/reference/reactors/pfr.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/reference/reactors/pfr.md b/doc/sphinx/reference/reactors/pfr.md index 33ae427e2b..ac1066c15b 100644 --- a/doc/sphinx/reference/reactors/pfr.md +++ b/doc/sphinx/reference/reactors/pfr.md @@ -97,7 +97,7 @@ To satisfy the constraint that the total surface coverage is 1, the conservation equation for the first surface species on each surface is replaced by this constraint: $$ -\dot{s}_{0,j} = 1 - \sum_{i\ne 0} \theta_{i,j} +\sum_{i} \theta_{i,j} = 1 $$ (pfr-surf-species-0) Without this constraint, the solver could find the trivial, non-physical solution From b2cb80bd8fa691a289da10420db530a197be3b13 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 15 Mar 2024 16:59:23 -0400 Subject: [PATCH 03/18] Migrate flame temperature example from Jupyter --- doc/sphinx/userguide/flame-temperature.md | 93 +++++++++++++++++++++++ doc/sphinx/userguide/index.md | 5 ++ 2 files changed, 98 insertions(+) create mode 100644 doc/sphinx/userguide/flame-temperature.md diff --git a/doc/sphinx/userguide/flame-temperature.md b/doc/sphinx/userguide/flame-temperature.md new file mode 100644 index 0000000000..acde064088 --- /dev/null +++ b/doc/sphinx/userguide/flame-temperature.md @@ -0,0 +1,93 @@ +--- +file_format: mystnb +kernelspec: + name: python3 +--- + +```{py:currentmodule} cantera +``` + +# Calculating Adiabatic Flame Temperature + +This guide demonstrates calculation of the adiabatic flame temperature for a methane/air +mixture, comparing calculations which assume either complete or incomplete combustion. + +## Complete Combustion + +The stoichiometric equation for complete combustion of a lean methane/air mixture ($\phi +< 1$) is: + +$$ +\t{\phi CH_4 + 2(O_2 + 3.76 N_2) \rightarrow + \phi CO_2 + 2\phi H_2O + 2 (1-\phi) O_2 + 7.52 N_2} +$$ + +For a rich mixture ($\phi > 1$), this becomes: + +$$ +\t{\phi CH_4 + 2(O_2 + 3.76 N_2) \rightarrow + CO_2 + 2 H_2O + (\phi - 1) CH_4 + 7.52 N_2} +$$ + +To find the flame temperature resulting from these reactions using Cantera, we create a +`Solution` object containing only the species in the above stoichiometric equations, and +then use the `Solution.equilibrate` function to find the resulting mixture composition +and temperature, taking advantage of the fact that equilibrium will strongly favor +conversion of the fuel molecule. + +```{code-cell} python +import cantera as ct +import numpy as np + +# Get all of the Species objects defined in the GRI 3.0 mechanism +species = {S.name: S for S in ct.Species.list_from_file("gri30.yaml")} + +# Create an ideal gas phase object with species representing complete combustion +complete_species = [species[S] for S in ("CH4", "O2", "N2", "CO2", "H2O")] +gas1 = ct.Solution(thermo="ideal-gas", species=complete_species) + +phi = np.linspace(0.5, 2.0, 100) +T_complete = np.zeros(phi.shape) +for i in range(len(phi)): + gas1.TP = 300, ct.one_atm + gas1.set_equivalence_ratio(phi[i], "CH4", "O2:1, N2:3.76") + gas1.equilibrate("HP") + T_complete[i] = gas1.T +``` + +## Incomplete Combustion + +In the case of incomplete combustion, the resulting mixture composition is not known in +advance, but must be found by calculating the equilibrium composition at constant +enthalpy and temperature: + +$$ +\t{\phi CH_4 + 2(O_2 + 3.76 N_2) \rightarrow + ? CO_2 + ? CO + ? H_2 + ? H_2O + ? O_2 + 7.52 N_2 + minor\ species} +$$ + +Now, we use a gas phase object containing all 53 species defined in the GRI 3.0 +mechanism, and compute the equilibrium composition as a function of equivalence ratio. + +```{code-cell} python +# Create an IdealGas object including incomplete combustion species +gas2 = ct.Solution(thermo="ideal-gas", species=species.values()) +T_incomplete = np.zeros(phi.shape) +for i in range(len(phi)): + gas2.TP = 300, ct.one_atm + gas2.set_equivalence_ratio(phi[i], "CH4", "O2:1, N2:3.76") + gas2.equilibrate("HP") + T_incomplete[i] = gas2.T +``` + +## Comparing Results + +```{code-cell} python +import matplotlib.pyplot as plt + +plt.plot(phi, T_complete, label="complete combustion", lw=2) +plt.plot(phi, T_incomplete, label="incomplete combustion", lw=2) +plt.grid(True) +plt.xlabel(r"Equivalence ratio, $\phi$") +_ = plt.ylabel("Temperature [K]") +``` diff --git a/doc/sphinx/userguide/index.md b/doc/sphinx/userguide/index.md index 340478a5cd..e8080c36a8 100644 --- a/doc/sphinx/userguide/index.md +++ b/doc/sphinx/userguide/index.md @@ -33,6 +33,9 @@ conditions, or calculating the voltage of a Lithium-ion battery as it is dischar - [](input-errors) - [](legacy2yaml-tutorial) +### Combustion Calculations +- [](flame-temperature) + ### Implementing Custom Models - [](extensible-reactor) @@ -66,5 +69,7 @@ thermobuild input-errors legacy2yaml-tutorial +flame-temperature + extensible-reactor ``` From 49c943e636964667251b447938c61a03d231c6fe Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 15 Mar 2024 20:41:38 -0400 Subject: [PATCH 04/18] Migrate heating value example from Jupyter --- doc/sphinx/userguide/heating-value.md | 103 ++++++++++++++++++++++++++ doc/sphinx/userguide/index.md | 2 + 2 files changed, 105 insertions(+) create mode 100644 doc/sphinx/userguide/heating-value.md diff --git a/doc/sphinx/userguide/heating-value.md b/doc/sphinx/userguide/heating-value.md new file mode 100644 index 0000000000..d6f0718337 --- /dev/null +++ b/doc/sphinx/userguide/heating-value.md @@ -0,0 +1,103 @@ +--- +file_format: mystnb +kernelspec: + name: python3 +--- + +```{py:currentmodule} cantera +``` + +# Calculating Heating Value of a Fuel + +This guide demonstrates how to use Cantera's thermodynamic data to calculate the lower +[heating value](https://en.wikipedia.org/wiki/Heat_of_combustion) (LHV) and higher +heating value (HHV) of methane and other fuels. + +## Heating value of Methane + +The complete reaction for heating methane is: + +$$\t{CH_4 + 2 O_2 \rightarrow CO_2 + 2 H_2O}$$ + +We compute the lower heating value (LHV) as the difference in enthalpy (per kg +*mixture*) between reactants and products at constant temperature and pressure, divided +by the mass fraction of fuel in the reactants. + +```{code-cell} python +import cantera as ct +gas = ct.Solution("gri30.yaml") + +# Set reactants state +gas.TPX = 298, 101325, "CH4:1, O2:2" +h1 = gas.enthalpy_mass +Y_CH4 = gas["CH4"].Y[0] # returns an array, of which we only want the first element + +# set state to complete combustion products without changing T or P +gas.TPX = None, None, "CO2:1, H2O:2" +h2 = gas.enthalpy_mass + +LHV = -(h2 - h1) / Y_CH4 / 1e6 +print(f"LHV = {LHV:.3f} MJ/kg") +``` + +The LHV is calculated assuming that water remains in the gas phase. However, more energy +can be extracted from the mixture if this water is condensed. This value is the higher +heating value (HHV). + +The ideal gas mixture model used here cannot calculate this contribution directly. +However, Cantera also has a non-ideal equation of state which can be used to compute +this contribution. + +```{code-cell} python +water = ct.Water() +# Set liquid water state, with vapor fraction x = 0 +water.TQ = 298, 0 +h_liquid = water.h +# Set gaseous water state, with vapor fraction x = 1 +water.TQ = 298, 1 +h_gas = water.h + +# Calculate higher heating value +Y_H2O = gas["H2O"].Y[0] +HHV = -(h2 - h1 + (h_liquid - h_gas) * Y_H2O) / Y_CH4 / 1e6 +print(f"HHV = {HHV:.3f} MJ/kg") +``` + +## Generalizing to arbitrary species + +We can generalize this calculation by determining the composition of the products +automatically rather than directly specifying the product composition. This can be done +by computing the *elemental mole fractions* of the reactants mixture and noting that for +complete combustion, all of the carbon ends up as $\t{CO_2}$, all of the hydrogen ends +up as $\t{H_2O}$, and all of the nitrogen ends up as $\t{N_2}$. From this, we can +compute the ratio of these species in the products. + +```{code-cell} python +def heating_value(fuel): + """Returns the LHV and HHV for the specified fuel""" + gas.TP = 298, ct.one_atm + gas.set_equivalence_ratio(1.0, fuel, "O2:1.0") + h1 = gas.enthalpy_mass + Y_fuel = gas[fuel].Y[0] + + # complete combustion products + X_products = { + "CO2": gas.elemental_mole_fraction("C"), + "H2O": 0.5 * gas.elemental_mole_fraction("H"), + "N2": 0.5 * gas.elemental_mole_fraction("N"), + } + + gas.TPX = None, None, X_products + Y_H2O = gas["H2O"].Y[0] + h2 = gas.enthalpy_mass + LHV = -(h2 - h1) / Y_fuel / 1e6 + HHV = -(h2 - h1 + (h_liquid - h_gas) * Y_H2O) / Y_fuel / 1e6 + return LHV, HHV + + +fuels = ["H2", "CH4", "C2H6", "C3H8", "NH3", "CH3OH"] +print("fuel LHV (MJ/kg) HHV (MJ/kg)") +for fuel in fuels: + LHV, HHV = heating_value(fuel) + print(f"{fuel:8s} {LHV:7.3f} {HHV:7.3f}") +``` diff --git a/doc/sphinx/userguide/index.md b/doc/sphinx/userguide/index.md index e8080c36a8..acf6fc6307 100644 --- a/doc/sphinx/userguide/index.md +++ b/doc/sphinx/userguide/index.md @@ -35,6 +35,7 @@ conditions, or calculating the voltage of a Lithium-ion battery as it is dischar ### Combustion Calculations - [](flame-temperature) +- [](heating-value) ### Implementing Custom Models @@ -70,6 +71,7 @@ input-errors legacy2yaml-tutorial flame-temperature +heating-value extensible-reactor ``` From 1baf95e6e941e5f5cb37ff574f626124818bcb63 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 17 Mar 2024 15:32:13 -0400 Subject: [PATCH 05/18] Revise lithium ion battery example based on Jupyter version Co-authored-by: Steven DeCaluwe --- .../samples/single-particle-battery.svg | 167 ++++++ doc/sphinx/userguide/index.md | 5 + .../surface_chemistry/lithium_ion_battery.py | 480 ++++++++++++------ 3 files changed, 510 insertions(+), 142 deletions(-) create mode 100644 doc/sphinx/_static/images/samples/single-particle-battery.svg diff --git a/doc/sphinx/_static/images/samples/single-particle-battery.svg b/doc/sphinx/_static/images/samples/single-particle-battery.svg new file mode 100644 index 0000000000..3ec1807d5b --- /dev/null +++ b/doc/sphinx/_static/images/samples/single-particle-battery.svg @@ -0,0 +1,167 @@ + + + + + + + Graphite anode + Li + + + + + + + + + + + + + + + + + + + + + LCO cathode + Liquid electrolyte + + Li + x + CoO + 2 + + Li + x + C + 6 + + + + + + + + + Li + Li + + + Li + + + Li + + + Li + + + Li + + + Li + + + e + - + + + + + + + + + + + + + Li + Li + + + e + - + + + + + i + Far + i + ionic + + + + + + \ No newline at end of file diff --git a/doc/sphinx/userguide/index.md b/doc/sphinx/userguide/index.md index acf6fc6307..6cbe83ea5f 100644 --- a/doc/sphinx/userguide/index.md +++ b/doc/sphinx/userguide/index.md @@ -37,6 +37,9 @@ conditions, or calculating the voltage of a Lithium-ion battery as it is dischar - [](flame-temperature) - [](heating-value) +### Electrochemistry Calculations +- [](/examples/python/surface_chemistry/lithium_ion_battery) + ### Implementing Custom Models - [](extensible-reactor) @@ -73,5 +76,7 @@ legacy2yaml-tutorial flame-temperature heating-value +/examples/python/surface_chemistry/lithium_ion_battery + extensible-reactor ``` diff --git a/samples/python/surface_chemistry/lithium_ion_battery.py b/samples/python/surface_chemistry/lithium_ion_battery.py index f5a463db77..65f36dd6ed 100644 --- a/samples/python/surface_chemistry/lithium_ion_battery.py +++ b/samples/python/surface_chemistry/lithium_ion_battery.py @@ -1,177 +1,373 @@ -""" -Lithium-ion battery -=================== +r""" +Lithium Ion Battery Discharge Curve +=================================== + +In this example, we illustrate how to calculate the open circuit voltage (voltage when +the external applied current is zero) for a lithium ion battery as a function of anode +and cathode lithium content. The open circuit voltage here is calculated via two means: +kinetically and thermodynamically. + +The thermodynamics are based on a graphite anode and a LiCoO₂ cathode (the typical +active materials for commercial batteries, as of 2019), and are modeled using the +:ct:`BinarySolutionTabulatedThermo` class. -This example calculates the cell voltage of a lithium-ion battery at -given temperature, pressure, current, and range of state of charge (SOC). +The system here can be thought of as consisting of two particles—anode and +cathode—connected by a liquid electrolyte: -The thermodynamics are based on a graphite anode and a LiCoO2 cathode, -modeled using the :ct:`BinarySolutionTabulatedThermo` class. -Further required cell parameters are the electrolyte ionic resistance, the -stoichiometry ranges of the active materials (electrode balancing), and the -surface area of the active materials. +.. image:: /_static/images/samples/single-particle-battery.svg + :alt: Cartoon of a Single Particle Battery Model + :width: 75% + :align: center -The functionality of this example is presented in greater detail in a jupyter -notebook as well as the reference (which also describes the derivation of the -:ct:`BinarySolutionTabulatedThermo` class): +For the sake of simplicity, we're going to assume that the anode and cathode capacities +are perfectly balanced. That is, if the cathode lithium content is X percent of its max +possible (its capacity), then we will assume that the anode is at 1-X percent. Without +loss of generality, we will define the anode composition. -Reference: +The routine below returns the steady-state cell voltage of a lithium-ion cell for a +given cell current, active material lithium stoichiometries, and electrolyte ionic +resistance. This functionality is presented in greater detail in the reference (which +also describes the derivation of the :ct:`BinarySolutionTabulatedThermo` class): - M. Mayur, S. C. DeCaluwe, B. L. Kee, W. G. Bessler, “Modeling and simulation - of the thermodynamics of lithium-ion battery intercalation materials in the - open-source software Cantera,” Electrochim. Acta 323, 134797 (2019), - https://doi.org/10.1016/j.electacta.2019.134797 + M. Mayur, S. C. DeCaluwe, B. L. Kee, W. G. Bessler, “Modeling and simulation of the + thermodynamics of lithium-ion battery intercalation materials in the open-source + software Cantera,” *Electrochim. Acta* 323, 134797 (2019), + https://doi.org/10.1016/j.electacta.2019.134797 -Requires: cantera >= 2.6.0, matplotlib >= 2.0 +In the future, this example may be developed further to demonstrate simulating +charge-discharge of the lithium ion battery. + +Requires: cantera >= 2.6.0, matplotlib >= 2.0, scipy >= 1.10.0 .. tags:: Python, surface chemistry, kinetics, electrochemistry, battery, plotting + +Kinetic equilibrium calculations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +In the kinetic equilibrium problem, steady state is achieved when the faradaic current +at each electrode (anode and cathode) surface and the ionic current through the +electrolyte are all zero. There are, essentially, four relevant electric potentials +which must be determined: + +- :math:`\phi_{\rm anode}`: electric potential of graphite anode. +- :math:`\phi_{\rm elyte,\, anode}`: electric potential of electrolyte + at anode interface. +- :math:`\phi_{\rm elyte,\, cathode}`: electric potential of + electrolyte at cathode interface. +- :math:`\phi_{\rm cathode}`: electric potential of LCO cathode. + +Setting one of these four to the reference potential of zero (because it is the +difference in electric potential which drives currents, the actual potential values are +irrelevant. Let's assume the anode electric potential is zero), there is only one +distribution of electric potentials across the cell such that the current is invariant +across the cell. I.e. we want the potentials such that: + +.. math:: + + i_{\rm Far,\, anode} = i_{\rm ionic} = i_{\rm Far,\,cathode}= i_{\rm app} + +where :math:`i_{\rm app}` is the user input for the applied current. For this example, +we assume an applied current of 0, to calculate the equilibrium voltage. + +Faradaic current +^^^^^^^^^^^^^^^^ + +The faradaic current for this model is calculated using Butler-Volmer +kinetics. For a Li-ion battery, this is: + +.. math:: + + i = S_{\rm elde}i_\circ\left[\exp\left(\frac{F\beta\eta}{RT}\right) + - \exp\left(-\frac{F(1-\beta)\eta}{RT}\right) \right] + +where :math:`S_{\rm elde}` is the specific surface area of the electrode in question, +:math:`F` is Faraday's constant, :math:`\beta` is the charge-transfer symmetry +parameter, :math:`R` the universal gas constant, :math:`T` the temperature, and +:math:`\eta` the overpotential, which is the electric potential difference between the +electrode and electrolyte, :math:`\Delta \phi = \phi_{\rm elde} - \phi_{\rm elyte}`, +relative to that at equilibrium, :math:`\Delta \phi_{\rm eq}`: + +.. math:: + + \eta = \Delta \phi - \Delta \phi_{\rm eq} + +:math:`i_\circ` is known as the “exchange current density,” which is equal to the rate +of the forward and reverse current at equilibrium (which are equal). :math:`i_\circ` and +:math:`\beta` are provided as user inputs in the YAML file. At any particular state, +(user-specified electric potentials, pressure, temperature, and chemical compositions), +Cantera calculates :math:`\eta` as part of the routine to evaluate reaction rates of +progress :math:`\left(\dot{q} = i_{\rm Far}/F \right)`. The user simply sets the state +values mentioned above. + +Ionic current +^^^^^^^^^^^^^ + +The electrolyte is modeled as a resistor with user-defined ionic resistance +:math:`R_{\rm io}`, and hence the ionic current is calculated as: + +.. math:: + + i_{\rm ionic} = \frac{\phi_{\rm elyte,\,ca} - \phi_{\rm elyte,\,an}}{R_{\rm io}} + +where positive current is defined as delivering Li\ :math:`^+` to the anode interface. +Given :math:`i_{\rm app}`, this equation can be inverted, to calculate the electric +potential of the electrolyte at the cathode interface, relative to that at the anode +interface: + +.. math:: + + \phi_{\rm elyte,\,ca} = \phi_{\rm elyte,\,an} + R_{\rm io}i_{\rm app} + +Again: in this example, :math:`i_{\rm app} = 0` and hence the two electric potential +values in the electrolyte are equal. + +Numerical routine +^^^^^^^^^^^^^^^^^ + +For the kinetic routine, there are three processes to determine the cell voltage +:math:`\phi_{\rm cathode} - \phi_{\rm anode}` which corresponds to the user-provided +:math:`i_{\rm app}`: + +1. Determine the :math:`\phi_{\rm elyte,\,anode}` value which corresponds to + :math:`i_{\rm app}`, given :math:`X_{\rm Li, anode}`, the percentage of Li in the + anode active material. + +2. Determine :math:`\phi_{\rm elyte,\,cathode}`, given :math:`\phi_{\rm elyte,\,anode}` + and :math:`i_{\rm app}`. + +3. Determine the :math:`\phi_{\rm cathode}` which corresponds to :math:`i_{\rm app}`, + given :math:`\phi_{\rm elyte,\,cathode}` and :math:`X_{\rm Li, anode}`, the + percentage of Li in the anode active material. + +The routines below are written generally such that an interested user may set +:math:`i_{\rm app}` to any value of interest. + +Import necessary packages +^^^^^^^^^^^^^^^^^^^^^^^^^ + """ -import cantera as ct import numpy as np +from scipy.optimize import fsolve +import time # Used for timing our calculations +import matplotlib.pyplot as plt +import cantera as ct +print(f"Runnning Cantera version: {ct.__version__}") + +# %% +# Define the phases +# ~~~~~~~~~~~~~~~~~ +# +# The phase thermodynamics are defined according to experimentally-measured open circuit +# voltage values, as described in the reference provided above -# Parameters -samples = 101 -soc = np.linspace(0., 1., samples) # [-] Input state of charge (0...1) -current = -1 # [A] Externally-applied current, negative for discharge -T = 293 # [K] Temperature -P = ct.one_atm # [Pa] Pressure - -# Cell properties -input_file = "lithium_ion_battery.yaml" # Cantera input file name -R_electrolyte = 0.0384 # [Ohm] Electrolyte resistance -area_cathode = 1.1167 # [m^2] Cathode total active material surface area -area_anode = 0.7824 # [m^2] Anode total active material surface area - -# Electrode balancing: The "balancing" of the electrodes relates the chemical -# composition (lithium mole fraction in the active materials) to the macroscopic -# cell-level state of charge. -X_Li_anode_0 = 0.01 # [-] anode Li mole fraction at SOC = 0 -X_Li_anode_1 = 0.75 # [-] anode Li mole fraction at SOC = 100 -X_Li_cathode_0 = 0.99 # [-] cathode Li mole fraction at SOC = 0 -X_Li_cathode_1 = 0.49 # [-] cathode Li mole fraction at SOC = 100 - -# Calculate mole fractions from SOC -X_Li_anode = (X_Li_anode_1 - X_Li_anode_0) * soc + X_Li_anode_0 -X_Li_cathode = (X_Li_cathode_0 - X_Li_cathode_1) * (1 - soc) + X_Li_cathode_1 - -# Import all Cantera phases +input_file = "../data/lithium_ion_battery.yaml" anode = ct.Solution(input_file, "anode") cathode = ct.Solution(input_file, "cathode") -metal = ct.Solution(input_file, "electron") -electrolyte = ct.Solution(input_file, "electrolyte") -anode_int = ct.Interface( - input_file, "edge_anode_electrolyte", adjacent=[anode, metal, electrolyte]) -cathode_int = ct.Interface( - input_file, "edge_cathode_electrolyte", adjacent=[cathode, metal, electrolyte]) - -# Set the temperatures and pressures of all phases -for phase in [anode, cathode, metal, electrolyte, anode_int, cathode_int]: - phase.TP = T, P - - -# Root finding function -def newton_solve(f, xstart, C=0.0): - """ - Solve f(x) = C by Newton iteration using initial guess *xstart* - """ - f0 = f(xstart) - C - x0 = xstart - dx = 1.0e-6 - n = 0 - while n < 200: - ff = f(x0 + dx) - C - dfdx = (ff - f0) / dx - step = - f0 / dfdx - - # avoid taking steps too large - if abs(step) > 0.1: - step = 0.1 * step / abs(step) - - x0 += step - emax = 0.00001 # 0.01 mV tolerance - if abs(f0) < emax and n > 8: - return x0 - f0 = f(x0) - C - n += 1 - raise Exception("no root!") - - -# This function returns the Cantera calculated anode current (in A) -def anode_current(phi_s, phi_l, X_Li_anode): - """ - Current from the anode as a function of anode potential relative to - electrolyte. - """ +# The 'elde' electrode phase is needed as a source/sink for electron +elde = ct.Solution(input_file, "electron") +elyte = ct.Solution(input_file, "electrolyte") +anode_interface = ct.Interface( + input_file, "edge_anode_electrolyte", [anode, elde, elyte] +) +cathode_interface = ct.Interface( + input_file, "edge_cathode_electrolyte", [cathode, elde, elyte] +) + +# %% +# Define battery conditions: temperature, pressure, stoichiometry, electrolyte resistance +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Inputs are: +# +# - Stoichiometries ``X_Li_ca`` and ``X_Li_an`` [-] (can be vectors) +# - Temperature ``T`` [K] +# - Pressure ``P`` [Pa] +# - Externally-applied current ``i_app`` [A] +# - Electrolyte resistance ``R_elyte`` [Ohm] +# - Anode total surface area ``S_an`` [m^2] +# - Cathode total surface area ``S_ca`` [m^2] + +# Array of lithium mole fractions in the anode +X_Li_an = np.arange(0.005, 0.995, 0.02) +# Assume that the cathode and anode capacities are balanced +X_Li_ca = 1.0 - X_Li_an + +# I_app = 0: Open circuit +I_app = 0.0 + +# At zero current, electrolyte resistance is irrelevant +R_elyte = 0.0 + +# Temperature and pressure +T = 300 # K +P = ct.one_atm + +F = ct.faraday + +# [m^2] Cathode total active material surface area +S_ca = 1.1167 + +S_an = 0.7824 # [m^2] Anode total active material surface area + +# %% +# Set phase temperatures and pressures +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +phases = [anode, elde, elyte, cathode, anode_interface, cathode_interface] +for ph in phases: + ph.TP = T, P + +# %% +# Helper Functions +# ~~~~~~~~~~~~~~~~ + +def anode_curr(phi_l, I_app, phi_s, X_Li_an): # Set the active material mole fraction - anode.X = {"Li[anode]": X_Li_anode, "V[anode]": 1 - X_Li_anode} + anode.X = {"Li[anode]": X_Li_an, "V[anode]": 1 - X_Li_an} # Set the electrode and electrolyte potential - metal.electric_potential = phi_s - electrolyte.electric_potential = phi_l + elde.electric_potential = phi_s + elyte.electric_potential = phi_l[0] - # Get the net reaction rate at the anode-side interface - # Reaction according to input file: - # Li+[electrolyte] + V[anode] + electron <=> Li[anode] - r = anode_int.net_rates_of_progress[0] # [kmol/m2/s] + # Get the net product rate of electrons in the anode (per m2^ interface) + r_elec = anode_interface.get_net_production_rates(elde) - # Calculate the current. Should be negative for cell discharge. - return r * ct.faraday * area_anode + anode_current = r_elec * ct.faraday * S_an + diff = I_app + anode_current + return diff -# This function returns the Cantera calculated cathode current (in A) -def cathode_current(phi_s, phi_l, X_Li_cathode): - """ - Current to the cathode as a function of cathode potential relative to electrolyte - """ + +def cathode_curr(phi_s, I_app, phi_l, X_Li_ca): # Set the active material mole fractions - cathode.X = {"Li[cathode]": X_Li_cathode, "V[cathode]": 1 - X_Li_cathode} + cathode.X = {"Li[cathode]": X_Li_ca, "V[cathode]": 1 - X_Li_ca} # Set the electrode and electrolyte potential - metal.electric_potential = phi_s - electrolyte.electric_potential = phi_l + elde.electric_potential = phi_s[0] + elyte.electric_potential = phi_l - # Get the net reaction rate at the cathode-side interface - # Reaction according to input file: - # Li+[electrolyte] + V[cathode] + electron <=> Li[cathode] - r = cathode_int.net_rates_of_progress[0] # [kmol/m2/s] + # Get the net product rate of electrons in the cathode (per m2^ interface) + r_elec = cathode_interface.get_net_production_rates(elde) - # Calculate the current. Should be negative for cell discharge. - return - r * ct.faraday * area_cathode + cathode_current = r_elec * ct.faraday * S_an + diff = I_app - cathode_current + return diff -# Calculate cell voltage, separately for each entry of the input vectors -V_cell = np.zeros_like(soc) -phi_l_anode = 0 -phi_s_cathode = 0 -for i in range(samples): - # Set anode electrode potential to 0 - phi_s_anode = 0 +# %% +# Run the calculations for all stoichiometries +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # Calculate anode electrolyte potential - phi_l_anode = newton_solve( - lambda E: anode_current(phi_s_anode, E, X_Li_anode[i]), - phi_l_anode, C=current) +# Tic +t0 = time.time() - # Calculate cathode electrolyte potential - phi_l_cathode = phi_l_anode + current * R_electrolyte +# Initialize array of OCVs +E_cell_kin = np.zeros_like(X_Li_ca) - # Calculate cathode electrode potential - phi_s_cathode = newton_solve( - lambda E: cathode_current(E, phi_l_cathode, X_Li_cathode[i]), - phi_s_cathode, C=current) +for i, X_an in enumerate(X_Li_an): + # Set anode electrode potential to 0 + phi_s_an = 0 + E_init = 3.0 - # Calculate cell voltage - V_cell[i] = phi_s_cathode - phi_s_anode + phi_l_an = fsolve(anode_curr, E_init, args=(I_app, phi_s_an, X_an)) -try: - import matplotlib.pyplot as plt + # Calculate electrolyte potential at cathode interface + phi_l_ca = phi_l_an[0] + I_app * R_elyte - # Plot the cell voltage, as a function of the state of charge - plt.plot(soc * 100, V_cell) - plt.xlabel("State of charge / %") - plt.ylabel("Cell voltage / V") - plt.show() + # Calculate cathode electrode potential + phi_s_ca = fsolve(cathode_curr, E_init, args=(I_app, phi_l_ca, X_Li_ca[i])) -except ImportError: - print("Install matplotlib to plot the outputs") + # Calculate cell voltage + E_cell_kin[i] = phi_s_ca[0] - phi_s_an + +# Toc +t1 = time.time() +print(f"{i:d} cell voltages calculated in {t1 - t0:3.2e} seconds.") + +# %% +# Plot cell voltage, as a function of the cathode stoichiometry +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +fig, ax = plt.subplots() +ax.plot(100 * X_Li_ca, E_cell_kin) +_ = ax.set(ylim=[2.5, 4.3], xlabel="Li Fraction in Cathode (%)", + ylabel="Open Circuit Potential (V)") + +# %% +# Thermodynamic Equilibrium Calculation +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# For the :math:`i_{app} = 0` case, we can also calculate the voltage using +# thermodynamics. At equilibrium, the net electrochemical potential change of the +# reaction must be zero: +# +# .. math:: +# \sum_k\nu_k\tilde{\mu}_k = 0 +# +# where :math:`\tilde{\mu}_k = \mu_k + z_kF\Phi_k`, where, in turn +# :math:`\mu_k = \frac{\partial g_k}{\partial n_k}` is the chemical potential, +# :math:`\nu_k` the net stoichiometric coefficient, :math:`z_k` the net elementary +# charge, and :math:`\Phi_k` the phase electric potential for species :math:`k`. +# +# From this, we can calculate the equilibrium electric potential difference +# :math:`\Delta \Phi_{\rm eq} = \left(\Phi_{\rm elde} - \Phi_{\rm elyte}\right)_{\rm eq}` +# as: +# +# .. math:: +# \Delta \Phi_{\rm eq} = -\frac{\Delta g_{\rm rxn}}{n_{\rm charge}F} +# +# where :math:`\Delta g_{\rm rxn} = \sum_k \nu_k\mu_k` is the chemical potential of the +# reaction and and :math:`n_{\rm charge} = \sum_{k,\,{\rm elde}} \nu_k z_k` is the net +# elementary charge transferred from the electrolyte to the electrode. + +# Tic +t0 = time.time() + +# Initialize array of OCVs +E_cell_therm = np.zeros_like(X_Li_ca) + +for i, X_an in enumerate(X_Li_an): + # Set anode electrode potential to 0 + anode.X = "Li[anode]:" + str(X_an) + ", V[anode]:" + str(1 - X_an) + dG_an = anode_interface.delta_gibbs[0] + n_charge = 1.0 + E_eq_an = -dG_an / n_charge / ct.faraday + + cathode.X = "Li[cathode]:" + str(1.0 - X_an) + ", V[cathode]:" + str(X_an) + dG_ca = cathode_interface.delta_gibbs[0] + n_charge = 1.0 + E_eq_ca = -dG_ca / n_charge / ct.faraday + + E_cell_therm[i] = E_eq_ca - E_eq_an + +# Toc +t1 = time.time() +print(f"{i:d} cell voltages calculated in {t1 - t0:3.2e} seconds.") + +# %% +# Plot thermodynamic OCV, and compare to results from kinetic method +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +fig, ax = plt.subplots() +ax.plot(100 * X_Li_ca, E_cell_kin, label="Kinetic") +ax.plot(100 * X_Li_ca, E_cell_therm, + linewidth=0.0, marker="o", markerfacecolor="none", label="Thermodynamic") +ax.set(ylim=[2.5, 4.3], xlabel="Li Fraction in Cathode (%)", + ylabel="Open Circuit Potential (V)") +_ = ax.legend() + +# %% +# As one would expect, the two approaches give identical results. While both methods are +# incredibly fast, the thermodynamic method is roughly 30 times faster. +# +# A large part of this is that the thermodynamic approach is an analytical approach +# (i.e. the answer is known from theory), while the kinetic approach relies on the +# root-finding fzero method to fit the correct voltage. Note also that the kinetic +# method, because of the use of Butler-Volmer kinetics, calculates the thermodynamic +# voltage, in order to calculate the overpotential :math:`\eta = \Delta \Phi - \Delta +# \Phi_{\rm eq}`. +# +# However, it is at last important to note that, while slower, the kinetic method is of +# course more robust, and can be used to find results away from equilibrium. The +# thermodynamic method is only applicable at equilibrium (zero current). From 9d7f34297a839a5d90903970d071f1747b2b4789 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 18 Mar 2024 23:49:46 -0400 Subject: [PATCH 06/18] Revise flame speed sensitivity example based on Jupyter version Co-authored-by: Santosh Shanbhogue --- .../_static/images/samples/flame-speed.svg | 120 +++++++++++++ .../python/onedim/flamespeed_sensitivity.py | 164 +++++++++++++++--- 2 files changed, 260 insertions(+), 24 deletions(-) create mode 100644 doc/sphinx/_static/images/samples/flame-speed.svg diff --git a/doc/sphinx/_static/images/samples/flame-speed.svg b/doc/sphinx/_static/images/samples/flame-speed.svg new file mode 100644 index 0000000000..38b5ef653f --- /dev/null +++ b/doc/sphinx/_static/images/samples/flame-speed.svg @@ -0,0 +1,120 @@ + + + + + + + + + + + + + + + + + + Reactantsρ + u + T + u + S + u + Productsρ + b + T + b + S + b + + + + + + + + + + + + + + + + Flame + \ No newline at end of file diff --git a/samples/python/onedim/flamespeed_sensitivity.py b/samples/python/onedim/flamespeed_sensitivity.py index bdbaea11ad..0d9433677b 100644 --- a/samples/python/onedim/flamespeed_sensitivity.py +++ b/samples/python/onedim/flamespeed_sensitivity.py @@ -1,41 +1,157 @@ -""" +r""" Laminar flame speed sensitivity analysis ======================================== -Sensitivity analysis for a freely-propagating, premixed methane-air -flame. Computes the sensitivity of the laminar flame speed with respect -to each reaction rate constant. +In this example we simulate a freely-propagating, adiabatic, premixed methane-air flame, +calculate its laminar burning velocity and perform a sensitivity analysis of its +kinetics with respect to each reaction rate constant. + +Requires: cantera >= 3.0.0, pandas + +.. tags:: Python, combustion, 1D flow, flame speed, premixed flame, + sensitivity analysis, plotting + +The figure below illustrates the setup, in a flame-fixed coordinate system. The +reactants enter with density :math:`\rho_u`, temperature :math:`T_u` and speed +:math:`S_u`. The products exit the flame at speed :math:`S_b`, density :math:`\rho_b` +and temperature :math:`T_b`. -Requires: cantera >= 2.5.0 +.. image:: /_static/images/samples/flame-speed.svg + :width: 50% + :alt: Freely Propagating Flame + :align: center -.. tags:: Python, combustion, 1D flow, flame speed, premixed flame, sensitivity analysis """ import cantera as ct +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +plt.rcParams['figure.constrained_layout.use'] = True + +# %% # Simulation parameters -p = ct.one_atm # pressure [Pa] -Tin = 300.0 # unburned gas temperature [K] -reactants = 'CH4:0.45, O2:1.0, N2:3.76' +# --------------------- + +# Define the gas mixture and kinetics used to compute mixture properties +# In this case, we are using the GRI 3.0 model with methane as the fuel +mech_file = "gri30.yaml" +fuel_comp = {"CH4": 1.0} + +# Inlet temperature in kelvin and inlet pressure in pascals +# In this case we are setting the inlet T and P to room temperature conditions +To = 300 +Po = ct.one_atm + +# Domain width in metres +width = 0.014 + +# %% +# Simulation setup +# ---------------- + +# Create the object representing the gas and set its state to match the inlet conditions +gas = ct.Solution(mech_file) +gas.TP = To, Po +gas.set_equivalence_ratio(1.0, fuel_comp, {"O2": 1.0, "N2": 3.76}) + +flame = ct.FreeFlame(gas, width=width) +flame.set_refine_criteria(ratio=3, slope=0.07, curve=0.14) + +# %% +# Solve the flame +# --------------- + +flame.solve(loglevel=1, auto=True) + +# %% +print(f"\nmixture-averaged flame speed = {flame.velocity[0]:7f} m/s\n") -width = 0.03 # m -# Solution object used to compute mixture properties -gas = ct.Solution('gri30.yaml') -gas.TPX = Tin, p, reactants +# %% +# Plot temperature and major species profiles +# ------------------------------------------- +# +# Check and see if all has gone well. -# Flame object -f = ct.FreeFlame(gas, width=width) -f.set_refine_criteria(ratio=3, slope=0.07, curve=0.14) +# Extract the spatial profiles as a SolutionArray to simplify plotting specific species +profile = flame.to_array() -f.solve(loglevel=1, auto=True) -print(f"\nmixture-averaged flamespeed = {f.velocity[0]:7f} m/s\n") +fig, ax = plt.subplots() + +ax.plot(profile.grid * 100, profile.T, ".-") +_ = ax.set(xlabel="Distance [cm]", ylabel="Temperature [K]") + +# %% + +fig, ax = plt.subplots() + +ax.plot(profile.grid * 100, profile("CH4").X, "--", label="CH$_4$") +ax.plot(profile.grid * 100, profile("O2").X, label="O$_2$") +ax.plot(profile.grid * 100, profile("CO2").X, "--", label="CO$_2$") +ax.plot(profile.grid * 100, profile("H2O").X, label="H$_2$O") + +ax.legend(loc="best") +plt.xlabel("Distance (cm)") +_ = ax.set(xlabel="Distance [cm]", ylabel="Mole fraction [-]") + +# %% +# Sensitivity Analysis +# -------------------- +# +# See which reactions affect the flame speed the most. The sensitivities can be +# efficiently computed using the adjoint method. In the general case, we consider a +# system :math:`f(x, p) = 0` where :math:`x` is the system's state vector and :math:`p` +# is a vector of parameters. To compute the sensitivities of a scalar function +# :math:`g(x, p)` to the parameters, we solve the system: +# +# .. math:: +# \left(\frac{\partial f}{\partial x}\right)^T \lambda = +# \left(\frac{\partial g}{\partial x}\right)^T +# +# for the Lagrange multiplier vector :math:`\lambda`. The sensitivities are then +# computed as +# +# .. math:: +# \left.\frac{dg}{dp}\right|_{f=0} = +# \frac{\partial g}{\partial p} - \lambda^T \frac{\partial f}{\partial p} +# +# In the case of flame speed sensitivity to the reaction rate coefficients, +# :math:`g = S_u` and :math:`\partial g/\partial p = 0`. Since +# :math:`\partial f/\partial x` is already computed as part of the normal solution +# process, computing the sensitivities for :math:`N` reactions only requires :math:`N` +# additional evaluations of the residual function to obtain +# :math:`\partial f/\partial p`, plus the solution of the linear system for +# :math:`\lambda`. +# +# Calculation of flame speed sensitivities with respect to rate coefficients is +# implemented by the method `FreeFlame.get_flame_speed_reaction_sensitivities`. For +# other sensitivities, the method `Sim1D.solve_adjoint` can be used. + +# Create a DataFrame to store sensitivity-analysis data +sens = pd.DataFrame(index=gas.reaction_equations(), columns=["sensitivity"]) # Use the adjoint method to calculate sensitivities -sens = f.get_flame_speed_reaction_sensitivities() +sens.sensitivity = flame.get_flame_speed_reaction_sensitivities() + +# Show the first 10 sensitivities +sens.head(10) + +# %% + +# Sort the sensitivities in order of descending magnitude +sens = sens.iloc[(-sens['sensitivity'].abs()).argsort()] + +fig, ax = plt.subplots() + +# Reaction mechanisms can contains thousands of elementary steps. Limit the plot +# to the top 15 +sens.head(15).plot.barh(ax=ax, legend=None) + +ax.invert_yaxis() # put the largest sensitivity on top +ax.set_title("Sensitivities for GRI 3.0") +ax.set_xlabel(r"Sensitivity: $\frac{\partial\:\ln S_u}{\partial\:\ln k}$") +ax.grid(axis='x') -print() -print('Rxn # k/S*dS/dk Reaction Equation') -print('----- ---------- ----------------------------------') -for m in range(gas.n_reactions): - print(f"{m: 5d} {sens[m]: 10.3e} {gas.reaction(m).equation}") +#sphinx_gallery_thumbnail_number = -1 From 82060c34872766ad69da447b4d43d80e8fb4fd1b Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 20 Mar 2024 09:24:51 -0400 Subject: [PATCH 07/18] [Doc] Automatically ignore unwanted Matplotlib output in gallery --- doc/sphinx/conf.py | 1 + samples/python/onedim/flamespeed_sensitivity.py | 7 ++++--- samples/python/surface_chemistry/lithium_ion_battery.py | 2 +- samples/python/thermo/equations_of_state.py | 2 +- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py index d595c8bdb4..07698ebce3 100644 --- a/doc/sphinx/conf.py +++ b/doc/sphinx/conf.py @@ -53,6 +53,7 @@ 'ignore_pattern': r'(__.*__\.py|test_examples\.m)', 'image_srcset': ["2x"], 'remove_config_comments': True, + 'ignore_repr_types': r'matplotlib\.(text|axes|legend)', 'examples_dirs': [ '../samples/python/', '../samples/cxx/', diff --git a/samples/python/onedim/flamespeed_sensitivity.py b/samples/python/onedim/flamespeed_sensitivity.py index 0d9433677b..558a690632 100644 --- a/samples/python/onedim/flamespeed_sensitivity.py +++ b/samples/python/onedim/flamespeed_sensitivity.py @@ -81,7 +81,8 @@ fig, ax = plt.subplots() ax.plot(profile.grid * 100, profile.T, ".-") -_ = ax.set(xlabel="Distance [cm]", ylabel="Temperature [K]") +ax.set_xlabel("Distance [cm]") +ax.set_ylabel("Temperature [K]") # %% @@ -93,8 +94,8 @@ ax.plot(profile.grid * 100, profile("H2O").X, label="H$_2$O") ax.legend(loc="best") -plt.xlabel("Distance (cm)") -_ = ax.set(xlabel="Distance [cm]", ylabel="Mole fraction [-]") +ax.set_xlabel("Distance [cm]") +ax.set_ylabel("Mole fraction [-]") # %% # Sensitivity Analysis diff --git a/samples/python/surface_chemistry/lithium_ion_battery.py b/samples/python/surface_chemistry/lithium_ion_battery.py index 65f36dd6ed..46a0511964 100644 --- a/samples/python/surface_chemistry/lithium_ion_battery.py +++ b/samples/python/surface_chemistry/lithium_ion_battery.py @@ -355,7 +355,7 @@ def cathode_curr(phi_s, I_app, phi_l, X_Li_ca): linewidth=0.0, marker="o", markerfacecolor="none", label="Thermodynamic") ax.set(ylim=[2.5, 4.3], xlabel="Li Fraction in Cathode (%)", ylabel="Open Circuit Potential (V)") -_ = ax.legend() +ax.legend() # %% # As one would expect, the two approaches give identical results. While both methods are diff --git a/samples/python/thermo/equations_of_state.py b/samples/python/thermo/equations_of_state.py index 3790be7d34..731cb5edf7 100644 --- a/samples/python/thermo/equations_of_state.py +++ b/samples/python/thermo/equations_of_state.py @@ -253,7 +253,7 @@ def plot(p, thermo_Ideal, thermo_RK, thermo_CoolProp, name): ax.text(960, 320, "p = 600 bar", color=color[9], rotation=-68) ax.set_xlabel("Density [kg/m$^3$]") ax.set_ylabel("Temperature [K]") -_ = ax.legend(handles=[ideal_line[0], RK_line[0], CP_line[0]]) +ax.legend(handles=[ideal_line[0], RK_line[0], CP_line[0]]) #%% # The figure compares :math:`T-\rho` plots for ideal, R-K, and Helmholtz EoS at From b54784df81cda7cce73049e9d29d78dca78c409d Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 20 Mar 2024 20:53:21 -0400 Subject: [PATCH 08/18] Migrate flame convergence analysis example from Jupyter Co-authored-by: Richard West --- .../flame_speed_convergence_analysis.py | 520 ++++++++++++++++++ 1 file changed, 520 insertions(+) create mode 100644 samples/python/onedim/flame_speed_convergence_analysis.py diff --git a/samples/python/onedim/flame_speed_convergence_analysis.py b/samples/python/onedim/flame_speed_convergence_analysis.py new file mode 100644 index 0000000000..5b1f41b7a1 --- /dev/null +++ b/samples/python/onedim/flame_speed_convergence_analysis.py @@ -0,0 +1,520 @@ +r""" +Flame Speed with Convergence Analysis +===================================== + +Requires: cantera >= 3.0.0, matplotlib >= 2.0, pandas + +In this example we simulate a freely-propagating, adiabatic, 1-D flame and + +* Calculate its laminar burning velocity +* Estimate the uncertainty in the laminar burning velocity calculation due to grid size. + +.. tags:: Python, combustion, 1D flow, flame speed, premixed flame, plotting + +The figure below illustrates the setup, in a flame-fixed co-ordinate system. The +reactants enter with density :math:`\rho_u`, temperature :math:`T_u` and speed +:math:`S_u`. The products exit the flame at speed :math:`S_b`, density :math:`\rho_b`, +and temperature :math:`T_b`. + +.. image:: /_static/images/samples/flame-speed.svg + :width: 50% + :alt: Freely Propagating Flame + :align: center +""" + +# %% +# Import Modules +# -------------- + +import cantera as ct +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt +import matplotlib +import scipy +import scipy.optimize + +# %% +# Define plotting preference +# -------------------------- + +plt.style.use("ggplot") +plt.style.use("seaborn-v0_8-deep") +plt.rcParams["figure.constrained_layout.use"] = True + +# %% +# Estimate uncertainty from grid size and speeds +# ---------------------------------------------- + +def extrapolate_uncertainty(grids, speeds, plot=True): + """ + Given a list of grid sizes and a corresponding list of flame speeds, + extrapolate and estimate the uncertainty in the final flame speed. + Also makes a plot, unless called with `plot=False`. + """ + grids = list(grids) + speeds = list(speeds) + + def speed_from_grid_size(grid_size, true_speed, error): + """ + Given a grid size (or an array or list of grid sizes) + return a prediction (or array of predictions) + of the computed flame speed, based on + the parameters `true_speed` and `error`. + + It seems, from experience, that error scales roughly with + 1/grid_size, so we assume that form. + """ + return true_speed + error / np.array(grid_size) + + # Fit the chosen form of speed_from_grid_size, to the last four + # speed and grid size values. + popt, pcov = scipy.optimize.curve_fit(speed_from_grid_size, grids[-4:], speeds[-4:]) + + # How bad the fit was gives you some error, `percent_error_in_true_speed`. + perr = np.sqrt(np.diag(pcov)) + true_speed_estimate = popt[0] + percent_error_in_true_speed = perr[0] / popt[0] + print( + f"Fitted true_speed is {popt[0] * 100:.4f} ± {perr[0] * 100:.4f} cm/s " + f"({percent_error_in_true_speed:.1%})" + ) + + # How far your extrapolated infinite grid value is from your extrapolated + # (or interpolated) final grid value, gives you some other error, `estimated_percent_error` + estimated_percent_error = ( + speed_from_grid_size(grids[-1], *popt) - true_speed_estimate + ) / true_speed_estimate + print(f"Estimated error in final calculation {estimated_percent_error:.1%}") + + # The total estimated error is the sum of these two errors. + total_percent_error_estimate = abs(percent_error_in_true_speed) + abs( + estimated_percent_error + ) + print(f"Estimated total error {total_percent_error_estimate:.1%}") + + if plot: + fig, ax = plt.subplots() + ax.semilogx(grids, speeds, "o-") + ax.set_ylim( + min(speeds[-5:] + [true_speed_estimate - perr[0]]) * 0.95, + max(speeds[-5:] + [true_speed_estimate + perr[0]]) * 1.05, + ) + ax.plot(grids[-4:], speeds[-4:], "or") + extrapolated_grids = grids + [grids[-1] * i for i in range(2, 8)] + ax.plot( + extrapolated_grids, speed_from_grid_size(extrapolated_grids, *popt), ":r" + ) + ax.set_xlim(*ax.get_xlim()) # Prevent automatic expansion of axis limits + ax.hlines(true_speed_estimate, *ax.get_xlim(), colors="r", linestyles="dashed") + ax.hlines( + true_speed_estimate + perr[0], + *ax.get_xlim(), + colors="r", + linestyles="dashed", + alpha=0.3, + ) + ax.hlines( + true_speed_estimate - perr[0], + *ax.get_xlim(), + colors="r", + linestyles="dashed", + alpha=0.3, + ) + ax.fill_between( + ax.get_xlim(), + true_speed_estimate - perr[0], + true_speed_estimate + perr[0], + facecolor="red", + alpha=0.1, + ) + + above = popt[1] / abs( + popt[1] + ) # will be +1 if approach from above or -1 if approach from below + + ax.annotate( + "", + xy=(grids[-1], true_speed_estimate), + xycoords="data", + xytext=(grids[-1], speed_from_grid_size(grids[-1], *popt)), + textcoords="data", + arrowprops=dict( + arrowstyle="|-|, widthA=0.5, widthB=0.5", + linewidth=1, + connectionstyle="arc3", + color="black", + shrinkA=0, + shrinkB=0, + ), + ) + + ax.annotate( + f"{abs(estimated_percent_error):.1%}", + xy=(grids[-1], speed_from_grid_size(grids[-1], *popt)), + xycoords="data", + xytext=(5, 15 * above), + va="center", + textcoords="offset points", + arrowprops=dict(arrowstyle="->", connectionstyle="arc3"), + ) + + ax.annotate( + "", + xy=(grids[-1] * 4, true_speed_estimate - (above * perr[0])), + xycoords="data", + xytext=(grids[-1] * 4, true_speed_estimate), + textcoords="data", + arrowprops=dict( + arrowstyle="|-|, widthA=0.5, widthB=0.5", + linewidth=1, + connectionstyle="arc3", + color="black", + shrinkA=0, + shrinkB=0, + ), + ) + ax.annotate( + f"{abs(percent_error_in_true_speed):.1%}", + xy=(grids[-1] * 4, true_speed_estimate - (above * perr[0])), + xycoords="data", + xytext=(5, -15 * above), + va="center", + textcoords="offset points", + arrowprops=dict(arrowstyle="->", connectionstyle="arc3"), + ) + + ax.set(xlabel="Grid size", ylabel="Flame speed (m/s)") + + return true_speed_estimate, total_percent_error_estimate + + +# %% +def make_callback(flame): + """ + Create and return a callback function that you will attach to + a flame solver. The reason we define a function to make the callback function, + instead of just defining the callback function, is so that it can store + a pair of lists that persist between function calls, to store the + values of grid size and flame speed. + + This factory returns the callback function, and the two lists: + (callback, speeds, grids) + """ + speeds = [] + grids = [] + + def callback(_): + speed = flame.velocity[0] + grid = len(flame.grid) + speeds.append(speed) + grids.append(grid) + print(f"Iteration {len(grids)}") + print(f"Current flame speed is is {speed * 100:.4f} cm/s") + if len(grids) < 5: + return 1.0 # + try: + extrapolate_uncertainty(grids, speeds) + except Exception as e: + print("Couldn't estimate uncertainty. " + str(e)) + return 1.0 # continue anyway + return 1.0 + + return callback, speeds, grids + + +# %% +# Define the reactant conditions, gas mixture and kinetic mechanism associated with the gas +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +# Inlet Temperature in Kelvin and Inlet Pressure in Pascals +# In this case we are setting the inlet T and P to room temperature conditions +To = 300 +Po = 101325 + +# Define the gas mixture and kinetics +# In this case, we are choosing a GRI3.0 gas +gas = ct.Solution("gri30.yaml") + +# Create a stoichiometric CH4/Air premixed mixture +gas.set_equivalence_ratio(1.0, "CH4", {"O2": 1.0, "N2": 3.76}) +gas.TP = To, Po + +# %% +# Define flame simulation conditions +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +# Domain width in metres +width = 0.014 + +# Create the flame object +flame = ct.FreeFlame(gas, width=width) + +# Define logging level +loglevel = 1 + +# Define tight tolerances for the solver +refine_criteria = {"ratio": 2, "slope": 0.01, "curve": 0.01} +flame.set_refine_criteria(**refine_criteria) + +# Set maximum number of grid points to be very high (otherwise default is 1000) +flame.set_max_grid_points(flame.domains[flame.domain_index("flame")], 1e4) + +# Set up the the callback function and lists of speeds and grids +callback, speeds, grids = make_callback(flame) +flame.set_steady_callback(callback) + +# %% +# Solve +# ~~~~~ +# +# After the first five iterations, it will start to estimate the uncertainty. + +flame.solve(loglevel=loglevel, auto=True) + +Su0 = flame.velocity[0] +print(f"Flame Speed is: {Su0 * 100:.2f} cm/s") + +# %% +# Use the final lists of grid sizes and flame speeds to make one final extrapolation +# "best guess" + +best_true_speed_estimate, best_total_percent_error_estimate = extrapolate_uncertainty( + grids, speeds +) + +best_true_speed_estimate + + +# %% +# Analyze the error predictions +# ----------------------------- +# +# Now let's see how good our error estimates were, with hindsight. +# +# If we assume that the final answer, with a very fine grid, has actually converged and +# is is the "truth", then we can find out how large the errors were in the previous +# values, and compare these with our estimated errors. This will show if our estimates +# are reasonable, or conservative, or too optimistic. + +def analyze_errors(grids, speeds, true_speed): + """ + If we assume that the final answer, with a very fine grid, + has actually converged and is is the "truth", then we can + find out how large the errors were in the previous values, + and compare these with our estimated errors. + This will show if our estimates are reasonable, or conservative, or too optimistic. + """ + true_speed_estimates = np.full_like(speeds, np.NaN) + total_percent_error_estimates = np.full_like(speeds, np.NaN) + actual_extrapolated_percent_errors = np.full_like(speeds, np.NaN) + actual_raw_percent_errors = np.full_like(speeds, np.NaN) + for i in range(3, len(grids)): + print(grids[: i + 1]) + true_speed_estimate, total_percent_error_estimate = extrapolate_uncertainty( + grids[: i + 1], speeds[: i + 1], plot=False + ) + actual_extrapolated_percent_error = ( + abs(true_speed_estimate - true_speed) / true_speed + ) + actual_raw_percent_error = abs(speeds[i] - true_speed) / true_speed + print( + "Actual extrapolated error (with hindsight) " + f"{actual_extrapolated_percent_error:.1%}" + ) + print(f"Actual raw error (with hindsight) {actual_raw_percent_error:.1%}") + + true_speed_estimates[i] = true_speed_estimate + total_percent_error_estimates[i] = total_percent_error_estimate + actual_extrapolated_percent_errors[i] = actual_extrapolated_percent_error + actual_raw_percent_errors[i] = actual_raw_percent_error + print() + + fig, ax = plt.subplots() + ax.loglog(grids, actual_raw_percent_errors * 100, "o-", label="raw error") + ax.loglog( + grids, + actual_extrapolated_percent_errors * 100, + "o-", + label="extrapolated error", + ) + ax.loglog( + grids, total_percent_error_estimates * 100, "o-", label="estimated error" + ) + ax.set(xlabel="Grid size", ylabel="Error in flame speed (%)") + ax.legend() + ax.set_title(flame.get_refine_criteria()) + ax.get_yaxis().set_major_formatter(matplotlib.ticker.PercentFormatter()) + flame.get_refine_criteria() + + data = pd.DataFrame( + data={ + "actual error in raw value": actual_raw_percent_errors * 100, + "actual error in extrapolated value": actual_extrapolated_percent_errors + * 100, + "estimated error": total_percent_error_estimates * 100, + }, + index=grids, + ) + return data + + +analyze_errors(grids, speeds, best_true_speed_estimate) + +# %% +# Repeat with less tight refine criteria +# -------------------------------------- + +refine_criteria = {"ratio": 3, "slope": 0.1, "curve": 0.1} + +# Reset the gas +gas.set_equivalence_ratio(1.0, "CH4", {"O2": 1.0, "N2": 3.76}) +gas.TP = To, Po + +# Create a new flame object +flame = ct.FreeFlame(gas, width=width) + +flame.set_refine_criteria(**refine_criteria) +flame.set_max_grid_points(flame.domains[flame.domain_index("flame")], 1e4) + +callback, speeds, grids = make_callback(flame) +flame.set_steady_callback(callback) + +# Define logging level +loglevel = 1 + +flame.solve(loglevel=loglevel, auto=True) + +Su0 = flame.velocity[0] +print(f"Flame Speed is: {Su0 * 100:.2f} cm/s") + +# Use the best true speed estimate from the fine grid tight criteria above +analyze_errors(grids, speeds, best_true_speed_estimate) + +# %% +# Default (loose) criteria +# ------------------------ + +flame = ct.FreeFlame(gas, width=width) +refine_criteria = flame.get_refine_criteria() +refine_criteria.update({"prune": 0}) +refine_criteria + +gas.set_equivalence_ratio(1.0, "CH4", {"O2": 1.0, "N2": 3.76}) +gas.TP = To, Po + +# Create a new flame object +flame = ct.FreeFlame(gas, width=width) + +flame.set_refine_criteria(**refine_criteria) +flame.set_max_grid_points(flame.domains[flame.domain_index("flame")], 1e4) + +callback, speeds, grids = make_callback(flame) +flame.set_steady_callback(callback) + +# Define logging level +loglevel = 1 + +flame.solve(loglevel=loglevel, auto=True) + +Su0 = flame.velocity[0] +print(f"Flame Speed is: {Su0 * 100:.2f} cm/s") + +analyze_errors(grids, speeds, best_true_speed_estimate) + +# %% +# Middling refine criteria +# ------------------------ + +refine_criteria = {"ratio": 3, "slope": 0.1, "curve": 0.1} + +# Reset the gas +gas.set_equivalence_ratio(1.0, "CH4", {"O2": 1.0, "N2": 3.76}) +gas.TP = To, Po + +# Create a new flame object +flame = ct.FreeFlame(gas, width=width) + +flame.set_refine_criteria(**refine_criteria) +flame.set_max_grid_points(flame.domains[flame.domain_index("flame")], 1e4) + +callback, speeds, grids = make_callback(flame) +flame.set_steady_callback(callback) + +# Define logging level +loglevel = 1 + +flame.solve(loglevel=loglevel, auto=True) + +Su0 = flame.velocity[0] +print(f"Flame Speed is: {Su0 * 100:.2f} cm/s") + +# %% +analyze_errors(grids, speeds, best_true_speed_estimate) + +# %% +# Try a Hydrogen flame (still with GRI mech) +# ------------------------------------------ + +# Tight criteria +refine_criteria = {"ratio": 2, "slope": 0.01, "curve": 0.01} + +# Reset the gas +gas.set_equivalence_ratio(1.0, "H2", {"O2": 1.0, "N2": 3.76}) +gas.TP = To, Po + +# Create a new flame object +flame = ct.FreeFlame(gas, width=width) + +flame.set_refine_criteria(**refine_criteria) +flame.set_max_grid_points(flame.domains[flame.domain_index("flame")], 1e4) + +callback, speeds, grids = make_callback(flame) +flame.set_steady_callback(callback) + +# Define logging level +loglevel = 1 + +flame.solve(loglevel=loglevel, auto=True) + +Su0 = flame.velocity[0] +print(f"Flame Speed is: {Su0 * 100:.2f} cm/s") + +# %% +# get a new best true speed estimate +best_true_speed_estimate, best_total_percent_error_estimate = extrapolate_uncertainty( + grids, speeds +) + +# %% +analyze_errors(grids, speeds, best_true_speed_estimate) + +# %% +# Middling refine criteria, Hydrogen flame +# ---------------------------------------- + +refine_criteria = {"ratio": 3, "slope": 0.1, "curve": 0.1} + +# Reset the gas +gas.set_equivalence_ratio(1.0, "H2", {"O2": 1.0, "N2": 3.76}) +gas.TP = To, Po + +# Create a new flame object +flame = ct.FreeFlame(gas, width=width) + +flame.set_refine_criteria(**refine_criteria) +flame.set_max_grid_points(flame.domains[flame.domain_index("flame")], 1e4) + +callback, speeds, grids = make_callback(flame) +flame.set_steady_callback(callback) + +# Define logging level +loglevel = 1 + +flame.solve(loglevel=loglevel, auto=True) + +Su0 = flame.velocity[0] +print(f"Flame Speed is: {Su0 * 100:.2f} cm/s") + +# %% +analyze_errors(grids, speeds, best_true_speed_estimate) From 234d91035121a1cab1823871fb727790aab8f2ce Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 22 Mar 2024 13:43:12 -0400 Subject: [PATCH 09/18] Revise twin flame example based on Jupyter Re-draw figure from Jupyter example as an SVG. Co-authored-by: Santosh Shanbhogue --- .../images/samples/twin-premixed-flame.svg | 153 ++++++++++++++++++ .../onedim/premixed_counterflow_twin_flame.py | 125 ++++++++------ 2 files changed, 229 insertions(+), 49 deletions(-) create mode 100644 doc/sphinx/_static/images/samples/twin-premixed-flame.svg diff --git a/doc/sphinx/_static/images/samples/twin-premixed-flame.svg b/doc/sphinx/_static/images/samples/twin-premixed-flame.svg new file mode 100644 index 0000000000..55b997d82d --- /dev/null +++ b/doc/sphinx/_static/images/samples/twin-premixed-flame.svg @@ -0,0 +1,153 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Flames + Products + reactants(fuel + oxidizer) + reactants(fuel + oxidizer) + \ No newline at end of file diff --git a/samples/python/onedim/premixed_counterflow_twin_flame.py b/samples/python/onedim/premixed_counterflow_twin_flame.py index df3630a494..51f68fe041 100644 --- a/samples/python/onedim/premixed_counterflow_twin_flame.py +++ b/samples/python/onedim/premixed_counterflow_twin_flame.py @@ -2,14 +2,22 @@ Symmetric premixed twin flame ============================= -Simulate two counter-flow jets of reactants shooting into each other. This simulation -differs from the similar :doc:`premixed_counterflow_flame.py +Simulate two counter-flow jets of premixed reactants shooting into each other. This +simulation differs from the similar :doc:`premixed_counterflow_flame.py ` example as the latter simulates a jet of reactants shooting into products. -Requires: cantera >= 3.0 +Requires: cantera >= 3.0, matplotlib >= 2.0 .. tags:: Python, combustion, 1D flow, premixed flame, strained flame, plotting + +An illustration of this configuration is shown in the figure below: + +.. image:: /_static/images/samples/twin-premixed-flame.svg + :width: 90% + :alt: Twin flame established by two opposed jets of fuel + oxidizer mixtures + :align: center + """ import sys @@ -32,58 +40,61 @@ def derivative(x, y): return dydx -def computeStrainRates(oppFlame): +def compute_strain_rates(opposed_flame): # Compute the derivative of axial velocity to obtain normal strain rate - strainRates = derivative(oppFlame.grid, oppFlame.velocity) + strain_rates = derivative(opposed_flame.grid, opposed_flame.velocity) # Obtain the location of the max. strain rate upstream of the pre-heat zone. # This is the characteristic strain rate - maxStrLocation = abs(strainRates).argmax() - minVelocityPoint = oppFlame.velocity[:maxStrLocation].argmin() + max_strain_location = abs(strain_rates).argmax() + min_velocity_point = opposed_flame.velocity[:max_strain_location].argmin() # Characteristic Strain Rate = K - strainRatePoint = abs(strainRates[:minVelocityPoint]).argmax() - K = abs(strainRates[strainRatePoint]) + strain_rate_point = abs(strain_rates[:min_velocity_point]).argmax() + K = abs(strain_rates[strain_rate_point]) - return strainRates, strainRatePoint, K + return strain_rates, strain_rate_point, K -def computeConsumptionSpeed(oppFlame): +def compute_consumption_speed(opposed_flame): + Tb = max(opposed_flame.T) + Tu = min(opposed_flame.T) + rho_u = max(opposed_flame.density) - Tb = max(oppFlame.T) - Tu = min(oppFlame.T) - rho_u = max(oppFlame.density) + integrand = opposed_flame.heat_release_rate / opposed_flame.cp - integrand = oppFlame.heat_release_rate/oppFlame.cp - - total_heat_release = np.trapz(integrand, oppFlame.grid) - Sc = total_heat_release/(Tb - Tu)/rho_u + total_heat_release = np.trapz(integrand, opposed_flame.grid) + Sc = total_heat_release / (Tb - Tu) / rho_u return Sc # This function is called to run the solver -def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1, - ratio=2, slope=0.3, curve=0.3, prune=0.05): +def solve_opposed_flame(opposed_flame, mass_flux=0.12, loglevel=1, + ratio=2, slope=0.3, curve=0.3, prune=0.05): """ - Execute this function to run the Opposed Flow Simulation This function + Execute this function to run the Opposed Flow Simulation. This function takes a CounterFlowTwinPremixedFlame object as the first argument """ - oppFlame.reactants.mdot = massFlux - oppFlame.set_refine_criteria(ratio=ratio, slope=slope, curve=curve, prune=prune) + opposed_flame.reactants.mdot = mass_flux + opposed_flame.set_refine_criteria(ratio=ratio, slope=slope, + curve=curve, prune=prune) - oppFlame.show() - oppFlame.solve(loglevel, auto=True) + opposed_flame.show() + opposed_flame.solve(loglevel, auto=True) # Compute the strain rate, just before the flame. This is not necessarily # the maximum We use the max. strain rate just upstream of the pre-heat zone # as this is the strain rate that computations compare against, like when # plotting Su vs. K - strainRates, strainRatePoint, K = computeStrainRates(oppFlame) + strain_rates, strain_rate_point, K = compute_strain_rates(opposed_flame) - return np.max(oppFlame.T), K, strainRatePoint + return np.max(opposed_flame.T), K, strain_rate_point +# %% +# Define parameters +# ----------------- # Select the reaction mechanism gas = ct.Solution('gri30.yaml') @@ -99,27 +110,39 @@ def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1, # Domain half-width of 2.5 cm, meaning the whole domain is 5 cm wide width = 0.025 -# Done with initial conditions +# %% +# Set up the simulation +# --------------------- + # Compute the mass flux, as this is what the Flame object requires -massFlux = gas.density * axial_velocity # units kg/m2/s +mass_flux = gas.density * axial_velocity # units kg/m2/s # Create the flame object -oppFlame = ct.CounterflowTwinPremixedFlame(gas, width=width) +opposed_flame = ct.CounterflowTwinPremixedFlame(gas, width=width) # Uncomment the following line to use a Multi-component formulation. Default is # mixture-averaged -# oppFlame.transport_model = 'multicomponent' +# opposed_flame.transport_model = 'multicomponent' -# Now run the solver. The solver returns the peak temperature, strain rate and -# the point which we ascribe to the characteristic strain rate. +# %% +# Run the solver +# -------------- +# +# The solver returns the peak temperature, strain rate and the point which we ascribe to +# the characteristic strain rate. -(T, K, strainRatePoint) = solveOpposedFlame(oppFlame, massFlux, loglevel=1) +T, K, strain_rate_point = solve_opposed_flame(opposed_flame, mass_flux, loglevel=1) -# You can plot/see all state space variables by calling oppFlame.foo where foo -# is T, Y[i], etc. The spatial variable (distance in meters) is in oppFlame.grid -# Thus to plot temperature vs distance, use oppFlame.grid and oppFlame.T +# %% +# Compute flame properties +# ------------------------ +# +# You can plot/see all state space variables by calling `opposed_flame.foo`, where +# `foo` is `T`, `Y[i]`, and so forth. The spatial variable (distance in meters) is in +# `opposed_flame.grid`. Thus to plot temperature vs distance, use `opposed_flame.grid`` +# and `opposed_flame.T`. -Sc = computeConsumptionSpeed(oppFlame) +Sc = compute_consumption_speed(opposed_flame) if "native" in ct.hdf_support(): output = Path() / "premixed_counterflow_twin_flame.h5" @@ -127,14 +150,18 @@ def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1, output = Path() / "premixed_counterflow_twin_flame.yaml" output.unlink(missing_ok=True) -oppFlame.save(output, name="mix") +opposed_flame.save(output, name="mix") print(f"Peak temperature: {T:.1f} K") print(f"Strain Rate: {K:.1f} 1/s") print(f"Consumption Speed: {Sc * 100:.2f} cm/s") -oppFlame.save("premixed_counterflow_twin_flame.csv", basis="mole", overwrite=True) +opposed_flame.save("premixed_counterflow_twin_flame.csv", basis="mole", overwrite=True) -# Generate plots to see results, if user desires +# %% +# Plot results +# ------------ +# +# Generate plots to see results, if desired. if '--plot' in sys.argv: import matplotlib.pyplot as plt @@ -143,24 +170,24 @@ def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1, # Axial Velocity Plot plt.subplot(1, 2, 1) - plt.plot(oppFlame.grid, oppFlame.velocity, 'r', lw=2) - plt.xlim(oppFlame.grid[0], oppFlame.grid[-1]) + plt.plot(opposed_flame.grid, opposed_flame.velocity, 'r', lw=2) + plt.xlim(opposed_flame.grid[0], opposed_flame.grid[-1]) plt.xlabel('Distance (m)') plt.ylabel('Axial Velocity (m/s)') # Identify the point where the strain rate is calculated - plt.plot(oppFlame.grid[strainRatePoint], - oppFlame.velocity[strainRatePoint], 'gs') + plt.plot(opposed_flame.grid[strain_rate_point], + opposed_flame.velocity[strain_rate_point], 'gs') plt.annotate('Strain-Rate point', - xy=(oppFlame.grid[strainRatePoint], - oppFlame.velocity[strainRatePoint]), + xy=(opposed_flame.grid[strain_rate_point], + opposed_flame.velocity[strain_rate_point]), xytext=(0.001, 0.1), arrowprops={'arrowstyle': '->'}) # Temperature Plot plt.subplot(1, 2, 2) - plt.plot(oppFlame.grid, oppFlame.T, 'b', lw=2) - plt.xlim(oppFlame.grid[0], oppFlame.grid[-1]) + plt.plot(opposed_flame.grid, opposed_flame.T, 'b', lw=2) + plt.xlim(opposed_flame.grid[0], opposed_flame.grid[-1]) plt.xlabel('Distance (m)') plt.ylabel('Temperature (K)') From 53c47b46222a40d4c4d75169a5350fbdf2271d2c Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 25 Mar 2024 11:35:36 -0400 Subject: [PATCH 10/18] Migrate 1D packed bed example from Jupyter Co-authored-by: Gandhali Kogekar --- samples/python/reactors/1D_packed_bed.py | 463 +++++++++++++++++++++++ 1 file changed, 463 insertions(+) create mode 100644 samples/python/reactors/1D_packed_bed.py diff --git a/samples/python/reactors/1D_packed_bed.py b/samples/python/reactors/1D_packed_bed.py new file mode 100644 index 0000000000..96d1d8d752 --- /dev/null +++ b/samples/python/reactors/1D_packed_bed.py @@ -0,0 +1,463 @@ +""" +One-dimensional packed-bed, catalytic-membrane reactor +====================================================== + +The model shown in this example simulates heterogeneous catalytic processes inside +packed-bed, catalytic membrane reactors. The gas-phase and surface-phase species +conservation equations are derived and the system of differential-algebraic equations +(DAE) is solved using the ``scikits.odes.dae`` IDA solver. + +Requires: cantera >= 3.0.0, matplotlib >= 2.0, scikits.odes >= 2.7.0 + +.. tags:: Python, surface chemistry, 1D flow, catalysis, packed bed reactor, + porous media, user-defined model +""" + + +# %% +# Methodology +# ----------- +# +# A one-dimensional, steady-state catalytic-membrane reactor model with surface +# chemistry is developed to analyze species profiles along the length of a packed-bed, +# catalytic membrane reactor. The same model can further be simplified to simulate a +# simple packed-bed reactor by excluding the membrane. The example here demonstrates the +# one-dimensional reactor model explained by G. Kogekar [1]_. + +# %% +# Governing equations +# ^^^^^^^^^^^^^^^^^^^ +# +# Assuming steady-state, one-dimensional flow within the packed bed, total-mass, species +# mass and energy conservation may be stated as [2]_: +# +# .. math:: +# +# \frac{d(\rho u)}{dz} &= \sum_{k=1}^{K_\t{g}} \dot {s}_k W_k A_\t{s} +# + \frac{P_\t{b}}{A_\t{b}} j_{k_\t{M}}, \\ +# +# \rho u \frac{dY_k}{dz} + A_\t{s} Y_k \sum_{k=1}^{K_\t{g}} \dot {s}_k W_k &= +# A_\t{s} \dot {s}_k W_k + \delta_{k, k_\t{M}} \frac{P_\t{b}}{A_\t{b}} j_{k_\t{M}}, \\ +# +# \rho u c_\t{p} \frac{dT}{dz} +# + \sum_{k=1}^{K_\t{g}} h_k (\phi_\t{g} \dot {\omega}_k + A_\t{s} \dot {s}_k) W_k +# &= \hat h \frac{P_\t{b}}{A_\t{b}}(T_\t{w} - T) +# + \delta_{k, k_\t{M}} \frac{P_\t{b}}{A_\t{b}} h_{k_\t{M}} j_{k_\t{M}}. +# +# The fractional coverages of the :math:`K_\t{s}` surface adsorbates :math:`\theta_k` +# must satisfy +# +# .. math:: +# \dot {s}_k = 0, \qquad (k = 1,\ldots, K_\t{s}), +# +# which, at steady state, requires no net production/consumption of surface species by +# the heterogeneous reactions [3]_. +# +# The pressure within the bed is calculated as: +# +# .. math:: +# \frac{dp}{dz} = - \left( \frac{\phi_\t{g} \mu}{\beta_\t{g}} \right) u, +# +# where :math:`\mu` is the gas-phase mixture viscosity. The packed-bed permeability +# :math:\beta_\t{g}` is evaluated using the Kozeny-Carman relationship as +# +# .. math:: +# \beta_\t{g} = \frac{\phi_\t{g}^3 D_\t{p}^2}{72 \tau_\t{g}(1 - \phi_\t{g})^2}, +# +# where :math:`\phi_\t{g}`, :math:`\tau_\t{g}`, and :math:`D_\t{p}` are the bed +# porosity, tortuosity, and particle diameter, respectively. +# +# The independent variable in these conservation equations is the position :math:`z` +# along the reactor length. The dependent variables include total mass flux :math:`\rho +# u`, pressure :math:`p`, temperature :math:`T`, gas-phase mass fractions :math:`Y_k`, +# and surfaces coverages :math:`\theta_k`. Gas-phase fluxes through the membrane are +# represented as :math:`j_{k, \t{M}}`. Geometric parameters :math:`A_\t{s}`, +# :math:`P_\t{b}`, and :math:`A_\t{b}` represent the catalyst specific surface area +# (dimensions of surface area per unit volume), reactor perimeter, and reactor +# cross-sectional flow area, respectively. Other parameters include bed porosity +# :math:`\phi_\t{g}` and gas-phase species molecular weights :math:`W_k`. The gas +# density :math:`\rho` is evaluated using the equation of state (such as ideal gas, +# Redlich-Kwong or Peng-Robinson). +# +# If a perm-selective membrane is present, then :math:`j_{k_\t{M}}` represents the +# gas-phase flux through the membrane and :math:`{k_\t{M}}` is the gas-phase species +# that permeates through the membrane. The Kronecker delta, :math:`\delta_{k, k_\t{M}} = +# 1` for the membrane-permeable species and :math:`\delta_{k, k_\t{M}} = 0` otherwise. +# The membrane flux is calculated using Sievert's law as +# +# .. math:: +# j_{k_\t{M}}^{\text{Mem}} = +# \frac{B_{k_\t{M}}}{t} \left( p_{k_\t{M} {\text{, mem}}}^\alpha - p_{k_\t{M} \text{, sweep}}^\alpha \right) W_{k_\t{M}} +# +# where :math:`B_{k_\t{M}}` is the membrane permeability, :math:`t` is the membrane +# thickness. :math:`p_{k_\t{M} \text{, mem}}` and :math:`p_{k_\t{M} \text{, sweep}}` +# represent perm-selective species partial pressures within the packed-bed and the +# exterior sweep channel. The present model takes the pressure exponent :math:`\alpha` +# to be unity. The membrane flux for all other species (:math:`k \neq k_\t{M}`) is zero. +# +# Chemistry mechanism +# ^^^^^^^^^^^^^^^^^^^ +# +# This example uses a detailed 12-step elementary micro-kinetic reaction mechanism that +# describes ammonia formation and decomposition kinetics over the Ru/Ba-YSZ catalyst. +# The reaction mechanism is developed and validated using measured performance in a +# laboratory-scale packed-bed reactor [4]_. This example also incorporates a Pd-based +# H₂ perm-selective membrane. +# +# Solver +# ^^^^^^ +# +# The above governing equations represent a complete solution for a steady-state +# packed-bed, catalytic membrane reactor model. The dependent variables are the +# mass-flux :math:`\rho u`, species mass-fractions :math:`Y_k`, pressure :math:`p`, +# temperature :math:`T`, and surface coverages :math:`\theta_k`. The equation of state +# is used to obtain the mass density, :math:`\rho`. +# +# The governing equations form an initial value problem (IVP) in a +# differential-algebraic (DAE) form as follows: +# +# .. math:: +# f(z,{\bf{y}}, {\bf y'}, {\bf c}) = {\bf 0}, +# +# where :math:`\bf y` and :math:`\bf y'` represent the state vector and its derivative +# vector, respectively. All other constants such as reference temperature, chemical +# constants, etc. are incorporated in vector :math:`\bf c` (Refer to [1]_ for more +# details). This type of DAE system in this example is solved using the +# ``scikits.odes.dae`` IDA solver. + +# %% +# Import Cantera and scikits +# -------------------------- + +import numpy as np +from scikits.odes import dae +import cantera as ct +import matplotlib.pyplot as plt + +# %% +# Define gas-phase and surface-phase species +# ------------------------------------------ + +# Import the reaction mechanism for Ammonia synthesis/decomposition on Ru-Ba/YSZ catalyst +mechfile = "example_data/ammonia-Ru-Ba-YSZ-CSM-2019.yaml" +# Import the models for gas-phase +gas = ct.Solution(mechfile, "gas") +# Import the model for surface-phase +surf = ct.Interface(mechfile, "Ru_surface", [gas]) + +# Other parameters +n_gas = gas.n_species # number of gas species +n_surf = surf.n_species # number of surface species +n_gas_reactions = gas.n_reactions # number of gas-phase reactions + +# Set offsets of dependent variables in the solution vector +offset_rhou = 0 +offset_p = 1 +offset_T = 2 +offset_Y = 3 +offset_Z = offset_Y + n_gas +n_var = offset_Z + n_surf # total number of variables (rhou, P, T, Yk and Zk) + +print("Number of gas-phase species = ", n_gas) +print("Number of surface-phase species = ", n_surf) +print("Number of variables = ", n_var) + +# %% +# Define reactor geometry and operating conditions +# ------------------------------------------------ + +# Reactor geometry +L = 5e-2 # length of the reactor (m) +R = 5e-3 # radius of the reactor channel (m) +phi = 0.5 # porosity of the bed (-) +tau = 2.0 # tortuosity of the bed (-) +dp = 3.37e-4 # particle diameter (m) +As = 3.5e6 # specific surface area (1/m) + +# Energy (adiabatic or isothermal) +solve_energy = True # True: Adiabatic, False: isothermal + +# Membrane (True: membrane, False: no membrane) +membrane_present = True +membrane_perm = 1e-15 # membrane permeability (kmol*m3/s/Pa) +thickness = 3e-6 # membrane thickness (m) +membrane_sp_name = "H2" # membrane-permeable species name +p_sweep = 1e5 # partial pressure of permeable species in the sweep channel (Pa) +permeance = membrane_perm / thickness # permeance of the membrane (kmol*m2/s/Pa) + +if membrane_present: + print("Modeling packed-bed, catalytic-membrane reactor...") + print(membrane_sp_name, "permeable membrane is present.") + +# Get required properties based on the geometry and mechanism +W_g = gas.molecular_weights # vector of molecular weight of gas species +vol_ratio = phi / (1 - phi) +eff_factor = phi / tau # effective factor for permeability calculation +# permeability based on Kozeny-Carman equation +B_g = B_g = vol_ratio**2 * dp**2 * eff_factor / 72 +area2vol = 2 / R # area to volume ratio assuming a cylindrical reactor +D_h = 2 * R # hydraulic diameter +membrane_sp_ind = gas.species_index(membrane_sp_name) + +# Inlet operating conditions +T_in = 673 # inlet temperature [K] +p_in = 5e5 # inlet pressure [Pa] +v_in = 0.001 # inlet velocity [m/s] +T_wall = 723 # wall temperature [K] +h_coeff = 1e2 # heat transfer coefficient [W/m2/K] + +# Set gas and surface states +gas.TPX = T_in, p_in, "NH3:0.99, AR:0.01" # inlet composition +surf.TP = T_in, p_in +Yk_0 = gas.Y +rhou0 = gas.density * v_in + +# Initial surface coverages +# advancing coverages over a long period of time to get the steady state. +surf.advance_coverages(1e10) +Zk_0 = surf.coverages + + +# %% +# Define residual function required for IDA solver +# ------------------------------------------------ + +def residual(z, y, yPrime, res): + """Solution vector for the model + y = [rho*u, p, T, Yk, Zk] + yPrime = [d(rho*u)dz, dpdz, dTdz, dYkdz, dZkdz] + """ + # Get current thermodynamic state from solution vector and save it to local variables. + rhou = y[offset_rhou] # mass flux (density * velocity) + Y = y[offset_Y : offset_Y + n_gas] # vector of mass fractions + Z = y[offset_Z : offset_Z + n_surf] # vector of site fractions + p = y[offset_p] # pressure + T = y[offset_T] # temperature + + # Get derivatives of dependent variables + drhoudz = yPrime[offset_rhou] + dYdz = yPrime[offset_Y : offset_Y + n_gas] + dZdz = yPrime[offset_Z : offset_Z + n_surf] + dpdz = yPrime[offset_p] + dTdz = yPrime[offset_T] + + # Set current thermodynamic state for the gas and surface phases + # Note: use unnormalized mass fractions and site fractions to avoid + # over-constraining the system + gas.set_unnormalized_mass_fractions(Y) + gas.TP = T, p + surf.set_unnormalized_coverages(Z) + surf.TP = T, p + + # Calculate required variables based on the current state + coverages = surf.coverages # surface site coverages + # heterogeneous production rate of gas species + sdot_g = surf.get_net_production_rates("gas") + # heterogeneous production rate of surface species + sdot_s = surf.get_net_production_rates("Ru_surface") + wdot_g = np.zeros(n_gas) + # specific heat of the mixture + cp = gas.cp_mass + # partial enthalpies of gas species + hk_g = gas.partial_molar_enthalpies + + if n_gas_reactions > 0: + # homogeneous production rate of gas species + wdot_g = gas.net_production_rates + mu = gas.viscosity # viscosity of the gas-phase + + # Calculate density using equation of state + rho = gas.density + + # Calculate flux term through the membrane + # partial pressure of membrane-permeable species + memsp_pres = p * gas.X[membrane_sp_ind] + # negative sign indicates the flux going out + membrane_flux = -permeance * (memsp_pres - p_sweep) * W_g[membrane_sp_ind] + + # Conservation of total-mass + # temporary variable + sum_continuity = As * np.sum(sdot_g * W_g) + phi * np.sum(wdot_g * W_g) + res[offset_rhou] = (drhoudz - sum_continuity + - area2vol * membrane_flux * membrane_present) + + # Conservation of gas-phase species + res[offset_Y:offset_Y+ n_gas] = (dYdz + (Y * sum_continuity + - phi * np.multiply(wdot_g,W_g) + - As * np.multiply(sdot_g,W_g)) / rhou) + res[offset_Y + membrane_sp_ind] -= area2vol * membrane_flux * membrane_present + + # Conservation of site fractions (algebraic constraints in this example) + res[offset_Z : offset_Z + n_surf] = sdot_s + + # For the species with largest site coverage (k_large), solve the constraint + # equation: sum(Zk) = 1. + # The residual function for 'k_large' would be 'res[k_large] = 1 - sum(Zk)' + # Note that here sum(Zk) will be the sum of coverages for all surface species, + # including the 'k_large' species. + ind_large = np.argmax(coverages) + res[offset_Z + ind_large] = 1 - np.sum(coverages) + + # Conservation of momentum + u = rhou / rho + res[offset_p] = dpdz + phi * mu * u / B_g + + # Conservation of energy + res[offset_T] = dTdz - 0 # isothermal condition + # Note: One can just not solve the energy equation by keeping temperature constant. + # But since 'T' is used as the dependent variable, the residual is res[T] = dTdz - 0 + # One can also write res[T] = 0 directly, but that leads to a solver failure due to + # singular jacobian + + if solve_energy: + conv_term = (4 / D_h) * h_coeff * (T_wall - T) * (2 * np.pi * R) + chem_term = np.sum(hk_g * (phi * wdot_g + As * sdot_g)) + res[offset_T] -= (conv_term - chem_term) / (rhou * cp) + +# %% +# Calculate the spatial derivatives at the inlet that will be used as the initial +# conditions for the IDA solver + +# Initialize yPrime to 0 and call residual to get initial derivatives +y0 = np.hstack((rhou0, p_in, T_in, Yk_0, Zk_0)) +yprime0 = np.zeros(n_var) +res = np.zeros(n_var) +residual(0, y0, yprime0, res) +yprime0 = -res + +# %% +# Solve the system of DAEs using IDA solver +# ----------------------------------------- + +solver = dae( + "ida", + residual, + first_step_size=1e-15, + atol=1e-14, # absolute tolerance for solution + rtol=1e-06, # relative tolerance for solution + algebraic_vars_idx=[np.arange(offset_Y + n_gas, offset_Z + n_surf, 1)], + max_steps=8000, + one_step_compute=True, + old_api=False, # forces use of new api (namedtuple) +) + +distance = [] +solution = [] +state = solver.init_step(0.0, y0, yprime0) + +# Note that here the variable t is an internal variable used in scikits. In this +# example, it represents the natural variable z, which corresponds to the axial distance +# inside the reactor. +while state.values.t < L: + distance.append(state.values.t) + solution.append(state.values.y) + state = solver.step(L) + +distance = np.array(distance) +solution = np.array(solution) +print(state) + +# %% +# Plot results +# ------------ + +plt.rcParams['figure.constrained_layout.use'] = True + +# %% +# Pressure and temperature +# ^^^^^^^^^^^^^^^^^^^^^^^^ + +# Plot gas pressure profile along the flow direction +f, ax = plt.subplots(figsize=(4,3)) +ax.plot(distance, solution[:, offset_p], color="C0") +ax.set_xlabel("Distance (m)") +ax.set_ylabel("Pressure (Pa)") + +# Plot gas temperature profile along the flow direction +f, ax = plt.subplots(figsize=(4,3)) +ax.plot(distance, solution[:, offset_T], color="C1") +ax.set_xlabel("Distance (m)") +ax.set_ylabel("Temperature (K)") + +# %% +# Mass fractions of the gas-phase species +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +# Plot major and minor gas species separately +minor_idx = [] +major_idx = [] +for j, name in enumerate(gas.species_names): + mean = np.mean(solution[:, offset_Y + j]) + if mean <= 0.1: + minor_idx.append(j) + else: + major_idx.append(j) + +# Major gas-phase species +f, ax = plt.subplots(figsize=(4,3)) +for j in major_idx: + ax.plot(distance, solution[:, offset_Y + j], label=gas.species_names[j]) +ax.legend(fontsize=12, loc="best") +ax.set_xlabel("Distance (m)") +ax.set_ylabel("Mass Fraction") + +# Minor gas-phase species +f, ax = plt.subplots(figsize=(4,3)) +for j in minor_idx: + ax.plot(distance, solution[:, offset_Y + j], label=gas.species_names[j]) +ax.legend(fontsize=12, loc="best") +ax.set_xlabel("Distance (m)") +ax.set_ylabel("Mass Fraction") + +# %% +# Site fractions of the surface-phase species +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +# Plot major and minor surface species separately +minor_idx = [] +major_idx = [] +for j, name in enumerate(surf.species_names): + mean = np.mean(solution[:, offset_Z + j]) + if mean <= 0.1: + minor_idx.append(j) + else: + major_idx.append(j) + +# Major surf-phase species +f, ax = plt.subplots(figsize=(4,3)) +for j in major_idx: + ax.plot(distance, solution[:, offset_Z + j], label=surf.species_names[j]) +ax.legend(fontsize=12, loc="best") +ax.set_xlabel("Distance (m)") +ax.set_ylabel("Site Fraction") + +# Minor surf-phase species +f, ax = plt.subplots(figsize=(4,3)) +for j in minor_idx: + ax.plot(distance, solution[:, offset_Z + j], label=surf.species_names[j]) +ax.legend(fontsize=12, loc="best") +ax.set_xlabel("Distance (m)") +ax.set_ylabel("Site Fraction") + +#%% +# References +# ---------- +# +# .. [1] G. Kogekar (2021). "Computationally efficient and robust models of non-ideal +# thermodynamics, gas-phase kinetics and heterogeneous catalysis in chemical +# reactors", Doctoral dissertation. +# +# .. [2] B. Kee, C. Karakaya, H. Zhu, S. DeCaluwe, and R.J. Kee (2017). "The Influence +# of Hydrogen-Permeable Membranes and Pressure on Methane Dehydroaromatization in +# Packed-Bed Catalytic Reactors", *Industrial & Engineering Chemistry Research* +# 56, 13:3551-3559. +# +# .. [3] R.J. Kee, M.E. Coltrin, P. Glarborg, and H. Zhu (2018). *Chemically Reacting +# Flow: Theory, Modeling and Simulation*, Wiley. +# +# .. [4] Z. Zhang, C. Karakaya, R.J. Kee, J. Douglas Way, C. Wolden (2019). +# "Barium-Promoted Ruthenium Catalysts on Yittria-Stabilized Zirconia Supports +# for Ammonia Synthesis", *ACS Sustainable Chemistry & Engineering* +# 7:18038-18047. + +# sphinx_gallery_thumbnail_number = -1 From c1cfe2619c4e92b8d1470cc79434b9dea0f6f134 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 26 Mar 2024 10:58:37 -0400 Subject: [PATCH 11/18] Revise nonideal shock tube example based on Jupyter Notebooks Brings in some updates based on nonideal_shock_tube.ipynb and batch_reactor_ignition_delay_NTC.ipynb. Co-authored-by: Santosh Shanbhogue Co-authored-by: Steven DeCaluwe --- samples/python/reactors/NonIdealShockTube.py | 301 +++++++++++-------- 1 file changed, 172 insertions(+), 129 deletions(-) diff --git a/samples/python/reactors/NonIdealShockTube.py b/samples/python/reactors/NonIdealShockTube.py index 9227b8a84a..bf3e684b65 100644 --- a/samples/python/reactors/NonIdealShockTube.py +++ b/samples/python/reactors/NonIdealShockTube.py @@ -2,27 +2,75 @@ Ignition delay time using the Redlich-Kwong real gas model ========================================================== -Ignition delay time computations in a high-pressure reflected shock tube -reactor, comparing ideal gas and Redlich-Kwong real gas models. - In this example we illustrate how to setup and use a constant volume, adiabatic reactor to simulate reflected shock tube experiments. This reactor -will then be used to compute the ignition delay of a gas at a specified -initial temperature and pressure. The example is written in a general way, -that is, no particular EoS is presumed and ideal and real gas EoS can be used -equally easily. - -The reactor (system) is simply an 'insulated box,' and can technically be used -for any number of equations of state and constant-volume, adiabatic reactors. +is then used to compute the ignition delay of a gas at a specified +initial temperature and pressure. -Other than the typical Cantera dependencies, plotting functions require that -you have matplotlib (https://matplotlib.org/) installed. +The example is written in a general way, that is, no particular equation of state (EoS) +is presumed and ideal and real gas EoS can be used equally easily. The example here +demonstrates the calculations carried out by G. Kogekar et al. [1]_. Requires: cantera >= 2.5.0, matplotlib >= 2.0 .. tags:: Python, combustion, reactor network, non-ideal fluid, ignition delay, plotting """ +# %% +# Methods +# ------- +# +# Reflected Shock Reactor +# ^^^^^^^^^^^^^^^^^^^^^^^ +# +# The reflected shock tube reactor is modeled as a closed, constant-volume, adiabatic +# reactor. The heat transfer and the work rates are therefore both zero. With no mass +# inlets or exits, the 1st law energy balance reduces to: +# +# .. math:: +# \frac{dU}{dt} = \dot{Q} - \dot{W} = 0 +# +# Because of the constant-mass and constant-volume assumptions, the density is also +# therefore constant: +# +# .. math:: +# \frac{d\rho}{dt} = 0 +# +# Along with the evolving gas composition, the thermodynamic state of the gas is defined +# by the initial total internal energy :math:`U = m u = m \sum_k Y_k u_k`, where +# :math:`u_k` and :math:`Y_k` are the specific internal energy (J/kg) and mass fraction +# of species :math:`k`, respectively. +# +# The species mass fractions evolve according to the net chemical production rates due +# to homogeneous gas-phase reactions: +# +# .. math:: +# \frac{dY_k}{dt} = \frac{W_k}{\rho}\dot{\omega}_k, +# +# where :math:`W_k` is the molecular weight of species :math:`k` (kg/kmol³), +# :math:`\rho` is the (constant) gas-phase density (kg/m³), and :math:`\dot{\omega}_k` +# is the net production rate of species :math:`k` (kmol/m³/s). +# +# Redlich-Kwong Parameters +# ^^^^^^^^^^^^^^^^^^^^^^^^ +# +# Redlich-Kwong constants for each species are calculated according to their critical +# temperature :math:`T_c` and pressure :math:`P_c`: +# +# .. math:: +# a = 0.4275 \frac{R^2 T_c^{2.5}}{P_c} +# +# and +# +# .. math:: +# b = 0.08664 \frac{R T_c}{P_c} +# +# where :math:`R` is the universal gas constant. +# +# For stable species, the critical properties are readily available. For radicals and +# other short-lived intermediates, the Joback method [3]_ is used to estimate critical +# properties. + # Dependencies: numpy, and matplotlib import numpy as np import matplotlib.pyplot as plt @@ -32,51 +80,36 @@ import cantera as ct print('Running Cantera version: ' + ct.__version__) - +# %% # Define the ignition delay time (IDT). This function computes the ignition -# delay from the occurrence of the peak concentration for the specified -# species. -def ignitionDelay(states, species): +# delay from the occurrence of the peak concentration for the specified species. +def ignition_delay(states, species): i_ign = states(species).Y.argmax() return states.t[i_ign] +# %% +# Define initial conditions and reaction mechanism +# ------------------------------------------------ +# +# In this example we will choose a stoichiometric mixture of n-dodecane and air as the +# gas. For a representative kinetic model, we use the one from Wang et al [2]_. # Define the reactor temperature and pressure: -reactorTemperature = 1000 # Kelvin -reactorPressure = 40.0*101325.0 # Pascals - -# Define the gas: In this example we will choose a stoichiometric mixture of -# n-dodecane and air as the gas. For a representative kinetic model, we use: -# -# H.Wang, Y.Ra, M.Jia, R.Reitz, Development of a reduced n-dodecane-PAH -# mechanism. and its application for n-dodecane soot predictions., Fuel 136 -# (2014) 25–36. doi:10.1016/j.fuel.2014.07.028 +reactor_temperature = 1000 # Kelvin +reactor_pressure = 40.0*101325.0 # Pascals -# R-K constants are calculated according to their critical temperature (Tc) and -# pressure (Pc): -# -# a = 0.4275*(R^2)*(Tc^2.5)/(Pc) -# -# and -# -# b = 0.08664*R*Tc/Pc -# -# where R is the gas constant. -# -# For stable species, the critical properties are readily available. For -# radicals and other short-lived intermediates, the Joback method is used to -# estimate critical properties. For details of the method, see: Joback and Reid, -# "Estimation of pure- component properties from group-contributions," Chem. -# Eng. Comm. 57 (1987) 233-243, doi: 10.1080/00986448708960487 +# Load the real gas mechanism: +real_gas = ct.Solution('nDodecane_Reitz.yaml', 'nDodecane_RK') +# Create the ideal gas object: +ideal_gas = ct.Solution('nDodecane_Reitz.yaml', 'nDodecane_IG') +# %% # Real gas IDT calculation - -# Load the real gas mechanism: -real_gas = ct.Solution('nDodecane_Reitz.yaml', 'nDodecane_RK') +# ------------------------ # Set the state of the gas object: -real_gas.TP = reactorTemperature, reactorPressure +real_gas.TP = reactor_temperature, reactor_pressure # Define the fuel, oxidizer and set the stoichiometry: real_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26', @@ -85,86 +118,75 @@ def ignitionDelay(states, species): # Create a reactor object and add it to a reactor network # In this example, this will be the only reactor in the network r = ct.Reactor(contents=real_gas) -reactorNetwork = ct.ReactorNet([r]) -timeHistory_RG = ct.SolutionArray(real_gas, extra=['t']) -# Tic -t0 = time.time() +reactor_network = ct.ReactorNet([r]) +time_history_RG = ct.SolutionArray(real_gas, extra=['t']) # This is a starting estimate. If you do not get an ignition within this time, # increase it -estimatedIgnitionDelayTime = 0.005 -t = 0 +estimated_ignition_delay_time = 0.005 +t = 0 +t0 = time.time() counter = 1 -while t < estimatedIgnitionDelayTime: - t = reactorNetwork.step() +while t < estimated_ignition_delay_time: + t = reactor_network.step() if counter % 20 == 0: # We will save only every 20th value. Otherwise, this takes too long # Note that the species concentrations are mass fractions - timeHistory_RG.append(r.thermo.state, t=t) + time_history_RG.append(r.thermo.state, t=t) counter += 1 # We will use the 'oh' species to compute the ignition delay -tau_RG = ignitionDelay(timeHistory_RG, 'oh') +tau_RG = ignition_delay(time_history_RG, 'oh') -# Toc t1 = time.time() print("Computed Real Gas Ignition Delay: {:.3e} seconds. " - "Took {:3.2f}s to compute".format(tau_RG, t1-t0)) - + "Took {:3.2f} s to compute".format(tau_RG, t1-t0)) +# %% # Ideal gas IDT calculation -# Create the ideal gas object: -ideal_gas = ct.Solution('nDodecane_Reitz.yaml', 'nDodecane_IG') +# ------------------------- # Set the state of the gas object: -ideal_gas.TP = reactorTemperature, reactorPressure +ideal_gas.TP = reactor_temperature, reactor_pressure # Define the fuel, oxidizer and set the stoichiometry: ideal_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26', oxidizer={'o2': 1.0, 'n2': 3.76}) r = ct.Reactor(contents=ideal_gas) -reactorNetwork = ct.ReactorNet([r]) -timeHistory_IG = ct.SolutionArray(ideal_gas, extra=['t']) +reactor_network = ct.ReactorNet([r]) +time_history_IG = ct.SolutionArray(ideal_gas, extra=['t']) -# Tic t0 = time.time() - t = 0 - counter = 1 -while t < estimatedIgnitionDelayTime: - t = reactorNetwork.step() +while t < estimated_ignition_delay_time: + t = reactor_network.step() if counter % 20 == 0: # We will save only every 20th value. Otherwise, this takes too long # Note that the species concentrations are mass fractions - timeHistory_IG.append(r.thermo.state, t=t) + time_history_IG.append(r.thermo.state, t=t) counter += 1 # We will use the 'oh' species to compute the ignition delay -tau_IG = ignitionDelay(timeHistory_IG, 'oh') - -# Toc +tau_IG = ignition_delay(time_history_IG, 'oh') t1 = time.time() print("Computed Ideal Gas Ignition Delay: {:.3e} seconds. " - "Took {:3.2f}s to compute".format(tau_IG, t1-t0)) + "Took {:3.2f} s to compute".format(tau_IG, t1-t0)) print('Ideal gas error: {:2.2f} %'.format(100*(tau_IG-tau_RG)/tau_RG)) -# Plot the result +# %% +# Plot the results +# ---------------- -plt.rcParams['xtick.labelsize'] = 12 -plt.rcParams['ytick.labelsize'] = 12 -plt.rcParams['figure.autolayout'] = True -plt.rcParams['axes.labelsize'] = 14 -plt.rcParams['font.family'] = 'serif' +plt.rcParams['figure.constrained_layout.use'] = True # Figure illustrating the definition of ignition delay time (IDT). - plt.figure() -plt.plot(timeHistory_RG.t, timeHistory_RG('oh').Y, '-o', color='b', markersize=4) -plt.plot(timeHistory_IG.t, timeHistory_IG('oh').Y, '-o', color='r', markersize=4) +plt.plot(time_history_RG.t, time_history_RG('oh').Y) +plt.plot(time_history_IG.t, time_history_IG('oh').Y, '-.') plt.xlabel('Time (s)') plt.ylabel(r'OH mass fraction, $\mathdefault{Y_{OH}}$') @@ -172,107 +194,109 @@ def ignitionDelay(states, species): plt.xlim([0, 0.00055]) ax = plt.gca() ax.annotate("", xy=(tau_RG, 0.005), xytext=(0, 0.005), - arrowprops=dict(arrowstyle="<|-|>", color='k', linewidth=2.0), - fontsize=14) -plt.annotate('Ignition Delay Time (IDT)', xy=(0, 0), xytext=(0.00008, 0.00525), - fontsize=16) + arrowprops=dict(arrowstyle="<|-|>", color='k', linewidth=2.0)) +plt.annotate('Ignition Delay Time (IDT)', xy=(0, 0), xytext=(0.00008, 0.00525)) -plt.legend(['Real Gas', 'Ideal Gas'], frameon=False) +plt.legend(['Real Gas', 'Ideal Gas']) # If you want to save the plot, uncomment this line (and edit as you see fit): # plt.savefig('IDT_nDodecane_1000K_40atm.pdf', dpi=350, format='pdf') - +# %% # Demonstration of NTC behavior -# Let us use the reactor model to demonstrate the impacts of non-ideal behavior on IDTs -# in the Negative Temperature Coefficient (NTC) region, where observed IDTs, counter -# to intuition, increase with increasing temperature. +# ----------------------------- +# +# A common benchmark for a reaction mechanism is its ability to reproduce the negative +# temperature coefficient (NTC) behavior. Intuitively, as the temperature of an explosive +# mixture increases, it should ignite faster. But, under certain conditions, we observe +# the opposite. This is referred to as NTC behavior. Reproducing experimentally observed +# NTC behavior is thus an important test for any mechanism. We will do this now by +# computing and visualizing the ignition delay for a wide range of temperatures # Make a list of all the temperatures at which we would like to run simulations: -T = np.array([1250, 1225, 1200, 1150, 1100, 1075, 1050, 1025, 1012.5, 1000, 987.5, - 975, 962.5, 950, 937.5, 925, 912.5, 900, 875, 850, 825, 800]) +T = np.array([1250, 1170, 1120, 1080, 1040, 1010, 990, 970, 950, 930, 910, 880, 850, + 820, 790, 760]) # If we desire, we can define different IDT starting guesses for each temperature: -estimatedIgnitionDelayTimes = np.ones(len(T)) +estimated_ignition_delay_times = np.ones(len(T)) # But we won't, at least in this example :) -estimatedIgnitionDelayTimes[:] = 0.005 +estimated_ignition_delay_times[:] = 0.005 # Now, we simply run the code above for each temperature. # Real Gas -ignitionDelays_RG = np.zeros(len(T)) +ignition_delays_RG = np.zeros(len(T)) for i, temperature in enumerate(T): # Setup the gas and reactor - reactorTemperature = temperature - real_gas.TP = reactorTemperature, reactorPressure + reactor_temperature = temperature + real_gas.TP = reactor_temperature, reactor_pressure real_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26', oxidizer={'o2': 1.0, 'n2': 3.76}) r = ct.Reactor(contents=real_gas) - reactorNetwork = ct.ReactorNet([r]) + reactor_network = ct.ReactorNet([r]) # create an array of solution states - timeHistory = ct.SolutionArray(real_gas, extra=['t']) + time_history = ct.SolutionArray(real_gas, extra=['t']) t0 = time.time() - t = 0 counter = 1 - while t < estimatedIgnitionDelayTimes[i]: - t = reactorNetwork.step() + while t < estimated_ignition_delay_times[i]: + t = reactor_network.step() if counter % 20 == 0: - timeHistory.append(r.thermo.state, t=t) + time_history.append(r.thermo.state, t=t) counter += 1 - tau = ignitionDelay(timeHistory, 'oh') + tau = ignition_delay(time_history, 'oh') t1 = time.time() - print("Computed Real Gas Ignition Delay: {:.3e} seconds for T={}K. " - "Took {:3.2f}s to compute".format(tau, temperature, t1-t0)) + print("Computed Real Gas Ignition Delay: {:.3e} seconds for T={:.1f} K. " + "Took {:3.2f} s to compute".format(tau, temperature, t1-t0)) - ignitionDelays_RG[i] = tau + ignition_delays_RG[i] = tau # Repeat for Ideal Gas -ignitionDelays_IG = np.zeros(len(T)) +ignition_delays_IG = np.zeros(len(T)) for i, temperature in enumerate(T): # Setup the gas and reactor - reactorTemperature = temperature - ideal_gas.TP = reactorTemperature, reactorPressure + reactor_temperature = temperature + ideal_gas.TP = reactor_temperature, reactor_pressure ideal_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26', oxidizer={'o2': 1.0, 'n2': 3.76}) r = ct.Reactor(contents=ideal_gas) - reactorNetwork = ct.ReactorNet([r]) + reactor_network = ct.ReactorNet([r]) # create an array of solution states - timeHistory = ct.SolutionArray(ideal_gas, extra=['t']) + time_history = ct.SolutionArray(ideal_gas, extra=['t']) t0 = time.time() t = 0 counter = 1 - while t < estimatedIgnitionDelayTimes[i]: - t = reactorNetwork.step() + while t < estimated_ignition_delay_times[i]: + t = reactor_network.step() if counter % 20 == 0: - timeHistory.append(r.thermo.state, t=t) + time_history.append(r.thermo.state, t=t) counter += 1 - tau = ignitionDelay(timeHistory, 'oh') + tau = ignition_delay(time_history, 'oh') t1 = time.time() - print("Computed Ideal Gas Ignition Delay: {:.3e} seconds for T={}K. " - "Took {:3.2f}s to compute".format(tau, temperature, t1-t0)) + print("Computed Ideal Gas Ignition Delay: {:.3e} seconds for T={:.1f} K. " + "Took {:3.2f} s to compute".format(tau, temperature, t1-t0)) - ignitionDelays_IG[i] = tau + ignition_delays_IG[i] = tau +# %% # Figure: ignition delay (tau) vs. the inverse of temperature (1000/T). -fig = plt.figure() -ax = fig.add_subplot(111) -ax.plot(1000/T, 1e6*ignitionDelays_RG, '-', linewidth=2.0, color='b') -ax.plot(1000/T, 1e6*ignitionDelays_IG, '-.', linewidth=2.0, color='r') -ax.set_ylabel(r'Ignition Delay ($\mathdefault{\mu s}$)', fontsize=14) -ax.set_xlabel(r'1000/T (K$^\mathdefault{-1}$)', fontsize=14) +fig, ax = plt.subplots() +ax.plot(1000/T, 1e6*ignition_delays_RG, '.-', linewidth=2.0) +ax.plot(1000/T, 1e6*ignition_delays_IG, '-.', linewidth=2.0) +ax.set_ylabel(r'Ignition Delay (μs)') +ax.set_xlabel(r'1000/T (K$^{-1}$)') -ax.set_xlim([0.8, 1.2]) +ax.set_xlim([0.8, 1.3]) # Add a second axis on top to plot the temperature for better readability ax2 = ax.twiny() @@ -280,12 +304,31 @@ def ignitionDelay(states, species): ax2.set_xticks(ticks) ax2.set_xticklabels((1000/ticks).round(1)) ax2.set_xlim(ax.get_xlim()) -ax2.set_xlabel('Temperature (K)', fontsize=14) +ax2.set_xlabel('Temperature (K)') -ax.legend(['Real Gas', 'Ideal Gas'], frameon=False, loc='upper left') +ax.legend(['Real Gas', 'Ideal Gas'], loc='upper left') # If you want to save the plot, uncomment this line (and edit as you see fit): # plt.savefig('NTC_nDodecane_40atm.pdf', dpi=350, format='pdf') # Show the plots. plt.show() + +# %% +# References +# ---------- +# +# .. [1] G. Kogekar, C. Karakaya, G. J. Liskovich, M. A. Oehlschlaeger, S. C. DeCaluwe, +# R. J. Kee (2018). "Impact of non-ideal behavior on ignition delay and chemical +# kinetics in high-pressure shock tube reactors," *Combust. Flame.* 189 1-11, +# https://doi.org/10.1016/j.combustflame.2017.10.014 +# +# .. [2] H. Wang, Y. Ra, M. Jia, R. Reitz (2014). "Development of a reduced +# n-dodecane-PAH mechanism and its application for n-dodecane soot predictions", +# *Fuel* 136, 25–36. https://doi.org/10.1016/j.fuel.2014.07.028. +# +# .. [3] K. G. Joback and R. C. Reid (1987). "Estimation of pure-component properties +# from group-contributions," *Chem. Eng. Comm.* 57, 233-243, +# https://doi.org/10.1080/00986448708960487. + +# sphinx_gallery_thumbnail_number = -1 From 1fb7b2927073b1dbe349a646f106a36b3474a453 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 26 Mar 2024 18:57:27 -0400 Subject: [PATCH 12/18] Migrate interactive path diagram example from Jupyter Co-authored-by: Sai Krishna --- AUTHORS | 1 + doc/sphinx/conf.py | 37 +++ .../reactors/interactive_path_diagram.py | 284 ++++++++++++++++++ 3 files changed, 322 insertions(+) create mode 100644 samples/python/reactors/interactive_path_diagram.py diff --git a/AUTHORS b/AUTHORS index 1f3cbf71b1..792a71a898 100644 --- a/AUTHORS +++ b/AUTHORS @@ -60,6 +60,7 @@ Ingmar Schoegl (@ischoegl), Louisiana State University Santosh Shanbhogue (@santoshshanbhogue), Massachusetts Institute of Technology Travis Sikes (@tsikes) Harsh Sinha (@sin-ha) +Sai Krishna Sirumalla (@skrsna), Northeastern University David Sondak Raymond Speth (@speth), Massachusetts Institute of Technology Su Sun (@ssun30), Northeastern University diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py index 07698ebce3..13e3081d72 100644 --- a/doc/sphinx/conf.py +++ b/doc/sphinx/conf.py @@ -14,6 +14,7 @@ import sys, os, re from pathlib import Path from sphinx_gallery.sorting import ExplicitOrder +from sphinx_gallery.scrapers import figure_rst # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the @@ -46,6 +47,41 @@ 'sphinx_copybutton', ] + +class GraphvizScraper(): + """ + Capture Graphviz objects that are assigned to variables in the global namespace. + """ + def __init__(self): + # IDs of graphviz objects that have already been seen and processed + self.processed = set() + + def __repr__(self): + return 'GraphvizScraper' + + def __call__(self, block, block_vars, gallery_conf): + import graphviz + # We use a list to collect references to image names + image_names = list() + + # The `image_path_iterator` is created by Sphinx-Gallery, it will yield + # a path to a file name that adheres to Sphinx-Gallery naming convention. + image_path_iterator = block_vars['image_path_iterator'] + + # Define a list of our already-created figure objects. + for obj in block_vars["example_globals"].values(): + if isinstance(obj, graphviz.Source) and id(obj) not in self.processed: + self.processed.add(id(obj)) + image_path = Path(next(image_path_iterator)).with_suffix(".svg") + obj.format = "svg" + obj.render(image_path.with_suffix("")) + image_names.append(image_path) + + # Use the `figure_rst` helper function to generate the reST for this + # code block's figures. + return figure_rst(image_names, gallery_conf['src_dir']) + + sphinx_gallery_conf = { 'filename_pattern': '\.py', 'example_extensions': {'.py', '.cpp', '.h', '.c', '.f', '.f90', '.m'}, @@ -54,6 +90,7 @@ 'image_srcset': ["2x"], 'remove_config_comments': True, 'ignore_repr_types': r'matplotlib\.(text|axes|legend)', + 'image_scrapers': ('matplotlib', GraphvizScraper()), 'examples_dirs': [ '../samples/python/', '../samples/cxx/', diff --git a/samples/python/reactors/interactive_path_diagram.py b/samples/python/reactors/interactive_path_diagram.py new file mode 100644 index 0000000000..915b25b7ba --- /dev/null +++ b/samples/python/reactors/interactive_path_diagram.py @@ -0,0 +1,284 @@ +""" +Interactive Reaction Path Diagrams +================================== + +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 + +.. tags:: Python, combustion, reactor network, plotting, reaction path analysis + +.. tip:: + To try the interactive features, download the Jupyter notebook version of this + example: :download:`interactive_path_diagram.ipynb`. +""" + +# %% +import numpy as np +from scipy import integrate +import graphviz +import os +from matplotlib import pyplot as plt +from collections import defaultdict +import cantera as ct + +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 import get_ipython + if "IPKernelApp" not in get_ipython().config: + raise ImportError("console") + if "VSCODE_PID" in os.environ: + raise ImportError("vscode") +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 widgets, interact + + +# %% +# 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: + +gas = ct.Solution("gri30.yaml") + +# %% +# Use Shock tube ignition delay measurement conditions corresponding to the experiments +# by Spadaccini and Colket [1]_. +# +# * CH₄-C₂H₆-O₂-Ar (3.29%-0.21%-7%-89.5%) +# * :math:`\phi` = 1.045 +# * P = 6.1 - 7.6 atm +# * T = 1356 - 1688 K + +# Set temperature, pressure, and composition +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 and Colket (1994). + +residence_time = 1e-3 + +# %% +# 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 + +# %% +# Store time, pressure, temperature and mole fractions + +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 +# --------------------------------- +# +# 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 = 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 + + diagram = ct.ReactionPathDiagram(gas, element) + diagram.threshold = threshold + diagram.title = f"time = {time:.2g} s" + + diagram.show_details = details + graph = graphviz.Source(diagram.get_dot()) + if is_interactive: + display(graph) + else: + 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, + ), + ) +else: + # For non-interactive use, just draw the diagram for a specified time step + diagram = plot_reaction_path_diagrams( + plot_step=100, + threshold=0.1, + details=False, + element="C" + ) + +# %% +# 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 = C2H6_stoichiometry.nonzero()[0] + +# %% +# The following cell calculates net rates of progress of reactions containing the +# species of interest and stores them. + +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] + ) + +# 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(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("<=>", "="), + 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.show() + +if is_interactive: + interact( + plot_instantaneous_fluxes, + annotation_cutoff=widgets.FloatSlider(value=0.1, min=1e-2, max=4, steps=10), + profiles=widgets.fixed(profiles) + ) +else: + plot_instantaneous_fluxes(annotation_cutoff=0.1, profiles=profiles) + +# %% +# 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(profiles["C2H6_production_rates"]), + np.array(profiles["time"]), + axis=0, + initial=0, +) + +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, + fontsize=8) + + plt.xlabel("Time (s)") + plt.ylabel("Integrated net rate of progress") + plt.title("Cumulative C₂H₆ formation") + 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 + ), + profiles=widgets.fixed(profiles), + integrated_fluxes=widgets.fixed(integrated_fluxes) + ) +else: + plot_integrated_fluxes( + profiles=profiles, + integrated_fluxes=integrated_fluxes, + annotation_cutoff=1e-5 + ) + +# %% +# References +# ---------- +# .. [1] L. J. Spadaccini and M. B. Colket (1994). "Ignition delay characteristics of +# methane fuels", *Progress in Energy and Combustion Science,* 20:5, 431-460. +# Prog. Energy Combust. Sci. 20, 431. +# https://doi.org/10.1016/0360-1285(94)90011-6. From e7d8d2e7533b0de102440b35cce360e98ae830f6 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 28 Mar 2024 09:34:22 -0400 Subject: [PATCH 13/18] [CI] Update packages needed for examples --- .github/workflows/main.yml | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 02d4f52e49..4d5ac69a11 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -315,7 +315,7 @@ jobs: libblas-dev liblapack-dev libhdf5-dev libfmt-dev libsundials-dev - name: Install conda dependencies run: | - mamba install -q doxygen=1.9.7 scons pip + mamba install -q doxygen=1.9.7 scons pip scikits.odes - name: Upgrade pip run: pip install -U pip setuptools wheel - name: Install Python dependencies @@ -323,7 +323,7 @@ jobs: pip install ruamel.yaml scons numpy cython 'sphinx>=7.2,<8' \ sphinxcontrib-matlabdomain sphinxcontrib-doxylink sphinxcontrib-bibtex \ pydata-sphinx-theme==0.14.1 sphinx-argparse sphinx_design myst-nb \ - sphinx-copybutton matplotlib pandas scipy pint + sphinx-copybutton matplotlib pandas scipy pint coolprop graphviz pip install "git+https://github.com/sphinx-gallery/sphinx-gallery.git@master" pip install "git+https://github.com/Cantera/sphinx-tags.git@main" - name: Build Cantera @@ -409,12 +409,15 @@ jobs: run: python3 -m pip install -U pip setuptools wheel - name: Install Python dependencies run: | - python3 -m pip install numpy ruamel.yaml pandas pyarrow matplotlib scipy pint graphviz + python3 -m pip install numpy ruamel.yaml pandas pyarrow matplotlib scipy \ + pint graphviz CoolProp python3 -m pip install --pre --no-index --find-links dist cantera - name: Run the examples - # See https://unix.stackexchange.com/a/392973 for an explanation of the -exec part + # See https://unix.stackexchange.com/a/392973 for an explanation of the -exec part. + # Skip 1D_packed_bed.py due to difficulty installing scikits.odes. run: | ln -s libcantera_shared.so build/lib/libcantera_shared.so.3 + rm samples/python/reactors/1D_packed_bed.py export LD_LIBRARY_PATH=build/lib find samples/python -type f -iname "*.py" \ -exec sh -c 'for n; do echo "$n" | tee -a results.txt && python3 "$n" >> results.txt || exit 1; done' sh {} + From 6952a5893917049d38fb1479803bad647752ebcd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 17 Apr 2024 09:55:39 -0400 Subject: [PATCH 14/18] Migrate CSTR example from Jupyter Co-authored-by: Santosh Shanbhogue --- .../samples/stirred-reactor-network.svg | 163 ++++++++++ .../images/samples/stirred-reactor-solo.svg | 124 ++++++++ samples/python/reactors/continuous_reactor.py | 295 ++++++++++++++++++ 3 files changed, 582 insertions(+) create mode 100644 doc/sphinx/_static/images/samples/stirred-reactor-network.svg create mode 100644 doc/sphinx/_static/images/samples/stirred-reactor-solo.svg create mode 100644 samples/python/reactors/continuous_reactor.py diff --git a/doc/sphinx/_static/images/samples/stirred-reactor-network.svg b/doc/sphinx/_static/images/samples/stirred-reactor-network.svg new file mode 100644 index 0000000000..e1dfd43226 --- /dev/null +++ b/doc/sphinx/_static/images/samples/stirred-reactor-network.svg @@ -0,0 +1,163 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Mixture tank + + + + + + + + + + + + + + + + + + + + + Stirred reactor + Pressure valve + Capture tank + + \ No newline at end of file diff --git a/doc/sphinx/_static/images/samples/stirred-reactor-solo.svg b/doc/sphinx/_static/images/samples/stirred-reactor-solo.svg new file mode 100644 index 0000000000..61808ff7d7 --- /dev/null +++ b/doc/sphinx/_static/images/samples/stirred-reactor-solo.svg @@ -0,0 +1,124 @@ + + + + + + + + + + + + + + + + + + + + + ṁ(t) + ṁ(t) + Inlet + Outlet + + + + + + + + + + + + P, V, T, Yᵢ + + + + + + + + + \ No newline at end of file diff --git a/samples/python/reactors/continuous_reactor.py b/samples/python/reactors/continuous_reactor.py new file mode 100644 index 0000000000..258237a552 --- /dev/null +++ b/samples/python/reactors/continuous_reactor.py @@ -0,0 +1,295 @@ +r""" +Continuously Stirred Tank Reactor +================================= + +In this example we will illustrate how Cantera can be used to simulate a +`continuously stirred tank reactor `__ +(CSTR), also interchangeably referred to as a perfectly stirred reactor (PSR), a well +stirred reactor (WSR), a jet stirred reactor (JSR), or a +`Longwell reactor `__ (and +there may well be more "aliases"). + +Requires: cantera >= 3.0, matplotlib >= 2.0, pandas + +.. tags:: Python, combustion, reactor network, well-stirred reactor + +Simulation of a CSTR/PSR/WSR +---------------------------- + +A diagram of a CSTR is shown below: + +.. image:: /_static/images/samples/stirred-reactor-solo.svg + :width: 50% + :alt: A continuously stirred tank reactor with inflow and outflow + :align: center + +As the figure illustrates, this is an open system (unlike a batch reactor, which is +isolated). *P*, *V* and *T* are the reactor's pressure, volume and temperature +respectively. The mass flow rate at which reactants come in is the same as that of the +products which exit, and on average this mass stays in the reactor for a characteristic +time :math:`\tau`, called the *residence time*. This is a key quantity in sizing the +reactor and is defined as follows: + +.. math:: + + \tau = \frac{m}{\dot{m}} + +where :math:`m` is the mass of the gas. +""" + +# %% +# Import modules and set plotting defaults +# ---------------------------------------- + +import time +from io import StringIO +import matplotlib.pyplot as plt +import pandas as pd +import cantera as ct + +print(f"Running Cantera version: {ct.__version__}") + + +# %% +# Define the gas +# -------------- +# +# In this example, we will work with n-C₇H₁₆/O₂/He mixtures, for which experimental data +# can be found in the paper by Zhang et al. [1]_ We will use the same mechanism reported +# in the paper. It consists of 1268 species and 5336 reactions. + +gas = ct.Solution("example_data/n-hexane-NUIG-2015.yaml") + +# %% +# Define initial conditions +# ------------------------- +# +# Inlet conditions for the gas and reactor parameters +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +reactor_temperature = 925 # Kelvin +reactor_pressure = 1.046138 * ct.one_atm # in atm. This equals 1.06 bars +inlet_concentrations = {"NC7H16": 0.005, "O2": 0.0275, "HE": 0.9675} +gas.TPX = reactor_temperature, reactor_pressure, inlet_concentrations + +residence_time = 2 # s +reactor_volume = 30.5 * (1e-2) ** 3 # m3 + +# %% +# Simulation parameters +# ^^^^^^^^^^^^^^^^^^^^^ + +# Simulation termination criterion +max_simulation_time = 50 # seconds + +# %% +# Reactor arrangement +# ------------------- +# +# We showed a cartoon of the reactor in the first figure in this notebook, but to +# actually simulate that, we need a few peripherals. A mass-flow controller upstream of +# the stirred reactor will allow us to flow gases in, and in-turn, a "reservoir" which +# simulates a gas tank is required to supply gases to the mass flow controller. +# Downstream of the reactor, we install a pressure regulator which allows the reactor +# pressure to stay within. Downstream of the regulator we will need another reservoir +# which acts like a "sink" or capture tank to capture all exhaust gases (even our +# simulations are environmentally friendly!). This arrangement is illustrated below: +# +# .. image:: /_static/images/samples/stirred-reactor-network.svg +# :width: 80% +# :alt: A complete reactor network representing a continuously stirred tank reactor +# :align: center +# +# Initialize the stirred reactor and connect all peripherals +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +fuel_air_mixture_tank = ct.Reservoir(gas) +exhaust = ct.Reservoir(gas) + +stirred_reactor = ct.IdealGasMoleReactor(gas, energy="off", volume=reactor_volume) + +mass_flow_controller = ct.MassFlowController( + upstream=fuel_air_mixture_tank, + downstream=stirred_reactor, + mdot=stirred_reactor.mass / residence_time, +) + +pressure_regulator = ct.PressureController( + upstream=stirred_reactor, + downstream=exhaust, + primary=mass_flow_controller, + K=1e-3, +) + +reactor_network = ct.ReactorNet([stirred_reactor]) + +# Create a SolutionArray to store the data +time_history = ct.SolutionArray(gas, extra=["t"]) + +# Set the maximum simulation time +max_simulation_time = 50 # seconds + +# Start the stopwatch +tic = time.time() + +# Set simulation start time to zero +t = 0 +counter = 1 +while t < max_simulation_time: + t = reactor_network.step() + + # We will store only every 10th value. Remember, we have 1200+ species, so there + # will be 1200+ columns for us to work with + if counter % 10 == 0: + # Extract the state of the reactor + time_history.append(stirred_reactor.thermo.state, t=t) + + counter += 1 + +# Stop the stopwatch +toc = time.time() + +print(f"Simulation Took {toc-tic:3.2f}s to compute, with {counter} steps") + +# %% +# Plot the results +# ^^^^^^^^^^^^^^^^ +# +# As a test, we plot the mole fraction of CO and see if the simulation has converged. If +# not, go back and adjust max. number of steps and/or simulation time. + +plt.figure() +plt.semilogx(time_history.t, time_history("CO").X, "-o") +plt.xlabel("Time (s)") +plt.ylabel("Mole Fraction : $X_{CO}$") + +# %% +# Illustration : Modeling experimental data +# ----------------------------------------- +# +# Let us see if the reactor can reproduce actual experimental measurements. +# +# We first load the data. This is also supplied in the paper by Zhang et al. [1]_ as an +# excel sheet + +experimental_data_csv = """ +T,NC7H16,O2,CO,CO2 +500,5.07E-03,2.93E-02,0.00E+00,0.00E+00 +525,4.92E-03,2.86E-02,0.00E+00,0.00E+00 +550,4.66E-03,2.85E-02,0.00E+00,0.00E+00 +575,4.16E-03,2.63E-02,2.43E-04,1.01E-04 +600,3.55E-03,2.33E-02,9.68E-04,2.51E-04 +625,3.36E-03,2.31E-02,1.42E-03,2.67E-04 +650,3.67E-03,2.45E-02,9.16E-04,1.46E-04 +675,4.38E-03,2.77E-02,2.25E-04,0.00E+00 +700,4.79E-03,2.87E-02,0.00E+00,0.00E+00 +725,4.89E-03,2.93E-02,0.00E+00,0.00E+00 +750,4.91E-03,2.84E-02,0.00E+00,0.00E+00 +775,4.93E-03,2.80E-02,0.00E+00,0.00E+00 +800,4.78E-03,2.82E-02,0.00E+00,0.00E+00 +825,4.41E-03,2.80E-02,1.49E-05,0.00E+00 +850,3.68E-03,2.80E-02,4.18E-04,1.66E-04 +875,2.13E-03,2.45E-02,1.65E-03,2.22E-04 +900,1.03E-03,2.05E-02,5.51E-03,3.69E-04 +925,5.82E-04,1.79E-02,8.59E-03,6.78E-04 +950,3.88E-04,1.47E-02,1.05E-02,1.07E-03 +975,2.35E-04,1.28E-02,1.19E-02,1.36E-03 +1000,1.14E-04,1.16E-02,1.34E-02,1.82E-03 +1025,4.83E-05,9.88E-03,1.52E-02,2.41E-03 +1050,1.64E-05,8.16E-03,1.83E-02,2.97E-03 +1075,1.22E-06,5.48E-03,1.95E-02,3.67E-03 +1100,0.00E+00,3.24E-03,2.14E-02,4.38E-03 +""" + +experimental_data = pd.read_csv(StringIO(experimental_data_csv)) +experimental_data.head() + +# %% + +# Define all the temperatures at which we will run simulations. These should overlap +# with the values reported in the paper as much as possible +T = [650, 700, 750, 775, 825, 850, 875, 925, 950, 1075, 1100] + +# Create a SolutionArray to store values for the above points +temp_dependence = ct.SolutionArray(gas) + +# %% +# Now we simply run the reactor code we used above for each temperature + +concentrations = inlet_concentrations + +for reactor_temperature in T: + # Use concentrations from the previous iteration to speed up convergence + gas.TPX = reactor_temperature, reactor_pressure, concentrations + + stirred_reactor = ct.IdealGasReactor(gas, energy="off", volume=reactor_volume) + fuel_air_mixture_tank = ct.Reservoir(gas) + mass_flow_controller = ct.MassFlowController( + upstream=fuel_air_mixture_tank, + downstream=stirred_reactor, + mdot=stirred_reactor.mass / residence_time, + ) + pressure_regulator = ct.PressureController( + upstream=stirred_reactor, downstream=exhaust, primary=mass_flow_controller + ) + reactor_network = ct.ReactorNet([stirred_reactor]) + + # Re-run the isothermal simulations + tic = time.time() + reactor_network.advance(max_simulation_time) + toc = time.time() + print(f"Simulation at T={reactor_temperature}K took {toc-tic:3.2f}s to compute") + + concentrations = stirred_reactor.thermo.X + temp_dependence.append(stirred_reactor.thermo.state) + +# %% +# Compare the model results with experimental data +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +plt.figure() +plt.plot( + temp_dependence.T, temp_dependence("NC7H16").X, color="C0", label="$nC_{7}H_{16}$" +) +plt.plot(temp_dependence.T, temp_dependence("CO").X, color="C1", label="CO") +plt.plot(temp_dependence.T, temp_dependence("O2").X, color="C2", label="O$_{2}$") + +plt.plot( + experimental_data["T"], + experimental_data["NC7H16"], + color="C0", + marker="o", + label="$nC_{7}H_{16}$ (exp)", +) +plt.plot( + experimental_data["T"], + experimental_data["CO"], + color="C1", + marker="^", + linestyle="none", + label="CO (exp)", +) +plt.plot( + experimental_data["T"], + experimental_data["O2"], + color="C2", + marker="s", + linestyle="none", + label="O$_2$ (exp)", +) + +plt.xlabel("Temperature (K)") +plt.ylabel(r"Mole Fractions") + +plt.xlim([650, 1100]) +plt.legend(loc=1) + +# %% +# References +# ---------- +# +# .. [1] K. Zhang, C. Banyon, C. Togbé, P. Dagaut, J. Bugler, H. J. Curran (2015). "An +# experimental and kinetic modeling study of n-hexane oxidation," *Combustion and +# Flame* 162:11, 4194-4207, https://doi.org/10.1016/j.combustflame.2015.08.001. + +# sphinx_gallery_thumbnail_number = -1 From a2dcfabe52401a24660907d6c2db08fb22fda319 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 18 Apr 2024 12:54:13 -0400 Subject: [PATCH 15/18] [CI] Fix failing Codecov uploads --- .github/workflows/main.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 4d5ac69a11..859f724474 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -280,9 +280,10 @@ jobs: path: interfaces/dotnet/coveragereport* retention-days: 5 - name: Upload Coverage to Codecov - uses: codecov/codecov-action@v3 + uses: codecov/codecov-action@v4 with: verbose: true + token: ${{secrets.CODECOV_TOKEN}} files: ./coverage.xml,./build/pycov.xml,./interfaces/dotnet/coverage.cobertura.xml fail_ci_if_error: true From 7bed9d4170110cede42af1223bd76d13c050b3a8 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 18 Apr 2024 14:29:57 -0400 Subject: [PATCH 16/18] [CI] Avoid errors from running certain examples --- .github/workflows/main.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 859f724474..97bc3a6298 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -416,9 +416,13 @@ jobs: - name: Run the examples # See https://unix.stackexchange.com/a/392973 for an explanation of the -exec part. # Skip 1D_packed_bed.py due to difficulty installing scikits.odes. + # Skip continuous_reactor.py due to thermo warnings + # Increase figure limit to handle flame_speed_convergence_analysis.py. run: | ln -s libcantera_shared.so build/lib/libcantera_shared.so.3 rm samples/python/reactors/1D_packed_bed.py + rm samples/python/reactors/continuous_reactor.py + echo "figure.max_open_warning: 100" > matplotlibrc export LD_LIBRARY_PATH=build/lib find samples/python -type f -iname "*.py" \ -exec sh -c 'for n; do echo "$n" | tee -a results.txt && python3 "$n" >> results.txt || exit 1; done' sh {} + From 993536a0569bb844a72077d43df2503e98e671dc Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 18 Apr 2024 15:43:37 -0400 Subject: [PATCH 17/18] [CI] Show Sphinx warnings separately --- .github/workflows/main.yml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 97bc3a6298..509eec4151 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -330,7 +330,13 @@ jobs: - name: Build Cantera run: scons build -j2 debug=n optimize=y use_pch=n - name: Build documentation - run: scons sphinx doxygen logging=debug + run: | + scons sphinx doxygen logging=debug \ + sphinx_options="-W --keep-going --warning-file=sphinx-warnings.txt" + - name: Show Sphinx warnings + run: | + cat sphinx-warnings.txt + if: failure() - name: Ensure 'scons help' options work run: | scons help --options From 923af6f74ba3169d3a14570b32136036a7488697 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 18 Apr 2024 19:46:02 -0400 Subject: [PATCH 18/18] Ignore warning about caching Sphinx config --- doc/sphinx/conf.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py index 13e3081d72..d03da1084d 100644 --- a/doc/sphinx/conf.py +++ b/doc/sphinx/conf.py @@ -119,6 +119,8 @@ def __call__(self, block, block_vars, gallery_conf): } } +suppress_warnings = ["config.cache"] # Triggered by objects in sphinx_gallery_conf + # Override sphinx-gallery's method for determining which examples should be executed. # There's really no way to achieve this with the `filename_pattern` option, and # `ignore_pattern` excludes the example entirely.