Aerodynamic study of the corrugated airfoil at ultra-low Reynolds number

Abstract

In this study, Aeshna cyanea dragonfly forewing mid-cross-section corrugated airfoil was simulated at ultra-low Reynolds number. The corrugated airfoil was compared with its smooth counterpart to study the effect of the corrugations upon the aerodynamic performance. Unsteady two-dimensional laminar flow was solved using FLUENT. This study was divided into gliding phase and flapping phase. In the gliding phase, the corrugated airfoil produced a higher lift force with respect to the profiled airfoil at both tested Reynolds numbers ( $1400$ , $200$ ) with comparable drag coefficient for all the tested angles of attack. In the flapping phase, both the corrugated airfoil and the flat-plate have a very similar flow behavior which yields a very similar aerodynamic performance at Re $=1400$ . A structural analysis was performed to compare the corrugated airfoil with the flat-plate. The analysis revealed the superiority of the corrugated airfoil over the flat-plate in decreasing the deflection under the applied load. The reduced frequency was varied to study its impact on the aerodynamic performance. By increasing the reduced frequency, the thrust and the lift forces increased by $82$ % and $75$ %, respectively. Any increase in the reduced frequency will increase lift and thrust forces, but the propulsive efficiency will deteriorate.

Introduction

In the 1990s, the Defense Advanced Research Projects Agency (DARPA) initiated the challenge to create small and lightweight flyers to mimic nature. Consequently, the micro aerial vehicles (MAVs) became a research frontier attracting the researchers all over the world.^1,2 The MAVs concept is based on being small and inexpensive unmanned aerial vehicles comparable in size to biological flyers which put them in a new category of aviation aircraft.^3,4

The aerodynamic challenge is represented by the low Reynolds number $(\mathrm{Re})$ where the MAVs fly. The Reynolds number can be defined as the ratio between the inertial forces to the viscous forces. Typically, for most of MAVs, the small chord length and the relatively low flight speed made them fly in a Reynolds number range from ${10}^2$ to ${10}^4$ . At this low Reynolds number, the viscous forces significantly dominate the flow. The aerodynamics at this range of Reynolds number is one of the most interesting and less understood aspects. ⁵ Based upon lift and thrust forces generation, we can classify the MAVs into fixed wings, rotary wings, and flapping wings. ⁶

In nature, the greatest optimizer, there are no examples of fixed wings or even rotary wings. Nature gave flapping wings the chance to be evolved over millions of years. All birds and insects have well-defined wings and fuselage that blend together. By flapping, birds and insects can produce both lift and thrust forces to make extreme maneuvers, fly effortlessly, and sustain flight for long durations. ⁶ Flapping wings are nature’s answer for the perfect MAVs. ⁷ The interest in flapping wings field has increased rapidly over the past few decades. Hence, a lot of efforts have been made to enhance our understanding of flapping aerodynamics. The more we understand, the better our MAVs design becomes. Insects such as dragonfly, locust, and hawkmoth represent a frontier for researches due to their unique wing shape and kinematics.

Several studies were concerned with dragonfly experimentally and numerically.^8–19 Dragonfly is basically a flapping insect which can flap its wing with approximately 40 beats/s.^20–23 The wingbeat frequency was found to be size-dependent as small species tend to use higher frequency than large species. ²¹ However, dragonflies can be considered to be also gliders. The Aeshna cyanea dragonfly species can glide up to $30\mathrm{s}$ without any noticeable loss in altitude. It has been hypothesized that dragonflies adopt the gliding flight mode in order to save energy during flapping flight. The gliding mode can be also used for cooling dragonfly wings during hot weather as they can run under the risk of overheating during different flapping flights.^11,24 The typical Reynolds number for the dragonfly can be within the range of ${10}^2<\mathrm{Re}<2400$ , which can be classified as ultra-low Reynolds number.

Okamoto et al. ⁹ were one of the very first researchers who studied the corrugated shape of the airfoil and the aerodynamic characteristics of a real dragonfly’s wings and body through three-dimensional (3D) experiments. All the aerodynamic characteristics have been studied within the Reynolds number range ${10}^3<\mathrm{Re}<{10}^4$ . The corrugations were found to have a slight effect on the lift slope curve and the minimum drag coefficient but have a great effect on the maximum lift coefficient and the drag coefficient at high angles of attack. The downward leading-edge was found to perform better than the upward leading-edge, which contradicts the results of Kesel. ¹⁰ The maximum lift coefficient was found to approach unity for the tested dragonfly wings, which was higher than that produced by the artificial airfoils. However, the minimum drag coefficient was higher. They referred these higher values to the high wing stiffness caused by the zigzag veins. Another two-dimensional (2D) experimental analysis was made by Kesel ¹⁰ to investigate the aerodynamic performance of a dragonfly forewing. The corrugated airfoils correspond to three different positions along the dragonfly wingspan. The tested Reynolds numbers were higher than nature reality, but low enough to treat the flow as a laminar flow. Due to the corrugations’ presence along the chord and the span, the air became trapped inside the valleys in the form of rotating vortices. These vortices produce a net negative pressure on the airfoil surfaces regardless the orientation of the valleys, which was agreed by Vargas et al., ¹¹ New et al., ²⁵ and Kim et al. ²⁴ and confirmed Rees ¹⁹ suggestion. Kesel ¹⁰ found that the leading-edge orientation has a minor effect on the lift force produced. In general, the corrugated airfoils showed a relatively low drag force with higher lift force compared to the flat-plate and the corrugated smooth counterpart. Kesel ¹⁰ concluded that neither a curved plate nor a flat-plate would fulfill the static demand of an insect wings unless they become thicker, which will increase the wings’ weight. The corrugations are very important in resisting bending loads, insect stability, and aerodynamic performance.

Vargas et al. ¹¹ used the immersed boundary method with non-body conformal Cartesian grids to solve the time-dependent, incompressible, viscous Navier–Stokes equations. The studied corrugated airfoil corresponds to profile 2 in Kesel. ¹⁰ They were able to clearly observe the trapped vortices inside the corrugations, which cause the airflow stream to be similar to the profiled airfoil. The trapped vortices produced a negative shear drag; therefore, the contribution of the shear drag in the total drag for the corrugated airfoil was lower than the other two tested airfoils. The corrugated airfoil was found to produce a comparable and even sometimes higher lift force with a comparable drag force to the profiled airfoil. Furthermore, Vargas et al. ¹¹ concluded that the corrugated airfoil has an aerodynamic performance benefit compared to the conventional airfoils. It provides a stiff structure with lightweight which can resist bending loads. They suggested that such a wing can be useful in the future MAVs’ wing design. Kim et al. ²⁴ conducted a 2D numerical analysis to study the effects of the corrugations upon the drag and the lift coefficients. Various corrugation shapes were solved using FLUENT solver at different Reynolds numbers for different angles of attack. They were able to capture the generated vortices inside the corrugation valleys which gave the dynamic smoothing effect. The numerical results showed that the corrugations found in the dragonfly forewing increased the lift coefficient for all the tested angles of attack compared to its smoothed version at $\mathrm{Re}=1400$ . In general, they found that profile 1 outperforms its smooth counterpart for all the tested Reynolds numbers and angles of attack. It was noted that the maximum lift coefficient increased as the Reynolds number increased. Kim et al. ²⁴ observed that the production of the lift force is mainly due to the pressure-side corrugations. Moreover, they reported the importance of the pressure side’s middle wide corrugation in the lift force production. It was noted that as the corrugations became shallower, both the lift and the drag forces were reduced. Another 2D numerical analysis was conducted to study the effect of the stroke plane angle for a flapping dragonfly corrugated airfoil upon thrust and lift generation. ²² FLUENT solver was used to solve the corrugated airfoil which corresponds to profile 2 in Kesel ¹⁰ in simple plunging and pitching motions (i.e. forward flight) within range of ${10}^3<\mathrm{Re}<5\times {10}^3$ . The simulation results showed that the difference between the tested corrugated airfoil and the flat-plate was very small based on the used flapping frequency range. They found that the thrust coefficient reaches its peak value at $\mathrm{Re}=3\times {10}^3$ and began to decrease as the Reynolds number increases. Moreover, as the Reynolds number was even higher, the thrust coefficient trend exhibits one peak in the upstroke instead of two peaks and starts to diminish. Chuang et al. ²² found that as the stroke plane angle $\mathrm{\beta }$ increased from $45°$ to $60°$ , the twin peak trend was only visible at $\mathrm{Re}<5\times {10}^3$ . The lift coefficient variation was found to be vital during the downstroke for $\mathrm{Re}<5\times {10}^3$ . Moreover, Chuang et al. ²² concluded that the increase in the stroke inclined angle and the Reynolds number has a positive effect on the lift coefficient. To gain maximum thrust force, we shall use $\mathrm{Re}=3\times {10}^3$ regardless of the stroke plane angle value. Another 2D simulations using FLUENT solver were carried out by Lavimi et al. ⁸ to study two different dragonfly’s cross sections in a flapping motion at Reynolds number of $5000$ and ${10}^4$ . They found that the airfoils with corrugated surface have a better aerodynamic performance than the flat-plate. They reported a rapid increase in the lift force caused by the pitching-up rotation. Moreover, the formations of new strong vorticity layers were responsible for the increase in the lift coefficient, not the stall absences. The corrugated airfoil caused the separation to be delayed in the pitching motion; hence, the lift coefficient was further increased. Lavimi et al. ⁸ found that by reducing the Reynolds number, strong vorticity layers were formed for all the tested airfoils during the flapping motion. These strong vorticity layers besides the pitching motion caused the lift force to increase. Finally, the second tested cross section was found to be a suitable alternative for the flat-plate in designing MAVs.

Based on the above-mentioned survey, it is of great importance to study the aerodynamic performance of the corrugated airfoil during different flight phases. This study is going to be as a step for a better understanding of nature for better MAVs performance. It is clear that the performed 2D analysis cannot predict the 3D effects of the actual dragonfly wing in nature. However, it gives a clear view about the flow behavior among the corrugations and the effect of the trapped vortices. On the current research, the corrugated airfoil based upon A. cyanea dragonfly forward wing mid-cross-section will be considered as the main airfoil. The corrugated airfoil will be analyzed in two different flight phases (i.e. gliding and flapping). For the gliding phase, different parameters will be varied to study their effect on the airflow stream, hence the aerodynamic performance. The results will be compared to the flat-plate as a base airfoil. For the flapping phase, the effect of reduced frequency variation will be considered to study its impact on the aerodynamic forces and the propulsive efficiency.

Numerical method

Numerical settings

The flow motion around the tested airfoils is governed by the so-called Navier–Stokes equations which consist of conservation of mass, momentum, and energy. The energy equation is ignored since the flow is assumed to be isothermal. In this study, FLUENT solver was used to solve the 2D, unsteady, incompressible, viscous flow. The flow is assumed to be laminar; hence, no turbulent models were needed to be coupled with the governing equations.

FLUENT is a finite volume solver that turns the Navier–Stokes equations into a set of algebraic equations to be solved numerically. The pressure-based solver was used since the flow is incompressible. Two algorithms can be found under the pressure-based solver, namely, segregated solver and coupled solver. In the segregated solver, FLUENT solves the pressure correction (i.e. the continuity equation) and the momentum equations sequentially. In the coupled solver, FLUENT solves the pressure correction and the momentum equations simultaneously. ²⁶ The pressure-based coupled solver was suitable for our calculations since we only deal with a single-phase flow (namely, air) with low speed.

The coupled solver requires more memory than the segregated solver; however, it produces more accurate solutions. The solutions of the discrete governing equations from FLUENT are stored at the center of the control volumes, so interpolation is needed to transfer the solutions to the control volume faces. For the gradient interpolation, the Least Squares Cell Based was used. Due to the presence of swirling flows trapped inside the corrugations’ valleys and the possible vortex shedding, it was recommended by FLUENT manual ²⁶ to use the Quadratic Upwind Interpolation (QUICK) scheme in the momentum interpolation. The Pressure Staggering Option (PRESTO) scheme was used for the pressure to accurately predict the high gradient over the corrugated airfoil surfaces.

The solutions obtained from FLUENT solver must satisfy our convergence criteria conditions, ²⁶ so we can clearly adopt the solution. The convergence criteria conditions are as follows:

• All the Navier–Stokes equations’ residuals must achieve the limit of $({10}^{-5})$ for every time step in case of unsteady flow or to the limit of $({10}^{-7})$ for steady flow.

• The overall mass imbalance must be less $(1\%)$ than the smallest flux through the domain boundary.

• The appearance of the requested phenomena (i.e. vortex shedding in our case).

The aerodynamic forces (i.e. lift and drag) are computed from the pressure and the shear stress integrations over the tested airfoil surfaces. Their corresponding coefficients are defined by

$${\mathrm{C}}_{\mathrm{l}}=\frac{\mathrm{L}}{\frac{1}{2}{\mathrm{\rho }\mathrm{U}}_{\infty }^2{\mathrm{S}}_{\mathrm{o}}}$$

(1)

$${\mathrm{C}}_{\mathrm{d}}=\frac{\mathrm{D}}{\frac{1}{2}{\mathrm{\rho }\mathrm{U}}_{\infty }^2{\mathrm{S}}_{\mathrm{o}}}$$

(2)

where the lift and the drag forces are denoted by $\mathrm{L}$ and $\mathrm{D}$ , respectively. ${\mathrm{U}}_{\infty }$ denotes the free stream velocity. ${\mathrm{S}}_{\mathrm{o}}$ is the wing area which equals the wing span times the mean chord length. In this study, the wing area is taken to be the mean chord length (i.e. per unit depth).

Airfoil geometry

Okamoto et al. ⁹ and Kesel ¹⁰ obtained different cross sections by slicing the dragonfly forewing into several elements. Many researchers including Vargas et al. ¹¹ used the same geometry. Vargas et al. ¹¹ constructed their study upon the mid-cross-section of the A. cyanea dragonfly forewing (profile 2 in Kesel ¹⁰ ). The corrugated airfoil used in this study is based on Vargas et al. ¹¹ The corrugated airfoil is shown in Figure 1(a).

Figure 1.

Tested airfoils: (a) corrugated airfoil, (b) profiled airfoil, and (c) flat-plate.

A lot of trials were carried out to get the points of the corrugated airfoil. The points were extracted digitally, but the corrugated shape was found always to lack the accuracy needed for the solution. The thickness of the airfoil and the corrugations shape (i.e. the sequence of edges and bends) were found to influence the solution significantly, ¹⁰ so the accuracy of the shape was of great importance. Finally, the airfoil points were taken from Vargas ²⁷ PhD thesis. The airfoil edges were kept sharp without smoothing as it was reported that the blunt corners and sharp corners experience almost the same results. ²⁸

All insects in nature with their different species were found to fly at ultra-low Reynolds number (i.e. $\mathrm{Re}<{10}^4$ ). So, to mimic nature as much as we can, the airfoils were solved at ultra-low Reynolds number. The corrugated airfoil and the corresponding profiled airfoil were compared together. The profiled airfoil was obtained by connecting the corner intersection points of the corrugated airfoil together to have a streamlined surface as shown in Figure 1(b). The process of connecting points was based on the solution of the corrugated airfoil at $\mathrm{\alpha }=0°$ and $\mathrm{Re}=1400$ . The solution of the corrugated airfoil showed that the air became trapped inside the valleys which gave a smooth streamline surface ¹¹ (i.e. profiled airfoil). A flat-plate in Figure 1(c) was solved with a constant thickness of $2.7\mathrm{mm}$ and flat edges on both the leading-edge and the trailing-edge as a base airfoil in the comparison.

The tested airfoils were solved under the same conditions. The study parameters were the same used by Vargas et al. ¹¹ The mean chord length $({\mathrm{c}}_{\mathrm{m}})$ of all the tested airfoils for the gliding phase was set to be $80\mathrm{mm}$ . Air properties are $\mathrm{\rho }=1.225\mathrm{kg}/{\mathrm{m}}^3$ and $\mu =1.7894\times {10}^{-5}\mathrm{kg}/(\mathrm{m}\mathrm{s})$ .

In seeking to complete the aerodynamic study, the flapping phase of the corrugated airfoil was necessary to be investigated. During a flapping motion, an effective angle of attack $({\mathrm{\alpha }}_{\mathrm{e}})$ was generated due to the flow wash out phenomena illustrated in Figure 2. The effective angle of attack can be calculated according to equation (3)

$${\mathrm{\alpha }}_{\mathrm{e}}\left(\mathrm{t}\right)=\mathrm{ta}{\mathrm{n}}^{-1}\left(\frac{-\mathrm{y}\prime \left(\mathrm{t}\right)}{{\mathrm{U}}_{\infty }}\right)-\mathrm{\theta }\left(\mathrm{t}\right)$$

(3)

where $\mathrm{y}\prime (\mathrm{t})$ is the velocity derivative with respect to time. $\mathrm{\theta }(\mathrm{t})$ is the instantaneous pitching angle of the airfoil. All the angles are measured with respect to the negative x-direction as the reference axis. In general, a complete flapping cycle consists of four main motions. The motions are the downstroke translation motion and the upstroke translation motion, besides the supination rotation (i.e. rotation of the wing at the end of the upstroke) and the pronation rotation (i.e. rotation of the wing at the end of the downstroke). During the supination and the pronation rotations, the angle of attack changes rapidly to ensure that the leading-edge always faces the airflow.

Figure 2.

Definition of the effective angle of attack. ²⁹

For all the flapping tested cases, the airfoil center of rotation was made at quarter the mean chord length $({\mathrm{c}}_{\mathrm{m}}/4)$ . The numerical investigation was performed at a fixed Reynolds number of $1400$ and a mean chord length of $10\mathrm{mm}$ .

Equations (4) and (5) are used to describe the coupled corrugated airfoil plunging and pitching motions

$$\mathrm{y}\left(\mathrm{t}\right)={\mathrm{y}}_{\mathrm{o}}\sin (2\mathrm{\pi }f\mathrm{t}+\mathrm{\phi })$$

(4)

$$\mathrm{\theta }\left(\mathrm{t}\right)={\mathrm{\theta }}_{\mathrm{o}}\sin \left(2\mathrm{\pi }f\mathrm{t}\right)$$

(5)

where ${\mathrm{y}}_{\mathrm{o}}$ is the plunging amplitude while ${\mathrm{\theta }}_{\mathrm{o}}$ is the pitching amplitude. $f$ denotes the motion frequency while $\mathrm{\phi }$ is the phase difference between the plunging motion and the pitching motion. Varying the phase difference yields advanced motion, synchronized motion, and delayed motion.

To be able to understand the aerodynamics behind the flapping mechanisms, several non-dimensional parameters have to be considered, one of them is the reduced frequency. The reduced frequency $(\mathrm{k})$ is one of the most important non-dimensional parameter used to describe the vortex shedding formed by the unsteady flapping airfoils. It can be defined to be the ratio of two characteristic velocities,^30,31 namely, the flapping speed and the free stream speed, equation (6)

$$\mathrm{k}=\frac{2\mathrm{\pi }f{\mathrm{c}}_{\mathrm{m}}}{{\mathrm{U}}_{\infty }}$$

(6)

The product of the frequency by the mean chord length (i.e. $f{\mathrm{c}}_{\mathrm{m}}$ ) is the flapping speed.

Power coefficient and thrust efficiency calculations

One of the current study parameters of interest is the power $({\mathrm{P}}_{\mathrm{o}})$ . The instantaneous input power is the amount of energy needed by the airfoil to overcome the aerodynamic forces opposing the motion during a steady flapping flight. The input power can be expressed as ³²

$${\mathrm{P}}_{\mathrm{o}}\left(\mathrm{t}\right)=-\mathrm{L}\left(\mathrm{t}\right)\times \frac{\mathrm{dy}\left(\mathrm{t}\right)}{\mathrm{dt}}-{\mathrm{M}}_{\mathrm{a}}\left(\mathrm{t}\right)\times \frac{\mathrm{d}\mathrm{\theta }\left(\mathrm{t}\right)}{\mathrm{dt}}$$

(7)

where $\mathrm{y}(\mathrm{t})$ is the vertical motion of the airfoil and $\mathrm{\theta }(\mathrm{t})$ is the pitching motion of the airfoil about the airfoil quarter chord. ${\mathrm{M}}_{\mathrm{a}}(\mathrm{t})$ is the moment created by the aerodynamic forces at the airfoil quarter chord. FLUENT solver is used to calculate the moment coefficient with respect to the origin for each time step. The moment has to be transferred to the airfoil center of rotation (i.e. airfoil quarter chord) using the following equation ³³

$${\mathrm{M}}_{\mathrm{a}}\left(\mathrm{t}\right)={\mathrm{M}}_{\mathrm{o}}\left(\mathrm{t}\right)+{\mathbf{R}}_{\mathbf{ao}}\times \mathrm{L}\left(\mathrm{t}\right)$$

(8)

where ${\mathrm{M}}_{\mathrm{o}}(\mathrm{t})$ is the moment with respect to the origin and ${\mathbf{R}}_{\mathbf{ao}}$ is the vector from the origin to the airfoil center of rotation. The simulation starts with the airfoil in its top position, which means that the position of the airfoil is given by ³³

$${\mathbf{R}}_{\mathbf{ao}}=-{\mathrm{y}}_{\mathrm{o}}[1-\sin (2\mathrm{\pi }f\mathrm{t}+\mathrm{\phi })]$$

(9)

The moment can be non-dimensionalized as

$${\mathrm{C}}_{\mathrm{m}}=\frac{{\mathrm{M}}_{\mathrm{o}}}{\frac{1}{2}{\mathrm{\rho }\mathrm{U}}_{\infty }^2{{\mathrm{c}}_{\mathrm{m}}}^2}$$

(10)

The input power can be also non-dimensionalized as follows

$${\mathrm{C}}_{\mathrm{power}}=\frac{{\mathrm{P}}_{\mathrm{o}}}{\frac{1}{2}{\mathrm{\rho }\mathrm{U}}_{\infty }^3{\mathrm{c}}_{\mathrm{m}}}$$

(11)

The propulsive efficiency is a measure of the energy lost in the wake versus the energy used in creating the necessary thrust, which can be expressed as ³²

$${\mathrm{\eta }}_{\mathrm{prop}}=\frac{{\mathrm{C}}_{\mathrm{t}}}{{\mathrm{C}}_{\mathrm{power}}}$$

(12)

where the thrust coefficient is defined as^33,34

$${\mathrm{C}}_{\mathrm{t}}=-{\mathrm{C}}_{\mathrm{d}}$$

(13)

The definition of the propulsive efficiency loses its meaning when the airfoil produces drag instead of thrust.

Computational domain

In this study, we nominated the computational domain to be a squared domain with structured elements to simulate reasonably the corrugated shape of the main airfoil. The schematic diagram of the domain is depicted in Figure 3.

Figure 3.

Computational domain.

The computational domain consists of four boundaries. All the boundaries were far enough to eliminate their effects on the geometry wall (i.e. airfoils). In the present gliding analysis, we were interested in a narrow range of angles of attack. Hence, the domain size was chosen to be relatively small. The upper and the lower boundaries were made at a distance of $5{\mathrm{c}}_{\mathrm{m}}$ from the airfoil; same was the upstream from the airfoil leading-edge. The downstream boundary is located $12{\mathrm{c}}_{\mathrm{m}}$ from the airfoil trailing-edge, where ${\mathrm{c}}_{\mathrm{m}}$ represents the mean airfoil chord length. The velocity inlet condition was applied to the upstream, the upper, and the lower boundaries, which corresponds to a Dirichlet boundary condition. The downstream boundary was set to be pressure outlet, which corresponds to a Neumann boundary condition. The airfoil geometry was set to be wall (no-slip condition).

For the flapping phase, all the independent studies’ simulations were conducted for the combined pitching and plunging motions of the corrugated airfoil at a flapping frequency of $25\mathrm{Hz}$ . The domain size in Figure 3 was first used; after then, the domain was enlarged step by step till the error between the new domain and the old domain was within acceptable range $(<8\%)$ . Hence, the new domain is suitable for the structured mesh quality and our computational resources.

Figure 4(a) and (b) illustrates one complete cycle for the forces coefficient values evaluated for three different domain sizes. The first number in the legend represents the upper, the lower, and the upstream distances, while the second number represents the downstream distance. The time axis was normalized for one complete flapping cycle. Finally, based on the obtained results, the domain size was updated so that the upper and the lower boundaries were positioned at a distance $12{\mathrm{c}}_{\mathrm{m}}$ from the maximum point the airfoil will reach during its motion. The upstream was located at a distance $12{\mathrm{c}}_{\mathrm{m}}$ from the airfoil leading-edge. The downstream was extended up to $24{\mathrm{c}}_{\mathrm{m}}$ far from the airfoil trailing-edge.

Figure 4.

Corrugated airfoil forces coefficient in a flapping motion for different domain sizes: (a) lift coefficient and (b) drag coefficient.

Mesh generation

Due to the geometry complexity, ICEM-CFD was needed. ICEM-CFD represents a powerful tool to control every single point along the domain which gave the ability to adjust accurately the position of the nodes along the chosen airfoils’ edges. The mesh consists of three grid regions with different cell sizes. The first region surrounds the airfoil with clustered mesh to capture the vortices with high fidelity. The growth rate was set to be $1.01$ . Care was taken while meshing to keep the maximum aspect ratio as minimum as possible and the minimum Eriksson Skewness above $0.7$ . Finally, a high number of nodes were always used along the computational domain to ensure high grid density around the studied airfoils as illustrated in Figure 5.

Figure 5.

Mesh around the corrugated airfoil.

For the flapping phase, a mesh-independent study was carried out to ensure that the mesh resolution was enough to capture the solution details (i.e. vortex shedding phenomena). The results shown in Figure 6(a) and (b) represent the forces coefficient values for one complete cycle at different mesh sizes, namely, fine mesh ( $\approx 1100000$ element), medium mesh ( $\approx 820000$ element), and course mesh ( $\approx 680000$ element). The results obtained indicated that the error between the medium mesh size and the fine mesh size is within the acceptable limit $(<5\%)$ . As a result, the medium mesh size was chosen to be used in our current simulations in order to save the computer CPU time.

Figure 6.

Corrugated airfoil forces coefficient in a flapping motion for different mesh sizes: (a) lift coefficient and (b) drag coefficient.

The need for dynamic mesh rises from the requirement of the airfoil to flap in a certain prescribed motion. The mesh must be able to absorb the deformation performed by the motion and simultaneously track the body. To be able to use dynamic mesh model in FLUENT, we must construct the very first volume mesh and provide a description of the boundary motion. FLUENT provides either boundary profiles, user-defined functions (UDFs), or 6-degree-of-freedom (6-DOF) solver to describe the motion of the boundary. ²⁶ In this study, the UDF was used to describe the boundary motion. The UDF is basically a C code with a lot of varieties of pre-written macros which can be linked to FLUENT. For FLUENT to recognize the UDF files, they must be compiled or either interpreted. All the dynamic meshing macros require compiling. ³⁵ The UDF macro used in this study based on the literature is the $\mathrm{DEFINE}\_\mathrm{CG}\_\mathrm{MOTION}$ macro. It is important to get the proper combination of the dynamic mesh methods to avoid the negative cell volume errors. Some of the options can be employed in only certain applications and several can be selected together.^26,36 For our application, the smoothing and the re-meshing methods were used together as they both produce better quality mesh and enable larger time steps. FLUENT uses all the activated dynamic mesh settings together every time step to optimize the updated mesh. Based on the enabled methods, FLUENT decides the best for the deforming zone and then applies the smoothing techniques followed by re-meshing the necessary zones. ³⁶

Time-independent study

For the gliding phase, the time step of $1\times {10}^{-4}\mathrm{s}$ was found to be sufficient for the current gliding analysis. ²⁴ However, a time-independent study was needed for the flapping phase to ensure the capability of the used time step in capturing the vortex developing and shedding phenomena accurately. To calculate the required time step, two rules of thumbs were used. The used time step must be less than the time step obtained from Courant number relation and always less than the minimum mesh length divided by the air free stream velocity. Figure 7 illustrates the conducted study which shows $\mathrm{\Delta }\mathrm{t}=1.5\times {10}^{-5}\mathrm{s}$ being sufficient enough for the present numerical investigation.

Figure 7.

Time-independent study for the corrugated airfoil in a flapping motion.

Validation

To be able to solve our main gliding and flapping problems, a validation case was required to validate the adapted mesh strategy and the used FLUENT setup. The validation problem can be summarized as a flat-plate with $30\mathrm{mm}$ chord length and a thickness of $1.5\mathrm{mm}$ that oscillates in a pure pitching motion at $\mathrm{Re}=8000$ . The airfoil starts the motion from rest with a pitch-up motion followed by a pitch-down motion at a reduced frequency of $0.31$ according to the reduced frequency definition in Okamoto and Azuma. ³⁷ The instantaneous angle of attack of the flat-plate can be described using equation (14)

$$\mathrm{\theta }\left(\mathrm{t}\right)={\mathrm{\theta }}_{\mathrm{o}}\sin \left(2\mathrm{\pi }f\mathrm{t}\right)$$

(14)

The pitching amplitude in our case equals $16°$ . The same domain and boundary conditions used in Lavimi et al. ⁸ were applied. The problem was found to be turbulent so, the Reynolds stress turbulence model with linear pressure-strain was used. The pressure-based SIMPLEC solver was found to be suitable in the present validation problem. Figure 8(a) and (b) shows the obtained lift and drag coefficients compared to the literature.^8,12,37 The trend of the present numerical study and the experimental results are close to each other. The hysteresis loop width was found to be higher in the present numerical study compared to Okamoto and Azuma. ³⁷ This difference can be related to the fact that the experiment was 3D, not 2D. The 3D flow along the span may affect the stability and the formation of the leading-edge vortex (LEV). Another contribution to the deviation between the present numerical study and the experiment is the low Reynolds number. The flow may not be fully turbulent and the boundary layer at the leading-edge could be laminar. ⁸ From the obtained results, we can say that the present numerical results exhibit an acceptable matching with the results found in the literature.

Figure 8.

Aerodynamic coefficients versus the instantaneous AOA: (a) lift coefficient and (b) drag coefficient.

Results and discussion

Gliding phase

During all the performed simulations, it was noted that the flow behaves as a steady flow at $\mathrm{Re}=200$ for all the tested angles of attack and behaves as unsteady flow beyond $\mathrm{\alpha }=15°$ . Moreover, at $\mathrm{Re}=1400$ and $\mathrm{\alpha }=0°$ , the flow was fully steady. As the angle of attack increases up to $\mathrm{\alpha }=5°$ , the unsteady flow pattern appears slightly. Furthermore, the flow was found to be fully unsteady beyond $\mathrm{\alpha }=5°$ . By increasing the Reynolds number up to ${10}^4$ , the flow becomes unsteady for all the angles of attack. The main reason behind the unsteadiness was the vortex field change around the airfoils. The change of the vortex field causes the lift and the drag forces’ fluctuation to increase by increasing the angle of attack.

Figure 9 illustrates the flow streamlines at the maximum lift coefficient for the corrugated airfoil at angles of attack of $0°,5°,10°$ , and $15°$ with $\mathrm{Re}=1400$ . Figure 9(a) clearly demonstrates the appearance of trapped recirculating flow inside the corrugated valleys at $\mathrm{\alpha }=0°$ . The trapped vortices gave the effect of dynamic smoothing of the flow streamlines around the corrugated airfoil to simulate the profiled airfoil in Figure 1(b). The trapped vortices generate a negative shear drag inside the valleys which helps in reducing the total corrugated airfoil drag force. Thus, the corrugated airfoil drag force became comparable and even sometimes lower than the profiled airfoil. As the angle of attack increased, the vortex along the corrugated airfoil upper surface increased in size till it occupied approximately half the corrugated airfoil chord length as clarified in Figure 9(d). Clearly, we can see the different stages on which the vortices on the airfoil upper surface became larger and began to attach to each other forming two large vortices as depicted in Figure 9(d). One of the vortices was found near the leading-edge, while the other one was found to occupy the airfoil aft portion. However, the trapped vortices can be still observed on the lower side of the airfoil and did not completely disappear.

Figure 9.

Streamlines of the corrugated airfoil at Re = 1400: (a) 0°, (b) 5°, (c) 10°, and (d) 15°.

The average lift and drag coefficients for the tested airfoils are depicted in Figures 10 and 11 when the Reynolds number equals $200$ and $1400$ , respectively. Clearly, we can see that the corrugated airfoil exhibits the highest lift forces among all the tested angles of attack compared to the profiled airfoil. These results indicate that the corrugations play an important role in producing lift force even if there is no recirculating flow inside its valleys.

Figure 10.

Coefficient values at Re = 200: (a) lift and drag coefficients versus AOA and (b) polar diagram.

Figure 11.

Coefficient values at Re = 1400: (a) lift and drag coefficients versus AOA and (b) polar diagram.

Figure 10 illustrates that the flat-plate produced almost the same lift force with respect to the corrugated airfoil with approximately the same drag coefficient. However, Figure 11 illustrates the superior performance of the flat-plate at $\mathrm{\alpha }=5°$ and $10°$ but with higher drag coefficient. At higher angles of attack, the flat-plate performance deteriorates and the corrugated airfoil performs better. The drag coefficient for all the tested airfoils was found to experience an increase in the values by decreasing the Reynolds number. The reason is that at lower Reynolds number values, the skin friction caused by the airflow viscosity becomes the main contributor to the drag coefficient. ¹¹

Figure 12(a) and (b) illustrates the gliding ratio for the tested airfoils at $\mathrm{Re}=200$ and $1400$ , respectively, with various angles of attack. The gliding ratio represents the wing aerodynamic performance during gliding flight (unpowered flight), which can be expressed as ${\mathrm{C}}_{\mathrm{l}}/{\mathrm{C}}_{\mathrm{d}}$ . The corrugated airfoil was found to exceed the profiled airfoil performance for the two tested Reynolds numbers for all the tested angles of attack. Compared to the flat-plate, the corrugated airfoil performs better at lower Reynolds number. However, at $\mathrm{Re}=1400$ and $\mathrm{\alpha }=5°$ , the flat-plate was slightly higher than the corrugated airfoil, but soon the flat-plate performance deteriorated rapidly to achieve the lowest value among the tested airfoils.

Figure 12.

Gliding ratio (C_l/C_d) versus AOA: (a) Re = 200 and (b) Re = 1400.

Figure 13 compares the static pressure coefficient $({\mathrm{C}}_{\mathrm{p}})$ distribution for the tested airfoils at $\mathrm{Re}=1400$ and $\mathrm{\alpha }=0°$ . Figure 13 demonstrates the role of corrugations in producing lift force higher than the profiled airfoil. At each corrugation valley, a sharp edge with high amplitude with respect to the negative sign was clearly observed; after then, the value began to decrease till the beginning of the next valley. The high amplitudes’ values are mainly due to the core strength of the recirculating flow inside each valley. The amplitude of the sharp edge was found to change from one valley to the other based on the strength of the recirculating flow. For the wide valley in the suction side, both airfoils have approximately the same values of pressure coefficient. However for the same wide valley in the pressure side, an increase in the pressure coefficient values was found mainly due to the trapped vortex. For the flat-plate, we can see that the pressure coefficients for the pressure side and the suction side are identical, which is clearly rational for symmetric airfoil at zero angle of attack. The flat-plate experienced the highest peak value at the leading-edge among the other tested airfoils, which can be attributed to the thick flat-plate leading-edge.

Figure 13.

Static pressure coefficient distribution for the tested airfoils at α = 0° and Re = 1400.

Figure 14 illustrates the streamlines formed around the tested airfoils at fixed angle of attack of $5°$ for $\mathrm{Re}=200$ and $1400$ , respectively. At $\mathrm{Re}=200$ , the flow sticks to the profiled airfoil surfaces with no sign for any separation reducing the pressure drag force coefficient of the airfoil but in-turn increasing the viscous drag force coefficient. Similar to the profiled airfoil, the flow was found to be fully attached to the upper and the lower surfaces of the flat-plate. Owing to the simple symmetric shape of the flat-plate (not streamlined shape), the pressure drag force coefficient was found to be higher than the profiled airfoil with lower value of the viscous drag force coefficient. For the corrugated airfoil, small trapped vortices can be noticed to be formed inside the valleys of the corrugated airfoil. Yet, strong enough to reduce the viscous drag force coefficient compared to the profiled airfoil. The corrugated airfoil was found to experience a higher pressure drag force coefficient almost two times the profiled airfoil value due to the corrugated airfoil shape. At $\mathrm{Re}=1400$ , moderate flow separation along all the tested airfoils can be clearly observed. The separated recirculating flow has increased the pressure drag force coefficient of the corrugated airfoil to be twice the profiled airfoil value. Yet, it decreased the viscous drag force coefficient to be half the value of the profiled airfoil. The area occupied by the generated recirculating flow can be seen to be bigger in case of the flat-plate compared to the profiled airfoil. Although the recirculating flow helps the flat-plate to produce less viscous drag force coefficient compared to the profiled airfoil, it produced higher pressure drag force coefficient.

Figure 14.

Streamlines of the tested airfoils at α = 5°: (a) Re = 200 and (b) Re = 1400.

Table 1 presents the components of the drag coefficient for the tested airfoils at $\mathrm{\alpha }=5°$ for two different Reynolds numbers. The pressure drag force coefficient denoted by ${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{p}}}$ and the viscous drag force coefficient denoted by ${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{v}}}$ .

Table 1.

Drag force coefficient’s components for the tested airfoils at α = 5° for two different Reynolds numbers.

	$\mathrm{Re}=200$			$\mathrm{Re}=1400$
	${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{p}}}$	${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{v}}}$	${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{total}}}$	${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{p}}}$	${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{v}}}$	${\mathrm{C}}_{{\mathrm{d}}_{\mathrm{total}}}$
Corrugated airfoil	$0.11998$	$0.15703$	$0.27701$	$0.07286$	$0.02951$	$0.10237$
Profiled airfoil	$0.06350$	$0.21972$	$0.28322$	$0.03992$	$0.06400$	$0.10392$
Flat-plate	$0.08539$	$0.19502$	$0.28042$	$0.06663$	$0.03909$	$0.10572$

The streamlines analysis showed the importance of the corrugated airfoil in reducing the net drag coefficient compared to the streamed airfoil (i.e. profiled) and the flat-plate by reducing the viscous drag force coefficient. From the performed analysis, we can say that the corrugated airfoil always performs better than its smoothed version (i.e. profiled airfoil) at low Reynolds number, and sometimes its performance is comparable to the flat-plate. Hence, the profiled airfoil will be excluded from the flapping analysis.

Flapping phase

The effect of the frequency was studied by fixing the other flapping parameters. The Reynolds number, plunging amplitude, pitching amplitude, and phase difference were kept fixed at $1400,1{\mathrm{c}}_{\mathrm{m}},30°$ , and $90°$ , respectively. Three frequencies were chosen as $25,40$ , and $60\mathrm{Hz}$ to yield reduced frequencies of $0.768,1.229$ , and $1.8435$ , respectively.

Figure 15 shows the instantaneous static pressure contour for the corrugated airfoil and the flat-plate at $\mathrm{Re}=1400,\mathrm{k}=1.229,{\mathrm{y}}_{\mathrm{o}}={\mathrm{c}}_{\mathrm{m}},{\mathrm{\theta }}_{\mathrm{o}}=30°$ , and $\mathrm{\phi }=90°$ . Figure 15 is mainly divided into two columns; the left column represents the static pressure contour for the corrugated airfoil with simultaneous lift and thrust coefficients. The right column represents the static pressure contour for the flat-plate with simultaneous lift and thrust coefficients. Figure 15 represents the downstroke and the upstroke motions of the tested airfoils divided into several intervals. During the beginning of the downstroke motion in Figure 15(a), a residual vortex can be clearly seen to occupy the lower aft portion of the airfoil from the previous motion. The residual vortex results in decreasing the pressure on the lower surface more than the pressure on the upper surface causing the lift force to be negative and to produce drag force for both airfoils. As the airfoil starts to move in the downstroke, the residual vortex starts to detach from the airfoil lower surface and moves further to the trailing-edge. The movement of the residual vortex results in enhancing the lift and the thrust coefficients. Moreover, by increasing the effective angle of attack in Figure 15(b), the trailing-edge vortex (TEV) was sheded downstream to form the well-known von Kármán vortex street but in a reversed manner, which is called, the reversed von Kármán vortex street. The reversed von Kármán vortex street is responsible for generating thrust force. At the same time, we can notice a LEV starts to form on the upper surface of the airfoil. As the time passes, the LEV grows in size and moves further downstream the airfoil upper surface. However, the LEV is still attached to the leading-edge, which causes the lift and the thrust forces to reach their maximum values during the flapping cycle. In Figure 15(c), the airfoil began to decelerate as it approaches the end of the downstroke motion. Here, we can see that LEV began to be detached from the airfoil upper surface and moves downstream to the trailing-edge. The deceleration of the airfoil is coupled with a rapid change in the effective angle of attack as the airfoil is changing the leading-edge direction to face the airflow stream (i.e. pronation rotation). The pressure difference between the upper and the lower surfaces is now decreased, which causes the lift and the thrust forces to further decrease.

Figure 15.

Static pressure contours for the corrugated airfoil and the flat-plate at Re = 1400, f = 40 Hz, y_o = c_m, θ_o = 90°, and φ = 90°: (a) t ≈ 0/10 T, (b) t ≈ 2/10 T, (c) t ≈ 4/10 T, (d) t ≈ 6/10 T, and (e) t ≈ 8/10 T.

Figure 15(d) illustrates the movement of the TEV downstream to the trailing-edge as the airfoil starts to accelerate in the upstroke motion. At the same time, the LEV began to be formed on the lower surface of the airfoil. As the airfoil continues in the upstroke motion, the effective angle of attack increases causing the LEV to further increase. At the same time, the TEV is sheded to form the reversed von Kármán vortex street hence, improving the thrust coefficient. Figure 15(e) represents the airfoil upstroke acceleration phase. During this phase, the LEV began to move downstream of the lower airfoil surface; however, the LEV is still attached to the airfoil surface. The lift and the thrust coefficients reach one more time almost their maximum values during the flapping cycle with respect to the negative sign for the lift coefficient. Thus, the thrust force is a twin positive peak during the complete cycle with a higher value during the downstroke than the upstroke. The airfoil now began to decelerate as the end of the upstroke is being approached. The end of the upstroke is coupled with a rapid rotation of the airfoil, called the supination rotation.

The two airfoils were found to have similar static pressure contours. Obviously, the corrugations did not have much influence on the flow field around the airfoil rather than their main influence was to change the thrust and the lift forces. The averaged values of the lift, the thrust, and the power coefficients together with the propulsive efficiency for both the corrugated airfoil and the flat-plate can be found in Table 2. The coefficient values are very similar for the two airfoils. The lift coefficient was found to be higher than that for the flat-plate, which can be attributed to the curvature presence in the corrugated airfoil due to the corrugations along the upper and the lower surfaces. The propulsive efficiency was found to be slightly higher in case of the corrugated airfoil than that for the flat-plate.

Table 2.

Performance parameters for the corrugated airfoil and the flat-pate (Re = 1400, k = 1.229, y_o = c_m, θ_o = 30°, and φ = 90°).

	Corrugated airfoil	Flat-plate
Cl	$0.01727$	$0.01304$
Ct	$0.7456$	$0.7512$
Cpower	$1.297$	$1.364$
ηprop%	$57.44$	$55$

Figures 16 and 17 represent the thrust and the lift coefficient variations of the corrugated airfoil for one complete flapping cycle resulting from the pressure variation caused by the vortices formation and shedding for the different tested frequencies. During the first few flapping cycles, the flow field around the airfoil seems to be unstable due to the start of the movement which affects the forces coefficient values. Eventually, the flow began to be stable causing the forces coefficient values to reach a quasi-steady condition. The forces coefficient values were calculated to be the average over one complete flapping cycle. The forces coefficient values need to be set for at least four cycles to reach an acceptable error between the cycles.

Figure 16.

Corrugated airfoil thrust coefficient for various flapping frequencies.

Figure 17.

Corrugated airfoil lift coefficient for various flapping frequencies.

In contrast to the lift force, the thrust force over one flapping cycle is being generated in both the upstroke and the downstroke motions as shown in Figures 16 and 17. This can be explained as a result of the aerodynamic force presences normal to the relative airflow around the airfoil, which has a positive thrust component in both the upstroke and the downstroke motions. In Figure 16, the first peak corresponds to the downstroke and the second peak corresponds to the upstroke.

The peak formed during the upstroke was smaller than that formed during the downstroke, which can be due to the corrugations effect on the forces coefficient values. By increasing the flapping frequency, the aerodynamic forces increased, which can be attributed to the effective angle of attack increase. Figure 17 illustrates that for $f=25$ and $40\mathrm{Hz}$ , the airfoil was not able to produce a positive lift force or at the beginning of the flapping cycle. The lift force generated during the downstroke was higher than that generated during the upstroke, which resulted in producing an average positive lift force for all the tested flapping frequencies. It can be observed that by increasing the flapping frequencies, the peaks’ values of the aerodynamic forces were slightly shifted to the left.

The effects of the reduced frequency on the thrust and the lift coefficients at ${\mathrm{R}}_{\mathrm{e}}=1400$ are shown in Figure 18. Figure 18 illustrates that increasing the reduced frequency has a great impact on the thrust and the lift coefficients. Increasing the reduced frequency from $0.768$ to $1.229$ enhances the thrust and the lift coefficients values by approximately $82\%$ and $75\%$ , respectively. Furthermore, by increasing the reduced frequency to $1.8435$ , the thrust and the lift coefficients further increased by approximately $62\%$ and $95\%$ , respectively. Although we can find that increasing the reduced frequency has a significant impact on both the thrust and the lift coefficients, it reduced the flapping propulsive efficiency from $57\%$ at $\mathrm{k}=1.229$ to $38\%$ at $\mathrm{k}=1.8435$ .

Figure 18.

Reduced frequency effect: (a) coefficient values and (b) propulsive efficiency.

Static structure analysis

As it was previously demonstrated that both the corrugated airfoil and the flat-plate have a very similar aerodynamic performance, it was necessary to investigate their structural performance. A wing with $70\mathrm{mm}$ span length was built using the corrugated airfoil and the flat-plate cross sections with a chord length of $10\mathrm{mm}$ . The chord length was fixed along the wing span. Both wings were tested under a uniform static pressure distribution along their upper surfaces. The static test was conducted using the commercial software ABAQUS. The uniform applied pressure magnitude is $5\mathrm{Pa}$ . One of the wing’s side surfaces was fixed and the other one was free to deform under the applied pressure. Care must be taken as the ABAQUS software is dimensionless software. SI units were employed in the performed analysis. To perform the static test, Young’s modulus and Poisson’s ratio were only needed. The material properties are given in Table 3.

Table 3.

Material properties. ³⁸

Young’s modulus (MPa)	$6100$
Poisson’s ratio	$0.25$

In the present simulation, we are only interested in the wing’s deflection under the applied load to compare the stiffness between the corrugated airfoil and the flat-plate. Figure 19 represents the corrugated airfoil deflection. The corrugated airfoil thickness varies along the chord length with a maximum value of $0.18\mathrm{mm}$ thick. Figure 20 represents the flat-plate deflection in $\mathrm{m}$ with different cross-section thickness values.

Figure 19.

Corrugated airfoil deflection in m.

Figure 20.

Flat-plate deflection in m for different flat-plate thicknesses: (a) 0.18 mm, (b) 0.5 mm, (c) 0.65 mm, (d) 0.7 mm, and (e) 0.75 mm.

Table 4 represents the weight of the different tested cross sections and their corresponding deflection. We can clearly see that the flat-plate thickness which corresponds to the thickest region in the corrugated airfoil (i.e. $0.18\mathrm{mm}$ thickness) deflects almost $60$ times the deflection of the corrugated airfoil. As the flat-plate cross-section thickness increased, the deflection began gradually to decrease. However, the weight of the flat-plate cross section increased till the weight of the flat-plate cross section became almost $\mathrm{four}$ times larger than the corrugated airfoil weight to produce almost the same deflection.

Table 4.

Weight and deflection variation of the different tested cross sections.

Airfoil	Deflection (mm)	Weight (g)
Corrugated airfoil	$1.66$	$0.152$
Flat-plate $(0.18\mathrm{mm})$	$98.6$	$0.152$
Flat-plate $(0.5\mathrm{mm})$	$20.7$	$0.420$
Flat-plate $(0.65\mathrm{mm})$	$9.48$	$0.546$
Flat-plate $(0.7\mathrm{mm})$	$7.60$	$0.588$
Flat-plate $(0.75\mathrm{mm})$	$0.0918$	$0.630$

Finally, although the corrugated airfoil and the flat-plate have almost the same aerodynamic performance, the flat-plate cannot be used to replace the corrugated airfoil unless its weight increased.

Conclusion

In this study, A. cyanea dragonfly forewing mid-cross-section corrugated airfoil was simulated at ultra-low Reynolds number. The corrugated airfoil was compared with its smooth counterpart to study the effect of corrugations upon the aerodynamic performance. A flat-plate was also added to the comparison as a base airfoil. This study was divided into two main phases, namely, gliding phase and flapping phase.

For the gliding phase analysis:

• At $\mathrm{Re}=1400$ and $\mathrm{\alpha }=0°$ , we were able to capture small trapped recirculating flow inside the corrugated valleys. Theses trapped vortices gave the effect of dynamic smoothing of the flow streamlines around the corrugated airfoil to simulate the profiled airfoil. The trapped vortices generate a negative shear drag inside the valleys, which help in reducing the total corrugated drag to become comparable to the profiled airfoil.

• The corrugated airfoil was found to produce an increasing lift force among all the tested angles of attack for the two tested Reynolds numbers. The results indicated the importance of the corrugations in producing the lift force even if there is no recirculating flow inside the valleys at high angles of attack. At $\mathrm{Re}=200$ , the flat-plate produced almost the same lift force as the corrugated airfoil with approximately the same drag force. However, at $\mathrm{Re}=1400$ , the flat-plate aerodynamic performance was improved among the corrugated airfoil at $\mathrm{\alpha }=5°$ and $10°$ with higher drag force. At higher angles of attack, the corrugated airfoil performs better than the flat-plate. The drag force for all the tested airfoils was found to experience an increase in the values by decreasing the Reynolds number. At lower Reynolds number, the skin friction caused by the flow viscous effect becomes the main contributor to the drag force. The profiled airfoil was found to behave relatively bad at the two tested Reynolds numbers compared to the corrugated airfoil and the flat-plate; hence, it was excluded from the flapping analysis.

For the flapping phase analysis:

• The static pressure contours for the corrugated airfoil and the flat-plate revealed the importance of the LEV formation and the shedding of the TEV in generating both thrust and lift forces by forming the reversed von Kármán vortex street. The results obtained for both the corrugated airfoil and the flat-plate under the same circumstances showed slight differences between their aerodynamic performances. The lift force obtained for the corrugated airfoil was slightly higher than that for the flat-plate, which can be explained due to the camber presences for the corrugated airfoil. The thrust trend obtained was a twin peak for both tested airfoils. The generation of the thrust force during the downstroke was higher than that for the upstroke for the corrugated airfoil, while it has almost the same value for the flat-plate during the upstroke and the downstroke.

• Based on the performed simulations, increasing the reduced frequency has a great impact on the thrust and the lift forces. Increasing the reduced frequency from $0.768$ to $1.229$ enhanced the thrust and the lift forces values by $82\%$ and $75\%$ , respectively. Furthermore, by increasing the reduced frequency to $1.8435$ , the thrust and the lift forces further increased by $62\%$ and $95\%$ , respectively. Although we can find that increasing the reduced frequency has a significant impact on both the thrust and the lift forces, it reduced the flapping propulsive efficiency from $57\%$ at $\mathrm{k}=1.229$ to $38\%$ at $\mathrm{k}=1.8435$ .

• Although the corrugated airfoil and the flat-plate have the same aerodynamic performance, the flat-plate cannot be used to replace the corrugated airfoil due to its low stiffness properties.

Finally, based on the conducted simulations, the corrugated airfoil is a preferable choice in MAVs due to the high performance at low Reynolds number and high strength with light structure weight.