Concentration Sensitivity Analysis of laminar premixed flames

[Fanne]

Hi all !

I am currently working with the Python interface of Cantera, and performing 1D premixed flame calculations of methane oxidation with the GRI30 (irreversible) mechanism.

To perform a scheme reduction, I would like to get the sensitivities of the concentrations of the solution of premixed flame calculation to the reaction rate coefficients of each reactions, in the same fashion as explained here : https://groups.google.com/forum/#!topic/cantera-users/qMOxAdolJiA
Now, since from my understanding -but I would reeaally like to stand corrected - there is no way of doing this directly as in CHEMKIN, I tried to systematically modify the multiplier of each reaction in order to recompute my flame (surprisingly enough, it doesn't take that much time to do), and to use the formula that can be found both in the previous post and in the premix.f file of CHEMKIN (II) to get the normalized sensitivities of each species.
My problem is that it does not confront well with CHEMKIN results. First of all, the orders of magnitude are way off (sensitivities in the order of 10e4 with cantera ...??). So I guess that my way of performing these calculations is wrong, but I cannot find a way to correct it without finding out :
what exactly is that "multiplier" in cantera, is it comparable to the result of CHEMKIN's subroutine CKRDEX ? (but even in CHEMKIN I always wondered what it did in reality)
why, when normalizing the sensitivities, do we use the multiplier value and not the reaction rate coef ?

Thank you for your help !


---------------------------------------------
Just in case, here is my sensitivity calculation routine (similar to the one given in the post I mentionned) :

fout = open('out.out','w')

T0 = f.T()
II=gas.nReactions()
KK=gas.nSpecies()
JJ=f.flame.nPoints()

#Sensitivity analysis------------------------------------------

#Similar to the DP of CHEMKIN
dk = 2.98e-8

#LOOP sur les reactions
for i in range(II):
    fout.write('         ' + str(i+1))
    gasbis.setMultiplier(1+dk, i)
    fbis.solve(loglevel, 0)
    T = fbis.T()
#LOOP sur les mesh points
    for n in range(JJ):
        f.setGasState(n)
        fbis.setGasState(n)
        G0 = gas.moleFractions()
        G = gasbis.moleFractions()
#LOOP sur Temp, especes, enrg
        SensiTemp=(T[n]-T0[n])/(T0[n]*dk)
        fout.write('  ' + str(SensiTemp))
        for k in range (KK):
            Su0=G0[k]
            Su=G[k]
            Sensi=(Su-Su0)/(Su0*dk)
            fout.write('  ' + str(Sensi))
    gasbis.setMultiplier(1, i)
    fout.write("\n")
fout.close()

[Ray Speth]

Fanne,

There are a couple of things that come to mind that might help you. 

First, the "multiplier" in Cantera is used to scale the forward rate of progress for the corresponding reaction. So if the rate of progress for a reaction is Rf based on the rate expression and the multiplier is M, the actual rate reported will be M*Rf. The Chemkin documentation for the CKRDEX function is not exactly enlightening, so I'm not sure how that function compares.

Second, You need to choose the size of the perturbation "dk" carefully. Obviously you don't want it to be so large that the you get a nonlinear response. However, you also need it to be large compared to the tolerances that you're specifying for the 1D solver (which aren't shown in the code you provided), otherwise the perturbation will be lost in the noise.

Finally, since it looks like you're trying to compute sensitivities at particular spatial locations in a freely propagating flame, it's important to consider how the constraint on the flame location will affect your computed sensitivities. In both Chemkin and Cantera, the flame position is fixed by setting the temperature to a fixed value at a single point along the domain. What value of the temperature is used for this will affect the resulting sensitivity coefficients. The fact that an artificial parameter introduced for the numerical solution changes your sensitivity coefficients suggests to me that these aren't really physically interesting sensitivities, and that you would be better off looking at other integral properties which do not depend on the choice of the temperature fixed point.

Regards,

Ray

[Fanne]

Ray,

Thank you for your reply, and for the definition of the multiplier, I thought that is what it did.
Concerning the sensitivities at particular spatial locations, if you are right in your explanation, it doesn't really matter; because when I post process I am always comparing two points with the same temperature, and not so much the same position. But in any way, such a small perturbation of the multiplier should not be giving me sensitivities of the order 10⁴, it's just too much !
You are right about the numerical value of the perturbation: it was smaller than the tolerances of my solver. Could it explain the order though ?

[Bryan W. Weber]

Ray,

You specify that the multiplier is used to scale the forward rates of progress. From the code in GasKinetics.cpp, the scaling of the forward rates is done before the calculation of the reverse rate of progress, and thus, for reversible reactions (at least for the gas, interface, and aqueous phases), the reverse rate is multiplied by the same factor as well. Correct?

Fanne,

If you grep (or find) in your CHEMKIN-II cklib source for "PAR(4" you can find all the places the perturbation factor is used. From my examination of the CHEMKIN-III cklib, it seems to serve a similar purpose as the multiplier in Cantera. However, I've never had any luck getting the CKRDEX keyword to work as I expected.Alternatively, you can use CKRAEX to get/set the Arrhenius coefficient of the reaction rates. For all the reactions except pressure dependent ones, this is equivalent to multiplying the rate of progress. Note that this will not be the case for the pressure dependent reaction formulations, such as Troe or SRI. You may also be able to modify the subroutine CKQC, which calculates reaction rates of progress in CKLIB.

Bryan

[Ray Speth]

Bryan,

Yes, the reverse rate of progress is also scaled by the multiplier. Doing otherwise would be inconsistent with the definition of the equilibrium constant.

Fanne,

An approximate way to think about the tolerance issue is this: With a solver tolerance of B and a perturbation dk, you should expect the change in a solution variable Y to be on the order of dY = Y*(B + dY/dk * dk). You're trying to find dY/dk by dividing this by dk, but if B is larger than dY/dk * dk, all you're going to get is B/dk.

For the spatial position issue, if your goal is to look at the sensitivity of the concentrations at a point where the temperature has a particular value, I think you have to do the transformation from Y(z) to Y(T) before you compute the sensitivity coefficient.

Regards,

Ray

[Fanne]

Thank you all for your help, it seems that changing the tolerance of my solver solved my current issues !
I am going to investigate this chemkin perturbation factor though, because i would love to find a way to optimise the calculation with cantera; in deed, it turns out that when the tolerance of the solver is smaller than the perturbation, it takes over two days to perform a 1D flame calulation with sensitivities (computed at all points for all the species, but still).
has anyone tried to perform sensitivities calculations ''on the fly'', with the use of the jacobian for example directly ?

[Katharina]

This thread was really helpful for me. Thanks a lot to everyone! What I still do not understand, is why you disable the "refine grid" option for solving the scaled reaction set. Will the solver still use the refined grid produced in the first, not manipulated, run? Will the scaling of forward rate of progress not also influence local gradients in the solution?

[Ray Speth]

Katharina,

With the refine grid option disabled, the solver will continue to use the same already-refined grid. The reason for disabling grid refinement is because we want to calculate the sensitivities for the particular discretized problem that we originally solved. Grid refinement introduces a lot of nonlinearities into the system, since you can only add or remove an integer number of grid points, and this is done on the various thresholds. If you left grid refinement enabled, you might find that perturbing a particular reaction rate would trigger the addition of a grid point at some location, which would change the computed flame speed, and what you would be calculating is the sensitivity with respect to adding that grid point rather than just the sensitivity to that reaction rate. As long as the perturbations are small, the original refined grid should still be valid for the perturbed problem.

The local gradients in the solution do of course change when the reaction rates are perturbed, which is why the perturbed system needs to be re-solved for each parameter.

Regards,

Ray

[Katharina]

Ray,

thanks for your reply, that was very helpful again. It is all working now.

Best regards,
Katharina