
[Sponsors] 
February 15, 2017, 05:23 
Bounding k and omega problems

#1 
New Member
Join Date: Mar 2015
Posts: 12
Rep Power: 11 
Hi,
I'm running a case in OpenFOAM (v1606+) where I have trouble with solution convergence and bounding of k and omega. The case is a oscillating cover (looks like a 2D cut of a bowler hat) in infinite fluid. There is no inflow or outflow, just still water and the cover oscillating in the ydirection. The simulation is obviously transient. For the mesh motion I'm using dynamicMotionSolverFvMesh with the displacementLaplacian solver and InverseDistance diffusivity model. For turbulence I'm using the komega SST model, without wall functions. For boundary conditions I have zeroGradient for U and P on all the domain sides (except for front and back which are empty). On the cover I have zeroGradient for P and movingWallVelocity with uniform value (0,0,0) for U. k and omega also have zeroGradient on all sides, but with fixedValue uniform 1e11; on the cover. nut have the same BC on the cover, while the other fields are set to calculated. I have been trying to have second order schemes, but with central schemes I get bounding problems with k and omega and the solution crashes. I just got the solution to run by switching to the following scheme: Code:
ddtSchemes { default Euler; } gradSchemes { default cellMDLimited Gauss linear 0.5; } divSchemes { default none; div(phi,U) Gauss inearUpwind cellMDLimited Gauss linear 0.5; div(phi,k) Gauss linearUpwind cellMDLimited Gauss linear 0.5; div(phi,omega) Gauss linearUpwind cellMDLimited Gauss linear 0.5; div((nuEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear limited corrected 0.5; laplacian(diffusivity,cellMotionU) Gauss linear limited corrected 0.5; } interpolationSchemes { default linear; } snGradSchemes { default limited corrected 0.5; } wallDist { method meshWave; } fluxRequired { default no; pcorr ; p ; } Code:
solvers { pcorr { solver GAMG; tolerance 0.02; relTol 0; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; cacheAgglomeration no; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } p { $pcorr; tolerance 1e7; relTol 0.01; } pFinal { $p; tolerance 1e07; relTol 0; } U { solver PBiCG; preconditioner DILU; tolerance 1e08; relTol 0;//0.1; minIter 3; } UFinal { $U; //tolerance 1e05; //relTol 0; } "cellDisplacement.*" { solver GAMG; tolerance 1e8; relTol 0; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } "(komega)" { solver PBiCG; preconditioner DILU; tolerance 1e08; reltol 0; minIter 3; } "(komega)Final" { solver PBiCG; preconditioner DILU; tolerance 1e08; reltol 0; minIter 3; } } PIMPLE { correctPhi yes; nOuterCorrectors 2; nCorrectors 1; nNonOrthogonalCorrectors 1; pRefCell 1001; pRefValue 0; } relaxationFactors { fields { } equations { "U.*" 1; "k.*" 1; “omega.*” 1; } } Code:
Time = 0 Mesh stats points: 64712 internal points: 0 faces: 127108 internal faces: 63212 cells: 31584 faces per cell: 6.02584 boundary patches: 6 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 30474 prisms: 202 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 908 Breakdown of polyhedra by number of faces: faces number of cells 7 798 8 110 Checking topology... Boundary definition OK. Cell to face addressing 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 inlet 40 82 ok (nonclosed singly connected) outlet 40 82 ok (nonclosed singly connected) topAndBottom 80 164 ok (nonclosed singly connected) front 31584 32356 ok (nonclosed singly connected) back 31584 32356 ok (nonclosed singly connected) cover 568 1136 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.625 0.625 0.00625) (0.625 0.625 0.02125) Mesh has 2 geometric (nonempty/wedge) directions (1 1 0) Mesh has 2 solution (nonempty) directions (1 1 0) All edges aligned with or perpendicular to nonempty directions. Boundary openness (7.57749e19 9.55767e20 7.24336e15) OK. Max cell openness = 4.09121e16 OK. Max aspect ratio = 24.2247 OK. Minimum face area = 4.86616e08. Maximum face area = 0.000985343. Face area magnitudes OK. Min volume = 7.29923e10. Max volume = 1.47801e05. Total volume = 0.0234301. Cell volumes OK. Mesh nonorthogonality Max: 35.522 average: 3.79143 Nonorthogonality check OK. Face pyramids OK. Max skewness = 1.90833 OK. Coupled point location match (average 0) OK. Mesh OK. Thanks! 

February 15, 2017, 10:42 

#2 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17 
Hi,
not still sure about the BCs although I read the description twice, but anyways... a] I am not very sure about the simulation with komega at the walls without wall functions employed. I am afraid you simply need to have a wall function involved by the definition of the model. Please refer to the list of functions here: http://openfoam.com/documentation/cp...Functions.html b] Second order  the linearUpwind is actually second order scheme, would you like to use more precise schemes? Your mesh is actually of a great quality so I would not worry much about the limitation of nonorthogonality. Not needed for your case. I wonder you'd be able to go to more precise schemes on turbulence in OpenFOAM. Actually what would happen if your omega would be upwinded only? Would the bounding problem remain? c] Instead of PBiCG you may try PBiCGStab  it is a new stabilised version d] general though  you are going to have a problem with k and omega if your U gradients are unphysical. Therefore your velocity field should be well converged for each time step. Is it the case? e] Maybe others will spot different problem. 

February 15, 2017, 11:11 

#3 
New Member
Join Date: Mar 2015
Posts: 12
Rep Power: 11 
Thanks for the answer Matej.
Sorry for being a bit unclear on the BCs, here are my BCs for k, nut, omega, p and U. edit: note that the left and right sides of the domain is called inlet and outlet even though there's now flow through them. Code:
FoamFile { version 2.0; format ascii; class volScalarField; location "0"; object k; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 2 0 0 0 0]; internalField uniform 1e11; boundaryField { inlet { type zeroGradient; } outlet { type zeroGradient; } topAndBottom { type zeroGradient; } front { type empty; } back { type empty; } cover { type fixedValue; value uniform 1e11; } } Code:
FoamFile { version 2.0; format ascii; class volScalarField; location "0"; object nut; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 1 0 0 0 0]; internalField uniform 1e11; boundaryField { inlet { type calculated; value uniform 1e11; } outlet { type calculated; value uniform 1e11; } topAndBottom { type calculated; value uniform 1e11; } front { type empty; } back { type empty; } cover { type fixedValue; value uniform 1e11; } } Code:
FoamFile { version 2.0; format ascii; class volScalarField; location "0"; object omega; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 0 1 0 0 0 0]; internalField uniform 1; boundaryField { inlet { type zeroGradient; } outlet { type zeroGradient; } topAndBottom { type zeroGradient; } front { type empty; } back { type empty; } cover { type fixedValue; value uniform 1e11; } } Code:
FoamFile { version 2.0; format ascii; class volScalarField; location "0"; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 2 0 0 0 0]; internalField uniform 0; boundaryField { inlet { type zeroGradient; } outlet { type zeroGradient; } topAndBottom { type zeroGradient; } front { type empty; } back { type empty; } cover { type zeroGradient; } } Code:
FoamFile { version 2.0; format ascii; class volVectorField; location "0"; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { inlet { type zeroGradient; } outlet { type zeroGradient; } topAndBottom { type zeroGradient; } front { type empty; } back { type empty; } cover { type movingWallVelocity; value uniform (0 0 0); } } b) I do know that linearUpwind is second order, so Im happy to settle with it. Ideally I would prefer central schemes, but I have not managed to get it to run, so the next step was to switch to linearUpwind and check if that helped. It did, as the simulations now are running without crashing, but I still have the boundedness problems as mentioned. I have not tried to use first order upwind schemes for omega, but if I nothing else helps I guess that's the only solution. Why just for omega btw? and not for both omega and k? c) I'll try PBiCGStab instead, thanks. d) I have no problems with convergence in U. The initial residual is reasonably low, and it drops down to ~10^15 in a couple of iterations. Last edited by BendikS; February 15, 2017 at 11:15. Reason: clarifications 

February 16, 2017, 03:43 

#4 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 182
Rep Power: 17 
Hi,
one more thing, I do not know where did you got these lines: div(phi,U) Gauss inearUpwind cellMDLimited Gauss linear 0.5; div(phi,k) Gauss linearUpwind cellMDLimited Gauss linear 0.5; div(phi,omega) Gauss linearUpwind cellMDLimited Gauss linear 0.5; but it does not make much sense any more. The current style is: div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss linearUpwind grad(k); div(phi,omega) Gauss uwind; If you wish to use purely linear scheme, as it is the holy grail of second order, you are always going to have problems with boundedness, as the scheme is bounded only in 1D situation. So if you manage to keep your Co<0.1 approx and Pe<=2 you should be OK. But it is hardly ever achievable. Even for LES we use LUST, which blends linear with linearUpwind. Second thing: if you insist on accuracy of your spacial discretisation you should be demanding in the same way on your time discretisation. Using Euler (first order  very dispersive on your oscilations), you should use CrankNicolson scheme. Last edited by matejfor; February 16, 2017 at 06:20. 

February 16, 2017, 05:03 

#5 
New Member
Join Date: Mar 2015
Posts: 12
Rep Power: 11 
Hi,
I believed those two statements were the same, but in the latest simulations I have been using the "good old" version. I switched to Upwind for omega now, and the bounding issue for omega stopped, but k is still having some "oscillatory" bounding issues, i.e. its fine for maybe ten iterations, but then k is negative for one iterations or two, before it's behaving well again. Regarding central schemes, I guess it will be to hard to get stable solutions with it. yup, switching to crank nicolson aswell, just wanted to figure out the spatial schemes first. 

February 17, 2017, 05:06 

#6 
New Member
Join Date: Mar 2015
Posts: 12
Rep Power: 11 
Ok so I switched to
ddt: CrankNicolson 0.9; grad: default cellMDLimited Gauss Upwind 0.5; div(U): linearUpwind grad(U); div(k): upwind; div(omega): upwind; and gauss linear corrected for laplacian schemes. The simulations runs, and I'm not having the bounding problems anymore, but suddenly around the t = 8s, the pressure residual increases and behaves a bit strangely considering that the cover is just continuing to oscillate. There is nothing "new" happening in the simulations physically. I attached the residual plot. Any thoughts? Thanks 

April 15, 2021, 16:29 

#7 
New Member
David Smith
Join Date: Jul 2013
Posts: 9
Rep Power: 13 
I got the same issue.


December 20, 2021, 04:38 

#8 
Member
Arthur
Join Date: Aug 2014
Location: Italy
Posts: 47
Rep Power: 12 
From my experience:
sometimes chosing the second order scheme for pressure can give some issued with k an omega too (even if upwind is used for turbulence). I tried a 3D case compressible, supersonic with shockWaves within the domain: using MUSCL for div(phid, p) in some cases, not always, I got some troubles with bounding k and omega. In some other cases also p was bounded in small regions of the domain. changing the scheme for div(phid, p) to LUST everything stabilised very quickly and work smoothly. I also got the impression that the choice of the numerical scheme is depend on the type of processor. Somebody can confirm it? 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
transsonic nozzle with rhoSimpleFoam  Unseen  OpenFOAM Running, Solving & CFD  8  July 1, 2022 06:54 
Behaviour of the kOmegaSST in a steadystate case  Max1234  OpenFOAM Running, Solving & CFD  18  October 31, 2018 08:03 
Bounding epsilon or bounding omega  Stylianos  OpenFOAM  8  February 23, 2018 13:41 
Unexpected deltaT decrease in pimpleFoam simulation  robyTKD  OpenFOAM Running, Solving & CFD  9  June 27, 2014 06:52 
Bounding OMEGA  barath.ezhilan  OpenFOAM  3  April 20, 2012 11:06 