## Abstract

This paper aims to modify the conventional one-coefficient instantaneous acceleration-based (IAB) model for better prediction of unsteady friction behavior. In this work, the energy dissipation caused by viscous stress during fluid volume compression–expansion (CE) was derived from the compressible Navier–Stokes equation. It is found that the energy dissipation term can be expressed by the product of the second-order partial derivative of velocity in space and the second viscosity coefficient. On this basis, a modified IAB-CE model was developed with the energy dissipation term and solved by the method of characteristic (MOC). The numerical results obtained from the modified model showed a good agreement with the four test cases, where the relative errors are improved by 0.26, 2.03, 9.56, and 36.67%, compared with the results from the original IAB model. The estimation for wave peak and valley is improved as well. Furthermore, the Bradley equation can be applied to establish the relationship between the dissipation coefficient and the Reynolds number. The modified model developed in this study takes into account the fluid CE effects and improves the prediction accuracy of wave amplitude of unsteady flow.

## HIGHLIGHTS

The numerical investigation of unsteady pipe flow with an efficient one-coefficient model.

The fluid compression–expansion effect was considered in the energy dissipation term.

The dissipation coefficient and the Reynolds number conform to the Bradley equation.

The wave amplitude prediction of friction flow was improved by the modified IAB-CE model.

## NOTATIONS

*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

## INTRODUCTION

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.

## METHODOLOGY

### One-dimensional unsteady friction 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:

*H*and flow velocity

*V*. The

*J*in Equation (2) stands for unsteady friction which could be calculated by the one-coefficient IAB model (Pezzinga 2000):where represents the signal of the instantaneous mean velocity, and

_{u}*k*is the empirical coefficient that could be determined by the trial-and-error method or the empirical formula (Reddy

*et al.*2012):where , .

### Derivation of the CE term

*F*is the mass force, and

*s*is the tensor of strain rate. When the fluid is incompressible,

*s*is zero in Equation (7); therefore, the momentum equation does not include the second viscosity coefficient representing the CE effects of the fluid.

_{kk}*x*is the coordinate along the pipe axis and ignoring momentum equations and velocity components in other directions as well as the mass force, the momentum equation in the

*x*-axis is 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.

*R*is the pipe radius. Integrating each term of Equation (11) in the control volume,

*V*and ignoring the small quantity, Equation (11) can be expressed as follows:

*J*and the constant coefficient as

_{d}*k*, and to eliminate the influence of the coordinate system, the CE effect term could be written as follows:

_{d}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^{5}, the coefficient *k* is at a constant order of 10^{−2}. 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*reaches an order of 10

_{d}^{2}

_{.}

### Method of solution

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.

### Experimental setup

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.

Parameters . | Bergant's experiment . | Adamkowski's experiment . | ||
---|---|---|---|---|

Upstream tank head H_{res} (m) | 32 | 128 | ||

Initial flow velocity V_{0} (m/s) | 0.1 | 0.2 | 0.3 | 0.94 |

Re | 1,870 | 3,750 | 5,600 | 15,843 |

Wave speed a (m/s) | 1,319 | 1,298.4 | ||

Valve closing time t_{c} (s) | 0.009 | 0.003 | ||

Density (kg/m^{3}) | 998.9794 | 997.5906 | ||

Dynamic viscosity () | 1.1227 × 10^{−3} | 0.947 × 10^{−3} |

Parameters . | Bergant's experiment . | Adamkowski's experiment . | ||
---|---|---|---|---|

Upstream tank head H_{res} (m) | 32 | 128 | ||

Initial flow velocity V_{0} (m/s) | 0.1 | 0.2 | 0.3 | 0.94 |

Re | 1,870 | 3,750 | 5,600 | 15,843 |

Wave speed a (m/s) | 1,319 | 1,298.4 | ||

Valve closing time t_{c} (s) | 0.009 | 0.003 | ||

Density (kg/m^{3}) | 998.9794 | 997.5906 | ||

Dynamic viscosity () | 1.1227 × 10^{−3} | 0.947 × 10^{−3} |

## RESULTS AND DISCUSSION

### Model calculation

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.

*et al.*2010) and was found applicable for multilayer sorption where most other equations fail (Walker

*et al.*1973). Some scholars also found that the diffusion constant and its affecting factors correlate to the Bradley equation (Barrer & Rideal 1939). In this study, further investigation establishes the relationship between the coefficient

*k*and the Reynolds number, which can be expressed by the Bradley equation as follows:

_{d}The *k _{d}* values fit well with Re using this empirical correlation (

*R*

^{2}= 0.998), suggesting that the Bradley equation can be applied to provide a referenced baseline for the coefficient

*k*.

_{d}### Relative error analysis

*H*

_{c}is the computed values obtained by the original IAB model or the modified IAB model,

*H*

_{m}is the experimental value, and

*H*

_{ini}is the upstream tank head which could be found in Table 1.

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.

Test cases . | Peak error . | Valley error . | ||
---|---|---|---|---|

IAB (%) . | IAB-CE (%) . | IAB (%) . | IAB-CE (%) . | |

Bergant (0.1 m/s) | 1.85 | 1.86 | 0.71 | 0.68 |

Bergant (0.2 m/s) | 1.77 | 1.16 | 1.52 | 0.36 |

Bergant (0.3 m/s) | 1.59 | 0.98 | 4.55 | 1.82 |

Adamkowski (0.94 m/s) | 5.70 | 3.60 | 3.46 | 2.03 |

Test cases . | Peak error . | Valley error . | ||
---|---|---|---|---|

IAB (%) . | IAB-CE (%) . | IAB (%) . | IAB-CE (%) . | |

Bergant (0.1 m/s) | 1.85 | 1.86 | 0.71 | 0.68 |

Bergant (0.2 m/s) | 1.77 | 1.16 | 1.52 | 0.36 |

Bergant (0.3 m/s) | 1.59 | 0.98 | 4.55 | 1.82 |

Adamkowski (0.94 m/s) | 5.70 | 3.60 | 3.46 | 2.03 |

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.

## CONCLUSIONS

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.

## ACKNOWLEDGEMENTS

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

## DATA AVAILABILITY STATEMENT

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