Fail to converge for zeroD simulation

[BangShiuh]

Problem description This problem usually occurs for a reaction mechanism with reverse reactions. The reverse reaction rates can become very large for some conditions which might be caused by using the mechanism outside of the temperature range.

Steps to reproduce error from running:

import cantera as ct
import numpy as np

ct.suppress_thermo_warnings()
gas = ct.Solution("https://raw.githubusercontent.com/BangShiuh/reaction-mechanisms/main/input/A2NOx.yaml", transport_model=None)

solutions = ct.SolutionArray(gas, extra=["t"])
air_composition = 'N2:3.76, O2:1.0'
fuel_composition = "POSF10325:1.0"

phi = 0.1
gas.set_equivalence_ratio(phi, fuel_composition, air_composition)
gas.TP = 300, ct.one_atm
inlet = ct.Reservoir(gas)

# initiate ideal gas reactor
gas.equilibrate('HP', estimate_equil=-1)
reactor = ct.IdealGasReactor(gas)
reactor.volume = 1.0

def mdot(t):
    return reactor.volume * reactor.density / 1e-3

inlet_mfc = ct.MassFlowController(inlet, reactor, mdot=mdot)
exhaust = ct.Reservoir(gas)
ct.PressureController(reactor, exhaust, master=inlet_mfc, K=1e-5)
sim = ct.ReactorNet([reactor])

rt = 0.001
solutions.append(reactor.thermo.state, t=sim.time)

while sim.time < rt * 10:
    # try:
        sim.step()
        solutions.append(reactor.thermo.state, t=sim.time)
    # except:
    #     break

# solutions.write_hdf("result-issue.h5", group='phi')

To see the large reverse reaction constant:

import cantera as ct
import matplotlib.pyplot as plt

ct.suppress_thermo_warnings()
gas = ct.Solution("/home/bschen/reaction-mechanisms/input/A2NOx.yaml")
gas.TP = 300, ct.one_atm

plt.plot(gas.reverse_rate_constants, ".")

print(gas.reaction(1153))
print(gas.reverse_rate_constants[1153])

After making this reaction irreversible, the simulation converged.

Behavior

Traceback (most recent call last): File "phi-issue.py", line 34, in sim.step() File "build/python/cantera/reactor.pyx", line 1158, in cantera.reactor.ReactorNet.step cantera._utils.CanteraError:

---

CanteraError thrown by CVodesIntegrator::step: CVodes error encountered. Error code: -3

At t = 0.0012841 and h = 2.1008e-11, the error test failed repeatedly or with |h| = hmin. Components with largest weighted error estimates: 196: -215.00762894680796 197: 210.9919998140404 21: 5.7857926262174235 2: -0.008801723314520293 23: 0.0030090242342010846 4: 1.0254549364550695e-06 6: 7.229548501420058e-07 14: 5.002437405334156e-07 12: 4.414005225560132e-07 8: 3.0764574908068557e-07

---

System information

• Cantera version: 3.0.0

Attachments

Additional context

[BangShiuh]

In GasKinetics.cpp,

void GasKinetics::updateKc()
{
    thermo().getStandardChemPotentials(m_grt.data());
    fill(m_delta_gibbs0.begin(), m_delta_gibbs0.end(), 0.0);

    // compute Delta G^0 for all reversible reactions
    getRevReactionDelta(m_grt.data(), m_delta_gibbs0.data());

    double rrt = 1.0 / thermo().RT();
    for (size_t i = 0; i < m_revindex.size(); i++) {
        size_t irxn = m_revindex[i];
        m_rkcn[irxn] = std::min(
            exp(m_delta_gibbs0[irxn] * rrt - m_dn[irxn] * m_logStandConc),
            BigNumber);
    }

    for (size_t i = 0; i != m_irrev.size(); ++i) {
        m_rkcn[ m_irrev[i] ] = 0.0;
    }
}

Can we use a "smaller" BigNumber to solve this problem? Or maybe there is a more physical way to set the limitation to the equilibrium constant.

[BangShiuh]

The reaction mechanism can be downloaded at https://raw.githubusercontent.com/BangShiuh/reaction-mechanisms/main/input/A2NOx.yaml