Unsteady Numerical Calculation of Oblique Submerged Jet

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.


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][2][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][9][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 Energies 2020, 13, 4728 3 of 13

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 L w = 55 D, a width of W w = 7 D, and a height of H w = 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 −4 s, and the total simulation time was 1 s.
Energies 2020, 13, x FOR PEER REVIEW 3 of 14 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 −4 s, and the total simulation time was 1 s. 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 ro1l 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 ro1l 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 × 10 6 , and the maximum number of mesh cells was 4.80 × 10 6 .  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: 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 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 o 2 xyz 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 o 2 xyz 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 × 10 6 , and the maximum number of mesh cells was 4.80 × 10 6 .  Figure 1 shows the comparison between the CFD and the experimental results of the axial velocity V/V b in the centerline of the jet. V b denotes the bulk velocity at the jet exit, which can be defined as: 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.   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.  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. 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: The turbulent eddy viscosity: where ρ is the density. S takes on the usual definition for mean strain: Wall blocking is accounted for by the damping function: where χ = R/ν and ν = µ/ρ. The switching function is: where d is the minimum distance to the nearest wall.

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/V max with the empirical formula for a fully developed circular jet. The calculation formula is as follows: where V max represents the maximum velocity at the jet exit, r represents the radial direction, and n is the empirical constant for the power law.
Energies 2020, 13, x FOR PEER REVIEW 6 of 14 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: 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. 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.  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/V b 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/V b 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.  Figure 5 shows the V/V b 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.

Analysis of the Oblique Submerged Jet Flow Field under Different Times
Energies 2020, 13, x FOR PEER REVIEW 7 of 14 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 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 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 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 7, under the condition of Re = 11,700, when 0 s ≤ t ≤ 0.05 s, the value of V/V b 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/V b 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/V b to peak and the time for reaching the state where the V/V b 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 Energies 2020, 13, 4728 9 of 13 value and then decreased slightly. With the increase of Re, the peak pressure on the axis increased, but its distribution was similar.    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 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 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 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 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 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 C pmax increased slowly, then rapidly increased, then slowly increased, and finally remained unchanged. With the increase of the Reynolds number, the change rule of the C pmax value with time was similar, but the increasing rate was larger, the maximum value of C pmax decreased slightly after remaining stable, but the time required to reach the maximum value decreased.

CPmax
Re=11700 Re=23400 Re=35100 t Figure 11. Variation of Cpmax with t at different Reynolds numbers.

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: (1) The jet fluid diffused gradually after leaving the nozzle, the ambient fluid was continuously Figure 11. Variation of C pmax with t at different Reynolds numbers.

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: (1) The jet fluid diffused gradually after leaving the nozzle, the ambient fluid was continuously sucked into the jet, and there were symmetrical vortices around the main body of the jet; when impinging on the plate, most of the jet fluid deflected in the forward flow direction, and the vortex range increased in the forward flow direction, while decreasing in the backward flow direction. (2) In the process of jet movement, there was a peak value of velocity on the axis, and the pressure at the peak value was the minimum value; with the increase of time, after a certain time, the velocity, flow angle, and pressure distribution along the axis remained unchanged: 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. (3) The variation of the maximum impact pressure coefficient on the impingement plate with time was affected by the Reynolds number, but its distribution law was the same.