custom reaction rate expression

[youk]

Hello all,

I'm trying to implement a custom single-step, global irreversible oxidation reaction (for now propane as a fuel) into cantera, so first let me give a little bit of context:

The single-step reaction uses temperature-dependant kinetic parameters (activation energy and pre-exponential factor) to re-create the (possibly non-monotonic) dependency of auto-ignition delay times on initial mixture temperature. Using a single-step reaction with this Arrhenius-type reaction rate expression:

Rate = (rho * Y_fuel / M_fuel)^(n_fuel) * (rho * Y_O2 / M_O2)^(n_o2) * A * exp(- E_a / (R * T))

rho is the mixture density and A is in 1/s in this expression (if reaction order is 1)

Analitically, it can be shown that using this reaction with reaction orders n_fuel + n_o2 = 1 (I use 0.5 for both), the autoignition behaviour can be recreated, given temperature dependent kinetic parameters. Ea and A are derived from autoignition simulations using detailed mechanisms and our own simple integrater code can reproduce the results from the detailed mechanism very well. (The methodology and applications in LES simulations  are explained in these papers, but don't waste your time on this :) 10.1088/1742-6596/2885/1/012039, https://doi.org/10.1016/j.proci.2024.105235). 

As far as I understand, because of these temperature-dependant kinteic parameters, I have to use a custom reaction rate expression. Based on the custom_reactions.py example and the documentation I have been able to make a working script.

Yet, the results are not  what I would expect. I need to multiply the reaction rate by 12 orders of magnitude to get similar (but not as accurate as it should) results to the detailed mechanism and apart from that, there is substantial deviation at high temperatures which makes me think it is not just an order of magnitude problem.

Sorry for rambling, but my question is:

Is there an obvious mistake in my implementation of the custom reaction and rate in my script? How should I handle potential unit conversion? I don't fully understand it from the documentation. 

Attached is my script, some results from a detailed mechanism (CRECk-polimi) for plotting and an example of my output.

 

Many thanks in advance!

[Ray Speth]

Hi,

What you’re seeing here is that the custom rate expression provides a value for the rate constant which is then (by default) multiplied by the reactant concentrations. Of course, since the concentration dependence is already embedded in the rate expression, that’s not desired here. You can override this default behavior by specifying the “order” of the reaction with respect to those reactants to be zero, before adding the reaction to the Solution:

custom_reaction.orders = {'C3H8': 0, 'O2': 0}

After doing this, you’ll see that both gas.forward_rate_constants and gas.forward_rates_of_progress return the same value for this reaction.

Regards,
Ray

​

[youk]

Thank you so much Ray,

setting the reaction orders to 0 gave errors during integration, but I removed the concentrations dependence from my rate expression and   set

custom_reaction.orders = {'C3H8': 0.5, 'O2': 0.5}
Which works just fine and is a bit more elegant imo.

Cheers!