Species Sensitivity Analysis for 1D Flames

AUTHORS.md

@@ -44,6 +44,7 @@ update, please report on Cantera's
 - **Cory Kinney** [@corykinney](https://github.com/corykinney)
 - **Gandhali Kogekar** [@gkogekar](https://github.com/gkogekar)
 - **Daniel Korff** [@korffdm](https://github.com/korffdm) - Colorado School of Mines
+- **Marina Kovaleva** [@marina8888](https://github.com/marina8888) - Tohoku University
 - **Ashwin Kumar** [@mgashwinkumar](https://github.com/mgashwinkumar) - Virginia Tech
 - **Jon Kristofer**
 - **Samesh Lakothia** [@sameshl](https://github.com/sameshl)

interfaces/cython/cantera/onedim.py

@@ -172,6 +172,34 @@ def set_profile(self, component, positions, values):
         """
         super().set_profile(self.flame, component, positions, values)
 
+    def get_species_reaction_sensitivities(self, species: str, ix: int) -> np.ndarray:
+        r"""
+        Compute the normalized sensitivities of a species,
+        :math:`X_i` with respect to the reaction rate constants :math:`k_i`:
+        .. math::
+            s_i = \frac{k_i}{X_i} \frac{dX_i}{dk_i}
+        :param species:
+            Species of interest for sensitivity calculation (i.e "NO")
+        :param distance:
+            Index of flame grid point for which the sensitivity analysis is undertaken
+        """
+
+        def g(sim):
+            return sim.X[sim.gas.species_index(species), ix]
+
+        n_vars = sum(dd.n_components * dd.n_points for dd in self.domains)
+
+        i_X = self.flame.component_index(species) + self.domains[1].n_components * ix
+
+        dgdx = np.zeros(n_vars)
+        dgdx[i_X] = 1
+        X0 = g(self)
+
+        def perturb(sim, i, dp):
+            sim.gas.set_multiplier(1 + dp, i)
+
+        return self.solve_adjoint(perturb, self.gas.n_reactions, dgdx) / X0
+
     @property
     def max_grid_points(self):
         """
@@ -773,12 +801,12 @@ def get_flame_speed_reaction_sensitivities(self):
         def g(sim):
             return sim.velocity[0]
 
-        Nvars = sum(D.n_components * D.n_points for D in self.domains)
+        n_vars = sum(dd.n_components * dd.n_points for dd in self.domains)
 
         # Index of u[0] in the global solution vector
         i_Su = self.inlet.n_components + self.flame.component_index('velocity')
 
-        dgdx = np.zeros(Nvars)
+        dgdx = np.zeros(n_vars)
         dgdx[i_Su] = 1
 
         Su0 = g(self)
@@ -1539,4 +1567,4 @@ def set_initial_guess(self, data=None, group=None):
 
         self.set_profile('velocity', [0.0, 1.0], [uu, 0])
         self.set_profile('spread_rate', [0.0, 1.0], [0.0, a])
-        self.set_profile("lambda", [0.0, 1.0], [L, L])
+        self.set_profile("lambda", [0.0, 1.0], [L, L])

samples/python/onedim/species_sensitivity.py

@@ -0,0 +1,90 @@
+r"""
+Sensitivity analysis
+========================================
+In this example we simulate an impinging jet flame, premixed ammonia/hydrogen-air flame,
+calculating the sensitivity of N2O to the A-factor reaction rate constants.
+Requires: cantera >= 3.0.0, pandas
+.. tags:: Python, combustion, 1D flow, species reaction, premixed flame,
+          sensitivity analysis, plotting
+"""
+import cantera as ct
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+# Define the gas mixture
+gas = ct.Solution("example_data/nakamura.yaml")  # Use the Nakamura mechanism for ammonia combustion blends
+gas.TP = 295, ct.one_atm
+air = "O2:0.21,N2:0.79"
+gas.set_equivalence_ratio(phi=0.6, fuel="NH3:0.7, H2:0.3", oxidizer=air)
+flame = ct.ImpingingJet(gas=gas, width=0.02)
+flame.inlet.mdot = 0.255 * gas.density
+flame.surface.T = 493.5
+flame.set_initial_guess("equil")
+
+
+# Refine grid to improve accuracy
+flame.set_refine_criteria(ratio=3, slope=0.025, curve=0.05)
+
+# Solve the flame
+flame.solve(loglevel=1, auto=True)  # Increase loglevel for more output
+
+# Plot temperature profile
+plt.figure(figsize=(8, 6))
+plt.plot(flame.grid * 1e3, flame.T, label="Flame Temperature", color="red")
+plt.xlabel("Distance (mm)")
+plt.ylabel("Temperature (K)")
+plt.title("Temperature Profile of a Flame")
+plt.grid(True, linestyle='--', alpha=0.7)
+plt.legend()
+plt.tight_layout()
+plt.show()
+
+# Create a DataFrame to store sensitivity-analysis data
+sens = pd.DataFrame(index=gas.reaction_equations(), columns=["sensitivity"])
+
+# Use the adjoint method to calculate species sensitivities at a set distance in the flame domain
+distance = 0.02
+species = "N2O"
+sens.sensitivity = flame.get_species_reaction_sensitivities(species, distance)
+
+sens = sens.iloc[(-sens['sensitivity'].abs()).argsort()]
+fig, ax = plt.subplots()
+
+sens.head(15).plot.barh(ax=ax, legend=None)
+ax.invert_yaxis()  # put the largest sensitivity on top
+ax.set_title(f"Sensitivities for {species} Using the Nakamura mechanism")
+ax.set_xlabel(r"Sensitivity: $\frac{\partial\:\ln X_{N2O}}{\partial\:\ln k}$")
+ax.grid(axis='x')
+plt.tight_layout()
+plt.show()
+
+
+# Forward sensitivities
+dk = 3e-4
+
+# get index in the grid at distance
+ix = np.argmin(np.abs(flame.grid - distance))
+
+Su0 = flame.X[gas.species_index(species), ix]
+fwd = []
+for m in range(flame.gas.n_reactions):
+    flame.gas.set_multiplier(1.0)  # reset all multipliers
+    flame.gas.set_multiplier(1 + dk, m)  # perturb reaction m
+    flame.solve(loglevel=0, refine_grid=False)
+    Suplus = flame.X[gas.species_index(species), ix]
+    flame.gas.set_multiplier(1 - dk, m)  # perturb reaction m
+    flame.solve(loglevel=0, refine_grid=False)
+    Suminus = flame.X[gas.species_index(species), ix]
+    fwd.append((Suplus - Suminus) / (2 * Su0 * dk))
+sens = pd.DataFrame(index=gas.reaction_equations(), columns=["sensitivity with forward"], data=fwd)
+sens = sens.iloc[(-sens['sensitivity with forward'].abs()).argsort()]
+fig, ax = plt.subplots()
+
+sens.head(15).plot.barh(ax=ax, legend=None)
+ax.invert_yaxis()  # put the largest sensitivity on top
+ax.set_title(f"Sensitivities for {species} Using Nakamura Mech")
+ax.set_xlabel(r"Sensitivity: $\frac{\partial\:\ln X_{i}}{\partial\:\ln k}$")
+ax.grid(axis='x')
+plt.tight_layout()
+plt.show()

test/python/test_onedim.py

@@ -530,6 +530,53 @@ def test_adjoint_sensitivities(self):
             fwd = (Suplus-Suminus)/(2*Su0*dk)
             assert fwd == approx(dSdk_adj[m], rel=5e-3, abs=1e-7)
 
+    def test_adjoint_sensitivities2(self):
+        self.gas = ct.Solution('gri30.yaml')
+        self.gas.TP = 300, 101325
+        self.gas.set_equivalence_ratio(1.0, 'CH4', {"O2": 1.0, "N2": 3.76})
+
+        self.flame = ct.FreeFlame(self.gas, width=0.014)
+        self.flame.set_refine_criteria(ratio=3, slope=0.07, curve=0.14)
+        self.flame.solve(loglevel=0, refine_grid=False)
+
+        # Adjoint sensitivities
+        ix = -1
+        species = "NO"
+        dk = 1e-4
+
+        adjoint_sensitivities = self.flame.get_species_reaction_sensitivities(species, ix)
+        reaction_equations = self.gas.reaction_equations()
+        adjoint_sens = [(reaction, sensitivity) for reaction, sensitivity in
+                        zip(reaction_equations, adjoint_sensitivities)]
+        adjoint_sens_sorted = sorted(adjoint_sens, key=lambda x: abs(x[1]), reverse=True)
+
+        Su0 = self.flame.X[self.gas.species_index(species), ix]
+        fwd_sensitivities = []
+        for m in range(self.flame.gas.n_reactions):
+            self.flame.gas.set_multiplier(1.0)  # reset all multipliers
+            self.flame.gas.set_multiplier(1 + dk, m)  # perturb reaction m
+            self.flame.solve(loglevel=0, refine_grid=False)
+            Suplus = self.flame.X[self.gas.species_index(species), ix]
+            self.flame.gas.set_multiplier(1 - dk, m)  # perturb reaction m
+            self.flame.solve(loglevel=0, refine_grid=False)
+            Suminus = self.flame.X[self.gas.species_index(species), ix]
+            fwd_sensitivities.append((Suplus - Suminus) / (2 * Su0 * dk))
+        fwd_sens = [(reaction, sensitivity) for reaction, sensitivity in zip(reaction_equations, fwd_sensitivities)]
+        fwd_sens_sorted = sorted(fwd_sens, key=lambda x: abs(x[1]), reverse=True)
+
+        # Extract top 10 reactions from both lists
+        top_reactions_adjoint = [reaction for reaction, _ in adjoint_sens_sorted[:10]]
+        top_reactions_fwd = [reaction for reaction, _ in fwd_sens_sorted[:10]]
+
+        # Count common reactions
+        common_reactions = list(set(top_reactions_fwd) & set(top_reactions_adjoint))
+
+        # Assert that at least 8 reactions match
+        assert len(
+            common_reactions) >= 8, f"Only {len(common_reactions)} Fwd and adjoint based species sensitivities do not match"
+
+
+
     def test_jacobian_options(self):
         reactants = {'H2': 0.65, 'O2': 0.5, 'AR': 2}
         self.create_sim(p=ct.one_atm, Tin=300, reactants=reactants, width=0.03)