Finite Volume Method for Transient Pipe Flow with an Air Cushion Surge Chamber Considering Unsteady Friction and Experimental Validation

: In various water transmission systems such as long-distance water transfer projects and hydropower stations, accurate simulation of water hammer is extremely important for safe and stable operation and the realization of intelligent operations. Previous water hammer calculations usually consider only steady-state friction, underestimating the decay of transient pressure. A second-order Finite Volume Method (FVM) considering the effect of unsteady friction factor is developed to simulate the water hammer and the dynamic behavior of air cushion surge chamber in a water pipeline system, while an experimental pipe system is conducted to validate the proposed numerical model. Two unsteady friction models, Brunone and TVB models, were incorporated into the water hammer equations, which are solved by the MUSCL–Hancock method. One virtual boundary method was proposed to realize the FVM simulation of Air Cushion Surge Chamber. Comparisons with water hammer experimental results show that, while the steady friction model only accurately predicts the ﬁrst pressure peak, it seriously underestimates pressure attenuation in later stages. Incorporating an unsteady friction factor can better predict the entire pressure attenuation process; in particular, the TVB unsteady friction model more accurately reproduces the pressure peaks and the whole pressure oscillation periods. For water pipeline systems with an air cushion surge chamber, energy attenuation of the elastic pipe water hammer is primarily due to pipe friction and the air cushion. The experimental results with the air cushion surge chamber demonstrate that the proposed FVM model with the TVB unsteady friction model and the air chamber polytropic exponent near 1.0 can well reproduce the experimental pressure oscillations.


Introduction
Water hammer often occurs in various water pipe systems, including long-distance water transfer projects and hydropower stations. An abnormal pressure surge may lead to a pipe burst, so many water hammer protection measures are introduced to reduce the water hammer intensity. The air cushion surge chamber is a closed chamber that is partially filled with water and compressed air [1]. As compared to an open-type pressure regulating chamber, it is rarely restricted by geological or topographical features, and offers numerous advantages, such as shorter construction supply, lower excavation volume, cost-effectiveness, and minimal ecological impact [2]. The air cushion surge chamber is widely used in hydroelectric power plants for water hammer protection, ensuring the safe hydraulic operation of the water pipe system. In order to realize the safe and stable operation of water systems and the realization of intelligent operations, it is extremely important to accurately simulate the transient flow of pipe systems with an air cushion surge chamber.
The Method of Characteristics (MOC) is currently a widely used simulation tool for modeling the hydraulic transition process of hydroelectric power plants. However, there are some complicated situations, such as short pipes, T-junctions, and series pipes, in actual water delivery systems. When using MOC, interpolation or wave speed regulation is required, which reduces computational efficiency and accuracy and introduces computational errors. Moreover, most existing water hammer calculation models only adopt a steady friction model, implying that the friction inside the pipeline remains the same as the steady-state during the transient process. However, in actual transient processes, the friction inside the pipeline is influenced by multiple factors, leading to significantly different calculation results from the actual results. Additionally, they cannot accurately describe the waveform distortion and peak attenuation of pressure waves [3].
The FVM discretized the calculation area in the pipeline into independent control units, and solved the differential equations in each unit separately. Godunov et al. [4] proposed a numerical scheme for solving nonlinear Riemannian problems. This scheme is very suitable for approximating smooth solutions and discontinuous solutions. Therefore, in recent years, a large number of researchers have gradually begun to construct the Godunov scheme to solve the hydraulic transient water hammer problem.
Yazdi et al. [5] pointed out in 2007 that when calculating hydraulic transients, if the Courant number condition is not met, the second-order Godunov scheme is more stable than the MOC method. Bi Sheng et al. [6] adopted the Godunov scheme to solve the two-dimensional flow-transport equation. This model can simultaneously solve the water flux and transport flux, which is highly efficient for simulating the dynamic characteristics of water flow in complex terrain, and effectively eliminates problems such as excessive numerical damping and unstable oscillation caused by the convection term in numerical calculation. Zheng Jieheng et al. [7] used the Godunov scheme to study the hydraulic transients in sequential transmission pipelines. Based on the finite volume method, Zhao Yue et al. [8] proposed a treatment method with double virtual boundary to numerically simulate the phenomenon of water hammer and water column separation in pipelines. Zhou et al. [9][10][11] proposed a method to simulate a liquid column separation-bridging water hammer using a second-order GODUNOV scheme. Hu et al. [12] proposed the application of a second-order GODUNOV scheme to simulate non-pressurized flow.
Currently, there are two main unsteady friction models extensively utilized. These models are the weighted function model represented by the Zielke [13] and the empirical correction model represented by the Brunone [14]. According to the Zielke unsteady friction model, the instantaneous shear stress in the pipe due to transient flow is composed of a constant term and an additional term. The additional term uses the weighted function to account for the impact of historical velocity and acceleration on the current flow state. However, this method has a long calculation time and requires a large storage space. Subsequently, Zielke's model was simplified by Trikha [15], Vardy [16], and other scholars, resulting in a weighted function class unsteady friction model with higher computational efficiency. The Brunone unsteady friction model links unsteady friction with instantaneous local acceleration and convective acceleration and proposes a new dynamic friction model. Vitkovsky [17] improved the Brunone model by predicting the direction of water flow and wave propagation, as well as the effects of specific acceleration and deceleration stages.
To simulate more accurately the hydraulic transient process of the pressurized delivery pipeline system with an air cushion surge chamber, this paper introduces the second-order Godunov format of FVM during the calculation process, incorporating the Trikha-Vardy-Brunone (TVB) and Brunone unsteady friction models. One virtual boundary method was proposed to realize the FVM simulation of an air cushion surge chamber. An experimental pipe system was designed and conducted to validate the proposed models in simulating the water hammer and dynamic behavior of an air cushion surge chamber.
The novelty of the paper is that a second-order FVM considering the effect of unsteady friction factor is developed to simulate the dynamic behavior of the air cushion surge chamber in a water pipeline system, while an experimental pipe system is conducted to validate the proposed numerical model. Pressure damping and energy dissipation of transient flow in a water pipeline system with air cushion surge chamber are carefully investigated and modeled, which have not been well considered in previous work.

Water Hammer Control Equations
Equations of motion and continuity for water hammer are [18]: The matrix forms of the above two equations can be expressed as follows: x is the distance along the pipe axis coordinate. H is the piezometric head; V is the flow velocity in the pipe; a is the wave velocity of the water hammer; D is the inner diameter of the pipe; J Q and J u represent steady friction and unsteady friction J u ; and t is time.
If V = 0, the convection term can be ignored and the classical water hammer governing equation can be changed.
Brunone used instantaneous local acceleration and convective acceleration to represent the unsteady frictional component of dynamic frictional resistance, resulting in an empirically modified model. Based on Brunone's work, Vitkovsky later added the identification of the flow direction of water and obtained an improved instantaneous acceleration model with higher calculation accuracy. The specific model forms are as follows: (5) in which the Brunone coefficient of friction k = √ C * /2; C* is the shear attenuation constant and the value depends on the Reynolds number. When the water flow in the pipe is laminar flow, C* = 0.00476; when the water flow in the pipe is turbulent flow, The Zielke model simplified by Vardy and Brown (TVB model) is as follows where τ u is the unsteady shear stress; ρ is density of fluid; ν is kinematic viscosity of the fluid; N is number of cells along the pipeline; Y ai (t) is weighting function; and R is pipe radius; When the water flow in the pipe is laminar The computational domain is discretized along both the x-axis and t-axis using the finite volume method. This results in the formation of multiple computational control volumes or cells, as shown in Figure 1. Calculations are then performed on these volumes.
where τu is the unsteady shear stress; ρ is density of fluid; ν is kinematic viscosity of the fluid; N is number of cells along the pipeline; Yai(t) is weighting function; and R is pipe radius; When the water flow in the pipe is laminar The computational domain is discretized along both the x-axis and t-axis using the finite volume method. This results in the formation of multiple computational control volumes or cells, as shown in Figure 1. Calculations are then performed on these volumes. Integrate Equation (4) from the control surface i − 1/2 to control surface i + 1/2. Additionally, since the control variable uniformly and continuously changes, substitute = 1 ∫ +1/2 −1/2 to obtain the integration expression for control variable U: where: Fi+1/2 and Fi−1/2 are the fluxes at i + 1/2 and i − 1/2, respectively; Δt is the time step; Δx is the length of control volume; the superscript n denotes the current time step; and n + 1 denotes the subsequent time step. Integrate Equation (4) from the control surface i − 1/2 to control surface i + 1/2. Additionally, since the control variable uniformly and continuously changes, substitute U i = 1 ∆x i+1/2 i−1/2 udx to obtain the integration expression for control variable U: where: F i+1/2 and F i−1/2 are the fluxes at i + 1/2 and i − 1/2, respectively; ∆t is the time step; ∆x is the length of control volume; the superscript n denotes the current time step; and n + 1 denotes the subsequent time step.

Control Equations of Air Cushion Surge Chamber
The continuity equation for node P at the bottom of the surge chamber is [18] where Q T is the flow rate at the upstream pipe outlet of the surge chamber; Q is the flow rate of the downstream pipeline inlet; and Q S is the flow into or out of the surge chamber ( Figure 2).

Control Equations of Air Cushion Surge Chamber
The continuity equation for node P at the bottom of the surge chamber is [18] = + where QT is the flow rate at the upstream pipe outlet of the surge chamber; Q is the flow rate of the downstream pipeline inlet; and QS is the flow into or out of the surge chamber ( Figure 2). In order to determine the pressure tube head at the bottom point P of the surge chamber, it is necessary to use the feature compatibility equation of the last calculation section of the upstream pipeline and the first calculation section of the downstream pipeline in the surge chamber.
where Hp is the piezometric head at the bottom of air cushion surge chamber.
Equations (10) and (11) are substituted into continuity Equation (9) at point P, and the variables QT and Q are eliminated where CP1, BP1, CM2, and BM2 are known variables at time t.
in which B is a function of the physical properties of the fluid and the pipeline, often called the pipeline characteristic impedance, and B = a/gA; R is the pipeline resistance coefficient R = f∆x/(2gDA 2 ); f is the Darcy-Weisbach friction factor; D is pipe diameter; and A is the cross-section area. Neglecting the water flow inertia and frictional head loss in the air cushion surge chamber, a relationship can be found between the pressure at the center point of the bottom and the water level in the air cushion surge chamber: In order to determine the pressure tube head at the bottom point P of the surge chamber, it is necessary to use the feature compatibility equation of the last calculation section of the upstream pipeline and the first calculation section of the downstream pipeline in the surge chamber.
C + : where H p is the piezometric head at the bottom of air cushion surge chamber.
Equations (10) and (11) are substituted into continuity Equation (9) at point P, and the variables Q T and Q are eliminated where C P1 , B P1 , C M2, and B M2 are known variables at time t.
in which B is a function of the physical properties of the fluid and the pipeline, often called the pipeline characteristic impedance, and B = a/gA; R is the pipeline resistance coefficient R = f ∆x/(2gDA 2 ); f is the Darcy-Weisbach friction factor; D is pipe diameter; and A is the cross-section area.
Neglecting the water flow inertia and frictional head loss in the air cushion surge chamber, a relationship can be found between the pressure at the center point of the bottom and the water level in the air cushion surge chamber: where H a is the absolute head equal to the gauge plus barometric pressure heads; Z s is the elevation of the air-water interface in air chamber; H atm is the absolute barometric pressure head; R s is the head loss coefficient of the impedance hole of the air chamber; and A s is the cross-section area of the air chamber.
The air was assumed to follow the reversible polytropic relation [18] PV n a = P 0 V n 0 = Constant (18) where V a is the volume of the air within the air chamber, and n is the polytropic exponent. After combining Equations (12)- (18), the flow rate, pressure head, and water level at the air chamber can be obtained.

Computation of Flux Term
The Riemann problem-solving method can be applied to obtain the flux values at the i + 1/2 and i + 1/2 boundary interfaces. The average value of the control variable U on the left side of the i + 1/2 interface is represented by U n L , and the average value of the control variable U on the right side of the i + 1/2 interface is represented by U n R .
The flux values at i + 1/2 can be calculated by: The MUSCL-Hancock method for second-order linear reconstruction is used to realize the second-order accuracy of computation results.
First step: Data Reconstruction.
Second step: Advance time calculation Third step: Solve the Riemann problem. To compute intercell flux f i+1/2 , the conventional Riemann problem with data can be calculated by By inserting the solved Equations (24) and (25) into the Equation (20), the fluxes in the second-order Godunov scheme at the boundary of each element can be obtained.

Time Integral
To advance the solution of the Godunov flux calculation at the n + 1 time step with second-order accuracy, time integration of Equation (8) is necessary. In the absence of friction, the following formula can be derived: The second-order Runge-Kutta method is used to obtain explicit results with a secondorder calculation accuracy when pipe friction is taken into account. The specific calculation process is shown as follows: First step: Second step: Last step:

Boundary Condition
The virtual boundary method is used to process the boundary. To implement this method, virtual cells numbered 0, −1, N + 1, and N + 2 are added to the upstream and downstream boundaries, respectively. The condition that the upstream boundary is and the downstream boundary is U −1 = U 0 = U 1/2 , U N+1 = U N+2 = U N+1/2 enable the determination of the flux value at the boundary using the Riemann invariant equation. Solving the Riemann invariant equation at the upstream boundary of the reservoir yields [10]: where V n 1 and H n 1 are the velocity and pressure head of the first control volume, and V 1/2 and H 1/2 are the velocity and pressure head at 1/2 interface at the upstream boundary.
Solving the Riemann invariant equation at the upstream boundary of the reservoir yields [10]: where V n N and H n N are the velocity and pressure head of the last control volume, and V N+1/2 and H N+1/2 are the velocity and pressure head at N + 1/2 interface at the downstream boundary.
The physical variable values at each transient moment of the air cushion surge chamber can be determined by calculating the physical variable values of the virtual units at the upstream and downstream sides of the chamber based on the control equation. Thus, in combination with the Riemann invariant equation, C P1 , B P1 , C M2 , and B M2 are known variables at time t, which can conclude that: B P1 = a gA 1 (33) where V n Nu and H n Nu represent the flow rate and the water head of the last control domain of the air cushion surge chamber on the upstream side pipeline, while V n 1d and H n 1d , respectively, represent the flow rate and water head of the first control domain of the air cushion surge chamber on the downstream side pipeline. A 1 and A 2 are the pipe cross-section areas at the upstream and downstream of the air chamber.
By substituting Equations (32)-(35) into Equation (8), the relation between the water level in the piezometric head at the bottom of the air cushion surge chamber and the flow rate into the air cushion surge chamber under virtual boundary conditions can be solved.

Experimental Setup
An experimental pipe system was designed and conducted to validate the proposed models in simulating the water hammer and dynamic behavior of the air cushion surge chamber. Figure 3 displays the experimental setup. The system consists of an upstream reservoir, variable frequency pump, constant pressure tank, upstream pipe, upstream electromagnetic flowmeter, and 1# ball valve. After the 1# ball valve, the pipeline serves as a return pipe. The water hammer experimental pipeline is 582 m long and is made of copper tubes, which have a wall thickness of 2 mm and an inner diameter of 21 mm. Along the pipeline, one 1/4 ball valve (1# ball valve) is arranged at the end, which is 582 m away from the upstream constant pressure tank. Additionally, five pressure sensors have been where and represent the flow rate and the water head of the last control domain of the air cushion surge chamber on the upstream side pipeline, while 1 and 1 , respectively, represent the flow rate and water head of the first control domain of the air cushion surge chamber on the downstream side pipeline. A1 and A2 are the pipe crosssection areas at the upstream and downstream of the air chamber.
By substituting Equations (32)-(35) into Equation (8), the relation between the water level in the piezometric head at the bottom of the air cushion surge chamber and the flow rate into the air cushion surge chamber under virtual boundary conditions can be solved.

Experimental Setup
An experimental pipe system was designed and conducted to validate the proposed models in simulating the water hammer and dynamic behavior of the air cushion surge chamber. Figure 3 displays the experimental setup. The system consists of an upstream reservoir, variable frequency pump, constant pressure tank, upstream pipe, upstream electromagnetic flowmeter, and 1# ball valve. After the 1# ball valve, the pipeline serves as a return pipe. The water hammer experimental pipeline is 582 m long and is made of copper tubes, which have a wall thickness of 2 mm and an inner diameter of 21 mm. Along the pipeline, one 1/4 ball valve (1# ball valve) is arranged at the end, which is 582 m away from the upstream constant pressure tank. Additionally, five pressure sensors have been

Water Hammer Problem in a Simple Reservoir-Pipe-Valve System
Four experimental conditions were conducted. The relevant parameters are shown in Table 1. H0 is constant pressure head in the upstream pressure tank.

Water Hammer Problem in a Simple Reservoir-Pipe-Valve System
Four experimental conditions were conducted. The relevant parameters are shown in Table 1. H 0 is constant pressure head in the upstream pressure tank. The wave velocity of the water hammer in this experiment is calculated according to the experimental data measured by the pressure sensor. According to the pressure experimental data, the time difference between any two adjacent wave peaks is recorded as 2T, and then, according to the pipeline length L, water hammer wave a = 2L/T can be obtained. In addition, many factors can affect the water hammer wave velocity. In order to eliminate the interference, repeated tests were carried out on different experimental conditions and many times. Finally, the average value of multiple tests was taken as the value of the water hammer wave velocity in the subsequent numerical simulation. The wave velocity of the water hammer measured by the experiment is between 1260 m/s and 1360 m/s. In the numerical simulations, the wave speed is a = 1290 m/s and the experimentally observed range of the average resistance coefficient for the pipeline system under stable flow conditions is 0.0312~0.0356. Valve closing times ranged between 0.012 s and 0.026 s. The experimental results show that the valve closing time is much less than half of the pressure fluctuation period. Therefore, it is believed that the valve is closed instantaneously.
As shown in Figure 4, Cases 1 and 2 can produce a similar trend for experimental pressure oscillations, although with different reference values. This is because Cases 1 and 2 have the same Reynolds number Re = 1284, but different driving pressure heads (H 0 = 31 and 57). When the initial steady condition is under laminar flow in Case 1 and Case 2, the steady friction water hammer model can only accurately predict the first pressure peak, but fail to calculate well the subsequent pressure oscillations. The reason should be that, during the fast transient flow event, the pressure damping is mainly caused by the dynamic shear force near the pipe wall. However, the traditional water hammer model assumes that the dynamic shear force coefficient is constant. In order to verify this, the unsteady friction factor is considered in the water hammer events. Figure 4a,b show that, compared to the steady friction water hammer model, the TVB and Brunone unsteady friction models can provide the much better simulation results. Meanwhile, it can be found that the TVB model shows the highest consistency with experimental results, while the Brunone unsteady friction model still produces the differences in simulating the later pressure peaks and the pressure oscillations. The main reason is that, as described above, the Brunone unsteady friction model is an empirically modified model, and the TVB model is a mechanism model which is derived from the physical equations.
Cases 1 and 2 exhibit laminar flow with Re = 1284, while cases 3 and 4 demonstrate low Reynolds number turbulence with Re = 4334.
The wave velocity of the water hammer in this experiment is calculated according to the experimental data measured by the pressure sensor. According to the pressure experimental data, the time difference between any two adjacent wave peaks is recorded as 2T, and then, according to the pipeline length L, water hammer wave a = 2L/T can be obtained. In addition, many factors can affect the water hammer wave velocity. In order to eliminate the interference, repeated tests were carried out on different experimental conditions and many times. Finally, the average value of multiple tests was taken as the value of the water hammer wave velocity in the subsequent numerical simulation. The wave velocity of the water hammer measured by the experiment is between 1260 m/s and 1360 m/s. In the numerical simulations, the wave speed is a = 1290 m/s and the experimentally observed range of the average resistance coefficient for the pipeline system under stable flow conditions is 0.0312~0.0356. Valve closing times ranged between 0.012 s and 0.026 s. The experimental results show that the valve closing time is much less than half of the pressure fluctuation period. Therefore, it is believed that the valve is closed instantaneously.
As shown in Figure 4, Cases 1 and 2 can produce a similar trend for experimental pressure oscillations, although with different reference values. This is because Cases 1 and 2 have the same Reynolds number Re = 1284, but different driving pressure heads (H0 = 31 and 57). When the initial steady condition is under laminar flow in Case 1 and Case 2, the steady friction water hammer model can only accurately predict the first pressure peak, but fail to calculate well the subsequent pressure oscillations. The reason should be that, during the fast transient flow event, the pressure damping is mainly caused by the dynamic shear force near the pipe wall. However, the traditional water hammer model assumes that the dynamic shear force coefficient is constant. In order to verify this, the unsteady friction factor is considered in the water hammer events. Figure 4a,b show that, compared to the steady friction water hammer model, the TVB and Brunone unsteady friction models can provide the much better simulation results. Meanwhile, it can be found that the TVB model shows the highest consistency with experimental results, while the Brunone unsteady friction model still produces the differences in simulating the later pressure peaks and the pressure oscillations. The main reason is that, as described above, the Brunone unsteady friction model is an empirically modified model, and the TVB model is a mechanism model which is derived from the physical equations.     mates the pressure damping in the later pressure oscillations. In contrast to the steady friction model, the TVB and Brunone unsteady friction models can accurately predict the entire pressure attenuation process. Among these models, the TVB model can most accurately reproduce the experimental results.
It can also be found from Figures 4 and 5 that, as the Reynolds number increases, the steady friction water hammer model can give better simulation results.
Overall, the TVB model can accurately reproduce the experimental results, regardless of the presence of laminar or turbulent flow, and is recommended to simulate the water hammer events.  Figure 1 demonstrates that, by opening ball valve #3 on the connection pipe of the air cushion surge chamber, the experimental device transforms into a pressurized pipeline water supply system equipped with an air cushion surge chamber. Furthermore, a pressure sensor PT-5# is placed on top of the air cushion chamber to monitor the gaseous pressure inside it, and the steady-state gas pressure of the cushion surge chamber is also measured by this pressure sensor. The velocity of pressure wave in the air cushion surge chamber experiment is consistent with that of the above water hammer experiment without air cushion surge chamber, i.e., the wave speed is a = 1290 m/s. Table 2 presents the experimental conditions of the air cushion surge chamber. The gas polytropic index, n, is varied to 1.0, 1.2 and 1.4, where n = 1.0 corresponds to an isothermal process, n = 1.4 corresponds to an adiabatic process, and n = 1.2 for the two in between, respectively, and the simulation results are depicted in Figure 6.  It can also be found from Figures 4 and 5 that, as the Reynolds number increases, the steady friction water hammer model can give better simulation results.

Water Hammer Problem with Air Cushion Surge Chamber
Overall, the TVB model can accurately reproduce the experimental results, regardless of the presence of laminar or turbulent flow, and is recommended to simulate the water hammer events. Figure 1 demonstrates that, by opening ball valve #3 on the connection pipe of the air cushion surge chamber, the experimental device transforms into a pressurized pipeline water supply system equipped with an air cushion surge chamber. Furthermore, a pressure sensor PT-5# is placed on top of the air cushion chamber to monitor the gaseous pressure inside it, and the steady-state gas pressure of the cushion surge chamber is also measured by this pressure sensor. The velocity of pressure wave in the air cushion surge chamber experiment is consistent with that of the above water hammer experiment without air cushion surge chamber, i.e., the wave speed is a = 1290 m/s. Table 2 presents the experimental conditions of the air cushion surge chamber. The gas polytropic index, n, is varied to 1.0, 1.2 and 1.4, where n = 1.0 corresponds to an isothermal process, n = 1.4 corresponds to an adiabatic process, and n = 1.2 for the two in between, respectively, and the simulation results are depicted in Figure 6. one is energy loss due to heat transfer during air compression and expansion of air chamber, and another is hydraulic loss caused by the pipe friction; (2) results in Figure 6 only consider the effect of the steady friction factor, neglecting the effect of unsteady friction. Therefore, in the following section, the effect of unsteady friction will be included in the numerical simulation to enhance the numerical accuracy, as well as to verify the above explanation.

Water Hammer Problem with Air Cushion Surge Chamber
(a) (b)  Figure 7 gives the results calculated by numerical models (steady and unsteady friction water hammer model) and experimental results in Case 1 (laminar flow) and Case 3 (turbulent flow). According to Section 4.2, the TVB unsteady friction model performs better simulations for water hammer pressure of the pipeline in this experimental system. Therefore, the TVB unsteady friction model is adopted for this section. Additionally, examination of the simulation results from steady friction models indicated that n = 1.0 provides the optimal simulation effect. Consequently, when using unsteady friction simulation, n is set to 1.0.
As shown in Figure 7, the results of Case 1 (laminar flow) and Case 3 (turbulent flow) both show that introduction of the unsteady friction model causes the first peak pressure to slightly increase, compared to the steady friction model. This is because air in the chamber experiences an expansion and deceleration process when the valve is quickly closed, and the unsteady friction model suppresses its deceleration, which results in the first pressure peak value increasing. However, in comparison to the steady friction model, the peak and cycle of the first period align more closely with experimental data. The peak value and cycle of the pressure decay process also exhibits better agreement with experimental data, indicating that the use of the unsteady friction model can more accurately simulate the hydraulic transient process of the pressurized pipeline water supply system. The difference between the peak pressure of the numerical simulation and the experimental results becomes larger, which is caused by the absence of wall heat exchange in the mathematical model. In future work, we will also take into account the energy loss caused by the heat exchange of the tube wall. (2) in the existing numerical models, the air was assumed to follow the reversible polytropic relation [18]; (3) the gas polytropic index, n, is varied to 1.0, 1.2 and 1.4, where n = 1.0 corresponds to an isothermal process, n = 1.4 corresponds to an adiabatic process, and n = 1.2 for middle process; (3) in this work, the heat transfer is serious, which makes the thermal process close to the isothermal process (n = 1.0).
Meanwhile, why, even for n = 1.0, do the calculated pressure oscillations attenuate lower than the measured ones? The reasons are that: (1) in the transient event in this work, the pressure damping (namely energy dissipation) is attributed to two aspects, in which one is energy loss due to heat transfer during air compression and expansion of air chamber, and another is hydraulic loss caused by the pipe friction; (2) results in Figure 6 only consider the effect of the steady friction factor, neglecting the effect of unsteady friction.
Therefore, in the following section, the effect of unsteady friction will be included in the numerical simulation to enhance the numerical accuracy, as well as to verify the above explanation. Figure 7 gives the results calculated by numerical models (steady and unsteady friction water hammer model) and experimental results in Case 1 (laminar flow) and Case 3 (turbulent flow). According to Section 4.2, the TVB unsteady friction model performs better simulations for water hammer pressure of the pipeline in this experimental system. Therefore, the TVB unsteady friction model is adopted for this section. Additionally, examination of the simulation results from steady friction models indicated that n = 1.0 provides the optimal simulation effect. Consequently, when using unsteady friction simulation, n is set to 1.0.
As shown in Figure 7, the results of Case 1 (laminar flow) and Case 3 (turbulent flow) both show that introduction of the unsteady friction model causes the first peak pressure to slightly increase, compared to the steady friction model. This is because air in the chamber experiences an expansion and deceleration process when the valve is quickly closed, and the unsteady friction model suppresses its deceleration, which results in the first pressure peak value increasing. However, in comparison to the steady friction model, the peak and cycle of the first period align more closely with experimental data. The peak value and cycle of the pressure decay process also exhibits better agreement with experimental data, indicating that the use of the unsteady friction model can more accurately simulate the hydraulic transient process of the pressurized pipeline water supply system. The difference between the peak pressure of the numerical simulation and the experimental results becomes larger, which is caused by the absence of wall heat exchange in the mathematical model. In future work, we will also take into account the energy loss caused by the heat exchange of the tube wall.

Conclusions
A second-order FVM considering the effect of unsteady friction factor is developed to simulate the water hammer and the dynamic behavior of air cushion surge chamber in a water pipeline system, while an experimental pipe system is conducted to validate the proposed numerical model. Two unsteady friction models, Brunone and TVB models, were incorporated into the water hammer equations, and the virtual boundary method was proposed to realize the FVM simulation of the air cushion surge chamber.
Comparisons with water hammer experimental results show that, while the steady friction model only accurately predicts the first pressure peak, it seriously underestimates pressure attenuation in later stages. Incorporating the unsteady friction factor can better predict the entire pressure attenuation process; in particular, the TVB unsteady friction model more accurately reproduces the pressure peaks and the whole pressure oscillation periods.
For water pipeline systems with an air cushion surge chamber, energy attenuation of the elastic pipe water hammer is primarily due to pipe friction and air cushion. The experimental results for the air cushion surge chamber demonstrate that the proposed FVM model with TVB unsteady friction model and the air chamber polytropic exponent near 1.0 can well reproduce the experimental pressure oscillations. However, the difference between the peak pressures of the numerical simulation and the experimental results becomes larger, which is caused by the absence of wall heat exchange in the mathematical model. In the future, we will also take into account the energy loss caused by the heat exchange of the tube wall.  Data Availability Statement: Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
A second-order FVM considering the effect of unsteady friction factor is developed to simulate the water hammer and the dynamic behavior of air cushion surge chamber in a water pipeline system, while an experimental pipe system is conducted to validate the proposed numerical model. Two unsteady friction models, Brunone and TVB models, were incorporated into the water hammer equations, and the virtual boundary method was proposed to realize the FVM simulation of the air cushion surge chamber.
Comparisons with water hammer experimental results show that, while the steady friction model only accurately predicts the first pressure peak, it seriously underestimates pressure attenuation in later stages. Incorporating the unsteady friction factor can better predict the entire pressure attenuation process; in particular, the TVB unsteady friction model more accurately reproduces the pressure peaks and the whole pressure oscillation periods.
For water pipeline systems with an air cushion surge chamber, energy attenuation of the elastic pipe water hammer is primarily due to pipe friction and air cushion. The experimental results for the air cushion surge chamber demonstrate that the proposed FVM model with TVB unsteady friction model and the air chamber polytropic exponent near 1.0 can well reproduce the experimental pressure oscillations. However, the difference between the peak pressures of the numerical simulation and the experimental results becomes larger, which is caused by the absence of wall heat exchange in the mathematical model. In the future, we will also take into account the energy loss caused by the heat exchange of the tube wall.  Data Availability Statement: Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.