Dynamic mesh in OpenFOAM

features in OpenFOAM is dynamic meshing. But there is almost no documentation. Here we build an example.

Case setup

Our case is a 100×35 mm rectangle with the 8x25mm obstacle on the top side. This obstacle will move left to right leaving 5mm space to vertical walls.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 0.001;

vertices
(
    (0.0 0 0) (0.0 0 1)
    (5.0 0 0) (5.0 0 1)
    (13.0 0 0) (13.0 0 1)
    (100.0 0 0) (100.0 0 1)
    (0.0 10.0 0) (0.0 10.0 1)
    (5.0 10.0 0) (5.0 10.0 1)
    (13.0 10.0 0) (13.0 10.0 1)
    (100.0 10.0 0) (100.0 10.0 1)
    (0.0 35.0 0) (0.0 35.0 1)
    (5.0 35.0 0) (5.0 35.0 1)
    (13.0 35.0 0) (13.0 35.0 1)
    (100.0 35.0 0) (100.0 35.0 1)
);
blocks
(
    hex (0 2 10 8 1 3 11 9) (5 10 1) simpleGrading (1 1 1)
    hex (2 4 12 10 3 5 13 11) (8 10 1) simpleGrading (1 1 1)
    hex (4 6 14 12 5 7 15 13) (87 10 1) simpleGrading (1 1 1)
    hex (8 10 18 16 9 11 19 17) (5 25 1) simpleGrading (1 1 1)
    hex (12 14 22 20 13 15 23 21) (87 25 1) simpleGrading (1 1 1)
);

edges
(
);

patches
(
    wall walls
    (
        (0 2 3 1)
        (2 4 5 3)
        (4 6 7 5)
        (6 14 15 7)
        (14 22 23 15)
        (16 8 9 17)
        (8 0 1 9)
    )
    patch obstacle
    (
        (20 12 13 21)
        (12 10 11 13)
        (10 18 19 11)
    )
    patch atmosphere
    (
        (22 20 21 23)
        (18 16 17 19)
    )
);

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

Initial water level is 20mm. This is achieved by system/setFieldsDict file

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      setFieldsDict;
}

defaultFieldValues
(
    volScalarFieldValue alpha.water 0
    volVectorFieldValue U ( 0 0 0 )
);

regions
(
    boxToCell
    {
        box ( 0 0 0 ) ( 0.100 0.020 0.001 );
        fieldValues
        (
            volScalarFieldValue alpha.water 1
        );
    }
 );

Note: here is no convertToMeters parameter, so all dimensions are in meters.

Other options are taken from “Breaking of a dam” tutorial.

Dynamic meshing setup

To use dynamic meshing features, two files should be created. The first is constant/dynamicMeshDict. It specifies dynamic mesh type and it’s parameters. In our case, it looks like below:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      dynamicMeshDict;

}
motionSolverLibs ( "libfvMotionSolvers.so"  );
dynamicFvMesh    movingConeTopoFvMesh;

movingConeTopoFvMeshCoeffs
{
  motionVelAmplitude (0.082 0 0); // 82mm/sec through X axis
  motionVelPeriod 1.5; // realy, its SEMIperiod

  leftEdge -0.020;
  leftObstacleEdge  0.030;
  rightObstacleEdge 0.050;

  left
  {
    minThickness 0.1e-3;
    maxThickness 0.2e-3;

  }

  right
  {
    minThickness 0.1e-3;
    maxThickness 0.2e-3;

  }
}

Here, we choose movingConeTopoFvMesh type with 82 mm/sec velocity amplitude and 1.5 sec semi-period (see details below).

The second file we should provide is constant/polyMesh/meshModifiers. It contains a list of mesh modifiers that will be applied to the mesh. Our file is:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      meshModifiers;

}
2
(
  right
  {
    type layerAdditionRemoval;
    faceZoneName rightExtrusionFaces;
    minLayerThickness 0.0008;
    maxLayerThickness 0.0012;
    oldLayerThickness -1;
    active on;
  }
   left
   {
    type layerAdditionRemoval;
    faceZoneName leftExtrusionFaces;
    minLayerThickness 0.0008;
    maxLayerThickness 0.0012;
    oldLayerThickness -1;
    active on;
  }
);

Here we define two modifiers with layerAdditionRemoval type. It’s parameters will be discussed below.

moivingConeTopoFvMesh

This modifier moves all points between two face zones named leftExtrusionFaces and rightExtrusionFaces. Points are moved according the law , where is velocity amplitude specified by motionVelAmplitude entry of movingConeTopoFvMeshCoeffs dictionary. Semi-period corresponds to motionVelPeriod value.

Other dictionary entries are used only when moivingConeTopoFvMesh is run in “automatic” mode. It happens when neither face zones, not point zones, cell zones or mesh modifiers is specified. In this case, moivingConeTopoFvMesh creates two layerAdditionRemoval modifiers, that are controlled via left and right sub-dictionaries. It also creates two face zones, one for each modifier. The first zone is named leftExtrusionFaces and includes all faces with centre’s X coordinate equal to . The second zone is named rightExtrusionFaces and includes all faces with centre’s X coordinate equal to . Note this values are hard-coded. It seems like originally leftObstacleEdge and rightObstacleEdge parameters were planned to set these values, but somebody forgot to use them.

layerAdditionRemoval

This modifier adds or removes cell layers near the face zone specified.

When is applied by the dynamic mesh engine (like moivingConeTopoFvMesh), it first checks the thickness of all cells touching the face in the face zone. If an average thickness increased relative to previous call and maximum thickness is greater than maxLayerThickness, than layerAdditionRemoval adds new layer of cells. In the same manner, if an average thickness decreased relative to previous call and maximum thickness is less than maxLayerThickness, than layerAdditionRemoval removes one layer of cells.

Entry oldLayerThickness specifies value for the first comparison with average layer thickness, when there is no value from the previous step.

Special note should be done about cell selection. Internal faces has two cells – owner and neighbour. It’s normal points from the cell with a smaller ID to the cells with the larger one. Boundary faces, obviously, have only owner cell, and their normals point out the cell. It seems like for both internal and boundary faces an algorithm chooses cell that is pointed to by the face normal. It means that boundary face normals should always be inverted.

Face zones

Face zone is a list of faces with normal flipping specified. List of all face zones is specified in constant/polyMesh/faceZones file.

FoamFile                            
{                                   
    version     2.0;                
    format      ascii;              
    class       regIOobject;        
    location    "constant/polyMesh";
    object      faceZones;          
}                                   

2
(
rightExtrusionFaces
{
    type faceZone;
    faceLabels List 5 (1 3 5 7 9);
    flipMap    List  5 (0 1 1 0 0);
}
leftExtrusionFaces
{
    type faceZone;
    faceLabels List 3 (2 4 6);
    flipMap    List  3 (0 1 1);
}
)

Each list item contains it’s type faceZone, list of face labels and list of boolean values that indicate if face normal should be flipped (relative to it’s original direction) or no.

The last thing is important, because layerAdditionRemoval uses normal direction for cell selection.

Face zone creation

To produce face zones, use topoSet utility. Here is the file system/topoSetDict:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      topoSetDict;

}

actions
(
{
        name    fLeft;
        type    faceSet;
        action  new;
        source  boxToFace;
        sourceInfo
        {
            box (0.0049 0 0) (0.0051 0.035 0.001);
        }
}
{
        name    fRight;
        type    faceSet;
        action  new;
        source  boxToFace;
        sourceInfo
        {
            box (0.0129 0 0) (0.0131 0.035 0.001);
        }
}
{
        name    cLeft;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (0.004 0 0) (0.005 0.035 0.001);
        }
}
{
        name    cRight;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (0.013 0 0) (0.014 0.035 0.001);
        }
}
{
        name    rightExtrusionFaces;
        type    faceZoneSet;
        action  new;
        source  setsToFaceZone;
        sourceInfo
        {
            faceSet fRight;
            cellSet cRight;
            flip    true;
        }
}
{
        name    leftExtrusionFaces;
        type    faceZoneSet;
        action  new;
        source  setsToFaceZone;
        sourceInfo
        {
            faceSet fLeft;
            cellSet cLeft;
            flip    true;
        }
}
);

In this file, we first create two face sets fLeft and fRight with all faces with X coordinate equal to 5mm and 13mm. Note that not only obstacle’s faces are used, but internal faces below them too.

Next, we define cell (yes, cell) sets cLeft and cRight. The first set includes cells left to the fLeft faces. In other words, it is the rightmost cell column before the obstacle. In the same manner, the second cell set is the leftmost column after the obstacle.

The last step is a face zones definition. It is done by setsToFaceZone source. As appears from it’s name, it creates face zone from the face set. But it has additional parameter cellSet. It seems like it refers to the cell set that face normals should point out from. More precisely, according to the setsToFaceZone.C:148, topoSet inverts normal if

• it is boundary face and it’s owner is not in cellSet
• it is internal face and it’s neighbour (but not owner) is in cellSet

The flip flag inverts the result of these rules.

These steps produce exactly what we need. Really, by default, all faces from fLeft have normals directed to the right (see fig.). All they also belongs to cells from cLeft. It means they aren’t touched by the rules above, but are then flipped by the flip flag.

As for faces from fRight set, they splits into two parts.

The first one is the obstacle’s right side (i.e. top part of fRight). It’s faces by default look to the left (as a boundary faces, that look outside domain). They also have owners from cRight cell set. As a result, they will be flipped only by the flip flag.

The second part of fRight includes faces below the obstacle right side. Initially, their normals look to the right. Note these faces have no owners in cRight cell set, but have neighbours. It implies that they will be flipped twice: by the second rule from the list above, and by the flip flag. So, they will remain the same as default, i.e. pointing to the right.

Resulting face zones are shown on the the figure right. Note that all normals look into cells that will be affected by the layer modifiers.

 

Case file

Here is a case archive: https://yadi.sk/d/LYyzg8HukjLqo