Optimization of lift force for a bio-inspired flapping wing model in hovering flight

Abstract

The lift force is an important component that affects the aerodynamic performance of bio-inspired flapping-wing micro aerial vehicles. However, there is a lack of endeavors in the optimization of the flapping wing parameters that affect the lift force of micro aerial vehicles. This research is therefore to establish a methodology that combines computational fluid dynamic simulation for the evaluation of the lift generation and wing parameters. The Taguchi’s framework for design of experiments, combined with the computational fluid dynamics simulations, is performed on a simplified robotic fly wing model used in the experiment found in Dickinson et al. (1999), under the hovering flight mode, to identify the most influential parameters on the lift generation of micro aerial vehicles. A commercial computational fluid dynamics code, ANSYS/FLUENT, along with a three-dimensional Navier–Stokes solver is used to simulate the unsteady flow field. With the optimization of the time-averaged lift force as the optimization objective, five typical parameters (flapping frequency, maximum angle of attacks during the upstroke and downstroke motions, stroke amplitude, and rotation type) in the flapping trajectory equation are selected as the input factors with each having four levels. Specific computational fluid dynamics cases are simulated in accordance with the chosen orthogonal arrays. By conducting statistical analyses with analyses of means and variance, the flapping frequency and the stroke amplitude are determined to be the two most influential parameters. The response surface of the time-averaged lift force and the power consumption contour are constructed with respect to these two parameters to determine the optimal combination for the generation of lift forces under a specific power constrain and provide a guideline for bio-inspired micro aerial vehicle designs.

Introduction

Our understanding of insect aerodynamics has advanced over the decades, thanks to numerous related research, which, in turn, have inspired and accelerated the development of flapping-wing micro aerial vehicles (MAVs). The interest in technologies that broaden the roles and capabilities of flapping-wing MAVs has rapidly increased over the last decade in the aeronautical engineering community. The important features of MAVs, such as their small size, lightness in weight, low cost, and good manoeuvrability and stealth, mean that they do well in intelligence, surveillance and reconnaissance (ISR) applications in urban and indoor environments.

Various approaches have been adopted in early related research, such as the direct observation of the flapping of real insects and birds, and three-dimensional experimental measurements and numerical simulations on models of flapping insect wings or at the full body scale of insects, to visualize and study the corresponding flow fields, aerodynamics, and flight mechanics. Notably, with the rapid development of high performance computers and computational fluid dynamic (CFD) algorithms over the past few decades, numerical simulation has become an essential method to better understand the aerodynamic performance of the flapping of insect wings. Weis-Fogh ¹ came up with the clap and fling mechanism. Ellington ² indicated that the conventional quasi-steady theory significantly underestimates the lift force in the flapping flights of insects. More recently, Ellington et al. ³ discovered that the leading-edge vortex is one of the important aspects for aerodynamic phenomena after performing flow-visualization of a hawkmoth (Manduca sexta). Later, Liu and Kawachi ⁴ investigated a hovering hawkmoth and numerically validated the flow over this flight mechanism.

In the subsequent experimental observations of Dickinson et al., ⁵ robotic models of double flapping wings were immersed in a tank of mineral oil. To match the real Reynolds number (Re) and reduced frequency in the real flight conditions of the fruit fly (Drosophila), a larger model than the real fly was flapped at a lower frequency in a higher viscosity fluid. They explained that the delayed stall mechanism can enhance the high lift forces during the translation phase while the mechanisms of wake capture and rotational circulation are responsible for the lift force peaks during the rotation phase. Sun and Tang ⁶ conducted a numerical simulation of a model fruit fly wing to confirm the experimental results of Dickinson et al. ⁵ However, they indicated that the generation of high lift force peaks is due to the fast acceleration of the wings during the beginning of a stroke and the fast pitching-up rotation of the wing near the end of a stroke, respectively, rather than the mechanisms of wake capture and rotational circulation. According to Dickinson et al. ⁵ and Sun and Tang, ⁷ the symmetrical rotation of the wings is used for long haul forward flights with less power requirement while advanced and delayed rotations are used for manoeuvring and flexible generation of lift forces. Furthermore, based on the experimental conditions and wing kinematics, Ramamurti and Sandberg ⁸ also performed three-dimensional numerical simulations to confirm and predict the flow structure in different types of rotations induced by the flapping motion of a fruit fly.

In summary, the above research on the flapping flight of various insects has focused on lift generation mechanisms, power expenditure, and the influence of different types of rotations. However, the significance of different flapping kinematic parameters has not been compared and discussed. In terms of engineering applications, especially the bio-inspired design of MAVs, not only is the lift generation mechanism an essential component, but the ways to identify the most influential parameters for the improvement of flight performance are also important.

Several efforts have been made in the optimization of flapping kinematics. Rakotomamonjy et al. ⁹ used both a neural network and a simulation model to maximize the time-averaged lift force. However, the simulation model adopted a typical local two-dimensional approach instead of a CFD-based method, so that the transversal components of the aerodynamic performance of the flow are neglected. Berman and Wang ¹⁰ combined a quasi-two-dimensional aerodynamic model with a hybrid algorithm of a genetic algorithm (GA) and a gradient-based (GB) algorithm to obtain optimal power minimizing kinematics.

On the other hand, Hamdaoui et al. ¹¹ considered a low-fidelity model coupled with a GA to conduct the multi-objective optimization of propulsive efficiency. The performance optimization of forward flights by solving two-dimensional unsteady Navier–Stokes equations along with a GB algorithm is showed in the work of Tuncer and Kaya. ¹² Milano and Gharib ¹³ combined a GA with water towing tank experiments to identify the maximum average lift force. Ghommem et al. ¹⁴ conducted kinematic optimization of a hovering flight to minimize the power consumption by linking two dimensional vortex lattice approach and the global-local optimization methods. The above approaches are restricted by two-dimensional due to low computational efficiency resultant of stochastic analysis. Culbreth et al. ¹⁵ found that there is the issue of high computational cost. They also investigated an optimization method for unconstrained maximum propulsive efficiency by using a three-dimensional Navier–Stokes numerical solver with a GB algorithm, which was conducted on 256 central processing units (CPUs). In order to avoid high computational cost, Jones and Yamaleev ¹⁶ carried out three-dimensional optimization with unsteady viscous flow by using the adjoint-based optimization method with a limited number of grids. Also, considering the high numerical cost encountered in CFD, Walker et al. ¹⁷ addressed the thrust optimization with an analytical aerodynamic theory of a deforming airfoil. Zheng et al. ¹⁸ investigated the optimal stroke amplitude of a hawkmoth with respect to lift force and power consumption by implementing a blade element model in a Navier–Stokes flow solver. The hovering kinematics of both rigid and flexible flapping wings were optimized using a surrogate-based approach. ¹⁹ The aerodynamic model considered no spanwise flow with the two dimensional strip theory manner. In summary, the numerous computations involved in optimization are the main problem to overcome, especially when extended to the three-dimensional flow. This work therefore aims to investigate the optimal kinematic combinations of the fruit fly for high lift forces by combining the CFD and Taguchi methods.

The Taguchi methods, developed by Genichi Taguchi in the 1950s, are statistical methods with the advantage of low cost and high efficiency. This method utilizes the specialized orthogonal to reduce the numbers of experiments and limit the arrangement. ²⁰ The Taguchi methods allow the investigation of how different input factors affect the response mean and variation, and allow estimation of the most influential parameters. ²¹ For the purpose of optimization, an adequate number of different combination cases are needed. By using rigorous and robust orthogonal arrays, the Taguchi methods could therefore help to find the influential parameters and reduce the number of simulation cases without a negative influence on the optimization results. Notably, the CFD simulation for each flapping case generally takes very long time to get accurate outcomes.

Materials and methods

A schematic plot of the methodology is shown in Figure 1, which is composed of a numerical modelling procedure and an optimization process. The CFD numerical procedure is first validated with the experimental data in previous literature. The validated numerical procedure is then used to calculate the lift forces for all cases in a Taguchi orthogonal array, followed by a statistical analysis to determine the most influential parameters and the corresponding response surface. Details on the methodology will be given in the following sections.

Figure 1.

Schematic plot of the research methodology.

Wing model

Wing model and wing kinematics

The wing configuration in this research, shown in Figure 2, is mainly constructed based on the robotic fruit fly in the experiments of Dickinson et al., ⁵ while the wing cross sections have been slightly modified from rectangles to ellipses in order to obtain a smooth flow field. The wing length R is 0.25 m, and the wing area S is 0.0167 m².

Figure 2.

Wing model for numerical simulation.

In order to validate the numerical method, this research uses two baseline cases for comparison purposes with the original experimental data. These two baseline cases only differ in the wing trajectories and follow those in the experiments of Dickinson et al. ⁵ and the numerical study of Sun and Tang. ⁶ Only two motions were considered: the translation motion (flapping motion) and the wing rotation motion (pitching motion), denoted by the stroke angle $\phi (t)$ and the angle of attack $\alpha (t)$ , as shown in Figure 3. The kinematic expressions with respect to time for $\phi (t)$ and $\alpha (t)$ are given by

$\phi (t)={{\phi }_m\sin -1K}_{\phi }\sin -1[K_{\phi }\sin (2\pi \mathrm{ft})]$

(1)

$\alpha (t)={\alpha }_m\tanh K_{\alpha }\tanh [K_{\alpha }\sin (2\pi (\mathrm{ft}+\Delta T))]$

(2)

where $K_{\phi }$ and $K_{\alpha }$ are the effect numbers for shape changes of the kinematics. ²² $K_{\phi }=0.963$ and $K_{\alpha }=1.911$ are selected for this research to accord the wing trajectory in the experiments of Dickinson et al. ⁵ As illustrated in Figure 3E from Dickinson et al., ⁵ three different rotation types, advanced, symmetrical and delayed rotations, are used to illustrate the phase differences between $\phi (t)$ and $\alpha (t)$ , $\Delta T=+0.08$ for the advanced baseline case (Case 1) and $\Delta T=0$ for the symmetrical baseline case (Case 2).

Figure 3.

Illustration of kinematic parameters and coordinate systems of a robotic flapping wing.

The stroke amplitude (or the flapping amplitude) ${\phi }_m$ , defined as the flapping angle from peak-to-peak, is 160° and the flapping frequency f is 0.145 Hz. The angle of attack (AoA) for Case 1 at the mid-stroke position is 40° for both the upstroke and downstroke motions. For Case 2, the AoAs at the mid-stroke position are 20° and 40° during the upstroke and downstroke motions, respectively.

 Scaling parameters

In this research, the wing is hovering in flight. Thus, the reference velocity is defined as the time-averaged velocity of the wingtip during one flapping cycle, as follows

$U_{\mathrm{ref}}=U_{\mathrm{tip}}=2{\phi }_m\mathrm{fR}$

(3)

The Reynolds number (Re) is defined as $\mathrm{Re}\equiv U_{\mathrm{ref}}{L}_{\mathrm{ref}}\nu =U_{\mathrm{ref}}{c}_m\nu$ . Here, the mean chord length c_m (defined as $c_{\text{m}}=S/R=0.0668m$ ) is used as the reference length $L_{\mathrm{ref}}$ , and v is the fluid kinematic viscosity ( $\nu =115\text{cSt}=1.15\times 10-5\text{m}2/\text{s}$ for the mineral oil used in the experiments of Dickinson et al. ⁵ ). Therefore, the Re for the baseline cases is about 117, which corresponds to the real Re of a fruit fly (Re = 136).

The instantaneous lift and drag coefficients are defined as $C_L\equiv L12\rho U_{\mathrm{ref}}^{2}S$ and $C_D\equiv D12\rho U_{\mathrm{ref}}^{2}S$ , respectively. L and D are the instantaneous lift and drag force, respectively, and ρ is the fluid density (for the mineral oil, $\rho =0.88\times 103\text{kg}/\text{m}3$ ).

Numerical methods

Flow solver

The commercial CFD solver ANSYS/FLUENT v14.5 is used to simulate the flow fields around the flapping wing in this research. To model and simulate the flapping wing motion, dynamic mesh algorithms were applied to the problem, which have the capability of solving moving mesh domains in accordance with arbitrary Lagrangian Eulerian (ALE) formulation.

The flow is considered to be laminar and incompressible. Thus, the governing equations are three-dimensional, unsteady and incompressible Navier–Stokes equations, such that

Mass conservation equation

$\nabla •u=0$

(4)

Momentum conservation equation

$\partial \overrightarrow{u}\partial t+(\overrightarrow{u}•\nabla )\overrightarrow{u}=-1\rho \nabla p+\nu \nabla 2\overrightarrow{u}$

(5)

where $\nabla =\partial ..\partial x+\partial ..\partial y+\partial ..\partial z$ is the gradient operator, $u$ the flow velocity, t the time, and p the static flow pressure.

The mass and momentum equations were solved by using finite volume methods approach in the inertial frame to adopt iterative solution strategies. A second-order upwind scheme was used for space discretization. Temporal discretization was limited to only the first-order accuracy due to the use of dynamic mesh algorithms. A standard SIMPLEC algorithm for pressure–velocity coupling was implemented.

Computational domain and boundary conditions

The generation of the computational domain in this simulation follows that in Bos et al. ²³ The spherical domain is divided into an inner and an outer domain, as schematically illustrated in Figure 4. The wing is located in the centre of the spherical domain. The far field boundary is located at a 30 mean chord length away from the wing in order to avoid influence. ²³ The inner domain uses O-type structural grids which rotate together with the flapping wing. Thus, good quality of the mesh is maintained around the wing. The outer domain uses unstructured and fixed meshes. Although the mesh quality may be decreased in comparison to the original generated mesh by the deforming and re-meshing during simulation, the outer domain is so far from the wing that negative effects on the numerical results could be negligible.

Figure 4.

The structured and unstructured meshes.

Two boundary conditions are used in this research: no-slip boundary conditions are applied on the wing body surface and the zero gradient boundary conditions are applied on the far-field boundary. Notably, because this research only focuses on hovering flights, there is no free stream in the simulation.

Optimization method (Taguchi method)

Optimization objective

The optimization objective for this research is to optimize the time-averaged lift force L. When the flapping wings are in hovering flight, the lift force should equal to its gravity weight. Nevertheless, increases in lift force could be beneficial for manoeuvring and flight control^5,7 and also for adding more devices to MAVs, such as micro-cameras, positioning systems, etc.

Identification of independent parameters

The optimization targets for this research are mainly the wing trajectories because their changes are relatively easier in progressing the design of MAVs. From the experience of the MAVs designed by our laboratory, the stroke amplitude and the flapping frequency could be simply changed by varying the frame structure and the battery voltage, respectively. In this research, five typical parameters that are usually used in aerodynamic performance analysis are selected to optimize the wing trajectory. The first four variables are the flapping frequency f, AoA at the midstroke position of the downstroke motion ${\alpha }_{\mathrm{downstroke}}$ , AoA at the midstroke position of the upstroke motion ${\alpha }_{\mathrm{upstroke}}$ , and the flapping amplitude ${\phi }_m$ . The last variable used is different rotation types: advanced, symmetrical, and delayed rotations.

Orthogonal array arrangement

Since the Taguchi method with five input parameters was used, the L16 (45) orthogonal array was chosen. ²⁴ Four levels for each parameter were selected and a total of 16 experiments were conducted in each Taguchi test.

Results and discussion

Validation of numerical method

Mesh resolution study

The numerical simulations of flapping wings have been characterized by low Reynolds numbers and demonstrated significant viscous forces and thick boundary layers. Three different dimensions of meshes are generated by the ICEM, three million for Mesh 1, four million for Mesh 2, and seven million for Mesh 3, respectively. The final comparison of the calculated lift coefficients for the three meshes is shown in Figure 5. The results for Meshes 2 and 3 are almost indistinguishable. The difference between the results for Meshes 1 and 2 is much less than 5%. In summary, a mesh dimension of three million is enough to obtain accurate results for this numerical model. The minimum spacing between the nodes in the boundary layer is 0.006 m.

Figure 5.

Comparison of the instantaneous lift coefficients for three different meshes. $T(=\mathrm{tf})$ is the dimensionless time.

Time step independence validation

As mentioned above, since dynamic mesh algorithms were utilized, FLUENT can be only limited to first order temporal discretization. Thus, proper time step sizes should be chosen to obtain good convergence behaviour and accurately capture the flow physics. Based on the criterion of CFL condition (CFL<1), three time-step sizes were chosen for a flapping frequency of 0.145 Hz for the baseline cases: 0.02 s (345 time steps/cycle), 0.03 s (300 time steps/cycle), and 0.05 s (138 time steps/cycle).

Form Figure 6, there is no apparent difference among the lift coefficient histories for the three different time steps. Therefore, 138 time steps/cycle was chosen for the current study for computational efficiency.

Figure 6.

Comparison of the instantaneous lift coefficients for different time steps.

Note that after the 10^th period, there were almost no changes in aerodynamic force coefficients of the periodic flapping and flow structures. Thus, the numerical results in the 10^th period were considered to be the results of the periodic convergence and chosen for the comparison.

Comparison of numerical results and experimental data

To validate the numerical method, the numerical results for Cases 1 and 2 were compared with the original experimental data of Dickinson et al. ⁵ in one flapping cycle, as shown in Figures 7 and 8, respectively. Notably, the lift forces in Figure 1E and Figure 3B of Dickinson et al. ⁵ were converted into the lift coefficients for comparison.

Figure 7.

Comparison of the numerical and the experimental data in Dickinson et al. ⁵ for Case 1.

Figure 8.

Comparison of the numerical and the experimental data in Dickinson et al. ⁵ for Case 2.

As seen, for both baseline cases, the lift coefficients show reasonable agreement with the experimental values in magnitude and variation trends. Moreover, the present numerical simulations capture the peak lift coefficients well. However, for Case 2, there are relatively large differences at the start of both the upstroke and downstroke motions, when the wing rotation began. These trends in the difference of the magnitude were also shown in the CFD studies by Sun and Tang ⁷ (Figure 4B) and Ramamurti and Sandberg ⁸ (Figure 3). Two possible reasons may be the cause of these differences. The first could be that the experiment used double wing models, while in the numerical simulation only a single wing was used. Thus, the numerical simulation lacks the wing–wing interaction. Secondly, the wing model for the numerical simulations is slightly different from that used in the experiments. The numerical model has an ellipse cross section while the experimental model has a rectangular one, which leads to much smoother flow structures in the numerical simulations than those of the experiments. Nevertheless, the simulated time-averaged lift force is $\overline{L}=0.103$ and is in good agreement with the experimental value of 0.104 for Case 2.

In hovering kinematic models, the time-averaged drag coefficient $\overline{C_D}$ should be zero in one flapping cycle so that the insect could keep hovering without moving forward or backward. Therefore, the simulated time-averaged drag coefficient is −0.009, approximate to zero. The wing motion is proven to be hovering in flight.

Figures 9 and 10 show the vorticity magnitude contours and the pressure contour with velocity streamlines, respectively, at 50% wingspan in the 10^th flapping cycle of Case 2. Compared with Figure 4 in Sun and Tang ⁷ and Figure 7 in Ramamurti and Sandberg, ⁸ the flow structures in this simulation well capture the typical mechanisms of dynamic stall.

Figure 9.

Vorticity magnitude contour at 50% wingspan in one flapping cycle for Case 2.

Figure 10.

Pressure magnitude contour with velocity streamlines at 50% wingspan in one flapping cycle for Case 2.

To further validate the numerical method, the result of Case 2 was compared with a real fruit fly. The data of a real fruit fly were obtained from Weis-Fogh, ¹ Sun and Tang,^6,7 and Shyy et al. ²⁵ According to the definition, the time-averaged lift coefficient of the real fruit fly is 0.374. The lift coefficient of Case 2 is 0.378. This is obviously larger than the real fruit fly, so it could well balance the weight. As explained by Sun and Tang ⁷ and observed by Ellington, ² the higher value of the time-averaged lift coefficient is because the AoA in this numerical study is greater than that of the real fruit fly ( $\text{AoA}\approx 35\circ$ at the midstroke position).

In summary, the numerical method is well validated and considered suitable for the following nonlinear optimization.

Optimization of wing trajectory (first Taguchi test)

Estimate levels of input parameters

Table 1 shows the L16 (45) orthogonal array for the first Taguchi test. As observed, four parameter levels are assigned to all five parameters. For the flapping frequency, the four levels are 0.145, 0.25, 0.402, and 0.5 Hz. The 0.145 and 0.402 Hz frequencies were used in the robotic fly experiments of Dickinson et al. ⁵ The flapping frequency f was extended to 0.5 Hz to match the Reynolds numbers for most hoverflies (Re ≈ 500, Dickinson et al. ⁵ ). The corresponding Re at f = 0.5 Hz is 528. f = 0.25 Hz in between 0.145 and 0.402 Hz was chosen as the fourth level. With regard to the AoA, the experiments of Dickinson et al. ⁵ used 20°, 40°, and 50° during the upstroke motion and 40° and 50° during the downstroke motion. Therefore, 20°, 30°, 40°, and 50° were chosen for the AoA at the upstroke motion, whereas 30°, 40°, 50°, and 60° were selected for the AoA at the downstroke motion to extend the upper and lower bounds used in the experiments. In terms of the stroke amplitude ${\phi }_m$ , Dickinson et al. ⁵ used 138°and 160° in their experiments. Four levels of 120°, 140°, 160°, and 180° were set to moderately extend the lower and upper bounds of the experimental ${\phi }_m$ , accordingly.

Table 1.

L16 (45) orthogonal array arrangement and numerical results of the first Taguchi test.

Case	Rotation type/ $\Delta T$	AoA upstroke (deg)	AoA downstroke (deg)	Stroke amplitude (deg)	Flapping frequency (Hz)	Time-averaged life force (N)
1	Ad/0.08	20	30	120	0.145	0.045
2	Ad/0.08	30	40	140	0.402	0.686
3	Ad/0.08	40	50	160	0.500	1.656
4	Ad/0.08	50	60	180	0.250	0.508
5	Sy/0	20	40	160	0.250	0.232
6	Sy/0	30	30	180	0.500	1.435
7	Sy/0	40	60	120	0.402	0.485
8	Sy/0	50	50	140	0.145	0.082
9	De/-0.08	20	50	180	0.402	0.569
10	De/-0.08	30	60	160	0.145	0.043
11	De/-0.08	40	30	140	0.250	0.118
12	De/-0.08	50	40	120	0.500	0.447
13	Ad/0.1	20	60	140	0.500	0.923
14	Ad/0.1	30	50	120	0.250	0.196
15	Ad/0.1	40	40	180	0.145	0.136
16	Ad/0.1	50	30	160	0.402	0.932

There are three different rotation types and advanced rotation obviously contributes to lift generation during hovering flights.^5–8,23 Hence, besides the three original rotation types in the robotic fly experiments of Dickinson et al. ⁵ ( $\Delta T=+0.08$ advanced-, $\Delta T=0$ symmetrical-, and $\Delta T=-0.08$ delayed rotation types), an advanced rotation of $\Delta T=+0.1$ was also selected for further investigation of its effects on aerodynamic performance.

CFD simulations

The wing model for the Taguchi study and CFD settings were the same as those of baseline cases except for the time step sizes. As discussed above, 138 time steps per flapping cycle is enough for the time step independency in the current numerical model. Thus, four time step sizes, 0.05, 0.033, 0.02, and 0.016 s, were chosen, which corresponded to the four different flapping frequencies of 0.145, 0.25, 0.402, and 0.5 Hz, respectively.

Due to the different stroke amplitudes and flapping frequencies, each experiment in the Taguchi test has its own specific reference velocity and Reynolds number. The calculated reference velocities and Reynolds numbers for all of the experiments in the L16 (45) orthogonal array are listed in Table 2.

Table 2.

Reference velocities and Reynolds numbers for Taguchi cases in Table 1.

Case	1	2	3	4	5	6	7	8
Uref (m/s)	0.159	0.491	0.698	0.393	0.349	0.785	0.735	0.177
Re	102	330	469	264	234	528	283	119
Case	9	10	11	12	13	14	15	16
Uref (m/s)	0.632	0.203	0.305	0.524	0.611	0.262	0.228	0.561
Re	424	136	205	352	410	176	153	377

Time-averaged lift force

The simulated time-averaged lift forces for all cases are listed in Table 1. In Figure 11, it is clearly observed that higher frequency and advanced rotation could yield higher lift forces. This agrees well with the results of previous research.^5–8,23 Furthermore, Taguchi statistical analyses—the analysis of means (ANOM) and the analysis of variance (ANOVA) were conducted.

Figure 11.

Comparisons of the instantaneous lift forces for all 16 cases in the first Taguchi test.

Analysis of means

Taguchi suggested the use of the ANOM to investigate the relative importance of each group of input parameters. The mean is the average response for each combination of input parameter levels. The ANOM lists out the means for each input parameter and calculates the difference between the maximum and minimum mean responses, denoted by delta R. That is, a larger delta value means a greater influence on the optimization objective. According to the orthogonal array, the delta R could be estimated as

$R=max\{\overline{L_i}\}-min\{\overline{L_i}\}$

(6)

Here, $\overline{L_i}$ is the average time-averaged lift force that corresponds to each level.

Figure 12 shows the graphs of the trends of the time-averaged lift forces and Table 3 is the response table for the means. The rank represents the ranking of the effects on the means. As can be observed, the most influential factor is the flapping frequency while other parameters have relatively small effects. The effect trends provide the optimal combination among all 16 cases, which is AoA = 40° for the upstroke motion and AoA = 50° for the downstroke motion, frequency of 0.5 Hz and stroke amplitude of 160° with advanced rotation of $\Delta T=+0.08$ . The means analysis was performed on the basis that other parameters were assumed to have same influence. This leads to its extreme sensitivity to outliers. Thus, the range is a very crude measurement and further analysis is needed for an accurate conclusion.

Figure 12.

Means plot for time-averaged lift forces in the first Taguchi test.

Table 3.

ANOM response of time-averaged lift forces for the first Taguchi test.

Level $\overline{L_i}$	Rotation type	AoA upstroke	AoA downstroke	Stroke amplitude	Frequency
1	0.724	0.442	0.632	0.293	0.076
2	0.547	0.590	0.375	0.452	0.264
3	0.294	0.599	0.626	0.716	0.668
4	0.559	0.492	0.490	0.662	1.115
R	0.429	0.157	0.257	0.423	1.039
Rank	2	5	4	3	1

AOA: angle of attack.

Analysis of variance

By comparing the variance between the group means and the variance within groups, the ANOVA can determine the statistical significance of the input parameters. ²¹ The results of the ANOVA are presented in Table 4.

Table 4.

Results of ANOVA analysis of the averaged-time lift force for the first Taguchi case.

Input parameter	DF	SS	MS	F-Test	p-Value
Rotation type	3	0.377	0.126	0.46	0.714
AoA upstroke	3	0.070	0.023	0.08	0.970
AoA downstroke	3	0.181	0.060	0.21	0.888
Stroke amplitude	3	0.456	0.152	0.57	0.643
Frequency	3	2.553	0.851	9.42	0.002

AoA: angle of attack; DF: degree of freedom; SS: sum of squares; MS: mean squares.

In Table 4, DF is the degree of freedom which equals to the number of input levels minus 1. SS, sum of squares, represents the variation from the mean. MS is mean squares, defined as follows

$\mathrm{MS}=\mathrm{SS}\mathrm{DF}$

(7)

The F-test is used to determine the significance of each input parameter on the response variable, in which F is a ratio of the mean squares within the groups and the mean squares between the groups, denoted by the mean squares of the effect ${\mathrm{MS}}_{\mathrm{effect}}$ and the mean squares of the error ${\mathrm{MS}}_{\mathrm{error}}$ , respectively.^20,21

$F={\mathrm{MS}}_{\mathrm{effect}}{\mathrm{MS}}_{\mathrm{error}}$

(8)

p-Value describes the contribution of each parameter to the influence on the response variable.

The ANOVA is performed at a 5% significance level in this study. The limitation for the F-value is $F_{\mathrm{critical}}=3.49$ , above which the corresponding parameter has significant effects on the time-averaged lift force.

From Table 4, the most significant parameter for the time-averaged lift force is the flapping frequency, which is the same conclusion reached by the means analysis based on the F-value. This conclusion is easy to understand and validated by previous CFD researches. It could also explain why birds and insects have high flapping frequencies. However, we can see that the F-values of other parameters are much less than the limitation due to the significant effect caused by the flapping frequency. In the meantime, the ANOVA shows that the second parameter of a significant effect is the stroke amplitude instead of different rotation type in the means effect analysis. Thus, we conducted another Taguchi test with the same flapping frequency to determine the second significant parameter.

Optimization of wing trajectory (extended Taguchi test)

Estimate levels of input parameters

For the extended Taguchi test, the L16 (45) orthogonal array was used again but the flapping frequency was set to 0.145 Hz and the same for each case to avoid influence. To further examine the relative importance of the remaining four parameters on the time-averaged lift force, the levels of the stroke amplitude and AoA for the upstroke and downstroke motions were set to be the same as those in the first Taguchi test. However, the levels of the different rotation types were changed. Since the analysis from the first Taguchi test and previous studies^5–8,23 show that advanced rotation could generate greater lift forces than symmetrical and delayed rotations, we focused on how advanced rotation could influence the lift force in the extended Taguchi test. Hence, all of the cases used advanced rotation with different starting times. They were $\Delta T=+0.06$ , +0.07, +0.08 and +0.09, respectively.

Time-averaged lift force

The final orthogonal array arrangement and numerical results for the extended Taguchi test are listed in Table 5.

Table 5.

Orthogonal array arrangement and numerical results of the extended Taguchi test.

Case	Rotation type/ $\Delta T$	AoA upstroke (deg)	AoA downstroke (deg)	Stroke amplitude (deg)
Ad/0.06	20	30	120	0.0443
Ad/0.06	30	40	140	0.0390
Ad/0.06	40	50	160	0.1225
Ad/0.06	50	60	180	0.1664
Ad/0.07	20	40	160	0.0811
Ad/0.07	30	30	180	0.1197
Ad/0.07	40	60	120	0.0623
Ad/0.07	50	50	140	0.0982
Ad/0.08	20	50	180	0.0916
Ad/0.08	30	60	160	0.0859
Ad/0.08	40	30	140	0.0811
Ad/0.08	50	40	120	0.0698
Ad/0.09	20	60	140	0.0438
Ad/0.09	30	50	120	0.0496
Ad/0.09	40	40	180	0.1392
Ad/0.09	50	30	160	0.1171

AoA: angle of attack; DF: degree of freedom; SS: sum of squares; MS: mean squares.

Thereafter, a statistical analysis was conducted to the results of the extended Taguchi test. Figure 13 and Table 6 show the means effect plot and the data of the ANOM for the time-averaged lift forces, respectively. Among the four input parameters, the most significantly influential parameter is stroke amplitude.

Figure 13.

Means plot of the time-averaged lift forces for the extended Taguchi test. AoA: angle of attack

Table 6.

ANOM response of time-averaged lift forces for the extended Taguchi test.

Level $\overline{L_i}$	Rotation type	AoA upstroke	AoA downstroke	Stroke amplitude
1	0.093	0.065	0.090	0.057
2	0.090	0.074	0.090	0.066
3	0.082	0.101	0.082	0.102
4	0.087	0.113	0.091	0.129
R	0.011	0.048	0.008	0.073
Rank	3	2	4	1

AoA: angle of attack.

Table 7 summarizes the final results of the ANOVA and F-test. From the ANOVA analysis, the F-value of the stroke amplitude is much larger than the limitation of the F-value ( $F_{\mathrm{critical}}=3.49$ ) while those of the other parameters are less than $F_{\mathrm{critical}}$ . The conclusion is the same as above. Consequently, for the extended Taguchi test, the most important influential parameter among the four parameters is shown to be the stroke amplitude.

Table 7.

Results of ANOVA analysis of the averaged-time lift force for the extended Taguchi test.

Input parameter	DF	SS	MS	F-Test	p-Value
Rotation type	3	0.00026	0.00009	0.05	0.984
AoA upstroke	3	0.00609	0.00203	1.66	0.228
AoA downstroke	3	0.00019	0.00006	0.04	0.990
Stroke amplitude	3	0.01353	0.004511	7.49	0.004

AoA: angle of attack; DF: degree of freedom; SS: sum of squares; MS: mean squares.

Three-dimensional response surface and power consumption contour

From the results of the first and the extended Taguchi tests, the two most effective parameters of the flapping trajectory were determined to be the flapping frequency and the stroke amplitude. Hence, the three-dimensional response surface of the time-averaged lift force was constructed with respect to these two parameters, as shown in Figure 14. Both parameters have four levels which is a total of 16 cases. The different levels for the flapping frequency and the stroke amplitude varied from 0.1 to 0.5 Hz and 120° to 180°, respectively. These levels were selected by summarizing those of previous studies on this robotic flapping wing model.^5–8,23

Figure 14.

Three-dimensional response surface of the time-averaged lift force with respect to the flapping frequency and the stroke amplitude.

Table 8 lists the different levels for the two parameters and the final CFD time-averaged lift force results for these 16 cases. From the response surface, the time-averaged lift force for this flapping wing model increases with an increase in flapping frequency and stroke amplitude. This response surface could be used as a design reference for a bio-inspired MAV. The time-averaged drag coefficients were also calculated for the 16 cases and shown in Table 8. The maximum deviation from 0 drag is less than 5%. The hovering flight condition is thus further verified.

Table 8.

Time-averaged lift forces for different flapping frequencies and stroke amplitudes used to construct the three-dimensional response surface.

	120	140	160	180
0.2	0.065	0.088	0.116	0.148
0.3	0.240	0.287	0.377	0.480
0.4	0.582	0.790	1.033	1.306
0.5	0.920	1.248	1.625	2.051

Implication of the optimization results to the design of a bio-inspired MAV

From the above optimization tests with the Taguchi method, advanced rotation could give relatively higher lift coefficients than symmetrical rotation while delayed rotation has obviously smaller lift coefficients than symmetrical rotation. Nevertheless, the frequency and the stroke amplitude are the most effective parameters for the generation of high lift forces. A strong relationship between the lift force and the flapping frequency has also been validated by both the four-wing and the two-wing MAVs designed by our laboratory. ²⁶ However, high flapping frequency is directly related to the high power consumption. Considering the aerodynamic power, P, only (neglecting the inertial power)

$P=M·\Omega =(M_t·\phi (t)•+M_r·\alpha (t)•)$

(9)

where M is the aerodynamic moment and Ω is the angular velocity. Notably, the translational moment and the rotational moment are combined in P.

Figure 15 shows the contour of the time averaged lift force and the power consumption with respect to the flapping frequency and the stroke amplitude. As expected, the power consumption will increase simultaneously with the lift. Nevertheless, Figure 15 provides a useful reference to MAV designers on the choice of proper power source and a rough estimation of flight endurance. Based on Figures 14 and 15, to provide a sufficient lift force with efficient power consumption for a MAV, the combination of lower flapping frequency and larger stroke amplitude should be considered. Therefore, further experiments will be conducted to validate the importance of stroke amplitude as shown in this study and the response surface.

Figure 15.

The contour of the time-averaged lift force and the power consumption with respect to the flapping frequency and the stroke amplitude.

Conclusion

In this research, three-dimensional CFD simulations of a robotic fly are conducted and validated with experimental data from the literature. The numerically simulated temporal variation of the lift coefficient is in reasonably good agreement with the experimental data. The numerical and experimental time-averaged lift coefficients are very similar. An optimization approach that combines the CFD and the Taguchi methods for the design of bio-inspired MAVs has been successfully established. Through two Taguchi tests, the optimization results are in accordance with previous research and indicate that for time-averaged lift forces, the two most influential parameters are the flapping frequency and the stroke amplitude, respectively. A three-dimensional response surface is then built with respect to the flapping frequency and the stroke amplitude. Consequently, designers of MAVs can focus on either one of the two parameters in to obtain a proper high time-averaged lift force. It is worth noting that this optimization approach could also be applied to other flapping wing MAV models under different flying modes. Multiple targets of optimization, such as high time-averaged lift force, lower power consumption, high lift-drag ratio and high efficiency will be used in studies in the near future.