Numerical Simulation of a Passive Control of the Flow Around an Aerofoil Using a Flexible, Self Adaptive Flaplet

Self-activated feathers are used by almost all birds to adapt their wing characteristics to delay stall or to moderate its adverse effects (e.g., during landing or sudden increase in angle of attack due to gusts). Some of the feathers are believed to pop up as a consequence of flow separation and to interact with the flow and produce beneficial modifications of the unsteady vorticity field. The use of self adaptive flaplets in aircrafts, inspired by birds feathers, requires the understanding of the physical mechanisms leading to the mentioned aerodynamic benefits and the determination of the characteristics of optimal flaps including their size, positioning and ideal fabrication material. In this framework, this numerical study is divided in two parts. Firstly, in a simplified scenario, we determine the main characteristics that render a flap mounted on an aerofoil at high angle of attack able to deliver increased lift and improved aerodynamic efficiency, by varying its length, position and its natural frequency. Later on, a detailed direct numerical simulation analysis is used to understand the origin of the aerodynamic benefits introduced by the flaplet movement induced by the interaction with the flow field. The parametric study that has been carried out, reveals that an optimal flap can deliver a mean lift increase of about 20% on a NACA0020 aerofoil at an incidence of 20o degrees. The results obtained from the direct numerical simulation of the flow field around the aerofoil equipped with the optimal flap at a chord Reynolds number of 2 × 104 shows that the flaplet movement is mainly induced by a cyclic passage of a large recirculation bubble on the aerofoil suction side. In turns, when the flap is pushed downward, the induced plane jet displaces the trailing edge vortices further downstream, away from the wing, moderating the downforce generated by those vortices and regularising the shedding cycle that appears to be much more organised when the optimal flaplet configuration is selected.


Introduction
The control of flow separation in wings at high angle of attack has been the focus of many research activities in the past. In particular, a number of biomimetic methodologies for separation control on wings in highly loaded conditions have been inspired by observing the flight or swimming characteristics of certain birds and fish [1][2][3][4]. In particular, the idea of reproducing the pop-up of birds feathers for stall delay and control is becoming increasingly popular because of its passive but still self-adaptive character: the feathers lift up is believed to be induced by the back-flow occurring when the flow separates as a consequence of the increased angle of attack [3,5,6].
The results reported by Schatz et al. [7] show that the use of self deploying flexible flaps, mounted on the suction side of a wing (HQ17 aerofoil), can deliver an increase in lift of about 10% in nominally stalled conditions at a chord Reynolds number, Re c 10 6 (Re c = U ∞ c/ν is the Reynolds number based on the magnitude of the free stream velocity U ∞ and the aerofoil chord c). Unsteady Reynolds-averaged numerical simulations were used together with experimental measurements to describe the mechanism that produce the added lift; however, due to the lack of information available from this kind of simulations, further study is necessary to fully undestrand this compex fluid-structure interaction problem. More recent experiments by Schluter [8] have also studied the effectiveness of an adaptive passive flap on a SD8020 aerofoil at moderate Reynolds number (Re c = 3 − 4 × 10 4 ) showing that its use promotes a lift increase in near stall conditions. Wang and Schluter [9] have extended the previous analysis to genuinely three-dimensional conditions considering the effects of a passive flap on a wing of finite span with the same aerofoil section. They found that the flap still deliver a substantial lift benefit if extending over the 80% of the wingspan leaving the tip clear. They also observed that the position and the length of the flap leading to improved aerodynamic performances were independent of the three dimensional character of the flow field. Bechert et al. [10] have extensively investigated the effects of wing mounted movable flaps in a series of wind tunnel experiments. Their results indicate that adaptive flaps show good aerodynamic performances on wings with a large aspect ratio by successfully suppressing flow separation that develops gradually upstream from the trailing edge. Traub and Jaybush [11] have systematically evaluated the effect of several self-actuated 3D spoiler geometries using wind tunnel experiments at Re c = 2.25 × 10 5 on a SD7062 profile. The best results, in terms of largest lift increase in quasi stalled conditions, were obtained when considering a square slotted spoiler. Bramesfeld and Maughmer [5] explored the effect of small, movable tabs mounted on the suction side of a S824 aerofoil in a low-speed wind tunnel experiment conducted at Re c 10 6 . From the surface pressure distributions they discovered that the effectors act as pressure dams that reduce the adverse effects of the separation, allowing higher pressures upstream of their location. Johnston et al. [12] and [13] made a comparison of the effectiveness of free-moving and fixed flaps mounted at different deployment angles over an angle of attack range from 12 o to 20 o , and found a similar behaviour in term of lift, with the maximum lift obtained for deployment angle less than 60 o . However, they found that the fixed flap produces more drag than the free-moving one. Recently, Bruecker and Weidner [14] used flexible flaps to delay the dynamic stall lift breakdown of a NACA0020 wing at moderate Reynolds number (i.e., Re c = 7.7 × 10 5 ) in ramp-up motion (α 0 = 0 and α s = 20 o ). The authors also offer a mechanistic explanation of the stall delay that would be due to a reduction of the backflow, and by a re-organisation of the shear layer roll-up process. In turns, the modified roll-up pattern would cause a delay in the onset of the non-linear growth of the shear layer via a mode-locking of the fundamental instability mode with the motion of the flaps. In disagreement with the majority of the research community, Kernstine et al. [15] found that the highest increase in lift, on separation onset, was obtained with a flap mounted in the first half of a NACA2412 aerofoil, slightly downstream of the leading edge. Very recently another parametric study on the geometry and location of the flap was performed by Altman and Allemand [16]. Their experiments could not confirm the best configuration suggested by Kernstine et al. [15]. More in general, the authors conjecture that it might not be possible to design a universal flap configuration improving post-stall performances.
Apart from the aerodynamic improvements offered by adaptive flaps in stall conditions, the use of similar devices has also been explored as a method for reducing structural vibrations in aerofoils. Liu et al. [17] and Montefort et al. [18] have investigated the effects of a single flexible, polymeric rectangular flap and of an array of small rectangular polymeric flaplets attached near the leading edge on the upper wing surface, considering a NACA0012 aerofoil and a flat-plate. They found that by manipulating the unsteady structure of the flow, these devices were able to reduce significantly wing vibrations particularly near the dominant first torsional mode.
The present contribution will just focus on the impact that self-adaptive very thin flaps have on the flow field structures around a wing at high angle of attack. In particular, the possibility of controlling the flow around a NACA0020 aerofoil using passive, selfadaptive, almost zero-thickness flaps attached to the suction side of the aerofoil will be explored performing a series of Direct Numerical Simulations. After having reported the results of a preliminary parametric study meant to bound the characteristics of the best performing geometries and locations, the attention will move on a detailed analysis of the three-dimensional flow field generated by the wing at α = 20 o degrees when a quasi optimal flaplet is mounted on its suction side at Re c = 2 × 10 4 . By carrying out an in-depth analysis of flow fields generated by direct numerical simulations, we will characterise the main effects induced by the presence of the flap and we will also propose a conceptual explanation of their effectiveness in delivering aerodynamic benefits in stalled configurations. This works differs from previous research on the topic, due to the presence of a torsional spring, holding the flap to the aerofoil surface. By adjusting the value of the torsional stiffness, the flap movement can be selectively locked-in a specific flow frequency, thus introducing an additional dynamic mean of manipulating the flow . The paper is structured as follows. Initially in Section 2 we will give a brief overview of the numerical methods and of the geometrical set-up that have been used. Then, in Section 4 we will illustrate the results of the preliminary two-dimensional parametric campaign that we have carried out to roughly identify the quasi optimal configuration and location of the flaplet. Finally, in Section 5 the results and the interpretation of the flow fields generated by a full direct numerical simulations are offered also by comparing the characteristics of the fields obtained with and without flaplet. Some conclusions will be drawn at the end of the paper in Section 6.

Baseline Numerical Formulation
To tackle the problem at hand, we consider an incompressible three-dimensional unsteady flow field, governed by the Navier-Stokes equations around a straight wing with an infinite spanwise dimension. The computational domain is shown in Fig. 1. The coordinate system is a Cartesian inertial one, with the x and y axis (x 1 and x 2 ) denoting the directions parallel and normal to the aerofoil chord, and z (x 3 ) being the axis normal to the paper. Also, u, v and w (u 1 , u 2 and u 3 ) denote the velocity components parallel and normal to the chord, and along the span respectively. With the given notations, using Einstein's summation convention, the dimensionless equations that govern the incompressible flow motion are: The equations have been made non-dimensional using the magnitude of the free stream velocity U ∞ and the aerofoil chord c. Also, in the momentum Eq. 1, Re c = U ∞ c/ν is the Reynolds number and f i represents a system of body forces used to keep into account the presence of the flap as it will be discussed later. The momentum and mass conservation Eqs. 1 and 2, are discretised on a cell-centred, co-located grid using a well-established curvilinear finite volume code [19][20][21][22]. The fluxes are approximated by a second-order central formulation, and the method of Rhie and Chow [23] is used to avoid spurious pressure oscillations. The equations are advanced in time by a second-order semi-implicit fractional-step procedure [24], where the implicit Crank-Nicolson scheme is used for the wall normal diffusive terms, and the explicit Adams-Bashforth scheme is employed for all the other terms. The pressure Poisson equation arising when imposing the solenoidal condition on the velocity field, is transformed into a series of two-dimensional Helmholtz equations in wave number space via Fast Fourier transform (FFT) in the spanwise direction. Each of the resultant elliptic 2D problem is then solved using a preconditioned Krylov method (PETSc library [25]). In particular, the iterative Biconjugate Gradient Stabilised (BiCGStab) method preconditioned by an algebraic multigrid preconditioner (boomerAMG) [26] revealed to be quite efficient. The code is parallelised using a streamwise domain decomposition via the MPI message passing library. Further details on the code, its parallelisation and the extensive validation campaign that has been carried out in the past can be found in Rosti et al. [21].
The aerofoil that has been selected for the present study is a symmetric NACA0020, which has been extensively studied at static and dynamic stalled conditions by the authours [21,27]. The flow domain around the aerofoil is meshed using a body fitted C grid arrangement (see Fig. 1b). The grid is adapted to the three dimensional case by repeating the baseline 2D grid uniformly in the spanwise direction. With this arrangement, the external surface that bounds the computational domain, contains both the inlet and the outlet (see Fig. 1a). To determine which portions of the external bounding surface act as an inlet (or an outlet), at each time step a local spanwise average of the fluid velocity is evaluated in a tiny region close to the boundary. When the averaged flow direction points outward, the corresponding portion of the boundary is assumed to be an outlet, and is treated using a convective boundary condition. Conversely, if the mean flow direction is directed inward, the corresponding boundary surface is considered to be an inlet, and a Dirichlet type condition based on an irrotational approximation is used. In particular, the values to be assigned to the velocity field on the Dirichlet portions of the boundary are determined by solving a companion potential equation discretised via the panel method of Hess and Smith [28].
The remaining boundary conditions are imposed as follows: impermeability and no slip conditions are set on the aerofoil wall, periodic conditions are assumed on the planes bounding the domain in the spanwise direction, and continuity of the flow variables is enforced through the top and bottom planes generated by the C-grid topology downstream of the trailing edge.
All the three-dimensional simulations that will be presented have been obtained at a chord Reynolds number Re c = 20000. Differently, the two-dimensional parametric study that will be presented in the next section has been carried out at Re c = 2000. For both the 3D and the 2D simulations, the angle of attack has been set to α = 20 o (stalled condition).
The grid system that has been chosen for the three dimensional simulations, has been determined after a number of trial and errors tests and companion grid convergence studies. Finally, we have found that a grid composed by 2785 × 626 × 97 nodes (in the x 1 , x 2 and x 3 directions, respectively), delivered a sound compromise between all the local resolution requirements set by the imposed high angle of attack: wing curvature, separation, attached turbulent boundary layers and shear layers embedded in the flow field. In terms of local wall units, in the attached turbulent layers, the corresponding mesh resolution verifies x + < 3.0, y + < 0.5 and z + < 7.5 (superscript + indicates standard local viscous units lengths: i.e., lengths made non dimensional using the kinematic viscosity ν, and the skin friction velocity u τ ). Finally, the spanwise size of the domain has been set equal to 0.9c, which guarantees a good velocity decorrelation between the periodic end planes [21].
Further details on the procedure that has been followed to generate the grid and the mesh refinement study campaign can be found in Rosti et al. [21]. Figure 2 shows the configuration that we have analysed in this study, it comprises a NACA0020 aerofoil with a rigid, nominally infinitely thin flaplet of length L mounted on the wing suction side. The flap is hinged to the surface via a torsional spring that constraints its motion to take place on the x − y plane. The evolution of the flap angular displacement, θ(t) can be modelled using the canonical second order differential equation:

Fluid-flap interaction model
In Eq. 3, I is the flap moment of inertia with respect to the rotation axis (i.e., I = mL 2 /3, m being the mass of the flap per unit spanwise length), C is an angular damping factor and K is the spring rotational stiffness (C and K are per unit spanwise length too). Finally, T is the total torque per unit spanwise length exerted by the fluid forces on the flap. When no damping is considered, a compact way to characterise the physical properties of the flap is based on specifying the spring stiffness K in terms of the moment of inertia I and its natural frequency f , obtained from the solution of the homogeneous equation associated to Eq. 3: The coupled motion of the flap and the surrounding fluid are enforced using an Immersed Bounday Method (IBM) [29][30][31][32]. In particular, at each time step, the presence of the flap is modelled by introducing a system of singular forces f i distributed along the flap and appearing as body forces in the momentum Eq. 1. In particular, this body force distribution is computed to impose the impermeability and the non-slip conditions on each instantaneous flap configuration determined by its angular position θ(t). On the other hand, at each time step, the integral of the elemental contributions of each force f i to the linear momentum balance about the hinge provides the total torque forcing Eq. 3. The coupled algorithm that we have briefly described is based on a particular version of the immersed boundary (IB) method (i.e., the Reproducing Kernel Particle Method -RKPM method developed by Pinelli et al. [32]) that will be shortly revised hereafter.
In common with many others IB algorithms, the first stage of the algorithm involves a discretisation of the immersed body by distributing N nodes X i , i = 1, . . . , N (termed as Lagrangian points) over the surface bounding the immersed object. Generally, this set of nodes do not coincide with the underlying body fitted grid x i,j,k used to discretise the domain around the baseline aerofoil. This geometric discrepancy introduces the necessity of having a tool able to transfer the body forces defined on an immersed surface into equivalent forces defined over a local volume surrounding the surface but belonging to the body fitted grid. A discussion on the effect of spreading forces from the immersed surface into a volume comprising nodes of the C-grid used to mesh the aerofoil is posticipated at the end of this section. In the general framework of a pressure correction method, the second stage of the IB procedure involves a preliminary time advancement of the momentum equations without considering the presence of the immersed surface. The obtained predicted velocity field u * (x i,j,k ) is then interpolated onto the embedded surface : U (X i ) = I(u * ), where the velocity corrections leading to the prescribed velocity distribution on the surface U are computed. These velocity corrections, per time unit, can be interpreted as system of local body forces that restore the desired boundary conditions on : In the final stage of the IB method, the previously obtained velocity field u * (x i,j,k ) is discarded and the momentum equations are advanced again using the boundary restoring forces obtained in the predictive stage (4). This force is evaluated on the fluid grid x i,j,k from the values at X i using a pseudo inverse of the operator I, indicated with C and termed as spread.
The spread operation formally allows to determine the singular forces on the fluid grid as: Aside from the flow field time advancement, also the position of the flap needs to be updated. Once the torque in Eq. 3 is computed by integrating each contribution of the singular forces F * (X i ) along the whole flap (4), the new angular position θ(t) is found by integrating (3) in time. Finally, all the flap Lagrangian coordinates, and their respective velocities are updated consistently with a rigid body rotation about the hinge. The global time advancement scheme finalises with the solution of a pressure Poisson equation and the final projection of the velocity field onto the consistent divergence free space. A summary of the basic steps involved in the algorithm used to advance in time the fully coupled flap-fluid system is provided in Fig. 3.
The distinguishing feature of any continuous forcing IB method is related to the way in which the two operators I and C are applied. In our case, we follow the RKPM approach, used by Liu et al. [33,34] and Zhang et al. [35], to construct a quasi Dirac's delta function that can be defined on arbitrary supports [32]. The derived mollifier shares a number of momentum properties with a genuine delta function and therefore can be used to approximate both the interpolation and the spreading (i.e., convolution) operators as it would be formally done with a delta function in a distribution sense. As an example, the approximation f a (x) of the value of a given smooth function at point x ∈ can be expressed as a convolution with a kernel function w d having the first moments of a genuine Dirac's function as: In particular, the RKPM method assembles the kernel w d by modifying a compact support weight function w with a correcting polynomial p which coefficients are found by matching moments with a real delta function: Because of the finite size of the support over which each regularised delta function acts to enforce the boundary conditions on the immersed object, the latter inherits an effective aerodynamic thickness. This thickness is related to the actual size of the support that in turns is determined by the local mesh size (typically the effective thickness is equivalent to the diagonal of a computational cell). In summary, whilst the flow around the baseline aerofoil is simulated using a classical body fitted, C-grid, the effects on the flow generated by the flaplet and the dynamic of the latter are kept into account via an immersed boundary method. In particular, the movement of the flap is determined via the integral of the fluid torques distributed along the flap itself. The local torques (per unit mass) are generated by the local accelerations computed by imposing the desired flap velocity at each time instant and the distance of each Lagrangian node to the flap hinge. The local accelation are obtained by interpolation from the body fitted grid into the immersed surface. On the other hand, the same singular force distribution is used as a set of body forces on the right hand side of the fluid momentum equations discretised on the body fitted grid. The operation to transform local forces distributed on a surface into a set of forces operating on a volume strip belonging to the body fitted mesh is carried out via a finite support compact pseudo Dirac's delta function. As a consequence the actual shape of the flap is not seen as a sharp object by the fluid flow, but rather as a diffused volume strip with a finite (i.e. a non-zero) aerodynamic thickness (seen by the flow). Further details on the RKPM IB method and its implementation in a finite volume context can be found in Pinelli et al. [32]. For an implementation in a Lattice Boltzmann framework including moving and deformable surfaces the reader can refer to Favier et al. [29].

Baseline Flow Characterisation
To introduce the main features of the flowfield that we wish to manipulate, we consider a NACA 0020 aerofoil at an angle of incidence of 20 o and at chord Reynolds number of 2 × 10 4 . In these conditions the flow is mainly characterised by a large recirculation zone covering almost the whole suction side as shown in Fig. 4a [21]. Moreover, both a secondary counter rotating vortex located by the trailing edge, and another very small recirculation bubble close to the aerofoil maximum thickness can also be observed. All the mentioned spanwise vortices are enclosed within a region bounded by the two shear layers originating at the leading and trailing edges. The leading edge shear layer is induced by the early separation of the free stream laminar flow approaching the wing (see Fig. 5) and by the subsequent  [21]. The unsteady behaviour of the spanwise vorticity field, determined by the shear layers instabilities and by the mutual interaction of the vortices embedded in the wake, is the ultimate responsible of the aerodynamic response of the aerofoil to stalled conditions. For this reason, any control strategy that aims at an overall improvement of the aerodynamic efficiency must tackle the direct manipulation of the vorticity field and its unsteadiness. Along this line of thought, this work investigates on the possibility of controlling the vorticity field generated by an aerofoil at high angle of attack using a self adaptive flaplet mounted on its suction side. In particular, the objective is to find a configuration that palliates the detrimental effects of stall by producing increased lift. To pursue such an objective, a parametric study covering a fully three-dimensional flow at the targeted chord Reynolds number would be computationally unrealistic. For this reason, a preliminary study on a low Reynolds number, fully laminar, two-dimensional flow has been Contours of instantaneous flow streamwise velocity and streamlines. Contour goes from −0.3U ∞ (blue) to 1.4U ∞ (red), and the snapshots cover a full shedding period of 1.87c/U ∞ carried out with the objective of bounding the parametric range that needs to be explored for achieving a good flap design in a realistic three dimensional scenario. Before describing the initial two-dimensional parametric study, a comparison between the two baseline cases (i.e., fully 3D case at higher Reynolds number versus the laminar case at lower Reynolds numbers) will be introduced to provide a conceptual justification of the procedure that has been followed. Figure 6 compares the character of the mean three dimensional x−wise velocity field at Re c = 2 × 10 4 and α = 20 o with the two dimensional field obtained at the same angle of attack but at one order of magnitude smaller Reynolds number, i.e., Re c = 2 × 10 3 . The two velocity fields show similar qualitative features: large recirculating regions of comparable magnitude covering the whole suction side of the aerofoil (i.e., the sizes of the recirculating regions are 0.5c and 0.35c in the 2D and 3D case, respectively). In both cases, the flow separates at the leading edge (x s ≈ 0.025) reattaching at x r ≈ 0.9 in the 2D case, whilst staying detached along all the rest of the suction side for the 3D case. The unsteadiness of both the 2D and the 3D stalled cases is mainly determined by the presence, the interaction and the shedding of the two large counter rotating vortices that characterise the region above the aerofoil (see Fig. 7). The dynamic of these two large vortices governing the lift oscillations, is mainly of 2D, laminar nature and basically involves only the interaction of the very large coherent structures embedded in the flow. Although the quantitative differences between the two-dimensional and the three-dimensional case are not negligible, the dominating effects and the events sequencing appear to be qualitatively similar. Moreover, since the self adaptive flap that we will use extends over the whole span of the wing, no significant 3D effects will be introduced by its presence as the flaplet will mainly interfere with the largest integral scales of the flow which are intrinsically two-dimensional in character.

Flaplet Design in 2D
Motivated by the aforementioned considerations, we have initially focused on the geometrical properties (i.e., size and location) and the flap dynamic response (i.e., its natural frequency) that deliver an optimal condition in a two dimensional, fully laminar flow at α = 20 o . Here, we define an optimal condition as the one that delivers the highest lift coefficient c l , preserving or improving the aerodynamic efficiency E = c l /c d .
We have started our analysis by considering the low Reynolds number (i.e., Re c = 2 × 10 3 ), 2D flow over a NACA0020 aerofoil at α = 20 o without any added flap. Figure 8 shows the time evolution of the lift and drag coefficients for the baseline configuration. Both coefficients are characterised by periodic oscillations: every period of lift coefficient    Fig. 9. The roll up of the trailing edge vortex, corresponds to a decrease in lift that gradually disappears as the vortex is shed into the wake. In the following time instants of the sequence Fig. 10, another pair of vortices is formed and shed away from the aerofoil. However, the newly generated lifting vortex quickly detaches from the wing surface, thus preventing the lift to raise. As the lift vortex is shed into the wake, it starts interacting with the trailing edge vortex that rolls up increasing its size. This interaction energises the trailing edge vortex with a consequent further decrease in lift, and with an impact in determining the structure of the near-wake (see Figs. 9 and 10). The final snapshots of the series, correspond to the end of the cycle with the generation of a new lifting vortex leading to the beginning of a new cycle.
Next, we have proceeded to perform a parametric study on the aerodynamic effects of the flaplet configuration. In particular, the flap reaction to the underlying unsteady flow field can be tuned by acting on various parameters: its length, position, inertia, spring stiffness and damping factor. The outcomes of the analysis conducted by varying the aforementioned parameters are summarised in Table 1 reporting some typical variations of the averaged aerodynamic coefficients (last three columns) when changing the flaplet characteristics (second to fifth columns). In particular, the length L of the flap was varied between 0.1c and 0.3c, the position of the flap hinge x F ranged between 0.6c and 0.8c (measured from the leading edge), the stiffness K of the spring was set such that its natural frequency f was between 1/4th and 4 times the shedding frequency f 0 of the baseline case without flap. The effects of the length and stiffness of the torsional spring on the value of the mean lift coefficient c l are also reported graphically in Fig. 12a. An optimum condition (i.e., maximum lift increase with respect to the baseline case) is achieved with a flaplet 0.2c long, resonating with the shedding frequency (flap natural frequency equal to the shedding one). Except for the cases of flaplets of very low natural frequency, if the latter doesn't match the baseline flow shedding frequency, the lift coefficient turns out to be almost unaffected by the length of the flap. On the other hand, when considering resonating conditions, the maximum lift and efficiency are achieved using a flaplet L = 0.2c long, a size roughly corresponding to half the dimensions of the recirculation region. Figure 12b shows how the lift coefficient  hinge location at 0.7c and a spring stiffness leading to a flaplet natural frequency matching the shedding one. For this specific flow condition and with the mentioned configuration, the flaplet interferes actively with the unsteady vorticity field delivering a 20% increase in the average aerodynamic efficiency. The corresponding time variations of the lift c l and drag c d coefficients are reported in Fig. 11 together with the elevation y of the tip of the flap from the surface of aerofoil. The time averaged c l is 35% higher than the case without flap (see Fig. 8), whilst the shedding frequency remains unchanged (i.e., f s = 0.555U ∞ /c). Performance of each configuration is evaluated by the lift coefficient c l , the drag coefficient c d , and the efficiency E = c l /c d . during a ramp-up manoeuvre. The aerofoil is NACA0020 and the Reynolds number is Re c = 2000. The flap parameters, i.e., the ratio between the spring natural frequency and the shedding frequency f/f 0 , the flap's length L, the hinge position x F , the spring rotational stiffness K, and the moment of inertia I is provided Also, the presence of the flap increases the r.m.s. of the lift coefficient by 15%. However, differently from the baseline case, the presence of the flap seems to regularise the shedding pattern, with all the lift extrema attaining almost the same value at each shedding period (see Fig. 8). The right columns in Figs. 9 and 10 show the spanwise vorticity over two shedding cycles at the times marked in Fig. 11a. In the initial snapshots Fig. 9, when the flap is almost laying on the aerofoil surface, a first vortex detaches from the trailing edge. Next, (Fig. 9f and h) the flap reaches its maximum elevation whilst a large lifting vortex is formed  Fig. 7). The circulation of the leading edge vortex which is responsible for the lift generation is only slightly increased by the presence of the flaplet (i.e., ≈ 3%), whilst the circulation of the trailing edge vortex, responsible for the generation of the downforce, is substantially reduced by a factor of ≈ 20%. Therefore, the increase in the average lift induced by the presence of the flaplet is mainly related with i) the regularisation of the shedding process and with ii) the reduction of the downward force induced by the trailing edge vortex.
This preliminary study conducted in a simplified 2D, laminar scenario has allowed to determine a point in the parameters space leading to a maximum increase in both lift and aerodynamic efficiency. The analysis has also characterised the features of the unsteady vorticity fields that develops when the optimal flap is used. The validity of our conjecture about the possibility of extending the results obtained with a simplified 2D scenario to a realistic 3D one will be discussed next.

Effect of the Adaptive Flaplet on a 3D Aerofoil
We now compare the three-dimensional flow fields around a NACA0020 at an angle of incidence of 20 o and at Re c = 2 × 10 4 , obtained when considering the unmodified aerofoil and when equipping the wing with a flaplet extending along its whole span, and featuring the quasi optimal configuration discussed in the previous section (flap length L = 0.2c, hinge location at x = 0.7c). Furthermore, inspired by the two-dimensional results, the stiffness of the torsional spring has been set to K = 0.150 leading to a natural frequency that matches the shedding one of the unmodified aerofoil.  Fig. 13b showing a net improvement when the flaplet is introduced with a mean efficiency growth from 1.8 (baseline case) to 2.0 (i.e., 11% increase with the flap). This improvement is in good agreement with the experimental results reported by Schatz et al. [7]. Furthermore, the time evolution of the aerodynamic coefficients clearly reveals the presence of a dominant frequency that corresponds to the shedding rate of the vortices into the wake. The introduction of the flaplet does not modify the value of the associated Strouhal number S t = f s c/U ∞ that remains fixed to S t = 0.534, a value that is almost the same as the one found in the 2D case at lower Reynolds number. When we compare the mean pressure coefficients C p of the two configuration, as shown in Fig. 14a, we notice that the pressure on the suction side of the aerofoil with flap is reduced upstream the flap position, thus generating a higher lift, in agreement with the results by Schatz et al. [7] and Bramesfeld and Maughmer [5], and then increases, downstream the location of its hinge. The friction coefficient C f , reported in Fig. 14, shows that the two aerofoils have a similar friction profile, with an early leading edge separation located at x ≈ 0.025c [21], except, in the leading edge peak which is enhanced by 10% in the case with flap.
Next, we analyse the effect of the flaplet on the average fields. We start by comparing the contours of the mean spanwise component of vorticity ω z in Fig. 15b. The figure shows that both the aerofoils are in a fully stalled condition with a large recirculation zone present on the whole suction side. Another smaller recirculation bubble is visible in both cases at about 0.25c from the leading edge, in proximity of the location of the aerofoil maximum thickness. The backflow region with positive vorticity (i.e., red: counter clockwise) on the suction side is clearly reduced when the flap is in use. Moreover, we can also notice that the presence of the flaplet reduces the size of the positive vorticity recirculating region by the trailing edge, also displacing the peak of positive vorticity further downstream, well beyond the trailing edge.
More information on the mean flow can be educed from the velocity profiles in Fig. 16 where the x and y velocity components are shown. Whilst the mean flow velocity on the pressure side is basically unaffected by the presence of the flaplet, on the suction side the velocity field changes in the region spanned by the flap movement. As compared to the baseline case, upstream of the flap location, at x = 0.6c, both velocity components are reduced in amplitude, with a corresponding overall reduction of reversed flow. Downstream  The effects of the flaplet on the flow become more visible when considering the distribution of higher order statistical quantities. Figure 17a, shows a comparison of the averaged turbulent kinetic energy k = 1/2 < u i u i >, in the controlled and uncontrolled cases. Consistently with the upstream laminar conditions, the kinetic energy is initially zero for both the aerofoils. Further downstream in the shear layer originated at the leading edge, k starts to grow similarly in both cases. On the other hand, the second shear layer formed past the trailing edge is influenced by the action of the flaplet. Its motion reduces the intensity of the velocity fluctuations. Downstream of the aerofoil, the two shear layers merge into the wake where the reduced levels of k, due to the flaplet action, are evident. This is clearly visible from Fig. 17b showing the turbulent kinetic energy profile at x ≈ 5.5c.
As previously mentioned, one of the consequences of the action of the flaplet on the flow field is the reduction in the intensity of the backflow on the aerofoil surface. To quantify this effect, in Fig. 18 we display the probability of finding a negative streamwise velocity component P (u < 0) in the two cases. In both situations, this probability is obviously zero in the outer flow where the u velocity is always positive, whilst its value increases in the recirculating region. In the reference case, the highest probability of backflow corresponds to the region close to the trailing edge, at x ≈ 0.8. In the case where the flaplet is active, the probability of having backflow is remarkably reduced not only in the region spanned by the flap movement but also upstream of it. To gain further insight on the effect of the flaplet-flow interaction we have used the classical Q-criterion proposed by Hunt et al. [36]. This technique assigns a vortex to all spatial regions that verify the condition where S = 1 2 ∇u + ∇u T is the strain rate tensor and = 1 2 ∇u − ∇u T is the vorticity tensor. Instantaneous Q iso-surfaces corresponding to the case without and with flaplet are shown in Figs. 19 and 20. From the first figure, it appears that the action of the flap contributes to the reductions of both the backflow and the generation of turbulent structures upstream of its location. Moreover, coherent structures in the wake dissipate faster in presence of the flaplet consistently with the drop observed in the profiles of turbulent kinetic energy Fig. 17. From the second figure, it is possible to recognise the principal flow features of the baseline case [21]. These are summarised hereafter to introduce the comparison with the flaplet case. Initially, the incoming laminar flow separates at the leading edge, forming a shear layer that rolls up into Kelvin-Helmholtz (KH) vortices [37][38][39][40][41]; this instability, locally, triggers the flow transition to turbulence; further downstream, the turbulent separated region appears to be characterised by fine texture, small-scale eddies, eventually merging into coherent larger structures; finally behind the aerofoil, a large turbulent wake is formed, the dynamics of which are similar to a von Karman vortex street typical of bluff body wakes. In contrast to classical vortex shedding process showing an alternately series of vortices of opposite sign and equal strength, here the wake is highly asymmetric presenting vortices of uneven strength. The loss of symmetry and the irregularity of the vortices pattern is related to the interaction between the two vortex generating mechanisms [42,43]:  To study the differences in the shear-layer characteristics between the baseline case and the optimal flap simulation, we carried out a Finite Time Lyapunov Exponent (FTLE) analysis. This technique is a Lagrangian coherent structures educing technique, see Haller [44] and Shadden et al. [45], that highlights the presence of strong shear layers in separated flows. The FTLE σ T (x, t) is a scalar function of space and time which measures the rate of separation of neighbouring particle trajectories initialised within a small ball centred at x at time t, and is defined as Here, λ max ( ) is the largest singular value of the Cauchy-Green deformation tensor computed over a finite time interval [t 0 , t 0 + T ] By looking at the time variation of the vorticity field another important effect of the interaction between the flow and the flaplet emerges. In particular, in Figs. 22e-f and 23, we compare the evolution of the spanwise vorticity ω z over two shedding cycles for both the cases, without (left column) and with flap (right column). The sequence of the reference case starts with the lifting vortex recently shed, and the trailing edge vortex being freshly formed and ready to be shed (Fig. 22a). As the lifting vortex detaches, another one is generated above the aerofoil (see Fig. 22), and eventually shed into the wake at a later stage (see Fig. 22) when the formation of the next trailing edge vortex takes place. The latter does not undergo a full evolution as it clearly appears from the following snapshots. In the following shedding cycle (see the left column in Fig. 23) the aforedescribed process almost repeats identically but with a remarkable difference: the trailing edge vortex is generated slightly more downstream than the previous one, thus allowing the new lifting vortex to expand more than its predecessor (see Figs. 22 and 23). The presence of the flaplet alters the previously described sequence. Here, the initial snapshot has been chosen to match the condition in which the flaplet lays on the aerofoil surface Fig. 22. In this situation, the trailing edge vortex has just been shed, and the lifting vortex is forming. As the flap lifts up Fig. 23 under the action of the pressure gradient induced by the passage of the lift vortex, a new trailing edge vortex is formed whilst the lifting vortex is shed away. As a consequence, the flap moves downward Fig. 23 under the action of the trailing edge vortex that is forming and subsequently, detaches from the trailing edge. The formation and roll up of the trailing edge vortex is conditioned by the movement of the flap that during its downward rotation generates a jet that pushes the vortex downstream. The displacement of the trailing edge vortex away from the aerofoil at every shedding cycle allows the incoming lifting vortex to grow and develop more freely without the constraint generated by the vicinity of a counter rotating vortex. The detachment of the trailing edge vortex induced by the flap generated jet has also a regularisation effect on the shedding cycle that now repeats identically with no difference between consecutive cycles. As the snapshots indicate, the position of the flap is strongly related with the passage of the lifting vortex. In particular, we have measured a correlation coefficient between the evolution of the lift and the flap position equal to 0.6. These findings are quite similar to the ones observed for the 2D laminar flow where the flaplet was regularising the lift/drag cycle with a movement characterised by the same value of the lift-flap position correlation coefficient.
To shed some more light on the mechanism driving the flap motion and the corresponding lift increase, in Fig. 24 we consider the instantaneous streamwise velocity profiles (displayed in the top panel) and the pressure coefficients (shown in the bottom one), sampled at two time instants corresponding to an ascending and descending flap movement. Note that all the profiles have been obtained after having averaged in the spanwise, homogeneous direction. When the flap is moving downward (blue line), the large recirculation region (negative velocity) is enclosed between the outer flow at the top, and a region characterised by positive velocities at the bottom. The fluid trapped in this region is displaced downstream under the action of the low pressure values associated with the core of the trailing edge vortex. This observation is confirmed by the pressure coefficient recorded at the same time instant showing a strong negative value at the trailing edge consistently with the   To provide a phenomenological explanation on the increased regularity of the shedding cycle, we have computed conditional averages of the flow fields. In particular, the averages have been conditioned by the value of the lift coefficient (i.e., ensemble averages between samples sharing the same phase in the shedding cycle). In particular, we averaged spanwise vorticity fields corresponding either to the maximum Fig. 25a and b or to the minimum displaced to the right when the flaplet is used. Moreover, in the case with the flap, the lift generating vortex in the maximum lift condition shows higher values of conditional averages of spanwise vorticity, as shown by a dense and compact region of saturated blue colour, thus indicating an increased level of coherence. Concerning the wake, the vortex street generated with the flap shows an almost uniform sequencing of the counter rotating vortices sharing the same spatial locations. The enhanced regularity of the cycle is also confirmed in Fig. 26a showing the time cross correlation ρ of the lift coefficient c l and the spanwise vorticity ω z at location (2.0c; 0.4c) (x coordinate measured from the leading edge, y from the profile chord). The cross correlation ρ is defined as follows where E [ ] and σ [ ] indicate the expected value and the standard deviation, respectively. In the case with flap, the evolution of the time cross correlation shows a clear periodic behaviour with high levels of correlations (0.35). Whilst in the case without flap, the correlation is much lower (0.05). Finally, the right panel of the figure shows the time cross correlations of the lift and drag coefficients with the flap movement defined by its elevation y. As already said, both the aerodynamic force coefficients are strongly correlated with the flap movement, especially the lift coefficient which has a value of cross correlation almost double the one for the drag coefficient. A phase shift is also apparent for the drag coefficient, whose peaks slightly precede the one of the lift coefficient. To summarise the last results, all this cross-correlation graphs show that the aerodynamic coefficients are strongly linked  As already done in the 2D case, to determine which mechanism is the main responsible for the increase in average lift obtained with the flap, we have computed the circulation over two closed surfaces C embedding the lift and the trailing edge vortices, respectively. The former is defined over the region x ∈ [0.5, 1.0], y ∈ [0.3, 0.6], the latter covers the area x ∈ [0.8, 1.3], y ∈ [0.0, 0.3] (see Fig. 7). Similarly to what we have observed for the 2D laminar case, the circulation of the leading edge vortex (the one that generates lift) is only slightly increased by the flap presence (≈ 2%), whilst the circulation of the trailing edge vortex (the one that reduces the lift, or increase the downforce) is substantially reduced by a factor of ≈ 15%.

Conclusion
This numerical study focused on the use of passive, self actuated flaps as lift enhancement devices in nominally stalled conditions. The main objective was to discover how the mutual interaction between these self deployable devices and the unsteady flow field generated by a foil at high angle of attack can improve the aerodynamic efficiency of stalled wings. Although the design of optimal flaps (i.e., delivering maximum lift increase) was not a primary objective of this work, we had to carry out a preliminary selection study to determine the characteristics (i.e., size, location and natural frequency) of a self-adaptive flaplet able to deliver substantial aerodynamic benefits in an otherwise stalled condition. This initial study has been conducted on a baseline NACA0020 aerofoil at 20 o degrees angle of attack at low (fully laminar) chord Reynolds number (i.e., Re c = 2 × 10 3 ). The impact on the aerodynamic performance of a rigid, thin flap hinged with a torsional spring on the aerofoil suction side has been analysed via a parametric study involving the size of the flap, the hinge location and the spring stiffness. It has been found that it is of fundamental importance to lock-in the flap oscillation frequency with the foil Strouhal number. In resonating conditions, the lift response become quite sensitive to the geometric properties of the flap. In particular, the quasi optimal performances (i.e., ≈ 20% increase in lift) are achieved with a flap length of one fifth of the chord hinged at about 70% of the aerofoil. Having determined the geometric and physical character of an aerodynamically efficient flaplet, we turned our attention to the understanding of the mechanisms responsible for the improved foil performances at high angle of attack. To this end, we have carried out Direct Numerical Simulations of the flow past a NACA0020 aerofoil at 20 o angle of attack at a chord Reynolds number of 2×10 4 considering both the baseline wing and the wing equipped with the quasi optimal flaplet determined in the 2D parametric campaign. Initially, considering the baseline wing, we have confirmed that the flow mechanisms taking place in the fully three-dimensional scenario, involving a laminar separation, a subsequent reattachment and a laminar-turbulent transition, determine a flow behaviour that is qualitatively similar to the two dimensional case used for the preliminary design study. The reasons for this similarity are related with the common laminar separation, and the convective inviscid instability of the leading edge shear layer responsible for the roll up of the large recirculation bubble on the aerofoil. In a second phase, we have systematically compared the flow fields generated with and without the flap. Although the mean velocity fields and the mean kinetic energy are very similar, the flaplet has a very strong impact in manipulating the unsteady character of the vorticity field. In particular, the flap is popped up by the passage of the lift vortex and when relaxing back to the equilibrium position generates a jet almost tangent to the wing surface, directed towards the trailing edge. This jet detaches the vortex street generated by the trailing edge shear layer instability away from the aerofoil. The displacement of the trailing edge vortices has a twofold effect. On one hand there is a net decrease in the downforce that is directly generated by these vortices leading to a global increase of the lift. On the other hand, the displacement of the trailing edge vortex allow for a complete evolution of the leading edge generated vortex that now does not interact directly with the trailing edge vortices. As a consequence, the periodic character of the wake is now mainly controlled by the shedding of the leading edge vorticity into the wake that regularises the shedding cycle also promoting a much more ordered wake topology.
Funding Information This study was funded by the European Commission, 7th framework project Pelskin (grant number 334954) and by the EPSRC, project Quiet Aerofoils of the Next Generation (grant number EP/N020413).

Compliance with Ethical Standards
Conflict of interests The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.