Setting up industrial simulations using Incompressible Navier-Stokes Solver

This blog explains how to set up industrial simulations using the Incompressible Navier-Stokes (IncNS) Solver in Nektar++. The geometry here is the Imperial Front Wing (IFW), an unraced front-wing configuration of the McLaren MP4-17D race car. Combined with a rolling wheel, it is a representative model for modern open-wheel motorsports, especially Formula 1 front-wing designs. We focus here on the simulation setup for this geometry, but the considerations apply to any large, 3D simulation for the IncNS solver.

Simulation setup

The major components for setting up:

Preprocessing: Geometry and Mesh, I/O considerations
Solver setup: Expansion, stabilisation, parameters, boundary conditions
Postprocessing: Dealing with large solution dumps and efficient use of FieldConvert filters.

The basic command to run the IncNS solver in parallel is:

$(which mpirun) IncNavierStokesSolver mesh.xml session.xml -v

Preprocessing

The starting point of the simulation is a valid CAD of the geometry (rig_CAD.stp). The meshing procedure comprises two steps: a) generating a linear mesh and b) converting to a high-order curvilinear mesh.

The first step, the linear mesh, is generated using commercial Finite Volume meshing software using one of the accepted input modules. NekMesh is used for the higher-order mesh generation. We use a conformal, unstructured, mixed element mesh with tetrahedron and prismatic elements and a single “macro” prism layer (rig.ccm) to obtain a linear mesh (rig_linear_mesh.xml).

NekMesh -v rig.ccm:ccm rig_linear_mesh.xml

In the second step, the linear mesh (rig_linear_mesh.xml) is projected on the CAD boundary (rig_CAD.stp) and higher-order nodes are added on edges, faces, and elemental interiors to obtain the projected mesh (rig_projected_mesh.xml). The projection order (order) should be at least the polynomial order of the simulation you plan to run.

NekMesh -m projectcad:file=rig_CAD.stp:order=4 rig_linear_mesh.xml rig_projected_mesh.xml -v

If there are any invalid elements after the projection step, they can be removed by the following step:

NekMesh -m linearise:invalid rig_projected_mesh.xml rig_valid_projected_mesh.xml -v

The macro prism layer is split using an iso-parametric approach to refine the boundary layer to the desired first cell height and distribution. The following command splits the macro prism layer into five layers with a growth rate of 1.2 on the selected surfaces.

NekMesh -m bl:surf=4,5,6:layers=5:r=1.20 rig_valid_projected_mesh.xml mesh.xml -v

The final, valid higher-order mesh obtained is mesh.xml. By default, it has two separate COMPOSITES for the two different types of elements in our mesh (tetrahedrons and prisms). We suggest using the HDF5 format for efficient I/O for large geometries.

NekMesh mesh.xml mesh.nekg

Solver setup

The next step is to set a valid session file for the large 3D simulations (session.xml). The basic structure of a valid Nektar++ session file can be found here. The section does not provide a full outline of the session file. Only things to consider for 3D simulations will be indicated. A sample session file for this simulation can be found here.

We set up a simulation here using Taylor-Hood elements, using a polynomial order of P=4,3 for velocity and pressure. The EXPANSIONS section is defined here as combining one-dimensional basis functions and increased quadrature points. Quadrature points should be defined considering PDE and geometrical aliasing; more details can be found in the paper here. Velocity and pressure are defined separately for different composites containing tetrahedrons and prisms. The snippet below shows how it is done for the velocity at P=4.

<EXPANSIONS>

<E COMPOSITE="C[0]" BASISTYPE="Modified_A,Modified_A,Modified_B" NUMMODES="5,5,5" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0" NUMPOINTS="7,7,6" FIELDS="u,v,w" />

<E COMPOSITE="C[2]" BASISTYPE="Modified_A,Modified_B,Modified_C" NUMMODES="5,5,5" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0,GaussRadauMAlpha2Beta0" NUMPOINTS="7,6,6" FIELDS="u,v,w" />

</EXPANSIONS>

The CONDITIONS section contains various tags detailed in the user guide. To stabilise the simulation, the following options are recommended:

<SOLVERINFO>

</SOLVERINFO>

<PARAMETERS>

</PARAMETERS>

A conjugate gradient is often the best choice for solving large matrix systems. We use a conjugate gradient solver with static condensation applied to the global matrix (IterativeStaticCond). The absolute tolerance approach has been tested to work well for these simulations (AbsoluteTolerance). Different preconditioners and settings are needed to set up the linear system for velocity and pressure due to the varying stiffness of the matrices arising from their formulation. The pressure system should be solved with at least a tolerance of 1e-4 to achieve an accurate solution. A reliable set of defaults is outlined below:

<GLOBALSYSSOLNINFO>

<V VAR="u,v,w">

</V>

<V VAR="p">

</V>

</GLOBALSYSSOLNINFO>

The definition of variables and boundary conditions have been skipped as they tend to be case-dependent. There are a few pointers to keep in mind while running large 3D simulations:

If possible, restart from an existing lower-fidelity solution. Restarting from a uniform field often imposes strong stability restrictions on the timestep.
Track the CFL value actively to adjust the timestep accordingly. The CFL value is just an indicator of the simulation’s stability —more often than not, the simulation will crash for a CFL > 1.
Start increasing the polynomial order from P=2,1 and keep increasing it as soon as it runs stably. The sooner you reach the desired polynomial order, the better.
Use HDF5 to manage the size of the output and provide I/O gains. To use this, Nektar++ needs to be compiled using HDF5 and MPI. The simulation can be run using HDF5 I/O with the following command.

$(which mpirun) IncNavierStokesSolver mesh.xml session.xml --io-format Hdf5 --use-hdf5-node-comm -v

Postprocessing

Managing the output files from large simulations is difficult. Factors like storage, memory requirements, file management, and interpreting the output further complicate this. 3D solution files (solution.fld) can be heavy, and we suggest the following ways to deal with this:

Use forces and point/plane probes (using the HistoryPoints filter) to track the evolving flow simulation.
Extract solutions on surfaces of interest using appropriate FieldConvert modules to deal with smaller surface files.
Use FieldConvert output as Filters while the simulation runs to have the desired output as FieldConvert checkpoints.

An example of how to extract surfaces from a 3D solution (instantaneous or average) and use it to calculate surface variables and wall shear stress is given below:

Extract the surface mesh of interest from the original 3D mesh using NekMesh

NekMesh -m extract:surf=1 mesh.xml surface1.xml

Extract the surface variables and wall shear stress on the same surface using the appropriate FieldConvert modules. Please note that the bnd can be different from the surf, investigate your mesh setup for this.

$(which mpirun) FieldConvert -f -m extract:bnd=0 mesh.xml session.xml solution.fld solution_surface1_surfvar.fld

$(which mpirun) FieldConvert -f -m wss:bnd=0:addnormals=1 mesh.xml session.xml solution.fld solution_surface1_wss.fld

Convert the extracted solution on the surface, using the surface mesh, to the VTU format for visualisation.

FieldConvert -f surface1.xml solution_surface1_surfvar.fld solution_surface1_surfvar.vtu:vtu:highorder