Numerical Analysis of a Flow over a Naca0015 Airfoil Containing the Flap and Gurney Flap Elements

The present work aimed at numerical simulation of a two-dimensional, incompressible and steady-state airflow over a NACA0015 airfoil with flap and Gurney flap using the standard turbulence models k-epsilon, k-omega and versions of these models with modified constants. All meshes used were structured. Comparison of the results of the turbulence models with experimental data from the literature shows differences between the results in the leading-edge region. The differences were minimized by adjusting the contents of the turbulence models. For validation, the k-epsilon turbulence model with modified constants was used in a new simulation of a profile (NACA0012), which agrees with several experimental studies; however, it does not show better results. present work aimed at numerical simulation of incompressible steady-state airflow over a NACA0015 hydrofoil with flap and Gurney flap using the standard k-epsilon and k-omega turbulence models. The meshes used for the simulation were structured meshes created using the ICEM CFD® module. The results were compared with experimental data obtained by [9], [10] and [11]. The commercial software ANSYS CFX® was used for the computational simulations.


I. Introduction
The numerical simulation of flows is called computational fluid dynamics (CFD). The use of this tool is based on the numerical solution of the conservation equations, since no general analytical answers for these equations are known, except for very simplified flows, which generally have limited application in engineering [1]. For permanent subsonic flows of Newtonian fluids, equations 1 and 2 represent the equation for conservation of mass and the equation for motion (in Cartesian coordinates -), respectively.
One of the major obstacles to understanding and predicting the behavior of flows, and thus simulating computations, lies in the random properties of turbulent flows [2]. To simulate a given turbulent flow, all turbulence scales must be discretized, resulting in the need for extremely sparse meshes and thus requiring many, sometimes prohibitively expensive, computational resources. This approach is known as DNS [3].
Other ways to simulate a flow with turbulent properties include mathematically modeling the smaller scales and solving for the larger turbulence scales (LES) or even fully modeling the turbulence scales, providing a way to circumvent the random properties of turbulence by relating them to average conditions and obtaining their time average after simulation (the technique is known as RANS, from Reynolds-Averaged Navier Stokes) [4].The use of turbulence models has its efficiency proved in several areas in engineering, having their results validated by experimental data. Some examples are: simulation in relief valves [5], totally turbulent flows in pipes [6], in natural circulation circuits of nuclear reactors [7] and in the prediction of the characteristics of aircraft [1].
Numerical simulations are an alternative or complement for testing geometries, especially in the aerospace field. The increasing memory and processing capacity of computers has enabled the use of more accurate codes and more complex models. Their great attractiveness lies in the reduction of the high cost of building models and their simulations in wind tunnels [8] and full-scale prototypes. However, it should be noted that numerical simulation cannot eliminate the need for practical testing, prototype construction, and physical studies, making it a complementary tool in the development of projects [1].
The present work aimed at numerical simulation of incompressible steady-state airflow over a NACA0015 hydrofoil with flap and Gurney flap using the standard k-epsilon and k-omega turbulence models. The meshes used for the simulation were structured meshes created using the ICEM CFD® module. The results were compared with experimental data obtained by [9], [10] and [11]. The commercial software ANSYS CFX® was used for the computational simulations.
Flows on airfoils are of great interest, mainly because of their numerous practical applications. They also exhibit special characteristics, such as abrupt changes in behavior with small changes in parameters, such as the angle of attack, with the possible sudden loss of lift and increase in drag force in a phenomenon known as stall.
After experimental investigations in search of a better understanding of stall [12], three stall configurations were registered: one in which a laminar separation bubble formed at the leading edge and later reconnected with the flow, another in which this bubble did not reconnect and caused a stall at the leading edge, and the stall at the trailing edge.
After its discovery, laminar detachment bubble was not paid much attention at first because it did not have a predominant influence in profiles with larger cross-section, although its presence in more obtuse profiles was known [13].
With the development and use of thinner airfoils, several studies were conducted to determine their effects on aerodynamics [14].
Eggert and Rumsey [15] performed simulations of a flow on a NACA 0018 aerofoil with a Reynolds number of . Their results showed that the shear stress transport model (SST) and the Spalart-Almaras model (AS) always respond well when the flow conditions are not transition conditions. In addition, better approximations to real flow were found when three-dimensional computational models were used instead of two-dimensional models and domains with the same dimensions of the wind tunnel in which the experimental data were collected. However, as reported, the experimental study was not able to ensure the minimization of the three-dimensional effects on the wind tunnel as well as the adequate treatment of the evolution of the boundary layer of the tunnel walls in the test section.
Balakumar [3] simulated the flow on a NACA 0012 airfoil using the DNS approach at Reynolds numbers from to and angles of attack from 5° to 15° and pointed out the formation of laminar separation bubbles in both cases. The dimensions of the laminar separation bubbles were related to the Reynolds numbers, indicating smaller bubbles at larger Reynolds numbers and vice versa. Static pressure peaks were found in the region where the bubbles formed.

II. Methodology 2.1. Geometry
For the proposed study, a geometry was created with the same dimensions as Lee and Su [9], which were strategically chosen to allow comparison of results. The design procedure started with the use of Javafoil® software and was completed with Autodesk Inventor® computer-aided design software. Since the latter is a three-dimensional CAD software, the profile was assigned an arbitrary thickness of 30 mm. Figure 1 shows the geometry of interest.

Discretization
Most of the errors in numerical analysis occur during the creation of the mesh and, in the case of external flows, also due to the insufficient determination of the size of the computational domain. In the simulations performed, the size of the computational domain and the distribution of elements in sufficient density and concentrated in regions of high property gradient were studied to obtain answers that are independent of the mesh.
ICEM CFD® software was used to build the domain and create the mesh, exporting geometry created with Inventor®.
Direct evidence for assuming the correct dimensions for the size of the computational domain for external flows has not been found by the author in the scientific literature. An approach for this determination is proposed.
The basis of the proposed methodology is based on an application involving an arbitrary airfoil with a span much larger than its chord, which allows the two-dimensional approximation of the problem. Other properties attributed to the flow are subsonic, incompressible, and adiabatic flow. Therefore, the Navier-Stokes equations begin to have elliptical mathematical properties, i.e., what happens downstream affects what happens upstream of the flow. The question, the answer to which was the starting point for the methodology, was: what are the distances from the blade surfaces at which the disturbances they cause are practically imperceptible? From this point on, the pressure and velocity gradients would be almost zero, since their values would then be constant. For the computational analysis, this means that the domain is dimensioned so that the velocities at a very small distance from the boundary conditions practically all have the same value, and that from this point on there is no more disturbance of the flow by the obstacle (the airfoil), so that an optimal domain is created. The parameter was used as a comparison value between the value measured at a given location ( ) and the distance value, where is a dimensionless quantity representing a deviation from the velocity of the free stream ( ).
A deviation was previously assumed to be acceptable, and an initial domain was established to serve as a starting point for successively increasing the dimensions of the computational domain. Figure 2 shows a simplified and general conceptual scheme of the variation of , the dimensions of the initial estimate, and the boundary conditions. Another hypothesis adopted for the methodology is that each turbulence model has its optimal domain, so the verification becomes necessary for each model. For the study, the k-omega [16] and k-epsilon [17] models were used, which employ the Boussinesq approach for the Reynolds stress tensor and define the turbulent viscosity by equations 3 and 4, the turbulent kinetic energy by equations 5 and 6, the specific dissipation rate by equations 7 and 8, and the Reynolds stress tensor by equation 9. (3) With being the turbulent dynamic viscosity, being the turbulent kinetic energy, being the dissipation by unit of turbulent kinetic energy, being the turbulent dissipation rate, being the Reynolds stress tensor and being the Kronecker delta. And the closure constants are for the k-omega model: ⁄ ⁄ ⁄ ⁄ ⁄ , and for the k-epsilon model: .
It is important to emphasize that these two-equation turbulence models were chosen because they are widely used in industry and because they are models that tend to converge more easily (stability and lower computational time). Of course, there are more complex and detailed turbulence models that could be more advantageous candidates for capturing transitional flows (e.g., separation bubbles), and we mention here the [18] or the physics-based k-kl-ω [19], but they are not the target of this study. The main idea is to minimize the errors by changing the constants of the standard turbulence models (k-epsilon and k-omega) and comparing the responses with the experimental data using the available and simple RANS models to obtain better average results.
ICEM CFD® was used to create the structured meshes and a computational domain with rectangular shape was chosen. Based on the first estimate, the velocities (along and at a very small distance from the boundary condition) were verified. Based on the values of V^*, the dimensions of the calculation area were completed. It was found that the distance between the boundary condition "ENTRANCE" and the profile had a deviation of about , so the original estimate was kept. For the other boundary conditions, it was found necessary to increase the original dimensions. For the distance between the wing and the boundary condition "EXIT", a larger deviation was observed in the velocity, so a fivefold increase in the chord was assigned for each change in the area. For the distance between the wing and the boundary condition "WALL", a twofold increase in the chord was assigned to the smallest variation. The values and properties assigned to the boundary conditions were: ENTRANCE velocity of 15.2 m/s, WALL free-slip condition, WING no-slip condition, EXIT relative pressure of 0 Pa, and SYMMETRY symmetry condition.
Other properties of the fluid for the simulation were: ideal gas behavior, incompressible and isothermal flow at a temperature of 25°C. The residual target for the simulations was for dimensionless quantities.       The software chosen for the analysis was ANSYS CFX®, which uses the finite volume method and performs only three-dimensional simulations, so the thickness of 30 mm generated with the software CAD Inventor® was maintained. The parallel symmetry condition shown in Figure 2 is a mathematical artifice to ensure that the simulation becomes two-dimensional.
The verification of the number of elements was based on the convergence of the values of the lift coefficient (CL) and the drag coefficient (CD) with the variation of the number of elements according to the methodology of [21]. The results shown in Figure 9 were obtained for successive refinements. The final mesh was selected, with approximately elements. Figure 10 presents the mesh and a detail in the trailing edge. For comparison with the experimental results of [9], curves of the pressure coefficient (CP) were plotted over the dimensionless chord for each model used ( Figure 11). Figure 11 shows the curves and the difference of the pressure coefficient near the leading edge of the numerical models and the experimental model. Figure 12 shows the results for the velocity field for the standard k-epsilon and k-omega models. In this figure, the differences in the regions marked A and B are highlighted. Figure 11. CP curves by dimensionless chord. In red color, the curve for the k-epsilon model; in green color, the one for the k-omega model; and in black color, the experimental [9]. Figure 12. Comparison between the color maps for the k-epsilon and k-omega models. The regions indicated by the -A‖ and -B‖ arrows present the main differences between the simulated flows.

Optimization of the Constants of the Turbulence Models
In the k-epsilon model, constants that presented greater variation, after a sensitivity analysis based on the chain rule ( and ) [22], were Cμ and Cε1. Based on the simulation with the default values, the ANSYS CFX® module, called Response Surface Optimization, was used, a software in which 40 exploration points were created. At these points, the closure constants were changed at the following intervals: and .

EC-2022-739
Experiments were designed using the optimal space-filling technique with maximum entropy (ANSYS, 2017) as the criterion. To check convergence, the software was configured to display the number of iterations needed to achieve convergence in addition to the CD and CL errors, with the goal of monitoring.
Using the results of these points, a response surface was created using the Kriging technique. Figures 13 and  14 show the response surfaces for the CL and CD errors, respectively. In these simulations, only two cases were found where the points did not converge.  Based on these plots, we looked for which combination of constants resulted in the best CL and CD values. The Pareto curve, shown in Figure 15, shows the different combinations of Cμ, Cε1 and their respective CD and CL error values. The software indicated three candidates with the best results, and Table 1 shows these values.
Since these points were determined based on the response curves obtained from the simulated points, it was necessary to simulate the selected condition to verify the error estimated by the optimization.  Solution 1 was selected from Table 1 because it yielded the smallest error for CL. Figure 16 shows the velocity field and the CP curve through the dimensionless chord for the selected solution. The calculated error values were 20.4 % for CD and 9.5% for CL. For the k-omega model, it was necessary to perform the sensitivity test of its constants. The results were, for CL: and for CD: .
The and were chosen for the optimization process because their influences were superior to the other constants. The interval chosen for the simulations was and .

EC-2022-741
The simulation results showed that only two points did not reach convergence. The response curves can be seen in Figures 17 and 18. Table 2 shows the three best values for and and their expected errors.   Various techniques and configurations were used for the values of the first solution, but convergence to an analysis on a permanent basis was not possible. When the transition regime option was used, it was observed (using the CD variables) that the flow did not exhibit the typical fluctuation caused by the release of vortices, but showed an almost random behavior that was exclusively numerical in origin and increased as the simulation progressed over time.
With the values of the second solution it was possible to achieve the convergence of the simulation with a permanent flow condition, but the errors of the coefficients CD and CL were 79.3% and 25.3%, respectively. For the last solution, the error values of the coefficients CD and CL were 26.6% and 53.7%, respectively.
Since there was no significant deviation in the drag value, an investigation of the stability of the constants was carried out, looking for limits at which convergence was not possible, in order to determine a safety interval for the constants. For this purpose, a separate analysis was performed for each constant, with a maximum deviation of ±50% of the original values. The result is shown in figures 19 to 24.     The results showed that the value of the -coefficient above 0.1035 caused the nonconvergence of the simulations. For the -closure coefficient, an interval was determined in which the change in the coefficients represented some change in the results; this interval is considered between 0.06825 and 0.08175.
Thus, for the k-omega model, the method of changing the constants did not have as positive effect as for the k-epsilon model, and for this reason the results for this model (k-omega) are not discussed in the verification/validation section.

Verification/Validation
In order to verify that the modifications to the k-epsilon model can produce responses with better quality, regardless of the geometry, the simulation of the NACA 0012 profile was performed, since it is a profile for which a lot of data are available in the literature, having been studied in detail in dozens of wind tunnels over a period of more than 50 years [20]. The experimental data were taken from the publication of Gregory [10] with free transition and Ladson [11] with fixed transition. An angle of attack of 10° and a Reynolds number of were defined in a model with 760 mm chord and air properties at 25°C for an undisturbed flow velocity of 55 m/s. The domain verification and the study of the mesh followed the methodology already presented. Figure 25 shows the velocity deviations at the extremities of the domain and the convergence curve using the coefficients CL and CD as parameters for monitoring according to the criterion of Stern and Wilson [21].  Figure 26 shows the comparison of the experimental data obtained by [10] and [11] with the k-epsilon and the optimized k-epsilon model, where the maximum lag obtained is 0.0001.
The comparison showed that the optimized model approximates the experimental results obtained from tests on airfoils where a device was used to control the transition from the laminar to the turbulent state (this feature was not present in the boundary conditions). Therefore, in the optimized model, either laminar separation bubbles did not form (this behavior occurs because the turbulence model considers the flow to be fully turbulent) in the leading edge region or the pressure fields were deformed to produce a CP curve through the dimensionless chord similar to the experiments. On the other hand, the standard K-epsilon model represented the pressure distribution curve (as shown in Figure 26) for a NACA0012 profile with a clean surface with reasonable accuracy in agreement with experiments. Figure 26. On the left, the velocity profile on the extremities of the domain, and on the right, the convergence study.

V. Conclusions
Numerical fluid dynamics was applied to an external, two-dimensional, incompressible subsonic airflow on an airfoil (NACA 0015 with flap and Gurney flap) in a structured mesh using k-epsilon and k-omega turbulence models.
The results of the simulations show that the k-epsilon model has better stability when its constants are changed. The changes to this model gave better results when the results of the numerical simulations were compared with the experimental data from [9]. However, when the changes were verified with another known geometry (NACA 0012), the trend was not sustained.
As expected, the separation bubble was not observed, but the behavior of the pressure distribution was determined with reasonable accuracy by the two-equation K-epsilon model.
Further studies should be conducted to determine the results of the turbulence models with the modified constants in other geometries. The modifications with the method proposed by Oguma [21] had no effect when applied to the k-omega model.