"Simple Reactor" is not isobaric - neglects dV/dt ?

[rwest]

As part of our solvation corrections project, @bslakman and I were trying to figure out how to change the reactor simulations to work for liquid phase (incompressible) equation of state instead of ideal gas, but it looks like SimpleReactor has been neglecting the volume change (due to stoichiometry differences of reactions) anyway. Did we overlook something, or is this currently missing?

Example: if A <=> 2B then at a conversion of 0.5 (NA/NA0 = 0.5) for isothermal isobaric gas phase reactor, CA/CA0 = 0.5/1.5=0.333. (rxn_rate = d(VC)/dt != (dC/dt)/V)

[connie]

If this is really the case this may solve our problem that the
concentration traces were not matching up with chemkin's solver.

[rwest]

I heard a rumour that you might want to look at this, @connie, as you have recently implemented analytical Jacobians here?...

[rwest]

I think how we resolve this may impact other parts of the code that we were beginning to look at - for example the calculation of conversion (for termination criteria and reporting) currently renormalizes by ideal gas law.
I guess my question is, should the ode solver track N concentrations (in which case we need to add the dilution term to each residual) or should it track N mole numbers and an N+1th scalar representing volume (in which case we need to calculate concentrations from mole number every time we want a rate) ?

[rwest]

Ah, I see @jwallen has just implemented something on his fork in 2224dff.
I should have looked at that first. Does this change @connie's Jacobians at all? (where are these anyway?)

[jwallen]

Yes, I haven't pushed to official yet because we wanted to compare the new implementation against Chemkin first. I think it shouldn't be too hard to swap in the incompressible liquid equation of state for the ideal gas law when the time comes.

@connie's jacobian implementation is in the SimpleReactor.jacobian() method, which currently occupies the majority of rmgpy/solver/simple.pyx.

[rwest]

The proposal is that the ode solver vector y now represents number of moles?
Normalized so that we start with one mole total?
This may make a difference in terms of absolute precision of ode solver? I guess it's not a bad way to normalize - at least its not too hard to explain what ATOL is this way.

[rwest]

I think this is closed by:

• 2224dff
• 36340fd
• 2e60641

(assuming the Jacobians are correct)
Should we add documentation somewhere to describe what y and y0 represent?