
[Sponsors] 
May 23, 2010, 07:37 
rhoSimpleFoam loss in total enthalpy

#1 
Senior Member
Join Date: Jan 2010
Location: Stuttgart
Posts: 146
Rep Power: 9 
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 

May 23, 2010, 12:03 

#2 
Senior Member
BastiL
Join Date: Mar 2009
Posts: 488
Rep Power: 13 
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 

May 23, 2010, 13:15 

#3 
Senior Member
Join Date: Jan 2010
Location: Stuttgart
Posts: 146
Rep Power: 9 
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 

May 23, 2010, 17:01 

#4 
Senior Member
BastiL
Join Date: Mar 2009
Posts: 488
Rep Power: 13 
Ok this looks ok. How bad is your mesh?
Regards Bastian 

May 23, 2010, 18:49 

#5 
Senior Member
Join Date: Jan 2010
Location: Stuttgart
Posts: 146
Rep Power: 9 
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 

May 31, 2010, 05:12 

#6  
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 798
Rep Power: 20 
Quote:
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() ); } 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 // 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() ); } 

June 2, 2010, 03:23 
how to implement this in OpenFoam

#7 
Senior Member
Join Date: Jan 2010
Location: Stuttgart
Posts: 146
Rep Power: 9 
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 

June 2, 2010, 03:42 

#8  
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 798
Rep Power: 20 
Quote:
Instead, simply place it after the OpenFOAM includes and above main(). Eg, Code:
tmp<volScalarField> thermalDissipationEff( ... ) { ... } int main(int argc, char *argv[]) { .... } Code:
... if (addDissipativeHeating) { hEqn = thermalDissipationEff(turbulence); } hEqn.relax(); ... 

June 2, 2010, 09:59 
Thank you

#9 
Senior Member
Join Date: Jan 2010
Location: Stuttgart
Posts: 146
Rep Power: 9 
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 

June 2, 2010, 10:09 

#10  
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 798
Rep Power: 20 
Quote:
See if there is a sign error or something else, or if particular bits of the domain are doing something strange. 

June 11, 2010, 03:27 

#11 
New Member
Johannes Kneer
Join Date: Mar 2009
Location: Germany, Karlsruhe
Posts: 13
Rep Power: 10 
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. 

June 11, 2010, 03:35 

#12 
Senior Member
Join Date: Jan 2010
Location: Stuttgart
Posts: 146
Rep Power: 9 
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 

July 12, 2010, 16:14 

#13 
Member
Sylvain Martel
Join Date: Apr 2009
Location: University of Sherbrooke/Quebec/Canada
Posts: 51
Rep Power: 9 
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 

July 13, 2010, 02:37 

#14 
Senior Member
Join Date: Jan 2010
Location: Stuttgart
Posts: 146
Rep Power: 9 
Hi,
Thank you for your advice. I yet had changed it to div(phi,U) Gauss vanLeerV; that gives good results. Regards Chris 

May 6, 2013, 03:17 

#15 
Senior Member
HECKMANN Frédéric
Join Date: Jul 2010
Posts: 237
Rep Power: 9 
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 ) 

Thread Tools  
Display Modes  


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  CDadapco  1  April 1, 2004 06:20 