rhoSimpleFoam loss in total enthalpy
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 
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 
1 Attachment(s)
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((rhoA(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*(1A(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 1e08; relTol 0.05; smoother GaussSeidel; cacheAgglomeration false; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; maxIter 500; } U { solver smoothSolver; smoother GaussSeidel; nSweeps 2; tolerance 1e06; relTol 0.1; } h { solver PBiCG; preconditioner DILU; tolerance 1e06; relTol 0.1; } k { solver smoothSolver; smoother GaussSeidel; nSweeps 2; tolerance 1e07; relTol 0.1; } epsilon { solver smoothSolver; smoother GaussSeidel; nSweeps 2; tolerance 1e07; 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 
Ok this looks ok. How bad is your mesh?
Regards Bastian 
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 (nonclosed singly connected) wout 2000 2050 ok (nonclosed singly connected) v 406 432 ok (nonclosed singly connected) win 2000 2050 ok (nonclosed singly connected) wr 9500 9550 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.267792 0.032787 0.0440929) (0.505012 0.0488568 0.302414) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (2.79925e19 6.46543e19 4.54471e19) OK. Max cell openness = 2.56061e16 OK. Max aspect ratio = 6.1821 OK. Minumum face area = 9.30065e07. Maximum face area = 2.2486e05. Face area magnitudes OK. Min volume = 2.33102e09. Max volume = 9.44257e08. Total volume = 0.00176285. Cell volumes OK. Mesh nonorthogonality Max: 56.5682 average: 8.99531 Nonorthogonality check OK. Face pyramids OK. Max skewness = 1.92253 OK. Mesh OK. Regards Chrisi 
Quote:
Eg, Code:
// The source for the enthalpy equation resulting from It could be, however, that the Fluent temperatures are still higher. From their docs, it appears that they may doing the same as cdadapco 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 
how to implement this in OpenFoam
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 
Quote:
Instead, simply place it after the OpenFOAM includes and above main(). Eg, Code:
tmp<volScalarField> thermalDissipationEff( ... ) Code:
... 
Thank you
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 
Quote:
Quote:

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. 
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 
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 
Hi,
Thank you for your advice. I yet had changed it to div(phi,U) Gauss vanLeerV; that gives good results. Regards Chris 
1 Attachment(s)
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) Code:
 fvc::div(phi, u2) 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.cfdonline.com/Forums/ope...mplecfoam.html ) 
All times are GMT 4. The time now is 06:37. 