High-Order Implicit Large-Eddy Simulations of Flow over a NACA0021 Aerofoil

A. Flux Reconstruction

For a detailed overview of the multidimensional FR approach, the reader is referred to the work of Witherden et al. [6] and the references therein. In this section, we will outline salient numerical aspects of FR within the context of solving the three-dimensional Navier–Stokes equations on unstructured hexahedral grids. In conservative form, the Navier–Stokes equations can be expressed as

To solve the system using the FR approach, we start by decomposing the physical domain $\mathrm{\Omega }$ into a set of $N_E$ conformal, nonoverlapping, hexahedral elements. Take ${\mathrm{\Omega }}_i$ to refer to the $i$th element in this set. It is convenient, for reasons of both mathematical simplicity and computational efficiency, to work in a transformed space. We accomplish this by introducing a standard element $\tilde{\mathrm{\Omega }}$, which exists in a transformed space $\tilde{\mathit{x}}$. Take ${\mathcal{J}}_i(\tilde{\mathit{x}})$ to be the relevant Jacobian matrix for the $i$th element and $J_i(\tilde{\mathit{x}})$ its determinant. We proceed to associate this standard element with a set of $N_u$ solution points denoted by $\{{\tilde{\mathit{x}}}_i^{(u)}\}\in \tilde{\mathrm{\Omega }}$ and a set of $N_f$ flux points denoted by $\{{\tilde{\mathit{x}}}_i^{(f)}\}\in \partial \tilde{\mathrm{\Omega }}$. At each flux point, we define a polynomial correction function ${\mathit{g}}_i^{(f)}(\tilde{\mathit{x}})$ with the property that ${\mathit{n}}_i^{(f)}\cdot {\mathit{g}}_j^{(f)}({\tilde{\mathit{x}}}_i^{(f)})={\delta }_{ij}$, where ${\mathit{n}}_i^{(f)}$ is the unit normal vector associated with the $i$th flux point. The specific form of these correction functions determines various properties of the resulting FR scheme, including stability and accuracy [2,14,15]. In particular, by the appropriate choice of correction functions, one can recover the nodal DG scheme [2].

Inside of each element, we represent the solution as a degree $k$ polynomial $u_{i\alpha }\in P_k(\tilde{\mathrm{\Omega }})$. This polynomial can be defined using a nodal basis as

This is accomplished as follows. First, the solution polynomial is evaluated at each of the flux points. As the flux points on the faces of adjacent elements pair up, this results in two distinct values of the solution at a single point in physical space. By employing a local discontinuous Galerkin (LDG)-type method, it is possible to correct these discontinuous solutions to obtain an approximation, which is $C^0$ continuous at the flux points and the gradient of which sits in $P_k^3(\tilde{\mathrm{\Omega }})$. The gradients of this approximate solution, along with the original discontinuous solution defined by Eq. (7), can be used to compute a discontinuous flux within each as

Furthermore, using the solution gradient, it is also possible to compute an approximation of the inviscid and viscous fluxes at each flux point pair. These can then be combined to compute a normal transformed flux, which is common to both flux points. For the inviscid portion, this involves application of an approximate Riemann solver, whereas for the viscous portion, a LDG-type scheme can be employed. This common approximation, however, is generally different from that obtained by evaluating the transformed discontinuous flux of Eq. (8) at the flux point and dotting it with the associated normal vector.

The final step involves scaling the vector correction function associated with each flux point by the difference between the common and and discontinuous flux at each flux point. By summing the scaled corrections from all of the flux points and adding these onto ${\tilde{\mathit{f}}}_{i\alpha }^D$, the result will be $C^0$ continuous at each of the flux points. The divergence of this polynomial can then be taken to give the desired semidiscretized form of the governing system.

B. Aliasing-Driven Instabilities

Inside of a FR step, it is necessary, at several points, to perform polynomial projections: specifically, the projection of the transformed flux into the polynomial space of the solution within an element, projection of the transformed normal common interface flux into the polynomial space of the correction function divergence on the face of an element, and the projection of the divergence of the corrected flux into the space of the solution within an element. In the prescription outlined previously, these projections are performed by simply evaluating the relevant functional at the appropriate set of nodal basis points, a so-called collocation projection.

Although extremely efficient, collocation-type projections can result in aliasing, whereby energy from modes that are unresolved is erroneously transferred into those that are. This aliasing can impact not just the accuracy of the solution but also the stability. Aliasing can be mitigated through modal filtering [3], adding spectral vanishing viscosity [10], switching to a skew-symmetric formulation [11,12], and antialiasing via an approximate $L^2$ projection with a suitable quadrature rule [13]. Two of these approaches are considered in this study: modal filtering and antialiasing via an approximate $L^2$ projection. Modal filtering is perhaps the most straightforward approach to suppressing aliasing errors. It is relatively simple to implement and relatively cheap to apply at a given time step. However, it often requires careful parameter tuning in order to achieve a balance between stability and accuracy and even then can be overly dissipative. Antialiasing via an approximate $L^2$ projection, on the other hand, can be relatively complex to implement, requiring significant changes to a codebase, and can be relatively expensive to apply at a given time step. However, if the $L^2$ projection is sufficiently accurate, it can suppress aliasing without additional tuning parameters and without being overly dissipative.

C. Modal Filtering

It is known that the effects of aliasing are most pronounced in the highest-frequency modes of the expansion [16–18]. Therefore, one means of stabilizing a simulation is to filter the higher modes of the expansion [3] every $N$ time steps. The modal expansion of the solution can be written as

D. Antialiasing via Approximate $L^2$ Projection

The principle behind antialiasing is to compute the modal expansion coefficients of the desired polynomial exactly. Taking the projection of the nonlinear volume flux as an example, we note that in modal form the polynomial inside of the $i$th element can be expressed as

When antialiasing via an approximate $L^2$ projection, it is necessary to evaluate the relevant functional at the chosen set of quadrature points. This requires all arguments of the functional, for example, the solution and its gradients, to be interpolated to these points. Next, the functional itself must be evaluated at the abscissa and a summation performed to obtain the desired modal coefficients. As a consequence, antialiasing via an approximate $L^2$ projection, especially with a large number of points, is more computationally expensive than a simple collocation projection.

A. Overview

In this study, we use PyFR to simulate unsteady separated turbulent flow over a NACA0021 aerofoil in deep stall. The aerofoil is set at a 60 deg angle of attack to the oncoming freestream flow. The Reynolds number $R{e}_{\infty }$ (based on the aerofoil chord $c$, the freestream velocity $u_{\infty }$, the freestream density ${\rho }_{\infty }$, and the fluid viscosity $\mu$) is set at 270,000. The Mach number $M{a}_{\infty }$ (based on the freestream velocity $u_{\infty }$, the freestream density ${\rho }_{\infty }$, and the freestream pressure $p_{\infty }$) is set at 0.1. The setup exhibits significant unsteady dynamics and is hence a good test case for scale-resolving approaches such as detached-eddy simulation (DES), LES, and DNS. Moreover, significant experimental [19] and computational (from the DESider project) [20,21] data exist for the test case, including time-span-averaged pressure coefficient distributions on the aerofoil, associated lift and drag coefficients, and lift and drag coefficient spectra.

B. Domains and Meshes

Four domains were considered. All had an extent of $30c$ in the stream- and crosswise directions. However, the spanwise extent $L_z$ was varied. Specifically, values of $L_z=c$, $2c$, $4c$, $7c$ were considered. Previous experimental data were obtained with a span of $7.2c$ [19], and previous computational results suggest that simulations should employ a span of at least $4c$ [20]. The leading edge of the aerofoil was located at the center of the cross-streamwise plane. Riemann-invariant boundary conditions were applied in the far field, a no-slip adiabatic wall boundary condition was applied on the aerofoil surface, and a periodic condition was imposed in the spanwise direction.

For each domain, two quadratically curved hexahedral meshes were made using Gmsh [22], one for simulation with third-order solution polynomials (nominally fourth-order accurate in space), henceforth referred to as P3 meshes, and one for simulation with fourth-order solution polynomials (nominally fifth-order accurate in space), henceforth referred to as P4 meshes. All meshes were unstructured in the stream- and crosswise directions but uniformly extruded in the spanwise direction. Elements adjacent to the aerofoil were sized such that the first (Gauss–Legendre) solution point sat a distance equivalent to $y^{+}\approx 1$ from the wall (where $y^{+}$ is based on a flat plate boundary layer at Reynolds number 270,000). Elements of the P3 meshes had a characteristic size $0.05c$ in the stream- and crosswise directions near the aerofoil and a uniform extent $0.05c$ in the span-wise direction throughout the domain. Elements of the P4 meshes had a characteristic size $0.0625c$ in the stream- and crosswise directions near the aerofoil and a uniform extent $0.0625c$ in the spanwise direction throughout the domain. For both meshes, the majority of elements was concentrated near the aerofoil, as shown in Fig. 1.

In total, the P3 meshes had 80,840 elements per spanwise extent $c$, and the $P4$ meshes had 51,632 elements per spanwise extent $c$. However, since third-order solution polynomials require 64 solution points within each hexahedral element and fourth-order solution polynomials require 125, the total number of degrees of freedom per spanwise extent $c$ was similar across P3 and P4 meshes.

C. Methodology

The compressible Navier–Stokes equations with constant viscosity were solved using PyFR version 0.8.0 [6]. A DG scheme was used for the spatial discretisation, solution points were located as a tensor product of Gauss–Legendre points within each hexahedra, flux points were located as a tensor product of Gauss–Legendre points on each face of each hexahedra, common inviscid fluxes were computed via a RoeM-type Riemann solver [23], and common viscous fluxes were computed using a LDG approach [24]. An explicit RK45[$2\mathrm{R}+$] scheme [25] with adaptive PI time-step control [26] was used to advance the solution in time. A wall-resolved ILES approach was adopted, and hence no turbulence model was employed.

Simulations on P3 and P4 meshes with $L_z=c$, $2c$, $4c$, $7c$ were started from a freestream initial condition with $R{e}_{\infty }=27,000$. Each simulation was then run “natively” without any antialiasing for $100c/u_{\infty }$ time units using first-order solution polynomials to facilitate the passing of initial transients. After this warmup period, various simulations were restarted on each mesh with $R{e}_{\infty }=270,000$. Specifically, P3 meshes with $L_z=c$, $2c$, $4c$, $7c$ were restarted with third-order solution polynomials and 1) antialiasing via modal filtering, henceforth referred to as P3-F simulations; 2) antialiasing via an $L^2$ projection, henceforth referred to as P3-AA simulations; and 3) natively without any antialiasing, henceforth referred to as P3-N simulations. P4 meshes with $L_z=c$, $2c$, $4c$, $7c$ were restarted with fourth-order solution polynomials and 1) antialiasing via an $L^2$ projection, henceforth referred to as P4-AA simulations, and 2) natively without any antialiasing, henceforth referred to as P4-N simulations. All the aforementioned simulations, which are summarized in Table 1, were run for a further $400c/u_{\infty }$ time units, over which period data were extracted for analysis.

In the previous description, antialiasing via modal filtering refers to the application of Eq. (11) with $\kappa =36$, $s=16$, and ${\eta }_c=1$ every $N=50$ time steps. A visual representation of the resulting filter coefficients is shown in Fig. 2. Antialiasing via an $L^2$ projection on the P3 and P4 meshes, with third-order and fourth-order solution polynomials, respectively, refers to use of 9th-order and 11th-order Gauss–Legendre quadrature rules, respectively, to perform an approximate $L^2$ projection of the transformed flux into the polynomial space of the solution (within an element) and an approximate $L^2$ projection of the transformed normal common interface flux into the polynomial space of the correction function divergence (on the face of an element).

We note that the filter parameters were selected based on the results of a series of numerical experiments using a P3 mesh with $L_z=c$ and third-order solution polynomials. First, experiments were performed varying filter strength $s$ and cutoff ${\eta }_c$, with fixed $\kappa =36$ and $N=50$. Outcomes are presented in Table 2. Looking at the table, it is clear that only the strongest filter ($s=16$ and ${\eta }_c=1$) is stable up to $t=400c/u_{\infty }$. Next, for this strongest filter, with fixed $\kappa =36$, $s=16$, and ${\eta }_c=1$, experiments were performed varying application frequency $N$. Outcomes are presented in Table 3. Looking at the table, it is clear that only the highest application frequency ($N=50$) is stable up to $t=400c/u_{\infty }$.

A. Stability

It was found that P3-N and P4-N simulations quickly diverged for all $L_z$ considered. However, P3-F, P3-AA, and P4-AA simulations remained stable for all $L_z$. Consequently, all subsequent results are from P3-F, P3-AA, and P4-AA simulations. This finding clearly demonstrates the importance of suppressing aliasing-driven instabilities when using FR schemes.

B. Comparison with Experimental Data

1. Pressure Distributions

Figure 3 plots time- and span-averaged distributions of pressure coefficient $C_p$ on the pressure and suction surfaces of the aerofoil. All simulations achieved good agreement with the experimental data on the pressure surface. P3-AA and P4-AA simulations achieved increasingly good agreement with the experimental data on the suction surface as $L_z$ was increased, visually converging to the experimental data when $L_z\ge 4c$. P3-F simulations exhibited the same trend. However, they failed to visually converge to the experimental data, even when $L_z=7c$.

2. Time-Averaged Force Coefficients

Figures 4, 5, and 6 compare the time-averaged lift coefficient $C_l$ and time-averaged drag coefficient $C_d$ with experimental data [19] and computational results from the DESider project [20,21], for P3-F, P3-AA, and P4-AA simulations, respectively Table 4. Figures 7 and 8 plot percentage errors in $C_l$ and $C_d$, relative to experimental data [19], for PyFR simulations and computational results from the DESider project [20,21], respectively.

As one would expect, the overall trends are similar to those observed for the $C_p$ distributions of Fig. 3. When $L_z=c$, all of the PyFR simulations overpredict both $C_l$ and $C_d$. This is especially evident for the P3-F simulation with $L_z=c$. However, as $L_z$ is increased, the PyFR simulations converge toward the experimental data, and when $L_z\ge 4c$, the P3-AA and P4-AA simulations exhibit less than 2.5% deviation in $C_l$ and $C_d$. Finally, we note that the P4-AA simulation with $L_z=7c$ achieves the best agreement with experimental data and that the resulting errors in $C_l$ and $C_d$ are lower than any of those obtained from the DESider project [20,21].

3. Force Spectra

Figure 9 compares the power spectral density (PSD) of $C_l$ with related experimental data of [19], where $St$ is a Strouhal number based on the freestream velocity $u_{\infty }$ and the chord $c$.

To obtain the PSD, $C_l$ was sampled every $0.025c/u_{\infty }$ time units, and the PSD was computed using Welch’s averaged periodogram method with windows of length 4096 samples and a shift between windows of 10 samples. It should be noted that, while the PyFR data were obtained from $C_l$ (the time- and surface-averaged lift coefficient), the experimental was obtained from a sectional lift coefficient at a fixed spanwise location. Consequently, it is only meaningful to compare frequencies for peak 1 and peak 2 (see Fig. 9), which are presented in Table 5. We see that the P3-AA and P4-AA simulations with $L_z=4c$ and $L_z=7c$ achieved very good agreement with the experimental results.

C. Discussion

1. Effect of $L_z$

Increasing $L_z$ toward the experimental value of $7.2c$ was found to give consistently better agreement with experimental data, in terms of both time-averaged statistics and spectral characteristics. This appears to be an improvement compared to previous results from the DESider project [20,21], which, as illustrated in Fig. 8, exhibited a degree of scatter in time-averaged statistics as $L_z$ was varied, including cases in which time-averaged statistics became worse with increasing $L_z$.

To better understand $L_z$ dependence, Fig. 10 plots time-averaged nondimensional pressure contours, with streamlines, on the midspan plane obtained from the P3-AA simulations for the four different $L_z$. In all cases, there are two primary vortices behind the suction surface of the aerofoil, which gives rise to a low-pressure region. As $L_z$ is increased, these vortices move farther away from the suction surface, and the extent of the pressure drop in the wake is also reduced. These observations suggest that as $L_z$ is increased the two primary vortices become weaker.

Further insight is gained by plotting isosurfaces of Q criterion, which define vortical structures [27], obtained from P3-AA simulations with $L_z=c$ and $L_z=4c$. Figure 11 plots such isosurfaces colored by nondimensional density and cut by the midspan plane. Small-scale vortical structures for the $L_z=c$ case are elongated and predominately perpendicular to the spanwise direction, whereas those for the $L_z=4c$ case appear more isotropic. Also, the large region of low density (indicating the presence of a large-scale vortical structure) behind the suction surface is more pronounced for the $L_z=c$ case compared with the $L_z=4c$ case. This is further illustrated by Fig. 12, which plots colour maps of nondimensional density on the midspan plane, and Fig. 13, which plots isosurfaces of Q criterion colored by nondimensional density and cut by a plane described in Fig. 12.

2. Effect of Antialiasing via Modal Filtering and Approximate $L^2$ Projection

Antialiasing via an $L^2$ projection is found to give consistently more accurate results than antialiasing via modal filtering. Figure 14 plots time-averaged nondimensional pressure contours, with streamlines, on the midspan plane obtained from P3-F and P3-AA simulations with $L_z=c$. For the P3-F case, the two primary vortices are closer to the surface of the aerofoil than those obtained by the P3-AA simulation. The P3-F simulation also has a larger low-pressure region, which indicates that the primary vortices of the P3-F simulation are stronger than those of the P3-AA simulation.

Further insight is gained by plotting isosurfaces of the Q criterion, obtained from P3-F and P3-AA simulations with $L_z=c$. Figure 15 plots isosurfaces of the Q criterion colored by nondimensional density and cut by the midspan plane. It is clear that the P3-AA simulation resolves a wider range of small-scale vortical structures compared to the P3-F simulation. The inability of the filtered P3-F simulations to accurately resolve small-scale structures likely underlies their observed lack of convergence to the experimental data as $L_z$ is increased.

D. Computational Cost

A total of 12 simulations are run on a variety of different GPU-accelerated machines: local clusters at Imperial College, London, Emerald at the Science and Technology Facilities Council Rutherford Appleton Laboratory, the Wilkes cluster at Cambridge University, and Piz Daint at the Swiss National Supercomputing Center. All these machines consist of NVIDIA GPUs and are targeted using the CUDA backend of PyFR.

To compare computational cost across all the simulations, we choose to conduct performance testing on single machine, Piz Daint, which is comprised of NVIDIA K20X GPUs. This is accomplished by running $0.1c/u_{\infty }$ extra time units after the final solution for each simulation. The number of NVIDIA K20X GPUs is selected to be $10L_z/c$ (effectively weak scaling). Hence, when $L_z=c$, a total of ten GPUs are employed, whereas at $L_z=7c$, a total of 70 GPUs are used. Depending on the order of accuracy and the type of antialiasing approach, GPU memory load varies from 5 to 10% of capacity. This setup mirrors how the actual simulations are performed in practice with a greater number of GPUs being used at larger spans and loading toward the limit of strong scaling.

Measurements of the wall clock time, excluding that required for file input/output, can be seen in Fig. 16. From the figure, it is clear that the P4-AA simulations are quickest for all $L_z$. This is attributed to the fact that the P4 meshes are significantly coarser than the P3 meshes. The increase in element size serves to offset the more severe time-step restriction that comes from increasing the solution polynomial order, so as to permit a larger overall time step to be taken as shown in Fig. 17. This in turn offsets the higher wall clock time per time step for the P4-AA simulations, shown in Fig. 18, leading to an overall reduction in time to solution.

Finally, we note that, while model filtering is generally considered a cheaper route to antialiasing than an approximate $L^2$ projection, in this case it leads to a higher wall clock time. This was because the adaptive PI time-step controller [26] selects a smaller time step for P3-F simulations than P3-AA simulations, as can be seen in Fig. 17. This reduction in time-step size offsets the reduced wall clock time per time step for the P3-F simulations, shown in Fig. 18, leading to an overall increase in time to solution. Further studies revealed that this effect is proportional to the filter application frequency $N$, as shown in Fig. 19. Thus, when employing an adaptive PI time-step controller, it is important to apply the filter as infrequently as possible. The interaction between modal filtering and adaptive time stepping remains a topic for further research.