CFD Investigations of Transient Cavitation Flows in Pipeline Based on Weakly-Compressible Model

: In hydraulic systems, transient ﬂow often occurs and may results in cavitation in pipelines. In this paper, the Computational Fluid Dynamics (CFD) method based on the Fluent software was used to investigate the cavitation ﬂow in pipeline; the density-pressure model was incorporated into the continuity equation by using further development of UDF (user deﬁned function), which reﬂects the variable wave speed of the transient cavitation ﬂow, and the related algorithms were established based on weakly compressible ﬂuid Reynolds Average Navier-Stokes (RANS) techniques. Firstly, the numerical simulations of the transient non-cavitation and cavitation ﬂows caused by the fast closing valve in the reservoir-pipe-valve system were carried out by using the grid slip technique. The simulation results can enrich the ﬂow ﬁeld information such as velocity, pressure and vapor volume fraction. Through the evolution process of the pressure ﬁeld, the propagation characteristics of pressure waves can be analyzed qualitatively and quantitatively. Through the evolution process of the velocity ﬁeld, it can be seen that the velocity distribution in the wall area changes rapidly and has a high gradient, which mainly depends on the viscosity. However, the change of the velocity distribution in the core region is related to the velocity distribution of the history of the past time, which mainly depends on the di ﬀ usion. The formation, development and collapse of the cavity can be successfully captured, and it can be clearly and visually observed that the uneven distribution of vapor cavity in the direction of pipe length and pipe diameter, and the vapor cavity move slowly along the top of the pipe wall. Rarefaction wave’s propagation into pressure decreasing region and pressure increasing region can lead to di ﬀ erent results of cavitation ﬂow. The accuracy and reliability of the weakly compressible ﬂuid RANS method were veriﬁed by comparing the calculated results with the experimental data.


Introduction
Hydraulic transients are caused by the rapid changes in pressurized pipelines and characterized by strong positive and negative pressures, which are often referred as water-hammer, and which probably cause device failures, system fatigues, leakages, or pipe ruptures, and can even be accompanied by cavitation when the liquid pressure in pipes drops below its corresponding vapor pressure, which may subsequently lead to more severe damages to the hydraulic system [1]. Such changes are normally triggered by valve opening/closing, pump startup/stop, load regulations of hydraulic turbine units, etc. Traditional models for one-dimensional (1D) hydraulic transients are widely used to estimate the hammer characteristics are revealed and some can demonstrate flow details in the transient events, cavitation is not considered in their simulations; thus, it is necessary to include the easily induced cavitation during water hammer in the calculation and study the mutual effect of cavitation and water hammer.
Unsteady friction is important in the prediction of water hammer decays. Ioriatti et al. [18] studied several unsteady friction models with the 1D and 2D finite volume schemes in transient flows, and concluded that the convolution integral models are significantly superior to instantaneous acceleration models concerning accuracy. Urbanowicz et al. [19] proposed a correction to the recursive formula of convolution integral model by Schohl [20] and the computation speed was improved. Vardy et al. [21] assessed the validity of frozen-viscosity conditions, which underpins the theoretical models of unsteady wall shear stress in transient flows in pipes and channels using detailed CFD simulations, and demonstrated that no frozen-viscosity distribution models perform well for large times after the commencement of acceleration while they perform well for short durations. This study used the CFD method to capture the axial velocity distribution during transient flows, and thus, the unsteady friction effect can be revealed sufficiently.
By using effective numerical simulation technology CFD to perform high-dimensional numerical simulation of transient cavitation flow, it can capture the detailed flow characteristics of complex pipeline hydraulic systems, and even related cavitation characteristics, and explore both the dynamic characteristics and internal laws of cavitation cavity. For calculations of transient flow using CFD, it was common to treat the problem with fluid medium of water as incompressible [22,23]. However, for fluid transient flow, especially with cavitation caused by pressure reduction, the propagation of variable pressure wave and cavitation simultaneously occur in the pipe; thus, the weak compressibility of water should be considered [10,24]. If the solution were directly solved based on the compressible method, the solution would be difficult to converge due to the huge difference between the speed of sound and the convection speed. Therefore, with the theory of weakly compressible fluids and the Mixture model for cavitation two-phase transient flows, based on the time-averaged continuity equation and the Reynolds time-averaged momentum equation, this paper introduced a fluid density-pressure model that can capture variable wave speed, and studied transient cavitation flow based on the Schnerr and Sauer cavitation model [25], which was used to calculate the mass transfer between the liquid phase and the vapor phase. In order to comprehensively and systematically investigate the applicability of aforementioned models and techniques, numerical simulations of transient non-cavitation flow and transient cavitation flow caused by the rapid valve closure of the reservoir-pipeline-valve system were both carried out. The time step and grid independence verifications, the analyses of pressure field, velocity field and cavitation cavity evolution process were carried out. The simulation results were analyzed and compared to the experimental data. This study was conducted by using the Fluent software (Version 14.5, ANSYS Inc., Pittsburgh, PA, USA), because it is a powerful CFD solver in which the density-pressure model can be incorporated using UDF so that the weakly compressible flow calculation method could be built.

Governing Equations of Weakly Compressible Fluid
In actual engineering, for the low-Mach number flow problem in which the fluid density change in the pressurized pipe is very small, if the solution were directly based on the compressible method, the solution would be difficult to converge due to the huge difference between the speed of sound and the convection speed. Thus, Song and Yuan [26] proposed a weakly compressible flow model for incompressible flow based on compressible models.
The continuity equation of the compressible model can be written as: where p is the pressure of the fluid and a is the pressure wave speed, thus, Equation (2) is a density-pressure model that is related to variable wave speed. Substituting Equation (2) into Equation (1) yields: , a * = a a 0 , then Equation (1) is converted to dimensionless as: where Strauhal number St = l 0 u 0 t 0 , Mach number Ma = u 0 a , p 0 is far from the flow field pressure, u 0 , t 0 , and l 0 are the reference flow rate, time, and length, respectively.
For weakly compressible fluid, Ma < 1, Ma 2 1, St 1, thus, the last term of Equation (4) can be omitted. Therefore, the continuous equation of weakly compressible fluid is: The momentum equation of weakly compressible fluid is: where µ is the dynamic viscosity of the fluid and f i is the i-direction component of body force. Equations (5) and (6) constitute the basic governing equations for weakly compressible fluids.

Solution of RANS Equation Based on the Weakly Compressible Model
RANS equations are used to describe the governing equations of weakly compressible fluids: where the over bar indicates the time average components and τ ij represents Reynolds stress and τ ij = −ρu i u j . There are 11 unknown variables in Equations (7) and (8), but there are only four equations. Therefore, the equations are not closed. Choosing different turbulence models to close the equations is the key to solving. The commonly used turbulence models are: standard k − ε model, RNG k − ε model, Realizable k − ε model, SST k − ω model etc. The standard two-equation eddy viscosity models (k − ε and k − ε types of model) based on the Boussinesq relationship τ ij = −ρu i u j = −2/3(ρk + µ t ∂u i /∂x i )δ ij + µ t (∂u i /∂x j + ∂u j /∂x i ) (where k is the turbulence kinetic energy, δ ij is the Kronecker delta, and µ t is the turbulent viscosity) were originally developed for single-phase non-cavitation flows. When simulating the sheet cavitation that appears on the hydrofoil with the above models, there is a tendency to overestimate the turbulent viscosity on the cavity collapse and on its downstream region. Thus, they are not suitable for the simulation of re-injection flow and cavity structure shedding process [27][28][29][30]. However, the simulation experience of transient cavitation flow in the pipeline is insufficient, and it is not clear whether the cavity inside the pipe has similar characteristic sheet cavitation. Therefore, the RNG k − ε model that can better simulate curvature is used, and the transport equations of the turbulence kinetic energy, k, and the turbulence dissipation rate, ε, are as follows: where µ t = ρC µ k 2 ε and G k is the generation of turbulence kinetic energy due to the mean velocity gradients, computed by ; α k and α ε are the turbulent Prandtl number for k and ε, and α k = α ε = 1.39. The other model constants are, respectively: C 1ε

Cavitation Model
At present, the single-fluid Mixture model based on the multi-phase flow Euler-Euler method [31] is often used to model the cavitation flow. The Schnerr and Sauer cavitation model was used in this paper [25]. The liquid-vapor mass transfer (evaporation and condensation) are governed by the vapor transfer equation: where subscript v represents the vapor phase, α is the vapor volume fraction, ρ v is the vapor density, u vi is the i-direction component velocity of vapor, R e is the phase change rate from liquid phase to vapor phase, and R c is the phase change rate from vapor phase to liquid phase. R e and R c are obtained based on the Rayleigh-Plesset equation as follows: When where subscript l represents liquid phase, subscript m represents mixture, and p is the local pressure, p v is the saturation vapor pressure in the bubble, R B is the bubble radius, and R B = ( α 3 , in which n is the bubble number density.

Computational Domains and Grid Model
The experimental data used in this paper are from the pressure fluctuation data of the transient flow experiments conducted by Simpson [32]. The experiment device consists of two closed water reservoirs, a slope copper pipe with a straight upward direction, a valve with adjustable closing time, and a pressure detection system, as is shown in Figure 1. In the experiment device, the water level and water head of the closed water reservoirs at both ends could be adjusted, and the rapid closing of the valve caused transient flow phenomena. There were pressure sensors installed at 1/4 of the pipe length to the upstream reservoir, 1/4 of the pipe length to the downstream reservoir, and at the valve to measure pressure fluctuations.
The experimental data used in this paper are from the pressure fluctuation data of the transient flow experiments conducted by Simpson [32]. The experiment device consists of two closed water reservoirs, a slope copper pipe with a straight upward direction, a valve with adjustable closing time, and a pressure detection system, as is shown in Figure 1. In the experiment device, the water level and water head of the closed water reservoirs at both ends could be adjusted, and the rapid closing of the valve caused transient flow phenomena. There were pressure sensors installed at 1/4 of the pipe length to the upstream reservoir, 1/4 of the pipe length to the downstream reservoir, and at the valve to measure pressure fluctuations. The diameter of the pipeline was 19.05 mm, the total length was 36.0 m, and the height difference of the upstream and downstream reservoirs was 1 m. The main parameters of different experimental schemes are shown in Table 1. Among them, experiment 1 is a transient non-cavitation flow (the wave speed in the pipe is 1280 m/s), experiment 2 is a transient cavitation flow, and both of them are turbulent flows. In this paper, the water hammers with initial steady state of turbulent flow were discussed, and the water hammer with initial Reynolds number smaller than 2320 of laminar flow will be studied in the future to fulfill the calculation model. Based on the experimental data provided by Simpson [32], the model was built using Proe software, and the geometry model sketch is shown in Figure 2. Since the pressure at the inlet and outlet of the pipeline was constant during the experiments, the reservoirs were not modeled but simplified with the given pressure boundary conditions. The positions of the pressure sensors at Sl and S2 of the pipeline were respectively located at 9.0 m and 27.0 m from upstream inlet, and the position of the pressure sensor S3 was 0.02 m upstream to the valve. The ball valve was modeled using a rotating cylinder with the same diameter of the pipe. The diameter of the pipeline was 19.05 mm, the total length was 36.0 m, and the height difference of the upstream and downstream reservoirs was 1 m. The main parameters of different experimental schemes are shown in Table 1. Among them, experiment 1 is a transient non-cavitation flow (the wave speed in the pipe is 1280 m/s), experiment 2 is a transient cavitation flow, and both of them are turbulent flows. In this paper, the water hammers with initial steady state of turbulent flow were discussed, and the water hammer with initial Reynolds number smaller than 2320 of laminar flow will be studied in the future to fulfill the calculation model. Based on the experimental data provided by Simpson [32], the model was built using Proe software, and the geometry model sketch is shown in Figure 2. Since the pressure at the inlet and outlet of the pipeline was constant during the experiments, the reservoirs were not modeled but simplified with the given pressure boundary conditions. The positions of the pressure sensors at S l and S 2 of the pipeline were respectively located at 9.0 m and 27.0 m from upstream inlet, and the position of the pressure sensor S 3 was 0.02 m upstream to the valve. The ball valve was modeled using a rotating cylinder with the same diameter of the pipe. Due to the simple structure of the pipeline, it was easy to implement structured meshing, thus, the structure grid was used in the model. The calculations were based on a 2D geometry model of the longitudinal section of the pipeline, and the grid models were also 2D, as shown in Figure 3. The quarter-turn ball valve was from the experiment of Simpson [32]. Due to the simple structure of the pipeline, it was easy to implement structured meshing, thus, the structure grid was used in the model. The calculations were based on a 2D geometry model of the longitudinal section of the pipeline, and the grid models were also 2D, as shown in Figure 3. The quarter-turn ball valve was from the experiment of Simpson [32]. Due to the simple structure of the pipeline, it was easy to implement structured meshing, thus, the structure grid was used in the model. The calculations were based on a 2D geometry model of the longitudinal section of the pipeline, and the grid models were also 2D, as shown in Figure 3. The quarter-turn ball valve was from the experiment of Simpson [32].

Computational Scheme and Boundary Conditions
In the calculations, the pressure-based algorithm was used, and the coupled method was used as the pressure-velocity coupling method, with the spatial discretization of pressure in PRESTO! scheme and spatial discretization of others in QUICK scheme. In the transient calculations, the first order implicit discrete scheme for the time term was used.
In correspondence to the experiments, the boundary conditions of the calculation cases are listed in Table 2. The inlet of the pipeline was set as the pressure inlet boundary, the outlet of the pipeline was set as the pressure outlet boundary, the wall surfaces were set to be stationary and no slip boundary, and the friction coefficient of the walls was default. Since the ball valve was closed quickly, it was considered that the rotational angular velocity of the ball valve was linear to time and the motion of the valve was realized by a user defined function (UDF). The contact spherical surfaces of the ball valve with the left and right pipes were defined as two pairs of interfaces. Gravity was considered in the calculations. By using the Fluent software, the steady state calculation was firstly carried out. After the iterative convergence, the mean flow velocity of each section in the pipeline was approximately equal to the experimental value, and the initial steady state of the flow field before the transient was obtained.
Then the transient calculation was carried out. In order to describe the compressibility of water, it is necessary to introduce a water pressure equation that characterizes the change of water density with pressure, which is, to load the following water physical property UDF:

Computational Scheme and Boundary Conditions
In the calculations, the pressure-based algorithm was used, and the coupled method was used as the pressure-velocity coupling method, with the spatial discretization of pressure in PRESTO! scheme and spatial discretization of others in QUICK scheme. In the transient calculations, the first order implicit discrete scheme for the time term was used.
In correspondence to the experiments, the boundary conditions of the calculation cases are listed in Table 2. The inlet of the pipeline was set as the pressure inlet boundary, the outlet of the pipeline was set as the pressure outlet boundary, the wall surfaces were set to be stationary and no slip boundary, and the friction coefficient of the walls was default. Since the ball valve was closed quickly, it was considered that the rotational angular velocity of the ball valve was linear to time and the motion of the valve was realized by a user defined function (UDF). The contact spherical surfaces of the ball valve with the left and right pipes were defined as two pairs of interfaces. Gravity was considered in the calculations. By using the Fluent software, the steady state calculation was firstly carried out. After the iterative convergence, the mean flow velocity of each section in the pipeline was approximately equal to the experimental value, and the initial steady state of the flow field before the transient was obtained.
Then the transient calculation was carried out. In order to describe the compressibility of water, it is necessary to introduce a water pressure equation that characterizes the change of water density with pressure, which is, to load the following water physical property UDF: where ρ l0 and K l are the density and bulk modulus of the liquid at absolute reference pressure p 0 , respectively, and ∆p = p − p 0 .

Grid Independence Verification
In order to explore the influence of the number of grids on the numerical simulation results, several sets of grids were used to carry out the grid independence trials. In the grid models, there was a contraction in the wall area, with the distance of the first grid node to the wall 0.5 mm, the grid spacing growth law in the radial direction from the wall to the core area was exponential with the ratio of 1.05, and the radial grid element number in this direction was 28. Since the grid on the radial direction was dense enough, the radial grid spacing of the pipeline stayed unchanged. Since the pipeline was long, the axial grid spacing varied with four numbers (1 mm, 3 mm, 5 mm, and 8 mm), and the time step of the transient simulation was uniformly set to be 0.0001 s. The cross-sectional predicted average pressure values at the S 3 monitoring position are shown in Figure 4.

Grid Independence Verification
In order to explore the influence of the number of grids on the numerical simulation results, several sets of grids were used to carry out the grid independence trials. In the grid models, there was a contraction in the wall area, with the distance of the first grid node to the wall 0.5 mm, the grid spacing growth law in the radial direction from the wall to the core area was exponential with the ratio of 1.05, and the radial grid element number in this direction was 28. Since the grid on the radial direction was dense enough, the radial grid spacing of the pipeline stayed unchanged. Since the pipeline was long, the axial grid spacing varied with four numbers (1 mm, 3 mm, 5 mm, and 8 mm), and the time step of the transient simulation was uniformly set to be 0.0001 s. The cross-sectional predicted average pressure values at the S3 monitoring position are shown in Figure 4. It can be seen from Figure 4 that the pressure fluctuation curves simulated by different grid models are almost identical. Moreover, the difference between them is less than 1% and they are in good agreement with the experimental data. Therefore, within this range, the number of grids has little influence on the calculation results of transient flow pressure fluctuation. In order to ensure the numerical calculation efficiency and accuracy, a grid model with an axial grid spacing of 5 mm was used, and the number of grids was 248,610. The final y+ values of the pipe wall were 5.6 and 4.4, respectively, in case 1 and case 2, the requirements were met. It can be seen from Figure 4 that the pressure fluctuation curves simulated by different grid models are almost identical. Moreover, the difference between them is less than 1% and they are in good agreement with the experimental data. Therefore, within this range, the number of grids has little influence on the calculation results of transient flow pressure fluctuation. In order to ensure the numerical calculation efficiency and accuracy, a grid model with an axial grid spacing of 5 mm was used, and the number of grids was 248,610. The final y+ values of the pipe wall were 5.6 and 4.4, respectively, in case 1 and case 2, the requirements were met.

Time Step Independence Verification
In addition to considering the effect of grid spacing on transient simulation results, the impact of time step on the simulation results was also studied. In the transient simulation, the number of grids was 248,610. The predicted cross-sectional average pressures at the S 3 monitoring position with different time steps (0.00001 s, 0.0005 s, 0.0001 s, and 0.001 s) are shown in Figure 5. It can be seen from Figure 5 that the maximum water hammer pressure and pressure fluctuation cycle are hardly affected by the time step. When the time step is smaller than 0.0001 s, the results are hardly affected by the time step. In order to ensure the numerical calculation efficiency and calculation accuracy, the final time step was chosen to be 0.0001 s. The C-F-L criterion was fulfilled with the selected time step size in this paper.
Water 2018, 10, x FOR PEER REVIEW 9 of 29

Time Step Independence Verification
In addition to considering the effect of grid spacing on transient simulation results, the impact of time step on the simulation results was also studied. In the transient simulation, the number of grids was 248,610. The predicted cross-sectional average pressures at the S3 monitoring position with different time steps (0.00001 s, 0.0005 s, 0.0001 s, and 0.001 s) are shown in Figure 5. It can be seen from Figure 5 that the maximum water hammer pressure and pressure fluctuation cycle are hardly affected by the time step. When the time step is smaller than 0.0001 s, the results are hardly affected by the time step. In order to ensure the numerical calculation efficiency and calculation accuracy, the final time step was chosen to be 0.0001 s. The C-F-L criterion was fulfilled with the selected time step size in this paper.

Analysis of the Pressure Fluctuation
Firstly, numerical simulation was carried out for experiment 1 of initial flow velocity of 0.239

Analysis of the Pressure Fluctuation
Firstly, numerical simulation was carried out for experiment 1 of initial flow velocity of 0.239 m/s. The transient flow process was caused by rapid closure of the valve, and the predicted cross-sectional average pressure values at S l , S 2 , and S 3 monitoring positions of each time point are shown in Figure 6, and compared to the experimental data.

Analysis of the Pressure Fluctuation
Firstly, numerical simulation was carried out for experiment 1 of initial flow velocity of 0.239 m/s. The transient flow process was caused by rapid closure of the valve, and the predicted crosssectional average pressure values at Sl, S2, and S3 monitoring positions of each time point are shown in Figure 6, and compared to the experimental data. It can be seen from Figure 6 that at the valve (the S3 monitoring point), the first peak value of the pressure fluctuation can be accurately predicted. However, the attenuation of the pressure in the simulations is lighter than the experimental result, after 0.465 s, the peak value of the experimental pressure head decreases from 55.33 m to 49.77 m due to the frictional resistance, and the calculated pressure head peak is 55.51 m. The discrepancy between the CFD results and experimental data after the first peak value of the pressure fluctuation is because the energy dissipation is underestimated during the transient process, which may be because the pipe friction effect in the transient process cannot be well simulated in the calculations or the turbulent dissipation is underestimated in the used turbulence model. One reason for the discrepancy may come from the unsteady friction models (i.e., the convolutional integral model and the instantaneous acceleration model) effects in the 1D model which are also not revealed well in the CFD simulations due to their complex mechanisms. Another reason for the discrepancy may be the fact that the 2D simulation cannot reveal well the 3D water hammer decay. The results at the other monitoring points are similar to those at the valve. And the simulation results can match well with the experimental data at each monitoring position. It can be seen from Figure 6 that at the valve (the S 3 monitoring point), the first peak value of the pressure fluctuation can be accurately predicted. However, the attenuation of the pressure in the simulations is lighter than the experimental result, after 0.465 s, the peak value of the experimental pressure head decreases from 55.33 m to 49.77 m due to the frictional resistance, and the calculated pressure head peak is 55.51 m. The discrepancy between the CFD results and experimental data after the first peak value of the pressure fluctuation is because the energy dissipation is underestimated during the transient process, which may be because the pipe friction effect in the transient process cannot be well simulated in the calculations or the turbulent dissipation is underestimated in the used turbulence model. One reason for the discrepancy may come from the unsteady friction models (i.e., the convolutional integral model and the instantaneous acceleration model) effects in the 1D model which are also not revealed well in the CFD simulations due to their complex mechanisms. Another reason for the discrepancy may be the fact that the 2D simulation cannot reveal well the 3D water hammer decay. The results at the other monitoring points are similar to those at the valve. And the simulation results can match well with the experimental data at each monitoring position.

Analysis of the Pressure and Velocity Fields
Figures 7 and 8 illustrate the complete evolutionary cycle of the pressure and velocity fields in the axial plane of the pipe (corresponding to time nodes of T 0 ∼ T 12 on Figure 6b). Since the transient non-cavitation flow of experiment 1 is single-phase, the pressure wave speed a is constant. The majority of the mean velocity reduction in the pipe occurs in the last one-third of time of the valve closure [32]; the initial closure of the valve does not induce much velocity reduction and the pressure increase is moderate, while the rapid change of velocity in the last part of the closure time induces a big pressure increase, and the big pressure and velocity change can be shown in the pressure field contour more clearly. Thus, the time nodes T 1 ∼ T 12 of the pressure and velocity fields can be expressed in terms of the valve closure time, Tc, and the one-way travel time of the pressure wave in the pipe, L/a, as shown in Table 3. After the above mentioned first fluctuation cycle, the propagation process will be repeated continuously and periodically, and finally, under the action of the frictional resistance, the fluctuation gradually weakens.   Figure 9 is the experimental data and numerical simulation results of pressure head fluctuation at S1 monitoring point, with marked time nodes 0~11 (different from time nodes in Figure 6b). The velocity distribution at S1 monitoring position and the corresponding pressure field on the axis plane of the pipeline at these time nodes are shown and analyzed in five stages, and the stages are divided

Analysis of the Velocity Distribution Change
The pictures in Figures 7 and 8 are magnified 186 times in radial direction of the slender pipe so that the results can be demonstrated clearly. During the transient flow in the long and slender pipeline, with the pipe inclination of 1.56 degrees, the biggest pressure caused by gravity is 1.5% compared to the pressure peak value due to water hammer, and the pressure gradient caused by gravity is also very small; thus, the pressure and pressure gradient caused by gravity are not apparently shown. The pressure change in the pipeline is mainly driven by the pressure wave propagation. Moreover, the pressure contours are seemingly to be with vertical lines. A fluctuation cycle is divided into four stages: (1) The first stage: the pressure increase wave propagates from the valve to the upstream reservoir. When the valve is quickly closed, the small range of flow region immediately adjacent to the valve stops first, and a pressure increase wave ∆p 1 (the pressure head of ∆p 1 is 31.39 m) is caused at that point. Figure 7a-d show that the pressure increase wave propagates at wave speed a, in which for transient non-cavitation flow the wave speed a is calculated with a = K l /ρ l0 and the result is a = 1280 m/s, and the mean flow velocity at the region where the wave propagates change to zero (Figure 8a-d).
At time t = L/a (L is the length of the pipeline), the pressure increase wave propagates to the upstream reservoir, and the entire pipeline is in a pressurized state.
(2) The second stage: at the end of the first stage, the fluid pressure in the pipeline is bigger than the pressure in the upstream reservoir, under the action of the pressure difference, a flow velocity opposite to the original one starts from the upstream reservoir, and causes pressure of the corresponding part of the pipe to decrease almost to the original value (considering the pressure decay due to dissipation, this value is smaller than the original one). That is, the pressure decrease wave propagates from the upstream reservoir to the valve at the wave speed a, as shown in Figure 7e-g. At time t = 2L/a, the pressure decrease wave propagates to the valve, and the pressure of the entire pipeline restores to a value a little smaller than the original one.
(3) The third stage: at time t = 2L/a, the entire pipeline has a flow velocity opposite to the original one, but at this time the valve is fully closed, and the velocity of the small range of flow region immediately adjacent to the valve becomes zero, which causes a pressure decrease wave ∆p 2 (the pressure head of ∆p 2 is −30.9 m, and the pressure head difference from the original steady state value at this point is less than the original pressure increase value ∆p 1 due to energy dissipation in the pipe), and Figure 7h-j show the process of the pressure decrease wave propagation from the valve to the upstream reservoir. At time t = 3L/a, the pressure decrease wave propagates to the upstream water reservoir, and the entire pipeline is in a decompressed state.
(4) The fourth stage: at time t = 3L/a, the fluid pressure in the entire pipeline is smaller than the pressure in the upstream water reservoir, under the action of the pressure difference, a flow velocity with the same direction to the original one starts from the upstream reservoir, and makes the pressure of corresponding part of the pipe return to a value a little smaller than the original one. That is, the pressure increase wave propagates from the upstream reservoir to the valve at the wave speed a, as shown in Figure 7k-m. At time t = 4L/a, the pressure increase wave propagates to the valve, the fluid pressure in the entire pipeline restores to a value a little smaller than the original one and the fluid velocity direction in the entire pipeline is from the reservoir to the valve, and the flow state in the entire pipeline is similar to that at the beginning of the first stage.
After the above mentioned first fluctuation cycle, the propagation process will be repeated continuously and periodically, and finally, under the action of the frictional resistance, the fluctuation gradually weakens.  Stage 1: Figures 10 and 11 are respectively the axial velocity distribution and local pressure distribution at the S1 monitoring position of different time, and Figure 12 is the pressure field evolution on the axis plane of the pipeline in this stage. Before the valve is closed, the velocity at 0 = 0 time at S1 monitoring position exhibits a fully developed steady-state velocity distribution, with ̅ = 0.239 m/s, the pressure field at this time in Figure 12a shows the steady-state pressure field, with = 24.24 m, and local pressure distribution at S1 monitoring position in Figure 11a shows that the pressure gradient on this line is almost 0. When the valve is closed, the pressure increase wave propagates from the valve to the upstream reservoir (Figure 12b, 1 = 0.044 s = + 0.71 × / , = 37.29 m), the flow velocity in the pipeline gradually decreases (Figure 10b, ̅ = 0.137 m/s); the local pressure distribution at S1 monitoring position (Figure 11b) shows that when the pressure wave passes, the pressure at the upper part of the pipe is bigger, while the pressure at the lower part of the pipe is smaller with the biggest pressure difference 27.7 Pa, but this value is still very small compared to the average pressure 365,126 Pa on S1 monitoring position and the biggest relative pressure difference is 27.7/365126 = 0.008%, thus, the 1D and quasi-2D model assumption that pressure on the cross-section is constant is still valid. After the pressure increase wave passes (Figures 10c and  12c, 2 = 0.05 = + 0.92 × / ), the pressure at S1 position reaches the maximum (pressure head = 55.1 m), the average flow velocity is zero, there is a significant reverse flow near the wall and the velocity gradient is large, the velocity distribution of the core region remains the same shape as the initial steady state, and the local pressure at S1 monitoring position is almost constant at this time (Figure 11c). Stage 1: Figures 10 and 11 are respectively the axial velocity distribution and local pressure distribution at the S 1 monitoring position of different time, and Figure 12 is the pressure field evolution on the axis plane of the pipeline in this stage. Before the valve is closed, the velocity at T 0 = 0 s time at S 1 monitoring position exhibits a fully developed steady-state velocity distribution, with U = 0.239 m/s, the pressure field at this time in Figure 12a shows the steady-state pressure field, with H = 24.24 m, and local pressure distribution at S 1 monitoring position in Figure 11a shows that the pressure gradient on this line is almost 0. When the valve is closed, the pressure increase wave propagates from the valve to the upstream reservoir (Figure 12b, T 1 = 0.044 s = Tc + 0.71 × L/a, H = 37.29 m), the flow velocity in the pipeline gradually decreases (Figure 10b, U = 0.137 m/s); the local pressure distribution at S 1 monitoring position (Figure 11b) shows that when the pressure wave passes, the pressure at the upper part of the pipe is bigger, while the pressure at the lower part of the pipe is smaller with the biggest pressure difference 27.7 Pa, but this value is still very small compared to the average pressure 365,126 Pa on S 1 monitoring position and the biggest relative pressure difference is 27.7/365126 = 0.008%, thus, the 1D and quasi-2D model assumption that pressure on the cross-section is constant is still valid. After the pressure increase wave passes (Figures 10c and 12c, T 2 = 0.05 s = Tc + 0.92 × L/a), the pressure at S 1 position reaches the maximum (pressure head H = 55.1 m), the average flow velocity is zero, there is a significant reverse flow near the wall and the velocity gradient is large, the velocity distribution of the core region remains the same shape as the initial steady state, and the local pressure at S 1 monitoring position is almost constant at this time (Figure 11c). cross-section is constant is still valid. After the pressure increase wave passes (Figures 10c and  12c, 2 = 0.05 = + 0.92 × / ), the pressure at S1 position reaches the maximum (pressure head = 55.1 m), the average flow velocity is zero, there is a significant reverse flow near the wall and the velocity gradient is large, the velocity distribution of the core region remains the same shape as the initial steady state, and the local pressure at S1 monitoring position is almost constant at this time (Figure 11c).  Stage 2: Figures 13 and 14 are respectively the axial velocity distribution and local pressure distribution at the S1 monitoring position of different time, and Figure 15 is the pressure field evolution on the axis plane of the pipeline in this stage. When the pressure decrease wave generated at the upstream reservoir propagates to the downstream valve (Figure 15a, 3 = 0.058 s = + 1.21 × / ), the pressure head at S1 position decreases to 41.29 m, the velocity direction is from the valve to the reservoir (Figure 13a, ̅ = −0.102 m/s), the reverse flow develops from the wall region and then diffuses to the core region, and the local pressure distribution in Figure 14a shows that the pressure at the upper part of the section is bigger while the pressure at the lower part of the section is smaller with a biggest pressure difference of 22.5 Pa, but the average pressure on S1 monitoring position at this time is 410,492 Pa and the biggest relative pressure difference is 22.5/410492 = 0.005% (still very small). After the pressure wave passes (Figure 15b (Figure 14b). In this process, the velocity distribution of the core region at S1 position stays the same shape as the initial state, and the velocity magnitude at the wall region is bigger than that of the core region. Stage 2: Figures 13 and 14 are respectively the axial velocity distribution and local pressure distribution at the S1 monitoring position of different time, and Figure 15 is the pressure field evolution on the axis plane of the pipeline in this stage. When the pressure decrease wave generated at the upstream reservoir propagates to the downstream valve (Figure 15a, 3 = 0.058 s = + 1.21 × / ), the pressure head at S1 position decreases to 41.29 m, the velocity direction is from the valve to the reservoir (Figure 13a, ̅ = −0.102 m/s), the reverse flow develops from the wall region and then diffuses to the core region, and the local pressure distribution in Figure 14a shows that the pressure at the upper part of the section is bigger while the pressure at the lower part of the section is smaller with a biggest pressure difference of 22.5 Pa, but the average pressure on S1 monitoring position at this time is 410,492 Pa and the biggest relative pressure difference is 22.5/410492 = 0.005% (still very small). After the pressure wave passes (Figure 15b (Figure 14b). In this process, the velocity distribution of the core region at S1 position stays the same shape as the initial state, and the velocity magnitude at the wall region is bigger than that of the core region. Stage 2: Figures 13 and 14 are respectively the axial velocity distribution and local pressure distribution at the S 1 monitoring position of different time, and Figure 15 is the pressure field evolution on the axis plane of the pipeline in this stage. When the pressure decrease wave generated at the upstream reservoir propagates to the downstream valve (Figure 15a, T 3 = 0.058 s = Tc + 1.21 × L/a), the pressure head at S 1 position decreases to 41.29 m, the velocity direction is from the valve to the reservoir (Figure 13a, U = −0.102 m/s), the reverse flow develops from the wall region and then diffuses to the core region, and the local pressure distribution in Figure 14a shows that the pressure at the upper part of the section is bigger while the pressure at the lower part of the section is smaller with a biggest pressure difference of 22.5 Pa, but the average pressure on S 1 monitoring position at this time is 410,492 Pa and the biggest relative pressure difference is 22.5/410492 = 0.005% (still very small). After the pressure wave passes (Figure 15b, T 4 = 0.064 s = Tc + 1.42 × L/a), the pressure head decreases to 24.46 m, the average axial velocity U = −0.236 m/s, and the pressure at S 1 monitoring position at this time is almost constant with a very small pressure gradient (Figure 14b). In this process, the velocity distribution of the core region at S 1 position stays the same shape as the initial state, and the velocity magnitude at the wall region is bigger than that of the core region.

Analysis of the Velocity Distribution Change
pressure at the upper part of the section is bigger while the pressure at the lower part of the section is smaller with a biggest pressure difference of 22.5 Pa, but the average pressure on S1 monitoring position at this time is 410,492 Pa and the biggest relative pressure difference is 22.5/410492 = 0.005% (still very small). After the pressure wave passes (Figure 15b, 4 = 0.064 s = + 1.42 × / ), the pressure head decreases to 24.46 m, the average axial velocity ̅ = −0.236 m/s, and the pressure at S1 monitoring position at this time is almost constant with a very small pressure gradient (Figure 14b). In this process, the velocity distribution of the core region at S1 position stays the same shape as the initial state, and the velocity magnitude at the wall region is bigger than that of the core region.  Stage 3: The pressure decrease wave generated at the valve propagates upstream, Figures 16 and  17 are respectively the axial velocity distribution and local pressure distribution at the S1 monitoring position of different time, and Figure 18 is the pressure field evolution on the axis plane of the pipeline in this stage. At 5 = 0.092 s = + 2.42 × / , the pressure head at S1 position is 24.08 m, and the velocity distribution in Figure 16a exhibits that the velocity at the wall region firstly begins to change (with cross-sectional average velocity of ̅ = −0.233 m/s ), and the velocity change gradually diffuses to the core region, at 6 = 0.1 s = + 2.7 × / , ̅ = −0.14 m/s. After the pressure wave passes (Figure 18c, 7 = 0.108 s = + 2.99 × / ), the pressure head at S1 monitoring position decreases to the minimum, −5.57 m, and the axial average velocity is 0 m/s, the velocity distribution in Figure 16c is similar to that in Figure 10c. In the process, the velocity distribution of the core region remains the same to the initial shape, but the velocity peak clearly spreads to the core region. The local pressure distribution shows that the pressure difference on the S1 monitoring position is very small before the pressure wave passes (Figure 17a), the pressure difference is bigger (15.3 Pa) when the pressure wave passes (Figure 17b) but the average pressure on S1 monitoring position at this time is 114,364 Pa and the biggest relative pressure difference is 15.3/114364 = 0.01% (still very small), and the pressure difference is also not big after the pressure wave passes (Figure 17c). Stage 3: The pressure decrease wave generated at the valve propagates upstream, Figures 16 and  17 are respectively the axial velocity distribution and local pressure distribution at the S1 monitoring position of different time, and Figure 18 is the pressure field evolution on the axis plane of the pipeline in this stage. At 5 = 0.092 s = + 2.42 × / , the pressure head at S1 position is 24.08 m, and the velocity distribution in Figure 16a exhibits that the velocity at the wall region firstly begins to change (with cross-sectional average velocity of ̅ = −0.233 m/s ), and the velocity change gradually diffuses to the core region, at 6 = 0.1 s = + 2.7 × / , ̅ = −0.14 m/s. After the pressure wave passes (Figure 18c, 7 = 0.108 s = + 2.99 × / ), the pressure head at S1 monitoring position decreases to the minimum, −5.57 m, and the axial average velocity is 0 m/s, the velocity distribution in Figure 16c is similar to that in Figure 10c. In the process, the velocity distribution of the core region remains the same to the initial shape, but the velocity peak clearly spreads to the core region. The local pressure distribution shows that the pressure difference on the S1 monitoring position is very small before the pressure wave passes (Figure 17a), the pressure difference is bigger (15.3 Pa) when the pressure wave passes (Figure 17b) but the average pressure on S1 monitoring position at this time is 114,364 Pa and the biggest relative pressure difference is 15.3/114364 = 0.01% (still very small), and the pressure difference is also not big after the pressure wave passes (Figure 17c).  Figure 18 is the pressure field evolution on the axis plane of the pipeline in this stage. At T 5 = 0.092 s = Tc + 2.42 × L/a, the pressure head at S 1 position is 24.08 m, and the velocity distribution in Figure 16a exhibits that the velocity at the wall region firstly begins to change (with cross-sectional average velocity of U = −0.233 m/s), and the velocity change gradually diffuses to the core region, at T 6 = 0.1 s = Tc + 2.7 × L/a, U = −0.14 m/s. After the pressure wave passes (Figure 18c, T 7 = 0.108 s = Tc + 2.99 × L/a), the pressure head at S 1 monitoring position decreases to the minimum, −5.57 m, and the axial average velocity is 0 m/s, the velocity distribution in Figure 16c is similar to that in Figure 10c. In the process, the velocity distribution of the core region remains the same to the initial shape, but the velocity peak clearly spreads to the core region. The local pressure distribution shows that the pressure difference on the S 1 monitoring position is very small before the pressure wave passes (Figure 17a), the pressure difference is bigger (15.3 Pa) when the pressure wave passes (Figure 17b) but the average pressure on S 1 monitoring position at this time is 114,364 Pa and the biggest relative pressure difference is 15.3/114364 = 0.01% (still very small), and the pressure difference is also not big after the pressure wave passes (Figure 17c). remains the same to the initial shape, but the velocity peak clearly spreads to the core region. The local pressure distribution shows that the pressure difference on the S1 monitoring position is very small before the pressure wave passes (Figure 17a), the pressure difference is bigger (15.3 Pa) when the pressure wave passes (Figure 17b) but the average pressure on S1 monitoring position at this time is 114,364 Pa and the biggest relative pressure difference is 15.3/114364 = 0.01% (still very small), and the pressure difference is also not big after the pressure wave passes (Figure 17c).  Stage 4: Figures 19 and 20 are respectively the axial velocity distribution and local pressure distribution at the S1 monitoring position of different time, and Figure 21 is the pressure field evolution on the axis plane of the pipeline in this stage. When the pressure increase wave generated at the upstream reservoir in the fourth stage reaches the S1 monitoring position (Figure 21a, 8 = 0.116 s = + 3.27 × / ), the pressure head at S1 position is 12.91 m, the velocity at the wall region is positive Figure 19a, with average axial velocity ̅ = 0.148 m/s; the local pressure difference is 14.5 Pa (Figure 20a), the average pressure on S1 monitoring position at this time is 126,411 Pa, and the biggest relative pressure difference is 14.5/126411 = 0.01% (still very small). After the pressure wave passes (Figure 21b, 9 = 0.124 s = + 3.56 × / ), the pressure head is 24.15 m, ̅ = 0.233 m/s, and the local pressure on S1 monitoring position is almost constant (Figure 20b). Finally, the velocity distribution in Figure 19b will be returned to the initial steady state. Stage 4: Figures 19 and 20 are respectively the axial velocity distribution and local pressure distribution at the S1 monitoring position of different time, and Figure 21 is the pressure field evolution on the axis plane of the pipeline in this stage. When the pressure increase wave generated at the upstream reservoir in the fourth stage reaches the S1 monitoring position (Figure 21a, 8 = 0.116 s = + 3.27 × / ), the pressure head at S1 position is 12.91 m, the velocity at the wall region is positive Figure 19a, with average axial velocity ̅ = 0.148 m/s; the local pressure difference is 14.5 Pa (Figure 20a), the average pressure on S1 monitoring position at this time is 126,411 Pa, and the biggest relative pressure difference is 14.5/126411 = 0.01% (still very small). After the pressure wave passes (Figure 21b, 9 = 0.124 s = + 3.56 × / ), the pressure head is 24.15 m, ̅ = 0.233 m/s, and the local pressure on S1 monitoring position is almost constant (Figure 20b). Finally, the velocity distribution in Figure 19b will be returned to the initial steady state.  (Figure 20b). Finally, the velocity distribution in Figure 19b will be returned to the initial steady state. Pa (Figure 20a), the average pressure on S1 monitoring position at this time is 126,411 Pa, and the biggest relative pressure difference is 14.5/126411 = 0.01% (still very small). After the pressure wave passes (Figure 21b, 9 = 0.124 s = + 3.56 × / ), the pressure head is 24.15 m, ̅ = 0.233 m/s, and the local pressure on S1 monitoring position is almost constant (Figure 20b). Finally, the velocity distribution in Figure 19b will be returned to the initial steady state.  Stage 5: The wave cycle will be repeated in the pipeline due to inertia, in this stage, the pressure increase wave generated at the valve propagates upstream. Figure 22 is the axial velocity distribution at S1 monitoring position and Figure 23 is the corresponding local pressure distribution. At( 10 = 0.148 s = + 4.41 × / ), the pressure head at S1 position is 24.85 m, a little bigger than that of the original steady state due to the effect of the first part of the wave speed, and the average axial velocity ̅ is 0.228 m/s. After the wave passes (Figure 24b Stage 5: The wave cycle will be repeated in the pipeline due to inertia, in this stage, the pressure increase wave generated at the valve propagates upstream. Figure 22 is the axial velocity distribution at S1 monitoring position and Figure 23 is the corresponding local pressure distribution. At( 10 = 0.148 s = + 4.41 × / ), the pressure head at S1 position is 24.85 m, a little bigger than that of the original steady state due to the effect of the first part of the wave speed, and the average axial velocity ̅ is 0.228 m/s. After the wave passes (Figure 24b, 11 = 0.164 s = + 4.98 × / ), the pressure head reaches the maximum, 51.89 m, and the average axial velocity ̅ = 0 m/s. The local pressure distribution in Figure 23a,b show that the pressure gradient at S1 monitoring position at both time are not big because both time nodes are not pressure wave passing time.  Stage 5: The wave cycle will be repeated in the pipeline due to inertia, in this stage, the pressure increase wave generated at the valve propagates upstream. Figure 22 is the axial velocity distribution at S 1 monitoring position and Figure 23 is the corresponding local pressure distribution. At (T 10 = 0.148 s = Tc + 4.41 × L/a), the pressure head at S 1 position is 24.85 m, a little bigger than that of the original steady state due to the effect of the first part of the wave speed, and the average axial velocity U is 0.228 m/s. After the wave passes (Figure 24b, T 11 = 0.164 s = Tc + 4.98 × L/a), the pressure head reaches the maximum, 51.89 m, and the average axial velocity U = 0 m/s. The local pressure distribution in Figure 23a,b show that the pressure gradient at S 1 monitoring position at both time are not big because both time nodes are not pressure wave passing time.
original steady state due to the effect of the first part of the wave speed, and the average axial velocity ̅ is 0.228 m/s. After the wave passes (Figure 24b, 11 = 0.164 s = + 4.98 × / ), the pressure head reaches the maximum, 51.89 m, and the average axial velocity ̅ = 0 m/s. The local pressure distribution in Figure 23a,b show that the pressure gradient at S1 monitoring position at both time are not big because both time nodes are not pressure wave passing time.  In the process of the change of the axial velocity distribution at S1 monitoring position, two different regions can be identified: a wall region with a fast flow change and a high gradient, and a core region related to the velocity distribution of the past time history. These regions behave differently, because they depend on different flow characteristics: the wall region depends on the viscosity, and the core region depends on the fluid inertia. During the propagation of pressure wave, the change of velocity distribution is mainly due to diffusion, and the distribution always spread from the wall region to the core region, which is also the main reason for the change of turbulent condition in the pipe. There is a local pressure difference when the pressure wave passes the monitoring position, but the pressure difference is very small and negligible compared to the local average pressure; at the time before and after the pressure wave passes, the pressure on the monitoring position is almost constant; thus, the 1D and quasi 2D model assumption that pressure on the cross-section during the pressure wave propagation is constant is still valid.
In this study, the Reynolds stress term − ′ ′ ̅̅̅̅̅̅ in Equation (8) always exists in the calculations.
One reason is that the initial states of the two calculated cases are all turbulent flows. Another reason is that during the transient process, the turbulence dissipation rate on the wall region (where there is large velocity gradient) is bigger than that on the core region; it can be seen from = 2 / that In the process of the change of the axial velocity distribution at S1 monitoring position, two different regions can be identified: a wall region with a fast flow change and a high gradient, and a core region related to the velocity distribution of the past time history. These regions behave differently, because they depend on different flow characteristics: the wall region depends on the viscosity, and the core region depends on the fluid inertia. During the propagation of pressure wave, the change of velocity distribution is mainly due to diffusion, and the distribution always spread from the wall region to the core region, which is also the main reason for the change of turbulent condition in the pipe. There is a local pressure difference when the pressure wave passes the monitoring position, but the pressure difference is very small and negligible compared to the local average pressure; at the time before and after the pressure wave passes, the pressure on the monitoring position is almost constant; thus, the 1D and quasi 2D model assumption that pressure on the cross-section during the pressure wave propagation is constant is still valid.
In this study, the Reynolds stress term − ′ ′ ̅̅̅̅̅̅ in Equation (8) always exists in the calculations.
One reason is that the initial states of the two calculated cases are all turbulent flows. Another reason is that during the transient process, the turbulence dissipation rate on the wall region (where there is large velocity gradient) is bigger than that on the core region; it can be seen from = 2 / that the turbulent viscosity is very small on the wall region during the transient process. Therefore, to some extent, the used turbulence model can capture the flow characteristics both when the local flow changes from turbulent to laminar and when the local flow changes from laminar to turbulent. Figure  25 is the contours of turbulence dissipation rate on S1 monitoring position local region at the initial In the process of the change of the axial velocity distribution at S 1 monitoring position, two different regions can be identified: a wall region with a fast flow change and a high gradient, and a core region related to the velocity distribution of the past time history. These regions behave differently, because they depend on different flow characteristics: the wall region depends on the viscosity, and the core region depends on the fluid inertia. During the propagation of pressure wave, the change of velocity distribution is mainly due to diffusion, and the distribution always spread from the wall region to the core region, which is also the main reason for the change of turbulent condition in the pipe. There is a local pressure difference when the pressure wave passes the monitoring position, but the pressure difference is very small and negligible compared to the local average pressure; at the time before and after the pressure wave passes, the pressure on the monitoring position is almost constant; thus, the 1D and quasi 2D model assumption that pressure on the cross-section during the pressure wave propagation is constant is still valid.
In this study, the Reynolds stress term −ρu i u j in Equation (8) always exists in the calculations. One reason is that the initial states of the two calculated cases are all turbulent flows. Another reason is that during the transient process, the turbulence dissipation rate on the wall region (where there is large velocity gradient) is bigger than that on the core region; it can be seen from µ t = ρC µ k 2 /ε that the turbulent viscosity is very small on the wall region during the transient process. Therefore, to some extent, the used turbulence model can capture the flow characteristics both when the local flow changes from turbulent to laminar and when the local flow changes from laminar to turbulent. Figure 25 is the contours of turbulence dissipation rate on S 1 monitoring position local region at the initial time of T = 0 s and at T = 0.082 s, and it shows that during the transient process, the turbulence dissipation rate has a relatively big value on the wall region than on the core region.

Analysis of the Pressure Fluctuation
The numerical simulation was carried out for the transient cavitation flow experiment 2 (initial flow velocity of 0.332 m/s). The predicted cross-sectional average pressure values of all time steps at Sl, S2, and S3 monitoring positions are compared with the experimental data in Figure 26. The wave speed in the first and second stage of the first pressure wave cycle is constant ( = 1280 m/s) due to the material in the pipeline is single phase water, while the wave speed in the third and fourth stage of the first pressure wave cycle is variable due to cavitation. The average pressure wave speed 1 in the third and fourth stage can be evaluated by the pressure head change at Sl monitoring position and with equation as follows: where = 36 m is the pipe length and is the arrival time of the pressure decrease wave at the S1 monitoring position.
The pressure increase wave travels the length of 2 in the wave speed of = 1280 m/s during the first two stages (no cavitation) of the first wave cycle, and then the pressure decrease wave travels the length of 0.75 from the valve to the S1 monitoring position in the wave speed of 1 during the third stage (with cavitation) of the first wave cycle. The arrival time of the pressure decrease wave at S1 monitoring position, , is identified with the rarefaction wave contours, as shown in Figures 27 and 28. If = 0.0774 s is identified as the rarefaction wave arrival time, the result is 1 = 1277 m/s, while if = 0.08 s is identified as the rarefaction wave arrival time, the

Analysis of the Pressure Fluctuation
The numerical simulation was carried out for the transient cavitation flow experiment 2 (initial flow velocity of 0.332 m/s). The predicted cross-sectional average pressure values of all time steps at S l , S 2 , and S 3 monitoring positions are compared with the experimental data in Figure 26. The wave speed a in the first and second stage of the first pressure wave cycle is constant (a = 1280 m/s) due to the material in the pipeline is single phase water, while the wave speed in the third and fourth stage of the first pressure wave cycle is variable due to cavitation. The average pressure wave speed a 1 in the third and fourth stage can be evaluated by the pressure head change at S l monitoring position and with equation as follows: where L = 36 m is the pipe length and Tv is the arrival time of the pressure decrease wave at the S 1 monitoring position. The pressure increase wave travels the length of 2L in the wave speed of a = 1280 m/s during the first two stages (no cavitation) of the first wave cycle, and then the pressure decrease wave travels the length of 0.75L from the valve to the S 1 monitoring position in the wave speed of a 1 during the third stage (with cavitation) of the first wave cycle. The arrival time of the pressure decrease wave at S 1 monitoring position, Tv, is identified with the rarefaction wave contours, as shown in Figures 27  and 28. If Tv = 0.0774 s is identified as the rarefaction wave arrival time, the result is a 1 = 1277 m/s, while if Tv = 0.08 s is identified as the rarefaction wave arrival time, the result is a 1 = 1136.84 m/s. Thus, there are different values of for the rarefaction wave average speed, and because the pressure gradient is bigger when T = 0.08 s than when T = 0.0774 s, T = 0.08 s is chosen as the rarefaction wave arrival time, and the value of a 1 = 1136.84 m/s is chosen as the wave speed of the rarefaction wave and it is used in the analysis of the pressure and velocity field in Section 3.4.2, which indicates the pressure wave speed in the cavitation flow is smaller than in the non-cavitation flow.
The simulation results demonstrate that the transient damping is increased in the transient cavitation flow [5]. The vapor cavity formation and duration, the instantaneous pressure peak caused by the cavity collapse, and the phase can match well with the experimental results. The simulation results demonstrate that the transient damping is increased in the transient cavitation flow [5]. The vapor cavity formation and duration, the instantaneous pressure peak caused by the cavity collapse, and the phase can match well with the experimental results.    The instance when maximum vapor volume fraction forms correspond to the moment of the minimum pressure head. The vapor volume fraction peak value reduces faster than the pressure head, and when the second cavitation occurs at T = 0.02359 s, the vapor volume fraction value increases again, but never reach the maximum value in the first vapor volume fraction peak. Later on, the vapor volume fraction diminishes, although there is still a pressure head change in the pipeline.   The instance when maximum vapor volume fraction forms correspond to the moment of the minimum pressure head. The vapor volume fraction peak value reduces faster than the pressure head, and when the second cavitation occurs at T = 0.02359 s, the vapor volume fraction value increases again, but never reach the maximum value in the first vapor volume fraction peak. Later on, the vapor volume fraction diminishes, although there is still a pressure head change in the pipeline.  The instance when maximum vapor volume fraction forms correspond to the moment of the minimum pressure head. The vapor volume fraction peak value reduces faster than the pressure head, and when the second cavitation occurs at T = 0.02359 s, the vapor volume fraction value increases again, but never reach the maximum value in the first vapor volume fraction peak. Later on, the vapor volume fraction diminishes, although there is still a pressure head change in the pipeline.   Figure 26c,d are the experimental data and numerical simulation results of pressure head fluctuation at S 3 monitoring position, with marked time nodes T 0 ∼ T 6 of several typical moments in the first cavity formation-development-crash evolution process. The vapor phase volume fraction of these typical moments at the valve upstream local region on the axial plane of the pipeline and the corresponding pressure and velocity fields are analyzed in two stages: the cavity region growth stage (Figures 30-32) and the cavity region reduction stage (Figures 33-35).

Analysis of the Vapor Volume Fraction and Pressure and Velocity Fields
Water 2018, 10, x FOR PEER REVIEW 23 of 29 during this time, the pressure decrease wave generated at the valve in the third stage of the pressure wave cycle propagates from the valve to the upstream reservoir, as can be seen in Figure 31, and velocity field in Figure 32 shows that the cross-section average velocity in the pipe gradually decreases to 0.  during this time, the pressure decrease wave generated at the valve in the third stage of the pressure wave cycle propagates from the valve to the upstream reservoir, as can be seen in Figure 31, and velocity field in Figure 32 shows that the cross-section average velocity in the pipe gradually decreases to 0.  The cavity region growth stage: at T 0 = 0.0757 s = Tc + 2L/a + 0.11L/a 1 , the pressure at the S 3 monitoring point is firstly detected to fall below the vapor pressure. As can be seen in Figure 33a, a local cavity is formed at the top of the pipe wall (the vapor phase volume fraction value is small), and the vapor zone develops upstream along the top of the pipe wall, which is different from the traditional DVCM and DGCM model assumption that "the vapor cavity occupies the entire calculated section and the position remains the same"; in the subsequent time, due to the pressure decrease wave propagates from the valve to the upstream reservoir, the local cavity volume fraction at the S 3 monitoring point gradually increases and the vapor region continues to develop along the top of the pipe wall toward the upstream reservoir (T 1 = 0.0761 s = Tc + 2L/a + 0.12L/a 1 ). The maximum length of the cavity region along the pipe is 3.73 m, and it occurs at T 2 = 0.078s = Tc + 2L/a + 0.18L/a 1 . The cavity region growth stage is from T 0 = 0.0757 s to T 2 = 0.078 s, it lasts 0.0023 s, during this time, the pressure decrease wave generated at the valve in the third stage of the pressure wave cycle propagates from the valve to the upstream reservoir, as can be seen in Figure 31, and velocity field in Figure 32 shows that the cross-section average velocity in the pipe gradually decreases to 0.  The cavity region reduction stage: although the pressure decrease wave continues to propagate upstream at 3 = 0.087 s = + 2 / + 0.47 / 1 and 4 = 0.097 s = + 2 / + 0.78 / 1 , the vapor region stops to grow and it begins to reduce, but the vapor phase fraction begins to grow, as can be seen in Figures 33a,b and 34a,b, and the corresponding velocity fields in Figure 35a,b shows that the cross-sectional average velocity in the pipeline gradually decreases to 0. When the pressure increase wave generated at the upstream reservoir propagates to the valve (Figure 34c    The total duration time of the cavity is 0.0597 s. When the rarefaction wave propagates into a zone of decreasing pressure, a vaporous cavitation region is formed; alternatively, when the rarefaction wave propagates into a zone of increasing pressure, no vaporous cavitation region is formed. The hydrostatic pressure difference from top to bottom of the pipe is small and negligible compared to the pressure peak value (1.5%) because the pipe is long and slender, the flow characteristics are mainly changed by the pressure wave propagation; while when cavitation occurs, the gravity effect is important to make the cavity region stay on the top of the pipe. Because the vapor density is much smaller than water density, and buoyancy is caused by the density gradient, the vapor formed at the top of the pipe is inclined to stay upward.

Analysis of the Variable Sound Speed Field in the Transient Cavitation Flow
where is the vapor bulk modulus and = 0 , is the specific heat ratio with = 1.4, and 0 is the reference pressure with 0 = 101,325 Pa.
In the transient cavitation flow, pressure wave speed is variable, and pressure wave in the cavitation region is smaller than in the non-cavitation region. Pressure wave field can match with the vapor volume fraction field at all these time nodes. The weakly compressible fluid calculation model in this paper can capture the variable wave speed in the transient cavitation flow.
As for the pressure wave speed contours in Figure 36, at T2 = 0.078 s, T3 = 0.087 s and T4 = 0.097 s, the average pressure wave speed values are all 1278 m/s, which is almost equal to the pressure The cavity region reduction stage: although the pressure decrease wave continues to propagate upstream at T 3 = 0.087 s = Tc + 2L/a + 0.47L/a 1 and T 4 = 0.097 s = Tc + 2L/a + 0.78L/a 1 , the vapor region stops to grow and it begins to reduce, but the vapor phase fraction begins to grow, as can be seen in Figure 33a,b and Figure 34a,b, and the corresponding velocity fields in Figure 35a,b shows that the cross-sectional average velocity in the pipeline gradually decreases to 0. When the pressure increase wave generated at the upstream reservoir propagates to the valve (Figure 34c, T 5 = 0.1333 s = Tc + 2L/a + 1.93L/a 1 and Figure 34d T 6 = 0.1354 s = Tc + 2L/a + 1.99L/a 1 ), the cavity is decompressed and the local cavity region reduces and finally disappears (Figure 33c,d); during this time, the velocity fields in Figure 35c,d show the velocity in the pipeline returns to the flow direction from the reservoir to the valve. The cavity region reduction stage is from T 3 = 0.087 s to T 6 = 0.1354 s, and it lasts for 0.0484 s.
The total duration time of the cavity is 0.0597 s. When the rarefaction wave propagates into a zone of decreasing pressure, a vaporous cavitation region is formed; alternatively, when the rarefaction wave propagates into a zone of increasing pressure, no vaporous cavitation region is formed. The hydrostatic pressure difference from top to bottom of the pipe is small and negligible compared to the pressure peak value (1.5%) because the pipe is long and slender, the flow characteristics are mainly changed by the pressure wave propagation; while when cavitation occurs, the gravity effect is important to make the cavity region stay on the top of the pipe. Because the vapor density is much smaller than water density, and buoyancy is caused by the density gradient, the vapor formed at the top of the pipe is inclined to stay upward. Figure 36 is the pressure wave speed field evolution in the transient cavitation flow. The pressure wave speed is calculated according to Equation (16) [33].

Analysis of the Variable Sound Speed Field in the Transient Cavitation Flow
where K v is the vapor bulk modulus and K v = γP 0 , γ is the specific heat ratio with γ = 1.4, and P 0 is the reference pressure with P 0 = 101,325 Pa.

Conclusions
Based on the weakly compressible fluid RANS method, it can make up for the shortcomings of the existing 1D MOC, capture the detailed flow characteristics of the complex pipeline hydraulic system, and even the related cavitation characteristics, and explore the dynamic characteristics of the cavity when transient cavitation occurs in the pipeline water delivery system and the inner rules.
In this paper, by considering the compressibility of water and introducing the fluid densitypressure equation that represented the pressure wave speed, a transient cavitation flow calculation model based on the weakly compressible fluid RANS method was proposed. The grid slip technique was used in the simulation of transient flow in reservoir-pipe-valve system, and the transient non- In the transient cavitation flow, pressure wave speed is variable, and pressure wave in the cavitation region is smaller than in the non-cavitation region. Pressure wave field can match with the vapor volume fraction field at all these time nodes. The weakly compressible fluid calculation model in this paper can capture the variable wave speed in the transient cavitation flow.