*Article* **Unsteady Aerodynamic Design of a Flapping Wing Combined with a Bionic Wavy Leading Edge**

**Xuan Bai, Hao Zhan and Baigang Mi \***

**\*** Correspondence: mibaigang@163.com

**Abstract:** Based on the bionic design of the humpback whale fin, a passive flow control method is proposed to obtain greater flapping lift by applying the wavy leading edge structure to the straight symmetrical flapping wing. The leading edge of the conventional flapping wing is replaced by the wavy shape represented by regular trigonometric function to form a special passive flow control configuration imitating the leading edge of the humpback whale fin. The dynamic aerodynamic performance and flow field characteristics of straight wing and wavy leading edge flapping wing with different parameters are compared and analyzed by CFD numerical simulation. The simulation results show that the wavy leading edge structure changes the flow field of the baseline flapping wing and reduces the pressure on the upper surface of the flapping wing during the process of downward flapping, thereby increasing the pressure difference between the upper and lower surfaces of the flapping wing and increasing the lift. The sensitivity analysis of the design parameters shows that in order to obtain the maximum lift coefficient while losing the least thrust, the smaller amplitude should be selected on the premise of selecting the smaller wavelength. Among the configurations of different design parameters calculated in this paper, the optimal wavy leading edge flapping wing configuration increases the time average lift coefficient by 32.86% and decreases the time average thrust coefficient by 14.28%. Compared with the straight wing, it has better low-speed flight and can withstand greater take-off weight.

**Keywords:** flapping wing; wavy leading edge; flow control; bionics; computational fluid dynamics (CFD)

### **1. Introduction**

Flapping wing aircraft generates lift to overcome gravity and thrust to fly forward by imitating birds flapping and twisting wings. This flight mode has the unique characteristics such as small size, high maneuverability, good flexibility, and strong concealment. It has shown broad development prospects in military, civilian, detection, earthquake relief and other industries. Therefore, flapping wing aircraft has become an important direction and hot issue in the research and development of UAV aircraft.

Because the geometric size of the flapping wing is relatively small, the chord length is usually less than 15 cm, and it is often in the low Reynolds number range (104–10<sup>6</sup> ), the flapping wing flight principle is complex. In addition, it exhibits a highly unsteady nonlinear flow phenomenon. In order to realize the efficient and stable flight of flapping wing aircraft, it is necessary to study on its aerodynamic characteristics and influencing factors based on the biological knowledge of bird flight. The main influencing factors include the static geometric shape of the flapping wing [1,2], multi-degree of freedom flapping [3], dynamic flexible deformation [4], wing tip slit [5,6], leading/trailing edge serrations [7,8], leading edge alula [9,10], and other small-scale flow control structures. In the process of the preliminary design and flow control of the flapping wing aircraft, it is very important and effective to improve its flow field and aerodynamic performance by referring to the biological structure of nature.

**Citation:** Bai, X.; Zhan, H.; Mi, B. Unsteady Aerodynamic Design of a Flapping Wing Combined with a Bionic Wavy Leading Edge. *Appl. Sci.* **2023**, *13*, 1519. https://doi.org/ 10.3390/app13031519

Academic Editors: Josep Maria Bergadà and Gabriel Bugeda Castelltort

Received: 27 December 2022 Revised: 18 January 2023 Accepted: 19 January 2023 Published: 24 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In 1995, Fish et al. [11] observed the cross section of the whale fins and studied the hydrodynamic performance of the convex structure of the fin limb of the humpback whale for the first time. They pointed out that the wavy leading edge can play a role like the lifting device. Since then, the leading edge wavy configuration has been tried to apply to the study of flow control under various scenarios and conditions. Sudhakar et al. [12] took a typical UAV as the research object and pointed out that compared with the straight wing, the wavy leading edge wing can increase lift and reduce drag and the lift-drag ratio increased by 25%. In addition, by using the oil flow method and PIV technique, they found that the wavy leading edge also makes the height and length of the laminar separation bubble on the straight wing change significantly under the condition of low angle of attack (α = 6◦ ) [13]. Wei et al. [14] studied the effect of wavy leading edge on the aerodynamic performance of swept conical SD7032 wing by experiment. The results show that the wavy leading edge can slightly improve the stall performance, and the stall angle of attack can be delayed by 2–4 degrees. However, when the angle of attack is in the range of 7–20 degrees, it makes the lift drag ratio of the original wing decrease by 0–40%. Abdelrahman et al. [15] adopted an infinite span rectangular wing with a high radian S1223 airfoil to numerically study the influence of leading edge wavys on its performance and concluded that the wavy structure can completely change the flow structure on the airfoil. The flow separation on the wing surface is limited to the trough of the wavy leading edge, which plays a role in flow control. Seyhan et al. [16] verified the feasibility of flow control of NACA0015 airfoil by wavy leading edge under low Reynolds number conditions. Gopinathan et al. [17] compared the effect of wavy leading edge structure on swept wings of two different airfoils by numerical and experimental methods, and pointed out that this method has great development potential. In addition to extensive researches on wings [12–20], as a flow control method, wavy leading edge has also been tried to apply in the design of rotor [21], hydraulic turbine cascade and steam turbine cascade [22–24], wind turbine blade [25–28], hydrofoil [29–32], etc., to improve mechanical efficiency. However, their motion mode and flow field characteristics are far different from those of flapping wings. Flapping wings have a complex flight principle and highly unsteady nonlinear flow phenomenon, which brings difficulty to flow control. Anwar et al. [33] numerically simulated the aerodynamic performance of six wavy leading edge flapping wing with different design parameters similar to insect flapping motion. The results show that the lift-to-drag ratio of the best performed wavy leading edge wing decreases by 0.34%, and the efficiency decreases by 0.35%. It concluded that the wavy leading edge wing is not better than the straight wing in insect flapping motion with azimuth rotation. However, the flapping mode of imitation bird flapping wings is different from that of insects, and the influence of wavy leading edge on its aerodynamic performance remains to be further studied. The studies on the bionic wavy leading edge are summarized as shown in Table 1.

In this paper, based on bionics, a passive flow control method is proposed to solve the problem that the average lift force of the straight symmetrical flapping wing is difficult to improve. The passive flow control method is applied to the flapping wing by applying the wavy leading edge which imitating the fin structure of the humpback whale. The aerodynamic performance parameters of the straight wing and the wavy leading edge flapping wing with different design parameters are calculated and compared, and the flow field characteristics and control mechanism are analyzed. At the same time, the sensitivity of the design parameters is analyzed, and the basis for selecting the optimal design parameters of the wavy leading edge flapping wing configuration is given. It brings new inspiration for flapping wing to obtain flapping lift and unsteady flow control under a low Reynolds number.


**Table 1.** Studies on flow control using bionic wavy leading edge.

### **2. Aerodynamic Characteristics Analysis of Baseline Flapping Wing**

### *2.1. Flapping Wing Motion Model*

The flapping wing motion can be decomposed into two basic motion forms: flapping and pitching, as shown in Figure 1. The cross-sections of bird's wing like thin and cambered airfoil. NACA0014 airfoil has been experimented [34] and numerically simulated [34,35] as a flapping wing. Therefore, the flapping wing model adopted in this paper is a rectangular straight wing with a finite span and the wing section of NACA0014 airfoil. The chord length of the wing is *c*, the span length is *L* and aspect ratio is *AR* = *<sup>L</sup> c* . The set flapping wing motion law is expressed as Equation (1) [35], where, *f* is the flapping frequency, *ψ*(*t*) represents the flapping motion, *ψ*<sup>1</sup> is the amplitude of the flapping motion, *α*(*t*) represents the pitching motion, *α*<sup>0</sup> is the pre-installation angle of the wing, *α*<sup>1</sup> is the amplitude of the pitching motion, *φ* is the phase difference, and the center of the pitching motion is

the quarter chord length of each section of the wing. Figure 2 shows the time history plot corresponding to Equation (1).

$$\begin{cases} \quad \psi(t) = \psi\_1 \cos(2\pi ft) \\\ a(t) = a\_0 + a\_1 \cos(2\pi ft + \phi) \end{cases} \tag{1}$$

center of the pitching motion is the quarter chord length of each section of the wing. Figure

is the phase difference, and the

 )

 1

0

 

 *ft*  )  )

 

(1)

(1)

 

 +

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 4 of 24

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 4 of 24

is the amplitude of the pitching motion,

2 shows the time history plot corresponding to Equation (1).

2 shows the time history plot corresponding to Equation (1).

**Figure 1.** Motion of flapping wing.**Figure 1.** Motion of flapping wing.

wing,

wing,

1

1

**Figure 2.** The time history plot corresponding to Equation (1). **Figure 2.** The time history plot corresponding to Equation (1).

the approaching velocity

, the reduced frequency is

, the reduced frequency is

 *m*

length

in Table 2.

in Table 2.

*c*

= 0.12

**Figure 2.** The time history plot corresponding to Equation (1). According to the flight speed and flapping frequency of small birds such as pigeons, According to the flight speed and flapping frequency of small birds such as pigeons, the approaching velocity *U* is 15 m/s, the flapping frequency *f* is 8 *Hz*, the chord length *c m* = 0.12 , the aspect ratio *AR* = 8, the Reynolds number is Re*U*=According to the flight speed and flapping frequency of small birds such as pigeons, the approaching velocity *U*<sup>∞</sup> is 15 m/s, the flapping frequency *f* is 8 Hz, the chord length *c* = 0.12m, the aspect ratio *AR* = 8, the Reynolds number is Re = *ρU*∞*c <sup>µ</sup>* <sup>=</sup> 1.2 <sup>×</sup> <sup>10</sup><sup>5</sup> , the reduced frequency is *k* = *π f c U*∞ = 0.2. The specific values of each parameter are shown in Table 2.

, the aspect ratio *AR* = 8, the Reynolds number is

**Parameter Value Parameter Value**

is 15 m/s, the flapping frequency

**Parameter Value Parameter Value**

*f*

Re

. The specific values of each parameter are shown

*U c*

=

. The specific values of each parameter are shown

is 8 *Hz*, the chord

 1.2 10

 

 *c*

 = 5

 1.2 10

 

 = 5

**Table 2.** Relevant parameters of flapping wing.

*U*

0.2 *fc <sup>k</sup> U* 

0.2 *fc <sup>k</sup> U* 

 = =

 = =


**Table 2.** Relevant parameters of flapping wing.

### *2.2. Calculation Method and Verification*

Flapping wing has flight characteristics of low speed and low Reynolds number. The numerical calculation in this paper is based on the three-dimensional unsteady Reynoldsaveraged N-S Equation, which based on the finite volume method. The second-order up-wind scheme is used in the convection direction, the central difference is used in the diffusion term, and the Coupled algorithm is used to solve [36].

The two-Equation shear stress (*k* − *ω* SST) turbulence model [37] modified by the low Reynolds number is used for calculation. The turbulence model simulates the transition process with the damping function in the low Reynolds number eddy viscosity model, which has wide applicability. The *k* − *ω* SST turbulence model Equation can be expressed as the following Equation (2):

$$\begin{array}{lcl} \mu\_{l} &= \frac{\rho \mathbf{a}\_{l} \mathbf{k}}{\max(\mu\_{l} \omega \epsilon \Omega \mathbf{F}\_{2})}\\ \frac{\partial(\rho \mathbf{k})}{\partial t} + u\_{l} \frac{\partial(\rho \mathbf{k})}{\partial \mathbf{x}\_{i}} &= P\_{\mathbf{k}} - \beta\_{k} \rho \mathbf{k} \omega + \frac{\partial}{\partial \mathbf{x}\_{i}} \left[ \left( \mu\_{l} + \frac{\mu\_{l}}{\sigma\_{\mathbf{k}}} \right) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{i}} \right] \\ \frac{\partial(\rho \omega)}{\partial t} + u\_{l} \frac{\partial(\rho \omega)}{\partial \mathbf{x}\_{i}} &= C\_{\omega} P\_{\omega} - \beta\_{k} \rho \omega^{2} + \frac{\partial}{\partial \mathbf{x}\_{i}} \left[ \left( \mu\_{l} + \frac{\mu\_{l}}{\sigma\_{\omega}} \right) \frac{\partial \omega}{\partial \mathbf{x}\_{i}} \right] + 2\rho (1 - F\_{1}) \frac{1}{\sigma\_{\omega 2}} \frac{1}{\omega} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{i}} \frac{\partial \omega}{\partial \mathbf{x}\_{i}} \end{array} \tag{2}$$

The low Reynolds number-corrected *k* − *ω* SST turbulence model constructs the mode coefficient *α* ∗ in the eddy viscosity coefficient *µ<sup>t</sup>* = *α* <sup>∗</sup>*ρκ*/*ω* of the *k* − *ω* Equation as Equation (3) [38], where Re*<sup>t</sup>* = *ρκ*/*µω* represents the turbulent Reynolds number:

$$\alpha^\* = \alpha^\*\_{\infty} \left[ \frac{\alpha^\*\_0 + \frac{R\varepsilon\_t}{R\_k}}{1 + \frac{R\varepsilon\_t}{R\_k}} \right] \tag{3}$$

In order to verify the reliability of the numerical simulation, the rectangular straight wing of NACA0014 airfoil with aspect ratio of *AR* = 8 is selected. The mode of motion is sinusoidal oscillation, as shown in Equation (4). *z*(*t*) is the wing displacement in the vertical direction at time *t*, amplitude *h* = 0.4*c*, reduction frequency *k* = 0.2, and Reynolds number *Re* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>6</sup> [34].

$$z(t) = h \cos(kt)\tag{4}$$

The boundary conditions of the calculation model are referred to [34]. Since the overset grid method is used in the subsequent calculation, the surrounding and far field of the wing are divided into two parts during grid division. Figures 3 and 4 show the comparison of the average lift coefficient and the average thrust coefficient of the flapping wing in one period calculated by using the overset grid method, and the results calculated by Jones [34]. The two have good consistency, indicating that the numerical simulation method used in this paper is reliable.

method used in this paper is reliable.

method used in this paper is reliable.

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 6 of 24

in one period calculated by using the overset grid method, and the results calculated by Jones [34]. The two have good consistency, indicating that the numerical simulation

in one period calculated by using the overset grid method, and the results calculated by Jones [34]. The two have good consistency, indicating that the numerical simulation

**Figure 3.** Comparison of average lift coefficient. **Figure 3.** Comparison of average lift coefficient. **Figure 3.** Comparison of average lift coefficient.

**Figure 4.** Comparison of average thrust coefficient. **Figure 4.** Comparison of average thrust coefficient. **Figure 4.** Comparison of average thrust coefficient.

### *2.3. Grid and Time Step Independence Verification*

line case.

line case.

*2.3. Grid and Time Step Independence Verification* Overset grid method is used in this paper, where the component grid is fluid domain grid including flapping wing. The boundary layer is set on the surface of the flapping wing to ensure y+ < 1. The background grid is a cylindrical area with a radius of 25*c* and *2.3. Grid and Time Step Independence Verification* Overset grid method is used in this paper, where the component grid is fluid domain grid including flapping wing. The boundary layer is set on the surface of the flapping wing to ensure y+ < 1. The background grid is a cylindrical area with a radius of 25*c* and a height of 10*c*, and it is locally refined to ensure that the mesh size at the junction of the Overset grid method is used in this paper, where the component grid is fluid domain grid including flapping wing. The boundary layer is set on the surface of the flapping wing to ensure y+ < 1. The background grid is a cylindrical area with a radius of 25*c* and a height of 10*c*, and it is locally refined to ensure that the mesh size at the junction of the background grid and the component grid is close. The component and background grid are divided by polyhedral unstructured grids. The computational domain grids and boundary conditions are shown in Figure 5. And Figure 6 shows the y+ contours of baseline case.

(**a**) (**b**)

(**a**) (**b**)

a height of 10*c*, and it is locally refined to ensure that the mesh size at the junction of the background grid and the component grid is close. The component and background grid are divided by polyhedral unstructured grids. The computational domain grids and boundary conditions are shown in Figure 5. And Figure 6 shows the y+ contours of base-

background grid and the component grid is close. The component and background grid are divided by polyhedral unstructured grids. The computational domain grids and boundary conditions are shown in Figure 5. And Figure 6 shows the y+ contours of base-

**Figure 5.** Grid and boundary conditions. (**a**) Background grid; (**b**) Component grid; (**c**) Boundary layer and surface mesh. **Figure 5.** Grid and boundary conditions. (**a**) Background grid; (**b**) Component grid; (**c**) Boundary layer and surface mesh. **Figure 5.** Grid and boundary conditions. (**a**) Background grid; (**b**) Component grid; (**c**) Boundary layer and surface mesh.

in one period calculated by using the overset grid method, and the results calculated by Jones [34]. The two have good consistency, indicating that the numerical simulation

Overset grid method is used in this paper, where the component grid is fluid domain grid including flapping wing. The boundary layer is set on the surface of the flapping wing to ensure y+ < 1. The background grid is a cylindrical area with a radius of 25*c* and a height of 10*c*, and it is locally refined to ensure that the mesh size at the junction of the background grid and the component grid is close. The component and background grid are divided by polyhedral unstructured grids. The computational domain grids and boundary conditions are shown in Figure 5. And Figure 6 shows the y+ contours of base-

method used in this paper is reliable.

**Figure 3.** Comparison of average lift coefficient.

**Figure 4.** Comparison of average thrust coefficient.

*2.3. Grid and Time Step Independence Verification*

line case.

**Figure 6.** The y+ contours of baseline case. **Figure 6.** The y+ contours of baseline case. **Figure 6.** The y+ contours of baseline case.

In order to ensure the accuracy of numerical calculation and save computational cost, the independence of computational mesh and time step is verified by the straight flapping wing model. The fifth flapping cycle is selected for comparison, and the calculation results are stable. The lift and trust coefficient time history of the fifth and sixth cycles is compared in Figure 7. Firstly, three groups of overset grids with different number of grids are selected as shown in Table 3 to calculate. The calculation results are shown in Figure 8. As can be seen, the grid solution results of Case2 and Case3 tend to coincide, which can be considered that the grid has reached the grid independence requirements. Therefore, the grid division scheme of Case2 is selected for subsequent calculation. Secondly, to verify the independence of the time step, three different time steps for calculation are taken; the unit is s, and the result is shown in Figure 9. When the time step is 0.0025 s and 0.001 s, the relative error between the time average lift coefficient, and the time average trust coefficient is less than 2%, which meets the requirement of independence. Therefore, in the subsequent calculation, the time step is set to 0.0025 s. In order to ensure the accuracy of numerical calculation and save computational cost, the independence of computational mesh and time step is verified by the straight flapping wing model. The fifth flapping cycle is selected for comparison, and the calculation results are stable. The lift and trust coefficient time history of the fifth and sixth cycles is compared in Figure 7. Firstly, three groups of overset grids with different number of grids are selected as shown in Table 3 to calculate. The calculation results are shown in Figure 8. As can be seen, the grid solution results of Case2 and Case3 tend to coincide, which can be considered that the grid has reached the grid independence requirements. Therefore, the grid division scheme of Case2 is selected for subsequent calculation. Secondly, to verify the independence of the time step, three different time steps for calculation are taken; the unit is s, and the result is shown in Figure 9. When the time step is 0.0025 s and 0.001 s, the relative error between the time average lift coefficient, and the time average trust coefficient is less than 2%, which meets the requirement of independence. Therefore, in the subsequent calculation, the time step is set to 0.0025 s. In order to ensure the accuracy of numerical calculation and save computational cost, the independence of computational mesh and time step is verified by the straight flapping wing model. The fifth flapping cycle is selected for comparison, and the calculation results are stable. The lift and trust coefficient time history of the fifth and sixth cycles is compared in Figure 7. Firstly, three groups of overset grids with different number of grids are selected as shown in Table 3 to calculate. The calculation results are shown in Figure 8. As can be seen, the grid solution results of Case2 and Case3 tend to coincide, which can be considered that the grid has reached the grid independence requirements. Therefore, the grid division scheme of Case2 is selected for subsequent calculation. Secondly, to verify the independence of the time step, three different time steps for calculation are taken; the unit is s, and the result is shown in Figure 9. When the time step is 0.0025 s and 0.001 s, the relative error between the time average lift coefficient, and the time average trust coefficient is less than 2%, which meets the requirement of independence. Therefore, in the subsequent calculation, the time step is set to 0.0025 s.

(**a**) (**b**)

(**a**) (**b**)

**Figure 7.** The comparison of time history of lift and thrust coefficient. (**a**) Lift coefficient; (**b**)

**Figure 7.** The comparison of time history of lift and thrust coefficient. (**a**) Lift coefficient; (**b**)

Trust coefficient.

Trust coefficient.

(**c**)

Boundary layer and surface mesh.

**Figure 6.** The y+ contours of baseline case.

subsequent calculation, the time step is set to 0.0025 s.

**Figure 5.** Grid and boundary conditions. (**a**) Background grid; (**b**) Component grid; (**c**)

In order to ensure the accuracy of numerical calculation and save computational cost, the independence of computational mesh and time step is verified by the straight flapping wing model. The fifth flapping cycle is selected for comparison, and the calculation results are stable. The lift and trust coefficient time history of the fifth and sixth cycles is compared in Figure 7. Firstly, three groups of overset grids with different number of grids are selected as shown in Table 3 to calculate. The calculation results are shown in Figure 8. As can be seen, the grid solution results of Case2 and Case3 tend to coincide, which can be considered that the grid has reached the grid independence requirements. Therefore, the grid division scheme of Case2 is selected for subsequent calculation. Secondly, to verify the independence of the time step, three different time steps for calculation are taken; the unit is s, and the result is shown in Figure 9. When the time step is 0.0025 s and 0.001 s, the relative error between the time average lift coefficient, and the time average trust coefficient is less than 2%, which meets the requirement of independence. Therefore, in the

**Figure 7.** The comparison of time history of lift and thrust coefficient. (**a**) Lift coefficient; (**b**) Trust coefficient. **Figure 7.** The comparison of time history of lift and thrust coefficient. (**a**) Lift coefficient; (**b**) Trust coefficient. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 8 of 24 *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 8 of 24



**Figure 8.** Grid independence verification. (**a**) Lift coefficient; (**b**) Trust coefficient. **Figure 8.** Grid independence verification. (**a**) Lift coefficient; (**b**) Trust coefficient. **Figure 8.** Grid independence verification. (**a**) Lift coefficient; (**b**) Trust coefficient.

**Figure 9.** Time step independence verification. (**a**) Lift coefficient; (**b**) Trust coefficient. **Figure 9.** Time step independence verification. (**a**) Lift coefficient; (**b**) Trust coefficient. leading edge vortex and the wing tip vortex begin to diffuse, the intensity is weakened. **Figure 9.** Time step independence verification. (**a**) Lift coefficient; (**b**) Trust coefficient.

lift during the upstroke. Many studies have shown that for this type of flapping wing, the flapping frequency increases can increase the peak value of positive and negative lift and increase the thrust, but the average lift in the cycle cannot increase much. The time-averaged lift coefficient is small and difficult to improve, and therefore, the takeoff weight of aircraft with straight symmetrical flapping wing is generally small, and the flight perfor-

wing. As shown in Figure 9, it produces positive lift during the downstroke, and negative lift during the upstroke. Many studies have shown that for this type of flapping wing, the flapping frequency increases can increase the peak value of positive and negative lift and increase the thrust, but the average lift in the cycle cannot increase much. The time-averaged lift coefficient is small and difficult to improve, and therefore, the takeoff weight of aircraft with straight symmetrical flapping wing is generally small, and the flight perfor-

The flapping wing flutters downward from the highest position (t/T = 0), forming the leading edge vortex at the leading edge of the upper surface and the wing tip vortex at the wing tip. When the flapping wing flaps down to the horizontal position (t/T = 0.25), the vortex intensity of leading edge vortex and tip vortex reaches the maximum, as shown in Figure 10; the left side is the wing tip direction. The inner part of the red circle is the vortex. As the flapping continues downward, the flapping wing begins to decelerate, the leading edge vortex and the wing tip vortex begin to diffuse, the intensity is weakened.

The flapping wing flutters downward from the highest position (t/T = 0), forming the leading edge vortex at the leading edge of the upper surface and the wing tip vortex at the wing tip. When the flapping wing flaps down to the horizontal position (t/T = 0.25), the vortex intensity of leading edge vortex and tip vortex reaches the maximum, as shown in Figure 10; the left side is the wing tip direction. The inner part of the red circle is the vortex. As the flapping continues downward, the flapping wing begins to decelerate, the

*2.4. Flow Field Structure Characteristics of Baseline Flapping Wing*

*2.4. Flow Field Structure Characteristics of Baseline Flapping Wing*

mance is hard to improve [39,40].

mance is hard to improve [39,40].

### *2.4. Flow Field Structure Characteristics of Baseline Flapping Wing*

The baseline flapping wing used in this paper is a straight symmetrical flapping wing. As shown in Figure 9, it produces positive lift during the downstroke, and negative lift during the upstroke. Many studies have shown that for this type of flapping wing, the flapping frequency increases can increase the peak value of positive and negative lift and increase the thrust, but the average lift in the cycle cannot increase much. The time-averaged lift coefficient is small and difficult to improve, and therefore, the takeoff weight of aircraft with straight symmetrical flapping wing is generally small, and the flight performance is hard to improve [39,40].

The flapping wing flutters downward from the highest position (t/T = 0), forming the leading edge vortex at the leading edge of the upper surface and the wing tip vortex at the wing tip. When the flapping wing flaps down to the horizontal position (t/T = 0.25), the vortex intensity of leading edge vortex and tip vortex reaches the maximum, as shown in Figure 10; the left side is the wing tip direction. The inner part of the red circle is the vortex. As the flapping continues downward, the flapping wing begins to decelerate, the leading edge vortex and the wing tip vortex begin to diffuse, the intensity is weakened. The asymmetry of the up and down flapping process and the strong leading edge vortex attached to the upper surface during the down flapping process are the important mechanisms of lift generation. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 9 of 24 The asymmetry of the up and down flapping process and the strong leading edge vortex attached to the upper surface during the down flapping process are the important mechanisms of lift generation.

**Figure 10.** Velocity vector diagram of 1% chordwise section at t/T = 0.25. **Figure 10.** Velocity vector diagram of 1% chordwise section at t/T = 0.25.

### **3. Design of Wavy Leading Edge Flapping Wing 3. Design of Wavy Leading Edge Flapping Wing**

### *3.1. Flow Control Mechanism of Wavy Leading Edge 3.1. Flow Control Mechanism of Wavy Leading Edge*

Figure 11 shows that the leading edge of the flippers of the humpback whale has special nodular structures. These wavy bulges can maintain lift at a high angle of attack and provide higher maneuverability during the predation process [11]. Figure 11 shows that the leading edge of the flippers of the humpback whale has special nodular structures. These wavy bulges can maintain lift at a high angle of attack and provide higher maneuverability during the predation process [11].

**Figure 11.** Flipper's front structure of humpback whale [41]. Many scholars have studied and analyzed the flow control mechanism of wavy leading edge structure on the wing or other blades. The numerical results of Abdelrahman et al. [15] show that the wavy leading edge changes the flow field, and the flow separation on the upper surface of the wing is suppressed at the wavy leading edge, which has the effect of delaying stall. Torró et al. [18] found that the unsteady aerodynamic force of the airfoil at the periodic vortex shedding frequency is reduced by the wavy leading edge through the large eddy simulation method. In the study of hydrofoil, Rostamzadeh et al. [30] pointed out that the effect of wavy leading edge structure on the flow field is similar to the principle of vortex generator. By generating shedding vortex, the pressure difference between the upper and lower surfaces of the wing is increased to achieve the effect of increasing lift and reducing drag. The longitudinal vortex also causes the downwash velocity at the peak position of the leading edge of the wing, which changes the flow state of the fluid. Under the condition of low angle of attack, flow separation will also occur at the leading edge of the wavy shape, and the vortex structure presents the "bi-periodic" state [13,31]. The bi-periodic state means that when the fluid flows through the wavy leading edge, a pair of backflows will be formed at the trough, as shown in Figure 12. Therefore, under the influence of unsteady effect and vortex generated by wavy leading edge, even if the

Many scholars have studied and analyzed the flow control mechanism of wavy leading edge structure on the wing or other blades. The numerical results of Abdelrahman et

on the upper surface of the wing is suppressed at the wavy leading edge, which has the effect of delaying stall. Torró et al. [18] found that the unsteady aerodynamic force of the airfoil at the periodic vortex shedding frequency is reduced by the wavy leading edge through the large eddy simulation method. In the study of hydrofoil, Rostamzadeh et al. [30] pointed out that the effect of wavy leading edge structure on the flow field is similar to the principle of vortex generator. By generating shedding vortex, the pressure difference between the upper and lower surfaces of the wing is increased to achieve the effect of increasing lift and reducing drag. The longitudinal vortex also causes the downwash velocity at the peak position of the leading edge of the wing, which changes the flow state of the fluid. Under the condition of low angle of attack, flow separation will also occur at

**Figure 10.** Velocity vector diagram of 1% chordwise section at t/T = 0.25.

**3. Design of Wavy Leading Edge Flapping Wing** *3.1. Flow Control Mechanism of Wavy Leading Edge*

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 9 of 24

anisms of lift generation.

flapping wing does not exceed the stall angle of attack, its flow field characteristics will change. Figure 11 shows that the leading edge of the flippers of the humpback whale has special nodular structures. These wavy bulges can maintain lift at a high angle of attack and provide higher maneuverability during the predation process [11].

The asymmetry of the up and down flapping process and the strong leading edge vortex attached to the upper surface during the down flapping process are the important mech-

**Figure 11.** Flipper's front structure of humpback whale [41]. **Figure 11.** Flipper's front structure of humpback whale [41]. istics will change. fore, under the influence of unsteady effect and vortex generated by wavy leading edge,

of the fluid. Under the condition of low angle of attack, flow separation will also occur at **Figure 12.** Upper surface limit streamline of wavy leading edge wing. **Figure 12.** Upper surface limit streamline of wavy leading edge wing.

### *3.2. Characteristic Parameters of Wavy Leading Edge 3.2. Characteristic Parameters of Wavy Leading Edge*

Based on the above analysis, the aerodynamic performance of the wavy leading edge flapping wing with NACA0014 as the base section in the unsteady flapping state is studied. The structure of the wavy leading edge wing is shown in Figure 13. The sinusoidal wavy is defined by the two parameters of amplitude *A* and wavelength *W* and varies along the wing wingspan according to Equation (5). In addition, *x* is the spanwise length from any section to the wing root. The average chord length and projection area of the wavy leading edge flapping wing are the same as those of the reference flapping wing. Based on the above analysis, the aerodynamic performance of the wavy leading edge flapping wing with NACA0014 as the base section in the unsteady flapping state is studied. The structure of the wavy leading edge wing is shown in Figure 13. The sinusoidal wavy is defined by the two parameters of amplitude *A* and wavelength *W* and varies along the wing wingspan according to Equation (5). In addition, *x* is the spanwise length from any section to the wing root. The average chord length and projection area of the wavy leading edge flapping wing are the same as those of the reference flapping wing. **Figure 12.** Upper surface limit streamline of wavy leading edge wing. *3.2. Characteristic Parameters of Wavy Leading Edge* Based on the above analysis, the aerodynamic performance of the wavy leading edge flapping wing with NACA0014 as the base section in the unsteady flapping state is studied. The structure of the wavy leading edge wing is shown in Figure 13. The sinusoidal wavy is defined by the two parameters of amplitude *A* and wavelength *W* and varies along the wing wingspan according to Equation (5). In addition, *x* is the spanwise length

$$Y(\mathbf{x}) = \frac{A}{2}\sin(\frac{2\pi}{W}\mathbf{x})\tag{5}$$

(5)

(5)

(6)

(6)

**Figure 13.** Wavy leading edge wing.**Figure 13.** Wavy leading edge wing.

max

*baseline*

max

*baseline*

*baseline*

*x*

*avy baseline*

*x x*

*avy baseline*

*baseline*

*x*

*wavy*

<sup>=</sup>

 <sup>=</sup>

*x x*

<sup>=</sup>

 <sup>=</sup>

*w*

*w*

*z z*

*z z*

*wavy*

*x x x*

*x x x*

*wavy*

*wavy*

max

+ − − −

max

+ − − −

airfoil at the middle is the same as the airfoil of the baseline wing.

[ ( ( ) )] [ ( ) ]

airfoil at the middle is the same as the airfoil of the baseline wing.

[ ( ( ) )] [ ( ) ]

( )

> ( )

*Y x*

*Y x* 2

*W*

2

*W*

 )

> )

 *x*

 *x*

*wavy*

*wavy*

sin(

sin(

A sinusoidal wave is generated by five sections along the spanwise direction, as shown in Figure 14a. The airfoil data for each section are defined by Equation (6), and the

2

=

2

*A*

*A*

=

max max

*x Y x c Y x c x x*

max max

*x Y x c Y x c x x*

The

**Configuration**

*A*

*W*

*baseline <sup>x</sup>*

and

A sinusoidal wave is generated by five sections along the spanwise direction, as shown in Figure 14a. The airfoil data for each section are defined by Equation (6), and the airfoil at the middle is the same as the airfoil of the baseline wing. leading edge wing, respectively, and max *<sup>x</sup>* represents the abscissa at the maximum thickness of the baseline wing. The *baseline z* and *wavy z* represent the ordinates of the baseline

represent the abscissas of the baseline wing and the wavy

$$\begin{cases} \mathbf{x}\_{\text{array}} = \begin{cases} \frac{\mathbf{x}\_{\text{label}}}{\mathbf{x}\_{\text{label}}} \left[ \mathbf{x}\_{\text{max}} + \left( Y(\mathbf{x}) - c \right) \right] - \left[ Y(\mathbf{x}) - c \right] \mathbf{x}\_{\text{row}} \times \mathbf{x}\_{\text{max}} \\\\ \mathbf{z}\_{\text{array}} = \begin{cases} \mathbf{x}\_{\text{basicline}} & \mathbf{x}\_{\text{new}} \ge \mathbf{x}\_{\text{max}} \end{cases} \\\\ \end{cases} \tag{6}$$

$$\begin{array}{c} \mathbf{x}\_{\text{label}} \\\\ \hline \\\\ \mathbf{x}\_{\text{new}} \end{array} \tag{7}$$

wing and the wavy leading edge wing, respectively. When generating spanwise section

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 11 of 24

*wavy*

*x*

**Figure 14.** Construction method of wavy leading edge wing. (**a**) Wing spanwise section of wavy leading edge wing; (**b**) Airfoil of wavy leading edge wing. **Figure 14.** Construction method of wavy leading edge wing. (**a**) Wing spanwise section of wavy leading edge wing; (**b**) Airfoil of wavy leading edge wing.

For purpose of further clarifying the influence of design parameters on the aerodynamic performance and flow control effect of the bionic wavy leading edge flapping wing, dimensionless processing was carried out for amplitude and wavelength, as shown in Equation (7). Twelve different flapping wing configurations composed of three amplitudes and four wavelengths are selected for calculation and comparison. The specific parameters of the configurations are shown in Table 4, and the values of these parameters are within the range of protuberances found on the humpback whale's pectoral fin in nature [42]. Figure 15 shows the schematic diagrams of all twelve configurations. The com-The *xbaseline* and *xwavy* represent the abscissas of the baseline wing and the wavy leading edge wing, respectively, and *x*max represents the abscissa at the maximum thickness of the baseline wing. The *zbaseline* and *zwavy* represent the ordinates of the baseline wing and the wavy leading edge wing, respectively. When generating spanwise section airfoil data, only the geometry from the leading edge to the maximum thickness is changed, while the geometry from the maximum thickness to the trailing edge is kept unchanged. The baseline airfoil is stretched or compressed according to the chord length at different spanwise positions, as shown in Figure 14b.

parison of different design parameters of the wavy leading edge wing configurations can be seen through Figure 15. *A A c W W L* = <sup>=</sup> (7) **Table 4.** Wavy leading edge flapping wing configuration. For purpose of further clarifying the influence of design parameters on the aerodynamic performance and flow control effect of the bionic wavy leading edge flapping wing, dimensionless processing was carried out for amplitude and wavelength, as shown in Equation (7). Twelve different flapping wing configurations composed of three amplitudes and four wavelengths are selected for calculation and comparison. The specific parameters of the configurations are shown in Table 4, and the values of these parameters are within the range of protuberances found on the humpback whale's pectoral fin in nature [42]. Figure 15 shows the schematic diagrams of all twelve configurations. The comparison of different design parameters of the wavy leading edge wing configurations can be seen through Figure 15.

$$\begin{cases} \frac{\overline{A}}{W} = \frac{\overline{A}}{\frac{\overline{W}}{L}}\\ \end{cases} \tag{7}$$

**Configuration**

1-1 0.05 0.1 2-1 0.1 0.1 3-1 0.2 0.1 1-2 0.05 0.05 2-2 0.1 0.05 3-2 0.2 0.05 1-3 0.05 0.02 2-3 0.1 0.02 3-3 0.2 0.02 1-4 0.05 0.01 2-4 0.1 0.01 3-4 0.2 0.01

*W* = 0.01

Trust coefficient.


**Table 4.** Wavy leading edge flapping wing configuration.

**Figure 15.** Wavy leading edge flapping wing configuration. (**a**) Configuration 1-1; (**b**) Configuration 1-2; (**c**) Configuration 1-3; (**d**) Configuration 1-4; (**e**) Configuration 2-1; (**f**) Configuration 2-4; (**g**) Configuration 1-3; (**h**) Configuration 2-4; (**i**) Configuration 3-1; (**j**) Configuration 3-2; (**k**) Configuration 3-3; (**l**) Configuration 3-4. **Figure 15.** Wavy leading edge flapping wing configuration. (**a**) Configuration 1-1; (**b**) Configuration 1-2; (**c**) Configuration 1-3; (**d**) Configuration 1-4; (**e**) Configuration 2-1; (**f**) Configuration 2-4; (**g**) Configuration 1-3; (**h**) Configuration 2-4; (**i**) Configuration 3-1; (**j**) Configuration 3-2; (**k**) Configuration 3-3; (**l**) Configuration 3-4.

### *3.3. Grid Independence Verification*

can meet the grid independence requirements.

*3.3. Grid Independence Verification* Due to the difference between the wavy leading edge flapping wing and the baseline wing in the model, Configuration 1-1 ( *A* = 0.05,*W* = 0.1 ) and Configuration 3-4 ( *A* = 0.2 , ) of the wavy leading edge are selected for grid independence verification to ensure the accuracy of the numerical simulation. The same background grid as the straight wing is used, and the number of grids is 1.35 million. The partitioning strategy of the component grid is also the same as the three cases in Section 2.3. Figures 16 and 17, respectively, show the calculation results of the lift coefficient and thrust coefficient of Con-Due to the difference between the wavy leading edge flapping wing and the baseline wing in the model, Configuration 1-1 (*A* = 0.05, *W* = 0.1) and Configuration 3-4 (*A* = 0.2, *W* = 0.01) of the wavy leading edge are selected for grid independence verification to ensure the accuracy of the numerical simulation. The same background grid as the straight wing is used, and the number of grids is 1.35 million. The partitioning strategy of the component grid is also the same as the three cases in Section 2.3. Figures 16 and 17, respectively, show the calculation results of the lift coefficient and thrust coefficient of Configuration 1-1 and Configuration 3-4 in one period with time. It can be seen that the Case2 can meet the grid independence requirements.

(**a**) (**b**)

**Figure 16.** The time history of lift and thrust coefficient (Configuration 11). (**a**) Lift coefficient; (**b**)

figuration 1-1 and Configuration 3-4 in one period with time. It can be seen that the Case2

(**a**) (**b**)

(**c**) (**d**)

(**e**) (**f**)

(**g**) (**h**)

(**i**) (**j**)

(**k**) (**l**) **Figure 15.** Wavy leading edge flapping wing configuration. (**a**) Configuration 1-1; (**b**) Configuration 1-2; (**c**) Configuration 1-3; (**d**) Configuration 1-4; (**e**) Configuration 2-1; (**f**) Configuration 2-4; (**g**) Configuration 1-3; (**h**) Configuration 2-4; (**i**) Configuration 3-1; (**j**) Configuration 3-2; (**k**) Configura-

Due to the difference between the wavy leading edge flapping wing and the baseline

sure the accuracy of the numerical simulation. The same background grid as the straight wing is used, and the number of grids is 1.35 million. The partitioning strategy of the component grid is also the same as the three cases in Section 2.3. Figures 16 and 17, respectively, show the calculation results of the lift coefficient and thrust coefficient of Configuration 1-1 and Configuration 3-4 in one period with time. It can be seen that the Case2

) of the wavy leading edge are selected for grid independence verification to en-

) and Configuration 3-4 (

*A* = 0.2,

*A* = 0.05,*W* = 0.1

**Figure 16.** The time history of lift and thrust coefficient (Configuration 11). (**a**) Lift coefficient; (**b**) Trust coefficient. **Figure 16.** The time history of lift and thrust coefficient (Configuration 11). (**a**) Lift coefficient; (**b**) Trust coefficient. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 13 of 24

**Figure 17.** The time history of lift and thrust coefficient (Configuration 3-4). (**a**) Lift coefficient; (**b**) Trust coefficient. **Figure 17.** The time history of lift and thrust coefficient (Configuration 3-4). (**a**) Lift coefficient; (**b**) Trust coefficient.

### **4. Numerical Results 4. Numerical Results**

aged thrust coefficient

tion (8).

decreases.

*CT*

*C*

=

tion 3-3; (**l**) Configuration 3-4.

*W* = 0.01

*3.3. Grid Independence Verification*

wing in the model, Configuration 1-1 (

can meet the grid independence requirements.

In order to compare the influence of wavy leading edge on the aerodynamic performance of flapping wing more directly, the calculation results are divided into two parts according to different wavelengths and different amplitudes. The aerodynamic performance of flapping wing is measured by time averaged lift coefficient *CL* and time aver-In order to compare the influence of wavy leading edge on the aerodynamic performance of flapping wing more directly, the calculation results are divided into two parts according to different wavelengths and different amplitudes. The aerodynamic performance of flapping wing is measured by time averaged lift coefficient *C<sup>L</sup>* and time averaged thrust coefficient *C<sup>T</sup>* in a flapping cycle. The two parameters are defined as Equation (8).

$$\begin{cases} \begin{array}{c} \overline{\mathcal{C}\_{L}} = \frac{1}{T} \int\_{0}^{T} \mathcal{C}\_{L}(t)dt\\ \overline{\mathcal{C}\_{T}} = \frac{1}{T} \int\_{0}^{T} -\mathcal{C}\_{D}(t)dt \end{array} = \frac{1}{T} \int\_{0}^{T} \mathcal{C}\_{T}(t)dt \end{cases} \tag{8}$$

in a flapping cycle. The two parameters are defined as Equa-

#### 0 ()*L L CT 4.1. Effects of Different Wavelengths on Aerodynamic Performance of Flapping Wing*

 *t dt*

0 0 1 1 ( ) ( ) *T T T D T C C t dt C t dt T T* = − = (8) *4.1. Effects of Different Wavelengths on Aerodynamic Performance of Flapping Wing* Figures 18 and 19, respectively, show the time history comparison of lift coefficient and thrust coefficient of straight flapping wing and wavy leading edge flapping wing with different wavelengths in one flapping cycle. It can be seen that the lift coefficient of the flapping wing with the wavy leading edge is higher than that of the straight wing in the Figures 18 and 19, respectively, show the time history comparison of lift coefficient and thrust coefficient of straight flapping wing and wavy leading edge flapping wing with different wavelengths in one flapping cycle. It can be seen that the lift coefficient of the flapping wing with the wavy leading edge is higher than that of the straight wing in the process of downward flapping (t/T = 0–0.5) and in the last quarter stage of upward flapping (t/T = 0.875–1). In the first three quarters of upward flapping (t/T = 0.5–0.875), the lift coefficient is similar to or slightly lower than that of the straight wing. The thrust coefficient is close to that of the straight wing in the whole process of the downward flapping (t/T = 0–0.5) and the first half of the upward flapping (t/T = 0.5–0.75), and it

process of downward flapping (t/T = 0–0.5) and in the last quarter stage of upward flapping (t/T = 0.875–1). In the first three quarters of upward flapping (t/T = 0.5–0.875), the lift

(t/T = 0–0.5) and the first half of the upward flapping (t/T = 0.5–0.75), and it decreases obviously in the latter half (t/T = 0.75–1), which has strong regularity. That is, the smaller the wavelength of the wavy leading edge flapping wing, the more the thrust coefficient

decreases obviously in the latter half (t/T = 0.75–1), which has strong regularity. That is, the smaller the wavelength of the wavy leading edge flapping wing, the more the thrust coefficient decreases. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 14 of 24 *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 14 of 24

**Figure 18.** Time history of lift coefficient in a flapping cycle. (**a**) *A* = 0.05 ; (**b**) *A* = 0.1 ; (**c**) *A* = 0.2 . **Figure 18.** Time history of lift coefficient in a flapping cycle. (**a**) *A* = 0.05; (**b**) *A* = 0.1; (**c**) *A* = 0.2. **Figure 18.** Time history of lift coefficient in a flapping cycle. (**a**) *A*=0.05; (**b**) *A*=0.1; (**c**) *A*=0.2.

**Figure 19.** Time history of trust coefficient in a flapping cycle. (**a**)

. **Figure 19.** Time history of trust coefficient in a flapping cycle. (**a**) *A* = 0.05 ; (**b**) *A* = 0.1 ; (**c**) *A* = 0.2 . **Figure 19.** Time history of trust coefficient in a flapping cycle. (**a**) *A* = 0.05; (**b**) *A* = 0.1; (**c**) *A* = 0.2.

*A* = 0.05

; (**b**)

*A* = 0.1

; (**c**)

*A* = 0.2

Figure 20 compares the time averaged lift and thrust coefficients of the straight wing and wavy leading edge wings. The wavy leading edge inevitably reduces the thrust coefficient while increasing the lift coefficient, and the smaller the wavelength, the larger the time averaged lift coefficient and the smaller the time averaged thrust coefficient. Figures 21 and 22, respectively, show the time-averaged lift coefficient and thrust coefficient during downstroke (t/T = 0–0.5) and upstroke (t/T = 0.5–1). The wavy leading edge wings mainly increases the lift in the downstroke, while the thrust decreases in both the up-anddown stroke of the flapping. Table 5 shows the average growth rate of time averaged lift coefficient relative to straight wing and the average reduction rate of time averaged thrust coefficient for flapping wing with different wavelengths. When the wavelength is *W* = 0.01 , the average lift coefficient increases the most, which is 41.04% of the straight wing, but at the same time, the average thrust coefficient decreases by 25.97%. Therefore, it is necessary to weigh the lift and thrust of the flapping wing during design, or to increase the thrust by increasing the flapping frequency. Wavy leading edge flapping wing can withstand greater take-off weight, and it has better low-speed flight ability. Figure 20 compares the time averaged lift and thrust coefficients of the straight wing and wavy leading edge wings. The wavy leading edge inevitably reduces the thrust coefficient while increasing the lift coefficient, and the smaller the wavelength, the larger the time averaged lift coefficient and the smaller the time averaged thrust coefficient. Figures 21 and 22, respectively, show the time-averaged lift coefficient and thrust coefficient during downstroke (t/T = 0–0.5) and upstroke (t/T = 0.5–1). The wavy leading edge wings mainly increases the lift in the downstroke, while the thrust decreases in both the up-anddown stroke of the flapping. Table 5 shows the average growth rate of time averaged lift coefficient relative to straight wing and the average reduction rate of time averaged thrust coefficient for flapping wing with different wavelengths. When the wavelength is *W* = 0.01 , the average lift coefficient increases the most, which is 41.04% of the straight wing, but at the same time, the average thrust coefficient decreases by 25.97%. Therefore, it is necessary to weigh the lift and thrust of the flapping wing during design, or to increase the thrust by increasing the flapping frequency. Wavy leading edge flapping wing can withstand greater take-off weight, and it has better low-speed flight ability. Figure 20 compares the time averaged lift and thrust coefficients of the straight wing and wavy leading edge wings. The wavy leading edge inevitably reduces the thrust coefficient while increasing the lift coefficient, and the smaller the wavelength, the larger the time averaged lift coefficient and the smaller the time averaged thrust coefficient. Figures 21 and 22, respectively, show the time-averaged lift coefficient and thrust coefficient during downstroke (t/T = 0–0.5) and upstroke (t/T = 0.5–1). The wavy leading edge wings mainly increases the lift in the downstroke, while the thrust decreases in both the up-anddown stroke of the flapping. Table 5 shows the average growth rate of time averaged lift coefficient relative to straight wing and the average reduction rate of time averaged thrust coefficient for flapping wing with different wavelengths. When the wavelength is *W* = 0.01, the average lift coefficient increases the most, which is 41.04% of the straight wing, but at the same time, the average thrust coefficient decreases by 25.97%. Therefore, it is necessary to weigh the lift and thrust of the flapping wing during design, or to increase the thrust by increasing the flapping frequency. Wavy leading edge flapping wing can withstand greater take-off weight, and it has better low-speed flight ability.

.

*W* = 0.01

(**a**) (**b**) (**c**) **Figure 18.** Time history of lift coefficient in a flapping cycle. (**a**)

(**a**) (**b**) (**c**) **Figure 19.** Time history of trust coefficient in a flapping cycle. (**a**)

**Figure 20.** Comparison of average lift coefficient and thrust coefficient. (**a**) Time averaged lift coefficient; (**b**) Time averaged trust coefficient. **Figure 20.** Comparison of average lift coefficient and thrust coefficient. (**a**) Time averaged lift coefficient; (**b**) Time averaged trust coefficient. **Figure 20.** Comparison of average lift coefficient and thrust coefficient. (**a**) Time averaged lift co-

*A* = 0.05

> *A* = 0.05

Figure 20 compares the time averaged lift and thrust coefficients of the straight wing and wavy leading edge wings. The wavy leading edge inevitably reduces the thrust coefficient while increasing the lift coefficient, and the smaller the wavelength, the larger the time averaged lift coefficient and the smaller the time averaged thrust coefficient. Figures 21 and 22, respectively, show the time-averaged lift coefficient and thrust coefficient during downstroke (t/T = 0–0.5) and upstroke (t/T = 0.5–1). The wavy leading edge wings mainly increases the lift in the downstroke, while the thrust decreases in both the up-anddown stroke of the flapping. Table 5 shows the average growth rate of time averaged lift coefficient relative to straight wing and the average reduction rate of time averaged thrust coefficient for flapping wing with different wavelengths. When the wavelength is

, the average lift coefficient increases the most, which is 41.04% of the straight

wing, but at the same time, the average thrust coefficient decreases by 25.97%. Therefore,

can withstand greater take-off weight, and it has better low-speed flight ability.

; (**b**)

; (**b**)

*A* = 0.1 ; (**c**)

*A* = 0.2

*A* = 0.1 ; (**c**)

*A* = 0.2.

**Figure 21.** Comparison of average lift coefficient of downstroke and upstroke. (**a**) Time averaged lift coefficient of downstroke; (**b**) Time averaged lift coefficient of upstroke. **Figure 21.** Comparison of average lift coefficient of downstroke and upstroke. (**a**) Time averaged lift coefficient of downstroke; (**b**) Time averaged lift coefficient of upstroke. **Figure 21.** Comparison of average lift coefficient of downstroke and upstroke. (**a**) Time averaged lift coefficient of downstroke; (**b**) Time averaged lift coefficient of upstroke.

**Figure 22.** Comparison of average thrust coefficient of downstroke and upstroke. (**a**) Time averaged

**Table 5.** Average change rates of flapping wings at different wavelengths.

**Table 5.** Average change rates of flapping wings at different wavelengths.

*CL*

*CL*

*W*

*W*

efficient; (**b**) Time averaged trust coefficient.

Time averaged trust coefficient

*CT*

Time averaged trust coefficient

*CT*

*W* = 0

*W* = 0

Time averaged lift coefficient

Time averaged lift coefficient

the straight wing (

the straight wing (

trust coefficient of downstroke; (**b**) Time averaged trust coefficient of upstroke. **Figure 22.** Comparison of average thrust coefficient of downstroke and upstroke. (**a**) Time averaged trust coefficient of downstroke; (**b**) Time averaged trust coefficient of upstroke. **Figure 22.** Comparison of average thrust coefficient of downstroke and upstroke. (**a**) Time averaged trust coefficient of downstroke; (**b**) Time averaged trust coefficient of upstroke.

**0.1 0.05 0.02 0.01**

16.01%↑ 28.89%↑ 34.32%↑ 41.04%↑

**0.1 0.05 0.02 0.01**

16.01%↑ 28.89%↑ 34.32%↑ 41.04%↑

11.66%↓ 14.04%↓ 20.71%↓ 25.97%↓

11.66%↓ 14.04%↓ 20.71%↓ 25.97%↓

*A* = 0.1

*A* = 0.1 as an

as an

moments (t/T = 0, 0.2, 0.5 and 0.8) in a flapping cycle are compared and analyzed by taking

With the purpose of further clarifying the mechanism of wavy leading edge changing the characteristics of flow field and the influence of different wavelengths, four typical moments (t/T = 0, 0.2, 0.5 and 0.8) in a flapping cycle are compared and analyzed by taking

example. Figure 23 shows the pressure contours of the upper surface of the flapping wing, and the unit of pressure is Pa. In the downward stage (t/T = 0 and 0.2), the pressure on the upper surface of the wavy leading edge flapping wing decreases, while in the upward stage (t/T = 0.5 and 0.8), the pressure on the upper surface increases. The change of pressure is mainly reflected in the position above 60% of the wingspan, and the change of the pressure difference between the upper and lower surfaces will cause the change of lift.

example. Figure 23 shows the pressure contours of the upper surface of the flapping wing, and the unit of pressure is Pa. In the downward stage (t/T = 0 and 0.2), the pressure on the upper surface of the wavy leading edge flapping wing decreases, while in the upward stage (t/T = 0.5 and 0.8), the pressure on the upper surface increases. The change of pressure is mainly reflected in the position above 60% of the wingspan, and the change of the pressure difference between the upper and lower surfaces will cause the change of lift.

) and the wavy leading edge flapping wing with

) and the wavy leading edge flapping wing with


**Table 5.** Average change rates of flapping wings at different wavelengths.

With the purpose of further clarifying the mechanism of wavy leading edge changing the characteristics of flow field and the influence of different wavelengths, four typical moments (t/T = 0, 0.2, 0.5 and 0.8) in a flapping cycle are compared and analyzed by taking the straight wing (*W* = 0) and the wavy leading edge flapping wing with *A* = 0.1 as an example. Figure 23 shows the pressure contours of the upper surface of the flapping wing, and the unit of pressure is Pa. In the downward stage (t/T = 0 and 0.2), the pressure on the upper surface of the wavy leading edge flapping wing decreases, while in the upward stage (t/T = 0.5 and 0.8), the pressure on the upper surface increases. The change of pressure is mainly reflected in the position above 60% of the wingspan, and the change of the pressure difference between the upper and lower surfaces will cause the change of lift. Figure 24 compares the pressure coefficients of the straight wing and the wavy leading edge wing with different wavelengths at 80% wingspan section. When t/T = 0 and 0.2, the pressure of the upper surface of the wavy leading edge flapping wing decreases compared with the straight wing, and the pressure of the lower surface increases. The pressure difference between the upper and lower surfaces increases, and the lift increases. Although Figure 24a is similar to Figure 24c, and Figure 24b is similar to Figure 24d, it can be seen from the pressure contours that the pressure on the upper surface of the flapping wing is greater than that on the lower surface. Therefore, after the wavy leading edge causes the pressure difference to increase, the downward force on the flapping wing is further increased. That is, the lift is reduced. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 16 of 24 Figure 24 compares the pressure coefficients of the straight wing and the wavy leading edge wing with different wavelengths at 80% wingspan section. When t/T = 0 and 0.2, the pressure of the upper surface of the wavy leading edge flapping wing decreases compared

Figure 25 shows the pressure coefficient contours and streamlines of the 50% spanwise section of the flapping wing when t/T = 0. As the wavelength decreases, the vortex generated by the wavy leading edge decreases the pressure coefficient on the upper surface of the wing section and increases the pressure coefficient on the lower surface, resulting in an increase in the lift coefficient. In addition, the pressure coefficient on the trailing edge decreases with the decrease of the wavelength, which leads to a decrease in the thrust coefficient of the flapping wing. with the straight wing, and the pressure of the lower surface increases. The pressure difference between the upper and lower surfaces increases, and the lift increases. Although Figure 24a is similar to Figure 24c, and Figure 24b is similar to Figure 24d, it can be seen from the pressure contours that the pressure on the upper surface of the flapping wing is greater than that on the lower surface. Therefore, after the wavy leading edge causes the pressure difference to increase, the downward force on the flapping wing is further increased. That is, the lift is reduced.

(**c**) (**d**)

(**c**) t/T = 0.5; (**d**) t/T = 0.8.

**Figure 23.** Pressure contours of the upper surface of flapping wing (Pa). (**a**) t/T = 0; (**b**) t/T = 0.2;

creased. That is, the lift is reduced.

**Figure 23.** Pressure contours of the upper surface of flapping wing (Pa). (**a**) t/T = 0; (**b**) t/T = 0.2; (**c**) t/T = 0.5; (**d**) t/T = 0.8. **Figure 23.** Pressure contours of the upper surface of flapping wing (Pa). (**a**) t/T = 0; (**b**) t/T = 0.2; (**c**) t/T = 0.5; (**d**) t/T = 0.8. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 17 of 24

Figure 24 compares the pressure coefficients of the straight wing and the wavy leading edge wing with different wavelengths at 80% wingspan section. When t/T = 0 and 0.2, the pressure of the upper surface of the wavy leading edge flapping wing decreases compared with the straight wing, and the pressure of the lower surface increases. The pressure difference between the upper and lower surfaces increases, and the lift increases. Although Figure 24a is similar to Figure 24c, and Figure 24b is similar to Figure 24d, it can be seen from the pressure contours that the pressure on the upper surface of the flapping wing is greater than that on the lower surface. Therefore, after the wavy leading edge causes the pressure difference to increase, the downward force on the flapping wing is further in-

**Figure 24.** Pressure coefficient of 80% wingspan section. (**a**) t/T = 1; (**b**) t/T = 0.2; (**c**) t/T = 0.5; (**d**) t/T = 0.8. **Figure 24.** Pressure coefficient of 80% wingspan section. (**a**) t/T = 1; (**b**) t/T = 0.2; (**c**) t/T = 0.5; (**d**) t/T = 0.8.

### Figure 25 shows the pressure coefficient contours and streamlines of the 50% *4.2. Effects of Different Amplitudes on Aerodynamic Performance of Flapping Wing*

**Figure 25.** Pressure contours and streamlines of 50% spanwise section of flapping wing at t/T = 0.

*W* = 0.02 ; (**e**)

*W* = 0.01.

; (**d**)

*4.2. Effects of Different Amplitudes on Aerodynamic Performance of Flapping Wing*

spanwise section of the flapping wing when t/T = 0. As the wavelength decreases, the vortex generated by the wavy leading edge decreases the pressure coefficient on the upper surface of the wing section and increases the pressure coefficient on the lower surface, resulting in an increase in the lift coefficient. In addition, the pressure coefficient on the trailing edge decreases with the decrease of the wavelength, which leads to a decrease in Figures 26 and 27, respectively, show the time history comparison of lift coefficient and thrust coefficient of straight flapping wing and wavy leading edge flapping wing with different amplitudes in one flapping cycle. When the wavelength is large, the change of amplitude has little effect on the lift coefficient, and at the same wavelength, the larger the amplitude, the more the thrust coefficient decreases.

(**a**) (**b**) (**c**)

(**d**) (**e**)

*W* = 0.05

; (**c**)

*W* = 0.1

(**a**)Straight Wing; (**b**)

the thrust coefficient of the flapping wing.

t/T = 0.8.

the thrust coefficient of the flapping wing.

amplitude, the more the thrust coefficient decreases.

**Figure 25.** Pressure contours and streamlines of 50% spanwise section of flapping wing at t/T = 0. (**a**)Straight Wing; (**b**) *W* = 0.1 ; (**c**) *W* = 0.05 ; (**d**) *W* = 0.02 ; (**e**) *W* = 0.01. **Figure 25.** Pressure contours and streamlines of 50% spanwise section of flapping wing at t/T = 0. (**a**) Straight Wing; (**b**) *W* = 0.1; (**c**) *W* = 0.05; (**d**) *W* = 0.02; (**e**) *W* = 0.01. different amplitudes in one flapping cycle. When the wavelength is large, the change of amplitude has little effect on the lift coefficient, and at the same wavelength, the larger the

*4.2. Effects of Different Amplitudes on Aerodynamic Performance of Flapping Wing*

(**a**) (**b**)

(**c**) (**d**) **Figure 24.** Pressure coefficient of 80% wingspan section. (**a**) t/T = 1; (**b**) t/T = 0.2; (**c**) t/T = 0.5; (**d**)

Figure 25 shows the pressure coefficient contours and streamlines of the 50% spanwise section of the flapping wing when t/T = 0. As the wavelength decreases, the vortex generated by the wavy leading edge decreases the pressure coefficient on the upper surface of the wing section and increases the pressure coefficient on the lower surface, resulting in an increase in the lift coefficient. In addition, the pressure coefficient on the trailing edge decreases with the decrease of the wavelength, which leads to a decrease in

**Figure 26.** Time history of lift coefficient in a flapping cycle. (**a**) *W* 0.1 ; (**b**) *W* 0.05 ; (**c**) *W* 0.02 ; (**d**) *W* 0.01. **Figure 26.** Time history of lift coefficient in a flapping cycle. (**a**) *W* = 0.1; (**b**) *W* = 0.05; (**c**) *W* = 0.02; (**d**) *W* = 0.01.

Figure 28 compares the time averaged lift and thrust coefficients of the straight wing and wavy leading edge wings with different amplitudes, and Figures 29 and 30, respectively, show the time-averaged lift coefficient and thrust coefficient during downstroke (t/T = 0–0.5) and upstroke (t/T = 0.5–1). Table 6 shows the average growth rate of time averaged lift coefficient and the average reduction rate of time averaged thrust coefficient of them. As can be seen, at the same wavelength, the larger the amplitude, the higher the average lift coefficient and the more the time-averaged thrust coefficient decreases.

(**c**) (**d**)

 = 0.01.

; (**d**)

Figures 26 and 27, respectively, show the time history comparison of lift coefficient

and thrust coefficient of straight flapping wing and wavy leading edge flapping wing with different amplitudes in one flapping cycle. When the wavelength is large, the change of amplitude has little effect on the lift coefficient, and at the same wavelength, the larger the

(**a**) (**b**)

(**c**) (**d**)

*W* = 0.1 ; (**b**)

*W* = 0.05 ; (**c**)

*W* = 0.02

; (**c**)

amplitude, the more the thrust coefficient decreases.

**Figure 27.** Time history of trust coefficient in a flapping cycle. (**a**) *W* = 0.1; (**b**) *W* = 0.05; (**c**) *W* = 0.02; (**d**) *W* = 0.01. of them. As can be seen, at the same wavelength, the larger the amplitude, the higher the average lift coefficient and the more the time-averaged thrust coefficient decreases.

**Figure 28.** Comparison of average lift coefficient and thrust coefficient. (**a**) Time averaged lift coefficient; (**b**) Time averaged trust coefficient. **Figure 28.** Comparison of average lift coefficient and thrust coefficient. (**a**) Time averaged lift coefficient; (**b**) Time averaged trust coefficient.

(**a**) (**b**) **Figure 29.** Comparison of average lift coefficient of downstroke and upstroke. (**a**)Time averaged

lift coefficient of downstroke; (**b**) Time averaged lift coefficient of upstroke.

(**a**) (**b**)

**Figure 28.** Comparison of average lift coefficient and thrust coefficient. (**a**) Time averaged lift

**Figure 27.** Time history of trust coefficient in a flapping cycle. (**a**)

Figure 28 compares the time averaged lift and thrust coefficients of the straight wing and wavy leading edge wings with different amplitudes, and Figures 29 and 30, respectively, show the time-averaged lift coefficient and thrust coefficient during downstroke (t/T = 0–0.5) and upstroke (t/T = 0.5–1). Table 6 shows the average growth rate of time averaged lift coefficient and the average reduction rate of time averaged thrust coefficient of them. As can be seen, at the same wavelength, the larger the amplitude, the higher the average lift coefficient and the more the time-averaged thrust coefficient decreases.

*W* = 0.02 ; (**d**)

*W* = 0.01.

coefficient; (**b**) Time averaged trust coefficient.

*W* = 0.1 ; (**b**)

*W* = 0.05 ; (**c**)

**Figure 29.** Comparison of average lift coefficient of downstroke and upstroke. (**a**)Time averaged lift coefficient of downstroke; (**b**) Time averaged lift coefficient of upstroke. **Figure 29.** Comparison of average lift coefficient of downstroke and upstroke. (**a**)Time averaged lift coefficient of downstroke; (**b**) Time averaged lift coefficient of upstroke. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 20 of 24

**Figure 30.** Comparison of average thrust coefficient of downstroke and upstroke. (**a**) Time averaged trust coefficient of downstroke; (**b**) Time averaged trust coefficient of upstroke. **Figure 30.** Comparison of average thrust coefficient of downstroke and upstroke. (**a**) Time averaged trust coefficient of downstroke; (**b**) Time averaged trust coefficient of upstroke.


**Table 6.** Average change rates of flapping wings at different amplitudes. **Table 6.** Average change rates of flapping wings at different amplitudes.

Taking the straight wing ( *A* = 0 ) and the wavy leading edge flapping wing with different amplitudes at wavelength *W* = 0.02 as an example, the four typical moments in a flapping cycle (t/T = 0, 0.2, 0.5, and 0.8) are compared and analyzed. The pressure contours of the upper surface of the flapping wing are shown in Figure 31; the unit of pressure is Pa. It can be seen that the wavy leading edge changes the pressure distribution on the flapping wing. In the downward flapping stage, the larger the amplitude, the smaller the pressure behind the trough, the smaller the pressure on the upper surface of the flapping Taking the straight wing (*A* = 0) and the wavy leading edge flapping wing with different amplitudes at wavelength *W* = 0.02 as an example, the four typical moments in a flapping cycle (t/T = 0, 0.2, 0.5, and 0.8) are compared and analyzed. The pressure contours of the upper surface of the flapping wing are shown in Figure 31; the unit of pressure is Pa. It can be seen that the wavy leading edge changes the pressure distribution on the flapping wing. In the downward flapping stage, the larger the amplitude, the smaller the pressure behind the trough, the smaller the pressure on the upper surface of the flapping wing, the larger the pressure difference between the upper and lower surfaces, and the higher the lift.

wing, the larger the pressure difference between the upper and lower surfaces, and the

(**a**) (**b**)

higher the lift.

**Figure 31.** Pressure contours of the upper surface of flapping wing. (**a**) t/T=0; (**b**) t/T = 0.2; (**c**) t/T = 0.5; (**d**) t/T = 0.8. **Figure 31.** Pressure contours of the upper surface of flapping wing. (**a**) t/T=0; (**b**) t/T = 0.2; (**c**) t/T = 0.5; (**d**) t/T = 0.8.

(**a**) (**b**) **Figure 30.** Comparison of average thrust coefficient of downstroke and upstroke. (**a**) Time averaged

flapping cycle (t/T = 0, 0.2, 0.5, and 0.8) are compared and analyzed. The pressure contours of the upper surface of the flapping wing are shown in Figure 31; the unit of pressure is Pa. It can be seen that the wavy leading edge changes the pressure distribution on the flapping wing. In the downward flapping stage, the larger the amplitude, the smaller the pressure behind the trough, the smaller the pressure on the upper surface of the flapping wing, the larger the pressure difference between the upper and lower surfaces, and the

**0.05 0.1 0.2**

24.34%↑ 29.33%↑ 36.52%↑

8.75%↓ 15.53%↓ 29.03%↓

) and the wavy leading edge flapping wing with dif-

as an example, the four typical moments in a

trust coefficient of downstroke; (**b**) Time averaged trust coefficient of upstroke.

**Table 6.** Average change rates of flapping wings at different amplitudes.

*CL*

*A* = 0

> *W* = 0.02

*A*

Time averaged trust coefficient

*CT*

Taking the straight wing (

ferent amplitudes at wavelength

higher the lift.

Time averaged lift coefficient

Figure 32 shows the Mach number contours at 50% spanwise section of the straight wing and the wavy leading edge flapping wing with different amplitudes when t/T = 0, where is the trough position. When the gas flows through the wavy leading edge, the airflow on the upper surface of the wing accelerates, and the airflow on the lower surface decelerates. Moreover, the larger the amplitude, the more obvious the change in velocity, so the greater the lift increase. Figure 32 shows the Mach number contours at 50% spanwise section of the straight wing and the wavy leading edge flapping wing with different amplitudes when t/T = 0, where is the trough position. When the gas flows through the wavy leading edge, the airflow on the upper surface of the wing accelerates, and the airflow on the lower surface decelerates. Moreover, the larger the amplitude, the more obvious the change in velocity, so the greater the lift increase.

**Figure 32.** Mach number contours of 50% spanwise section of flapping wing at t/T = 0. (**a**) Straight

In order to further analyze the influence of design parameters of wavy leading edge flapping wing on its aerodynamic performance, the sensitivities of time averaged lift coefficient and thrust coefficient to wavelength and amplitude are calculated, respectively, as shown in Table 7. At the same amplitude, the larger the wavelength of the wavy leading edge flapping wing, the smaller the time averaged lift coefficient and the larger the time averaged thrust coefficient. At the same wavelength, the larger the amplitude, the larger the time averaged lift coefficient and the smaller the time averaged thrust coefficient.

**Wavelength**

*W*

−77.4% 22.6%

**Amplitude**

*A*

; (**d**)

*A* = 0.2.

*CL*

Wing; (**b**)

*A* = 0.05 ; (**c**)

**Table 7.** Sensitivity of design parameters.

Time averaged lift coefficient

*4.3. Sensitivity Analysis of Design Parameters*

*A* = 0.1

(**c**) (**d**)

= 0.5; (**d**) t/T = 0.8.

so the greater the lift increase.

**Figure 32.** Mach number contours of 50% spanwise section of flapping wing at t/T = 0. (**a**) Straight Wing; (**b**) *A* = 0.05 ; (**c**) *A* = 0.1 ; (**d**) *A* = 0.2 . **Figure 32.** Mach number contours of 50% spanwise section of flapping wing at t/T = 0. (**a**) Straight Wing; (**b**) *A* = 0.05; (**c**) *A* = 0.1; (**d**) *A* = 0.2.

**Figure 31.** Pressure contours of the upper surface of flapping wing. (**a**) t/T=0; (**b**) t/T = 0.2; (**c**) t/T

Figure 32 shows the Mach number contours at 50% spanwise section of the straight wing and the wavy leading edge flapping wing with different amplitudes when t/T = 0, where is the trough position. When the gas flows through the wavy leading edge, the airflow on the upper surface of the wing accelerates, and the airflow on the lower surface decelerates. Moreover, the larger the amplitude, the more obvious the change in velocity,

#### *4.3. Sensitivity Analysis of Design Parameters 4.3. Sensitivity Analysis of Design Parameters*

In order to further analyze the influence of design parameters of wavy leading edge flapping wing on its aerodynamic performance, the sensitivities of time averaged lift coefficient and thrust coefficient to wavelength and amplitude are calculated, respectively, as shown in Table 7. At the same amplitude, the larger the wavelength of the wavy leading edge flapping wing, the smaller the time averaged lift coefficient and the larger the time averaged thrust coefficient. At the same wavelength, the larger the amplitude, the larger the time averaged lift coefficient and the smaller the time averaged thrust coefficient. In order to further analyze the influence of design parameters of wavy leading edge flapping wing on its aerodynamic performance, the sensitivities of time averaged lift coefficient and thrust coefficient to wavelength and amplitude are calculated, respectively,as shown in Table 7. At the same amplitude, the larger the wavelength of the wavy leading edge flapping wing, the smaller the time averaged lift coefficient and the larger the time averaged thrust coefficient. At the same wavelength, the larger the amplitude, the larger the time averaged lift coefficient and the smaller the time averaged thrust coefficient.

**Table 7.** Sensitivity of design parameters. **Table 7.** Sensitivity of design parameters.


To obtain the maximum lift coefficient while losing the least thrust, the minimum wavelength and the smallest amplitude should be selected. That is, in the configuration calculated in this paper, the configuration with wavelength *W* = 0.01 and amplitude *A* = 0.05 is optimal. Compared with the straight wing, the time averaged lift coefficient is increased by 32.86%, and the time averaged thrust coefficient is reduced by 14.28%.

### **5. Conclusions**

In this paper, a flow control method is proposed to apply the wavy leading edge structure to the straight symmetrical flapping wing to improve the average flapping lift. The aerodynamic performance and flow field characteristics of straight flapping wing and wavy leading edge flapping wing with different design parameters are calculated and compared by numerical simulation method. The mechanism of improving average lift is analyzed, and the sensitivity of the design parameters of wavy leading edge is carried out. The calculation results show that:


the configuration calculated in this paper, the configuration 1–4 with wavelength *W* = 0.01 and amplitude *A* = 0.05 is the best. Compared with the straight wing, this configuration can increase the average lift coefficient by 32.86% and reduce the average thrust coefficient by 14.28%.

**Author Contributions:** Conceptualization, B.M.; methodology, X.B.; data curation, X.B.; writing original draft, X.B.; writing—review & editing, B.M.; supervision, H.Z.; project administration, H.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** National Natural Science Foundation of China (Grant No. 12202363); Natural Science Basic Research Program of Shaanxi (Program No. 2021JQ-084).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The relevant information of numerical calculation has been reflected in the manuscript in the form of pictures, charts, etc.

**Acknowledgments:** The authors would like to acknowledge the support of National Natural Science Foundation of China (Grant No. 12202363) and the support of Natural Science Basic Research Program of Shaanxi (Program No. 2021JQ-084).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**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.
