#### Abstract

Accurate simulation of cavitating flows in pipeline systems is important for cost-effective surge protection. However, this is still a challenge due to the complex nature of the problem. This paper presents a numerical model that combines the discrete vapor cavity model (DVCM) with the quasi-two-dimensional (quasi-2D) friction model to simulate transient cavitating flows in pipeline systems. The proposed model is solved by the method of characteristics (MOC), and the performance is investigated through a numerical case study formulated based on a laboratory pipeline reported in the literature. The results obtained by the proposed model are compared with those calculated by the classic one-dimensional (1D) friction model with the DVCM and the corresponding experimental results provided by the literature, respectively. The comparison shows that the pressure peak, waveform, and phase of pressure pulsations predicted by the proposed model are closer to the experimental results than those obtained by the classic 1D model. This demonstrates that the proposed model that combines the quasi-2D friction model with the DVCM has provided a solution to more accurately simulate transient cavitating flows in pipeline systems.

#### 1. Introduction

Significant hydraulic transients caused by pump trips or rapid valve closure in a pipeline system can induce pipe failure [1, 2]. These damaging transient events typically involve not only very high positive pressure but also very low negative pressure. Sometimes, the negative pressure can produce more severe damage relative to the positive pressure. In particular, when the pressure instantaneously drops to the saturated vapor pressure at a given liquid temperature, column separation occurs in transient pipe flows. This causes the appearance of vapor cavities to separate liquid columns at some specific locations such as closed ends and high points, and this is known as column separation. When a vapor cavity collapses, it can generate significant positive pressure peaks and the magnitude sometimes is higher than the Joukowsky pressure [3]. These activities could occur repeatedly in a transient event and seriously damage hydraulic pipeline systems.

Hydraulic transients with column separation have been explored in some aspects by many researchers [3]. A number of vaporous cavitation models have been developed for simulating transient cavitating flows in pipeline systems [4–6]. Among these models, the discrete vapor cavity model (DVCM) combined with the one-dimensional (1D) quasi-steady friction model is the classic modelling approach and is widely used [7]. It is relatively easy to implement and is solved by the method of characteristics (MOC) [8]. Although the results provide insights into the real phenomenon, the classic 1D model may generate numerical oscillations and unrealistic pressure spikes associated with multicavity collapsing [4].

In order to improve these issues, it is necessary to calculate the frictional effect more accurately. Some previous studies have proposed improved friction models such as 1D unsteady friction models and 1D quasi-steady friction models with correction factor. Those unsteady friction models include the constant coefficient instantaneous acceleration-based models [9] and the convolution-based models [10, 11]. Results show that pressure pulsations of cavitating transient flows simulated by the DVCM with 1D unsteady friction models have a good agreement with those of the experimental ones. Mosharaf-Dehkordi and Firoozabadi [7] investigated the effect of quasi-steady friction term on the accuracy of the classic DVCM model and found that the friction term affects the calculation of the timing of pressure pulses.

Generally, the 1D model is relatively easy to program and has a high computational efficiency, but it cannot accurately model the frictional effect, especially for complex flow conditions such as cavitating flow. Taking into account the velocity profile, the two- or three-dimensional friction models are able to simulate transient flows in a pipe more accurately than the 1D models [12]. Some studies [13, 14] make use of computational fluid dynamic (CFD) methods for simulation of transient cavitating flows in pipeline systems. However, the CFD methods require much more computational costs than the 1D method. Recently, quasi-two-dimensional (quasi-2D) models have been introduced to simulate transient cavitating pipe flows. Pezzinga and Cannizzaro [15] developed a transient cavitating 2D model combining a distributed vaporous cavitation model with the quasi-2D model. The predictor-corrector method is used to solve the transient cavitating 2D model. Using the 2D model proposed by Pezzinga and Cannizzaro, Gao et al. [16] compared the influence of using two different turbulence models, including the two-region turbulence model and the five-region turbulence model, in the quasi-2D model. Their results showed that numerical results simulated using both turbulence models are similar. Santoro et al. [17] compared the performance of various forms of the 1D model and analyzed the difference of features between the 1D model and 2D model such as the selection of grid size and weighting parameters. Their results showed both 1D and 2D models are weakly sensitive to grid size and the 2D model is sensitive to the selection of weighting parameters.

The early studies have demonstrated that the quasi-2D friction model, combined with the vaporous cavitating model, is capable of simulating transient cavitating flows more accurately. This paper proposes a numerical water hammer model coupling the quasi-2D friction model with the DVCM. For the quasi-2D model, Vardy and Hwang model [18] is employed, which is accurate and stable. And Jang et al. [19] improved an efficient calculation method of the Vardy and Hwang model and coupled the quasi-2D friction model with the 1D model. Based on the Jang et al. method [19], the proposed model is solved by the MOC with staggered grid for transient cavitating flows in pipeline systems. Numerical validation is conducted on a pipeline system that is established based on a laboratory system reported in the literature [4]. The results obtained from the proposed approach are compared with those from the 1D quasi-steady friction model with the DVCM and with the experimental results provided by the literature [4], respectively.

#### 2. Mathematical Model

##### 2.1. Governing Equations

When the pressure is above the vapor pressure, the two-dimensional governing equations for transient pipe flows with the assumption of the axially symmetric flows can be expressed as [18]where *x* is the distance in axial direction of the pipe, *r* is the distance in radial direction from the pipe axis, *t* is the time, *ρ* is the liquid density, *H* is the pressure head, *u* and are the axial and radial velocities, respectively, a is the wave speed, is the gravitational acceleration, and *τ* is the shear stress. Through integrating equations (1) and (2) across the pipe cross section, the 1D water hammer equations are derived as [20]withwhere *Q* is the discharge, *A* is the cross-sectional area of the pipe, *R* and *D* are the radius and the diameter of the pipe, respectively, *ν*_{T} is the total viscosity, *π* is the circular constant, and is the wall shear stress.

##### 2.2. Discrete Vapor Cavity Model (DVCM)

The DVCM is the most commonly used model for simulation of transient cavitating pipe flows at the present time. It assumes that vapor cavities are allowed to form at any of the computational grid points if the pressure is computed to be below the vapor pressure [8]. The vapor cavities are thus confined to computational grid points and a pure liquid column with a constant pressure wave speed is assumed in between the adjacent grid points. The pressure head in a cavity is set equal to the vapor pressure head:where is relative vapor pressure head at the reference temperature.

A continuity equation for the volume of a discrete vapor cavity is given by [8]where *Q* and are the discharges at the downstream and upstream sides of the cavity, respectively. The numerical integration of equation (6) using the staggered MOC grid can be expressed aswhere *ψ* is the weighting coefficient (0.5 ≤ *ψ* ≤ 1), *t* is the computational time, is the time interval, and the subscript is the upstream sides at a cavity.

#### 3. Numerical Scheme

The solution for transient flows in the homogeneous liquid zone is calculated by a quasi-2D model using the method of characteristics (MOC). The characteristic forms of 2D transient flow equations are given aswhere the shear stress . The detailed derivation process of equation (8) from equations (1) and (2) is given in Appendix. And the corresponding 1D transient flow equations are expressed aswhere *B* = *a/gA*. The computational grids for the discretization of characteristic equations are shown in Figure 1. As shown in Figure 1(a), the pipe is discretized into *N*_{r} cylinders with varying thickness in the radial direction. The pipe length *L* is divided into *N*_{x} reaches of constant length Δ*x* in the axial direction. The axial velocity *u* is located in the middle of each reach in the radial direction, and both radial velocity and turbulent viscosity are located at the boundaries of each reach in the radial direction. The time internal Δ*t* is computed by the equation Δ*t* = Δ*x/a.*

**(a)**

**(b)**

Integrating equation (8) on the characteristic lines between time *n*Δ*t* and (*n* + 1)Δ*t* (Figure 1(b)), the discretized forms of equations (8) can be expressed aswithwhere *ε* and *θ* are the weighting coefficients of the viscous term and radial flux term in the characteristic equations, respectively (0.5 ≤ *ε* ≤ 1, 0.5 ≤ *θ* ≤ 1), *q* (=*rv*) is the radial flux, the subscripts *i* and *j* indicate the axial and radial step indexes, respectively, the superscript *n* denotes the time step index, and *r*_{j} and *r*_{cj} are the coordinates of the boundary and middle points of reaches in the radial direction, respectively. The source terms *K*_{pi,j} and *K*_{ni,j} are the known values, whose elements depend on *H*, *u*, and *q* at the previous time level.

Through adding and subtracting algebraic operation in equations (9)–(11), the equations for axial velocity, pressure head, and radial velocity can be expressed as [21]

Both equations are the fully implicit equations in terms of the unknowns (*u*, *H*, and *q*), and these solutions are calculated using the Thomas algorithm [22].

The discretized forms of 1D characteristic equation (10) are expressed as follows:with

When the pressure is below the saturated vapor pressure, the vapor cavity is formed at the fixed grid points. The influence of radial pipe displacements on the cavity volume is neglected. Meanwhile, the order of the radial velocity is very small relative to the axial velocity. Thus, the radial velocity and its derivatives can be neglected between the liquid column and the vapor cavity. The assumption of neglecting the radical velocity is used by solving the 2D transient cavitating flow equations. When the vapor cavity is formed, the characteristic forms of 2D transient cavitating flow equations are expressed aswith

The pressure head at the cavity is equal to the saturated vapor pressure, shown as follows:

Hence, the upstream and downstream axial velocities can be solved by using equations (18) and (19) using the Thomas algorithm.

Then, the wall shear stress is computed as follows:

Finally, the upstream and downstream discharges are computed using the discretized forms of 1D characteristic equation (9) in staggered grid, as shown in Figure 2:with

The computational procedure for the proposed quasi-2D model can be summarized as follows:(1)Calculate initial steady-state conditions including discharge, pressure head, and velocity profiles(2)Calculate the transient-state conditions in the next time step(3)Calculate the axial velocity profiles by using equation (13)(4)Calculate the discharge and pressure head by using equations (15) and (16)(5)If the pressure is below the vapor pressure, the pressure head is equal to the relative vapor pressure head; if the pressure is above the vapor pressure, go to Point 10(6)Calculate the upstream and downstream axial velocity profiles by using equations (18) and (19)(7)Calculate the wall shear stress by using equation (20)(8)Calculate the upstream and downstream discharges by using equations (21) and (22)(9)Calculate the volume of the cavity by equation (8); if the volume is smaller than zero, it is set to zero(10)Print results(11)Go to Point 2 until the maximum computational time

#### 4. Numerical Results and Discussion

In order to validate the accuracy of the proposed model, the numerical results computed by both the classic 1D model and the proposed quasi-2D model are compared with experimental results obtained by Bergant and Simpson [4]. The experimental system was a typical reservoir-pipe-valve system, which was composed of a copper straight and sloping pipe (length 37.23 m, internal diameter 22.1 mm, and pipe slope 3.2°). The transient flows were caused by a rapid downstream valve closure. The wave speed was 1319 m/s. Experimental results on the pressure oscillation at the valve for the initial mean velocities of *V*_{0} = 0.30 m/s, 0.71 m/s, and 1.40 m/s are used in the current research. Through testing different grid sizes, the *N*_{x} value of 16 and the *N*_{r} value of 90 are used in the proposed model, which are capable of simulating transient cavitating pipe flows effectively and accurately.

##### 4.1. Effects of the Weighting Coefficients *ε*, *θ*, and *ψ* in the Proposed Quasi-2D Model

The pressure heads at the downstream valve of pipe predicted by the proposed quasi-2D model with different weighting coefficients *ε* for three cases of the initial mean velocities are presented in Figure 3. The weighting coefficients *θ* and *ψ* values are specified as 1.00 and 0.50 for three cases of the initial mean velocities in Figures 3(a)–3(c), respectively. It is clear that the pressure heads computed by the proposed quasi-2D model with different weighting coefficients *ε* for the cases of *V*_{0} = 0.30 m/s and 1.40 m/s have very slight differences in Figures 3(a) and 3(c). In Figure 3(b), the first pressure peak and the first and second cavity duration time calculated using different weighting coefficients *ε* for the cases of *V*_{0} = 0.71 m/s are almost the same and show good agreement with the experimental results. But the numerical results after the first three pressure pulses are quite different. And the numerical results calculated using the weighting coefficient *ε* = 1.00 are closer to the experimental results, as shown in Figure 3(b).

**(a)**

**(b)**

**(c)**

The pressure heads at the downstream valve of pipe predicted by the proposed quasi-2D model with different weighting coefficients *θ* for three cases of the initial mean velocities are presented in Figure 4. The weighting coefficients *ε* and *ψ* values are both set to 1.00 for three cases of the initial mean velocities in Figures 4(a)–4(c). As observed in Figure 4(a), the magnitude of pressure heads computed using different weighting coefficients *θ* for the case of *V*_{0} = 0.30 m/s has some slight differences. But the pressure pulses after the short duration pressure peak obtained using the weighting coefficient *θ* = 0.50 occur in advance, comparing with the pressure pulses obtained using the weighting coefficients *θ* = 1.00 and 0.75. And the similar phenomenon also appears in Figures 4(b) and 4(c). It can be seen that the pressure pulses obtained using the weighting coefficient *θ* = 1.00 are closer to the experimental results.

**(a)**

**(b)**

**(c)**

The pressure heads at the downstream valve of pipe predicted by the proposed quasi-2D model with different weighting coefficients *ψ* for three cases of the initial mean velocities are presented in Figure 5. The weighting coefficients *ε* and *θ* values are both set to 1.00 for three cases of the initial mean velocities in Figures 5(a)–5(c). It can be seen that the timing of pressure pluses obtained using different weighting coefficients *ψ* is quite different for the cases of *V*_{0} = 0.30 and 0.71 m/s. It is clear that the magnitude and timing of pressure pluses computed using different weighting coefficients *ψ* for the case of *V*_{0} = 1.40 m/s have some slight difference. In Figure 5(a), the pressure pluses calculated using the weighting coefficient *ψ* = 0.50 are closer to the experimental results. But, in Figure 5(b), the numerical results obtained using the weighting coefficients *ψ* = 1.00 and 0.50 are more accurate than those obtained using the weighting coefficient *ψ* = 0.75.

**(a)**

**(b)**

**(c)**

Through the abovementioned analysis, the weighting coefficients *ε* and *θ* values of 1 are more suitable for simulating transient cavitating pipe flows in the proposed quasi-2D model. Then, the trial-and-error processes are required to confirm the weighting coefficient *ψ.* Therefore, the weighting coefficients *ε* and *θ* values are set to 1 in this study. Through trial-and-error processes, for the cases of *V*_{0} = 0.30 m/s and 0.71 m/s, the weighting coefficient *ψ* value is 0.50 in the proposed quasi-2D model. But, for the case of *V*_{0} = 1.40 m/s, the weighting coefficient *ψ* value is set to 1.00 in the proposed quasi-2D model. And the weighting coefficient *ψ* is set to 1.00 in the classic 1D model for the three cases of the initial mean velocities.

##### 4.2. Comparison of Numerical Results Obtained by the 1D Model and Quasi-2D Model

Comparison of the numerical and experimental results with *V*_{0} = 0.30 m/s is shown in Figure 6. The results indicate that both the classic 1D model and the quasi-2D model can predict the short-duration pressure peak, *H*_{sd}, after the collapse of the first cavity at the valve, and both models can calculate the duration of the first cavity accurately. *H*_{sd} is 95.60 m in the experiment, 98.94 m calculated by the quasi-2D model, and 100.22 m by the 1D model. In addition, the pressure head traces after that high pressure peak computed by the classical 1D model have some unrealistic pressure spikes. It is noted that the pressure head traces obtained by the proposed quasi-2D model better match the experimental ones than the results from the 1D model. The waveform and phase of the pressure pulses predicted by the quasi-2D model are more consistent with those of the experimental results.

In Figure 7, the pressure responses at valve calculated by the classic 1D model and the proposed quasi-2D model are compared with the experimental results for *V*_{0} = 0.71 m/s. Both the classic 1D model and the quasi-2D model can predict the first cavity duration after the Joukowsky pressure rise at the downstream valve accurately. But it can be observed that the pressure simulated by the classic 1D model has some unrealistic pressure spikes. Moreover, the unrealistic pressure head peaks obtained by the classic 1D model are much higher than the Joukowsky pressure rise. On the contrary, the waveform and phase of the pressure oscillations obtained by the quasi-2D model are in better agreement with those of the experimental ones. Meanwhile, the time of vapor cavities formation and collapse calculated by the quasi-2D model is closer to the experimental ones than the one obtained by the classic 1D model.

The pressure heads calculated by both the classic 1D model and the proposed quasi-2D model in comparison to the experimental results for *V*_{0} = 1.40 m/s are shown in Figure 8. It can be seen that both the classic 1D model and quasi-2D model cannot accurately predict the pressure head peaks after the collapse of the cavities. However, the quasi-2D model is able to reproduce the waveform and attenuation of pressure fluctuations closer to the experimental results. In addition, it is noted that the duration of the cavities and the time of pressure rise after the collapse of the cavities computed by the quasi-2D model better match the experimental ones, relative to the ones calculated by the classic 1D model.

Figure 9 shows the comparison of wall shear stress at midpoint of the pipe for *V*_{0} = 1.40 m/s computed by the classic 1D model and the proposed quasi-2D model. It can be seen that the behaviour of the wall shear stress calculated by the quasi-2D model in comparison to those obtained by the classic 1D model is quite different. The former shows a shorter and more intense pulsation process, whereas the latter changes gradually as the time increases. This is because the classic 1D model uses the average velocity in the cross-section of the pipe to calculate the friction dissipation term. However, the quasi-2D model calculates the wall shear stress using the instantaneous velocity profile along the axis of the pipe. Therefore, the quasi-2D model computed the dissipation caused by friction more accurately than the classic 1D model.

Figure 10 shows the behaviour of velocity profiles at the midpoint of the pipe calculated by the quasi-2D model. It can be observed that different behaviors of velocity profiles along the axis of the pipe are presented at different phases. When the flow is at the steady state, all components of the axial velocity are positive and there is a uniform variation in the centre of the pipe and a large velocity gradient near the pipe wall. When the pressure wave passes through the midpoint of the pipe at time *t* = *L*/*a*, the reverse flow appears near the pipe wall. At time *t* = 2*L*/*a*, all components of the axial velocity are negative and there exists a quite large velocity gradient near the pipe wall. However, the variation is gradual in the centre of the pipe. As the time increases, the velocity gradient decreases near the pipe wall. Therefore, the quasi-2D model, taking into account the velocity profile in the cross-section of the pipe, calculates the dissipation caused by friction more accurately than the classic 1D model.

The variations of vapor cavity volumes at valve for *V*_{0} = 1.40 m/s computed by the classic 1D model and the proposed quasi-2D model are shown in Figure 11. It can be seen that vapor volumes computed by the quasi-2D model are relatively smaller than the ones obtained by the classic 1D model. The size of vapor cavities calculated by both models decreases as the time increases. However, the rate of decrease as predicted by the quasi-2D model is faster than the one obtained by the classic 1D model. That is because different calculation results of discharges at the cavity are obtained by both models using the difference physical descriptions for the wall shear stress.

The above analysis shows that both the classic 1D model and the proposed quasi-2D model can predict the Joukowsky pressure rise accurately when the valve closes rapidly, and both models cannot predict the pressure peak after the collapse of the first cavity at the valve exactly. However, the result calculated by the quasi-2D model is better than the one obtained by the classic 1D model, and the quasi-2D model can reduce the unrealistic numerical oscillation of pressure with respect to the classic 1D model. The waveform and phase of the pressure fluctuations predicted by the quasi-2D model are closer to the experimental results with respect to those obtained by the classic 1D model. Meanwhile, the duration of vapor cavity formation and collapse calculated by the quasi-2D model is in better agreement with the experimental result than the one computed by the classic 1D model.

However, computational results obtained by the proposed quasi-2D model are sensitive to the values of the weighting coefficients *ε* and *θ* in the characteristic equations, and the value of the weighting coefficient *ψ* in DVCM. The selection of the value of the weighting coefficients *ε* and *θ* affects the phase of the pressure fluctuations. And the influence of the choice of the weighting coefficient *ψ* value is the waveform of the pressure fluctuations. Therefore, when the proposed quasi-2D model is used, the *ε* and *θ* values are tested firstly based on the experimental results. And, then, the *ψ* value is confirmed in order that the simulation results are closer to the experimental ones.

#### 5. Conclusions

This paper presents a quasi-2D approach for simulating transient cavitating flows in pipeline systems. The research combines the DCVM with the quasi-2D friction model and solves the resulting equations by MOC. The proposed model improves the accuracy of friction dissipation calculation in the classic DVCM model by using the instantaneous velocity profile in the cross-section of the pipe.

To validate the new approach, the numerical results calculated by the classic 1D model and the proposed quasi-2D model are compared with the corresponding experimental results. The comparison shows that the pressure peaks, the waveform, and the phase of the pressure fluctuations predicted by the proposed model are closer to the experimental results with respect to those obtained by the classic 1D model. Meanwhile, the duration of vapor cavity formation and collapse calculated by the proposed model is in better agreement with the experimental result than the one computed by the classic 1D model. The proposed quasi-2D approach provides a solution for more accurate simulation of predicting transient cavitating flows in pipeline systems.

#### Appendix

#### Derivation Process of Equation (8) from Equations (1) and (2)

The two-dimensional governing equations for transient pipe flows are identified as *L*_{1} and *L*_{2} from equations (1) and (2):

These equations are combined linearly using an unknown multiplier *λ*:

By rearranging equation (A2), the equation may be expressed as

If the unknown multiplier *λ* is written asthen equation (A3) becomes the ordinary differential equation as follows:

The solution of equation (A5) yields the two particular values of *λ*:

Substituting these values of *λ* into equation (A5), the characteristic equations could be expressed as

#### Nomenclature

A: | Cross-sectional area of the pipe (m^{2}) |

a: | Wave speed (m/s) |

D: | Pipe diameter (m) |

: | Acceleration of gravity (m/s^{2}) |

H: | Pressure head (m) |

: | Relative vapor pressure head at reference temperature (m) |

K_{n}: | Known values |

K_{p}: | Known values |

L: | Pipe length (m) |

n: | Time step index |

Q: | Discharge (m^{3}/s) |

q: | Radial flux (m^{2}/s) |

R: | Radius of the pipe (m) |

r: | Distance in radial direction from the pipe axis (m) |

r_{j}: | Coordinates of boundary points of reaches in radial direction (m) |

r_{cj}: | Coordinates of middle points of reaches in radial direction (m) |

t: | Time (s) |

: | Time interval (s) |

u: | Longitudinal velocity components (m/s) |

V_{cav}: | Cavity volume (m^{3}) |

: | Radial velocity components (m/s) |

x: | Distance in axial direction of pipe (m) |

*Greek symbols*

ε: | Weighting coefficient of the viscous term |

θ: | Weighting coefficient of the radial flux term |

: | Total viscosity (m/s^{2}) |

π: | Circular constant |

ρ: | Liquid density (kg/m^{3}) |

: | Wall shear stress (Pa) |

τ: | Shear stress (Pa) |

ψ: | Weighting coefficient. |

*Superscripts*

n: | Time step index. |

*Subscripts*

: | Upstream sides at a cavity |

i: | Axial step index |

j: | Radial step index |

*Abbreviations*

CFD: | Computational fluid dynamic |

1D: | One-dimensional |

DVCM: | Discrete vapor cavity model |

MOC: | Method of characteristics |

Quasi-2D: | Quasi-two-dimensional. |

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was financially supported by the Fundamental Research Funds for the Central Universities (Grant no. 2572018BJ14), National Natural Science Foundation of China (Grant nos. 51808102 and 51978202), Science Foundation of Harbin University of Commerce (18XN068), and the Natural Science Fund in Heilongjiang Province (LH2019E111).