CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Full energy equation for enthalpy (http://www.cfd-online.com/Forums/openfoam-solving/59468-full-energy-equation-enthalpy.html)

hartinger June 8, 2006 06:22

Hi, i've implemented the fu
 
Hi,

i've implemented the full energy equation for enthalpy without neglecting anything and put it into sonicLiquidFoam.

I would like to post it here to check whether i got the physics right, and to discuss whether the discretisation is correct and appropriate.
I've tested it and the results look alright. Code and testcase is attached.
I'm sure it can be improved, any comments welcome.

http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif testCaseAndCode.tar

So, here we go:

energy equation:
rho Dh/Dt = - grad(k div(T)) - tau && grad(U) + DpDt

enthalpy:
h = u + p/rho // u - internal energy
// enthalpy - temperature relation for liquids
h = Cp(T-T0) // T0 - ambient temperature | Cp - heat capacity

substantial derivatives:
('d' being partial differential)
Dh/Dt = dh/dt + U & grad(h)

Dp/Dt = dp/dt + U & grad(p)
with U & grad(p) = div(pU) - p div(U)
Dp/Dt = dp/dt + div(pU) - p div(U)

conduction term:
grad(k div(T)) = grad(k/Cp div(h))

that translated into foamSpeak:

volTensorField gradU = fvc::grad(U);

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

volScalarField tauGradU = tau && gradU;

solve
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(kl/Cp, h)
==
fvc::ddt(p)
+ fvc::div(p*U)
- p * fvc::div(U)
- tauGradU
);

T = T0 + h / Cp;

thank you
markus

david_h June 8, 2006 09:03

Markus, I see 3 things whic
 
Markus,

I see 3 things which you might look at:

1) Is the enthalpy of formation included in the definition of enthalpy(i.e.):
h = h_0 + cp*(T - T0)
Certain thermodynamic classes include this in
the definition of enthalpy.

2) There is a function to fvc:DDt() to calculate
the equivalent of:
fvc::ddt(p)
+ fvc::div(p*U)
- p * fvc::div(U)
(look at pEqn.H in rhoTurbFoam)
should the last term be "-U & grad(p)" ?

3) You might also try discretizing tauGradU as:
div( phi/interplate(rho) & interpolate(tau) ) - U & div(tau)
(I believe this is similar to procedure used for calculating the U & grad(p) term in the DDt() )

David H.

david_h June 8, 2006 09:14

Ignore my comment after point
 
Ignore my comment after point 2
Dave H.

hartinger June 8, 2006 12:35

David, thanks for your remarks
 
David, thanks for your remarks.

to 1)
personally i won't do any reaction, so should be alright to assume the enthalpy of formation to be zero.

to 2)
i wasn't aware of that function, thanks

> should the last term be "-U & grad(p)" ?
what do you mean? i replaced "U & grad(p)" with "div(pU) - p div(U)" to have a more conservative form. that's correct, isn't?

thank you
Markus

david_h June 12, 2006 10:20

Marcus, Below I derived an
 
Marcus,

Below I derived an internal energy equation in terms of pressure and internal enthalpy by subtracting the kinetic energy equation from the total energy equation (Kuo,Combustion). I get a "non-conservative" form for the pressure:
ddt(p) + U & grad(p)

and a "conservative" form for the stress:
div( phi/rho & tau ) - U & div(tau)

derivation:

total energy equation:
ddt(rho,h_T) + div(phi,h_T) - ddt(p) =
div(-q) + div(phi/rho & tau )

h_T = h + K
K ~ mean kinetic energy
h_T ~ total enthalpy

1) ddt(rho,h) + div(phi,h) - ddt(p) -
ddt(rho,K) + div(phi,K) = div(-q) + div(phi/rho & tau )

kinetic energy equation
ddt(rho,K) = U & ddt(rhoU) + K ddt(rho)

2) ddt(rho,K) + U & div(phi,U) + U & grad(p)
+ K & div(phi) = U & div(tau)

subtracting 2) from 1) give an internal energy:

3) ddt(rho,h) + div(phi,h)
- ( ddt(p) + U & grad(p) ) = div(-q) +
( div( phi/rho & tau) - U & div(tau) )

Dave

hartinger June 13, 2006 08:46

yep, looks good. I guess it i
 
yep, looks good.
I guess it is a good idea to replace 'tau && grad(U)' with 'div( phi/rho & tau ) - U & div(tau)' as you suggested

markus

sradl July 17, 2007 07:19

I've got a question regarding
 
I've got a question regarding the flux calculation in Markus' programm (which is the same as in the sonicLiquidFoam solver).

Obviously phi is the surface mass flux. in the sonicLiquidEnergyFoam however, a new surface flux is introduced: phid. This should be something like a volume flux as it is multiplied with a pressure - howevever only something like a volume...

So can anybode explain what does is the meaning of phid in the code section:

surfaceScalarField phid = psi*
(
(fvc::interpolate(rho*U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)/fvc::interpolate(rho);

phi = (rho0/psi - p0)*phid;

cheers
Stefan

hartinger July 17, 2007 08:45

Hi Stefan, have a look here
 
Hi Stefan,

have a look here, might help
http://www.cfd-online.com/OpenFOAM_D...ges/1/104.html

regards
markus

isabelle July 31, 2007 09:55

Hi, I have a question about
 
Hi,

I have a question about the enthalpy equation implemented in various solvers of openFOAM-1.4,
such as buoyantFoam, sonicTurbFoam or dieselEngineFoam for instance. If I understand correctly, it seems that the viscous heating term "tauGradU", that was in Markus test (1st message above) is not accounted for in OpenFoam-1.4. Is it because viscous heating can often be assumed negligible or is there any other reason ?

Best regards,
Isabelle

hartinger July 31, 2007 10:09

Hi Isabelle, yes, for gases
 
Hi Isabelle,

yes, for gases usually negligible.

regards
markus

sradl July 31, 2007 10:21

Dear Foamers, I've combined
 
Dear Foamers,

I've combined Markus' code with a new Non-Newtonian models and Heat Transfer in the surrounding solid. If someone is interested for testing, I can send him the package via email (however, it would be a quite big package of libraries, solvers and boundary conditions).

The heat of dissipation in my problem (polymer flow) plays a significant role: the temperature profile has quite significant gradients with DT = 7 [K].

br
Stefan Radl

cosimobianchini August 28, 2007 12:32

Hi David, I would like an ex
 
Hi David,
I would like an explanation about something you wrote sometime ago in this thread:

tau && grad(U) = div( phi/interplate(rho) & interpolate(tau) ) - U & div(tau);

I would have written something like:
tau && grad(U) = div( U & tau ) - U & div(tau)

In other terms how can you do (scalar & tensor)? what does that mean? (I tried something like that but I didn't manage even to compile it)

I'm asking you this because I'm implementing a solver for total enthalpy so I directly need to model the div( U & tau ) term. I would prefer to make it dependent on mass flux instead of velocity. How can I do it?
Thanks
Cosimo

david_h August 30, 2007 09:20

Cosimo, I made a mistake in
 
Cosimo,

I made a mistake in writing the first term of the equation. I used phi as though it represented the "mass flow rate vector" at the face center instead of the "mass flow rate" (scalar) at the face center normal to the face. The second term in the equation comes from subtracting the relevant term from momentum equation dotted with the velocity vector.

I would revise the term as:

uFace = interpolate(U,phi,"methodString");
uFace = uFace + ( phi/interpolate(rho) - uFace & mesh.Sf() ) * mesh.Sf()/mesh.magSf()/mesh.magSf();

viscDiss = div( uFace & interpolate(tau) & mesh.Sf() ) - U & div(tau)

The reasoning of the above was to derive the viscous dissipation term using the discretized forms of the "total energy" and "momentum" equations as opposed to discretizing the viscous dissipation term in the "internal energy" equation from its continuous representation (tau && grad(U) ).

The method as outlined above implies an inconsistency in how the stress is calculated in the momentum and energy equations. The stress used in the momentum equation is calculated directly (nearly) at the face. Whereas the above expression interpolates the stress at the face from the stress evaluated at the cell centers.

Dave

cosimobianchini August 31, 2007 10:45

Hi David, thanks for the answe
 
Hi David, thanks for the answer.

Just a couple of questions more:

1) Markus said viscous heating = tau && gradU is usually negligible for gases.
However I'm solving for total enthalpy where viscous dissipation = div(U & tau) = tau && gradU + U & div(tau). Is the error for neglecting such term of the same order of magnitude as using the static enthalpy formulation or not?

2) I'm not sure wether or not turbulent stress R should be added to tau in case of turbulent flows?

Thanks a lot,
Cosimo

hartinger September 1, 2007 05:40

Hey Cosimo, 1) whether visc
 
Hey Cosimo,

1) whether viscous heating is negligible is determined by the physics. A dimensional analysis should tell you. The way you implement it doesn't really matter.

2) look at sonicTurbFoam, the turbulent fluxes need to be modeled, which is done there using the effective turbulent thermal diffusivity. If you consider viscous heating you should use the effective turbulent viscosity.

regards
markus

cosimobianchini September 3, 2007 08:21

Hi Markus and thanks for the r
 
Hi Markus and thanks for the reply but I need further clarification.

1) Maybe my question was not very clear. I was asking if on the same case (same physics) the error for neglecting viscous work was bigger in case of static enthalpy equation or total enthalpy equation.
The source terms for these two different conservation equations are not the same:
staticEnthalpy = tau && gradU
totalEnthalpy = div(U & tau) = tau && gradU + U & div(tau)
Is there a general rule to establish which one is bigger?

2) Thanks you confirm what I have implemented was correct.

regards
Cosimo

mehulkumar April 22, 2009 06:24

nonNewtonianFoam with energy eqn
 
Hey..
I have a case wich require nonNewtonianFoam solver with energy equation.
As mentioned here..I am intrested to test this solver.
Can any body provide me code for this???

Thanks
Mehulkumar
Email: patel.mehulkumar.l@gmail.com

Chrisi1984 May 21, 2010 07:43

how to implement this in OpenFoam
 
Quote:

Originally Posted by david_h (Post 197236)
Cosimo,

I made a mistake in writing the first term of the equation. I used phi as though it represented the "mass flow rate vector" at the face center instead of the "mass flow rate" (scalar) at the face center normal to the face. The second term in the equation comes from subtracting the relevant term from momentum equation dotted with the velocity vector.

I would revise the term as:

uFace = interpolate(U,phi,"methodString");
uFace = uFace + ( phi/interpolate(rho) - uFace & mesh.Sf() ) * mesh.Sf()/mesh.magSf()/mesh.magSf();

viscDiss = div( uFace & interpolate(tau) & mesh.Sf() ) - U & div(tau)

The reasoning of the above was to derive the viscous dissipation term using the discretized forms of the "total energy" and "momentum" equations as opposed to discretizing the viscous dissipation term in the "internal energy" equation from its continuous representation (tau && grad(U) ).

The method as outlined above implies an inconsistency in how the stress is calculated in the momentum and energy equations. The stress used in the momentum equation is calculated directly (nearly) at the face. Whereas the above expression interpolates the stress at the face from the stress evaluated at the cell centers.

Dave

Hi,

I now that this read is old, but i would like to know how I should implement this

uFace = interpolate(U,phi,"methodString");
uFace = uFace + ( phi/interpolate(rho) - uFace & mesh.Sf() ) * mesh.Sf()/mesh.magSf()/mesh.magSf();

viscDiss = div( uFace & interpolate(tau) & mesh.Sf() ) - U & div(tau)

in OpenFoam (1.6) speech?

I didn't succed.

Regards Chrisi

SiCaRiUs April 23, 2012 07:31

Hi,

is it possible to integrate these terms into the simpleFoam solver to get the heat due to friction in an moving liquid?

Thanks!


All times are GMT -4. The time now is 00:45.