Unsteady friction model modi ﬁ ed with compression – expansion effects in transient pipe ﬂ ow

This paper aims to modify the conventional one-coef ﬁ cient instantaneous acceleration-based (IAB) model for better prediction of unsteady friction behavior. In this work, the energy dissipation caused by viscous stress during ﬂ uid 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 coef ﬁ cient. On this basis, a modi ﬁ ed IAB-CE model was developed with the energy dissipation term and solved by the method of characteristic (MOC). The numerical results obtained from the modi ﬁ ed 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 coef ﬁ cient and the Reynolds number. The modi ﬁ ed model developed in this study takes into account the ﬂ uid CE effects and improves the prediction accuracy of wave amplitude of unsteady ﬂ ow.

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

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. 2017aDuan et al. , 2017bDuan et al. , 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 proved 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 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 exist 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) furtherly 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 for 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 largescale and complicated pipeline systems (Duan et al. 2017a(Duan et al. , 2017b(Duan et al. , 2020Negharchi & 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.

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.
AQUA -Water Infrastructure, Ecosystems and Society Vol 00 No 0, 3

Uncorrected Proof
(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: Continuity equation: Momentum equation: where D is the pipe diameter, a is the wave speed, g is the gravitational acceleration, x is the coordinate along the pipe axis, t is the temporal coordinate, and f is the Darcy-Weisbach friction factor that could be computed by the following equations. Based on the critical Re for flow transition (Adamkowski & Lewandowski 2006;Urbanowicz & Zarzycki 2008), the Darcy-Weisbach friction factor or laminar flows (Re 2,320) depends only on flow characteristics (Re), according to the Hagen-Poiseuille law: For turbulent flows (Re . 2,320), it depends on flow characteristics (Re) and absolute pipe wall roughness (K/D) according to the Colebrook-White equation: Equations (1) and (2) belong to hyperbolic partial differential equations with unknown variates of the piezometric head H and flow velocity V. The J u in Equation (2) stand for unsteady friction which could be calculated by the one-coefficient IAB model (Pezzinga 2000): where sign(V) represents the signal of the instantaneous mean velocity, and k is the empirical coefficient that could be determined by the trial-and-error method or the empirical formula (Reddy et al. 2012): where C Ã ¼ 12:86=Re k , k ¼ log 10 (15:29=Re 0:0567 ).

Derivation of the CE term
The tensor expression of the Navier-Stokes equation is as follows: where m is the dynamic viscosity, m 0 is the the second viscosity, F is the mass force, and s is the tensor of strain rate. When the AQUA -Water Infrastructure, Ecosystems and Society Vol 00 No 0, 4 Uncorrected Proof fluid is incompressible, s kk is zero in Equation (7); therefore, the momentum equation does not include the second viscosity coefficient representing the CE effects of the fluid. As indicated by Bratland (1986), the transient friction can be analyzed with more detail by considering the radial variation of flow velocity characteristics. In order to reflect the energy dissipation caused by the CE effects in the one-dimensional unsteady friction model, the momentum equation is re-derived from Equation (7) under the compressible condition. Assuming 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: Equation (8) needs to be transformed into the form in the cylindrical coordinate system, and the conversion formula is as follows: where u is the circumferential coordinate, and r is the radial coordinate. Substituting Equations (9) and (10) into the Equation (8), and then, 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.
Taking into account the cross-sectional average velocity and pressure, where R is the pipe radius. Integrating each term of Equation (11) in the control volume, In Equation (16), the second viscous coefficient of flow with pressure wave is much larger than the dynamic viscous coefficient, so the dynamic viscous coefficient is ignored. The wall friction t w in Equation (17) is as follows: Substituting Equations (14)- (17) into Equation (11), and multiplying both sides by 1=rgA and replacing the cross-sectional average pressure P with the hydraulic head H. The result is as follows: Noticing that the second term on the right side of Equation (19) is the central difference form of the diffusion term, so for the whole pipe, the result is as follows: After expanding the full derivative of the velocity V and ignoring the small quantity, Equation (11) can be expressed as follows: Comparing Equation (21) and Equation (2), it could be found that, except for the constant friction term, Equation (21) includes the term related to the second viscosity coefficient, which is the dissipation caused by the CE effects. To define this term as J d and the constant coefficient as k d , and to eliminate the influence of the coordinate system, the CE effect term could be written as follows: Incorporating Equation (22) into Equation (5), the modified one-coefficient IAB model applied in this paper is obtained as follows: AQUA -Water Infrastructure, Ecosystems and Society Vol 00 No 0, 6

Uncorrected Proof
In order to identify the main parameters governing the hydraulic friction loss, dimensional analysis is carried out to estimate the order of magnitude of each term in the momentum equation. The variables in Equation (2) are normalized by the Joukowsky head (Wahba 2017), wave propagation time, and pipe length: And the dimensionless form of the momentum equation is expressed as follows: where the modified unsteady friction term is as follows: It shows that the term ( k d g V 0 L 2 ) 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 than other unsteady frictional terms in turbulent flow when k d reaches to an order of 10 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 segment equally, Dx and Dt 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.
According to the MOC, Equations (1) and (2) can be transformed into ordinary differential equations along with the characteristics lines with known quantity at points A and B. For the one-coefficient IAB unsteady friction model, the contained convective acceleration term @V=@x and the local acceleration term @V=@t could be discretized using the following equations, respectively: 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.
And the discretization at points P, B, and C can be expressed as follows: (Dx) 2 (33)

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;   Table 1.

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 hereinbefore 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 more emphasized than that of the oscillating loss (Brunone et al. 2004). In such cases, the CE effect may  Uncorrected Proof 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.
As a common mathematical form for expressing empirical correlation, the Bradley equation was originally employed to describe the adsorption behavior (Krupadam 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 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
In order to allow for a precise evaluation, a quantitative measure of the difference between the pressure head obtained by the original IAB model and by the modified IAB model has been made. The calculation method is as follows: where 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      Uncorrected Proof significantly better results, and its relative errors are smaller than that 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 device, 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 request in 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.