Numerical Investigation of a High-Pressure Submerged Jet Using a Cavitation Model Considering Effects of Shear Stress

Abstract

In the current research, a high-pressure submerged cavitation jet is investigated numerically. A cavitation model is created considering the effect of shear stress on cavitation formation. As such, this model is developed to predict the cavitation jet, and then the numerical results are validated by high-speed photography experiment. The turbulence viscosity of the renormalization group (RNG) k-ε turbulence model is used to provide a flow field for the cavitation model. Furthermore, this model is modified using a filter-based density correction model (FBDCM). The characteristics of the convergent-divergent cavitation nozzle are investigated in detail using the current CFD simulation method. It is found that shear stress plays an important role in the cavitation formation in the high-pressure submerged jet. In the result predicted by the Zwart-Gerber-Belamri (ZGB) cavitation model, where critical static pressure is used for the threshold of cavitation inception, the cavitation bubble only appears at the nozzle outlet and the length of the cavity is much shorter than the actual length captured by the high-speed photography experiment. When the shear stress term is added to the critical pressure, the length of the predicted cavity is close to the experimental result and three phenomena of the jet are captured, namely, growth, shedding, and collapsing, which agrees well with the experimental high-speed image. According to the orthogonal analysis based on the simulation result, when the jet power is unchanged, the main geometry parameter of the divergent-convergent nozzle that affects the jet performance is the divergent angle. For the nozzle with three different divergent angles of 40°, 60°, and 80°, the one with the medium angle generates the most intensive cavitation cloud, while the small one shows the weakest cavitation performance. The obtained simulation result is confirmed by cavitation erosion tests of the Al1060 plate using these three nozzles.

1. Introduction

Cavitation is a well-known phenomenon in the field of fluid machinery. It can damage the hydraulic parts of the machines accompanied by noise and vibration [1,2,3]. Efforts have been made by researchers to avoid or reduce cavitation impact on pumps and turbines [4,5,6]. When the cavitation phenomenon is controlled and used properly, it will have enormous advantages in the engineering fields [7,8,9,10,11,12]. One of the commonly used techniques to generate intensive cavitation is a high-pressure submerged jet.

Research about the characteristics of high-pressure submerged cavitation jets has been mainly conducted by experimental methods. The effect of nozzle geometry, as well as jet parameters such as pressure and stand-off distance to the target, is investigated, aiming to improve the cavitation impact and promote the efficiency of the cleaning and cutting devices [13,14,15]. For the usage of cavitation peening, the residual stress, surface roughness and the micro-hardness of the peened metals are used to evaluate the performance of the jet. Since the jet velocity usually reaches more than 200 m/s, it is difficult to be detected by sensors. In addition, the flow region of the jet is filled with cavitation bubbles, which makes the Particle Image Velocimetry (PIV) difficult in process, where the bubbles may refract the laser light. Consequently, few experiments to detect the velocity field of high-pressure cavitation jets have been reported.

With the development of the CFD technique, the numerical simulation is becoming a common way to reveal the flow characteristic, which has advantages over experiments for investigating the theoretic phenomena in fluid dynamics. CFD codes have been widely used for predicting the performance of hydraulic machines [16,17]. Cavitating flow around hydrofoil or pumps is also simulated using multiphase flow models coupled with cavitation models based on the Rayleigh-Plesset equation. Qiang Guo et al. investigated the tip leakage vortex cavitation of a hydrofoil by coupling a cavitation model with a modified RANS model. It is found that the turbulence structure and the strength of the vortex have a great impact on the cavitating flow [18]. F. Khatami et al. [19] found that the compressibility of the gas-liquid mixture should be considered during the simulation of transient cavitating flow, and a set of URANS equations are used for the prediction of vortex cavitation in an elliptic hydrofoil. When the transient characteristics like vortex and shedding phenomenon are concerned, the turbulence model affects the simulation result to a large degree. The main reason is the compressibility and the large density gradient in the region that is close to cavitation clouds, and an accurate prediction on eddy viscosity is important for the simulation. Large Eddy Simulation (LES) has the ability to provide good results for cavitating simulation [20,21]. However, it is time-consuming and has a seriously high requirement for mesh quality. Consequently, different ways are developed to solve the transient cavitating flow using RANS models mainly by adjusting eddy viscosity which is usually over-predicted by RANS models [22,23,24].

The commonly used cavitation models that are available in most of the commercial algorithms are the Singhal model [25], Sauer model [26], and Zwart-Gerber-Belamri (ZGB) model [27]. The Singhal model is the so-called full cavitation model, which takes different parameters of the flow field and materials into consideration: mass transfer between phases and turbulent fluctuation of pressure effect of incondensable gases. The ZGB model assumes that the bubbles are mono-dispersed and the mass transfer rate during the cavitation can be acquired from the bubble density numbers and the mass change rate of a single bubble that is deduced from the Rayleigh-Plesset equation. The Sauer model was created in the same way as the Singhal and ZGB models, where the expression for the net mass transfer between liquid and vapor is derived. Guoyi P. et al. [28] conducted simulation on a submerged jet using a compressible multiphase flow model coupled with a cavitation model deduced from the Rayleigh-Plesset equation and the periodical shedding process was predicted and compared with experimental results. The period for the growth and shedding of the simulation result was close to the experimental data, while the cavity length of the CFD result was a little shorter than the high-speed photography images. This indicates that the mass transfer rate of the cavitation is underestimated. The cavitation models are also used for predicting the diesel fuel atomization characteristics to optimize the orifice geometry, which reduces the time required for nozzle design [29,30,31,32]. To take the effect of heat transfer on the cavitation process into consideration, cavitation models that consider the thermal effect are built while simulating the cryogenic cavitation flow in liquid natural gas or liquid nitrogen. Sun et al. [33] investigated unsteady sheet cavitation over a hydrofoil in a thermo-sensitive fluid using a modified ZGB model that contains the thermal effect. The cavitation models mentioned above assume that cavitation happens when the local pressure drops below a critical threshold value and the saturation pressure of the liquid is used as the threshold pressure for the cavitation inception. According to the thermal dynamic theory, the saturation pressure is deduced under the condition of a steady and equilibrium state for fluid where the shear stress is ignored during the movement of the liquid. For low-speed flows, the shear stress is relatively low, and the cavitation formation can be predicted using models without considering it. For cavitation jets, the exit velocity at the nozzle outlet is more than 200 m/s while the fluid of the ambient environment is almost steady, and the local shear stress is large.

In the current research, a high-pressure submerged cavitation jet is simulated using a cavitation model considering the effect of the shear stress and RNG turbulent model with FBDCM adjustment for turbulent viscosity. The transient characteristic of the cavitation cloud is captured and validated using a high-speed photography experiment. Eventually, the erosion experiment on Al 6061 material is used for testing the optimization of the nozzle geometry.

2. Method for Experiment and Simulation

2.1. Experimental Method

Figure 1 shows the cavitation jet experiment platform and test system. The system uses an Italian AR high-pressure piston pump to provide high pressure for the jet. The maximum working pressure of the piston pump is 50 MPa, the rated rotating speed is 1450 rpm, and the flow rate is 15 L/min. The water used in the experiment is pure water and the temperature is 25 °C. The impurities are removed by a Y-type filter before being delivered to the plunger pump. The pressure relief valve and the pressure gauge are connected downstream. The upstream pressure of the nozzle is controlled by adjusting the speed of the plunger pump. In order to prevent the influence of the upstream pipe bending on the nozzle pressure, a 300 mm stainless steel straight pipe is connected upstream to the nozzle. The jet platform is divided into two parts, a test water tank, and a water storage tank. The water storage tank is used to supply the water to the plunger pump and is located at the bottom of the test water tank. In order to facilitate the visual study of the submerged jet, the test tank is made of a transparent material polymethyl methacrylate, which has a refractive index close to that of water. It can effectively avoid positional errors caused by photography. In order to ensure that the outlet can return to the same position after each nozzle change, the nozzle is fixed on the three-degree-of-freedom moving slide by the clamp, and the repeated positioning accuracy of the slide is 0.01–0.02 mm.

2.2. Numerical Method

2.2.1. Multiphase Model

In this study, mixture model is used to obtain the velocity and pressure field of each phase. In the simulation, the flow is assumed to be isothermal, incompressible and the fluids are Newtonian. The conservation equation for the mass and momentum of the mixture is as follows:

$$\frac{\partial }{\partial t}\left({\rho }_m\right)+\nabla ·\left({\rho }_m{\overrightarrow{v}}_m\right)=0$$

(1)

$$\begin{array}{c}\frac{\partial }{\partial t}\left({\rho }_m{\overrightarrow{v}}_m\right)+\nabla ·\left({\rho }_m{\overrightarrow{v}}_m{\overrightarrow{v}}_m\right)=-\nabla p+\nabla ·\left[{\mu }_m\left(\nabla {\overrightarrow{v}}_m+{\overrightarrow{v}}_m^T\right)\right]\\ +{\rho }_m\overrightarrow{g}+\overrightarrow{F}+\nabla ·({\displaystyle {\displaystyle \sum }_{k=1}^n}{\alpha }_k{\rho }_k{\overrightarrow{v}}_{dr,k}{\overrightarrow{v}}_{dr,k})\end{array}$$

(2)

where ${\overrightarrow{v}}_m$ is the mass-averaged velocity and ${\rho }_m$ is the mixture density:

$${\rho }_m={\displaystyle \sum }_{k=1}^n{\alpha }_k{\rho }_k$$

(3)

${\mu }_m$ is the viscosity of the mixture, which is defined as:

$${\mu }_m={\displaystyle \sum }_{k=1}^n{\alpha }_k{\mu }_k$$

(4)

n is the number of phases, $\overrightarrow{F}$ is a body force and ${\overrightarrow{v}}_{dr,k}$ is the drift velocity for the secondary phase k.

2.2.2. Turbulence Model

Since the velocity of the high-pressure jet is extremely high around 200 m/s at the nozzle outlet, the Reynolds number of the main flow is high, and the turbulent model is necessary for the simulation. In the current research, the RNG k-ε turbulence model is employed. RANS models (e.g., k-ε and k-ω models) are widely used in cavitating two-phase flow, which tends to over-predict the turbulent viscosity in the cavitating region. For this reason, turbulence viscosity is modified in order to capture the shedding phenomenon of the cavitation jet in the current research. The transport equations of the turbulent model proposed by YAKHOT [34] are expressed as follows:

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}[{\alpha }_k\left(\mu +{\mu }_t\right)\frac{\partial k}{\partial x_j}]+G_k+G_b-\rho \epsilon -Y_M+S_k$$

(5)

$$\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\frac{\partial }{\partial x_i}\left(\rho \epsilon u_i\right)=\frac{\partial }{\partial x_j}[{\alpha }_{\epsilon }\left(\mu +{\mu }_t\right)\frac{\partial \epsilon }{\partial x_j}]+C_{1\epsilon }\frac{\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}-R_{\epsilon }+S_{\epsilon }$$

(6)

$${\mu }_t=\frac{C_{\mu }\rho k^2}{\epsilon }$$

(7)

where k and ε are the turbulent energy and turbulent dissipation rate, G_k is turbulent energy generation term, the default values for the empirical constants are: $C_{1\epsilon }=1.42,C_{2\epsilon }=1.68$. C_μ is coefficient of the model and the default value is 0.09. Since the shedding phenomenon has an important effect on the development of cavitating flows, a reasonable turbulence production is necessary for the simulation. Two main modification methods are used in cavitation flow simulation based on RANS models, namely, the filter-based model (FBM) [35] and density corrected model (DCM) [36]. The main idea of the DCM is to modify the eddy viscosity in the region filled with cavitation clouds, where high-density gradient exists. For the FBM approach, turbulence viscosity is modified where the turbulence length scale is larger than the local mesh size. Recently, a new method named the filter-based density correction model (FBDCM) was proposed, which combines DCM and FBM approaches using a bridging function [37]. The advantages of FBDCM over FBM and DCM are confirmed by Huang et al. [38] and Yu An et al. [39]. The adjusted turbulence viscosity is as follows:

$${\mu }_{t\_FBDCM}=\frac{C_{\mu }{\rho }_m{k}^2}{\epsilon }{f}_{FBDCM}$$

(8)

$$f_{FBDCM}=\psi \left({\rho }_m/{\rho }_l\right)f_{FBM}+\left[1-\psi \left({\rho }_m/{\rho }_l\right)\right]f_{DCM}$$

(9)

$$f_{FBM}=\min \left(1,\frac{\lambda ·\epsilon }{k^{3/2}}\right)$$

(10)

$$f_{DCM}=\frac{{\rho }_v+{\left(1-{\alpha }_v\right)}^n\left({\rho }_l-{\rho }_v\right)}{{\rho }_v+\left(1-{\alpha }_v\right)\left({\rho }_l-{\rho }_v\right)}$$

(11)

$$\psi \left({\rho }_m/{\rho }_l\right)=0.5+\tan \mathrm{h}\left[\frac{C_1\left(C_3{\rho }_m/{\rho }_l-C_2\right)}{C_4\left(1-2C_2\right)+C_2}]/[2\tan \mathrm{h}\left(C_1\right)\right]$$

(12)

Here, the constants are set as suggested, C₁ = 4, C₂ = 0.2, C₃ = 0.6, C₄ = 0.2, and n = 10 according to the literature. The filter scale λ is set as 1.1 times the mesh size in the refined region.

2.2.3. Cavitation Model

When the mixture model is used to describe the multiphase flow with cavitation, the transport equation of the vapor mass fraction with mass transfer is in the following form:

$$\frac{\partial }{\partial t}\left(\alpha {\rho }_v\right)+\nabla ·\left(\alpha {\rho }_v{\overrightarrow{v}}_v\right)=R_e-R_c$$

(13)

where R_e and R_c are evaporation and condensation rates. To provide the value for the source term, different models are deduced based on Rayleigh-Plesset equation. In real cavitating flow, bubbles are poly-dispersed and evaluate with coalescence and breakage besides growth and collapse. If all the factors are considered, the model is too complicated with a highly nonlinear feature, which increases the CPU costs and makes the simulation difficult to converge. One of the reduced models is proposed by Zwart-Gerber-Belamri [27], where the bubbles are assumed to be mono-dispersed and the mass transfer rate is calculated using bubble density numbers. The mass transfer term of the ZGB model is as follows:

If $\mathrm{P}\le P_v$

$$R_e=F_{vap}\frac{3{\alpha }_{nuc}\left(1-{\alpha }_v\right){\rho }_v}{R_B}\sqrt{\frac{2}{3}\frac{P_v-P}{{\rho }_l}}$$

(14)

If P $>P_v$

$$R_{\mathrm{c}}=F_{cond}\frac{3{\alpha }_v{\rho }_v}{R_B}\sqrt{\frac{2}{3}\frac{P-P_v}{{\rho }_l}}$$

(15)

where R_B is the bubble radius, α_nuc is the nucleation site volume fraction, $F_{vap}$ is the evaporation coefficient and $F_{cond}$ is the condensation coefficient. The default values for the constants are R_B = 10⁻⁶ m, α_nuc = $5\times {10}^{-4}$, $F_{vap}=50$, $F_{cond}=0.01$.The influence of turbulence on the threshold pressure is modeled in the way proposed in Singhal et al. Model [25]:

$$P_v=P_{sat}+\frac{1}{2}\left(coeff\right){\rho }_l{k}_l$$

(16)

where ${\rho }_l$ and $k_l$ are the liquid phase density and turbulence kinetic energy, respectively. The recommended value for coeff is 0.39 and this value is used by default.

The Zwart-Gerber-Belamri model is widely used to predict the cavitating flow fields in hydraulic machines like pumps and hydrofoils. However, the threshold for the pressure that drives evaporation and condensation is static pressure, which suits the cavitation flow without or only with slight shear stress. In high-pressure submerged cavitation jet flow, the velocity at the nozzle exit reaches more than 200 m/s, and the velocity gradient between the jet and the ambient water is extremely high, which generates high local shear stress. The shear stress layer is proposed to have a great impact on the cavitation generation and development of the jet [40]. In the current research, the threshold pressure is adjusted, and the shear stress is considered.

$${\tau }_{ij}=-p{\delta }_{ij}+S_{ij}=\left[\begin{array}{ccc}-p+S_{11}& S_{12}& S_{13}\\ S_{21}& -p+S_{22}& S_{23}\\ S_{31}& S_{32}& -p+S_{33}\end{array}\right]$$

(17)

The critical pressure threshold for the cavitation inception should be the maximum eigenvalue of the stress tensor in Equation (17). To reduce the difficulty to calculate the eigenvalue, the shear strain rate is used to model the maximum viscous stress tensor:

$$\dot{\gamma }=\sqrt{2D_{ij}{D}_{ij}}$$

(18)

where D_ij is the deformation rate, defined as $D_{ij}=\left(\partial u_i/\partial x_j+\partial u_j/\partial x_i\right)$.

The magnitude of the stress tensor can be expressed by Equation (19):

$$\left|{\tau }_{ij}\right|=\left|-p{\delta }_{ij}+\mu \dot{\gamma }{\delta }_{ij}\right|$$

(19)

The pressure threshold for the onset of cavitation in the flow field can be expressed as:

$$p-\mu \dot{\gamma }\le p_{saturation},p_{threshold}=\mu \dot{\gamma }+p_{saturation}$$

(20)

2.2.4. Geometry and Numerical Scheme

In this research, water jet from a submerged convergent-divergent nozzle is investigated. Figure 2 and Figure 3 illustrate the schematic configuration and the real photo of the simulated and tested nozzles. The key parameters include the throat diameter d, the convergent angle β, the divergent angle α, and the corresponding length of the three stages.

Figure 4 shows the simulation domain for the jet flow field. A structure mesh is created using the commercial software ANSYS ICEM 17.1. The inlet of the nozzle is set as a pressure inlet boundary, and the static pressure is inserted as 20 MPa, which is consistent with the experimental condition. The free surface between liquid and air is simplified to an asymmetrical boundary using rigid-lid hypothesis, where the normal velocity and the gradient of the variables are set as zero. This hypothesis has been verified and is used in research about the open Chanel flows [41]. The pressure outlet boundary condition is used for the out boundary of the domain and the pressure is set as 1 bar. The SIMPLE algorithm is used for the velocity pressure coupling, and the time step size is set as 5 μs. The iteration is conducted until the flow field is fully developed.

3. Results and Discussion

3.1. Validation of the Simulation Method

To describe the flow characteristic inside the nozzles, global parameters such as the discharge coefficient and local parameters such as vapor distribution and velocity field can be used. In the current research, the flow discharge coefficient is used for the validation of the numerical accuracy, which is defined as:

$$C_d=\frac{Q_{eff}}{Q_{ideal}}$$

(21)

The ideal flow rate through the nozzle, which is deduced from the Bernoulli equation, is as follows:

$$Q_{ideal}=A\sqrt{\frac{2\left(p_1-p_2\right)}{\rho }}$$

(22)

where A is the area of the throat section, p₁ and p₂ are the pressure upstream and downstream of the nozzle.

The effective flow rate predicted by the simulation is calculated using the following integral equation:

$$Q_{eff}=\underset{0}{\overset{R}{{\displaystyle \int }}}2\pi ru\left(r,z\right)dr$$

(23)

where r is the radial position and u is the local axial velocity at exit of the nozzle throat.

Figure 5 shows the discharge coefficient by experimental testing and numerical prediction. It can be found that the simulation result has a good agreement with the tested value, especially under a low-pressure condition. The difference between the CFD and experimental value increases slightly as the pressure rises, while the largest deviation is around 5%. The discharge coefficient has a slightly increasing tendency with the increase of the upstream pressure of the nozzle. The difference between the simulation and the experiment result is larger under the low-pressure condition because the relative test error is smaller when the flow rate is larger. In general, it can be concluded that the current simulation method provides a credible result.

Figure 6 shows the vapor volume fraction distribution predicted by the ZGB cavitation model and the modified cavitation model which is compared with the experimental result. In the ZGB model, the generation of cavitation is triggered by the pressure drop at the nozzle throat, where the pressure must be below the saturated vapor pressure. Thus, the bubbles collapse quickly when they drift to the high-pressure regions. As such, only a short cavity appears right at the nozzle exit, which is quite different from the experimental result. Comparing the cavitation bubble region predicted by ZGB cavitation model with the experimental image, it is obvious that this model underestimates the cavitation rate of the jet and the shedding while the collapsing process of the cavitation cloud cannot be captured. When the effect of shear stress is considered, the simulated cavitation cloud developed in a similar way to that shown in the experiment. In the central region of the jet at the nozzle outlet, the vapor volume fraction is almost zero. This is caused by the sudden change of dynamic pressure to static pressure for the high momentum jet stream. At the same time, a ring-type cavitation zone surrounds this central region, which develops from the high-speed shear layer. Strong and concentrated vortexes are formed under the effect of the intensive shear stress, which is verified via experiments carried out by researchers [42,43]. When the cavity grows and reaches a special length, the downstream part sheds off from the main cloud. From the UDC result shown in Figure 6, three different zones, namely, growing, shedding, and collapsing zone can be seen clearly. This phenomenon agrees well with experimental images taken by high-speed photography. The shedding process is proposed to have a great effect on the performance of cavitating jet used for cavitation peening, cleaning, and rock breaking since the collapsing of the cloud happens mainly after it sheds off from the cloud [44]. After the collapse, the cavitation bubbles with small diameter diffuse under the effect of turbulent flow. The cavitating area increases while the vapor volume fraction becomes lower. Comparing the simulation result of the new cavitation model with the experimental image, it can be judged that the current modified cavitation model predicts high-pressure cavitation jets with a more reasonable result.

To reveal the transient characteristic of the cavitation jet, the evolution process of the jet cavity is captured by numerical simulation and experiment. Figure 7 and Figure 8 show the evolution of the cavitation cloud of cavitation jet from a convergent-divergent nozzle with 60° divergent angle at 20 Mpa. One period for the development that includes growing, shedding and collapsing is also illustrated. From Figure 6 it can be found that the time between the two shedding is 1250 μs, which indicates that the period for the development of the cavitation jet is around 1250 μs. After the first shedding, the upstream part starts to collapse and diffuse, while the downstream part grows continuously. From Figure 8, it can be seen that the period for the evolution of the cavitation jet is in agreement with the experimental images. As an advantage of numerical simulation, more details about the velocity profile, vapor distribution, and information about the pressure of the jet can be obtained from the simulation. At the moment of t = t₀, the part ahead of the cavitation cloud has just shed off and starts to collapse, while the cavitation cloud at the nozzle outlet is thin and short. At t = t₀ + 250 μs, the part near the nozzle outlet extends in both axial and radial direction. The profile of this cavity is kept regular and the boundary is smooth. At the same time, under the effect of the turbulent disturbance and large dimensional vortex, the shed cavity rolls up to a ring shape, which disappears at t= t₀ + 750 μs. The cavity reaches the maximum length at t = t₀ + 1000 μs, after which the second shedding happens. Before shedding, the profile of the cavitating cloud is changed. At the nozzle outlet, the diameter of the cavity is reduced again, and the boundary becomes curved under the effect of the surrounding flows. By comparing the developing process of cavitation jet illustrated in Figure 7 and Figure 8, it can be found that the currently developed cavitation model can also predict the transient characteristics of high pressure submerged cavitation jets reasonably.

3.2. Influence of Geometric Parameters on Nozzle Performance

3.2.1. Nozzle Geometry Analysis Using Design of Experiments

The validated simulation method is used for the nozzle optimization in combination with Design of Experiments (DoE) using an orthogonal design, which reduces the calculation and experiment efforts in a large degree. According to former research, the cavitation performance of convergent-divergent nozzle is affected by four main parameters which are shown

Table 1. Factor level.

Number	A L2 (mm)	B d (mm)	C L3 (mm)	D β (°)
1	2	0.5	2	40
2	4	1.0	4	60
3	6	1.5	6	80

, namely, throat length L₂, throat diameter d, divergent length L₃ and divergent angle β. To evaluate the effect of these parameters on the cavitation jet, orthogonal analysis is carried out on the nozzle based on the result of numerical simulation. Three different levels are proposed in each factor, as shown in Table 1. The simulation plane is designed according to the L₉ orthogonal table, which is shown in

Table 2. Orthogonal test schemes.

Number	A	B	C	D
1	2	0.5	2	40
2	2	1.0	4	60
3	2	1.5	6	80
4	4	0.5	4	80
5	4	1.0	6	40
6	4	1.5	2	60
7	6	0.5	6	60
8	6	1.0	2	80
9	6	1.5	4	40

.

The most effective way to evaluate the performance of the cavitation nozzle is to compare the erosion capacity of the jet on the impinged material surface. Most of the cavitation models available do not contain such functions because the shock wave impact is difficult to be captured by CFD simulation. Also, the fluid-structure interaction related to the cavitation erosion process is difficult and consumes an enormous effort of CPU. Here, the mean value of vapor volume fraction on a radial distributed line at x = 0.1 m is used for evaluating the nozzle cavitation performance where the calculated result for the 9 cases is shown in

Table 3. Orthogonal calculation result.

Number	1	2	3	4	5	6	7	8	9
VOF	0.039	0.049	0.091	0.023	0.069	0.199	0.019	0.064	0.179

.

The main function of orthogonal analysis is to determine the factor that affects the objective function mostly with the minimum amount of data, which saves resource and time for experiments or calculation. Another advantage of orthogonal analysis is the ability to optimize the performance of the products that are affected by multiple factors. In the current research, this method is used to find the vital parameter that affects the cavitation of the high pressure submerged jet. The result of range analysis is listed in

Table 4. Range analysis of vapor volume fraction.

Number	A	B	C	D
K1	1.027146	0.264775	1.071998	1.169487
K2	1.058102	1.089607	0.96127	0.953527
K3	0.979006	1.709872	1.030986	0.94124
k1	0.342382	0.088258	0.357333	0.389829
k2	0.352701	0.363202	0.320423	0.317842
k3	0.326335	0.569957	0.343662	0.313747
R	0.026365	0.481699	0.036909	0.076082

. It can be found from the table that the range arrangement is R_B > R_D > R_C > R_A, which means that the main factor for the cavitation capacity of the nozzle is the throat diameter. In all calculated cases of orthogonal analysis, the upstream pressure is kept constant and as a result, more power should be consumed to supply the pressure as the throat diameter of the nozzle is increased. The second important parameter is the divergent angle of the nozzle, which affects the cavitation performance mostly when the supplied power is not changed. Consequently, the effect of the nozzle divergent angle on the cavitation jet is analyzed in detail.

3.2.2. Effect of Divergent Angle on the Nozzle Performance

Since the vapor volume fraction of the jet varies periodically, it makes no sense to compare the vapor distribution of different nozzles at an arbitrary moment. To make the volume fraction comparable, the time-averaged value for one period is calculated. Figure 9 shows the time-averaged volume fraction of cavitation jet generated by convergent-divergent nozzles with three different divergent angles, namely, 40°, 60°, and 80°. It can be seen from the figures that, in the time-averaged volume fraction contour there is a region with high vapor concentration which is located right at the nozzle outlet. The length of the high vapor volume fraction region varies with the increase of the nozzle divergent angle. The nozzle with medium divergent angle (60°) generates the longest cavitation cloud with high vapor concentration. In the nozzle with the smallest divergent angle (40°), the cavitation cloud is much thinner and shorter than that of the other two nozzles.

To analyze the cavitation growing and collapsing characteristics of jets with nozzles with different divergent angles, the vapor volume fraction along the radial distributed lines at different axial position, namely, x = 0.04 m, x = 0.06 m and x = 0.08 m is extracted and compared, which is shown in Figure 10. It can be seen that the highest values for the volume fraction of 60° and 80° nozzles are similar, which are close to 0.85 at x = 0.04 m, while the highest volume fraction of 40° nozzle is only around 0.7. This means that in the region close to the nozzle outlet, the cavitation generation rate of nozzle with a large divergent angle is faster than that with a small divergent angle. This can be explained by the discharge coefficient. The only difference between these three nozzles is the divergent angle. When the divergent angle is extremely small (closed to 0°) the throat length is increased, which increases the resistance of the fluids to flow past the nozzle and thus reduces the jet velocity when the upstream pressure is kept constant. The shear stress that promotes cavitation generation is reduced when the velocity gradient becomes smaller. At position x = 0.06 m, the volume fraction of 80° nozzle drops quickly to around 0.3 which is close to the 40° nozzle, while for the 60° nozzle it remains close to 0.5. This indicates that the flow field of a nozzle with a medium divergent angle can generate a cavitation cloud with larger size, namely, more cavitation bubbles.

To see the effect of divergent angles on the velocity field, velocity at four different positions of the nozzles is plotted in Figure 11a–c. According to the empirical theory, the velocity at the nozzle outlet follows the similarity law when the axial velocity and the radial position are normalized by dividing u_m and y_0.5, where u_m is the maximum velocity along the radial distributed line at special axial position, y_0.5 is the radial position where the axial velocity is half of the maximum velocity. From the first three figures, it can be found that the velocity at the center of the jet is the highest, which decreases gradually along the radial direction. The center velocity also decreases along the axial direction. The axial velocity decreases quickly at x = 0.02 m, and the decrease rate slows down along the direction of flow. At the position of x = 0.02 m, the maximum axial velocity increases from around 100 m/s to 120 m/s when the divergent angle of the nozzle is increased from 40° to 80°, because the increased angle provides less constraint on the jet and thus the resistance is reduced. Figure 11d–f show the normalized velocity distribution of jets from these three nozzles. It can be found that the jet with different nozzles follows the similarity law and the distribution curves of the velocity are almost coincident with each other.

To verify the optimization result, the nozzle with three different angles which is simulated above is manufactured and tested. The cavitation impact is evaluated using the erosion capacity of the Al1060 material. As the divergent angle affects the flow field and the velocity distribution, the optimal stand-off distance varies for different nozzles. Figure 12 shows the eroded aluminum surface by the three nozzles at the optimal stand-off distance. It can be found that the eroded area and the degree of the metal by the 40° nozzle are the slightest. The surface is smooth and only small shallow pits appear after 60 min impingement. The 60° nozzle creates the most severe damage on the aluminum surface where deep pits appear on a large area of the impinged surface and parts of the material are removed from the plate after the surface fatigue. For the 80° nozzle, the damage on the metal surface is slightly weaker than that of the 60° nozzle, the roughness of the surface is increased, and the fatigue starts to happen, while little material is removed.

Figure 13 shows the mass loss of the eroded Al1060 by the three different nozzles with different impingement times. The mass loss increases gradually with time, while the nozzle geometry affects the mass loss rate considerably. After impingement for 90 min using nozzles of the three different configurations, the mass loss are 9.6 mg, 900 mg and 144 mg, respectively. The 60° nozzle has the best cavitation impact performance and the 40° one is the worst. This confirms the result concluded from the numerical simulation, and it proves that the simulation method can be used to evaluate the cavitation performance of the nozzles.

4. Conclusions

In this study, a high-pressure submerged cavitation jet is numerically investigated based on a cavitation model considering the effect of shear stress and RNG turbulent model with FBDCM adjustment for turbulent viscosity. The proposed method for simulation of cavitation jet is then used for the geometry optimization for a convergent-divergent nozzle based on DoE using an orthogonal design. The main obtained results are as follows:

(1)

The shear stress plays an important role in the simulation of the high-pressure cavitation jet. In the result of the ZGB cavitation model, the cavitation cloud of the submerged jet only appears close to the nozzle outlet and the length of it is quite different from the experimental result. When considering the shear stress in the model, the simulated result is matched with the experiment data, which contains the growing, shedding and collapsing parts of the jet. The period of the jet evolution predicted by the current cavitation model also agrees well with the high-speed photograph experimental result.

(2)

According to the orthogonal analysis, the main geometric parameter of the divergent-convergent nozzle that affects the jet performance is the nozzle diameter, and the divergent angle is the second most important factor. When the pressure is constant, the jet power increases with the increase of the nozzle diameter. Thus, the divergent angle of the nozzle plays the most important role in the cavitation impact of the jet, while the jet power is not changed.

(3)

The velocity field of the jet from different nozzles follows the similarity law, and the normalized distribution curves of the velocity are coincident for each case. The divergent angle of the nozzle affects the magnitude of the central velocity of the jet. The maximum velocity at the position 0.02m downstream of the nozzle outlet decreased by about 20m/s when the nozzle divergent angle was reduced from 80° to 40°.

(4)

According to the experimental test, the area of the metal eroded by the 40° nozzle is small and the surface roughness is still low, the 60° nozzle creates severe damage on the aluminum surface and deep pits appear on a large area of the impinged surface; the 80° nozzle shows a medium cavitation impact performance between the other two nozzles. The result using numerical simulation agrees well with the experimental erosion test, which concludes that the convergent-divergent nozzle has optimal cavitation performance when the divergent angle is close to 60°.