A NUMERICAL STUDY OF EFFECTS OF THE AXIAL STRESS ON UNSTEADY LIQUID PIPELINE FLOWS

We study the effects of the axial component of the shear stress on unsteady pipeline flows. We show that the axial component of the shear stress should be introduced in the modeling of unsteady flows, and as a numerical model, we propose a one-dimensional momentum equation in which a term containing the second derivative of the velocity with respect to space is introduced. The momentum equation and the continuity equation are converted into a system suitable for the application of upstream difference approximations. Numerical results are presented, and their correspondence with experimental results is examined to see how our model captures phenomena observed experimentally.


Introduction.
A small disturbance in a liquid pipeline flow may result in an unexpectedly rapid increase of the pressure.Such a burst of the pressure is called waterhammer, and it is one of phenomena peculiar to unsteady pipeline flows.On the other hand, the generation of the high pressure due to waterhammer becomes a pressure wave, and it propagates with the sound speed.The rapid transmission of a pressure wave is another phenomenon characteristic of unsteady pipeline flows.In this paper, we show how these phenomena can be treated mathematically.In particular, we introduce a new model to capture phenomena observed experimentally and show how it can be analyzed numerically.
In Section 2, we introduce some experimental results to illustrate unsteady pipeline flows.A pipeline was connected to a water tank at one end and a valve was set at the other end.The valve was set open initially to allow a uniform flow to form.Then the valve was closed suddenly, and waterhammer was generated.We present the temporal change of the pressure which was measured in the transition.In Section 3, we propose a new model to analyze unsteady liquid pipeline flows numerically.We show that an expression derived from the Darcy-Weisbach equation corresponds to the radial component of the shear stress in a steady radially symmetric laminar flow, and that the axial component should be retained in the modeling of unsteady flows.In Section 4, we turn to a numerical aspect of the problem.We show how upstream difference approximations can be applied to the momentum equation and the continuity equation.We also present some numerical results and examine their correspondence with experimental results introduced in Section 2.

Description of unsteady pipeline flows.
In this section, we present some results obtained in an experimental study of unsteady pipeline flows.Three pipes of lengths L 1 = 26.87m, L 2 = 16.9 m, and L 3 = 10.6 m were connected to form a pipeline.The total length L = L 1 + L 2 + L 3 of the pipeline was 54.37 m.The diameter D and the wall thickness e of each of those pipes were equal to 0.05 m and 0.003 m, respectively.It was connected to a water tank at one end, which we call the upstream end, and a valve was set at the other end, which we call the downstream end.The first pipe of length L 1 was set horizontally, and its height h 1 , relative to the valve, was 3.3 m.The height h 2 of the water level of the tank, relative to the upstream end, was 2.99 m.The total height h = h 1 + h 2 of the water level of the tank, relative to the valve, was 6.29 m.An outline of the pipeline is shown in Figure 2.1.There was actually a part of the pipeline between the first pipe and the tank.It stretches vertically from the end of the first pipe, and then bends perpendicularly, as shown in the figure.However, its effect is negligible as far as unsteady flows are concerned, and we may regard it as a part of the tank.
Initially, the valve was set open to allow a uniform flow to form.Then the valve was suddenly closed, and the initial uniform flow turned to an unsteady flow.In the meantime, the temporal change of the pressure was measured at the point 0.1 m upstream from the valve, that is, x = 54.27 m.Figures 2.2 and 2.3 show two experimental results with different values of the velocity V 0 (m/s) of the initial uniform flow: (a) V 0 = 0.089 m/s (Figure 2.2) and (b) V 0 = 0.122 m/s (Figure 2.3).We set the moment the pressure head starts rising to be the initial time t = 0 s.In those figures, the pressure p (Pa) is represented in terms of the pressure head H (m). The relationship between the pressure and the pressure head is represented by the following equation: (2.1) Here ρ (kg/m 3 ) is the density of the fluid, g (m/s 2 ) is the gravitational acceleration, and z (m) is the height of the point on the pipeline where the pressure is measured.
In both cases, the pressure head rapidly increases from the initial level which approximately equals the height h = 6.29 m of the water level of the tank relative to the valve after the valve was closed.By t = 0.02 s, the pressure head becomes as high as 18 m and 22 m in cases (a) and (b), respectively.In general, the relationship between the change of the pressure head ∆H (m) and the change of the velocity ∆V (m/s) is represented by the following equation [15]: Here a (m/s) denotes the wave speed, that is, the velocity of the pressure wave.The wave speed is theoretically given by where K (Pa) is the bulk modulus of the elasticity of the fluid, e (m) is the pipe wall thickness, E (Pa) is Young's modulus of elasticity, and µ is Poisson ratio; −µ is the ratio of the lateral unit strain to the axial unit strain of the pipe.Theoretical values of the parameters are K = 2.03×10 9 Pa, E = 196×10 9  We can interpret what happened inside the pipeline as follows.Immediately after the valve was closed, the pressure wave starts propagating upstream.As soon as the pressure wave reaches the upstream end, the high pressure developed over the entire stretch of the pipeline exerts the force on the fluid toward the tank, and then the reverse flow is generated.The generation of the reverse flow at the upstream end propagates downstream towards the valve.As the reverse flow is generated, the pressure starts decreasing.Thus the pressure over the entire stretch of the pipe returns to the original level, when the generation of the reverse flow reaches the valve.Our experimental results show that it took 0.07 s to reach this stage.Meanwhile, the reverse flow still continues even after the pressure over the entire stretch of the pipeline returns to the original level, and the pressure starts decreasing further at the valve.The low pressure built at the valve also becomes a wave and propagates upstream towards the tank.As soon as the wave reaches the upstream end, the low pressure developed over the entire stretch exerts the force on the fluid towards the valve, and the fluid starts moving towards the valve.The generation of the flow towards the downstream end starts at the upstream end, and propagates towards the valve, and the as soon as the flow is generated the pressure returns to the original level.The generation and the propagation of the pressure waves are repeated and the state of the flow approaches the stationary state.In the following sections, we illustrate how such unsteady pipeline flows can be analyzed numerically.

Introduction of the turbulent effect into the modeling of liquid pipeline flows.
We turn to an analytical aspect in studies of unsteady pipeline flows.In particular, we propose a new model for a turbulent flow.In studies of unsteady liquid pipeline flows, the momentum equation and the continuity equation are often analyzed numerically.Here x (m) represents the distance from the origin of a one-dimensional coordinate set along the pipeline, and t (s) represents time.The velocity and the pressure head of the flow are denoted by V (m/s) and H (m), respectively.The subscripts x and t denote the partial differentiation with respect to x and t, respectively.The Darcy-Weisbach friction factor is denoted by f .The last term in the right-hand side of the momentum equation represents the shear stress, which was developed from the Darcy-Weisbach equation for steady flows [15].In a three-dimensional xyz-space, the shear stress exerted on a Newtonian fluid in the x-direction is given in terms of the expression where u is the x-component of the velocity, and ν (m 2 /s) is the kinematic viscosity whose value for water is 10 −6 m 2 /s.Suppose that the flow is radially symmetric.In the cylindrical coordinates (x,r ,θ), where the shear stress can be represented in terms of the following expression: ( Moreover, when the velocity is independent of x, the first term is dropped, and it becomes (3.6) On the other hand, the x-component of the velocity of the Poiseuille flow, a steady radially symmetric laminar flow, is given by where u 0 is a constant.Its average ū over a cross section is given by Now the shear stress can be represented by Thus, in a steady radially symmetric laminar flow, the shear stress is proportional to the average velocity.Meantime, for the experimental results shown in Section 2, the velocity of the initial flows are greater than 0.08 m/s.If we set ν = 10 −6 m 2 /s, the Reynolds number is at least 4 × 10 4 , for the diameter D of the pipeline was 0.05 m.In general, a pipeline flow whose Reynolds number is greater than 1/3×10 4 is a turbulent flow, and the initial flows in the experiment should be turbulent.When the flow is turbulent, the shear stress is known to be proportional to the square of the average velocity [12].
The foregoing discussion suggests that the fourth term on the left-hand side of (3.1), f V |V |/2D, corresponds to the radial component of the shear stress of a steady flow.For an unsteady flow, the term proportional to the second derivative of the velocity component with respect to x must be retained, and the momentum equation should become Here ν 0 is a constant that represents a generalized kinematic viscosity called the eddy viscosity.In the following section, we illustrate how the system of partial differential equations (3.2) and (3.11) can be analyzed numerically to evaluate the velocity and the pressure head.

Numerical analysis of unsteady pipeline flows.
In this section, we show how the system of the partial differential equations (3.2) and (3.11) can be solved numerically.We also present some numerical results to show how the effects of the axial component of the shear stress appear.In order to solve the system of partial differential equations (3.2) and (3.11) numerically, we convert it to a form for which upstream difference approximations are suitable [7].The change of variables converts (3.2) and (3.11) to the following system which is suitable for upstream difference approximations: Here the functions Q ± (V ) are defined by For a typical pipeline flow, the wave speed a is greater than 10 3 m/s, whereas the magnitude of the velocity is at most 10.0 m/s, and it may be assumed that V − a < 0 and V + a > 0. Then a forward difference approximation can be applied to (4.2) and a backward difference approximation can be applied to (4.3) efficiently [2].Divide the interval [0,L] into N equally spaced intervals, and let ∆x denote the length of the intervals: Let x i = i∆x (i = 0, 1, 2,...,N).Choose ∆t > 0 and let t j = j∆t (j = 0, 1, 2,...).Denote approximate solutions V and H of the system (3.2) and (3.11) at (x, t) = (x i ,t j ) by V i,j and H i,j , respectively.Suppose that Y i,j and Z i,j are given by Then Y i,j and Z i,j are approximate solutions of the system (4.2) and ( 4.3) at (x, t) = (x i ,t j ).Then we obtain the following difference equations: The functions Q − (V i,j ) and Q + (V i,j ) now become where In order to obtain our numerical results, we introduced the effect of pressure gradient into the friction term, and set Here f s denotes the Darcy-Weisbach friction factor, and f d denotes the factor of the friction due to the divergence of the gradient of the pressure head from the slope of the pipe.For Q − (V i,j ), (H x ) i is given by and for Q + (V i,j ), it is given by The height z of a point on the pipeline is given as the following function of x: where Our interest is to see how the axial component of the shear stress affects numerical results.On the other hand, the eddy viscosity varies from one problem to another, and horizontal eddy viscosities range up to 10 11 times the kinematic viscosity [3].Here we tried our model setting ν 0 = 5.0 m 2 /s and V 0 = 0.122 m/s.Then, in order to find an appropriate value of f d , we tested 201 evenly spaced points in the interval [10,50] when the other parameter values were fixed, and selected the one in which the error between the numerical solution and the experimental data was minimized.Thus, we set the following values of the parameters: Here T is the time it takes for the valve to close completely.The values of the other parameters were as described before.Approximate solutions were computed with the initial value   In the following section, we discuss the significance of those numerical results.

Conclusion.
.3 shows that both the numerically computed pressure heads increased almost to 21 m by 0.02 s.Meanwhile, the pressure head based on the system (3.2) and (3.11) increased to 22 m by 0.05 s and stayed at this level for a while, whereas the pressure head based on the system (3.1) and (3.2) stayed below the other one.Then they started decreasing.Figure 4.4 shows that the pressure head based on the system (3.2) and (3.11) decreased to −7.5 m by 0.12 s and that it stayed at that level until it started increasing again approximately at t = 0.13 s.On the other hand, the pressure head based on the system (3.1) and (3.2) stayed above the other one.Note that the pressure head based on the system (3.2) and (3.11) was closer to the experimental result when it stayed above or close to 21 m between t = 0.02 s and t = 0.07 s, or when it stayed below or close to −6.5 m between t = 0.11 s and t = 0.15 s.In summary, the system (3.2) and (3.11) captured the characteristics of an unsteady pipeline flow observed experimentally.
In addition to waterhammer and propagation of pressure waves, another phenomenon characteristic of unsteady pipeline flows is a formation of two-phase flow of the mixture of liquid and liquid-vapor due to an extreme decrease of the pressure.When the pressure decreases to a certain level, called vapor pressure, liquid-vapor can be generated in a liquid flow, and the two-phase flow can form.Vapor pressure is known to be 32 hPa, and it corresponds to the pressure head H = −10 m.Numerical techniques to capture the phenomena have been developed for the past two decades [1,4,5,6,8,9,10,11,13,14]. In analysis of the two-phase flows which arise from liquid flows, it is important to know precisely when the pressure decreases to vapor pressure, for that is when a liquid flow starts turning to a two-phase flow.We have demonstrated a value of our new model also in this respect.
Pa, and µ = 0.3.For ρ = 1000 kg/m 3 , we find that the theoretical value of a is approximately equal to 1324.54 m/s, and for this value of a (m/s), theoretical values of ∆H obtained by substituting ∆V = −0.089m/s and ∆V = −0.122m/s in (2.2) are approximately ∆H = 12.02 m and ∆H = 16.48 m, respectively.Note that these predicted values of the increase of the pressure head based on (2.2) match well with the values obtained experimentally: (a) ∆H ≈ 18−6.29 = 11.71m and (b) ∆H ≈ 22 − 6.29 = 15.71 m.After the waterhammer occurs, the pressure head remains almost constant approximately until t = 0.07 s, and then it decreases rapidly.The pressure head decreases to −5 m in case (a), and −10 m in case (b), and the pressure head starts increasing after it reaches the lowest level.

Figures 4. 1
Figures 4.1, 4.2, 4.3, and 4.4 show the numerical results based on the system (3.1) and (3.2), and the system (3.2) and (3.11) for V 0 = 0.089 m/s and V 0 = 0.122 m/s, respectively.They show the temporal change of the pressure head at x = 54.27m which corresponds to the position where the pressure was measured experimentally to obtain the results shown in Figures 2.2 and 2.3.In those figures, numerical solutions of the system without the axial component of the shear stress, the system (3.1) and (3.2), and

, 4 . 2 , 4 . 3 ,Figure 4 . 3 .Figure 4 . 4 .
Figures 4.1, 4.2, 4.3, and 4.4 show the numerical results based on the system (3.1) and (3.2), and the system (3.2) and (3.11) for V 0 = 0.089 m/s and V 0 = 0.122 m/s, respectively.They show the temporal change of the pressure head at x = 54.27m which corresponds to the position where the pressure was measured experimentally to obtain the results shown in Figures 2.2 and 2.3.In those figures, numerical solutions of the system without the axial component of the shear stress, the system (3.1) and (3.2), and