Camera-based deformation measurement of an insect-like micro air vehicle flapping wing

Abstract

Wing kinematics and wing deformation have a strong effect on the lift production and aerodynamic efficiency of bio-inspired flapping wing aerial vehicles. The goal of this study is to automatically track and classify wing kinematics components and wing deformation components of an insect-like flapping wing. Insight into wing motion and deformation characteristics may allow for improved performance of flapping wing drones. This work uses a stereoscopic high-speed camera setup and feature detection algorithm to measure the kinematics and deformation components of an insect-like flapping wing system with passive wing inclination. The wing membrane is painted with a black and white speckle pattern. Tracking and tracing of features on the wing surface is used to recover rigid wing kinematics and deformation of the entire wing surface throughout the wing stroke. Measurement results correspond with the characteristics of natural and artificial flapping wing systems in literature. The wing performs a figure-of-eight motion. The wing exhibits large deformations that vary throughout its stroke, related to both aerodynamic forces and inertial effects. Stroke velocity appears to influence angle of attack as well as spanwise bending, camber and twist.

Keywords

flapping wing micro air vehicles FWMAV bio-inspired flexible wings high-speed camera feature detection

Introduction

Flapping wing fliers in nature, such as hummingbirds and insects, show impressive flight controllability at high aerodynamic efficiency.^1–3 These promising aerodynamic characteristics have motivated researchers to develop artificial flapping wing micro and nano aerial vehicles (FWMAV and FWNAV).^4–7 Several of these FWAVs are capable of lifting off the ground and performing controlled flight maneuvers. Flight times vary from a few seconds to several minutes.

Numerical simulations^2,4,8,9 and experimental studies^1,10–12 on (artificial) flapping insect wings indicate that both wing kinematics and wing deformation play an important role in aerodynamic force production and aerodynamic efficiency. Investigations into elastic flapping wing structures often study the effect of wing deformation by comparing force production and power consumption to a simulation of an identical rigid wing structure undergoing the same wing motion. The consensus is that spanwise bending, camber and twist have a beneficial effect on aerodynamic performance. Elastic wings outperform rigid wings in situations characterized by high angles of attack, such as insect-like flapping flight.^1,9,13 The improved performance is attributed to a more stable attachment of the leading edge vortex (LEV), which prevents flow separation.^4,9,14 Percin et al.¹¹ found that wing deformation is necessary for the wings to benefit from lift enhancing phenomena such as wake capture.

Several experiments have captured the wing motion and deformation of natural and artificial flapping wings. Camera experiments on natural fliers often involve manually tracking wing features frame by frame.^15–17 This is a time consuming process, difficult for use in sensitivity studies.^6,16 Other studies draw markers, usually a few dozen,^6,10 onto the wing surface and track them either manually or by (semi-)automatic computer tracking algorithms.¹⁸

This study applies a speckle pattern to the wing surface and tracks wing features through feature detection. Use of a stereoscopic camera setup allows for a three dimensional reconstruction of wing kinematics and deformation. Visualizing the wing kinematics and wing deformation may inform improvements in driveline design and wing design. Improved design may result in higher lift capacity, increased flight endurance and flight control.

Materials and methods

This work uses an experimental investigation of wing kinematics and wing deformation on an artificial insect-like flapping wing. Figure 1 shows the artificial wing. The wing membrane is a $10 μ$ m Mylar sheet. A carbon fiber skeleton provides structural stiffness. The carbon fiber skeleton consists of four veins: a leading edge (LE), root edge (RE) and two inner veins (R1 and R2). Table 1 lists the wing parameters. Figure 2 shows the flapping wing driveline. The wing is actuated through a direct drive.⁷ A gear transmission connects the wing to a brushless DC motor that drives the stroke motion. A rubber band positioned at the wing root generates a restoring torque that tries to rotate the wing towards a zero inclination angle. The inclination motion is passive and is the result of an equilibrium between aerodynamic forces and torque from the rubber band. The wing rotates around its first carbon fiber vein. The rubber band at the wing base ensures a positive angle of attack, so this first carbon vein always corresponds to the leading edge. Torsion springs are attached at the wing base. When the wing stroke angle increases, the torsion springs store elastic energy, which is converted and released as kinetic energy at the start of the next stroke. This way the torsion springs assist in stroke reversal. The spring stiffness is chosen such that the natural frequency of the wing driveline is close to the stroke frequency of 20Hz. The resulting resonance is beneficial to driveline efficiency.^5,7,15 Stroke frequency and stroke amplitude are determined through a control circuit that consists of a ’Seeed Studio XIAO nRF52840 Sense’ processing unit and a printed circuit board that was designed in-house. A LiPo battery supplies power to the control circuit and motor.

Figure 1.

The wing has a skeleton with four carbon fiber rods attached at the wing base. Two outer rods form the wing leading edge (LE) and root edge (RE). Two inner rods (R1 and R2) add stiffness to the wing. A Mylar sheet covers the carbon fibre skeleton and forms the wing membrane. L is the length of the wing membrane, measured along the leading edge. The chord length c is defined as the distance from the leading edge to the trailing edge (TE), measured perpendicular to the leading edge.

Figure 2.

Experimental setup (left), with close-up of driveline and speckle pattern wing (right).

Table 1.

Wing parameters.

Parameter	Symbol	Value	Unit
Wing length	$L$	73.5	mm
Single wing span	$R$	83.5	mm
Mean chord length	$\bar{c}$	20.0	mm
Aspect ratio		4.17	–
Single wing surface area	$S$	1670	mm $^{2}$
LE diameter	$d_{L E}$	0.64	mm
R1 diameter	$d_{R 1}$	0.28	mm
R2 diameter	$d_{R 2}$	0.28	mm
RE diameter	$d_{R E}$	0.42	mm
Membrane thickness	$t$	10	$μ$ m
Wing mass	$m$	0.225	g

Experimental setup

Figure 2 shows the experimental setup. The setup uses two cameras (Ximea xiB-64 BC120RGCM-X8G3 and JAI SP-12000M-CXP4) that capture wing motion. The transparent wing membrane is covered with a layer of white spray paint and a black speckle pattern of approximately 800 dots with a diameter of approximately 1 mm. Speckle patterns have been used in the past to determine insect-like flapping wing characteristics. Wu et al.¹⁸ used a speckle pattern flapping wing in vacuum to distinguish between inertial and aerodynamic contributions to wing deformation. Doan et al.¹⁹ used a speckle pattern wing to perform modal analysis on an insect-like wing structure. The effect of spray paint on wing inertia and wing structural properties is negligible. The painting process increases wing mass by 0.012 grams or approximately 5%.

Motion tracking

The camera image acquisition is synchronized by a trigger unit. Frame rate is set to 800 Hz. At a stroke frequency of 20Hz, 40 images are taken per stroke cycle or 20 images per stroke. Exposure time is set to 100 $μ$ s to limit motion blur. The measurement is repeated with the driveline in various orientations with respect to the cameras to capture the full wing motion pattern. The two camera setup triangulates the coordinates of feature points on the wing membrane throughout wing stroke in a fixed reference frame $x_{0} y_{0} z_{0}$ connected to the driveline (cfr. Figure 3a). Image processing is done in MATLAB R2022a.

Figure 3.

(a) Illustration of wing kinematic parameters. Wing kinematics consist of three rotations, expressed with respect to a fixed coordinate system $x_{0} y_{0} z_{0}$ . The center of $x_{0} y_{0} z_{0}$ corresponds to the wing center of rotation. $x_{0}$ is parallel to the wing when it is at midstroke. $x_{0} y_{0}$ spans the body horizontal plane. $x_{c} y_{c} z_{c}$ follows the wing motion. The stroke angle $ϕ$ is a rotation of the wing inside the body horizontal plane, around $z_{0}$ . The deviation angle $δ$ is a rotation around the rotated y-axis. Deviation is the wing out-of-plane motion. The inclination angle $α$ is a rotation around the $x_{c}$ -axis. The inclination angle determines the angle of incidence. At zero inclination angle, the wing surface is parallel to $z_{0}$ . R is the single wing span, measured along the leading edge from the wing center of rotation to the wing tip. (b) Illustration of wing deformation parameters. $x_{c} z_{c}$ corresponds to the rigid wing plane. Spanwise bending is defined as the leading edge deformation perpendicular to the wing plane, along $y_{c}$ . The twist angle $θ_{t w}$ is the angle between $x_{c} z_{c}$ and a chord line that connects the local leading edge and trailing edge. Camber $q_{c}$ is the distance to the local wing chord, measured perpendicular to the chord line going from leading edge to trailing edge.

Wing kinematics

Rigid wing motion requires three angular parameters: stroke $ϕ$ , deviation $δ$ and inclination $α$ . Figure 3a illustrates the wing kinematics. Wing stroke is a periodic rotation around $z_{0}$ . In studies of birds and insects, these two strokes are classed as an upstroke and downstroke respectively.¹⁵ Previous in-house experiments show that wing motion is symmetric.⁷ Therefore there is no need for a distinction between upstroke or downstroke, so only one stroke is analyzed. Wing deviation is the wing motion out of the body horizontal plane $x_{0} y_{0}$ . The deviation motion determines whether the wing follows a planar,^2,8,9,14,20 elliptical,^10,14,21 figure-of-eight^6,7 or double figure-of-eight⁷ profile. All rotations are defined as rotations around positive coordinate axes, so a negative deviation angle corresponds to a wing that is tilted upward. Finally inclination is the wing rotation around its leading edge.

Determining stroke, deviation and inclination from tracked point data is relatively straightforward. The wing rotates around its leading edge, so points p on the leading edge have the highest value $p_{z 0}$ relative to their distance from the center of rotation (CoR), according to Figure 3 and with $p = (p_{x 0}, p_{y 0}, p_{z 0})^{T}$ . A regression line is drawn through a selection of detected points with the highest value for $p_{z 0}$ relative to their distance from the CoR. This line is identified as the leading edge. During operation, points on the leading edge undergo stroke motion and deviation motion. Inclination motion does not affect their position. After stroke and deviation, a point at a distance $r$ along the leading edge has coordinates:

R_{z} (ϕ) R_{y} (δ) [\begin{matrix} r \\ 0 \\ 0 \end{matrix}] = [\begin{matrix} r \cos ϕ \cos δ \\ r \sin ϕ \cos δ \\ - r \sin δ \end{matrix}]

(1)

$R_{z} (ϕ)$ and $R_{y} (δ)$ are the rotation matrices related to stroke and deviation respectively. From coordinates of a point p on the leading edge, the stroke angle and deviation angle are found:

{\begin{matrix} ϕ = \arctan (\frac{p_{y 0}}{p_{x 0}}) \\ δ = \arcsin (\frac{- p_{z 0}}{\sqrt{p_{x 0}^{2} + p_{y 0}^{2} + p_{z 0}^{2}}}) \end{matrix}

(2)

To determine inclination, a rigid wing plane is fitted through the marker coordinate data. The rigid wing plane corresponds to the $x_{c} z_{c}$ -plane, which contains the wing membrane when no wing deformation is present. A subset of tracked markers closest to the wing CoR is used to fit the rigid wing plane. Points close to the CoR are best suited for this fit, since close to the wing center the effect of wing deformation is expected to be minor. The inclination angle is then the angle between the rigid wing plane and a vertical plane that contains $x_{c}$ .

Wing deformation

Three components together describe wing deformation: spanwise bending, camber (chordwise bending) and twist (torsion). Figure 3b illustrates each deformation component.

The wing membrane is divided into wing chords. Chords are thin strips of wing membrane oriented perpendicular from the LE to the TE, as shown in Figure 3b. Tracked points are assigned to wing chords based on their distance from the CoR, measured along the LE. For each chord, the LE and TE are computed based on the average of a subset of the three highest and three lowest detected points within that chord.

Spanwise bending $q_{s b}$ is the deformation of the leading edge perpendicular to the rigid wing plane. For each chord, spanwise bending corresponds to the $y_{c}$ -value of its LE coordinates. Positive spanwise bending means the LE lags behind the rigid wing plane. Negative spanwise bending means the LE is ahead of the rigid wing plane.

Camber $q_{c}$ is the distance between a point on the wing membrane and the wing chord line. Positive camber means the point on the wing membrane lags behind the rigid wing plane. Negative camber means the point is ahead of the rigid wing plane.

Twist $θ_{t w}$ is the angle between the rigid wing plane and a straight line drawn from LE to TE. Positive twist means the wing chord is pitched down with respect to the rigid wing plane, so the angle of attack is lower compared to the case of a rigid wing. Negative twist means that the wing chord is pitched up with respect to the rigid wing plane, so the angle of attack is higher compared to the case of a rigid wing.

Per wing chord there is a single value for spanwise bending, a single value for twist and camber varies along the chord.

Results

Kinematics measurements

Figures 4 and 5 show the measured wing kinematics. Figure 4 shows the evolution of the deviation angle and the inclination angle throughout wing stroke. Figure 5 shows the wing tip path. The wing motion has a stroke amplitude of 120°. The deviation angle ranges from $- 10$ °to 8°. Based on a polynomial fit, the mean deviation angle is equal to $-$ 3°. So on average the wing is tilted upward. Mean inclination angle is equal to 29°and maximum inclination angle is equal to 50°.

At the start of wing stroke, the wing has a high angle of attack and moves downward. The high angle of attack is approximately constant until midstroke ( $ϕ =$ 0°).

As the wing approaches midstroke, the wing moves upward.

After midstroke the wing pitches down, so the angle of attack decreases. The decrease in angle of attack is likely due to stronger aerodynamic forces acting on the wing membrane, caused by high stroke velocity.

Near the end of stroke, the wing stroke velocity decreases. The wing pitches back towards zero inclination angle. The wing also moves back to its deviation angle from the beginning of stroke.

Figure 4.

Evolution of (a) deviation angle and (b) inclination angle throughout stroke. The wing moves from left to right.

Figure 5.

Bold dots show the wing tip path, moving from left to right. The data points are mirrored (transparent dots) with respect to $y_{0} z_{0}$ to clarify the figure-eight motion profile.

Figure 6 shows a side view of the wing motion from start of stroke ( $\hat{t} = 0$ ) to end of stroke ( $\hat{t} = 1$ ). This figure confirms that the inclination angle is close to zero at the start of stroke, then remains small during the first half of stroke. The inclination angle is largest after midstroke and returns to zero at end of stroke. The wing initially dips down, then tilts upward as it nears midstroke and comes back down near the end of stroke.

Figure 6.

Side view of wing motion.

Deformation measurements

Figure 7 shows the evolution of each wing deformation component throughout wing stroke.

At the start of stroke ( $ϕ = - 47.9$ ° $\dots$ $- 23.2$ °), inertial effects dominate wing deformation. While the wing does not yet experience a high stroke velocity, there are significant spanwise bending, negative camber and high twist values close to the wing tip. The leading edge is set in motion, while the rest of the wing membrane appears to lag behind, causing the strong deformation that is also visible in Figure 6 at $\hat{t} = 0.1$ . The wing moves at low velocity, yet spanwise bending is significant. Over the wing span, bending increases gradually from 0 mm at the wing root to 3.7 mm at the wing tip. There is simultaneously a large positive camber present between the root edge and inner vein R1, and a large negative camber between R1 and the leading edge. The maximum negative camber is $- 3.2$ mm. Relative to the mean chord length, this equates to $q_{c} / \bar{c} = - 16 %$ . The wing experiences a high twist that increases along the wing span up to 26°.

Towards midstroke ( $ϕ = - 16.2$ ° $\dots$ 10.5°), the effect of aerodynamic forces acting on the wing membrane becomes more prominent. This is reflected by the strong positive camber, up to 2.7 mm or $q_{c} / \bar{c} = 13 %$ . The largest camber is situated between R1 and the leading edge. Spanwise bending appears rather steady and remains concentrated toward the wing tip. Twist values are lower than at the start of stroke. Maximum twist is at the wing tip, at 19°.

Deformation after midstroke ( $ϕ =$ 24.0° $\dots$ 33.2°) is smaller than before midstroke. Both spanwise bending and camber values decrease, likely because of the increased inclination angle that lowers the wing angle of attack and thereby aerodynamic forces acting on the wing membrane. Twist values increase on the outer half of the wing, with a maximum of 23°at the wing tip. The increase in wing tip twist is expected. Near the wing tip structural stiffness is low, so the wing has a low resistance to bending. At the same time stroke velocity is highest at the wing tip, causing high aerodynamic forces. The combination of high aerodynamic forces and low stiffness should indeed result in an increased wing tip deformation.

As the wing nears the end of stroke ( $ϕ =$ 47.9° $\dots$ 54.9°), stroke velocity decreases. At the same time inclination angle is high, so angle of attack is relatively low compared to the first half stroke (cfr. Figure 4). Lower stroke velocity and lower angle of attack result in small deformation.

At the very end of stroke ( $ϕ \approx$ 63°), the wing rapidly decelerates and inertial effects again take over. This is reflected by negative spanwise bending and a large region of negative camber. The leading edge and wing membrane are ahead of the rigid wing plane.

Figure 7.

Evolution of deformation components throughout stroke: (a) Nondimensional spanwise bending; (b) Nondimensional camber; (c) Twist, in chronological order from top to bottom.

Discussion

The measured wing kinematics follow a logical evolution that agrees with findings from literature.^6,7,16 Roshanbin et al.⁶ measured a similar wing tip trajectory for their flapping wing mechanism, where the wing initially appears to strongly deviate down, then comes back up toward midstroke and returns to its mean value near end of stroke. Maeda et al.¹⁶ measured the same motion profile on a hummingbird in hovering flight. Previous in-house measurements by Timmermans⁷ show a double figure-of-eight pattern when using the same driveline mechanism at a lower stroke amplitude. The wing exhibits a large initial oscillation followed by a second oscillation at lower amplitude, creating a double figure-of-eight pattern. Aerodynamic forces acting on the wing membrane may explain this change in deviation motion profile. At higher stroke velocity, the increase in aerodynamic force acting on the wing membrane prevents the wing from deviating downward a second time.

Wu et al.¹⁸ measured total structural deformation of an artificial flapping wing reinforced with internal veins, which is closest to the wing design that is used in this study. Although they did not separate deformation contributions into spanwise bending, camber and twist, a comparison can still be made. The wing deformation in this work follows a very similar evolution to that measured by Wu et al. Along the leading edge, the out-of-plane deformation increases gradually on the part of the wing that is closest to the wing root and then increases more rapidly towards the wing tip. As stroke velocity increases, deformation related to camber and twist is highest at the trailing edge near the wing tip, where the structural stiffness is lowest. Wu et al. measured a total deformation of 18% compared to chord length. The overshoot related to inertial effects at the end of stroke also agrees with the findings of Wu et al. Camber values between 10% and 25% of mean chord length are common in elastic insect-like flapping wing systems.^13,16,17,21 Camber increases along the wing span. A maximum in positive camber is observed near midstroke. All these observations in literature agree with the findings from the previous section. Hummingbirds and artificial flapping wing systems without internal venation experience higher twist values, with maxima up to 60°.^10,15,16

The measurements from the previous section successfully characterize individual wing kinematic components and deformation components. Measurement values and profiles correspond to findings from literature. The experimental setup and processing methodology are suited for additional investigations into wing kinematics and deformation. Future research efforts should explore the effects of wing actuation parameters such as stroke amplitude and stroke frequency and wing design parameters such as vein diameter and venation pattern on wing deformation and aerodynamic performance. Also, the wing structural resonance frequency should be considered during the wing design process, as flapping wings appear more efficient when actuated close to their structural resonance frequency.^4,8 Insights from such extensive investigations should prove useful for the wing design process and lead to improved lift production and power requirements of FWAVs. Experiments should be performed under hovering conditions as well as conditions that mimic forward flight, in order to reveal the effects of an additional incoming airflow on wing kinematics and wing deformation. Elastic wing models offer opportunities in the construction of flight controllers that account for the effects of wing deformation. Such models may be especially useful to FWAVs with passive inclination motion, as changes in driving parameters such as stroke amplitude influence or maneuvers such as forward flight also affect the inclination motion. This can affect aerodynamic force production in a non-obvious way.³ The influence of a control action or flight maneuver on force production is therefore especially difficult to predict in FWAVs with passive wing motion.

Conclusion

This work uses a stereoscopic high-speed camera setup to track wing motion and deformation of an artificial insect-like flapping wing mechanism. The wing tip undergoes a periodic figure-of-eight motion. Wing inclination is largest after midstroke, when stroke velocity is highest. The wing exhibits large deformation throughout the stroke. When stroke velocity is high, aerodynamic forces likely dominate the wing deformation. At and around stroke reversal, inertial effects are more prominent. Both the kinematics and deformation correspond well with findings from literature.

The large deformations are expected to have a strong effect on aerodynamic force production and aerodynamic efficiency. Follow-up experiments should be performed that include measurements of lift and power consumption. A sensitivity analysis may be performed into the effect of actuation and wing design on wing kinematics and wing deformation and linked with lift production and power consumption. Insights gathered during additional experimental campaigns may lead to improved wing designs with better aerodynamic performance and to a low-computational cost aerodynamic force production model that accounts for wing deformation.

Footnotes

Acknowledgments

The control circuit that drives wing stroke motion was designed by Tom Henskens, KU Leuven, Celestijnenlaan 300 box 2420, Leuven, Belgium.

ORCID iDs

Thomas Roelandt

Thijs Willems

Frank Naets

Dirk Vandepitte

Funding

The author(s) received no financial support for the research, authorship and/or publication of this article: The Research Foundation - Flanders (FWO) is gratefully acknowledged for its support through research grant numbers G095120N and 1S60222N. Internal Funds KU Leuven are gratefully acknowledged for their support.

Conflicting interests

The author(s) declare that there are no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Young

Walker

Bomphrey

, et al. Details of insect wing design and deformation enhance function and flight efficiency. Science 2009; 325: 1549–1552.

Engels

Wehmann

Lehmann

. Three-dimensional wing structure attenuates aerodynamic efficiency in flapping fly wings. J R Soc Interface 2020; 17: 20190804.

Menzer

Ren

Guo

, et al. Wing kinematics and unsteady aerodynamics of a hummingbird pure yawing maneuver. Biomimetics 2022; 7: 115.

Liu

Ravi

Kolomenskiy

, et al. Biomechanics and biomimetics in insect-inspired flight systems. Philos Trans R Soc Lond B Biol Sci 2016; 371: 20150390.

Phan

Park

. Insect-inspired, tailless, hover-capable flapping-wing robots: recent progress, challenges, and future directions. Prog Aerosp Sci 2019; 111: 100573.

Roshanbin

Abad

Preumont

. Kinematic and aerodynamic enhancement of a robotic hummingbird. AIAA J 2019; 57: 4599–4607.

Timmermans

. Flight performance of nano aerial vehicles with flapping wings: experiment-based system modelling. KU Leuven, Belgium: PhD Dissertation, 2021.

Kumar

Singh

Mohite

, et al. Role of flexibility on the aerodynamic performance of a resonating hummingbird-inspired wing. Aeronaut J 2023; 127: 193–212.

Tobing

Young

Lai

JCS

. Effects of wing flexibility on bumblebee propulsion. J Fluids Struct 2017; 68: 141–157.

10.

Tanaka

Suzuki

Kitamura

, et al. Lift generation of hummingbird wing models with flexible loosened membranes. In: IEEE/RSJ international conference on intelligent robots and systems (IROS), Tokyo, Japan, November 3–7, 2013.

11.

Percin

van Oudheusden

, et al. Wing flexibility effects in clap-and-fling. Int J Micro Air Veh 2011; 3: 217–227.

12.

Zhao

Huang

Deng

, et al. Aerodynamic effects of flexibility in flapping wings. J R Soc Interface 2009; 7: 485–497.

13.

Gehrke

Richeux

Uksul

, et al. Aeroelastic characterisation of a bio-inspired flapping membrane wing. Bioinspir Biomim 2022; 17: 065004.

14.

Shyy

Aono

Chimakurthi

, et al. Recent progress in flapping wing aerodynamics and aeroelasticity. Prog Aerosp Sci 2010; 46: 284–327.

15.

Ingersoll

Lentink

. How the hummingbird wing is tuned for efficient hovering. J Exp Biol 2018; 221: jeb178228.

16.

Maeda

Nakata

Kitamura

, et al. Quantifying the dynamic wing morphing of hovering hummingbird. R Soc Open Sci 2017; 4: 170307.

17.

Walker

Thomas

ALR

Taylor

. Deformable wing kinematics in the desert locust: how and why do camber, twist and topography vary through the stroke? J R Soc Interface 2009; 6: 735–747.

18.

Stanford

Bowman

, et al. Digital image correlation techniques for full-field displacement measurements of micro air vehicle flapping wings. Exp Tech 2009; 33: 53–58.

19.

Doan

Park

, et al. Modal analysis using a virtual speckle pattern based digital image correlation method: an application for an artificial flapping wing. Exp Mech 2022; 62: 253–270.

20.

Choi

Gong

Lee

, et al. Simulation of the flapping wing aerial vehicle using flexible multibody dynamics. Int J Micro Air Veh 2021; 13: 1–12.

21.

Dong

Song

Xue

, et al. 3D numerical simulation of a hovering hummingbird-inspired flapping wing with dynamic morphing. J Appl Fluid Mech 2022; 15: 873–888.