Heat transfer measurements and CFD simulations of an impinging jet

Heat transport in impinging jets makes a part of many experimental and numerical studies because some similarities can be identified between a pure impingement jet and industrial processes like, for example, the heat transfer at the bottom of an agitated vessel. In this paper, experimental results based on measuring the response to heat flux oscillations applied to the heat transfer surface are compared with CFD simulations. The computational cost of a LES-based approach is usually too high therefore a comparison with less computationally expensive RANS-based turbulence models is made in this paper and a possible improvement of implementing an anisotropic explicit algebraic model for the turbulent heat flux model is evaluated.


Introduction
Impinging jets are frequently used in industrial processes because high heat transfer intensities can be achieved at relatively small surfaces. A lot of experimental, numerical and analytical data concerning the hydrodynamics and heat or mass transfer in impinging jets can be found in the literature, mainly in a configuration where the fluid (liquid) jet impinges a smooth plane surface.
In CFD simulations of an impinging jet, the computational cost of LES-based or DNS-based approach is usually too high therefore less computationally expensive RANSbased (Reynolds averaged Navier-Stokes equation) turbulence models [1] are more frequently used in practise. With a RANS-based approach, the Reynolds stress tensor in the time-averaged Navier-Stokes equation must be described by some turbulence model. The simplest definition is based on the Boussinesq eddy-viscosity approximation [2] which assumes that the Reynolds (turbulent) stress tensor is proportional to mean strain-rate tensor S i j and the proportional factor is the eddy (turbulent) viscosity.
Turbulence models like, for example, Spalart-Allmaras, k− ε and k − ω are based on this approach and using a set of another (transport) equations, quantities like turbulence kinetic energy k and energy dissipation ε are described and the turbulent viscosity is then expressed as a scalar value.
In some situations when the turbulence anisotropy is important (highly swirling flows, stress induced secondary flows in rectangular ducts, impinging jets, etc), this approach may not give correct results and a second-order closure Reynolds-Stress Model [3] could be more useful.
is frequently described by the scalar value of the eddy viscosity and turbulent Prandtl number (which is frequently taken as constant 0.85 [4]).
This approach is not able to reflect the possible anisotropy in the turbulent heat flux as pointed out by [5] and also other authors. A full second-order closure transport equation for the turbulent heat flux, similar to the transport equation for the Reynolds stress tensor [6] could resolve the anisotropy of the turbulent heat flux, nevertheless with many other unknown terms which would have to be modeled somehow. Therefore, simpler explicit algebraic heat flux models can be used as proposed by Daly and Harlow [7] − and other authors [8,9]. Such models can improve results in complex cases where the turbulence anisotropy is important and they have relatively small impact on the total computational requirements.

Experimental methods
Experimental data in the literature describing the heat transfer for impinging jets on a flat plate can show relatively large scatter (20 − 30%), see Figure  data [10][11][12][13][14][15][16] . The discrepancies in the available data can lie in slightly different geometrical configurations, boundary conditions, as well as in the experimental method itself and its accuracy. Some experimental techniques are based on measuring the velocity or temperature profiles in the fluid and the quantities like the shear stress or heat transfer coefficients at the impinging surface are then indirectly evaluated using the momentum or heat transfer equations. Large group of experimental methods measures the momentum and/or heat transfer characteristics directly at the impinged surface. The shear stress (shear rate) at the impinged surface can be measured by the electrodiffusion method which is based on measuring the limiting diffusion current between two electrodes. The small working electrode is situated at the impinged surface and the size of the diffusion current is proportional to the molar flux of ion exchange and, at some conditions, it is also proportional to the heat transfer coefficient. The limiting diffusion current can be used to determine the shear stress (rate) at the location of the working electrode, see for example [17][18][19][20]. Using the analogy between mass and heat transfer, the heat transfer characteristics can be determined [21].
The heat transfer coefficients are frequently measured by the infrared-thermography methods. A stream of fluid impinging the plate heated by Joule's heat leads to a nonuniform temperature field at the plate surface which can be measured by the infrared thermography and consequently, the local heat transfer coefficients can be determined, see [16], for example. The infrared thermography was also used by Freund in his work [22] to measure the temperature response of the plate induced by an oscillating heat flux impinging the plate. This was an extension of work by [23] and it is briefly described below. See [24][25][26] for applications of this method (TOIRT -Temperature Oscillation Infra Red Thermography).
The principle of the TOIRT method by Wandelt and Roetzel [23] is depicted on Figure 2. It is based on measuring the time dependency of T | z=δ (or other corresponding characteristic) at the surface plate of width δ induced by oscillating heat flux q| z=δ . The temperature response de-pends on thermophysical properties of the plate and heat transfer coefficients on both sides of the plate, i.e. α 0 and α δ . Transient temperature field T in a homogeneous plate with heat diffusivity a is described by the Fourier equation Neglecting the lateral heat conduction described by the first two terms on the right-hand side and adding boundary conditions on both sides of the plate, it is possible to find an analytical solution for the case of a periodically oscillating impinging heat flux. The boundary conditions of the third kind can be applied on both sides, with a periodic oscillation of the heat flux on one side. They can be described as where λ is the thermal conductivity of the plate, ω is the frequency of the oscillating heat flux with amplitudeq, and α 0 , α δ are the already mentioned heat transfer coefficients on both sides of the plate. Using the Laplace transform, a temperature solution on both sides of the plate for amplitude A and phase shift ϕ can be found [23] T Then, the phase shift at the surface with induced temperature oscillations (z = δ) can be described [23] as where the following dimensionless variables and parameters are used It is maybe not clear from the complicated solution (10), but the phase shift of the temperature response depends on the size of unknown heat transfer coefficient α 0 , and, on the contrary, it does not depend on the amplitude of the heat flux and, furthermore, it does not depend on the accuracy of the oscillating temperature field measurements. In this case with the neglected lateral heat conduction, it is not necessary to measure the surface temperature precisely, only the exact time dependency between the impinging heat flux and the measured temperature response is required. Figure 3 shows an example of relation between the heat transfer coefficient α 0 and phase shift ϕ| z=δ for our experimental parameters and values mentioned in the following text. Figure 2. The principle of the TOIRT method by Wandelt and Roetzel [23]. The left side of the wall is impinged by a stream of air, the right side is heated by an oscillating heat flux and the temperature response T (t, z = δ) is then measured by an infrared camera.

Experimental results
The heat transfer coefficient in an impinging jet was measured on a circular plate with diameter 500 mm and thickness 1 mm. The plate was made of steel with heat conductivity 15 W m −1 K −1 , density 7628 kg m −3 and specific heat capacity 460 J kg −1 K −1 , that is thermal diffusivity 4.275 × 10 −6 m 2 s −1 . The impinging jet was produced by a stream of air leaving the circular jet nozzle of inner diameter D = 45.8 mm at distance H = 91.6 mm from the impinged plate, that is H/D = 2, see Figure 2. The temperature of the jet air and surroundings was 25.9 • C, and the corresponding properties are: kinematic viscosity 15.85×10 −6 m 2 s −1 , thermal conductivity 0.026 W m −1 K −1 , and the Prandtl number 0.71. Corresponding mean volumetric (bulk) velocity u b = 7.82 m s −1 used in our measurements was calculated for the Reynolds number 22596.
The oscillating heat flux impinging the other side of the plate was produced by 2 × 500 W halogen lamp. The light intensity modulation, that is the size of the impinging heat flux was maintained within range 44% -100% of the maximum value using the voltage variation. The modulation proceeded with frequency 0.1 Hz and it was controlled by a precise sine-impulse generator. The temperature response of this plate side was measured by the infrared camera thermoIMAGER TIM 160 with sampling frequency 10 Hz (camera resolution 160 x 120 points, the real distance between the measured points corresponded to approximately 1.56 mm, the maximum sampling frequency is 120 Hz which we have not used because the measured data are stored in a proprietary file format at this frequency and cannot be processed easily). As mentioned above, the approximate solution (10) neglects the lateral heat conduction. Freund et al [25] wrote that if ξ < 0.5, the lateral conduction is negligible compared to the heat transfer rate through the wall and the spatial resolution is then expected to be limited by the resolution of the camera. In our case, the value the dimensionless thickness was ξ = 0.27 therefore we were in accordance with this condition. Figure 3 presents the solution of Equation (10) for the thermophysical parameters mentioned above. The corresponding values of the Nusselt number were then expressed as where λ is the thermal conductivity of the fluid (air) and they are presented on Figure 4 in comparison with other literature data [12,14,16]. The measured data do not agree very well with the literature data. The main reasons lie in relatively small values of the heat transfer coefficients (below 10 W m −2 K −1 ) therefore the measured dependency between the phase shift ϕ and the heat transfer coefficient (see Figure 3) was at the accuracy limit of our measurements.

Numerical simulations
We have performed several numerical simulations to get a qualitative as well as quantitative comparison of several turbulence models available in ANSYS Fluent. diameter D at distance H from a plate which is heated by constant heat flux q. We have solved the case as symmetric (axisymmetric), and a fully developed velocity profile obtained in a separate computation for the given Reynolds number Re = u b D/ν and corresponding turbulence model was applied at the inlet. The pressure outlet boundary condition was used at the top and the right boundary of the solution domain. The maximum radius of the outer boundary related to the jet diameter was r/D = 10 in order to minimize the influence of reversed (back-flow) quantities. At the bottom wall, constant heat flux q = 1000 W/m 2 were used in heat transfer computations. The mesh used in our simulations satisfied condition y + ≈ 1 at the impinged wall boundary. A grid independence study was performed for three different mesh sizes (13000, 23000, 50000), and the grid convergency index [27] for the finest mesh was evaluated as 3.3% for the mean velocity field at r/D = 1.0, and 0.03% for the temperature value at the stagnation point. Figure 6 shows relatively small differences of numerical results on a full 3-D geometry (700 000 mesh elements) compared to the 2-D axisymmetric case therefore we made most of our simulations on the 2-D finest mesh with 50 000 mesh elements.    Figure 6. The Nusselt number comparison of 2-D axisymmetric case (50 000 mesh elements) and a full 3-D case (700 000 mesh elements) with experimental data [12,14]. Height of the jet discharge above the plate is H/D = 2, Reynolds number Re = 23000. The Intermittency Transition turbulence model was used here.
In our previous work [28], we have described the implementation of anisotropic models for the turbulent heat flux in ANSYS Fluent. The anisotropic turbulent heat flux was implemented through a user-defined-scalar (UDS) transport equation of scalar φ It provides a possibility to use anisotropic diffusivity Γ i j .
The following definition based on Daly-Harlow model (5) for the turbulent heat flux and φ = c P T can give us the following equation which, divided by specific heat capacity c P , corresponds to Equation (3) for Reynolds averaged temperature T (steadystate case with no internal heat sources, constant density and specific enthalpy). Formally same equation as (16) and same definition (15) for diffusivity coefficient Γ i j can be obtained for temperature as the user-defined scalar, that is φ = T . Some differences should appear only in the definition of boundary conditions, see [28] for more details. Implementing an anisotropic turbulent heat flux with some available RSM turbulence models in ANSYS Fluent unfortunately did not yield in a better agreement with the experimental data, the values of the Nusselt number were overpredicted especially in the stagnation region, see Figure 7. Daly-Harlow explicit algebraic heat flux model [7] was used here and qualitatively similar results were obtained with other models [8,9]. The following code shows an example of UDF function implementing diffusivity coefficient Γ i j in the user-defined scalar equation for the Daly-Harlow model in 2-D. The overprediction of the Nusselt values near the stagnation point could be resolved by a production limiter like Kato-Launder [29] which limits the excessive production of turbulence kinetic energy caused by a very high level of shear stress rate in the stagnation regions [4]. See Figure 8 comparing the results for SST k − ω turbulence model with and without the Kato-Launder limiter.
The Kato-Launder limiter [29] is also part of the AN-SYS Fluent implementation of the SST k −ω Intermittency Transition model which provided the best agreement with the experimental data. It is a modified version of the Transition SST model which is also known as γ − Re θ transition model [30]. The Intermittency Transition model solves only one additional transport equation for intermittency γ in addition to the kinetic turbulence energy and dissipation equations. It avoids the dependency of Re θ on the velocity [4] and makes the Intermittency Transition model Galilean invariant in contrast to γ − Re θ transition model which has another transport equation for momentum-thickness Reynolds number Re θ . The coupling of the transition model with the SST turbulence model is done through the modified production and destruction terms in the equation for turbulent kinetic energy k, see [4,[30][31][32] for more details.

Conclusions
Results of numerical simulations of a round impinging jet were compared with our own experimental data as well as with the literature experimental data. The measured heat transfer coefficients were relatively small therefore the eval- uated phase shift with the TOIRT method was at the accuracy limit of our measurements. A different fluid media, water instead of air, for example, could help to resolve these problems and get better agreement with the literature data.
We have implemented an anisotropic model for the turbulent heat flux through a user-defined-scalar equation in ANSYS Fluent to get better prediction of the heat transfer characteristics in an impinging jet. Nevertheless, the best agreement was provided with the SST k − ω Intermittency Transition model which belongs to the RANS turbulence models based on the scalar eddy viscosity (Boussinesq approximation). The second-order closure Reynolds Stress models with our implementation of the explicit algebraic model for the anisotropic turbulent heat flux were not able to get better agreement with experimental data. The location of the secondary peak in the Nusselt dependency was not correct and the values of the Nusselt number were overpredicted in the stagnation region substantially. A further investigation is necessary to identify the main reasons of these discrepancies.

Acknowledgment
This work has been supported by research project of the Grant Agency of Czech Republic No. 14-18955S.