Water hammer analysis using an implicit finite-difference method

The Implicit Finite-Difference Method (IFDM) for the solution of water hammer in pipe networks is presented. All the equations necessary to calculate the flow and pressure in each node of the network are shown in detail. Section-by-section coupling through the Karney equation allows obtain a system of equations for each pipe section which is easy to solve applying the Thomas’ algorithm. Also, there are presented the original expressions in the form of finite differences are presented for: (1) the frictional term RQP|Q| of the dynamics equation, (2) the transient friction factor proposed by Brunone-Vítkovský, and (3) the short pipe replacement elements which allow increase the time step. It is demonstrated that the proposed methodology allows modelling the transient flow with higher level of stability and numerical accuracy in comparison to the Method of Characteristics (MOC), especially when the Courant number (Cn) is less than 1. However, because of IFDM works with weighting coefficients (θ1 and θ2) which must adopt values generally close to 0.5 depending on the analyzed problem, the achievement of the best near-to-exact solution requires to analyze each case separately, being obligatory to apply a trial/error procedure that can make the analysis cumbersome and time consuming.


INTRODUCTION
There is little literature dealing with the water hammer simulation using finite difference methods of an implicit type, where the solution of the basic equations of transient flow requires solving a system of equations.This may be due to two reasons: the IFDM implementation has a certain level of complexity because it is necessary to solve a coupled set of equations (linear or nonlinear).The boundary conditions under the IFDM context are difficult to handle because it requires adding fictitious grid points located beyond the pipe boundaries, being also necessary to establish additional systems of equations that allow solve the flow and pressure in each boundary node of the network.In general, IFDM is unconditionally stable and is used for solving the transient flow in those cases where it is necessary to adopt larger time steps (Δt) without taking into account the limitations given by the Courant number (C n ).Nevertheless, IFDM is not exempted dispersion and numerical attenuation even when the transient flow is solved in a very simple pipe network [6].
Other authors [28] proposed solving the water hammer equations in a network system using the implicit central difference method to permit large time steps, where resulting nonlinear difference equations are organized in a sparse matrix and are solved using the Newton-Raphson procedure.The implicit centered method gives a solution which is marginally stable, where time increments may be used in the solution [21].Implicit computer programs are somewhat more difficult to write and debug, in that the simultaneous solution hides any source of error in programming.Hybrid methods (HM) based on IFDM are more efficient to model the unsteady flow in complex pipe networks, with greater flexibility and robustness than Method of Characteristics (MOC), particularly for large Courant numbers [16,17,22].This is because the HM avoids the interpolations required by the MOC, thus lessening dispersion problems.Samani and Khayatzadeh [18] propose a method in which the implicit finite difference was coupled with the Method of Characteristics to obtain the discretized equations for the disproportionate reaches.The developed numerical method has a good fit-level when it is compared with the exact solution available for many test examples.Some authors [9] reviewing the work of Wylie and Streeter [29], conclude that major advantage of the IFDM is that they are stable for large time steps.Computationally, however, implicit schemes increase both the execution time and the storage requirement and they need a matrix inversion solver because a large system of equations has to be solved.Other authors [1,14] recognize that HM based on IFDM is a sophisticated method which has proven to be efficient to model water hammer in complex pipe networks, with greater flexibility than MOC and other methods when the Courant number is greater than 1.IFDM is unconditionally stable; that is to say that its behavior is the same regardless of time step size or adopted Courant number.[15].Sepehran and Badri [19] propose a full implicit finite difference method that works using a non-symmetrical staggered grid.The comparison between experimental data and numerical results shows agreement with MOC approach.However, the scheme uses pseudo parameters to calculate the state variables in the pipe boundaries that tend to complicate the analysis unnecessarily.

GOVERNING EQUATIONS OF THE TRANSIENT FLOW
When analyzing a volume control it is possible to obtain a set of non-linear partial differential equations of hyperbolic type valid for describing the one-dimensional (1-D) transient flow in pipes with circular cross-section [7]: Where: equations ( 1) and (2) correspond to the continuity and momentum (dynamics), respectively.Besides, ∂ = partial derivative, H = piezometric head, a = wave speed, c = (gA/a), g = gravity constant, A = pipe cross−section, Q = fluid flow and R = f/2DA with f = friction factor (Darcy−Weisbach) and D = pipe diameter.The subscripts x and t denote space and time dimensions, respectively.Equations ( 1) and (2), in conjunction with the equations related the boundary conditions of specific devices, describe the phenomenon of wave propagation for a water hammer event.

WAVE SPEED
For water (without the presence of free air or gas) the more general equation to calculate the water hammer wave speed magnitude in one-dimensional (1-D) flows is [23,24,26]: With K = water bulk modulus; ρ = liquid density; E = pipe elasticity modulus (Young); e = pipe wall thickness and ψ = factor related with the pipe supporting condition.

METHOD OF THE CHARACTERISTICS (MOC)
The Method of Characteristics (MOC) is an Eulerian numerical scheme [27] very used for solving the equations which govern the transient flow.
Because it works with "a" constant and, unlike other methodologies based on a finite difference or finite element, it can easily model wave fronts generated by very fast transient flows.MOC works converting the computational space (x) -time (t) grid (or rectangular mesh, see Figure 1) in accordance with the Courant condition.It is useful for modelling the wave propagation phenomena in water distribution systems due to its facility for introducing the hydraulic behavior of different devices and boundary conditions such as valves, pumps, reservoirs, etc. [16].Some main advantages of MOC are its ease of use, speed and explicit nature, which allows calculate the variables Q and H directly from previously known values [5,30].The main disadvantage of the MOC is that it must to fulfil with the Courant stability criterion that can limit the time step magnitude (Δt) common for the entire network.The MOC stability criterion states that [26]: With Δx = reach length, N = number of reaches and L = pipe length.In order to get C n = 1.0, some initial pipe properties can be modified (length and/ or wave speed).Another way is to keep the initial conditions and apply numerical interpolations with risk a of generating errors (numerical dissipation and dispersion) in the solution [10].In general, MOC gives exact numerical results when C n = 1.0; otherwise, it generates erroneous results in the way of attenuations when C n < 1.0 or numerical instability when C n > 1.0 [10].

IMPLICIT FINITE DIFERENCE METHOD
According to Figure 2, the state variables at time t have been calculated, being necessary to calculate new values at time t = t +Δt.To achieve this purpose, the partial derivatives in equations ( 1) and ( 2) may be approximated by finite differences, as follows [22]: On the other hand: Taking into account the equations ( 5) to ( 9), the equations ( 1) and ( 2) can be represented as follows [22]: Where the coefficients are the following [22]: Where θ 1 and θ 2 are weighting coefficients, ε = linearization constant and f av corresponds to an average value

CALCULATION OF THE BOUNDARY CONDITIONS
Karney [12] proposes a comprehensive and systematic approach to model different boundary conditions (with constant or variable consumption, reservoirs, etc.).and which it allows linking the hydraulic behavior of all piping, consumptions or tanks connected to each network node by: Where C c and B c are known constants and Q ext is the nodal consumption or flow rate that may be constant or a function of time.The importance of equation ( 20) is that it enables to separate or decouple the pipelines of complex networks in each node, restoring the flow continuity and the piezometric head in the node (when there is not storage nor singular losses).

COUPLING EQUATIONS
Equations ( 10) to ( 19) and ( 20) can be used to construct a diagonal band system of dimensions 2X(N+1) -see Figure 3, which can be efficiently converted into a tri-diagonal system and then solved using the Thomas algorithm [17].Figure 3 shows the system of equations for a pipe divided into N reaches, where the first and last rows correspond to H P t+Δt value (equation 20) of each boundary node.The remaining rows are generated from application of the equations ( 10) and (11) in each pipe reach.All coefficients have a superscript that identifies the reach number (for example, in c 1 2 the superscript 2 indicates that c 1 is calculated using the reach 2 data).This indication is relevant because the state variables may vary from one reach to another.For this reason, the matrix will not be always simmetrical although its structure has a form of bands, as shown in Figure 3.
Each system of equations will have a different size according to the number of reaches assigned to each pipe.Once the network is disconnected, it generates a system of linear equations for each pipe of the network.The calculation of the state variables is realized in each pipe independently from the rest.The algorithm works as follows for every pipe of the system: For each internal node i: • Compute H i t+Δt and Q i t+Δt for i = 2 to N using equations 10 to 19 arranged according to Figure 3. • Compute H P t+Δt for each network node using equation 20.
• Solve the system of equations (Figure 3) and reassign doing b) Back to a) until complete the simulation time.

TERM RQ|Q|
Friction always acts against the direction of flow (reversal) and it must change its sign accordingly.The effect of friction is expressed in terms of Q 2 in the basic differential equations, and it does not fulfill this requirement.Such difficulty can be eliminated by writing Q 2 as Q|Q| which has the same magnitude as just mentioned and it changes its sign automatically as flow reverses [20].A number of approaches have been proposed for the friction term of the dynamics equation, standing out the following approximation due to its general form [13]: Where dx = variation of x.The approximation given by equation ( 21) can be represented in the dynamics equation ( 10) by the IFDM as follows [22]: Equation ( 21) provides an excellent way for analyzing the effect of the friction factor on the transient simulation [13].It is suitable to adopt ε = 1.0 in equation (21) to obtain greater numerical stability, because when f value is high, instabilities can be produced whenever C n = 1.0 [11].Hence the approximation of RQ|Q| using ε = 1.0 ensures a more stable value which is independent of f [11].Wylie [28] has recommended the application of equation ( 21) with ε = 1.0 to obtain an improved alternative, with only a minor penalty in computational effort.
Nevertheless, some authors [13] indicate that values of ε near 0.85 are almost optimal for most applications.The main advantage of equation ( 21) is which weighting term ε influences the friction factor without requiring the discretization (Δx, Δt or a) to be changed [13].Hence, equation ( 21) provides an excellent way of assessing the sensitivity of a transient simulation to friction values.For example, if two values of ε produce significantly different results, the MOC (or IFDM) grid is too coarse, and smaller time step is required [13].Finally, is convenient to state that the friction term analyzed in this section comes from stationary state that introduces progressive distortion in the modeling of the phenomenon towards the middle and end cycles of the transient wave, both in amplitude and phase.

TRANSIENT FRICTION FACTOR
Another related problem with dynamics equation referrers to the determination of f, since the values of the friction factor of each pipe could affect the magnitude of the water hammer pressures.In this sense some authors have stated a number of ways to determine f, going from theoretical approaches up to others of practical nature.Among theoretical methods highlights the following equation for modelling the transient friction factor [4]: Where f t = transient friction factor; x = axial coordinate and k = decay coefficient which is dependent on maximum piezometric heads (respect to the steady−state value).It must be noticed that whenever k = 0 in equation ( 26), the friction factor corresponding to steady state (f st ) is obtained.In equation ( 26) three parameters have to be estimated: the characteristic roughness size η (Colebrook-White equation for f st ), a and k [4].Is important to note that Colebrook-White's formula is limited to the range of the Reynolds' number (R e ) in the transition regime (where both, R e and the pipe's relative roughness have influence on the friction factor).In flow with R e > 2105 and rough pipes is suitable to apply the second formula of Karmann-Prandtl.The estimation of f s must be carried out taking into account both the pressure and discharge data in a steady-state condition, the periodicity of pressure waves and the damping of pressure peaks [4].Brunone's approximation has significantly improved modeling of unsteady friction, and the decay coefficient k can be analytically deducted from Vardy and Brown's shear decay coefficient C * , obtaining with this procedure an unsteady friction model using a variable k based on instantaneous local conditions [25].In general, equation ( 26) can be represented in the dynamics equation ( 2) by the IFDM as follows [22]: Where:

MODELING OF SHORT PIPES
Water distribution systems (WDS) contain many pipes with a great variety of lengths and wave speeds.MOC requires that each pipe has at least one sub-section (N = 1) as a precondition to discretize the pipe network from a common Δt to all pipes.However, the adopted discretization can affect stability and/or convergence of MOC by not ensuring that all sections of the network have C n = 1.0.This can worsen when there are relatively "short" pipes limiting the magnitude of Δt and N, so that the most pipes are discretized with C n < 1.0, forcing the application of interpolation processes that can degrade the solution or to make adjustments to the wave speed which means modify the initial conditions of the problem.Alternatively, some authors [12][13] propose to replace the shortest pipe of the network for a mathematical representation called "pipe replacement element" PRE (Figure 4) for obtaining a greater Δt without altering the optimum conditions of stability and convergence of MOC.PRE can be represented as lumped inertia element (PRE-LIE) or as finite difference (PRE-FD) approximation [22].This last option works with the following equations taking into account the equation 20 [22] -see Figure 4: With: Where E i (i = 1, …, 6) are explicit coefficients which are dependant of d i (i = 1, …, 4) -see equations ( 12) to ( 15), as follows: Known Q ext1 and Q ext2 is possible to obtain the piezometric heads H P1 t+Δt and H P2 t+Δt for Δt applying equations ( 33) and (34) in ( 31) and (32).

EXAMPLE OF APPLICATION 1
The single pipe apparatus [3] used for investigating water hammer waveforms comprises a metal (copper) pipeline of length 37.2 m, 22 mm internal diameter and 1.6 mm wall thickness that is upward sloping (Figure 5).The transient event is generated by a rapid closure of the downstream end valve.The apparatus is installed in Robin hydraulic laboratory of the Department of Civil and Environmental Engineering at the University of Adelaide [3].The initial flow velocity is V 0 = 0.2 (m/s), static head in the tank 2, H T = 32 (m), valve closure time T c = 0.009 (s), and water hammer wave speed a = 1319 (m/s).The number of reaches for each computational run is N = 32, and the time step (Δt) = 0.00088135 (s).The average error in Figures 6 and 7 between the results obtained by MOC and IFDM is just +0.4% and +0.8%, respectively.In this case the computation runtimes were: MOC = 3.3 (s) and IFDM = 15.9 (s).The software was carried out a PC with CPU N280 @ 1.66 GHz.

EXAMPLE OF APPLICATION 2
Figure 14 shows the single pipeline which will be used to verify the application of the PRE [30].The pipe network in Figure 14 consists in a single pipe leading from a reservoir to a valve.(s) and its closure curve is defined by [30]: For analysis purposes, the initial pipe length was divided into three parts, and Table 3 shows the adopted discretization in this case when Δt = 0.010 (s).Because pipe 2 is too short relative to pipes 1 and 3, is convenient to replace it by a PRE in order to avoid solving the problem using a too small time step.

CONCLUSIONS
IFDM constitutes a useful tool for simulating the water hammer in pipe networks.Nevertheless, it spends more computational memory and takes longer to complete the execution simulation time.
A main disadvantage of IFDM is that it must be applied using different combinations of θ 1 and θ 2 to achieve the best result.In example 1, when C n = 1 the IFDM gets a near-to-exact solution when θ 1 = θ 2 = 0.501 (Figures 6 and 7).When C n = 0.5, IFDM obtains a smaller numerical attenuation with respect to the MOC when θ 1 = θ 2 = 0.550 (Figures 8 to 11).IFDM obtains a result with transient f similar to that reported in the literature [3] when θ 1 = θ 2 = 0.515 (Figures 12 and 13).In this case the transient f was applied in tandem with the formulation for variable f [8].It is found that transient friction factor generates an attenuation of pressure, being anyway a less conservative but closer to reality solution according to the results reported by some authors [3,4,25].In other words, the simulation results show that the steady friction model underestimates the damping observed in the experimental results.Additionally, the steady friction model does not predict the evolution of the shape of the pressure oscillations, demonstrating steady friction's inability to model frequency-dependent attenuation [3].MOC use within a rectangular space-time grid can create instabilities that can increase artificially the numerical attenuation of the transient friction [25].For this reason, the results obtained by the transient f should be considered with caution.In the case of the frictional term RQ|Q|, the IFDM obtains the best solution when θ 1 = θ 2 = 0.510 (tables 1 and 2).Analyzing the results obtained using equation (21) with ε = 0.0 and ε = 1.0 it verifies that the friction factor influence on the results is irrelevant because the differences in each and other case are very small.For that reason, in this case the adopted grid (Δx, Δt) is adequate and a smaller time step is not required.In example 2, the use of the two-node boundary condition called PRE allows reaching a more efficient and near-to-exact solution when θ 1 = θ 2 = 0.500 (Figure 15).Some authors [30] have recognized that the PRE can be used to know the transient response of a system in which excessively short tubes slow the program execution because of the time step becomes smaller and the simulation requires a greater amount of computational memory.The general conclusion is that the IFDM must adopt different values for θ 1 and θ 2 depending on the case, which forces to analyze each case separately applying a trial and error procedure that can make the analysis slow and cumbersome.The need to adopt different values of θ 1 and θ 2 to obtain the best solution in each analyzed case allows the conclusions of Arfaie and Anderson [2], who affirm that better results are obtained when the IFDM is applied using fixed values for the weighting factors equal to 0.500 and 0.515.Although these last values allow control the spurious numerical peaks that appear in the IFDM when very fast transients are simulated, this article shows that a more exhaustive analysis is required in each case.

Figure 3 .
Figure 3. Band diagonal system for a pipe with N reaches.

Figures 8 , 9 ,Figure 7 .
Figures 8, 9, 10 and 11 compare the pressure vs. time curve at the pipe midpoint and in the valve when a transient flow is analyzed using both MOC and IFDM (θ 1 = θ 2 = 0.550) with C n = 0.5.In all cases, the results are compared with exact solution (MOC with C n = 1.0).Comparison between f steady and f transientFigures 12 and 13 compare the pressure vs. time curves when a transient flow is analyzed using

Figure 15 .
Figure 15.Results obtained by MOC (without PRE) and IFDM (with PRE) at the valve.

Table 3 .
Discretization for the single pipeline.

table 4
) when Δt = 0.045 (s).Replacement of pipe 2 by a PRE allowed reduce the size of N and increase the magnitude of Δt, allowing reduce the execution time of the program from 26.4 (s) to 5.8 (s).

Table 4 .
Re-discretization for the system.