High Fidelity Numerical Characterization of NACA 23012 at Re = 300, 000 through a Discontinuous Galerkin LES Solver

1. Introduction

Rotary wing aircraft are characterised by a very complex aerodynamics. In fact, the motion of the lifting surfaces produces complex wakes, which shape is strongly dependent on the flying condition, see e.g. [1,2]. The presence of wakes with complicated geometries generates a plenty of aerodynamic interactions between the various components of the aircraft. According to the flying condition, a relevant interaction that may affects all these possible configurations is the result of the collision of the tip vortex generated by a blade against another blade of the rotor. This phenomenon, said blade vortex interaction (BVI), plays an important role in modern rotary wing aircraft design. In fact, BVI occurrences both generate meaningful acoustic waves and sudden local variations of the lift produced by the stricken blade [3]. Among the several different ways a BVI event can occur a meaningful one is the so–called parallel BVI, where the axis of the tip vortex is locally parallel to the spanwise direction of the blade. While other kinds of BVI events involve almost a single slice of the blade, the parallel BVI acts on a larger portion of the blade due to the peculiar alignment between the vortex and the blade axes. Then, parallel BVI potentially induces larger variations in the overall aerodynamic loads. Recently, there appear two interesting studies concerning parallel BVI events for a NACA 23012 airfoil at $Re=3·{10}^5$: a high fidelity numerical study at $\alpha ={10}^{\circ }$ [4] and an extensive experimental analysis [5]. In particular, this experimental study highlights how the parallel BVI is capable to generate a meaningful variation of the pressure distribution on the suction side, that may triggers the stall of the airfoil when it operates very near to the stall condition.

This experimental investigation provides an interesting global description of the phenomenon. However, the full detail of the intrinsic phenomenology that rules the blade vortex interaction is hardly ever achievable by means of a wind tunnel test. Of course, experimental data are essential for the comprehension of this aerodynamic interaction as for the validation of any numerical tool, but it is evident that an accurate numerically reconstructed flow field easily provides information that are otherwise barely obtainable experimentally. For example, this is the case of the skin friction coefficient distribution over the airfoil and of the details of the interaction between the impinging vortex and the turbulent structures inside the suction side boundary layer. This surely motivates the interest towards high fidelity numerical simulations, but it is not the only reason. In fact, even if namely $2D$, every wind tunnel experiment is not completely two-dimensional, especially near the stall. As important implication data are, up to a certain level, wind tunnel dependent. This is rather evident by the collection of $2D$$C_L\left(\alpha \right)$ curves in Figure 1 for a NACA 23012 airfoil in clean configuration and $Re$ in the range $3·{10}^5$ – $5.5·{10}^5$. The data scattering near the stall condition surely exceeds whatever Reynolds dependence, that is expected small inside the considered range. Likely, after proper validation, high fidelity numerical studies are also more portable between different research groups.

The present paper provides a high fidelity numerical characterization of the clean configuration of NACA 23012 airfoil at $Re=3·{10}^5$ for $5^{\circ }\le \alpha \le {17}^{\circ }$, with a peculiar focus for $\alpha \ge {10}^{\circ }$, i.e., around the stall regime. This study is intended as the first step towards the numerical simulation of parallel BVI events, that would be a valuable addition to the work [5]. Note that such future perspective is not redundant with the study [4], since it would be performed at higher and, in light of [5], more meaningful angles of attack. This final goal determines the selected airfoil and operative conditions. However, this work hopefully is more general, providing a valid numerical reference for the clean NACA 23012 for other researchers interested to any other possible study on this airfoil. Reference that, up to the authors knowledge, is not already available for NACA 23012 at these, or comparable, operative conditions.

2. Materials and Methods

In order to ensure high accurate results, a LES approach has been selected. The present simulations are run using the local discontinuous Galerkin (LDG) LES solver implemented inside the open source FEMilaro library [11], written according to the latest Fortran and MPI standards. This library has been selected by virtue of the highly parallelizable nature of LDG schemes [12], that allow very efficient implementations.

2.1. Mathematical Model

The selected solver approaches the Navier–Stokes equations in their compressible form. The mathematical formulation behind this solver is well reported in [13] but, for the sake of completeness, its main aspects are briefly recalled here. In particular, according to a common practice in compressible LES, equations are filtered supposing the commutation between the LES filter and the derivatives. Indicating $\overline{·}$ as the LES filter and $\widetilde{·}$ the Favre decomposition associated to $\overline{·}$ and defined as $\tilde{f}=\overline{\rho f}/\overline{\rho }$ [14], the compact form (1) for the filtered Navier–Stokes equation is obtained.

$$\begin{array}{cc}\hfill & {\partial }_t\mathbf{U}+\nabla ·{\mathbf{F}}^e\left(\mathbf{U}\right)=\nabla ·{\mathbf{F}}^v(\mathbf{U},\mathcal{G})-\nabla ·{\mathbf{F}}^{SGS}(\mathbf{U},\mathcal{G})\hfill \\ \hfill & \mathcal{G}-\nabla \Phi =0\hfill \end{array}$$

(1)

Here, $\mathbf{U}=[\overline{\rho },\overline{\rho }{\tilde{e}}^{tot},\overline{\rho }\tilde{\mathbf{u}}]$, $e^{tot}=e+\left|\right|\tilde{\mathbf{u}}{\left|\right|}^2/2$ and $\Phi =[\tilde{T},\tilde{\mathbf{u}}]$, where $\rho$ is the density, $\mathbf{u}$ the velocity, e the internal energy and T the temperature. The notation ${\partial }_t$ and ${\partial }_i$ means $\partial /\partial t$ and $\partial /\partial x_i$ respectively. Note that the system (1) is written as a system of first order PDEs, since the LDG discretization requires this kind of rewriting [15,16].

The terms ${\mathbf{F}}^e$ and ${\mathbf{F}}^v$ are the convective (or eulerian) and the diffusive (or viscous) fluxes classically defined in the Navier-Stokes equations, but evaluated with filtered variables. Instead, ${\mathbf{F}}^{SGS}$ is the subgrid flux and contains all the terms that requires a model. In particular, its expression is reported in equation (2).

$${\mathbf{F}}^{SGS}=\left[\begin{array}{c}0\\ {\tau }^{SGS}\\ \frac{1}{k}{\mathbf{Q}}^{SGS}+\frac{\gamma M^2}{2}({\mathbf{J}}^{SGS}-{\tau }_{kk}^{SGS}\tilde{\mathbf{u}})\end{array}\right]$$

(2)

In equation (2), ${\tau }^{SGS}$ is the subgrid stress tensor, ${\mathbf{Q}}^{SGS}$ the subgrid heat flux and ${\mathbf{J}}^{SGS}$ the subgrid turbulent diffusion flux.

Among the different closures proposed inside FEMilaro, the anisotropic dynamic model originally formulated in [17] and extended to the compressible case in [13] is adopted to model the term (2). The details of the mathematical formulation are not reported here since they are not the main focus of the present work, but a couple of comments on the properties of this model are appropriate. In particular, this turbulent closure is both anisotropic and based on a revisitation of Germano’s dynamic procedure. Then, it is capable to capture small scales anisotropicity and it is self–tailored on the local characteristics of the flow, being able to describe also backscatter phenomena. In light of these theoretic features and of the validation with FEMilaro library in [13,18], this model is assumed accurate enough for the present study. Note that the current implementation of the model presents a small change with respect to its original formulation. Indeed, the filter size is not explicitly defined but it is included in the model coefficients determined by the dynamic procedure, as in [19].

2.2. Simulation Setup

The computational domain is obtained from the extrusion in the z direction for a distance $d_z$ of a $2D$ C–type domain defined in the $xy$ plane. The two–dimensional geometry is reported in Figure 2 and it is parameterized in R frontal radius and L wake length. The airfoil trailing edge is located at $(1,0)$ and it is put at the angle of attack $\alpha$ through a rotation around this point. The airfoil chord is $c=1$. Parameters R, L and $d_z$ are selected combining the results of some fast preliminary independence tests, which results are reported in Section 3.1, with available literature concerning high fidelity numerical simulations of $2D$ airfoils, e.g., [20,21,22,23]. The final choice is reported in Table 1.

On this domain, a hybrid structured–unstructured mesh of tetrahedra is build through the open source mesher GMSH. The structure portion around the airfoil is created building a structured $3D$ of hexahedra, composed by 25 layers of 225 elements in the $xy$ plane and repeated 8 times in the z direction. Each of these hexahedra is divided into 6 tetrahedra based on the hexahedron’s vertices. The resulting overall grid consists in 365,938 thetraedra. The construction strategy is identical for all the angles of attack considered in this work. A sketch of the structured portion is reported in Figure 3 for $\alpha ={15}^{\circ }$.

Concerning boundary conditions, a Dirichlet datum equal to the free stream values is enforced on all the outer boundaries except on the sides perpendicular to the spanwise direction, where periodic conditions are imposed. The airfoil surface is modeled as an isothermal wall with no slip condition on the velocity. Instead, initial conditions are the default ones of FEMilaro solver, i.e., a uniform flow equal to ${\mathbf{u}}_{\infty }$ and damped around the airfoil according to a Gaussian function, in order to ensure the no slip condition still from the first time steps.

2.3. Simulation Strategy

The simulation strategy consist in letting the simulation to evolve from the initial conditions to a statistically steady state regime. Then, this state is observed over a statistically meaningful time interval and the main quantities of interest are acquired. In this way, data mean tangential and fluctuating velocity profiles, skin friction coefficient $C_F$ and aerodynamic coefficients are collected.

By a technical viewpoint, all the simulations of this work are performed with a uniform polynomial distribution of degree $q=2$. This results into an overall number of degree of freedom of $3,659,380$ per simulation. Moreover, a time step $dt=9·{10}^{-6}$ is employed to satisfy the $CFL$ condition in the whole domain for the entire duration of the simulations. This choice guarantees a maximum Courant number of $0.85$.

Simulations are run exploiting 2048 parallel threads on AMD® Zen 2 EPYC™ 7H12 CPUs on Karolina, MeluXina and Discoverer supercomputers belonging to the EuroHPC JU consortium. In details, runs on Karolina and MeluXina clusters are executed using 2048 physical cores, while computations on Discoverer cluster take advantage of hyperthreading, running 2048 processes on 1024 physical cores without any meaningful time overhead.

3. Results

3.1. Preliminary Tests

The preliminary tests consists in the comparison between the first time instant of a set of ILES simulation at $\alpha ={15}^{\circ }$ and variable R. This peculiar angle of attack is selected since it is one of the points associated to the largest uncertainty in terms of $C_L$ and $C_D$, as emerged looking at Figure 1. Then, this condition is the one requiring the best possible simulation.

These tests are performed on the time interval $t\in [0,5]$, that it is not enough to reach a statistically steady state regime, as evidently reported from Figure 4.

This is a conscious choice, since the aim is to perform computationally cheap preliminary tests, but requires a careful comparison between the data. The absence of a final state implies that there is no regime value to compare across the tests. To avoid the comparison of pointwise values, that are rather meaningless during a transient with turbulence, the equivalence of the time histories is evaluated comparing the mean value of the difference between the time histories. Moreover, since a compressible solver is used, data must be compared after the pressure wave generated by the boundary conditions at $t=0$ reaches the airfoil. Of course, the change in the grid size implies boundary conditions are felt for the first time at different time instants in the different simulations. This fact is well depicted in Figure 4: after the violent oscillations in the lift coefficient, expected from the inviscid theory and that represents the first phase of the airfoil startup [24], there appear a negative peak generated by the boundary pressure wave at a location that is R dependent. Then, data must be compared at least after $t\ge 2$. Results of the comparison are reported in Table 2 for $t\ge 2.5$.

Note that in these tests $L=2R$ is supposed. The comparison with literature data [20,21,22,23] suggests a a shortening, then $L=17$ is selected.

3.2. Simulations

In the present work six different angles of attack are considered: five around the stall regime, that are ${10}^{\circ }$, ${13}^{\circ }$, ${15}^{\circ }$, ${16}^{\circ }$, ${17}^{\circ }$, and one in the linear range of the $C_L\left(\alpha \right)$ curve, at $5^{\circ }$. This latter case might appear out of the scope of this work. However, the purpose of this low $\alpha$ simulation consists in offering a more robust validation of the results.

All the simulations are run according to the simulation strategy proposed in Section 2.3. Note that each single run requires a small ad hoc customization. In fact, the growing flow unsteadiness due to the increase in the angle of attack makes the observation interval of the final statistically steady state simulation-dependent. Simulations at the lowest values of $\alpha$ require a much shorted acquisition window. Instead, the settling time to reach this final regime is roughly constant for all the tested angles.

The instantaneous flow fields in the steady state regime are reported in Figure 5 according to the angle of attack. At $\alpha =5^{\circ }$ the flow is completely attached and exhibits a thin boundary layer with a low turbulence level. Increasing the angle of attack, a larger and larger low velocity region appears at the trailing edge from $\alpha ={15}^{\circ }$, indicating an incoming trailing edge stall. The same information is deduced by the progressive thickening of the boundary layer. Moreover, as evident by the leading edge detail in Figure 6, the transition point moves progressively towards the front portion of the airfoil. There, the transition is triggered by a laminar separation bubble (LSB), especially at high values of $\alpha$.

The qualitative description deduced from the the flow fields is supported from the analysis of the velocity profiles in Figure 7. Figure 7 reports the mean tangential velocity profiles extracted along lines normal to the airfoil surface. The chord position of these lines is reported in Table 3, while Figure 8 offers a their sketch at $\alpha ={15}^{\circ }$.

The analysis of Figure 7 confirms a progressively larger mean trailing edge separation starting for $\alpha \ge {15}^{\circ }$. Moreover, also the first velocity profile exhibits a recirculation region, confirming the presence of a laminar separation bubble near the leading edge. As already observed in Figure 6, this LSB is stronger higher the angle of attack is. The presence of such aerodynamic structure is expected at these operative conditions [25] and also low fidelity simulations through solvers as XFOIL exhibit a LSB for NACA 23012 at $Re=3·{10}^5$. The presence of the laminar separation bubble, as well as of the trailing edge separation, is confirmed by the skin friction coefficient $C_F$ reported in Figure 9. In fact, the are negative regions of $C_F$ both on the front and at the rear of the airfoil. These regions are larger higher $\alpha$ is, confirming what observed about profiles 1 and 10.

It is also of interest the analysis of some turbulent quantities, as the resolved turbulent fluctuations. Then, some profiles for the fluctuating quantities $<u_i^{\prime }{u}_j^{\prime }>$ are investigated, where $u_i^{\prime }={\overline{u}}_i-<{\overline{u}}_i>$ and $<·>$ denotes the time average. With reference to the same chord positions used for the mean velocity profiles and reported in Table 3, the profiles for $<u_1^{\prime }{u}_1^{\prime }>$, $<u_2^{\prime }{u}_2^{\prime }>$ and $<u_1^{\prime }{u}_2^{\prime }>$ are reported in Figure 10 for the angles of attack in Figure 9. Also the lower value of $\alpha ={10}^{\circ }$ is added to provide a better comparison.

At the position of the laminar separation bubble the highest fluctuations are observed both in terms of $<u_1^{\prime }{u}_1^{\prime }>$ and $<u_2^{\prime }{u}_2^{\prime }>$, even if they are limited in a small region near the wall. Also the largest negative value of $<u_1^{\prime }{u}_2^{\prime }>$ is observed at this location. Note that a negative value of $<u_1^{\prime }{u}_2^{\prime }>$ is connected to the production of turbulent kinetic energy. Then, the LSB is effectively triggering transition from the laminar to turbulent boundary layer, generating the most part of the turbulent kinetic energy. Moreover, note how this production is larger higher the angle of attack is. Comparing this information with Figure 9, it is clear that a stronger LSB introduces more turbulent kinetic energy in the boundary layer. Moving along the suction side, the absolute values of the fluctuations are smaller and smaller near the wall, but it is evident how, especially near the trailing edge, these fluctuations are no longer limited to the near wall region. After the LSB the turbulent production progressively decreases, while larger structures progressively appears. In practice, small turbulent structures near the leading edge are transported by the flow towards the leading edge and during this transport redistribution mechanisms make them grow. As expected, higher the angle of attack is stronger this phenomenon is.

Concerning the characterization of the global behavior of NACA 23012, it is provided by means of lift and drag coefficients. Figure 11 presents $C_L\left(\alpha \right)$ and $C_D\left(\alpha \right)$ numerical curves compared with experimental data from [6,7,8,9,10]. Numerical coefficients are averaged in time over an appropriate time interval. Data are summarized in Table 4.

Remarkably, the lift coefficient data demonstrate the presence of a very smooth stall. The stall is a trailing edge separation and it is not triggered by the blow up of the laminar separation bubble. The LSB exhibits a stable behaviour across all the considered angles of attack, suggesting that its presence is only detrimental for the maximum lift coefficient [25] but does not play a role in the stall mechanism in the light stall region.

4. Discussion

In order to validate the present results, the comparison with experimental data in Figure 11 is essential. However, the numerical $C_L\left(\alpha \right)$ curve in Figure 11 agrees only with some experimental sets of measurements. Since the evident scattering of the experimental data this is rather expected, but complicates the validation of the current results.

In the linear range, all experiments almost agree one with the other except the one of Dent et al. [7]. Looking at this peculiar set of data, the $C_L\left(\alpha \right)$ slope is smaller. This reduction likely is related to a $3D$ aerodynamic effect as non negligible tip vortices. Moreover, the presence of tip vortices that reduce the effective angle of attack of the airfoil justifies also the high stall angle, around ${18}^{\circ }$, higher than in all the other measures.

Excluding this set of data, the numerical result at $5^{\circ }$ is in good agreement with all the measurements. Instead, for higher angles, i.e., $\alpha \ge {10}^{\circ }$, the numerical estimations are confirmed only by Ashender et al. [6] and Tenzin [10]. In this case, note that the deviation from Jacobs et al. [8] and Rivero et al. [9] is mainly due to a different behavior outside the linear range. Observe that the deviation from the linear range usually takes place at the point where viscous effects start playing an important role and separation begins occurring. In these conditions, experiments can be affected by the peculiar characteristics of the wind tunnel where they are performed. The larger scales of the flow structures that appear approaching the stall potentially are more sensitive to the finite size of the test chamber. Moreover, the interaction between the boundary layer on the wind tunnel walls and the one on the airfoil alters the flow on the model tips in a way that surely is wind tunnel dependent. Furthermore, this interaction is reasonably larger thicker the airfoil boundary layer is, then near the stall condition. The lower predictions in terms of $C_L\left(\alpha \right)$ and stall angle of Jacobs et al. and Rivero et al. are attributed to wind tunnel effect as the just discussed ones. Since the full detail concerning the test strategies is not available this is only an hypothesis, but a reasonable one.

Then, consider the two remaining measurements. Tenzin data confirm the LES maximum lift coefficient, since the experimental datum is $C_{L,max}^{Tenzin}=1.49$ and the numerical prediction is $C_{L,max}^{LES}=1.51$. However, Tenzin suggests an earlier stall, around $\alpha ={14}^{\circ }$. On the other hand, Ashender reports a similar stall, both in terms of smoothness and angle of maximum lift coefficient, but a lower $C_{L,max}$ of $1.393$. However, this data is provided with a confidence interval of $\pm 0.048$ and the LES estimation is quite near to the upper bound of this interval. In addition, the previous discussion on Jacobs and Rivero data suggests that a slightly lower experimental prediction of $C_{L,max}$ with respect to the perfectly $2D$ and ideal LES environment is reasonable.

In light of these considerations, the good agreement of numerical data with Tenzin and, particularly, with Ashender measures suggests the correctness of the found results. Then, a numerical baseline is available for future studies involving a NACA 23012 airfoil at the operative conditions of $Re=3·{10}^5$. The presence of these accurate numerical flow fields is a valuable addition to experiments, since it provides a reference not only for lift and drag coefficients, but also for the details of the flow physics.