# Solving the total energy equation

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 LinkBack Thread Tools Display Modes
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: 9
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:
 volScalarField eTot ( IOobject ( "eTot", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), e+0.5*magSqr(U), e.boundaryField().types() );
and solve a transport equation for eTot.

Can anybody tell me how I have to reformulate the energy transport equation
Quote:
 { fvScalarMatrix eEqn ( fvm::ddt(rho, e) + mvConvection->fvmDiv(phi, e) == DpDt ); eEqn.relax(); eEqn.solve(); thermo.correct(); }
?

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: 9 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: 9 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();``` dongchao yang, mm.abdollahzadeh and babala like this. 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: 9 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: 9 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()) );``` (here it is actually solved for rhoE, but you might do this for eTot as well 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())``` Hope that works for you! mm.abdollahzadeh and Pagoda like this.

 October 5, 2010, 13:44 #6 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 9 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: 9 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. mm.abdollahzadeh and elham usefi like this.

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
Attached Files
 hEqn.H (565 Bytes, 106 views)

 October 14, 2010, 15:29 #10 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 9 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 non-standard 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. mm.abdollahzadeh likes this.

 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: 9 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)); }``` where rhoE corresponds to rho*(e+K) and the div(p*U) is included in the phiEp. Can anybody explain why sigmaDotU appears in the inviscid part? you have to use the difference between the fvm and fvc operator for the laplacian? Thanks, 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: energy equation in rhoCentralFoam

 November 11, 2010, 17:19 #15 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 9 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. mm.abdollahzadeh likes this.

 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 Ranque-Hilsch 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 high-speed 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. charmc likes this.

 December 11, 2012, 13:34 #19 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 9 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 kEpsilon::devRhoReff() const { return tmp ( new volSymmTensorField ( IOobject ( "devRhoReff", runTime_.timeName(), mesh_, IOobject::NO_READ, IOobject::NO_WRITE ), -muEff()*dev(twoSymm(fvc::grad(U_))) ) ); } tmp kEpsilon::divDevRhoReff(volVectorField& U) const { return ( - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T())) ); }``` Note that divDevRhoReff is nothing more than the implicit matrix form of div(devRhoReff), i.e. devRhoReff is the stress tensor. 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 Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post xiuying OpenFOAM Running, Solving & CFD 8 August 27, 2013 15:33 vw.cfd OpenFOAM 6 August 7, 2009 05:44 carsten OpenFOAM Bugs 11 September 12, 2008 11:16 sivakumar OpenFOAM Pre-Processing 9 September 9, 2008 12:53 Zhong Lei Main CFD Forum 23 May 14, 1999 13:22

All times are GMT -4. The time now is 20:40.

 Contact Us - CFD Online - Top