Next Article in Journal
Diatom Species Richness in Swiss Springs Increases with Habitat Complexity and Elevation
Next Article in Special Issue
Effects of Dam Regulation on the Hydrological Alteration and Morphological Evolution of the Volta River Delta
Previous Article in Journal
Slope Stability of a Scree Slope Based on Integrated Characterisation and Monitoring
Previous Article in Special Issue
Numerical Prediction of the Short-Term Trajectory of Microplastic Particles in Laizhou Bay
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
College of Water Resources and Civil Engineering, China Agricultural University, Beijing 100083, China
2
Beijing Engineering Research Centre of Safety and Energy Saving Technology for Water Supply Network System, China Agricultural University, Beijing 100083, China
*
Author to whom correspondence should be addressed.
Water 2020, 12(2), 448; https://doi.org/10.3390/w12020448
Submission received: 24 November 2019 / Revised: 10 January 2020 / Accepted: 31 January 2020 / Published: 7 February 2020

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:
ρ t + ( ρ u i ) x i = 0
where ρ 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 = p ρ
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:
p t + ρ a 2 u i x i + u i p x i = 0
Let ρ * = ρ ρ 0 , t * = t t 0 , x i = l 0 x i * , u * = u u 0 , p * = p p 0 ρ 0 u 0 2 , a * = a a 0 , then Equation (1) is converted to dimensionless as:
( S t ) ( M a ) 2 ( p * t * ) + ρ * a * 2 u i * x i * + ( M a ) 2 u i * p * x i * = 0
where Strauhal number S t = l 0 u 0 t 0 , Mach number M a = 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, M a < 1 , M a 2 1 , S t 1 , thus, the last term of Equation (4) can be omitted. Therefore, the continuous equation of weakly compressible fluid is:
p t + ρ a 2 u i x i = 0
The momentum equation of weakly compressible fluid is:
t ( ρ u i ) + x j ( ρ u j u i ) = p x i + x j [ μ ( u i x j + u j x i ) ] 2 3 x i ( μ u k x k ) + ρ f i
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.

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:
p ¯ t + ρ a 2 u ¯ i x i = 0
t ( ρ u ¯ i ) + x j ( ρ u ¯ j u ¯ i ) = p ¯ x i + x j ( μ u ¯ i x j + τ i j ) + ρ f i
where the over bar indicates the time average components and τ i j represents Reynolds stress and τ i j = ρ 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 τ i j = ρ u i u j ¯ = 2 / 3 ( ρ k + μ t u ¯ i / x i ) δ i j + μ t ( u ¯ i / x j + u ¯ j / x i ) (where k is the turbulence kinetic energy, δ i j 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:
( ρ k ) t + ( ρ k u i ) x i = x j [ α k ( μ + μ t ) k x j ] + G k ρ ε
( ρ ε ) t + ( ρ ε u i ) x i = x j [ α ε ( μ + μ t ) ε x j ] + C 1 ε * ε k G k C 2 ε ρ ε 2 k
where μ t = ρ C μ k 2 ε and G k is the generation of turbulence kinetic energy due to the mean velocity gradients, computed by G k = μ t ( u i x j + u j x i ) u i x j ; α k and α ε are the turbulent Prandtl number for k and ε , and α k = α ε = 1.39 . The other model constants are, respectively: C 1 ε * = C 1 ε η ( 1 η / η 0 ) 1 + β η 3 ; C 1 ε = 1.42 ; C 2 ε = 1.68 ; C μ = 0.0845 ; η 0 = 4.377 ; β = 0.012 ; η = ( 2 E i j E i j ) 1 / 2 k / ε ; E i j = 0.5 ( u i / x j + u j / 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:
t ( α ρ v ) + x i ( α ρ v u v i ) = R e R c
where subscript v represents the vapor phase, α is the vapor volume fraction, ρ v is the vapor density, u v i 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   p v p :   R e = ρ l ρ v ρ m α ( 1 α ) 3 R B 2 3 p v p ρ l
When   p v < p :   R c = ρ l ρ v ρ m α ( 1 α ) 3 R B 2 3 p p v ρ l
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 = ( α 1 α 3 4 π 1 n ) 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. 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.
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. 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:
ρ l = ρ l 0 ( 1 + Δ p K l )
where ρ l 0 and K l are the density and bulk modulus of the liquid at absolute reference pressure p 0 , respectively, and Δ 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 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.

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

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

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 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, T c , and the one-way travel time of the pressure wave in the pipe, L / a , as shown in Table 3.
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 = K l / ρ l 0 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 = 2 L / 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 = 2 L / 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 = 3 L / 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 = 3 L / 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 = 4 L / 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 S1 monitoring point, with marked time nodes T 0 T 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 according to the pressure wave propagation directions. The change of average axial velocity U ¯ and pressure head H at S1 position are recorded and analyzed.
Stage 1: Figure 10 and Figure 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 T 0 = 0 s time at S1 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 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, T 1 = 0.044   s = T c + 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 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 (Figure 10c and Figure 12c, T 2 = 0.05   s = T c + 0.92 × L / a ), the pressure at S1 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 S1 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 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, T 3 = 0.058   s = T c + 1.21 × L / a ), the pressure head at S1 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 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, T 4 = 0.064   s = T c + 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 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, Figure 16 and Figure 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 T 5 = 0.092   s = T c + 2.42 × L / a , 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 U ¯ = −0.233 m/s), and the velocity change gradually diffuses to the core region, at T 6 = 0.1   s = T c + 2.7 × L / a , U ¯ = −0.14 m/s. After the pressure wave passes (Figure 18c, T 7 = 0.108   s = T c + 2.99 × L / a ), 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 4: Figure 19 and Figure 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, T 8 = 0.116   s = T c + 3.27 × L / a ), the pressure head at S1 position is 12.91 m, the velocity at the wall region is positive Figure 19a, with average axial velocity U ¯ = 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, T 9 = 0.124   s = T c + 3.56 × L / a ), the pressure head is 24.15 m, U ¯ = 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 ( T 10 = 0.148   s = T c + 4.41 × L / a ), 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 U ¯ is 0.228 m/s. After the wave passes (Figure 24b, T 11 = 0.164   s = T c + 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 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 ρ 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 S1 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 Sl, S2, and S3 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 Sl monitoring position and with equation as follows:
2 L a + 3 4 L a 1 = T v
where L = 36 m is the pipe length and T v is the arrival time of the pressure decrease wave at the S1 monitoring position.
The pressure increase wave travels the length of 2 L 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.75 L from the valve to the S1 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 S1 monitoring position, T v , is identified with the rarefaction wave contours, as shown in Figure 27 and Figure 28. If T v = 0.0774 s is identified as the rarefaction wave arrival time, the result is a 1 = 1277 m/s, while if T v = 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 S3 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 S3 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 (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   s = T c + 2 L / a + 0.11 L / a 1 , the pressure at the S3 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 S3 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 = T c + 2 L / a + 0.12 L / a 1 ). The maximum length of the cavity region along the pipe is 3.73 m, and it occurs at T 2 = 0.078 s = T c + 2 L / a + 0.18 L / a 1 . The cavity region growth stage is from T0 = 0.0757 s to T2 = 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   s = T c + 2 L / a + 0.47 L / a 1 and T 4 = 0.097   s = T c + 2 L / a + 0.78 L / 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 = T c + 2 L / a + 1.93 L / a 1 and Figure 34d T 6 = 0.1354   s = T c + 2 L / a + 1.99 L / 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 T3 = 0.087 s to T6 = 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 = K l [ ρ l ( ρ l ρ v ) α v ] [ 1 + ( K l K v 1 ) α v ]
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.
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 wave speed value of 1277 m/s, when the rarefaction wave arrives at the S1 monitoring position at T v = 0.0774 s (illustrated in Section 3.4.1). However, the average wave speed from Figure 36 at T0 = 0.0757 s, T1 = 0.0761 s and T6 = 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.

Author Contributions

Conceptualization, X.T. and X.D.; methodology, X.D. and H.G.; software, H.G.; validation, X.D., X.T. and H.G.; formal analysis, X.D.; investigation, X.D., X.L. and X.S.; resources, H.G., X.L. and X.S.; data curation, X.D. and X.S.; writing—original draft preparation, X.D.; writing—review and editing, X.D.; visualization, H.G. and X.L.; supervision, X.T.; project administration, X.T.; funding acquisition, X.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 51779257, 51479196, 51179192. The APC was funded by the National Natural Science Foundation of China, grant number 51779257, 51479196, 51179192.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Chaudhry, M.H. Applied Hydraulic Transients; Springer: Cham, Switzerland, 2014. [Google Scholar]
  2. Simpson, A.R.; Bergant, A. Numerical Comparison of Pipe-Column-Separation Models. J. Hydraul. Eng. 1994, 120, 361–377. [Google Scholar] [CrossRef] [Green Version]
  3. Bergant, A.; Simpson, A.R. Pipeline Column Separation Flow Regimes. J. Hydraul. Eng. 1999, 125, 835–848. [Google Scholar] [CrossRef] [Green Version]
  4. Urbanowicz, K.; Zarzycki, Z.; Kudzma, S. Universal Weighting Function in Modeling Transient Cavitating Pipe Flow. J. Theor. Appl. Mech. 2012, 50, 889–902. [Google Scholar]
  5. Zhu, Y.; Duan, H.F.; Li, F.; Wu, C.G.; Yuan, Y.X.; Shi, Z.F. Experimental and Numerical Study on Transient Air–Water Mixing Flows in Viscoelastic Pipes. J. Hydraul. Res. 2018, 56, 877–887. [Google Scholar] [CrossRef]
  6. Gao, H.; Tang, X.; Li, X.; Shi, X. Analyses of 2D Transient Models for Vaporous Cavitating Flows in Reservoir-Pipeline-Valve Systems. J. Hydroinform. 2018, 20, 934–945. [Google Scholar] [CrossRef] [Green Version]
  7. Jang, T.U.; Wu, Y.; Xu, Y.; Newman, J.; Sun, Q. Efficient Quasi-Two-Dimensional Water Hammer Model on A Characteristic Grid. J. Hydraul. Eng. 2016, 142, 06016019. [Google Scholar] [CrossRef]
  8. Martins, N.M.C.; Delgado, J.N.; Ramos, H.M.; Newman, J. Maximum transient pressures in a rapidly filling pipeline with entrapped air using a CFD model. J. Hydraul. Res. 2017, 55, 506–519. [Google Scholar] [CrossRef]
  9. Martins, N.M.C.; Soares, A.K.; Ramos, H.M.; Covas, D.I.C. CFD modeling of transient flow in pressurized pipes. Comput. Fluids 2016, 126, 129–140. [Google Scholar] [CrossRef]
  10. Wang, C.; Nilsson, H.; Yang, J.; Petit, O. 1D–3D coupling for hydraulic system transient simulations. Comput. Phys. Commun. 2017, 210, 1–9. [Google Scholar] [CrossRef]
  11. Wu, D.; Yang, S.; Wu, P.; Wang, L. MOC-CFD Coupled Approach for the Analysis of the Fluid Dynamic Interaction between Water Hammer and Pump. J. Hydraul. Eng. 2015, 141, 06015003. [Google Scholar] [CrossRef]
  12. Yang, S.; Wu, D.; Lai, Z.; Du, T. Three-Dimensional Computational Fluid Dynamics Simulation of Valve-Induced Water Hammer. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2017, 231, 2263–2274. [Google Scholar] [CrossRef]
  13. Saemi, S.; Raisee, M.; Cervantes, M.J.; Nourbakhsh, A. Computation of Two- and Three-Dimensional Water Hammer Flows. J. Hydraul. Res. 2019, 57, 386–404. [Google Scholar] [CrossRef]
  14. Mandair, S.; Karney, B.; Magnan, R.; Morissette, J.F. Comparing Pure CFD and 1-D Solvers for the Classic Water Hammer Models of a Pipe-Reservoir System. In Proceedings of the WDSA/CCWI Joint Conference, Kingston, ON, Canada, 23–25 July 2018; Volume 1. [Google Scholar]
  15. Su, H.; Luo, Y.; Qu, L.; Feng, G. 3D CFD Simulation of Water Hammer in High Pressure Common Rail System of Diesel Engine. CSCBD 2019, 2019. [Google Scholar] [CrossRef] [Green Version]
  16. Heng, Z.; Zhijie, Z.; Wenbo, P. Research on Water Hammer Phenomenon during Stop Valve Closing Process Based on CFD. In Proceedings of the International Symposium on Big Data and Artificial Intelligence, Hong Kong, China, 29–30 December 2018. [Google Scholar] [CrossRef]
  17. Wang, J.; Wei, Z.; Gou, J.; Xiao, Q.; Li, S.; Li, Y. CFD Numerical Simulation of Water Hammer in a Vortex Diode. In Proceedings of the 2018 26th International Conference on Nuclear Engineering, London, UK, 22–26 July 2018. [Google Scholar] [CrossRef]
  18. Ioriatti, M.; Dumbser, M.; Iben, U. A Comparison of Explicit and Semi-Implicit Finite Volume Schemes for Viscous Compressible Flows in Elastic Pipes in Fast Transient Regime. ZAMM J. Appl. Math. Mech. 2017. [Google Scholar] [CrossRef]
  19. Urbanowicz, K. Fast and Accurate Modelling of Frictional Transient Pipe Flow. ZAMM J. Appl. Math. Mech. 2018, 98, 802–823. [Google Scholar] [CrossRef]
  20. Schohl, G.A. Improved Approximate Method for Simulating Frequency-Dependent Friction in Transient Laminar Flow. J. Fluids Eng. 1993, 115, 420–424. [Google Scholar] [CrossRef]
  21. Vardy, A.E.; Brown, J.M.B.; He, S.; Ariyaratne, C.; Gorji, S. Applicability of Frozen-Viscosity Models of Unsteady Wall Shear Stress. J. Hydraul. Eng. 2014, 141, 04014064. [Google Scholar] [CrossRef] [Green Version]
  22. Li, Z.; Bi, H.; Karney, B.; Wang, Z. Three-dimensional transient simulation of a prototype pump-turbine during normal turbine shutdown. J. Hydraul. Res. 2017, 55, 520–537. [Google Scholar] [CrossRef]
  23. Zhang, X.; Cheng, Y.; Xia, L.; Yang, J. CFD simulation of reverse water-hammer induced by collapse of draft-tube cavity in a model pump-turbine during runaway process. IOP Conference Series: Earth and Environmental Science. IOP Publ. 2016, 49, 052017. [Google Scholar]
  24. Zhang, J.; Zhang, Y.W. Analysis on hydrodynamic characteristic and cavity form of high-speed projectile with small angle of attack. In Proceedings of the Second International Conference on Mechanic Automation & Control Engineering, Hohhot, China, 15–17 July 2011. [Google Scholar]
  25. Schnerr, G.H.; Sauer, J. Physical and numerical modeling of unsteady cavitation dynamics. In Proceedings of the Fourth International Conference on Multiphase Flow, New Orleans, LA, USA, 27 May–1 June 2001. [Google Scholar]
  26. Song, C.C.S.; Yuan, M. A weakly compressible flow model and rapid convergence methods. J. Fluids Eng. 1988, 110, 441–445. [Google Scholar] [CrossRef]
  27. Coutier-Delgosha, O.; Fortes-Patella, R.; Reboud, J.L. Simulation of unsteady cavitation with a two-equation turbulence model including compressibility effects. J. Turbul. 2002, 3, 58–65. [Google Scholar] [CrossRef] [Green Version]
  28. Coutier-Delgosha, O.; Reboud, J.L.; Delannoy, Y. Numerical simulation of the unsteady behaviour of cavitating flows. Int. J. Numer. Methods Fluids 2003, 42, 527–548. [Google Scholar]
  29. Li, D.Q.; Grekula, M.; Lindell, P. A modified SST k-ω turbulence model to predict the steady and unsteady sheet cavitation on 2D and 3D hydrofoils. In Proceedings of the 7th International Symposium on Cavitation, Ann Arbor, MI, USA, 16–20 August 2009. [Google Scholar]
  30. Li, Z.; Pourquie, M.; Van Terwisga, T.J.C. A numerical study of steady and unsteady cavitation on a 2D hydrofoil. J. Hydrodyn. 2010, 22, 728–735. [Google Scholar] [CrossRef]
  31. Tullis, J.P. Fundamentals of Cavitation; Springer: Dordrecht, The Netherlands, 2005. [Google Scholar]
  32. Simpson, A.R. Large Water Hammer Pressure Due to Column Separation in Sloping Pipes; University of Michigan: Ann Arbor, MI, USA, 1986. [Google Scholar]
  33. Wylie, E.B.; Streeter, V.L. Fluid Transients; McGraw-Hill Book Co.: New York, NY, USA, 1983. [Google Scholar]
Figure 1. Sketch of reservoir-pipe-valve experiment device.
Figure 1. Sketch of reservoir-pipe-valve experiment device.
Water 12 00448 g001
Figure 2. Sketch of the geometry model.
Figure 2. Sketch of the geometry model.
Water 12 00448 g002
Figure 3. Sketch of the grid.
Figure 3. Sketch of the grid.
Water 12 00448 g003
Figure 4. Comparison of the cross-sectional average pressure fluctuations at the S3 position in different grids.
Figure 4. Comparison of the cross-sectional average pressure fluctuations at the S3 position in different grids.
Water 12 00448 g004
Figure 5. Comparison of cross-sectional average pressure fluctuations at the S3 position of different time steps.
Figure 5. Comparison of cross-sectional average pressure fluctuations at the S3 position of different time steps.
Water 12 00448 g005
Figure 6. Comparison of experimental data and numerical simulation results of pressure head fluctuation at each monitoring point. (a) S1, (b) S2, (c) S3.
Figure 6. Comparison of experimental data and numerical simulation results of pressure head fluctuation at each monitoring point. (a) S1, (b) S2, (c) S3.
Water 12 00448 g006aWater 12 00448 g006b
Figure 7. Pressure field evolution during the first wave cycle on the axis plane of the pipeline.
Figure 7. Pressure field evolution during the first wave cycle on the axis plane of the pipeline.
Water 12 00448 g007
Figure 8. Velocity field evolution during the first wave cycle on the axis plane of the pipeline.
Figure 8. Velocity field evolution during the first wave cycle on the axis plane of the pipeline.
Water 12 00448 g008
Figure 9. Experimental data and numerical simulation results of pressure head fluctuation at S1 monitoring point.
Figure 9. Experimental data and numerical simulation results of pressure head fluctuation at S1 monitoring point.
Water 12 00448 g009
Figure 10. Axial velocity distribution at different time ( T 0 = 0 s, T 1 = 0.044 s, and T 2 = 0.05 s) of the S1 monitoring point.
Figure 10. Axial velocity distribution at different time ( T 0 = 0 s, T 1 = 0.044 s, and T 2 = 0.05 s) of the S1 monitoring point.
Water 12 00448 g010
Figure 11. Local pressure distribution at S1 monitoring position of different time ( T 0 = 0 s, T 1 = 0.044 s, and T 2 = 0.05 s).
Figure 11. Local pressure distribution at S1 monitoring position of different time ( T 0 = 0 s, T 1 = 0.044 s, and T 2 = 0.05 s).
Water 12 00448 g011
Figure 12. Pressure field evolution on the axis plane of the pipeline ( T 0 = 0 s, T 1 = 0.044 s, and T 2 = 0.05 s).
Figure 12. Pressure field evolution on the axis plane of the pipeline ( T 0 = 0 s, T 1 = 0.044 s, and T 2 = 0.05 s).
Water 12 00448 g012
Figure 13. Axial velocity distribution at different time ( T 3 = 0.058 s and T 4 = 0.064 s) of the S1 monitoring point.
Figure 13. Axial velocity distribution at different time ( T 3 = 0.058 s and T 4 = 0.064 s) of the S1 monitoring point.
Water 12 00448 g013
Figure 14. Local pressure distribution at S1 monitoring position of different time ( T 3 = 0.058 s and T 4 = 0.064 s).
Figure 14. Local pressure distribution at S1 monitoring position of different time ( T 3 = 0.058 s and T 4 = 0.064 s).
Water 12 00448 g014
Figure 15. Pressure field evolution on the axis plane of the pipeline ( T 3 = 0.058 s and T 4 = 0.064 s).
Figure 15. Pressure field evolution on the axis plane of the pipeline ( T 3 = 0.058 s and T 4 = 0.064 s).
Water 12 00448 g015
Figure 16. Axial velocity distribution at different time ( T 5 = 0.092 s, T 6 = 0.1 s, and T 7 = 0.108 s) of the S1 monitoring point.
Figure 16. Axial velocity distribution at different time ( T 5 = 0.092 s, T 6 = 0.1 s, and T 7 = 0.108 s) of the S1 monitoring point.
Water 12 00448 g016
Figure 17. Local pressure distribution at S1 monitoring position of different time ( T 5 = 0.092 s, T 6 = 0.1 s, and T 7 = 0.108 s).
Figure 17. Local pressure distribution at S1 monitoring position of different time ( T 5 = 0.092 s, T 6 = 0.1 s, and T 7 = 0.108 s).
Water 12 00448 g017
Figure 18. Pressure field evolution on the axis plane of the pipeline ( T 5 = 0.092 s, T 6 = 0.1 s, and T 7 = 0.108 s).
Figure 18. Pressure field evolution on the axis plane of the pipeline ( T 5 = 0.092 s, T 6 = 0.1 s, and T 7 = 0.108 s).
Water 12 00448 g018
Figure 19. Axial velocity distribution at different time ( T 8 = 0.116 s and T 9 = 0.124 s) of S1 monitoring point.
Figure 19. Axial velocity distribution at different time ( T 8 = 0.116 s and T 9 = 0.124 s) of S1 monitoring point.
Water 12 00448 g019
Figure 20. Local pressure distribution at S1 monitoring position of different time ( T 8 = 0.116 s and T 9 = 0.124 s).
Figure 20. Local pressure distribution at S1 monitoring position of different time ( T 8 = 0.116 s and T 9 = 0.124 s).
Water 12 00448 g020
Figure 21. Pressure field evolution on the axis plane of the pipeline ( T 8 = 0.116 s and T 9 = 0.124 s).
Figure 21. Pressure field evolution on the axis plane of the pipeline ( T 8 = 0.116 s and T 9 = 0.124 s).
Water 12 00448 g021
Figure 22. Axial velocity distribution at different time ( T 10 = 0.148 s and T 11 = 0.164 s) of the S1 monitoring point.
Figure 22. Axial velocity distribution at different time ( T 10 = 0.148 s and T 11 = 0.164 s) of the S1 monitoring point.
Water 12 00448 g022
Figure 23. Local pressure distribution at S1 monitoring position of different time ( T 10 = 0.148 s and T 11 = 0.164 s).
Figure 23. Local pressure distribution at S1 monitoring position of different time ( T 10 = 0.148 s and T 11 = 0.164 s).
Water 12 00448 g023
Figure 24. Pressure field evolution on the axis plane of the pipeline ( T 10 = 0.148 s and T 11 = 0.164 s).
Figure 24. Pressure field evolution on the axis plane of the pipeline ( T 10 = 0.148 s and T 11 = 0.164 s).
Water 12 00448 g024
Figure 25. Turbulence dissipation on S1 monitoring position local region at T = 0 s and T = 0.082 s.
Figure 25. Turbulence dissipation on S1 monitoring position local region at T = 0 s and T = 0.082 s.
Water 12 00448 g025
Figure 26. Comparison between Computational Fluid Dynamics (CFD) simulation results and experimental data.
Figure 26. Comparison between Computational Fluid Dynamics (CFD) simulation results and experimental data.
Water 12 00448 g026
Figure 27. T = 0.0774 s Pressure field on the axis plane of the pipeline (a) and local pressure distribution around S1 monitoring position (b).
Figure 27. T = 0.0774 s Pressure field on the axis plane of the pipeline (a) and local pressure distribution around S1 monitoring position (b).
Water 12 00448 g027
Figure 28. T = 0.08 s Pressure field on the axis plane of the pipeline (a) and local pressure distribution around S1 monitoring position (b).
Figure 28. T = 0.08 s Pressure field on the axis plane of the pipeline (a) and local pressure distribution around S1 monitoring position (b).
Water 12 00448 g028
Figure 29. Pressure head change and vapor volume fraction change at S3 monitoring position (0.02 m upstream of the valve).
Figure 29. Pressure head change and vapor volume fraction change at S3 monitoring position (0.02 m upstream of the valve).
Water 12 00448 g029
Figure 30. Evolution of vapor volume fraction on the axis plane of the pipeline in the cavity growth stage.
Figure 30. Evolution of vapor volume fraction on the axis plane of the pipeline in the cavity growth stage.
Water 12 00448 g030
Figure 31. Pressure field evolution on the axis plane of the pipeline in the cavity growth stage.
Figure 31. Pressure field evolution on the axis plane of the pipeline in the cavity growth stage.
Water 12 00448 g031
Figure 32. Velocity field evolution on the axis plane of the pipeline in the cavity growth stage.
Figure 32. Velocity field evolution on the axis plane of the pipeline in the cavity growth stage.
Water 12 00448 g032
Figure 33. Evolution of vapor volume fraction on the axis plane of the pipeline in the cavity collapse stage.
Figure 33. Evolution of vapor volume fraction on the axis plane of the pipeline in the cavity collapse stage.
Water 12 00448 g033
Figure 34. Pressure field evolution on the axis plane of the pipeline in the cavity collapse stage.
Figure 34. Pressure field evolution on the axis plane of the pipeline in the cavity collapse stage.
Water 12 00448 g034
Figure 35. Velocity field evolution on the axis plane of the pipeline in the cavity collapse stage.
Figure 35. Velocity field evolution on the axis plane of the pipeline in the cavity collapse stage.
Water 12 00448 g035
Figure 36. Pressure wave speed field evolution in the transient cavitation flow.
Figure 36. Pressure wave speed field evolution in the transient cavitation flow.
Water 12 00448 g036
Table 1. Main parameters of different experimental schemes.
Table 1. Main parameters of different experimental schemes.
Experiment No.Initial Velocity (m/s)Initial Reynolds NumberWater Level of Upstream Reservoir (m)
Experiment 10.239453124.30
Experiment 20.332629423.41
Table 2. Boundary conditions of the calculation cases.
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 10.239237,739.7235,960.424
Case 20.332229,060.1225,776.916
Table 3. Time nodes of T 1 T 12 in terms of T c and L / a .
Table 3. Time nodes of T 1 T 12 in terms of T c and L / a .
Time NodeTime(s)In Terms of T c and L / a
T 1 0.024 T c + 0 × L / a
T 2 0.036 T c + 0 × L / a
T 3 0.052 T c + 1 × L / a
T 4 0.06 T c + 1.28 × L / a
T 5 0.072 T c + 1.71 × L / a
T 6 0.08 T c + 1.99 × L / a
T 7 0.086 T c + 2.2 × L / a
T 8 0.098 T c + 2.63 × L / a
T 9 0.108 T c + 2.99 × L / a
T 10 0.12 T c + 3.41 × L / a
T 11 0.13 T c + 3.77 × L / a
T 12 0.136 T c + 3.98 × L / a

Share and Cite

MDPI and ACS Style

Tang, X.; Duan, X.; Gao, H.; Li, X.; Shi, X. CFD Investigations of Transient Cavitation Flows in Pipeline Based on Weakly-Compressible Model. Water 2020, 12, 448. https://doi.org/10.3390/w12020448

AMA Style

Tang X, Duan X, Gao H, Li X, Shi X. CFD Investigations of Transient Cavitation Flows in Pipeline Based on Weakly-Compressible Model. Water. 2020; 12(2):448. https://doi.org/10.3390/w12020448

Chicago/Turabian Style

Tang, Xuelin, Xiangyu Duan, Hui Gao, Xiaoqin Li, and Xiaoyan Shi. 2020. "CFD Investigations of Transient Cavitation Flows in Pipeline Based on Weakly-Compressible Model" Water 12, no. 2: 448. https://doi.org/10.3390/w12020448

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop