TURBULENCE MODELS APPLIED TO CLASSICAL FLUID MECHANICS AND HEAT TRANSFER PROBLEMS: THE PERFORMANCE EVALUATION OF THE OPEN SOURCE CFD PACKAGE OPENFOAM MODELOS DE TURBULÊNCIA APLICADOS A PROBLEMAS DE MECÂNICA DOS FLUIDOS E TRANSFERÊNCIA DE CALOR: AVALIAÇÃO DO DESEMPENHO DO PACOTE DE DFC DE CÓDIGO ABERTO OPENFOAM

Topics related to the modeling of turbulent flow feature significant relevance in several areas, especially in engineering, since the vast majority of flows present in the design of devices and systems are characterized to be turbulent. A vastly applied tool for the analysis of such flows is the use of numerical simulations based on turbulence models. Thus, this work aims to evaluate the performance of several turbulence models when applied to classic problems of fluid mechanics and heat transfer, already extensively validated by empirical procedures. The OpenFOAM opensource software was used, being highly suitable for obtaining numerical solutions to problems of fluid mechanics involving complex geometries. The problems for the evaluation of turbulence models selected were: twodimensional cavity, Pitz-Daily, air flow over an airfoil, air flow over the Ahmed blunt body and the problem of natural convection between parallel plates. The solution to such problems was achieved by utilizing several Reynolds Averaged Equations (RANS) turbulence models, namely: k-ε, k-ω, LamBremhorst k-ε, k-ω SST, Lien-Leschziner k-ε, Spalart-Allmaras, Launder-Sharma k-ε, renormalization group (RNG) k-ε. The results obtained were compared to those found in the literature which were empirically obtained, thus allowing the assessment of the strengths and weaknesses of the turbulence modeling applied in each problem.


INTRODUCTION
The use of numerical simulations as an integral part of some engineering systems projects, such as in the power generation, automotive, aerospace etc. is already a reality that brings several benefits, such as reducing the cost and time for product development. The increasing use of computational procedures is due to the advancement of the processors' performance as well as the numerical simulation software and the decreasing costs of hardware equipment (ROCHA; SILVEIRA, 2012;WEHMANN, et al., 2018). Given that most of the flows, which are observed in the practical context of designing devices and systems, are characterized by being turbulent (VERSTEEG; MALALASEKERA, 2007), the study on the modeling of such problems is of utmost importance in the training of professionals that will work with this subject, which occurs mainly in technological and exact sciences. In many cases, the evaluation of turbulent flows is made by the application of turbulence models, which are solved by the use of appropriate numerical methods. In this sense, the evaluation of the performance of turbulence models in the solution of problems commonly encountered in the engineering context becomes very important.
Numerical techniques for modeling and simulation are an efficient tool to analyze turbulent flows. Direct Numerical Simulation (DNS) of Navier-Stokes equations can solve the calculation of incompressible turbulent flows, being limited to low and moderate Reynolds numbers, because it requires a large computational effort (KUROKAWA; CORRÊA; DE QUEIROZ, 2018). Unsteady Reynolds Average Navier-Stokes (URANS) equations with high Reynolds k-ε turbulence model may be a useful alternative to Direct Numerical Simulation (DNS) in computational simulations (OZMEN-CAGATAY;KOCAMAN, 2010;QUECEDO, et al., 2005), however its modeling performance needs an efficient numerical method to result in an oscillation-free solution without introducing excessive artificial dissipation (KUROKAWA; CORRÊA; DE QUEIROZ, 2018).
The turbulence model chosen to run the simulation usually has a direct influence in the results. Each turbulent model can be appropriate to specific computational problems, reaching better results to specific cases. The most know turbulence models are the k-ε and k-ω, being successfully applied for industrial applications. The k-ε model is based in two equations which are used in the evaluation of the turbulent viscosity. This viscosity is mathematically modeled as a function of k, the turbulent kinetic energy, and ε, the viscous dissipation rate of turbulent kinetic energy. This model is recommended for fully turbulent flows. The k-ω model is also a two equations model, and it is most v.1, n.5, 2021 recommended for near-wall regions, achieving better results in this kind of problem. Beyond the boundary layer region, the ω equation behaves largely sensitive, so an upgrade model called Shear Stress Transport (SST) k-ω was developed. This model is a combination/transition between the k-ω model in near-wall regions and the k-ε model in regions far from the boundary layer. This approach is known to be successful in cases with boundary layer separation (JONES; LAUNDER, 1972;MAZARBHUIYA;BISWAS;SHARMA, 2020;MENTER;KUNTZ;LANGTRY, 2003).
Over the years, turbulence models have been applied by many researchers to solve complex engineering problems that simulate real-world applications in CFD softwares. Khavaran (1999) consolidated the use of an acoustic analogy with a compressible Reynols-Avaraged Navier-Stokes (RANS) solution and k-ε turbulence model to predict jet noise. Menter, Kuntz and Langtry (2003) described the SST turbulence model, as well as this model enhancement, modifying near wall treatment of the equations. Quecedo et al., (2005) described and compared two mathematical methods for solving a dam-break problem using the finite element method, one by solving Navier-Stokes equations and other by solving Shallow Water equations. Considering the computational efforts required, the results showed that the Shallow Water approach should be used for computations involving large domains, and Navier-Stokes approach for fine calculations in which knowledge of the 3D structure of the flow is required. Similarly, Ozmen-Cagatay and Kocaman (2010) also got results to dam-break flows, comparing experimental and numerical results using a commercially available CFD program, solving RANS equations with k-ε turbulence model involving the Shallow Water equations. Sacomano Filho, Fukumasu and Krieger (2013) applied numerical simulations based on RANS k-ε turbulence model within the stochastic separated flow formulation to analyze the physics of ethanol turbulent spray flames . Dos Santos et al. (2014) compared the use of Large Eddy Simulation (LES) and RANS with k-ε turbulent model in combined convective and radiative heat transfer for nonreactive turbulent channel flows. Yahiaoui et al. (2016) developed a numerical study by the commercial CFD code, ANSYS Fluent, using a finite volume method to solve the steady-state RANS equations, and the turbulence models Spalart-Allmaras, realizable k-ε and k-ω SST were used to produce a closed system of solvable equations. Rosa et al. (2017) compared four models to predict jet mixing noise, two derived from the Lighthill Acoustic Analogy (LAA) and other two from the Linearized Euler Equations (LEE) with added source terms. RANS equations, DNS and LES were used to the computational solution of the mean flow, and empirical functions were used to model turbulent correlations. Kurokawa, Corrêa and de Queiroz (2018) analyzed the effectiveness of a numerical technique based on the finite-difference GENSMAC methodology coupled with high Reynolds k-ε turbulence model and the ADBQUICKEST upwind scheme to deal with the advective terms for a dam-break flow and a turbulent jet impinging orthogonally onto a flat surface. Coimbra and da Silva (2020) used RANS equations to solve the turbulent closure problem. The impact of mesh refinement levels and boundary conditions on the swirl number and overall flow structure was investigated. The realizable  In front of the presented, our work aims to apply and test the turbulence models that use the Reynolds Averaged Equations (RANS), namely: k-ε, k-ω, Lam-Bremhorst k-ε, k-ω SST, Lien-Leschziner k-ε, Spalart-Allmaras, Launder-Sharma k-ε and renormalization group (RNG) k-ε; with the subsequent comparison of the numerical results obtained against those empirically known for the classic problems under study: two-dimensional cavity, Pitz-Daily, air flow over an aerodynamic profile, air flow over the Ahmed body and the natural convection problem between parallel plates. The OpenFOAM opensource package was used in obtaining the numerical solutions, being highly efficient when working on problems involving complex geometries. The proposed approach allows the evaluation of the strengths and weaknesses of the turbulence models performance for each problem.

METHODOLOGY
To perform the numerical simulations, OpenFOAM opensource software was used, which In this work, the approach used to assess the problems was as follows: presentation of each problem and the information needed to solve them (still in the methodology section), comparison of the results obtained by applying each turbulence model for each problem against those found in the selected literature (in the results section).

Two-dimensional cavity
The first problem approached is that of the two-dimensional cavity consisting of an isothermal incompressible flow of a newtonian fluid, in a two-dimensional domain with the shape of a square, whose upper face presents a movement to the right, as illustrated in Figure 1. The outline and initial conditions applied are as follows:

REVISTA CIENTÍFICA ACERTTE ISSN 2763-8928 MODELOS DE TURBULÊNCIA APLICADOS À PROBLEMAS DE MECÂNICA DOS FLUIDOS E TRANSFERÊNCIA DE CALOR: AVALIAÇÃO DO DESEMPENHO DO PACOTE DE DFC DE CÓDIGO ABERTO OPENFOAM Paulo Alexandre Costa Rocha, Felipe Pinto Marinho, Victor Oliveira Santos Stéphano Praxedes Mendonça, Maria Eugênia Vieira da
-It is assumed that the initial pressure field is uniform with an intensity of 0 m²/s²; -The velocity field is assumed to be uniform with magnitude 0 m/s; -The configuration of a movable upper contour (solid wall with constant velocity) and the remaining boundaries being stationary (solid walls) is applied; -In the stationary boundaries, zero gradient conditions are used for pressure, and non-slip with null velocity; -The upper wall received the non-slip boundary condition, with a fixed velocity equal to 1 m/s on the x-direction.

Transport properties
The only transport property required for this problem is the kinematic viscosity of the fluid, ν, where it was used the hypothesis that it is equivalent to 10-5 m²/s, so it is possible to determine the value of the Reynolds number for the flow, Re, by Equation (1).
Obtaining then a Re = 10000, which characterizes a fully turbulent regime for internal flows.

Mesh generation
Initially, a uniform mesh of 10 x 10 elements (10 mm x 10 mm) was used. Then the mesh went through a gradual refinement until the variation of the results was small enough for the numerical solution to stabilize. The mesh generator tool provided by OpenFOAM, blockMesh, was used for the generation of the uniform mesh. The mesh used to obtain the results was 100 x 100 elements (1 mm x 1 mm), being the most refined, as well as the one that provided the most consistent results with the literature. The mesh used is illustrated in Figure 2.  In a more detailed way, the following boundary and initial conditions were adjusted as: ▪ The non-slip condition on the walls was considered; ▪ A uniform velocity field was defined at the inlet. As carried out in the Pitz-Daily experiment, the cases for which the entry velocity were 9.1 m/s were considered for the case where Re = 1.5 x 10 4 , 13.3 m/s and 22.2 m/s for Reynolds values of 2.2 x 10 4 and 3.7 x 10 4 respectively; ▪ A null pressure value (manometric) was considered at the outlet; ▪ For the initial value of the turbulent kinetic energy, k, a value of 5% of the velocity was considered. Thus, for the case in which U = 13.3 m/s, for example, k = 0.663 J by Equation (2) and by the hypothesis that the components of the random fluctuations of velocity U'x, U'y, U'z, are equivalent:  Table 1 illustrates the initial values of the turbulent flow properties used in the models applied to the problem. The method of calculating these is well detailed by Versteeg and Malalasekera (2007).

Mesh generation
The check of the mesh quality was performed using different degrees of refinement and evaluating the change in the results for the velocity fields along the sections. The results considered in this comparison were obtained after a sufficient period for the flow simulation to be steady state.  able to produce results with a considerable degree of precision, making the use of a more refined mesh not necessary for the purposes of this work.   and perpendicular to the flow stream. A more detailed study on aerodynamic profiles is presented by Abbott and Von Doenhoff (1959) and Anderson Jr (2017).
With that in mind, the problem consisted in determining the Cl and Cd curves for the chosen profile, analyzing the influence of the Reynolds number and the angle of attack. The hypotheses applied to the flow are that it was incompressible, turbulent, two-dimensional, viscous and stationary.
The applied boundary and initial conditions in the simulation are: ▪ Uniform pressure field with zero gradient on the walls; ▪ For the flow velocity, a variable field was defined according to the direction of the considered angle of attack and the speed module defined for each Reynolds number. In addition, the non-slip condition was defined on the walls; ▪ The Reynolds number values considered are 2 x 10 5 , 5 x 10 5 and 1 x 10 6 ; ▪ The fluid properties were defined as air at atmospheric pressure and at 0 ºC temperature.
Furthermore, a density of 1.2928 kg/m³ and kinematic viscosity of 1.3304 x 10 -5 m²/s were used.
Being μ the dynamic viscosity of the fluid and l the characteristic length of the airfoil, in this case the value of the chord length.

Air flow over the Ahmed body
The Ahmed body is a classic and extremely important problem for the automobile industry.
Because of its geometry, proposed by Ahmed, Ramm and Faltin (1984), it generates the essential characteristics of the flow field of a real vehicle, with the exception of more specific effects, as seen on the lower part of the vehicle and surface protrusions such as rear-view mirrors. Figure 6 illustrates the dimensions of the geometry used in the simulation. The dimensions are in mm and the angle adopted
Therefore, the problem of air flow over the Ahmed body consisted in determining its drag coefficients. The hypotheses used are of incompressible, turbulent, three-dimensional, viscous and stationary flow.
For the simulations, the boundary and initial conditions are: ▪ The non-slip condition over the walls is considered; ▪ A uniform velocity field was defined at the inlet with an intensity of 30 m/s; ▪ A null pressure value was considered at the outlet and null pressure gradient at the inlet and over the contact regions with solid surfaces, such as the ground floor and the body; ▪ The reference area for calculating the nondimensional coefficients was assumed to be equivalent to the cross-sectional area and the reference length as the length of the body, in the case: Aref = 0.112032 m², lref = 1.044 m.

Mesh generation
To choose the level of mesh refinement, six different meshes were analyzed, whose results obtained are represented in Table 3. According to Table 3, it is perceived that the influence on the results due to the mesh refinement becomes small even for a significant increase in the number of elements, starting from the mesh with 13500 elements. Because of that, this mesh was chosen for the resolution and analysis of the problem.

Natural convection between parallel plates
This problem aims to study natural convection in a closed space. The geometry of interest consists of a cavity with a face of higher temperature, and another of lower temperature, separated by  Table 4.
In addition, the hypotheses applied to the flow are that it is three-dimensional, steady-state, turbulent, compressible and viscous regime. The boundary and initial conditions used in the simulation are: ▪ The velocity field was initially considered null and the non-slip condition was applied over the 6 faces. ▪ The pressure was defined as being uniform of value 1 x 10 5 Pa.

Mesh generation
The k -ω SST method was used to select the mesh and the effect of the refinement on the temperature values at an intermediate point (y = 0.5 m) was evaluated, and this effect can be seen in Figure 8. The different refinement levels applied are illustrated in Table 5.   Figure 8, it is clear that from mesh 2 on, the effect on the temperature is no longer significant. Thus, this level of refinement was chosen for the other simulations of this problem.

Two-dimensional cavity
The results for the problem of the two-dimensional cavity are represented in Figures 9 and 10, where the vertical component of the velocity is evaluated along the horizontal center axis of the cavity, while the horizontal velocity component is evaluated along the vertical center axis. The numerical results obtained show a behavior very close to those found experimentally in the work of Ghia, Ghia and Shin (1982), as can be observed by Figures 9 and 10. However, this similarity occurred only when the most refined mesh was used, 100 x 100. Thus, it can be affirmed that all turbulence models presented a good performance in the numerical prediction of the flow field for this problem. This fact is expected since the problem of the two-dimensional cavity is considered a standard case to test the accuracy and efficiency of numerical algorithms. In this case, since all models used are considered classic and have already been extensively validated, it is expected that they presented a good behavior for this problem. For this computational domain, the greatest deviations occur at the corners of the cavity where there is a tendency to the formation of recirculating zones, mainly affecting the performance of the models k -ε, Lam-Bremhorst k -ε, Launder-Sharma k -ε, Lien-Leschziner k -ε. The best fit is achieved by k -ω and k -ω SST models. In the first, the great advantage is that its integration to the solid boundaries does not require wall functions, since it is a low Reynolds number model. Nevertheless, there is a problem related to the specification of the boundary conditions for the turbulent frequency, ω, in the free stream, because the same tends to zero causing the dynamic viscosity turbulent, μt, tend to infinity. Because of this, small non-null values of ω must be specified in these regions. The k-ω SST model can adjust its behavior according to the region to be analyzed, since near the walls it presents a similar behavior to the k-ω model, while in the free stream it resembles the k -ε model, presenting the main advantages of each model. A better description of these characteristics is presented by Versteeg and Malalasekera (2007), Wilcox (1988) and Menter (1994).

Pitz-Daily
As already presented, one of the main goals of this problem is the determination of the reattachment point of the stream after the step. Thus, Table 6 shows the percentage errors, differences between the numerical results obtained and those found by Pitz and Daily where they are divided by the empirical results. These errors were found considering all the velocities addressed in the problem, 9.1, 13.3 and 22.2 m/s. It is noteworthy that when the Spalart-Allmaras model was applied, no numerical convergence could be achieved. Lam-Bremhorst k -ε 0.612% -3.233% Launder-Sharma k -ε 5.544% -8.800% Renormalization Group k -ε 4.892% -9.806% The analysis of Table 6 shows that, in general, all turbulence models show good performance, highlighting the k-ε and Lam-Bremhorst k-ε models. The former is widely used in industrial applications, because it presents good performance in areas of free stream with small pressure gradients and for being extensively well validated. However, its performance becomes poor in regions close to walls where there is a decrease in the Reynolds number and in situations in which there are high gradients of adverse pressure, as well as in some unconfined flows. The Lam-Bremhorst k-ε (LB), Launder-Sharma k-ε (LS) and Lien-Leschziner k-ε (LL) models are modified versions of the k-ε model for cases where the flow presents regions of low Reynolds numbers. The difference between the first two is in the boundary conditions of the rate of dissipation of turbulent energy, ε. This difference is well explained by Garg (1998). The latter was developed to replicate the behavior of the one equation Wolfshtein, Norris and Reynolds model, as well as providing the correct asymptotic response in the vicinity of the wall, as illustrated by Leschziner (2015). The Spalart-Allmaras model was designed for applications of external aerodynamic flows, in addition to presenting a good performance in boundary layers with adverse pressure gradients, which are important for the prediction of the stall phenomenon (VERSTEEG; MALALASEKERA, 2007). In general, the model is inappropriate for many internal flows, which may justify its numerical non-convergence.  Figure 11 shows the velocity fields plotted for a plane located at x = 0.1347 m for the case where Re = 2.2 x 10 4 , thus illustrating, by the similarity of the velocity profiles in the x direction, that the model applied does not directly impact the behavior of the field in this direction. Figure 11: Velocity field plotted for a plane located at x = 0.1347 m for the case Re = 2.2 x 10 4

Air flow over an aerodynamic profile
The values used as reference were provided by the database of the Research Group on Applied Aerodynamics of the University of Ilinois (UIUC AIRFOIL DATA SITE, 2019), and are presented in Figures 12 and 13. The green curve corresponds to Re = 2 x 10 5 , purple to Re = 5 x 10 5 and orange to Re = 1 x 10 6 . As it can be seen, Cl lies between 1 and 1.3 and Cd ranges from 0.012 to 0.022. It is worth to note that the model which presented the best performance was Spalart-Allmaras, and the results obtained by it are shown in Figure 14. This good concordance of the prediction provided by the Spalart-Allmaras model was already expected, since, as already discussed in this work, this model, which is characterized by being of a transport equation over one auxiliary variable that is related to the turbulent viscosity, is designed for applications involving external aerodynamic flows. For a better detailing on the mathematical description, one can recommend the references Versteeg and Malalasekera (2007) and Spalart and Allmaras (1992).

Air flow around the Ahmed body
The values of the drag coefficient for all the models applied and their empirical values are shown in Table 7. the application of the renormalization group theory, mathematical formalism from the statistical mechanics to the standard k -ε model. The main advantage is related to the presence of a deformation-dependent correction term, since one of the main limitations of the k -ε model is in the transport equation for the rate of dissipation of turbulent energy for flows that experience high deformation rates. For more details, it is recommended the reference Yakhot et al., (1992). The k -ω model also showed good performance, mainly due to its good characteristics in the evaluation of the boundary layer near the boundaries of the profile.

Natural convection between parallel plates
The results for the temperature distribution and the velocity field in a position that is half the height of the gap between plates are represented in Figures 15 and 16, where the percentage of 50% serves only to indicate that the position corresponds to half of the distance between the plates. The experimental results were obtained by Betts and Bokhari (2000).