Numerical Analysis of Turbulent Heat Transfer in the Case of Minijets Array

The presented numerical investigations show an analysis of the turbulent single-phase array of ten minijets impinging a heated surface, which lead to the intensification of heat transfer between the fluid and the surface. Attention was devoted to the comparison between phenomena occurring for the heated flat and concave surface geometry. The selection of the shapes was based on the impinging jets applications. From the numerical point of view, the focus was placed on a comparison of the Reynolds Averaged Navier–Stokes (RANS) turbulence model implementations in ANSYS Fluent software, and their impact on the modeling precision of the thermal and hydrodynamic boundary layers phenomena. The 3D numerical model was based on the continuity, momentum, and energy transport equations, together with three various RANS turbulence models: k-ω SST Kato-Launder, k-ε RNG Kato-Launder, and Intermittency Transition. The water submerged minijets, characterized by three various values of Reynolds number, were considered. Average surface Nusselt number values for all analyzed cases were compared with the experimental correlations and exhibited the same tendency but differed in detail. Numerically obtained average Nusselt number values agreed with the results of two from three correlations in the range of 10–20%. The flat surface was characterized by higher heat transfer than the concave one and an influence of the cross flow, changing the symmetrical distribution of the Nusselt number, was more visible for it. A cross flow impact was found in fuzzy hexagonal or tetragonal symmetry of this distribution. Additionally, the areas of high temperature gradient values were identified in the region of the strongest jets’ interactions, which can be important for mechanical strength analysis.


Introduction
An efficient heat transfer relates to better performance of the cooling/heating equipment, lower costs of operation, protection of the environment, and may lead to the manufacturing of smaller devices. The utilization of jet impingement phenomena as the heat transfer enhancement technology is well known and commonly applied (glass and metal industry, gas turbine cooling, cryosurgery, paper drying, etc. [1][2][3][4][5][6][7]), but there is still increasing interest in its development and new applications.
Considering the heat rate values that need to be withdrawn, the mentioned technology can be reduced to a single jet or to be extended toward the multijets array. Numerical modeling of the first one is relatively easy and can be a source of data for a deeper understanding of the heat transfer processes, occurring within the boundary layer. From the application point of view, only the jets array is able to transfer high heat fluxes, from real devices. One of the examples is connected with the heat exchangers, small in size, but very efficient. There are many ways to intensify the heat transfer processes in such devices, for example, through the application of nanofluids [8], application of jets [9,10], increasing the heat transfer surface and modifying the flow structure [11], or increasing the heat transfer surface Experimental investigations are very important for prototype testing and for numerical analysis as a source of crucial information (reference data). Unfortunately, experimental investigations are time-consuming and can be very expensive. That is why numerical studies are the cornerstone of modern engineering. The utilization of Computational Fluid Dynamics (CFD) can improve the design process and speed up the overall development time. Moreover, numerical analysis can provide detailed information about the studied phenomenon. However, both approaches should not be considered separately, because both are important: experimental approach as a reference case and numerical one as a case for development. Additionally, an accurate and fast numerical model is essential for the future optimization process of the real device. Therefore, the preparation of a suitable numerical model is a fundamental process.
The issue of the turbulent jet impingement phenomenon has been studied for quite a long time, but its numerical representation is still challenging. The modeling of the hydrodynamic phenomena seems to be well known, since in this aspect the numerical and experimental results exhibit good agreement [15]. On the other hand, the analysis of heat transfer processes still requires improvement, because the Nusselt number distribution shows some discrepancies. Additionally, jet impingement in a majority of the cases is a turbulent phenomenon, which brings to light another challenge for the numerical analysis-the selection of the proper turbulence treatment. One of the examples of such work is the publication of Yang et al. [16]. They analyzed the phenomenon of jets array impingement on the non-flat surface and included the comparison between the performance of various turbulence models. Additionally, they conducted an experimental analysis to verify the numerical data. The authors claimed that the best choice would be the SST k-ω model, but its experimental validation was not fully successful. Another publication, worth to mention, is the paper by Deng et al. [17]. The numerical analysis of high-power laser diode cooling, using the water jets array impingement, was also based on the k-ω model. It allowed them to obtain a detailed analysis of the phenomena, nevertheless, the validation part did not include a comprehensive discussion on the turbulence models and their impact on the results. Similarly, Sabato et al. [18] focused on the performance of the SST k-ω model to represent the turbulent flow and heat transfer in the case of jets. In comparison to previous articles the most interesting is the outcome of the work by Draksler [19]. In his numerous articles, he focused on the Large Eddy Simulations, but the outcome of work [19] gives many advices on the turbulence model evaluation methods. Additionally, the LES based analysis of the jets array impingement might be a valuable verification database to check the feasibility of RANS based models.
The presented paper describes the continuation of research work started by Kura [20], in which the focus was placed on the phenomena occurring in the hydrodynamic and thermal boundary layers Experimental investigations are very important for prototype testing and for numerical analysis as a source of crucial information (reference data). Unfortunately, experimental investigations are time-consuming and can be very expensive. That is why numerical studies are the cornerstone of modern engineering. The utilization of Computational Fluid Dynamics (CFD) can improve the design process and speed up the overall development time. Moreover, numerical analysis can provide detailed information about the studied phenomenon. However, both approaches should not be considered separately, because both are important: experimental approach as a reference case and numerical one as a case for development. Additionally, an accurate and fast numerical model is essential for the future optimization process of the real device. Therefore, the preparation of a suitable numerical model is a fundamental process.
The issue of the turbulent jet impingement phenomenon has been studied for quite a long time, but its numerical representation is still challenging. The modeling of the hydrodynamic phenomena seems to be well known, since in this aspect the numerical and experimental results exhibit good agreement [15]. On the other hand, the analysis of heat transfer processes still requires improvement, because the Nusselt number distribution shows some discrepancies. Additionally, jet impingement in a majority of the cases is a turbulent phenomenon, which brings to light another challenge for the numerical analysis-the selection of the proper turbulence treatment. One of the examples of such work is the publication of Yang et al. [16]. They analyzed the phenomenon of jets array impingement on the non-flat surface and included the comparison between the performance of various turbulence models. Additionally, they conducted an experimental analysis to verify the numerical data. The authors claimed that the best choice would be the SST k-ω model, but its experimental validation was not fully successful. Another publication, worth to mention, is the paper by Deng et al. [17]. The numerical analysis of high-power laser diode cooling, using the water jets array impingement, was also based on the k-ω model. It allowed them to obtain a detailed analysis of the phenomena, nevertheless, the validation part did not include a comprehensive discussion on the turbulence models and their impact on the results. Similarly, Sabato et al. [18] focused on the performance of the SST k-ω model to represent the turbulent flow and heat transfer in the case of jets. In comparison to previous articles the most interesting is the outcome of the work by Draksler [19]. In his numerous articles, he focused on the Large Eddy Simulations, but the outcome of work [19] gives many advices on the turbulence model evaluation methods. Additionally, the LES based analysis of the jets array impingement might be a valuable verification database to check the feasibility of RANS based models. The presented paper describes the continuation of research work started by Kura [20], in which the focus was placed on the phenomena occurring in the hydrodynamic and thermal boundary layers of single jet and jets array. All analyses were conducted with utilization of three various RANS turbulence models, namely, k-ω SST Kato-Launder (k-ω SST KL), k-ε RNG Kato-Launder (k-ε RNG KL), and Intermittency Transition and three values of Reynolds number. An additional aim of the investigations was connected with evaluation of their suitability for the jet impingement modeling. At the current state of investigations, the local values of Nusselt number at flat and concave surfaces were calculated. Moreover, taking into account temperature gradient pattern, potentially dangerous regions, from a mechanical point of view, were identified. The average values of the Nusselt number was used as the measure of flat and concave jets impingements performance.

Mathematical Model
Assuming that the flow was Newtonian, incompressible, without additional mass, momentum, and energy sources, the mathematical model consisted of three transport equations: mass, momentum, and energy, defined as follows [21]: where U i,j are the mean velocity components, m/s; ρ is the density, kg/m 3 ; P is the pressure, Pa; µ is the dynamic viscosity, Pa·s; S ij is the strain rate tensor, 1/s; u i u j is the Reynolds stress term, m 2 /s 2 ; Θ is the mean temperature, K; a is the thermal diffusivity, and m 2 /s and u j θ is the turbulent heat flux. In order to solve turbulent phenomena, the additional equations of turbulence model must be considered. The preliminary stage was devoted to selection of the most suitable model.
The numerical results of preliminary study of this research, devoted to a single jet analysis, show that the perfect turbulence model does not exist. The numerical calculations have been performed for 2D axisymmetric model at various values of the Reynolds numbers and non-dimensional geometrical parameters H/D. The Reynolds number was defined as follows [21]: where ρ is the fluid density, kg/m 3 ; U is the average orifice velocity, m/s; D is the orifice diameter, m; and µ is the fluid dynamic viscosity, Pa·s. The geometric concept for the preliminary study is presented in Figure 2.
The numerical results of the Nusselt number distribution along the heated wall (as a function of non-dimensional parameter x/D) for various turbulence models have been compared with the experimental data. Local Nusselt number was defined as follows [21]: where α is the heat transfer coefficient, W/m 2 K; D is the orifice diameter, m; λ is the fluid thermal conductivity, W/mK; q is the wall heat flux, W/m 2 ; T w is the local wall temperature, K; and T in is the reference inlet temperature, K.  An example of a comparison between numerical and experimental [2,3,5,6] investigations can be found in Figure 3.  For all the analyzed models, discrepancies between numerical and experimental data were observed. The k-ω models family (especially Intermittency Transition) gave the results, which matched the Nusselt number values better than the results of other models for low values of H/D (where H is the distance between orifice inlet and the impinged surface, D is the orifice diameter, e.g., equal to 2). For higher values of H/D parameter (e.g., equal to 6), the k-ε model was the one, which results were the closest to the experimental data (the secondary maximum in Nusselt number does not occur). Preliminary investigations confirmed that it is very difficult to select just one turbulence model that is adequate for all analyzed cases in all the parameter ranges. Therefore, three turbulence models were chosen for jets array analysis: k-ω SST KL, k-ε RNG KL, and Intermittency Transition Models. The first two models require two additional equations to be solved (k and ω/ε), the last one needs three other equations (k, ω, and intermittency) [21].

Jets Array
Calculations were conducted for three volumetric flow rates (100, 75, 50 L/h) taken from [10,13] and represented by a Reynolds number value.
Two geometric configurations were investigated: jets array impinging flat surface (Variant A) and concave surface (Variant B). The second one is based on the existing heat exchanger [10,13]. An example of a comparison between numerical and experimental [2,3,5,6] investigations can be found in Figure 3. An example of a comparison between numerical and experimental [2,3,5,6] investigations can be found in Figure 3.  For all the analyzed models, discrepancies between numerical and experimental data were observed. The k-ω models family (especially Intermittency Transition) gave the results, which matched the Nusselt number values better than the results of other models for low values of H/D (where H is the distance between orifice inlet and the impinged surface, D is the orifice diameter, e.g., equal to 2). For higher values of H/D parameter (e.g., equal to 6), the k-ε model was the one, which results were the closest to the experimental data (the secondary maximum in Nusselt number does not occur). Preliminary investigations confirmed that it is very difficult to select just one turbulence model that is adequate for all analyzed cases in all the parameter ranges. Therefore, three turbulence models were chosen for jets array analysis: k-ω SST KL, k-ε RNG KL, and Intermittency Transition Models. The first two models require two additional equations to be solved (k and ω/ε), the last one needs three other equations (k, ω, and intermittency) [21].

Jets Array
Calculations were conducted for three volumetric flow rates (100, 75, 50 L/h) taken from [10,13] and represented by a Reynolds number value.
Two geometric configurations were investigated: jets array impinging flat surface (Variant A) and concave surface (Variant B). The second one is based on the existing heat exchanger [10,13]. For all the analyzed models, discrepancies between numerical and experimental data were observed. The k-ω models family (especially Intermittency Transition) gave the results, which matched the Nusselt number values better than the results of other models for low values of H/D (where H is the distance between orifice inlet and the impinged surface, D is the orifice diameter, e.g., equal to 2). For higher values of H/D parameter (e.g., equal to 6), the k-ε model was the one, which results were the closest to the experimental data (the secondary maximum in Nusselt number does not occur). Preliminary investigations confirmed that it is very difficult to select just one turbulence model that is adequate for all analyzed cases in all the parameter ranges. Therefore, three turbulence models were chosen for jets array analysis: k-ω SST KL, k-ε RNG KL, and Intermittency Transition Models. The first two models require two additional equations to be solved (k and ω/ε), the last one needs three other equations (k, ω, and intermittency) [21].

Jets Array
Calculations were conducted for three volumetric flow rates (100, 75, 50 L/h) taken from [10,13] and represented by a Reynolds number value.
Two geometric configurations were investigated: jets array impinging flat surface (Variant A) and concave surface (Variant B). The second one is based on the existing heat exchanger [10,13]. Geometrical Geometrical configurations, together with boundary conditions, are presented in Figure 4, while Figure 5 shows all 18 studied cases.
(a) (b)  The geometrical configurations ( Figure 4) represent the computational domains, which are only 1/4 of entire model. Water was considered as a working fluid for both geometries. In both variants (A and B), the inlet of working fluid is located in the bottom part of the geometry. Water flows through presented orifices creating a jet, which impinge heated flat or concave surfaces, and then is directed toward one outlet, since the second one is blocked. Such geometry causes the appearance of cross flow. One computational case considered the system without cross flow occurrence and in this case, there were two outlets. Both types of heat transfer surfaces have a length of 21 mm. All orifices have a diameter equal to 1 mm. The distance between orifices' exits and heated surfaces is equal to 2 mm (H/D = 2) for both Variants, while distance between the orifices' centers is about 4.2 mm (S/D = 4.2). The geometry of Variant B has a ratio of surface curvature to the orifice diameter (R/D) equal to 8. Constant thermal properties of water: heat capacity Cp = 4182 J/kgK, density ρ = 998.2 kg/m 3 , thermal conductivity λ = 0.6 W/mK, dynamic viscosity μ = 0.001003 Pa·s, were applied in conducted analyses.
A fully developed profile (controlled velocity and turbulence parameters: k, ε, ω, intermittency) at the inlet was used for all studied cases. Detailed information about boundary conditions is presented in Table 1.    The geometrical configurations ( Figure 4) represent the computational domains, which are only 1/4 of entire model. Water was considered as a working fluid for both geometries. In both variants (A and B), the inlet of working fluid is located in the bottom part of the geometry. Water flows through presented orifices creating a jet, which impinge heated flat or concave surfaces, and then is directed toward one outlet, since the second one is blocked. Such geometry causes the appearance of cross flow. One computational case considered the system without cross flow occurrence and in this case, there were two outlets. Both types of heat transfer surfaces have a length of 21 mm. All orifices have a diameter equal to 1 mm. The distance between orifices' exits and heated surfaces is equal to 2 mm (H/D = 2) for both Variants, while distance between the orifices' centers is about 4.2 mm (S/D = 4.2). The geometry of Variant B has a ratio of surface curvature to the orifice diameter (R/D) equal to 8. Constant thermal properties of water: heat capacity Cp = 4182 J/kgK, density ρ = 998.2 kg/m 3 , thermal conductivity λ = 0.6 W/mK, dynamic viscosity μ = 0.001003 Pa·s, were applied in conducted analyses.
A fully developed profile (controlled velocity and turbulence parameters: k, ε, ω, intermittency) at the inlet was used for all studied cases. Detailed information about boundary conditions is presented in Table 1.  The geometrical configurations ( Figure 4) represent the computational domains, which are only 1/4 of entire model. Water was considered as a working fluid for both geometries. In both variants (A and B), the inlet of working fluid is located in the bottom part of the geometry. Water flows through presented orifices creating a jet, which impinge heated flat or concave surfaces, and then is directed toward one outlet, since the second one is blocked. Such geometry causes the appearance of cross flow. One computational case considered the system without cross flow occurrence and in this case, there were two outlets. Both types of heat transfer surfaces have a length of 21 mm. All orifices have a diameter equal to 1 mm. The distance between orifices' exits and heated surfaces is equal to 2 mm (H/D = 2) for both Variants, while distance between the orifices' centers is about 4.2 mm (S/D = 4.2). The geometry of Variant B has a ratio of surface curvature to the orifice diameter (R/D) equal to 8. Constant thermal properties of water: heat capacity C p = 4182 J/kgK, density ρ = 998.2 kg/m 3 , thermal conductivity λ = 0.6 W/mK, dynamic viscosity µ = 0.001003 Pa·s, were applied in conducted analyses.
A fully developed profile (controlled velocity and turbulence parameters: k, ε, ω, intermittency) at the inlet was used for all studied cases. Detailed information about boundary conditions is presented in Table 1. The list of all assumptions used in the described investigation is as follows: For all the analyzed cases the average Nusselt number was calculated and compared with the experimental correlations available in the literature [22][23][24]. The average Nusselt number was calculated by area-weighted averaging method, it ensured that the average value was mesh cell size independent.
Numerical calculations were carried out in ANSYS Fluent software, in a high performance computer in CYFRONET AGH.

Mesh
For Variant A the block-structured (about 4.9 M elements) and for Variant B hybrid structured-unstructured (about 4.9 M elements) meshes have been generated. Figure 6 presents mesh views for both geometries. The list of all assumptions used in the described investigation is as follows: • incompressible flow • constant material properties • symmetry boundary condition • all walls adiabatic except the heated surface • fully developed inlet velocity profile with constant temperature For all the analyzed cases the average Nusselt number was calculated and compared with the experimental correlations available in the literature [22][23][24]. The average Nusselt number was calculated by area-weighted averaging method, it ensured that the average value was mesh cell size independent.
Numerical calculations were carried out in ANSYS Fluent software, in a high performance computer in CYFRONET AGH.

Mesh
For Variant A the block-structured (about 4.9 M elements) and for Variant B hybrid structuredunstructured (about 4.9 M elements) meshes have been generated. Figure 6 presents mesh views for both geometries. Additionally, for all walls in Variants A and B, the inflation layer has been applied to ensure proper modeling of all near-wall phenomena. Information about y+ can be found in Table 2. The average Nusselt number on the heated surface was selected as a parameter for mesh sensitivity study. An example of such a study is presented in Table 3. Version 3 (v3) of the mesh was selected for further calculations. Additionally, for all walls in Variants A and B, the inflation layer has been applied to ensure proper modeling of all near-wall phenomena. Information about y+ can be found in Table 2. The average Nusselt number on the heated surface was selected as a parameter for mesh sensitivity study. An example of such a study is presented in Table 3. Version 3 (v3) of the mesh was selected for further calculations.

Numerical Procedure
The numerical procedure was the same for all cases (Figure 7). First, the steady-state calculation was carried out, and then the transient one. A second-order upwind scheme was used for all equations. For temporal discretization, first-order scheme was used, together with Coupled Pressure-Velocity Algorithm, which solves the momentum and pressure-based continuity equations together [21]. The time step was equal to 10 −6 s for all cases. All residuals below 10 −6 and stabilization of average and maximum Nusselt number values on the heated surface were used as the convergence criterions.

Numerical Procedure
The numerical procedure was the same for all cases (Figure 7). First, the steady-state calculation was carried out, and then the transient one. A second-order upwind scheme was used for all equations. For temporal discretization, first-order scheme was used, together with Coupled Pressure-Velocity Algorithm, which solves the momentum and pressure-based continuity equations together [21]. The time step was equal to 10 −6 s for all cases. All residuals below 10 −6 and stabilization of average and maximum Nusselt number values on the heated surface were used as the convergence criterions.

Results
The local Nusselt number distribution for all studied cases is presented in Figures 8-10. The results for Variant B are projected on the plane, therefore its surface is larger.

Results
The local Nusselt number distribution for all studied cases is presented in Figures 8-10. The results for Variant B are projected on the plane, therefore its surface is larger.

Numerical Procedure
The numerical procedure was the same for all cases (Figure 7). First, the steady-state calculation was carried out, and then the transient one. A second-order upwind scheme was used for all equations. For temporal discretization, first-order scheme was used, together with Coupled Pressure-Velocity Algorithm, which solves the momentum and pressure-based continuity equations together [21]. The time step was equal to 10 −6 s for all cases. All residuals below 10 −6 and stabilization of average and maximum Nusselt number values on the heated surface were used as the convergence criterions.

Results
The local Nusselt number distribution for all studied cases is presented in Figures 8-10. The results for Variant B are projected on the plane, therefore its surface is larger. (e) (f)  The areas of jet impingement are clearly visible. The Nusselt number distribution slightly differs between Variants A and B. Cross flow influence on the heat transport phenomena can be easily identified through the elongated shapes of areas of intensified heat transfer. It looks like the k-ω SST KL and Intermittency Transition models exhibit similar tendency and their utilization leads to the higher values of the Nusselt number. In the case of single jet (Figure 3), they were in a better agreement with experimental data in the stagnation region, while the k-ε RNG KL model's results were characterized by slightly lower values of the Nusselt number. The same trend can be observed for the jets array. A very convenient way of analyzing real devices is to use average values, therefore, numerically calculated average Nusselt number values on the heated surface were compared with the experimental correlations presented in the literature [22][23][24]. The correlations' equations are collected in Table 4, in which S is the distance between orifices, H is the distance between orifice exit and heated surface, D is the orifice diameter, Pr is the Prandtl number, and K, G, F, and f are the equation coefficients.

Reference Correlation Limitations
Martin [22] Nu = Pr 0.42 KGF A comparison between numerical analysis and correlation data is presented in Figure 11. and heated surface, D is the orifice diameter, Pr is the Prandtl number, and K, G, F, and f are the equation coefficients. A comparison between numerical analysis and correlation data is presented in Figure 11.

Reference Correlation Limitations
Martin [22] KGF 42 . 0 Pr = Nu 05 . 0 6 ) ) / 6 . 0 Robinson [23]  An analysis of the results shows that the average Nusselt number values are similar for k-ω SST KL and Intermittency Transition models, but for k-ε RNG KL the discrepancies are higher than for kω model family of about 20%. Another observation is that even CFD results differ from correlation data, all numerical Nusselt number values are within the empirical values range. The heat transfer enhancement is higher for Variant A than for Variant B of about 15%. The average values of Nusselt number show that all presented turbulence models can be used for this type of analysis, when the An analysis of the results shows that the average Nusselt number values are similar for k-ω SST KL and Intermittency Transition models, but for k-ε RNG KL the discrepancies are higher than for k-ω model family of about 20%. Another observation is that even CFD results differ from correlation data, all numerical Nusselt number values are within the empirical values range. The heat transfer enhancement is higher for Variant A than for Variant B of about 15%. The average values of Nusselt number show that all presented turbulence models can be used for this type of analysis, when the amount of transferred heat rate is important. Considering the empirical correlations and their agreement with the numerical results, the best correspondence can be found for the Fabbri's one.
More in-depth analysis reveals a strong impact of the crossflow phenomenon on the local Nusselt number distribution. The deformation of local Nusselt number distribution in jets array in reference to single jet case is presented in Figure 12.
amount of transferred heat rate is important. Considering the empirical correlations and their agreement with the numerical results, the best correspondence can be found for the Fabbri's one.
More in-depth analysis reveals a strong impact of the crossflow phenomenon on the local Nusselt number distribution. The deformation of local Nusselt number distribution in jets array in reference to single jet case is presented in Figure 12. A comparison between Variant A and Variant B for all studied turbulence models is presented in Figure 13. In general, locations of the maxima for all turbulence models are clearly visible and more or less are the same. However, the agreement between k-ω SST KL and Intermittency Transition models is particularly good both in value and shape. The k-ε RNG KL model exhibits some discrepancies, especially in the case of flat surface, for which it predicts additional maximum.
One of the main purposes of heat transfer numerical investigations, besides optimization and performance analysis, is information about temperature gradients. This information is crucial for future structure analysis, like thermal stresses calculations. The prediction of thermal stresses in real device is important for design phase, manufacturing, and for the further operation of the heat exchanger. A comparison of temperature gradient values between Variant A and Variant B for all studied turbulence models is presented in Figure 14. The values represent two components of temperature gradient, tangential to the heated surface. On the flat surface the temperature gradient takes higher values than in the case of concave one. The decreasing value of flow rate causes increase in the temperature gradient value. It suggests that an influence of the crossflow, in the case of higher flow rates values, reduces the temperature gradient values. All analyzed turbulence models predict the maxima of the temperature gradient in similar locations. A comparison between Variant A and Variant B for all studied turbulence models is presented in Figure 13. In general, locations of the maxima for all turbulence models are clearly visible and more or less are the same. However, the agreement between k-ω SST KL and Intermittency Transition models is particularly good both in value and shape. The k-ε RNG KL model exhibits some discrepancies, especially in the case of flat surface, for which it predicts additional maximum.
One of the main purposes of heat transfer numerical investigations, besides optimization and performance analysis, is information about temperature gradients. This information is crucial for future structure analysis, like thermal stresses calculations. The prediction of thermal stresses in real device is important for design phase, manufacturing, and for the further operation of the heat exchanger. A comparison of temperature gradient values between Variant A and Variant B for all studied turbulence models is presented in Figure 14. The values represent two components of temperature gradient, tangential to the heated surface. On the flat surface the temperature gradient takes higher values than in the case of concave one. The decreasing value of flow rate causes increase in the temperature gradient value. It suggests that an influence of the crossflow, in the case of higher flow rates values, reduces the temperature gradient values. All analyzed turbulence models predict the maxima of the temperature gradient in similar locations. Symmetry 2020, 12, x FOR PEER REVIEW 12 of 19      Temperature gradient analysis leads to an identification of potentially dangerous regions, from mechanical strength point of view. They are connected with the areas of the most intense interaction between jets, and which creates particular pattern on the heated wall. In the case of flat surface, it looks like broken/fuzzy symmetry of hexagonal cells, while in the case of concave one, of tetragonal cells. In these areas a special treatment should be applied. Their identification can also aid in the prediction of expected life of the device, which is important information for the customer and for environmental protection.
To evaluate the impact of cross flow on the heat transfer processes, one additional case was considered. It represented the system without appearance of crossflow, it means that the outlets were at both sides of the geometry. The results of local Nusselt number calculations are shown in Figure 18, while the temperature gradient tangential components magnitude in Figure 19. In the central part of Figure 18, the hexagonal cells symmetrical distribution is visible in the case of flat surface, and tetragonal one in the case of concave surface. The temperature gradient magnitude (Figure 19) exhibits the areas of the strongest jet interaction and at the same time accentuates the regular and symmetrical shape of the cells. These results confirm the influence of the cross flow and its significance for heat transfer processes.
prediction of expected life of the device, which is important information for the customer and for environmental protection.
To evaluate the impact of cross flow on the heat transfer processes, one additional case was considered. It represented the system without appearance of crossflow, it means that the outlets were at both sides of the geometry. The results of local Nusselt number calculations are shown in Figure  18, while the temperature gradient tangential components magnitude in Figure 19. In the central part of Figure 18, the hexagonal cells symmetrical distribution is visible in the case of flat surface, and tetragonal one in the case of concave surface. The temperature gradient magnitude (Figure 19) exhibits the areas of the strongest jet interaction and at the same time accentuates the regular and symmetrical shape of the cells. These results confirm the influence of the cross flow and its significance for heat transfer processes.

Summary
The presented research shows the numerical analyses of impinging jets array in the case of flat and concave surfaces. The studies concentrated on heat transfer aspects of this phenomenon and the comparison of average Nusselt number for jets array between numerical analysis and empirical correlations. The CFD results were presented for two different geometrical configurations, and three different volumetric flow rates. Three RANS turbulence models were analyzed. Based on the results from all preliminary studies and presented major cases, it can be concluded that all three turbulence models predict the average Nusselt number correctly within the range of values coming from the empirical correlations. However, the best agreement was found with the Fabbri's correlation. The flat surface was characterized by higher heat transfer performance than the concave one of about 15%. Moreover, the deformation of local Nusselt number distribution due to the cross flow in jets array is presented and compared with a single jet analysis. This deformation was found through the broken symmetry of hexagonal and tetragonal cells in the case of flat and concave surfaces, respectively. The potentially risky areas, from the mechanical strength point of view, were identified and they correspond to the mutual interaction between the jets. The difference in the tangential components of the temperature gradient characteristics of flat and concave surfaces was found. In the case of the flat surface, they took higher values than in the case of concave one. For both surfaces, they increased with decreasing values of the flow rate. To find the correspondence between the hydrodynamics, heat transfer, and mechanical strength of mentioned previously prototype heat exchanger, further analyses are necessary, and they should also be focused on the turbulence characteristics. In such analyses, k-ω SST KL and Intermittency Transition turbulence models should be considered, due to correspondence in prediction of thermal processes occurring in the vicinity of the heated wall.