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

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 (φe0≤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.


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 leadingedge 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 (Re ≤ 130) [17,25].Moreover, insects with high Reynolds numbers (Re ≈ 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 figureeights) 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.

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.

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. shown in Figure 1 is the wing area.Moreover, the Reynolds number used in this study was defined as follows: Further, the analysis in Sections 3.1 and 3.3 was performed with  ≈ 134 based on the values shown in Table 1.In addition, the mean velocity of the wing tip ( = 0.645, 1.29, 2.58, 5.16, and 10.32) was varied and analyzed as  ≈ 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 ( = 54.5, 109, 218, 436, and 872) was similarly varied when the mean velocity of the wing tip ( = 0.645, 1.29, 2.58, 5.16, and 10.32) was varied.

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 -axis, wing axis  , and -axis, respectively.Note that the lapping motion of an insect Moreover, the Reynolds number used in this study was defined as follows: Further, the analysis in Sections 3.1 and 3.3 was performed with Re ≈ 134 based on the values shown in Table 1.In addition, the mean velocity of the wing tip (u tip = 0.645, 1.29, 2.58, 5.16, and 10.32) was varied and analyzed as Re ≈ 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 tip = 0.645, 1.29, 2.58, 5.16, and 10.32) was varied.

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 < φ e0 < 180 and 180 < φ e0 < 360), except for the three line-shaped (θ e = 0 • ) and U-shaped (φ e0 = 0 • and φ e0 = 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".
where θ p is the flapping angle, θ f is the feathering angle, and θ e is the elevation angle; θ p,amp , θ f ,amp , and θ e,amp are each amplitude; φ p , φ f , and φ e are each initial position; φ e0 is the initial phase of the elevation angle and f is the flapping frequency.In this study, θ p,amp , θ f ,amp , θ e,amp , φ p , φ f , φ e , φ e0 , and f were set to 70      , and 218 Hz, 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.

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  in the lattice Boltzmann equation is expressed using a discrete velocity vector  .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: ≈   +   +    , (5) where  is the weight coefficient and  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.

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: 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,  ,  , and  in Equation ( 5) are constants that satisfy the following relations: In addition, a 0 , b i , and c ij in Equation ( 5) are constants that satisfy the following relations: where ρ is the fluid density, ρu i is the moment, and Π neq ij 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: e αi e αj c 2 − where the first term on the right is the equilibrium distribution function f eq α and the second term on the right is the nonequilibrium part of the distribution function f neq α .Thus, they can be replaced as follows: e αi e αj c 2 − Further, the time evolution equation of the lattice Boltzmann equation can be expressed as follows: where τ is the relaxation time, which is defined as follows: 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 α : where c s is the speed of sound, which is defined as follows: Thus, the pressure p, velocity component u i , and nonequilibrium part of the stress tensor Π neq ij can be expressed as follows: Furthermore, the equilibrium distribution function p eq α can be expressed as follows: Thus, the time evolution equation of the lattice Boltzmann equation with the incompressibility formulation is defined as follows:

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.

𝑐 = √3
. (16) Thus, the pressure , velocity component  , and nonequilibrium part of the stress tensor  can be expressed as follows: Furthermore, the equilibrium distribution function  can be expressed as follows: Thus, the time evolution equation of the lattice Boltzmann equation with the incompressibility formulation is defined as follows:

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: 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: where u wall 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.
where  is the velocity of the object on the wall and  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  ∶  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  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  was expressed by pressures  and  at points ℎ and ℎ away from the boundary point in the direction normal to the wall surface, which is presented as follows: where pressures  and  were interpolated with weights from the surrounding four grid points.In this case, ℎ and ℎ 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  was calculated by substituting the physical quantities  and  at the virtual boundary point into Equation (20) as follows: The virtual equilibrium distribution function  * (,  ) at point E was linearly extrapolated from the equilibrium distribution function  (,  ), calculated in Equation (25), and the equilibrium distribution function  (,  ) at point D from the internal ratio :  shown in Figure 6.The nonequilibrium part of the virtual distribution function  * (,  ) at point E was interpolated with the nonequilibrium part of the distribution function  (,  ) at point D as follows: Pressure p vb 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 vb 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: 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 eq α was calculated by substituting the physical quantities u vb and p vb at the virtual boundary point into Equation (20) as follows: The virtual equilibrium distribution function p eq * α (t, x E ) at point E was linearly extrap- olated from the equilibrium distribution function p eq α (t, x vb ), calculated in Equation (25), and the equilibrium distribution function p eq α (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 neq * α (t, x E ) at point E was interpolated with the nonequilibrium part of the distribution function p neq α (t, x D ) at point D as follows: 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 eq * α (t, x E ) and the nonequi- librium part of the distribution function p neq * α (t, x E ) at point E as follows:

Computational Model
Figure 7 shows a schematic view of the computational model used in this study.The computational domain has a size of 40L × 40L × 40L, 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 tip .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 = 40L, whereas the Neumann boundary condition was used in other planes for pressure and in all planes for velocity.

Computational Model
Figure 7 shows a schematic view of the computational model used in this study.The computational domain has a size of 40 40 40, where the representative length  is the mean chord length  and the representative speed  is the mean velocity of the wing tip  .Note that a grid model consisting of four blocks of different grid sizes (fourtiered 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  was set to 1 3 ⁄ and velocity  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 - plane of  = 40, whereas the Neumann boundary condition was used in other planes for pressure and in all planes for velocity.

Evaluation Parameters
We calculated the lift coefficient  , thrust coefficient  , and power coefficient  as the evaluation indices of the fluid in the flapping flight, as described in Section

Evaluation Parameters
We calculated the lift coefficient C L , thrust coefficient C T , and power coefficient C PWR 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: where F y and F z are the forces in the y and z directions, f local is the force, and u local 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: where u is the velocity vector and ω is the vorticity vector.Moreover, Ω ij is the vorticity tensor and S ij is the deformation velocity tensor, which are expressed as follows: where D is the velocity gradient tensor, which is calculated as follows: In this study, the Q value and helicity density h d are nondimensionalized as follows:

Grid Independence and Flapping Cycle Convergence Tests
A grid independence test was performed for three patterns with grid sizes of 16, 32, and 64 cells/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.
Biomimetics 2024, 9, x FOR PEER REVIEW movement; however, it converged immediately, and little change was found in t history.In addition, not much difference was observed after the fifth cycle (Table 3) fore, the analysis conducted in this study was performed using the data from th cycle.For the grid independence test, Figure 9 shows the time history of the lift coe  , the pressure component of the lift coefficient  , and the viscous stress com the of lift coefficient  over one cycle at each resolution, and Table 4 shows th averaged values of them  ,  , and  .Based on these results, the difference i each resolution decreased as the resolution increased, and not much difference w served between 32 and 64 cells/.Therefore, a resolution of 32 cells/ was sele ter considering the computational cost and accuracy.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 cycleaveraged 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 cells/L.Therefore, a resolution of 32 cells/L was selected after considering the computational cost and accuracy.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): 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.

Example Test for Analyzing the Three-Dimensional Flapping Motion
Fluid forces generated by a rectangular rigid plate in simple motion were compar with those in previous studies to verify the physical validity of the analysis results fo moving rigid plate in a three-dimensional analysis.The motion is oscillating, which expressed through the following translational displacement (heaving motion) and ro tion (flapping motion): where () is the displacement of the -coordinate of the center of gravity and () the oscillatory angle.Figure 10 shows a schematic view of the movement of an oscillati plate.Figure 11 shows the time history of lift coefficient generated from an oscillatory pla As shown in Figure 11, the simulated values were similar to the results of other stud 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.[36,37].In addition, the cycle-averaged value of the lift was 0.223 in this study, w value was 0.22 in a previous study [36], with a difference of 1.3% and the results analysis were nearly consistent with the previous study.

Effect of Figure-Eight Motion
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.
analysis were nearly consistent with the previous study.Figure 13 shows the isosurfaces of the vortex structure  * = 0.8 formed via each motion at  = 0.25, 0.50, 0.75, and 1.00.The isosurfaces were colored on the basis of the normalized helicity density ℎ * .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    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  , the thrust coefficient  , and the power coefficient  over one stroke cycle for each motion, and Table 5 shows the cycle-averaged values of coefficients  ,  , and  and the ratio of the lift coefficient to the power coefficient   ⁄ for each motion.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 PWR over one stroke cycle for each motion, and Table 5 shows the cycle-averaged values of coefficients C L , C T , and C PWR and the ratio of the lift coefficient to the power coefficient C L /C PWR for each motion.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  = 0.1, 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 ( * = 0.8) formed by each motion.The isosurfaces were colored on the basis of the normalized pressure coefficient  .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 -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  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.1T, 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 figureeight 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.

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 -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.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 cycleaveraged thrust coefficient in each was small (Table 5).Therefore, the effect of the figureeight 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.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.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 φ e0 = 0 • and φ e0 = 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 φ e0 = 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 Figure 16, the motions with  = 0° and  = 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  = 45° at any Reynolds number.Therefore, in this study, the figure-eight motion was defined as the motion with  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 ( ≤ 90°) had a higher ratio of the lift coefficient to the power coefficient, and the motion could generate lift more efficiently.The motion with  = 90° at  = 33.5 and  = 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  = 10° by increasing the vertical motion of the wings in each stroke.On the contrary, the motion with  = 90° had a larger vertical motion and a larger angle of attack compared with the motion with  = 45°; however, the upward motion was larger.Therefore, the motion was less efficient than the motion with  = 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  = 90° at  = 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   ⁄ 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.As shown in Table 6, the motion with a smaller initial phase of the elevation motion angle (φ e0 ≤ 90 • ) had a higher ratio of the lift coefficient to the power coefficient, and the motion could generate lift more efficiently.The motion with φ e0 = 90 • at Re = 33.5 and φ e0 = 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 φ e0 = 10 • by increasing the vertical motion of the wings in each stroke.On the contrary, the motion with φ e0 = 90 • had a larger vertical motion and a larger angle of attack compared with the motion with φ e0 = 45 • ; however, the upward motion was larger.Therefore, the motion was less efficient than the motion with φ e0 = 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 φ e0 = 90 • at Re = 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 PWR 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.As shown in Figure 16, the motions with  = 0° and  = 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  = 45° at any Reynolds number.Therefore, in this study, the figure-eight motion was defined as the motion with  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 ( ≤ 90°) had a higher ratio of the lift coefficient to the power coefficient, and the motion could generate lift more efficiently.The motion with  = 90° at  = 33.5 and 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 AR can be expressed as AR = l l /l s , where AR 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 θ el 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.
= 90° at other Reynolds numbers was the most efficient in generating lift.As sho in Figure 16, this result could be attributed to the fact that the angle of attack was lar than the original figure-eight motion with  = 10° by increasing the vertical motio the wings in each stroke.On the contrary, the motion with  = 90° had a larger vert motion and a larger angle of attack compared with the motion with  = 45°; howe the upward motion was larger.Therefore, the motion was less efficient than the mo with  = 45° when it was averaged over the entire stroke.However, horizontal fo dominate rather than vertical forces toward motion at low Reynolds numbers; theref increasing the vertical motion of the flapping motion and the angle of attack are impor for generating more lift.Therefore, the motion with  = 90° at  = 33.5 was the m efficient in generating lift.Moreover, the motion patterns of the most efficient in gene ing lift at each Reynolds number were investigated in detail.Figure 18 shows the rela between the initial phase of the elevation angle of the most efficient motion in genera lift and Reynolds number.Note that the most efficient motion in generating lift was culated from the tangential point where   ⁄ was the greatest slope on the ellipt approximation shown in Figure 17.As shown in Figure 18, the initial phase of the ele tion angle of the most efficient motion in generating lift increased as the Reynolds num decreased, whereas the angle decreased and approached zero as the Reynolds num increased.In addition, as shown in Table 6 and Figure 17f, all motions generated more lift as Reynolds number increased, and more lift was generated compared with the power.above discussion is based on nondimensionalized values without considering differen in the speed of the wings (flapping frequency) to simplify the comparison of the lift power generated by each motion.Then, the physical values of the lift and power co cients were summarized to investigate the effect of the flapping frequency on aero namic characteristics.Table 7 shows the cycle-averaged values of the force in the  di tion  , power  , and the ratio of the lift coefficient to the power coefficient     17, the eight types of figure-eight motions showe elliptical distribution in the relation between the lift coefficient and power coeffic Therefore, we discuss the relation between the lift coefficient and power coefficien making an elliptical approximation along the distribution.As shown in Figure 17a-e area surrounded by the ellipse and the W8 line (light blue shown in the figure) decre as the Reynolds number increased, and the number of motions that generate lift m efficiently than without the figure-eight motion also decreased.This decrease indi that the effect of figure-eight motion decreased as the Reynolds number increased elliptical shape, which is primarily related to the decrease in area, was focused on to cuss this difference in effect in more detail.Therefore, the effect of the change in the R olds number on the figure-eight motion was determined by calculating the aspect and inclination angle of the ellipse at each Reynolds number.As shown in Figure 19 long and short sides of the ellipse for each Reynolds number were defined as  an respectively, and the aspect ratio  can be expressed as  =   ⁄ , where  i ratio of the long side to the short side of the ellipse.In addition, the inclination ang the ellipse and the location of the center of the ellipse in Figure 17 were defined as  (, ).Table 8 shows the values of the elliptic approximation, and Figure 20 show 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 cha remarkably; the aspect ratio increased with increasing Reynolds number, and the as ratio converged when the Reynolds number increased above a certain value.The incr in the aspect ratio indicated that the elliptical approximation of various figure-eight tions approached the straight "without-8-line" shown in Figure 17, which indicates the difference between the motions with and without a figure-eight motion decreased addition, the increase in the Reynolds number was synonymous with the increase in in size.Therefore, varying the elevation motion angle that represents the figure-eight mo was not very effective for insects under high Reynolds number conditions (large s Particularly, the figure-eight motion may be a vital mechanism for insects to generate more efficiently under low Reynolds number conditions (small size).
In this study, insect wings were analyzed as rectangular rigid plates without th ness.Regarding the shape, Kirishna et al. [38] reported that flight efficiency did not cha much depending on the shape of the wings between a rectangular model with a mo based on the actual shape of a blowfly.However, the changes in the shape of the w may affect the efficiency of each figure-eight motion pattern.Therefore, in future stu of insects, it is necessary to evaluate figure-eight motion by considering not only the sh of the wings but also the flexibility of the wings and other parameters.

Conclusions
In this study, we investigated the more efficient figure-eight motion patterns in g erating lift during the hovering motion related to the Reynolds number.For this purp we compared the ratios of the cycle-averaged lift coefficient to the power coefficient g erated via each motion by varying the elevation motion angle and Reynolds number.C sequently, the motion with a smaller initial phase of the elevation motion angle ( 90°) could generate lift more efficiently at all Reynolds numbers.In addition, the in phase of the elevation angle of the most efficient motion in generating lift increased as Reynolds number decreased, whereas the angle decreased and approached  = 0 the Reynolds number increased.Moreover, varying the elevation motion angle that resents the figure-eight motion was effective for insects under low 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.

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 (φ e0 ≤ 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 φ e0 = 0 as the Reynolds number increased.Moreover, varying the elevation motion angle that represents the figureeight motion was effective for insects under low Reynolds number conditions and was not very effective for insects under high Reynolds number conditions.Therefore, the motion may be a vital mechanism for insects under low Reynolds number conditions to generate lift more efficiently.

Figure 1 .
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 .
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 2024, 9 ,
x FOR PEER REVIEW 4 of 26 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  , 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 <  < 180 and 180 <  <360), except for the three line-shaped ( = 0° ) and U-shaped ( = 0° and  = 180° ) motions.Therefore, in summarizing all the motions, we defined the motion with  = 0° as a motion "without figure-eight motion" and the motion with  0° as a motion "with figure-eight mo-

Figure 2 .
Figure 2. Definition of the flapping motion of the insect in (a) bird's eye view, (b) - plane, and (c) - 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 .
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 3 .
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 3(b),The gray solid and dotted lines represent the trajectory of the wings tips and the body of the insect, respectively.

Figure 3 .
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 4 .
Figure 4. 3D27V model for the three-dimensional lattice Boltzmann method.The numbers (0~26) represent the 27 directions of the 3D27Vmodel.

Figure 4 .
Figure 4. 3D27V model for the three-dimensional lattice Boltzmann method.The numbers (0~26) represent the 27 directions of the 3D27Vmodel.

Figure 5 .
Figure 5. Virtual boundary points.The gray and white areas represent the interior of the object and fluid, respectively.

Figure 5 .
Figure 5. Virtual boundary points.The gray and white areas represent the interior of the object and fluid, respectively.

Figure 6 .
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 .
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 7 .
Figure 7. Schematic view of the computational model: (a) bird's eye view, (b) - plane, and (c) - plane.The representative length  is the mean chord length  shown in Section 2.1.In Figure 7(a), the wings and grid model are shown in green and black or gray lines, respectively.In Figure 7(b) and (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 .
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 8 .
Figure 8.Time history of the lift coefficient in seven cycles with  = 0.04.

Figure 8 .
Figure 8.Time history of the lift coefficient in seven cycles with U = 0.04.

Figure 9 .
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 .
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 .
Figure 9.Time history of the (a) lift coefficient, (b) pressure component of lift coefficient, and viscous stress component of lift coefficient for three resolutions.

Figure 10 .
Figure 10.Schematic view of the movement of an oscillating plate.

Figure 10 .
Figure 10.Schematic view of the movement of an oscillating plate.

Figure 11 .
Figure 11.Time history of the lift coefficient of an oscillating plate.The result was compare those of Trizila[36] and Wang et al.[37].

A
comparison of the vortex structure and fluid forces generated by each moti conducted to investigate the effect of the figure-eight motion on the aerodynamic motion was represented by varying the elevation motion angle  which is a ro motion that represents the figure-eight motion shown in Equation (4).Motion with eight motion was represented by the elevation motion angle with the parameters in Section 2.2, and motion without figure-eight motion was represented by  = 0 tions (41) and (42) show the angles with and without figure-eight motions as follo  = 70 (2(2) − 10) + 10 ,  = 0.

Figure 12
Figure 12 shows a schematic view of each motion represented by the flapping angle  and the feathering motion angle  shown in Section 2.2 in addition to vation motion angle.

Figure 11 .
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].

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:

Figure 11 .
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].

AFigure 12
Figure 12 shows a schematic view of each motion represented by the flapping motion angle  and the feathering motion angle  shown in Section 2.2 in addition to the elevation motion angle.

Figure 12 .
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 .
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 13
Figure13shows the isosurfaces of the vortex structure Q * = 0.8 formed via each motion at t = 0.25T, 0.50T, 0.75T, and 1.00T.The isosurfaces were colored on the basis of the normalized helicity density h * d .As shown in Figure13, 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.

Figure 13 .
Figure 13.Vortex structures and normalized helicity density over one stroke cycle for each motion.(a) With and (b) without figure-eight motions at  = 0.25; (c) with and (d) without figure-eight motions at  = 0.50; (e) with and (f) without figure-eight motions at  = 0.75; and (g) with and (h) without figure-eight motions at  = 1.00.

Figure 13 .
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.25T; (c) with and (d) without figure-eight motions at t = 0.50T; (e) with and (f) without figure-eight motions at t = 0.75T; and (g) with and (h) without figure-eight motions at t = 1.00T.

Figure 14 .
Figure 14.Time histories of (a) lift coefficient, (b) thrust coefficient, and (c) power coefficient over one stroke cycle for each motion.

Figure 14 .
Figure 14.Time histories of (a) lift coefficient, (b) thrust coefficient, and (c) power coefficient over one stroke cycle for each motion.

Figure 15 .Figure 15 .
Figure 15.Vortex structures ( * = 0.8) and normalized pressure coefficient for each motion in the downstroke at  = 0.1.Vortex structures ( * = 0.8) and normalized pressure coefficients from (a) Figure 15.Vortex structures (Q * = 0.8) and normalized pressure coefficient for each motion in the downstroke at t = 0.1T.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.

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 figureeight motion and Reynolds number.These motions were represented as eight types of figure-eight motions by varying the initial phase φ e0 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.Biomimetics 2024, 9, x FOR PEER REVIEW 19 of 26 with and (b) without figure-eight motions.Normalized pressure coefficient at the wing tip from (c) with and (d) without figure-eight motions.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  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.(a)  = 0°.(b)  = 45°.(c)  = 90°.(d)  = 135°.(e)  = 180°.(f)  = 225°.(g)  = 270°.(h)  = 315°.

Figure 18 .
Figure 18.Relation between the initial phase of elevation angle of the most efficient motion in generating lift and Reynolds number.

imetics 2024, 9 ,
x FOR PEER REVIEW 21 of 26 546; (f) all Reynolds numbers.Light blue areas shown in the figure ( ) indicate more efficiency than without figure-eight motions.

)
indicate more efficiency than without figure-eight motions.

Figure 18 .
Figure 18.Relation between the initial phase of elevation angle of the most efficient motion in erating lift and Reynolds number.

Figure 18 .
Figure 18.Relation between the initial phase of elevation angle of the most efficient motion in generating lift and 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.Re [−] F z [µN] P [µW] F z /P [−]

Figure 19 .
Figure 19.Schematic view of the elliptic approximation.

Figure 19 .
Figure 19.Schematic view of the elliptic approximation.

Figure 20 .
Figure 20.Relation between the aspect ratio and Reynolds number.

Figure 20 .
Figure 20.Relation between the aspect ratio and Reynolds number.

Table 1 .
Physical properties of an insect.

Table 1 .
Physical properties of an insect.

Table 2 .
Advection parameters of the particles.

Table 3 .
Cycle-averaged lift coefficient in seven cycles.

Table 4 .
Cycle -averaged lift coefficient, pressure component of lift coefficient, and viscous stress component of lift coefficient for three resolutions.

Table 3 .
Cycle-averaged lift coefficient in seven cycles.

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 str component of lift coefficient for three resolutions.

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 figureeight motion.

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 ( ) 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 ( (f) all Reynolds numbers.Light blue areas shown in the figure ( ) indicate more efficiency than without figure-eight motions.

Table 7
shows the cycle-averaged values of the force in the  direction  , power  , and the ratio of the lift coefficient to the power coefficient   ⁄

Table 7 .
Cycle-averaged values of the force in the  direction  , power , and the ratio of the in the  direction to the power   ⁄ without figure-eight motion at each Reynolds number.

Table 8 .
Respective parameters of the elliptic approximation for each Reynolds number.