Flow around a NACA 0012 airfoil
In this tutorial, the simulation around a NACA 0012 airfoil at $Re=5000$ and $Ma=0.4$ is considered. First, we explain how to set the main flow parameters. We then describe the evaluation of lift and drag and visualization of the flow field. Finally, we show how to use the sponge zone to remove artificial reflections from the outflow boundary, so that a clean acoustic field is retained.
Contours of velocity magnitude for the flow around a NACA 0012 airfoil.
Copy the naca0012 tutorial folder to your desired working directory.
cp -r $FLEXI_TUTORIALS/naca0012 .
Compiler options
Make sure that FLEXI is compiled with the cmake options listed in the following table.
| Option | Value | Comment |
|---|---|---|
| CMAKE_BUILD_TYPE | Release | |
| FLEXI_EQYNSYSNAME | navierstokes | |
| FLEXI_PARABOLIC | ON | |
| FLEXI_MPI | ON | optional |
Table: Cmake options for the NACA 0012 simulation.
The standard settings are sufficient for this example. To check whether they are set, change to your build folder and open the cmake GUI
ccmake [flexi root directory]
If necessary, set the above options and then compile the code by issuing
make
Mesh Generation with HOPR
The mesh file used by FLEXI is created from the external linear mesh NACA0012_652.cgns and a 3rd order boundary description NACA0012_652_splitNg2.cgns using HOPR
hopr parameter_hopr.ini
This creates the mesh file NACA0012_652_Ng2_mesh.h5 in HDF5 format. If HOPR is not available, the mesh file is supplied in this tutorial.
Flow Simulation with FLEXI
The simulation setup is defined in parameter_flexi.ini. The initial condition is selected via the variable vector RefState=(/1.,0.990268069,0.139173101,0.,4.4642857/) which represents the vector of primitive solution variables $(\rho, u, v, w, p)^T$. The chosen velocity vector $(u,v)^T$ yields an angle of attack of $\alpha=8^\circ$ and a velocity magnitude of 1. The first supplied RefState is given the number 1, the second is given the number 2 and so on.
IniRefState = 1 : the initial condition uses RefState 1 for the initial flow field solution.
IniExactFunc = 1 : exact function routine for initialization, case 1 initializes a freestream state based on IniRefState.
Material properties are given in the following. Based on the ideal gas law, we get
$$Ma=1/\sqrt{\kappa p/\rho}=0.4$$
Note that in this non-dimensional setup the mesh is scaled such that the chord length is unity, i.e. $C=1$. Then to arrive at $Re=\rho u C / \mu = 5000$, the viscosity is set to
$$ \mu = \rho u C / Re = 1/Re = 0.0002 $$
| Property | Variable | Value |
|---|---|---|
| dynamic viscosity $\mu$ | mu0 | 0.0002 |
| ideal gas constant $R$ | R | 4.4642857 |
| Prandtl number | Pr | 0.72 |
| isentropic coefficient $\kappa$ | kappa | 1.4 |
Table: Material properties set in the parameter file.
Numerical settings
The DG solution on the mesh is represented by piecewise polynomials and the polynomial degree in this tutorial is chosen as $N=3$.
The main code settings are displayed in the following table:
| Variable | Description | Value |
|---|---|---|
| N | Polynomial degree | 3 |
| MeshFile | Mesh file to be used | NACA0012_652_Ng2_mesh.h5 |
| tend | end time of the simulation | 10 |
| Analyze_dt | time interval for analysis | 0.01 |
| nWriteData | dump solution every n’th Analyze_dt | 10 |
| CFLscale | 0.9 | |
| DFLscale | 0.9 |
Table: Numerical settings.
Boundary conditions
The boundary conditions were already set in the mesh file by hopr. Thus, the simulation runs without specifying the boundary conditions in the FLEXI parameter file. The freestream boundaries of the mesh are Dirichlet boundaries using the same state as the initialization, e.g. IniRefState by default.
The boundary types and states used by the simulations can always be checked in the initial console output of FLEXI. In order to run a quick test run, first set
maxiter=1
in your parameter file to enforce the simulation to stop after one time step. Then run the code
flexi parameter_flexi.ini
The code will finish after one time step. The initialization output contains the boundary conditions
Name Type State Alpha
BC_wall 3 0 0
BC_inflow 2 0 0
BC_outflow 2 0 0
BC_zminus 1 0 1
BC_zplus 1 0 -1
Type 2 identifies the weak Dirichlet boundary condition, type 3 the adiabatic wall. The boundary conditions in $z$ direction are not relevant for this 2D example and are realized as periodic boundaries.
Suppose we want to change the wall boundary condition from adiabatic to isothermal.
- This is done by first defining a second RefState. Add to the parameter file the second of the following lines
RefState=(/1.,0.990268069,0.139173101,0.,4.4642857/)
RefState=(/1.,0.990268069,0.139173101,0.,4.9107143/)
Note that in the second RefState the pressure was increased by 10\%. From the ideal gas law, we then know that the temperature also is raised by 10\%, since the density has been kept constant. Also note that the isothermal wall ignores the velocity components but only uses the thermodynamic quantities.
- Second, overwrite the initial wall boundary condition
BoundaryName=BC_wall
BoundaryType=(/4,2/)
Here, 4 implies the use of the isothermal boundary condition and 2 indicates the RefState to be used for the wall. Note that the temperature is set indirectly by the RefState given in terms of $\rho$ and $p$.
To check the newly set boundary condition, run the code once more for one iteration. The output should now look like this
Boundary in HDF file found | BC_wall
was | 3 0
is set to | 4 2
..................................................
Name Type State Alpha
BC_wall 4 2 0
BC_inflow 2 0 0
BC_outflow 2 0 0
BC_zminus 1 0 1
BC_zplus 1 0 -1
Now the wall boundary condition is isothermal using the temperature defined by the second RefState listed in your parameter file. A complete overview over the boundary conditions and how to use them is given in Section 4.3.2 of the user guide.
Before starting your simulation, remember to disable the maxiter line again by setting it to -1.
Running the code
We proceed by running the code in parallel. For example using 4 processors, use the following command
mpirun -np 4 flexi parameter_flexi.ini
On a 2012 laptop with core i5 processor, this simulation takes about 40 minutes.
Evaluation of the lift and drag forces
The forces acting on the airfoil are one of the main desired output quantities from the simulation. They are calculated on the fly during runtime. The associated flags in the parameter file are
CalcBodyForces=T
WriteBodyForces=T
The first line activates the calculation of the forces at each Analyze_dt, the second line enforces output of the forces to a .dat file.
The body forces are a good measure for convergence. In the context of time-dependent flows this determines whether the solution has reached a quasi steady state. The following figure shows the $x$ and $y$ components of the force acting on the airfoil until TEnd=10.
Resulting forces on the airfoil up to $t=10$.
The lift and drag coefficients can be easily calculated by rotating these forces from the computational reference frame to the one of the freestream.
From the forces, it is clear that the steady state has not yet been reached and the simulation must be run further. Before we proceed with the simulation, we will nonetheless examine the preliminary results to check the quality of the simulation.
Wall velocities
Due to the weak coupling between the grid cells and to boundaries, boundary conditions are enforced weakly, e.g. by applying a specific flux. This adds largely to the stability of the scheme. However, as a result the no-slip condition at the wall is not exactly fulfilled by the numerical solution. Rather, it is approximated as far as the resolution allows.
Evaluation of the velocity vector near the wall helps quantifying this error, which can be seen as a quality measure for the near wall resolution. In the parameter file, the computation and output to a .dat file of the average and extreme values of the wall velocity at every Analyze_dt are activated by setting
CalcWallVelocity=T
WriteWallVelocity=T
During the computation, we get output like the following:
Wall Velocities (mean/min/max) :
BC_wall 2.737875251E-02 3.639096145E-04 6.228999317E-01
In our case, the wall velocity is on average at about 3\% of the freestream velocity, reaching a peak of 60\%. This peak typically occurs at the quasi-singularity at the trailing edge. To decrease this deviation from the theoretical no-slip condition, either the wall-normal mesh size must be decreased or the polynomial degree increased.
It is important to note that both of the above measures will, besides increasing the number of degrees of freedom, decrease the time step, which directly affects the computational time. Thus, it is important to achieve an acceptable trade-off between the acceptable error and the computational time.
In this tutorial, the observed slip velocity is deemed uncritical and we proceed with the same resolution.
Visualization
To visualize the solution, the State-files must be converted into a format suitable for Paraview. Issue the command
mpirun -np 4 posti_visu parameter_postiVisu.ini parameter_flexi_navierstokes.ini NACA0012_Re5000_AoA8_State_0000000.0*
to generate the corresponding vtu-files using supersampling with $(NVisu+1)^3$ points in each cell. The vtu file can then be loaded into Paraview. The following figure shows a visualization of the density at $t=10$. The levels of $0.99<\rho<1.01$ are chosen to reveal the acoustic radiation from the airfoil.
Contours of density with $0.99<\rho<1.01$ at $t=10$ for the flow around a NACA 0012 airfoil.
The large scale vortex shedding of the wake due to the high angle of attack is clearly visible. Acoustic radiation from the airfoil can also be observed.
Now, a problem becomes apparent: the vortex street propagating towards the outflow boundary results in a second, artificial acoustic source at the outflow boundary. This is one of the fundamental problems in direct aeroacoustic computations. Before we proceed with the simulation, we will now make use of the sponge zone functionality of FLEXI to remove this artificial source.
Remove outflow reflections using the sponge zone
The sponge zone introduces a dissipative source term to the discrete operator, which is only active in a user-specified region, typically upstream of the outflow boundary. We use the sponge zone to dampen the vortices convected downstream before they hit the outflow boundary. See Flad et al. 2014 for the background of our sponge zone implementation.
In order to activate the sponge zone, set
SpongeLayer=T
in the parameter file.
Most of the other parameters are already set for the current simulation:
| Variable | Description | Value |
|---|---|---|
| SpongeShape | Set shape of sponge: 1: cartesian | 1 |
| damping | Damping factor of sponge | 1. |
| xStart | Coordinates of start position of sponge | (/2.0,0,0/) |
| ramp (for SpongeShape=1) | ||
| SpongeDistance | Set shape of sponge: 1: Cartesian | 3.0 |
| sponge region. | ||
| SpongeDir | Direction vector of the sponge ramp | (/1,0,0/) |
| (for SpongeShape=1) | ||
| SpongeBaseFlow | Type of baseflow to be used for sponge | 4 |
| 4: moving average (Pruett baseflow) | ||
| tempFilterWidth | Temporal filter width used to advance | 2.0 |
| Pruett baseflow in time. | ||
| SpongeViz | Write a visualization file of the | T |
Table: Sponge settings
The source term is of the form
\begin{equation}
\tilde{U}_t=U_t – d\sigma(\vec{x}) \left( U-U_B \right)
\end{equation}
First, damping determines the strength of the source term, i.e. $d$ in the above equation. It is dependent on the mean convection velocity, the desired amount of amplitude reduction and the thickness of the sponge zone. Typically, some trial and error is necessary to obtain an appropriate value. In non-dimensional calculations, i.e. velocity and length scale are of $\mathcal{O}(1)$, $d=0.1 … 2$.
Ramping of the source term from 0 is necessary to avoid reflections at the sponge interface. If such reflections occur, it is necessary to choose a thicker sponge ramp so that the source term is ramped up less steeply. We choose a parallel ramp by setting SpongeShape=1. The ramp’s starting position, thickness and direction are set by xStart, SpongeDistance and SpongeDir, respectively. These parameters govern the shape function $\sigma(\vec{x})$ which smoothly ramps the source term from 0 to 1.
It is important to emphasize that the regions where the source term is active should not be interpreted physically. They should be considered a part of the boundary condition. Likewise, they must be placed sufficiently far downstream of the airfoil such that they do not influence the near field solution. With the chosen settings, the sponge zone starts one chord behind the airfoil and is ramped up to 1 at the outflow boundary, located 4 chords behind the airfoil.
In order to visualize the ramping function $d\sigma(\vec{x})$, set SpongeViz=T. For that purpose, change the project name to e.g. “test” and set maxiter=1 to immediately stop the simulation after one iteration. This will output a .vtu file ready for visualization. After that, don’t forget to change back the project name and reset maxiter=-1.
Now, we have to choose the desired baseflow ($U_B$). For the current configuration, the moving average (SpongeBaseFlow=4) is appropriate. It produces a mean field slowly progressing in time, which adapts to the potential flow around the airfoil. The parameter tempFilterWidth determines the effective time window for the moving average. It should be chosen somewhat larger than the largest time scales to be damped. Its value of 2.0 here is chosen based on the frequency of the oscillations in the body forces seen in Figure \ref{fig:naca0012_bodyforces}.
The moving average base flow needs an initial value. One option is to provide an initial flow field from a file using the SpongeBaseflowFile parameter. If this parameter is not set, the code will initialize the base flow with the same values as the solution itself. Thus, in the case of a fresh computation, the base flow will be initialized with the IniExactFunc used to initialize the solution. In the case of a restart, like in the present tutorial, the base flow will be initialized with the state file provided for the restart of the simulation.
Note that using the moving average baseflow, the code will dump additional *baseflow*.h5 files to the hard drive, which are necessary to restart the simulation with the corresponding base flow from a given time. If these files are available and exhibit the same project name as the current simulation, the code automatically detects the right base flow file to initialize the base flow at the restart time. However, if in that situation SpongeBaseflowFile is still provided in the parameter file, the moving average will be restarted using the flow field provided in the parameter file, which may not be desirable.
Restarting the simulation
If the simulation run is interrupted or if you decide to run the simulation further than TEnd=10, FLEXI can easily be restarted. In the current setting, the code dumps the solution in $0.1$ time units intervals. Since we want to progress our simulation further, set TEnd=25 in the parameter file. Since we have now turned on the sponge zone, it is also advisable to modify the project name, i.e.
ProjectName=NACA0012_Re5000_AoA8_SP
If the state files are in the current folder, restart the simulation as follows
mpirun -np 4 flexi parameter_flexi.ini NACA0012_Re5000_AoA8_State_0000010.000000000.h5
Note that it is also possible to change the polynomial degree of the simulation during restart. This can be used to generate a rough initial solution at low $N$ and then switch to a high quality approximation at higher $N$ later. The code will automatically project the solution onto the new polynomial basis at startup. Restart is however only possible using the same mesh file.
Once the simulation is done, visualize the results again following the instructions in the visualization section above.
The following two figures show the resulting instantaneous density contours with and without sponge zone. Clearly, the source at the outflow is damped to uncritical levels. Fine-tuning of the sponge parameters may improve the results.
Contours of density with $0.99<\rho<1.01$ at $t=25$ for the flow around a NACA 0012 airfoil with moving average sponge zone.
Contours of density with $0.99<\rho<1.01$ at $t=20$ for the flow around a NACA 0012 airfoil without sponge zone.
Another look at the body forces indicates that a quasi-steady state has been reached.
Resulting forces on the airfoil up to $t=25$.
Two-dimensional computation
The laminar flow around this airfoil is in nature two-dimensional. Nonetheless, we ran this simulation using a three-dimensional code by imposing periodic boundary conditions and only one mesh element in the spanwise direction. This means we
- Compute one unecessary variable (momentum in spanwise direction),
- Compute three-dimensional fluxes for all variables,
- Compute one unecessary gradient,
- Use several degrees of freedom in the spanwise direction due to the high order ansatz in each element.
To circumvent this, FLEXI provides the option to perform truely two-dimensional calculations. To enable this, we need to set the flag FLEXI_2D to ON during configuration. FLEXI will then perform two-dimensional computations, but you will need to provide a
mesh that consists of only one element in the third dimension. Only then will you be able to start the computation. Since this is already the case for the mesh used in this tutorial,
we can start a two-dimensional computation right away after we compile our code with FLEXI_2D set to on. So switch to your build directory, set the CMAKE option FLEXI_2D to on, configure and compile your code once more. Then switch back to the tutorial folder.
You can now immediately run the code again using the same parameter file, but your new executable. All the options are the same for two-dimensional and three-dimensional computations. For compatibility reasons, you will have to specify vectorial parameters using
three dimensions (e.g. the momentum in the RefState), but the third dimension will be ignored. Pay attention to how much faster your code will run when using the two-dimensional version.
By default the State-files will be written like three-dimensional data by extruding the solution in the third dimension so the files will be compatible with all already existing post-processing tools. You can set the option Output2D to TRUE if you want to
write two-dimensional files, thus reducing the size of the output considerably. Also keep in mind that no matter if two- or three-dimensional, the number of dimensions of arrays will always be the same (one dimension simply beeing of size 1 for two-dimensional
computations) and also the number of variables will be the same. This means we will still have a momentum in the third dimension, which is simply set to zero (no calculations will be done for this momentum).
Using record points
Let’s assume we are interested in the temporal evolution of the flow variables in the vicinity of the profile’s upper and lower surfaces. Using the State-files, we will always have the whole flow field. If we want to have a very fine temporal resolution or our
field is very large, the cost of writing and storing these files can get prohibitive. As a remedy, FLEXI provides an option to sample specific points in the flow field at a very fine temporal resolution (down to a single time step) and only output those
specific points. We call these points record points.
If record points should be used, three distinctive steps have to be taken:
- Create the record points themselves
- Run a simulation with active record points
- Visualize the results
The first and third step are performed with pre- and postprocessing tools respectively, while a simulation using record points simply requires to set the specific options in FLEXI.
Prepare the record points
The coordinates of the record points are defined using a POSTI tool called preparerecordpoints, which is build when the respective CMake option is set. The tool takes a single parameter file as an input, so it can be run like
posti_preparerecordpoints parameter_recordpoints.ini
A sample parameter file is included in the NACA0012 tutorial folder. You need to specify the mesh that should be used to create the record points and a project name. The options NSuper and maxTolerance control the algorithm that tries to find
the record points in the mesh – a good default for NSuper is to set it to at least twice of the polynomial degree of your mesh. The main part of the parameter file is the actual definition of the record points. They are created in so called sets,
and each set is part of a named group. In the example, we have two sets – one for the suction side and one for the pressure side of the airfoil – each assigned to a separate group. There are several different types of sets available, which can be used
to create single points, lines or planes of record points. For the NACA example, we are using the boundary layer plane type of set. This special type of plane is defined by a spline which is projected onto the nearest boundary. Along this projected line,
a user-specified number of points is distributed in an equidistant way. The plane is then created by following the normal vector of the boundary at each of the points for a specified distance. Along the normal direction, the record points can be stretched to cluster
them near a wall.
The definition of one of the groups looks like this:
GroupName = suctionSide
BLPlane_GroupID = 1
BLPlane_nRP = (/20,30/)
BLPlane_nCP = 2
BLPlane_CP = (/0.9,0.014,0.5/)
BLPlane_height = 0.05
BLPlane_CP = (/0.999,0.001,0.5/)
BLPlane_height = 0.05
BLPlane_fac = 1.04
The first parameter creates a group which is used for identification purposes – each set needs to be assigned to a group. Then, the new set is defined by setting the respective GroupID option, in this case BLPlane_GroupID.
The ID is set to 1 to assign it to the first group. Following this general definitions, the actual parameters of the set follow. For the boundary layer plane, we first define the number of the points used in wall-tangential and
wall-normal direction using the BLPlane_nRP option. The next six options define the spline that is projected onto the boundary. First we set the number of control points in the boundary, then define the coordinates of each
of the points and the height in the wall-normal direction at this point. The last option defines the stretching that is used in the wall-normal direction.
To find out how the other types of sets are defined, have a look at the descriptions of their parameters.
If you now run the preparerecordpoints tool using the syntax given above, the algorithm will first calculate the physical coordinates of the recordpoints following your definitions. Then, it will attempt to identify the element in the mesh that contains each
record point and the coordinates in the reference element therein. This is done by inverting the mapping from reference to physical space using Newton’s method. This information is needed to later interpolate the polynomial solution to exactly the
positions of the record points. The results are then written to a file called ‘ProjectName_RPSet.h5’. If you set the doVisuRP option to true, you will also get a visualization of the record points that can be viewed in ParaView.
Using record points
To actually use the recordpoints during a simulation, you simply need to set a few options in the FLEXI parameter file. The import ones are the following:
RP_inUse = T
RP_DefFile = NACA0012_RPSet.h5
RP_SamplingOffset = 1
With the first option, the general usage of the record points system can be turned on or off. The RP_DefFile option must then be set to the name of the file that has been created using the preparerecordpoints tool. You can specify that the solution should
be sampled at every RP_SamplingOffset timestep.
If you now run the simulation, besides the State-files additional RP-files will be created, containing the value of the conservative variables at each record point at every sampling time step.
Visualizing the solution at the record points
The visualizerecordpoints tool is used to post-process the recordpoints data. It can be used to merge several of the RPSet.h5 files together and create a longer timeseries. Besides simply visualizing the timeseries, there are a multitude of
post-processing options available. You can use the time series to determine the mean and fluctuating part of the solution or perform spectral analysis using FFTs. For boundary layer planes, different types of turbulent quantities like the
skin friction can be directly computed.
As an example, we want to calculate the mean flow $\bar{U}$ and the temporal fluctutations $U’$ at our two record points planes. Again, a sample parameter file for the tool is given in the tutorial folder. The options set are
ProjectName = NACA0012
GroupName = suctionSide
GroupName = pressureSide
RP_DefFile = NACA0012_RPSet.h5
OutputTimeAverage = T
doFluctuations = T
where besides setting a name for the project and giving the path to the file that contains the record point definition, we define the name of the groups we want to visualize. They must be the ones we defined when using the preparerecordpoints tool.
Regarding actual evaluation, we want to calculate the temporal average and the fluctuations, simply calculated as $U’=U-\bar{U}$. We did not specify any variables that we want to visualize, in that case all conservative variables that are stored in the
RPSet.h5 files will be used. If you want to visualize a specific or a derived variable (e.g. the pressure), it can simply be set in the parameter file like this:
VarName = Pressure
You can run the tool using
posti_visualizerecordpoints parameter_visualizeRecordpoints.ini NACA0012_Re5000_AoA8_RP_*
This will take all the time samples recorded during our simulation as input. For each of the planes, a separate .vts file contaning the temporal average will be written. The fluctuations are time-resolved data, and due to limitations in the VTK file format
every time step has to be written to a single file. To avoid creating a lot of files in the working directory, those will be stored in a subfolder called timeseries. In the working directory, .pvd files for the fluctuations are created.
These files can be opened with ParaView and will contain the complete time series with the correct time value.