Field bounding in combustion solvers

rhythm · November 2, 2016, 11:13

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:

Code:

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 · November 2, 2016, 11:19

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

$\frac{d (\rho Y_i)}{d t} = W_i R_i$

$R_i =\mp K_f \prod_j C_j$

For a trivial reaction $A \rightarrow B$ and assuming density is constant this can be solved analytically:

$\frac{d Y_i}{d t} = \mp K_f \frac{W_i}{W_A} Y_A$

$Y_A = Y_{A}(0) \exp(- K_f t)$

$Y_B = Y_{A}(0) \frac{W_B}{W_A} \Big[1-\exp(-K_{f} t)\Big]$

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.: $W_i; K_{f};$ 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: $A=1; \beta=0 ; Ta = 0$
- 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

November 2, 2016, 11:13	Field bounding in combustion solvers	#1
rhythm Member Ben Jankauskas Join Date: Jun 2013 Location: Exeter Posts: 34 Rep Power: 12	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: Code: template<class ThermoType> void Foam::multiComponentMixture<ThermoType>::correctMassFractions() { volScalarField Yt("Yt", 1.0Y_[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

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Foam::error::PrintStack	almir	OpenFOAM Running, Solving & CFD	91	December 21, 2022 04:50
Schools for Combustion (CFD) Ph.D.?	gdbb89	Main CFD Forum	5	July 21, 2016 20:36
XiFoam lean-premixed combustion instabilities	remir	OpenFOAM Running, Solving & CFD	8	April 15, 2016 08:31
Bounding k and epsilon & coldEngineFoam	sasanghomi	OpenFOAM	1	September 13, 2013 12:12
Velocity Field in Combustion	Guosheng Kang	Main CFD Forum	5	October 11, 1999 20:55

November 2, 2016, 11:19		#2
rhythm Member Ben Jankauskas Join Date: Jun 2013 Location: Exeter Posts: 34 Rep Power: 12	As an example, I broke the reactingFoam code to solve the following equations: $\frac{d (\rho Y_i)}{d t} = W_i R_i$ $R_i =\mp K_f \prod_j C_j$ For a trivial reaction $A \rightarrow B$ and assuming density is constant this can be solved analytically: $\frac{d Y_i}{d t} = \mp K_f \frac{W_i}{W_A} Y_A$ $Y_A = Y_{A}(0) \exp(- K_f t)$ $Y_B = Y_{A}(0) \frac{W_B}{W_A} \Big[1-\exp(-K_{f} t)\Big]$ 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.: $W_i; K_{f};$ 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: $A=1; \beta=0 ; Ta = 0$ - 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