Implicit large eddy simulations for NACA0012 airfoils using compressible and incompressible discontinuous Galerkin solvers

We present implicit Large Eddy Simulations for NACA0012 airfoils at various Reynolds numbers (Re = 1× 10, Re = 1× 10 and Re = 1× 10) and Angles of Attack (0 ◦ ≤ AoA ≤ 10 ◦) using two discontinuous Galerkin formulations. On the one hand, we use a compressible solver based on a nodal DGSEM formulation and supplemented with a stabilising split-form formulation (Pirozzoli) and Roe interface fluxes. On the other hand, we use an incompressible DG-Fourier formulation that uses the interior penalty parameter to provide localised dissipation. Both solvers enable high order computations by varying the polynomial order inside mesh elements, which are here set to P=3 and P=4. We provide results of aerodynamic coefficients and pressure distributions using both solvers to show how they are able to provide underresolved flows that agree well with experimental data and well established solvers (Xfoil or Ansys-Fluent).


Introduction
High order Discontinuous Galerkin (DG) methods provide accurate solutions by enabling arbitrarily high polynomial approximations inside each grid element. For high order polynomials, the numerical errors are not distributed along all wave-numbers but localised at high wave-numbers [23,32,18,33,29]. This characteristic of high order methods results in very accurate simulations with low dissipative and dispersive errors. Although this characteristic seems a-priori beneficial for well resolved simulations, when computing under-resolved Large Eddy Simulations (LES), it can prove difficult to obtain stable simulations. In implicit (or under-resolved) Large Eddy Simulations (iLES), the smallest numerical eddies are larger than would have been in a finer mesh, leading to numerical under-resolution (i.e. coarse grid or low polynomial order) and aliasing [7]. Various methods have been proposed to stabilise under-resolved computations with aliasing. Among others, split forms or skew symmetric variants [35,6]), localised interior penalty fluxes [12], over-integration [26,25,31] or filtering [5] may be incorporated into the solver to stabilize the computations and remove or alleviate the aliasing.
Contrarily to low order methods, high order methods do not have enough inherent numerical dissipation in under-resolved simulations, to dissipate large flow structures (when compared to Kolmogorov scales). Therefore, computation of iLES flows using high order DG solvers require localised dissipative mechanisms to dissipate flow structures close to cut-off size. In what follows, we compare two dissipative stabilising mechanisms that enable the simulation of turbulent under-resolved flows. On the one hand, we use a compressible formulation with an energy conserving split-form and dissipation through Roe fluxes [24]. On the other hand, the incompressible solver uses the viscous discretisation through interior penalty formulation to enhance stability [12]. We challenge both formulations with a NACA0012 airfoil at various angles of attack in turbulent regimes, to explore both accuracy and stability. We compare simulated results to experimental data and simulations using low order methods (Xfoil and Ansys-Fluent).

Methodologies
We first introduce the two different mechanisms used to stabilise both compressible and incompressible high order DG formulations. The explanation included here is brief and aims only at introducing the fundamental concepts and motivating ideas. Further details can be found in the following references by the authors [24,12].
The 3D Navier-Stokes equations can be written as: where u is the vector of conservative variables u = (ρ, ρv 1 , ρv 2 , ρv 3 , ρe) T in compressible solvers. For incompressible solvers u = (v 1 , v 2 , v 3 ) T and Eq.(1) is complemented with ∇ · u. Details on the definition of inviscid and viscous solvers can be found in [24,12]. To derive discontinuous Galerkin schemes, we consider Eq. (1) for one mesh element el, multiply by a locally smooth test function φ j , for 0 ≤ j ≤ P , where P is the polynomial degree, and integrate on el: We can now integrate by parts the inviscid fluxes, F e , integral to obtain a local weak form of the equations (one per mesh element): where n is the normal vector at element boundaries ∂el. We replace discontinuous fluxes at inter-element faces by a numerical inviscid flux, F * e , to obtain a weak form for the equations for each element, where, we have omitted the fluxes at external boundaries, for simplicity. This set of equations for each element is coupled through the inviscid fluxes F * e and governs flow behaviour. Note that one can proceed similarly and integrate by parts the viscous terms (see [3,12]), but here for simplicity we retain the volume integral.
The non-linear inviscid and viscous terms that can be discretised to control dissipation in the numerical scheme have been underlined. Riemman solvers are the classic option to include numerical dissipation in DG schemes [34,4], since they naturally arise when discretising the nonlinear terms. Comparison of different fluxes for homogeneous turbulence can be found in [17,24]. A different option is to modify the viscous terms to enhance its dissipative properties. The latter has been proposed in [12] using an increased penalty parameter (compared to the minimum required to ensure coercivity of the scheme) when discretising the viscous terms using a interior penalty formulation.

Compressible DGSEM solver
The compressible solver uses conservative variables to solve the Navier-Stokes equations. We use a particular nodal variant of DG methods: the Discontinuous Galerkin Spectral Element Method (DGSEM), see for exam-ple [8]. In addition, the compressible formulation is modified to be energy preserving [20]. The required split-form necessitate Gauss-Lobatto points to cancel out boundary terms using the summation-by-parts simultaneousapproximation-term property (SBP-SAT). The interested reader is referred to [19,9,29,20]. These energy conserving schemes are designed to remain stable and energy conserving and consequently do not necessitate additional localised numerical dissipation. Nonetheless, in this work we introduce dissipation through Roe fluxes, to enhance robustness at high Reynolds numbers. Additionally, viscous terms are discretised using the Bassi-Rebay 1 (BR1) scheme, which is equivalent to the interior penalty formulation when using Gauss-Lobatto points and hexahedral elements [30]. Let us note that this formulation for the viscous fluxes is neutrally stable [21] and adds the minimum dissipation required to achieve a stable scheme, whilst others may introduce some extra dissipation. Other techniques are available to discretise second order derivatives and can be found in the classic review by Arnold et al. [3].

Incompressible DG-Fourier solver
Flow solutions of the incompressible Navier-Stokes equations, are obtained from the 3D unsteady high order h/p Discontinuous Galerkin-Fourier solver [16,11,15,14,12]. The solver uses a second order stiffly stable approach to discretise the NS equations in time whilst spatial discretisation is provided by the discontinuous Galerkin-Symmetric Interior Penalty formulation with modal basis functions in the x-y plane. Here, x represents the streamwise flow direction and y is the normal direction. Spatial discretisation in the z-direction (here defining the spanwise airfoil length) is provided by a purely spectral method that uses Fourier series and allows computation of spanwise periodic three-dimensional flows. Since high order methods (e.g. discontinuous Galerkin and Fourier) are unable to provide enough numerical dissipation to enable under-resolved high Reynolds computations (e.g. as necessary in Large Eddy Simulations), we have adapted the original laminar version of the solver to increase (controllably) the dissipation and enhance the stability in under-resolved simulations [12]. This dissipative formulation has minimal impact on well resolved flow regions and its implicit treatment does not restrict the use of relatively large time steps, thus providing an efficient stabilization mechanism for Large Eddy Simulations. The solver has been widely validated for a variety of flows, including bluff body flows, airfoil and blade aerodynamics and vertical axis turbines under static and rotating conditions [16,11,15,14,12,28,13].

Re = 1 × 10 6 and various Angles of Attack
We start by illustrating the highest Reynolds number case, which is the most challenging in terms of stability and robustness. To illustrate the range of the flow behaviour at various AoAs, we show in figure 2, velocity contours for AoA: 0 • ,5 • and 10 • , computed using the incompressible DG solver. It can be seen that at Re = 1 × 10 6 the flow remains attached for all angles, and that only mild separation is seen near the trailing edge. We will see in the next section that at lower Reynolds numbers this is not necessarily the case.

AoA=5 • and various Reynolds numbers
Having shown the overall good performance in terms of aerodynamic quantities at the most challenging Reynolds numbers, we now focus our attention on the angle AoA=5 • and compare the usability of the solvers to study the NACA0012 boundary layer evolution.
First, we compare the aerodynamic coefficients for AoA=5 • , and Reynolds numbers Re = 1 × 10 5 and Re = 1 × 10 6 , using the incompressible and compressible solvers, both with polynomial order P=3 and P=4, in table 1. We observe good agreement for the highest polynomial order. Small discrepancies  Table 1 NACA0012 airfoil at AoA=5 • for Re = 1 × 10 5 and Re = 1 × 10 6 . Comparison of Lift and Drag using the DG compressible and DG incompressible solvers and two polynomial orders P=3 and P=4.
boundary layer using both solvers in figure 4. It can be seen that detachment near the trailing edge is similar for both solvers. Regarding transition to turbulence (represented by fluctuations in velocity contour), both solvers capture transition on the suction side. The compressible solver shows a transition location near the maximum thickness (x/c ≈ 0.4), whilst the incompressible solver shows transition closer to the leading edge (x/c ≈ 0.2). We have observed significant variations of the transition location for the compressible solver when varying the polynomial order, that we have not seen in the incompressible solver. Further studies are necessary to completely assess the influence of discretisation in the transition location for the two solvers.
Second, we explore the pressure coefficient distribution along the airfoil profile when varying the Reynolds number. We only depict results for the incompressible DG solver since these are very similar to the results provided by the compressible solver. Note that this is not surprising, since the lift coefficients at Re = 1 × 10 5 and Re = 1 × 10 6 are very similar for P=4 at AoA = 5 • , see table 1. Figure 5 shows velocity contours for Re = 1 × 10 4 , Re = 1 × 10 5 and Re = 1 × 10 6 at AoA = 5 • . It can be seen that for the lowest Reynolds, the boundary layer remains laminar until it detaches after the maximum thickness, showing a highly unsteady wake. When the Reynolds number increases, the boundary layer shows transition to turbulence before the maximum thickness, as appreciated by the fluctuations and small scales appearing in figure 5.
To quantify these results, we depict in figure 6, the pressure distribution (Cp) for the three Reynolds numbers. In the top row, we show instantaneous Cp against averaged for incompressible DG solver. In the bottom row, we compare mean Cp distributions against Xfoil [10] (with critical N-factor N cr = 1) and Fluent SST (fully turbulent simulation) [1]. At Re = 1 × 10 4 , the top figure shows that the boundary layer detaches before transition occurs and after the maximum thickness, as shown by the velocity contours in figure  5. Since the flow detaches leading to a highly unsteady wake, there is little hope that the averaged Cp captures the actual behaviour of the boundary layer. This is why, in the bottom figure, the mean values obtained using the incompressible DG solver do not agree with the mean Xfoil and Fluent values that assume steady turbulent flow. At Re = 1 × 10 5 and At Re = 1 × 10 6 , the instantaneous Cp values (top row) show scattering in the data associated to transition. This occurs close to the leading edge on the suction side, whilst it is delayed towards the trailing edge on the pressure side. The bottom row shows that the DG results compare very well to Xfoil when using a critical N=1 (to set the transition point close to the leading edge), whilst Fluent SST (fully turbulent) shows lower Cp values associated to simulating the complete boundary layer as turbulent (no laminar region). This results suggest that DG solvers using iLES approaches (compressible and incompressible) can capture transitional behaviour in boundary layers even when relatively coarse meshes are selected.

Conclusions
In this contribution, we have presented results for turbulent flows over a NACA0012 airfoil. High order discontinuous Galerkin formulations require localised dissipation to remain stable for under-resolved turbulent flow conditions, often referred to as implicit Large Eddy Simulations. Here we have presented compressible and an incompressible DG formulations (with different stabilising mechanisms) that are able to cope with high Reynolds number flows. Both DG formulations provide aerodynamic coefficients and boundary layer information that compare favorably to experimental data and well established low order solvers. We conclude that the compressible and incompressible formulations included in this work can be very useful in aeronautical applications.