Sensitivity analysis

[tjabbo]

Hi people,

I want to do a sensitivity analysis for the combustion of ammonia with respect to temperature, pressure and equivalence ratio. That means that I am looking for the most sensitive 10 reactions with respect to temperature,pressure and equivalence ratio. For the beginning I used the sensitivity script of cantera but there are some lines which I dont get what they are doing. For example in the script they define tolerances for the solution and for the sensitivty coefficient:

sim.rtol = 1.0e-6
sim.atol = 1.0e-15
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6

What do these tolerances mean?

Furthermore I dont quite understand how I can set the range of temperature, pressure or equivalence ratio. My understanding is that in order to calculate a sensitivity with respect to temperature one has to vary the temperature in order to see how the different reaction rates will change. Does anyone know where I can set these ranges?

This is the sensitivty script of cantera:

gas = ct.Solution('gri30.yaml')
temp = 1500.0
pres = ct.one_atm

gas.TPX = temp, pres, 'CH4:0.1, O2:2, N2:7.52'
r = ct.IdealGasConstPressureReactor(gas, name='R1')
sim = ct.ReactorNet([r])

# enable sensitivity with respect to the rates of the first 10
# reactions (reactions 0 through 9)
for i in range(10):
    r.add_sensitivity_reaction(i)

# set the tolerances for the solution and for the sensitivity coefficients
sim.rtol = 1.0e-6
sim.atol = 1.0e-15
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6

states = ct.SolutionArray(gas, extra=['t', 's2', 's3'])

for t in np.arange(0, 2e-3, 5e-6):
    sim.advance(t)
    s2 = sim.sensitivity('OH', 2)  # sensitivity of OH to reaction 2
    s3 = sim.sensitivity('OH', 3)  # sensitivity of OH to reaction 3
    states.append(r.thermo.state, t=1000*t, s2=s2, s3=s3)

    print('{:10.3e} {:10.3f} {:10.3f} {:14.6e} {:10.3f} {:10.3f}'.format(
        sim.time, r.T, r.thermo.P, r.thermo.u, s2, s3))

Thank you for your answer!

Kind regards,

Tjabbo

[Bryan Weber]

Hi,

Sorry for the delayed response! The tolerances you note in the beginning of your post control the numerical accuracy required by the integrator before moving to the next time step.

The default sensitivity analysis is a local sensitivity calculated by numerically perturbing the parameters of the model and finding the change in the specified model solution variable. So, for example, this line:

sim.sensitivity('OH', 2)

finds the sensitivity of the concentration of OH (one of the solution variables) with respect to perturbing the reaction rate of reaction 2 (one of the model parameters) at the current time step.

To answer your question, we need to know what you're trying to calculate to find the sensitivity. For example, if you're calculating the ignition delay, computing the sensitivity of the ignition delay with respect to the reaction rates cannot be done directly, because the ignition delay is not one of the solution variables. Instead, you'll need to define a separate formula for how you want to compute that sensitivity. Please reply back with any questions and with some more details about what you want to do, and we'll be happy to help more.

Hope that helps!

Bryan

[tjabbo]

Hello Bryan,

Thank you for your answer.

I want to compute the sensitivity of the reaction rates with respect to the temperature, equivalence ratio and pressure. So for example I want to know which reaction rates are strongly influenced by a change in the temperature. I have added the script which I am currently using in the attachment (largely based on the script found in https://skill-lync.com/student-projects/week-9-senstivity-analysis-assignment-18).

However when I am running the script it gives me following error:

TypeError: sensitivity() takes at least 2 positional arguments (1 given)

This is due to the line:

sens  = sim.sensitivity(l[range_val-1]),r.component_name(param_index) #sensitivity of the different reactions with respect to the temperature

The only thing which I have changed from the original script is that I have swapped the parameters because i want to know how the temperature influences the reaction rates and not the other way round. Do you know what I am doing wrong and how I could fix this?

Furthermore, I dont really know how I could change the parameter temperature to another parameter like pressure or equivalence ratio.

Thank you in advance!

Tjabbo Reinhardus

[Steven DeCaluwe]

Dear Tjabbo,

From the sensitivity method documentation, it looks like the function only does sensitivities in one direction: it will calculate the sensitivities of the solution variables with respect to an input parameter (e.g. reaction rate coefficient, species enthalpy, etc.)

To find the production rate sensitivities to those variables, there would be two possible paths:

1. Calculate them directly from theory.  

Given that the reaction rate is computed according to one of the formulas listed in Cantera’s science section, you don’t actually need to run the reactor to find the sensitivity (if I am understanding you correctly), but rather could derive them by hand.

2. Define your own custom sensitivity routine.

That said, I can see the potential convenience of letting the software calculate the sensitivities. In this case, you would want to use finite differencing and calculate the sensitivities yourself, using the formula given here, for example.

As an example, calculating sensitivity w/r/t temperature might look something like this:

import cantera as ct

gas = ct.Solution('gri30.yaml')

T_o = 1500 # Default temperature, K

dT = 1 # Temperature perturbation, K

gas.TP = T_o, ct.one_atm

rop_o = gas.net_rates_of_progress # Rxn rates at default conditions

gas.TP = T_o+dT, None # perturb the temperature

rop_new = gas.net_rates_of_progress # Rxn rates at perturbed state.

# Add a small number to avoid possible division by zero)

rxn_sens_T = T_o*(rop_new - rop_o)/dT/(rop_o + 1e-24)

Note this script returns mostly zeros since H2 is the only species present when I load gri 3.0 like this. 

I do have broader questions about why you want these sensitivities in the context of a reactor simulation, since T, P, and equivalence ratios are results / outcomes, not parameters / inputs, and the reaction rates are process variables, not observable results.  In other words: finding the sensitivity of a reaction rate to the temperature at a given time in the reactor has limited use, since the temperature is a dependent variable in this case, not something one can control / select.

Hope this helps,

Steven

——————————————————

Steven C. DeCaluwe, Ph.D | Associate Professor of Mechanical Engineering

COLORADOSCHOOLOFMINES

Brown Building W410B

Golden, CO 80401

Email: deca...@mines.edu

Phone: (303) 273-3666

Twitter: @CORESresearch

https://people.mines.edu/decaluwe/

https://cores-research.mines.edu

He / Him / His

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/60f04252-3557-4ab7-9295-c82a703acfabn%40googlegroups.com.
<Sensitivty Analysis Script.txt>

[tjabbo]

Hello Steven,

Thank you for your answer! To come back to your question on why I am interested in the sensitivity of the reaction rate:

I am currently working on an engine running on ammonia and I am analysing the emission of this engine. In our case we can control the combustion temperature by for example increasing the oxygen concentration (more ammonia is burned). This temperature increase results in a reduction of the N2O emissions. Now I want to further research this reduction of N2O by looking at reaction paths and their dependence on the temperature.  So for example if I conduct a sensitivity analysis and I see that the reaction rates which will form N2O are reduced with increasing temperature, this would confirm the experimental results. Hope this clarifies on why I want to conduct this sensitivity analysis.

Do you think this is a good approach?

Kind regards,

Tjabbo