Parametric study of a corrugated airfoil in a forward flight at ultra-low Reynolds number

In the current study, the mid cross section of the dragonfly forewing was simulated at ultra-low Reynolds number. The study aims to understand better the contribution of corrugations found along the wing on the aerodynamic performance during a forward flight. Different flapping parameters were employed. FLUENT solver was used to solve unsteady, two-dimensional, laminar, incompressible Navier–Stokes equations. The results revealed that any stroke amplitude less than 1cm generated no thrust force. The stroke amplitude had to be increased to form the reversed Kármán vortices responsible for generating thrust force. The highest propulsive efficiency was found in the Strouhal number range 0.2 < St < 0.4 with a peak efficiency of 57% at St = 0.39. Changing the phase difference between pitching and plunging motions from advanced to synchronized caused lift force to drop 91% and thrust force to increase by 15%. On the other hand, changing the phase difference from synchronized to delayed caused lift and thrust forces to increase by 89% and 36%, respectively, and propulsive efficiency to deteriorate significantly. In all performed simulations, the airfoil was assumed to start motion from rest with no initial angle of attack. The increase in initial angle of attack generates a very high lift force with a fair loss for both thrust force and propulsive efficiency. The decomposition of flapping motion into its elementary motions revealed that the aerodynamic forces generated are a non-linear superposition from both pure pitching and pure plunging aerodynamic forces. This can be attributed to the non-linear interaction between unsteady vortices generated from these decomposed motions.


Introduction
Interest in the flapping wings field has increased rapidly over the past decades. Hence, a lot of efforts have been undertaken to enhance our understanding of flapping aerodynamics. The more we understand the better our micro aerial vehicles (MAVs) designs become. Due to the rapidly increasing interest in the flapping wings field over the past decades, a lot of effort has been put into broadening the understanding of said field. Consequently, this interest has brought upon an advantage: the deeper we delve into flapping aerodynamics, the better our micro aerial vehicles (MAVs) designs.
Compared to fixed wing theories, the flapping wing theory is quite complex. On one hand, fixed wings require the generation of lift force from the airflow around the wings, and the production of thrust force from the engines. On the other hand, birds and insects simultaneously produce these two forces just by flapping their wings [1]. In this case, flapping wings in nature can be classified into two main categories: bird-like and insect-like, while the latter can be subcategorized into a tandem wing (e.g. dragonflies and hawkmoths) and a single wing (e.g. bees and houseflies). Nature gave flapping wings the chance to be evolved over millions of years. Favored by nature-induced evolution, all birds and insects have well-defined and intertwined wings and abdomen. By flapping, birds and insects can produce both lift force and thrust force to make extreme maneuvers, fly effortlessly and sustain flight for a long duration [2]. However, neither fixed wings nor rotary wings exist naturally-they are purely man-made.
Both birds and insects use aerodynamic forces differently. Birds benefit from their high aspect ratio wings and high Reynolds number to generate lift force to sustain their weight [3]. The high aspect ratio is more suitable for steady state forward flight or soaring flight [4]. The higher the aspect ratio is, the better the birds' ability to maximize flight duration and range. Meanwhile, the aspect ratio of insects is very low compared to that of birds. Thus, they tend to use a high flapping frequency with varying amplitude to maximize the duration of flight and range [3]. The low aspect ratio helps insects improve their maneuvering [4].
Moreover, wing stroke kinematics was found to vary widely between natural flyers; it can even vary in the same flight phase depending on flight speed. For instance, as the flight speed decreases, the wing stroke plane becomes more horizontal, leading them to the hovering phase. However, the wing stroke plane becomes more vertical when they tend to fly faster. Stroke kinematics was also found to be more complicated as the flyer gets smaller in size; e.g. larger birds exhibit a relatively simple wingtip motion compared to the small flyers like the blowfly. Natural flyers can work their body, legs, tail in a perfect harmony enabling them to perform many sophisticated maneuvers such as takeoff, landing, and track a target precisely [1].
According to the classical theories of fixed wings, biological flyers' angles of attack can reach high values that can cause an aerodynamic stall. However, biological flyers were found to exhibit favorable aerodynamic forces at these high angles of attack [5] which was quite interesting to the researchers studying the flapping aerodynamic mechanisms.
Several experimental and numerical studies were made on dragonflies [6][7][8][9][10][11][12]. Many researchers were concerned with the flapping phenomena and the vortex structure generated rather than studying the wing airfoil shape. It was equally important to know how the vortices interact with each other and how these interactions affect the performance of the flapping creatures. Young et al. [13] studied the wake of a 2D NACA0012 airfoil oscillating in pure vertical sinusoidal motion at Reynolds number equals 2 × 10 4 . The analysis was made using a compressible Reynolds averaged Navier-Stokes flow solver. The flow was assumed to be both fully laminar and fully turbulent. They showed that for pure pitching airfoil and pure plunging airfoil, the Mach number has an important effect on the predicted forces, but it does not significantly influence the wake structure. The results indicate the importance of the leading-edge on the aerodynamic forces generated from the plunging motion at the tested Reynolds number: the trailing-edge was found to affect both the wake structure and the aerodynamic forces. The laminar wake structure results were found to reproduce the experimental results in the literature. Young et al. [13] showed that it cannot be a single parameter controlling the flow at the tested Reynolds number. However, the reduced frequency and its multiplication with the non-dimensional plunging amplitude have to be treated separately for a reduced frequency of less than 4.
Another numerical and experimental study was proposed on the wake structure of a 2D NACA0012 airfoil in hovering motion [14]. The airfoil's total displacement is fixed to 6c where cis the airfoil chord length. The angular displacement is fixed to range from − 45° to 45°. The numerical study was conducted to solve unsteady and laminar flow using a direct numerical simulation code. Based on the airfoil chord and the maximum motion velocity, the Reynolds number was 10 3 . The experimental analysis was performed in a water tunnel using the particle image velocimetry (PIV) technique. During the analysis over a one-time period, three types of vortices were observed and classified. These vortices are namely the leading-edge vortex (LEV) at the airfoil leading-edge, the translational vortex (TV) at the airfoil trailing-edge and the rotational stopping vortex (RSV), also at the airfoil trailingedge. During the airfoil translation motion, both LEV and TV were generated. However, the RSV was only generated during the airfoil rotation motion. These generated vortices were found to interact with the generated vortices from the previous stroke, affecting the aerodynamic forces. The experimental results gave a more complicated vortex topology compared to those of the numerical simulation.
The response surface methodology (RSM) optimization technique was used to maximize the propulsive efficiency of a 2D plunging NACA0012 airfoil [15]. The unsteady laminar flow was solved using FLUENT for various Reynolds numbers. A user-defined function (UDF) was utilized to account for the airfoil motion. The reduced frequency and the plunging amplitude were considered the control parameters. The results showed that both the thrust coefficient and the input power efficiency increase with the increase of the reduced frequency. However, for high values of reduced frequency, the input power coefficient increases faster than the thrust coefficient. Therefore, a lower propulsive efficiency was obtained. Mekadem et al. [15] found that at any Re, high propulsive efficiency can be achieved by applying a high frequency which corresponds to a large thrust force and reduces the input power by adopting low plunging amplitude.
Additionally, a 2D NACA0012 airfoil was used to model the forward flight of a flapping insect wing at Re = 150. The plunging and the pitching motions' amplitudes were fixed at 0.4c and 45°, respectively. Both the mean angle of attack and the phase difference were kept at 0°.  [16] reported that the ground effect either helps the flapping airfoil to reduce the mean drag force at high frequencies or enhances the mean lift force at low frequencies. They reported the strong influence of the ground on the vortex shedding process. This influence became stronger by increasing the flapping frequency, which in-turn can form the oblate vortices in the wake.
Furthermore, Lua et al. [17] conducted a 2D experimental and numerical investigation to study the vortex evolving with Strouhal number (St) upon thrust coefficient. The experiment was made using a re-circulating water channel using the digital particle image velocimetry (DPIV) technique. The numerical analysis was conducted using FLUENT at Re = 5000. The flow was treated as a laminar flow. To simulate the insect wing, a 2D elliptic airfoil was used within 0.1 < St < 0.8 with different induced effective angles of attack. The elliptic airfoil has an aspect ratio (span/chord) of 10 and 5 mm as a maximum thickness. The obtained results showed that regardless of the effective angle of attack's values, the thrust coefficient ( C t ) increased with the increase in St until it reached a peak value at a critical Strouhal number value ( St cr. ) then deteriorates beyond this value. A higher value of the effective angle of attack would cause a higher St cr. , and as a result, a higher C t peak. There are three classified types of C t over one flapping stroke: type-1, characterized by a single peak, type-2, similar to type-1 except the flat edge after the stroke reversal and type-3, which contains multiple peaks. Below the St cr. , type-1 and type-2 were found while type-3 occurred mainly beyond the St cr. . For type-3, the airfoil produced drag force most of the flapping cycle time. The drag force production increased with the increase in the St. The deterioration of type-3 can be attributed to the adverse domination effect of the airfoil rotation rate on the near wake region. However, the favorable effect of the near structure wake was dominant for low St. The three types' generations were influenced by the combined effects of the adjacent wake structure and the rate of wing rotation. Evidently, regardless of the effective angle of attack value, the maximum propulsive efficiency occurs within 0.25 < St < 0.35.
The present study aims to gain the necessary knowledge needed to design MAVs by studying the performance of Aeshna Cyanea's dragonfly flapping insect. The mid cross section corrugated airfoil of its front wing will be considered the main airfoil. The numerical approach is listed first dealing with the conducted independent studies and the adopted flapping equations to simulate the forward flight. Different parameters will be altered during the forward flapping motion to study their effect on the airstream flow, hence on the aerodynamic performance. Although the performed 2D analysis will not be able to illustrate the actual performance of the corrugated airfoil due to the spanwise flow, it will give a closer look at the behavior of the flow along the corrugated surfaces. In Sect. 2, the numerical method is represented beside the adopted flapping motion equations in the current study supported by all the necessary independent studies. Section 3 illustrate the results for the varied parameters during the flapping motion. This section is divided into four subsections, each subsection contains the tested parameter with the obtained results and any necessary comparison. Finally, the conclusion is depicted in Sect. 4. Figure 1a illustrates the corrugated airfoil corresponding to the forewing middle cross section of the Aeshna Cyanea dragonfly while Fig. 1b illustrates the simplified computational shape used in the present numerical study. The airflow over the corrugated airfoil was solved using the FLUENT solver.

Numerical method
The validation problem was conducted for a flat-plate in a pure pitching motion with pre-specified motion parameters. The validation problem was solved and compared to the published numerical and experimental solution of the problem [19][20][21]. The present numerical results show a good agreement with the literature [18].
To obtain the numerical parameters used in the present study, several independent studies were performed including the domain independent study, the time independent study and the mesh independent study. The interested readers can refer to Ref. [18] for more details about the used numerical parameters, the independent studies' results and the validation problem.

Flapping phase
In order to understand the aerodynamics behind the flapping mechanisms, the Strouhal number has to be considered. The Strouhal number (St) is a non-dimensional parameter used to describe the wake patterns behind a stationary object such as Kármán vortices or to describe the unsteady flow behind a 2D airfoil undergoing flapping motion [1]. The Strouhal number can be defined as: where A o denotes the wake width approximated by the peak-to-peak maximum excursion of the airfoil's trailing-edge and can be approximated also to be twice the plunging amplitude. f denotes the flapping frequency while U ∞ denotes the free stream velocity. The Strouhal number can be related to the reduced frequency as [15]: where y o denotes the plunging amplitude, k denotes the reduced frequency and c m denotes the mean airfoil chord length. For a certain number of Strouhal numbers, the airfoil produces a forward thrust force rather than drag force. The thrust force is produced due to the formation of vortices that are similar to the Kármán vortices' structure but with a reversed vortices direction; therefore, thrust force vortices are named reverse Kármán vortices [1]. The Strouhal number for the swimmers is in the range of 0.2 < St < 0.4 where the propulsive efficiency is quite high. The natural flyers are in the same range, but only during the cruising phase [22]. The Strouhal number is not applicable in the hovering phase as the forward flight speed approaches zero.
For all the tested cases, the airfoil center of rotation was made at quarter the airfoil mean chord length (c m /4). The numerical investigation was performed at a fixed Reynolds number of 1400 and a mean chord length of 10 mm.
Equations (4) and (5) were used to describe the coupled corrugated airfoil plunging and pitching motions.
where o is the pitching amplitude and is the phase difference between the plunging motion and the pitching motion. Varying the phase difference yields advanced motion, synchronized motion, and delayed motion. For the power and the thrust efficiency calculations' details, interested readers can refer to Ref. [18] for the equations.  The lift and the thrust coefficients for the corrugated airfoil at different plunging amplitudes during one complete flapping cycle are shown in Fig. 2. Figure 2a illustrates that by increasing the plunging amplitude, the lift coefficients' peaks increase. Although the flapping motion is considered symmetrical, the aerodynamic forces generated during the upstroke and the downstroke are different. These differences can be attributed to the corrugation pattern along the upper and lower surfaces of the airfoil as well as the non-linear interaction between the generated vortices. Fortunately, these differences resulted in producing average positive lift forces for all the tested plunging amplitudes.

Strouhal number effect
As shown in Fig. 2b, the thrust coefficient is a twin peak curve as the thrust force is produced in both upstroke and downstroke. The second peak is slightly lower than the first peak, which is evident in y o = 1c m . As the plunging amplitude decreased, the peak value for the thrust force was more inclined to the right. It is clear that based on the strength of the generated vortices, the amplitude of the coefficient forces varies. It was noted that for y o = 0.25c m the airfoil was not able to produce thrust force during any phase of the flapping cycle. However, for y o = 0.5c m , the airfoil was only able to produce small amounts of thrust force at some phases during the flapping cycle. For y o = 1c m or even higher, the airfoil was able to produce thrust force during all the flapping cycle phases.
Static pressure contours were needed to gain better insight into the production of thrust and lift forces. The static pressure contours will enhance our understanding by visualizing the vortices formation and shedding process. Figure 3 represents the instantaneous static pressure contours of the corrugated airfoil for various plunging amplitudes (i.e.,y o = 0.25c m , 0.5c m , 1c m , 1.5c m ) during one complete flapping cycle. Note that the color legend has been truncated to be in the range from −22 to 5 to easily spot the differences between the plunging amplitudes' contours. Figure 3 shows that as the plunging amplitude increases, as the strength of the developed vortices along the airfoil surface also increases. In Fig. 3a, we can recognize two different sized residual vortices for the case of y o = 1.5c m . One of the vortices is positioned at the aft portion of the airfoil and is fully detached from its surface with a dual vortex core, while the other vortex is positioned near the airfoil leading-edge with a relatively weak strength. These vortices' strength and position contributed in producing a net positive lift force at the beginning of the flapping cycle in compliance to the other plunging amplitudes. It can be pointed out that by increasing the plunging amplitude, the residual vortices became more detached from the airfoil surface. Thus, the pressure along the lower airfoil surface was higher than that on the upper surface. It was found that at the beginning of the flapping cycle all the tested plunging amplitudes produced no thrust force. In Fig. 3b, as the airfoil accelerates in the downstroke phase, the leading-edge vortex (LEV) was developed for all the tested plunging amplitudes. As the LEV began to separate from the airfoil surface, the aerodynamic forces attained their maximum values. The LEV strength depends heavily on the plunging amplitude; however, the structure of the vortex was not affected [23]. By increasing the plunging amplitude, the strength of the LEV increased which in-turn enhanced the lift and the thrust forces for the airfoil. As the airfoil continued the downstroke, the LEV moved further downstream the upper surface, which caused the airfoil to stall and the aerodynamic forces to deteriorate. As the airfoil began to While no airfoil at the wing pronation had the ability to produce thrust force, they were able to produce lift force. The motion is symmetrical so, what applied for the downstroke motion and wing pronation, also applied for the upstroke and wing supination. The airfoil corrugations were the main reason for aerodynamic forces differences' between the upstroke and downstroke motions. The stronger vortices generated by increasing the plunging amplitude can be attributed to the increase in the effective angle of attack during the flapping cycle. However, it is believed that the bigger the size of the generated LEV, the less effective the corrugations will be on the aerodynamic forces.
In Fig. 4, one can see that increasing the stroke amplitude above a certain limit results in producing a thrust force instead of a drag force. The results indicate that by increasing the stroke amplitude up to 1c m , the wake shedding forms high strength reversed Kármán vortices instead of regular Kármán vortices, thus producing drag force. The lift coefficient was found to be enhanced by 55% by increasing the stroke amplitude from 0.25c m to By further increasing the stroke amplitude to 1.5c m , both thrust and lift forces were enhanced. The lift force was found to be almost 23 times the lift force produced with the stroke amplitude of 1c m . Moreover, the thrust force was found to be three times larger than the thrust force produced by the stroke amplitude of 1c m . On the other hand, a loss in the propulsive efficiency was noticed for increasing the stroke amplitude to 1.5c m by more than 1.5 times the value for the 1c m stroke amplitude. The highest propulsive efficiency was found in the Strouhal number range 0.2 < St < 0.4 , with a peak efficiency of approximately 57% at St = 0.39 which is in a good agreement with the literature.

Phase difference effect
The wing rotation's timing and duration were found to have a substantial effect on the resulting aerodynamic coefficients as well as on the power producing efficiency of most flying insects. At the end of each half cycle, the wing rotates rapidly (i.e. supination and pronation). Different wing rotation can be achieved by changing the phase difference angle between the plunging motion and the pitching motion. Three types of rotation were investigated: advanced, symmetrical, and delayed rotations, as illustrated in Fig. 5. In advanced rotation, the wing rotates before the end of the stroke, while in delayed rotation, the wing rotates after the beginning of the stroke. In symmetrical rotation, the wing rotates at the end of the stroke.
In the present numerical simulation, three phase difference values were chosen as 60°, 90° and 120° to yield advanced rotation, synchronized rotation and delayed rotation. The Reynolds number, frequency, stroke amplitude and pitching amplitude were kept fixed at 1400, 40Hz, 1c m and 30°, respectively. Figure 6 represents the lift and the thrust coefficients for the corrugated airfoil at three different rotations during one complete flapping cycle. It is obviously illustrated, for both the lift and thrust forces, that delayed rotation was the first to reach its peak values, followed by the synchronized rotation and advanced rotation, respectively. Figure 6a shows that the lift force amplitude for advanced rotation is the highest among the other two rotations. Nevertheless, the force amplitude for synchronized rotation and delayed rotation is almost the same, with the maximum value shifted from each other. Figure 6b indicates the twin peak behavior of the thrust coefficient as the thrust force is being produced in both the upstroke and the downstroke, with the second peak lower than the first peak. Advanced rotation was found to produce the highest amount of thrust force, while delayed rotation produced the lowest amount. We can explain the higher performance of advanced rotation by the increase in the circulation around the airfoil. This circulation will increase the total airflow velocity on one surface of the airfoil and decrease it on the other surface. Furthermore, airfoil advanced rotation increases the strength of the developed vortex on the airfoil surface hence, the strength of the reversed Kármán vortices at the end of each half stroke. Figure 7 represents the static pressure contours for the corrugated airfoil for various phases' difference. In Fig. 7a, during the beginning of the downstroke phase, there are two well-shaped, differently-sized residual vortices on the lower surface of the corrugated airfoil during the advanced rotation. The smaller vortex can be seen clearly near the corrugated airfoil leading-edge, while the larger one was found to occupy the mid-lower surface of the airfoil. These residual vortices were stronger than those developed for synchronized rotation, thus resulting in decreasing the lift force by almost more than three times the value of that of synchronized rotation. Contrarily, for delayed rotation, there is a shed residual trailing-edge vortex at the beginning of the downstroke. The pressure on the lower airfoil surface was higher than that of the upper surface, resulting in the production of a positive lift force.
Unlike advanced and delayed rotation, synchronized rotation produces drag force instead of thrust force at the beginning of the cycle. As the airfoil continued the downstroke movement and accelerates as the residual vortices began to completely detach and shed from the airfoil surface. Figure 7b illustrates that for synchronized rotation, the residual vortex had already been shed from the airfoil lower surface, and a LEV began to form at the upper surface of the airfoil.
However, for advanced rotation, the residual vortex was almost shed from the trailing-edge to form the reversed Kármán vortices. For delayed rotation, the LEV was formed on the airfoil's upper surface and moved downstream to occupy the aft portion of the airfoil's surface. The aerodynamic forces were enhanced as long as the formed LEV was attached to the airfoil's surface near the leading-edge, and began to deteriorate due to the movement of the LEV downstream. Figure 7c reveals the airfoils deceleration for wing pronation. Clearly, one can see the beginning of the LEV formation for advanced rotation, while for synchronized rotation, the LEV has already been formed and moved downstream from the airfoil leading-edge. For delayed rotation, the vortex on the upper surface of the airfoil can be seen in the form of weak dual core above each other at the trailing-edge. These vortices produced negative lift force and began to produce drag force instead of thrust force. Due to the motion symmetry, same explanation can be used for the upstroke motion and the supination.
We can clearly see that the same sequence of vortices forming and shedding from either airfoil surfaces was the same for the three chosen rotations. Synchronized rotation can be seen to be ahead of advanced rotation by approximately one step, which in turn, lagged the delayed rotation by approximately one step. This caused the coefficients peak values to be shifted from each other. Figure 8 represents the thrust coefficient, the lift coefficient and the propulsive efficiency variation against the phase difference. Advanced rotation generated the highest lift force compared to the other two motions in accordance with Ref. [24]. Changing the phase difference from advanced rotation to synchronized rotation causes the lift force to drop about 91% of its value. However, this reduction in-turn causes the thrust force to increase by approximately 15%, with an increase in the propulsive efficiency from 53% to 57%. On the other hand, when the phase difference changed from synchronized rotation to delayed rotation, lift force increased once again by about 89% but it caused a 36% loss in thrust force. Moreover, the propulsive efficiency was found to deteriorate significantly when the phase difference became delayed rotation.

Initial angle of attack effect
In all conducted simulations, we assumed that airfoil will start motion from rest with zero initial angle of attack (i.e. the airfoil is completely horizontal). However, in nature, this assumption is not completely true. Insects start their wing motion with very high initial angles of attack which can approach 60°, going higher or lower depending on the required maneuver and the wing's morphological limit. For the current performed simulation, the high initial angles  of attack will cause a high mesh deformation based on the adopted mesh strategy which consequently requires a high computational effort.
In order to study the effect of the initial angle of attack on the generated aerodynamic forces, a 15° initial angle of attack was chosen [9]. Figure 9 illustrates the lift and the thrust coefficients for the corrugated airfoil at both 0° and 15° initial angles of attack. Figure 9a clearly shows the superior performance of the corrugated airfoil with the initial angle of attack compared to the one with zero initial angle of attack. The initial angle of attack contributes to generating a higher lift force through the downstroke phase compared to the zero-initial angle of attack case. Through the upstroke phase, the zero-initial angle of attack generates a negative lift force (i.e. the second peak during the flapping cycle) which can reverse the effect of the first peak, thus lowering the average lift force. In compliance, the corrugated airfoil with the initial angle of attack produces a lift force value which almost fluctuates around zero, hence increasing the total average lift force during the flapping cycle. Compared to the case with zero initial angle of attack, the initial angle of attack helps in producing a positive lift force during the start of the flapping cycle.
For the case with a 15° initial angle of attack, Fig. 9b shows that the thrust coefficient is produced once during the flapping cycle. The thrust force was found to be generated only during the downstroke, while the drag force was produced during the upstroke in compliance with the case with a zero initial angle of attack where the thrust is being produced during both phases. The reason lies in the fact that during the upstroke phase, increasing the initial angle of attack change the position of the airfoil to a more vertical one facing the airflow, which in-turn produced drag force instead of thrust force. It was noted that the thrust force produced during the downstroke for the case with an initial angle of attack is less than that produced for the case with zero initial angle of attack and more shifted to the left. Table 1 illustrates the coefficient values and the propulsive efficiency produced during the flapping cycle for the two different initial angles of attack. The thrust force generated for the case with zero initial of attack was found to be three times more than that produced with an initial of attack with a higher propulsive efficiency. However, the case with initial angles of attack was found to compensate the loss in thrust force by generating very high lift force (i.e. 126 times the lift force) compared to that with zero initial angle of attack. The propulsive efficiency for the case with zero initial angle of attack was 1.5 times higher than that of the case with the initial angle of attack.

Flapping motions effect
The flapping wing motion can be classified into three categories based on the wing motion kinematics and the forces generation mechanism: the vertical stroke plane, the inclined stroke plane and the horizontal stroke plane.  In the present simulation, we are interested only in the vertical stroke plane. In general, flapping motion can be decomposed to its elementary motions: pitching motion and plunging motion as shown in Fig. 10. Although the flapping motion (i.e. combined motion) is a linear superposition of pitching motion and plunging motion, aerodynamic forces generated from the combined motion are not a linear superposition of the aerodynamic forces generated from the pitching motion and the plunging motion [25]. To understand the contribution of the basic elementary motions' kinematics to the unsteady aerodynamic forces generated for the combined motion, a pure pitching motion and a pure plunging motion were investigated separately. Figure 11 represents the lift and the thrust coefficients generated over the corrugated airfoil during one complete flapping cycle for different flapping motions. In Fig. 11a, one can see that the lift coefficient amplitude for the pure plunging motion is the highest among the other lift coefficients generated from pure pitching motion and combined motion. The pure pitching airfoil was found to generate the least amount of lifting force during the cycle which gave an average negative lift force. The pure pitching airfoil was found to reach its peak value first, while the combined motion reached its peak value last.
Thrust force can be categorized by twin peak in all three different motions. Figure 11b indicates that pure pitching motion is the only motion that was not able to produce any thrust force during the flapping cycle. The Strouhal number for pure pitching was calculated based on the relation given by Eq. 6 to be 0.13 [25].
The value of the Strouhal number is well below the critical value for producing a thrust force, as the shedding vortices formed are the regular Kármán vortices. The pure plunging motion produced thrust force during the downstroke motion, while the ability to produce thrust force decreased during the upstroke. Hence, the average thrust force produced was low. The combined motion was found to outperform the other two implemented motions in generating thrust force.
From Fig. 11a, the pure pitching motion was the first to reach its peak value, followed by the pure plunging motion and the combined motion, respectively. However, in Fig. 11b, pure plunging motion was the first to reach its peak value followed by pure pitching motion with respect to the negative sign and the combined motion, respectively. These aerodynamic forces' behavior can be attributed to the formation and shedding of the LEV. Figure 12 illustrates the static pressure contours for the three different motions implemented in the current analysis. Obviously, we can see that the vortices developed over the airfoil during the plunging motion are much stronger and bigger than those developed for the two other motions. For pure plunging motion, a high strength residual vortex near the trailing-edge in Fig. 12a  Meanwhile, for the pure pitching motion, a weak strength residual vortex is attached to the airfoil lower surface near the trailing-edge. The pressure difference between the upper and the lower surfaces results in producing a negative lift force at the beginning of the flapping cycle. No airfoil was able to produce thrust force at the cycle's beginning. As the periodic time passed, the airfoil began to accelerate in the three different implemented motions. Figure 12b illustrates the formation of a high strength LEV over the upper surface of the airfoil for the pure plunging motion. This LEV causes the lift force to attain its almost peak value during the flapping cycle. However, at the same time, the thrust force had already reached its peak value and began to deteriorate. For the pure pitching motion, the LEV results in increasing the lift force to reach its peak value during the flapping cycle. Although the pure pitching motion was able to produce a positive lift force during the downstroke, it was not able to produce a positive thrust force. Figure 12c represents the airfoil deceleration which, in return, causes the LEV to move further downstream the airfoil's upper surface. For the pure plunging motion, the LEV was found to be enormous enough in size to occupy almost the aft portion of the airfoil's upper surface. The LEV is now completely detached from the airfoil upper surface, causing the lift and the thrust forces to deteriorate. For the pure pitching motion, the LEV is still relatively attached to the airfoil upper surface. The same conclusion can be observed during the upstroke motion. Figure 13 illustrates the effective angle of attack variation during the flapping cycle for the three different motions. For pure plunging motion, we can see that the effective angle of attack oscillates between 51° and − 51°, which results in producing a high resultant unsteady aerodynamic force. Although the unsteady aerodynamic resultant force was able to produce high lift force amplitude, it produced low thrust force amplitude. For pure pitching motion, the effective angle of attack was found to oscillate between 30° and − 30°. Although the effective angle of attack was relatively high, the pure pitching motion could not produce a positive lift force. Furthermore, the pure pitching motion produced drag force instead of thrust force. On the other hand, the combined motion's effective angle of attack was found to oscillate between 21° and − 21°. It can be noticed that the effective angle of attack for the combined differed slightly in pattern from those obtained for the other two motions. The combined motion was found to outperform the other two motions in producing thrust force regardless the low effective angle of attack. Despite the fact that the combined motion could not produce as much lift force as the pure plunging motion could, it produced an average positive lift force during the complete flapping cycle.
It is obviously depicted in Figs. 11a and 11b that the aerodynamic forces generated for the combined motion are a non-linear superposition from the pure pitching aerodynamic forces and pure plunging aerodynamic forces. This can be attributed to the non-linear interaction between the unsteady vortices generated from the pure pitching motion and the pure plunging motion, as can be seen in the pressure contours illustrated in Fig. 12.

Conclusions
In the current study, the corrugated airfoil was employed numerically to investigate the effect of corrugations on the aerodynamic performance at ultra-low Reynolds number by altering different parameters during the forward flight. The corrugated airfoil used corresponds to the mid cross section of Aeshna Cyanea dragonfly forward wing [10]. The analysis of Strouhal number indicates that the stroke amplitude has to be above a certain limit to produce thrust force. The lift force was enhanced by 55% by increasing the stroke amplitude from 0.25 to 0.5 c m , but it then declined by 76% as the stroke amplitude increased to 1 c m . By further increasing the stroke amplitude to 1.5 c m , lift and thrust forces were almost 23 times and 3 times larger than those produced for 1 c m stroke amplitude. On the other hand, a loss by more than 1.5 in the propulsive efficiency was noticed. The highest propulsive efficiency was in the range 0.2 < St < 0.4 with a peak efficiency of 57% at St = 0.39. In the same manner, changing the phase difference has a great impact on the flapping performance. Utilizing synchronized rotation instead of advanced rotation caused lift force to drop 91% of its value. However, this change in-turn caused thrust force to increase by 15%, with 3% increase in the propulsive efficiency. On the other hand, as the phase difference changed from synchronized to delayed rotation, lift force increased once again by 89%, but it caused a 36% loss in thrust force. Moreover, the propulsive efficiency deteriorated significantly. Another key point that has a greater impact on the performance is the initial angle of attack. Adjusting the initial angle of attack to be 15° at the beginning of the flapping cycle caused lift and thrust forces to be 126 times higher and 3 times lower compared to those with no initial angle of attack, respectively. Furthermore, the propulsive efficiency was 1.5 times higher for the case with no initial angle of attack. Finally, to be able to understand the flapping motion better, the motion was decomposed to its elementary motions, namely, the pitching motion and the plunging motion. The pure plunging motion produced the highest lift force while the combined motion produced the highest thrust force. In contrast, the pure pitching motion was incapable of producing neither positive lift force nor thrust force. The aerodynamic forces generated for combined motion are a non-linear superposition from pure pitching and pure plunging aerodynamic forces. This can be attributed to the non-linear interaction between the unsteady vortices generated from the two motions.
Optimizing the corrugated airfoil will be our priority in future studies to be able to compensate the small details that we cannot model or yet fully understood to reach the superior performance for future MAVs. Optimizing the airfoil shape is relatively easier than optimizing the flapping parameters due to the imposed mechanical limitations. Moreover, the effect of wing flexibility is very important to be considered to help us understand the passive control used during incest's flight.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.