A Simulation Calculation Method of a Water Hammer with Multpoint Collapsing

The traditional discrete vapor cavity model (DVCM) is widely used in water hammer simulation. Water column separation in pipelines is usually predicted with this model. Nevertheless, the main weaknesses of this model consist of numerical instability and nonconvergence. Regarding the weaknesses of the traditional model, this paper discusses an improved method. The new method uses a new water hammer velocity formula, a new cavity model, and a floating grid method. Through simulations to test the effects of the new model, an experimental platform can be established to realize a water hammer with multipoint collapsing. The numerical simulation was programmed in C++ and the test was carried out with an actual pipeline model built in the laboratory. After certain modelling and calibration, the parameters in the simulation calculation were consistent with the measured parameters in the test. The numerical simulation results were compared with the experimental results. For the hydraulic transient system with multipoint collapsing, the superposition effect of the wave crest of the pseudo-water hammer in the traditional calculation model was obvious. The pressure of the water hammer in the simulation calculation was significantly higher than the actual value and the convergence effect of the water hammer wave was not good. Compared with the results of the traditional model, the simulation results of the new model were closer to the measured values. Therefore, the new model has better numerical solution accuracy, stability, and convergence, which is worth further study and promotion.


Introduction
Water hammers with multipoint collapsing are the main cause of pipeline burst accidents in water transportation projects. When the fluid pressure in pipelines temporarily drops to the vapor pressure of the fluid, a vaporization steam chamber is formed. When the steam chamber collapses, the water hammer pressure boost is tens of times that of an ordinary water hammer. The formation mechanism and simulation calculation of water hammers with multipoint collapsing are popular research topics [1]. Previous studies have focused on water hammers with one-point collapse, but the superposition or reflection of the pressure waves of multipoint collapse has seldom been reported. The calculation of multipoint separated liquid columns in hydraulic transient processes remains a challenge [2].
Water column separation may lead to severe industrial accidents, which is why engineers should estimate its effects in pipelines. Some scholars have studied water hammers with one-point collapsing. Valuable information about various models used to simulate water column separation can be found in Hatcher [3] and Adamowski [4]. The traditional discrete vapor cavity model (DVCM) is used for most engineering transient simulation software packages because it is easily implemented and reproduces many features of the physical events of column separation. The traditional DVCM can provide accurate results if the number of reaches is restricted; unrealistic spikes may be generated due to multicavity collapses. Many scholars have studied gas release. Autrique et al. studied the theoretical model of the water hammer of cavities collapsing with interruption and experimentally verified the model [5]. Soares et al. pointed out that when the gas release was ignored, the liquid viscosity significantly affected the friction [6]. Kranenburg et al. [7] proposed dividing the flow into three areas, namely, the water column separation area, cavity area, and conventional water hammer area, and indicated that the flow was accompanied by gas release and water column separation. Apollonio et al. considered the phase transition of gas release and liquid vaporization at the separation position of the water column and believed that the wave velocity of water hammer varied with the gas content [8]. A large number of research results show that the gas holdup in liquid will affect the wave velocity of the water hammer, which is constantly changing in the process of hydraulic transience. Many scholars are interested in the study of cavity models. Daily et al. [9] analyzed and studied the basic equations of transient gas-liquid two-phase flow and proposed the basic assumptions of the models as follows. First, in the one-dimensional model, it is assumed that the velocity, pressure, and water density are the same at any position in the pipeline section. In other words, the average values of the velocity, pressure, and water density in the pipeline section can be adopted [10]. Second, the momentum can be ignored in the gas-liquid two-phase exchange. Third, the relative motion between the gas and liquid phases can be ignored. Fourth, compared with the liquid phase, the mass of the gas phase can be ignored. Fifth, the bubbles are uniformly distributed in the liquid pipeline. According to the above assumptions, a uniform bubble distribution model can be derived. The research results of Daily et al. are the theoretical basis of the traditional DVCM model. Kiyama et al. studied the water hammer model of cavities collapsing with multiple interruptions in pipeline systems through theoretical simulations and experiments [11]. The bubble uniform distribution model of the gas-liquid two-phase flow mixture is based on the single-phase mixture. Bubbles are assumed to be uniformly distributed in the flowing liquid, whereas the surface tension is not considered. The relative motion between the two phases is ignored and the pressures inside and outside bubbles are the same. Therefore, the velocities of the gas and liquid phases are the same all the time and are equal to the mixing speed, and the mass of the gas phase can be neglected compared to the mass of the liquid phase. However, when bubbles are not evenly distributed in water, the model can be used to calculate the velocity of the water hammer with gas, but the calculation results contain significant errors. It is a long-term research topic for scholars to establish appropriate calculation methods. Kranenburg et al. [6,12] indicated that the finite difference method allowed similar accuracy to that of the characteristic method and could avoid the limitation of the time step, but it was difficult to guarantee its convergence rate. With the finite element method, Weyler et al. [13] and Watt et al. [14] found that the calculation result of the pressure jump of transient flow had difficulty converging. Amara et al. [15] used the Mac Cormack finite difference method to solve the water hammer equation. H Kim, S Kim et al. solved the control equation of the quasi two-dimensional model by the mixed scheme of the characteristic method and the finite difference method, taking into account the relationship between head reduction and cavitation [16,17]. S Malavasi et al. studied the dependence of incipient cavitation number on the most relevant parameters through experiments, and gave an estimation formula of incipient cavitation number [18]. A Ghodhbani et al. used the characteristic method and wave velocity adjustment method in the calculation.
In contrast, the traditional discrete vapor cavity model (DVCM) is the main method of predicting the water hammer process with liquid column separation since it has the advantages of simple operation and simple convergence conditions [19]. The basic theory of this method is that the steam region is concentrated at the calculated cross-section of the pipeline. When the pressure in the region drops to or below the steam pressure at the corresponding temperature, a cavity occurs. The fluid pressure in the cross-section with the cavity is equal to the saturated vapor pressure and remains constant. The cavity is filled with steam, thus leading to discontinuity in the original fluid, which is divided into two independent and continuous fluid regions. The water hammer wave propagates with a constant wave velocity only in the continuous fluid. It can be understood theoretically that only the vaporization of fluid is considered in the cavity, whereas the release of dissolved gas is not considered when the pressure is reduced. Even though the actual cavity area is composed of steam and released gas, when the fluid pressure is lower than the saturated steam pressure, the gas release in the fluid is a slow process. Therefore, it seems reasonable to ignore the gas release to simplify the model. However, when the fluid pressure in the pipeline decreases, gas release and liquid vaporization occur and the gas content in the fluid varies. The water hammer wave speed is greatly affected by the gas content in the fluid, and the water hammer wave speed directly affects the pressure boosting water hammer. Therefore, the influence of the gas content on the water hammer wave speed cannot be ignored; otherwise, the accuracy of the calculation results is greatly affected [17,20]. The main assumptions of the traditional discrete vapor cavity model (DVCM) are analyzed below to demonstrate its disadvantages [20][21][22]: The wave velocity of the water hammer is considered to be constant over the whole hydraulic transient process. The wave velocity of water hammer in a steady state is not in line with the actual conditions. In the process of hydraulic transition, the wave velocity of the water hammer is constantly changing. The change in pressure leads to a change in the gas content of the liquid, which affects the change in the water hammer wave velocity. According to the formula of Joukowski, the wave velocity of the water hammer directly affects the pressure boosting of the water hammer.

2.
Usually, with the node taken as the research object, it is assumed that when the flow into the node is not equal to that out of the node, a cavity is generated immediately. This assumption is not in line with the actual conditions. The fluid is elastic, and when the water column is stretched within the range of fluid elasticity, no cavity is generated.

3.
It is assumed that the steam cavity fills the whole cross-section of the pipeline without considering the length of the cavity. The model is a one-dimensional model describing the flow of an unstable pipeline. This assumption is not in line with the actual conditions. The parameters of the cavity boundary are simply replaced with the calculation parameters of fixed nodes. The replaced results are not consistent with the actual conditions. In view of the shortcomings of the traditional DVCM model, the improved method was studied.
A new wave velocity formula for a water hammer and a cavity model considering liquid elasticity were established, and the calculation was carried out by the floating grid method. The disadvantage of the traditional DVCM is the existence of unrealistic oscillations that suddenly increase with the number of reaches (nonconvergence). The improved model should be able to overcome this problem. The simulation results show that the new model is superior to the traditional DVCM model in accuracy, stability, and convergence.

Wave Velocity of the Water Hammer
In the process of hydraulic transition, the wave velocity of the water hammer is constantly changing. The change in pressure in the transient process leads to the change in the gas content of the fluid, which affects the change in the water hammer wave velocity. According to the formula of Joukowski, the wave velocity of the water hammer directly affects the pressure boosting of the water hammer. The accuracy of the water hammer wave velocity formula largely determines the accuracy of the simulation calculation.
According to traditional theories of water hammer, the formula for direct water hammer pressure in a horizontal pipeline can be derived based on the momentum theorem [23]: where a is the water-hammer wave speed (m/s), ρ is the density of the fluid(kg/m 3 ), K is the bulk modulus of the fluid (Pa), E is the Young's modulus of elasticity of the pipe material (Pa), e is pipe wall thickness (mm), and D is the inner pipe diameter (mm).
Equation (1) is widely used to calculate the wave velocity of a water hammer. The propagation velocity of the pressure wave is an important parameter in a hydraulic transition process. In the transient process without gas release and liquid column separation, the wave velocity is a fixed value. In a hydraulic transient process, when the pressure declines to a certain value the gas nucleus expands to form smaller bubbles and the dissolved gas is released. In such a process, the liquid compressibility varies significantly and the propagation velocity of the pressure wave is reduced. If the variation in the wave velocity of the water hammer is not considered in the simulation computation, a significant deviation is generated in the calculation results.
Kobori et al. [21] conducted experiments and theoretical discussion on the relationship between the gas content and wave velocity and analyzed the influences of gas content on wave velocity. Peassall et al. [24] obtained the theoretical formula of the relationship between gas content and wave velocity and drew a chart. The previous study was based on the steady-state process, but the data in the steady-state process were not applicable to describe the transient process and the pressure change in the transient process led to the unsteady gas content in the pipeline. Wylie et al. believed that the gas release process was a one-way process [25]. Adamkowski et al. analyzed the gas release phenomenon and indicated that the gas release might occur anywhere in the pipeline [4].
The influences of the wave velocity of the water hammer on the pressure boost have been widely studied. Various calculation formulas of the wave velocity of a water hammer have been proposed based on the consideration of released gas as follows.
Later, Crepo et al. [26] derived the calculation formula of the wave velocity of a water hammer below: where γ is gaseous polytropic index, P is the liquid pressure (Pa), α is the gas content in unit volume liquid, ρ e is the density of liquid phase (kg/m 3 ), and ρ g is the density of gas phase (kg/m 3 ). For the convenience of comparison, kg is used to substitute for γP(kg = γP) and then Equation (3) is obtained as Henry et al. [27] proposed the calculation formula of the wave velocity of a water hammer below: k g . Then, substituting these terms into Equation (4) gives: When ρ g /ρ e << 1, (a g /a e ) 2 << 1 and Equation (5) can be simplified as [8]: Keiji Miyazaki Yoichi Fujiion [28] obtained the calculation formula of the wave velocity of a water hammer below: a = 1 where dρ g dP = ρ g K g and dρ e dP = ρ e K e . Substituting these terms into Equation (7) gives: Martin, Padmanabahn, and Wiggert [29] proposed the calculation formula of the wave velocity of a water hammer below: Equation (9) is the calculation formula of the wave velocity derived from the three-group characteristic line equation [10], where u is the nonuniformity coefficient of the pipe wall material, K g , K e are the elastic moduli of the gas phase, and liquid phase, D is the pipe diameter (m).
Tullis, J.P., V.L. Streeter, and E.B. Wylie [30] proposed the following formula for the wave velocity of a water hammer: Since where V is equal to the gas proportion in unit volume. Substituting α = MPT P into Equation (10) gives: Wylie, E.B et al. [31] gave the calculation formula of the wave velocity of the water hammer below: where , and H is the piezometer height in the pipes.
After substituting the formulas mentioned above into Equation (13), based on P = ρgH, we obtained The proposal of water hammer wave velocity is as follows: The representative six calculation formulas of water hammer wave velocity in previous studies are summarized below. The six formulas are compared and analyzed, and their applicable conditions are also provided. In addition, a new calculation formula for the water hammer wave velocity is derived.
Equation (8) neglects the influence of pipe wall deformation on the wave velocity, whereas Equation (9) does not consider the influence of gas elasticity. Therefore, Equation (8) can be corrected as: where α is the gas content in unit volume liquid; ρ e is the density of liquid phase (kg/m 3 ), ρ g is the density of gas phase (kg/m 3 ), and K g , K e are the elastic moduli of the gas phase and liquid phase, respectively. D is the inner pipe diameter (m), E is the pipe wall Young's modulus of elasticity (Pa), e is the pipe wall thickness (mm), u is the nonuniformity coefficient of the pipe wall material.
Where if α = 0, the wave velocity of a single-phase fluid can be expressed as: In Equation (15), if α = 1, the wave velocity of a single-phase fluid can be expressed as: Thus, in Equation (15), α K g , 1−α K e , and Du Ee indicate the influences of gas deformation, liquid deformation, and pipe wall deformation on wave velocity in the two-phase gas-liquid fluid, respectively. The abovementioned six wave velocity formulas proposed by Henry et al. are all simplifications of the unified wave velocity formula Equation (15) under certain conditions. However, these formulas have some application limitations. Therefore, the applicable conditions of these equations should be considered.
The applicable conditions of the above six wave velocity formulas are discussed below. Compared with Equations (3) and (15), Equation (3) only considers the influence of gas deformation on the wave velocity and has a coefficient of correction. Equation (3) is basically the same as Equation (4). Therefore, Equation (3) is only applicable to the standard gas-fluid two-phase flow with α, which is neither close to 1 nor close to 0. Henry et al. believed that Equation (3) was applicable only when 0.01 ≤ α ≤ 0.5.
Compared with Equations (5) and (15), Equation (5) only neglects the influence of pipe wall deformation on the wave velocity and is more applicable than Equation (2). However, this equation is not suitable for the two-phase flow with a low gas content, small wall thickness, and small elasticity, especially for the hydraulic transition process of the common water distribution system.
Equation (9) has a wider application scope than Equations (3) and (5). Compared with Equation (15), Equation (9) only omits one factor, ρ g × α. For a two-phase flow of air and water, even if α = 0.9, which is impossible in a homogeneous two-phase flow of gas-water, the omitted ρ e α is only 1.2% of ρ e (1 − α). Therefore, Equation (9) has a certain calculation accuracy in the ordinary two-phase flow of gas and water and in the common hydraulic transient process.
Equations (10) and (13) are applicable to the hydraulic transient process when α is smaller, but Equation (10) is not applicable to the cavity of liquid column separation.
Equation (15) has been widely applied in the process of hydraulic transience since this equation comprehensively considers the effects of gas deformation, liquid deformation, and pipe wall deformation on the wave velocity.

Model of Collapsing Cavity
There are many causes for the formation of liquid column separation. It is generally believed that a collapsing cavity is composed of expanded bubbles released from water and vapors vaporized from liquid. The pressure of the first liquid column separation could be reduced to the vaporization pressure of the liquid, and then the continuous peak of each subsequent wave gradually decreases [32].
When the depressurization wave is not enough to reduce the pressure in the tube below the vaporization pressure of the liquid, the main cause for the increase in the liquid porosity is gas release and expansion. This increase in the liquid porosity is slight, and the collapsing phenomenon is not obvious. However, when the pipeline is flat and long enough, the boost wave cannot return within a short time and gas release and bubble expansion occur in a long pipeline. The continuous accumulation of bubbles forms a cavity, which hinders the fluid flow and greatly reduces the wave velocity of the water hammer. When the air content in the fluid is high, even if the pressure is not reduced below the vaporization pressure of water, the collapsing phenomenon occurs due to the gas release and bubble expansion under the state of low pressure. Most of the cavities formed in such a case are air. When the boost wave returns, in the closing process of the water columns, air absorption by the fluid is slow and air that is not dissolved in time acts as a cushion, thus reducing the pressurization of the water column. The lower the gas content of the fluid was, the closer the collapsing pressure was to the vaporization pressure of water [32]. Nevertheless, more bubbles in the pipeline also cause the risks of air mass interception and a pump-starting hammer caused by water impact.
When the depressurization wave decreases the pressure in the tube below the vaporization pressure of liquid and the critical pressure of bubbles (p* < p ≤ PV), the pressure range is small. p(v) = 0.024 kg/cm 2 when the water temperature is 20 • C. When the temperature increases, this range may be expanded. For example, when the water temperature is 50 • C, p(v) = 0.125 kg/cm 2 . The inhomogeneity of bubbles in fluids and the temperature instability in the transient process affect the solution of critical pressure p*. It seems impossible to solve for p* accurately. Usually, the value of p* can only be roughly estimated [33]. Normally, the gas accumulated in the pipeline can be eliminated after long-term steady-state operation, and the pipeline system is relatively tight. The gas content in the fluid is low, and the effects of released gas or bubble expansion are not obvious in the hydraulic transient process. The preliminary calculation can be neglected. The collapsing cavity is believed to be a steam chamber produced by liquid vaporization and the pressure of the collapsing cavity maintains the pressure of liquid vaporization [31].
Taking the pipe segment as the calculation unit and a water column as an elastic body, the water column is stretched rather than being vaporized immediately to produce a cavity when the flow rate of the latter node is slightly larger than that of the former node. To solve this problem, the statistical data of the water quantity in each node of the pipe segment are added. ∆t(Q a − Q c ) represents the variation in the water quantity in a pipe segment, where Q a indicates the flow rate of the former node; Q c indicates the flow rate of the latter node. The vaporization pressure of water is approximately the pressure generated by a -10 m high water column. The elastic modulus of water at the normal temperature is K = 2.18 × 10 9 Pa. Hence, the critical stretched water volume required for the formation of a cavity can be easily calculated according to the volume formula of the elastic modulus [34]: K = − dPV 0 dV , where V 0 indicates the volume of unstretched water; dV indicates the volume variation and can be expressed as the difference between the pipe segment volume and the calculated water volume; and dP indicates the vaporization pressure of water and is approximately the pressure generated by −10 m high water column, depending on the actual conditions. According to the calculated critical value, the cavity state can be determined. When the water volume decreases below the critical value, a cavity is formed. When the water volume is greater than the critical value, the elasticity of the fluid and pipe wall compensate for the volume reduction. Based on the consideration of the fluid elasticity, the shortcomings of the traditional calculation method can be avoided to obtain a more precise calculation result.
If the water volume in the pipe segment decreases below the critical value, a cavity is formed. The void fraction can be calculated with the flow difference between two successive pipe segments.

Method of the Floating Grid
Before collapsing occurs, the pressure is calculated with Equations (18) and (19) [35]: If collapsing occurs at a point, this means that a boundary condition of constant pressure is given at this point. If this point is at the headend, it means that the control conditions for the process of liquid column separation cannot work, at least before the water column is rebridged at this point. If this point is at a certain point in the pipeline, then the pipeline system is divided into two segments. Therefore, the wave of the water hammer in the respective water column is diffused, reflected, and superposed in a way that can be expressed with Equations (18) and (19). If the length of the broken water column is short and the cavitation area is short, it can be neglected. Moreover, if the depressurization time is not long and the gas release is not significant in the transient state, each segment after the collapsing can be calculated with Equations (18) and (19). Therefore, the grids of the characteristic equation are regular, as shown in Figure 1a.
Energies 2020, 13, x FOR PEER REVIEW 9 of 17 If collapsing occurs at a point, this means that a boundary condition of constant pressure is given at this point. If this point is at the headend, it means that the control conditions for the process of liquid column separation cannot work, at least before the water column is rebridged at this point. If this point is at a certain point in the pipeline, then the pipeline system is divided into two segments. Therefore, the wave of the water hammer in the respective water column is diffused, reflected, and superposed in a way that can be expressed with Equations (18) and (19). If the length of the broken water column is short and the cavitation area is short, it can be neglected. Moreover, if the depressurization time is not long and the gas release is not significant in the transient state, each segment after the collapsing can be calculated with Equations (18) and (19). Therefore, the grids of the characteristic equation are regular, as shown in Figure 1a. In Figure 1, △…△… indicates the spot and time range where the collapsing occurs. However, if the broken water column and the cavitation area are long and the gas release causes a significant change in the wave velocity, the calculation and the grids of the characteristic line are much more complicated [20,36].
When collapsing occurs in a long horizontal pipe segment, a long liquid column separation area often appears. The pressure in this area remains PV. In a previous study on the calculation of liquid column separation [17,31], the length of the breaking water column was neglected and treated as a fixed point, thus generating calculation errors in the flow velocity of the collapsing cavity node. Therefore, the flow velocity at the boundary point of the cavity should be used to indicate the flow change of the whole cavity. In the traditional fixed characteristic grid method, the flow rate at each fixed node is calculated. When the length of the collapsing cavity is too long to be neglected, the velocity calculated at the fixed node is significantly different from that calculated at the cavity boundary. With the extension of computing time, the cumulative errors in the repeated calculation at the node cannot be ignored. To reduce the errors, the following methods combining virtual In Figure 1, . . . . . . indicates the spot and time range where the collapsing occurs. However, if the broken water column and the cavitation area are long and the gas release causes a significant change in the wave velocity, the calculation and the grids of the characteristic line are much more complicated [20,36].
When collapsing occurs in a long horizontal pipe segment, a long liquid column separation area often appears. The pressure in this area remains P V . In a previous study on the calculation of liquid column separation [17,31], the length of the breaking water column was neglected and treated as a fixed point, thus generating calculation errors in the flow velocity of the collapsing cavity node. Therefore, the flow velocity at the boundary point of the cavity should be used to indicate the flow change of the whole cavity. In the traditional fixed characteristic grid method, the flow rate at each fixed node is calculated. When the length of the collapsing cavity is too long to be neglected, the velocity calculated at the fixed node is significantly different from that calculated at the cavity boundary. With the extension of computing time, the cumulative errors in the repeated calculation at the node cannot be ignored. To reduce the errors, the following methods combining virtual variables with linear interpolation should be adopted. For the sake of simplicity, the characteristic line method of fixed grids can be adopted at other points. This method is illustrated below (Figure 1b).
It is assumed that liquid column separation occurs at the headend of the pipeline at t = t1. The flow parameters at Point B are known. At t = t2, assuming that the boundary point of liquid column separation reaches R, then the pressure at Point R is P = Pv and the velocity V at Point R is unknown. The characteristic line from Point B does not intersect at Point R, so the flow velocity at Point R cannot be obtained directly. First, a virtual Point C with a velocity of Vc and a pressure of Pc = Pv is set. Vc can be obtained easily with the characteristic line of Point B: Then, V R can be obtained with the linear interpolation method. According to Figure 1b, we obtain: There are two variables, namely, x R and V R . To solve Equation (22), an additional equation is provided below: Substituting Equation (23) into Equation (22) gives: or In the calculation for the operation period of t2 = t3, the pressure and velocity of point P can be determined first. When V C = V R , the boundary value of each period can be calculated with the water hammer equations. Figure 2a shows the schematic diagram of the experimental device. The water transportation system was designed to transport water from the water tank on the first floor (simulated water source) to the water tank on the roof. A circulation loop was formed from overflow pipes on the top of the water tank on the roof and the water tank on the first floor. Other components in the pipeline system include valves for control and repair, pressure sensors, electromagnetic flow meters, and pressure gauges. A special quick-closing valve was set at the outlet of the pump for flow control. Four observation segments made of plexiglass tubes were set in the pipeline. The system was 400 m long and allowed the simultaneous occurrence of water hammers of cavities collapsing with multiple interruptions. Eight pressure measurement points were set in the pipeline system to measure the pressure in real time. Figure 2b shows the diagram of the experimental pipeline system. The pipe segments were connected with clamps. The main pipe was composed of DN100 galvanized steel pipe and organic glass pipe. Valves and float ball type fluid level gauges were installed in the tank. Manual gate valves and butterfly valves were installed in pipe segments to facilitate maintenance and control for stabilizing the flow rate. In this experiment [32], there were two cavities and data recording and phenomenon observation were realized. The two cavities occurred at the Number 1 and Number 2 pressure measuring points (Figure 1a). The elevations of the Number 1 and Number 2 pressure measuring points were 0.75 m and 7.21 m, respectively.

Instrument Calibration
An electromagnetic flowmeter was used in the experiment, and it had been calibrated before delivery. The range of the omega pressure sensor was −1 bar ~ 15 bar. The pressure sensor can be calibrated with the aid of a pressure gauge. During steady-state operation, the consistency between the displayed values of the electromagnetic flowmeter and pressure sensor and the stored values of In this experiment [32], there were two cavities and data recording and phenomenon observation were realized. The two cavities occurred at the Number 1 and Number 2 pressure measuring points (Figure 1a). The elevations of the Number 1 and Number 2 pressure measuring points were 0.75 m and 7.21 m, respectively.

Instrument Calibration
An electromagnetic flowmeter was used in the experiment, and it had been calibrated before delivery. The range of the omega pressure sensor was −1 bar~15 bar. The pressure sensor can be calibrated with the aid of a pressure gauge. During steady-state operation, the consistency between the displayed values of the electromagnetic flowmeter and pressure sensor and the stored values of LabVIEW was checked to ensure the transmission accuracy. The parameters in the experiment were accurately measured to ensure that the parameters of the simulation calculation were consistent with the parameters of the experimental platform.

Experimental Methods
After the water pump started and entered the normal operation state, the manual gate valve near the water pump outlet was adjusted to obtain a stable initial flow rate V 0 in the whole pipeline system. After stable operation for an hour, the accumulated gas in the pipeline was removed. After the initial state was determined, the valve closing speed was adjusted as required, and the switch of the speed governing box of the electric control valve was opened. The remote control valve was remotely closed to open the electric control valve. Finally, the electric control valve was closed by remote control to generate the water hammer of cavities collapsing with multiple interruptions. During the experimental period, experimental data and images were recorded in real time and stored through LabVIEW.

Results and Discussion
According to the experimental conditions, the hydraulic transition process at the Number 1 and Number 2 pressure measuring points was simulated. The initial flow rate was 77 m 3 /h, and the valve closing time was 0.9 s. Under the same working conditions, the operation was repeated three times and the difference among the measurement results should be within five ten thousandths. Only one group of data was selected to analyze the cavity.
The original model was used to simulate the experimental conditions, and the simulation results were compared with the measured pressure ( Figure 3). Taking the 2.8 s after pump stop as an example, the pressure calculated with the original model was 6.5 bar higher than the measured pressure, and the convergence rate of the original model was slow. At 12 s after pump stop, the pressure still oscillated obviously and still presented the state of a cut-off cavity, which was far from the actual conditions. Energies 2020, 13, x FOR PEER REVIEW 12 of 17 accurately measured to ensure that the parameters of the simulation calculation were consistent with the parameters of the experimental platform.

Experimental Methods
After the water pump started and entered the normal operation state, the manual gate valve near the water pump outlet was adjusted to obtain a stable initial flow rate V0 in the whole pipeline system. After stable operation for an hour, the accumulated gas in the pipeline was removed. After the initial state was determined, the valve closing speed was adjusted as required, and the switch of the speed governing box of the electric control valve was opened. The remote control valve was remotely closed to open the electric control valve. Finally, the electric control valve was closed by remote control to generate the water hammer of cavities collapsing with multiple interruptions. During the experimental period, experimental data and images were recorded in real time and stored through LabVIEW.

Results and Discussion
According to the experimental conditions, the hydraulic transition process at the Number 1 and Number 2 pressure measuring points was simulated. The initial flow rate was 77 m 3 /h, and the valve closing time was 0.9 s. Under the same working conditions, the operation was repeated three times and the difference among the measurement results should be within five ten thousandths. Only one group of data was selected to analyze the cavity.
The original model was used to simulate the experimental conditions, and the simulation results were compared with the measured pressure ( Figure 3). Taking the 2.8 s after pump stop as an example, the pressure calculated with the original model was 6.5 bar higher than the measured pressure, and the convergence rate of the original model was slow. At 12 s after pump stop, the pressure still oscillated obviously and still presented the state of a cut-off cavity, which was far from the actual conditions. The optimized model was used to simulate the experimental conditions and the simulation results were compared with the measured pressure. A comparison between the measured pressure and the simulated pressure at the Number 1 pressure measuring point is shown in Figure 4. We take 2.8 s after pump stop as an example. The simulated pressure was 0.8 bar higher than the measured pressure. Compared with the original model optimization, the improved model showed faster The optimized model was used to simulate the experimental conditions and the simulation results were compared with the measured pressure. A comparison between the measured pressure and the simulated pressure at the Number 1 pressure measuring point is shown in Figure 4. We take 2.8 s after pump stop as an example. The simulated pressure was 0.8 bar higher than the measured pressure. Compared with the original model optimization, the improved model showed faster convergence. After 6 s, no cavity existed, and the pressure gradually returned to a stable state. The pressure data of the Number 1 pressure measuring point in the process of hydraulic transition were analyzed by one-way ANOVA. The three groups of data were the simulated value before optimization (time1), the simulated value after optimization (time2), and the measured value of the experiment (time3).
The results of analysis of these three data groups ( Figure S1 in the Supplementary Information) showed that F(2, 33) = 15.950 (p < 0.001). The results also showed that there were significant differences among the preoptimization simulated value (time1), the simulated postoptimization value (time2), and the measured value (time3). The differences were then analyzed by two-to-two comparisons of the three sets of data. The results of the analysis ( Figure S2) showed a significant difference between the preoptimization simulated value (time1) and the measured value (time3) (p < 0.001). A comparison of the optimized simulated value (time2) with the measured value (time3) showed that p = 0.418 (p > 0.05); there was no significant difference between the optimized simulated value (time2) and the measured value (time3). Figure 5 shows a comparison between the measured pressure and the simulated calculated pressure obtained with the original model at the Number 2 pressure measuring point. We take 4.1 s after pump stop as an example to illustrate. The pressure calculated with the original model was 6.7 bar higher than the measured maximum pressure, and the convergence rate of the original model was slow. After 12 s of pump stop, the pressure still oscillated obviously and still presented the state of a cut-off cavity, which was far from the actual conditions. The pressure data of the Number 1 pressure measuring point in the process of hydraulic transition were analyzed by one-way ANOVA. The three groups of data were the simulated value before optimization (time1), the simulated value after optimization (time2), and the measured value of the experiment (time3).
The results of analysis of these three data groups ( Figure S1 in the Supplementary Information) showed that F(2, 33) = 15.950 (p < 0.001). The results also showed that there were significant differences among the preoptimization simulated value (time1), the simulated postoptimization value (time2), and the measured value (time3). The differences were then analyzed by two-to-two comparisons of the three sets of data. The results of the analysis ( Figure S2) showed a significant difference between the preoptimization simulated value (time1) and the measured value (time3) (p < 0.001). A comparison of the optimized simulated value (time2) with the measured value (time3) showed that p = 0.418 (p > 0.05); there was no significant difference between the optimized simulated value (time2) and the measured value (time3). Figure 5 shows a comparison between the measured pressure and the simulated calculated pressure obtained with the original model at the Number 2 pressure measuring point. We take 4.1 s after pump stop as an example to illustrate. The pressure calculated with the original model was 6.7 bar higher than the measured maximum pressure, and the convergence rate of the original model was slow. After 12 s of pump stop, the pressure still oscillated obviously and still presented the state of a cut-off cavity, which was far from the actual conditions.  Figure 6 shows a comparison between the measured pressure and the simulated calculated pressure obtained with the improved model at the Number 2 pressure measuring point. The maximum pressure calculated by the optimized model was close to the measured maximum pressure and the peak pressure occurred 4.1 s after the pump stopped. The simulated pressure was 0.3 bar higher than the measured pressure. Compared with the original model optimization, the improved model showed faster convergence. After 6 s, no cavity existed and the pressure gradually returned to a stable state. After 12 s, the pressure gradually recovered. The pressure data of the Number 2 pressure measuring point in the process of hydraulic transition were analyzed by one-way ANOVA. The three groups of data were the simulated value before optimization (time1), the simulated value after optimization (time2), and the measured value of the experiment (time3).
The results of analysis of these three data groups ( Figure S3) showed that F(2, 21) = 13.893, (p < 0.001). The results also showed that there were significant differences among the preoptimization simulated value (time1), the simulated postoptimization value (time2), and the measured value (time3). The differences were then analyzed by two-to-two comparisons of the three sets of data. The results of the analysis ( Figure S4) showed a significant difference between the preoptimization simulated value (time1) and the measured value (time3) (p < 0.001). A comparison of the optimized simulated value (time2) with the measured value (time3) showed that p = 0.934 (p > 0.05); there was  Figure 6 shows a comparison between the measured pressure and the simulated calculated pressure obtained with the improved model at the Number 2 pressure measuring point. The maximum pressure calculated by the optimized model was close to the measured maximum pressure and the peak pressure occurred 4.1 s after the pump stopped. The simulated pressure was 0.3 bar higher than the measured pressure. Compared with the original model optimization, the improved model showed faster convergence. After 6 s, no cavity existed and the pressure gradually returned to a stable state. After 12 s, the pressure gradually recovered.  Figure 6 shows a comparison between the measured pressure and the simulated calculated pressure obtained with the improved model at the Number 2 pressure measuring point. The maximum pressure calculated by the optimized model was close to the measured maximum pressure and the peak pressure occurred 4.1 s after the pump stopped. The simulated pressure was 0.3 bar higher than the measured pressure. Compared with the original model optimization, the improved model showed faster convergence. After 6 s, no cavity existed and the pressure gradually returned to a stable state. After 12 s, the pressure gradually recovered. The pressure data of the Number 2 pressure measuring point in the process of hydraulic transition were analyzed by one-way ANOVA. The three groups of data were the simulated value before optimization (time1), the simulated value after optimization (time2), and the measured value of the experiment (time3).
The results of analysis of these three data groups ( Figure S3) showed that F(2, 21) = 13.893, (p < 0.001). The results also showed that there were significant differences among the preoptimization simulated value (time1), the simulated postoptimization value (time2), and the measured value (time3). The differences were then analyzed by two-to-two comparisons of the three sets of data. The results of the analysis ( Figure S4) showed a significant difference between the preoptimization simulated value (time1) and the measured value (time3) (p < 0.001). A comparison of the optimized simulated value (time2) with the measured value (time3) showed that p = 0.934 (p > 0.05); there was The pressure data of the Number 2 pressure measuring point in the process of hydraulic transition were analyzed by one-way ANOVA. The three groups of data were the simulated value before optimization (time1), the simulated value after optimization (time2), and the measured value of the experiment (time3).
The results of analysis of these three data groups ( Figure S3) showed that F(2, 21) = 13.893, (p < 0.001). The results also showed that there were significant differences among the preoptimization simulated value (time1), the simulated postoptimization value (time2), and the measured value (time3). The differences were then analyzed by two-to-two comparisons of the three sets of data. The results of the analysis ( Figure S4) showed a significant difference between the preoptimization simulated value (time1) and the measured value (time3) (p < 0.001). A comparison of the optimized simulated value (time2) with the measured value (time3) showed that p = 0.934 (p > 0.05); there was no significant difference between the optimized simulated value (time2) and the measured value (time3).
The calculated results of the original model were significantly higher than the measured values and did not converge. The results of the optimized model were closer to the actual measurement values. There were still some gaps between the calculated results of the optimized model and the measured values (Figures 4 and 6). The simulation calculation values were higher than the measured values and the convergence rate was slow. The main reasons are provided as follows. (1) The vibration of the laboratory pipeline reduced the water hammer wave. The pipes in the laboratory were supported and fixed by angle steel frames. In the process of the water hammer of the cavity collapsing with multiple interruptions, there was a large vibration. Part of the energy of the water hammer wave was converted into vibration energy [17], thus causing the energy reduction of water hammer wave. (2) In the pipeline system, there were more turning elbows, which increased the head loss. (3) The sensitivity of the sensor was not high enough, and the small pressure fluctuation in the actual measurement was not recorded. Therefore, the measured pressure line converged faster than the simulated values.

Conclusions
Aiming at the shortcomings of the traditional DVCM model, research was carried out and results were achieved in this paper. The new water hammer wave velocity formula considers the influence of gas holdup, liquid deformation, and pipe wall deformation. The water hammer wave velocity changes with the real-time change in gas holdup. In the new cavity model, taking the pipe section as the calculation unit, when judging whether there is a cutoff cavity, the elasticity of the liquid is considered. When the amount of water reduction in the pipe section exceeds the elasticity of the water body and the compensation ability of the pipe wall, the cutoff cavity is considered. This can effectively reduce the occurrence of pseudo water hammers in simulation calculations. The floating grid method takes the influence of cavity length into account, which can bring the calculation results closer to the actual situation. It can be seen from the test, compared with the traditional DVCM model before optimization, the calculated results of the optimized model were closer to the experimental values. Additionally, the simulated value of the new model was higher than the measured value. One of the reasons for this result is the vibration of the experimental pipeline and the connection of the experimental pipeline by multiple elbows, which also leads to faster pressure wave attenuation in the measured value than in the simulated value. The second reason may be that the sensitivity of the sensor was not high enough to record the small pressure fluctuation in the actual measurement.
Supplementary Materials: The following are available online at http://www.mdpi.com/1996-1073/13/5/1103/s1, Figure S1: The comparison of the pre-optimization simulated value, post-optimization simulated value, and the measured value (No. 1 manometer point), Figure S2: The two-to-two comparisons of the pre-optimization simulated value, optimized simulated value and the measured value (No. 1 manometer point), Figure S3: The comparison of the pre-optimization simulated value, post-optimization simulated value, and the measured value (No. 2 manometer point), Figure S4: The two-to-two comparisons of the pre-optimization simulated value, optimized simulated value and the measured value (No. 2 manometer point).  pipe wall thickness (m) E pipe wall elastic Young's modulus K g , K e elastic moduli of gas phase and liquid phase M mass content of gas per unit volume P pressure t time coordinate u nonuniformity coefficient of pipe wall material V movement velocity of fluid in the pipe X position coordinate Greek letters α gas content in unit volume liquid γ gaseous polytropic index ρ e density of liquid phase ρ g density of gas phase ρ m gas-liquid mixed fluid density R gas constant