CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Solving the total energy equation

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree11Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   January 26, 2010, 13:09
Default Solving the total energy equation
  #1
Member
 
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 40
Rep Power: 8
dohnie is on a distinguished road
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
mm.abdollahzadeh likes this.
dohnie is offline   Reply With Quote

Old   October 5, 2010, 10:30
Default
  #2
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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
cliffoi is offline   Reply With Quote

Old   October 5, 2010, 10:40
Default
  #3
Member
 
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 40
Rep Power: 8
dohnie is on a distinguished road
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.
dohnie is offline   Reply With Quote

Old   October 5, 2010, 10:44
Default
  #4
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   October 5, 2010, 11:27
Default
  #5
Member
 
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 40
Rep Power: 8
dohnie is on a distinguished road
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.
dohnie is offline   Reply With Quote

Old   October 5, 2010, 13:44
Default
  #6
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
Thanks Florian, this will make a very good starting point for my work.

Ivor
cliffoi is offline   Reply With Quote

Old   October 11, 2010, 04:37
Default
  #7
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 198
Rep Power: 7
Chris Lucas is on a distinguished road
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
Chris Lucas is offline   Reply With Quote

Old   October 12, 2010, 17:36
Default
  #8
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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 likes this.
cliffoi is offline   Reply With Quote

Old   October 14, 2010, 13:46
Default
  #9
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 198
Rep Power: 7
Chris Lucas is on a distinguished road
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
File Type: h hEqn.H (565 Bytes, 96 views)
Chris Lucas is offline   Reply With Quote

Old   October 14, 2010, 15:29
Default
  #10
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   October 18, 2010, 03:23
Default
  #11
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 198
Rep Power: 7
Chris Lucas is on a distinguished road
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
Chris Lucas is offline   Reply With Quote

Old   October 21, 2010, 04:13
Default
  #12
Member
 
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 40
Rep Power: 8
dohnie is on a distinguished road
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
dohnie is offline   Reply With Quote

Old   November 8, 2010, 08:58
Default
  #13
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 198
Rep Power: 7
Chris Lucas is on a distinguished road
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
Chris Lucas is offline   Reply With Quote

Old   November 9, 2010, 01:53
Default
  #14
Senior Member
 
Nakul
Join Date: Apr 2010
Location: India
Posts: 147
Rep Power: 7
nakul is on a distinguished road
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
nakul is offline   Reply With Quote

Old   November 11, 2010, 17:19
Default
  #15
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   November 16, 2010, 02:26
Question
  #16
New Member
 
Hait Anatoly
Join Date: Sep 2010
Posts: 3
Rep Power: 6
HaitAV is on a distinguished road
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!
HaitAV is offline   Reply With Quote

Old   November 23, 2010, 05:50
Default
  #17
New Member
 
Hait Anatoly
Join Date: Sep 2010
Posts: 3
Rep Power: 6
HaitAV is on a distinguished road
Hi

Can anybody help me?

All my attempts to create good solver based on total enthalpy for supersonic flows were unsuccessful.

Thanks again.
HaitAV is offline   Reply With Quote

Old   December 11, 2012, 10:30
Default
  #18
New Member
 
Traib
Join Date: Sep 2012
Posts: 27
Rep Power: 4
Traib is on a distinguished road
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.
Traib is offline   Reply With Quote

Old   December 11, 2012, 13:34
Default
  #19
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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()))
    );
}
Note that divDevRhoReff is nothing more than the implicit matrix form of div(devRhoReff), i.e. devRhoReff is the stress tensor.

Regards
Ivor
cliffoi is offline   Reply With Quote

Old   December 11, 2012, 14:46
Default
  #20
New Member
 
Traib
Join Date: Sep 2012
Posts: 27
Rep Power: 4
Traib is on a distinguished road
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
Traib is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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 Pre-Processing 9 September 9, 2008 12:53
Why FVM for high-Re flows? Zhong Lei Main CFD Forum 23 May 14, 1999 13:22


All times are GMT -4. The time now is 05:29.