CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Heat transfer in liquid water suggestions for chioce of solver (

hkb September 17, 2008 10:49

Hi all I'm new to OpenFOAM as
Hi all
I'm new to OpenFOAM as well as CFD, and would like to ask for some advice on how to choose solvers. I'm interested in thermal fatigue in pipes, and hence I would like to model thermal transport in such systems. The medium is liquid water, and hence rho=const (or in some cases rho=rho(T)). I would eventually like to do RANS as well as LES. As far as I can see, the standard library offers either incompressible solvers with no heat transport, or compressible solver utilizing a perfect gas equation of state. Based on that it seems that my choices are:

1. Adding temperature transport to a (or several) incompressible solver.

2. creating a new equation of state, where rho=const (or rho=rho(T)), and use existing compressible solvers.

I am kindly asking for suggestion or comments in this matter. Maybe someone has done something similar? I should also add that apart from being a CFD beginner, my knowledge in c++ is also limited.

Best Regards

dpalko September 17, 2008 11:14

Hi Henrik, I suggest you st
Hi Henrik,

I suggest you start with incompressible solvers. The following is an interesting wiki link, how to add the energy equation to icoFoam.


grtabor September 17, 2008 11:15

Hi. I would choose option 1

I would choose option 1.; nonlinear density just adds to the complexity of solution, whereas an additional transport equation should be easy to include into turbFoam (or something similar).


hkb September 17, 2008 14:51

Hi Thanks for your rapid repl
Thanks for your rapid reply. I have started adding a temperature transport to turbFoam (based on this thread: I figured that since the flow affects the temperature but not the other way around, the temperature equation should be solved after the PISO loop, is that correct?

Best Regards

ameya September 17, 2008 15:21

Hi Henrik Even I am working
Hi Henrik

Even I am working on heat transfer in liquids and my work comprises of determining temperature fluctuations occuring due to mixing of non isothermal jets.

I am using heat transfer codes available in OpenFOAM. If you want to know which solvers to use go to : OpenFOAM<version>/applications/solvers/heattransfer

If you want rho = rho(T) i suggest you use OpenFOAM-1.5 . For 1.5 version you can download code for density variation using boussinesq approximation. you can use that code as it is or easily modify it for piecewise linear and piecewise polynomial variation of density

In case you need any further help feel free to ask

grtabor September 17, 2008 17:52

I had a look at the code at:
I had a look at the code at:

Am I right in thinking it misses off the viscous heating term??


hkb September 18, 2008 12:50

Yes, I think the viscous heati
Yes, I think the viscous heating is missing.
Based on 0Fall-2004/1CE85465-E605-4582-841B-24C3966E332A/0/energy.pdf I would write (ignoring div(u) terms)

ddt(T)+div(phi,T)-laplacian(alphaEff,T)- 1/(rho*Cv) tau && grad(U) == 0

This is essentially what is proposed in this thread: but for compressible flow. In that thread I also found this relation:

tau = mu * (gradU + gradU.T()) - (2.0/3.0 * mu * fvc::div(U)) * I

where I guess the second term drops out due to incompressibility.

Any thoughts on this?

Ameya: Are you referring to the Boussinesq solver? I looked at it a while back (from the OF-1.4.1-dev distribution) and as I recall it was only for laminar flow, while I need to consider turbulence.


grtabor September 19, 2008 06:52

Henrik, I came to the concl

I came to the conclusion that the term should not be there if the fluid is incompressible. Temam and Miranville "Mathematical Modelling in Continuum Mechanics" (good book, btw) gives the equation for this case as

rho C(dT/dt + u.del T) - k laplacian(T) = r

where r is the source term. Anyway, I put this and the Boussinesq term into icoFoam yesterday (isn't OF great!!) and got some rather pretty results out.


hkb September 19, 2008 08:40

Hi Gavin have you given any t
Hi Gavin
have you given any thought to the reason for the absence of the viscous heating term in the incompressible case? It seems to be fairly common to omit it in the compressible case (I can't find it in any of the standard solvers in OF I have checked), so maybe it has a neglible impact? Or is there anything in the physics that invalidates the term in the incompressible case? I guess that (with a better knowledge of c++) I could ask the solver to print the different contributions seperately, to check their relative size.


grtabor September 19, 2008 09:31

Temem and Miranville give the
Temem and Miranville give the general energy equation as

rho De/Dt + q_i,i - \sigma_ij u_i,j = r

u_i,j = div(u) so this term = 0 in incompressible flow. q is the heat transfer by conduction which is modeled by Fourier's law, and for a perfect gas, the internal energy e=CT.

Looking through the link that you provided, it includes a viscous heating term as well, \phi, but then states that it is small for almost every case, and proceeds to ignore it. They don't ignore div(u) until the very last paragraph, where they state that the Boussinesq approximation makes this zero.


hkb September 19, 2008 09:45

Hi Gavin Forgive my ignorance
Hi Gavin
Forgive my ignorance, but what is \sigma in the formula above? I actually ordered the book this morning, but I guess it won't arrive in a week or so.
Another question (also from a beginners perspective): from energy conservation, does the viscous heating equal the viscous dissipation? And in that case, could the viscous heating be estimated from the turbulent flow? Would the property of interest be epsilon?


grtabor September 19, 2008 10:16

\sigma=stress tensor (I think.
\sigma=stress tensor (I think...)

I'm guessing the viscous heating should roughly equal the viscous dissipation in turbulent flow. Or rather; any shear in the flow will dissipate energy through a viscous mechanism; your source gives the formula

2\mu e:e - 2/3 \mu del.u^2

for this (e here is the strain I think). In turbulent flow the main contribution to this is likely to be the turbulent one.


hkb September 19, 2008 12:11

Hi Gavin thanks for your pati
Hi Gavin
thanks for your patience. I have one last question:
In one of the above links I picked up this expression for the heat transfer coefficient:
This seem reasonable as this means that the heat transfer is multiplied by the relative increase of the viscosity due to turbulence:


Would you believe in this?


grtabor September 19, 2008 15:35

Hi Henrik, Yes, I can easil
Hi Henrik,

Yes, I can easilly believe that - the turbulence will enhance heat transfer. Worth checking this dimensionally though; the Prandtl number should be dimensionless, shouldn't it, so that would give [alphaEff] = [nut] - is that correct? (Sorry, end of long day and I can't be bothered to work it out properly.)


dpalko September 19, 2008 15:53

Henrik, alphaEff is not the

alphaEff is not the heat transfer coefficient, it is the heat diffusivity. If you take the law of conservation of energy, what you get after some simplifications and assumptions is div(T)=(k/cp*rho)*laplace(T), where k/(rho*cp)=alpha is thermal diffusivity. This just assumes that to the control volume the energy is transfered by conduction. The second term, nut/Pr accounts for increase of heat transport due to turbulence.


hkb September 19, 2008 16:43

I'm sorry, this is all new to
I'm sorry, this is all new to me, and it seems like I messed it up in the naming.
So taking a step backwards:
based on Fourier's law I would expect the coefficient in front of the laplacian to be k/(rho*C), at least in the laminar case. Do we agree on that? Next, it is feasible that heat will be transported faster in a turbulent case, and by adding nut/Prt we preserve dimensions, so maybe it is ok. Still with me? Finally, in the compressible cases, the coefficient in front of the laplacian is called alphaEff, however, in those cases the coefficient don't include rho and C (since they are solving for h and we are solving for T).
I tried to look into the origin of alphaEff in the compressible codes, but I didn't understand very much. Maybe someone who knows a little more about these things can add something about the similarities (or differences) between alphaEff/(rho*C) and nut/Pr+k/(C*rho)?


dpalko September 19, 2008 17:06

Maybe i should do it more clea
Maybe i should do it more clear:
When you average NS equations, you end up with one additional term, which is the Reynolds stress. This term is then modeled by turbulence model and this is when nut appears (sorry i don't know how to write equations in this forum).
The same thing as for the momentum equation, you do for the equation of transport of energy. Here you end up with one additional term. This term you then put equal to nut/Prt, where Prt is the turbulent Prandtl number, which is defined as the ratio between the turbulent transport of momentum and energy, this number is usually set up by your turbulence model and is Prt=0.9. This is the standard approach when using eddy viscosity turbulence models.


hkb October 6, 2008 11:10

Thanks a lot David so k/C*rho
Thanks a lot David
so k/C*rho takes care of conduction, while nut/Prt accounts for turbulent convection, is that the way to understand this? I looked a little further into the compressible codes, and it seems that the situation is similar, but in those cases a turbulent Prandtl number of 1.0 is used.
I realize I have to catch up on the theory in order to get further. Can you (or anyone) recommend a textbook that explains this in a complete but still comprehensible way?


christian October 16, 2008 09:30

When running the boussinesqBuo
When running the boussinesqBuoyantSimpleFoam solver found at the OpenFOAM wiki on my own case, the following error appear:

Starting time loop

Time = 1

incompatible dimensions for operation
[U[0 1 -2 0 0 0 0] ] == [-grad(p)[1 -2 -2 0 0 0 0] ]#0 Foam::error::printStack(Foam: in "/home/csvs/OpenFOAM/OpenFOAM-1.5/lib/linux64GccDPOpt/"
#1 Foam::error::abort() in "/home/csvs/OpenFOAM/OpenFOAM-1.5/lib/linux64GccDPOpt/"
#2 void Foam::checkMethod<foam::vector<double> >(Foam::fvMatrix<foam::vector<double> > const&, Foam::DimensionedField<foam::vector<double>, Foam::volMesh> const&, char const*) in "/home/csvs/OpenFOAM/OpenFOAM-1.5/applications/bin/linux64GccDPOpt/boussinesqBuo yantSimpleFoam"
#3 Foam::tmp<foam::fvmatrix<foam::vector<double> > > Foam::operator==<foam::vector<double> >(Foam::fvMatrix<foam::vector<double> > const&, Foam::tmp<foam::geometricfield<foam::vector<double >, Foam::fvPatchField, Foam::volMesh> > const&) in "/home/csvs/OpenFOAM/OpenFOAM-1.5/applications/bin/linux64GccDPOpt/boussinesqBuo yantSimpleFoam"
#4 main in "/home/csvs/OpenFOAM/OpenFOAM-1.5/applications/bin/linux64GccDPOpt/boussinesqBuo yantSimpleFoam"
#5 __libc_start_main in "/lib64/"
#6 Foam::regIOobject::readIfModified() in "/home/csvs/OpenFOAM/OpenFOAM-1.5/applications/bin/linux64GccDPOpt/boussinesqBuo yantSimpleFoam"

From function checkMethod(const fvMatrix<type>&, const GeometricField<type,>&)
in file /home/csvs/OpenFOAM/OpenFOAM-1.5/src/finiteVolume/lnInclude/fvMatrix.C at line 1208.

FOAM aborting

Anyone with an idea? I can assure you that the dimension of my U field is m/s.

gschaider October 17, 2008 09:21

Hi Christian! If you have a
Hi Christian!

If you have a closer look you might noticed that the difference in dimensions is a density. See discussions about different dimensions in pressure for compressible and incompressible solvers


All times are GMT -4. The time now is 15:22.