*Article* **Unsteady Numerical Calculation of Oblique Submerged Jet**

#### **Weixuan Jiao 1, Di Zhang 1, Chuan Wang 1,\*, Li Cheng 1,\* and Tao Wang <sup>2</sup>**


Received: 31 July 2020; Accepted: 8 September 2020; Published: 11 September 2020

**Abstract:** A water jet is a kind of high-speed dynamic fluid with high energy, which is widely used in the engineering field. In order to analyze the characteristics of the flow field and the change of law of the bottom impact pressure of the oblique submerged impinging jet at different times, its unsteady characteristics at different Reynolds numbers were studied by using the Wray–Agarwal (W-A) turbulence model. It can be seen from the results that in the process of jet movement, the pressure at the peak of velocity on the axis was the smallest, and the velocity, flow angle, and pressure distribution remain unchanged after a certain time. In the free jet region, the velocity, flow angle, and pressure remained unchanged. In the impingement region, the velocity and flow angle decreased rapidly, while the pressure increased rapidly. The maximum pressure coefficient of the impingement plate changed with time and was affected by the Reynolds number, but the distribution trend remained the same. In this paper, the characteristics of the flow field and the law of the impact pressure changing with time are described.

**Keywords:** unsteady flow; numerical simulation; submerged jet

#### **1. Introduction**

The main advantages of an impinging jet are that it can effectively strengthen local heat and mass transfer, and thus has been widely used in many engineering fields [1–3] such as water conservancy and hydropower discharge engineering, in which the plunging nappe impacts the plunge pool, the cooling of electronic components, the chemical mixing jet, the micro-sprinkler irrigation in water-saving irrigation, etc. The impinging jet depends on the shape of the nozzle, the jet outlet speed, the impinging angle, the impinging height, and other parameters. Previous studies have shown that the flow field of the impinging jet mainly includes a free jet zone, an impingement zone, and a wall jet zone [4,5]. Moreover, with the change of parameters, the flow field in different regions presents different laws. Compared with a vertical impinging jet, the flow field and pressure distribution of the oblique impinging jet in the wall jet region are not uniform and are more complex [6,7].

With the improvement of computer performance, computational fluid dynamics technology provides a low-cost and high-efficiency solution to the complex problems in water conservancy engineering [8–10], municipal engineering [11,12], and other practical projects. Many scholars use Computational Fluid Dynamics (CFD) technology to study the flow field of an impinging water jet. Ghiti et al. [13] studied the jet impingement on a horizontal plate at different Reynolds numbers. The results showed that there is a positive correlation between vorticity and the Reynolds number. The turbulent characteristics near the wall plate are greatly affected by the Reynolds number. Based on the large eddy simulation (LES) turbulence technique, So et al. [14] carried out three-dimensional

analysis of the flow mechanism of the unconfined plane impinging jet and provided a detailed flow structure. Gopalakrishnan et al. [15] studied round jet impingement at multiple impingement angles. The results showed that the larger the jet angle, the higher the jet growth that was predicted. Baghel et al. [16] presented the hydrodynamic structure and heat transfer performance of a vertical free surface jet impinging on a horizontal plate. By comparing the simulation results of the free surface jet velocity field with different turbulence models, it was found that the v2f turbulence model can predict the flow field better than other turbulence models. There is also much research on the impact pressure of a water impinging jet. Wei et al. [17] studied the dynamic pressure caused by a submerged inclined jet by forced ventilation at the bottom of a deep plunge pool in detail, and proposed that the combination of a reduced jet thickness and an enhanced jet aeration can effectively reduce the impact on the bottom of a deep plunge pool. Tian et al. [18] investigated the impact pressure characteristics of jet equipment that can reach a maximum jet velocity of 50 m/s. Through numerical simulation, the reason for a scale effect of the impact pressure was elucidated. Shao et al. [19] improved the traditional smoothed particle hydrodynamics (SPH) method and verified its effectiveness by simulating the water jet impingement problem. The results showed that the improved SPH method can accurately simulate the shape and pressure distribution of an impinging jet. Based on the W-A turbulence model, Wang et al. [20] used CFD technology to study the flow characteristics of a submerged impinging water jet at different impact heights. The results show that there is a negative correlation between impact pressure, impact pressure coefficient and impact height. Through CFD technology, scholars can also carry out research on the unsteady flow process of water jets. Based on the incompressible Reynolds-averaged Navier–Stokes (RANS) equations with the W-A turbulence model, Zhang et al. [21] carried out numerical simulation of single-jet impingement. The results showed that the W-A model has good applicability in numerical simulation of single-jet impingement. Zhang et al. [22] studied the time-dependent initial flow structure of a subsonic unsteady elliptic jet using the large eddy simulation (LES) method. Based on the unsteady FLUENT(ANSYS Inc., Pittsburgh, PA, USA) numerical simulation technology, Yang et al. [23] carried out numerical simulation on the internal flow characteristics of different impinging jets. The results showed that the local Nusselt number oscillates over time. Chung et al. [24] studied the momentum and heat transfer characteristics of an unsteady impinging jet by using the numerical simulation method. The results showed that the impact heat transfer is very unstable due to the influence of the primary vortex generated by the jet nozzle. Zhang et al. [25] carried out numerical simulation of the internal flow field of the self-excited pulsed cavitating jet nozzle, analyzed the evolution process of the initial generation, energy concentration, and release of the cavitation in the chamber in one cycle, and explained the occurrence mechanism of the cavitating jet. Zaafori et al. [26] used the finite difference method to analyze the influence of the initial perturbation of the co-current flow on the jet dynamics and thermal behavior over time.

An oblique submerged jet is a complex flow with a practical engineering application background and important theoretical research value. The systematic study of it not only helps to deepen the understanding of the interaction mechanism between jet and wall, but also provides a scientific basis for engineering applications. In order to reveal the evolution law of the hydraulic characteristics of an oblique submerged jet flow field with time and the influence of the Reynolds number on its flow characteristics, the UDF (User Define Function) of FLUENT software was used to introduce the W-A model to simulate the unsteady flow field and pressure distribution of an oblique submerged impinging jet under different Reynolds numbers in this paper. In this study, a general circular jet, i.e., a fully developed circular jet, was selected for unsteady calculation.

#### **2. Numerical Simulation**

#### *2.1. Geometric Model and Boundary Conditions*

A three-dimensional non-closed oblique submerged impinging jet model was built, as shown in Figure 1. The jet was obliquely sprayed into a regular sink through a circular nozzle. The bottom surface was an impingement plate, the top surface was a free surface, the left and right sides were outlet boundaries, and the front and back sides were fixed wall boundaries. The inner diameter of the nozzle was *D* = 20 mm, and the outer diameter was *d* = 22 mm. All lengths in this paper were dimensionless with *D*; in order to ensure that the jet outlet flow was a fully developed turbulent pipe flow, the nozzle length was 50 *D*; the calculation domain of the sink was a length of *Lw* = 55 *D*, a width of *Ww* = 7 *D*, and a height of *Hw* = 12 *D*; the vertical distance between the center of the nozzle outlet and the impingement plate was *H* = 3 *D*. The calculation time step was 10−<sup>4</sup> s, and the total simulation time was 1 s.

**Figure 1.** 3-D oblique submerged jet.

In order to better analyze the flow field of the oblique submerged impinging jet, as shown in Figure 1, the oblique submerged impinging jet was described in a double coordinate system. The rectangular coordinate system *ro*1*l* was used to describe the impinging jet flow in the normal central plane before the impingement. The rectangular coordinate system *o2xyz* was used to describe the impinging jet flow in the normal central plane after the impingement. The origin of *ro*1*l* was set at the center of the nozzle outlet; r and l represent the radial and axial (along the jet) directions, respectively. The *o2xyz* coincides with the geometric center (GC) of the jet on the impingement plane, with x, y, and z representing the longitudinal, transverse, and vertical directions, respectively.

The grid quality has a crucial influence on the numerical calculation. In this paper, ICEM software was used to divide the calculation domain into hexahedral grids. A grid sensitivity analysis was carried out to assess the required grid density. As shown in Table 1, several cases with different mesh cell numbers were compared. The minimum number of mesh cells was 2.65 <sup>×</sup> 106, and the maximum number of mesh cells was 4.80 <sup>×</sup> 106.

**Table 1.** Grid information of the different cases.


Figure 1 shows the comparison between the CFD and the experimental results of the axial velocity *V*/*Vb* in the centerline of the jet. *Vb* denotes the bulk velocity at the jet exit, which can be defined as:

$$V\_b = 4\text{Q/}\pi\text{D}^2\tag{1}$$

where *V* represents the velocity at any position of the jet exit; Q is the flow rate of the jet.

As shown in Figure 2, when the number of mesh cells of the overall calculation domain reached 3.73 million, the difference between the CFD numerical simulation results and the experimental results was very small. After grid sensitivity analysis, the total number of mesh cells in the calculation domain was finally determined to be 4.26 million, as shown in Figure 3. The drawn grid can effectively reduce the false diffusion caused by poor grid quality, and can also improve the computing speed.

**Figure 2.** (**a**) Definition of the bulk mean velocity (*Vb*); (**b**) grid independence test.

**Figure 3.** Grid of the oblique submerged impinging jet: (**a**) Computational domain; (**b**) regular sink.

In the calculation, based on the pressure solver, the SIMPLE algorithm for pressure–velocity coupling was adopted. The pressure equation was discretized by the second-order scheme, and the momentum and turbulence model equations were discretized by the second-order upwind scheme. The turbulent viscosity ratio was used to define the turbulence on the boundary of the flow field. The inlet turbulent viscosity ratio was 10. The inlet adopted velocity inlet, the speed was evenly distributed, and the inlet working condition is shown in Table 2. The outlet was set as a pressure outlet, and the gauge pressure value was 0 Pa. This paper simulated the three-dimensional impinging jet in an open channel water environment, and its free surface remained almost unchanged, so the rigid-lid (RL) assumption method was used to process the free surface, and it was set as a symmetrical surface. A fixed wall with no slip was selected for the wall surface.



#### *2.2. Wray–Agarwal Turbulence Model*

#### 2.2.1. Introduction of Wray–Agarwal Turbulence Model

The Wray–Agarwal turbulence model is a low Reynolds number one-equation model developed using the *k-*ω model. In this paper, the maximum y+ was less than 1 to ensure that the near wall treatment for the W-A model was accurate. Wang et al. [20] and Han et al. [27] carried out a comparison between numerical simulations and experimented on various turbulence models, including the W-A model. The results showed that the W-A model simulates the separation characteristics of the boundary layer more accurately and has good applicability in a numerical simulation of a submerged jet. The W-A model includes the cross-diffusion in the ω equation, and *R* = κ/ω is determined by the following transport equation:

$$\frac{\partial R}{\partial \mathbf{t}} + \frac{\partial u\_{j}R}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} \Big[ (\sigma\_{R}R + \nu) \frac{\partial R}{\partial \mathbf{x}\_{j}} \Big] + C\_{1}RS + f\_{1}C\_{2\text{kio}}\frac{R}{S} \frac{\partial R}{\partial \mathbf{x}\_{j}} \frac{\partial S}{\partial \mathbf{x}\_{j}} - (1 - f\_{1})C\_{2\text{k}c}R^{2} \left( \frac{\frac{\partial S}{\partial \mathbf{x}\_{j}} \frac{\partial S}{\partial \mathbf{x}\_{j}}}{S^{2}} \right) \tag{2}$$

The turbulent eddy viscosity:

$$
\mu = \rho f\_u \mathbb{R} \tag{3}
$$

where ρ is the density. *S* takes on the usual definition for mean strain:

$$S = \sqrt{2S\_{ij}S\_{ij}} \ S\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{4}$$

Wall blocking is accounted for by the damping function:

$$f\_{\mu} = \frac{\chi^3}{\chi^3 + \mathcal{C}\_W^3},$$

where χ = *R*/ν and ν = μ/ρ. The switching function is:

$$\mathcal{G}\_1 = \min\{\tanh(\arg\_1^4), 0.9\} \\ \arg\_1 = \frac{1 + \frac{d\sqrt{\mathbb{R}\mathbb{S}}}{\nu}}{1 + \left[\frac{\max\{d\sqrt{\mathbb{R}\mathbb{S}}, 1.5\mathbb{R}\}}{20\nu}\right]^2} \tag{6}$$

where *d* is the minimum distance to the nearest wall.

The constants are:

$$\begin{array}{l} \text{C}\_{1k\nu} = 0.0829 \,\text{C}\_{1k\varepsilon} = 0.1127\\ \text{C}\_{1} = f\_{1}(\text{C}\_{1k\omega} - \text{C}\_{1k\varepsilon}) + \text{C}\_{1k\varepsilon} \\ \sigma\_{k\omega} = 0.72 \,\sigma\_{k\varepsilon} = 1.0 \\ \sigma\_{R} = f\_{1}(\sigma\_{k\omega} - \sigma\_{k\varepsilon}) + \sigma\_{k\varepsilon} \\ \text{C}\_{2k\omega} = \frac{\text{C}\_{1k\omega}}{\text{x}^{2}} + \sigma\_{k\omega} \,\text{C}\_{2k\varepsilon} = \frac{\text{C}\_{1k\varepsilon}}{\text{x}^{2}} + \sigma\_{k\varepsilon} \\ \kappa = 0.41 \,\text{C}\_{\omega\nu} = 8.54 \end{array} \tag{7}$$

#### 2.2.2. Verification of the Wray–Agarwal Turbulence Model

The Wray–Agarwal, Standard *k-*ε, RNG *k-*ε (Renormalization Group *k-*ε), Realizable *k-*ε, Standard *k-*ω, and SST *k-*ω (Shear Stress Transfer *k-*ω) turbulence models were used to simulate the oblique submerged impinging jet with an impact angle of θ = 45◦ and an impact height of H/D = 3. The numerical results were compared with the experimental data of PIV [28].

Figure 4 shows a comparison between the numerical velocity *V*/*Vmax* with the empirical formula for a fully developed circular jet. The calculation formula is as follows:

$$\text{V/V}\_{\text{max}} = (1 - 2\text{r/D})^{1/\text{n}} \tag{8}$$

where *Vmax* represents the maximum velocity at the jet exit, r represents the radial direction, and n is the empirical constant for the power law.

**Figure 4.** Comparison of calculation results with different turbulence models: Profiles of the normalized axial velocity *V*/*Vmax* near the jet exit (l/D = 0.5).

When the value of n was 7, the velocity distribution calculated by Formula (8) was approximately consistent with the velocity distribution of a fully developed circular jet [29]. As shown in Figure 4, the numerical calculation results with different turbulence models are in good agreement with the empirical formula. This also shows that the flow at the jet exit was a fully developed jet flow. At the same time, the calculated results of the W-A model are closest to the empirical formula.

Table 3 illustrates the percentage difference between the CFD and the experimental results of the axial velocity *V*/*Vb* in the centerline of the jet. The validity of the turbulence model was verified by comparing the numerical simulation and experimental results of the axial velocity *V*/*Vb* in the centerline of the jet. As can be seen from Table 3, the calculated results of the W-A model are closest to the experimental results. Hence, the W-A model can accurately simulate the unsteady characteristics of an oblique submerged impinging jet.

**Table 3.** Numerical and experimental results with different turbulent models.


#### **3. Results**

#### *3.1. Analysis of the Oblique Submerged Jet Flow Field under Di*ff*erent Times*

Figure 5 shows the *V*/*Vb* velocity contours and streamline (*Re* = 35,100) of the mid-section (*y*/*D* = 0) of the jet at different times (0.01~1 s). When *t* = 0.01 s, the jet presented a circular arc shape distribution; when *t* = 0.02 s, the jet fluid moved axially, there was a vortex around the front of the jet, and the vortex had good symmetry with respect to the central axis; with the flow of the jet, the core position of the vortex gradually approached the impingement plate. When *t* = 0.1 s, most of the jets were deflected in the forward flow direction, the vortex in the forward flow direction gradually increased, and the vortex in the backward flow direction gradually decreased. With the increase of time *t*, the jet flowed along the wall, the vortex core remained near the front of the wall jet region in the forward flow direction, and the vortex gradually increased.

**Figure 5.** Normalized mean velocity *V*/*Vb* contour and streamline of the mid-section of the oblique submerged impinging jet at different times: (**a**) *t* = 0.01 s, (**b**) *t* = 0.02 s, (**c**) *t* = 0.05 s, (**d**) *t* = 0.06 s, (**e**) *t* = 0.1 s, (**f**) *t* = 0.2 s, (**g**) *t* = 0.5 s, (**h**) *t* = 1.0 s.

Figure 6 shows the vorticity contours of the mid-section of the oblique submerged impinging jet at different times (*Re* = 35,100). It can be seen that in the stage where the jet ejected the nozzle until it hit the impingement plate (0 ≤ *t* ≤ 0.5 s), the vorticity intensity was small. When the jet fluid collided with the impingement plate, the vorticity intensity increased in the backward flow direction. When comparing with Figure 4, it can be observed that when the jet deflected and flowed along the wall jet region, the vortex mainly occurred in the tongue part before the jet, and the vortex intensity was large. With the increase of time, the flow of the tongue part before the jet became complicated, and the vorticity distribution was irregular.

**Figure 6.** Vorticity of mid-section of the oblique submerged impinging jet at different times: (**a**) *t* = 0.01 s, (**b**) *t* = 0.02 s, (**c**) *t* = 0.05 s, (**d**) *t* = 0.06 s, (**e**) *t* = 0.1 s, (**f**) *t* = 0.2 s, (**g**) *t* = 0.5 s, (**h**) *t* = 1.0 s.

#### *3.2. Comparison of the Oblique Submerged Jet Flow Field under Di*ff*erent Reynolds Numbers*

Figures 7 and 8 show the dimensionless velocity *V*/*Vb* and flow angle ϕ distribution on the axis of the submerged jet at different times under different Reynolds numbers. As shown in Figure 7, under the condition of *Re* = 11,700, when 0 s ≤ *t* ≤ 0.05 s, the value of *V*/*Vb* remained unchanged, then decreased rapidly, and the influence range of jet increased; when *t* = 0.1 s, the jet velocity first remained constant, then rose, and finally dropped rapidly; when *t* ≥ 0.5 s, the distribution of the *V*/*Vb* value was as follows: In the region of 0 ≤ *l*/*D* ≤ 3, it remained unchanged, and in the region of *l*/*D* ≥ 3, it decreased rapidly. When *Re* was increased, that is, when the speed of the jet nozzle was increased, the time for *V*/*Vb* to peak and the time for reaching the state where the *V*/*Vb* value did not change became shorter. It can be seen from Figure 8 that with the increase of time, the flow angle in the free jet region gradually increased to the impinging angle and then remained unchanged, and the decrease rate of the flow angle in the impingement region increased.

Figure 9 shows the pressure distribution on the axis of the submerged jet at different times under different Reynolds numbers. As shown in the figure, under the condition of *Re* = 11,700, when *t* = 0.01 s, the pressure of the jet along the axis line increased rapidly from a negative value to a positive value and then decreased to 0; with the increase of *t*, the pressure was maintained at the value of 0, then decreased, then increased, and finally decreased to 0; compared with Figure 7a, it can be seen that the pressure was the smallest when the velocity was maximum. When the jet velocity distribution remained unchanged, the pressure distribution along the axis of the jet also remained unchanged: The pressure in the free jet region remained near 0, and the pressure in the impingement region increased rapidly to the maximum

value and then decreased slightly. With the increase of *Re*, the peak pressure on the axis increased, but its distribution was similar.

**Figure 7.** Normalized mean axial velocity *V*/*Vb* along the jet centerline. (**a**) *Re* = 11,700, (**b**) *Re* = 23,400, (**c**) *Re* = 35,100.

**Figure 8.** *Cont.*

**Figure 8.** Definition of flow angle ϕ along the jet centerline. (**a**) *Re* = 11,700, (**b**) *Re* = 23,400, (**c**) *Re* = 35,100.

**Figure 9.** Definition of pressure along the jet centerline. (**a**) *Re* = 11,700, (**b**) *Re* = 23,400, (**c**) *Re* = 35,100.

*3.3. Pressure Distribution on the Impingement Plate at Di*ff*erent Times*

Figure 10 shows the pressure contours (*Re* = 23,400) of the jet impinging on the plate (*oxy* plane) at different times. It can be seen from the figure that when *t* = 0.06 s, the pressure on the impingement plate was almost 0; when *t* = 0.07 s, the pressure on the impingement plate was elliptical and had good symmetry; when *t* = 0.08 s, the pressure increased near the impinging origin. With the increase of time, the maximum pressure on the impingement plate increased gradually. When *t* ≥ 0.5 s, the pressure distribution near the impinging origin was similar.

**Figure 10.** Pressure contours of jet impinging on the plate (*oxy* plane) at different times: (**a**) *t* = 0.06 s, (**b**) *t* = 0.07 s, (**c**) *t* = 0.08 s, (**d**) *t* = 0.09 s, (**e**) *t* = 0.10 s, (**f**) *t* = 0.15 s, (**g**) *t* = 0.20 s, (**h**) *t* = 0.20 s, (**i**) *t* = 1.00 s.

Figure 11 shows the change in the maximum pressure coefficient on the impingement plate with time under different Reynolds numbers. It can be seen from the figure that when *Re* = 11,700, with the increase of time, the maximum pressure coefficient *Cpmax* increased slowly, then rapidly increased, then slowly increased, and finally remained unchanged. With the increase of the Reynolds number, the change rule of the *Cpmax* value with time was similar, but the increasing rate was larger, the maximum value of *Cpmax* decreased slightly after remaining stable, but the time required to reach the maximum value decreased.

**Figure 11.** Variation of *Cpmax* with *t* at different Reynolds numbers.

#### **4. Conclusions**

Unsteady numerical calculation of an oblique submerged impinging jet was carried out by using the W-A model, and the following conclusions were obtained:


**Author Contributions:** Data curation, C.W. and L.C.; formal analysis, W.J. and D.Z.; methodology, L.C. and W.J.; writing—original draft, W.J. and D.Z.; writing—review and editing, T.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 51979240 and 51609105), the Priority Academic Program Development of Jiangsu Higher Education Institutions (Grant number PAPD), the Open Research Subject of Key Laboratory of Fluid and Power Machinery (Xihua University), and the Ministry of Education (szjj2016-061).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
