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

Derivation of the Temperature Equation

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

Like Tree15Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 7, 2017, 12:10
Default Derivation of the Temperature Equation
  #1
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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 S_h 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
Kummi likes this.
__________________
Keep foaming,
Tobias Holzmann

Last edited by Tobi; March 9, 2017 at 06:10. Reason: Wrong continuity equation / change #2
Tobi is offline   Reply With Quote

Old   March 7, 2017, 13:15
Default
  #2
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,278
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
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.
arjun is offline   Reply With Quote

Old   March 7, 2017, 15:39
Default
  #3
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21
Zeppo will become famous soon enough
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 likes this.
Zeppo is offline   Reply With Quote

Old   March 7, 2017, 16:27
Default
  #4
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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 03:31. Reason: Updated based on wrong continuity equation in the first post
Tobi is offline   Reply With Quote

Old   March 8, 2017, 09:17
Default
  #5
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21
Zeppo will become famous soon enough
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.
Zeppo is offline   Reply With Quote

Old   March 8, 2017, 09:40
Default
  #6
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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 07:03.
Tobi is offline   Reply With Quote

Old   March 8, 2017, 10:20
Default
  #7
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,278
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
Why can not you use the very first equation you wrote, the enthalpy equation?
arjun is offline   Reply With Quote

Old   March 8, 2017, 10:40
Default
  #8
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 8, 2017, 11:10
Default
  #9
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,278
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
Quote:
Originally Posted by Tobi View Post
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 likes this.
arjun is offline   Reply With Quote

Old   March 8, 2017, 12:35
Default
  #10
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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.


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.
Attached Images
File Type: png Residuals.png (18.2 KB, 77 views)
File Type: png Temperature.png (10.9 KB, 83 views)
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 8, 2017, 12:50
Default
  #11
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,278
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
Quote:
Originally Posted by Tobi View Post
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.
Tobi likes this.
arjun is offline   Reply With Quote

Old   March 8, 2017, 14:00
Default
  #12
Senior Member
 
Zeppo's Avatar
 
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21
Zeppo will become famous soon enough
Quote:
Originally Posted by Tobi View Post
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.
Zeppo is offline   Reply With Quote

Old   March 8, 2017, 15:14
Default
  #13
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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
Tobi is offline   Reply With Quote

Old   March 9, 2017, 02:19
Default
  #14
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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 05:02. Reason: Used constant values - wrong results
Tobi is offline   Reply With Quote

Old   March 9, 2017, 05:34
Default
  #15
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Alcorcón
Posts: 313
Rep Power: 15
agustinvo is on a distinguished road
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.
agustinvo is offline   Reply With Quote

Old   March 9, 2017, 06:05
Default
  #16
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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 07:54. Reason: Remark might be not be decribed well
Tobi is offline   Reply With Quote

Old   March 9, 2017, 06:57
Default
  #17
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,278
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
Quote:
Originally Posted by Tobi View Post
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
arjun is offline   Reply With Quote

Old   March 9, 2017, 07:30
Default
  #18
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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
Tobi is offline   Reply With Quote

Old   March 9, 2017, 07:39
Default
  #19
New Member
 
Hasan Shetabivash
Join Date: Jan 2017
Location: Montreal
Posts: 17
Rep Power: 12
hasan_shetabivash is on a distinguished road
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.
hasan_shetabivash is offline   Reply With Quote

Old   April 10, 2018, 05:48
Default cloud u pls share your solution here
  #20
Member
 
Karelke Yu
Join Date: Dec 2014
Posts: 96
Rep Power: 11
cfdopenfoam is on a distinguished road
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
cfdopenfoam is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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

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


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


All times are GMT -4. The time now is 17:14.