mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
158 lines
5.0 KiB
Python
158 lines
5.0 KiB
Python
r"""
|
|
Laminar flame speed sensitivity analysis
|
|
========================================
|
|
|
|
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`.
|
|
|
|
.. image:: /_static/images/samples/flame-speed.svg
|
|
:width: 50%
|
|
:alt: Freely Propagating Flame
|
|
:align: center
|
|
|
|
"""
|
|
|
|
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
|
|
# ---------------------
|
|
|
|
# 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")
|
|
|
|
|
|
# %%
|
|
# 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.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')
|
|
|
|
#sphinx_gallery_thumbnail_number = -1
|