From 32a6d3da4efa066fee1931ee9714a395eb79cc91 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 18 Mar 2024 23:49:46 -0400 Subject: [PATCH] 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 000000000..38b5ef653 --- /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 bdbaea11a..0d9433677 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 >= 2.5.0 +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`. + +.. 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' +# --------------------- -width = 0.03 # m +# 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} -# Solution object used to compute mixture properties -gas = ct.Solution('gri30.yaml') -gas.TPX = Tin, p, reactants +# 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 -# Flame object -f = ct.FreeFlame(gas, width=width) -f.set_refine_criteria(ratio=3, slope=0.07, curve=0.14) +# Domain width in metres +width = 0.014 -f.solve(loglevel=1, auto=True) -print(f"\nmixture-averaged flamespeed = {f.velocity[0]:7f} m/s\n") +# %% +# 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") + + +# %% +# Plot temperature and major species profiles +# ------------------------------------------- +# +# Check and see if all has gone well. + +# Extract the spatial profiles as a SolutionArray to simplify plotting specific species +profile = flame.to_array() + +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() -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}") +# 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') + +#sphinx_gallery_thumbnail_number = -1