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
Henrik

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.
http://www.openfoamwiki.net/index.ph...ure_to_icoFoam

David

Hi. I would choose option 1
 

Hi.

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).

Gavin

Hi Thanks for your rapid repl
 

Hi
Thanks for your rapid reply. I have started adding a temperature transport to turbFoam (based on this thread: http://www.cfd-online.com/OpenFOAM_D....html#POST6623). 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
Henrik

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

I had a look at the code at:
 

I had a look at the code at:

http://www.openfoamwiki.net/index.ph...ure_to_icoFoam

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

Gavin

Yes, I think the viscous heati
 

Yes, I think the viscous heating is missing.
Based on http://ocw.mit.edu/NR/rdonlyres/Earth--Atmospheric--and-Planetary-Sciences/12-80 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: http://www.cfd-online.com/OpenFOAM_D...es/1/2513.html 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.

Henrik

Henrik, I came to the concl
 

Henrik,

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.

Gavin

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.

Henrik

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.

Gavin

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?

Henrik

\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.

Gavin

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:
alphaEff=nut/Pr+k/(C*rho)
This seem reasonable as this means that the heat transfer is multiplied by the relative increase of the viscosity due to turbulence:

nut/Pr+k/(C*rho)=k/(C*rho)*(mut/mu+1)=k/(C*rho)*(mut+mu)/mu

Would you believe in this?

Henrik

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.)

Gavin

Henrik, alphaEff is not the
 

Henrik,

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.

David

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)?

Henrik

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.

David

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?

Henrik

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:http://www.cfd-online.com/OpenFOAM_D...part/proud.gifstream&) in "/home/csvs/OpenFOAM/OpenFOAM-1.5/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/csvs/OpenFOAM/OpenFOAM-1.5/lib/linux64GccDPOpt/libOpenFOAM.so"
#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/libc.so.6"
#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.

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

Bernhard

Hi Bernard, I realised an h
 

Hi Bernard,

I realised an hour ago that pressure in compressible solvers is of dimension Pascal and that pressure in incompressible solvers is of dimension pressure/density. The solver is running now. Sorry, I should have posted a message that everything is well, for now http://www.cfd-online.com/OpenFOAM_D...part/happy.gif

Hi, dear foamers, When I run the boussinesqBuoyantSimpleFoam solver on my own case, the following error appears:
Time = 0.05
DILUPBiCG: Solving for Ux, Initial residual = 0.989258, Final residual = 0.0103564, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.973745, Final residual = 0.0111817, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.975119, Final residual = 0.00821149, No Iterations 1
GAMG: Solving for p, Initial residual = 9.53872e-06, Final residual = 5.94514e-07, No Iterations 1
time step continuity errors : sum local = 1.07053e+17, global = -1.94455e+13, cumulative = -1.94455e+13
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 0.00430228, No Iterations 1
bounding epsilon, min: -1.3712e+37 max: 4.32279e+38 average: 7.36798e+32
#0 _ZN4Foam5error10printStackERNS_7OstreamE-0xb34150
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libOpenFOAM.so"
#1 _ZN4Foam6sigFpe13sigFpeHandlerEi-0xad1820
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libOpenFOAM.so"
#2 Uninterpreted: [0xa0000000000107e0]
#3 _ZNK4Foam13LimitedSchemeIdNS_14vanLeerLimiterINS_6 NVDTVDEEENS_10limitFuncs6magSqrEE7limiterERKNS_14G eometricFieldIdNS_12fvPatchFieldENS_7volMeshEEE-0x1d0538e
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libfiniteVolume.so"
#4 _ZNK4Foam33limitedSurfaceInterpolationSchemeIdE7we ightsERKNS_14GeometricFieldIdNS_12fvPatchFieldENS_ 7volMeshEEE-0x1f1c680
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libfiniteVolume.so"
#5 _ZNK4Foam2fv21gaussConvectionSchemeIdE6fvmDivERKNS _14GeometricFieldIdNS_13fvsPatchFieldENS_11surface MeshEEERNS3_IdNS_12fvPatchFieldENS_7volMeshEEE-0x1f192a0
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libfiniteVolume.so"
#6 Foam::tmp<Foam::fvMatrix<double> > Foam::fvm::div<double>(Foam::GeometricField<double , Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::word const&) in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
#7 Foam::tmp<Foam::fvMatrix<double> > Foam::fvm::div<double>(Foam::GeometricField<double , Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&) in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
#8 _ZN4Foam10boussinesq9RASModels8kEpsilon7correctEv-0x28bce90
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libboussinesqRASModels.so"
#9 main in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
#10 __libc_start_main-0x734df0
in "/lib/tls/libc.so.6.1"
#11 _start in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
Floating exception
Any comments? thank you very much in advance.

Several questions arise from your output:
- why do you have a timestep!=1 for a steady solver (not that it matters much, but it gives me the suspicion that the settings you're using were meant for a transient solver)
- is this the first time-step where continuity explodes?
- do you have relaxation?
- was the vanLeer your idea?
- what are the boundary conditions?

Bernhard

Hello, Bernhard, thank you very much for your precious comments.
About the relaxation, it is shown as follows.
relaxationFactors
{
p 0.15;
U 0.3;
k 0.3;
epsilon 0.3;
R 0.7;
nuTilda 0.7;
T 0.3;
}


Sorry, I have to post the results in the follows for the reason of too many characters.

I changed the deltaT=1, and do calculation again, the results are shown in follows.

Create time

Create mesh for time = 0

Reading transportProperties


Reading environmentalProperties
Reading field p

Reading field T

Reading field Q

Reading field U

Creating field alphaEff

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
Cb 1.44;
alphaEps 0.76923;
}


Starting time loop

Convergence criterion for U = 0.001
Convergence criterion for p = 0.01
Convergence criterion for T = 0.001

Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.00959259, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.0060194, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.00859174, No Iterations 1
GAMG: Solving for p, Initial residual = 1, Final residual = 0.00880705, No Iterations 9
time step continuity errors : sum local = 0.000515452, global = -1.32363e-05, cumulative = -1.32363e-05
DILUPBiCG: Solving for epsilon, Initial residual = 0.0516045, Final residual = 0.000326048, No Iterations 1
bounding epsilon, min: -1.05947 max: 85.5521 average: 2.55034
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 0.0173683, No Iterations 1
DILUPBiCG: Solving for T, Initial residual = 0.000886871, Final residual = 2.99323e-06, No Iterations 1
ExecutionTime = 142.679 s ClockTime = 143 s

Initial residual for U = 1
Initial residual for p = 1
Initial residual for T = 0.000886871


Time = 3

DILUPBiCG: Solving for Ux, Initial residual = 0.998546, Final residual = 0.0174537, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.999007, Final residual = 0.0185219, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.999519, Final residual = 0.0171578, No Iterations 1
GAMG: Solving for p, Initial residual = 0.984913, Final residual = 0.00683274, No Iterations 3
time step continuity errors : sum local = 904905, global = -7597.79, cumulative = -7597.66
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 0.00633617, No Iterations 1
bounding epsilon, min: -2.0693e+12 max: 3.04074e+13 average: 2.12447e+08
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 0.00279179, No Iterations 1
bounding k, min: -2.59096e+10 max: 1.58311e+12 average: 4.51853e+06
DILUPBiCG: Solving for T, Initial residual = 0.0518603, Final residual = 0.00156515, No Iterations 1
ExecutionTime = 357.442 s ClockTime = 358 s

Initial residual for U = 0.999519
Initial residual for p = 0.984913
Initial residual for T = 0.0518603


Time = 5

DILUPBiCG: Solving for Ux, Initial residual = 0.989258, Final residual = 0.0103564, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.973745, Final residual = 0.0111817, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.975119, Final residual = 0.00821149, No Iterations 1
GAMG: Solving for p, Initial residual = 9.53872e-06, Final residual = 5.94514e-07, No Iterations 1
time step continuity errors : sum local = 1.07053e+19, global = -1.94455e+15, cumulative = -1.94455e+15
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 0.00430228, No Iterations 1
bounding epsilon, min: -1.3712e+37 max: 4.32279e+38 average: 7.36798e+32
#0 _ZN4Foam5error10printStackERNS_7OstreamE-0xb34150
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libOpenFOAM.so"
#1 _ZN4Foam6sigFpe13sigFpeHandlerEi-0xad1820
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libOpenFOAM.so"
#2 Uninterpreted: [0xa0000000000107e0]
#3 _ZNK4Foam13LimitedSchemeIdNS_14vanLeerLimiterINS_6 NVDTVDEEENS_10limitFuncs6magSqrEE7limiterERKNS_14G eometricFieldIdNS_12fvPatchFieldENS_7volMeshEEE-0x1d0538e
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libfiniteVolume.so"
#4 _ZNK4Foam33limitedSurfaceInterpolationSchemeIdE7we ightsERKNS_14GeometricFieldIdNS_12fvPatchFieldENS_ 7volMeshEEE-0x1f1c680
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libfiniteVolume.so"
#5 _ZNK4Foam2fv21gaussConvectionSchemeIdE6fvmDivERKNS _14GeometricFieldIdNS_13fvsPatchFieldENS_11surface MeshEEERNS3_IdNS_12fvPatchFieldENS_7volMeshEEE-0x1f192a0
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libfiniteVolume.so"
#6 Foam::tmp<Foam::fvMatrix<double> > Foam::fvm::div<double>(Foam::GeometricField<double , Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::word const&) in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
#7 Foam::tmp<Foam::fvMatrix<double> > Foam::fvm::div<double>(Foam::GeometricField<double , Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&) in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
#8 _ZN4Foam10boussinesq9RASModels8kEpsilon7correctEv-0x28bce90
in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/lib/linuxIA64GccDPOpt/libboussinesqRASModels.so"
#9 main in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
#10 __libc_start_main-0x734df0
in "/lib/tls/libc.so.6.1"
#11 _start in "/home/g2/e090012/OpenFOAM/OpenFOAM-1.5/applications/bin/linuxIA64GccDPOpt/boussinesqBuoyantSimpleFoam"
Floating exception

No Idea what could be the problem. Two recommendations:

- when starting a new case base the settings on a stable calculation with the same solver. Usually such settings can be found in $FOAM_TUTORIALS. So what I'd recommend is to use fvSchemes, fvSolution and controlDict from the buoyantSimple(!)Foam-hotroom as a first guess for your case. The references to vanLeer indicate that this is not the case for your calculation
- for outflows that have zeroGradient for U don't use zeroGradient for transported quantities like k, epsilon and T. Use inletOutlet. Otherwise things might explode when you get a backflow there

Bernhard