Reduced-Order Modeling and Experimental Studies of Bilaterally Coupled Fluid–Structure Interaction in Single-Degree-of-Freedom Flapping Wings

Flapping insect wings continue to inspire several emerging technologies, such as flapping wing micro air vehicles (FWMAVs) and elastic airfoil energy harvesting devices. FWMAVs are a robotic platform [1–3] that could enable low-cost remote sensing with unprecedented spatial resolution. Airfoil-based energy harvesters have the potential for highly efficient energy extraction from ambient flows [4–6] and could power the extensive sensor networks employed in many “Internet of Things” applications. However, the mathematical models necessary to design and optimize such technologies are often inefficient and require significant computational resources. As a result, many flapping wing models are challenged by the high-dimensional parametric studies essential for engineering design.

As an artificial or biological wing flaps, it deforms from both fluid and structural forces. This fluid–structure interaction (FSI) plays a critical role in flapping wing dynamics and has been studied extensively. Many flapping wing FSI models rely on direct computational methods, such as finite element analysis (FEA) coupled to computational fluid dynamics (CFD) [7–13]. However, both CFD and FEA require considerable computational resources to estimate flapping wing dynamics, and these inefficiencies are compounded when the two computational methods are coupled together. From the structural standpoint, large flapping rotations lead to periodic centrifugal forces that cause FEA stiffness matrices to become time-varying [14]. If direct FEA is used to calculate wing deformation, the stiffness matrix must be updated at each interval of analysis. The result is a huge number of degrees-of-freedom (DOFs), and the time required to evaluate the response of all DOFs is extensive. From the fluid standpoint, CFD must resolve the flow field over an entire control volume in order to estimate the pressure distribution over the wing surface [15]. This often requires solving several thousands of equations which makes CFD computationally intensive. To reduce the computational demand of the fluid dynamic solver, others have employed unsteady vortex lattice methods (UVLMs) to investigate flapping wing FSI [16,17]. Despite lower computation times compared to CFD, UVLM is a numerical method that requires considerable resources depending on the desired solution accuracy. Thus, direct computational methods are not well suited for efficient calculations of flapping wing FSI.

To reduce computational complexity, many researchers leverage quasi-static methods rooted in blade element theory (BET) [18–21]. BET discretizes a wing into airfoils (blade elements) that run along the wing’s chord. The elemental aerodynamic forces are estimated over each individual blade using 2D airfoil theory and are then integrated over the wing to calculate net aerodynamic forces. While this is an efficient method to estimate aerodynamics, BET is generally limited to rigid wings. It has been used only a handful of times to address the effects of wing flexibility. Wang et al. developed a flapping wing FSI model based upon BET, and assumed the wing’s leading edge was rigid [22]. As a result, the model was best suited to estimate torsional deformation rather than bending deformation. Stanford et al. developed an FSI model that accounted for bending, however their structural solver required each physical DOF be solved for [23]—they did not leverage modal truncation to reduce the order of the structural model. Jankauski developed a reduced-order aeroelastic framework for flapping wings using modal truncation and BET, but this framework was used only to study one-way coupled FSI where the fluid was able to affect the structure but not vice versa [24]. It is possible that bilateral coupling between fluid and structure non-trivially affects flapping wing dynamics in some circumstances. While differences in unilateral and bilateral FSI models have not been studied extensively in flapping wings, bilateral FSI models of wind turbines predict stresses greater than those predicted by unilateral models [25]. Furthermore, in simplified cases such as a flexible plate in an ambient flow, bilateral FSI models predict physical phenomena missed entirely by unilateral models [26]. Based on these findings in other systems, bilateral coupling must be investigated more thoroughly in the context of flexible flapping wings.

Based upon this literature review, there remains a need for a reduced-order, bilaterally coupled FSI model for flapping, flexible wings. In this work, we develop this FSI model for a rectangular wing undergoing single-degree-of-freedom (SDOF) flapping. Though real insect wings have complex, three-dimensional rotational kinematics [27], these kinematics are challenging to replicate experimentally at high frequencies observed in insect flight. Further, multi-degrees-of-freedom (MDOF) kinematics give rise to complex fluid dynamic phenomena, such as rotational lift and damping [28]. Before formulating a reduced-order FSI model which considers MDOF kinematics, it is sensible to first develop and experimentally study a simpler SDOF model in order to demonstrate the feasibility of the approach. Moving forward, we will generalize this framework in order to account for the more complex flapping kinematics and flight conditions observed in real insects. The work presented here is a necessary first step toward accomplishing this goal.

The remainder of the paper is organized as follows. First, we derive the FSI model using the Lagrangian approach for the structural equation of motion (EoM) and BET for the fluid model. BET is a quasi-steady method which neglects unsteady fluid dynamic phenomena such as dynamic stall or vortex shedding. Next, we detail a simple SDOF flapping experiment used to verify our model both in-air and in-vacuum. We then compare experimental measurements to simulation results, and conclude with a discussion on how the fluid environment affects wing deformation. This paper extends the findings originally published in our ASME IDETC 2020 conference paper [29].

Here, we derive a reduced-order bilaterally coupled FSI model for flexible wings subject to SDOF flapping. We begin by determining the structural EoM via the Lagrangian method. We then identify aerodynamic forces and coupling through a BET approach. Aerodynamic terms are included in the EoM using the principle of virtual work.

The FSI framework in this section originated in Refs. [24,30] for a wing rotating in three dimensions, though these previous studies considered only unilateral fluid–structure coupling. We now consider SDOF rotation but with bilateral fluid–structure coupling. The model is summarized briefly to provide clarity to this paper. For a more thorough treatment, the reader is directed to these references.

Now, we determine the aerodynamic modal force Q_k using a BET formulation. We assume the following:

• Unsteady fluid dynamic forces are negligible.

• Drag is the only aerodynamic force that contributes to wing deformation because the induced fluid velocity is always normal to the wing’s surface.

• Inertial and aerodynamic forces do not vary along the wing chord and are assumed to act at the chord’s centroid.

• The displacement of each vibration mode retained does not vary with respect to the chord. This implies that the wing cannot exhibit a torsional response.

In this section, we describe a simple experiment designed to study our SDOF FSI model. We construct an SDOF rotation stage to prescribe flapping kinematics to a rectangular paper wing. Mode shapes and natural frequencies of the paper wing are estimated via FEA and are subsequently verified using a scanning vibrometer. During flapping experiments, we measure the spanwise strain at a point near the base of the wing using a uniaxial strain gage. We conduct flapping experiments both in-vacuum and in-air to verify the isolated structural and FSI models independently.

The SDOF rotation stage is pictured in Fig. 2, and a video of the stage during operation is included in the Supplementary Material on the ASME Digital Collection. All mounting brackets are 3D printed with FormLabs durable resin. A 60 W DC motor (Maxon Motors, 310007) drives the motion of the wing. The motor is equipped with an optical encoder that provides position feedback to a motor controller/driver (Maxon Motors, EPOS 24/5). The motor controller uses a proportional-integral-derivative framework to maintain prescribed flapping kinematics and minimize overshoot. All motion profiles are prescribed through a laptop computer running LabVIEW. In this work, we consider discrete flapping frequencies ranging from 5 to 15 Hz and a rotation amplitude of 45 deg. All rotations are sinusoidal. Each trial at a particular flapping frequency is conducted three times and the measurements from each trial are averaged in the frequency domain.

The motor connects to a wing clamp through a shaft coupler. The clamp secures the wing edge. A 350 ohm strain gage (Omega Engineering, GD-2/350-DY11) is adhered near the wing base. We use a National Instruments NI 9236 cDAQ module to provide excitation voltage to the gage as well as to record the temporal strain during experiments. The wing clamp is terminated with a low friction flange mount ball bearing. A female-end quantized analog encoder (US Digital, MAE3-A10-250-220-7-B) records the angular position of the terminated shaft end. The entire rotation stage is housed in an acrylic vacuum chamber (Sanatron, Fig. 3) capable of operating at pressures as low as 500 milliTorr. At this pressure, the medium density is roughly 0.05% of ambient air. All vacuum feed-through components are provided by Kurt J. Lesker company. The ability to conduct experiments in-vacuum allows us to evaluate the accuracy of the structural model prior to investigating the FSI model.

We use a simple rectangular paper wing in all flapping experiments. The wing is made of thick card stock and is cut with a shear. All material and geometric properties of the wing and the strain gage mounted to the wing are shown in Table 1. We model the experimental wing in abaqus fea to determine its natural frequencies and mode shapes. The FEA model assumes the wing is clamped at its base edge (Fig. 4) which implies no rotation or translation in this clamped region. We include the strain gage in the FEA model because it has a thickness on the same order of magnitude as that of paper. According to the manufacturer, the gage is composed primarily of polyimide film. As a result, the gage locally stiffens the wing in a way that cannot be neglected. The model is discretized into 250 elements, which we found was sufficient for convergence of the first natural frequency. For this work, we retain only a single vibration mode. Across the experimental parameters considered, higher-order modes had a negligible contribution to the wing’s dynamic response. The first natural frequency predicted via FEA is ω₁ = 31.5 Hz and corresponds to a bending mode (Fig. 6).

Next, we verify FEA-predicted mode shapes and natural frequencies experimentally. Because the wing is lightweight and has a large surface area, we measure these parameters in-air as well as in-vacuum to remove added mass effects. We secure the paper wing to a modal shaker (Modal Shop, K2007E007) using a metal clamp. The shaker excites the wing at its base via a linear swept sine signal ranging from 10 to 1000 Hz over 3.2 s. We measure basal excitation with a piezoelectric accelerometer (PCB Piezotronics, 352A21) and the response velocity of the wing at several points using a planar scanning vibrometer (Polytec PSV-400). We acquire data at 2.56 kHz, which results in a spectral resolution of 3200 FFT lines over the frequency range considered. We average the frequency response function over three trials at each measurement point to reduce spectral noise. The FRF averaged over the surface of the wing is shown in Fig. 5. Measured responses are reconstructed to identify the first vibration mode shape. This mode shape agrees well with that determined via FEA (Fig. 6). We then calculate the frequency response function averaged over the wing surface G(ω) and use FEMTools modal parameter extractor to estimate the first natural frequency and damping ratio from this averaged frequency response (Table 2).

Overall, the agreement between natural frequencies calculated via FEA (ω₁ = 31.5 Hz) and measured experimentally (ω₁ = 30.2 Hz) is good. We consider the natural frequency measured in-vacuum for this comparison. The discrepancy is likely due to uncertainty in Young’s modulus of the paper wing. We assume a Young’s modulus of 9.5 GPa for paper (Table 1), whereas reported values range between roughly 7 and 12 GPa [32]. To minimize the uncertainty in material properties, we use measured natural frequencies rather than those determined via FEA for all simulations that follow. This effectively adjusts Young’s modulus of the FEA model to reconcile differences between it and the physical structure.

The natural frequency in-air is slightly lower than that measured in-vacuum due to the added mass. We also observe that the damping ratio is greater in-air, which suggests that aerodynamic damping may affect the structural response during flapping experiments. We found no notable differences between the mode shape measured in-air and in-vacuum. For the simulations that follow, we use experimentally measured natural frequencies and damping ratios rather than those determined via FEA.

In this section, we compare model predictions to measurements taken from a wing flapping in-vacuum and in-air. We then use the FSI model to gain insight into the mechanisms responsible for wing deformation.

We first evaluate the accuracy our structural model (Eq. (5)) without aerodynamics (Q_k = 0). We compare model predictions of strain to those measured from the wing flapping in-vacuum. Equation (5) is solved numerically over 50 periods (where each period is discretized into 100 uniform time-steps) to estimate the wing’s modal response. We selected this temporal discretization such that the effective sampling frequency is 20 times greater than the fifth harmonic of the flapping frequency. Increasing the time-steps per period from 100 to 500 maximally affected peak wing strains by less than 0.5%, so we maintained the 100 time-steps per period in order to facilitate efficient parametric studies. We consider 100 evenly spaced flapping frequencies from 5 to 15 Hz, which is the same flapping frequency range considered in the experiment. Strain at the location of the gage is determined through Eq. (6). We calculate the Fourier transform of strain numerically and identify peak-to-peak magnitude at the driving frequency and each significant harmonic thereof. Across the range of flap frequencies ω considered, we observe appreciable response components at ω and 3ω. We see a large 5ω response for flapping frequencies between 5 and 7 Hz, however we do not focus on these responses here given that most insects flap at 30% or higher of their wing’s first natural frequency [33].

Then, the magnitude of ω and 3ω strain components as a function of flapping frequency ω is shown in Fig. 7. In general, the model predicts experimental findings fairly well. However, the predicted strain at ω is slightly higher than what is measured experimentally (particularly when flapping frequencies exceed 10 Hz), and the peak 3ω response occurs at a lower frequency for the model than for the experiment. We conjecture these small discrepancies stem a weak hardening nonlinearity, which may occur in thin plates undergoing large displacements [34].

Note the significant 3ω response when ω = 10 Hz. The strain response at 3ω has nearly the same magnitude as the primary ω response. For the model, the 3ω response is a result of the inertial force (which has frequency content only at ω) interacting with the time-varying wing stiffness (Eq. (5)). Because 3ω is near the fundamental frequency of the wing, the dynamic response at 3ω is large. Thus, even though the wing model is linear, the time-variance of the stiffness coefficient may still elicit superharmonic resonances [35].

Now that we have verified that the structural model predicts in vacuo dynamics with reasonable accuracy, we repeat the flapping experiment in air. We include the aerodynamic modal forces given by Eq. (14) into the EoM assuming the aerodynamic drag coefficient is C = 3.0 and the density of air is ρ_f = 1.22 kg/m³. The aerodynamic drag coefficient and air density were not measured explicitly and are instead approximated from Refs. [19,32]. Please note that inaccuracies in these parameters may introduce systematic bias into the model, however potential errors appear to be small in the context of this work. All other parameters are identical to those presented for in vacuo simulations. FSI model predictions and experimental results are compared in Fig. 8. In general, the model-theory agreement is very good. We find the primary ω strain response increases modestly in air, while the 3ω strain response is substantially attenuated. We believe the increased ω response is due to aerodynamic loading, whereas the decreased 3ω response is due to aerodynamic damping. This is discussed in greater detail in Sec. 4.2.

In terms of computational efficiency, we are able to predict the response over a flapping cycle in about 50 ms with matlab’s ODE45 solver. Unfortunately, we were unable to find any reported computation times of direct FSI methods for comparison. However, in the past we have used CFD to calculate pressure distributions over a rigid flapping wing [30]. Using an equivalent time-step and coarse surface mesh, it required approximately one hour per wingbeat to resolve to the flow field without considering fluid–structure coupling. Both computations are made using the same custom workstation with the following hardware: Intel Core i9-9900K Coffee Lake 8-Core, 16-Thread, 3.6 GHz processor, CORSAIR Vengeance LPX 32GB (2 × 16GB) 288-Pin DDR4 SDRAM, Gigabyte Z390 Aorus PRO WIFI LGA 1151 (300 Series) Intel Z390 HDMI SATA 6 Gb/s USB 3.1 ATX Intel Motherboard, Seagate BarraCuda ST2000DM008 2TB 7200 RPM 256MB Cache SATA 6.0Gb/s 3.5″ Hard Drive, and Samsung 970 Evo PLUS 2TB Internal Solid State Drive. Based on these evaluation times, we estimate our new model predicts the wing response at least four orders of magnitude faster than direct FSI methods. In reality, the computational savings are likely much greater.

In Sec. 4.1, we observed a large 3ω dynamic response when the wing flapped at roughly one-third of its first natural frequency. The 3ω response was pronounced in vacuum but significantly attenuated in air. Here, we aim to identify the aerodynamic mechanism responsible for attenuating the 3ω response. There are two plausible explanations by which aerodynamics will reduce the 3ω response. The first mechanism is aerodynamic loading, Q_A, which does not rely on structural deformation. If Q_A has a 3ω component, it is possible that it will either constructively or destructively interfere with the in vacuo 3ω response depending on their relative phase. The second mechanism is through aerodynamic damping, Q_ζ, which relies on bilateral coupling of fluid and structure. This term is less straightforward to analyze because the structural response must be known in order to calculate it.

In order to identify which mechanism is responsible for attenuating the 3ω response in air, we simulate both unilaterally and bilaterally coupled FSI models. The unilateral model includes only aerodynamic loading, Q_A, whereas the bilateral model considers both aerodynamic loading Q_A and aerodynamic damping Q_ζ. We use the same parameters and flapping frequency that were used in previous simulations. Strain magnitude as a function of flap frequency for both unilateral and bilateral FSI models is shown in Fig. 9. Both the unilateral and bilateral FSI models suggest that aerodynamic loading Q_A increases the ω response magnitude relative to the in vacuo case. The unilateral model also predicts that strain magnitude at 3ω will approximately double from the in vacuo case as a result of aerodynamic loading, which is not consistent to what we observed experimentally. On the other hand, the bilateral FSI accurately predicts the experimentally measured 3ω response which suggests that aerodynamic damping Q_ζ is responsible for attenuating the 3ω response.

To explore this further, we calculate the magnitude of Q_A and Q_ζ individually, as well as their sum, as a function of ω (Fig. 10). Q_A has a nontrivial 3ω component that leads to the large 3ω response predicted by the unilateral FSI model in Fig. 9. In contrast, when Q_A and Q_ζ are considered together, their 3ω components interact destructively due to a difference in phase. When the wing is flapping at approximately one-third of its natural frequency, the net aerodynamic force has no 3ω component. We note that Q_ζ is dependent on the modal response and cannot be treated as a simple time-dependent input, and thus the coupling of the fluid and structure is ultimately what attenuates the 3ω component of the net aerodynamic force.

Additionally, Fig. 10 indicates that over the range of flapping frequencies considered, Q_A is always greater than Q_ζ and that the two forces interfere destructively. As a result, the magnitude of Q_A by itself is always greater than the magnitude of Q_A and Q_ζ summed together. At flapping frequencies of approximately 12 Hz and higher, Q_ζ grows faster than Q_A at ω. Consequently, the net aerodynamic forces at ω shrink as flapping frequencies exceed 12 Hz. The net aerodynamic forces at 3ω shrink as flapping frequencies exceed 13 Hz. However, it is difficult to see the reduction in aerodynamic loads at high frequency in the experimental results because inertial forces tend to dominate the structural response for the wing considered. Inertial forces increase monotonically with respect to flapping frequency, which implies wing deformation in general increases with flapping frequency despite reductions in aerodynamic forces.

The two primary insights identified by our flapping FSI model and experiment are (1) a 3ω superharmonic resonance when the wing flaps at one-third of its natural frequency in vacuo, and (2) significant attenuation of this superharmonic resonance in air due to aerodynamic damping. Superharmonic resonance and aerodynamic damping have been addressed using other analytic [36,37] and computational [38,39] FSI models, however we believe the framework presented here provides additional understanding into the mechanisms that underlie them.

The 3ω superharmonic resonance in flapping wings is often attributed to a weak cubic stiffness term in the wing’s structural equation of motion [36,39]. Weak nonlinear cubic stiffness can generate odd harmonics for a system excited by a harmonic force at ω, and if one of these harmonics coincides with the system’s resonant frequency, a superharmonic resonance occurs [40]. However, linear time-varying stiffness, a characteristic of flapping wings much less discussed in the literature, contributes to superharmonic resonance as well. Time-varying stiffness arises from centrifugal softening, where centrifugal softening relies on coupling between the flapping angular velocity and wing elastic deformation [14]. Physically, centrifugal softening implies the perceived stiffness of the wing is at a minimum as it passes its mid-stroke (maximum angular velocity) and at a maximum upon stroke reversal (zero angular velocity). Like a system with cubic stiffness, a system with linear time-varying stiffness excited at a single input frequency will exhibit a response at that input frequency and odd harmonics thereof. In the present work, the linear time-varying approximation yields good agreement between the experiment and theory which leads us to believe time-variance is primarily responsible for the superharmonic resonance observed here. Nonetheless, we saw secondary effects of strain hardening behavior during in vacuo experiments. In reality, it is likely that both time-variance and structural nonlinearity contribute to flapping wing dynamics—additional efforts must be made to identify under which conditions one is dominant over the other.

Our FSI model also provides an approximate expression for aerodynamic damping, which relies on the interplay of the wing’s rotational velocity and structural deformation. Because wing kinematics are prescribed, the aerodynamic damping term is also linear and time-varying. Kang and Shyy utilized a similar linear time-varying analytic expression for aerodynamic damping, though the flapping kinematics considered in their study are different than those considered here [37]. In contrast to this linear time-varying form, Ramananarivo et al. employed a nonlinear aerodynamic damping model that does not rely on the wing’s rigid body motion [36]. Despite the differences in the two models, both Kang and Ramananarivo found in their respective studies that aerodynamic damping caused a phase lag in the wing’s structural response that was beneficial to aerodynamic thrust production. In the present study, we identify another critical function of aerodynamic damping: to attenuate superharmonic resonances. This may be important in the context of insect flight, where several insects flap around one-third the fundamental frequency of their wings [33]. While a modest 3ω wing response is thought to be beneficial to flapping wing flight in terms of energy efficiency [35] and aerodynamic performance [39], an extremely large 3ω resonant response would eventually compromise lift production or damage the wing. It is plausible that aerodynamic damping, which scales linearly with rotation amplitude according to our model, is responsible for maintaining beneficial deformation amplitudes at higher order harmonics of the flapping frequency. This must be further investigated via FSI models that can account for MDOF flapping kinematics.

Artificial and biological wings deform under both aerodynamic and inertial forces. This deformation is believed to play an important role in flight and has a notable effect on aerodynamic force production and energy efficiency. However, most flapping wing FSI models are high-order and require long solution times, which hinders their ability to carry out parameter studies efficiently. Lower-order models may employ assumptions that limit their accuracy, for example assuming unilateral coupling between fluid and structure.

Here, we formulate a reduced-order bilaterally coupled FSI model of a wing undergoing SDOF flapping. The structural model is derived via the Lagrangian approach and general fluid loading is accounted for via the principle of virtual work. We then estimate specific fluid forces using a quasi-steady BET approach. We find that the fluid forces can be represented as two terms: (1) a time-dependent fluid loading term that is a function only of rigid body flapping kinematics and (2) a fluid damping term that is a function both of the wing’s rigid body motion and rate of elastic deformation. The fluid damping term effectively couples fluid and structure.

To study this model experimentally, we construct a mechanical flapper and measure the basal strain of a paper wing flapping at 45 deg with frequencies from 5 to 15 Hz. We initially flap the wing in-vacuo to benchmark the structural model. We find that when flapping at one-third the wing’s natural frequency, the wing experiences a large superharmonic resonance. In air, this superharmonic response is attenuated substantially. Through simulation, we determine that aerodynamic damping is responsible for attenuating this superharmonic resonant response. The developed models predict the strain response magnitude with good accuracy for both in-air and in-vacuo experimental studies. In addition to its accuracy, the model is very efficient and can be solved in milliseconds.

However, the model and experiment here are limited to SDOF flapping. We also neglected unsteady fluid dynamic phenomena such as vortex shedding or flow separation. It is unclear if these assumptions will be justified in the context of real insect flight, where the wing structure and MDOF flapping kinematics are significantly more complex. Nonetheless, we have demonstrated here that a reduced-order, bilaterally coupled flapping wing FSI model performs well in a simplified case. This motivates possible extensions of this general modeling approach to consider more realistic flight conditions.

This research was supported by the National Science Foundation under Award No. CBET-1855383. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. We would like to thank Sanatron LLC for providing us with the vacuum chamber used in our experiments.