Numerical Simulation of Laminar Water Hammer Flow

Title: Numerical Simulation of Laminar Water Hammer Flow in a Pipe with Varying Cross-section

PRANAB NARAYAN JHA, Clearview Subsea, Houston, Texas, USA 77084 Email: [email protected] JIANFENG ZHANG, Clearview Subsea, Houston, Texas, USA 77084 Email: [email protected] (author for correspondence) CHARLES DALTON, University of Houston, Houston, Texas, USA 77204 Email: [email protected]

Running Head (i.e., the text that appears in the top margin of the published page): Laminar water hammer in varying cross-section pipe.

This is the Authors' Original Manuscript of the article published online by Taylor & Francis in Journal of Applied Water Engineering and Research on 16 Feb 2017, available online: http://iahr.tandfonline.com/doi/abs/10.1080/23249676.2017.1287016.

1

Title: Numerical Simulation of Laminar Water Hammer Flow in a Pipe with Varying Cross-section 1. Abstract Numerical simulation of water hammer effect in a pipe with varying cross-section is presented for an initially laminar flow by solving Navier-Stokes equations in time using Ansys Fluent software. Results for pressure variation with time are compared to available literature to validate the numerical method. The pressure signal is presented at various locations in the pipe and compared with existing studies qualitatively. Pressure and velocity profiles are shown at various times during the pressure wave along the length of the domain which indicates fluctuation of these variables as the wave passes through the changing geometry. Attenuation of the pressure amplitude was found to be slower in the current geometry compared with experimental results from the literature due to a more gradual change in cross-sectional area. The initial pressure wave and its interaction with the geometry changes in the domain is discussed in detail and comparisons are made with existing literature. Keywords: Water hammer; numerical simulation; varying cross-section; laminar flow; pressure wave. 2. Introduction Water hammer is a common occurrence in piping systems. Due to sudden closure of valves, partial blocking of lines, etc. the pressure wave generated in the system can potentially damage the pipes and fittings. Many experimental studies have been conducted in the past and one-dimensional computational models have been developed to predict the effect of pressure transients in liquid filled pipes and piping systems. But very few studies exist that have simulated the transient flow in a pipe and even fewer concerning changing cross-sectional area. One of the earlier results are presented by Holmboe & Rouleau (1967). They conducted experiments on the problem of water hammer using a 1 inch diameter pipe with laminar flow (𝑅𝑒 ≈ 82) and used a quick-closing valve at one end of the pipe. They provided pressure transient data which has served as a very good reference and validation data set for studies carried out by other researchers. Studies such as the one conducted by Zielke (1968) have presented one-dimensional models of transient flow using numerical techniques like the method of characteristics. Anton et al. (2001) reviewed the literature for numerical models used for predicting pressure transients in single pipe and piping systems. Using a 1D model Brunone et al. (2004) studied the pressure decay in laminar transient flow and concluded that their model needed to be modified in order to account for the pipe material. This suggests 2

that more detailed calculations need to be carried out to understand the water hammer phenomenon, even for laminar flow. Riasi et al. (2009) are among the few researchers who have conducted numerical simulation of laminar and turbulent flows through a pipe by solving the unsteady Navier-Stokes equations. They compared their laminar flow results in a straight pipe with Holmboe & Rouleau (1967) and found a good agreement for transient pressure data. Their results described the region of reverse flow and strong gradients near the pipe wall during flow deceleration. They also observed the change in point of zero velocity and inflection point in the velocity profile over the transient pressure change. Jinping et al. (2010) also conducted CFD simulations of turbulent flow in a long straight pipe and their geometry was very similar to that used by Riasi et al. (2009). Meniconi et al. (2012) presented both experimental and numerical data of the interaction of a pressure wave in a pipe with changing cross section. They conducted experiments using pipes of different diameters in series and joined using sudden contraction-expansion connectors as well as gradual conical connectors. This study will only qualitatively compare the results with their data because their experiments were carried out at turbulent flow speeds. This study presents a numerical simulation of laminar flow in a pipe with varying cross-section. The present arrangement is much shorter than previous studies conducted by most of research groups. This arrangement is part of a sensor and has a measurement section (hereafter referred to as MS) of rectangular cross section between two circular pipe segments. The MS is connected to the pipes using nozzle-diffuser connectors. The software package Ansys Fluent (Ansys Inc., 2015) was used for the current simulation. 3. Problem Description 3.1. Geometry Figure 1 shows the geometry used in the simulation. Two pipes of 25.4 mm (1-inch) inside diameter and given lengths are connected to a measurement section (MS) which has a rectangular cross section with rounded edges. The cross-sectional area ratio of the pipe to the MS is about 10. A nozzle on the upstream end and a diffuser on the downstream end ensure a smooth transition from the circular cross section of the pipe to the MS. Both the nozzle and the diffuser have a half angle of 5o. This angle was chosen based on a separate study conducted to find the optimum angle for the flow through the current geometry. This design minimized the overall pressure drop across the axial length (from inlet to outlet) and eliminated any recirculation region in the diffuser section.

3

The MS is 101.6 mm (4 inches) in length. The coordinate system origin for this geometry lies at the center of the MS with the z-axis being the flow direction. Thus, the MS spans +/- 1 mm in the normal, +/- 50.8 mm in stream-wise (axial) and +/- 12.7 mm in the span-wise direction relative to the coordinate system origin. Related dimensions are also indicated in Figure 1. This arrangement is part of a sensor that will be deployed on a subsea pressure vessel.

Figure 1. Geometry and dimensions (in mm) of pipe with gradual changing cross-section. Front view (top) and top view (bottom).

This geometry is somewhat different than the one used by Meniconi et al. (2012). They conducted experiments on ‘in-series’ visco-elastic pipe system of 93.3, 38.8 and 93.3 mm diameter (area ratio 5.8) and lengths were 54.49, 100 and 110.44 m, respectively. Both connections used between the pipes were sudden contraction – expansion and gradual conical connectors (half angle 36o). Another set of experiments was conducted using a short pipe of 38.8 mm diameter and 0.48 to 0.06 m length connected to two pipes of 93.3 mm diameter and greater than 75 m length. This configuration was referred to as a ‘partial blockage’ system. They only tested their system with sudden contraction – expansion connections. Apart from having a varying cross-section, the arrangement used in the present study might have significant industrial applications in systems where quick closing valves are used in conjunction with pipes having short axial dimensions.

4

Figure 2. Top view of the mesh on the boundaries in the vicinity of upstream pipe, nozzle and measurement section. 3.2. Mesh Figure 2 shows the mesh used in the present study. The upstream-downstream pipes were meshed using hexahedral elements and the nozzle-diffuser connectors using tetrahedral elements. A seven cell high boundary layer mesh was used near all the pipe walls. It can be observed that the MS had finer mesh to resolve the higher velocity expected in this region due to smaller flow area. The mesh was progressively refined from the upstream-downstream pipe towards the MS. The total mesh cell count for the present study was 0.69 million. A separate mesh dependence study was conducted using mesh resolutions of 0.31 million and 1.68 million cells. It was concluded that the present mesh was sufficient to resolve the laminar flow in the system and results are not reported here. 4. Numerical Method 4.1. Governing equations The flow in consideration was modelled as compressible laminar flow. The mass and momentum balance equations are 𝜕𝜌 𝜕(𝜌𝑢𝑗 ) + = 0 and 𝜕𝑡 𝜕𝑥𝑗

(1)

𝜕(𝜌𝑢𝑖 ) 𝜕(𝜌𝑢𝑗 𝑢𝑖 ) 𝜕𝑃 𝜕 𝜕𝑢𝑖 + = − + (𝜇 ) 𝜕𝑡 𝜕𝑥𝑗 𝜕𝑥𝑖 𝜕𝑥𝑗 𝜕𝑥𝑗

(2)

respectively. Here, ρ is the fluid density, uj is the fluid velocity in the jth direction, P is the static pressure, and μ is the fluid viscosity. The pipe material is considered to be made of stainless steel with a yield strength, E = 180 GPa. To account for the effect of pipe elasticity on the fluid a modified bulk modulus Kf’ of elasticity for the fluid is introduced,

5

𝐾𝑓′ =

𝐾𝑓 . 𝐾𝑓 𝐷 1 + 𝑒𝐸

(3)

Here Kf = 2.2 GPa is the bulk modulus of elasticity of water, D is the pipe diameter, and e is the pipe thickness. A pipe thickness of 4.5 mm (0.179 in) was assumed for modeling pipe elasticity (this had no effect on the geometric dimensions used in the simulation). Pressure wave speed depends on the bulk modulus. For a rigid pipe without taking elasticity into account, a higher value of bulk modulus would be obtained, leading to higher wave speed which would be unrealistic. The pressure wave speed a is given as 𝐾𝑓′ 𝑎=√ . 𝜌

(4)

The flow is assumed to be compressible and isothermal. The density variation in the fluid domain depends on the pressure change and is computed at every cell from 𝑑𝜌 =

𝜌𝑑𝑃 𝐾𝑓′

(5)

Equations for speed of sound and density of fluid were implemented in Fluent by a user-defined function. 4.2. Boundary Conditions and Discretization schemes The governing equations were solved over the discretized mesh domain using Ansys Fluent (Ansys Inc., 2015). A total pressure inlet boundary was prescribed on the upstream left end of the domain and a static pressure outlet condition on the downstream right end (Figure 1). After establishing an initial steady flow field in the domain, the downstream end was changed to wall boundary condition and an unsteady simulation was conducted starting at t = 0. All other surfaces were prescribed as wall boundaries. The Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm was used for pressure-velocity coupling. The pressure equation was solved using a second-order scheme and the nonlinear convective terms in the transport equations were modeled using the third order MUSCL (Monotone Upstream-Centered Schemes for Conservation Laws) scheme. The convergence criteria set for pressure and velocity components were to obtain normalized residuals below 1e-5 at each time step. A first-order implicit method was used for time derivative terms. A time-step of 1e-5 seconds was chosen to resolve the pressure wave moving from one end of the domain to the other. This choice of time-step was based on the fact that the current domain was smaller in axial dimension than reference numerical work discussed here. Riasi et al. (2009) used a time-step of 1e-4 s because their domain was longer. They could resolve the pressure wave as it moved across the domain extent and simulate the flow 6

field accurately. The current time-step is based on a similar logic considering the current geometric dimensions. 5. Results and Discussion 5.1. Validation The numerical scheme was validated using data from Holmboe & Rouleau (1967) for laminar flow in a straight pipe (no cross section area change). A two-dimensional axisymmetric pipe domain was used for this purpose with an inside diameter 0.0254 m (1 inch) and length 36.09 m. The domain was discretized using 400 divisions in axial and 50 divisions in radial direction. This mesh resolution was chosen based on the mesh resolution study presented by Riasi et al. (2009) and similar numerical method was adopted. A total pressure at the inlet and atmospheric (static) pressure at the outlet boundaries were specified. The kinematic viscosity of the liquid was 𝜈 = 3.967 × 10−5 𝑚2 /𝑠. The inlet pressure was chosen such that a Reynolds number of 𝑅𝑒 = 82 was obtained, the same as both Holmboe and Rouleau (1967), and Riasi et al. (2009). Figure 3 shows the pressure transient for this validation case compared favorably with results from the reference publication (Holmboe & Rouleau, 1967). This study suggested that the time-step (10-4 s) was small enough to capture the pressure wave in the long pipe and that the bulk modulus was modeled correctly. Also inferred is that the numerical discretization of equations (1) and (2) was done accurately. The pressures near the valve and the pipe centerline match with the experimental results. Given good agreement between the two sets of data, this numerical scheme was applied to simulate flow in a pipe with changing cross-sectional area and results are discussed in next section.

7

Figure 3. Transient pressure comparison for validation. “Exp.” indicates Experimental data of Holmboe & Rouleau (1967) 5.2. Varying Cross-Section Geometry

Figure 4. Transient pressure recorded at various locations in the domain.

In this section results from simulation conducted of flow through the geometry discussed in section 3.1 are presented and discussed. The inlet was prescribed a total pressure boundary condition such that Re = 82 was maintained in the pipe. This Reynolds number is based on the inlet pipe diameter 0.0254 m, average velocity 0.0032 m/s, sea-water density (ρ = 1027 kg/m3) and viscosity (μ = 0.00175 kg/m-s). As with the validation case, a steady flow was established in the pipe before valve closure. The valve closure action was simulated by closing the outlet instantly at t = 0. Figure 4 shows the pressure

8

head (P/ρg) in the pipe at three different locations – near the valve, at the center of the measurement section, and in the nozzle upstream of the measurement section. It can be observed that the pressure wave behaves in a similar trend at the three locations, although the amplitude of the wave is lower in the MS center (center) and even lower in the nozzle located upstream of the MS. This reduction in amplitude of the wave was not noticed in the straight pipe calculation carried out by Holmboe & Rouleau (1967) (and discussed in the previous section). This decrease can be attributed to the change in cross sectional area ratio and the associated change in velocity between these locations. The first wave near the valve showed multiple peaks and troughs, especially during the initial pressure rise. Though quantitative comparison cannot be made with Meniconi et al. (2012) data, a qualitative analysis suggests that the relative decrease in pressure signals was more magnified for their experiments. By the fifth pressure cycle the amplitude of the pressure signal had decreased by more than 60% for their case, while in the present simulation this decrease was approximately 35%. This is quite likely due to the gradual change in area for the present case (5o half angle) compared with either a sudden contraction – expansion or a 35o half angle nozzle – diffuser. Figure 5 shows the axial velocity variation in time at the same three locations as shown in Figure 4. The reason for decrease in the pressure amplitude as the wave moves from the valve towards the inlet becomes clearer. The velocity at the MS center is much higher than the nozzle or the valve due to the cross-sectional change in area. This leads to the energy carried with the pressure wave to be converted to kinetic energy. Moreover, the mean velocity for all curves in Figure 5 is positive because this plot shows velocity computed at the centerline. A velocity reversal was observed near the wall and discussed later in Figure 7. The fluctuations in crosssectional average velocity (not shown) at any given axial location had a mean value of zero. Such velocity plots and behavior were not reported by Meniconi et al. (2012) and would certainly be difficult to obtain using an experimental setup.

Figure 5. Axial velocity variation at different locations. 9

Figure 6 shows pressure head and axial velocity along the axis of the pipe at different times. No significant variation was observed initially at t = 0.0001 s in the pressure profile along the length of the pipe. This was because, when compared with the profiles at 0.001 s and 0.003 s, the magnitude was quite small. Times of t = 0.001 s and 0.003 s correspond to the high and low amplitudes of the pressure wave, respectively, as observed from Figure 4. The high and low values of static pressure can be seen near the valve (z/D = 17). The reversal of pressure gradient along the pipe can also be observed. Velocity along the axis was initially high in the measurement section and low in the upstream and downstream pipe sections. This was due to the small area of the measurement section compared with the upstream / downstream pipe. During the pressure wave the velocity decreased at 0.001 s, a negative velocity was observed at 0.002 s and switched back to positive at 0.003 s. This indicated outflow from the inlet boundary at 0.002 s. Velocity oscillations also decreased in amplitude with further decreases in the pressure wave (not shown in the figures). This information, which was obtained from the simulation about the pressure rise and velocity reversal, is crucial in designing the piping system, developing the valve operational procedure, and taking physical measurements at the MS.

10

Figure 6. Pressure head (top) and axial velocity (bottom) along the center line.

Figure 7. Axial velocity profile on a spanwise line in the center of MS.

11

Figure 7 shows the velocity profiles at a span-wise location in the center of the MS. The curves correspond to different times during the first pressure cycle. It can be observed that the velocity profile changed from a positive value to a negative value over the cycle, and then back to a positive value. The shape of the profile toward the center axis remained the same, but the biggest change was observed near the walls. Similar phenomenon was observed by Riasi et al. (2009) for their simulations of an initially turbulent flow case. For the present case if L/a is chosen as a characteristic time scale (L being the axial length of the domain), the reversal of velocity profile to a negative value takes place at 3L/a (0.002 s); whereas the profiles showed reversal at 1.6L/a for the turbulent flow simulations of Riasi et al. (2009). This is most likely due to the change in the cross-sectional area.

Figure 8. Pressure signal during the first wave.

Figure 8 shows details of the first pressure wave recorded at four different locations. As expected, the pressure signal near the valve exhibits higher amplitude than in the MS center, upstream nozzle, and the inlet of the computational domain. The diminishing amplitudes can be attributed to the losses due to change in geometry. The ‘pipes in-series’ system used by Meniconi et al. (2012) showed very sharp variation in the pressure signals when the wave crossed a measurement point. In the present simulation, the variations were much smoother but more subdued. For example, at the MS center pressure wave decreased significantly, reaching only about 65% (at t = 0.0004 s) of the pressure that was noticed near the valve (at t = 0.0001 s). This observation was unlike the data presented by Meniconi et al. (2012), in which the relative change in pressure was shown to be about the same or even more in the smaller diameter pipe. But their data did not show such crests in pressure signals for the ‘partial blockage case’. Moreover, their results showed very little difference in pressure signals between the cases with conical (half angle 36o) and sudden 12

changes in cross-sectional area. This behavior can be attributed to two differences from the present simulation – abrupt changes in area and the initial turbulent flow in their experiments. The second and third crests (at t = 0.0006 s and 0.0011 s) in Figure 8 resulted from the wave reflecting from the downstream diffuser and then again from the upstream nozzle of the domain. Subsequently, we see a decrease in pressure signal leading to a negative head which was due to the pressure reflecting back from the outlet section.

Figure 9. Pressure (top) and velocity profile (bottom) during the initial wave.

Figure 9 shows a normalized pressure gradient ((dP/dz)/ρg) and axial velocity along the axis of the computational domain during the initial pressure wave. Right after valve closure (t = 0.0001 s) the sharp increase in gradient can be observed near the outlet. The positive and negative gradients of the pressure wave can be seen in the figure at 0.0003 and 0.0006 s. At 0.0003 s the gradient inside the MS itself is positive allowing flow reversal to take place, as observed in Figure 7. At 0.0006 s a reflected

13

pressure wave can be seen traveling back toward the outlet. But during this cycle the velocity remains positive although we see it decreasing significantly in the MS between 0.0001 s to 0.001 s. The shape of the profile remains almost similar during this time, with some variation in the profile at 0.0003 s in the MS. 6. Conclusions This study has shown that geometrical changes in the flow path affect the pressure waves and velocity changes generated due to sudden valve closure. Simulation of water hammer in a pipe with initially laminar flow was conducted (𝑅𝑒 ≈ 82) by solving the Navier-Stokes equations using Ansys Fluent software. The pipe had a varying cross sectional area with a rectangular cross section channel (indicated as MS) in the center connected to two pipes, one at each end, using nozzle-diffuser connectors with 5o half angle. The numerical method was initially validated by comparing with results of Holmboe & Rouleau (1967) for water hammer in a straight pipe. This method was then used to simulate flow through the current design. Pressure transients at different locations were reported, including near the valve, center of the MS and upstream of the MS, over the simulated time period. Pressure and axial velocity profiles at various times during the wave cycle were also presented and discussed. The first wave cycle was discussed in detail, including transients corresponding to different peaks in the pressure signal. At different times during the initial cycle, pressure gradient and velocity profiles were presented. These profiles showed how the pressure wave travelled along the domain length and affected the velocity profile in the MS. These results were qualitatively compared with experimental data of Meniconi et al. (2012). It was found that the amplitude of the pressure wave attenuated faster for their experiments where they employed sudden expanding-contracting connectors. But some features of the pressure transients obtained in the current study, especially in the first wave cycle, clearly showed the presence of peaks in the signal that could be related to the pipe geometry. This analysis is useful in designing pipe systems with valves and other components in close proximity that might get affected due to water hammer. 7. Acknowledgements The authors would like to thank Nick Bourgeois (Clearview Subsea) for his help with creating the geometry and Ansys, Inc., for the technical support on the simulations.

14

8. Funding The authors gratefully acknowledge the funding of the study by National Energy Technology Laboratory, Department of Energy. Funding for the project (Project No. 12121-6301-03) is provided through the “Ultra-Deepwater and Unconventional Natural Gas and Other Petroleum Resources Research and Development Program” authorized by the Energy Policy Act of 2005. This program—funded from lease bonuses and royalties paid by industry to produce oil and gas on federal lands—is designed to assess and mitigate risk enhancing the environmental sustainability of oil and gas exploration and production activities. RPSEA is under contract with the U.S. Department of Energy’s National Energy Technology Laboratory to administer three areas of research. RPSEA is a 501(c) (3) nonprofit consortium with more than 120 members, including 16 of the nation's premier research universities, five national laboratories, other major research institutions, large and small energy producers and energy consumers. The mission of RPSEA, headquartered in Sugar Land, Texas, is to provide a stewardship role in ensuring the focused research, development and deployment of safe and environmentally responsible technology that can effectively deliver hydrocarbons from domestic resources to the citizens of the United States. Additional information can be found at www.rpsea.org. 9. References Ansys Inc. (2015). ANSYS® Fluent, Release 16.0, User's Guide. PA 15317, Canonsburg. Anton, B., Simpson, A. R., & Vìtkovsk, J. (2001). Developments in unsteady pipe flow friction modelling. Journal of Hydraulic Research, 249-257. Brunone, B., Ferrante, M., & Cacciamani, M. (2004). Decay of pressure and energy dissipation in laminar transient flow. Journal of Fluids Engineering, 126(6), 928-934. Holmboe, E. L., & Rouleau, W. T. (1967). The Effect of Viscous Shear on Transients in Liquid Lines. ASME Journal of Basic Engineering, 89(1), 174-180. doi:10.1115/1.3609549 Jinping, L., Peng, W., & Jiandong, Y. (2010). CFD Numerical Simulation of Water Hammer in Pipeline Based on the Navier-Stokes Equation. European Conference on Computational Fluid Dynamics ECCOMAS CFD 2010. Lisbon. Meniconi, S., Brunone, B., & Ferrante, M. (2012). Water-hammer pressure waves interaction at crosssection changes. Journal of Fluids and Structures, 33, 44-58. Riasi, A., Nourbakhsh, A., & Raisee, M. (2009). Unsteady Velocity Profiles in Laminar and Turbulent Water Hammer Flows. Journal of Fluids Engineering, 131(12), 121202. doi:10.1115/1.4000557 Zielke, W. (1968). Frequency-dependent friction in transient pipe flow. ASME Journal of Basic Engineering, 109-115.

15

16