cylinder case
In this tutorial is presented a quite simple case of a cylinder mesh made with snappyHexMesh utility. You can realize it in many different ways (with open or close STL, with big, small or right dimensions of blockMesh and so on). All this different ways are presented below.
First of all a picture that describes the case settings (each aspect will be explained later). You can copy the files contained in motorBike tutorial contained in $FOAM_TUTORIALS/incompressible/simpleFoam/motorBike and modify the case.
At this link you can download the whole case. It may contain errors so check all the files while reading this guide and then try to modify some parameters to see what changes in results.
STL
In this tutorial, for educational purposes, two different STL are created:
- Close: the close STL represent the pipe with all surfaces meshed with triangles
- Open: the open STL represent the pipe only by the lateral surface (inlet and outlet faces are missing)
We will use both in the tutorial with different blockMesh.
The dimension of the cylinder are:
- base: with radius of 0.5
- height of 6 in z direction
In STL files we set just one layer name, including all surfaces, called suplat
Directory system
See cube case tutorial for controlDict, fvSchemes, fvSolution and decomposeParDict.
blockMeshDict
First of all we want to create a mesh of cubes (a classic with snappy), so also the blockMesh has to be made of cubes (we will create a mesh of elongated elements after this case).
We analyze three type of blockMesh:
Big blockMesh
Is a blockMesh greater than the STL geometry especially in axial direction (z)
Just the important part of blockMeshDict is shown
convertToMeters 1;
a 20;
b 20;
c 80;
yi -1;
yf 1;
zi -1;
zf 7;
xi -1;
xf 1;
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)
);
Right blockMesh
Is a blockMesh long as the cylinder (in axial direction)
To change the blockMesh just change the bold parameters
a 20;
b 20;
c 60; //to maintain the cubes dimensions unchanged
yi -1;
yf 1;
zi 0;
zf 6;
xi -1;
xf 1;
Small blockMesh
Is a blockMesh smaller than the cylinder (in axial direction)
To change the blockMesh just change the bold parameters
a 20;
b 20;
c 40; //to maintain the cubes dimensions unchanged
yi -1;
yf 1;
zi 1;
zf 4;
xi -1;
xf 1;
The results is shown in this pictures, all the blockMesh little cubes have the same dimension, the total number of these cubes (especially the number in axial direction) changes from one to another blockMesh
Snap problems
As usual we encounter some problems in snappyHexMesh, especially in the snap part, but with some arrangements regarding the STL and blockMesh choice we can obtain a good result.
We want to analyze the snap problems, so deactivate the layer part in activation rows.
castellatedMesh true;
snap true;
addLayers false;
and try to run snappyHexMesh with these (semi)default values:
- we define the surface name as cylinder and recall the STL region suplat, named mysuplat;
- we set a single refinement level with 6 cells in it for cylinder surface;
- we set a quite high value of featureAngle to avoid high refinement zones;
- we select an internal point to obtain internal mesh
- we set quite improved values of snap controls (improved compared to default values)
geometry
{
cylinder_close.stl; //or cylinder_open.stl;
{
type triSurfaceMesh;
name cylinder;
regions
{
suplat
{
name mysuplat;
}
}
}
};
castellatedMeshControls
{
maxLocalCells 1500000;
maxGlobalCells 2000000;
minRefinementCells 10
nCellsBetweenLevels 6;
features
(
);
refinementSurfaces
{
cylinder
{
level (1 1);
}
}
resolveFeatureAngle 95;
refinementRegions
{
}
locationInMesh (0 0 3);
}
snapControls
{
nSmoothPatch 8;
tolerance 2.0;
nSolveIter 100;
nRelaxIter 10;
}
Close STL
At first we try to run snappyHexMesh with the close STL with different blockMesh, the snappyHexMeshDict file doesn't change, only change the files into constant/polyMesh directory (the blockMesh files)
Big blockMesh
Two problems:
- snappy doesn't recognize the sharp edges and can't snap well the mesh points to the surface
- we get a mesh thickening also towards the inlet and outlet faces, while often, in pipe simulations, you just want a lateral walls thickening.
The first problem is very serious and difficult to resolve (by now I don't know what to do), the second can be solved defining two regions in STL file (one for lateral walls and one for inlet and outlet faces), setting a 0 refinement level for the whole surface and a single refinement for the wall region.
refinementSurfaces
{
cylinder
{
level (0 0);
regions
{
mysuplat
{
level (1 1);
}
}
}
}
The result is shown in this picture
You can try to improve snap controls (decreasing tolerance and increasing all other options) but the results doesn't change
Right blockMesh
With the right blockMesh the mesh quality improves on sharp edges, but we still have some imperfections (honestly they could be attributed to the bad mesh quality in STL file near the edge).
We still have the wrong thickening problem, but it is not a big deal and we can solve it. This is a slice that shows the thickening problem.
Small blockMesh
With a blockMesh smaller than the cylinder we resolve the sharp edges problems and the wrong thickening problems, but (obviously) we obtain a smaller mesh. So if you want to simulate the flow in a 30 meters pipe, you have to make an STL file of 40 meters and then the blockMesh will cut your geometry. It is not a disaster, but not so convenient.
By the way the results are great (but the cylinder is smaller)
Conclusion
The best option with a close STL is to use a small blockMesh (smaller than the geometry), but you have to draw a bigger STL that will be cut by the blockMesh in snappyHexMesh process. Also the right blockMesh could work well, but you have to draw a very refined STL especially near sharp edges.
Open STL
With a simple geometry like our cylinder, we can try to run snappyHexMesh using an open STL, to avoid the previous (close STL) problems. This STL has no inlet and outlet faces, but only the lateral surface (called suplat in STL file)
Big blockMesh
We can't use this blockMesh because it is greater than the geometry and the STL surfaces can't subdivide it in two distinct regions, so we can't obtain the internal mesh.
Right blockMesh
The mesh obtained with a right blockMesh is perfect. No evident problems are found: sharp edges are fine, the thickening (a problem in close STL + right blockMesh) is correct and the cylinder has the proper length.
Small blockMesh
The mesh is very similar to the one obtained with close STL + small blockMesh, with the same pro and con.
Conclusion
We have found the perfect mesh, with no problem, is the open STL + right blockMesh. Now we can analyze the layer addition.
Layers
After finding the perfect combination to obtain a good mesh we can proceed adding layers.
We activate the addLayer controls in activation rows and set the snappyHexMeshDict file:
- we decide to add three layers on mysuplat sub-region with relative sizes
- final layer thickness is 0.7 times the near cells thickness, with an expansion ratio of 1.1
- we increase the maxFaceThicknessRatio and the maxThicknessToMedialRatio to help layer addition also in presence of highly warped cells
addLayersControls
{
relativeSizes true;
layers
{
mysuplat
{
nSurfaceLayers 3;
}
}
expansionRatio 1.1;
finalLayerThickness 0.7;
minThickness 0.1;
nGrow 1;
// Advanced settings
featureAngle 110;
nRelaxIter 3;
nSmoothSurfaceNormals 1;
nSmoothNormals 3;
nSmoothThickness 10;
maxFaceThicknessRatio 0.9;
maxThicknessToMedialRatio 0.9;
minMedianAxisAngle 130;
nBufferCellsNoExtrude 0;
nLayerIter 50;
}
Close STL
We analyze now the layer addition on mesh obtained with close STL
Big blockMesh
Near the sharp edges (badly meshed) the layers are unacceptable. This problem stems from the bad mesh previously created.
Right blockMesh
Also in this case the layer addition doesn't run very well, especially in proximity of inlet outlet faces (because of the bad previously mesh)
Small blockMesh
We obtain good layers, but the cylinder is smaller than the STL geometry
Conclusion
"With great mesh comes great layers." [Uncle Snappy cit.]
Open STL
We analyze the layer addition on mesh obtained with open STL. We don't consider big blockMesh because we didn't obtain any mesh with open STL + big blockMesh
Right blockMesh
Perfect layer addition. Nothing to say
Small blockMesh
Perfect layer addition, but smaller cylinder (have to consider that).
Conclusion
"With great mesh comes great layers." [Uncle Snappy cit.]
Axial case (open STL)
A common problem of the previous mesh is the high number of cells. Our case is very simple and grossly meshed, but in real cases, using cubes, we will obtain a very high number of elements, hence slow simulations.
To solve this problem we can create a blockMesh with elongated elements (parallelepipeds), snappyHexMesh will halve all dimensions obtaining four similar parallelepipeds by the first one.
This is a common solution in pipe flow analysis, but with too elongated elements we may obtain elements with nonOrthoFaces.
In order to create a case with elongated elements (called axial case) we use the best couple obtained in previous simulations: open STL + right blockMesh with two different blockMesh
- Axial : slightly elongated elements
- Axial extreme : heavily elongated elements, minor number of elements, but presence of nonOrthoFaces
blockMeshDict
You can recycle the old blockMeshDict and just change the head parameters.
Axial
a 20;
b 20;
c 15;
yi -1;
yf 1;
zi 0;
zf 6;
xi -1;
xf 1;
Extreme axial
a 20;
b 20;
c 7;
yi -1;
yf 1;
zi 0;
zf 6;
xi -1;
xf 1;
snappyHexMeshDict
The snappyHexMeshDict file is similar to the one used for previous case. Just change a control parameter:
- we disable the control on nonOrthoFaces because if snappyHexMesh find some of them it tries to repair them and fail, so doesn't add layers (5.4.7 OpenFOAM UG). We ignore the control on them but we have to keep in mind that we could have nonOrthoFaces if we encounter problems solving a case with this mesh.
geometry
{
cylinder_wiki_open.stl;
{
type triSurfaceMesh;
name cylinder;
{
suplat
{
name mysuplat;
}
}
}
};
castellatedMeshControls
{
maxLocalCells 1500000;
maxGlobalCells 2000000;
minRefinementCells 10
nCellsBetweenLevels 6;
features
(
);
refinementSurfaces
{
cylinder
{
level (1 1);
}
}
resolveFeatureAngle 95;
refinementRegions
{
}
locationInMesh (0 0 3);
}
snapControls
{
nSmoothPatch 8;
tolerance 2.0;
nSolveIter 100;
nRelaxIter 10;
}
addLayersControls
{
relativeSizes true;
layers
{
mysuplat
{
nSurfaceLayers 3;
}
}
expansionRatio 1.1;
finalLayerThickness 0.7;
minThickness 0.1;
nGrow 1;
featureAngle 110;
nRelaxIter 3;
nSmoothSurfaceNormals 1;
nSmoothNormals 3;
nSmoothThickness 10;
maxFaceThicknessRatio 0.9;
maxThicknessToMedialRatio 0.9;
minMedianAxisAngle 130;
nBufferCellsNoExtrude 0;
nLayerIter 50;
}
meshQualityControls
{
//- Maximum non-orthogonality allowed. Set to 180 to disable.
maxNonOrtho 180; //65;
...
...
}
Axial
The result is very nice, the layer addition worked properly and we managed to keep a lower number of elements
In this slice you can see the parallelepipeds elements split
Extreme axial
In this case the number of cells decreases , the mesh is very handy but we risk getting nonOrthoFaces. Moreover the number of cells between levels isn't six (as we set) but snappy mess a little in boundary refinement zones. We should check the mesh.
checkMesh
The checkMesh is an useful utility which finds any errors in our mesh. You can launch it by typing checkMesh in a command line, you can also add the -latestTime options to check just the latest timestep created
Cube case (open STL)
We should control also the previous mesh, which we called ''perfect'', to see if it is ready for a simulation process or not. just type checkMesh -latestTime in open STL + right blockMesh folder and read the output. It should be like this:
Create time
Create polyMesh for time = 3
Time = 3
Mesh stats
points: 56933
faces: 162928
internal faces: 155312
cells: 53040
boundary patches: 2
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 49440
prisms: 2880
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 720
Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface topology
defaultFaces 896 946 ok (non-closed singly connected)
mysuplat 6720 6776 ok (non-closed singly connected)
Checking geometry...
Overall domain bounding box (-0.499997 -0.499997 0) (0.499997 0.499997 6)
Mesh (non-empty, non-wedge) directions (1 1 0)
Mesh (non-empty) directions (1 1 0)
***Number of edges not aligned with or perpendicular to non-empty directions: 34175
<<Writing 35158 points on non-aligned edges to set nonAlignedEdges
Boundary openness (-2.14369e-16 -1.32609e-15 5.96795e-18) OK.
Max cell openness = 3.09416e-16 OK.
Max aspect ratio = 3.44995 OK.
Minumum face area = 0.000775603. Maximum face area = 0.01. Face area magnitudes OK.
Min volume = 3.87981e-05. Max volume = 0.000748594. Total volume = 4.70017. Cell volumes OK.
Mesh non-orthogonality Max: 31.4371 average: 5.60554
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.371812 OK.
Failed 1 mesh checks.
End
The checkMesh failed, because of the presence of nonAlignedEdges. This is because snappy creates some default patch (called defaultFaces, you can see it in output) and set to empty.
boundary and autoPatch
We have two ways to fix this problem.
- defaultFaces is a patch containing inlet and outlet faces. If we don't care to distinguish these faces we can simply edit 3/polyMesh/boundary file and set defaultFaces as a patch, not empty
2
(
defaultFaces
{
type patch; // this was empty before our change
nFaces 896
startFace 155312;
}
mysuplat
{
type wall;
nFaces 6720;
startFace 156208;
}
)
If we run checkMesh now it should worked well.
- If we care to distinguish from inlet and outlet faces (as usual) we should run autoPatch utility on third step (latest). First of all modify system/controlDict file and set
startFrom latestTime; // It was startTime before our change
so autoPatch will work on mesh contained in 3 folder and will create a 4 folder with patched mesh. To run autoPatch type autoPatch 75 in a command line. 75 is the featureAngle used for patching (play on this if autoPatch doesn't work properly). The utility will create a new boundary file in 4/polyMesh folder:
5
(
defaultFaces
{
type empty;
nFaces 0;
startFace 155312;
}
mysuplat
{
type wall;
nFaces 0;
startFace 155312;
}
auto0
{
type patch;
nFaces 448;
startFace 155312;
}
auto1
{
type patch;
nFaces 448;
startFace 155760;
}
auto2
{
type patch;
nFaces 6720;
startFace 156208;
}
)
autoPatch keeps old patches and add the new; auto0 and auto1 are inlet and outlet, while auto2 represents the walls. We should fix it
3
(
inlet
{
type patch;
nFaces 448;
startFace 155312;
}
outlet
{
type patch;
nFaces 448;
startFace 155760;
}
wall
{
type wall;
nFaces 6720;
startFace 156208;
}
)
You could set wall type for wall patch if you will use wall functions on this surfaces.
Now the checkMesh should run properly.
Axial case (open STL)
If we run checkMesh on axial case, we obtain an output very similar to the open STL case; checkMesh warns us of the presence of nonAlignedEdges (is an intrinsic problem of snappyHexMesh) we could fix this as shown before.
If we try to run checkMesh on axial_extreme case the output is:
Create time
Create polyMesh for time = 3
Time = 3
Mesh stats
points: 6420
faces: 17388
internal faces: 15792
cells: 5502
boundary patches: 2
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 5026
prisms: 336
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 140
Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface topology
defaultFaces 812 870 ok (non-closed singly connected)
mysuplat 784 840 ok (non-closed singly connected)
Checking geometry...
Overall domain bounding box (-0.499997 -0.499997 0) (0.499997 0.499997 6)
Mesh (non-empty, non-wedge) directions (1 1 0)
Mesh (non-empty) directions (1 1 0)
***Number of edges not aligned with or perpendicular to non-empty directions: 5373
<<Writing 5855 points on non-aligned edges to set nonAlignedEdges
Boundary openness (2.10851e-17 9.47129e-17 5.47847e-18) OK.
Max cell openness = 2.60068e-16 OK.
Max aspect ratio = 28.184 OK.
Minumum face area = 0.000800175. Maximum face area = 0.0857147. Face area magnitudes OK.
Min volume = 0.000347927. Max volume = 0.00647358. Total volume = 4.70014. Cell volumes OK.
Mesh non-orthogonality Max: 75.6873 average: 14.825
*Number of severely non-orthogonal faces: 672.
Non-orthogonality check OK.
<<Writing 672 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 0.378278 OK.
Failed 1 mesh checks.
End
checkMesh warn us of the presence of 672 nonOrthoFaces besides the nonAlignedEdges.