On the Flow Past a Three‑Element Wing: Mean Flow and Turbulent Statistics

Large eddy simulations (LES) on the flow past the 30P30N three-element high-lift wing at a moderate Reynolds number Re c = 750, 000 and three different angles of attack 𝛼 = 5 , 9 and 23 ◦ are conducted. The main focus is on the time-averaged statistics of the turbulent flow. The form drag noticeably increases with the angle of attack, while viscous drag remains roughly constant and contributes minimally to the total drag. This is associated with the significant pressure peaks found in the main element with increasing angles of attack and hence, the development of stronger adverse pressure gradients. At 𝛼 = 23 ◦ , this leads to the development of a prominent wake downstream this element that eventually evolves into a visible recirculation region above the flap, indicating the onset of stall conditions. In the flap, strong adverse pressure gradients are observed at small angles of attack instead, i.e., 𝛼 = 5 and 9 ◦ . This is attributed to the flap’s deflection angle with respect to the main wing, which causes a small separation of the boundary layer as the flow approaches the trailing edge. At the stall angle of attack, i.e., 𝛼 = 23 ◦ , the spread of the main element wake maintains attached the flow near the flap wall, thus mitigating the pressure gradient there and preventing the flow to undergo separation. The shear layers developed on the slat and main coves are also analysed, with the slat shear layer showing more prominence. In the slat, its size and intensity noticeably decrease with the angle of attack as the stagnation point moves towards the slat cusp. Conversely, the size of the shear layer developed in the main element cavity remains approximately constant regardless of the angle of attack. At the lower angles of attack, i.e., 𝛼 = 5 and 9 ◦ , the development of the shear layer is anticipated by the turbulent separation of the flow along the pressure side of the main wing, leading to increased levels of turbulence downstream. At the higher angle of attack, i.e., 𝛼 = 23 ◦ , the shear layer is originated by the cavity separation and transition to turbulence occurs within the cavity.


Introduction
In civil aviation, the use of high-lift devices is indispensable during takeoff and landing operations, when the aircraft speed is considerably lower than during cruise.The deployment of these devices gives rise to elevated noise levels, representing a key design factor to reduce the emitted airframe noise.Moreover, the optimization of these elements can also entail other benefits such as lower drag coefficients, enhanced payload capacities, shorter airfield distances, or even smaller operational empty weights (OEW) (van Dam 2002).Consequently, understanding the flow complexities taking place in multi-element wings is crucial for designing a more efficient aviation, with computational fluid dynamics (CFD) playing an important role in achieving so.
In this regard, various workshops can be found in the literature that aimed to establish CFD capabilities in dealing with multi-element wings.This is the case of the high-lift CFD challenge workshop held by NASA (Klausmeyer and Lin 1997) or the third workshop on benchmark problems for air-frame noise computations (BANC-III) organized by AIAA (Choudhari and Lockard 2015), both focusing on the benchmark configuration referred as 30P30N.
In the former, all computations were based on Reynolds-averaged Navier-Stokes (RANS) models or potential/Euler solvers and compared the predicted aerodynamic coefficients against experimental results to assess the accuracy of the computations.In general, all investigations showed good results, despite the analysis being limited to general aerodynamics parameters and mean statistics.
The interest of the latter was placed on accurately computing the noise emitted by the slat cove, identified as the main source of noise in multi-element wings.Therefore, the dynamics of the flow had to be analysed by means of unsteady CFD approaches.Both experimentalists (Jenkins et al. 2004;Pascioni et al. 2014;Pascioni and Cattafesta 2016;Murayama et al. 2018) and computational investigators (Bodart et al. 2013;Lockard et al. 2014;Terracol et al. 2015;Ashton et al. 2016;Housman et al. 2019) contributed to provide detailed flow measurements to the workshop.In those studies, it was pointed out that the cavity-like geometry of the slat cove excites acoustic modes similar to those in open cavity flows and hence, generates a series of narrowband peaks in the low to middle frequency range due to a feedback mechanism.They also identified a broadband noise at a higher frequency associated with the turbulent fluctuations in the slat cove and the vortex shedding occurring downstream in the trailing edge.Most of the computational contributions at BANC-III employed hybrid RANS-LES approaches, such as detached eddy simulations (DES), delayed detached eddy simulations (DDES) or other variations.Only one wall-modelled large eddy simulation (WMLES) was presented (Bodart et al. 2013).From here, several subsequent works have used the results of BANC-III to evaluate the accuracy of their novel methodologies and thus, did not provide any additional insights into the physics of the flow.Among these investigations, higher order (Shi et al. 2018;Gao et al. 2020;Jin et al. 2020) and reduced dissipation schemes (Sakai et al. 2019), adaptive mesh refinements algorithms (Ueno and Ochi 2019), and variations of the classical hybrid RANS-LES turbulence models (Kojima et al. 2020;Shur et al. 2023) were tested.Additionally, another WMLES was published (Zhang et al. 2017), which represented the second work applying this advanced turbulence modelling approach.In the experimental side, Li et al. (2017) validated a novel concept for the test sections of wind tunnels using tensioned Kevlar cloth and a perforated plate.This was aimed at accurately obtaining aerodynamic and aeroacoustic measurements with minimal interference of the sidewalls.Later, Li et al. (2018) extended their measurements on this wing model to put more emphasis on the physical analysis, effectively capturing the mid-frequency peaks reported in the previous studies.
Besides the 30P30N configuration, other high-lift wing geometries have been investigated in the literature.In this context, ONERA (the French Aerospace Lab) and DRL (the German Aerospace Center) collaborated on a joint project called LEISA, aiming to develop silent take-off and landing devices.Within this project, Deck and Laraufie (2013) performed zonal detached eddy simulations (ZDES) on the DLR-F15 three-element high-lift wing model, providing a thorough description of the flow around the slat and flap regions and extending the aeroacoustic analysis to the main cove cavity as well.Nevertheless, the lack of validation data motivated the follow-on LEISA2 (Manoha and Pott-Pollenske 2015) project, in which the two institutions built an extensive experimental database to validate the numerical and aeroacoustics computations of high-lift wings.The DLR-F16 configuration was selected as the reference wing model for this project.Later, Terracol and Manoha (2020) performed on this wing configuration the only wall-resolved large eddy simulations (LES) documented in the literature for this type of geometries, as far as the author is concerned.In their work, they provided a very detailed analysis of the flow around the three elements along with aeroacoustic measurements.The DLR-F16 model was also used to explore the capabilities of the Coanda effect to delay or remove the separation over the flap at high deflection angles, and hence enhance the lift coefficient of high-lift wings.This work was performed employing zonal overset large eddy simulations (Ananthan et al. 2020).
All the studies presented until now are performed at moderate Reynolds numbers ( Re c ∼ 10 6 ).Studies at a lower Reynolds number regime ( Re c ∼ 10 3 − 10 4 ) can also be found in the literature in both the experimental (Wang et al. 2018(Wang et al. , 2019) ) and computational (Vadsola et al. 2021) fields.These studies showed a distinct flow behaviour and observed the presence of Görtler vortices in the slat wake, which appears below a critical Reynolds number when no roll-ups of the slat shear layer are observed.In the moderate Reynolds number regime, the slat wake is dominated by streamwise and spanwise vortices instead.
Although the numerous studies performed in the past regarding the flow over high-lift wings, there is still a gap in the literature concerning scale-resolving simulations on this type of flow.This is particularly relevant considering the potential of LES to analyse complex and unsteady flows, where flow separations, shear layers and wakes might develop.Moreover, the existing literature on the 30P30N wing configuration is rather limited, primarily focusing on the dynamics of the slat at a fixed angle of attack.The present work aims at studying the effects brought by the angle of attack on the turbulent flow statistics, thereby extending the analysis to include previously unexplored parameters, such as the boundary layer development in high-lift wings.To this end, LES of the flow over the 30P30N wing at a Reynolds number of Re c = 750, 000 and three different angles of attack, namely = 5 , 9 and 23 • , are performed.Following up on a preceding publication (Montalà et al. 2024) that was primarily centred on the examination of the unsteady flow characteristics, the focus herein is shifted towards the mean features of the flow.The paper has been structured as follows: Sec. 2 introduces the mathematical framework used for the simulation (Sec.2.1), as well as the computational domain employed (Sec.2.2).In Sec. 3 the results obtained for the different angles of attack are discussed, first analysing the aerodynamic coefficients (Sec.3.1) and then studying different regions of the flow in detail, i.e., the boundary layer development along the main and flap elements (Sec.3.2), the shear layers formed in the slat and main cavities (Sec.3.3) and the wakes originated downstream each element (Sec.3.4).Finally, conclusions are presented in Sec. 4.

Numerical Configuration
As mentioned in the introduction, large eddy simulations (LES) of the flow past the 30P30N high-lift three-element wing configuration are performed at a Reynolds number Re c = 750, 000 and the angles of attack = 5 , 9, and 23 • .Note that the Reynolds number is based on the nested chord c of the wing and is defined as Re c = U ∞ c∕ ; with being the fluid kinematic viscosity and U ∞ the freestream velocity.

Mathematical and Numerical Modelling
LES are conducted by solving the incompressible spatially filtered Navier-Stokes equations.These can be written as where u i (or u, v and w) and p are the filtered (or resolved) velocity and pressure fields, the fluid density, and ij the subgrid-scale (SGS) stress tensor.The deviatoric of this term is modelled via an eddy viscosity ( sgs ) model, and hence the resulting stresses can be written as 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 ij is the Kronecker delta.In this case, the eddy viscosity sgs is evaluated through the Vreman (2004) model.
The set of equations described above are solved with Alya (Vázquez et al. 2016), an inhouse finite-element (FE) code developed at the Barcelona Supercomputing Center (BSC).In this code, the convective operator of the equations is approximated by a low-dissipation scheme that preserves energy, momentum and angular momentum at the discrete level (Lehmkuhl et al. 2019).Numerical dissipation is only introduced through a fractional-step method (Codina 2001) to stabilize the pressure, allowing the use of equal order interpolations for the velocity and pressure.Regarding the time integration, this is advanced using an energy-conserving 2nd-order Runge-Kutta explicit method (Capuano et al. 2017) combined with an eigenvalue-based time-step estimator (Trias and Lehmkuhl 2011), which allows maximizing the time step by adapting the linear stability domain of the dynamical system, always ensuring that the CFL condition is fulfilled.

Computational Domain
The geometrical configuration of the case is depicted in Fig. 1.This also describes the global (xyz) and local ( z ) coordinates system, the latter representing a rotation of the global coordinate system around the z direction to align the axis towards the wall (1) tangential and normal directions, with the origin on the wing surface.The origin of the global coordinate system is located at the leading edge of the slat in the stowed configuration.In the deployed configuration studied here, the origin is maintained at the same location as in the stowed one and, as illustrated in Fig. 1, this is roughly at the trailing edge of the slat.
Solutions are obtained in a circular domain (x-y plane) of radius r xy = 10c that is extruded spanwise (z-direction) a distance L z = 0.1c .As shown in Fig. 2, this length ensures that the two-point correlations of the x-and y-velocities along the z direction drop to nearly zero before reaching the mid-span.In other words, this length is large enough to contain the largest turbulent scales of the flow.This is in agreement with the suggestions done at the BANC-III workshop (Choudhari and Lockard 2015).
For the boundary conditions, a uniform velocity profile is applied at the inlet, which depends on the angle of attack as (u, v, w) = U ∞ (cos , sin , 0) .At the outlet, a pres- sure-based condition is imposed, i.e., p − 0.5 u 2 n = 0 , with u n being the normal veloc- ity at the boundary and an expansion coefficient set to = 1 in this work.The no-slip condition is prescribed at the wing walls, i.e., (u, v, w) = (0, 0, 0) , while periodic bound- ary conditions are set in the spanwise direction.
Concerning the grid, the present computations are obtained employing a hybrid unstructured mesh of about 58 million grid points.Structured-like inflation layers around the wing walls are considered, which allows adapting the near-wall spacings according to the LES requirements, namely Δ + < 2 , Δ + ≈ 50 − 150 and Δz + ≈ 15 − 40 (Piomelli and Chas- nov 1996), without notably increasing the number of grid points.In the vicinity of the walls, maximum non-dimensional distances of Δ + < 1 and Δ + < 80 are forced in the wall-normal and streamwise directions, respectively.This results in a two-dimensional mesh of 449, 000 grid points.This is then extruded in the z-direction using 129 equal planes, which leads to a near-wall spanwise spacing of Δz + < 50 .Note that the superscript + denotes that quantities have been non-dimensionalized with the scaling factor u ∕ , where u = √ w ∕ represents the friction velocity and is computed with the wall shear stress w .A second level of mesh refinement was further considered, as shown in Table 1.In Montalà et al. (2024), a thorough grid independence analysis was conducted, proving that the baseline mesh was a reasonable compromise between accuracy and computational effort.Therefore, the results shown in the present work are obtained with this mesh.

Results and Discussion
To initialise the simulations, the velocity field resulting from running 20 U ∞ ∕c in a coarser mesh is interpolated into the mesh described in Sec.2.2.This coarser simulation is initialised with a homogeneous velocity field and is run until reaching the statistically steady state.Then, after initializing the finer mesh, flow is advanced in time for about 3 U ∞ ∕c before statistics are collected during approximately 14 U ∞ ∕c .To obtain the statistics of the flow, the instantaneous field is averaged both in time and in the spanwise direction.Note that all quantities presented throughout the whole work are non-dimensionalized using the reference length, velocity, and time as c, U ∞ , and c∕U ∞ .
The instantaneous vortical structures identified through the Q-criterion iso-contours are visualized in Fig. 3.This figure highlights the unsteady nature of the flow.Cavity flow types are observed in the slat and main coves, characterized by large recirculation bubbles bounded by the respective shear layers.Laminar-to-turbulent transitions occur at the  leading edge of each element.For the main and flap elements, this transition is observed at all the angles of attack, whereas for the slat, the transition is only triggered at the highest angle.Downstream of these transitions, the development of wall-bounded turbulence is observed.Wake turbulence can also be identified downstream of each element.The wakes of the slat and main elements become apparently more prominent as the angle of attack increases.Concerning the flap, flow separation can be glimpsed at the smaller angles of attack, leading to wake formation behind this element.Nevertheless, this is not observed at the highest angle of attack.Instead, the flow downstream of the wing is dominated by the main element wake.The analysis of all these transient phenomena was conducted in Montalà et al. (2024).In the present work, the interest is placed on the time-averaged statistics of the turbulent flow.
A general overview of the predicted mean flow field is visualized in Fig. 4.This shows how the streamlines evolve with the angle of attack.The following flow features can be  concentration of the streamlines along thin layers.These represent the location of the mixing layers formed by the interaction between the flow from the upper and lower surfaces of the wing when they converge at the trailing edge of the element.The mixing layers are then convected downstream along the wake.It is detected that the wakes of the slat and main element separate from the wing surface as the angle of attack increases.This expansion is especially noticeable at = 23 • , where a spread of the streamlines is identified, yielding a low-velocity recirculation area above the flap.All these effects are addressed in more detail in the sections below.

Aerodynamic Coefficients
The pressure C p = (p − p ∞ )∕(0.5 U 2 ∞ ) and skin friction C f = w ∕(0.5 U 2 ∞ ) coefficients along the wing walls for the three angles of attack studied are depicted in Fig. 5.These quantities are validated against experimental studies available in the literature.Concretely, the investigations carried out in the FSAT (Pascioni et al. 2014), JAXA (Murayama et al. 2018) and LTPT (Klausmeyer and Lin 1994) facilities.For a more thorough validation, including a grid independence analysis, the reader is referred to Montalà et al. (2024).
The pressure coefficient (see Fig. 5a) shows a good agreement with the experimental data.At increasing angles of attack, higher suction peaks are developed on the slat and main elements.As discussed in Fig. 4, this is related to the acceleration of the flow along the leading edge of these elements as the angle of attack increases and the stagnation point is displaced downstream.In the flap, the pressure coefficient shows slight variations with the angle of attack as the flow around this element is highly constrained by the large dimensions of the main element.A small reduction of the C p peak is detected.This decrease might be understood by the incidence angle of the freestream flow and the loss of momentum linked to this.From the streamlines in Fig. 4, one may notice that the flow on the flap leading edge is mainly transported from the main wing pressure side through the gap between both elements.As the angle of attack is increased, the main element pressure side exhibits lower velocities, which in turn results in smaller velocities in the flap leading edge and hence, smaller suction peaks.
As for the skin friction (see Fig. 5b), higher discrepancies are found with the experiments.This might be associated with the differences in the Reynolds number between the present computations and the experiments since the C f is more sensitive to variations in this parameter.Nonetheless, both the computations and the experiments show the same trend and relatively small differences are detected as the angle of attack is modified.In the present computations, the evolution of C f with the angle of attack can be explained by differences in the transitional regions near the leading edges and by the adverse pressure gradients developed downstream in the wall-bounded turbulence region, i.e., as the flow approaches the trailing edge.
The integration of the C p and C f distributions (local coefficients) over the wing sur- face yields the aerodynamic force acting on it.By decomposing the force in the direction parallel and perpendicular to the freestream, the lift C L = L∕(0.5SU 2 ∞ ) and drag C D = D∕(0.5SU 2 ∞ ) coefficients can be computed, respectively; with S = c L z denoting the reference surface.Therefore, these quantites represent the integral values of the C p and C f distributions.The obtained lift and drag coefficients for the three angles of attack are illustrated in Fig. 6.The partial contribution of each element to the total C L and C D coefficients is also included.
As observed in Fig. 6a, the primary contributor to the lift coefficient is the main element, characterized by its larger surface and higher suction peaks in the pressure coefficient.These are noticeably increased with the angle of attack (see Fig. 5a), resulting in an increased C L coefficient for this element.In the flap, the C p peaks are kept relatively modest and slightly decrease with the angle of attack, causing the lift coefficient to follow the same trend.Regarding the slat, higher lift coefficients are observed, which are linked to the suction peaks developed as the angle of attack increases.Nevertheless, the small surface of this element keeps its contribution low.
The drag coefficient is depicted in Fig. 6b.The contribution from the flap remains nearly constant for all the angles of attack considered, with a slight increase as the angle increases.Generally, this element shows the higher drag values.Larger variations are found in the contributions from the main and slat elements.Additionally, they both exhibit opposite behaviours.As the C D of the main element increases with the angle of attack, that of the slat decreases.Nonetheless, the net effect is an increase in the total drag coefficient.
It is worth pointing out that, conversely to the lift coefficient, the absolute C D contri- bution of both the slat and main elements is nearly equivalent despite the difference in the wet surface.This is related to the inclination of the slat, which is s = 30 • deflected Fig. 6 Aerodynamic coefficients at the different angles of attack together with the partial contribution from each wing element.a Lift coefficient and b drag coefficient with respect to the main element.This results in a significant portion of the pressure forces being directed upstream, i.e., in the opposite direction of the drag force.
While pressure forces are primarily responsible for the lift coefficient, viscous forces might have a greater impact on the drag coefficient.These act in the streamwise direction of the flow and predominantly contribute to increasing the drag.In Fig. 7, the form (or pressure) C D,press. and viscous C D,visc.drag coefficients generated by each element are computed.
As a result of the small variations of C f with the angle of attack (see Fig. 5b), it is observed that C D,visc. is approximately constant for the different angles of attack considered.As expected, the largest contributor is the main element, as it involves higher values of C f integrated over a larger surface, whereas the contribution of the remaining elements is relatively small.
On the other hand, the C D,press. shows higher variations with the angle of attack, which are linked to the visible changes in the C p distribution (see Fig. 5a).Furthermore, the pres- sure contribution to the total drag is considerably higher than the viscous one, about one order of magnitude higher.This explains why the evolution of C D illustrated in Fig. 6b is closely related to the C D,press. shown in Fig. 7a.

Boundary Layer Development
As evidenced by the sudden increase in the C f shown in Fig. 5b, a laminar-to-turbulent transition takes place at the leading edge of the main and flap elements ( x∕c ≈ 0.10 and x∕c ≈ 0.93 , respectively).After this, turbulent coherent structures shed from the slat and main wakes interact with the main and flap boundary layers, respectively, leading to a reduction in the C f .This phenomenon was extensively addressed in Montalà et al. (2024).Following this transitional region, a turbulent boundary layer (TBL) is developed along the main and flap suction sides, i.e., x∕c > 0.2 and x∕c > 0.95 , respectively.In the slat suction side, a TBL is also observed at = 23 • (see Fig. 3c).However, its analysis is omitted in the present work.
In Fig. 8, the TBL along the main element suction side is analysed.Figure 8a shows the tangential velocity profiles at different chord locations; together with the boundary layer thickness , computed as the height where the tangential flow velocity reaches the 95% of the local edge velocity U e .In Fig. 8b, the evolution of the momentum thickness Reyn- olds number Re is depicted.This is calculated in terms of the edge velocity U e and the momentum thickness .Finally, the shape factor H and Clauser pressure gradient parameter along the main suction side are presented in Fig. 8c and Fig. 8d, respectively.The shape factor is defined as the ratio between the displacement and momentum thicknesses H = * ∕ .These two thicknesses are obtained through the integration of the quantities 1 − u ∕U e and u ∕U e (1 − u ∕U e ) , respectively, up to the boundary layer edge .Regard- ing the Clauser pressure gradient parameter, this is computed as = * ∕ w dP e ∕dx , where dP e ∕dx represents the tangential pressure gradient at the boundary layer height.
From Fig. 8, it can be detected that predictions at the lower angles of attack, i.e., = 5 and 9 • , show a modest boundary layer growth compared to the case at = 23 • .This is evi- denced in the evolution of and Re depicted in Fig. 8a and Fig. 8b, respectively.As mov- ing downstream, at = 5 and 9 • , the growth rate of these two quantities is progressively reduced and, as was shown in the preceding investigations (Montalà et al. 2024), tends to the typical evolution observed in a flat plate with a zero pressure gradient (ZPG).Conversely, at = 23 • , the growth of the boundary layer gets more pronounced near the trailing edge.These effects can be understood by analysing H and in Fig. 8c and Fig. 8d, respectively, which are indicators of the pressure gradient.At = 5 and 9 • , the Clauser pressure gradient parameter tends to zero as the trailing edge is approached (see Fig. 8d), while at the same time, the shape factor gets close to the range of H = 1.3 − 1.4 (see Fig. 8c), the expected values in a ZPG flat plate.Despite some theoretical correlations pointing out that this is H = 1.3 (White 2006), recent computations have shown that this value would be closer to H = 1.4 (Schlatter and Örlü 2010; Kumar and Dey 2019).Therefore, the TBL resembles the behaviour of a ZPG, which can also be observed in the pressure coefficient plateau in the rear part of the main element (see Fig. 5a).On the other hand, at = 23 • , the pressure coefficient keeps decreasing as the trailing edge is approached and both H and increase considerably instead, as illustrated in Fig 8c and Fig. 8d, respectively.Hence, due to the augmented angle of attack, flow experiences a stronger adverse pressure gradient (APG), resulting in a significant flow deceleration in the streamwise direction and an increase of the wall-normal convection.Nevertheless, despite the APG, no separation is detected near the trailing edge and, as shown in Fig. 5b, the skin friction coefficient does not reach C f = 0 at any chordwise location.Finally, the abrupt increase in (Fig. 8a), Re (Fig. 8b) and H (Fig. 8c) observed in the final section of all plots, as well as the sudden reduction in (Fig. 8d), corresponds to the small lift-up of the boundary layer induced by the flow passing through the gap between the main and flap elements.
In Fig. 9, the TBL along the flap suction side is presented in the same fashion as for the main element.Opposite to what is observed in the main element, the flap exhibits a distinct behaviour, with the lower angles of attack displaying more pronounced TBL growth rates (see and Re in Fig. 9a and Fig. 9b, respectively).Despite being at a low angle of attack, i.e., = 5 and 9 • , the divergent-like geometry of the flap, this being located at a deflection angle of f = 30 • with respect to the main element, increases the APG and so as the Clauser pressure gradient parameter and the shape factor, as observed in Fig. 9c and Fig. 9d, respectively.The intensified APG, with and H growing exponentially, induces flow separation near the flap trailing edge.This phenomenon can be identified through both, the streamlines depicted in Fig. 4a, b and the skin friction coefficient in Fig. 5b, which reaches C f = 0 at approximately x∕c = 1.12 .For the predictions at = 23 • , despite the C p initially following the same trend near the leading edge as for the other angles of attack, it stops decreasing at a given chord position (see Fig. 5a), resulting in a diminished APG.Consequently, for this angle of attack, it can be detected that H (Fig. 9c) and (Fig. 9d) exhibit lower values and remain nearly constant at the different chordwise positions, leading to no flow separation downstream.Note that the flow over the flap is highly affected by both the main TBL and the jet formed by the flow passing through the main-flap gap.At = 23 • , the flow from the main element TBL exhibits a higher vertical momentum and diminished streamwise velocities.As visualized in Fig. 4c, this induces a low-pressure zone over the flap which leads to the separation of the mixing layer from the main element wake.Precisely, this expansion of the streamlines maintains the jet flow coming from the gap between the main and the flap elements attached to the flap walls, preventing the flow from undergoing a strong APG.Nevertheless, despite the more favourable flow conditions developed in the flap TBL, a recirculation zone is created away from the wall, associated with an increase in the total wing drag.
In Fig. 10, the Reynolds stresses in wall units (+) are displayed at different chord positions of the main ( x∕c = 0.4 , 0.6, 0.8) and flap ( x∕c = 0.975 , 1.025, 1.075) suction sides, as highlighted by the dashed lines in Fig. 8 and Fig. 9, respectively.On the main wing surface (see Fig. 10a, b and c), higher magnitudes of tangential Reynolds stresses are observed compared to the other stresses, characterized by the discernible inner peak typically found in wall-bounded flows ( y + < 100 ).The magnitude of this peak shows to be approximately the same for all the angles of attack and located at the same y + height, highlighting the uni- versality of the flow in the near-wall region.The magnitudes of the other stresses, i.e., the wall-normal and the shear stresses, are nearly zero in this region and get higher as moving further away from the wall ( 100 < y + < 1000 ) due to the damping of the vertical velocity fluctuations near the wall.Additionally, this outer region is highly affected by the outer conditions of the flow, i.e., the APG.As moving downstream along the different stages, the impact of the APG present at = 23 • can be observed.At this angle of attack, the wall- normal fluctuations become more significant in the outer region of the TBL as approaching the trailing edge, accompanied by the development of a prominent outer peak in the streamwise stresses and an increase of the shear stresses as well.This was already observed in a flat plate (Lee and Sung 2008) and in a NACA0012 airfoil (Vinuesa et al. 2017).The APG contributes to reducing the streamwise momentum and enhances the wall-normal convection of the flow, thereby intensifying the turbulent activity in the TBL outer region and leading to more energetic large-scale motions in this region (Vinuesa et al. 2018).For the lower angles of attack, i.e., = 5 and 9 • , where nearly a ZPG TBL is obtained, the Reynolds stresses exhibit little variation along the different chord positions.
On the flap surface (see Fig. 10d, e and f), the Reynolds stresses also display similar structures as in the main wing: the inner peaks in the streamwise stresses are identified, while the wall-normal and shear stresses contain most of their energy in the outer region of the boundary layer.In this case but, as shown in Fig. 9d, a stronger APG is detected at the angles of attack = 5 and 9 • , and not at = 23 • .Hence, the development of an outer peak in the streamwise Reynolds stresses, along with a general increase of the wall-normal and shear stresses in this region, is observed at these angles of attack.Indeed, the elevated values of observed in the flap (see Fig. 9d) compared to the main element (see Fig. 8d), promote the formation of an outer peak larger than the inner one (see Fig. 10f).Additionally, the stronger APG on the flap surface results in significantly higher velocity fluctuations compared to the main wing.Further away from the wall, i.e., at higher values of y + , an increase of the Reynolds stresses can be anticipated (not shown here).This is analysed later on Sec.3.4 and is associated with the signature of the main and slat wakes.Note that on the flap, the boundary layer thickness is shifted towards smaller heights ( y + < 500 ) due to the confinement by the presence of the main wake.

Shear Layers
As illustrated by the streamlines depicted in Fig. 11, the cavities located in the slat and main elements promote the formation of recirculation bubbles.In the case of the slat cove, its shape is directly related to the position of the stagnation point, which moves towards the slat cusp as the angle of attack increases, resulting in a smaller recirculation bubble.This is visualized in Fig. 11a, c and e, where the streamlines in the slat cove region are displayed.In the main wing cavity, the shape is less dependent on the angle of attack, and its size remains relatively constant across the different angles of attack studied (see Fig. 11b, d  and f).The large dimensions of the main element constrain the flow downstream, and as approaching the cavity, flow is mainly streamlined along its pressure side.Nevertheless, despite showing similar flow directions, the flow conditions in the shear layer are not the same for all the angles of attack, i.e., the local velocities and fluctuations.Indeed, larger Moffat eddies (Moffatt 1964) are observed at increasing angles of attack.These are a series of vortices with decreasing size and intensity that are generated in corners by a source of motion away from the wall.As getting close to the corner these vortices start to develop independently when the Stokes flow is reached, i.e., Re << 1 .However, their absolute size is determined by the conditions far from the wall, where the stirring force is produced.While at = 5 and 9 • only a sequence of two eddies can be detected, at = 23 • three consecutive eddies can be identified.More eddies should be developed as approaching the corner, but the mesh resolution should be finer to visualize them.
A more quantitative analysis is performed in Figs. 12 and 13.These figures show the velocity magnitude U and turbulent kinetic energy TKE profiles at different positions along the slat and main shear layers, respectively.Namely, the positions are defined to be the 20, 40, 60 and 80% of the shear layer total length in each case and are referred to as S1, S2, S3 and S4 for the slat; and M1, M2, M3 and M4 for the main shear layers, respectively.
As observed in Fig 12a, the initial stage of the slat shear layer (S1) exhibits higher velocities in the region outside the cavity ( d∕c > 0 ) as the angle of attack is reduced.Due to the location of the stagnation point at low angles of attack, flow undergoes a higher acceleration as it moves along the slat surface towards the cusp.The steep velocity gradient near the cusp gives rise to natural Kelvin-Helmholtz (K-H) instabilities, with lower intensities observed as the angle of attack increases.Precisely, the excited frequencies in a mixing layer are related to the velocity magnitude on both sides ( U 1 and U 2 ), and the velocity gradient across it (dU/dn) as f 0 = (1∕2)(U 1 + U 2 )∕(7 0 ) , where 0 = |U 1 − U 2 |∕max(dU∕dn) .Considering the Strouhal number to be the non-dimen- sional frequency St = fc∕U ∞ , the previous equation leads to approximately St ≈ 60 at = 5 and 9 • , and St ≈ 37 for the predictions at = 23 • .These values almost match with the frequencies observed.For more details on the shear layer frequencies, the reader is referred to Montalà et al. (2024).
Progressing along the shear layer (S2, S3 and S4), the steep gradients begin to vanish, denoting the thickening of the shear layer (see Fig. 12b, c and d).This is associated with the appearance of three-dimensionalities on the flow and the onset of turbulent conditions, which increases the mixing across the shear layer.These effects are also showcased in the turbulent kinetic profiles.In Fig. 12e, very narrow peaks are observed at the beginning of the shear layer (S1), while wider profiles are obtained as progressively moving downstream (S2, S3 and S4), i.e., Fig. 12f, g and h.Also note that the magnitudes of turbulent kinetic energy peaks are amplified downstream, which is associated with the turbulent structures that are developed.
For the shear layer located in the main element cove, its profiles are displayed in Fig. 13.Compared to the slat (see Fig. 12), reduced velocity and TKE magnitudes are obtained here.It is worth pointing out that the M1, M2, M3 and M4 stages are not in the same relative position of the shear layer for all the angles of attack considered.These stages are obtained by considering that all shear layers start from the cusp of the cavity while this is not entirely true.As discussed in Montalà et al. (2024), at low angles of attack, there is a turbulent separation of the flow along the pressure side of the main element, which anticipates the formation of the shear layer at = 5 and 9 • .On the other hand, at = 23 • , flow remains attached and separation is produced when flow reaches the cavity.This explains the steeper velocity gradients observed in Fig. 13a at = 23 • , compared to the other angles of attack, as well as the higher TKE peaks detected in Fig. 13e at = 5 and 9 • in the initial stage of the shear layer (M1).In the lower angles of attack, due to the early separation of the flow, the levels of TKE are already high when flow reaches the cove.At = 23 • , the flow undergoes a turbulent transition within the shear layer inside the cavity, leading to progressively increasing TKE magnitudes along it, as shown in Fig. 13e, f, g and h.
As mentioned earlier, larger Moffat eddies are identified when the angle of attack increases (see Fig. 11b, d and f).The earlier development of K-H at low angles of attack, due to the flow separation, leads to the earlier onset of turbulent conditions.This, combined with the higher velocities of the outer stream ( d∕c > 0 ), enhance the momentum mixing from the outer to the inner regions of the cavity, promoting the formation of vortical structures inside it and resulting in larger velocity magnitudes ( d∕c < 0 ).This phenomenon is evidenced in Fig. 13a and e, where the velocity and TKE profiles for the initial stage of the shear layer (M1) are visualized.Therefore, this higher momentum observed at = 5 and 9 • may disfavour the formation of the Moffet eddies (delays the appearance of the Stokes condition), in contraposition to the lower velocities observed at = 23 • inside the cavity.

Wakes
In this section, the wakes behind the slat and main elements are analysed.In Fig. 14, the velocity and TKE profiles in the slat wake are inspected six different locations ( x∕c = 0.02 , 0.04, 0.06, 0.08, 0.010 and 0.12).The location of the wake can be traced by the velocity deficit and the TKE peaks detected in the middle of the profiles.It is observed that increasing the angle of attack deflects the slat wake away from the main element wall.As explained at the beginning of Sec.3.2, the turbulent coherent structures from the slat wake sweep the boundary layer of the main element, locally reducing the skin friction coefficient.The deflection of the slat wake leads to a lower reduction in the C f distribution (see Fig. 5b, x∕c ≈ 0.170 ) as fewer impingement events are produced.
The slat wake results from the interaction between the flow from the slat lower and upper sides as they meet at the trailing edge.Regarding the flow convected from the slat cove (lower part of the profiles in Fig. 14), the first stage ( x∕c = 0.2 ) exhibits increasing velocities with the angle of attack, alongside smaller values of TKE.This is associated with the location of the main element stagnation point (see Fig. 4), where the flow experiences higher acceleration through the gap betweeen the slat and main elements as the angle of attack is increased.This, combined with the smaller recirculation regions formed in the slat cavity, leads to smaller TKE values (see discussion in Sec.3.3 Fig. 12e, f, g and h).
Focusing on the region of flow convected from the suction side (upper part of the profiles in Fig. 14) at x∕c = 0.2 , the velocity and TKE magnitudes increase with the angle of attack.The increased TKE values are related to the development of a TBL along the slat suction side.In this element, the larger velocities at higher angles of attack translates into higher local Reynolds numbers in the suction side at = 23 • , which early triggers the tran- sition from laminar-to-turbulent flow (see Fig. 3c).
Further downstream along the different stages in Fig. 14, wake rapidly dissipates from the velocity field at the lower angles of attack, while it is still visible at x∕c = 0.12 for the case at = 23 • .Regarding the TKE, a two-lobbed profile is obtained at = 23 • , probably due to the footprint of both separated flows coming from the slat upper and lower sides.At the lower angles of attack, i.e. = 5 and 9 • , the TKE distribution exhibits a single peak with a similar magnitude than at = 23 • .This single bump is attributed to the velocity fluctuations primarily occurring in the slat cove, since no significant turbulent conditions are identified on the slat suction side at these angles of attack.
Similarly to the slat wake, Fig. 15 shows different velocity and TKE profiles of the flow around the flap at six different stages ( x∕c = 0.9 , 1.0, 1.1, 1.2, 1.3 and 1.4).These profiles allow the visualization of the slat, main and flap wakes.
The footprint of the slat wake is only visible at = 23 • ; its signature in the velocity fluctuations is still discernible despite the TKE levels being relatively small.In the other cases, i.e., = 5 and 9 • , the TKE levels are nearly zero and the slat wake has completely vanished from the mean flow.
Also from Fig. 15, increased levels of turbulent intensities are detected in the region where the main element wake is found.This can be identified for all the angles of attack, despite being more prominent at = 23 • .As commented in Sec.3.2, the stronger APG Fig. 14 a Velocity and b turbulent kinetic energy profiles at x∕c = 0.02 , 0.04, 0.06, 0.08, 0.10 and 0.12 leads to a higher momentum loss along the main wing suction side.As in the slat wake (see Fig. 14), two TKE peaks are detected in the main element wake due to the confluence of the flows from the main suction and pressure sides.At x∕c = 0.9 , the lower peak linked to the main element wake is considerably higher due to the accelerated fluid from the main-flap gap.As moving downstream, the turbulent mixing promotes the merging of both peaks, leading to a wider zone of increased TKE levels.This effect can barely be seen at lower angles of attack, where TKE levels are attenuated as the flow approaches the flap trailing edge.At = 23 • , off the flap surface, a large velocity deficit and increased TKE levels are observed.This velocity deficit leads to a visible recirculation area, which leaves a prominent wake downstream the flap (see Fig. 4c).
At the lower angles of attack, i.e., = 5 and 9 • , the flow separation near the flap trailing edge also induces a velocity deficit and increased TKE levels (see for instance x∕c = 1.2 and 1.3 in Fig. 15).However, the flap wake rapidly dissipates downstream.

Conclusions
Large eddy simulations of the flow past the three-element wing 30P30N are performed at a constant Reynolds number ( Re c = 750, 000 ) and three different angles of attack ( = 5 , 9 and 23 • ).Interesting features regarding the aerodynamic coefficients have been observed.The lift coefficient variations with the angle of attack are closely related to the changes in the pressure coefficient distribution along the wing surface.The main element, characterized by a larger wet area and higher suction peaks in the pressure coefficient, constitutes the predominant contributor to the total lift.Moreover, more pronounced suction peaks are identified as the angle of attack increases, leading to higher lift coefficients.Similar effects are observed in the slat, whereas the flap shows small variations of the lift coefficient as the angle of attack is modified.At the highest angle of attack, i.e., = 23 • , the contribution of the pressure coefficient to the lift coefficient on the main element is reduced, while its impact gets larger on the drag coefficient.This is associated with the onset of stall conditions, which is also evidenced by the spread of the streamlines above the flap.
Regarding the drag coefficient, the pressure and viscous contributions are analysed separately, demonstrating that the drag is mainly caused by pressure forces, with the main wing being the principal contributor to it.The total drag increases with the angle of attack, this being associated with the growth of the pressure drag, rather than with variations in the viscous drag, which is kept nearly constant along the different angles of attack.This is also noticeable from the skin friction distribution, which shows small differences with the angle of attack.In the slat, due to the increase of the suction peaks at higher angles of attack and its inclination with respect to the inflow, a negative pressure drag coefficient is obtained (pointing towards the upstream direction).Nevertheless, the total drag coefficient is still positive due to the higher contributions from the main and flap elements.At = 23 • , important flow features are brought about as a consequence of the increase in the APG along the main wing surface.This increase is translated into a more pronounced growth of both the TBL thickness and the momentum thickness, representing an increase of the wall-normal momentum to the detriment of the streamwise velocities.This effect can also be identified by the appearance of an outer peak in the streamwise Reynolds stresses as approaching the trailing edge, together with an increase of the wall-normal and shear stresses, denoting larger turbulent structures in the outer region of the TBL.
The development of a more prominent main element wake influences the evolution of the flap TBL.Due to the inclination angle of the flap relative to the main wing, higher values of the Clauser pressure gradient parameter and shape factor are observed at the lower angles of attack, i.e., stronger APG.Nonetheless, at the highest angle of attack studied, i.e., = 23 • , the spread of the main element wake above the flap attaches the flow transported from the gap between the main and flap elements to the flap wall, reducing the APG at this angle of attack and preventing the flow separation seen at the other angles of attack.
Due to the cavity-like shape of the slat and main coves, two recirculation bubbles are identified in these regions, bounded by the respective shear layers.In the slat cove, reduced recirculation areas are observed as the angle of attack increases, accompanied by a reduction of the TKE levels along the shear layer.In the main cove, the size of the recirculation bubble remains practically the same for all the angles of attack considered.Nevertheless, at the lower angles of attack, i.e., = 5 and 9 • , the formation of the shear layer is anticipated due to the flow separation on the main element pressure side as it approaches the cavity.This leads to higher TKE values downstream.At = 23 • , the transition to turbulence occurs inside the cavity and hence, the TKE levels increase as moving downstream along the shear layer.At this angle of attack, the low turbulent conditions found at the initial stages of the shear layer favour the formation of Moffat eddies in the corner of the cove.
It should be noted that, despite an anticipated trend between angles of attack, the above conclusions hold true only for the values studied in this work, and estimations made for other angles of attack within this range should be treated with due respect.

Fig. 1
Fig. 1 Computational domain.Green dashed lines represent inflow conditions; blue dashed lines represent outflow conditions.a General view and b Zoomed view

Fig. 2
Fig. 2 Two-point correlations along the spanwise direction at = 9 • .a Location of the probes coloured accordingly, b x-velocity (u) correlations and c y-velocity (v) correlations

Fig. 5 a
Fig. 5 a Pressure coefficient and b skin friction coefficient (main and flap suction sides only)

Fig. 7
Fig. 7 Drag coefficient at different angles of attack.a Pressure drag coefficient and b viscous drag coefficient

Fig. 8
Fig. 8 Boundary layer development along the main suction side.a Tangential velocity profiles (solid lines) and boundary layer thickness (dotted lines), b momentum thickness Reynolds number, c shape factor and d Clauser pressure gradient parameter

Fig. 9
Fig. 9 Boundary layer development along the flap suction side.a Tangential velocity profiles (solid lines) and boundary layer thickness (dotted lines), b momentum thickness Reynolds number, c shape factor and d Clauser pressure gradient parameter

Fig. 10
Fig. 10 Tangential (solid lines), wall-normal (dashed lines) and shear (dash-dotted lines) Reynolds stresses at different locations along the a, b, c main and d, e, f flap suction sides

Fig. 12
Fig. 12 Velocity (top row) and turbulent kinetic energy (bottom row) profiles at different locations along the slat shear layer path.a, e S1; b, f S2; c, g S3; and d, h S4

Fig. 13
Fig. 13 Velocity (top row) and turbulent kinetic energy (bottom row) profiles at different locations along the main shear layer path.a, e M1; b, f M2; c, g M3; and d, h M4