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

Abstract

In hydraulic systems, transient flow 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 flow in pipeline; the density-pressure model was incorporated into the continuity equation by using further development of UDF (user defined function), which reflects the variable wave speed of the transient cavitation flow, and the related algorithms were established based on weakly compressible fluid Reynolds Average Navier-Stokes (RANS) techniques. Firstly, the numerical simulations of the transient non-cavitation and cavitation flows 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 flow field information such as velocity, pressure and vapor volume fraction. Through the evolution process of the pressure field, the propagation characteristics of pressure waves can be analyzed qualitatively and quantitatively. Through the evolution process of the velocity field, 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 diffusion. 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 different results of cavitation flow. The accuracy and reliability of the weakly compressible fluid RANS method were verified by comparing the calculated results with the experimental data.

1. 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 extreme pressures typically observed in the first cycles. Many mathematical models for 1D transient cavitation flows have been proposed, including Discrete Vapor Cavity Model (DVCM) [2], Discrete Gas Cavity Model (DGCM) [3], and the relevant modification models [4]. In these models, the transient cavitation flow phenomenon is replayed, and some of the simulation results are basically consistent with the experimental data, but they assume that the vapor cavity occupies the entire calculated section and the position remains the same, which is not the real condition, and they are impossible to reflect the essence of the cavity volume fraction transport. Zhu et al. [5] used two different methods of DVCM and DGCM to calculate the wave speed in homogeneous air-water mixing flows in viscoelastic pipes and the results demonstrate that the existence of air content in transient water flows may reduce greatly the influence of viscoelasticity effect and meanwhile increase relatively the importance of unsteady friction effect on transient damping. As efficient techniques to get more information of hydraulic transients, the more accurate quasi two-dimensional (2D) unsteady friction model, which takes the velocity distribution of the cross-section into consideration attract more and more attention [6,7]. However, the quasi-2D model is based on the assumption that the flow velocity is axisymmetric, the radial momentum equation is neglected, the radial velocity component and its derivatives in the axis momentum equation are neglected, and the pressure in the cross-section is assumed to be constant, which are not the real conditions.

With rapid developments in computer capabilities and Computational Fluid Dynamics (CFD), CFD is a mature technology and reliable engineering tool for analyses of hydraulic engineering problems, including hydraulic transient in pipes [8]. Martins et al. [9] used CFD to analyze the hydraulic transient flows in pressurized pipes, and the calculated velocity profiles demonstrated two regions with different behaviors: the wall region dominated by the fluid viscosity showing flow changes are faster with high gradients and the pipe core region strongly dependent on the fluid inertial forces; but the downstream boundary condition used was the flowrate-time function with the flowrate change from the transient flow experiment, which cannot represent the ball valve motion effect, and the flowrate change time step was not verified; moreover, the transient cavitation flow was not modeled. Wang et al. [10] combined 1D method of characteristics and the three-dimensional (3D) finite volume method of CFD based on OpenFOAM open source CFD software to carry out transient simulations of hydraulic system caused by the valve closing/opening in a pipe, but the data transfer between the 1D Method of Characteristics (MOC) part and the 3D CFD part may have errors and the accuracy of the data exchange scheme was not mathematically verified, besides, the transient cavitation flow was not modeled, and the wave speed was constant. Wu et al. [11] also used 1D MOC and 3D CFD to study the interaction between valve-induced water hammer and pump during the rapid closure of the valve and demonstrate that the H-Q curve in transient operation evidently deviate from the steady-state value.

Yang et al. [12] calculated the valve closure induced water hammer using the 3D CFD method and showed advantage over MOC, but the mesh size on the wall region of the pipe was not fine enough and the flow characteristics at different region of the pipe during the transient process was not analyzed. Saemi et al. [13] used 2D and 3D CFD method to simulate water hammer flows with initial laminar and turbulent steady state, and the results demonstrate that at the instance of wave passage, the inertia and pressure gradient contributions are considerable away from the wall, while in the near-wall region the effects of viscosity, pressure gradient, and inertia are important. Mandair et al. [14] compared the MOC method and 3D CFD method in calculation of the water hammer in Reservoir-Pipeline-Valve system and demonstrated that the CFD method can provide more flow information than the MOC method. Su et al. [15] used the 3D CFD method to calculate the water hammer pressure in a high-pressure common rail system and revealed a reverse flow, secondary flow, and flow separation region during the water hammer process. Heng et al. [16] used the CFD method to calculate the flow field of the stop valve in the closing process, and the water hammer pressure affected by the valve closing method was analyzed. Wang et al. [17] used CFD simulation to describe the water hammer in a pipeline with a vortex diode and found that a vortex diode used as a leaky check-valve in a pipeline system acts significantly in suppressing pressure fluctuation of water hammer. In these papers, although water 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.

2. Methodology

2.1. 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:

$$\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u_i)}{\partial x_i}=0$$

(1)

where $\rho$ is the density of the fluid and $u_i$ is the i-direction component of velocity of the fluid.

For compressible fluid considering the transient flow, there is:

$$a^2=\frac{\partial p}{\partial \rho }$$

(2)

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:

$$\frac{\partial p}{\partial t}+\rho a^2\frac{\partial u_i}{\partial x_i}+u_i\frac{\partial p}{\partial x_i}=0$$

(3)

Let ${\rho }^{*}=\frac{\rho }{{\rho }_0}$, $t^{*}=\frac{t}{t_0}$, $\frac{\partial }{\partial x_i}=l_0\frac{\partial }{\partial x_i^{*}}$, ${\overrightarrow{u}}^{*}=\frac{\overrightarrow{u}}{u_0}$, $p^{*}=\frac{p-p_0}{{\rho }_0{u}_0^2}$, $a^{*}=\frac{a}{a_0}$, then Equation (1) is converted to dimensionless as:

$$(St){(Ma)}^2\left(\frac{\partial p^{*}}{\partial t^{*}}\right)+{\rho }^{*}{a}^{*}{}^2\frac{\partial u_i^{*}}{\partial x_i^{*}}+{(Ma)}^2{u}_i^{*}\frac{\partial p^{*}}{\partial x_i^{*}}=0$$

(4)

where Strauhal number $St=\frac{l_0}{u_0{t}_0}$, Mach number $Ma=\frac{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$, $M{a}^2\ll 1$, $St\gg 1$, thus, the last term of Equation (4) can be omitted. Therefore, the continuous equation of weakly compressible fluid is:

$$\frac{\partial p}{\partial t}+\rho a^2\frac{\partial u_i}{\partial x_i}=0$$

(5)

The momentum equation of weakly compressible fluid is:

$$\frac{\partial }{\partial t}(\rho u_i)+\frac{\partial }{\partial x_j}(\rho u_j{u}_i)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right]-\frac{2}{3}\frac{\partial }{\partial x_i}\left(\mu \frac{\partial u_k}{x_k}\right)+\rho f_i$$

(6)

where $\mu$ 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.

2.2. Solution of RANS Equation Based on the Weakly Compressible Model

RANS equations are used to describe the governing equations of weakly compressible fluids:

$$\frac{\partial \overline{p}}{\partial t}+\rho a^2\frac{\partial {\overline{u}}_i}{\partial x_i}=0$$

(7)

$$\frac{\partial }{\partial t}(\rho {\overline{u}}_i)+\frac{\partial }{\partial x_j}(\rho {\overline{u}}_j{\overline{u}}_i)=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}(\mu \frac{\partial {\overline{u}}_i}{\partial x_j}+{\tau }_{ij})+\rho f_i$$

(8)

where the over bar indicates the time average components and ${\tau }_{ij}$ represents Reynolds stress and ${\tau }_{ij}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$.

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-\epsilon$ model, RNG $k-\epsilon$ model, Realizable $k-\epsilon$ model, SST $k-\omega$ model etc. The standard two-equation eddy viscosity models ($k-\epsilon$ and $k-\epsilon$ types of model) based on the Boussinesq relationship ${\tau }_{ij}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}=-2/3(\rho k+{\mu }_t\partial {\overline{u}}_i/\partial x_i){\delta }_{ij}+{\mu }_t(\partial {\overline{u}}_i/\partial x_j+\partial {\overline{u}}_j/\partial x_i)$ (where $k$ is the turbulence kinetic energy, ${\delta }_{ij}$ is the Kronecker delta, and ${\mu }_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-\epsilon$ model that can better simulate curvature is used, and the transport equations of the turbulence kinetic energy, $k$, and the turbulence dissipation rate, $\epsilon$, are as follows:

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

(9)

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

(10)

where ${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$ and $G_k$ is the generation of turbulence kinetic energy due to the mean velocity gradients, computed by $G_k={\mu }_t(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})\frac{\partial u_i}{\partial x_j}$; ${\alpha }_k$ and ${\alpha }_{\epsilon }$ are the turbulent Prandtl number for $k$ and $\epsilon$, and ${\alpha }_k={\alpha }_{\epsilon }=1.39$. The other model constants are, respectively: $C_{1\epsilon }{}^{*}=C_{1\epsilon }-\frac{\eta (1-\eta /{\eta }_0)}{1+\beta {\eta }^3}$; $C_{1\epsilon }=1.42$; $C_{2\epsilon }=1.68$; $C_{\mu }=0.0845$; ${\eta }_0=4.377$; $\beta =0.012$; $\eta ={(2E_{ij}{E}_{ij})}^{1/2}k/\epsilon$; $E_{ij}=0.5\cdot (\partial u_i/\partial x_j+\partial u_j/\partial x_i)$.

2.3. 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:

$$\frac{\partial }{\partial t}(\alpha {\rho }_v)+\frac{\partial }{\partial x_i}(\alpha {\rho }_v{u}_{vi})=R_e-R_c$$

(11)

where subscript $v$ represents the vapor phase, $\alpha$ is the vapor volume fraction, ${\rho }_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:

$$\mathrm{When}{p}_v\ge p:R_e=\frac{{\rho }_l{\rho }_v}{{\rho }_m}\alpha (1-\alpha )\frac{3}{R_B}\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_l}}$$

(12)

$$\mathrm{When}{p}_v<p:R_c=\frac{{\rho }_l{\rho }_v}{{\rho }_m}\alpha (1-\alpha )\frac{3}{R_B}\sqrt{\frac{2}{3}\frac{p-p_v}{{\rho }_l}}$$

(13)

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={(\frac{\alpha }{1-\alpha }\frac{3}{4\pi }\frac{1}{n})}^{\frac{1}{3}}$, in which $n$ is the bubble number density.

2.4. 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 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. Main parameters of different experimental schemes.

Experiment No.	Initial Velocity (m/s)	Initial Reynolds Number	Water Level of Upstream Reservoir (m)
Experiment 1	0.239	4531	24.30
Experiment 2	0.332	6294	23.41

. 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₂ 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₃ 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].

2.5. 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. Boundary conditions of the calculation cases.

Case No.	Initial Velocity (m/s)	Pressure of Inlet (Pa)	Pressure of Outlet (Pa)	Valve Closing Time (ms)
Case 1	0.239	237,739.7	235,960.4	24
Case 2	0.332	229,060.1	225,776.9	16

. 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:

$${\rho }_l={\rho }_{l0}(1+\frac{\Delta p}{K_l})$$

(14)

where ${\rho }_{l0}$ and $K_l$ are the density and bulk modulus of the liquid at absolute reference pressure $p_0$, respectively, and $\Delta p=p-p_0$.

3. Results and Discussion

3.1. 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₃ 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.

3.2. 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₃ 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.

3.3. Analysis of Transient Non-Cavitation Flow Results

3.3.1. 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₂, and S₃ 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 S₃ 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.

3.3.2. Analysis of the Pressure and Velocity Fields

Figure 7 and Figure 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\sim 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\sim 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. Time nodes of $T_1\sim T_{12}$ in terms of $Tc$ and $L/a$.

Time Node	Time(s)	In Terms of $\mathit{T}\mathit{c}$ and $\mathit{L}\mathbf{/}\mathit{a}$
$T_1$	0.024	$Tc+0\times L/a$
$T_2$	0.036	$Tc+0\times L/a$
$T_3$	0.052	$Tc+1\times L/a$
$T_4$	0.06	$Tc+1.28\times L/a$
$T_5$	0.072	$Tc+1.71\times L/a$
$T_6$	0.08	$Tc+1.99\times L/a$
$T_7$	0.086	$Tc+2.2\times L/a$
$T_8$	0.098	$Tc+2.63\times L/a$
$T_9$	0.108	$Tc+2.99\times L/a$
$T_{10}$	0.12	$Tc+3.41\times L/a$
$T_{11}$	0.13	$Tc+3.77\times L/a$
$T_{12}$	0.136	$Tc+3.98\times L/a$

.

The pictures in Figure 7 and Figure 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=\sqrt{K_l/{\rho }_{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.

3.3.3. Analysis of the Velocity Distribution Change

Figure 9 is the experimental data and numerical simulation results of pressure head fluctuation at S₁ monitoring point, with marked time nodes $T_0\sim T_{11}$ (different from time nodes in Figure 6b). The velocity distribution at S₁ 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 according to the pressure wave propagation directions. The change of average axial velocity $\overline{U}$ and pressure head $H$ at S₁ position are recorded and analyzed.

Stage 1: Figure 10 and Figure 11 are respectively the axial velocity distribution and local pressure distribution at the S₁ 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₁ monitoring position exhibits a fully developed steady-state velocity distribution, with $\overline{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₁ 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\mathrm{s}=Tc+0.71\times L/a$, $H$ = 37.29 m), the flow velocity in the pipeline gradually decreases (Figure 10b, $\overline{U}$ = 0.137 m/s); the local pressure distribution at S₁ 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₁ 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 (Figure 10c and Figure 12c, $T_2=0.05\mathrm{s}=Tc+0.92\times L/a$), the pressure at S₁ 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₁ monitoring position is almost constant at this time (Figure 11c).

Stage 2: Figure 13 and Figure 14 are respectively the axial velocity distribution and local pressure distribution at the S₁ 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\mathrm{s}=Tc+1.21\times L/a$), the pressure head at S₁ position decreases to 41.29 m, the velocity direction is from the valve to the reservoir (Figure 13a, $\overline{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₁ 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\mathrm{s}=Tc+1.42\times L/a$), the pressure head decreases to 24.46 m, the average axial velocity $\overline{U}$ = −0.236 m/s, and the pressure at S₁ 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₁ 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, Figure 16 and Figure 17 are respectively the axial velocity distribution and local pressure distribution at the S₁ monitoring position of different time, and Figure 18 is the pressure field evolution on the axis plane of the pipeline in this stage. At $T_5=0.092\mathrm{s}=Tc+2.42\times L/a$, the pressure head at S₁ 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 $\overline{U}$ = −0.233 m/s), and the velocity change gradually diffuses to the core region, at $T_6=0.1\mathrm{s}=Tc+2.7\times L/a$, $\overline{U}$ = −0.14 m/s. After the pressure wave passes (Figure 18c, $T_7=0.108\mathrm{s}=Tc+2.99\times L/a$), the pressure head at S₁ 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₁ 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₁ 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: Figure 19 and Figure 20 are respectively the axial velocity distribution and local pressure distribution at the S₁ 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 S₁ monitoring position (Figure 21a, $T_8=0.116\mathrm{s}=Tc+3.27\times L/a$), the pressure head at S₁ position is 12.91 m, the velocity at the wall region is positive Figure 19a, with average axial velocity $\overline{U}$ = 0.148 m/s; the local pressure difference is 14.5 Pa (Figure 20a), the average pressure on S₁ 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, $T_9=0.124\mathrm{s}=Tc+3.56\times L/a$), the pressure head is 24.15 m, $\overline{U}$ = 0.233 m/s, and the local pressure on S₁ 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 S₁ monitoring position and Figure 23 is the corresponding local pressure distribution. At ($T_{10}=0.148\mathrm{s}=Tc+4.41\times L/a$), the pressure head at S₁ 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 $\overline{U}$ is 0.228 m/s. After the wave passes (Figure 24b, $T_{11}=0.164\mathrm{s}=Tc+4.98\times L/a$), the pressure head reaches the maximum, 51.89 m, and the average axial velocity $\overline{U}$ = 0 m/s. The local pressure distribution in Figure 23a,b show that the pressure gradient at S₁ 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 S₁ 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 $-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$ 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 ${\mu }_t=\rho C_{\mu }{k}^2/\epsilon$ 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₁ 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.

3.4. Analysis of Transient Cavitation Flow Results

3.4.1. 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₂, and S₃ 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:

$$\frac{2L}{a}+\frac{\frac{3}{4}L}{a_1}=Tv$$

(15)

where $L$ = 36 m is the pipe length and $Tv$ is the arrival time of the pressure decrease wave at the S₁ 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₁ 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₁ monitoring position, $Tv$, is identified with the rarefaction wave contours, as shown in Figure 27 and Figure 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.

Figure 29 is the pressure head change and vapor volume fraction change at S₃ monitoring position (0.02 m upstream of the valve) during the transient event. 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.

3.4.2. Analysis of the Vapor Volume Fraction and Pressure and Velocity Fields

Figure 26c,d are the experimental data and numerical simulation results of pressure head fluctuation at S₃ monitoring position, with marked time nodes $T_0\sim 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 (Figure 30, Figure 31 and Figure 32) and the cavity region reduction stage (Figure 33, Figure 34 and Figure 35).

The cavity region growth stage: at $T_0=0.0757\mathrm{s}=Tc+2L/a+0.11L/a_1$, the pressure at the S₃ 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₃ 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\mathrm{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.0757 s to T₂ = 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 $T_3=0.087\mathrm{s}=Tc+2L/a+0.47L/a_1$ and $T_4=0.097\mathrm{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\mathrm{s}=Tc+2L/a+1.93L/a_1$ and Figure 34d $T_6=0.1354\mathrm{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₃ = 0.087 s to T₆ = 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.

3.4.3. Analysis of the Variable Sound Speed Field in the Transient Cavitation Flow

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].

$$a=\sqrt{\frac{K_l}{[{\rho }_l-({\rho }_l-{\rho }_v){\alpha }_v][1+(\frac{K_l}{K_v}-1){\alpha }_v]}}$$

(16)

where $K_v$ is the vapor bulk modulus and $K_v=\gamma P_0$, $\gamma$ is the specific heat ratio with $\gamma$ = 1.4, and $P_0$ is the reference pressure with $P_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 T₂ = 0.078 s, T₃ = 0.087 s and T₄ = 0.097 s, the average pressure wave speed values are all 1278 m/s, which is almost equal to the pressure wave speed value of 1277 m/s, when the rarefaction wave arrives at the S₁ monitoring position at $Tv$ = 0.0774 s (illustrated in Section 3.4.1). However, the average wave speed from Figure 36 at T₀ = 0.0757 s, T₁ = 0.0761 s and T₆ = 0.1354 s are 1280 m/s because of the small vapor volume fraction at these times. Therefore, the pressure wave speed in Figure 36 calculated with Equation (16) are almost identical to the values calculated with the arrival time of the rarefaction wave.

4. 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 density-pressure 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-cavitation flow and transient cavitation flow caused by fast closure of the valve were numerically simulated. The main conclusions are:

(1) The weakly compressible fluid RANS method transient cavitation calculation model predicted a high degree of agreement between the calculated pressure head and the experimental data. Moreover, the accuracy and reliability of the calculation model of the transient cavitation flow of the weakly compressible fluid RANS method and related boundary conditions were verified.

(2) The grid scale has little effect on the calculation results of the pressure head. The time step has little effect on the maximum water hammer pressure and pressure fluctuation cycle, but it has a greater influence on the attenuation of the pressure head: the smaller the time step, the slower the decay; the larger the time step, the faster the decay.

(3) Based on the weakly compressible fluid RANS method, it can make up for the deficiency of the existing 1D MOC, and the propagation characteristics of the pressure wave can be analyzed qualitatively and quantitatively through the evolution process of the pressure field. The evolution process of the velocity field showed that the velocity distribution in the wall region changed rapidly and had a high gradient, and the change depended mainly on the viscosity. The change of the velocity distribution in the core region was related to the velocity distribution in the past time history, mainly depending on the diffusion. When the pressure wave passes, the local pressure difference on the monitoring cross-section is small and negligible compared with the local average pressure.

(4) Based on the weak compressible fluid RANS method, the development process and distribution of the vapor cavity can be clearly and intuitively observed. The vapor cavity was unevenly distributed in the pipe axis and pipe diameter direction and moved slowly along the top of the pipe wall, and the entire process of cavity formation-development-cracking was captured. The non-symmetry vapor distribution during the transient process demonstrate the advantage of the used 2D CFD method over the quasi-2D MOC based method, which assume that the pressure is constant on the cross-section and cannot capture the non-symmetry cavity distribution. When the rarefaction wave propagated into a zone of decreasing pressure, a vaporous cavitation region was formed; alternatively, when the rarefaction wave propagated into a zone of increasing pressure, no vaporous cavitation region was formed.

(5) The pressure wave speed in the cavitation region is smaller than in the non-cavitation region, the pressure wave speed field can match with the vapor volume fraction field, and the presented calculation model can capture the variable wave speed in the transient cavitation flow. With the transient cavitation flow pressure wave calculation method, the variable wave speed essence during the transient cavitation flows is revealed, and the results can match with the pressure wave speed calculated from the rarefaction wave travel time.