[rhythm]

Hi folks,

I have been looking at bounding constraint implementation in combustion solvers. As known the species scalarFields for these solvers are supposed to be bounded to range [0:1]. The lower bound implementation is clear as it is only needed to clean up numerical error build up, whereas the upper bound is slightly more involved and raises some questions.

As far as I can tell there is only one part in the code that takes care of upper species bound. It can be seen in multiComponentMixture class. There is a member function called correctMassFractions(), which in summary does the following:

template<class ThermoType>
void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
{
    volScalarField Yt("Yt", 1.0*Y_[0]);

    for (label n=1; n<Y_.size(); n++)
    {
        Yt += Y_[n];
    }
<...>
    forAll(Y_, n)
    {
        Y_[n] /= Yt;
    }
}

This member function is invoked during construction of multiComponentMixture object. Which leads me to my questions:

1) As it can be seen the normalisation constraint is only applied to the internalField values of Y_. This does not prevent silly users (like me ) from specifying boundary condition values higher than 1. Wouldn't it make sense to include some sort of fail safe check to prevent users from specifying unphysical values at the boundaries?

2) correctMassFractions() function, as far as I can tell, is only called during construction of multiComponentMixture, which probably makes sense as we wouldn't want to additionally fiddle with the field values at runtime apart from what the solver calculates. Nonetheless, I would expect some sort of procedure telling the code that there is an upper constraint for field values. Unfortunately, I can't seem to find any of that in the code. Does anyone know if I am missing something or are the field values kept in range [0:1] by a correct choice of reaction, thermo parameters?

Cheers,
Ben

[rhythm]

As an example, I broke the reactingFoam code to solve the following equations:





For a trivial reaction and assuming density is constant this can be solved analytically:







NB.: I am aware that reactingFoam takes density as well as temperature changes into account, so the equations that it would be solving, as well as the results would be slightly different. Nonetheless, they wouldn't (or at least they shouldn't) be drastically different

NBB.: are molecular weights of the species and Arrhenius reaction rate.

I ran the test on a simple cavity case with no velocity field and uniform species field values and plotted the values at the cell[0], also plotted the analytical solution.
Other params:
- MolWeights:
- Arrhenius:
- thermodynamics and transport coeffs I picked from reactingFoam tutorial

reactingFoam:
massFractionEvo.png

analytical:
analyticalEvo.png

In addition I broke the code even more and made density const as well as set the temperature as constant which in the end gives identical solution to the analytical:
massFractionEvo_constCase.png

I am aware that breaking the code like that isn't always the best choice when one wants to understand how the algorithm works, but I thought that it is quite interesting to see whether it can be reduced to give an analytical solution.

Anyways, if anyone has some insight as to how species scalarField bounding is implement I would highly appreciate your input.

Cheers,
Ben