Creating OpenFOAM meshes with Gmsh

Following on from Alex's post, I'd like to demonstrate how Gmsh can be used to create a mesh for OpenFOAM, using a 2D bifurcating network as a simple example. Meshes can be created interactively using a GUI or by writing a .geo file using Gmsh's own scripting language, which will often be more convenient. For easy modification of the geometry, it is useful to start with a definition of the relevant parameters:

//--------------------------------------------------------------------
// Geometry parameters - inputs
//--------------------------------------------------------------------
hp  =   1;              // parent channel height
hd1 = 0.5;              // 1st daughter channel height
hd2 = 0.5;              // 2nd daughter channel height
L   =  2;               // channel lengths
L2  =   1;              // Daughter channel contraction lengths
nCells =  15;           // number of cells in transverse direction
xCells = 30;            // number of cells in streamwise direction
//--------------------------------------------------------------------
// Geometry parameters - calculated
//--------------------------------------------------------------------
midp = hp*Sqrt(3)/2;
lp = hp/nCells;
ld1 = hd1/nCells;
ld2 = hd2/nCells;
dm = 2*midp/hp;
dx = (L2^2/(1+dm^2))^(1/2);
dy = dm*dx;
dx1 = 0.5*((hp-hd1)^2/(1+(1/dm)^2))^(1/2);
dy1 = dx1/dm;
dx2 = 0.5*((hp-hd2)^2/(1+(1/dm)^2))^(1/2);
dy2 = dx2/dm;

Next, we specify the grid points that define the geometry, based on the parameters above. The expression inside the parentheses is the point's ID number; the first three columns inside the braces are the x, y, z coordinates, and the 4th column denotes the prescribed mesh element size near that point.

//--------------------------------------------------------------------
// Points
//--------------------------------------------------------------------
// Junction - triangle
Point(1) = {0, hp/2, 0, lp};
Point(2) = {0, -hp/2, 0, lp};
Point(3) = {midp, 0, 0, lp};
// Junction - contractions
Point(4) = {dx+dx1, hp/2+dy-dy1, 0, ld1};
Point(5) = {midp+dx-dx1, dy+dy1, 0, ld1};
Point(6) = {dx+dx2, -(hp/2+dy-dy2), 0, ld2};
Point(7) = {midp+dx-dx2, -(dy+dy2), 0, ld2};

Between points, the mesh element size will interpolate, so grading could be introduced by placing points at the centre of the channels. The channels can be created by an alternative method, so we'll only define points for the junction/contraction part of the geometry. Points are joined together by lines, with the numbers in braces specifying the two points you wish to connect:

//--------------------------------------------------------------------
// Lines
//--------------------------------------------------------------------
// Junction - triangle
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 1};
// Junction - contractions
Line(4) = {1, 4};
Line(5) = {4, 5};
Line(6) = {5, 3};
Line(7) = {3, 7};
Line(8) = {7, 6};
Line(9) = {6, 2};

Arcs and splines can also be easily created by specifying a centre of rotation and a number of control points, respectively. Plane surfaces are created using a line loop. For the line loop, the numbers in braces specify, in order, the lines which constitute the perimeter of the surface. The line loop must be closed (i.e. it ends where it begins) and continuous (i.e. a negative line ID indicates the opposite direction).

//--------------------------------------------------------------------
// Surfaces
//--------------------------------------------------------------------
Line Loop(1) = {1, 2, 3};
Plane Surface(1) = {1};
Line Loop(2) = {4, 5, 6, 3};
Plane Surface(2) = {2};
Line Loop(3) = {7, 8, 9, 2};
Plane Surface(3) = {3};

Strictly speaking, lines 2 and 3 aren't needed, and surfaces 1, 2, and 3 could be combined into a single surface. However, by default, Gmsh uses an unstructured mesh and it is often a good idea to partition the geometry to constrain the mesh generation. As the channel sections are straight, it is convenient to use a structured mesh for them; for this, we use the extrude command to translate a line:

//--------------------------------------------------------------------
// Channels
//--------------------------------------------------------------------
pS[] = Extrude {-L, 0, 0} {
  Line{1};
  Layers{xCells};
  Recombine;
};
d1S[] = Extrude {dx*L/L2, dy*L/L2, 0} {
  Line{5};
  Layers{xCells};
  Recombine;
};

d2S[] = Extrude {dx*L/L2, -dy*L/L2, 0} {
  Line{8};
  Layers{xCells};
  Recombine;
};

The extrude command automatically creates all the necessary points, lines, and surface between the chosen line(s) and its translated counterpart. The first set of braces gives the x, y, z components of the extrusion; the second set of braces specifies a) the line(s) you wish to extrude, b) how many cells you wish to divide the extruded surface into, and c) an optional command to recombine the triangular cells into rectangles. The variables pS, d1S, and d2S contain the IDs of the newly extruded line (in e.g. pS[0]), the surface it creates (in e.g. pS[1]), and the joining lines. OpenFOAM requires meshes to be 3D, so we next have to extrude the entire surface by one layer of arbitrary thickness (in the z-direction). This is achieved using the same extrude command, but specifying surfaces as an argument rather than lines:

//--------------------------------------------------------------------
// Unit depth
//--------------------------------------------------------------------
zV[] = Extrude {0, 0, -1} {
  Surface{1,2,3,pS[1],d1S[1],d2S[1]};
  Layers{1};
  Recombine;
};

Similarly to the line extrusions, the variable zV contains all the extruded surfaces and volumes. Finally, we create "physical surfaces" that OpenFOAM will recognise as patches for boundary conditions:

//--------------------------------------------------------------------
// Physical surfaces
//--------------------------------------------------------------------
Physical Surface("inlet") = {zV[21]};
Physical Surface("outlet1") = {zV[27]};
Physical Surface("outlet2") = {zV[33]};
Physical Surface("walls") = {zV[7],zV[9],zV[13],zV[15],zV[20],zV[22],zV[26],zV[28],zV[32],zV[34]};
Physical Surface("topAndBottom") = {1,2,3,pS[1],d1S[1],d2S[1],zV[0],zV[5],zV[11],zV[17],zV[23],zV[29]};
Physical Volume("internalMesh") = {zV[1],zV[6],zV[12],zV[18],zV[24],zV[30]};

Note, all physical surfaces appear in the "boundary" file as type "patch", so you will need to change "topAndBottom" to type "empty". Naming this file "bifurcation2d.geo", we can generate the mesh using the following command:

gmsh bifurcation2d.geo -3

If the mesh is of poor quality, the flag -optimize is useful for optimising the mesh element quality. The flags -clmin float and -clmax float are also useful for constraining the minimum and maximum element sizes, respectively. This will create the file "bifurcating2d.msh". If this file is in the same directory as the OpenFOAM case, then to convert the resulting .msh file into an OpenFOAM mesh, run

gmshToFoam bifurcation2d.msh

If the .msh file is not in the same directory as the OpenFOAM case, the flag -case DIR needs to be used. Below, shows the resulting mesh in Paraview, and the solution for the steady-state velocity profile using the icoFoam solver.

We can easily change the geometry of the network by altering the inputs. For example, changing the size of the daughter channels to hd1=0.7 and hd2=0.3 creates the following mesh and velocity profile solution: