Numerical Study of Airfoil Selection and Analysis of 3D Flow Phenomenon past Finite-Span Wings for Small UAVs

Small unmanned aerial vehicle (SUAV) is an unmanned aircraft vehicle (UAV) that flies at an altitude of lower than 1,100 m from the ground, has a maximum gross takeoff weight of 10 kg, and a flight speed of less than 50 m/s. One of the design factors for the small UAV design with a fixed-wing propeller is the airfoil selection. The selection of an airfoil profile using aerodynamic concepts leads to a performance coefficient that determines the selected airfoil’s sustainability and efficiency. The coefficients used are CL, CD, and CM. Numerical studies were carried out using Computational Fluid Dynamics using XFLR5 and ANSYS Fluent 19.1 software to evaluate airfoils in 2D and evaluate the phenomenon of induced drag on the wings in 3D. Airfoil selection was made on five types of airfoils: AH 83-150 Q, E399, E431, E715, and E662. The coefficients of CL, CD, and CM were obtained by varying α. 3D analysis of selected airfoil geometry with finite span. Simulation of steady conditions using Reynolds-Averaged-Navier-Stokes (RANS) in the Spalart-Allmaras turbulent model with variations of α = 0 ◦ , 8◦ , 12◦ , and 16◦ . The post-processor visualized the flow around the wing with pressure contours, velocity pathlines, and tip vortices. The analysis was carried out on the aerodynamic coefficients of CL, CD, CM, and CMr with α variation on the finite span wing. Based on the research, the results showed that the selected airfoil was E431, the aerodynamic performance of the CL, CD, CL/CD, CM, and CMr wings. In addition, information was also obtained regarding a decrease in the pressure difference between the upper surface and lower surface of the wing with an increasing span, 3D streamline, the extent of the contour of the vorticity magnitude, and a streamline on the wingtip on the upper surface and lower surface of the wing.


Introduction
Small Unmanned Aerial Vehicle (SUAV) is an unmanned aircraft vehicle (UAV) that flies at an altitude of lower than 1,100 m from the ground, has a maximum gross takeoff weight of 10 kg, and a flight speed of less than 50 m/s. The uses of small UAVs include conducting investigative surveys by taking pictures of the surrounding area, geomagnetic surveys, and collecting meteorological data.
One of the important design factors in designing a UAV with a fixed-wing and propeller propulsion is the airfoil selection. Airfoil selection takes into account the characteristics of the flight behavior and the type of UAV model selected. Airfoil profile selection uses basic aerodynamic concepts that refer to a performance coefficient that determines the selected airfoil's sustainability and efficiency. The coefficients used are C L (lift coefficient), C D (drag coefficient), and C M (pitching moment coefficient). The three coefficients are related to flight speed, wing area, wing chord, angle of attack, aspect ratio, profile shape, and Reynolds number [1].
Airfoil selection was done by comparing airfoil shape parameters, value parameters, and curves from the required performance coefficients C L , C D , and C M . The airfoil shape parameters compared include airfoil thickness and trailing edge shape. The parameter values and curves of the performance coefficients of C L , C D , and Cm against the variation of α obtained from running computational fluid dynamics (CFD) software are plotted in several graphs. The main performance analyzed is presented in the graph C L /C D against α, C L against α, and C M 1 4 chord against α. The airfoil selection was carried out in research on UAV concept design for an exploration mission in Colombia by Rocha and Solaque using the CFD method with XFLR5 software. The selection was made by comparing the graphs of the results of running XFLR5 on airfoils E168, E197, and E214. Airfoil E197 was used because it has a low C D , high C L with a slow stall, but a high C M value is obtained, which results in negative rotation [1]. Airfoil selection was also carried out in UAV design research to improve reliability. Airfoil selection uses VSAERO and PANUKL 2002 software on airfoils ILL415, ILL517, ILL515, DoA5, and NACA 2415. Airfoil NACA 2415 was used because it has better efficiency and pitching moment from airfoil charts [2].
Airplanes fly in subsonic flow. The drag on the wing with infinite span is caused by the profile drag and induced drag. The profile drag is the drag caused by the wing-body. The profile drag is the amount of surface friction resistance caused by the shear stress acting on the surface coupled with the pressure drag caused by the pressure imbalance in the flow direction, which causes the boundary layer to be separated. Induced drag is a type of drag caused by a pressure imbalance on the tip of the wing with a finite span between the top of the surface (suction area) and the bottom of the surface (pressure area). The pressure imbalance is necessary to generate positive lift, but near the wingtips on the lower surface (high air pressure) towards the upper surface (lower pressure), resulting in a streamline to curl. The three-dimensional flow results in the formation of a vortex, which alters the flow and causes a velocity component in the wing's downward direction called downwash. The forced flow pattern causes the velocity to be non-downward relative to any section of the wing's airfoil section, which reduces the initial α. The lift vector is tilted backward, and a component of the force that opposites the thrust force appears, which is called induced drag. Reducing the size of the vortex tip and minimizing induced drag is very important for aeronautical engineers [3].
The phenomena of tip vortex and induced drag occur on airplane wings. After 2D analysis for airfoil selection, 3D analysis was carried out to determine the characteristics of the induced drag phenomenon on wings with a certain span length. A 3D study examining the effect of induced drag was carried out by Panagiotou [3]. Panagiotou made six types of winglet configurations and then compared the six configurations to obtain the best aerodynamic efficiency using CFD. One of the six wing configurations is a wing without winglets, and the other five configurations use a blended winglet with variations in cant angle. The study used an airfoil PSU 94-097. Design and 3D computation were carried out with the ANSYS CFX flow solver on a small mesh with 3,000,000 nodes with the Spalart and Allmaras turbulent models. Inlet conditions with a speed of 140 km/h, flying altitude 2000 m. The computations were carried out at Re = 1, 2 × 10 6 and 4, 8 × 10 5 at the root and tip of the winglets with α = 8 • to 16 • in 4 • increments. The optimal configuration selection was done by comparing the C L , C D , and C M r coefficients on the different wing α. Panagiotou, et al. also compared the vortex in the downstream wing tip at x/c = 0.2 where C w,t is the root chord winglet, and x is the freestream velocity axis.
Based on the studies that have been carried out in previous studies, this research will focus on selecting the appropriate type of airfoil for SUAV. The first step in the 2D selection of airfoil types for wings with a thickness of 14.5% to 15.5% with performance parameters C L , C D , and C M in α variation. The next step, perform a 3D analysis of the wing with one of the selected airfoils. Wings with certain span lengths were analyzed to study the induced drag phenomenon through flow visualization and evaluation of C L , C D , C M , and C M r performance coefficients at certain α variations.

Research Description
Numerical studies were carried out using Computational Fluid Dynamics using XFLR5 software to evaluate airfoils in 2D and ANSYS Fluent 19.1 to evaluate trailing edge and induced drag phenomena on the wings in 3D. The stages of the research carried out were divided into two stages, namely airfoil selection and aerodynamic performance analysis and visualization of the flow through the UAV wing at variations of α = 0 • , 8 • , 12 • , and 16 • .

Airfoil Selection
The stages of research on 2D airfoil selection using XFLR5 software by downloading airfoil data from the University of Illinois at Urbana-Champaign (UIUC) airfoil data site website [4]. The airfoil that will be selected in this study is shown in Table 1 and Figure 1. The airfoil that will be selected is the result of the research by Perdana [5]. Enter airfoil coordinate data to XFLR5. Refining globally changes the number of panels to 250 points for each airfoil analyzed [6]. Then, entering the simulation conditions and running parameters (Table 2) into the XFLR5 software. The simulation conditions were obtained from Kim [7].
The results of running XFLR5 software are graphs of the performance coefficients of C L , C D , and C M with α variation. Then, evaluating airfoil selection by comparing the graphs of The airfoil criteria which will be selected namely has a high value of

Pre-processor stage of UAV Wing Aerodynamic Performance Analysis
This was done by importing airfoil profile coordinates into Autodesk Inventor software. The airfoil profile coordinates used are obtained from the results of the airfoil selection research phase. Then, making a wing geometry model according to the criteria in Table 3 and the boundary box as shown in Figure 2.
Create a meshing was started with a splitting face and then generating elements using ANSYS Meshing. The meshing was done with structured Hexa meshing. The mesh was then refined in the areas around the wing and estimated the y + value. Lastly, quality checks were conducted and meshing quality was improved.  The boundary condition settings are in Table 4. The simulation conditions were steady with the turbulence of the Spalart-Allmaras model. The materials used in the simulation were air with a density of 1.223 kg/m 3 and dynamic viscosity of 1.7887x 10 -5 N.s/m 2 . The fluid operated under the operating condition absolute pressure of 101050 Pa following Kim's research conditions [7].
The parameters shown in Table 4 on the inlet side include a turbulent intensity (I) value of 1% [8] and a length scale (l) = 1.020 × 10 −3 m. Reference value to determine the reference value for the calculation of the drag coefficient, lift coefficient, and moment coefficient, obtained from the inlet velocity condition.
The solution method used to simulate pressurevelocity coupling with the SIMPLE scheme. Spatial discretization using Second-Order Upwind. The convergence criteria used in the iteration process is 10 -6 . Simulations were carried out at α variations at 0 • , 8 • , 12 • , and 16 • .

Post Processor Stage Analysis of UAV Wing Aerodynamic Performance
The data collection procedure that will be used in the post-processor stage in the three-dimensional simulation of the UAV wing using ANSYS Fluent 19.1 software. The data taken in this simulation include aerodynamic parameters, added flow visualization in the form of vorticity magnitude contours, pressure contours, and velocity magnitudes. The post-processor stage displayed the aerodynamic coefficients:

2D Airfoil Simulation
2D airfoil simulations were performed using XFLR5 software on five types of airfoils, namely AH 83 150 Q, E399, E431, E715, and E662 to evaluate the aerodynamic coefficients of C L , C D , and C M with α of -10 • to 20 • . The simulation results are presented in graphs of C L /C D vs α, C L vs α, C D vs α, and C M vs α. The simulation was carried out at Re = 1.15 × 10 5 with M = 0.03. The airfoil selection was done by taking a high aerodynamic coefficient value of C L /C D , high C L max at high α, and C M close to a positive value. Airfoil used was E431 which has high C L /C D at α = 0 • and 4 • , namely 11.37 and 19.94, C Lmax = 1.413 at α = 12 • , and C M low at α = 0 • , namely -0.088 (Table 5).

2D Airfoil Simulation
2D airfoil simulations were performed using XFLR5 software on five types of airfoils, namely AH 83 150 Q, E399, E431, E715, and E662 to evaluate the aerodynamic coefficients of C L , C D , and C M with α of -10 • to 20 • . The simulation results are presented in graphs of C L /C D vs α, C L vs α, C D vs α, and C M vs α. The simulation was carried out at Re = 1.15 × 10 5 with M = 0.03.

Meshing Wing dan Grid Independency Test
Meshing was done by tightening the grid close to the wing and in areas with high turbulence phenomena, namely the wingtip and trailing edge. The meshing grid is shown in Figures 3 and 4. Aerodynamic Coefficient Airfoil   Grid independency test is the step to determine the optimal number of mesh for 3D wing simulation. Table  6 shows that Grid D is the optimal meshing because the increasing number of cells, the value of C L , shows a small error, namely 0.2% with a value of y + less than 1. Threedimensional simulation of wings using Grid D with the number of cells 4473748.

Validation of 3D Wing Simulation
The results of the 3D numerical modeling of the wings with the ANSYS Fluent 19.1 software were compared with the numerical modeling results with a higher degree of accuracy, namely the XFLR5 software. Comparisons were made to take the pressure coefficient (C p ) data on the wing mid-span α = 0 • with Cp airfoil data E431 α = 0 • obtained from the results of running XFLR5. This was done to find out how accurate the numerical modeling of the UAV wing was.
Validation by comparing C p is shown in Figure 5. The C p graph of the E431 airfoil has the same trendline as the C p graph on the wing midspan with the two graph lines sometimes coinciding. This shows that the 3D simulation results with Grid D have good accuracy.

Pressure Coefficient (C p ) in the Spanwise Direction
In the simulation the wings were taken at Z/s=0.25, 0.50, 0.75 and 0.95. The purpose of collecting C p data was to see how the difference and distribution of pressure coefficients along the span, especially the area near the tip, in the three-dimensional wing analysis. Figure 6 shows the distribution of C p in the variation of α with C p on the upper surface had a lower graph value than the lower surface. The C p of the lower surface for Z/s = 0.25, 0.50, 0.75, and 0.95 tends to have a slightly smaller value with the increase in the position of the span towards the tip. In contrast, the C p of the upper surface had a value that increased drastically with increasing the span tip position. Increasing α results in a more negative change in the C p of the upper surface, while the C p of the lower surface had a slightly.
The C L value is shown in Figure 7(a) which increased with increasing the value of α to α = 16 • . The maximum C L value was at α = 16 • with a value of 1.49.
The wing C D value is shown in Figure 7(b) with a trendline that increased with increasing α. The maximum CD value was at α = 16 • with a value of 0.143.
The C L /C D values indicate the aerodynamic efficiency shown in Figure 7(c). Figure 7(c) shows the C L /C D trendline which increased then decreased with increasing α. The maximum C L /C D value was at α = 8 • with a value of 17.91.
The C M value is the wing pitching moment coefficient shown in Figure 7(c). The pitching moment causes the nose to pitch down on the wing. The figure 7(d) shows the C M trendline which increased then greatly decreased and increased with increasing α. The maximum C M values was at α = 16 • with a value of -0.101.
The C M r value is the root bending moment coefficient of the wings which is shown in Figure 7(d). Root bending moment is caused by flow moving from the high pressure of the lower surface of the wing to low pressure on the upper surface of the wing. Figure 7(e) shows the C M r trendline which increased then decreased with increasing α. The maximum C M r value was at α = 16 • with a value of 3.13.   Figure 8. The pressure contour comparison was observed to know the influence of the vortex tip on the UAV wing. Air fluid flows at freestream velocity across the wing and results in a pressure difference in the area around the wing.
Fluid flows through the wing so that a stagnation point is formed in the leading edge area, which is marked with a red contour in Figure 8. The wing shape with an E431 airfoil caused the flow to cross the wing's upper sur-face faster (lower pressure) than the lower surface (higher pressure). The difference in pressure on the upper and lower surface causes the wing to have a lift. The difference in the pressure contour between the upper and lower surface was getting smaller with the increase of the Z/s value towards the wingtip due to the existence of the tip vortex and the increasing adverse pressure gradient. The smaller the difference in the pressure contour towards the wingtip makes the lift in the midspan bigger than the wingtip.

Contours of 3D Streamline and Velocity Magnitude
Streamline contours and velocity magnitude were used to visualize the fluid velocity through the wing on an E431 airfoil. This flow visualization was taken at midspan (Z/s = 0.50) and at the wingtip (Z/s = 0.95) with α = 0 • , 8 • , 12 • , and 16 • . Figure 9 shows the velocity magnitude contour on the wing's upper surface, which had a higher velocity than the lower surface of the wing. The streamline at Z/s = 0.50 had a trendline that follows the wing-body while at Z/s = 0.95 the streamline moved away from the upper surface of the wing-body at the rear and tends to be 3D streamlined in the form of a curl. This happened because, at Z/s = 0.95, it was close to the wingtip where there was a flow on the lower surface with high pressure towards the upper surface with lower pressure.
The wake area behind the wing-body at Z/s = 0.50 increased from α = 0 • to α = 16 • . The wake area behind the wing-body at Z/s = 0.95 had a small area at α=0 • , 8 • , 12 • , and 16 • .  Figure 10 shows the line of vorticity magnitude after the flow passes through the wing body α = 0 • , namely at X/c = 1.0, X/c = 1.5, X/c = 2.0, and X/c = 2.5. The vorticity magnitude value shows the vortex tip that occurred after the flow passed through the wing body. The value of vorticity magnitude had a high value at the center of rotation.
The vorticity magnitude of the center of rotation increased at X/c = 1.0 to X/c = 1.5 but continued to decrease at X/c = 2.0 and X/c = 2.5. The area of vorticity magnitude increased with increasing X/c value. Figure 11 shows the contours of the vorticity magnitude at X/c = 1.   The upper surface of the wing body α = 16 • shows the pressure contour of the leading edge along the span, which had a low pressure of about -188 Pa and the area's pressure contour along the trailing edge that had an adverse pressure gradient of about -36 Pa. The streamlines on the wingtip were different in orientation and tended to be curl-shaped. The lower surface of the wing body shows the contour on the leading edge had a high pressure of 21 Pa, with the trailing edge having a low pressure of about -17 Pa. The streamline on the wingtip was different in orientation to the direction of the freestream velocity.

Conclusions
Two-dimensional (2D) numerical analysis of airfoils with XFLR5 software and 3D wings with ANSYS Fluent 19.1 software has been performed. The airfoil selection result was E431 which had high C L /C D at α = 0 • and 4 • , namely 11.37 and 19.94, respectively, C L max = 1.413 at α = 12 • , and low C M at α = 0 • which is -0.088. The aerodynamic performance of the wing had a maximum C L at α = 16 • with a value of 1.49, the maximum C D at α = 16 • with a value of 0.143, high C L /C D at α = 8 • with a value of 17.91, C M max = 0.101, and C M r max = 3.13 which was high at α = 16 • .
The upper and lower surface pressure contours in the wingtip area had a small difference in pressure due to the tip vortex. A decrease in the pressure difference with increasing span indicated a reduction in wing lift in the area around the wingtip. Velocity magnitude at the intersection of the wing bodies shows that the wake area was greater at Z/s = 0.50 than at Z/s = 0.95. The 3D streamline shows at Z/s = 0.95 the flow tended to move from the lower surface to the upper surface of the wing body.
The extent of the contour of the vorticity magnitude increased with the increase in the value of X/c behind the wing body and α. Vorticity magnitude shows the phenomenon of trailing vortex on the wing body. Streamlines on the wingtips on the upper surface and lower surface of the wing body experienced different orientations with the flow of freestream velocity and tended to be curl-shaped.