# Temperature calculation from total enthalpy in OpenFOAM (XiFOAM)

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

 August 25, 2017, 13:10 Temperature calculation from total enthalpy in OpenFOAM (XiFOAM) #1 New Member   Jim KIT Join Date: Aug 2012 Location: Germany Posts: 25 Rep Power: 12 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 void Foam::heheuPsiThermo::correct() { if (debug) { Info<< "entering heheuPsiThermo::correct()" << endl; } // force the saving of the old-time values this->psi_.oldTime(); calculate(); if (debug) { Info<< "exiting heheuPsiThermo::correct()" << endl; } }``` The function calculate can be found in the same h-data: Code: ```template void Foam::heheuPsiThermo::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 Type> inline Foam::scalar Foam::species::thermo::THE ( const scalar he, const scalar p, const scalar T0 ) const { return 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 Type> inline Foam::scalar Foam::species::thermo::THa ( const scalar ha, const scalar p, const scalar T0 ) const { return T ( ha, p, T0, &thermo::Ha, &thermo::Cp, &thermo::limit ); }``` and from the same h-data the temperature would be calculated using Newton-Raphson method: Code: ```template class Type> inline Foam::scalar Foam::species::thermo::T ( scalar f, scalar p, scalar T0, scalar (thermo::*F)(const scalar, const scalar) const, scalar (thermo::*dFdT)(const scalar, const scalar) const, scalar (thermo::*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::T(scalar f, scalar T0, " "scalar (thermo::*F)" "(const scalar) const, " "scalar (thermo::*dFdT)" "(const scalar) const, " "scalar (thermo::*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 inline Foam::scalar Foam::janafThermo::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 04:48.

 October 8, 2020, 10:16 #2 New Member   Join Date: Sep 2020 Posts: 11 Rep Power: 4 Did you get the solution for this? Can anyone give a comment on this topic?