Experimental and Numerical Simulation of Water Hammer in Gravitational Pipe Flow with Continuous Air Entrainment

<jats:p>Water hammer is an undesired hydraulic shock phenomenon in water supply pipe systems. It is very important to simulate water hammer for preventing the hazard of over pressure. In order to predict the transient pressure caused by a valve closing in a gravitational pipe with continuous air entrainment, a numerical model based on the Lax-Wendroff format is established, and the matched boundary model is provided. Compared with the traditional methods, this study provides another access by considering the influence of the pipe flow velocity on the wave propagation to simulate transient processes. A corresponding experiment is conducted to optimize the numerical model. Based on the experimental result, an additional friction function is proposed to evaluate the influence of the air content on the attenuation. The result shows that the energy dissipation of the shock waves may be underestimated in air-water mixture flow using the common steady friction. By introducing the additional friction function, the improved model can more accurately simulate the attenuation of the water hammer in the gravitational pipe with continuous air entrainment. As there are plenty of practical water supply systems running with air content, the improved Lax-Wendroff Method (LWM) is valued in accurately predicting water hammer processes especially in those conditions.</jats:p>

section_content: As an undesired hydraulic phenomenon in pipe systems, water hammer is usually caused by a valve closing and opening, a pump shutting and starting, as well as misoperation in pipe systems [1, 2]. Various water hammer events commonly occur in water pipe systems, pump systems and hydro-power systems [3]. They can bring many risks, such as overpressures, vibrations, and bursts [4, 5]. It is very important to predict water hammer for protecting hydraulic system [6] [7] [8]. In the past it has been difficult to simulate water hammer pressure by mathematical algorithm. With the development of the computer, mathematical algorithms become more popular and improved greatly. 

Numerical simulation has now become the main approach for transient analysis [9] [10] [11]. The key of the water hammer calculation is to solve the hyperbolic partial differential equations. Although the water hammer equation is a closed-form expression, no theoretical analytical solution is established at present. In order to obtain the numerical solution, various methods are applied to water hammer simulation, typically including the Finite Difference Method (FDM) [12] and the Finite Element Method (FEM) [13]. According to the existing literatures, FDM is the most favored method. Especially the method of characteristic (MOC), as an FDM, is nowadays the most popular method. It is widely used to simulate the water hammers in various engineering practices, such as valve closing [14] [15] [16] [17] [18], and pump failure [19]. Not only that, it has also been applied to special system and boundary conditions. Sibetheros [20] improved the flexibility of numerical simulations by spline polynomials interpolations in MOC. Gong [21] proposed a reconstructed MOC to estimate the pipe parameters. Compared with conducted experiment on a single pipeline with a thin wall thickness, the reconstructed MOC is validated to be of good reliability. Travas [22] proposed a mixed MOC and FDM to consider unsteady friction. Equipped with a feature of relatively low complexity, the mixed method is easier to numerically implement. The Kelvin-Voigt model is introduced to MOC for considering pipe-wall viscoelasticity [23, 24]. Nevertheless, based on wave propagation principles, Wood [25] proposed a wave characteristic method (WCM) to simplify the analysis of a water network. According to the comparisons with experimental data and other traditional methods, WCM is confirmed to be of good accuracy except for a significant improvement on convenience. Guo [26] used a sliding mesh method to simulate a water hammer induced by ball-valve closure. They found that Fluent can better simulate the water hammer than the traditional method of characteristics, especially when considering the compressibility of water. In actual water supply systems, there are plenty of pipelines running with air content. In previous research, it has already been ascertained that the numerical solutions for two-phase flow is familiar in format to the model for one-phase liquid flow [27, 28]. Bergant et al. [29, 30] ascertained that water hammer attenuation may be significantly affected by unsteady friction, cavitation, and trapped air pockets. In their researches, various cases with different conditions are studied and analyzed. It is found that the water hammer attenuation has always been over estimated in traditional simplified methods. Zhou [31] did numerical simulations and conducted an experiment on air leakage effect in a water hammer process, which shows that air leakage has a significant cushioning effect which reduces the intensity of the water hammer. 

As a conventional numerical method, the Lax-Wendroff Method (LWM) was widely used in computational fluid dynamics (CFD) [32]. However, only a few studies have briefly reported LWM in a water hammer simulation [12, 33, 34], without considering the convective acceleration terms, and it has hardly been applied in water hammer practices. Considering the real working conditions of water supply systems, it is more important to improve the accuracy of predictions on a water hammer in a pipeline with air content inside. In this paper, water hammer was simulated by LWM for gravitational pipe flow with continuous air entrainment. To analyze the influence of air content on the wave propagation and attenuation, the general numerical model was established by introducing LWM to complete governing equations, and the matched boundary conditions were provided. Then an experiment was conducted to validate and improve the proposed model. Compared with the experiment, the maximum pressure and attenuation were simulated by the proposed model. However, the attenuation is heavily unevaluated in the numerical simulation. Based on the experimental result, an additional friction function was added to the model for approximating the water hammer attenuation. The application shows that the improved model can more accurately simulate the attenuation of water hammer. In water supply systems, especially under low pressure conditions, different air content may occur in the pipelines. It makes the wave speed and energy loss much different from the flows without air. The proposed improved LWM model can well simulate the transient processes of pipe flow with continuous air entrainment. It is valuable in accurately predicting extreme pressure in water hammers with air entrainment, and is helpful in guiding suitable operations in water supply systems for extreme pressure control. Further research of this study will be focused on the physical interactions between entrained air and water mass in air liquid mixture flow. Based on the equipped experimental platform, experiments in an inner visible pipeline will be conducted to trace the air entrainment, delivering, bubble forming and breaking processes by observing the actual inside transient process. 

Water 2018, 10, 928 3 of 16
section_name: Introduction

section_name: Numerical Modeling and Meshing

section_content: Many factors can affect water hammer, including the flow velocity, the pressure and the properties of the pipe wall and the fluid. In hydraulic transient, the water hammer is usually considered as an isothermal process, consequently, continuity and momentum equations can describe the dynamics characteristics of the water hammer. For a closed-pipe system, the complete dynamic governing equations of the water hammer can be written as follows [1]

where x and t are independent variables denoting location and time respectively. Correspondingly, h and v are dependent variables denoting hydraulic pressure and flow velocity respectively.
section_name: Governing Equations

section_content: The dynamic equations system consists of two hyperbolic partial differential equations. It is difficult to establish the general analytical solution [1]. MOC is usually used to establish the numerical solutions, without considering the convective acceleration terms. Considering the influence of convective acceleration terms on water hammer, LWM is used to solve the complete equations system without omitting convective terms. The Lax-Wendroff Method (LWM) is widely applied in CFD [32], but it is hardly used in water hammer simulations. Previously, Chaudhry [12] briefly reported the LWM solution of simplified water hammer equations. Equations ( 1 ) and (2) meet LWM in time and space, using a second-order Taylor expansion, then the discretized FDM can be written as follows

In these equations, time and space can be described by discrete nodes and mesh, as seen in Figure 1, where i is a series to denote the location of the nodes. When a pipe is divided as n segments, i ∈ {0, 1, 2,. . ., n}. i = 0 and i = n refer to the first and the last node respectively, and the corresponding locations are x = 0 and x = l. ∆t is the time step. Equations ( 3 ) and (4) are the discretized second-order Taylor expansion for velocity and hydraulic pressure. The left terms are the variables at time t + ∆t, and the right terms are the variables at time t. 

The unknown variables at t + ∆t can be solved after the first-order and second-order derivative terms being determined in Equations ( 3 ) and (4). According to Equations ( 1 ) and ( 2 ), the first-order derivative terms can be converted as

Then, the second-order derivative terms can be derived according to Equations ( 5 ) and (6). Equations ( 5 )- (8) establish the relationship between time derivatives and spatial derivatives. 

The first-order and second-order space derivatives can be approximated by the central differences as Equations ( 9 )- (12) [32]
section_name: Discretization Model

section_content: A time-space differences mesh is needed to simulate the transient process. Figure 1 shows a general space and time mesh. In the mesh, x denotes one-dimensional space, representing the distance from the inlet of the pipe, and ∆x is the space step. t is the time variable and ∆t is the time step. There are three kinds of nodes according to time and space: general interior nodes, initial nodes, and side nodes. In Figure 1, hollow circles denote the internal nodes, solid circles denote the initial nodes and rectangles denote the side nodes. At the initial state t = 0, the beginning of the first step, all the variables of initial nodes are known or can be determined according to the initial flow conditions. Then the internal nodes at the time t = ∆t, the ending of the first step, can be determined based on the parameters in the beginning of the current step. Then the side nodes will be determined according to the relative boundary conditions. Consequently, all the nodes were solved at the time t = ∆t, which is the ending of the current step and will be considered as the beginning of the next step. Analogously, the water hammer process can be determined by time advancement, if the initial boundary conditions are provided.
section_name: Meshing and Advancement

section_content: As shown in Figure 2, the gravitational pipe system consists of a long pipe, an upstream reservoir, and a downstream valve. Keeping the upstream water level constant, the flow is steady at the initial state. The water hammer will occur when the downstream valve close suddenly. In order to simulate the transient process by the proposed model, the matching boundary condition model is established as follows. It includes the upstream reservoir and the downstream control valve. 

Meshes and node types.
section_name: Boundary Conditions

section_content: As shown in Figure 2, the gravitational pipe system consists of a long pipe, an upstream reservoir, and a downstream valve. Keeping the upstream water level constant, the flow is steady at the initial state. The water hammer will occur when the downstream valve close suddenly. In order to simulate the transient process by the proposed model, the matching boundary condition model is established as follows. It includes the upstream reservoir and the downstream control valve. In this experiment, the detailed parameters of devices and equipment are listed as follows. The pipe connecting the upstream and downstream reservoirs is a steel pipe 29 m in length. Its internal diameter is 0. 107 m, and the thickness of the pipe wall is 0. 005 m. The elasticity modulus of the pipe material is 210 GPa. The control valve downstream is a butterfly valve, and its characteristics are illustrated in Section 3. 2. Along the pipeline, there are 26 test nipple joints distributed uniformly at one side of the pipe, available for connecting pressure sensors. The pressure sensors are of high precision, with the time interval of two captured pressure values being 0. 001 s. The water hammer experimental data is obtained by recording the pressure data of the sensors during the transient processes. In this experiment, the detailed parameters of devices and equipment are listed as follows. The pipe connecting the upstream and downstream reservoirs is a steel pipe 29 m in length. Its internal diameter is 0. 107 m, and the thickness of the pipe wall is 0. 005 m. The elasticity modulus of the pipe material is 210 GPa. The control valve downstream is a butterfly valve, and its characteristics are illustrated in Section 3. 2. Along the pipeline, there are 26 test nipple joints distributed uniformly at one side of the pipe, available for connecting pressure sensors. The pressure sensors are of high precision, with the time interval of two captured pressure values being 0. 001 s. The water hammer experimental data is obtained by recording the pressure data of the sensors during the transient processes.
section_name: Boundary Conditions

section_content: Figure 3 shows the logic boundary model of the pipe ends. For a large reservoir, the water level can be considered as a constant. Then the boundary model needs to include the constant water level. Thus, the pressure head can be determined by the boundary. 

With ∂h(0, t)/∂t = 0, combining Equation ( 6 ), the flow velocity of the upstream side nodes can be solved as

Water 2018, 10, 928
section_name: Upstream Boundary Model

section_content: With (0, ) 0 h t t ∂ ∂ =, combining Equation ( 6 ), the flow velocity of the upstream side nodes can be solved as

Boundary model of the pipe ends.
section_name: of 16

section_content: The valve is a flow control device, and its rapidly closing will cause water hammer. The boundary equation of the valve can be written as

where d C is the discharge coefficient, which depends on the valve performance curve and opening ratio at a specific time. The opening ratio is usually known according to the valve-closing process. Thus, the discharge coefficient may also be determined at any time. 

As for the head coefficient, convert Equation ( 6 ) to the different form of (
section_name: Downstream Boundary Model

section_name: Experimental Measurement

section_content: The hydraulic head loss is an important factor for water hammer. Traditional water hammer simulations commonly approximate friction resistance based on the steady status. In order to
section_name: Steady Hydraulic Friction

section_content: The valve is a flow control device, and its rapidly closing will cause water hammer. The boundary equation of the valve can be written as

where C d is the discharge coefficient, which depends on the valve performance curve and opening ratio τ at a specific time. The opening ratio is usually known according to the valve-closing process. Thus, the discharge coefficient may also be determined at any time. 

As for the head coefficient, convert Equation ( 6 ) to the different form of
section_name: Downstream Boundary Model

section_content: The hydraulic head loss is an important factor for water hammer. Traditional water hammer simulations commonly approximate friction resistance based on the steady status. In order to consider total head loss, including minor loss, various velocity experiments are conducted. Various steady flow conditions of different flow velocities are accessed by controlling the upstream reservoir water level. 

Two pressure sensors are set at the upstream side and downstream side of the pipeline respectively. The pressure difference of those two sensors represents the total head loss under the corresponding flow velocity. Figure 4 shows the total head loss change with the velocities, and the fitting friction curve is provided with an equivalent friction factor 0. 269. It shows that the steady friction can be obtained by the second power law. 

Water 2018, 10, x FOR PEER REVIEW 7 of 17 consider total head loss, including minor loss, various velocity experiments are conducted. Various steady flow conditions of different flow velocities are accessed by controlling the upstream reservoir water level. Two pressure sensors are set at the upstream side and downstream side of the pipeline respectively. The pressure difference of those two sensors represents the total head loss under the corresponding flow velocity. Figure 4 shows the total head loss change with the velocities, and the fitting friction curve is provided with an equivalent friction factor 0. 269. It shows that the steady friction can be obtained by the second power law.
section_name: Steady Hydraulic Friction

section_content: Discharge properties of the valve can greatly affect the water hammer process. A relatively smooth valve-closing process can effectively reduce the water hammer strength. To numerically simulate the experimental transient process, it is always necessary preparatory work to measure the valve's real discharge coefficient curve. To obtain the characteristics of the valve, the experimental data of various steady flow conditions are recorded. Adjusting the open ratio of the valve, according to the corresponding water head difference at two sides and the flow discharge through
section_name: Valve Performance

section_content: Discharge properties of the valve can greatly affect the water hammer process. A relatively smooth valve-closing process can effectively reduce the water hammer strength. To numerically simulate the experimental transient process, it is always necessary preparatory work to measure the valve's real discharge coefficient curve. To obtain the characteristics of the valve, the experimental data of various steady flow conditions are recorded. Adjusting the open ratio of the valve, according to the corresponding water head difference at two sides and the flow discharge through the valve, several mapping relationship dots are recorded based on Equation (15). After a least square method fitting, the performance curve of the valve in the experiment is shown in Figure 5.
section_name: Valve Performance

section_content: A related experiment was conducted in order to validate and improve the model. As shown in Figure 2, the experiment consists of a straight pipe, valves, and upstream reservoir. In the system, the length of the pipe is 29 m and the internal diameter is 0. 107 m. A water hammer is generated when the valve is closed in 2. 2 s. Figure 6 shows the water hammer pressure process. As seen in the figure, the wave-shaped curve shows the stable frequency and periods. 

Water 2018, 10, x FOR PEER REVIEW 8 of 17
section_name: Water Hammer Measurements

section_content: A related experiment was conducted in order to validate and improve the model. As shown in Figure 2, the experiment consists of a straight pipe, valves, and upstream reservoir. In the system, the length of the pipe is 29 m and the internal diameter is 0. 107 m. A water hammer is generated when the valve is closed in 2. 2 s. Figure 6 shows the water hammer pressure process. As seen in the figure, the wave-shaped curve shows the stable frequency and periods. According to the periods and length of the travel path, the experimental wave speed can be determined as follows The experimental wave speed is only about 151. 0 m/s. However, without taking into account the air, the intrinsic wave speed in the pipe should be about 1340 m/s according to the water hammer wave speed formula. Obviously, the experimental wave speed is far lower than for a pipe without air, because of air entrainment. For pipe flow with air content [1], the relationship between air content and the shock wave speed is K ρ According to the periods and length of the travel path, the experimental wave speed can be determined as follows

The experimental wave speed is only about 151. 0 m/s. However, without taking into account the air, the intrinsic wave speed in the pipe should be about 1340 m/s according to the water hammer wave speed formula. Obviously, the experimental wave speed is far lower than for a pipe without air, because of air entrainment. For pipe flow with air content [1], the relationship between air content and the shock wave speed is

where a is the actual wave velocity in a pipe system, K is the volume modulus, ρ is the density, E is the elasticity modulus, and m is the amount of substance of the gas in unit volume. Provided with the actual experimental wave speed, the air content can be evaluated as

For the case shown in Figure 6, the amount of air per unit volume is about 0. 205, and air volume content is about 0. 459% according to this measurement. 

Thus, all the necessary coefficients are established, which provide the basic parameters required to numerically simulate the water hammer process.
section_name: Water Hammer Measurements

section_name: Simulation and Analysis

section_content: Above we presented the LWM model for general internal nodes, and a matching boundary conditions model was also established. In order to debug the proposed model, the pipe system in Figure 2 is simulated and compared with MOC. The performance of the valve is shown in Figure 5. Figure 7 shows the numerical solution of the water hammer process produced in the condition that the valve is closed in 2. 2 s. As shown in the figure, the proposed model is in agreement with the MOC; indeed, they can almost produce the same numerical simulation. However, a small difference can be found in the phase of the wave. Obviously, the phase difference increases with time, because the proposed model considers the influence of the pipe flow velocity on the wave propagation. With a lower speed, the flow speed can be considered by the convective term. As shown in Figure 7, both LWM and MOC can give out reasonable maximum pressure, which is the key to water hammer prediction. 

Water 2018, 10, x FOR PEER REVIEW 9 of 17
section_name: Simulation and Comparison

section_content: Above we presented the LWM model for general internal nodes, and a matching boundary conditions model was also established. In order to debug the proposed model, the pipe system in Figure 2 is simulated and compared with MOC. The performance of the valve is shown in Figure 5. Figure 7 shows the numerical solution of the water hammer process produced in the condition that the valve is closed in 2. 2 s. As shown in the figure, the proposed model is in agreement with the MOC; indeed, they can almost produce the same numerical simulation. However, a small difference can be found in the phase of the wave. Obviously, the phase difference increases with time, because the proposed model considers the influence of the pipe flow velocity on the wave propagation. With a lower speed, the flow speed can be considered by the convective term. As shown in Figure 7, both LWM and MOC can give out reasonable maximum pressure, which is the key to water hammer prediction.
section_name: Simulation and Comparison

section_content: As seen in Figure 7, the LWM can simulate the actual maximum water hammer pressure. It may meet the common requirement for predicting and preventing water hammer overpressure. However, the numerical results heavily distorted the properties of the attenuation, compared to the experimental result. After the complete closing of the valve, the peak pressure of the water hammer wave decreases rapidly in the experiment, but the water hammer wave decreases slowly in the numerical simulation. Obviously, the attenuation is heavily underestimated in the numerical
section_name: Improvement to Approximate the Attenuation

section_content: As seen in Figure 7, the LWM can simulate the actual maximum water hammer pressure. It may meet the common requirement for predicting and preventing water hammer overpressure. However, the numerical results heavily distorted the properties of the attenuation, compared to the experimental result. After the complete closing of the valve, the peak pressure of the water hammer wave decreases rapidly in the experiment, but the water hammer wave decreases slowly in the numerical simulation. Obviously, the attenuation is heavily underestimated in the numerical simulation. In order to modify the deviation, an additional friction function is introduced according to the flow patterns. 

In Figure 7, the numerical model evaluates the unsteady friction resistance according to the experimental measurement shown in Figure 4, which denotes the steady friction in pipe flow. The model can simulate a water hammer well before the valve has closed completely, and during this period the pipe flow maintained positive velocity. However, after the valve is completely closed, the fluid will have a great impact on the pipeline, and then the steady friction cannot actually approximate the practical energy dissipation. The experiment can actually measure the attenuation process of the water hammer wave, but it cannot measure the impact flow rate and the actual energy dissipation. In order to evaluate the attenuation, an additional friction function is introduced to the numerical model. The friction term is proposed as follows

where f is the proposed additional function after the valve's complete closing, which can be described

. where t c is the valve closing time. 

Considering a, m, g and L as the main physical variables, combining the Darcy-Weisbach friction factor, the additional function can be described as

where λ is the Darcy-Weisbach friction factor. According to Rayleigh's method of dimensional analysis, the exponential function is proposed as

where C is a dimensionless coefficient that needs to be calibrated by experiment. After a complete valve closing, the flow is only impacted at low velocity. According to the flow patterns, the λ can separately be described as λ = 64ν vD in laminar model (LM) and λ = 8gn 2 3 √ R p in turbulent model (TM) [35]. 

Figure 8 shows that the flow velocity may cross laminar and turbulent zones during hydraulic transient. The compound model (CM) can be described as

where v c is the critical flow velocity between laminar and turbulent flow. 

where c v is the critical flow velocity between laminar and turbulent flow. In order to calibrate the dimensionless coefficient C, the positive and negative peaks are identified and analyzed in Figure 9. 

Water 2018, 10, x FOR PEER REVIEW 11 of 17

In order to calibrate the dimensionless coefficient C, the positive and negative peaks are identified and analyzed in Figure 9. Based on the identified peaks, the method of least squares is used to obtain the optimal C. The squared residuals between experimental and simulated peaks are defined as

Then, the sum of squared residuals is ( )

The optimal C can be determined when the sum of squared residuals is the least. In order to optimize the dimensionless coefficient C, the water hammer processes are simulated by various dimensionless coefficients and three kinds of friction models. The relationship between the sum of squared residuals and the dimensionless coefficient C in each model are shown in Figure 10. Table 1 shows the first 10th peaks of the optimized simulations in each model. The result shows that the CM model can give the least sum of squared residuals, in which the matched coefficient is about 8. 48. Then the additional function is written as

Figure 11 compares the water hammer simulations. Obviously, the improved model can more actually simulate the attenuation. Based on the identified peaks, the method of least squares is used to obtain the optimal C. The squared residuals between experimental and simulated peaks are defined as

Then, the sum of squared residuals is

The optimal C can be determined when the sum of squared residuals is the least. In order to optimize the dimensionless coefficient C, the water hammer processes are simulated by various dimensionless coefficients and three kinds of friction models. The relationship between the sum of squared residuals and the dimensionless coefficient C in each model are shown in Figure 10. Table 1 shows the first 10th peaks of the optimized simulations in each model. The result shows that the CM model can give the least sum of squared residuals, in which the matched coefficient is about 8. 48. Then the additional function is written as

Figure 11 compares the water hammer simulations. Obviously, the improved model can more actually simulate the attenuation. By introducing the additional friction function, the improved model can simulate water hammer more actually, not only on the maximum pressure, but also the attenuation of the water hammer waves. 

For validation, the improved model is applied to various air-water mixtures conditions as shown in Table 2. Changing the valve closing time and the air mass entrainment, various water hammer processes are produced. These water hammers are experimentally measured and numerically simulated. Although the actual air content is hard to measure, it can be estimated by the actual wave speed according to previous research [1]. In the numerical model, the air content is taken into account by adjusting m in Equation ( 26 ), the amount of substance of the gas in unit volume, which affects the energy loss degree greatly. This is inferred from the measured actual wave speed in experiment using Equation (19). As shown in Figure 12, the numerical simulation results using the improved LWM are in agreement with the experimental results both on the maximum pressure and the attenuation trends. By introducing the additional friction function, the improved model can simulate water hammer more actually, not only on the maximum pressure, but also the attenuation of the water hammer waves. 

For validation, the improved model is applied to various air-water mixtures conditions as shown in Table 2. Changing the valve closing time and the air mass entrainment, various water hammer processes are produced. These water hammers are experimentally measured and numerically simulated. Although the actual air content is hard to measure, it can be estimated by the actual wave speed according to previous research [1]. In the numerical model, the air content is taken into account by adjusting m in Equation ( 26 ), the amount of substance of the gas in unit volume, which affects the energy loss degree greatly. This is inferred from the measured actual wave speed in experiment using Equation (19). As shown in Figure 12, the numerical simulation results using the improved LWM are in agreement with the experimental results both on the maximum pressure and the attenuation trends.
section_name: Improvement to Approximate the Attenuation

section_name: Discussion

section_content: The LWM is widely used in computational fluid dynamics for fluid field simulation. However, it is uncommonly used in numerical water hammer-simulation for pipe systems. As an alternative approach, LWM can be applied to the water hammer hyperbolic partial differential equations system. In particular, the convective acceleration terms are completely considered in the model, although they only affect the water hammer slightly. Considering the applicability in numerical computation, the LWM can yield the water hammer maximum pressure as the MOC, which shows that the LWM can suit the water hammer equations. However, considering the water hammer experiment in the pipe, LWM may heavily distort the attenuation of the water hammer wave. Based on the experiment, an additional friction function is proposed to simulate the actual wave attenuation. Traditional numerical simulations greatly underestimate air resistance when the valve is closed completely. By introducing the additional resistance function, the proposed model can simulate the water hammer well, both on the maximum pressure and the attenuation process. In particular, it is difficult to reveal the mechanism of the resistance in the shock flow after completely closing valve. Many studies have reported the difference of the water hammer attenuation process in theory and experiment. Traditional models underestimate the resistance, which makes the fluctuation remain longer in the numerical result. Considering the restriction of the experimental scale of the pipe used in this research on the water hammer in gravitational pipe flow system with air entrainment, it may require a different friction coefficient in flow with air content after complete valve closing.
section_name: Discussion

section_content: Based on the Lax-Wendroff difference Method, a numerical water hammer model is established for gravitational flow pipe systems. The proposed model considers the influence of convective acceleration terms, and it also supports flexibly setting the time steps and space steps without the limitation of the characteristics in MOC. With a second-order accuracy, the model is validated by MOC and used in an experiment to predict the maximum transient pressures. The experiment shows that the water hammer processes in water supply systems with air entrainment differ much from those in one phase of pure water, especially on the attenuation processes. The improved model with a modified friction function using a compound model can simulate the water hammer processes well, both on the maximum pressures and the attenuation process. This research reveals that we used to underestimate the energy loss process of a water hammer in pipe systems, especially in conditions with air content. Proposed improved LWM can simulate water hammers much better in such conditions. As plenty of practical water supply systems are running with air content, this research can be helpful on predicting water hammers in actual engineering. 

In further research, the physical interaction of air entrainment in water hammer processes and the air-delivery process in a water hammer will be researched. To figure out the air inlet conditions and its physical interactions, a further experiment will be conducted using a tube with a visible interior to observe the air-delivery trail and its shape-changing process in water hammers.
section_name: Conclusions

section_content: Funding: This research work was funded by the National Natural Science Foundation of China (Grant Nos. 51779216, 51279175 ) and Zhejiang Provincial Natural Science Foundation of China (Grant No. LZ16E090001 ).