[philippose]

Hello,
Could someone give me more information regarding the implication of the "Time Step Continuity Error" output which is displayed while solving a steady state problem using simpleFoam ?

What exactly does it mean, and what kind of values can be tolerated? And also, what factors does it depend on ?

Thank you

Philippose

[pierre]

Hi Philippose,

I'll try to help, the code goes:

volScalarField contErr = fvc::div(phi);

scalar sumLocalContErr = runTime.deltaT().value()*
mag(contErr)().weightedAverage(mesh.V()).value();


Info<< "time step continuity errors : sum local = " << sumLocalContErr


since typically runTime.deltaT().value() = 1 for simpleFoam calculations (i.e. iteration counter), div(phi) is your face flux (interpolated velocity field) divergence and ideally is zero (if continuity is something you're looking for). The magnitude of div(phi) volume averaged over each cell should hence also tend to zero. If you want you can normalise it further to give you an idea of magnitude (maybe by dividing by average inlet face flux).

Hope this helps

Pierre

[philippose]

Hello Pierre,
Thanks a lot for the information!! I had not looked into the code to try and figure out the meaning of the term. I assumed it was a standard term used in CFD solvers, which one was expected to know.

Thank you for the explanation :-) Have a nice day!

Philippose

[roger]

Hi Foamers!

I'm new in OF and I'm having problems in running a small mesh (around 400.000 tetra cells, 20 MB) with simpleFoam. Time step continuity error rises too much, and it happens suddently. It takes so long to simulate just 20 seconds (about 1 hour). Anyone can help me having a look at this? I've tried changing tolerances, solver types, ... Thanks in advance!


Time = 9

DILUPBiCG: Solving for Ux, Initial residual = 0.458639, Final residual = 0.0233617, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.379868, Final residual = 0.0171841, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.658214, Final residual = 0.0397693, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.727959, Final residual = 9.87705e-07, No Iterations 45
time step continuity errors : sum local = 135.971, global = -1.13269e-06, cumulative = -1.30202e-06
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 4.32485e-15, No Iterations 1
bounding epsilon, min: -4.46707e+08 max: 2.69515e+25 average: 6.1936e+19
DILUPBiCG: Solving for k, Initial residual = 0.999987, Final residual = 8.72274e-06, No Iterations 1
ExecutionTime = 1349.42 s ClockTime = 2229 s

Time = 10

DILUPBiCG: Solving for Ux, Initial residual = 0.0742623, Final residual = 8.17034e-09, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0742549, Final residual = 5.84205e-09, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0742894, Final residual = 5.0114e-09, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.511445, Final residual = 9.97221e-07, No Iterations 487
time step continuity errors : sum local = 102.865, global = 2.62581, cumulative = 2.62581
DILUPBiCG: Solving for epsilon, Initial residual = 0.522913, Final residual = 1.72356e-17, No Iterations 1
bounding epsilon, min: -3.41487e+17 max: 8.29693e+24 average: 1.98496e+19
DILUPBiCG: Solving for k, Initial residual = 0.887151, Final residual = 0.0281178, No Iterations 1
bounding k, min: -2.08013e+12 max: 5.53603e+18 average: 4.26882e+13
ExecutionTime = 1491.63 s ClockTime = 2487 s

Time = 11

DILUPBiCG: Solving for Ux, Initial residual = 0.999971, Final residual = 0.0517584, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.99998, Final residual = 0.0517602, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.999913, Final residual = 0.0517559, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.859263, Final residual = 3.05651e-08, No Iterations 1
time step continuity errors : sum local = 5.63857e+13, global = 0.469436, cumulative = 3.09525
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 8.3548e-18, No Iterations 1
bounding epsilon, min: -1.881e+17 max: 1.25296e+50 average: 2.87937e+44
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 4.56399e-08, No Iterations 1
bounding k, min: -4.21771e+12 max: 2.1231e+34 average: 9.18412e+28
ExecutionTime = 1511.3 s ClockTime = 2559 s

Time = 12

DILUPBiCG: Solving for Ux, Initial residual = 0.0744755, Final residual = 1.84446e-11, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0744755, Final residual = 1.90422e-11, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0744755, Final residual = 1.69176e-11, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.481839, Final residual = 438.302, No Iterations 1001
time step continuity errors : sum local = 7.48027e+23, global = 6.2141e+20, cumulative = 6.2141e+20
DILUPBiCG: Solving for epsilon, Initial residual = 0.163347, Final residual = 9.82397e-16, No Iterations 1
bounding epsilon, min: -2.10117e+37 max: 4.04364e+51 average: 4.20323e+47
DILUPBiCG: Solving for k, Initial residual = 6.9407e-13, Final residual = 6.9407e-13, No Iterations 0
ExecutionTime = 1761.29 s ClockTime = 2919 s

Time = 13

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.000907592, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.000907587, No Iterations 2
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.000907558, No Iterations 2
DICPCG: Solving for p, Initial residual = 1, Final residual = 5.12291e-10, No Iterations 1
time step continuity errors : sum local = 3.26533e+37, global = -3.23767e+21, cumulative = -2.61627e+21
DILUPBiCG: Solving for epsilon, Initial residual = 1, Final residual = 2.1039e-20, No Iterations 1
bounding epsilon, min: -1.985e+21 max: 1.16326e+89 average: 2.67323e+83
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 0.000185863, No Iterations 1
bounding k, min: 2.99871e-16 max: 2.54158e+68 average: 5.8407e+62
ExecutionTime = 1780.3 s ClockTime = 3165 s

[paulo]

Hello Mayol and all,

As simpleFoam will look for the steady state solution, I usually apply a big epsilon starting value in the domain (about 100 times the value of a boundary) and wait for the solver to "change" it. You'll see that the values go down quickly, and the solution almost always gets stabilized.

Give a feedback if it worked for you.

Any comments are appreciated.

Regards,

Paulo Rocha.

[naish]

Hello Paulo ;
I have the same problem and just trying your hint,
do you mean the epsilon not the relaxationsfactor of epsilon right ??
I read in other threads to keep the relax small.
Do you mean with 100 times the boundary value inlet velocity ??

To mayol It takes such a long time because of the itterations of P go and change the solver from PBiCK to GAMG

[paulo]

Hello Naish and all,

Answers:

1 - Yes, it is the starting value of epsilon in the domain (the value contained in the 'epsilon' file in the '0' folder). About the relaxation factor I use 0.3 (also for k), and don't have problems. You can try a smaller if you need;

2 - Yes, 100 times the value of epsilon at the inlet (beware,not for 'k').

Regards,

Paulo Rocha.

[roger]

Hey guys!!

I've tried changing epsilon value (I have velocity inlet 5m/s, so I consider epsilon = 500), reducing relaxation factors to 0.3 (U,p,k,E,...), and also changing the solver for avoiding to much iterations on pressure (PBICK to GAMG). It doesn't work. And it doesn't allow me to change the solver for pressure, I have the following report:

Starting time loop

Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.0118552, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.0130391, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.0149314, No Iterations 1


--> FOAM FATAL IO ERROR : expected 'true' or 'false', found yes

file: /home/roger/OpenFOAM/roger-1.4.1/run/tutorials/simpleFoam//Senzilla/system/fvSol ution::p::scaleCorrection at line 41.

From function operator>>(Istream&, bool&)
in file primitives/bool/boolIO.C at line 72.

Any idea?? Thanks for your support!!

[roger]

Hey Naish!

I have a problem trying to change PBICK to GAMG, because on "caceAgglomeration" I can only choose between yes/no instead of true/false, that's why it appears the error shown above.
I run checkMesh and it looks like everything is fine, maybe you can have a look and tell me what you think:

Create polyMesh for time = constant

Time = constant

Mesh stats
points: 86486
edges: 546809
faces: 895472
internal faces: 845132
cells: 435151
boundary patches: 6
point zones: 0
face zones: 0
cell zones: 1

Number of cells of each type:
hexahedra: 0
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 435151
polyhedra: 0

Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Topological cell zip-up check OK.
Face vertices OK.
Face-face connectivity OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface
CONJUNT 16156 8297 multiply connected surface (shared edge)
TRIMSRF 7074 3834 ok (not multiply connected)
IN2_SUP 20 24 ok (not multiply connected)
IN4_SUP 27 31 ok (not multiply connected)
OUT_SUP 55 44 ok (not multiply connected)
HABITACLE_SUP 27008 14038 ok (not multiply connected)
<<Writing 8 conflicting points to set nonManifoldPoints

Checking geometry...
Domain bounding box: (8016.26 0 11534.5) (11376.3 1533 12814.5)
Boundary openness (-1.52724e-16 5.65707e-17 4.50097e-16) OK.
Max cell openness = 1.63021e-16 OK.
Max aspect ratio = 49.6919 OK.
Minumum face area = 7.26316. Maximum face area = 2520.11. Face area magnitudes OK.
Min volume = 1.67593. Max volume = 39882.8. Total volume = 3.8029e+09. Cell volumes OK.
Mesh non-orthogonality Max: 84.7936 average: 14.6168
*Number of severely non-orthogonal faces: 278.
Non-orthogonality check OK.
<<Writing 278 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 2.42082 OK.
Min/max edge length = 2.83426 88.9473 OK.
All angles in faces OK.
All face flatness OK.

Mesh OK.

Thanks again!!

[tomf]

Hi Mayol,

I actually had these kind of problems as well and it seems like you may have made the same mistake I did. When I look at the dimensions of your domain bounding box, it says it is about 3000m long (x-direction) and same order of magnitude in the other dirctions. I guess that should be meters instead. If you imported the mesh from a fluent mesh, you should import it again, but scaled (fluentMeshToFoam -scale 0.001). This will bring down your Re by a factor of 1000, which will get better y+ values and physically more stable equations. Than your calculations should also be more stable. I am not sure about the non-orthogonality warning, but you could consult other topics for advise on that aspect.

Good luck!

[naish]

Hi Tom
As I mentioned I have the same problem, and your hit with the large domain sounds good. But what to do if the domanin has to be very large 100 - 200 meters ? I have about 10 mio cells a smaler version with 4 mio cells worked fine.

Thanks Naish

[matejfor]

Hi Mayol,

As Tom has pointed out, what I typically do, is to forget to downscale the mesh. It usually helps. Another thing, which might help you, looking at your mesh report and iteration progress is some limiting in discretisation schemes. Not sure what schemes you're using, but what helps me is cellLimiters for convective terms:
div(phi,U) Gauss linearUpwindV cellLimited leastSquares 1.0;

(this is laminar flow, so you might want to change linearUpwind with some other scheme.)

hope this helps...

[naish]

Hi
now my case is running, a I think that I have understand what was going wrong.

Please correct me if a I m wrong instead of downscaling I also can increase the viscosity, with the same effect (Re decreases). After a couple of timesteps when the model is more stable I can increase Re by changing viscosity.

Naish

[roger]

Hi!!

Tom, now it works quite better. I had a problem with the scale, thanks. Matej, I can show you my fvSchemes file, maybe you can suggest me some changes (although I haven't found "cellLimited" inside of any of the several options). Thanks guys, every step I make gets me closer to a reliable solution

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,R) Gauss upwind;
div(R) Gauss linear;
div(phi,nuTilda) Gauss upwind;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
default none;
laplacian(nuEff,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
interpolate(U) linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p;
}

[naish]

Hi mayol
Hope your case is running well, I plaied the whole weekend with different combination of parameters but not really successful in some cases it rises slower but it is sill diverging.

Have you alredy tried to run potentialFoam first?
With it you can create a startingfield, it makes the beginning of the calculation much more stable.

Please let me know if you found a solution, I am sill searching

[roger]

Hey Naish & cia.

I have tried several things but doesn't converge. Do you think it could work changing the values of "nNonOrthogonalCorrectors" , "pRefCell" , "pRefValue" ? And modifying relaxationFactors (and which values)? Thanks again.

[matejfor]

Mayol,
looking at your mesh orthogonality, it does not look very pretty, so it is a good idea to have number of nonorthogonal correctors around 2-4 I guess. Pressure reference in isothermal, incompressible case does not have a meaning.

Relaxation factors should be fine in default, I think the problem is elsewhere. What do you prescribe for pressure at the boundaries?
I think, you should check the epsilon once again. What is your inlet size to have this big epsilon?
I always use for the first hit this little script:

---------------------------------------------------
#!/bin/sh
echo ' ==========================================='
echo ' give me Uref, Turbulent Intensity and L'
echo ' I will give you TKE and epsilon for inlet'
echo ' ==========================================='
###
Uref=$1
ti=$2
L=$3
tke=$(echo "scale=4;1.5 *($Uref * $ti)^2"|bc -l )
epsilon=$(echo "scale=4;0.09^(3/4)*$tke^(3/2)/$L"|bc -l)
echo 'tke =' $tke
echo 'eps =' $epsilon
--------------------------------------------------
Does your gemoetry size match the inlet values?

For the schemes, I would suggest you to try this version:

divSchemes
{
default none;
div(phi,U) Gauss QUICKV cellLimited leastSquares 1.0;
div(phi,k) Gauss QUICK cellLimited leastSquares 1.0;
div(phi,epsilon) Gauss QUICK cellLimited leastSquares 1.0;
div(phi,R) Gauss linear;
div(R) Gauss linear;
div(phi,nuTilda) Gauss QUICK cellLimited leastSquares 1.0;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

Of course you do should try upwind, not a quick for a start and maybe even then go for some TVD scheme like VanLeer, but it is not that important in this moment.

hope this helps
good luck
matej

[eduardasilva]

Hello!

I am starting to implement a simple case in OpenFoam just to compare accuracy and stability with CFX. I imported a mesh from Gambit using fluentMeshToFoam. The mesh seems good after running checkMesh:

Overall domain bounding box (-1.22461e-18 -0.02 -0.0198552) (0.2 0.02 0.0198552)
Mesh (non-empty) directions (1 1 1)
Mesh (non-empty, non-wedge) dimensions 3
Boundary openness (-2.30502e-20 1.33809e-18 7.99073e-19) OK.
Max cell openness = 2.68493e-16 OK.
Max aspect ratio = 11.3925 OK.
Minumum face area = 2.2802e-06. Maximum face area = 2.92195e-05. Face area magnitudes OK.
Min volume = 1.10766e-08. Max volume = 1.21731e-07. Total volume = 0.000248888. Cell volumes OK.
Mesh non-orthogonality Max: 13.8762 average: 3.72945
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.381761 OK.

I am trying to run a laminar, steady and incompressible case in a straight tube which I was expecting to be easily feasible. I am using simpleFoam as solver and I cannot get the solution to converge. I checked all boundary conditions and tried different relaxation factors but nothing diferent happens. It seems to me that the residuals for U and p are going down quickly but for some reason the cumulative continuity errors are not.


Time = 252

DILUPBiCG: Solving for Ux, Initial residual = 4.23184e-07, Final residual = 4.23184e-07, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 5.03488e-05, Final residual = 5.03488e-05, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 1.3677e-05, Final residual = 1.3677e-05, No Iterations 0
DICPCG: Solving for p, Initial residual = 8.46934e-05, Final residual = 8.46934e-05, No Iterations 0
time step continuity errors : sum local = 0.000289333, global = 8.75706e-05, cumulative = 0.00393069
ExecutionTime = 5.61 s ClockTime = 5 s

What do you think I am doing wrong?

Thanks for your time.
Eduarda

[bastil]

Eduarda,

I dont see a reason why this should not be converged?

Regards

[eduardasilva]

Hello BastiL,

I am still learning OpenFoam so I don't exactly understand why this is happening. It is not running any iterations for U and p but the run doesn't converge. I am sending the outcome of two different times. Pressure is not varying and the time step continuity errors are not going down for some reason. Any ideas?

Time = 1025

DILUPBiCG: Solving for Ux, Initial residual = 1.66682e-14, Final residual = 1.66682e-14, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 3.82132e-14, Final residual = 3.82132e-14, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 3.94543e-14, Final residual = 3.94543e-14, No Iterations 0
DICPCG: Solving for p, Initial residual = 9.72482e-05, Final residual = 9.72482e-05, No Iterations 0
time step continuity errors : sum local = 0.000332199, global = 2.07969e-06, cumulative = -0.00607836
ExecutionTime = 20.03 s ClockTime = 20 s



Time = 3211

DILUPBiCG: Solving for Ux, Initial residual = 8.23629e-17, Final residual = 8.23629e-17, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 1.74586e-16, Final residual = 1.74586e-16, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 2.00532e-16, Final residual = 2.00532e-16, No Iterations 0
DICPCG: Solving for p, Initial residual = 9.72482e-05, Final residual = 9.72482e-05, No Iterations 0
time step continuity errors : sum local = 0.000332199, global = 2.07969e-06, cumulative = -0.00153215
ExecutionTime = 60.56 s ClockTime = 61 s