temperature calculation from enthalpy

stockzahn · December 16, 2019, 11:36

Dear Foamers,

I'm trying to implement tables and/or functions for the calculation of state properties in OF6. I adapt the rhoPimpleFoam solver and unfortunately I'm stuck trying to change the temperature calculation as function of enthalpy and pressure.

That's my assumption:

1) Using "heRhoThermo" the member function "calculate" in "heRhoThermo.C" calculates the temperature field from solved state variables (the part where the temperature is calculated is colored red).

Code:

template<class BasicPsiThermo, class MixtureType>
void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
{
    const scalarField& hCells = this->he();
    const scalarField& pCells = this->p_;

    scalarField& TCells = this->T_.primitiveFieldRef();
    scalarField& psiCells = this->psi_.primitiveFieldRef();
    scalarField& rhoCells = this->rho_.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]);
        rhoCells[celli] = mixture_.rho(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& rhoBf =
        this->rho_.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& prho = rhoBf[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]);
                prho[facei] = mixture_.rho(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]);
                prho[facei] = mixture_.rho(pp[facei], pT[facei]);
                pmu[facei] = mixture_.mu(pp[facei], pT[facei]);
                palpha[facei] = mixture_.alphah(pp[facei], pT[facei]);
            }
        }
    }
}

2) I've found several member functions called "THE", the most promising though leading to the file "absoluteEnthalpy.H" (I couldn't find one in the class "basicMixture" or in a derived class), which leads to a function called "THa":

Code:


            scalar THE
            (
                const Thermo& thermo,
                const scalar h,
                const scalar p,
                const scalar T0
            ) const
            {
                return thermo.THa(h, p, T0);
            }

3) The function "THa" seems to be located in the file "thermoI.H":

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

4) Which finally (?) returns a volScalarField "T".

But until now I couldn't find the piece of code calculating the temperature field from the given input variables.

Does anybody know, where to find that or am I on the wrong track?

Thanks in advance and have a nice evening (or whatever),
stockzahn

clapointe · December 16, 2019, 11:54

Check out this blog post : https://caefn.com/openfoam/temperature-calculation.

Caelan

stockzahn · December 19, 2019, 07:00

Quote:

Originally Posted by clapointe

Check out this blog post : https://caefn.com/openfoam/temperature-calculation.

Caelan

Dear Caelan,

thank you very much for this link - the information indeed was very helpful and (I think) I was able to follow almost all the steps how the enthalpy/internal energy field is created and, after solving the LES, the temperature is obtained from the resulting field.

However, unfortunately I have a missing link I couldn't find and maybe you or somebody else know where to find it. Examplarily I use the rhoPimpleFoam solver to explain my problem:

In the file createFields.H the object thermo of the class fluidThermo is created.

Code:

 autoPtr<fluidThermo> pThermo 
 (
 
     fluidThermo::New(mesh)
 
 );
 
 fluidThermo& thermo = pThermo();

fluidThermo inherits the methods of its parents basicThermo and compressibleTransportModel, whereas the method he() to calculate and retrieve the energy field is declared in the basicThermo class.

Code:

virtual tmp<volScalarField> he 
             (
 
                 const volScalarField& p,
 
                 const volScalarField& T
 
             ) const = 0;

The calculation itself is done by heThermo, which is a template class.

Code:

template<class BasicThermo, class MixtureType>
 
 class heThermo
 
 :
 
     public BasicThermo,
 
     public MixtureType
 
 {
  ...


             virtual tmp<volScalarField> he 
             (
 
                 const volScalarField& p,
 
                 const volScalarField& T
 
             ) const;     



  ...
 }

My problem is that I don't understand how the object thermo (class fluidThermo) "knows" to use the method defined in the template class heThermo, wenn calling the method .he() in EEqn.H.

Code:


volScalarField& he = thermo.he();

How/where a heThermo object is created or how does the object thermo knows to call the function from heThermo.C? Do you maybe know the answer to that?

stockzahn

fetc95 · April 20, 2020, 05:15

Hello,

I am concerned with the same problem, but the link provided above does not work. Could you explain how the temperature is obtained from the resulting field if you understand it?

The meaning of "he" is more intricated for me be because of diffusivity term in the EEq.H of buoyantPimpleFoam. In the simplest case the energy equation in the buoyantPimpleFoam can simplified as

Code:

fvm::ddt(rho, he) 
      + fvm::div(phi, he)
      - fvm::laplacian(alpha, he)
     ==
      + fvOptions(rho, he)

Let "he" is an energy which is calculated in some way. The term

Code:

fvm::laplacian(alpha, he)

is seemed to be a heat conductivity, but why the "he" is under Laplacian and not the temperature? The heat diffusivity is proportional to the gradient of the temperature, not energy.

Regards

stockzahn · April 20, 2020, 11:49

Dear fetc95,

I currently dealt with exactly the same question. Unfortunately I wasn't able to find an answer I'm sure about, so I only can tell you my thoughts. Maybe you can verify or falsify them.

- The diffusion term normally should be fvc::laplacian(kappa, T)
- The disadvantage is that the diffusion cannot be calculated implicitly
- The term fvm::laplacian(alpha, he) works if the enthalpy/internal energy is proportional to the temperature; only the factor Cp^-1/Cv^-1 (specific heat capacity), which is contained in alpha, must be taken into account (in OF alpha = kappa/Cp)
- This works for simple models, e.g. perfect gas
- In my case I integrated a property library (of helium); at some states helium has an inverse correlation between temperature and enthalpy/internal energy, therefore I think I must not use the implicit term, but have to use fvc::laplacian(kappa, T) which increase instability

To sum up: In my understanding the implicit term fvm::laplacian(alpha, he) can be used, when

1) T proportional to he
2) Cp * T = h | Cv * T = e

otherwise the code has to be changed to calculate the diffusion explicitly from fvc::laplacian(kappa, T)

What do you think?

fetc95 · April 20, 2020, 19:45

Quote:

Originally Posted by stockzahn

Dear fetc95,

I currently dealt with exactly the same question. Unfortunately I wasn't able to find an answer I'm sure about, so I only can tell you my thoughts. Maybe you can verify or falsify them.

- The diffusion term normally should be fvc::laplacian(kappa, T)
- The disadvantage is that the diffusion cannot be calculated implicitly
- The term fvm::laplacian(alpha, he) works if the enthalpy/internal energy is proportional to the temperature; only the factor Cp^-1/Cv^-1 (specific heat capacity), which is contained in alpha, must be taken into account (in OF alpha = kappa/Cp)
- This works for simple models, e.g. perfect gas
- In my case I integrated a property library (of helium); at some states helium has an inverse correlation between temperature and enthalpy/internal energy, therefore I think I must not use the implicit term, but have to use fvc::laplacian(kappa, T) which increase instability

To sum up: In my understanding the implicit term fvm::laplacian(alpha, he) can be used, when

1) T proportional to he
2) Cp * T = h | Cv * T = e

otherwise the code has to be changed to calculate the diffusion explicitly from fvc::laplacian(kappa, T)

What do you think?

Hello again,

I thought about this issue yesterday and today and that what I think.

It means that all compressible, combustion, heat transfer, etc. solvers are not valid for cases of temperature-dependent heat capacity and are valid only when hConst or eConst is choosen. I believe such important fact would be mentioned somewhere in openfoam if it was so. I also found the following in CFD direct https://cfd.direct/openfoam/energy-equation/. I did not find there any information about requirement of proportionality between

and

.

I suggest the following. Let us want to solve some heat equation with temperature-dependent heat capacity, e.g., hPolynomial. So in LHS of this equation we have a term
$c\left(T\right)\frac{\partial T}{\partial t}$
If we want to solve it in conservative form we can reformulate it by locating c\left(T\right) under derivative as follows
$c\left(T\right)\frac{\partial T}{\partial t} = \frac{\partial}{\partial t}\left( \int_{TRef}^T c\left(\theta \right) d\theta \right) = \frac{\partial he}{\partial t}$
Similarly we can do when formulating the heat flux law
$-\vec{q} = \lambda grad T$
by multiplying and dividing it by $c\left(T\right)$ and locating $c\left(T\right)$ under grad:
$-\vec{q} = \frac{\lambda}{c\left(T\right)} c\left(T\right) grad T = \frac{\lambda}{c\left(T\right)} grad \left( \int_{TRef}^T c\left(\theta \right) d\theta \right) = \alpha grad (he)$

Check me am I right mathematically.

So, I suggest that there is no requirement of proportionality between

and

, and the only requirement is the proportionality between

and

. Thus, the above formulation of the heat equation is valid for cases of temperature-dependent heat capacity. This is to be valid for any function c(T) which satisfies the following:

1) Finite number of points of jump discontinuities in range [TRef,TMax]
2) No essential discontinuities in range [TRef,TMax]

Another question is that how the temperature is calculated in case of non-constant heat capacity. E.g. if the c(T) is a polynomial of some order n
$c\left(T\right) = P_n\left(T\right)$
he is a polynomial of n+1 order
$he\left(\vec{x},t\right) = P_{n+1}\left(T(\vec{x},t)\right)$
The above algebraic polynomial equation is to be solved for temperature for all x, y, z and t. The advantage of this equation is that we know that for T>0 the he as function of T is always monotonical. It means that for fixed x, y, z and t it has only one solution. It seems that there is some iteration procedure is needed in OpenFOAM to obtain the temperature. It would be good to find it in the code.

stockzahn · April 21, 2020, 10:56

First of all thanks for this discussion. I also would like to start with the plausibility check:

Quote:

Originally Posted by fetc95

It means that all compressible, combustion, heat transfer, etc. solvers are not valid for cases of temperature-dependent heat capacity and are valid only when hConst or eConst is choosen. I believe such important fact would be mentioned somewhere in openfoam if it was so. I also found the following in CFD direct https://cfd.direct/openfoam/energy-equation/. I did not find there any information about requirement of proportionality between

and

.

That's a valid point, it would be strange, if they'd implement a term which only can be used in very rare and unrealistic cases (I have to check the combustion solver etc., I actually don't now if they use the same term).

On the other hand, the idea to treat the diffusion term implicitly to increase stability seems so obvious that I would be surprised, if I wouldn't have read it in a book until now - but actually I haven't. Maybe I've got the wrong books ...

Regarding the math, just to be sure about the last line: Locating the specific heat capacity inside the gradient only is allowed, if it is not a function of a spatial coordinate (so c != f(x,y,z)), so if it is constant. Or am I wrong? But for the discretized equation we can use (different) constant values for each cell, therefore I would agree using discrete values. Then the temperature gradient only would be multiplied by certain factor, which is "corrected" by using alpha instead of kappa. Considering it like that, it maybe could work anyway for discretized functions, independent of the model to calculate the specific property.

But for sure the code would be helpful, since normally in the discretized diffusion term the face values of the diffusive parameter are used, which would lead to a (small) error, if theay are corrected by the cell value of the diffusive parameter (alpha in this case). E.g. for 1D:

$D = \frac{\Delta t}{\Delta x^2 c_C} \left[\kappa_w c_w T_W - \left(\kappa_w c_w + \kappa_e c_e \right) T_C + \kappa_e c_e T_E \right]$

Still, somehow for me this only works if he is more or less an "image" of the temperature - multiplied by a factor c. I have to admit that until now I tried to find the code for "laplacian" a few times, but there are many different files and for me the structure of this basic files functions is not clear yet, so I never spent the time to dig into it until now - although I now that sooner or later that will be necessar anyway ...

What do you think? And thanks again for starting this discussion

stockzahn

fetc95 · April 21, 2020, 17:08

Hello,

I am also grateful that you participate this discussion. I firstly qoute this:

Quote:

Originally Posted by stockzahn

On the other hand, the idea to treat the diffusion term implicitly to increase stability seems so obvious that I would be surprised, if I wouldn't have read it in a book until now - but actually I haven't. Maybe I've got the wrong books ...
Regarding the math, just to be sure about the last line: Locating the specific heat capacity inside the gradient only is allowed, if it is not a function of a spatial coordinate (so c != f(x,y,z)), so if it is constant. Or am I wrong? But for the discretized equation we can use (different) constant values for each cell, therefore I would agree using discrete values. Then the temperature gradient only would be multiplied by certain factor, which is "corrected" by using alpha instead of kappa. Considering it like that, it maybe could work anyway for discretized functions, independent of the model to calculate the specific property.

stockzahn

I offer to remember the properties of complicated functions and their derivatives. Concretely, we have an enthalpy or internal energy as function of temperature which is function of x,y,z, i.e. "he" depends on coordinates only by means of temperature. Consider the derivative of "he"

$\frac{\partial he\left(T\left(x,y,z\right)\right)}{\partial x} = \frac{\partial he\left(T\left(x,y,z\right)\right)}{\partial T\left(x,y,z\right)} \frac{\partial T\left(x,y,z\right)}{\partial x}$

The specific heat capacity in our case is temperature derivative of "he". For example see the paper 10.1016/j.applthermaleng.2009.04.018 figures 1 and 2. The blue line is internal energy, the black line is specific heat capacity. You can see that the black line is derivative of the blue line.
So, on the contrast, you can back locate this term under derivative. Note that when I located the specific heat capacity under grad I located not directly heat capacity but its integral over temperature. In the simplest case of constant specific heat capacity this integral equals just c*T.

I think we must to find code of obtaining temperature from "he". If that code is constructed for c(T), not only constant, it is obviously that this heat equation is valid for any function c(T).

But while we will look for the code I offer to turn back to your concrete case to which the current thread is devoted. Did I understand you correctly that when you implemented fvm::laplacian(alpha, he) you obtained the unexpected results and when you implemented fvc::laplacian(kappa,T) you obtained the results you are expected but the stability decreased?

stockzahn · April 22, 2020, 13:57

Quote:

Originally Posted by fetc95

I offer to remember the properties of complicated functions and their derivatives. Concretely, we have an enthalpy or internal energy as function of temperature which is function of x,y,z, i.e. "he" depends on coordinates only by means of temperature.

Yes, of course. My bad. I'm thinking too much by means of discretized functions.

Quote:

Originally Posted by fetc95

Note that when I located the specific heat capacity under grad I located not directly heat capacity but its integral over temperature. In the simplest case of constant specific heat capacity this integral equals just c*T.

Analytically you sum up all the infinitely small products of dt*c, sure. For the numerical treatment in OF we would need different laplacian algorithms, depending on the model for c. Assuming e.g. a polynomial approach: There is no entry for a reference temperature - by the way in thermoI.H the function T calculates the temperature from "he" and p. If the calculation of c for the entire temperature range starting from absolute zero has to be done every iteration step (which would be computationally costly), one polynom must cover the entire temperatire range of the property. There is the option to define limit temperatures (limit), but in the polynomial approach (hPolynomialThermoI.H) this limit cannot be used as reference temperatrue, but is used by the T() in thermoI.H.

Quote:

Originally Posted by fetc95

I think we must to find code of obtaining temperature from "he". If that code is constructed for c(T), not only constant, it is obviously that this heat equation is valid for any function c(T).

So your (analytical) consideration, at least as far as I understand, is not implemented in the thermo models. I also can't image that there are different laplacian algorithms depending on the thermo model, but of course it is possible. Using the normal numerical laplacian operator for finite volumes (as in my last post), I cannot imagine how the code should now when using the face values and when the cell values in the discretized function, which would result in a small mistake since the face values of the diffusion coefficient differ from the cell value of the calculated cell. But as you suggested: We need to find the code used for solving the laplacian termin the energy equation - looking at the structure and the information given in fvSchemes I doubt there are different algorithms apart from the typical schemes.

I'm not sure if I will be able to start looking for the laplacian code before the weekend, but maybe I will try it then.

Quote:

Originally Posted by fetc95

But while we will look for the code I offer to turn back to your concrete case to which the current thread is devoted. Did I understand you correctly that when you implemented fvm::laplacian(alpha, he) you obtained the unexpected results and when you implemented fvc::laplacian(kappa,T) you obtained the results you are expected but the stability decreased?

Unexpected results, I don't want to give a concrete statement right now. But, the version using the explicit discretizartion with fvc::laplacian(kappa, T) is less stable.

stockzahn

fetc95 · April 24, 2020, 19:35

Dear stockzahn

Could you hint me what procedure do you use in your solver to get the temperature for using it in the equation?

stockzahn · April 26, 2020, 09:22

Quote:

Originally Posted by fetc95

Dear stockzahn

Could you hint me what procedure do you use in your solver to get the temperature for using it in the equation?

Dear fetc95,

when I made my first thermodynamical model, I tried to write down some procedures and explanations - unfortunately the document still contains some typos. However, in the end there are some graphics showing the calling paths for getting the properties. Does that answer your question? Otherwise please let me know which temperature at which point of the calculation you are interested in. Maybe I have an idea...

thermophysicalModel_small.pdf

sincerely,
stockzahn

fetc95 · April 26, 2020, 16:10

Quote:

Originally Posted by stockzahn

Dear fetc95,

when I made my first thermodynamical model, I tried to write down some procedures and explanations - unfortunately the document still contains some typos. However, in the end there are some graphics showing the calling paths for getting the properties. Does that answer your question? Otherwise please let me know which temperature at which point of the calculation you are interested in. Maybe I have an idea...

Attachment 76963

sincerely,
stockzahn

I was referring to which line of code you use to call the temperature in the solver. Is it thermo.T() or something else?

About your ideas which you tell earlier. I look at the of function T in thermoI.H. As i tell it is iteration procedure for obtaining temperature from some function F(T, p) which depends on case. Concretely it is the Newton method, which usually has good convergence rate so the temperature calculation has not be expencive.

When we consider the term laplacian(kappa/c(T), he(T)) we do not need any special implementation of the laplacian(*,*) for this case because if we consider this term on some time/iteration step (n+1), we consider the first argument of the function on step n. So from the point of view of (n+1) th step it is just a known volume field. All we need is to interpolate it on the cell faces.

stockzahn · April 27, 2020, 16:31

Quote:

Originally Posted by fetc95

I was referring to which line of code you use to call the temperature in the solver. Is it thermo.T() or something else?

yes, thermo.T() calls the temperature field which is stored in the thermo class.

Quote:

Originally Posted by fetc95

About your ideas which you tell earlier. I look at the of function T in thermoI.H. As i tell it is iteration procedure for obtaining temperature from some function F(T, p) which depends on case. Concretely it is the Newton method, which usually has good convergence rate so the temperature calculation has not be expencive.

I don't even use this function, since I retrieve my temperature values from a library function (T = f(p,h)). But just if I understand the code correctly: The function T in thermoI.H varies the temperature until it results in the correct enthalpy value (for a fixed pressure)? If I try to "translate" this line:

Code:

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

I get (ignoring the limit)

$T_{new} = T_{old} - \frac{h\left(p, T_{old}\right) - h_{known}}{c_p\left(p, T_{old}\right)}$

That again, if I understand correctly, can only work if $h = c_p \cdot T$ , or what am I missing? Unless c_p is just used as a mathematical factor without a physical meaning ...

Quote:

Originally Posted by fetc95

When we consider the term laplacian(kappa/c(T), he(T)) we do not need any special implementation of the laplacian(*,*) for this case because if we consider this term on some time/iteration step (n+1), we consider the first argument of the function on step n. So from the point of view of (n+1) th step it is just a known volume field. All we need is to interpolate it on the cell faces.

Still we use the enthalpy gradient instead of the temperature gradient. Can we use it, since in the small temperature range of a cell and its neighbours we can assume the specific heat capacity more or less constant? We do not need the absolute temperature, just the difference, therefore it doesn't really matter what the value of c is. But we'll get an error which increases with increasing variation of c with changing T.

In your previous posts you pointed out, that the approach would work if the enthalpy would be calculated as integral of c(T)dT. But I can't find any function doing that.

fetc95 · April 28, 2020, 21:42

Hello!

Quote:

Originally Posted by stockzahn

$T_{new} = T_{old} - \frac{h\left(p, T_{old}\right) - h_{known}}{c_p\left(p, T_{old}\right)}$

That again, if I understand correctly, can only work if $h = c_p \cdot T$ , or what am I missing? Unless c_p is just used as a mathematical factor without a physical meaning ...

I believe that if it would work only when h = cp*T, there would no need to realize that iteration procedure.

Let's try to gather information from your first post with your previous posts. As you rote the temperature is calculated

Code:

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]);
        rhoCells[celli] = mixture_.rho(pCells[celli], TCells[celli]);

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

The function THE is seemed to be called from the file which depends on chosen enthalpy model in the thermoPhysicalProperties dictionary of the considered case. This may be sinsibleEnthalpy.H, sensibleInternalEnthalpy.H, absoluteEnthalpy.H, etc. However in all this files function THE is implemented as

Code:

scalar THE
            (
                const Thermo& thermo,
                const scalar h,
                const scalar p,
                const scalar T0
            ) const
            {
                return thermo."FUNCTION"(h, p, T0);
            }

Where THa,THs, TEs or others are substituted instead FUNCTION depending on which of sensible, absolute, internal Sensible or some other have been chosen. Anyway this function are all implemented in the thermol.H in similar way. For example 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);
 
 }

This function calls the other function T which is also implemented in thermol.H. You wrote its code above. As it is written in the description of this function it returns the temperature corresponding to the value of the thermodynamic property f, given the function f = F(p, T) and dF(p, T)/dT. So this function solves the equation:

$F\left(T, p\right) - f = 0$

using the Newton method:

$T^{k+1} = T^k - \frac{F\left(T^k, p\right) - f}{\partial F\left(T^k, p\right)/\partial T}$

If we gather all of this we see that to calculate temperature we go through all cell and for every cell we call the function:

$T_{cell}^{presentTime} = T\left( h_{cell}^{presentTime}, p, T_{cell}^{previousTime}, h(T, p), c_p(T, p)\right)$

Since in the 5th argument of the function T() is as said in the description a derivative of the 4th argument, we see that cp is passed as derivative of h. And it does not depend on that cp is constant or not. We already due to this can believe that

$c_p = \frac{\partial h}{\partial T}$

and correspondingly

$h = \int_0^T {c_p dT}$

Now let's try to talk about integration of cp. We do not find an implementation of this integration in general case because it would be very costly numerically. So the integration is to depend on concrete chosen model hConst, hPolynimial, etc. For example when cp is constant this integral is

$h = \int_0^T {c_p dT} = c_p T$

and you can find in hConstThermol.H in line 109 that h is calculated as cp times T.

Now let's consider the case of hPolynomial when cp is polynomial. In this case is

$c_p = \sum\limits_{k=0}^N a_k T^k$

The integration of this gives

$h = b + \sum\limits_{k=0}^{N} \frac{a_k}{k+1} T^{k+1}$

Since we know that enthalpy or internal energy equals zero at zero temperature we set b=0. So, summarizing we can say that in implementation to set the polynomial we just need to store the order of the polynomial and array of its coefficients. Suppose we have an array V of size M=N+1 of polynomial coefficients, where

Such information we obtain from thermoPhysicalProperties when setting the cp-polynomial coefficients. So now all we need to set the integral of cp is set an array newV of size M+1 where

Code:

newV[0] = 0;
newV[j+1] = V[j]/(j+1), where j = 0,1, …, M

In hPolynimialThermo.H hPolynimialThermol.H hPolynimialThermo.C you can find that hCoeffs and cpCoeffs are declared as special type polyType which is implemented in polynomial.H and polynomial.C. Herewith the hCoeffs are calculated as

Code:

hCoeffs  = cpCpoeffs.integral().

where method integral() is implemented in polynomial.C and when applied for some array V of size N returns the array newV of size N+1 where newV[0] = 0, and other newV[j] are calculated as newV[j+1] = V[j]/(j+1).

We could check my suggestion by solving the following model problem. Consider the equation:

$\rho c_p\frac{\partial T}{\partial t} - \frac{\partial}{\partial x}\left( \kappa \frac{\partial T}{\partial x} \right) = S$

This equation can be reformulated as

$\rho \frac{\partial h}{\partial t} - \frac{\partial}{\partial x}\left( \alpha \frac{\partial h}{\partial x} \right) = S$

Setting the boundary condition this equation can be easily solved with finite-different method which is easy to implement. From the other hand this equation can be solved with thermoFoam. If my suggestion is right these solutions must concise.

stockzahn · May 4, 2020, 17:02

First of all: Sorry for the late reply - it was a strange week. Plus I obviously need much more time to process this mathematically. You are really good at this - I'm impressed. Thanks to your last post I think I've achieved a better understanding and thanks to your formulation I finally get how the calculation of

works.

If I get it right, for the polynomial approach, the calculation of the polynomial coefficents for the enthalpy is implemented based on the polynoms of the specific heat capacity. I assume similar pieces of code exist for the other models. Let's have a look at the calling path to the function T in thermoI.H (using the Newton 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
 {
     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;
 }

In the enthalpy case it is called by the function THa (or THs) also defined in thermoI.H:

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

The fifth argument is the specific heat capacity, whose function is called (thermo<Thermo, Type>::*dFdT)(const scalar, const scalar), defined in the hPolynomialThermoI.H file:

Code:

 template<class EquationOfState, int PolySize>
 inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::Cp
 (
     const scalar p, const scalar T
 ) const
 {
     return CpCoeffs_.value(T) + EquationOfState::Cp(p, T);
 }

Which in my interpretation returns a specific heat capacity value calculated by the polynomial model plus the summand "EquationOfState::Cp(p, T)", of which I'm not sure where to find it. I should test it, but my assumption for now is that it is defined in rPolynomialI.H:

Code:

 template<class Specie>
 inline Foam::scalar Foam::rPolynomial<Specie>::Cp(scalar p, scalar T) const
 {
     return 0;
 }

So I think I was able to understand the calculation of the temperature given an enthalpy (and pressure), but I still cannot see, how the discretized enthalpy diffusion could work

Code:

 fvm::laplacian(turbulence->alphaEff(), he)

. When the he-field is generated, at the end the function Ha (or Hs) in hPolynomialThermoI.H is called:

Code:

 template<class EquationOfState, int PolySize>
 inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::Ha
 (
     const scalar p, const scalar T
 ) const
 {
     return hCoeffs_.value(T) + EquationOfState::H(p, T);
 }

As you explained very clearly, the

Code:

 hCoeffs_

can be calculated from the

Code:

 cpCoeffs_

, but in my understanding it still only holds a current field value. It does not contain the polynomial coefficients and I can't find a function to retrieve them. There is a similar issue with your suggestion to check with FDM: From an analytical point of view that all makes sense, but the discretized form (using just values for the diffusion coefficients kappa and alpha) seems to have the same problem like I thought the OF method should have.

For the 1D-case the (explicit) FDM discretization of your last two equations (for

) should be

$\rho c_p\frac{\partial T}{\partial t} - \frac{\partial}{\partial x}\left( \kappa \frac{\partial T}{\partial x} \right) = S$
$T^{n+1} = \frac{\Delta t}{\rho} \underbrace{\frac{\kappa}{c_p}}_{\alpha} \frac{T^n_{i+1} - 2 T^n_{i} + T^n_{i-1}}{\Delta x ^2} + T^n_i$

$\rho \frac{\partial h}{\partial t} - \frac{\partial}{\partial x}\left( \alpha \frac{\partial h}{\partial x} \right) = S$
$h^{n+1} = \frac{\Delta t}{\rho} \alpha \frac{h^n_{i+1} - 2 h^n_{i} + h^n_{i-1}}{\Delta x ^2} + h^n_i$

Using just "values" of the diffusion coefficients, the equations for temperature and enthalpy are identical they differ only by the source terms at the boundaries. That works if we assume the difference of the specific heat capacities for two neighboured cells negligible. Or did you have in mind something else?

fetc95 · May 4, 2020, 19:41

Hi, stockzahn

I start the reply from the end.

Quote:

Originally Posted by stockzahn

There is a similar issue with your suggestion to check with FDM: From an analytical point of view that all makes sense, but the discretized form (using just values for the diffusion coefficients kappa and alpha) seems to have the same problem like I thought the OF method should have.

For the 1D-case the (explicit) FDM discretization of your last two equations (for

) should be

$\rho c_p\frac{\partial T}{\partial t} - \frac{\partial}{\partial x}\left( \kappa \frac{\partial T}{\partial x} \right) = S$
$T^{n+1} = \frac{\Delta t}{\rho} \underbrace{\frac{\kappa}{c_p}}_{\alpha} \frac{T^n_{i+1} - 2 T^n_{i} + T^n_{i-1}}{\Delta x ^2} + T^n_i$

$\rho \frac{\partial h}{\partial t} - \frac{\partial}{\partial x}\left( \alpha \frac{\partial h}{\partial x} \right) = S$
$h^{n+1} = \frac{\Delta t}{\rho} \alpha \frac{h^n_{i+1} - 2 h^n_{i} + h^n_{i-1}}{\Delta x ^2} + h^n_i$

Using just "values" of the diffusion coefficients, the equations for temperature and enthalpy are identical they differ only by the source terms at the boundaries. That works if we assume the difference of the specific heat capacities for two neighboured cells negligible. Or did you have in mind something else?

If I correctly understand your question is how to obtain the boundary conditions for h when solving the equation for it. These conditions can be easily calculated using the conditions for T.

I also cannot understand why you retrieve cp from derivative in equation for h.

I will consider the case of temperature dependent cp. Let's for example

The h in this case is

$h = \int_0^T{c_p(g)dg} = \int_0^T{2TdT} = T^2$

So the equation is (let $\rho = 1$ )

$\frac{\partial h}{\partial t} - \frac{\partial}{\partial x}\left( \frac{\kappa}{2T^2} \frac{\partial h}{\partial x} \right) = S$

The boundary conditions for this equation can be calculated using the conditions for temperature and the dependencies between h, cp and T. Let in the following

$x \in \left[ 0,1 \right], \left.T\right|_{x=0} = T_0, \left.\frac{\partial T}{\partial x}\right|_{x=1} = 0$

The condition for the function h at x=0 you can directly calculate from known function h = h(T):

$\left.h\right|_{x=0} = h\left( \left.T\right|_{x=0} \right) = h\left( T_0 \right) = T_0^2,$

The gradient condition at x = 1 can be obtained by multiplying and dividing it by cp and doing the same procedure as I used several posts above:

$\left.\frac{c_p}{c_p}\frac{\partial T}{\partial x}\right|_{x=1} = \left.\frac{1}{c_p}\frac{\partial}{\partial x} \left( \int_0^T{c_p(g)dg} \right) \right|_{x=1} = \left.\frac{1}{c_p}\frac{\partial h}{\partial x} \right|_{x=1}$

Finally

$\left.\frac{\partial h}{\partial x}\right|_{x=1} = 0$

In OF or in implicit FVM we will actually solve the following equation on every moment of time:

$- \Delta t \frac{\partial}{\partial x}\left( f\left( x \right) \frac{\partial h^{n+1}}{\partial x} \right) + h^{n+1} = \Delta t S^{n+1} + h^n$

with conditions:
$\left.h\right|_{x=0} = T_0^2, \left.\frac{\partial h}{\partial x}\right|_{x=1} = 0$

where

$f\left( x \right) = \frac{\kappa}{c_p \left( T^n \right)}$

f(x) is calculated from temperature obtained on the previous time step and on the current step represents just a known field. It seems some decoupling, but it can be compensated by doing several iterations on every time step as it is in OF. And there is know need to assume the difference of the specific heat capacities for two neighboured cells negligible. You just need to solve the Poison equation on every step with coordinate-dependent diffusion coefficient.

Quote:

but in my understanding it still only holds a current field value. It does not contain the polynomial coefficients and I can't find a function to retrieve them.

Honestly, I cannot understand what you mean. cpCoeffs_ is an object which has different methods. It stores the coefficients and the order of the polynomial and returns the value of field when is called as

Code:

cpCoeffs_.value(T)

This routine returns the value at point T. But you are right when say that in solvers cpCoeffs returns the "current field" since it is calculated from the "current" field T which means the field obtained on the previous time/iteration step.

Quote:

When the he-field is generated, at the end the function Ha (or Hs) in hPolynomialThermoI.H is called:

The he-field is not generated from Ha(). He-field is obtained by you when you solve the equation in the solver. As I understand the function Ha() is used only to say openfoam what is the relationship between enthalpy and temperature which it has to use when obtaining the temperature from the known value of the enthalpy by Newton's method.

stockzahn · May 5, 2020, 09:39

Quote:

Originally Posted by fetc95

If I correctly understand your question is how to obtain the boundary conditions for h when solving the equation for it. These conditions can be easily calculated using the conditions for T.

I also cannot understand why you retrieve cp from derivative in equation for h.

I will consider the case of temperature dependent cp. Let's for example

No, sorry, that's a misunderstanding - I'm aware of how to apply the boundary condition. I try to formulate my thoughts more clearly: How is it done in OF? That's the line from the energy equation:

Code:

 fvm::laplacian(turbulence->alphaEff(), he)

he (e.g. in rhoPimpleFoam) is a reference to the energy field which is called in each SIMPLE iteration step when including the file EEqn.H:

Code:

 volScalarField& he = thermo.he();

It calls the energy field he_, which is a member of the heThermo class (heThermo.H).

Code:

virtual volScalarField& he()
{
    return he_;
}

The boundary field of he_ seems to be recalculated (h[celli] = f(p[celli],T[celli])) according to the calling path illustrated in figure 1 of the pdf from a previous post (e.g. when he is used to set up LES or when calculating T = f(p, h), maybe by calling the member function boundaryFieldRef() of a GeometricField). The function Ha (or Hs) then returns a value for the enthalpy - nothing else, as well as values are stored in he_ which can be seen since he_ is of the class volScalarField.

turbulence->alphaEff() is a function in the turbulence class, which calls the function alphaEff(alphat) in heThermo.C:

Code:

template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::volScalarField>
Foam::heThermo<BasicThermo, MixtureType>::alphaEff
(
    const volScalarField& alphat
) const
{
    tmp<Foam::volScalarField> alphaEff(this->CpByCpv()*(this->alpha_ + alphat));
    alphaEff.ref().rename("alphaEff");
    return alphaEff;
}

Also here, as far as I can see, a field is generated consisting only of values without any other information (class volScalarField). So if there is no special treatment the procedure looks like fvm::laplacian(volScalarField diffusionCoefficent, volScalarField diffusedProperty).

So I cannot see, where a function is used, for me it seems that (constant) values are passed to the solver. So I'm not saying it is not possible to set up a FDM simulation with cp = f(T) - of course we can. I don't get how OF does it and as I just tried to explain: What I can see from the code an according small FDM simulation would be run with a temperature dependent specific heat capacity, but not passing the function, just the value:

$c_{p,num} = c_p(p,T)$
$T^{n+1} = \frac{\Delta t}{\rho} \underbrace{\frac{\kappa}{c_{p,num}}}_{\alpha} \frac{T^n_{i+1} - 2 T^n_{i} + T^n_{i-1}}{\Delta x ^2} + T^n_i$

$\alpha_{num} = \alpha(p,T)$
$h^{n+1} = \frac{\Delta t}{\rho} \alpha_{num} \frac{h^n_{i+1} - 2 h^n_{i} + h^n_{i-1}}{\Delta x ^2} + h^n_i$

But maybe I have to dig a little bit deeper and calling paths contain more steps than I think right now.

fetc95 · May 5, 2020, 19:05

Unfortunately, I cannot understand where you see that the properties passed to the solvers must be constant. In what part of the code do you see that the field which is passed as diffusion property (the 1th argument of the laplacian) cannot be varied from cell to cell? The alpha (or alphaEff, it is not important) is a volScalarField as you say. But you also say "with no any information". The type vol<Type>Field already means that this property can be varied from cell to cell. On the current solution step this field is already known and can be used in matrix coefficients of the linear equation system which is solved for defining the sought field. What also information is needed? May be it is needed to find how the diffusion coefficient alpha is calculated.

In OpenFOAM descriptions, e.g., in Extended Code Guide in openfoam.com or in Programmers Guide from openfoam.org the term

Code:

laplacian(Gamma, G)

means

$\nabla \left( \Gamma \nabla G \right)$

or more exactly

$\nabla \left( \Gamma \left( x,y,z \right) \nabla G \right)$

so the diffusion coefficient in general can be the function of coordinates, i.e., volField or surfaceField and the passed properties must not be constant. In addition, I am sure that there is no point for the developers to implement the diffusion term only for constant diffusion coefficient. Anyway, I believe that if it was so it would be mentioned in the programmers guide. But if you are doubt we need to analyze the code of fvm::laplacian operator in similarity with that we done for thermo models. However, I think that it is obviously. Concretely, I am sure that after the Gauss integration over the finite volume of the above term and discretization of the surface gradients in 1D case the OF gets:

$\frac{1}{\Delta x^2} \left( \Gamma_{j+\frac{1}{2}}G_{j+1} - \left( \Gamma_{j+\frac{1}{2}} + \Gamma_{j-\frac{1}{2}} \right) G_j + \Gamma_{j-\frac{1}{2}}G_{j-1} \right)$

and not

$\frac{\Gamma_j}{\Delta x^2} \left(G_{j+1} - 2 G_j + G_{j-1} \right)$

From the information provided by you I cannot see that there is some contradiction and that it is impossible. The implementation of the first way of discretization is not programmatically expensive, so I do not think that OF developers ignored that.

Your second scheme is to be like this

$\alpha^{num}_i = \alpha(p^n_i,T^n_i)$
$h^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} h^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) h^n_{i} + \alpha^{num}_{i-1/2}h^n_{i-1}}{\Delta x ^2} + h^n_i$

where $\alpha^{num}_{i\pm1/2}$ are values of field $\alpha^{num}$ interpolated onto cell faces. Here is no any need to know how the $\alpha^{num}$ was calculated. All we need is to know the cell values of the field $\alpha^{num}$ which were obtained when thermo.correct() was called on the previous step. The Laplacian(alpha,h) just will take this known field and calculate the diffusion term.

stockzahn · May 6, 2020, 12:17

[QUOTE=fetc95;768930]
Unfortunately, I cannot understand where you see that the properties passed to the solvers must be constant.
[QUOTE/]

I don't know why I can't express my self adequately. I didn't mean that the entire field has the same value, but the single values of the cells (which of course can differ) are passed as numerical value not as functions. But I see that you can think that, since I wrote the equation in my last post as I did (and not like in my (overall) post #7). I try to be more precise.

I try again: If values are passed and not functions then

$T^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} T^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) T^n_{i} + \alpha^{num}_{i-1/2}T^n_{i-1}}{\Delta x ^2} + T^n_i$

corresponds to

$h^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} h^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) h^n_{i} + \alpha^{num}_{i-1/2}h^n_{i-1}}{\Delta x ^2} + h^n_i$

which only works, if $T \propto h$ or $\Delta T \propto \Delta h$ and that's only valid, if the change of

over the temperature range $T_{i-1}$ - $T_{i+1}$ including $T^{n+1}$ is negligible. Else the enthalpy (evolution) does not represent the temperature (evolution).

fetc95 · May 6, 2020, 20:06

[QUOTE=stockzahn;769086]

Quote:

Originally Posted by fetc95

Unfortunately, I cannot understand where you see that the properties passed to the solvers must be constant.
[QUOTE/]

I don't know why I can't express my self adequately. I didn't mean that the entire field has the same value, but the single values of the cells (which of course can differ) are passed as numerical value not as functions. But I see that you can think that, since I wrote the equation in my last post as I did (and not like in my (overall) post #7). I try to be more precise.

I try again: If values are passed and not functions then

$T^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} T^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) T^n_{i} + \alpha^{num}_{i-1/2}T^n_{i-1}}{\Delta x ^2} + T^n_i$

corresponds to

$h^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} h^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) h^n_{i} + \alpha^{num}_{i-1/2}h^n_{i-1}}{\Delta x ^2} + h^n_i$

which only works, if $T \propto h$ or $\Delta T \propto \Delta h$ and that's only valid, if the change of

over the temperature range $T_{i-1}$ - $T_{i+1}$ including $T^{n+1}$ is negligible. Else the enthalpy (evolution) does not represent the temperature (evolution).

Dear stockzahn,

It seems that you ignore some moments we already discussed above. The first scheme in your last post is valid only when cp = const, i.e., it is equivalent to (for non-constant kappa)

$T^{n+1} = \frac{\Delta t}{\rho c_p} \frac{\kappa^{num}_{i+1/2} T^n_{i+1} - \left( \kappa^{num}_{i+1/2} + \kappa^{num}_{i-1/2} \right) T^n_{i} + \kappa^{num}_{i-1/2}T^n_{i-1}}{\Delta x ^2} + T^n_i$

while the second scheme is valid for any variable cp. So, I cannot understand what you mean when you say that your two schemes "correspond" to each other. The first scheme does not correspond to the second one. You are right that your first scheme in your last post works only when cp=const. But it is not right for the second scheme, because it works for any cp.

Let us try again. You have an equation (rho = 1)

$c_p \frac{\partial T}{\partial t} - \frac{\partial}{\partial x}\left( \kappa \frac{\partial T}{\partial x} \right) = 0$

Let us in the following believe that cp and kappa depend on T

$\kappa = \kappa(T)$

This equation you also can write as follows

$c_p \frac{\partial T}{\partial t} - \frac{\partial}{\partial x}\left( \frac{c_p\kappa}{c_p} \frac{\partial T}{\partial x} \right) = 0$

Now we use that the changing of enthalpy is proportional to the changing of the temperature

which gives

$\frac{dh}{dT} = c_p(T)$

which gives

$h = \int_0^T c_p(\theta)d\theta$

For constant cp:

$h = \int_0^T c_pd\theta = c_pT$

For polynomial cp:

$h = \int_0^T \left( \sum\limits_{k=0}^N a_k \theta^k \right)d\theta = \sum\limits_{k=0}^{N} \frac{a_k}{k+1} T^{k+1}$

and for any function cp(T) which can be analytically integrated over T you can obtain the analytical expression for h=h(T).

As I wrote above you can use the following transformations

$c_p\frac{\partial T}{\partial t} = \frac{\partial}{\partial t}\left( \int_0^T c_p\left(\theta \right) d\theta \right) = \frac{\partial h}{\partial t}$

$\frac{\partial T}{\partial x} = \frac{c_p}{c_p} \frac{\partial T}{\partial x} =\frac{1}{c_p} \frac{\partial}{\partial x}\left( \int_0^T c_p\left(\theta \right) d\theta \right) = \frac{1}{c_p}\frac{\partial h}{\partial x}$

Substituting these in the original equation we get the fully mathematically equivalent equation

$\frac{\partial h}{\partial t} - \frac{\partial}{\partial x}\left( \alpha \frac{\partial h}{\partial x} \right) = 0$

where

$\alpha = \frac{\kappa}{c_p} = \frac{\kappa\left(T\right)}{c_p\left(T\right)}=\alpha\left(T\right)$

The only difference between this equation and the first equation is that the last is written in "conservative" form while the first one does not. The first equation is conservative only when cp=const. In other issues these equations are equivalent.

The procedure of the numerical solving of the first equation (in temperature terms) is as follows:
Suppose you calculated the temperature on some step n which is T^n. Now in the end of the nth step you calculate:

$c_{pj} = c_p(T^n_j)$

$\kappa_j = \kappa(T^n_j)$

After this you can pass them to the (n+1) step and calculate the temperature T^(n+1)

$T^{n+1}_j = \frac{\Delta t}{c_{pj}} \frac{\kappa_{j+1/2} T^n_{j+1} - \left( \kappa_{j+1/2} + \kappa_{j-1/2} \right) T^n_{j} + \kappa_{j-1/2}T^n_{j-1}}{\Delta x ^2} + T^n_j$

Note that our numerical method does not know and does not remember that cp_j and kappa_j were calculated as functions of temperature on the previous step. On the step (n+1) it senses them just as some known grid fields in FDM or volume fields in FVM. In other words the last scheme will not be changed if you will calculate cp_j and kappa_j on the previous step in some other way. On what step of this procedure do you need to pass them as "functions"? You just calculate this volume fields on previous step in some given way (in the present case - as functions of temperature, but it is not essential) and then pass them to the following step.

Similarly it is when you solve the equivalent equation for the enthalpy h. You calculated the temperature on some step n which is T^n. Now in the end of the nth step you calculate:

$\alpha_j = \alpha(T^n_j)$

After this you can pass alpha to the (n+1) step and calculate the enthalpy h^(n+1)

$h^{n+1}_j = \Delta t \frac{\alpha_{j+1/2} h^n_{j+1} - \left( \alpha_{j+1/2} + \alpha_{j-1/2} \right) h^n_{j} + \alpha_{j-1/2}h^n_{j-1}}{\Delta x ^2} + h^n_j$

After this you calculate the temperature T^(n+1) solving the equation:

$h\left(T^{n+1}_j\right) - h^{n+1}_j = 0$

for every j using the Newton method as we discussed above.

Similarly, we see that our method does not know and does not remember that alpha_j was calculated as function of temperature. And it does not need to know that. On the step (n+1) it senses alpha just as some known volume field. You just calculate this volume field on previous step in some given way and then pass it to the following step.
This works for any function h=h(T) and cp=cp(T) and there is no need of $h \propto T$ . The only need is $dh \propto dT$ . However, you should not intricate these two relations. These are different relations. The first one is an algebraic relation. The second one is a differential relation

December 16, 2019, 11:54		#2
clapointe Senior Member Join Date: Aug 2015 Posts: 494 Rep Power: 14	Check out this blog post : https://caefn.com/openfoam/temperature-calculation. Caelan stockzahn likes this.

April 20, 2020, 05:15	Temperature calculation	#4
fetc95 New Member Russian Federation Join Date: Apr 2020 Posts: 18 Rep Power: 6	Hello, I am concerned with the same problem, but the link provided above does not work. Could you explain how the temperature is obtained from the resulting field if you understand it? The meaning of "he" is more intricated for me be because of diffusivity term in the EEq.H of buoyantPimpleFoam. In the simplest case the energy equation in the buoyantPimpleFoam can simplified as Code: fvm::ddt(rho, he) + fvm::div(phi, he) - fvm::laplacian(alpha, he) == + fvOptions(rho, he) Let "he" is an energy which is calculated in some way. The term Code: fvm::laplacian(alpha, he) is seemed to be a heat conductivity, but why the "he" is under Laplacian and not the temperature? The heat diffusivity is proportional to the gradient of the temperature, not energy. Regards

May 4, 2020, 17:02		#15
stockzahn Member Join Date: Oct 2018 Location: France Posts: 34 Rep Power: 7	First of all: Sorry for the late reply - it was a strange week. Plus I obviously need much more time to process this mathematically. You are really good at this - I'm impressed. Thanks to your last post I think I've achieved a better understanding and thanks to your formulation I finally get how the calculation of works. If I get it right, for the polynomial approach, the calculation of the polynomial coefficents for the enthalpy is implemented based on the polynoms of the specific heat capacity. I assume similar pieces of code exist for the other models. Let's have a look at the calling path to the function T in thermoI.H (using the Newton 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 { if (T0 < 0) { FatalErrorInFunction << "Negative initial temperature T0: " << T0 << abort(FatalError); } scalar Test = T0; scalar Tnew = T0; scalar Ttol = T0tol_; 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; } In the enthalpy case it is called by the function THa (or THs) also defined in thermoI.H: 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 ); } The fifth argument is the specific heat capacity, whose function is called (thermo<Thermo, Type>::dFdT)(const scalar, const scalar), defined in the hPolynomialThermoI.H file: Code: template<class EquationOfState, int PolySize> inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::Cp ( const scalar p, const scalar T ) const { return CpCoeffs_.value(T) + EquationOfState::Cp(p, T); } Which in my interpretation returns a specific heat capacity value calculated by the polynomial model plus the summand "EquationOfState::Cp(p, T)", of which I'm not sure where to find it. I should test it, but my assumption for now is that it is defined in rPolynomialI.H: Code: template<class Specie> inline Foam::scalar Foam::rPolynomial<Specie>::Cp(scalar p, scalar T) const { return 0; } So I think I was able to understand the calculation of the temperature given an enthalpy (and pressure), but I still cannot see, how the discretized enthalpy diffusion could work Code: fvm::laplacian(turbulence->alphaEff(), he) . When the he-field is generated, at the end the function Ha (or Hs) in hPolynomialThermoI.H is called: Code: template<class EquationOfState, int PolySize> inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::Ha ( const scalar p, const scalar T ) const { return hCoeffs_.value(T) + EquationOfState::H(p, T); } As you explained very clearly, the Code: hCoeffs_ can be calculated from the Code: cpCoeffs_ , but in my understanding it still only holds a current field value. It does not contain the polynomial coefficients and I can't find a function to retrieve them. There is a similar issue with your suggestion to check with FDM: From an analytical point of view that all makes sense, but the discretized form (using just values for the diffusion coefficients kappa and alpha) seems to have the same problem like I thought the OF method should have. For the 1D-case the (explicit) FDM discretization of your last two equations (for ) should be $\rho c_p\frac{\partial T}{\partial t} - \frac{\partial}{\partial x}\left( \kappa \frac{\partial T}{\partial x} \right) = S$ $T^{n+1} = \frac{\Delta t}{\rho} \underbrace{\frac{\kappa}{c_p}}_{\alpha} \frac{T^n_{i+1} - 2 T^n_{i} + T^n_{i-1}}{\Delta x ^2} + T^n_i$ $\rho \frac{\partial h}{\partial t} - \frac{\partial}{\partial x}\left( \alpha \frac{\partial h}{\partial x} \right) = S$ $h^{n+1} = \frac{\Delta t}{\rho} \alpha \frac{h^n_{i+1} - 2 h^n_{i} + h^n_{i-1}}{\Delta x ^2} + h^n_i$ Using just "values" of the diffusion coefficients, the equations for temperature and enthalpy are identical they differ only by the source terms at the boundaries. That works if we assume the difference of the specific heat capacities for two neighboured cells negligible. Or did you have in mind something else?

May 5, 2020, 19:05		#18
fetc95 New Member Russian Federation Join Date: Apr 2020 Posts: 18 Rep Power: 6	Unfortunately, I cannot understand where you see that the properties passed to the solvers must be constant. In what part of the code do you see that the field which is passed as diffusion property (the 1th argument of the laplacian) cannot be varied from cell to cell? The alpha (or alphaEff, it is not important) is a volScalarField as you say. But you also say "with no any information". The type vol<Type>Field already means that this property can be varied from cell to cell. On the current solution step this field is already known and can be used in matrix coefficients of the linear equation system which is solved for defining the sought field. What also information is needed? May be it is needed to find how the diffusion coefficient alpha is calculated. In OpenFOAM descriptions, e.g., in Extended Code Guide in openfoam.com or in Programmers Guide from openfoam.org the term Code: laplacian(Gamma, G) means $\nabla \left( \Gamma \nabla G \right)$ or more exactly $\nabla \left( \Gamma \left( x,y,z \right) \nabla G \right)$ so the diffusion coefficient in general can be the function of coordinates, i.e., volField or surfaceField and the passed properties must not be constant. In addition, I am sure that there is no point for the developers to implement the diffusion term only for constant diffusion coefficient. Anyway, I believe that if it was so it would be mentioned in the programmers guide. But if you are doubt we need to analyze the code of fvm::laplacian operator in similarity with that we done for thermo models. However, I think that it is obviously. Concretely, I am sure that after the Gauss integration over the finite volume of the above term and discretization of the surface gradients in 1D case the OF gets: $\frac{1}{\Delta x^2} \left( \Gamma_{j+\frac{1}{2}}G_{j+1} - \left( \Gamma_{j+\frac{1}{2}} + \Gamma_{j-\frac{1}{2}} \right) G_j + \Gamma_{j-\frac{1}{2}}G_{j-1} \right)$ and not $\frac{\Gamma_j}{\Delta x^2} \left(G_{j+1} - 2 G_j + G_{j-1} \right)$ From the information provided by you I cannot see that there is some contradiction and that it is impossible. The implementation of the first way of discretization is not programmatically expensive, so I do not think that OF developers ignored that. Your second scheme is to be like this $\alpha^{num}_i = \alpha(p^n_i,T^n_i)$ $h^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} h^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) h^n_{i} + \alpha^{num}_{i-1/2}h^n_{i-1}}{\Delta x ^2} + h^n_i$ where $\alpha^{num}_{i\pm1/2}$ are values of field $\alpha^{num}$ interpolated onto cell faces. Here is no any need to know how the $\alpha^{num}$ was calculated. All we need is to know the cell values of the field $\alpha^{num}$ which were obtained when thermo.correct() was called on the previous step. The Laplacian(alpha,h) just will take this known field and calculate the diffusion term. Last edited by fetc95; May 5, 2020 at 19:52. Reason: adding

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
efficiency calculation using reference temperature	lostking18	CFX	0	September 10, 2019 13:08
Warning...Extrapolating on the temperature calculation	-nan	CONVERGE	7	April 21, 2017 11:10
Dynamic mesh adaption in parallel calculation	xh110120	FLUENT	1	March 12, 2016 08:05
where is the calculation of the temperature field	Tobi	OpenFOAM	1	July 30, 2012 10:40
temperature / enthalpy fields depending on type of fvPatchField	astein	OpenFOAM Programming & Development	0	June 28, 2010 07:10

April 20, 2020, 11:49		#5
stockzahn Member Join Date: Oct 2018 Location: France Posts: 34 Rep Power: 7	Dear fetc95, I currently dealt with exactly the same question. Unfortunately I wasn't able to find an answer I'm sure about, so I only can tell you my thoughts. Maybe you can verify or falsify them. - The diffusion term normally should be fvc::laplacian(kappa, T) - The disadvantage is that the diffusion cannot be calculated implicitly - The term fvm::laplacian(alpha, he) works if the enthalpy/internal energy is proportional to the temperature; only the factor Cp^-1/Cv^-1 (specific heat capacity), which is contained in alpha, must be taken into account (in OF alpha = kappa/Cp) - This works for simple models, e.g. perfect gas - In my case I integrated a property library (of helium); at some states helium has an inverse correlation between temperature and enthalpy/internal energy, therefore I think I must not use the implicit term, but have to use fvc::laplacian(kappa, T) which increase instability To sum up: In my understanding the implicit term fvm::laplacian(alpha, he) can be used, when 1) T proportional to he 2) Cp * T = h \| Cv * T = e otherwise the code has to be changed to calculate the diffusion explicitly from fvc::laplacian(kappa, T) What do you think?

April 24, 2020, 19:35		#10
fetc95 New Member Russian Federation Join Date: Apr 2020 Posts: 18 Rep Power: 6	Dear stockzahn Could you hint me what procedure do you use in your solver to get the temperature for using it in the equation?

May 6, 2020, 12:17		#19
stockzahn Member Join Date: Oct 2018 Location: France Posts: 34 Rep Power: 7	[QUOTE=fetc95;768930] Unfortunately, I cannot understand where you see that the properties passed to the solvers must be constant. [QUOTE/] I don't know why I can't express my self adequately. I didn't mean that the entire field has the same value, but the single values of the cells (which of course can differ) are passed as numerical value not as functions. But I see that you can think that, since I wrote the equation in my last post as I did (and not like in my (overall) post #7). I try to be more precise. I try again: If values are passed and not functions then $T^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} T^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) T^n_{i} + \alpha^{num}_{i-1/2}T^n_{i-1}}{\Delta x ^2} + T^n_i$ corresponds to $h^{n+1} = \frac{\Delta t}{\rho} \frac{\alpha^{num}_{i+1/2} h^n_{i+1} - \left( \alpha^{num}_{i+1/2} + \alpha^{num}_{i-1/2} \right) h^n_{i} + \alpha^{num}_{i-1/2}h^n_{i-1}}{\Delta x ^2} + h^n_i$ which only works, if $T \propto h$ or $\Delta T \propto \Delta h$ and that's only valid, if the change of over the temperature range $T_{i-1}$ - $T_{i+1}$ including $T^{n+1}$ is negligible. Else the enthalpy (evolution) does not represent the temperature (evolution).