Aerodynamics of Owl-like Wing Model at Low Reynolds Numbers

1)Department of Mechanical Engineering, Tokyo University of Science, Tokyo 125–8585, Japan 2)Department of Aerospace Engineering, Tohoku University, Sendai, Miyagi 980–8579, Japan 3)Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga, Fukuoka 816–8580, Japan 4)Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency, Sagamihara, Kanagawa 252–5210, Japan 5)Department of Information and Computer Technology, Tokyo University of Science, Tokyo 125–8585, Japan


Introduction
In recent years, engineering communities are paying attention to small unmanned air vehicles (UAVs). The UAVs are potentially used for civilian tasks 1) including transportation, communication, agriculture, disaster mitigation, environment preservation and planetary aerial exploration. 2) The UAVs typically operate in the low Reynolds number flow regime Re=Oð10 3À5 Þ because of their size and flight velocity. Under the low Re flow conditions, viscous effects of the flow around the airfoil in the aerodynamics of the UAVs 3-5) play a more important role than those in the aerodynamics of commercial aircrafts.
Recently, Nagai et al. 6,7) presented a preliminary UAV design for the aerial exploration of Mars. Their design was obtained through a multi-objective design optimization process with certain constraints. The configuration and propulsion of the UAV are a fixed-wing and direct-current motor-driven propellers, respectively. The chord and span length of the main wing and flight velocity of the model are 0.63 m, 2.45 m, and 60 m/sec, respectively. Thus, the estimated chord-based Re is approximately 2:7 Â 10 4 . The airfoil shape of the main wing is called the Ishii airfoil. 8) Although decent maximum lift coefficient C l and lift-to-drag ratio l=d have been presented using large-eddy simulations (LES) and force measurements in a previous study, 8) Nagai et al. 7,9) emphasized further improvements in the aerodynamic performance of the main wing in order to built a complete UAV system with sufficient robustness for use on Mars.
Previous investigations [10][11][12][13][14][15][16] have suggested a variety of cross-sectional airfoil shapes for aircraft with low flight speed. Particularly, Selig and Guglielmo 16) developed a high-lift, low-Re airfoil (S1223) at the Re of 2 Â 10 5 . Shyy et al. 17) presented the S1223 with a thinner thickness than the original one, improving the aerodynamic performance at the Re of 7:5 Â 10 4 . Interestingly, Liu et al. 18) mentioned that the S1223 16) has the same maximum camber line and thickness coordinates of bird airfoils (i.e., airfoils shaped like those of the seagull and merganser wing 18,19) ). Birds are experts of gliding flight in the Re range Oð10 4À5 Þ 20) in the Earth's atmosphere. Owls, [21][22][23][24] swifts, 25,26) and eagles 27) are good examples. They fly in the Re range Oð10 4À5 Þ based on the mean chord length c of the wing and the gliding flight velocity u 1 . Through adaptation to their living environment, the wings of some birds possess unique characteristics, such as leading-edge serration, velvet-like surface and feathers, a trailing-edge fringe, structural flexibility, and unique airfoil shape and planform. Until now, the roles and mechanisms associated with each of the above-mentioned features in bird flight have been actively studied. 24) Most studies have reported that leadingedge serrations, [28][29][30] feathers, 25,26) fringes, 23,24,31) the effects of wing flexibility, 24,32,33) and airfoil shape, 27,[34][35][36][37] have a positive influence on aerodynamic performance.
Very recently, Ananda and Selig 37) designed bird-like airfoils at the Re of 6:0 Â 10 4 , 1:0 Â 10 5 , and 1:5 Â 10 5 for the main wings of UAVs. Even though Ananda and Selig 37) mainly focused on the airfoil shape, their results have shown the best l=d and C l use in UAVs requiring high wing loading. Moreover, Okamoto et al. 38) studied the effects of thickness, camber, surface roughness, and leading-edge sharpness on aerodynamic performance at the Re of 1:1 Â 10 4 and 1:4Â 10 4 . They found that the roughness of the airfoil surface can increase the maximum C l and l=d. Sunada et al. 39,40) conducted a systematic comparison of the aerodynamic characteristics of many airfoil shapes at the Re of 4:0 Â 10 3 . They revealed that either an under-camber of 5% and sharp leadingedge or proper corrugation could improve the aerodynamics of a wing with a rectangular cross-section. Recently, Winslow et al. 41) investigated various airfoil characteristics at a Re of 10 4À5 using a two-dimensional (2D) Reynolds-averaged Navier-Stokes (NS) solver with a Sparat-Allmaras turbulence model and a correlation-based laminar-turbulent boundarylayer transition model. They showed that cambered flat-plate airfoils have better C l and C d characteristics than thick conventional airfoils with a rounded leading-edge when the Re becomes 10 4 . Moreover, our previous study 42) examined the aerodynamics of many airfoil shapes at the Re of 2:3 Â 10 4 using 2D laminar simulations. It has been found that the owl-like airfoil has a higher maximum l=d than that of the Ishii airfoil 8) and other airfoil shapes considered. Subsequently, we constructed an owl-like wing model based on the cross-sectional profile of the owl wing at a 40% wingspan from the root and investigated aerodynamic characteristics at the Re of 2:3 Â 10 4 using LES. 35) The results proved that the owl-like wing model reached a higher l=d where compared with conventional thin and thick symmetrical airfoils at the Re of 2:3 Â 10 4 . However, the aerodynamics associated with the owl-like wing model at other Re has remained unclear.
Based on the above-mentioned studies, the current study focuses on the aerodynamics associated with the owl-like wing model enabling a wider Re range compared to those previously studied. Specifically, the Re range of the UAV cruise flights in the Earth and Mars atmospheres and the gliding flight of birds in the Earth's atmosphere (i.e. the Re ranging from 1:0 Â 10 4 to 5:0 Â 10 4 ). This study utilizes the LES approach in order to accurately handle unsteady flow phenomena (i.e., separation, transition, and reattachment) in this low-Re flow regime. Experimental force measurements are also conducted to validate the results of LES. Although the simplified owl-based airfoil shape is analyzed this study aims to obtain a detailed understanding regarding flow structures around the owl-like wing model and the interplay between flow structures and aerodynamic performance at a low Re.

Flow conditions
This study considers that chord and freestream-based Reynolds numbers (Re) vary from 1:0 Â 10 4 to 5:0 Â 10 4 and the angle of attack (¡) changes from 0.0°to 14.0°, repsectively. This Re range covers the level flights of UAVs and the gliding flights of owls. 23) For computation, the freestream Mach number is set to 0.2. This value is low enough that the compressibility of fluid is negligible. Therefore, the flow field obtained using our simulation is considered to be similar to that with a low freestream velocity. The specific heat ratio and Prandtl number are set at 1.4 and 0.72, respectively.

Owl-like wing model construction
This study considers the same owl-like wing model as that used in our previous study. 35) The airfoil geometry is decomposed into upper and lower outlines that are constructed based on the mathematical expressions reported in Liu et al. 18) The upper (z upper ) and lower (z lower ) outlines are calculated based on the camber line (z c ) and thickness distribution (z t ) as follows: where, z upper ðxÞ, z lower ðxÞ, z c ðxÞ, and z t ðxÞ are the upper and lower surface coordinates and a function of the camber line and airfoil thickness, respectively. The Birnbaum-Glauert camber line 18) is used to extract the mean camber line (z c ) from wing surface measurements and is expressed as follows: where, c, ¼ x=c, and z c;max are the local wing chord length, the normalized chord-wise coordinate, and the maximum camber coordinate, respectively. The thickness distribution (z t ) is given as follows: Trans. Japan Soc. Aero. Space Sci., Vol. Based on Eqs. (1), (2), (3), (4), and (5), the cross-sectional shape based on the owl wing is constructed. According to previous reports 18,23,24) and assuming that an owl wing nearly possesses an elliptic lift distribution in the span-wise direction, the inner 40% of the wing from the root provides half of the lift and is selected. Figure 1(a) shows the crosssectional owl-like wing model at $ ¼ 2y=b ¼ 0:4. The maximum thickness ratio is approximately 0.055 at x=c ¼ 0:11. This value is similar to that reported in Ananda and Selig. 37) It should be noted that the thickness of the computational model near the trailing-edge (i.e. 0:9 < x=c < 1:0) is assumed to be zero, while the thickness of the experimental model is 1.25% c.

Fluid dynamic simulation 2.3.1. Governing equations
This study uses an in-house fluid dynamic solver 43) to solve unsteady flow around the wing. Three-dimensional (3D) spatially-filtered NS equations, non-dimensionalized by the freestream density, freestream velocity, and chord length of the wing, are employed as the governing equations in this study. In the non-dimensional form, the governing equations are represented as follows: where, x i , u i , q i , μ, p, e, ( ij , ij , and t denote the nondimensional forms of the positional vector, velocity vector, heat flux vector, density, static pressure, total energy per unit volume, stress tensor, Kronecker delta, and time, respectively. The Reynolds (Re), Prandtl (P r ), and Mach (Ma) numbers are non-dimensional parameters and defined as follows: where, " 1 , U 1 , a 1 , C p , and 1 represent viscosity, velocity, speed of sound, constant pressure specific heat, and thermal conductivity, respectively. The subscript 1 denotes the quantity under freestream conditions. Here, the viscosity is calculated using Sutherland's law.

Numerical methods
In NS equations, the spatial derivatives of convective and viscous terms, metrics, and Jacobian are evaluated using the sixth-order compact differencing scheme. 44,45) At the first and second points off the wing surface, a second-order explicit differencing scheme is used. This study utilizes the LES approach in order to prevent the uncertainty of the results that arise because of modeling boundary layer turbulence. Resolving the boundary layer turbulence is important to accurately simulate unsteady flow phenomena. While additional stress and heat flux terms are appended in an ordinary LES approach, they are not implemented in an implicit LES approach. 46) In this study, based on the supposition that a high-order, low-pass filter selectively damps only unresolved high-frequency waves, the implicit LES is adopted. A tenth-order filter 44,47) is utilized with a filter coefficient of 0.495. This implicit LES approach has been well validated by Visbal et al. 48,49) for several problems, and the results of the implicit LES model have agreed with the experimental data and numerical results with standard subgrid-scale models. For time integration, the second-order backward difference is used, and it is converged through five subiterations 50) of the alternating directional implicit symmetric Gauss-Seidel implicit method 51) in each time step. The computational time is 2:5 Â 10 À4 in non-dimensional time corresponding to the maximum Courant number of approximately 1.4. All computations are conducted using either the supercomputer at the Japan Aerospace Exploration Agency (JAXA) or K supercomputer at RIKEN Advanced Institute for Computational Science (AICS) in Japan.

Computational grids and boundary conditions
The computational grid around the wing is shown in Fig. 1(b). A C-type structure mesh is utilized. The grid coordinates are oriented such that ² traverses clockwise around the wing, © follows a span-wise direction, and ¦ is normal  Trans. Japan Soc. Aero. Space Sci., Vol. 63, No. 1, 2020 to the surface. The outer boundary is positioned 30c away from the wing to reduce its influence on the solution near the wing. The first grid points away from the wing surface are 1:4 Â 10 À4 c. There are 615 grid points traversing around the wing in the clockwise (²) direction, 201 points for the span-wise direction (©), and 101 points for the direction normal to the surface (¦), so that the total number is 12,485,115 points. The grid resolution of the suction side in terms of the maximum wall units in the case of a Re of 46,000 and an angle of attack of 6°are Á$ þ max 25, Á þ max 15, and Á þ max 1; where, Á$ þ is the circumferential grid spacing, Á þ is the span-wise grid spacing, and Á þ is the wall-normal grid spacing. The superscript + denotes the normalized value based on the wall unit. With these criteria, 52) near wall fully-developed turbulent structures are sufficiently resolved. At the outer boundaries, all variables are extrapolated from one point inside of the outer boundary. A no-slip and adiabatic boundary condition is adopted for the wing surface. For the span-wise direction, 20% chord length is computed with a periodic boundary condition. This boundary condition is imposed using a 10 point overlap.

Experimental setup and flow conditions
Experiments are conducted using a low-speed wind tunnel at Kyushu University in order to support and validate the numerical results of this study. The wind tunnel used is of the circuit type, and the closed test section has a cross-section of 180 mm Â 360 mm. The turbulence level is approximately 0.3% for the wind velocity of 5 m/sec. The chord and span lengths of the experimental owl-like wing model are 80 mm and 180 mm, respectively. It is made of balsa wood and has smooth contours. It should be noted that part of the airfoil shape is different from the original geometry shown in Fig. 1(a). Although the trailing-edge part from around x=c ¼ 0:9 to x=c ¼ 1:0 corresponds to feathers, and its thickness is almost zero for the computational model, the experimental model has a 1-mm thickness in the trailingedge part due to the difficulty of forming. The difference is only in the thickness and the camber line is almost the same. However, the difference in thickness is expected to influence the trailing-edge with a drag force at low ¡. A three-component force balance (model: LMC-3531-5N, NISSHO-ELECTRIC-WORKS Co., Ltd.) is used to measure the aerodynamic forces. The Re are set to be 2:3 Â 10 4 , 4:0 Â 10 4 and 5:0 Â 10 4 , and the ¡ is varied from 0.0 deg to 14.0 deg. Note that pretests have shown that measurements taken at the Re of 1:0 Â 10 4 reach the limits of the experimental set-up due to given limitations in the precision of force balance used in this study. The repeatability of the experiment has been evaluated and confirmed. Further details can be found in our previous report. 53)

Verification and validation analysis
Grid sensitivity analysis is performed in the case of a Re of 4:6 Â 10 4 and an ¡ of 6.0 deg. Table 1 shows the number of grid points in each direction. Table 1 also includes a comparison of mean lift and drag coefficients (C l and C d ) among the three grids and shows that the difference among the grids is less than 3%. Moreover, the distribution of mean surface pressure (C p ) and friction coefficients (C f ) among three grids are compared and confirmed (Fig. 2). Therefore, the intermediate grid is selected for all simulations in this study.
Next, the mean C l and C d of numerical simulations as a function of ¡ are compared with the experimental data. The bars shown in Fig. 3 indicate standard deviations in the temporal variation of the aerodynamic coefficients under a quasi-steady-state, and filled symbols with line and open symbols indicate numerical and experimental results, respectively. This shows that the computational results are in acceptable agreement with the experimental data, and depict a qualitative trend with respect to the change in ¡. However, a quantitative difference is observed in ¡ when there is a nonlinear change in C l . This discrepancy most likely comes from the difference in flow conditions at the inlet when comparing the numerical simulations and experiments. Moreover, the C d of experimental results are slightly higher than those of numerical results. This is most likely because of the difference in the thickness at the trailing-edge when comparing the computational and experimental models, in addition to the different inlet flow conditions. Those points should be addressed in the near future. It is noted that the fluid dynamic solver used in this study has been validated for flows around several wings 8,54,55) and the flat plate 56) in similar Re ranges.

Results and Discussion
3.1. Aerodynamic performance and fluid dynamics Figure 3 shows temporal and span-averaged C l and C d as a function of ¡. It is seen that C l and C d increase as ¡, except for C d at low angles of attack (i.e., < 3 deg). Since the airfoil shape is asymmetric, C l is not zero at ¡ of 0 deg. A nonlinear change in the C l curve with respect to the change in ¡ is presented for all Re. The change in Re affects ¡ where the nonlinear change in the C l curve appears (i.e., ¼ 4:5 deg for Re ¼ 1:0 Â 10 4 , ¼ 3:0 deg for Re ¼ 2:3 Â 10 4 , and ¼ 1:5 deg for Re ¼ 4:6 Â 10 4 , respectively). This nonli-  Trans. Japan Soc. Aero. Space Sci., Vol. 63, No. 1, 2020 nearity in the C l curve slope has also been observed in previous reports 8,55,57) that studied other airfoils. Mueller and Batill 57) mentioned that a leading-edge separation bubble induces such transition, which accounts for a sudden increase in C l . It should be noted that such nonlinear C l phenomenon is a well-known factor of low Re airfoils. [3][4][5] Figure 4 plots locations of the separation and reattachment points as a function of ¡. In Fig. 4, L.E. and T.E. denote the leading-edge (x=c ¼ 0) and trailing-edge (x=c ¼ 1). A distance in the x=c between the separation and reattachment points indicates the length of the LSB. For the suction side of the wing (Fig. 4(a)), the flow starts to separate from the trailing-edge side (i.e., trailing-edge separation). As ¡ increases, the separation point moves from the trailing-edge to the leading-edge. When ¡ becomes larger than 6 deg, the flow separates near the leading-edge (i.e., leading-edge separation). The separated shear layers reattach on the surface when ¡ becomes larger than 4.5 and 6 deg at the Re of 2:3 Â 10 4 and 4:6 Â 10 4 respectively, resulting in formation of the LSB. The separated shear layers at the Re of 1:0 Â 10 4 do not reattach on the surface. While for the pressure side of the wing, the separated shear layers reattach on the surface at all Re. The separation points move from the leading-edge to the mid-chord of the wing, while reattachment points move from the trailing-edge to the mid-chord of the wing. The length of the LSB in the chord-direction reduces as ¡ increases. At all Re, the flow past the pressure side is fully attached when ¡ is higher than 4.5 deg. Moreover, the separation point at the same ¡ locates at almost the same positions for all Re. Thus, the length of the LSB in the chord-wise direction on the pressure side mainly depends on the position of the reattachment point, and the length of the LSB decreases as Re increases. Now, the possible reasons of the nonlinearity in the C l curve slope for the owl-like wing model are analyzed and discussed based on the computational results corresponding to ¡ according to Figs. 3 and 4 for each Re. Figure 5 visualizes instantaneous flow structures around the wing at two ¡ for the Re of 1:0 Â 10 4 , 2:3 Â 10 4 , and 4:6 Â 10 4 . In Fig. 5, the iso-surface represents the second invariant of the velocity gradient tensor Q, 58) its color is the vorticity vector in the chord-direction velocity direction (! x ), the background contours indicate the magnitude of chord-wise velocity, and T.E. denotes the trailing-edge of the wing. The second invariant of the velocity gradient tensor Q and chord-wise direction velocity u are normalized by the freestream velocity u 1 and chord length c. The shear layers are separated, and these separated layers induce a vortex due to Kelvin-Helmholtz instability in the region of upper wing surface. The vortex induced breaks down after it travels a few percent of c. It seems that the flow transits from laminar to a turbulent state. At the ¡ corresponding to the nonlinear change in the C l curve, the location of vortex shedding from the separated shear layers is found near the upper surface of the wing. Moreover, it is seen that the turbulent vortices travel along the wing surface in the snapshot of the flow around the wing. At the same time, an increase in chorddirection velocity over the suction side of the wing is observed due to the change in the flow structures. Figure 6 shows a comparison of temporal and span-averaged C p distribution on the wing surface at two ¡ for each Re. It is found that the suction peaks are enhanced by either the leadingedge separation in the case of Re of 1:0 Â 10 4 (Fig. 5(a) and Fig. 6(a)) or generation of the LSB on the suction side in the cases of Re of 2:3 Â 10 4 and 4:6 Â 10 4 (Fig. 5(b), (c) and Fig. 6(b), (c)) when the ¡ corresponding to the nonlinearity in the C l curve slope appears. At the same time, the C p distributions over the pressure side change and increase. Therefore, the significant increase in suction peak due to either the change of type of separation or reattachment of flow near the trailing-edge and the recovery of C p on the  Trans. Japan Soc. Aero. Space Sci., Vol. 63, No. 1, 2020 pressure side are most likely responsible for the nonlinearity in the C l curve slope observed in Fig. 3. The maximum C l is approximately 1.0 for all Re considered and favorable. After ¡ corresponding to the maximum C l , the mean C l still maintains a high value although the wing is stalled (see Fig. 3). For most ¡, C l and C d become large and small, respectively, as Re increases. Moreover, due to the variation in time-averaged C l and C d , l=d also shows a nonlinear curve with respect to the change in ¡. Table 2 and Fig. 7 summarize the data of maximum l=d reported in relevant studies. 8,13,[36][37][38]41,42,55,[59][60][61][62] The results indicate that the owl-like wing model has favorable aerodynamic performance in the Re range region of interest. It should be noted that at the same Re, the experimental l=d max tends to be smaller than the computational one. This is most likely due to the slight difference in drag force data. Furthermore, the values of l=d max for the Ishii airfoil, 8) SD7003, 63) and NACA0012 55) listed in Table 2 were obtained using the current LES solver. Thus, we think that the values of l=d max in Table 2 have a certain consistency. Figure 8 shows temporal-and span-averaged chord-wise direction velocity u and surface pressure coefficient C p distribution for the maximum l=d conditions at each Re. In Fig. 8, the symbols of S and R denote the locations of flow separation and reattachment points on the wing surface, respectively. It is observed that the laminar boundary layers separate from the wing surface and separation occurs at different locations among the cases. At the Re of 2:3 Â 10 4 and 4:6 Â 10 4 , flow reattaches on the suction side, resulting in the formation of LSB. High suction peaks of C p are seen for all Re conditions. Especially, the suction peak of C p at the Re of 2:3 Â 10 4 is remarkably high. For the pressure side, C p keeps a positive value over most of the surface at all Re, hence increasing C l .
In addition, in order to explore the reasons that the owllike wing model has better performance than other Reynolds number airfoils, we conduct a comparison of flowfields around wings and C p distribution. Here we select the Ishii airfoil 8) and SD7003 airfoil 13,63) from many existing airfoils for comparison. This is because the Ishii airfoil 8) and SD7003 airfoil 13,63) are known as the airfoil shapes that have good aerodynamic performance under the low Re conditions considered and computational results 8,63) are available. Moreover, we focus on and perform this analysis for the maximum l=d condition at the Re of 2:3 Â 10 4 . Figure 9 illustrates the comparison of cross-sectional airfoil shapes. The owl-like wing model is thinner than the Ishii airfoil 8) and SD7003 airfoil. 13,63) The outline of the pressure side of the owl-like wing model shows a deep concave shape in comparison with that of the Ishii airfoil 8) and SD7003 air-  Trans. Japan Soc. Aero. Space Sci., Vol. 63, No. 1, 2020 foil, 13,63) while the outline of the suction side of the owl-like wing model is relatively similar to that of two wings. Figure 10 shows the comparison of instantaneous and temporal-and span-averaged flow structures around the owl-like wing model, around the Ishii airfoil, 8) and around the SD7003 airfoil, 63) respectively. Figure 11 plots the comparison of temporal-and span-averaged C p distribution along the airfoil surfaces. Since we consider the flow conditions corresponding to the maximum l=d condition of each wing, all wings are not at the same ¡. It is found that fundamental structures of instantaneous and averaged flow structures over the suction side of all wings look similar (see Fig. 10). In terms of the instantaneous flow structures, laminar separation occurs, the spanwise vortex is shed from the separated shear layer, and a hairpin vortex is developed due to breakdown of the spanwise vortex. Additionally, in terms of the temporaland span-averaged flow structures, LSB is observed on the suction side of all wings, the length in the chord-wise direction of the LSB is similar, and the location of LSB is different due to the different location of separation and reattachment points on the suction side. As seen in Fig. 11, a considerable difference in C p distribution is evident. In comparison with that of the Ishii and SD7003 airfoils, the owl-like wing model shows a lower suction peak and surface pressure in the plateau region and a higher pressure region on the pressure side. This is most likely due to the different airfoil shapes and value of ¡. For similar ¡ conditions, the main characteristics of C p distribution are almost the same (results are not shown

Re
Owl-like wing model(CFD, 2D) 36,42) Owl wing model with feathers(EXP) 61) Owl-like wing model(CFD, current) Other airfoils(CFD) 8,36,37,41,42,55,60,62) Other airfoils(EXP) 13,38,59) (×10 4 ) Trans. Japan Soc. Aero. Space Sci., Vol. 63, No. 1, 2020 here). We find that the shape of the pressure side of the owllike wing model has a positive influence on lift generation. Furthermore, we believe that the above-mentioned characteristics in C p distribution of the airfoil surface due to its geometry are possible reasons that the owl-like wing airfoil has better performance than other low-Re airfoils (here, the Ishii airfoil and SD7003 are considered) under the maximum l=d condition and at a Re of 2:3 Â 10 4 . 3.2. Effects of pressure-side flow separation on wing aerodynamics In this subsection, the effects of pressure-side flow separation on aerodynamics are discussed. Here, flow structures around the wing at the ¡ of 0.0 deg are adopted for analysis. As seen in Figs. 4 and 12, at this ¡, flow over the pressure side presents laminar separation, laminar-to-turbulent transition, reattachment, and generation of LSB due to a deeply concave lower surface for all Re, while there is no LSB generation on the suction side. In Fig. 12, the turbulent kinetic energy TKE is defined as ðu 02 þ v 02 þ w 02 Þ=ð2u 2 1 Þ, and TKE max indicates the maximum TKE along the pressure side of the wing surface.
For all Re, the C p on the pressure side near the reattachment point shows a positive value and gradually decreases toward the trailing-edge. Those positive C p regions generated on the pressure side enhance lift generation. Moreover, the magnitude of C p corresponding to the inside of the LSB is similar, while the magnitude of C p near the reattachment point and its position vary according to the Re. The distribution of TKE max correlates with that of C p on the pressure side. Figure 13 shows the instantaneous flow structures around the pressure side of the wing and the power spectral density (PSD) of the chord-wise direction velocity u at each Re. Note that the positions of PSD sampling points for chord-wise direction velocity u correspond to the location of TKE max with variation in the chord-wise direction. The data are ensemble-averaged in the span-wise direction. Results present the unsteady flow features including the vortex shedding from the separated shear layer and the formation of 3D turbulent vortices. The u PSD results shown in Fig. 13 show that ¹5/3 power decay can be observed at x=c ¼ 0:70 at the Re of 1:0 Â 10 4 , x=c ¼ 0:50 at the Re of 2:3 Â 10 4 , and x=c ¼ 0:35 at the Re of 4:6 Â 10 4 , and the flow around the wing from the chord-wise position is judged to be turbulent. It is clear that the flow separation over the pressure side has an impact on the time histories of C l (see Fig. 14) and the C l fluctuates at all Re. The amplitude and frequency of C l fluctuation become smaller and higher as Re increases because of the change in the vortical structures around the wing (see Fig. 13).

Conclusion
This paper investigates the flow and aerodynamic characteristics associated with the owl-like wing model under low Owl-like wing model Ishii airfoil [8] SD7003 airfoil [13,63]   Owl-like wing model Ishii airfoil [8] SD7003 airfoil [63] Suction side Fig. 11. Comparison of temporal-and span-averaged C p distribution along the wings in the case of maximum l=d at Re ¼ 23;000. Solid red, dashed black, and broken blue lines correspond to owl-like wing, Ishii, and SD7003 airfoils, respectively. Re conditions (1 Â 10 4 -5 Â 10 4 ) based on 3D LES and force measurements. The computational time-averaged C l and C d are in reasonable agreement with the experimental results.
The computational results show that the owl-like wing model attains a relatively high maximum C l and higher l=d when compared to common low-Re airfoils, such as SD7003, 63) the Ishii airfoil, 8) and other airfoils in the Re range between 1:0 Â 10 4 and 5:0 Â 10 4 . Although the flow separation and movements and formation of LSB occur on the pressure side, the concave lower surface improves the time-averaged aerodynamic performance at all ¡ in this low-Re flow regime. Moreover, the analyses reveal that the nonlinearity in the C l curve slope with respect to the change of ¡ is most likely due to the increase in the suction peaks in conjunction with either the change of type of separation or reattachment of the flow near the trailing-edge and the recovery of pressure on the pressure side. Those findings could be the basic characteristics of the owl-like wing model and used for improving the aerodynamic design of fixed-wing UAVs flying in the Re range from 1:0 Â 10 4 to 5:0 Â 10 4 . Trans. Japan Soc. Aero. Space Sci., Vol. 63, No. 1, 2020