The energy conservation equation is expressed in terms of internal energy \(e\) or enthalpy \(h\) in OpenFOAM and the temperature field is calculated from these solution variables. More precisely, we can specify the energy variable in the energy entry in the thermophysicalProperties file and the available options are the followings:
- specific sensible enthalpy \(h {\rm [J/kg]}\) (sensibleEnthalpy)
- specific sensible internal energy \(e {\rm [J/kg]}\) (sensibleInternalEnergy)
If we select the hePsiThermo thermophysical model, the temperature field \(T\) is calculated from the solved energy variable in the following function, where the compressibility \(\psi\), dynamic viscosity \(\mu\) and thermal diffusivity \(\alpha\) are also calculated.
Source Code
- hePsiThermo: hePsiThermo.C
- (heRhoThermo: heRhoThermo.C)
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 |
template<class BasicPsiThermo, class MixtureType> void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate() { const scalarField& hCells = this->he_; const scalarField& pCells = this->p_; scalarField& TCells = this->T_.primitiveFieldRef(); scalarField& psiCells = this->psi_.primitiveFieldRef(); scalarField& muCells = this->mu_.primitiveFieldRef(); scalarField& alphaCells = this->alpha_.primitiveFieldRef(); 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]); } volScalarField::Boundary& pBf = this->p_.boundaryFieldRef(); volScalarField::Boundary& TBf = this->T_.boundaryFieldRef(); volScalarField::Boundary& psiBf = this->psi_.boundaryFieldRef(); volScalarField::Boundary& heBf = this->he().boundaryFieldRef(); volScalarField::Boundary& muBf = this->mu_.boundaryFieldRef(); volScalarField::Boundary& alphaBf = this->alpha_.boundaryFieldRef(); forAll(this->T_.boundaryField(), patchi) { fvPatchScalarField& pp = pBf[patchi]; fvPatchScalarField& pT = TBf[patchi]; fvPatchScalarField& ppsi = psiBf[patchi]; fvPatchScalarField& phe = heBf[patchi]; fvPatchScalarField& pmu = muBf[patchi]; fvPatchScalarField& palpha = alphaBf[patchi]; if (pT.fixesValue()) { forAll(pT, facei) { const typename MixtureType::thermoType& mixture_ = this->patchFaceMixture(patchi, facei); phe[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(phe[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]); } } } } |
The calculation procedure of the temperature field depends on the selected energy variable, so the following THE function accordingly switches the called method.
364 365 366 367 368 369 370 371 372 373 |
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); } |
If we choose the sensible enthalpy as the energy variable, the following THs function is called to calculate the temperature from the sensible enthalpy.
117 118 119 120 121 122 123 124 125 126 127 128 |
//- Temperature from sensible enthalpy // given an initial temperature T0 scalar THE ( const Thermo& thermo, const scalar h, const scalar p, const scalar T0 ) const { return thermo.THs(h, p, T0); } |
376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 |
template<class Thermo, template<class> class Type> inline Foam::scalar Foam::species::thermo<Thermo, Type>::THs ( const scalar hs, const scalar p, const scalar T0 ) const { return T ( hs, p, T0, &thermo<Thermo, Type>::Hs, &thermo<Thermo, Type>::Cp, &thermo<Thermo, Type>::limit ); } |
The calculation of the temperature is done iteratively using the Newton-Raphson method.
376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 |
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 { if (T0 < 0) { FatalErrorInFunction << "Negative initial temperature T0: " << T0 << abort(FatalError); } 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_) { FatalErrorInFunction << "Maximum number of iterations exceeded: " << maxIter_ << abort(FatalError); } } while (mag(Tnew - Test) > Ttol); return Tnew; } |
If the specific heat capacity at constant pressure \(c_p\) is expressed in the form of temperature polynomial function (hPolynomial)
\begin{equation}
c_p(T) = \sum_{i=0}^7 c_i T^i, \tag{1}
\end{equation}
the temperature in the j-th cell \(T_j\) is calculated from the following equation
\begin{equation}
\displaystyle \int_{T_{std}}^{T_j} \left( \sum_{i=0}^7 c_i T^i \right) dT = h_j, \tag{2}
\end{equation}
where \(h_j\) is the sensible enthalpy value in the j-th cell. In general the equation will be nonlinear, the iterative solution technique is implemented.
Compressible Flow Solvers |
The above function calculate() is called by rhoPimpleFoam, rhoSimpleFoam and sonicFoam etc. from the line “thermo.correct()” after solving the energy conservation equation EEqn:
27 28 29 30 31 |
EEqn.solve(); fvOptions.correct(he); thermo.correct(); |
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 |
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class BasicPsiThermo, class MixtureType> void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct() { if (debug) { InfoInFunction << endl; } // force the saving of the old-time values this->psi_.oldTime(); calculate(); if (debug) { Info<< " Finished" << endl; } } |
After updating the compressibility field \(\psi\), the pressure Poisson equation pEqn is constructed and solved.
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 |
while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvc::div(phiHbyA) - fvm::laplacian(rhorAUf, p) == fvOptions(psi, p, rho.name()) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA + pEqn.flux(); } } |
The density field \(\rho\) is calculated with updated pressure and compressibility fields.
81 82 83 84 85 86 87 88 89 90 91 |
// Explicitly relax pressure for momentum corrector p.relax(); U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); pressureControl.limit(p); p.correctBoundaryConditions(); rho = thermo.rho(); |
93 94 95 96 |
Foam::tmp<Foam::volScalarField> Foam::psiThermo::rho() const { return p_*psi_; } |
98 99 100 101 102 |
template<class Specie> inline Foam::scalar Foam::perfectGas<Specie>::psi(scalar p, scalar T) const { return 1.0/(this->R()*T); } |
Summary |
Energy variable | Function used to calculate temperature |
sensible enthalpy | THs |
absolute enthalpy | THa |
sensible internal energy | TEs |
absolute internal energy | TEa |
Study with great curiosity and let’s see the physical phenomena from a HIGHER point of view!
Thanks a lot. This is a very useful content.
thanks for your useful information.
how to add temperature equation to viscoelasticFluidFoam solver.