1. Introduction
Insects fly by moving their wings rapidly back and forth in a complex pattern. With this wing motion, they generate the aerodynamic forces required for flight. Due to the unsteady nature of the involved flow dynamics and the complexity of the wing motion, the interaction between the wings and the surrounding air is difficult to capture accurately with simple models.
A common approach to describe and predict the effects of this interaction is by using the concept of aerodynamic mechanisms (Ellington Reference Ellington1984b), also known as phenomenology-based models (Zhang, Hedrick & Mittal Reference Zhang, Hedrick and Mittal2015). These semi-empirical models relate the aerodynamic forces exerted on a flapping wing directly to the instantaneous motion of the wing and its geometry. Each individual aerodynamic mechanism provides a prediction of the aerodynamic forces that are related to specific aspects of the motion, orientation and geometric characteristics of the wing (e.g. the wing speed, angle of attack and surface area). Several such aerodynamic mechanisms have been identified, including translational lift (stroke-rate-related forces) (Ellington et al. Reference Ellington, van den Berg, Willmott and Thomas1996; Dickinson, Lehmann & Sane Reference Dickinson, Lehmann and Sane1999; Birch & Dickinson Reference Birch and Dickinson2001; Birch, Dickson & Dickinson Reference Birch, Dickson and Dickinson2004), rotational lift and drag (wing-pitch-rate-related forces) (Dickinson et al. Reference Dickinson, Lehmann and Sane1999; Nakata, Liu & Bomphrey Reference Nakata, Liu and Bomphrey2015; Bomphrey et al. Reference Bomphrey, Nakata, Phillips and Walker2017; van Veen, van Leeuwen & Muijres Reference van Veen, van Leeuwen and Muijres2019) and added-mass forces (acceleration-reaction forces) (Sedov Reference Sedov1965; Walker & Westneat Reference Walker and Westneat2000; Sane & Dickinson Reference Sane and Dickinson2001; Whitney & Wood Reference Whitney and Wood2010; Kang et al. Reference Kang, Aono, Cesnik and Shyy2011). These effects are associated with the instantaneous motion of the wing, and are therefore also commonly included in so-called quasi-steady aerodynamic models to describe the aerodynamics of insect flight (Ellington Reference Ellington1984a; Sane & Dickinson Reference Sane and Dickinson2002; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Lee et al. Reference Lee, Lua, Lim and Yeo2016; Wang, Goosen & Van Keulen Reference Wang, Goosen and Van Keulen2016; Cai et al. Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021).
Additional effects, which depend on the time history of the wing motion, are wake capture (interaction between the wing and the previously generated wake) (Dickinson Reference Dickinson1994; Sane & Dickinson Reference Sane and Dickinson2001; Birch & Dickinson Reference Birch and Dickinson2003; Lehman Reference Lehman2008) and the Wagner effect (delayed build-up of bound circulation as a result of inertia of the fluid) (Wagner Reference Wagner1925; Dudley & Ellington Reference Dudley and Ellington1990; Dickinson & Götz Reference Dickinson and Götz1996; Walker & Westneat Reference Walker and Westneat2000; Sane Reference Sane2003). Due to their dependence on the time history of airflow development and wing-beat kinematics, these mechanisms are generally considered incompatible with a quasi-steady modelling approach.
Here, we combined numerical simulations with aerodynamic modelling to study the unsteady aerodynamics of insect flight. In our proposed quasi-steady aerodynamic model, we included an extended model of the added-mass mechanism, and a novel model that captures the Wagner effect on aerodynamic force production (table 1). For this study, we used two-winged Diptera as the model insect.
Overview of several quasi-steady aerodynamic models used to study the aerodynamics of flapping insect wings, and our new model for accelerating wings. For each model, we highlight which aerodynamic mechanism is included. These are the translational-force mechanism, the rotational-lift mechanism, added mass, the Wagner effect and wake capture. For the added-mass mechanism, we also highlight whether the model includes tangential forces, and whether the added fluid volume is independently estimated. For the wake capture, no quasi-steady aerodynamic force model exists; instead, here we highlight whether the wing-beat kinematics used for model development included stroke reversal, and thus for which wake-capture effects were present in the data. 
$^{*}$Added-mass is included only for the power requirement model, not for the force model. 
$^{\S }$ The model was developed using kinematics without wake capture (no stroke reversal), but it was tested on wings with stroke reversal. 
$^{+}$The Wagner effect was applied as a correction only for an accelerating wing, not for a flapping wing. 
$^{{\dagger} }$The current model focuses on the wing-stroke motion alone.
1.1. Wing-beat kinematics and aerodynamics of flying Dipteran insects
Among insects, the order Diptera is a highly diverse group of flyers that exhibit a large range of wing-beat patterns (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014, Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a; Wiegmann & Yeates Reference Wiegmann and Yeates2017). Common among Diptera is that they all have two wings that they beat back and forth at high angles of attack (of the order of 45
$^{\circ }$) (figure 1a,b).
(a,b) Schematic representation of a hovering fruit fly (Drosophila hydei), including the wing motion during the forward wing stroke (a) and the backward wing stroke (b). The lollypop symbols show the wing chord, whereby the circle depicts the leading edge of the wing. (c–e) Wing-beat kinematics of a hovering fruit fly (blue) (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014) and malaria mosquito (Anopheles coluzzii) (red) (Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a), including angle of attack (c), stroke rate (d) and stroke acceleration (e). The dot indicates stroke reversal. ( f) Definition of the world reference frame, wing reference frame, the stroke angle and angle of attack. (g,h) Design of the kinematics used in this study, showing the stroke rate (g) and stroke acceleration (h) throughout the simulation. The colour of each line in panels g and h indicates the stroke acceleration at the end of the simulation (see colour bar). Here, green data show accelerating wings, purple show decelerating wings and grey show that of a non-accelerating wing.
The sweeping wing-stroke motion of these animals can be divided into two phases: the forward stroke and backward stroke (figures 1a and 1b, respectively). During the forward stroke, the wing is moving from the dorsal side of the insect towards the ventral side, and in the opposite direction for the backward stroke. Stroke reversal is defined as the phase where the wing transitions from the forward stroke to the backward stroke, and vice versa. Around stroke reversal, the wing has to first decelerate and then accelerate rapidly in the opposite direction. Furthermore, the wing rotates around its spanwise axis over an angle of around 90
$^{\circ }$, such that the angle of attack is again close to 45
$^{\circ }$ during the next stroke (figure 1c). Due to the high flapping frequency (of the order of 100 Hz), this stroke reversal has to occur in a very short time and is thus accompanied with very high wing-stroke and wing-pitch accelerations (of the order of 
$10^{6}$ rad s
$^{-2}$) (Dickinson & Muijres Reference Dickinson and Muijres2016; Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a), see also figure 1(e).
These wing accelerations result in added-mass force production, while accelerating wings can further experience a Wagner-effect-related phenomenon (Sane Reference Sane2003; Wang et al. Reference Wang, Goosen and Van Keulen2016). Currently, several models have been developed for flapping insect wings that include the added-mass effect (Sane & Dickinson Reference Sane and Dickinson2002; Andersen, Pesavento & Wang Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Lee et al. Reference Lee, Lua, Lim and Yeo2016; Wang et al. Reference Wang, Goosen and Van Keulen2016; Cai et al. Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021), but none of these take the Wagner effect into account (table 1). Without a proper Wagner-effect model the description of the aerodynamics of rapidly oscillating insect wings is still incomplete.
As a result, the aerodynamics of flapping insect wings is often modelled using simplified versions of these acceleration-related models, or these effects are ignored altogether (table 1). The most simplified quasi-steady aerodynamic models for insect flight ignore the acceleration-based mechanisms and are thus based on only on the velocity-dependent mechanisms. These are the mechanism based on stroke rate of the wing (translational lift) and the mechanism based on wing-pitch rate (rotational lift) (Ellington Reference Ellington1984a; Dudley & Ellington Reference Dudley and Ellington1990; Dickinson & Muijres Reference Dickinson and Muijres2016).
Most quasi-steady aerodynamic models that do include acceleration-based aerodynamics use simplified models of the added-mass mechanism. These models are based on potential-flow theory, assume a constant added fluid volume and produce only wing-normal pressure forces as the wing is considered very thin (Sane & Dickinson Reference Sane and Dickinson2002; Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Wang et al. Reference Wang, Goosen and Van Keulen2016). Only a few studies estimated the added volume for flapping insect wings independently (Lee et al. Reference Lee, Lua, Lim and Yeo2016; Cai et al. Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021), or modelled tangential added-mass forces based on potential-flow theory (Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007). None of these studies, however, tested the assumptions used to estimate tangential added-mass force production.
Because of the lack of a usable Wagner-effect model, all these quasi-steady models for insect flight have ignored the Wagner effect on aerodynamic force production, except for Wang et al. (Reference Wang, Goosen and Van Keulen2016) (table 1). In their study, the authors added a Wagner-effect-based correction to the translational (circulatory) forces produced by a wing starting from rest, based on the classic Wagner theory, but this model could not be applied to flapping wings.
1.2. Unsteady aerodynamics of insect flight
In this section, we shortly review the characteristics of the different unsteady aerodynamic mechanisms identified above and how they can potentially be treated with a (quasi-steady) modelling approach. The focus of our study is on the stroke-related effects, and therefore these will receive particular attention. In the subsequent analysis, we will accordingly consider wings that perform only a stroke motion, without wing-pitch movements. In view of this, the particular scope of the present study is to reconsider the stroke-related components of the quasi-steady modelling of flapping wings. The definition of the coordinate system used in the analysis is illustrated in figure 1( f).
The stroke-rate-based forces (translational lift) depend on the angular stroke velocity of the wing and the angle of attack. Because insect wings tend to operate at high angles of attack, the stroke-rate-based forces in insect flight are often associated with the occurrence of a leading-edge vortex (LEV). The LEV constitutes a low-pressure region of separated flow that is located above the wing. The flow separates at the leading edge and reattaches before the trailing edge. The LEV is commonly identified as the mechanism by which insect wings can operate at such high angles of attack without experiencing stall, and consequently produce the high lift force coefficients supporting flight (Maxworthy Reference Maxworthy1979; Ellington et al. Reference Ellington, van den Berg, Willmott and Thomas1996; Dickinson et al. Reference Dickinson, Lehmann and Sane1999; Lentink & Dickinson Reference Lentink and Dickinson2009).
Because the circulation of this LEV and associated forces depends on the translational movement of the air across the wing surface, we opt here to refer to the stroke-rate-based forces as the commonly used translational forces (
$\boldsymbol {F}_{transl}$). Note, however, that the flapping insect wing itself makes a (curvilinear) revolving stroke motion and thus does not translate in the strict rectilinear sense.
Translational forces can be modelled in relation to the wing geometry and kinematics as (Dickinson et al. Reference Dickinson, Lehmann and Sane1999)
\begin{equation} \boldsymbol{F}_{transl} = \tfrac{1}{2} \rho \omega^{2} S_{yy} \boldsymbol{C}_{F_{transl}}(\alpha), \end{equation}
where 
$\rho$ is the air density, 
$\omega$ is the stroke rate of the wing, 
$\boldsymbol {C}_{F_{transl}}$ is the translational-force coefficient vector as a function of the angle of attack 
$\alpha$ and 
$S_{yy}$ is the spanwise second moment of area of the wing. It expresses the scaling of translational forces with wing geometry, and is defined as
\begin{equation} S_{yy} = \int_{0}^{R} c(y) y^{2} \, {{\rm d} y} , \end{equation}
where 
$y$ is the position along the spanwise axis, 
$c$ is the wing chord length as a function of 
$y$ and 
$R$ is the span of the wing.
This translational-force model, when appropriately calibrated with empirical data, was found to be able to predict the lift forces of a hovering fruit fly rather well (Sane & Dickinson Reference Sane and Dickinson2001; Dickinson & Muijres Reference Dickinson and Muijres2016), but falls short in predicting the drag forces (Sane & Dickinson Reference Sane and Dickinson2001). The translational forces are most prominent during mid-stroke, where the stroke rate is highest and consequently also the stroke-rate-induced airflow over the wing is maximum (figure 1d). The translational aerodynamic force vector is often defined in terms of the lift and drag components 
$\boldsymbol {F} = (L, D)$, while the spanwise force is mostly ignored as it contributes little to powering flight (Dickinson et al. Reference Dickinson, Lehmann and Sane1999). Here, we will also ignore spanwise aerodynamic forces, but instead of using lift and drag components we will decompose all force vectors into the chord-wise wing-tangential force component and the force component normal to the wing surface 
$\boldsymbol {F} = (F_x, F_z)$ (figure 1f).
The pitch-rate-based forces are related to the pitching motion of the wing (i.e. rotation around the spanwise axis of the wing), and are a combination of a pitch–stroke interaction (Magnus-effect-based rotational lift) (Dickinson et al. Reference Dickinson, Lehmann and Sane1999) and a pure pitch-related effect (rotational drag) (Nakata et al. Reference Nakata, Liu and Bomphrey2015; Bomphrey et al. Reference Bomphrey, Nakata, Phillips and Walker2017). They are most prominent around stroke reversal, when wing-pitch rates are highest. These rotational forces have been modelled in a quasi-steady approach, based on pitch rate, stroke rate and wing geometry, in a previous study of the authors with a similar computational-fluid-dynamics-based methodology as used in the present investigation (van Veen et al. Reference van Veen, van Leeuwen and Muijres2019).
The wake-capture aerodynamic mechanism is most prominent at the start of each new stroke, after stroke reversal (Dickinson Reference Dickinson1994; Sane & Dickinson Reference Sane and Dickinson2001; Birch & Dickinson Reference Birch and Dickinson2003). This mechanism relates to the interaction of the wing with the wake generated in the previous stroke. Because the latter is dependent on the time history of the flow field, the wake-capture aerodynamic mechanism cannot be modelled based on the (instantaneous) wing kinematics and morphology alone. Despite this complication, still several studies have developed models for wake capture in both revolving and flapping wings (Usherwood & Ellington Reference Usherwood and Ellington2002; Nabawy & Crowther Reference Nabawy and Crowther2014; Oh et al. Reference Oh, Lee, Park, Choi and Kim2019).
The added-mass forces are caused by a change of the stroke rate, and are also referred to as acceleration-reaction forces or virtual-mass forces. They arise because when the wing accelerates, a volume of fluid around the wing (the ‘added mass’) is accelerated together with the wing. This accelerating fluid exerts a reaction force on the wing, which can be modelled as (Daniel Reference Daniel1984)
\begin{equation} \boldsymbol{F}_{AM} = \rho V \frac{{\rm d}}{{\rm d}t} (\boldsymbol{v} ), \end{equation}
where 
$V$ is the fluid volume that accelerates with the wing, 
$\boldsymbol {v}$ the velocity vector of the wing and 
$t$ is time. For insect flight, the added-mass forces are often approximated by the potential-flow model of a two-dimensional thin flat plate with an elliptical cross-section (Sedov Reference Sedov1965). For such a model, the forces normal to the wing surface depend on a cylindrical added volume with a cylinder diameter equal to the chord of the wing. The chord-wise tangential forces depend on a cylindrical added volume with a diameter equal to the thickness 
$\tau$ of the wing (Sedov Reference Sedov1965).
Based on these assumptions and using blade-element theory, the added mass for a revolving wing is modelled as
$$\begin{gather} \boldsymbol{F}_{AM} = (F_{x_{AM}}, F_{z_{AM}}) , \end{gather}$$
$$\begin{gather}F_{x_{AM}} = \frac{\rm \pi}{4} \rho \dot{\omega} S_{\tau y} \cos (\alpha) , \end{gather}$$
$$\begin{gather}F_{z_{AM}} = \frac{\rm \pi}{4} \rho \dot{\omega} S_{cy} \sin (\alpha) , \end{gather}$$
where 
$\dot {\omega }$ is the angular stroke acceleration of the wing; 
$F_{x_{AM}}$ and 
$F_{z_{AM}}$ are the chord-wise tangential and wing-normal added-mass forces, respectively (figure 1f). The wing geometry parameters 
$S_{\tau y}$ and 
$S_{cy}$ are the corresponding wing thickness-based second moment of area and chord-based second moment of area, respectively. These are defined as
$$\begin{gather} S_{\tau y} = \int_0^{R} \tau (y)^{2} y \,{{\rm d} y} , \end{gather}$$
$$\begin{gather}S_{cy} = \int_0^{R} c(y)^{2} y \,{{\rm d} y} , \end{gather}$$
where 
$\tau$ is the wing thickness.
Because the added-mass forces depend linearly on the angular acceleration of the wing, they are positive during the acceleration phase of the wing beat, and negative during the deceleration phase. As a result, wing-beat-average added-mass forces are often assumed to be close to zero, and are even exactly zero for a fully symmetrical motion (Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007; Zhang et al. Reference Zhang, Hedrick and Mittal2015). For this reason, added-mass forces are regularly ignored when studying the aerodynamics of insect flight (table 1) (Ellington Reference Ellington1984a; Dudley & Ellington Reference Dudley and Ellington1990; Dickinson & Muijres Reference Dickinson and Muijres2016).
Studies that included added-mass forces in their insect flight studies have used several variants of the potential-flow model (table 1). The majority of these studies assume a constant cylindrical added volume, according to the flat-plate potential-flow model (Sane & Dickinson Reference Sane and Dickinson2002; Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Wang et al. Reference Wang, Goosen and Van Keulen2016). This assumption was only tested in two studies by fitting an added-mass model with variable fluid volume to numerical simulation data (Lee et al. Reference Lee, Lua, Lim and Yeo2016; Cai et al. Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021). For both cases, the estimated added volume was found to be less than the cylindrical volume assumed according to potential theory ((1.4)–(1.6)).
Because the relative thickness of insect wings tends to be low, the wing thickness-based second moment of area is very small. On that ground, most studies ignore the tangential added-mass forces component (Sane & Dickinson Reference Sane and Dickinson2002; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Lee et al. Reference Lee, Lua, Lim and Yeo2016; Wang et al. Reference Wang, Goosen and Van Keulen2016), while only few include it for completeness (Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007).
As these tangential added-mass force models are based on inviscid flow theory, they only take the pressure forces into account, while viscous forces are ignored. Whether one can ignore the effect of viscous forces in added-mass models for low-Reynolds-number insect wings still remains to be tested. Viscous added-mass forces on insect wings have been studied in a different type of study (Zhang et al. Reference Zhang, Hedrick and Mittal2015). Using force decomposition and a panel method applied to flapping flight of a fruit fly and a hawk moth, they showed that viscosity-based centripetal added-mass forces significantly contribute to lift force production. This suggests that the quasi-steady added-mass model should be extended to include such a viscous force component.
The Wagner effect describes how the build-up of bound circulation around a wing starting from rest is delayed, with a corresponding delay in aerodynamic force production (Wagner Reference Wagner1925; Dickinson & Götz Reference Dickinson and Götz1996; Walker & Westneat Reference Walker and Westneat2000; Sane Reference Sane2003; Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015). This mechanism was first described by Wagner (Reference Wagner1925), who developed a general analysis of the potential flow that is generated by a wing starting from rest. In his study, he further elaborated two cases explicitly, namely the classical impulsively started wing as well as that of a wing undergoing a constant acceleration (velocity increasing linearly with time). Both cases similarly displayed the effect of delayed circulation build-up when compared with the quasi-steady value based on the instantaneous wing velocity.
The analysis of the Wagner effect was subsequently re-examined for wings operating at high angles of attack, as is the case for flying insects (Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015). This showed that the Wagner effect might significantly affect aerodynamic force production in flying insects. This effect is expected to be most prominent around stroke reversal, when wing accelerations are highest (Dudley & Ellington Reference Dudley and Ellington1990; Dickinson & Götz Reference Dickinson and Götz1996). Due to this, the Wagner effect might be especially relevant for high-frequency flappers such as mosquitoes. These animals have wing-stroke amplitudes of only several wing chords and thus the Wagner effect can significantly reduce aerodynamic forces throughout their complete wing stroke (Liu, Du & Sun Reference Liu, Du and Sun2020).
Because the Wagner effect depends on the time history of the wing movement, it is commonly expressed in terms of the distance that the wing has travelled from its initial position. Therefore, it does not lend itself easily to incorporation in quasi-steady modelling, as this type of modelling inherently relies on properties of the instantaneous motion. At least partly due to this complication, most studies on the aerodynamics of insect flight have ignored the Wagner effect in their force model (Ellington Reference Ellington1984a; Dudley & Ellington Reference Dudley and Ellington1990; Sane & Dickinson Reference Sane and Dickinson2002; Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007; Liu Reference Liu2009; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Dickinson & Muijres Reference Dickinson and Muijres2016; Wang et al. Reference Wang, Goosen and Van Keulen2016).
To our knowledge, only one application of a quasi-steady aerodynamic model has been reported that includes the Wagner effect (Wang et al. Reference Wang, Goosen and Van Keulen2016). Here, the Wagner effect was modelled as a direct correction on the translational lift, based on the classic Wagner function (Wang et al. Reference Wang, Goosen and Van Keulen2016). This approach proved reasonably successful in predicting aerodynamic forces for a wing starting from rest, but it was discarded for the analysis of flapping wings.
A more sophisticated analysis of the Wagner effect shows that it is, in essence, related to the inertia of the fluid around the accelerating wing (Wagner Reference Wagner1925; Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015). Therefore, we propose that the Wagner effect could potentially be modelled in a more general sense based on instantaneous wing accelerations, but as yet there exists no such model that could be applied to conditions of flapping (insect) wings. As one of the innovative contributions in the current investigation, we intend to fill this gap by proposing a quasi-steady model of the Wagner effect for flapping wings, in the form of a stroke-acceleration-dependent modification to the translational force on these wings (see methods § 2.1.2 for details).
All together, the individual aerodynamic mechanisms outlined above can be combined in an overall quasi-steady aerodynamic force model (Ellington Reference Ellington1984b). Such a model can then be used to describe the aerodynamic forces produced by a flying insect, based only on the motion and shape of its wings. Simplified versions of these models have been employed, which are restricted to a linear summation of the wing-velocity-based aerodynamic mechanisms, i.e. the stroke-rate-dependent translational-force mechanism and the pitch-rate-dependent rotational-lift mechanism (e.g. Dickinson & Muijres Reference Dickinson and Muijres2016). The aerodynamic mechanisms that depend on wing acceleration, the added-mass mechanism and Wagner effect, are often ignored, for different reasons. As mentioned above, the Wagner effect is difficult to capture in a quasi-steady model, while added-mass forces are often ignored because the wing-beat-average inviscid component reduces to zero for fully symmetric cyclic wing beats.
However, many insects tend to move their wings further to the back than to the front and the angle of attack is often different during the forward and backward wing strokes (figure 1a,b). As a result, added-mass forces do not fully cancel within a wing beat. Moreover, even during fully cyclic wing beats, added-mass forces still affect instantaneous force production, and as such influence the power requirement of flight (Ellington Reference Ellington1984a). Secondly, viscous added-mass forces also do not cancel during a cyclic wing beat (Zhang et al. Reference Zhang, Hedrick and Mittal2015).
In addition, because insects flap their wings at very high frequencies, the wing accelerations can become very high. For example, a hovering fruit fly (Drosophila hydei) flaps its wings at 200 Hz (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014), resulting in peak wing-stroke accelerations of two million radians per second squared (figure 1e). A similarly sized malaria mosquito (Anopheles coluzzii) flaps its wings at 600 Hz during hovering flight (Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a), resulting in peak wing-stroke accelerations of four million radians per second squared (figure 1e). These high wing accelerations suggest that acceleration-based aerodynamic mechanisms, such as the added mass and the Wagner effect, cannot be ignored when studying the aerodynamics of insect flight.
1.3. Modelling the stroke-related aerodynamic forces produced by accelerating wings
In this investigation, we studied how revolving wing-stroke accelerations affect aerodynamic forces on the wings of average-sized insects, operating at Reynolds numbers of the order of 100 (Dickinson et al. Reference Dickinson, Lehmann and Sane1999). As the particular focus here is on the wing-stroke-related effects, we considered wings that perform only a pure (revolving) wing-stroke motion, without rotational motions around the wing-pitch axis. In this way, we did consider variations in wing-pitch angle (i.e. angle of attack), but not wing-pitch rates or pitch accelerations. We developed and tested novel models that describe the acceleration-based aerodynamic mechanisms, i.e. added mass and Wagner effects (table 1).
To achieve this, we used computational-fluid-dynamics (CFD) simulations to systematically explore, over a wide parameter space, how the aerodynamic forces on a rigid wing vary with wing shape, angle of attack, stroke rate and stroke acceleration. We subsequently compared the architecture and parameters of our newly developed aerodynamic model with those of previously developed models (table 1).
Finally, we applied our aerodynamic model for accelerating wings to the wing-beat kinematics of a hovering fruit fly (Drosophila hydei) and a malaria mosquito (Anopheles coluzzii). These insects function here as model species for average-sized Dipteran insects, where the fruit fly represents insects with relatively low wing-beat frequencies (Dickinson & Muijres Reference Dickinson and Muijres2016), and the malaria mosquito represents insects with exceptionally high wing-beat frequencies (Bomphrey et al. Reference Bomphrey, Nakata, Phillips and Walker2017; Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a; Liu et al. Reference Liu, Du and Sun2020).
This study contributes to the understanding of the highly unsteady fluid dynamics of rapidly oscillating and accelerating wings, and thereby provides further insight into the evolutionary trade-offs in the flight apparatus of natural flyers. Furthermore, the study may assist in improving the flapping flight motor systems of bio-inspired flying robots (Wang et al. Reference Wang, Goosen and Van Keulen2016; Karasek et al. Reference Karasek, Muijres, De Wager, Remes and Croon2018).
2. Material and methods
2.1. A generic model for wing-stroke-based aerodynamic force production in insect wings
The main goal of our study is to improve our understanding of the wing-stroke-related aerodynamic mechanisms in flapping insect flight. In particular, we investigate systematically how stroke accelerations affect aerodynamic force production, and how this can be interpreted through the unsteady aerodynamic mechanisms of added mass and the Wagner effect. For this reason, we exclude the effect of wing-pitch rates, which has been given separate attention in a similar numerical study (van Veen et al. Reference van Veen, van Leeuwen and Muijres2019). We did, however, take variations in wing-pitch angle (i.e. angle of attack) into account. Accordingly, we analysed how the variations of wing-stroke rate 
$\omega$, wing-stroke acceleration 
$\dot {\omega }$, angle of attack 
$\alpha$ and wing shape affect the aerodynamic forces on a wing. To differentiate between the various force mechanisms, we introduce for convenience a nomenclature linked to the particular aerodynamic mechanism, and state their dependence on the kinematics and geometry parameters in brackets. For example, we will refer to the aerodynamic forces caused by the translational air movements over the wing as translational forces, and denote them as 
$\boldsymbol {F}_{transl}(\alpha, \omega, S_{yy} )$, indicating that they depend on the angle of attack 
$\alpha$, the stroke rate 
$\omega$ and the second moment of area 
$S_{yy}$ of the wing.
Adopting this notation, we define a generic model for the wing-stroke-related aerodynamic force production as
\begin{equation} \boldsymbol{F}_{total} = \boldsymbol{F}_{transl} (\alpha, \omega, S_{yy} ) + \boldsymbol{F}_{AM} (\alpha, \dot{\omega}, S_{cy} ) + \boldsymbol{F}_{WE} . \end{equation}
This model consists of three terms defined together with their parametric dependencies. These are the translational forces (
$\boldsymbol {F}_{transl}$), the added-mass-based forces (
$\boldsymbol {F}_{AM}$) and the Wagner-effect-based forces (
$\boldsymbol {F}_{WE}$). The models for the translational forces and added-mass-based forces are analogous to those defined in (1.1) and (1.5 and 1.6), respectively, while for the Wagner effect we will develop a new model.
2.1.1. The added-mass force model
In our model, we include added-mass forces both normal and tangential to the accelerating wing surface, as was also done by Andersen et al. (Reference Andersen, Pesavento and Wang2005) and Berman & Wang (Reference Berman and Wang2007). Compared with these previous approaches, we propose two important modifications:
(i) Firstly, we do not assume a priori that the added volume around the accelerating wing necessarily equals that of a cylinder, as suggested by potential-flow theory. Instead, we incorporate an adjustable added-mass force coefficient into our model that allows us to independently estimate the added fluid volume from the numerical simulation data. This approach is similar to that used by Lee et al. (Reference Lee, Lua, Lim and Yeo2016) and Cai et al. (Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021), for only the wing-normal forces.
(ii) Secondly, we do not assume that tangential added-mass forces are governed by pressure forces, in accordance with potential-flow theory, as done in all previous models (Sedov Reference Sedov1965; Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007). Instead, we propose that viscous forces are significant, if not dominant, for the tangential added-mass force dynamics, as is also the case for the translational forces (Sane & Dickinson Reference Sane and Dickinson2002). Given this assumption, the tangential forces should not scale with the wing thickness, but instead with wing chord, similar as for the wing-normal added-mass forces.
The resulting added-mass force model, comprising both the wing-normal and -tangential components, therefore becomes
\begin{equation} \boldsymbol{F}_{AM} = \rho \dot{\omega} S_{cy} \boldsymbol{C}_{F_{AM}}(\alpha) , \end{equation}
where 
$\boldsymbol {C}_{F_{AM}} = (C_{Fx_{AM}}, C_{Fz_{AM}})$ is the added-mass force coefficient vector with its wing-tangential and wing-normal components.
2.1.2. The Wagner-effect-based force model
In contrast to the models for translational forces and added mass, it is as yet unclear how the Wagner-effect-based force can be modelled based on wing kinematics and geometry. Wang et al. (Reference Wang, Goosen and Van Keulen2016) applied it as a Wagner-effect-based reduction of the translational forces using the Wagner function (Jones Reference Jones1940). This model predicted the aerodynamic forces on a impulsively starting wing well, but it did not apply well to flapping (insect) wings (Wang et al. Reference Wang, Goosen and Van Keulen2016). One important reason for this is that the Wagner function is expressed in terms of the distance travelled after start-up, and moreover depends on the nature of the acceleration. Both aspects are problematic for any quasi-steady model of a flapping wing that continuously accelerates or decelerates.
Detailed analyses of the Wagner effect on impulsively accelerating wings show that the Wagner solution comprises three terms, representing the translational lift, the added mass and a third term that represents the non-equilibrium effect of the acceleration on the circulatory lift (Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015). This third term has since then been named after Wagner as the Wagner effect.
In view of this, we here propose to model the Wagner effect of accelerating insect wings using a similar physics-based approach. Here, the Wagner effect thus reflects the reduction in instantaneous translational-force production on an accelerating wing, as a result of a delay in the build-up of bound circulation. A valuable insight, in this respect, is that the Wagner effect can therefore essentially be interpreted as resulting from a specific non-equilibrium condition, where the instantaneous translational force deviates from its equilibrium value corresponding to the instantaneous wing-stroke rate. Accordingly, when considering the case of an accelerating wing, we propose to interpret the interaction between stroke rate and stroke acceleration as a generalization of the Wagner effect, and consequently model it as an acceleration-based modification of the translational forces.
To determine how this modification might scale with wing kinematics and geometry, we introduce a non-dimensional measure that expresses the relative importance of the stroke acceleration with respect to the stroke rate. For this purpose, we define the ratio between the time scales of translational air movement across the wing surface (
$T_{transl}$) and acceleration of the moving air (
$T_{accel}$). For a linearly translating wing this ratio reads
\begin{equation} \frac{T_{transl}} {T_{accel}} = \frac{c/v}{\sqrt{c/\dot{v}}} = \frac{\sqrt{c\dot{v}}}{v} , \end{equation}
where 
$v$ is the wing speed, 
$\dot {v}$ the wing acceleration and 
$c$ the wing chord. As the Wagner effect increases with the relative acceleration effect, we can likewise expect that its strength will increase with the magnitude of the above parameter. We tentatively propose, furthermore, that this relation is (approximately) linear, which is an assumption that will be justified by the final results upon fitting this model to the simulation data base. Based on these considerations, we model the acceleration-related modification of the translational forces as a result of the Wagner effect as
\begin{equation} \boldsymbol{F'}_{WE} = \frac{1}{2} \rho v^{2} c \left(\frac{\sqrt{c \dot{v}}}{v}\right) \boldsymbol{C}_{F_{WE}}(\alpha) . \end{equation}
Here, 
$\boldsymbol {F'}_{WE}$ is the spanwise-sectional Wagner-effect-based translational-force modification, and 
$\boldsymbol {C}_{F_{WE}}(\alpha )$ is the corresponding force coefficient. When converting to the reference frame of a revolving wing and integrating across the span, the Wagner-effect model becomes
\begin{equation} \boldsymbol{F}_{WE} = \tfrac{1}{2} \rho \omega \sqrt{\dot{\omega}} S_{WE} \boldsymbol{C}_{F_{WE}}(\alpha) . \end{equation}
Here, 
$S_{WE}$ is the wing geometry scaling parameter of the Wagner effect, which follows from the preceding analysis as
\begin{equation} S_{WE} = \int_0^{R} \sqrt{c(y)^{3} y^{3}} \,{{\rm d} y} . \end{equation}Note that the here-proposed Wagner-effect model is a linear approximation of a complex aerodynamic mechanism that is inherently nonlinear (Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015). However, this simplification allows one to apply the model to flapping insect wing kinematics using a quasi-steady approach.
2.1.3. Our newly proposed aerodynamic force model for insect wings with stroke accelerations
In conclusion, the complete model for the stroke-based aerodynamic force produced by an accelerating wing is summarized as
\begin{equation} \boldsymbol{F}_{total} = \boldsymbol{F}_{transl} (\alpha, \omega, S_{yy} ) + \boldsymbol{F}_{AM} (\alpha, \dot{\omega}, S_{cy} ) + \boldsymbol{F}_{WE} (\alpha, \omega, \dot{\omega}, S_{WE}) , \end{equation}where the translational forces, added-mass forces and Wagner-effect-based forces scale with wing geometry and kinematics as
$$\begin{gather} \boldsymbol{F}_{transl} = 1/2 \rho {\omega}^{2} {S}_{yy} \boldsymbol{C}_{F_{transl}}(\alpha) , \end{gather}$$
$$\begin{gather}\boldsymbol{F}_{AM} = \rho \dot{\omega} S_{cy} \boldsymbol{C}_{F_{AM}}(\alpha) , \end{gather}$$
$$\begin{gather}\boldsymbol{F}_{WE} = 1/2 \rho \omega \sqrt{\dot{\omega}} S_{WE} \boldsymbol{C}_{F_{WE}}(\alpha) . \end{gather}$$
Here, both the added-mass forces and Wagner-effect-based aerodynamic forces depend on wing stroke acceleration, and therefore we will refer to them combined as the acceleration-based forces (
$\boldsymbol {F}_{accel} = \boldsymbol {F}_{AM} + \boldsymbol {F}_{WE}$).
This model is different from previous aerodynamic models for insect flight in that it combines an extended added-mass model, including variable added volume and a viscous-flow-based tangential force component, with a novel Wagner-effect model that can be applied using the quasi-steady modelling approach (table 1).
It may be recalled here that, in view of the scope of the current study, our analysis considers only the effect of the revolving stroke motion on aerodynamic forces, and thus we have left the effect of wing-pitch rate and pitch accelerations on aerodynamic force production out of consideration. We did, however, model how aerodynamic forces vary with the wing-pitch angle (i.e. angle of attack).
We tested the validity of this newly proposed model for wing-stroke-based aerodynamic force production, by performing a model fit to an extensive series of numerical simulations with systematically varying wing-beat kinematics and wing geometry. The analysis consisted of the following consecutive steps:
(i) We first used the forces from our numerical simulations of accelerating fruit fly wings in a model-fitting procedure to determine the force coefficients
$\boldsymbol {C}_{F_{transl}}$, $\boldsymbol {C}_{F_{AM}}$ and $\boldsymbol {C}_{F_{WE}}$, and how they vary with the angle of attack $\alpha$. Here, we also tested the validity of our Wagner-effect model (2.5) and our added-mass force model with viscous-flow-based tangential forces (2.2).(ii) Secondly, we used the simulation data of accelerating and decelerating fruit fly wings to test how well the model captures the acceleration-based aerodynamic forces during both wing-stroke accelerations and decelerations.
(iii) Thirdly, we tested how well the model predicts forces produced by wings with variable wing-planform geometries;
Following this, we studied the underlying aerodynamics of stroke-acceleration-based force production using the CFD flow visualizations. Then, using a quasi-steady implementation of our model, we investigated how stroke-acceleration-based forces affect flight force production of hovering fruit flies and malaria mosquitoes, as model species for low-frequency and high-frequency flapping flyers, respectively. Finally, in the discussion section, we compared the architecture and parameters of our model with those of quasi-steady aerodynamic models previously reported in the literature (table 1).
2.2. Numerical set-up
All CFD simulations were conducted with the immersed-boundary adaptive-mesh-refinement viscous-flow solver IBAMR (Bhalla et al. Reference Bhalla, Bale, Griffith and Patankar2013). A single rigid wing with a span of approximately 
$R = 3$ mm and a thickness of 
$\tau = 18\ \mathrm {\mu } {\rm m}$ was placed with its root in the centre of a domain with a size of 
$5\ {\rm cm}\times 5\ {\rm cm}\times 5\ {\rm cm}$. A periodic boundary condition was used on the sides of the domain, and a no-slip boundary condition was used at the surfaces of the wing. We determined the aerodynamic forces acting on the wing directly from the constraint force term of the no-slip boundary condition at the wing surfaces (Bhalla et al. Reference Bhalla, Bale, Griffith and Patankar2013).
The solver has been tested and validated in a previous study by comparing the results of a simulation of the hovering kinematics of a fruit fly with forces measured in an experiment with a robotic flapper (van Veen et al. Reference van Veen, van Leeuwen and Muijres2019). This validation approach allowed us to test the model accuracy during the complete range of wing rotation rates and accelerations relevant for flying insects (figure 1). Based on this validation study, we chose the temporal and baseline spatial resolutions as 
$\Delta t = 1\times 10^{-7}$ s and 
$\Delta x =0.01$ mm, respectively. The solver uses adaptive mesh refinement based on the vorticity in the flow field; in total, we used three mesh refinement levels with vorticity thresholds at 50, 500, 5000 s
$^{-1}$ from coarsest to finest, respectively. This resulted in a mesh size of approximately 4 million cells and a total of 10 000 time steps per simulation.
2.3. Wing geometries
We used 14 different wing shapes to systematically study the effect of wing geometry on the stroke-based forces (2.7). Each wing was modelled as a rigid plate with a thickness of 18 
$\mathrm {\mu }$m, defined as a set of approximately 100 000 equidistantly arranged Lagrangian boundary points. Two wing planforms were based on the wing morphology of a fruit fly and malaria mosquito, respectively. The other 12 were elliptical wings with systematically varying geometries, consisting of three sets (figure 2a–e):
(i) The first set of four elliptical wings had a fixed wing span
$R =3$ mm and a varying maximum chord length $c_{max} = 0.75$ to 1.5 mm with steps of 0.25 mm (figure 2d, squares).(ii) The second set had a fixed maximum chord length
$c_{max} = 1$ mm and a varying wing span $R = 2$ to 4 mm in steps of 0.5 mm (figure 2d, circles).(iii) The third set of elliptical wings were designed such that the chord-based second moment of area
$S_{cy}$ was kept constant (figure 2e, diamonds). This set of wings allows testing of how added-mass force production depends on wing geometry. For elliptical wings, the geometric scaling parameter of the stroke-acceleration model as described in (1.4) is defined as (Sedov Reference Sedov1965)
(2.11)where\begin{align} S_{cy} = \int_0^{R} c(y)^{2} y \, {{\rm d} y} = \int_{-({1}/{2}){\rm \pi}}^{({1}/{2}) {\rm \pi}} (c_{max} \cos(\varTheta))^{2} \left(\frac{1}{2} R\sin(\varTheta) + \frac{1}{2}R \right) \frac{1}{2}R \, {\rm d}\varTheta , \end{align}$\varTheta$ is the polar coordinate defined in figure 2c. This shows that, for elliptical-shaped wings, $S_{cy}$ remains the same if the product $c_{max}R$ is kept constant. Therefore, for the third set of elliptical wings, we increased $c_{max}$ in four steps from 0.75 to 1.5 mm, while simultaneously decreasing $R$ from 4 to 2 mm (figure 2d, diamonds)).
Description of the investigated parameter space, including variations in wing-planform geometry and wing kinematics. (a–c) Geometry of the fruit fly wing, malaria mosquito wing and ellipse-shaped wings, respectively. (d,e) Variations in geometry of all tested wings, including variations in maximum chord length vs wing span (d), and the spanwise second moment of area 
$S_{yy}$ vs the chord-based second moment of area 
$S_{cy}$ (e). ( f–h) The kinematics parameter space as defined by variations in angle of attack, stroke rate and stroke acceleration. Each dot highlights a set of simulations at those values. Here, blue dots indicate the simulations with a fruit fly wing, red dots indicate the simulations with the fly wing, mosquito wing and wings with constant 
$S_{cy}$, and green dots shows simulations that are performed on all wing geometries (table 2). The box highlights the set of simulation that we use as an example in § 2.5. The blue and red surfaces show the kinematics spanned by a hovering fruit fly (blue) (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014) and malaria mosquito (red) (Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a).
The values of maximum chord length, wing span, 
$S_{yy}$ and 
$S_{cy}$ of both the fruit fly and the malaria mosquito fall within the parameter space spanned by the sets of elliptical wings (figure 2d,e).
2.4. Reference frames
The position and orientation of the wing and the forces computed by the numerical solver were expressed in a right-handed world reference frame with its origin at the root of the wing (figure 1f). The orientation of the wing was expressed in this reference frame using Euler angles by first rotating around the 
$z_{world}$-axis of the world reference frame with the stroke angle 
$\gamma$, and then by rotating it around its spanwise axis (y-axis of the wing reference frame) with a wing-pitch angle (
$\phi$). Because in the simulations there is no difference between the pitch angle and the angle of attack, we will always refer to this angle as the angle of attack (
$\alpha$). If there were an additional wing rotation (i.e. a deviation motion component), or a non-zero free-stream velocity field, the pitch angle and the angle of attack would not coincide.
The aerodynamic forces on the wings are expressed in the wing reference frame, which is a right-handed coordinate frame with the origin at the wing hinge location, the x-axis parallel to the wing surface pointing towards the trailing edge, the y-axis parallel to the surface pointing towards the wing tip and the z-axis perpendicular to the wing surface (figure 1f).
2.5. Kinematics set-up
We defined a set of simulations that allowed us to isolate the effect of rotational stroke acceleration on aerodynamic force production from the other aerodynamic mechanisms. Each set of wing kinematics was controlled by three main parameters: stroke rate, stroke acceleration and angle of attack. The angle of attack was kept constant throughout each simulation, but was varied between simulations. The stroke rate and stroke acceleration varied throughout the simulation. The wing was not given an additional translational motion, and as such we simulated the rotational wing stroke motion as occurring in hovering flight. We extracted the aerodynamic forces for our analysis at the last time step of each simulation, and accordingly classified each simulation with the combination of stroke rate and stroke acceleration at the same end of the simulation.
To explain how the wing kinematics study was designed, we consider as an example sub-set of all simulations (figure 1g,h). This sub-set consists of simulations with the fruit fly wing operating at a single (constant) angle of attack 
$\alpha = 36^{\circ }$, a single stroke rate at the end of the simulation (
$\omega = 1000$ rad s
$^{-1}$) and a range of simulation-end stroke accelerations, varying from 
$\dot {\omega } = -4.5\times 10^{6}$ to 
$\dot {\omega } = 4.5\times 10^{6}$ rad s
$^{-2}$. These simulations can be further subdivided into three types: simulations with accelerating wings (
$\dot {\omega } > 0$ rad s
$^{-2}$), decelerating wings (
$\dot {\omega } < 0$ rad s
$^{-2}$) and a simulation with a non-accelerating wing (
$\dot {\omega } = 0$ rad s
$^{-2}$). All simulations of this sub-set are highlighted by the black squares in figure 2( f–h), while the corresponding kinematics is shown in figure 1(g,h).
All wings started at a zero stroke rate and smoothly accelerated to a constant baseline stroke-rate value during the first 0.5 ms of the simulation (not shown in figure 1g,h). For the non-accelerating wing, this baseline value was equal to the final stroke rate of 
$\omega = 1000$ rad s
$^{-1}$, which was then kept constant during the rest of the simulation (grey data in figure 1g,h). For the wings with positive stroke accelerations, the baseline value was equal to 60 % of the stroke rate at the end of the simulation (
$\omega = 600$ rad s
$^{-1}$; green data in figure 1g). For the decelerating wings, the baseline value was equal to 140 % of the stroke rate at the end of the simulation (
$\omega = 1400$ rad s
$^{-1}$; purple data in figure 1g). For all accelerating and decelerating wings, we then smoothly changed the stroke acceleration from zero to the required stroke acceleration, which was then kept constant until the end of the simulation (figure 1h).
To control for the effect of time on the development of the flow, we designed the kinematics such that each simulation covered exactly the same duration (1 mm). The time it takes to accelerate or decelerate to the required stroke rate at the end of the simulation varies with the stroke-acceleration magnitude for that simulation. Therefore, the acceleration and deceleration of each wing were started at different moments in time, such that the required stroke rate (in this case, 
$\omega = 1000$ rad s
$^{-1}$) was always reached at 
$t = 1$ mm (figure 1g,h).
At the end of all the simulations in this set, the stroke rate and angle of attack were the same (figure 1g,h); only the stroke acceleration was different between the cases. This means that only the stroke-acceleration parameter is responsible for the variations in aerodynamic forces between the simulations.
2.6. Kinematics parameter space
The kinematics described in § 2.5 has three kinematic input parameters: stroke rate, stroke acceleration and angle of attack. In this study, these parameters were systematically varied to assess their influence on aerodynamic force production (figure 2f–h). The range of these parameters was determined based on the typical range encountered during hovering flight of a fruit fly (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014) and a malaria mosquito (Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a) (figure 2f–h, blue and red data, respectively). All these simulations can roughly be divided into two sets: a set of simulations with the fruit fly wing, and a set of simulations with all wing geometries (table 2).
Overview of the range of kinematics of all simulations with the fruit fly wing (Drosophila hydei, D.h.), the malaria mosquito wing (Anopheles coluzzii, An.c.) and the elliptic wings with constant chord-based second moment of area (
$\text {e}_{Scy}$), constant wing span (
$\text {e}_R$) and constant maximum chord length (
$\text {e}_c$). For each parameter, we give range and step size (
$\Delta$).
The set of simulations with the fruit fly wing was primarily used to quantify in detail the relation between the wing kinematics and the resulting total aerodynamic force 
$F_{total}$, as defined in (2.7). For this, we performed 568 fruit fly wing simulations, consisting of 40 simulations with non-accelerating wings, 480 simulations with accelerating wings and 48 simulations with decelerating wings (blue dots in figure 2f–h, table 2):
(i) The simulations with non-accelerating fruit fly wings were performed at a stroke-rate range of 250–1750 rad s
$^{-1}$ with steps of 500 rad s$^{-1}$, and angles of attack ranging from 9$^{\circ }$ to 90$^{\circ }$ in steps of 9$^{\circ }$. At the angle of attack of 36$^{\circ }$, we performed simulations at a stroke-rate step size of 250 rad s$^{-1}$.(ii) The accelerating fruit fly wing simulations were performed at a stroke-acceleration range of
$1\times 10^{6}$ to $4.5\times 10^{6}$ rad s$^{-2}$ with steps of $0.5\times 10^{6}$ rad s$^{-2}$, a stroke-rate range of 500–1750 rad s$^{-1}$ with steps of 250 rad s$^{-1}$ and a range of angles of attack of 9$^{\circ }$ to 90$^{\circ }$ in steps of 9$^{\circ }$.(iii) The decelerating fruit fly wing simulations were performed at an angle of attack of 36
$^{\circ }$, a range of stroke accelerations of $-1\times 10^{6}$ to $-4.5\times 10^{6}$ rad s$^{-2}$ with steps of $0.5\times 10^{6}$ rad s$^{-2}$ and a stroke-rate range of 500 to 1750 rad s$^{-1}$ with steps of 250 rad s$^{-1}$.
The set of simulations with all wing types was primarily used to quantify the dependence of the total stroke-based forces 
$F_{total}$ on wing-planform geometry. For this, we ran 153 additional simulations with the mosquito wing and the 12 elliptical wings (figure 2). All of these simulations were performed at a constant angle of attack 
$\alpha = 36^{\circ }$, at both positive stroke accelerations and zero acceleration (green and red dots in figure 2f–h, table 2).
The simulations at zero stroke acceleration were performed on all thirteen additional wings, and at stroke rates of 
$\omega = 500$, 1000 and 1500 rad s
$^{-1}$. This resulted in 39 additional simulations at zero stroke acceleration.
The simulations with positive stroke accelerations consisted of two sub-sets:
(i) With the mosquito wing and the five elliptical wings with constant
$S_{cy}$ we performed simulations at stroke accelerations of $1\times 10^{6}$ and $4.5\times 10^{6}$ rad s$^{-2}$, and at stroke rates of $\omega = 500$ to 1750 rad s$^{-1}$ in steps of 250 rad s$^{-1}$. This resulted in 72 simulations with these six wings.(ii) For the seven elliptical wings with varying wing span and varying chord length, we performed simulations at stroke accelerations of
$1\times 10^{6}$ and $4.5\times 10^{6}$ rad s$^{-2}$, and stroke rates of 500, 1000 and 1500 rad s$^{-1}$. This resulted in 42 simulations with these seven additional wings.
The 721 CFD simulations described above (figure 2, table 2) were used to systematically study how wing-stroke accelerations and wing geometry affect aerodynamic force production, and to study the underlying aerodynamics (see § 2.7 below). The stroke rate and stroke accelerations at the end of these simulations ranged from 
${\omega } = 250$ to 1750 rad s
$^{-1}$ and 
$\dot {\omega } =-4.5\times 10^{6}$ to 
$4.5\times 10^{6}$ rad s
$^{-1}$, respectively. For these stroke rates, the average tested wing was operating at a mid-span Reynolds-number range of 36–252. This range is equivalent to the Reynolds numbers of flapping wings of average-sized flying insects (Dickinson et al. Reference Dickinson, Lehmann and Sane1999). The stroke angle travelled for all simulations ranged from 
$18^{\circ}$ to 130
$^{\circ }$ for the combination of lowest stroke rate and highest stroke acceleration, and highest stroke rate and lowest stroke acceleration, respectively. Thus none of the wings moved through their own wake, avoiding wake capture.
2.7. Aerodynamics of acceleration-based force production
We studied the aerodynamics of the acceleration-based force production, by extracting the airflow dynamics that can be associated with the added-mass forces 
$F_{AM}$ and Wagner-effect forces 
$F_{WE}$ (cf. (2.7)–(2.10)). Note that this cannot be obtained directly from the simulations, as for an accelerating wing the three effects (translational, added-mass and Wagner effect) are always occurring in combination. Hence, a special data processing procedure was followed to separate them.
We quantified the aerodynamics using pressure distributions on the top and bottom surfaces of the wing, as well as in a plane perpendicular to the span of the wing (figure 7a). We extracted the pressure fields on the wing surface from surfaces 0.06 mm above and below the wing surface. This offset between wing surface and pressure field was needed because the cell height of the numerical mesh was of the same order of magnitude as the wing thickness; for our numerical set-up, the offset of 0.06 mm was proven to result in the minimal numerical error (van Veen et al. Reference van Veen, van Leeuwen and Muijres2019). Note that we used these pressure fields to study the aerodynamics of the accelerating wings only. The actual forces acting on the wing were calculated from the constraint force term of the immersed-boundary method functions (Bhalla et al. Reference Bhalla, Bale, Griffith and Patankar2013).
2.7.1. Aerodynamics of the added-mass forces
We approximated the pressure distribution around an accelerating wing that is due to the added-mass mechanism, following three consecutive steps:
(i) We first extracted from the simulation the (instantaneous) pressure distribution
$p_{total}$ around the wing, for the case of a given stroke acceleration ($\dot {\omega }$) and the lowest stroke rate ($\omega = 500$ rad s$^{-1}$); this represents the total force $F_{total}$ in (2.7).(ii) Secondly, we extracted the equivalent pressure distribution
$p_{transl}$ around the same wing when moving at the same stroke rate ($\omega = 500$ rad s$^{-1}$), but without a stroke acceleration ($\dot {\omega } = 0$ rad s$^{-2}$). This pressure distribution therefore corresponds to the translational-based aerodynamic force $F_{transl}$.(iii) We then estimated the pressure distribution corresponding to the added-mass force at the given stroke acceleration as
(2.12)\begin{equation} p_{AM} \approx p_{total} - p_{transl} . \end{equation}
Here, we assume that, at this lowest stroke rate, the Wagner effect on this pressure distribution is limited, and that thus the extracted pressure dynamics is primarily caused by the added-mass effect. Using this approach, we systematically studied how the pressure varied with changes in wing stroke acceleration, angle of attack and wing geometry.
We also used these pressure distributions to estimate the equivalent thickness of the fluid layer that is accelerated with the wing (
$\delta _{AM}$). Assuming that the 
$p_{AM}$ distribution on the wing surface can be associated with the local acceleration of a fluid volume of thickness 
$\delta _{AM}$, this implies that
\begin{equation} \delta_{AM}(x,y) = \frac{p_{AM}(x,y)}{\rho \dot{\omega} y} . \end{equation}We determined this added-mass fluid-layer distribution for the set of elliptical wings with variable chord length, in order to test how the fluid-layer thickness varies along the span for this range of wings. Based on this, we tested how the resulting spanwise distribution of the added volume deviates from that of the previously assumed cylindrical added-volume model (1.4) (Sedov Reference Sedov1965). For the visualization of the spanwise fluid layer thickness distribution we assumed an elliptical chord-wise distribution.
2.7.2. Aerodynamics of the Wagner-effect forces
We estimated the pressure distribution around an accelerating wing due to the Wagner effect 
$p_{WE}$ as the difference in pressure distribution between two wings moving at the same stroke acceleration 
$\dot {\omega }$, but with different stroke rates 
$\omega$ as
\begin{align} p_{WE} \approx [p_{total} (\dot{\omega}, \omega_1) - p_{transl} (\dot{\omega}=0, \omega_1 )] - [p_{total} (\dot{\omega}, \omega_2) - p_{transl} (\dot{\omega}=0, \omega_2 )] , \end{align}
where the stroke-rate difference 
$\Delta \omega = \omega _1 - \omega _2$ is kept constant. Thus, the Wagner-effect- based pressure 
$p_{WE}$ is obtained as the difference in the acceleration-based pressure (
$p_{total} - p_{transl}$) between two wings moving with a stroke-rate difference of 
$\Delta \omega$. Note that, as the acceleration is the same for both wings, the added-mass contributions cancel. The result should therefore correspond to the pressure distribution that causes the Wagner-effect-based forces for a stroke rate equal to 
$\Delta \omega$ and stroke acceleration 
$\dot {\omega }$ (2.10). Using this approach, we determined the 
$p_{WE}$ distribution on the top and bottom surfaces of a fruit fly wing, operating at various stroke-rate differences, with and without a stroke acceleration or stroke deceleration, and at a range of angles of attack.
Next to this, we studied how the Wagner effect affects the vortex system around the wing, and particularly how it affects the LEV strength. For this, we compared the spanwise vorticity around fruit fly wings that moved at various levels of wing-stroke acceleration. This analysis consisted of two sets of tests:
(i) Firstly, we tested how the LEV strength differed between accelerating wings and non-accelerating wings, operating at a range of angles of attack (
$\alpha =27^{\circ }$, 45$^{\circ }$, 63$^{\circ }$ and 81$^{\circ }$). Here, all wings were moving at a stroke rate of $\omega = 1250$ rad s$^{-1}$, and the accelerating wings had a stroke acceleration of $\dot {\omega } = 4.5\times 10^{6}$ rad s$^{-2}$.(ii) Secondly, we compared the spanwise vortex system around accelerating and decelerating wings, at the stroke accelerations
$\dot {\omega } = -4\times 10^{6}$, $-2\times 10^{6}$, 0, $2\times 10^{6}$ and $4\times 10^{6}$ rad s$^{-2}$. All these wings were operating at the same stroke rate $\omega = 750$ rad s$^{-1}$ and angle of attack $\alpha = 36^{\circ }$, and thus only the stroke accelerations differed between the cases.
2.8. Applying our new accelerating wing models to the wing beat of flying insects
To test how wing-stroke accelerations affect lift and drag forces during the wing beat of flying insects, we formulated a quasi-steady version of our basic aerodynamic force model for accelerating wings ((2.7)–(2.10)). We then applied this quasi-steady model to the actual wing-stroke kinematics of a hovering fruit fly (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014) and a hovering malaria mosquito (Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a). This allowed us to test how the different aerodynamic mechanisms (translational, added mass and Wagner effect) affect lift and drag forces during the forward and backward wing stroke of these flying animals. Note that we did not include a (wing-pitch-rate-related) rotational-lift mechanisms in our model, and thus we here only test how added-mass forces and the Wagner-effect-based forces compare with the translational forces produced by a flapping insect wing. For more details about the quasi-steady application of the new model, see supplementary material S1.1 available at https://doi.org/10.1017/jfm.2022.31.
3. Results
The results of the numerical simulations of accelerating wings are processed and analysed according to the following steps:
(i) In § 3.1, we analyse how the aerodynamic forces depend on wing kinematics and morphology, and use these data to test the validity of our proposed new aerodynamic force model ((2.7)–(2.10)).
(ii) In § 3.2, we use the airflow visualization results to discuss the aerodynamics that underlies the observed stroke-acceleration-based force production.
(iii) Finally, in § 3.3, we apply the newly developed and validated force model to analyse the contribution of acceleration effects in the wing-beat kinematics of a hovering fruit fly and a hovering malaria mosquito.
3.1. Acceleration-based aerodynamic forces produced by insect wings
We used the CFD simulations with the fruit fly wing at zero and positive stroke accelerations to determine how stroke-acceleration-based aerodynamic forces vary with stroke rate, stroke acceleration and angle of attack (supplementary figures S1 and S2). For this, we focused on the chord-wise tangential forces and the forces acting normal to the wing, as the spanwise viscous forces contribute very little to lift and drag production. All forces, including the spanwise tangential forces, are provided in the supplementary figures S1 and S2 and supplementary database S1.
3.1.1. Normal and tangential aerodynamic forces produced by accelerating insect wings
By focusing on the normal and chord-wise tangential forces, our force-vector-based aerodynamic model can be converted into a set of scalar-valued models, one for the wing-normal forces (acting in the z-direction) and one for the chord-wise tangential forces (acting in the x-direction, figure 1f). These scalar-based models are defined as
$$\begin{gather} F_{z_{total}} = F_{z_{transl}} (\alpha, \omega, S_{yy} ) + F_{z_{AM}} (\alpha, \dot{\omega}, S_{cy} ) + F_{z_{WE}} (\alpha, \omega, \dot{\omega}, S_{WE}) \end{gather}$$
$$\begin{gather}F_{z_{transl}} = 1/2 \rho {\omega}^{2} {S}_{yy} C_{Fz_{transl}}(\alpha) \end{gather}$$
$$\begin{gather}F_{z_{AM}} = \rho \dot{\omega} {S}_{cy} C_{Fz_{AM}}(\alpha) \end{gather}$$
$$\begin{gather}F_{z_{WE}} = 1/2 \rho \omega \sqrt{\dot{\omega}} S_{WE} C_{Fz_{WE}}(\alpha) , \end{gather}$$and
$$\begin{gather} F_{x_{total}} = F_{x_{transl}} (\alpha, \omega, S_{yy} ) + F_{x_{AM}} (\alpha, \dot{\omega}, S_{cy} ) + F_{x_{WE}} (\alpha, \omega, \dot{\omega}, S_{WE}), \end{gather}$$
$$\begin{gather}F_{x_{transl}} = 1/2 \rho {\omega}^{2} {S}_{yy} C_{Fx_{transl}}(\alpha), \end{gather}$$
$$\begin{gather}F_{x_{AM}} = \rho \dot{\omega} S_{cy} C_{Fx_{AM}}(\alpha) , \end{gather}$$
$$\begin{gather}F_{x_{WE}} = 1/2 \rho \omega \sqrt{\dot{\omega}} S_{WE} C_{Fx_{WE}}(\alpha) , \end{gather}$$
where 
$F_z$ and 
$F_x$ are the wing-normal and chord-wise tangential force scalars, respectively; 
$C_{Fz}$ and 
$C_{Fx}$ are the corresponding normal and tangential force coefficients.
The normal and tangential aerodynamic forces produced by the fruit fly wing for all CFD simulations with zero or positive stroke accelerations are shown in figure 3. Here, all forces are normalized as 
$F^{*} = F / W$, with weight 
$W=10$ 
$\mathrm {\mu }$N as the approximate weight of our two model species, the fruit fly and malaria mosquito.
Weight-normalized total wing-normal forces and chord-wise tangential forces on an accelerating fruit fly wing as a function of the stroke acceleration 
$\dot {\omega }$, stroke rate 
$\omega$ and angle of attack 
$\alpha$. The left column (a–f) shows the normal forces (
$F_z$), and the right column (g–i) shows the equivalent chord-wise tangential forces (
$F_x$), at the same scale. Each panel shows the results for a single angle of attack (see title); each data point shows the results of a single simulation, as the normalized aerodynamic force (ordinate) for a given stroke acceleration (abscissa) and stroke rate (colour coded, see colour bar). The lines in each panel show the nonlinear least-squares fitting results at that angle of attack of the model defined in ((3.1)–(3.8)).
We grouped all results according to the separate panels in figure 3, and then fitted our aerodynamic model ((3.1)–(3.8)) to these results, for each combination of angle of attack and force direction separately. For each fit, we used a nonlinear least-squares fitting routine in Matlab (Mathworks Inc). The root mean square fitting errors for all fits are 
$0.02 \pm 0.01 \ (\text {mean} \pm \text {standard deviation},\ n = 20 \ \text {fits})$.
The normal forces are seen to increase in magnitude with both stroke rate and stroke accelerations, except for the cases with a combination of high stroke rates and low stroke accelerations. Here, the wing-normal forces slightly decrease with increasing stroke accelerations (figure 3a–j). The fitting results, as shown as curved lines in each panel of figure 3, capture the variations in aerodynamic force production well, including the above-mentioned nonlinear dynamics in pressure forces at high stroke rates and low stroke accelerations. The model is capable of capturing this because the Wagner-effect model includes a square root dependence on stroke acceleration.
From these fitting results at each angle of attack, we extracted the wing-normal and -tangential force coefficients (
$C_{Fz}$ and 
$C_{Fx}$, respectively) of all three force components in our model (figure 4a,b). All wing-normal force coefficients increase in magnitude with angle of attack, whereby the translational-force coefficients and added-mass coefficients have positive values and the Wagner-effect coefficients are negative (figure 4a). The positive added-mass force coefficient shows that added-mass forces increase with stroke accelerations, whereas the negative Wagner-effect force coefficients show that the Wagner effect causes a reduction in aerodynamic force produced when a wing accelerates. These trends are in line with the physical mechanism for each of these components, as discussed in § 1.2.
(a) Wing-normal pressure force coefficients of the translational-based force model, the added-mass force model and the Wagner-effect model (in blue, red and green, respectively) vs the angle of attack for an accelerating fruit fly wing. The data points show separate fits with their 95 % confidence interval, and the curves are least-squares sine fits through these data. (b) Wing chord-wise tangential viscous force coefficients of the translational-based force model, the added-mass force model and the Wagner-effect model (in blue, red and green, respectively) vs the angle of attack for an accelerating fruit fly wing. The data points show separate fits with their 
$95\,\%$ confidence interval, and the curves are least-squares fits through these data. (c) The tangential-to-normal force coefficient ratio of the translational-based force model, the added-mass force model and the Wagner-effect model (in blue, red and green, respectively) vs the angle of attack, determined from the data in panels a–b as 
$RC_{Fx-z} = C_{Fx} / C_{Fz}$. (d) The weight-normalized aerodynamic forces acting normal to a fruit fly wing with an angle of attack 
$\alpha = 36^{\circ }$, and that is both accelerating and decelerating. Colours indicate results at variable stroke rates (see colour bar). Each dot shows the results of a single CFD simulation and the lines show nonlinear least-squares fitting results of the model defined in ((3.1)–(3.4)). The coloured lines show the fit based on the complete data set, and the black lines show the fit based on the steady and accelerating wings only. (e) Wing-normal force coefficients of the translational-force model, added-mass force model and Wagner-effect model (in blue, red and green, respectively) from the nonlinear least-squares fits in panel b, for steady and accelerating wings only (left) and for the complete data set (right). For each data point, we show the fit and the 
$95\,\%$ confidence interval.
The tangential forces are small compared with the normal forces, except for the lowest angle of attack of 9
$^{\circ }$ where they are of similar magnitude (figure 3). This is also the case for the tangential and normal force coefficients (figure 4a,b), showing that, for all three aerodynamic mechanisms, the tangential forces are on average relatively small compared with the pressure forces. To quantify this difference, we determined the tangential-to-normal force coefficient ratio as 
$RC_{Fx-z} = C_{Fx} / C_{Fz}$, for each combination of aerodynamic mechanism and angle of attack (figure 4c). This ratio is smaller than one for all three aerodynamic mechanisms and at all simulated angles of attack, showing that the wing-normal forces dominate over the tangential forces, notably at higher angles of attack.
For the Wagner effect, the ratio remains near zero throughout the complete range of simulated angles of attack, with an average value of 
$RC_{Fx-z_{WE}} =-0.04 \pm 0.03$ (
$n = 10$ angles of attack). For both the translational force and added-mass mechanism, the tangential-to-normal force coefficient ratios vary with angle of attack. As a result, this ratio is of the order of one at low angles of attack, showing that, here, tangential forces contribute significantly to total force production. In contrast, at angles of attack larger than 36
$^{\circ }$ the tangential force coefficients are on average two orders of magnitude smaller than the wing-normal force coefficients (
$RC_{Fx-z_{transl}} =-0.01 \pm 0.01$ and 
$RC_{Fx-z_{\text {AM}}} =0.04 \pm 0.04$, 
$n = 6 \ \text {angles-of-attack})$ (figure 4c).
3.1.2. The force coefficients for insect wings moving with a wing-stroke acceleration
For each modelled aerodynamic mechanism, the variations in normal force coefficients with angle of attack were captured well using a simple sine function (figure 4a), and thus the pressure force coefficients of each sub-model can be expressed as
$$\begin{gather} C_{Fz_{transl}}(\alpha) = C_{Fz\alpha_{transl}} \sin{\alpha}, \end{gather}$$
$$\begin{gather}C_{Fz_{AM}}(\alpha) = C_{Fz\alpha_{AM}} \sin{\alpha}, \end{gather}$$
$$\begin{gather}C_{Fz_{WE}}(\alpha) = C_{Fz\alpha_{WE}} \sin{\alpha} , \end{gather}$$
where 
$C_{Fz\alpha }$ are the angle-of-attack-specific wing-normal force coefficients, equal to 
$C_{Fz\alpha _{transl}} = 3.13 \, (3.10, 3.15)$, 
$C_{Fz\alpha _{AM}} = 0.96 \, (0.95, 0.97)$ and 
$C_{Fz\alpha _{WE}} = -1.02 \, (-1.06, -0.98)$, as the fit result and 95 % confidence interval.
Because the tangential force coefficients of the Wagner effect remain close to zero for the complete range of angles of attack, we here propose that the Wagner effect on tangential forces is negligible (figure 4b). For both the translational model and added-mass model, the tangential force coefficients tend to decrease in magnitude with angle of attack (figure 4b). Therefore, we modelled the variations in these force coefficients with angle of attack using a second-order polynomial and a cosine function, respectively (figure 4b). Note that the latter dependence on the angle of attack is similar to the potential-flow added-mass model, (1.5). The resulting tangential force coefficients of the three sub-models are
$$\begin{gather} C_{Fx_{transl}}(\alpha) = A_{Fx\alpha_{transl}} {\alpha^{2}} + B_{Fx\alpha_{transl}} {\alpha} + C_{Fx\alpha_{transl}}, \end{gather}$$
$$\begin{gather}C_{Fx_{AM}}(\alpha) = C_{Fx\alpha_{AM}} \cos{\alpha}, \end{gather}$$
$$\begin{gather}C_{Fx_{WE}}(\alpha) = 0 . \end{gather}$$
Here, the fit result and 95 % confidence interval of the angle-of-attack-specific tangential force coefficients of the translational model are 
$A_{Fx\alpha _{transl}} = 8.5\times 10^{-5} \, (6.8\times 10^{-5}, 10.2\times 10^{-5})$, 
$B_{Fx\alpha _{transl}} = -1.2\times 10^{-2} \, (-1.4\times 10^{-2}, -1.0\times 10^{-2})$ and 
$C_{Fx\alpha _{transl}} = 0.41 \, (0.37, 0.44)$. The angle-of-attack-specific tangential force coefficient of the added-mass model is 
$C_{Fx\alpha _{AM}} = 0.104 \, (0.101, 0.107)$.
The tangential added-mass force coefficients estimated using our viscous-flow-based scaling (2.2) are therefore of the same order of magnitude as the tangential translational-force coefficients (figure S3b). Alternatively, we can also use these tangential forces to test the added-mass model based on potential-flow theory (4.3), as was used in several previous studies (Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007). Doing so results in tangential pressure force coefficients that are unrealistically high, with an average value of 206 over the entire angle-of-attack range (figure S3a). This reveals that our viscous force scaling captures the tangential added-mass forces much better than the corresponding pressure force scaling (figure S3). This confirms that the tangential added-mass forces are not dominated by pressure forces, but instead by viscous added-mass forces acting on the accelerating wing surface.
In summary, the tangential force component of the Wagner effect can be assumed negligible, and the tangential forces of both the translational and added-mass mechanisms are dominated by viscous forces and therefore most prominent at low angles of attack. At angles of attack larger than 36
$^{\circ }$, these viscous forces comprise only 1 % and 4 % of the total aerodynamic forces, for the translational and added-mass mechanisms, respectively (figure 4c). The wings of a flying insect mostly operate at angles of attack larger than this 36
$^{\circ }$ (figure 1c), and therefore the aerodynamics of accelerating insect wings is dominated by the wing-normal pressure forces. Because of this, we will from here on primarily focus on these pressure forces.
3.1.3. Aerodynamic forces produced by accelerating and decelerating insect wings
We continued to test whether our model based on CFD simulations of fruit fly wings with only zero and positive stroke accelerations is also valid for decelerating wings. We did this by analysing the results of all CFD simulations with fruit fly wings operating at an angle of attack of 36
$^{\circ }$ (figure 4d). This set consists of the above-used simulations with zero and positive stroke accelerations, and additional simulations with negative stroke accelerations. We compared the wing-normal pressure force coefficients resulting from the above-used fit based on simulations at zero and positive stroke accelerations only, with the pressure force coefficients resulting from all simulations, including those with decelerating wings (figure 4d–e). For all three pressure force coefficients, the 95 % confidence intervals of the complete set (all 
$\dot {\omega }$) and reduced set (
$\dot {\omega } \geq 0$ rad s
$^{-2}$) overlap, showing that they do not differ significantly. This indicates that the above-developed and validated acceleration-based model also applies to decelerating wings.
Note that, because the Wagner-effect model includes a square root of stroke acceleration, we needed to adjust it to include negative stroke accelerations as
\begin{equation} F_{z_{WE}} = 1/2 \rho \omega \, \text{sign}(\dot{\omega}) \sqrt{|\dot{\omega}|} S_{WE} C_{Fz_{WE}}(\alpha). \end{equation}3.1.4. Aerodynamic forces produced by accelerating insect wings with various wing planforms
Next, we tested whether the geometric scaling of our aerodynamic model is valid for the range of wing planforms that we used. We did so by fitting our pressure force model ((3.1)–(3.4)) to the results of all CFD simulations of each wing separately, operating at 
$\alpha = 36^{\circ }$ (figure 5). We then tested, for all three aerodynamic mechanisms, whether the pressure force coefficients of the different wings scale predominantly with the corresponding wing geometry parameter. We did so by fitting linear regression models to the data, and testing using 
$t$ statistics whether the regression slope was significantly different from zero (using fitlm in Matlab). For all of the three aerodynamic mechanisms, the pressure force coefficient did not significantly vary with the corresponding wing-planform geometry, assuming a 5 % significance level (translational-force model: 
$p\text {-value}=0.24$; added-mass model: 
$p\text {-value}=0.35$; Wagner-effect model: 
$p\text {-value}=0.16$). This shows that our model captures the effect of wing geometry on aerodynamic pressure force production well.
Estimates of the wing-normal pressure force coefficients vs wing geometry parameters for all tested wings operating at angle of attack 
$\alpha = 36^{\circ }$. For each wing, we show (a) the translational-force coefficient 
$C_{Fz_{transl}}$ vs the spanwise second moment of area 
$S_{yy}$, (b) the added-mass force coefficient 
$C_{Fz_{AM}}$ vs the chord-based second moment of area 
$S_{cy}$ and (c) the Wagner-effect-based force coefficient 
$C_{Fz_{WE}}$ vs the second moment of area of the Wagner-effect model 
$S_{WE}$. Each data point shows the fit and 
$95\,\%$ confidence interval for a single wing type. As shown in the legend on the bottom, the tested wings are the fruit fly wing (blue), the malaria mosquito wing (red), the elliptical wings with variable chord length (green squares), the elliptical wings with variable span (green circles) and the elliptical wings with constant 
$S_{cy}$ (green diamonds). The curves show the linear regression fits through the data as fit result (dashed line) the and 
$95\,\%$ confidence interval of the linear fit (dotted lines).
3.1.5. Aerodynamic forces produced by insect wings within the tested parameter space of stroke rates and stroke accelerations
Finally, we determined how the three different aerodynamic mechanisms contribute relatively to the total aerodynamic pressure force production within the parameter space of stroke rates and stroke accelerations of our simulations (figure 6). This shows that, as expected, the total aerodynamic pressure force 
$F_{total}$ is highest at high stroke rates and high stroke accelerations, which result in high translational and added-mass forces, respectively. The aerodynamic forces become negative for wings that decelerate at low stroke rates (figure 6a). For flapping wings, this occurs directly prior to stroke reversal.
The distribution of aerodynamic pressure forces produced by accelerating and decelerating fruit fly wings throughout the tested parameter space of stroke rates and stroke accelerations. (a) The total weight-normalized aerodynamic force 
${F_z}^{*}_{total}$, as estimated using our aerodynamic model for accelerating wings ((3.1)–(3.8)). (b–d) The aerodynamic forces resulting from the sub-components in our aerodynamic model, including (b) the normalized translational forces 
${F_z}^{*}_{transl}$, (c) the added-mass forces 
${F_z}^{*}_{\text {AM}}$ and (d) the Wagner-effect-based aerodynamic forces 
${F_z}^{*}_{WE}$. In all panels, the parameter space of stroke rates and stroke accelerations is spanned by the axes, and weight-normalized aerodynamic forces are shown in colours, as defined by the colour map on the right.
Within the visualized parameter space, translational pressure forces 
$F_{transl}$ are responsible for the majority of aerodynamic force production, especially at high stroke rates (figure 6b). But added-mass forces 
$F_{AM}$ and Wagner-effect-based forces 
$F_{WE}$ cannot be ignored (figure 6c,d). The added-mass forces are dominant at low stroke rates, and the Wagner effect is highest when both stroke rates and stroke accelerations are large. Note that the added-mass mechanism and Wagner effect tend to produce pressure forces of opposite signs: for an accelerating wing, the added-mass forces are positive but the Wagner-effect-based forces are negative. For a decelerating wing, this is reversed.
3.2. Aerodynamics of acceleration-based force production in insect wings
We studied the aerodynamics of acceleration-based wing-normal force production by separately extracting the flow dynamics that results in added-mass forces 
$F_{AM}$ and Wagner-effect-based forces 
$F_{WE}$, using the procedure described in § 2.7.
3.2.1. Aerodynamics of the added-mass-based force production
We estimated the pressure distribution that would be caused by the added-mass mechanism using (2.12). By doing so, we assumed that, at the combination of lowest simulated stroke rate and highest stroke acceleration, added-mass forces are higher than Wagner-effect-based forces. We tested this using our here-developed aerodynamic model, as shown in figure 6. At the combination of lowest simulated stroke rate and highest simulated stroke acceleration, the added-mass forces are 6.6 times higher than the Wagner-effect-based forces. This shows that, although the Wagner effect is not negligible, the majority of pressure forces determined using (2.12) result in added-mass force production.
The added-mass-based pressure distribution 
$p_{AM}$ around an accelerating fruit fly wing (figure 7) displays a low-pressure region above the wing and a high-pressure region on the bottom side of the wing. The peak in low pressure above the wing appears to be smaller in magnitude than the high-pressure peak below the wing. Thus, the added-mass forces are primarily caused by increased pressure below the wing, and less so by reduced pressures above the wing.
The effect of stroke acceleration on the pressure distribution resulting from the added-mass mechanism (
$p_{AM}$) around a fruit fly wing, as determined using (2.12). All wings were moving at a stroke rate of 
$\omega = 500$ rad s
$^{-1}$ and angle of attack of 
$\alpha = 36^{\circ }$. The first to last columns show the stroke-acceleration-induced pressure distributions at stroke accelerations of 
$1\times 10^{6}$, 
$2\times 10^{6}$, 
$3\times 10^{6}$ and 
$4\times 10^{6}$ rad s
$^{-2}$, respectively. (a) Schematic representation of a fruit fly wing, including the extraction plane used for panels b–e. (b–m) Air-pressure distribution 
$p_{AM}$ around the wing (b–e), on the wing top surface ( f–i) and on the wing bottom surface (j–m).
The pressure region on top of the wing shows a small asymmetry along the chord, as the lower pressure region is located more towards the trailing edge. This might be due to an interaction with the LEV on the top surface of the wing. The pressure region below the wing is symmetrical with respect to the chord and has its maximum at around two third of the wing span. The pressures on both sides of the wing increase approximately linearly with stroke accelerations (figure 7), which agrees with the observed linear dependence of 
$F_{AM}$ on the stroke acceleration (2.2).
The change in added-mass-based pressure with angle of attack 
$\alpha$ shows a similar pattern as the change with stroke acceleration 
$\dot {\omega }$ (figure 8). With increasing angle of attack, both the high- and low-pressure regions increase in magnitude. But these variations are not linear, because the pressure increase with angle of attack is high at low angles, and then levels off at angles near 90 degrees (figure 8a–d). Interestingly, also with increasing angle of attack, the peak in low pressure above the wing moves further away from the root (figure 8i–l).
The effect of the angle of attack on the pressure distribution resulting from the added-mass mechanism (
$p_{AM}$, (2.12)) around a fruit fly wing moving at a stroke rate of 
$\omega = 500$ rad s
$^{-1}$ and a stroke acceleration of 
$\dot {\omega } = 4\times 10^{6}$ rad s
$^{-2}$. The first to last columns show the stroke-acceleration-induced pressure distributions at an angle of attack of 27
$^{\circ }$, 45
$^{\circ }$, 63
$^{\circ }$ and 81
$^{\circ }$, respectively. (a–l) Air-pressure distribution 
$p_{AM}$ around the wing (a–d), on the wing top surface (e–h) and on the wing bottom surface (i–l). The pressure field planes in (a–d) are defined in figure 7(a).
We used elliptical wings with variable chord length to study how the thickness of the fluid layer that is accelerated with the wing (
$\delta _{AM}$) varies along the span and with wing geometry (figure 9). We estimated the thickness of this fluid layer from the pressure distribution 
$p_{AM}$ on the wing surface using the method described in § 2.7.1 (figure 9).
The effect of wing-planform geometry on the thickness of the accelerated fluid layer 
$\delta _{AM}$ (2.13) for the set of elliptical wings with varying maximum chord length. (a,b) Pressure distribution on the top surface (a) and bottom surface (b) of an elliptical wing with maximum chord length 
$c = 1$ mm moving at an angle of attack of 
$\alpha = 36^{\circ }$, stroke rate 
$\omega = 500$ rad s
$^{-1}$ and stroke acceleration 
$\dot {\omega } = 4.5\times 10^{6}$ rad s
$^{-2}$. (c,d) Pressure distribution on the top surface (c) and bottom surface (d) of the same wing as in (a,b), but now moving at zero stroke acceleration. (e, f) The stroke-acceleration-induced pressure field 
$p_{AM}$ on the top (e) and bottom wing surfaces ( f), defined as the difference between the pressure fields in panels (a,b) and (c,d). (g) Schematic representation of the thickness of the fluid layer accelerated by the elliptical wing. (h–l) Fluid-layer thickness distribution across the wing span for the set of five elliptical wings with varying maximum chord lengths operating at an angle of attack of 
$\alpha = 36^{\circ }$, a stroke rate of 
$\omega = 500$ rad s
$^{-1}$ and a stroke acceleration of 
$\dot {\omega } = 4.5\times 10^{6}$ rad s
$^{-2}$.
This shows that the thickness of the fluid layer increases with wing chord length, and it decreases along the span (figure 9h–l): the fluid-layer thickness is highest near the wing root, and lowest at the wing tip. The conventional added-mass model assumes that the fluid-layer thickness is proportional to the local wing chord length throughout the span (1.4). Our CFD computations show that this assumption is not correct for a wing with a rotational stroke acceleration.
3.2.2. Aerodynamics of Wagner-effect-based force production
We estimated the pressure distribution that would be caused by the Wagner-effect mechanism using (2.14). Here, we test this method by comparing the Wagner-effect forces from our model (3.4) with the model forces calculated using the same approach as in (2.14). Both calculations result in exactly the same force estimation, highlighting the validity of this analysis approach.
The pressure distribution due to the Wagner effect 
$p_{WE}$ around a fruit fly wing moving at variable stroke accelerations is shown in figure 10. At all tested stroke accelerations, the Wagner-effect-based pressure on the bottom of the wing remains near zero, and thus the air pressure below the wing is almost unaffected by the Wagner effect (figure 10e–h). In contrast, the pressure distribution on top of the wing is strongly affected, and this effect increases with the stroke-acceleration level (figure 10a–d).
The effect of stroke acceleration on the pressure distribution resulting from the Wagner effect (
$p_{WE}$), for the fruit fly wing operating at an angle of attack 
$\alpha = 36^{\circ }$. The Wagner-effect air pressures were determined using (2.14) from simulations at stroke rates 
$\omega _1 = 1250$ and 
$\omega _2 = 750$ rad s
$^{-1}$, and variable stroke accelerations. Hereby, the first to last columns show the results for wings moving at a stroke acceleration of 
$1\times 10^{6}$, 
$2\times 10^{6}$, 
$3\times 10^{6}$ and 
$4\times 10^{6}$ rad s
$^{-2}$, respectively. (a–l) The Wagner-effect pressure distribution 
$p_{WE}$ on the top surface of a wing (a–d), the bottom surface of a wing (e–h) and the pressure difference across the wing surface (i–l).
The Wagner-effect-based pressure distribution on top of the wing shows a complex pattern: the pressure is elevated within a large band running diagonally across the wing surface, and near the wing tip a small region with negative pressure is formed. The increase in pressure on top of the wing causes a reduction of the pressure difference across the wing surface (
$\Delta p_{WE} = p_{{WE}_{bottom}} - p_{{WE}_{top}}$) (figure 10i–l), explaining why the Wagner-effect-based forces are negative for an accelerating wing.
When varying the angle of attack, we find a similar behaviour (figure 11a–d). The Wagner-effect-based pressure difference across the wing 
$\Delta p_{WE}$ is mostly negative, and becomes more negative with increasing angle of attack. At the same time, the area of highest negative pressure moves from the tip region more towards the root of the wing.
The Wagner-effect-based pressure difference across the wing surface (a–d), and the spanwise vorticity around the fruit fly wing (e–l), at various angles of attack. The first to last columns show data for fruit fly wings operating at 
$\alpha = 27^{\circ }$, 
$\alpha = 45^{\circ }$, 
$\alpha = 63^{\circ }$ and 
$\alpha = 81^{\circ }$, respectively. (a–d) The Wagner-effect pressure difference across the wing (
$\Delta p_{WE}$) was determined using (2.14) from simulations at stroke rates 
$\omega _1 = 1250$ rad s
$^{-1}$ and 
$\omega _2 = 750$ rad s
$^{-1}$, and stroke acceleration 
$\dot {\omega } = 4.5\times 10^{6}$ rad s
$^{-2}$. (e–l) The spanwise vorticity around steady moving and accelerating fruit fly wings (top and bottom rows, respectively). (e–h) The steady moving wing has a stroke rate 
$\omega = 1250$ rad s
$^{-1}$ and zero stroke acceleration. (i–l) The accelerating fruit fly wing is moving at the same stroke rate but also accelerates at a stroke acceleration of 
$\dot {\omega } = 4.5\times 10^{6}$ rad s
$^{-2}$. The vorticity field planes in (e–l) are defined in figure 7(a).
In summary, these results show that the Wagner effect primarily causes an increase in air pressure on top of the wing (
$p_{{WE}_{top}}$), resulting in a reduction in the pressure difference across the wing surface (
$\Delta p_{WE}$). This effect increases with both stroke acceleration and angle of attack. The reduction in 
$\Delta p_{WE}$ explains why the Wagner-effect-based forces reduce the total aerodynamic forces on the accelerating wing (§ 3.1).
To visualize the mechanism behind the Wagner-effect-based increase in pressure on top of the wing (figure 10a–d), we compared the spanwise vorticity around wings that move without and with a stroke acceleration (figures 11e–h and 11i–l, respectively). Here, all wings moved at a stroke rate of 
$\omega = 1250$ rad s
$^{-1}$, the highest modelled stroke acceleration of 
$\dot {\omega } = 4.5\times 10^{6}$ rad s
$^{-2}$ and variable angle of attack. The results show that the wings moving with a stroke acceleration have a smaller LEV compared with the non-accelerating wing, for all tested angles of attack. This reveals that the reduced aerodynamic forces in the accelerating wing are a consequence of an acceleration-induced reduction in the growth of the LEV.
Next, we tested whether the same mechanism of delayed LEV development also explained the Wagner-effect-based increase in aerodynamic forces in decelerating wings. For this, we compared the spanwise vorticity fields around accelerating and decelerating wings (figure 12). All these fruit fly wings were operating at the same stroke rate (
$\omega = 750$ rad s
$^{-1}$) and angle of attack 
$\alpha = 36^{\circ }$, and thus the differences in vorticity fields are due to the difference in stroke acceleration. This comparison reveals that the LEVs on accelerating wings are smaller than for the steady moving wing, whereas the LEVs on decelerating wings are larger than for the steady moving wing. The effect of stroke acceleration is most striking when directly comparing the LEV strength between the wings accelerating and decelerating at the same values (right and left columns in figure 12, respectively). Also, the change in LEV strength is larger between the steady moving wing (figure 12a,d) and the wings acceleration and decelerating at 
$2\times 4$ rad s
$^{-2}$ (figure 12b,e) than between the 
$2\times 4$ and 
$4\times 10^{4}$ rad s
$^{-2}$ cases (figure 12b,e vs 12c, f). This can be explained by the fact that at the here simulated stroke rate (
$\omega = 750$ rad s
$^{-1}$), the change in Wagner effect is largest at low stroke accelerations, and levels off at higher values (figure 6d).
The spanwise vorticity around decelerating and accelerating fruit fly wings shows that the Wagner effect causes an increase in LEV strength in decelerating wings, equivalent to the decrease in LEV strength of accelerating wings. (a–c) Spanwise vorticity around fruit fly wings decelerating at stroke accelerations 
$\dot {\omega } = 0$, 
$-2\times 10^{6}$ and 
$-4\times 10^{6}$ rad s
$^{-2}$, respectively. (d–f ) Spanwise vorticity around wings accelerating at the equivalent stroke accelerations 
$\dot {\omega } = 0$, 
$2\times 10^{6}$ and 
$4\times 10^{6}$ rad s
$^{-2}$, respectively. All wings were moving at the same stroke rate and angle of attack (
$\omega = 750$ rad s
$^{-1}$ and 
$\alpha = 36^{\circ }$, respectively), and thus only the stroke accelerations differed. The vorticity field planes are defined in figure 7(a).
The results of figures 11 and 12 combined confirm that the Wagner effect causes not only a delay in LEV growth on accelerating wings, but that it also causes a delay in the decrease of LEV strength for decelerating wings. This supports the validity of our proposed physical modelling of the Wagner-effect mechanism, for both accelerating and decelerating wings, and at the large range of tested angles of attack.
3.3. Application of the new model: acceleration-based forces throughout the wing-beat kinematics of a hovering fruit fly and malaria mosquito
Based on our newly developed aerodynamic force model for accelerating wings ((3.1)–(3.14)), we formulated an equivalent quasi-steady aerodynamic model for predicting the translational forces, added-mass forces and Wagner-effect-based forces produced by flapping insect wings (see materials and methods § 2.8). We then applied this quasi-steady model to the kinematics of a hovering fruit fly (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014) and a hovering malaria mosquito (Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a). By doing so, we tested how the different acceleration-based force components affect lift, drag and pressure forces during the forward and backward wing strokes of these flying insects (figure 13a–d).
Application of the quasi-steady aerodynamic force model for accelerating wings to the wing motion of (a,b) a hovering fruit fly (Muijres et al. Reference Muijres, Elzinga, Melis and Dickinson2014), and (c,d) a hovering malaria mosquito (Muijres et al. Reference Muijres, Chang, Veen, Spitzen, Biemans, Koehl and Dudley2017a). The total weight-normalized pressure forces are in grey; the separate force components translational force, added-mass force and Wagner-effect force are in blue, green and red, respectively (see legend at top). (e,g) Weight-normalized pressure forces on the fly and mosquito wing, respectively. ( f,h) Average pressure forces on the fly and mosquito wing during the forward and backward strokes. (i,k) Weight-normalized lift forces perpendicular to the wing velocity vector for the fly and mosquito, respectively. (j,l) Average lift force during the forward and backward strokes for the fly and mosquito, respectively. (m,o) Weight-normalized drag forces, parallel to the wing velocity vector, for the fly and mosquito, respectively. (n,p) Average drag forces during the forward and backward strokes for the fly and mosquito, respectively. The vertical dashed lines indicate stroke reversal.
For the hovering fruit fly, the translational forces dominate aerodynamic force production, as the translational-force mechanism produces 96 % and 98 % of the wing-beat-average lift and drag forces, respectively (figures 13j and 13n, respectively). For the hovering mosquito, the translational forces are also dominant, but this is less pronounced as here the translational forces produce 92 % and 77 % of the wing-beat-average lift and drag, respectively (figures 13l and 13p, respectively).
For the fruit fly, added-mass forces explain 6 % of the total lift and less than 1 % of the total drag, and the Wagner effect explains –2 % of lift and 2 % of drag. For the mosquito, added-mass forces explain 15 % of the wing-beat-average lift and 31 % of the wing-beat-average drag. The Wagner effect explains –7 % and –8 % of the wing-beat-average lift and drag, respectively.
This shows that acceleration-based forces are primarily relevant for the mosquito with its high-frequency and low-amplitude flapping wing beat. For these high-frequency flappers, the added-mass mechanism tends to contribute positively to both lift and drag production, whereas the Wagner effect has on average a negative effect on both.
The temporal dynamics of the acceleration-based forces shows that the instantaneous acceleration-based forces are quite significant for both the fly and mosquito, as these forces are sometimes even higher than the translational forces (figure 13). For both the fly and mosquito, the translational forces are close to zero at stroke reversal, and maximum around mid wing stroke where the stroke rates are highest. In contrast, both the added-mass forces and Wagner-effect forces have two peaks in force production during a single wing stroke (forward stroke and backward stroke): one force peak occurs during the acceleration phase at the start of the wing stroke, and a second peak occurs during the deceleration phase at the end of the wing stroke (figure 13). For the added-mass mechanism, the acceleration force peak is positive, and the deceleration force peak is negative. For the Wagner effect, this dynamics is opposite as Wagner-effect-based forces scale negatively with stroke acceleration (2.10). As a result, the instantaneous added-mass forces and Wagner-effect-based forces tend to partly cancel each other throughout the wing beat.
Despite the fact that the instantaneous added-mass forces and Wagner-effect forces are opposite in sign throughout the wing beat, they do not fully cancel each other and thus they still produce significant instantaneous and wing-beat-average forces, especially for high-frequency flapping mosquitoes.
4. Discussion
4.1. An aerodynamic model for accelerating insect wings
In this study, we showed that the aerodynamic forces resulting from a rotational wing-stroke motion produced by an insect wing (
$F_{total}$) consist of three terms: the translational forces (
$F_{transl}$), the added-mass forces (
$F_{AM}$) and the Wagner-effect-based forces (
$F_{WE}$). The last two force components depend on wing-stroke accelerations, and thus we call them combined the acceleration-based forces (
$F_{accel}$). Here, we developed novel aerodynamic force models for the translational-force mechanism, the added-mass mechanism and the Wagner effect. We did this by linking the stroke rate 
$\omega$, stroke acceleration 
$\dot {\omega }$, the angle of attack 
$\alpha$ and the wing-planform geometry to the generation of these unsteady aerodynamic forces on accelerating and decelerating insect wings.
Our model differs from previous aerodynamic models for insect wings in several aspects (table 1). The three most important ones are:
(i) Our added-mass model is among the most complete models as it includes both tangential viscous forces and a variable virtual mass (table 1). Both aspects have been modelled separately (Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007; Lee et al. Reference Lee, Lua, Lim and Yeo2016; Cai et al. Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021), but no previous study has combined both tangential viscous forces and a variable virtual mass in a single added-mass model for flapping insect wings (table 1).
(ii) We combined our added-mass model with a novel generalized Wagner-effect model, which models the Wagner effect as a stroke-acceleration-based reduction in translational-force production. This model is defined in the time domain, instead of the spatial domain as in a previous study (Wang et al. Reference Wang, Goosen and Van Keulen2016), and therefore can be implemented in a quasi-steady approach.
(iii) We based our model on 721 CFD simulations of 14 different wing geometries accelerating with variable stroke kinematics. This allowed us to explicitly test the dependence of the unsteady aerodynamics on wing-planform shape. In contrast, most other models are based on a single wing, or only a few different wing types (table 1).
Next to these novel aspects of our study, our systematic numerical modelling approach also has its limitations. The most important ones are:
(i) To isolate the stroke-acceleration-based forces from other aerodynamic mechanisms such as wake capture and rotational lift, we kept the angle of attack
$\alpha$ constant throughout each simulation, and did not model stroke reversal. Therefore, we did not test how added-mass forces and the Wagner effect relate to wing-pitch motions or wake capture. Note that we did test, however, how the stroke-acceleration-based forces depend on pitch angle (i.e. the angle of attack).(ii) Our novel generalized Wagner-effect model is a linear approximation of a complex aerodynamic mechanism that is inherently nonlinear (Wagner Reference Wagner1925; Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015). Therefore, this simplified model is intended to produce a first-order estimate of Wagner-effect forces in accelerating insect wings; more sophisticated analyses are needed to study the detailed Wagner effects on flapping insect wings.
(iii) We designed our numerical study such that all CFD simulations took place in exactly 1 ms. This allowed us to systematically vary stroke rates and stroke accelerations without varying the temporal dynamics of the start-up effect. An alternative approach would be to keep distance travelled constant between simulations. A future study with such alternative approach would be useful to test the robustness of our aerodynamic model.
(iv) In this study, we considered rigid wings. However, the wings of insects are flexible and can therefore deform under aerodynamic and inertial loads (Combes & Daniel Reference Combes and Daniel2003). These deformations might influence the generation of forces. Based on the pressure distribution of the added-mass forces (figure 7) the centre of pressure can be approximated at half the chord length and two thirds of the span. The location of the centre of pressure lies close to the rotation axes, and therefore we expect that the wing deformations due to stroke acceleration will not be as large as for instance the wing deformation due to pitch-based motions (Dickinson et al. Reference Dickinson, Lehmann and Sane1999; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Bomphrey et al. Reference Bomphrey, Nakata, Phillips and Walker2017; van Veen et al. Reference van Veen, van Leeuwen and Muijres2019). Moreover, although wing deformations might affect the forces produced by accelerating wings, the underlying aerodynamic mechanisms would of course not change.
4.2. The Wagner-effect forces for accelerating insect wings
As part of our general aerodynamic model for accelerating insect wings, we developed a sub-model for the Wagner-effect-based forces that are produced by these accelerating wings. Most previous efforts to model the Wagner effect on accelerating wings have been based on the (classic) case of a wing impulsively started from rest, and related to the distance travelled after the start of the wing movement (Wagner Reference Wagner1925; Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015; Wang et al. Reference Wang, Goosen and Van Keulen2016). These types of models do not lend themselves well to application to flapping insect wings, because such flapping wings continuously accelerate and decelerate (Wang et al. Reference Wang, Goosen and Van Keulen2016). Furthermore, the concept of an impulsive start is physically unrealistic, as they imply infinite acceleration forces at the start. Because of these limitations, Wagner effects are commonly ignored for insect flight, which has been mentioned as an important caveat of previous models (Lee et al. Reference Lee, Lua, Lim and Yeo2016).
We solved this issue by modelling the aerodynamic forces resulting from the Wagner effect as an acceleration-induced modification of the translational forces, defined in the temporal domain instead of the spatial domain (chord lengths travelled). In our study, we showed that this time-dependent Wagner-effect model captures the effect of wing-stroke accelerations on translational-force production well, for both accelerating and decelerating wings (figure 4).
The chord-wise tangential force coefficients were found to be very small for the complete range of tested angles of attack (figure 4a–c), and therefore we assumed the tangential Wagner-effect forces to be negligible. The wing-normal pressure force coefficients of the Wagner-effect model were of similar magnitude as the pressure coefficients of the other modelled aerodynamic mechanisms, but with negative values (figure 4a). This confirms that, as predicted, the Wagner effect results in an acceleration-based reduction in the translation-based forces on an accelerating wing. Our proposed Wagner-effect-based aerodynamic model is therefore defined as
$$\begin{gather} F_{z_{WE}} = 1/2 \rho \omega \, \text{sign}(\dot{\omega}) \sqrt{|\dot{\omega}|} S_{WE} C_{Fz\alpha_{WE}} \sin{\alpha}, \end{gather}$$
$$\begin{gather}F_{x_{WE}} = 0 , \end{gather}$$
where the angle-of-attack-specific wing-normal force coefficient for the Wagner effect is 
$C_{Fz\alpha _{WE}} = -1.02 \, (-1.06, -0.98)$.
The model predicts that, due to inertia of the fluid, the translational forces on an accelerating wing are lower than assumed based on the stroke rate alone (figure 6). In contrast, translational forces on a decelerating wing are higher than for the steady case, as a result of an inversion of the Wagner effect. The fact that the same model (without adapting the model coefficients) applies for these two conditions further justifies the new model (figure 4d). The contribution of this effect is strongest at large stroke accelerations in combination with large stroke rates (figure 6d).
Comparing the structure of our simplified Wagner-effect model with previous studies (Wagner Reference Wagner1925; Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015; Wang et al. Reference Wang, Goosen and Van Keulen2016) shows a striking similarity. This is most apparent when comparing our model with the discussion in Pullin & Wang (Reference Pullin and Wang2004). Here, the Wagner solution is shown to consist of three terms, representing the translational lift, the added mass and a third term that represents the non-equilibrium effect of the acceleration on the circulatory lift, i.e. the Wagner effect. Our data-driven model contains the same three terms, whereby the Wagner effect is modelled as resulting from a similar non-equilibrium effect. This confirms that our linear Wagner-effect model captures a similar physical mechanism as described in the original studies on the topic (Wagner Reference Wagner1925; Pullin & Wang Reference Pullin and Wang2004).
Our analysis of the airflow dynamics around acceleration and decelerating wings highlights the aerodynamic mechanism that causes Wagner-effect-based forces (figures 11 and 12). The LEV is important for aerodynamic force production by flapping insect wings (van den Berg & Ellington Reference van den Berg and Ellington1997; Sane & Dickinson Reference Sane and Dickinson2002; Sane Reference Sane2003). For accelerating wings, this LEV is observed to be smaller than for wings moving with a constant stroke rate. In contrast, decelerating wings display a larger LEV than for the steady case. These dynamics can directly be attributed to the Wagner effect as, due to inertia of the fluid, the development of the LEV is delayed. As a result, the LEV on an accelerating wing grows less fast than would be expected from the increase in stroke rate alone, and LEVs on decelerating wings reduce in size more slowly than based on the decrease in stroke rate alone. Our Wagner-effect model captures this dynamics using the stroke-acceleration-dependent correction on translational pressure force production.
4.3. The added-mass forces produced by accelerating insect wings
Here, we developed an aerodynamic model for added-mass force production by an insect wing with a stroke acceleration. The resulting model includes both a pressure force component normal to the wing surface, and a chord-wise tangential force component.
We showed that these tangential added-mass forces are best captured by a viscous added-mass force model dependent on the wing surface area instead of a tangential pressure force model dependent on the frontal area of the wing (figure S3). This suggests that the tangential added-mass forces on an accelerating insect wing are governed by viscous forces, and not by pressure forces as often implicitly assumed (Sedov Reference Sedov1965; Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007). Our findings confirm the results of a previous study that suggests that viscosity-based (centripetal) added-mass forces are not negligible in insect flight (Zhang et al. Reference Zhang, Hedrick and Mittal2015).
Our newly developed added-mass force model for accelerating insect wings is defined as
$$\begin{gather} F_{z_{AM}} = \rho \dot{\omega} {S}_{cy} C_{Fz\alpha_{AM}} \sin{\alpha} , \end{gather}$$
$$\begin{gather}F_{x_{AM}} = \rho \dot{\omega} {S}_{cy} C_{Fx\alpha_{AM}} \cos{\alpha} , \end{gather}$$
where 
$C_{Fz\alpha _{AM}} = 0.96 \, (0.95, 0.97)$ and 
$C_{Fx\alpha _{AM}} = 0.104 \, (0.101, 0.107)$. The angle-of- attack-specific tangential viscous force coefficient is thus an order of magnitude smaller than the angle-of-attack-specific wing-normal pressure force coefficient. This shows that the pressure forces dictate added-mass force production in the wing-normal direction, especially at the high angles of attack at which flapping insect wings tend to operate (figure 1c). This confirms that the potential-flow-based added-mass model as often used for insect flight is mostly valid (Sedov Reference Sedov1965; Sane & Dickinson Reference Sane and Dickinson2002; Lee et al. Reference Lee, Lua, Lim and Yeo2016), although viscous added-mass forces cannot be ignored. Moreover, viscous forces dominate the tangential-force component at all angles of attack, and they dominate total aerodynamic forces at low angles of attack.
Our model for the added-mass-based pressure forces (4.3) resembles several models based on two-dimensional potential flow (table 1). Most of these previous models assumed a cylindrical added volume with a diameter equal to the local wing chord length (Sedov Reference Sedov1965; Sane & Dickinson Reference Sane and Dickinson2002; Andersen et al. Reference Andersen, Pesavento and Wang2005; Berman & Wang Reference Berman and Wang2007; Nakata et al. Reference Nakata, Liu and Bomphrey2015; Wang et al. Reference Wang, Goosen and Van Keulen2016). This results in an added-mass coefficient of 
${\rm \pi} /4 \approx 0.79$ (1.4). Two studies tested this assumption of a cylindrical added volume by adding a variable virtual mass to their insect flight model (Lee et al. Reference Lee, Lua, Lim and Yeo2016; Cai et al. Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021). Using this approach, both studies predicted that the added volume in insect flight is approximately 20 % lower than the cylindrical added volume. In contrast, here we estimated the added volume as 22 % larger than the 
${\rm \pi} /4$ of a cylindrical added volume.
This discrepancy between our added-volume estimate and that of Lee et al. (Reference Lee, Lua, Lim and Yeo2016) and Cai et al. (Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021) can be traced back to the absence of a Wagner-effect model in these previous studies. Throughout a wing beat, the added mass and Wagner effect tend to produce forces of opposite sign (figure 13), and thus the Wagner effect often cancels part of the added-mass forces. When these Wagner effects are not modelled, the acceleration-based forces are assumed to result from the added-mass mechanism only. Because the Wagner effect counteracts the effect of added mass, this might cause an under-prediction of the added-mass-based force coefficient.
We tested this hypothesis and how large this effect might be by removing the Wagner-effect component from our wing-normal aerodynamic force model for accelerating wings (3.1). We then again fitted this reduced model to all CFD simulations with the fruit fly wing at zero and positive stroke accelerations. This is equivalent to the analysis shown in figures 3 and 4(a,b), but then without the Wagner-effect model included. The resulting angle-of-attack-specific pressure force coefficients of the translational and added-mass mechanisms are 
$C_{Fz\alpha _{transl}} = 2.68 \, (2.65, 2.71)$ and 
$C_{Fz\alpha _{AM}} = 0.73 \, (0.73, 0.74)$, respectively. These values are 14 % and 24 % lower than the translational and added-mass coefficients when the Wagner-effect component is included, respectively. Most importantly, the added-mass coefficients that we estimated without modelling the Wagner effect are in reasonable agreement with the result of Lee et al. (Reference Lee, Lua, Lim and Yeo2016) (0.73 vs 0.65). This confirms that, by ignoring the Wagner effect, the contribution of added mass to aerodynamic force production for flying insects is under-estimated.
Our modelling approach thus showed that the effective added volume of air around the accelerating wing is 22 % larger than predicted using potential-flow theory (Sedov Reference Sedov1965; Sane & Dickinson Reference Sane and Dickinson2002; Lee et al. Reference Lee, Lua, Lim and Yeo2016), and 42 % larger than estimated previously (Lee et al. Reference Lee, Lua, Lim and Yeo2016; Cai et al. Reference Cai, Kolomenskiy, Nakata, Liu and Cai2021). This is confirmed by our pressure field analysis, which indicates that the shape of the accelerated fluid was far from cylindrical (figure 9). Most importantly, our simulations show that the thickness of the accelerated fluid layer varies along the span, irrespective of the spanwise variations in chord length. In fact, the fluid layer was thickest near the wing root and thinnest near the wing tip. As a result, the simplified cylindrical model underestimates the volume of added fluid near the root, and overestimates this volume near the wing tip, resulting in a 22 % underestimation of the added volume for the complete wing.
The relatively large fluid-layer thickness near the wing root does not result in large added-mass forces near the root, because translational accelerations resulting from the revolving stroke motion reduce to zero at the root. This is most apparent from the pressure distributions resulting from the added-mass mechanism (figures 7–9). Note that, here, we estimated the spanwise fluid-layer thickness distribution, but not the chord-wise distribution. In our visualization, we assumed the chord-wise distribution to be elliptical (figure 9). Additional research is needed to also estimate this chord-wise distribution.
Finally, the added-mass forces on flapping insect wings have been shown in the literature to vary with Reynolds number (Lee et al. Reference Lee, Lua, Lim and Yeo2016). To model this effect, Lee et al. (Reference Lee, Lua, Lim and Yeo2016) used a Reynolds-number-dependent correction factor. Additional Reynolds-number-specific simulations are needed to test how the Reynolds number would affect our added-mass model, and potentially also the Wagner effect.
4.4. Translational forces produced by revolving insect wings with a stroke acceleration
Although it was not the primary focus of our research, our modelling approach also allowed us to develop an aerodynamic model of the translational forces on revolving insect wings. The resulting translational-force model is very similar to those previously developed (table 1), with the most important difference that, in the present discussion, we chose to decompose the total aerodynamic force into a wing-normal (pressure) force component and a chord-wise (viscous) force. Most previous models use a decomposition in lift and drag instead. Our resulting translational-force model is
$$\begin{gather} F_{z_{transl}} = 1/2 \rho {\omega}^{2} {S}_{yy} C_{Fz\alpha_{transl}} \sin{\alpha} , \end{gather}$$
$$\begin{gather}F_{x_{transl}} = 1/2 \rho {\omega}^{2} {S}_{yy} C_{Fx_{transl}}(\alpha) , \end{gather}$$
where the angle-of-attack-specific pressure force coefficient is 
$C_{Fz\alpha _{transl}} = 3.13 \, (3.10, 3.15)$ and the viscous force coefficient is
\begin{equation} C_{Fx_{transl}}(\alpha) = A_{Fx\alpha_{transl}} {\alpha^{2}} + B_{Fx\alpha_{transl}} {\alpha} + C_{Fx\alpha_{transl}} . \end{equation}
The coefficients of the viscous force coefficients are 
$A_{Fx\alpha _{transl}} = 1.0\times 10^{-4} \, (0.7\times 10^{-4}, 1.3\times 10^{-4})$, 
$B_{Fx\alpha _{transl}} = -1.4\times 10^{-2} \, (-1.6\times 10^{-2}, -1.1\times 10^{-2})$ and 
$C_{Fx\alpha _{transl}} = 0.43 \, (0.37, 0.48)$.
We compared the force coefficients of our model with those of the fruit fly model of Sane & Dickinson (Reference Sane and Dickinson2002). In that study, the maximum pressure force coefficient equals 3.5, and the maximum viscous force coefficient is 0.39. These are 11 % higher and 5 % lower than the maximum pressure and viscous force coefficients of our model, respectively. This similarity between the two models suggests that the translational-force model developed using wings rotating at a constant stroke rate (Sane & Dickinson Reference Sane and Dickinson2002) is largely also valid for accelerating wings.
4.5. Importance for insect flight
Wing-acceleration-based aerodynamic mechanisms are often not taken into account when modelling the aerodynamics of insect flight. To test whether this is a valid simplification, we applied a quasi-steady version of our aerodynamic wing-stroke model to the wing-beat kinematics of a hovering fruit fly and malaria mosquito (figure 13).
When focusing on the results for the hovering fruit fly (figure 13, left column) we see that, on average, translational forces produce the majority of lift and drag. Thus, when one is only interested in wing-beat-average force production, then it is reasonable to neglect acceleration-based aerodynamic forces for these insects. In contrast, the instantaneous acceleration-based aerodynamic forces throughout the wing beat can certainly not be ignored, as they are often of the same order of magnitude as the translational forces. Hence, they are also relevant for the power requirement of flight (Ellington Reference Ellington1984a).
This effect is even more explicit for the hovering mosquito that has a high-frequency, low-amplitude wing beat (figure 13, right column). For the mosquito, the added-mass forces are the largest force component during both the start and end of the wing beat, when the wing rapidly accelerates and decelerates, respectively. As a consequence, for these high-frequency flapping flyers also wing-stroke-average forces are, for the large part, governed by acceleration-based aerodynamic forces. This confirms the conclusions of a recent study that showed that added-mass forces are particularly important for high-frequency flyers such as mosquitoes (Liu et al. Reference Liu, Du and Sun2020). This study did not model the Wagner effect, although they do mention that the Wagner-effect-based delay in LEV development is particularly relevant for these high-frequency flyers.
For the fruit fly and malaria mosquito, the wing-beat-average lift forces estimated using our model equal to 64 % and 48 % of their average body weights, respectively. These percentages are similar to the relative lift forces when ignoring rotational lift (Sane & Dickinson Reference Sane and Dickinson2002; Dickinson & Muijres Reference Dickinson and Muijres2016; Bomphrey et al. Reference Bomphrey, Nakata, Phillips and Walker2017; Liu et al. Reference Liu, Du and Sun2020). The fact that this percentage is smaller for the mosquito supports the suggestion that rotational lift and drag forces near stroke reversal are also particularly important in these high-frequency flapping flyers (Bomphrey et al. Reference Bomphrey, Nakata, Phillips and Walker2017; Liu et al. Reference Liu, Du and Sun2020). A further study that combines our models of the Wagner effect and added mass with wing-pitch rate and pitch acceleration models would allow us to determine how these mechanisms together govern mosquito flight aerodynamics.
When comparing all force components produced during the hovering wing beat of a fruit fly and malaria mosquito (figure 13), we find that the largest stroke-acceleration-related forces are the added-mass forces on flying mosquitoes. These added-mass forces augment wing-beat-average lift production by 15 %, but this comes at the cost of a steep 30 % increase in drag (figures 13l and 13p, respectively). The instantaneous added-mass forces are even higher, as they dominate aerodynamic force production during the start and end of each wing stroke (figure 13k,o). During the start of the wing stroke, the high wing acceleration results in large added-mass-based drag forces, which thus increase the energetic cost of flight. In contrast, during the deceleration phase of the wing the added-mass forces are negative, resulting in a negative total drag force.
To counteract this negative added-mass-based drag (i.e. thrust) during the end of the wing stroke, the animal should use an active wing braking mechanism. Alternatively, the animal can use an elastic system that stores elastic energy during the deceleration phase prior to stroke reversal, and releases this elastic energy during the consecutive acceleration phase after stroke reversal. This elastic storage could occur in the flexible thorax that also operates as the flight motor system (Dickinson & Lighton Reference Dickinson and Lighton1995) and in the deformable wings (Reid et al. Reference Reid, Schwab, Maxcer, Peterson, Johnson and Jankauski2019).
The Wagner effect has the opposite effect on instantaneous drag production during a wing beat: the Wagner effect reduces drag during the wing acceleration phase, and increases drag in the deceleration phase. As a result, the Wagner effect augments both wing braking before stroke reversal and wing accelerating after stroke reversal. This shows that the Wagner effect directly causes a reduction in the energetic cost of flight.
Intriguing is also that as a result of their opposite dynamics, the instantaneous added-mass forces and Wagner effect tend to partly cancel each other. Therefore, although the instantaneous forces of both mechanisms separately are quite high, the instantaneous forces of them combined are much smaller. Despite this partial cancellation, the acceleration-based forces near stroke reversal are still significant, for both the fruit fly and the malaria mosquito. Here, the drag component of these combined acceleration-based forces increase the power requirement for flight, suggesting that an elastic storage mechanism remains important for reducing the cost of flight (Dickinson & Lighton Reference Dickinson and Lighton1995; Reid et al. Reference Reid, Schwab, Maxcer, Peterson, Johnson and Jankauski2019).
The high acceleration-based forces at stroke reversal might also play a prominent role in flight manoeuvrability and control. Many adjustments of the wing-beat kinematics for manoeuvring occur during stroke reversal (Beatus, Guckenheimer & Cohen Reference Beatus, Guckenheimer and Cohen2015; Dickinson & Muijres Reference Dickinson and Muijres2016; Muijres et al. Reference Muijres, Iwasaki, Elzinga, Melis and Dickinson2017b), which is the phase where both the added-mass forces and Wagner effect appear prominent in both the fruit fly and the mosquito. Thus, insects can manipulate these forces to generate the aerodynamic torques required for flight control.
5. Conclusions
Using a systematic parametric study based on CFD simulations, we developed an aerodynamic force model for insect wings with a revolving stroke acceleration. This model captures the effect of the unsteady aerodynamic mechanisms, added mass and Wagner effect, on the aerodynamic forces produced by the revolving insect wing.
Using this model, we showed that added-mass forces on beating insect wings are higher than previously assumed due to two aspects: first, the fluid volume that accelerates with the flapping insect wing is larger than previously thought, resulting in high air-pressure-based added-mass forces; second, the added-mass forces oriented tangential to an insect wing are governed by viscous forces, and are consequently orders of magnitude higher than assumed using potential theory. Furthermore, our novel Wagner-effect model defined in the time domain can be applied to a quasi-steady modelling approach, allowing to study the Wagner effect in flapping insect flight.
We applied our model to the hovering wing beats of two Dipteran insect model species, the fruit fly and the malaria mosquito. This shows that stroke accelerations play an important role in the transient aerodynamic force production of these flying animals. Especially for high-frequency flappers, such as the malaria mosquito, the stroke-acceleration-based forces significantly influence the forces related to powering flight, required for weight support, and used for flight control. Due to its generality, whereby wing geometry scaling is explicitly included, our model can be applied to a wide range of flying animals, and can thus be used to generate a broader and more universal understanding in flapping flight systems.
Supplementary material
Supplementary material are available at https://doi.org/10.1017/jfm.2022.31.
Acknowledgements
We thank H. Schipper and L. Honfi Camilo for their IT support, and the PhD thesis committee members of W.G.v.V. for their constructive feedback on the topic of this paper.
Funding
F.T.M. was supported by a grant from The Netherlands Organization for Scientific Research (NWO/VI.Vidi.193.054).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available at https://doi.org/10.1017/jfm.2022.31.
Author contributions
W.G.v.V., J.L.v.L. and F.T.M. conceived and designed the study. W.G.v.V. developed the computational system, carried out the simulations and drafted the manuscript. W.G.v.V. and F.T.M. performed the data analysis. F.T.M. provided the fruit fly and mosquito data, coordinated the study and helped to draft the manuscript. All authors contributed critically to the development of the theoretical models, the writing of the manuscript, and gave final approval for publication.
Open access
