[hartinger]

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.

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]

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]

Ignore my comment after point 2
Dave H.

[hartinger]

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]

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]

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]

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]

Hi Stefan,

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

regards
markus

[isabelle]

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]

Hi Isabelle,

yes, for gases usually negligible.

regards
markus

[sradl]

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]

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]

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]

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]

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]

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]

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]

Originally Posted by david_h

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]

Hi,

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

Thanks!