Aerodynamic Shape Optimization of a Wavy Airfoil for Ultra-Low Reynolds Number Regime in Gliding Flight

The effect of the number of waves and the width of the ridge and valley in chord direction for a wavy airfoil was investigated at the angle of attack of 0 ∘ and Reynolds number of 10 3 through using the two-dimensional direct numerical simulation for four kinds of wavy airfoil shapes. A new method for parameterizing a wavy airfoil was proposed. In comparison with the original corrugated airfoil profile, the wavy airfoils that have more distinct waves show a lower aerodynamic efficiency and the wavy airfoils that have less distinct waves show higher aerodynamic performance. For the breakdown of the lift and drag concerning the pressure stress and friction stress contributions, the pressure stress component is significantly dominant for all wavy airfoil shapes concerning the lift. Concerning the drag, the pressure stress component is about 75 % for the wavy airfoils that have more distinct waves, while the frictional stress component is about 70 % for the wavy airfoils that have less distinct waves. From the distribution of pressure isoline and streamlines around wavy airfoils, it is confirmed that the pressure contributions of the drag are dominant due to high pressure on the upstream side and low pressure on the downside; the frictional contribution of the drag is dominant due to large surface areas of the airfoil facing the external flow. The effect of the angle of attack on the aerodynamic efficiency for various wavy airfoil geometries was studied as well. Aerodynamic shape optimization based on the continuous adjoint approach was applied to obtain as much as possible the highest global aerodynamic efficiency wavy airfoil shape. The optimal airfoil shape corresponds to an increase of 60 % and 62 % over the aerodynamic efficiency and the lift from the initial geometry, respectively, when optimal airfoil has an approximate drag coefficient compared to the initial geometry. Concerning an fixed angle of attack, the optimal airfoil is statically unstable in the range of the angle of attack from − 1 ∘ to 6 ∘ , statically quasi-stable from − 6 ∘ to − 2 ∘ , where the vortex is shedding at the optimal airfoil leading edge. Concerning an angle of attack passively varied due to the fluid force, the optimal airfoil keeps the initial angle of attack value with an initial disturbance, then quickly increases the angle of attack and diverges in the positive direction.


Introduction
Recently, micro-air vehicles (MAVs) [1], nano-air vehicles (NAVs) [2], and pico-air vehicles (PAVs) [3] have been widely used by researchers, security and law enforcement agencies, search and rescue operators, firefighters, farmers, filmmakers, photographers, and delivery companies. Motivated by the enormous global appetite for applying MAVs, NAVs, and PAVs, the design of high aerodynamic performance airfoils for these applications has received considerable attention in recent years. Micro-, nano-, and pico-air vehicles are working at low Reynolds number (Re) [4,5]. In particular, it is worth mentioning Defense Advanced Research Projects Agency (DARPA) specifications for NAVs with an extremely small wingspan <7.5 cm and working at Re < 15,000, and PAVs working at Re < 3000 [6]. Thus, the study of aerodynamics and shape optimization for airfoils at ultra-low Reynolds numbers becomes pertinent for engineered flying objects. For the Reynolds number region from 20,000 to 5 × 10 5 , various investigations of aerodynamic phenomena have appeared [5,[7][8][9]. However, the databases of airfoils at Reynolds number less than this range airfoil data below this range are very limited [10][11][12][13].
Of particular interest in this study is the level of Reynold number at 10 3 . Here, this Reynolds number regime is called ultra-low. Although many MAVs, NAVs, and PAVs have been designed over the last few years, the typically applied airfoil platform is still the contractible version of that applied for the giant aircraft at high Reynolds number. Due to the difference between the high Reynold number regime and the ultra-low Reynolds number regime, such as the dominant inertial force for the high Reynold number regime but dominant viscous force for the ultra-low Reynold number regime, the optimal airfoil design for the giant aircraft might be distinguished momentously from that for MAVs, NAVs, and PAVs. Therefore, there is an important requirement to readjust the typical streamlines airfoil platform for ultra-low Reynolds numbers regime.
One way of designing effective airfoils at ultra-low Reynolds number is by biomimetics [5,[14][15][16]. The dragonfly is considered as a high-performance flyer. The typical range of the Reynolds number for dragonflies is from 100 to 10,000 [17], which can be categorized as the ultra-low Reynolds number flow regime. Dragonflies have both gliding and flapping flight modes. When flapping, they can suddenly accelerate and stop, hover, and turn [18]. In gliding flight, the dragonfly elevates into the air using powered (flapping) flight and makes use of potential energy to move horizontally above the ground [19,20]. It is generally known that the dragonfly has highly corrugated wings where the wing's cross-section varies along the chordwise direction, as seen in Figure 1, unlike that of a typical engineered airfoil which is streamlined and smooth cambered. Extensive biological and aerodynamic studies of corrugated airfoils have been carried out by many researchers [18,[21][22][23][24][25][26]. Visualization of the two-dimensional flow field around the blade cross section using a bent airfoil has been made by Komine [27]. Our previous study has conducted numerical analysis of the twoand three-dimensional flow fields around the corrugated airfoil and quantitative evaluation of the aerodynamic characteristics [28]. It has been found that the corrugated airfoil performs significantly better than the streamlined airfoil at ultra-low Reynolds number with the low angle of attacks, and performs as well at the high angle of attacks. There are several attempts to explain the reason of the unexpected improvement in the aerodynamic performance for the corrugated airfoil at ultra-low Reynolds number. One explanation is that air flow over the airfoil corrugations is trapped between the valleys and essentially forms a virtual streamlined airfoil [20,23]. Another explanation is that the delay of the flow separation or the earlier reattachment of the flow separation for the corrugated airfoil corresponds to the improved performance [18,25]. In spite of different explanations for the improved aerodynamic performance, these studies unanimously agree that the corrugated airfoil performs well in ultra-low Reynolds number regimes. Thus, such corrugated airfoils could be candidates for MAVs, NAVs, and PAVs. Since the dragonfly's flight forms include flapping and glide, the corrugated airfoil shape can be considered to have been optimized for multiple purposes. Research on flapping using corrugated wing has been conducted by experiments and numerical calculations from the viewpoint of developing MAVs, NAVs, and PAVs [29][30][31]. Flapping wings are said to be able to obtain higher thrust than gliding wings, but application issues such as complicated mechanisms results in increased weight and less flight time. For fixed wing MAV, NAVs, and PAVs, wings that simultaneously provide a superior aerodynamic performance and structural robustness are critical. However, there is no guarantee that wing sections like dragonflies are optimized for aerodynamic characteristics during gliding. Therefore, in order to obtain an airfoil shape that exhibits higher aerodynamic performance in gliding at ultra-low Reynolds number, aerodynamic shape optimization for a corrugated airfoil of dragonflies is necessary.
Over the past several decades, aerodynamic shape optimization has been successfully applied for two-dimensional and simple three-dimensional configurations. General aerodynamic shape optimization includes a global method based on a heuristic algorithm and gradient method. The global optimization method is a probabilistic method that simulates biological phenomena and physical phenomena, as represented by genetic algorithms [32], artificial neural networks [33], and response surface method [34]. This method can search for solutions without convergence to local solutions even when the objective function has multimodality. However, it takes up huge amounts of computational resources and is time-consuming to obtain an optimal solution. On the other hand, researchers [35] found efficient global optimization of expensive black-box functions by using response surfaces for global optimization lies in balancing the need to exploit the approximating surface with the need to improve the approximation. For the gradient-based aerodynamic shape optimization, through simplifying the description of real-life problems, the optimal solution is searched by analytically finding the gradient vector and Hessian matrix related to the design variables concerning the objective function. In particular, adjoint-based methods [36][37][38][39][40] seem to be an attractive alternative, since the sensitivity analysis is independent of the number of design variables and proportional to the number of aerodynamic cost functions. In addition, the cost function gradient concerning the design variables is calculated indirectly by solving the flow governing equations and the adjoint equation, in which the cost of calculating the adjoint equation is almost the same as that of solving the flow equations. Based on the order of discretization and derivation of the adjoint equations, the adjoint-based method can be further classified as continuous and discrete. The advantage of the continuous adjoint method is that it is independent of the method used to solve the flow equations. Adjoint-based methods have been utilized in diverse areas such as aerospace [41,42], marine [43], and biomedical engineering [44]. In addition, they have been used to design optimal airfoil [6,39,43,[45][46][47]. In this study, an aerodynamic shape optimization procedure for a corrugated airfoil of dragonflies using a continuous adjoint method is implemented. A time marching finite difference method is utilized to solve the flow and adjoint equations. The properties of the adjoint equations are very similar to the flow equations, thus many numerical schemes for flow solvers can be applied to the adjoint solvers as well. The simple steepest descent algorithm is employed to minimize the objective function.
By the way, many of these studies treat wings as rigid bodies. Naka and Hashimoto have conducted experiments using a wing model of an elastic film and a rigid rod that reproduces the deformation of a dragonfly wing when flapping and gliding flight [48,49]. Numerical calculations using this wing model are also performed for glide, and the flow is left by deforming the wing so that it twists around the span direction axis. It is clear that the drag coefficient reduction rate is increased more than that of the lift coefficient, resulting in a higher lift-drag ratio than the rigid wing. Numerical calculation of the corrugated airfoil considering deformation is also important for evaluating the strength of wing for receiving force from fluid. In particular, when handling a rigid three-dimensional corrugated wing, it is predicted that the moment of force around the root of the wing will be larger than that of an elastic wing. At present, few studies have been conducted to evaluate and compare the hydrodynamic force and moment of force acting in each case where the corrugated wing during flight is an elastic body and a rigid body. Therefore, in this study, the hydrodynamic moment of force and the attitude stability of the optimal airfoil for ultra-low Reynolds number regime in gliding flight is investigated as well.
One objective of this study is to conduct a systematic investigation to study the improvement in aerodynamic performance for a series of wavy airfoil profiles generated by modifying the number of waves and the width of the ridge and valley in the chord direction. A novel scheme is proposed to represent corrugations, that is, the width of the ridge and valley of a wavy airfoil in the chord direction is expressed by the wavenumber of the sine wave. The height of the wave in the chord direction is limited by the envelope of the aircraft airfoil NACA20408. First, the effect of the number of waves and the width of the ridge and valley on the lift coefficient, drag coefficient, and the aerodynamic efficiency are investigated through intentionally changing waves for the airfoil. Various cases are considered. The number, width, and height of waves are varied on the whole airfoil, while the thickness of the wavy airfoil is fixed at 1% of the chord length. The effect of the angle of attack on the aerodynamic efficiency for various wavy airfoil geometries is studied as well.
Another objective of this study is to optimize the waves for high aerodynamic performance (in terms of lift-to-drag ratio) using a continuous adjoint method. The aerodynamic shape optimization is conducted on a high aerodynamic efficiency wavy airfoil, in which the NACA2408 airfoil is used as the envelope for the height of corrugations. Performance of the optimized wavy airfoil, the original corrugated airfoil, and NACA2408 airfoil are compared. Due to our application pointing to ultra-low Reynolds number regime in gliding flight, the Reynolds number and angle of attack for the optimization is Re = 1 × 10 3 and α = 0 • , respectively. The flow governing equations and the adjoint equations are calculated through a time-marching finite difference method. The minimizations of the objective function are carried out by using the steepest descent algorithm.
The third objective of the present paper is to evaluate the hydrodynamic moment of force and the attitude stability of the optimal airfoil for ultra-low Reynolds number regime in gliding flight is investigated as well. There are two applied problem settings. One is that the angle of attack of the optimal wavy airfoil, the original corrugated airfoil, and NACA2408 airfoil are fixed to discuss the hydrodynamic moment of force and statically stability. Another is that the angle of attack of the optimal wavy airfoil is passively changed by the fluid force to study the attitude stability.

Computational Methodology
In this study, there are two applied problem settings. One is that the angle of attack of the airfoil is fixed to investigated the improvement in aerodynamic performance for a series of wavy airfoil profiles generated by modifying the number of waves and the width of the ridge and valley in the chord direction and optimize the waves for high aerodynamic efficiency using a continuous adjoint method. The computational theory for the first problem setting is introduced in Sections 2.1-2.3. Another is that the angle of attack of the airfoil is passively changed by the fluid force to analyze the attitude stability. For the latter, a moving grid technical adapted to the moving airfoil surface is applied. The computational theory for the second problem setting is introduced in Section 2.4. More details of computation of the attitude stability of airfoil are described in our previous studies [28].

Flow Equations
This section presents the flow governing equations for numerically calculating the flow field over the airfoil. Consider a domain Λ, with boundary Θ that is spatial and temporal and occupied by a fluid of density ρ and dynamic viscosity µ. The spatial and temporal coordinates are denoted by x and t. For the incompressible Newtonian fluid, the continuity equation and the Navier-Stokes equations can be represented as where u means the velocity, p denotes the pressure, and E is identity tensor. The boundary conditions are either on the flow velocity or stress. Both Dirichlet and Neumann-type boundary conditions are considered: Here, n is the unit normal vector on the boundary Θ. Θ g and Θ h are the subsets of the boundary Θ. The boundary Θ can be further split, such as Θ U , Θ D , and Θ S are the upstream, downstream, and lateral boundaries, respectively, and Θ B is the body surface. The initial velocity condition can be written as Here, u 0 indicates divergence free. The drag and lift force coefficient, (C d , C l ), on the body are calculated using the following expression: The time-averaged drag and lift coefficients can be computed from: Note that the time should start at t 0 until the flow is fully developed and the time interval for averaging T should be large enough so that the mean value is stationary.

The Continuous Adjoint-Based Aerodynamic Optimization Method
The continuous adjoint-based aerodynamic optimization method is presented in this section. Let Θ B be the segment of the boundary Θ, whose shape is to be determined by the set of parameters ω = (ω 1 , . . . , ω m ). Let I c (Ω, ω) be the objective, in which the independent variables includes the flow variables Ω = (u, p) and shape parameters ω. Thus, minimizing the objective function I c (Ω, ω), the shape parameters ω are determined in the optimization procedure. When ω makes a small change, a corresponding small change of I c appears. Namely, The flow Equations (1) and (2) can be written as = ( u , p ) = 0, which expresses the dependence on Ω and ω and seems to exist as the constraint conditions for the objective function I c = I c (Ω, ω). Then, (Ω, ω) = 0 can be written. Thus, δΩ can be determined from An augmented objective function is constructed to convert the constrained problem into an unconstrained one. The flow equations are augmented to the original objective function by introducing a set of Lagrange multipliers or adjoint variables ψ = (ψ u , ψ p ): If the flow variables, Ω, exactly satisfy the flow Equations (1) and (2), the augmented objective function (11) degenerates to the objective function I c . Then, the first variation of the augmented objective function (9) can be written as follows: where ∂I ∂Ω ∂I ∂ω ∂I ∂ψ Moreover, since the variation in δ is zero, it can be multiplied by ψ and subtracted from the variation δI without changing the result. Thus, Equation (12) can be written as Here, ψ should be carefully selected and satisfy the adjoint equations: Substitution of the adjoint Equation (17) into Equation (16) results in the following equation for I and ω. where From the above equation, an optimal solution of shape parameters can be obtained as the gradient of the augmented objective function reaches 0, i.e., δI = 0. Thus, the gradient G given by Equation (19) is utilized to find the optimal shape parameters. Since the gradient G is independent of δΩ, the gradient of I with respect to an arbitrary number of design variable can be determined with the need for additional flow field evaluations. The main cost is in solving the adjoint Equation (17). In general, the adjoint problem is about as complex as a flow solution. Once Equation (19) is obtained, the gradient G can be fed into any numerical optimization algorithm to obtain an improved design.

The Adjoint Equations
The equation and the boundary conditions for adjoint variables can be carried out when the representation written in Equation (17) is set up. Thus, the adjoint equations can be written as: The adjoint Equations (20) and (21) are a set of coupled linear partial differential equations. Unlike the flow Equations (1) and (2), the equations for the adjoint variables are posed backward in time. The boundary conditions on the adjoint variables are where s = {uψ u − Eψ p + µ[∇ψ u + (∇ψ u ) T ]} · n, Θ U means the upstream boundary, Θ D means the downstream boundary, Θ S means the lateral boundaries, and Θ B is the body surface. The terminal condition on the adjoint velocity is given by:

Attitude Stability
The basic equations are the continuity and Navier-Stokes equations for an incompressible Newtonian fluid. A non-dimensionalization is applied to all variables by means of the streamwise velocity U uni and chord length L. Since the boundary-fitted-grid is employed in our computation, a general curvilinear coordinates has to be applied and is represented as ξ, η, ς. Based on the references [50,51] where J represents the Jacobian of the coordinate transformation, U k means the contravariant velocity, Re means the Reynolds number, as in ρU uni L/ν, V k denotes the moving speed component of the grid point in the general coordinate. In the pitching direction, the angle of attack variation because of the fluid force can be written as: Here, I m means the inertia moment, C M and C M0 are the moment coefficient of the fluid force and the support moment coefficient in the pitching direction, where the support point of the airfoil is at the 25% chord length from the airfoil leading edge, S denotes the effective area, and α means the angle of attack, which is opposite to the direction of the fluid moment. There is therefore a negative sign is applied on the left-hand side of Equation (29).

Shape Parameterization
The original corrugated airfoil used in this study is the same as our previous study used in numerical simulation [28]; see Figure 1. The position x c treated in this study is a coordinate system along the chord direction, and the position of the leading edge is defined as the origin x c = 0, and the position of the trailing edge is standardized as x c = 1. In order to try to explore a new form of the corrugated airfoil for gliding at ultra-low Reynolds number, as well as reduce the design variables for shape optimization, the width of the ridge or valley in the chord direction for the corrugated airfoil is expressed by the wavenumber of the sine wave. This kind of corrugated airfoil (i.e., a wavy airfoil) has different ridge and valley widths on the leading and trailing edge regions; see an example in Figure 2. Therefore, the wavenumber n(x c ) of the wavy airfoil is a model with the position x c (0 ≤ x c ≤ 1 ) in the chord direction as the independent variable. In this study, the wavenumber of the wavy wing is expressed by a quadratic function: Here, A, B, and C are constants and are determined to satisfy the following conditions: a is the wave number at the leading edge, b is the wave number at the trailing edge, and c is the position of the quadratic function axis. From the above, the wavenumber of the corrugated wing can be expressed using the parameters a, b, and c as follows: Next, the height of the wave for the wavy airfoil remains to be decided. Since the original corrugated airfoil has a shape with the envelope of NACA2408 airfoil as shown in Figure 3, the height of the wave is also determined by using NACA2408 airfoil as the envelope, see Figure 2. Let y e be the thickness of NACA2408 airfoil at the position of x c in the chord direction. The thickness of the wavy airfoil is fixed at 1% of the chord length in this study. Finally, the height of the wavy airfoil y w can be expressed as: When y w has an inflection point, the wavy airfoil will not use the NACA2408 airfoil as the envelope, as shown in Figure 4. Thus, this case of existing inflection point has to be rejected. Substitute Equation (32) into Equation (33) and take the first derivative of y w /y e identically equalling to zero, leading to the condition expression for the parameters a, b, and c that must satisfy as follows: As a result, based on the above discussion, the shape of the wavy airfoil handed in this study is determined.

The Objective Function
When a relatively large lift coefficient and low drag coefficient are both achieved, we term this airfoil as high performance airfoil. Thus, the time-averaged aerodynamic efficiency is employed as the objective function for our optimization:

Computational Domain, Mesh, and Boundary Conditions
The computational domain and boundary conditions for corrugated airfoil and NACA2408 airfoil are shown in Figure 5, where L is the chord length of the airfoil. A Cartesian coordinate system is used to define x in the mainstream direction and y in the vertical direction. Due to the H-type boundary-fitted grid employed, a general curvilinear coordinate system has to be applied and represented as ξ and η for all computations, in which the direction following the mainstream direction is denoted as ξ, and the direction away from the surface of airfoil is defined as η. The computational mesh near the airfoil surface is shown in Figure 6. Through using the mesh moving scheme on the initial grid points, the different geometries meshes are created. In this scheme, the mesh in the flow domain is modeled as a linear elastic solid and is deformed to accommodate the new geometry of the airfoil [52][53][54]. When the angle of attack of the airfoil is passively changed by the fluid force, transfinite interpolation method is used to regenerate the grids on the moving airfoil surface [55]. The validation of the mesh convergence for evaluating the adequacy of the spatial resolution was carried out in our previous studies [28].  The computational domain and boundary condition are showed in Figure 5. More details are introduced in our previous article [28]. The computational parameters are summarized in Table 1. Firstly, in order to to conduct a systematic investigation to study the improvement in aerodynamic performance for a series of wavy airfoil profiles generated by modifying the number of waves and the width of the ridge and valley in the chord direction, the direct numerical simulations of the flow around different wavy airfoils are conducted with the computational parameters of Case_1 in Table 1. Next, for optimizing the corrugations for high aerodynamic efficiency using a continuous adjoint method, the computational parameters of Case_2 in Table 1 are used. Finally, to evaluate the attitude stability of the optimal wavy airfoil generated by using our shape optimization procedure, we perform the simulations with the computational parameters of Case_3 given in Table 1.

Implementation of the Numerical Method and Optimization Procedure
A finite difference method is utilized to solve the flow and adjoint equations. Since the properties of the adjoint equations are very similar to the flow equations, numerical schemes for flow solvers can be applied to the adjoint solvers as well. In this study, applying a cell-centered, collocated arrangement the flow equations and adjoint equations are discretized, in which all physical variables and the corresponding contravariant components are located at the cell center and cell-face center, respectively. A 2-order spatial central finite difference discretization technique is employed. For coupling the pressure field and the continuity equation, the fractional method is applied [56]. For the time marching of the flow and adjoint equations, the 2-order Adams-Bashforth method is applied to the convective terms, the 2-order Crank-Nicolson method is used to the viscous terms in order to eliminate the viscous stability constraint, and the backward Euler method is utilized for the pressure term. For the calculation, the Poisson equation is the most time-consuming procedure and calculated through the residual cutting method [57]. This numerical technique used in this study has been validated extensively in several turbulent flows [28,[58][59][60][61]. For the problem setting where the angle of attack of airfoil is passively changed by the fluid force, the 2-order Runge-Kutta method is utilized for the time marching of Equation (29).
The cycle of the aerodynamic optimization starts with an initial geometry. Then, the structured meshes for the initial geometry is produced. The numerical calculation of the flow field is conducted enough times until the flow fields were fully developed and the aerodynamic coefficients exhibit adequate convergence. The adjoint variables are calculated using the date from the last 1000 steps solution of the flow field. Then, the gradients G in Equation (19) for the objective function (36) can be obtained using the solutions of the flow governing equations and adjoint equations. Through a numerical optimization algorithm, the gradient G is employed to obtain an improved design. In order to accommodate a new airfoil profile, the structured mesh is adjusted through a mesh moving scheme. In Figures 7 and 8, the details of an aerodynamic optimization procedure and an example of the iteration history of the present objective function are showed, where the angle of attack is α = 0 • and the values of the design variables for the initial geometry are ω * = (0.60, 3.0, 0.25). Table 2 shows the change in the value of C L /C D before and after the shape optimization procedure. An improvement in the lift-to-drag ratio of approximately 18% is confirmed after a shape optimization procedure.  Table 2. Change in time-averaged lift-to-drag ratio value by shape optimization procedure.

Number of Iterations Step
Time-Averaged Lift-to-Drag Ratio Value 0 1.38 13 1.63

Results and Discussion
All calculations were computed on an NEC SX-8R supercomputer of Cybermedia Center, Osaka University. The validations of our numerical solvers for the flow-through two-dimensional corrugated airfoil at a fixed angle of attack or the angle of attack passively changed by the fluid force were conducted in our previous article [28]. All numerical calculations are conducted enough times until the flow fields were fully developed and the aerodynamic coefficients exhibit adequate convergence. All results are collected by time-averaging.

Effect of Waves
The improvement in aerodynamic performance for a series of wavy airfoil profiles generated by modifying the number of waves and the width of the ridge and valley in the chord direction is explored in this section. The Reynolds number and the angle of attack applied for this investigation is 1000 and ranges from −2 • to 5 • with the increments of 1 • , respectively. The details of the computational parameters are shown in Case_1 of Table 1.
Performances and flow characteristics of various wavy airfoil profiles are showed. For comparison, the performances and flow characteristics for the original corrugated and NACA2408 airfoil profiles are listed as well. The aerodynamic performance and profile of the various airfoil geometries are summarized in Table 3 for the angle of attack of 0 • . Figures 9-14 correspondingly show the flow characteristics such as distribution of mean pressure isolines around the various geometries at the angle of attack of 0 • . Cases are labeled as Wavy_x for the explored wavy airfoil profiles, where x varies from 1 to 4. Waves are mainly placed on the different regions of the airfoil, such as on the front part with respect to the case Wavy_1, on the rear part with respect to the case Wavy_2. The width of the ridge and valley in the chord direction and the number of waves are varied, such as the small ridge and valley and many waves for the case Wavy_3, the large ridge and valley and few waves for the case Wavy_4. All of the wavy airfoil shapes and the original corrugated airfoil profile not only have a higher lift but also lower drag compared to the NACA2408 airfoil profile, which is the same as conclusions of several investigations of the corrugated dragonfly wings [18,[21][22][23][24][25][26]. Interestingly, in comparison with the original corrugated airfoil profile, Wavy_1 and Wavy_3 have more distinct waves and show a lower aerodynamic performance (in terms of lift-to-drag ratio) and Wavy_2 and Wavy_4 have less distinct waves and show higher aerodynamic performance. In particular, Wavy_3 has the maximum waves and shows the lowest aerodynamic performance, Wavy_4 has the minimum waves and shows the highest aerodynamic performance. For the drag, from the original corrugated airfoil profile to the all wavy airfoil shapes, even if the corrugations and waves on the airfoil are different and changed, there is almost nothing changed, i.e., the value at about 0.109 except Wavy_4. Thus, the lift coefficient greatly affects the aerodynamic performance compared to the drag coefficient for wavy airfoil shapes. Table 3. The time-averaged aerodynamic coefficients and its breakdown for the pressure and friction stress contributions of the various geometries at Re = 10 3 and α = 0 • . The superscripts p and f denote the pressure and friction stress components, respectively.   Table 3 also shows the breakdown of lift and drag coefficients for the contributions of the pressure and friction stress in various geometries. The superscripts p and f denote pressure and friction stress contributions, respectively. For instance, the pressure stress contribution of the lift is represented as C p L . It can be seen that the breakdown components of the lift coefficient with the different wavy airfoil shapes are almost the same, while the breakdown components of the drag coefficient concerning the different wavy airfoil shapes is different. The first, for the lift coefficient, the pressure stress component concerning all wavy airfoil shapes, the original corrugated and NACA2408 airfoil profiles are significantly dominant, whereas the frictional components are extremely minor. Through investigations of the pressure distribution of the wavy airfoils in Figures 11a, 12a, 13a, and 14a, It can be confirmed that the negative pressure generated at the valleys of the corrugated airfoil contributes to the increased lift. This also corresponds to the large flow acceleration on the blade upper surface in the streamline distribution, as can be seen in Figures 11b, 12b, 13b, and 14b. Secondly, for the drag coefficient, the pressure stress component accounts for about 75% for Wavy_1, Wavy_3, and the original corrugated airfoil, while the frictional component accounts for about 70% for Wavy_2, and Wavy_4, and the NACA2408 airfoil. That is, Wavy_1 and Wavy_3 have more distinct waves and are akin to appear as a corrugated airfoil for the drag, see Figures 9, 11 and 13. However, Wavy_2 and Wavy_4 have less distinct waves and are akin to appear as a smoothed surface NACA2408 airfoil for the drag, see Figures 10, 12 and 14. From Figures 11a and 13a, we seem to confirm that the reason why the pressure stress component of the drag coefficient for the wavy airfoils (Wavy_1 and Wavy_3) is dominant is that the downstream side and upstream side are low pressure and high pressure, respectively. On the other hand, the reason why the frictional stress component of drag for the wavy airfoils (Wavy_2 and Wavy_4) is large is that large airfoil surface areas face the external flow, see Figures 12a and 14b. In Figures 11b and 13b, it can be seen that there is a distinct trapped vortex in each cavity of the wavy airfoil which causes the main external flow to pass the airfoils without facing the surface area of the airfoil. However, compared to the streamlines around Wavy_1 and Wavy_3 airfoils, the large surface area faces the external flow for Wavy_2 and Wavy_4 airfoil In order to discuss the effect of the angle of attack significantly on the aerodynamic performance (in terms of lift-to-drag ratio) with respect to various airfoil geometries, Figure 15 shows the aerodynamic performance of various wavy airfoil shapes, the original corrugated and NACA2408 airfoil profiles at Re = 1000 as the angle of attack is varied from −2 • to 5 • . Above all, we could confirm that the angle of attacks has almost no effect on the aerodynamic performance of various geometries, apart from a small region of negative angle of attacks for Wavy_3 and Wavy_4. Therefore, it is completely acceptable that antecedent computations are performed at an angle of attack α = 0 • to evaluate the influence of airfoil geometries on the aerodynamic performance. From the results of Figure 15 and Table 3, we can confirm that the geometry shape of Wavy_4 exhibits the highest aerodynamic performance. Thus, in order to obtain as much as possible the highest global aerodynamic performance wavy airfoil shape, the subsequent aerodynamic optimization procedure begins with the highest aerodynamic performance airfoil shape Wavy_4 as the initial geometry.  Figure 9. Flow past the NACA2408 airfoil (Table 3) Figure 11. Flow past the shape Wavy_1 (Table 3) Figure 12. Flow past the shape Wavy_2 (Table 3) Figure 13. Flow past the shape Wavy_3 (Table 3) Figure 14. Flow past the shape Wavy_4 (Table 3)

Shape Optimization
After an investigation of the effect of waves on an airfoil with the envelope of NACA2408 airfoil, we carry out shape optimization for high aerodynamic performance (in terms of time-averaged lift-to-drag ratio) using a continuous adjoint method. The Wavy_4 airfoil is selected as the initial profile guess for aerodynamic shape optimization in order to obtain as much as possible the highest global aerodynamic efficiency wavy airfoil shape. The Reynolds number is 10 3 and the angle of attack is 0 • . Table 4 presents the aerodynamic performance and profile of the optimal airfoil obtained applying our optimization procedure for the aerodynamic efficiency maximization. For comparison, the data for the initial geometry Wavy_4 and the original corrugated airfoil are also listed. It can be found that, from Table 4, the shape of the optimal airfoil is similar to a curved plate wing without waves on the airfoil even though the initial guess has waves on the airfoil. Compared to the initial geometry and the corrugated airfoil, the optimal airfoil is associated with higher aerodynamic efficiency, as well as higher lift coefficient. The optimal airfoil shape corresponds to an increase of 60% over the aerodynamic efficiency (lift-to-drag ratio) from the initial geometry, an increase of 132% concerning the corrugated airfoil. The corresponding increase in lift is 62% with respect to the initial geometry and 120% with respect to the corrugated airfoil. However, in comparison with the initial airfoil geometry, only a very little improvement of the drag for the optimal airfoil appears. Figure 16 shows the distribution of pressure isolines and streamlines for the flow past the optimal airfoil shape. Here, for each color of the streamline and pressure isolines, the threshold values are the same as those in Figures 9-14 for the color bar. From the results in Table 4, the breakdown of the lift coefficient for the optimal airfoil is dominated by the pressure component, while the breakdown of the drag coefficient is dominated by the frictional contribution. These reasons can be found from the distribution of pressure isolines and streamlines in Figure 16. The optimal airfoil is similar to a curved plate and represents the upward convex. The acceleration of the flow near the upper surface of the optimal airfoil is therefore remarkable, it is considered that the pressure distribution on the upper surface of the blade is reduced and the lift coefficient is greatly earned. Since the optimal airfoil does not have distinct corrugations on the airfoil, there is no trapped vortex in each cavity of the corrugations which limits the surface of the airfoil directly facing the external flow. From this point of view, it can be seen that the shape of the airfoil obtained by optimizing the wavy airfoil model does not have the aerodynamic characteristics of a corrugated airfoil but have the aerodynamic characteristics of a smoothed-face airfoil.  Figure 16. Flow past the optimal airfoil (Table 4) at Re = 10 3 and α = 0 • : (a) time-averaged pressure isolines; and (b) time-averaged streamlines. Table 4. Mean aerodynamic coefficients and its breakdown for the pressure and friction stress components of the initial airfoil geometry and optimal airfoil for aerodynamic efficiency maximization at Re = 10 3 and α = 0 • .

Hydrodynamic Moment and Attitude Stability
Using the airfoil shape obtained by the shape optimization for the aerodynamic efficiency maximization performed above, numerical calculations related to hydrodynamic moment of force and attitude stability are conducted. The details of this computation are described in our previous article [28].

For the Problem Setting of Fixed Angle of Attack
First, for the problem setting of a fixed angle of attack, the moment coefficient CM because of the fluid force for the optimal airfoil is evaluated, in which the support point is at 25% chord length from the airfoil leading edge. The time-averaged data of the static moment coefficient are collected after the value exhibits adequate convergence. The mean static moment coefficient of the optimal airfoil at Re = 4000 with varying angle of attack is showed in Figure 17, where −6 • ≤ α ≤ 6 • with increment of 1 • . We know that, if an airfoil is stable in the range of α, as α increases the hydrodynamic moment of force for the airfoil should increase as well, moreover, the direction of the increase of the hydrodynamic moment of force should be the opposite to the direction of increase of α. In Figure 17, it can be seen that the slope of the moment coefficient C M of the fluid force concerning the angle of attack α is negative in the range of −1 • ≤ α ≤ 6 • , i.e., the angle of attack in this region is statically unstable for the optimal airfoil.  Figure 18 shows the distribution of the lift coefficient with respect to dimensionless time at Reynolds number 4000, the angle of attack −2 • . After the dimensionless time t * = 20, the evolution of the lift coefficient becomes quasi-stationary, but it can be seen that it oscillates periodically. In Figure 19, the distribution of vorticity around the optimal airfoil at Re = 4000 and α = −2 • , and dimensionless time t * = 28 is showed. From Figure 19, it can be seen that vortices are periodically discharged from the leading edge of the optimal airfoil, which may indicate the reason why the lift coefficient oscillates greatly and periodically. After further computation, it has been confirmed that the vortex shedding at the optimal airfoil leading edge is in the range of the angle of attack −6 • ≤ α ≤ −2 • . Actually, the range of α having the vortex shedding exactly corresponds with the range of the angle of attack with the hydrodynamic moment fluctuating. Based on the above discussion, as a conclusion, the optimal airfoil is either statically unstable or statically quasi-stable.

For the Problem Setting of Angle of Attack Passively Changed by the Fluid Force
In order to balance the initial hydrodynamic moment of force because of the fluid force acting on the airfoil at an initial α, a constant moment of force M 0 is used, where the rotating support point is at 25% chord length from the leading edge and the initial α is at 0 • . We observe the time evolution of the objective airfoil's α after an initial disturbance, i.e., ∆α = +0.001 • , is imposed. More details of the computational setup are introduced in our previous article [28]. Figure 20 shows the time evolution of α for the optimal airfoil, the original corrugated and NACA2408 airfoil at Re = 4000 with the initial α = 0 • . From this figure, the time evolution of the angle of attack for the corrugated and NACA2408 airfoil is relatively stable and its amplitude is increasingly large. For the optimal airfoil, the angle of attack keeps the initial angle of attack with an initial disturbance, which indicates it is neutral with respect to the moment of fluid force. However, as the dimensionless time pass, the angle of attack increases quickly and diverges in a positive direction. From the above results, coupled with the distribution of the moment coefficient CM of the fluid force concerning the fixed angle of attack in Figure 17, it can be confirmed that the optimal airfoil is destabilized.

Conclusions
In this study, the effect of the number of waves and the width of the ridge and valley in chord direction for a wavy airfoil was investigated at the angle of attack of 0 • and Reynolds number of 10 3 . All these wavy airfoil geometries have a higher lift, as well as, lower drag in comparison with that of NACA2408 airfoil profile. Interestingly, compared to the original corrugated airfoil profile, the wavy airfoils that have more distinct waves show a lower aerodynamic performance (in terms of lift-to-drag ratio) and the wavy airfoils that have less distinct waves show higher aerodynamic performance. In particular, the wavy airfoil that has the maximum waves shows the lowest aerodynamic performance and the wavy airfoil that has the minimum waves shows the highest aerodynamic performance. For the breakdown of the lift and drag coefficient concerning the pressure and friction stress contributions, the pressure stress contributions are significantly dominant for all wavy airfoil shapes concerning the lift, whereas the drag coefficient concerning the different wavy airfoil shapes are different, such as the pressure stress component accounts for about 75% for the wavy airfoils that have more distinct waves, while the frictional stress component accounts for about 70% for the wavy airfoils that have less distinct waves. From the distribution of pressure isoline and streamlines around wavy airfoils, it is confirmed that the pressure stress contribution of the drag are dominant due to high pressure on the upstream side and low pressure on the downside; the frictional stress contribution of the drag is dominant because of the large surface area facing the external flow. An adjoint-based aerodynamic shape optimization method was utilized to find as much as possible the highest global aerodynamic efficiency wavy airfoil shape at Reynolds number 10 3 and angle of attack 0 • . Interestingly again, the shape of the optimal airfoil is similar to a curved plate without distinct waves on the airfoil even though the initial guess has waves on the airfoil. The optimal airfoil shape corresponds to an increase of 60% over the aerodynamic efficiency from the initial geometry and increase of 132% concerning the original corrugated airfoil. The corresponding increase in lift is 62% concerning the initial geometry and 120% concerning the original corrugated airfoil. The optimal airfoil has an approximate drag coefficient compared to the initial geometry and the original corrugated airfoil. The hydrodynamic moment of force and attitude stability of the optimal airfoil was investigated as well. Concerning an fixed angle of attack, the optimal airfoil is statically unstable in the range of the angle of attack from −1 • to 6 • , statically quasi-stable from −6 • to −2 • , where the vortex is shedding at the optimal airfoil leading edge. Concerning an angle of attack passively varied due to the fluid force, the optimal airfoil keeps the initial angle of attack value with an initial disturbance, then quickly increases the angle of attack and diverges in the positive direction. From the above, it is confirmed that the optimal airfoil is destabilized.