Coupling Quasi-Two-Dimensional Friction Model and Discrete Vapor Cavity Model for Simulation of Transient Cavitating Flows in Pipeline Systems

School of Civil Engineering, Institute of Artificial Environment Control and Energy Application, Northeast Forestry University, Harbin, China School of Architecture, Harbin Institute of Technology, Harbin, China Laboratory of Cold Region Urban and Rural Human Settlement Environment Science and Technology, Ministry of Industry and Information Technology, Harbin, China Harbin University of Commerce, Harbin 150028, China Department of Mechanics, Kim Il Sung University, Pyong Yang 999093, Democratic People’s Republic of Korea School of Municipal and Environmental Engineering, Harbin Institute of Technology, Harbin 150090, China


Introduction
Significant hydraulic transients caused by pump trips or rapid valve closure in a pipeline system can induce pipe failure [1,2]. ese 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. is 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]. ese 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][5][6]. Among these models, the discrete vapor cavity model (DVCM) combined with the one-dimensional (1D) quasisteady 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. ose 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 quasisteady 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, quasitwo-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. e 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. eir 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. eir 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. e 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.
is 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]. e 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.

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] g a 2 zH zt 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 v are the axial and radial velocities, respectively, a is the wave speed, g is the gravitational acceleration, and τ is the shear stress. rough integrating equations (1) and (2) across the pipe cross section, the 1D water hammer equations are derived as [20] gA with where 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 τ w is the wall shear stress.

Discrete Vapor Cavity Model (DVCM).
e 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]. e 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. e pressure head in a cavity is set equal to the vapor pressure head: where H v is relative vapor pressure head at the reference temperature. A continuity equation for the volume of a discrete vapor cavity is given by [8] 2 Mathematical Problems in Engineering where Q and Q u * are the discharges at the downstream and upstream sides of the cavity, respectively. e numerical integration of equation (6) using the staggered MOC grid can be expressed as where ψ is the weighting coefficient (0.5 ≤ ψ ≤ 1), t is the computational time, Δt is the time interval, and the subscript u * is the upstream sides at a cavity.

Numerical Scheme
e solution for transient flows in the homogeneous liquid zone is calculated by a quasi-2D model using the method of characteristics (MOC). e characteristic forms of 2D transient flow equations are given as where the shear stress τ � − ρ] T zu/zr. e detailed derivation process of equation (8) from equations (1) and (2) is given in Appendix. And the corresponding 1D transient flow equations are expressed as where B � a/gA. e 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. e pipe length L is divided into N x reaches of constant length Δx in the axial direction. e axial velocity u is located in the middle of each reach in the radial direction, and both radial velocity v and turbulent viscosity v T are located at the boundaries of each reach in the radial direction. e time internal Δt is computed by the equation Δt � Δx/a.
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 as with Mathematical Problems in Engineering where ε 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. e 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. rough adding and subtracting algebraic operation in equations (9)-(11), the equations for axial velocity, pressure head, and radial velocity can be expressed as [21] εC u1j u n+1 Both equations are the fully implicit equations in terms of the unknowns (u, H, and q), and these solutions are calculated using the omas algorithm [22]. e 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. e 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.
us, the radial velocity and its derivatives can be neglected between the liquid column and the vapor cavity.
e 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 as with e 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 omas algorithm. en, 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 e 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)  e 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. e 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 e 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. e 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.   Figure 2: MOC method with staggered grid [8].
weighting coefficients ψ for three cases of the initial mean velocities are presented in Figure 5. e 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. rough 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. en, the trial-and-error processes are required to confirm the weighting coefficient ψ. erefore, the weighting coefficients ε and θ values are set to 1 in this study.
rough trial-and-error processes, for the cases of

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. e 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. e 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  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. e former shows a shorter and more intense pulsation process, whereas the latter changes gradually as the time increases. is is because the classic 1D model uses the average velocity in the crosssection 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. erefore, 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 � 2L/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. erefore, the quasi-2D model, taking into account the velocity profile in the cross-section of the pipe, calculates the dissipation   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. e 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. at 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. e 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. e 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. e 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. erefore, 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.

Conclusions
is paper presents a quasi-2D approach for simulating transient cavitating flows in pipeline systems. e research combines the DCVM with the quasi-2D friction model and solves the resulting equations by MOC. e 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. e 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. e proposed quasi-2D approach provides a solution for more accurate simulation of predicting transient cavitating flows in pipeline systems.