Net Reaction Sensitivities per Species for 1D

[marina8888]

Abstract

Tutorial or example on extracting (net) species sensitivities in 1D Flames with respect to certain reactions.

Motivation

I've been struggling to find a way to extract the net sensitivity of a specific species to certain reactions in 1D flames. This is used for analysing and improving mechanisms. For example, if experimental data predicts a different production of a certain species, it is possible to find the main reactions that impact the production of that species, and hence focus on improving those specific reactions.

Possible Solutions

Is there a way this is possible with the solve_adjoint method? I can see that this method was already implemented with the FlameBase.get_flame_speed_reaction_sensitivities method, but for laminar burning velocity instead of for a specific species. I would try to contribute the code myself, but I am still very new to using Cantera.

References

The question has also been asked by a different user below, and is especially useful for those transitioning from CHEMKIN:
https://groups.google.com/g/cantera-users/c/wd171ne9DVk

Thank you in advance for your help.

[bryanwweber]

This would be really awesome! One approach is to take the existing code for flame speed sensitivities:

https://github.com/Cantera/cantera/blob/7142a7b0c4d5ae06b8129a835f8b84bf3b9842aa/interfaces/cython/cantera/onedim.py#L528-L554

and copy-paste it into your file where you're running the flame. Then, change all the references to self to refer to the flame object, and change all the references to u to be to the species of interest. I haven't actually tried that, but that should at least get you some error messages you can start debugging.

The easiest way to work on this will most likely be in a Jupyter Notebook. If you're familiar with the GitHub workflow and want to make a pull request to the https://github.com/cantera/cantera-jupyter repository, I'd be happy to help review and debug the code!

[marina8888]

Ok, I will have a go, and see how I get on. Thank you for the links and pointing me in the right direction. I haven't actually used Jupyter Notebook before, so if it works with you, I will start with some python code to inherit from the FreeFlame class and debug from there.

[bryanwweber]

Sure! I don't think you need to inherit from anything, to be clear. I think starting from the flame speed example (https://cantera.org/examples/python/onedim/adiabatic_flame.py.html) and literally copy-pasting that function into that example, and calling it with the appropriate arguments replaced, is going to be the easiest way to start.

[marina8888]

HI @bryanwweber, just to check, are there any instructions on how to build Cantera from soure?

[bryanwweber]

Yes, on our website at https://cantera.org/install or more specifically https://cantera.org/compiling/installation-reqs.html#sec-installation-reqs

[marina8888]

Hello @bryanwweber. I have replaced references in the get_flame_speed_reaction_sensitivities function from flame speed u, to a species of a user's choice (at distance Z along the domain). However, after comparing between known species sensitivity results, the function seems to give the incorrect values. Would you have any suggestions on what might be going wrong under the hood deeper in cantera?

def get_flame_speed_reaction_sensitivities(self, species_name, distance_Z):
        r"""
        Compute the normalized sensitivities of the laminar flame speed
        :math:`S_u` with respect to the reaction rate constants :math:`k_i`:

        .. math::

            s_i = \frac{k_i}{S_u} \frac{dS_u}{dk_i}
        """

        def g(sim, species_name:str, distance):
            # finds grid point closest to distance and set at that gas value:
            index = min(range(len(self.flame.grid)), key=lambda i: abs(self.flame.grid[i]-distance))
            self.set_gas_state(index)

            # sets type of gas and converts it from list
            index = self.gas.species_names.index(species_name)
            val = list(getattr(sim.gas[index], 'X'))
            val = val[0]

            return val

        # sum across grid:
        Nvars = sum(D.n_components * D.n_points for D in self.domains)

        # Index of gas[0] in the global solution vector
        i_Su = self.inlet.n_components + self.flame.component_index(species_name)

        dgdx = np.zeros(Nvars)
        dgdx[i_Su] = 1

        Su0 = g(self, species_name = species_name, distance = distance_Z)

        def perturb(sim, i, dp):
            sim.gas.set_multiplier(1+dp, i)

        return self.solve_adjoint(perturb, self.gas.n_reactions, dgdx) / Su0

[bryanwweber]

@marina8888 Are you trying to solve for the sensitivity of the flame speed to perturbing species concentrations, or the sensitivity of a species concentration to perturbing reaction rates? Right now, you're solving for the latter, but I just want to make sure since I wasn't quite sure from your initial post.

Assuming that is the case, I created this modified version of the adiabatic_flame.py example, based on your code above:

import cantera as ct
import numpy as np

# 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

# Solution object used to compute mixture properties
gas = ct.Solution("gri30.yaml")
gas.TPX = Tin, p, reactants

# Flame object
f = ct.FreeFlame(gas, width=width)
f.set_refine_criteria(ratio=3, slope=0.07, curve=0.14)

f.solve(loglevel=1, auto=True)
print("\nmixture-averaged flamespeed = {:7f} m/s\n".format(f.velocity[0]))


# Use the adjoint method to calculate sensitivities
# sens = f.get_flame_speed_reaction_sensitivities()
def get_species_reaction_sensitivities(sim, species_name, distance):
    r"""
    Compute the normalized sensitivities of the laminar flame speed
    :math:`S_u` with respect to the reaction rate constants :math:`k_i`:

    .. math::

        s_i = \frac{k_i}{S_u} \frac{dS_u}{dk_i}
    """

    def g(sim, species_name: str, distance):
        # finds grid point closest to distance and set at that gas value:
        index = np.argmin(np.abs(sim.flame.grid - distance))
        sim.set_gas_state(index)

        # sets type of gas and converts it from list
        index = sim.gas.species_index(species_name)
        val = sim.gas.X[index]
        return val

    # sum across grid:
    Nvars = sum(D.n_components * D.n_points for D in sim.domains)

    # Index of gas[0] in the global solution vector
    i_Su = sim.inlet.n_components + sim.flame.component_index(species_name)

    dgdx = np.zeros(Nvars)
    dgdx[i_Su] = 1

    Su0 = g(sim, species_name=species_name, distance=distance)

    def perturb(sim, i, dp):
        sim.gas.set_multiplier(1 + dp, i)

    return sim.solve_adjoint(perturb, sim.gas.n_reactions, dgdx) / Su0


sens = get_species_reaction_sensitivities(f, "CH4", 0.015)


print()
print("Rxn #   k/S*dS/dk    Reaction Equation")
print("-----   ----------   ----------------------------------")
for m in range(gas.n_reactions):
    print("{: 5d}   {: 10.3e}   {}".format(m, sens[m], gas.reaction_equation(m)))

I'm not sure what values are appropriate, so I'm not sure what to look for here. How are you generating the "known" sensitivities?

[marina8888]

Thank you very much for sending the improved implementation of the code. I am seeing some strange behaviour which I still need to look into before making a PR. Please could it be possible to check that you are seeing the same behaviour when you run the script?

• when using FreeFlame as the model:
  Values printed are actually the sensitivities for flame speed and not for the input species. For example, when using 'HCO' as the species at distance = end of the reactor, the top reactions should be as listed below (graph found in literature):
  

But instead, the results produced are almost identical to the flame speed sensitivity, as shown below (graph found in literature): :

- when using the ImpingingJet model:
This calculates the correct results with respect to the input species (validated for NO and NH3), but only when the mdot input value is 1/4 of what I believe to be the true value of the mdot. Running at the true mdot value results in the sensitivities showing as mostly 0, or just incorrect. Equally running with mdot too low also gives incorrect sensitivities.

Impinging Jet flame experiments are usually conducted at around the laminar burning velocity as the input velocity to prevent flashback. For example, methane has 35 cm/s laminar burning velocity at stoichiometry (and multiplied by the approximate density of 1kg/m^3 for the premixed air/fuel solution) should give mdot value of 0.35kg/m^2/s. However, the example script below uses below 1/4 of this value:
stagnation flame example script
Would you happen to know why this might be the case, or where I may be miscalculating?

[bryanwweber]

Hi! I'm getting values around 1E-20 to 1E-15 using the script I posted in my last comment. That's clearly not the same order of magnitude as the graphs you're showing from the literature. Unfortunately, this is outside my area of expertise (adjoint sensitivity specifically) so I'm not sure where to go from here. If you understand how adjoint sensitivity should be implemented, maybe you could try to create your own version in pure Python? It might be slow, but a slow + correct answer is better than a fast + wrong answer.

[marina8888]

Hi @jiweiqi apologies for the late reply - I have honestly lost track of this since realising that this may not be as easy of an issue to solve as it may appear on the surface.

In the case of my research work, I am using the approach of finding reactions with low sensitivity to laminar burning velocity (for a mechanism tuned to laminar burning velocity) and high sensitivity to species where there are inaccuracies compared to experimental data.

Thank you for sending across the paper - I will read it and consider its implications for my work as well.

I will close this issue, given that this might not be something I can solve directly by tweaking the source code. Thank you @bryanwweber for all your help thus far.