[piotr.mecht]

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

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:

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.

 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

[HPE]

mesh+epsilon

[piotr.mecht]

Thank you for your feedback, but could you be a tiny bit more specific about these two words?

[JNSN]

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

[piotr.mecht]

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

[HPE]

>> 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]

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

[piotr.mecht]

Originally Posted by HPE

>> 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):

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):

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:

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)

[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:

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:

    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.

[piotr.mecht]

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.

[JNSN]

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