Consolidate reactor models and extend sparse solver applicability

[speth]

Abstract

The aim of this work is to consolidate the different implementations of governing equations used for zero-dimensional reactors, and to allow the sparse/preconditioned solver to be used with non-ideal equations of state.

Motivation

The inheritance diagram for reactor models currently looks like this:

Many of these models are conceptually identical in terms of the physical system they represent, and only differ in terms of the state variables used in the model. Specifically, we have two dimensions where we have two choices of state variable:

1. Implementation of the energy equation in terms of energy (internal energy or enthalpy depending on the physical constraint) or temperature
2. Use of mass fractions or (total) moles for species

By making two specific improvements to a couple of these models, we can eliminate the other nearly-redundant model options, leaving us with:

• Reactor implementing a model with volume, temperature, and species moles as the state variables, supporting any compressible equation of state
• ConstPressureReactor implementing a model with temperature, pressure, and species moles as the state variables, supporting any equation of state
• ReactorSurface implementing a model with surface species moles as the state variables, supporting any surface equation of state.
• FlowReactor / FlowReactorSurface - TBD

Description

Generalize the temperature-dependent energy equations to account for non-ideal fluids

Currently, the temperature-dependent reactor formulations such as IdealGasReactor assume an ideal gas equation of state when rewriting the energy equation to change the state variable from $U$ to $T$, starting with the definition:

$$U = m \sum_k Y_k u_k(T)$$

and differentiating it to obtain:

$$\frac{dU}{dt} = u \frac{dm}{dt} + m c_v \frac{dT}{dt} + m \sum_k u_k \frac{dY_k}{dt}$$

However, it is also possible to expand a generic definition of $U = U(T, V, m_1, \ldots, m_K)$ where $m_k = mY_k$:

$$\frac{dU}{dt} = \frac{\partial U}{\partial T}\Bigg|_{V,m_k}\frac{dT}{dt} + \left.\frac{\partial U}{\partial V}\right|_{T,m_k}\frac{dV}{dt} + \sum_k \frac{\partial U}{\partial m_k}\Bigg|_{T,V,m_{j\ne k}}\frac{dm_k}{dt}$$

Here we can make the following substitutions:

$$ \left(\partial U / \partial T \right)_{V,m_k} = m c_v $$

$$\pi_T \equiv (\partial U / \partial V)_{T,m_k}$$

where $\pi_T$ is the internal pressure, and

$$(\partial U / \partial m_k)_{T,V,m_{j\ne k}} = \tilde{u}_k /W_k $$

where $\tilde{u}_k$ are the partial molar internal energies at constant temperature and volume. This gives:

$$ \frac{dU}{dt} = m c_v \frac{dT}{dt} + \pi_T \frac{dV}{dt} + \sum_k \frac{\tilde{u}_k}{W_k}\frac{dm_k}{dt} $$

This modification has been implemented in Cantera/cantera#2094, covering the IdealGasMoleReactor and IdealGasReactor classes, along with implementations of these new thermo properties for relevant phase types.

Something similar (and I believe even a bit simpler) can be done for constant pressure reactors.

Improve scaling / tolerances for mole-based reactors

Generally speaking, the mole-based reactors perform well compared to the mass-fraction based variants, often requiring somewhat fewer timesteps during integration. These are also the reactor types where we have preconditioners implemented (for the ideal gas variants) due to the simpler form of the derivatives when using moles as the state variable. As such, there doesn't seem to be a strong reason to keep the mass fraction-based versions around.

One minor issue that I think should be addressed before making a direct substitution of reactor types is to provide a better (dynamic) default for the integrator absolute tolerances for these reactors. In the mass fraction-based case, the state variables are already scaled, and the integrator absolute tolerance is already a nondimensional value. However, in the case of the mole-based reactors, the scale of the species moles varies depending on the reactor volume, temperature, and pressure. As such, for a very small reactor or at very low pressure, an absolute tolerance of 1.0e-15 may actually be fairly large compared to the species mole numbers for some species that actually do need to be resolved. Conversely, for large volumes or high pressures, this could be an unnecessarily tight tolerance.

As such, I think what we should do for these reactors is set the default absolute tolerance based on the initial conditions of the reactor so it represents an absolute tolerance of 1.0e-15 on the species mole fractions. Some open questions are (a) whether we should update this value at all as the simulation progresses, in case there is a large change in the problem's scale and (b) whether we should likewise scale user-specified values for the absolute tolerance in such a case. My initial inclination is to keep it simple and go with "no" on both of these fronts.

Rename models and delete other versions

• Rename IdealGasMoleReactor to Reactor
• Rename IdealGasConstPressureMoleReactor to ConstPressureReactor
• Rename MoleReactorSurface to ReactorSurface
• Delete the remaining redundant reactor types

This is obviously a breaking change, so the goal is to complete this as part of the Cantera 4.0 major version jump.

Alternatives

The other option I considered for unifying some of the different models was to switch to having the equation of state be one of the governing equations, creating a DAE system where both $U$ and $T$ are state variables. I thought this might turn out to have some nice features, but unfortunately the integration was substantially slower due to additional Jacobian evaluations that seem to be required by the Sundials IDAS solver.

References

• Reconsider composition state variables for reactors, flames, and other systems #154
• Addition of Jacobian Preconditioning Functionality to Cantera. #47