Numerical simulation of jet impingement relaminarization using nonlinear eddy viscosity turbulence models

Abstract

This numerical study investigated a single air jet impingement at a constant temperature. Different RANS turbulence models, namely the v²-f model and four Reynolds stress turbulence models (the linear pressure-strain model, the quadratic pressure-strain model, the stress-omega model and the stress-BSL model) were employed to simulate and compare the numerical results with other published measurement data. Here, we show that the v²-f model has better numerical precision in simulating the local Nusselt number and velocity profiles and is in good agreement with a variety of Reynold numbers. These findings can assist scientists or engineers to improve mathematical turbulence modeling and industrial equipment for further development and to enable innovation.

KEYWORDS:

• Relaminarization
• nonlinear eddy viscosity
• pressure strain term

1. Introduction

Jet impingement is important in many industrial applications, such as cooling internal combustion engines, gas turbines, heat treatment processes, defogging of optical surfaces and semiconductor components. At present, many scientists and engineers are striving to make the detailed calculations involved more efficient and accurate. An issue that is likely to cause an impasse in jet impingement studies is the flow at the stagnation point and a combination between laminar-turbulent transitions on the hot surface. The transition from a laminar to a turbulent regime is frequently triggered by instability (Ichimiya, 1999). In jet impingement, secondary vortices may occur due to the turbulent flow field along the wall and can be caused by turbulent variations in lateral/radial velocity related to pressure gradient changes, as seen in Figure 1, which can lead to separation or relainarization zones. Creating a separation or relainarization zone may increase localized heat transfer rates (Cousteix, 2003; Dewan et al., 2012).

Figure 1. The vortex motion in the impinging jet.

The main reason for the poor rate of heat transfer happens the jet flow characteristic approaching the surface-induced stagnation or laminar deceleration region. The flow builds up a higher static pressure on the surface and loses the flow momentum in the radial direction, substantially impacting local turbulence and heat transfer properties (Zuckerman & Lior, 2006).

To better understand the previous scientific findings, this paper is organized as follows: section 1 describes the experimental research, numerical simulation study and exciting applications; section 2 explains the mathematical and numerical modeling; section 3 presents the computational investigation methodology; section 4 presents the computational investigation results; and section 5 provides the discussion of the results. Finally, section 6 discusses the main conclusions and findings.

1.1. Experiment research

Here, we have reviewed the experimental research related to jet impingement, which includes a hypothesis, a variable that the researcher can manipulate, and variables that can be measured. Initially, Kalifa et al. (2016) investigated the dynamics of the turbulent behavior of jets impinging on a flat plate. The experimental investigations of the characteristics on a flat surface can be re depicted and conducted using particle image velocimetry. An impinging jet's behavior clearly shows that the co-flow stream imposed noticeable restrictions on the spreading of the impinging jet and considerably reduced the air jet's entrainment. Bader et al. (2017) experimentally investigated relaminarization in accelerated flat-plate boundary-layer flows using laser Doppler anemometry. To understand the processes in the relaminarization mechanism, the mean value of velocity fluctuation was investigated and showed that strong acceleration leads to a suppression of turbulence production. Chen et al. (2018) evaluated the local flow and average heat transfer characteristics of jet impingement cooling with supercritical pressure fluids using carbon dioxide. The microelectromechanical system (MEMS) method and a low thermal conductivity substrate were integrated with a thermal sensor chip to read and collect the hot surface temperature. The test data showed that the local heat transfer coefficient at the stagnation point improved with increasing heat flux. Singh et al. (2020) determined the effect of the Reynolds number, nozzle-to-plate distance, and loading fraction. The results identified that raising the Reynolds number and loading fraction increased the Nusselt number but increasing the H/D ratio decreased the cooling performance. In addition, they proposed a correlation function based on the experimental data that a proper correlation function might be the best option in jet impingement studies. Issac et al. (2020) applied the impingement of a round jet on a heated flat plate with constant heat flux, which was investigated experimentally and numerically at four different nozzle-to-plate spacings. The numerical results showed that the RNG k-ε turbulence model accurately predicted the outcomes except at stagnation, the k-$\omega$ shear stress transport (SST) turbulence model was in good agreement at various H/D ratios only for low Reynolds numbers, and the v²-f turbulence model predicted a more accurate Nusselt number at the stagnation results. Liu et al. (2020) utilized the volume of fluid with pulse inlet pressure combined with the large eddy simulation (LES) turbulence model to understand the erosion area. They concluded that the frequency of the vortices shedding from the jet nozzle dominated pressure pulsation in the first ring area. Leocadio and Passos (2021) experimentally and computationally investigated heat transport during water jet impingement cooling of a heated steel plate. The heat transfer coefficient (HTC) increased as the surface temperature decreased, and the correlations were used to predict the HTC and the rewetting front radius. This finding with a new correlation should enhance cooling system efficiency in hot strip mills.

1.2. Numerical simulation study

Many real-world engineering problems often encounter a design device with complicated geometry. Jet impingement has limitations in solving complex engineering problems, and the best-in-class turbulence model is still an open question. Large eddy simulation approach shows as a promising technique to simulate turbulent flows which work with many engineering applications (Zahiri & Roohi, 2021; Zahiri & Roohi, 2019; Lar-Kermani et al., 2018). For LES jet impingements, Shademan et al. (2016) used an LES to analyze a turbulent impingement jet with a Reynold number of 28,000 compared with the particle image velocimeter technique. The results suggested that the mean normalized centerline velocity in the free jet, impingement areas, and pressure distribution were consistent between the models and tests. Also, they found that the vortices caused form the Kelvin-Helmholtz instabilities near the nozzle exit were combined, broken down and changed into large-scale structure configuration. Aillaud et al. (2016) focused on a wall-resolved large eddy simulation (LES) of an adiabatic, circular, and submerged impinging jet on a hot flat plate, using a Reynolds number of 23,000 and a nozzle-to-plate distance (H/D) of 2.0. Two conclusions based on the temporal mean Nusselt number were described. First, the LES data provided an excellent level of accuracy compared to existing experimental data. Second, via high-order statistics used to discover the fluid characteristics near the second peak of the Nusselt number, the time variation of pressure distribution revealed the rebound of the main vortices close to the location of the secondary peak and permitted identification of the beginning of unstable separation. Grenson and Deniau (2017) concluded that LES provided the most accurate predictions of the secondary maximum in the Nusselt number of the impinging jet when the hot inlet jet exits have fully developed profiles.

Similarly, Bader et al. (2017) concluded that the viscous effects were significant for a narrow boundary layer flow over a solid body surface. The wall friction and heat transfer characteristics can affect laminar or turbulent flow conditions, causing a reverse transition or relaminarization from turbulent to laminar. Next, Dey et al. (2018) investigated the flow with the pre-attachment region and found that the velocity profiles follow an exponential function. Huang et al. (2018) studied numerical heat transfer with slot jet impingement by employing the combination turbulence model. The Kato-Launder model, the intermittency transition model, and the crossflow transition model were merged with the primarily original SST turbulence model. When compared to both experimental data and accessible numerical results from previous works, it was discovered that the proposed model performs well regarding various parameters.

In 2020, Pakhomov et al. numerically studied the effect of the pulse form on heat transfer in a pulsed impinging jet using the unsteady RANS model and second-moment closure. The study results found that the heat transfer at the stagnation point varies according to the pulse frequency and the Reynolds number. Singh et al. (2020) investigated the fluid flow properties of a turbulent offset jet impinging on a wavy wall surface. The results provide that pressure near the reattachment point decreases and the magnitude of the decay in the normalized streamwise maximum velocity depends on the wavy surface's amplitude. Similarly, Zhang et al. (2020) investigated ultrafast cooling equipment using the self-excited pulsed nozzle. The results reveal that the self-excited pulsed jet can enhance the perturbations of the boundary layer and significantly increase the turbulence intensity of the plate surface.

1.3. Current and future interesting applications

According to the jet impingement phenomenon, heat transfer intensification methods are vital in many industrial applications (Chen & Chen, 2022; Steinberg et al., 2022; Ul Haq et al., 2022; Xing et al., 2022; Zhou et al., 2022). Regarding the importance of a wide range of applications, such as aerospace, automotive and semiconductor components, the following sections discuss the applications that require a fidelity standard of accuracy. Gao et al. (2018) investigated the safety impact of gas jets on a jet blast deflector of aircraft engines that significantly interferes with the flight deck environment. This work confirmed that those predicted results from SST k-ω turbulence model was better for compressible viscous gas jet simulations. Zhou et al. (2018) proposed a new jet impingement chip-scale cooling structure to improve silicon wafer manufacturing. This unique design was estimated to provide 70% higher cooling performance for pumping power and more uniform cooling. Radmard et al. (2020) explored the heat transfer from directly manufactured micro pins undergoing jet impingement. Their study showed that the outlet port position significantly affected the hydraulic and thermal resistances. Rane and He (2020) analyzed a two-phase turbine used for power generation in thermal power plants based on a Ranz-Marshall heat transfer correlation and concluded there was good agreement in the BNL nozzle.

However, as late as 2021, the aero-acoustic properties of dual impinging jets remained less understood. Mehta et al. (2021) adopted the shadowgraph technique to investigate the aero-acoustic properties of dual impinging jets. Their findings indicated that changing the nozzle spacing affected the fountain flow strength. The aero-acoustic broadband and discrete amplitude tones were also affected by changing the spacing, with the overall sound pressure level being robust to the gap between nozzles. Erasmus et al. (2021) studied the application of impinging jets in thermal receivers and bayonet tube heat exchangers for concentrating solar power using the transition SST RANS turbulence model and proposed an empirical correlation of the pressure loss factor. Their numerical investigation clearly showed that the dominating total pressure loss occurred when there was rapid jet expansion. Kura et al. (2021) applied a jet impingement on a non-planar surface for a wide range of geometries designed for heat exchangers. The numerical simulations were carried out using an advanced ζ-f turbulence model implemented in OpenFOAM. The results suggested that the geometrical shape influenced the secondary peak in the Nusselt number distribution. Cademartori et al. (2021) employed CFD tools to demonstrate a novel cooling technique through slot jet impingement and compared those of predicted results with an available vertical galvanizing industrial line primarily intended using a round jet shape. From their investigation, slot jet impingement was made more efficient and favorable by increasing the jet wall distance and thus improving the uniformity of transversal temperature distributions. Most recently, Das et al. (2022) applied jet impingement on the absorber plate to enhance the effectiveness of a solar air heater. Their results showed that superior improvement of elongated streamlines and higher turbulent kinetic energy could improve the heat transfer performance.

The above literature review demonstrates that the reverse transition from turbulent to laminar flow has been studied experimentally by many investigators since it not only has a very complex nature but also offers fundamental knowledge into how to control the turbulence mechanism, as well as expand this capability from academic research to practical industrial work. Additionally, it is not easy to understand the flow mechanism in any field with a general turbulence model. Therefore, to increase the understanding of the jet impingement behavior due to the variation of the anisotropic turbulence mechanism, the current study presents all the hidden turbulence equations in the calculation framework and focuses on a new turbulence model that is implemented in commercial software by including separation-induced transition correction such as v²-f, and Reynolds Stress based on $\epsilon$-approach and the $\omega$-approach. Furthermore, these numerical investigation results are compared with previous experimental results (Baughn et al., 1991; Cooper et al., 1993).

2. Mathematical and numerical modeling

This section provides turbulence model descriptions considered in the current study. Initially the v²-f model is proposed, which includes the effect of an anisotropy turbulent structure and changes to appropriate velocity scale near walls. Then, the Reynold stress model (RSM) model is presented, which considers the effects of streamline curvature and rapid changes in strain rate tensors. In addition, the RSM in this study is considered using $\epsilon$-based and $\omega$-based approaches, because both models are still questionable regarding their capability in an impingement jet numerical study. All the equation systems are shown in the next section.

2.1. The v²-f turbulence model

The v²-f model is a RANS-based turbulence model which does not require a wall function. It includes additional terms which take account of near-wall anisotropy turbulence and nonlocal pressure-strain in the model. For the standard k – $\epsilon$ model, the primitive velocity and time scales are dependent on the turbulent kinetic energy, $k^{1/2}$ and $k/\epsilon$ respectively. Since those of previous scales failed to predict wall-bounded flow applications, the v²-f model applies the velocity scale ${\stackrel{\leftharpoonup }{v}}^2$ instead of $k$ for an improved scale near walls, which can match with the boundary condition for the normal component. As a result, the v²-f model has an advantage over the previous models because it can provide proper scaling in representing the damping of turbulent transport close to the wall.

The v²-f model contains four transport equations. Clearly, both the turbulence kinetic energy and the turbulence dissipation equations are employed with those of the k-$\epsilon$ model but then they use different eddy viscosity assumptions and another turbulence time scale. Also, the ${\stackrel{\leftharpoonup }{v}}^2$ transport equation and the elliptic-relaxation equation are included. The effects of the anisotropic are included and the Helmholtz elliptic equation are solved to define the elliptic relaxation function. For the computational performance improvement of the v²-f model, a variant of this eddy-viscosity model is derived when the change in variables is introduced. This turbulence variable can be regarded as the ratio of the two-time scales: scalar $k/\epsilon$ (isotropic) and lateral $\overline{{\upsilon }^2}/\epsilon$ (anisotropic). More information on the v²-f model can be found in ANSYS Inc. (2021)

2.1.1 Transport equation for the v²-f model

The turbulence kinetic energy, turbulence dissipation rate, velocity variance scale, and elliptic relaxation function are derived as follows: (1) $\begin{array}{rl}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\rho k{u}_i)& =P-\rho \epsilon +\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[(\mu +{\mu }_t)\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\right]\end{array}$(1) (2) $\begin{array}{rl}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\rho \epsilon u_i)& =\frac{C_{\epsilon 1}^{\mathrm{\prime }}P-1.9\rho \epsilon }{T}\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[(\mu +0.7692{\mu }_t)\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{x}_j}\right]\end{array}$(2) (3) $\begin{array}{rl}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\rho \overline{v^2}{u}_i)& =\rho kf-6\rho \overline{v^2}\frac{\epsilon }{k}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[(\mu +{\mu }_t)\frac{\mathrm{\partial }\overline{v^2}}{\mathrm{\partial }{x}_j}\right]\end{array}$(3) (4) $\begin{array}{rl}f-L^2\frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}_j^2}& =\frac{0.2667-0.4\overline{\frac{v^2}{k}}}{T}+\frac{3P}{10\rho k}+\frac{5\overline{v^2}/k}{T}\end{array}$(4) where $P=2{\mu }_t{S}^2$, $S^2\equiv S_{ij}{S}_{ij}$, and $S_{ij}=\frac{1}{2}\left(\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}\right)$. The turbulent length (L) and time (T) scales are expressed as follow: (5) $\begin{array}{rl}{T}^{\mathrm{\prime }}& =\text{max}\left[\frac{k}{\epsilon },6\sqrt{\frac{\mu }{\rho \epsilon }}\right]\end{array}$(5) (6) $\begin{array}{rl}T& =\text{min}\left[T^{\mathrm{\prime }},\frac{1.5755k}{\overline{v^2}\sqrt{2S^2}}\right]\end{array}$(6) (7) $\begin{array}{rl}{L}^{\mathrm{\prime }}& =\text{min}\left[\frac{k^{1.5}}{\epsilon },\frac{2.6244k^{1.5}}{\overline{v^2}\sqrt{2S^2}}\right]\end{array}$(7) (8) $\begin{array}{rl}L& =0.23\text{max}\left[L^{\mathrm{\prime }},70{\left(\frac{{\mu }^3}{{\rho }^3\epsilon }\right)}^{0.25}\right]\end{array}$(8) (9) $\begin{array}{rl}C{}_{\epsilon 1}^{\mathrm{\prime }}& =1.4+0.063\sqrt{k/\overline{v^2}}\end{array}$(9) (10) $\begin{array}{rl}{\mu }_t& =0.22\rho \overline{v^2}T\end{array}$(10)

The variable $f$ is solved from the elliptic relaxation equation (4(4) $\begin{array}{rl}f-L^2\frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}_j^2}& =\frac{0.2667-0.4\overline{\frac{v^2}{k}}}{T}+\frac{3P}{10\rho k}+\frac{5\overline{v^2}/k}{T}\end{array}$(4) ). Here, the model utilizes an elliptic operator to calculate a term comparable to the pressure-strain correlation of the RSM. Ellipticity is described by a modified Helmholtz operator, introducing wall effects via a linear differential equation (ANSYS Inc., 2021).

2.2. The reynolds stress model

The Reynolds stress model (RSM) is the most decorative RANS turbulence model. In particular, the model of RSM has higher-level turbulence closures and represents a complete classical turbulence model. The method of closure applied is usually called second-order closure. In the RSM, the stress tensor components of the Reynolds stress are directly determined while the eddy viscosity assumption is avoided. This model relies on the exact Reynolds stress transport equation by couple solving between the transport equations of Reynolds stresses and the dissipation rate; an additional set of seven transport equations are required for complex three-dimensional problems. The advantage of the RSM is that it includes the effects from swirling motion, rotational motion, rapid change in strain rate, and streamline curvature, so that it could be expected to predict jet impingement behavior more accurately. In particular, the $\epsilon$-based (linear and quadratic pressure-strain model) and $\omega$-based (stress-omega and stress-BSL model) RSMs are used to model the pressure strain and dissipation rate terms, respectively. They are accountable for compromising the accurate prediction of RSM results.

2.2.1 Reynolds stress transport equations

The complete exact transport equations for the transport of the Reynolds stress, $\rho \overline{u{}_j^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}}$, are written as: (11) $\begin{array}{rl}\underset{C_{ij}\equiv Convection}{\underbrace{\frac{\mathrm{\partial }}{{\mathrm{\partial }}_k}(\rho u_k\overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}})}}& =\underset{D_{T,ij}\equiv TurbulentDiffusion}{\underbrace{-\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_k}[\rho \overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}+\overline{p^{\mathrm{\prime }}({\delta }_{kj}u{}_i^{\mathrm{\prime }}+{\delta }_{ik}u{}_j^{\mathrm{\prime }})}]}}\\ & \underset{M_{L,ij}\equiv MolecularDiffusion}{\underbrace{+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_k}\left[\mu \frac{\mathrm{\partial }}{{\mathrm{\partial }}_k}(\overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}})\right]}}\\ & +-\underset{P_{ij}\equiv StressProduction}{\underbrace{\rho \left[\overline{u{}_i^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_k}+\overline{u{}_j^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_k}\right]}}\\ & +\underset{{\psi }_{ij}\equiv PressureStrain}{\underbrace{p^{\mathrm{\prime }}\left(\frac{\mathrm{\partial }u{}_i^{\mathrm{\prime }}}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }u{}_j^{\mathrm{\prime }}}{\mathrm{\partial }{x}_i}\right)}}-\underset{{\epsilon }_{ij}\equiv Dissipation}{\underbrace{\overline{2\mu \frac{\mathrm{\partial }u{}_i^{\mathrm{\prime }}}{\mathrm{\partial }{x}_k}\frac{\mathrm{\partial }u{}_j^{\mathrm{\prime }}}{\mathrm{\partial }{x}_k}}}}\end{array}$(11) In equation (11(11) $\begin{array}{rl}\underset{C_{ij}\equiv Convection}{\underbrace{\frac{\mathrm{\partial }}{{\mathrm{\partial }}_k}(\rho u_k\overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}})}}& =\underset{D_{T,ij}\equiv TurbulentDiffusion}{\underbrace{-\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_k}[\rho \overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}+\overline{p^{\mathrm{\prime }}({\delta }_{kj}u{}_i^{\mathrm{\prime }}+{\delta }_{ik}u{}_j^{\mathrm{\prime }})}]}}\\ & \underset{M_{L,ij}\equiv MolecularDiffusion}{\underbrace{+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_k}\left[\mu \frac{\mathrm{\partial }}{{\mathrm{\partial }}_k}(\overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}})\right]}}\\ & +-\underset{P_{ij}\equiv StressProduction}{\underbrace{\rho \left[\overline{u{}_i^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_k}+\overline{u{}_j^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_k}\right]}}\\ & +\underset{{\psi }_{ij}\equiv PressureStrain}{\underbrace{p^{\mathrm{\prime }}\left(\frac{\mathrm{\partial }u{}_i^{\mathrm{\prime }}}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }u{}_j^{\mathrm{\prime }}}{\mathrm{\partial }{x}_i}\right)}}-\underset{{\epsilon }_{ij}\equiv Dissipation}{\underbrace{\overline{2\mu \frac{\mathrm{\partial }u{}_i^{\mathrm{\prime }}}{\mathrm{\partial }{x}_k}\frac{\mathrm{\partial }u{}_j^{\mathrm{\prime }}}{\mathrm{\partial }{x}_k}}}}\end{array}$(11) ), $C_{ij}$, $M_{L,ij}$, and $P_{ij}$ do not require additional modeling. However, $D_{T,ij}$, ${\psi }_{ij}$, and ${\epsilon }_{ij}$ need additional modeled equations to describe the assumption to close the equation set.

2.2.2 Modeling the dissipation rate

The scalar dissipation rate computed from the steady transport equation is: (12) $\begin{array}{rl}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\rho \epsilon u_i)& =0.72\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\mu +\frac{9\rho k^2}{100\epsilon }\right)\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{x}_j}\right]P_{ii}\frac{\epsilon }{k}\\ & -1.92\rho \frac{{\epsilon }^2}{k}\end{array}$(12) For the RSM based on the $\epsilon$-equation, the dissipation tensor ${\epsilon }_{ij}$ is written as: (13) ${\epsilon }_{ij}=\frac{2}{3}{\delta }_{ij}(\rho \epsilon )$(13)

However, where the RSM is coupled based on the $\omega$ or $BSL$ equation, the dissipation tensor is modeled as: (14) ${\epsilon }_{ij}=\frac{2}{3}{\delta }_{ij}\rho {\beta }_{RSM}^{\ast }k\omega$(14) where ${\beta }_{RSM}^{\ast }$ is computed from modeling the pressure-strain term as described in the next section. For the stress-omega model and stress-$BSL$ model, the specific dissipation rate is calculated similarly as the standard k-$\omega$ model and the baseline k-$\omega$ model, respectively (ANSYS Inc., 2021).

2.2.3 Modeling the turbulence kinetic energy

The trace of the Reynolds stress tensor is introduced to model the turbulence kinetic energy (TKE) term as follow: (15) $k=\frac{1}{2}\overline{u{}_i^{\mathrm{\prime }}u{}_i^{\mathrm{\prime }}}$(15) For the $\epsilon$-equation based RSMs, a transport equation for the TKE is solved to obtain boundary conditions for the Reynolds stresses using the following model equation: (16) $\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\rho k{u}_i)=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\mu +0.09\rho \frac{k^2}{\epsilon }\right)\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\right]+0.5P_{ii}$(16)

2.2.4 Modeling the pressure-strain term

In this study, we use two $\epsilon$-based RSMs (linear and quadratic pressure-strain models) and $\omega$-based RSMs (stress-omega and stress-$BSL$ models), respectively. The elaborate equation system details are as follows:

2.2.4.1. Linear Pressure-Strain Model (RSM-LPS)

The classical approach to modeling the pressure-strain term (${\psi }_{ij}$) uses the following decomposition according to proposals by Gibson and Launder (1978) and Launder (1989): (17) ${\psi }_{ij}={\psi }_{ij,1}+{\psi }_{ij,2}+{\psi }_{ij,w}$(17) The slow pressure-strain term, ${\psi }_{ij,1}$, is modeled by: (18) ${\psi }_{ij,1}=-1.8\rho \frac{\epsilon }{k}\left[\overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}}-\frac{2}{3}{\delta }_{ij}k\right]$(18)

The rapid pressure-strain term, ${\psi }_{ij,2}$, is modeled by: (19) ${\psi }_{ij,2}=-0.60[(P_{ij}-C_{ij})]-\frac{1}{3}{\delta }_{ij}(P_{kk}-C_{kk})$(19) The wall-reflection term, ${\psi }_{ij,w}$, represents the redistribution of near wall-normal stresses. It tends to dampen the normal stress perpendicular to the wall while enhancing the stresses parallel to the wall. The calculation model is: (20) $\begin{array}{rl}{\psi }_{ij,w}& \equiv \frac{\epsilon }{2k}(\overline{u{}_k^{\mathrm{\prime }}u{}_m^{\mathrm{\prime }}}{n}_k{n}_m{\delta }_{ij}-\frac{3}{2}\overline{u{}_i^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}{n}_j{n}_k\\ & -\frac{3}{2}\overline{u{}_j^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}{n}_i{n}_k)\varphi \\ & +0.3({\psi }_{km,2}{n}_k{n}_m{\delta }_{ij}-\frac{3}{2}{\psi }_{ik,2}{n}_j{n}_k\\ & -\frac{3}{2}{\psi }_{jk,2}{n}_i{n}_k)\varphi \end{array}$(20)

In equation (20(20) $\begin{array}{rl}{\psi }_{ij,w}& \equiv \frac{\epsilon }{2k}(\overline{u{}_k^{\mathrm{\prime }}u{}_m^{\mathrm{\prime }}}{n}_k{n}_m{\delta }_{ij}-\frac{3}{2}\overline{u{}_i^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}{n}_j{n}_k\\ & -\frac{3}{2}\overline{u{}_j^{\mathrm{\prime }}u{}_k^{\mathrm{\prime }}}{n}_i{n}_k)\varphi \\ & +0.3({\psi }_{km,2}{n}_k{n}_m{\delta }_{ij}-\frac{3}{2}{\psi }_{ik,2}{n}_j{n}_k\\ & -\frac{3}{2}{\psi }_{jk,2}{n}_i{n}_k)\varphi \end{array}$(20) ), $n_k$ is the $x_k$unit normal vector in each component, $d$ is the normal distance to the wall, and $\varphi =0.3924k^{3/2}/\epsilon d$.

2.2.4.2. Quadratic Pressure-Strain Model (RSM-QPS)

This model has been revealed to provide advanced performance for a range of basic shear flows, including plane strain, rotating plane shear, and asymmetrical expansion/contraction. It leads to produce a high level of accuracy for complex engineering flows, especially in streamline curvature, as proposed by Speziale et al. (1991) and written as: (21) $\begin{array}{rl}{\psi }_{ij}& =-(3.4\rho \epsilon +1.8P){\tau }_{ij}\\ & +4.2\rho \epsilon \left({\tau }_{ik}{\tau }_{kj}-\frac{1}{3}{\tau }_{mn}{\tau }_{mn}{\delta }_{ij}\right)\\ & +\left(0.8-1.3\sqrt{{\tau }_{ij}{\tau }_{ij}}\right)\rho k{S}_{ij}\\ & +1.25\rho k\left({\tau }_{ik}{S}_{jk}+{\tau }_{jk}{S}_{ik}-\frac{2}{3}{\tau }_{mn}{S}_{mn}{\delta }_{ij}\right)\\ & +0.4\rho k({\tau }_{ik}{\mathrm{\Omega }}_{jk}+{\tau }_{jk}{\mathrm{\Omega }}_{ik})\end{array}$(21) where ${\tau }_{ij}$ represent the Reynolds stress anisotropy tensor and ${\mathrm{\Omega }}_{ij}$ is the mean rate of rotation tensor defined by: (22) $\begin{array}{rl}{\tau }_{ij}& =-\left(\frac{-\rho \overline{u{}_i^{\mathrm{\prime }}u{}_j^{\mathrm{\prime }}}+\frac{2}{3}\rho k{\delta }_{ij}}{2\rho k}\right)\end{array}$(22) (23) $\begin{array}{rl}{\mathrm{\Omega }}_{ij}& =\frac{1}{2}\left(\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}-\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}\right)\end{array}$(23)

2.2.4.3. Stress-Omega Model (RSM-SO)

The stress-omega model is a stress transport model based on the ω-equations and the LRR model. It resembles the κ−ω model due to its excellent predictions for a wide range of turbulent flows requiring wall reflection treatments. The relation of pressure-strain can be written by: (24) $\begin{array}{rl}{\psi }_{ij}& ={\psi }_{ij,1}+{\psi }_{ij,2}\end{array}$(24) Therefore, (25) $\begin{array}{rl}{\psi }_{ij}& =-C_1\rho {\beta }_{RSM}^{\ast }\omega \left[\overline{u_i^{\mathrm{\prime }}{u}_j^{\mathrm{\prime }}}-\frac{2}{3}{\delta }_{ij}k\right]\\ & -\overline{\alpha }\left[P_{ij}-\frac{1}{3}{P}_{kk}{\delta }_{ij}\right]\\ & -\overline{\beta }\left[{\mathrm{\Gamma }}_{ij}-\frac{1}{3}{P}_{kk}{\delta }_{ij}\right]\\ & -\overline{\gamma }k\left[S_{ij}-\frac{1}{3}{S}_{kk}{\delta }_{ij}\right]\end{array}$(25) where $C_1$, ${\beta }_{RSM}^{\ast }$, $\overline{\alpha }$, $\overline{\beta }$, ${\mathrm{\Gamma }}_{ij}$,$f_{\beta }^{\ast }$, and $\overline{\gamma }$ are defined by: (26) $\begin{array}{rl}{C}_1& =1.8\left[\frac{(5/3)+(R{e}_t/12)}{1+(R{e}_t/12)}\right]\end{array}$(26) (27) $\begin{array}{rl}{\beta }_{RSM}^{\ast }& =0.09\left[\frac{0.2667+{(R{e}_t/8)}^4}{1+{(R{e}_t/8)}^4}\right]f_{\beta }^{\ast }\end{array}$(27) (28) $\begin{array}{rl}{f}_{\beta }^{\ast }& =\{\begin{array}{ll}1& {\chi }_k\le 0\\ \frac{1+640{\chi }_k^2}{1+400{\chi }_k^2}& {\chi }_k>0\end{array}\end{array}$(28) (29) $\begin{array}{rl}\overline{\alpha }& =\frac{1+0.0645R{e}_t}{1+(R{e}_t/12)}\end{array}$(29) (30) $\begin{array}{rl}\overline{\beta }& =0.1964\left[\frac{R{e}_t/12}{1+(R{e}_t/12)}\right]\end{array}$(30) (31) $\begin{array}{rl}{\mathrm{\Gamma }}_{ij}& =-\rho \left[\overline{u{}_i^{\mathrm{\prime }}u{}_m^{\mathrm{\prime }}}\frac{\mathrm{\partial }{u}_m}{\mathrm{\partial }{x}_j}+\overline{u{}_j^{\mathrm{\prime }}u{}_m^{\mathrm{\prime }}}\frac{\mathrm{\partial }{u}_m}{\mathrm{\partial }{x}_i}\right]\end{array}$(31) (32) $\begin{array}{rl}\overline{\gamma }& =0.4945\left[\frac{0.007+(R{e}_t/12)}{1+(R{e}_t/12)}\right]\end{array}$(32)

Where (33) $\begin{array}{rl}{\chi }_k& \equiv \frac{1}{{\omega }^3}\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\frac{\mathrm{\partial }\omega }{\mathrm{\partial }{x}_j}\end{array}$(33) (34) $\begin{array}{rl}R{e}_t& =\frac{\rho k}{\mu \omega }\end{array}$(34)

2.2.4.4. Stress-BSL Model (RSM-S-BSL)

More specifically, the stress-BSL model is the $\omega$-based Reynold stress model, as presented by Speziale et al. (1991), which solves the scale equation from the baseline (BSL) k-$\omega$ model and adopts a blending function for the transport of turbulent dissipation in the same way as the SST k-$\omega$ turbulence model. It uses the same relation in equation (25(25) $\begin{array}{rl}{\psi }_{ij}& =-C_1\rho {\beta }_{RSM}^{\ast }\omega \left[\overline{u_i^{\mathrm{\prime }}{u}_j^{\mathrm{\prime }}}-\frac{2}{3}{\delta }_{ij}k\right]\\ & -\overline{\alpha }\left[P_{ij}-\frac{1}{3}{P}_{kk}{\delta }_{ij}\right]\\ & -\overline{\beta }\left[{\mathrm{\Gamma }}_{ij}-\frac{1}{3}{P}_{kk}{\delta }_{ij}\right]\\ & -\overline{\gamma }k\left[S_{ij}-\frac{1}{3}{S}_{kk}{\delta }_{ij}\right]\end{array}$(25) ) for the pressure-strain correlation in the Reynolds stress equations as the stress-omega model. The coefficient ${\beta }_{RSM}^{\ast }$ is equal to ${\beta }_{RSM}^{\ast }/f_{\beta }^{\ast }$.

2.2.5 Modeling turbulent diffusive transport

For the $\epsilon$-based RSM, $D_{T,ij}$ can be modeled by a scalar turbulent diffusivity equation: (35) $D_{T,ij}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_k}\left(\frac{{\mu }_t}{{\sigma }_k}\frac{\mathrm{\partial }\overline{u_i^{,}{u}_j^{\mathrm{\prime }}}}{\mathrm{\partial }{x}_k}\right)$(35) The turbulent viscosity,${\mu }_t$, for the $\epsilon$-based model is computed using: (36) ${\mu }_t=0.09\frac{\rho k^2}{\epsilon }$(36)

The $\omega$-equation (stress-omega) and BSL-equation (stress-BSL) calculated from: (37) ${\mu }_t=\frac{0.144\mu \omega +\rho k}{6\mu \omega +\rho k}$(37) In equation (35(35) $D_{T,ij}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_k}\left(\frac{{\mu }_t}{{\sigma }_k}\frac{\mathrm{\partial }\overline{u_i^{,}{u}_j^{\mathrm{\prime }}}}{\mathrm{\partial }{x}_k}\right)$(35) ), the model constant is ${\sigma }_k=0.82$. However, for the $\omega$-based approach, the stress-omega model sets ${\sigma }_{k,1}=2.0$ and ${\sigma }_{k,2}=1.0$, whereas the stress-BSL models use the blending function $F_1$ from: (38) $\begin{array}{rl}{\sigma }_k& =F_1{\sigma }_{k,1}+(1-F_1){\sigma }_{k,2}\end{array}$(38) (39) $\begin{array}{rl}{F}_1& =tanh({\phi }_1^4)\end{array}$(39) (40) $\begin{array}{rl}{\phi }_1& =min\left[max\left(\frac{\sqrt{k}}{0.09\omega y},\frac{500\mu }{\rho y^2\omega }\right),\frac{4\rho k}{1.168D_{\omega }^{+}{y}^2}\right]\end{array}$(40) (41) $\begin{array}{rl}{D}_{\omega }^{+}& =max\left[\frac{2\rho }{1.168}\frac{1}{\omega }\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\frac{\mathrm{\partial }\omega }{\mathrm{\partial }{x}_j},{10}^{-10}\right]\end{array}$(41) where $D_{\omega }^{+}$ represents the positive portion of the cross-diffusion term and $y$ is the distance to the next surface.

3. Computational investigation

The summary of the computational domain is shown in Figure 2. The configuration includes the circular orifice and the target constant heat flux plate ($\dot{q}$ = 300 $W/m^2$). The jet Reynolds number of 23,000 based on the exit nozzle diameter is ejected out from a nozzle. According to previous experimental literature reviews (Baughn et al., 1991; Fenot et al., 2005; Hee Lee et al., 2002; Attalla, 2015), the golden ratio of H/D = 2 is capable of clearly generating the second peak of the Nusselt number, which is the focus parameter for the investigation of the turbulence model in the current paper, where ANSYS Fluent 2020 R1 is used. The inlet velocity profiles are fully developed flow, with no-slip and adiabatic wall boundary conditions applied at the pipe surfaces (Grenson & Deniau, 2017).

Figure 2. Computational domain for air impinging jet on a flat plate Re = 23,000 and H/D = 2.

3.1. Grid independent test

The computational grid in Figure 3 was generated using the ANSYS meshing software with the multizone method. The rest of the computational grid uses a pure structural grid. The hot flat plate surface was created with $1$ micrometer first layer cells and 140 circumferences of grid points to obtain $y^{+}<1$ (for boundary layer relaminarization or separation recommended by ANSYS 2020R1, 2020). The grid-independent analysis was carried out over five successively finer grids, as shown in Table 1. All grids shared the same stretching ratio of 1.10 for the structural grid while loosely refined at far outlet boundaries where to save computational cost, the gradients of the flow variables were expected to be negligible.

Figure 3. Computational grid.

Table 1. Computational grid.

Case	circumferential division	radial division	nodes	elements
Grid-1	70	60	479,656	457,800
Grid-2	90	80	807,512	773,505
Grid-3	100	90	2,578,204	2,537,190
Grid-4	120	100	5,337,469	5,272,320
Grid-5	140	140	8,159,200	8,245,961

The comparison was performed using the Nusselt number on the grid size, with various Nusselt profiles plotted against the $r/D$ ratio, as shown in Figure 4, compared to previous publication (Baughn et al., 1991). According to the findings, the Nusselt profile for Grid-5 was consistent with the experiment. Another important finding in that study was that the grid distribution in the radial direction must be dented sufficiently to obtain the second peak of the Nusselt number prediction. As a result, the current study will investigate all cases using the same grid number.

Figure 4. Comparison of Nusselt numbers for grid-independent analysis.

3.2 Numerical setting and initial conditions

The initialization of the incoming jet into the domain required velocity profiles, the TKE, and the specific dissipation rate at the inflow boundary, solved from a long pipe until fully developed. This condition is matched with those of available experiment. Furthermore, the boundary conditions and numerical settings are described in Table 2.

• Inlet: inlet boundary with fully developed flow characteristics, such as velocity, TKE, and the specific dissipation rate

• Top auxiliary inlet: pressure inlet with the inflow direction

• Outlet: pressure: outlet boundary with a constant reference pressure

• Pipe wall: nonslip adiabatic wall

• Heat bottom wall: no-slip wall specified constant heat flux, $\dot{q}$= 300 W/m²

Table 2. Numerical setup.

Numerical parameter	settings	turbulence model
Solver	pressure-based	v^2-f, all RSMs
Pressure-Velocity Coupling	Coupled	v^2-f, all RSMs
Multigrid	F-Cycle	v^2-f, all RSMs
Under-relaxation factors
	0.45 – momentum	v^2-f, all RSMs
	0.45 – pressure	v^2-f, all RSMs
	0.65 – tke	v^2-f, all RSMs
	0.65 – specific dissipation rate	v^2-f, all RSMs
Special discretization	$2^{\text{nd}}$ order – pressure	v^2-f, all RSMs
	$2^{\text{nd}}$ order upwind – momentum	v^2-f, all RSMs
	$2^{\text{nd}}$ order upwind – tke	v^2-f, all RSMs
	$2^{\text{nd}}$ order upwind – momentum	v^2-f, all RSMs
	$2^{\text{nd}}$ order upwind – energy	v^2-f, all RSMs
	$2^{\text{nd}}$ order upwind – velocity variance scale	$v^2-f$
	$2^{\text{nd}}$ order upwind – elliptic relaxation function	$v^2-f$
Gradient	Green-gauss node-based

4. Results

This section is divided into three sections. The first part shows the predicted results of the average Nusselt number on a flat plate with a constant heat flux thermal boundary condition. The available experimental data from Baughn et al. (1991) are also plotted for comparison. The second part demonstrates the turbulence quantities versus the axial vertical distance. Finally, the contour plots of the Nusselt number and velocity profiles near the transition wall zones are shown.

4.1. Secondary peak occurrence due to variation in anisotropic turbulence model

To capture the secondary peak phenomenon, a simulation with a Reynolds number of 23,000 is performed and tested with common turbulence models such as v²-f, RSM-LPS, RSM-QPS, RSM-SO, and RSM-S-BSL. These turbulence model variations are applied with the constant heat flux at the bottom plate. The numerical investigation, shown in Figure 5, describes the distribution of the local Nusselt number magnitude against the $r/D$ ratio. The secondary peak of the Nusselt number magnitude grows in the distribution at $r/D$ = 2.0 and a value of 110. It can be deduced that all anisotropic models, except for the v²-f turbulence model, fail to predict this second peak. Although the RSM turbulence model is more detailed than those of the general turbulence models (and even though it includes all Reynolds stress transport for all directions), it still cannot capture relaminarization due to a lack of separation-induced transition correction. However, the second peak Nusselt number in RSM-LPS and RSM-QPS can still be observed. For the pressure-strain terms in equations (17(17) ${\psi }_{ij}={\psi }_{ij,1}+{\psi }_{ij,2}+{\psi }_{ij,w}$(17) ) and (21(21) $\begin{array}{rl}{\psi }_{ij}& =-(3.4\rho \epsilon +1.8P){\tau }_{ij}\\ & +4.2\rho \epsilon \left({\tau }_{ik}{\tau }_{kj}-\frac{1}{3}{\tau }_{mn}{\tau }_{mn}{\delta }_{ij}\right)\\ & +\left(0.8-1.3\sqrt{{\tau }_{ij}{\tau }_{ij}}\right)\rho k{S}_{ij}\\ & +1.25\rho k\left({\tau }_{ik}{S}_{jk}+{\tau }_{jk}{S}_{ik}-\frac{2}{3}{\tau }_{mn}{S}_{mn}{\delta }_{ij}\right)\\ & +0.4\rho k({\tau }_{ik}{\mathrm{\Omega }}_{jk}+{\tau }_{jk}{\mathrm{\Omega }}_{ik})\end{array}$(21) ), they are $\epsilon$-based Reynolds stress models with a slow pressure-strain term (also known as the return-to-isotropy), a rapid pressure-strain term, and a wall-reflection term. This relationship may not precisely explain the direct effect on the second peak, but it is expected to calculate the trend in the magnitude of Nusselt number. As shown in equations (21(21) $\begin{array}{rl}{\psi }_{ij}& =-(3.4\rho \epsilon +1.8P){\tau }_{ij}\\ & +4.2\rho \epsilon \left({\tau }_{ik}{\tau }_{kj}-\frac{1}{3}{\tau }_{mn}{\tau }_{mn}{\delta }_{ij}\right)\\ & +\left(0.8-1.3\sqrt{{\tau }_{ij}{\tau }_{ij}}\right)\rho k{S}_{ij}\\ & +1.25\rho k\left({\tau }_{ik}{S}_{jk}+{\tau }_{jk}{S}_{ik}-\frac{2}{3}{\tau }_{mn}{S}_{mn}{\delta }_{ij}\right)\\ & +0.4\rho k({\tau }_{ik}{\mathrm{\Omega }}_{jk}+{\tau }_{jk}{\mathrm{\Omega }}_{ik})\end{array}$(21) )–(23(23) $\begin{array}{rl}{\mathrm{\Omega }}_{ij}& =\frac{1}{2}\left(\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}-\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}\right)\end{array}$(23) ), RSM-QPS has more complex equations than RSM-LPS, most notably the Reynolds stress anisotropy tensor, mean strain rate, and rotation rate. These equations are related to the logarithmic region of a turbulent boundary layer, and we may be able to obtain more accurate results if we improve this calculation algorithm or develop a new unique correlation approach in the boundary layer region. Finally, we still believe that the RSM of these models can help future researchers improve their performance.

Figure. 5. Comparison of Nusselt numbers between experiments.

The v²-f model is like the standard k-ε model; however, it includes near-wall turbulence anisotropy and nonlocal pressure-strain effects. It is a general low-Reynolds number turbulence model valid up to solid walls and thus does not require wall functions. Although the model was designed to simulate attached or mildly separated boundary layers, it also accurately simulates separation-dominated flows. The v²-f model is distinguished by its use of the velocity scale rather than TKE to calculate eddy viscosity, defined as the velocity fluctuation normal to the streamlines by the variable $f$, and is the solution to the elliptic relaxation equation. In this case, the v²-f model uses an elliptical operator to compute a term like the RSM’s pressure strain correlation. Ellipticity is defined by a modified Helmholtz operator that introduces wall effects via a linear differential equation. From the above explanation, it is advantageous to discover that the near-wall relation has more chance for state-of-the-art improvement to separation and relaminarization behavior.

4.2. Velocity profiles

The axial velocity component profiles of the impinging jet (onto a constant hot plate H/D = 2.0, Re = 23,000) are shown in Figure 6 and compared to experimental data from Cooper (1993). The results show that the effect of near-wall velocity was consistent with all turbulence models in the range of 0.5 $\le r/D\le$ 1.0. Furthermore, the error increases significantly as the $r/D$ ratio increases. Interestingly, Figure 6(c) at $r/d$ = 1.5 depicts the velocity profiles at the local Nusselt number minimum in the negative first peak that all the turbulence models are not capable of predicting, especially RSM. When we look more closely at the RSM turbulence model, it is discovered that the $\epsilon$-based models could predict a slight trend of local Nusselt second peak trends faster than experimental data. Comparing the resulting velocity to the experiment, the RSM $\epsilon$-based (RSM-LPS, RSM-QPS) model had consistently less error than the $\omega$-based (RSM-SO, RSM-S-BSL) approach. As a result, these findings clearly show that the prediction of the second peak may be related to the $\epsilon$-transport equation; it is highly recommended that to improve the accuracy of the RSM inner sublayer, we need some mathematical modifications, such as an inner new $\epsilon$-transport equation or a new correlation of the internal sublayer (Panda, 2019). In addition, the stress-omega model, which is $\omega$-based stress-transport model related to the $\omega$-transport equations and does not require any wall reflection treatments. The stress-$\omega$ model eliminates the free-stream sensitivity that is examined by the baseline (BSL) model. The near-wall model equation was integrated into both models. As a result, it may be difficult to predict the occurrence of the Nusselt second peak behavior. In contrast, in Figure 6(e–f), the velocity profiles negatively deviate more from the RSM, and v²-f is overpredicted while moving away from the jet core region. We suspect that the energy dissipation of turbulence is relatively high. However, considering the increased $X/D$ in vertical distance, v²-f still predicted axial velocity well compared to the experiment. It can be concluded that v²-f and all RSM turbulence models are not practical for predicting inside near-wall boundary velocity, although v²-f has improved near-wall turbulence anisotropy and nonlocal pressure-strain effects. Nevertheless, the RSM has also modified the wall function but is still not as accurate as we expected. This finding is like the numerical results of Zuckerman et al. (2006)

Figure 6. Comparisons of simulated and measured mean velocity with H/D = 2 and Re = 23,000: (a) r/D = 0.5, (b) r/D = 1.0, (c) r/D = 1.5, (d) r/D = 2.0, (e) r/D = 2.5, and (f) r/D = 3.0.

When $r/D$ = 2.0 is indicated as the second peak of the local Nusselt number described in Figure 6(d), we discovered that the v²-f model accurately predicts velocity profiles in all regimes. However, all RSM turbulence models fail to predict the entire range. In the detailed analyses, the $\epsilon$-based approach can detect the second local Nusselt number faster than experimental data, but it has more error in the velocity profiles than the $\omega$-based approach. We can conclude that the faster second peak causes increasingly faster turbulent energy dissipation, resulting in low momentum in the entire region. The v²-f model uses a velocity variance scale and an elliptic relaxation function to detect the relaminarization behavior in the region of separation-induced relaminarization.

4.3. Turbulence intensity and turbulence kinetic energy profiles

Since the flow of the jet impingement is anisotropic, in Figure 7(a–d), the mentioned parameters are shown for turbulence intensity (TI) and turbulence kinetic energy (TKE) for the cases $H/D$ = 2.0 and $r/D$ = 1.5 and 2.0. Physically, the TI or unresolved unsteadiness is calculated from the mean square fluctuation components divided by the mean absolute velocity, while TKE is characterized by measured root-mean-square velocity fluctuations or can be calculated based on the closure method in the RANS equations. Generally, TKE can be produced by fluid shear, friction, or integral scale and then transferred down the turbulence energy cascade and dissipated by viscous forces. This production, transport, and dissipation process can be expressed as separation-induced relaminarization (Ting, 2016). Therefore, if we find that the Nusselt number and velocity profiles of v²-f can produce the best results compared to experimental data, we can assume that in this study, v²-f is the best turbulence model even before comparison with the RSM; hence, the results in this section will focus on v²-f computational data to represent the baseline of accuracy (Issac et al. (2020); Zuckerman et al (2006)).

Figure 7. Turbulence intensity and turbulence kinetic energy distributions before and after the second peak: (a) TI at r/D = 1.5, (b) TI at r/D = 2.0, (c) TKE at r/D = 1.5, and (d) TKE at r/D = 2.0.

Considering all RSM turbulence models in terms of TI and TKE estimation the RSM-LPS and RSM-QPS values are nearly identical across all regions, RSM-S-BSL values are in the same order if $X/D\le 0.1$, and RSM-SO is a significant blunder across all regimes. Despite these challenges, when we move from the jet core to the second peak region ($r/D$ = 2.0), we suspect a very complex mechanism due to relaminarization behavior is occurring. However, we believe that the Nusselt number and average velocity profile prediction results from v²-f model provides the best mapping with an experiment in the relaminarization region. Hence, we hypothesize that TI and TKE must be good enough to be a key for this study in terms of analyses for another turbulence model. When considered close to the wall, the $\epsilon$-based method is superior to the $\omega$-based method. As a result, we can infer that the RSM $\omega$-based approach is not suitable for separation-induced relaminarization processes; if it is used in the future, more modification will be needed to improve the near-wall region. However, the TI is often referred to as the turbulence level, which is correlated with the TKE. Here, we can show that the improvement may modify the TKE equation, TI equation, or a related equation (Khalaji et al., 2016; Huang et al., 2018; Zhang et al., 2019).

4.4 Velocity and nusselt contour

The predicted velocity field with H/D = 2.0 is depicted in Figure 8. Once the grid is independent and fine enough with $y^{+}<1$, the velocity vector described in Figure 8(a) shows parabolic profiles in the inner layer of the boundary layer. Next, for the local Nusselt number contour of v²-f, as depicted in Figure 8(b), given the details at the center of jet impingement with the so-called stagnation point, it increases the local Nusselt number and decreases monotonically at r/D = 2.0, which produces a simulated result in the second peak. These results are consistent with the mechanism described in other research (Zhang et al., 2019; Ekkad & Singh, 2021), giving us confidence that all our settings are correct and that we have obtained numerical fidelity results. The region of highest TKE immediately adjacent to the wall corresponds to the part of the secondary peak in the Nusselt number; so, this peak can be attributed to the development of turbulent flow in the wall jet (Zuckerman et al., 2006). Figure 8(c) shows a region where the deceleration of the flow and decay of the wall-jet boundary layer caused a reduction in turbulence and a decrease in the Nusselt number. The other anisotropic turbulence models could not predict the secondary peak because they produced different (and likely incorrect) velocity fields, directly influencing the thermal energy transfer calculations. The four-equation v²-f model resolves the complete boundary layer and includes an elliptic relaxation function, which tracks the damping effects of walls on the velocity variance scale. This finding emphasizes that if we need greater accuracy regarding the Nusselt number, we should not only improve the elaborate anisotropy transport equation model but also the near-wall modeling must be considered, such as the $\epsilon$-based approach and the $\omega$-based approach.

Figure 8. Results of characteristic flow distributions for H/D = 2 with Re = 23,000 of v²-f: (a) velocity vector, (b) local Nusselt number, and (c) velocity contour.

4.5. The effect of the Reynolds Number

For comparison, the results shown in Figure 9 focus only on the v²-f turbulence model to investigate the effect of the Reynolds number. The formation of the secondary peak of the Nusselt number can be captured by the v²-f turbulence model for all Reynolds numbers in this study. The location of the secondary peak removes downward in the radial direction as the Reynolds number increases. Based on these findings, we hypothesize that the secondary peak is caused by the spreading wall jet’s transition from a laminar to a turbulent boundary layer flow. Furthermore, as the Reynolds number increases, the effect of the toroidal flow phenomenon, which forms in the shear region around the circumference of the jet, striking the impinging plate appears to be enhanced (Issac et al., 2020).

Figure 9. Variation of local Nusselt numbers for H/D = 2: (a) Re = 30,000, (b) Re = 23,000 (c) Re = 15,000, and (d) Re = 10,000.

5. Conclusions

Numerical study of jet impingement using a variety of relevant anisotropic turbulence models were investigated and compared those of results with available experimental data. Surprisingly, better agreement with experimental data was obtained with the turbulence models based on some near-wall turbulence anisotropy and nonlocal pressure-strain effects. The advantage of the v²-f model is using the velocity scale instead of the turbulence kinetic energy to evaluate the eddy viscosity. The anisotropic wall effects are modeled through the elliptic relaxation function by solving an additional elliptic equation of the Helmholtz type. Consequently, the secondary peak in the Nusselt number dependency on a radius could be predicted with reasonable accuracy compared to the second-order closure Reynolds stress turbulence models. The secondary peak location was found at r/D = 2.05 with the experiment and produced similar results to other study. At H/D = 2.0, the golden ratio that clearly shows the secondary peak of the Nusselt number decreases in the radial direction between the stagnation point; at a radial distance of approximately r/D = 2.05, a second peak is observed. Initially, the first peak of Nusselt number prediction is at r/D = 0.5 due to the change in the radial flow velocity as it travels radially downstream from the stagnation point. When the jet deflects after hitting the surface, it accelerates rapidly but at a greater distance, as the spreading of the jet causes its velocity to diminish. Finally, we found that the transport equation based on v²-f produces better results for all flow characteristics compared to the RSM in all approaches. The critical point for turbulence modeling improvement is identified by focusing on the near-wall transport theorem. However, this modification or improvement is challenging and may require considerably more relevant referenced tests.

Disclosure statement

No potential conflict of interest was reported by the author(s).