Optimization of a Twistable Hovering Flapping Wing Inspired by Giant Hummingbirds Using the Unsteady Blade Element Theory

Abstract

Due to the complexity of tailoring the wing flexibility and selecting favorable kinematics, the design of flapping wings is a considerably challenging problem. Therefore, there is an urgent need to investigate methods that can be used to design wings with high energy efficiency. In this study, an optimization model was developed to improve energy efficiency by optimizing wing geometric and kinematic parameters. Then, surrogate optimization was used to solve the design optimization model. Finally, the optimal design parameters and the associated sensitivity were provided. The optimized flapping wing, inspired by hummingbirds, features large geometrical parameters, a moderate amplitude of the flapping angle, and low frequency. With the spanwise twisting deformation considered in the parameterization model, the optimization solver gave an optimized wing with a pitching amplitude of approximately 39 deg at the root and 76 deg at the tip. According to the sensitivity analysis, the length of the wing, flapping frequency, and flapping amplitude are the three critical parameters that determine both force generation and power consumption. The amplitude of the pitching motion at the wing root contributes to lowering power consumption. These results provide some guidance for the optimal design of flapping wings.

1. Introduction

Flapping-wing micro air vehicles (FWMAVs) exploit unsteady aerodynamic mechanisms to enhance their flight performances, attracting the attention of scientists and engineers in bio-inspired engineering [1,2]. The design parameters of the flapping wing as the core component to provide lift mainly include geometric and kinematic parameters. If a mechanical structure is developed to enable controllable wing configuration and additional degrees of freedom are added to the multi-body system, this will increase the demand for onboard power supply, resulting in an increased load on FWMAVs. In contrast, a flexible wing design that combines reciprocating flapping and passive rotation can replicate the primary flight features observed in small insects [3] and hummingbirds [4,5]. In addition, this design provides engineers with a simple approach to developing an energy-efficient flapping wing, as evidenced by previous research [6,7]. However, due to the complexity of tailoring the wing flexibility and selecting favorable kinematics [8], the design of flapping wings is a considerably challenging problem [6,9]. As a result, exploring how to efficiently control wing beat patterns and devise the most power-efficient wings that enable the FWMAV to maintain aloft [10], especially to achieve remarkable hovering flight endurance, has become a new hot topic.

During the past decades, researchers have made many attempts to reveal the influence of wing shape and kinematics on aerodynamic forces and power consumption [11,12,13,14]. The calculation is simplified by treating the flapping wings as idealized rigid plates. With this design, it becomes more feasible to analyze the effects of shape and motion parameters on the performance of FWMAVs. Some studies have explored the impact of wing flexibility on lift and energy requirements by modeling the intricate fluid-structure interaction between the wing and the surrounding air. This approach has enabled researchers to uncover the contribution of both flexibility and inertia to FWMAV performance [15,16,17]. However, to effectively increase the hovering time, finding the best balance between these factors is still an important task from the perspective of developing a FWMAV but remains a largely trial-and-error challenge due to the lack of valid design principles [8]. Therefore, it is important for the design of the FWMAV to determine the design parameters of the flapping wings as soon as possible. In the initial stages of design, some methods can guide the selection of design parameters. For example, we often rely on the scaling law and empirical formula established for flying creatures or the micro vehicles developed previously to determine some parameters [18]. In addition, another method is engineering design through surrogate modeling [19], which has been applied successfully to some specific optimization problems in various engineering fields. Recently, numerical simulations and phenomenology-based aerodynamic modeling [20,21,22,23] have been employed to directly relate the aerodynamic forces exerted on a flapping wing to its time-varying motion and geometry. The effects of the design parameters on the wing performance can be predicted quickly, which allows for combining it with the surrogate optimization method [24] to advance the design of the wing. Using a quasi-steady model instead of the more computationally costly method, such as the direct numerical simulation [15,25,26], allows us to perform optimization programs and detailed sensitivity analysis.

Previously, some attempts were made to find the optimal design of a hovering flapping wing. Using the quasi-steady aerodynamic model, Berman and Wang [27] explored the optimal wing kinematics pattern to minimize the power density of the flapping wing during the hovering flight. They found that this optimization yielded kinematics that are qualitatively and quantitatively similar to the data observed previously on natural insects. Their further research [28,29] concluded that the aerodynamic efficiency of hovering flapping flight with optimal kinematic parameters exceeds that of steady wing flight with the optimal angle of attack. The analytical model accounting for only the translational mechanism was also adopted by Nabawy and Crowther [8] to acquire the minimum power constrained by a given lift, and they concluded that the flapping angle with a triangular profile and the pitching angle with a rectangular profile consume the lowest power with the given wing geometry parameters. X. Ke and W. Zhang [30] adopted a revised quasi-steady aerodynamic model to evaluate the wing performance of dynamically scaled fruit flies. Optimizations of wing geometry and kinematic parameters were conducted, revealing that the optimal flapping angle followed a harmonious profile, while the optimal pitching angle exhibited a rounded trapezoidal profile with a faster time scale of pitching reversal. Although previous studies have explored the optimization design of hovering wings to some extent using various aerodynamic models, the studies mentioned above either focused solely on optimizing kinematic parameters under a given shape or primarily investigated hovering fliers operating at relatively low Reynolds numbers. The wing was modeled as a single rigid body, and its flexibility or morphing capabilities were ignored. Compared to smaller hovering insects, most vertebrates, such as bats and hummingbirds, fly at Reynolds numbers beyond ~10,000. However, flying creatures and robots that operate at higher Reynolds numbers have received less attention in research. Therefore, a broader parameter space, which includes shape and kinematic parameters, still holds significant potential for exploration. Wing performance can be significantly impacted by wing twisting. For example, Hoang Vu Phan et al. [31] conducted a parametric study on a three-dimensional passively twistable wing and reported that the vertical aerodynamic force and the aerodynamic power consumption for optimal flat and twisted wings are almost identical. However, another work [32] showed that the power loading of wings with torsion is higher than flat wings. The function of wing twist in biological hovering flapping flight is not yet fully understood, although twisting has been observed in vertebrate species. Therefore, disregarding spanwise twisting may lead to suboptimal optimization results when accounting for the highly coupled kinematic parameters. Considering that the rotational motion of the flapping wings of insects is not always actively controlled, the wings of FWMAVs are also designed to pitch passively. Some quasi-steady models [33,34,35] were developed previously to predict the passive pitching motion over an entire stroke. Structural dynamics modeling techniques have been utilized to model flexible wings and investigate the impacts of flexibility and inertia on wing performances. The wings were simplified as isotropic or anisotropic plates, reinforced with stiffeners, to study the interaction between passively twisted wings and air. Similarly, comparable studies are observed in the research on wind turbine blades [36,37,38]. Research on the optimization of passively twisted wings can be found in references [39,40].

In this paper, we performed an optimization study to provide a pair of energy-efficient wings for a FWMAV in the initial design stage. Based on the findings of previously reported studies on wing twist effects, we recognized wing-spanwise twisting as a means to improve hovering flight efficiency. We specified a linear distribution of rotation angle along the wingspan on a hummingbird-scale model to introduce dynamic twisting. Accordingly, a set of optimal geometric and kinematic parameters was numerically determined. Next, we performed a sensitivity analysis of design parameters at the optimal point for the established wing model. That helped us identify some key parameters that may affect wing performance. The obtained results will assist in guiding the stiffness design and the deformation constraint of flexible wings widely used today.

The rest of this paper is organized as follows. Section 2.1 develops the parametric characterization of a hummingbird-scale flapping wing with spanwise twisting deformation. Section 2.2 briefly introduces the aerodynamic modeling method used in the current research: unsteady blade element theory. Section 2.3 establishes a mathematical model for the combined optimization and introduces the surrogate optimization used to solve the problem. Section 3.1 gives the optimization results and the aerodynamic and inertia quantities calculated for the optimal wing. Section 3.2 expounds on the sensitivity analysis results, and some discussions are given in Section 3.3. Lastly, the conclusions are drawn at the end of this paper.

2. Materials and Methods

2.1. Morphological and Kinematic Parameterization of a Wing

Hummingbird wings are relatively small and narrow, and the plane shape near the wing root is wide. We think that a quarter ellipse can reflect most features of hummingbird wings, so the planar geometry of the ellipse is established, as shown in Figure 1. More importantly, it also simplifies the process of parameterization of wing geometry. The morphological parameters, such as wing length ($L$) and wing root chord ($C_R$), are annotated there. An offset ($r_0$) exists between the wing root and the flapping axis ($Z-$axis). The chord length ($c(r)$) of any local cross-section at the radial radius of $r$ is assumed to be of the following formula:

$$c(r)=\{\begin{array}{ll}0& (0\le r<r_0)\\ \frac{4\overline{c}}{\pi }\sqrt{1-\frac{{(r-r_0)}^2}{L^2}}& (r_0\le r\le r_0+L)\end{array},$$

(1)

where $\overline{c}$ is the mean chord length and is related to the root chord by $\overline{c}=\pi C_R/4$. The wing area ($S$) is expressed as $S=L\overline{c}$, and the ratio of the wing length ($L$) to the mean chord length ($\overline{c}$) defines the wing aspect ratio ($AR$) as $AR=L/\overline{c}$. The radius of the second moment ($R_2$) of the airfoil about the flapping axis can be calculated according to the following formula:

$$R_2=\sqrt{1/S{\displaystyle {\int }_0^{r_0+L}{r}^{2}c(r)\mathrm{d}r}}.$$

(2)

Consequently, it is necessary to determine three parameters, $L$, $C_R$, and $r_0$, to determine the plane shape of the wing.

Generally, the reciprocating motion of any wing section can be divided into two progressive movements: flapping and pitching motions. These two degrees of freedom without deviation motions of the wing section are considered in the current design. As shown in Figure 2, two orthogonal coordinate systems (the inertial frame $O-XYZ$ and the co-rotating frame $O^{\prime }-\xi \eta \zeta$) and two Euler angles (the flapping angle $\phi$ and the rotation angle ${\theta }_r$) are introduced to describe the spatial location of any wing section with a spanwise distance of $r$ from the flapping axis. With the defined framework in Figure 2, the $\xi$ axis coincides with the straight leading edge at each instant time, and it rotates back and forth along with the wing. The flapping motion is assumed to be a cosinusoidal oscillation. The flapping angle as a function of time is expressed as

$$\phi (t)=-{\phi }_m\cos (2\pi ft)$$

(3)

where ${\phi }_m$ is the flapping amplitude, $f$ is the flapping frequency, and $t$ is time. We abstracted the flapping motion as a cosinusoidal waveform because it has been produced widely in many bionic systems and robot experiments [2,6,7,8]. Furthermore, in engineering, the flapping motion of a wing is usually driven by a linkage mechanism that converts the rotational motion output by the motor into the reciprocating motion of the rocker. A time-varying profile approaching cosinusoidal is usually generated by that. Considering that the leading-edge stiffness is usually large enough to limit the variation of the flapping angle along the wingspan, we make this angle independent of the radial radius. At the same instant, this angle has an identical value in all wing cross-sections. In this case, the angular speed, as expressed in $\dot{\phi }(t)=2\pi f{\phi }_m\sin (2\pi ft),$ can indicate the stroke phase; $\dot{\phi }>0$ and $\dot{\phi }<0$ represent the upstroke and downstroke phases, respectively.

As for the other motion, the pitching axis is the leading edge of this wing. The rotation angle, defined as the angle between the wing chord and the positive direction of the $\eta$-axis, can describe the pitching motion of a two-dimensional wing cross-section. The profile of the rotation angle displays a round-trapezoidal waveform. That was proposed by Berman and Wang [17] and is mathematically represented as

$${\theta }_r(r,t)=\frac{\pi }{2}-\frac{{\theta }_{rm}(r)}{\tanh C_{\theta }}\tanh [C_{\theta }\sin (2\pi ft)],$$

(4)

where ${\theta }_{rm}$ is the amplitude of the pitching angle, and $C_{\theta }$ is a parameter for regulating the duration of stroke reversals. If $C_{\theta }\to 0$, the profile will approach a sinusoidal pattern; if $C_{\theta }\to \infty$, the profile will appear step-shaped. In this work, we limited its value to the range of 0–5. Too-short time intervals are impractical. Different from the invariant flapping angle along the span, the rotation angle ${\theta }_r$ is a function of two variables, time $t$ and radial radius $r$, as seen in Equation (4). So, torsional deformation is introduced due to fact that the amplitude of the pitching angle ${\theta }_{rm}(r)$ is the dependent variable of radial radius $r$. We can also see from Equation (4) that there is no phase difference of the rotation angle relative to the flapping angle for each wing section. At the beginning of each half-stroke, each one is placed vertically at an angle of attack of 90 deg. Here, we assume that the pitching amplitude is distributed linearly along the wingspan, and its value is determined by three variables: the wing length and the amplitude of the pitching angle at the root (${\theta }_{rm}{}_{,r}$) and tip (${\theta }_{rm}{}_{,t}$). Therefore, this distribution can be expressed as

$$\begin{array}{l}{\theta }_{rm}(r)={\theta }_{rm,r}+k(r-r_0)\\ \mathrm{and}k=({\theta }_{rm}{}_{,r}-{\theta }_{rm}{}_{,t})/L,\end{array}$$

(5)

where the slope ($k$) is a quantity independent of time, and the twist angle (${\theta }_{rm}{}_{,r}-{\theta }_{rm}{}_{,t}$) is the difference between ${\theta }_{rm}{}_{,r}$ and ${\theta }_{rm}{}_{,t}$. In order to show that the approximation of the profile and distribution of the rotation angle is reasonable, we modeled the dynamic rotation angles at different radial sections of a deformable wing [41]. Figure 3 shows experimental measurements (shown in dashed lines) and approximate modeling results (drawn in solid lines). Note that the functional form used here is the same as that in Equation (4) except for the start time and the phase difference relative to the flapping angle. As shown in Figure 3, the difference here is that during the plateau stage, there are some oscillations observed in the time history curves obtained through experimental measurements. The wing flexibility results in short phase delays of the rotation angle in various wing sections. However, those effects are hard to assess accurately due to the complex fluid-structure interaction involved, and, therefore, we ignored those in our studies. Linear twist deformations during the flapping of the wings have been achieved in some developed FWMAVs [41,42] and are essential for economical flight based on those studies. Even in the future, the design of flapping wings might also be revolutionized by actively controlling the twist angle of the wings. The morphing pattern of wing twist at the middle of the upstroke is shown in Figure 4. At this specific moment, the angles of attack of the airfoil cross-sections vary along the wingspan.

In summary, a total of four unknown parameters need to be determined for a wing cross-section to define its stroke: $f$ and ${\phi }_m$ in Equation (3), and ${\theta }_{rm}$ and $C_{\theta }$ in Equation (4). Based on the assumption that the rotation angle is distributed linearly along the wingspan, the local value of ${\theta }_{rm}$ depends only on the two parameters: ${\theta }_{rm}{}_{,r}$ and ${\theta }_{rm}{}_{,t}$ in Equation (5). As a result, the parameters that control the wing kinematics become five. Eventually, three morphological parameters and five kinematic parameters constitute the design variables in the optimization design, and those are summarized as [$r_0$, $L$, $C_R$, $f$, ${\phi }_m$, ${\theta }_{rm,r}$, ${\theta }_{rm,t}$, $C_{\theta }$].

2.2. Estimation of Force and Power Consumption

A wing theoretically consists of an infinite number of wing cross-sections, but it can be discretized into finite blade elements for numerical analysis. Since there are distinctive rotation angles in different wing sections even at the same instant, we discretized the wing model into infinitesimal strips along the spanwise directions. An aerodynamic model, named the unsteady blade element theory (UBET) in Refs. [22,43], was then employed to estimate the force generated by any two-dimensional wing section and the associated power requirement. In this model, the resultant aerodynamic load acting on a blade element is decomposed into translation-induced load, rotation-induced load, and added-mass load. In addition, the inertial load can also be obtained with the known mass distribution. The power components required to overcome aerodynamic and inertial moments around the flapping and rotation axes are readily available. Note that in the optimization problem it is difficult to parameterize the mass distribution of a wing in a case where its configuration is unknown. Therefore, in our study, the density was assumed to be homogeneous in order to account for inertial effects. The UBET allows us to quickly evaluate the wing performance with the given morphological and kinematic parameters and to excavate the influence law of parameters, which is extremely important for optimization. Force and power calculations are introduced in the following paragraphs.

For any wing cross-section, the translational force usually depends on its angular stroke velocity $\dot{\phi }$ and geometric angle of attack $\theta$. As shown in Figure 2, $\theta$ can be determined according to the rotation angle ${\theta }_r$. $\theta$ is identical to ${\theta }_r$ during the upstroke, and $\theta$ equals $\pi -{\theta }_r$ during the downstroke. The symbol $V$ stands for inflow velocity, $V_T$ denotes the translational velocity due to the flapping motion, and $V_i$ is the induced velocity. The angle between the inflow velocity and the chord is the effective angle of attack $\alpha$, while $\psi$ is the induced angle of attack, the angle between the inflow velocity and the translational velocity. The equations about the induced velocity and the effective angle of attack can be established based on the blade element theory and the momentum theory [44], as follows:

$$V_T=\left|\dot{\phi }r\right|,$$

(6)

$$\psi ={\tan }^{-1}(\frac{V_i}{V_T}),$$

(7)

$$\alpha =\theta -\psi ,$$

(8)

$$U_{ref}=4{\phi }_{m}f{R}_2,$$

(9)

$$\mathrm{Re}=\frac{{\rho }_{air}{U}_{ref}\overline{c}}{\mu },$$

(10)

$$C_L=(1.966-3.94{\mathrm{Re}}^{-}{}^{0.429})\sin (2\alpha ),$$

(11)

$$C_D={\mathrm{Re}}^{-0.764}+(1.873-3.14{\mathrm{Re}}^{-0.369})(1-\cos (2\alpha )),$$

(12)

$$\mathrm{and}\frac{c}{8{\phi }_{m}r}\sqrt{{(\frac{V_i}{V_T})}^2+1}\left[C_L(\alpha )-C_D(\alpha )\frac{V_i}{V_T}\right]-{(\frac{V_i}{V_T})}^2=0,$$

(13)

where $U_{ref}$ is the reference velocity, and $\mathrm{Re}$ is the Reynolds number. Once the morphological and kinematic parameters ($R_2$, $\overline{c}$, ${\phi }_m$, $f$, and $\dot{\phi }$), the density of air (${\rho }_{air}$), and its kinematic viscosity ($\mu$) are determined, then, for the spanwise blade element with a radial radius of $r$, these equations will be reduced to a set of non-linear equations with two unknowns: $V_i$ and $\alpha$. The corresponding numerical solutions can be obtained by using the Newton–Raphson method. After that, the translational force generated by a wing section can be determined according to

$$\mathrm{d}{F}_{T\eta }(r,t)=-\frac{1}{2}\mathrm{sign}(\dot{\phi }){\rho }_{air}(V_i{}^2+V_T{}^2)[C_L(\alpha )\sin \psi +C_D(\alpha )\cos \psi ]c(r)\mathrm{d}r$$

(14)

$$\mathrm{and}\mathrm{d}{F}_{T\zeta }(r,t)=\frac{1}{2}{\rho }_{air}(V_i{}^2+V_T{}^2)[C_L(\alpha )\cos \psi -C_D(\alpha )\sin \psi ]c(r)\mathrm{d}r,$$

(15)

where $\mathrm{d}{F}_{T\eta }$ and $\mathrm{d}{F}_{T\zeta }$ are the $\eta$- and $\zeta$-axis components of the translational force produced by the blade element with a width of $\mathrm{d}r$, and ‘sign[x]’ is a function of x, defined as

$$\mathrm{sign}[x]=\{\begin{array}{l}1\hspace{1em}\hspace{1em}(x>0)\\ -1\hspace{1em}\hspace{0.17em}(x<0)\end{array}.$$

(16)

The added mass forces stem from a change in the flapping and pitching speeds. In this model, only the air within an elliptical cylinder that is determined by the chord length $c(r)$ as the major diameter, the projected chord length in the motion direction $c(r)\sin \theta$ as the minor diameter, and the thickness $\mathrm{d}r$ as the height are called the ‘added mass’. The mass element $\mathrm{d}{m}_f$ can be calculated as

$$\mathrm{d}{m}_f={\rho }_{air}\frac{\pi }{4}c{(r)}^2\sin \theta \mathrm{d}r.$$

(17)

A point at the half chord $c(r)/2$ is selected as a reference point, and it is considered to be the point of action of the resultant force of the added-mass forces. The normal acceleration $a_n$ of this point is expressed as

$$a_n=r\ddot{\phi }\sin {\theta }_r+\frac{c(r)}{2}({\dot{\phi }}^2\cos {\theta }_r\sin {\theta }_r+{\ddot{\theta }}_r),$$

(18)

where $\ddot{\phi }$ and ${\ddot{\theta }}_r$ are the angular accelerations of flapping and rotation motions, respectively. Therefore, the added-mass force perpendicular to the wing can be decomposed into $\eta$- and $\zeta$-axis components as

$$\mathrm{d}{F}_{\mathrm{A}\eta }(r,t)=-\mathrm{d}{m}_f{a}_n\sin (\pi -{\theta }_r)=-\frac{\pi }{4}{\rho }_{air}c{(r)}^2{a}_n\sin \theta \sin {\theta }_r\mathrm{d}r$$

(19)

$$\mathrm{and}\mathrm{d}{F}_{\mathrm{A}\zeta }(r,t)=-\mathrm{d}{m}_f{a}_n\cos (\pi -{\theta }_r)=\frac{\pi }{4}{\rho }_{air}c{(r)}^2{a}_n\sin \theta \cos {\theta }_r\mathrm{d}r.$$

(20)

The rotational force is related to the pitching motion. The expression for the rotational force acting on a radial cross-section is

$$\mathrm{d}{F}_{\mathrm{R}}=\rho V_T\mathrm{d}{\Gamma }_{\mathrm{R}},$$

(21)

where $\mathrm{d}{\Gamma }_{\mathrm{R}}$ is the circulation induced by the rotation motion and is determined by

$$\mathrm{d}{\Gamma }_{\mathrm{R}}=C_{\mathrm{Rot}}{\dot{\theta }}_{r}c{(r)}^2\mathrm{d}r.$$

(22)

In Equation (22), $C_{\mathrm{Rot}}$ is the rotational force coefficient and can be expressed as

$$C_{\mathrm{Rot}}=\pi [\frac{0.75-x_f}{c(r)}],$$

(23)

where $x_f$ is the distance from the pitching axis to the leading edge. The straight leading edge is the pitching axis in this work, and $x_f$ is thus zero. Eventually, the rotational force can be decomposed into $\eta$- and $\zeta$-axis components as

$$\mathrm{d}{F}_{\mathrm{R}\eta }(r,t)=\mathrm{d}{F}_{\mathrm{R}}\cos (\frac{\pi }{2}+{\theta }_r)=-{\rho }_{air}{V}_T{C}_{\mathrm{Rot}}{\dot{\theta }}_{r}c{(r)}^2\sin {\theta }_r\mathrm{d}r$$

(24)

$$\mathrm{and}\mathrm{d}{F}_{\mathrm{R}\zeta }(r,t)=\mathrm{d}{F}_{\mathrm{R}}\cos {\theta }_r={\rho }_{air}{V}_T{C}_{\mathrm{Rot}}{\dot{\theta }}_{r}c{(r)}^2\cos {\theta }_r\mathrm{d}r.$$

(25)

The inertial force generated by a wing section depends on the chord-wise distributions of the mass and angular acceleration. For a material point with a distance of $x_m$ from the pitching axis, its accelerations in the $\eta$- and $\zeta$-directions can be determined according to

$$\ddot{\eta }=r\ddot{\phi }+x_m({\dot{\phi }}^2\cos {\theta }_r+{\ddot{\theta }}_r\sin {\theta }_r+{\dot{\theta }}_r{}^2\cos {\theta }_r)$$

(26)

$$\mathrm{and}\ddot{\zeta }=-x_m({\ddot{\theta }}_r\cos {\theta }_r-{\dot{\theta }}_r{}^2\sin {\theta }_r).$$

(27)

Thus, the inertial force can be obtained by integrating those components generated by all chord-wise infinitesimal elements from the leading edge (LE) to the trailing edge (TE). Generally, it is easier to obtain the inertial components of a chord-wise element in the two directions. They are calculated by multiplying the mass ($\mathrm{d}{m}_w$) by the acceleration in the corresponding direction. Because the direction of the force is opposite to the direction of the acceleration, they can be expressed as

$$\mathrm{d}{F}_{\mathrm{I}\eta }(r,t)={\displaystyle {\int }_{\mathrm{LE}}^{\mathrm{TE}}-\mathrm{d}{m}_w\ddot{\eta }}={\displaystyle {\int }_{\mathrm{LE}}^{\mathrm{TE}}-{\rho }_S\mathrm{d}{x}_m\mathrm{d}r\ddot{\eta }}={\rho }_S({\displaystyle {\int }_{\mathrm{TE}}^{\mathrm{LE}}\ddot{\eta }\mathrm{d}{x}_m})\mathrm{d}r$$

(28)

$$\mathrm{and}\mathrm{d}{F}_{\mathrm{I}\zeta }(r,t)={\displaystyle {\int }_{\mathrm{LE}}^{\mathrm{TE}}-\mathrm{d}{m}_w\ddot{\zeta }}={\displaystyle {\int }_{\mathrm{LE}}^{\mathrm{TE}}-{\rho }_S\mathrm{d}{x}_m\mathrm{d}r\ddot{\zeta }}={\rho }_S({\displaystyle {\int }_{\mathrm{TE}}^{\mathrm{LE}}\ddot{\zeta }\mathrm{d}{x}_m})\mathrm{d}r,$$

(29)

where the constant ${\rho }_S$ is the surface density of wing mass and $\mathrm{d}{F}_{\mathrm{I}\eta }$ and $\mathrm{d}{F}_{\mathrm{I}\zeta }$ are the inertial force components in the $\eta$- and $\zeta$-directions, respectively.

The power components required by a wing section to overcome the moment around the flapping ($\mathrm{d}{P}_f$) and pitching axes ($\mathrm{d}{P}_r$) at an instant time $t$ can be determined as follows:

$$\mathrm{d}{P}_{f,j}(r,t)=\hspace{0.17em}\dot{\phi }r\cdot \mathrm{d}{F}_{j\eta }(r,t)$$

(30)

$$\mathrm{and}\mathrm{d}{P}_{r,j}(r,t)={\dot{\theta }}_r[x_j\cdot (d{F}_{j\eta }\sin \theta -d{F}_{j\zeta }\cos \theta )],$$

(31)

where the subscript $j$ can take one of ‘$T$’, ‘$A$’, ‘$R$’, and ‘$I$’, $\mathrm{d}{F}_{j\eta }$ and $\mathrm{d}{F}_{j\zeta }$ represent the corresponding force components in the $\eta$- and $\zeta$-directions, $\mathrm{d}{P}_{f,j}$ and $\mathrm{d}{P}_{r,j}$ are the power components associated with the force $\mathrm{d}{F}_{j\eta }$ and $\mathrm{d}{F}_{j\zeta }$, and $x_j$ denotes the distance between the action point of the related force and the pitching axis. The distance from the centers of the translational and rotational forces to the leading edge is determined as a function of the geometric angle of attack, as expressed in

$$x_T(r,\theta )=x_R(r,\theta )=c(r)(\frac{0.82}{\pi }\left|\theta \right|+0.05).$$

(32)

In contrast, the action points of the added-mass and inertial forces are assumed to be the center of the wing chord. The associated expression can be written as

$$x_A(r)=x_I(r)=\frac{c(r)}{2}$$

(33)

A slight difference in the chordwise locations of the action point of resultant forces is estimated to have a weak influence on the total power consumption, as discussed in Ref. [23].

2.3. Establishment of the Optimization Problem

By parameterizing the morphology and kinematics of the wing model, a total of eight decision variables remain to be determined for this optimization: the three parameters related to the geometry of the wing model, $r_0$, $L$, $C_R$, and the five associated kinematic parameters, $f$, ${\phi }_m$, ${\theta }_{rm,r}$, ${\theta }_{rm,t}$, and $C_{\theta }$. Biological scaling law was introduced additionally to limit the range of geometric variables within which the expected wing model is more like a real hummingbird wing in shape. Consequently, some parameters of the hummingbird-inspired flapping wing model are explored only in a limited sample space. Next, we will elaborate on this. According to the relation between body weights and parameters of wing performance [45], the predicted values of wing length and wing area are expressed as follows:

$$\widehat{L}=21.07m^{0.62}$$

(34)

$$\mathrm{and}\widehat{S}=133.42m^{1.15}.$$

(35)

where $m$ is the body mass in grams (g), $\widehat{L}$ is the estimated value of wing length in millimeters (mm), and $\widehat{S}$ is the estimated wing area in square millimeters (mm²). With the assumption that the flat wing has a quarter elliptic shape, the predicted value of the root chord can be simplified as

$${\widehat{C}}_R=\frac{4\widehat{S}}{\pi \widehat{L}}=8.06m^{0.53}.$$

(36)

Given a scaling factor $\sigma$ ($0<\sigma <1$) and the target mass $m$ of the FWMAV to be designed, the wing length and root chord will be confined to a range close to the predicted value obtained in Equations (42) and (43). As a result, we wrote those as

$$L\in \left[\widehat{L}(1-\sigma ),\widehat{L}(1+\sigma )\right]$$

(37)

$$\mathrm{and}{C}_R\in \left[{\widehat{C}}_R(1-\sigma ),{\widehat{C}}_R(1+\sigma )\right].$$

(38)

Under the conditions of $\sigma =0.2$ and $m=30\mathrm{g}$, the constraint boundaries for these two variables are displayed in

Table 1. The decision variables and their boundary values.

Variables	Description	Min	Max
$r_0$	Wing root offset	0.0 mm	100 mm
$L$	Wing length	78.9 mm	188.3 mm
$C_R$	Root chord	39.1 mm	58.7 mm
$f$	Frequency	0.0 Hz	30.0 Hz
${\phi }_m$	Flapping Amplitude	0.0 deg	90.0 deg
${\theta }_{rm,r}$	Pitching Amplitude at the root	0.0 deg	90.0 deg
${\theta }_{rm,t}$	Pitching Amplitude at the wingtip	0.0 deg	90.0 deg
$C_{\theta }$	Regulating parameter of rotation angles	0.0 deg	5.0

. In addition, the root offset $r_0$ is specially set aside for a length of up to 100 mm to meet the needs of the robotic arm design. We can also see from Table 1 that the kinematic parameters are explored in a relatively larger sample space than the morphological parameters.

$$\mathrm{Re}=\frac{4{\rho }_f{\phi }_{m}f{R}_2\overline{c}}{\mu }\le \mathrm{13,000}$$

(39)

Although Table 1 lists the boundaries of $r_0$, $L$, $C_R$, ${\phi }_m$, and $f$, these parameters are bounded additionally by the two conditions. One of them is the applicable condition of the used aerodynamic model. As expressed in Equation (39), the permissible Reynolds number cannot exceed 13,000 [27]. Thus, $r_0$ and $C_R$ are restricted further because $R_2$ and $\overline{c}$ depend on these two parameters in that formula. The other is that the vertical force generated by a couple of flapping wings must be equal to or greater than the total body weight. That means the optimized wing can have enough ability to support the actual MAV weight or go upward against gravity. Given a set of design parameters, including the eight decision variables, the wing performance can be evaluated quickly using the UBET introduced in Section 2.2. Following the earlier discussion, the forces in the $\zeta$-direction produced by a spanwise blade element with the width of $\mathrm{d}r$ at an instant time can be summed up as follows:

$$\mathrm{d}{F}_{\zeta }(r,t)=\mathrm{d}{F}_{T\zeta }(r,t)+\mathrm{d}{F}_{R\zeta }(r,t)+\mathrm{d}{F}_{A\zeta }(r,t)+\mathrm{d}{F}_{I\zeta }(r,t).$$

(40)

The resultant vertical force generated by a wing cross-section over one stroke period can be determined according to

$$\mathrm{d}{F}_{\zeta }(r)=f{\displaystyle {\int }_0^{\frac{1}{f}}\mathrm{d}}{F}_{\zeta }(r,t).$$

(41)

Then, the time-averaged vertical force for one stroke cycle is obtained by integrating the vertical forces across all divided wing sections, as shown in Equation (42):

$$F_{\zeta }={\displaystyle {\int }_0^{r_0+L}\mathrm{d}{F}_{\zeta }}(r)$$

(42)

Eventually, the constraint on the generation of vertical force is written as

$$2F_{\zeta }\ge mg,$$

(43)

where g is the acceleration of gravity and is numerically equal to 9.8 N/kg.

In this study, we focus on the optimal design of energy-minimizing wings in hover. Therefore, the first task is to assess the energy cost of flapping wings with different configurations. Various power components consumed by a wing section, including aerodynamic and inertial power, have been discussed in Section 2.2. The power components required by an entire wing to create flapping and pitching motions can be expressed as

$$P_{f,j}(t)=-{\displaystyle {\int }_0^{r_0+L}\mathrm{d}{P}_{f,j}(r,t)}=-{\displaystyle {\int }_0^{r_0+L}\dot{\phi }r\cdot \mathrm{d}{F}_{j\eta }(r,t)}$$

(44)

$$\mathrm{and}{P}_{r,j}(t)=-{\displaystyle {\int }_0^{r_0+L}\mathrm{d}{P}_{r,j}(r,t)}=-{\displaystyle {\int }_0^{r_0+L}{\dot{\theta }}_r[x_j\cdot (d{F}_{j\eta }\sin \theta -d{F}_{j\zeta }\cos \theta )]},$$

(45)

where $P_{f,T}$($P_{r,T}$), $P_{f,A}$($P_{r,A}$), $P_{f,R}$($P_{r,R}$), and $P_{f,I}$($P_{r,I}$) denote the translational, added mass, rotational, and inertial power components about the flapping (rotation) axis, respectively. However, the mechanical power required to move each wing should be treated with caution. In our calculation, the true cycle-averaged aerodynamic (${\overline{P}}_{m,a}$) and inertial power (${\overline{P}}_{m,i}$) are obtained by averaging only the corresponding positive components, accounting for

$${\overline{P}}_{m,a}=f{\displaystyle {\int }_0^{\frac{1}{f}}\{\Xi [P_{f,T}(t)]+\Xi [P_{f,R}(t)]+\Xi [P_{f,A}(t)]+}\Xi [P_{r,T}(t)]+\Xi [P_{r,R}(t)]+\Xi [P_{r,A}(t)]\}\mathrm{d}t,$$

(46)

$$\mathrm{and}{\overline{P}}_{m,i}=f{\displaystyle {\int }_0^{\frac{1}{f}}\{\Xi [P_{f,I}(t)]+\Xi [P_{r,I}(t)]}\text{}}\mathrm{d}t,$$

(47)

where $P_m$ is the mechanical power, and $\Xi$ is a function, defined as

$$\Xi [x]=\{\begin{array}{l}x(x>0)\\ 0(x\le 0)\end{array}.$$

(48)

The cycle-averaged mechanical power ${\overline{P}}_m$ is the sum of these two powers, expressed as

$${\overline{P}}_m={\overline{P}}_{m,a}+{\overline{P}}_{m,i}.$$

(49)

The strategy of obtaining the energy cost is consistent with that employed in another study [38]. After comparing the cycle-averaged mechanical power obtained using the UBET estimation and the experimental measurement of the robot wing model, the effectiveness of this method was proved. It is also reasonable in our study because the wing system has no elastic storage, and the twisting of the wing model is active instead of relying on resilient energy storage or release. No energy cost is needed to dissipate negative power. This treatment ensures that the negative parts of each power component are excluded from the net cycle-average power:

$$\begin{array}{l}\underset{\chi }{\min }{P}^{*}\\ \mathrm{subject}\mathrm{to}\\ \hspace{1em}\hspace{1em}2F_{\zeta }-mg\ge 0\\ \hspace{1em}\hspace{1em}\mathrm{Re}=\frac{4{\rho }_f{\varphi }_{m}f{R}_2\overline{c}}{\mu }\le \mathrm{13,000}\\ \hspace{1em}\hspace{1em}\chi \in ({\chi }_{LB},{\chi }_{UB})\end{array}$$

(50)

where the vector $\chi$ of design variables includes the eight parameters $\chi$ = [$r_0$, $L$, $C_R$, $f$, ${\phi }_m$, ${\theta }_{rm,r}$, ${\theta }_{rm,t}$, $C_{\theta }$], and ${\chi }_{LB}$ and ${\chi }_{UP}$ are its lower and upper bounds, respectively.

The mechanical power normalized by the MAV mass $P^{*}={\overline{P}}_m/m$ was used as the evaluation index of the energy economy [46,47,48]. Thus, the constrained optimization problem was finally established, as expressed in Equation (51). In order to find the optimal solution of this optimization problem, we converted it into a single-objective unconstrained optimization problem after introducing some penalty terms. Thus, $P^{*}$ becomes the main term of the objective function to be minimized in the current optimization problem. Here, the exterior penalty function method was adopted to build the objective function $F$, as follows:

$$F=P^{*}+r\min {\{(\frac{2F_{\zeta }}{mg}-1),0\}}^2+s\min {\{(\mathrm{13,000}\mu -4{\rho }_f{\phi }_{m}f{R}_2\overline{c}),0\}}^2,$$

(51)

where $\min \{x,0\}$ will return the smaller value between the independent variable $x$ and zero, and the nonlinear inequality constraints about the vertical force $F_{\zeta }$ and the Reynolds number $\mathrm{Re}$ are rewritten in form. In addition, the values of real penalty factors, such as $r$ and $s$, are chosen as 1000. As seen in Equation (50), the corresponding ‘penalty’ terms are added to the objective function for the sample points that violate the constraint conditions. Those terms will return larger values. Thus, the points farther from the feasible domains will be punished more severely. By contrast, the sample points that meet the two constraint conditions will not be punished. As a result, the actual power density can be returned as the objective function value. Eventually, this optimization problem transfers to that displayed as follows:

$$\begin{array}{l}\underset{\chi }{\min }F=P^{*}+r\min \{(\frac{2F_{\zeta }}{mg}-1,0)\}+s\min \{(\mathrm{13,000}\mu -4{\rho }_f{\varphi }_{m}f{R}_2\overline{c}),0\}\\ \mathrm{subject}\mathrm{to}\\ \hspace{1em}\hspace{1em}\chi \in ({\chi }_{LB},{\chi }_{UB})\end{array}$$

(52)

2.4. Optimization Techniques

Surrogate optimization is suitable for those problems with computationally expensive function evaluations. Using the surrogate model instead of the actual simulation model can reduce the computation time, save computational efforts, and preserve acceptable approximation accuracy. The general flowchart of a surrogate optimization program is plotted in Figure 5. The surrogate model algorithms consist of the steps described next. First, an initial experimental design is generated to determine some sample points at which the function evaluation is to be done. Then, the program will spend time obtaining the objective function values at all these points. After that, there is going to be a surrogate model (response surface) built through the data fitting for the initial design points and adopted to predict the objective function values at unsampled points in the variable domain to decide at which point(s) to conduct the subsequent expensive function evaluation. the function evaluation will be conducted at that (those) selected in the previous step. The surrogate model can be updated using the new data set. All of the above steps will be repeated until the user-defined stop criterion is met.

Surrogate model algorithms in the reported literature differ mainly concerning the generation of the initial experimental design, the chosen surrogate model, and the strategy for selecting the sample point in each iteration. In this paper, we used the symmetric Latin hypercube design (‘SLHD’) [49] to generate 30 initial sample points in the design spaces. The return value of the objective function is calculated for these sample points based on this output analysis of the UBET. Next, a surrogate model is built to find a globally good approximation at an unobserved location for this optimization problem. A mixture surrogate model (‘MIX_RcKg’) containing cubic radial basis function interpolation [50], the Kriging model with a Gaussian correlation function (first order regression polynomial) [19], and a multivariate adaptive regression spline [51], are employed to prevent accidentally selecting the worst surrogate model, especially in the case where nothing about the behavior of the black-box objective function is known. For this problem with only continuous variables, we used the candidate point approach (‘CAND’) [52] for selecting the sample point in each iteration. In this approach, two groups of candidates for the next sample point are created. Those in the first group are generated by perturbing the best point found so far. Those in the second group are combined by uniformly selecting sample points from the whole box-constrained variable domain. Then, the objective function values at the candidate points are predicted using the established surrogate model. The distance from each candidate sample point to the set of already sampled points is computed. A weighted sum of these two results is used to determine the best candidate point. That avoids selecting the one that is locally optimal rather than globally optimal. Computational experiments also show that the candidate point approach leads to the best results. In the current problem, the iteration generally starts outside the feasible domain. The sample point that violates the constraint receives the corresponding ‘penalty’ and is excluded from the candidate points. As the number of iterations increases, the increasing penalty forces the iteration point to approach the feasible domain. A surrogate model optimization toolbox [53] and self-developed UBET program were used to solve the optimization problem, and the flowchart of the optimization program is shown in Figure 6.

3. Results and Discussion

3.1. Optimization Results and Optimal Wing Performance

A total of six hundred design points were explored in the sample space for the combined optimization of wing geometric and kinematic parameters. Figure 8a shows the objective function values at these design points on a logarithmic scale. For those points where the obtained value of $F$ is greater than 10, we marked them with hollow red diamonds and identified them as points that violate the constraints and suffer penalties. In contrast, those marked with green circles were recognized to be within feasible regions because the smaller objective function values were returned there. They account for about 52.6% of the design points, more than half of the total. Therefore, it is reasonable to believe that the sample space has been explored well and that the global optimum has been obtained so far.

Table 2. Optimal parameters obtained for the optimization problem.

Variables	Description	Value
$r_0$	Wing root offset	72.27 mm
$L$	Wing length	179.81 mm
$C_R$	Root chord	51.95 mm
$f$	Frequency	6.54 Hz
${\phi }_m$	Flapping Amplitude	65.88 deg
${\theta }_{rm,r}$	Pitching Amplitude at the root	39.01 deg
${\theta }_{rm,t}$	Pitching Amplitude at the wingtip	75.96 deg
$C_{\theta }$	Regulating parameter of rotation angles	3.18

lists the combined optimization results, including the optimal combination of those eight design variables. For optimal wing kinematics, the associated flapping and rotation angles at the root and tip are plotted as a function of normalized time in Figure 8b, where various phases of each half-stroke, such as translation, pronation, and supination, are annotated. The angle of attack gradually increases and decreases during supination and pronation, respectively. It is constant during the translational phase. In other words, this can be indicated by the rate of change of the geometric angle of attack over time. $\dot{\theta }>0$ and $\dot{\theta }<0$ represent supination and pronation, respectively. $\dot{\theta }=0$ represents the translation phase. As shown in this figure, the rotation angle is characterized by a round-trapezoidal profile, with a relatively larger amplitude at the wingtip than at the root. The optimal wing twist angle is approximately 37 deg, and both supination and pronation take a fifth of the time of the half cycle. In terms of wing morphology, all wing geometric parameters are slightly smaller than their respective upper bounds. According to the statistical analysis of the aspect ratio of real hummingbirds [39], $\widehat{AR}=6.9m^{0.09}$, the estimated aspect ratio of a single wing ($\widehat{AR}/2$) was calculated to be 4.68. That is very close to the current result, 4.41. A comparison study for wing models obtained through three different approaches was carried out, and the results were plotted in Figure 7. Compared with the simplified wing from a tiny hummingbird (see Figure 7a), the optimized wing in Figure 7c in this study has a larger area due to a higher loading requirement. The wing planform depicted in Figure 7b,c share similarities, and the former was derived based on a statistical analysis prediction for a specific target load.

Moreover, the Reynolds number that controls the flow around the optimal wing was also examined. Its value is about 12,740, approaching the upper limit of the permissible Reynolds number (13,000). The coupling between geometric and kinematics parameters through Re plays a significant role in the search for the optimal value of each design parameter in the combined optimization. After these analyses, it was found that the optimal wing has the feature of large geometrical parameters, medium amplitude of the flapping angle, and low frequency. In a previous study on a hovering fruit fly’s wing, Xijun Ke et al. [47] used a modified quasi-steady aerodynamic model and hybrid genetic algorithm to obtain a similar result, although the dynamic twisting of the wing was ignored. Different from the current result, the optimal wing used a much higher frequency to produce enough lift at a lower Reynolds number. In a study of a hummingbird-like wing with passive twisting [39], the results from the optimization with morphological and kinematic parameters considered simultaneously also have a similar feature to that found in our study. Due to the lower requirement for the lift there, the optimal single wing length is much shorter, and the flapping amplitude reaches the maximum, up to 90 deg. However, the optimal aspect ratio is close to that obtained in this work. The differences in these results also fully illustrate the existence of a strong coupling between these parameters.

The performance of the optimized wing was also evaluated using the UBET, and all obtained results were plotted in Figure 8c–f. Wherein, Figure 8c shows the time-varying curve of the total vertical force and the respective force components generated by a single wing. It can be seen that inertial force has a significant negative contribution to that in the vertical direction, which is especially obvious when the wing rotates rapidly due to pitching acceleration or deceleration. In addition, during pronation, negative and downward rotational forces are produced by pushing air upward, away from the wing surface. During the whole stroke cycle, other force components promote the generation of positive vertical forces. Among them, the translational force component has the most significant effect. After weighted averaging this force over one stroke period, the total generated by a single wing is approximately 15.05 g. So, the constraint on the vertical force is satisfied. To further quantify the time-varying power, we plotted the power output of the optimal wing along the flapping axis and the pitching axis in Figure 8d and Figure 8e, respectively. These two power outputs include aerodynamic power consisting of various components caused by translational, rotational, and added-mass loads, and inertia power of the wing model itself (displayed in Figure 8). Additionally, the positive values of these power components were summed to obtain the corresponding instantaneous contribution of aerodynamic power to total power. The associated inertial contribution was obtained by truncating the negative portion of the inertial power curve in Figure 8d. These two quantities were then added to calculate the total mechanical power utilized, to estimate the cycle-averaged mechanical power. All components are displayed in Figure 8f. The time-history curves of vertical and horizontal forces are also plotted in Figure 9, and those can be better understood by combining the distribution of ${\theta }_r$ along the span.

For the power around the flapping axis in a stroke cycle (see Figure 8d), the translational component is more dominant than the inertial power, which is positive when the wing accelerates and negative when it decelerates. Drag (the $\eta$-axis force component), always in the opposite direction of the translational velocity, is generated in the forward or backward reciprocating cycle. The translational power is thus positive during the forward and backward strokes to overcome the drag. It peaks in the middle of each half-stroke when the angular speed of the flapping angle is at its maximum. However, the rotational and added-mass components contribute little to flapping power consumption during the supination and pronation of each half-stroke, respectively. Due to the deceleration of pronation and the acceleration of supination, during these two fast rotations, the added-mass forces have different directions (see Figure 8c), always opposite to the normal acceleration of the point used for reference in each spanwise element. The same is true for the rotational force, as it is in the opposite direction to each rotation. These are why the flapping power associated with them has opposite signs during these two rotations. Thus, the added-mass component is dominant during pronation, but the rotational component, by contrast, plays a dominant role in the flapping power consumption during supination. As for the power along the pitching axis in Figure 8e, only the rotational component has positive portions. Because the normal velocity and the rotational force component of a wing section induced by rotation motion are opposite, two positive local peaks occur during fast rotations. However, all other power components have negative values over stroke duration. These results mean that the pitching of this wing does not require much power consumption to overcome the moment around the pitching axis. Compared with the flapping power consumption, it has a relatively smaller magnitude. These are also why subtle differences in the chord-wise location of the centers of these forces have a weak influence on total power consumption, as described earlier. As seen in Figure 8f, the contribution of the aerodynamic component to the mechanical power output has three local peaks similar to those in the curve of the total vertical force (Figure 8g). When the wing decelerates, the negative portion can be dissipated. As a result, the mechanical power out is concentrated mainly during the earlier phase of each half-stroke. The synergistic effect between aerodynamic and inertial contributions leads to three positive peaks in the total mechanical power output. The cycle-averaged power is calculated to be approximately 0.99 W, and the power required to overcome the aerodynamic load is about 1.8 times that of the inertial load. The inertial contribution is not as significant as the aerodynamic contribution, due to the small wing density in this study.

3.2. Sensitivity Analysis of Design Parameters

To figure out how total vertical force and mechanical power change with the design parameters at the optimal point, we evaluated wing performance at a range of design points by perturbing each design parameter while keeping the others constant and equal to their optimized values. It is necessary to understand the sensitivity of design parameters in the early stages of flapping wing design. The lower and upper limits of the design variable are limited to about 0.8 and 1.2 times its optimal value, respectively. Those values of the design parameters are obtained at equal intervals between their boundary values, forming a series of design points. For each design parameter, the results obtained at these points are plotted in Figure 10. All geometric parameters, and those parameters related to the flapping motion, consistently and monotonically increase vertical force and mechanical power. However, the parameters related to the rotation motion reduce them synchronously, except for $C_{\theta }$. These two quantities first decrease and then increase slightly with increasing $C_{\theta }$. When it is improved to a certain extent, the curve in power becomes almost flat. Both of them, by contrast, are minimally affected by $C_{\theta }$. The vertical force and the mechanical power have the feature of simultaneous increase and decrease, and thus the optimal values for these parameters are loosely restricted by the constraint about the vertical force. The opposite effects between these design parameters indicate that the rotation motion cannot be ignored for the optimal design of energy-efficient wings.

For geometric parameters, the single wing length has the most significant influence. Secondly, the significance of wing-root offset and root chord is almost coincident. With a 20% increase relative to the optimum, these three design variables increase the vertical force by approximately 47%, 19%, and 21%, respectively. The corresponding power consumption also increases by about 55%, 23%, and 25%, respectively. Compared with the wing with no offset between the root and the flapping axis, that with a wing-root offset has an additional velocity. It is added to the translational speed that varies linearly along the chord; this effect also leads to similar results with the same trend. Although the increase in the root chord and wing length results in the changes in the wing aspect ratio, the wing area, and the second moment radius, the force and power changes caused by the wing length effect are nonlinear. Those are similar to the changes arising from the parameters related to the flapping motion, such as $f$ and ${\varphi }_m$. When these two parameters are increased by 20% based on their optimal values, the associated force and power are raised by approximately 43% and 73%, and 44% and 61%, respectively. The mechanical power is most affected by the flapping frequency, followed by the amplitude of the flapping angle. However, it should be noted that wing length has the most significant effect on the vertical force among all the parameters of concern. The pitching amplitude in the wing-tip section is one of the most influential parameters related to the rotation motion. The second one is the pitching amplitude at the root of the wing, and the regulating parameter of rotation angle ($C_{\theta }$) has minimal impact on both force and power. When ${\theta }_{rm,t}$, ${\theta }_{rm,r}$, and $C_{\theta }$ are 0.8 times their respective optimal values, the vertical force is approximately 24%, 7%, and 1% higher than that obtained for the optimized wing. The associated power is increased by 38%, 14%, and 4%, respectively. When the pitching amplitude at the tip remains constant, the increase in the wing-root pitching amplitude causes the angle of attack across all wing sections to decrease. The same is true when the wing-tip pitching amplitude decreases while keeping the amplitude at the root unchanged. When the angle of attack decreases, the vertical force generated by the wing reduces. That is the same situation as the change in lift generated by a non-twisting wing at a low angle of attack (<45 deg). How the spanwise distributions of wing rotation angle outside the currently considered range affect wing performance will be discussed further in the next section. The optimal value of $C_{\theta }$ determines the profile of the rotation angle approaching the round-trapezoidal shape. As its value gradually increases, the duration of rotation motions decreases. Due to the increase of the angular acceleration, the negative contribution of the inertial component to the total vertical force becomes greater, and thus the resultant force decreases. However, once the value of $C_{\theta }$ corresponding to the minimum vertical force is exceeded, the contribution of the translational component becomes more dominant. The vertical force gradually increases due to the increased time spent in the translational stage. That also causes a slight increase in the associated power. However, other components, by contrast, are relatively unimportant because the opposite force component is generated during pronation and supination. That makes little contribution to the cycle-averaged value. Further, the power is less sensitive because negative power components are excluded from the mechanical power.

3.3. Effect of Wing Rotation Angle Distribution

The wing-root and wing-tip pitching amplitudes (${\theta }_{rm,r}$ and ${\theta }_{rm,t}$) were selected to analyze their effects on vertical force and mechanical power in the optimized wing and to investigate how spanwise torsion affects wing performance in a broader design space. In Figure 11a, the vertical force contours have an elliptical shape, with higher forces concentrated in smaller ellipses. The force gradually decreases from inside to outside, beginning with the center of the innermost ellipse. However, the contour line of the mechanical power shows an inclined curve. The power consumption is less when both the root and tip sections of the wing model have larger pitching amplitudes, as shown in Figure 11b. When the wing model has no pitching motion and is at an angle of attack of 90 deg, there is almost no vertical force, but only the force component inside the stroke plane is produced. Therefore, a large amount of power is required to overcome the aerodynamic drag. That corresponds to the case in the lower left corner of the contour plots. Compared to that case, the wing model with an angle of attack of 0 deg also produces almost no vertical force. That is because the wing surface is parallel to the airflow direction. Thus, the contribution of the aerodynamic component to the mechanical power is reduced, while only the inertial contribution is retained. The upper-right corner of the contour plots indicates this case. For the mechanical power, those results imply a gradual increase from the lower left corner to the upper right corner of the contour plots. In addition, for a non-twisting wing where ${\theta }_{rm,r}={\theta }_{rm,t}=45\mathrm{deg}$, the maximum vertical force is generated while maintaining moderate power consumption. However, a lower power consumption and a comparable force were obtained for a twisted wing with the combination of a smaller ${\theta }_{rm,r}$ and a larger ${\theta }_{rm,t}$, such as the case where ${\theta }_{rm,r}=20\mathrm{deg}$ and ${\theta }_{rm,t}=60\mathrm{deg}.$ By contrast, a higher power was consumed in the case where ${\theta }_{rm,r}=60\mathrm{deg}$ and ${\theta }_{rm,t}=20\mathrm{deg}$, although the same vertical force is generated. This comparison study shows that an energy-efficient wing has the feature of a larger pitching amplitude at its tip but a lower one at its root. The optimal model obtained in the current optimization has just such characteristics. However, because the requirement for vertical force is not the highest here, the pitching amplitudes at these two spanwise positions vary slightly in value. In addition, the wing-root section of a hummingbird wing is usually at a high angle of attack in nature due to the limitation of the shoulder joint on the rotation motion. That forces them to improve efficiency by lowering the angle of attack in the wing-tip section. This point is qualitatively consistent with the current result.

In a previous comparison study [32], the power loading ($F_{\zeta }/\overline{P}$) of the optimal twisted wing was approximately 5.9% higher than that of the optimal flat wing. The twisted wing with the largest power loading has an average angle of attack of 44 deg at the root and 25 deg at the tip. However, in the current work, the optimized wing has an angle of attack of 50 deg, slightly higher than that obtained in Ref. [32]. This difference comes from the distribution of the rotation angle along the wingspan. A linear distribution was specified on the current wing in our study. Despite that, the average geometric angle of attack for all wing cross-sections is approximately 30 deg in that study, approaching the 32 deg found in the currently obtained wing.

4. Conclusions

A method for finding the optimal combination of wing morphological and kinematic parameters that determine the minimum energy cost during the hovering flight was proposed. The optimization problem was established by parameterizing the geometric shape and three-dimensional kinematics through appropriate simplification. We considered a combination of three-dimensional flapping and pitching motions here. Additionally, the pitching amplitude was assumed to vary linearly along the wingspan to introduce the spanwise twisting deformation pattern to the wing model. That provided an effective way of simplification for the torsional deformation of hummingbird wings and the feasibility of the optimization involving the spanwise twist of the wing model. After parameterization, eight design parameters were obtained to establish a flapping wing configuration. The wing model was discretized into finite blade elements along the wingspan, which were analyzed using the unsteady blade element theory. It allowed for rapid evaluation of wing performance across various combinations of design parameters. Considering the design goal of a high-efficiency flapping wing, we transformed this engineering design problem into a mathematical optimization problem. In order to minimize energy costs, we formulated an optimization problem for the wing’s geometry and kinematics parameters based on mechanical power consumption. An optimization algorithm based on surrogate modeling was adopted to implement the combined optimization of design parameters. Ultimately, we identified the optimal combination of parameters that satisfied the constraint for vertical force. The current study proposed a design method for passively twisted wings, which can provide an initial design of hover-capable wings for FWMAVs according to the desired target load.

To assess the impact of individual design parameters on wing performance, we conducted a sensitivity analysis by traversing one parameter near the optimum point while keeping the other parameter values constant. The results indicate that wing length has the most significant influence on the vertical force, followed by frequency and flapping amplitude, which have similar effects and are ranked second and third, respectively. The remaining parameters showed a relatively lower sensitivity. Mechanical power is most significantly influenced by frequency, followed by flapping amplitude and wing length. It was also found that the root bending angle has a significant impact. A decrease in this parameter improves the efficiency of vertical force generation. These results provided a basis for adjusting the wing design parameters to improve hover efficiency and are significant in optimizing wing performance.

Finally, we further investigated the effect of wing torsion. Our findings suggest that an energy-efficient wing exhibits a larger pitching amplitude at the wingtip and a lower amplitude at the root. The optimized wing obtained in this work has just such characteristics. That confirms the effectiveness of the optimization results and the role of a twisted wing in reducing the energy cost of hovering. It is constructive for adjusting mass and stiffness distribution to achieve an optimal pitch amplitude distribution.