Next Article in Journal
Strontium and Copper Co-Doped Multifunctional Calcium Phosphates: Biomimetic and Antibacterial Materials for Bone Implants
Previous Article in Journal
Morphological Reconstruction for Variable Wing Leading Edge Based on the Node Curvature Vectors
Previous Article in Special Issue
A Photosensitivity-Enhanced Plant Growth Algorithm for UAV Path Planning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Simulation of the Advantages of the Figure-Eight Flapping Motion of an Insect on Aerodynamics under Low Reynolds Number Conditions

1
Department of Master’s Program of Mechanophysics, Kyoto Institute of Technology Matsugasaki Goshokaido-cho, Sakyo-ku, Kyoto 606-8585, Japan
2
Department of Mechanical Engineering, Kyoto Institute of Technology Matsugasaki Goshokaido-cho, Sakyo-ku, Kyoto 606-8585, Japan
*
Author to whom correspondence should be addressed.
Biomimetics 2024, 9(4), 249; https://doi.org/10.3390/biomimetics9040249
Submission received: 25 February 2024 / Revised: 17 April 2024 / Accepted: 18 April 2024 / Published: 20 April 2024
(This article belongs to the Special Issue Bio-Inspired Design and Control of Unmanned Aerial Vehicles (UAVs))

Abstract

:
In proceeding with the advanced development of small unmanned aerial vehicles (UAVs), which are small flying machines, understanding the flight of insects is important because UAVs that use flight are attracting attention. The figure-eight trajectory of the wing tips is often observed in the flight of insects. In this study, we investigated the more efficient figure-eight motion patterns in generating lift during the hovering motion and the relationship between figure-eight motion and Reynolds number. For this purpose, we compared the ratios of the cycle-averaged lift coefficient to the power coefficient generated from each motion by varying the elevation motion angle, which is the rotational motion that represents the figure-eight motion, and the Reynolds number. The result showed that the motion with a smaller initial phase of the elevation motion angle ( φ e 0 90 ° ) could generate lift more efficiently at all Reynolds numbers. In addition, the figure-eight motion was more effective when the Reynolds number was low.

1. Introduction

Research on small unmanned aerial vehicles (UAVs) is attracting attention in military and civilian applications [1]. There are three main types of small UAVs: micro air vehicles (MAVs), nano air vehicles (NAVs), and pico air vehicles (PAVs). The size of each UAV is defined as follows. For MAVs, the Defense Advanced Research Projects Agency (DARPA) defined UAVs as having a size within 15 cm, a mass of 50–100 g, a flight speed of 30–60 km/h, and a flight duration of 10 km as the goal of developing UAVs by Broad Agency Announcement (BAA 97-29) in 1997 [2]. In addition, for NAVs, the DARPA defined UAVs as having a size within 7.5 cm, a mass of 10 g, a flight speed of 18–36 km/h, and a flight duration of 1 km as the goal of developing UAVs by Broad Agency Announcement (BAA 06-06) in 2005 [3]. Furthermore, for PAVs, Wood et al. [4] defined them as UAVs with a size of within 5 cm and a mass of 500 mg. The performance of the small UAVs makes it possible for UAVs to enter spaces where people cannot enter and tight spaces in the case of a natural disaster, accident, or war, and to easily check the conditions of these areas using cameras and sensors. The aforementioned application potential is the development purpose of UAVs. In achieving this purpose, small UAVs are required not only to downsize but also to have high flight performance such as the ability to fly freely, high hovering power, and high flight stability to withstand wind in outdoor environments. Therefore, considerable research and development have been conducted to achieve downsizing and motion control [5,6,7]. These types of UAVs have been developed for research purposes: fixed-wing, rotary-wing, and flapping-wing. Among these, various studies on fixed-wing and rotary-wing UAVs have been conducted because they are similar to airplanes and helicopters and are easy to design. Fixed-wing UAVs can efficiently fly long distances; however, they are difficult to hover. On the contrary, rotary-wing UAVs have high hovering power and maneuverability, making them suitable for small UAV development [6]. However, fixed-wing UAVs cannot efficiently fly as the wing size is increased until they generate more lift [8], and rotary-wing UAVs do not show high flight performance because of their high power consumption [9] if they are downsized. On the contrary, flapping-wing UAVs can fly more efficiently than rotary-wing UAVs under low Reynolds number conditions [10], and these UAVs have been attracting attention recently. Flapping-wing UAVs were developed based on the flapping motion of birds and insects, which can not only generate enough lift to support their own weight but also have high flight performance, exhibited in such actions as making sharp turns and taking off swiftly in nature [11]. However, the development of flapping-wing UAVs remains in its nascent stages [6], and to advance the development of advanced small UAVs, deepening our understanding of the characteristics of the fluid surrounding insects and their flight is necessary, which have the same Reynolds number as small UAVs.
Insects have a high flight performance because of the high frequency and unsteady motion of their wings. However, the motion is completely different from the mechanism of lift generation by steady flow in airplanes and helicopters, and much has not been clarified because such a mechanism is difficult to explain using conventional aerodynamic principles. Therefore, considerable research on insect flight has been conducted to clarify the mechanism of lift generation in insects. Ellington et al. [12] indicated that a leading-edge vortex is a vortex that is generated from the leading edge of the wings, and more lift is generated by attaching this vortex to the wings. This factor primarily explains the mechanism of lift generation in insects, and it is common to all insects, including small and large insects [13,14]. In addition, regarding the flapping motion, Weis-Fogh [15] reported that clap and fling are the motions of opening and closing the wings, which generate more lift. Moreover, Cheng et al. [16] indicated that smaller insects perform a U-shaped upstroke motion, which generates more lift. Thus, we can understand the mechanism of lift generation in insects by investigating their flapping motion.
In the flapping motion of insects, simple motions such as line-shaped, oval-shaped, and the figure-eight trajectory of the wing tips (i.e., figure-eight motion) are often observed [17,18]. Considerable research on these motions has been conducted, and Aghav [19] and Lehmann et al. [20] indicated that figure-eight motion generated the most lift in these motions. Therefore, understanding the figure-eight motion will provide new knowledge for the development of UAVs that refer to the flapping motion of insects. Regarding the figure-8 motion, Galinski et al. [21] achieved an electromechanical device incorporating the figure-8 motion via experimentation. In addition, Ishihara [22] indicated that the figure-8 motion can be created via fluid–structure interaction (FSI) and that the lift increased significantly as the wing tip path mode shifted to the figure-eight mode. Moreover, many studies [20,23,24] have compared different strokes of insect wings, such as line-shaped, oval-shaped, figure-eight-shaped, and pear-shaped, and reported that the differences in their motions have a significant effect on the generation of lift. Therefore, we focused on the differences in the motion mode of the figure-eight motions in particular with the self-intersection of the wing-tip trajectories in detail. In addition, the motion is observed in insects of various sizes; however, it is especially observed in insects with low Reynolds numbers ( R e 130 ) [17,25]. Moreover, insects with high Reynolds numbers ( R e 720 ) move their wings horizontally (long and thin figure-eights) relative to the flapping plane, whereas insects with low Reynolds numbers move their wings vertically (thick figure-eights) with a large stroke deviation [17]. However, few studies have summarized the relationship between figure-eight motion and Reynolds number in detail. Therefore, this study aims to investigate the more efficient figure-eight motion patterns in generating lift during hovering motion related to the Reynolds number.

2. Materials and Methods

2.1. Model of an Insect

In this study, fruit fry was used as a model for analyzing the flapping of small insects. The parameters of the insect were based on the physical properties of fruit fry used by Aono et al. [26], and these values are shown in Table 1. The flight mode adopted hovering, which is the simplest flapping flight mode. In a previous study, Yamauchi et al. [27] reported that the body’s presence was not important in hovering because the presence has little effect on the cycle-averaged lift. Therefore, the analysis was performed with only two wings, without the reproduction of the insect body (Figure 1). Regarding the insect wing, it may be possible to perform a simplified analysis with only one wing using symmetric boundary conditions. However, in this study, it is necessary to capture the fluid dynamic interactions generated by each wing to consider the effects of the different motion modes of the wings. Therefore, the analysis was performed with two wings for higher accuracy of the results. The wings were reproduced with rectangular rigid plates without thickness to capture the shape of the wings on Cartesian grids as easily as possible. S W shown in Figure 1 is the wing area.
Moreover, the Reynolds number used in this study was defined as follows:
R e = u t i p ¯ C m ν a i r .
Further, the analysis in Section 3.1 and Section 3.3 was performed with R e 134 based on the values shown in Table 1. In addition, the mean velocity of the wing tip ( u t i p ¯ = 0.645 , 1.29 , 2.58 , 5.16 , and 10.32 ) was varied and analyzed as R e 33.5 , 67 , 134 , 268 , and 536 when the Reynolds number dependence was evaluated in Section 3.4. Note that the frequency of the wings ( f = 54.5 , 109 , 218 , 436 , and 872 ) was similarly varied when the mean velocity of the wing tip ( u t i p ¯ = 0.645 , 1.29 , 2.58 , 5.16 , and 10.32 ) was varied.

2.2. Motion of the Insect

As shown in Figure 2, the flapping motion of an insect is represented by combining the rotational motion of the three axes. Each rotation angle is represented using Equations (2)–(4), where the positional angle, feathering angle, and elevation angle are around the z -axis, wing axis x W , and y -axis, respectively. Note that the lapping motion of an insect was based on the study of fruit fry used by Aono et al. [26] as described in Section 2.1. In this study, a figure-eight motion was represented by varying the elevation motion angle θ e , which is a rotational motion that represents the figure-eight motion. The motions observed by the change of this angle were motions with the self-intersection of the wing-tip trajectory in the other motions ( 0 < φ e 0 < 180 and 180 < φ e 0 < 360), except for the three line-shaped ( θ e = 0 ° ) and U-shaped ( φ e 0 = 0 ° and φ e 0 = 180 ° ) motions. Therefore, in summarizing all the motions, we defined the motion with θ e = 0 ° as a motion “without figure-eight motion” and the motion with θ e 0 ° as a motion “with figure-eight motion”.
θ p = θ p , a m p cos 2 π f t + φ p ,
θ f = θ f , a m p sin 2 π f t + φ f ,
θ e = θ e , a m p cos 2 π 2 f t φ e 0 + φ e ,
where θ p is the flapping angle, θ f is the feathering angle, and θ e is the elevation angle; θ p , a m p , θ f , a m p , and θ e , a m p are each amplitude; φ p , φ f , and φ e are each initial position; φ e 0 is the initial phase of the elevation angle and f is the flapping frequency. In this study, θ p , a m p , θ f , a m p , θ e , a m p , φ p , φ f , φ e , φ e 0 , and f were set to 70 ° , 70 ° , 10 ° , 10 ° , 0 ° , 10 ° , 10 ° , and 218   H z , respectively. Note that these parameters were also based on the properties of fruit fry used by Aono et al. [26] as described in Section 2.1. Based on the abovementioned conditions, Figure 3 shows the time history of the represented rotational motion angles θ p , θ f , and θ e , and the trajectories of the wing tips.

2.3. Governing Equation

In this study, the three-dimensional normalized lattice Boltzmann method [28,29] (3D27V model) with incompressible formulation was used as the governing equation for fluid resulting from three-dimensional flapping motion. The method was developed on the basis of the lattice Boltzmann method [30,31], which has the advantages of a simple algorithm and high computational efficiency. Thus, the proposed method can reduce memory usage and improve computational stability while maintaining the advantages of the lattice Boltzmann method.
The distribution function f α in the lattice Boltzmann equation is expressed using a discrete velocity vector e α . In the normalized lattice Boltzmann equation, the incompressible Navier–Stokes equation can be represented using the second-order moments of the distribution function as follows:
f α ω α a 0 + b i e α i + c i j e α i e α j ,
where ω α is the weight coefficient and e α is the velocity vector, which are given in Table 2 for the advection direction α of the 3D27V model shown in Figure 4. Note that the 3D27V model was used in this study instead of the simple 3D19V model to capture more accurately complex flow fields in three dimensions.
In addition, a 0 , b i , and c i j in Equation (5) are constants that satisfy the following relations:
ρ = α f α ,
ρ u i = α e α i f α ,
Π i j n e q = α e α i e α j f α c 2 3 ρ δ i j ρ u i u j ,
where ρ is the fluid density, ρ u i is the moment, and Π i j n e q is the nonequilibrium part of the stress tensor. Thus, by substituting Equations (6)–(8) into Equation (5), the distribution function f α can be expressed as follows:
f α = ω α ρ 1 + 3 e α i u i c 2 + 9 e α i u i 2 2 c 4 3 u i u i 2 c 2 + 9 ω α 2 c 2 e α i e α j c 2 1 3 δ i j Π i j n e q ,
where the first term on the right is the equilibrium distribution function f α e q and the second term on the right is the nonequilibrium part of the distribution function f α n e q . Thus, they can be replaced as follows:
f α = f α e q + f α n e q ,
f α e q = ω α ρ 1 + 3 e α i u i c 2 + 9 e α i u i 2 2 c 4 3 u i u i 2 c 2 ,
f α n e q = 9 ω α 2 c 2 e α i e α j c 2 1 3 δ i j Π i j n e q .
Further, the time evolution equation of the lattice Boltzmann equation can be expressed as follows:
f α t + δ t , x + e α δ t = f α t ,   x + 1 1 τ f α n e q t , x ,
where τ is the relaxation time, which is defined as follows:
τ = 3 ν c δ x + 1 2 ,
where ν is the kinematic viscosity.
As previously stated, the incompressibility formulation is applied to reduce the error caused by compressibility. In this case, the pressure distribution function p α is defined using the following density distribution function f α :
p α = c s 2 f α ,
where c s is the speed of sound, which is defined as follows:
c s = c 3 .
Thus, the pressure p , velocity component u i , and nonequilibrium part of the stress tensor Π i j n e q can be expressed as follows:
p = α p α ,
u i = 1 ρ c s 2 α e α i p α ,
Π i j n e q = 1 c s 2 α e α i e α j p α c 2 3 ρ δ i j c 2 3 ρ u i u j .
Furthermore, the equilibrium distribution function p α e q can be expressed as follows:
p α e q = ω α p + ρ 0 e α · u + 3 e α · u 2 2 c 2 u 2 2 .
Thus, the time evolution equation of the lattice Boltzmann equation with the incompressibility formulation is defined as follows:
p α t + δ t , x + e α δ t = p α e q t ,   x + 1 1 τ p α n e q t , x .

2.4. Virtual Flux Method

In this study, the virtual flux method (VFM) [32,33,34] was used to represent insect wings on a Cartesian grid. This method has several advantages. It has a simple algorithm that is easy to implement, the physical quantities around the objects are accurately calculated, and the computational efficiency is higher than that of the immersed boundary method. As shown in Figure 5, the calculation of the VFM requires the placement of a virtual boundary point at the intersection of the boundary of the virtual object and the discrete velocity vector.
In this study, the no-slip condition for velocity and the Neumann boundary condition for pressure were used as the boundary conditions of the virtual object, and each boundary condition can be expressed as follows:
u v b = u w a l l ,
p n = 0 ,
where u w a l l is the velocity of the object on the wall and n is the normal vector to the virtual boundary surface.
The calculation procedure for the VFM is shown in Figure 6. As shown in Figure 6, we consider the case in which an object at point E moves from point E to point D when calculating the distribution function for the next time step at point D. The object at point E cannot pass over the virtual boundary because the virtual boundary point is located at the point that divides the grid width into a   : b between points D and E. Thus, we need to find the distribution function for the pressure at the next time step of point D using the distribution function at the virtual boundary point.
Pressure p v b at the virtual boundary can be calculated with extrapolation, and the numerical precision of the pressure is variable. In this study, we used extrapolation with second-order accuracy and pressure p v b was expressed by pressures p 1 and p 2 at points h 1 and h 2 away from the boundary point in the direction normal to the wall surface, which is presented as follows:
p v b = h 2 2 p 1 h 1 2 p 2 h 2 2 h 1 2 ,
where pressures p 1 and p 2 were interpolated with weights from the surrounding four grid points. In this case, h 1 and h 2 were set to 3 and 2 3 times the grid width, respectively, to prevent the grid points from straddling the virtual boundary surface. The equilibrium distribution function p α e q was calculated by substituting the physical quantities u v b and p v b at the virtual boundary point into Equation (20) as follows:
p α e q t , x v b = ω α p v b + ρ 0 e α ·   u v b + 3 e α ·   u v b 2 2 c 2   u v b 2 2 .
The virtual equilibrium distribution function p α e q t , x E at point E was linearly extrapolated from the equilibrium distribution function p α e q t , x v b , calculated in Equation (25), and the equilibrium distribution function p α e q t , x D at point D from the internal ratio a : b shown in Figure 6. The nonequilibrium part of the virtual distribution function p α n e q t , x E at point E was interpolated with the nonequilibrium part of the distribution function p α n e q t , x D at point D as follows:
p α e q t ,   x E = a + b a p α e q t ,   x v b b a p α e q t ,   x D ,
p α n e q t ,   x E = p α n e q t ,   x D .
However, when the internal ratio a was small, the denominator was small and the calculation would possibly diverge. Therefore, in the case of a < 0.5 , the distribution function is calculated using the equilibrium distribution function and the nonequilibrium part of the distribution function at point C, which is the next point.
Thus, the distribution function p α ( t + δ t , x D ) at point D of the next time step can be calculated using the virtual equilibrium distribution function p α e q t ,   x E and the nonequilibrium part of the distribution function p α n e q t ,   x E at point E as follows:
p α ( t + δ t ,   x D ) = p α e q t ,   x E + 1 1 τ p α n e q t ,   x E .

2.5. Computational Model

Figure 7 shows a schematic view of the computational model used in this study. The computational domain has a size of 40 L × 40 L × 40 L , where the representative length L is the mean chord length C m and the representative speed U is the mean velocity of the wing tip u t i p ¯ . Note that a grid model consisting of four blocks of different grid sizes (four-tiered multiblock method [35]) was used to reduce computation time. In the case of the four-tiered block, each grid width was two times larger than the inner one and the coarsest grid width was eight times larger than the finest one. For the initial conditions of the computational domain, pressure p was set to 1 / 3 and velocity u was set to 0 . For the boundary conditions of the computational domain (outside of block1), the pressure was fixed similarly to the initial condition in the x z plane of y = 40 L , whereas the Neumann boundary condition was used in other planes for pressure and in all planes for velocity.

2.6. Evaluation Parameters

We calculated the lift coefficient C L , thrust coefficient C T , and power coefficient C P W R as the evaluation indices of the fluid in the flapping flight, as described in Section 2.1, to investigate the effects of the figure-eight motion on aerodynamics. They are defined as follows:
C L = F z 1 2 ρ a i r u t i p ¯ 2 S W ,
C T = F y 1 2 ρ a i r u t i p ¯ 2 S W ,
C P W R = f l o c a l · u l o c a l 1 2 ρ a i r u t i p ¯ 3 S W ,
where F y and F z are the forces in the y and z directions, f l o c a l is the force, and u l o c a l is the velocity at a certain point of the wing.
Moreover, we used the Q value, which expresses the vortex structure, and the helicity density h d , which indicates the degree of twisting structure of the swirling flow by the vortex, as evaluation indices of the flow field. Such indices are defined as follows:
Q = 1 2 Ω i j 2 S i j 2 ,
h d = u · ω ,
where u is the velocity vector and ω is the vorticity vector. Moreover, Ω i j is the vorticity tensor and S i j is the deformation velocity tensor, which are expressed as follows:
Ω i j = 1 2 D i j D j i   ,
S i j = 1 2 D i j + D j i   ,
where D is the velocity gradient tensor, which is calculated as follows:
D = u x u y u z v x v y v z w x w y w z = D 11 D 12 D 13 D 21 D 22 D 23 D 31 D 32 D 33 .
In this study, the Q value and helicity density h d are nondimensionalized as follows:
Q = Q U / L 2   ,
h d = h d u · ω   .

3. Results and Discussion

3.1. Grid Independence and Flapping Cycle Convergence Tests

A grid independence test was performed for three patterns with grid sizes of 16 , 32 , and 64   c e l l s / L for the representative length of the finest grid size in the multiblock to investigate the number of grid sizes with numerical reliability in this analysis. The representative speed of the flapping flight U was set to 0.04 , and the resolution was verified when the lift coefficient was stable in the seven flapping cycles of the flapping cycle convergence test. Figure 8 shows the time history of the lift coefficient in seven cycles, and Table 3 shows the cycle-averaged lift coefficient per cycle. As shown in Figure 8, a disturbance was observed at the beginning of the cycle because of the beginning of the wing movement; however, it converged immediately, and little change was found in the time history. In addition, not much difference was observed after the fifth cycle (Table 3). Therefore, the analysis conducted in this study was performed using the data from the sixth cycle.
For the grid independence test, Figure 9 shows the time history of the lift coefficient C L , the pressure component of the lift coefficient C L p , and the viscous stress component the of lift coefficient C L τ over one cycle at each resolution, and Table 4 shows the cycle-averaged values of them C L ¯ , C L p ¯ , and C L τ ¯ . Based on these results, the difference in lift at each resolution decreased as the resolution increased, and not much difference was observed between 32 and 64   c e l l s / L . Therefore, a resolution of 32   c e l l s / L was selected after considering the computational cost and accuracy.

3.2. Example Test for Analyzing the Three-Dimensional Flapping Motion

Fluid forces generated by a rectangular rigid plate in simple motion were compared with those in previous studies to verify the physical validity of the analysis results for a moving rigid plate in a three-dimensional analysis. The motion is oscillating, which is expressed through the following translational displacement (heaving motion) and rotation (flapping motion):
x ( t ) = L sin 2 π f t ,
θ t = π 2 π 4 sin 2 π f t + π 3 ,
where x t is the displacement of the x -coordinate of the center of gravity and θ t is the oscillatory angle. Figure 10 shows a schematic view of the movement of an oscillating plate.
Figure 11 shows the time history of lift coefficient generated from an oscillatory plate. As shown in Figure 11, the simulated values were similar to the results of other studies [36,37]. In addition, the cycle-averaged value of the lift was 0.223 in this study, while the value was 0.22 in a previous study [36], with a difference of 1.3 % and the results of this analysis were nearly consistent with the previous study.

3.3. Effect of Figure-Eight Motion

A comparison of the vortex structure and fluid forces generated by each motion was conducted to investigate the effect of the figure-eight motion on the aerodynamics. Each motion was represented by varying the elevation motion angle θ e which is a rotational motion that represents the figure-eight motion shown in Equation (4). Motion with figure-eight motion was represented by the elevation motion angle with the parameters shown in Section 2.2, and motion without figure-eight motion was represented by θ e = 0 . Equations (41) and (42) show the angles with and without figure-eight motions as follows:
θ e = 70 cos 2 π 2 f t 10 + 10 ,
θ e = 0 .
Figure 12 shows a schematic view of each motion represented by the flapping motion angle θ p and the feathering motion angle θ f shown in Section 2.2 in addition to the elevation motion angle.
Figure 13 shows the isosurfaces of the vortex structure Q = 0.8 formed via each motion at t = 0.25 T , 0.50 T , 0.75 T , and 1.00 T . The isosurfaces were colored on the basis of the normalized helicity density h d . As shown in Figure 13, the normalized helicity density was similar for each motion; however, the lengths of the vortex structures were different and the vortices with a figure-eight motion were longer than those without a figure-eight motion. This is because the vortices with a figure-eight motion stay longer on the top surface of the wings than those without a figure-eight motion, whereas the vortices without a figure-eight motion disappear immediately.
This difference in the vortex structure indicates a difference in the flow field, which indicates a difference in the relative velocity vector between the wings and the fluid. The difference in velocity vectors affects forces such as lift and thrust generated by the wings. Therefore, the time histories of the fluid forces generated via each motion were investigated. Figure 14 shows the time histories of the lift coefficient C L , the thrust coefficient C T , and the power coefficient C P W R over one stroke cycle for each motion, and Table 5 shows the cycle-averaged values of coefficients C L ¯ , C T ¯ , and C P W R ¯ and the ratio of the lift coefficient to the power coefficient C L ¯ / C P W R ¯ for each motion.
As shown in Figure 14a, a comparison between the lift coefficient with and without figure-eight motions showed a different trend throughout the entire cycle and an essentially different trend in the lift coefficient in the first half of each stroke. We examined the vortex structure and pressure coefficient distribution at t = 0.1 T , which corresponds to the initial stage of the downward motion of the wings (downstroke), to further discuss the cause of this large difference in the lift coefficient. Figure 15 shows the isosurfaces of the vortex structure ( Q = 0.8 ) formed by each motion. The isosurfaces were colored on the basis of the normalized pressure coefficient C p . As shown in Figure 15, the motion with a figure-eight motion produced a larger leading-edge vortex and wing tip vortex on the upper surface of the wings than the motion without a figure-eight motion. This difference could be attributed to the fact that the z -directional motion of the figure-eight motion increased the angle relative to the motion direction (angle of attack) of the wings by adding a downward motion in the first half of the downstroke. Consequently, a larger negative pressure was produced on the upper surface of the wings, and the pressure difference between the positive pressure on the lower surface and the negative pressure on the upper surface of the wings generated a larger lift force than the motion without a figure-eight motion. In addition, as shown in Figure 14a, the relation of the lift coefficient in each motion reversed after the middle of the downstroke. This result was based on the abovementioned fact and was considered to be due to the fact that the figure-eight motion in the z -direction increases the angle of attack of the wings by adding upward motion after the middle of the downstroke. As shown in Table 5, the average lift coefficient with a figure-eight motion was approximately 26 % larger than that without a figure-eight motion.
As shown in Figure 14b, a comparison between the thrust coefficient with and without figure-eight motions showed a similar trend throughout the entire cycle, and the cycle-averaged thrust coefficient in each was small (Table 5). Therefore, the effect of the figure-eight motion on the thrust coefficient was small. This effect could be attributed to the fact that the thrust coefficient was offset by the upward motion of the wings (upstroke) and downstroke.
As shown in Figure 14c, a comparison between the power coefficient with and without figure-eight motions showed a different trend throughout the entire cycle, and this difference could be attributed to the difference in the lift coefficient for each motion. As shown in Table 5, a comparison of the cycle-averaged power coefficient in each motion confirmed that the motion with figure-eight motion required more power. This finding could be attributed to the motion in the z -direction.
As shown in Table 5, the ratio of the lift coefficient to the power coefficient with a figure-eight motion was approximately 17% higher than that without a figure-eight motion. Therefore, these results indicated that a flight with a figure-eight motion consumed more power than a flight without a figure-eight motion but generated more lift, and it may be a more efficient way to fly while hovering.

3.4. Effect of Various Figure-Eight Motions and Reynolds Number

We compared the fluid forces generated via various figure-eight motions in addition to the motions shown in Section 3.3 to investigate the more efficient figure-eight motion patterns in generating lift during the hovering motion and the relationship between figure-eight motion and Reynolds number. These motions were represented as eight types of figure-eight motions by varying the initial phase φ e 0 of the elevation motion angle shown in Equation (4) from 0 ° to 315 ° in increments of 45 ° . The flapping motion angle and the feathering motion angle were the same as those shown in Section 2.2. Figure 16 shows a schematic view of each motion. Moreover, we investigated the dependence of the Reynolds number on the figure-eight motion.
Table 6 shows the cycle-averaged values of the lift coefficient and the ratio of the lift coefficient to the power coefficient for each motion at each Reynolds number. Figure 17 shows the relation between the cycle-averaged power coefficient and the cycle-averaged lift coefficient. Figure 17 shows the cycle-averaged power coefficient on the horizontal axis and the cycle-averaged lift coefficient on the vertical axis. Hence, comparing each type of motion, the motion that generates more lift for a certain amount of power (more efficient motion) is located in the upper left corner. This positional relationship through the comparison of motions can be described similarly by the slope of a line connecting the origin and one point of the results. Thus, the motion with a greater slope of the line is more efficient in generating lift. The dotted line shown in Figure 17 was used to compare each figure-eight motion and without a figure-eight motion, which was the line connecting the value of without-figure-eight motion and the origin of the graph. In this study, the line was known as “without-8-line”. The motion located to the upper left of “without-8-line” generated lift more efficiently than without a figure-eight motion, whereas the motion located to the lower right of “without-8-line” was less efficient than without a figure-eight motion.
As shown in Figure 16, the motions with φ e 0 = 0 ° and φ e 0 = 180 ° were recognized as U-shaped motions. However, as shown in Table 6, the most efficient motions were not the U-shaped motions but the figure-eight motions with φ e 0 = 45 ° at any Reynolds number. Therefore, in this study, the figure-eight motion was defined as the motion with θ e 0 ° , without discussing the U-shaped motion independently.
As shown in Table 6, the motion with a smaller initial phase of the elevation motion angle ( φ e 0 90 ° ) had a higher ratio of the lift coefficient to the power coefficient, and the motion could generate lift more efficiently. The motion with φ e 0 = 90 ° at R e = 33.5 and φ e 0 = 90 ° at other Reynolds numbers was the most efficient in generating lift. As shown in Figure 16, this result could be attributed to the fact that the angle of attack was larger than the original figure-eight motion with φ e 0 = 10 ° by increasing the vertical motion of the wings in each stroke. On the contrary, the motion with φ e 0 = 90 ° had a larger vertical motion and a larger angle of attack compared with the motion with φ e 0 = 45 ° ; however, the upward motion was larger. Therefore, the motion was less efficient than the motion with φ e 0 = 45 ° when it was averaged over the entire stroke. However, horizontal forces dominate rather than vertical forces toward motion at low Reynolds numbers; therefore, increasing the vertical motion of the flapping motion and the angle of attack are important for generating more lift. Therefore, the motion with φ e 0 = 90 ° at R e = 33.5 was the most efficient in generating lift. Moreover, the motion patterns of the most efficient in generating lift at each Reynolds number were investigated in detail. Figure 18 shows the relation between the initial phase of the elevation angle of the most efficient motion in generating lift and Reynolds number. Note that the most efficient motion in generating lift was calculated from the tangential point where C L   ¯ / C P W R ¯ was the greatest slope on the elliptical approximation shown in Figure 17. As shown in Figure 18, the initial phase of the elevation angle of the most efficient motion in generating lift increased as the Reynolds number decreased, whereas the angle decreased and approached zero as the Reynolds number increased.
In addition, as shown in Table 6 and Figure 17f, all motions generated more lift as the Reynolds number increased, and more lift was generated compared with the power. The above discussion is based on nondimensionalized values without considering differences in the speed of the wings (flapping frequency) to simplify the comparison of the lift and power generated by each motion. Then, the physical values of the lift and power coefficients were summarized to investigate the effect of the flapping frequency on aerodynamic characteristics. Table 7 shows the cycle-averaged values of the force in the z direction F z , power P , and the ratio of the lift coefficient to the power coefficient F z / P generated via the motion without a figure-eight motion at each Reynolds number. As shown in Table 7, a higher flapping frequency generated more power than lift.
Moreover, as shown in Figure 17, the eight types of figure-eight motions showed an elliptical distribution in the relation between the lift coefficient and power coefficient. Therefore, we discuss the relation between the lift coefficient and power coefficient by making an elliptical approximation along the distribution. As shown in Figure 17a–e, the area surrounded by the ellipse and the W8 line (light blue shown in the figure) decreased as the Reynolds number increased, and the number of motions that generate lift more efficiently than without the figure-eight motion also decreased. This decrease indicates that the effect of figure-eight motion decreased as the Reynolds number increased. The elliptical shape, which is primarily related to the decrease in area, was focused on to discuss this difference in effect in more detail. Therefore, the effect of the change in the Reynolds number on the figure-eight motion was determined by calculating the aspect ratio and inclination angle of the ellipse at each Reynolds number. As shown in Figure 19, the long and short sides of the ellipse for each Reynolds number were defined as l l and l s , respectively, and the aspect ratio A R can be expressed as A R = l l / l s , where A R is the ratio of the long side to the short side of the ellipse. In addition, the inclination angle of the ellipse and the location of the center of the ellipse in Figure 17 were defined as θ e l and O ( x ,   y ) . Table 8 shows the values of the elliptic approximation, and Figure 20 shows the relation between the aspect ratio of the ellipse and the Reynolds number.
As shown in Table 8 and Figure 20, the inclination angle of the ellipse did not change remarkably; the aspect ratio increased with increasing Reynolds number, and the aspect ratio converged when the Reynolds number increased above a certain value. The increase in the aspect ratio indicated that the elliptical approximation of various figure-eight motions approached the straight “without-8-line” shown in Figure 17, which indicates that the difference between the motions with and without a figure-eight motion decreased. In addition, the increase in the Reynolds number was synonymous with the increase in insect size. Therefore, varying the elevation motion angle that represents the figure-eight motion was not very effective for insects under high Reynolds number conditions (large size). Particularly, the figure-eight motion may be a vital mechanism for insects to generate lift more efficiently under low Reynolds number conditions (small size).
In this study, insect wings were analyzed as rectangular rigid plates without thickness. Regarding the shape, Kirishna et al. [38] reported that flight efficiency did not change much depending on the shape of the wings between a rectangular model with a model based on the actual shape of a blowfly. However, the changes in the shape of the wings may affect the efficiency of each figure-eight motion pattern. Therefore, in future studies of insects, it is necessary to evaluate figure-eight motion by considering not only the shape of the wings but also the flexibility of the wings and other parameters.

4. Conclusions

In this study, we investigated the more efficient figure-eight motion patterns in generating lift during the hovering motion related to the Reynolds number. For this purpose, we compared the ratios of the cycle-averaged lift coefficient to the power coefficient generated via each motion by varying the elevation motion angle and Reynolds number. Consequently, the motion with a smaller initial phase of the elevation motion angle ( φ e 0 90 ° ) could generate lift more efficiently at all Reynolds numbers. In addition, the initial phase of the elevation angle of the most efficient motion in generating lift increased as the Reynolds number decreased, whereas the angle decreased and approached φ e 0 = 0 as the Reynolds number increased. Moreover, varying the elevation motion angle that represents the figure-eight motion was effective for insects under low Reynolds number conditions and was not very effective for insects under high Reynolds number conditions. Therefore, the figure-eight motion may be a vital mechanism for insects under low Reynolds number conditions to generate lift more efficiently.

Author Contributions

Conceptualization, M.Y. and T.F.; methodology, M.Y. and T.F.; software, M.Y. and T.F.; validation, M.Y.; formal analysis, M.Y.; investigation, M.Y.; resources, T.F.; data curation, M.Y.; writing—original draft preparation, M.Y.; writing—review and editing, M.Y. and T.F.; visualization, M.Y.; supervision, T.F.; project administration, T.F.; funding acquisition, T.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

Dataset available on request from the authors.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Yahyanejad, S.; Rinner, B. A fast and mobile system for registration of low-altitude visual and thermal aerial images using multiple small-scale UAVs. ISPRS J. Photogramm. Remote Sens. 2015, 104, 189–202. [Google Scholar] [CrossRef]
  2. Xiao, S.; Hu, K.; Huang, B.; Huichao, D.; Ding, X. A Review of Research on the Mechanical Design of Hoverable Flapping Wing Micro-Air Vehicle. J. Bionic Eng. 2021, 18, 1235–1254. [Google Scholar] [CrossRef]
  3. Keennon, M.; Klingebiel, K.; Won, H.; Andriukov, A. Development of the Nano Hummingbird: A Tailless Flapping Wing Micro Air Vehicle. In Proceedings of the 50th AIAI Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Nashville, TN, USA, 9–12 January 2012. Paper No. AIAI 2012-0588. [Google Scholar]
  4. Wood, R.J.; Finio, B.; Karpelson, M.; Ma, K.; Perez-Arancibia, N.O.; Sreetharan, P.S.; Tanaka, H.; Whitney, J.P. Progress on ‘Pico’ air vehicles. Int. J. Robot. Res. 2012, 31, 1292–1302. [Google Scholar] [CrossRef]
  5. Ang, H.S.; Xiao, T.H.; Duan, W.B. Flight mechanism and design of biomimetic micro air vehicles. Sci. China Ser. E 2009, 52, 3722–3728. [Google Scholar] [CrossRef]
  6. Cai, G.; Dias, C.; Seneviratne, L. A Survey of Small-Scale Unmanned Aerial Vehicles: Recent Advances and Future Development Trends. Unmanned Syst. 2014, 2, 175–199. [Google Scholar] [CrossRef]
  7. Hassanalian, M.; Abdelkefi, A. Classifications, applications, and design challenges of drones: A review. Prog. Aerosp. Sci. 2017, 91, 99–131. [Google Scholar] [CrossRef]
  8. Mueller, T.J.; Delaurier, J.D. An overview of micro air vehicle aerodynamics. Fixed Flapping Wings Aerodyn. Micro Air Veh. Appl. 2001, 195, 1–10. [Google Scholar]
  9. Bohorquez, F.; Samuel, P.; Sirohi, J.; Pines, D.; Pudd, L. Design, Analysis, and Hover Performance of a Rotary Wing Micro Air Vehicle. J. Am. Helicopter Soc. 2003, 48, 80–90. [Google Scholar] [CrossRef]
  10. Floreano, D.; Wood, R.J. Science, technology, and the future of small autonomous drones. Nature 2015, 521, 460–466. [Google Scholar] [CrossRef]
  11. Abas, M.F.B.; Rafic, A.S.B.M.; Yusoft, H.B.; Ahmad, K.A.B. Flapping wing micro-aerial-vehicle: Kinematic, membranes, and flapping mechanism of ornithopter and insect flight. Chin. J. Aeronaut. 2016, 29, 1159–1177. [Google Scholar] [CrossRef]
  12. Ellington, C.P.; Berg, C.V.D.; Willmott, A.P.; Thomas, A.L.R. Leading-edge vortices in insect flight. Nature 1996, 384, 626–630. [Google Scholar] [CrossRef]
  13. Bomphrey, R.J.; Nakata, T.; Henningsson, P.; Lin, H. Flight of the dragonflies and damselflies, philosophical transaction of the royal society b. Biol. Sci. 2016, 371, 20150389. [Google Scholar] [CrossRef] [PubMed]
  14. Liu, H.; Aono, H. Size effect on insect hovering aerodynamics: An integrated computational study. Bioinspir. Biomim. 2009, 4, 015002. [Google Scholar] [CrossRef] [PubMed]
  15. Weis-Fogh, T. Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production. J. Exp. Biol. 1973, 59, 169–230. [Google Scholar] [CrossRef]
  16. Cheng, X.; Sum, M. Very small insects use novel wing flapping and drag principle to generate the weight-supporting vertical force. J. Fluid Mech. 2018, 855, 646–670. [Google Scholar] [CrossRef]
  17. Lyu, Y.Z.; Zhu, H.J.; Sun, M. Flapping-mode changes aerodynamic mechanisms in miniature insects. Phys. Rev. E 2019, 99, 012419. [Google Scholar] [CrossRef] [PubMed]
  18. Conn, A.T.; Ling, C.S.; Burgess, S.C. Biomimetic Analysis of Insect Wing Kinematics for Flapping MAVs. Int. J. Micro Air Veh. 2011, 3, 1–11. [Google Scholar] [CrossRef]
  19. Aghav, H. Effects of stroke deviation on the aerodynamics of the smallest flying insects. J. Eng. Math. 2022, 137, 4. [Google Scholar] [CrossRef]
  20. Lehmann, F.O.; Pick, S. The aerodynamic benefit of wing–wing interaction depends on stroke trajectory in flapping insect wings. J. Exp. Biol. 2007, 210, 1362. [Google Scholar] [CrossRef]
  21. Galinski, C.; Zbikowski, R. Insect-like flapping wing mechanism based on a double spherical Scotch yoke. J. R. Soc. Interface 2005, 2, 223. [Google Scholar] [CrossRef]
  22. Ishihara, D. Role of fluid-structure interaction in generating the characteristic tip path of a flapping flexible wing. Phys. Rev. E 2018, 98, 032411. [Google Scholar] [CrossRef]
  23. Hu, F.; Liu, X. Effects of stroke deviation on hovering aerodynamic performance of flapping wings. Phys. Fluids 2019, 31, 111901. [Google Scholar] [CrossRef]
  24. Hu, F.; Wang, Y.; Li, D.; Liu, X. Effects of asymmetric stroke deviation on the aerodynamic performance of flapping wing. Proc. Inst. Mech. Eng. G J. Aerosp. Eng. 2023, 237, 480–499. [Google Scholar] [CrossRef]
  25. Bomphrey, R.J.; Nakata, T.; Phillips, N.; Walker, S.M. Smart wing rotation and trailing-edge vortices enable high frequency mosquito flight. Nature 2017, 544, 92–95. [Google Scholar] [CrossRef] [PubMed]
  26. Aono, H.; Liang, F.; Liu, H. Near and far-field aerodynamics in insect hovering flight: An integrated computational study. J. Exp. Biol. 2008, 211, 239–257. [Google Scholar] [CrossRef] [PubMed]
  27. Yamauchi, K.; Fukui, T.; Morinishi, K. Numerical simulation of influences of the body’s presence on flow around the wings in insect flapping flight. In Proceedings of the AEME-JSME-KSME 2019 Joint Fluids Engineering Conference, San Francisco, CA, USA, 28 July–1 August 2019. Paper No. AJKFLUIDS2019-5176. [Google Scholar]
  28. Izham, M.; Fukui, T.; Morinishi, K. Application of Regularized Lattice Boltzmann Method for Incompressible Flow Simulation at High Reynolds Number and Flow with Curved Boundary. J. Fluid Sci. Technol. 2011, 6, 812–822. [Google Scholar] [CrossRef]
  29. Morinishi, K.; Fukui, T. Parallel computational of turbulent flows using moment base lattice Boltzmann method. Int. J. Comput. Fluid Dyn. 2016, 30, 363–369. [Google Scholar] [CrossRef]
  30. Ladd, A.J.C. Numerical simulations of particulate suspensions via a discretized Boltzmann equation part 1 theoretical foundation. J. Fluid Mech. 1994, 271, 285–309. [Google Scholar] [CrossRef]
  31. Chen, S.; Doolen, G.D. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 1998, 30, 329–364. [Google Scholar] [CrossRef]
  32. Morinishi, K.; Fukui, T. An Eulerian approach for fluid structure interaction problems. Comput. Fluids 2012, 65, 92–98. [Google Scholar] [CrossRef]
  33. Tanno, I.; Morinishi, K.; Matsuno, K.; Nishida, H. Validation of virtual flux method for forced convection flow. JSME Int. J. Ser. B 2006, 49, 1141–1148. [Google Scholar] [CrossRef]
  34. Kawaguchi, M.; Fukui, T.; Morinishi, K. Comparative study of the virtual flux method and immersed boundary method coupled with regularized lattice Boltzmann method for suspension flow simulations. Comput. Fluids 2022, 246, 105615. [Google Scholar] [CrossRef]
  35. Yu, D.; Mei, R.; Shyy, W. A multi-block lattice Boltzmann method for viscous fluid flows. Int. J. Numer. Methods Fluids 2002, 39, 99–120. [Google Scholar] [CrossRef]
  36. Trizila, P.; Kang, C.K.; Aono, H.; Shyy, W. Low-Reynolds-Number Aerodynamics of a Flapping Rigid Flat Plate. AIAA J. 2011, 49, 806–823. [Google Scholar] [CrossRef]
  37. Wang, S.; He, G.; Zhang, X. Lift enhancement on spanwise oscillating flat-plates in low-Reynolds-number flow. Phys. Fluids 2015, 27, 1–19. [Google Scholar] [CrossRef]
  38. Krishna, S.; Cho, M.; Wehmann, H.N.; Engels, T.; Lehmann, F.O. Wing Design in Flies: Properties and Aerodynamic Function. Insects 2020, 11, 466. [Google Scholar] [CrossRef]
Figure 1. Computational model of an (a) insect and (b) wing. Only the two wings shown in gray were analyzed without body reproduction. The red line represents the wing length.
Figure 1. Computational model of an (a) insect and (b) wing. Only the two wings shown in gray were analyzed without body reproduction. The red line represents the wing length.
Biomimetics 09 00249 g001
Figure 2. Definition of the flapping motion of the insect in (a) bird’s eye view, (b) x z plane, and (c) y z plane. The wings, which are the analysis objects, are shown in gray. On the other hand, the body, which is not the analysis object, is shown by the dotted line.
Figure 2. Definition of the flapping motion of the insect in (a) bird’s eye view, (b) x z plane, and (c) y z plane. The wings, which are the analysis objects, are shown in gray. On the other hand, the body, which is not the analysis object, is shown by the dotted line.
Biomimetics 09 00249 g002
Figure 3. Representation of insect flapping motion. (a) Time history of the positional angle, feathering angle, and elevation angle. (b) Trajectory of wing tips. In Figure 3b,The gray solid and dotted lines represent the trajectory of the wings tips and the body of the insect, respectively.
Figure 3. Representation of insect flapping motion. (a) Time history of the positional angle, feathering angle, and elevation angle. (b) Trajectory of wing tips. In Figure 3b,The gray solid and dotted lines represent the trajectory of the wings tips and the body of the insect, respectively.
Biomimetics 09 00249 g003
Figure 4. 3D27V model for the three-dimensional lattice Boltzmann method. The numbers (0~26) represent the 27 directions of the 3D27Vmodel.
Figure 4. 3D27V model for the three-dimensional lattice Boltzmann method. The numbers (0~26) represent the 27 directions of the 3D27Vmodel.
Biomimetics 09 00249 g004
Figure 5. Virtual boundary points. The gray and white areas represent the interior of the object and fluid, respectively.
Figure 5. Virtual boundary points. The gray and white areas represent the interior of the object and fluid, respectively.
Biomimetics 09 00249 g005
Figure 6. Schematic view of the physical quantity calculation method at the virtual boundary point. The gray and white areas represent the interior of the object and fluid, respectively.
Figure 6. Schematic view of the physical quantity calculation method at the virtual boundary point. The gray and white areas represent the interior of the object and fluid, respectively.
Biomimetics 09 00249 g006
Figure 7. Schematic view of the computational model: (a) bird’s eye view, (b) x z plane, and (c) x y plane. The representative length L is the mean chord length C m shown in Section 2.1. In Figure 7a, the wings and grid model are shown in green and black or gray lines, respectively. In Figure 7b,c, the wings are shown in gray areas, and block areas of grid model consisting of different grid sizes are represented by black, green, red, and blue lines.
Figure 7. Schematic view of the computational model: (a) bird’s eye view, (b) x z plane, and (c) x y plane. The representative length L is the mean chord length C m shown in Section 2.1. In Figure 7a, the wings and grid model are shown in green and black or gray lines, respectively. In Figure 7b,c, the wings are shown in gray areas, and block areas of grid model consisting of different grid sizes are represented by black, green, red, and blue lines.
Biomimetics 09 00249 g007
Figure 8. Time history of the lift coefficient in seven cycles with U = 0.04 .
Figure 8. Time history of the lift coefficient in seven cycles with U = 0.04 .
Biomimetics 09 00249 g008
Figure 9. Time history of the (a) lift coefficient, (b) pressure component of lift coefficient, and (c) viscous stress component of lift coefficient for three resolutions.
Figure 9. Time history of the (a) lift coefficient, (b) pressure component of lift coefficient, and (c) viscous stress component of lift coefficient for three resolutions.
Biomimetics 09 00249 g009
Figure 10. Schematic view of the movement of an oscillating plate.
Figure 10. Schematic view of the movement of an oscillating plate.
Biomimetics 09 00249 g010
Figure 11. Time history of the lift coefficient of an oscillating plate. The result was compared with those of Trizila [36] and Wang et al. [37].
Figure 11. Time history of the lift coefficient of an oscillating plate. The result was compared with those of Trizila [36] and Wang et al. [37].
Biomimetics 09 00249 g011
Figure 12. Trajectory of the wing tip in each motion. (a) With figure-eight motion and (b) without figure-eight motion. The gray line represents the trajectory of the wing tip.
Figure 12. Trajectory of the wing tip in each motion. (a) With figure-eight motion and (b) without figure-eight motion. The gray line represents the trajectory of the wing tip.
Biomimetics 09 00249 g012
Figure 13. Vortex structures and normalized helicity density over one stroke cycle for each motion. (a) With and (b) without figure-eight motions at t = 0.25 T ; (c) with and (d) without figure-eight motions at t = 0.50 T ; (e) with and (f) without figure-eight motions at t = 0.75 T ; and (g) with and (h) without figure-eight motions at t = 1.00 T .
Figure 13. Vortex structures and normalized helicity density over one stroke cycle for each motion. (a) With and (b) without figure-eight motions at t = 0.25 T ; (c) with and (d) without figure-eight motions at t = 0.50 T ; (e) with and (f) without figure-eight motions at t = 0.75 T ; and (g) with and (h) without figure-eight motions at t = 1.00 T .
Biomimetics 09 00249 g013
Figure 14. Time histories of (a) lift coefficient, (b) thrust coefficient, and (c) power coefficient over one stroke cycle for each motion.
Figure 14. Time histories of (a) lift coefficient, (b) thrust coefficient, and (c) power coefficient over one stroke cycle for each motion.
Biomimetics 09 00249 g014
Figure 15. Vortex structures ( Q = 0.8 ) and normalized pressure coefficient for each motion in the downstroke at t = 0.1 T . Vortex structures ( Q = 0.8 ) and normalized pressure coefficients from (a) with and (b) without figure-eight motions. Normalized pressure coefficient at the wing tip from (c) with and (d) without figure-eight motions.
Figure 15. Vortex structures ( Q = 0.8 ) and normalized pressure coefficient for each motion in the downstroke at t = 0.1 T . Vortex structures ( Q = 0.8 ) and normalized pressure coefficients from (a) with and (b) without figure-eight motions. Normalized pressure coefficient at the wing tip from (c) with and (d) without figure-eight motions.
Biomimetics 09 00249 g015
Figure 16. Trajectory of the wing tip in each motion: (a) φ e 0 = 0 ° , (b) φ e 0 = 45 ° ; (c) φ e 0 = 90 ° ; (d) φ e 0 = 135 ° ; (e) φ e 0 = 180 ° ; (f) φ e 0 = 225 ° ; (g) φ e 0 = 270 ° ; and (h) φ e 0 = 315 ° . The gray line represents the trajectory of the wing tip.
Figure 16. Trajectory of the wing tip in each motion: (a) φ e 0 = 0 ° , (b) φ e 0 = 45 ° ; (c) φ e 0 = 90 ° ; (d) φ e 0 = 135 ° ; (e) φ e 0 = 180 ° ; (f) φ e 0 = 225 ° ; (g) φ e 0 = 270 ° ; and (h) φ e 0 = 315 ° . The gray line represents the trajectory of the wing tip.
Biomimetics 09 00249 g016
Figure 17. Relation between the cycle-averaged lift coefficient and power coefficient for each motion at each Reynolds number. (a) R e = 33.5 ; (b) R e = 67 ; (c) R e = 134 ; (d) R e = 268 ; (e) R e = 546 ; (f) all Reynolds numbers. Light blue areas shown in the figure (Biomimetics 09 00249 i001) indicate more efficiency than without figure-eight motions.
Figure 17. Relation between the cycle-averaged lift coefficient and power coefficient for each motion at each Reynolds number. (a) R e = 33.5 ; (b) R e = 67 ; (c) R e = 134 ; (d) R e = 268 ; (e) R e = 546 ; (f) all Reynolds numbers. Light blue areas shown in the figure (Biomimetics 09 00249 i001) indicate more efficiency than without figure-eight motions.
Biomimetics 09 00249 g017
Figure 18. Relation between the initial phase of elevation angle of the most efficient motion in generating lift and Reynolds number.
Figure 18. Relation between the initial phase of elevation angle of the most efficient motion in generating lift and Reynolds number.
Biomimetics 09 00249 g018
Figure 19. Schematic view of the elliptic approximation.
Figure 19. Schematic view of the elliptic approximation.
Biomimetics 09 00249 g019
Figure 20. Relation between the aspect ratio and Reynolds number.
Figure 20. Relation between the aspect ratio and Reynolds number.
Biomimetics 09 00249 g020
Table 1. Physical properties of an insect.
Table 1. Physical properties of an insect.
NameSymbolValue
Body length l b 2.78   [ m m ]
Mean chord length C m 0.78   [ m m ]
Wing length R 2.39   [ m m ]
Mean velocity of wing tip u t i p ¯ 2.58   [ m / s ]
Wingbeat frequency f 218   [ H z ]
Kinematic viscosity of air ν a i r 1.5 × 10 5   [ m 2 / s ]
Density of air ρ a i r 1.2   [ k g / m 3 ]
Table 2. Advection parameters of the particles.
Table 2. Advection parameters of the particles.
α e α e α ω α
0 0 , 0 , 0 0 8 / 27
1 6 ± 1 , 0 , 0 , 0 , ± 1 , 0 , 0 , 0 , ± 1 c 2 / 27
7 18 ± 1 , ± 1,0 , ± 1,0 , ± 1 , 0 , ± 1 , ± 1 2 c 1 / 54
19 26 ± 1 , ± 1 , ± 1 3 c 1 / 216
Table 3. Cycle-averaged lift coefficient in seven cycles.
Table 3. Cycle-averaged lift coefficient in seven cycles.
Cycle C L ¯
10.484
20.379
3 0.390
4 0.387
50.383
60.382
70.383
Table 4. Cycle-averaged lift coefficient, pressure component of lift coefficient, and viscous stress component of lift coefficient for three resolutions.
Table 4. Cycle-averaged lift coefficient, pressure component of lift coefficient, and viscous stress component of lift coefficient for three resolutions.
Resolution C L ¯ C L p ¯ C L τ ¯
16   c e l l s / L 0.3730.4080.038
32   c e l l s / L 0.3820.4270.048
64   c e l l s / L 0.3920.4420.052
Table 5. Cycle-averaged values of lift coefficient, thrust coefficient, power coefficient, and the ratio of the lift coefficient to the power coefficient. (A) With figure-eight motion and (B) without figure-eight motion.
Table 5. Cycle-averaged values of lift coefficient, thrust coefficient, power coefficient, and the ratio of the lift coefficient to the power coefficient. (A) With figure-eight motion and (B) without figure-eight motion.
Case C L ¯ C T ¯ C P W R ¯ C L ¯ / C P W R ¯
(A)0.382 3.45 × 10 3 0.3801.005
(B)0.302 2.54 × 10 2 0.3520.857
Table 6. Cycle-averaged values of lift coefficient and the ratio of the lift coefficient to the power coefficient for each motion at each Reynolds number. (A) With figure-eight motion; (B) without figure-eight motion. Cycle-averaged values with a figure-eight motion (Biomimetics 09 00249 i001) indicate more efficiency than those without figure-eight motions.
Table 6. Cycle-averaged values of lift coefficient and the ratio of the lift coefficient to the power coefficient for each motion at each Reynolds number. (A) With figure-eight motion; (B) without figure-eight motion. Cycle-averaged values with a figure-eight motion (Biomimetics 09 00249 i001) indicate more efficiency than those without figure-eight motions.
C L ¯ C L ¯ / C P W R ¯
Case Re33.56713426853633.567134268536
φ e 0
(A)00.2700.3030.3540.3960.4210.5090.7260.9611.1651.335
450.3500.3980.4630.5220.5580.6130.8491.0751.2461.358
900.4010.4580.5200.5810.6130.6240.8381.0191.1591.223
1350.3830.4480.5160.5720.5890.5420.7330.8951.0141.069
1800.2960.3630.4170.4600.4830.4130.5940.7410.8510.922
2250.2030.2520.3120.3410.3360.3060.4630.6370.7570.803
2700.1770.2020.2460.2740.2780.3060.4470.6230.7660.864
3150.2020.2260.2650.2940.3100.3840.5580.7680.9561.116
(B)-0.2200.2510.3020.3500.3790.4290.6230.8571.0651.207
Table 7. Cycle-averaged values of the force in the z direction F z , power P , and the ratio of the force in the z direction to the power F z / P without figure-eight motion at each Reynolds number.
Table 7. Cycle-averaged values of the force in the z direction F z , power P , and the ratio of the force in the z direction to the power F z / P without figure-eight motion at each Reynolds number.
R e [−] F z [μN] P [μW] F z / P [−]
33.50.2050.3080.665
670.9331.932 0.483
1344.49513.53 0.332
26820.84101.0 0.206
53690.32772.30.117
Table 8. Respective parameters of the elliptic approximation for each Reynolds number.
Table 8. Respective parameters of the elliptic approximation for each Reynolds number.
R e O ( x , y ) l l l s θ e l A R
33.5(0.623, 0.284)0.1310.079541.658
67(0.511, 0.329)0.154 0.073 542.110
134(0.468, 0.385)0.176 0.068 532.588
268(0.437, 0.429)0.197 0.068 532.897
536(0.415, 0.446)0.2080.071532.930
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yoshida, M.; Fukui, T. Numerical Simulation of the Advantages of the Figure-Eight Flapping Motion of an Insect on Aerodynamics under Low Reynolds Number Conditions. Biomimetics 2024, 9, 249. https://doi.org/10.3390/biomimetics9040249

AMA Style

Yoshida M, Fukui T. Numerical Simulation of the Advantages of the Figure-Eight Flapping Motion of an Insect on Aerodynamics under Low Reynolds Number Conditions. Biomimetics. 2024; 9(4):249. https://doi.org/10.3390/biomimetics9040249

Chicago/Turabian Style

Yoshida, Masato, and Tomohiro Fukui. 2024. "Numerical Simulation of the Advantages of the Figure-Eight Flapping Motion of an Insect on Aerodynamics under Low Reynolds Number Conditions" Biomimetics 9, no. 4: 249. https://doi.org/10.3390/biomimetics9040249

APA Style

Yoshida, M., & Fukui, T. (2024). Numerical Simulation of the Advantages of the Figure-Eight Flapping Motion of an Insect on Aerodynamics under Low Reynolds Number Conditions. Biomimetics, 9(4), 249. https://doi.org/10.3390/biomimetics9040249

Article Metrics

Back to TopTop