Comparing different CFD software with NACA 2412 airfoil

The engineering application’s design process starts with a concept, based on our theoretical knowledge and continues with a numerical simulation. In our paper, we review the finite volume method (FVM) which is used generally for heat and fluid dynamic simulations. We compare three different computational fluid dynamics (CFD) software (based in the fine volume method) for validating a NACA airfoil, which can be used for example in the aerospace industry for an airplane’s wing profile, and it can be used for example in the renewable industry for a wind turbine’s blade or a water turbine’s impeller profile. At the end of this paper, the result of our simulations will be compared with a validation case and the difference between the CFD software and the measured data will be presented.


INTRODUCTION
In agricultural areas which are far from the electric grid, establishing a connection with the energy system is expensive and if the investment returns, it will be in a long time. For this region, an independent, off-grid renewable energy source is a good option, i.e., the wind and hydropower. The development process of this power generation equipment is the same as with other products. It starts with a preconcept which is based on our theoretical background or a previous similar product. After the preconcept, the design engineer creates a model which is tested with simulation products.
The simulation method can be different depending on the simulated physics, i.e., for stress simulation the engineer usually uses a finite element software, for fluid simulation they use a finite volume (FV) software. Besides, there are the meshless, finite difference, domain decomposition and other methods. With the simulation tools, the engineers can iteratively improve product performance and eliminate the possible errors.
For fluid dynamics collectively we called the simulation software computational fluid dynamics (CFD). In our article, we will use three different CFD software, which are based on the finite volume method (FVM). We will validate an airfoil (which is also called wing profile) with these tools. The airfoils are the principal parts of the modern lift force-based equipment, i.e., the previously indicated wind turbines' blade and water turbines' impeller. We will reflect the difference with this validation process within the same method with different software.

Briefly about airfoils
Nowadays the geometry of wind turbines blades is based on wing profiles used in the aerospace industry and the impellers use the lift force favorably. For this process, the turbines' blade is made of airfoils, which are known as wing profiles in the aerospace industry. If the working medium is water, the airfoil's name changes to hydrofoil referring to the liquid medium.
Perhaps one of the most well-known airfoil families is the National Advisory Committee for Aeronautics (NACAs) wing profiles, for which a large amount of measurement results are available. The NACA profiles have a four-or five-digit code, which represents the geometry. In the four-digit system, 1. The first digit shows the maximum asymmetry (called chamber) between the two acting surfaces of an airfoil in the percentage of the cord line length. 2. The second digit shows where the point defined with the first digit can be found on the cord line. Also, this point is where the cord line and the chamber mean line (green line on the next figure) have the greatest distance. 3. The third and fourth show the percentage of the length and thickness (Abbott & von Doenhoff, 1959). Fig. 1 presents a NACA 2412, which we chose for simulation. If we put the airfoil in the position seen in the next figures (Figs. 2 and 3) and the direction of the free flow is horizontal, the pressure on the upper part of the airfoil decreases, and the velocity increases. In the lower part, the pressure increases while wind velocity drops. The pressure difference between the sides generates the lift. The lift force can be calculated with the following equation: In the equation F L is the lift force, c L is the lift coefficient, r is the density, A is the area of the blade's cross-section, v 2 ∞ is the velocity in free stream.

Shortly about numerical fluid dynamics
On the CFD software's market there are a large number of CFD software with a different mathematical background. We used three CFD software (CFX from ANSYS, OpenFOAM from OpenFOAM Foundation, and FloEFD from Mentor Graphics) based on FVMs.
The FVM based software divides the computational domain into FVs. Pressure, velocity, and temperature can be defined with the help of the conservation laws (mass, momentum, and energy). This calculation typically is based on the following transport equation: v vt In the equation the vU=vt is the time-dependent member (vU=vt ¼ 0 the steady-state), U is a conserved quantity (e.g., mass, impulse, energy), F is the same quantity's flux, S v is the U quantity's volumetric source, S A is the U quantity's surface source, V is the arbitrary enclosed control volume, A is the surface of this control volume and R is the residual (error of the equation).
For the computational domain's discretized part (each mesh cells), the solver starts to compute with an initial value and iteratively recalculates the flow field. In this calculation, the software analyses the residual quantity (R), which should decrease and go towards 0 or a prescribed zero close value. After the computation, the solver calculates each physical value, in Reynold-Averaged Navier-Stokes method (RANS), the main quantities are the velocity, pressure, and temperature.

Turbulence
In the engineering application, the flow characteristics are usually turbulence phenomena. The turbulence is a transient, random, chaotic, diffusive, and spatial (3D) phenomenon. It can occur due to near the wall's boundary layer, in the mixing flow shear layer and the unstable density stratification. The development of a turbulence flow is shown in Fig. 4.
In the CFD software, the turbulence can be modeled in several ways. The most computational resource-intensive method is not to use any turbulence model and wall function (which describe the flow in the near-wall region). In this case, the software calculates the viscous sublayer, the buffer layer and the turbulent region, where the mesh size should be smaller than the smallest eddy that influences the flow field and its properties.
There is a less computational resource method, which is to use a turbulent model. The most common RANS transport equation-based turbulence models are the two-equation models, for example, the k-« and k-u models. Each model has advantages and disadvantages, therefore Menter in 1993 created the SST (Shear Stress Transport) model, which includes the benefits of both models.
When we use a wall function to describe the flow field in the near-wall regions, the computational requirements are much less than in the previous case.

k-« turbulence model
The k-« turbulence model explains the flow field variability with two equations. The first equation is the k, which describes the kinetic energy, the second is the «, for the turbulence dissipation energy.
The benefit of the model is the need for low computing resources, and it was one of the first turbulence models, so it was validated with a lot of experiments and its limits are well known.
The disadvantage of the model is that it shows poor performance for fully developed flows in noncircular pipes, for rotating flows, for flow with large strains, and for mixing layers (Jakub ık & Feszty, 2016; Krist of, 2019).
The most common k-« models are the Standard, the Realisable, and the RNG versions. We used the Standard k-« model for our simulation.

k-u turbulence model
The k-u model is a two-equation eddy-viscosity turbulence model like the k-«. In this case, instead of «, the solver computes u, which is the specific rate of dissipation.
The model shows better performance near the walls, but its results highly depend on the turbulence intensity of the boundary conditions (BCs), which is a disadvantage compared to the k-« (Jakub ık & Feszty, 2016).

SST turbulence model
The SST turbulence model is a two-equation eddy-viscosity turbulence model, which used k-u near the walls and k-« for free stream. The model used the benefits of each model, therefore it describes the fluid field more accurately than the other two-equation models. Its known disadvantage is the overestimation of the turbulence intensity in the large strain regions (e.g., at the stagnation point which is shown in Fig. 5), but it is less like the overestimation of the k-« (CFD Online).

SIMULATIONS PROPERTIES Computational domain
We used a 60 3 61 3 1 m domain with plain flow (2D) simplification for the simulation. Fig. 6 shows the whole domain with the airfoil.

Mesh for Mentor Graphics FloEFD
The FloEFD uses its own Cartesian mesh process, which generates structured hexahedron elements, whose edges are aligned with the global X, Y, and Z coordinate axes.
The mesher is using the cell mating technology for finer mesh size and it uses the cut-cell method for the elements near the boundary of the walls. The cell mating process finer the grid by splitting the cells into the half, in 3D it generates 8 cells from 1 cell (halved by the X, Y, and Z axis). The cut cell method divides the mesh elements near-wall in two pieces and generates polyhedral elements on the fluid and the solid sides, using Mentor Graphics' SmartCells technology to describe the flow field for this control volumes (Mentor Graphics, 2019). Fig. 7 shows these methods (the colors from the blue to red represent the cell mating's steps).  The FloEFDs mesh grid is in each case (3D and 2D) a volumetric mesh, so in our study, the initial mesh (represented as blue in the previous figure) thickness has only one element.
The simulation may give different results depending on the mesh density, therefore for 08 angle of attack (AoA) we ran mesh dependency simulations, where we followed the convergence of the lift coefficient (shown in Fig. 8).
As a result of the analysis, we chose the 1,520,510 number element settings. The final mesh can be found in the next three figures (Figs. 9, 10 and 11).

Mesh for ANSYS CFX and OpenFOAM
For the CFX and OpenFOAM, we were able to use the same mesh. The mesh can be converted from the ANSYS products to the Cfx4ToFoam, fluentMeshToFoam and fluent3DMeshToFoam commands to OpenFOAM, and with the FoamMeshToFluent we convert the grid from OpenFOAM to ANSYS Fluent, which can be used for CFX. For this grid, we used an unstructured, deformable surface followed mesh, which contained hexahedron and prism elements. Each software uses a volumetric mesh for 3D and 2D simulations, so the final mesh was a one element thick mesh.
Like previously, we started the validation process with a mesh dependency analysis for the AoA 5 08 case, whose results are shown in Fig. 12. The lift coefficient started converging very quickly, and in the last five cases, the difference was below 1% compared to the previous simulation. We did not see this fast convergence in the case of the drag coefficient, therefore we continue to finer the cells near the airfoil's region. For the grid which contain 152,578 elements the drag coefficient change was decreased to 1.7% (compared with the previous grid), and its value does not change significantly for the last grid, which contains 377,644 cells. As a result of these simulations, we chose the 152,578 elements mesh, which can be seen in the next three figures (Figs. 13, 14 and 15). In the three different software, the open condition was prescribed in two different ways, in FloEFD, we defined the absolute environment pressure (101,325 Pa), for the other software we defined a relative pressure (0 Pa). In the FloEFD and CFX, the symmetry was made with the symmetric condition. For the OpenFOAM, in 2D we have to use the "empty" BC instead of symmetry, which is in this CFD software a 3D BC.
The inlet velocity was 42.89 m/s. We changed the inlet vector between 0 and 108. The Reynold number was RE 5 2.85310 6 .
In FloEFD and CFX, we used the "Air" fluid from the software material database. In the OpenFOAM, there is no material database, so we defined the material with the following properties: r 5 1.2041 kg/m 3 , v 5 1.511083 3 10 À5 m 2 /s.

Stopping criteria
The software solves Eq. (2) for each mesh elements and monitors the U quantity. Based on the continuity laws, the entering flow should leave the flow fields and no residual fluid should remain in the system. The FVM based CFD software iteratively monitor this residual U value and the calculations run again and again until the residual value decreases near the required 0 value.
For each CFD software, we use the default quantities (velocity [u, v, and w], pressure, turbulence [k and « or u] and the temperature) and the lift and drag coefficients and forces. For the lift coefficient, we used the following equation (which is the (1) equation rearranged).
In Eq. (3) c L is the lift coefficient, F L is the lift force, r is the density, A is the area of the wing's cross-section (which is, in this case, the cord length multiplied by the wings' thickness (Guly as, 2016)), v ∞ is the air velocity in the free-steam. For the AoA 5 68, the changing of the lift coefficient is shown in the next figure (Fig. 16) during the first iteration to the last iteration.
Having compared the simulations, we found that the change of the c L s value decreased below 10 À5 between 400th and 500th iterations, below 10 À6 between 550th and 600th iterations and below 10 À7 between 700th and 750th iterations. Therefore, we set the maximum iteration to 1,100.

Comparing steady-state with transient simulation
Validating the NACA 2412 airfoil is a time-dependent process, so it should be modeled with a time-dependent simulation. We ran a steady-state and a transient simulation for 0 and 108 AoA. The results of the simulations are shown in the next figures (Figs. 17,18,19 and 20) and in Table 1.
In Table 1, F x is the force on the wing's surface (X component), F y is the force on the wing's surface (Y component), F L is the lift force and c L is the lift coefficient.
Comparing the flow field and the results in Table 1 (especially the lift coefficient), the differences between the steady-state and transient results are negligible. The consequence of running the simulation in steady-state is that we are unable to examine the turbulent phenomenon's time dependency, only the average value of the eddies can be examined.

SIMULATIONS' RESULTS
We ran our simulations between 08 and 108 with 28 steps. The next figures  show the AoA 5 08 and the AoA 5 88 the cases' velocity field. The lift coefficients can be found in Fig. 25. Examining the diagram, the lift coefficients are approximately the same at 0-68, above the CFXs, and OpenFOAMs results, remain the same, while the Mentor Graphics FloEFDs c L values were a bit higher than the others.

COMPARISON OF RESULTS
We would like to compare our results with the measured data. We chose three different measurements and we made one averaged value about them. These values are shown in Fig. 26. The three c L series were measured with different free stream velocities, but they could be compared without any correction because the lift coefficient was independent of Reynold's  (Abbott & von Doenhoff, 1959;Matsson et al., 2016;Seetharam et al., 1977) Fig. 25. Lift coefficients versus AoA number (Abbott & von Doenhoff, 1959). We can conclude that the measured value is similar to Abbot's and Seetharam's measurements. There was only a significant difference over 108.
Matsson's results below 108 (except for AoA 5 88) are similar, but the values are less than in the two other cases. In Fig. 27, we can see the comparison of the simulations' results and the average of the measured lift coefficients.
In Fig. 28, we can see the difference between the simulations' results and the measured average lift coefficients.
In the two previous figures, the lift coefficient deviated with a couple of percentage in the 0-68 AoA range for the CFX and the OpenFOAM, while for the FloEFD the difference in this range was between 2.47 and 7.44%. In the 8 and 108 for each software and each turbulence model, the deviation increased up to 11.17%.
In the case of FloEFD, the lift value was higher, for OpenFOAM with SST, the lift value was always lower in each case than the averaged measured values. In the case of other software and turbulence models, we could not find any regularities.

CONCLUSIONS
In this article, we reviewed how to calculate lift coefficients for wing profiles. Then, we briefly described the FVM and our chosen software's meshers and then we presented the results of our simulations.
Our simulation was about to run validation for a NACA 2412 profile in the range of 0-108. Comparing the results in steady-state (stationery) and in transient (time-dependent, unsteady state) for AoA 5 08 and AoA 5 108, the difference was negligible, therefore we ran our simulation in steady-state.
With our simulation, we established that with the chosen method and software in the case of the lower AoA, the lifting coefficient can be simulated with lower differences, in the case of the larger AoAs (in this case 88 and above) transient simulation is necessary for accurate results.