Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support Chemkin parameterization for multiple bath gas species #193

Open
ljxsbest opened this issue Jan 9, 2024 · 3 comments
Open

Support Chemkin parameterization for multiple bath gas species #193

ljxsbest opened this issue Jan 9, 2024 · 3 comments
Labels
feature-request New feature request

Comments

@ljxsbest
Copy link

ljxsbest commented Jan 9, 2024

The Chemkin format is not supported using the latest Cantera version. For example, in our mechanism, the format is:

H+O2(+M)=HO2(+M)    4.66E12 0.44 0.0E0 
    HE/0.57/
    AR/0.0/
    N2/1.0/
    O2/1.0/
    H2/2.0/
    CH4/2.0/
    CO2/3.25/
    H2O/17.6/
    CO/4.0/
    LOWMX  /4.0662E19 -1.4E0 -1.80537E2/
    TROEMX /5.0E-1 1.0E0 1.0E10 1.0E30/
    LOWSP  /N2 1.91E+20 -1.5568 253.86/
    TROESP /N2 5.0E-1 1.0E0 1.0E10 1.0E30/

However, the format is not supported in the latest Cantera version, and we must extract the N2 component and change the format into two separate reactions:

H+O2(+M)=HO2(+M)    4.66E12 0.44 0.0E0
    HE   /0.57/
    AR   /0.0 /
    N2   /0.0 /
    O2   /1.0 /
    H2   /2.0 /
    CH4  /2.0 /
    CO2  /3.25/
    H2O  /17.6/
    CO   /4.0 /
    LOW  / 4.0662E19 -1.4E0 -1.80537E2 /
    TROE / 5.0E-1 1.0E0 1.0E10 1.0E30  /

H+O2(+N2)=HO2(+N2)    4.66E12 0.44 0.0E0   
    LOW  / 1.91E+20 -1.5568 253.86    / 
    TROE / 5.0E-1 1.0E0 1.0E10 1.0E30 /

When we try to adjust the format to fit the Cantera, it’s not working by Chemkin.
Could you please try to fix this issue in your side to support the Chemkin format in the next version?

@ljxsbest ljxsbest added the feature-request New feature request label Jan 9, 2024
@speth
Copy link
Member

speth commented Jan 9, 2024

I haven't seen these Chemkin keywords before: LOWMX, TROEMX, LOWSP, TROESP. I have several questions here:

  • Is there any public documentation of what these keywords mean and how the resulting rate coefficient is supposed to be calculated?
  • When were these keywords introduced?
  • Is the second input form actually equivalent, or is the formula for calculating the rate coefficient more involved? I don't understand why additional syntax would have been added if the second form is equivalent.
  • What do you mean that the second input form isn't working? That looks like syntax that would have been parseable all the way back to Chemkin-II.

I'm wondering if this ends up being a special form of Mike Burke's mixing rules (see #157), which @pjsingal is currently working on implementing in Cantera.

@rwest
Copy link
Member

rwest commented Jan 9, 2024

It's new to me too. I'm pretty sure fair-use doctrine allows me to post this excerpt of a Release 2021R2 technical manual that I just found online:

3.6.4. Usage of Multiple Bath-gas Species
The typical usage of the reaction types outlined in Unimolecular/Recombination Fall-off Reactions (p. 36) through General Pressure Dependence Using Logarithmic Interpolation (p. 40) is with a specific species present as the "bath" gas. Many times the reactions for different bath-gas species are all included in the same mechanism. Note that the reaction rate computed at the high or low pressure limit would, in theory, be identical regardless of which species acts as the bath gas. However, if multiple instances of the same reaction with different bath-gas species are present in a given mechanism, then the net rate reaction rate is computed by summing up the individual rates. If each rate expression gives the reaction rate at the limiting pressure as $R$ and suppose that there are $n$ such rate expressions (for $n$ different bath-gas species), then the total rate is $n*R$ (due to the summation). Note, however, that the expected value is $R$. To accommodate such cases, an enhanced way of specifying the reaction with multiple bath-gas species is provided. In particular, the rate expressions for different bath-gas species may be written together in the same reaction. For reactions written in this manner, the net rate is computed by summing the individual rates multiplied by the mole fraction of corresponding bath-gas species (or the mixture). Thus, the correct limiting rate is obtained.

Consider the following example:

H +O2 (+M) = HO2 +(M) 4.650e+012 0.440 0.0
LowMX/1.737e+019 -1.230 0.0/
TroeMX/ 6.700e-001 1.000e-030 1.000e+030 1.000e+030/ LowSP/ AR 6.810e+018 -1.200 0.0/
TroeSP/AR 7.000e-001 1.000e-030 1.000e+030 1.000e+030/ LowSP/ HE 9.192e+018 -1.200 0.0/
TroeSP/HE 5.900e-001 1.000e-030 1.000e+030 1.000e+030/ HE / 1.0 / AR / 1.0 / H2/ 1.30/ H2O/ 10.00 /

Note that the keywords LOW and TROE (see the Chemkin-Pro Input Manual) here are augmented by SP and MX; which indicate the rate expressions for a particular species (such as Argon (AR) and Helium (HE) in this case) and the rest of the mixture. The net rate for this reaction is then
$R = R_{AR}X_{AR} + R_{HE}X_{HE} + R_{MIX}X_{MIX}$ (3.38)
where $R$ and $X$ indicate rate and mole fraction, respectively, and
$X_{MIX} = 1 - (X_{AR} + X_{HE})$
is the mole fraction of the "rest" of the mixture. In the example above, each rate expression makes use of the third-body efficiencies. For the specific bath-gas species (that is 'AR' and 'HE' in this example) non-unity third-body efficiency cannot be specified.

This formulation is also available when (a) the reaction line input is for the low- pressure limit, that is, for individual species and for a mixture, the rates are indicated by the keyword HIGHMX and HIGHSP and (b) for general logarithmic pressure dependence use keywords PLOGMX and PLOGSP as shown in the example below:

H + O2 = HO2         1.0 0. 0.
   PLOGMX/1.      2.48E+17 -2.226  36. /
   PLOGMX/10.     1.90E+18 -2.194   6. /
   PLOGMX/100.    6.05E+19 -2.330 635. /
   PLOGSP/HE 1.   2.48E+17 -2.226  36. /
   PLOGSP/HE 10.  1.90E+18 -2.194   6. /
   PLOGSP/HE 100. 6.05E+19 -2.330 635. /

For PLOG input, third-body enhancements are not allowed. In the example above, the rates for the specific species (HE in this case) and for the rest of the mixture are first computed by appropriate interpolation for the specified pressure. The net rate is then computed by multiplying by the corresponding mole fractions as shown in Equation 3.38 (p. 42).

@speth speth changed the title Mechanism compile issue Support Chemkin parameterization for multiple bath gas species Mar 19, 2024
@speth
Copy link
Member

speth commented Nov 13, 2024

Thanks, @rwest. Now that the linear-Burke reaction rate type is implemented, which handles multiple bath gas species, I think we have a possible way of handling this, with one major caveat. Taking the first example from the Chemkin docs, we could use the same parameters to write a linear-Burke reaction:

reactions:
- units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol}
- equation: H + O2 (+M) <=> HO2 (+M)
  type: linear-Burke
  duplicate: true
  colliders:
  - name: M
    type: falloff
    low-P-rate-constant: [1.737e+019, -1.230, 0.0]
    high-P-rate-constant: [4.650e+012, 0.440, 0.0]
    Troe: {A: 6.700e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
  - name: AR
    type: falloff
    low-P-rate-constant: [6.810e+018, -1.200, 0.0]
    high-P-rate-constant: [4.650e+012, 0.440, 0.0]
    Troe: {A: 7.000e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
    efficiency: {A: 1, b: 0, Ea: 0}
  - name: He
    type: falloff
    low-P-rate-constant: [9.192e+018, -1.200, 0.0]
    high-P-rate-constant: [4.650e+012, 0.440, 0.0]
    Troe: {A: 5.900e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
    efficiency: {A: 1, b: 0, Ea: 0}
  - name: H2
    efficiency: {A: 1.3, b: 0, Ea: 0}
  - name: H2O
    efficiency: {A: 10.0, b: 0, Ea: 0}

(@pjsingal, let me know if you agree this interpretation is correct; I think the translation is all species use the same high-pressure rate constant, and the efficiency for those species with a specific rate parameterization is 1)

The caveat is that the linear-Burke mixing rule is different from the mole fraction weighted sum that is implemented by Chemkin, and almost universally more accurate than that ad hoc approach for combining the rates (see Burke and Song, 2017). In case the kinetic data is based on best available data for each collider, then the linear-Burke form should give better results, and there's no reason not to use it. However, if someone has done some global optimization of their mechanism and rate coefficients based on the mole fraction weighted sum approach, then switching the way the rates are evaluated could lead to worse results. On the other hand, I'm not really keen on implementing another reaction rate type just to provide a less accurate way of representing such reactions.

For a more complete example, I've put together a YAML file and small script for comparing the three different approaches for the above reaction, and consistent with the previous literature result, you can see that using the more accurate mixing rule can have a fairly large effect on the rate constant, here compared as the ratio of the linear-Burke rate to either the "separate reactions" approach or the "mole fraction weighted sum" approach.

ckpdep.yaml
phases:
- name: gas
  thermo: ideal-gas
  kinetics: bulk
  species:
  - {gri30.yaml/species: [O2, H, AR, N2, H2O, HO2, H2]}
  - {nasa_gas.yaml/species: [He]}
  state: {T: 300, P: 1 atm, X: {O2: 1.0, H: 1.0e-5, AR: 4.4, HE: 2.0, N2: 4.0, H2O: 0.1}}

# H +O2 (+M) = HO2 +(M) 4.650e+012 0.440 0.0
# LowMX/1.737e+019 -1.230 0.0/
# TroeMX/ 6.700e-001 1.000e-030 1.000e+030 1.000e+030/ LowSP/ AR 6.810e+018 -1.200 0.0/
# TroeSP/AR 7.000e-001 1.000e-030 1.000e+030 1.000e+030/ LowSP/ HE 9.192e+018 -1.200 0.0/
# TroeSP/HE 5.900e-001 1.000e-030 1.000e+030 1.000e+030/ HE / 1.0 / AR / 1.0 / H2/ 1.30/ H2O/ 10.00 /

reactions:
- units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol}
# Reaction 0: Combined rate using Burke reduced pressure linear mixing rule
- equation: H + O2 (+M) <=> HO2 (+M)
  type: linear-Burke
  duplicate: true
  colliders:
  - name: M
    type: falloff
    low-P-rate-constant: [1.737e+019, -1.230, 0.0]
    high-P-rate-constant: [4.650e+012, 0.440, 0.0]
    Troe: {A: 6.700e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
  - name: AR
    type: falloff
    low-P-rate-constant: [6.810e+018, -1.200, 0.0]
    high-P-rate-constant: [4.650e+012, 0.440, 0.0]
    Troe: {A: 7.000e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
    efficiency: {A: 1, b: 0, Ea: 0}
  - name: He
    type: falloff
    low-P-rate-constant: [9.192e+018, -1.200, 0.0]
    high-P-rate-constant: [4.650e+012, 0.440, 0.0]
    Troe: {A: 5.900e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
    efficiency: {A: 1, b: 0, Ea: 0}
  - name: H2
    efficiency: {A: 1.3, b: 0, Ea: 0}
  - name: H2O
    efficiency: {A: 10.0, b: 0, Ea: 0}

# Reaction 1-3: Separate reactions for colliders with unique rate constant expressions
- equation: H + O2 (+M) <=> HO2 (+M)
  type: falloff
  duplicate: true
  low-P-rate-constant: [1.737e+019, -1.230, 0.0]
  high-P-rate-constant: [4.650e+012, 0.440, 0.0]
  Troe: {A: 6.700e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
  efficiencies: {AR: 0, He: 0, H2: 1.3, H2O: 10.0}
- equation: H + O2 (+AR) <=> HO2 (+AR)
  type: falloff
  duplicate: true
  low-P-rate-constant: [6.810e+018, -1.200, 0.0]
  high-P-rate-constant: [4.650e+012, 0.440, 0.0]
  Troe: {A: 7.000e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
- equation: H + O2 (+He) <=> HO2 (+He)
  type: falloff
  duplicate: true
  low-P-rate-constant: [9.192e+018, -1.200, 0.0]
  high-P-rate-constant: [4.650e+012, 0.440, 0.0]
  Troe: {A: 5.900e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}

# Reaction 4-6: Separate reactions, but calculated on the full mixture concentration.
# These are the rates needed to compute the Chemkin-style mole-fraction based weighting
- equation: H + O2 (+M) <=> HO2 (+M)
  type: falloff
  duplicate: true
  low-P-rate-constant: [1.737e+019, -1.230, 0.0]
  high-P-rate-constant: [4.650e+012, 0.440, 0.0]
  Troe: {A: 6.700e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
  efficiencies: {H2: 1.3, H2O: 10.0}
- equation: H + O2 (+M) <=> HO2 (+M)
  type: falloff
  duplicate: true
  low-P-rate-constant: [6.810e+018, -1.200, 0.0]
  high-P-rate-constant: [4.650e+012, 0.440, 0.0]
  Troe: {A: 7.000e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
- equation: H + O2 (+M) <=> HO2 (+M)
  type: falloff
  duplicate: true
  low-P-rate-constant: [9.192e+018, -1.200, 0.0]
  high-P-rate-constant: [4.650e+012, 0.440, 0.0]
  Troe: {A: 5.900e-001, T3: 1.000e-030, T1: 1.000e+030, T2: 1.000e+030}
ckpdep.py
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.constrained_layout.use'] = True

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


def kf_burke(gas):
    return gas.forward_rate_constants[0]

def kf_separate(gas):
    return sum(gas.forward_rate_constants[1:3])

def kf_chemkin(gas):
    XM = 1 - sum(gas['AR','HE'].X)
    XAR, XHE = gas['AR','HE'].X
    return gas.forward_rate_constants[4:] @ [XM, XAR, XHE]


burke = []
separate = []
chemkin = []

for i in range(10000):
    T = 300 + 1000 * np.random.random()
    P = 10**(-2 + 10 * np.random.random())
    X = np.random.random(gas.n_species)
    gas.TPX = T, P, X
    burke.append(kf_burke(gas))
    separate.append(kf_separate(gas))
    chemkin.append(kf_chemkin(gas))

burke = np.array(burke)
separate = np.array(separate)
chemkin = np.array(chemkin)

fig, ax = plt.subplots(1, 2, figsize=(9,5))
ax[0].plot(burke / separate, '.')
ax[1].plot(burke / chemkin, '.')
ax[0].set(ylabel=r'$k_{\text{Burke}} / k_{\text{separate}}$')
ax[1].set(ylabel=r'$k_{\text{Burke}} / k_{X\text{-weighted}}$')
plt.show()

burke-ratios

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
feature-request New feature request
Projects
None yet
Development

No branches or pull requests

3 participants