JET 21 ANALYTICAL AND NUMERICAL ANALYSIS OF TRAPPED AIR POCKET DYNAMIC RESPONSE DUE TO PRESSURE CHANGE IN LIQUID-FILLED PIPELINES ANALITIČNA IN NUMERIČNA ANALIZA DINAMIČNEGA ODZIVA ZRAČNEGA MEHURJA NA TLAČNO MOTNJO V CEVOVODIH, NAPOLNJENIH Z VODO Jošt Pekolj 1 , Anton Bergant 1,2ℜ , Matjaž Perpar 1 Keywords: pipelines, air pocket, transient pressures, analytical solution, method of characteris- tics Abstract This paper serves as an overview of the different approaches to modelling the filling and empty- ing of liquid-filled pipelines. Specifically, it compares an analytical and numerical analysis of the start-up of a liquid column in a one-dimensional pipeline with a trapped air pocket. The analytical analysis is based on the premise that there are no pressure waves travelling in the water column. For this reason, it is not suitable for pipeline systems with very small initial trapped air pockets JET Volume 15 (2022) p.p. 21-30 Issue 2, 2022 Type of article: 1.01 www.fe.um.si/si/jet.html ℜ Corresponding author: Anton Bergant, PhD, Full time: Litostroj Power d.o.o., Litostrojska 50, 1000 Ljubljana, Slo- venia, anton.bergant@litostrojpower.eu; Part time: Faculty of Mechanical Engineering, University of Ljubljana, Ašk- erčeva 6, 1000 Ljubljana, Slovenia 1 Faculty of Mechanical Engineering, University of Ljubljana, Aškerčeva 6, 1000 Ljubljana, Slovenia 2 Litostroj Power d.o.o., Litostrojska 50, 1000 Ljubljana, Slovenia Analytical and numerical analysis of trapped air pocket dynamic response due to pressure change in liquid-filled pipe- lines Analitična in numerična analiza dinamičnega odziva zračnega mehurja na tlačno motnjo v cevovodih, napolnjenih z vodo Jošt Pekolj, Anton Bergant, Matjaž Perpar 22 JET 22 JET Jošt Pekolj, Anton Bergant, Matjaž Perpar JET Vol. 15 (2022) Issue 2 or very high reservoir pressures. In the pipelines with parameters, as shown in this paper, the results of the analytical analysis correlate sufficiently with those of the numerical analysis. Povzetek Prispevek predstavlja pregled dveh različnih pristopov k modeliranju polnjenja in praznjenja cevovodov s kapljevino. Primerja analitično in numerično analizo polnjenja v enodimenzijskem cevovodu z ujetim zračnim mehurjem na enem koncu. Analitična analiza temelji na predpostav- ki, da po vodnem stolpcu ne potujejo tlačni valovi. Zaradi tega ta pristop ni primeren za cevne sisteme z zelo majhnimi začetnimi zračnimi mehurji ali zelo visokimi tlaki rezervoarja. Za predsta- vljene primere v tem delu se analitični in numerični rezultati ujemajo. 1 INTRODUCTION During filling or emptying of liquid-filled pipelines an air pocket, or pockets, of various sizes can get trapped at distinct positions (high point, pipe end). The presence of an air pocket influences the system’s response to a pressure disturbance. Pressure in an air pocket could increase by multiples of the original pressure change. The response of an air pocket can be simulated using different computational models. In this paper, the authors examine the one-dimensional analyt- ical and numerical models of air pocket dynamic response. Air pockets become trapped at the downstream end of the gravitational type of pipeline. The basis of the analytical analysis is to describe the behaviour of the system based on differential equations. The method of character- istics will be used for the numerical analysis, as it is convenient for tracking pressure waves in the water filled region of the pipe. The analytical and numerical results will be compared for different parameters of the hydraulic pipeline system during partial filling, as shown in Figure 1. At the leftmost part of the pipeline there is a reservoir with a pressure head H rez , which is attached to the pipeline with a water column of length L(t) and an air pocket of length L gas (t). The water and air regions are initially separated by a valve, which is located at x 1 . After the rapid valve opening, the water column moves towards the air pocket at velocity v(t). The response of an air pocket in similar pipelines has also been studied experimentally in [1], [2] and [3]. When comparing the two types of analysis, the numerical one shows better agreement with the experimental results, because it also includes the interaction of pressure waves and unsteady friction [2]. Figure 1: A hydraulic pipeline system with a trapped air pocket JET 23 Analytical and numerical analysis of trapped air pocket dynamic response due to pressure change in liquid-filled pipelines 2 THEORETICAL MODELLING The following subchapters describe the most important equations used in analytical and numer- ical analysis. The equations are valid for the system shown in Figure 1. 2.1 Water region The water column in the water region of the pipe can be modelled as rigid or elastic. The one-di- mensional motion of the rigid water column is governed by the momentum conservation, as shown in equation (2.1): dv dt ( t ) = g ∙ ( H a , r ez − H a , gas ( t ) ) L ( t ) + g ∙ sin θ − ¿ − λ 2 D ∙ v ( t ) | v ( t ) | − K v 2 L ( t ) ∙ v ( t ) | v ( t ) | − K e + 1 2 L ( t ) ∙ v 2 ( t ) ∙ H ( v ( t ) ) . (2.1) dv dt ( t ) = g ∙ ( H a , r ez − H a , gas ( t ) ) L ( t ) + g ∙ sin θ − ¿ − λ 2 D ∙ v ( t ) | v ( t ) | − K v 2 L ( t ) ∙ v ( t ) | v ( t ) | − K e + 1 2 L ( t ) ∙ v 2 ( t ) ∙ H ( v ( t ) ) . (2.1) dv dt ( t ) = g ∙ ( H a , r ez − H a , gas ( t ) ) L ( t ) + g ∙ sin θ − ¿ − λ 2 D ∙ v ( t ) | v ( t ) | − K v 2 L ( t ) ∙ v ( t ) | v ( t ) | − K e + 1 2 L ( t ) ∙ v 2 ( t ) ∙ H ( v ( t ) ) . (2.1) dv dt ( t ) = g ∙ ( H a , r ez − H a , gas ( t ) ) L ( t ) + g ∙ sin θ − ¿ − λ 2 D ∙ v ( t ) | v ( t ) | − K v 2 L ( t ) ∙ v ( t ) | v ( t ) | − K e + 1 2 L ( t ) ∙ v 2 ( t ) ∙ H ( v ( t ) ) . (2.1) (2.1) Similar equations can be seen in [4] and [5]. Equation (2.1) is simplified by neglecting some of the terms; namely, friction factor λ, coefficients of valve head losses K v and reservoir exit losses K e , while the slope term of the pipe will be set to zero (horizontal pipe). For the elastic water region, it is necessary to consider the motion of pressure waves in the water. This is based on the assumption that the waves are one-dimensional, which is justified, since the diameter of the pipe is substantially smaller than its length. The motion of the water is governed by the continuity and momentum conservation, as shown in equations (2.2) and (2.3) [6]: ∂ H ∂ t + v ∙ ∂ H ∂ x − v ∙ s in θ + a 2 g ∂ v ∂ x = 0, ( 2. 2) g ∙ ∂ H ∂ x + ∂ v ∂ t + v ∙ ∂ v ∂ x + λ ∙ v ∙ | v | 2 D = 0 . ( 2. 3) H i = C p B m + C m B p B p + B m , ( 2. 4) Q i = C p − C m B p + B m . ( 2. 5) H a , gas ( t ) ∙ L gas n ( t ) = C A . ( 2. 6) (2.2) ∂ H ∂ t + v ∙ ∂ H ∂ x − v ∙ s in θ + a 2 g ∂ v ∂ x = 0, ( 2. 2) g ∙ ∂ H ∂ x + ∂ v ∂ t + v ∙ ∂ v ∂ x + λ ∙ v ∙| v | 2 D = 0 . ( 2. 3) H i = C p B m + C m B p B p + B m , ( 2. 4) Q i = C p − C m B p + B m . ( 2. 5) H a , gas ( t ) ∙ L gas n ( t ) = C A . ( 2. 6) (2.3) It is assumed that the velocity of pressure waves a is constant. Thereafter, the method of charac- teristics is used to transform equations (2.2) and (2.3) into four ordinary differential equations, which are solved using the finite difference method [6]. They are combined to derive a solution for the piezometric head H i and volume flow rate Q i at a node in a characteristic grid, as shown in equations (2.4) and (2.5), where the friction head losses are neglected: ∂ H ∂ t + v ∙ ∂ H ∂ x − v ∙ s in θ + a 2 g ∂ v ∂ x = 0, ( 2. 2) g ∙ ∂ H ∂ x + ∂ v ∂ t + v ∙ ∂ v ∂ x + λ ∙ v ∙ | v | 2 D = 0 . ( 2. 3) H i = C p B m + C m B p B p + B m , ( 2. 4) Q i = C p − C m B p + B m . ( 2. 5) H a , gas ( t ) ∙ L gas n ( t ) = C A . ( 2. 6) (2.4) ∂ H ∂ t + v ∙ ∂ H ∂ x − v ∙ s in θ + a 2 g ∂ v ∂ x = 0, ( 2. 2) g ∙ ∂ H ∂ x + ∂ v ∂ t + v ∙ ∂ v ∂ x + λ ∙ v ∙ | v | 2 D = 0 . ( 2. 3) H i = C p B m + C m B p B p + B m , ( 2. 4) Q i = C p − C m B p + B m . ( 2. 5) H a , gas ( t ) ∙ L gas n ( t ) = C A . ( 2. 6) (2.5) In order to calculate variables H i and Q i at the boundary, a device-specific equation, or equations, replaces one of the water hammer compatibility equations (the reservoir and air pocket in this case) [2], [6]. 24 JET 24 JET Jošt Pekolj, Anton Bergant, Matjaž Perpar JET Vol. 15 (2022) Issue 2 2.2 Air region The air pocket is modelled as a lumped body. Therefore, the presence of pressure waves is not included in the air region in the model. A similar description is illustrated in [6]. During the filling of the pipeline, the water column compresses the air pocket, which undergoes a polytropic pro- cess. The changes in pressure and volume of the air pocket are described using equation (2.6): ∂ H ∂ t + v ∙ ∂ H ∂ x − v ∙ s in θ + a 2 g ∂ v ∂ x = 0, ( 2. 2) g ∙ ∂ H ∂ x + ∂ v ∂ t + v ∙ ∂ v ∂ x + λ ∙ v ∙ | v | 2 D = 0 . ( 2. 3) H i = C p B m + C m B p B p + B m , ( 2. 4) Q i = C p − C m B p + B m . ( 2. 5) H a , gas ( t ) ∙ L gas n ( t ) = C A . ( 2. 6) (2.6) The gas constant CA is calculated by using the known initial conditions (pressure head and air pocket length) in the equation (2.6). From previous experimental investigations, such as [1-3], it was determined that the numerical model most closely matches the experimental data when the polytropic coefficient for the process in air equals 1.4. 2.3 Gas-water interface When using the rigid water column approach, equations (2.1) and (2.6) are coupled for the air and water regions. The result is a linear first-order ordinary differential equation. Equation (2.7) describes the motion of a gas-water interface [5]: d v 2 dL + ( K + 1 L + λ D ) v 2 = 2 g ( H a , r ez L − C A ( x L − L ) n + s in θ ) . ( 2. 7) v = L − ( K + 1 2 ) ∙ √ 2 g ∙ √ H a , r ez ∫ L 0 L L ¿ K ∙ d L ¿ − C A ∫ L 0 L L ¿ K ( x L − L ¿ ) n ∙ d L ¿ + ¿ + sin θ ∫ L 0 L L ¿ K + 1 ∙ d L ¿ . ( 2. 8) (2.7) It is assumed that the gas-water interface is always vertical (both in the rigid and elastic water approach). Equation (2.7) can be solved analytically (in this case) or numerically [5]. By using the elastic water column approach, the interface changes with the volume of the air pocket. 3 ANALYTICAL AND NUMERICAL SOLUTIONS 3.1 Analytical solution The analytical method solves the differential equation for the gas-water interface (2.7) [5]. This equation was solved by using an integration factor. The solution is an expression for the speed of the water column in terms of its length (2.8): Analytical and numerical analysis of trapped air pocket dynamic response due to pressure change in liquid-filled pipelines 5 ---------- 𝑣𝑣 = 𝐿𝐿 � � ��� � � ∙ � 2𝑔𝑔 ∙ � 𝐻𝐻 �,��� ∫ 𝐿𝐿 ∗ � ∙ 𝑑𝑑𝐿𝐿 ∗ � � � − 𝐶𝐶 � ∫ � ∗ � (� � �� ∗ ) � ∙ 𝑑𝑑𝐿𝐿 ∗ � � � + + sin 𝜃𝜃 ∫ 𝐿𝐿 ∗ ��� ∙ 𝑑𝑑𝐿𝐿 ∗ � � � . (2.8) With the help of equation (2.6) for the air region, it was also possible to calculate how the pressure in the air pocket changes over time. 3.2 Numerical solution For the numerical analysis, the numerical model based on the method of characteristics [6] was used. The nodes in the numerical grid are fixed and are distributed over the length of the initial water column – from the leftmost node at the reservoir to the rightmost node at the valve. For the interior nodes of the water region characteristic grid, the piezometric head and volume flow rate are calculated using (2.4) and (2.5). At the leftmost node, at the boundary with the reservoir, the modified characteristic equations (2.4) and (2.5) are coupled with the boundary equation of the constant head reservoir [6]. For the node at the gas-water interface, the air region equation coupled with the modified characteristic equations was used. When the water-air interface moves, it is assumed that the variables Hi and Qi at the interface and at the fixed node on the valve are the same. 4 COMPARISONS OF ANALYTICAL AND NUMERICAL RESULTS The analytical and numerical solutions to the equations were calculated using different parameters of the hydraulic pipeline system, as depicted in Figure 1. The parameters for the different pipelines are shown in Tables 1, 2 and 3. For all of the cases illustrated, the pressure loss coefficient K and slope of the pipe θ are equal to zero. Table 1: Pipeline system parameters for Case 1 [4], [5] Hrez [m] Ha, gas0 [m] L0 [m] Lgas0 [m] D [m] n [-] ρv [kg/m 3 ] a [m/s] 31 10.3 100 15 0.3 1.4 1000 1000 Table 2: Pipeline system parameters for Case 2 [1] Hrez [m] Ha, gas0 [m] L0 [m] Lgas0 [m] D [m] n [-] ρv [kg/m 3 ] a [m/s] 8.15 10.32 5.57 3.25 0.04 1.4 1000 400 Table 3: Pipeline system parameters for Case 3 [1] Hrez [m] Ha, gas0 [m] L0 [m] Lgas0 [m] D [m] n [-] ρv [kg/m 3 ] a [m/s] 12.23 10.28 5.57 3.25 0.04 1.4 1000 400 (2.8) With the help of equation (2.6) for the air region, it was also possible to calculate how the pres- sure in the air pocket changes over time. JET 25 Analytical and numerical analysis of trapped air pocket dynamic response due to pressure change in liquid-filled pipelines 3.2 Numerical solution For the numerical analysis, the numerical model based on the method of characteristics [6] was used. The nodes in the numerical grid are fixed and are distributed over the length of the initial water column – from the leftmost node at the reservoir to the rightmost node at the valve. For the interior nodes of the water region characteristic grid, the piezometric head and volume flow rate are calculated using (2.4) and (2.5). At the leftmost node, at the boundary with the reservoir, the modified characteristic equations (2.4) and (2.5) are coupled with the boundary equation of the constant head reservoir [6]. For the node at the gas-water interface, the air region equa- tion coupled with the modified characteristic equations was used. When the water-air interface moves, it is assumed that the variables H i and Q i at the interface and at the fixed node on the valve are the same. 4 COMPARISONS OF ANALYTICAL AND NUMERICAL RESULTS The analytical and numerical solutions to the equations were calculated using different param- eters of the hydraulic pipeline system, as depicted in Figure 1. The parameters for the different pipelines are shown in Tables 1, 2 and 3. For all of the cases illustrated, the pressure loss coeffi- cient K and slope of the pipe θ are equal to zero. Table 1: Pipeline system parameters for Case 1 [4], [5] H rez [m] H a, gas0 [m] L 0 [m] L gas0 [m] D [m] n [-] ρ v [kg/m 3 ] a [m/s] 31 10.3 100 15 0.3 1.4 1000 1000 Table 2: Pipeline system parameters for Case 2 [1] H rez [m] H a, gas0 [m] L 0 [m] L gas0 [m] D [m] n [-] ρ v [kg/m 3 ] a [m/s] 8.15 10.32 5.57 3.25 0.04 1.4 1000 400 Table 3: Pipeline system parameters for Case 3 [1] H rez [m] H a, gas0 [m] L 0 [m] L gas0 [m] D [m] n [-] ρ v [kg/m 3 ] a [m/s] 12.23 10.28 5.57 3.25 0.04 1.4 1000 400 Figures 2 to 4 show the analytical solutions for the velocity of the water column in terms of its length for all the three cases, as illustrated in Tables 1, 2 and 3. The results shown in Figure 2 agree with those presented by Tijsseling et al. [5]. 26 JET 26 JET Jošt Pekolj, Anton Bergant, Matjaž Perpar JET Vol. 15 (2022) Issue 2 Figure 2: Water column velocity in terms of column length for Case 1 Figure 3: Water column velocity in terms of column length for Case 2 JET 27 Analytical and numerical analysis of trapped air pocket dynamic response due to pressure change in liquid-filled pipelines Figure 4: Water column velocity in terms of column length for Case 3 Figures 5 to 7 show the comparison of calculated pressure in the air pocket by using the analyt- ical and numerical methods (Chapter 3) for all three cases in Tables 1, 2 and 3. As can be seen from the figures, the air pocket pressure response for the discussed cases is similar to that of the analytical and numerical analysis. Figure 5: A comparison of analytical and numerical air pocket pressure for Case 1 28 JET 28 JET Jošt Pekolj, Anton Bergant, Matjaž Perpar JET Vol. 15 (2022) Issue 2 Figure 6: A comparison of analytical and numerical air pocket pressure for Case 2 Figure 7: A comparison of analytical and numerical air pocket pressure for Case 3 JET 29 Analytical and numerical analysis of trapped air pocket dynamic response due to pressure change in liquid-filled pipelines 5 CONCLUSION The analytical and numerical results match quite well, however there are still some differences. The reason for this is that the physical models are based on slightly different approaches, e.g. the rigid and elastic water column, and a different consideration of the water region between the valve and air pocket. The models match for those cases with large initial air pockets and low reservoir pressures. Since the analytical analysis does not include pressure waves in the water column and unsteady friction [2], it is not suitable for cases with small initial air pockets and high reservoir pressures. The numerical method could be improved by providing a better description of the gas-water interface. The part of the water column that extends over the valve will be in- cluded in the elastic water region in future research. References [1] L. Zhou, D. Liu, B. Karney, P . Wang: Phenomenon of White Mist in Pipelines Rapidly Fill- ing with Water with Entrapped Air Pockets, Journal of Hydraulic Engineering, Vol.139, Iss. 10, p.p.1041-1051, 2013 [2] A. Bergant, A. Tijsseling, Y. Kim, U. Karadžić, L. Zhou, M.F. Lambert, A.R. Simpson: Un- steady Pressures Influenced by Trapped Air Pockets in Water-Filled Pipelines, Strojniški vestnik – Journal of Mechanical Engineering, Vol.64, Iss.6, pp.501-512, 2018 [3] L. Zhou, Y. Cao, B. Karney, A. Bergant, A. S. Tijsseling, D. Liu, P.Wang: Expulsion of Entrapped Air in a Rapidly Filling Horizontal Pipe, Journal of Hydraulic Engineering, Vol.146, Iss.7, p.p.1-16, 2020 [4] E. Cabrera, J. Abreu, R. Perez, A. Vela: Influence of Liquid Length Variation in Hydraulic Transients, Journal of Hydraulic Engineering, Vol.118, Iss.12, p.p.1639-1650, 1992 [5] A. S. Tijsseling, Q. Hou, Z. Bozkuş: Analytical expressions for liquid-column velocities in pipelines with entrapped air, In: Proceedings of the ASME 2015 Pressure Vessels & Piping Division Conference, ASME, Paper PVP2015-45184, 2015 [6] F. B. Wylie, V. L. Streeter: Fluid Transients in Systems, Prentice Hall, 1993 Nomenclature A Area a Velocity of pressure waves B m Proportionality constant B p Proportionality constant C A Gas constant C m Proportionality constant C p Proportionality constant D Internal diameter of the pipe 30 JET 30 JET Jošt Pekolj, Anton Bergant, Matjaž Perpar JET Vol. 15 (2022) Issue 2 g Gravitational acceleration (g = 9.81 m/s 2 ) H Heaviside step function H a,gas Air pocket absolute pressure head H a,gas0 Air pocket initial absolute pressure head H a,rez Reservoir absolute pressure head H i Pressure head at a node H rez Reservoir pressure head K Combined pressure head loss coefficient K e Pressure head loss coefficient for reservoir exit K v Pressure head loss coefficient for the valve L Water column length L * Dummy length in integral L 0 Initial water column length L gas Air pocket length L gas0 Initial air pocket length n Polytropic coefficient Q i Volumetric flow rate at a node t Time v Water column velocity x Space coordinate in direction of the pipeline x 1 Location of the valve x L Full length of the pipeline λ Friction factor θ Slope of the pipe ρ v Water density Acknowledgments Acknowledgments The authors gratefully acknowledge the support of the Slovenian Research Agency (ARRS) con- ducted through research project L2-1825.