1. Introduction
Flapping wings are used to generate lift and thrust for birds and insects. In contrast to fixed wings used in traditional air vehicles, which are more likely to stall in the low Reynolds number (Re) regime (
$Re < 10^5$) due to the low resistance to the adverse pressure gradient for the laminar boundary layer, flapping wings could produce a higher lift efficiently (Maxworthy Reference Maxworthy1981; Wilga & Lauder Reference Wilga and Lauder1999; Sane Reference Sane2003; Fish & Lauder Reference Fish and Lauder2006; Platzer et al. Reference Platzer, Jones, Young and Lai2008; Shyy et al. Reference Shyy, Lian, Tang, Viieru and Liu2008; Pesavento & Wang Reference Pesavento and Wang2009; Liu, Wang & Liu Reference Liu, Wang and Liu2024). Due to their superior force production in the low Reynolds number regime, bio-inspired flapping wings have been used in micro-air vehicles (MAVs), which are typically operated within the same low Reynolds number regime as experienced by natural flyers (Jones, Duggan & Platzer Reference Jones, Duggan and Platzer2001; Pines & Bohorquez Reference Pines and Bohorquez2006; Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010). Therefore, understanding force production mechanisms of flapping wings and promoting their applications in designing bio-inspired MAVs are of significant interest (Ellington Reference Ellington1999).
Bio-inspired flapping wings generate high lift through a few sophisticated unsteady mechanisms, of which examples are the complicated vortex system, rapid pitching, wake capture, passive deformation and synchronisations of kinematics, wing flow and powers (Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010; Eldredge & Jones Reference Eldredge and Jones2019; Huang et al. Reference Huang, Bhat, Yeo, Young, Lai, Tian and Ravi2023; Liu et al. Reference Liu, Wang and Liu2024). These mechanisms are associated with sophisticated kinematics such as stroke, pitching, deviation and other components. Stroke is the flapping motion of the spanwise axis of the wing; pitching is the rotation of the wing about the spanwise axis; and deviation is the movement of the spanwise axis away from the stroke plane defined by the wing root and the wing tip of the maximum and minimum positions. The stroke motion can be decomposed into translation and rotation, the former is also called heaving motion and the latter is known as revolving motion if the rotational direction is not changing. The vortex system consists of leading-edge vortex (LEV), trailing-edge vortex (TEV) and tip vortex (TV), as shown in figure 1. Among these unsteady mechanisms, the LEV has been discovered to contribute to force production significantly (Maxworthy Reference Maxworthy1979; Ellington et al. Reference Ellington, Van Den Berg, Willmott and Thomas1996; Dickinson, Lehmann & Sane Reference Dickinson, Lehmann and Sane1999; Videler, Stamhuis & Povel Reference Videler, Stamhuis and Povel2004; Muijres et al. Reference Muijres, Johansson, Barfield, Wolf, Spedding and Hedenström2008). The LEV has also been found in the dynamic stall of helicopter blades and aerofoils undergoing a transient change in the effective angle of attack (McCroskey Reference McCroskey1982; Mulleners & Raffel Reference Mulleners and Raffel2012; Choudhry et al. Reference Choudhry, Leknys, Arjomandi and Kelso2014; Sangwan, Sengupta & Suchandra Reference Sangwan, Sengupta and Suchandra2017). Due to its importance, great effort has been made to study the LEV of flapping wings/foils during the last three decades (Shyy & Liu Reference Shyy and Liu2007; Eldredge & Jones Reference Eldredge and Jones2019; Liu et al. Reference Liu, Wang and Liu2024). For example, flow visualisation experiments conducted by van den Berg & Ellington (Reference van den Berg and Ellington1997) showed the LEV formation on a hawkmoth wing model, and exposed that the LEV in insect hovering flight can be generated by either the wing rotation or stroke. By using a simple linear stability analysis, Triantafyllou, Triantafyllou & Grosenbaugh (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993) found that the optimal Strouhal number for efficient propulsion is in the range of 
$0.25$–
$0.35$, which is the frequency range of the maximum amplification of the convectively unstable wake. This finding was confirmed by an experimental study (Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998), and observed in natural flyers (Taylor, Nudds & Thomas Reference Taylor, Nudds and Thomas2003). In addition, the conditions of high efficiency are associated with the formation of moderately strong LEVs, which interact with the TEVs, and consequently result in a reverse Karmán street (Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998). Actually, the LEV evolution is mainly determined by the reduced frequency and the effective angle of attack rather than the Strouhal number, as shown by an experimental study on the force generation and the LEV evolution of pitching and heaving aerofoils (Baik et al. Reference Baik, Bernal, Granlund and Ol2012). By considering an accelerated flat-plate aerofoil experimentally, Pitt Ford & Babinsky (Reference Pitt Ford and Babinsky2013) found that the circulation build-up is contained in the LEV and the TEV instead of the vortex sheet bounded on the wing in the conventional thin-aerofoil theory. The detachment of the LEV could occur during its evolution due to its continuing growth. The leading-edge shear layer and the secondary vortex play an important role in the growth and detachment of the LEV, as shown by Rival et al. (Reference Rival, Kriegseis, Schaub, Widmann and Tropea2014), Widmann & Tropea (Reference Widmann and Tropea2015) and Li et al. (Reference Li, Feng, Kissing, Tropea and Wang2020), who considered the length scales and other dimensionless parameters. At very low or very high amplitude-based Strouhal numbers, the three-dimensional transition of the LEV would occur, leading to the strong three-dimensional effects (Zurman-Nasution, Ganapathisubramani & Weymouth Reference Zurman-Nasution, Ganapathisubramani and Weymouth2020). A linear stability analysis was conducted by Gao et al. (Reference Gao, Cantwell, Son and Sherwin2023) to understand the spanwise instabilities of the LEV and their effects on the force of a heaving aerofoil. The spanwise flow is able to stabilise the LEV for natural flyers (Maxworthy Reference Maxworthy1979; Ellington et al. Reference Ellington, Van Den Berg, Willmott and Thomas1996; Liu et al. Reference Liu, Ellington, Kawachi, Van Den Berg and Willmott1998; Jardin & David Reference Jardin and David2014), which was further explained by Xia & Mohseni (Reference Xia and Mohseni2023) with a finite-area sink model.
Figure 1. Definition sketch of vortices around a flapping wing.
The three-dimensional effects and the TVs are the major features triggered by low-aspect-ratio flapping wings (Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010). Many aerial and aquatic animals have wings or fins of low aspect ratios that allow them to fly or swim in tight spaces at low Reynolds numbers (Azuma Reference Azuma2012). In contrast to the large-aspect-ratio wings or foils used in most experiments and two-dimensional numerical simulations, the low-aspect-ratio wings/foils are more suitable for the MAVs due to the size requirements (Pines & Bohorquez Reference Pines and Bohorquez2006). Therefore, studying force production and the flow structures of low-aspect-ratio flapping wings is essential for advancing flapping-wing MAV design. There have been many studies on low-aspect-ratio wings, uncovering a few important findings. The TVs contribute to the total force exerted on the wing in three major ways: direct vortex force, stabilising the LEV and enhancing downwash (Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010; Sengupta Reference Sengupta2015). For example, the instantaneous dynamics and the history of the hovering wing dominate the aerodynamic force (Liu et al. Reference Liu, Ellington, Kawachi, Van Den Berg and Willmott1998). The axial flow detected within the LEV can stabilise the LEV (Liu et al. Reference Liu, Ellington, Kawachi, Van Den Berg and Willmott1998; Shyy et al. Reference Shyy, Kang, Chirarattananon, Ravi and Liu2016). The peak in the thrust forces is higher when the wing rotation is ahead of the stroke reversal compared with other cases (Ramamurti & Sandberg Reference Ramamurti and Sandberg2002). The propulsive performance of low-aspect-ratio flapping foils is determined by their vortical wake, which is characterised by two sets of interconnected vortex loops (Dong, Mittal & Najjar Reference Dong, Mittal and Najjar2006). The wake dynamics is dominated by the aspect ratios due to the interactions with the TV, as shown by a numerical study focusing on the effects of aspect ratio, angle of attack and geometry (Taira & Colonius Reference Taira and Colonius2009b). In contrast to the fixed finite wing, the TV could not only generate the lift for flapping wings but also delay the shedding of the LEV (Shyy et al. Reference Shyy, Trizila, Kang and Aono2009). In addition, as an important contribution to the forces generated by the low-aspect-ratio pitching plates, the wake circulation shed in one stroke scales with the kinematic parameters and the geometric parameters (Buchholz & Smits Reference Buchholz and Smits2008; Buchholz, Green & Smits Reference Buchholz, Green and Smits2011). Garmann & Visbal (Reference Garmann and Visbal2014) performed direct numerical simulations to study the vortex dynamics of low-aspect-ratio revolving wings and found that the Coriolis force does not contribute to the attachment of the LEV. The vorticity transport for heaving wings with aspect ratio (AR) of 
$AR = 2$ and 4 was analysed by Akkala & Buchholz (Reference Akkala and Buchholz2017) who found that the opposite-sign vortex between the LEV and the boundary layer regulates the growth of the LEV. Menon, Kumar & Mittal (Reference Menon, Kumar and Mittal2022) found that, in the wake region, the streamwise vortices contribute to positive lift production, in contrast to the spanwise vortices producing a negative lift on the wing. As shown by the above discussion, the evolution of flow structures for a low-aspect-ratio wing is much more complex than that of two-dimensional or large-aspect-ratio wings due to the three-dimensional effects and the TVs.
As a new type of air vehicle, bio-inspired MAVs powered by low-aspect-ratio wings encounter many unique challenges due to their small size, low weight and complex flight environments (Pines & Bohorquez Reference Pines and Bohorquez2006; Moulin & Karpel Reference Moulin and Karpel2007). Although the unsteady vortices generated by flapping wings can bring many benefits, such as a high transient aerodynamic force and energy harvesting, they can also bring unfavourable effects, such as flutter, buffeting and dynamic response (Bisplinghoff, Ashley & Halfman Reference Bisplinghoff, Ashley and Halfman2013). Therefore, it is desired to develop novel control strategies to maintain manoeuvrability and aerodynamic efficiency under different flying conditions, which natural flyers have been doing well (Othman et al. Reference Othman, Zekry, Saro-Cortes, Lee and Wissa2023). Inspired by the natural flyers, a few control strategies have been developed (Ho et al. Reference Ho, Nassef, Pornsinsirirak, Tai and Ho2003; Fish & Lauder Reference Fish and Lauder2006; Jones & Platzer Reference Jones and Platzer2010; Rose, Natarajan & Gopinathan Reference Rose, Natarajan and Gopinathan2021). Flow control on traditional wings/aerofoils to mitigate or delay the formation of the dynamic stall vortex by suppressing flow separation has been widely studied (McCroskey Reference McCroskey1982; Carr Reference Carr1988; Greenblatt & Wygnanski Reference Greenblatt and Wygnanski2000), including suction, blowing and plasma controls (see e.g. Huang et al. Reference Huang, Huang, LeBeau and Hauser2004; Corke & Post Reference Corke and Post2005; Chng et al. Reference Chng, Rachman, Tsai and Zha2009; Chen, Seele & Wygnanski Reference Chen, Seele and Wygnanski2013; Fahland et al. Reference Fahland, Stroh, Frohnapfel, Atzori, Vinuesa, Schlatter and Gatti2021). Compared with fixed wings/foils used in traditional vehicles, flapping wings are operated not only at a much higher Strouhal number that is larger than the lower limit for obtaining the dynamic stall (Baik Reference Baik2011), but also in a low Reynolds number range where the flow transition and reattachment can occur more often (Pines & Bohorquez Reference Pines and Bohorquez2006). Despite the importance, active flow control strategies in flapping wings have not received much attention. There have been a few studies on flapping-wing flow controls, including both passive and active control strategies, on unsteady wings or aerofoils. For passive control, great effort has been made towards the effects of the serrations (see e.g. Rao & Liu Reference Rao and Liu2020; Ji et al. Reference Ji, Wang, Ravi, Tian, Young and Lai2022), corrugations (see e.g. Rees Reference Rees1975; Vargas, Mittal & Dong Reference Vargas, Mittal and Dong2008; Bomphrey et al. Reference Bomphrey, Nakata, Henningsson and Lin2016), alula (see e.g. Walker, Thomas & Taylor Reference Walker, Thomas and Taylor2012; Tian et al. Reference Tian, Tobing, Young, Lai, Walker, Taylor and Thomas2019) and covert feathers (see e.g. Carruthers, Thomas & Taylor Reference Carruthers, Thomas and Taylor2007; Wang & Tian Reference Wang and Tian2022) on the force production of wings. These strategies affect the force production by mainly controlling the vortex generation or the boundary layer transition. Few studies have been devoted to active control. For example, Taira & Colonius (Reference Taira and Colonius2009a) investigated the post-stall flow control for an impulsively started translating wing and found that the lift enhancement is due to the roll-up of the trailing-edge sheet into the TV, which presses the LEV onto the wing surface. Synthetic jets are found to be able to delay or promote the shedding of the LEV and suppress the formation of the TEV to improve the aerodynamic performance of a two-dimensional heaving foil at 
$St=0.3$ and 
$Re=100$ (Wang & Tang Reference Wang and Tang2018). The lift enhancement can be achieved by increasing the negative pressure on the suction side of an aerofoil and the circulation of the LEV through suction control, as shown in a two-dimensional numerical study (Feng, Li & Chen Reference Feng, Li and Chen2020). The peak circulation of the LEV is increased by mitigating the formation of the secondary vortex under the LEV through plasma control for a wing at a Strouhal number (St) of 
$St= 0.1$ in the experiments of Kissing et al. (Reference Kissing, Stumpf, Kriegseis, Hussong and Tropea2021). As we can see, the synthetic jet, through blow and/or suction control, could improve the lift force significantly through a variety of mechanisms. This work employs a steady-blowing jet as the control strategy due to its simplicity and proven effectiveness in a stationary foil at very low Reynolds numbers (Taira & Colonius Reference Taira and Colonius2009a). In addition, using steady-blowing jets reduces the parameter space significantly.
Despite the studies introduced above, the active flow control of a three-dimensional low-aspect-ratio flapping wing to enhance the lift generation has not been well explored. A systematic study on using flow controls to enhance the lift generation is desired. This work numerically studies the flow control of a low-aspect-ratio flat-plate flapping wing at an average angle of attack of 
$10^{\circ }$ with a steady-blowing jet. In addition, the effects of Reynolds number and aspect ratio on the effectiveness of the control are considered. The focus is on the control-enhanced lift production and its mechanism through analysis of the vortices and the flow fields around the wing.
The rest of this paper is organised as follows. Section 2 introduces the physical model, the numerical method and the simulation set-up. In § 3, a physical consideration for the flapping-wing flow control is discussed, which is the basis of applying the control strategy. The results and discussions are presented in § 4. Finally, concluding remarks of this work are provided in § 5.
2. Physical model and numerical methods
2.1. Physical model and governing equation
Here, a three-dimensional low-aspect-ratio thin flat-plate wing with a sharp leading edge undergoing a heaving motion with a steady-blowing jet actuator is considered, as illustrated in figure 2, where 
$c$ and 
$b$ are respectively the chord and spanwise length of the wing, 
$U_\infty$ is the incoming flow velocity, 
$\alpha$ is the angle of attack relative to 
$U_\infty$ and 
$\theta$ is the angle between the blowing direction and the positive chord-wise direction of the wing. It is noted that the leading edge is sharp, so that the separation point is fixed at the leading edge, as used in a few previous studies (Shyy et al. Reference Shyy, Trizila, Kang and Aono2009; Thielicke & Stamhuis Reference Thielicke and Stamhuis2015; Eldredge & Jones Reference Eldredge and Jones2019). A thin wing is used because the wing thickness and the leading-edge geometry do not affect the aerodynamics and vortex structures significantly for flapping wings with significant flow separation (also referred to as vortex dominated flows) (Rival et al. Reference Rival, Kriegseis, Schaub, Widmann and Tropea2014; Thielicke & Stamhuis Reference Thielicke and Stamhuis2015). As one of the fundamental kinematics of a complex flapping motion, studying the flow control on a heaving wing helps us to build the required knowledge for designing flow control strategies to enhance the performance of more complicated flapping foils/wings. In order to facilitate the discussions later, two reference frames are introduced: the global inertial frame 
$XYZ$ and the body-fixed frame 
$xyz$, as shown in figure 2. The spanwise length of the wing (
$b$) is fixed to 
$2c$ first to study the flow control effects under an 
$AR=b/c$ of 2, which is around the typical values of small birds adopting flapping flight like sparrows, magpies and pigeons (Tobalske & Dial Reference Tobalske and Dial1996; Azuma Reference Azuma2012). The effects of 
$AR$ will be investigated after the discussions of the cases with 
$AR = 2$. Here, 
$\alpha$ is set to be 
$10^{\circ }$ to promote the formation of the vortical structures during the heaving-down motion. The following equation describes the heaving motion:
\begin{equation} h(t) = h_0^* \cos(2 k t), \end{equation}
where 
$h_0^*=h_0/c=0.5$ and 
$k= {\rm \pi}f c / U_\infty$ are, respectively, the dimensionless heaving amplitude and the reduced frequency with 
$h_0$ being the heaving amplitude and 
$f$ being the heaving frequency. Although the flapping amplitude varies in species, the peak-to-peak flapping amplitude is usually of the same order as the wing chord (Drucker & Lauder Reference Drucker and Lauder2002; Azuma Reference Azuma2012). Therefore, 
$h_0^* = 0.5$ has been widely used in previous studies on flapping wings (Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998; Von Ellenrieder, Parker & Soria Reference Von Ellenrieder, Parker and Soria2003; Dong et al. Reference Dong, Mittal and Najjar2006). The Strouhal number, 
$St=2h_0f/U_\infty$, is fixed at 0.2, within the range of efficient flying or swimming of animals in cruise conditions (Taylor et al. Reference Taylor, Nudds and Thomas2003). Once 
$h_0^*$ and 
$St$ are determined, the reduced frequency can be evaluated to be 
$0.2{\rm \pi}$, which is within the range of the operating reduced frequency of common MAVs (Jones et al. Reference Jones, Duggan and Platzer2001; Jones & Babinsky Reference Jones and Babinsky2010; Decroon et al. Reference Decroon, Percin, Remes, Ruijsink and de Wagter2016).
Figure 2. The schematic of the heaving flat-plate wing with a steady-blowing jet.
The incompressible Navier–Stokes equations govern the flow dynamics over the flapping wing, as the typical operating Mach number of flapping wings is much less than 0.3 for animals and MAVs
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 , \end{gather}$$
$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{ u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{ u} ={-}\boldsymbol{\nabla} p + \frac{1}{Re} \nabla^2 \boldsymbol{u}+ \boldsymbol{ f}, \end{gather}$$
where 
$p$, 
$\boldsymbol {u}$ and 
$\boldsymbol {f}$ are the non-dimensionalised pressure, velocity and body force density, respectively, and 
$Re$ is the Reynolds number defined as 
$Re = U_\infty c / \nu$, with 
$\nu$ being the kinematic viscosity of the fluid.
The steady-blowing jet is applied uniformly along the spanwise direction of the wing, covering the whole span of the wing. It is implemented as a uniform body force along a ghost structure line element (also known as the actuator line) positioned four cells away from the wing, similar to the treatment used by Taira & Colonius (Reference Taira and Colonius2009a). This treatment is to avoid the interference between the boundary of the flapping wing and the blowing jet. It should be noted that specifying the amplitude of the body force corresponding to the jet is much more convenient than specifying the velocity in the feedback immersed boundary method (IBM). To make sure it is comparable to the velocity-based method, we calibrate the amplitude of the body force by conducting the control force in the quiescent flow to achieve the prescribed jet velocity (or strength) (Taira & Colonius Reference Taira and Colonius2009a), which is measured by the momentum coefficient according to
\begin{equation} C_{\mu} = \frac{\rho {U_{act}}^2b\sigma}{1/2\rho U_{\infty}^2bc} = 2 \left(\frac{U_{act}}{U_{\infty}}\right)^2\left(\frac{\sigma}{c}\right), \end{equation}
where 
$\sigma$ is the width of the jet slot, the same as the mesh size near the wing, and 
$U_{act}$ is the maximum flow velocity recorded at the steady state.
2.2. Numerical method
The flow dynamics of the flapping wing described above is solved by the recursive regularised lattice Boltzmann method (LBM) (Chen & Doolen Reference Chen and Doolen1998; Aidun & Clausen Reference Aidun and Clausen2010; Coreixas et al. Reference Coreixas, Wissocq, Puigt, Boussuge and Sagaut2017; Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017) with a feedback IBM (Peskin Reference Peskin2002; Mittal & Iaccarino Reference Mittal and Iaccarino2005; Huang & Tian Reference Huang and Tian2019; Griffith & Patankar Reference Griffith and Patankar2020). In the LBM, the discrete lattice Boltzmann equation describes the evolution of the particle distribution function
\begin{equation} f_i( \boldsymbol{x}+\boldsymbol{e}_i \Delta t,t+\Delta t) - f_i( \boldsymbol{x},t) = \varOmega_i( \boldsymbol{x},t) + \Delta t G_i, \end{equation}
where 
$f_i( \boldsymbol {x},t)$ is the particle distribution function, 
$\boldsymbol {e}_{i}$ is the discrete lattice velocity, 
$\varOmega _i( \boldsymbol {x},t)$ is the collision function, 
$\Delta t$ is the lattice time step and 
$G_i$ denotes the discrete lattice effects due to the body force density 
$\boldsymbol {f}$. Here, the model proposed by Guo, Zheng & Shi (Reference Guo, Zheng and Shi2002a) is used to compute 
$G_i$
\begin{equation} G_i=\left(1-\frac{1}{2\tau}\right)w_i \left(\frac{\boldsymbol{e}_i-\boldsymbol{u}}{c_s^2} + \frac{\boldsymbol{e}_i \boldsymbol{\cdot} \boldsymbol{u}}{c_s^4}\boldsymbol{e}_i\right) \boldsymbol{\cdot} \boldsymbol{f}, \end{equation}
where 
$w_i$ is the weight of the velocity sets of the lattice model, and 
$c_s$ is the ideal speed of sound of the LBM model in lattice units, defined as 
$c_s=\Delta x / (\Delta t \sqrt {3})$, with 
$\Delta x$ being the lattice spacing. Further details of 
$c_s$ and the actual speed of sound in numerical simulations can be found in Krüger et al. (Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017) and Sengupta et al. (Reference Sengupta, Jha, Sengupta, Joshi and Sundaram2023). For a good balance between the accuracy and the computational cost, the three-dimensional nineteen discrete velocity (D3Q19) lattice model is adopted in this work. To enhance the numerical stability at moderate and high Reynolds numbers, the collision model used here is the recursive regularised Bhatnagar–Gross–Krook collision model (Coreixas et al. Reference Coreixas, Wissocq, Puigt, Boussuge and Sagaut2017)
\begin{equation} \varOmega_i(\boldsymbol{x},t) = \frac{1}{\tau}f_i^{(1)}, \quad \tau = \frac{\nu}{c^2_{s} \Delta t} + \frac{1}{2}, \end{equation}
where 
$\tau$ is relaxation time, and 
$f^{(1)}_i$ is the non-equilibrium part of the particle distribution function. The non-equilibrium part 
$f^{(1)}_i$ is evaluated through the recursive regularised steps
\begin{equation} f_i^{(1)} = \frac{w_i}{2c_s^4}\boldsymbol{H}_i^{(2)} : \boldsymbol{\varPi}, \quad \boldsymbol{\varPi} =\sum_{i}\boldsymbol{e}_{i} \boldsymbol{e}_{i} (\,f_i - f_i^{(eq)}), \end{equation}
where 
$\boldsymbol {\varPi }$ is the deviatoric stress tensor, 
$f_i^{(eq)}$ is the equilibrium distribution function and 
$\boldsymbol {H}_i^{(2)}$ is the second-order Hermite tensor (Latt & Chopard Reference Latt and Chopard2006). Including certain third-order Hermite terms into 
$f_i^{(eq)}$ and 
$f^{(1)}_i$ by the recursive procedure can further improve the numerical stability and accuracy. To avoid repetition, readers are referred to Coreixas et al. (Reference Coreixas, Wissocq, Puigt, Boussuge and Sagaut2017) for evaluating the 
$f_i^{(eq)}$ and 
$f^{(1)}_i$ through recursive procedures. The fluid density, pressure and velocity can be evaluated from 
$f_i$ according to the following equations:
\begin{equation} \rho=\sum_{i}f_i,\quad p=\rho c_s^2, \quad \mathrm{and} \quad \rho\boldsymbol{u}=\sum_{i}\boldsymbol{e}_i f_i+\frac{1}{2}\boldsymbol{f}\Delta t.\end{equation}The boundary condition at the fluid–wing interface and the body force used to generate the steady-blowing jet are implemented by the feedback IBM (Goldstein, Handler & Sirovich Reference Goldstein, Handler and Sirovich1993). Owing to its convenience in mesh treatment, the IBM has been extensively applied in flows involving complex geometries and fluid–structure interactions (Mittal & Iaccarino Reference Mittal and Iaccarino2005; Huang & Tian Reference Huang and Tian2019; Griffith & Patankar Reference Griffith and Patankar2020). In the IBM, the body force is used to achieve the boundary conditions
\begin{equation} \boldsymbol{f} ={-}\int \boldsymbol{F}(s,t)\delta( \boldsymbol{x}-\boldsymbol{X}(s,t))\, {\rm d} A,\end{equation}
where 
$\boldsymbol { F}(s,t)$ is the Lagrangian force on the structure exerted by the fluid, 
$dA$ is the surface area of structure elements and 
$\delta$ is the smoothed Dirac delta function. Here, 
$\boldsymbol {F}(s,t)$ at the fluid–wing interface is determined by the feedback scheme
$$\begin{gather} \boldsymbol{F} = \beta(\boldsymbol{U}_{ib} (s,t) - \boldsymbol{U}(s,t)), \end{gather}$$
$$\begin{gather}\boldsymbol{U}_{ib}(s,t) = \int \boldsymbol{u}(\boldsymbol{x},t) \delta( \boldsymbol{x}-\boldsymbol{X}(s,t))\, {\rm d}\kern0.07em\boldsymbol{x}, \end{gather}$$
where 
$\beta$ is the feedback coefficient, 
$\boldsymbol {U}_{ib}(s,t)$ is the fluid velocity on the structure element interpolated from the surrounding fluid and 
$\boldsymbol {U}(s,t)$ is the moving velocity of the structure element. In this work, the feedback coefficient 
$\beta$ is set to be 2.6 in lattice units, as suggested by Huang et al. (Reference Huang, Liu, Wang, Ravi, Young, Lai and Tian2022). For the steady-blowing jet, the magnitude of the body force 
$|\boldsymbol {f}_{act}|$ is defined as
\begin{equation} |\,\boldsymbol{f}_{act}| = \tfrac{1}{1000} \tfrac{1}{2} \rho U_\infty^2, \end{equation}
leading to 
$C_{\mu }=2.03\,\%$, which is large enough to affect the evolution of the vortical structures because the purpose of this work is to explore the control strategies and explain the mechanisms. In the simulation, the Dirac delta function 
$\delta (\boldsymbol {x})$ is replaced by the smoothed Dirac delta function 
$\delta _h (\boldsymbol {x})$ proposed by Peskin (Reference Peskin2002) (see (6.27) in the reference).
At 
$Re = 5000$, the turbulent flow is modelled with large eddy simulation (LES). In the LBM framework, turbulence modelling can be implemented by adding the eddy viscosity 
$\nu _{sgs}$ to the kinematic viscosity to get the effective relaxation time (Malaspinas & Sagaut Reference Malaspinas and Sagaut2012)
\begin{equation} \tau_{eff} = \frac{\nu + \nu_{sgs}}{c^2_s \Delta t} + \frac{1}{2}. \end{equation}
In this work, the dynamic Smagorinsky model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991) is adopted to evaluate 
$\nu _{sgs}$ for its higher accuracy in dealing with a few types of three-dimensional turbulent flows lacking a homogeneous direction according to 
$\nu _{sgs} = (C \Delta x)^2 \sqrt { 2 \bar {\boldsymbol {S}} : \bar {\boldsymbol {S}}}$, where 
$C$ is the Smagorinsky constant, and 
$\bar {\boldsymbol {S}}$ is the instantaneous strain rate tensor evaluated by the second-order central difference scheme. To determine 
$C$ dynamically at every time step, the test-filtered quantities are evaluated by a Gaussian box filter mapped to the D3Q19 lattice model (Touil, Ricot & Lévêque Reference Touil, Ricot and Lévêque2014) due to its simplicity and efficiency in the LBM framework. For example, the test-filtered strain rate tensor can be evaluated according to 
$\widetilde { \overline {\boldsymbol {S}(\boldsymbol {x},t)} } = \sum _{i=0}^{19} w_i \overline {\boldsymbol {S}(\boldsymbol {x} + \boldsymbol {e}_i \Delta t,t)}$. Once all the required grid-filtered and test-filtered quantities are evaluated, 
$C$ is evaluated according to the modified dynamic Smagorinsky model proposed by Lilly (Reference Lilly1992).
2.3. Simulation set-up
Simulations are conducted in a rectangular domain spanning from 
$(-15c, -15c, -15c)$ to 
$(45c, 15c, 15c)$, as shown in figure 3. The multi-block technique (Yu, Mei & Shyy Reference Yu, Mei and Shyy2002; Xu et al. Reference Xu, Tian, Young and Lai2018) is adopted to improve the computational efficiency. The mesh with a high resolution is applied near the wing, and the coarse mesh with a coarsening ratio of 2 is applied away from the wing. Six blocks are used, where the finest mesh of 
$\Delta x = c/160$ and the coarsest mesh of 
$\Delta x = c/5$ are respectively used for the innermost and the outermost blocks. The Dirichlet-type velocity boundary condition and the Neuman-type pressure boundary condition are used for the external boundaries except for the outlet where the Dirichlet-type pressure boundary condition and the Neuman-type velocity boundary condition are applied. The non-equilibrium extrapolation method (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002b) is employed to determine 
$f_i$ at boundary nodes. The absorbing layers (Xu & Sagaut Reference Xu and Sagaut2013) are applied at the lattices near the external boundaries to absorb the outgoing acoustics and minimise the acoustic wave reflections, where the thickness of the absorbing layers is 
$10c$ near the outlet and 
$5c$ near other external boundaries. All cases at 
$Re=1000$ are simulated for four heaving cycles, and the results in the last cycle are used for analysis. At 
$Re=5000$, the results in the last eight cycles are used to evaluate phase average aerodynamic coefficients and flow structures.
Figure 3. A multiblock computation domain except for the outermost block showing the top-port side of the wake around a heaving flat-plate wing of 
$AR = 2$ at 
$Re=5000$. Each block contains a Cartesian mesh grid.
The lift and drag coefficients used to describe the aerodynamic performance of the flapping wing are defined as
\begin{equation} C_L =\frac{L}{\dfrac{1}{2} \rho U_\infty^2 A} , \quad C_D =\frac{D}{\dfrac{1}{2} \rho U_\infty^2 A}, \end{equation}
where 
$L$ and 
$D$ are the lift and drag on the wing, respectively, and 
$A = cb$ is the area of the wing. The force on the wing is evaluated from the Lagrangian force in IBM, as shown in (2.10). The pressure distribution on the wing is calculated at a point 
$3.5 \Delta x$ away from the wing along the surface normal direction by using the linear inverse-distance-weighted interpolation, similar to the treatment in Huang et al. (Reference Huang, Bhat, Yeo, Young, Lai, Tian and Ravi2023). Validations, mesh-independent studies and the cycle-to-cycle variations of the mean drag and lift coefficients have been conducted to make sure that the results and their analysis presented here are reliable, as given in Appendices A and B.
3. Physical considerations for flapping-wing flow control strategy
To consider the flapping-wing flow control strategy in simulations, the force-production mechanism on flapping wings is discussed first. According to conservation of momentum within a control volume 
$V_f$ surrounding the wing, the aerodynamic force on the flapping wing can be written as (Wu, Ma & Zhou Reference Wu, Ma and Zhou2007; Wang et al. Reference Wang, Zhang, He and Liu2015)
\begin{equation} \boldsymbol{F}_W = \rho \int_{V_f} \boldsymbol{u} \times \boldsymbol{\omega} \, {\rm d} {V} - \rho \int_{V_f} \frac{\partial \boldsymbol{u}}{\partial t} \, {\rm d} {V} - \oint_\varSigma p \boldsymbol{n} \, {\rm d} S - \oint_\varSigma \frac{\rho q^2}{2} \boldsymbol{n} \, {\rm d} S + \oint_\varSigma \boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\tau} \, {\rm d} S, \end{equation}
where 
$\varSigma$ is the control surface bounding 
$V_f$, the first term is the vortex force (Saffman Reference Saffman1993) denoted by 
$\boldsymbol {F}_\omega$, the second term is the force due to the fluid acceleration in 
$V_f$ and the remaining terms are forces due to the pressure, kinetic energy and viscosity on the external surface, which are negligible if 
$V_f$ is properly selected (Wang et al. Reference Wang, Zhang, He and Liu2015). The inertia force due to the acceleration of the wing-occupied air is ignored as the thickness of the wing is considered as zero here. As the control strength 
$C_\mu$ is only 2.03 %, the fluid acceleration in 
$V_f$ is not expected to be affected by the control significantly (Taira & Colonius Reference Taira and Colonius2009b). According to previous studies (Wang et al. Reference Wang, Zhang, He and Liu2015; Moriche, Flores & García-Villalba Reference Moriche, Flores and García-Villalba2017), the vortex force contributes to the total force of a flapping wing significantly for a moving body in unsteady viscous flows (Wang et al. Reference Wang, Zhang, He and Liu2015; Moriche et al. Reference Moriche, Flores and García-Villalba2017). Therefore, we only focus on the vortex force
\begin{equation} \boldsymbol{F}_\omega = \rho \int_{V_f} \boldsymbol{u} \times \boldsymbol{\omega} \, {\rm d} {V}, \end{equation}
where 
$\boldsymbol {u} \times \boldsymbol {\omega }$ is the Lamb vector.
The control strategies can be designed based on (3.2). Specifically, the lift due to the vorticity is generated by the vertical (i.e. 
$Y-$) component of the Lamb vector, 
$-u\omega _Z + w\omega _X$. Therefore, increasing the vorticity in the spanwise and streamwise directions is an effective way to increase the lift for a three-dimensional flapping wing. Firstly, viscous flows separate from the flapping wing due to the high angle of attack, where most of the vorticity is concentrated in the LEV instead of being distributed in the boundary layer (Pitt Ford & Babinsky Reference Pitt Ford and Babinsky2013). Therefore, enhancing the circulation of the LEV caused by flapping motions should enhance lift production. To qualitatively explain how to enhance the circulation of the LEV, the two-dimensional assumption, which can be understood for low Reynolds numbers or flows with coherent LEVs/TEVs, is used. The strength of an LEV is determined by the mass flux of the shear layer and the secondary flow induced by the interaction between the LEV and the wing surface, of which the former is considered to be the main contributor of the LEV formation (Wojcik & Buchholz Reference Wojcik and Buchholz2014; Eldredge & Jones Reference Eldredge and Jones2019). Therefore, increasing the mass influx of the shear layer at the leading edge might enhance the LEV circulation for flapping wings, leading to an enhancement of the lift generation. Secondly, TVs of low-aspect-ratio flapping wings have been found to be able to generate lift and delay the shedding of the LEV for flapping wings (Shyy et al. Reference Shyy, Trizila, Kang and Aono2009), which benefits the aerodynamic performance. This is different from the fixed wings, of which the TVs are considered detrimental to lift generation (Anderson Reference Anderson2011) because they induce the downwash near the wing to reduce the effective angle of attack. The TVs provide a positive contribution to the lift of the wing (Shyy et al. Reference Shyy, Trizila, Kang and Aono2009). Unlike the flow past the stationary wing, the circulation distribution near the wing tips is enhanced, leading to an increase in the lift. Therefore, enhancing the TVs of low-aspect-ratio flapping wings could trigger two force generation mechanisms: a direct vortex force and stabilisation of the LEV. Similar to the LEV enhancement strategy mentioned above, enhancing the mass flux of the shear layer at the wing tips might improve the TV circulation.
Unlike the flow separation in the dynamic stall and the post-stall, the boundary layer separates from the leading edge and rolls into unorganised vortices. The forcing frequency must be selected carefully to promote the boundary layer transition or reorganise the small vortices into a large distinct vortex (Wu et al. Reference Wu, Lu, Denny, Fan and Wu1998). The frequency of the shear layers and the vortex shedding is locked to the flapping frequency for insect-bird-like flapping wings due to the high flapping frequency and large flapping amplitude, where the shear layer immediately rolls into a single large vortex. Therefore, steady-blowing jets are used in this work to reduce the parameter space.
Three types of steady-blowing jets are designed. The leading-edge blowing jet is designed to enhance the shear layer at the leading edge to improve the LEV strength. The mid-chord blowing jet is designed to interact with the reattachment point of the LEV and the secondary vortex between the LEV and the wing to delay the detachment of the LEV (Kissing et al. Reference Kissing, Stumpf, Kriegseis, Hussong and Tropea2021). The trailing-edge blowing jet is inspired by Taira & Colonius (Reference Taira and Colonius2009a), which increases the strength of TVs in post-stall flow control for an impulsively started translating low-aspect-ratio wing. The blowing directions are restricted to the directions parallel to the 
$xy$ plane to reduce the complexity and the parameter space.
The control configurations are abbreviated as 
$LE/MD/TE-\theta$ for clarity, where 
$LE$, 
$MD$ and 
$TE$ refer to leading-edge, mid-chord and trailing-edge blowing jets respectively, and 
$\theta$ is the angle between the blowing direction and the positive chord-wise direction (
$x+$) in the body-fixed frame. Simulations with three blowing locations and five blowing angles, 
$0^\circ$ to 
$180^\circ$ in steps of 
$45^\circ$, are conducted, but only the results with relatively large force variation are presented and discussed.
4. Results and discussions
4.1. Trailing-edge downstream blowing jet
4.1.1. Force generation
The trailing-edge downstream blowing (
$TE-0^{\circ }$) jet is discussed first due to its maximum lift enhancement among all control configurations considered in this work. The time histories of the lift and drag coefficients of the baseline cases and the corresponding variations in 
$C_D$ and 
$C_L$ are shown in figure 4, where a few interesting observations are obtained. Firstly, the maximum lift occurs at around 
$t/T=0.2$ within the heaving-down period, which is increased by approximately 
$15.4\,\%$ at 
$Re=1000$ and 
$7.6\,\%$ at 
$Re=5000$ compared with the baseline cases. The relative enhancement in the lift at 
$Re = 5000$ is smaller than that at 
$Re = 1000$. Secondly, during the heaving-up period, the 
$TE-0^{\circ }$ jet also increases the lift coefficient, which is smaller than that in the heaving-down period, especially at 
$Re=5000$. The 
$TE-0^{\circ }$ jet accelerates the flow near the trailing edge, leading to the enhancement in the TV and the pressure drop near the trailing edge, of which the latter can be qualitatively explained by the Bernoulli equation. Thirdly, at both 
$Re = 1000$ and 5000, the variations of 
$C_L$ and 
$C_D$ with the actuation of the 
$TE-0^{\circ }$ jet are similar to those of the corresponding baseline cases. The lift and drag on the wing increase to the maximum at approximately 
$t/T = 0.2$ and decrease to the minimum at approximately 
$t/T=0.7$, which is not affected by the Reynolds number and the 
$TE-0^{\circ }$ jet, indicating that the time when the circulations of the primary vortices reach the maximum does not change under the blowing control. Because of the smaller viscous friction at 
$Re = 5000$, the lift coefficient is slightly higher and the drag coefficient is slightly lower than those at 
$Re=1000$ for both baseline and controlled cases. Because of the positive geometry angle of attack, 
$10^{\circ }$, the viscous friction always contributes to a negative lift. At 
$Re = 5000$, this negative contribution decreases, leading to a slightly larger lift than those at 
$Re = 1000$. The increase in the drag coefficient can be attributed to two components: the increase in the pressure-induced force due to the stronger suction of the lift-production vortices and the increase in the skin friction because of the larger 
$\partial u / \partial y$ near the wing surface caused by the 
$TE-0^{\circ }$ jet. The increase in the drag coefficient due to the blowing effects is also smaller at 
$Re = 5000$ due to the smaller viscosity than that at 
$Re = 1000$. For the two baseline cases at both Reynolds numbers, the difference between the drag coefficients is small during the heaving-down period (
$0\leq t/T \leq 0.5$) because the pressure-induced force from the lift-production vortices dominates the total drag. Finally, during the heaving-up period (
$0.5\leq t/T \leq 1.0$), the flow separation is not as strong as that in the heaving-down period, because the effective angle of attack is smaller. Therefore, the friction drag component is larger than that during the heaving-down period.
Figure 4. The time histories of the aerodynamic coefficients at 
$Re = 1000$ and 5000: (a) lift coefficient of the baseline cases, (b) variations in drag coefficient due to the 
$TE-0^{\circ }$ cases, (c) lift coefficients of the baseline cases and (d) are variations in lift coefficients of the 
$TE-0^{\circ }$ cases.
The average aerodynamic performance coefficients are shown in table 1. From table 1, it is found that the mean lift coefficient is increased by 50.0 % at 
$Re = 1000$ and 22.9 % at 
$Re=5000$. The price for the lift enhancement is the increase in the drag coefficient, of which the mean value is increased by 41.7 % at 
$Re = 1000$ and 23.8 % at 
$Re=5000$. Consequently, there is a slight increase in the lift-to-drag ratio at 
$Re = 1000$ (i.e. 6.0 %) and a slight decrease of 3.5 % at 
$Re=5000$. The observation of the increase in the lift-to-drag ratio at 
$Re = 1000$ is consistent with the results reported by Taira & Colonius (Reference Taira and Colonius2009a), where a more significant increase in the lift-to-drag ratio was observed at a much lower Reynolds number (i.e. 300) rather than the high Reynolds numbers. Therefore, this control strategy is particularly beneficial for the lift enhancement, and thus can be used to improve the maximum lift of the wing to enhance its manoeuvrability. The control strategy further benefits the lift-to-drag ratio at lower Reynolds numbers, and thus can be used to improve the flight efficiency during cruise flight.
Table 1. The average lift, drag and lift-to-drag ratio coefficients for the baseline cases and 
$TE-0^\circ$ cases at 
$Re = 1000$ and 5000.
4.1.2. Flow structures and pressure distributions
To uncover the mechanism of the enhanced lift generation of the controlled cases, the flow structures and the pressure distributions on the wing are shown in figure 5, where the vortical structures are identified by the iso-surface of the Q-criterion (Jeong & Hussain Reference Jeong and Hussain1995; Tian et al. Reference Tian, Dai, Luo, Doyle and Rousseau2014; Sengupta et al. Reference Sengupta, Sharma, Sengupta and Suman2019), and the contours show the pressure coefficient, 
$C_p = (p-p_\infty ) / (0.5\rho U_\infty ^2)$ to highlight the low-pressure cores of the LEV and the TVs. For clarity, the values selected for the iso-surface of the Q-criterion are 10 at 
$Re = 1000$ and 50 at 
$Re = 5000$. The most significant feature of the controlled cases lies in the strengthened TVs which cause the extension of the low-pressure regions near the wing tips towards the trailing edge, making the major contribution to the lift enhancement. This observation is more obvious at 
$t/T=0.4$ than that at 
$t/T=0.2$, and is detailed as follows.
Figure 5. The three-dimensional flow structures coloured by the pressure coefficient, 
$C_p$, and the pressure distribution on the suction side at 
$t/T = 0.2$ (a–d) and 0.4 (e–h). Panels (a) and (e) are the baseline cases at 
$Re = 1000$. Panels (b) and (f) are the 
$TE-0^{\circ }$ cases at 
$Re = 1000$. Panels (c) and (g) are the baseline cases at 
$Re=5000$. Panels (d) and (h) are the 
$TE-0^{\circ }$ cases at 
$Re = 5000$. In each panel, the left image shows the vortical structure, and the right image shows the pressure coefficient on the wing surface.
At the initial stage of the heaving-down period (
$t/T = 0.2$), the flow separates from the leading edge and wing tips. Then, it rolls up into the LEV and TVs immediately for both the baseline cases and the 
$TE-0^{\circ }$ cases at both 
$Re = 1000$ and 5000. The difference in the pressure distribution on the wing of the baseline cases and the 
$TE-0^{\circ }$ cases is not very significant at this time instant, where the pressure near the wing corners is slightly lower for the 
$TE-0^{\circ }$ cases at both 
$Re = 1000$ and 5000. The 
$TE-0^{\circ }$ effects on the vortical structures are demonstrated by the fact that the TEV sheet and the TVs are slightly extended. The suction effects from the LEV and TVs are concentrated at the two corners of the wing because of the detachment of the middle section of the LEV. At 
$Re = 5000$, the differences in the pressure distribution and the vortical structures between the baseline case and the 
$TE-0^{\circ }$ case are smaller than that at 
$Re = 1000$. In addition, the increase in the size of TVs at 
$Re = 5000$ is smaller than that at 
$Re=1000$, as shown in figures 5(f) and 5(h). These factors result in a smaller increase in the lift coefficient for 
$Re = 5000$, as shown in figure 4(c) and table 1.
After the initial development of the LEV and TVs, at 
$t/T = 0.4$, the LEV transitions to an arch shape and travels downstream and inward on the suction side of the wing for all cases, as shown in those snapshots at 
$t/T=0.4$ (see figure 5e–h). Due to the anchoring effects induced by the TVs, the LEV remains attached to the wing, which is similar to the observations from previous studies (Shyy et al. Reference Shyy, Trizila, Kang and Aono2009; Taira & Colonius Reference Taira and Colonius2009b; Visbal, Yilmaz & Rockwell Reference Visbal, Yilmaz and Rockwell2013). The two low-pressure regions directly connected with the arch-like LEV are observed. The TVs also provide a suction force to the wing, which is consistent with the observation that the TVs can provide lift force for flapping wings (Shyy et al. Reference Shyy, Trizila, Kang and Aono2009). When the Reynolds number increases from 1000 to 5000, the vortex structures remain similar to those at 
$Re = 1000$ with smaller scales developed around the LEV and TVs due to the flow instabilities. At 
$Re = 5000$, a secondary low-pressure region appears in front of the low-pressure region caused by the LEV, which is probably caused by the non-continuous LEV ejections due to the shear layer instabilities. Meanwhile, for the baseline case at 
$Re=5000$, the pressure under the LEV and the TVs are relatively lower than that at 
$Re=1000$, resulting in a larger lift.
Under the actuation of the 
$TE-0^{\circ }$ jet, as shown in figures 5(f) and 5(h), an elongated TEV sheet can be observed at both 
$Re = 1000$ and 5000. The interactions between the TEV and the TVs are delayed at 
$Re = 1000$ (figure 5f), similar to that reported in the post-stall flow control (Taira & Colonius Reference Taira and Colonius2009a). The loop-like vortical structure in the wake region has not formed at the same downstream location as the baseline case (figure 5e). The LEV remains similar to that in the baseline case, and the TVs are enhanced clearly, leading to a larger wing area subjected to the suction force at both 
$Re = 1000$ and 5000. The pressure on the regions connected to the arch-like LEV is also lower than that of the baseline case with the aid of the trailing-edge downstream blowing jet at both Reynolds numbers. Careful comparison between figures 5(e) and 5(f) or figures 5(g) and 5(h) shows that the LEVs of controlled cases are closer to the wing, with TVs as the stabiliser. Therefore, it is clear that trailing-edge downstream blowing jet stabilises the LEV, and enhances the strength of the LEV and the TV at 
$Re = 1000$ and 5000, leading to a larger lift on the wing. At 
$Re = 5000$, the control mechanisms are similar to those at 
$Re = 1000$, but the variations in the flow structures and the pressure distribution on the wing are smaller than those at 
$Re = 1000$. Therefore, the enhancement in the lift is smaller than that at 
$Re=1000$.
To study the effects of the steady-blowing jet on the TVs, the circulation and the core centres of the TV are shown in figure 6, where the 
$\varGamma _2$ and 
$\varGamma _1$ criteria proposed by Graftieaux, Michard & Grosjean (Reference Graftieaux, Michard and Grosjean2001) are adopted to estimate the circulation and the vortex centres at different chord-wise locations (in the body-fixed frame) along the TV. The results at 
$Re=5000$ are averaged over eight flapping cycles to reduce the influence of velocity fluctuations due to the small scales in determining the TV boundary. Multiple slices normal to the wing are set across the TV to calculate the circulation and the vortex centres at different chord-wise locations, as shown in figure 6(a). The 
$TE-0^{\circ }$ enhances circulation along the TV at 
$Re = 1000$ and 5000, as shown in figure 6(b). Under the actuation, at 
$Re=1000$, the circulation variation along the TV tube is similar to the baseline case but with a much larger amplitude. At 
$Re=5000$, the circulation enhancement along the TVs is much smaller than that at 
$Re=1000$, and the circulation peaks are slightly shifted downstream under the actuation. Owing to the smaller improvement in the circulation of the TV at 
$Re=5000$, the lift enhancement is not as high as that for 
$Re=1000$. The variations in the TV core centres in figure 6(c) show that the TV is slightly closer to the wing under the actuation due to the traction effects of the 
$TE-0^{\circ }$ jet, acting like a sink to the upstream flow. Still, at 
$Re = 5000$, the variation in the TV centre location is not as large as that at 
$Re = 1000$. To qualitatively estimate the effects of the variation of the TV core centres, we decompose the velocity field 
$\boldsymbol {u}$ as 
$\boldsymbol {U} + \boldsymbol {u}^{\prime }$, where 
$\boldsymbol {U}$ is curl free (irrotational), while 
$\boldsymbol {u}^{\prime }$ is divergence free (solenoidal) and induced by the vorticity. With the decomposition, (3.2) can be rewritten as 
$\boldsymbol {F}_\omega = \rho \boldsymbol {U} \times \int _{V_f} \boldsymbol {\omega } \, {\rm d} \boldsymbol {V} + \rho \int _{V_f} \boldsymbol {u}^{\prime } \times \boldsymbol {\omega } \, {\rm d} \boldsymbol {V}$, which is the linear–nonlinear splitting of the Lamb-vector integral (Wu, Ma & Zhou Reference Wu, Ma and Zhou2015). In the region near TVs, 
$\boldsymbol {U}$ is approximately aligned with the wing, 
$\boldsymbol {u}^{\prime }$ near the wing is pointing inwards and 
$\boldsymbol {w}$ is tilted outwards. Both terms in 
$\boldsymbol {F}_\omega$ make a positive contribution to the lift if 
$V_f$ is selected to cover any one of the TVs. The intention here is not to accurately determine 
$\boldsymbol {U}$ and 
$\boldsymbol {u}^{\prime }$. Instead, we qualitatively estimate 
$\boldsymbol {u}^{\prime }\propto 1/r^2$ according to the Biot–Savart law, where 
$r$ is the distance between the vortex filament and the wing centre. Therefore, the 
$TE-0^{\circ }$ jet strengthens the direct vortex force by enhancing TVs, and drawing the TVs closer to the wing.
Figure 6. The circulation and the core centres of the TV at 
$t/T = 0.4$: (a) the slices in the chord-wise direction along the TV coloured by 
$\omega _x$, (b) the circulation of 
$\omega _x$ and (c) the vertical distance between the TV core centres and the wing surface.
The time histories of the variation in tip suction force coefficient, 
$\Delta C_{s,tip} = \int \Delta C_P \, {\rm d} S$, at 
$Re = 1000$ and 5000 are shown in figure 7(b), where 
$\Delta C_P$ is the pressure difference on the top surface compared with the corresponding baseline case, and the integral region is the one enclosed by the dashed line, as shown in figure 7(a). Firstly, the suction force at the wing-tip region increases during the heaving-down period and recovers to almost zero in the heaving-up period. This observation further shows that the trailing-edge blowing jet increases the wing-tip suction force by enhancing the TV. Secondly, the maximum value of 
$\Delta C_{s,tip}$ at 
$Re = 1000$ is larger than that at 
$Re = 5000$, which is consistent with the previous observations in 
$C_L$ and 
$\varGamma$ along the TV. Finally, the effective operation time of the trailing-edge blowing jet is approximately in the interval of 
$0\leq t/T \leq 0.5$. The jet could be switched off during the heaving-up period to save energy input.
Figure 7. The wing-tip suction force variations at 
$Re = 1000$ and 5000: (a) the schematic showing the region for evaluating the surface integral of the pressure, and (b) time histories of the variation in the suction coefficient compared with the corresponding baseline cases.
The reasons for the enhancement of the TV are analysed here. As discussed in § 3, the growth rate of a vortex is primarily determined by the mass flux of the shear layer feeding into the vortex, and the induced flow between the vortex and the wing (Widmann & Tropea Reference Widmann and Tropea2015; Eldredge & Jones Reference Eldredge and Jones2019). To show this, figure 8 shows the vertical component of the velocity in the body-fixed frame across the shear layers near the wing tips at 
$x/c = 0.5$, 0.7 and 0.9 (as shown in figure 8a) at four instants: 
$t/T =0$, 0.1, 0.3 and 0.4. At 
$t/T = 0$, the beginning of the heaving-down period, the 
$TE-0^{\circ }$ jet already increases the shear layer flux from the mid-chord position to the trailing edge for both 
$Re = 1000$ and 5000, as shown in figures 8(b) and 8(c). The increase at 
$Re = 5000$ is slightly smaller than that at 
$Re = 1000$. At 
$t/T = 0.1$, as the wing starts to heave down, the effects of the 
$TE-0^{\circ }$ jet on the shear layer flux is quite small for both 
$Re = 1000$ and 5000, as shown in figures 8(d) and 8(e), which indicates that the shear layer strength is determined by the heaving motion when the wing start to accelerate. Later, at 
$t/T = 0.3$, as shown in figures 8(f) and 8(g), the gradient of the velocity profiles are slightly smaller than those at 
$t/T = 0.2$ due to the wing starting to decelerate and the increase in the shear layer flux due to the 
$TE-0^{\circ }$ jet is also increased slightly. At 
$t/T = 0.4$, close to the end of the heaving-down period, the shear layer flux is increased at all three chord-wise locations at both 
$Re = 1000$ and 5000, especially at the position near the trailing edge, as shown in figures 8(h) and 8(i). When the heaving velocity is small, the shear layers flux at the wing tips are affected significantly by the 
$TE-0^{\circ }$ jet. Under the actuation, the maximum velocity and velocity gradient of the shear layer are increased due to the larger spatial pressure gradient than the baseline case, as the pressure gradient is the primary source of the generation of vortices on the wing (see e.g. Rosenhead Reference Rosenhead1963; Wu & Wu Reference Wu and Wu1996; Buchholz & Smits Reference Buchholz and Smits2008). Therefore, the TV enhancement due to the larger shear layer mass flux is because the 
$TE-0^{\circ }$ jet increases the spanwise pressure gradient near the wing tip. The 
$TE-0^{\circ }$ jet increases the mass flux across the rear surface of the control volume enclosing the fluid and vortices on the top of the suction side. In order to satisfy the mass conservation within the control volume, the mass flux from other surfaces, especially the shear layer near the wing tip, is enhanced. At 
$Re=5000$, the difference in the shear layer profiles between the baseline case and the 
$TE-0^{\circ }$ case is much smaller than that at 
$Re=1000$, resulting in a smaller improvement in the circulation of the TV.
Figure 8. Shear layer profiles at three chord-wise locations near the wing tip: (a) schematic of three chord-wise locations near the wing tip, (b) 
$t/T=0.0$ at 
$Re = 1000$, (c) 
$t/T=0.0$ at 
$Re = 5000$, (d) 
$t/T=0.1$ at 
$Re = 1000$, (e) 
$t/T=0.1$ at 
$Re = 5000$, (f) 
$t/T=0.3$ at 
$Re = 1000$, (g) 
$t/T=0.3$ at 
$Re = 5000$, (h) 
$t/T=0.4$ at 
$Re = 1000$ and (i) 
$t/T=0.4$ at 
$Re = 5000$.
The mass flux of the shear layer is a function of the shear layer thickness and its velocity distribution. Widmann & Tropea (Reference Widmann and Tropea2015) approximated the velocity distribution of the shear layer feeding into the LEV using a first-order method and found that the temporal evolutions of the LEV circulation with different wing chord lengths collapse onto each other if normalised by the volume flux of the leading-edge shear layer, 
$\delta _{SL} U_\infty$ instead of using 
$cU_\infty$ as the scale, where 
$\delta _{SL}$ is the shear layer thickness. This study shows that the mass flux of the shear layer determines the growth of the vortex. Since the flat-plate wing used in the simulations is infinitely thin, the separation of the shear layer should be independent of the geometry of the edges, as a sharp edge cannot withstand any adverse pressure gradient. Therefore, it is plausible to consider that the properties of the shear layer emanating from the wing edges should be tightly related to the laminar boundary layer before the separation. The thickness of this boundary layer is roughly scaling with the Falkner–Skan solution, 
$\delta _{SL}=(\nu / x_0 U_{SL})^{1/2}$ (Schlichting, Kestin & Shapiro Reference Schlichting, Kestin and Shapiro1961), where 
$x_0$ is a constant determined by the local fluid condition near the stagnation point, and 
$U_{SL}$ is the velocity at the 99 % thickness of the boundary layer. When 
$U_{SL}$ is increased, the thickness of the boundary layer will decrease as 
$\delta _{SL} \sim 1 / U_{SL}^{1/2}$. As we can see, increasing the velocity of the shear layer tends to decrease its thickness, resulting in a self-limiting mechanism of increasing the mass flux of the shear layer. Therefore, at a higher 
$Re$, the effectiveness of the trailing-edge downstream blowing jet on the lift enhancement is reduced.
It should be noted that the same control strategy can also benefit the stationary foil in the lift generation, by an increase of approximately 20 % at an angle of attack (AoA) of 
$40^{\circ }$, which is much smaller than that of the flapping foil in this work and that reported by Taira & Colonius (Reference Taira and Colonius2009b) for a stationary wing at 
$Re=300$. For stationary foils, the stabilisation effect of the steady-blowing jet on LEV and TEV is more obvious at low Reynolds number, while the TV enhancement is weaker than the flapping wing.
4.2. Leading-edge and mid-chord blowing jets
4.2.1. Force generation
The control effects of leading-edge and mid-chord blowing jets are considered in this section. As the LEV contributes to the force significantly in flapping wings, the leading-edge blowing (
$LE-\theta$) and the mid-chord blowing (
$MD-\theta$) jets are expected to modify the shear layer feeding into the LEV directly and interact with the LEV near the reattachment point to affect its detachment, respectively. Because § 4.1 shows that the results at 
$Re = 1000$ and 5000 are similar for both baseline and controlled cases, only the results at 
$Re = 1000$ are discussed here. Three controlled cases are selected to be discussed in this section, according to the variations in aerodynamic performance. The time histories of the aerodynamic force coefficients of other blowing configurations are included in Appendix B due to the relatively small variations compared with the baseline case.
The time histories of the variations in lift and drag coefficients are shown in figure 9, and the mean lift, drag and lift-to-drag coefficients are given in table 2. Based on our simulations, two major observations are obtained. Firstly, the lift production of these cases is not as fruitful as the 
$TE-0^{\circ }$ jet discussed in § 4.1. Specifically, the 
$MD-90^{\circ }$ jet enhances the lift during the whole flapping cycle, and the trend of 
$C_D$ and 
$C_L$ is similar to that of the baseline case. The lift enhancement is smaller than that of 
$TE-0^{\circ }$ with an increase of 34.3 % in the average lift coefficient and an increase of 9.4 % in the instantaneous lift at 
$t/T=0.2$. Secondly, some cases could enhance the lift-to-drag ratio significantly. For example, the 
$MD-180^{\circ }$ jet decreases the average drag by 36.0 %, resulting in a 50.2 % higher 
$\overline {C_L} / \overline {C_D}$ than the baseline case, and the variation in the 
$C_L$ is small, with an increase of 10.0 % in the average lift and 26.0 % in the instantaneous lift at 
$t/T=0.7$. The 
$LE-180^{\circ }$ jet not only increases the negative drag (thrust) significantly but also extends the thrust-production period at the beginning of the heaving-up period, leading to an 18.0 % lower average drag and a 22.1 % higher 
$\overline {C_L} / \overline {C_D}$ than the baseline case. This might be due to the 
$LE-180^{\circ }$ jet enhancing the formation of the LEV effectively during the heaving-up period.
Figure 9. The time histories of the variations in lift and drag coefficients with respect to the baseline at 
$Re = 1000$: (a) variation in the lift coefficient, and (b) variation in the drag coefficient.
Table 2. The average lift, drag and lift-to-drag ratio coefficients for the baseline case and those under leading-edge, mid-chord and trailing-edge blowing jets at 
$Re = 1000$.
4.2.2. Flow structures and pressure distributions
The three-dimensional flow structures and the pressure distributions for the baseline case and the 
$LE-180^{\circ }$ case are shown in figure 10, while the vorticity contours at the mid-span section of the wing are shown in figure 11. At 
$t/T=0.2$, a larger LEV tube is initiated at the leading edge under the 
$LE-180^{\circ }$, as shown in figures 10(a) and 10(b). However, this larger LEV provides less suction force for the middle section of the wing compared with that of the baseline case, resulting in a small variation in the lift and drag coefficients during the heaving-down period (see figure 9). As shown in figures 10(b) and 11(b), not only is the LEV advected away from the wing surface under the actuation, it is also not as compact as that in the baseline case. At 
$t/T = 0.7$, in contrast to the heaving-down period, the 
$LE-180^{\circ }$ jet enhances the formation of the LEV, rendering larger negative lift and drag on the wing than that of the baseline case, as shown in figures 10(c) and 10(d). The increase in the negative lift is relatively small because the angle of attack is 
$10^{\circ }$, where the drag is primarily from the horizontal component of the vortical suction force. Although the middle section of the LEV starts to detach from the wing surface, the overall suction force and the wing area subjected to the LEV suction are still larger than that in the baseline case due to the anchoring effects caused by the TVs. In figures 11(c) and 11(d), it can be seen that the LEV identified by vorticity contours at the middle spanwise location of the wing is much larger than that of the baseline case at the middle spanwise location of the wing, and the LEV keeps its coherent shape well, in contrast to that at 
$t/T=0.2$. Therefore, the effectiveness of the 
$LE-180^{\circ }$ jet depends on the effective angle of attack at the initial LEV development stage. The effective angle of attack would affect the angle between the shear layer and the wing surface, which has been known to relate to the growth and detachment state of the LEV (Li et al. Reference Li, Feng, Kissing, Tropea and Wang2020). The leading-edge upstream blowing jet can enhance the LEV formation while keeping its coherent shape well if the effective angle of attack is small at the LEV initial development stage.
Figure 10. The three-dimensional flow structures coloured by the pressure coefficient and the pressure distribution on the suction side of the wing at 
$t/T = 0.2$ and 0.7 at 
$Re=1000$: (a) baseline case at 
$t/T = 0.2$, (b) 
$LE-180^{\circ }$ case at 
$t/T = 0.2$, (c) baseline case at 
$t/T = 0.7$ and (d) 
$LE-180^{\circ }$ case at 
$t/T = 0.7$.
Figure 11. The vorticity contours at the mid-span section of the wing at 
$Re=1000$: (a) baseline case at 
$t/T = 0.2$, (b) 
$LE-180^{\circ }$ case at 
$t/T = 0.2$, (c) baseline case at 
$t/T = 0.7$ and (d) 
$LE-180^{\circ }$ case at 
$t/T = 0.7$.
The three-dimensional flow structures and the pressure distributions on the wing with the actuation of the mid-chord blowing jets are shown in figure 12. As shown in figure 12(a), the 
$MD-90^{\circ }$ jet induces a vortical structure covering the entire wing span with a smaller size approaching the wing tip, which is similar to the LEV. This blowing-induced middle vortex causes an additional suction force on the wing, leading to a larger lift than the baseline case. Regarding the 
$MD-180^{\circ }$ jet, shown in figure 12(b), the blowing jet from the downstream of the LEV destroys the coherent shape of the LEV during the initial growth of the LEV. The low-pressure region under the LEV disappears, leading to a smaller drag and lift on the wing than the baseline case.
Figure 12. The three-dimensional flow structures and the pressure distribution on the suction side of the wing at 
$t/T = 0.2$ at 
$Re=1000$ with the actuation of the mid-chord blowing jets: (a) 
$MD-90^{\circ }$ case, (b) 
$MD-180^{\circ }$ case.
The vorticity contours superposed by the streamlines with the actuation of the mid-chord blowing jet at the mid-span section are shown in figure 13. As shown in figure 13(a), the 
$MD-90^{\circ }$ jet changes the topological characteristics of the vortex–wing system, which has been known to affect the formation and detachment of the vortex on flapping wings (Perry & Chong Reference Perry and Chong1987; Rival et al. Reference Rival, Kriegseis, Schaub, Widmann and Tropea2014). With the actuation, two additional nodes (the two green nodes), one full saddle point (the red node) and two half-saddle points (the two blue nodes) are formed. As shown in figure 13(b), the 
$MD-180^{\circ }$ jet drives the reversed flow from the trailing edge at the early heaving-down period, which triggers the ‘bluff body detachment’ mechanism (Widmann & Tropea Reference Widmann and Tropea2015) of the LEV during the initial growth stage of the LEV. When the LEV detaches due to the initialisation of the reversal flow, it is called the ‘bluff body detachment’ type. Another one is the ‘vorticity layer eruption’ type, where the LEV detaches due to the vorticity eruption in the boundary layer (Doligalski, Smith & Walker Reference Doligalski, Smith and Walker1994).
Figure 13. The vorticity contour superposed by the streamlines at the mid-span section of the wing at 
$t/T = 0.2$ at 
$Re=1000$: (a) 
$MD-90^{\circ }$ case, and (b) 
$MD-180^{\circ }$ case. In (a), the green, blue and red points represent the node, half-saddle and full saddle.
4.3. Effects of aspect ratios
4.3.1. Force generation
The effects of the aspect ratio (
$AR$) on flapping-wing flow control are explored in this section by considering 
$AR = 1$, 2 and 4, as 
$AR$ plays a significant role in force generation and vortex dynamics of low-aspect-ratio flapping wings (Dong et al. Reference Dong, Mittal and Najjar2006). Four blowing configurations are selected according to the results of 
$AR = 2$. The first case is the leading-edge upstream blowing jet because it enhances the LEV formation during the heaving-up period effectively to increase the thrust; the second case is the mid-chord vertical blowing jet because it improves the lift by forming a blowing-induced middle vortex; the third case is the mid-chord upstream blowing jet because it promotes the ‘bluff body detachment’ of the LEV to decrease the drag, resulting in an increase in the lift-to-drag ratio; and the last case is the trailing-edge downstream blowing jet because of the highest lift increase in all cases of 
$AR = 2$.
The time histories and the average lift and drag coefficients for 
$AR = 1$ and 4 are given in figure 14 and table 3. With the actuation of the 
$LE-180^{\circ }$ jet, the drag coefficient decreases quickly at the beginning of the heaving-up period for 
$AR = 1$, similar to that for 
$AR = 2$. For 
$AR =4$, however, the drag reduction magnitude and period are smaller and shorter than those for 
$AR = 1$ and 2, indicating that the thrust enhancement of the 
$LE-180^{\circ }$ jet diminishes with increasing 
$AR$. This may imply that the variations in the LEV evolution during the heaving-up period decrease when increasing 
$AR$ from 1 to 4, resulting in a smaller force variation for 
$AR = 4$. There are 24.0 %, 17.9 % and 1.2 % reductions in 
$\overline {C_D}$ for 
$AR = 1$, 2 and 4, respectively. Similarly, with the 
$MD-180^{\circ }$ jet, the control effects diminish with increasing 
$AR$ from 1 to 4 regarding the variations in the lift and drag coefficients. The enhancement on 
$\overline {C_L} / \overline {C_D}$ decreases when increasing the 
$AR$, where they are 77.8 %, 50.1 % and 44.6 % increases in 
$\overline {C_L} / \overline {C_D}$ for 
$AR = 1$, 2 and 4, respectively. Regarding the 
$MD-90^{\circ }$ and 
$TE-0^{\circ }$ control, 
$AR$ does not affect flow control effects much in terms of the variations in the magnitude and the trend of the lift and drag coefficients. For the 
$MD-90^{\circ }$ blowing configuration, the increases in 
$\overline {C_L}$ are 36.2 %, 34.6 % and 28.4 % for 
$AR = 1$, 2 and 4, respectively, where the lift enhancement decreases with increasing 
$AR$ but not much. For the 
$TE-0^{\circ }$ blowing configuration, the increases in 
$\overline {C_L}$ are 50.8 %, 50.0 % and 45.3 % for 
$AR = 1$, 2 and 4, indicating that the lift enhancement is not very sensitive to 
$AR$ for low-aspect-ratio wings.
Figure 14. The time histories of the drag and lift coefficients of baseline cases and their changes due to the control for 
$AR = 1$ and 4 (with markers) at 
$Re=1000$: (a) drag coefficient of baseline case, (b) changes of drag coefficient due to the control, (c) lift coefficient of baseline case and (d) changes of lift coefficient due to the control.
Table 3. The average lift, drag and lift-to-drag ratio coefficients for the baseline and controlled cases of 
$AR = 1$ and 4 at 
$Re = 1000$.
4.3.2. Flow structures and pressure distributions
The flow structures of the baseline and 
$LE-180^{\circ }$ cases for 
$AR = 1$ and 4 are shown in figure 15. From figure 15(a), a much more apparent vortex loop is formed for the wing of 
$AR =1$ compared with that of 
$AR = 4$ (figure 15c) due to the stronger interaction between the TVs and the TEV for a low-
$AR$ wing. With the actuation of the 
$LE-180^{\circ }$ jet, the LEV is enhanced during the heaving-up period for both 
$AR = 1$ and 4, as shown in figures 15(b) and 15(d), where the size of the LEV and the core pressure are larger and lower than those of the corresponding baseline cases. Therefore, thrust production is enhanced, as shown in figures 14(b) and 14(d). For the wing with 
$AR = 4$, only the middle section of the LEV is enhanced, resulting in a relatively smaller thrust enhancement, where the relative increase in the LEV suction area to that of the baseline case is much smaller than that for 
$AR = 1$. Overall, it can be seen that the leading-edge upstream blowing jet is more effective for a heaving wing with a lower 
$AR$.
Figure 15. The three-dimensional flow structures coloured by the pressure coefficient and the pressure distributions on the suction side of the wing at 
$t/T = 0.6$: (a) baseline case for 
$AR = 1$, (b) 
$LE-180^{\circ }$ for 
$AR = 1$, (c) baseline case for 
$AR = 4$ and (d) 
$LE-180^{\circ }$ for 
$AR = 4$.
The vortical structures of the baseline and the 
$MD-180^{\circ }$ cases of 
$AR = 1$ and 4 are shown in figure 16. As 
$AR$ increases from 1 to 4, the middle portion of the LEV shows two-dimensional behaviour due to weaker three-dimensional effects, as shown in figures 16(a) and 16(c). The LEV for 
$AR = 4$ is more robust to the 
$MD-180^{\circ }$ jet than that for 
$AR = 1$, maintaining low pressure at its vortex core and hence the suction force to the wing despite the formation of the small streamwise vortices. For 
$AR = 1$, the LEV cannot maintain its low pressure in the core, and the suction force on the wing is significantly reduced. With the actuation of the 
$MD-180^{\circ }$ jet, the size of the streamwise vortices is comparable to the LEV for 
$AR = 1$, which is roughly 
$0.5c$ in the spanwise direction. Based on the observations for 
$AR = 1$, 2 and 4, the LEV for a large 
$AR$ appears to be more stable than that of a small 
$AR$ in the presence of the streamwise vortices caused by the 
$MD-180^{\circ }$ jet because of the larger ratio between the size of the LEV and the streamwise vortices. Therefore, the enhancement of 
$\overline {C_L} / \overline {C_D}$ decreases as 
$AR$ increases from 1 to 4.
Figure 16. The three-dimensional flow structures coloured by the pressure coefficient and the pressure distributions on the suction side of the wing at 
$t/T = 0.2$: baseline case for 
$AR = 1$, (b) 
$MD-180^{\circ }$ case for 
$AR = 1$, (c) baseline case for 
$AR = 4$ and (d) 
$MD-180^{\circ }$ case for 
$AR = 4$.
5. Conclusions
The flow control of a low-aspect-ratio heaving wing through a steady-blowing jet has been investigated numerically by a feedback immersed boundary–LBM. It is found that the trailing-edge downstream blowing jet increases lift production most among all configurations at both 
$Re = 1000$ and 5000 by increasing the net circulation on the wing. An extended and stronger TV is formed by the increase in the mass flux of the shear layer near the wing tip, resulting in a larger suction area with lower pressure. However, this mass flux increase is reduced at 
$Re=5000$ due to its self-limiting mechanism, resulting in a lower TV circulation enhancement, and thus a lower lift generation. In addition, the pressure on the part of the wing connected to the arch-like LEV is lower, indicating that the LEV is strengthened, despite the non-spanwise orientation of the LEV.
The leading-edge upstream blowing jet enhances the LEV at the beginning of the heaving-up period, resulting in higher thrust production. However, in the heaving-down period, it cannot effectively enhance thrust production because the blowing jet advects the LEV away from the wing, and the LEV is not as compact as that of the baseline case. The mid-chord vertical blowing jet increases the lift by forming a blowing-induced middle vortex due to the modification of the topological characteristics of the vortex–wing system. The mid-chord upstream jet destroys the LEV coherent structure and promotes its detachment by triggering the ‘bluff body detachment’ mechanism early at the growth stage. Although the drag and lift are reduced due to the detachment of the LEV, the 
$\overline {C_L} / \overline {C_D}$ is significantly enhanced.
The aspect ratio (
$AR$) does not significantly influence the flow control effects except for the leading-edge upstream blowing and the mid-chord upstream blowing jets. With the actuation of the leading-edge upstream blowing jet, as 
$AR$ increases from 1 to 4, the thrust enhancement during the initial heaving-up period is reduced because only the middle portion of the LEV is enhanced for 
$AR = 4$. For 
$AR = 4$, the LEV is more robust against the mid-chord upstream blowing jet, probably due to the larger ratio between the size of the LEV and the blowing-induced streamwise vortices.
Acknowledgements
We would like to thank Dr J. Kissing at TU Darmstadt for fruitful discussions.
Funding
This work was supported by the Australian Research Council Discovery Project (grant no. DP200101500) and conducted with the assistance of computational resources from the National Computational Infrastructure (NCI), which is supported by the Australian Government.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Solver validation and mesh independence study
A.1. Solver validation at 
$Re=10\ 000$
The immersed boundary(IB)–LBM has been extensively validated and applied in several biological flows (Xu et al. Reference Xu, Tian, Young and Lai2018; Ma et al. Reference Ma, Wang, Young, Lai, Sui and Tian2020; Huang et al. Reference Huang, Liu, Wang, Ravi, Young, Lai and Tian2022). Here, a heaving SD7003 aerofoil section with spanwise length of 
$0.2c$ in a uniform flow is conducted to validate the capability of the current solver to model three-dimensional flapping wings at a moderate Reynolds number. The aerofoil is heaving according to (2.1) at an angle of attack of 
$4^{\circ }$, 
$h_0=0.05c$, 
$k=3.93$ and 
$Re=10\ 000$. The computational domain is similar to that defined in § 2.3, except the spanwise direction 
$(z)$ for which 
$0.2c$ is used, the same as the spanwise length of the aerofoil. The boundary conditions are also similar to those defined in § 2.3, except for the spanwise boundaries where the periodic boundary condition is applied. Figure 17 shows the drag and lift coefficients from the current solver agree well with the result in the high fidelity simulations of Visbal (Reference Visbal2009). Further validations of LES can be found in Wang et al. (Reference Wang, Liu, Jin, Huang, Young and Tian2024).
Figure 17. The time histories of the drag and lift coefficients for a heaving SD7003 aerofoil in a uniform flow: (a) drag coefficient, and (b) lift coefficient.
A.2. Mesh independence study
To study the mesh convergence, the baseline cases with three different grids around the wing, shown in figure 2, at 
$Re=10\ 000$ are performed. The finest grid around the wing is 
$c/80$, 
$c/160$ and 
$c/320$ for the coarse, medium and fine grids, respectively. The phase average time histories of drag and lift coefficients are shown in figure 18. As shown in table 4, the relative variation in the mean drag and lift coefficients for the medium grid agrees with those for the fine grid. Therefore, the medium grid with medium mesh resolution is considered to be sufficient for this work.
Table 4. Summary of the mesh independence study: 
$\Delta x_{min}$ is the finest grid of the inner block surrounding the wing; 
$N$ is the number of lattices; and 
$\epsilon$ is the fractional change of the aerodynamic coefficients after refining.
Here, the periodicity of the aerodynamic force is examined. As shown in figure 19, the simulations can be considered to achieve the statistically steady state from the third heaving cycle at both 
$Re = 1000$ and 5000.
Figure 18. The time histories of the phase average drag and lift coefficients of the baseline case at 
$Re=10\ 000$: (a) drag coefficient, and (b) lift coefficient.
Figure 19. The relative variations of the cycle average drag and lift coefficient to the previous heaving cycle: (a) drag coefficient at 
$Re = 1000$, (b) lift coefficient at 
$Re = 1000$, (c) drag coefficient at 
$Re = 5000$ and (b) lift coefficient at 
$Re = 5000$.
Figure 20. The time histories of the drag and lift coefficients of the wing of 
$AR = 2$ with the actuation of all blowing configurations at 
$Re = 1000$. Panels (a,c,e) are the drag coefficients for the leading-edge, mid-chord and trailing-edge blowing, respectively. Panels (b,d,f) are the lift coefficients for the leading-edge, mid-chord and trailing-edge blowing, respectively.
Appendix B. Aerodynamic force variations at $Re=1000$
The time histories of the aerodynamic coefficient of the wing of 
$AR = 2$ under the actuation of all blowing configurations at 
$Re = 1000$ are shown in figure 20. For the leading-edge blowing, as shown in figures 20(a) and 20(b), the drag reduction effect at the beginning of the heaving-up period increases with increasing the blowing angle from 
$0^{\circ }$ to 
$180^{\circ }$, suggesting that the LEV strength is improved due to increasing the blowing angle. Compared with 
$LE-135^{\circ }$ and 
$LE-180^{\circ }$ jets, other leading-edge blowing configurations do not affect the aerodynamic coefficients significantly during the whole flapping cycle. Regarding the mid-chord blowing, as shown in figures 20(c) and 20(d), all the blowing configurations do not change the trend of 
$C_D$ of 
$C_L$ much. On increasing the blowing angle, 
$C_D$ is reduced, leading to increased 
$C_L / C_D$. For the trailing-edge blowing, as shown in figures 20(e) and 20(f), the 
$TE-0^{\circ }$ jet could increase 
$C_L$ significantly. On increasing the blowing angle, the lift improvement degrades. Regarding 
$C_D$, it decreases with increasing the blowing angle because the suction force from the low-pressure region of fluids decreases, consistent with the variations in 
$C_L$.
