Heat transfer intensification by jet impingement – numerical analysis using RANS approach

. Jet impingement is a method of the heat transfer enhancement applied in the engineering systems. The idea is to generate fast-flow fluid jet which impinge on the heated (or cooled) surface, causing significantly higher heat transfer rate. Although some flat surface jet impingement cases are described in the literature, the validated data is still limited. The reason is coming from the fact, that these flows are hydrodynamically complex. Therefore the numerical analysis is necessary to understand the phenomena, especially in the range of turbulent flow. The well-known and accurate method is Reynolds Averaged Navier-Stokes (RANS) approach. However, depending on the applied turbulence model, various results can be obtained. The reason is that the jet impingement strongly depends on the complex boundary layer effects and their resolving is still challenging for RANS models and until now it is their weakest point. In the paper, the hydrodynamic and thermal, numerical results of jet impingement are presented, depending on selected RANS based models. The aim was to indicate their ability to anticipate the turbulence parameters.


Introduction
The topics of pollution and energy waste, related to energy generation and utilization are ones of the biggest scientific issues of modern days. Novel methods of energy generation are being investigated, therefore the ideas, such as the fuel cells or hydrogen based cars appear more and more common. Moreover, old-fashioned facilities are being refurbished or replaced. This approach, to convert already existing object into modern ones is even more important. The new techniques, which can reduce energy consumption and/or increase the efficiency are searched for.
Together with many other ideas of how to improve the performance of existing devices, the novel concept of heat exchanger was proposed at the Gdansk University of Technology [1][2]. It offers interesting combination of increased heat transfer rates between two media and almost negligible pressure losses. The innovation is connected with the implementation of partitions around the cylindrical heat exchanging wall, having many small orifices drilled through it. As a consequence, the orifices cause generation of the jets, which are influenced by the crossflow. When jets hit the surface, they cause increase of the transferred heat flux. Experimental investigations are the source of data regarding the macroscale phenomena and device performance. However, there is lack of their understanding at microscale. Turbulence generated by the jets near the wall is expected to be the main reason of heat transfer intensification. Therefore the proof is needed, mainly to get the ability of further optimization of the original design.
The fluid mechanics, gas dynamics, materials durability or heat and mass transfer processes -they all are very demanding for scientists. Amongst them however, the most problematic is the issue of turbulence. Three approaches are present in the world of science: analytical, experimental and numerical. For the case of turbulence, there is no universal analytical solution of governing equations (Navier-Stokes), experimental facilities are problematic to be constructed, due to existence of many scales characterising turbulence, and, lastly, the numerical approaches exist in hundreds of different types, each one able to simulate just particular case. Nevertheless, in the Authors opinion, the last approach is the most optimistic for complex systems -it is also a subject of the biggest concern of the scientific world within fluid dynamics. The reason is rather simple: each year an increasing power of the computers allows the users to perform more and more demanding simulations in reasonable time, what is described by the Moore's law [3], coming from the observation. With an access to computer clusters it is possible to obtain detailed results faster, than with the experimental approach.
Three different types of turbulence numerical analyses are commonly used today [4]:  Reynolds Averaged Navier-Stokes (RANS),  Large Eddy Simulations (LES),  Direct Numerical Simulation (DNS). First method is the least demanding from the computational power point of view, while the last one is the most demanding and still almost impossible to be used for typical applications. LES methods can be stated as the compromise between two other approaches. The difference between them is mostly due to the resolved scale of turbulence. As a result, both RANS and LES methods include some validated analytical formulas.
For the purposes of detailed analysis of turbulence phenomenon, occurring in abovementioned heat exchanger the Authors have chosen the RANS approach. According to suggestions of Zuckerman [5], the 2 v f  [6] family of models was considered, as those, which are expected to provide the best results, when analysing jet impingement. The idea behind the research was to start with analyses of single jet in various configurations to gather enough data to model whole orifices' system taken from the heat exchanger. Some results concerning this topic were already published in [7]. In the following article, the focus was placed on description of the advantages and disadvantages of particular RANS methods. In parallel, an attention was placed on the most important parameters characterising the phenomena occurring in the boundary layer and close to it as the result of surface jet impingement.

Mathematical model
The system of equations describing phenomenon of single phase jet impingement consists of conservation of mass (Equation 1), momentum (Equation 2) and energy (Equation 3). Single impinging jet of the constant velocity tends to reach steady-state, therefore the steadystate simulations, with the use of SIMPLE algorithm, were performed. In the RANS approach, variables such as for example velocity or temperature are split into two terms, constant mean value, and its fluctuation. Terms including the fluctuations are placed in the above system of equations, and to close it, the turbulence model equations have to be included. The analysed fluid was the air. Compressibility effect is negligible with small temperature differences, such as in analysed cases, therefore density was kept as a constant value: where ui,j are the velocity components, 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; ' ' i j u u is the Reynolds stress term, m 2 /s 2 ; Θ is the mean temperature, K;  is the temperature fluctuation, K; α is the thermal diffusivity, m 2 /s and ' j u  is the turbulent heat flux, K•m/s. The symbol "¯" represents averaged quantity.

Numerical model
Analyses were managed with OpenFOAM, the powerful open-source software, developed to perform the finite volume operations on the system of differential equations [8]. For single jet, it was possible to perform the stationary analysis. Also, instead of full 3D geometry, the axisymmetric cases could be prepared. Both assumptions did not impacted quality of the solution. The second-order, upwind discretization schemes were used for the advection terms, to obtain the high accuracy of the results.
In Figure 1, the schematic view of the geometrically important parameters, as well as the placement of boundary conditions is presented. Indicated parameters are: D which is the orifice diameter, H, which is the height between impinged wall and orifice exit and x, being the distance of particular point from the stagnation point. At the inlet the constant bulk velocity was kept, with the turbulent flow profile maintained along the orifice diameter. Neumann boundary condition (zero gradient of the pressure) was set at the outlet. In Figure 2(a), the close-up of the mesh near wall can be seen. Total number of mesh elements reached up to 63 000, and the characteristic near-wall mesh parameter y+ was always kept below 0.5, which was equal to the first near-wall element height of approximately 4e-6 m. The results of the mesh independence test are shown in Figure 2(b), differences between local Nusselt number values for applied meshes were negligible. For the validation purposes, the reference case from ERCOFTAC database was chosen [9], with boundary conditions listed in Table 1. As mentioned in Section 1, the 2 v f  family of RANS models is considered to be the most suitable for the jet impingement. Its superior performance over k-ε model, especially in the stagnation zone, can be seen in Figure 3. The original version of 2 v f  model was proposed by Durbin [9], and extended for example by Kenjeres et al. [10]. Nowadays, more than 10 different types of it exist, according to Billard [11]. The Authors focused on two of them, proposed by Lien and Kalitzin in 2001 [6] and Hanjalic in 2004 [12]. The first one is available in the most common CFD software, while the second one, as it will be presented in the following article, shows noticeably better performance, but was not previously implemented. Therefore Hanjalic's model was implemented and verified in OpenFOAM by Authors for the purposes of the described research. The models are characterized by the system of four complementary equations of the following variables:  k (turbulence kinetic energy),  ε (dissipation of turbulence kinetic energy, k), According to [11], while being very accurate, the original 2 v f  model of Durbin presented highly unstable behavior. The reason was mostly due to the formulation of wall boundary condition for f: where N is the model constant; ν is the kinematic viscosity, m 2 /s; y is the distance between the wall and the first mesh node near it, m. Such models require dense mesh near the wall for a proper calculation of the boundary layer processes, therefore possible value of y can be very small. As it can be indicated in the denominator of Equation 4, fourth power of y suggests that the resulting value of the fwall might cause computational problems. Lien and Kalitzin [6] reformulated the whole model, with the main modification proposed for the problematic term: 0. wall f  (5) Their model turned out to be very robust, but it lacked accuracy of the original one. That is the reason of the Hanjalic's research, who besides another reformulation of fwall: with noticeable difference being the y second power in the denominator, changed also one variable in the equations system. Instead of 2 v , ζ was proposed, which is 2 v normalized by k. Therefore the Hanjalic's model is called ζ-f.

Results
In the following Section, the brief explanation of important parameters concerning jet impingement and their relations is presented and explained. In Figure 3, the very important feature of local heat transfer distribution, the secondary peak along the heated surface, can be seen. It occurs under particular conditions, which was confirmed also by the experiments [13], and is very difficult to obtain using the numerical models, as it can be seen through Behnia results. The reason lies in overproduction of turbulence kinetic energy, k, in the stagnation zone, which is typical for the most common RANS models, as shown in Figure  4(a) and 4(b).  [9] and Yan [13] with the results obtained by Authors (k-ε and ζ-f). Figures 5 and 6 concern different parameters, local turbulence kinetic energy k, defined in [7] and vorticity, ω, defined in [8] and obtained with the application of the ζ-f model. It is important to notice, that in both figures the values from two different probed lines are shown, one inside and one outside the boundary layer (H/D = 4e-4 and H/D = 0.1). Line inside was placed according to the location of the first mesh node near the wall, while the one outside was placed above the boundary layer height. As it turned out, only values investigated inside boundary layer, just above the heated surface are meaningful for the purposes of presented research. Both distributions of k and ω for that location can be compared with location of the Yan [13] local Nusselt number from Figure 3. Inflections of all curves occur nearly at the same x/D points, as shown by dashed lines in all figures. However, the order of maximal values is reversed. Another conclusion regards the fact that vorticity values near the wall turned out to be significantly higher, than outside it. For turbulence kinetic energy, the situation was opposite.
(a) k-ε (b) ζ-f  In Figure 7 another turbulence kinetic energy k profile is shown, obtained for both k-ε and ζ-f models -perpendicular to the heated surface, at the distance from stagnation point equal to x/D ≈ 1.78, where the secondary peak of Nusselt number occurs. The results correspond with the results presented in Figure 4(a) and 4(b). In addition, the velocity U profile is also shown. The scale for both parameters remains the same. The visible peaks of turbulence kinetic energy is located above the peak of velocity. They occur in the zone of the k maxima area, noticeable in Figures 4(a), 4(b) and 5. As presented earlier, the situation outside the boundary layer does not have the direct impact on the thermal performance of the jet impingement heat transfer. The differences between both turbulence models are visible, especially in the near-wall zone, where the characteristic skew of turbulence kinetic energy k profile can be noticed for ζ-f model. In Figure 8 the profile of vorticity ω is shown, obtained for both k-ε and ζ-f models, located at the same distance x/D ≈ 1.78 as the profile of turbulence kinetic energy k, shown in Figure 7. In contrary to the analysis of Figure 7, the profiles of vorticity and velocity exhibit similarities outside the boundary layer. Highly turbulent character of the flow in the area of the boundary layer is also visible, with noticeably high values of vorticity visible at the lowest part of the Figure 8. Discrepancies between both utilized models can be noticed.

Summary
The scope of presented research was to provide information related to the selection of proper turbulence model, when analysing the jet impingement. In addition, some suggestion were provided about the appropriate post-processing and analysis of the results -as suggested in Section 1, the phenomena causing the heat transfer intensification during the jet impingement occur mostly close to the wall. The focus should be placed on investigation of the boundary layer, as fluid behaviour outside it not necessarily impacts the intensification. Presented own results proved, that the chosen ζ-f model can be accurate enough to represent impingement heat transfer.