CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

Poor convergence: just mesh or numerics?

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By piotr.mecht
  • 1 Post By HPE

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 3, 2020, 15:57
Default Poor convergence: just mesh or numerics?
  #1
Member
 
Piotr Ładyński
Join Date: Apr 2017
Posts: 51
Rep Power: 5
piotr.mecht is on a distinguished road
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
to verify my meshes. Totally clean output can be achieved only with simple meshes with blockMesh or gmsh, or some axisymmetric/2D, but when you use mesh with grading then you have at least concave cells on transition between mesh levels. From what Iíve read it doesnít matter much in modern OpenFOAM versions, so I can ignore it. Other check fails also supposedly donít matter if my mesh isnít highly skew or non-orhogonal (itís never stressed in the texts how bad is not bad, so it is also confusing), but I always try to find the best output anyway.


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 100-500 iterations, but often it doesnít change much. Simulation doesnít crash, but residuals wonít reduce below 1E-3, 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 k-omega-espsilon ,
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 hexa-blocks 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.1385e-06 m^2/s) , the flow is 4.7778E-4 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)";
Outlet
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 3-4(+) failed checks, mostly in large numbers underdetermined cells and low quality faces/negative volume test and it it is very hard to avoid highly non-orthogonal 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 5E-3 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 1E-4 for p and U, 1E-3 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 p-U resuiduals arenít falling below 3E-3: 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.745998e-09, No Iterations 13
 DILUPBiCG:  Solving for Uy, Initial residual = 0.0003662707, Final residual = 4.305333e-09, No Iterations 11
 DILUPBiCG:  Solving for Uz, Initial residual = 0.001522483, Final residual = 3.394068e-09, No Iterations 14
 GAMG:  Solving for p, Initial residual = 0.00307961, Final residual = 8.509195e-08, No Iterations 16
 GAMG:  Solving for p, Initial residual = 0.0002415345, Final residual = 8.896399e-08, No Iterations 7
 GAMG:  Solving for p, Initial residual = 2.354338e-05, Final residual = 9.480076e-08, No Iterations 4
 GAMG:  Solving for p, Initial residual = 3.640671e-06, Final residual = 6.430156e-08, No Iterations 3
 GAMG:  Solving for p, Initial residual = 8.784028e-07, Final residual = 4.785081e-08, No Iterations 2
 GAMG:  Solving for p, Initial residual = 2.556547e-07, Final residual = 4.384352e-08, No Iterations 1
 GAMG:  Solving for p, Initial residual = 8.643125e-08, Final residual = 8.643125e-08, No Iterations 0
 time step continuity errors : sum local = 3.549486e-07, global = 1.949192e-10, cumulative = 3.79238e-08
 DILUPBiCG:  Solving for epsilon, Initial residual = 0.0003076259, Final residual = 1.124087e-10, 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.475334e-09, No Iterations 5
 ExecutionTime = 14769.84 s  ClockTime = 14892 s
Last saved results (from 9700th, 9800th, 9900th , 10000th) looked a bit oscilating. Iíd say you could read pressure drop of 520-530 Pa from it, which was very close to my previous SolidWorks Flow Simulation results, but it doesnít seem to be reliable if you read residuals so high. Computations took a long time tho and it wouldnít be worth an effort in a proffessional situation.


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 k-Epsilon with y+ >30 and the second one with k-OmegaSST for y+<1 with wall function disabled. The solution converged with upwind div schemes in both cases and the resulting pressure drop was similar (570-580Pa), which was higher than you could calculate from Darcy-Weisbach + 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
Attached Images
File Type: png pressuredropMeshView.png (56.8 KB, 19 views)
File Type: png pressureDropBlockMesh.png (47.1 KB, 20 views)
File Type: png pressureDropCfMesh.png (45.5 KB, 14 views)
Attached Files
File Type: zip checkMeshLogs.zip (5.7 KB, 3 views)
File Type: zip myFiles.zip (17.4 KB, 5 views)
piotr.mecht is offline   Reply With Quote

Old   May 3, 2020, 17:04
Default
  #2
HPE
Senior Member
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 583
Rep Power: 5
HPE is on a distinguished road
mesh+epsilon
HPE is online now   Reply With Quote

Old   May 4, 2020, 02:02
Default
  #3
Member
 
Piotr Ładyński
Join Date: Apr 2017
Posts: 51
Rep Power: 5
piotr.mecht is on a distinguished road
Thank you for your feedback, but could you be a tiny bit more specific about these two words?
Tobi likes this.
piotr.mecht is offline   Reply With Quote

Old   May 7, 2020, 14:15
Default
  #4
Member
 
JNSN's Avatar
 
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 80
Rep Power: 15
JNSN is on a distinguished road
Hi,

with upwind div-scheme 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
JNSN is offline   Reply With Quote

Old   May 7, 2020, 16:28
Default
  #5
Member
 
Piotr Ładyński
Join Date: Apr 2017
Posts: 51
Rep Power: 5
piotr.mecht is on a distinguished road
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
piotr.mecht is offline   Reply With Quote

Old   May 7, 2020, 17:26
Default
  #6
HPE
Senior Member
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 583
Rep Power: 5
HPE is on a distinguished road
>> 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 min-max 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.
JNSN likes this.
HPE is online now   Reply With Quote

Old   May 8, 2020, 01:38
Default
  #7
Member
 
JNSN's Avatar
 
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 80
Rep Power: 15
JNSN is on a distinguished road
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 non-orthogonality 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
JNSN is offline   Reply With Quote

Old   May 9, 2020, 17:42
Default
  #8
Member
 
Piotr Ładyński
Join Date: Apr 2017
Posts: 51
Rep Power: 5
piotr.mecht is on a distinguished road
Quote:
Originally Posted by HPE View Post
>> Would you mind to run potentialFoam with very high number of correctors, and monitor min-max of velocity components by using the postProcess tool? It may provide something useful.
Than you for your answers.


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
potentialFoam postProcess output for my (images would look more or less the same):
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
Max values are on the edges of expansion and contraction.



I tried some scenarios with k-OmegaSST, 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 351-9899 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 non-zero 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)
Attached Images
File Type: png PotentialFoamFullScale.png (19.4 KB, 3 views)
File Type: png PotentialFoamMinMax.png (41.7 KB, 4 views)
File Type: png PotentialFoamRescaled.png (41.2 KB, 4 views)
Attached Files
File Type: zip someLogs.zip (169.6 KB, 1 views)
piotr.mecht is offline   Reply With Quote

Old   May 10, 2020, 13:04
Default
  #9
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,383
Blog Entries: 6
Rep Power: 41
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
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 pressure-momentum coupling:


- SIMPLE is non-consistent 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);                                    
    }
There is nothing related to relaxation. Thus, as your case does not have any pressure field relaxation, the turbulence properties are still bounding, I expect that you do have fluctuations that do not allow your residuals to go down. Additionally, the residuals are not always identical. It is always a matter of choice. Foam uses the L2 norm if I am not wrong here. There are different norms. And, as HPE already said, the residual output is not always a instrument of measurement.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   May 11, 2020, 13:18
Default
  #10
Member
 
Piotr Ładyński
Join Date: Apr 2017
Posts: 51
Rep Power: 5
piotr.mecht is on a distinguished road
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 SIMPLE-C differences before. My comment on this in fvSolution was literal copy-pasted sentence from this pdf at wolfdynamics.com (slide 665 and 666... and 667 with advice to put a bit more conservative underrelaxation on SIMPLE-C). When JNSN told me to put tight relaxation I switched to SIMPLE for my last simulations.
piotr.mecht is offline   Reply With Quote

Old   May 12, 2020, 11:05
Default
  #11
Member
 
JNSN's Avatar
 
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 80
Rep Power: 15
JNSN is on a distinguished road
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
JNSN is offline   Reply With Quote

Reply

Tags
convergence criteria, mesh qualty, numerical scheme

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[ICEM] Bad Quality Hybrid Mesh External Flow tim13 ANSYS Meshing & Geometry 0 March 8, 2020 02:22
[snappyHexMesh] High quality mesh for wind in complex urban environment ziboaa OpenFOAM Meshing & Mesh Conversion 0 January 14, 2020 12:15
decomposePar problem: Cell 0contains face labels out of range vaina74 OpenFOAM Pre-Processing 32 February 7, 2018 08:26
sliding mesh problem in CFX Saima CFX 45 September 22, 2015 10:53
[ICEM] Problem making structured mesh on a surface froztbear ANSYS Meshing & Geometry 4 November 10, 2011 08:52


All times are GMT -4. The time now is 08:48.