Correct handling of third-body concentration in falloff reactions [issue #579]

[yeanment]

This pull request aims to resolve the issues reported in #579 in handling of the three-body falloff reactions of type A + B + C <=> AB + C, where C represents a species like N2.

Specifically, I have made the following updates in my recent two commits:

1. Add an if not (reaction.rate.type == "falloff"): in the block to avoid update of the reactant or product for normal three-body Arrhenius reactions.

if reaction.third_body:        
    if len(reaction.third_body.efficiencies) == 1:
        if isclose(reaction.third_body.default_efficiency, 0.0):
            # Donot update threebody falloff reactions
            if not (reaction.rate.type == "falloff"):
               all_reactants = ...
               all_products = ...
               reagents = copy.deepcopy(...

2. Update the three-body efficiencies for three-body falloff reactions of type A + B + C <=> AB + C.

# three-body-Arrhenius + falloff
efficiencies = reaction.third_body.efficiencies
if (
    falloff 
    and len(reaction.third_body.efficiencies) == 1
    and isclose(reaction.third_body.default_efficiency, 0.0)
):
    for symbol in species_info.all_species_list:
        if symbol not in efficiencies:
            efficiencies[symbol] = 0.0

3. Add statements like cw.writer(fstream, f"alpha = std::max(alpha, 1e-16);") to avoid float-point exceptions. However, I am not sure whether this will have a big influence on the convergencies.

It is noted that I would like to reuse the existing code as much as possible. Therefore, the generated mechanism.H is not optimized, where a better implementation may be available. For instance, in the present implementation, the generated Corr writes like

Corr = mixture + (-1) * sc[0] + (-1) * sc[1] + (-1) * sc[2] + (-1) * sc[3] +
         (-1) * sc[4] + (-1) * sc[5] + (-1) * sc[6] + (-1) * sc[8] +
         (-1) * sc[9] + (-1) * sc[10] + (-1) * sc[11] + (-1) * sc[12] +
         (-1) * sc[13] + (-1) * sc[14] + (-1) * sc[15] + (-1) * sc[16] +
         (-1) * sc[17] + (-1) * sc[18] + (-1) * sc[19] + (-1) * sc[20] +
         (-1) * sc[21] + (-1) * sc[22] + (-1) * sc[23] + (-1) * sc[24] +
         (-1) * sc[25] + (-1) * sc[26] + (-1) * sc[27] + (-1) * sc[28] +
         (-1) * sc[29] + (-1) * sc[30] + (-1) * sc[31] + (-1) * sc[32] +
         (-1) * sc[33] + (-1) * sc[34] + (-1) * sc[35] + (-1) * sc[36] +
         (-1) * sc[37] + (-1) * sc[38] + (-1) * sc[39] + (-1) * sc[40] +
         (-1) * sc[41] + (-1) * sc[42] + (-1) * sc[43] + (-1) * sc[44] +
         (-1) * sc[45] + (-1) * sc[46] + (-1) * sc[47] + (-1) * sc[48] +
         (-1) * sc[49] + (-1) * sc[50] + (-1) * sc[51] + (-1) * sc[52] +
         (-1) * sc[53] + (-1) * sc[54] + (-1) * sc[55] + (-1) * sc[56] +
         (-1) * sc[57] + (-1) * sc[58] + (-1) * sc[59] + (-1) * sc[60];

which is equivalent to

Corr = sc[7];

I'm not familiar with parsing simplification, so I just left it as is.

[baperry2]

Closes #579.

Thanks for submitting this! A bit swamped right now, but having this as a PR will make considering the fixes much easier when time allows.

[yeanment]

Hi, @baperry2, I have revised the implementation for the code generations to specifically address third-body fall-off reactions with a single third-body species.

• Targeted Reactions: Pressure-dependent reactions of the form A + B (+C) = D (+C), such as H + O2 (+N2) = HO2 (+N2).
• Third-body Efficiency: Implemented a logic shift where the default third-body efficiency is set to 0 for all species except the explicitly defined third-body.
• Test Kinetic Model: The recent NH3/H2 model by Zhu et al. published in Combustion and Flame (doi:10.1016/j.combustflame.2023.113239). The mechanism file used for validation is attached as Zhu2024NUIG.zip

---

Implementation Detail

The commit addresses the issues raised above and the code now generates cleaner Corr for mechanism.H compared to previous implementations, as shown below.

---

Validation & Results

To verify the update, I conducted simulations of a planar unstretched flame with PeleLMeX, and the initial profile is computed using a Cantera solution. Both the original and updated implementations were tested under identical setups for a finite time interval.

• Black Line: Initial flame front location.
• Blue Line: Original/previous implementation.
• Red Line: New implementation (this PR).

Conclusion: The new model provides more reliable results, showing a distinct improvement in predicting the flame's progression compared to the original code in dealing with such reactions.