# Prediction Method of Cavitation Jet Wave Attenuation Based on Five-Equation Two-Fluid Model

Authors: You Fu; Zhaopeng Xie; Weiguo Zhao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, the seven-equation two-fluid flow model is streamlined to a five-equation numerical calculation method. This method is applied to predict the wave attenuation of the cavitation jet. Compared with 9 different experimental schemes in two separate laboratories, it is found that the velocity flow characteristics are different from normal transient flow when the cavitation jet is formed in the pipeline. Such a difference in velocity flow characteristics will cause changes in the prediction deviation of pressure response time and amplitude. A numerical method for predicting the attenuation of pressure fluctuation with the cavitation jet is presented taken into account these flow characteristics. In this method, the unsteady friction model is opened after the jet wave transmits 3–5 cycles, and the cavitation jet wave attenuation is well simulated. This calculation method could provide a research basis for pipeline leakage, vibration, noise, and other fields that need accurate pressure signals.


---

## Body

## 1. Introduction

Water hammer with column separation in a pipeline system can cause cavitation jet [1–3]. The formation of the cavitation jet will cause a sharp rise in pressure, which has a serious impact on the safety of the pipeline system [4–8]. In recent years, scholars have begun to use the two-fluid model to study the cavitation jet phenomenon both in time-effectiveness and calculation accuracy in the pipeline system [9–14]. A typical two-fluid two-phase flow model under one-dimensional direction requires the solution of 7–9 equations [15–23]. At the same time, the cavitation jet which is a typical transient flow phenomenon usually calculated in an explicit format. The complexity of the calculation model and solution format takes more computation time. Therefore, this method sometimes is not appropriate for engineering applications.There are two main factors influencing the calculation accuracy of water hammer cavitation jet in the pipeline system: (1) residual void of gas and (2) wave attenuation. Residual void of gas means that the volume of initial gas increases with the progress of time iteration. The principle of this phenomenon is that the initial calculation needs a certain amount of insoluble gas to satisfy the appeared second phase (water vapor) in the two-fluid model. The energy loss is not considered in the interaction between the two phases by the two-fluid model in theory. The expansion value of bubbles is often larger than the compression, which is expected to result in the void of gas gradually increasing with the advance of time. So the expansion and compression stages of the bubbles are unbalanced under the condition of constant relaxation coefficient in real calculation [19–21, 23], resulting in the residual void of gas. The residual void of gas will cause an obvious influence on the velocity of wave transmission. It will continuously affect the wave, which will lead to obvious errors in the prediction of water hammer wave.Wave attenuation refers to the attenuation of the pressure wave and the degree of energy dissipation with the advance of time, which is also an important part of transient research [24–32]. Accurate prediction of water hammer cavitation jet in a pipeline system is the basis of prediction of pipeline obstruction and leakage. The main problem of wave attenuation lies in the combination of unsteady friction model and cavitation jet calculation method. The two decay coefficient models added into the two-fluid model will take such solution process more complicated. When the large intensity of the cavitation jet is generated, the applicability of the unsteady friction model with the two-fluid model has been rarely reported in the literature. Therefore, it is an urgent problem to predict the wave attenuation of the cavitation jet.Another factor is the timeliness of the calculations. The calculation of the cavitation jet in the pipeline system is generally isentropic flow. The energy equation in this calculation model is used to represent the energy change between phases when the relaxation phenomenon happened. However, the change of energy between phases is mainly reflected in the correction of two-phase pressure, which results in the waste of computing resources.This paper simplified the two-phase flow model with seven equations to a numerical calculation method with five equations. A numerical method for predicting the attenuation of pressure fluctuation with the cavitation jet is presented. This method is applied to predict the wave attenuation of the cavitation jet. The main contents are as follows:(1)
The calculation method of velocity and pressure relaxation coefficient is rederived strictly through the interphase action mechanism of the seven-equation two-fluid model.(2)
A five-equation calculation method is constructed.(3)
The internal flow parameters were analyzed to reveal the internal flow trend of the cavitation jet.(4)
By setting the internal switch, the two-coefficient decay model was added into the five-equation calculation method. The calculation method of wave attenuation suitable for water jet cavitation was constructed.


## 2. Two-Phase Flow Model and Properties

### 2.1. Two-Fluid Model

The two-fluid model usually consists of 7-8 equations. Each phase has its density, velocity, and pressure. The two-phase equations are coupled by the interphase forces and relaxation coefficients. In the macroscopic Euler-Euler model, the microscopic parameters in the macroscopic model can be described by adding the void transport equation and the number of bubble concentration equation. The numerical calculation model selected in this paper is shown in the following equations:(1)∂α1∂t+uI∇α1=μP1−P2,(2)∂αkρk∂t+∇⋅αkρkuk=0,(3)∂αkρkuk∂t+∇⋅αkρkuk2+αkPk=PI∇αk±λu2−u1,(4)∂αkρkEk∂t+∇⋅αkρkukEk+αkPkuk=PIuI∇αk±μPI¯P2−P1±λuI¯u2−u1−Kkgfuk,where the subscript k is 1 and 2, in which 1 is the gas phase and 2 is the liquid phase; I is the interphase; α is the void fraction; u is the velocity, m/s; P is the pressure, Pa; ρ is the density, kg/m3; E is the total energy E=e+1/2u2, J; μ is the pressure relaxation coefficient; λ is the velocity relaxation coefficient; and Kkgf is the gravity and resistance term.The calculation method of the interphase boundary was referred to the literature [15, 16]. The expression equation is as follows:(5)PI=PI¯=α1P1+α2P2,(6)uI=uI¯=α1ρ1u1+α2ρ2u2α1ρ1+α2ρ2.


### 2.2. Conventional Closure Relation

The equation of state (EOS) was introduced for the closure of the equation in this model. The function of the EOS is to create a relationship between internal energy, pressure, and density to transform the model.This paper cited a modified form of the stiffened gas in Suarel [18] and Crouzet [12]. The expression for this equation of state is as follows:(7)ekpk,ρk=pk+γkπkρkγk−1+qk,k=1,2,where γk is the specific heat ratio, which is γk=Cp,k/Cv,k. π1 usually takes the value of 0. The expression for π2 is as follows:(8)π2=ρ2Trefγ2−1Cv,2Tref−psatTref,where Tref is the experimental ambient temperature. psat is the saturated vapor pressure.


### 2.3. Improved Relaxation Calculation Method

Relaxation refers to a process in the transient flow where two phases tend to balance between each other [19–21, 23]. In transient flow, as the parameters of two-phase are usually different, the velocity of the wave will be greatly deviated. After each time of transients, the two phases need an equilibrium process to end the difference caused by the transient flow. The relaxation coefficient is an empirical coefficient used to calculate the relaxation process.The construction methods of relaxation coefficient are divided into finite and infinite relaxation coefficient construction methods. In this paper, a new construction method of relaxation coefficients is established by combining the two kinds of relaxation coefficients.During the relaxation process, it is considered that the development of this process is very short, and some literatures indicate that the development process is nearly 1 ms. In the process of interphase relaxation, the influence of spatial variables on the relaxation process can be ignored, and only the change of the partial derivative of the time term is reckoned. Therefore, equations (1)–(4) were sorted out, ignoring the space term and the source term independent of the relaxation process. The equations are modified into the following form:(9)∂α1∂t=μP1−P2,(10)∂P1∂t=μP2−P1PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1,(11)∂P2∂t=μP1−P2PI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2.Equations (9)–(11) could be converted to (26)–(28):(12)α1∗−α10=μP1−P2Δt,(13)P1∗−P10=μP2−P1PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1Δt,(14)P2∗−P20=μP1−P2PI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2Δt,where the variables marked with ∗ indicate the relaxed ones and the variables marked with 0 are the calculated ones. Relaxation process is an interphase equilibrium process. Thus at the end of the relaxation process, there should be P1∗−P2∗=−2S/R or be a tendency to −2S/R. By subtracting equation (13) from equation (14) and substituting P1∗−P2∗=−2S/R into the original system of equations, the following expression can be obtained:(15)P1∗−P2∗=P10−P201+ΔtμPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1,(16)P1∗−P2∗=−2SR,(17)P10−P201+ΔtμPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1=−2SR.The relaxation process of pressure is caused byP20−P10≠0. So the following formulas are required to achieve the entire relaxation process:(18)1+ΔtμPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1=−2SRP10−P20,(19)μPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1=−2SRP10−P20−11Δt.Further finishing could be obtained as follows:(20)μ=−2S/RP10−P20+1PI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1Δt,(21)PI=α1P1+α2P2,(22)2SR≈1∼5×104,Δt≈1×10−5.For the velocity relaxation method, as the relaxation phenomenon is independent of the spatial variables and is only a function in the time direction, the original equations were transformed into the following equations:(23)∂α1ρ1u1∂t=λu2−u1,(24)∂α2ρ2u2∂t=−λu2−u1,(25)∂u1∂t=λu2−u1α1ρ1,(26)∂u2∂t=−λu2−u1α2ρ2.The equation could be transformed as follows:(27)u1∗−u10=Δtλu2−u1α1ρ1,(28)u2∗−u20=−Δtλu2−u1α2ρ2.Further finishing could be obtained as follows:(29)u1∗=Δtλu2−u1α1ρ1+u10,(30)u2∗=−Δtλu2−u1α2ρ2+u20.It is desirable to the relaxation stateu1∗−u2∗=0 or tends to 0 during the transient process. Subtracting the above two equations gives the following expressions:(31)u1∗−u2∗=Δtλu2−u1α1ρ1+Δtλu2−u1α2ρ2+u10−u20,(32)Δtλu2−u1α1ρ1+Δtλu2−u1α2ρ2+u10−u20=0,(33)λ=Δt1α1ρ1+1α2ρ2−1.Therefore, equation (33) is an expression for the improved calculation of the velocity relaxation method.


## 2.1. Two-Fluid Model

The two-fluid model usually consists of 7-8 equations. Each phase has its density, velocity, and pressure. The two-phase equations are coupled by the interphase forces and relaxation coefficients. In the macroscopic Euler-Euler model, the microscopic parameters in the macroscopic model can be described by adding the void transport equation and the number of bubble concentration equation. The numerical calculation model selected in this paper is shown in the following equations:(1)∂α1∂t+uI∇α1=μP1−P2,(2)∂αkρk∂t+∇⋅αkρkuk=0,(3)∂αkρkuk∂t+∇⋅αkρkuk2+αkPk=PI∇αk±λu2−u1,(4)∂αkρkEk∂t+∇⋅αkρkukEk+αkPkuk=PIuI∇αk±μPI¯P2−P1±λuI¯u2−u1−Kkgfuk,where the subscript k is 1 and 2, in which 1 is the gas phase and 2 is the liquid phase; I is the interphase; α is the void fraction; u is the velocity, m/s; P is the pressure, Pa; ρ is the density, kg/m3; E is the total energy E=e+1/2u2, J; μ is the pressure relaxation coefficient; λ is the velocity relaxation coefficient; and Kkgf is the gravity and resistance term.The calculation method of the interphase boundary was referred to the literature [15, 16]. The expression equation is as follows:(5)PI=PI¯=α1P1+α2P2,(6)uI=uI¯=α1ρ1u1+α2ρ2u2α1ρ1+α2ρ2.


## 2.2. Conventional Closure Relation

The equation of state (EOS) was introduced for the closure of the equation in this model. The function of the EOS is to create a relationship between internal energy, pressure, and density to transform the model.This paper cited a modified form of the stiffened gas in Suarel [18] and Crouzet [12]. The expression for this equation of state is as follows:(7)ekpk,ρk=pk+γkπkρkγk−1+qk,k=1,2,where γk is the specific heat ratio, which is γk=Cp,k/Cv,k. π1 usually takes the value of 0. The expression for π2 is as follows:(8)π2=ρ2Trefγ2−1Cv,2Tref−psatTref,where Tref is the experimental ambient temperature. psat is the saturated vapor pressure.


## 2.3. Improved Relaxation Calculation Method

Relaxation refers to a process in the transient flow where two phases tend to balance between each other [19–21, 23]. In transient flow, as the parameters of two-phase are usually different, the velocity of the wave will be greatly deviated. After each time of transients, the two phases need an equilibrium process to end the difference caused by the transient flow. The relaxation coefficient is an empirical coefficient used to calculate the relaxation process.The construction methods of relaxation coefficient are divided into finite and infinite relaxation coefficient construction methods. In this paper, a new construction method of relaxation coefficients is established by combining the two kinds of relaxation coefficients.During the relaxation process, it is considered that the development of this process is very short, and some literatures indicate that the development process is nearly 1 ms. In the process of interphase relaxation, the influence of spatial variables on the relaxation process can be ignored, and only the change of the partial derivative of the time term is reckoned. Therefore, equations (1)–(4) were sorted out, ignoring the space term and the source term independent of the relaxation process. The equations are modified into the following form:(9)∂α1∂t=μP1−P2,(10)∂P1∂t=μP2−P1PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1,(11)∂P2∂t=μP1−P2PI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2.Equations (9)–(11) could be converted to (26)–(28):(12)α1∗−α10=μP1−P2Δt,(13)P1∗−P10=μP2−P1PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1Δt,(14)P2∗−P20=μP1−P2PI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2Δt,where the variables marked with ∗ indicate the relaxed ones and the variables marked with 0 are the calculated ones. Relaxation process is an interphase equilibrium process. Thus at the end of the relaxation process, there should be P1∗−P2∗=−2S/R or be a tendency to −2S/R. By subtracting equation (13) from equation (14) and substituting P1∗−P2∗=−2S/R into the original system of equations, the following expression can be obtained:(15)P1∗−P2∗=P10−P201+ΔtμPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1,(16)P1∗−P2∗=−2SR,(17)P10−P201+ΔtμPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1=−2SR.The relaxation process of pressure is caused byP20−P10≠0. So the following formulas are required to achieve the entire relaxation process:(18)1+ΔtμPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1=−2SRP10−P20,(19)μPI−ρ22∂e2/∂ρ2α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1α1ρ1∂e1/∂P1=−2SRP10−P20−11Δt.Further finishing could be obtained as follows:(20)μ=−2S/RP10−P20+1PI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1Δt,(21)PI=α1P1+α2P2,(22)2SR≈1∼5×104,Δt≈1×10−5.For the velocity relaxation method, as the relaxation phenomenon is independent of the spatial variables and is only a function in the time direction, the original equations were transformed into the following equations:(23)∂α1ρ1u1∂t=λu2−u1,(24)∂α2ρ2u2∂t=−λu2−u1,(25)∂u1∂t=λu2−u1α1ρ1,(26)∂u2∂t=−λu2−u1α2ρ2.The equation could be transformed as follows:(27)u1∗−u10=Δtλu2−u1α1ρ1,(28)u2∗−u20=−Δtλu2−u1α2ρ2.Further finishing could be obtained as follows:(29)u1∗=Δtλu2−u1α1ρ1+u10,(30)u2∗=−Δtλu2−u1α2ρ2+u20.It is desirable to the relaxation stateu1∗−u2∗=0 or tends to 0 during the transient process. Subtracting the above two equations gives the following expressions:(31)u1∗−u2∗=Δtλu2−u1α1ρ1+Δtλu2−u1α2ρ2+u10−u20,(32)Δtλu2−u1α1ρ1+Δtλu2−u1α2ρ2+u10−u20=0,(33)λ=Δt1α1ρ1+1α2ρ2−1.Therefore, equation (33) is an expression for the improved calculation of the velocity relaxation method.


## 3. Simplified Calculation Model

When the model (1)–(4) is used for calculation, it is usually settled as isentropic flow. The equation of energy is primarily used to correct the imbalance of pressure between two phases. When the expression of pressure relaxation coefficient is obtained, the energy equation for the gas and liquid phases can be omitted. The model can be solved only by the void transport equation, continuity equations, and momentum equations. Some literature has pointed out that the velocity of two phases in the two-fluid model can be considered as negligible when the void of the fraction of the gas phase is not large [33]. So the equation can be further simplified as follows:(34)∂α1∂t+u∇α1=μP1−P2,(35)∂αkρk∂t+∇⋅αkρku=0,(36)∂αkρku∂t+∇⋅αkρku2+αkPk=PI∇αk.The density term in the equation can be transformed to the pressure through the relation between density and pressure (37) and (38):(37)dρkdt=rkck2dPkdt,(38)ρk∗=ρk0+rkck2Pk∗−Pk0.The equations could be simplified into the following one-dimensional form to solve the cavitation jet in the pipeline:(39)∂α1∂t+u∂α1∂x=μP1−P2,(40)∂P1∂t+u∂P∂t+ρ1c12r1∂u∂x=c12ρ1μP2−P1r1α1,(41)∂P2∂t+u∂P2∂t+ρ2c22∂u∂x=c22ρ2μP1−P2r2α2,(42)∂u∂t+u∂u∂t+1ρ1∂P1∂x+P1−PIα1ρ1∂α1∂x=−K1gfα1ρ1,(43)∂u∂t+u∂u∂t+1ρ2∂P2∂x+PI−P2α2ρ2∂α1∂x=−K2gfα2ρ2.The system of equations has five distinct eigenvalues, and the equation is unconditional hyperbolic:(44)λ1=uλ2=u+c22r2λ3=u−c22r2λ4=u+c12r1λ5=u−c12r1,where ck represents the wave velocity, and the expression is shown as follows:(45)ck2=Pk/ρk2−∂ek/∂ρk∂ek/∂Pk.


## 4. Numerical Resolution

### 4.1. Hyperbolic Operator

The equations can be written as follows:(46)∂ϕ∂t+∂fϕ∂x=S.The parameter of (46) could be expressed as follows:(47)∂f∂ϕ=a,a=a++a−a+=12a+aa−=12a−a,(48)f=a+ϕ+a−ϕ=f++f−∂ϕ∂t+∂f+∂x+∂f−∂x=S.The numerical calculation method in NND format (49) can be obtained by further sorting out these equations.The NND format is a nonoscillatory, nonfree parameter and dissipative finite difference scheme proposed by Zhang Hanxin [34, 35]. This scheme can effectively suppress the oscillation of numerical solutions. The NND format corrects the equation by appropriately adding a third-order term and adopts different calculation formats for the upstream and downstream regions of the shock wave:(49)∂ϕ∂tjn+fj+1/2L+n+fj+1/2R−n∂x+fj−1/2L+n+fj−1/2R−n∂x=S.Ashj+1/2n=fj+1/2L+n+fj+1/2R−nandhj−1/2n=fj−1/2L+n+fj−1/2R−n, the changed difference format can be expressed as follows:(50)∂ϕ∂tjn=−1Δxhj+12n−hj−12n+S,(51)fj+1/2L+n=fj+n+12minmodΔfj−1/2+n,Δfj−1/2+n,(52)Δfj−1/2L+n=fj−1+n+12minmodΔfj−3/2+n,Δfj−1/2+n,(53)Δfj+1/2R−n=fj+1−n−12minmodΔfj+1/2−n,Δfj+3/2−n,(54)fj−1/2R−n=fj−n−12minmodΔfj−1/2−n,Δfj+1/2−n.The definition of min mod (x, y) has the following meaning: when x and y are the same number, min mod (x, y) takes the value of x, the element with the smallest absolute value in y, and when x and y are different, min mod (x, y) takes a value of zero.


### 4.2. Numerical Resolution for Nonconservative Equation

As the two-fluid model is a typical nonconservative equation, the NND format needs to be corrected for it. For the improving NND format, the Godunov scheme [18, 36–38] (Saurel and Abgrall 1999a; Saurel and Abgrall 1999b; Saurel and Lemetayer 2001; and Saurel et al. 2003) can be used as follows:(55)∂α1∂x=α1j+1/2∗−α1j−1/2∗2⋅Δx,α1j+1/2∗=α1j+1∗−α1j∗,α1j−1/2∗=α1j∗−α1j−1∗.


### 4.3. Boundary Condition

Boundary conditions were constructed by using the method of characteristics. Taking∂ϕ/∂t+B∂ϕ/∂x=S into the form ∂ϕ/∂t+TλT−1∂ϕ/∂x=S and then using the feature matrix left multiplied the equation, the following equation could be obtained:(56)T−1DϕDt=T−1S,∂x∂t=λ.The equation of the boundary condition could be changed to the following form:(57)T−1ϕP−ϕLΔt=T−1SLΔt.The compatibility equation used the following form for interpolation calculation:(58)ϕR=ϕC+ΔtΔxc−uϕC−ϕB,ϕL=ϕC−ΔtΔxc+uϕC−ϕA,ϕL=ϕC−ΔtΔxuϕC−ϕA.


## 4.1. Hyperbolic Operator

The equations can be written as follows:(46)∂ϕ∂t+∂fϕ∂x=S.The parameter of (46) could be expressed as follows:(47)∂f∂ϕ=a,a=a++a−a+=12a+aa−=12a−a,(48)f=a+ϕ+a−ϕ=f++f−∂ϕ∂t+∂f+∂x+∂f−∂x=S.The numerical calculation method in NND format (49) can be obtained by further sorting out these equations.The NND format is a nonoscillatory, nonfree parameter and dissipative finite difference scheme proposed by Zhang Hanxin [34, 35]. This scheme can effectively suppress the oscillation of numerical solutions. The NND format corrects the equation by appropriately adding a third-order term and adopts different calculation formats for the upstream and downstream regions of the shock wave:(49)∂ϕ∂tjn+fj+1/2L+n+fj+1/2R−n∂x+fj−1/2L+n+fj−1/2R−n∂x=S.Ashj+1/2n=fj+1/2L+n+fj+1/2R−nandhj−1/2n=fj−1/2L+n+fj−1/2R−n, the changed difference format can be expressed as follows:(50)∂ϕ∂tjn=−1Δxhj+12n−hj−12n+S,(51)fj+1/2L+n=fj+n+12minmodΔfj−1/2+n,Δfj−1/2+n,(52)Δfj−1/2L+n=fj−1+n+12minmodΔfj−3/2+n,Δfj−1/2+n,(53)Δfj+1/2R−n=fj+1−n−12minmodΔfj+1/2−n,Δfj+3/2−n,(54)fj−1/2R−n=fj−n−12minmodΔfj−1/2−n,Δfj+1/2−n.The definition of min mod (x, y) has the following meaning: when x and y are the same number, min mod (x, y) takes the value of x, the element with the smallest absolute value in y, and when x and y are different, min mod (x, y) takes a value of zero.


## 4.2. Numerical Resolution for Nonconservative Equation

As the two-fluid model is a typical nonconservative equation, the NND format needs to be corrected for it. For the improving NND format, the Godunov scheme [18, 36–38] (Saurel and Abgrall 1999a; Saurel and Abgrall 1999b; Saurel and Lemetayer 2001; and Saurel et al. 2003) can be used as follows:(55)∂α1∂x=α1j+1/2∗−α1j−1/2∗2⋅Δx,α1j+1/2∗=α1j+1∗−α1j∗,α1j−1/2∗=α1j∗−α1j−1∗.


## 4.3. Boundary Condition

Boundary conditions were constructed by using the method of characteristics. Taking∂ϕ/∂t+B∂ϕ/∂x=S into the form ∂ϕ/∂t+TλT−1∂ϕ/∂x=S and then using the feature matrix left multiplied the equation, the following equation could be obtained:(56)T−1DϕDt=T−1S,∂x∂t=λ.The equation of the boundary condition could be changed to the following form:(57)T−1ϕP−ϕLΔt=T−1SLΔt.The compatibility equation used the following form for interpolation calculation:(58)ϕR=ϕC+ΔtΔxc−uϕC−ϕB,ϕL=ϕC−ΔtΔxc+uϕC−ϕA,ϕL=ϕC−ΔtΔxuϕC−ϕA.


## 5. Unsteady Friction Model

Accurate prediction of the water hammer cavitation jet in a pipeline system is the basis of prediction in pipeline obstruction and leakage. In a one-dimensional pipeline system, the unsteady friction model and steady friction model are used to predict the reduction of pressure waves.A one-dimensional pipeline flow resistance model is usually expressed as follows:(59)J=Js+Ju,where Js is the steady friction model. The steady friction can be expressed by the Darcy–Weisbach relation:(60)Ju=fVV2gD.The value of steady-state flow resistance coefficient is related to the Reynolds number. For pipeline flow, when the Reynolds number Re < 2320,f=64/Re, and when Re > 2320, according to the Colebrook–White formula, the following expression can be obtained:(61)1f=−2lg2.51Ref+K/D3.71.In general, the steady friction model cannot well simulate the wave attenuation process in transient flow. Therefore, some scholars put forward the unsteady friction model. Among these unsteady friction models, the two-coefficient decay model has the best application as follows:(62)Ju=Kut∂u∂t+KuxSignuc∂u∂x.For the values ofKut and Kux, one can refer to the literature [31]. The expression of ∂u/∂t and ∂u/∂x in the formula is related to the transfer direction of the characteristic line (63) and (64).Along the positive characteristic,(63)∂u∂x=uin−ui−1nΔx,∂u∂t=uin−uin−1Δt.Along the negative characteristic,(64)∂u∂x=uin−ui+1nΔx∂u∂t=uin−uin−1Δt.As there are few successful cases in the numerical calculation of wave attenuation with cavitation jet. In this paper, the unsteady friction model is partially integrated with the calculation method of the cavitation jet. When the wave experiences 3–5 phases, the friction model gets switch on to the unsteady friction model to guarantee accurate capture of the pressure wave decay.


## 6. Experimental Comparison

### 6.1. Experimental System 1

The length of the experimental system is 37.2 m. The pipeline has a slope of 2 m. The pipe inner diameter is 22 mm. The pipe wall thickness is 1.6 mm. Poisson’s ratio 0.34, Young’s modulus isE = 120 ± 5 GPa, and the wave velocity in the pipe is 1321 m/s. The experiment temperature was 20. The saturated steam pressure was 2340 Pa. The material of the pipe was copper [39].There are three experimental schemes:Experimental scheme 1: the flow rate was 0.28 m/s, the flow direction was p1–p5, the valve closed time was 0.009 s, the absolute head upstream was 28.332 m, and the test time was 0.7 s.Experimental scheme 2: the flow rate was 0.40 m/s, the flow direction was p5–p1, the valve closed time was 0.0095 s, the absolute head upstream was 23.22 m, and the test time was 1.0 s.Experimental scheme 3: the flow rate was 1.50 m/s, the flow direction was p5–p1, the valve closed time was 0.009 s, the absolute head upstream was 23.22 m, and the test time was 1.2 s. The experimental program is shown in Table1, and the schematic diagram of experimental system 1 is shown in Figure 1:(65)c=K/ρ1+KD/Eec′=2.19×109/9981+2.19×109/120×1090.0221/0.0016×1.02=1321m/s.Table 1
Experimental program.
Experimental scheme
Velocity (m/s)
Initial water head (m)
Theoretical head rise (m)
1
0.280
28.33
37.69
2
0.400
23.22
53.84
3
1.500
22.66
201.90Figure 1
Schematic diagram of experimental system 1.


### 6.2. Experimental System 2

The length of the experimental system is 36 m. The pipeline has a slope of 1 m. The pipe inner diameter is 19.5 mm. The pipe wall thickness is 1.588 mm. Poisson’s ratio 0.34, Young’s modulus isE = 119 ± 5 GPa, and the wave velocity in the pipe is 1280 m/s. The experiment temperature was 23.9. The saturated steam pressure was 2320 Pa. The material of the pipe was copper [7]. The experimental program is shown in Table 2, and the schematic diagram of experimental system 1 is shown in Figure 2:(66)c=K/ρ1+KD/Eec′=2.19×109/997.31+2.19×109/119×1090.01905/0.0016×0.88=1280m/s.Table 2
Experimental program.
Experimental scheme
Velocity (m/s)
Initial water head (m)
Valve closed time (s)
Theoretical head rise (m)
4
0.332
33.55
0.03608
43.33
5
0.401
33.54
0.01848
52.34
6
0.466
33.32
0.03520
60.82
7
0.696
33.05
0.02640
90.84
8
0.938
36.68
0.03608
122.43
9
1.125
31.88
0.04400
146.83Figure 2
Schematic diagram of experimental system 2.


## 6.1. Experimental System 1

The length of the experimental system is 37.2 m. The pipeline has a slope of 2 m. The pipe inner diameter is 22 mm. The pipe wall thickness is 1.6 mm. Poisson’s ratio 0.34, Young’s modulus isE = 120 ± 5 GPa, and the wave velocity in the pipe is 1321 m/s. The experiment temperature was 20. The saturated steam pressure was 2340 Pa. The material of the pipe was copper [39].There are three experimental schemes:Experimental scheme 1: the flow rate was 0.28 m/s, the flow direction was p1–p5, the valve closed time was 0.009 s, the absolute head upstream was 28.332 m, and the test time was 0.7 s.Experimental scheme 2: the flow rate was 0.40 m/s, the flow direction was p5–p1, the valve closed time was 0.0095 s, the absolute head upstream was 23.22 m, and the test time was 1.0 s.Experimental scheme 3: the flow rate was 1.50 m/s, the flow direction was p5–p1, the valve closed time was 0.009 s, the absolute head upstream was 23.22 m, and the test time was 1.2 s. The experimental program is shown in Table1, and the schematic diagram of experimental system 1 is shown in Figure 1:(65)c=K/ρ1+KD/Eec′=2.19×109/9981+2.19×109/120×1090.0221/0.0016×1.02=1321m/s.Table 1
Experimental program.
Experimental scheme
Velocity (m/s)
Initial water head (m)
Theoretical head rise (m)
1
0.280
28.33
37.69
2
0.400
23.22
53.84
3
1.500
22.66
201.90Figure 1
Schematic diagram of experimental system 1.


## 6.2. Experimental System 2

The length of the experimental system is 36 m. The pipeline has a slope of 1 m. The pipe inner diameter is 19.5 mm. The pipe wall thickness is 1.588 mm. Poisson’s ratio 0.34, Young’s modulus isE = 119 ± 5 GPa, and the wave velocity in the pipe is 1280 m/s. The experiment temperature was 23.9. The saturated steam pressure was 2320 Pa. The material of the pipe was copper [7]. The experimental program is shown in Table 2, and the schematic diagram of experimental system 1 is shown in Figure 2:(66)c=K/ρ1+KD/Eec′=2.19×109/997.31+2.19×109/119×1090.01905/0.0016×0.88=1280m/s.Table 2
Experimental program.
Experimental scheme
Velocity (m/s)
Initial water head (m)
Valve closed time (s)
Theoretical head rise (m)
4
0.332
33.55
0.03608
43.33
5
0.401
33.54
0.01848
52.34
6
0.466
33.32
0.03520
60.82
7
0.696
33.05
0.02640
90.84
8
0.938
36.68
0.03608
122.43
9
1.125
31.88
0.04400
146.83Figure 2
Schematic diagram of experimental system 2.


## 7. Discussion

The results of the 7- and 5-equation models are compared using these nine experiments. The time-effectiveness, calculation accuracy, residual void of gas, and internal velocity characteristics of the two models are analyzed.


### 7.1. Time-Effectiveness

The 7-equation model, 5-equation model, and method of characteristic (MOC) [4, 6, 7] are selected for the comparison for time-effectiveness. The time step and the number of grids are guaranteed equality in the calculation. Nine experiments were chosen for comparative verification. The calculated time is statistically calculated and plotted as shown in Figure 3. In Figure 3, MOC is used as the measurement standard of unit 1, and the data in the 7-equation model and the 5-equation model are multiples of the calculation time of the MOC.Figure 3
Comparison diagram of computing time-effectiveness of the three models.


### 7.2. Calculation Accuracy

According to the calculation principle of the MOC method, it does not have the capacity to carry out accurate numerical simulation of the two-phase flow when cavitation jet occurs. Therefore, the calculation accuracy of the MOC is low. However, the MOC is easily programmable, so it has high time-effectiveness.The calculation results of the three models are compared, as shown in Figures4–7. In the simulation of 1–3 wave phases, both 7-equation and 5-equation models could well capture the wave changes of jet cavitation.Figure 4
Comparison of scheme 3 calculation scheme (endpoint of the system).Figure 5
Comparison of scheme 3 calculation scheme (midpoint of the system).Figure 6
Comparison of scheme 9 calculation scheme (endpoint of the system).Figure 7
Comparison of scheme 9 calculation scheme (a quarter of the system).It can be seen that the 5-equation model has a good calculation time-effectiveness on the basis of guaranteeing the calculation accuracy.


### 7.3. Residual Void of Gas

The main reason for the inaccurate calculation of the model is the residual void of the gas. So the distribution void of fraction in the 7-equation model and 5-equation model is presented. As shown in Figures8 and 9, the 5-equation model can better eliminate residual void of the gas after each iteration. From equations (13), (14) and (40), (41), it can be seen that the 5-equation model reduce the interference between phase parameters and thermodynamic parameters, making the void of fraction simulation more accurate.Figure 8
Comparison of the cavitation rate of scheme 3 (endpoint of the system).Figure 9
Comparison of the cavitation rate of scheme 9 (endpoint of the system).


### 7.4. Lag Time Characteristic Analysis

The internal response curves of pressure, velocity, and void of fraction of schemes 3 and 9 were taken for analysis. So the wave phase and period of normal launching cavitation jet are taken for theoretical analysis. In Table3, the wave velocity is considered as a constant value.Table 3
Wave phase and cycle of scheme 3 and scheme 9
Scheme
Wave velocityC (m/s)
Wave phase 2 L/c (s)
Wave cycle 4 L/c (s)
The actual cycle (s)
Deviation (s)
3
1320
0.0564
0.1128
0.5713
0.4585
9
1280
0.0563
0.1125
0.3133
0.2008It can be seen in Table3 that the larger the intensity of the cavitation jet is, the longer the lag time (deviation) of the wave is. The variation trend of velocity and pressure near the tail position is displayed (under the 5-equation model), as shown in Figures 10–13.Figure 10
Change of internal velocity in scheme 3.Figure 11
Change of internal velocity in scheme 9.Figure 12
The internal velocity of the calculation scheme with no cavitation jet.Figure 13
The internal pressure of the calculation scheme with no cavitation jet.As shown in Figures12 and 13, the internal velocity distribution forming is divided into six stages. The six calculation schemes are, respectively, divided into the following five processes for the calculation scheme of cavitation-free jet:(1)
In the reverses stage of 1-2, the liquid starts to move in the opposite direction of the flow, and the pressure starts to gradually decrease(2)
In the recoveries stage of 2-3, the backflow starts to slow down gradually, and the velocity starts to gradually reduce and stops after reaching 0(3)
In the stabilizations stage of 3-4, the velocity is 0(4)
In the acceleration stage of 4-5, it begins to accelerate in the same direction with the initial flow, and the pressure starts to gradually recover(5)
In the deceleration stage of 5-6, the velocity gradually decreases, and pressure gradually increasesThe above five stages are repeated until the water hammer wave attenuation disappears.As shown in Figures10 and 11, the internal velocity distribution of the formation of the cavitation jet is divided into seven stages, and the six processes are as follows:(1)
In the stage of 1-2 flow reverses, the liquid starts to move in the opposite direction of the flow direction, and the pressure starts to gradually decrease.(2)
In the stage of 2-3 flow recoveries, the backflow began to slow down gradually, and the velocity gradually decreased but does not return to 0. The void of the cavitation jet is still forming.(3)
In the stage of 3-4 flow recoveries and acceleration, the velocity in the first half slows down, while in the second half accelerates. In the whole process, the slope is greater than in stage 2-3. At this point, the void generated by the cavitation jet starts to disappear rapidly, and the pressure is assured to remain unchanged.(4)
In the gentle flow stage of 4-5, the slope of the forward accelerated flow is lower and the duration of time is shorter than that in the 3-4 stage.(5)
In the acceleration stage of 5-6, it continues to accelerate in a positive direction, and the pressure is guaranteed to remain unchanged.(6)
In the flow deceleration stage of 6-7, the velocity gradually decreases and pressure gradually increases.It can be seen that the main differences lie in the 3-4 stage and 4-5 stage, in which the velocity increases significantly, leading to the generation of voids and the lag of time characteristics in the calculation.


### 7.5. Analysis of Pressure Amplitude Characteristics

It can be seen from the difference in the distribution of 3-4, 4-5, and 5-6 stages on the curve that these three stages play relatively different roles in the change of cavitation. The two-fluid model is taken as an example, where the transport equation of the void can be written as follows:(67)∂α1∂t+uI∇⋅α1=μP1−P2+m˙ρI.When only considering the mathematical properties of two-phase coupling (ignoring variables such as viscosity and surface tension), the generation of the void is generated from the following two aspects:μP1−P2: gas volume changes due to two-phase pressure changesm˙/ρI: the change of cavity volume caused by phase transitionIt is considered that at the end of each calculation, the two-phase pressure tends to be equal, which isP1∗−P2∗=0, the relaxation coefficient of the pressure changes as follows:(68)μ=1PI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1Δt,(69)μP1−P2=P1−P2PI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1Δt.The above equations are expressed as the change of gas void caused by the change of pressure.The change of cavity volume caused by phase transition is expressed as follows:(70)m˙/ρI=2S/RPI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1/α1ρ1∂e1/∂P1Δt,where 2S/R is usually greater than P1−P2. According to the slope of the curve, it can be concluded that the 4-5 stage is usually the stage of cavitation condensation and collapse, while the 5-6 stage is the gas volume change caused by the two-phase pressure.It can be found that in the experimental scheme 1 with low cavitation intensity (theoretical head change 43.33 m), the actual pressure head change reaches 67.58 m, but in the experimental scheme 6 with low cavitation intensity (theoretical head change 146.83 m), the actual pressure head change is only 144.98 m. Therefore, the internal flow characteristics of the two curves were analyzed separately, as shown in Figures14 and 15.It can be seen in Figure14 that the 4-5 stage of change in the internal characteristics of velocity is steeper. Stage 4-5 corresponds to the elimination stage of cavitation. The main reason for high pressure is that the cavitation formed in a small size, so it quickly collapses under high pressure. After cavitation collapse, the high-pressure liquid forms jet flow, which rapidly impacts the wall and causes the pressure to increase instantly.Figure 14
Scheme 1 calculation of internal velocity change of the scheme.However, it is different from the velocity distribution trend in Figure14 (scheme 1). The change rate of stage 4-5 more flattens than that in stages 3-4 and 5-6 in Figure 15. In scheme 6, the cavitation volume is larger, which results in longer elimination time. The slope is smaller, and the pressure increase is not as large as in scheme 1. This explains the reasons why the intensity of the cavitation jet is low and the pressure rise is high. The pressure difference between phases in Figures 16 and 17 also confirms the phenomenon that after cavitation collapse, high-pressure liquid forms jet flow and rapidly impacts the wall, leading to an instantaneous increase in pressure.Figure 15
Scheme 6 calculation of internal velocity change of the scheme.Figure 16
Comparison of scheme 1 pressure calculation scheme and pressure difference between two phases.Figure 17
Comparison of scheme 6 pressure calculation scheme and pressure difference between two phases.Stage 5-6 corresponds to the gas compression stage. Since the gas modulus is larger than that of the water vapor, the slope of the 5-6 phase is lower.


### 7.6. Unsteady Friction Model

When the steady flow resistance model is adopted, only the first 2-3 waves of cavitation jet pressure can be predicted accurately. The whole attenuation and propagation process of the pressure wave cannot be analyzed and demonstrated in detail. Therefore, it is necessary to point out the unsteady flow resistance model to the calculation method for analysis. Calculation examples of scheme 7 and scheme 8 are adopted for verification.In calculation scheme 7, the following are observed:(1)
In the calculation method without the addition of an unsteady flow friction model, the simulation accuracy of first and second pressure waves is good, but the subsequent pressure attenuation process is far from the experiment(2)
In the calculation method with the addition of an unsteady flow friction model, the amplitude deviation of the first and second pressure wave simulation is relatively large, but the subsequent pressure attenuation process is in good agreement with the experimentIn calculation scheme 8, the following are observed:(1)
In the calculation method without the addition of an unsteady flow friction model, the simulation accuracy of the first and second pressure waves is good, but the subsequent pressure attenuation process is more different from the experiment(2)
In the calculation method with an unsteady flow friction model, the simulated amplitude and phase of the first and secondary pressure waves have deviated, but the subsequent pressure attenuation process is higher than the experimental valueIn conclusion, for the process of pipeline cavitation jet, the applicability of the two-coefficient unsteady flow resistance model decreases as the severity of cavitation jet increases. Therefore, the phase deviation and the calculation error of the pressure attenuation process are caused. To simulate the full attenuation process of cavitation jet pressure wave reasonably, the double-coefficient flow resistance model needs to be re-integrated.


### 7.7. Application of Improved Transient Resistance Model

It can be seen from Figures18 and 19 and Figures 20 and 21 that to achieve accurate simulation of the process with cavitation jet, a switch needs to be set in the calculation. This switch does not turn on the unsteady flow resistance model after 3-4∗ (2 L/c) (wave phase). In this case, the model can not only meet the pressure and phase accuracy requirements of the first-two peaks but also ensure a good simulation of pressure wave attenuation. The simulation results are shown in Figures 22–25.Figure 18
Comparison of calculation results of scheme 7 with and without unsteady flow friction model (endpoint of the system).Figure 19
Comparison of calculation results of scheme 7 with and without unsteady flow friction model (midpoint).Figure 20
Comparison of calculation results of scheme 8 with and without unsteady flow friction model (endpoint of the system).Figure 21
Comparison of calculation results of scheme 8 with and without unsteady flow friction model (midpoint).Figure 22
Comparison of calculation results of scheme 7 improved unsteady flow friction model (endpoint of the system).Figure 23
Comparison of calculation results of scheme 7 improved unsteady flow friction model (midpoint).Figure 24
Comparison of calculation results of scheme 8 improved unsteady flow friction model (endpoint of the system).Figure 25
Comparison of calculation results of scheme 8 improved unsteady flow friction model (midpoint).


## 7.1. Time-Effectiveness

The 7-equation model, 5-equation model, and method of characteristic (MOC) [4, 6, 7] are selected for the comparison for time-effectiveness. The time step and the number of grids are guaranteed equality in the calculation. Nine experiments were chosen for comparative verification. The calculated time is statistically calculated and plotted as shown in Figure 3. In Figure 3, MOC is used as the measurement standard of unit 1, and the data in the 7-equation model and the 5-equation model are multiples of the calculation time of the MOC.Figure 3
Comparison diagram of computing time-effectiveness of the three models.


## 7.2. Calculation Accuracy

According to the calculation principle of the MOC method, it does not have the capacity to carry out accurate numerical simulation of the two-phase flow when cavitation jet occurs. Therefore, the calculation accuracy of the MOC is low. However, the MOC is easily programmable, so it has high time-effectiveness.The calculation results of the three models are compared, as shown in Figures4–7. In the simulation of 1–3 wave phases, both 7-equation and 5-equation models could well capture the wave changes of jet cavitation.Figure 4
Comparison of scheme 3 calculation scheme (endpoint of the system).Figure 5
Comparison of scheme 3 calculation scheme (midpoint of the system).Figure 6
Comparison of scheme 9 calculation scheme (endpoint of the system).Figure 7
Comparison of scheme 9 calculation scheme (a quarter of the system).It can be seen that the 5-equation model has a good calculation time-effectiveness on the basis of guaranteeing the calculation accuracy.


## 7.3. Residual Void of Gas

The main reason for the inaccurate calculation of the model is the residual void of the gas. So the distribution void of fraction in the 7-equation model and 5-equation model is presented. As shown in Figures8 and 9, the 5-equation model can better eliminate residual void of the gas after each iteration. From equations (13), (14) and (40), (41), it can be seen that the 5-equation model reduce the interference between phase parameters and thermodynamic parameters, making the void of fraction simulation more accurate.Figure 8
Comparison of the cavitation rate of scheme 3 (endpoint of the system).Figure 9
Comparison of the cavitation rate of scheme 9 (endpoint of the system).


## 7.4. Lag Time Characteristic Analysis

The internal response curves of pressure, velocity, and void of fraction of schemes 3 and 9 were taken for analysis. So the wave phase and period of normal launching cavitation jet are taken for theoretical analysis. In Table3, the wave velocity is considered as a constant value.Table 3
Wave phase and cycle of scheme 3 and scheme 9
Scheme
Wave velocityC (m/s)
Wave phase 2 L/c (s)
Wave cycle 4 L/c (s)
The actual cycle (s)
Deviation (s)
3
1320
0.0564
0.1128
0.5713
0.4585
9
1280
0.0563
0.1125
0.3133
0.2008It can be seen in Table3 that the larger the intensity of the cavitation jet is, the longer the lag time (deviation) of the wave is. The variation trend of velocity and pressure near the tail position is displayed (under the 5-equation model), as shown in Figures 10–13.Figure 10
Change of internal velocity in scheme 3.Figure 11
Change of internal velocity in scheme 9.Figure 12
The internal velocity of the calculation scheme with no cavitation jet.Figure 13
The internal pressure of the calculation scheme with no cavitation jet.As shown in Figures12 and 13, the internal velocity distribution forming is divided into six stages. The six calculation schemes are, respectively, divided into the following five processes for the calculation scheme of cavitation-free jet:(1)
In the reverses stage of 1-2, the liquid starts to move in the opposite direction of the flow, and the pressure starts to gradually decrease(2)
In the recoveries stage of 2-3, the backflow starts to slow down gradually, and the velocity starts to gradually reduce and stops after reaching 0(3)
In the stabilizations stage of 3-4, the velocity is 0(4)
In the acceleration stage of 4-5, it begins to accelerate in the same direction with the initial flow, and the pressure starts to gradually recover(5)
In the deceleration stage of 5-6, the velocity gradually decreases, and pressure gradually increasesThe above five stages are repeated until the water hammer wave attenuation disappears.As shown in Figures10 and 11, the internal velocity distribution of the formation of the cavitation jet is divided into seven stages, and the six processes are as follows:(1)
In the stage of 1-2 flow reverses, the liquid starts to move in the opposite direction of the flow direction, and the pressure starts to gradually decrease.(2)
In the stage of 2-3 flow recoveries, the backflow began to slow down gradually, and the velocity gradually decreased but does not return to 0. The void of the cavitation jet is still forming.(3)
In the stage of 3-4 flow recoveries and acceleration, the velocity in the first half slows down, while in the second half accelerates. In the whole process, the slope is greater than in stage 2-3. At this point, the void generated by the cavitation jet starts to disappear rapidly, and the pressure is assured to remain unchanged.(4)
In the gentle flow stage of 4-5, the slope of the forward accelerated flow is lower and the duration of time is shorter than that in the 3-4 stage.(5)
In the acceleration stage of 5-6, it continues to accelerate in a positive direction, and the pressure is guaranteed to remain unchanged.(6)
In the flow deceleration stage of 6-7, the velocity gradually decreases and pressure gradually increases.It can be seen that the main differences lie in the 3-4 stage and 4-5 stage, in which the velocity increases significantly, leading to the generation of voids and the lag of time characteristics in the calculation.


## 7.5. Analysis of Pressure Amplitude Characteristics

It can be seen from the difference in the distribution of 3-4, 4-5, and 5-6 stages on the curve that these three stages play relatively different roles in the change of cavitation. The two-fluid model is taken as an example, where the transport equation of the void can be written as follows:(67)∂α1∂t+uI∇⋅α1=μP1−P2+m˙ρI.When only considering the mathematical properties of two-phase coupling (ignoring variables such as viscosity and surface tension), the generation of the void is generated from the following two aspects:μP1−P2: gas volume changes due to two-phase pressure changesm˙/ρI: the change of cavity volume caused by phase transitionIt is considered that at the end of each calculation, the two-phase pressure tends to be equal, which isP1∗−P2∗=0, the relaxation coefficient of the pressure changes as follows:(68)μ=1PI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1Δt,(69)μP1−P2=P1−P2PI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1Δt.The above equations are expressed as the change of gas void caused by the change of pressure.The change of cavity volume caused by phase transition is expressed as follows:(70)m˙/ρI=2S/RPI−ρ22∂e2/∂ρ2/α2ρ2∂e2/∂P2−PI−ρ12∂e1/∂ρ1/α1ρ1∂e1/∂P1/α1ρ1∂e1/∂P1Δt,where 2S/R is usually greater than P1−P2. According to the slope of the curve, it can be concluded that the 4-5 stage is usually the stage of cavitation condensation and collapse, while the 5-6 stage is the gas volume change caused by the two-phase pressure.It can be found that in the experimental scheme 1 with low cavitation intensity (theoretical head change 43.33 m), the actual pressure head change reaches 67.58 m, but in the experimental scheme 6 with low cavitation intensity (theoretical head change 146.83 m), the actual pressure head change is only 144.98 m. Therefore, the internal flow characteristics of the two curves were analyzed separately, as shown in Figures14 and 15.It can be seen in Figure14 that the 4-5 stage of change in the internal characteristics of velocity is steeper. Stage 4-5 corresponds to the elimination stage of cavitation. The main reason for high pressure is that the cavitation formed in a small size, so it quickly collapses under high pressure. After cavitation collapse, the high-pressure liquid forms jet flow, which rapidly impacts the wall and causes the pressure to increase instantly.Figure 14
Scheme 1 calculation of internal velocity change of the scheme.However, it is different from the velocity distribution trend in Figure14 (scheme 1). The change rate of stage 4-5 more flattens than that in stages 3-4 and 5-6 in Figure 15. In scheme 6, the cavitation volume is larger, which results in longer elimination time. The slope is smaller, and the pressure increase is not as large as in scheme 1. This explains the reasons why the intensity of the cavitation jet is low and the pressure rise is high. The pressure difference between phases in Figures 16 and 17 also confirms the phenomenon that after cavitation collapse, high-pressure liquid forms jet flow and rapidly impacts the wall, leading to an instantaneous increase in pressure.Figure 15
Scheme 6 calculation of internal velocity change of the scheme.Figure 16
Comparison of scheme 1 pressure calculation scheme and pressure difference between two phases.Figure 17
Comparison of scheme 6 pressure calculation scheme and pressure difference between two phases.Stage 5-6 corresponds to the gas compression stage. Since the gas modulus is larger than that of the water vapor, the slope of the 5-6 phase is lower.


## 7.6. Unsteady Friction Model

When the steady flow resistance model is adopted, only the first 2-3 waves of cavitation jet pressure can be predicted accurately. The whole attenuation and propagation process of the pressure wave cannot be analyzed and demonstrated in detail. Therefore, it is necessary to point out the unsteady flow resistance model to the calculation method for analysis. Calculation examples of scheme 7 and scheme 8 are adopted for verification.In calculation scheme 7, the following are observed:(1)
In the calculation method without the addition of an unsteady flow friction model, the simulation accuracy of first and second pressure waves is good, but the subsequent pressure attenuation process is far from the experiment(2)
In the calculation method with the addition of an unsteady flow friction model, the amplitude deviation of the first and second pressure wave simulation is relatively large, but the subsequent pressure attenuation process is in good agreement with the experimentIn calculation scheme 8, the following are observed:(1)
In the calculation method without the addition of an unsteady flow friction model, the simulation accuracy of the first and second pressure waves is good, but the subsequent pressure attenuation process is more different from the experiment(2)
In the calculation method with an unsteady flow friction model, the simulated amplitude and phase of the first and secondary pressure waves have deviated, but the subsequent pressure attenuation process is higher than the experimental valueIn conclusion, for the process of pipeline cavitation jet, the applicability of the two-coefficient unsteady flow resistance model decreases as the severity of cavitation jet increases. Therefore, the phase deviation and the calculation error of the pressure attenuation process are caused. To simulate the full attenuation process of cavitation jet pressure wave reasonably, the double-coefficient flow resistance model needs to be re-integrated.


## 7.7. Application of Improved Transient Resistance Model

It can be seen from Figures18 and 19 and Figures 20 and 21 that to achieve accurate simulation of the process with cavitation jet, a switch needs to be set in the calculation. This switch does not turn on the unsteady flow resistance model after 3-4∗ (2 L/c) (wave phase). In this case, the model can not only meet the pressure and phase accuracy requirements of the first-two peaks but also ensure a good simulation of pressure wave attenuation. The simulation results are shown in Figures 22–25.Figure 18
Comparison of calculation results of scheme 7 with and without unsteady flow friction model (endpoint of the system).Figure 19
Comparison of calculation results of scheme 7 with and without unsteady flow friction model (midpoint).Figure 20
Comparison of calculation results of scheme 8 with and without unsteady flow friction model (endpoint of the system).Figure 21
Comparison of calculation results of scheme 8 with and without unsteady flow friction model (midpoint).Figure 22
Comparison of calculation results of scheme 7 improved unsteady flow friction model (endpoint of the system).Figure 23
Comparison of calculation results of scheme 7 improved unsteady flow friction model (midpoint).Figure 24
Comparison of calculation results of scheme 8 improved unsteady flow friction model (endpoint of the system).Figure 25
Comparison of calculation results of scheme 8 improved unsteady flow friction model (midpoint).


## 8. Conclusion

(1)
In this paper, a more efficient simplified numerical calculation method is established by simplifying the seven-equation two-fluid model, which not only has high computational efficiency but also can further eliminate the problems existing in the two-fluid model.
(2)
It is found that when a cavitation jet is formed, the velocity flow characteristics are different from the traditional ones. These differences result in different response time and response amplitude of water hammer pressure.
(3)
Coupled with the two-coefficient transient flow friction model, a numerical calculation method suitable for predicting the attenuation of the cavitation jet pressure wave is found. By opening the transient friction model when the water hammer wave transmits 3–5 phases, this method can simulate the reduction of wave well. It provides a reference for the research of pressure pipeline leakage, vibration, and noise.



---

*Source: 4150691-2020-02-11.xml*