Unsteady friction model modified with compression–expansion effects in transient pipe flow

a

Wave speed

A

Cross-sectional area

D

Pipe diameter

f

Darcy–Weisbach friction factor

g

Gravitational acceleration

H

Pressure head

k

Empirical coefficient

L

Pipe length

P

Cross-sectional average pressure

Q

Volumetric flow rate

r

Radial coordinate

R

Pipe radius

Re

Reynolds number

V

Cross-sectional average velocity

t

Time

x

Coordinate along the pipe axis

ρ

Fluid density

μ

Dynamic viscosity

θ

Circumferential coordinate

τ

Wall friction

Hydraulic transient is an intermediate stage from one steady flow state to another. The transient pipe flow is such a critical flow phenomenon, which can be characterized by a series of positive and negative pressure waves moving along pipes at wave speed. It is also known as water hammer, which is usually caused by improper operation of the pump or rapid motion of the valve in the pipeline. Considering the fast propagation speed and varying amplitude of the pressure wave in unsteady pipe flow, the accurate modeling and prediction of the pressure wave is essential for pipeline design and defect detection (Colombo et al. 2009; Duan et al. 2017a, 2017b, 2020).

In order to simulate the pressure wave in transient pipe flow accurately, unsteady friction models were established, which contain the instantaneous acceleration-based (IAB) model and the convolution-based (CB) model. Daily et al. (1956) first proposed the IAB model, and they found that unsteady wall shear stress was positive in accelerated flow and negative in decelerated flow. With the assumption that the damping was attributable to unsteady friction caused by instantaneous accelerations, the accelerations were calculated from the average cross-sectional values without considering consideration of the influence of the velocity distribution at the cross-section in the IAB model. This modification method of the IAB model became the main research direction afterward. Axworthy et al. (2000) verified the rationality of the IAB model by theoretical analysis and then determined the applicable scope of the model. Brunone et al. (1995) established the single coefficient IAB model, which considered the wave velocity and convection acceleration to the unsteady friction based on the IAB model. This developed model was then proven to be compatible in obtaining results consistent with the experimental values (Bergant et al. 2001). However, the Brunone model neglected the symbol's convection acceleration term, which may lead to unreasonable results when convective acceleration is with the opposite sign of instantaneous acceleration. Considering the coordinate system is very important for the Brunone model, Pezzinga (2000) modified the Brunone model by adding a term estimating the velocity symbol in the equation, so that the applicability limitation and the influence of the coordinate system could be eliminated.

In the IAB model, local acceleration and convection acceleration have different physical meanings. Local acceleration can only affect the pressure velocity, while convection acceleration can only produce friction. In order to distinguish these two effects, Ramos et al. (2004) proposed a two-coefficient IAB model. According to Ramos' research, the phase of the pressure wave is determined by the local acceleration, and the amplitude of the pressure wave is determined by the convective acceleration. The selection of empirical parameters in the IAB model has a great influence on its calculation accuracy, which is generally determined by the trial-and-error method. In order to determine the empirical parameters more accurately, Vardy & Brown (2003) proposed the empirical calculation formula related to the Reynolds number based on the concept of global acceleration, and the Reynolds number could be determined according to the initial stable flow state or the instantaneous value, but the calculation accuracy of the latter is relatively high. For the two-coefficient IAB model, Storli & Nielsen (2011) proposed the empirical calculation formula that changes with the change of pipeline length. Compared with the constant empirical parameters, the results obtained from the proposed model are more consistent with the measured values. However, the empirical formula is only applicable to the simple water tank–pipeline–valve experimental system, and its applicability to other cases needs to be verified. The IAB model successfully connects the friction loss with the flow characteristics and only contains the first-order partial derivatives of the velocity, which is an important feature for efficient program-solving.

The convolution-based model is another unsteady friction model, which was proposed by Zielke (1968). This model was obtained by solving the Laplace transform of Navier–Stokes equations in the axial direction of the pipeline. Therefore, the CB model has a solid theoretical basis and, to some extent, reveals the essence of transient flow friction. The CB model requires a lot of computing time and storage space to calculate the convolution of local acceleration and weight function at all times of each node. Vardy & Brown (2004) introduced the Reynolds number to the weight function and divided the flow into two regions. One is the outer annular regions where the eddy viscosity was assumed to vary linearly with distance from the pipe wall. Another is the inner core regions where the eddy viscosity is assumed to be constant. On this basis, a weight function model suitable for hydraulic smooth and rough pipelines was then established. The Zarzycki model is similar to the Vardy model, and the difference lies in the fact that the Zarzycki model adopts the assumption of four-layer viscosity coefficient distribution and gives different Reynolds number formulas to distinguish laminar flow and turbulent flow (Zarzycki & Kudzma 2007). In stable flow, Meniconi et al. (2014) found that the proportion of unsteady friction loss to total friction loss in the CB model was inversely proportional to the Reynolds number. Hence, they modified the model with the instantaneous Reynolds number. This CB model has a strong theoretical basis, and the accuracy of this model has been verified by many experiments. However, the model includes the calculation of historical acceleration, and even the simplified model still depends on much computational power. Meanwhile, due to the absence of the convection term in the model, errors may easily accumulate in the calculation of transient flows with strong convection.

For the IAB model and the CB model, the unsteady friction term is usually constructed in the momentum equation based on the shear dissipation effect, and the calculation of amplitude and phase of the pressure wave has a good agreement at the initial stage of the fluctuation, but the error increases with the calculation time. In addition, as pressure also exists in laminar flow and static fluid, and the turbulent dissipation effect is not applicable, it is necessary to consider other energy dissipation mechanisms in the one-dimensional model apart from the shear dissipation effect and the turbulent dissipation effect.

According to Brunone & Golia (2006), different transition phases of accelerating and decelerating can be classified by the acceleration parameter. On this basis, a proposed velocity coefficient was applied to characterize the different behaviors of the velocity profile. In order to better investigate the unsteady energy dissipation behavior, Brunone & Berni (2010) experimentally measured the instantaneous velocity profiles through the ultrasound method. They observed different variations of the wall shear stress for the accelerating and decelerating flow phase, while the distinction was less obvious after wave propagation and reflection at a large time scale. Knowing that the unsteady friction component in the water hammer equation is closely related to fluid acceleration, Duan et al. (2012) further carried out a dimensional analysis and found that the unsteady friction effect is more important for small-scale water supply and transmission systems.

Apart from the analytical results from the one-dimensional model, Martins et al. (2017) studied the transient dissipation behavior by the CFD method using a realizable k–ε model. Waveform comparison suggests that the CB model better fits the CFD model as it considers the velocity time history in calculating the unsteady friction term. It should be noted that despite a higher numerical resolution in detailed flow field analysis, CFD models may still have limitations in terms of the neglected fluid-structure energy exchange, simplified laminar assumption of wall shear law, and obviously high computational cost, when compared with eligible IAB or CB models in practical application for large-scale hydraulic systems (Mandair et al. 2020; Marsili et al. 2022).

According to Riedelmeier et al. (2014), the pressure wave damping obtained by the three-dimensional model was higher, while the wall friction was small. This indicates that the factors affecting the pressure wave energy dissipation include shear dissipation and other effects. Dörfler (2011) and Landry et al. (2012) pointed out that the unsteady friction could be represented by a friction term dependent on local inertia and wall friction or by a dissipation term dependent on the expansion viscosity. Mitosek & Szymkiewicz (2012) found that the transient flow model with the diffusive term, which represents the friction generated by compression or expansion of the fluid, could reproduce the irreversible processes at the pressure wavefront. In elastic pipes, a similar attenuation of the pressure wave magnitude was observed and was attributed to the radial compression or expansion of the pipe wall (Chahardah-Cherik et al. 2021; Essaidi & Triki 2021).

Due to the high computational efficiency of the IAB model, the application of the IAB model is more beneficial in large-scale and complicated pipeline systems (Duan et al. 2017a, 2017b, 2020; Negharchi & Shafaghat 2021). However, the deviation of wave damping during energy dissipation makes it a necessity for increasing the accuracy of the IAB model. Since the compressibility is the fundamental cause of pressure waves in liquid, the energy dissipation caused by fluid volume compression or expansion needs to be well considered. The main purpose of this paper is to modify the IAB model with energy dissipation generated by compression–expansion (CE) of the fluid. First, the development of the CE effect term was derived from compressible Navier–Stokes equations. Second, the modified one-coefficient IAB-CE model was calculated by the method of characteristic (MOC), and then the simulation results were compared with the published experimental data. Finally, a quantitative error analysis was carried out to evaluate the prediction accuracy of the developed model.

According to the following assumption (Chaudhry 2014), the one-dimensional model of transient pipe flow is established:

• (1)

  The flow is one-dimensional with cross-sectional average quantities.

• (2)

  The flow velocity is small compared with the pressure wave speed.

• (3)

  The pressure wave speed is used to represent the compressibility of the fluid.

• (4)

  The pipe wall is elastic to allow for pipe deformation.

The one-dimensional governing equations for transient pipe flow can be expressed as follows:

Integrating Equation (11) over a control volume of length dx with a constant cross-section A, the control volume can be presented in Figure 1.

It shows that the term (⁠⁠) acts as a governing factor for the CE effect on the unsteady friction. It is also expected that the CE effect is more significant for shorter pipes under a large Reynolds number. The approximate order magnitude of each individual term can be analyzed using their scales. For most pipe flow cases when the Reynolds number is lower than 10⁵, the coefficient k is at a constant order of 10⁻². In product with other non-dimensional terms, it suggests that the effect of the CE becomes insignificant when k_d < 10 at low Reynolds number, while it obtains a 10 times order magnitude compared to other unsteady frictional terms in turbulent flow when k_d reaches an order of 10²_.

Because a general solution of the partial differential Equations (1) and (2) is not available, numerical procedures must be used. In this paper, the numerical integration is formed by using the MOC. Dividing the pipe into N segments equally, and are the spatial step and the time step for wave propagation, respectively. In this way, with the pipeline axis as the horizontal axis and time as the vertical axis, the grid for solving the characteristic equation is constituted. In order to facilitate the subsequent discretization of the unsteady friction term, the rectangular grid will be used for the MOC calculation, as is shown in Figure 2.

The modified one-coefficient IAB model in Equation (23) includes a diffusion term, which can be discretized according to the discrete grid, as shown in Figure 3.

As the purpose of this paper is to investigate the modified one-coefficient IAB model, a simple reservoir–pipe–valve system, which gives an extremely strong transient event after the valve is suddenly closed, is adopted.

The schematic diagram of the experimental device is shown in Figure 4, and the initial velocity of water in the pipe could be changed by adjusting the water level of the tank.

During the experiment, the pressure fluctuations at the middle of the pipe and the valve, namely the monitoring points 1 and 2, are recorded. In this paper, four different sets of experimental data from published literature (Bergant et al. 2001; Adamkowski & Lewandowski 2006) are used to analyze the modified one-coefficient IAB model. System parameters for each test case are listed in Table 1.

For Bergant's experiment, after grid independence verification, the number of nodes is 29, so the space and time steps are set as 1.3286 m and 0.001 s, respectively. When the initial flow velocity is 0.1 m/s, the comparison between the simulation results of the one-coefficient IAB model and the modified one-coefficient IAB-CE model with the experimental data is as shown in Figure 5.

For the modified one-coefficient IAB model, the value of k_d needs to be determined. The determination of the related second viscosity coefficient is still difficult involving complicated empirical methods. Therefore, this paper adopts the commonly used trial-and-error method in the IAB unsteady friction model to determine the value of the coefficient, which is taken as 10 after repeated comparison calculation. From Figure 5, it is found that the results of the modified one-coefficient IAB model are slightly more consistent with the experimental data, and the waveform change curve near the equilibrium pressure also agrees with the waveform measured in the experiment to an acceptable degree.

When the initial flow velocity is 0.2 m/s, the pipe flow is turbulent. The trial-and-error method is also used to determine the coefficient k_d, which is 80. The results are shown in Figure 6, and it is found that the amplitude and phase of the pressure wave calculated by the modified one-coefficient IAB model are in better agreement with the experimental data than the original model.

It is also noted that the improvement in wave prediction is more obvious than that in the laminar flow. This can be explained by the dimensional analysis mentioned above, that the CE effect becomes small at a low Reynolds number. For the laminar transient flow, greater shear stress caused by a large velocity gradient enhances the friction loss which should be emphasized more than that of the oscillating loss (Brunone et al. 2004). In such cases, the CE effect may not be a dominating factor compared with the modification of the decay coefficient in the unsteady friction term. Since the estimation of decay coefficient draws no conclusive formula, a variable coefficient might be an option for the laminar transient model (Pezzinga 1999).

When the initial flow velocity is 0.3 m/s and the coefficient is 110, the corresponding results are presented in Figure 7. It is seen that except for the deviation between the minimum pressure value of the first cycle and the experimental data, the simulation results of the modified IAB model are better than the original model in terms of the phase waveform distortion of the pressure wave.

For Adamkowski's experiment, after grid independence verification, the number of nodes is 33, so the space step is 3.07 m, and the time step is 0.0024 s. According to the trial-and-error method, the coefficient k_d is 190, and the results are given in Figure 8. It is found that the waveform distortion near the equilibrium pressure is more consistent with the measured waveform, which is similar to the results under other conditions.

It is seen that the calculation results from the modified model show more obvious improvement in turbulent unsteady flow. It is also found that the coefficient k_d in the modified dissipation term increases with the increasing Reynolds number.

The k_d values fit well with Re using this empirical correlation (R² = 0.998), suggesting that the Bradley equation can be applied to provide a referenced baseline for the coefficient k_d.

The results of the calculated relative error are shown in Figures 9–12. In each figure, (a) and (b) represent the result from the original IAB model and the modified IAB model, respectively. It is found that the modified IAB model with CE effects gives significantly better results, and its relative errors are smaller than those of the original IAB model generally. For the four cases with different initial flow velocities, the average relative errors of the modified IAB model are 7.63, 16.38, 35.02, and 10.69%, respectively. These results are, respectively, improved by 0.26, 2.03, 9.56, and 36.67%, compared with the calculation results using the original IAB model. Therefore, a more satisfying computational accuracy can be obtained by considering the energy dissipation caused by CE effects in the momentum equation.

For the unsteady pipe flow, the estimation of peak pressure and valley pressure is important in process design, as they may largely affect the safety of devices, such as pumps, valves, and connected vessels. The average relative errors of each pressure head peak and valley are listed in Table 2. The maximal test head is used as the denominator herein to avoid distortion when compared with data that deviate largely from the initial value. The comparative results show that the modified IAB-CE model can provide satisfied estimations for both head peaks and head valleys, especially in test cases with a large flow velocity. This is similar to the previous results found in the whole head curve analysis. The improvement in prediction for peak and valley proves that the modified IAB-CE model is beneficial for investigating transient wave characteristics.

It should be noted that the continuous effort in IAB model development is still expected for the laminar pipe flow. A more comprehensive evaluation of the CE effect can be made by considering possible synchronous effects on other constant and unsteady frictional terms, as well as taking into account the separated accelerating or decelerating phase in future work. However, this may require the multiple coefficients method with enhanced model complexity.

In this paper, the energy dissipation caused by CE effects was derived from the compressible Navier–Stokes equations. It is found that the product of the second-order partial derivative of velocity in space and the second viscosity coefficient could express the energy dissipation term. According to the presented results, the mathematical model of transient pipe flow as the original one-coefficient IAB model with the CE energy dissipation term could obtain a more satisfying agreement between the results of calculation and the experimental data. The relative errors are improved by 0.26, 2.03, 9.56, and 36.67%, respectively, and the prediction for wave peak and valley is also improved. Furthermore, it is also found that the coefficient k_d of the modified dissipation term increases with the increasing Reynolds number, and their relationship can be well established by the Bradley equation, which is beneficial to provide a referenced baseline for obtaining the dissipation coefficient.

This work has been supported by the National Natural Science Foundation of China (Grant No. 21978227), the Fundamental Research Funds for the Central Universities (Grant No. xjh012019022), the Natural Science Basic Research Program of Shaanxi Province (Grant No. 2021JQ052), and the Hong Kong Scholars Program (Grant No. XJ2020045).

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