Derivation of the Temperature Equation

Tobi · March 7, 2017, 13:10

Dear all,

I have a question to you all and hope that someone can clarify everything. What I am calculating is nothing special. It is just the temperature distribution inside a solid. Basically it is just the Laplace equation with the time derivative. However, I will show you how I got to the equation and what I figured out which I do not really understand.

Enthalpy equation
I started the derivation of the temperature equation by using the common enthalpy equation:

$\frac{\partial(\rho h)}{\partial t} + \nabla \bullet (\rho \textbf{U} h) = -\nabla \textbf{q} + S_h$

Based on the fact that I have a solid, there is no convection. In addition, further source terms

are not considered. Only the thermal conduction $\textbf{q}=-\lambda \nabla T$ is considered. Thus we get to the simple equation:

$\frac{\partial(\rho h)}{\partial t} = \nabla (\lambda \nabla T)$

The problem that I had is to replace the enthalpy with the temperature. We know that:

$h = \int_{T_{ref}}^{T} c_p \mathrm{d}T$

Okay, now I split the time derivative using the product rule:

$h\frac{\partial\rho}{\partial t} + \rho\frac{\partial h}{\partial t} = \nabla (\lambda \nabla T)$

Now I extended the second term as:

$\rho\frac{\partial h}{\partial t} = \rho\frac{\partial h}{\partial T}\frac{\partial T}{\partial t}$

I think we can change the partial derivative to finite differences - that should not be a big deal ( $\lim_{\partial t \to \Delta t}$ ) but please correct me if I am wrong here:

$\rho\frac{\partial h}{\partial T}\frac{\partial T}{\partial t} \longrightarrow \rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t}$

inserting the definition of the enthalpy we get:

$\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$

Summing up I got:

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

Inserting it to the enthalpy equation we get the following:

$h\frac{\partial\rho}{\partial t} + \rho c_p \frac{\partial T}{\partial t} = \nabla (\lambda \nabla T)$

The first term vanishes based on the continuity equation which should hold always (right?):
(Note... I was writing the continuity equation wrong - corrected now)

$\frac{\partial\rho}{\partial t} + \nabla \bullet (\rho \textbf{U}) = 0$

Based on the fact that there is no velocity $\textbf{U} = 0$ , all terms vanishes. To be more precise, the convective term vanishes and we get:

$\frac{\partial \rho}{\partial t} = 0$

Therefore, the first term in the new enthalpy equation vanishes and we get:

$\rho c_p \frac{\partial T}{\partial t} = \nabla (\lambda \nabla T)$

Till now, I did not assume that the density neither the heat capacity or the thermal conductivity is constant. So the first question is, is that derivation correct?

Implementation
I made simple checks in order to check out the temperature distributions. What I did is the following. I solved the three different kind of equations:

Code:

//- No convergence after 10 Outer loops
//  After the third outer loop, the initial residual stays constant (~ 0.0042)

fvScalarMatrix TEqn
(
    rho*cp*fvm::ddt(T)
    ==
    fvm::laplacian(k, T)
);

TEqn.solve()
Info<< "T[10] -> " << T[10] << endl;

Code:

//- No convergence after 10 Outer loops
//  After the third outer loop, the initial residual stays constant (~ 0.0042)
//  Same as above

fvScalarMatrix TEqn
(
    fvm::ddt(rho*cp, T)
    ==
    fvm::laplacian(k, T)
);

TEqn.solve()
Info<< "T[10] -> " << T[10] << endl;

Code:

//- Convergence after 4 iterations to 1e-10

fvScalarMatrix TEqn
(
    fvm::ddt(T)
    ==
    fvm::laplacian(k/rho/cp, T)
);

TEqn.solve()
Info<< "T[10] -> " << T[10] << endl;

More informations are:

Simple simulation (cavity case, left fixed T (300), right fixed T (850), bottom fixed T (350)
2D Mesh
10 Outer loops
the material properties are temperature depended (polynomials)

The second question that I have is, I am allowed to divide the whole equation by $c_p, \rho$ and put it inside the Laplace term? From a mathematical point I am not sure. Okay, these quantities are temperature depended and the Laplace operator (or gradient/divergence) is space operator. However, based on the temperature distribution these quantities are somehow space depended too.

The third question is, why are the first and second implementation not converging while the third one does?

Any help would be highly appreciated and it would be nice if someone can check my derivation. I am really not sure about that.

Thank you in advance,
Kind regards,

Tobi

arjun · March 7, 2017, 14:15

A . B = 0 and A = 0

does not imply that B = 0 too.

so i do not understand how you get d rho / d t = 0 from u = 0.

Zeppo · March 7, 2017, 16:39

You didn't post the results you've gotten from the code runs

Code:

Info<< "T[10] -> " << T[10] << endl;

They should be different in case rho and/or Cp are temperature dependent. Are the results different from each other?

Look at the source code.

Code:

rho*cp*fvm::ddt(T)

is equivalent to

Code:

new GeometricField<Type, fvPatchField, volMesh>
(
    io,
    rho*cp * 1.0/mesh().time().deltaT() *(T - T.oldTime())
)

And

Code:

fvm::ddt(rho*cp, T)

is equivalent to

Code:

new GeometricField<Type, fvPatchField, volMesh>
(
    io,
    1.0/mesh().time().deltaT() * (rho*cp*T - (rho*cp).oldTime()*T.oldTime())
)

These two snippets result in two different fvMatrixes.

Quote:

The third question is, why are the first and second implementation not converging while the third one does?

Why do you think your first and secondexamples are not converged? A low absolute level of residuals doesn't imply that the solution is diverged. watching some monitors should help better to say if it's true or not.

Tobi · March 7, 2017, 17:27

Hey Arjun,

you are right and I was too stupid to write the correct continuity equation. I corrected it now and you can see that the derivation should be correct. Thanks for mentioning it (indirectly).

Edit based on the reply from Zeppo

You are right, again I made a mistake. I was using wrong words. I did not want to say converged and not converged but the matrix system (the residual error) does not get smaller anymore (comparing the first two implementations and the last one). That 's what I wanted to say. Tomorrow I will make an analysis of the absolute values and post it. As far as I have it in mind, #1 and #2 were similar where #3 was different (but only ~ 0.3K).

Zeppo · March 8, 2017, 10:17

Tobias, your final equation

$\rho c_p \frac{\partial T}{\partial t} = \nabla (\lambda \nabla T)$

is not correct. The correct one is:

$c_p\frac{\partial (\rho T)}{\partial t} = \nabla \cdot (\lambda \nabla T)$

which is equivalent to

$\frac{\partial (\rho c_p T)}{\partial t} = \nabla \cdot (\lambda \nabla T) + \rho T \frac{\partial c_p}{\partial t}$

The proof is given in the book by Moukalled at alias on the pages 62-63.

Tobi · March 8, 2017, 10:40

Dear Sergei,

thanks for the hint. I will check it out, perfect! Thank you again.

arjun · March 8, 2017, 11:20

Why can not you use the very first equation you wrote, the enthalpy equation?

Tobi · March 8, 2017, 11:40

Hi Arjun,

I could do but I am just interested in the temperature distribution in my solid body and I do not want to recalculate the temperature based on the enthalpy. Well FOAM offers the recalculation but I wanted to keep my solver as simple as possible. Never checked out the code how FOAM recalculate T from h and vice versa.

arjun · March 8, 2017, 12:10

Quote:

Originally Posted by Tobi

Hi Arjun,

I could do but I am just interested in the temperature distribution in my solid body and I do not want to recalculate the temperature based on the enthalpy. Well FOAM offers the recalculation but I wanted to keep my solver as simple as possible. Never checked out the code how FOAM recalculate T from h and vice versa.

Its alright.

I use the equation based on total energy and total enthalpy in wildkatze and directly calculate the temperature from it without recalculating from enthalpy.

What you do is cheaper though, no doubt. I did that way because it can accommodate most of the things ie solids u, p etc terms gone, in case of fluids taken care. The same could be used in combustion (i am in process of coding now), same could be used with user defined gas law for example.

PS: Last 1 month I am learning and trying to code combustion it is fun :-) So I am getting into that stiff ODE problem now.

Tobi · March 8, 2017, 13:35

Hi Arjun,

Offtopic
I am dealing with combustion since 3 years. Maybe you know that I am developing a flamelet creator. The only missing part is solving the Stiff ODE. Unfortunately I made my own C++ project to learn much more in OOP but now I cannot use the ODE solvers from FOAM like, SEULEX, RODOS etc. Well, I hope that one day I can implement into my code. If you want to contribute, you are welcomed.

Topic
I agree that in fluids it would be nice calculating the enthalpy. As you already said, I just calculate the temperature because it is cheaper and I do not solve other stuff (only thermal stress). However, if someone is interested, I made a simple test case (attached - you have to modify it to e.g. laplacianFoam).

However, I checked the different temperature implementations named 1, 2, 3 , 4 and 5 in the charts. The implementation was as follows:

Code:

// Calculated marked with 1 
rho*cp*fvm::ddt(T) == laplacian(k, T)

// Calculated marked with 2
cp*fvm::ddt(rho, T) == laplacian(k, T)

// Calculated marked with 3
fvm::ddt(cp*rho, T) == laplacian(k, T)

// Calculated marked with 4
fvm::ddt(cp*rho, T) - rho*T*fvc::ddt(cp) == laplacian(k, T)

// Calculated marked with 5
fvm::ddt(T) == laplacian(k/rho/cp, T)

The result is given in the two plots. The first plot gives the maximum temperature in the domain while the other one shows the residuals during the outer corrector loops (just a cut of the calculation). It is interesting that using the product rule (as Zappo mentioned) let the simulation diverge. However, I also checked SuSpand Sp for that term (based on the fact that it is always negative, there is no need for SuSp or Sp. However, another fact is, that 1, 2, 3 are similar in the values while 3 does not converge based on the residuals after 10 outer corrections and the residuals are even larger.

Thank you Zappo for the hint. I will use implementation (2) for now.

arjun · March 8, 2017, 13:50

Quote:

Originally Posted by Tobi

Hi Arjun,

Offtopic
I am dealing with combustion since 3 years. Maybe you know that I am developing a flamelet creator. The only missing part is solving the Stiff ODE. Unfortunately I made my own C++ project to learn much more in OOP but now I cannot use the ODE solvers from FOAM like, SEULEX, RODOS etc. Well, I hope that one day I can implement into my code. If you want to contribute, you are welcomed.

.

I intend to make the VODE based solver free for others to use as library, it is C++ implementation from cvode that was available. What would be different is that in year 2013 I created a faster LU solver for dense matrices. So wish to see if that gives speed up here. (Compared to starccm+ internal direct solver it was 30 times faster for 300x300 matrices and efficiency gape increases as the size increases).

I am at the moment working on finite rate species combustion and now I start to work on ERENA algorithm and VODE based stiff solver (after ERENA) (will compare the two on cost of calculation).

I created a project called firekatze around these things like flamelets and other combustion related things (I added chemkin import and transport properties import etc).
Shall provide you then species properties from temperature and pressure.

At the moment the focus is to get this combustion thing work on finite rate and it shall be faster. Made lots of change to accomodate it in FVUS (like explicit navier stokes for example)

PS: I learned all that in last 30 days, i remember i started on 7th feb and its 8th now. So huge amount of code i did in this month, still going on.

Zeppo · March 8, 2017, 15:00

Quote:

Originally Posted by Tobi

However, I checked the different temperature implementations named 1, 2, 3 , 4 and 5 in the charts. The implementation was as follows:

Code:

// Calculated marked with 1 
rho*cp*fvm::ddt(T) == laplacian(k, T)

// Calculated marked with 2
cp*fvm::ddt(rho, T) == laplacian(k, T)

// Calculated marked with 3
fvm::ddt(cp*rho, T) == laplacian(k, T)

// Calculated marked with 4
fvm::ddt(cp*rho, T) - rho*T*fvc::ddt(cp) == laplacian(k, T)

// Calculated marked with 5
fvm::ddt(T) == laplacian(k/rho/cp, T)

The result is given in the two plots. The first plot gives the maximum temperature in the domain while the other one shows the residuals during the outer corrector loops (just a cut of the calculation). It is interesting that using the product rule (as Zappo mentioned) let the simulation diverge.

I am surprised by the fact that #4 is so far off the mark. How about decreasing the timestep (say, by the factor of 2) and addidng some more outer iteration (say, 2 times more) on every timestep. Just to be sure the residuals get plain and solution don't evolve at the timestep anymore.

Tobi · March 8, 2017, 16:14

I was also thinking about decreasing the time step and increasing the outer loops. But I wanted to have Same settings. Well, tomorrow I will Check it out.

Sent from my HTC One mini using CFD Online Forum mobile app

Tobi · March 9, 2017, 03:19

I posted wrong results here because the last check yesterday was with constant values (cp, rho, k). Thats the reason why 4 was not blowing up. New results are in preparation.

agustinvo · March 9, 2017, 06:34

Hello Tobi,
some time ago I was considering the same sitation but for fluids, if I should include rho and Cp inside or outside the derivatives. It's a long time I don't think about it, but I should!

In the other hand, I think you had a typo here:
$\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\partial t} = c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow c_p \frac{\partial T}{\partial t}$

rho disappeared misteriously:
$\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\partial t} =\rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$

Tobi · March 9, 2017, 07:05

Yes you are right. The rho should be considered. Typing all the equations in here is sometimes a bit ... wired. I will change it. However, in the literature Zappo is mentioning, you will find the derivation. However, I don't know where the mistake in my derivation is. The difference that I can point out at the moment is that Moukalled et al. use the total derivative which then lead the the density inside. What I mean is, that I was doing that:

$\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\partial t} =\rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$

While Moukalled was getting the total derivative:

$\rho c_p \frac{\mathrm{D} T}{\mathrm{D} t}$
It is clear why the density is inside the derivative if we rearrange the equation without the total derivative. I don't know, maybe I was missing something or my assumption I made with the lim is not correct. I might miss some important thing.

By the way there is a typo on page 62 in equation 3.67. DP -> Dp.

Summing up:

Moukalled et al. replaced:

$\mathrm{d}h = c_p \mathrm{d}T$

into the total derivative of the enthalpy:

$\rho \frac{\mathrm{D}h}{\mathrm{D t}} = \rho c_p \frac{\mathrm{D}T}{\mathrm{D t}}$

and I used the lim to go from finite differences to partial differential stuff. So I think there I made a mistake. To be honest, I am not sure if am allowed to do what I did. Would be nice if someone can check it and tell the mistake.

Note:
I changed the equations above related to the mistake I made (rho missing).

Remark: (not correct - see my next post. Everything is fine)
After checking the first post, the difference comes from this line:

$\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$

Especially this one:

$\rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$

Which should be wrong based on Moukalled et al. because he related it to the total Derivative:

$c_p \mathrm{d} T \longrightarrow \rho c_p \frac{\mathrm{D} T}{\mathrm{D} t}$

arjun · March 9, 2017, 07:57

Quote:

Originally Posted by Tobi

Hi Arjun,

Offtopic
I am dealing with combustion since 3 years. Maybe you know that I am developing a flamelet creator. The only missing part is solving the Stiff ODE. Unfortunately I made my own C++ project to learn much more in OOP but now I cannot use the ODE solvers from FOAM like, SEULEX, RODOS etc. Well, I hope that one day I can implement into my code. If you want to contribute, you are welcomed.

.

Forgive me for hijacking your thread for 1 or two more posts. When you get time, please let me know the interface that you would want the ode solver and other things that you may want to use.

If possible just name the functions something like this
void Function( arguements ) ... for the things that you would want so I would create them into library. That way it would save you lot of headache.

Also know that firekatze can load chemkin, so if would want it would give you other things like diffusion coefficients of species, enthalphy , conductivity, viscosity etc of species etc etc.
So when you get time, just prepare a note what things you wish and library shall have it.

PS: At the moment i am using chemkin and transport data produced from
http://akrmys.com/KUCRS/index.htm.en

Tobi · March 9, 2017, 08:30

Offtopic:

Hi Arjun, I made my own code that loads chemistry, transport and thermodynamic data. The code reads a dictionary where the user can specify the fuel and oxidizer stream and other solver related stuff. After that, I construct the polynomials for fitting (viscosity, thermal conductivity, binary diffusion etc.). However, the user can specify the calculation of the gas kinetic equation and if he want to use fitting procedure or direct calculation. Then I build my flamelet class with all information which then should solve everything. Everything is already implemented, TROE, LOW, ENHANCED, SRI ... but the only thing that is missing is the solving of the chemistry. I was thinking of using sub-time stepping, like flamelet equation dt and then solve chemistry with an ODE in order to reach that dt. After that go to the next time step (I think it is a common approach). At the moment (or 3 weeks ago) I started to implement simple Euler Integration which is easy but not really powerful. However, doing that, I would have the basis of other solver implementations. But like always - no time. So the project is stlil at the moment. If you are interested: http://afc.holzmann-cfd.de/ but it is outdated. If you want, I can upload the latest doxygen stuff today evening.

Topic

Just want to inform you that everything is fine and the equation I derived at the beginning is correct within my assumptions. As Zeppo mentioned, the temperature equation is given in Moukalled et al. However, as I pointed out above, I got a different equation which is in fact correct because I do not have any velocities. Therefore, taking the equation from Moukalled et al. and applying my assumptions (no motion -> no velocity):

$c_p \left[ \frac{\partial \rho T}{\partial t} + \nabla \bullet (\rho \textbf{U} T) \right] = \nabla \bullet (k \nabla T)$

$c_p \left[ \rho \frac{\partial T}{\partial t} + T \frac{\partial \rho}{\partial t} + \rho \textbf{U} \bullet \nabla T + T \nabla \bullet (\rho \textbf{U}) \right] = \nabla \bullet (k \nabla T)$

$c_p \left[ \rho \frac{\partial T}{\partial t} + \rho \textbf{U} \bullet \nabla T + T \left\{ \underbrace{\frac{\partial \rho}{\partial t} + \nabla \bullet (\rho \textbf{U})}_{\text{continuity}}\right\} \right] = \nabla \bullet (k \nabla T)$

$c_p \rho \frac{\partial T}{\partial t} + \underbrace{c_p \rho \textbf{U} \bullet \nabla T}_{\text{\textbf{U} is zero - no velocity}} = \nabla \bullet (k \nabla T)$

$c_p \rho \frac{\partial T}{\partial t} = \nabla \bullet (k \nabla T)$

Thats the reason why the calculations of the temperature of 1 and 3 are almost similar. I hope that this is the proof and everybody would agree.

hasan_shetabivash · March 9, 2017, 08:39

Dear Tobias,
First of all assuming constant density and heat capacity for solids is reasonable and the energy equation you've obtained for solid is acceptable (see the book by Moukalled page 65 Eq. 3.80). The difference between your mathematical operations with those of Moukalled is the fact that zero velocity was your foremost assumption and you derived all the equations based on the zero velocity in the domain. However, in the book the zero velocity was the last assumption which lead to your final equation. The only thing that you need to consider is non-linearity that is added to the equation by dependence of Cp to temperature. If Cp(T) is changing dramatically you need an appropriate method to deal with the nonlinear problem.

Cheers,

cfdopenfoam · April 10, 2018, 06:48

Hi Tobi,

I am facing the same situation here, planning to solve a temperature equation as yours but with flow U.

Is it possible for u to share your solution/code here setting us a valuable example?

I think foamer struggling with the same situation will appreciate great thanks to you.

Best,
Karekle

March 7, 2017, 13:10	Derivation of the Temperature Equation	#1
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Dear all, I have a question to you all and hope that someone can clarify everything. What I am calculating is nothing special. It is just the temperature distribution inside a solid. Basically it is just the Laplace equation with the time derivative. However, I will show you how I got to the equation and what I figured out which I do not really understand. Enthalpy equation I started the derivation of the temperature equation by using the common enthalpy equation: $\frac{\partial(\rho h)}{\partial t} + \nabla \bullet (\rho \textbf{U} h) = -\nabla \textbf{q} + S_h$ Based on the fact that I have a solid, there is no convection. In addition, further source terms are not considered. Only the thermal conduction $\textbf{q}=-\lambda \nabla T$ is considered. Thus we get to the simple equation: $\frac{\partial(\rho h)}{\partial t} = \nabla (\lambda \nabla T)$ The problem that I had is to replace the enthalpy with the temperature. We know that: $h = \int_{T_{ref}}^{T} c_p \mathrm{d}T$ Okay, now I split the time derivative using the product rule: $h\frac{\partial\rho}{\partial t} + \rho\frac{\partial h}{\partial t} = \nabla (\lambda \nabla T)$ Now I extended the second term as: $\rho\frac{\partial h}{\partial t} = \rho\frac{\partial h}{\partial T}\frac{\partial T}{\partial t}$ I think we can change the partial derivative to finite differences - that should not be a big deal ( $\lim_{\partial t \to \Delta t}$ ) but please correct me if I am wrong here: $\rho\frac{\partial h}{\partial T}\frac{\partial T}{\partial t} \longrightarrow \rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t}$ inserting the definition of the enthalpy we get: $\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$ Summing up I got: $\rho\frac{\partial h}{\partial T}\frac{\partial T}{\partial t} = \rho c_p \frac{\partial T}{\partial t}$ Inserting it to the enthalpy equation we get the following: $h\frac{\partial\rho}{\partial t} + \rho c_p \frac{\partial T}{\partial t} = \nabla (\lambda \nabla T)$ The first term vanishes based on the continuity equation which should hold always (right?): (Note... I was writing the continuity equation wrong - corrected now) $\frac{\partial\rho}{\partial t} + \nabla \bullet (\rho \textbf{U}) = 0$ Based on the fact that there is no velocity $\textbf{U} = 0$ , all terms vanishes. To be more precise, the convective term vanishes and we get: $\frac{\partial \rho}{\partial t} = 0$ Therefore, the first term in the new enthalpy equation vanishes and we get: $\rho c_p \frac{\partial T}{\partial t} = \nabla (\lambda \nabla T)$ Till now, I did not assume that the density neither the heat capacity or the thermal conductivity is constant. So the first question is, is that derivation correct? Implementation I made simple checks in order to check out the temperature distributions. What I did is the following. I solved the three different kind of equations: Code: //- No convergence after 10 Outer loops // After the third outer loop, the initial residual stays constant (~ 0.0042) fvScalarMatrix TEqn ( rhocpfvm::ddt(T) == fvm::laplacian(k, T) ); TEqn.solve() Info<< "T[10] -> " << T[10] << endl; Code: //- No convergence after 10 Outer loops // After the third outer loop, the initial residual stays constant (~ 0.0042) // Same as above fvScalarMatrix TEqn ( fvm::ddt(rhocp, T) == fvm::laplacian(k, T) ); TEqn.solve() Info<< "T[10] -> " << T[10] << endl; Code: //- Convergence after 4 iterations to 1e-10 fvScalarMatrix TEqn ( fvm::ddt(T) == fvm::laplacian(k/rho/cp, T) ); TEqn.solve() Info<< "T[10] -> " << T[10] << endl; More informations are: Simple simulation (cavity case, left fixed T (300), right fixed T (850), bottom fixed T (350) 2D Mesh 10 Outer loops the material properties are temperature depended (polynomials) The second question that I have is, I am allowed to divide the whole equation by $c_p, \rho$ and put it inside the Laplace term? From a mathematical point I am not sure. Okay, these quantities are temperature depended and the Laplace operator (or gradient/divergence) is space operator. However, based on the temperature distribution these quantities are somehow space depended too. The third question is, why are the first and second implementation not converging while the third one does? Any help would be highly appreciated and it would be nice if someone can check my derivation. I am really not sure about that. Thank you in advance, Kind regards, Tobi Kummi likes this. __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; March 9, 2017 at 07:10. Reason: Wrong continuity equation / change #2*

March 8, 2017, 10:17		#5
Zeppo Senior Member Sergei Join Date: Dec 2009 Posts: 261 Rep Power: 22	Tobias, your final equation $\rho c_p \frac{\partial T}{\partial t} = \nabla (\lambda \nabla T)$ is not correct. The correct one is: $c_p\frac{\partial (\rho T)}{\partial t} = \nabla \cdot (\lambda \nabla T)$ which is equivalent to $\frac{\partial (\rho c_p T)}{\partial t} = \nabla \cdot (\lambda \nabla T) + \rho T \frac{\partial c_p}{\partial t}$ The proof is given in the book by Moukalled at alias on the pages 62-63. Tobi, Zhiheng Wang, Kummi and 1 others like this.

March 8, 2017, 10:40		#6
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Dear Sergei, thanks for the hint. I will check it out, perfect! Thank you again. __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; October 14, 2017 at 08:03.

March 8, 2017, 11:40		#8
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Hi Arjun, I could do but I am just interested in the temperature distribution in my solid body and I do not want to recalculate the temperature based on the enthalpy. Well FOAM offers the recalculation but I wanted to keep my solver as simple as possible. Never checked out the code how FOAM recalculate T from h and vice versa. __________________ Keep foaming, Tobias Holzmann

March 8, 2017, 16:14		#13
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	I was also thinking about decreasing the time step and increasing the outer loops. But I wanted to have Same settings. Well, tomorrow I will Check it out. Sent from my HTC One mini using CFD Online Forum mobile app __________________ Keep foaming, Tobias Holzmann

March 7, 2017, 14:15		#2
arjun Senior Member Arjun Join Date: Mar 2009 Location: Nurenberg, Germany Posts: 1,291 Rep Power: 35	A . B = 0 and A = 0 does not imply that B = 0 too. so i do not understand how you get d rho / d t = 0 from u = 0. Tobi likes this.

March 7, 2017, 17:27		#4
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Hey Arjun, you are right and I was too stupid to write the correct continuity equation. I corrected it now and you can see that the derivation should be correct. Thanks for mentioning it (indirectly). Edit based on the reply from Zeppo You are right, again I made a mistake. I was using wrong words. I did not want to say converged and not converged but the matrix system (the residual error) does not get smaller anymore (comparing the first two implementations and the last one). That 's what I wanted to say. Tomorrow I will make an analysis of the absolute values and post it. As far as I have it in mind, #1 and #2 were similar where #3 was different (but only ~ 0.3K). __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; March 8, 2017 at 04:31. Reason: Updated based on wrong continuity equation in the first post

March 8, 2017, 11:20		#7
arjun Senior Member Arjun Join Date: Mar 2009 Location: Nurenberg, Germany Posts: 1,291 Rep Power: 35	Why can not you use the very first equation you wrote, the enthalpy equation?

March 9, 2017, 03:19		#14
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	I posted wrong results here because the last check yesterday was with constant values (cp, rho, k). Thats the reason why 4 was not blowing up. New results are in preparation. __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; March 9, 2017 at 06:02. Reason: Used constant values - wrong results

March 9, 2017, 06:34		#15
agustinvo Senior Member Agustín Villa Join Date: Apr 2013 Location: Alcorcón Posts: 314 Rep Power: 15	Hello Tobi, some time ago I was considering the same sitation but for fluids, if I should include rho and Cp inside or outside the derivatives. It's a long time I don't think about it, but I should! In the other hand, I think you had a typo here: $\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\partial t} = c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow c_p \frac{\partial T}{\partial t}$ rho disappeared misteriously: $\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\partial t} =\rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$ Tobi likes this.

March 9, 2017, 07:05		#16
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Yes you are right. The rho should be considered. Typing all the equations in here is sometimes a bit ... wired. I will change it. However, in the literature Zappo is mentioning, you will find the derivation. However, I don't know where the mistake in my derivation is. The difference that I can point out at the moment is that Moukalled et al. use the total derivative which then lead the the density inside. What I mean is, that I was doing that: $\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\partial t} =\rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$ While Moukalled was getting the total derivative: $\rho c_p \frac{\mathrm{D} T}{\mathrm{D} t}$ It is clear why the density is inside the derivative if we rearrange the equation without the total derivative. I don't know, maybe I was missing something or my assumption I made with the lim is not correct. I might miss some important thing. By the way there is a typo on page 62 in equation 3.67. DP -> Dp. Summing up: Moukalled et al. replaced: $\mathrm{d}h = c_p \mathrm{d}T$ into the total derivative of the enthalpy: $\rho \frac{\mathrm{D}h}{\mathrm{D t}} = \rho c_p \frac{\mathrm{D}T}{\mathrm{D t}}$ and I used the lim to go from finite differences to partial differential stuff. So I think there I made a mistake. To be honest, I am not sure if am allowed to do what I did. Would be nice if someone can check it and tell the mistake. Note: I changed the equations above related to the mistake I made (rho missing). Remark: (not correct - see my next post. Everything is fine) After checking the first post, the difference comes from this line: $\rho\frac{\mathrm{d} h}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho\frac{\mathrm{d} \left[\int_{T_{ref}}^{T} c_p \mathrm{d}T\right]}{\mathrm{d} T}\frac{\mathrm{d} T}{\mathrm{d} t} = \rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$ Especially this one: $\rho c_p \frac{\mathrm{d} T}{\mathrm{d} t} \longrightarrow \rho c_p \frac{\partial T}{\partial t}$ Which should be wrong based on Moukalled et al. because he related it to the total Derivative: $c_p \mathrm{d} T \longrightarrow \rho c_p \frac{\mathrm{D} T}{\mathrm{D} t}$ __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; February 16, 2018 at 08:54. Reason: Remark might be not be decribed well

March 9, 2017, 08:30		#18
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Offtopic: Hi Arjun, I made my own code that loads chemistry, transport and thermodynamic data. The code reads a dictionary where the user can specify the fuel and oxidizer stream and other solver related stuff. After that, I construct the polynomials for fitting (viscosity, thermal conductivity, binary diffusion etc.). However, the user can specify the calculation of the gas kinetic equation and if he want to use fitting procedure or direct calculation. Then I build my flamelet class with all information which then should solve everything. Everything is already implemented, TROE, LOW, ENHANCED, SRI ... but the only thing that is missing is the solving of the chemistry. I was thinking of using sub-time stepping, like flamelet equation dt and then solve chemistry with an ODE in order to reach that dt. After that go to the next time step (I think it is a common approach). At the moment (or 3 weeks ago) I started to implement simple Euler Integration which is easy but not really powerful. However, doing that, I would have the basis of other solver implementations. But like always - no time. So the project is stlil at the moment. If you are interested: http://afc.holzmann-cfd.de/ but it is outdated. If you want, I can upload the latest doxygen stuff today evening. Topic Just want to inform you that everything is fine and the equation I derived at the beginning is correct within my assumptions. As Zeppo mentioned, the temperature equation is given in Moukalled et al. However, as I pointed out above, I got a different equation which is in fact correct because I do not have any velocities. Therefore, taking the equation from Moukalled et al. and applying my assumptions (no motion -> no velocity): $c_p \left[ \frac{\partial \rho T}{\partial t} + \nabla \bullet (\rho \textbf{U} T) \right] = \nabla \bullet (k \nabla T)$ $c_p \left[ \rho \frac{\partial T}{\partial t} + T \frac{\partial \rho}{\partial t} + \rho \textbf{U} \bullet \nabla T + T \nabla \bullet (\rho \textbf{U}) \right] = \nabla \bullet (k \nabla T)$ $c_p \left[ \rho \frac{\partial T}{\partial t} + \rho \textbf{U} \bullet \nabla T + T \left\{ \underbrace{\frac{\partial \rho}{\partial t} + \nabla \bullet (\rho \textbf{U})}_{\text{continuity}}\right\} \right] = \nabla \bullet (k \nabla T)$ $c_p \rho \frac{\partial T}{\partial t} + \underbrace{c_p \rho \textbf{U} \bullet \nabla T}_{\text{\textbf{U} is zero - no velocity}} = \nabla \bullet (k \nabla T)$ $c_p \rho \frac{\partial T}{\partial t} = \nabla \bullet (k \nabla T)$ Thats the reason why the calculations of the temperature of 1 and 3 are almost similar. I hope that this is the proof and everybody would agree. hwangpo, Thamali, mbookin and 1 others like this. __________________ Keep foaming, Tobias Holzmann

March 9, 2017, 08:39		#19
hasan_shetabivash New Member Hasan Shetabivash Join Date: Jan 2017 Location: Montreal Posts: 17 Rep Power: 12	Dear Tobias, First of all assuming constant density and heat capacity for solids is reasonable and the energy equation you've obtained for solid is acceptable (see the book by Moukalled page 65 Eq. 3.80). The difference between your mathematical operations with those of Moukalled is the fact that zero velocity was your foremost assumption and you derived all the equations based on the zero velocity in the domain. However, in the book the zero velocity was the last assumption which lead to your final equation. The only thing that you need to consider is non-linearity that is added to the equation by dependence of Cp to temperature. If Cp(T) is changing dramatically you need an appropriate method to deal with the nonlinear problem. Cheers, Tobi likes this.

April 10, 2018, 06:48	cloud u pls share your solution here	#20
cfdopenfoam Member Karelke Yu Join Date: Dec 2014 Posts: 96 Rep Power: 12	Hi Tobi, I am facing the same situation here, planning to solve a temperature equation as yours but with flow U. Is it possible for u to share your solution/code here setting us a valuable example? I think foamer struggling with the same situation will appreciate great thanks to you. Best, Karekle

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Static Temperature / Opening Temperature	JulianP	CFX	12	April 10, 2019 19:00
Multiphase flow - incorrect velocity on inlet	Mike_Tom	CFX	6	September 29, 2016 02:27
Constant velocity of the material	Sas	CFX	15	July 13, 2010 09:56
temperature equation for one phase in multiphase solver	romant	OpenFOAM Programming & Development	0	April 13, 2010 06:06
Adding temperature equation in settlingFoam	sachin	OpenFOAM Running, Solving & CFD	2	March 31, 2010 04:21