Curved surface minijet impingement phenomena analysed with ζ-f turbulence model

: The jet impingement phenomenon plays an important role among the heat transfer intensiﬁcation methods. Very often its application and analyses refer to simple ﬂat surfaces, while there is a lack of information in the literature for cases addressing curved surfaces. In the present work, the single jet impingement on the non-ﬂat (concave and convex) surface is studied for a wide range of geometries, which originate from the mini-jet heat-exchanger design. The numerical simulations were performed by an advanced ζ - f turbulence model implemented in the open-source OpenFOAM (ESI-OpenCFD Ltd, Bracknell, United Kingdom) code. Noticeable differences in the phenomena occurring on the convex and concave surfaces were identiﬁed in the stagnation zone. Besides, the existence and location of the secondary peak in the Nusselt number distribution differed between the cases. These distributions were inﬂuenced by the shape of geometry, which determined ﬂow characteristics and resulting heat transfer performance. The origins of these differences were looked at in the turbulence characteristics close to the impinged surface of the stagnations zone and its vicinity, where turbulence kinetic energy and enstrophy were analysed. It was stated that the differences are already noticeable for the single jet impingement case, but they might sum up when multiple jets are considered. Therefore, here presented results would be important for the analysis of the overall unit of mentioned mini-jets heat-exchanger.


Flat Surface Jet Impingement
Between numerous turbulence-related cases that are considered numerically, the jet impingement is an example of the phenomenon that despite the development of various computational methods and many scientists involved, still causes problems, even though its analyses started as early as just after the Second World War [1]. Uncertain accuracy of available computational methods, especially related to the thermal results, was commented by Zuckerman and Lior [2] or Dewan et al. [3].
Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS)-based methods, previously applied for analyses of particular, simple cases of low Reynolds number, have already become available also for complex flows, including for example additional heat transfer or crossflow phenomena. However, it does not mean that the Reynolds Averaged Navier-Stokes (RANS) methods should be forgotten. They still are important due to significantly lower computational time, sufficient accuracy and common availability [2,3]. Moreover their advanced examples are step by step introduced with an ultimate goal of reaching the highest accuracy and possible universality.  [35][36][37]. Red arrows indicate the hot fluid flow, while blue arrows indicate a flow of cold fluid. Orifices, in which the jets are generated, are visible on the right side of the image. The solid wall, which constricts the cold fluid, is not shown, for clarity.
The presented results [35][36][37], significant from the engineering point of view, refer to the overall output. However, the description of the small-scale turbulent phenomena occurring near the heat transfer curved wall due to the jets' impingement, being the source of the heat transfer intensification, are not yet available. Such information will be of great importance, for the purposes of further optimisation of the similar devices. It requires distinction of numerical model, which will provide the best possible accuracy of the results. Previous studies, which have been already published [34], were based on the Intermittency Transition Model described by Menter et al. [38].
The results exhibited some issues referring to an accuracy and further research was necessary to make a proper choice. It was already performed and published in [39,40]. In both articles, which concerned single jet impingement on both flat and non-flat surfaces, the superior performance of the Hanjalic's ζ-f model over others investigated was proven. Additionally this model allows one to predict the turbulent-laminar transition effects that might occur when the jets array is considered [41]. The main conclusion coming from [39,40] was, that the hydrodynamic phenomena occurring in both cases (flat and non-flat geometries impingement) are very similar, so the turbulence model predicting the heat and mass transfer well for the flat case would also be suitable for the non-flat ones.
Presented paper concentrates on analysis of non-flat, namely concave and convex, impinging jet, from the hydrodynamic and heat transfer point of view. Special attention is devoted towards their relation in the impinged surface proximity. Following scheme, The presented results [35][36][37], significant from the engineering point of view, refer to the overall output. However, the description of the small-scale turbulent phenomena occurring near the heat transfer curved wall due to the jets' impingement, being the source of the heat transfer intensification, are not yet available. Such information will be of great importance, for the purposes of further optimisation of the similar devices. It requires distinction of numerical model, which will provide the best possible accuracy of the results. Previous studies, which have been already published [34], were based on the Intermittency Transition Model described by Menter et al. [38].
The results exhibited some issues referring to an accuracy and further research was necessary to make a proper choice. It was already performed and published in [39,40]. In both articles, which concerned single jet impingement on both flat and non-flat surfaces, the superior performance of the Hanjalic's ζ-f model over others investigated was proven. Additionally this model allows one to predict the turbulent-laminar transition effects that might occur when the jets array is considered [41]. The main conclusion coming from [39,40] was, that the hydrodynamic phenomena occurring in both cases (flat and non-flat geometries impingement) are very similar, so the turbulence model predicting the heat and mass transfer well for the flat case would also be suitable for the non-flat ones.
Presented paper concentrates on analysis of non-flat, namely concave and convex, impinging jet, from the hydrodynamic and heat transfer point of view. Special attention is devoted towards their relation in the impinged surface proximity. Following scheme, presented in Figure 2, simply defines the problem, which authors tried to solve and the main outcomes of this paper. presented in Figure 2, simply defines the problem, which authors tried to solve and the main outcomes of this paper.

Geometrical Model and Boundary Conditions
Two non-flat surface types were considered, concave and convex. Both of particular curvature radius R. The schematic view of important geometrical parameters is presented in Figure 3. Other necessary variables are as follows: the diameter D of the jet generating nozzle and distance H between the nozzle exit and stagnation point (the central point of the impinged surface). The analyses regarded the axisymmetric cases, so geometry is shown in such manner.
Applied boundary conditions were as follows: • at the inlet: the constant bulk velocity, determined from the Reynolds number values (namely Re = 5000; Re = 10,000; Re = 23,000) and the constant bulk temperature equal to 293 K, • at impinged wall: the constant heat flux equal to 1000 W/m 2 , • at the upper wall: the constant heat flux equal to 0 W/m 2 , • at the outlet: the gauge pressure equal to 0 Pa.

Geometrical Model and Boundary Conditions
Two non-flat surface types were considered, concave and convex. Both of particular curvature radius R. The schematic view of important geometrical parameters is presented in Figure 3. Other necessary variables are as follows: the diameter D of the jet generating nozzle and distance H between the nozzle exit and stagnation point (the central point of the impinged surface). The analyses regarded the axisymmetric cases, so geometry is shown in such manner.
Applied boundary conditions were as follows: presented in Figure 2, simply defines the problem, which authors tried to solve and the main outcomes of this paper.

Geometrical Model and Boundary Conditions
Two non-flat surface types were considered, concave and convex. Both of particular curvature radius R. The schematic view of important geometrical parameters is presented in Figure 3. Other necessary variables are as follows: the diameter D of the jet generating nozzle and distance H between the nozzle exit and stagnation point (the central point of the impinged surface). The analyses regarded the axisymmetric cases, so geometry is shown in such manner.
Applied boundary conditions were as follows: • at the inlet: the constant bulk velocity, determined from the Reynolds number values (namely Re = 5000; Re = 10,000; Re = 23,000) and the constant bulk temperature equal to 293 K, • at impinged wall: the constant heat flux equal to 1000 W/m 2 , • at the upper wall: the constant heat flux equal to 0 W/m 2 , • at the outlet: the gauge pressure equal to 0 Pa.  In Table 1, the summary of all investigated cases is presented. Variables are shown in a non-dimensional way, an approach typical in the literature regarding the jet-impingement. The inlet Reynolds number based on inlet bulk velocity u b : where D is the jet issuing part diameter, m and ν is the kinematic viscosity, m 2 /s. As the inlet diameter size might strongly impact the results of the heat transfer at the impinged surface, which was proven in [42][43][44], it was kept constant throughout the research and equal to 0.01 m. This information is very often missed in published sources, which can rise misunderstandings. Moreover without it, particular results should not be used as the data validation.
Combining all the above cases, the 54 simulations were performed. The Reynolds number equal to 23,000 was selected as a reference case, widely used as the jet impingement benchmark [14,[18][19][20][21]. Other values of the Reynolds number were adapted in the basis of [35][36][37] and conditions occurring in the minijet heat exchanger. The values of H/D and R/D were chosen arbitrarily. Nevertheless they reflect the possible geometries. For example curvatures of R/D ≤ 4 were not possible to be applied for the jet impingement cases, as the geometry would be too narrow. Also the H/D ratios higher than 4 would significantly increase the overall size, which is unwelcomed in the case of the heat exchangers.
In Figure 3, the placement and definition of all boundary conditions are also shown. At the inlet, the fully developed profile was applied, using the mapping method [45], which required prescription of the bulk velocity u b , and the bulk temperature T b . Turbulent inlet conditions were then concluded with every consecutive iteration by mapping the downstream flow back to the inlet until the flow was fully developed. The final turbulence intensity was afterwards equal to 5.8%, while the ratio between the centreline velocity and bulk velocity u c /u b was equal to 1.2. It is important to note, that such intensity was relatively high, and strongly impacted the results in the stagnation zone. As shown in [22,23,44], assumed inlet conditions (both regarding velocity profile and turbulence parameters) affect significantly the heat transfer rates at the impinged surface. Such information very often is not presented in the literature, which precludes exact representation of boundary conditions, because behind the statement fully-developed profile many possible configurations can be hidden. For the purposes of the following article, aforementioned parameters were therefore kept unchanged. In the outlet, the zero gauge pressure boundary condition was applied in all cases. At the surface, the constant heat flux q was applied, equal to 1000 W/m 2 , being more physical boundary condition than fixed temperature [23].
The working fluid was the air, in an incompressible form. All its properties were kept constant, with values listed in Table 2. The reason was mostly due to the chosen heat flux applied at the impinged surface, equal to 1000 W/m 2 . The resulted temperature gradient was relatively small, and temperature dependent parameters of the working fluid would not have impacted the accuracy of the results.

Numerical Model and Solver Configuration
OpenFOAM was chosen to perform the numerical analyses, mainly due to its opensource character and possibility to extend its functionalities [45]. To investigate the jet impingement phenomenon, it was necessary to include the conservation laws of mass, momentum and energy in the model: where u i,j are the components of velocity vector, m/s; ρ is the fluid density, kg/m 3 ; p is the pressure, Pa; µ is the fluid dynamic viscosity, Pa·s; S ij is the tensor of strain rate, 1/s; u i u j is the term for the Reynolds stress tensor, m 2 /s 2 ; Θ is the mean fluid temperature, K; θ is the fluctuation of fluid mean temperature, K; a is the thermal diffusivity, m 2 /s and u j θ is the turbulent heat flux vector, K·m/s. Turbulent component of the heat flux were obtained with the use of eddy diffusivity concept [45]. As the RANS based turbulence modelling approach was applied, averaged quantities were employed as the variables. The single jet impingement was analysed, and axisymmetric cases were considered, with symmetry axis presented in Figure 3. Single jet impingement can be thought as the phenomenon that with time-constant boundary conditions reaches steady state. Therefore, the steady-state SIMPLE algorithm was applied. Second order numerical schemes were used [45].

Selection of Turbulence Model
Based on the conclusions from [39,40,44] the ζ-f model was chosen to be used in the following paper, after being successfully implemented in the OpenFOAM software. The governing equation for the model including the turbulence kinetic energy is presented in Equation (6) and the governing equation for its dissipation is shown in Equation (7): where G is the turbulence kinetic energy production term, m 2 /s 3 , C ε1 , C ε2 , C 1 , C' 2 , σ k , σ ε , σ ζ are the variables specific for this particular model, τ is the model time scale and L is its length scale, which are used to limit the model and enhance the quality of the results obtained with it. Being a low-Reynolds type model, it avoids any wall-functions that might perturb the performance. Its main feature, in comparison to the precursor, the v 2 − f model [17], has better robustness, achieved due to reformulation of the elliptic relaxation function f, defined by Equation (9), especially by modifying its wall-boundary value and assuring the better reproduction of the turbulence kinetic energy production in the near wall region. This region is especially important when investigating the jet impingement, because there the heat transfer strongly depends on the near-wall turbulence. Finally, instead of using the v 2 as the variable, which is the wall-normal velocity scale, ζ is introduced, defined by Equations (5) and (8), and obtained by normalization of the velocity scale with turbulence kinetic energy k.
As it was mentioned in the Motivation section, the numerical model was validated [39,40] and exemplary results are shown in Figure 4, presenting the comparison of the local heat transfer results for various RANS turbulence models and experimental data obtained for the flat surface impinging jet. The green lines in Figure 4 indicate the results obtained with ζ-f model [19], their difference lies in the formulation of the turbulence intensity at the inlet [44] and indicates the borderlines of the Nusselt number results that can be obtained with this model. The high accuracy of the ζ-f model, in comparison with the reference data, can be seen.
where G is the turbulence kinetic energy production term, m 2 /s 3 , Cε1, Cε2, C1, C'2, σk, σε, σζ are the variables specific for this particular model, τ is the model time scale and L is its length scale, which are used to limit the model and enhance the quality of the results obtained with it. Being a low-Reynolds type model, it avoids any wall-functions that might perturb the performance. Its main feature, in comparison to the precursor, the 2 v f − model [17], has better robustness, achieved due to reformulation of the elliptic relaxation function f, defined by Equation (9), especially by modifying its wall-boundary value and assuring the better reproduction of the turbulence kinetic energy production in the near wall region. This region is especially important when investigating the jet impingement, because there the heat transfer strongly depends on the near-wall turbulence. Finally, instead of using the 2 v as the variable, which is the wall-normal velocity scale, ζ is introduced, defined by Equations (5) and (8), and obtained by normalization of the velocity scale with turbulence kinetic energy k.
As it was mentioned in the Motivation section, the numerical model was validated [39,40] and exemplary results are shown in Figure 4, presenting the comparison of the local heat transfer results for various RANS turbulence models and experimental data obtained for the flat surface impinging jet. The green lines in Figure 4 indicate the results obtained with ζ-f model [19], their difference lies in the formulation of the turbulence intensity at the inlet [44] and indicates the borderlines of the Nusselt number results that can be obtained with this model. The high accuracy of the ζ-f model, in comparison with the reference data, can be seen.
. In Figure 5 the exemplary results of the mesh independence tests are presented, for two of the cases mentioned in Table 1, represented by Nusselt number values. The mesh was generated in the same manner as in [39,40,44]. The most important was to maintain dense enough mesh near the heat transfer wall, as well as near the nozzle wall. While typically it is recommended for the low-Reynolds RANS turbulence models to keep the Nu In Figure 5 the exemplary results of the mesh independence tests are presented, for two of the cases mentioned in Table 1, represented by Nusselt number values. The mesh was generated in the same manner as in [39,40,44]. The most important was to maintain dense enough mesh near the heat transfer wall, as well as near the nozzle wall. While typically it is recommended for the low-Reynolds RANS turbulence models to keep the non-dimensional wall normal distance y + at a value of approximately 1, for the purposes of this research that value was below 0.5, with decent number of mesh elements in the viscous sublayer. As can be seen, despite the mesh size, obtained results remained almost the same, as long as the indicated requirements were maintained. It is worth mentioning, that the designed mesh which fulfilled the criteria for the highest Reynolds number cases was also used in less demanding cases, where Reynolds number had lower values. For all analysed cases, the mesh consisted only of hexahedral elements to minimise the impact of non-orthogonality and skewness on the results. Also the boundary layer mesh was optimised to limit the aspect ratio of its elements to avoid the results distortion.
non-dimensional wall normal distance y + at a value of approximately 1, for the purposes of this research that value was below 0.5, with decent number of mesh elements in the viscous sublayer. As can be seen, despite the mesh size, obtained results remained almost the same, as long as the indicated requirements were maintained. It is worth mentioning, that the designed mesh which fulfilled the criteria for the highest Reynolds number cases was also used in less demanding cases, where Reynolds number had lower values.
For all analysed cases, the mesh consisted only of hexahedral elements to minimise the impact of non-orthogonality and skewness on the results. Also the boundary layer mesh was optimised to limit the aspect ratio of its elements to avoid the results distortion.

Results
Prior to the results description, an information regarding the way of their presentation is necessary to be provided. The location of some parameters will be shown in a nondimensional way, by using the distance x/D, as indicated in Figure 6, where x is the length of the arc between the stagnation point and any selected point on the surface and D is the nozzle diameter. It may be argued, that instead of the arc length, the relative length of the straight segment connecting the points, also shown in Figure 6, should be used. However, for the range of the distances analysed in the following article, the difference between those two, was negligible. Due to the shape of the surfaces, difference in the distance between them and the nozzle exit occurs (particular H/D value is maintained exactly constant only at the symmetry axis). However, the discrepancies between these distances in the whole stagnation

Results
Prior to the results description, an information regarding the way of their presentation is necessary to be provided. The location of some parameters will be shown in a nondimensional way, by using the distance x/D, as indicated in Figure 6, where x is the length of the arc between the stagnation point and any selected point on the surface and D is the nozzle diameter. It may be argued, that instead of the arc length, the relative length of the straight segment connecting the points, also shown in Figure 6, should be used. However, for the range of the distances analysed in the following article, the difference between those two, was negligible.
non-dimensional wall normal distance y + at a value of approximately 1, for the purposes of this research that value was below 0.5, with decent number of mesh elements in the viscous sublayer. As can be seen, despite the mesh size, obtained results remained almost the same, as long as the indicated requirements were maintained. It is worth mentioning, that the designed mesh which fulfilled the criteria for the highest Reynolds number cases was also used in less demanding cases, where Reynolds number had lower values.
For all analysed cases, the mesh consisted only of hexahedral elements to minimise the impact of non-orthogonality and skewness on the results. Also the boundary layer mesh was optimised to limit the aspect ratio of its elements to avoid the results distortion.

Results
Prior to the results description, an information regarding the way of their presentation is necessary to be provided. The location of some parameters will be shown in a nondimensional way, by using the distance x/D, as indicated in Figure 6, where x is the length of the arc between the stagnation point and any selected point on the surface and D is the nozzle diameter. It may be argued, that instead of the arc length, the relative length of the straight segment connecting the points, also shown in Figure 6, should be used. However, for the range of the distances analysed in the following article, the difference between those two, was negligible. Due to the shape of the surfaces, difference in the distance between them and the nozzle exit occurs (particular H/D value is maintained exactly constant only at the symmetry axis). However, the discrepancies between these distances in the whole stagnation Due to the shape of the surfaces, difference in the distance between them and the nozzle exit occurs (particular H/D value is maintained exactly constant only at the symmetry axis). However, the discrepancies between these distances in the whole stagnation zone are very small. In addition, it makes the comparison of results between flat and non-flat surface jet impingement possible.
Lastly, the h variable is also shown in Figure 6, which would be used in Section 4.3 to provide results in particular, various non-dimensioned distance h/H from the impinged surface.
The Results section is divided into three main parts: • overall thermal results of the jet impingement for all 54 cases, • influence of the impinged surface type on the flow behaviour near it, • impact of the impinged surface type on the impinging jet.

Thermal Performance of Convex and Concave Cases
In Figures 7-9 the local values of the Nusselt number at the impinged surface for all 54 cases are plotted. Results related with the convex geometry are indicated by the dashed line and these related to the concave with the solid line. First conclusion that can be stated after analysis of the charts from Figures 7-9 is related to the difference in the values of Nusselt number in the stagnation point and close to it. For all convex cases, those values were noticeably higher than for concave ones. As seen in Figures 7-9, the reason lies in the stagnation zone, where the visible differences of Nusselt number values occur. Therefore, it is the area that has to be analysed in detail to explain the nature of such phenomenon.
With increasing H/D ratio, the characteristic feature of the local Nusselt number curve for some jet impingement cases is vanishing-the secondary peak. Its origin was described in detail in [46]. They related the skewness of the pressure temporal distribution with the rebound of primary vortices observed for a single vortex ring operating near the location of the secondary peak. They have stated, that the rebound is the event causing the mean heat transfer enhancement to stop. While the RANS based turbulence models do not allow the extraction of so detailed data, still their ability to predict the secondary peak on the Nusselt number curve indicates their suitability for trustworthy jet impingement analyses-which was another reason to choose the ζ-f model for the presented research [23]. This ability in reference to particular RANS models was also reviewed in [2].
zone are very small. In addition, it makes the comparison of results between flat and nonflat surface jet impingement possible.
Lastly, the h variable is also shown in Figure 6, which would be used in Section 4.3 to provide results in particular, various non-dimensioned distance h/H from the impinged surface.
The Results section is divided into three main parts: • overall thermal results of the jet impingement for all 54 cases, • influence of the impinged surface type on the flow behaviour near it, • impact of the impinged surface type on the impinging jet.

Thermal Performance of Convex and Concave Cases
In Figures 7-9 the local values of the Nusselt number at the impinged surface for all 54 cases are plotted. Results related with the convex geometry are indicated by the dashed line and these related to the concave with the solid line. First conclusion that can be stated after analysis of the charts from Figures 7-9 is related to the difference in the values of Nusselt number in the stagnation point and close to it. For all convex cases, those values were noticeably higher than for concave ones. As seen in Figures 7-9, the reason lies in the stagnation zone, where the visible differences of Nusselt number values occur. Therefore, it is the area that has to be analysed in detail to explain the nature of such phenomenon.
With increasing H/D ratio, the characteristic feature of the local Nusselt number curve for some jet impingement cases is vanishing-the secondary peak. Its origin was described in detail in [46]. They related the skewness of the pressure temporal distribution with the rebound of primary vortices observed for a single vortex ring operating near the location of the secondary peak. They have stated, that the rebound is the event causing the mean heat transfer enhancement to stop. While the RANS based turbulence models do not allow the extraction of so detailed data, still their ability to predict the secondary peak on the Nusselt number curve indicates their suitability for trustworthy jet impingement analyses-which was another reason to choose the ζ-f model for the presented research [23]. This ability in reference to particular RANS models was also reviewed in [2].    In Figure 10, selected Nusselt number curves are plotted also in two-dimensional plot. Chosen geometry configurations reflects two outermost cases. As indicated in Figures 7-9, differences between Nusselt number occurs due to the both hydrodynamic and geometric variations. In Figure 10, selected Nusselt number curves are plotted also in two-dimensional plot. Chosen geometry configurations reflects two outermost cases. As indicated in Figures 7-9, differences between Nusselt number occurs due to the both hydrodynamic and geometric variations.
.   In general, the shape of all the curves, with mentioned secondary peak less visible with increasing H/D, is similar to findings from experimental research of Lee and Lee [47], and O'Donovan and Murray [48]. While they focused on the flat surface jet impingement, similarities between their results and the ones presented here can be observed. When it comes to shape of particular Nusselt number curves in relation to specific H/D ratios and chosen Reynolds number values-it is another suggestion, that the phenomena behind the jet impingement on flat and curved surfaces share the same origin and the turbulence model adequate for flat cases would also perform well for non-flat types analyses. An influence of the inlet Reynolds number on the Nusselt number is also visible in Figures 7-9. For lower value of H/D, the Reynolds number required for the secondary peak to occur was also lower. This statement is true for both, convex and concave geometries. Nevertheless, the secondary peak, if occurs, is more visible for convex geometries. This secondary maximum of Nusselt number tends to move away from the stagnation point with increasing Reynolds number.
In Figure 11 the stagnation point Nusselt number value are presented, while in Figure 12 the Nusselt number values averaged over the surface, up to x/D = 2, are presented. This distance was chosen, as further away from the stagnation zone the stream start to change the behaviour to the typical wall jet, independent from the jet impingement phenomenon [18,19,22,25,26]. Generally, for the convex geometry with the increased curvature ratio R/D the stagnation point Nusselt number decreases, while for the concave geometry the behaviour is opposite. Such results are in agreement with the preliminary research presented in [41]. The source of such behaviour is related to the size and characteristics of the stagnation zone, varying between convex and concave impinged geometries. This is more visible for local values from Figure 11. Averaged values, shown in Figure 12, tend to differ less. It can be stated, that the averaged values remain more or less similar. Such tendency, however, occurs not for all cases and less clearly-the reason is due to averaging, which might diminish the differences. Similar situation happened, when analysing the area averaged Nusselt number values of jets array impingement in [41]. Despite the impinged surface type, discrepancies were small. Nevertheless, from Figures 11 and 12 it can be concluded, that jet impingement in the convex geometry environment assures slightly higher heat transfer rates. ary maximum of Nusselt number tends to move away from the stagnation point with increasing Reynolds number.
In Figure 11 the stagnation point Nusselt number value are presented, while in Figure 12 the Nusselt number values averaged over the surface, up to x/D = 2, are presented. This distance was chosen, as further away from the stagnation zone the stream start to change the behaviour to the typical wall jet, independent from the jet impingement phenomenon [18,19,22,25,26]. Generally, for the convex geometry with the increased curvature ratio R/D the stagnation point Nusselt number decreases, while for the concave geometry the behaviour is opposite. Such results are in agreement with the preliminary research presented in [41]. The source of such behaviour is related to the size and characteristics of the stagnation zone, varying between convex and concave impinged geometries. This is more visible for local values from Figure 11. Averaged values, shown in Figure 12, tend to differ less. It can be stated, that the averaged values remain more or less similar. Such tendency, however, occurs not for all cases and less clearly-the reason is due to averaging, which might diminish the differences. Similar situation happened, when analysing the area averaged Nusselt number values of jets array impingement in [41]. Despite the impinged surface type, discrepancies were small. Nevertheless, from Figures 11 and  12 it can be concluded, that jet impingement in the convex geometry environment assures slightly higher heat transfer rates. In Table 3, the location x/D of the secondary Nusselt number maximum is presented, for the cases in which its occurrence was the most distinguishable. The differences of about 5% occurred between convex and concave geometries. For other cases the maxima In Table 3, the location x/D of the secondary Nusselt number maximum is presented, for the cases in which its occurrence was the most distinguishable. The differences of about 5% occurred between convex and concave geometries. For other cases the maxima were less visible therefore they were omitted in Table 3.

Influence of Geometry Type on Stagnation Zone Flow Behaviour
In Figure 13, the illustrative close-up of the stagnation zones for selected geometry case is presented, for both convex and concave surfaces, for H/D = 2 and R/D = 4. Reynolds number was equal to 23,000. The zoom includes also small part of the nozzle. Analysis of this area, where relatively low velocity occurred, confirms already stated conclusions, that the stagnation zone is slightly larger in concave geometry. The reason is related strictly to the geometry and the flow distribution. In convex geometry it is easier for the stream to leave the impinging area and move further. In concave geometry, due to its shape, the outflow is being slightly limited by the shape of curvature. Stagnation zone is also characterised by the high-pressure gradient, a result of the velocity decrease. In Figure 14, local values of turbulence kinetic energy, k, are shown, for both concave and convex geometries, with an example case, where H/D = 2 and R/D = 4. The values were obtained at the first mesh nodes near the heat transferring wall. The reason for choosing these particular values was shown in the authors' publication [49]. As Nusselt number distribution is influenced mostly by the phenomena in the near-wall region, it should be referred to corresponding values of other variables from this zone, despite the fact, that for example global maxima of turbulence kinetic energy appear significantly above that region.
As can be seen in Figure 14, with increasing Reynolds number, the second peak in the turbulence kinetic energy distribution becomes more evident, and eventually higher, than the first one. Values of turbulence kinetic energy for convex cases are higher until reaching the first peak, afterwards they get lower than for concave cases. After reaching the first minima, they again exceed the values of concave cases. Eventually, after the second peak, they reach lower values. The distributions for convex and concave geometries intersperse.
The location of both inflection points is strongly related to the location of peaks in the In Figure 14, local values of turbulence kinetic energy, k, are shown, for both concave and convex geometries, with an example case, where H/D = 2 and R/D = 4. The values were obtained at the first mesh nodes near the heat transferring wall. The reason for choosing these particular values was shown in the authors' publication [49]. As Nusselt number distribution is influenced mostly by the phenomena in the near-wall region, it should be referred to corresponding values of other variables from this zone, despite the fact, that for example global maxima of turbulence kinetic energy appear significantly above that region.
As can be seen in Figure 14, with increasing Reynolds number, the second peak in the turbulence kinetic energy distribution becomes more evident, and eventually higher, than the first one. Values of turbulence kinetic energy for convex cases are higher until reaching the first peak, afterwards they get lower than for concave cases. After reaching the first minima, they again exceed the values of concave cases. Eventually, after the second peak, they reach lower values. The distributions for convex and concave geometries intersperse.
The location of both inflection points is strongly related to the location of peaks in the Nusselt number values of the same cases, as indicated in Figure 14c by the green lines. For all concave cases, those peak points are moved slightly away from the stagnation points. Away from the stagnation zone, since the distance of about x/D ≥ 2.5, the wall jet creation can be noticed from the plots of Figure 14. The relaminarization processes can be assumed there, as the turbulence kinetic energy is decreasing and reaching relatively small values. In Figure 15, the distribution of turbulence kinetic energy dissipation, ε, is compared with the distribution of local enstrophy e. It is described by Equation (10), on the basis of vorticity (velocity curl), with its unit being 1/s 2 [45]: Both variables, e and ε are shown for both concave and convex geometries, using an example case, where H/D = 2, R/D = 4 and Reynolds number was equal to 23,000. One noticeable difference is related with the values in the stagnation zone, where the turbulence dissipation ε reached relatively high values, which is related to the rapid flow suppression. Local enstrophy, on the other hand, related to the change in velocity field, according to its definition from Equation (10), does not reach high values in the stagnation region.
The visible connection between those two variables can be concluded from the Figure  15. It confirms the findings of Zhu and Antonia [50]. Away from the stagnation point regions of high local enstrophy are associated with those of high turbulence kinetic energy dissipation. It is a valuable information, especially, if the analyses of the jets array impingement would be performed, since enstrophy can possibly be used for the optimisation In Figure 15, the distribution of turbulence kinetic energy dissipation, ε, is compared with the distribution of local enstrophy e. It is described by Equation (10), on the basis of vorticity (velocity curl), with its unit being 1/s 2 [45]: Both variables, e and ε are shown for both concave and convex geometries, using an example case, where H/D = 2, R/D = 4 and Reynolds number was equal to 23,000. One noticeable difference is related with the values in the stagnation zone, where the turbulence dissipation ε reached relatively high values, which is related to the rapid flow suppression. Local enstrophy, on the other hand, related to the change in velocity field, according to its definition from Equation (10), does not reach high values in the stagnation region.
The visible connection between those two variables can be concluded from the Figure 15. It confirms the findings of Zhu and Antonia [50]. Away from the stagna-tion point regions of high local enstrophy are associated with those of high turbulence kinetic energy dissipation. It is a valuable information, especially, if the analyses of the jets array impingement would be performed, since enstrophy can possibly be used for the optimisation process. Moreover, both variables reach high values in near-wall region, indicating the importance of analyses of that region and its impact on thermal performance of various jet impingement cases. In Figure 16, the plots showing distribution of the turbulence kinetic energy dissipation, ε, and distribution of the local enstrophy, e, are presented. Both values are again shown for concave as well as convex geometries, with an example cases, where H/D = 2, R/D = 4 and Reynolds number was equal to 23,000. The values were obtained at the first mesh nodes near the heat transfer wall. The curves related to the dissipation, ε, share the same tendency of maxima location and shape, as the results of turbulence kinetic energy, shown in Figure 14c. However, analysis of the enstrophy curve leads to the conclusion, that slight difference between it and the dissipation can be seen. The first peaks of both concave and convex cases, noticeably higher, than the second, can be related to the end of the stagnation zone and start of the wall jet formulation. Even more distinguished change of location occurs for the second maxima. This time, however, they occur before the peak of turbulence kinetic energy, and appear close to the location of the secondary peak of Nusselt number. Therefore, while qualitative comparison of images from Figure 15 suggest close relation between the distribution of ε and e, plots shown in Figure 16 reveal more detailed comparison and some discrepancies. Finally, in Figure 17 the values of the near the impinged surface is different with different surface types-for convex geometry they are closer to the wall than for concave one, which is in agreement with the conclusions coming from Figure 13. In Figure 16, the plots showing distribution of the turbulence kinetic energy dissipation, ε, and distribution of the local enstrophy, e, are presented. Both values are again shown for concave as well as convex geometries, with an example cases, where H/D = 2, R/D = 4 and Reynolds number was equal to 23,000. The values were obtained at the first mesh nodes near the heat transfer wall. The curves related to the dissipation, ε, share the same tendency of maxima location and shape, as the results of turbulence kinetic energy, shown in Figure 14c. However, analysis of the enstrophy curve leads to the conclusion, that slight difference between it and the dissipation can be seen. The first peaks of both concave and convex cases, noticeably higher, than the second, can be related to the end of the stagnation zone and start of the wall jet formulation. Even more distinguished change of location occurs for the second maxima. This time, however, they occur before the peak of turbulence kinetic energy, and appear close to the location of the secondary peak of Nusselt number. Therefore, while qualitative comparison of images from Figure 15 suggest close relation between the distribution of ε and e, plots shown in Figure 16 reveal more detailed comparison and some discrepancies. Finally, in Figure 17 the values of thev 2 , calculated with Equation (5), are shown. They were obtained along the centreline of the impinging jet and h/D = 0 corresponds then with the stagnation point. The "wall-normal" velocity scale v 2 is a good indicator of the turbulence transported with the impinging jet. With higher values of the Reynolds number, the jet core is preserved longer, which is indicated by more evident curvatures close to the wall. The location of the v 2 maxima near the impinged surface is different with different surface types-for convex geometry they are closer to the wall than for concave one, which is in agreement with the conclusions coming from Figure 13.

Impact of Surface Shape on Impinging Jet Hydrodynamics
The source of the heat transfer intensification away from the stagnation point is attributed to the generation of the boundary layer due to the wall jet formulation and its interference with the big eddies, which are generated at a distance x/D corresponding with the secondary Nusselt number maxima. It is related with the inlet conditions, as shown in

Impact of Surface Shape on Impinging Jet Hydrodynamics
The source of the heat transfer intensification away from the stagnation point is attributed to the generation of the boundary layer due to the wall jet formulation and its interference with the big eddies, which are generated at a distance x/D corresponding with the secondary Nusselt number maxima. It is related with the inlet conditions, as shown in

Impact of Surface Shape on Impinging Jet Hydrodynamics
The source of the heat transfer intensification away from the stagnation point is attributed to the generation of the boundary layer due to the wall jet formulation and its interference with the big eddies, which are generated at a distance x/D corresponding with the secondary Nusselt number maxima. It is related with the inlet conditions, as shown in Figures 7-9. But even stronger influence of those occurs with the heat transfer intensification at the stagnation point and near it. As shown previously, the highest values of Nusselt number are obtained exactly there. However, as indicated in Figure 13, the stagnation zone is an area with negligible flow and high-pressure gradients. Thus the energy transfer occurs there, caused directly by the dissipation of eddies, which are transported with the inflow. It is assumed [5,50], that they are the main reason of such high heat transfer. In Figure 18, the velocity profiles in particular distance values from the heated and impinged surface, are shown, for the case, in which Re = 5000, H/D = 1 and R/D = 4. As it can be seen, while close to the nozzle exit conditions remain the same, the jet core (in which the velocity is stable along the jet normal direction) downwards starts to differ between the cases. The black arrow placed near the plots (Figures 18-23) shows the flow direction, as the plots are oriented according to it. The ratio h/H indicates the particular distance above the impinged surface.  Figure 13, the stagnation zone is an area with negligible flow and high-pressure gradients. Thus the energy transfer occurs there, caused directly by the dissipation of eddies, which are transported with the inflow. It is assumed [5,50], that they are the main reason of such high heat transfer. In Figure 18, the velocity profiles in particular distance values from the heated and impinged surface, are shown, for the case, in which Re = 5000, H/D = 1 and R/D = 4. As it can be seen, while close to the nozzle exit conditions remain the same, the jet core (in which the velocity is stable along the jet normal direction) downwards starts to differ between the cases. The black arrow placed near the plots (Figures 18-23) shows the flow direction, as the plots are oriented according to it. The ratio h/H indicates the particular distance above the impinged surface. In Figure 19, the values of turbulence kinetic energy, corresponding with the locations and cases from Figure 18, are shown. The maxima, visible in the first 3 plots, indicate the regions, in which the shear occurs due to the jet hitting initially still fluid. It can be noticed, that the values of turbulence kinetic energy at the stagnation zone and close to it increase on consecutive figures-which proves the statement about the energy transfer due to the dissipation of eddies, which are transported with the inflow. In Figure 20, similarly to the situation from Figure 19, the values of enstrophy flux are presented. As for the case of turbulence kinetic energy, the shear regions are evident. However, when the jet approach the heated surface, vorticity values are low and increase while moving away from the stagnation point. It can be attributed with the already mentioned fact that almost no direct flow occurs in the stagnation zone. In Figure 19, the values of turbulence kinetic energy, corresponding with the locations and cases from Figure 18, are shown. The maxima, visible in the first 3 plots, indicate the regions, in which the shear occurs due to the jet hitting initially still fluid. It can be noticed, that the values of turbulence kinetic energy at the stagnation zone and close to it increase on consecutive figures-which proves the statement about the energy transfer due to the dissipation of eddies, which are transported with the inflow. In Figure 20, similarly to the situation from Figure 19, the values of enstrophy flux are presented. As for the case of turbulence kinetic energy, the shear regions are evident. However, when the jet approach the heated surface, vorticity values are low and increase while moving away from the stagnation point. It can be attributed with the already mentioned fact that almost no direct flow occurs in the stagnation zone.          As can be seen, the differences between the cases are less visible now. The reason is due to relatively high R/D ratio. As proven by authors in [34] such convex and concave geometries tend to differ less from the flat case impingement, and from themselves. Analysis of the velocity profiles indicates, that the jet core is preserved for a longer period of the jet-it leads then to the heat transfer intensification in the stagnation zone, compared to cases where Reynolds number was equal to 5000. The black arrows also show the flow direction.

Conclusions
In the article, the thermal and hydrodynamic results of non-flat surface single air jet impingement were presented. They were obtained with the use of tRANS-based, turbulence model of Hanjalic, ζ-f, [23], which in the authors' opinion is capable of predicting properly the phenomena occurring in the analysed cases. The ζ-f turbulence model performs well in the stagnation region, where many RANS-based turbulence models tend to overpredict the generation of the turbulence kinetic energy [2]. According to the author's knowledge, such cases were analysed for the first time with that model. It was implemented and validated in the open-source software-OpenFOAM. The model's main features, which are relative robustness and convergence and more importantly the ability to indicate the secondary peak in the local Nusselt number distribution at the heated and impinged surface made it possible to calculate 54 various cases, with variable non-dimensional geometry constraints: R/D and H/D, as well as the inlet Reynolds number.
Noticeable differences between convex and concave surfaces impingement were identified, regarding not only stagnation point or, more generally, stagnation zone thermal performance, but also the existence of the secondary peak in relation to the mentioned variables. The shape of particular geometry impacts the size of stagnation zone, which results in the flow characteristics when the wall jet starts to formulate after impingement. As can be seen, the differences between the cases are less visible now. The reason is due to relatively high R/D ratio. As proven by authors in [34] such convex and concave geometries tend to differ less from the flat case impingement, and from themselves. Analysis of the velocity profiles indicates, that the jet core is preserved for a longer period of the jet-it leads then to the heat transfer intensification in the stagnation zone, compared to cases where Reynolds number was equal to 5000. The black arrows also show the flow direction.

Conclusions
In the article, the thermal and hydrodynamic results of non-flat surface single air jet impingement were presented. They were obtained with the use of tRANS-based, turbulence model of Hanjalic, ζ-f, [23], which in the authors' opinion is capable of predicting properly the phenomena occurring in the analysed cases. The ζ-f turbulence model performs well in the stagnation region, where many RANS-based turbulence models tend to overpredict the generation of the turbulence kinetic energy [2]. According to the author's knowledge, such cases were analysed for the first time with that model. It was implemented and validated in the open-source software-OpenFOAM. The model's main features, which are relative robustness and convergence and more importantly the ability to indicate the secondary peak in the local Nusselt number distribution at the heated and impinged surface made it possible to calculate 54 various cases, with variable non-dimensional geometry constraints: R/D and H/D, as well as the inlet Reynolds number.
Noticeable differences between convex and concave surfaces impingement were identified, regarding not only stagnation point or, more generally, stagnation zone thermal performance, but also the existence of the secondary peak in relation to the mentioned variables. The shape of particular geometry impacts the size of stagnation zone, which results in the flow characteristics when the wall jet starts to formulate after impingement. Distribution of local Nusselt number was also related to the turbulence kinetic energy, k, values near the heated wall, where differences, especially connected with the Reynolds number at the inlet, were shown.
In addition, the trial to interconnect the rate of turbulence dissipation, ε, and enstrophy, e, in the analysed models were done and presented. Enstrophy, related to the vorticity and velocity distribution can be taken as a measure of flow energy dissipation. While the big-scale comparison of ε and e showed similarities, the near-wall comparison indicated differences, related mostly with the impact of wall jet formulation on the enstrophy field.
It is worth to mention, that while enstrophy can provide a meaningful added value, it is not a common practice to use it when simulating jet impingement with the use of RANS models. Authors did not find any, related to the results above, examples of such approach. It is assumed, that this variable might be even more useful, when applied to the analyses of jets array impingement and used to identify the zones of particular thermal behaviour.
To sum up, only selected results were presented. The number of investigated cases made it impossible to include all possible variations. Nevertheless, the presented cases indicate the most important outcome from the research-altogether illustrating the most important differences between the convex and concave surfaces impingement.
The presented results will be of great importance after the next step of the research project-investigation of the jets arrays impingement-which is already started, to be able to optimise the occurring flow patterns.