Turbulent heat transport and its anisotropy in an impinging jet

. The turbulent heat transport is anisotropic in many cases as reported by several researchers. RANS-based turbulence models use the turbulent viscosity when expressing the turbulent heat ﬂux in the energy balance (analogy of the Reynolds stresses in the momentum balance). The turbulent (eddy) viscosity calculation comes from the Boussinesq analogy mainly and it represents just a scalar value, hence a possible anisotropy in the turbulent ﬂow ﬁeld cannot be simply transferred to the temperature ﬁeld. The computational cost of a LES-based approach can be too prohibitive in complex cases, therefore simpler explicit algebraic heat ﬂux models describing the turbulent heat ﬂux in the time-averaged energy equation could be used to get more accurate CFD results. This paper compares several turbulence models for the case of a turbulent impinging jet and deals with a methodology of implementing a user-deﬁned function describing the anisotropic turbulent heat ﬂux in a CFD code.


Introduction
Several approaches like RANS, LES or DNS can be used in numerical solutions of turbulent flows. The RANS approach is based on the Reynolds averaged Navier-Stokes equations (RANS) [1]. Substituting the instantaneous velocity as a sum of its mean and fluctuating components u = u + u to the Navier-Stokes equation and averaging with respect to time, we can get the following equation (simplified to incompressible case with no body forces) where the last term on the right-hand side represents the Reynolds stress tensor. It is a symmetric tensor with 6 independent components which must be described by some other equations in order to close Equation (1). The simplest approach to describe the Reynolds stress tensor 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.
This approach is used in turbulence models like Spalart-Allmaras, k − ε and k − ω models which describe the new unknown variables, turbulence kinetic energy k and energy dissipation ε, for example, by a set of another (transport) equations. The turbulent viscosity is then usually expressed as a scalar value. These models (RANS) and their variants perform very well in many situations but their main deficiency is that the eddy viscosity is assumed to be an isotropic scalar quantity a Corresponding author: karel.petera@fs.cvut.cz hence they are not able to give correct results in situations when the turbulence anisotropy is important (highly swirling flows, stress induced secondary flows in rectangular ducts, impinging jets, etc). A second-order closure Reynolds-Stress Model [3] can be used to solve such problems, but the problem is that the turbulent heat flux representing similar term in the time-averaged form of the energy equation as the Reynolds stress tensor in the momentum equation (for the sake of simplicity, S h referring to other internal sources and incompressible system is assumed here) is usually described using 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 LES-based approach is computationally very expensive due to large requirements on grid resolution as well as the transient basis of solution which has to be averaged over a sufficiently long flow-time interval. A full second-order closure transport equation for the turbulent heat flux, similar to the transport equation for the Reynolds stress tensor [6] could provide a possibility to describe the anisotropy of the turbulent heat flux, nevertheless with many other unknown terms which would have to be modeled somehow. The general form of this transport equation shows an explicit dependency on the Reynolds stresses u i u j which is not the case of Equation (5). If we wanted to avoid solving another transport equation, simpler explicit algebraic heat flux models can be used as pro- or Abe and Suga [8] − Such models can improve results in complex cases where the turbulence anisotropy is important with a relatively small impact on computational requirements. This paper compares some simulation results of a single round impinging jet with experimental data and gives focus on methodology of implementing a user-defined function in ANSYS Fluent solver describing the anisotropic turbulent heat flux on the basis of an explicit algebraic heat flux model.

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. Figure 1 shows our domain of interest, that is a single round jet of 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.

Flow field simulation results
First, we focused on mean velocity profiles at several dimensionless distances r/D from the axis, that is on lines parallel with the symmetry axis (see the three vertical black lines in Figure 1). Figure 2 describes the dependency of mean velocity u m = √ u 2 + v 2 on distance from the plate y/D at these lines. It is normalized by the bulk (average) velocity u b at the jet nozzle. Experimental data [9] are compared here with turbulence models based on the Boussinesq eddy-viscosity approximation (EVM), namely Realizable k − ε model with Enhanced Wall Treatment, SST k − ω model [1], Intermittency Transition Model (new in AN-SYS Fluent v15 [4]), k-kl-ω model [10], and Transition SST model [11]. One can see here that RKE and k-kl-ω models do not give very good results compared to other models. The Intermittency Transition model as well as Transition SST models produced almost same results, and the two-equation SST k − ω model is relatively close to them. Full second-order RSM closures with linear pressurestrain model, quadratic pressure strain model, and Low-Re Stress-Omega model (see [4] for complete description) are presented on Figure 3. Surprisingly, their results compared to the EVM results are much worse.
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 [12] 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. We used the finest mesh with 50 000 mesh elements in all our simulations. Figure 4 shows a comparison of experimental data [13,14] with our simulation results based on "eddy-viscosity" models. The best agreement with experimental data was obtained for the Intermittency Transition Model and Transition SST models -there is practically no difference between these models hence the Transition SST model is omitted in the figure. They are the only models of this category which are able to predict the secondary peak in Nusselt number dependency on radius r/D.

Heat transfer simulation results
Concerning the RSM models, the simulation results are not satisfactory which comes from the incorrect prediction of velocity profiles. We were not able to get better agreement with experimental data in contrast to [15]. Figure 5 shows that only the quadratic pressure-strain RSM model predicts the secondary peak in the Nusselt number however the position is wrong and the Nusselt value at the stagnation point is underpredicted.
Nevertheless, experimental data in literature [16,17,13,18,14,19] describing the heat transfer for impinging jets on a flat plate show relatively large scatter (20 − 30%), see for example Figure 6 which compares the literature data with our simulation results based on some turbulence models. It is the case with larger distance of the jet nozzle from the plate, H/D = 6, and it demonstrates here that a comparison with simulation results is not always as clear as one would wish to be.

UDF methodology and implementation
As mentioned in the introduction of this paper, the possible anisotropy in the turbulent flux is not reflected by current CFD turbulence models of ANSYS Fluent even with the RSM models. ANSYS Fluent provides a way to create user-defined functions (UDF) and scalar equations which can be used to modify or add user-defined thermophysi- cal properties, boundary profiles of some quantities, source terms in transport equations, etc. This can be used to incorporate a model like (6) or (7) in order to describe an anisotropic behavior of the turbulent heat flux and possibly get more accurate results in heat transfer simulations. The first and quick attempt to implement an anisotropic model of the turbulent heat flux could focus on redefinition of terms appearing in Equation (5), for example c P , μ t , or Pr t . But the ANSYS Fluent programming interface for these quantities is designed to return scalar values only so  it is not possible to use them in definitions of vector or tensor quantities which are necessary with equations like (6) or (7).
Another possible approach to implement a user-defined function describing the turbulent heat flux can be based on source term S h in Equation (4). Subtracting Equation (5), which represents the standard definition of the turbulent heat flux, and adding the definition based on an explicit algebraic model, for example Equation (6) following The ANSYS Fluent interface for the user defined functions provides necessary function declarations (macros) to access values of temperature gradients (C T G), turbulent viscosity (C MU T), turbulence kinetic energy (C K), turbulence energy dissipation (C D), specific enthalpy (C CP), and the Reynolds-stress components (C RUU, C RVV, . . . ). They can be used in DEFINE SOURCE macro to describe such source term. The problem with the approach based on the source term lies in the necessity to express the spatial derivatives of the term in parenthesis in Equation (8) which is not an easy task in ANSYS Fluent UDF.
A user-defined-scalar (UDS) transport equation seems to be a more straightforward approach in the ANSYS Fluent. In general, an arbitrary transport equation for scalar φ of the following form can be solved simultaneously with and Reynolds number Re = 23000, see references [16,17,13,18,14,19]. The literature data are compared with our simulation results for some turbulence models.
the basic momentum, energy and mass transport equations.
In our case, we can focus just on an enthalpy balance, that is φ = c P T . Diffusivity coefficient Γ can be an anisotropic tensor in general case, so using for example Daly's model (6) for the turbulent heat flux, we can derive the following definition where λ stands here for thermal conductivity. Assuming a steady-state case with no internal heat (enthalpy) sources, constant density and specific enthalpy, substitution to Equation (9) gives which corresponds to Equation (4) for Reynolds averaged temperature T . Formally the same equation as (11) can be obtained for temperature as the user-defined scalar, that is φ = T . The differences should arise in boundary conditions which are summarized in Table 1. The definition of diffusivity coefficient Γ should stay same in both cases as described by Equation (10). Figure 7 presents simulation results of a test case in the laminar flow regime (Re = 270) obtained by a solution Table 1. Two possible definitions of the user-defined scalar φ and consequent differences in boundary conditions. scalar definition boundary flux boundary value of the standard energy equation. It is compared with three variants based on the solution of the user-defined scalar equation. The first one represents the case of isotropic diffusivity coefficient Γ = λ/c P which exactly match the solution of energy equation. The second variant shows an anisotropic situation with 50% larger diffusivity in y-direction, that is in the direction normal to the impinged plate (wall). The third variant represent the test case of 20 times larger diffusivity coefficient in the direction of r-coordinate, that is in the direction from the center to the outer edge of the plate. This figure just clearly illustrates the impact of the anisotropic diffusivity on the Nusselt number and represents a verification that we are able to produce same results with UDS equation as well as with the standard energy equation. A similar test case was done for the turbulent flow regime, Re = 23000, see Figure 8. The flow-field solution is based on the Intermittency Transition turbulence model [4]. The diffusivity in this case is not equal to the molecular thermal conductivity λ only but it is complemented by the term coming from the definition of turbulent heat flux (5) with Pr t = 0.85.
Due to larger velocities along the impinged plate in the turbulent flow regime, the anisotropic diffusivity in r-direction was artificially increased 2000 times in the user-defined function to demonstrate the impact on the heat transfer clearly.

Conclusions
Results of numerical simulations using several turbulence models were compared with experimental data from literature for a round impinging jet. Surprisingly, better agree- The anisotropic diffusivity in r-direction is 2000 times larger, and the flow-field solution is based on Intermittency Transition turbulence model [4]. Experimental data [13] are included for comparison. ment with experimental data was obtained with the turbulence models based on the scalar eddy viscosity (Boussinesq approximation), namely the Transition SST model [11] or Intermittency Transition model [4]. They were able to predict the secondary peak in the Nusselt number dependency on radius with reasonable accuracy in contrast to the second-order closure Reynolds Stress turbulence models. The location of the secondary peak was evaluated as r/D = 2.08 with the Intermittency Transition model which is close to value r/D = 1.97 reported by [20]. With the Reynolds Stress turbulence models, we were not able to get sufficiently correct flow-field predictions compared to [15] which had unfavorable impact on the heat transfer results. LES based simulations have the potential to provide better results and are more universal, of course. But it is questionable if the increase in computational costs by several orders compared to RANS based simulations are compensated by not much better agreement with experimental data [21,20].
The second-order closure Reynolds Stress models are able to overcome the problem with isotropic value of the turbulent (eddy) viscosity, nevertheless the anisotropy of the turbulent heat flux is usually not considered in current CFD codes where only the scalar value of the turbulent viscosity is commonly used in its calculation. Instead of using another complex transport equation for the turbulent heat flux, simpler explicit algebraic heat flux models can be used to describe its dependency on the Reynolds stress components and possibly reflect the anisotropy of the turbulent heat flux itself. The methodology of implementing a user-defined model of the anisotropic turbulent heat flux in ANSYS Fluent through a user-defined scalar (UDS) equation is described in this paper. The UDS simulation results have been validated against the results based on solution of the standard energy transport equation in laminar and turbulent flow regimes and perfect match has been obtained. A test case of the anisotropy impact on the heat transfer (Nusselt number) has been illustrated.