Equilibrium chemistry calculations
 

Hi,

I'm working on a custom solver based on buoyantBoussinesqSimpleFoam.
I'd like to introduce a very simple equilibrium chemistry module, mostly as a preliminary/benchmarking activity for future coupling with a dedicate Gibbs Energy Minimization tool.

At the moment, I'd like to simulate precipitation/dissolution between two transported species.
The specific equilibrium model does not matter for now, as I'm more interested in how to perform the actual coupling between the OF solver and a generic equilibrium chemistry model.
Let's assume that equilibrium and mass conservation are prescribed as the algebraic constraint



where both species concentrations and are also subject to separate transport equations.
In the literature, I've found several examples of coupling between CFD and equilibrium chemistry codes.
The common approach is to perform sequential segregated transport and reaction steps, and iterate until proper convergence is achieved.

What I'm failing to understand is how is it possible to transport separately different species and at the same time prescribe some algebraic equilibrium relation without over-constraining the overall numerical system.
Since only equilibrium is involved, reaction rates disappear from transport equations.
This means that each chemistry iteration redistributes somehow local mass between the species to enforce equilibrium, causing the "transport" part alone to never reach numerical convergence (e.g. based on residuals).

What am I missing? Is there some trivial aspect I am not considering?
If someone can give some advice or share some useful reference it would be greatly appreciated.

Thanks,
Andrea