A finite volume study for pressure waves propagation in a straight section of pipeline with caviation

The main objective of this research was to study the pressure waves propagation generated by a sudden closure of a valve in a straight pipe. The physical model consisted of a head tank that can be pressurized with air, and a copper pipe with a fast-closing ball valve on the downstream end of the line. The cavitation and fluid-structure interaction phenomena were integrated analytically into the one-dimensional continuity and momentum equations, by assuming that the fluid density and the flow area vary with pressure. These equations were solved through a high resolution finite volume method, in combination with others numerical methods such as Taylor series expansion, Newton method, Simpson’s Rule and quadratic interpolation. Due to the complexity of the solution procedure, a computational code in FORTRAN 95 language was developed in order to obtain numerical solutions. Several discretizations of the computational grid were achieved to assess their impact on the solution. The model was validated with experimental data and analytic results obtained by other researchers. Several pressure values, in different points of pipe, were compared, and an excellent agreement was found for both cases. NOMENCLATURE A = pipe flow area Ao = pipe flow area for the unpressurized state A’ = derivate of A with respect to pressure a = wave speed through homogeneous fluid contained in pipe a = wave speed through homogeneous fluid contained in pipe for initial conditions c1 = pipe restrain factor Do = pipe inside diameter E = modulus of elasticity for the pipe material e = pipe wall thickness F = pressure continuous function k = bulk modulus of water L = pipe length Nz = number of computational cells n = time level P = fluid pressure Psat = saturation condition of pressure Int. Jnl. of Multiphysics Volume 7 · Number 4 · 2013 259


INTRODUCTION
The pressure waves caused by hydraulic transient in a pipeline, like a sudden closure of a valve, can become in a stress increase and generated noise, vibrations and damages in the pipe materials.For these reasons it is consider an undesirable phenomenon in pipelines and is necessary to know its effects by including them into the design process of any hydraulic system.It also can cause localized and distributed vaporization which may have a devastating effect on a pipeline system.A famous example is the Oigawa hydropower accident in Japan, were a fast valve-closure, due to the draining of an oil control system during maintenance, caused an extremely high-pressure water hammer wave that split the penstock open.The resultant release of water generated a low-pressure wave resulting in substantial column separation that caused crushing (pipe collapse) of a significant portion of the upstream pipeline due to the external atmospheric pressure [1].This research studies the pressure waves propagation along a fluid, generated by hydraulic transient in a straight pipe section.It also includes the phenomena of fluid structure interaction and cavitation analytically, by considering that the fluid density and the flow area vary with pressure.The most widely used approach for modeling cavitation in pipelines, which is not used herein, is the Discrete Vapor Cavity Model (DVCM) [2], in which the cavity is assumed to occupy the entire cross section of the pipe on specific positions.This model is usually applied in combination with a Finite Difference Method, being its major drawback the fact that solutions tend to become artificially spiky as the mesh is refined.This problem is addressed in this work through the implementation of the MUSCL-Hancock method, which is a High Resolution Finite Volume Method that gives second-order accuracy on smooth solutions and integrates a slope limiting technique to work near to the discontinuities to prevent nonphysical oscillations in the solution [3].
Second order accurate methods give a very good accuracy on smooth solutions, but fail near discontinuities, where oscillations are generated.Upwind Methods have the advantage of keeping the solution monotonically varying in regions where the solution should be monotone, even though the accuracy is not very good.The idea with High Resolution methods is to combine the best features of both methods.Second order of accuracy is obtained where possible, but without insist on it in regions where the solution is not behaving smoothly [4].
The application of this method has important advantages: solutions do not develop artificial spikiness when refined computational meshes are employed, the predicted pressure amplitudes do not exhibit important overshooting, and since the fluid is treated as homogeneous, no special considerations are required for the treatment of boundary conditions adjacent to dynamic components such as closing valves or pumps [3].
Validation of the mathematical model and the finite-volume solution method is provided by comparing results against experimental data obtained from the water hammer experiment of Simpson and Wylie [5], the analytic results obtained by Chaiko [3] and the analytic results obtained by Da Silva [2].

PHYSICAL MODEL
The physical model for this analytical study is defined by the experimental apparatus used by Simpson and Wylie [5] to study liquid-column separation on the upstream side of a suddenly closed valve.The pipeline is considered to be horizontal.The system consists of a head tank that can be pressurized with air, and a copper pipe with a fast-closing ball valve on the downstream end of the line.The pipe inside diameter is 19.05 mm and the wall thickness is 1.588 mm.Pipe length from the head tank to the valve is 36 m.Water temperature was maintained at approximately 24°C.The wave speed for the liquid filled line is determined 1298 m/s.The initial fluid velocity is 0.332 m/s and the reservoir pressure head is 23.41 m.A simplified schematic of the physical model is shown in Figure 1.

MATHEMATICAL MODEL
The mathematical model is based on a one dimensional, homogeneous-fluid continuity and momentum equations (1) (2) where P(z,t) is the pressure, ρ(P) is the fluid density, u(z,t) is the fluid velocity , A(P) is the pipe flow area, t is the time, z is the axial distance along the pipe and A'=dA/dP (which must be constant for Eqn.(2) to applies).For the integration of the fluid-structure interaction phenomena, the flow area can vary with fluid pressure due to the next stress-strain relation [6] (3) where D o and A o denote the inside diameter and the pipe flow area or the unpressurized state, e is the thickness of the pipe, E is the modulus of elasticity for the pipe material (for copper ) and c 1 is the piping restrain factor for a pipe anchored against longitudinal movements (4) where v is the Poisson's ratio, which is 0.35 for the copper pipe [7].
The density variations depend only on the pressure and are assumed to occur at constant entropy.For the sub-cooled liquid region, where P is greater than P sat (5) where T o is the initial temperature of the fluid, k is the bulk modulus of the fluid (2.07GPa for water at ordinary temperatures [24]).For the saturated region, where P is less or equal to P sat : (6) where υ is the fluid specific volume, the subscripts f y g refers to saturated liquid and vapor, respectively, and X is the quality of the mixture.

BOUNDARY CONDITIONS
Once the transient begins, caused by the sudden closure of the valve, the pressure is considered to remain constant at the upstream end of the pipe (z=0), therefore, in this position the pressure is constant and equal to the tank pressure P T (7) On the downstream end (z=L), the fluid velocity is an specific function of time during the valve closure stroke and zero for the remaining part of the transient.(8) where t 1 is the time interval during which the data are applied as a boundary condition (closure time of the valve) and is a specific velocity boundary condition derived from the experimental data [5] by curve fitting (Figure 2) resulting for t < 0.02s: (9) For 0.02375 s < t ≥ 0.02 s: (10) For 0.03 s < t ≥ 0.02375 s: (11) V(t)=-44326t +3767,5t -106,82t+1,0107 3 2 (1 )

ANALYTICAL MODEL
The governing continuity and momentum equations are solved by the MUSCL-Hancock method, which is a high-resolution finite volume method and prevent numerical oscillations by using slope limiting.The solution domain is divided into Nz computational cells of equal length ∆z.To ensure the numerical stability of the method, ∆t is calculated according to the CFL condition (12)

EQUATIONS INTEGRATION
The governing equations in their conservation form are numerically integrated on the finite volume grid.Central differences in both space and time are used to obtain the next finite difference formulas where n represents the time level.

DATA RECONSTRUCTION
The first step for applying this method is the linear reconstruction of the data.The variables values ρA and ρuA are found in the cells neighbor for the current time level n.For interior cells (j = 2,…, Nz-1), slope limiting is accomplished through the application of the minmod slope limiter.The slope m of ρA and ρuA respectively, on interior cells at time n is where the definition of the minmod function is (17) For the first and last cells (j=1 y j=Nz), slope limiting is not used since they only have one neighbor cell.Inward sloping difference formulas are employed to obtain their slopes Using the previous calculated slopes, the quantities ρA and ρuA are obtained on the left and right side, respectively, of each cell boundary according to

TAYLOR SERIES EXPANSION
The second step of the solution method is the Taylor Series Expansion of the linear reconstructed data on n time level where the slopes obtained in the previous step are used again.The fluid velocity u is calculated in t n+1/2 with the quantities ρA and ρuA.In the subcooled liquid region, direct calculation of pressure is accomplished by the solution of a quadratic equation.In the saturated two-phase region, an iterative numerical procedure is used to calculate the pressure.A guess is made for the pressure, and then the flow area and density are computed from Eqn.
(3) and the known value of the product ρA.Newton's method is applied to Eqn. ( 6) in order to recalculate pressure from the density.The calculated pressure is then compared against the guessed value to determine if convergence has been achieved, or if another iteration is required [3].

RIEMANN PROBLEM
With the quantities of u and P on every side of the cell boundary, the Riemann problem is solved to have only one value of each variable in each cell boundary.In the subcooled liquid region, where P is greater than P sat (T o ), the weak dependency of density and wave speed pressure for F(P) is neglected to obtain [3] (30) For interior cell boundary (j = 2,…,Nz-1), where rightward and leftward propagation characteristics intersect [3] (31) (32) On the left boundary (z = 0), the pressure is equal to the tank pressure P T .This boundary

RESULTS
When the transient begins with the closure of the valve, it causes a high pressure wave which propagates to the upstream of the pipe.Initially the pressure at the valve gradually increases to its maximum value, just when the valve is completely closed to 0.03 s.Then, the pressure remains constant until the first pressure wave is reflected back from the tank to produce the first pressure drop.From that moment the pressure begins to drop sharply until reaching the saturation conditions of the liquid and produce the appearance of the first steam cavities at about 0.077 s.
Due to the constant entropy condition implemented in the model, once it reaches the saturation conditions, the pressure continues decreasing to a certain point and then begins to increase again until reaching the saturation pressure.This variation is so small that cannot be appreciated in Figure 3, which suggests that the vapor fraction in the mixture is quite small.This behavior matches the physical description of the cavitation phenomenon; therefore, it is considered that the proposed model provides satisfactory results in this regard.( ) As time increases, results are getting away from the ones obtained by Chaiko [3], who incorporates friction in his one-dimensional model.Therefore, he gets values of pressure that slightly decrease as time increases, due to the dissipative effect of friction.For this reason, his results are closer to the experimental ones for higher times.Figure 4 compares calculated pressure at z = 9 m, near from the tank, using 60 computational cells.Again, it shows an excellent agreement with the experimental data.Initially, pressure remains constant from initial value until 0.02s approximately, when the first high pressure wave arrives to this point and causes pressure increase.Evidently, the effect of cavitation is less severe in this position of the pipeline than it is in the valve.The saturation state is reach by short period of time.Comparing the results using the finite-volume method and the ones obtained by Da Silva [2] using DVCM, as shown in Figure 5, total agreement is achieved between them until 0.15 s.From this point, several waves start to propagate at the same time in the system, which increases the discontinuities and produces numerical oscillations for the DVCM.For shorter lengths of pipe, more time the pressure waves spread for the same time level, as shown in Figures 6, 7 and 8.This is because the length of the line is shorter, which suggest that a shorter time is also required for the wave traveling towards one end of the pipe and returns.This turns in a decrease of the cavitation severity, by not allowing the vapor fraction in the mixture to increase.The fluid reaches the saturation pressure, however, does not remain in that state for a long time, when the high pressure wave coming back reaches again.

CONCLUSIONS
The finite volume method implemented, which incorporated slope limiters, allowed to obtain a stable solution and also with a low computational cost.The solutions for all the cases studied were obtained in seconds and all of them showed a very good agreement with the experimental data, even with as few as 20 computational cells.The most severe conditions of cavitation were found at the valve.The friction effect on the pressure consider by Chaiko [3] is not significant for shorts periods of time.
for the unpressurized state A' = derivate of A with respect to pressure a = wave speed through homogeneous fluid contained in pipe aᐉ = wave speed through homogeneous fluid contained in pipe for initial conditions

Figure 1 .
Figure 1.Simplified schematic of the physical model

Figure 2 .
Figure 2. Curve Fitting for fluid velocity on valve neighbor from experimental data.
velocity of the fluid on the boundary z = L is specify by Eqn.(8), and the fluid pressure is given implicitly by the relation (34)

Figure 3
Figure3compares calculated pressure at the valve (z = L = 36 m) against those obtained by Chaiko[3] and the experimental data of Simpson and Wylie[5].It is done using 40 computational cells.At this point of the pipeline, the most severe condition of cavitation is expected.The results show very good agreement with both comparisons.The mathematical model used does not consider the slope of the pipeline of the physical model, and that's the reason why the values of pressure obtained are lightly higher than the experimental ones, especially in the initial time intervals.

Figure 4 .
Figure 4. Comparison of finite-volume calculation of pressure at z = 9 m using 60 computational cells against experimental data and Chaiko's results [3].