sensitivity analysis of ignition delay

[Shrabanti Roy]

Hi

I am using cantera 2.3.0 with python interface.

Can anyone please suggest me a good example to start with for calculating sensitivity analysis of ignition delay time.

I checked examples of mechanism_reduction.py, sensitivity.py and reaction_path.py. I am confused which one of these examples is related with ignition delay sensitivity. Or do I need to modify or use some kind of combination of these example?

Thanks.

[Nick Curtis]

Hi Shrabanti,

None of those examples directly calculates the ignition delay sensitivity, however Bryan Webber has a code that emulates Chemkin's ignition delay sensitivity that may be of use to you,

Best,

Nick

[Weiqi Ji]

Hi Shrabanti,

Recently we have shown that the sensitivity vector of ignition delay time is parallel to the local temperature sensitivity vectors at ignition. For instantly, you can simply use the local temperature sensitivity at ignition to represent the relative sensitivity of the ignition delay time. For the details, you can find a pre-print in the Github repo https://github.com/jiweiqi/IgnSens. I also uploaded a demo using Jupyter notebook. The work will appear in the comping Proceedings of the Combustion Institute.

Here is a quich link to the pdf file: https://github.com/jiweiqi/IgnSens/blob/master/Manuscript/PCI_Sensitivity_R2_June_24_Clean.pdf

Regarding computing the local sensitvities, you can find details in the cantera menu, it should be something like:

sim = ct.ReactorNet([r])
	for i in range(gas.n_reactions):
		r.add_sensitivity_reaction(i)

Best,

Weiqi

[Bryan W. Weber]

Shrabanti,

Actually, the code that Nick linked to uses Chemkin directly to compute the information of interest.  However, you need to define for yourself first what do you mean by ignition delay sensitivity? How do you want to calculate that? You'll need to go look through the literature for ways of calculating it and which would be best for you. The one that I've found most useful is

S = ln(tau_+/tau_-)/ln(k_+/k_-)

where the + and - indicate where you've increased or decreased the reaction rate of the reaction, respectively. This is more useful than, say, a percentage change, because ignition delays are inherently non-linear (they follow an exponential trend). Cantera includes the useful set_multiplier method to increase or decrease the reaction rate of a given reaction; then you can just run a for loop to calculate the ignition delay (again, you need to define that) and compute S as defined above. Make sure you reset the multiplier to 1 after each loop iteration.

The sensitivity1.py example computes effectively the partial derivatives of the target (in this case, OH) with respect to a given reaction. This is useful, but "ignition delay" is not one of the targets we can put in there (in part because the definition varies by application).

Best,

Bryan

[Weiqi Ji]

Hi Bryan,

We have recently shown the correlation between the sensitivity of ignition delay and the local composition sensitivity at ignition. Details can be found in my previous reply. I am wondering if I could make it a Cantera example. Do you think it is necessary to write an API like:

def GetIgnSens(gas, T, p, Composition, 'UV/HP'):

     # first, compute the local temperature sensitivity at ignition point.

     # second, figure the sign of the sensitivity coefficient.

     return sensitivity vector

Or there are better ways for a demo.

Best,

Weiqi

[Shrabanti Roy]

Thanks everyone for your valuable replies.

Based on the information I tried to calculate sensitivity for both ways (following the equation given by Weiqi Ji in his manuscript and using just the temperature). In both cases I am having some troubles. 

1. If I use the equation stated in the manuscript, is that mean, I have to calculate ignition delay time at each iteration?

In my case I am failing to get the updated ignition delay for each iteration. 

2. If I use sensitivity analysis based on temperature using r.add_sensitivity_reaction(i)  , Do I have to use sim.sesnsitivity() to get the sensitivity value?

I tried this approach where I am getting error with CVOde integral error. I might doing something wrong with the code.

I attached both the codes that I tried. Can I get some guideline with these codes? 

Thanks

Shrabanti

[Weiqi Ji]

Hi Shrabanti,

I have written a function to compute the sensitivity in brute force way and another one with local sensitivity. You can find it in https://github.com/jiweiqi/IgnSens/blob/master/demo/ign_sens_excel.ipynb.

Sorry for the poor comment in the code.

I am out of the office. The quick observations of your second python file is that:

1) Do you need to fix the time step, if not, you can switch it to sim.step().

2) You can obtain the temperature sensitivity during computing the ignition delay. You don't have to call sim.sesnsitivity() explicitly.

Weiqi

[Weiqi Ji]

For the question 1, you have to re-calculate the ignition delay after perturbing the rate constant. Instead of the way you are doing in the attached code. Actually, you only have to perturb the top two sensitive reactions according to the manuscript.

Best,

Weiqi

On Sunday, July 8, 2018 at 5:32:12 AM UTC+8, Shrabanti Roy wrote:

[Weiqi Ji]

There is an example code for computing the local sensitivity during autoignition: http://cantera.org/docs/sphinx/html/cython/examples/reactors_sensitivity1.html?highlight=sim%20sensitivity. 

i.e., the sensitvity1.py from www.cantera.org

The code is for OH sensitivity, you can change it to temperature sensitivity. Although they are in the same direction according to the manuscript.

On Sunday, July 8, 2018 at 5:32:12 AM UTC+8, Shrabanti Roy wrote:

[Weiqi Ji]

Sorry, you have to call sim.sensitvity() to get the sensitivity.

[Nick Curtis]

>Bryan W. Weber

>Actually, the code that Nick linked to uses Chemkin directly to compute the information of interest.

whoops! my mistake, didn't read close enough

Nick

[Bryan W. Weber]

Hi Weiqi,

If you already have a Jupyter Notebook, you can make a PR for the Cantera/cantera-jupyter repository on GitHub. The reason we include functions like that in the example scripts is so that we can use the if __name__ == '__main__': ... pattern to avoid running costly functions if the module is imported.

Best,

Bryan

[VISHAL GAWLI]

Hello Shrabanti,

Did you find out the final solution to the problem which you faced while writing this code??

[John D]

Bryan:

In the formula for S (S = ln(tau_+/tau_-)/ln(k_+/k_-)) wouldn't the denominator be a fixed value based on the value of dk chosen? For instance if dk = 0.01 (1 %) then k+ would be 1.01 k_nominal and k- would be 0.99 k_nominal (if k_nominal is the unperturbed rate constant). Am I missing something basic?

Thanks

JD

[Bryan Weber]

Hi JD,

You're exactly right, the denominator only needs to be calculated once.

Best,

Bryan

[John D]

Thanks, Bryan.

Perhaps it would be better to have something like a "central differencing" formula, S = (Tau+ - Tau-)/(2 * Tau_o * dk).  This is similar to the formula used in the flame-speed sensitivity example (which is one-sided differencing formula).

JD

[Bryan Weber]

Hi JD,

The reason I've used the logarithmic formula in the past is because ignition delay times vary over several orders of magnitude in general (that is, they follow the Arrhenius law). As such, a linear formulation of the sensitivity will tend to overemphasize the reactions that make ignition delay longer. Another way to look at it is that the logarithmic formulation makes it so that a change from 1 -> 10 s in ignition delay to have the same sensitivity value as 0.1 -> 1s. The same isn't true for flame speeds, which tend to vary in a fairly narrow range. That said, you could certainly apply a logarithmic formula to the central difference. I'm not sure that it would meaningfully change your analysis, though.

Incidentally, I have to give credit to Bill Pitz at LLNL who first explained this to me :-)

Best,

Bryan

[John D]

Thanks, Bryan.  Your note is helpful.

Is there a quick way to verify that gas.set_multiplier(1 + dk, m) for reaction 'm' is working.

I tried to create an Arrhenius object and print it out as

print("A = ",arr.pre_exponential_factor) # A

before and after applying the gas.set_multiplier factor

but got the same value.

I was just trying to check the change in ignition delay by perturbing one important reaction manually as a test and I saw no change.

I thought I would try out your formula and see what I get.

Just wanted to ensure that I was not doing anything silly

Thanks

JD

[Bryan Weber]

Hi JD,

The multiplier is considered as part of the total reaction rate calculation, it does not modify the pre-exponential factor of the reaction object in memory. You can think of it as a vector of factors that's multiplied with the vector of computed rates. If you share your script, I can help debug what might be going on. Depending on the value of dk, you might not be making much of a change. Can you start a new thread to focus discussion on your script?

Thanks!

Bryan

[John D]

Thanks, Bryan.  I did realize that the set_multiplier does what you mentioned in your post after my own post.  I do see the expected behavior in the results with a change in the dk value (same as hardwiring A*(1+dk)).

Thanks for all your posts.

JD