Air Pocket Modeling in Water Mains With an Air Valve by Bernardo Carvalho Trindade A thesis submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Master of Science Auburn, Alabama May 7, 2012 Keywords: Water, Air, Water Mains Copyright 2011 by Bernardo Carvalho Trindade Approved by Jos e Goes Vasconcelos Neto, Chair, Professor of Civil Engineering Prabakhar Clement, Professor of Civil Engineering Xing Fang, Professor of Civil Engineering Abstract Water mains lling events follow certain maintenance operations that require the com- plete or partial emptying of pipelines. These operational procedures are performed carefully in order to prevent air pocket formation, as such pockets may trigger damaging surges, as well impose blockages that reduce pipe conveyance. While previous research addressed the modeling of the pipeline lling events, the air phase within the conduits is in most cases either not included in the simulation, or its formulation is very simpli ed, not including air pressure gradients. This paper presents two modeling frameworks which couple a Two-Component Pressure Approach model to simulate water ow with models for the air phase ow within the conduits. The rst is a discretized model for air phase based on the isothermal Euler equations and the second is based on the ideal gas law (as presented in [Zhou et al., 2002a]). Also, the results of an experimental investigation with an apparatus based on the experi- mental setup presented by [Trajkovic et al., 1999] without ideal ventilation which includes the essential features of a pipeline lling event are presented and discussed. Relevant and ow features characteristics during the lling event are analyzed such as the occurrence of interface breakdown. These experimental results as well as eld data of the lling of an ac- tual water main are used to calibrate and assess the model. A feasible modeling framework to simulate the lling of water mains accounting for air e ects was presented. Also, it was demonstrated that lumped approach for air phase modeling has comparable accuracy with discretized model based on Euler equations at a much reduced computational e ort. ii Acknowledgments I would like to thank my advisor and friend, Professor Jose Goes Vasconcelos, for all the patience and time spent in the development of my numerical and experimental research, as well as for all the funny talks about all sorts of non professional topics which made work more pleasant. I also would like to thank Carmen Chosie for the help in the laboratory with the experiments (and with the squeegee) and my dear friend Paulo Ferreira for the help with programming. I would like also to thank my parents, Aluisio (dad) and Regina (mom) Trindade, whose examples were the most important inspiration source for me to move on and want to always achieve greater things in life, both professionally and personally. My thanks to my dear sister, Helena Trindade, for all the funny and sad moments shared up until now, for the constant support during easy and hard times here so far away from home, and for making everything in life more colorful to me. Finally, I am also grateful to all my friends in the civil engineering department who made my time spent in the university very enjoyable. Those are: Gabriel Leite, Zheng Min, Luana Ozelim, Jessica Looper, Tom Hatcher, the nepalian crew (Sushban Shrestha, Nirajan Dhakal and Manoj KC). I would like also to thank Audrey Blackmon for the nice and calming resting breaks and for being with me at late nights while I wrote the the second half of this document. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Di erent types of ows and their classi cation . . . . . . . . . . . . . . . . . 3 1.2 Mathematical models for 1-D water . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Derivation of unsteady fully pressurized ows equations . . . . . . . . 4 1.2.2 Free surface ows equations: the Saint-Venant equations . . . . . . . 10 1.3 Mathematical models for air phase ows . . . . . . . . . . . . . . . . . . . . 12 1.3.1 Euler equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.2 Ideal gas law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4 Traditional solution methods for pressurized pipe model . . . . . . . . . . . . 16 1.4.1 Characteristics form of unsteady equations for closed pipes and its numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4.2 Lumped Inertia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.1 Flow Regime Transition Modeling . . . . . . . . . . . . . . . . . . . . . . . . 21 2.1.1 Interface-tracking Models . . . . . . . . . . . . . . . . . . . . . . . . 22 2.1.2 Shock-Capturing models . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2 Two-phase ow studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 iv 2.2.1 Numerical models and experiments considering the e ects of entrapped air in closed conduits . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.3 Numerical schemes to solve hyperbolic PDEs . . . . . . . . . . . . . . . . . . 43 3 Knowledge gap and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1 Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1.1 Water phase modeling . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.1.2 Water phase source terms . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1.3 Water phase boundary conditions . . . . . . . . . . . . . . . . . . . . 58 4.1.4 Air phase modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Air phase modeling - Uniform Air Pressure Head (UAPH) approach . . . . . 63 4.3 Air phase modeling - Euler equations approach . . . . . . . . . . . . . . . . . 63 4.3.1 Air phase source terms . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.3.2 Air phase boundary conditions . . . . . . . . . . . . . . . . . . . . . . 65 4.4 Experimental program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.4.1 Experimental apparatus setup . . . . . . . . . . . . . . . . . . . . . . 68 4.4.2 Experimental procedure . . . . . . . . . . . . . . . . . . . . . . . . . 69 5 Results and analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.1 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.2 First version of TPAir . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2.1 Assessment of the Euler equation approach . . . . . . . . . . . . . . . 77 5.2.2 Comparison between both approaches . . . . . . . . . . . . . . . . . . 81 5.3 Second version of TPAir . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.3.1 Comparison between experimental results and numerical model pre- dictions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.3.2 Model comparison with actual pipeline lling event . . . . . . . . . . 97 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 v Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 vi List of Figures 1.1 Control volume for continuity equation derivation. Adapted from [Wylie and Streeter, 1993] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Control volume for momentum equation derivation. Adapted from [Wylie and Streeter, 1993] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Control volume for momentum equation derivation. Adapted from [Sturm, 2001] 10 1.4 Control volume for momentum equation derivation. Adapted from [Sturm, 2001] 11 1.5 Example of a typical ow characteristic grid . . . . . . . . . . . . . . . . . . . . 18 1.6 Control volume for lumped inertia models. Adapted from [Zhou et al., 2002a] . 20 2.1 Scheme of shock tting for free-surface ows presented in [Cunge et al., 1980] . 22 2.2 Characteristic grid x-t near the ow regime transition. Adapted from [Song et al., 1983] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3 Bores in the model by [Le on et al., 2010] . . . . . . . . . . . . . . . . . . . . . . 27 2.4 Model conceptualization. Adapted from [Liou and Hunt, 1996] . . . . . . . . . . 28 2.5 Model conceptualization. Adapted from [Li and McCorquodale, 1999] . . . . . . 29 2.6 Preissmann slot representation . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.7 Post shock oscillations (red dotted line). . . . . . . . . . . . . . . . . . . . . . . 31 vii 2.8 Modi ed Preissmann slot geometry. Extracted from [Le on et al., 2009] . . . . . 32 2.9 Bubble and air pockets motion in closed pipe ows. Extracted from [Falvey, 1980] 38 4.1 Flow chart for the model calculation procedures . . . . . . . . . . . . . . . . . . 53 4.2 Model physical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3 Average Courant number for di erent portions of the ow when entrapped air pocket is present. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.4 Pocket nd graphical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.5 Test representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.1 Air phase heads close to the ventilation ori ce for all tested conditions. . . . . . 74 5.2 Air phase heads close to the ventilation ori ce for all tested conditions. . . . . . 75 5.3 Trajectories of moving bore for recirculation ow rates of 2.53 l/s (a) and 5.05 l/s (b) and slope of 2% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.4 Trajectories of moving bore (thick line) and pressurization interface (thin line) for recirculation ow rate of 2.53 l/s and slope of 1% . . . . . . . . . . . . . . . 76 5.5 Model?s rst version layout. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.6 Pressures comparison at the ventilation ori ce for a) initial depth of 0:5D and b) 0:8D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.7 Pressure history at a) x = 7:7m and y0 = 0:5D, b) x = 5:7m and y0 = 0:5D, c) x = 7:7m and y0 = 0:8D, and d) x = 5:7m and y0 = 0:8D . . . . . . . . . . . . 81 5.8 Air velocity at the upstream end of the pipe for tested ow rates . . . . . . . . 82 viii 5.9 Head pro les along the pipe at selected instants . . . . . . . . . . . . . . . . . . 82 5.10 Hypothetical real sized pipeline sketch . . . . . . . . . . . . . . . . . . . . . . . 83 5.11 Comparison of heads at the ventilation ori ce for Q = 0:125 . . . . . . . . . . . 84 5.12 Comparison of heads at the ventilation ori ce for Q = 0:250 . . . . . . . . . . . 85 5.13 Comparison of heads at the ventilation ori ce for Q = 0:500 . . . . . . . . . . . 86 5.14 Air velocities predicted using Euler equations model for all three Q and d orif = 0:15 86 5.15 Pressure history at x = 500m for d orif = 0:05 . . . . . . . . . . . . . . . . . . . 87 5.16 Experimental and predicted pressures at the ventilation ori ce for all tested cases with slope of 0.5% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.17 Experimental and predicted pressures at the ventilation ori ce for all tested cases with slope of 1.0% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.18 Experimental and predicted pressures at the ventilation ori ce for all tested cases with slope of 2.0% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.19 Experimental and predicted pressures at x = 0:39 (from downstream valve) for all tested cases with slope of 0.5% . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.20 Experimental and predicted reservoir heads for all tested cases with slope of 0.5% 94 5.21 Real water main?s sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.22 Field measurements and predicted heads at the upstream ventilation valve of the water main . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.23 Field measurements and predicted heads at the downstream ventilation valve of the water main . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 ix 5.24 Field measurements and ow rates at the upstream ventilation valve of the water main . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 x List of Tables 4.1 Experimental variables. Flow rate was normalized as Q = Q=pgD5 and venti- lation diameter as d = dorif=D where D is the pipe diameter . . . . . . . . . . 70 5.1 Tested values for Euler equation approach assessment (ori ce sizes based on the work by [Zhou et al., 2002a]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.2 Pipeline con gurations parameters . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.3 Long pipe tested values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.4 Summary of results of both models . . . . . . . . . . . . . . . . . . . . . . . . . 96 xi List of Abbreviations a celerity the acoustic waves in the pressurized ow A water ow cross sectional area Aorif ori ce area Aair air ow cross sectional area 4D2 A bvel pressurizing bore velocity Cd discharge coe cient that is assumed as Cd = 0:65 dorif ventilation ori ce diameter d normalized ventilation diameter dorif=D D pipeline diameter E is the total energy per unit volume f friction head loss in the short pipe portion inside the boundary condition right after the inlet Fg Gravitaty force acting on a control volume FShear Shear force acting on a control volume F(U) vector with the ux of conserved variables g acceleration of gravity hc distance between the free surface and the centroid of the ow cross section (limited to D=2) xii hs surcharge head hair extra head due to entrapped air pocket pressurization Hres reservoir water level K elasticity modulus of the material Keq overall local loss coe cient in the inlet Mair mass of air within the pocket Mairout air mass that escapes through the ventilation ori ce MOC Method of characteristics n time step index Pwet water ow wet perimeter Q water ow rate Q normalized ow rate Q=pgD5 Qrec water ow rate which is admitted into the reservoir from the recirculation system Qin water ow rate which enters the upstream end of the pipe S vector of source terms S1 source term for a continuity equation S2 source term for a momentum equation Sdisp;i source term that accounts for air displacement Sf source term that accounts for shear stress between air and pipe walls t time xiii u ow velocity U vector of the conserved variables U vector conserved variables prior to air displacement correction UAPH Uniform air pressure head Vp air pocket volume wdepth local water depth weti boolean variable that is true if the cell is considered wet and false is considered dry bed. Being dry implies that the ow in the cell will not be calculated x x axis direction y y axis direction z z axis direction or pipe elevation relative to the datum 0 shear tension between water and the pipe walls celerity of the acoustic waves in the air x control volume?s length water speci c weight speci c mass of water air speci c mass of air correction factor for air displacement source terms angle formed by free surface ow width and the pipe centerline or between the channel bottom the a horizontal line xiv Chapter 1 Introduction Water mains are important components of the urban infrastructure designed to operate in pressurized ow mode, but under abnormal operational conditions air phase may be present within the conduits. Such conditions are commonly observed during the lling and emptying of water mains, but may also occur following accidental water main ruptures. Those conditions can also happen when water demand exceeds supply, in which case the pipes act as reservoirs. In the process of lling, the air phase present in the pipeline may become entrapped between masses of water, creating air pockets, which should be eliminated by air valves located at selected points along the pipeline. The appearance of these entrapped air pockets in pipelines is well documented in several experimental investigations, such as [Falvey, 1980], [Zhou et al., 2002b], [Vasconcelos et al., 2006]. Entrapped air pockets can cause loss of pipe cross-section, localized head losses, high transient pressures, and other problems if not properly expelled. Yet, even with air valves, the re lling process is often performed cautiously. The reason for this are the potentially damaging surges that may occur due to the air-water interactions during the pipe lling and to ensure that no air pocket remains in the pipe. Because of these uncertainties, the lling is generally performed slowly, and the restoring of normal pipeline conditions may take several hours, or even days, depending on the event that caused disruption of the distribution system. Because this ow interruption has major impacts in water consumers, the ability to quickly and safely restore operational conditions in water mains is highly desirable. Therefore, it is important to understand how to predict the appearance of an air pocket and how to mitigate their e ects. 1 Numerical models represent a feasible way to minimize the uncertainties related to the air-water interactions during the lling process. Three dimensional, CFD packages are a possible solution, but considering some characteristics of the problem, these solutions are not practical. Some of those characteristics are the length of the pipeline systems, the complex geometry of some pipeline elements (e.g. partially closed valves), the ow unsteadiness, and the time scale of the lling events (up to several hours of unsteady ow conditions). One dimensional, unsteady ow models would be an alternative, but they also have limitations. Unsteady ow models in use by water work authorities, usually applied to assess waterhammer e ects and closed pipe surges, assume homogeneous, single-phase ow conditions. By contrast, the re lling process is characterized by two-phase, strati ed ow conditions, invalidating the assumed hypothesis in such single phase models, as demonstrated by [Vasconcelos, 2008]. Therefore, engineers still lack a model which can simulate air pockets in a precise and realistic way. Progress in the area still have not resulted in a tool to help engineers to simulate pipeline priming, particularly considering the e ects of entrapped air. Similarly to other events that involve the transition between pressurized and free surface ows (also referred to as mixed ows), there are certain characteristics on water pipeline lling events that are challenging to numerical models: 1. Even though there are a variety of models which are able to handle ow regime transi- tion, such models have limitations and di culties which include the complexity of the model itself, post-shock oscillations, limitation on the interaction between bores, and others; 2. This is a two-phase (air/water) problem, and models handling the two separate phases need to be appropriately linked. The handling of the interface between air and water is particularly challenging; 3. Due to the formation of bores and the di erences in the celerity magnitudes among di erent portions/phases of the ow, non-linear numerical schemes should be used. 2 This is required to account for both free-surface and water pressurized at the same time limiting issues related to di usion and oscillations at bores and shocks; 4. At certain regions of the ow, particularly at the vicinity of curved air-water interfaces, shallow water assumptions are not applicable due to strong vertical acceleration, and the problem is intrinsically three-dimensional [Benjamin, 1968]; 5. Several di erent mechanisms may result in the entrapment of air pockets during lling events, and conditions leading to such entrapments still not fully understood, thus cannot be numerically predicted. 1.1 Di erent types of ows and their classi cation Water ows can occur in several di erent applications in engineering. However, the characteristics of those ows change from application to application. All of them are strictly three dimensional in space and also varies with time, even if this variation is very gradual. Also, they can have di erent phases, mildly variable speci c mass according to the pressure, can be a free surface or pressurized, and turbulent or laminar. However, it?s very di cult to take all relevant characteristics of water ows into consid- eration when a study is required on a particular case. As one attempts to create a predictive model of a ow in particular conditions, hypothesis are introduced to minimize the com- plexity of the analysis and are related to how the model handles some of the aforementioned characteristics. For instance, regarding the time dependency, ows can be considered steady or transient. The later are more di cult to model, many times requiring the solution of a set of partial di erential equations (PDE) such as the Saint-Venant equations showed in [Cunge et al., 1980], [Sturm, 2001] and [Chaudhry, 2008]. In other applications, when the temporal variation is negligible, ows are assumed to be steady, which eliminates time from analysis, making it much simpler. 3 A ow with phases that do not mix completely is named a multi-phase ow. For those cases, each ow may be analyzed using a di erent set of coupled partial di erential equations. This increases the problem di culty due to the increased number and complexity of those equations that needs to be solved. However, in some cases the properties of both phases can averaged for a certain volume of ow, as in the work presented by [Issa and Kempf, 2003], simplifying the model solution. Another important factor is whether the ow is opened to the atmosphere or not. In the case of rivers, channels, lakes and other ows in that are not con ned, the ow is referred to as free-surface ow. The biggest challenge when modeling this type of ow appears when dealing with the water surface variation in a one, two or three-dimensional transient analysis, so that very sophisticated numerical techniques such as nite volumes or volume of uid (VOF) must be applied. On the other hand, in case of pressurized water pipes, con ned aquifers, drainage pipes, and other ows, the ow is referred to as pressurized ow. Those have known cross section, and techniques to simulate them are more consolidated up to date, such as the Method of Characteristics for closed pipes. There is still the case in which both free-surface and pressurized ow regimes coexist in the same ow, such as in drainage systems, water main lling events, gas/oil transmission lines, and others. 1.2 Mathematical models for 1-D water The model presented in subsection 1.2.1 is used to simulate unsteady ows in closed conduits, while the model presented in subsection 1.2.2 is used for unsteady free-surface ows. The methods for the solution of these equations are presented in section 1.4. 1.2.1 Derivation of unsteady fully pressurized ows equations In this section, the equations for mass and momentum conservation in pressurized ows are presented. Those equations are the basis for the unsteady ows in pipe models. 4 Piezometric line Datum Figure 1.1: Control volume for continuity equation derivation. Adapted from [Wylie and Streeter, 1993] Continuity Equation Figure 1.1 shows a control volume inside a pipe where water would enter from the left. It is possible to have some amount of water stored inside this control volume due to an stretch of the pipe cross-section or length or the to compression of the water mass. Applying the mass continuity the control volume of gure 1.1 yields: Au Au+ @@x( Au) x = @@t A x (1.1) Au@ @x + u@A@x +A@u@x = A@ @t @A@t (1.2) where A is the pipe?s cross-section area, p is the local pressure, x is the control volume?s length, is the water speci c mass, t is the time, u is the ow velocity. Dividing both sides by A yields: 5 1 @ @t + 1 @ @xu+ 1 A @A @t + 1 A @A @xu+ @u @x = 0 (1.3) 1 d dt + 1 A dA dt + @u @x = 0 (1.4) For convenience, the rst term of Equation 1.4, which represent the compressibility e ects, should be rewritten in terms of pressure. For this, the bulk modulus of elasticity concept described in equation 1.5 is used : 1 d dt = 1 K dp dt (1.5) in which K is the elasticity modulus of the material. This concept assumes negligibility of the thermodynamic e ects, which is a good hypothesis due to the low compressibility of water [Vasconcelos, 2007]. For a prismatic pipe, which is the case in most practical applications [Wylie and Streeter, 1993], the gradient @A@x is equal to 0. Hence, the variation of the pipe cross-section area is related only to the variation of the pressure, as in equation 1.6 dA dt = @A @p dp dt (1.6) Replacing the two rst terms in equation 1.4 by equations 1.5 and 1.6: 1 K dp dt + 1 A dA dp dp dt + @u @x = 0 (1.7) 1 K dp dt 1 + KAdAdp + @u@x = 0 (1.8) According to [Wylie and Streeter, 1993], the wave propagation velocity (celerity of the acoustic wave) in a pressurized pipe is given by the expression: 6 a2 = K= 1 + (K=A)(dA=dp) (1.9) which can be written as: K A dA dp = K a2 1 (1.10) and then inserted in equation 1.8 yielding: dp dt 1 a2 + @u @x = 0 (1.11) @p @t 1 a2 + @p @x u a2 + @u @x = 0 (1.12) According to [Wylie and Streeter, 1993], for low Mach numbers (u a) u=a2 0. With this approximation the second term of equation 1.12 can be dropped, yielding: @p @t 1 a2 + @u @x = 0 (1.13) The pressure in the pipe is related to the piezometric head H according to the expression: p = g(H z) (1.14) in which z is the pipe elevation relative to the datum. Applying this expression in equation 1.13, multiplying both sides by a2 and dividing by g, and considering that @z@t = 0 yields: @H @t + a2 g @u @x = 0 (1.15) Note that in case of lateral in/out ow considered in the control volume of gure 1.1, the mass conservation equation 1.15 would show a third term (source term), as shown in [McInnis and Karney, 1995] and [Wylie and Streeter, 1993]: 7 @H @t + a2 g @u @x + a2 gAq = 0 (1.16) in which q is the volume of water gained/lost per length per time. Momentum Equation Piezometric line Datum Figure 1.2: Control volume for momentum equation derivation. Adapted from [Wylie and Streeter, 1993] The application of Newton?s second law to the control volume in gure 1.2, where the summation of all the external forces applied to the control volume is equal to the variation of momentum, yields: pA pA+ @@xpA x +p@A@x x A xsin 0 D x = Adudt x (1.17) in which is the water speci c weight, 0 is the shear tension between water and the pipe walls, D is the pipe diameter, and is the pipe slope. Expanding the derivative of pA, dividing both sides by x and canceling terms yields: 8 A@p@x p@A@x +p@A@x Asin 0 D = Adudt (1.18) which leads to: A@p@x + Asin 0 D = A@u@t (1.19) Knowing that @z@x sin and using the hydrostatic hypothesis represented by equation 1.14, equation 1.19 can be rewritten as: gA @H @x sin gAsin 0 D = A@u@t (1.20) Canceling terms and dividing both sides by A yields: g@H@x 4 0 D @u@t = 0 (1.21) According to [Wylie and Streeter, 1993] the shear stress in pipe ows 0 can be written as: 0 = fu 2 8 (1.22) with f as the friction factor. Writing u2 as ujuj, in order to avoid gain of energy in case of a negative discharge, the momentum equation is obtained: g@H@x + @u@t +fujuj2D = 0 (1.23) With this, the set of PDEs that describe unsteady, pressurized ow in pipes becomes: 8 >>< >>: @H @t + a2 g @u @x = 0 g@H@x + @u@t +fujuj2D = 0 (1.24) 9 1.2.2 Free surface ows equations: the Saint-Venant equations The Saint-Venant equations are also a mathematical model to express mass and momen- tum conservations, but unlike the previous case, the Saint-Venant equations were derived for one-dimensional free-surface ows such as rivers and channels. The hypothesis assumed in the derivation are discussed in [Sturm, 2001], [Vasconcelos, 2007], and others, and are: The ow is unidimensional, resulting in a hydrostatic pressure distribution along the vertical axis of the channel cross-section hydrostatic; The slope is small, thus sin( ) = ; The uid (water) is incompressible, thus is constant; The water depth is small when compared to the wave length; The channel?s bed is stable (e.g. no movable bed conditions); The shear tensions are constant and do not depend on the transient nature of the ow. Mass conservation equation Lateral in ow, qL x Pro le Cross section Figure 1.3: Control volume for momentum equation derivation. Adapted from [Sturm, 2001] Mass balance on gure 1.2.2 yields: 10 Q Q+ @Q@x x t +qL x t = @A@t x t (1.25) Which when simpli ed yields: @A @t + @Q @x = qL (1.26) Momentum conservation Lateral in ow, qL x Centroid Cross sectionPro le Figure 1.4: Control volume for momentum equation derivation. Adapted from [Sturm, 2001] To derive the momentum equation, second Newton?s law is applied to the system in the gure 1.4. The net momentum ux and external forces balance is: @( A x u) @t + @( A xu u) @x = hcA ( hcA+ @( hcA) @x x) +Fg Fshear (1.27) where Fg = A xsin( ) A xS0 (1.28) Fshear = 0Pwet x (1.29) 11 With further mathematical development and combining the equations, one obtains: @( A u) @t + @( hcA) @x + @( Au u) @x = AS0 0Pwet (1.30) where hc is the distance from the water surface to the centroid of the ow cross-section, is the channel slope in degrees, S0 is the channel slope in percentage and Pwet is the wet perimeter. Assuming that Au = Q, is constant, Sf is the friction slope (boundary friction force for each unit weight of water per x), so that Sf = 0Pwet A , the channel is prismatic, and there is no lateral in ow, is canceled in both sides yielding: @Q @t + @ @x Q2 A +ghcA = gA(S0 Sf) (1.31) Therefore, the nal set of equations which are called the Saint-Venant equations is: 8 >>< >>: @A @t + @Q @x = qL @Q @t + @ @x Q2 A +ghcA = gA(S0 Sf) (1.32) which can also be written in conservative format: (U)t + Fx(U) = S(U) (1.33) where U = 2 64A Q 3 75; F(U) = 2 64 Q Q2 A +gAhc 3 75; S(U) = 2 64 0 gA(S0 Sf) 3 75 (1.34) 1.3 Mathematical models for air phase ows This section presents the background for the two alternatives of air-phase modeling implemented in this work. For the discretized alternative, the largely documented Euler 12 equations were chosen. A simpler implementation of the ideal gas law was chosen as the non- discretized air model because of its simplicity, easiness to implement, and low computational requirements. For the air ow in the air pockets, the following hypothesis were considered: Viscosity is negligible: shear forces due to viscosity are negligible when compared to the impelling forces caused by the air-water interaction; The ow is uni-dimensional (1-D): the length of the pocket is much bigger than its lateral dimensions; Isothermal ow: The temperature variation of the air ow during is small. 1.3.1 Euler equations The most complete set of equations for the study of ows is the Navier-Stokes equation, which can be found in [Potter, 2004] and other books. However, they are still very di cult to solve, so that simpli cation hypothesis are of great advantage. With the no shear hypothesis, the viscosity terms of the Navier-Stokes equations disapear and it turns into the Euler equations (equation 1.35), which can be found in [LeVeque, 1992], [Toro, 2009] and other books. (U)t + Fx(U) + Gy(U) + Hz(U) = S(U) (1.35) 13 U = 2 66 66 66 66 66 4 u v w E 3 77 77 77 77 77 5 ; F(U) = 2 66 66 66 66 66 4 u u2 +p uv uw u(E +p) 3 77 77 77 77 77 5 ;G(U) = 2 66 66 66 66 66 4 v uv v2 +p vw v(E +p) 3 77 77 77 77 77 5 ; H(U) = 2 66 66 66 66 66 4 w uw vw w2 +p w(E +p) 3 77 77 77 77 77 5 ; S(U) = 2 66 66 66 66 66 4 S1 S2 S3 S4 S5 3 77 77 77 77 77 5 9 >>> >>>> >>> >>> >>>> >>> >>>> >>> >= >>> >>>> >>> >>>> >>> >>>> >>> >>> >; (1.36) where is the uid speci c mass, u, v and w are the uid velocity in the x, y an z axis, p is the pressure, and E is the total energy per unit volume. The source terms Sn can represent several features of the ow, such as shear stresses and ow area variations, and will be discussed later. The full version of the Euler equations is a three-dimensional, hyperbolic and non-linear conservative set of partial di erential equations, relatively complex and costly to solve. With the application of the second hypothesis (1-D ow), the equations are simpli ed to equation 1.37: (U)t + Fx(U) = S(U) (1.37) where U = 2 66 66 4 u E 3 77 77 5 ; F(U) = 2 66 66 4 u u2 +p u(E +p) 3 77 77 5 ; S(U) = 2 66 66 4 S1 S2 S3 3 77 77 5 (1.38) Further, according to [LeVeque, 1992], for a isothermal problem the energy is constant, making the energy derivatives equal to zero. This leads to equation 1.39: 14 U = 2 64 u 3 75; F(U) = 2 64 u u2 +p 3 75; S(U) = 2 64S1 S2 3 75 (1.39) with p = a2 (1.40) where a is the sound speed in the air. Since, according to [Tran, 2011], the ow condition in pipeline llings may be approximated to isothermal conditions, equation 1.39 is the set of equations used in the discretized model. 1.3.2 Ideal gas law The ideal gas law, which was rst presented in [Kr onig, 1856], states that the state of a mass of a gas depends on its pressure, volume, temperature, the ideal gas constant, and the amount of gas in the studied volume. The relation is: PV = nRT (1.41) where P is the pressure, V is the volume, n is the amount of gas (in moles), R is the ideal gas constant, and T is the temperature. Since for the studied case it was assumed that the event happens under isothermal conditions, the right side of equation 1.41 becomes constant: PiVi = PfVf (1.42) where the the subindexes i and f denotes initial and nal values. This mathematical model is usually supplemented with an ori ce-type equation to relate air pressure with discharge in presence of ventilation ori ces. 15 1.4 Traditional solution methods for pressurized pipe model In this section, two methods for solving the fully pressurized pipe model derived in section 1.2.1 are presented. The rst one is the method of characteristics (also referred to as MOC), which consists in nding two ODEs for head ow rate that are valid along characteristic lines. The second one assumes certain simpli cations which simplify the system of PDEs showed in equation 1.24 to a single ODE, making it easier to implement and compute, but with more restrict applicability. 1.4.1 Characteristics form of unsteady equations for closed pipes and its nu- merical solution A popular way to model transient ows in water mains is the method of characteristics. Among the advantages of this method are the simple implementation and computational e ciency. The limitation is the fact that it assumes that the ow is homogeneous and single phase, which is not true for pipe lling problems. The idea of the method is to apply a characteristic transformation on the mass and momentum equations for pressurized ows and integrate them. Therefore, the goal is to nd a constant that renders the linear combination of the PDE system as an ODE system, as described below: 8 >>< >>: @H @t + a2 g @u @x + a2 g q A = L1 @u @t +g @H @x +f ujuj 2D = L2 (1.43) L2 + L1 = 0 (1.44) or @u @t +g @H @x +f ujuj 2D + @H @t + a2 g @u @x + a2 g q A = 0 (1.45) 16 As it can be seen in [Wylie and Streeter, 1993], the that satis es equation 1.45 is: = ga (1.46) Substituting in equations 1.44 by the expression in equation 1.46 and multiplying everything by A yields the following characteristic form of mass and momentum conservation equations: 8 >>> >>>> >>< >>>> >>> >>: @Q @t +A g a @H @t +f QjQj 2DA = 0 dx dt = a @Q @t A g a @H @t +f QjQj 2DA = 0 dx dt = a (1.47) The next step towards a numerical solution of equation 1.24 is to integrate this set of equations in order to nd analytic expressions for both pairs of equations (positive and negative signs of a). The steps of the integration can be seen in [Wylie and Streeter, 1993] and result in: 8 >< >: HP = HA B(QP QA) RjQAjQP HP = HB +B(QP QB) +RjQBjQP (1.48) Some constants are de ned in order to simplify the previous equations: CP = Hn 1i 1 +BQn 1i 1 CM = Hn 1i+1 BQn 1i+1 BP = B +RjQn 1i 1j BM = B +RjQn 1i+1j (1.49) 17 t x x t Figure 1.5: Example of a typical ow characteristic grid With this, equations 1.48 turn to be: 8 >< >: Hi = CP BPQi Hi = CM +BMQi (1.50) which are valid for the intersection between the characteristic lines, as in gure 1.5: dx dt = a (1.51) Equations 1.50 are the nal version of the numerical integration of characteristics form of mass and momentum equation. This form is the one present computational codes, such as Hammer [Bentley, 2010] and Pipe2011 [KYPipe, 2010]. 1.4.2 Lumped Inertia Lumped inertial model is a class of models with assumes that the whole water mass can be treated as a rigid column. This simpli cation is made in order to transform equations 1.24 into a single ordinary di erential equation. The hypothesis for the lumped inertia approach are: The water phase is incompressible (d dt = 0); 18 The pipe walls are perfectly rigid (this combined with the rst hypothesis implies that dH dx = const); Considering those hypothesis, the calculation for the water phase can be reduced from a system of two PDEs (equation 1.24) into one ODE. Considering this control volume and applying the hypothesis? to equation 1.23 with localized head loses, it follows that: g@H@x + @u@t + f +KeqDL ujuj 2D = 0 (1.52) Isolating @u@t yields: @u @t = g @H @x f ujuj 2D Keq u2 2L (1.53) Knowing that dHdx is constant over the pipe length (second hypothesis), the nal form of the equation becomes: @u @t = H H0 L f ujuj 2D Keq u2 2L (1.54) According to [Vasconcelos, 2007], this equation can be solved analytically, with the solution as follows: u(t) = r 2g(Hu Hd) R tanh t 2L p 2g(Hu Hd)R (1.55) in which R = fLD + 1 and the subindexes u and d mean upstream and downstream. Another option is to solve the equation numerically for each point using Runge-Kutta or Forward Euler methods [Press, 1989]. [Martin, 1976], [Liou and Hunt, 1996], [Zhou et al., 2002a] and [Fuertes et al., 2000] applied lumped inertia models to mixed ows in pipes. Their models in addition to using the hypothesis of the lumped inertia, assumed that the contact interface between air and water 19 is well-de ned and perpendicular to the pipe walls. This second hypothesis is illustrated in gure 1.6. Well-de ned interfaceWater Air Figure 1.6: Control volume for lumped inertia models. Adapted from [Zhou et al., 2002a] 20 Chapter 2 Literature Review This section presents a summary of the current body of knowledge on ow regime transition modeling, two phase ows, and some numerical schemes. It starts by presenting interface-tracking and shock-capturing models for ow regime transition with perfect venti- lation, followed later by studies in which the e ect air alongside water ow is accounted for. This chapter complements the theoretical basis of the present work, showing more recent ideas that were applied in the model proposed in this dissertation, and similar alternatives that have been proposed along the past few decades. 2.1 Flow Regime Transition Modeling Flow regime transition is a type of liquid ow in which there is the transition from free-surface to pressurized ow and vice-versa, so that the existence of both regimes need to be taken into account by the model. Some of the inherently di culties with this problem are: No single set of equations is able to account for both regimes; Wave celerities may be up to three orders of magnitude di erent between two ow regimes; Di culties on how to handle the interface between both regimes. In this section, two di erent modeling approaches for handling ow regime transition are presented. One is the interface tracking approach, which tracks the position of the interface between both phases and divides the ow into two di erent domains (pressurized 21 (n+ 1) t n t C B1,B2 B A xA xB A1 A2 B1 B2 (n+ 1) tn t x t t x S R Figure 2.1: Scheme of shock tting for free-surface ows presented in [Cunge et al., 1980] and the free-surface). The other approach is the shock-capturing, which uses a single set of equations to calculate both portions by introducint a conceptual model, eliminating the need of tracking the interface position. 2.1.1 Interface-tracking Models Basic concept of Shock-Fitting A general concept behind interface-tracking is the Shock-Fitting technique, described in [Cunge et al., 1980]. In this approach, a discontinuity in the solution in the form of a hydraulic jump or a pipe- lling front (pressurizing bore) is tracked explicitly over time. This tracking represent the calculation of the discontinuity position and velocity and other relevant ow variables. Those variables are both depths and discharges (up and downstream), the position and the velocity of the discontinuity. In the shock tting method presented in [Cunge et al., 1980] for river ow, the bore motion is solved by means of the Method of Characteristics applied to the Saint-Venant equations, as presented in gure 2.1. 22 Therefore, the system of equations to solve the ow in both sides of the discontinuity are the characteristic equations for both sides of the discontinuity. The usage of the method of characteristics adds three more unknowns for the problem, which are the origin points of the characteristic lines upstream and downstream the bore, as shown in gure 2.1. Note that the supercritical side provides only one characteristic line. Also, mass and momentum equations across the bore are necessary to solve the problem. Hence, according to [Cunge et al., 1980], the four characteristic equations for the sub- critical (downstream) side are: xB2 xS = (tB2 tL) (Q=A) B2 +cB2 2 + (Q=A)S +cS 2 (2.1) (Q=A)B2 + 2cB2 = (Q=A)S + 2cS +g(tB2 tS) S 0B2 +SfB2 2 + S0S +SfS 2 (2.2) xB2 xR = (tB2 tR) (Q=A) B2 cB2 2 + (Q=A)R cR 2 (2.3) (Q=A)B2 2cB2 = (Q=A)R 2cR +g(tB2 tL) S 0B2 +SfB2 2 + S0R +SfR 2 (2.4) where the subindexes represent the positions in gure 2.1. Also, according to [Cunge et al., 1980], the characteristic equations for the supercritical side (upstream) are: xB1 xL = (tB1 tL) (Q=A) B1 +cB1 2 + (Q=A)L +cL 2 (2.5) 23 (Q=A)B1 + 2cB1 = (Q=A)L + 2cL +g(tB1 tL) S 0B1 +SfB1 2 + S0L +SfL 2 (2.6) Another two equations are the ones required to nd the velocity of the shock. S = AB1(uB1 uB2)A B1 AB2 +uB2 (2.7) uB1 uB2 = gAB1 AB2A B1AB2 (AB1hcB1 AB2hcB2) (2.8) The last equations and also the simplest one, it the di erential expression of the bore front velocity: dx dt = S (2.9) Also, according to gure 2.1, xB1 = xB2 = xB The system behind the track of bores has nine equations (equations 2.1 to 2.9) for nine unknowns (xL, xS, xR, x, S, QA, AA, QB and AB). It is also important to notice, as shown in gure 2.1, that the interface between the two regimes (the bore) is considered to be perfectly de ned and vertical. Models based on Shock-Fitting theory This idea of applying shock tting to simulate ow regime transitions was rst used by [Wiggert, 1972]. This model solves the free surface ow using the method of characteristics (similar to the development in 1.4.1) for equation 1.34 and calculates the front parameters as presented above. On the other hand, a lumped inertia approach as the discussed in subsection 1.4.2 for the pressurized portion of the ow was used. To validate the model, it was used a relatively long channel (30 meters long) with rectangular cross-section and an intermediate portion partially covered forming a rectangular pipe, having a still water level as initial 24 Figure 2.2: Characteristic grid x-t near the ow regime transition. Adapted from [Song et al., 1983] conditon. As boundary conditions the apparatus had a downstream weir, and a upstream reservoir with a gate which would be opened to generate a ow that would pressurize the channel. The in ow front attained its maximum speed in a short amount of time after the gate opening, decelerating later along the channel and then slowly reaccelerating again when close to the weir and becoming steeper along the whole path. The main goal of those experiments was to calibrate and validate the proposed model. [Song et al., 1983] proposed a model that would use the method of characteristics also for both portions of the ow. The practical problem of this alternative is the discrepancy of the time steps between both portions. In order to overcome this problem, the model calculates several time steps of the pressurized portion for each time step of the free surface portion. The resulting characteristic grid for the vicinity of the bore would be as show in gure 2.2. Other examples of models that use MOC are [Fuamba, 2002], [G omez and Achiaga, 2001]. Even though those models can simulate relatively complex interactions between the free surface ow and the pressurizing bore, they have certain limitations. One of them comes from the MOC for the free surface ow, which incurs in a inability to simulate non-pressurizing bores, which are common in practical applications. Another limitation is the handling of the interface between free surface and pressurized ow, which is considered to be always 25 in the form of a bore. The model presented by [Politano et al., 2007] addressed the last limitation allowing gradual ow regime transition, calculating mass and momentum balance for a control volume comprising the last pressurized and the rst free-surface nodes. A model that partially overcomes some of those limitations is presented in [Le on et al., 2008]. The authors used the nite volume non-linear scheme (Riemann solver) HLL to a circular geometry developed in [Le on et al., 2006] in order to solve the Saint-Venant equations for the free surface portion of the ow in a shock capturing fashion, which makes the predictions of non pressurizing bores very precise without needing to track them ( gure 2.3c). For the pressurized portion of the ow a compressible water hammer formulation which is able to handle distributed cavitation was used, which is also solved with the HLL solver. The model agrees well with experimental data for systems with small ventilation, however, some of the inconveniences of the other shock- tting approaches still persists. One of them is the complexity of the model, which in this case comes from the calculations of the air water interface and also from the variation of the HLL scheme that was used. Also, a common di culty for all the Shock-Fitting models is the tracking of multi- ple pressurizing fronts, which happens often in practical applications, and the interactions between them, and other ow features. Other interface tracking models The interface-tracking model present in [Liou and Hunt, 1996] proposes a solution for the ow regime transition problem in a water main lling event without going through the complexity of the shock- tting theory. This model use the lumped inertia approach assuming a vertical interface between air and water and a very rapid ow so that as the water reaches a pipe cross-section this cross-section goes instantaneously from wet bed to fully pressurized without the occurrence of a free-surface ow, which makes this model not to be formally a ow regime transition model, but still dealing with ow regime transition problems. This model solves equation 1.54 for each pipe, as shown in gure 2.4 with a Runge-Kutta method 26 Figure 2.3: Bores in the model by [Le on et al., 2010] with adaptive step size, showed in [Press, 1989]. One important limitation of the model is that air intrusion is not considered, otherwise free-surface ow would occur and the lumped inertia hypothesis become invalid. The model was used to simulate the lling of a theoretical water main with variable slope and of a laboratory experimental apparatus, which had good agreement with the model. Also, the experimental data was used to validate the model only at the early stages of the ow, when the velocity of the front is still very high, not representative of the typical velocities at the rest of the event. Another lumped inertia model was proposed by [Li and McCorquodale, 1999]. This model considers the water portion of each ow regime as a rigid column, calculated with the lumped inertia approach, as seen in gure 2.5. The two rigid columns (pressurized and free-surface) are simulated with data from the moving bore by using the following continuity and momentum relations: (V1 +Vbore)A1 = (V2 +Vbore)A2 (2.10) 27 Figure 2.4: Model conceptualization. Adapted from [Liou and Hunt, 1996] and w(y1A1) +pA2 ( wy2A2 + wZA2) = wA1(V1 +Vbore)(V2 V1) (2.11) where w is the speci c weight of water and Z is the pressure head on the surcharge side of the hydraulic jump. Air e ects are also considered in this model as it will be explained later, as this approach have similarities with the model proposed in this work. [Malekpour et al., 2011] simulated the rapid lling of water mains with the method of Characteristics for a full dynamic model showing good agreement with experimental data. However, the applicability of the approach is limited considering that the rapid lling of water mains is not desired this may damage the pipe. As it was shown, currently there are di erent options of interface-tracking models for ow regime transition with varying degrees of complexity and able of simulating relevant ow features. However, this class of models still needs to overcome some limitations, such as the di culty in simulating ows in more complex systems with several pressurizing bores and depression waves which may interact with each other. Another drawback of such models is the complexity of the calculations for the interface between free-surface and pressurized ows that exceeds the alternative provided by shock capturing models. 28 Figure 2.5: Model conceptualization. Adapted from [Li and McCorquodale, 1999] 2.1.2 Shock-Capturing models Another family of ow regime transition models is called shock-capturing. The name is due to the fact that those models don?t need to explicitly track a moving bore. This is accomplished by using a single set of equations and an appropriate choice of conserved variables (e.q. Q and A instead of y and u) for both ow regimes, which eliminates the need of changing the way the ow is calculated before upstream and downstream the pressurizing bore. Preissmann Slot models The rst Shock-Capturing model handling ow regime transition was presented by [Cunge and Wegner, 1964] using an idea presented in [Preissmann, 1961] named the Preiss- mann Slot. This model assumes the existence of a ctitious slot on the crown of the pipe which simulates a pressurized ow as a free surface ow with the water lling the slot, as shown in gure 2.6, so that it can be calculated with the Saint-Venant equations (equation 29 hslot tslot D Figure 2.6: Preissmann slot representation 1.34). In gure 2.6 D represents the pipe diameter, tf the width of the slot and hslot the water level inside the slot (from the crown of the pipe until the surface). The key of this strategy is to set two di erent ways to calculate the static momentum of the geometric cross-section (A hc) in the Saint-Venant momentum equation, one for free-surface ow and the other for a pressurized ow. For the free-surface ow this static momentum for a circular cross-section is calculated using the formula shown in [Akan, 2006]: 8 >< >: D3 24 [3sin( ) sin 3( ) 3 cos( )] = arccos[(y D=2)(D=2)] (2.12) where D is the pipe diameter, is the angle formed by free surface ow width and the pipe centerline, and y is the water level (which is smaller than D), while for the pressurized portion of the ow the static momentum is calculated by: Ahc = tsloth 2 slot 2 + D2 4 hslot + D2 (2.13) To consider the desired celerity of the pipe when a pressurized ow is present in the simulations the slot width must be properly chosen. For this, equation 2.14 is used: 30 Pressurized owFree-surface ow withentraped air pocket Figure 2.7: Post shock oscillations (red dotted line). tslot = gA 2 pipe a2 (2.14) Because of the use of a single set of equations to solve both ow regimes, this model is signi cantly simpler than interface-tracking models. However, the Preissmann slot model has two disadvantages. The rst one is shown in [Cunge and Wegner, 1964] and is about negative pressure ows. In the Preissmann slot model, if the pressure, which is seen as the water level, drops bellow the crown of the pipe for any portion inside the pressurized ow zone the free surface ow will be restored for this low pressure zone, which means an unrealistic assumption of perfect ventilation along the whole pipeline. The second limitation is regarding spurious post-shock oscillations, which will always happen in the pressurized portion of the ow right next to the pressurizing bore, as shown in gure 2.7. Those oscillations can lead to peaks which are more than one order of magnitude higher than the actual pressure, distorting the results and making a computational code crash. These oscillations happen due to the di erence in celerity between both portions of the ow, and are discussed in more details by [Vasconcelos et al., 2009a]. To overcome this limitation, a common strategy is to use an arti cially low celerity values in between 50 to 200 m/s (as in [Trajkovic et al., 1999] and [Vasconcelos et al., 2006]) instead of using the real pipe celerity given by equation 1.9. 31 Figure 2.8: Modi ed Preissmann slot geometry. Extracted from [Le on et al., 2009] However, when very low celerity are adopted (e.q less than 50 m/s) this can lead to other problems. When high pressures happen in the ow, signi cant continuity errors can happen due to the accumulation of water in the slot, which is not physical. With this, arti cially low pressure peaks [Vasconcelos et al., 2006] and wave propagation velocity can be anticipated by the model. In order to avoid the use of a low celerity for the whole pipe, [Sj oberg and CTH., 1982] proposed a modi cation in Preissmann slot model in which the width of the slot varies with the ow head and the pipe diameter according to equation 2.15, so that it would avoid a sudden transition between a low celerity to a high celerity in the moment the ow touches the crown of the pipe. Ts = D 10 6 + 0:05423e (h=D)24 when h 0:9999D (2.15) The model presented by [Le on et al., 2009] applied a similar idea but instead of using a equation such as equation 2.15 a criteria based on the head and diameter with linear approximation between each combination of head and a percentage of the diameter was used, as shown in gure 2.8. 32 Another improvement in the Preissmann slot was made by [Kerger et al., 2011]. This model considers a slot with negative depth for a pressure drop to a negative value in the pressurized portion of the ow instead of restoring a free-surface ow. However, this model still su ers from the limitation of the arti cially low pipe celerity to avoid the post shock oscillations. Also, the hypothesis of the negative depth of the slot is not at rst intuitive, making it di cult to understand its physical behavior and formulate other hypothesis for improvements in the model. Other works which use the Preissmann slot model are [Capart et al., 1997], and [Tra- jkovic et al., 1999]. The latter performed experiments to calibrate and successfully validate the model with an apparatus consisting of a 10 m long pipe with 10 cm of inner diameter, eight ventilation pipes to eliminate the e ects of air compression, a upstream reservoir two gates being one in the end and the other 1.5 m away from the tank. This experimental set up was adapted for the experiments performed in this dissertation. A detailed explanation about the Preissmann Slot can be found in [Cunge et al., 1980]. Two-component Pressure Approach (TPA) In order to face the limitation of negative pressure of the Preissmann Slot model and the di culty of handling a large number of pressurizing bores and their interactions with the shock- tting models, [Vasconcelos et al., 2006] proposed the TPA model. This model makes use of the similarities between equations 1.24 and 1.34, with hypothesis on how the celerity is adjusted for both ow regimes add on a surcharge term enabling both coexisting ow regimes to be modeled with a single set of equations. As shown in [Vasconcelos, 2005], assuming sin S0 and fDujuj2g = Sf, equation 1.24 can be written with some manipulation as: @A @t + @Q @x = 0 @Q @t + @ @x Q2=A+gAh c = gA(S 0 Sf) (2.16) 33 If rewritten in quasi-linear form by adding ghAx to both sides of the equation, the above equation becomes: Ut + ~AUx = S ~A= 2 64 0 1 Q2A2 + A @P@A 2QA 3 75 (2.17) Equation 1.34 (Saint Venant equations) expressed in terms of total pressure and written in quasi-linear form yields exactly the same expression as equation 2.17. With this similar- ity observed, the only di erence between both equations in the interpretation of the term (A= )(@p=@A). For open channel ows, this term is described as the square of celerity, as shown in equation 2.18 A(h) = Z h 0 Ts(y)dy @A(h) = Ts(h)@h , and couped with @p = g@h A @p @A = g A(h) Ts(h) = c 2 (2.18) where c is the ow celerity, y is the depth for the integration and h is the actual ow depth. For a incompressible pressurized ow, assuming a elastic behavior for the pipe walls, (A= )(@P=@A) represents the acoustic wave speed, which is the equivalent of the celerity of the free surface ow for pressurized ows. With this, one has: A @P @A = c 2 = a2 = P A A + with = 0 ! a2 = A P A (2.19) Considering that P = g h as this P is the extra pressure due to the ow pressur- ization, equation 2.19 can be modi ed so that it yields: hs = a 2 g A Apipe (2.20) 34 where hs is the extra head due to the pressurization and A is the ow area the exceeded the original cross-section area of the pipe, making this original cross-section stretched or contracted. This term is then inserted in the Saint-Venant equations (equation 1.34) as a head that is added to the free-surface head which goes beyond the centroid of the cross- section in pressurized ows (hc). This yields the nal format of the Saint Venant equations modi ed by the TPA model: U = 2 64A Q 3 75; F(U) = 2 64 Q Q2 A +gA(hc +hs) 3 75; S(U) = 2 64 0 gA(S0 Sf) 3 75 (2.21) hs = 8 >>< >>: 0 !free-surface ow a2 g A Apipe !pressurized ow (2.22) hc = 8 >>> >>> < >>> >>> : D 3 3sin( ) sin3( ) 3 cos( ) 2 sin(2 ) !free-surface ow where = arccos[(y D=2)(D=2)] D 2 !pressurized ow (2.23) where = arccos[(y D=2)(D=2)]. In the case where the ow is pressurized and the pressure is negative A in equation 2.21 becomes negative instead of regenerating the free- surface ow. This gives a negative value of hs. Equation 2.21 is also reached if a surcharge pressure is considered in the derivation of equation 1.34. Unlike Preissmann slot, hc is limited by the value D=2. With this, the TPA model overcame the inability of the conventional Preissmann slot model to simulate pressurized ows with negative pressure. Since its introduction, other shock-capturing models have also overcome the limitation of the traditional Preissmann slot model. A work by [Bourdarias and Gerbi, 2007] describes a model which is very similar 35 to the TPA. However, as described in [Vasconcelos et al., 2009a], all those models shock- capturing models based on the Saint-Venant equations still su er from the limitation of the occurrence of post-shock oscillations, leading to errors due to a arti cially low celerity as described for the Preissmann slot models. Alternatives to mitigate those oscillations based on numerical ltering and a new ux function called hybrid ux are described in [Vasconcelos et al., 2009a]. 2.2 Two-phase ow studies Air in water mains and slug ows are a concern for engineers that has been studied both numerically and experimentally by several researchers. Air in water mains is characterized by entrapped air, normally during the lling process, which engineers want to eliminate by means of air valves, drag forces caused by water, among others. Another related subject, according to [Fabre and Lin e, 1992], are slug ows, which are a ow pattern with sequences of long air pockets almost lling the pipe, followed by liquid slugs which may contain small bubbles. Those ows are observed in systems with steady injections of liquid/gas phase and occur for a range of ow rates in gas/water ows, as described in [Falvey, 1980]. An important di erence between slug ows and problems on air present in water main is that in the slug ow the gas portion of the ow is forced into the conduits and are transported along the pipe with the liquid phase, such as in oil/natural gas pipelines. With water mains, air is initially present and is gradually expelled by means of lling. Regarding the numerical approach to solve each case, a feasible way of trying to model the problem of air and water ow in pipes is by coupling ow regime transition models presented in the last section with air phase models whereas for slug ow a statistical approach is one alternative to describe ow characteristics on both phases. As presented in this subsection, several researchers have conducted experiments about entrapped air in water ows. While some of those experiment were performed with ex- perimental focus, several others were performed to assess the application and/or calibrate 36 numerical models. Many of those works about entrapped air were done focusing on the dragging of air pockets and bubbles. However, since this work focus on ventilation of air pockets those other works will not be detailed here. 2.2.1 Numerical models and experiments considering the e ects of entrapped air in closed conduits Studies about entrapped air in water pipe lines started as early as [Kalinske and Bliss, 1943], when the authors performed experiments on dragging of air pockets by air ows and developed equation 2.24 to determine the air ow rate given a water ow rate based on the experimental data. The proposed expression is valid as long as the water ow rate is high enough drag forces remove the air, even when buoyancy forces oppose the air motion. The work presented [Pothof and Clemens, 2010] and [Pozos et al., 2010] tackled the same problem, presenting expressions that relate air discharge to water discharge obtained by correlation of experimental data and conservation laws. Q a Qw max = 0:0066(Fr 1)1:4 (2.24) Until the beginning of the last decade, most one-dimensional two-phase ow models were of the interface-tracking type. Possibly, the rst such work was [Martin, 1976], which presented a lumped inertia model, assuming a well de ned interface between air and water that advances over the dry portion of the pipe expelling the air through a ori ce modeled with a simple head discharge relation. However, the well de ned interface hypothesis is unrealistic for water main lling events. An important technical report presented by [Falvey, 1980] has a chapter focusing on air- water interactions in closed pipe ows. In the rst section on pipe ows, a comprehensive compilation of previous works is presented, including results, charts, and design criteria for air valves. An air-water ow classi cation is presented in terms of the air/water interaction pattern (small bubbles, big pockets, etc) based on the relative ow rates of both phases. After 37 Figure 2.9: Bubble and air pockets motion in closed pipe ows. Extracted from [Falvey, 1980] presenting air/water ow classi cation, a discussion is presented on ows in partially lled conduits, which is basically an open channel ow inside a pipe with a moving air layer on the top. The focus of that section is to present several models that try to predict the amount of air that would be transported along the pipe from upstream driven by the drag forces in between the two phases. The third section of that chapter presents a discussion about ows having a hydraulic jump that lls the conduit, which includes a review on several works that established empirical relations for the air pockets and bubbles drag, such as [Kalinske and Bliss, 1943] and [WISNER and Bucarest, 1967]. A chart presents the conditions anticipated for air pockets and bubbles to be dragged (hydraulic removal of air) according to pipe geometric characteristics and water ow rate. Air ventilation in conduits is also discussed in the report, with suggested criteria for location of ventilations valves, precautions against freeze, cavitation, and harm to personal in the vicinity. In the sixth and seventh sections, design criteria for the ventilation in pipelines and pump systems are presented, along with several types of ventilation valves, charts and head vs: discharge relations are presented. 38 The work presented by [Hamam and McCorquodale, 1982] aimed to study mechanisms for air pocket formation using a apparatus with rectangular and circular cross-section pipes. The studied mechanism was the formation of waves in the water surface caused by the relative air-water motion that eventually touch the crown of the pipe. This requires a certain amount of shear stress between air and water, which was accomplished by blocking downstream the free surface ow in the pipe. This blockage creates a pressurizing bore that the air ahead of it with a velocity in the opposite direction of the water ow, creating a high relative velocity between phases. With their results they concluded that the ow regime transition associated with those pockets can cause signi cant pressure peaks associated with the air expelling. The model presented in [McCorquodale and Hamam, 1983] is based in the lumped inertia approach for the water phase a simple compressible ow theory for the air phase. The model presented in [Li and McCorquodale, 1999] is a re nement of [Hamam and McCorquodale, 1982] in the sense that the rst one considers the movement of the air pockets. A similar model was presented in [Zhou et al., 2002a], in the sense that this model also used a lumped inertia approach with a moving, well de ned vertical interface which expels the air ahead of it. A di erence from the model presented in [Martin, 1976] is the ori ce formulation, which considers a chocked ori ce. The focus of this work was to measure maximum surges under di erent in ow/ventilation con gurations. The model was compared to experimental results which showed good agreement except for the oscillation pattern for substantial air release. Regarding the experimental part, the apparatus consisted of a reservoir linked to a horizontal pipe with a control valve separating the initially lled and empty portions of the pipe; and interchangeable nozzles with di erent diameters providing ventilation on the other end of the pipe. Three typical outcomes were noticed: negligible water hammer e ect, in which the air pocket "absorbs" the water hammer e ect; mitigated water hammer e ect, in which a the water hammer e ect is again absorbed but a high pressure peak is observed and the end of the lling event when the water hits the ventilation 39 ori ce; and water hammer dominant, in which the air content is rapidly expelled and the water slams the ventilation ori ce generating a pressure peak almost four times higher as the mitigated one. They found that the ratio between the ventilation ori ce and pipe diameter is determinant to the intensity of the transient pressures and oscillations pattern. However, the very rapid lling nature of the experimental conditions as well as the extreme high pressures when compared to the pipe diameter do not represent actual water main lling conditions. [De Martino et al., 2008] also studied those pressure peaks due to the expel of air pockets through ventilation ori ces in a purely experimental fashion, deriving an expression for the maximum expected pressure peak. Other models dealing with lumped inertia-type approach for water phase are [Kabiri Samani et al., 2006] and [Izquierdo et al., 1999], this last one proposing a model for the lling of whole sloping water pipe lines with air pockets. The model was tested with experimental data in [Fuertes et al., 2000] in order to validate the model again under extreme conditions. [Chaiko and Brinckman, 2002] presented a comparative study of three models to simulate air/water interactions in a pipe lling problem. The rst model presented was a full dynamic approach with both phases being solved by the method of characteristics, so that water characteristic lines need to be interpolated to match the grid; the second was the same full dynamic for water and lumped for air; and the third applied lumped inertia approach for air MOC for water phase, however calculating only the unperturbed portion of the water ow so that the characteristic lines have constant slope and match the grid without the need of interpolation. The author runs tests for a vertical set up which consisted of a a cylinder with a air pocket on the upper part which is compressed by the water phase due to a increase in the water pressure at the bottom of the cylinder. The authors showed that the second model (full dynamic for water and lumped for air) capture all the relevant events as well as the rst model, even though the second didn?t capture small oscillations due to the re ection of the pressure wave in the air, which has no practical importance. However, the problem that 40 was proposed by the authors is too idealized, given the experimental set up and given the well de ned interface between phases. [Zhou et al., 2011] developed a model based on the method characteristics for the water phase and a variation of the ideal gas law with a polytropic coe cient of 1.4 instead of 1.0 for the air phase, with simple mass and momentum conservations in the interface between phases for the interface tracking. Also, a series of experiments were performed to show the e ect of small air pockets (0% to 8.02% of void fraction in the pipe) in the pressure peak and validate the model. The authors concluded that instead of having "the maximum peak pressure of air pocket increases with the decrease of the initial void fraction of air pocket" as it was stated in their previous works (such as [Zhou et al., 2002a]), the maximum peak happened for a void fraction of 6.18%. A recent work with the interface-tracking approach is presented in [Le on et al., 2010] as a enhancement of the model proposed in [Le on et al., 2008]. The model presented in the work handles the water phase using the HLL Riemann solver for the Saint-Venant equations (free surface ow) and for the compressible water hammer formulation presented rst in [Guinot, 2003] (pressurized ow). For the air phase the approach was the same presented in [Martin, 1976] and [Zhou et al., 2002a]. Two conditions were considered for the air phase: with and without air release. For both of them, the model matched well with the experimental data from [Vasconcelos et al., 2006] and [Zhou, 2000]. On the shock-capturing front, [Arai and Yamamoto, 2003] developed a model based on the Preissmann slot which assume the existence of a cap over the slot to avoid ventilation. A modi ed version of the Saint-Venant equations was used for the water phase as well as a structurally similar version of those equations were used for the air-phase. This particular set of equations was used in order to apply the four-point Preissmann implicit scheme presented in [Cunge et al., 1980] to solve the equations. In order to calibrate a model, their experimental studies were made with a scale model of an underground drainage system with 122.08 meters of pipe, being one of the longest reported in literature. With their experiments they showed 41 that the existence of entrapped air a ects the time for the ow regime transition to happen, suppresses inertial oscillation, and rises the maximum observed pressure in the system, also increasing the pressure variation in the inlet. Their two-phase full dynamic ow model showed good agreement with the experimental data. Another shock capturing model was developed by [Vasconcelos and Wright, 2009] as a enhancement of the TPA model presented in [Vasconcelos et al., 2006]. In this model, a term to account for air pressure is added to the momentum equation of the TPA model, changing the vector F in equation 2.21 for the one in equation 2.25. F(U) = 2 64 Q Q2 A +gA(hc +hs) +gApipehair 3 75 (2.25) where hair is the term that accounts for the air pressure, making part of the link between both phases. While the air ow rate through the ventilation ori ce is given by the expression presented in [Zhou et al., 2002a], the air head is calculated by the continuity of the air phase, which leads to the following expression: hair = 12g a w u S(Apipe Afs) CdAorifY (2.26) where uS is the velocity of the moving bore, Afs is the area of the free-surface ow, and Y is an expansion factor showed in [Zhou et al., 2002a]. The work by [Tran, 2011] presented a lumped model for bubbly ows which accounts for the e ects of liquid compressibility, pipe elasticity and temperature rise across a pressure wave. The author found that below a certain air content the e ects of pipe elasticity and liquid compressibility are signi cant, while for a bubbly ow with a reasonably high amount of air (above 10% or 20%) the transient ow may be considered isothermal in part due to the high thermal capacity of water. 42 In [Fabre and Lin e, 1992] the authors presented a model for slug ow. The problem of long pockets as well as small bubbles generated by slug ow is dealt with in a statistical fashion given the virtually random behavior and intermittency of the details of this type of ow. Also on slug ows, the model presented in [Issa and Kempf, 2003] is a shock-capturing technique coupled with a formulation that accounts for the Kelvin-Helmholtz instabilities which form a slug ow, as described in [Falvey, 1980]. However, the focus of this work was to only predict the frequency and the occurrence of slug ow due to the Kelvin-Helmholtz instabilities, not calculating the e ects of the air pressure of the entrapped air in the slug ow. Since slug ows is outside the scope of this work, this subject will not be further discussed here. 2.3 Numerical schemes to solve hyperbolic PDEs There are basically two types of numerical schemes for hyperbolic PDEs (such as Saint- Venant equations) solution, which are linear and the non-linear numerical schemes. The linear schemes, as the name indicates, use a linear combination of the previous time step values to calculate the uxes of conserved variables, while the non-linear schemes are based on the Riemann problem. Some linear schemes and its characteristics are [Toro, 2001]: Lax-Friedrichs (LxF): this rst order accurate (takes into consideration until the second derivative in a Taylor expansion) scheme is easy to implement, fast to calculate, and stable, however it is notorious for the high amount of numerical di usion in case of a discontinuity in the solution of the equations, such as in a supersonic air ow or a moving bore in a open channel ow. The algebraic expression for this scheme in nite volume is: 43 8> < >: ~Fn+1i = 1 2[ ~Fni 1 + ~Fni+1] + 1 2 t x( ~Uni+1 ~Uni 1) ~Un+1i = 1 2[ ~Uni+1 + ~Uni 1] + 1 2 t x( ~Fni 1 ~Fni+1) + tSni (2.27) in which i is the node number, n is the time step number, ~U is the vector of conserved variables and ~F is the uxes vector. Lax-Wendro (LxW): this scheme has second order accuracy (depends on the two last time steps values) and is more complex to calculate than the Lax-Fridriechs scheme because of the second order precision, being still considerably fast to calculate and easy to implement. A problem with LxF is the occurrence of very strong spurious oscillations on the vicinity of discontinuities in the solution (bores), which may easily make the code crash by making the water depth goes negative. In order to overcome this problem some researchers make use of ux limiters such as the TVD methods (Total Variation Diminishing) to mitigate those oscillations, as shown in [Toro, 2001]. The algebraic expression for McCormack variation of the LxW scheme (without TVD) is divided into two steps, being the rst the predictor and the second the corrector: (Predictor step) 8 >>< >>: ~Fn+12 i+12 = 1 2[~F n i+1 + ~F n i ] + 1 2 x t(~U n i ~U n i+1) + 1 2 tS n i ~Un+12 i+12 = 1 2[~U n i+1 + ~U n i ] + 1 2 t x(~F n i ~F n i+1) + 1 2 tS n i (Corrector step) 8 >>< >>: ~Fn+1i = ~Fni + t x(~Fn+12 i+12 ~Fn+12 i 12 ) ~Un+1i = ~Uni+1=2 + t x(~Fn+1i ~Fn+1i+1 ) + tSni (2.28) FORCE: This is a hybrid scheme which, in its original form, is calculate by averaging the uxes of Lax-Friedrichs and Lax-Wendro : FFORCEi+1 2 = (1 + )12(FLxF + FLxWi+1 2 ) (2.29) 44 with = 0:5. This scheme is less di use than Lax-Friedrichs and presents no oscilla- tions, however it demands a high computing e ort when compared to the other two alternatives due to the fact that it calculates both Lax-Friedrichs and Lax-Wendro for then calculate the actual uxes. For linear schemes as the ones presented, the degree of numerical di usion and spurious oscillations in the results, as well as the stability, are controlled by a non-dimansional pa- rameter named the Courant number, de ned by equation 2.30 [Sturm, 2001]. This condition related the propagation velocity of a ow feature predicted by a numerical scheme with the discretization-based velocity x= t. Explicit schemes require a Courant number below the unity for stability. Cr = juj+c x= t (2.30) where Cr is the Courant number. The open channel celerity c is substituted by the acoustic wave speed a a pipe ow and by the sound speed the the gas for a gas ow. According to [Godunov, 1959], the Godunov?s theorem states that when Courant num- ber is below unity, linear schemes results present either numerical di usion ( rst order schemes) or oscillations (higher order schemes) at the vicinity of bores. Those problems are worsened as the Courant number becomes smaller then unity, rendering the use of those linear schemes particularly problematic for problems containing several discontinuities, where part of the ow have a lower Courant number smaller than unity, considering that Cr 1 must be enforced in the whole solution domain. In the model proposed in this dissertation, the air and both free-surface and pressurized water ows must be calculated together for every time step, which creates the necessity of using the same time step for both the phases. This may lead to a low Courant number for the free surface water ow because the air speed of sound and the acoustic wave speed in pressurized ows can be two three orders of magnitude higher than the open channel 45 celerity c. As mentioned, Godunov?s theorem states that the air ow would su er from very high di usivity or very high spurious oscillations (sometimes even with a TVD method implemented) at the vicinity of bores/shocks if calculated with a linear scheme, however shocks are not anticipated in air ows in this work. In order to overcome the anticipated limitations of linear schemes in the simulation of hyperbolic partial di erential equations, non-linear schemes were developed. The rst one was by Godunov and solves the problem?s discontinuities not by performing a linear combination, buy instead by solving exactly the initial value problem proposed by Riemann. This solution is generally iterative and time-consuming, which lead to the development of several other non-linear schemes performing an approximate solution for the Riemann problem with the objective to derive alternative expressions for Finite Volume uxes across cell interfaces. Among such schemes one includes HLL [Harten et al., 1983], HLLC [Toro et al., 1994] and Roe?s scheme [Roe, 1981]. Both Roe and HLL rst order accurate schemes would be good options for implementa- tion in the proposed model. The HLL scheme is generally simpler than Roe scheme, however its uxes formulation is highly dependent on the ow cross-section geometry, which makes it cumbersome to be implemented for the case of a circular pipe. Since the cross-section geom- etry to be studied with in the present work is circular and Roe scheme does not depend on the cross-section geometry on the calculation of inter-cell uxes, the implementation of Roe scheme suggested by [Macchione and Morelli, 2003] was adopted. The algebraic expressions for the update of the conserved variables vector ~U (as in equation 2.21) is: ~Un+1i = ~Uni 2 (" (~Fni + ~Fni+1) X j j jj( w(j))i+1=2 r(j)i+1=2 # " (~Fni 1 + ~Fni ) X j j jj( w(j))i 1=2 r(j)i1=2 #) + t~Sni (2.31) where = t= x and F is the uxes vector. This equation was derived from the implemen- tation of Roe scheme in a 1-D Finite Volume framework. Also as shown in [Macchione and 46 Morelli, 2003], for the Jacobian matrix right eigenvectors r(j) the following inter-cell ow averaged variables, named Roe averages, should be considered for each interface (i;i+ 1) as: Ai+1=2 = pAiAi+1 (2.32) Qi+1=2 = pA iQi+1 + pA i+1Qip Ai +pAi+1 (2.33) For the celerity the Roe average is: c = s gI1i+1 I1iA i+1 Ai when Ai+1 6= Ai (2.34) c = s 1 2g(Ai+1 +Ai) 1 2g(bi+1 +bi) when Ai+1 = Ai or (I1i+1 I1i)(Ai+1 Ai) < 0 (2.35) so that the approximate eigenvalues matrix is characterized by the following eigenvalues and eigenvectors: eigenvalues: 1 = Q A + c 2 = Q A c (2.36) eigenvectors: r(1) = 12 c[1; 1]T r(2) = 12 c[1; 2]T (2.37) The variations w(1)(2) (strength of the wave crossing the (i,i+1) interface) at the point i+ 1=2 are expressed as follows: w(1)(2) = (Qi+1 Qi) + Qi+1=2 Ai+1=2 ci+1=2 (Ai+1 Ai) (2.38) 47 Chapter 3 Knowledge gap and objectives The problem of the ventilation of entrapped air in water mains is a complex and relevant problem, which remains to date poorly understood. Many experimental investigations, such as [Falvey, 1980], [Benjamin, 1968], [Pothof and Clemens, 2010] and others, focused on se- lected features of this problem, but the current knowledge is still very limited especially with regards to numerical modeling attempts, as shown in the literature review. Therefore, more numerical research and experimental investigations on water ows with entrapped air pock- ets still need to be done in order to provide engineers with a more complete understanding of this phenomena. The goal of this work is to obtain further insight on air-water interactions during water pipeline lling operations, with the overarching objective of developing a numerical model that may be used to simulate a priori lling operations in pipelines and detect operational issues related to the entrapment of air pockets. To achieve this objective, a numerical model is proposed to simulate the lling of pipelines and it applies the TPA approach presented by [Vasconcelos and Wright, 2009] to describe the water phase. Air phase modeling is performed either by using a discretized framework that applies the Euler equation or by using a type of UAPH model. A secondary objective was to assess the bene ts of using a discretized framework to simulate air phase. An experimental investigation was also conducted using a scale model apparatus that includes the essential features of a water pipeline. Key parameters in the problem are sys- tematically varied, including in ow rate, pipeline slope and ventilation degree. Experimental measurements included pressure, pressurization interface trajectories and in ow rates. Both 48 modeling alternatives for air phase were compared to experimental data and to the eld data of an actual water main lling event presented by [Vasconcelos et al., 2009b]. 49 Chapter 4 Methodology This chapter describes the methodology used to accomplish the stated objectives. In order to create a model able to simulate air pressurization e ects in water mains undergo- ing re lling operations, the model based on the TPA approach with air head described in [Vasconcelos and Wright, 2009] was coupled with two di erent approaches for the air phase calculation. This model aims to simulate water ows with air phase in pockets shrinking in volume due to the re lling, which means that they are not dragged and do not oat. Only the strictly necessary boundary conditions (for both phases) were developed, which consists on the ones required by the assessment of the model, experimental program, and eld conditions, in order to test the model. A systematic experimental investigation was conducted for this work due to the limited published data about this problem. A set of eld data about the lling process of a real water main in the city of Bras lia - Brazil used for model validation. This process of validation required a calibration of the energy dissipation based on the Manning?s n friction factor and on the contraction coe cient of the ventilation ori ce due to the lack of information on those two parameters and to the fact that the authors needed to approximate the behavior of an actual air valve by a single ori ce. 4.1 Numerical model Certain ow features of the water main pipeline lling problem were determinant in the model?s formulation so that it could describe the lling process adequately. With regards to the water phase, these features include: 50 Flow regime transition: addressed by using a ow regime transition model [Vasconcelos and Wright, 2009]; Post-shock oscillations at pipe- lling bore fronts: use of a numerical ltering scheme [Vasconcelos et al., 2009a]; Air pocket entrapment and pressurization: used either Euler equation or uniform air pressure head (UAPH) model; Free-surface and pipe- lling bores: Use of the approximate Riemann solver presented by Roe [Macchione and Morelli, 2003]; Dry water bed: assumed thin water layer (depth of 0.001 m) present in the whole dry portion of the pipe; and Solution stationarity: use approach presented by [Sanders et al., 2011]. Air phase in the model is represented by a well-de ned air pocket that is not signi cantly fractured. This pockets shrinks due to compression by the water phase that gradually occu- pies the lowest points in the pipeline pro le. Air is displaced and escapes through ventilation ori ces located at selected locations. According to [Tran, 2011] for such ow conditions air compression process may be considered isothermal. This assumption is used in both models used to simulate air phase during the lling process. The air phase is calculated as if the only outside connections (with atmosphere) occur at ventilation points, which are treated as ori ces for simplicity. Ideal ventilation with negligible air phase pressure head is assumed to exist prior to the formation of an entrapped air pocket, as it will be discussed later. When a pocket forms, it is delimited by the ventilation ori ce and a ow regime transition interface, either abrupt or gradual. In the proposed model, an air pocket is formed by the closure of a downstream valve or by the pressurization interface generated as water lls the lowest points of the pipeline, creating a pressurization interface. 51 Figure 4.2 presents a sketch of a typical application, whereas Figure 4.1 presents the overall structure of the proposed model. 4.1.1 Water phase modeling The TPA model, used in the water phase simulation, modi es the Saint-Venant equa- tions, enabling them to simulate both pressurized ows and free-surface ow regimes. This model has been improved in the past years and the alternative used here was presented in [Vasconcelos and Wright, 2009]. This alternative has a term that accounts for air phase pressure head, so that the modi ed St. Venant equations are, in divergence format: (U)t + Fx(U) = S(U) (4.1) where U = 2 64A Q 3 75; F(U) = 2 64 Q Q2 A +gA(hc +hs) +gApipehair 3 75; S(U) = 2 64 0 gA(S0 Sf) 3 75 (4.2) hair 8 >< >: = 0!free-surface ow without entrapped air pocket / Pressurized ow 6= 0!free-surface ow with entrapped air pocket (4.3) hs = 8> >< >>: 0 !free-surface ow a2 g A Apipe !pressurized ow (4.4) hc = 8> >>> >>< >>> >>> : D 3 3sin( ) sin3( ) 3 cos( ) 2 sin(2 ) !free-surface ow where = arccos[(y D=2)(D=2)] D 2 !pressurized ow (4.5) 52 Figure 4.1: Flow chart for the model calculation procedures 53 Ventilation ori ces Air pockets Water ow Energy grade line (EGL) Figure 4.2: Model physical scheme where U = [A;Q]T is the vector of the conserved variables, A is the ow cross sectional area, Q is the ow rate, F(U) is the vector with the ux of conserved variables, g is the acceleration of gravity, hc is the distance between the free surface and the centroid of the ow cross section (limited to D=2), hs is the surcharge head, hair is the extra head due to entrapped air pocket pressurization, is the angle formed by free surface ow width and the pipe centerline, D is the pipeline diameter, and a is the celerity the acoustic waves in the pressurized ow. The numerical scheme used in the implementation of the water phase model used the Finite Volume Method and the approximate Riemann solver of Roe, as presented in [Mac- chione and Morelli, 2003]. This choice was motivated by the signi cant discrepancy in the celerity values between the free-surface and pressurized ows, and between air and water ows. This discrepancy may be in the order of 2 or 3 orders of magnitude and yields an extremely low Courant number for the free-surface water ow, as shown in gure 4.3 and equation 4.6. Cr = juj+c x= t (4.6) 54 Pressurized owFree-surface ow a 200-1000 m/s c 300 m/s c 1 m/s Figure 4.3: Average Courant number for di erent portions of the ow when entrapped air pocket is present. In such conditions, free surface bores could be simulated under extremely low Courant numbers with no signi cant di usion oscillations. The mathematical formulation for the Roe scheme is presented in [Macchione and Morelli, 2003] and can be seen in equations 2.31 to 2.38. For dry bed regions of the ow, it was assumed that the ow depth would start as a minimum water depth of 1 mm. In such cases, the model would then predict the existence of a non-physical ow of this thin layer down the pipeline slope. To deal with this problem it was used the approach presented in [Sanders et al., 2011]. In this formulation, in order keep a minimum water layer with no motion and at the same time keep the stationarity of the solution, two criteria were followed in order to calculate the ow at a certain nite volume cell: one is based on the ratio between friction forces and the other based on the minimum submerged area of the cell. After computing these criteria to all cells in the domain, only the cells in which at least one of the two criteria is met have the ow calculated. The rst criteria is s>s where s was chosen to be 0.1 as suggested in [Sanders et al., 2011] and is based on the following equations: s = cDu 2 gRhjdz=dxj (4.7) 55 where Rh is the hydraulic radius and: cD = gn2Rh 1=3 (4.8) The second criteria follows the same logic, in which the ow is calculated if A > A , with A de ned in equation 4.9: A i = D 3 96j zij 9sin 2 6cos 2 +sin 3 2 (4.9) where = 2cos 1(1 2j zij=D) (4.10) and zi is the vertical elevation di erence between both cell sides. With this, the criteria to consider if a cell is a dry bed (weti = false) or is wet is: 8 >< >: weti = true; if Ai >A i and si >s i weti = false; otherwise (4.11) where weti is a boolean variable used to synthesize those two criterion in a single vari- able. The momentum equation 2.31 is solved for Qn+1i only if weti = true, with Qn+1i = 0 otherwise. 4.1.2 Water phase source terms Two source terms were considered for the water phase modeling, one accounting for pipe walls friction and another one accounting for pipe slope, both presented in [Sanders et al., 2011]. For the pipe walls friction source term, a semi-implicit formulation based on the Manning?s equation was used, while for the pipe slope a formulation which preserves stationarity of the solution was used. 56 The source term for gravity forces in a sloped pipe presented in [Sanders et al., 2011] follows the formulation presented in [Capart et al., 2003], which preserves stationarity of the solution, avoiding non-physical oscillations. In this formulation, two variables hs and wlevel represent the linear change of hs and wlevel in the cell. These are computed as: If vi = false then wlevel = 0 and: hsi = 8 >< >: 0 if si s zi2 [1 +cos( si=s )] if si < >: 0 if si s zi2 [1 +cos( si=s )] if si >>> >>> >>> < >>> >>> >>>> : r = dt=dx us = u n 1 +r(u n 2c n 1 u n 1c n 2 ) 1 +r( un1 +un2 +cn1 cn2 ) cs = c n 1 +rus(c n 1 c n 2 ) 1 +r(cn1 cn2 ) ys = wleveln1 +rjus csj (wleveln2 wleveln1 ) (4.21) and nally wleveln+11 = wlevelni + (un1 us)csg (4.22) where us, cs and ys are the interpolated velocity, celerity and water depth in between the rst and the second cells. If ow depth at inlet is less than the pipe diameter D then the surcharge head hs is set to zero. On the other hand, if wleveln+11 > D, ow at inlet is pressurized, hair is set to zero and hs is recalculated to match the piezometric head at the upstream end, calculated with the energy equation. With the depth and the pressure head updated, the ow area An+11 is updated and a new ow rate is then calculated with un+11 and An+11 . 59 The downstream boundary condition can be a fully opened or closed valve. For the case in which the valve is fully opened, the approach called transmissive condition presented in [Toro, 2001] is used: 8 >>> >>< >>> >>: wlevelnNo = wlevelnNo 1 hsnNo = hsnNo 1 unNo = unNo 1 (4.23) where the subindex No means the number of nodes used in the discretization. For the case in which the downstream boundary condition is a closed valve the boundary condition is calculated enforcing the relevant characteristic equation and zero velocity at the downstream end. The MOC Hartree is then used as: 8 >>>> >>> >>>> < >>> >>> >>> >>: ur = u n No +r( u n Noc n No 1 +c n Nou n No 1) 1 +r(unNo unNo 1 +cnNo cnNo 1) cr = c n No +rur(c n No 1 c n No) 1 +r(cnNo cnNo 1) Sfr = n2 urjurj(Rhn No)1:333 Kr = ur +gyrc r g(Sfr S0No) t; (4.24) and nally, with the velocity unNo set to zero: wleveln+1No = crg (Kr) (4.25) where ur, cr, Sfr and Kr are the interpolated velocity, celerity, friction slope, and a constant accounting for the interpolated previous time step values in between the last and the second last cells. If wleveln+11 > D, the ow depth at the downstream end becomes pressurized. In such case, wleveln+11 is set to the value of the pipe diameter and hs is set as the extra head of the cross section minus D, and hair is set to zero. 60 4.1.4 Air phase modeling In the proposed model, air is initially considered as a continuous layer over the water layer (strati ed water in free-surface ow mode). During the simulation, air is handled in one of two following manners: 1. At the rst stages of the lling, when the air within a given pipe reach consists in an entire layer that is connected to atmosphere at the ventilation ori ce and at its lowest point within the pipe (ideal ventilation), it is assumed that the air pressure in the entire layer is approximately atmospheric, and air velocity is assumed negligible. This condition persists until a pocket is formed at the lowest point during the lling process. 2. At the second stage, once an air pocket is formed, the connection to the atmosphere at its lowest point is lost, and air phase pressure is expected to be varying above atmospheric values, requiring calculations with either one of the two presented air phase models to determine its pressure and in uence in the water ow. An algorithm was developed to track air pocket volume, start and end nodes as it shrinks in order to simulate its behavior with any of the two models. For this, the mechanism considered for air pocket formation is the isolation of an air mass due to the development of a ow regime transition interface or a closed downstream valve. As mentioned, it is assumed that during a pipeline lling event this air pocket will be delimited by a ventilation ori ce and a ow regime transition interface (or alternatively a closed valve). This interface will move mainly towards the air pocket ventilation point, compressing the air pocket and forcing its elimination through the ventilation ori ce, as it is sketched in Figure 4.2. The rst step for tracking the air pocket is to check if the air mass has contact from down or upstream with sources of perfect ventilation. For this, a routine scans the pipe twice, leaving from the rst cell (upstream) towards the last cell (downstream) and from last towards the rst, to check the availability of ideal ventilation from the left and/or from 61 Downstream to upstream scan Upstream to downstream scan False False False False False False False False False False False False False True True True Figure 4.4: Pocket nd graphical example the right sides of the air masses. Figure 4.4 shows both scans and their return true or false for di erent cells in the pipe. The scan from upstream to downstream will set all the nodes of the pipe since its beginning as false for contact with ideal ventilation sources from the upstream until it reaches a source of perfect ventilation, which means until the downstream end of the pipe. On the other hand, the scan routine which starts from the last node (downstream) will set all the pipe nodes as true for contact with perfect ventilation sources from the downstream sides until it nds a cell which has water touching the crown of the pipe, so that from this point the contact with the downstream source of perfect ventilation is blocked. After the pipe is scanned from both sides, the cells with contact to ideal ventilation sources from both sides set to false and water level wdepth below the crown of the pipe will be set as true for air pocket presence. Those cells will require the calculation of the air phase for them, having a non-atmospheric pressure associated with them in each time step. Another important part of the algorithm is to determine the ventilation status of the cell. A ventilated cell means that it has direct contact with air, being this air mass is a air pocket or not. Therefore, the variable for ventilation vi is set to true as long as wdepthi < >: nNoair = nNoair 1 uairnNoair = uairnNoair 1 + 2bvel (4.37) where bvel is the pressurizing front velocity and Noair is the number of cells in the air pocket. This boundary condition accounts for the pressurizing front velocity, so that when it moves from a cell to another the loss of air due to the shrinking of the pocket was previously com- pensated. However, since the source terms which account for the cross-section area shrinking already accounts for the movement of this mass, both cannot be used simultaneously or else the air present in the last cell would be counted twice. If the area variation source terms are used, bvel equals zero. On the other hand, if the wall movement is considered in the boundary condition, the mentioned source terms are not used and bvel needs to be calculated. For this case, two formulations were considered: bvel = (uA)i 1 k (uA)i+kA i 1 k Ai+k (4.38a) bvel = xi m imlast t (4.38b) where i is the cell number of the pressurization front, k is the number of cells before and after the pressurizing front to be considered in the control volume surrounding the pressurizing front, m is the last time step when the pocket?s last cell moved from one cell to another, mlast 66 is the time step when the pressurizing front moved considered for the last bvel calculation, and t is time interval between both instants. For instance, one can choose to recalculate the pressurizing front velocity when the front moves two cells, to that m would be the current time step and mlast would equal m 2. It is important to notice that the rst formulation in equation 4.38 can be used for the whole simulation time while the second one needs a velocity to start with since no previous values exist then. For that, the rst formulation is used until the pressurizing front moves a certain distance and from then on the second formulation can be used. At the uppermost point, the ventilation ori ce boundary condition for the Euler equation approach is similar to the UAPH model in the sense that both apply a continuity equation along with the ori ce equation. The continuity equation for this boundary condition is: tAn+1air n+11 un+11 = An+1air x( n+11 n1 ) +Mairn+1out (4.39) where Mairout is the air mass discharged through the ventilation ori ce, calculated using equation 4.35. Equation 4.39 is solved for un+11 using the Riemann invariants for the isother- mal, one-dimensional, primitive version of the Euler equation [Pulliam, 1994]: uair1 = uair2 (log10 2 log10 1) (4.40) 4.4 Experimental program An experimental investigation was conducted to gather insights on the characteristics of the pipeline lling problem, and also to validate the proposed numerical model. The experimental apparatus was set up and works started in May, lasting until late June, 2011. The current apparatus was inspired in the one presented in [Trajkovic et al., 1999]. 67 Ventilation ori ce Downstream valve Water inlet from a reservoir 6.69 m 10.96 m Pressure transducers Recirculation pumpADV Figure 4.5: Test representation 4.4.1 Experimental apparatus setup A sketch of the experimental apparatus is presented in Figure 4.5. The experimen- tal apparatus was a 10.96-m long, 101.6-mm diameter clear PVC pipeline with adjustable slope. At the upstream end, a 0.66 m3 capacity water reservoir supplied ow to the pipeline through a 50 mm ball valve; at the downstream end ow discharged freely through a 101.6 mm knife gate valve into a 0.62 m3 reservoir, and ow was subsequently recirculated with pumps. Right after the inlet control valve, a T junction was installed in the pipe so that di erent caps with ventilation ori ces could be installed. Initial, steady ow conditions were such that free surface ows exist at the whole pipeline, as the downstream gate was fully opened. A sudden closure maneuver (within 0.3 seconds) of the knife gate valve at the down- stream end of the pipeline blocked the downstream ventilation, triggered a backward-moving pressurizing interface, and resulted in the entrapment of an air pocket. These air pockets became pressurized as water accumulated at downstream end of the pipe pushed the air mass through the ventilation ori ce in the beginning of the pipe. Two pressure transducers MEGGIT-ENDEVCO 8510B-5 were installed at selected locations along the pipe (upstream end and 39% of the pipeline length from the knife gate valve). Transducer results were calibrated each experimental run with the aid of four digital manometers, with of 3.5-m H20 maximum pressure head and 0.3% accuracy. Flow rates were measured with a MicroADV 68 positioned in the recirculation system, and con rmed by a paddle-wheel ow meter. The presented apparatus was inspired in the one presented in work by [Trajkovic et al., 1999]; among the main modi cations was the absence of intermediate ventilation points. 4.4.2 Experimental procedure 1. With the desired slope set in the pipeline, the pumps were started; valves near the pump were opened enough to provide selected steady ow rate to the system; 2. The desired ventilation ori ce was installed; 3. When water level at upstream reservoir attained steady level, it was performed readings at all manometers, as well the upstream reservoir head; 4. The data logging was started for the pressure transducers, for the ow rate measured with the MicroADV, and for the head at the upstream reservoir measured with the manometer; 5. The downstream knife gate valve was rapidly maneuvered and closed entrapping an air pocket and creating a backward-moving pressurization front; 6. Digital cameras (30 FPS) recorded the whole pipe lling process, one of them tracking the bore and another one tracking the pressurization interface; 7. When the pressurization interface approached the ventilation ori ce, the ventilation ori ce was shut to avoid water leakage; 8. Pump was shut down and control valves closed so that pressure could attain a static level; 9. Manometers were read and data collection with logging software stopped; 69 The use of two cameras to track the in ow/pressurization front was particularly neces- sary in the case when interface breakdown [Vasconcelos and Wright, 2005] was noticed. Oth- erwise, just one camera tracking the pipe- lling bore front was used. The described exper- imental program varied systematically ow rates, ventilation ori ce diameters and pipeline slopes. Table 4.1 presents the ranges of the tested experiment variables, with a total of 36 conditions tested. A minimum of two repetitions were performed to ensure consistency of the data collected. Variables Tested values Normalized values Flow rate 2.53, 3.79 and 5.05 l/s 0.245, 0.368 and 0.490 Slope 0.5, 1 and 2% N/A Ventilation ori ce diameter 0.63, 0.95, 1.27, and 5.06 cm 0.0625, 0.09375, 0.125, 0.5 Table 4.1: Experimental variables. Flow rate was normalized as Q = Q=pgD5 and venti- lation diameter as d = dorif=D where D is the pipe diameter 70 Chapter 5 Results and analysis In this chapter, a set of numerical and experimental results were presented. It starts with a discussion focused on the experimental results. This discussion will focus on quantitative results, such as pressure and pressurization front trajectories, and also in the ow features, such as interface breakdown and depression waves. Later, comparisons between the modeling alternative for air phase and between models and experimental/ eld data are presented and discussed. The tests were created with three main purposes: 1. Evaluate the feasibility of running parallel discretized models for both air and water phases; 2. compare the di erent model approaches to simulate air phase during pipeline lling events; and 3. validate the model with experimental and eld data. The development of the model in its current version led to the creation of intermediate functional versions were tested, a process by which the strengths and weaknesses of di erent proposed ideas are evaluated. In the following two sections, intermediate and nal versions of the model are presented with a evaluation of each version. Only the nal version of the numerical model was compared with experimental and elds results. 5.1 Experimental results Figure 5.1 shows the pressure history close to the ventilation ori ce for all tested cases in the experimental program with normalized ori ce diameter d orif = dorif=D smaller or 71 equal to 0.125. The transducer at that station was located at the pipe crown, so it measured air phase pressures for most of the lling processes. As anticipated, higher pressurization levels were observed for smaller ventilation ori ces and that the lling time was smaller for higher ow rates. Air phase pressure head results were not signi cantly di erent for varying pipeline slopes. On the other hand, there was a slight di erence in the lling time between di erent pipeline slopes for a given ventilation ori ce and ow rate. This di erence is attributed to the di erent initial water levels in the apparatus prior to the closing of the knife gate valve at the downstream end. Also in Figure 5.1, it can be noticed that for the smallest ventilation the air phase pressure head kept increasing during the lling process, indicating steady ow for these cases was not attained. Figure 5.2 presents the pressure head histories for a point x = x=L = 0:39 (measured from the knife gate valve) for experimental runs with d orif 0:125 with d orif = dorifD . These pressures were also measured at the pipeline crown. A sudden step up in the pressure values the moment in which the ow regime transition interface reached the transducer. As in the case of pressure measurements at the ventilation ori ce, these pressures kept increasing due to the increase in the air pressure for the smallest ventilation case. The magnitude of the jump in the pressure readings was an indication of the strength of the pipe- lling bore front, and it increased for larger in ow rates and ventilation ori ces. The absence of this discontinuity was a sign of either gradual pressurization interface and/or the occurrence of interface breakdown feature due to the interaction of the in ow front and the depression wave generated at the inlet caused by air pressurization. The relevance of this ow feature is that its occurrence may pose di culties to pipeline lling models that use well-de ned in ow interfaces as a modeling hypothesis, such as the models by [Liou and Hunt, 1996] and [Izquierdo et al., 1999]. To further illustrate the impact of the interface breakdown feature, Figure 5.3 presents two set of trajectories of moving pressurization interfaces for normalized ow rates of Q = 72 Q=pgD5 = 0:245 and 0:490 and 2% slope, measured for all tested ventilation diameters. All these interfaces start as pipe- lling bore fronts at x = 0 as the gate valve is closed, and such bores lasted until x 0:28. For both ow rates interface breakdown was noticed when the smallest ventilation ori ce was used. Upon occurrence of interface breakdown, the pipe- lling bore becomes an open-channel bore that moves more slowly toward the ventilation ori ce, leaving an air intrusion on its top. Interestingly, as the backward-moving bores approached the ventilation ori ce there seems to be an acceleration on their motion. The cause for this change in bore velocity is not determined at this point. When interface breakdown occurred, interface measurements included both the position of the open-channel bore and the pressurization front. The latter was a gradual transition, and immediately following the interface breakdown it was noticed that the pressurization front retreated. Soon afterwards, the pressurization front resumed the motion toward the ventilation ori ce, trailing the open-channel bore. Figure 5.4 presents the trajectories for the condition Q = 0:245 and 1% slope, for all four tested d orif. The largest one (d orif = 0:50) has not generated any sign of air pressurization, and the pipe- lling bore kept its shape as it propagated toward the ventilation point. However, for smaller ventilation ori ces there was the occurrence of interface breakdown, and the the trajectory of the pressurization fronts (thin lines) are plotted along the trajectories of the pressurization bores. For the cases when d orif = 0:125 and d orif = 0:0938 the trajectory of the pressurization front was approximately parallel to the backward moving bore, which was moderately slowed by the interface breakdown. For the smallest ventilation (d orif = 0:0625), on the other hand, the velocity of the pressurization bore was signi cantly reduced, with a much larger separation between the pressurization front and interface breakdown. 73 0 1 2 3 4 5 6 7 Q = 0:245;d orif = 0:125 0 1 2 3 4 5 6 7 h air = hair =D Q = 0:245;d orif = 0:09375 0 1 2 3 4 5 6 7 0 1 2 3 Q = 0:245;d orif = 0:0625 Q = 0:368;d orif = 0:125 Q = 0:368;d orif = 0:09375 0 1 2 3 t = t=(L=(gD)0:5) Q = 0:368;d orif = 0:0625 Q = 0:490;d orif = 0:125 Q = 0:490;d orif = 0:09375 0 1 2 3 Q = 0:490;d orif = 0:0625 S = 0:5% S = 1:0% S = 2:0% Figure 5.1: Air phase heads close to the ventilation ori ce for all tested conditions. 74 0 1 2 3 4 5 6 7 0 1 2 3 Q = 0:245;d orif = 0:0625 0 1 2 3 4 5 6 7 h air = hair =D Q = 0:245;d orif = 0:09375 0 1 2 3 4 5 6 7 Q = 0:245;d orif = 0:125 0 1 2 3 t = t=(L=pgD) Q = 0:368;d orif = 0:0625 Q = 0:368;d orif = 0:09375 Q = 0:368;d orif = 0:125 0 1 2 3 Q = 0:490;d orif = 0:0625 Q = 0:490;d orif = 0:09375 Q = 0:490;d orif = 0:0125 S = 0:5% S = 1:0% S = 2:0% Figure 5.2: Air phase heads close to the ventilation ori ce for all tested conditions. 75 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 0.2 0.4 0.6 0.8 1 t = t=( L= pg D) x = x=L a) Interface breakdown 0 0.2 0.4 0.6 0.8 1 x = x=L b) Interface breakdown d orif = 0:0625 d orif = 0:0938 d orif = 0:125 d orif = 0:5 Figure 5.3: Trajectories of moving bore for recirculation ow rates of 2.53 l/s (a) and 5.05 l/s (b) and slope of 2% 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 0.2 0.4 0.6 0.8 1 t = t=( L= pg D) x = x=L Interface breakdown d orif = 0:0625d orif = 0:0938 d orif = 0:125d orif = 0:5 Figure 5.4: Trajectories of moving bore (thick line) and pressurization interface (thin line) for recirculation ow rate of 2.53 l/s and slope of 1% 76 5.2 First version of TPAir 5.2.1 Assessment of the Euler equation approach Ventilation ori ce Downstream valve Water inlet 10.83 m Device providing constant head and ow rate Figure 5.5: Model?s rst version layout. The rst test presented in [Trindade and Vasconcelos, 2011] was thought as a means to assess the idea of running parallel models for both phases, which means the coupling of the Euler equations for air phase with the Saint Venant equations for the water phase. The test case for this rst model was comprised by a 10 cm diameter, 8.3 m-long pipeline slope was 2.7% downward with the ow, and free discharge conditions existed at the downstream end. At the upstream end ow rate and depth were kept constant due to the supercritical ow condition at that location. The general layout is presented in 5.5. In this rst stage, the rst functional version of the model with Euler equations approach was developed, using: Water source terms presented in [Sanders et al., 2011]; Single phase Euler equation test model for air phase in order to develop the boundary conditions; Upstream ventilation ori ce air boundary condition; Downstream moving bore (bvel) calculated as in section 4.3.2 using air mass balance (equation 4.38a) for the whole calculation; and Air phase described by the isothermal, one-dimensional Euler equations coupled with the Saint Venant equations to model water phase. 77 The initial set up for the model consists in setting a ow rate and a initial depth predicted with Manning equation for the whole pipe. In order to get a better prediction for steady state, the model runs for 60 seconds with a transmissive boundary condition (equation 4.23). Unsteady ow conditions arise from abruptly closing the downstream valve at the time 60 seconds, after steady state is reached, blocking the ow and instantly generating an upstream- moving pipe lling bore, which expelled the air phase through the ori ce at the upstream end. The main di erences from the proposed example to the original con guration by [Trajkovic et al., 1999] are the ow rates tested (larger ows) and the fact that ventilation in the system is limited. Ideal ventilation condition is assumed to occur only at the downstream end prior to the valve closure. The parameters varied for the test cases were the diameter of the ventilation ori ce and the initial steady ow rate in the system, which resulted in varying initial water level. The choice of air ori ce diameters was based on the work of [Zhou et al., 2002a] plus a larger ori ce which provides near-ideal ventilation at the upstream and, yielding a total of three di erent ventilations. The smallest ventilation ori ce corresponded to a case in which the air cushioning e ects were dominant while the second larger ori ce corresponds to the case in which air pressurization was minimal. Notice that albeit the imposition of Q in the pipeline is not a realistic type of boundary condition, its still useful for the model assessment that is the main goal of this research phase. Variables Tested values Q Normalized ow rates ( QpgD5 ) 0.67 and 1.32 d orif Normalized ventilation ori ce diameter (dorifD ) 0.028, 0.114 and 0.171 Table 5.1: Tested values for Euler equation approach assessment (ori ce sizes based on the work by [Zhou et al., 2002a]). 78 Figure 5.6 shows the variation of the air pressure heads in the discharge ori ce over time for both tested initial depths y0 = 0:5D and y0 = 0:8D ( gures 5.6a and 5.6b, respec- tively). This gure indicates that the pipe lling time was smaller for y0 = 0:8D, which was anticipated considering the faster motion of the upstream propagating bore upon valve closure. Since the in ow rate (and therefore uA in equation 4.38) is higher for the case when y0 = 0:8D, while the area di erence A2 A1 is smaller, the bore moves much faster it larger depth. The resulting conditions for the air ow is that at the pressurizing bore the air phase pushed with a higher velocity, resulting in higher pressurization heads calculated at the air discharge cell with equation 4.35. Figure 5.6: Pressures comparison at the ventilation ori ce for a) initial depth of 0:5D and b) 0:8D For the two larger ori ce sizes in gure 5.6a it can be noticed that the air ow reaches the steady state at about 60% of the simulation time. At this point, the air pressure gradient is uniform throughout the air phase. One notices the large di erence in the nal pressure values, as the smallest ori ce yields a pressure head over two orders of magnitude higher than the one with the largest ori ce. The discrepancy in the lling time, which should be the same because the in ow rate and pipe volume are xed, comes from the relatively small celerity used in the computations, around 85 m/s. This situation is analogous to running simulations with a wide Preissmann slot value, as presented by [Vasconcelos and Wright, 2004] and discussed in section 2.1.2. While this version of the model is still unstable for higher celerity values (e.g. 200-300 m/s), this limitation was addressed in later versions of 79 the proposed model. Interestingly, the maximum air pressure predicted for higher initial ow rate conditions increased with the initial ow depth, except for the smallest ori ce diameter. For the smallest ventilation ori ce case in gure 5.6a and for all cases in gure 5.6b, small pressure oscillations can be noticed. Those appear due to the air pressure waves re ecting at the edges of an smaller air cavity, which occurs in a much smaller period than the total lling time. Such an e ect could not be represented in a model which considers uniform pressure along throughout air pocket. While the magnitude of these oscillations is small at this particular application, it is thought that for longer air pockets this may become a more signi cant issue, which could justify the use of Euler equations model in order to obtain accurate results. Figure 5.7 presents two pairs of pressure histories illustrating the piezometric pressure evolution with time at a point close to the downstream valve (x = 7:7m) and another point further upstream (x = 5:7m), as in [Trajkovic et al., 1999], for the di erent tested depths. Figures 5.7a and 5.7b indicate slight di erences in the pressure histories, with an abrupt pressure rise in Figure 5.7a at t = 1:5s caused by the bore arrival at that point. A similar trend was also noticed for the largest initial water depth, with the bore sweeping at x = 5:7m at the instant t = 0:23s, albeit with higher pressures and much more rapid lling. Figure 5.8 presents the velocity history at the upstream station for both simulated ows. Generally, the proposed model predicts an increase of the velocity magnitude at the discharge point, with a value that approaches the celerity of the upstream propagating bore. As anticipated, this increase is more pronounced for larger discharge ori ces. For the largest water depth, there is a nal drop in the ow velocity, which is possibly due to the decrease in the upstream ow depth caused by air pressurization. Scatter in the data that appears for the smallest ori ce results is due to the model di culty in determining precisely the exact location of the bore interface at each time step, particularly for high air pressures conditions. Figure 5.9 presents the piezometric pro le at selected instants for both initial ow rates tested. One notices that the total piezometric pro les are much higher than the piezometric 80 Figure 5.7: Pressure history at a) x = 7:7m and y0 = 0:5D, b) x = 5:7m and y0 = 0:5D, c) x = 7:7m and y0 = 0:8D, and d) x = 5:7m and y0 = 0:8D pro le solely due to water phase because of the strong air pressurization, and as expected this pressurization increases for smaller discharge ori ces. At t = 3s, for the case when y0 = 0:5D, one notices a fairly uniform pressure gradient for the air phase (at left) followed by a small pressure increase, which corresponds to the arrival of the bore at the location. On the other hand, for the case when y0 = 0:8D, the air pressure is not uniform ahead of the bores, particularly for the two smaller ori ce results. This non-uniformity will be re-evaluated in future versions of this model, after the inclusion of friction terms between air and water phases and between air and the pipe walls. 5.2.2 Comparison between both approaches [Vasconcelos and Trindade, 2011] introduced the UAPH approach to the model frame- work and compared both approaches to simulate the air phase in the context of large-scale pipelines by means of a numerical investigation. To perform the comparison between both approaches, a two-reach pipeline con guration with an intermediate low point was proposed. 81 Figure 5.8: Air velocity at the upstream end of the pipe for tested ow rates Figure 5.9: Head pro les along the pipe at selected instants Both con gurations parameters are presented in table 5.2. Friction factor for air phase cal- culation was based on Manning roughness values. Figure 5.10 presents a sketch of the hypothetical pipeline. L1 (m) L2 (m) S1 (%) S2 (%) D (m) n Conf. 1 1000 600 4 2 1.0 0.010 Conf. 2 4000 2000 Table 5.2: Pipeline con gurations parameters The rst reach was simulated using the TPA approach for the water phase and Euler equations or uniform pressure approach for the air phase. The second reach was treated as the downstream boundary condition of the rst reach, being calculated using method of characteristics coupled with a variable depth until the whole connection is lled; afterwards the second pipe is computed using rigid column theory. 82 Figure 5.10: Hypothetical real sized pipeline sketch At the upstream end a xed head reservoir, which was calculated using equations 4.20 with xed reservoir head and water depth and 4.21, provides water in ow (three di erent pairs of ow rates and heads considered) but does not allow for air escape; the air valve im- mediately downstream from the reservoir valve provides limited ventilation for the upstream pipe reach. Three ventilation sizes for the air valve were considered in terms of the diameter of the pipe: 5%, 10% and 15% of the pipe diameter. Overall nine conditions were tested for each one of the modeling approaches, and the variables tested with respective tested range are presented in table 5.3. Variables Tested values Q Initial normalized ow rates ( QpgD5 ) 0.125, 0.250, 0.500 d orif Normalized ventilation ori ce diameters (dorifD ) 0.05, 0.1, 0.15 Table 5.3: Long pipe tested values. The qualitative behavior of the ow was generally the same for all tested cases. Some time after the ow was initialized (with a chosen constant reservoir head able to supply the speci ed initial ow rate), water reached the bottom of the system (connection between pipes) and started to ll this point as a pool. At this point, the system is still considered su ciently ventilated so that the air pressure inside the upstream pipe equals the atmospheric pressure. After a certain time the whole cross-section at this bottom junction was lled and a water column formed at the downstream upward sloped pipe. When this happened a 83 Figure 5.11: Comparison of heads at the ventilation ori ce for Q = 0:125 pressurizing bore started to move upstream in the upstream pipe starting to pressurize the air, which was expelled through the ventilation ori ce. The pressurized water ow in both pipes then oscillated in a way that resembles a U-pipe oscillation until the air got expelled. Figures 5.11, 5.12 and 5.13 show head and water ow rates for the upstream end of the pipe, where the ventilation ori ce and water inlet are located. It can be seen that the ow rate decreased as the pressurizing bore moves upstream while the air pressure increased. While air pressure continued to rise for the smallest ventilation tested, sometimes to signi cant levels, it becomes stable for ori ces diameters of 10% and 15% of the water main diameter. Because ow rate decreases proportionally to the air pressure, the lling time may increase signi cantly (up to 25%) when compared with a condition with larger ventilation size, as it would be expected. Figures 5.11 to 5.13 also show that both models pressure predictions at the ventilation ori ce were highly similar when air pressurization was low, which means below or around 1 m of water column. 84 Figure 5.12: Comparison of heads at the ventilation ori ce for Q = 0:250 Figure 5.14 shows the air velocity at the upstream end of the pipe (inside the pocket, not through the ventilation ori ce) for all three ow rates tested and d orif = 0:15. In those gures it can be noticed that for higher ow rates, and therefore higher pressures as seen in gures 5.11 to 5.13, the air velocity was higher, as it would be expected. Values in gure 5.14 are negative because of the adopted referential, which is positive in the downstream direction and negative to upstream. The magnitudes of the observed velocity oscillations seem exaggerated and are possibly caused by a limitation on how this phase is represented in the numerical code. The pressure history at x = 500m for d orif = 0:05 is shown in gure 5.15. In this gure it can be seen that while the connection was not lled the pressure at that point is zero, raising to a small value after air pressurization starts (connection is lled). After that the pressure rose abruptly when the returning bore reaches that point, rising steadily later until the pipeline is lled, being this rising more pronounced for larger rates. 85 Figure 5.13: Comparison of heads at the ventilation ori ce for Q = 0:500 Figure 5.14: Air velocities predicted using Euler equations model for all three Q and d orif = 0:15 86 Figure 5.15: Pressure history at x = 500m for d orif = 0:05 87 5.3 Second version of TPAir In the second and nal version of the model, the calculations of the air phase downstream boundary condition in the Euler equations model are performed by assuming bvel = 0 (section 4.3.2, equations 4.37 and 4.38) and using the source terms described in equation 4.32. With this modi cation the continuity errors decreased considerably and the code became more stable. 5.3.1 Comparison between experimental results and numerical model predic- tions The comparison between experimental results presented in section 5.1 and corresponding numerical predictions is presented in Figures 5.16, 5.19 and 5.20. All comparison were performed assuming that the PVC pipeline Manning roughness coe cient was n=0.0085, and the wave celerity assumed for the PVC pipe was 200 m/s. Pressure head predictions by the model for all the tested cases that resulted in air pres- surization for 0.5% slope are compared with experimental results in Figure 5.16. Both Euler and UAPH approaches models showed good agreement with experimental data for most of the cases, specially for the lower ow rates (Q = 0:245 and Q = 0:368). In some of the numerical predictions there were non-physical, high-frequency oscillations on the air pres- surization results as the pocket volume shrank to zero. Results obtained in the intermediate point (x = 0:39 measured from the knife gate valve) for the same slope are presented in Figure 5.19, and indicate fair agreement between numerical and experimental results. There is a tendency of the numerical model to anticipate the arrival of the pressurization front at the station, which results in the early jump in the pressure results. However, in general the predicted pressure increase over time matched what was measured by the transducers. Another issue with the numerical predictions was linked to interface breakdown occur- rences. The proposed numerical model (both air phase model implementations) was able to predict the onset of the interface breakdown as the interaction of the depression wave 88 0 1 2 3 4 5 6 7 0 1 2 3 Q = 0:245;d orif = 0:0625 0 1 2 3 4 5 6 7 h r es = hr es =D Q = 0:245;d orif = 0:09375 0 1 2 3 4 5 6 7 Q = 0:245;d orif = 0:125 0 1 2 3 t = t=(L=pgD) Q = 0:368;d orif = 0:0625 Q = 0:368;d orif = 0:09375 Q = 0:368;d orif = 0:125 0 1 2 3 Q = 0:490;d orif = 0:0625 Q = 0:490;d orif = 0:09375 Q = 0:490;d orif = 0:125 Euler equation UAPH Exp rep 1 Exp rep 2 Figure 5.16: Experimental and predicted pressures at the ventilation ori ce for all tested cases with slope of 0.5% 89 0 1 2 3 4 5 6 7 0 1 2 3 Q = 0:245;d orif = 0:0625 0 1 2 3 4 5 6 7 h r es = hr es =D Q = 0:245;d orif = 0:09375 0 1 2 3 4 5 6 7 Q = 0:245;d orif = 0:125 0 1 2 3 t = t=(L=pgD) Q = 0:368;d orif = 0:0625 Q = 0:368;d orif = 0:09375 Q = 0:368;d orif = 0:125 0 1 2 3 Q = 0:490;d orif = 0:0625 Q = 0:490;d orif = 0:09375 Q = 0:490;d orif = 0:125 Euler equation UAPH Exp rep 1 Exp rep 2 Figure 5.17: Experimental and predicted pressures at the ventilation ori ce for all tested cases with slope of 1.0% 90 0 1 2 3 4 5 6 7 0 1 2 3 Q = 0:245;d orif = 0:0625 0 1 2 3 4 5 6 7 h r es = hr es =D Q = 0:245;d orif = 0:0938 0 1 2 3 4 5 6 7 Q = 0:245;d orif = 0:125 0 1 2 3 t = t=(L=pgD) Q = 0:368;d orif = 0:0625 Q = 0:368;d orif = 0:0938 Q = 0:368;d orif = 0:125 0 1 2 3 Q = 0:490;d orif = 0:0625 Q = 0:490;d orif = 0:0938 Q = 0:490;d orif = 0:125 Euler equation UAPH Exp rep 1 Exp rep 2 Figure 5.18: Experimental and predicted pressures at the ventilation ori ce for all tested cases with slope of 2.0% 91 0 1 2 3 4 5 6 7 0 1 2 3 Q = 0:245;d orif = 0:0625 0 1 2 3 4 5 6 7 h r es = hr es =D Q = 0:245;d orif = 0:09375 0 1 2 3 4 5 6 7 Q = 0:245;d orif = 0:125 0 1 2 3 t = t=(L=pgD) Q = 0:368;d orif = 0:0625 Q = 0:368;d orif = 0:09375 Q = 0:368;d orif = 0:125 0 1 2 3 Q = 0:490;d orif = 0:0625 Q = 0:490;d orif = 0:09375 Q = 0:490;d orif = 0:125 Euler equation UAPH Exp rep 1 Figure 5.19: Experimental and predicted pressures at x = 0:39 (from downstream valve) for all tested cases with slope of 0.5% 92 and pipe lling bore resulted in an open channel bore. Unlike the experiments, the pre- dicted pressurization front does not retreat following the breakdown. There was a slight over-prediction of the pressure head in cases when interface breakdown occurred; results obtained with Euler equation indicate the instant of the breakdown by a second, smaller increase in pressure head at t 0:2. This discrepancy, however, was not signi cant and has not compromised the general accuracy of the numerical model. Figure 5.20 contains a comparison between experimental results and numerical predic- tions of the pressure head variation at the upstream reservoir during the lling events for cases with 0.5% slope. One notices generally good agreement between experimental and numerical results. The reservoir discharge point within the pipeline was in free surface ow regime until the arrival of the backward moving pressurization interface. One recalls that prior to the knife gate valve closure, the reservoir head was steady. Considering that the in ow rate into the reservoir was constant, the increase in reservoir head following the knife gate valve maneuver indicates a drop in the in ow rate admitted into the pipeline due to the almost instantaneous air pressurization. Steeper reservoir pressure head increase is linked to stronger air pressurization, and the drop in ow rate resulted in the depression wave which trigged interface breakdown events. Error bars plotted in Figure 5.20 re ect the accuracy of the manometer used in the experiments ( 0:01m). It can be observed that both models yielded similarly accurate results when compared to experimental data despite the two considerably di erent approaches to simulate air phase. An aspect to be considered is the computational e ort involved in each air phase modeling alternative. In general, the simulation time using the Euler equation model was over 9 times larger than one required by the UAPH model approach for the comparison with the experimental results. Not only due to the additional model complexity, but the enforcement of the Courant condition for the air phase simulation using resulted in even smaller time steps as the celerity of the air phase was in the order of 300 m/s. 93 0 1 2 3 4 5 6 7 0 1 2 3 Q = 0:245;d orif = 0:0625 0 1 2 3 4 5 6 7 h r es = hr es =D Q = 0:245;d orif = 0:09375 0 1 2 3 4 5 6 7 Q = 0:245;d orif = 0:125 0 1 2 3 t = t=(L=pgD) Q = 0:368;d orif = 0:0625 Q = 0:368;d orif = 0:09375 Q = 0:368;d orif = 0:125 0 1 2 3 Q = 0:490;d orif = 0:0625 Q = 0:490;d orif = 0:09375 Q = 0:490;d orif = 0:125 Experimental Euler equations UAPH Figure 5.20: Experimental and predicted reservoir heads for all tested cases with slope of 0.5% 94 Another aspect to be considered is the continuity error for the air phase. The Euler equation model had an average continuity error of 7.97%, while the UAPH model had 0%. The likely source of the continuity error for the Euler equation model are the ori ce boundary condition and the limitation of the correction factor (equation 4.32) to a certain value, which distorts the actual required air compression due to pocket vertical shrinking. However, this continuity error does no seem to a ect the nal pressure results, as a link between this error and a higher discrepancy between model and experimental results could not be established in table 5.4 and gures 5.16 to 5.19. Both models had close average percentages of the total initial air volume expelled before calculations stopped (stop criteria of an one cell air pocket or crash), which were 98.3% of the initial air for the Euler equations model and 96.5% for the UAPH. However, this value was less than 90% for three cases with the UAPH model, which was considered as an issue of early simulation stop, while Euler equation model had the index above 90% for all the cases. Also, 50% of the UAPH model simulations stopped due to crash close to the simulation end, while only 11 % of the cases for the Euler equations had this problem. However, the oscillations at the end of the simulations with the Euler equations model might have caused the second water cell to touch the crown of the pipe, causing an air pocket of one cell length and thus causing the stop criteria to be prematurely triggered. 95 P arameters Euler UA PH S (%) Q (GPM) d or if (in) Ph ysical time (s) Computational time (s) Filled Air v ol. (%) Con t. error (%) Jump P o- sition (m) Ph ysical time (s) Computational time (s) Filled Air v ol. (%) Con t. error (%) Jump P o- sition (m) 0.5 40 0.250 1 00.8 656 99.3 1.913 1.90 89.4 46 66.3 0 7.91 0.5 40 0.375 8 4.9 363 97.8 3.656 3.95 84.2 34 92.6 0 3.74 0.5 40 0.500 7 9.1 254 98.3 2.706 1.03 79.1 28 99.1 0 0.22 0.5 40 2.000 7 8.0 248 97.6 0.416 0.60 78.1 27 99.2 0 0.32 0.5 60 0.250 7 9.8 297 93.3 6.609 5.58 78.3 29 82.1 0 5.42 0.5 60 0.375 7 2.8 187 97.1 3.306 0.54 72.9 21 98.6 0 0.22 0.5 60 0.500 7 0.7 136 98.2 6.13 0.32 70.7 20 98.8 0 0.16 0.5 60 2.000 6 9.8 120 98.4 8.26 0.27 69.7 17 98.9 0 0.16 0.5 80 0.250 7 3.3 165 92.4 8.768 4.44 74.0 25 92.4 0 2.33 0.5 80 0.375 6 8.5 156 98.0 8.808 0.27 68.6 18 98.6 0 0.16 0.5 80 0.500 6 6.9 93 98.3 12.29 6 0.22 66.9 17 98.5 0 0.16 0.5 80 2.000 6 6.3 81 98.6 12.35 2 0.16 66.2 13 98.5 0 0.16 1.0 40 0.250 1 07.0 652 100.0 4.953 0.22 103 .3 57 88.1 0 3.84 1.0 40 0.375 8 9.7 417 98.7 6.141 2.17 90.0 41 96.9 0 1.95 1.0 40 0.500 8 5.1 339 98.7 5.837 1.08 84.8 35 98.9 0 0.70 1.0 40 2.000 8 3.2 308 97.6 2.053 0.65 83.4 33 98.9 0 0.76 1.0 60 0.250 8 7.5 357 98.7 7.586 2.33 87.1 43 92.6 0 2.60 1.0 60 0.375 7 5.9 245 98.9 7.786 0.32 76.4 25 99.3 0 0.16 1.0 60 0.500 7 3.9 172 98.7 5.574 0.32 73.8 21 99.2 0 0.16 1.0 60 2.000 7 2.1 149 98.8 7.977 0.27 72.0 18 99.2 0 0.16 1.0 80 0.250 7 8.2 226 99.2 11.17 1.68 79.0 29 97.7 0 0.70 1.0 80 0.375 7 1.5 195 98.6 8.68 0.27 71.7 22 98.9 0 0.16 1.0 80 0.500 6 9.8 127 98.7 11.804 0.22 69.9 19 98.9 0 0.16 1.0 80 2.000 6 8.8 110 98.9 12.723 0.16 68.7 16 98.8 0 0.16 2.0 40 0.250 1 11.1 653 100.0 4.956 0.22 110 .0 63 94.8 0 1.84 2.0 40 0.375 9 1.6 418 98.2 6.796 0.97 92.1 41 98.5 0 0.97 2.0 40 0.500 8 5.2 317 98.4 6.551 0.92 85.6 33 99.5 0 0.16 2.0 40 2.000 8 3.6 299 98.6 3.698 0.43 83.6 35 100.0 0 0.27 2.0 60 0.250 8 9.9 372 99.1 10.585 1.41 91.2 42 97.1 0 1.14 2.0 60 0.375 7 7.3 272 99.5 11.902 0.16 78.0 28 99.4 0 0.16 2.0 60 0.500 7 5.3 196 98.9 9.017 0.27 75.4 23 99.3 0 0.16 2.0 60 2.000 7 3.6 169 99.0 11.426 0.22 73.5 21 99.3 0 0.16 2.0 80 0.250 8 1.2 415 97.4 12.85 1.14 83.0 34 98.5 0 0.49 2.0 80 0.375 7 3.9 239 98.8 10.599 0.27 74.3 24 99.1 0 0.16 2.0 80 0.500 7 1.2 145 98.9 15.379 0.22 71.4 20 99.0 0 0.16 2.0 80 2.000 7 0.6 131 99.1 15.669 0.16 70.5 18 99.0 0 0.16 T able 5.4: Summary of results of b oth m o dels 96 5.3.2 Model comparison with actual pipeline lling event The comparisons between the eld data and the numerical predictions for the lling of CAESB ductile iron pipeline are presented below in Figures 5.22 to 5.24. This 350 mm diameter transmission main has a pump station, and the lling process occurs in two steps. In the rst step, the initial 4.4-km extension line is lled by gravity, throttling the upstream butter y valve so that the in ow rate is limited to Q = 0:18. In the second step, pumps are turned on and the remainder 2.8 km of the pipeline is lled. The analysis below focuses in simulation of the initial 1,700 meters of the gravity lling. The air valves positioned at x=400 m correspond to a couple of 50-mm, spherical shutter air release valves. The actual discharge area of these valves was not measured, and was calibrated in the numerical model so that the observed air pressure at the discharge point was approximately similar to the measurements. The pipeline pro le is shown 5.21, and the assumed values for the Manning roughness was n = 0:011 and for the celerity was 100 m/s. While the anticipated celerity is much larger, the adopted value is adequate considering that the modeling is not focusing on transient pressure issues but instead on pipeline lling. Moreover, the larger celerity helped limit the computational e ort for the simulations. The work by [Vasconcelos et al., 2009b] applied the TPA model that did not incorporate e ects of air pressurization to simulate pipeline lling. Figure 5.22 presents a comparison of the pressure head hydrograph measured downstream from the pump station (x 380m), and the sample frequency was 4 seconds. One notices that the eld measurement signal an increase in the pressure head at about 200 seconds into the simulation and attains a stable level. At about t > 1100 s the pressure begins to steadily rise again and will arrive at 8 m when t > 2700 s. The results obtained with the traditional TPA model indicate a small pressure (corresponding to the water depth) until about t = 2300 s, when the pressure rapidly climbs achieving levels over 7.5 m after t = 3500 s. The results obtained with both the proposed model better approach the eld measurements. Pressure begins to climb at 97 t=500 seconds, and will arrive at 8.0 m for T=2900 seconds (Euler equation model). The UAPH model presents fairly good agreement too, but at t = 1700 s it begins to diverge from the solution obtained with the Euler equation, and pressure will attain the 8.0 m only for t = 3500 s. An analogous comparison, this time however focusing on the measured and predicted in ow rate admitted into the pipeline, is presented in Figure 5.24. Flow measurements in the water main were performed with an electromagnetic ow meter, with a sampling frequency of 1 minute. The butter y valve opening was gradual, and took approximately 4 minutes. The simulation performed with the traditional TPA model (presented in [Vasconcelos et al., 2009b]) reproduced this gradual opening; the results presented here have skipped this for simplicity, assuming the nal opening right on the onset of simulation. Flow measurement indicate an initial ow rate slightly above 40 L/s, which will start declining for t> 1600 s, stabilizing in 29 L/s when t=2700 seconds. Assuming that the ow rate drop is caused by air pressurization (as in the case of the experiments performed in this study), there seems to be a slight inconsistency with the pressure measurements which indicate that pressure begins to climb when for t> 1200s. The cause for this possible inconsistency is not determined. The numerical prediction by the TPA model indicate that the ow rate drop will occur much later, whereas the proposed model indicate the ow rate drop occurring much sooner, as soon as air pressure begins to climb in the pipeline. 98 Figure 5.21: Real water main?s sketch 0 2 4 6 8 10 0 500 1000 1500 2000 2500 3000 3500 4000 head (m) time (s) Measured Traditional TPA UAPH Euler equations Figure 5.22: Field measurements and predicted heads at the upstream ventilation valve of the water main 99 0 5 10 15 20 25 0 500 1000 1500 2000 2500 3000 3500 4000 head (m) time (s) Measured Traditional TPA UAPH Euler equations Figure 5.23: Field measurements and predicted heads at the downstream ventilation valve of the water main 0 0.01 0.02 0.03 0.04 0.05 0.06 0 500 1000 1500 2000 2500 3000 3500 4000 o w rate (m 3=s ) time (s) Measured Traditional TPA UAPH Euler Figure 5.24: Field measurements and ow rates at the upstream ventilation valve of the water main 100 5.4 Conclusion This work presented two model alternatives for the modeling of the lling of a water main, together with experimental investigations on the subject with a scale model. Several numerical models have been developed to simulate this problem, however each of them with certain limitations which prevent most of them to be used for practical applications. Also, published data for the expelling of air pockets through an air valve during a re lling process is limited. This work presented results from a proposed model that couples the Saint-Venant equa- tions and the Euler equations to simulate the lling of water mains accounting for air pres- surization. The idea was to better simulate pressure gradients in the air phase, and with this be able to predict the lling events with greater accuracy. While examples presented for the rst version of the model were not a precise depiction of the conditions anticipated during water main lling, it served to illustrate the importance of incorporating air pressurization in the computations of water main lling with limited ventilation. Also, it was noted that for the simulation of this type of problem, especially when interface breakdown occurs, the model must account variations in the cross sectional area of the air ow. With this, the second version of the models incorporated a more sophisticated mathe- matical model based the Euler equations for the air phase, allowing for variable cross sectional area as well as for air frictional losses. In the second version, both presented numerical mod- els (Euler equations and UAPH approaches) successfully predicted the pressures heads and ow rates for laboratory and eld data, which showed both model?s capacity of simulating typical engineering actual conditions. Both models were also able to predict the occurrence of di erent ow features such as interface breakdown and interaction between bores and depression waves. Also, a comparison between both approaches with actual experimental results showed that, for the range of tested cases in the laboratory as well as for the actual water main case, the UAPH approach can be used instead of the Euler equations model 101 with no signi cant loss of precision and with a gain in implementation simplicity and com- putational e ort. Still, the UAPH model had more crashes and early stops than the Euler equations model, which may prevent certain cases to be properly simulated with UAPH. Im- provements still need to be done to account better for the occurrence of interface breakdown and drown inlets, which still cause non physical results. Also functionalities to account for air pocket movement (dragging and otation) need to be incorporated. An experimental investigation was performed varying common actual operation param- eters, which are the inlet ow rate, pipe slope, and ventilation size. With the results of this investigation, data was provided to calibrate/assess the proposed models and also other researcher?s models still to be developed. A comprehensive analysis of what was observed in the experimental data was also included, clarifying events and features of this type of ow, such as interface breakdowns. Experimental and numerical results show the importance that ventilation design has on the maximum pressures observed in a system. A considerable increase in the maximum pressure in the system together with a increase in the lling time was observed if ventilation is smaller than adequate. Experiments also showed the sequence of events leading to an interface breakdown in a supercritical lling event, which is: 1. Air pressure increase due to the formation of an air pocket; 2. In ow rate decreases generating a depression wave; and 3. This depression wave intercepts the pressuring bore and, if strong enough, causes the interface breakdown to happen. The two main scienti c contributions of this work were: 1. Development of a feasible modeling framework to simulate the lling of water mains accounting for air e ects; and 102 2. Demonstration that lumped approach for air phase modeling (UAPH) has comparable accuracy with discretized model based on Euler equations at a much reduced compu- tational e ort. The suggested follow ups working on this topic are: Modi cation on the modeling frameworks to account for air pocket otation and drag- ging; More generic model implementation to make it possible the simulation more realistic problems with complex geometries and several air pockets; Consideration of a wider range of air pocket formation mechanisms; and Development of wider range of anticipated boundary conditions in water main lling events; and Development of a coupled numerical solution for Euler and Saint Venant equations, possibly based on the HLL Riemann solver. 103 Bibliography [Akan, 2006] Akan, A. (2006). Open channel hydraulics. Butterworth-Heinemann. [Arai and Yamamoto, 2003] Arai, K. and Yamamoto, K. (2003). Transient analysis of mixed free-surface-pressurized ows with modi ed slot model 1: Computational model and ex- periment. In Proc. FEDSM03 4th ASME-JSME Joint Fluids Engrg. Conf. Honolulu, Hawaii, Paper, volume 45266. [Benjamin, 1968] Benjamin, T. B. (1968). Gravity currents and related phenomena. Journal of Fluid Mechanics, 31(02):209{248. [Bentley, 2010] Bentley, S. (2010). Hammer. [Bourdarias and Gerbi, 2007] Bourdarias, C. and Gerbi, S. (2007). A nite volume scheme for a model coupling free surface and pressurised ows in pipes. Journal of Computational and Applied Mathematics, 209(1):109{131. [Capart et al., 2003] Capart, H., Eldho, T., Huang, S., Young, D., and Zech, Y. (2003). Treatment of natural geometry in nite volume river ow computations. Journal of Hy- draulic Engineering, 129:385. [Capart et al., 1997] Capart, H., Sillen, X., and Zech, Y. (1997). Numerical and experimen- tal water transients in sewer pipes. Journal of hydraulic research, 35(5):659{672. [Chaiko and Brinckman, 2002] Chaiko, M. and Brinckman, K. (2002). Models for analysis of water hammer in piping with entrapped air. Journal of uids engineering, 124:194. [Chaudhry, 2008] Chaudhry, M. (2008). Open-channel ow. Springer Verlag. [Cunge and Wegner, 1964] Cunge, J. and Wegner, M. (1964). Int egration num erique des equations d? ecoulement de barr e de saint-venant par un sch ema implicite de di erences nies. La Houille Blanche, (1):33{39. [Cunge et al., 1980] Cunge, J. A. et al. (1980). Practical Aspects of Computational River Hydraulics. Pitman Advanced Pub. Program. [De Martino et al., 2008] De Martino, G., Fontana, N., and Giugni, M. (2008). Tran- sient ow caused by air expulsion through an ori ce. Journal of Hydraulic Engineering, 134:1395. [Fabre and Lin e, 1992] Fabre, J. and Lin e, A. (1992). Modeling of two-phase slug ow. Annual review of uid mechanics, 24(1):21{46. 104 [Falvey, 1980] Falvey, H. (1980). Air-water ow in hydraulic structures. NASA STI/Recon Technical Report N, 81:26429. [Fuamba, 2002] Fuamba, M. (2002). Contribution on transient ow modelling in storm sew- ers contribution sur la mod elisation de l^a[U+0080][U+0099] ecoulement transitoire dans les collecteurs. Journal of hydraulic research, 40(6):685. [Fuertes et al., 2000] Fuertes, V., Arregui, F., Cabrera, E., and Iglesias, P. (2000). Experi- mental setup of entrapped air pockets model validation. In BHR GROUP CONFERENCE SERIES PUBLICATION, volume 39, pages 133{146. Bury St. Edmunds; Professional En- gineering Publishing; 1998. [Godunov, 1959] Godunov, S. (1959). A di erence method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. Matematicheskii Sbornik, 89(3):271{306. [G omez and Achiaga, 2001] G omez, M. and Achiaga, V. (2001). Mixed ow modelling pro- duced by pressure fronts from upstream and downstream extremes. In Urban Drainage Modeling, pages 461{470. ASCE. [Guinot, 2003] Guinot, V. (2003). Godunov-type schemes: an introduction for engineers, volume 1. Elsevier Science. [Hamam and McCorquodale, 1982] Hamam, M. and McCorquodale, J. (1982). Transient conditions in the transition from gravity to surcharged sewer ow. Canadian Journal of Civil Engineering, 9(2):189{196. [Harten et al., 1983] Harten, A., Lax, P., and Van Leer, B. (1983). On upstream di erencing and godunov-type schemes for hyperbolic conservation laws. SIAM review, pages 35{61. [Issa and Kempf, 2003] Issa, R. and Kempf, M. (2003). Simulation of slug ow in horizontal and nearly horizontal pipes with the two- uid model. International journal of multiphase ow, 29(1):69{95. [Izquierdo et al., 1999] Izquierdo, J., Fuertes, V., Cabrera, E., Iglesias, P., and Garcia-Serra, J. (1999). Pipeline start-up with entrapped air. Journal of hydraulic research, 37(5):579{ 590. [Kabiri Samani et al., 2006] Kabiri Samani, A., Borghei, S., and SAEIDI, M. (2006). En- trapped air in long water tunnels during transition from a pressurized to free-surface ow regime. Scientia Iranica, 13(2):174{186. [Kalinske and Bliss, 1943] Kalinske, A. and Bliss, P. (1943). Removal of air from pipe lines by owing water. Proceedings of the American Society of Civil Engineers (ASCE), 13(10):3. [Kerger et al., 2011] Kerger, F., Erpicum, S., Dewals, B., Archambeau, P., and Pirotton, M. (2011). 1d uni ed mathematical model for environmental ow applicated to aerated mixed ows. Advances in Engineering Software, 42(9):660{670. 105 [Kr onig, 1856] Kr onig, A. (1856). Grundz uge einer theorie der gase. Annalen der Physik, 175(10):315{322. [KYPipe, 2010] KYPipe (2010). Pipe2010: Surge. [Le on et al., 2009] Le on, A., Ghidaoui, M., Schmidt, A., and Garc a, M. (2009). Application of godunov-type schemes to transient mixed ows. J. Hydraulic Res, 47(2):147{156. [Le on et al., 2010] Le on, A., Ghidaoui, M., Schmidt, A., and Garc a, M. (2010). A robust two-equation model for transient-mixed ows. [Le on et al., 2006] Le on, A., Ghidaoui, M., Schmidt, A., Garcia, M., et al. (2006). Godunov- type solutions for transient ows in sewers. Journal of Hydraulic Engineering, 132:800. [Le on et al., 2008] Le on, A., Ghidaoui, M., Schmidt, A., Garc a, M., et al. (2008). E cient second-order accurate shock-capturing scheme for modeling one-and two-phase water ham- mer ows. Journal of Hydraulic Engineering, 134:970. [LeVeque, 1992] LeVeque, R. (1992). Numerical methods for conservation laws. Birkhauser. [Li and McCorquodale, 1999] Li, J. and McCorquodale, A. (1999). Modeling mixed ow in storm sewers. Journal of Hydraulic Engineering, 125:1170. [Liou and Hunt, 1996] Liou, C. P. and Hunt, W. A. (1996). Filling of pipelines with undu- lating elevation pro les. J. Hydr. Engrg., 122(10):534{539. [Macchione and Morelli, 2003] Macchione, F. and Morelli, M. (2003). Practical aspects in comparing shock-capturing schemes for dam break problems. Journal of Hydraulic Engi- neering, 129:187. [Malekpour et al., 2011] Malekpour, A., BW, K., et al. (2011). Rapid lling analysis of pipelines with undulating pro les by the method of characteristics. ISRN Applied Math- ematics, 2011. [Martin, 1976] Martin, C. S. (1976). Entrapped air in pipelines. Second BHRA International conference on pressure surges. [McCorquodale and Hamam, 1983] McCorquodale, J. and Hamam, M. (1983). Modeling surcharged ow in sewers. In Proc., Int?l Symposium on Urban Hydrology, Hydraulics and Sediment Control, University of Kentucky, Lexington, Kentucky, pages 331{338. [McInnis and Karney, 1995] McInnis, D. and Karney, B. W. (1995). Transients in distribu- tion networks: Field tests and demand models. Journal of Hydraulic Engineering. [Politano et al., 2007] Politano, M., Odgaard, A., and Klecan, W. (2007). Case study: Nu- merical evaluation of hydraulic transients in a combined sewer over ow tunnel system. Journal of Hydraulic Engineering, 133:1103. 106 [Pothof and Clemens, 2010] Pothof, I. and Clemens, F. (2010). Experimental study of air- water ow in downward sloping pipes air-water ow in downward sloping pipes. Interna- tional Journal of Multiphase Flow. [Potter, 2004] Potter, M. C. (2004). Mec^anica dos Fluidos. Thomson, 3 edition. [Pozos et al., 2010] Pozos, O., Sanchez, A., Rodal, E., and Fairuzov, Y. (2010). E ects of water-air mixtures on hydraulic transients. Canadian Journal of Civil Engineering, 37(9):1189{1200. [Preissmann, 1961] Preissmann, A. (1961). Propagation des intumescences dans les canaux et rivieres. In First Congress of the French Association for Computation, Grenoble, France, pages 433{442. [Press, 1989] Press, W. (1989). Numerical recipes in Pascal: the art of scienti c computing. Cambridge Univ Pr. [Pulliam, 1994] Pulliam, T. H. (1994). The euler equations. Technical report, NASA Ames Research Center. [Roe, 1981] Roe, P. (1981). Approximate riemann solvers, parameter vectors, and di erence schemes. Journal of computational physics, 43(2):357{372. [Sanders et al., 2011] Sanders, B., Bradford, S., et al. (2011). Network implementation of the two-component pressure approach for transient ow in storm sewers. Journal of Hydraulic Engineering, 137:158. [Sj oberg and CTH., 1982] Sj oberg, A. and CTH. (1982). The Sewer Network Models DAGVL-A and DAGVL-DIFF. Water Resource Publications. [Song et al., 1983] Song, C., Cardle, J., and Leung, K. (1983). Transient mixed- ow models for storm sewers. Journal of hydraulic engineering, 109(11):1487{1504. [Sturm, 2001] Sturm, T. W. (2001). Open Channel Hydraulics. McGraw-Hill, 1 edition. [Toro, 2009] Toro, E. (2009). Riemann solvers and numerical methods for uid dynamics: a practical introduction. Springer Verlag. [Toro et al., 1994] Toro, E., Spruce, M., and Speares, W. (1994). Restoration of the contact surface in the hll-riemann solver. Shock waves, 4(1):25{34. [Toro, 2001] Toro, E. F. (2001). Shock-Capturing Methods for Free-Surface Shallow Flows. WILEY. [Trajkovic et al., 1999] Trajkovic, B., Ivetic, M., Calomino, F., and D?Ippolito, A. (1999). Investigation of transition from free surface to pressurized ow in a circular pipe. Water science and technology, 39(9):105{112. [Tran, 2011] Tran, P. (2011). Propagation of pressure waves in two-component bubbly ow in horizontal pipes. Journal of Hydraulic Engineering, 137:668. 107 [Trindade and Vasconcelos, 2011] Trindade, B. and Vasconcelos, J. (2011). Numerical sim- ulation of water pipeline lling events with limited ventilation. ASCE. [Vasconcelos, 2005] Vasconcelos, J. (2005). Dynamic approach to the description of ow regime transition in stormwater systems. PhD thesis, University of Michigan. [Vasconcelos and Wright, 2004] Vasconcelos, J. and Wright, S. (2004). Numerical modeling of the transition between free surface and pressurized ow in storm sewers. Innovative modeling of urban water systems, Monograph, 12:189{214. [Vasconcelos and Wright, 2009] Vasconcelos, J. and Wright, S. (2009). Investigation of rapid lling of poorly ventilated stormwater storage tunnels. Journal of hydraulic research, 47(5):547{558. [Vasconcelos et al., 2006] Vasconcelos, J., Wright, S., and Roe, P. (2006). Improved simula- tion of ow regime transition in sewers: Two-component pressure approach. Journal of Hydraulic Engineering, 132:553. [Vasconcelos et al., 2009a] Vasconcelos, J., Wright, S., and Roe, P. (2009a). Numerical os- cillations in pipe- lling bore predictions by shock-capturing models. Journal of Hydraulic Engineering, 135:296. [Vasconcelos, 2007] Vasconcelos, J. G. (2007). M etodos Num ericos Aplicados a Escoamentos Hidr aulicos. Faculdade de Tecnologia. [Vasconcelos et al., 2009b] Vasconcelos, J. G., Moraes, J. R. S., and Gebrim, D. V. B. (2009b). Field measurements and numerical modeling of a water pipeline lling events. In Event Proc. 33rd IAHR Congress. [Vasconcelos and Trindade, 2011] Vasconcelos, J. G. and Trindade, B. C. (2011). Discretized air phase modeling in the simulation of pipeline re lling operations. [Vasconcelos and Wright, 2005] Vasconcelos, J. G. and Wright, S. J. (2005). Experimental investigation of surges in a stormwater storage tunnel. Journal of Hydraulic Engineering, 131(10):853{861. [Vasconcelos, 2008] Vasconcelos, J. G., M. D. T. B. (2008). Applicability of transient ow models to water mains subjected to ow regime transition. In Proc. XXXI Interamerican Congress in Sanitary Engineering, Santiago, Chile. [Wiggert, 1972] Wiggert, D. (1972). Transient ow in free-surface, pressurized systems. Journal of the Hydraulics division, 98(1):11{27. [WISNER and Bucarest, 1967] WISNER, P. and Bucarest, R. (1967). Air demand and pul- satory pressures in bottom outlets. Rapport sur la r eunion, 1:46. [Wylie and Streeter, 1993] Wylie, E. B. and Streeter, V. L. (1993). Fluid Tansients. McGraw-Hill, 9 edition. 108 [Zhou, 2000] Zhou, F. (2000). E ects of trapped air on ow transients in rapidly lling sewers/{by Fayi Zhou. PhD thesis, University of Alberta. [Zhou et al., 2002a] Zhou, F. et al. (2002a). Transient ow in a rapidly lling horizontal pipe containing trapped air. Journal of Hydraulic Engineering, 128(6). [Zhou et al., 2002b] Zhou, F., Hicks, F. E., and Ste er, P. M. (2002b). Observations of air{ water interaction in a rapidly lling horizontal pipe. Journal of Hydraulic Engineering, 128(6):635{639. [Zhou et al., 2011] Zhou, L., Liu, D., and Zhang, Q. (2011). The in uence of entrapped air pockets on hydraulic transients in water pipelines. Journal of Hydraulic Engineering, 1:264. 109