Structure optimization of submerged water jet cavitating nozzle with a hybrid algorithm

ABSTRACT It is critical to obtain a suitable geometry for a cavitating nozzle of submerged water jet, since a suitable geometry determines its cavitation intensity and cavitation distribution. This paper puts forward a water jet cavitating nozzle optimization platform aiming to improve its axial maximum vapor volume fraction. Considering the influence of all eight geometric parameters of the nozzle itself, a hybrid optimization algorithm is developed to optimize the nozzle by combining with the CFD technique. After the optimization, the optimal geometry results to a 9.41% increment in the axial maximum vapor volume fraction at 50 m underwater depth, and improves the cavitation effect obviously especially in deeper water. From the results, it is found that the cylindrical section diameter, contraction angle and diffusion angle have a great influence on the axial maximum vapor volume fraction, while the inlet diameter have the least effect on the same. The optimization presented in this study can lay the foundation for further study on water jet cavitating nozzle.


Introduction
Cavitation is a common but harmful phenomenon in hydraulic machinery field. Since it can cause strong erosion to equipment, it should be avoided or controlled as much as possible. The occurrence of cavitation is due to the sudden pressure drop of the flowing liquid, and the local static pressure is lower than the saturated vapor pressure, generating cavitation bubbles. These bubbles then undergo a period of development, migration and collapse. When the cavitation collapses in a high-pressure zone near the wall, it will generate huge energy and damage the surface, eventually forming cavitation erosion. In recent years, the erosion effect caused by cavitation has been gradually utilized in surface treatment, material cutting (Chen, Nie, Wu, & Li, 2006), cleaning (Chen & Hu, 2017a, 2017b, rock breaking and other fields. The cavitating water jet technology is one of the important means to achieve cavitation utilization and enhance the ability of material crushing and cutting. Because the cavitation effect of submerged water jet is better than that of non-submerged water jet (Kang, Zhou, Yang, & Wang, 2013;Wang, 2009;Yanaida, Nakaya, Eda, & Nishida, 1985;Yang & Lei, 2014), this paper focuses on the submerged water jet cavitation study. Cavitating nozzles are the most important actuators for generating cavitating jets. Different internal structure of a nozzle CONTACT Yanzhen Chen 634497026@qq.com has different flow state, which could increase in the flow velocity, and decrease in pressure. Also, the vortex is formed in the shear layer by the entrainment between the high-speed cavitating jet and the surrounding low-speed fluid. Then, vortex cavitation occurs in low-pressure region of the vortex core, shedding as a periodical phenomenon. In order to improve the cavitation effect of the cavitating jet and meet ideal application requirements, it is important to optimize the nozzle structure reasonably. Some previous studies focused only on the internal flow field of the nozzle. Yang and Lei (2014) studied the characteristics of the internal flow field of the artificial submerged jet nozzle and found that the jet and the surrounding water can produce strong shearing action, which causes uneven speed and pulsation, forms a negative pressure zone, and provides conditions for cavitation. Yang, Tan, Zhang, Liu, and Jin (2010) analyzed the pressure characteristics of the internal flow field of an anglepipe and an organ-pipe cavitating nozzle. It is concluded that the cavitation effect of the angle-pipe nozzle is superior to the organ-pipe nozzle under the specific structural parameters. Yao et al. (2015) compared the axial distribution of vapor volume fraction among four kinds of nozzles with the same size, including a conical nozzle, an angle nozzle, a contraction nozzle and a contractionexpansion nozzle. Among them, the angle nozzle was found to have the best cavitation effect. Guan, Deng, Guo, and Hua (2012) used a two-dimensional model of the axisymmetric angle nozzle instead of three-dimensional model to study the internal flow field of nozzle. Considering the radial and axial distribution of vapor volume fraction and velocity as the criterion for evaluating the cavitation intensity, the influence of cylindrical section diameter and the inlet or outlet pressure on the cavitation was analyzed.
In addition, most of the studies focused on the impact erosion of cavitating nozzle, ignoring the optimization of the nozzle's structure. Kang et al. (2013) studied the flow characteristics of the central-body nozzle under nonsubmerged and submerged conditions by numerical simulation, and verified the cavitation erosion ability of the nozzle by erosion test. It is found that the erosion ability of the submerged jet is significantly better than that of the non-submerged jet within a certain stand-off distance. The cavitating jet effect is affected by the stand-off distance, which is closely related to the distribution of the bubbles. This shows that the cavitation effect plays a key role in jet erosion. Li et al. (2009) studied the erosion performance of a conical and an angle nozzle. By comparing the simulated cavitation results with the experimental erosion patterns, the axisymmetric annular distribution of the bubbles corresponds to the annular erosion region formed on the target surface. Lei, Deng, Guan, Chen, and Zhang (2015) studied the cavitation process and flow field characteristics of a submerged water jet with the vapor volume fraction as the evaluation index of cavitation intensity. This also indicates that the vapor volume fraction is reliable as an indicator. Liu, Gao, Jiao, and Tian (2016) studied the effects of different diameter nozzles on the submerged jet dynamic pressure under different environmental pressures and jet pressure conditions, but did not involve cavitation models. Zhang, Liang, and Gao (2000) used a two-dimensional model to study the flow field characteristics of an axisymmetric ultrahigh pressure water jet impacting rock under submerged state. Chahine, Kapahi, Choi, and Hsiao (2016) used a twodimensional model to study the flow field characteristics of surface cleaning by cavitation bubble dynamics and collapse, which can remove the dirty particle near material wall. Hutli, Nedeljkovic, Bonyár, and Légrády (2017), Hutli et al. (2016a), Hutli, Nedeljkovic, Radovic, andBonyár (2016b) and Kang, Liu, Zhang, and Li (2017) studied the influences of different cavitation numbers, jet velocity, jet pressure, nozzle diameter, stand-off distance, impact angle and impact time on the cavitation erosion characteristics of the submerged cavitating jet impacting target surface. Liu, Shao, Kang, and Gong (2015) and Peng, Li, and Tian (2017) also analyzed the mechanism of cavitation incubation, development, shedding and collapse by erosion tests. The distribution of cavitation bubbles affected the cavitation area and stand-off distance, while the intensity of cavitation determined the number of bubbles and erosion pits. These just show that the degree and pattern of erosion damage have a direct and close relationship with the distribution and quantity of cavitation bubbles in the flow field. Cavitation erosion can be predicted approximately based on vapor volume fraction and cavity profiles by 2D Fluent model (Li, Pourquie, & Terwisga, 2014). The cavitation effect affects greatly the results of erosion tests, and the vapor volume fraction can be used as an indicator for evaluating the cavitation effect of nozzle, which is difficult to get in advance due to complex experiments and sophisticated test facilities.
In order to improve the cavitation effect of nozzle, some optimization work has been done. Cao, Meng, Sun, and Wang (2017) presented a new cleaning nozzle in their study. The length of nozzle's cylinder section was reduced. The optimization results show that the new model improved the cavitation performance. Deng, Shen, Li, Feng, and Tong (2008) presented a numerical simulation for the angle nozzle using a SIMPLEC and κ-ε turbulence model. The results showed that the diffuser angle of the nozzle had great influence on the internal flow field and had the optimal value of 30 degree. Tan (2011) used a Fluent software to simulate the angle nozzle and the shrinkable tubing nozzle. Then the diffusion angle, the diffusion segment and cylindrical section length were studied in order to obtain the optimal structural parameters. Fu et al. (2013) studied the different structure of nozzle influencing jet characteristics with orthogonal experimental design method. The results showed that the exit section shape of nozzle had the biggest influence, while the nozzle exit section length had the least influence.
As studies mentioned above, it is confirmed that the vapor volume fraction can indicate both the cavitation intensity and its distribution. It is suitable to use the vapor volume fraction as a criterion for evaluating the cavitation performance of nozzle. But these studies did not pay attention to the cavitation intensity of a specific location for engineering applications. Therefore, this paper focuses on the cavitation effect on the center axis of the jet, which should be as strong as possible, and uses the axial maximum vapor volume fraction as the objective value to represent the cavitation effect of the nozzle. Still, some of studies only considered the internal flow field of the nozzle, and did not involve the flow characteristics outside the nozzle. Therefore, in this paper, the external flow field is considered to fully observe the flow field evolution process between water jet and ambient water after the nozzle exit. This makes the simulation be more consistent with the actual submerged flow field. In addition, Some of these studies hardly considered the influence of the structural factors of the cavitation nozzle itself on the flow field, which does exist. These various geometric parameters can change the internal structure of the nozzle and have a certain influence on the flow characteristics and the intensity and distribution of the cavitation. In order to obtain better cavitation performance, the optimization of the nozzle structure is analyzed. It is worth noting that, in the previous researches, the cavitation performance of the nozzle was compared between only several representative values, which were not comprehensive, to get a better cavitation effect. The optimization method was not used to find the optimal structural parameters in a global position. Therefore, in this paper, a nozzle optimization loop is built to study the effects of different structure papameters in order to achieve better cavitation erosion effect.
The optimization algorithm is a method that utilizes a mathematical approach to obtain an optimal solution in a linear or non-linear research space. It has been widely used to find an optimal solution in various fields. Many improved optimization algorithms have been put forward in order to improve the search capability (Ghorbani, Kasaeian, Toopshekan, Bahrami, & Maghami, 2017;Maryam & Behrouz, 2018;Tezdogan et al., 2018;Zhang, Zhang, Tezdogan, Xu, & Lai, 2017). For the gradient algorithm, it is highly dependent on the initial point, and it is also likely to fall into a local solution. Zhang, Tezdogan, Zhang, Xu, and Lai (2018) used an experimental design method to find an initial value of the gradient algorithm, achieving higher evaluation accuracy. Thus, a Latin hypercube design (LHD) (McKay, Beckman, & Conover, 1979), an experimental design method, is selected in this study to obtain an optimal initial value, and a MMFD algorithm (Vanderplaats, 1984), a gradient algorithm, is used to search the local part in detail in the water jet cavitating nozzle optimization space.
With the rapid development of the approximation technique, many methods can be used to calculate the data approximately besides Computational Fluid Dynamics (CFD) tools (Taormina, Chau, & Sivakumar, 2015;Wu & Chau, 2011;Yaseen, Sulaiman, Deo, & Chau, 2019). However, the CFD tools can obtain the results more accurately. Therefore, the aim of this study is to present a water jet cavitating nozzle optimization loop with the CFD method and the hybrid algorithm. The model presented by Yao et al. (2014) is selected as the validation model. The ANSYS Fluent software is used to simulate the performance for the water jet cavitating nozzle, and the CFD results are validated against data by Yao et al. (2014). Then, eight parameters are selected as design variables, and the axial maximum vapor volume fraction is used as objective function. The optimization loop is developed to optimize the nozzle.

Governing equations
The submerged water jet is a turbulent flow in the vaporliquid two-phase flow field (Carlomagno & Ianiro, 2014;Hutli et al., 2016aHutli et al., , 2016bLei et al., 2015;Lin, Zhang, Qin, & Dang, 2015;Liu et al., 2016). Multiphase model, turbulence model and cavitation model should be considered for the water jet cavitating nozzle numerical simulation. In the literature of Lin et al. (2015), it is pointed out that the volume of fluid (VOF) model separates the vapor phase and the liquid phase in cavitating jet and is often used to track the interface of two or more incompatible fluids. Since the object of this study is submerged cavitating jet, in which the water jet is the same fluid as the surrounding medium in the outer flow field, and the characteristics of the whole jet are mainly explored, the VOF model is not suitable. Tan (2011) pointed out that the mixture model and the Euler model are the most commonly used multiphase flow models for cavitating jets. Both two models are suitable for the case where the volume fraction of the dispersed phase exceeds 10%. The Euler model has a long calculation time and its solution is easy to diverge, while the mixture model is suitable for studying the wide distribution of dispersed phases and is often used to study the location and intensity of cavitation. In general, the mixture model is selected to deal with similar issues by Celik, Arikan Ozden, and Bal (2014), Chen, Li, Gong, Chen, and Lu (2019), Ji, Luo, Wu, Zhang, and Xu (2010), Lei et al. (2015), Li, Zheng, Hong, and Ni (2017), Yang and Lei (2014), Zhao, Wei, Zhang, Dong, and Jin (2009), and Zhou and Wang (2008), and hence, the mixture model is selected as the multiphase model in this study.
The mixture continuity equation can be expressed as The momentum equation for the mixture model can be written as where ρ m is the mixture density, p is the mixture pressure, u i and u j are the velocity in the i and j direction, respectively, τ ij denotes the viscous tensor, which can be expressed as where μ m and μ t are the molecular and turbulence viscosities obtained from the turbulence model, respectively.
where β v is vapor volume fraction, ρ l and ρ v represent the density of liquid and vapor, respectively. μ l and μ v are the liquid and vapor dynamic viscosities, respectively.

Turbulence model
According to standard k-ε model, Yakhot and Orszag (1986) (2012), and Zhou and Wang (2008), it is verified that the RNG k-ε model can get satisfactory results in the CFD simulation of cavitation flow. In addition, because two-equation turbulence models are the most widely used in engineering calculations (Akbarian et al., 2018;Han, Xiong, & Chen, 2009;Mou, He, Zhao, & Chau, 2017;Ramezanizadeh, Nazari, Ahmadi, & Chau, 2019), a comparison has been made with the same Schnerr-Sauer cavitation model between different turbulence models, which are Standard k-ε, RNG k-ε, Realizable k-ε, Standard k-ω and SST k-ω, as shown in Table 1. It can be seen that the result of the RNG k-ε model has the smallest error with the Yao's results and is very suitable for cavitation simulation. Therefore, the RNG k-ε model is used to simulate the performance of water jet cavitating nozzle. The turbulent kinetic energy k and its dissipation rate ε are obtained from the following transport equations: where and , G k is the generation of turbulent kinetic energy by the mean velocity gradients, C 1ε and C 2ε are the empirical constants. α k and α ε are the reciprocal of the effective turbulent Prandtl number of the turbulent kinetic energy k and the dissipation rate ε, respectively.

Cavitation model
The axial maximum vapor volume fraction of different cavitation models, Zwart-Gerber-Belamri (ZGB) and Schnerr-Sauer (S-S), is shown in Table 2. It can be seen that the result of S-S model has a smaller error with the Yao's results. Moreover, the S-S model has more advantages in the calculation of the pressure coefficient, and the calculation is more stable and easy to converge (Lei et al., 2015). Therefore, the Schnerr-Sauer cavitation model is used as a cavitation model in this study. The conservation equation for vapor volume fraction has a general form as where the source terms S e and S c refer to the evaporation and condensation of vapor bubbles, respectively, accounting for the mass transfer between the vapor and liquid phases in cavitation.
The source terms are derived from the Rayleigh-Plesset equation and are defined as where P is local far-field pressure, P v is saturation vapor pressure. The bubble radius R can be related to the vapor volume fraction β v and the bubble number density N b as follows: where N b is the only parameter to be provided as input for FLUENT with a default value of N b = 1.0e + 13, which is also the optimal value proven by the research of Liu et al. (2012) and Li, Kelecy, Egeljamaruszewski, and Vasquez (2008a). The bubble number density is related to bubble nuclei distribution in the fluid. It is observed from the research of Li et al. (2014) that a larger bubble number density than the default value did not have significant effects on the vapor volume, however, a smaller bubble number density predicted a larger vapor volume. And it is explained that this phenomenon could be caused by an overshoot in the source strength fluctuations, which is proportional to bubble number density. This overshoot seems to be an artifact of this cavitation model and is not associated with physics. Therefore, the default value of N b = 1.0e + 13 is used in this study.

Hybrid algorithm
A hybrid algorithm is developed in this paper to obtain the global optimal solution. In this algorithm, a Latin hypercube design (LHD), developed by McKay et al. (1979), is selected to obtain an optimal global position for a non-linear optimization space, and a Modified Method of Feasible Directions (MMFD) algorithm (Vanderplaats, 1984) is used to search the local part in detail.
Latin hypercube design (LHD) is an effective method to explore the process characteristics for a research space with fewer samples. In LHD method, the design space for each factor is divided uniformly (the same number of divisions, n, for all factors). These levels are randomly combined to specify n points defining the matrix design (each level of a factor is studied only once) . The space-filling capacity of LHD method is good Figure 1. Samples designed using the LHD method (Zhao & Cui, 2007). and it can fit non-linear response easily. Figure 1 shows the samples designed with LHD method for two factors with nine points. The X1 and X2 are the two dimensions in the research space.
Modified Method of Feasible Directions (MMFD) method is a gradient optimization method which has been widely used for different optimization applications. This method exploits the local area around the initial design point and rapidly obtains a local optimum design. It is well-suited for highly non-linear design spaces (Velden & Koch, 2010). Thus, MMFD method is used as an optimization method in this study. Velden and Koch (2010) pointed out that MMFD method can be used to optimize the non-linear problems when starting from a feasible design point. In addition, Zhang et al. (2018) assumed that a LHD is a good method to obtain an optimal initial point for a gradient optimization method. Therefore, this study uses a LHD algorithm to choose a feasible design initial point. Figure 2 shows the flow chart for water jet cavitating nozzle optimization combining LHD and MMFD methods. The basic steps can be summarized as follows: (1) LHD algorithm is used to design some samples in a given space.
(2) The CFD technique is applied to estimate the objective function for each sample. (3) Compare each result and output the best sample. (4) The best sample is selected as an initial design point for a MMFD method. (5) Finally, MMFD method is employed to optimize this mathematical space. (6) Output the optimal solution and optimal geometry.

Optimization process
The optimization is carried out on a personal computer (Inter(R) Core(TM) i7-3770 K CPU @ 3.50 GHz), and each CFD case is computed for approximately 12 CPU hours. The total runs for the optimization is 40 times  by hybrid algorithm and only MMFD algorithm. The ANSYS Fluent software package is applied to simulate the flow field for the water jet cavitating nozzle. The LHD and the MMFD methods are employed to optimize the geometry model for the water jet cavitating nozzle. The complete optimization flow chart is presented in Figure 3.
(2) Change design variables using a MMFD method.
(3) Change the geometry using the full parametric modeling algorithm. (4) Build computational domain and calculate the objective function using CFD technique. (5) Output objective value. (6) Output the optimal geometry if the MMFD method is convergent. Otherwise, return to Step 2.

Geometry model
Conical nozzle and angle nozzle have been widely used in engineering practice, considering the jet characteristics and the difficult degree of manufacturing. Compared with the angle nozzle, the structure of the conical nozzle has no diffusion section, but other structural parts are the same. Since the cavitation effect for the angle nozzle is better than conical nozzle (Lu et al., 2009), in order to obtain the optimal cavitation effect, the angle nozzle is used in this study. A 2D parameterized model is established to express variables of the model structure. Figure 4 shows the geometry of angle nozzle. The cylindrical section length is 1 mm, the contraction angle is 6.75 degrees, and the diffusion angle is 30 degrees. Guan et al. (2012), Li et al. (2014) and Zhang et al. (2000) used the CFD technique to simulate the performance of two-dimensional nozzles. These results show that the two-dimensional model could descript the flow field characteristics of the actual model when the velocity in the z-direction is ignored. In the literatures of Carlomagno and Ianiro (2014), Hutli et al. (2016b), Lei et al. (2015), Li and Shen (2008), and Peng, Tian, Li, and Alehossein (2018), it is explained that the cylindrical symmetric nozzle can produce an average jet axisymmetric velocity profile in a two-dimensional plane, and their simulation results are well verified by experimental results. Therefore, for the simulation of cylindrical symmetric nozzles, the two-dimensional rotating section of the actual model can be used instead of the actual model to improve the calculation efficiency. Since the cavitating nozzle is also a cylindrical symmetric nozzle, the two-dimensional model is selected as the research objective without considering the velocity in the z-direction. Figure 5 shows the two-dimension model used in this study. The length of the computational domain is set as 150 mm, and the width of this domain is 50 mm, which can adequately satisfy the flow field boundary of submerged free jet (Zhang et al., 2000), and eliminate the influence of the flow field boundary on the simulation results.

Grid-independency study
A grid-independency study is undertaken to estimate the impact of the current CFD model on the axial maximum vapor volume fraction. Table 3 clearly shows the mesh configurations of the current CFD model, as well as the volume fraction. As can be seen from the table, the volume fraction value obtained by fine mesh is closer to the Yao's result, hence the fine mesh is used in this study.

Grid quality
Compared with the structured mesh, the computational process of the unstructured quadrilateral mesh is more complicated. However, the unstructured quadrilateral mesh shows the strong adaptability and is more accurate to capture the flow field (Wang, 2004). Therefore, the unstructured quadrilateral mesh is used as mesh generation method. In order to get a higher resolution and clearer grid figure, a grid that is five times larger than the grid size used in current actual simulation calculation, as shown in Figure 6. Figure 7 shows the refined mesh after the gird size is magnified five times near the water jet cavitating nozzle. The actual overall elements number is equal to approximately 60 thousand, and the average element quality reaches 0.971 by the local refinement of 0.05 mm element size. The average skewness is 0.031, the average orthogonal quality is 0.996, the average Jacobian ratio is 1.043, and the average aspect ratio is 1.047. According to Wang and Hu (2015), these evaluation parameters can show that the grid used in this paper is of good quality for the application of CFD numerical simulation.

Simulation setup
For all simulations, The pressure inlet condition is selected at the nozzle's inlet (inlet1) and the environmental inlet (inlet2). The pressure outlet condition is applied at the outlet boundary. The top boundary and nozzle are set as wall, and the bottom boundary is set as axis. The basic pressure parameters of this simulation are listed in Table 4 (Yao et al., 2014). These simulations are carried out under different water depths with different pressure values at inlet2 and outlet boundaries. Since the pressure difference between the inlet2 and outlet is very small, it can provide a certain kinetic energy to ensure that the ambient fluid can flow out through the outlet. The turbulence intensity is set according to the formula I = 0.16Re −0.125 , and the hydraulic diameter, which is 6 mm at the nozzle's inlet, are used as the turbulent boundary conditions. The nozzle and wall are regarded as no-slip walls. The materials of two phases are selected as that the vaporization pressure of water at 300 K is P v = 3540 Pa, where the density and dynamic viscosity of the waterliquid are ρ l = 998.2 kg/m 3 and μ l = 1.003e − 03 Pa s, respectively, and the vapor density is ρ v = 0.5542 kg/m 3 and the vapor viscosity is μ v = 1.34e − 05 Pa s. For the numerical strategy, the monitor is set to track the vertex of the maximum vapor volume fraction on the axis. Several time steps as 0.0001, 0.00001, and 0.000001 s, were tested with t = 0.00001 s, which can give reasonable results with relatively short calculational times. Different iteration per time step will result in different stages of the periodic cavitation clouds at the end of each time step, affecting the development of cavitation cloud. A small number of iteration may result in non convergence at each time step. The more the number of iteration is selected, the more development of cavitation will be at each time step. But if the iteration number is too many, the calculation time will be too long. For a periodic solution, an appropriate number of iteration per time step can obtain a more similar solution, as shown in Table 5. Therefore, the iteration number of 20 is selected. And the SIMPLEC algorithm is used to solve the pressure and velocity equations. The pressure is discretized by using the PRESTO! Scheme, and the first order upwind scheme was used for the momentum equations, vapor volume fraction equation and the convection terms in the turbulence equations.

Results and discussions
The simulations are carried out at different underwater depths. Figures 8-15 show the velocity distribution and    Table 6 shows the axial maximum velocity and the axial maximum vapor volume fraction at different underwater depths, in which the theoretical maximum jet velocity u max = 44.7 √ P t (Xue, 2006), where P t is the total pump pressure, and error values σ = (u − u Yao /u Yao ), where u is the theoretical value u max for σ 1 , and the current simulation value for σ 2 , respectively.
As can be seen from Figures 8-15, there are some deviations in the local part of curves. The Yao's curve fluctuates greatly because the value is obtained at the time when cavitation clouds are shedding, and current   simulation value is taken at a stable periodic interval, so the curve is smoother. As shown in Figure 13, the bulge on Yao's curve is caused by new cavitation cloud just shedding off, causing the whole curve to stretch forward. In Figure 14, the depression is caused by 'necking'. In Figure 15, the bulge indicates a state in which the cavitation cloud is about to disappear after escaping, and is far away from the overall bubble at the outlet. These deviations are caused by the periodic shedding of the vortex in   the shear layer, which causes the cavitation cloud to have four stages of incubation -developing -shedding -collapse. The first stage is the incubation of cavitation. When the pressure in the fluid suddenly drops below the vapor pressure, cavitation occurs and initial bubbles appear. There are two low-pressure region in the jet flow field. One is that cavitation may occur first in the cylindrical section, according to Bernoulli's theory, the high-speed jet will reduce the pressure in the cylindrical section, and the other is the low-pressure region of the vortex core in the shear layer around the jet, which is also prone to cavitation, and promotes the expansion and development of the previous bubbles, thus forming a big cavitation cloud. The growth process of these bubbles is the second stage of development. Subsequently, the cavitation cloud is affected by the re-entrant jet, chock effect and the vortex shedding, then forming a periodic shedding phenomenon. This is the third stage of shedding. Finally, in the fourth stage of collapse, when bubbles move with the jet to the downstream environmental pressure region or the high-pressure region near the wall, the collapse of bubbles will occur, generating shock waves and micro jets, which could hit the material surface and damage it. The periodic variation of cavitation clouds lead to pressure changes with different nozzle open area at different stages, as a result of the cavitation appearance and its 'chock effect', causing errors in the results (Peng et al., 2017). The shedding frequency is influenced significantly by the number of iteration per time step. Moreover, the interaction between the cavitating jet and the surrounding fluid is also random, which causes the target value to fluctuate a little. The error between this simulation results and Yao's results about the axial maximum vapor volume fraction and the axial maximum velocity are both within 1% deviation, as shown in the Tables 6 and 7, which is within acceptable range (Pan, Zhang, & Zhou, 2012). Overall, the trend of the curve under each working condition is the same, which can truly and accurately represent the changed trend of variable and reflect the optimized effect in the following structural optimization. Therefore, the simulation method is feasible. It can be seen from Figures 8, 10, 12, and 14 that the axial velocity distribution conforms to the jet distribution law. The velocity increases in the contraction section and then reaches a maximum within the cylindrical section. Also, a core constant velocity zone is formed. After the constant velocity section, as the distance on the axis increases, the velocity gradually decreases due to the turbulent diffusion caused by the action between the jet and the surrounding fluid. When the nozzle pressure is constant, as the water depth increases, the environmental pressure also increases, so that the length of the core constant velocity section is significantly reduced. The maximum velocity on the axis is related to the nozzle inlet pressure, so it varies little with water depth.
Focusing on the distribution of the axial vapor volume fraction in Figures 9, 11, 13, and 15, the cavitation bubbles are produced from the exit of the cylindrical section on the axis, then the vapor volume fraction is significantly increased in the diffusion section. The peak point of the curve indicates that there is a significant maximum vapor volume fraction on the axis. As the water depth increases, the axial vapor volume fraction gradually decreases, and this trend becomes more and more obvious with the increasing underwater depth. In shallower places, such as 10 and 25 m, the axial vapor volume fraction and the length of the cavitation bubble region on the axis are larger and less affected by the depth change. In deeper places, such as 50 and 100 m, the axial vapor volume fraction and the length of cavitation bubble region on the axis are obviously reduced, and will continue to decrease significantly as the water depth increases. The pressure of the surrounding environment increases with the increasing underwater depth, which inhibits the formation of cavitation bubbles and reduces the cavitation effect.
It can be confirmed that when the nozzle's structure at a certain working condition is optimized, the working effect of the nozzle on other depths can also be improved. The water depth has little effect on the axial maximum velocity. In the shallow water, the axial maximum vapor volume fraction does not change significantly with the water depth. However, if the depth is too deep, the environmental pressure is too large, which is not conducive for cavitation to occur. Combined with the actual depth of underwater operations, this paper selects environmental conditions at 50 m underwater depth, which also has a minimum error in Table 7, as an optimized working condition.

Optimization
For different structures and varying local parameters, the recognized best values will also be biased, therefore, it is necessary to find the best combination of geometric parameters to obtain an optimal cavitation effect.

Object function
The damage caused by the submerged cavitating jet consists of two parts, one is the impact force formed by the continuous jet, and the other is the microjet produced by the cavitation bubbles collapse. Among them, cavitation erosion by microjet accounts for the dominant position of destruction. For the impact force F d of the continuous jet, it is mainly affected by the jet velocity u max and the volume flow rate q v , according to the empirical formula F d = ρq v u max (Xue, 2006). When the inlet1 pressure is constant at working conditions, the axial maximum velocity almost unchanged, and when the nozzle diameter changes a little, the flow rate change is also not obvious. Most importantly, the microjet formed by the collapse of the bubbles is highly destructive. The fatigue damage it causes is the core of cavitation erosion. The vapor volume fraction indicates the amount of bubbles generated by cavitation. It can reflect the cavitation intensity. Since the nascent, growth and collapse of cavitation bubbles are obviously affected by the environment, incomplete development and early collapse may occur. Therefore, the more intense cavitation occurs, the greater the number of bubbles produced. This can lead to a possible increase in cavitation erosion, which increases damage to the surface. For the application of cavitation technology, it is often required to form a suitable degree of damage at a specified part of the target. Therefore, in addition to considering the cavitation intensity, the cavitation effect should also pay attention to the location of cavitation. The different distribution of cavitation bubbles in the flow field determines the pattern of different shapes formed by the destruction. At present, the patterns of cavitation erosion generated by nozzles in many experiments are mostly ring-shaped. However, in many applications, higher working accuracy is required, which requires a more concentrated range of erosion pits, causing the destruction energy to gather together, increasing the damage strength and improving working efficiency. Therefore, this paper does not consider the maximum vapor volume fraction in the entire flow field, commonly generated by the nozzle, but focuses on the maximum volume fraction on the center axis of the jet. The larger the axial maximum volume fraction, the greater the potential for strong central erosion will be. Thus, the axial maximum volume fraction is selected as the objective function in this study.

Design variables
The cylindrical section length and the diameter of the nozzle influence the flow resistance and flow coefficient, s well as flow state of thin-walled and thick-walled holes. A large cylindrical section will reduce the velocity at the outlet boundary. The contraction angle will influence the flow resistance and jet velocity. The diffusion section will improve the cavitation effect. However, a small diffusion section will reduce the number of bubbles, while a large diffusion section will increase the resistance and decrease the velocity at the outlet boundary. To obtain the optimal structure, eight parameters are selected in the optimization, as shown in Figure 4 where a is the inlet section length, b is the contraction section length, c is the cylindrical section length, d is the diffusion segment length, e is the inlet diameter, f is the cylindrical section diameter, g is the contraction angle, h is the diffusion angle. The value ranges of these parameters are listed in Table 8. Table 9 shows the 100 samples with different design parameters and corresponding axial maximum vapor volume fraction. Figure 16 shows the three-dimensional space distributions of 100 samples. The optimal sample nozzle is No.56 (a = 2.222 mm, b = 4.131 mm, c = 4.111 mm, d = 6.101 mm, e = 1.364 mm, f = 0.222 mm, g = 9.94°, h = 13.64°) with the corresponding axial maximum vapor volume fraction 0.96143.

Results and discussions
Next, an optimal sample (No.56) is selected as the initial value of the MMFD method in hybrid algorithm. Figure 17 shows the time history of axis maximum vapor volume fraction using hybrid algorithm. Figure 18 shows the time history obtained by using only MMFD method, and the initial value is the axis maximum vapor volume fraction for the original model. After the optimization, the axis maximum vapor volume fraction for the optimal model increases by 9.41% and by 3.63% using hybrid method and only MMFD method respectively. As can be seen from the results, the hybrid algorithm can obtain a better result than only MMFD method. Therefore, the hybrid algorithm presented in this study is a promising method for the nozzle optimization. Table 10 shows the optimal solution obtained using the hybrid method. As shown in Table 10, the inlet section length, the cylindrical section length, the diffusion segment length and the contraction angle increase. The Table 9. Design parameters with corresponding axial maximum vapor volume fraction. contraction section length, the inlet diameter, the cylindrical section diameter and the diffusion angle decrease. Figure 19 shows a comparison of the geometry for the original model and the optimal model. Figure 20 shows the comparison of velocity distribution contours between original and optimal models. It can be seen that, due to the Bernoulli effect, both of them form a core high-speed region from the entrance of the cylindrical section in the red region, and the maximum speeds are not much different. In the radial direction, the speed at the axis is the largest and gradually decreases toward the periphery. In the axial direction, a conical core velocity zone is formed from the cylindrical section, wrapped in the peripheral diffusion zone. The diffusion zone interacts with the surrounding low-velocity fluid to form the turbulent boundary layer, where energy is dissipated and the velocity is reduced. The core velocity region of the original nozzle is relatively concentrated, narrow and longer, so the free space in the diffusion section is large. However, the core velocity region of the optimized nozzle has a larger coverage area in the diffusion section, resulting in a relatively faster diffusion. Figure 21 shows the comparison of the axial distribution of velocity between original and optimal models. It can be found that the axial maximum velocity is equal to the maximum velocity of the entire flow field, indicating  that the position where the flow velocity is largest appears on the axis, and a significant difference in the trends can be observed. For the axial velocity distribution of the original model, the velocity reaches a peak in the cylindrical section and forms the core constant velocity zone. Then the velocity drops slightly outside the nozzle outlet. Subsequently, there is a significant decrease as the  axial distance increases. However, the axial velocity distribution of the optimized model is very different. There is a distinct two-order constant velocity section. In the cylindrical section, the increase of speed is significantly lower than before, and the maximum speed of the axis is reached at the exit of the cylindrical section. The first segment of the constant velocity region is then formed within the diffusion segment. The first significant velocity drop occurs at the exit of the diffuser section, followed by the second constant velocity zone. The total length of the two constant velocity segments is similar to that of the original model. Subsequently, as the x-axial distance increases, a significant speed reduction occurs. The impact kinetic energy of the optimized nozzle is reduced faster, but the fatigue action of the micro jet formed by the bubble collapse is the main destructive factor. Figure 22 shows the comparison of the radial velocity profiles at different axial positions between original and optimal models. Four positions at the axis are chosen as x = 20 mm, x = 22 mm, x = 25 mm and x = 30 mm. It can be seen that for both the original and optimal models, the velocity has a similar trend. And the velocity of the   optimal model is significantly lower than the velocity of the original model. As the x-axial distance increases, the velocity at the axis, where the radial distance is 0, gradually decays and the decrease becomes larger and larger, which is consistent with the above analysis of the axial velocity. There is very little speed reduction in the core constant velocity zone. After leaving the core area, the speed is significantly reduced by the existence of shear layer, such as the curve of x = 30 mm, which is obviously lower than the curves of other positions closer to the outlet. At different axial positions, the radial velocity also has the same trend. The closer to the central axis, the higher the speed. As the radial distance increases, a constant velocity section in the core zone is also generated, and the speed reduction is small. Subsequently, the jet enters the boundary layer with a faster speed drop. As the axial distance increases, the radial constant velocity section becomes shorter and shorter. This also corresponds to the above description of the conical core velocity zone. Figure 23 shows the comparison of vapor volume fraction distribution contours between original and optimal models. It can be seen that both nozzles form a conical cavitation zone. The bubble region of the original model is concentrated near the axis, forming three layers of weak-strong-weak, while the cavitation region of the optimized model covers almost the entire diffusion segment, the red area is large and the zone is wide. The conical cavitation region converges to the center of the axis outside the nozzle. The periphery also produces a transition layer with a lower volume fraction. The maximum values of the entire flow field of the two models are 0.966 and 0.984 near the wall of the diffusion section, respectively, which are significantly improved. Figure 24 shows the comparison of the axial distribution of vapor volume fraction between original and optimal models. It can be seen that the axial vapor volume fraction of both models are generated from the inlet of the diffusion section, and gradually disappears near x = 35 mm position. The lengths of the cavitation zone on the axis are substantially equal. For the original model, the axial maximum axis vapor volume fraction is reached outside the nozzle, then as the distance along the axis increases, the value begins to show a significant drop. In contrast, the optimized model reaches the axial maximum vapor volume fraction near the nozzle exit, which is 0.965, significantly larger than the value of 0.882 in the original model, indicating that the cavitation intensity is improved. Subsequently, the volume fraction drops very slowly, and within a long axial distance in the outer flow field, the axial vapor volume fraction remains at higher values, significantly better than the value of the original model that has been significantly reduced. When the position at the axis is near x = 30 mm, it begins to drop rapidly until it disappears. The vapor volume fraction of the optimized axis distribution above 0.8 has a longer length. Figure 25 shows the comparison of axial maximum vapor volume fraction at different water depths for original and optimal models. As clearly seen from Figure 24, the optimal model shows a good cavitation performance, not only at 50 m but also at 10, 25 and 100 m water depth, where the axial maximum vapor volume fraction is increased by 0.829%, 1.36%, and 17.6%, respectively. For the optimized nozzle, the change of the axial maximum vapor volume fraction at different water depths is the same as that of the non-optimized one, which decreases with the increasing water depth, and shows a small decrease first and then a larger decrease. However, the difference is that the target value is slightly decreased after optimization, and the reduction is significantly reduced compared with that before optimization, which also makes it possible to obtain a similar and higher axial maximum vapor volume fraction at different water depths. For the optimization effect, it is significantly increased as underwater depth increases. This is because in the shallow water, the cavitation effect is less affected by the environmental pressure, and the original target value is relatively large, so the optimization effect is not obvious. However, in deep water, due to the large environmental pressure, the cavitation effect of the non-optimized nozzle is largely suppressed. When it is optimized, the improvement is large and the optimization effect is obvious. From another point of view, the optimization of nozzle structure can reduce the suppression of cavitation due to the increase of environmental pressure to a certain extent, which significantly enhance the cavitation effect of the nozzle in deep water, and improve the working efficiency in deep water.

Regression analysis
Regression analysis is a statistical analysis method that determines the interdependence of the quantitative relationship between two or more variables. According to the number of variables involved, regression analysis is divided into one and multiple regression analysis. In order to design the optimal nozzle, it is very important to analyze the relationship between the geometry size and the vapor volume fraction. Therefore, the regression analysis method is used to analyze the data from Table 9.
At present, most of the researches only consider the maximum vapor volume fraction in the entire flow field, and the effect of each geometric parameter variable on this target is shown in Table 11. The Multiple R is 0.793 for the maximum vapor volume fraction, showing a moderate correlation with the maximum vapor volume fraction. Among them, the influence of the diffusion angle h is the largest, which is consistent with previous studies (Deng et al., 2008;Fu et al., 2013). It can be seen that the LHD method is credible for choosing samples in a fixed space.
The regression analysis for the axial maximum vapor volume fraction is shown in Table 12. The results obtained from Table 12 suggest that the Multiple R is 0.928. It can be seen that eight parameters are highly relevant with the axial maximum vapor volume fraction. The cylindrical section diameter, contraction angle and diffusion angle are beneficial for improving the axial maximum vapor volume fraction, in which the cylindrical section diameter has the greatest effect. The influence of the length of each part is generally less than the influence of the angle on the axial maximum vapor volume fraction, wherein the inlet section length has the greatest influence, and the effect of cylindrical section length is greater than contraction section length and diffusion segment length. However, the inlet diameter has the least effect on the axial maximum vapor volume fraction. It is worth noting that cylindrical section diameter, diffusion angle and contraction section have a large influence on both target values, in which the contraction section length affects the maximum vapor volume fraction, while its angle significantly affects the axial maximum vapor volume fraction.

Conclusion
This paper presents an cavitating nozzle structure optimization method. A hybrid algorithm is proposed to optimize the angle nozzle geometry. For the hybrid algorithm, a LHD method is applied to obtain an optimal global position for a non-linear optimization space, and a MMFD is used to search the local part in detail. After the optimization, the optimal model is obtained, which can increase the axial maximum vapor volume fraction by 9.41% at 50 m underwater depth. The cavitation effect between original and optimal models is compared at different water depths. The results indicate that the optimal model has a good cavitation effect at different water depths. The improvement of the cavitation effect after optimization at 100 m is the most obvious. The optimization of the nozzle structure can reduce the suppression of cavitation due to the increase of environmental pressure to a certain extent, significantly enhance the cavitation effect of the nozzle in deep water, and improve the working efficiency in deep water operation.
The relationship between the geometry size and the vapor volume fraction are discussed with regression analysis. The cylindrical section diameter, contraction angle and diffusion angle are beneficial for improving the axial maximum vapor volume fraction, in which the cylindrical section diameter has the greatest effect. The influence of the length of each part is generally less than the influence of angle on the axial maximum vapor volume fraction, wherein the inlet section length has the greatest influence, and the effect of cylindrical section length is larger than contraction section length and diffusion segment length. However, the inlet diameter has the least effect on the axial maximum vapor volume fraction. These results make the foundation for the future study on the optimal design of a cavitating nozzle structure and the development of water jet technology.
The optimization method presented in this paper only pays attention to the single target of the axial maximum volume fraction. In further research, multiple target values, including velocity, pressure, core constant velocity section length, and axial cavitation section length, will be considered simultaneously as a multi-objective optimization, so that the method can meet more complex application requirements. Combined with the actual working conditions, more external constraints, such as flow, pump pressure, target distance, can be added, not just considering their own structural parameters, so that the optimization results can better guide the actual operation and improve the operating efficiency.

Disclosure statement
No potential conflict of interest was reported by the authors.