1D Flame sensitivity analysis for mass fractions

[Antoine]

Hello,

I would like to know if there is a way in Cantera to perform a sensitivity analysis similar to Chemkin to obtain sensitivity coefficients for the concentration of a given species?

I was previously using a brute-force approach for small mechanisms, but I am now looking at larger detailed mechanisms, and the time to perform a single brute-force analysis is prohibitive. I saw the get_flame_speed_reaction_sensitivities and was wondering if an analogous method exists for a different target than flame speed

Thanks,

Antoine

[Nick Curtis]

Hi Antoine,

Is this python sensitivity example is what you're looking for?

Best,

Nick

[Antoine]

Hi Nick,

I am looking to obtain sensitivity indeces in 1D flames. From my understanding sensitivity and sensitivities are only applicable to 0D reactors.

But yes, this is what I would be looking for in freely-propagating flames.

Thanks,

Antoine

[Nick Curtis]

Hi Antoine,

It appears that 

get_flame_speed_reaction_sensitivities

is merely a thin-wrapper around the solve_adjoint function, which as far as I can tell _should_ be defined for all 1-D flames.

I think you could probably write your own version fairly easily, as long as your target was in the solution vector.

Nick

[Antoine]

Hi Nick,

Thanks for the answer. This is what I figured out after experimenting a bit with solve_adjoint.

Antoine

[Rodolfo Rocha]

Hi everyone,

I have already asked it in another topic, but maybe you could help me with that.

I created a function based on get_flame_speed_reaction_sensitivities, replacing the 'u' references with a species variable in a separate function (for species a string like 'CO', 'CO2', etc.). My code is this:

def get_species_reaction_sensitivities(self, species, grid_point):

     r"""

   Compute the normalized sensitivities of the species production

 :math:`s_{i, spec}` with respect to the reaction rate constants :math:`k_i`:

   .. math::

              s_{i, spec} = \frac{k_i}{[X]} \frac{d[X]}{dk_i}

"""

        def g(sim):

            return sim.X[self.gas.species_index(species), grid_point]

        Nvars = sum(D.n_components * D.n_points for D in self.domains)

        # Index of species in the global solution vector

  i_spec = self.inlet.n_components + self.flame.component_index(species)

        dgdx = np.zeros(Nvars)

 dgdx[i_spec] = 1

        spec_0 = g(self)

        def perturb(sim, i, dp):

               sim.gas.set_multiplier(1+dp, i)

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

For 'self' being the flame object. It gives me the same value divided by the species concentration in each point, which means that the value of

self.solve_adjoint(perturb, self.gas.n_reactions, dgdx)

is always the same.

For you who figured it out, how did you solve it? What else should I change in the function?

Thank you very much,

[Antoine]

Hello Rodolfo,

From what I can see here, you are not changing the solution node position in the index i_spec. The solution is stored in the following way

    #     point 1                       point 2                      point 3        ...

    #  ______\_______     ______\_______     ______\_______

    # |All components\     |All components\      |All components\

    # |u,V,T,lambda,Xi\

So your index should corresponds to the species at the solution point you are looking at and look something like this

i_spec = f.domains[1].n_components*pointNumber + i_component

where you change the pointNumber to match your node and the i_component to the index of the component of interest.

Hope this helps,

Antoine

[Rodolfo Rocha]

Hello Antoine,

Thank you very much for your response.

I have tried this approach, but it gave me very strange values. For you to see, this sensitivity values I got were way higher than 1. In my case, I'm looking for mole fraction sensitivity. My function is the following:

def get_species_reaction_sensitivities(self, species, grid_point):

r"""

Compute the normalized sensitivities of the species production

:math:`s_{i, spec}` with respect to the reaction rate constants :math:`k_i`:

.. math::

s_{i, spec} = \frac{k_i}{[X]} \frac{d[X]}{dk_i}

"""

def g(sim):

return sim.X[self.gas.species_index(species), grid_point]

Nvars = sum(D.n_components * D.n_points for D in self.domains)

pos = int(Nvars / len(self.grid))

# Index of spec in the global solution vector

# i_spec = self.inlet.n_components + self.flame.component_index(species) + grid_point * pos

i_spec = self.inlet.n_components + self.flame.component_index(species) + self.domains[1].n_components*grid_point

dgdx = np.zeros(Nvars)

dgdx[i_spec] = 1

spec_0 = g(self)

def perturb(sim, i, dp):

sim.gas.set_multiplier(1+dp, i)

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

Which for the CO molar fraction sensitivity of the first 10 reactions of GRI-mech 3.0 results in the graphic I've attached.

What do you think is wrong now? Do you think that maybe dividing by the molar fraction is what is giving me those values? How can I fix it?

Thanks again for the help,

[Antoine]

Hello Rodolfo,

What type of burner is this? I suspect that the flame is not lit at 0mm from the burner. In this context, you have concentrations close to zero before the flame so it might be why it is happening.

In your case, if you are looking to do sensitivity analysis of the entire profile, I would suggest that you look into brute force sensitivity analysis. Especially with GRI3.0 as it is pretty fast. You can validate your results using the adjoint method.

Antoine

[Rodolfo]

Hello Antoine,

I'm using, for this case, the FreeFlame simulation for 1D premixed flames. Actually, I want to make a bar plot of the sensitivity in the flame region, the full profile was taken in order to see if things had worked well, such as in the example shown below, from another topic. I'm also going to do it for BurnerFlame, for a longer tube with an experimental temperature profile.

You were right, the flame did not ignite at 0 mm from the burner. As you can see in the temperature profile, it ignites after 20 mm. Do I have to do a brute force sensitivity for this case? How can I do it?

--
You received this message because you are subscribed to a topic in the Google Groups "Cantera Users' Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/cantera-users/vfejZRDKKuw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cantera-user...@googlegroups.com.
To post to this group, send email to canter...@googlegroups.com.
Visit this group at https://groups.google.com/group/cantera-users.
For more options, visit https://groups.google.com/d/optout.

--

Rodolfo Cavaliere da Rocha

[Antoine]

Hello Rodolfo,

Here is a link to a detailed example of brute-force sensitivity analysis  with Cantera (https://charlesreid1.com/wiki/Cantera/Sensitivity_Analysis). I suggest that you try with the two methods to ensure that your results are consistent. 

If you have an example from a different setup, maybe you should try to reproduce the SA with the adjoint method before looking at your problem. 

Antoine

[Rodolfo Rocha]

Hello Antoine,

I have managed to solve it. Sensitivity values in my code are good in the flame region, but not before the flame. From the link you sent me I have noticed that they use a different value for the perturbation. In my case, dgdx[i_spec] = 1, and in theirs, dgdx[i_spec] = 0.05. It does not change the shape of the plots, but gives a fraction of the values.

For a flame of BurnerFlame, the only change I needed was to put self.burner.n_components instead of self.inlet.n_components. It worked like a charm.

Thank you very much for your help!

Best regards,