
[Sponsors] 
May 3, 2020, 15:57 
Poor convergence: just mesh or numerics?

#1 
Member
Piotr Ładyński
Join Date: Apr 2017
Posts: 55
Rep Power: 8 
General problem
I’m having troubles with getting solution convergence using simpleFoam on quite simple channel flow problems and I can’t identify the reason behind this. I’ve read tons of pdf tutorials and I often encounter two statements, which confuse me: 1. Mesh has to be flawless in order to produce accurate solution and other factors don’t matter that much 2. Mesh can have some bad cell criteria, but the convergence might be slower, accuracy might be lower and it is okay, when it doesn’t blow up I always use Code:
checkMesh allGeometry allTopology If I try to use divSchemes other than upwind (which is not recommended option in general) usually everything crashes really fast (in case of bounded limited linearUpwind or limitedLinear it is bounding for k/epsilon/omega and the error accumulates up to infinity), even with very small underrelaxation factors. I try to increase them later if my simulation doesn’t crash after 100500 iterations, but often it doesn’t change much. Simulation doesn’t crash, but residuals won’t reduce below 1E3, which is way too high to agree on convergence. I tried some tips and schemes from [LINK], but it didn’t help much. I was expecting, that there is something wrong with my boundary conditions, but most of the time I use the most basic ones for (inlet→fixed velocity, fixed komegaespsilon , outlet→ fixed pressure, walls→ wallFunctions/noSlip) and I can’t find anything wrong. My simulations often work fine with 2d/axisymmetric domain. Could someone help me to identify my problem or make this clear? I struggle a lot with simple 3d geometries, but I want to put more complex shapes inside my flow channel. I really need a comment on that. Particular Problem I had some experience in succesfull transient buoyancy simulations with buoyantPimpleFoam before, but now I wanted to solve some steady flow problems and i’ve found it’s harder than I thought. This weekend I created this simple geometry with sudden expansion to a wide channel and sudden reduction with direction change (I somehow couldn’t figure out how to make it with hexablocks in gmsh/blockMesh, despite its simplicity). I knowingly ignored possibility of symmetry plane, because it doesn’t matter at this moment It is a pipe of ~Ø106mm diameter with Ø20mm connectors perpendicular to this pipe, the medium is cold water (nu= 1.1385e06 m^2/s) , the flow is 4.7778E4 m3/s (~0.76 m/s average on inlet, Re=~13356). BC’s I’m pretty sure that there is nothing special about these, and I haven’t messed it up. I used groovyBC to specify exponential flow profile, but my problem stays the same with the uniform fixed velocity vector Inlet U: exponential profile with groovy BC or uniform fixed value p: zeroGradient k: estimated fixedValue epsilon/omega: estimated fixedValue nut: calculated Velocity profile described with: Code:
valueExpression "vector(0,0.9505*pow((1(sqrt(pow(pos().z,2)+pow(pos().x,2))/0.01)),0.14286),0)"; U: inletOutlet (0 0 0) p: fixedValue of 0 k: zeroGradient epsilon/omega: zeroGradient nut: calculated Walls (Connector/Pipe) U: noSlip p: zeroGradient k: kqRWallFunction epsilon/omega: epsilonWallFunction/omegaWallFunction nut: nutkWallFuntion Results with cfMesh (cartesianMesh) The meshes generated with cfMesh always look great, but when you display checkMesh allGeometry allTopology it is almost sure, that there will be 34(+) failed checks, mostly in large numbers underdetermined cells and low quality faces/negative volume test and it it is very hard to avoid highly nonorthogonal faces or skewness. Poeple on the Internet told me that it is all acceptable if you can produce results (can I trust this?). The problem is that my initial residuals won’t fall below 5E3 or something like that (k, epsilon – upwind, U linearUpwindV grad(U), they seem too high to assume convergence. I switched velocity to upwind and after some randomly applied tweaks I had some results after >6300 iterations (the goal was set to only 1E4 for p and U, 1E3 for turbulence). I obtained a pressure drop of ~660Pa, which was ~120Pa more than I obtained using Solidworks Flow Simulation (dont’ know if it’s any better, because it’s not verbose on its numerics or exact mesh details). Results with snappyHexMesh Obtaining good meshes (good in terms of checkMesh quality criteria) with sMH is really time consuming even for simple geometries (and it is more frustrating, the simplier my channel shape seem to be), but you can end up with single underdetermined cells or low quality cell, which seem to be okay (I don’t count ‘concave cells’, because it is related to transition cells between mesh levels). In my example I had this problem, that pU resuiduals aren’t falling below 3E3: I let the solver flip matrices for 10k iterations (on upwind scheme for velocity, k, epsilon) and residuals weren’t changing much. Code:
Time = 10000 DILUPBiCG: Solving for Ux, Initial residual = 0.001280595, Final residual = 7.745998e09, No Iterations 13 DILUPBiCG: Solving for Uy, Initial residual = 0.0003662707, Final residual = 4.305333e09, No Iterations 11 DILUPBiCG: Solving for Uz, Initial residual = 0.001522483, Final residual = 3.394068e09, No Iterations 14 GAMG: Solving for p, Initial residual = 0.00307961, Final residual = 8.509195e08, No Iterations 16 GAMG: Solving for p, Initial residual = 0.0002415345, Final residual = 8.896399e08, No Iterations 7 GAMG: Solving for p, Initial residual = 2.354338e05, Final residual = 9.480076e08, No Iterations 4 GAMG: Solving for p, Initial residual = 3.640671e06, Final residual = 6.430156e08, No Iterations 3 GAMG: Solving for p, Initial residual = 8.784028e07, Final residual = 4.785081e08, No Iterations 2 GAMG: Solving for p, Initial residual = 2.556547e07, Final residual = 4.384352e08, No Iterations 1 GAMG: Solving for p, Initial residual = 8.643125e08, Final residual = 8.643125e08, No Iterations 0 time step continuity errors : sum local = 3.549486e07, global = 1.949192e10, cumulative = 3.79238e08 DILUPBiCG: Solving for epsilon, Initial residual = 0.0003076259, Final residual = 1.124087e10, No Iterations 5 bounding epsilon, min: 0.001938776 max: 5.461397 average: 0.1375029 DILUPBiCG: Solving for k, Initial residual = 0.0003674363, Final residual = 1.475334e09, No Iterations 5 ExecutionTime = 14769.84 s ClockTime = 14892 s blockMesh results I simplified my geometry, so the channel was straight (no direction change) and I could waste some time on blockMesh. I wanted to get ANY reliable results. I prepared one case with kEpsilon with y+ >30 and the second one with kOmegaSST for y+<1 with wall function disabled. The solution converged with upwind div schemes in both cases and the resulting pressure drop was similar (570580Pa), which was higher than you could calculate from DarcyWeisbach + local resistance coefficients (something like 495 Pa). Final thoughts What's wrong with my reasoning? Why does it take so many attempts to prepare ANY results for such a simple example (not speaking of more accurate ones)? Is it the mesh or some fvSolution options I should consider? Thanks for any feedback in advance. *i couldn't include stl or fms due to a size limit 

May 3, 2020, 17:04 

#2 
Senior Member
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 932
Rep Power: 12 
mesh+epsilon
__________________
The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki1/Wiki2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESIOpenCFDTrademark) Bugs/Features: FOAMExtend (WikkiFSB) Bugs: OpenFOAM.org How to create a MWE New: Forkable OpenFOAM mirror 

May 4, 2020, 02:02 

#3 
Member
Piotr Ładyński
Join Date: Apr 2017
Posts: 55
Rep Power: 8 
Thank you for your feedback, but could you be a tiny bit more specific about these two words?


May 7, 2020, 14:15 

#4 
Senior Member
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 126
Rep Power: 19 
Hi,
with upwind divscheme you can not expect any good results. Try to use something like linearUpwind, which is a quite good compromise between stability and accuracy. And by the way, residulas does not help much about judging the quality of the results. Compare with experimental or values from theory (what you have done). Best regards, Jan 

May 7, 2020, 16:28 

#5 
Member
Piotr Ładyński
Join Date: Apr 2017
Posts: 55
Rep Power: 8 
I already know, that upwind scheme is not very accurate, but it was the only scheme to produce any results: I couldn't get much stability with other schemes in my steady state simulations and this made me sure, that i'm doing something wrong.
If it is the mesh, that is not flawless enough (in terms of checkMesh allGeometry allTopology output), then I'd say it is close to impossible to have totally clean checkMesh report unless it is a very simple, smooth or purely cuboid shape without mesh level transitions (or just some 2D mesh). How people deal with this? In case of cfMesh even the tutorials examples have many bad cells. Some margin of error should be somehow acceptable I suppose (?). This needs more explaination to me. I try to apply general rules of thumb to my fv* files or BCs, based on tutorials, pdfs, forum search or user guide, but each time I think I finally understand the parameters I end up wasting my time to reconstruct oversimplified shapes with extensive blockMesh scripts or just making transient simulations 

May 7, 2020, 17:26 

#6 
Senior Member
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 932
Rep Power: 12 
>> I’m having troubles with getting solution convergence using simpleFoam on quite simple channel flow problems and I can’t identify the reason behind this.
As a first step, I would evaluate the convergence based on measuring a physical quantity within the numerical domain instead of monitoring the residuals. For example, probing the dominant velocity component at points within the field of interest. This paradigm shift may lead to the potential fact that your simulation does convergence actually. Second matter is whether the physical phenomenon you try to numerically model is actually steady, or involving some sort of periodicity or transient at some spatial/temporal location. I will continue to read your post with a fresh eye, but so far I haven't seen anything disturbing in either your setup or numerical modelling. Would you mind to run potentialFoam with very high number of correctors, and monitor minmax of velocity components by using the postProcess tool? It may provide something useful by indicating a hidden problem even though it is cheap enough to run. EDIT: Also can you try to use kOmegaSST instead of kEpsilon? Epsilon is notoriously bitchy near walls (what a scientific language I have just used!  in terms of boring language=because of the distance to wall term in the denominator of one of the governing equation terms which I find lazy to show explicitly ). Anyways, please do try kOmegaSST which is a robust generic closure model.
__________________
The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki1/Wiki2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESIOpenCFDTrademark) Bugs/Features: FOAMExtend (WikkiFSB) Bugs: OpenFOAM.org How to create a MWE New: Forkable OpenFOAM mirror 

May 8, 2020, 01:38 

#7 
Senior Member
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 126
Rep Power: 19 
Good point with kEps, kOmegaSST is definitely worth to try.
I checked your checkMesh results (at least for SHM). It is absolutely fine. The additional output due to options allGeometry and allTopology is normally not neccessary. In terms of stability nonorthogonality is most important (in general). Below 65 is absolute sufficient. Some comments on your fvSolution: set relTol for p to something like 0.01, nonOrthoCorr to 0, remove residual control and start with tight relaxation, e.g. 0.3 for p and 0.7 for U,k,omega and eps. And of course, change upwind Scheme to linearUpwind for U Best 

May 9, 2020, 17:42 

#8  
Member
Piotr Ładyński
Join Date: Apr 2017
Posts: 55
Rep Power: 8 
Quote:
potentialFoam postProcess output for my cfMesh (pictures): Code:
Time = 0 Reading fields: volVectorField: U Executing functionObjects fieldMinMax minMaxComponents(U) write: min(U) = (0.239974 0.974195 0.31967) in cell 177304 at location (0.0065529 0.052477 0.378724) max(U) = (0.368088 0.607389 0.738501) in cell 272015 at location (0.0043907 0.0523493 0.00863062) End Code:
Time = 0 Reading fields: volVectorField: U Executing functionObjects fieldMinMax minMaxComponents(U) write: min(U) = (0.716378 0.408556 1.18003) in cell 63109 at location (0.00641546 0.0525871 0.379634) max(U) = (0.220886 1.25981 0.186272) in cell 15145 at location (0.00896128 0.0527133 0.00362175) End I tried some scenarios with kOmegaSST, but without expected success. I attached some logs from my various modifications. Applied scenarios description: Code:
logOmegaSST_v01  turbulence model set to kOmegaSST,  relTol for p set to 0.01  velocity divScheme changed to linearUpwindV grad(U)  k, omega left with upwind (linearUpwind crashes in ~15 iterations)  snappyHexMesh mesh with y+ targeted at ~32 (one thick layer) the one I used before  forgot to log, then I saved some iterations from t = 6200, no convergence logOmegaSST_v02  same as v01, but I realised my estimation for inlet omega way too low, so i changed it (in equation: omega ~= epsilon/(Cmu*k) i forgot about dividing by Cmu)  forgot to log again, some iterations from t = 2500, no convergence stopped it, because not much changed logOmegaSST_v03  Switched to cartesianMesh (cfMesh)  log with 10k iterations (deleted time 3519899 to fit the file size limit), no convergence logOmegaSST_v04  remeshed with cfMesh to a target y+ of ~1 (11 layers)  resolved boundary layer conditions (fixedValue ~0 k and nut)  log of ~4k iterations (deleted time 351 tp 4049 to fit the file size limit)  i gave up calculations, because output was similar to previous attempt logOmegaSST_v05  changed relTol on velocity to a nonzero value (it was 0 before)  not much difference to initial residuals (relative to previous results), I abandonned further computations after ~3.5 iterations (removed time 301  3499) logOmegaSST_v06  switched k and omega divSchemes to linearUpwind  crashed after few iterations logOmegaSST_v07  remeshed with snappyHexMesh for better mesh quality control to y+ ~=1  linearUpwind divSchemes to k and omega  crashed after few iterations (omega explodes into infinity) logOmegaSST_v08  switched divSchemes for k and omega to upwind  crashed after few iterations again (omega again) 

May 10, 2020, 13:04 

#9 
Super Moderator
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,694
Blog Entries: 6
Rep Power: 50 
Dear Piotr,
well written post. Thumbs up. All information are given and in a clear and understandable way. Before I want to answer, is it possible to get a case we can investigate? In the first post, the 'STL' files are missing, so we cannot check your case. However, my first comment... I do not agree regarding not using upwind. It is the only scheme that produces unbounded physical values with the bottleneck of high diffusion (which smooths your solution); Taylor series expansion bla blub. I always start with upwind and switch to higher order schemes later. A common way people proceed in science and industry. Furthermore, your statement in the fvSolutions: Code:
The SIMPLECmethod relaxes the pressure in a consistent manner and additional relaxation of the pressure is not generally necessary Since when is the SIMPLEC relaxing the pressure? For my understanding regarding the pressuremomentum coupling:  SIMPLE is nonconsistent as a pressure term is neglected during the derivation.  SIMPLEC (c = consistent) includes this pressure term and makes it consistent in terms of mathematics I never made the check but using the keyword consistent activates a precalculation: Code:
if (simple.consistent()) { rAtU = 1.0/(1.0/rAU  UEqn.H1()); phiHbyA += fvc::interpolate(rAtU()  rAU)*fvc::snGrad(p)*mesh.magSf(); HbyA = (rAU  rAtU())*fvc::grad(p); }
__________________
Keep foaming, Tobias Holzmann 

May 11, 2020, 13:18 

#10 
Member
Piotr Ładyński
Join Date: Apr 2017
Posts: 55
Rep Power: 8 
I couldn't fit into attachment size limt at first so I didn't put my geometry, but here is a link from my google drive (my triSurface folder for this case setup plus optional .fms file for cfMesh).
I must admit I haven't looked into the source code for SIMPLE vs SIMPLEC differences before. My comment on this in fvSolution was literal copypasted sentence from this pdf at wolfdynamics.com (slide 665 and 666... and 667 with advice to put a bit more conservative underrelaxation on SIMPLEC). When JNSN told me to put tight relaxation I switched to SIMPLE for my last simulations. 

May 12, 2020, 11:05 

#11 
Senior Member
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 126
Rep Power: 19 
Hi Piotr,
I just run your case, with proper relaxation it runs absolutely fine (at least in terms of stability) with linearUpwind, for kEps and kOmega, both SIMPLE and SIMPLEC. I agree, convergence is not good. But I think it is not due to bad mesh quality or numerics, the flow is not stationary. May after running steady state, switch to transient. Best, Jan 

Tags 
convergence criteria, mesh qualty, numerical scheme 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
sliding mesh problem in CFX  Saima  CFX  46  September 11, 2021 07:38 
[snappyHexMesh] High quality mesh for wind in complex urban environment  ziboaa  OpenFOAM Meshing & Mesh Conversion  1  January 12, 2021 15:33 
decomposePar problem: Cell 0contains face labels out of range  vaina74  OpenFOAM PreProcessing  37  July 20, 2020 05:38 
[ICEM] Bad Quality Hybrid Mesh External Flow  tim13  ANSYS Meshing & Geometry  0  March 8, 2020 02:22 
[ICEM] Problem making structured mesh on a surface  froztbear  ANSYS Meshing & Geometry  4  November 10, 2011 08:52 