[jns-v]

Dear Folks,

I have absolutely no clue as to what's wrong here, thus I'm not sure if my issue
is actually related to the meshing process or if it is a more generic problem. I'll
just try to give you the information I think may be neccessary in order to find the
root cause of my issues, so please bear with me.

I'm working a project as part of a CFD course at university. I'm trying to model
natural convection in a cylinder which is heated at the bottom (const. heat flux
of 40 kW/m²) but otherwise has adiabatic walls.Initially the cylinder is filled with 20 L of water at 300 K. I'm also using Boussinesq approximation for the buoyancy effects.

The OpenFOAM release I'm working with is 8.20201114-1 (openfoam.org)

In order to solve this in a reasonable time I'm using a wedge of 10° (I know it's probably
a bit large, but that's not the culprit). I followed https://cfd.direct/openfoam/user-guide/v8-blockMesh/ to setup the vertices etc.

So everything basically went well, mesh generation, solving with (buoyantPimpleFOAM),
only the results weren't exactly physical and also didn't match too well the reference calculations conducted with Ansys Fluent. Anyway that's how I know the blockMeshDict follows the syntax and topology etc.

One of the options I came up with was using buoyantBoussinesqPimpleFoam (the incompressible version of buoyantPimpleFoam as far as I understand things, which I deemed to be better suited for this case anyway) which is only available in the openfoam.com version, so I went ahead and tried to install that. Installing that from source actually worked kind of even though the linker reported some errors.

Anyway the mesh generation using blockMesh (com version) failed which is why I was unable to start solving anything.

As both versions org and com are not available for my distro (Arch Linux) my tutor suggested to go with a docker image. Same behaviour.

I tried to go back to openfoam-org v8 which at least generated a proper mesh before but now that's no longer working as well and which is where I finally get to the point of this post. I hope you're still with me.

So my blockMeshDict is:

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

convertToMeters 1;      // all values given in m // I understand thats "scale" in .comVersion

// geometry parameters
V 20.0;                 // Volume in l
V #calc "$V * 1e-3";    // convert to SI units
dhRatio 2.335;          // diameter-to-height ratio

// spatial discreization
dr 2.0e-3;                 // spacial step size in radial direction in mm
dz $dr;                 // spacial step size in axial direction in mm
pi #calc "acos(-1.0)";  // PI
angle 5.0;              // angle in degrees (10 degree wedge means +-5° from neutral plane
angleRad #calc "$angle * $pi / 180.0"; // angle in radians

// calculate diameter and height from volume an D/H-ratio
D #calc "pow(4 * $V / $pi * $dhRatio, 1.0/3.0)"; // diameter in m
R #calc "$D / 2.0";      // radius in m
H #calc "$D / $dhRatio"; // height in m

// Cell counts
Nr #calc "round($R / $dr)"; // Number of cells in radial direction
Nz #calc "round($H / $dz)"; // Number of cells in z-direction (i.e. top to bot)

// domain boundaries
xMin 0;
xMax #calc "$R * cos($angleRad)";
yMax #calc "$R * sin($angleRad)";
yMin #calc "$R * sin($angleRad) * -1.0";


vertices
(
    ($xMin $xMin $xMin) // 0,A
    ($xMax $yMin $xMin) // 1,B
    ($xMax $yMax $xMin) // 2,C
    ($xMin $xMin $H)    // 3,D
    ($xMax $yMin $H)    // 4,E
    ($xMax $yMax $H)    // 5,F
);

blocks
(
    hex (0 1 2 0 3 4 5 3) ($Nr 1 $Nz) simpleGrading (1 1 1)
);

edges
(
    arc 1 2 (92.67 0 0)
    arc 4 5 (92.67 0 370.67)
);

boundary
(
    wall
    { // cylinder outer wall
        type    wall;
        faces
        ( (1 2 4 5) );
    } // end wall
    top
    { // top
        type    wall;
        faces
        ( (3 4 5 3) );
    }; // end top
    bottom
    { // top
        type    wall;
        faces
        ( (0 1 2 0) );
    }; // end top
    front
    { // face pointing clockwise around z-axis
        type    wedge;
        faces
        ( (0 1 4 3) );
    } // end front
    back
    { // face pointing counterclockwise around z-axis
        type    wedge;
        faces
        ( (0 2 5 3) );
    } // end front
    axis
    { // face pointing counterclockwise
        type    empty;
        faces
        ( (0 3 3 0) );
    } // end front
);

Which initially produced a proper Mesh (took it from a former working test case). Note it creates prisms at the symmetry axis which makes sense...:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 8-15ae41985af0
Exec   : checkMesh
Date   : Jun 30 2021
Time   : 12:02:41
Host   : "Arch"
PID    : 65248
I/O    : uncollated
Case   : /home/jonas/Studium/tuw/2021S/CFD/Projekt/simulation/validationCase300K_leastSquares
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           16745
    internal points:  0
    faces:            33026
    internal faces:   16282
    cells:            8232
    faces per cell:   5.9898
    boundary patches: 6
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     8148
    prisms:        84
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing 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                  
    wall                84       170      ok (non-closed singly connected)  
    top                 98       197      ok (non-closed singly connected)  
    bottom              98       197      ok (non-closed singly connected)  
    front               8232     8415     ok (non-closed singly connected)  
    back                8232     8415     ok (non-closed singly connected)  
    axis                0        0        ok (empty)                        

Checking geometry...
    Overall domain bounding box (0 -0.0170089 0) (0.194412 0.0170089 0.167156)
    Mesh has 2 geometric (non-empty/wedge) directions (1 0 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Wedge front with angle 5.00002 degrees
    Wedge back with angle 5.00002 degrees
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (1.02257e-14 -5.18071e-15 -7.09525e-18) OK.
    Max cell openness = 3.06561e-16 OK.
    Max aspect ratio = 2.00621 OK.
    Minimum face area = 3.44308e-07. Maximum face area = 6.76938e-05.  Face area magnitudes OK.
    Min volume = 6.85157e-10. Max volume = 1.33606e-07.  Total volume = 0.00055274.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.323205 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

Now using the .com version (this is actually on a Windows host using WSL2 but I got the same for generic install on Arch and using docker on Arch). It doesn't create any prisms and I guess that's why I get those 2 failed mesh checks:

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _7bdb509494-20201222 OPENFOAM=2012
Arch   : "LSB;label=32;scalar=64"
Exec   : checkMesh
Date   : Jun 30 2021
Time   : 12:05:26
Host   : DESKTOP-JD8R1M6
PID    : 2799
I/O    : uncollated
Case   : /home/jonas/sim/testMesh
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Time = 0

Mesh stats
    points:           16830
    internal points:  0
    faces:            33110
    internal faces:   16282
    cells:            8232
    faces per cell:   6
    boundary patches: 6
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     8232
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
 ***Total number of faces on empty patches is not divisible by the number of cells in the mesh. Hence this mesh is not 1D or 2D.
    Cell to face addressing 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
    wall                84       170      ok (non-closed singly connected)
    top                 98       198      ok (non-closed singly connected)
    bottom              98       198      ok (non-closed singly connected)
    front               8232     8415     ok (non-closed singly connected)
    back                8232     8415     ok (non-closed singly connected)
    axis                84       170      ok (non-closed singly connected)

Checking faceZone topology for multiply connected surfaces...
    No faceZones found.

Checking basic cellZone addressing...
    No cellZones found.

Checking geometry...
    Overall domain bounding box (0 -0.0170088 0) (0.194412 0.0170088 0.167156)
    Mesh has 2 geometric (non-empty/wedge) directions (1 0 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Wedge front with angle 4.99999 degrees
    Wedge back with angle 4.99999 degrees
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-1.42018e-14 -2.82449e-16 1.44999e-17) OK.
    Max cell openness = 2.82981e-16 OK.
    Max aspect ratio = 2.00621 OK.
 ***Zero or negative face area detected.  Minimum area: 0
  <<Writing 84 zero area faces to set zeroAreaFaces
    Min volume = 6.85153e-10. Max volume = 1.33605e-07.  Total volume = 0.000552737.  Cell volumes OK.
    Mesh non-orthogonality Max: 0 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
 ***Max skewness = 1.32253e+147, 84 highly skew faces detected which may impair the quality of the results
  <<Writing 84 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

Failed 2 mesh checks.

End

My plan thus was to generate the mesh using the working .org-Version 8 and solve using the .com version with the incompressible solver. However after reinstalling v8 I can't even generate the mesh anymore (same blockMeshDict, same everything) but I get a floating point exception. (output shortend due to character count restriction see full log in attachment).

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 8-15ae41985af0
Exec   : blockMesh
Date   : Jun 30 2021
Time   : 11:32:48
Host   : "Arch"
PID    : 50013
I/O    : uncollated
Case   : /home/jonas/Studium/tuw/2021S/CFD/Projekt/simulation/testMesh
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Creating block mesh from
    "system/blockMeshDict"
[...]
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/usr/lib/libc.so.6"
#3  Foam::face::centre(Foam::Field<Foam::Vector<double> > const&) const at ??:?
#4  Foam::blockDescriptor::check(Foam::Istream const&) at ??:?
#5  Foam::blockDescriptor::blockDescriptor(Foam::dictionary const&, int, Foam::Field<Foam::Vector<double> > const&, Foam::PtrList<Foam::blockEdge> const&, Foam::PtrList<Foam::blockFace> const&, Foam::Istream&) at ??:?
#6  Foam::block::block(Foam::dictionary const&, int, Foam::Field<Foam::Vector<double> > const&, Foam::PtrList<Foam::blockEdge> const&, Foam::PtrList<Foam::blockFace> const&, Foam::Istream&) at ??:?
#7  Foam::block::New(Foam::dictionary const&, int, Foam::Field<Foam::Vector<double> > const&, Foam::PtrList<Foam::blockEdge> const&, Foam::PtrList<Foam::blockFace> const&, Foam::Istream&) at ??:?
#8  void Foam::PtrList<Foam::block>::read<Foam::block::iNew>(Foam::Istream&, Foam::block::iNew const&) at ??:?
#9  Foam::blockMesh::createTopology(Foam::IOdictionary const&, Foam::word const&) at ??:?
#10  Foam::blockMesh::blockMesh(Foam::IOdictionary const&, Foam::word const&) at ??:?
#11  ? in "/opt/OpenFOAM/OpenFOAM-8/platforms/linux64GccDPInt32Opt/bin/blockMesh"
#12  __libc_start_main in "/usr/lib/libc.so.6"
#13  ? in "/opt/OpenFOAM/OpenFOAM-8/platforms/linux64GccDPInt32Opt/bin/blockMesh"
[1]    50013 floating point exception (core dumped)  blockMesh

Obviously i can omit the checkMesh in this case I suppose.
I don't understand why this fails after reinstallation when I didn't change the blockMeshDict at all (in fact I setup a base case and only always copy that to create a new one).

Also note that I get the same behaviour runnun blockMesh on two tutorial cases "freePiston" and "DLR_A_LTS" which also use wedge geometries so I think my blockMeshDict is fine but I have a more serious issue here.

Please if you need more info I'm happy to provide it. I hope you made it this far anyway so any help is much appreciated as this is a bit beyond me...

Best regards
Jonas

[jns-v]

Okay so after I've grown a some fresh grey hair I figured out that the solver buoyantBoussinesqPimpleFoam doesn't let me specify externalWallheatFlux BC so it's unusable for me anyway.
The remaining issue I'd like to be able to resolve is the floating point exception when running blockMeshDict.

[HPE]

couldn't reproduce the issue with v2106.

might be useful to report the issue into .org bug tracker.

[jns-v]

Originally Posted by HPE

couldn't reproduce the issue with v2106.

might be useful to report the issue into .org bug tracker.

May I ask, did you check the mesh? i.e. did it produce prisms for you at the axis? as you are on v2106.


Anyway thanks for looking at it. I'm reluctant to file a bug report because as mentioned the blockMeshDict used to compile just fine but after reinstalling .org version I got the fpe so I thought it hast to have sth to do with my system setup rather than OF itself. But I'll reconsider this and may setup a proper reproducable case in order to file a bug report.


Thanks.

[HPE]

no, I haven't checked the mesh. But checking in terms of what?

[Ship Designer]

Hello Jonas,

I tried your blockMeshDict in OF v8 and it works without errors, I've noticed a few things though.

scale 1;

works in both OF .org and .com.

The edges of type arc are superfluous if in your y-direction the cell density is only 1, because the edge of a single cell is straight anyway.






I tried entering a density of 10 and got an odd looking block:






I then changed the arc definition as follows:

edges
(
    //arc 1 2 (92.67 0 0)
    //arc 4 5 (92.67 0 370.67)
    arc 1 2 ($R 0 0)
    arc 4 5 ($R 0 $H)
);

and this looks better I think:




Then I've noticed that the vertex ordering in two faces, bottom and back, had wrong ordering, although blockMesh didn't complain about it. To make definition of faces easier and less prone to errors, there's a simpler way I discovered in one of the blockMesh tutorial cases. Faces for patches can be defined either with a list of vertices, ordered counterclockwise for the face normal pointing outward (right thumb rule), or, by the block name and face number. For the latter blockMesh takes care of correct vertex ordering. So in your case, give the block a name:

blocks
(
    name Wedge hex (0 1 2 0 3 4 5 3) ($Nr 1 $Nz) simpleGrading (1 1 1)
);

and then replace the vertex lists in the faces definition with the name and face number pair:

boundary
(
    wall
    { // cylinder outer wall
        type    wall;
        faces
        ( (Wedge 1) );
    } // end wall
    top
    { // top
        type    wall;
        faces
        ( (Wedge 5) );
    }; // end top
    bottom
    { // top
        type    wall;
        faces
        ( (Wedge 4) );
    }; // end top
    front
    { // face pointing clockwise around z-axis
        type    wedge;
        faces
        ( (Wedge 2) );
    } // end front
    back
    { // face pointing counterclockwise around z-axis
        type    wedge;
        faces
        ( (Wedge 3) );
    } // end front
    axis
    { // face pointing counterclockwise
        type    empty;
        faces
        ( (Wedge 0) );
    } // end front
);

I run blockMesh with your case with both OF v8 and v1912 and it works okay, see attached log files. Interesting to note is the discrepancy between the wedge angles. I'm afraid I have no idea why your blockMesh crashes.

[jns-v]

Originally Posted by HPE

no, I haven't checked the mesh. But checking in terms of what?

The .com version of blockMesh does not produce prisms at the axis but highly skewed hex blocks whereas the .org version correctly implements prisms at the axis.


So my question was if you also recognized this.

[jns-v]

Originally Posted by Ship Designer

Hello Jonas,

I tried your blockMeshDict in OF v8 and it works without errors, I've noticed a few things though.

scale 1;

works in both OF .org and .com.

The edges of type arc are superfluous if in your y-direction the cell density is only 1, because the edge of a single cell is straight anyway.






I tried entering a density of 10 and got an odd looking block:






I then changed the arc definition as follows:

edges
(
    //arc 1 2 (92.67 0 0)
    //arc 4 5 (92.67 0 370.67)
    arc 1 2 ($R 0 0)
    arc 4 5 ($R 0 $H)
);

and this looks better I think:




Then I've noticed that the vertex ordering in two faces, bottom and back, had wrong ordering, although blockMesh didn't complain about it. To make definition of faces easier and less prone to errors, there's a simpler way I discovered in one of the blockMesh tutorial cases. Faces for patches can be defined either with a list of vertices, ordered counterclockwise for the face normal pointing outward (right thumb rule), or, by the block name and face number. For the latter blockMesh takes care of correct vertex ordering. So in your case, give the block a name:

blocks
(
    name Wedge hex (0 1 2 0 3 4 5 3) ($Nr 1 $Nz) simpleGrading (1 1 1)
);

and then replace the vertex lists in the faces definition with the name and face number pair:

boundary
(
    wall
    { // cylinder outer wall
        type    wall;
        faces
        ( (Wedge 1) );
    } // end wall
    top
    { // top
        type    wall;
        faces
        ( (Wedge 5) );
    }; // end top
    bottom
    { // top
        type    wall;
        faces
        ( (Wedge 4) );
    }; // end top
    front
    { // face pointing clockwise around z-axis
        type    wedge;
        faces
        ( (Wedge 2) );
    } // end front
    back
    { // face pointing counterclockwise around z-axis
        type    wedge;
        faces
        ( (Wedge 3) );
    } // end front
    axis
    { // face pointing counterclockwise
        type    empty;
        faces
        ( (Wedge 0) );
    } // end front
);

I run blockMesh with your case with both OF v8 and v1912 and it works okay, see attached log files. Interesting to note is the discrepancy between the wedge angles. I'm afraid I have no idea why your blockMesh crashes.

Hey Ship Designer,


First of all thanks for your dedicated answer.
The naming of the block really makes the syntax much more readable.
I wasn't really aware of how to use this properly.


As for my OpenFOAM v8 installation, I've got a feeling that I messed up with some systemwide installed libraries during the process of compiling installing uninstalling etc so I'm kinda stuck with the docker image for now


Best regards,
Jonas

[olesen]

Originally Posted by jns-v

The .com version of blockMesh does not produce prisms at the axis but highly skewed hex blocks whereas the .org version correctly implements prisms at the axis.


So my question was if you also recognized this.

You may want the older merge points instead of topological merge if you indeed have wedge geometries. Can specify in the blockMeshDict or as a blockMesh option (see blockMesh -help-full).

[olesen]

Originally Posted by Ship Designer

Hello Jonas,

I tried your blockMeshDict in OF v8 and it works without errors, I've noticed a few things though.

scale 1;

works in both OF .org and .com.

Indeed, the scale keyword has been in OpenFOAM since 2008: https://github.com/olesenm/OpenFOAM-...d6ec06b1ead41c

This precedes the refactoring of blockMesh as a library: https://github.com/olesenm/OpenFOAM-...47100e83f77a52

which was largely done in ramp-up for adding a paraview module to display blockMesh layouts: https://github.com/olesenm/OpenFOAM-...b9c370adf38aee