Transient flow analysis using the method of characteristics MOC with five-point interpolation scheme

An original 2nd-order Method of Characteristics (MOC) which it works with a five-point interpolating scheme valid for solving the hyperbolic and quasilinear partial differential equations that describe the transient flow phenomenon in pipelines is shown. The results obtained by both MOC 2nd-order and exact solution (MOC 1st-order) are compared. It is shown that MOC 2nd-order allows obtain near-to-exact results within a wide range of Courant numbers.


Introduction
Water hammer is a hydraulic phenomenon that it appears in pipe networks when flow is affected by some disturbance that it alters the flow velocity, such as the valve opening or closure, the relief devices failure, the pump start-up or shutdown, etc. Pressure fluctuations can cause serious damage in piping and hydraulic devices, which it can result in accidents and other consequences such as costs increase associated with water supply and water quality. For this reason, in transient flow analysis it is very important to know about the pressure waves and about the prevention or minimization of their undesirable effects. Wave analysis is a vital task that it must be included in the water supply systems design to ensure their safety and reliability in emergency conditions, and also to detect obstructions, defects, leaks (losses) and pipe wall damages (Duan et al., 2014). Ignoring the water hammer effects may lead to erroneous pipe wall thickness' dimensioning during the design stage. Transient flow phenomenon is difficult to understand intuitively, although the current computers have allowed its automatic calculation -by means of numerical schemes-in systems with different complexity levels. The Method of Characteristics (MOC) has been widely used for transient analysis (Chaudhry and Hussaini, 1985;Chaudhry, 2013), where the equations governing the phenomenon are transformed into ordinary differential equations that can be solved through characteristic lines using finite-difference approximations. For numerical stability reasons, the Courant number (C n ) must always be less than or equal to 1.0; C n = aΔt/Δx , where a is the wave speed, Δt the time step and Δx = L P /N is the reach length with L P the pipe length and N the number of reaches. 1 storder methods yield satisfactory results when the friction factor is small and C n = 1.0. In systems transporting highly viscous liquids or in systems with small diameter pipes, the simulation may become unstable even with C n < 1.0. On the other hand, the interpolations to be applied whenever solved using the MOC with C n < 1.0 tend to soften wave fronts, and disturbances due to flow changes travel at a greater speed than the correct one (Chaudhry and Hussaini, 1985), Twyman, J. (2018). 24, 62-70 being necessary to apply alternative methodologies to overcome these disadvantages. Trikha (1975) recommends using different time steps for each pipe in order to obtain C n = 1.0. Despite its advantages, this methodology requires interpolating at the pipe boundary nodes, which it is an error source when modelling complex boundary conditions. Wiggert and Sundquist (1977) propose an interpolation algorithm in the space (x) axis, which it has been widely used in water hammer programs because of its efficiency and ease of programming. Holly and Preissman (1977) present a more complex formulation based on a 4 thorder two-point method that is significantly more accurate than other higher-order methods using interpolations based on the four or five-point calculation. Goldberg and Wylie (1983) present two interpolation schemes (one implicit and the other one explicit) that operate on the time (t) axis of the space-time grid, both capable of reducing the numerical dispersion when MOC is applied with C n < 1.0. Shimada and Okushima (1984) present two 2 nd -order methods based on the solution in series and Newton-Raphson that they deliver results quickly and efficiently by omitting trivial terms calculated within the truncation error. Swaffield and Maxwell-Standing (1986) cite the Everett method as appropriate to reduce the numerical attenuation associated with C n < 1.0. This method must be applied in tandem with Newton-Gregory method to interpolate in the nodes nearto-pipe boundaries. Lai (1989) presents three multimodaltype schemes that combine implicit schemes with explicit ones in order to attenuate the numerical dispersion. Hsu (1990, 1991) propose an improved version of the Holly-Preissman scheme, which however it presents a greater cost in terms of computational time. Sibetheros et al. (1991) investigate MOC applicability in tandem with spline-type interpolation polynomials for the water hammer's numerical analysis in a frictionless horizontal pipe, finding that spline curves application significantly improves overall accuracy respect to both the MOC (linear interpolation) and the explicit finite-difference scheme (2 nd -order). However, the method becomes complicated when it is required to represent the spline boundary conditions, and it presents problems when the transient modelling in short pipes is required. Ghidaoui and Karney (1994) investigate interpolation methods and the Holly-Preissmann scheme. They use the equivalent hyperbolic differential equation concept that it allows evaluate the numerical scheme consistency, providing a mathematical description of both dissipation and numerical dispersion independently of Courant. In addition, the analysis allows comparing alternative approaches, and it explains why high-order methods should be avoided. Karney and Ghidaoui (1997) propose a flexible discretization scheme where several interpolation schemes are combined with the method of wave-speed adjustments. This hybrid scheme includes interpolations in a secondary characteristic line that it minimizes the distance between the interpolated point and the original characteristic line.
The current article introduces an original MOC 2 nd -order scheme with linear interpolation in the spatial axis valid for the water hammer analysis in pipe networks. The results obtained by both MOC 2 nd -order and exact result with MOC 1 st -order with C n = 1.0, are compared. Equations governing the transient flow along with wave speed and MOC 1 st -order formulation are extensively discussed in the classic books by Wylie and Streeter (1978), Chaudhry (1979) and Watters (1984). These topics can also be studied in recent articles by Twyman (2016aTwyman ( , 2016bTwyman ( , 2016cTwyman ( , 2017. Finally, it is possible to study how to pose and solve boundary conditions within the MOC's context in Karney (1984) and Karney and McInnis (1992), so no further details will be given here.

Basic equations of the transient flow
When analyzing a volume control is possible to obtain a set of non-linear partial differential equations of hyperbolic type valid for describing the one-dimensional transient flow in pipes with circular cross-sectional area (Chaudhry, 1979;Wylie, 1984): where: (1) and (2) correspond to continuity and momentum (dynamics) equations, respectively. Besides, is the partial derivative, H is the hydraulic grade-line elevation, a is the wave speed, g is the gravity constant, A is the pipe crosssectional area, Q is the fluid flow, f is the friction factor (Darcy-Weisbach) and D is the inner pipe diameter. The subscripts x and t denote spatial and time dimension, . Transient flow analysis using the method of characteristics MOC with five-point interpolation scheme. 24, 62-70 respectively. Equations (1) and (2), in conjunction with the equations related with the boundary conditions of specific devices, describe the phenomenon of wave propagation for a water hammer event.

Wave speed
The more general equation to calculate the wave speed in fluids without air is (Twyman, 2016b;Wan and Mao, 2016;Wylie and Streeter, 1978): where K is the bulk modulus of the water, ρ is the water density, e is the pipe wall thickness, K is the volumetric compressibility modulus of water, E is the elasticity modulus of the pipe, c 1 is a factor related with the pipe support condition, generally equal to 1 -u 2 (u is the Poisson's ratio of the wall material), which corresponds to pipeline anchored against longitudinal movement (Wylie and Streeter, 1978).

Method of characteristics MOC
MOC is very used for solving the transient flow equations because it works with a constant wave speed and, unlike other methodologies based on 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) in accordance with the Courant condition. It is useful for modelling the wave propagation phenomena due to its facility for introducing the hydraulic behaviour of different devices and boundary conditions (valves, pumps, reservoirs, etc.). According to Karney (1984) and Karney and McInnis (1992) the MOC proceeds by combining the dynamic and continuity equations together with unknown multiplier. Suitable chosen values of this multiplier allows the partial differential terms to be combined together and replaced by total differentials. When using the simplified governing equations, the result of this process is: The two equations associated with the positive value are usually termed the C + equations and the remaining two relations associated with the negative value are called the Cequations ( Figure 1). The head and flow values are known at time t and it is desired to know these values at time t + Δt. A typical such point is P, with unknowns H P and Q P . The known values at time t can correspond to initial state values or, instead, they can correspond to computed values at previous Δt. From Figure 1, Q and H are known in the nodes (i -1) and (i + 1), and their new values must be calculated at node P. Due to this the characteristic lines are projected from P to the x-axis so as to be able to intercept it at points L and R. Because of x P and t P are specified by the analyst, L and R coordinates can be written as: In order to calculate P the values in L and R must be known. However, only the grid point values in (i -1) and (i + 1) are known. The values in L and R can be computed using linear (or higher order) interpolation schemes from the known conditions in (i -1) and (i + 1). It is known that interpolations cause numerical dispersion and attenuation, which it has forced to apply more accurate high-order interpolation schemes whenever the pipe has C n < 1.0 and when it is desired to leave unchanged the initial conditions of the problem (lengths, wave speed, etc.), or when there is a lack of a more stable numerical scheme than MOC. In MOC context, some analytic expressions can be obtained for the state variables at L and R nodes using numerical schemes with different interpolation orders (usually one or two). If U represents some state variable (H or Q), then the following expressions can be deduced according to Newton-Gregory interpolation scheme (i = 2,..., N) . 24, 62-70 (Twyman, 2004): (8) and (9) are general since when C n →0, U L →U R → U i . When C n = 1.0 the interpolation does not exist and U L becomes U i-1 and U R becomes U i+1 . On the other hand, when C n = 2.0, U L becomes U i-2 and U R becomes U i+2 . For this reason, the application of 2 nd -order arrangements is sufficient to interpolate when the Courant number varies between 0 and 2.0. In this case the values for the variables U 0 and U N+2 , both located beyond the nodes 1 and N + 1 respectively (Figure 2), they can be calculated by an extrapolation procedure defined by the following expressions: Figure 2: Pipe nodes (adapted from Chaudhry and Hussaini, 1985) where U 1 and U N+1 are the values corresponding to upstream and downstream pipe sections, respectively. Once known U L and U R is possible to calculate H P and Q P from the following formulas when U is equal to Q and H according to the convenience: where c = (gA)/a and R f = f/2DA. Equations (11) and (12) can be rewritten as: with: It should be noted that (16) is valid along the positive characteristic line C + and (17) is valid along the negative characteristic line C -. By simultaneously solving (14) and (15) it is possible to obtain H P and Q P for all interior points P in the characteristic mesh: once Q P is known H P , value can be calculated from (14) or (15). For boundary sections an additional formula is required, which it must be solved together with the characteristic equation (positive or negative) depending on the position of the first or last reach, respectively. The calculation starts with an initial condition, usually steadystate flow.

Artificial viscosity
Artificial viscosity (AV) is a fictitious mathematical term that it is introduced in the equations when working with finite differences, and it generates a similar effect to the real viscosity. The AV is a numerical technique that generates a fictitious attenuation over the spurious instabilities that frequently appear when fast transient flow is solved using 2 nd -order numerical schemes. These instabilities generally appear in peaks form without physical meaning, and they may lead to destabilize the calculation due to equations nonlinearity (Brufau and García-Navarro, 2000). In other words, the phenomenon occurs when an approximate solution has an oscillatory and unrealistic behaviour respect to the analytical solution. In general, the precise cause of this effect is unknown, although it is known that they are generated when the first derivatives of the basic equations have more influence than the second derivatives or when the space-time grid is thick, with significant spacing between the nodes. Another cause could be the impossibility of the numerical schemes to capture all the discontinuities features, so there would be certain motion scales without numerical solution (Malekpour and Karney, 2016). The main difficulty of applying the AV lies in the dispersion amount determination needed to smooth (or . Transient flow analysis using the method of characteristics MOC with five-point interpolation scheme. 24, 62-70 eliminate) the numerical oscillations without sacrificing the method's accuracy level, and without unnecessarily affect the wave shape (Brufau and García-Navarro, 2000). The AV proposed in this paper recalculates H and Q every two time steps using the dissipation constant γ when the node i varies between pipe sections 2 and N, being the formula as follows (Abbott and Basco, 1989): With U equal to Q or H as the case may be. In general, AV is easy to program and apply, although it has the disadvantage that its coefficient γ is difficult to estimate, being necessary to apply a trial/error procedure in order to determine its most suitable value. AV, independently of its shape, is an indispensable tool to suppress numerical oscillations (Hou et al., 2012) generated by high-order numerical schemes.

Results: example 1
Water hammer is analyzed in a simple pipe network ( Figure  3) using MOC 1 st -order and 2 nd -order. System includes a constant head reservoir (upstream end, node 1) with H 0 = 100 m, a steel pipe and a valve (downstream end, node 2), where H 0 = 98.1 m, being the pipe head loss due to friction Δh f equal to 1.9 m. The pipe data are: length L P = 4800 m, wave speed a = 1200 m/s, pipe cross−sectional area A = 3.14 m 2 , initial flow Q 0 = 2.632 m 3 /s. Pipe network was discretized with N = 10. The friction factor (Darcy) is f = 0.022. Transient flow is generated by the valve closure in T c = 35 s (Figure 4), and it can be considered as slow because the ratio 2(L P / a) = 8 s is much smaller than T c . For all effects the exact result is given by MOC 1 st -order with C n = 1.0. The maximum simulation time is T max = 60 s. Figures 5 and 6 show the envelopes for the maximum and minimum pressure heads obtained by MOC 1 st -order and 2 nd -order when the pipe is discretized in order to obtain the following Courant numbers: 0.2, 0.4, 0.6, 0.8 and 1.0, which is achieved with the time steps Δt = 0.08, 0.16, 0.24, 0.32 and 0.40, respectively. For example, the value C n = 0.2 results from the expression: aΔtN/L P = 1200 · 0.08 · 10/4800. Figure 5 shows that MOC 1 st -order generates attenuations in the maximum pressure head that tend to decrease as Courant number grows toward the critical value C n = 1.0. Figure 6 shows how MOC 2 nd -order also generates attenuations, but on a smaller scale compared to 1 st -order MOC as C n moves away from 1.0. When C n < 1.0, Figures 7 and 8 show that MOC 2 nd -order is more accurate than MOC 1 st -order. Regarding CPU . 24, 62-70 execution times, the speed of both schemes is almost similar, with no significant differences (Figure 9). The examples were carried out in a standard PC with processing speed of 1.48 GHz.

Results: example 2
System analyzed (Amara et al., 2013) is shown in Figure  10, which consists of a steel pipe of length L P = 10000 m, wave speed a = 1000 m/s and diameter D = 1 m. The steady-state flow is Q 0 = 2.0 m 3 /s, and the head in the constant level reservoir (node 1) is H 0 = 400 m. The friction factor is f = 0.01976 and the pipe head loss due to friction Δh f is equal to 65.3 m. Transient conditions are generated by instantaneous closure of the valve located at the pipe downstream end (node 2). In this case transient flow can be considered as fast because 2(L P / a) = 20 s is much larger than T c = 0 s. For all effects the exact result is given by MOC 1 st -order with C n = 1.0. The pipe was divided using N = 30. With the Courant number C n = 1.0, the time step is Δt = L P /(aN) = 10000 / (1000 · 30) = 0.333 s. The maximum simulation time is T max = 120 s. Figures 11 and 12 show the pressure versus head plot at the valve when the transient is simulated using MOC 1 st and 2 nd -order, where it is noticed in the second case the appearance of the typical spurious peaks related to the 2 nd -order schemes. These fictitious pressure peaks can be easily attenuated by applying the AV. Figure 13 shows the pressure head versus time plot when the MOC 2 nd -order is applied using the AV with γ = 0.01. . Transient flow analysis using the method of characteristics MOC with five-point interpolation scheme. 24, 62-70 Figure 13: Pressure versus head plot when MOC (1 st and 2 ndorder) is applied Depending on the system complexity, the detection of the most suitable γ value may mean performing a trial/error procedure. One interesting point about MOC 2 nd -order is that it can work with C n > 1.0 without losing accuracy or numerical stability. For example, when C n is greater than 1.0 (for example, 1.08 or 1.80), the results shown in Figure  14 are obtained. As expected, the MOC 1 st -order shows an unstable numerical behaviour, rapidly tending to a fictitious state of water column separation near pipe midpoint. On the other hand, MOC 2 nd -order was able to conclude the simulation without numerical problems, although it was necessary to adjust γ value up to 0.10 (when C n = 1.08) and up to 0.20 (when C n = 1.80), which allowed to obtain near exact result, especially when extreme pressures values must be shown. Table 1 shows MOC 1 st and 2 nd -order execution time when time step can be increased thanks to γ attenuation effect, which helps to eliminate the typical instability effect of the explicit schemes when they are applied with C n > 1.0. In Table 1 it can be verified that MOC 2 nd -order allows the increase of time step magnitude up to 82%, reducing the execution time by approximately 14%, being in this case faster than MOC 1 st -order.

Discussion and conclusions
For the numerical examples presented in this paper, it is evident that MOC 1 st -order works properly only when C n = 1.0, showing attenuations and numerical instability when C n < 1.0 or when C n > 1.0, respectively, which obviously limits its application field. This turns evident in case of attempting to solve the water hammer in complex pipe networks with great number of different pipes, where the accurate solution it can be impossible to reach without altering or modifying some initial condition. A solution way is to leave unchanged the initial conditions and try to apply MOC by solving the state variables using interpolation techniques. MOC 1 st -order with linear interpolation works fine only when C n is closer to 1.0, showing significant errors as C n moves away from 1.0. In contrast, MOC 2 nd -order with quadratic interpolation registers smaller magnitude errors when C n < 1.0, regardless of whether it is a slow or fast transient, although in the second case it must resort to a numeric filter (e.g. artificial viscosity AV) to soften the fictitious numerical oscillations. The AV is simple and easy of programming, although the optimal value determination may take some time due to the trial/ error procedure implementation, which may mean running the computational program as many times as necessary. On the other hand, MOC 2 nd -order formulation slightly affects the transient simulation performance with respect to the calculation speed and memory capacity use. In addition, the 2 nd -order interpolation approach has a good numerical behaviour when C n > 1.0, even though in certain cases the γ value must be adjusted. As general conclusion, MOC 2 ndorder is as fast and efficient as traditional MOC 1 st -order, and it is less sensitive to numerical effect (attenuation, instability) which is generated when transient flow is solved in systems where C n is different than 1.0, which makes a useful method for solving the transient flow in . 24, 62-70 pipe networks where its shape and configuration makes very difficult or impossible to fulfil with the Courant condition in all pipe sections.