A Modified Quasisteady Aerodynamic Model for a Sub-100 mg Insect-Inspired Flapping-Wing Robot

This study proposes a modified quasisteady aerodynamic model for the sub-100-milligram insect-inspired flapping-wing robot presented by the authors in a previous paper. The model, which is based on blade-element theory, considers the aerodynamic mechanisms of circulation, dissipation, and added-mass, as well as the inertial effect. The aerodynamic force and moment acting on the wing are calculated based on the two-degree-of-freedom (2-DOF) wing kinematics of flapping and rotating. In order to validate the model, we used a binocular high-speed photography system and a customized lift measurement system to perform simultaneous measurements of the wing kinematics and the lift of the robot under different input voltages. The results of these measurements were all in close agreement with the estimates generated by the proposed model. In addition, based on the model, this study analyzes the 2-DOF flapping-wing dynamics of the robot and provides an estimate of the passive rotation—the main factor in generating lift—from the measured flapping kinematics. The analysis also reveals that the calculated rotating kinematics of the wing under different input voltages accord well with the measured rotating kinematics. We expect that the model presented here will be useful in developing a control strategy for our sub-100 mg insect-inspired flapping-wing robot.

The flapping flight of insects is characterized by a high angle of attack and a high rotational component, causing a large number of separations of the boundary layer and potentially generating vortices attached to the leading and trailing edges of the wings. These complex boundary conditions make it difficult to build an accurate steady aerodynamic model. Instead, researchers have attempted to develop quasisteady or unsteady aerodynamic models both for the flapping flight of insects and for insect-inspired robots [10][11][12][13][14][15][16][17][18][19][20][21]. For example, by studying a dynamically scaled model of the fruit fly (Drosophila melanogaster) [10], Dickinson et al. identify three mechanisms-delayed stall, rotational circulation, and wake capture-to explain how insects produce high lift at low Reynolds numbers. In a subsequent study, Sane and Dickinson [11] use blade element theory (BET) to develop a quasisteady aerodynamic model for insect flapping. From studies of free-falling cards [22][23][24], Bergou et al. [12] propose a 2-DOF quasisteady aerodynamic model for insect flapping that considers the circulation, dissipation, added mass, and inertial effect. In recent years, several improved models have been presented. For example, Nabawy and Crowther [16] propose a quasisteady aerodynamic model for the evaluation of the steady translational force coefficients of flapping wings in normal hover. The model is generic in that it can be applied to wings of arbitrary morphology and kinematics without the use of experimental data, and the aerodynamic components of the model are linked directly to morphology and kinematics via physical relationships. On this basis, Nabawy and Crowthe [17] propose a novel lifting-line theory introducing the concept of the equivalent angle of attack, which enables the capture of the steady nonlinear aerodynamics at high angles of attack and allows accurate estimation of aerodynamic forces from geometry and kinematic information alone. For detailed studies of the aerodynamic moment and the movement of the center of pressure on an insect wing, Han et al. [18] present an improved quasisteady aerodynamic model of a hawkmoth-scale flapping wing that considers the movement of the center of pressure. They find that the model becomes more suitable when the center of pressure is assumed to be at the half-chord rather than the quarter-chord. Nakata et al. [19] propose a novel CFD-informed quasisteady model which assumes that the aerodynamic forces on a flapping wing can be decomposed into quasisteady forces and parameterized based on CFD results and which is capable of predicting flappingwing aerodynamic forces and power with higher accuracy. Lee et al. [20] present an improved quasisteady aerodynamic force and power model for rigid flapping by considering the effects of Reynolds numbers, Rossby numbers, wing aspect ratios, and taper ratios. Their model has the advantage of being applicable over a wider range of flow conditions without prior tuning or calibration. Wang et al. [21] propose a predictive quasisteady model by considering the wing's translation, rotation, and coupling, as well as the addedmass effect. This model shows high accuracy in predicting the center of pressure, the aerodynamic loads, and the passive pitching motion for various Reynolds numbers without any empirical parameters.
Based on the nonrigid connection of a flexible hinge, the artificial wing used in the robot we designed generates lift through its passive rotation while flapping. The kinematics of the wing can be considered as 2-DOF of flapping and rotating. In addition, since the observed deformation of the wing during flapping is negligible, it can be regarded as a rigid plate, similar to the approximation of a twodimensional rigid plate (or wing) in Wang et al. [12,23]. However, the latter assumes that the axis of rotation of the plate is on its midline and that the chord length is maintained-characteristics not found in insect wings or in insectinspired robots.
The purpose of the present study is to provide a more accurate estimate of the lift and passive rotation generated by the robot's artificial wing. We first develop a quasisteady aerodynamic model based on BET by modifying the terms of dissipation and added-mass from [12]. We then report the findings of experiments performed to assess the model.
The wing kinematics and the lift generated by the robot were captured simultaneously by a binocular high-speed photography system and a customized lift measurement system, respectively. The lift and passive rotation were calculated both on the basis of the model and from the experimentally captured wing kinematics under different input voltages. The results predicted by our model are in close accord with the experimental measurements.
2. Description of the Artificial Wing 2.1. Morphology. The artificial wing studied in this paper is morphologically similar to the Eristalis tenax wing [25], as shown in Figure 2.
To reduce the complexity of the artificial wing, the wing root, veins, and thickness of the insect wing are ignored. The plane profile and the geometric parameters of the artificial wing are shown in Figure 3.
The x w -axis and y w -axis represent the span and chord direction of the wing, respectively, and the direction of the z w -axis corresponds to the right-hand rule. R is the length of the wing (measured as 13 mm). x r is the distance from the wing root of the wing O ′ to the center of rotation O, which is small compared to R and is taken to be 0 in this paper. The green part shown in Figure 3 is the spanwise strip of the wing, r is the distance from the spanwise strip to the root of the wing,r = r/R is defined as the normalized distance, d r is the width of the spanwise strip, and d y is the width of a chord element on the spanwise strip. P com is the center of mass of the wing, P cog is the geometric center of the wing, P A is any point on the spanwise strip, P M is the midpoint of the spanwise strip, and P r is the intersection of the spanwise strip and the x w -axis. The leading edge y l and the trailing edge y t of the wing are parameterized into functions of y l ðrÞ and y t ðrÞ, respectively, as shown in Table 1, in which the coefficients are obtained by polynomial fitting of the morphology of the Eristalis tenax wing in [25] based on MATLAB (Math-Works Inc.).
The chord length c can then be expressed as cðrÞ = y l ðrÞ − y t ðrÞ, and the area of the single wing can be integrated as S = Ð R 0 cðrÞdr (calculated as 49.27 mm 2 ). The mean chord is c = S/R (calculated as 3.79 mm), the normalized chord is defined asĉ = c/ c, and the aspect ratio is AR = R/ c (calculated as 3.43). A nondimensional radius of the kth moment of the wing area can then be expressed asr k k ðSÞ = Ð 1 0ĉ r∧ k dr, as outlined in [25]. Table 2 compares the morphological parameters of the artificial wing and the Eristalis tenax wing.

Kinematics.
The coordinate systems and the angles in 3-DOF of a wing are defined in Figure 4. The x 3 y 3 z 3 -coordinate system is defined as the experimental reference coordinate system, in which the x 3 z 3 -plane is the flapping stroke plane and the y 3 -axis is parallel to the average lift of the robot while hovering. The x w y w z w -coordinate system is defined to describe the fixed-wing plane. The flapping angle φ is defined as the opposite number of the rotation angle of the x 2 y 2 z 2coordinate system relative to the x 3 y 3 z 3 -coordinate system along the −y 3 -axis. The deviation angle θ is defined as the  where _ φ and _ ψ are the first derivative with respect to time of φ and ψ, respectively. Note that the linear velocity of each spanwise strip along the x w -axis is ignored based on BET; that is, we assume w v A,x = 0. In the x w y w z w -coordinate system, the wing inertia matrix after ignoring the thickness of the wing can be expressed as where w I zz = w I xx + w I yy and w I xy = w I yx .
The mass of the single artificial wing (m w ) is measured as 0.5 mg, and the mass distribution of the wing in the x w y w z w -coordinate system is calculated by Inventor (Autodesk Inc.) as shown in Table 3, in which w P com is the center of mass of the wing.

Modeling Based on BET
In this section, the circulation, dissipation, added mass, and inertial effect are considered in the development of a quasisteady aerodynamics model based on BET.
3.1. Circulation. As shown in Figure 5 Figure 2: Shape of (a) the artificial wing and (b) the Eristalis tenax wing from [25].
Trailing edge y t ⁎ Figure 3: The plane profile and the geometric parameters of the artificial wing. 3 Applied Bionics and Biomechanics dynamics (CFD) and experimental results of Pesavento and Wang [22], the aerodynamic force and aerodynamic moment generated by the circulation mechanism on a single spanwise strip can be expressed as If ρ air is the air density (taken as 1:29 × 10 −6 g/mm 3 ), w v M,x = 0, and Γ is the circulation, composed of the translational Γ trans and rotational Γ rot circulation, where C T and C R are dimensionless coefficients, respectively, taken as 1.8 and π [15,22], and α = arccos ð w v M,z /j w v M | Þ is the angle of attack, defined as the angle between w v M and the y w -axis of each spanwise strip. Note that since the angle of attack of each spanwise strip is different, the angle of attack is defined for each strip rather than for the entire wing.

Dissipation.
Based on [24,26], for a chord element d y d r on the spanwise strip shown in Figure 3, the aerodynamic force generated by the dissipation mechanism can be expressed as , ð7Þ     Figure 5: A single spanwise is taken as an example to define the circulation and angle of attack. The angle of attack is the angle between the velocity of the midpoint and the y w -axis of the spanwise strip; it varies between −π/2 and π/2. The angle of attack in this case is negative. 4 Applied Bionics and Biomechanics where the dimensionless coefficient C D is taken as 3.4. The aerodynamic force and aerodynamic moment generated by the dissipation mechanism on a single spanwise strip can then be integrated as 3.3. Added Mass. During wing flapping, the acceleration and deceleration of the wing cause the surrounding air to accelerate and decelerate. In this process, the surrounding air pushes against the wings, creating the added-mass effect. For the translation and rotation of thin two-dimensional rigid wings, the aerodynamic force and aerodynamic moment generated by the added-mass effect on a single spanwise strip can be expressed (based on [27]) as where w _ v r,i and w _ ω x are the first derivatives with respect to time of the linear velocity w v r of w P r and of w ω r in the x w y w z w -coordinate system, respectively. λ i and λ ij are the added-mass coefficients, defined as 3.4. Inertial Effect. Part of the aerodynamic force and moment generated by the wing is used for the linear and rotational acceleration of the wing itself. Although the mass of the wing is small, its inertial effect is not negligible. In the x w y w z w -coordinate system, the inertial force and the inertial moment of a single wing can be expressed as where w _ v com is the first derivative with respect to time of the linear velocity w v com of w P com in the x w y w z w -coordinate system.
3.5. Total Aerodynamic Force and Moment. From the above, the total aerodynamic force and the total aerodynamic moment acting on a single wing can be integrated in the x w y w z w -coordinate system as It is clear from this that the aerodynamic force and moment are completely determined in the proposed model by the morphology and 2-DOF transient kinematics (φ and ψ) of the artificial wing. In addition, some unsteady aerodynamic mechanisms, such as the start-up vortices [28], the spanwise flow [29], the wake capture [10], and the tip vortices [30], are partly included in the dimensionless aerodynamic coefficients used in this study. In this way, the aerodynamics generated by some unsteady-state mechanisms are also accounted for in the model.

Measurement of Wing Kinematics and Lift.
With the fixed-shape artificial wing designed for our robot, the flapping aerodynamics are determined only by the 2-DOF wing kinematics of φ and ψ. However, the lift generated by the robot is also valuable in assessing the model. Our experiments measured the kinematics of the wing and the lift of the robot simultaneously. As illustrated in Figure 6(a), a binocular high-speed photography system was used to measure the flapping and rotational kinematics of the left wing of the robot. At the same time, the lift generated by the robot was measured using a customized lift measurement system described in [31] and illustrated in Figures 6(b) and 6(c). A trigger was used to ensure that the photography system and the lift measurement system worked synchronously.
More specifically, the binocular high-speed photography system consisted of two high-speed cameras (Phantom LC111, working synchronously at 3000 fps), with the tips and intersections of the wing veins serving as mark points. Based on the binocular ranging method and threedimensional motion reconstruction, the spatial position of these mark points in each frame is obtained by the direct linear transformation (DLT) method [32], and the flapping and rotating kinematics are then derived from the geometric relationship of these mark points.

Applied Bionics and Biomechanics
For the lift measurement, the robot is mounted vertically on the support plate of the lift measurement system, as shown in Figure 6(c). The real-time lift generated by the robot is transformed by the deformable Invar-made doublecantilever beam into a slight displacement in the vertical direction of the target plate, and the real-time lift of the robot is then extracted by measuring the real-time displacement of the target plate using a capacitive displacement sensor (CS005, Micro-Epsilon). Since a magnet is used as part of the electromagnetic actuator in the robot, a long truss made of carbon fiber was introduced to extend the distance between the robot and the Invar-made support plate so as to prevent electromagnetic interference. After calibration, the lift measurement system had a dynamic resolution of 0.457 μN (0.82 mg), and a sensitivity of 2.19 μm/mN.
The sub-100 mg insect-inspired flapping-wing robot presented in [6] was used as the prototype for the measurement. Sinusoidal voltages with a fixed frequency of 80 Hz were applied to the robot. The measured flapping kinematics φ meas and rotating kinematics ψ meas of the left wing under different voltage amplitudes are shown in Figures 7(a) and 7(b), in which each curve is obtained from these measurement points by sixth-order Fourier fitting. The rotating kinematics are delayed overall with respect to the flapping kinematics by a quarter of a stroke cycle. As the amplitude of the input voltage is increased, the flapping and rotating kinematics of the wing change more dramatically. The relationship between the peak-to-peak amplitude of the measured wing kinematics and the input voltage amplitude is shown in Figure 7(c). It is clear that the amplitudes of the wing kinematics can be modulated by the input voltage amplitude, which is helpful in controlling the lift of the robot. Table 1 and the measured flapping and rotating kinematics of the wing into Equations (12) and (13), the aerodynamics of the single artificial wing under different input voltage amplitudes were calculated based on the proposed quasisteady aerodynamic model, as shown in Figure 8. Obviously, the magnitudes of the aerodynamic force and moment are determined by the amplitude of the input voltage. For the calculated aerodynamic force 3 F Aero,y along the y 3 -axis (i.e., the lift direction of the robot while hov-ering), the total force and circulation force almost coincide, indicating that the circulation mechanism makes the dominant contribution to the aerodynamic force. The two peaks of the total aerodynamic force are located in the middle of the upstroke and the middle of the downstroke, with the translational circulation force making the main contribution. However, at the reversal of the stroke, the total aerodynamic force is close to zero or even negative. Note that the curve should be symmetric around the stroke reversal; the difference between the two peaks of force shown in Figure 8 is mainly due to the asymmetry of the wing kinematics caused by assembly and measurement errors. Since the flapping frequency of the robot is 80 Hz, the mean aerodynamic force is more valuable for flight control than the transient force. The distribution of the mean aerodynamic force after cycle averaging is shown in Table 4. The mean force is almost entirely determined by the circulation and the added mass, while the contributions of the dissipation and the inertial effect are negligible. Unlike the aerodynamic force, the calculated aerodynamic moment is mainly determined by the dissipation, the added mass, and the inertial effect rather than by the circulation. Since both the added-mass moment and the inertial moment are proportional to the linear and angular acceleration of segments of the wing at specific points, the trends of the added-mass moment and inertial moment are similar. Both the dissipation moment and the circulation moment are greater in the middle of the upstroke and downstroke, while they are almost zero at the reversal of the stroke.

Experimental Verification.
To verify the applicability of our quasisteady aerodynamic model, we analyzed the lift and the passive rotation generated by the robot based both on the model and on experimental measurements.

Lift Estimated by the Aerodynamic
Force. The robot we designed consists of two wings, the kinematics of which are assumed to be perfectly s1ymmetric while hovering. Thus, the calculated lift generated by the robot can be expressed as To reduce the interference from mechanical and electrical noise at high frequencies, the measured instantaneous lift generated by the robot was filtered by a digital low-pass filter  Applied Bionics and Biomechanics with a cutoff frequency of 300 Hz. Figure 9 compares the measured and the calculated lift forces generated by the robot under different input voltages during one flapping cycle. The measured lift L meas and the calculated lift L calc are in good accord with regard to overall amplitude and trend, but unlike the calculated results, the measured amplitude of the upstroke is close to that of the downstroke. This indicates that the lift generated by the robot in the actual process is to some extent insensitive to the imperfect symmetry of the wing kinematics. Considering the unavoidable simplification in the mathematical model, the transient deviations between the measured and calculated results are acceptable. Both the measured curve and the theoretical curve show that the lift generated by the flapping wings fluctuates strongly, especially near the stroke reversal point, where the transient lift may even be negative. Figure 9(e) compares the calculated mean lift L calc to the measured mean lift L meas of the robot under different input voltage amplitudes after cycle averaging. The comparison demonstrates the accuracy of the proposed quasisteady aerodynamic model in estimating the lift generated by the robot.
Although the values of the aerodynamic constants C T , C R , and C D in the model are assumed from experience in our calculations, it is interesting to study the changes in the calculated lift of the robot when these three constants fluctuate separately. As described in Section 3, the total aerodynamic force acting on a single wing w F Aero is linear with C T , C R , and C D , respectively. With the fixed morphology and captured kinematics of the wing, the change in the mean lift of the robot is therefore directly proportional to the change in each aerodynamic constant. That is, where C i refers to C T , C R , and C D , respectively. In order to quantify the relationship between the relative change in the 7 Applied Bionics and Biomechanics mean lift of the robot and the relative change in each aerodynamic constant, we designate the sensitivity of the former to the latter as k i (i.e., k T , k R , and k D , respectively). Then, where C i0 is the assumed value used in Section 3 ( C T0 = 1:8, C R0 = π, and C D0 = 3:4), and L 0 is the calculated mean lift of the robot corresponding to C i0 . Considering that the range of ΔC i /C i0 is [-30%, 30%], we now use the model to calculate the corresponding relative changes in the mean lift of the robot under the input voltages from 0.1 V to 1.1 V. For simplicity, Figure 10 shows just a part of the calculated results. The calculated k T varies mostly between 0.55 and 0.65, which is 2-3 times the corresponding k R , while k D is close to zero under any input voltage. These results indicate that changes in C T and C R are more likely to cause changes in the calculated mean lift and that changes in C D have almost no effect on the lift. This finding is similar to the results illustrated in Figure 8 and Table 4; that is, the circulation (especially the translational circulation) is the dominant contributor to the calculated lift, while the contribution of the dissipation is negligible.

Passive Rotation Estimated by the Aerodynamic
Moment. As shown in Figure 2, each side of the artificial wing is connected to the transmission by a flexible hinge. Due to the nonrigid connection of the hinge, the flapping of the wing produces passive rotations, which are the main factor for the generation of lift [33,34]. In this study, the rotating kinematics of the wing are accurately obtained by the binocular highspeed photography system based on three-dimensional motion reconstruction. However, further control tests of the robot require that the rotating kinematics of the wing be estimated quickly and easily rather than measured accurately. The binocular high-speed photography system requires at least two high-speed cameras as well as nonreal-time image processing due to a large amount of data and therefore is more complex and less efficient than estimations generated by a mathematical model. The flapping kinematics of the wing can be easily measured in real time by only one high-speed camera. Based on the aerodynamic moment The wing hinge used in the robot is a sandwich structure consisting of a flexible polymer middle layer between thin, rigid outer layers, as shown in Figure 11. The torsional stiffness can be expressed as where E h (3.5 GPa) is the elastic modulus of the flexible polymer layer and l h (120 μm), w h (1800 μm), and t h (7.5 μm) are, respectively, the length, width, and thickness of this layer. From the analysis in Section 3, in the x w y w z w -coordinate system, the total moment that causes the passive rotation of the wing can be expressed as

Applied Bionics and Biomechanics
According to [35], the dynamics of the wing in 2-DOF (flapping and rotating) can be modeled as follows, based on the principle of virtual power: Substituting Equations (4), (13), and (19) into Equation (20) yields where w M A,x is (in accordance with Equation (13)) solely a function of φ and ψ given the fixed morphology of the artificial wing. Therefore, by substituting the measured flapping kinematics φ meas into Equation (21), the calculated rotating kinematics ψ calc of the wing are obtained as shown in Figure 12. The measured rotating kinematics ψ meas under the same input conditions are also plotted to verify the accuracy of the calculation results. It is clear that the calculated rotating kinematics ψ calc are close to the measured results ψ meas in terms of overall amplitude and trend, especially when the amplitude of the input voltage is relatively large. However, the amplitude of rotation is slightly overestimated by the proposed model at local positions such as the peaks in Figure 12, and the overestimation gradually increases as the input voltage (i.e., the amplitude of the flapping kinematics) decreases. The main reason is that when the flapping frequency is fixed and the flapping amplitude is reduced, the contribution of the translational aerodynamics to the drag force is further reduced due to the low flapping velocity.  10 Applied Bionics and Biomechanics Furthermore, the passive rotation of the wing is slightly restricted by the unavoidable minor shaking and bending of the wings, and these unfavorable factors become significant as the flapping amplitude of the wings is further decreased. As a result, the accuracy of the model's estimation of the passive rotation declines with the decrease in the flapping amplitude. Like the lift results, the measured rotating kinematics are more symmetric about the stroke reversal point than the calculated rotating kinematics because the mathematical model is slightly sensitive to the asymmetry of input wing kinematics. Since the flapping amplitude of the robot is generally above ±55°in actual flight tests, and the model's estimates of the passive rotation in this range are relatively accurate, the estimations are a useful resource for the control of the robot.

Conclusion
In this paper, we present a modified quasisteady aerodynamic model for an electromagnetically driven sub-100-milligram insect-inspired flapping-wing robot whose design was presented in [6]. Based on the blade-element theory, the aerodynamic mechanisms of the circulation, dissipation, added mass, and inertial effect are considered in developing the model. Since the deviation of the artificial wing is negligible for our robot, the aerodynamics generated on the wing are completely determined, in this model, by the wing morphology and the 2-DOF wing kinematics of flapping and rotating.
To verify the applicability of the model, the wing kinematics and the lift generated by the robot were measured synchronously by a binocular high-speed photography system and a customized lift measurement system. Based on the measured flapping and rotating kinematics of the artificial wing, we calculated the aerodynamic forces and moments acting on the wing under different input voltages. Our results show that the circulation mechanism makes the dominant contribution to the transient aerodynamic lift, while the contributions of the dissipation and the inertial force are minimal. However, the aerodynamic moment is mainly determined by the dissipation, added mass, and inertial effect rather than by the circulation. Combined with the synchronous lift measurement of the robot, the transient lift and the mean lift estimated by the model are all in good agreement with the measured results under different input voltages. In addition, the robot generates lift by the passive rotation of the wings. Although the passive rotation can be captured by the binocular high-speed photography system used in this study, a more convenient way to estimate it is needed. To supply this need, we analyzed the 2-DOF flapping-wing dynamics of the robot based on the proposed model and obtained estimates of the rotating kinematics from the measured flapping kinematics alone. The calculated rotating kinematics of the wing under different input voltages are all in good accord with the measured values.
We therefore conclude that the modified quasisteady aerodynamic model developed in this study is applicable with high accuracy to the sub-100 mg insect-inspired flapping-wing robot presented previously and that the model will provide theoretical support for the development of control strategy.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.

11
Applied Bionics and Biomechanics