1. Introduction
Flapping wings represent an intriguing, biologically inspired concept for aero/hydrodynamic force generation, essentially devoted to exploit the optimal capabilities of flying/swimming animals. Such capabilities include, in particular, high propulsive efficiency, agile manoeuvring and hovering at a fixed position (Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010; Haider et al. Reference Haider, Shahzad, Qadri and Shah2021). Furthermore, when compared with other concepts (i.e. based on fixed or rotating wings), flapping wings have the potential to overcome certain challenges appearing at relatively low Reynolds number, such as the drop in performance when viscous effects are not confined to thin boundary layers (Li, Dong & Cheng Reference Li, Dong and Cheng2020). From an engineering viewpoint, flapping wings can be especially employed in small-scale devices working at relatively low Reynolds numbers, such as microaerial vehicles (MAVs) (Kumar & Michael Reference Kumar and Michael2012) as well as autonomous underwater vehicles (AUVs) (Triantafyllou, Techet & Hover Reference Triantafyllou, Techet and Hover2004). Such devices find relevant applications in, for example environmental monitoring, surveillance or search-and-rescue activities. Moreover, another potential employment of flapping wings is represented by (actively driven) flow energy harvesting (Young, Lai & Platzer Reference Young, Lai and Platzer2014).
During the last decades, research efforts have been devoted to the development of mature flapping-wing MAV/AUV prototypes (Licht, Hover & Triantafyllou Reference Licht, Hover and Triantafyllou2004; de Croon Reference de Croon2020). In parallel, fundamental research has been carried out to better understand the flow physics involved in these complex unsteady aero/hydrodynamic problems, such as the characteristic formation, development and breakup process of leading-edge vortices (LEVs), which are typically the main things responsible for the generated thrust and lift force (Platzer et al. Reference Platzer, Jones, Young and Lai2008; Moriche, Flores & García-Villalba Reference Moriche, Flores and García-Villalba2016; Eldredge & Jones Reference Eldredge and Jones2019; Wu et al. Reference Wu, Zhang, Tian, Li and Lu2020). Despite these advances, however, it can be noted how the current knowledge is essentially limited to rather idealised conditions, i.e. assuming that the free stream interacting with the body is a laminar and unperturbed uniform flow. In realistic configurations, on the other hand, the environmental flow surrounding the body is generally perturbed and, in many cases, properly turbulent (Watkins et al. Reference Watkins, Milbank, Loxton and Melbourne2006). How environmental perturbations (e.g. atmospheric or oceanic turbulence) affect the aero/hydrodynamic mechanisms of a biologically inspired propulsion system thus remains a largely unaddressed topic and constitutes the basis for the present work.
Microaerial vehicles pose peculiar issues related to their characteristic flight envelope, which involves low flight speeds and low operational altitudes. The latter typically lie in the so-called atmospheric boundary layer (ABL), i.e. the region directly influenced by the surface frictional effect and characterised by strong fluctuations and mixing (Mohamed et al. Reference Mohamed, Clothier, Watkins, Sabatini and Abdulrahim2014a ). As a result, they can experience remarkably high turbulence intensities (here, the turbulence intensity can be defined as the ratio between the root mean square of the fluid velocity fluctuations and the locomotion speed). A similar scenario can be generally expected also for AUVs operating in realistic aquatic environments. In fact, in this case wave-induced coherent perturbations can be present along with substantially random ones (Dantas, da Cruz & de Barros Reference Dantas, da Cruz and de Barros2013). In any case, the design of robust and accurate controllers is key for achieving satisfactory navigation capabilities, especially in the case of severe external disturbances (Guerrero et al. Reference Guerrero, Torres, Creuze and Chemori2020; Panda, Mitra & Warrior Reference Panda, Mitra and Warrior2020). Focusing on flapping-wing devices, it is also important to point out that their size and flapping period can be often comparable to the characteristic spatial and temporal scales of the environmental flow perturbations, respectively (Watkins et al. Reference Watkins, Milbank, Loxton and Melbourne2006; Fisher et al. Reference Fisher, Ravi, Watkins, Watmuff, Wang, Liu and Petersen2016).
The importance of aero/hydrodynamic performance and control of MAVs/AUVs (not necessarily biologically inspired, but also of more conventional kinds) in gusty or perturbed flow environments has been highlighted by some studies in the literature (Mohamed et al. Reference Mohamed, Clothier, Watkins, Sabatini and Abdulrahim2014a , Reference Mohamed, Watkins, Clothier, Abdulrahim, Massey and Sabatinib ; Ortega-Jimenez et al. Reference Ortega-Jimenez, Sapir, Wolf, Variano and Dudley2014; Shyy et al. Reference Shyy, Kang, Chirarattananon, Ravi and Liu2016; Chirarattananon et al. Reference Chirarattananon, Chen, Helbling, Ma, Cheng and Wood2017; Watkins et al. Reference Watkins, Milbank, Loxton and Melbourne2006, Reference Watkins, Burry, Mohamed, Marino, Prudden, Fisher, Kloet, Jakobi and Clothier2020). Nevertheless, several models that have been proposed to describe the arising effects are essentially phenomenological, and only a few seminal contributions addressed the problem from a more fundamental fluid-dynamical perspective (Poudel & Yu Reference Poudel and Yu2018; Poudel, Yu & Hrynuk Reference Poudel, Yu and Hrynuk2021; SureshBabu et al. Reference SureshBabu, Medina, Rockwood, Bryant and Gopalarathnam2021; Wei & Shende Reference Wei and Shende2023). In this regard, it has to be pointed out that many of the available studies on flapping wings in perturbed free stream conditions have some relevant limitations. Firstly, they often rely on two-dimensional numerical simulations of gusts or chaotic flows over airfoils rather than properly three-dimensional and fully developed turbulence over wings (Lian Reference Lian2009; Poudel & Yu Reference Poudel and Yu2018; Poudel et al. Reference Poudel, Yu and Hrynuk2021). Secondly, they typically employ Reynolds-averaged Navier–Stokes or large-eddy simulations approaches, hence not resolving the full range of scales and with uncertainties in the modelling assumptions (Viswanath & Tafti Reference Viswanath and Tafti2010; Jones & Yamaleev Reference Jones and Yamaleev2012). A remarkable exception in this regard is represented by the work of Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019) based on direct numerical simulations (DNS) of realistic insect flight configurations in turbulent environments. Specifically, such studies highlighted for the first time that the mean aerodynamic forces and moments are minimally influenced by free stream turbulence (FST), also for very strong perturbations, whereas their fluctuations increase with the turbulence intensity. Moreover, Engels et al. (Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019) suggested that the integral scale is another key factor in determining the influence of FST on flapping-wing aerodynamics.
A naturally arising question is whether the minimal influence of FST observed for a specific geometry and kinematics by Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019) holds in a more canonical configuration. Following this line of reasoning, in the present work we focus on the impact of FST on the aerodynamics of a flapping wing with a simplified geometry and a prescribed sinusoidal heaving and pitching motion. In particular, we aim at addressing the following questions: (i) How do the characteristic vortical patterns, force (i.e. thrust and lift) and moment generated by the flapping wing get altered by the presence of external flow perturbations? (ii) What are the statistical trends with respect to relevant properties of the impinging FST (i.e. its intensity and characteristic size)? By means of (three-dimensional) DNSs employing a synthetic turbulence approach and an immersed boundary method, we perform a parametric study varying both the intensity and integral length scale of the incoming flow disturbances, providing a statistical characterisation of the main aerodynamic features.
The paper is structured as follows. First, the methodology is introduced in § 2. Then, § 3 presents and discusses the most relevant findings. The section begins with a qualitative description of the instantaneous flow features observed in a representative case. Then, to provide a statistical perspective, the mean vorticity patterns are analysed using phase- and spanwise-averaged flow fields at various stages of the flapping cycle. This is followed by an analysis of the evolution of mean aerodynamic coefficients and load distribution over the cycle to complement the flow field analysis. Subsequently, the study investigates vortex breakup and aerodynamic fluctuations arising from the interaction between the wing and external perturbations, with a particular focus on variations relative to the unperturbed case. Finally, a spectral analysis is conducted to identify the dominant time scales associated with the aerodynamic fluctuations, offering deeper insights into the temporal dynamics of the wing’s response to turbulence. Conclusions are offered in § 4.
Figure 1. Sketch and set-up of the investigated problem: a flapping wing of infinite aspect ratio (in green) is immersed in a free stream of mean velocity
$U$
with external flow perturbations generated in a region of influence (in light blue) using a synthetic turbulence approach. The streamwise, spanwise and transverse directions are denoted as
$x$
,
$y$
and
$z$
, respectively. Here
$h(t)$
and
$\theta (t)$
schematically indicate the heaving and pitching motion, with
${\textit{PP}}$
being the pivot point. Here
$L_x$
,
$L_y$
and
$L_z$
represent the sizes of the computational domain box. The boundaries of the domain are coloured and labelled according to the applied boundary conditions.
2. Methods
We consider a rigid wing of infinite aspect ratio to exclude tip effects, whose cross-section is made of a NACA 0012 airfoil (figure 1). The choice of this kind of airfoil is motivated by the general purpose of investigating the fundamental mechanisms, noting, however, that real-world animals can have remarkably different geometrical features. The Reynolds number, based on the wing chord
$c$
, the free stream (mean) velocity
$U$
and the kinematic viscosity
$\nu$
, is chosen to be
${\mathit{Re}} = U c / \nu = 1000$
. The flapping wing’s motion is set as a combination of heave
$h(t)$
(i.e. translation in the transverse direction) and pitch
$\theta (t)$
(i.e. rotation around an axis along the spanwise direction and passing through a pivot point
${\textit{PP}}$
), both evolving sinusoidally in time,
\begin{align} h(t) &= h_0 \cos (\omega t), \end{align}
\begin{align} \theta (t) &= \theta _0 \cos (\omega t + \phi ). \end{align}
Here,
$h_0$
and
$\theta _0$
are the heaving and pitching amplitude, respectively,
$\omega$
is the angular frequency (which defines the flapping period
$T_{{w}} = 2\pi /\omega$
) and
$\phi$
is the phase lag between pitch and heave. In the case of an unperturbed free stream, this configuration has been investigated in a previous work by Moriche, Flores & García-Villalba (Reference Moriche, Flores and García-Villalba2017), who found some parametric combinations that are optimal in terms of the generated thrust (and eventually also lift) force. In particular, we focus on the ‘A090’ case reported in Moriche et al. (Reference Moriche, Flores and García-Villalba2017), which was found to give the highest thrust coefficient and propulsive efficiency. Hence, the parameters for the wing kinematics are chosen and retained throughout the present study as
$h_0=c$
,
$\theta _0=30^\circ$
,
$\omega = 1.41\, U/c$
,
$\phi =90^\circ$
, with the pivot point located at a quarter of the chord, i.e.
$x_{\textit{PP}}/c = 1/4$
. Note that, in this configuration, the mean pitching angle is zero and, as a result, the mean lift and pitching moment coefficient are negligible (Moriche et al. Reference Moriche, Flores and García-Villalba2017). Also, it is relevant to highlight that the resulting flow field is strictly periodic and two-dimensional (Moriche et al. Reference Moriche, Flores and García-Villalba2016). The two-dimensional nature of the unperturbed free stream case allows direct comparison with the results in the presence of FST in terms of both instantaneous and statistical quantities.
The wing is immersed in an incompressible flow governed by the Navier–Stokes equations,
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u} &= 0, \end{align}
\begin{align} \frac {\partial \boldsymbol {u}}{\partial t} + {(\boldsymbol {u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}} &= - \frac {1}{\rho }\boldsymbol{\nabla }p + \nu \boldsymbol{\nabla} ^2\boldsymbol {u} + \boldsymbol {f}_{\textit{IBM}} + \boldsymbol {f}_{\textit{ST}}, \end{align}
where
$\boldsymbol {u} = (u_x, u_y, u_z)$
and
$p$
are the velocity and pressure field, respectively,
$\rho$
is the fluid density and
$\nu$
is the kinematic viscosity. Boundary conditions are set as follows (see also figure 1): a uniform velocity
$U$
is imposed at the inlet, an advective boundary condition is prescribed at the outlet and periodic conditions apply at the lateral boundaries (i.e. normal to
$y$
or to
$z$
). The infinite span of the wing is modelled by making the wing span equal to
$L_y$
. The immersed boundary method proposed by Uhlmann (Reference Uhlmann2005) is used to enforce the no-slip condition on the surface of the moving wing, which reflects in the presence of the volumetric forcing
$\boldsymbol {f}_{\textit{IBM}}$
in the momentum equation. On the other hand, the free stream perturbations are injected using a synthetic turbulence inflow generator (STIG) (Klein, Sadiki & Janicka Reference Klein, Sadiki and Janicka2003). Specifically, we use a volumetric forcing
$\boldsymbol {f}_{{ST}}$
acting within a ‘region of influence’ located upstream of the wing (light blue box in figure 1), as proposed by Schmidt & Breuer (Reference Schmidt and Breuer2017). In this method, the generation of FST is controlled by means of two input parameters: (i) the turbulence intensity
$\mathit{TI}_0$
and (ii) the integral length scale
$\varLambda _0$
. A characterisation of the generated free stream turbulent perturbations (in the absence of the flapping wing) is reported in Appendix A.
Equations (2.2a ) and (2.2b ) are solved numerically using the in-house code TUCAN (Moriche Reference Moriche2017). It is based on the fractional-step algorithm of Brown, Cortez & Minion (Reference Brown, Cortez and Minion2001), (second-order) central finite differences for spatial discretisation and a low-storage, semi-implicit three-stage Runge–Kutta method for time integration. The code has been extensively validated and employed in previous studies on aerodynamic problems (Arranz, Flores & Garcia-Villalba Reference Arranz, Flores and Garcia-Villalba2020; Moriche et al. Reference Moriche, Sedky, Jones, Flores and García-Villalba2021), also including recent work on the generation and aerodynamic impact of free stream perturbations (Catalán et al. Reference Catalán, Moriche, Flores and Garcia-Villalba2024a , Reference Catalán, Olivieri, García-Villalba and Floresb ). Notably, a version of TUCAN has recently been developed for running on graphical processing units (Guerrero-Hurtado et al. Reference Guerrero-Hurtado, Catalán, Moriche, Gonzalo and Flores2025), which is also employed for the present investigation.
The origin of the reference frame is at the centre of the STIG’s region of influence, with the inlet of the domain lying at
$x/c = -4$
. The wing’s leading edge (LE) is initially (i.e. at
$t=0$
, coinciding with the beginning of the downstroke) set at
$x_{{LE}}/c = 1$
. After testing the influence of varying its size in each direction, the computational domain is set to
$-4 \leqslant x/c \leqslant 4$
,
$-2.5 \leqslant y/c \leqslant 2.5$
and
$-2.5 \leqslant z/c \leqslant 2.5$
(i.e.
$L_x = 8c$
and
$L_y=L_z=5c$
). Concerning the spatial resolution, using a uniform grid spacing
$\varDelta \approx 0.02 c$
in all directions was found to be the best compromise between accuracy and computational cost. Both aspects are further discussed in Appendix B. Lastly, the (constant) time step is chosen in order to have a maximum Courant–Friedrichs–Lewy number of approximately 0.4.
In performing all the simulations, we run the first two flapping cycles (equivalent to approximately two flowthrough times) as a transient that is discarded; then, we run at least other 400 flapping cycles, which are considered for the statistical postprocessing. We perform both a phase average (i.e. between different cycles) and a spanwise average (i.e. between different cross-sections along the
$y$
direction). In fact, the two averages are interchangeable and their combination will be uniquely denoted in the following by
$\langle \boldsymbol{\cdot }\rangle$
. Instead, a temporal average will be denoted by
$\overline {\boldsymbol{\cdot }}$
.
We define the thrust, lift and pitching moment coefficient as follows:
\begin{align} C_T &= \frac {- 2 \, F_x}{\rho U^2 S}, \end{align}
\begin{align} C_L &= \frac {2 \, F_z}{\rho U^2 S}, \end{align}
\begin{align} C_M &= \frac {2 \, M_y^{\textit{PP}}}{\rho U^2 S c}, \end{align}
where
$F_x$
and
$F_z$
are the streamwise and transverse components of the total aerodynamic force acting on the flapping wing, respectively,
$M_y^{\textit{PP}}$
is the aerodynamic pitching moment with respect to the pivot point and
$S = c \, L_y$
is the planform area of the wing. In a similar way, we can introduce the sectional aerodynamic coefficients
$C_t$
,
$C_l$
,
$C_m$
(which will generally be a function of the spanwise coordinate
$y$
, other than time). Finally, we define the propulsive efficiency as
\begin{equation} \eta _{\textit{p}} = \frac {1}{T_{\textit{a}v\textit{g}}} \int _0^{T_{\textit{a}v\textit{g}}} \frac {-F_x U} {(F_z \dot {h} + M_y^{\textit{PP}} \dot {\theta })} \, {\textrm{d}}t, \end{equation}
where
$T_{\textit{a}v\textit{g}} \geqslant 400 T_{{w}}$
.
3. Results and discussion
A series of numerical simulations has been carried out in which the two key quantities characterising the turbulent free stream, i.e. the (nominal) turbulence intensity
$\mathit{TI}_0$
and integral length scale
$\varLambda _0$
, are varied. Overall, the set of considered values for the turbulence intensity is
$\mathit{TI}_0 = 10\,\%, 20\,\%, 30\,\%, 40\,\%$
and 50 % (while fixing
$\varLambda _0/c = 0.6$
), whereas that for the integral length scale is
$\varLambda _0/c = 0.2, 0.4, 0.6, 0.8$
and 1 (while fixing
$\mathit{TI}_0 = 30\,\%$
). Additionally, for the sake of comparison, a simulation for the unperturbed case, i.e.
$\mathit{TI}_0=0$
, has been carried out. In total, 10 different configurations are considered.
3.1. Main features of a flapping wing in a turbulent free stream
To start, we provide a qualitative description of the instantaneous flow obtained in a representative case where the flapping wing is immersed in a turbulent free stream with
$\mathit{TI}=30\,\%$
and
$\varLambda _0/c=0.6$
. Figure 2 shows, for different instants over an entire flapping cycle, the vortical structures detected by means of the well-known
$Q$
-criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988). For comparison, the same snapshots for the unperturbed free stream case are reported alongside. For the unperturbed case, as already anticipated, the instantaneous flow is two-dimensional and the generated vortices are only diffusing while transported downstream (Moriche et al. Reference Moriche, Flores and García-Villalba2016). For the turbulent free stream case, the vortical pattern is much more complex and clearly three-dimensional. Here, vortices are undergoing an instability process, as it is especially evident once sufficiently far from the surface of the wing. Indeed, such instability could be expected given the chaotic character of the free stream perturbations. (The reader is also referred to supplementary videos available online (Olivieri et al. Reference Olivieri, Catalan, Garcia-Villalba and Flores2024), showing the evolution of the flow on a longer time span and for several of the investigated cases.)
Figure 2. Flow visualisations at different instants of a flapping cycle, for (i) a turbulent free stream case with
$\mathit{TI}_0 = 30\,\%$
and
$\varLambda _0/c = 0.6$
and (ii) the unperturbed free stream case (i.e.
$\mathit{TI}_0=0$
). The left- and right-hand sides report four instants of the downstroke (from top to bottom,
$t/T_{{w}} = 0, 1/8, 2/8, 3/8$
) and upstroke (from top to bottom,
$t/T_{{w}} = 4/8, 5/8, 6/8, 7/8$
), respectively. Vortical structures are identified by isosurfaces at
$Q=10 (U/c)^2$
(transparent) and
$30 (U/c)^2$
(opaque), coloured with the spanwise vorticity.
Figure 3. Time histories of the aerodynamic coefficients for a selected turbulent free stream case with
$\mathit{TI}_0 = 30\,\%$
and
$\varLambda _0/c = 0.6$
(in green) and for the unperturbed case (in grey): (a) thrust; (b) lift; (c) pitching moment. In order to highlight the amplitude modulation caused by the turbulent free stream, the minimum and maximum of the unperturbed case solution are marked by horizontal dashed lines.
Looking at the time history of the aerodynamic force coefficients, shown in figure 3, the scenario is consistent with the flow visualisations of figure 2. Note that here both the turbulent free stream case under consideration (green curves) and the unperturbed case (grey curves) are reported. For the latter, the horizontal dashed lines indicate the peak-to-peak amplitude of such periodic signals. Comparing the two cases, it is clear how the turbulent free stream case departs from the fully periodic behaviour of the unperturbed one. Essentially, the most relevant effect that can be highlighted is a modulation in the amplitude of the generated force and moment, particularly noticeable in the thrust and moment coefficient (figure 3 a,c).
Figure 4. Effect of turbulence intensity on the averaged spanwise-vorticity field
$\langle \omega _y \rangle$
, for four stages of the downstroke (from left to right,
$t/T_{{w}} = 0, 1/8, 2/8, 3/8$
). The reported cases are, from top to bottom: (a)
$\mathit{TI}_0 = 10\,\%$
, (b)
$30\,\%$
and (c)
$50\,\%$
, while fixing the same integral length scale
$\varLambda _0/c = 0.6$
. The same colormap is used for all cases, which ranges from
$-30 U/c$
(in blue and anticlockwise) to
$30 U/c$
(in red and clockwise). Additionally, isolines from the unperturbed free stream case are depicted at
$\pm 30 U/c$
(in orange/blue, with white border for visual clarity).
Figure 5. Effect of integral length scale on the averaged spanwise-vorticity field
$\langle \omega _y \rangle$
, for four stages of the downstroke (from left to right,
$t/T_{{w}} = 0, 1/8, 2/8, 3/8$
). The reported cases are, from top to bottom: (a)
$\varLambda _0/c = 0.2$
, (b)
$0.6$
and (c)
$1$
, while fixing the same turbulence intensity
$\mathit{TI}_0 = 30\,\%$
. The same colormap is used for all cases, which ranges from
$-30 U/c$
(blue and anticlockwise) to
$30 U/c$
(red and clockwise). Additionally, isolines from the unperturbed free stream case are depicted at
$\pm 30 U/c$
(in orange/blue, with white border for visual clarity).
3.2. Mean vorticity patterns
To obtain a clearer view on a statistical basis, we analyse the phase- and spanwise-averaged flow fields at different stages of the flapping cycle. For the sake of brevity, only the downstroke phase is shown, since the upstroke presents mirror symmetry. Let us begin with considering the spanwise component of the averaged vorticity field, i.e.
$\langle \omega _y \rangle$
(where
$\boldsymbol{\omega } = \boldsymbol{\nabla }\times \boldsymbol {u}$
denotes the fluid vorticity vector field), shown in figures 4 and 5. The resulting flow fields are depicted for a series of selected cases from the full set of performed simulations. In particular, figure 4 highlights the effect of varying the turbulence intensity
$\mathit{TI}_0$
while figure 5 concerns the variation of the integral length scale
$\varLambda _0$
. For comparison, the unperturbed free stream case is also represented in the figures by means of superimposed vorticity isolines.
Starting with the effect of the turbulence intensity (figure 4) it can be observed that, as
$\mathit{TI}_0$
is increased (from figure 4
a to figure 4
c), the characteristic vortex patterns get progressively smoothed. This feature is clear when looking at the rear part of the airfoil (i.e. approximately from one quarter of the chord) or downstream, and is consistent with the enhanced mixing expected for a turbulent flow. Here, FST strongly modifies the intensity and size of the vortices generated by the wing, which on average appear more diffused. It can be noted, however, that this effect is not so evident in the proximity of the LE, with a much more limited impact of FST on the boundary layer development. In other words, the formation region of the LEV is rather robust to the impinging perturbations. Note that this different outcome between the front and rear part of the airfoil is reflected in the resulting aerodynamic load distribution (as later discussed in § 3.3).
Qualitatively similar observations can be made when looking at the effect of the integral length scale (figure 5), with a promoted diffusion of the generated coherent structures for disturbances of larger size (i.e. from figure 5
a to figure 5
c). However, differently from what observed for the turbulence intensity (figure 4), here it appears that the variation is somehow saturating as
$\varLambda _0$
is increased (i.e. small differences between figures 5
b and 5
c).
A first comparison can be drawn between the phenomenology observed here and that described by Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019) who considered a realistic insect model flapping in a turbulent free stream. Remarkably, in Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019) the generated LEVs are stable, i.e. they never detach from the wing. Instead, in the present case the wing motion is two-dimensional and the resulting flow does not have a net centrifugal component, which is known to be essential to have a stable LEV (Chen et al. Reference Chen, Zhou, Werner, Cheng and Wu2023). Therefore, there is a similarity for what it concerns the robust LEV formation stage; however, the detachment of the vortices and its subsequent enhanced diffusion is not occurring in the case considered by Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019). The aerodynamic impact of these alterations in the vortical patterns are discussed in the next section.
Figure 6. Evolution over the flapping cycle of the mean sectional aerodynamic coefficients in several FST conditions: (a,b) thrust; (c,d) lift; (e,f) pitching moment; (a,c,e) effect of turbulence intensity (red curves, with
$\mathit{TI}_0 = 10, 30, 50 \,\%$
increasing with darkness) while fixing
$\varLambda _0/c=0.6$
; (b,d,f): effect of integral length scale (blue curves, with
$\varLambda _0/c = 0.2, 0.6, 1$
increasing with darkness) while fixing
$\mathit{TI}_0 = 30\,\%$
. Additionally, the unperturbed case (i.e.
$\mathit{TI}_0=0$
) is depicted by the dashed grey curves. Downstroke and upstroke are filled in white and light grey, respectively.
Figure 7. Chordwise distribution of the average pressure coefficient on the upper (solid line) and lower (dashed line) surface at different stages of the downstroke (from top to bottom row): (a,c,e,g) effect of turbulence intensity (red curves, with
$\mathit{TI}_0 = 10, 30, 50 \,\%$
increasing with darkness) while fixing
$\varLambda _0/c=0.6$
; (b,d,f,h) effect of integral length scale (blue curves, with
$\varLambda _0/c = 0.2, 0.6, 1$
increasing with darkness) while fixing
$\mathit{TI}_0 = 30\,\%$
. Additionally, the unperturbed case (i.e.
$\mathit{TI}_0=0$
) is depicted by the grey curves.
3.3. Mean aerodynamic coefficients and load distribution
Along with the inspection of the average flow field, a complementary observation is provided by the evolution over the flapping cycle of the most relevant aerodynamic coefficients, i.e. the thrust, lift and pitching moment. These quantities are shown in figure 6, separating the effect of varying the turbulence intensity (figure 6 a,c,e, red curves) from that of varying the integral length scale (figure 6 b,d, f, blue curves). In each plot, the unperturbed case (dashed grey curves) is also reported as a reference.
Overall, at least for the investigated ranges, the effect of
$\mathit{TI}_0$
appears more pronounced than that of
$\varLambda _0$
. For the mean thrust coefficient
$\langle C_t \rangle$
(figure 6
a,b), whether the turbulence intensity or the integral length scale is increased, we notice that the peaks (at
$t/T_{{w}} \approx 0.25$
and 0.75) are smeared while a mild force increment emerges for
$0 \lesssim t/T_{{w}} \lesssim 0.25$
and
$0.5 \lesssim t/T_{{w}} \lesssim 0.75$
, i.e. in the beginning of the upstroke and downstroke. Similar evidence is found for the mean lift coefficient
$\langle C_l \rangle$
(figure 6
c,d). In the first part of the downstroke, the external perturbations are responsible for moderately increasing the generated lift. However, FST is decreasing the peak of
$\langle C_l \rangle$
at
$t/T_{{w}} \approx 0.25$
. Note that in the upstroke, for the aforementioned symmetry, the same occurs but with opposite signs, with a larger negative lift in
$0.5 \lesssim t/T_{{w}} \lesssim 0.75$
and a smeared negative peak at
$t/T_{{w}} \approx 0.75$
(not shown). Lastly, we look at the time evolution of the mean pitching moment coefficient
$\langle C_m \rangle$
(figure 6
e,f). Compared with the thrust or lift coefficient, in the unperturbed case (grey dashed curves) this quantity has a shape that is richer and more structured. In fact, it appears to be more sensitive to the effect of FST, as it can be appreciated when comparing, e.g. the unperturbed case with the lowest
$\mathit{TI}_0 = 10\,\%$
case, or considering the impact of varying
$\varLambda _0$
(which produces a substantially more contained variation for the thrust and lift).
The variation of the resulting aerodynamic load is consistent with the vorticity patterns previously discussed in § 3.2 and can be illustrated by looking at the (average) pressure coefficient on the wing, reported in figure 7 for several stages of the downstroke (note that, as usual, the vertical axis of the plots is inverted). At first, let us focus on the unperturbed free stream case (grey curves). At the beginning of the downstroke (first column of figures 4 and 5), the (anticlockwise) vortex is in the proximity of the lower surface of the wing and close to its trailing edge. Its influence on the pressure coefficient distribution can be indeed detected in figure 7 (dashed lines in figure 7
a and figure 7
b), essentially limited to the rear part of the wing, i.e.
$0.6 \lesssim s/c \lt 1$
. Due to the orientation of the wing, the contribution of the vortex to the thrust and lift would be essentially detrimental. Since this vortex is quickly damped by the FST (see figures 4 and 5) for
$t/T_{{w}}$
approximately between 0 and 1/8 we consistently observe a mild increase in both
$\langle C_t \rangle$
and
$\langle C_l \rangle$
with
$\mathit{TI}_0$
and
$\varLambda _0$
when compared with the unperturbed case (figure 6). On the other hand, around mid-downstroke (i.e.
$t/T_{{w}}$
approximately between 1/4 and 3/8) the new (clockwise) vortex generated at the LE and developing along the upper surface (see last two columns of figures 4 and 5) is the main responsible for the positive peaks in the thrust and lift coefficient. However, in this stage of the flapping cycle the same effect of FST (i.e. dissipating the vortex more rapidly as it moves downstream) turns out to be detrimental for the aerodynamic performance, eventually smearing the peaks of thrust and lift. This feature is also reflected in the different distributions of the pressure coefficient (solid lines in figure 7
e,f,g,h). Also here, note that the impact of FST is concentrated in the rear part of the wing (i.e.
$0.4 \lesssim s/c \lt 1$
) while the region upstream is significantly less altered by the disturbances. Lastly, it should be noted that the stronger effect of FST on the intensity of the LEV and on the extension of separated regions on the rear part of the airfoil explains the strongest and more complex effect of FST on the moment coefficient
$\langle C_m \rangle$
, mentioned above while discussing figure 6(e,f).
Figure 8. Effect of turbulence intensity on the averaged turbulent kinetic energy field
$\langle \mathit{TKE} \rangle$
, for four stages of the downstroke (from left to right,
$t/T_{{w}} = 0, 1/8, 2/8, 3/8$
). The reported cases are (a)
$\mathit{TI}_0 = 10\,\%$
, (b)
$30\,\%$
and (c)
$50\,\%$
, while fixing the same integral length scale
$\varLambda _0/c = 0.6$
. The same colormap, ranging from
$0$
(blue) to
$U^2$
(red), is used for all cases.
Figure 9. Effect of integral length scale on the averaged turbulent kinetic energy field
$\langle \mathit{TKE} \rangle$
, for four stages of the downstroke (from left to right,
$t/T_{{w}} = 0, 1/8, 2/8, 3/8$
). The reported cases are (a)
$\varLambda _0/c = 0.2$
, (b)
$0.6$
and (c)
$1$
, while fixing the same turbulence intensity
$\mathit{TI}_0 = 30\,\%$
. The same colormap, ranging from
$0$
(blue) to
$U^2$
(red), is used for all cases.
Figure 10. Evolution over the flapping cycle of the standard deviation of the sectional aerodynamic coefficients in several FST conditions: (a,b) thrust; (c,d) lift; (e,f) pitching moment; (a,c,e) effect of turbulence intensity (red curves, with
$\mathit{TI}_0 = 10, 30, 50 \,\%$
increasing with darkness) while fixing
$\varLambda _0/c=0.6$
; (b,d,f): effect of integral length scale (blue curves, with
$\varLambda _0/c = 0.2, 0.6, 1$
increasing with darkness) while fixing
$\mathit{TI}_0 = 30\,\%$
. Downstroke and upstroke are filled in white and light grey, respectively.
3.4. Vortex breakup and aerodynamic fluctuations
After the characterisation of the mean flow and mean aerodynamic coefficients, we analyse the relevant fluctuations that emerge from the interaction between the flapping wing and the external perturbations. To this aim, figures 8 and 9 show the phase- and spanwise-averaged turbulent kinetic energy field
$\langle \mathit{TKE} \rangle$
during the downstroke, for varying turbulence intensity (figure 8) and integral length scale (figure 9). Note that
$\langle \mathit{TKE} \rangle$
is defined by considering the fluctuations with respect to the phase- and spanwise-averaged velocity field. At a glance, it is quite evident that the region where
$\langle \mathit{TKE} \rangle$
is higher (coloured in red) corresponds to the location of the detached LEVs (figures 4 and 5). The most intense turbulent spot is found at
$t/T_{{w}}=0$
(beginning of the downstroke), located below the wing and approximately between the midchord and the trailing edge. Note that at the beginning of the upstroke (
$t/T_{{w}}=4/8$
, not shown) the equivalent turbulent spot is located above the wing. This would follow the sequence shown in the figure. At other stages of the flapping cycle, this region of higher
$\langle \mathit{TKE} \rangle$
changes its intensity, size and position, in agreement with the displacement and dissipation of the separated LEVs. Combined with the observations on the instantaneous fields (figure 2), this suggests that the LEV undergoes an instability and its mean kinetic energy is eventually converted into
$\mathit{TKE}$
. Also, note that the region where such highly turbulent spots are produced becomes wider as
$\mathit{TI}_0$
or
$\varLambda _0$
is increased. However, the trend while increasing the integral length scale appears to saturate, as it was observed for the vorticity field (figure 5).
Next, we focus on how FST affects the fluctuation of the aerodynamic coefficients. Figure 10 shows the evolution over the flapping cycle of the standard deviation of the thrust, lift and pitching moment coefficient. For the thrust coefficient (figure 10
a,b), as the turbulence intensity or integral length scale is increased, a clear pattern emerges towards an approximately triangular profile with a peak at midstroke (i.e.
$t/T_{{w}} = 0.25$
and 0.75). Recalling the shape of the mean thrust coefficient (figure 6
a,b), this is consistent with the idea that the vortex breakup provides the major contribution in the resulting
$\mathit{TKE}$
, and in turn in the aerodynamic force and moment fluctuations. In fact, the same trend is found for the lift coefficient (figure 10
c,d), although this quantity shows a more complex shape. For the pitching moment coefficient (figure 10
e,f), the evolution of its standard deviation at sufficiently large
$\mathit{TI}_0$
or
$\varLambda _0$
approximately resembles a sawtooth profile. Overall, note that the fluctuation in the force and moment is governed not only by the intensity of the vortex breakup, but also by its degree of proximity to the wing surface. For this reason the peaks of the standard deviations shown in figure 10 are found at
$t/T_{{w}} \approx 0.25$
and
$0.75$
(or slightly delayed for the pitching moment).
3.5. Fluctuations with respect to the unperturbed case
To further understand how the aerodynamic coefficients are affected by the presence of FST, it can be useful to introduce their (instantaneous) variation with respect to the unperturbed case (i.e.
$\mathit{TI}_0=0$
),
\begin{align} \Delta C_T &= C_T - C_T|_{\mathit{TI}_0=0}, \end{align}
\begin{align} \Delta C_L &= C_L - C_L|_{\mathit{TI}_0=0}, \end{align}
\begin{align} \Delta C_M &= C_M - C_M|_{\mathit{TI}_0=0}. \end{align}
Figure 11. Time histories of the instantaneous variation with respect to the unperturbed case of the aerodynamic coefficients, for the selected turbulent free stream case with
$\mathit{TI}_0 = 30\,\%$
and
$\varLambda _0/c = 0.6$
: (a) thrust; (b) lift; (c) pitching moment.
A representative sample of these quantities is shown in figure 11 for the same case already considered in § 3.1 (i.e.
$\mathit{TI}_0=30\,\%$
and
$\varLambda _0 = 0.6 c$
). At a glance, it can be noted the erratic nature of these signals, suggesting a wide frequency spectrum (later shown in figure 15). Also, it is apparent that
$\Delta C_L$
and
$\Delta C_M$
have negligible mean values, although it is less clear for
$\Delta C_T$
. Indeed,
$\overline {\Delta C_T}$
is shown for all the simulated cases in figure 12, focusing on its dependency with both the turbulence intensity and integral length scale. For the unperturbed case, the net (i.e. time-averaged) aerodynamic coefficients are found to be
$\overline {C_T}|_{\mathit{TI}_0=0} \approx 0.92$
,
$\overline {C_L}|_{\mathit{TI}_0=0} \approx 0.00$
and
$\overline {C_M}|_{\mathit{TI}_0=0} \approx 0.00$
. For the thrust coefficient, since its value in the unperturbed case is around unity, it can be estimated that the relative variation due to FST is in any case no larger than approximately
$5\,\%$
. In other words, FST does not produce a considerable net variation in terms of the mean thrust. On the other hand, given the symmetry of the kinematics,
$\overline {\Delta C_L} = \overline {\Delta C_M} = 0$
(not shown). Furthermore, negligible variations are also found for the propulsive efficiency
$\eta _{{p}}$
, which always lies around 0.35, as in the unperturbed case (Moriche et al. Reference Moriche, Flores and García-Villalba2017).
Figure 12. Mean (i.e. time-averaged) variation of the thrust coefficient with respect to the unperturbed case, as a function of (a) turbulence intensity (while fixing
$\varLambda _0/c=0.6$
) and (b) integral length scale (while fixing
$\mathit{TI}_0 = 30\,\%$
). The bars indicate the standard error of the sample mean.
Figure 13. Standard deviation of the variation of the aerodynamic coefficients with respect to the unperturbed case: (a,b) thrust; (c,d) lift; (e,f) pitching moment; (a,c,e) effect of turbulence intensity (while fixing
$\varLambda _0/c=0.6$
); (b,d,f): effect of integral length scale (while fixing
$\mathit{TI}_0 = 30\,\%$
).
Figure 13 reports the (time-averaged) standard deviations of
$\Delta C_T$
,
$\Delta C_L$
and
$\Delta C_M$
, which are found to follow some monotonic trends in both
$\mathit{TI}_0$
and
$\varLambda _0$
. For all the aerodynamic coefficients, the standard deviation appears to grow linearly with the turbulence intensity (figure 13
a,c,e). Conversely, it tends to saturate with the integral length scale, in a more or less accentuated way depending on the aerodynamic coefficient considered (figure 13
b,d,f). Note that the quantities shown in figure 13 are, in general, slightly different from the temporal average of the standard deviations that are reported in figure 10. Nevertheless, this subtle distinction does not affect the resulting trends.
Figure 14. Probability density functions (PDFs) of the variation of the aerodynamic coefficients with respect to the unperturbed case: (a,b) thrust; (c,d) lift; (e,f) pitching moment; (a,c,e) effect of turbulence intensity (red curves, with
$\mathit{TI}_0 = 10, 20, 30, 40, 50 \,\%$
increasing with darkness) while fixing
$\varLambda _0/c=0.6$
; (b,d,f) effect of integral length scale (blue curves, with
$\varLambda _0/c = 0.2, 0.4, 0.6, 0.8, 1$
increasing with darkness) while fixing
$\mathit{TI}_0 = 30\,\%$
.
The PDF of the aerodynamic coefficient fluctuations is shown in figure 14. Increasing the turbulence intensity (figure 14
a,c,e) or the integral length scale (figure 14
b,d,f) the distributions get wider, as expected from the trends of the standard deviation (figure 13). For the effect of the turbulence intensity, in particular, the curves can be easily collapsed by normalising with a linear scaling in
$\mathit{TI}_0$
(not shown). In passing, it can be noted a certain difference between the thrust coefficient (figure 14
a,b) exhibiting a mild degree of non-Gaussianity (i.e. slightly negative skewness and kurtosis of the PDF typically larger than 4), and the other (lift and moment) aerodynamic coefficients resembling more a normal distribution (i.e. no skewness and a kurtosis of around 3).
It can be noted that these findings overall agree with, and generalise those of Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019) who reported a negligible sensitivity of the mean aerodynamic coefficients, a linear trend of their fluctuations with respect to the turbulence intensity and a stronger effect for sufficiently large integral length scales. Remarkably, the statistical aerodynamical response looks similar also when considering a general heaving-and-pitching motion as in the present case, where the LEVs are not stable.
Figure 15. Premultiplied power spectral densities (PSDs) of the variation of the aerodynamic coefficients with respect to the unperturbed case: (a,b) thrust; (c,d) lift; (e,f) pitching moment; (a,c,e) effect of turbulence intensity (red curves, with
$\mathit{TI}_0 = 10, 30, 50 \,\%$
increasing with darkness) while fixing
$\varLambda _0/c=0.6$
; (b,d,f) effect of integral length scale (blue curves, with
$\varLambda _0/c = 0.2, 0.6, 1$
increasing with darkness) while fixing
$\mathit{TI}_0 = 30\,\%$
. In (a) and (b), the dashed grey line indicates twice the flapping motion’s frequency
$2 f_{{w}} = 2/T_{{w}}$
, whereas in the other panels it indicates the flapping motion frequency
$f_{{w}} = 1/T_{{w}}$
; the dash–dotted grey lines indicate its higher harmonics at which secondary peaks are found. The coloured dotted lines indicate the characteristic turbulent frequency
$f_0 = 1/\tau _0$
(see also figure 16).
3.6. Spectral analysis
As the last step of this study, we look at the characteristic time scales that can be identified in the aerodynamic fluctuations under consideration. Figure 15 shows the PSD of the force and moment variation for the same cases analysed in §§ 3.3 and 3.4, again separating the effect of the turbulence intensity and integral length scale (note that the PSD are premultiplied by the frequency to better appreciate the most energetic frequencies on a log–linear plot). In the plots, we have also indicated with a dashed vertical line the main frequency related to the flapping motion,
$f = n f_{{w}} = n/T_{{w}}$
, with
$n=2$
for
$\Delta C_T$
and
$n=1$
for
$\Delta C_L$
or
$\Delta C_M$
. Note that this is essentially the characteristic frequency of the unperturbed case solution, which has been subtracted in the definitions of
$\overline {C_T}$
,
$\overline {C_L}$
and
$\overline {C_M}$
(see (3.1)). Nevertheless, from the PSD of the thrust variation for different turbulence intensities (figure 15
a) it clearly appears that the peaks are still in correspondence of
$2 f_{{w}}$
. Similarly, the peaks of both
$\Delta C_L$
and
$\Delta C_M$
(figure 15
c,d,e,f) are in the proximity of
$f_{{w}}$
. Additionally, it can be noted the presence of secondary peaks at higher harmonics of
$f_{{w}}$
(dash–dotted lines), particularly evident for the pitching moment.
Figure 16. Premultiplied PSDs of the free stream velocity fluctuations (normalised with the variance): (a) effect of turbulence intensity (red curves, with
$\mathit{TI}_0 = 10, 30, 50 \,\%$
increasing with darkness) while fixing
$\varLambda _0/c=0.6$
; (b) effect of integral length scale (blue curves, with
$\varLambda _0/c = 0.2, 0.6, 1$
increasing with darkness) while fixing
$\mathit{TI}_0 = 30\,\%$
. For each case, the dotted coloured lines indicate the estimated frequency
$f_0 = 1/\tau _0$
.
When looking at the effect of the integral length scale on the thrust (figure 15
b), focusing on the largest value of
$\varLambda _0$
considered (darkest blue curve), however, we have the emergence of a significant widebanded region at lower frequencies (centred around
$fc/U \approx 0.3$
). Conversely, for the lift and moment a similar trend is not evident. Such evidence might be explained by the combination of different time scales. Indeed, along with that of the flapping motion, we can introduce another time scale that is related to the FST. In particular, we consider a characteristic time scale based on the integral length scale and the free stream mean velocity,
$\tau _0 \sim \varLambda _0 / U$
. To better explain its meaning, figure 16 shows the PSD of the turbulent velocity fluctuations
$u_{\mathrm{ST}}$
that are injected within the region of influence upstream of the wing. In both figures 15 and 16, the frequency
$f_0 = 1/\tau _0$
is indicated by vertical dotted lines, coloured according to the corresponding case (note that in figures 15
a,c,e and 16
a,c,e
$f_0$
is the same for all cases, being
$\varLambda _0$
the same). As expected, in figure 16 the peaks are found in the proximity of this characteristic frequency.
Considering again the thrust variation of figure 15(b), it can be noticed that the widebanded region for the largest
$\varLambda _0/c=0.6$
is somehow located between
$f_0$
and
$f_{{w}}$
. The fact that the peaks in the spectra of
$\Delta C_T$
,
$\Delta C_L$
and
$\Delta C_M$
do not appear on
$f_0$
suggests that the effect of FST on the aerodynamic forces is not direct, but mediated through the LEV dynamics. Although this observation is limited to the considered range of
$\mathit{TI}_0$
and
$\varLambda _0$
, it looks consistent with the results presented above. Nevertheless, a deeper analysis and an extension of the investigated parametric space would be needed in order to better conclude on this point, which we leave open for future work.
4. Conclusions
This work has been motivated by improving our understanding on how biologically inspired flapping-wing devices (e.g. MAVs or AUVs) can perform in perturbed or turbulent flow environments. Tackling the issue from a fundamental viewpoint, we have focused on the impact of FST on the aerodynamics of a flapping wing of infinite planform aspect ratio with an imposed heaving and pitching motion that, in unperturbed free stream conditions, generates thrust with an optimal propulsive efficiency (Moriche et al. Reference Moriche, Flores and García-Villalba2017). We have presented results from DNS at a relatively low Reynolds number of 1000 employing a STIG to introduce free stream perturbations of controlled intensity and integral length scale. These two key parameters have been spanned over a representative range to highlight some peculiar features and trends. Furthermore, the simulations were conducted over long time intervals – typically exceeding 400 flapping cycles – to ensure statistical convergence, with durations generally longer than those in previous studies (Engels et al. Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019).
As expected, introducing the turbulent perturbations causes the qualitative alteration of the vortical patterns, which transition from the two-dimensional and periodic character of the reference unperturbed case to a fully three-dimensional and chaotic one in which the vortices are unstable. This is accompanied by a significant modulation in the amplitude of the aerodynamic coefficients. From the phase- and spanwise-averaged flow fields it can be deduced that FST is able to effectively enhance the dissipation of the generated LEVs in the rear part of the wing or downstream. However, the boundary layer in the proximity of the LE is found to be particularly robust to the action of the external disturbances, even in the most severe turbulence conditions that were tested. The shape of the phase- and spanwise-averaged aerodynamic coefficients can be altered by the impinging perturbations, resembling what is observed in previous studies with flapping motions involving centrifugal effects (Fisher et al. Reference Fisher, Ravi, Watkins, Watmuff, Wang, Liu and Petersen2016; Engels et al. Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019). Yet, FST does not induce an appreciable variation of the net thrust. Of particular interest is the comparison with the results of Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019) for a model bumblebee undergoing a proper flapping motion generating stable LEVs, since the present results suggest that the stability of the LEV is not a necessary condition in order to have a minimal influence of FST on the aerodynamic performance, at least in a statistical sense.
The vortex breakup (triggered by FST) appears to be the main source of flow disturbances that the wing is most likely to experience in terms of fluctuations of the aerodynamic force and moment, with a clear variability during the flapping cycle. The standard deviations of the aerodynamic coefficient fluctuations (with respect to the unperturbed case) appear to scale linearly with the turbulence intensity, as also reported by Engels et al. (Reference Engels, Kolomenskiy, Schneider, Lehmann and Sesterhenn2016, Reference Engels, Kolomenskiy, Schneider, Farge, Lehmann and Sesterhenn2019), while they tend to saturate with the integral length scale. Overall, the aerodynamic fluctuations appear to be symmetrically distributed and generally close to Gaussian, only with a certain departure found for the thrust coefficient. Lastly, from a spectral analysis it can be outlined that the dominant frequency in the thrust coefficient fluctuations is given by a non-trivial combination of the flapping frequency and the characteristic time scale of the incoming flow disturbances.
Our findings provide statistical information on the interaction between a turbulent free stream and an isolated flapping wing with fully prescribed kinematics. Directions for future work include: (i) assessing the influence of airfoil and planform geometry on the interaction between FST and aerodynamic load generation; (ii) releasing certain degrees of freedom to model aspects of flight or swimming dynamics; (iii) testing active or passive control techniques to mitigate the aerodynamic disturbances induced by FST; (iv) achieving higher Reynolds number flows and wider scale separation.
Acknowledgements
The authors thankfully acknowledge the computer resources at Picasso and the technical support provided by the Supercomputing and Bioinnovation Center (SCBI) of the University of Malaga (RES-IM-2023-1-0021 and RES-IM-2023-2-0017). S.O. acknowledges CINECA and INFN for the availability of high-performance computing resources and support.
Funding
This work is supported by grant PID2022-142135NA-I00 by MCIN/AEI/10.13039/501100011033 and grants FJC2021-047652-I and TED2021-131282B-I00 by MCIN/AEI/10.13039/501100011033 and European Union NextGenerationEU/PRTR.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Characterisation of FST generation
Here, we provide some additional information on the turbulent free stream generation. To this aim, we consider the situation in which the flapping wing is absent, in order to isolate and measure the properties of the free stream perturbations. Note that the same problem has been recently surveyed in Catalán et al. (Reference Catalán, Olivieri, García-Villalba and Flores2024b ).
Figure 17. Streamwise evolution of the effective turbulence intensity (a,b) and large-scale anisotropy ratio (c,d) for cases at different (nominal) turbulence intensity
$\mathit{TI}_0 = 10, 30, 50 \,\%$
((a,c) red curves) and integral length scale
$\varLambda _0/c = 0.2, 0.6, 1$
((b,d) blue curves). Both
$\mathit{TI}_0$
and
$\varLambda _0$
increase with darkness. The shaded area indicates the extension of the STIG’s influence region in the baseline case.
Figure 17 shows how the (effective) turbulence intensity
$\mathit{TI}=\sqrt {2/3\,\mathit{TKE}} / U$
and large-scale anisotropy ratio
$A = 2 u_{\textit{rms}} / (v_{\textit{rms}} + w_{\textit{rms}})$
vary along the streamwise direction, for various combinations of the nominal turbulence intensity
$\mathit{TI}_0$
(figure 17
a,c) and integral length scale
$\varLambda _0$
(figure 17
b,d). Looking at the turbulence intensity (figure 17
a,b), it can be noticed how the evolution of the flow is somehow different depending on the value of
$\mathit{TI}_0$
and
$\varLambda _0$
. Specifically, a noticeable decay of turbulence is found when
$\mathit{TI}_0$
is the largest or
$\varLambda _0$
is the smallest. For the purpose of the present study, however, the flapping wing is always set with its LE at
$x/c = 1$
, so that the properties of the turbulent fluctuations are always close to the nominal ones. In terms of flow isotropy (figure 17
c,d) we observe that, once downstream of the influence region (shaded area), the large-scale anisotropy ratio
$A$
approximately follows a plateau at values close to unity. A mild deviation from
$A=1$
(i.e. streamwise and transverse fluctuations with same intensity) is observed for cases with larger
$\varLambda _0$
, which can be ascribed to finite-size effects and could be further limited by increasing the domain extension.
Next, focusing on the baseline case, we examine the relevant length scales associated with the generated fluctuations, as reported in figure 18(a). The plot shows the streamwise evolution of (i) the integral longitudinal length scale
$\varLambda _f$
(solid), (ii) the longitudinal Taylor length scale
$\lambda _f$
(dashed) and (iii) Kolmogorov length scale
$\eta$
(dotted line). Following the standard approach and definitions, the first two quantities are obtained from the velocity autocorrelation coefficient whereas the latter is obtained by the turbulent energy dissipation rate (Pope Reference Pope2000). While the Taylor length scale undergoes some variations,
$\varLambda _f$
remains almost constant and close to the nominal value
$\varLambda _0$
in the whole range of interest. Similarly, we observe a slow decrease for
$\eta$
. Finally, figure 18(a) shows the resulting energy spectrum at various streamwise locations. It is clear the tendency to fill the spectrum at higher wavenumbers while moving downstream, consistent with the full development of the turbulent flow (Catalán et al. Reference Catalán, Olivieri, García-Villalba and Flores2024b
).
Figure 18. Streamwise evolution of characteristic length scales and energy distribution for the baseline case (
$\mathit{TI}_0 = 30\,\%$
and
$\varLambda _0/c=0.6$
): (a) longitudinal integral length scale
$\varLambda _f$
(solid), longitudinal Taylor length scale
$\lambda _f$
(dashed) and Kolmogorov length scale
$\eta$
(dotted line) as a function of the streamwise coordinate (the shaded area indicates the STIG’s influence region); (b) energy spectrum function evaluated at
$x/c \approx 1$
, 2, 3 and 4 (corresponding to curves of increasing brightness, respectively). In (b), the inset shows the energy spectrum premultiplied by the wavenumber in a log–linear scale.
Figure 19. Sensitivity of the results with respect to the size of the computational domain. The series of cases where
$\varLambda _0$
is varied (as shown in figures 6
b,d,f and 10
b,d,f) is here enriched by additional simulations at
$\varLambda _0/c = 1$
(denoted by different symbols, as explained in the main text) where the size of the domain is varied. Here (a,c,e) and (b,d,f) report the evolution of the mean and standard deviation of the sectional aerodynamic coefficients over the flapping cycle, respectively. (a,b) thrust; (c,d) lift; (e,f) pitching moment. Downstroke and upstroke are filled in white and light grey, respectively.
Appendix B. Convergence studies
This section provides additional information on the choice of the numerical settings (i.e. extent of the computational domain and grid resolution) adopted for the present investigation.
To test the sensitivity of the solution to the domain size, we have focused on the case at
$\mathit{TI}_0 = 30\,\%$
and
$\varLambda _0/c = 1$
, i.e. that with the largest integral length scale considered in this study, for which finite-size effects are expected to be most critical. Departing from the baseline setting (i.e.
$-4 \leqslant x/c \leqslant 4$
,
$-2.5 \leqslant y/c \leqslant 2.5$
and
$-2.5 \leqslant z/c \leqslant 2.5$
), we have repeated the simulation varying the domain extension. The results are presented in figure 19, where these additional cases are denoted by symbols and reported along with the other cases where the integral length scale is varied. The tests have been carried out as follows. First, we have assessed the sensitivity with respect to the streamwise direction (while keeping the same spanwise and transverse extension), extending the domain downstream (i.e.
$-4 \leqslant x/c \leqslant 12$
; the case is denoted by circles in figure 19) or upstream (i.e.
$-8 \leqslant x/c \leqslant 4$
; squares in figure 19) to better accommodate the wake generated by the wing or the region of influence of the synthetic turbulence generator, respectively. Looking at the resulting aerodynamic coefficients, the effect of such variations are found to be minimal. Then, keeping the initial streamwise extension
$-4 \leqslant x/c \leqslant 4$
, we have tested the effect of doubling the spanwise extension (i.e.
$-5 \leqslant y/c \leqslant 5$
; triangles in figure 19), observing negligible variations. Last, we have looked at the impact of doubling the transverse extension (i.e.
$-5 \leqslant y/c \leqslant 5$
; diamonds in figure 19). In this case, a non-negligible, yet contained, effect can be appreciated. However, these variations turn out to be small when compared with those caused by the changes in
$\varLambda _0$
considered in the study.
Moving to the assessment of the adequate grid resolution
$\varDelta$
, it is relevant to decouple the requirements posed by the generation of the turbulent free stream, on one hand, from those of solving the flow around the flapping wing in the unperturbed free stream case, on the other one. Both issues have been separately addressed, at the same Reynolds number of this work, in previous works by Catalán et al. (Reference Catalán, Olivieri, García-Villalba and Flores2024b
) (for the turbulent free stream) and Moriche et al. (Reference Moriche, Flores and García-Villalba2016, Reference Moriche, Flores and García-Villalba2017) (for the flapping wing in unperturbed free stream). In Catalán et al. (Reference Catalán, Olivieri, García-Villalba and Flores2024b
), it has been shown that a uniform grid spacing
$\varDelta \approx 0.04 c$
is sufficient to accurately resolve a turbulent flow with similar properties as in the present case. On the other hand, Moriche et al. (Reference Moriche, Flores and García-Villalba2016, Reference Moriche, Flores and García-Villalba2017) showed that a finer resolution of
$\varDelta \approx 0.01 c$
is recommended to accurately simulate a flapping wing in a uniform free stream. Therefore, the latter turns out to give the strictest requirement on the adequate
$\varDelta$
. Note, however, that for the present study, unlike the case of a uniform free stream, there is a need to achieve a satisfactory statistical convergence by running
${O}(10^2)$
flapping cycles. Consequently, after testing we have chosen a grid spacing of
$\varDelta \approx 0.02 c$
which, although being slightly under-resolved compared with that used by Moriche et al. (Reference Moriche, Flores and García-Villalba2016), represents the most feasible compromise for carrying out the present statistical analysis.
Open access
