## Abstract

In the analysis of flexible flapping wings of insects, the aerodynamic outcome depends on the combined structural dynamics and unsteady fluid physics. Because the wing shape and hence the resulting effective angle of attack are *a priori* unknown, predicting aerodynamic performance is challenging. Here, we show that a coupled aerodynamics/structural dynamics model can be established for hovering, based on a linear beam equation with the Morison equation to account for both added mass and aerodynamic damping effects. Lift strongly depends on the instantaneous angle of attack, resulting from passive pitch associated with wing deformation. We show that both instantaneous wing deformation and lift can be predicted in a much simplified framework. Moreover, our analysis suggests that resulting wing kinematics can be explained by the interplay between acceleration-related and aerodynamic damping forces. Interestingly, while both forces combine to create a high angle of attack resulting in high lift around the midstroke, they offset each other for phase control at the end of the stroke.

## 1. Introduction

Flapping motions are common to biological flyers. Insects flap their wings to generate lift to stay aloft. For decades, the aerodynamics of insect flight remained inexplicable. A well-known example is the myth that bumblebees cannot fly according to the classical, aerodynamics theories [1]. Since then, findings on unsteady lift enhancing mechanisms, such as the clap and fling [2], dynamic stall via leading-edge (LE) vortices [3], wing–wake interaction via wake-capture [4,5] and rotational forces owing to combined translation and rotation [4,6], have uncovered the flapping wing aerodynamics and contributed to our understanding of how an insect may stay aloft [7].

One intriguing feature of biological flight that has not been adequately explored is the role of the wing flexibility on insect aerodynamics. Insect wings are flexible and can substantially deform during flight [8–10]. Several studies have suggested that insect wing rotations may be passive [11,12], meaning that the resulting rotation is owing to a dynamic balance between the wing inertial force, elastic restoring force and fluid dynamic force.

For flexible flapping wings, wing shape and motion are *a priori* unknown and complex as the unsteady aerodynamics and structural dynamics are coupled and a satisfactory analysis needs to take both effects into account. The instantaneous shape and rates of shape change affect force generation. Simultaneously, fluid forces influence wing deformations. To date, most studies on insect flight have analysed the resulting aerodynamics based on measured or mimicked wing kinematics [4,13–15]. For example, by measuring the wing shape deformations of a locust wing, better aerodynamic performance could be correlated to the flexibility of the wing [13]. However, such approaches are based on the measured, resulting wing kinematics rather than by modelling how the fluid and structure interact.

While high-fidelity simulations [16–19] and experimental work [20,21] can yield instantaneous information of the flow field and wing deformations of flexible flapping wings in hover flight, analytical formulation of the lift and wing deformation as a function of time has not been developed for flexible flapping wings in hover. Such an analytical procedure not only offers a faster way to analyse the aerodynamics compared with computational or experimental methods, but also directly evinces the relations between different physical mechanisms and offers a framework to better understand the nature of fluid–structure interactions (FSIs).

Recently, several scaling analyses have successfully modelled the lift generation of a flexible wing, e.g. [21,22], by capturing the first-order, important mechanism in the FSI. For example, we introduced scaling relationships that provide time-averaged estimates for the lift, power input and efficiency [22]. The proposed scaling was consistent with relationships based on experiments in air [21] and in water (appendix A), covering a wide range of fluid densities [22].

In general, the fluid dynamic force acting on a moving wing can be decomposed into two terms: (i) the acceleration–reaction force, often linearized by the added mass [23] and (ii) the force induced by the vorticity in the flow field [22,24], or the aerodynamic damping term. For example, the rigid wing models by Theodorsen [25] for forward flight and Dickinson *et al.* [4] for hover flight also include both forces.

However, the aforementioned scaling analyses account only for one of the two components. For example, the scaling that we derived [22] mainly accounts for the acceleration-related forces, i.e. added mass and wing inertial force, and neglects the aerodynamic damping; furthermore, the effect of shape deformation on the fluid dynamics was not explicitly addressed. Consequently, the scaling analysis only gives an estimate of the time-averaged quantities.

The instantaneous information of both wing deformation and lift can be critical for the understanding of biological flight, especially in regard to offering input to flight stability and control. For insects, with the notable exception of butterflies, the instantaneous information from flapping may not be as critical for simple flight dynamics models, because flapping time is much faster than that of the entire flyer [7]. Regardless, better instantaneous aerodynamic force models are needed to improve the time-averaged estimate and our understanding of insect flight. Moreover, the dynamic flight stability is an inherent property of a flapping insect [26] with three types of flight control in general [26] (i) stabilization control to stabilize inherently unstable flight by applying time-varying control forces and moments; (ii) manoeuvre control to produce fast turns by applying large, time-varying control forces and moments and (iii) steady-state control, for example, to change from hover to forward flight. Because the dynamics of a flapping flyer is the result of coupling between unsteady aerodynamics and time-varying body and wing motions, an accurate analysis of stability and dynamics requires instantaneous information of how wing deformation and lift change in time.

Here, we show that under the assumption of small wing deformation, lift and drag can be decoupled in the two-way FSI. The combined fluid–structure coupling simplifies to a dynamic balance between chordwise wing deformation and drag. We model the drag on the wing with the Morison equation [27], accounting for both added mass and aerodynamic damping. The resulting partial differential equation can be solved for the instantaneous wing deformation. Wing deformations produce a passive pitch with an angle of attack *α*, a key aerodynamic metric [16], which directly influences the generation of aerodynamics force. Specifically, instantaneous lift follows from the quasi-steady model [4,16] based on the corresponding angle of attack owing to wing deformations.

Figure 1 highlights the outcome of the model proposed by showing three representative angles of attack owing to passive pitch *α* and lift coefficients *C*_{L} as a function of time. Both *α* and *C*_{L} agree excellently with high-fidelity numerical results [16]. The high-fidelity solutions were obtained from a well-validated Navier–Stokes equation solver, fully coupled with a structural dynamics solver [16,22]. The Reynolds number is relevant for fruit flies: *Re* = *Uc*/*ν* = 100, based on the midstroke velocity *U*, wing chord length *c* and kinematic viscosity of fluid *ν*. A sinusoidal plunge motion *h* with amplitude *h*_{a} and frequency *f* is imposed on the LE of the wing as *h*(*t*) = *h*_{a}cos(2π*ft*) as a function of time *t* (figure 1*h*). In the absence of freestream for hovering flight, we set the maximum translational velocity *U* of the flat plate at LE as the reference velocity, such that *U* = 2π*fh*_{a} [15,16,28]. The reduced frequency becomes a geometrical quantity: *k* = π*fc*/*U* = *c/*(2*h*_{a}) and is inversely proportional to the Keulegan–Carpenter number *K*_{C}, often used in discussion of the Morison equation [29], *K*_{C} = π/*k*.

Suffice it to say that numerous other important questions remain, for example, how insects can optimize energy consumption using flexible wings. These issues are beyond the scope of this study and are considered to be unresolved issues.

## 2. Material and methods

### 2.1. Dimensionless parameters for the fluid–structure interaction and case set-up

The dimensionless parameter for the FSI is the frequency ratio *f*/*f*_{1}, where *f*_{1} is the first natural frequency of the wing in chordwise direction. For the high-fidelity numerical solution, the ratio between the structural density *ρ*_{s} and fluid density *ρ*_{f} was set to *ρ** = *ρ*_{s}/*ρ*_{f} = 7.8, and the thickness of the wing *h*_{s} was 2% of the chord, i.e. the thickness ratio is *h*_{s}* = *h*_{s}/*c* = 0.02. The density ratio is similar to steel in water or a light material in air. The case set-up in the high-fidelity simulations was motivated by rigid wing experiments in mineral oil [4–6,30]. For insects, wing bending is mainly owing to wing inertia as the density ratio *ρ** is typically of the order *O*(10^{3}) [22,31]. The mineral oil has a density that is 7 × 10^{2} times greater than air, such that the acceleration–reaction force, which is often linearized by the added mass, becomes the dominant force that deforms the wing [22]. Still, the study of the aeroelastic response of insect wings in such a low density ratio system is justified as long as the resulting aeroelasticity matches the actual motion [30] as illustrated by a rigid wing with *ρ** = 1.4 [4], crane fly [12] or shown by the γ-scaling analysis [22]. More importantly, the presented analytical model is not restricted to specific values of density ratios. In addition, a recent study suggests that flexible bending for animal locomotion is universal in air and water [32].

In the high-fidelity analysis, we varied the reduced frequency *k* and the frequency ratio *f*/*f*_{1} by changing the plunge amplitude *h*_{a} and Young's modulus of the wing *E*, respectively, in order to probe their influence on the resulting aerodynamics and the structural deformations [16]. The range selection of *h*_{a} was motivated from biological flyers [33], such that 0.25 < *k* < 2.35.

For most natural flyers, the first bending mode is the most dominant mode for the deformation of a flapping wing [21,22,34]. In this study, we consider *f*/*f*_{1} ≤ 0.4, such that the resulting phase lag covers advanced rotation and symmetric rotation, which are critical for rigid and flexible wings [4,16]. Biological relevance of these parameters and case set-up were discussed in our previous work [16]. A table of data points can be found in the electronic supplementary material.

While three-dimensional effects, such as spanwise flow that seem to stabilize the LEVs [35] or LEV–tip–vortex interaction [36] are notable in general, the chordwise flexibility is essential and warrants an independent study [16,37]. Furthermore, at *Re* = *O*(10^{2}), the effects of spanwise flow on the overall aerodynamics are less important than they are at higher Reynolds numbers [35,38]. In addition, the characteristics of the LEVs in two-dimensions for plunging motions are representative of three-dimensional flapping wings as long as the stroke-to-chord ratio is within the range of typical insects, i.e. around 4 to 5 [33,39,40], which we consider in this study.

### 2.2. Fluid–structure interaction of flexible wings in hover

We consider hovering motions. In a hover flight, the flyer maintains its position and, therefore, the forward flight velocity is zero. In other words, the freestream velocity is zero. The main purpose of plunging motions in a hover flight is to generate a force, which is lift, to sustain the weight of the flyer, see for example the seminal work by Dickinson *et al.* [4] on the aerodynamics of hovering insects. The projection of the resulting fluid dynamic force in the direction of the plunge motion is called drag (figure 1*h*).

We consider a non-dimensional flow field with unit density initiated by a two-dimensional flat plate with unit chord and flat edges [16]. The wing is a linear elastic flat plate of uniform thickness *h*_{s}, density *ρ*_{s} and Young's modulus *E*. The wing is oriented vertically. As the flat plate follows the imposed horizontal motion *h*(*t*) at the LE, the resulting fluid dynamic force dynamically balances with the wing inertia and the elastic bending forces, modelled locally as a linear Euler–Bernoulli beam as
2.1where *v* is the displacement owing to bending motion, *Π*_{0} = *ρ***h*_{s}*(*k*/π)^{2} is the effective inertia, the inertia of the wing normalized by the fluid dynamic variables [22], *c*_{d} is the non-dimensional structural damping coefficient, *Π*_{1} = *Eh*_{s}*^{3}/(12*ρ*_{f}*U*^{2}) is the effective stiffness, the wing stiffness normalized by the fluid dynamic variables [22] and *F* is the uniformly distributed transverse fluid force per unit length on the wing, such that *F** = *F*/(*ρ*_{f}*U*^{2}). The superscript (asterisk) indicates non-dimensional variables. The dimensional variables are non-dimensionalized with *U*, *c* and the period of a motion stroke *T* = 1/*f*, respectively.

The FSI in flexible flapping wings is a two-way coupled system: the fluid dynamic forces influence the wing deformations and, simultaneously, the wing deformations affect fluid dynamic forces. The force *F* on a real insect wing is governed by the unsteady Navier–Stokes equations. However, Navier–Stokes equations are nonlinear and a closed, analytical expression does not exist for *F* in general. We discuss in §2.3, how we can approximate *F* by identifying the first-order mechanisms, maintaining the two-way coupling.

Zero initial conditions are assumed. Remaining boundary conditions are that the angle at the LE is zero and that the trailing edge (TE) is a free end. The inhomogeneous boundary condition at the LE can be homogenized by considering the wing deformation relative to the LE, *w* = *v* – *h* [22]. The resulting wing camber deformations *w* can also be regarded as an angle between the TE and LE leading to the angle of attack *α*(*t**) (figure 1*h*). The acceleration term owing to *h* becomes the inertial force.

### 2.3. Small deformation assumption and decoupling of lift and drag

One of the assumptions of the Euler–Bernoulli beam model is that d*w*/d*s* is assumed to be small, where *s* is the coordinate along the wing (figure 1*i*). Therefore, the angle *θ* between TE and LE also remains small, because . As a consequence, the fluid dynamic force can be assumed to act normal to the wing, i.e. only the drag *D* has the principal contribution to the wing deformation. The contribution from lift is much smaller than from drag and, therefore, neglected in the dynamic balance of the FSI system.

Based on scaling arguments, we found that, for the deformation of flexible flapping wings, the fluid dynamics force term could be well approximated by the added mass force for various cases at *k* > *O*(1) and *Re* > *O*(10^{2}) [22]. Moreover, based on equation (2.1), we derived scaling relationships between the aerodynamic performance, such as the time-averaged propulsive force and the propulsive efficiency, and the shape deformation parameter *γ* defined as
2.2where *γ* can be obtained from *a priori* determined parameters. The Strouhal number, *St* = *fh*_{a}/*U*, another important parameter in flapping wing aerodynamics, becomes a constant for the chosen velocity scale *U* in hovering flight. A physical interpretation of *γ* is a length scale for the resulting TE deformation relative to LE, normalized by the chord [22]. Therefore, the small wing deformation assumption translates to , which corresponds to .

The drag on a moving body, immersed in viscous fluid can be described by the semi-empirical Morison equation [27,29]. The Morison equation includes both added mass and aerodynamic damping. The added mass is proportional to wing acceleration. The aerodynamic damping term is proportional to the square of the wing velocity. Because we assume small wing deformations, the drag on the wing is dominated by the imposed plunge motion *h*. As a result, we neglect the higher-order contributions from *w* on drag *D.* The resulting fluid dynamic force *F* in equation (2.1) in terms of the resulting Morison equation becomes
2.3The added mass was modelled with the classic solution for the force acting normal to a thin flat plate with chord *c* moving with the position *h* [22]. The aerodynamic damping coefficient *C*_{v} depends on *k*, *Re*, etc., in general [29] and may also be a function of space and time. Experimental measurements on an oscillating rigid flat plate in a water tank suggest that *C*_{v} = 15(*k*/π)^{0.5} exp(1.88/*Re*^{0.547}), applicable for 1.01 < *Re* < 1057 and 0.67 < *k* < 2 [29]. In the original formulation, the period number was used, which is equal to the Keulegan–Carpenter number, instead of the reduced frequency, i.e. *P*_{e} = *K*_{C} = π/*k.* The Reynolds number of *Re* = 100 is in the range of applicability, and the majority of the high-fidelity solutions are in the applicable range of *k*. The magnitude of *C*_{v} increases with *k* and when the relative importance of the viscous effects becomes more important, i.e. as the Reynolds number decreases.

The quadratic term of the aerodynamic damping yields the additional damping term *ḣẇ*, because the square of the local wing velocity is . The structural damping of an insect wing is a material and complex property, which is still unknown. Thus, an analytical treatment of structural damping of insect wings is difficult. To enable analysis, we assume the damping coefficient *c*_{d} in equation (2.1) as a constant, combining the structural damping and additional aerodynamic damping *ḣẇ*. We assume *c*_{d} = *x*_{1}*k*^{2} + *x*_{2}*k* + *x*_{3}, where we have included the reduced frequency terms to account for *ḣ*. The coefficients *x*_{1}, *x*_{2} and *x*_{3} are determined by minimizing the difference of the phase lag *φ* between the prediction from the current model and the high-fidelity numerical result [16] using a least-squares method. The shape deformation parameter *γ*, derived based on added mass, was shown to scale the time-averaged lift and midstroke angle of attack [16]. However, added mass alone could not predict the phase lag, and we will show that aerodynamic damping should also be considered. Note that the high-fidelity results were parametrized with *k* and *f*/*f*_{1}, see also the table in the electronic supplementary material for the values of *k*. The resulting coefficients are *x*_{1} = 0.83, *x*_{2} = −0.21 and *x*_{3} = 0.16. Dependence of *c*_{d} on other parameters, such as *f*/*f*_{1} or *Re*, and the use of other empirical fit models will be investigated in the future.

The resulting Euler–Bernoulli beam equation for the structure coupled with the Morison equation for the fluid models the two-way coupled FSI system for flexible flapping wing in hover flight. The main assumption is that the wing deformations are small.

## 3. Results and discussion

### 3.1. Passive pitch angles of the flexible flapping wing

The solution of the coupled system equation (2.1) with the fluid dynamic force term given by the Morison equation (2.3) can be obtained using the method of separation for the wing deformation relative to the imposed motion *w*(*x**, *t**) = *v*(*x**,*t**) − *h*(*t**) = *χ*(*x**)*τ*(*t**) for the spatial component *χ* and the temporal component *τ*. The solution procedure will not be detailed in this paper as it is straightforward and similar to our previous investigation [22]. The solution of the spatial component is the same as that of the classical vibration of a clamped beam with a free end [22]. The temporal part results in an ordinary differential equation for an oscillator for the first half-stroke 0 < *t** ≤ 0.5 as
3.1where the non-dimensional damping factor is *ζ* = *c*_{d}/*П*_{0} and the non-dimensional natural frequency is with *k*_{1} ≈ 1.875 being the eigenvalue belonging to the first spatial mode *χ*_{1}. We assume a symmetry for the second half-stroke 0.5 < *t** ≤ 1. The first term in the RHS of equation (3.1) models the effects of the added mass and the inertial force of the wing, which are both acceleration-related. The second term represents the aerodynamic damping term, and the third term is owing to the structural damping. The coefficients *δ* and *β* are defined as and *β* = *C*_{v}/*Π*_{0}.

The solution of equation (3.1) is
3.2with
where *A* and *B* are coefficients that can be determined from the initial conditions. Because the transient term *T*_{transient} of equation (3.2) will be exponentially damped out, we discuss only the effects of the remaining steady-state terms. The acceleration-related term *T*_{acceleration} has the same origin as *γ* and accounts for the added mass and the wing inertial force. For a system without structural damping, i.e. *ζ* = 0, the proportionality coefficient of this term becomes *γ* = *δ*/(4*π*^{2} + *ω*_{f}^{2}). The effects of this term were discussed in our previous work [22].

Analysis of the time history of the wing deformations from the high-fidelity solution showed that wing deformations can be approximated by a first-order harmonic [16]. The departure from the first-order sinusoidal motion increased with increasing wing deformations. Because we assume small deformations, we expect the resulting wing motion to be close to a first-order harmonic. Furthermore, the imposed zero initial condition is not the same as the condition in the high-fidelity solutions: high-fidelity solutions showed nearly periodic passive pitch motions [16]. To account for the difference in the initial condition, we first determine the wing deformations at the middle and at the end of the strokes, *w*_{m} = *w*(0.25) and *w*_{e} = *w*(0.5), which are converted to angles of attack owing to passive pitch *α*_{m} = arctan(*w*_{m}) and *α*_{e} = arctan(*w*_{e}), see also figure 1*h*. Based on these two angles, a first-order harmonic approximation can be constructed for the angle of attack *α*_{FH}(*t**) = 90° − *α*_{a}cos(2*πt** + *φ*) by solving for the phase lag *φ* and the angular amplitude *α*_{a} with *α*_{FH}(0.25) = *α*_{m} and *α*_{FH}(0.5) = *α*_{e} [16].

For *f*/*f*_{1} < 0.2, the angle of attack owing to passive pitch at the end-of-the stroke *α*_{e} increases with *f*/*f*_{1}, after which *α*_{e} reduces (figure 2*a*). On the other hand, the midstroke angle *α*_{m} decreases with *f*/*f*_{1} (figure 2*b*). The resulting *φ* strongly correlates to *f*/*f*_{1}, yielding the advanced, symmetric and delayed rotation modes with increasing *f*/*f*_{1} (figure 2*c*). The rotational mode becomes symmetric at around *f*/*f*_{1} = 0.35, corresponding to the maximum lift [16]. Despite its complicated behaviour in *f*/*f*_{1}, the predicted *α*_{a} is in reasonable agreement with the high-fidelity numerical results (figure 2*d*). Overall, the predicted angles owing to deformations show strong agreement with the high-fidelity solutions [16] and are able to capture the critical trend of the relation between the different rotational modes and the frequency ratio.

The angles *α*_{m} and *α*_{a} in figure 2*b,d* show a complicated behaviour relative to *f*/*f*_{1}. The scattered trend is an indication that *f*/*f*_{1} is not a good scaling parameter for these two angles. A better scaling parameter is *γ* [16] as shown in figure 3, illustrating a better correlation. The angles *α*_{m} and *α*_{e} are obtained from the dynamic balance of the wing and the fluid under the small deformation assumption. As a result, the error, measured as the absolute difference between the estimated and high-fidelity solutions, Δ*α*_{m} and Δ*α*_{e}, increases with *γ* (figure 3*c*,*d*).

While the midstroke angle *α*_{m} and amplitude *α*_{a} scale with *γ*, the correct scaling for the angle owing to passive pitch at the end-of-the stroke *α*_{e} and the phase lag *φ* were unknown [16]. A potential reason that the shape deformation parameter *γ* was not able to scale *φ* is that *γ* was derived from the added mass and inertial force, not accounting for the aerodynamic damping. As discussed further in §3.3, the good agreement in *α*_{e} and the phase lag *φ* (figure 2*a*,*c*) suggests that contribution of aerodynamic damping is essential for the estimation of *φ* and the instantaneous wing deformation.

### 3.2. Lift on the flexible flapping wing

In the model for the dynamic balance between the wing deformations and the fluid dynamic forces, higher-order contributions of the relative wing deformation *w* and d*w*/d*s* on drag were neglected. However, the relative wing deformation *w* and the corresponding passive pitch angle are essential for the generation of lift. In addition to the imposed plunge motion *h*, instantaneous lift depends on the instantaneous wing shape and the rates of wing deformation: passive pitch velocity and acceleration.

For rigid wings, the lift on a vertically oriented wing in hover is expected to be zero in general. Only under certain conditions, a pure plunge motion without pitch in fluid in rest can yield positive lift force [41]. The resulting vortex dynamics around the wing breaks the imposed symmetry and generates a non-zero force that is perpendicular to the motion, the lift. For flexible wings, the symmetry is broken by the relative wing deformations.

In principle, a Navier–Stokes equation solver can be used with the obtained angle of attack due to passive pitch *α*_{FH} and the imposed plunge *h* to obtain a high-fidelity instantaneous lift. In this study, we apply the quasi-steady model for lift on a rigid wing [4,16] based on the plunge motion (*h*) and the pitch motion (*α*_{FH}) including their first- and second-order time derivatives to provide an analytical model for instantaneous lift generation on a flexible wing in hover. The differences in the quasi-steady analysis and high-fidelity numerical solution were discussed in our previous work [16].

A quasi-steady model can be valuable in predicting the lift generation on a flexible flapping wing if a proper account can be given to address the coupling between unsteady aerodynamics and structural dynamics. This is the case even though such a model assumes that the fluid dynamic forces depend on the instantaneous wing velocities and accelerations, and can under-represent, for example, the influence from the nonlinear wing–wake interaction relevant in certain cases [4,16]. Because a flexible wing can streamline its shape to mitigate the lift-degrading wing–wake interaction, the lack of an explicit model for the wing–wake interaction is not as severe as for the rigid wings [16]. However, the commonly employed quasi-steady model neglects the influence from the wing camber deformations [16].

On balance, the present model can predict both the instantaneous angle of attack *α*_{FH}(*t**) and the corresponding lift *C*_{L}(*t**), matching well the high-fidelity solutions (figure 1). As a result, the overall level of agreement of the time-averaged lift coefficient is favourable (figure 4*a*). The normalization is [16].

It is noted that as the wing deformations increase, the effects of *w* on the Morison equation become non-negligible and the predictive capability of this model worsens as expected. The error in the prediction of lift as a function of time, indicated with the mean-squared error of lift, MSE(*C*_{L}) = Σ* ^{N}*(

*C*

_{L,num}−

*C*

_{L,pred})

^{2}/

*N*initially increases with

*γ*(figure 4

*b*). There were 240 time steps per half-stroke in the high-fidelity computation [16], i.e.

*N*= 240. However, MSE(

*C*

_{L}) does not show a linear trend with

*γ*. The resulting error also depends on the error introduced in the quasi-steady model for lift as evident from the difference in resulting lift (figure 1

*f*), for the case where the wing deformations, predicted from the current model, matches well with the high-fidelity numerical results (figure 1

*c*). Furthermore, we empirically determined the coefficient for the structural and aerodynamic damping terms. In addition, as the first-order harmonic representation of the passive pitch neglects the higher-order transient response of the wing structure, which can contribute to the observed differences in the resulting lift [16].

### 3.3. Interplay between acceleration-related and aerodynamic damping forces

For dynamically scaled, rigid, fruit fly wings with active rotations, a midstroke angle of attack of 40° < *α*_{m} < 50° generated the optimal lift when the wing rotation is slightly advanced relative to the translation [4,6]. For these advanced rotation modes, the lift was slightly higher than at symmetric rotation owing to the rotational force [4,6], whereas the delayed rotation performed the worst [4]. For passive pitch of a deformed wing, the rotational force is not as dominant owing to the compliant nature of the wing [16,42]. As a result, the advanced rotation mode does not necessarily produce a higher lift (figure 4*c*). For the case considered, the highest lift was obtained when the passive pitch is nearly symmetric; this finding is very consistent with measurements of hovering fruit flies [14], honeybees [15], hoverflies [43] and hymenopterans in general [9].

Analysis of individual contributions of the acceleration-related force, i.e. wing inertial force and added mass, and aerodynamic damping evince the FSI of the passive pitch rotation for the considered cases. The acceleration-related force has the greatest contribution at lower *f*/*f*_{1} for *α*_{e} (figure 5*a*). However, as *f*/*f*_{1} increases, the aerodynamic damping starts to play a more dominant role. Around *f*/*f*_{1} = 0.35, corresponding to symmetric rotation (figure 4*c*) and optimal lift [16] with a lift coefficient of 1.8, both forces offset each other. The observation that the highest lift is obtained at a subresonance frequency ratio is consistent with numerous experimental measurements [21,44,45]. At the midstroke, where the translational velocity is at its maximum, both acceleration and aerodynamic damping forces work in the same direction, causing larger wing deformations and, therefore, higher angle of attack owing to passive pitch *α*_{m} (figure 5*b*), resulting in a high lift coefficient. While the translational acceleration is zero at the midstroke, the acceleration from the resulting passive wing rotation can help generating higher lift.

## 4. Concluding remarks

The fundamental limitation of the simplified model is the lack of information of the key input parameters, such as the lift coefficient as a function of angle of attack in the quasi-steady model [4], the aerodynamic damping coefficient *C*_{v} in the Morison equation [29] or the damping coefficient *c*_{d}. These coefficients are empirically determined from various experiments. The major merit of this study is that based on the insight learned from the previously conducted high-fidelity simulations, along with the empirical observations made in the literature, we have substantially enhanced the usefulness of analytical models by providing much improved predictive capability. Our results help to understand the complex, and intriguing physics of flexible flapping wings of insects.

Of course, a real insect wing is three-dimensional and intriguingly complex. For example, some insects show figure-8, O-shaped or U-shaped flapping patterns [4,9,14,46]. Moreover, the wing motions are three-dimensional with wing rotation at the wing root [14]. In this study, we focus on the role of chord wise flexibility and passive pitch in two-dimensional plunge motion. It is noted that the structural damping of insects is difficult to measure owing to their complicated wing structures, and an empirical fit for the structural damping term needs to be obtained from high-fidelity numerical simulations or physical experiments. Development of analytical models towards more realistic wing kinematics [46] with nonlinear beam models with higher accuracies for larger wing deformations is planned as future work.

The present model highlights the principal mechanisms of flapping flexible wing aerodynamics. A key feature is that we account for both acceleration-related and aerodynamic damping forces with the Morison equation, resulting in an analytical model for instantaneous wing deformation and lift. Our results suggest that aerodynamic damping is important for an accurate prediction of phase lag of passive pitch as both forces work in opposite directions at stroke ends.

By assuming small deformations, lift and drag can be decoupled, enabling an analytical analysis. Small deformation and small angle approximations have been widely used in different branches of mechanics. Euler–Bernoulli beam model in structural dynamics and thin-aerofoil theories in aerodynamics are well-known examples that have made great contribution to our understanding of mechanics. The solution of the governing partial differential equation provides needed information of both amplitude and phase of instantaneous wing deformation, leading to an analytical model for instantaneous lift. While errors of the model increase with wing deformations, in practice, even for a large passive wing rotation amplitude of 45°, which is beyond the small angle regime, the qualitative trend matches the high-fidelity solution. In addition, it is well documented in the literature [22,37] that if the wing becomes too flexible, the performance deteriorates.

The proposed model can be used as a fast, accurate model providing information needed to understand the flight dynamics of insects and flapping micro flying robots.

## Data accessibility

A list of data points is uploaded as the electronic supplementary material.

## Acknowledgements

C.K. thanks Dr Keith Hollingsworth of the University of Alabama in Huntsville for the fruitful discussions on the Morison equation.

## Appendix A

The γ-scaling relationship [22] between the time-averaged force on the flexible flapping wing and *γ* was shown to be consistent with the relationship that was based on an experiment in air [20,21]. The scaling is also consistent with recent experiments performed in water [47]. The consistency between the γ-scaling and the water experiments can be shown as follows. The γ-scaling was
A1where *γ* becomes equivalent to the elastoinertial number for air [22]. More recently, a similar experiment has been performed in water [47] from which a scaling relationship for forward flight was observed as
A2where AR is the aspect ratio of the wing. The frequency ratio *f** was empirically determined by locating the maximum wing tip deflection, which includes the effects of added mass. It can be shown that the empirically determined scaling for flapping wings in water, equation (A 2), is also consistent with the scaling that we proposed, equation (A 1), because
for the LHS of equation (A 1), where *S* = AR *c*^{2} is the area of a rectangular wing.

In our study, we considered the added mass as an external fluid dynamic force, which is included in the definition of *γ*. Alternatively, because added mass is proportional to the acceleration, a modified natural frequency can be established that includes the effects of added mass, i.e. the empirical frequency ratio *f** [47].

- Received August 19, 2014.
- Accepted September 17, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.