Establishment and Solution of Four Variable Water Hammer Mathematical Model for Conveying Pipe

Abstract

Transient flow in pipe is a much debated topic in the field of hydrodynamics. The water hammer effect caused by instantaneous valve closing is an important branch of transient flow. At present, the fluid density is regarded as a constant in the study of the water hammer effect in pipe. When there is gas in the pipe, the variation range of density is large, and the pressure-wave velocity should also change continuously along the pipe. This study considers the interaction between pipeline fluid motion and water hammer wave propagation based on the essence of water hammer, with the pressure, velocity, density and overflow area set as variables. A new set of water hammer calculation equations was deduced and solved numerically. The effects of different valve closing time, flow rate and gas content on pressure distribution and the water hammer effect were studied. It was found that with the increase in valve closing time, the maximum fluctuating pressure at the pipe end decreased, and the time of peak value also lagged behind. When the valve closing time increased from 5 s to 25 s, the difference in water hammer pressure was 0.72 MPa, and the difference in velocity fluctuation amplitude was 0.076 m/s. The findings confirm: the greater the flow, the greater the pressure change at the pipe end; the faster the speed change, the more obvious the water hammer effect. High-volume flows were greatly disturbed by instantaneous obstacles such as valve closing. With the increase of time, the pressure fluctuation gradually attenuated along the pipe length. The place with the greatest water hammer effect was near the valve. Under the coupling effect of time and tube length, the shorter the time and the shorter the tube length, the more obvious the pressure fluctuation. Findings also confirm: the larger the gas content, the smaller the fluctuation peak of pipe end pressure; the longer the water hammer cycle, the smaller the pressure-wave velocity. The actual pressure fluctuation value was obviously lower than that without gas, and the size of the pressure wave mainly depended on the gas content. When the gas content increased from 1% to 9%, the difference of water hammer pressure was 0.41 MPa.

1. Introduction

Since the mid nineteenth century, the unsteady flow of fluid has been a much discussed topic for scholars. In steady flow, the parameter at any point does not change with time; while in unsteady flow, the parameter at any point is a function of time and changes with it [1,2]. When boundary conditions at the end of a pipeline change, the steady flow in the pipeline will also change accordingly, manifested as a sudden change in the flow rate and pressure of the fluid, which causes momentum conversion. The occurrence time of this phenomenon is generally very short, so it is called transient flow, also known as water hammer or water hammer [3,4]. At present, the most representative work on water hammer is “Transient Flow” written by Professor Streeter, an American hydraulicsist [5]. This monograph systematically explains the types of water hammer process, the mechanism and derivation of water hammer, the applicable conditions of the model and the corresponding solution method, etc. Today, the method of characteristics (MOC), first proposed in the monograph, are widely used to characterize the water hammer process.

The classical transient flow models of pressurized pipe flow all use a reservoir-pipe-valve system as the research object and close the valve as the research condition [6,7]. The mathematical model includes the motion equation and the continuity equation. The pipeline element is selected as the control body. The momentum equation is established by the momentum change of the fluid entering and exiting the control body, and the continuity equation is established by the mass conservation of the fluid entering and exiting the control body. When the valve is closed, a short section of fluid at the valve stops flowing instantaneously, and the next section of fluid next to this section of fluid and in the upstream direction is blocked, decelerating to a standstill. By analogy, the fluid in the tube gradually stops moving along the upstream direction. Since the kinetic energy of the fluid is reduced and converted into pressure energy, the pressure in the tube gradually increases along the upstream direction. That is, a pressure wave propagates upstream at a wave velocity a [8,9]. Figure 1 shows the movement of the liquid in a horizontal pipeline after the valve is suddenly closed at the end of a pipeline, ignoring the friction and local loss.

In recent years, many scholars have performed research on the transient flow in fluid pipelines. Baba [11] considered the numerical calculation of the transient flow model of a hydrogen-natural gas mixture in pipelines with different types of valves. The influence of the valve function was analyzed, and the calculation results under different mass ratios and gas mixture flow parameters were given. Wahba [12] used one-dimensional and two-dimensional water hammer models to numerically study the attenuation of turbulence in pipelines. The one-dimensional model was solved by the characteristic line method, and the two-dimensional model was solved by the semidiscretization method and the fourth-order Runge—Kutta method. Kandil [13] used the characteristic line method to numerically solve the water hammer equation and studied pressure fluctuations by changing the elastic modulus of the pipeline and Poisson’s ratio. Hossam [14] first used the wave characteristic line method to numerically simulate the transient flow in a viscoelastic pipe, re-derived the Joukowsky equation of the elastic pipe, and derived the quadratic equation of energy conservation in the viscoelastic pipe. Khamoushi [15] proposed a one-dimensional model of transient flow in non-Newtonian fluids, using Zielke’s unsteady friction solution for power law fluids and cross fluids. The time increment updated the weight function of the Zielke model. Abdeldayem [16] studied and explored the performance of different unsteady friction models in terms of accuracy, efficiency, and reliability to determine the most suitable model for engineering practice. Through comparison, it was found that Vítkovský’s unsteady friction model is the most suitable. Taking the reservoir-pipe-valve system as an example, the effectiveness of the method was verified by experimental data. Lashkarbolok [17] used radial-basis functions and least-squares optimization technology to solve the numerical solution of one-dimensional equations governing transient pipe flow. The method can deal with geometrically nonuniform pipes by employing an arbitrary distribution of scattering nodes in the space domain. Aliabadi [18] studied fluid-structure interaction (FSI) during water hammer in a viscoelastic (VE) pipe in the frequency domain. The main aim was to investigate pressure and stress waves using the transfer matrix method (TMM) in a typical reservoir-VE pipe-valve (RPV) system. Both major coupling mechanisms, namely Poisson and junction, are taken into account. Firouzi [19] proposed a time-varying reliability method to predict the failure risk of pipeline strength reduction due to the combination of water hammer pressure and corrosion and proposed the coupling model of pitting corrosion and water hammer pressure. It was found that with an increase in the average incidence of water hammer, the risk of pipeline failure increases, and the sudden closure of pipeline valves will have a destructive impact on the safety of a pipe network. Henclik [20] proposed an original concept of applying a shock response spectrum (SRS) method to water hammer events. The method is shortly presented, its numerical approach developed, and its application to water hammer loads described. SRS computations and analysis of water hammer experimental results measured at a laboratory pipeline are performed and concluded. Chen [21] constructed a fluid-solid coupling 4 equation model applied to high-temperature and high-pressure oil and gas based on the fluid-solid coupling theory and selected the best coupling method and combined the characteristic line method to obtain the pressure-wave velocity. Zhu [22,23] conducted a transient flow experiment using a sudden valve closing of a gas-liquid two-phase pressurized pipe flow and used the experiment to verify that the pressure-wave velocity of a pressurized pipe flow transient flow model was corrected by the gas content rate. At the same time, considering the nonconstant friction and the viscoelasticity of the pipe, the influence on gas flow and pressure-wave dissipation was studied, and the relevant parameters were checked. Fu [24,25] used the split-phase flow model to simulate the liquid column separation caused by transient flow in the pressure pipe flow and also considered the influence of the pipe stress wave on the pressure wave. Due to the existence of nonconservative terms in the model and to ensure the convergence of the results, this model needed to be solved separately while relaxation parameters were introduced. Han [26] examined bridging the water hammer in long-distance pipelines. He based his study on classic water hammer theory, considering the growth and reduction of cavities, the air content in the pipeline, and the water hammer after the instantaneous change in the pressure in the corrugated tube, among other factors. The process was described in semianalytical mathematics, and the influence of the growth and reduction time of cavities, the size of the gas content, and the relative elevation of the highest point of the pipeline to the horizontal section on bridging the water hammer pressure was analyzed.

After investigation, it was found that a large number of domestic and foreign literatures regard the fluid density in the conveying pipe as a constant. When the fluid in the tube is all liquid, the error caused by the fluid density as a constant is within the engineering allowable range. However, when there is gas (natural gas, etc.) in the tube, the gas has different degrees of compressibility and the density varies widely. Treating it as a constant violates the actual process of fluid flow in the tube. The density is closely related to the pressure-wave speed, which should also be constantly changing along the pipeline.

Furthermore, when calculating water hammer pressure, earlier researchers did not consider the change and inconsistency of the research object at the beginning and end of the period due to fluid inflow and outflow in the dt period. Based on the essence of water hammer, the interaction between pipeline fluid motion and water hammer wave propagation is considered in this paper, and the pressure, velocity, density and overflow area are set as variables. A new set of water hammer calculation equations is deduced and solved numerically. The effects of different valve closing time, flow rate and gas content on pressure distribution and water hammer effect are studied. Such calculations of the propagation process of the water hammer wave are more reasonable and reliable.

2. Mathematical Model

2.1. Fluid Motion Equation

Taking the water body in the flow section of length x as the research object, this paper tracks it to study the water body in the flow section. The analysis diagram of the interaction between the pipe water movement and the water hammer wave propagation is shown in Figure 2. The initial length of the stream segment is set to x. The liquid in the flow section between section 1-1 and section 2-2 at the initial time is taken as the research object. Taking into account the original flow of the liquid, after the dt period, the flow segment 1-2 moves to a new position 1′-2′, and the crest of the water hammer wave also propagates to the position of section 1′-1′.

As shown in Figure 2, after the dt period, section 1-1 of the studied control body moves to the position of section 1′-1. The movement distance is ds. The water hammer wave moves from section 2-2 to Section 1′-1′ at wave speed c. The movement distance is y, and section 2-2 moves to Section 2′-2′.

At the beginning of dt, the momentum of the liquid in the 1-2 flow section is:

$$\begin{array}{l}\left(\rho +\frac{\partial \rho }{\partial s}\frac{x}{2}\right)\left(A+\frac{\partial A}{\partial s}\frac{x}{2}\right)\left(v+\frac{\partial v}{\partial s}\frac{x}{2}\right)x=\rho Avx+\frac{1}{2}\frac{\partial \left(\rho Av\right)}{\partial s}{x}^2\\ +\frac{1}{4}\left(\rho \frac{\partial A}{\partial s}\frac{\partial v}{\partial s}+A\frac{\partial \rho }{\partial s}\frac{\partial v}{\partial s}+v\frac{\partial \rho }{\partial s}\frac{\partial A}{\partial s}\right)x^3+\frac{1}{8}\frac{\partial \rho }{\partial s}\frac{\partial A}{\partial s}\frac{\partial v}{\partial s}{x}^4\approx \rho Avx+\frac{1}{2}\frac{\partial \left(\rho Av\right)}{\partial s}{x}^2\end{array}$$

(1)

At the end of time dt, the average density of the fluid in the 1′-2′ flow section is:

$$\begin{array}{l}\tilde{\rho }=\frac{1}{2}\left\{\rho +\frac{\partial \rho }{\partial s}\mathrm{d}s+\frac{\partial }{\partial t}\left(\rho +\frac{\partial \rho }{\partial s}\mathrm{d}s\right)\mathrm{d}t+\rho +\frac{\partial \rho }{\partial s}\left(x+\mathrm{d}l\right)+\frac{\partial }{\partial t}\left[\rho +\frac{\partial \rho }{\partial s}\left(x+\mathrm{d}l\right)\right]\mathrm{d}t\right\}\\ \approx \rho +\frac{\partial \rho }{\partial s}\mathrm{d}s+\frac{\partial \rho }{\partial t}\mathrm{d}t+\frac{1}{2}\frac{\partial \rho }{\partial s}x\end{array}$$

(2)

In the same way, at the end of time dt, the average area of the fluid in the 1′-2′ flow section is:

$$\tilde{A}\approx A+\frac{\partial A}{\partial s}\mathrm{d}s+\frac{\partial A}{\partial t}\mathrm{d}t+\frac{1}{2}\frac{\partial A}{\partial s}x$$

(3)

The average velocity of the fluid in the 1′-2′ stream section is:

$$\tilde{v}\approx v+\frac{\partial v}{\partial s}\mathrm{d}s+\frac{\partial v}{\partial t}\mathrm{d}t+\frac{1}{2}\frac{\partial v}{\partial s}x$$

(4)

Therefore, at the end of time dt, the fluid momentum in the 1′-2′ flow section is:

$$\begin{array}{ll}mo& =\left(\rho +\frac{\partial \rho }{\partial s}\mathrm{d}s+\frac{\partial \rho }{\partial t}\mathrm{d}t+\frac{1}{2}\frac{\partial \rho }{\partial s}x\right)\left(A+\frac{\partial A}{\partial s}\mathrm{d}s+\frac{\partial A}{\partial t}\mathrm{d}t+\frac{1}{2}\frac{\partial A}{\partial s}x\right)\left(v+\frac{\partial v}{\partial s}\mathrm{d}s+\frac{\partial v}{\partial t}\mathrm{d}t+\frac{1}{2}\frac{\partial v}{\partial s}x\right)\left(y+\mathrm{d}l\right)\\ & =\rho Avx+\frac{\partial \left(\rho Av\right)}{\partial s}x\mathrm{d}s+\frac{\partial \left(\rho Av\right)}{\partial t}x\mathrm{d}t+\frac{1}{2}\frac{\partial \left(\rho Av\right)}{\partial s}{x}^2+\rho Avx\frac{\partial v}{\partial s}\mathrm{d}t\end{array}$$

(5)

At the time of dt, the increase of momentum in the studied control body is:

$$\begin{array}{ll}mo=& \left[\rho Avx+\frac{\partial \left(\rho Av\right)}{\partial s}x\mathrm{d}s+\frac{\partial \left(\rho Av\right)}{\partial t}x\mathrm{d}t+\frac{1}{2}\frac{\partial \left(\rho Av\right)}{\partial s}{x}^2+\rho Avx\frac{\partial v}{\partial s}\mathrm{d}t\right]-\\ & \left[\rho Avx+\frac{1}{2}\frac{\partial \left(\rho Av\right)}{\partial s}{x}^2\right]=\frac{\partial \left(\rho Av\right)}{\partial s}x\mathrm{d}s+\frac{\partial \left(\rho Av\right)}{\partial t}x\mathrm{d}t+\frac{1}{2}\frac{\partial \left(\rho Av\right)}{\partial s}{x}^2\end{array}$$

(6)

The external force acting on the control body includes the hydrodynamic pressure of the cross-section at both ends, the hydrodynamic pressure of the side, gravity and the frictional resistance of the pipe wall. Because the distance of the control body movement in the dt period is very small, when analyzing the combined external force it receives, the analysis of the 1-2 flow section is equivalent to the 1′-2′ analysis. Therefore, at the time of dt, the impulse of the external force received by the control body is:

$$\frac{\partial \left(\rho Av\right)}{\partial s}x\mathrm{d}s+\frac{\partial \left(\rho Av\right)}{\partial t}x\mathrm{d}t+\rho Avx\frac{\partial v}{\partial s}\mathrm{d}t=\left[-\frac{\partial \left(pA\right)}{\partial s}x+p\frac{\partial A}{\partial s}x+\rho gAx\sin \theta -{\tau }_0\chi x\right]\mathrm{d}t$$

(7)

After sorting, the formula can be obtained:

$$\left(\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}p}+\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}p}\right)\frac{\mathrm{d}p}{\mathrm{d}t}+\frac{1}{v}\frac{\mathrm{d}v}{\mathrm{d}t}+\frac{\partial v}{\partial s}+\frac{1}{\rho v}\frac{\partial p}{\partial s}-\frac{g}{v}\sin \theta +\frac{{\tau }_0\chi }{\rho Av}=0$$

(8)

Through the simultaneous establishment of a one-dimensional unsteady continuity equation:

$$\frac{1}{v}\frac{\partial v}{\partial t}+\frac{\partial v}{\partial s}+\frac{1}{\rho v}\frac{\partial p}{\partial s}-\frac{g}{v}\sin \theta +\frac{{\tau }_0\chi }{\rho Av}=0$$

(9)

Because of $\frac{\partial z}{\partial s}=-\sin \theta$, substitute the above formula to obtain:

$$\frac{1}{g}\frac{\partial v}{\partial t}+\frac{v}{g}\frac{\partial v}{\partial s}+\frac{1}{\gamma }\frac{\partial p}{\partial s}+\frac{\partial z}{\partial s}+\frac{{\tau }_0\chi }{\gamma A}=0$$

(10)

Equation (10) is the water hammer motion equation.

2.2. Fluid Continuity Equation

The derivation of a water hammer continuity equation uses the law of mass conservation and takes the fluid in the flow section with length x as the research object, as shown in Figure 2. Both ends of the pipe section are section 1-1 and section 2-2, respectively.

The difference in the quality of the liquid flowing into and out of the control body during the dt period is:

$$\frac{1}{g}\frac{\partial v}{\partial t}+\frac{v}{g}\frac{\partial v}{\partial s}+\frac{1}{\gamma }\frac{\partial p}{\partial s}+\frac{\partial z}{\partial s}+\frac{{\tau }_0\chi }{\gamma A}=0$$

(11)

At the beginning of the dt period, the fluid quality in flow section 1-2 is:

$$\left(\rho +\frac{\partial \rho }{\partial s}\frac{x}{2}\right)\left(A+\frac{\partial A}{\partial s}\frac{x}{2}\right)x\approx \rho Ax+\frac{1}{2}\frac{\partial \left(\rho A\right)}{\partial s}{x}^2$$

(12)

At the end of the dt period, the fluid quality in flow section 1-2 is:

$$\left(\rho +\frac{\partial \rho }{\partial s}\frac{x}{2}+\frac{\partial \rho }{\partial t}\mathrm{d}t\right)\left(A+\frac{\partial A}{\partial s}\frac{x}{2}+\frac{\partial A}{\partial t}\mathrm{d}t\right)x\approx \rho Ax+\frac{\partial \left(\rho A\right)}{\partial s}\frac{x^2}{2}+\frac{\partial \left(\rho A\right)}{\partial t}x\mathrm{d}t$$

(13)

Therefore, the fluid mass increase in flow segment 1-2 during the dt period is:

$$\left[\rho Ax+\frac{\partial \left(\rho A\right)}{\partial s}\frac{x^2}{2}+\frac{\partial \left(\rho A\right)}{\partial t}x\mathrm{d}t\right]-\left[\rho Ax+\frac{\partial \left(\rho A\right)}{\partial s}\frac{x^2}{2}\right]=\frac{\partial \left(\rho A\right)}{\partial t}x\mathrm{d}t$$

(14)

According to the law of conservation of mass, during the dt period, the difference between the mass of fluid flowing into and out of the control body should be equal to the increase in the mass of the fluid in the control body during the same period, that is:

$$\frac{\partial \left(\rho A\right)}{\partial t}+\frac{\partial \left(\rho Av\right)}{\partial s}=0$$

(15)

Equation (15) is the water hammer continuity equation derived in this paper. Therefore, the water hammer calculation equations can be obtained as:

$$\left\{\begin{array}{l}\frac{\partial \left(\rho A\right)}{\partial t}+\frac{\partial \left(\rho Av\right)}{\partial s}=0\\ \frac{1}{g}\frac{\partial v}{\partial t}+\frac{v}{g}\frac{\partial v}{\partial s}+\frac{1}{\gamma }\frac{\partial p}{\partial s}+\frac{\partial z}{\partial s}+\frac{{\tau }_0\chi }{\gamma A}=0\end{array}\right.$$

(16)

2.3. Complete Equations of Water Hammer

Generally in engineering, it is customary to use the piezometer head H to reflect the pressure p. Therefore, substituting $p=\rho g\left(H-z\right)$ into the equation of motion, we can get:

$$\frac{\mathrm{d}p}{\mathrm{d}t}=\frac{\partial p}{\partial t}+v\frac{\partial p}{\partial s}=\frac{\partial \left[\rho g\left(H-z\right)\right]}{\partial t}+v\frac{\partial \left[\rho g\left(H-z\right)\right]}{\partial s}=\rho g\frac{\mathrm{d}\left(H-z\right)}{\mathrm{d}t}+g\left(H-z\right)\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

(17)

where, $\frac{\partial z}{\partial t}=0$, $\frac{\partial z}{\partial s}=-\sin \theta$, θ is the included angle between the pipe axis and the horizontal direction.

The above formula considers the variation of liquid density ρ with time t and distance s during water hammer, which is different from the traditional water hammer equations.

Since:

$$\frac{\partial p}{\partial s}=\frac{\partial \left[\rho g\left(H-z\right)\right]}{\partial s}=\gamma \frac{\partial H}{\partial s}-\gamma \frac{\partial z}{\partial s}+g\left(H-z\right)\frac{\partial \rho }{\partial s}$$

(18)

Then substitute $\gamma =\rho g$, ${\tau }_0=\frac{1}{8}f\rho v^2$ and ${\chi }_0=\pi D=2\sqrt{\pi A}$ into the equation of motion to obtain:

$$\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial s}+g\frac{\partial H}{\partial s}+\frac{\left(H-z\right)g}{\rho }\frac{\partial \rho }{\partial s}+\frac{fv\left|v\right|}{4}\sqrt{\frac{\pi }{A}}=0$$

(19)

Since the fluid elasticity equation is:

$$\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}p}=\frac{1}{K}$$

(20)

The elastic equation of the pipe is:

$$\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}p}=\frac{D}{Ee}$$

(21)

The simultaneous hydroelasticity Equations (20) and (19) can be obtained:

$$\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial s}-\frac{\rho }{\left[\frac{K}{\rho g}-\left(H-z\right)\right]}\left(\frac{\partial H}{\partial t}+v\frac{\partial H}{\partial s}+v\sin \theta \right)=0$$

(22)

The simultaneous elastic equation of pipe (21) and Equation (19) can be obtained:

$$\frac{\partial A}{\partial t}+v\frac{\partial A}{\partial s}-\frac{\frac{2AK}{Ee}\sqrt{\frac{A}{\pi }}}{\left[\frac{K}{\rho g}-\left(H-z\right)\right]}\left(\frac{\partial H}{\partial t}+v\frac{\partial H}{\partial s}+v\sin \theta \right)=0$$

(23)

After sorting, the formula can be obtained:

$$\frac{\partial H}{\partial t}+v\frac{\partial H}{\partial s}+v\sin \theta +\frac{\frac{1}{\rho g}-\frac{\left(H-z\right)}{K}}{\frac{2}{Ee}\sqrt{\frac{A}{\pi }}+\frac{1}{K}}\frac{\partial v}{\partial s}=0$$

(24)

Thus, the water hammer closure equations are:

$$\left\{\begin{array}{l}\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial s}+g\frac{\partial H}{\partial s}+\frac{\left(H-z\right)g}{\rho }\frac{\partial \rho }{\partial s}+\frac{fv\left|v\right|}{4}\sqrt{\frac{\pi }{A}}=0\\ \frac{\partial H}{\partial t}+v\frac{\partial H}{\partial s}+v\sin \theta +\frac{\frac{1}{\rho g}-\frac{\left(H-z\right)}{K}}{\frac{2}{Ee}\sqrt{\frac{A}{\pi }}+\frac{1}{K}}\frac{\partial v}{\partial s}=0\\ \frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial s}-\frac{\rho }{\left[\frac{K}{\rho g}-\left(H-z\right)\right]}\left(\frac{\partial H}{\partial t}+v\frac{\partial H}{\partial s}+v\sin \theta \right)=0\\ \frac{\partial A}{\partial t}+v\frac{\partial A}{\partial s}-\frac{\frac{2AK}{Ee}\sqrt{\frac{A}{\pi }}}{\left[\frac{K}{\rho g}-\left(H-z\right)\right]}\left(\frac{\partial H}{\partial t}+v\frac{\partial H}{\partial s}+v\sin \theta \right)=0\end{array}\right.$$

(25)

3. Model Solving and Boundary Conditions

3.1. Difference Equation

In the calculation of water hammer equations, the more common methods are the characteristic line method and the direct difference method based on the finite difference method [27,28]. Because their deformation form cannot be transformed into characteristic equations, the characteristic line method cannot be used for calculation of the Gestalt equations proposed in this paper. Instead, the direct difference method is selected for calculation.

Compatibility, convergence and stability must be considered in difference calculations [28,29,30]. Change the differential equation into a difference equation. If the computational grid is reduced, the difference equation tends to be consistent with the differential equation. For example, the solution of the difference equation converges to the solution of the differential equation, which is convergence. If the rounding error and initial error in the calculation are always controlled within a limited range, and the calculated value is close to the true solution, then this scheme is stable. For the diffusion difference scheme adopted in this paper, the stability condition is the Courant condition. Since the water hammer equation is a hyperbolic equation, its stability condition should still meet the following requirements even though the concept of water hammer wave velocity is not introduced:

$$\frac{\Delta s}{\Delta t}\ge v \mathrm{OR} \frac{\Delta s}{\Delta t}\ge {\left|v\pm \frac{1}{\sqrt{\rho \left(\frac{1}{K}+\frac{D}{Ee}\right)}}\right|}_{\max }$$

(26)

In order to solve this model, the basic equation needs to be discretized and the differential equation needs to be transformed into a differential algebraic equation. Here, the diffusion difference scheme is used to solve the equation. Figure 3 shows the L-t plane of characteristic grid of the four variable water hammer model. Use X to represent the objective function head H, velocity V and density ρ and area A. Then the partial differential of X to time and position can be expressed as:

$$\frac{\partial X}{\partial t}=\frac{X_i^k-\left[\xi X_i^{k-1}+\frac{1-\xi }{2}\left(X_{i+1}^{k-1}+X_{i-1}^{k-1}\right)\right]}{\Delta t}$$

(27)

$$\frac{\partial X}{\partial z}=\frac{X_{i+1}^{k-1}-X_{i-1}^{k-1}}{2\Delta z}$$

(28)

By introducing Equations (27) and (28) into the water hammer closure calculation model, the following difference equations can be obtained:

$$\begin{array}{l}\frac{v_i^k-\left[\xi v_{i }^{k-1}+\frac{1-\xi }{2}\left(v_{i+1 }^{k-1}+v_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{v_{i+1 }^{k-1}-v_{i-1 }^{k-1}}{2\Delta z}\\ +g\frac{H_{i+1 }^{k-1}-H_{i-1 }^{k-1}}{2\Delta z}+\frac{H_{i }^{k-1}g}{{\rho }_{i }^{k-1}}\frac{{\rho }_{i+1 }^{k-1}-{\rho }_{i-1 }^{k-1}}{2\Delta z}+\frac{f{v}_{i }^{k-1}\left|v_{i }^{k-1}\right|}{2D}=0\end{array}$$

(29)

$$\begin{array}{l}\frac{H_i^k-\left[\xi H_{i }^{k-1}+\frac{1-\xi }{2}\left(H_{i+1 }^{k-1}+H_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{H_{i+1 }^{k-1}-H_{i-1 }^{k-1}}{2\Delta z}\\ +v_{i }^{k-1}\sin \theta +\frac{\frac{1}{\rho g}-\frac{\left(H-z\right)}{K}}{\frac{2}{Ee}\sqrt{\frac{A}{\pi }}+\frac{1}{K}}\frac{v_{i+1 }^{k-1}-v_{i-1 }^{k-1}}{2\Delta z}=0\end{array}$$

(30)

$$\begin{array}{l}\frac{{\rho }_i^k-\left[\xi {\rho }_{i }^{k-1}+\frac{1-\xi }{2}\left({\rho }_{i+1 }^{k-1}+{\rho }_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{{\rho }_{i+1 }^{k-1}-{\rho }_{i-1 }^{k-1}}{2\Delta z}-\frac{{\rho }_{i }^{k-1}}{\frac{K}{{\rho }_{i }^{k-1}g}-\left(H_{i }^{k-1}-z_i\right)}\cdot \\ \left\{\frac{H_i^k-\left[\xi H_{i }^{k-1}+\frac{1-\xi }{2}\left(H_{i+1 }^{k-1}+H_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\left(\frac{H_{i+1 }^{k-1}-H_{i-1 }^{k-1}}{2\Delta z}+\sin \theta \right)\right\}=0\end{array}$$

(31)

$$\begin{array}{l}\frac{A_i^k-\left[\xi A_{i }^{k-1}+\frac{1-\xi }{2}\left(A_{i+1 }^{k-1}+A_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{A_{i+1 }^{k-1}-A_{i-1 }^{k-1}}{2\Delta z}-\frac{2K{A}_{i }^{k-1}}{{\rho }_{i }^{k-1}Ee}\sqrt{\frac{A_{i }^{k-1}}{\pi }}\cdot \\ \left\{\frac{{\rho }_i^k-\left[\xi {\rho }_{i }^{k-1}+\frac{1-\xi }{2}\left({\rho }_{i+1 }^{k-1}+{\rho }_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{{\rho }_{i+1 }^{k-1}-{\rho }_{i-1 }^{k-1}}{2\Delta z}\right\}=0\end{array}$$

(32)

Make the following assumptions:

$$\left\{\begin{array}{l}{X}_{H1}=\xi H_i^{k-1}+\frac{1-\xi }{2}\left(H_{i+1}^{k-1}+H_{i-1}^{k-1}\right), X_{H2}=\left(H_{i+1}^{k-1}-H_{i-1}^{k-1}\right)\frac{\Delta t}{\Delta z}\\ X_{v1}=\xi v_i^{k-1}+\frac{1-\xi }{2}\left(v_{i+1}^{k-1}+v_{i-1}^{k-1}\right), X_{v2}=\left(v_{i+1}^{k-1}-v_{i-1}^{k-1}\right)\frac{\Delta t}{\Delta z}\\ X_{\rho 1}=\xi {\rho }_i^{k-1}+\frac{1-\xi }{2}\left({\rho }_{i+1}^{k-1}+{\rho }_{i-1}^{k-1}\right), X_{\rho 2}=\left({\rho }_{i+1}^{k-1}-{\rho }_{i-1}^{k-1}\right)\frac{\Delta t}{\Delta z}\\ X_{A1}=\xi A_i^{k-1}+\frac{1-\xi }{2}\left(A_{i+1}^{k-1}+A_{i-1}^{k-1}\right), X_{A2}=\left(A_{i+1}^{k-1}-A_{i-1}^{k-1}\right)\frac{\Delta t}{\Delta z}\end{array}\right.$$

(33)

After sorting, the expressions of the four unknowns are as follows:

$$\left\{\begin{array}{l}{H}_i^k=X_{H1}-v_i^{k-1}\left(\frac{X_{H2}}{2}+\Delta tsin\theta \right)-X_{v2}\left[\frac{\frac{1}{{\rho }_i^{k-1}g}-\frac{\left(H_i^{k-1}-z_i\right)}{K}}{\frac{4}{Ee}\sqrt{\frac{A_i^{k-1}}{\pi }}+\frac{1}{K}}\right]\\ v_i^k=X_{v1}-\frac{v_i^{k-1}}{2}{X}_{v2}-\frac{g}{2}{X}_{H2}-\frac{\left(H_i^{k-1}-z_i\right)g}{2{\rho }_i^{k-1}}{X}_{\rho 2}-\frac{f\Delta t}{2D}{v}_i^{k-1}\left|v_i^{k-1}\right|=0\\ {\rho }_i^k=X_{\rho 1}-\frac{v_i^{k-1}}{2}{X}_{\rho 2}+\frac{{\rho }_i^{k-1}}{\left[\frac{K}{{\rho }_i^{k-1}g}-\left(H_i^{k-1}-z_i\right)\right]}\left[H_i^k-X_{H1}+v_i^{k-1}\left(\frac{X_{H2}}{2}+\Delta tsin\theta \right)\right]\\ A_i^k=X_{A1}-\frac{v_i^{k-1}}{2}{X}_{A2}+\frac{2A_i^{k-1}K}{{\rho }_i^{k-1}Ee}\sqrt{\frac{A_i^{k-1}}{\pi }}\left({\rho }_i^k-X_{\rho 1}+\frac{v_i^{k-1}}{2}{X}_{\rho 2}\right)\end{array}\right.$$

(34)

The fluid in the oil pipeline is not all liquid in the field. What happens when there is a small amount of gas in the oil pipeline? Assuming that there are only gas-liquid two phases in the pipe, the density of the gas-liquid mixture is [31]:

$${\rho }_{\mathrm{m}}={\rho }_{\mathrm{g}}{C}_{\mathrm{g}}+{\rho }_{\mathrm{l}}\left(1-C_{\mathrm{g}}\right)$$

(35)

where, C_g is the gas content.

Since the density in this paper is considered a variable, according to the difference principle, we can get:

$${\rho }_{\mathrm{m}}{}_i^k={\rho }_{\mathrm{g}}{}_i^k{C}_{\mathrm{g}}{}_i^k+{\rho }_{\mathrm{l}}{}_i^k\left(1-C_{\mathrm{g}}{}_i^k\right)$$

(36)

$${\rho }_{\mathrm{m}}{}_i^{k-1}={\rho }_{\mathrm{g}}{}_i^{k-1}{C}_{\mathrm{g}}{}_i^{k-1}+{\rho }_{\mathrm{l}}{}_i^{k-1}\left(1-C_{\mathrm{g}}{}_i^{k-1}\right)$$

(37)

$${\rho }_{\mathrm{m}}{}_{i+1}^{k-1}={\rho }_{\mathrm{g}}{}_{i+1}^{k-1}{C}_{\mathrm{g}}{}_{i+1}^{k-1}+{\rho }_{\mathrm{l}}{}_{i+1}^{k-1}\left(1-C_{\mathrm{g}}{}_{i+1}^{k-1}\right)$$

(38)

$${\rho }_{\mathrm{m}}{}_{i-1}^{k-1}={\rho }_{\mathrm{g}}{}_{i-1}^{k-1}{C}_{\mathrm{g}}{}_{i-1}^{k-1}+{\rho }_{\mathrm{l}}{}_{i-1}^{k-1}\left(1-C_{\mathrm{g}}{}_{i-1}^{k-1}\right)$$

(39)

Accordingly, the flow rate also becomes the mixing flow rate, and its expression is:

$$v_{\mathrm{m}}=\frac{Q_{\mathrm{g}}+Q_{\mathrm{l}}}{A}$$

(40)

Then according to the difference principle, we can get:

$$v_{\mathrm{m}}{}_i^k=\frac{Q_{\mathrm{g}}{}_i^k+Q_{\mathrm{l}}{}_i^k}{A_i^k}$$

(41)

$$v_{\mathrm{m}}{}_{i }^{k-1}=\frac{Q_{\mathrm{g}}{}_{i }^{k-1}+Q_{\mathrm{l}}{}_{i }^{k-1}}{A_{i }^{k-1}}$$

(42)

$$v_{\mathrm{m}}{}_{i+1 }^{k-1}=\frac{Q_{\mathrm{g}}{}_{i+1 }^{k-1}+Q_{\mathrm{l}}{}_{i+1 }^{k-1}}{A_{i+1 }^{k-1}}$$

(43)

$$v_{\mathrm{m}}{}_{i-1 }^{k-1}=\frac{Q_{\mathrm{g}}{}_{i-1 }^{k-1}+Q_{\mathrm{l}}{}_{i-1 }^{k-1}}{A_{i-1 }^{k-1}}$$

(44)

Thus, formulas (29)–(32) are deformed as:

$$\begin{array}{l}\frac{v_{\mathrm{m}}{}_i^k-\left[\xi v_{\mathrm{m}}{}_{i }^{k-1}+\frac{1-\xi }{2}\left(v_{\mathrm{m}}{}_{i+1 }^{k-1}+v_{\mathrm{m}}{}_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{\mathrm{m}}{}_{i }^{k-1}\frac{v_{\mathrm{m}}{}_{i+1 }^{k-1}-v_{\mathrm{m}}{}_{i-1 }^{k-1}}{2\Delta z}\\ +g\frac{H_{i+1 }^{k-1}-H_{i-1 }^{k-1}}{2\Delta z}+\frac{H_{i }^{k-1}g}{{\rho }_{\mathrm{m}}{}_{i }^{k-1}}\frac{{\rho }_{\mathrm{m}}{}_{i+1 }^{k-1}-{\rho }_{\mathrm{m}}{}_{i-1 }^{k-1}}{2\Delta z}+\frac{f{v}_{\mathrm{m}}{}_{i }^{k-1}\left|v_{\mathrm{m}}{}_{i }^{k-1}\right|}{2D}=0\end{array}$$

(45)

$$\begin{array}{l}\frac{H_i^k-\left[\xi H_{i }^{k-1}+\frac{1-\xi }{2}\left(H_{i+1 }^{k-1}+H_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{\mathrm{m}}{}_{i }^{k-1}\frac{H_{i+1 }^{k-1}-H_{i-1 }^{k-1}}{2\Delta z}\\ +v_{\mathrm{m}}{}_{i }^{k-1}\sin \theta +\frac{\frac{1}{{\rho }_{\mathrm{m}}{}_{i }^{k-1}g}-\frac{\left(H_{i }^{k-1}-z\right)}{K}}{\frac{2}{Ee}\sqrt{\frac{A_{i }^{k-1}}{\pi }}+\frac{1}{K}}\frac{v_{\mathrm{m}}{}_{i+1 }^{k-1}-v_{\mathrm{m}}{}_{i-1 }^{k-1}}{2\Delta z}=0\end{array}$$

(46)

$$\begin{array}{l}\frac{{\rho }_{\mathrm{m}}{}_i^k-\left[\xi {\rho }_{\mathrm{m}}{}_{i }^{k-1}+\frac{1-\xi }{2}\left({\rho }_{\mathrm{m}}{}_{i+1 }^{k-1}+{\rho }_{\mathrm{m}}{}_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{{\rho }_{\mathrm{m}}{}_{i+1 }^{k-1}-{\rho }_{\mathrm{m}}{}_{i-1 }^{k-1}}{2\Delta z}-\frac{{\rho }_{\mathrm{m}}{}_{i }^{k-1}}{\frac{K}{{\rho }_{i }^{k-1}g}-\left(H_{i }^{k-1}-z_i\right)}\cdot \\ \left\{\frac{H_i^k-\left[\xi H_{i }^{k-1}+\frac{1-\xi }{2}\left(H_{i+1 }^{k-1}+H_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\left(\frac{H_{i+1 }^{k-1}-H_{i-1 }^{k-1}}{2\Delta z}+\sin \theta \right)\right\}=0\end{array}$$

(47)

$$\begin{array}{l}\frac{A_i^k-\left[\xi A_{i }^{k-1}+\frac{1-\xi }{2}\left(A_{i+1 }^{k-1}+A_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{A_{i+1 }^{k-1}-A_{i-1 }^{k-1}}{2\Delta z}-\frac{2K{A}_{i }^{k-1}}{{\rho }_{i }^{k-1}Ee}\sqrt{\frac{A_{i }^{k-1}}{\pi }}\cdot \\ \left\{\frac{{\rho }_{\mathrm{m}}{}_i^k-\left[\xi {\rho }_{\mathrm{m}}{}_{i }^{k-1}+\frac{1-\xi }{2}\left({\rho }_{\mathrm{m}}{}_{i+1 }^{k-1}+{\rho }_{\mathrm{m}}{}_{i-1 }^{k-1}\right)\right]}{\Delta t}+v_{i }^{k-1}\frac{{\rho }_{\mathrm{m}}{}_{i+1 }^{k-1}-{\rho }_{\mathrm{m}}{}_{i-1 }^{k-1}}{2\Delta z}\right\}=0\end{array}$$

(48)

The corresponding assumptions will also be changed accordingly:

$$\left\{\begin{array}{l}{X}_{v1}^{\mathrm{m}}=\xi v_{\mathrm{m}}{}_{i }^{k-1}+\frac{1-\xi }{2}\left(v_{\mathrm{m}}{}_{i+1 }^{k-1}+v_{\mathrm{m}}{}_{i-1 }^{k-1}\right),X_{v2}^{\mathrm{m}}=\left(v_{\mathrm{m}}{}_{i+1 }^{k-1}-v_{\mathrm{m}}{}_{i-1 }^{k-1}\right)\frac{\Delta t}{\Delta z}\\ X_{\rho 1}^{\mathrm{m}}=\xi {\rho }_{\mathrm{m}}{}_{i }^{k-1}+\frac{1-\xi }{2}\left({\rho }_{\mathrm{m}}{}_{i+1 }^{k-1}+{\rho }_{\mathrm{m}}{}_{i-1 }^{k-1}\right),X_{\rho 2}^{\mathrm{m}}=\left({\rho }_{\mathrm{m}}{}_{i+1 }^{k-1}-{\rho }_{\mathrm{m}}{}_{i-1 }^{k-1}\right)\frac{\Delta t}{\Delta z}\end{array}\right.$$

(49)

After sorting, the expressions of four unknowns considering gas-liquid two-phase are obtained as follows:

$$\left\{\begin{array}{l}{H}_i^k=X_{H1}-{v_{\mathrm{m}}}_i^{k-1}\left(\frac{X_{H2}}{2}+\Delta tsin\theta \right)-X_{v2}^{\mathrm{m}}\left[\frac{\frac{1}{{{\rho }_{\mathrm{m}}}_i^{k-1}g}-\frac{\left(H_i^{k-1}-z_i\right)}{K}}{\frac{4}{Ee}\sqrt{\frac{A_i^{k-1}}{\pi }}+\frac{1}{K}}\right]\\ {v_{\mathrm{m}}}_i^k=X_{v1}^{\mathrm{m}}-\frac{{v_{\mathrm{m}}}_i^{k-1}}{2}{X}_{v2}^{\mathrm{m}}-\frac{g}{2}{X}_{H2}-\frac{\left(H_i^{k-1}-z_i\right)g}{2{{\rho }_{\mathrm{m}}}_i^{k-1}}{X}_{\rho 2}^{\mathrm{m}}-\frac{f\Delta t}{2D}{v_{\mathrm{m}}}_i^{k-1}\left|{v_{\mathrm{m}}}_i^{k-1}\right|=0\\ {{\rho }_{\mathrm{m}}}_i^k=X_{\rho 1}^{\mathrm{m}}-\frac{{v_{\mathrm{m}}}_i^{k-1}}{2}{X}_{\rho 2}^{\mathrm{m}}+\frac{{{\rho }_{\mathrm{m}}}_i^{k-1}}{\left[\frac{K}{{{\rho }_{\mathrm{m}}}_i^{k-1}g}-\left(H_i^{k-1}-z_i\right)\right]}\left[H_i^k-X_{H1}+{v_{\mathrm{m}}}_i^{k-1}\left(\frac{X_{H2}}{2}+\Delta tsin\theta \right)\right]\\ A_i^k=X_{A1}-\frac{{v_{\mathrm{m}}}_i^{k-1}}{2}{X}_{A2}+\frac{2A_i^{k-1}K}{{{\rho }_{\mathrm{m}}}_i^{k-1}Ee}\sqrt{\frac{A_i^{k-1}}{\pi }}\left({{\rho }_{\mathrm{m}}}_i^k-X_{\rho 1}^{\mathrm{m}}+\frac{{v_{\mathrm{m}}}_i^{k-1}}{2}{X}_{\rho 2}^{\mathrm{m}}\right)\end{array}\right.$$

(50)

3.2. Boundary Conditions

The boundary point is different from the inner node. It has an equation with boundary conditions, and the number of equations is more than the number of unknowns. For the previous solutions, one is to synthesize the continuity equation and motion equation into a partial differential equation and then transform it into a difference equation. The second is not to establish the continuity equation and the motion equation, but to introduce the characteristic equation. For example, the inverse characteristic line equation can be introduced for the upstream and the along characteristic line equation can be introduced for the downstream [32,33]. For the equation form of four variables in this paper, neither method is applicable. The solution obtained by method one does not converge and cannot meet the stability condition. However, if the continuity Equation (24) is abandoned, the ideal results can be obtained by using the equation of motion (19) and Equations (20) and (21) together with the upstream boundary conditions. The downstream boundary of a simple pipeline can be treated by method one, and the results are stable and convergent.

The equations of motion (19)–(21) are used to solve the upstream boundary conditions. When the upstream is a reservoir, the reservoir water level can be regarded as constant, and the boundary condition is a constant, that is $\frac{\partial H}{\partial t}=0$, from the motion equation:

$$v_0^k=v_{0 }^{k-1}\left(1-{X^{\prime }}_{v2}\right)-g\left[{X^{\prime }}_{H2}+\frac{\left(H_{0 }^{k-1}-z_0\right)}{{\rho }_{0 }^{k-1}}{X^{\prime }}_{\rho 2}-\frac{f\Delta t}{2gD}{v}_{0 }^{k-1}\left|v_{0 }^{k-1}\right|\right]$$

(51)

$${\rho }_0^k={\rho }_{0 }^{k-1}-v_{0 }^{k-1}{X^{\prime }}_{\rho 2}+\frac{{\rho }_{0 }^{k-1}}{\frac{K}{{\rho }_{0 }^{k-1}g}-\left(H_{0 }^{k-1}-z_0\right)}{v}_{0 }^{k-1}\left({X^{\prime }}_{A2}+\Delta t\sin \theta \right)$$

(52)

Combine formulas (20) and (21) and substitute the calculated value to obtain:

$$A_0^k=A_0^{k-1}-v_{0 }^{k-1}{X^{\prime }}_{A2}+\frac{2A_{0 }^{k-1}K}{{\rho }_{0 }^{k-1}Ee}\sqrt{\frac{A_{0 }^{k-1}}{\pi }}\left({\rho }_0^k-{\rho }_0^{k-1}+v_{0 }^{k-1}{X^{\prime }}_{\rho 2}\right)$$

(53)

where in:

$$\left\{\begin{array}{l}{X^{\prime }}_{H2}=\left(H_{1 }^{k-1}-H_{0 }^{k-1}\right)\frac{\Delta t}{\Delta z} {X^{\prime }}_{v2}=\left(v_{1 }^{k-1}-v_{0 }^{k-1}\right)\frac{\Delta t}{\Delta z}\\ {X^{\prime }}_{\rho 2}=\left({\rho }_{1 }^{k-1}-{\rho }_{0 }^{k-1}\right)\frac{\Delta t}{\Delta z} {X^{\prime }}_{A2}=\left(A_{1 }^{k-1}-A_{0 }^{k-1}\right)\frac{\Delta t}{\Delta z}\end{array}\right.$$

(54)

Here, the equation of motion (19) and the continuity Equation (24) are integrated into a partial differential equation and then transformed into a difference equation. The specific form is:

$$\frac{\partial H}{\partial t}+\left(v+v_0\right)\frac{\partial H}{\partial s}+\frac{v_0}{g}\frac{\partial v}{\partial t}+\left[\frac{\frac{1}{\rho g}-\frac{\left(H-z\right)}{K}}{\frac{2}{Ee}\sqrt{\frac{A}{\pi }}+\frac{1}{K}}+\frac{v_{0}v}{g}\right]\frac{\partial v}{\partial s}+\frac{\left(H-z\right)v_0}{\rho }\frac{\partial \rho }{\partial s}+v\sin \theta +v_0\frac{fv\left|v\right|}{2gD}=0$$

(55)

When the downstream end of the pipeline is a valve, it can be regarded as an orifice to obtain [34]:

$$v_t^A=v_{\mathrm{m}}{\tau }_k\sqrt{\frac{H_t^A}{H_0}}$$

(56)

where H₀ is the initial head pressure of the downstream section. For downstream N section, rewrite it as:

$$v_N^k=v_{\mathrm{m}}{\tau }_k\sqrt{\frac{H_N^k}{H_N^0}}$$

(57)

Solve the above equation with Equations (20), (21) and (39):

$$v_N^k=-Y_1{Y}_2+\sqrt{{\left(Y_1{Y}_2\right)}^2-2Y_2{Y}_3}$$

(58)

When the flow rate is obtained, the following formula (37) can be obtained:

$$H_N^k=-Y_1{v}_N^k-Y_3$$

(59)

Substitute into Equations (20) and (21) to obtain:

$${\rho }_N^k={\rho }_{N }^{k-1}-v_{N }^{k-1}{X^{″}}_{\rho 2}+\frac{{\rho }_{N }^{k-1}}{\frac{K}{{\rho }_{N }^{k-1}g}-\left(H_{N }^{k-1}-z_0\right)}\left[H_{N }^k-H_{N }^{k-1}+v_{N }^{k-1}\left({X^{″}}_{A2}+\Delta t\sin \theta \right)\right]$$

(60)

$$A_N^k=A_N^{k-1}-v_N^{k-1}{X^{″}}_{A2}+\frac{2A_N^{k-1}K}{{\rho }_N^{k-1}Ee}\sqrt{\frac{A_N^{k-1}}{\pi }}\left({\rho }_N^k-{\rho }_N^{k-1}+v_N^{k-1}{X^{″}}_{\rho 2}\right)$$

(61)

where v_m is the initial flow velocity at the downstream section of the pipeline and τ_k is the opening of the downstream valve.

$$Y_1=\frac{v_0}{g}$$

(62)

$$Y_2=\frac{{\left(v_{\mathrm{m}}{\tau }_k\right)}^2}{2H_N^0}$$

(63)

$$\begin{array}{ll}{Y}_3=& -H_N^{k-1}+\left(v_N^{k-1}+v_0\right){X^{″}}_{H2}-v_0\left[\frac{v_N^{k-1}}{g}-\frac{\left(H{v}_N^{k-1}-z_0\right)}{\rho v_N^{k-1}}{X^{″}}_{\rho 2}-\frac{f\Delta t}{2gD}{v}_N^{k-1}\left|v_N^{k-1}\right|\right]\\ & +\left[\frac{\frac{1}{{\rho }_N^{k-1}g}-\frac{\left(H_N^{k-1}-z_N\right)}{K}}{\frac{2}{Ee}\sqrt{\frac{A_N^{k-1}}{\pi }}+\frac{1}{K}}+\frac{v_0{v}_N^{k-1}}{g}\right]{X^{″}}_{v2}+v_N^{k-1}\cdot \Delta t\sin \theta \end{array}$$

(64)

$$\left\{\begin{array}{l}{X^{″}}_{H2}=\left(H_{N }^{k-1}-H_{N-1 }^{k-1}\right)\frac{\Delta t}{\Delta z} {X^{″}}_{v2}=\left(v_{N }^{k-1}-v_{N-1 }^{k-1}\right)\frac{\Delta t}{\Delta z}\\ {X^{″}}_{\rho 2}=\left({\rho }_{N }^{k-1}-{\rho }_{N-1 }^{k-1}\right)\frac{\Delta t}{\Delta z} {X^{″}}_{A2}=\left(A_{N }^{k-1}-A_{N-1 }^{k-1}\right)\frac{\Delta t}{\Delta z}\end{array}\right.$$

(65)

The above is the treatment of the boundary problem of a simple pipeline. Figure 4 shows the calculation flowchart.

3.3. Model Validation

To verify the accuracy of the water hammer model, a water hammer calculation was carried out on a case of a fluid pipeline in the literature, and the calculation results were compared with those in the literature [35].

As shown in Figure 5, the calculation results of water hammer are compared. The red is the calculation result of the literature, and the blue is the calculation result of the article model. It can be seen from the figure that there are certain deviations in the results, which may be caused by different mathematical models since the model in this paper considers the four variables of pressure, flow rate, density and flow area.

4. Field Application and Result Discussion

4.1. Sample Pipeline Foundation Parameters

Taking an oil pipeline X as an example, this paper simulates the pressure fluctuation in the pipeline after the valve at the end is closed and the effects of different factors such as valve closing time, flow and gas content on water hammer pressure. The buried depth of the pipeline is 1.8 m, the pipe diameter is 300 mm, and the wall thickness is 6 mm. The physical parameters of oil, soil, anticorrosion coating and pipeline are shown in

Table 1. Physical parameters of oil, soil, anticorrosion coating and pipe.

Parameter	Value	Parameter	Value
Density of oil	841 kg/cm3	Specific heat of soil	2225 J/(K·kg)
Thermal diffusivity of oil	5.91 mm2/s	Specific heat of pipe	465 J/(K·kg)
Thermal conductivity of oil	0.162 W/(m·K)	Specific heat of anticorrosive coating	2000 J/(K·kg)
Coefficient of thermal expansion of oil	8.3 × 10−4 K−1	Specific heat of oil	1900 J/(K·kg)
Density of soil	1455 kg/cm3	Thermal conductivity of soil	1.4 W/(m·K)
Density of pipe	7850 kg/cm3	Thermal conductivity of pipe	48 W/(m·K)
Density of anticorrosive coating	1200 kg/cm3	Thermal conductivity of anticorrosive coating	0.15 W/(m·K)

. The atmospheric temperature is 20 °C, the soil temperature is 5 °C, and the oil temperature is 22 °C.

4.2. Result and Discussion

Assuming that the valve closing time is about 5 s, the valve closing index m is 1.0 (linear valve closing), and the calculation time is 200 s. Figure 6 is the change of pressure at the pipe end after valve closing. It can be seen from the figure that at the moment of valve closing, the fluid flow rate at the pipe end suddenly changes to 0, the pressure increases rapidly, reaches the peak value of 4.95 MPa at 5.08 s, then starts to decline, reaching the valley value of 3.37 MPa at 12.80 s. Due to friction along the way, the fluctuating pressure will decay quickly and finally reach the equilibrium pressure.

Take the section close to the pipe end as the research plane and calculate the fluctuation of velocity in this plane as shown in Figure 7. Before closing the valve, the initial flow velocity of the section is 2.85 m/s. When the valve is closed, the flow rate drops rapidly and fluctuates around 0. According to the principle of the Bernoulli equation, the fluctuation will become weaker and weaker over time. Finally, the fluid pressure in the pipe will be a certain value and the flow rate will be 0.

The effects of different valve closing time on the fluctuation of pressure and velocity at the pipe end are studied by the model. As shown in Figure 8a, the simulation results of fluctuating pressure at the pipe end under five different valve closing times are shown. According to the calculation, with the gradual increase of valve closing time, the maximum fluctuating pressure at the pipe end decreases, and the time of peak value lags behind. Figure 8b shows the simulation results of fluctuation velocity near the pipe end under five different valve closing times. The shorter the valve closing time, the earlier the flow rate drops to 0 and fluctuates around 0, and the greater the fluctuation range.

Figure 9 shows the influence of different flow rates on pipe end pressure. It can be seen from Figure 9a that the pressure at the pipe end first decreases and then increases at the moment of valve closing. When the flow in the pipe is 0.4 kg/s, the decline amplitude is the largest, but the rise speed is also the fastest and the peak value is the highest. Similarly, when the flow rate is 0.4 kg/s, the decline rate of pipe end pressure from peak to valley is the fastest and the valley is the smallest. That is, the greater the flow, the greater the change of pipe end pressure, the faster the speed change, and the more obvious the water hammer effect.

As Figure 9b shows, the greater the flow rate, the greater the initial flow velocity in the pipe and the fastest the decline rate of flow velocity. According to the enlarged view shown in Figure 9b, when the flow rate is large (e.g., flow rate is 0.4 kg/s), flow rate fluctuation is more intense. Moreover, there will be small fluctuations at inflection points such as falling from peak value and rising from valley value. The larger the flow, the more obvious the small fluctuations. It shows that the flow with large flow is greatly disturbed by instantaneous obstacles such as valve closing. However, with the increase of time, the impact of flow will be less and less obvious, whether it is pressure or flow rate.

After the fluctuating pressure is generated at the pipe end, it will propagate upstream along the pipeline to form a pressurization wave surface. To intuitively obtain the propagation of pressure wave in all calculation time, the three-dimensional diagram of pressure-wave propagation with pipe length and time after valve closing is drawn (Figure 10).

The effects of different pipe length and time on the fluctuating pressure in the pipe are studied, as shown in Figure 11. Figure 11a shows the pressure change chart in a 500 m section at 40 s. Pressure changes with pipe length when the blue line is 40 s. The red line shows the change of pressure with time at 500 m. By analogy, Figure 11b–d are the pressure change charts of a 1000 m section at 80 s, 1500 m section at 120 s and 2000 m section at 160 s after valve closing. The meaning of the blue line and the red line is the same as that in Figure 11a. It can be seen from the figure that the pressure fluctuation will gradually attenuate along the pipe length with the increase of time. The place with the greatest water hammer effect is near the valve. Under the coupling effect of time and tube length, the shorter the time and the shorter the tube length, the more obvious the pressure fluctuation; the longer the time and the longer the tube, the weaker the pressure fluctuation.

Through the water hammer model established in this paper, the changes of pipe end pressure and flow velocity under five conditions with gas content of 1%, 3%, 5%, 7% and 9% were calculated, respectively (Figure 12). Figure 12a shows that the greater the gas content, the smaller the fluctuation peak of pipe end pressure. The time of the first peak of pressure wave is basically the same, but the time of the first valley value begins to be different. The greater the gas content, the later the minimum pressure appears. With the increase of time, this difference becomes more and more obvious, which is reflected in a longer water hammer cycle. Since the water hammer period is the ratio of pipe length to pressure-wave velocity, the greater the void fraction, the smaller the pressure-wave velocity. Slightly different is that there is a difference in the time when the first fluctuation value of velocity wave appears. The greater the void fraction, the longer the time lag of flow rate fluctuation, as shown in Figure 12b. According to Bernoulli’s theory, pressure and velocity should be in a trade-off relationship and as the former increases, the latter decreases. Therefore, when the gas content is large, although the velocity fluctuation lags, the absolute value of velocity is large.

When there is free gas in the oil pipe, the gas will form uneven bubbles in the pipe. When the bubble is free in the oil, it can be approximately regarded as an elastic deformable body. When there is pressure fluctuation in the oil, the bubbles will be forced to compress and increase their internal energy. At this time, the bubbles are not stable, and work will be done to the surrounding oil to release energy, so that the momentum of the oil increases and compresses the surrounding bubbles. Due to energy dissipation, the increase of oil momentum is not equal to the increase of internal energy after bubble forcing, and the energy of bubble forcing comes from pressure fluctuation. Therefore, the actual pressure fluctuation value is obviously lower than that without gas, and the size of the pressure wave mainly depends on the gas content [36,37,38].

5. Conclusions

Based on the essence of water hammer, the interaction between pipeline fluid motion and water hammer wave propagation is considered, and the pressure, velocity, density and overflow area are set as variables. A new set of water hammer calculation equations is deduced and solved numerically. The effects of different valve closing time, flow rate and gas content on pressure distribution and water hammer effect were studied. The following conclusions can be drawn:

With the increase of valve closing time, the maximum fluctuating pressure at the pipe end decreases, and the time of peak value also lags behind. The shorter the valve closing time, the earlier the flow rate drops to zero and fluctuates around zero, and the greater the fluctuation range. When the valve closing time increases from 5 s to 25 s, the difference of water hammer pressure is 0.72 MPa and the difference of velocity fluctuation amplitude is 0.076 m/s.

With greater flow, the pressure change at the pipe end is greater. With faster speed changes, the water hammer effect becomes more obvious. With larger flow, flow velocity fluctuations increase, and there will be a small amplitude fluctuation. Our calculations show that the flow with large flow is greatly disturbed by instantaneous obstacles such as valve closing. However, with the increase of time, the impact of flow will be less and less obvious, whether it is pressure or flow rate.

With the increase of time, the pressure fluctuation will gradually attenuate along the pipe length. The place with the greatest water hammer effect is near the valve. Under the coupling effect of time and tube length, the shorter the time and the shorter the tube length, the more obvious the pressure fluctuation.

Our calculations show that the larger the gas content is, the smaller the fluctuation peak of pipe end pressure is, the longer the water hammer cycle is, and the smaller the pressure-wave velocity is. The actual pressure fluctuation value is obviously lower than that without gas, and the size of the pressure wave mainly depends on the gas content. When the gas content increases from 1% to 9%, the difference in water hammer pressure is 0.41 MPa. When the gas content is large, although the velocity fluctuation lags, the absolute value of velocity is large.