Abstract
On the premise of water hammer theory, a numerical model is proposed for simulating the filling process in an initially empty water conveyance pipeline with an undulating profile. Assuming that the pipeline remains full and ignoring air and water interactions in the already filled pipeline, the ongoing filling is simulated using the method of characteristics on an adaptive computational grid. The performance of the model is verified using previously published experimental and rigid column data. The model nicely replicates published experimental data. The model shows that the movement of the filling front into the system can be assumed as a rigid column as long as the flow away from the filling front is undisturbed elsewhere. Furthermore, applying the model to a hypothetical pipe system with an inline-partially open valve shows that the proposed model is robust enough to capture the transient events initiated within the moving column, a vital capability that the existing rigid water column models lack.
1. Introduction
Empty water conveyance pipelines are usually filled cautiously to prevent the onset of harmful transient pressures. Lack of sufficient experimental data forces operators to take a conservative approach to filling. As Watters [1] recommended, a safe operation requires that the filling be carried out with a discharge equivalent to that produced if the pipe were to run full with the velocity 0.3 m/s or less. This velocity criterion practically assures safe operation provided that designated air valves function properly.
Such a controlled operation over an undulating terrain creates alternating segments of open channel and pressurized flow along the line with the number and arrangement of these segments depending on both the line’s profile and the filling protocol. These segments are connected via a series of moving hydraulic jumps which are in turn responsible for slowly pushing the air out of the pipeline through installed and functioning air valves.
In practice, filling may not be carried out as smoothly as required by the dictates of good design. Poor operation caused by either human error or lack of proper control will result in rapid filling with its consequent drawbacks. Izquierdo et al. [2] show that severe overpressures may occur when the extending water column is suddenly arrested by a stagnant water pool remaining in the pipeline. Liou and Hunt [3] show that negative pressures can be induced as the moving water column passes through knees in the profile; such events might result in difficulties including the intrusion of contamination into the supply system, structural buckling of the pipe wall, or column separation with its associated high risk of excessive column rejoinder pressures.
Two important features of the flow during rapid filling are the bulk motion of the water column and any water hammer pressures waves propagated through the water column. The second feature is of great interest only once the moving water column is locally disturbed such as if the water column meets a partially open valve or a stagnant water column. However, in the absence of abrupt changes, the bulk motion of the water column dominates flow.
Depending on what flow features are of primary concern, different mathematical models are available ranging from simplified to highly complex. If the bulk motion of water is of primary interest, the rigid column model is usually the numerically cheapest solution in terms of computational time and the ease of implementation. This simplified model assumes that the whole water column accelerates and decelerates as a unit. Liou and Hunt [3] used a rigid column model to obtain the time history of the length and velocity of the moving water column and to calculate the instantaneous pressure distribution along the water column. This approach predicted the velocity history and evaluated the likelihood of column separation along the pipeline. Izquierdo et al. [2] employed a rigid column model to simulate the filling of a pumping line with an undulating elevation profile containing two static water columns separated by a large air pocket. They reported a huge pressure spike as the result of the rejoining water columns. Zhou et al. [4] applied a rigid column model to a simple system capped with an adjustable size orifice at its end to investigate how the degree of ventilation affects the transient pressures during filling.
However, if the water column is internally disturbed, a rigid column model can no longer capture the resulting water hammer pressures and a more complete dynamic model is then required to simulate the key flow features. Only a few more complete dynamic models specifically focusing on undulating pipelines have thus far been presented, though more have been proposed to cover filling in sewer systems. Sewer systems have likely received more interest and care because filling events are more frequent over a sewer’s lifetime, and thus these systems must be carefully protected from the drawbacks of filling. By contrast, filling pressures in water conveyance lines can usually be controlled through a well chosen operational strategy. Nevertheless the filling processes in sewer systems and in water conveyance pipelines are physically similar except for the dominant pace of the transient. This implies that sewer models used can be used for water conveyance applications provided that they can capture any water hammer events that occur.
Shock capturing and shock fitting methods are two available techniques that have been proposed to solve the full dynamic 1D equations (momentum and continuity equations). In typical shock capturing methods, a hypothetical slot at the top of the pipe allows simulation of both pressurized and free flow by a single set of equations. Since a single set of equations is employed in the simulation, the track of the pressurized front is automatically included. The method was originally suggested in honour of Preissmann in 1961 [5] and is thus named the Preissmann slot method or slot method for short. Garcia-Navarro et al. [6], Capart et al. [7], Ji [8], Trajkovic et al. [9], and Vasconcelos and Wright [10] among others used the slot method to simulate the flow regime transitions in some sewer systems.
However, this method possesses some important numerical limitations. First, to avoid nonlinear numerical instability, an artificially low wave velocity is required computationally. As a result, the slot method not only underestimates transient pressures (which depend strongly on the wave velocity magnitude) but also fails to properly track the timing of the pressure waves. In addition, the method cannot predict negative pressures in the pressurized zone of the real system; whenever negative pressures occur in the model, the flow is automatically switched back to open channel flow. The later limitation is addressed by Vasconcelos et al. [11] who improved the slot method by permitting a more realistic approach to negative pressures in the pressurized zone. Nevertheless, this method still suffers from the first limitation.
Shock fitting is another popular method which is extensively employed in transient analysis of the sewer systems [12–15]. In this method, the pressurization front or “the interface” is traced during each time step and the flow in both free flow and pressurized sections is treated by their corresponding theories. This method is able to maintain negative pressure and high wave speeds, but the interface instability remains an important computational challenge. Recently, Leon et al. [16] proposed a shock-fitting-based method which benefits from employing the second-order Godunov method in both its open channel and pressurized flow sections. This paper illustrates that the model can handle wave velocities up to about 1000 m/s, a significant achievement over previous shock fitting methods. Although the examples covered in Leon et al. show that interface stability is achieved even with wave velocities as high as 1000 m/s, it is not known if interface concerns might still arise if high and low water hammer pressure waves both interact with the interface. Moreover, manipulating boundary condition in the Godunov method is challenging particularly in the early stages when the first computational cell is being filled. Though not yet fully tested, this reasoning induced the authors to avoid the Leon approach in the context of the current research.
Another shock fitting family method was presented by Malekpour and Karney [17]. They used an implicit finite difference method with an expandable computational grid to simulate the rapid filling in the water conveyance systems. Although the proposed method succeeded in capturing the bulk motion of flow, it could not efficiently simulate any resulting water hammer pressures.
The first intent of the current paper is to develop and test a robust and flexible numerical model for simulating rapid filling in water conveyance pipelines having an undulating profile. The proposed model is a type of shock fitting method which applies the method of characteristics along with an expandable computational cell. As will be seen, the model robustly captures extreme water hammer spikes during the filling without stability issues. Although beyond the scope of the current work, it can be shown that this method can also capture column separation during rapid filling events. The secondary aim is to provide physical insight into the filling process through evaluating the proposed model results in comparison with results obtained from existing rigid column models.
2. Theoretical Background
2.1. Governing Equations and Assumptions
The transient flow or water hammer in the closed conduits is governed by the following two partial differential equations known as the momentum and continuity equations, respectively, [18]: where is velocity, is piezometric head, is pipe diameter, is elastic wave velocity, is gravitational acceleration, is friction factor in the Darcy-Weisbach equation, and , are distance and time independent variables, respectively.
To solve the equations in the context of rapid filling, the following assumptions are made: (1) the pipeline remains full with a well-defined vertical front during the filling, (2) the pipe system is sufficiently vented so that the air pressure remains essentially atmospheric and imposes little retarding force against the filling water column, (3) except in the immediate vicinity of the front (which is assumed incompressible), the water-pipe systems are slightly compressible everywhere, and (4) the frictional flow-resistance relationship for steady flow is a good approximation for the transient flow.
To evaluate whether the first assumption is reasonable during filling, a velocity-based criterion was proposed by Liou and Hunt [3]. They suggested that if the velocity of the moving water column exceeds a particular reference value, the pipeline remains full with a well-defined filling front. Otherwise, air will intrude into the moving water column and open channel flow is developed. The reference velocity was taken as that of an air cavity intruding along a filled horizontal pipeline with initial zero velocity. Using previous work, they suggested that, for horizontal and nearly horizontal pipes with the diameters between 4 and 18 cm, the air cavity velocity can be calculated by , where and are gravitational acceleration and pipe diameter, respectively. They also recommended that, for larger pipes with horizontal and downward angle approaching to 45°, the velocity can be computed by , and respectively. Vasconcelos and Wright [19], however, argued that, in some conditions, even once the velocity criterion is met, the shape of the front resembles a dam break front running over a dry bed with resistance rather, and not a well-defined vertical front. Although this might be an important issue if the air and water interaction is also of main concern, its distortion of the overall solution can reasonably be assumed small in the following research, as the goal is to describe essential events in the pressurized zone.
2.2. Numerical Strategy
Inspired by the traditional model used in furrow and border irrigation modeling [20, 21], Malekpour and Karney [17] proposed a numerical model based on the water hammer theory to simulate rapid pipeline filling. To handle the numerical issue analogous to the dried bed zone, they proposed an implicit finite difference with an adaptive computational mesh just like those were initially developed to simulate flow in irrigation furrows and borders. As shown in Figure 1, to represent the position of advancing front, the number of computational cells increases by one at each time step. Knowing the position of the water column front, the current time step, the discharge and piezometric heads in the intermediate and the first nodes were simultaneously calculated using an implicit finite difference method.
Although this model nicely replicates experimental data, convergence was a challenging issue. To assure convergence, the employed finite difference must be fully implicit. Yet, the resulting scheme provided accurate results only if the Courant number is kept close to 1 or when there are slower transients in systems with large effective storage, through either low wave speeds or additional storage elements in the system. Otherwise, the dissipative nature of the scheme distorts the accuracy of the results. However, in the model proposed here, the Courant numbers are enforced computationally rather than by setting a prior constant value. Since the time step is calculated based on the velocity of the bulk motion of water column () rather than the elastic wave velocity, the model has to be run with a high Courant number which in some stages may temporarily be as high as 500. This obviously makes the proposed model inefficient in capturing and tracing any water hammer pressures that might occur.
To overcome this problem, an improved version of the aforementioned model is proposed herein. In this model, the time step is fixed during the simulation and the water column front position is sought at the current time. This implies that, as shown in Figure 2, the water column traverses the front cell in several time steps rather than one. Under such a condition, any explicit method could be incorporated in modelling to make sure that no important pressures event is overlooked during the filling. The method of characteristics is used herein because of its long application in water hammer modelling, its proven accuracy, its physical compatibility, and its established superiority in handling complex boundary conditions.
2.3. The Method of Characteristics
The method of characteristics is a powerful tool for solving partial differential equations governing transient flow in closed conduits. This method transforms the governing equations into four ordinary differential equations which can be integrated through the paths called the characteristic curves to find the dependent variables at specific distances and times.
The characteristic equations for the water hammer equations are as follows [18]: These equations show that along the positive characteristic () and the negative characteristic () lines, (2.2) and (2.4), the so-called compatibility equations, govern the flow and head variations. At the intersection where the two characteristic lines meet, the independent variables and can then be calculated by simultaneous solution of the algebraic equations. The material presented here provides the basic formulation of the method of characteristic for the water hammer equations. Other issues, such as choosing the time step, stability issues, and dealing with the boundary conditions, are discussed later in the context of the current research.
Though analytical solutions are possible in trivial cases, the numerical methods are inevitable in practical systems with complex boundary conditions. The solution of the compatibility equations can be easily made based on a finite difference computational mesh depicted in Figure 3. This figure shows that the unknown variables at the current time for a given grid point () can be calculated using the information of known points and . For the sake of generality, it is assumed that the characteristic lines do not meet exactly at grid points where the information is stored. The discrete form of the compatibility equations can be then written as follows [18]: where In these equations, is discharge while and are space and time increments, respectively. As can be seen, (2.6) is linear and can be simultaneously solved to determine and at point . , , , and are linearly interpolated from the grid points information by using the following equations:
2.4. Numerical Structure and Detailed Formulation
The simulation starts with a computational grid of just one cell. Since the water column is elongating with time in this cell, it is impossible to directly use (2.1). Instead, an integrated form of the equation between two subsequent time lines can be easily derived. Assuming that the pressure at the front cells remains atmospheric during the filling and considering Figure 4(a), the discretized form of both continuity and momentum equations can be written as follows: where is discharge, is piezometric head, is cross sectional area of the pipe, is pipe diameter, is water column length in the front cell in the previous time line, is water column advancement in current time step, is computational time step, is time weighing factor, is pipeline slope, is energy slope, is gravitational acceleration, and , are time sub indexes for the previous and current time lines.
(a)
(b)
(c)
To calculate the three unknowns, , , and , in (2.9), an extra equation is needed. This equation is easily developed by considering energy conversions at the upstream boundary. Assuming that the pipeline is fed by a constant water level reservoir (see Figure 5), the boundary equation can be written as follows: where is valve head loss coefficient, is stagnant water column length, is head loss coefficient accounting for the friction loss in the stagnant water column, and is reservoir water height.
To prevent dealing with short pipe segments which can greatly increase simulation times, the water found between the valve and the reservoir is treated in the boundary equation rather than as a separate pipe. As can be seen in (2.10), the third term calculates the head required to accelerate water in this initially stagnant water column. This term is crucial during the early phases of filling because it slows down the water column by absorbing most part of the driving head of the reservoir.
No analytical solution is available for (2.9) and (2.10), so the Newton-Raphson iterative method is employed. Since in this stage the unknowns are calculated implicitly, the time step is not restricted by numerical instability issues. However, to obtain greater resolution in the early stage of the filling and also for the sake of consistency with other stages of calculations, the time step is limited through using the well-known Courant-Friedrich-Lewy (CFL) stability criterion.
For a system with a single pipe, the time step can be easily calculated by where the Courant number is set to 1. In systems with multiple pipes, each pipe may require its own time step as indicated by the following relation: where is the length of the pipe line, is wave velocity, and is an integer number showing the number of numerical divisions in the pipe.
In practice, however, it is impractical to use different time steps and care is needed to retain similar steps in each pipe. One solution is to apply the least calculated time step for all pipes. Since in this case the characteristic lines do not meet the computational grid points for many pipes, linear interpolation is required to extract information of the previous time line. The second alternative is to slightly adjust wave velocities to make all time steps the same, though the adjustment should not exceed 15% or the original value [18]. Nevertheless this later alternative has proven preferable, because the numerical solution quickly loses accuracy when large linear interpolations are used.
Once the time step is determined, the solution can be established step by step until the water column length equals or slightly exceeds the computational distance interval . When the water column length exceeds , the time increment and corresponding upstream head and discharge are linearly interpolated. It is important to know that the interpolation does not produce significant error particularly because column advancement in each time step is much shorter than the cell length so that the water column may only exceed the cell length by an infinitesimal amount.
Once the first cell has been filled, as shown in Figure 4(b), the computation advances by adding a new cell to the computational grid. In this stage, the calculation of the unknowns at the boundary point () is somewhat different from that in the last stage. To calculate the two unknowns, the head and discharge, in the upstream boundary, the energy equation (2.10) should be combined with the following equation, the negative characteristic equation, where
Simultaneous solution of (2.10) and (2.12) result in the following equation through which the discharge at the upstream boundary node can be easily calculated: where Having calculated discharge at upstream boundary, the corresponding head can be calculated at this point by substituting the discharge in (2.12).
To calculate the remaining unknowns, consider Figure 4(b). For sake of generalization, as shown in the figure, it is assumed that there would be a minor loss at the second node and that the pressure heads at the upstream and downstream side of the node are different. The unknowns in this stage are the discharge, the pressure heads at the upstream and downstream of the node, and the flow advancement length in the front cell. Considering the conservation of mass and momentum in the front cell, two equations can be easily formed through substituting the , , and in (2.9) by the , , and , respectively. The two other equations can be achieved by applying the energy equation across the node and through using the positive characteristic equations as follows: where is the head loss coefficient of the node, , and The resulting nonlinear equations are then solved by the Newton-Raphson method. The solution then proceeds to the new time level, and the same procedure is repeated until the front reaches the end of the second cell. At this time, another cell is added and the solution will go to the third stage. The solution procedure for this stage is exactly similar to the last one except that the unknowns at the intermediate nodes, as shown in Figure 4(c), should be calculated through simultaneous solution of the positive and negative characteristic equations and the energy equation across the node: Simplifying the above equations results in the following equation calculating the discharge at the associated node: where , .
The heads at upstream and downstream of the node can be then calculated by substituting the discharge into (2.18) and (2.19).
The numerical procedure presented for stage 3 can be employed for the rest of the simulation till the water column reaches the end of the pipe system.
3. Numerical Results
To evaluate the model results, the experimental data presented by Liou and Hunt [3] are first used. Their test apparatus consists of an upstream constant-head water supply reservoir and a PVC pipe with inside diameter = 2.29 cm, outside diameter = 2.67 cm, and length = 666 cm. The pipe extended 19.5 cm into the reservoir from a side wall near the tank bottom. The centre of the pipe inlet is set 3.43 cm above the reservoir bottom. A quarter-turn ball valve with an inside diameter of 2.44 cm is installed on the pipe at 43.82 cm from the inlet. The average slope for the first 335 cm of the pipe and the remaining portion of the pipe are 2.66° and 2.25° downward, respectively.
Steady flow tests showed that the average values for the entrance loss coefficient and the Darcy-Weisbach friction factor are 0.8 and 0.0242, respectively. Twenty filling tests were carried out, and the advance times were recorded in eight timing sections which are located at 57.6, 72.9, 118.6, 179.3, 301.6, 423.5, 545.5, and 661.7 cm from the inlet.
Figure 6 compares the water column velocity versus the water column length during the filling events. As can be seen, the water hammer model results accurately replicate the experimental data but the rigid column model proposed by Liou overestimates the velocities. Liou and Hunt attribute the discrepancy to the fact that the head required to accelerate the water in the reservoir was not considered in the calculation so that the water column received more deriving head than that in reality. However, the argument seems suspect since the particle velocity in the reservoir so the acceleration head does not seem to be significant except right at the inlet, so the associated acceleration head cannot absorb significant driving head. Instead, the discrepancy might be attributed to the velocity head which was neglected in Liou’s model. To test this idea, a different version of the water hammer model results which is obtained without considering the velocity head is also presented in Figure 6. As shown, the close similarity between the results signifies that the main culprit must be the velocity head.
Indeed, the understanding of the dynamic of the flow seems to be challenging because, in the different stages of the filling, the contribution of the driving and the resistive forces would be different. To examine how these forces come to action, consider Figure 7 in which the contribution of the acceleration head, velocity head, minor loss, and pipe head loss is demonstrated. This figure clearly shows that the quick changes only occur in the early stage of the filling, but as the water column increases, the pipe head loss mostly controls the pace of the motion. Putting aside the fully dynamic behaviour of the flow during the early stage of the filling, it is obvious that short the acceleration of the water column quickly reaches a low value and the quasisteady flow dominates the water column motion afterward.
Figure 6 also indicates that, even in the early stage, the rigid column motion dominates the flow in the rapid filling process. This can be physically explained through considering the energy relation during the rapid filling. As Karney [22] explained, as long as the compressibility index () which is the relation of the total change in internal energy to the total change in kinetic energy did not exceed 0.1, the rigid column theory dominates transient flow. Putting this argument in the context of the following research, one can conclude that the internal energy of the water column must not be changed significantly during the filling. This is indeed true because the water column has a greater tendency to occupy the empty pipe than to get compressed. In the absence of the significant compression, the compressibility index tends to zero because the internal energy of the water column experiences almost no changes during filling.
Finally, to evaluate if the proposed model is capable of capturing huge water hammer pressure spikes, the model was applied on a simple pipe system (see Figure 8) with a partially open valve (head loss coefficient, ) in its middle. The model results show that, once the water column reaches the valve, a high pressure spike would be created. The pressure spike travel back and forth between the reservoir and the valve until damped. Some numerical snapshots of the event captured by the model are shown in Figure 9. It should be noted that the damping is accomplished quickly even before the water column fills the immediate cell after the valve. This is because the high internal energy stored in the system would be quickly released though the partially open valve.
The water column velocity versus the length of the column is also presented in Figure 10 for two different cases, partially open valve and full open valve. It can be seen how the water column is retarded by the valve action. The figure also confirms that the water column velocity receives no fluctuations after the onset of the pressures spikes.
All these numerical evidences demonstrate the robustness of the proposed model in simulating the rapid filling in the pipe system, significantly even when the modelled system experiences strong water hammer pressures.
4. Summary and Conclusions
Based on water hammer theory, a numerical model is proposed for simulating rapid filling in the pipe system with undulating profile. Assuming the pipe remains full during the filling and neglecting the air and water interactions, the model is formulated through using the method of the characteristics along with an adaptable computational grid. The model is then verified by the experimental data presented in the literature. The model is also tested for the capability in capturing high water hammer pressures through a hypothetical example.
According to the results obtained from the model, the following conclusions are made. First, the water column moves in the system as a rigid column during the rapid filling provided that the water column is not locally disturbed along its length. Second, rapid filling can have dramatic effects at its early phase when the water column experiences high acceleration. As the water column lengthens, the acceleration reduces and the pipe head losses control the pace of the filling front. Under such a circumstance, a quasisteady model can also be utilized to capture the bulk motion of water. However, if the water column be disturbed along its length, neither the rigid column nor quasisteady models are capable of capturing the main feature of the flow. Third, the interaction of the moving column with the inline-partially open valves results in the onset of overpressures whose magnitude depends on the water column velocity at the moment of contact with the valve and also the initial rate of the opening of the valve. Though the overpressures can be severe, they are not persistent. In fact, the overpressures are quickly damped because the valve efficiently releases the internal energy accumulated in the system.
Finally, the numerical evidences presented in this paper demonstrate that the proposed model is robust enough to capture the water hammer pressures (if any) while simulating bulk motion of the water column. Capturing such water hammer pressures, as is attempted here, is indeed a vital component of any serious study of column separation in systems having an undulating profile. The interesting conclusion of the work is that applied modelling of complex physical process is possible but great care must be used to match the numerical aspects to the physical ones and that some subtlety is needed to ensure a good agreement between predicted and measured responses.