[sharifi]

Hi all,

I'm currently working with XiFoam-Solver (OF 2.4) and I'm trying to understand how exactly the temperature is calculated from the total enthalpy. Actually I'm looking for a reference where it is explained clearly. It would be great If you could share such a reference with me if you are aware of any.
As far as I understood, in XiFoam we are dealing with one step chemistry
reactants --> products
So it is possible to calculate the Nasa polynomial constants using the data for let's say methan-air combustion using equivalence ratio and "thermoData" of OpenFOAM.

The temperature is calculated by calling
thermo.correct(); from EaEqn.H

Because the XiFoam uses:

thermoType
{
    type            heheuPsiThermo;
    mixture         homogeneousMixture;
    transport       sutherland;
    thermo          janaf;
    equationOfState perfectGas;
    specie          specie;
    energy          absoluteEnthalpy;
}

The function correct(); can be found in "heheuPsiThermo.C" as:

template<class BasicPsiThermo, class MixtureType>
void Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::correct()
{
    if (debug)
    {
        Info<< "entering heheuPsiThermo<BasicPsiThermo, MixtureType>::correct()"
            << endl;
    }

    // force the saving of the old-time values
    this->psi_.oldTime();

    calculate();

    if (debug)
    {
        Info<< "exiting heheuPsiThermo<BasicPsiThermo, MixtureType>::correct()"
            << endl;
    }
}

The function calculate can be found in the same h-data:

template<class BasicPsiThermo, class MixtureType>
void Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::calculate()
{
    const scalarField& hCells = this->he_.internalField();
    const scalarField& heuCells = this->heu_.internalField();
    const scalarField& pCells = this->p_.internalField();

    scalarField& TCells = this->T_.internalField();
    scalarField& TuCells = this->Tu_.internalField();
    scalarField& psiCells = this->psi_.internalField();
    scalarField& muCells = this->mu_.internalField();
    scalarField& alphaCells = this->alpha_.internalField();

    forAll(TCells, celli)
    {
        const typename MixtureType::thermoType& mixture_ =
            this->cellMixture(celli);

        TCells[celli] = mixture_.THE
        (
            hCells[celli],
            pCells[celli],
            TCells[celli]
        );

        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);

        muCells[celli] = mixture_.mu(pCells[celli], TCells[celli]);
        alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);

        TuCells[celli] = this->cellReactants(celli).THE
        (
            heuCells[celli],
            pCells[celli],
            TuCells[celli]
        );
    }

    forAll(this->T_.boundaryField(), patchi)
    {
        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
        fvPatchScalarField& pTu = this->Tu_.boundaryField()[patchi];
        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];

        fvPatchScalarField& ph = this->he_.boundaryField()[patchi];
        fvPatchScalarField& pheu = this->heu_.boundaryField()[patchi];

        fvPatchScalarField& pmu_ = this->mu_.boundaryField()[patchi];
        fvPatchScalarField& palpha_ = this->alpha_.boundaryField()[patchi];

        if (pT.fixesValue())
        {
            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);

                ph[facei] = mixture_.HE(pp[facei], pT[facei]);

                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                pmu_[facei] = mixture_.mu(pp[facei], pT[facei]);
                palpha_[facei] = mixture_.alphah(pp[facei], pT[facei]);
            }
        }
        else
        {
            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);

                pT[facei] = mixture_.THE(ph[facei], pp[facei], pT[facei]);

                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                pmu_[facei] = mixture_.mu(pp[facei], pT[facei]);
                palpha_[facei] = mixture_.alphah(pp[facei], pT[facei]);

                pTu[facei] =
                    this->patchFaceReactants(patchi, facei)
                   .THE(pheu[facei], pp[facei], pTu[facei]);
            }
        }
    }
}

Now, to calculate the temperature the function THE is called from mixture_ as mixture_.THE
After some search I think I have to go to "specie/thermo/thermo/thermoI.H". I'm not sure why. It would be great if you could tell me if I am right? and if yes why?.


So in
"specie/thermo/thermo/thermoI.H":

 template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::THE
(
    const scalar he,
    const scalar p,
    const scalar T0
) const
{
    return Type<thermo<Thermo, Type> >::THE(*this, he, p, T0);
}

Then it seems it goes to "specie/thermo/absoluteEnthalpy/absoluteEnthalpy.H". But I don’t know, where it is defined for the code that we are using "absolute enthalpy". With other word, where in the code the energy form will be read?

//- Temperature from absolute enthalpy
            //  given an initial temperature T0
            scalar THE
            (
                const Thermo& thermo,
                const scalar h,
                const scalar p,
                const scalar T0
            ) const
            {
                return thermo.THa(h, p, T0);
}

and again back to "specie/thermo/thermo/thermoI.H" to calculate THa

template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::THa
(
    const scalar ha,
    const scalar p,
    const scalar T0
) const
{
    return T
    (
        ha,
        p,
        T0,
        &thermo<Thermo, Type>::Ha,
        &thermo<Thermo, Type>::Cp,
        &thermo<Thermo, Type>::limit
    );
}

and from the same h-data the temperature would be calculated using Newton-Raphson method:

template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::T
(
    scalar f,
    scalar p,
    scalar T0,
    scalar (thermo<Thermo, Type>::*F)(const scalar, const scalar) const,
    scalar (thermo<Thermo, Type>::*dFdT)(const scalar, const scalar)
        const,
    scalar (thermo<Thermo, Type>::*limit)(const scalar) const
) const
{
    scalar Test = T0;
    scalar Tnew = T0;
    scalar Ttol = T0*tol_;
    int    iter = 0;

    do
    {
        Test = Tnew;
        Tnew =
            (this->*limit)
            (Test - ((this->*F)(p, Test) - f)/(this->*dFdT)(p, Test));

        if (iter++ > maxIter_)
        {
            FatalErrorIn
            (
                "thermo<Thermo, Type>::T(scalar f, scalar T0, "
                "scalar (thermo<Thermo, Type>::*F)"
                "(const scalar) const, "
                "scalar (thermo<Thermo, Type>::*dFdT)"
                "(const scalar) const, "
                "scalar (thermo<Thermo, Type>::*limit)"
                "(const scalar) const"
                ") const"
            )   << "Maximum number of iterations exceeded"
                << abort(FatalError);
        }

    } while (mag(Tnew - Test) > Ttol);

    return Tnew;
}

So, F is probably the absolute Enthalpy as a function of temperature and the constants of Nasa polynomial from "src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H"

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::ha
(
    const scalar p,
    const scalar T
) const
{
    const coeffArray& a = coeffs(T);
    return this->RR*
    (
        ((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T
      + a[5]
    );
}

and dFdT is its derivative with respect to T.

So my questions:

1: How are F and dFdT defined?
2: Do we need pressure (p) for absolute enthalpy?
3: Don't we need velocity to calculate the kinetic energy and substract it from the absoulute enthalpy to calculate temperature?
4:Is the Nasa polynomial (see http://combustion.berkeley.edu/gri_m...nasa_plnm.html) written for absolute or sensible enthalpy?

Any help will be appreciated.

[fluid2020]

Did you get the solution for this? Can anyone give a comment on this topic?