Ignition Delay Time

[Jiawei Wan]

Hi there,

When I tried to calculate IDT, I noticed that the source code given by the official website (https://cantera.org/examples/jupyter/reactors/batch_reactor_ignition_delay_NTC.ipynb.html) used the ".set_equivalence_ratio()" method to define the composition of the gas. The ".X()" method should also be able to achieve the same result.

When I compared the two methods, the results were quite different. At the same time, I changed other yaml files and calculated the IDT of other fuels, and there were also many differences.

Has anyone encountered a similar situation or knows about this?

".set_equivalence_ratio()" method: gas.set_equivalence_ratio(phi=1.0, fuel="nc7h16", oxidizer={"o2": 1.0, "n2": 3.76})

".X()" method: gas.X = {'nc7h16': 0.018740629685157422, 'o2': 0.20614692653673164, 'n2': 0.775112443778111}

Thank you so much for your time!

Best regards,

Jay

[Anant Girdhar]

Hi Jay, 

Based on a quick experiment I just did, it seems like using using gas.X changes the state (the pressure caught my attention) and that's what's causing the difference in IDT. Here are three cases that I tried running with the results:

# TRIAL 1: Initialize the Solution object using set_equivalence_ratio():

gas = ct.Solution("seiser.yaml")
reactor_temperature = 1000  # Kelvin
reactor_pressure = 101325  # Pascals
gas.TP = reactor_temperature, reactor_pressure
# Define the fuel, oxidizer and set the stoichiometry

gas.set_equivalence_ratio(phi=1.0, fuel="nc7h16", oxidizer={"o2": 1.0, "n2": 3.76})

Result: Computed Ignition Delay: 3.248e-02 seconds. Took 0.44s to compute

# TRIAL 2: Initialize the Solution object using X():

gas = ct.Solution("seiser.yaml")
reactor_temperature = 1000  # Kelvin
reactor_pressure = 101325  # Pascals
gas.TP = reactor_temperature, reactor_pressure
# Define the fuel, oxidizer and set the stoichiometry

gas.X = {"nc7h16": 0.018740629685157422, "o2": 0.20614692653673164, "n2": 0.775112443778111}

Result: Computed Ignition Delay: 3.461e-02 seconds. Took 0.46s to compute

# TRIAL 3: Use set_equivalence_ratio() but set TP after that:

gas = ct.Solution("seiser.yaml")
reactor_temperature = 1000  # Kelvin
reactor_pressure = 101325  # Pascals
# Define the fuel, oxidizer and set the stoichiometry

gas.set_equivalence_ratio(phi=1.0, fuel="nc7h16", oxidizer={"o2": 1.0, "n2": 3.76})

gas.TP = reactor_temperature, reactor_pressure  # DOING THIS STEP AFTER SETTING PHI

Result: Computed Ignition Delay: 3.248e-02 seconds. Took 0.44s to compute

# TRIAL 4: Use X() but set TP after that:

gas = ct.Solution("seiser.yaml")

reactor_temperature = 1000  # Kelvin
reactor_pressure = 101325  # Pascals
# Define the fuel, oxidizer and set the stoichiometry

gas.X = {"nc7h16": 0.018740629685157422, "o2": 0.20614692653673164, "n2": 0.775112443778111}

gas.TP = reactor_temperature, reactor_pressure  # DOING THIS STEP AFTER SETTING COMPOSITION

Result: Computed Ignition Delay: 3.248e-02 seconds. Took 0.45s to compute

So Trials 3 and 4 give the same results. You can also see that by printing out the state of the gas (by executing gas()). Here's the output you get from there:

TRIAL 1: Set TP and then use set_equivalence_ratio()

  gas:

       temperature   1000 K
          pressure   1.0132e+05 Pa
           density   0.36789 kg/m^3
  mean mol. weight   30.188 kg/kmol
   phase of matter   gas

                          1 kg             1 kmol
                     ---------------   ---------------
          enthalpy        7.1647e+05        2.1629e+07  J
   internal energy        4.4105e+05        1.3314e+07  J
           entropy              8152        2.4609e+05  J/K
    Gibbs function       -7.4355e+06       -2.2447e+08  J
 heat capacity c_p            1316.7             39749  J/K
 heat capacity c_v            1041.3             31434  J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
                n2           0.71929           0.77511           -25.105
                o2           0.21851           0.20615           -28.134
            nc7h16          0.062207          0.018741           -92.735
     [ +157 minor]                 0                 0

TRIAL 2: Set TP and then use X()

  gas:

       temperature   1000 K
          pressure   94027 Pa
           density   0.3414 kg/m^3
  mean mol. weight   30.188 kg/kmol
   phase of matter   gas

                          1 kg             1 kmol
                     ---------------   ---------------
          enthalpy        7.1647e+05        2.1629e+07  J
   internal energy        4.4105e+05        1.3314e+07  J
           entropy            8172.6        2.4672e+05  J/K
    Gibbs function       -7.4561e+06       -2.2509e+08  J
 heat capacity c_p            1316.7             39749  J/K
 heat capacity c_v            1041.3             31434  J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
                n2           0.71929           0.77511            -25.18
                o2           0.21851           0.20615           -28.209
            nc7h16          0.062207          0.018741           -92.809
     [ +157 minor]                 0                 0

TRIAL 4: Use X() and then set TP

  gas:

       temperature   1000 K
          pressure   1.0132e+05 Pa
           density   0.36789 kg/m^3
  mean mol. weight   30.188 kg/kmol
   phase of matter   gas

                          1 kg             1 kmol
                     ---------------   ---------------
          enthalpy        7.1647e+05        2.1629e+07  J
   internal energy        4.4105e+05        1.3314e+07  J
           entropy              8152        2.4609e+05  J/K
    Gibbs function       -7.4355e+06       -2.2447e+08  J
 heat capacity c_p            1316.7             39749  J/K
 heat capacity c_v            1041.3             31434  J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
                n2           0.71929           0.77511           -25.105
                o2           0.21851           0.20615           -28.134
            nc7h16          0.062207          0.018741           -92.735
     [ +157 minor]                 0                 0

So you can see that TRIAL 2 produces a different state than the other cases which is probably what's affecting the IDT calculation.

I hope this helps at least a little! 

anant

--
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/a427cd56-53c2-4db1-b6c2-3345a3d76f6fn%40googlegroups.com.

[Jiawei Wan]

Dear Anant,

thank you so much for your wonderful suggestions. Based on your idea, I made a more comprehensive test and found some kind of strange results. I would like to share with you and hear your opinions or views.

The Reaction_mechanism comes from

https://www.sciencedirect.com/science/article/pii/S0010218021003965

The yaml format I have attached.

There are some major conclusion:

The IDT employed by the method .set_equivalence will not be influenced by the type of mechanism or initial temperature.

While, the IDT calculated by method .X() is affected  by initial T or mechanisms. Even if the gas statue showed the same, the calculated IDTs were different.

Thank you for your time!

Best Regards,

Jay

[Ray Speth]

Hi Jay and anant,

Hopefully I can demystify this a bit.

As the documentation for the set_equivalence_ratio method notes, it sets the equivalence ratio while holding T and P constant. In contrast, just setting the mole fractions alone with gas.X = ... hold the other intrinsic state variables of the phase constant. For the ideal gas phase, these are temperature and density, which means that setting the mole fractions will change the pressure. If you want to hold T and P constant while setting mole fractions, you can use syntax like the following:

gas.TPX = None, None, 'CH4:0.8, O2:0.2'

Setting gas.TP = ... after using gas.X = ... also works, of course.

Regards,
Ray

​

[Jiawei Wan]

Hello Ray,

Thank you so much for your explanation about the set_equivalence_ratio method, it makes me very clear. Thanks again.

 

While I attempted to try whether I could obtain the same results by using the set_equivalence_ratio method and set gas.X first and then setting gas.TP method, the results were not consistent.

Taking the well-known gri30 as example, when the temperature is 1000K the results is consistent, while increasing the temperature to 1600K or 2000K, it starts showing the discrepancy. (set the phi=1)

Taking the heptane mechanism ( the 160 species mechanism by Seier et al. 2000, Proc. Comb. Inst.) as example. (Similarity T = 1000K, the two methods are the same, while T = 1600K, the two results difference)

What is puzzling is that the outputs of gas statue were the same, the calculated IDTs were different.

Thanks and Best regards,

Jay

[Ray Speth]

Hi Jay,

Could you please provide the code used for these examples? There's not enough information here to tell why you're getting different results.

Regards,

Ray

[Jiawei Wan]

Hello Ray,

Thanks for the reply. And the attached is the test code.

Best,

Jay

[Ray Speth]

Hi Jay,

What’s happening here is that the results of setting X, then T and P versus setting T and P then setting the equivalence ratio result in very small differences in the state, on the order of the rounding error for floating point calculations (relative differences of about 2e-16). However, during the time integration of the ignition delay problem, there are many threshold conditions, such that the solver may not take identical time steps even for conditions that are so close together. Further, the definition of the ignition delay used here (from some of the Cantera samples) is somewhat coarse — it can only be one of the discrete output times of the integrator.

One thing you can do to make this result more consistent would be to at least keep all the output from the integrator by changing the integration loop and removing the “if” condition here:

while t < estimated_ignition_delay_time: t = reactor_network.step() if not counter % 10: # We will save only every 10th value. Otherwise, this takes too long # Note that the species concentrations are mass fractions time_history.append(r.thermo.state, t=t) counter += 1

I think that unless you have a very large mechanism, the admonition here about saving all the state data "taking too long" is probably not a significant concern.

Regards,
Ray