3D CFD simulation and analysis of transient flow in a water pipeline

The following symbols are used in this paper:

a

wave speed (m/s)

C_d

discharge coefficient in MOC

D

inner diameter of the pipe (m)

E

Young modulus of the elastic pipe (Pa)

e

pipe wall thickness (m)

external body forces (N)

G_k

the production of turbulence kinetic energy (kg/(m·s³))

G_ω

the generation of ω (kg/(m·s³))

gravitational acceleration (m/s²)

ΔH

head loss (m)

K_f

Bulk modulus of the fluid elasticity (Pa)

Bulk modulus of the pipe (Pa)

K_v

head loss coefficient

k

turbulent kinetic energy (m²/s²)

p

pressure (Pa)

t

time (s)

t₁

ball valve's closure time (s)

V

mean velocity (m/s)

v_t

turbulent kinematic viscosity (m²/s)

velocity components (m/s)

x_i, x_j

Cartesian coordinate numbers

Γ_k

the effective diffusivity of k

Kronecker delta

Γ_ω

the effective diffusivity of ω

θ

rotation angle (°)

θ_max

maximum closure angle (°)

λ_a

sweep/D, mesh size ratio in the axial direction to pipe diameter

λ_c

MMS/πD, mesh size ratio in circumferential direction

λ_r

the ratio of the first-layer thickness (FLT) to the thickness of the turbulent-flow viscous sublayer in the radial direction

ρ

fluid density (kg/m³)

stress tensor

ω

specific dissipation rate (1/s)

ω_r

angular velocity (rad/s)

1D

one-dimensional

3D

three-dimensional

CFD

computational fluid dynamics

FLT

first-layer thickness

MOC

method of characteristics

MMS

mesh maximum size

SIMPLE

Semi-Implicit Method for Pressure-Linked Equations

SST

shear–stress transport

Hydraulic transients often occur in a variety of pressurized pipe systems, including basin long-distance water transfer projects, hydropower and pump stations, urban water supply systems, and agricultural irrigation systems. Water-hammer events caused by accidental pump or turbine shutdowns or quick valve actions, likely produce abnormally high pressures, which may induce pipe rupture and damage other hydraulic devices (Ghidaoui & Kolyshkin 2001; Saemi et al. 2019). Therefore, the accurate numerical simulation of the water-hammer events is crucial for the proper design and safe operation of the pipe systems.

Traditional one-dimensional (1D) water-hammer models, in which shear stress along the pipeline is approximated by using the steady-state dissipation formulas, tend to underestimate the energy dissipations. To improve the numerical accuracy of water-hammer pressure, various 1D unsteady pipe wall shear stress formulas have been developed. According to the review of Ghidaoui et al. (2005), wall shear stress models can be classified into three types: (i) quasi-steady wall shear models. Wall shear stress expressions (such as Darcy–Weisbach and Hazen–Williams formulas) are assumed to remain valid under unsteady conditions. (ii) Empirical-based corrections to quasi-steady wall shear models. This type of 1D unsteady friction model is based on the instantaneous acceleration of the fluid (Brunone et al. 1991). Daily et al. (1956) developed a model in which the unsteady friction is dependent on instantaneous mean flow velocity and instantaneous local acceleration. Carstens & Roller (1959) presented the influence of instantaneous mean flow velocity and instantaneous local acceleration on friction terms. The instantaneous-acceleration unsteady friction model is simple to implement and computationally efficient but is difficult to determine the empirical parameter. (iii) Physically-based wall shear models. This type of 1D unsteady friction model is the convolution-based approach based on the model of Zielke (1968), who proposed a physically based approach for the derivation of unsteady wall shear in transient laminar flows by determining the convolution of previous fluid accelerations and a weighting function. Vardy & Brown (1995, 2003, 2004b) extended Zielke's approach to turbulent water-hammer flows. Original convolution-based models, Zielke's model, and the Vardy–Brown model are inefficient but have been shown to be able to accurately match computed to measured results. Based on the approximation method of Trikha (1975), Vardy & Brown (2004a), Vitkovsky et al. (2004), and Urbanowiez & Zarzycki (2012) successively developed various efficient and accurate calculations of Zielke and Vardy–Brown unsteady friction in pipe transients. Recently, Urbanowiez et al. (2016) developed formulas to enable the determination of shear stresses on the pipe walls together with a formula for coefficients of non-stationary friction losses. Urbanowiez (2017) presented certain analytical formulas to determine the coefficients of simplified effective weighting functions in a simple direct way. Urbanowiez (2018) presented a correction to the erratic recursive formula, which was used to calculate the unsteady wall shear stresses.

For transient pipe flow, the generation of vorticity at the pipe wall and its subsequent diffusion have been observed experimentally by Das & Arakeri (1998) and Brunone et al. (2000) and found numerically by Vardy & Hwang (1991), Silva-Araya & Chaudhry (1997), Pezzinga (1999), and Ghidaoui & Mansour (2002). However, it is nearly impossible for the 1D approach to accurately simulate vorticity at the pipe wall and describe the relations between dynamic vorticity and energy dissipation.

Moreover, the standard 1D approach in dealing with the dynamic hydraulic losses caused by the closing of the valve in unsteady pipe fluid flow is to assume the energy losses to be the same as if the flows were steady at the same velocity. Importantly, transient dynamic characteristics of the closing ball valve, and the unstable flow characteristics of the local vortex and backflow during the closing process of the valve are complicated three-dimensional (3D) flow problem, which is difficult to be determined by the 1D methods or the experimental measures.

Martins et al. (2014) introduced Reynolds-averaged Navier–Stokes (RANS) equations and a realizable k–ε turbulence model to realize 3D unsteady pipe-flow simulations. Subsequently, Martins et al. (2016) investigated the velocity distribution and wall shear stress by using the 3D computational fluid dynamics (CFD) model. Martins et al. (2017, 2018) explored the complex nature of the transient energy dissipation in the low Reynolds number turbulent flow. The authors (Wang et al. 2016; Zhou et al. 2018) have developed CFD models to investigate liquid column separation and transient flow with an entrapped air pocket in a pipeline, respectively. Ferriera et al. (2018) used the 3D CFD model to analyze the behavior of the ball valve under steady and unsteady conditions, which was compared with the measured results. Results show the 3D CFD methods are feasible to accurately predict the transient pressures, and its significant advantage is to vividly reveal the physical process.

With respect to previous papers (Martins et al. 2014; Ferriera et al. 2018), the primary novelties of this paper are as follows: (1) one alternative 3D CFD model is presented to systematically explore the 3D dynamic flow characteristics in pipe transients; (2) the detailed analysis of the ball-valve characteristics under static and dynamic conditions are carried out by the proposed 3D CFD approach; and (3) the effects of closure time on head loss coefficient and discharge coefficient are carefully investigated. The 3D CFD water-hammer approach is developed to investigate the transient pipe flow and the 3D detailed dynamic characteristics of the closing ball valve. As discussed above, existing studies primarily concern the 1D pipe transient simulations or the valve flow characteristics based on the 3D steady-flow incompressible-water simulations. Here, the proposed 3D CFD approach accounts for the water compressibility and the viscous sublayer, which are sometimes neglected in previous studies. 3D CFD pressure results are compared with both the measured data and 1D results (which are coded and run in Visual Studio 2010). The 3D transient flow characteristics of the ball valve during the closing process are further investigated by using the validated 3D CFD model. Moreover, the effects of the valve closing time on the valve characteristics are carefully discussed.

The experimental apparatus in the work of Bergant et al. (2001) is used here to validate the proposed 3D model and to further investigate the water-hammer events in the pipeline.

As shown in Figure 1, the experimental system consists of an upstream pressure tank with a constant pressure head of 32 m, a 37.23 m long copper pipe with an inner diameter of 22.1 mm, a downstream one-quarter-turn ball valve, and a downstream pressure tank. The wave speed was measured at 1,319 m/s. Three experimental cases are considered here. In these three tests, initially, the ball valve is totally open, and various steady flows with different velocities V₀ are formed by adjusting the pressure of the downstream pressure tank. Water-hammer events were induced by suddenly closing the ball valve (closure time 0.009 s). Here, Case 1: laminar flow, initial pipe velocity V₀ = 0.1 m/s, and Reynolds number Re = 1,870; Case 2: transitional flow, V₀ = 0.2 m/s, and Re = 3,750; and Case 3: turbulent flow, V₀ = 0.3 m/s, and Re = 5,600.

The standard two-equation eddy-viscosity models (k–ε and k–ω class) based on the Boussinesq relation are originally developed for the pipe flows. For the fast transient compressible-fluid flow, they often overpredict the turbulence viscosity, in particular for the more complicated water-hammer problems (i.e., cavitating pipe flow) (Wang et al. 2016). Finally, the shear-stress transport (SST) k–ω model is implemented in this work, since it can better predict the transient pressure during the valve closing process, which is also demonstrated in the authors’ previous work (Wang et al. 2016). Importantly, the SST k–ω model is robust and easy to include the refined mesh and the influence of the boundary layer, which have important effects on shear stress (discussed later in this paper).

The solution of these governing equations could be efficiently implemented in the general CFD software (i.e., OpenFOAM, ANSYS FLUENT, etc). In this work, ANSYS FLUENT 18.0 is used to realize the 3D simulations. The pressure–velocity coupling algorithm, Semi-Implicit Method for Pressure-Linked Equations (SIMPLE), is used to solve the pressure, the velocity, and the density equations in compressible flows. The discretization is based on a finite volume approach. The pressure, the density, and the momentum equations are solved by the second-order method, while the turbulent kinetic energy and specific dissipation rate are computed by the first-order upwind method. Moreover, the first-order implicit transient formulation is selected to realize time-dependent calculations.

Based on previous studies (Zielke 1968; Martins et al. 2014), the unsteady friction factor (or transient shear stress) is a crucial influence factor for precisely predicting pressure oscillations. Therefore, in order to accurately elaborate the shear stress near the pipe wall, the boundary layer and its refined mesh near the pipe wall are carefully treated in this work. In the current 3D simulations, the structured hex cells are adopted. Three grid density parameters are defined here, λ_a= sweep/D, the ratio of the mesh size in the axial direction to pipe diameter, and sweep is ‘the mesh size in the axial direction of the pipe’; λ_c= MMS/πD, the ratio of the mesh size in the circumferential direction to the perimeter, and MMS is ‘mesh maximum size’; λ_r= FLT/y₅, the ratio of the first layer thickness (FLT) to the thickness of the turbulent-flow viscous sublayer in the radial direction.

To ensure both the high accuracy and the economical computation, λ_a = 2.09, λ_c = 0.02083, and λ_r = 0.18 are recommended in the 3D water-hammer simulations (which is recommended by Martins et al. 2014). Moreover, it is necessary to refine the grid within the complicated flow region of the moving ball valve. Finally, total of 780,000 cells are suggested for the current 3D simulations (in Figure 2). The grid density and grid number are discussed later.

The accuracy and efficiency of 3D CFD simulation results depend on the mesh quality. Here, a series of grid independency studies are performed to determine the most suitable mesh size. Importantly, the viscous sublayer is cautiously considered in this work to explore and evaluate the effects of pipe shear stress on other physical variables.

As discussed in the section ‘3D grids and boundary conditions’, some definitions of grid parameters (e.g., λ_a, λ_c, and λ_r) are similar to those proposed by Martins et al. (2014), are introduced to describe the grid density in different dimensional directions.

A great number of meshing schemes were conducted to investigate the effects of grid density (λ_a, λ_c, and λ_r). It was found that (a) as the grid densities (λ_a, λ_c, and λ_r) increase, both the numerical accuracy and computational cost increase; (b) when grid densities are near the suggested values (λ_a = 2.09, λ_c = 0.02083, and λ_r = 0.18), the 3D CFD simulations can reach to the high accuracy and the economical computation; (c) when being far larger than the recommended values, higher grid density cannot further improve the accuracy; on the contrary, it increases the computations and even leads to numerical instability; and (d) compared with the mesh size in the axial direction, the variation of the mesh size in cross section is more sensitive to the numerical accuracy. Considering the page space limitation, only three kinds of meshing schemes are chosen to present the above related conclusions. The detailed mesh information and computation results are shown in Table 1 and Figure 3.

This section has several goals: (1) to verify the proposed 3D CFD model; (2) to elaborate on 3D dynamic flow characteristics of transient pipe flow including the 3D distributions and internal relations of shear stress, velocity, and pressure; (3) to compare the pressure results predicted by the 3D model and the 1D model; (4) to explore the influences of the 3D meshing method and the mesh density; and (5) to investigate the effects of the closing-valve time on the flow characteristics of the closing valve. First, comparisons of the 3D CFD and experimental pressure oscillations demonstrate that the proposed 3D CFD approach is capable of accurately predicting the transient pressure of the water hammers. After that, the 3D meshing method is discussed. Finally, the 3D CFD model is used to explore the nature and mechanism of water-hammer results and the 3D dynamic characteristics of closing the ball valve.

Figure 4 displays the comparisons of the 3D pressure results and the experimental pressures in Case 1, Case 2, and Case 3. Note that the pressures were measured at 0.02 m upstream of the closing ball valve in all three experimental cases.

Figure 4 demonstrates that the 3D simulation pressure results are in excellent visual agreement with the experimental data. However, compared with 3D CFD pressure oscillations, the experimental results in all three cases display a slight time lag after the fifth cycle (the maximum lag time is: case 1 0.0084 s; case 2 0.0094 s; case 3 0.0091 s). The possible reasons are that the wave cycles keep constant in the numerical simulations, however, it can be found from the three pressure oscillations in Figure 4 that the actual wave speed probably decreases a little bit as the transient even progresses.

Overall, the proposed 3D CFD model can accurately reproduce the maximum pressure, the subsequent pressure peaks and periods, and the pressure damping and attenuation. The most important benefit is that the 3D CFD model can clearly visualize the changes of the flow structure during the closing process of the valve and is importantly helpful to explore the 3D complicated flow characteristics, which are not well understood by using the experimental and 1D methods. Consequently, the proposed 3D CFD model is a robust method to simulate the transient pipe flow.

It is difficult to accurately obtain transient parameter distributions (velocity, pressure, and shear stress) on pipe cross section by using both the experimental measurement and the existing 1D simulation. Fortunately, the validated 3D CFD approach developed in this paper is convenient to show all transient flow information of any point at any time including the shear stress, velocity, and pressure. Taking Case 3 as an example, Figure 6(a) and 6(b) display the pressure and velocity oscillations based on the mean values of pipe across the section in the midsection of the pipeline. As shown in Figure 5, five feature points (FP 1–5) near the turning points of both pressure and velocity oscillations are chosen to further analyze the velocity, pressure, and shear stress profiles along the pipe across the section displayed in Figure 6. Moreover, Figure 7 gives the wall shear stress time history during the transient process in laminar flow (Case 1) and turbulent flow (Case 3).

First, the 3D results in Figure 5 confirm that the variations between mean pressure head and mean velocity completely accord with energy conservation, which is consistent with the 1D analysis of the water-hammer event in the classic reference book (Wylie et al. 1993).

Figure 6(a)–6(c) reveals 3D dynamic flow characteristics in pipe transients at five feature points of Figure 5 in the midpoint of Case 3. As shown in Figure 6(a), the pressure distribution in the radial direction indicates that at any transient time point, the pressure always remains consistent from the pipe wall to the pipe core. However, the velocity and shear stress distributions in the radial direction significantly change as the pipe transient progresses (in Figure 6(b) and 6(c)). Importantly, a strong relationship between the velocity and the shear stress distributions was found during the water-hammer wave propagation. As portrayed in Figure 6, the sequence of transient events is as follows:

At feature point 1, before the high-pressure wave caused by instantaneous valve closure reaches (in Figure 5), the pipe flow at the pipeline midpoint is still identical to initial steady flow. It can be found from Figure 6(b) that the mean velocity distribution from the pipe wall to the pipe core is parabolic, and the initial velocity is zero at the pipe wall and reaches the maximum at the pipe core. That is because at the moment, as shown in Figure 6(c), shear stress reaches the maximum at the pipe wall, and it drops to zero quickly, at about 0.98 R (R: pipe radius). Shear stress is resistance to the fluid flow near the pipe wall.

At feature point 2, once the high-pressure wave caused by instantaneous valve closure firstly arrives at the pipeline midpoint, the mean pressure head rapidly increases to the peak, while the mean velocity drops to zero (in Figure 5). However, as shown in Figure 6(b), unlike the mean velocity, transient velocity is varied along the radius direction. The main reasons are that (1) before the high-pressure wave arrives at the pipeline midpoint, the water velocity near the pipe center is higher than that near the pipe wall, (2) although the high-pressure wave trends to stop the pipe flow, due to both inertia effect and water deformation, the water near the pipe centerline still moves downstream with a small velocity, and (3) meanwhile, the water near pipe wall begins to move upstream under the extrusion of the water near the pipe centerline. Correspondingly, the shear stress between 0.98 R and R changes into a positive value from a negative value (Figure 6(c)), which just resists the water backflow near the pipe wall.

At feature point 3, as the decompression pressure wave from the upstream inlet (reservoir) firstly travels to the pipeline midpoint, the mean pressure head returns to the initial condition, but the backflow occurs with about initial mean velocity. However, as shown in Figure 6(b), transient velocity along the radius direction is still inconsistent. Interestingly, the water velocity near the pipe wall is higher than that at the pipe center. The reason is that before the decompression pressure wave arrives at the pipeline midpoint, the water near the pipe center still moves downstream with a small velocity, and the water near the pipe wall moves upstream, where the mean velocity distribution is similar to that at feature point 2. Correspondingly, the shear stress between 0.98 R and R becomes larger (Figure 6(c)), in order to resist the enhanced water backflow near the pipe wall.

At feature point 4, the negative pressure wave from the downstream closed valve firstly reaches the pipeline midpoint, the mean pressure head drops to the lowest value, and the backflow nearly stops (mean velocity becomes zero). More interestingly, the mean velocity distribution consists of three distinct regions: positive flow near the wall from R to 0.96 R, negative flow in the transition region from 0.96 R to 0.7 R, and positive flow near the pipe center from 0.7 R to 0. This attributes to (1) before the negative pressure wave arrives at the pipeline midpoint, the water velocity of backflow near the pipe wall is higher than that at the pipe center, and the mean velocity distribution is nearly identical to that at the feature point 3; (2) under the negative pressure wave, the water parts near wall from R to 0.96 R and near pipe center from 0.7 R to 0, first stop the backflow, then continue to move downstream with a small positive velocity due to both inertia effect and water compressibility; (3) meanwhile, the water in the transition region from 0.96 R to 0.7 R still keeps a little backflow. Correspondingly, the shear stress between 0.98 R to R turns negative (Figure 6(c)), in order to resist the positive flow near the pipe wall.

At feature point 5, as the compression pressure wave from the upstream inlet (reservoir) first travels to the pipeline midpoint, the pressure head, and the flow velocity both almost return to the initial condition. As the transient event progresses, the behaviors from feature point 1 to feature point 5 occur periodically.

Remarkably, it is observed that shear stress is the resistance to the fluid flow near the pipe wall. At any transient time, the largest shear stress always occurs at the pipe wall, and the shear stress quickly decreases along the radius direction, reaching almost zero at about 0.98 R. Therefore, the velocity characteristics near the pipe wall can be primitively judged based on the variations of transient wall shear stress in Figure 7 (see the detailed analysis at the five feature points). Moreover, as shown in Figure 7, compared with the turbulent flow (Case 3), the laminar flow (Case 1) has the identical variation trend of wall shear stress, but in which the amplitude becomes smaller.

Figure 8 displays the comparisons of the 3D pressure results and the 1D pressure results at the 0.02 m upstream of the closing ball valve in Case 1 to Case 3. In the 1D simulations, the first-order method of characteristics (MOC) scheme is used to solve the 1D models, and the 1D grid number is 32 (Bergant et al. 2001).

Figure 8(a)–8(c) demonstrates that the 1D models can get the basically identical pressure oscillations with the 3D CFD pressure results. However, due to the 1D inherent assumptions and limitations, it is difficult for the existing 1D models to further explore the transient parameter distributions along the radius direction.

Generally, the flow characteristics of the valves are measured from physical model tests under static closing conditions. In this work, flow characteristics of both static closing conditions and dynamic closing conditions are carefully investigated by using the proposed 3D CFD approach. First, the effects of valve closure time on the velocity and pressure are analyzed; second, the velocity and pressure contours are presented to further explore the 3D characteristics of water-hammer flow events; finally, two coefficients, head loss coefficient K_v and discharge coefficient C_d, are introduced to comprehensively investigate the effects of closure time.

The initial condition and system parameters of Case 1 are used in current 3D CFD simulations. When the ball valve is fully opened, the closure angle is defined as 0°; when the ball valve is fully closed, the closure angle is 90°. The closure percentage τ is defined as: τ=θ/θ_max, where θ is closure angle and θ_max is maximum closure angle (θ_max = 90°). The maximum effective closure percentage is 94.5%. Effective closure time (closure percentage from 0 to 94.5%) is 0.0085 s.

Four tests with different valve closure time are calculated in this section. As shown in Table 2, Serial 1 is the static closing condition (valve stops at a given angle) and Serials 2–6 are dynamic closing conditions (closing the ball valve without any stop).

Velocity and pressure variations at the ball valve are displayed in Figures 9 and 10, respectively.

As shown in Figure 9, the variations of velocity undergo three distinct stages during the closing process: a slight change in Stage 1, then a drastic reduction in Stage 2, and finally slowing down in Stage 3. Correspondingly, the differential pressure at the valve has similar three stages (in Figure 10). The ranges of closure percentage at three stages are associated with the valve closure time.

For the static closing condition, Stage 1 with a closure percentage of 0–28%, the velocity decreases slightly along with the closing process and the differential pressure increase slowly. During this period, since the closure percentage is relatively small, it has little effect on the pipe flow. As the closure percentage increases from 28 to 94.5% in Stage 2, the velocity at the valve quickly reduces and differential pressure is increasing at the same time. For Stage 3 with a closure percentage of 94.5–100%, the velocity nearly reduces to zero and differential pressure decreases and then keeps constant. Moreover, Figures 9 and 10 indicate that for the relatively long closure process (63 s in Serial 4), the velocity and pressure results are almost close to those in the static closing condition. However, when the valve is closed faster (9 and 0.009 s), some of the following obvious differences are found: (1) the range of Stage 1 becomes larger since the inertia effect of the water column is more dominated by shorter closure process; (2) the velocity and pressure change more drastically at Stage 2, attributing to that the valve resistance dominates in this stage; and (3) at Stage 3, residual velocity becomes higher and the pressure fluctuates more intensely, due to the combined effects of the inertia and compressibility of the water column.

Moreover, compared with the static closing conditions, velocity and pressure fields are more complex and unstable in dynamic closing conditions. For example, Figures 11 and 12 display the velocity and pressure contours of the closure percentage of 77.8% (Stage 2) in Serial 1 and Serial 3.

As discussed in the Introduction, Ferriera et al. (2018) carefully compared the 3D CFD simulation results and the experimental data of the valve static behavior. Their work validates the 3D CFD approach and provides an important reference for the current research. Based on the work of Ferriera et al. (2018), the current work more carefully explores the difference between static and dynamic characteristics of the valve. The current results indicate that the ball valve's transient characteristics are different from static characteristics, and the characteristics also change with different closure time. As described above, the velocity changes during valve closing are divided into three stages, and further analysis is presented on the specific characteristics of simulation results at different stages.

To comprehensively study the characteristics of velocity and pressure during the closing process, the pipe flow through a valve is evaluated by the following two coefficients: head loss coefficient K_v and discharge coefficient C_d.

As shown in Figures 13 and 14, the values of K_v and C_d under dynamic closing conditions are different from those in static closing conditions, especially for the abrupt closing process (t_c = 0.009 s). It can be found that (1) for the closure percentage of 0–70%, closure time obviously influences the values of K_v and C_d; (2) when the valve is closed in 0.009 s, the value of K_v becomes higher and the corresponding value of C_d decreases; and (3) for the closure percentage 70–94.5%, the influence of closure time on K_v and C_d is slight.

The 3D CFD model is developed to investigate the water-hammer pipe flow and 3D detailed dynamic characteristics of the closing ball valve. The calculated pressure results are compared both to those from the experimental data and a 1D model. The velocity fields, shear stress distributions, and transient pressures are analyzed simultaneously. Moreover, head loss coefficient and discharge coefficient were investigated under both the static and dynamic closing-valve conditions. Based on the results, the conclusions may be drawn as follows:

• (1)

  The proposed 3D CFD model can accurately reproduce the experimental pressure oscillations including the pressure peaks and periods and the pressure damping. The most important benefit is that the 3D CFD model can clearly visualize the changes in the flow structure during the closing process of the valve and is importantly helpful to explore the 3D complicated flow characteristics, which are not well understood by using the experimental and 1D methods.

• (2)

  The meshing approach and the boundary layer treatment proposed in this work are feasible to realize both the high accuracy and the economical computation of transient pipe flow. For the relatively coarse grid, as the grid densities (λ_a, λ_c, and λ_r) increase, both the numerical accuracy and computational cost increase; however, when being far larger than the recommended values (λ_a = 2.09, λ_c = 0.02083, and λ_r = 0.18), higher grid density cannot further improve the accuracy; on the contrary, it increases the computations and even leads to numerical instability.

• (3)

  During the water-hammer events, the pressure distribution always remains consistent from the pipe wall to the pipe core; however, the velocity and shear stress distributions in the radial direction significantly change as the pipe transient progresses. The velocity characteristics near the pipe wall can be primitively judged based on the variations of transient wall shear stress. The shear stress is the resistance to the fluid flow near the pipe wall. The largest shear stress always occurs at the pipe wall, and the shear stress quickly decreases along the radius direction, reaching zero at about 0.98 R.

• (4)

  During the closing-valve process, the variations of velocity through the valve undergo three distinct stages: a slight change in Stage 1, then a drastic reduction in Stage 2, and finally slowing down in Stage 3. Correspondingly, the differential pressure at the valve has the similar three stages. The ranges of closure percentage at three stages are associated with the valve closure time. As the valve is closed faster, the range of Stage 1 becomes larger; the velocity and pressure change more drastically at Stage 2; and at Stage 3, residual velocity becomes higher and the pressure fluctuates more intensely. Moreover, the head loss coefficient K_v and discharge coefficient C_d under dynamic closing conditions are different from those in static closing conditions.

The authors gratefully acknowledge the financial support for this research from the National Natural Science Foundation of China (Grant Nos 51679066 and 51839008), the Fundamental Research Funds for the Central Universities (Grant No. 2018B43114), the Fok Ying Tong Education Foundation for Young Teachers in the Higher Education Institutions of China (Grant No. 161068), and the Open Research Fund Program of State Key Laboratory of Water Resources and Hydropower Engineering Science (Wuhan University) (Grant No. 2016SDG01).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.