Generating a 2D mesh from CAD

Learning Objectives

This tutorial covers the basic functionalities of the curvilinear mesh generator that is provided with NekMesh, and is designed to show the main features of the two-dimensional mesh generator. Upon completion of this tutorial, the user will be able to

• create a simple boundary layer 2D mesh of a NACA aerofoil from a supplied CAD file in STEP format;
• add a refinement region to increase resolution locally;
• run a quick simulation with the IncNavierStokesSolver;
• visualise the mesh and the flow results.

The mesh configuration file

The design of the mesh generation system focuses on automation. As such, once provided with a CAD description of the domain only 4 numerical parameters are required to produce a mesh. The system uses a curvature based refinement and smoothness algorithms to obtain mesh sizings automatically. In this tutorial you will be introduced to using manual refinement fields to achieve more control over the mesh sizing.

NekMesh works on input and output files and both need to be specified as command line arguments. In the case of mesh generation the input is a .mcf, mesh configuration file, which is formatted similarly to the Nektar++ XML format. This file contains the instructions and parameters the mesh generation system needs, such as the input geometry name, the type of mesh to be generated, and the polynomial order at which to generate.

A barebones example of an .mcf file is:

<NEKTAR>
    <MESHING>
        <INFORMATION>
            <I PROPERTY="CADFile" VALUE="cyl.stp" />
            <I PROPERTY="MeshType" VALUE="2D" />
        </INFORMATION>
        <PARAMETERS>
            <P PARAM="MinDelta" VALUE="0.5" />
            <P PARAM="MaxDelta" VALUE="3.0" />
            <P PARAM="EPS" VALUE="0.1" />
            <P PARAM="Order" VALUE="4" />
        </PARAMETERS>
    </MESHING>
</NEKTAR>

Six key pieces of information are required to run the mesh generator:

• CADFile specifies the geometry file which will be meshed. It should typically be in STEP format, but some additional formats are supported: see the other tutorials for examples of this.
• MeshType tells the mesher what type of elements to produce. Depending on the dimension of the CAD, this can be one of 3D or 2D, or when a boundary layer is required, 3DBndLayer or 2DBndLayer.
• MinDelta, MaxDelta and EPS determine the minimum and maximum element sizes.
• Order defines the order of the high-order mesh.

The system automatically determines element sizing, $\delta$ as \(\delta = 2R\sqrt{\varepsilon(2-\varepsilon)},\) where $R$ is the radius of curvature of the CAD entity. Since curvature is typically unbounded, MinDelta and MaxDelta place limits on the value of $\delta$, so that $\delta_\min \leq \delta \leq \delta_\max$. A somewhat counterintuitive behaviour of this format is that, in the case of constant radius geometries, such as cylinders, changing the value of MaxDelta and MinDelta may not alter the mesh. In this setting, the EPS parameter can be used to control the sizing.

Creating a simple mesh of a cylinder

As a first step, let’s create a simple 2D cylinder in a rectangular domain. We do this from a .stp file (cyl.stp) which has been provided. We therefore run the command:

!NekMesh cyl.mcf cyl.xml

Notice there is no output. However if we supply the -v parameter, we can see information about the CAD file, and th mesh generation steps.

!NekMesh -v cyl.mcf cyl.xml

This produces a Nektar++ XML file, which you can then visualise. If you have a visualisation package such as ParaView or VisIt, you can convert the mesh to VTK format using FieldConvert.

!FieldConvert -f cyl.xml cyl.vtu

Alternatively, you might consider using the Tecplot format, either binary or plain ASCII.

!FieldConvert -f cyl.xml cyl.dat
!FieldConvert -f cyl.xml cyl.plt

Since these meshes are converted to a linear mesh, since this is what the visualisation software supports, the curved elements can look quite faceted. To counteract this, we can add an optional -n X command line option to FieldConvert, which increases the number of subdivisions used to view the elements through interpolation. For example, to subdivide each high-order element into 7 linear elements:

!FieldConvert -f -n 8 cyl.xml cyl-8.vtu

Visualising the mesh

In the above commands, the cyl.xml file is converted to a VTK file, cyl.vtu. For visualisation, there are a number of options:

• interactively, in the notebook: you can test this by executing the cell directly below this text;
• download the files generated by the commands in this notebook, and then use software such as ParaView or VisItto display them locally on your computer.

In this instance we suggest executing the cell below, which runs some basic Python code to read and display the mesh. You can use the mouse scroller to zoom in and out, and click/drag to rotate.

• If you click on the Z plane button , this will snap the view into two dimensions.
• The Surface with edges button will show cell boundaries and allow you to see the underlying mesh.

import pyvista as pv

# First read the VTK file
mesh = pv.read('cyl.vtu')

# Now create an itkwidget to visualise this in-browser.
pl = pv.PlotterITK()

# Add the mesh to the plotter.
pl.add_mesh(mesh, smooth_shading=True)

# Show the viewer.
pl.show(False)

Some important points to note in the visualisation of high-order meshes and fields:

• This visualisation shows a linear representation of the mesh. That is, inside each element, we are sampling the coordinates at discrete points and then performing linear interpolation between these. Therefore lower-order simulations may look jagged when the underlying error is actually rather small.
• Owing to this, in the visualisation above, you are viewing all integration points of the element. In triangular elements in particular, this can be somewhat visually odd: in Nektar++, we utilise a collapsed coordinate system to represent triangles, where two corners of a quadrilateral are merged into one. Thus the distribution of points within a triangular element is ‘squashed’ in one direction.

Adding a boundary layer

Since we want to perform a simulation of flow over the cylinder, we now need to do some additional refinement, in order to produce a suitable mesh for our solvers. The first step is to change the configuration to generate a boundary layer mesh. To do this, first change the mesh type in the configuration file cyl.mcf to read

<I PROPERTY="MeshType" VALUE="2DBndLayer" />

In addition to this, we need to provide information to specify the boundary to be refined and the level of refinement required. Add the following to the configuration file cyl.mcf:

<P PARAM="BndLayerSurfaces"    VALUE="5,6,7,8" />
<P PARAM="BndLayerThickness"   VALUE="0.25"    />
<P PARAM="BndLayerLayers"      VALUE="2"       />
<P PARAM="BndLayerProgression" VALUE="1.5"     />

These parameters control:

• the IDs of the CAD surfaces (or in this 2D case, the CAD curves) onto which the boundary layer will be added;
• the total thickness of the boundary layer;
• the number of layers in the boundary layer mesh; and
• the geometric progression used to define the reducing height of the boundary layer elements as they approach the surface.

A version that includes these modifications can be found in cyl-bl.mcf

!NekMesh -v cyl.mcf cyl.xml

Running a simulation

We can now try to run some simulations on this mesh. A session file, session_cyl.xml for the cylinder mesh is provided, which will execute a low Reynolds number incompressible simulation yielding classical vortex shedding. It is a relatively quick simulation and can run a large number of convective lengths in a short time. But, using the mesh configuration described above, the simulation will be quite under-resolved. You will notice that the vortices will shed at a significant angle to the freestream.

!IncNavierStokesSolver cyl.xml session_cyl.xml

Now lets visualize the results of the simulaiton. The simulation is run for 5000 timestep with $\Delta t=1\times10^{-2}$ and the results are dumpt each 50 steps. Hence, at the end of the simulation we will have 100 solution file. To visualize the results we first need to convert these files to an appropriate format for visulaization which can be done via the following python code:

for n in range(100):
    !FieldConvert -f -m vorticity cyl.xml session_cyl.xml cyl_$n\.chk out_cyl_$n\.vtu

The above command converts the .chk, i.e. check point files to .vtu files that can be visualize using Paraview, or inside the notebook. More information on the post-processing and visulaization can be found in the FieldConvert tutorials. Now, that all the solution files are converted, we can use the following script to visualize the animation of the flow over the time.

import pyvista as pv
import time
import ipywidgets as widgets

# Here we read the output data. Note that in the FieldConvert command above,
# we set the output file names to be cyl_0.vtu, cyl_1.vtu ...
# the " 'cyl_%d.vtu' % i " will produce the file names.
#
# In the subsequent sections, you should replace the "cyl_" part by the 
# name you give for the outputs.
meshes = [ pv.read('out_ud_%d.vtu' % i) for i in range(1,100) ]

plotter = pv.PlotterITK(notebook=True)
plotter.add_mesh(meshes[0], scalars='W_z')
viewer = plotter.show(False)

# we have 100 output files
nfiles = 100

i = 1
while True:
    viewer.geometries = meshes[i-1]
    i = (i+1) % nfiles
    time.sleep(0.25)

Adding a refinement region

The mesh which has been generated is very coarse and the simulation is lower order at $P = 2$. As such the simulation is highly under-resolved. A useful feature for this task is the ability to manually add refinement controls to the mesh generation, which can be used to refine the wake of the cylinder so that vortices can be captured accurately. This is done by adding a refinement tag to the cyl.mcf file, inside the MESHING XML tag, which for this example takes the form:

<REFINEMENT>
    <LINE>
        <X1> −1.0 </X1>
        <Y1> 0.0 </Y1>
        <Z1> 0.0 </Z1>
        <X2> 15.0 </X2>
        <Y2> 0.0 </Y2>
        <Z2> 0.0 </Z2>
        <R> 1.5 </R>
        <D> 0.4 </D>
    </LINE>
</REFINEMENT>

In this configuration, we define a line using two points, with coordinates (X1,Y1,Z1) and (X2,Y2,Z2). We also define the parameter R which is the radius of influence of the refinement around this line. D then defines the size of the elements within the refinement region. Add the aove settings to the cyl.mcf file to have a refined region.

If under your conditions the mesh has invalid elements you can also activate the variational mesh generation tool, which will attempt to perform interior deformation by adding the following tag inside the MESHING tag:

<BOOLPARAMETERS>
    <P VALUE="VariationalOptimiser" />
</BOOLPARAMETERS>

You should now be in a position to create a new mesh with your specifications and alter the simulation accordingly. You should note that increasing the resolution, you may affect the CFL condition and need to reduce the timestep. Change the mesh and solver parameters to produce different meshes and flow results, so that vortices can be captured accurately. Suggested parameters can be seen in the completed solution below.

Extra tasks

Task 1

Now that you have the refinement to the cyl.mcf file, generate a mesh using these new setting.

Hint the mesh generation command is !NekMesh -v cyl.mcf cyl.xml

Task 2

Run the simulation using the new mesh, convert and visulaize the results.

!FieldConvert -f -n 8 -m vorticity final/cyl-final.xml session_cyl.xml final/cyl-final.fld cyl-final.vtu

Take-a-ways:

It is worthy to keep in mind the following points:

• NekMesh can read CAD files in STEP format.

• In the Finite Difference Method, the boundary layer mesh usually has a progression rate smaller than $1.2$. However, Spectral/hp Element Method can be applied on meshes whose progression rates are larger than $1.5$.

End of the tutorial