A three-dimensional (3D) computational fluid dynamics (CFD) approach is developed to elaborate the water-hammer pipe flow and 3D detailed dynamic characteristics of a closing ball valve. The proposed CFD approach considers the water compressibility and the viscous sublayer, which are sometimes neglected in previous studies. Comparisons of the CFD results, the measured pressures and the one-dimensional results, demonstrate that the current 3D CFD approach better reproduces the experimental pressure oscillations while helping to visualize the associated physical processes and to further explore the 3D transient characteristics. The mean velocity distributions in the radial direction significantly change as the pipe transient progresses, which is closely associated with transient shear stress. Mean velocity variations at the valve during the closing process undergo three distinct stages: slight change, then drastic reduction, and finally slowing down. Head loss coefficient and discharge coefficient of the valve change as the valve closing time shortens.

  • The 3D dynamic flow characteristics in pipe transients are studied systematically.

  • Further analysis of the ball-valve characteristics under static and dynamic conditions carried out by 3D CFD simulations.

  • The effects of closure time on head loss coefficient and discharge coefficient are investigated.

Graphical Abstract

Graphical Abstract
Graphical Abstract

The following symbols are used in this paper:

     
  • a

    wave speed (m/s)

  •  
  • Cd

    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)

  •  
  • Gk

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

  •  
  • Gω

    the generation of ω (kg/(m·s3))

  •  
  • gravitational acceleration (m/s2)

  •  
  • ΔH

    head loss (m)

  •  
  • Kf

    Bulk modulus of the fluid elasticity (Pa)

  •  
  • Bulk modulus of the pipe (Pa)

  •  
  • Kv

    head loss coefficient

  •  
  • k

    turbulent kinetic energy (m2/s2)

  •  
  • p

    pressure (Pa)

  •  
  • t

    time (s)

  •  
  • t1

    ball valve's closure time (s)

  •  
  • V

    mean velocity (m/s)

  •  
  • vt

    turbulent kinematic viscosity (m2/s)

  •  
  • velocity components (m/s)

  •  
  • xi, xj

    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/m3)

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

Physical and experimental model

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 V0 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 V0 = 0.1 m/s, and Reynolds number Re = 1,870; Case 2: transitional flow, V0 = 0.2 m/s, and Re = 3,750; and Case 3: turbulent flow, V0 = 0.3 m/s, and Re = 5,600.

Figure 1

Schematic diagram of experimental devices.

Figure 1

Schematic diagram of experimental devices.

Close modal

Governing equations

For the water-hammer problem, the mass and momentum equations can be written as
formula
(1)
formula
(2)
where ρ is the density; t is the time; v, , and are the velocity components; p is the pressure; is the stress tensor; g is the gravitational acceleration; and F is the external body forces.
In this paper, the compressibility of the water and the elasticity of the pipe are considered to simulate the water-hammer effects. To stabilize the pressure solution for compressible flows, an extra term related to the wave speed is introduced in the pressure correction equation. The density variation is considered by coupling density and pressure,
formula
(3)
formula
(4)
where Kf is the bulk modulus of the fluid elasticity; is the bulk modulus of the pipe; D is the inner diameter of the pipe; E is the Young modulus of elastic pipe; and e is the pipe wall thickness.
The wave speed in liquid, a, can be calculated by
formula
(5)

Turbulence modeling

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

In the SST kω model, the term in Equation (2) is defined via the Boussinesq approximation,
formula
(6)
where k is the turbulent kinetic energy; vt is the turbulent kinematic viscosity; and is the Kronecker delta.
The additional transport equations, k and ω equations of the SST kω model are,
formula
(7)
formula
(8)
where Gk represents the production of turbulence kinetic energy; Gω represents the generation of ω; Γk and Γω represent the effective diffusivity of k and ω. and represent the dissipation of k and ω; represents the cross diffusion term; and are source terms.

Numerical solution method

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.

3D grids and boundary conditions

In the current simulations, as shown in Figure 1, the upstream tank and downstream tank are constant pressure inlet and outlet, and the pipe wall is a no-slip ‘Wall’. The closing of the ball valve is realized by using the moving mesh method. The initial steady pipe flow before each transient flow simulation is realized by performing 24 s steady-flow calculations. Therefore, in the whole simulations, the angular velocity of the ball valve, ωr, is
formula
(9)
where t is time; θ is rotation angle; t1 is the closure time of the ball valve.

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/y5, 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.

Figure 2

The grids in cross section and axial direction.

Figure 2

The grids in cross section and axial direction.

Close modal

Mesh independence analysis

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.

Table 1

Details of different mesh sizes

RankingCellsMesh of circular cross sectionMesh of axial direction
1 525,000 420 800 
2 780,000 624 800 
3 1,900,800 1,200 1,152 
RankingCellsMesh of circular cross sectionMesh of axial direction
1 525,000 420 800 
2 780,000 624 800 
3 1,900,800 1,200 1,152 
Figure 3

Grid independency studies.

Figure 3

Grid independency studies.

Close modal

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.

Verification of the 3D CFD model by comparing experimental results

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

Comparison between experimental data with 3D CFD results and analytical solutions: (a) Case 1; (b) Case 2; and (c) Case 3.

Figure 4

Comparison between experimental data with 3D CFD results and analytical solutions: (a) Case 1; (b) Case 2; and (c) Case 3.

Close modal

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.

3D Dynamic flow characteristics in pipe transients

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

Figure 5

Pressure and velocity oscillations, feature points at the midpoint of the pipe in Case 3.

Figure 5

Pressure and velocity oscillations, feature points at the midpoint of the pipe in Case 3.

Close modal
Figure 6

3D dynamic flow characteristics at five feature points of Figure 5 at the midpoint of Case 3, (a) pressure, (b) velocity, and (c) shear stress.

Figure 6

3D dynamic flow characteristics at five feature points of Figure 5 at the midpoint of Case 3, (a) pressure, (b) velocity, and (c) shear stress.

Close modal
Figure 7

Wall shear stress time history in laminar flow (Case 1) and turbulent flow (Case 3).

Figure 7

Wall shear stress time history in laminar flow (Case 1) and turbulent flow (Case 3).

Close modal

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.

Comparison of pressures simulated by 3D model and 1D models

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

Comparison between 1D (MOC) and 3D (FLUENT) results: (a) Case 1; (b) Case 2; and (c) Case 3.

Figure 8

Comparison between 1D (MOC) and 3D (FLUENT) results: (a) Case 1; (b) Case 2; and (c) Case 3.

Close modal

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.

Comparison of static and dynamic ball-valve characteristics

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 Kv and discharge coefficient Cd, are introduced to comprehensively investigate the effects of closure time.

Head loss coefficient Kv is a dimensionless coefficient that characterizes the head loss of pipe flow through a valve:
formula
(10)
where ΔH is head loss; V is the mean velocity in the pipe; g is the gravitational acceleration.
In 1D or quasi 2D simulations of the valve, an important parameter CdAg should be calculated based on previous experimental data. CdAg means the discharge coefficient times the area of valve opening, and Cd is defined as:
formula
(11)
The relationship between Kv and Cd is:
formula
(12)

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

Table 2

Details of closure time

Serial number 1 2 3 4 5 6 
Closure time Static 0.009 s 0.1 s 1 s 9 s 63 s 
Serial number 1 2 3 4 5 6 
Closure time Static 0.009 s 0.1 s 1 s 9 s 63 s 

Effects of closure time on the velocity and pressure

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

Figure 9

3D simulation results of the velocity history.

Figure 9

3D simulation results of the velocity history.

Close modal
Figure 10

3D simulation results of differential pressure history.

Figure 10

3D simulation results of differential pressure history.

Close modal

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.

Figure 11

Velocity contours at different closure angles: (a) static closing condition and (b) dynamic closing condition (9 s).

Figure 11

Velocity contours at different closure angles: (a) static closing condition and (b) dynamic closing condition (9 s).

Close modal
Figure 12

Pressure contours at different closure angles: (a) static closing condition and (b) dynamic closing condition (9 s).

Figure 12

Pressure contours at different closure angles: (a) static closing condition and (b) dynamic closing condition (9 s).

Close modal

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.

Effects of closure time on head loss coefficient Kv and discharge coefficient Cd

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 Kv and discharge coefficient Cd.

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

Figure 13

Simulation results of the head loss coefficient Kv.

Figure 13

Simulation results of the head loss coefficient Kv.

Close modal
Figure 14

Simulation results of the discharge coefficient Cd.

Figure 14

Simulation results of the discharge coefficient Cd.

Close modal

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 Kv and discharge coefficient Cd 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.

ANSYS FLUENT [Computer software]. ANSYS 18.0, Canonsburg, PA
.
Bergant
A.
,
Simpson
A. R.
&
Vitkovsky
J. P.
2001
Developments in unsteady pipe flow friction modeling
.
J. Hydraul. Res.
39
(
3
),
249
257
.
Brunone
B.
,
Golia
U. M.
&
Groco
M.
1991
Some remarks on the momentum equations for fast transients
. In:
Hydraulic Transients with Column Separation (9th and Last Round Table of the IAHR Group)
.
IAHR
,
Valencia
,
Spain
, pp.
201
209
.
Brunone
B.
,
Karney
B. W.
,
Mecarelli
M.
&
Ferrante
M.
2000
Velocity profiles and unsteady pipe friction in transient flow
.
J. Water Resour. Plan. Manage.
126
(
4
),
236
244
.
Carstens
M. R.
&
Roller
J. E.
1959
Boundary-shear stress in unsteady turbulent pipe flow
.
J. Hydr. Div.
85
(
2
),
67
81
.
Daily
J. W.
,
Hankey
W. L.
,
Olive
R. W.
&
Jordaan
J. M.
1956
Resistance coefficient for accelerated and decelerated flows through smooth tubes and orifices
.
Trans., ASME
78
,
1071
1077
.
Das
D.
&
Arakeri
J. H.
1998
Transition of unsteady velocity profiles with reverse flow
.
J. Fluid Mech.
374
,
251
283
.
Ferriera
J. P. B. C. C.
,
Martins
N. M. C.
&
Covas
D. I. C.
2018
Ball valve behavior under steady and unsteady conditions
.
J. Hydraul. Eng.
144
(
4
),
04018005
.
Ghidaoui
M. S.
&
Kolyshkin
A. A.
2001
Stability analysis of velocity profiles in water-hammer flows
.
J. Hydraul. Eng.
127
(
6
),
499
512
.
Ghidaoui
M. S.
&
Mansour
S. G. S.
2002
Efficient treatment of the Vardy-Brown unstrady shear in pipe transients
.
J. Hydraul. Eng.
128
(
1
),
102
112
.
Ghidaoui
M. S.
,
Zhao
M.
,
Mclnnis
D. A.
&
Axworthy
D. H.
2005
A review of water hammer theory and practice
.
Appl. Mech. Rev.
58
,
49
76
.
Martins
N. M. C.
,
Carrico
N. J. G.
,
Ramos
H. M.
&
Covas
D. I. C.
2014
Velocity-distribution in pressurized pipe flow and mesh analysis using CFD: accuracy and mesh analysis
.
Comput. Fluids
105
,
218
230
.
Martins
N. M. C.
,
Soares
A. K.
,
Ramos
H. M.
&
Covas
D. I. C.
2016
CFD modeling of transient flow in pressurized pipes
.
Comput. Fluids.
126
(
2016
),
129
140
.
Martins
N. M. C.
,
Brunone
B.
,
Meniconi
S.
,
Ramos
H. M.
&
Covas
D. I. C.
2017
CFD and 1D approaches for the unsteady friction analysis of low Reynolds number turbulent flows
.
J. Hydraul. Eng.
143
(
12
),
04017050
.
Martins
N. M. C.
,
Brunone
B.
,
Meniconi
S.
,
Ramos
H. M.
&
Covas
D. I. C.
2018
Efficient computational fluid dynamics model for transient laminar flow modeling: pressure wave propagation and velocity profile changes
.
ASME. J. Fluids Eng.
140
(
1
),
011102
.
Pezzinga
G.
1999
Quasi-2D model for unsteady flow in pipe networks
.
J. Hydraul. Eng.
125
(
7
),
676
685
.
Saemi
S.
,
Raisee
M.
,
Cervantes
M. J.
&
Nourbakhsh
A.
2019
Computation of two- and three-dimensional water hammer flows
.
J. Hydraul. Res.
57
(
3
),
386
404
.
Silva-Araya
W. F.
&
Chaudhry
M. H.
1997
Computation of energy dissipation in transient flow
.
J. Hydraul. Eng.
123
(
2
),
108
115
.
Urbanowiez
K.
2018
Fast and accurate modelling of frictional transient pipe flow
.
Z. Angew. Math. Mech.
98
(
5
),
802
823
.
Urbanowiez
K.
&
Zarzycki
Z.
2012
New efficient approximation of weighting functions for simulations of unsteady friction losses in liquid pipe flow
.
J. Theor. Appl. Mech.
50
(
2
),
487
508
.
Urbanowiez
K.
,
Tijsseling
A. S.
&
Firkowski
M.
2016
Comparing convolution-integral models with analytical pipe-flow solutions
.
J. Phys.: Conf. Ser.
760
,
012036
.
Vardy
A. E.
&
Brown
J. M. B.
1995
Transient, turbulent, smooth pipe friction
.
J. Hydraul. Res.
33
(
4
),
435
456
.
Vardy
A. E.
&
Brown
J. M. B.
2003
Transient turbulent friction in smooth pipe flows
.
J. Sound Vib.
259
(
5
),
1011
1036
.
Vardy
A. E.
&
Brown
J. M. B.
2004a
Efficient approximation of unsteady friction weighting functions
.
J. Hydraul. Eng.
130
(
11
),
1097
1107
.
Vardy
A. E.
&
Brown
J. M. B.
2004b
Transient turbulent friction in fully rough pipe flows
.
J. Sound Vibr.
270
(
1–2
),
233
257
.
Vardy
A. E.
&
Hwang
K.
1991
A characteristics model of transient friction
.
J. Hydr. Res.
29
(
5
),
669
684
.
Visual Studio [Computer software]. Visual Studio 2010
.
Vitkovsky
J. P.
,
Stephens
M.
,
Bergant
A.
,
Martin
L.
&
Simpson
A.
2004
Efficient and accurate calculation of Zielke and Vardy-Brown unsteady friciton in pipe transients
. In
Proc., 9th Int. Conf. on Pressure Surges
.
BHR Group
,
Cranfield
,
UK
, pp.
405
419
.
Wang
H.
,
Zhou
L.
,
Liu
D. Y.
,
Karney
B.
,
Wang
P.
,
Xia
L.
,
Ma
J. J.
&
Xu
C.
2016
CFD approach for column separation in water pipelines
.
J. Hydraul. Eng.
142
(
10
),
04016036
.
Wylie
E. B.
,
Streeter
V. L.
&
Suo
L.
1993
Fluid Transient in Systems
.
Prentice Hall
,
Englewood Cliffs, NJ
.
Zhou
L.
,
Wang
H.
,
Karney
B.
,
Liu
D. Y.
,
Wang
P.
&
Guo
S.
2018
Dynamic behavior of entrapped air pocket in a water filling pipeline
.
J. Hydraul. Eng.
144
(
8
),
04018045
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).