rhoSimpleFoam loss in total enthalpy
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.
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.
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 | |
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
default Gauss linear;
grad(U) Gauss linear;
grad(p) Gauss linear;
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
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;
div(U,p) upwind phi;
// ************************************************** *********************** //
pMin pMin [ 1 -1 -2 0 0 0 0 ] 100;
p 0.01; // 0.3 -> 0.01
rho 0.9; // 0.05 -> 0.9
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?
Ok this looks ok. How bad is your mesh?
I think my mesh is alright:
internal faces: 321704
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
tet wedges: 0
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)
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.
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:
how to implement this in OpenFoam
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)
const compressible::turbulenceModel& turb
tmp<volTensorField> tgradU = fvc::grad(turb.U());
- ((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?
Instead, simply place it after the OpenFOAM includes and above main().
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.
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.
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);
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
eqnResidual = hEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
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.
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.
Thank you for your advice.
I yet had changed it to
div(phi,U) Gauss vanLeerV;
that gives good results.
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:
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 )
|All times are GMT -4. The time now is 17:32.|