Adsorption/Desorption reaction rates

[Bjarne Kreitz]

Hi all,

I'm using Cantera to run a transient simulation of a CSTR for a heterogeneously catalyzed reaction, where I alternate the feed between two different mixtures. I use the advance() command to move the reactor forward in time. Everything works just fine and the gas phase concentration and coverage profiles look absolutely reasonable. However, when I compute the net rates of progress, I can see some weird oscillations in the adsorption/desorption steps. These oscillations do not occur for the reactions, which include only adsorbates (see the pdf).
The oscillations can also be observed in the net production rates of the gaseous products.

Does anybody have an idea why these oscillations occur and how to get rid of them? I don't understand the cause of the oscillation. I'm also not sure whether this is an issue since the rest of the results look great. Still, it makes the evaluation of transient simulations with sensitivity analyses or degree of rate control quite challenging. I played with the tolerances, but it did not reduce the noise.

An example with the corresponding .yaml file is attached.

Best,

Bjarne

[Ingmar Schoegl]

Hi Bjarne,

This reminds me of one of the MATLAB examples, which also exhibits periodic behavior: https://github.com/Cantera/cantera/blob/main/samples/matlab/periodic_cstr.m

So I believe steady state behavior is not necessarily to be expected.

-ingmar-

[Ingmar Schoegl]

PS: The corresponding Python example is here https://cantera.org/examples/python/reactors/periodic_cstr.py.html

[Richard West]

Hi Ingmar, 

I don’t think Bjarne’s issue is that it doesn’t reach steady state, but that the net rates of progress are not being calculated precisely, and show a lot of numerical noise - the orange and blue lines on the third and fourth plots that Bjarne shared are bouncing up and down with a lot of imprecision, compared to the smooth green and red lines that are about the same overall magnitude.

My assumption is that the reactions are fast but close to equilibrium and so there is a high forward rate of progress and an almost equally high reverse rate of progress so the net rate of progress is small but imprecise because it’s the difference between two large numbers. But where is the imprecision?  None of the reactants or products have particularly negligible concentrations or coverages (the orange one is H2 reacting with vacant sites to give adsorbed H atoms, all three of which things have non-negligible amounts on the first two plots).

Picking through the code….

InterfaceKinetics::updateMu0()

gets the standard chemical potentials, then 

InterfaceKinetics::updateKc()

gets the equilibrium constant Kc by doing exp(-∆G/RT), (with a warning in the source code about overflow) 

// WARNING this may overflow HKM

m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);

thenceforthe m_rkcn is the reciprocal of Kc.

this is called by  InterfaceKinetics::_update_rates_T()

which is called by InterfaceKinetics::updateROP()

which then calculates the reverse rate 

kr = kf / Kc

m_ropr[i] = m_ropf[i] * m_rkcn[i];

these kf and kr are then multiplied by the activity concentrations to get forward and reverse rates of progress

    // multiply ropf by the activity concentration reaction orders to obtain
    // the forward rates of progress.
    m_reactantStoich.multiply(m_actConc.data(), m_ropf.data());

    // For reversible reactions, multiply ropr by the activity concentration
    // products
    m_revProductStoich.multiply(m_actConc.data(), m_ropr.data());

and finally the net is calculated as the forward minus the reverse

m_ropnet[j] = m_ropf[j] - m_ropr[j];

if the forward and reverse are large and almost equal, but not very precise, you’d get a noisy net rate.

If Kc >> 1   and activity concentration of products >> activity concentration of reactants, or vice versa.

It’s occurring only in the examples of reactions that are adsorption/desorption, where the number of gas phase and surface species changes. 

Either increase the precision at which all these intermediate values are stored, or normalize them (redefine activity concentration for adsorbates?) such that they’re more nearly equal to the activity concentration of gas phase species?   I started playing around seeing if there’s a cunning approximation for (1-x) for small x, or something, but didn’t see an obvious quick win.

Or just give up and say it can’t be solved to the desired precision.

-- 
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/1518d353-422f-4aef-95ae-e5a856342743n%40googlegroups.com.

[Ingmar Schoegl]

Richard,

I agree with your assessment - the oscillations are actually noise. To illustrate the point, I replotted only the forward direction based of Bjarne's example, and everything is smooth. It is still surprising that the data are that noisy, as the relative errors are quite large - there are two and four orders of magnitude separation between what appear to be smooth directional data and the noisy difference.

I did notice a rather unusual absolute tolerance for the simulations; I didn't probe deeper though.

-ingmar-

[Ray Speth]

Hi all,

This ends up being a fairly subtle error, but one where at least related versions have come up before. The issue is that when a Solution object is used inside a ReactorNet, there is no guarantee that the state of the Solution after a call to ReactorNet.advance or ReactorNet.step corresponds exactly to the state of the system at the time indicated by the ReactorNet. The CVODES integrator underlying the ReactorNet just uses the attached Solution objects as “calculators”, and can leave them in arbitrary states that do not correspond to the state of the system at the ReactorNet‘s current time.

One case where this occurs is when using the advance method, where CVODES takes arbitrary internal timesteps, governed only by the absolute and relative tolerances on the solution, and then uses interpolation to determine the state at the specified output time. In this case, the Solution objects will never be automatically set to the state corresponding to the output time, but instead will have a state related to the last time reached by the integrator.

Cantera provides accessors for the Solution objects that will automatically synchronize the state to the correct values determined by CVODES, specifically the Reactor.thermo property for bulk phases and the ReactorSurface.kinetics property for interfaces. This problem is exacerbated in the case of surface phases, where the reaction rates are a function of both the bulk and surface states, which may end up in an inconsistent state if only one of the two is correctly synchronized. The recommendation I would make that once a Solution object has been used to construct a Reactor or ReactorSurface, it should not be accessed directly anymore, but only accessed via these properties of the reactor classes, to ensure that the state is set correctly.

A more general solution that would prevent such an error from being made would possibly require changing a feature of the current reactor network model, which is that a single Solution object may be used to represent multiple reactors within a reactor network. However, given the number of times that problems like this have arisen, re-evaluating this concept may be worth it.

Regards,
Ray

​

[Ingmar Schoegl]

Dear Ray,

Thanks for your insights here. This makes perfect sense!

-ingmar-

[Bjarne Kreitz]

Hi Ray,

Thanks again for solving this puzzle.

Best,

Bjarne