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

rhoSimpleFoam loss in total enthalpy

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

Like Tree5Likes
  • 1 Post By olesen
  • 2 Post By olesen
  • 2 Post By Chrisi1984

Reply
 
LinkBack Thread Tools Display Modes
Old   May 23, 2010, 07:37
Default rhoSimpleFoam loss in total enthalpy
  #1
Senior Member
 
Join Date: Jan 2010
Location: Stuttgart
Posts: 129
Rep Power: 7
Chrisi1984 is on a distinguished road
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
Chrisi1984 is offline   Reply With Quote

Old   May 23, 2010, 12:03
Default
  #2
Senior Member
 
BastiL
Join Date: Mar 2009
Posts: 471
Rep Power: 11
bastil is on a distinguished road
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
bastil is offline   Reply With Quote

Old   May 23, 2010, 13:15
Default
  #3
Senior Member
 
Join Date: Jan 2010
Location: Stuttgart
Posts: 129
Rep Power: 7
Chrisi1984 is on a distinguished road
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
Attached Files
File Type: gz 0.tar.gz (1.4 KB, 18 views)
Chrisi1984 is offline   Reply With Quote

Old   May 23, 2010, 17:01
Default
  #4
Senior Member
 
BastiL
Join Date: Mar 2009
Posts: 471
Rep Power: 11
bastil is on a distinguished road
Ok this looks ok. How bad is your mesh?

Regards Bastian
bastil is offline   Reply With Quote

Old   May 23, 2010, 18:49
Default
  #5
Senior Member
 
Join Date: Jan 2010
Location: Stuttgart
Posts: 129
Rep Power: 7
Chrisi1984 is on a distinguished road
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
Chrisi1984 is offline   Reply With Quote

Old   May 31, 2010, 05:12
Default
  #6
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18
olesen will become famous soon enough
Quote:
Originally Posted by Chrisi1984 View Post
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,
Code:
//- 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:
Code:
//- 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.
mm.abdollahzadeh likes this.
olesen is offline   Reply With Quote

Old   June 2, 2010, 03:23
Default how to implement this in OpenFoam
  #7
Senior Member
 
Join Date: Jan 2010
Location: Stuttgart
Posts: 129
Rep Power: 7
Chrisi1984 is on a distinguished road
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
Chrisi1984 is offline   Reply With Quote

Old   June 2, 2010, 03:42
Default
  #8
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18
olesen will become famous soon enough
Quote:
Originally Posted by Chrisi1984 View Post
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,

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

int main(int argc, char *argv[])
{
....
}
Within hEqn.H you can then use it
Code:
...
    if (addDissipativeHeating)
    {
         hEqn -= thermalDissipationEff(turbulence);
    }

    hEqn.relax();
...
Presumably you'll have a switch in your fvSolution::SIMPLE dictionary to turn it on/off.
fredo490 and mm.abdollahzadeh like this.
olesen is offline   Reply With Quote

Old   June 2, 2010, 09:59
Default Thank you
  #9
Senior Member
 
Join Date: Jan 2010
Location: Stuttgart
Posts: 129
Rep Power: 7
Chrisi1984 is on a distinguished road
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
Chrisi1984 is offline   Reply With Quote

Old   June 2, 2010, 10:09
Default
  #10
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18
olesen will become famous soon enough
Quote:
Originally Posted by Chrisi1984 View Post
Hi Mark,

Thank you for your help!!!

I got it running.

But the solution is not like expected.
As I mentioned previously:
Quote:
Originally Posted by olesen View Post
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.
olesen is offline   Reply With Quote

Old   June 11, 2010, 03:27
Default
  #11
New Member
 
Johannes Kneer
Join Date: Mar 2009
Location: Germany, Karlsruhe
Posts: 13
Rep Power: 8
johannesk is on a distinguished road
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.
johannesk is offline   Reply With Quote

Old   June 11, 2010, 03:35
Default
  #12
Senior Member
 
Join Date: Jan 2010
Location: Stuttgart
Posts: 129
Rep Power: 7
Chrisi1984 is on a distinguished road
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
Chrisi1984 is offline   Reply With Quote

Old   July 12, 2010, 16:14
Default
  #13
Member
 
Sylvain Martel
Join Date: Apr 2009
Location: University of Sherbrooke/Quebec/Canada
Posts: 51
Rep Power: 8
smart is on a distinguished road
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
smart is offline   Reply With Quote

Old   July 13, 2010, 02:37
Default
  #14
Senior Member
 
Join Date: Jan 2010
Location: Stuttgart
Posts: 129
Rep Power: 7
Chrisi1984 is on a distinguished road
Hi,

Thank you for your advice.

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

that gives good results.

Regards
Chris
Chrisi1984 is offline   Reply With Quote

Old   May 6, 2013, 03:17
Default
  #15
Senior Member
 
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 236
Rep Power: 7
fredo490 is on a distinguished road
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:
Code:
- fvm::Sp(fvc::div(phi), h)
2) what is:
Code:
- 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 ( Big problem predicting adiabatic wall temperature (wing, mach 0.33, rhoSimplecFoam) )
Attached Images
File Type: png GoverningEquations.png (43.7 KB, 34 views)
fredo490 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
Problem with zeroGradient wall BC for temperature - Total temperature loss cboss OpenFOAM 10 March 5, 2015 07:57
Problem with rhoSimpleFoam : exploding enthalpy and density at the walls david39 OpenFOAM Running, Solving & CFD 6 January 18, 2011 12:49
Would you help explain the drop of total enthalpy? samxsl FLUENT 2 May 6, 2009 20:45
Fluent VOF Method - At a total loss advice required please LSF Main CFD Forum 5 April 13, 2009 21:56
meaning static enthalpy and total enthalpy mspark CD-adapco 1 April 1, 2004 06:20


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