Investigations of hydraulic transient flows in pressurized pipeline based on 1D traditional and 3D weakly compressible models

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.

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.

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.

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:

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.

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.

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

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.

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

The pressure variations at the S₁ 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.

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.

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.

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.

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.

In the first stage, the pressure increase wave propagates from the valve to the upstream reservoir. At T₀, the pressure and velocity fields in the pipeline are in the steady-state condition. At T₁, 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 T₂, the pressure wave reaches the midpoint of the pipeline, and the velocity from the valve to the midsection reduces to zero. At T₃, 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 T₁, T₂ and T₃, respectively.

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

In the third stage, the pressure decrease wave propagates from the valve to the upstream reservoir. At T₆, the pressure wave reaches the midsection of the pipe, and the reversed velocity diminishes from the valve to the midpoint of the pipe. At T₇, 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.

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 (R₀ is the pipe inner radius and R/R₀ is a dimensionless parameter, where R is positive for y > 0 and negative for y < 0).

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.5U₀) occurs at −0.8R and the velocity profile is smooth, while at 0.5D, the maximum adverse velocity (−1.1U₀) 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.7U₀ at −0.4R₀ and −0.9U₀ at −0.9R₀. At t = 3T, the maximum adverse velocity for 0.25D and 0.5D occurs at −0.6R₀ and −0.9R₀, 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.4R₀. 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.2R₀. 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.

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.

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

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.

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

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.

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.

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

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