
[Sponsors] 
January 26, 2010, 13:09 
Solving the total energy equation

#1  
Member
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 41
Rep Power: 10 
Dear FOAMers,
I was wondering why the compressible solvers in OF solve a transport equation for static energy (or enthalpy)? My idea was to keep the internal energy e but solve for the total energy. So I created a field Quote:
Can anybody tell me how I have to reformulate the energy transport equation Quote:
How do I make sure that the kinetic part of eTot is reevaluated after every time step? Your help is greatly appreciated. Florian 

October 5, 2010, 10:30 

#2 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
Hi Florian,
Did you have any success implementing the total energy equation? This is something I need to look into over the next few months. Ivor 

October 5, 2010, 10:40 

#3 
Member
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 41
Rep Power: 10 
Hi Ivor,
yes, I got it worked out, it seems to run fine: Code:
fvScalarMatrix eTotEqn ( fvm::ddt(rho, eTot) + mvConvection>fvmDiv(phi, eTot)  fvm::laplacian(turbulence>alphaEff(), eTot) + fvc::div(p*U) ); eTotEqn.relax(); eTotEqn.solve(); e = eTot  0.5*magSqr(U); thermo.correct(); Last edited by dohnie; October 5, 2010 at 11:28. 

October 5, 2010, 10:44 

#4 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
Florian,
It looks right. Thanks for sharing it. Did you run into any problems with boundary conditions? I've noticed issues in the past when copying boundary types from one field to another. 

October 5, 2010, 11:27 

#5 
Member
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 41
Rep Power: 10 
Ivor,
I defined new boundary conditions for eTot which was a little tricky. A more elegant way is shown in rhoCentralFoam: Code:
e = rhoE/rho  0.5*magSqr(U); e.correctBoundaryConditions(); thermo.correct(); rhoE.boundaryField() = rho.boundaryField()* ( e.boundaryField() + 0.5*magSqr(U.boundaryField()) ); Code:
// solve the eTot equation as shown above e = eTot  0.5*magSqr(U); e.correctBoundaryConditions(); thermo.correct(); eTot.boundaryField() = e.boundaryField() + 0.5*magSqr(U.boundaryField()) 

October 5, 2010, 13:44 

#6 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
Thanks Florian, this will make a very good starting point for my work.
Ivor 

October 11, 2010, 04:37 

#7 
Senior Member
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 200
Rep Power: 10 
Hi,
I'm not sure, but I think there might be an error in the equation above. fvm::laplacian(turbulence>alphaEff(), eTot) This part of the equation is the heat transfer term (or in this case the internal energy diffusion) plus the viscous work term due to viscous stress. However, I think that the viscous work term should be multiplied with the effective viscosity instead of the effective alpha. Please correct me if I'm wrong! Regards, Christian 

October 12, 2010, 17:36 

#8 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
Yes, I also noticed there is something not quite right with the conduction term. I have looked at a few texts on the subject and in all cases the total energy equation is written with the laplacian term as a function of temperature, i.e. laplacian(k,T), so they aren't much help. If you follow the definition for Etot=e+0.5*(U & U), the derivative becomes dE=de+(U & dU).
When applied to the heat conduction term the resulting expression is fvm::laplacian(alpha, e)  fvc::div(alpha*(U & grad(U))) I haven't looked at the exact discretisation of the second term but there are a few options here. fvc::div((fvc::interpolate(alpha*U) & fvc::snGrad(U)) * mesh.magSf()) fvc::div(fvc::interpolate(alpha*((U & fvc::grad(U)) & mesh.Sf())) At first glance, however, this looks like it could be a source of significant error if treated incorrectly. I'm looking to do conjugate heat transfer calculations and this so far suggests to me that I should maybe forget about trying to use total energy. 

October 14, 2010, 13:46 

#9 
Senior Member
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 200
Rep Power: 10 
Hi,
I think that the total enthalpy equation should be like this volScalarField K=0.5*magSqr(U); // we want to solve d/dt{rho*(h + K)} +div{phi*(h+K)}= d/dt{p} + div{ alpha_eff * grad{h}} + div{ U&tau} // tau > stress tensor NS equation // so tau=mu_eff*(grad(U)+grad(U).T()2/3*I*trace(grad(U)) { solve ( fvm::ddt(rho, h) + fvm::div(phi, h)  fvm::laplacian(turbulence>alphaEff(), h) == fvc::div(phi, K) fvc::ddt(rho, K) +fvc::ddt(p) +fvc::div(turbulence>muEff()*(fvc::grad(U)&U)) +fvc::div(turbulence>muEff()*(dev2(fvc::grad(U)().T())&U)) ); thermo.correct(); } I haven't validated the code, but I don't get a compiler error and the code runs in a test case. If anyone as any suggestions or sees any error, please tell me Regards, Christian 

October 14, 2010, 15:29 

#10 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
I agree with Christian's approach. This bypasses the issues I pointed out in my last post and removes the need to try construct boundaries for the total energy based on temperature and/or internal energy (This can be problematic for nonstandard boundaries). The discretization seems correct as well. I also noticed that the shear stress tensor can be replaced with turbulence>devRhoReff() to make sure it is consistent with the turbulence model we're using.
Using the same approach, the total energy equation becomes. solve ( fvm::ddt(rho, e) + fvm::div(phi, e)  fvm::laplacian(turbulence>alphaEff(), e) ==  fvc::div(phi, K)  fvc::ddt(rho, K) + fvc::div(p*U)  fvc::div(turbulence>devRhoReff() & U) ); I have a question about the diffusion term... sonicFoam writes it using turbulence>alphaEff() as above. Should there be some form of gamma=Cp/Cv correction to this? fvm::laplacian(gamma*turbulence>alphaEff(), e) This may not be correct for turbulent flows though. 

October 18, 2010, 03:23 

#11 
Senior Member
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 200
Rep Power: 10 
Hi,
I have to agree with you that the gamma is missing for the laminar heat transfer. However, I'm not sure if this is true for the turbulent part. Regards, Christian 

October 21, 2010, 04:13 

#12 
Member
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 41
Rep Power: 10 
In the density based solver rhoCentralFoam the inviscid part is solved first, the (optional) viscous terms are added afterwards:
Code:
solve ( fvm::ddt(rhoE) + fvc::div(phiEp)  fvc::div(sigmaDotU) ); e = rhoE/rho  0.5*magSqr(U); e.correctBoundaryConditions(); thermo.correct(); rhoE.boundaryField() = rho.boundaryField()* ( e.boundaryField() + 0.5*magSqr(U.boundaryField()) ); if (!inviscid) { volScalarField k("k", thermo.Cp()*mu/Pr); solve ( fvm::ddt(rho, e)  fvc::ddt(rho, e)  fvm::laplacian(thermo.alpha(), e) + fvc::laplacian(thermo.alpha(), e)  fvc::laplacian(k, T) ); thermo.correct(); rhoE = rho*(e + 0.5*magSqr(U)); } and the div(p*U) is included in the phiEp. Can anybody explain why
Florian 

November 8, 2010, 08:58 

#13 
Senior Member
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 200
Rep Power: 10 
Hi,
have you tested the total energy equation yet. In my simulation, I get a very unphysical result for my temperature? The temperature has sudden jumps. Also, the temperature given at the inlet suddenly changes to another value. Have you encountered similar problems? Regards, Christian 

November 9, 2010, 01:53 

#14 
Senior Member
Nakul
Join Date: Apr 2010
Location: India
Posts: 147
Rep Power: 8 
Hi Florian,
Did you find any answer to the questions you posted in your last post ? I am also confused about this formulation of rhoCentralFoam!! I If you have found any solution, please do share!! My specific queries are posted here: http://www.cfdonline.com/Forums/ope...ntralfoam.html 

November 11, 2010, 17:19 

#15 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
Hi all,
I have spent some time looking at this and finally implemented the total enthalpy equation as follows: volScalarField Ek = 0.5*magSqr(U); htot = h + Ek; htot.correctBoundaryConditions(); { solve ( fvm::ddt(rho, htot) + fvm::div(phi, htot)  fvm::laplacian(turbulence>alphaEff(), htot) + fvc::laplacian(turbulence>alphaEff(), Ek) == fvc::ddt(p) + fvc::div(tau & U, "div((muEff*dev(twoSymm(fvc::grad(U))))&U)") ); } h = htot  Ek; h.correctBoundaryConditions(); Here tau is the stress tensor. Even though it's small I have included it in my solver. Note that the total enthalpy requires customized boundaries, like the gradientEnthalpyFvPatchScalarField and similar classes supplied in thermophysicalModels to get consistent solutions. The solutions from this implementation are looking pretty good. Going back to my previous comment, if we use enthalpy we don't have to worry about whether turbulence>alphaEff() is applicable. 

November 16, 2010, 02:26 

#16 
New Member
Hait Anatoly
Join Date: Sep 2010
Posts: 3
Rep Power: 8 
Hi all
I am a postgraduate student and I research air flow in RanqueHilsch vortex tubes. I can not make adequate model for energy separation in OpenFOAM, though ANSYS CFX gives not bad results. If I understand you correct, this problem may be in Energy Equation. That is why I try to correct my sonicFoam solver in agreement with your proposals, but I still have wrong solution. In some cases temperature in vortex chamber decrease lower 80 K. This is incorrect result of course. Cliffoi, can you explain me how to use your algorithm? 1) How to create htot field. 2) How to declare tau tensor. etc. I am not good C++ programmer. If it is possible, can you give me total source code of your solver? May be my problem are not in energy equation? Thank you very much! 

November 23, 2010, 05:50 

#17 
New Member
Hait Anatoly
Join Date: Sep 2010
Posts: 3
Rep Power: 8 
Hi
Can anybody help me? All my attempts to create good solver based on total enthalpy for supersonic flows were unsuccessful. Thanks again. 

December 11, 2012, 10:30 

#18 
New Member
Traib
Join Date: Sep 2012
Posts: 27
Rep Power: 6 
Hi,
I have written the hsEqn for highspeed reactive flow, including thermal conduction and diffusion, as follows. volScalarField k("k", thermo.Cp()*turbulence>muEff()/Pr); // thermal conductivity { fvScalarMatrix hsEqn ( fvm::ddt(rho, hs) + mvConvection>fvmDiv(phi, hs)  fvm::laplacian(turbulence>alphaEff(), hs) // thermal Diffusion  fvc::laplacian(k,T) // thermal Conduction == dpdt + fvc::div((turbulence>muEff()*(fvc::grad(U)) + turbulence>muEff()*dev2(Foam::T(fvc::grad(U))))&U) // rate of work done by stresses  (fvc::ddt(rho, K) + fvc::div(phi, K)) // kinetic energy terms + combustion>Sh() // Source for enthalpy equation in reactive flows ); hsEqn.relax(); hsEqn.solve(); thermo.correct(); } The solver compiles fine and I checked with a test case and it looks good. Unlike some earlier posts here, I think the stress term tau is actually the whole expression: turbulence>muEff()*(fvc::grad(U)) + turbulence>muEff()*dev2(Foam::T(fvc::grad(U)))) , rather than just the second term (devrhoReff(U)). I did the matrix manipulation with hand and the results confirm this expression. However, I found the second term alone is defined in rhoCentralFoam as tauMC, which is later combined with the earlier term in UEqn and energy equation (using sigmaDotU function) Please correct me if I have missed something. Regards. 

December 11, 2012, 13:34 

#19 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
Hi Traib,
The stress tensor is defined in each turbulence model by the member function devRhoReff. There should be no need for you to manually discretize this. e.g. from kEpsilon.C Code:
tmp<volSymmTensorField> kEpsilon::devRhoReff() const { return tmp<volSymmTensorField> ( new volSymmTensorField ( IOobject ( "devRhoReff", runTime_.timeName(), mesh_, IOobject::NO_READ, IOobject::NO_WRITE ), muEff()*dev(twoSymm(fvc::grad(U_))) ) ); } tmp<fvVectorMatrix> kEpsilon::divDevRhoReff(volVectorField& U) const { return (  fvm::laplacian(muEff(), U)  fvc::div(muEff()*dev2(fvc::grad(U)().T())) ); } Regards Ivor 

December 11, 2012, 14:46 

#20 
New Member
Traib
Join Date: Sep 2012
Posts: 27
Rep Power: 6 
Hi Ivor,
Actually it was from the keModel codes that I had found the definition of the function devDivRhoReff() but I did not know it was resolved within the turbulence equations I don't fully understand the function calls inside these codes due to my limited programming knowledge. I was about to use the code with my new case but now I will change it first, and I should pay more attention to how the codes are implemented. Thanks a lot for letting me know about this. Traib 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
How to write k and epsilon before the abnormal end  xiuying  OpenFOAM Running, Solving & CFD  8  August 27, 2013 15:33 
Error log  vw.cfd  OpenFOAM  6  August 7, 2009 05:44 
Differences between serial and parallel runs  carsten  OpenFOAM Bugs  11  September 12, 2008 11:16 
Unknown error  sivakumar  OpenFOAM PreProcessing  9  September 9, 2008 12:53 
Why FVM for highRe flows?  Zhong Lei  Main CFD Forum  23  May 14, 1999 13:22 