Turbulence models impact on the flow and thermal analyses of jet impingement

. Accurate numerical reconstruction of heat and mass transfer processes in particular applications, such a jet impingement, is difficult to obtain even with the use of modern computational methods. In the proposed paper, the flow and thermal phenomena occurring during single minijet impingement on the flat, concave and convex, heated surfaces were considered. Problem of impingement on non-flat surface, still not common and purely described in the literature, can be of big importance in engineering applications, such as the heat exchangers. Numerical analyses, based on the mass, momentum and energy conservation laws, were conducted with the OpenFOAM software. Focus was placed on the proper model construction, in which turbulence and boundary layer modelling was crucial, due to their significance in the heat transfer processes. Analysis of results obtained by RANS models focused mostly on the comparison of turbulent and hydrodynamics parameters.


Introduction
Modern world and its technical development is connected strictly with the trials to improve the ways of energy resources utilization. One of the trends is to apply methods known from one branch of engineering in other. An example can be a jet impingement. It was being used in the cooling systems of electronics or in some metallurgical processes. Recently, however, it was proposed in the novel construction of cylindrical heat exchanger [1], which resulted in very promising values of transferred heat rates and lower than in other similar devices hydraulic resistance [2][3].
In abovementioned device, about 1000 orifices were generating minijets with the core diameter of ~1 mm. Due to complexity of the system, experimental investigation of it can be very difficult to perform, especially, when the detailed flow phenomena in the small scale is the main goal of research. Until now only general data was available, then, concerning macroscale parameters, such as already mentioned pressure losses or overall heat transfer efficiency. Understanding of the phenomena occurring in the device requires more detailed analysis, concerning also microscale.
For such cases numerical methods are very helpful. However, even now they demonstrate lack of accuracy in some scientific problems. Jet impingement unfortunately is one of the examples. According to Zuckerman [4] proper predicting of flow and thermal parameters in the case of jet impingement and its numerical modelling depends significantly on the type of turbulence model chosen for the simulation. In his paper it can be found, that the best results for RANS simulations can be obtained with v 2 -f model, which is the four-equation, enhanced k-ε type, model. While hydrodynamics of the phenomena is predicted with decent level of accuracy by RANS methods, heat transfer causes problems. Hadžiabdić [5] in his doctoral thesis confirmed this statement, moreover he concluded, that even LES methods exhibited lack of accuracy when used to predict heat transfer in jet impingement.
Before analyzing complex systems, it is essential to identify the advantages and disadvantages of possible numerical methods to be applied. The following paper regards the analysis of single minijet that impinge various surfaces -flat, convex and concave. Their geometry was chosen in the basis of previous research [6], as well as data from [1][2][3]. While the available information for the flat case is generally broad, non-flat impingement is still not a popular and widely described scientific topic. Correlations do not exist, which can be used in the problems regarded in [1][2][3]. The goal of paper was to identify the turbulence model impact on the results and describe the differences between impingement on various surfaces -since such knowledge would be essential to correctly analyze arrays of minijets, existing for example in heat exchangers [1][2][3]].

Mathematical model and geometry
Steady-state, single-phase, two-dimensional axisymmetric analyses were performed. High resolution, second order discretization schemes were applied to provide sufficient results. The conservation laws of mass (Equation 1), momentum (Equation 2) and energy (Equation 3) were applied, using Reynolds averaging approach. All variables marked with overline, such as (), would represent time-averaged value. They were coupled with various turbulence models, both high-(k-ε and v 2 -f) and low-Reynolds (SST k-ω) types: where u is the velocity, m/s; ρ is the density, kg/m 3 ; p is the pressure, Pa; μ is the dynamic viscosity, Pa·s; Sij is the strain rate tensor, 1/s, '' ij uu is the Reynolds stress term, m 2 /s 2 ; E is the total energy, J; T is the temperature, K; λef is the effective thermal conductivity, W/(m·K).    As it can be seen, all cases were axisymmetric. Moreover, values of H, which was the height between orifice exit and stagnation point, and R, which was the radius of curvature, were chosen to be multiplications of orifice diameter D. Multiplication factors used in the paper are listed in Table 1.

Flat surface, validation case
Numerical analyses of jet impingement could be compared for example with the ERCOFTAC benchmark cases, describing the experimental work by Cooper et al. [7] and Yan [8]. They represented the situation, in which air is impinging the flat surface, heated with constant heat flux. Orifice to surface distance H was equal to two diameters of the orifice, 2D. Heat flux at the surface was equal to 1000 W/m 2 . Reynolds number, defined as: at the exit of the orifice was equal to 23000, and, in addition, the flow there was fully developed, which was achieved with the mapped inlet boundary condition. Parameter D in Equation 4 is the jet impingement characteristic length, denoting the orifice diameter. Table 1 presents the boundary conditions for validation case. They were also used for non-flat cases, from Section 4. Numerous publications were based on those results, among which the papers by Behnia et al. [9][10] were chosen as the reference case. Both papers confirmed, that v 2 -f model presents the best performance when analyzing jet impingement. Moreover, the main drawback of standard k-ε model, which is overproduction of turbulence kinetic energy in the stagnation zone, was also confirmed. Results obtained with v 2 -f model were compared with standard k-ε and SST k-ω models. First two were used in low-Reynolds mode, while SST k-ω was used in the combination with wall functions, to check their performance.
Initially, for the validation case, it was important to fulfil the mesh requirements regarding jet impingement simulations. Mesh construction process was time consuming, required not only knowledge of phenomena but also the trial-error method. Regular mesh independence checks were performed as well. The final space division was chosen in the basis of the comparison of obtained results with [7][8][9][10]. Table 2 presents number of mesh elements that were chosen to analyze the process.
where α is the convective heat transfer coefficient, W/(K·m 2 ) and λ is the thermal conductivity, W/(K·m).
As can be seen, depending on the mesh size and settings, very different results were obtained. Figure 5 shows selected results from Figure 4, compared with benchmark data, experimental results by Yan [8] and numerical results by Behnia et al. [9][10]. Also results obtained with SST k-ω are presented there, to check the performance of wall-functions. Only the values calculated with Mesh 2 were included, as they exhibited the best agreement with [8][9][10]. Moreover, results presented in Section 4 were also obtained with utilization of Mesh 2 construction process.
As it can be noticed, far from the stagnation region all presented results are almost the same. However, large discrepancies occur in the stagnation zone. The k-ε model overpredicted the heat transfer significantlybut this effect was expected. On the other hand results obtained with OpenFOAM v 2 -f model also did not reflect the ones obtained by Behnia et al. [9][10]. Especially, in the region, where the distance from the stagnation point is slightly higher than the orifice radius. The answer for that discrepancy can be found in the paper of Billard and Lawrence [11]. From theoretical point of view they analyzed evolution of the v 2 -f models, because, depending on particular implementations, various results could be obtained. They described, that some models were adjusted to use very robust Dirichlet boundary condition at the wall for the elliptic relaxation function, f. It led to the possibility of their usage in the segregated solvers, commonly applied in many commercial codes or OpenFOAM. However, such method led to omitting one term in the relaxation equation and, as a result, another overpredictionof the velocity scale, v 2 . It explains the issue with Nusselt number values, presented in Figure 5. In the papers by Behnia et al. [9][10], their own solver was used, so they were able to use the original v 2 -f model, not having such drawbacks. Still, though, the difference between their results and the experiment can be noticed that may possibly be never fully avoided, as described in [4].
Another reason of presented differences may lay in the characteristics of the fully developed flow. Behnia et al. in [9] included Figure 14, in which the impact of different orifice-exit velocity and turbulence profiles on the Nusselt number distribution at the impinged surface was very significant.
Another conclusion can be related with results obtained with SST k-ω model. In that case the heat transfer was underpredicted in the stagnation region and proper far from it. However, because of the wall functions consideration, which are just a simplification of actual heat and mass transfer processes, this model would not be analyzed and described in the next sections. Apart from the thermal parameters, also a flow prediction by particular numerical models is very important for results validation. Figure 6 presents obtained numerically and compared with experimental [5] profiles of the velocity in various locations far from the stagnation point, x/D. At the bulk Reynolds number in the orifice equal to 23000, bulk air velocity ub was equal to 34.5 m/s. Analysis of Figure 6 leads to conclusion, that in contrary to thermal parameters, flow behavior was predicted well by the v 2 -f model implemented in OpenFOAM and used in presented studies. Differences between it and standard k-ε are also clearly visible.
Profiles of velocity obtained in [9][10] were very similar, including the noticeable difference between the experimental and numerical results starting at ratio H/D higher than 0.2, for x/D = 1 and 2.5. In Figure 7, the turbulence kinetic energy distribution for validation case is presented. Its budget can be written as follows: where μt is the eddy viscosity; Pa·s, ε is turbulence dissipation rate, m 2 /s 3 . Two terms, production and dissipation, are emphasized, as they were used in Section 4 for data presentation. The k-ε model was characterized by overproduction of turbulence kinetic energy in the stagnation region. Its maximum was located there. On the other hand, usage of v 2 -f model caused the maximum to move outside this region. It reflects the real-life situation that can be observed in [5,[9][10]. As mentioned at the beginning of Section 3, the most common implementation of v 2 -f model is not able to limit the extensive production of turbulence in stagnation zone properly. That is the reason of relatively high values of turbulence kinetic energy, visible in Figure 7(b), which do not occur in [5,[9][10]. Nevertheless, the v 2 -f model was chosen for the next analyses, presented in Section 4. Before analyzing arrays of jets, that impinge non-flat surfaces, it is important to define the influence of the surface shape on the flow behavior. In [6] Authors tried to define the critical curvature radius to orifice diameter ratio R/D, for which the curvature effect plays a role. It was concluded, that depending on the type of surface and mentioned ratio, the difference in heat transfer occurs between the particular non-flat and flat surface impingement. Moreover, with increasing ratio, the values tended to vary lesssince the curvature of stagnation zone was almost negligible. However, [6] did not contained hydrodynamic data, which also is very important. In [6] it was noticed, that for ratio of surface curvature and orifice diameter equal to 4, the curvature effect is noticeable. Therefore, this ratio was selected for representation in the following studies. In Figure 8, comparison of normalized (in the same manner, as in Figure 5) velocity profiles, depending on the type of surface: flat, concave or convex, is presented. For the radial distance x/D higher than 1 only slight differences occur. However, when analyzing results in Figure 8(a) and 8(b), it can be seen, that depending on the type of surface, the flow behaves in different way. Therefore, the stagnation zone is the area, where the most noticeable discrepancies takes place. It is especially visible in Figure 8(b), for height values H/D ≥ 0.25. To be able to provide data described above, it was important to propose the method of comparison between flat and non-flat surfaces. For that purpose, the distance x was measured as the straight chord connecting the stagnation point and particular point on the curve. Moreover it was measured, that difference between its length and the length of the arc connecting both points was negligible.  Table 3. In Figure 9, distribution of the turbulence kinetic energy during impingement on concave, Figure 9(a) and convex, Figure 9(b), surfaces is shown. Its maxima are located outside the stagnation zone, as in the flat case. In Table 3, locations of turbulence kinetic energy maxima are listed, in relation with orifice diameter D. They were established for flat and non-flat cases in the same manner, as data from Figure 8. Differences between each case existit is an important aspect to mention, because it may influence the results when the jets array impingement occurs, such as in [1][2][3]. Multiple jets might strongly interfere with each other. In Figure 10, two terms of turbulence kinetic energy budget are presented, along the curvature or segment distant from stagnation point, at height normal to impinged surface, where maxima from Table 3 occurred. In the stagnation zone, the highest production exists in the convex case, followed by flat and concave ones. For the dissipation, situation is opposite. However, in the regions close to the highest values of turbulence kinetic energy, located at the distances x/D mentioned in Table  3, for both production and dissipation highest values were obtained for convex surface. In general, behavior presented in plots for each situation is quite similar.

Summary
In this paper, thermal and hydrodynamic analyses of jet impingement on flat and non-flat surfaces were presented. Both boundary conditions and scope of interest were based on [1][2][3][5][6][7][8]. Mesh and software configurations were determined by ERCOFTAC data [7][8]. Velocity profiles and turbulence kinetic energy budgets comparison revealed important differences. In Authors opinion those variations would matter, when the whole jets array impinging flat and non-flat surfaces would be analyzed. Different turbulence model were considered, however, as it was proved in [4], v 2 -f has given the best results. While RANS models, presented in the following paper, can reveal important information, they should be verified and extended also by more extensive methods, such as LES approach.