Numerical calculation interFoam never finishes
 

This might be a beginners question, but even though I've searched for an answer on this forum I didn't find any solution or tutorial explaining about this:

My problem is that I've imported a mesh from an stl. file, both using netgen (exported to fluent format) and tetgen. Converted to foam mesh with foam mesh conversion utilities. Nothing works:

Running with interfoam solver I can see that the time step, delta t is getting smaller and smaller (after perhaps 1-2 minutes: 1.5E-7 and getting smaller).

At the same time the max Courant number is growing with time. I suspect the reason for this behaviour is bad mesh and checkMesh also gives a some information about "skewness" and "severe non-orthogonality".

I don't know what either means, but I assume it's bad. Can anyone give an explanation and perhaps tell how to correct mesh errors if this is the problem?

With netgen I also tried to refine the mesh but nothing works.

Some output from checkMesh:
-------------------------

Create polyMesh for time = constant

Time = constant
Boundary definition OK.

Number of points: 1337
edges: 6656
faces: 9576
internal faces: 7448
cells: 4256
boundary patches: 3
point zones: 0
face zones: 0
cell zones: 0

Checking topology and geometry ...
Point usage check OK.

Upper triangular ordering OK.

Topological cell zip-up check OK.

Face vertices OK.

Face-face connectivity OK.


Basic topo ok ...

Checking patch topology for multiply connected surfaces ...

Patch Faces Points Surface
inlet 8 9 ok (not multiply connected)
atmosphere 52 45 ok (not multiply connected)
defaultFaces 2068 1054 ok (not multiply connected)


Patch topo ok ...
Topology check done.

Domain bounding box: min = (-4740.54 -5959.13 0) max = (4740.54 0 381) meters.

Checking geometry...
Boundary openness in x-direction = 5.52973e-10
Boundary openness in y-direction = 1.43336e-09
Boundary openness in z-direction = 2.80852e-09
Boundary closed (OK).
Max cell openness = 2.91038e-11 Max aspect ratio = 2.41902. All cells OK.

Minumum face area = 1493.91. Maximum face area = 126001. Face area magnitudes OK.

Min volume = 31884.7. Max volume = 7.52037e+06. Total volume = 7.10995e+09. Cell volumes OK.

Severe non-orthogonality for face 344 between cells 109 and 1298: Angle = 71.3136 deg.
Severe non-orthogonality for face 1437 between cells 505 and 4070: Angle = 71.4588 deg.
Severe non-orthogonality for face 2858 between cells 1204 and 1208: Angle = 71.327 deg.
Severe non-orthogonality for face 2865 between cells 1206 and 1209: Angle = 71.6745 deg.
Severe non-orthogonality for face 3041 between cells 1296 and 1301: Angle = 72.2016 deg.
Severe non-orthogonality for face 3052 between cells 1299 and 2726: Angle = 70.4236 deg.
Severe non-orthogonality for face 6739 between cells 3675 and 3679: Angle = 71.3936 deg.
Severe non-orthogonality for face 6741 between cells 3676 and 3681: Angle = 70.0194 deg.
Severe non-orthogonality for face 7208 between cells 4067 and 4071: Angle = 71.7692 deg.
Severe non-orthogonality for face 7209 between cells 4068 and 4073: Angle = 71.8987 deg.
Severe non-orthogonality for face 7210 between cells 4069 and 4072: Angle = 70.9372 deg.
Number of non-orthogonality errors: 0. Number of severely non-orthogonal faces: 11.
Mesh non-orthogonality Max: 72.2016 average: 33.2238
Non-orthogonality check OK.

Writing 11 non-orthogonal faces to set nonOrthoFaces

Face pyramids OK.

Max skewness = 173.084 percent. Face skewness OK.

Minumum edge length = 56.7614. Maximum edge length = 1058.37.

All angles in faces are convex or less than 10 degrees concave.

Geometry check done.

Number of cells by type:
hexahedra: 0
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 4256
polyhedra: 0
Number of regions: 1 (OK).
Mesh OK.

----
I was completely sure that Foam would be able to handle this kind of problem?

Ok, according to: http://www.c
 

Ok, according to: http://www.cfd-online.com/cgi-bin/Op...st=401#POST401 , then I'll have to play with the fvSchemes, fvSolution, maybe solver tolerances, non-orth correctors, number of PISO correctors and max Co number?

Any hints on doing this? This is very new to me...

Still I would like to hear about the above mesh (see output from checkMesh) - is it too bad or is it reasonable? Those angles don't tell me anything - I have absolutely no idea about "how bad" this mesh is :-)

Hi, The mesh is OK, nothing
 

Hi,

The mesh is OK, nothing wrong with it. The geometry seems a bit big: 4 kilometers - did you forget to scale it by any chance. Orthogonality is a bit big but not excessively so and FOAM will deal with it without trouble. If you want ot check the mesh a bit more, try doing a potential flow solution or trying an incompressible flow solver on the same geometry. My guess is that you've messed up the initial or boundary conditions or chosen inappropriate discretisation parameters.

Good luck,

Hrv

Hi Hrvoje, It seems like yo
 

Hi Hrvoje,

It seems like you must have a lot of experience here - I think I remember that I read something from a report about error estimation by you.

1a) Yep, I must scale it :-)
1b) I also had problems with the icoFoam solver, but I'll try to experiment again...

2) How can you tell that the mesh is ok and that the mesh doesn't have to be repaired or rebuild/refined? How coarse can I make the mesh?

3) Is there any (online) documentation available (somewhere) about orthogonality, skewness or something appropriate within this topic for beginners?

Hi, Martin, Can you explain
 

Hi, Martin,

Can you explain how the case was initialized? For interFoam, you need to initialize few layer of cells to the same liquid as the inlet liquid (for example, water coming into an empty tube, you must initialize few cells close to the inlet to water - gamma = 1). This used to be a big problem prior to version 1.2. But, I am still having the same problem in version 1.2.

Also, I never had good luck with tet elements when doing free surface flow. Try changing to Hex elements if you can.

Pei

Hi Pei, Why do you need to
 

Hi Pei,

Why do you need to initialize the cells near the entrance with gamma = 1?

I have done simulations with the entrance of a fluid and usually setting the boundary face (gamma = 1) is enough. What type of flow are simulating?

Billy.

Hi Pei (and Billy). Thanks
 

Hi Pei (and Billy).

Thanks for the comments. I think (at least I hope) that the gamma value is initialized correctly since I've got at test-case running which is a square box made from a modified blockMeshdict from the icoFoam / cavity-tutorial (this works correctly with interFoam).

So the only difference now, is that I'm trying to convert the box-case (which works with mesh made by blockMesh) and getting the stl-file case working...

Various parameters:

--------------------
/0/U:
--------------------
internalField uniform (0 0 0);

boundaryField
{
inlet
{
type fixedValue;
value uniform (0 -5 0);
}
atmosphere
{
type fixedValue;
value uniform (0 0 0);
}
defaultFaces
{
type fixedValue;
value uniform (0 0 0);
}

--------------------
/0/gamma
--------------------
internalField uniform 0;

boundaryField
{
inlet
{
type fixedValue;
value uniform 1;
}
atmosphere
{
type zeroGradient;
// type fixedValue;
// value uniform 0;
}
defaultFaces
{
type zeroGradient;
}
}

--------------------
/0/pd
--------------------
internalField uniform 0;

boundaryField
{
inlet
{
type zeroGradient;
}
atmosphere
{
type fixedValue;
value uniform 0;
}
defaultFaces
{
type zeroGradient;
}
}
--------------------

Also a big concern of mine is that I don't know much about discretizing methods but I'll look in the users guide again.

--------------------
/system/fvSchemes
--------------------
ddtSchemes
{
default Euler;
}

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

divSchemes
{
div(rho*phi,U) Gauss upwind;
div(phi,gamma) Gauss limitedLinear01 1;
div(phirb,gamma) Gauss limitedLinear01 1;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
pd;
pcorr;
gamma;
}

--------------------
/system/fvSolution
--------------------
solvers
{
pcorr ICCG 1e-10 0;
pd ICCG 1e-10 0;
gamma BICCG 1e-08 0;
U BICCG 1e-06 0;
}

PISO
{
momentumPredictor no;
nCorrectors 3;
nNonOrthogonalCorrectors 1;
nGammaCorr 1;
nGammaSubCycles 4;
cGamma 2;
}

I've tried to raise the nNonOrth..correctors to around 20 I think. Actually I don't know where to start changing values :-)

Hi, Martin and Billy, This
 

Hi, Martin and Billy,

This is the problem description of the problem I worked on before:

a simple tube (and/or a rectangle channel), initially filled with air (gamma = 0). At t=0, water (gamma =1) flows into the domain. In my problem, both surface tension and wall contact angle were important. Prior to version 1.2, delta t quickly dropped to 10e-8, 10e-9 or smaller and may diverge. The only cure was to initialize few layer of cells to gamma =1 near the inlet. Henry made improvements to interFoam in version 1.2. However, I am still having problem if I do not initialize the cells near the inlet to the inlet fluid.

Martin,

for divSchemes, I usually use:
div(rho*phi,U) Gauss limitedLinearV01 1;
div(phi,gamma) Gauss limitedLinear01 1;
div(phirb,gamma) Gauss gammaCompression 1;

But, I do not think this is the cause of your problem though.

If you do not mind, edit setFieldDict -> pick a box that will cover 2 -3 layers of cells next to inlet and set them to gamma = 1. This will modify 0/gamma. Plot the flow field to check if gamma field is correct using paraFoam. Then run the case to see you still have the same problem. If this still fails, I am willing to take a quick look of you case if you do not mind (although I cannot guaranttee you I can find the root cause of your problem).

pei

Hi Pei, I tried both your s
 

Hi Pei,

I tried both your suggestions. No luck.

AFAIR I can't post the case here since it's exceeding the upload limit (the tar.gz file is roughly 100 kb), so I've sent you an e-mail. Let me know if you didn't get it.

If anyone else wants to look at the case, let me know.

Also: Perhaps Hrvoje (or somebody) else will explain how to distinguish a bad mesh from a good mesh with checkMesh, since I have tried other imports which have had far more "severe non-orthog/skewness" than here.

The warnings in primitiveMesh
 

The warnings in primitiveMesh are just that. They might cause problems, especially with less damped numerics, or might not. It depends on the algorithm (interFoam, simpleFoam etc) and settings (divSchemes, turbulenceModel etc.) you are trying to run it with.

The errors though (non-orthogonality errors, wrong oriented faces) will cause problems and should be fixed in the mesh.

In general checkMesh is the first thing to run and will give you an indication of mesh related problems. As Hrvoje says run a simple solver (e.g. potentialFoam) on it to see if solving a basic elliptic equation (the basis of almost all codes) is possible. Then try running your code, first with stable numerics.