## Abstract

The quasi-steady friction model is generally adopted in water hammer simulation in pipe network systems, which cannot accurately reflect the attenuation of pressure, while the existing unsteady friction model is challenging to use in complex pipe network systems. In this study, a convenient method for treating the friction term is proposed based on the Moody diagram. The attenuation process of water hammer pressure can be accurately reflected by reading the relationship curve between Reynolds number and the Darcy friction factor in the pipeline transient process. Combined with the classical water hammer experiment and the long pipe valve closing experiment in our laboratory, the accuracy of this model is verified, and the influence of absolute roughness (*e*) and Reynolds number (*Re*) on the model was analyzed as well. The results show that the pressure attenuation using the Method of Characteristics (MOC) and the proposed friction model has a good agreement with the experimental data. The absolute roughness has little influence on the results in hydraulically smooth pipe, while the minimum Reynolds number has a significant influence. When selecting the minimum Reynolds number, 2% ∼ 5% of the initial flow rate is recommended for calculation.

## HIGHLIGHT

In this study, a convenient method for treating the friction term is proposed based on the Moody diagram. The attenuation process of water hammer pressure can be accurately reflected. The sensitivity of the model was analyzed. It can can not only maintain the convenience of the standard MOC but also reflect the attenuation of pressure in the pipe.

## INTRODUCTION

At present, unsteady friction models are used in the calculation of water hammer in pipelines, and it is believed that the Darcy friction factor is changing during the transition process. The attenuation simulation error of water hammer pressure is small. The accuracy of water hammer pressure attenuation is very important in some applications. It is needed to correctly evaluate the overlapping of pressure waves. It is also used for fault detection, such as pipeline leakage detection models based on inverse analysis of water hammer (Wang *et al.* 2002; Capponi *et al.* 2020). In the process of hydraulic transients, with the change of water flow velocity inside the pipe, the shear stress of the pipe is complicated, and the actual frictional resistance is also changing constantly. Brunone measured the instantaneous velocity profiles by an ultrasonic Doppler velocimeter in a long PEHD pipe. A simple stress model is combined to local velocity measurements to evaluate the instantaneous wall shear stress (Brunone & Berni 2010).

A lot of work has been done on the study of unsteady friction in pipelines. Based on the results of practical pipe research and empirical formula, Moody proposed a Moody chart similar to the experimental results of an artificial rough pipe (Moody 1944). Daily proposed an empirical correction model for unsteady frictional resistance calculation (Daily *et al.* 1956). It was believed that the instantaneous wall shear stress was directly proportional to the instantaneous acceleration. According to the experimental results, it was found that the unsteady shear stress was positive when the fluid was accelerating and negative when the flow was decelerating. For a short time, there is no wave superposition. For accelerating flows, the instantaneous wall shear stress is larger than the steady, whereas for decelerating flows it is smaller or almost equal to the steady. Regardless of acceleration or deceleration, the instantaneous wall shear stress is larger than the steady wall shear stress, when the overlapping of more and more reflected waves (Brunone & Golia 2008; Brunone & Berni 2010). Zielke first proposed a weighted function model based on convolution and gave the analytic expression of the friction term of pipeline laminar flow (Zielke 1968). However, the model ignored the convection term and saved the value of each step in the calculation, so the memory consumption was large, and the calculation was slow. Considering the influence of boundary layer and recent energy dissipation, Wood and Funk studied the influence of frequency-dependent effects of transients and friction on system reactions (Wood & Funk 1970). Suo and Wylie proposed a frequency-dependent friction factor model using Fourier transform (Suo & Wylie 1989). Brunone, based on Daily's research and combined with experiments, introduced instantaneous local acceleration and convective acceleration and constructed the classic IAB model (Brunone *et al.* 1995), but it was calculated by an empirical formula, and interpolation calculation was required by the Method of Characteristics (MOC). Vardy and Hwang derived the quasi-two-dimensional model by using the concentric cylindrical ring model (Vardy & Hwang 1991). Vardy proposed a new weighting model for turbulence based on Zielke's model (Vardy *et al.* 1993). Silva-Araya and Chaudhry defined the energy dissipation factor to modify the friction term in the governing equation (Silva-Araya & Chaudhry 1997). Pezzinga made a similar diagram and applied dimensional analysis to determine the dimensionless parameters of unsteady friction (Pezzinga 2000). Bergant and Vitkovsky modified the convection term of the IAB model and proposed the MIAB model (Bergant *et al.* 2001; Vitkovsky *et al.* 2006a, 2006b). Brunone studied pressure attenuation and energy dissipation in laminar flow (Brunone *et al.* 2004). Duan and Ghidaoui showed that the viscoelastic effect plays an important role in the late damping and phase shift (Duan *et al.* 2010). Duan studied the influence of pipe length and diameter on unsteady friction (Duan *et al.* 2012). Duan verified the weighting function-based model and the instantaneous acceleration-based model (Duan *et al.* 2017). The differences between models and data are measured using the local transient analysis norm and the integral total energy norm along the pipeline. Sahad Khilqa established an empirical formula for damping ratio by means of quantity-value analysis and nonlinear regression formula, and analyzed damping uncertainty (Khilqa *et al.* 2019a, 2019b).

The unsteady friction model can accurately reflect the attenuation of water hammer pressure. However, due to the complexity of the unsteady friction model, interpolation is needed in the simulation using the method of characteristics, and it is also difficult to deal with complex boundaries. The existing results are only applied in simple pipelines in the laboratory, which are difficult to use in actual water supply projects. Motivated by the assumption of turbine characteristics that the turbine characteristic curve obtained from the steady-state model experiment was also suitable for the transient state (Chaudhry 2014). The importance of unsteady friction in rapidly decelerating flows diminishes with the initial Reynolds number. The accuracy of existing convolutional unsteady friction model are a function of the initial and not a time-dependent Reynolds number (Meniconi *et al.* 2014). In this study, the relationship between friction coefficient and Reynolds number in the Moody diagram is used to reflect the change of friction following the treatment of pump and turbine characteristics in the transition process. Coupled with the standard MOC, the water hammer pressure of pipeline can be calculated, which can not only maintain the convenience of the standard MOC but also reflect the attenuation of pressure in the pipe. Under the premise of ensuring the accuracy of pipeline pressure attennuation, the presented method has the potential to be used in practical engineering calculation, while the physical mechanism is not discussed in this study.

The rest of this paper's format is as follows: in the second part, the mathematical model of quasi-steady friction and the Moody diagram based approach friction are introduced. In the third part, Bergant's classical experiment was respectively adopted to verify the model proposed in this study and the influences of different absolute roughness(*e*) and different minimum Reynolds number (*Re _{m}*) were discussed respectively. The fourth part is the conclusion.

## MATHEMATICAL MODEL

### Mathematical model of unsteady flow in pressurized pipelines

*V*is the flow velocity, m/s;

*t*is time, s;

*H*is the piezometric head, m;

*x*is the distance along the centerline of the pipe, m;

*g*is the acceleration of gravity, m/s

^{2};

*f*is the Darcy friction factor, in quasi-steady friction,

*f*= ;

*D*is the pipe diameter, m;

*a*is the wave velocity of water hammer, m/s. At present, MOC is a commonly used numerical calculation method for hydraulic transient solution of the pipeline system. Compared with other methods, MOC has better accuracy, easy programming, and can handle more complex systems.

*λ*can obtain two equations with

*H*and

*V*as dependent variables, which can be transformed into two sets of ordinary differential equations. The difference equation is obtained by integrating the equation:where,

*Q*is introduced to replace

*V*in the formula.

### Friction model based on the Moody diagram

*e*is the absolute roughness of the pipeline and

*D*is the diameter of the pipeline. As can be seen from Figure 1, when the relative roughness of the pipeline remains unchanged, the loss coefficient (

*f*) along the pipeline increases with the decrease of Reynolds number. However, in a conventional pipeline water hammer calculation,

*f*is a constant, inconsistent with reality. By using the MOC, the Reynolds number of the intersection point of any internal grid is:

The Darcy friction factor (*f)* is read through the Moody diagram and according to the flow rate and Reynolds number in the transient process. Then the friction term in the transition process is changed.

The Moody diagram is the experimental result obtained under the condition of steady flow. When the *Re* is large, the inertial force plays a leading role in the pipeline. When the relative roughness is constant, with the decrease of the velocity *V* in the pipeline, the smaller the Reynolds number is, the larger the Darcy friction factor (*f*) is. However, when the *Re* approaches zero infinitely, *f* tends to infinity.

*f*) was obtained by reading the Moody diagram, and the error was large. When the flow velocity decreases to a certain value, the calculated value of Equation (5) cannot truly reflect the relationship between the real viscous force and the inertia force of the pipeline. In the momentum equation the friction force per unit weight has to be computed. Although

*f*tends to infinity for Re approaching zero, the friction force for laminar flow is:

*f*goes to infinity. To solve this problem, a minimum value for the

*Re*or the speed value should be specified such that, if the computed value is less than the specified minimum, then

*V*is assumed equal to

*V*

_{min.}*Re*

_{m}is the minimum allowable Reynolds number. Use Equation (5) to find

*Re*

_{m}and read

*f*

_{c}, then continue calculating the next pipe section. When the instantaneous Reynolds number is less than

*Re*

_{m},

*f*is a constant

*, f*=

*f*

_{c}. Silva-Araya and Chaudhry used 2% ∼ 5% of the initial steady-state velocity as the basis to determine the energy dissipation factor (Silva-Araya & Chaudhry 1997). For low-speed laminar flow, the error caused by the friction term is small. In this paper, 2% ∼ 5% of the initial steady-state velocity is also used to obtain

*Re*according to Equation (5),

_{m}*f*

_{c}can be obtained according to Equation (7):

## RESULTS AND DISCUSSION

### The experiment of Bergant

In this section, the friction treatment method proposed in this paper is studied based on the classic valve closing water hammer experiment of Bergant. In this experiment, the water level of the upstream water tank was 32 m, and the valve closing time was 0.009 s. The attenuation of the water hammer in the pipeline was observed after the closure of the downstream valve. Specific experimental parameters are shown in Table 1:

Length (m) . | Diameter (m) . | Friction . | Wave velocity (m/s) . | Velocity (m/s) . | Kinematic viscosity (m^{2}/s)
. | Gravitational acceleration (m/s^{2})
. | Closing valve time (s) . |
---|---|---|---|---|---|---|---|

37.2 | 0.022 | 0.037 | 1,319 | 0.3 | 1.184 × 10^{−6} | 9.8 | 0.009 |

Length (m) . | Diameter (m) . | Friction . | Wave velocity (m/s) . | Velocity (m/s) . | Kinematic viscosity (m^{2}/s)
. | Gravitational acceleration (m/s^{2})
. | Closing valve time (s) . |
---|---|---|---|---|---|---|---|

37.2 | 0.022 | 0.037 | 1,319 | 0.3 | 1.184 × 10^{−6} | 9.8 | 0.009 |

According to Equation (8), the initial Reynolds number can be obtained as 5,574, the pipeline loss along the pipeline can be obtained through the measured valve pressure, and the Darcy friction factor can be calculated as 0.037. By reading the Moody diagram, the corresponding curve *e/D* is found to be 0.00045, where *e* = 0.01 mm.The initial velocity of the experiment was 0.3 m/s, ranging from 0.006 m/s to 0.015 m/s. According to Equation (8), the assumed minimum Reynolds number range of the experiment was 111 ∼ 279. Taking *Re _{m}* = 200 (

*V*= 0.01 m/s) as the specified minimum value, put it into Equation (7) to obtain the critical Darcy friction factor (

*f*

_{c}) as 0.32. When

*Re*>200, read the Moody diagram to obtain

*f*for calculation; when

*Re <*200, calculate as

*f*= 0.32. The curve of the Darcy friction factor used in this calculation is shown in Figure 2. When

*Re*reaches 2,300, the flow pattern changes, and

*f*= 0.4951.

Figure 3(a) and 3(b) show the pressure of the valve and the middle part of the pipeline calculated by using quasi-steady and Moody diagram-based approach friction models, respectively. The horizontal axis represents the time after the valve is closed, and the vertical axis represents the pressure at the measured point.

Because the upstream water tank area is much larger than the cross-sectional area of the pipeline, the water level at the upstream water tank changes less, which can be regarded as a constant value. The maximum and minimum values of pressure in the quasi-steady friction model, the Moody diagram-based approach friction model, and the measured data in the process of pressure change all appear in the first phase.The pressure changes of the two models fit well with the measured data in the first half of the phase length, and the maximum value of both models is 72.33 m, 0.57 m less than the measured data. Therefore, the quasi-steady friction model can be used to approximate the maximum value of water hammer in practical engineering. The minimum value calculated by the quasi-steady friction model is −8.09 m, while the measured data is −6.15 m. The minimum value calculated by the Moody diagram-based approach friction model is −6.31 m. There is a certain safety margin to calculate the minimum value of water hammer using the quasi-steady friction model in practical engineering. Over time, the damping effect of friction on the water hammer wave becomes more and more evident. The quasi-steady friction model ignores the change of friction in the transition process, and the loss coefficient along the passage in the calculation is smaller than the actual situation. The difference between the model and the experimental result is gradually magnified over time. The pressure decays slowly under the quasi-steady friction model. The Moody diagram-based approach friction simulation method presented in this paper is in good agreement with the experimental data. Although there is a certain phase difference, it can reflect the attenuation of water hammer pressure better.

In order to eliminate the judgment of an individual experiment on the results of this method, the water hammer valve closing experiment was carried out again on this basis. A wider range of valve closing times is included, and the pipe length and pipe radius are increased to further verify the feasibility of this method.

### Influence of absolute roughness on the Moody diagram-based approach friction

*Re*,

*f*and

*e/D*. When

*e*is changed, the curve read on the Moody diagram will change accordingly. In this section,

*e*= 0.001 mm,

*e*= 0.01 mm,

*e*= 0.1 mm, and

*e*= 5 mm are selected respectively. Bergant's experiment is still used for other boundary conditions. The pressure attenuation of water hammer with different pipeline roughnesses was compared. Due to the different roughnesses, the initial pressure of the valve is different, at 31.7252 m, 31.7207 m, 31.6804 m, 31.4257 m, and 30.6476 m, respectively. When , even if the order of magnitude changes, the friction along the path is not much different when it is stable. Figure 5 shows the attenuation range of the maximum pressure in each period compared with the maximum pressure in the previous period under the Moody diagram-based approach friction model. Under the same calculation model, the attenuation rate is calculated according to the maximum pressure values

*H*

_{T1max}and

*H*

_{T2max}in the period occurring in two adjacent periods

*T1*and

*T2*.

As shown in Figure 4, the maximum water hammer pressure at the end of the first phase under five types of *e* is 72.3 ± 0.1 m, indicating that *e* has little influence on the maximum water hammer pressure. With the increase of *e*, the water hammer pressure decreases faster. Figure 5 shows that the greater *e* is in the same period, the greater *R _{a}* is. In the same case of

*e*,

*R*becomes smaller and smaller with the increase of the attenuation period. Under different

_{a}*e*, Δ

*R*becomes smaller and smaller with the increase of the attenuation period. When

_{a}*e*changes between 0.001 mm and 0.1 mm, the values of relative roughness give rise to friction Reynolds numbers

*Re** less than 5, and Δ

*R*< 0.05%, indicating that absolute roughness has little effect on the attenuation of water hammer in hydraulically smooth pipe. It should be noted that the relative roughness used is the ratio of the average height of the surface roughness to the diameter of the pipe and is not equivalent to the instantaneous friction of the pipe.

_{a}### Influence of Re on the Moody diagram-based approach friction

When the velocity is relatively large; that is, when the Reynolds number calculated according to is large, the *f* along the passage is calculated by the Moody diagram. When *Re* is reduced to a certain value, the critical friction(*f*_{c}) is selected for calculation. In order to study the influence of minimum Reynolds number on the calculation results, this section takes *Re _{m}* = 10,

*Re*= 200, and

_{m}*Re*= 2,000 as nodes for calculation (the corresponding critical friction along the path is 6.4, 0.32, and 0.032, respectively). The attenuation of water hammer pressure is analyzed and compared; that is, when the Reynolds number calculated in the pipeline is greater than the minimum Reynolds value, the

_{m}*f*along the pipeline is read according to the Moody diagram; when the Reynolds number is less than the minimum value, the

*f*

_{c}corresponding to the minimum Reynolds number is used. The calculated results and experimental data are shown in Figure 6.

Select *Re _{m}* = 2,000. If the Reynolds number calculated according to the velocity is greater than 2,000, read the Darcy friction factor,

*f,*along the pipeline according to the Moody diagram; if it is less than 2,000, use

*f*= 0.032, the corresponding

*Re*is 2,000. In this case, the attenuation of the water hammer is relatively slow; that is, when the Reynolds number node is large, the effect of viscous force in the attenuation process of the water hammer is underestimated. According to the Moody diagram, when

*Re*= 10, the loss coefficient along the pipe,

*f*= 6.4. When

*Re*= 10 is selected, it is found that the attenuation rate of the water hammer wave is much faster than the measured value, indicating that the Darcy friction factor,

_{m}*f,*used in the calculation process is larger than the actual situation. In the process of hydraulic transition, the effect of the viscous force is amplified, and the effect of the inertia force is reduced. According to 2% ∼ 5% of the initial steady-state velocity as

*V*, the Re

_{min}*can be obtained. In this experiment,*

_{m}*V*= 3.3%

_{min}*V*was used, corresponding to

*Re*= 200. As shown in Figure 6, when the minimum Reynolds number is 200, the numerical calculation results are consistent with the measured results.

_{m}### Pipeline closing valve experiment

_{2}O ∼ 50 mH

_{2}O, the accuracy class is 0.5%F.S, the frequency response is 20 kHz, and the temperature range is 0 °C ∼ 50 °C. A flow regulating valve is set at 1 m behind the butterfly valve to adjust the flow of the pipeline. The distance between the flow regulating valve and the downstream water tank is 1 m. In the experiment, a triangular thin-walled weir was used to measure the flow. The pipe flow of the right triangle thin wall weir can be calculated by Equations (10) and (11):where

*Q*is the pipeline flow rate, m

^{3}/s;

*C*

_{0}is the flow coefficient of the triangular thin-walled weir;

*H*is the head of the weir, m;

*P*

_{1}is the distance from the barrier to the bottom of the weir, m;

*B*is weir width, m. According to Thompson's experiment, when the barrier angle is 90°, the approximate value is 1.4.

DTFXP-B16A-NN ultrasonic flowmeter and weighing method were used to check the flow. The working principle of the ultrasonic flowmeter is that two transducers are used as ultrasonic transmitter and receiver at the same time, transmitting and receiving frequency modulated pulse sound energy alternately between the two transducers. According to the time difference, the velocity of liquid in the pipe can be obtained. The weighing method calculates the flow rate by weighing the water quality and the required time in an iron bucket. The error range of flow values at each flow rate was required to be within ±5%. The pipe and instrument are left to rest for 30 minutes after adjusting a set of flows. If the flow does not change, the flow is considered to be stable.

*T'*of two adjacent wave peaks. The measured wave velocity is obtained as 1,317 m/s through Equation (12). The water level of the upstream water tank was 6.06 m, and the initial flow rate was 0.12 ± 0.005 m/s. The water hammer experiment was repeated ten times to avoid an experimental error. Compared with the numerical simulation, the experimental parameters are shown in Table 2.

Length (m) . | Diameter (m) . | Friction . | Wave velocity (m/s) . | Velocity (m/s) . | Kinematic viscosity (m^{2}/s)
. | Gravitational acceleration (m/s^{2})
. | Closing valve time (s) . |
---|---|---|---|---|---|---|---|

150 | 0.05 | 0.0455 | 1,317 | 0.118 | 1 × 10^{−6} | 9.8 | 0.039 |

Length (m) . | Diameter (m) . | Friction . | Wave velocity (m/s) . | Velocity (m/s) . | Kinematic viscosity (m^{2}/s)
. | Gravitational acceleration (m/s^{2})
. | Closing valve time (s) . |
---|---|---|---|---|---|---|---|

150 | 0.05 | 0.0455 | 1,317 | 0.118 | 1 × 10^{−6} | 9.8 | 0.039 |

The closing time of the valve is 0.039 s. According to the common equivalent roughness table of the pipeline, *e* = 0.05 is selected. According to the corresponding *f*-*Re* curve read in the Moody diagram, the initial Darcy friction factor, *f,* along the pipeline is 0.0455. The initial steady flow of the experiment was 0.118 m/s, and *V _{min}* was between 0.0024 m/s and 0.0059 m/s. According to Equation (5), the minimum Reynolds number was assumed to be in the range of 120 ∼ 295.

*Re*= 200 (

*V*= 0.004 m/s) was taken as the specified minimum value, and the corresponding critical friction,

*f*

_{c,}along the pipeline was 0.32. Comparison between the measured pressure in front of the valve and the numerical simulation results is shown in Figure 8.

The numerical simulation results of the unsteady friction model based on the Moody diagram can be well matched with the measured data. Both the experimental data and the Moody diagram based approach friction model show that the water hammer wave tends to decay with time. In the first cycle, the maximum value of water hammer pressure is 21.1578 m, and the value calculated by numerical simulation is 21.8441 m, and the relative error between the two is 3.14%. The minimum value of water hammer pressure is −7.2893 m, and the value calculated by numerical simulation is −7.1352 m. The relative error between the two is 2.11%. The measured wave peak values in the following periods have certain oscillations. The other measured data are very close to the numerical results of the Moody diagram based approach friction model. It shows that the method of the Moody diagram based approach friction is effective.

## CONCLUSION

Although CFD can also accurately simulate pressure attenuation, it depends on the number of computational cells (Martins *et al.* 2017). For large-scale engineering, it takes up more memory and computing time. The method proposed in this paper has small amount of calculation and good accuracy. In this paper, the friction processing method for the transition process based on the Moody diagram is proposed. The curve of the relationship between Reynolds number and the Darcy friction factor in the transient process of the pipeline is read in real-time, and the water hammer pressure of pipeline is calculated by coupling with the standard characteristic line method. It can not only keep the convenience of the standard characteristic line method but also reflect the attenuation of pipe water hammer. The simulation results are compared with the classical water hammer experiment and the long pipe valve closing experiment in the laboratory. The effects of absolute roughness and Reynolds number on the model are analyzed. The results show that the quasi-steady friction model fitted well with the measured data in the first half of the phase length but could not reflect the attenuation of the subsequent water hammer. Using the method proposed in this paper, the simulation results are consistent with the experimental results, and the relative error is less than 5%, which can better reflect the attenuation of water hammer. The absolute roughness has little influence on the extreme value of water hammer pressure in hydraulically smooth pipe, while the minimum Reynolds number node has a more significant influence on the result. When selecting the nodes of the Reynolds number calculated by the model, the range of Reynolds number calculated by 2% ∼ 5% of the initial steady velocity can be used as the basis for judging the minimum Reynolds point.

In the future, the physical mechanism of this method will be further studied and analyzed, as well as the applicability for indirect water hammer, which is more common in practical engineering.

## ACKNOWLEDGEMENT

This work was supported by the National Natural Science Foundation of China (Grant Nos.: 52179062 and 51879087).

## DATA AVAILABILITY STATEMENT

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