CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Temperature calculation from total enthalpy in OpenFOAM (XiFOAM)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By sharifi

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 25, 2017, 12:10
Default Temperature calculation from total enthalpy in OpenFOAM (XiFOAM)
  #1
New Member
 
Jim KIT
Join Date: Aug 2012
Location: Germany
Posts: 25
Rep Power: 14
sharifi is on a distinguished road
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:

Code:
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:

Code:
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:


Code:
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":

Code:

Code:
 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);
}
Code:



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?

Code:
//- 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);
}
Code:



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

Code:
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:

Code:
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"



Code:
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.

jennyjian likes this.

Last edited by sharifi; August 26, 2017 at 03:48.
sharifi is offline   Reply With Quote

Old   October 8, 2020, 09:16
Default
  #2
New Member
 
Join Date: Sep 2020
Posts: 11
Rep Power: 6
fluid2020 is on a distinguished road
Did you get the solution for this? Can anyone give a comment on this topic?
fluid2020 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
whats the cause of error? immortality OpenFOAM Running, Solving & CFD 13 March 24, 2021 07:15
Sensible Enthalpy vs. Total enthalpy Sethi OpenFOAM Programming & Development 3 March 23, 2017 21:48
How does fluent calculate the total enthalpy? bap FLUENT 0 June 22, 2016 10:55
Error in run Batch file saba1366 CFX 4 February 10, 2013 01:15
temperature / enthalpy fields depending on type of fvPatchField astein OpenFOAM Programming & Development 0 June 28, 2010 07:10


All times are GMT -4. The time now is 19:11.