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

Abstract

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 20^o 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 × 10⁴ 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.

1 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⁶ (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⁴) 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⁵ 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⁶. 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⁵) in ramp-up motion (α₀ = 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, self-adaptive, 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⁴. 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.

2 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₁ and x₂) denoting the directions parallel and normal to the aerofoil chord, and z (x₃) being the axis normal to the paper. Also, u, v and w (u₁, u₂ and u₃) 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:

$$\begin{array}{@{}rcl@{}} &&\frac{\partial u_{i}}{\partial t} +\frac{\partial u_{i} u_{j}}{\partial x_{j}} = -\frac{\partial P}{\partial x_{i}} + \frac{1}{Re_{c}} \frac{\partial^{2} u_{i}}{\partial x_{j} \partial x_{j}} + f_{i}, \end{array} $$

(1)

$$\begin{array}{@{}rcl@{}} &&\frac{\partial{u_{i}}}{\partial x_{i}} = 0, \end{array} $$

(2)

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.

Fig. 1

a Sketch of the computational domain. b Grid in the proximity of the aerofoil (nodes are plotted with a skip index of six). The inserted figure is an enlargement of the area surrounding the trailing edge

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₁, x₂ and x₃ 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].

2.1 Fluid-flap interaction model

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:

$$ I \ddot{\theta} + C \dot{\theta} + K \theta = \mathcal{T}, $$

(3)

In Eq. 3, I is the flap moment of inertia with respect to the rotation axis (i.e., I = mL²/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, $\phantom {\dot {i}\!}\mathcal {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: K = (2πf)²I.

Fig. 2

Sketch of the flap hinged on the suction side of the aerofoil

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 Γ: $\phantom {\dot {i}\!}\boldsymbol {U}(\boldsymbol {X}_{i})=\mathcal {I}(\boldsymbol {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 Γ:

$$ \boldsymbol{F}^{*}(\boldsymbol{X}_{i})=\frac{\boldsymbol{U}^{{\Gamma}}(\boldsymbol{X}_{i}) - \boldsymbol{U}^{*}(\boldsymbol{X}_{i})}{{\Delta} t}. $$

(4)

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 $\phantom {\dot {i}\!}\mathcal {I}$, indicated with $\phantom {\dot {i}\!}\mathcal {C}$ and termed as spread. The spread operation formally allows to determine the singular forces on the fluid grid as:

$$ \boldsymbol{f}^{*}(\boldsymbol{x}_{i,j,k})=\mathcal{C} \left( \boldsymbol{F}^{*}(\boldsymbol{X}_{i}) \right). $$

(5)

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.

Fig. 3

Schematic of a full time step for the coupled fluid-wing-flap system

The distinguishing feature of any continuous forcing IB method is related to the way in which the two operators $\phantom {\dot {i}\!}\mathcal {I}$ and $\phantom {\dot {i}\!}\mathcal {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:

$$ f_{a} \left( x \right) \approx {\int}_{{\Omega}} w_{d} \left( x-s \right) f \left( s \right) ds, $$

(6)

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:

$$ {\int}_{{\Omega}} (x-s)^{i} \, (y-t)^{j} \, (z-v)^{k} \; p(x-s,y-t,z-v) w(x-s,y-t,z-v) ds dt dv=\left\{\begin{array}{ll} 1 & \text{if}\; \; i=j=k = 0 \\ 0 & \text{otherwise} \end{array} \right. $$

(7)

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].

3 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⁴. 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 convective Kelvin-Helmholtz (KH) instability that determines its downstream development eventually leading to turbulent transition. A similar behaviour is observed for the trailing edge shear layer that undergoes a KH instability too with a consequent roll up responsible for the formation of the trailing edge vortex street Fig. 5. Further downstream, past the aerofoil, a large wake is formed by the joint contribution of the vorticity generated from both the leading and trailing edges. The uneven vorticity contributions from the two layers is ultimately responsible for the lack of symmetry characterising the wake topology. The global effect of the wake unsteadiness can be evinced from Fig. 4b showing the time evolutions of the lift and drag coefficients obtained by integrating the wall pressure and the shear stress at each time step (mean values: c _l = 0.64 and c _d = 0.35). From the figure, one can observe the presence of a dominant oscillation period clearly associated to the alternating vortex shedding in the wake, with a corresponding non dimensional frequency, in terms of Strouhal number, equal to St = f _s c/U _∞ ≈ 0.534 [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 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⁴ 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³. 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.

Fig. 4

a Contours of mean flow streamwise velocity and streamlines. Contour goes from $-0.3U_{\infty }$ (blue) to $1.4U_{\infty }$ (red). b Instantaneous lift c _l (solid line) and drag c _d (dashed line) coefficients as a function of time

Fig. 5

Contours of instantaneous flow streamwise velocity and streamlines. Contour goes from $-0.3U_{\infty }$ (blue) to $1.4U_{\infty }$ (red), and the snapshots cover a full shedding period of $1.87c/U_{\infty }$

Fig. 6

Contours of the mean flow x-component velocity u. The colour contours are used for the 2D case, and goes from $-0.4U_{\infty }$ (blue) to $1.2U_{\infty }$ (red), whilst the contour lines (with the same levels separated by $0.6U_{\infty }$) are used 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.

Fig. 7

Contours of the instantaneous spanwise component of vorticity ω _z , corresponding to the minima (a, b) and maxima (c, d) locations of the lift coefficient over one shedding period for the 2D (a, c) and 3D (b, d) cases. Blue lines used for negative, clockwise vorticity, red ones for positive values ($\pm 5U_{\infty }/c$)

4 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³), 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 corresponds to the shedding of a vortex, at a shedding frequency equal to f _s = 0.555U _∞ /c. The lift coefficient evolution also shows the presence of a lower frequency f = 0.308U _∞ /c (almost half the shedding frequency). The instantaneous vorticity fields ω _z over this two shedding periods are shown in Figs. 9 and 10 (left column). The presence of two dominant vortices formed as a consequence of the leading and trailing edge shear layer instabilities characterises all the time series. In particular, their opposite circulations are responsible for the lift and downforce generated by the clockwise rotating vortex (blue), and the counter clockwise rotating one (red), respectively. The first few snapshots of the reported vorticity time series correspond to a maximum lift condition in which the leading edge vortex has already formed whilst the trailing edge one is rolling up, on the verge of being shed from the aerofoil 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.

Fig. 8

Instantaneous lift c _l (solid line) and drag c _d (dashed line) coefficients as a function of time. The dots indicate the selected time snapshots shown in Figs. 9 and 10

Fig. 9

Contours of the instantaneous spanwise component of vorticity ω _z during two shedding cycles (corresponding to $3.247c/U_{\infty }$ non-dimensional time units) for the baseline (left column) and optimal flaplet configuration (right column). The snapshots correspond with the time instants marked in Fig. 8. Blue negative (clockwise) vorticity, red positive (counter clockwise) in the range $\pm 5U_{\infty }/c$

Fig. 10

Continuation of Fig. 9

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₀ 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 changes as a function of the hinge position when considering a L = 0.2c long flaplet in resonating conditions. The optimal position, in terms of maximum lift, is at about 0.7c, where, when unlifted, the flaplet end almost reaches the trailing edge (Figs. 11 and 12).

Table 1 Flap configurations considered in the 2D parametric study

Fig. 11

Optimal flaplet configuration: instantaneous lift c _l (solid line) and drag c _d (dashed line) coefficients as a function of time. The thin solid line represents the elevation y of the tip of the flap. The set of bullets on the graphs indicates the instants in time where the vorticity snapshots have been sampled, see Figs. 9 and 10

Fig. 12

Mean lift c _l as a function of a the spring natural frequency f and of b the position of the hinge x _F . The dashed, solid and dash-dot lines are used for the cases with L = 0.1c, L = 0.2c and L = 0.3c, respectively

In summary, when a low Reynolds number, 2D case at an angle of attack of α = 20^o is considered, the flaplet configuration that maximises the mean lift features a length of 0.2c, a 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). 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 above the aerofoil inducing a maximum lift force. The cycle is closed by the formation of a new trailing edge vortex Figs. 9j and 10b. The mutual interaction of the flow field with the flaplet has a strong impact on the shedding process and therefore on the structure of the wake (see Figs. 10b and j). The importance of the interaction is further stressed by the high correlation between the lift oscillations and the flap motion (correlation coefficient is ≈ 0.6), and in particular by the fact that the maximum lift is reached when the flap is almost at its maximum elevation (the time lag between the two functions is ≈ 0.2c/U _∞ ).

As a further measure of the effect of the flaplet on the vorticity field, we have quantified the circulations of the velocity field along two closed rectangular loops bounding the lifting vortex (x ∈ [0.5,1.0], y ∈ [0.3,0.6]), and the trailing edge vortex (x ∈ [0.8,1.3], y ∈ [0.0,0.3]), respectively (see 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.

5 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⁴, 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. 13

a Lift c _l (black) and drag c _d (grey) coefficients as a function of time. b Evolution of the aerodynamic efficiency E = c _l /c _d . Solid lines are used for the aerofoil with flap, whilst dashed lines for the reference values

Figure 13a compares the time evolution of the lift and drag coefficients of the reference case versus the ones obtained when using the flaplet. Their time averaged values are c _l = 0.64 and c _d = 0.35, for the baseline case, increasing to c _l = 0.74 and c _d = 0.37 with the flap, thus obtaining a 16% improvement in lift and a slightly augmented drag (6%). Also, the r.m.s. of the lift coefficient increases from 0.15 to 0.17 (14%). The aerodynamic efficiency, E = c _l /c _d is reported in 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.

Fig. 14

a Pressure C _P and b friction C _f coefficient distributions. Solid and dashed lines are used for the aerofoil with and without flap, respectively

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.

Fig. 15

Contours of the mean (i.e., time and z-averaged) spanwise component of vorticity ω _z and mean streamlines of the NACA 0020 aerofoil at α = 20^o and Re _c = 2 × 10⁴. Left panel a results for the baseline wing; right panel b wing equipped with a flaplet (L = 0.2c, x _F = 0.7c, K = 0.150). Blue negative vorticity (clockwise), red positive ($\pm 7U_{\infty }/c$)

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 of the flap, at x = 0.9c, in the region traversed by the flap oscillations, the velocity intensity is reduced because of the no-slip and no-penetration boundary condition on the flap solid surface. Finally, in the near wake region, the velocity defect is slightly enhanced in the case with flap.

Fig. 16

Mean x-wise (a) and y-wise (b) velocity components profiles over the aerofoil and in the near wake. Lines are used for the aerofoil without flap; symbols refer to aerofoil with flap

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 $\phantom {\dot {i}\!}k = 1/2 <u^{\prime }_{i} u^{\prime }_{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.

Fig. 17

a Mean turbulent kinetic energy profiles over the aerofoil and in the near wake, and b further downstream at x ≈ 5.5c. Lines are used for the aerofoil without flap; symbols refer to aerofoil with flap

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 $\phantom {\dot {i}\!}\mathcal {P}\left (u<0\right )$ 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.

Fig. 18

Intermittency factor $\mathcal {I}_{-}=\mathcal {P}\left (u<0\right )$ for the reference (a) and flap (b) cases. Contour levels go from blue ($\mathcal {I}_{-}= 0$) to red ($\mathcal {I}_{-}= 1$)

To gain further insight on the effect of the flaplet-flow interaction we have used the classical $\phantom {\dot {i}\!}\mathcal {Q}$-criterion proposed by Hunt et al. [36]. This technique assigns a vortex to all spatial regions that verify the condition

$$ \mathcal{Q}=\frac{1}{2} \left( \vert \boldsymbol{{\Omega}} \vert^{2} - \vert \boldsymbol{S} \vert^{2} \right) > 0, $$

(8)

where $\phantom {\dot {i}\!}\boldsymbol {S} = \frac {1}{2} \left ( \boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla } \boldsymbol {u}^{T} \right )$ is the strain rate tensor and $\phantom {\dot {i}\!}\boldsymbol {{\Omega }} = \frac {1}{2} \left ( \boldsymbol {\nabla } \boldsymbol {u} - \boldsymbol {\nabla } \boldsymbol {u}^{T} \right )$ is the vorticity tensor. Instantaneous $\phantom {\dot {i}\!}\mathcal {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]: the vortices rolling up under the action of the KH leading edge shear layer instability and the street of vortices shedding from the trailing edge. The main features of this flow process, largely present also in the 2D, laminar case, remain practically unaffected by the presence of the flap, except in the leading edge area, where the KH instability is delayed further downstream Fig. 20. The regularisation of the shedding of vortices from the leading edge is responsible of the increased pressure coefficient Fig. 14a, that ultimately produce the increased lift coefficient.

Fig. 19

Visualisation of instantaneous vorticity field by means of $\mathcal {Q}$-iso-surfaces ($\mathcal {Q}= 450U_{\infty }^{2}/c^{2}$) coloured by the spanwise vorticity. a and b are the cases without and with the flap, respectively

Fig. 20

Visualisation of instantaneous vorticity field by means of $\mathcal {Q}$-iso-surfaces ($\mathcal {Q}= 450U_{\infty }^{2}/c^{2}$) coloured by the spanwise vorticity. a and b are the cases without and with the flap, respectively

Fig. 21

Istantaneous contour plot of the FTLE σ^T during a shedding period for the case with flap. The contour levels go from 0 (white) to $7U_{\infty }/c$ (red). The black contour lines are used for the baseline case without flap

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

$$ \sigma^{T} \left( \boldsymbol{x}, t \right) = \frac{1}{T} ln \sqrt{\lambda_{max} \left( {\Delta} \right)}. $$

(9)

Here, λ _{m a x} (Δ) is the largest singular value of the Cauchy-Green deformation tensor computed over a finite time interval [t₀, t₀ + T]

$$ {\Delta} = \frac{\partial \boldsymbol{x} \left( t_{0} + T, \boldsymbol{x}_{0}, t_{0} \right)}{\partial \boldsymbol{x}_{0}}. $$

(10)

Figure 21a is the beginning of the shedding cycle from the trailing edge, with no vortex at the trailing edge, whilst at the leading edge shear layer rolls up under the action of a KH instability. Figure 21 show how the trailing edge shear layer undergoes a KH instability and a vortex is generated. The detachment of the vortex is illustrated in Fig. 21. In the two cases, the instantaneous shapes of the leading edge shear layers are very similar. However, the locations of the trailing edge vortex cores at these two time instants (panels b and d of Fig. 21) do not correspond with the controlled case showing a streamwise shift of the actual position. This relative displacement is induced by the downward movement of the flap and the consequent appeareance of a positive streamwise velocity generated by the relative movement of the flap and the aerofoil which will be analysed later on. The opposite effect is also visible during the flap lifting phase (panels c and d), when fluid momentum is entrained upstream by the flap movement.

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.

Fig. 22

Baseline aerofoil (left column) and aerofoil equipped with the flaplet (right column): contours of the instantaneous spanwise component of vorticity ω _z , over two shedding periods. Blue negative vorticity (i.e., clockwise), red positive ($\pm 5U_{\infty }/c$)

Fig. 23

Continuation of figure 22.

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 incipient formation of the trailing edge vortex (see Fig. 23h) and a decrease of the torque $\phantom {\dot {i}\!}\mathcal {T}$ acting on the flap. On the other hand, during the upward motion of the flap (red line) the sign of the recorded streamwise velocity in the region between the flap and aerofoil is negative. In this condition, the full detachment of the trailing edge vortex Fig. 23d reduces the suction also allowing for an increase of the torque on the flap.

Fig. 24

(Top) Spanwise averaged streamwise velocity profile u as a function of the y-coordinate at two time instants. Three velocity profiles are provided: at x = 0.85c, 0.9c and 0.95c. Note that, the profiles have been shifted for clarity indicating with the vertical dashed line the respectives zeros. (Bottom) Spanwise averaged pressure coefficient C _p recorded at the same time instants as above. In particular, in both panels, the blue and red lines correspond to the downward and upwards rotation of the flap, Figs. 23d and h, respectively

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 Fig. 25c and d lift force for both the cases. For both situations of minimum and maximum lift, it is possible to notice that the positive rollers (red ones, generating downforce) are 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

$$ \rho \left( \tau \right) = \frac{E\left[ c_{l} \left( t \right) \omega_{z} \left( t+\tau \right) \right]}{\sigma \left[ c_{l} \right] \; \sigma\left[ \omega_{z} \right]}, $$

(11)

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 with the flap movement and with the vorticity field, proving that the flap movement is determined and linked with the vortex dynamic which ultimately determines the aerodynamic behaviour.

Fig. 25

Contours of the conditional averaged spanwise component of vorticity ω _z , for the case with (b, d) and without (a, c) flap. Blue negative vorticity (i.e., clockwise), red positive ($\pm 5U_{\infty }/c$). The top and bottom rows correspond to the times of maximum and minimum lift, respectively

Fig. 26

a Time cross correlation ρ of the lift coefficient c _l and spanwise vorticity ω _z at the location $\left ( 2.0c; 0.4c \right )$. Solid and dashed lines are used for the case with and without flap, respectively. b Time cross correlation ρ of the lift coefficient c _l (solid line) and drag coefficient c _d (dashed line) with the flap elevation

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%.

6 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³). 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⁴ 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.