Power Modeling and Experiment Study of Large Flapping-Wing Flying Robot during Forward Flight

A power estimation approach for calculating the power of a flapping-wing air vehicle (FWAV) in forward flight is proposed in this paper. One of the challenges and essential points of FWAVs is endurance. In order to optimize FWAVs, it is necessary to analyze power required for flight in addition to kinematic and aerodynamic analyses of the prototype. Previously, calculating the power of birds was limited to calculating their average power, which assumed the lift was usually the same as the gravitational force. However, the lift varies with the flapping angle during flight. As a result, the power required for forward flight of FWAVs is determined in this work by using a kinematic model of the drive element and wing flapping, along with the aerodynamic model, which varies with the flapping angle during the flapping cycle. Experiments were performed with two prototypes with wingspans of 1.6 and 1.8 m, utilizing a wind tunnel platform. The correlations between power and angle of attack, flapping frequency, and incoming flow velocity were discovered, and recommendations for FWAVs and flying mode design were provided. However, several challenges are highlighted in the application of the model to practical design efforts.


Introduction
Flapping-wing air vehicles (FWAVs) have caught the interest of numerous academics in recent years. Endurance is a main challenge for FWAVs. Compared with birds in nature, the aerodynamic efficiency of FWAVs is not high and its loading capacity is limited, making it difficult to achieve exceptional endurance. Therefore, this paper aims to reveal the power consumption of FWAVs during forward flight by developing a power consumption model, with the goal of optimizing FWAVs. Researchers have built and constructed many prototypes in the field of flapping-wing micro air vehicles (FWMAVs). Typical systems include Robobee with a mass of 60 mg [1] from Harvard University; Hummingbird [2] from Purdue University with a span of 17 cm and a mass of 12.1 g; Nano Hummingbird [3] with a span of 16 cm and a mass of 19 g developed by AeroVironment funded by DARPA; Delfly with a span of 33 cm and a mass of 28.2 g [4] from Delft University of Technology in the Netherlands; and KUBeetle [5] with a span of 16 cm and a mass of 21.4 g from KAIST in Korea. However, due to the small size of FWMAVs, their endurance is severely limited, making outdoor free flying problematic.
In the development of flapping-wing air vehicles (FWAVs) with a large wingspan, Festo designed a two-segment FWAV-SmartBird [6], with a wingspan of 1.96 m and a mass of 480 g. The most important feature of Smartbird is that the wings of Smartbird are able to fold and retract during the upstroke, which reduces the effective area of wing-air interaction, thereby reducing the negative lift generated on the wings and improving flight 1. The fuselage skeleton is made of a high-strength carbon fiber frame with a thickness of 2 mm. It is utilized to mount and connect all structural components and provide structural support for all functional components. As with the skeletal structure of birds, the fuselage skeleton of the HIT-Hawk is constructed with a large number of hollow holes to provide strength while minimizing weight. 2. The wing unit is manufactured entirely of high-strength, lightweight carbon fiber.
The single-segment wing is designed to imitate the size and shape of a large bird's wing, and nylon polyester fabric with exceptional toughness qualities is used to create HIT-Hawk's feathers. Through the reciprocal flapping of the wings, the lift and thrust required for flight can be generated. 3. The tail unit is utilized to maintain the balance of the vehicle and modify its attitude during flight. The tail unit of HIT-Hawk contains two servos, which are used for pitch control and yaw control, and is designed with an outer circular tail shape, mimicking the tail of an eagle. 4. The drive unit of HIT-Hawk consists of a high-speed brushless motor and a reduction gear box with a gear ratio of 81:1. The reduction gear box converts the high speed and low torque motor rotation into a low speed and high torque circular rotation, which is used to drive the reciprocating flapping of the wing unit. 5. Space linkage unit-the linkage unit of HIT Hawk is a four-link RSSR space mechanism coupled to the drive unit and the wing unit. Lift and thrust can be generated by the wing through converting the circular motion of the drive unit into the reciprocating flapping of the wing unit by the space linkage unit. Figure 2 shows the power output and transmission mechanism of the HIT-Hawk. The drive source of the flapping-wing air vehicle is a brushless DC motor, which drives the output shaft to make rotational movement by a gear set with a gear ratio of 81:1. It drives the wing of the FWAV to flap periodically through the space RSSR four-link mechanism.

1.
The fuselage skeleton is made of a high-strength carbon fiber frame with a thickness of 2 mm. It is utilized to mount and connect all structural components and provide structural support for all functional components. As with the skeletal structure of birds, the fuselage skeleton of the HIT-Hawk is constructed with a large number of hollow holes to provide strength while minimizing weight.

2.
The wing unit is manufactured entirely of high-strength, lightweight carbon fiber. The single-segment wing is designed to imitate the size and shape of a large bird's wing, and nylon polyester fabric with exceptional toughness qualities is used to create HIT-Hawk's feathers. Through the reciprocal flapping of the wings, the lift and thrust required for flight can be generated.

3.
The tail unit is utilized to maintain the balance of the vehicle and modify its attitude during flight. The tail unit of HIT-Hawk contains two servos, which are used for pitch control and yaw control, and is designed with an outer circular tail shape, mimicking the tail of an eagle.

4.
The drive unit of HIT-Hawk consists of a high-speed brushless motor and a reduction gear box with a gear ratio of 81:1. The reduction gear box converts the high speed and low torque motor rotation into a low speed and high torque circular rotation, which is used to drive the reciprocating flapping of the wing unit.

5.
Space linkage unit-the linkage unit of HIT Hawk is a four-link RSSR space mechanism coupled to the drive unit and the wing unit. Lift and thrust can be generated by the wing through converting the circular motion of the drive unit into the reciprocating flapping of the wing unit by the space linkage unit. Figure 2 shows the power output and transmission mechanism of the HIT-Hawk. The drive source of the flapping-wing air vehicle is a brushless DC motor, which drives the output shaft to make rotational movement by a gear set with a gear ratio of 81:1. It drives the wing of the FWAV to flap periodically through the space RSSR four-link mechanism.

Kinematics of the Flapping-Wing Air Vehicle Flight
The first and last two spinning pairs and the middle two ball pairs form the spatial RSSR mechanism (see Figure 2).
The kinematic equation, detailed in Appendix A, is obtained to establish the relationship between input angle θ and flapping angle ϕ. The angular displacement of the pendulum CD is described as: where The angular velocity and angular acceleration of the pendulum CD can be calculated by determining the first-order and second-order derivatives from ϕ to time. Here, l 1 = 9 mm, l 2 = 54.5 mm, l 3 = 27 mm and s = 1.97 mm, h = 48.97 mm, and d = 18 mm. The first and last two spinning pairs and the middle two ball pairs form the spatial RSSR mechanism (see Figure 2).
The kinematic equation, detailed in Appendix A, is obtained to establish the relationship between input angle θ and flapping angle  . The angular displacement of the pendulum CD is described as: The flapping motion of birds can be divided into up-and-down flapping, forward and backward swinging, twisting around the spanwise, and folding of the wing surface during upward flapping. During forward flight, forward and backward swinging and folding of the wing are ignored in this work.
The variation law of flapping angle  and twist angle  with time t can be established using this basic geometric model of FWAV: where y is the coordinate of wing in spanwise, f is the flapping frequency, and b is the wingspan of a single wing; max  is the maximum twist angle and 0  is the initial twist angle, which are assumed to be 15° and 5°, respectively, in this paper.
The object of investigation is a tiny strip along the wing span. When the FWAV flies forward horizontally at speed  v , the tiny strip of the wing has three speeds at this time: horizontal forward flight speed  v , flapping speed  v , and twisting speed  v . The The flapping motion of birds can be divided into up-and-down flapping, forward and backward swinging, twisting around the spanwise, and folding of the wing surface during upward flapping. During forward flight, forward and backward swinging and folding of the wing are ignored in this work.
The variation law of flapping angle ϕ and twist angle γ with time t can be established using this basic geometric model of FWAV: where y is the coordinate of wing in spanwise, f is the flapping frequency, and b is the wingspan of a single wing; γ max is the maximum twist angle and γ 0 is the initial twist angle, which are assumed to be 15 • and 5 • , respectively, in this paper. The object of investigation is a tiny strip along the wing span. When the FWAV flies forward horizontally at speed v ∞ , the tiny strip of the wing has three speeds at this time: horizontal forward flight speed v ∞ , flapping speed v ϕ , and twisting speed v γ . The motion model of the wing is shown in Figure 3. If the wing is relatively stationary, the airflow makes relative motion to the wing, and the aerodynamic center of the wing, i.e., the onequarter chord position, is used as the reference point for calculation; the wing motion model is shown in Figure 3.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 5 of 17 motion model of the wing is shown in Figure 3. If the wing is relatively stationary, the airflow makes relative motion to the wing, and the aerodynamic center of the wing, i.e., the one-quarter chord position, is used as the reference point for calculation; the wing motion model is shown in Figure 3. The following formulae can be used to calculate the strip flapping velocity and torsional velocity: The following formulae can be used to calculate the strip flapping velocity and torsional velocity: where c(y) is the chord length of the strip at position y, determined according to the wing shape.
According to the analysis of Rayner [20] and Kuethe [21] on the unsteady aerodynamic mechanism of flying organisms, the downwash velocity can be estimated by the following equation: where α 0 is the zero-lift angle of the wing (assumed to be 0 • here); λ is the aspect ratio of a wing; and θ is the static pitch angle, which is equal to the sum of the wing mounting angle θ 0 and the twist angle γ. Therefore, the relative velocity of the airflow after synthesis at the aerodynamic center of the strip can be further calculated: Because of the combined effect of flapping and torsional motion, the dynamic angle of attack is: Further, the unsteady normal force coefficient method of [22] can be used to correct and obtain the dynamic angle of attack α'. Finally, the total angle of attack of the strip is obtained:

Power Model of the Flapping-Wing Air Vehicle in Forward Flight
Flapping-wing air vehicles, similar to birds, need energy to maintain flight. A FWAV uses its energy to generate lift and overcome drag in flight. The power consumption of the FWAV during horizontal forward flight consists mostly of two components: aerodynamic power and inertial power. The resistance of the air movement used causes the aerodynamic drag that acts on the flapping-wing vehicle. The power required for wing motion during flight is known as inertial power, and it is very low in the event of tiny wing-flapping amplitude and low flapping frequency. During the flight of the HIT-Hawk designed in this study, however, the inertial power is not negligible.

Aerodynamic Power
The aerodynamic drag acting on the wings of a flapping-wing air vehicle during steady-state forward flight can be divided into two categories [13]: Profile drag (D pro ): the drag caused by wing-air friction.
The parasitic drag (D par ), is defined as the drag caused by friction between the frame and the air, and it only affects the airframe. It is included in the total aerodynamic drag. As a result, overall aerodynamic drag may be written as: The total aerodynamic power can then be estimated for the steady-state flight process by multiplying the aerodynamic drag by the forward flight speed. The following equation can be used to calculate the different power components, and the three power components follow: Appl. Sci. 2022, 12, 3176 6 of 16 Different power components can be obtained from Equation (10). The three power components are induced power (P ind ), profile power (P pro ), and parasite power (P par ). The total aerodynamic power required for horizontal flight can be estimated as follows, similar to the drag component: 1.

Induced power
One of the main focuses of this article is the calculation of induced power. The induced power is primarily due to two factors: the tilt of the left and right lift vectors in the course of the flapping cycle, and the continual change in lift vector amplitude over the flapping cycle [12].
The induced power of each wing strip is calculated as the product of the strip's induced drag and velocity.
The induced drag caused by the formation of lift is determined by the wing shape, velocity, and air mass, and is computed as follows: where ρ is the air density, equal to 1.29Kg/m 3 ; V is the relative air velocity; C D,ind is the induced drag coefficient; c(y) is the chord length of the strip; and dy is the width of the strip. The induced drag coefficient is related to the lift coefficient and the wing parameters [23]: where C L is the lift coefficient; k is the induced drag factor, usually taken as 1.1-1.2; and λ is the aspect ratio. The total lift of a segment of the wing can be given by the circulatory normal force and an additional normal force contribution, which comes from added mass effect [24]: where C n (y) = 2π(θ 0 + γ + α ) and The overall lift generated by the flapping-wing vehicle should be more than the force in the vertical direction, which is employed to resist gravity in horizontal level flight due to the presence of the wing-flap angle and angle of attack. This is shown in Figure 4.
where CL is the lift coefficient; k is the induced drag factor, usually taken as 1.1-1.2; and  is the aspect ratio.
The total lift of a segment of the wing can be given by the circulatory normal force and an additional normal force contribution, which comes from added mass effect [24]: where The overall lift generated by the flapping-wing vehicle should be more than the force in the vertical direction, which is employed to resist gravity in horizontal level flight due to the presence of the wing-flap angle and angle of attack. This is shown in Figure 4. As a result, the wing-lift coefficient can be calculated as: In turn, the induced drag coefficient can be derived from the equation above. Thus, the induced power of the wings is equal to the total of the induced powers of each wing  As a result, the wing-lift coefficient can be calculated as: In turn, the induced drag coefficient can be derived from the equation above. Thus, the induced power of the wings is equal to the total of the induced powers of each wing strip, as shown in the equation below: where l is the length of the wing.

Profile power
The profile power component is the most difficult to calculate of all the power components. The reason for this is that the profile drag coefficient C D,pro is dependent on the shape and roughness of the wing surface, and free-stream circumstances and Reynolds number dependence [23]. Reynolds number can be expressed as Re = V·c ν , where c is the average chord length and ν is the kinematic viscosity of air.
To simplify the calculation, the basic expression for the profile drag of a wing strip is: Thus, the profile power can be expressed as: where C D,pro is the profile drag coefficient.

Parasite power
The parasitic power is only associated with the airframe design, not the wing, and is calculated as: where C D,par is the parasitic drag coefficient, S b is the frontal projected area of the body frame, and S b is related to the angle between the body frame and incoming flow.

Inertial Power
The inertial power is the power required to flap the wings of the FWAV, and is zero when the wings remain stationary. Calculating the inertia of the wings around the flapping axis is the most important component for determining the inertia power.
The inertia of a wing contains two parts; one is the inertia I w of the wing itself flapping reciprocally around the flapping axis. As depicted in Figure 5, the wing is decomposed into n units of width a, and the distance of the centerline of the wing unit from the flapping axis is defined as d i , and the mass of each unit is m i .
The inertial power is the power required to flap the wings of the FWAV, and is zero when the wings remain stationary. Calculating the inertia of the wings around the flapping axis is the most important component for determining the inertia power.
The inertia of a wing contains two parts; one is the inertia Iw of the wing itself flapping reciprocally around the flapping axis. As depicted in Figure.5, the wing is decomposed into n units of width a, and the distance of the centerline of the wing unit from the flapping axis is defined as di, and the mass of each unit is mi. Obviously, Iw is obtained by using the parallel axis principle to calculate the inertia of each strip relative to the flapping axis, if the number of units is large.
In addition, the air motion caused by wing flapping must be taken into consideration, which means that the virtual mass of air mv must be attached to the wings. The virtual mass of each segment of the wing can be calculated as: Then, the total virtual mass is described as: As a result, the total inertia of the wings around the flapping axis is: Obviously, I w is obtained by using the parallel axis principle to calculate the inertia of each strip relative to the flapping axis, if the number of units is large.
In addition, the air motion caused by wing flapping must be taken into consideration, which means that the virtual mass of air m v must be attached to the wings. The virtual mass of each segment of the wing can be calculated as: Then, the total virtual mass is described as: As a result, the total inertia of the wings around the flapping axis is: It is worth noting that birds in nature can considerably reduce the wing inertia during upward flapping by folding their wings. However, because this article solely considers a single segment wing of FWAV, the case of folded wings is not addressed here.
The power required when the wings flap at a specific angular velocity equals the rate of change of kinetic energy 0.5Iω 2 against time, hence the inertial power is expressed as: In summary, since the mechanical efficiency of the FWAV is not 100%, then the required power can be calculated as: P tot = (P ind + P pro + P par + P iner )/η (25) where η is the mechanical efficiency of the FWAVs. The mechanical efficiency η is defined as a lumped efficiency, which includes the drive unit transmission efficiency and conversion efficiency. Since the load applied on the wing varies during the flapping cycle, for this brushless DC motor, its mechanical efficiency also varies with the load on the wing during the flapping cycle, but it is difficult for us to obtain its actual value through experiment. Therefore, the lumped efficiency is regarded as a constant value [17,25]. In this paper, the mechanical efficiency η is set at 0.85.

Experiment Setup
This paper analyzes and compares model data with actual flapping data from the wind tunnel experiment, which provides an environment that is steadier than outdoor flight. The experimental setup of the FWAV is shown in Figure 6 and the test platform system in Figure 7. The actual power consumption of the flapping-wing system is calculated by measuring the input current and voltage, and the experimental conditions are changed by controlling the flapping frequency, angle of attack, and wind speed. The test platform system includes the angle of attack adjustment platform, sensors, mounting bracket, HIT-Hawk, and optical platform. The wind tunnel at Harbin Institute of Technology (Shenzhen) was applied to test the HIT-Hawk prototype. The wind tunnel presents as a 24 m × 6 m × 3.6 m test section with a maximum wind speed of 35 m/s. The wind tunnel provides a 1.84 × 10 −5 Pa·s kinematic viscosity of air and the Reynolds number is 5.73 × 10 4~1 .39 × 10 5 .
ciency also varies with the load on the wing during the flapping cycle, but it is difficul for us to obtain its actual value through experiment. Therefore, the lumped efficiency i regarded as a constant value [17,25]. In this paper, the mechanical efficiency  is set a 0.85.

Experiment Setup
This paper analyzes and compares model data with actual flapping data from the wind tunnel experiment, which provides an environment that is steadier than outdoo flight. The experimental setup of the FWAV is shown in Figure 6 and the test platform system in Figure 7. The actual power consumption of the flapping-wing system is calcu lated by measuring the input current and voltage, and the experimental conditions are changed by controlling the flapping frequency, angle of attack, and wind speed. The tes platform system includes the angle of attack adjustment platform, sensors, mounting bracket, HIT-Hawk, and optical platform. The wind tunnel at Harbin Institute of Technol ogy (Shenzhen) was applied to test the HIT-Hawk prototype. The wind tunnel present as a 24 m  6 m  3.6 m test section with a maximum wind speed of 35 m/s. The wind tunnel provides a kinematic viscosity of air and the Reynolds number i . Figure 6. Block diagram of signal acquisition and motion control system principles for the wind tunnel experiment platform. In the experiment, the six-dimensional force/moment, wing flapping angle, current and voltage of the FWAV under incoming flow (5 and 10 m/s) in the wind tunnel wer recorded by LabVIEW to calculate the actual power consumed for the power model veri fication. The M3816BH six-dimensional force sensor from SRI was used for measuring th lift and thrust; the AS5048 module was used as an angle sensor to measure the position o the wings and to calculate the flap frequency; the current sensor and voltage sensor with a resolution of 12 bit were ACS711KLCTR-12AB-T and MCP6L01, respectively. From th In the experiment, the six-dimensional force/moment, wing flapping angle, current, and voltage of the FWAV under incoming flow (5 and 10 m/s) in the wind tunnel were recorded by LabVIEW to calculate the actual power consumed for the power model verification. The M3816BH six-dimensional force sensor from SRI was used for measuring the lift and thrust; the AS5048 module was used as an angle sensor to measure the position of the wings and to calculate the flap frequency; the current sensor and voltage sensor with a resolution of 12 bit were ACS711KLCTR-12AB-T and MCP6L01, respectively. From the product of voltage and current, it is easy to obtain the overall power consumption during the experiment. The NI-6343 acquisition card with highest sampling frequency of 250 kS/s is utilized to acquire current and voltage together with force data. Since the maximum flapping frequency of the designed FWAV is 5 Hz, the sampling frequency used for the experiment is set to 1000 Hz.
Wingspans of 1.8 and 1.6 m were used in the FWAV, see detailed parameters in Table 1. The average chord length of the FWAV was designed to be 0.29 m. The flapping frequency varied in the range 1.5 to 3.5 Hz with a step size of 0.5 Hz. The angles of attack were set to be 0 • , 5 • , 10 • , 15 • , and 20 • . The motion control was achieved by the STM32 platform. To ensure the consistency of the experiments, the wings were tested with the same fuselage skeleton and drive unit. The fuselage skeleton was fastened on a mounting bracket that was mounted on an optical platform. The sensors in the experiment were reset after installing the wings to be tested before turning on the wind tunnel.

Results
First, the forces in the x-axis direction in Figure 8 were measured in the wind tunnel using an experimental prototype with or without wings. The purpose of this experiment was to measure the drags acting only on the wings or only on the fuselage, then to further calculate the parasitic drag and profile drag coefficients. In the wind tunnel, the FWAV with wings was first mounted in the experimental platform and the drag D1 was measured at different wind speeds using the six-dimensional force sensor. Note that the wings were kept fully expanded and static during the experiment. Secondly, the wings were removed and only the fuselage was kept; the drag D2 on the fuselage was measured in the same way. The drag D3 acting only on the wings can be calculated as D3 = D1 − D2. Using Equations (17) and (19), we can calculate cd1 and cd2 obtained from the experiment. The results are shown in Table 2.
Wind tunnel experiments on the wing #1 prototype with a wingspan of 1.8 m and the wing #2 prototype with a wingspan of 1.6 m were done in this study. Except for the wing to be analyzed, all components in the experiments were the same, including the power unit, fuselage skeleton, etc. The initial twist angle γ 0 of a single wing segment was assumed to be 5 • , with a maximum twist angle γ max of 15 • . ured at different wind speeds using the six-dimensional force sensor. Note that the wings were kept fully expanded and static during the experiment. Secondly, the wings were removed and only the fuselage was kept; the drag D2 on the fuselage was measured in the same way. The drag D3 acting only on the wings can be calculated as D3 = D1 − D2. Using Equations (17) and (19), we can calculate cd1 and cd2 obtained from the experiment. The results are shown in Table 2.  Wind tunnel experiments on the wing #1 prototype with a wingspan of 1.8 m and the wing #2 prototype with a wingspan of 1.6 m were done in this study. Except for the wing to be analyzed, all components in the experiments were the same, including the power unit, fuselage skeleton, etc. The initial twist angle 0  of a single wing segment was assumed to be 5°, with a maximum twist angle max  of 15°.
It should be noted that due to the presence of noise, the experimental data were processed with low-pass filtering with a cutoff frequency of 50 Hz. Lift and thrust are very important for the flight of a flapping-wing vehicle. A graph of the lift (vertical direction)

Wing #1
Wing #2  It should be noted that due to the presence of noise, the experimental data were processed with low-pass filtering with a cutoff frequency of 50 Hz. Lift and thrust are very important for the flight of a flapping-wing vehicle. A graph of the lift (vertical direction) and thrust (horizontal direction) measured by experiment in the proposed model is in Figure 8a. As shown in Figure 8a, the lift results measured by model and experiment of wing #1 at angle of attack of 15 • , flapping frequency of 3.5 Hz, and wind speed of 5 m/s are compared. The trend of lift with time calculated by the model is basically the same as the experimental results, but the peak value of the model is slightly larger than the experimental results. The power comparison between the model and the experiment for wing #1 is illustrated in Figure 8b. Because aerodynamic power is cubically related to velocity, the model data can be seen clearly in two peaks in one flapping cycle, which is also corroborated in the experimental data.
To better show the model and experimental results, the average power comparison between the model and experimental results under each condition and the difference percentage between the model and experimental results are plotted.
As shown in Figure 9, with a wind speed of 5 m/s, the experimental data of wing #1 agrees well with the model calculations. At a given wind speed, when the flapping frequency rises, the velocity produced by wing flapping rises as well, and the closer to the wingtip, the greater the v ϕ , which leads to an increase in the sum velocity of the wing relative to the air, thus increasing the power. The angular velocity of the wings, on the other hand, increases as the flapping frequency increases, increasing the inertial power required for wing flapping.
agrees well with the model calculations. At a given wind speed, when the flapping frequency rises, the velocity produced by wing flapping rises as well, and the closer to the wingtip, the greater the v  , which leads to an increase in the sum velocity of the wing relative to the air, thus increasing the power. The angular velocity of the wings, on the other hand, increases as the flapping frequency increases, increasing the inertial power required for wing flapping.  Note that the experimental data are much smaller than the model results at low flapping frequency, i.e., f = 1.5 Hz. This part of the gap could be the fact that in the model, we default the motor to rotate at constant speed. In reality, because of the soft mechanical characteristics of the brushless DC motor, it cannot output enough torque at low speed because the aerodynamic damping torque is variable during the flapping cycle.
Furthermore, the power slightly increases with a greater angle of attack for the same wind speed and flapping frequency. On the one hand, the increase in the angle of attack leads to an increase in the frontal projected area Sb of the fuselage along the incoming flow direction, thus increasing the parasite power. On the other hand, increasing the angle of attack also causes an increase in the lift coefficient C L , so that increases the induced coefficient C D,ind with increasing induced power.
As shown in Figure 10, the results of different wind speeds and flapping frequencies are shown for the wing #1 and wing #2 prototypes at an angle of attack of 5 • .
In Figure 10a,b, the model trends well with the experimental results. However, it can still be seen in Figure 10c,d that at low flapping frequencies, the model results are similar to the experimental results in Figure 9, which are analyzed in detail in the conclusion section. When the flapping frequency is greater than 3 Hz, as the flapping frequency increases, as shown in Figure 10d, the flexible deformation of the wing is somewhat different from the deformation assumed in the paper, and the flexible deformation is difficult to measure in the experiment. Therefore, the difference percentage increases if the flapping frequency is higher than 3.0 Hz. Figure 10, the results of different wind speeds and flapping frequencies are shown for the wing #1 and wing #2 prototypes at an angle of attack of 5°.

As shown in
In Figure 10a,b, the model trends well with the experimental results. However, it can still be seen in Figure 10c,d that at low flapping frequencies, the model results are similar to the experimental results in Figure 9, which are analyzed in detail in the conclusion section. When the flapping frequency is greater than 3 Hz, as the flapping frequency increases, as shown in Figure 10d, the flexible deformation of the wing is somewhat different from the deformation assumed in the paper, and the flexible deformation is difficult to measure in the experiment. Therefore, the difference percentage increases if the flapping frequency is higher than 3.0 Hz.

Discussion
To discuss the feasibility of the model, this section discusses the model and experimental results separately in detail and explains the reasons for the differences between the model and the experimental results.
In Figure 9a,b, from the proposed model in this work, according to Equation (7), the total angle of attack α of the wing strip increases as the static angle of attack θ increases with the same wind speed and flapping frequency in the model, which results in an increase in the lift coefficient C n (y) in Equation (15), where C n (y) is positively associated with the overall angle of attack. This in turn leads to an increase in the induced drag coefficient C D,ind from Equation (14) and the induced power P ind in Equation (16). On the other hand, according to Equation (7), an increase in the static angle of attack θ causes an increase in the induced velocity w, which increases the induced drag power P ind . In addition, the induced power is proportional to the velocity V in Equation (16) and the profile power is proportional to the cube of velocity V 3 in Equation (18), so that the increase in velocity increased by flapping frequency causes an increase in power. Similarly, as the wind speed increases, the parasite power, which is proportional to the cube of horizontal velocity, increases as well, according to Equation (19), in addition to the induced power and profile power. This result is shown in Figure 10.
In the experiments, the average power measured experimentally was obtained by controlling the FWAV to work at different wind speeds, different angles of attack, and different frequencies for the same time with the number of flapping cycles. The trend of the experimental results, together with the values, can be fitted with the model calculations, and increasing any of the conditions of wind speed, angle of attack, or frequency causes an increase in the average power consumption of the FWAV.
However, there are discrepancies between the experimental results and the model when the flapping frequency is 1.5 Hz or the flapping frequency is larger than 3 Hz.
In Figures 9c and 10c, the reason for the large difference between the average power calculated by the model and the average power measured experimentally at very low flapping frequencies comes mainly from the prototype platform. The brushless DC motor (SunnySky X2208) that we used is a motor that is mainly used in quadrotors. The alternating load applied to the motor is constantly changing as the wing flaps periodically, and the soft mechanical characteristics of such motors mean that motor has difficulty with constant speed output when the motor rotates at low speeds. The model assumes a constant speed output of the drive unit at 2πf; the angular velocity and angular acceleration in reality are smaller compared to the model when the motor rotates at low speeds, thus producing an obvious percentage difference at low speeds. Especially for wing #1, which has a larger wingspan and area than wing #2, it is more difficult for the motor to drive the wing #1 flapping at low speed, thus leading to the difference between the average power calculated by the model and the average power measured experimentally at 1.5 Hz flapping frequency, as shown in Figure 9a.
In Figure 10d, the percentage difference is more significant when the flapping frequency is larger than 3.0 Hz. In Equation (2), to simplify the calculation, we simply set the maximum wing twist angle (which occurs at the wing tip) and the initial twist angle to 15 • and 5 • , respectively. However, this is inaccurate because as the flapping frequency changes, so does the flexible deformation of the wing, so that at higher flapping frequencies the assumed twist angle value deviates from the actual twist, resulting in a difference between the model and experimental results.
In this work, the flexible deformation of the wings could not be measured and the stiffness of the wings was difficult to determine, so we took the approach of using assumed values. However, it should be emphasized that, considering the reliability of the drive unit, the FWAV in this paper usually flaps at about 3 Hz or less in level flight to prevent wearing of the gear set.

Conclusions
Unlike previous studies, in this paper, the lift is not directly and simply equal to gravity when calculating the induced drag and induced drag coefficient. Because the lift generated by the wings should change during the flapping cycle, we believe that the induced drag associated with it should also vary. However, such a simulation necessitates a more precise lift calculation, because the induced drag coefficient and induced power can only be determined if the lift calculation is accurate.
However, the limitation of this method is that it relies on the accuracy of the flapping kinematic and aerodynamic models, for example, in the experiments, when the flapping frequency is low and the actual rotational speed deviates from the kinematic model, the calculation results will have deviations. Moreover, the variation of flexible wing deformation with flapping frequency is not considered comprehensively enough in this work, which is a subject we will address, including establishing the complete deformation equations of wing span and axial direction, and verifying them with the assistance of experimental equipment. However, the advantage of this method is to calculate the accurate power consumption of the FWAV during level flight, instead of simply estimating it by gravity. Extended efforts should concentrate on improving the accuracy of aerodynamic models and then using the power model to solve for the optimal flight speed of FWAVs and even to optimize the wing characteristics.
In any case, the proposed model provides some guidance for the prototype's stability and flight strategy design. The prototype's center of gravity is normally situated closer to the rear end of the body to ensure flight stability, so that the prototype is always kept in a "head-up" posture with a large angle of attack. This design, however, is not conducive to expanding the endurance of FWAVs. As a result, during the cruising flight, the angle of attack should be reduced to improve endurance performance. Second, enhancing thrust and wind resistance performance is another technique to extend the range while maintaining enough lift for flight. According to the experimental data, the total power consumption of the FWAVs decreases slightly with increasing forward flight speed within the permitted wind speed. Third, the induced power and inertial power of FWAVs will be much greater if the wings keep flapping than if the wings remain stationary. Thus, the controller should be designed to manipulate the FWAVs to be able to fly similar to birds in nature, with a fully expanded wingspan and capable of performing flap-glide flight. (x C − x B ) 2 + (y C − y B ) 2 + (z C − z B ) 2 = l 2 2 (A1) Therefore, the angular displacement of the pendulum CD is described as: where A = 2l 3 (h − l 1 sin θ) B = −2dl 3 C = l 2 1 − l 2 2 + l 2 3 + h 2 + d 2 + s 2 − 2sl 1 cos θ − 2hl 1 sin θ (A3)