Abstract

Transient flow characteristics and dissipation mechanism in pressurized pipeline were investigated based on 1D friction models and 3D turbulence models, where the pressure–density model was combined into the 3D continuity equation allowing for the elasticity of the fluid and the pipes. The applicability of 3D realizable kε and 3D SST (shear stress transport) kω turbulence models was verified with comparison to 1D traditional water hammer models and the experimental data for fast closing of the valve in the reservoir–pipe–valve system. The valve closure rule was instantaneously carried out using the grid slip CFD (computational fluid dynamics) technique. The SST kω turbulence model has the highest accuracy in predicting the pressure attenuation of transient flows. The 3D detailed flow field confirms that the asymmetric flows induced by the change of valve opening within approximately three-fourths of the pipe inner diameter before the valve are captured. In the pressure wave cycles, the unsteady inertia, axial pressure gradient, viscous shear stress and turbulent shear stress mainly influence the velocity variations. During the pressure wave propagation, the viscous and turbulent dissipation are critical in the pressure attenuation in the wall region; the viscous dissipation is mainly concentrated in the viscous sublayer, while the turbulent dissipation increases to the maximum values at y+ = 13–23.

HIGHLIGHTS

  • The 3D weakly compressible and dissipation mechanism models are introduced to investigate transient flow characteristics.

  • The asymmetric flows induced by the change of valve opening are captured well within approximately three-fourths of the pipe inner diameter before the valve.

  • The viscous dissipation is mainly concentrated in the viscous sublayer, while the turbulent dissipation has maximum influence at y+ = 13–23.

INTRODUCTION

Hydraulic transients in pressurized pipeline may result in substantial pressure variations of flows due to the sudden change of fluid velocity. The propagation and reflection of pressure waves in the pipeline is the main cause of damages to hydraulic systems, such as device failures, leakages and pipe ruptures (Chaudhry 2014). The pressure change and flow characteristics in the pipeline were predicted to optimize the design of a hydraulic system, enhance protection and thus reduce the harm of water hammer events (Wang & Yang 2015). So far, a great deal of successful investigations were obtained in various analytical, experimental and numerical approaches.

The numerical prediction of hydraulic transients has been developed, and many models of transient flows have been proposed (Simpson 1986; Bergant et al. 2006). The 1D method of characteristics is widely used due to the fast and high-efficient prediction of transient flows. It solves the continuity and momentum equations governing the 1D transient flows. There are mainly two kinds of transient flow models considering friction resistance: the quasi-steady model (Ghidaoui et al. 2002) and the unsteady friction model (Vítkovsky et al. 2006). The 1D quasi-steady friction (1D-SF) model can precisely predict the first pressure peak of the water hammer period, but the accuracy of prediction decreases gradually as time evolves. The main reason is that the 1D-SF solver is based on the eddy viscosity distribution assumption of quasi-steady flow, which underestimates the observed energy dissipation of transient flows (Ferràs et al. 2016). To solve this problem, the 1D unsteady friction (1D-UF) models proposed by Vardy & Hwang (1993) and Brunone et al. (2004) introduced an independent unsteady friction to approximate the instantaneous energy dissipation, where the extra friction losses are assumed to be the function of the fluid local instantaneous acceleration. The amplitude and phase of pressure wave predicted by the 1D-UF model are in better agreement with the experimental data in the whole process of transient flows, and the energy dissipation can be better predicted by the 1D-UF model than the 1D-SF model. The 1D-UF models can effectively calculate the pressure and velocity variation during the transient process with relatively high accuracy, but it cannot provide detailed flow field information.

The quasi-two-dimensional (2D) unsteady friction model, which takes the velocity distribution of the cross section into consideration, was developed to obtain more information on the flow characteristics during hydraulic transients (Jang et al. 2016; Martins et al. 2017). Gao et al. (2018) investigated transient flows in the pipeline system with the 2D axial-symmetric model, and the MacCormack scheme was used to solve the governing equations, while the five region turbulence model was utilized for predicting turbulent characteristics of the transient flow in the study. The 2D model can better predict the magnitude and phase of pressure wave than the 1D model in most cases, but it is based on the axisymmetric flow velocity distribution assumption, which neglects the radial velocity component and its derivatives in the axis momentum equation. By using the 3D CFD (computational fluid dynamics) method, the velocity components in all three directions in the pipeline can be considered, and the transient characteristics of the pipeline system as well as the detailed flow field information can be obtained. So, the 3D CFD model is recommended to study the global flow field information, including the asymmetrical flow characteristics in the pipeline. But the 3D CFD method also has some drawbacks. For example, the 3D CFD method consumes more computing resources than the 1D and 2D methods; moreover, the accuracy of the 3D CFD simulation has great independence on the choice of turbulence models. Generally, the cost and run time of the 1D models are much smaller than that of the 3D CFD models, but the accuracy of the 3D CFD models is not totally higher than that of the 1D models in the pressure variation prediction, and it depends on the choice of the specific 1D friction or 3D turbulence models (Vítkovsky et al. 2006; Hooff et al. 2017).

The traditional 3D CFD method obtains the flow field information by solving the Navier–Stokes equation, in which the fluids are generally considered to be incompressible (Zhang et al. 2016; Li et al. 2017). However, when this method is used in calculating the direct water hammer events induced by the rapid valve closure, the propagation process of pressure wave in the pipeline cannot be captured. The weakly compressible model proposed by Song & Yuan (1988) can be used to solve this problem. A further study was also conducted to simulate transient flows based on the weakly compressible model. Martins et al. (2016) conducted 3D CFD simulation of hydraulic transient flows in pressurized pipes using the weakly compressible model, and the pressure wave propagation and the velocity profiles were analyzed, but the outlet boundary condition was set as flow rate varying with time, which is not accurate enough because there may be errors during the process to convert the pressure-head values into the discharge–time function, and there exists discrepancies between the fitted curve and the experimental flow rate curve. A more accurate method is to model the valve and consider its motion during the valve closure (Tang et al. 2020). Saem et al. (2019) carried out 2D and 3D simulations for transient flows, and valve modeling was considered in 3D simulation. The unsteady inertia, the pressure gradient and the viscous shear stress were found to be in balance in the near-wall region when the pressure wave passes. The 1D–3D coupling method is also used by some researchers to simulate transient flows. Wang et al. (2016) calculated the transient flow in a reservoir–pipe–valve system using the 1D–3D coupling method, and the position of coupling interface between the 1D and 3D regions is recommended in the study. Li et al. (2019) implemented the 1D–3D coupling method to predict both the transient process in the pipeline system and turbulent flows inside the pump, and the relationship between the external characteristics and the internal flow field in the pump is investigated. But different time steps and the data transfer between different computing domains may cause numerical errors.

In 3D CFD simulation, the Reynolds-averaged Navier–Stokes (RANS) method is often used to obtain the turbulent characteristics of unsteady transient flows, and it is very important to select a suitable turbulence model to accurately predict the energy dissipation when using RANS methods. A series of turbulence models were proposed in order to close the RANS equations, such as the standard kε model, RNG (Re-Normalisation Group) kε model, realizable kε model and SST (shear stress transport) kω model. The realizable kε model differs from the standard kε and RNG kε models, and it contains an alternative formulation for the turbulent kinetic energy k and a modified transport equation for the dissipation rate ɛ, resulting in the satisfaction of certain mathematical constraints (i.e. the positive of normal Reynolds stresses and Schwarz's inequality for turbulent shear stress) on the Reynolds stresses (Tsan et al. 1995), consistent with the physics of turbulent flows, and the three mentioned kε turbulence models here are only valid for fully developed turbulence, so they can only be used to solve the flow in the turbulent core region. However, a low Reynolds number turbulence model is needed for the flow in the near-wall region of pipes, and the low Reynolds SST kω turbulence model was then chosen for the simulation of the turbulent water hammer predictions. It employs the original kω model of Wilcox in the near-wall region and the standard kε model in the turbulent core region (Menter 1994).

The numerical hydraulic transient flows in the reservoir–pipe–valve system are investigated comprehensively and systematically based on 1D traditional and 3D CFD models in this study. The compressibility of fluid in the transient process is considered by introducing a weakly compressible model in 3D CFD simulation, and the pressure wave propagation is revealed by this model. The valve is modeled in 3D CFD simulation, and the instantaneous valve closure rule is carried out using the grid slip technique instead of the flow rate varying with the time outlet boundary condition. First, the numerical results predicted separately by 1D traditional and 3D CFD mathematical models are compared to a series of experimental data, respectively. Second, the comparisons are made between the numerical results by the 3D CFD weakly compressible model and the 1D traditional model. Lastly, the further comparisons are made between different 3D CFD turbulence models, and the applicability of the turbulence models is verified. Pressure and velocity evolutions in the axial section of the pipeline are analyzed in an entire water hammer cycle. Velocity profiles in several axial and radial positions are selected to study the asymmetrical flow characteristics induced by the valve closure, and the range of asymmetric flow is further captured. The cross-sectional velocity distributions during the pressure wave propagation are investigated, and the mechanism of velocity change is illustrated by analyzing the various terms in the axial momentum equation. The viscous and turbulent dissipation in the near-wall region of the pipe are revealed, and the applicability of the turbulence model for energy dissipation prediction is verified.

MATHEMATICAL MODEL AND NUMERICAL SCHEME

The 1D water hammer model

The traditional 1D-SF model uses the eddy viscosity distribution assumption of steady flows. While the 1D-UF model introduces an independent unsteady friction to approximate the energy dissipation and improve the accuracy of prediction. The method of characteristics can be utilized in both the models to solve the water hammer equations.

1D water hammer equations

The 1D water hammer equations comprise the momentum and continuity equations governing the transient flows (Chaudhry 2014).
formula
(1)
formula
(2)
where V is the mean velocity, H is the pressure head, t is the time, x is the distance along the pipe, g is the gravitational acceleration, a is the wave speed, α is the inclination angle of pipeline and J is the hydraulic losses.
In the 1D-SF model, only the steady friction Jq is considered in Equation (3):
formula
(3)
While in the 1D-UF model, both the steady friction Jq and the unsteady friction Ju are included in the hydraulic losses J (Vardy & Hwang 1993; Brunone et al. 2004):
formula
(4)
in which Ju is calculated as:
formula
(5)
where k′ is Brunone's friction coefficient, and Sign is the sign function as follows:
formula
(6)
Equations (1) and (2) constitute the integral equation of both quasi-steady and unsteady 1D water hammer flows.

Method of characteristics

The partial differential equations of the 1D-SF model were transformed into the ordinary differential equations as Equations (7) and (8) based on the method of characteristics (Pezzinga 2000).
formula
(7)
formula
(8)
The corresponding ordinary differential equations of the 1D-UF model were expressed as Equations (9) and (10):
formula
(9)
formula
(10)

Equations (7)–(10) were transformed into difference equations and solved to obtain the variation of pressure head H and flow rate Q with time in 1D transient flows.

3D CFD simulation

Governing equations of the weakly compressible model

The fluid was usually treated as incompressible in traditional CFD simulation, which is not applicable for transient flows, where the pressure wave propagation exists. However, if the solution is directly based on the compressible method, it is difficult to converge due to the huge difference between speed of sound and convection speed. So a weakly compressible model proposed by Song & Yuan (1988) was introduced in this paper.

The RANS method was utilized in the present study to predict the turbulent characteristics of unsteady transient flows. The RANS equations of continuity and momentum are as follows:

Continuity:
formula
(11)
Momentum:
formula
(12)
where ui is the fluid velocity in the xi spatial coordinate (with i = 1, 2, 3), p is the pressure, the overbar indicates the time average components, ρ is the density of the fluid, μ is the dynamic viscosity of the fluid, fi is the body force and τij denotes the Reynolds stress tensor:
formula
(13)
To consider the compressibility effects of transient flows, the pressure—density equation is added into the continuity equation by coupling density and pressure of fluids:
formula
(14)
where a is the pressure wave velocity, computed as follows (Bergant et al. 2006):
formula
(15)
in which K is the bulk modulus of elasticity of the fluid, E is the pipe Young modulus of elasticity, e is the pipe wall thickness and D is the inner diameter of the pipe. Elasticity of fluids and pipes are taken into consideration using Equation (15), resulting in a pressure variation more identical to the actual engineering.

Turbulence models

In the current study, the realizable kε and SST kω models were used as comparisons to verify the applicability of different turbulent models in solving transient flows in the pipeline. The transport equations of the turbulence kinetic energy k and the dissipation rate ε in the realizable kε model are as follows:
formula
(16)
formula
(17)
where . Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients, calculated as described in Modeling Turbulent Production in the kε models. Gb is the generation of turbulence kinetic energy due to buoyancy, calculated as described in Effects of Buoyancy on Turbulence in the kε models. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, calculated as described in Effects of Compressibility on Turbulence in the kε models. C2 = 1.9 and C1ε = 1.44 are the constants. σk = 1.0 and σε = 1.2 are the turbulent Prandtl numbers for k and ε, respectively. Sk and Sε are user-defined source terms.
The transport equations of the turbulence kinetic energy k and the specific dissipation rate ω in the SST kω turbulence model are as follows:
formula
(18)
formula
(19)
where represents the generation of turbulence kinetic energy due to mean velocity gradients, represents the generation of ω, and represent the effective diffusivity of k and ω, respectively, and represent the dissipation of k and ω due to turbulence, represents the cross-diffusion term, and and are user-defined source terms.

Computational techniques and methods

Computational domains and the grid model

The experimental data from Martins et al. (2016) are used to verify the accuracy of the above-mentioned numerical methods. As shown in Figure 1, the test system includes an upstream reservoir with constant head, a copper pipe and a ball valve located at the end section of the pipe. Transient flows in the pipe were induced by rapid closing of the valve, and the valve is closed from fully open to fully close within 14 ms. Two pressure transducers located at immediately upstream of the valve (PT1) and at the midsection of the entire pipe (PT2) were used to measure the pressure variation with a data acquisition frequency of 3 kHz. The X, Y and Z coordinates in Figure 1 represent the axial, vertical and horizontal directions of the flow, respectively.

Figure 1

Schematic of geometric models (Martins et al. 2016).

Figure 1

Schematic of geometric models (Martins et al. 2016).

The main parameters of pipe are as follows: pipe length L = 15.22 m, pipe thickness e = 1 mm and pipe inner diameter D = 20 mm. The wave speed in the pipeline is 1,250 m/s. Three test cases involving different initial conditions are calculated and analyzed in the paper. The parameters are presented in Table 1.

Table 1

Parameters of the test cases

Case IDInitial velocity, U0 (m/s)Re0 (dimensionless)Water level of the upstream reservoir, H0 (m)
Case 1 0.228 4,540 46.14 
Case 2 0.306 6,089 45.36 
Case 3 0.384 7,638 44.66 
Case IDInitial velocity, U0 (m/s)Re0 (dimensionless)Water level of the upstream reservoir, H0 (m)
Case 1 0.228 4,540 46.14 
Case 2 0.306 6,089 45.36 
Case 3 0.384 7,638 44.66 

The computational domains are established based on the parameters given here. During the calculation, the S1 monitor point is set 0.02 m upstream of the valve to monitor the pressure variation of PT1, and the S2 monitor point is set at the midsection of the pipe to record pressure variation of PT2.

Computational scheme and boundary conditions

The 1D water hammer equations are transformed into difference equations and solved by utilizing C# program. The variations of fluid pressure head H and flow rate Q with time in the 1D transient flow process are obtained.

Meanwhile, the 3D CFD simulations are conducted using the FLUENT solver of ANSYS platform. The pressure-based solver is used for the simulation of weakly compressible fluids, and the SIMPLE algorithm is chosen for pressure–velocity coupling in the flow field. The spatial discretization of pressure is a standard scheme (interpolating the pressure values at the faces using momentum equation coefficients; Rhie 1983), all the convective terms in the governing equations are discretized by the second-order upwind scheme, and the diffusion terms in the governing equations are discretized by the central differencing scheme. In addition, the first-order implicit scheme is employed for the discretization of the time derivatives. The reservoir was not modeled, but the constant pressure head provided by the reservoir is simulated by setting a constant pressure boundary condition at the inlet of pipeline. The slip grid technique is used to simulate the valve closure.

Grid independence verification

The whole flow domains are modeled using hexahedral structured meshes, and grid distributions in the cross section of pipes are shown in Figure 2. The grid distribution of the realizable kε and SST kω models is the same in the axial direction of the pipe and is different in the cross section of the pipe. In the turbulent core region, the grid distribution is the same, but in the near-wall region, the grid distribution for the SST kω model is more concentrated than that for the realizable kε model. The distance of the first grid node to the wall for the realizable kε model is set to be 0.25 mm in order to make sure that the corresponding y+ value is around 11.63. While the distance of the first grid node to the wall for the SST kω model is refined to be 0.064 mm to make sure that the corresponding y+ value is less than 1. The independence verification is conducted for both the realizable kε and SST kω models, and the results are presented here only for the SST kω model because they provide similar effect. Three sets of grids used in this study are with grid numbers of 2188442 (coarse-scale mesh), 2681578 (medium-scale mesh) and 4163058 (fine-scale mesh).

Figure 2

Grid distributions in the cross section of the pipeline: (a) the realizable kε model and (b) the SST kω model.

Figure 2

Grid distributions in the cross section of the pipeline: (a) the realizable kε model and (b) the SST kω model.

The pressure variations at the S1 monitor point for Case 1 are calculated using the three sets of grids, and the results are compared with the experimental data. The predicted results match well with the experimental data, and with the increase of grid number, the difference between the predicted data and the experimental results decreases and remains almost unchanged (around 1% by the coarse-scale mesh, less than 0.5% by the medium-scale mesh and fine-scale mesh). The calculations were performed using the 12 cores 3.7 GHz Intel Xeon CPU. The time consumption of the coarse-, medium- and fine-scale meshes is around 12, 16 and 36 h, respectively. Therefore, from the perspective of computational accuracy and computing resources, the medium-scale mesh (grid number 2681578) is selected for the CFD calculation. The 1D–3D coupling method will be used by the authors in the future, and the efficiency and accuracy of this method will be compared to that of the 1D water hammer and 3D CFD models, respectively.

RESULTS AND DISCUSSIONS

The pressure propagation and flow field characteristics in the pipeline during the hydraulic transient process will be analyzed and compared by different 1D and 3D water hammer models. The accuracy of all the predictions will be further compared and discussed.

Analysis of predicted results by 1D and 3D models

Pressure variations obtained from different models against the experimental data for Case 1 are shown in Figure 3. The 1D traditional models used here include the quasi-steady friction model (1D-SF), the unsteady friction model (1D-UF), the 3D weakly compressible models that allow for the realizable kε turbulence model (3D kε) and the SST kω turbulence model (3D kω). The pressure predicted by 3D models is obtained based on the cross-sectional area-weighted average method.

Figure 3

Pressure variations versus time at PT1 obtained by different models for Case 1: (a) pressure variation during 0–0.25s and (b) zoomed sections of pressure variations in (a) during 0–0.05 s.

Figure 3

Pressure variations versus time at PT1 obtained by different models for Case 1: (a) pressure variation during 0–0.25s and (b) zoomed sections of pressure variations in (a) during 0–0.05 s.

It can be seen from Figure 3 that in the first period of pressure wave evolutions, the predicted amplitudes and phases by all models are close to the experimental data except for the 1D-SF model. The differences of the first pressure peak values between the numerical results by all the models and the experimental data are no more than 1%. The 1D-SF model only predicts the first pressure peak value precisely, and the relevant difference of pressure peak values between the numerical and the experimental data increases after the first pressure peak. In contrast, since the 1D-UF model considers the local instantaneous acceleration, its prediction of pressure attenuation is more accurate than that of the 1D-SF model. The 3D realizable kε turbulence model predicts the pressure attenuation better than the 1D-SF model but worse than the 1D-UF model. The realizable kε turbulence model uses the normal wall functions in the near-wall region, and it cannot capture the energy dissipation well, so its prediction of pressure attenuation is relatively inaccurate. The 3D SST kω turbulence model utilizes the original kω model with high grid resolution solving the flow near the wall (Menter 1994), and it can capture the energy dissipation more accurately than the realizable kε turbulence model, so it is adopted. It shows in Figure 3 that the predicted pressure variations by the 3D SST kω turbulence model are the closest to the experimental data in all used models, and the first and second peak values of pressure head are 74.5 and 72.7 m, respectively. The pressure peak relative errors around 0.22 s of four models are 14.7, 9.4, 4.4 and 4.1% for the 1D-SF model, the 3D Realizable kε model, the 1D-UF model and the 3D SST kω model, respectively.

Pressure variations predicted by the 3D realizable kε turbulence and 3D SST kω turbulence models for Cases 2 and 3 are shown in Figure 4. The predicted pressure variations by the 3D SST kω turbulence model are still closer to the experimental data than the 3D realizable kε turbulence model. Based on above analysis, the 3D SST kω turbulence model was used to study the flow characteristics in the pipeline. The transient flows under higher initial Reynolds number conditions (especially cavitation flows) will be further studied in the future works.

Figure 4

Pressure variations versus time at PT1 obtained from different 3D turbulence models for Case 2 (a) and Case 3 (b).

Figure 4

Pressure variations versus time at PT1 obtained from different 3D turbulence models for Case 2 (a) and Case 3 (b).

Analysis of the flow characteristics in the pipeline

Pressure and velocity evolutions in the axial section of the pipeline

The pressure and velocity evolutions for Case 1 in axial and vertical sections of the pipe in a whole period of water hammer transient flow are shown in Figures 5 and 7, respectively (T is the period of the wave cycle). The whole pressure wave propagation process is divided into four stages. The calculated pressure wave speed is 1,247.5 m/s, with respect to the experimental data 1,250 m/s.

Figure 5

Pressure evolutions during one entire water hammer period.

Figure 5

Pressure evolutions during one entire water hammer period.

In the first stage, the pressure increase wave propagates from the valve to the upstream reservoir. At T0, the pressure and velocity fields in the pipeline are in the steady-state condition. At T1, the pressure increase wave induced by the valve closure travels upstream, and the velocity near the valve decreases to zero. The velocity in the near-wall region first decreases to zero due to the small inertia, and the velocity profile at the core region tends to stay the original shape due to the large inertia. At T2, the pressure wave reaches the midpoint of the pipeline, and the velocity from the valve to the midsection reduces to zero. At T3, the pressure wave reaches the upstream reservoir, the velocity in the entire pipe is almost zero and the reversed velocity appears near the upstream reservoir. The pressure distribution at cross sections of the pipe is almost uniform at any moment of the transient process. As shown in Figure 6, the area-weighted average pressure values at the pipe midsection are 462,232, 712,216 and 745,058 Pa at moments T1, T2 and T3, respectively.

Figure 6

Cross-sectional pressure distribution at the midsection of the pipeline: (a) T1 = 1/8T, (b) T2 = 2/8T, and (c) T3 = 3/8T.

Figure 6

Cross-sectional pressure distribution at the midsection of the pipeline: (a) T1 = 1/8T, (b) T2 = 2/8T, and (c) T3 = 3/8T.

In the second stage, the pressure decrease wave travels back from the upstream reservoir to the valve. At T4, the pressure wave and the reversed velocity reach the midpoint of the pipe and, at T5, arrive at the downstream valve.

In the third stage, the pressure decrease wave propagates from the valve to the upstream reservoir. At T6, the pressure wave reaches the midsection of the pipe, and the reversed velocity diminishes from the valve to the midpoint of the pipe. At T7, the pressure wave arrives at the valve, the velocity restores to zero in the entire pipe and positive velocity appears near the reservoir.

In the fourth stage, the pressure increase wave propagates from the upstream reservoir to the valve and restores the pressure and velocity fields in the pipeline to almost the initial steady state.

It also shows in Figure 5 and Figure 7 that the pressure contour is perpendicular to the pipe axis, while the velocity contour is not uniform in the pipeline. The pressure wave propagates continuously and periodically in the pipeline and weakens gradually due to the frictional resistance.

Figure 7

Velocity distributions during one entire water hammer cycle.

Figure 7

Velocity distributions during one entire water hammer cycle.

Asymmetrical characteristics of 3D transient flow

In steady-state conditions, the flow in pipes can be regarded as axial symmetry; however, transient flows caused by the change of valve opening may result in asymmetrical flow characteristics near the valve. Five positions in the axial direction, namely 0.25D, 0.5D, 0.75D, 1D and 1.25D upstream of the valve as shown in Figure 8(a) and three circumferential positions α = 0, 45 and 90° from the vertical plane to the horizontal plane as shown in Figure 8(b), are selected to investigate the velocity profiles (Figure 9) at different moments t = T, 2 and 3T (R0 is the pipe inner radius and R/R0 is a dimensionless parameter, where R is positive for y > 0 and negative for y < 0).

Figure 8

Velocity profiles at different locations from the valve: (a) the axial direction and (b) the radial direction.

Figure 8

Velocity profiles at different locations from the valve: (a) the axial direction and (b) the radial direction.

Figure 9

Axial velocity profiles in different radial and axial directions at different moments.

Figure 9

Axial velocity profiles in different radial and axial directions at different moments.

As shown in Figure 9(a), at α = 0°and t = T, the velocity profiles at 0.25D and 0.5D demonstrate apparent asymmetry; When R < 0, at 0.25D, the maximum adverse velocity (−0.5U0) occurs at −0.8R and the velocity profile is smooth, while at 0.5D, the maximum adverse velocity (−1.1U0) is almost near the wall and the velocity distribution changes rapidly. At t = 2T, the velocity profiles only demonstrate complex asymmetry characteristics at 0.25D, and when R < 0, there are two adverse velocity peaks: −0.7U0 at −0.4R0 and −0.9U0 at −0.9R0. At t = 3T, the maximum adverse velocity for 0.25D and 0.5D occurs at −0.6R0 and −0.9R0, respectively, with the different shapes of variation.

As shown in Figure 9(b), at α = 45° and t = T, the velocity profile at 0.25D demonstrates asymmetric, and there is an offset in the velocity profile at 0 to −0.4R0. At t = 2T and 3T, the velocity profiles are symmetric at all axial positions except for 0.25D.

As demonstrated in Figure 9(c), at α = 90° and t = T, the entire velocity profiles are symmetric in all axial positions, and at 0.25D position, there exists slight offset at −0.2 to 0.2R0. At t = 2T and 3T, all the velocity profiles are almost symmetric.

The asymmetric flows induced by the valve closure are well captured within approximately 3/4D upstream of the valve, and it is most obvious at the vertical plane. The main reason is that with the reduction of valve opening, the flow area of the valve decreases, which leads to the fluid velocity at the pipeline's lower part (R < 0) higher than that at the pipeline's upper part (R > 0), and finally leads to the asymmetric characteristics of the flow field near the valve. But the asymmetric flow at the valve does not affect the symmetric distribution of pressure in the cross section during pressure wave propagation.

Analysis of cross-sectional velocity distribution in the pipeline

Velocity profiles at the midsection of the pipe are analyzed at different moments in Figure 10. The whole process can be divided into five stages based on the velocity changes.

Figure 10

Pressure-head evolutions at the midsection of the pipeline.

Figure 10

Pressure-head evolutions at the midsection of the pipeline.

As shown in Figure 11(a), at t = 0.0 s, the velocity in the midsection of the pipeline is in a fully developed distribution before the valve is closed, with the mean velocity U0. At t = 0.014 s, the velocity profile decreases entirely due to the passing of pressure increase wave, which is induced by the closure of the downstream valve as shown in Figure 11(b). At t = 0.019 s, the pressure reaches its maximum value and maintains in a short period after the pressure increase wave passes, and in the near-wall region, a significant reverse flow occurs, with a large velocity gradient. As shown in Figure 12, during the second stage of water hammer cycle, as the pressure decrease wave propagates from the upstream reservoir to the downstream valve, the adverse velocity increases gradually, while the profiles keep the same shape in the core region as the steady state. Velocity variations during the third stage are shown in Figure 13, as the pressure decrease wave propagates upstream, the mean velocity decreases to zero; the profile in Figure 13(b) is similar to that in Figure 11(c). As it can be seen from Figure 14, in the fourth stage, as the pressure increase wave propagates downstream, the velocity restores gradually to the steady state, and the velocity profile has a spike in the near-wall region. The fine grid distributions and the SST kω turbulence model used are efficient to capture detailed transient flows in the near-wall regions. Figure 15 shows the pressure increase wave propagation in the early stage of the second water hammer cycle. The performance is almost the same as that of the first water hammer cycle, and a slight velocity variation is captured in the near-wall region compared to Figure 11. The pressure wave propagation affects the entire velocity distribution in the cross section of the pipeline, while the unsteady inertia and viscous shear stress deicide the velocity gradient between the near-wall region and the turbulent core region. When the pressure wave passes, the velocity in the near-wall region first decreases to zero due to the small inertia, and after the pressure wave passes, there is a high gradient due to the viscous effect. But in the core region, the velocity profile tends to stay the original shape due to the large inertia.

Figure 11

Velocity profile variations during the first stage of the water hammer period: (a) t = 0.0s, (b) t = 0.014s, and (c) t = 0.019s.

Figure 11

Velocity profile variations during the first stage of the water hammer period: (a) t = 0.0s, (b) t = 0.014s, and (c) t = 0.019s.

Figure 12

Velocity profile variations during the second stage of the water hammer cycle: (a) t = 0.27s, (b) t = 0.030s, and (c) t = 0.033s.

Figure 12

Velocity profile variations during the second stage of the water hammer cycle: (a) t = 0.27s, (b) t = 0.030s, and (c) t = 0.033s.

Figure 13

Velocity profile variations during the third stage of the water hammer cycle: (a) t = 0.40s and (b) t = 0.044s.

Figure 13

Velocity profile variations during the third stage of the water hammer cycle: (a) t = 0.40s and (b) t = 0.044s.

Figure 14

Velocity profile variations during the fourth stage of the water hammer cycle: (a) t = 0.51s, (b) t = 0.055s, and (c) t = 0.057s.

Figure 14

Velocity profile variations during the fourth stage of the water hammer cycle: (a) t = 0.51s, (b) t = 0.055s, and (c) t = 0.057s.

Figure 15

Velocity profile variations during the fifth stage of the water hammer cycle: (a) t = 0.064s and (b) t = 0.068s.

Figure 15

Velocity profile variations during the fifth stage of the water hammer cycle: (a) t = 0.064s and (b) t = 0.068s.

Mechanism of velocity change during the water hammer cycle

To study the mechanism of velocity changes during the water hammer cycle in cross sections of pipes, variations of different terms in the x-momentum (Equation (20)) in the cylindrical coordinate during the water hammer event are analyzed here. Terms 1, 2, 3 and 4 represent unsteady inertia, axial pressure gradient, viscous shear stress gradient and turbulent shear stress gradient in Equation (20), respectively, while the other three terms are smaller in the order of magnitude than these four terms based on the calculated results.
formula
(20)
As it can be seen from Figure 16, the unsteady inertia term is zero at the initial steady state. The axial pressure gradient, the viscous shear stress gradient and the turbulent shear stress gradient are in fully developed balance state, where the sum of those three terms is almost zero. And the axial pressure gradient maintains constant from pipe core to pipe wall.
Figure 16

Variations of different gradient terms along the radial direction at t = 0.0s.

Figure 16

Variations of different gradient terms along the radial direction at t = 0.0s.

The pressure wave passes through the midsection of the pipeline at t = 0.014 s. As shown in Figure 17, in the core region, the axial pressure gradient (50,354.0 Pa/m) and the unsteady inertia (−53,955.2 Pa/m) dominate the flow and make the velocity reduce, while the absolute values of the other two terms are less than 61.0. In the near-wall region, the four terms keep balance at y+ ≈ 9 (r ≈ 0.0095 m), and the viscous shear stress gradient and the turbulent shear stress gradient are −252.4 and 228.4 Pa/m respectively; in the region closer to the wall, the turbulent shear stress effect disappears gradually; meanwhile, the unsteady inertia reduces sharply, and the viscous shear stress gradient reaches a maximum value of 442.1 Pa/m at y+ ≈ 1.1 (r ≈ 0.00994 m). In the region where y+ < 1.1, the unsteady inertia, the viscous force, and the pressure are in dynamic balance.

Figure 17

Comparison of various terms in the axial momentum equation at the midsection of the pipe at time t = 0.014s: (a) overall display and (b) zoomed section of (a).

Figure 17

Comparison of various terms in the axial momentum equation at the midsection of the pipe at time t = 0.014s: (a) overall display and (b) zoomed section of (a).

At t = 0.019 s as shown in Figure 18, in the core region, the four terms are relatively small compared to that at t = 0.014 s. In the near-wall region, the viscous force (−6,024.7 Pa/m) and the unsteady inertia (5,834.1 Pa/m) dominate the flow at y+ ≈ 2.5 (r ≈ 0.00986 m), which corresponds to the high-velocity gradient zone in Figure 11(c).

Figure 18

Variations of different gradient terms along the radial direction at t = 0.019s.

Figure 18

Variations of different gradient terms along the radial direction at t = 0.019s.

Based on the above analysis, the turbulent shear stress has great influence on the near-wall region around y+ = 9. So, the turbulent and viscous dissipation are further analyzed. The dissipation function (White 1991) is introduced here to reveal the viscous and turbulent dissipation in pipeline flows:
formula
(21)
in which represents the viscous dissipation, represents the turbulent dissipation.

It can be seen from Figure 19 that in the viscous sublayer, the viscous dissipation is very high, while turbulent dissipation is particularly small. The values of viscous dissipation at the wall are 63.5, 479.9 and 1,491.8 Pa/s for t = 0.0, 0.014 (when the pressure wave passes) and 0.019 s (after the pressure wave passes), respectively. As y+ increases, the viscous dissipation reduces, while the turbulent dissipation increases and reaches the maximum values of 16.0, 15.4 and 16.0 Pa/s at y+ = 13, 23 and 13 for t = 0.0, 0.014 and 0.019 s. When approaching the core region, the viscous dissipation is almost zero, and the turbulent dissipation decreases gradually from the maximum value to zero. The maximum reverse velocity appears at y+ = 2.5 in Figure 11(c), and in the same region, the viscous shear stress reaches the maximum value and is in dynamic balance with the unsteady inertia (Figure 18); there is a high-velocity gradient near this region, and the viscous dissipation is also very large. At y+ = 37, the velocity is zero in Figure 11(c), and the turbulent dissipation appears mainly around y+ = 13–37. The energy dissipation is accurately captured by the SST kω model, so the pressure attenuation is predicted well.

Figure 19

Viscous and turbulent dissipation at the midsection of the pipe at time t = 0.0, 0.014, 0.019s: (a) overall display and (b) zoomed section of (a).

Figure 19

Viscous and turbulent dissipation at the midsection of the pipe at time t = 0.0, 0.014, 0.019s: (a) overall display and (b) zoomed section of (a).

CONCLUSIONS

In this paper, the 1D traditional and 3D weakly compressible models have been employed to investigate transient flows in pressurized pipeline, and the applicability of the models is verified by comparing the predicted results. The cost and run time of the 1D models (less than 30 s on single core 3.7 GHz Intel Xeon CPU) are much smaller than those of the 3D CFD models (more than 10 h on 12 cores 3.7 GHz Intel Xeon CPU). The 1D simulations have been implemented by considering different types of hydraulic losses in water hammer equations, while in the 3D simulations, various turbulence models have been used to take the turbulent characteristics of transient flows into consideration. The pressure variation predicted by different models and the multidimensional inner flow characteristics have been analyzed. The dissipation mechanism model is used to analyze the velocity distribution near the wall during the pressure wave propagation. The main conclusions are as follows:

  1. The 1D-SF model can only predict the first pressure peak value precisely, and the difference of pressure peak values between the numerical and the experimental data increases after the first pressure peak. The 1D-UF model predicts the pressure attenuation more accurately than those by the 1D-SF model. The predicted results of 1D model is used to verify the 3D model. The 3D realizable kε turbulence model predicts the pressure attenuation better than the 1D-SF model but worse than the 1D-UF model. The predicted pressure variations by the 3D SST kω turbulence model are the closest to the experimental data in all used models. The pressure peak relative errors around 0.22 s of the 3D SST kω model, the 1D-UF model, the 3D realizable kε model and the 1D-SF model are 4.1, 4.4, 9.4 and 14.7%, respectively.

  2. By using the weakly compressible model and the 3D SST kω turbulence model, the pressure propagation evolution is captured. The pressure contour is perpendicular to the pipe axis, while the velocity contour is not uniform during the pressure wave propagation. It was observed that the 3D asymmetric flows induced by the change of valve opening are within approximately three-fourths of the pipe inner diameter upstream of the valve.

  3. The velocity variations are mainly affected by the unsteady inertia, axial pressure gradient, viscous shear stress and turbulent shear stress. When the pressure wave passes, in the core region, the unsteady inertia and axial pressure gradient dominate the flow; in the region around y+ = 9, the four terms keep balance; in the region where y+ < 1.1, the unsteady inertia, the viscous force and the pressure are in balance. After the pressure wave passes, the unsteady inertia and viscous shear stress are the essential terms around y+ = 2.5.

  4. When the pressure wave passes the midsection of the pipe, the viscous dissipation is very high in the viscous sublayer, and it reaches a higher value after the pressure wave passes. As y+ increases, the viscous dissipation reduces, while the turbulent dissipation increases to the maximum value. When approaching the core region, the viscous dissipation is almost zero, and the turbulent dissipation decreases gradually from the maximum value to zero.

  5. After the pressure wave passes, the maximum reverse velocity appears at y+ = 2.5, and in the same region, the viscous shear stress reaches the maximum value and is in dynamic balance with the unsteady inertia. At y+ = 37, the velocity is zero, and the turbulent dissipation appears mainly around y+ = 13–37. The energy dissipation is accurately captured by the SST kω model, so the pressure attenuation is predicted well.

ACKNOWLEDGMENT

This work was supported by the National Natural Science Foundation of China (Grant Nos. 51779257, 52079140 and 51479196).

DATA AVAILABILITY STATEMENT

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

REFERENCES

REFERENCES
Bergant
A.
,
Simpson
A. R.
&
Tijsseling
A. S.
2006
Water hammer with column separation: a historical review
.
Journal of Fluids & Structures
22
(
2
),
135
171
.
Brunone
B.
,
Ferrante
M.
&
Cacciamani
M.
2004
Decay of pressure and energy dissipation in laminar transient flow
.
Journal of Fluids Engineering
126
(
6
),
928
934
.
Chaudhry
H. M.
2014
Applied Hydraulic Transients
.
Springer Verlag
,
New York, NY
,
USA
.
Ferràs
D.
,
Manso
P. A.
,
Schleiss
A. J.
&
Covas
D. I. C.
2016
Experimental distinction of damping mechanisms during hydraulic transients in pipe flow
.
Journal of Fluids and Structures
66
,
424
446
.
Gao
H.
,
Tang
X.
,
Li
X.
&
Shi
X.
2018
Analyses of 2D transient models for vaporous cavitating flows in reservoir-pipeline-valve systems
.
Journal of Hydroinformatics
20
(
3–4
),
934
945
.
Ghidaoui
M. S.
,
Mansour
S. G. S.
&
Zhao
M.
2002
Applicability of quasisteady and axisymmetric turbulence models in water hammer
.
Journal of Hydraulic Engineering
128
(
10
),
917
924
.
Jang
T. U.
,
Wu
Y.
,
Xu
Y.
,
Newman
J.
&
Sun
Q.
2016
Efficient quasi-two-dimensional water hammer model on a characteristic grid
.
Journal of Hydraulic Engineering
142
(
12
),
060160191
.
Li
Z.
,
Bi
H.
,
Karney
B.
,
Wang
Z.
&
Yao
Z.
2017
Three-dimensional transient simulation of a prototype pump-turbine during normal turbine shutdown
.
Journal of Hydraulic Research
55
(
4
),
520
537
.
Li
X.
,
Tang
X.
,
Zhu
M.
&
Shi
X.
2019
1D-3D coupling investigation of hydraulic transient for power-supply failure of centrifugal pump-pipe system
.
Journal of Hydroinformatics
21
(
5
),
708
726
.
Martins
N. M. C.
,
Soares
A. K.
,
Ramos
H. M.
&
Covas
D. I. C.
2016
CFD modeling of transient flow in pressurized pipes
.
Computers & Fluids
126
(
1
),
129
140
.
Martins
N. M. C.
,
Delgado
J. N.
,
Ramos
H. M.
&
Covas
D. I. C.
2017
Maximum transient pressures in a rapidly filling pipeline with entrapped air using a CFD model
.
Journal of Hydraulic Research
55
(
4
),
506
519
.
Pezzinga
G.
2000
Evaluation of unsteady flow resistances by quasi-2D or 1D models. Technical Note
.
Journal of Hydraulic Engineering
126
(
10
),
778
785
.
Saem
S. I.
,
Raisee
M.
,
Cervantes
M. J.
&
Nourbakhsh
A.
2019
Computation of two- and three-dimensional water hammer flows
.
Journal of Hydraulic Research
57
(
3
),
386
404
.
Simpson
A. R.
1986
Large Water Hammer Pressure due to Column Separation in Sloping Pipes
.
University of Michigan
,
Ann Arbor, MI
.
Song
C. C. S.
&
Yuan
M.
1988
A weakly compressible flow model and rapid convergence methods
.
Journal of Fluids Engineering
110
(
4
),
441
445
.
Tsan-Hsing
S.
,
William
W. L.
&
Aamir
S.
1995
A new k-ε eddy viscosity model for high Reynolds number turbulent flows: model development and validation
.
Computers & Fluids
24
(
3
),
227
238
.
Vardy
A. E.
&
Hwang
K. L.
1993
A weighting function model of transient turbulent pipe friction
.
Journal of Hydraulic Research
31
(
4
),
533
548
.
Vítkovsky
J. P.
,
Bergant
A.
,
Simpson
A. R.
&
Lambert
M. F.
2006
Systematic evaluation of one-dimensional unsteady friction models in simple pipelines
.
Journal of Hydraulic Engineering
132
(
7
),
696
708
.
Wang
C.
&
Yang
J. D.
2015
Water hammer simulation using explicit-implicit coupling methods
.
Journal of Hydraulic Engineering
141
(
4
),
04014086
.
Wang
C.
,
Nilsson
H.
,
Yang
J.
&
Petit
O.
2016
1D-3D coupling for hydraulic system transient simulations
.
Computer Physics Communications
210
,
1
9
.
White
F. M.
1991
Viscous Fluid Flow
.
McGraw-Hill
,
New York, NY
,
USA
.
Zhang
X.
,
Cheng
Y.
,
Xia
L.
&
Yang
J.
2016
CFD simulation of reverse water-hammer induced by collapse of draft-tube cavity in a model pump-turbine during runaway process
.
IOP Conference
49
(
5
),
052017
.