[nicholas.jones]

All,

I am trying to complete the meshing of a cylinder tightly packed with particles. On few occasions, I get a stable mesh, which has been attached.

However, if I were to use the same snappyHexMesh dict, while doubling the number of cells in the blockMeshDict (a=b=32, c=48), the simulation crashes before a single iteration, giving the following error (when not run in parallel, to keep it smaller):

Starting time loop

Time = 0.0001

Courant Number mean: 6.80645e-05 max: 0.0282073
smoothSolver: Solving for Ux, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver: Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver: Solving for Uz, Initial residual = 1, Final residual = 9.89261e-06, No Iterations 39
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam:ICPreconditioner::calcReciprocalD(Foam::Fie ld<double>&, Foam::lduMatrix const&) at ??:?
#4 Foam:ICPreconditioner:ICPreconditioner(Foam::l duMatrix::solver const&, Foam::dictionary const&) at ??:?
#5 Foam::lduMatrix:reconditioner::addsymMatrixConst ructorToTable<Foam:ICPreconditioner>::New(Foam:: lduMatrix::solver const&, Foam::dictionary const&) at ??:?
#6 Foam::lduMatrix:reconditioner::New(Foam::lduMatr ix::solver const&, Foam::dictionary const&) at ??:?
#7 Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#8 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ??:?
#9 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#10
at ??:?
#11
at ??:?
#12 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#13
at ??:?
Floating point exception (core dumped)

My blockmeshdict/snappyhexmeshdict are below:

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

a 16; // Number of elements in x direction
b 16; // Number of elements in y direction
c 24; // Number of elements in z direction
xi -0.01; // Minimum x coordinate
xf 0.07; // Maximum x coordinate
yi -0.01; // Minimum y coordinate
yf 0.07; // Maximum y coordinate
zi -0.01; // Minimum z coordinate
zf 0.11; // Maximum z coordinate


vertices
(
($xi $yi $zi) //0
($xf $yi $zi) //1
($xf $yf $zi) //2
($xi $yf $zi) //3
($xi $yi $zf) //4
($xf $yi $zf) //5
($xf $yf $zf) //6
($xi $yf $zf) //7
);


blocks
(
hex (0 1 2 3 4 5 6 7) ($a $b $c) simpleGrading (1 1 1)
);

edges
(
);

patches
(
);

mergePatchPairs
(
);

// ************************************************** *********************** //

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object autoHexMeshDict;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

castellatedMesh true;
snap true;
addLayers false;

geometry
{
Inlet.stl {type triSurfaceMesh; name Inlet;}
Outlet.stl {type triSurfaceMesh; name Outlet;}
Balls.stl {type triSurfaceMesh; name Balls;}
Walls.stl {type triSurfaceMesh; name Walls;}
};


castellatedMeshControls
{
maxLocalCells 10000000;
maxGlobalCells 25000000;
minRefinementCells 100;
nCellsBetweenLevels 2;

features
(
// {file "Inlet.eMesh"; level 2;}
// {file "Outlet.eMesh"; level 2;}
// {file "Balls.eMesh"; level 2;}
// {file "Walls.eMesh"; level 2;}
);

refinementSurfaces
{
Inlet {level (1 3);}
Outlet {level (1 3);}
Balls {level (2 5);}
Walls {level (2 5);}
}

resolveFeatureAngle 60;

refinementRegions
{
}

locationInMesh (0.01519 0.01008 0.0982);

allowFreeStandingZoneFaces true;
}


snapControls
{
nSmoothPatch 2;
tolerance 4;
nSolveIter 60;
nRelaxIter 20;
}





addLayersControls
{
relativeSizes true;

layers
{
Balls
{
nSurfaceLayers 6;
}
}

expansionRatio 1.1;
finalLayerThickness 0.7;
minThickness 0.1;
nGrow 1;


// Advanced settings
featureAngle 60;
nRelaxIter 5;
nSmoothSurfaceNormals 1;
nSmoothNormals 3;
nSmoothThickness 10;
maxFaceThicknessRatio 0.9;
maxThicknessToMedialRatio 0.9;
minMedianAxisAngle 130;
nBufferCellsNoExtrude 0;
nLayerIter 50;
nRelaxedIter 20;
}



meshQualityControls
{

maxNonOrtho 65;
maxBoundarySkewness 20;
maxInternalSkewness 4;
maxConcave 80;
minFlatness 0.5;
minVol 1e-13;
minTetQuality 1e-30;
minArea -1;
minTwist 0.05;
minDeterminant 0.001;
minFaceWeight 0.05;
minVolRatio 0.01;
minTriangleTwist -1;

// Advanced

nSmoothScale 4;
errorReduction 0.75;


relaxed
{
maxNonOrtho 75;
}
}


// Advanced

// Flags for optional output
// 0 : only write final meshes
// 1 : write intermediate meshes
// 2 : write volScalarField with cellLevel for postprocessing
// 4 : write current intersections as .obj files
debug 0;

mergeTolerance 1E-6;


// ************************************************** *********************** //

I have tried running unset FOAM_SIGFPE, although the simulation still crashes, citing a nan input.

I can typically generate meshes for fewer particles, but once I increase beyond 20, the majority of my meshes crash on first iteration. I have tried toying with snappyHexMeshDict and blockMeshDict, including fine grids to start, with few refinement iterations, or coarse grids to start, with many refinement iterations.

Any ideas how to solve this issue? .stl files are available at the following link:
http://s000.tinyupload.com/index.php...75582890393339

controlDict is below:

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application icoFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 0.5;

deltaT 0.0001;

writeControl timeStep;

writeInterval 10;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;


// ************************************************** *********************** //

Any help/suggestions are appreciated!

[danielpiaget]

Hello Nicolas,

The crash seems to be coming when you try to solve the matrix to solve the pressure equation. The error seems to indicate that the solver is excepting an integer but is getting something else. The first step would be to check the mesh quality. Regarding the mesh, you can check the quality of the mesh with the following command:

checkMesh -allGeometry -allTopology

And compare between the two meshes (20 particles and more that 20 particles). Compare the logs.
After that check your boundary conditions. Are the boundary conditions dependent on the number of particles ?

Good luck,

Daniel

[nicholas.jones]

Thanks for the tip!

I have been doing a standard checkMesh, the added parameters should be a useful tool. I will be comparing the stable and unstable results tomorrow. If you are interested/willing, I have attached the boundary conditions and checkMesh results. In both tests, boundary conditions are unchanged, and are listed below. I have tried to keep things simple.

FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
Inlet
{
type zeroGradient;
}

Outlet
{
type fixedValue;
value uniform 0;
}

Walls
{
type zeroGradient;
}

Balls
{
type zeroGradient;
}

defaultFaces
{
type empty;
}
}

// ************************************************** *********************** //

FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
Inlet
{
type fixedValue;
value uniform (0 0 0.1);
}

Outlet
{
type zeroGradient;
}

Balls
{
type fixedValue;
value uniform (0 0 0);
}

Walls
{
type fixedValue;
value uniform (0 0 0);
}

defaultFaces
{
type empty;
}
}

// ************************************************** *********************** //

The checkMesh -allGeometry -allTopology for the stable mesh is shown below, followed by the unstable mesh. The only thing that changes is the number of cells in the blockMeshDict (doubles in all directions)

Enabling all (cell, face, edge, point) topology checks.

Enabling all geometry checks.

Time = 0

Mesh stats
points: 118931
faces: 285221
internal faces: 249936
cells: 86998
faces per cell: 6.15137
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 62326
prisms: 4661
wedges: 0
pyramids: 0
tet wedges: 406
tetrahedra: 15
polyhedra: 19590
Breakdown of polyhedra by number of faces:
faces number of cells
4 4834
5 4326
6 3510
7 270
8 296
9 3316
10 184
11 179
12 1477
13 119
14 92
15 680
16 16
17 11
18 269
21 11

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Topological cell zip-up check OK.
Face-face connectivity OK.
<<Writing 2 cells with two non-boundary faces to set twoInternalFacesCells
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology Bounding box
defaultFaces 0 0 ok (empty)
Inlet 1590 1923 ok (non-closed singly connected) (1.50183e-05 1.50183e-05 0) (0.0634863 0.0634863 0.000367069)
Outlet 1590 1923 ok (non-closed singly connected) (1.50183e-05 1.50183e-05 0.0996329) (0.0634863 0.0634863 0.1)
Balls 25489 33163 ok (non-closed singly connected) (5.99532e-05 4.09208e-05 0.0100567) (0.0634428 0.0634402 0.0641586)
Walls 6616 8453 ok (non-closed singly connected) (1.29906e-05 1.29906e-05 0) (0.063494 0.063494 0.1)

Checking geometry...
Overall domain bounding box (1.29906e-05 1.29906e-05 0) (0.063494 0.063494 0.1)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-1.27511e-16 6.44614e-18 -3.47273e-16) OK.
Max cell openness = 2.96904e-16 OK.
Max aspect ratio = 8.85437 OK.
Minimum face area = 5.13905e-08. Maximum face area = 0.000104055. Face area magnitudes OK.
Min volume = 2.05956e-11. Max volume = 1.04459e-06. Total volume = 0.000225262. Cell volumes OK.
Mesh non-orthogonality Max: 64.7458 average: 13.4013
Non-orthogonality check OK.
Face pyramids OK.
***Max skewness = 5.05559, 1 highly skew faces detected which may impair the quality of the results
<<Writing 1 skew faces to set skewFaces
Coupled point location match (average 0) OK.
Face tets OK.
Min/max edge length = 9.57176e-05 0.0103234 OK.
*There are 1782 faces with concave angles between consecutive edges. Max concave angle = 74.4631 degrees.
<<Writing 1782 faces with concave angles to set concaveFaces
Face flatness (1 = flat, 0 = butterfly) : min = 0.741794 average = 0.99559
*There are 4 faces with ratio between projected and actual area < 0.8
Minimum ratio (minimum flatness, maximum warpage) = 0.741794
<<Writing 4 warped faces to set warpedFaces
Cell determinant (wellposedness) : minimum: 0 average: 15.9207
***Cells with small determinant (< 0.001) found, number of cells: 2
<<Writing 2 under-determined cells to set underdeterminedCells
***Concave cells (using face planes) found, number of cells: 5782
<<Writing 5782 concave cells to set concaveCells
Face interpolation weight : minimum: 0.0604884 average: 0.442852
Face interpolation weight check OK.
Face volume ratio : minimum: 0.011866 average: 0.707596
Face volume ratio check OK.

Failed 3 mesh checks.

End

Enabling all (cell, face, edge, point) topology checks.

Enabling all geometry checks.

Time = 0

Mesh stats
points: 273198
faces: 660154
internal faces: 580315
cells: 201648
faces per cell: 6.15166
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 143290
prisms: 13162
wedges: 0
pyramids: 0
tet wedges: 695
tetrahedra: 5
polyhedra: 44496
Breakdown of polyhedra by number of faces:
faces number of cells
4 10372
5 8823
6 9153
7 718
8 666
9 7889
10 368
11 284
12 3654
13 216
14 180
15 1595
16 37
17 31
18 497
21 13

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Topological cell zip-up check OK.
Face-face connectivity OK.
<<Writing 4 cells with zero or one non-boundary face to set oneInternalFaceCells
*Number of regions: 5
The mesh has multiple regions which are not connected by any face.
<<Writing region information to "0/cellToRegion"
<<Writing region 0 with 201644 cells to cellSet region0
<<Writing region 1 with 1 cells to cellSet region1
<<Writing region 2 with 1 cells to cellSet region2
<<Writing region 3 with 1 cells to cellSet region3
<<Writing region 4 with 1 cells to cellSet region4

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology Bounding box
defaultFaces 0 0 ok (empty)
Inlet 3763 4426 ok (non-closed singly connected) (1.36051e-05 1.36051e-05 0) (0.0634248 0.0634248 0.000182237)
Outlet 3763 4426 ok (non-closed singly connected) (1.36051e-05 1.36051e-05 0.0998178) (0.0634248 0.0634248 0.1)
Balls 55595 72259 ok (non-closed singly connected) (0.000624071 0.000537626 0.0100199) (0.0633189 0.063386 0.0642066)
Walls 16718 18977 ok (non-closed singly connected) (5.69448e-06 1.25903e-05 0) (0.0634887 0.0634887 0.1)

Checking geometry...
Overall domain bounding box (5.69448e-06 1.25903e-05 0) (0.0634887 0.0634887 0.1)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-1.83914e-16 -5.00593e-17 1.42976e-16) OK.
Max cell openness = 3.15507e-16 OK.
Max aspect ratio = 5.72821 OK.
Minimum face area = 7.623e-09. Maximum face area = 2.59008e-05. Face area magnitudes OK.
Min volume = 7.62194e-13. Max volume = 1.29755e-07. Total volume = 0.000229602. Cell volumes OK.
Mesh non-orthogonality Max: 59.9526 average: 13.2447
Non-orthogonality check OK.
Face pyramids OK.
***Max skewness = 5.41298, 2 highly skew faces detected which may impair the quality of the results
<<Writing 2 skew faces to set skewFaces
Coupled point location match (average 0) OK.
Face tets OK.
Min/max edge length = 6.6954e-05 0.00512291 OK.
*There are 3518 faces with concave angles between consecutive edges. Max concave angle = 68.4873 degrees.
<<Writing 3518 faces with concave angles to set concaveFaces
Face flatness (1 = flat, 0 = butterfly) : min = 0.722171 average = 0.996185
*There are 17 faces with ratio between projected and actual area < 0.8
Minimum ratio (minimum flatness, maximum warpage) = 0.722171
<<Writing 17 warped faces to set warpedFaces
Cell determinant (wellposedness) : minimum: 0 average: 16.0598
***Cells with small determinant (< 0.001) found, number of cells: 4
<<Writing 4 under-determined cells to set underdeterminedCells
***Concave cells (using face planes) found, number of cells: 13524
<<Writing 13524 concave cells to set concaveCells
Face interpolation weight : minimum: 0.0946927 average: 0.444179
Face interpolation weight check OK.
Face volume ratio : minimum: 0.0189728 average: 0.713364
Face volume ratio check OK.

Failed 3 mesh checks.

End

[danielpiaget]

Good morning Nicolas,

What version of openFOAM are you using ? If you are using version 2.3.0, your dictionary file seems to have missing parameters.....Check the tutorials in your folder

openFOAM/openFOAM-2.3.0/tutorials/mesh/snappyHexMesh/

I have notice that the number of regions is different between the two mesh.You have one region for the first mesh and 5 regions in the second mesh ? Also, there seems to have cells that are not connected. See the message from the second mesh.

"The mesh has multiple regions which are not connected by any face."

Daniel

[nicholas.jones]

Hey Daniel,

Thanks for the analysis! The missing entries have been added to the snappyHexMeshDict.

Do you have any idea why there are multiple regions being created? How can I control this beyond trial and error?

Edit: solved by using splitMeshRegions -largestOnly. However, any insight you can provide as to why this is happening is appreciated!

[danielpiaget]

Hello,

One cannot be sure, but I would speculate that when snappyHexMesh tried
to determine which part of the initial mesh (blockMes) was to be kept, some sections of the mesh where not connect, so it created different regions instead of one region (created by default if not specified).

Does your simulation work now ?

Thanks,

Daniel

[nicholas.jones]

Hey Daniel,

The splitMeshRegions -largestOnly worked, the resulting polymesh is used with no further modifications to any dict/boundary files.

Thanks for the help!