Cantera / cantera

Chemical kinetics, thermodynamics, and transport tool suite
https://cantera.org
Other
615 stars 348 forks source link

Bug with the computation of [M] with real gases #967

Closed rhouim closed 3 years ago

rhouim commented 3 years ago

We have been looking through the Cantera source code to examine how the activity concentrations are implemented in preparation for us to add a real-gas equation of state model.

When looking through GasKinetics.cpp and ThirdBodyCalc.h we noticed what could be a bug when computing the real-gas third body concentration. On line 63 of GasKinetics.cpp the activity concentrations are stored in m_conc and on line 64 the normal total concentration is stored in ctot. The activity concentrations, m_conc, and total molar density, ctot, is fed into 3-body update function.

On lines 36-44 of ThirdBodyCalc.h, [M] is calculated in the update function. Note: m_conc is named conc in this function. These lines of code perform:

[M] = m_default*ctot + sum(m_eff*conc_i),

where (presumably) m_default = 1 and m_eff = actual efficiency factor - 1. The portion of [M] computed using ctot does not include real gas effects, whereas the other portion includes the real gas effects. Written another way using c_i for molar concentration and a_i for the activity concentration for species i:

[M] = sum(c_i + (efficiency_factor_i - 1)*a_i)

This seems inconsistent with real gas thermodynamics. I am also pretty sure that, in general, sum(c_i) != sum(a_i).

It seems to me that ctot needs to be redefined in the computation of [M] as sum(a_i).

Perhaps this is a bug that is left over from optimizing the computation of [M] for ideal gas mixtures?

speth commented 3 years ago

Using the NonIdealShockTube.py example as a point of comparison, the modification in #968 changes the calculated ignition delay from 4.093e-04 s to 4.081e-04 s (-0.3%), compared to an ignition delay of 4.333e-4 s calculated using the ideal gas model. So the impact of this correction is pretty modest, at least in this case.

decaluwe commented 3 years ago

Thanks, @rhouim for bringing this up.

The implementation certainly seems inconsistent--one term ignores the activity coefficients, while the other considers them--but I'm wondering if the fix should be to remove activity coefficients from the latter term, not to add them to the former (just to be clear, neither term completely ignores real gas effects, since the EoS is used to calculate the molar density).

I've never really used or looked into third-body reactions, but my naive understanding is that the "third body" [M] is simply a physical object, providing energy via collision. Is this correct? If so, I would think that the collision rate, which is a function of concentration and temperature, would be the relevant term, and that the activities, which represent energy/thermodynamic interactions, would be irrelevant.

speth commented 3 years ago

That's an interesting thought, @decaluwe. In that case, there would still be an odd asymmetry (if not really an inconsistency) in calculating the rate as the product of two activities and one concentration.

I guess part of the issue here is that assuming that any third body has an equivalent effect is already a pretty significant simplification, so it's hard to think about what "should" be used here.

rhouim commented 3 years ago

These a very good points.

In our own codes, we have had the same internal debate: wether to base the computation of [M] with the activity concentrations or the normal molar concentrations? It is clear for what to do with the other concentration terms to get proper equilibrium. However, [M] is a bit ambiguous since it doesn't affect equilibrium.

Nevertheless, it might be more consistent to compute [M] with the activity concentrations. For example, there are occasionally reactions where the third body is explicitly listed in the reaction string itself. E.g., in GRI-3.0

H+O2+N2<=>HO2+N2
H+O2+AR<=>HO2+AR

In these reactions the concentrations of the third bodies, N2 and AR, would likely be computed with the activity concentrations. However, the following reaction

H+O2+M<=>HO2+M

reaction could have the third body concentration computed using normal the molar concentrations. This may be a little inconsistent.

decaluwe commented 3 years ago

So I chatted with @cfgoldsmith and @rwest offline about this, and their take does comport with my naïve understanding, above: the third body [M] provides energy via inelastic collisions to thermalize the reaction. The efficiency with which a collision imparts this energy certainly varies from species to species, but this dependence varies more with things like geometry and molecular structure, which are really separate from the factors that influence the species activity.

You are correct that the species-specific third-body reactions you write above will use the activity concentrations. This seems incorrect, given the paragraph above, although it is academic for ideal gas mechanisms like GRI-3.0 (the activity concentrations and the molar concentrations are the same).

But thinking about it some more, for real gases, if I change the state in a way that changes the activity for some species k but without changing the concentration of that species (for example, by changing the concentrations of other species for which k has a thermodynamic affinity), I don't really see how that would impact the energy imparted by collisions with k.

I suppose the "correct" way to handle this would be to write it as a standard 3-body with a generic [M], but give an efficiency of 1.0 for the species of interest, and a default of 0.0 for all others. Is this currently possible to do in a convenient manner?

Regardless, I don't think we should let the tail wag the dog, here. I would recommend editing #968 so that the molar concentration, rather than the activity concentration, is used for both terms.

rhouim commented 3 years ago

Thanks for the info. Based on this discussion, I agree that from a physical standpoint that treating [M] with normal concentrations may be more "correct".

However, there are some practical issues from a user standpoint. Some reaction mechanisms used in combustion can be quite large with many many species. Users working with supercritical combustion would have to go through all of those reactions and reclassify some of the reactions as third-body reactions and manually setting the all but one of the many enhancement factors (sometimes ~50 to 100 species) to zero. Most users, including myself, would likely (and perhaps unknowingly) eat the inconsistency.

Perhaps a way to handle this, assuming the you folks opt to switch to molar concentrations, that doesn't require users to manually go through and alter mechanism files would be to re-classify reactions with an explicit third-body in the reaction string as a third body reaction. This could, in principle, be performed while the mechanism file is being read and searching for species in normal reactions where deltaN_i (change in moles is zero between the products and reactants) is zero. If there is a species where deltaN_i is zero, then reclassify the reaction as a third body with all of the enhancement factors zero, except for the one species. However, this could potentially cause some confusion if Cantera functions are used to query the properties of a reaction and the classification in Cantera doesn't match what is in the input file.

In general, supercritical and high-pressure combustion is kind of in the Wild West right now. There are lots of open questions and inconsistencies such as using chemical reaction mechanisms developed under ideal gas conditions because nothing else is available. The third body thing is quite a small issue in comparison to these. Thus, either way: normal concentrations or activity concentrations to compute [M] doesn't matter to me. Just so that it is documented and consistent.

To be honest, I only noticed this when I started looking through the supercritical combustion literature for preparation for a course in multiphase combustion and adding a real-gas combustion module into our in-house code. No papers (at least ones that I found) state how they computed [M] for real-gas combustion. So I figured that I would take a look at Cantera and see what you folks did.

We will be using Cantera to perform mechanism reaction reduction to simulate solid-detonation fireballs using the Becker-Kistiakowsky-Wilson EoS model. I will be adding this EoS to Cantera. The compressibility factor is initially quite high (Z ~3 to 5) and we need to reduce large chemical mechanisms to ones that are more practical for multidimensional CFD for our rather unique conditions.

speth commented 3 years ago

Per the suggestion of @decaluwe, @rwest, and @cfgoldsmith, I've changed the implementation in #968 to use the physical concentration when calculating [M] for three body and falloff reactions.

To the question of how to treat reactions like H+O2+N2<=>HO2+N2, @decaluwe's suggestion is possible. In YAML, one would replace:

- equation: H + 2 O2 <=> HO2 + O2
  rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}

with

- equation: H + O2 + M <=> HO2 + M
  type: three-body
  default-efficiency: 0
  efficiencies: {O2: 1.0}

This is maybe not ideal, since now if you print all of the reaction equations, you'll see a bunch of reactions with the equation H + O2 + M <=> HO2 + M, and no indication of why they have different rates. Alternatively, we could do something like what we do for falloff reactions that name a specific third body, which would be to automatically identify the third body based on the stoichiometric coefficients and set up the efficiencies as above when given a definition like:

- equation: H + 2 O2 <=> HO2 + O2
  type: three-body

Then, the last step would be to also identify such reactions when converting from Chemkin input files and label them as three-body reactions.

ischoegl commented 3 years ago

All - I am posting here also as I just pushed a (friendly) amendment to @speth's PR on #968. I believe the question on how to treat concentration has been settled already, but a process to convert reactions that explicitly specify a collision partner is still missing (per @decaluwe's suggestion). As an alternative to what is posted above, it is possible to leave the existing syntax

- equation: H + 2 O2 <=> HO2 + O2
  rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}

and use simple rules to identify this reaction as a three-body reaction in a pre-processing step (similar to what is already done for the detection of electrochemical reactions). I.e. O2 would be detected as the collision partner, and treated as if one would manually specify

- equation: H + O2 + M <=> HO2 + M
  rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}
  type: three-body
  default-efficiency: 0
  efficiencies: {O2: 1.0}

Specifically, I am suggesting the following three rules for converting an elementary to a three-body reaction:

  1. reaction type must not be specified
  2. integer stoichiometric coefficients for all participating species
  3. one side has to have exactly 3 participating species to qualify

There is a draft version in commits to #968 that creates the machinery that makes detection possible, but I wanted to make sure that I am not overlooking anything (i.e. rules are not finalized). Any comments would be appreciated (@speth, @decaluwe, @rwest, and @cfgoldsmith).

rwest commented 3 years ago

I've been following along the various discussions and haven't yet thought of any Gotchyas that haven't already been mentioned. (Hence my silence)

speth commented 3 years ago

For posterity, I wanted to share (with permission) a few comments that Mike Burke made in an email conversation:

I read through the discussion in #967 and #968 and I have little to add beyond (and agree with essentially all of) the main points that eventually found consensus:

  • The role of the third body in both "three-body" and "falloff" reactions is to provide or take away energy via collision.
  • The species dependence of various third bodies comes from their respective efficiencies in transferring energy in a collision and their collision rates (which I suspect correlates more closely to physical concentrations than real-gas activity concentrations, though others on this email chain are better versed in real gas mixtures than I).
  • Based on the above, it probably makes the most sense (to me) to simply use the physical concentrations for the calculation -- at least for now because...
  • Real-gas kinetics is indeed relatively uncharted territory and I would guess that none of the proposed expressions in those issue discussions will quite work because they are still largely premised on ideal-gas-like behavior -- but until more is known about real gas kinetics, there's not much else I can think of doing for real gases (@cfgoldsmith may have more thoughts here)...