[Chrisi1984]

Hi Foamers,

I am simulating a turbulent compressible air flow through a curved pipe with adiabatic walls. I am using the solver rhoSimpleFoam. The air speed is quite high with a maximum of about 260m/s. The pressure and velocity fields looks well. It is similar to the outcome of Fluent.

My problem is that I have a loss in total enthalpy during the flow of about 3.8 percent. So my outlet Temperature is too low.

Viscous Heating is not imlemented in the solver, but this should be neglectible in gas fluxes, right?

So what else can be the reason for this problem?

I have read in this forum that others had similar problems, but I could not find a solution for it.


Didn't someone solved the problem ?

I am doing my master thesis with OpenFoam, and it would be really nice if you can help me.

If you need more information of my case, I can give you.

Regards
Chrisi

[bastil]

Hi,

could you please post your solver setting? If possible it would be best to have the whole case. I have now idea right now but I coult hake a look into it.

Regards Bastian

[Chrisi1984]

Ok,

very nice of you!!

Here is my case:

These are my schemes:

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
grad(U) Gauss linear;
grad(p) Gauss linear;
}

divSchemes
{
div(phi,U) Gauss upwind;
div((muEff*dev2(grad(U).T()))) Gauss linear;
div(phi,h) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,k) Gauss upwind;
div(phid,p) Gauss linear;
div(U,p) Gauss linear;

div(U) Gauss linear;//zusaetzlich fuer Energy
}

laplacianSchemes
{
laplacian(muEff,U) Gauss linear corrected;
laplacian(alphaEff,h) Gauss linear corrected;
laplacian((rho|A(U)),p) Gauss linear corrected;
laplacian((rho*rAU),p) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
laplacian((rho*(1|A(U))),p) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
div(U,p) upwind phi;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p ;
}


// ************************************************** *********************** //

my fvSolution:

solvers
{
p
{
solver GAMG;
tolerance 1e-08;
relTol 0.05;
smoother GaussSeidel;
cacheAgglomeration false;
nCellsInCoarsestLevel 10;
agglomerator faceAreaPair;
mergeLevels 1;
maxIter 500;
}

U
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-06;
relTol 0.1;
}

h
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-06;
relTol 0.1;
}

k
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.1;
}

epsilon
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.1;
}
}

SIMPLE
{
nNonOrthogonalCorrectors 0;
pMin pMin [ 1 -1 -2 0 0 0 0 ] 100;
}

relaxationFactors
{
p 0.01; // 0.3 -> 0.01
rho 0.9; // 0.05 -> 0.9
U 0.7;
k 0.7;
epsilon 0.7;
h 0.5;
}

and in the attatched file you will find my boundary and initial conditons.

Tell me if you need more information.

One more thing:
Is it possible that the high velocity magnitude is a problem for the solver?

I read that rhoSimpleFoam would only be for small Ma.
But my maximum Ma is near 1.

What can be effects of the high velocity which are neglected by rhoSimpleFoam?

Best regards
Chrisi

[bastil]

Ok this looks ok. How bad is your mesh?

Regards Bastian

[Chrisi1984]

I think my mesh is alright:

Mesh stats
points: 117072
faces: 336016
internal faces: 321704
cells: 109620
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 109620
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0

Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface topology
p 406 432 ok (non-closed singly connected)
w-out 2000 2050 ok (non-closed singly connected)
v 406 432 ok (non-closed singly connected)
w-in 2000 2050 ok (non-closed singly connected)
w-r 9500 9550 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-0.267792 -0.032787 0.0440929) (0.505012 0.0488568 0.302414)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (2.79925e-19 6.46543e-19 -4.54471e-19) OK.
Max cell openness = 2.56061e-16 OK.
Max aspect ratio = 6.1821 OK.
Minumum face area = 9.30065e-07. Maximum face area = 2.2486e-05. Face area magnitudes OK.
Min volume = 2.33102e-09. Max volume = 9.44257e-08. Total volume = 0.00176285. Cell volumes OK.
Mesh non-orthogonality Max: 56.5682 average: 8.99531
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 1.92253 OK.

Mesh OK.

Regards Chrisi

[olesen]

Originally Posted by Chrisi1984

Viscous Heating is not imlemented in the solver, but this should be neglectible in gas fluxes, right?

In general I don't think it's safe to assume that viscous heating is negligible for compressible flows, but rather than trusting an opinion, why not simply add viscous heating and see its influence for your flow?
Eg,

//- The source for the enthalpy equation resulting from
//  viscous and turbulent dissipation
tmp<volScalarField> thermalDissipationEff
(
    const compressible::turbulenceModel& turb
)
{
    tmp<volTensorField> tgradU = fvc::grad(turb.U());

    return
    (
            turb.mu()* (tgradU() && dev(twoSymm(tgradU())))
        + turb.rho() * turb.epsilon()
    );
}

For debugging purposes, you can also output this field to see how the power is distributed in your domain.

It could be, however, that the Fluent temperatures are still higher. From their docs, it appears that they may doing the same as cd-adapco and using an equilibrium condition for the thermal dissipation.

This would correspond to something like the following:

//- The source for the enthalpy equation resulting from
//  the effective viscous dissipation
//  (ie, when turbulent production and dissipation are in equilibrium)
tmp<volScalarField> thermalDissipationEquilibrium
(
    const compressible::turbulenceModel& turb
)
{
    tmp<volTensorField> tgradU = fvc::grad(turb.U());

    return
    (
        (
            turb.muEff()*dev(twoSymm(tgradU()))
          - ((2.0/3.0)*I) * turb.rho() * turb.k()
        ) && tgradU()
    );
}

I can't see a good justification to assume that the energy entering at the largest scales should appear immediately as a thermal contribution and justify an equilibrium assumption.

[Chrisi1984]

Hi,

Thank you for your help.

But I am sorry I am not able to implement this:

//- The source for the enthalpy equation resulting from
// the effective viscous dissipation
// (ie, when turbulent production and dissipation are in equilibrium)
tmp<volScalarField> thermalDissipationEquilibrium
(
const compressible::turbulenceModel& turb
)
{
tmp<volTensorField> tgradU = fvc::grad(turb.U());

return
(
(
turb.muEff()*dev(twoSymm(tgradU()))
- ((2.0/3.0)*I) * turb.rho() * turb.k()
) && tgradU()
);
}


correctly in my solver. It can not be compiled.
How exactly should the source code of hEqn look like and what must be changed in createFields?

Regards
Chrisi

[olesen]

Originally Posted by Chrisi1984

Hi,

Thank you for your help.

But I am sorry I am not able to implement this:
...
correctly in my solver. It can not be compiled.
How exactly should the source code of hEqn look like and what must be changed in createFields?

You are defining a function, so it cannot be defined within another function - if you place it in createFields.H or hEqn.H etc, you'll be trying to define a function within the 'main' function - that is not valid C/C++!!

Instead, simply place it after the OpenFOAM includes and above main().
Eg,

tmp<volScalarField> thermalDissipationEff( ... )
{
    ...
}

int main(int argc, char *argv[])
{
....
}

Within hEqn.H you can then use it

...
    if (addDissipativeHeating)
    {
         hEqn -= thermalDissipationEff(turbulence);
    }

    hEqn.relax();
...

Presumably you'll have a switch in your fvSolution::SIMPLE dictionary to turn it on/off.

[Chrisi1984]

Hi Mark,

Thank you for your help!!!

I got it running.

But the solution is not like expected.

No my outlet temerature is even lower and
The inlet pressure has increased form 105000 to 110000.

So now my result is more different compared to the Fluent solution.

But I did not change somthing in fvSolution with Simple.

What did you mean with turn on and off. Perhabs I should really change it for a better result.

Regards
Chrisi

[olesen]

Originally Posted by Chrisi1984

Hi Mark,

Thank you for your help!!!

I got it running.

But the solution is not like expected.

As I mentioned previously:

Originally Posted by olesen

For debugging purposes, you can also output this field to see how the power is distributed in your domain.

See if there is a sign error or something else, or if particular bits of the domain are doing something strange.

[johannesk]

Hi Chrisi1984,

have you further studied the problem and verified that the change to the solver is correct or if there is a sign error? I have patched rhoPorousSimpleFoam (runs much more stable for me than rhoSimpleFoam), but not yet created a simple testcase to verify everything.

I ran it in my original complex simulation though, but I can only see a minor change in the solution. Yet I do not have a basis to judge wether the solution is now better or worse.

[Chrisi1984]

Hi,

I want to give a update.

I analysed the hEqn and the equations for internal energy as well as for kinetic energy (page 68 in Jasaks Dissertation).
Because I want to make up the balance for the total enthalpy I added both.

Then I have imlemented the resulting equation for total enthalpie in the hEqn.

That is the result:

{


volScalarField u2 = 0.5 * magSqr(U);
volTensorField gradU = fvc::grad(U);
volTensorField kin1 = turbulence->muEff() * (gradU + gradU.T());
volScalarField divU = fvc::div(U);
volVectorField kin2 = 2.0/3.0 *turbulence->muEff() * divU * U;
volScalarField kin3 = fvc::div(kin2);
volVectorField kin4 = kin1 & U;
volScalarField kin5 = fvc::div(kin4);


fvScalarMatrix hEqn
(
fvm::div(phi, h)
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)

==

- kin3
+ kin5
-fvc::div(phi, u2)
);

hEqn.relax();

eqnResidual = hEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);

thermo.correct();
}


And I yet calculated with the new solver an easy case.

For me it seems that it deliver a perfect result. But the solver is not so stable.

Perhabs one should use phi instead of U to get a more conservative implementation.
It would be nice if someone have an idea how to do it, or can tell me if it is even necessary.

I would be glad about other suggestions to improve the solver.

Regards
Chrisi

[smart]

Hi Chrisi, try to change:

div(phi,U) Gauss upwind; to
div(phi,U) Gauss linear;

in the fvSchemes and let me know if you see improvement.

SM

[Chrisi1984]

Hi,

Thank you for your advice.

I yet had changed it to
div(phi,U) Gauss vanLeerV;

that gives good results.

Regards
Chris

[fredo490]

I know this post is really old (3 years) but it is still up to date.

I have some questions about the code from the post number 12. I found the original equation given by Jasaks (see enclosed) but I cannot find exactly what is included in the equation. I mean, what are the assumptions and what is considered / neglected ?

Moreover, about the code, I understand that the source term (rho*Q) is not included and that the gravity effect is neglected (rho*g*u). But there are two lines I don't understant:
1) Where does this term come from:

- fvm::Sp(fvc::div(phi), h)

2) what is:

- fvc::div(phi, u2)

3) Why is there no term about pressure in the code while the second RHS term of the equation mention the pressure ?

4) The equation seems to include the viscous effect (as it uses the effective viscousity) but when I look at my computation results, it looks like there is almost no viscous heating near the wall while other software like Fluent shows relatively important heating ( http://www.cfd-online.com/Forums/ope...mplecfoam.html )