A Ayvazian

I Cuenca

E R Gowree

FP92-AERO2022-Gowree On bio-inspired induced drag reduction techniques

The current study was aimed at identifying wing planform geometries which favours low induced drag. The geometries were defined by mathematical functions and following a subsequent parametric study, with a focus on sweep angle and taper ratio, inspired by birds wings and aquatic animals flukes or caudal fins planforms. The parametric studies were carried out using a vortex lattice method (VLM) which was validated with experimental data from literature. The main finding is a confirmation of the fact that the maximum Oswald efficiency factor, ε is not exclusive to an elliptical wing planform. Despite its apparent simplicity, in fact a swept-back-tapered wing with a linear chord distribution can reach ε ≈ 1, the limit or minimum induced drag established from elliptical lift distribution. Moreover, the crescent wing geometry, such as the one displayed by the swift, can even surpass this limit.

INTRODUCTION

The implementation of new techniques for reducing drag in commercial aircraft will help in reducing fuel consumption and thus emission from aircraft which has detrimental impact on the environment. In the search for new technologies to reduce drag, thanks to its millions of years of evolution based on optimal energy consumption, nature is an exceptional source of inspiration for potential solutions. By taking inspiration from nature here we attempt to identify and propose potential planform geometries for minimum induced drag, seen in migratory birds which display remarkable endurance. The major parameters that defines the planform geometry are the taper and sweep when viewed from the top. Most of the man made systems will consist of a simplified uniform or linear distribution of taper or sweep angle, however when compared to natural flyers, in Fig. 1, the taper, λ and sweep angle Λ are seen to vary in a more non-linear fashion for instance in the crescent-moon shape of the swift's wing, where the from first sight the leading edge can be expected to follow an elliptical profile. In case of the falcon the sweep even changes direction from forward in the inner region to backward in the outer region constituting of the primary feathers, similar variation are also observed on frigates not shown here and is slightly less pronounced on albatross and this configuration will be referred as the "M-configuration" here onwards. In Fig. 1, the corresponding optimised platform geometries emanating from the current study has been presented and interesting point to consider here is that λ → 0 in most natural flyers.

The lifting line theory provides a simple model for the induced drag, C D i , for an arbitrary spanwise lift distribution along the span and it can be expressed as

C D i = C 2 L πεAR (1)
C D i is therefore directly proportional to the lift coefficient, C 2 L , and inversely proportional to the aspect ratio AR, and the Oswald efficiency factor, ε which describes how efficiently the lift is distributed along the span, and for an elliptical lift distribution it is equal to 1, thus the minimum induced drag. On a finite wing the spanwise lift distribution is compromised by the spillage at the tip, leading to the formation of tip-vortices which produces the downwash effect leading to an induced drag component and a non-uniform spanwise distribution of circulation. The efficiency of lift generation and the induced drag can be reduced by re-organising the spanwise distribution in lift, through the introduction of taper, twist, sweep, dihedral and local change in aerofoil geometry to quote a few parameters. As a consequence, these parameters were targeted throughout the current study. However, it was demonstrated by Zimmer [START_REF] Zimmer | The significance of wing end configuration in airfoil design for civil aviation aircraft[END_REF] and later by van Dam [START_REF] Van Dam | Induced-drag characteristics of crescent-moon-shaped wings[END_REF] that crescent-moon planform shapes inspired by nature can significantly reduce the induced drag with the Oswald efficiency factor exceed 1, which violates the limit proposed by lifting line theory for elliptical lift distribution.

Amongst the first, van Dam [START_REF] Van Dam | Efficiency characteristics of crescent-shaped wings and caudal fins[END_REF] investigated the benefits of the crescent wing geometry inspired by nature. The efficiency of an elliptically-loaded crescent shaped wing was compared with those of an elliptical wing with the same spanwise chord distribution.The crescent wing showed an improvement in the Oswald efficiency factor of nearly 9%. He concluded that the aerodynamic efficiency is improved when the backward curvature was increased. By combining this with a large wingspan, it is possible to greatly reduce induced drag. In this research it was demonstrated also that this type of geometries generated non-planar trailing vortex sheets that produce this important improvement in terms of induced drag. Later BURKETT [START_REF] Burkett | Reductions in induced drag by the use of aft swept wing tips[END_REF] studied the benefits of aft swept wing tips, similar to that on the swift and the albatross wing. The rearward curvature on these type of wings can reduce the induced drag factor to values below those obtained from unswept elliptical wing once again considered the optimal geometry in classical finite wing theory. These studies showed the potential of crescent wing, like those of swifts, has the possibility to offer significant reductions in induced drag specially at moderate angles of attack. These crescent wings were compared with elliptic geometries during wind tunnel test campaigns at low Reynolds numbers again by van Dam and co-workers [START_REF] Van Dam | Experimental investigation on the effect of crescent planform on lift and drag[END_REF] [2] confirming the significant reduction in the induced drag, but here the leading-edge separation-induced vortex flow over the highly swept and tapered tips of this type of wings was postulated to be responsible for the benefits.

The current study was intended to revisit the above mentioned induced drag reduction technique and provide further insight into the mechanism which can be exploited for unmanned aerial systems (UAS) and (MAS) applications which operated in similar Reynolds number regime as natural flyers. A low order vortex lattice method (VLM) was first employed to generate the optimal planform geometry tending towards the natural flyers and higher order CFD was then employed to analyse the behaviour of the flow.

METHODS

For the geometrical optimisation process the XFLR5 software was used due to its low computational cost, which permits a the simulation and analysis a many cases, of the order of hundreds in this case due to the low order modelling approach. The geometries where obtained by varying the local sweep and taper and this was either done using linear or non-linear functions defined in section 2.1. Once an optimal geometry was identified CFD simula-tions were conducted using a commercial package to further understand the mechanism of induced drag reduction through the re-organisation of the spanwise flow field.

Wing planform geometry generation

A major step involved the planform geometry definition based on the bird species. Despite an extensive review of the literature, it was still difficult to obtain the real wing geometry, so we proceeded by analysing images of birds in flight and reconstructed the most representative wing planforms using mathematical functions described below, this allowed for the local change in sweep angle and taper ratio. A series of planform were thus generated using the mathematical functions and was analysed using XFLR5. From equation 1, the induced drag is a strong function of the aspect ratio, AR and C L , in order to compare the efficiency of different planforms in reducing C D i , these two parameters had to be kept constant.

Crescent wing

The mathematical functions to define the crescent wing is defined below and also plotted in Fig. 2.2.2 and 2.1.1. Since the elliptical wing planform is known to minimise induced drag, two ellipses were used to define the leading and trailing edges. This helped in defining wing planform similar to that studied experimentally by van Dam [START_REF] Van Dam | Experimental investigation on the effect of crescent planform on lift and drag[END_REF] shown in Fig. 2.2.2 which also shows the equations and the geometrical parameters used to define this geometry.

Mathematically, an ellipse is defined by the semi-axes, which means that here there are 4 unknowns. In addition to these, an offset for one ellipse has been introduced, in order to have a non-zero wing area in case the leading and trailing edges are the same, like in the case where, λ = 1 is desired, resulting into 5 unknowns. However, the major axes are considered to be pre-defined parameters, because the wingspan, b, is kept constant throughout the entire study, thus reducing the system of equations to 3 unknowns. Since the, taper ratio, λ is known to affect the spanwise loading and thus induced drag it was kept as a variable. To achieve a crescent wing planform sweep, α, was introduced which connected the leading edge at the tip and the root would also impact the planform shape of the wing and therefore the induced drag. The parameters y 1 and y 2 are not related to the wing planform, they instead represent the semi-minor axes of the ellipses. h represents the chord length at the tip and is directly linked to the offset of the equation representing the leading edge. After some mathematical manipulations, the resulting constraints are expressed in Eqs. 2, 3 and 4 below.

y 1 = 1 2 bsin(α) (2)
Airbus A320

Swift Peregrine Falcon stooping Albatross

Table 1: Illustration of some wing planforms observed in nature and from man-made systems. Images from open access sources.

Figure 1: The leading and trailing edge functions to define the geometry of the elliptical crescent wing on the left and the parabolic crescent wing on the right

h = b AR(1 -(λ -1) π 4 ) (3)
y 2 = y1 + (λ -1)h (4)
In this study we also studied crescent wings defined by parabolic functions, which can be defined by the functions shown in Figure 2.1.1. Again, there are three constants that need to be determined in order to define the wing planform: a, c and d. The parameters chosen to act as constraints are the same as in the case of the crescent wing defined with elliptical functions: λ , AR and α. The wingspan, b, is kept constant throughout all the studies, therefore it was treated as a known parameter. Following some mathematical manipulations, the constraints imposed could be expressed, equations Eqs. 5, 6 and 7.

d = 0.5b tan(α) (5)
a = ARφ (d 2 ) 4d 2 + bARφ (6) c = ( 1 a 2 - 1 d 2 ) b 2 4(λ -1) (7)
where φ is a parameter introduced for mathematical convenience:

φ = 1 3 + 1 λ -1 (8)

M-configuration wing

The albatross is also a high endurance flyer [START_REF] Stempeck | Aerodynamic performance of albatross-inspired wing shape for marine unmanned air vehicles[END_REF], but their wing planform as shown on the top of Fig.

c root = 2b AR(1 + λ ) (9)
c tip = λ c root (10)
c slope = 2(c tip -c root ) b (11)

Numerical methods

CFD simulations were also conducted, for seeing the behaviour of the flow around the shape of crescent wings. For validating the CFD methods used, the wing geometry of [START_REF] Van Dam | Experimental investigation on the effect of crescent planform on lift and drag[END_REF] experimental studies, Figure ??, was simulated with STAR-CCM+. The models used for the simulation were Kω turbulence model with the submodel turbulence supression for simulating the boundary layer transition at the 5% of the chord. The computational domain used was defined by 20 chords downstream, 5 chords upstream and 7 chords for top and bot margins. A mesh independence study was carried on resulting finally in a mesh of around 5.300.000 cells.

Low order modelling 2.2.2 High order modelling

Following the optimisation using XFLR5, more detailed RANS simulations were undertaken to generate further insight into the behaviour of the flow, limited in XFLR5, we would like to emphasis on the fact this approach was not used to calculate the induced drag directly. Here, the commercial package STAR-CCM+ was employed. The computational domain was defined by 20 chords downstream, 5 chords upstream and 7 chords for top and the bottom of the wing placed at the centre. In order to resolve the flow in the wake downstream a cubic region of very fine mesh was introduced as shown in Fig. at least part of the evolution of the vortical flow downstream. For very accurate calculations of the induced drag from the quantities in the wake, such a mesh needs to be extended until the Treftz plane which can be assumed to start at approximately 50 chords downstream at least which would result in a very large number of cell. Calculations with such as number of cell was not feasible within the scope of this project, but as shown in section 2.3, the RANS simulation was reliable enough to provide the spanwise loading and the total lift which could in agreement with experimental data of van Dam and co-workers [START_REF] Van Dam | Experimental investigation on the effect of crescent planform on lift and drag[END_REF]. The kω turbulence model with a turbulence supression submodel for simulating the boundary layer transition at 5% chord was employed so as to simulate the tripping position during the experiment. A mesh independence study was carried-out on the unstructured mesh resulting in a total number of mesh of around 5.3 million cell showing the best compromise.

Figure 3: The mesh refinement in the wake of the wing

Validation of numerical methods

Although, XFLR5 has shown to be very reliable for 2D aerofoils over a large range of Reynolds number, the 3D calculation over finite wings has not received a thorough validation in the literature. In addition, due to the highly dissipative characteristics of unstructured RANS simulation which also handles the flow in the far wake rather poorly, due to the expansion in the cell size downstream for a better compromise in terms of simulation time and cost. Here, we proceeded by a verification of the numerical method with the experimental results of van Dam and co-workers [START_REF] Van Dam | Experimental investigation on the effect of crescent planform on lift and drag[END_REF] in the subsonic wind tunnel at NASA Langley. The experiment was conducted over a crescent wing shown on the top left of Fig. 4, with an AR = 6.9 at a Reynolds number Re = 4.8 × 10 5 based on the root chord and Mach number, M = 0.205 for a range of angle of incidence. A similar geometry and flow conditions were simulated in XFLR5 using both the LLT and VLM methods with transition fixed at 5% chord during the 2D batch calculations. In RANS formulation the Kω turbulence model and the turbulence surpression sub-model was employed for simulating the boundary layer transition at 5% of chord as in the experiment.

From the total lift coefficient, C L shown on the top right of Fig. 2.3, an expected linear variation due to the Reynolds number being sufficiently high and in the presence of a boundary layer trip is observed. While the results from LLT methods showed noticeable overprediction with respect to the experimental results, the VLM methods showed closer agreement. No, noticeable difference was found between the number of panels and therefore the rest of the optimisation process using XFLR later was undertaken for 2000 panels. The RANS results also showed acceptable agreement while being slightly under-predicted. The total drag coefficient, C D showed the usual U-shape behaviour where drag increases with increasing angle of attack. At zero lift the total drag constitute purely of the profile drag which can be split into skin friction drag and form drag. But, with increasing angle of incidence the contribution of form drag and now induced drag starts to increase and cannot be easily separated from the total drag. Once again the VLM method showed reasonable agreement with respect to the experimental results, while the LLT was significantly under-predicting the drag. Except at zero or very low angle of attack or lift condition, the results from the RANS simulation was significantly over-predicted. This is due to the inability of RANS to account for the increase in form drag, but more importantly the induced drag as the quantifies in the far wake suffers from dissipation potentially due to the meshing strategy, resulting in large errors in calculating the total drag, from integration of the far wake quantities.

When comparing the spanwise lift distribution on the bottom right of figure Fig. 2.3 and elliptical C L distribution is observed, where the agreement between the experiment, VLM method and the RANS simulations are satisfactory up to about a 80% along the span. Further outboard the numerical simulation results showed slightly better agreement as opposed to the VLM method. From this comparison, it is fair to suggest that he RANS simulation can predict the near-field with reasonable accuracy and could be therefore employed to study the spanwise re-arrangement in the flow field over the wing to generate further insight into the potential induced drag reduction mechanism, despite the overall induced drag prediction being unreliable due to the meshing strategy. Since, the current validation was undertaken over a crescent wing the analysis that follows could be carried-out with reasonable fidelity.

PARAMETRIC STUDIES

For all the parametric optimisation studies, in order to have comparable geometries, the Aspect Ratio was kept constant at 17.664. This is in the range of the average values of the adults albatross [START_REF] Richardson | How do albatrosses fly around the world without flapping their wings?[END_REF], and they were chosen because they are representative of animals and UAV geometries. The airspeed was selected to 20 m/s, which is the normal cruise speed of the albatross corresponding to an average Re ≈ 250 × 10 3

Tapered Wing

From the literature the two influential geometric parameters for induced drag were identified as the taper taper ratio and the sweep. In conventional subsonic, transonic or supersonic wing design, the later is usually favoured for the reduction of wave drag or increase in drag divergence Mach number and improvement of stability but little consideration has been given for induced drag. Here, we revisited the effect of uniform sweep and taper on induced drag starting from a rectangular planform. Following the parametric optimisation process, from the top Fig. 5, even for a simple uniform sweep and taper the Oswald efficiency factor can be increased to a value of 1 (shown by the red dot), which corresponds to the minimum achievable from lifting line theory based on an elliptic lift or circulation distribution. This is achieved for the case of λ = 0.25 and sweep angle of 25 • . The current finding has not been reported previously on uniform swept-tapered wing, however most of the previous studies have been focussed on high Reynolds number and taper ratio. The effect of low freestream Reynolds number which reduces even further the local Reynolds number towards the tip can affect the local lift and this has not been exploited to great length.

M-configuration

Following on with uniform sweep and taper configurations, the effect of changing the direction of the sweep angle in the in-board section of the wing can be implemented through the M-configuration defined above. This configuration is expected to approach the planform geometries of albatross, also known for very high endurance. A similar optimisation process to that in the section above was adopted and in this case 5 parameters can be varied, the back or forward sweep angles, α 1 , α 2 and α 3 , and two breakpoint position along the wingspan where these angles were varied, S 1 and S 2 . In order, to Following the optimisation process, it was rather surprising to see that the M-configuration which was qualitatively representative of the albatross wing did not yield a high Oswald efficiency factor, in fact it was well below the optimised uniform swept-back-tapered configuration shown in Fig. 5, with ε = 0.981. The efficiency factor of such configuration could be further improved by starting the breakpoint position, where the second aft sweep is applied, α 3 , further inboard as shown in the contour plots in Fig. 6. This resulted in a planform tending towards the uniform swept wing in Fig. 5 and for the case of breakpoint position at 2.5% of span as shown at the bottom of Fig. 6 the Oswald efficiency factor tended towards a value of 1.

Crescent Wing

The crescent wing can be considered as a swept-tapered wing with a non-linear spanwise profile of the leading and trailing edge. Using the same optimisation process as earlier, from the top of Fig. 7 the effect of varying the sweep angle as a function of taper is shown for the planform with a parabolic and elliptic profile. The region where the Oswald efficiency factor ε > 0.95 expanded over a larger portion of the sweep and taper map as opposed to the uniform swept-tapered and non-uniform swept-tapered Mconfiguration. In fact both types of crescent wings could easily achieve Oswald efficiency factor, ε ≥ 1, for sweep angle of 10 • , where the parabolic planform attained a maximum of ε = 1.005 for a taper ratio, λ = 0.25 and the elliptic planform was marginally higher with, ε = 1.008, for a much smaller taper ratio, λ = 0.1.

For both parabolic and elliptic planforms, the sweep angle of 10 • quoted here is only applicable very close to the root of the wing as it was introduced in a way to vary the local change in curvature. However, this leads to a local change in sweep angle which can increase very rapidly along the span of the wing, whereby exceeding values of 60 • in the out-board region of the elliptic planform. On the other hand, for the elliptic wing the very low taper ratio would result in relatively low Reynolds number in the outboard region and would be prone to laminar separation. The viscous-coupled VLM method in XFLR5 could potentially account for the local change in Reynolds number, since the 2D viscous calculation from XFLR5 has been demonstrated to deal with low Reynolds number cases quite reliably. However, threedimensionality in the viscous flow due to large change in local sweep angle and the formation of large crossflow components cannot be modelled with the 2D approach of XFLR5. The combination of high sweep and low Reynolds number towards the tip would favour the formation of leading edge vortices which would generate extra lift due to the contribution from vortex lift. This prompted further studies of the flow field over the wing.

THE SPANWISE FLOW FIELD OVER THE WING

During the optimisation exercise in the previous section only the overall performance of the wing was characterised using the VLM method in XFLR5. It was demonstrated that if the sweep and taper is varied optimally the minimum induced drag condition, ε = 1.0, predicted from lifting line theory can be achieved or even surpassed. From the spanwise evolution of the local lift normalised by the total lift presented in Fig. 8, a major difference can be observed between the rectangular planform, which shows the traditional elliptical lift distribution and the rest of the configurations emanating from the current optimisation study. The swept-tapered and the M-configuration showed similar trend in spanwise loading with low inboard loading and higher out-board loading compared to the rectangular case, as seen on conven- tional swept-tapered transonic wing. However, from Tab 2, the conventional swept-tapered civil aircraft achieves an Oswald efficiency factor lower than the optimal swepttapered wing defined here which could attain a value of 1. However, the spanwise loading from the crescent wing showed a more intriguing characteristics, for the case of the parabolic planform very small or almost negligible variation in the lift was observed starting from the inboard region up to almost 60% along the span, where it increased and started dropping at approximately 95% of span with a trend similar to the swept-tapered wing. On the other hand the crescent wing with elliptic planform experienced larger extent of the constant spanwise lift, upto to 80% and unlike the parabolic planform underwent a marginal reduction upto around 95% of span where it increased rapidly again. This behaviour was predicted by the viscous-VLM method which is in fact a representation of the 2D viscous flow calculated from XFLR5 as a function of Reynolds number which reduces along the span in this case.

Along the elliptic crescent wing the sectional Reynolds number, Re, ranged from 290 × 10 3 at the root to 29 × 10 3 at the tip, for Re < 100 × 10 3 a large laminar separation bubble was predicted from the 2D XFLR5 viscous calculation at an incidence of 6 • and therefore in the presence a high local sweep the strong spanwise velocity component would generate a vortex which would contribute to the additional lift. In order, to support this hypothesis further RANS simulation was conducted on the elliptic crescent wing and a rectangular wing with similar aspect ratio for comparison. Even if the overall drag prediction was not in agreement with the experimental results, here we were only interested in the flow over the near-field of the wing and from the bottom right of Fig. 4 the spanwise loading from the RANS simulations is seen to predicted with a reasonable accuracy when compared with the experimental results. Therefore, we proceeded with RANS simulation with a mesh of approximately 13 millions cells on a rectangular wing and the optimal elliptic crescent wing. In this case the γ transition model was used in addition to the kω turbulence model. The topology of the flow over these two configurations is presented in terms of Q-criterion iso-surfaces coloured by vorticity in Fig. 9, from these two configuration the tip vortices are clearly put into evidence, however the most noticeable differences is in the region right before the tip where an additional pair of vortices can be observed in the case of the elliptic planform. The presence of this pair of vortices are clearer when comparing a plane orthogonal to the streamwise flow immediately downstrea the trailing edge in Fig. ??. This is further evidence of the presence of the vortex lift prevailing from the vortical struture similar to a leading vortex postulated earlier.

CONCLUSION

Here we have demonstrated that the VLM method in XFRL5 can provide a reliable estimate of the spanwise loading and thus induced drag when compared with existing experimental results in the literature. This was an important step towards validating the low-order modelling which was then employed to optimise wing planform geometry with lower induced drag characteristics. From a simple uniform variation of sweep and taper we showed that an Oswald efficiency factor of 1 which defines the minimum induced drag from lifting line theory can be attained. The efficiency factor could be further increased for the variation in sweep and taper applied in a nonlinear fashion to approach planform geometries found in nature. Here the crescent moon shape observed on the wings of swifts and swallows and, the caudal fins and flukes in thunniforms and caetaceans, achieved the highest efficiency factor when fitted with an elliptic leading and trailing edge chord distribution, inspired from existing knowledge of the benefit from the elliptic loading. Preliminary CFD simulations have shown that in addition to the tip vortices, for the crescent wing a combination of low Reynolds number and high sweep towards the tip favours the formation of another vortex which provides additional loading towards the tip and thus improving the efficiency of the overall wing in producing lift and thus reducing the induced drag.

  2 is very different to the crescent wing. Unlike in the case of the crescent wing the planform geometry is defined as a linear variation of the sweep and taper, as shown at the bottom of Fig 2, in this case both a backward and a forward sweep was required. The parameters s1 and s2 represent the positions where the sweep angle α1, α2 and α3 are defined, but α1 and α3 were always kept over 270 • during the study, so they were effectively backward sweep angles. The parameters imposed on the M wing geometry were the sweep and the taper ratio, a linear chord distribution along the wingspan. These yielded constraints on the root chord length (c root ), the tip chord length (c tip ) and the slope of the linear function defining the chord distribution along the wingspan (c slope ) given by equation 9, 10 and 11

Figure 2 :

 2 Figure 2: The parameters used to define the M wing geometry.

Figure 4 :

 4 Figure 4: The top left shows the crescent wing tested experimentally by van Dam [11]. The top right and bottom left shows the total lift and total drag coefficient from the experiment, LLT, VLM and CFD calculations. Bottom right, spanwise distribution of lift from the experiment of van Dam and from the VLM and CFD calculations.

Figure 5 :

 5 Figure 5: Effect of sweep angle and taper ratio on the Oswald efficiency factor, ε, on the crescent wing with parabolic (middle) and elliptic (bottom) wing planform profile.

Figure 6 :

 6 Figure 6: Effect of sweep angle and taper ratio on the Oswald efficiency factor, ε on the M-configuration

Figure 7 :

 7 Figure 7: The results of the optimisation process for the crescent with a parabolic and elliptical planform profile on the left and right respectively. The figures at the bottom shown the optimal geometry.

Figure 8 :

 8 Figure 8: Spanwise lift distribution of the main geometries studied at AoA=6 o with VLM

Figure 9 :

 9 Figure 9: Q-criterion iso-surface coloured by vorticity over rectangular wing on the left and elliptic crescent wing on the right.

Table 2 :

 2 Geometrical characteristics and efficiency of the optimized elliptical functions based crescent wing. Aircraft's data taken from[START_REF] Nita | Estimating the oswald factor from basic aircraft geometrical parameters[END_REF]


A Domínguez

J Tomé

E Gowree

FP68-AERO2022-Bauerheim Bio-inspired vortex lift for enhanced manoeuvrability

Inspired by highly manoeuvrable species of birds like the peregrine falcon and the swift, static and dynamic computational fluid dynamics (CFD) simulations were conducted to investigate vortex lift in unsteady flows. The configuration corresponds to a 50 • sweep delta wing with sharp leading edge at Re= 5.0 × 10 4 . CFD simulations were performed using a Direct Numerical Simulation (DNS) approach with a Lattice-Boltzmann Method as well as Unsteady Reynolds Averaged Navier-Stokes (URANS) simulations. Aerodynamic forces as well as the overall structure of the leading edge vortices were compared with existing literature. The evolution of the flow structures was studied when the wing performs a pitching manoeuvre from 0 • to 20 • angle of attack. Close agreement between both methods was found for the static and pitching lift curves, with the URANS solver presenting substantial limitations to capture complex unsteady phenomena such as the vortex breakdown. A time lag was observed in the flow dynamics during the manoeuvre, with the vortex breakdown delayed during pitch-up resulting in an improved aerodynamics performance, but more present and intense when pitching down. A sinusoidal motion was tested with the URANS solver and compared with the linear ramp case, showing performance advantages as well as higher similarity to real manoeuvres.

INTRODUCTION

Birds, such as the peregrine falcon and the swift, are known for their high manoeuvrability, associated with their morphological transformation capability during flight [START_REF] Re | The leading-edge vortex of swift wing-shaped delta wings[END_REF]. In particular, the peregrine falcon wing's sweep angle (Λ) can vary from Λ min = 40 • (fig. 1) to ap-proximately 90 • , which allows them to adapt to different flight phases, being able to increase their maximum speed or pull out from a dive. Figure 1: Wing shape adopted by the falcon during the final phase of the pull-out manoeuvre [START_REF] Selim | The peregrine falcon's dive: On the pull-out maneuver and flight control through wing-morphing[END_REF] Previous experimental and high-fidelity simulations by Gowree et al. [START_REF] Gowree | Vortices enable the complex aerobatics of peregrine falcons[END_REF] showed that the flow over a peregrine falcon was dominated by large vortical structures. In their subsequent study [START_REF] Selim | The peregrine falcon's dive: On the pull-out maneuver and flight control through wing-morphing[END_REF] they demonstrated that in fact the overall lift force of the peregrine falcon's wings agrees well with the delta wing lift theory, which suggests the similitude between these two configurations.

Delta wings exhibit a leading-edge vortex (LEV) at low Reynolds number, even at low incidence [START_REF] Gursul | Unsteady aerodynamics of nonslender delta wings[END_REF]. Due to the LEV, these planforms reach higher stall angles while generating higher lift than conventional wings, thus enhancing manoeuvrability. While extensive research for slender delta wings (typically with a sweep angle greater than 60 • ) has been carried out, both numerically and experimentally, considerably less investigation has been undertaken for nonslender wings (Λ < 55 • ). Depending on the Reynolds number and incidence, nonslender wings will exhibit a dual or trial vortex structure, contrasting with the less energetic secondary vortex seen on slender configurations. In addition, nonslender delta wings will exhibit early vortex breakdown, which will lead to a de-crease in lift as well as more intense unsteadiness.

Consequently, the early vortex breakdown and increased unsteadiness, characteristic of nonslender planforms, can affect the manoeuvrability and stability of these wings. Recent numerical and experimental studies by Yi et al. [START_REF] Yi | Dynamic lift characteristics of nonslender delta wing in large-amplitude-pitching[END_REF] have shown that, for a sinusoidal pitching movement, the unsteady effect becomes dominant with substantial differences when pitching up and pitching down. On one hand, a strong and coherent LEV is found when pitching up, but vortex breakdown is found close to the apex when pitching down even at angles of incidence corresponding to pre-stall on the static case. In addition, lift hysteresis is observed and found to be enhanced by greater reduced frequencies. At high pitch-up rates, the strong LEV and its delayed breakdown phenomenon create a local suction peak, which translates into an increase in lift. On the other hand, when pitching down, the lower surface of the wing becomes a suction surface, resulting in a lift loss when comparing to the static case.

Yavuz [START_REF] Yavuz | Transformation of flow structure on a delta wing of moderate sweep angle during pitch-up maneuver[END_REF] investigated the transformation of the flow structure of a 50 • swept wing undergoing a ramp pitchup manoeuvre at low Reynolds numbers of 15 000 and 20 000. One of the vortex in the dual vortex system was observed to be absent, as well as the appearance of vortex breakdown at higher angle of incidence. Due to the pitch up movement, a time lag is observed in the development of the flow structures.

The present investigation focuses on the evolution of the force and vortex structures when performing not only a pitch-up, but also a pitch-down manoeuvre, as the results and flow dynamics between both of them will differ substantially, as well as when comparing them to the static case. The scope of this study is to establish a baseline tool which can be further developed towards its application to future bio-inspired MAV's and UAV's, hence the choice of a nonslender planform and low Reynolds number configuration.

NUMERICAL METHODOLOGY

The configuration considered consists of a 50º sweep non-slender delta wing with sharp leading edge (fig. 2), at a Reynolds number of 5 × 10 4 . Two approaches are employed: (1) Unsteady Reynolds Averaged Navier-Stokes (URANS) simulations using the commercial software STAR-CCM+, and (2) Direct Numerical Simulations (DNS) with a low-dispersion Lattice-Boltzmann method (LBM).

Static numerical simulations for angles of attack of 5 • and 15 • were performed with LBM, whereas 5 , 15 • and 20 • were simulated with the URANS approach. For the pitching movement, a [0 • -20 • -0 • ] linear ramp pitchup and pitch-down movement (as shown in fig. 3) was considered, with a constant pitch rate of ±10 • /s.

Lattice-Boltzmann Method

The numerical simulations using the Lattice-Boltzmann method were carried out using the open source code Palabos [START_REF] Latt | Palabos: Parallel lattice boltzmann solver[END_REF]. Rather than focusing on the macroscopic fluid quantities, this low-dispersive method tracks the time and space evolution of a discrete particle distribution function (F i (x,t))) at a discretized velocity c i on a lattice grid. The Lattice-Boltzmann advection equation reads:

F i (x + c i ∆t,t + ∆t) -F i (x,t) = - 1 τ (F i (x,t) -F eq i (x,t)) (1)
where the right hand side of the previous equation is the collision operator, here approximated by the BGK model. This term drives the particle distribution in relation to the equilibrium state with τ as the relaxation time, associated with the viscosity of the flow ν.

It uses a cartesian mesh with a immersed boundary method that allows to have a static mesh when performing pitch manoeuvres, provided that the geometry stays in the finest grid level. The computational domain was defined as a cube with a side length equal to 40c, and the delta wing was placed in the center.

The cartesian grid contains 7 refinement levels with a total of approximately 95 million cells. The minimum and maximum cell size and time-step considered for the simulations are given in table 1. Time and mesh Independence studies were carried out to ensure correct resolution.

∆x min /c r ∆x max /c r ∆t •U 0 /c r 0.003 0.192 1.92 × 10 4
Table 1: Time and space characteristics of the LBM grid

Grid spacing influence

An initial assessment of the influence of the grid spacing was performed. In the LBM solver, due to the nature of the cartesian grid, the spacing was constant in all three coordinates. The results will be presented in relation to c/∆x, the number of cells along the chord of the delta wing. Figure 4 depicts the axial vorticity for 3 different meshes, plotted on a plane perpendicular to the vortex core. As observed, in the coarser mesh (c r /∆x = 154), the structures do not resemble a classical dual-vortex structure found in these configurations at 5 • of angle of attack [START_REF] Gordnier | Compact difference scheme applied to simulation of low-sweep delta wing flow[END_REF]. In fact, the expected behaviour is only found for the higher resolution (c r /∆x = 333).

URANS

A set of unsteady RANS simulations in STAR-CCM+ were conducted in parallel with the LBM counterparts. The main objective is to use the force evolution data and main flow behaviour if validated with the DNS due to the lack of experimental results in the literature. Although it is expected that this approach will lack the resolution of the latter, with difficulties capturing complex unsteady flow dynamics like the vortex breakdown and its evolution during manoeuvres.

The meshing domain consists of a moving region containing the wing geometry and space where the main vortical activity is expected, placed inside an outer static domain. Both meshes interact through an overset interface, where the moving region is contained within the smallest refinement level of the static one throughout the rotating motion. For better consistency and coherence when comparing the results, both static and dynamic simulations have been carried out using this meshing approach, even though for the static case it would not be required.

An unstructured trimmer mesh with no prism layer has been used for the case, as it will consist essentially on a separated flow. The total number of cells is around 8.4 million, with only 1.4 million corresponding to the outer static region and the remaining 7 million to the inner overset mesh.

The outer static domain applies a coarse mesh with 6 levels of refinement and 4 layers per level. The two finest levels have been stretched downstream for a better capturing of the wake, with the cell size of the most refined matching the one at the outer boundaries of the inner overset mesh (figure 5). The moving overset mesh follows the same approach, although with a considerably lower cell size. This region presents 4 refinement levels ranging between 8 and 12 layers each, the finest being located around the sharp leading edge, with a cell size dx min = 0.488 mm. An extensive mesh independence study has been carried out to arrive to this point, where different domain and refinement shapes, meshing models, growth rates and cell size values have been tested. A time-step of 0.01 s has been used for the implicit unsteady solver of the static simulations, as it has been proven to be sufficient to capture the unsteadiness associated to the large structures captured by the URANS solver. This value was reduced to 0.005 s for the rotation in the dynamic cases after an initial static convergence, as to increase the number of time-steps per degree of pitch to 20.

A second order scheme in space and time has been chosen, and the turbulence model used is the Shear-Stress-Transport (SST) kω, as it provided closer results to both the literature and LBM simulations, while also better capturing the unsteady behaviour of the flow, unlike other tested options like the kε model.

RESULTS AND DISCUSSION

Static simulations 3.1.1 Lift coefficient

Reliable experimental or high fidelity numerical simulation data is scarce in literature for nonslender delta wings, which limited a thorough verification of the results. According to Polhamus [START_REF] Ec | Predictions of vortex-lift characteristics by a leading-edge suction analogy[END_REF], using the leading-edge suction analogy, the vortex-lift has a contribution to the total lift which decreases for lower sweep angles. Nevertheless, even for low angles of attack and low Reynolds number, a dual and trial vortex structure is found in these configurations, resulting in a non-negligible contribution to the overall lift which contradicts the basic assumptions of the theory. In addition, a Reynolds dependency has been observed, suggesting that flow reattachment has an influence on stall [START_REF] Gursul | Unsteady aerodynamics of nonslender delta wings[END_REF]. [START_REF] Yi | Dynamic lift characteristics of nonslender delta wing in large-amplitude-pitching[END_REF]. By comparing the two experimental curves for Λ = 50 • , one can observe that the maximum lift coefficient C L max increases for smaller Reynolds numbers. However, this also results in a lower α max and a more sudden stall, with a rapid post-stall decrease in C L . According to Gordnier et al. [START_REF] Gordnier | Compact difference scheme applied to simulation of low-sweep delta wing flow[END_REF], the increase in Re leads to the strengthening of the dual or trial vortex structure (depending on the incidence), which may lead to higher suction on the surface and therefore higher C L .

Comparing the curves corresponding to different sweep angles allows to deduce the effect of this parameter on the lift coefficient, where a highly swept delta wing will present increased C L max and α max . However, this will also result in a lower slope and therefore lower C L at smaller angles of attack, which could influence the planform choice depending on the desired performance when aiming at a potential application.

The results obtained from the URANS and LBM solvers show a good agreement with the experimental data corresponding to the Λ = 50 • , Re=100 000 case [START_REF] Gs | Passive flow control over flexible nonslender delta wings[END_REF], even though performed at different Reynolds number. This comparison is justifiable as both Reynolds numbers are outside the sensitive range, presented by Gursul et al. [START_REF] Gs | Passive flow control over flexible nonslender delta wings[END_REF] as a region below Re=30 000 where the results are strongly influenced by variations in Re. Close to the apex of the wing, a thin shear layer develops, which will later constitute the primary vortex with the greatest intensity. Before reaching half of the chord, the development of a secondary vortex, with the same sign as the primary but less intense, is seen. The existence of this dual vortex structure has been confirmed both numerically and experimentally for nonslender planforms [START_REF] Gursul | Unsteady aerodynamics of nonslender delta wings[END_REF]. The vortex breakdown can be seen in the aft section of the wing, where the loss of coherence of the primary leading-edge vortex is observed, resulting in a highly unsteady behaviour but without causing the wing to stall at this stage.

Mean flow structure

The streamwise vorticity at two sections along the chord has been plotted for both LBM and URANS approaches (figure 9), in order to compare the capacity to capture the complex flow structures of the latter with the DNS case. Clear differences can be observed for both x/c = 0.5 and x/c = 0.8 planes. The first one shows a double vortex structure for the LBM case, while only one large primary vortex can be appreciated with URANS. As seen in figure 9a, the vortex breakdown has already occurred at the aft section, however, the flow structures at this angle of attack are not large enough to be captured by the URANS solver, which now shows an averaged double vortex structure with no signs of breakdown.

Figure 10 presents a side-by-side comparison between literature data [START_REF] Gordnier | Compact difference scheme applied to simulation of low-sweep delta wing flow[END_REF] and one of the static URANS simulations, both corresponding to an angle of attack α = 15 • . The surface streamlines are plotted over the wing surface, coloured with the pressure coefficient C p , which allows to observe the relation between the vortical structures the suction they produce over the body. Despite the Reynolds number differing by a factor of 2, both cases are considered as low Re regime, thus it can be observed that the level of agreement with the literature is quite high, with the same vortical structures appearing and similar overall flow behaviour. The primary leading edge vortex can be clearly distinguished, with very close results between both cases in terms of vortex size and angle of the attachment line (PA) with respect to the x-axis. Also, a smaller tertiary vortex can be observed between the primary one and the leading edge of the wing, which breaks down at approximately half of the chord for the Re=26 000 case, and around one third of the chord for the Re=50 000 simulation, suggesting that indeed the vortex breakdown position is affected by the Reynolds number, advancing to-wards the apex as Re increases. Figure 11: Lift coefficient obtained for the manoeuvre as well as for static experimental and numerical results [START_REF] Gs | Passive flow control over flexible nonslender delta wings[END_REF].

Following an initial transient state when initiating the movement, the lift coefficient follows an almost linear trend. In fact, comparing with the static data, the slope of the curve seems to be similar, with a translation ∆C L ≈ 0.2.

When performing a constant pitch up or pitch down manoeuvre, the delta wing sees a change in the effective angle of attack due to the introduction of a vertical velocity component. In this case, denoting Ω [ • /s] as the pitch rate with axis on the centre of the wing and positive in the y-direction, and α as the angle of attack, the change in local incidence would be of the order of magnitude of α + 0.5Ω. This is consistent with the data found for C L in figure 11, where the manoeuvre curves at a certain α would match static curves at α + 5 • (pitch-up) or α -5 • (pitch-down) for the current case. Nevertheless, a substantially higher value is found for C L at α = 20 • when pitching up (C L ≈ 1.3) when comparing to the static case (C L ≈ 0.9), suggesting that vortex breakdown and stall are delayed when the wing is manoeuvring. This will be explained in the following section, while analysing the evolution of the vortical structures.

An extra step has been taken following the previous linear ramp manoeuvre, replacing the constant pitch rate by a sinusoidal movement. This could not only ease the simulation convergence at the transition between phases of the manoeuvre, but also improve the performance of the rotating wing. To make the comparison between both approaches possible, the frequency of the movement f=0.25 Hz has been chosen to provide the same average pitch rate along the manoeuvre as the previous case (Ω avg = 10 • /s). This results in a reduced frequency k = π f c/U ∞ = 0.785, far from the quasi-unsteady range as k > 0.1 thus allowing the dynamic effects on the flow with respect to the static cases. The angle of incidence along the manoeuvre is

α = 10π 180 × 1 -cos π 2 (t -t 0 ) [rad], (2)
which gives

Ω = 10π 180 × π 2 sin π 2 (t -t 0 ) [s -1 ]. (3)
Figure 13 shows the sinusoidal C L curve following a smooth hysteresis loop which matches the linear ramp case within the range of α = 0 • and α = 20 • . The main difference between the two approaches lies in the middle part of the manoeuvre, where the sinusoidal shows a much better pitch-up performance, in exchange for a reduced C Lmax and a larger lift drop at the beginning of the pitch-down phase, even reaching negative values towards the end of the cycle. Studies by Yi et al. [START_REF] Yi | Dynamic lift characteristics of nonslender delta wing in large-amplitude-pitching[END_REF] confirm this behaviour for large reduced frequencies in the unsteady range (k > 0.1).

Due to software limitations, it has not been possible to reproduce this non-constant pitch rate manoeuvre for the LBM approach, but the close agreement observed for the linear ramp case, seen in figure 11 has led to believe the URANS results are reliable enough for this preliminary analysis.

Evolution of the vortical structures

To understand the behaviour of the vortical structures over the wing when performing the manoeuvre, the position of the primary LEV is tracked. For that purpose, two positions have been analysed for the ramp manoeuvre: x/c = 0.4 and x/c = 0.6, where b is defined as the local semi-span at each location.

In figure 14, the position of the primary vortex in relation to the centre chord, measured in y/b, is depicted for the pitch-up/pitch-down manoeuvre, where y is the vortex core position in the y-axis at each section. As observed for the static case, when α increases, the vortex core moves closer to the centre chord of the delta wing, as expected. During the pitch-up manoeuvre, the trends captured with the LBM solver and the URANS solver are identical, however, the values obtained differ by around 10%. This evolution is depicted in figures 17 and 18, which illustrate the LEV development with Q-criterion iso-surfaces, where the primary vortex moves inward and When pitching down, the vortex core does not move closer to the leading edge, and the flow becomes substantially different from what is observed in the static case. Due to its early vortex breakdown, a consistent primary vortex is only seen up to x/c = 0.5, reason why the position of the vortex core is only plotted for x/c = 0.4. In this case, the differences between the two solvers become clear. For the URANS solver, the primary vortex position in figure 14 is nearly constant, as for the LBM solver, the primary vortex moves closer to the symmetry plane of the wing, showing a difference in spanwise position ranging from 20% to 40% between both solvers.

In addition, as seen in figure 17, vortex breakdown advances to the apex as the wing is pitching down. This highly unsteady phenomenon is not seen as intensively for URANS as for the DNS solver, clearly observed when comparing the vortex development in figures 17 and 18.

The height of the primary vortex in relation to the wing's surface is also tracked and analysed, and its evolution is depicted in figure 16. Contrary to what was observed for its y-position, the overall trend captured for both solvers is in good agreement. For the pitch-up phase, the primary vortex moves away from the wing surface as its dimension increases with the higher incidence. As it starts to pitch down, the initial height at maximum incidence is conserved and it moves closer to the surface as the angle of attack decreases.

It is evident that the LBM solver is able to resolve the unsteady features of vortex breakdown, characteristic of the pitch-down phase of the delta wing, as opposed to the URANS solver. However, an excellent level of agreement is seen for the evolution of the lift coefficient, which provides confidence for further analyses using purely the URANS approach.

In order to further understand how the vortical structures evolve along the manoeuvre, the streamwise vorticity has been analysed and compared with the static data at x/c = 0.4 and x/c = 0.6, for α = 5 • and α = 15 • . At α = 5 • , a thin shear layer can be seen at x/c = 0.4 when pitching up as the primary leading-edge vortex is being formed, as seen in figure 15. At x/c = 0.6 the primary vortex starts to be clearer. During pitch-up, a delay between the static and dynamic case is observed. This causes the LEV to adhere to the suction surface of the wing until higher angles of attack, providing higher lift overshoot, confirmed experimentally by Yi et al. [START_REF] Gs | Passive flow control over flexible nonslender delta wings[END_REF].

For that reason, a dual vortex structure is not seen at α = 5 • when pitching up, compared with the static case. A small weaker vortex appears inward of the primary one for the LBM solver, seen in figure 17 for α = 5 • , which is rapidly dissipated.

During the pitch-down phase, a triple vortex structure is seen at α = 5 • at x/c = 0.4. This type of structure, although closer to the leading edge of the wing, is typically spotted at higher static angles of attack. This suggests that there is also a delay in the recovery of the vortical structures when pitching down. At x/c = 0.6, vortex breakdown has occurred and the coherent vortical structure is no longer retrieved downstream.

At α = 15 • , depicted in figure 15, the same trend is observed. When pitching up, the vortical system starts to resemble a triple vortex structure described by Raymond et al. [START_REF] Gordnier | Compact difference scheme applied to simulation of low-sweep delta wing flow[END_REF], with lower intensity. Once again, this suggests that there is a time delay in the formation of the LEV due to the manoeuvre. During the pitch-down phase, as seen in figure 17, vortex breakdown moves very close to the apex of the wing. For that reason, a loss of coherence of the LEV downstream is expected, which can be confirmed in figure 15. At x/c = 0.4, three co-rotating vortex cores can be identified, which are rapidly dissipated as seen at x/c = 0.6.

CONCLUSION

The vortical flow over a 50 • sweep delta wing with sharp leading edge was simulated with an unsteady RANS solver in STAR-CCM+ and a DNS approach with a low-dispersion LBM. Static simulations at Re=50 000 and different angles of attack showed very close agreement between both methods in terms of force calculation, but the less computationally expensive URANS solver presented limitations when attempting to capture highly unsteady phenomena like the vortex breakdown, especially at low angles of attack due to the lower size of the unsteady flow structures.

Dynamic simulations for a 0 • -20 • linear pitchup/pitch-down manoeuvre were performed with both solvers at a constant pitch rate Ω = ±10 • /s. an excellent agreement was observed between the lift curves corresponding to the two approaches, where these were found to follow a similar trend as the static one but with a fixed C L translation. A time lag was observed in the development of the flow structures, where the vortex breakdown was delayed during pitch-up, effectively increasing the maximum lift coefficient and its corresponding angle of attack. During the pitch-down phase, the vortex breakdown became more intense and advanced towards the apex, resulting in a lift drop and a lower performance of the URANS solver when measuring the vortex position along the manoeuvre. A sinusoidal approach to the pitching manoeuvre, with a reduced frequency k = 0.785, was tested with the URANS solver and compared with the linear ramp, resulting in a C L curve following a smoother hysteresis loop which eased simulation convergence and improved manoeuvring performance during most of the pitch-up motion.

Figure 2 :Figure 3 :

 23 Figure 2: Top view of the geometry (left) and detailed view of the leading edge (right)

Figure 4 :

 4 Figure 4: Axial vorticity obtained with LBM, at x/c = 0.5, plotted on a plane perpendicular to the vortex core for 3 different meshes: a) -c r /∆x = 154, b) -c r /∆x = 205 and c) -c r /∆x = 333 .

Figure 5 :

 5 Figure 5: Rotating overset mesh inside most refined level of static domain in the x-z plane. The wing is represented at α = 20 • .

Figure 6 :

 6 Figure 6: Close-up view of the computational domain in the x-z plane. The delta wing is represented at its maximal angle of attack at α = 20 • . Only 6 refinement level are represented in this figure.

Figure 7 :

 7 Figure 7: Lift coefficient obtained with URANS and LBM simulations at Re=50 000 compared with experimental results at Re=20 000 and Re=100 000 [1][9][11].

Figure 7

 7 Figure 7 depicts the lift coefficient C L as a function

Figure 8

 8 Figure8depicts the axial vorticity at several planes perpendicular to the vortex core. At 5 • angle of attack, the vortex structure found using the LBM solver is consistent with the literature.

Figure 8 :

 8 Figure 8: Mean vortex structure over the wing for α = 5 •

Figure 10 :

 10 Figure 10: Surface streamline pattern and pressure coefficient at α = 15 • from [2] at Re=26 000 (up) and STAR-CCM+ at Re=50 000 (down).

3. 2 Figure 9 :

 29 Figure 9: Streamwise vorticity with LBM and URANS for α = 5 • at axial positions a) x/c = 0.5 and b) x/c = 0.8

Figure 13 :

 13 Figure 13: Lift coefficient obtained for the sinusoidal and linear ramp manoeuvres obtained with URANS

Figure 14 :

 14 Figure 14: Evolution of the primary vortex position y/b at axial positions x/c = 0.4 (■) and 0.6 (+) for LBM (orange) and URANS (black) during pitch-up (-) and pitchdown (--) manoeuvres.

Figure 12 :

 12 Figure 12: Streamwise vorticity obtained with LBM at two axial positions x/c for α = 5 • at a) -pitch-up, b) -static and c) -pitch-down

Figure 15 :

 15 Figure 15: Streamwise vorticity obtained with LBM at two axial positions x/c for α = 15 • at a) -pitch-up and b) -pitchdown

Figure 17 :

 17 Figure 17: LEV development over the delta wing shown by Q-criterion iso-surfaces coloured by the velocity magnitude obtained for the LBM solver. The arrows indicate whether the delta wing is pitching up or down.

Figure 18 :

 18 Figure 18: LEV development over the delta wing shown by Q-criterion iso-surfaces coloured by the velocity magnitude obtained in STAR-CCM+. The arrows indicate whether the delta wing is pitching up or down.


Francesco Scarano
email: francesco.scarano@isae-supaero.fr

Marc C Jacob

Xavier Carbonneau

Erwin R Gowree

Characterisation of bursts in a turbulent boundary layer over circular cavities

Boundary layer surveys are performed to characterise the turbulent boundary layer grazing over a a surface with flush-mounted circular cavities disposed in a staggered arrangement (staggered angle of 45 deg). The momentum thickness based Reynolds number varies from Re θ = 1830 to 3380 and the cavity diameter and spacing in wall units from d + = 130 to 250 and L + = 587 to 1075 respectively. A decrease in the local skin friction drag with respect to a smooth baseline is evidenced as well as a thickening of the viscous sublayer. The bursts frequency profile for the perforated case is shifted away from the wall while the intensity of the bursts is reduced. The production term of the turbulent kinetic energy budget is obtained from Particle Image Velocimetry data. An upward shift and a decrease in magnitude of the peak associated to the streaks suggests that the cavities modify the near wall turbulent cycle.

INTRODUCTION

Skin friction drag accounts for more than (50%) of the overall drag for an aircraft in cruise condition. Reducing the skin friction drag would then have a strong impact on the CO 2 emissions reductions since 1% of drag reduction is converted into a 0.75% in fuel-burn savings [START_REF] Ricco | A review of turbulent skin-friction drag reduction by near-wall transverse forcing[END_REF][START_REF] Michael | Friction-drag reduction by transverse wall motion -a review[END_REF]. The skin friction drag is generated by the turbulent boundary layer grazing over the surfaces of the aircraft. In a turbulent boundary layer the effect of the coherent structures, namely the hairpin vortices largely documented in literature [START_REF] Adrian | Hairpin vortex organization in wall turbulence[END_REF], is responsible for the majority of the turbulent kinetic production and transport in the near wall region and the skin friction drag. In particular, the turbulent activity, sometimes called bursting, can be decoupled in sweep (u > 0, v < 0) and ejections (u < 0, v > 0) which give a positive contribution to the Reynolds shear stress and to the production of turbulent kinetic energy. The control of this near-wall activity would lead to a potential skin friction reduction [START_REF] Michael | Friction-drag reduction by transverse wall motion -a review[END_REF].

Several techniques have been proposed for skin friction drag reduction; the most studied are the riblets [START_REF] Walsh | Drag characteristics of v-groove and transverse curvature riblets[END_REF][START_REF] Choi | Near-wall structure of a turbulent boundary layer with riblets[END_REF] and the spanwise wall oscillations [START_REF] Ricco | A review of turbulent skin-friction drag reduction by near-wall transverse forcing[END_REF][START_REF] Choi | Turbulent boundary-layer control by means of spanwise-wall oscillation[END_REF][START_REF] Quadrio | Streamwise-travelling waves of spanwise wall velocity for turbulent drag reduction[END_REF]. Despite the different mechanism involved (the riblets inhibit the spanwise movement of the structures while the spanwise wall oscillations impose a spanwise movement) and the charateristic size in wall units (15-20 viscous lengths for the riblets sizes and up to 8000 viscous lengths for the spanwise wall oscillations wavelengths) these techniques have proven to be highly effective in turbulent skin friction drag reduction. Difficulties in the implementations of these techniques lay in the sophisticated manufacturing for the riblets and in the high complexity realisation of spanwise wall oscillations in real life conditions as well as in laboratories [START_REF] Auteri | Experimental assessment of drag reduction by traveling waves in a turbulent pipe flow[END_REF][START_REF] Marusic | An energy-efficient pathway to turbulent drag reduction[END_REF].

More recently, studies on large dimples [START_REF] Van Nesselrooij | Drag reduction by means of dimpled surfaces in turbulent boundary layers[END_REF]) (diameters up to four times the boundary layer thickness) and small cavities [START_REF] Gowree | Skin friction drag reduction over staggered three dimensional cavities[END_REF][START_REF] Silvestri | Attenuation of sweep events in a turbulent boundary layer using micro-cavities[END_REF] (diameter in wall units 20 < d + < 145) have boosted the interest in the application of these wall patterns to reduce the turbulent activity and the skin friction.

The aim of this study is to further investigate the modification of the turbulent activity of a turbulent boundary layer grazing over a staggered array of flush-mounted cavities. This is done by the combination of hot wire surveys and stereoscopic PIV. The diameter and depth of the cavities have been chosen in agreement with the previous studies [START_REF] Gowree | Skin friction drag reduction over staggered three dimensional cavities[END_REF][START_REF] Silvestri | Attenuation of sweep events in a turbulent boundary layer using micro-cavities[END_REF].

METHODS

Hot wire setup

The experimental campaign has been conducted in the SaBRe wind tunnel at ISAE-SUPAERO. The length of the test section is 2.5 m in the streamwise direction (X), its width is 1.2 m in the horizontal (spanwise) direction (Z), whereas its height is 0.8 m (Y ). Inside the test sec-tion, a 2.5 m long, 1.2 m wide flat plate was mounted horizontally; the plate has an elliptic leading edge and it is equipped with a flap at the trailing edge. The flap was set at an angle of 3.7°during the entire experimental campaign in order to ensure a zero streamwise pressure gradient along the measurement domain. A Dantec P5517 hot wire probe for boundary layer connected to a two axes traversing system coupled with an optical probe to wall alignment technique proved to be a reliable method for measurements close to the wall.

PIV setup

Stereoscopic PIV was used in the current setup. The measurement was in the direction of the flow (XY ) at the trailing edge of the cavities. A commercial PIV system from Dantec Dynamics is used for the measurements. The laser is placed on top of the test section and it is connected to a laser arm which allows a precise movement and alignment of the laser sheet along the domain. The laser sheet is shone from the top and passes through a glass ceiling of the wind tunnel. The thickness of the laser sheet is approximately 0.5 mm. To avoid reflections on the models, the bottom part (which is transparent because of the transparent plexiglass) is covered by a black matte film. Black velcro was applied to the surfaces of the tunnel to avoid reflections during the acquisition process. The flow is seeded with a TOPAS smoke generator located downstream of the divergence of the tunnel. Before acquiring each measurements. The seeding system was activated for some time prior to the measurement in order to ensure that the boundary layer was sufficiently seeded.

The images are taken from the sides through the glass wall by two CCD cameras disposed in a steroscopic arrangement equipped with a NIKON 300 mm f/4 lens. DANTEC dynamic studio 6.1 was used for the calibration, for synchronising image acquisition to the laser, and for the calculation of the PIV vectors. The field of view (FOV) is 8.5 mm×4.5 mm. 1000 image pairs are acquired at 50 Hz. The two laser pulses have a ∆t which depends on the freesetream velocity and the measurement location, and it was such that a maximum particle displacement in a range of 10 pixels is obtained. For the streamwise planes the ∆t used are 15 µs for 10 m/s, 20 µs for 15 m/s and 30 µs for 20 m/s. The image pairs are processed using the 2D3C crosscorrelation PIV algorithm of Dantec Dynamics studio. The interrogation window size is iteratvely changed passing from 32 × 32 pixel to 16 × 16 with an overlap of 50%. This leads to a grid of 40400 vectors with a spatial resolution of 0.21 mm (∆X + , ∆Y + = 4). Vectors close to the wall Y + are discarded due to the reflections at the wall.

Models and measurements locations

The test samples consist of two square Plexiglas plates of 400 mm by 400 mm and 4 mm thickness, one completely smooth and the other with a staggered array of circular cavities sealed at the bottom. The diameter d of the cavities is 5 mm. The spacing between two cavities on a same row is L =22 mm (4.4d), with a stagger angle of 45°. The characteristic dimensions of the cavities and the three flow conditions are reported in Table 1. The hot wire measurements have been performed in the mid-span plane of the model, upstream of the cavities at the middle and downstream of the cavities where the PIV plane has been acquired. The measurement at the leading edge and at trailing edge are at one cavity distance L upstream/downstream of the last row of cavities. In Fig 3 the mean and RMS velocity profiles downstream the last row of cavities are compared with the corresponding smooth results. The values are reported in wall units and they are obtained using the actual value of the friction velocity calculated using the Clauser chart technique [START_REF] Wei | Comment on the clauser chart method for determining the friction velocity[END_REF]. For the sake of brevity only the results for the 10 m/s case (namely d + = 130) are reported in Fig 3. Numb

RESULTS AND DISCUSSIONS

Y + a) b) c) d) e) f) g) h) i) l) m) n) 0 1 2 3 4
•10 4 Y + smooth perforated

Number of bur

Y + a) b) c) d) e) f) g) h) i) l) m) n)
10 2 10 3 Y + smooth perforated a) b) c) d) e) f) g) h) i) l) m) n) 0 1 2 3 4
•10 4 Y + smooth perforated

10 2 10 3 Y + smooth perforated a) b) c) d) e) f) g) h) i) l) m) n) 0 1 2 3 4
•10 4 Y + smooth perforated A thickening of the viscous sublayer can be evidenced as well as a decrease in the peak of the RMS of the streamwise velocity. These results are in agreement with [START_REF] Gowree | Skin friction drag reduction over staggered three dimensional cavities[END_REF] and [START_REF] Silvestri | Attenuation of sweep events in a turbulent boundary layer using micro-cavities[END_REF] and suggest a local skin friction reduction. In Fig 2 the local skin friction reduction percentage between the perforated and the baseline surface is reported. A local skin friction reduction can be observed for all the freestream conditions but a decrease of the skin friction reduction with d + can be evidenced; similar trend are re-ported in the studies by [START_REF] Gowree | Skin friction drag reduction over staggered three dimensional cavities[END_REF] and [START_REF] Silvestri | Attenuation of sweep events in a turbulent boundary layer using micro-cavities[END_REF].

10 2 10 3 Y + s p a) b) c) d) e) f) g) h) i) l) m) n) 0 1 2 10 1
To better understand the turbulent behaviour, the Variable Interval Time Averaging (VITA) technique has been applied to the profiles [START_REF] Blackwelder | On the wall structure of the turbulent boundary layer[END_REF]. This technique, largely used for turbulent boundary layer research, is a simple and powerful way to detect bursts in a turbulent boundary layer. It consists in comparing the short-term RMS to the long-time RMS of the entire signal. If the short-term RMS is larger than the long-term RMS multiplied by a constant k [START_REF] Sullivan | Coherent structure identification from the analysis of hot-wire data[END_REF] a burst is detected. The short-term RMS is computed considering a time window which has to be of the order of magnitude of the characteristic time of the passage of a structure T + ∼ = 15. Given a fluctuating quantity Q(x i ,t) the short term average of the variable in the interval time T can be defined as

Q(x i ,t, T ) = t+ T 2 t-T 2 Q(x i , s)ds (1)
and the short term variance is

var(x i ,t, T ) = U 2 (x i ,t, T ) -[ U(x i ,t, T )] 2 . ( 2
)
The detection is then based on a Heaviside function

D(t) = 1 if var(x i ,t, T ) > k •U 2 RMS ; 0 otherwise (3)
The VITA technique allow to detect, for each location of the boundary layer, namely for each Y + , the number of burst and to compute the average burst intensity. This leads to a profile of the bursts frequency and intensity.

In Fig 4 the results for the d + = 130 and the comparison with the smooth baseline are presented (similar consideration can be made for the other flow conditions). In Fig 4 a) one can see that the burst frequency profile for the perforated case is shifted upward by ∆Y + = 5. A similar shift has been documented previously for spanwise wall oscillations [START_REF] Choi | Turbulent boundary-layer control by means of spanwise-wall oscillation[END_REF][START_REF] Maria | Particle image velocimetry investigation of a turbulent boundary layer manipulated by spanwise wall oscillations[END_REF] and riblets [START_REF] Choi | Near-wall structure of a turbulent boundary layer with riblets[END_REF]. This is in line with the thickening of the boundary layer and suggest an upward shift of the turbulent activity.

When looking at the burst intensity profile, Fig 4 b), a decrease in the average intensity with respect to the smooth baseline can be observed. This is very clear in the region which corresponds to the lobe at 20 < Y + < 80.

The percentage intensity decrease with respect to the smooth baseline is reported in Fig 5 for the three flow conditions. As for the friction coefficient, the intensity reduction is more pronounced for the lower d + and decreases when increasing the d + namely the Reynolds number. This results is in agreement with what found by Silvestri et al. [START_REF] Silvestri | Attenuation of sweep events in a turbulent boundary layer using micro-cavities[END_REF] who applied the technique only for Y + = 100 (see red dots in In Fig 6 a) a snapshot of the streamwise velocity contour is reported, again only for d + = 130. In the region downstream of the last row of cavities, the turbulent kinetic energy production term is evaluated

P = uv dU dY . (4)
The production term is averaged in the region downstream of the cavity and the average profile is reported in Fig 6 b). An upward shift and an attenuation of the first peak of the turbulent kinetic energy production term has to be highlighted. The first peak is associated with the near wall cycle and the kinetic energy produced by the streaks close to the wall. Interestingly, the second peak which is associated with the shear layer of the boundary layer is increased in magnitude suggesting an overall upward shift of the turbulent kinetic energy production. These results are in agreement with the VITA results presented in Fig 4.

Quadrant Analysis in the downstream region (results not reported here) revealed that the attenuation of the inner peak is due to an attenuation of the sweeps closer to the wall. Conversely the ejections increases in the logarithmic layer justifying the increase of the outer peak.

CONCLUSION & FURTHER WORKS

Turbulent boundary layer measurements have been conducted on a surface with flush-mounted circular cavities disposed in a staggered arrangement. The measurements have been conducted by means of hot wire surveys and stereoscopic PIV at three different flow conditions resulting in d + between 130 and 250. The cavities modify the near wall turbulent cycle, shifting upwardly and attenuating the production term of the turbulent kinetic energy

Figure 1 :

 1 Figure 1: Hot wire setup a), PIV setup b), test models and measurement locations c) (red dots for hot wire measurements, green line for the trace of the PIV plane).

Figure 2 :

 2 Figure 2: Percentage skin friction reduction obtained using the Clauser chart technique.

Figure 3 :

 3 Figure 3: Mean velocity and root mean square of the streamwise velocity in wall units at the trailing edge of the model, d + = 130.

Figure 4 :

 4 Figure 4: Number of bursts profile a), average burst signature for smooth and perforated b), average burst intensity profile for smooth and perforated c), results for d + = 130.

Figure 5 :

 5 Figure 5: Intensity reduction profile for the three flow conditions.

Fig 5 ).

 5 Their study highlighted the loss of benefit in terms of intensity reduction for values of d + above 145. On contrary, the results reported in Fig 5 show that even for higher values of d + (190 and 250) there is an entire region of the boundary layer (20 < Y + < 80 highlighted in blue) where the intensity reduction is positive. These results are in agreement with the local skin friction measurements reported in Fig 2.

Figure 6 :

 6 Figure 6: Contour of a snapshot of the velocity field a), average turbulent kinetic energy production term b).

  budget. This is due to an overall attenuation of burst intensity in the region 20 < Y + 80 and an upward shift of the burst frequency. The modification of the turbulence is translated into a thickening of the viscous sublayer and a local skin friction reduction. Further works are encouraged to 1. verify the spanwise uniformity of the aforementioned behaviour 2. study separately the effect of geometrical parameters (d + , L + ) and the flow condition.

Table 1 :

 1 Flow and geometrical parameters for the three different test conditions, downstream measurement location.

	CASE U ∞ [m/s] Re θ	d +	h + Re τ	L +	δ /d δ /h	u τ 0
	a)	10	1830 130 104 580 587	4.3	3.4 0.42
	b)	15	2710 190 152 745 838	3.9	3.1 0.60
	c)	20	3380 250 200 950 1075 3.7	3.0 0.77

ACKNOWLEDGEMENTS

The authors thank the AID (Agence Innovation Défense) for their financial support to the FriDA project.