Reference pressures: bars or atmospheres?

[Nicholaswogan]

I have built a photochemical model which has a design partially inspired by Cantera: https://github.com/Nicholaswogan/photochem

I've been benchmarking my kinetics against Cantera's. I load in the same reactions file into Cantera, and my photochemical model, then I compute rates and reverse rates for reactions .Everything looks identical, except the reverse rates of progress.

For example

# Already created a cantera `gas` object with identical composition to my photochem model
i = 0
print(gas.reverse_rates_of_progress[i]*(ct.avogadro/1e6))

The result from Cantera is 2.4172383036990714 and the result from my model is 2.3856287231176796. The ratio of these two numbers is 2.4172383036990714/2.3856287231176796 = 1.01325. In other words, the number of bars in a standard atmosphere.

The difference ultimately arises when reverse rates are computed. The formula is

k_r = (k_f/exp(-DG_r / (R*T))) * (k_B*T/P_ref)**(n_p - n_r)

k_r is the reverse rate constant, k_f is the forward rate constant, DG_r is the standard gibbs energy of the forward reaction, k_B is boltzman constant, T is temperature, P_ref is the reference pressure, n_p is the number of products in the forward reaction, and n_r is the number of reactants in the forward reaction.

In my model I have been assuming that P_ref = 1 bar. I believe in Cantera, P_ref = 1 atm. I believe this because when I change P_ref in my code to 1 atm, then my reverse rates of progress from my model perfectly match Canteras (to with machine precision).

I used P_ref = 1 bar because of the text above Equation 9 in this paper: https://iopscience.iop.org/article/10.1088/0004-637X/738/1/72/meta . Perhaps this is wrong?

My question: Should P_ref 1 bar or 1 atm? Can you give a citation to justify this? Note that I'm using the Shomate polynomials from NIST for gibbs free energies. Perhaps this is relevant.

[ischoegl]

@Nicholaswogan ... Thanks for submitting the issue. Could you please share the YAML input file used for Cantera? (or at least the portion that specifies the reaction you tested?)

[Nicholaswogan]

Attached is the yaml file. Also attached is a pickle file that contains temperature, pressure and a composition. You can then run the following script to print the forward and reverse rate of the first reaction (falloff), which reproduces what I describe in the original post.

import numpy as np
import cantera as ct
import pickle

with open('TPX.pkl','rb') as f:
    TPX = pickle.load(f)
    
gas = ct.Solution('zahnle_earth_1_ct.yaml')
T = TPX['T']
P = TPX['P']
X = TPX['X']
gas.TPX = T, P/1.0e1, X

i = 0
print(gas.reaction_equations()[i])
print(gas.forward_rates_of_progress[i]*(ct.avogadro/1e6))
print(gas.reverse_rates_of_progress[i]*(ct.avogadro/1e6))

# My model gets identical forward rate,
# but the reverse rate is a factor of 1.01325 smaller.

TPX.pkl.zip
zahnle_earth_1_ct.yaml.zip

[speth]

Cantera's default is to use 1 atm as the reference pressure, though you can change it. This is documented for the YAML format here. Also, you can check from Python:

>>> import cantera as ct
>>> gas = ct.Solution('gri30.yaml')
>>> print(gas.reference_pressure)
101325.0
>>> print(gas.species('H2').thermo.reference_pressure)
101325.0

NIST appears to specify 1 bar as the reference pressure for the data they provide in the Chemistry WebBook (e.g. CH4), so if you're using data from there, you should include setting the reference-pressure in your YAML species definitions.

I would hesitate to change the default be 1 bar for species that use the Shomate model, but this caveat should probably be mentioned where we suggest NIST as a source of data.

Lack of clarity about reference pressures in various data sources has been a thorny issue for a while. See #336 for reference.

[Nicholaswogan]

@speth thanks. This resolves my question!

[Nicholaswogan]

@speth This does not quite resolve my issue. I tried changing the reference pressure in my Cantera input yaml file by setting reference-pressure to 1 bar for all species. The result of the following is 1 bar.

print(gas.species('He').thermo.reference_pressure)

But the result of the following is one atmosphere

print(gas.reference_pressure)

When I compute reverse rates, it seems like nothing has changed - the result still appears to assume a reference pressure of 1 atm for all species.

[speth]

Thanks, @Nicholaswogan, for the additional testing. It seems that the user-specified reference pressures are indeed being ignored. I've prepared a change that should fix this problem.

[Nicholaswogan]

I see this issue should be resolve in v3.0.0. But I'm still having issues setting the reference pressure to 1 bar. How should this be done?

[speth]

You should add reference-pressure: 1 bar to the thermo section for each species in the YAML input file. If that's not working for you, please provide a full example, including your input file. Everything I've tried seems fine, including the fact that the reverse rate constants change by a factor of 1.01325 for reactions with a change in the net number of moles.

It wouldn't hurt to double check that you're actually using Cantera 3.0.0 by examining the value of ct.__version__, too.