On the dynamics of the turbulent flow past a three-element wing

A comprehensive analysis of the unsteady flow dynamics past the 30P30N three-element high lift wing is performed by means of large eddy simulations at different angles of attack ( a ¼ 5 (cid:2) , 9 (cid:2) , and 23 (cid:2) ) and at a Reynolds number of Re c ¼ 750 000 (based on the nested chord). Results are compared with experimental and numerical investigations, showing a quantitatively good agreement and, thus, proving the reliability and accuracy of the present simulations. Within the slat and main coves, large recirculation bubbles are bounded by shear layers, where the onset of turbulence is triggered by Kelvin – Helmholtz instabilities. In the energy spectrum of the velocity fluctuations, the footprint of these instabilities is detected as a broadband peak; its frequency being moved toward lower values as the angle of attack increases. Kelvin – Helmholtz vortices roll-up and break down into small scales that eventually impinge into the slat and main coves lower surfaces. The slat impingement shows to be more prominent, and hence, larger velocity and pressure fluctuations are observed. The impingement strength diminishes with the angle of attack in both coves, while higher fluctuations are originated on the slat and main respective suction sides, leading to larger boundary layers. This is associated with the displacement of the stagnation point with the angle of attack. Another salient feature observed is the laminar-to-turbulent flow transition in the main and flap leading edges although the average location of this transition seems to not be affected by the angle of attack. Tollmien – Schlichting instabilities precede this transition, with the disturbances amplified by the inviscid mode at low angles of attack, while at a ¼ 23 (cid:2) , the local Reynolds number on the main suction side is incremented and the viscous mode becomes important. The analysis shows that the turbulent wake formed at the trailing edge of all elements dominates the dynamics downstream. This is especially true at the higher angle of attack, where a large region of velocity deficit above the flap is observed, thus indicating the onset of stall conditions.


I. INTRODUCTION
During takeoff and landing operations, aircraft noticeably reduce their air velocity, and the deployment of high-lift devices, i.e., flaps and slats, becomes necessary to compensate for the loss of lift.In the early stages, the aeronautical industry was mainly focused on the aerodynamic performance during the cruise, while the design of high-lift devices only aimed to maintain reasonable takeoff and landing distances.Through the years, the attention has progressively turned to increasing the performance of these elements, especially due to the development of modern numerical methods and advances in highperformance computing.The unsteady flow generated by high-lift devices is very complex, and its analysis was out of reach for the computational capabilities of the early years.Some benefits of optimizing these devices are the increase in the total payload 1 or the reduction in noise emissions, 2 leading to larger operations periods and a lower environmental impact.Therefore, understanding the underlying physics of multi-element wings is fundamental for the design of future modern aircraft.
One of the most common research models to carry out these investigations is the 30P30N configuration.This high-lift wing model, designed by McDonnell-Douglas, was first used for the high-lift CFD challenge workshop held at NASA Langley Research Center, which entailed experimental 3,4 and Reynolds-averaged Navier-Stokes (RANS) numerical investigations. 5The primary objectives of these early studies were to analyze some key aerodynamic parameters and assess the capabilities of the computational tools of the time in predicting the flow around multi-element wings.With the years, some subsequent experiments 6,7 and RANS contributions [8][9][10] were published.However, it was a few years later when Lockard and Choudhari 11 applied more advanced turbulence models for the first time and compared the previous RANS computations with their delayed detached eddy simulations (DDES).Results showed slight differences, with the greater discrepancies occurring in the flap separation region.They also evaluated how the flow evolves as the Mach (M 1 ¼ U 1 =a) and Reynolds (Re c ¼ U 1 c=) numbers are modified, obtaining almost no dependency on the Reynolds number; U 1 being the freestream velocity, a is the speed of sound, c is the nested chord, and is the kinematic viscosity.
During these years, the attention turned to the flow dynamics at the slat cove and the 30P30N geometry became popular among the scientific community for the study of the aeroacoustics in high-lift wings.Investigators were aware that these devices contribute to increasing the airframe noise, in particular the recirculation bubble at the slat cove.This generates a prominent shear layer that emanates from the slat cusp and reattaches further downstream, impinging on the slat's lower surface and producing a series of acoustic peaks.Due to a feedback mechanism analogous to the Rossiter modes 12 in open cavity flows, narrowband peaks at the mid-frequency range appear, while broadband spectral peaks associated with the slat trailing edge vortex shedding are observed at higher frequencies.All these findings were reported in the third workshop on Benchmark problems for Air-frame Noise Computations (BANC-III) held by AIAA, 13 which studied the 30P30N wing model at Re c ¼ 1:7 Â 10 6 , angle of attack a ¼ 5:5 , and M 1 ¼ 0: 17.In this workshop, some experimentalists contributed to increasing the validation databases.This is the case of Jenkins et al., 7 Pascioni et al., 14 and Murayama et al. 15 Later on, Pascioni and Cattafesta 16,17 continued their research and performed more experiments on this wing configuration.In their works, they observed that the strength of the narrowband peaks decreases with the angle of attack and linked them with the shear layer instabilities, matching the spatial wavelength of the coherent structures with the respective frequencies.They also associated a low-frequency hump that was previously observed in other experiments with the shear layer flapping motion and applied a modification of the Rossiter's equations 18 to predict the cavity tones, finding a good agreement with their experiments.Other experimentalists tried to consider slat cove fillers (SCF) to reduce the emitted noise.This is the case of as Jawahar et al. 19 and Zhang et al., 20 who also evaluated the effects of slat extensions.In both studies, SCFs totally removed the narrowband peaks, slightly increasing the overall broadband noise levels.Concerning the slat extensions, they reduced the tonal peaks and displaced them to higher frequencies, but their effectiveness showed to be lower in reducing the emitted noise.
Plenty of numerical studies were presented at BANC-III, the majority employing hybrid RANS/LES models, such as detached eddy simulations (DES), delayed detached eddy simulations (DDES), or improved delayed detached eddy simulations (IDDES).Within the framework of this workshop, several participants published their results, [21][22][23][24][25] while other authors started using the 30P30N configuration to address their investigations, most of them also considering hybrid models.Among them, Shi et al. 26 and Gao et al. 27 employed finite element (FE) high-order schemes to tackle the problem by increasing the degree of the polynomial inside the element, i.e., flux reconstruction and spectral-differences methods, respectively.Both were focused on the slat cove dynamics and employed this geometry as a benchmark case to test the reliability of high-order methods in accurately predicting near-and far-field aeroacoustics, achieving both a qualitatively good agreement with the previous results.Jin et al. 28 also employed higher-order methods using a sixth-order cell-centered finite difference method, extending the aeroacoustics analysis to the main and flap elements as well.In spite of the interesting features contributing to the aeroacoustic noise in the different elements, they found that the major contributor to the emitted noise is the slat, in agreement with previous investigations.Sakai et al. 29 used a reduced dissipation scheme that introduced a c coefficient to control the numerical dissipation and also compared the capabilities of unstructured grids vs structured ones.They detected that the turbulent kinetic energy becomes larger at earlier stages of the shear layer, which could yield better resolved velocity fluctuations.To test the capabilities of their methods, Ueno et al. 30 considered an adaptive mesh refinement method to increase grid density in places of interest, and Kojima et al. 31 proposed an embedded-LES (ELES) approach as a possible substitute for the classical hybrid RANS/LES methods.In their approach, pseudoturbulence inside the turbulent boundary layer is generated, which promotes three-dimensionalities at earlier stages in the shear layer compared to DDES, although no significant differences were encountered in the slat wake.Shur et al. 32 also tested a new experimental hybrid approach named advanced detached eddy simulation (ADES) characterized by user-specified RANS and LES regions.The goal was to explore the sensibility of the method to the user's selection and, despite the method showed to be robust, the automatization of the selection seemed difficult and far from being prepared for industrial applications.
Away from these hybrid approaches, only Bodart et al. 21and Zhang et al. 33 performed wall-model large eddy simulations (LES).The former conducted simulations at sixteen different angles of attack up to stall conditions and mainly focused on the aerodynamic curves, i.e., lift and drag coefficients at the different angles.They pointed out that stall is not originated from the separation of the boundary layer from the surface of the wing but from the separation of the shear layer between the main and flap boundary layers.The latter focused on the slat aeroacoustics.In their research, they confirmed that the slat cove is the primary sound source, while the main element dominates the lowfrequency spectra and the noise of the flap is negligible.They also showed that the acoustic intensity scales as M 5 1 .All previous investigations are conducted at moderate to high Reynolds number and mainly focused on the aeroacoustics.Studies at lower Reynolds numbers can also be found in the literature.In the experimental field, Wang et al. 34,35 studied the formation of G€ ortler vortices on the slat wake in the low-Reynolds regime.They found that these vortices are originated at medium angles of attack, between a ¼ 2 and 12 , due to the recirculation bubble formed in the slat cove.They detected that G€ ortler vortices appear below a critical Reynolds number when no roll-ups of the slat shear layer occur.Above the critical value, roll-ups are formed in the shear layer and the slat wake is then dominated by streaks and spanwise vortices instead.Similar findings were reported by Vadsola et al., 36 who performed direct numerical simulations (DNS) at similar Reynolds numbers.Long et al. 37 also carried out experiments at low Reynolds numbers and observed that the entrainment flux in the transitional boundary layer of the main element is controlled by the length of the interface, suggesting that this could be a universal phenomenon.
In addition to the 30P30N wing, numerical simulations on the flow past other three-element high-lift wings are available in the literature.Examples of these are the investigations of Deck and Laraufie 38 who, by means of a zonal detached eddy simulation (ZDES) at Re c ¼ 2:09 Â 10 6 on the DLR-F15 wing configuration, presented a thorough description of the flow field around all elements.Similar slat cove flow dynamics as those of the 30P30N configuration were reported, and the aeroacoustics analysis was extended to the main cove cavity as well.In addition to the self-sustained mechanism in each cavity, they further pointed out the existence of a large-scale coupling between both coves which may lead to a more intense sound.Another noteworthy contribution was carried out by Terracol and Manoha 39 who conducted large eddy simulations (LES) on the flow past the DLR-F16 wing at Re c ¼ 1:23 Â 10 6 .The simulation further confirmed the existence of a large-scale coupling mechanism between the slat and flap cavities and provided some insights into the main laminar-to-turbulent transition process.
A summary of the most relevant experimental and numerical contributions found in the literature at moderate Reynolds numbers is provided in Tables I and II, respectively.All in all, although a wealth of literature exists on the subject, most studies either concentrate on general aerodynamic coefficients, on the slat cove dynamics and aeroacoustics, or serve as test cases for novel numerical methodologies.Moreover, these studies employ hybrid RANS/LES models or wallmodeled LES, and no wall-resolved LES can be found in the literature concerning the 30P30N configuration, which clearly manifests the necessity of high-fidelity predictions.Thus, considering the state-ofthe-art, the present work aims to comprehensively investigate and analyze the flow dynamics in a high-lift wing at different flow conditions so as to deepen into the underlying physics, which could serve to identify potential strategies to optimize the performance of wings.To achieve so, wall-resolved LES of the flow past the 30P30N high-lift airfoil are performed at a Reynolds number of Re c ¼ 750 000 and different angles of attack a ¼ 5 , 9 , and 23 .

II. COMPUTATIONAL CONFIGURATION
The geometry considered in the present study is the 30P30N wing (Fig. 1), which embodies an unswept high-lift three-element wing with the slat (d s ) and flap (d f ) deflected 30 .All measurements are non-dimensionalized with the stowed chord (c) of the airfoil; the slat (s) and flap (f) chords being 0:15c and 0:3c, respectively.For the flow properties, the Reynolds number based on the nested chord is set to Re c ¼ 750 000, and three angles of attack a ¼ 5 , 9 , and 23 are contemplated.The first two angles correspond to pre-stall conditions, while the latter represents a post-stall flow regime. 21As the focus of the present research is to perform wall-resolved LES, the Reynolds number used here is lower than in most of the experiments.Based on Choi and Moin 40 estimates, the computational requirements to capture all boundary layer turbulent structures scales as Re 13=7 .For this reason, to make simulations computationally more affordable, the Reynolds number has been slightly scaled down compared to the previous experiments without loss of generality in the results.

A. Computational methods
The incompressible spatially filtered Navier-Stokes equations in Einstein's notation summation can be written as where u i with i ¼ 1, 2, 3 (or u, v, w) and p are the filtered velocity and pressure fields, and q are the fluid kinematic viscosity and density, and s ij is the subgrid-scale (SGS) stress tensor.This last term encloses the effect of the non-resolved scales and, hence, requires to be modeled.Its deviatoric part is modeled via an eddy viscosity ( SGS ) model, where S ij ¼ 1=2ð@u i =@x j þ @u j =@x i Þ is the rate-of-strain tensor of the resolved velocity field and d ij is the Kronecker delta.sgs is here modeled using the Vreman model. 41he above system of equations is solved using Alya, 42 a lowdissipation finite element simulation code developed at the Barcelona Supercomputing Center (BSC).The convective term is discretized using an energy, momentum, and angular momentum conserving (EMAC) scheme initially developed by Charnyi et al. 43 and later extended for fully turbulent flows by Lehmkuhl et al. 44 Neither upwinding nor any equivalent momentum stabilization is employed.Only a non-incremental fractional-step method 45 is used to stabilize pressure, which allows for the use of finite element pairs that do not satisfy the inf-sup conditions, such as equal order interpolations for the velocity and pressure.For the linear finite elements considered here, the spatial discretization is second-order, which is equivalent to the error obtained for the pressure-velocity coupling in unstructured collocated finite volume codes.For the time integration, an explicit energy-conserving second-order Runge-Kutta method is applied that, despite being of second order, preserves kinetic energy to a higher order through the application of pseudo-symplectic methods. 46This is combined with an eigenvalue-based time step estimator. 47In this approach, the eigenvalues of the dynamical system are analytically bounded and the linear stability domain of the time-integration scheme is adapted in order to maximize the time step.In any case, in the present computations, the time step at each iteration is kept below CFL < 1, this being the Courant-Friedrichs-Lewy number.For more details about the present methodology and numerical tests, the reader is referred to Lehmkuhl et al. 44 The computational setup presented in this section has been successfully applied in scale-resolving simulations with similar configurations, e.g., the flow past a SD7003 airfoil considering active flow control (AFC) devices, 48 the flow past a circular cylinder at high Reynolds numbers, 49 or the flow past a sphere. 50

B. Computational mesh
The computational domain employed in the present simulations is presented in Fig. 2. In the x-y plane, the two-dimensional geometry of the airfoil is positioned at the center of a circular domain that extends radially a distance of r ¼ 10c away from the wing walls.Similar domain sizes are successfully applied by other authors in the literature. 9,27,33Along the spanwise direction z, the model is extruded a distance of L z ¼ 0:1c.This spanwise size is adequate to obtain a proper spanwise decorrelation of the slat cove flow with periodic boundary conditions, as demonstrated by Lockard and Choudhari. 10egarding the boundary conditions, a uniform velocity profile is imposed at the inlet boundary, with the velocity components depending on the angle of attack as ½u; v; w 1 ¼ U 1 ½cosðaÞ; sinðaÞ; 0; zerogradients are imposed at the outlet regions; the no-slip condition is applied on the airfoil walls; and periodic conditions are considered in the z direction.These boundary conditions are also summarized in Fig. 2.
For all the simulations, a hybrid unstructured mesh is generated with structured-like inflation layers around the airfoil walls (see Fig. 3).To correctly represent the structures in the near-wall region, a nondimensional wall-normal distance of Dg þ < 1 is considered in the wall-bounded turbulent zones, while the wall-tangential distance is kept below Dn þ < 80.Only in the coves, where recirculation bubbles are found, the wall-normal spacing reaches values up to Dg þ ¼ 30.On the other hand, the spanwise direction z is discretized using 129 grid points, which translates into a non-dimensional distance of Dz þ < 50.Note that Dg is defined as the wall-normal distance of the first off-wall node, whereas Dn and Dz are the wall-tangential and spanwise spacings of the wall nodes, respectively.According to Piomelli and Chavnov, 51 the spacings employed in the present computations are enough to represent the near-wall structures.All these quantities are non-dimensionalized using the scaling factor u s = and are denoted with the superscript þ, where u s ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi js w j=q p is the friction velocity and s w represents the wall shear stress.Using the local coordinate system n À g, the wall shear stress is computed as the wallnormal derivative of the tangential velocity evaluated in the wall, s w ¼ l ðdu n =dgÞ g¼0 .In regions of separated flow, the mean spacings are estimated to be 15g K , where g K is the Kolmogorov length scale.According to Pope, 52 motions responsible for dissipation are in the range of 8 < g K < 60, the peak being at g K % 24.Considering this, the mesh used in the present computations should be capable of capturing the main flow turbulent features.However, since the present work does not aim to resolve the far-field wake, downstream the flap (at about 1c) the mesh is rapidly coarsened.All this yields a grid count of 449, 000 nodes within the x-y plane and 58 Â 10 6 nodes throughout the whole domain.To further assess the mesh resolution, a grid independence analysis is performed in Sec.III A. Generally, it is shown that the current mesh is capable of capturing the main flow parameters with a good compromise between accuracy and computational time to solve the current simulations.

III. RESULTS
To save computational time, all simulations are started from an initial map obtained with a coarser mesh of about 6 Â 10 6 grid points.These computations are performed for about 20 U 1 =c.Once the statistical steady state is reached, the resulting velocity field is interpolated into the mesh described in Sec.II B and simulations are relaunched.After a transient period of approximately 3 U 1 =c, time-average statistics are obtained by the continuous integration of the instantaneous flow over up to 14 U 1 =c.A similar integration time was used by Terracol and Manoha 39 in a similar three-element wing configuration.For the presentation of the results, all quantities are nondimensionalized using the reference length, velocity, and time as c, U 1 , and c=U 1 .

A. Grid independence assessment
To assess the grid independence over the results, a finer mesh is created.In this mesh, while Dg þ is maintained practically the same, the streamwise spacings on the upper surfaces are reduced at the expense of coarsening the lower ones.This allows to reach Dn þ < 50 without practically changing the number of grid points in the x-y plane.The discretization along the z direction is increased from 129 to 225 points, which gives a spanwise resolution of Dz þ < 30.All in all, the total grid count is increased from 58 to 98 Â 10 6 nodes.Note that the finer mesh is only employed for validation purposes at Re c ¼ 750 000 and a ¼ 5 .All the mesh details can be consulted in Table III, where N refers to the number of grid points.In Table IV, the lift C l ¼ l=ð0:5qSU 2 1 Þ and drag C d ¼ d= ð0:5qSU 2 1 Þ coefficients predicted by each mesh are compared, as well as the root mean square (rms) values of these parameters; l is the lift force, d is the drag force, and S is the reference surface (S ¼ L z c).This table shows that the difference in C l and C d is roughly 1% and 3%, respectively, with a relative error of approximately 1% and 2% in the rms values.These small differences are a good indicator of gridindependent results.To further assess this, Fig. 4 displays the pressure C p ¼ ðp À p 1 Þ=ð0:5qU 2 1 Þ and skin friction C f ¼ s w =ð0:5qU 2 1 Þ coefficients at the wing surface; p 1 being the freestream pressure.As can be seen from the figure, both meshes exhibit a good agreement and only some discrepancies are found in the transition zone at the main leading edge (LE).In the remaining areas, the collapse of both graphs is almost perfect.

Designation
In Fig. 5, the mean tangential velocity profiles (u n =U 1 ) at different chordwise locations are compared.In general, a quite good match is observed at all stations.Along the main suction side, slight discrepancies are found within the boundary layer, while the edge values are kept the same in the two meshes.In the flap, the fine mesh exhibits reduced velocities in the wake of the main and also shows an earlier onset of flow separation although the discrepancies are roughly negligible.At the same locations as for the velocity, in Fig. 6, the turbulent kinetic energy (TKE=U 2 1 ) profiles are evaluated.Since this is a second-order statistic, higher sensitivity on the mesh resolution is expected.Despite this, the results predicted by both meshes collapse reasonably well.The main differences are found in the peak magnitudes, where the fine mesh exhibits lower values.Again, an anticipated separation of the flow is observed on the TKE near the flap trailing edge (TE).The sensitivity of the flow over the flap to the numerical conditions, especially on the separation phenomena, was also highlighted by other authors 38,39 with similar high-lift wing configurations.
Overall, the quantitatively small differences encountered in the two meshes yield the conclusion that both predict similar results.Hence, considering the computational cost for the finer grid, all the results presented in this work have been obtained with the baseline mesh.

B. Validation
In this section, the present LES computations are compared against experimental results from FSAT, 14 JAXA, 15 and LTPT 3,4 facilities.In those experiments, the Reynolds number was higher (Re c ¼ 1:71 Â 10 6 and Re c ¼ 5 Â 10 6 ) than the value under consideration in this work (Re c ¼ 0:75 Â 10 6 ).Nevertheless, as argued by Lissaman, 53 the aerodynamic behavior of airfoils above a given Reynolds number starts to depend weakly on this parameter when the transition to turbulence due to the incremented Reynolds number prevents the formation of laminar separation bubbles.This is confirmed by the small differences encountered by Lockard and Choudhari 11 in their numerical simulations considering different Reynolds numbers, or comparing the results in Fig. 7(c) of JAXA and LTPT experiments.In the case of the C f , this can be more sensible to the variations of the Reynolds number, as shown by Vinuesa et al. 54 Nevertheless, despite variations in the Reynolds number, a similar trend in the behavior of the skin friction should be anticipated.Considering the lack of data available in the literature at similar Reynolds numbers, the comparisons presented in this section, together with the grid independence analysis, should provide enough insights to validate the present study.
These comparisons are depicted in Fig. 7, which shows the lift coefficient C l , along with the pressure coefficient distribution C p at the three angles of attack studied.In light of the differences in the Reynolds number between the present computations and the experiments, the observed agreement is rather good.Indeed, when comparing the contributions of the three elements to the overall lift, there is a good collapse between the LES predictions and the experiments in both the linear, i.e., a ¼ 5 and 9 , and stall regimes, i.e., a ¼ 23 .These observations are consistent with those of Lockard and  Choudhari, 11 who did not detect noticeable differences in the C p when the Reynolds number varied from 1:4 Â 10 6 to 2:4 Â 10 6 .However, when reproducing the experimental results, some differences deserve attention.In these experiments, certain blockage effects arising from the wind tunnel walls are reported.For example, the angles of attack at FSAT (with hard walls) and JAXA differ, with a FSAT À a JAXA % À2 in order to achieve similar C p distributions, and corrections can exceed À3 for high angles of attack. 16Some of the small discrepancies among these experiments, as well as with the present computations, can be attributed to these effects.By applying corrections of a LES % a JAXA þ 1 and a LES % a FSAT À 1 to the JAXA and FSAT measurements, respectively, the computed C l would be closer to the values reported in the experiments.These corrections displace the present C l À a curve to the left and right, respectively.Nevertheless, they were not applied to Fig. 7(a).
For the C p distributions depicted in Figs.7(b)-7(d), it is observed that at a ¼ 5 , the conditions most similar to the experimental setups are those of FSAT at a ¼ 5:5 and JAXA at a ¼ 5 .One might expect that closer C p distributions would be obtained for a ¼ 4 and 6 , respectively.Unfortunately, these data points are not available for comparison.For a ¼ 9 , the experiments of JAXA at a ¼ 8 seem to agree well with the present predictions, while there is a slight disagreement for LTPT measurements at a ¼ 8 .Instead, it is estimated that a ¼ 7:5 would yield a better match, but again, these measurements are not reported.For a ¼ 23 , the only available experiments are those carried out at LTPT considering a ¼ 21 .Here, the effects derived from the differences in the Reynolds number might become important, especially at high angles of attack, where flow conditions are away from the C l À a linear regime.
In Fig. 8, the computed skin friction coefficient is compared with the experimental data from LTPT. 4 Analyzing the C f distributions along the main and flap walls, it can be seen that this parameter is less sensible to a; only small variations are detected at different angles.However, it shows a higher dependence on the Reynolds number, and hence, deviations from experimental values may become significant.It is worth noting that the experiments were conducted at Re c ¼ 5 Â 10 6 , which is considerably higher than the Reynolds number under study.The dependence of the skin friction with the Reynolds number was also reported by Vinuesa et al. 54 for a NACA 4412 airfoil, who observed lower skin friction magnitudes at increasing Reynolds numbers.Nevertheless, the general behavior seems to be well captured by the simulations, and the order of magnitude appears to be in alignment with the experiments.
At this point, results are validated with other numerical simulations.The most computationally studied case in the literature is the one considered at BANC-III, 13 i.e., a ¼ 5:5 and Re c ¼ 1:71 Â 10 6 , which is close to the flow conditions of the present predictions at a ¼ 5 .Figure 9 compares the predicted C l , the root mean square lift coefficient (C l;rms ), and the slat shear layer impingement distance (d imp: ) at a ¼ 5 with the numerical results presented at BANC-III.All three quantities fall within the stipulated confidence intervals and are even close to the mean values of the workshop.Only, the C l;rms exhibits a little deviation from the mean value, which is also the parameter showing the higher scattering within the workshop results.This confluence of results further underscores the reliability of the present simulations.
Finally, Fig. 10 provides a comparative analysis of the velocity magnitude (U=U 1 ) and turbulent kinetic energy (TKE =U 2 1 ) profiles against the experimental results at the FSAT wind tunnel (a ¼ 5:5 and Re c ¼ 1:71 Â 10 6 ) along the L2, L4, and L7 lines, which are established to intersect the slat shear layer at specific locations (see the Appendix).The value d used in the graphs refers to the distance from the origin of the line and is non-dimensionalized with its total distance (d max ).Although these measurements were obtained in Pascioni et al., 14 the experimental dataset employed here has been communicated through personal correspondence.As can be seen from the figure, U=U 1 profiles show good agreement with the experiments, and only some notable differences are found in the extreme locations, i.e., d % 0 and d % 1.The TKE is also in line with the experimental data despite some minor discrepancies found in the peak magnitudes.In general, results collapse well with the experiments, and the small variations might be attributed to different facts.First, differences in flow conditions, particularly in a, might induce a slightly different geometrical distribution of the shear layer and, hence, different intersections.Second, although a strong dependence on the Reynolds number is not expected, it is known that the transition to turbulence in the shear layers and its location might be affected by this parameter. 55,56Thus, if the inception of shear layer instabilities occurs at a different location depending on the Reynolds number, this will, in turn, result in small variations in the profiles of the TKE along the shear layers.Last but  not least, limitations inherent in the experimental particle image velocimetry (PIV) techniques, alongside potential computational errors, can further contribute to the observed variations.In summary, the numerical results presented in this section demonstrate a consistent agreement with both existing experimental and computational data, proving the reliability of the present LES.

C. Results overview
To provide an overview of the flow past the 30P30N wing, the vortical structures at a ¼ 5 are visualized by means of Q-criterion isocontours in Fig. 11.This representation allows for the identification of various flow dynamics inherent to high-lift three-element airfoils.Large recirculation bubbles are observed in the slat and main coves, which are bounded by the respective shear layers (SLs).Also, laminarto-turbulent transitions at the main and flap LEs are detected, which are followed by the development of turbulent boundary layers (TBLs) along the respective suction sides.Finally, turbulent wakes are generated at the TEs of all elements.In the case of the slat and main wakes, these interact with the TBLs of the downstream elements.
Upon identifying the primary characteristics of the flow, Fig. 12 shows the evolution of the vortical structures with a.In general, augmenting a results in a reduction in the size of the slat recirculation bubble, while, at the same time, leading to the development of a TBL on the suction side and a more pronounced slat wake.In the main LE, the laminar-to-turbulent transition is identified at nearly constant locations, but the peak velocities are noticeably incremented with a, yielding higher suction peaks (see also Fig. 7).The vortical structures on the TBL prove to be larger as well, which is associated with the stronger adverse pressure gradient (APG) present at increasing a.Then, as in the slat, a more prominent main wake is developed.Regarding the flap, the turbulent transition is also localized at approximately constant locations on the LE.On the TE, flow separation is detected at a ¼ 5 and 9 along the suction side, while flow remains attached at a ¼ 23 .In this case, only the flow separation produced by the blunt TE is observed.More insights into some of the previously described flow phenomena are discussed in the following sections.In Sec.III D, the slat and main cove shear layers are addressed; the turbulent transitions are studied in more detail in Sec.III E; and finally, the wake topology is tackled in Sec.III F.

D. Shear layers
The formation of SLs in the slat and main coves can be devised in Fig. 13, which displays the instantaneous spanwise vorticity contours (x z c=U 1 ) in these zones.Their location can also be clearly visualized by the turbulent kinetic energy (TKE=U 2 1 ) contours illustrated in Fig. 14, since the maximum velocity fluctuations are found within the SL width.Similar flow dynamics are manifested in both, the slat and main coves, despite some differences encountered.The flow features addressed in this section can be visualized in Fig. 15 (Multimedia view), where a video of the instantaneous spanwise vorticity contours at a ¼ 5 is provided.
In the slat, SLs originate from the cusp in all cases, whereas in the main cove, their origin is slightly different depending on a.As shown from the streamlines in Fig. 16, flow along the main pressure side exhibits an incipient separation close to the TE at a ¼ 5 and 9 , thereby initiating the development of their respective SLs from these locations.As a is increased, the separation point is displaced toward the TE and, at a ¼ 23 , the flow remains attached, resulting in the formation of the SL further downstream when the flow reaches the cove cavity.These separations are also detected in the localized peaks devised in Fig. 17(c), which shows the root mean square of the fluctuating pressure coefficient (C 0 p;rms ).After the formation of the SLs, the abrupt velocity gradients give rise to natural Kelvin-Helmholtz (K-H) instabilities at a certain distance from the SL origin, where the vortex roll-ups are clearly identified in the x z c=U 1 contours of Fig. 13.Due to the inception of these instabilities, the two-dimensional structures are broken down into complex three-dimensional vortices and the transition to turbulence takes place.This transition process is evidenced in Fig. 14, wherein a localized increase in the TKE is observed.
Finally, the turbulent coherent structures originated along the SLs impinge into the lower surfaces of the slat and main coves, respectively.The location of the impingement can be traced by the C 0 p;rms distribution in Fig. 17.The impingement locally increases the pressure fluctuations and, hence, is identified as a local maximum, represented by circles in the plots.Note that the relative widths of the maxima can be associated with the left-to-right flapping motion of the impingement points, transmitted by the dynamics of the K-H vortices.Compared to the slat, the impingement in the main cove exhibits a wider motion, with the maximum fluctuations in the rightmost segment of the displacement, near the LE of the flap.In contrast, in the slat, the maximum values are observed at the midpoint of the displacement.In this region, similar C 0 p;rms distributions were reported in the literature 13,22,23,33 with comparable magnitudes in the impingement region, namely, in the range of 0.16-0.30.Although these numerical investigations were conducted at a higher Reynolds number, it all points out that the main differences stem from the turbulence model employed, as variability in this value is obtained among the literature depending on this choice.In those studies, a discrete peak near the TE was also observed, which is not shown in the present work.As demonstrated by Lockard et al., 22 this is associated with the blunt TE edge considered in those investigations in contrast to the sharp TE considered here.The blunt geometry promotes the formation of vortex shedding and the appearance of a high-frequency hump in the spectra.
The footprint of the impingement is also identified by the higher TKE values in Fig. 14.The increase in the TKE is solely visible in the slat, indicating a less prominent impingement in the main cove, as can also be observed in the C 0 p;rms peak magnitudes (Fig. 17).Indeed, the slat SL shows a greater prominence compared to the main SL, as substantiated by the higher levels of x z (Fig. 13), C 0 p;rms (Fig. 17), and TKE (Fig. 14) along the SL path.
As the angle of attack is progressively incremented, the TKE contours in Fig. 14 also reveal the mitigation of the turbulence levels in both the slat and main SLs.In the slat, it can be seen as a displacement of the stagnation point toward the slat cusp, reducing the recirculation bubble and creating a lower perturbation of the freestream in the cove, hence reducing the velocity fluctuations within this zone.Due to the shift of the stagnation point, a TBL is developed on the suction side of the slat at a ¼ 23 .The boundary layer transition can be localized by the sudden increase in the C 0 p;rms values (Fig. 17).In the main cove, variations of a do not significantly affect the size of the recirculation bubble, most likely due to the large surface of the main element, which constrains the development of the flow downstream.Instead, the observed differences stem from the separation of the SL along the pressure side, which advances the development of the SL upstream.As the separation is displaced toward the main cove cavity at higher a values, the velocity fluctuations within the cove are progressively reduced, while these are noticeably incremented on the main suction side.As in the slat, this is associated with the location of the stagnation point on the main LE, which moves downstream along the pressure side as a is increased.This increases the C p suction peaks and, hence, generates stronger streamwise APG, leading to a more pronounced development of the TBL on the suction side.
In Fig. 14, the location of some numerical probes, denoted as P1, P2, P3, and P4, is also depicted.These have been placed to capture the pressure signals and are positioned at specific locations along the SLs, namely, at the 20%, 40%, 60%, and 80% of the SL total length.The pressure spectra of these signals can be visualized in Fig. 18.As discerned from the graphs, probes located near the slat cusp, i.e., P1 and P2, capture a broadband peak.This peak is related to the K-H instabilities and is attenuated downstream, leading to a significantly broader spectrum.The transition to turbulence takes place, and hence, the energy is spread into wider ranges of frequencies, i.e., P3 and P4.As can be seen, the peak frequency is shifted toward lower frequencies and a broadband peak is still visible at P4, which is associated with the horizontal flapping of the impingement.Thus, the oscillatory nature of the K-H instabilities is transmitted downstream along the SL.
Strouhal number is defined as St ¼ fc=U 1 , with f being the frequency.Therefore, the frequency at which K-H instabilities occur diminishes with a since the driven shear force is reduced, i.e., the velocity gradient across the SL.

E. Laminar-to-turbulent transitions
A laminar-to-turbulent transition of the boundary layer is observed on the main and flap suction sides.These transitions are identified near the LEs of the respective elements (see Fig. 11).As discussed in Sec.III D, the slat boundary layer also exhibits a turbulent transition at a ¼ 23 .Nevertheless, the analysis of the slat is omitted in this section, and only the transition in the main and flap is assessed.
Not only transition to turbulence can be captured qualitatively by means of Q-criterion visualizations, but also the location of the transition can be detected in the C 0 p;rms distributions as shown in Fig. 17.However, to precisely identify the location of the transition, here, the non-dimensional skin friction distribution is used (see Fig. 19).This transition is manifested as a sudden reduction of this coefficient, and similar transition mechanisms are found in both elements.On the main LE, it appears that the mean C f vanishes in two locations, at about x=c ¼ 0:068 (T1) and x=c ¼ 0:170 (T2); whereas on the flap LE, this occurs at x=c ¼ 0:920 (T3) and x=c ¼ 0:945 (T4).
The first troughs on each element (T1 and T3), where the mean skin friction coefficient reaches almost C f ¼ 0, are associated with the natural transition processes.The emergence of the Tollmien-Schlichting (T-S) instabilities might be linked to the APG since the transition occurs immediately after the peak in the C p distributions (see Fig. 7).These peaks are located at nearly identical x/c locations (maximum curvature point), which leads to the development of the first disturbances of the flow at approximately constant locations for the different a considered.The frequencies of the T-S in the main LE can be clearly identified as peaks in Fig. 20, which shows the spectra of the velocity components (u, v, w) in T1.It is observed that as a is increased the frequencies of these peaks increase as well (St % 255, 285, and 340).Therefore, the stronger APGs at increasing a induce higher frequencies of the instabilities, and thereby, the more rapid growth of the fluctuations and further reductions of the local mean C f .
Once the T-S instabilities develop, the flow undergoes the transition process which entails the growth of the disturbances until reaching the fully turbulent state.This is identified by the Q-criterion isocontours in Fig. 21.The initial two-dimensional T-S waves are convected downstream and evolve into three-dimensional vortices, which experience breakdown in regions of high localized shear and eventually create turbulent spots to finally coalesce into fully turbulent flow.The growth of the disturbances in the transitional region is further assessed in Fig. 22, which shows the TKE contours together with the locus of the inflection point of the mean velocity profiles.This is defined as the height where the second derivative of the mean tangential velocity along the wall-normal direction changes sign.8][59] As substantiated by Fig. 22, the initial near-wall fluctuations appear along the locus of the inflection points, which denotes that the disturbances growth is dominated by the outer-layer inviscid instabilities.The viscous mode gains only importance in the main element transition at a ¼ 23 , where the initial disturbances are also observed in the near-wall region.A similar behavior is reported by Tamaki and Kawai, 60 wherein they associate the prevalence of the inviscid mode with the relatively low freestream Reynolds number and  observed that, at higher Reynolds numbers, the transition process shifts toward the viscous mode.In the present LES, the local Reynolds number would be incremented by the higher velocities observed at higher angles of attack.
After the natural transition, a second trough (T2 and T4) can be identified in the main and flap C f , which entails a lower reduction of the coefficient (see Fig. 19).In this case, this is caused by the interaction with the slat and main wakes, respectively.As illustrated in Fig. 23, some turbulent coherent structures shed from the slat and main TEs eventually reach the surface of the main and flap LEs, respectively.These structures swept the boundary layer, yielding a localized reduction of the C f .In the flap, the C f reduction due to this phenomenon is much lower compared to the effect in the main.This is believed to be linked with two facts: the divergent geometry of the flap, which tends to prevent the direct impingement of the coherent structures onto the flap surface, and the less prominent activity of the main cove compared to the slat cove.Wakes are partially formed by the turbulent coherent structures shed from the coves.However, in the main cove, most of them are dissipated before reaching the TE and, thus, are not convected toward the flap.All this translates into fewer impingement episodes, and hence, a lower affection of the C f .This is in alignment with the findings of Terracol and Manoha 39 in the main LE of the F16 high-lift wing at Re c ¼ 1:23 Â 10 6 , who also reported the existence of these two flow phenomena, i.e., the natural transition and the slat wake interaction.In their case, the localized turbulent spots originated by the slat wake preceded the natural  transition, inverting the order of the troughs in the C f distribution, which might be simply attributed to the geometrical differences in the wing configuration.
Skin friction predictions are then compared with the theoretical correlations of a flat plate. 61Figure 24 shows the C f nondimensionalized with the local edge velocity of the boundary layer (U e ), i.e., C f ¼ s w =0:5qU 2 e .Furthermore, the comparison is also made with the boundary layer thickness (d) in Fig. 25.In all these plots, the distance s is used in the x axis to allow the comparison with a flat plate, being this the distance from the stagnation point following the curvature of the surface.As can be observed from the graphs, C f and d values within the laminar regions deviate from the theoretical solutions.This might be attributed to historical effects and the high curvature of the main and flap LEs.On the other hand, further downstream and after  the turbulent transition, computations resemble pretty well the results of the TBL in a flat plate.In the main element, the major discrepancies appear for the case at a ¼ 23 ; whereas in the flap, discrepancies are found for the predictions at a ¼ 5 and 9 .These differences can be attributed to the stronger APGs.Moreover, near the TEs more differences arise.In this case, these are related to other flow phenomena: the layer lift-up in the main and the TE flow separation in the flap.Comparable C f distributions were predicted by Terracol and Manoha 39 along the suction side of the DLR-F16 main element at a ¼ 6:15 .Although the quantitative differences in the fully turbulent region can also be linked to the different geometry or angle of attack, the slight offset toward lower values is consistent with the explanations given in Sec.III B regarding the higher Reynolds number effects on the C f .

F. Wake turbulence
As shown in Fig. 11, wake turbulence is exhibited downstream of each of the airfoil elements.In Fig. 26, the slat wake can be visualized through the representation of velocity magnitude (U=U 1 ) and turbulent kinetic energy (TKE =U 2 1 ) contours.This originates from the confluence of the flow on the slat suction and pressure sides.On the suction side, the boundary layer thickness grows with a, and larger turbulent structures can be detected, as discussed in Sec.III D. On the pressure side, following the impingement of the SL, some turbulent structures are re-entrained into the cove, whereas others are convected downstream and noticeably accelerated through the gap between the slat and the main element, resulting in an elongation of these structures; all in an intermittent fashion due to the flapping motion of the shear layer.From this figure, it can be seen that the width of the wake increases with a, especially due to the turbulent structures emanating from the thicker TBL.
The flow over the flap results from the confluence of the slat wake and main wakes, the main-flap gap jet, and the flap boundary layer, as described in Fig. 27, which depicts the U=U 1 and TKE/ U 2 1 contours in this region.It is then detected that, as the a is increased, the APG on the suction side of the main element becomes stronger, leading to the development of a more prominent TBL and, consequently, more discernible wake structures downstream.The main wake is apparent at a ¼ 5 and 9 but becomes particularly pronounced at a ¼ 23 , suggesting that the wing is in stall conditions.Similar flow features were described by Bodart et al. 21   Regarding the flap wake, this remains relatively modest in scale compared to the main wake and it is reduced as a is incremented.At a ¼ 5 and 9 , the wake is initiated upon the separation of the flow along the flap suction side.This separation is also described by Jin et al. 28 at a ¼ 5:5 and Re c ¼ 1:7 Â 10 6 at a qualitatively similar location.Conversely, at a ¼ 23 , no separation is discerned due to the noticeable spread of the main wake, which maintains the flow attached along the surface of the flap.Note that in all cases, the initial spread of the main wake is limited by the entrance of the jet, which dissipates downstream and allows the development of the wake.At a ¼ 23 , the growth of the main wake pushes the jet toward the flap surface.Irrespective of the angle of attack, K-H instabilities are identified immediately downstream of the flap TE due to the blunt geometry, which creates a localized small recirculation area here and hence triggers these instabilities.

IV. CONCLUSIONS
Large eddy simulations (LES) of the flow past a 30P30N wing are conducted at a Reynolds number of Re c ¼ 750 000 and three different angles of attack a ¼ 5 , 9 , and 23 .A mesh grid independence analysis is first carried out and two meshes are compared: the baseline mesh and a finer one.Overall, small differences are found between the two meshes, which showcases the independence of the results with respect to the selected resolution.For the sake of lowering the computational time, the baseline mesh is employed throughout the whole work.Results are then validated against experimental and numerical investigations available in the literature, proving the reliability of the present predictions.The validated computational framework is then employed to examine the flow dynamics and assess how they evolve as a is progressively increased.The different features and vortical structures identified throughout the flow field are studied.
Large recirculation bubbles in the slat and main coves are observed.These structures are bounded by shear layers and exhibit a complex three-dimensional nature.Initially, Kelvin-Helmholtz two-dimensional instabilities are formed within the shear layers, which are then broken into smaller turbulent structures.The impingement of these structures into the coves inner surface increases pressure fluctuations and turbulent kinetic energy levels noticeably, with the slat cove showing a more prominent impingement compared to the main cove.The strength of the impingement is considerably reduced at higher angles of attack.In the slat cove, the shift of the stagnation point toward the cusp as the angle of attack is increased leads to a reduction of the recirculation bubble size and diminishes the turbulent activity along the shear layer.This, in turn, results in a less prominent impingement.Opposed to what is observed in the slat cove, where the shear layer originates from the cusp, in the main cove, the shear layer originates from the separation point of the flow along the pressure side.As the angle of attack is incremented, the separation is displaced toward the cove cavity, resulting in a reduction of the turbulent kinetic energy levels inside the cove and a corresponding decrease in the impingement strength.In the main cove, the size of the recirculation bubble is maintained practically constant.
Laminar-to-turbulent transitions are detected at the main and flap suction sides, near the leading edges.These transitions are induced by the Tollmien-Schlichting instabilities that arise due to streamwise adverse pressure gradients in the zone.In this manner, the transition occurs near the pressure suction peak, where the higher suction levels at increasing angles of attack provoke more local separation events (lower mean C f ) and increase the frequencies of the Tollmien-Schlichting instabilities.In all cases, the growth of the disturbances is dominated by the inviscid mode.Only at the main leading edge at a ¼ 23 , the appearance of the viscous mode starts to be visible, most likely due to the higher velocities encountered in the region, which increases the local Reynolds number.
Wake turbulence downstream of all elements plays a significant role in the flow characteristics.The slat and main wakes interact with the turbulent boundary layers of downstream elements, affecting the flow patterns.As the angle of attack is increased, more prominent wakes are detected.At a ¼ 23 , the main wake becomes notably large, exhibiting von K arm an-like vortices and indicating the entrance to stall conditions, thus limiting the maximum lift coefficient.The frequencies of these oscillations are also detected in the lift and drag coefficient signals.The flap wake, on the other hand, is reduced as a increases.At lower angles, the flap wake originates from the separation of flow along the suction side, while at higher angles, the absence of separation is attributed to the spread of the main wake.Kelvin-Helmholtz instabilities are formed downstream of the flap trailing edge due to the blunt geometry.ing (equal).Ivette Rodriguez: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal);

FIG. 3 .
FIG. 3. Computational mesh near the airfoil walls, together with the global and local coordinate systems.

FIG. 15 .
FIG.15.Video of the instantaneous spanwise vorticity contours (x z c=U 1 ) in the slat and main coves at a ¼ 5 .Multimedia available online.

FIG. 20 .
FIG. 20.Energy spectrum of the (a) u, (b) v, and (c) w velocity component signals (E uu ; E vv , and E ww , respectively) at T1.

FIG. 19 .
FIG. 19.Zoomed view of the skin friction coefficient (C f ) along the (a) main and (b) flap suction sides.

FIG. 21 .
FIG. 21.Main leading edge instantaneous vortical structures identified by means of Q-criterion isocontours colored by the velocity magnitude at a ¼ 5 .
FIG. 29.Instantaneous velocity magnitude contours (U=U 1 ) in the flap region at a ¼ 23 along a vortex shedding period.

TABLE I .
Experiments on the flow past the 30P30N wing model at moderate to high Reynolds numbers.

TABLE II .
Computations on the flow past the 30P30N wing model at moderate to high Reynolds numbers.
FIG.1.Schematic of the 30P30N wing in stowed (above) and deployed (below) positions.FIG.2.Schematics of the computational domain and boundary conditions.

TABLE III .
Fine and baseline mesh properties.

TABLE IV .
Predicted aerodynamic coefficients by the fine and baseline grids at a ¼ 5 .

TABLE V .
Points ðx; yÞ=c defining the location of the L2, L4, and L7 lines.