[dohnie]

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

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

{
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

[cliffoi]

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

[dohnie]

Hi Ivor,
yes, I got it worked out, it seems to run fine:

    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();

[cliffoi]

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.

[dohnie]

Ivor,
I defined new boundary conditions for eTot which was a little tricky. A more elegant way is shown in rhoCentralFoam:

        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

   // 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!

[cliffoi]

Thanks Florian, this will make a very good starting point for my work.

Ivor

[Chris Lucas]

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

[cliffoi]

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.

[Chris Lucas]

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

[cliffoi]

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.

[Chris Lucas]

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

[dohnie]

In the density based solver rhoCentralFoam the inviscid part is solved first, the (optional) viscous terms are added afterwards:

        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

[Chris Lucas]

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

[nakul]

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.cfd-online.com/Forums/ope...ntralfoam.html

[cliffoi]

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.

[HaitAV]

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]

Hi

Can anybody help me?

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

Thanks again.

[Traib]

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.

[cliffoi]

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

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

[Traib]

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