An Overview of the Numerical Approaches to Water Hammer Modelling: The Ongoing Quest for Practical and Accurate Numerical Approaches

: Here, recent developments in the key numerical approaches to water hammer modelling are summarized and critiqued. This paper summarizes one-dimensional modelling using the ﬁnite difference method (FDM), the method of characteristics (MOC), and especially the more recent ﬁnite volume method (FVM). The discussion is brieﬂy extended to two-dimensional modelling, as well as to computational ﬂuid dynamics (CFD) approaches. Finite volume methods are of particular note, since they approximate the governing partial differential equations (PDEs) in a volume integral form, thus intrinsically conserving mass and momentum ﬂuxes. Accuracy in transient modelling is particularly important in certain (typically more nuanced) applications, including fault (leakage and blockage) detection. The FVM, ﬁrst advanced using Godunov’s scheme, is preferred in cases where wave celerity evolves over time (e.g., due to the release of air) or due to spatial changes (e.g., due to changes in wall thickness). Both numerical and experimental studies demonstrate that the ﬁrst-order Godunov’s scheme compares favourably with the MOC in terms of accuracy and computational speed; with further advances in the FVM schemes, it progressively achieves faster and more accurate codes. The current range of numerical methods is discussed and illustrated, including highlighting both their limitations and their advantages.


Introduction
Increases in computational power create an opportunity for analysts to better understand many complex phenomena, including water hammer. In essence, water hammer events involve hydraulic shock waves generated in a pipeline by rapid changes in one or more boundary values, such as valve closure. Such actions can have a devastating effect on system operation if improperly managed. As per Toro [1], shock waves of strong character are a very rapid transition of a physical quantity (i.e., pressure, density, or temperature). Water hammer problems with a rapid boundary change characterized, for example, by a valve closure time that is less than the wave return time over the entire system length, produce such rapid changes of pressure and velocity that the motion is typically characterized as a mathematical discontinuity. However, even with significant increases in computational speed, complex and multidimensional fluid problems still require a daunting numerical budget, whether using computational fluid dynamics (CFD) or considering the detailed transient modelling of a complex pipe network. This paper reviews advances in the numerical modelling of transient pressurized flow, highlighting the use of the more recently developed finite volume method (FVM) for water hammer work.
What makes the finite volume (FV) approach special is that it evaluates the conservative forms of the governing partial differential equations (PDEs) in algebraic form [1,2]. Indeed, the FVM is developed by considering a small volume surrounding each nodal point in a geometric mesh, and then relating the volume integral to the surface integrals using the divergence theorem. The resulting surface integrals act as fluxes at the boundaries of each finite volume, and the whole formulations enable the tracking of the evolution of a system's state.
In addition to a general overview of the traditional method of characteristics (MOC) and finite difference method (FDM) computational approaches, this paper highlights some newer developments in the numerical approach to water hammer analysis. In particular, the paper summarizes the inherent limitations of traditional methods-limitations that naturally lead to the introduction of the various FV approaches. The paper concludes by suggesting several directions for future research.

Governing Equations
Numerical estimates of water hammer pressure and flow are obtained by approximating the solutions of the governing continuity and momentum equations. To that end, several numerical methods are available.

Full-Form 1D Equations
The 1D governing equations for water hammer can be derived from considerations of conservation of mass and momentum, and are often written as follows [2,3]: where P is the pressure, V is the velocity, ρ is the fluid density, g is the gravitational acceleration, θ is the pipe slope, D is the pipe diameter, a is the celerity, f is the Darcy-Weisbach friction factor, x is the spatial coordinate along the pipeline, and t is the time.
The one-dimensional flow equations can also be written as:

∂U ∂t
where U is the vector of flow variables, F is the flux vector, and S is the source term vector. By comparison, in the one-dimensional Euler equation used in gas dynamics, the terms U and F are represented as: where total energy E is ρ 1 2 V 2 + e , the first term specifies the kinetic energy, and the second term specifies the internal energy of the gas. Note that the governing equation for water hammer varies from that of gas dynamics (Euler equation), primarily in terms of the energy conservation equation. In water hammer, the energy conservation equation is dropped based on the (usually reasonable) assumption of a negligible change in water temperature over a transient event. Equation (2) can be further simplified as: where J is the vector of time dependent flow variables in the x direction. Riemann Invariant: The Riemann invariant is a key conceptualization in various formulations, and is essentially a mathematical transformation of the governing conservation forms. Riemann invariants assist in the solution of the governing equations along the characteristic curves. Considering only the hyperbolic portion of Equation (4)-that is, assuming the source term S = 0-the governing equation can be expressed as K −1 ∂U ∂t + K −1 JKK −1 ∂U ∂x = 0, and given that dW = K −1 dU produces: where Λ is a diagonal that contains the eigenvalues of matrix J, hence for the pth term: where λ (p) is the pth eigenvalue.
Hence, the total derivative is given as DW (p) Dt = 0 along dx dt = λ (p) . Thus, this equation effectively defines the invariance of function W(p) along the characteristic line x = λ (p) t + x 0 , where W(p) is the pth Riemann invariant. Now, for the more general case involving m − 1 waves, the Riemann invariants can be expressed as: where K (p) 1 is the pth eigenvector of Λ. Equation (7) is known as a generalized Riemann invariant. A schematic of the generalized Riemann problem is graphically conceptualized in Figure 1. Considering only the hyperbolic portion of Equation (4)-that is, assuming the source term = 0-the governing equation can be expressed as K + K JKK = 0, and given that dW = K dU produces: where Λ is a diagonal that contains the eigenvalues of matrix , hence for the pth term: where λ ( ) is the pth eigenvalue.
Hence, the total derivative is given as ( ) = 0 along = λ ( ) Thus, this equation effectively defines the invariance of function W(p) along the characteristic line x = λ ( ) t + x , where W(p) is the pth Riemann invariant. Now, for the more general case involving m − 1 waves, the Riemann invariants can be expressed as: where ( ) is the pth eigenvector of Λ. Equation (7) is known as a generalized Riemann invariant. A schematic of the generalized Riemann problem is graphically conceptualized in Figure 1.

Simplified Equations
An early fundamental form for an analytical solution of the simplified water hammer equations was presented by Joukowsky [4] as follows: is the celerity or acoustic wave speed, is the pressure head = ( − ), is the elevation of the pipe centreline from a given datum, is the piezometric head, is the fluid density, is the cross-sectional average velocity, is the local longitudinal velocity, is the cross-sectional area of the pipe, and is the gravitational acceleration.

Simplified Equations
An early fundamental form for an analytical solution of the simplified water hammer equations was presented by Joukowsky [4] as follows: where a is the celerity or acoustic wave speed, P is the pressure head = ρg(H − z), Z is the elevation of the pipe centreline from a given datum, H is the piezometric head, ρ is the fluid density, V is the cross-sectional average velocity, u is the local longitudinal velocity, A is the cross-sectional area of the pipe, and g is the gravitational acceleration. The ± sign typically represents whether the water hammer wave is moving upstream or downstream. Allievi [5] developed the theory in the form of a differential equation relating the pressure gradient to the total derivative − ∂P ∂x = ρ g dV dt . The underlying assumptions in Equations (1), (2), and (4) include slightly deformable pipes and slightly compressible fluid. If we were to suppose that the fluid is incompressible and that the pipe walls are rigid, the transient flow then becomes a so-called rigid column motion-a motion in which the fluid velocity and pressure changes propagate instantaneously throughout the pipe system. The rigid column theory is generally applicable to relatively slow transients when the action of the boundary conditions gives rise to inertial effects but is too slow to evoke notable compressibility effects.
Solving Equations (1) and (4) numerically presents several challenges. Streeter and Lai [6] numerically formulated the governing Equations (1) and (4), initially accounting for their intrinsic nonlinearity. However, they also demonstrated that the nonlinear convective terms V ∂V ∂x and V ∂P ∂x are relatively small and, thus, justified their neglect. Moreover, the body force due to gravity is also often considered negligible. The simplified conservation PDEs, given by Chaudhry [3], Fox [7], and Wylie and Streeter [8] are often written as follows: a 2 g ∂V ∂x where H is the piezometric pressure head = P/γ, and τ w is the shear stress at the pipe wall.

Finite Volume Method
The equations presented so far are considered the classical 1D water hammer equations, and full derivations are found in the works of Chaudhry [3] and Wylie and Streeter [8]. Interestingly, alternate approaches usually utilize different presentations of the equations. Thus, for the FVM approach, the following form similar to Equation (2) is usually preferred [9,10] for its intrinsic conservation of mass and momentum:

∂U ∂t
, whereǍ is the average cross-sectional area of the pipe, M is the mass flowrate, µ is the mass of fluid per unit length of pipe, and M is the fluid's momentum M = µV = ρǍV . The wave celerity can be defined as a = d(Ǎ P) dµ . Variables µ and M are conservative; they can be converted to non-conservative pressure and density forms as required by using the following equation, which is obtained by integrating the wave celerity equation: where µ ref =Ǎρ ref , ρ ref is the reference fluid density, and a is the wave velocity. Equation (11) is classified as a hyperbolic system of conservation laws, and expresses its characteristic form assuming that U and F are continuous in space and time, and F is only a function of U [10]: Water 2021, 13, 1597 5 of 39 In this case, the Jacobian matrix of a flux vector is: where J is the Jacobian matrix of F in Equation (11). The solution methodology associated with Equation (11) is summarized here in the section on the finite volume approach.

Two-Dimensional Flows
If a more accurate estimate of shear stress is required, a 2D version of the flow can be used. The general governing equation of two-dimensional models after dropping the source term can take the following forms:

∂U ∂t
The constants J F and J G are the Jacobian matrices of F and G with respect to U [11]. The most common semi-two-dimensional, non-conservative administering conditions were given by Vardy and Hwang [12]. Although these conditions were produced using different methodologies, and are clearly structured differently, they can be described by the following pair of continuity and momentum equations in radial coordinates: where r is the radial distance, ν is the kinematic viscosity of the fluid, t is the time, and x is the distance along the pipe, H(x, t) is the piezometric head, V(x, r, t) is the radial velocity, and τ is the shear stress.
Owing to the turbulent nature of pipe flow, Reynolds averaging is performed in order to obtain the time-averaged solutions of Navier-Stokes equations. The time-averaging produces a new term known as Reynolds shear stress −ρu v , where u and v are velocity fluctuations in the x and r directions, respectively. Using the well-known Boussinesq approximation, the Reynolds shear stress can be expressed as −ρu v = µ T ∂V ∂r , where µ T is the eddy viscosity. Therefore, the total shear stress consists of two components-namely, viscous shear stress, and Reynolds shear stress. The total shear stress as a function of the velocity gradient is written as Neglecting nonlinear terms, Equations (17) and (18) can be rewritten as: Furthermore, Equation (20) can be simplified using dA = 2πrdr, which results in the following: Equations (19) and (21) [13]. One-equation models usually solve an equation for the turbulent kinetic energy, which is obtained during the averaging process. One-equation models include Prandtl's one-equation model, the Baldwin-Barth model, and the Spalart-Allmaras model. Although the Spalart-Allmaras model is a simple one-equation model based on the kinematic turbulent viscosity, it is not popular for water hammer analysis. More advanced models such as k − ε and k − ω models solve two equations-one for kinetic energy, and the other for the rate of dissipation of turbulent kinetic energy. An extended list of the turbulence models is given in Section 6.1.
Equations (19) and (21) are typically solved using one of a variety of available Reynolds-averaged Navier-Stokes (RANS)-equation-based turbulence models, such as two-layer or three-layer eddy viscosity models (e.g., [14] and [15]). Detailed implementation strategies for the eddy viscosity turbulence closure model for a system experiencing sudden valve closure are summarized in the works of Hanmaiahgari and Maji [16] and Silva-Araya and Chaudhry [17].

MOC Schemes
Although MOC schemes are well known by both researchers and practitioners, to facilitate comparison, it is valuable to quickly summarize the key MOC approaches here. In essence, the MOC approach converts the governing pair of PDEs into a set of four coupled ODEs [3,8]. This transformation preserves many of the attributes of the original equations, and lends itself to fast, accurate, and intuitive codes. However, one fundamental challenge is that the MOC specifies a fixed relationship between space and time steps (see Figure 2)-a requirement that can create numerical challenges in complex systems with many pipe segments, or in systems with a variable wave speed. Equations (19) and (21) are RANS (Reynolds-averaged Navier-Stokes) equations. Solution of RANS equations requires turbulence closure models. There are various turbulence closure models: zero-equation (mixing-length model), one-equation, two-equation, and seven-equation models. Zero-equation models are widely used: for example, the Cebeci-Smith model, Baldwin-Lomax model, Johnson-King model, and a roughness-dependent model [13]. One-equation models usually solve an equation for the turbulent kinetic energy, which is obtained during the averaging process. One-equation models include Prandtl's one-equation model, the Baldwin-Barth model, and the Spalart-Allmaras model. Although the Spalart-Allmaras model is a simple one-equation model based on the kinematic turbulent viscosity, it is not popular for water hammer analysis. More advanced models such as − and − models solve two equations-one for kinetic energy, and the other for the rate of dissipation of turbulent kinetic energy. An extended list of the turbulence models is given in Section 6.1.
Equations (19) and (21) are typically solved using one of a variety of available Reynolds-averaged Navier-Stokes (RANS)-equation-based turbulence models, such as twolayer or three-layer eddy viscosity models (e.g., [14] and [15]). Detailed implementation strategies for the eddy viscosity turbulence closure model for a system experiencing sudden valve closure are summarized in the works of Hanmaiahgari and Maji [16] and Silva-Araya and Chaudhry [17].

MOC Schemes
Although MOC schemes are well known by both researchers and practitioners, to facilitate comparison, it is valuable to quickly summarize the key MOC approaches here. In essence, the MOC approach converts the governing pair of PDEs into a set of four coupled ODEs [3,8]. This transformation preserves many of the attributes of the original equations, and lends itself to fast, accurate, and intuitive codes. However, one fundamental challenge is that the MOC specifies a fixed relationship between space and time steps (see Figure 2)-a requirement that can create numerical challenges in complex systems with many pipe segments, or in systems with a variable wave speed.  Resolving the integration of properties along the characteristic lines can be done either explicitly or implicitly. In the explicit method of characteristics, unknown values at time t + ∆t are obtained point-by-point by integrating characteristic lines from time t to time t + ∆t, in which the source term is computed using the known values at time t. In the implicit method of characteristics, unknown values at time t + ∆t are obtained simultaneously by integrating characteristic lines; however, for computational simplicity, the source (friction) term is computed by blending the known values at time t and the unknown values at time t + ∆t.
Gray [19] first used the MOC to compute pressure and velocity in a frictionless pipe system using fixed ∆t and ∆x intervals. Lister [20] investigated the fixed-grid method and the method of specified time intervals, and concluded that the fixed-grid method facilitates programming, leading directly to its popularity. Lai [21] included the frictional losses in the MOC. Streeter [22] developed the boundary conditions in the MOC for opening and closing the valves. Streeter and Lai [6] experimentally validated the MOC for different valve closure times. Perkins [23] proposed that the fixed-grid method is numerically stable only for Courant numbers less than unity. However, in order to efficiently resolve device behaviour with the fixed-grid scheme, a common time step for all pipes is essential, creating a challenge to satisfy the Courant condition C r = a∆t ∆x ≤ 1. To address this problem, space-line interpolations were suggested to approximate the head and flow at the foot of the characteristic line [20]. This approach allowed the use of different time steps for specific pipes, thus permitting faster codes with still manageable errors. Vardy [24] proposed reachout-space-line interpolation for applying the MOC to pipe networks. This method also requires interpolation at the boundary, causing generation of error and propagation of the error over time. To overcome the limitations of space-line interpolation, Wiggert and Sandquist [25] combined space-line interpolation and space reachout interpolation. In contrast, Goldberg and Wylie [26] developed a relatively more accurate time-line interpolation method, which uses the results from previous time step data to carry out the interpolation. A combination of all of the available interpolation techniques was presented by Lai [27] under the name of the multimode scheme. This method provides a vast range of choices for ∆t and ∆x, thus offering considerable flexibility. Sibetheros et al. [28] demonstrated that a spline strategy is often well suited to adapting MOC grids to more complex domains, and specifically to anticipating transient conditions at unknown points when the general transient behaviour is known. Figure 3 depicts various interpolation techniques.  Ghidaoui and Karney [30] analyzed a wave speed adjustment method; therefore, characteristic lines pass through the known-time-level grid points to unknown-time-level grid points. Karney and Ghidaoui [31] introduced a hybrid interpolation technique in which they considered a combination of numerous approaches for interpolation. The It is worth mentioning that interpolation is fundamentally non-physical, and inevitably introduces numerical errors in the form of dissipation and dispersion. Here, the term nonphysical is used, since it produces false or artificial diffusion, resulting from the numerical approximation of the convection term in the conservation equations [29]. Though dissipation and dispersion errors are normally presented, the first-order phase error is usually less apparent due to the numerical dissipation of the method. Moreover, the phase error is small; thus, a travelling wave proceeds (on average) at the correct speed a.
Moreover, interpolated methods have no intrinsic difficulty in accounting for the convection term in the momentum equations. In such an approach, the characteristic speeds change to u ± a. Thus, there is a slight asymmetry in the left and right running characteristics. For a user of MOC, there are two choices: using interpolation and accepting the numerical error but avoiding limitations to the physics and geometry; or, alternatively, to set the C r = 1 case, where the convection term must be dropped, and accepting some limitation or distortion on the system's geometry. The latter is a minor issue today, since the computing power permits the use of very short spatial increments ∆x.
Ghidaoui and Karney [30] analyzed a wave speed adjustment method; therefore, characteristic lines pass through the known-time-level grid points to unknown-time-level grid points. Karney and Ghidaoui [31] introduced a hybrid interpolation technique in which they considered a combination of numerous approaches for interpolation. The composite calculation [31] can be executed as a preprocessor step, and consequently utilizes memory effectively, executes rapidly, and provides an adaptable tool for discretization in pipe networks. Ghidaoui et al. [32] proposed an integrated-energy approach in order to comprehend and control the discretization errors associated with space-line and time-line interpolation methods. Shimada et al. [33][34][35] introduced an exact method to compute the numerical errors of both space-line and time-line interpolation methods. According to Ghidaoui et al., [36], interpolation distorts the physics, and it is particularly problematic in short pipes due to the requirement of involving several reaches, implying the need for tiny time steps and greater computational times.
Tian et al. [37] worked to moderate valve-induced water hammer effects on a parallel pipe system using the MOC method. Adopting an alternate startup process, they found that using a damping torque could effectively reduce the induced pressure vibration due to valve closure. Bergant et al. [38,39] highlighted several additional transient processes, including fluid-structure interaction (FSI), and showed how they could be included in the MOC scheme. They investigated various parameters affecting water hammer, and suggested terms for unsteady friction, cavitation, and viscoelasticity, as well as how to represent both leakages and blockages. The first of their papers deals with mathematical models to incorporate those effects, while the second presents case studies along with comparisons of the MOC to the rigid-pipe approach.
Afshar and Rohani [18] proposed an implicit MOC model to address the shortcomings of the widely used MOC model. They simulated two transient flow problems, and favourably compared the predicted results with an explicit MOC scheme.
Column separation generally occurs when the calculated pressure becomes less than the saturated vapor pressure-a remarkably common complication during certain transient flows. One of the most well-known and utilized methods for modelling cavitation is to use a discrete vapor cavity model (DVCM), as developed by Wylie and Streeter [8]. Bergant and Simpson [40] modified the standard DVCM model by incorporating shock waves and other interactions. This new model is known as the generalized interface vaporous cavitation model (GIVCM), and provides a simpler and more numerically stable approach to column separation. Bergant and Simpson [41] provided a detailed comparison between the experimental and numerical results. Sadafi et al. [42], using the GIVCM, investigated column separation caused by the rapid closure of a valve-an outcome that led to pipe failure.
When rapidly filling an empty pipe, the resulting momentum of the water column when arrested can lead to severe pressures-an outcome widely investigated. For example, Wang and Yang [43] combined an MOC approach with an implicit solution to the filling problem, and showed that a coupled explicit-implicit method can provide accurate and stable numerical results.
Jang et al. [44] developed a new model with a quasi-2D approach based on the MOC, in order to reduce the computational nodes by proposing 1D characteristic equations for pressure and discharge. They showed that this method can give equivalent results to those obtained with an explicit 2D solution.

Limitations of MOC
Historically, the most popular approaches for solving water hammer problems have been the MOC-based methods. The MOC initially gained popularity due to its explicit nature, accuracy, computational speed, and comparative simplicity. Such methods directly track the movement of hydraulic disturbances. In fact, the MOC preserves wave behaviour, which is useful for modelling high-frequency transient events. In the MOC, the convective terms are typically neglected, because the Riemann invariants in the water hammer problem are only weakly interdependent [9]. As long as the wave speed is constant and a V, dropping the convective term makes sense, and the MOC is well suited to describe transient events. In any case, the finite volume method tends to be as sensitive to the wave damping as the MOC if a first-order method is used. However, the MOC approaches have some notable disadvantages:

1.
The MOC approaches must be artificially modified if the Courant number is not set to unity [25,30]. Furthermore, mass and momentum conservation are not assured for non-uniform Courant numbers [45]. Interpolation is nonphysical, and tends to cause artificial damping of the solution [25,46].

2.
The MOC makes it difficult to include the convective acceleration term-a reality that can occasionally induce errors, particularly when the objective is to find a solution with C r = 1 [47].

3.
Conservation of momentum and mass cannot be ensured even if interpolation is not required.

4.
Overall, the MOC approach is best suited to systems with invariant celerity, since a condition of C r = 1 can be met.

5.
For C r = 1, explicit schemes are easy to program, and provide simpler solutions, whereas implicit schemes are more complex and computationally demanding to both code and execute [43]. The application of the explicit scheme is conditionally stable because of the source term, as discussed in Section 4.1, whereas the implicit scheme is unconditionally stable. Despite the popularity of explicit and implicit schemes, such approaches are often computationally inaccurate when used with larger spatial steps or conduits with variable cross-sectional areas [43].
Besides these concerns, complex device characteristics (e.g., for pumps and turbines), and even unsteady friction-not to mention certain air-water interactions-can be thought of as hints that simplified models themselves have limitations that can confound 1D analysis.

Need for Improved Numerical Methods
Thus, considerable challenges in water hammer calculations remain when using traditional MOC solvers. For example, when drawing water from a lake or reservoir, transient events can cause air to be released or to be entrained in the flow, causing the wave celerity to vary both temporally and spatially. Air release and reabsorption may require the modelling of heat transfer between air and water, which affects the dissipation of air pressure oscillations. In such cases, better numerical approaches may be both necessary and justified. Specific cases include multidimensional domains, cases where the wave speed is variable (e.g, associated with the gas release or readsorption with strong pressure changes), or occasionally in cases involving higher velocities within highly flexible pipes. Such applications tend to induce significant errors in traditional schemes, and thus justify a search for improved schemes.
In particular, the so-called inverse methods change the usual approach of predicting the response given by the system to one where the response is measured, and consistent system attributes are sought. According to Duan et al. [15], inverse methods use the entire transient response in the time domain for calibration, thereby utilizing both the leak-induced damping and reflection information. For that reason, inverse methods have become understandably more popular for finding anomalies in water networks and for calibrating systems. Inverse methods are typically more sensitive than typical design calculations to parameter values. Indeed, inverse problems are not well posed mathematically and, thus, demand more than an estimate of the maximum and minimum pressuresvalues that are often central to design questions. In fact, inverse methods perform well only when there is an accurate match between computed and measured time series. Unsteady pipe flows are sometimes characterized by temporal changes of wave speed and unsteady friction, influencing experimental outcomes and timings. In this context, the computational effort is still important, but achieving high accuracy, though of great value, has tended to be elusive to achieve in practice.

The Finite Volume Method (FVM)
The FVM is an example of a Godunov-type [48] algorithm, and is a technique for solving hyperbolic equations. Guinot [9,11,49] recommended this method because of its desirable attribute of conserving mass, momentum, and energy within a control volume for solving partial differential equations. This can solve either exact or approximate Riemann problems at each intercell boundary, primarily with the attribute of first-order accuracy in both space and time. The usually small control volume used in this context is often referred to as a cell ( Figure 4). Unlike finite difference methods, the FVM does not usually retain state values at grid points, but rather stores average values within a cell. Finite volume methods update these average values by accounting for the mass and momentum flux contributions at the cell boundaries within a single time step [50]. However, both the MOC and FVM are based on the Riemann invariant. If is derivable, then this characteristic form will be continuous with respect to space. If space is discretized in divisions, then for the th cell at time + 1 the discretization of the source term in Equation (24) is as follows: for −2 ≤ Guinot [49] considered only the hyperbolic term of the equation, and solved ( ) for [ , ] in between − and + cells, and later source term was incorporated into the equation, as presented in Equation (25). The Riemann [54] concept of wave propagation was applied to find the fluxes; / / and / / at the cell interfaces, for the first-order solution, become: In the MOC, non-conservative governing equations are integrated along the characteristic lines, whereas in the FVM, integration is done for an individual cell. In fact, there are many approaches based on the method chosen for time integration of the derivative with respect to space in the FVM [10]. Generally, integration in the FVM scheme is carried out in two steps: In the first step, non-conservative governing equations are integrated along characteristic lines, assuming a piecewise linear variation of the dependent variables, to find numerical fluxes at control volume (cell) boundaries. In the second step, conservative governing equations are integrated over time space within the control volume, and the computed values are represented at the cell centroid. Either only the interpolation is made on the boundaries, and a time integration is done separately using-for example-the Crank-Nicolson method, or a three-level implicit method; or else the second-order accuracy in time is obtained by adding the first step to compute values at time level t + ∆t 2 . There is an analogy between this method and the classic Lax-Wendroff method, or a two-step Lax-Wendroff method. Naturally, this approach is explicit, and is well suited for wave propagation problems. The first (implicit) approach is used by commercial CFD codes, as well as by thermal-hydraulic codes.

Finite Volume Method (FVM) Numerical Procedure
The finite volume method is broadly utilized in hyperbolic problems, including gas dynamics, shallow water flows, and many other contexts [1]. However, using this method for solving water hammer problems is fairly recent. This method has a clear advantage of better capturing shocks and discontinuities [51]. Mass conservation with the FVM approach is exact, but occurs in the MOC only when the solution is sufficiently accurate. However, in the case of two-phase flow, conservation never occurs with the MOC approach. For example, MOC methods usually deal with variable celerity problems by either interpolating parameters that are inaccurate and may cause mass conservation problems, or by forcing the celerity to be a constant using some procedure, such as that used in the discrete gas cavity model (DGCM). Such resolutions can cause serious numerical errors [52]. Although initially neglected, Guinot [9] applied the finite volume method to water hammer modelling by assuming a compressible pipe with its walls deforming elastically under pressure changes; that the conduit remains full; that friction losses and velocity head are negligible; that velocity is uniform over the cross-section; and that the reservoir level remains constant. He disregarded the advective terms, and built up a Riemann-type scheme, thus creating a first-order accurate FV scheme-an approach that was quite similar to the MOC with direct space-line interpolation. Guinot [9] expressed the governing equations for the water hammer in conservative form. Guinot's [9] approach has recently been implemented by Pal et al. [53].
The crucial takeaway from this is that this paper contains important concepts of the methodology and, thus, assists comprehension. In this method, Equation (13) is solved in two steps: The first step solves the homogeneous part of Equation (13), i.e., ∂U ∂t + J ∂U ∂x = 0. The next step solves the non-hyperbolic source term ∂U ∂t = S. Thus, the nonlinear hyperbolic problem is converted to an initial boundary value problem, and the partial differential equation will be of the form: ∂U ∂t having specified initial and boundary conditions as per problem requirements. Considering = i∆x, and the interface ∆x-and simplifying, the final form of Godunov's scheme can be obtained as: where U is the vector of flow variables, i is the number of a spatial cell, n is the number of temporal cells, ∆t is the time step, ∆x is the grid size, and F n+1/2 i−1/2 and F n+1/2 i+1/2 are the fluxes at the cell interfaces of the ith cell [49].
If U is derivable, then this characteristic form will be continuous with respect to space. If space is discretized in N divisions, then for the ith cell at time n + 1 the discretization of the source term in Equation (24) is as follows: is the solution after solution of the hyperbolic part.
According to Guinot [49], U n+1, * i obtained at the end of Godunov's scheme is taken as a starting point to produce the final solution U n+1 i . Thus, at the maximum permissible time step a , ∆t max, source , and for the uniform grid ∆t max = min ∆x a , ∆t max, source . Guinot [49] considered only the hyperbolic term of the equation, and solved F(U) for t n , t n+1 in between i − 1 2 and i + 1 2 cells, and later source term S was incorporated into the equation, as presented in Equation (25).
The Riemann [54] concept of wave propagation was applied to find the fluxes; F n+1/2 i−1/2 and F n+1/2 i+1/2 at the cell interfaces, for the first-order solution, become: or in more generalized form: which can be used for MUSCL reconstruction. Here, λ is an eigenvalue, as represented in Equation (7). Here, U L and U R represent the values in the left and right cells, respectively, while U * represents the intermediate state of the two solutions at the left and right edges. Guinot [49] suggested solving the Riemann problem at the interface by using a number of exact or approximate solvers. If the location of initial discontinuity is taken at x 0 , then the solution of the problem can be expressed as: where dx dt = −a and dx dt = a in Figure 5 represent the two contact discontinuities separating those three regions [9].   The basic (first-order) scheme uses a piecewise constant approximation for each cell, resulting in a first-order upwind discretization with a cellcentered indexing system, which is indicated as i. For the first-order Godunov's scheme, the above initial interface treatment reduces to U L = U n i and U R = U n i+1 , as represented in Figure 5.
The Riemann problem in its elementary form can be represented as:

∂U ∂t
The solution of Equation (29) consists of m + 1 constant states separated by m number of waves, each of which is associated with a particular eigenvalue λ i . The λ i represents wave speed in case of a water hammer problem ( Figure 6). There is a possibility of three types of elementary wave solutions as shock waves, including contact discontinuity and rarefaction waves, as shown in Figure 7a-c; however, as far as the water hammer is concerned, the shock wave solution dominates.   The underlying challenge associated with the simple first-order scheme is its smearing effect, which is inefficient in shock capturing; thus, specific reconstruction techniques are commonly adopted for specifying initial cell values. Those higher order reconstruction techniques enable higher order resolution of the discontinuity and, hence, provide better results, although sometimes at the cost of higher computational time [53]. Guinot [49] provided insight on the application of those higher order schemes for solving Riemann problems, as the first-order schemes tend to be inaccurate due to strong numerical diffusion. Due to the problem of handling the exact solution in higher order schemes, a number of approximations based on the exact solver were presented. Some notable solvers are the two-rarefaction Riemann solver, two-shock-based Riemann solver, and adaptive-type Riemann solver. The adaptive-type Riemann solver uses a simple Riemann solver for smooth flow, and near the shock location it uses a complex solver in an adaptive way [1,55.56]. Although those schemes were found to be effective in gas dynamics problems, they have thus far received little attention from the water hammer viewpoint.
Apart from the cumbersome exact Riemann solver, Guinot [49] proposed the use of a variety of solvers-including the Harten-Lax-van Leer (HLL) solver, Roe solver, approximate-state solver, and the rarefaction wave solver-in order to convert the hyperbolic problem to a solvable linear problem. Application of the MUSCL (monotonic upstream-centered scheme for conservation laws) technique based on van Leer reconstruction reduced this type of numerical error to a certain extent. A suitable limiter was required for eliminating numerical oscillations in the higher order methods. Since the value of a limiter depends on the solution itself, the solution methods utilizing limiters can be said to be nonlinear. Appropriate reconstruction methods such as MUSCL, the PPM (piecewise parabolic method), and the DPM (discontinuous profile method) were proposed for converting the generalized Riemann problem into a calculation-friendly Shock Wave: Guinot [9] suggested the application of the Rankine-Hugoniot conditions for capturing the shock wave characteristics due to the jump conditions. For the application of the MOC or FVM across the discontinuity generated in shock waves, correction needs to be undertaken, which is known as the equal area rule [2]. A general relationship can be developed based on this principle between the speed of shock and the values of U behind and ahead of it. This relationship can be expressed as follows: where U R , and U L are constant values of the shock at left and right side cell velocity, over time δt the amount of U contained in control volume increased by δU, and a = shock speed. This can further be represented in the form of fluxes, as follows: Equation (32) is known as the Rankine-Hugoniot jump relationship (Figure 7a). The underlying challenge associated with the simple first-order scheme is its smearing effect, which is inefficient in shock capturing; thus, specific reconstruction techniques are commonly adopted for specifying initial cell values. Those higher order reconstruction techniques enable higher order resolution of the discontinuity and, hence, provide better results, although sometimes at the cost of higher computational time [53]. Guinot [49] provided insight on the application of those higher order schemes for solving Riemann problems, as the first-order schemes tend to be inaccurate due to strong numerical diffusion. Due to the problem of handling the exact solution in higher order schemes, a number of approximations based on the exact solver were presented. Some notable solvers are the tworarefaction Riemann solver, two-shock-based Riemann solver, and adaptive-type Riemann solver. The adaptive-type Riemann solver uses a simple Riemann solver for smooth flow, and near the shock location it uses a complex solver in an adaptive way [1,55,56]. Although those schemes were found to be effective in gas dynamics problems, they have thus far received little attention from the water hammer viewpoint.
Apart from the cumbersome exact Riemann solver, Guinot [49] proposed the use of a variety of solvers-including the Harten-Lax-van Leer (HLL) solver, Roe solver, approximate-state solver, and the rarefaction wave solver-in order to convert the hyperbolic problem to a solvable linear problem. Application of the MUSCL (monotonic upstream-centered scheme for conservation laws) technique based on van Leer reconstruction reduced this type of numerical error to a certain extent. A suitable limiter was required for eliminating numerical oscillations in the higher order methods. Since the value of a limiter depends on the solution itself, the solution methods utilizing limiters can be said to be nonlinear. Appropriate reconstruction methods such as MUSCL, the PPM (piecewise parabolic method), and the DPM (discontinuous profile method) were proposed for converting the generalized Riemann problem into a calculation-friendly equivalent Riemann problem. In the following few paragraphs, a summary of various approximations of the exact solution is presented.  Here, the approximate Riemann solver at the interface can be represented as: where * = , and the intercell flux along is obtained as:  Here, the approximate Riemann solver at the interface can be represented as: , and the intercell flux along t is obtained as: Thus, the corresponding HLL intercell flux will be: The Intercell flux F * The main problem of this approximation arises from the underlying two-wave assumption, which becomes unacceptable for a larger system. Toro [1] presented a modification of this technique. Toro et al. [58] named it the HLLC model, where instead of a two-wave model a three-wave model was used. This model results in two-star regions and, hence, can tackle a large number of equations ( Figure 9): and the corresponding HLLC numerical flux:  Roe Method: This method was first documented by Roe [59], and later modified by Roe and Pike [60]. Here, the Jacobian term of the equation has been replaced by a constant Jacobian matrix to approximate the solution ̅ = ̅ ( , ), and then substituted into the original PDE to produce the following approximate Riemann problem: where the cell values are: Figure 9. Waveforms for an approximate HLLC solver.
Thus, upon integration: Roe Method: This method was first documented by Roe [59], and later modified by Roe and Pike [60]. Here, the Jacobian term J of the equation has been replaced by a constant Jacobian matrix to approximate the solution J = J(U L , U R ), and then substituted into the original PDE to produce the following approximate Riemann problem: where the cell values are: Thus, U n+1 i = U n i − ∆t ∆x J U n i − U n i−1 for backward marching, and Applying the vector splitting technique: Thus, the upwinding scheme becomes: and then the intercell fluxes can be represented as: Applying the flux splitting technique to F produces: which then produces the final upwind scheme in conservative form as follows: and the intercell fluxes are: Assuming hyperbolicity, F(U) = J(U)U; thus: Finally, the intercell fluxes for this problem are given as: where the approximation for J is given as J = 1 0 J(θ)dθ, and θ is a parameter varying linearly between 0 and 1 along a straight path connecting U L and U R .
This process involves the conversion of a non-linearized problem to a simple linear matrix algebra-type problem. Hence this method is also known as the Roe linearization technique. In Roe and Pike's [60] approach, the requirement of finding out J was ruled out by applying the direct averaging technique to a set of scalar quantities in order to discover the required eigenvalues and eigenvectors, as given below: , where U is average velocity and given as U = while ρ L and ρ R denote the left and right values of mass density. Thus, the corresponding eigenvalues (λ) and eigen vectors (K) are: Osher Solver: The nonlinear Jacobian matrix J presented in Equation (13) is linearized by splitting as follows, as per Engquist and Osher [61] and Osher and Solomon [62]: where the first term represents positive or zero eigenvalues, and the second term represents negative or zero eigenvalues. Thus, the flux at interface F(U) = F + (U) + F − (U), and ∂F + ∂U = J + (U), ∂F − ∂U = J − (U). The initial data of the problem are denoted by: hence, the corresponding fluxes can be written as: Performing conservation integration and simplifying the Osher fluxes can be given as: These integrals depend on the integral path chosen. A simplistic representation of the integration path is shown in Figure 10. Performing conservation integration and simplifying the Osher fluxes can be given as: These integrals depend on the integral path chosen. A simplistic representation of the integration path is shown in Figure 10. Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL): Instead of having the constant cell value employed in the first-order Godunov's scheme, this higher order scheme assumes a linear cell value using a reconstruction technique. This method was first presented by van Leer [63][64][65] where he introduced a higher order solution of Godunov's scheme using a TVD (total variation diminishing) scheme. The MUSCL scheme is one of the more popular higher order techniques for the finite volume method due to its higher accuracy and its somewhat simpler approach than other available higher accuracy schemes. In this method, Godunov's scheme's constant approximation of cell value is replaced by the reconstructed and slope-limited left and right states. Numerical fluxes are calculated at each of the cell boundaries. Those calculated fluxes act as the inputs of the Riemann solver to provide an averaged solution at the next time step. Hence, the most important part of this solver is the linear reconstruction of the cell values. Typically, the piecewise linear approximation is carried out using a central difference scheme, given as follows: The modified data will then be used to solve the generalized Riemann problem (GRP) to compute the intercell flux F i± 1 2 .
The formation of GRP is represented as:

∂U ∂t
MUSCL Hancock Scheme: A time integration of the hyperbolic portion of the governing equation is now represented as: A modification to this scheme was presented by van Leer [66]. In this method, the evolution of U L i and U R i , which are the flux vectors at the cell interfaces, by time 0.5∆t is given as: The piecewise Riemann problem is represented as: Limiters A suitable limiter function is used to keep the cell values within the cell limit. Therefore, proper limiter functions are implemented as per the requirement. The limiter function limits the spatial derivative values to realistic ones, in order to protect the solution from unreliable convergence and keep the oscillations within acceptable limits. The TVD fluxes with the application of a limiter function can be represented as: , and φ(r) are low-and high-resolution flux limiter functions and the slope of the limiter function, respectively. In principle, the flux is calculated by Equation (63); however, in the case of several equations, the limiter evaluation is not unique. For this reason, the variables are usually limited, and the fluxes are calculated based on the limited variables.
The ratio of successive gradients r i = and φ(r) ≥ 0, which represent that for φ(r) = 0 flux value changes to a low-resolution scheme, and for φ(r) = 1 flux value changes to a high-resolution scheme. Various limiters have been presented in the works of Toro [1] and Guinot [10], which can be used to reconstruct the cell value for water hammer problems.
Pal et al. [53] used MINMOD, since it produced the best results compared to other methods. In the MINMOD method, numerical flux in the predictor step is given as: The cell interface fluxes are obtained using Equations (64) and (65). Figure 11 shows different types of limiter configurations with slopes in the Sweby diagram [67]. Figure 12a,b represents MUSCL reconstruction for a single cell and three cells, respectively. It is wort mentioning that in the case of Riemann solvers, the amount of dissipation depends on the difference between the left and right states. The figure clearly shows that a first-order method produces more dissipation.
The cell interface fluxes are obtained using Equations (64) and (65). Figure 11 shows different types of limiter configurations with slopes in the Sweby diagram [67]. Figure  12a,b represents MUSCL reconstruction for a single cell and three cells, respectively. It is wort mentioning that in the case of Riemann solvers, the amount of dissipation depends on the difference between the left and right states. The figure clearly shows that a firstorder method produces more dissipation.  A piecewise quadratic reconstruction of the equation as given below was presented by Woodward and Colella [68] and Woodward and Colella [69]: where and are associated with first-and second-space derivatives of ( )and = . This is a third-order accurate space discretization, and is popularly known as the piecewise parabolic method (PPM). Another novel method, named the piecewise linear interpolation method (PLIM), was also presented in connection with the analysis of the study of water withdrawal and return after volcanic activity, but has found no use beyond nuclear safety applications [70]. It should be noted that numerical smearing cannot be completely removed in higher order schemes. A nonlinear method with adaptive multiresolution analysis is required to remove the oscillations [71].
Among other approaches for improving problems associated with the non-conservative solution of the hyperbolic problem, the flux-corrected transport (FCT) technique was introduced by Boris and David [72]. This approach essentially employs a finite difference algorithm, and finds its primary application in the conservative solution of the Euler equation in gas dynamics when seeking to capture shocks and discontinuities. The Eulerian FCT algorithm was nicknamed SHASTA (sharp and smooth transport algorithm) and consisted of two stages: (a) the transport stage (predictor stage), and (b) the anti-diffusion stage (corrector stage). This kind of higher order method requires an artificial viscosity term in the momentum equation to prevent, or at least limit, numerical dispersion. However, this method has little application within the water hammer community.
Hwang and Chung [73] utilized the FV technique for water hammer, without neglecting the advective term, and using the conservative form with the density as unknown instead of the pressure head. Szydlowski [74,75] used a second-order accurate FVM scheme for one-dimensional transient flow in a pipeline using a Riemann solver with the Roe linearization technique, in order to analyze the effect of water hammer in a single pipeline system. Szydlowski [75] represented the numerical flux at the interface in a similar way to Equations (49)(50)(51). Increments in pressure and mass flow between boundaries were represented by ∆ , which is calculated by averaging their values at the left and right states. Estimation of ∆ is carried out by using the cell-centered values of and as left and right state cell values, and using MUSCL scheme extrapolation for ensur- A piecewise quadratic reconstruction of the equation as given below was presented by Woodward and Colella [68] and Woodward and Colella [69]: where χ 1 i and χ 2 i are associated with first-and second-space derivatives of U i (x) and ψ = 1 3 . This is a third-order accurate space discretization, and is popularly known as the piecewise parabolic method (PPM). Another novel method, named the piecewise linear interpolation method (PLIM), was also presented in connection with the analysis of the study of water withdrawal and return after volcanic activity, but has found no use beyond nuclear safety applications [70]. It should be noted that numerical smearing cannot be completely removed in higher order schemes. A nonlinear method with adaptive multiresolution analysis is required to remove the oscillations [71].
Among other approaches for improving problems associated with the non-conservative solution of the hyperbolic problem, the flux-corrected transport (FCT) technique was introduced by Boris and David [72]. This approach essentially employs a finite difference algorithm, and finds its primary application in the conservative solution of the Euler equation in gas dynamics when seeking to capture shocks and discontinuities. The Eulerian FCT algorithm was nicknamed SHASTA (sharp and smooth transport algorithm) and consisted of two stages: (a) the transport stage (predictor stage), and (b) the anti-diffusion stage (corrector stage). This kind of higher order method requires an artificial viscosity term in the momentum equation to prevent, or at least limit, numerical dispersion. However, this method has little application within the water hammer community.
Hwang and Chung [73] utilized the FV technique for water hammer, without neglecting the advective term, and using the conservative form with the density as unknown instead of the pressure head. Szydlowski [74,75] used a second-order accurate FVM scheme for one-dimensional transient flow in a pipeline using a Riemann solver with the Roe linearization technique, in order to analyze the effect of water hammer in a single pipeline system. Szydlowski [75] represented the numerical flux at the interface in a similar way to Equations (49)- (51). Increments in pressure and mass flow between boundaries were represented by ∆U, which is calculated by averaging their values at the left and right states. Estimation of ∆U is carried out by using the cell-centered values of U i and U i+1 as left and right state cell values, and using MUSCL scheme extrapolation for ensuring second-order accuracy. Although this method produced predictions in close agreement with analytical results, its large deviation from experimental results demonstrated its inability to capture the additional energy dissipation processes that often influence transient events in real pipe systems.

Flux Vector Splitting (Boltzmann Approach)
Along with Godunov's scheme, which uses the Riemann solver (flux difference splitting) technique, the flux vector splitting method also gained popularity due to its lower computational effort requirements and simpler, more efficient application. Though popular, it suffers from some poorer discontinuity resolution, which makes it a poor candidate for high-shock-wave problems, such as water hammer. Most notable work on developing this method was undertaken by Sanders and Prendergast [76], Steger and Warming [77], and van Leer [78,79]. The upwind differencing is carried out, and the final forms of the flux vectors and corresponding intercell fluxes are represented as: and: For more details about this method, readers are advised to consult Sanders and Prendergast [76], Steger and Warming [77], and van Leer [78,79].
Boundary treatment: For the computation of boundary fluxes at x = 0 and x = L, most popular techniques extend the computational domain beyond the physical domain. These additional cells are called "ghost cells", "fictitious cells", or "virtual cells" [1]. They exist solely for handling boundaries, and do not necessitate extra equations at the boundaries. Two cells, I 0 and I m+1 , are taken as the two extreme point fictitious nodes adjacent to real nodes I 1 and I m , and the corresponding Riemann problem is solved for finding the Godunov fluxes F1 2 and F m+ 1 2 . Some potential physical boundary conditions that can arise for the proposed case are: (1) transmissive boundary, which allows the passage of waves without affecting the boundary used as the far-field boundary; and (2) reflective boundary, or valve boundary, which simulates a reflective impermeable boundary at x = L. The valve boundary condition represents the specific relation of discharge rate to time due to variable valve closure. A detailed description of those boundaries and their treatments is provided in Toro [1], and hence, is excluded here.

Finite Volume Method (FVM) Applications
Zhao and Ghidaoui [80] further investigated applying the FV technique based on mass and momentum conservation principles for water hammer modelling. They considered an MOC-type system of two pipes connected in series with a valve and a reservoir. They solved the numerical flux by the exact solution of the Riemann problem using both firstand second-order Godunov's schemes. The advective term was calculated by the arithmetic mean in order to solve the Riemann problem. They considered the conventional case as zero, and then compared the solutions. Finally, they concluded that although the first-order MOC, like the Godunov's scheme, showed high numerical dissipation with reduced accuracy, the second-order Godunov's scheme with a MINMOD limiter significantly increased the accuracy and efficiency. Although a second-order formulation was used for interior nodes, a first-order formulation was applied at the boundaries. For wave propagation, the firstorder interpolation on the boundaries cannot be of vital importance, although it certainly affects the results. Later, the same authors [51] sought to improve the scheme by using the MUSCL-Hancock method with a virtual cell technique at the boundary, thus maintaining a second-order accuracy throughout the computational domain. In this scheme, the left and right states of discontinuity at the ith cell are expressed as follows: where ∆t and ∆x are temporal and spatial steps. This method showed improved accuracy and efficiency for large pipe networks compared to the previous method [51,80,81]. Yazdi et al. [82] used a similar approach for the second-order Godunov's scheme in order to develop a one-dimensional coupled model based on Riemann solutions with a similar type of arrangement, and compared the results with analytical solutions. They concluded that the second-order scheme agreed closely with the analytical method when the Courant number (C r ) was equal to one, but showed some numerical dissipation for lesser Courant numbers. Leon et al. [83] proposed a virtual cell technique for the treatment of boundaries. Later, Leon et al. [84] further continued their work and introduced an FV shock-capturing scheme for both one-and two-phase flows. They compared the results between the FV model and the available MOC model for single-phase flow, and concluded that the MOC is suitable for C r = 1, but for C r < 0.95 the efficiency and accuracy of the shock-capturing FV scheme proved to be significantly better. According to Leon et al. [84], a comparison of numerical results between the FV scheme and the fixed-grid MOC for two-phase flow demonstrated significant advancement in the efficiency of the FV method.
Kerger et al. [85] developed a unified numerical model for all types of flows, such as pressurized flow, open-channel flow, and mixed flow. In their method, Kerger et al. [85] used the Preissmann slot technique [86] with a first-order Godunov's scheme based on the exact Riemann solver to solve the pressurized flow. The comparison between analytical, experimental, and numerical results was good, but numerical oscillations were present when the wave celerity was large. A 2 × 2 system for water hammer was solved by Guinot [46] using a higher order piecewise parabolic scheme along with Godunov's scheme for a simple single-pipe system. A six-point-in-space model was proposed along with a number of boundary condition treatments with a virtual cell technique. Although it showed improved accuracy over the MOC, it proved that the ill-treated boundary condition could also result in degraded accuracy for highly accurate and efficient schemes [46]. Daude et al. [87] used the FVM model for simulating water hammer events with column separation. This experimental case corresponds to water hammer with moderate cavitation.
The good agreement between the numerical pressure history at the valve and its experimental counterpart is clearly shown in Figure 13. The primary pressure peak was generated by the instantaneous closure of the valve, and has a value equal to the Joukowsky value ≈ 6.83 bar. Afterwards, cavitation occurs at the valve due to the reflection of the rarefaction wave at t = 57 ms. The collapse of the vapor pocket at t = 120 ms generates the secondary water hammer wave. The superimposition of the primary and the secondary water hammer transients produces a pressure peak with a magnitude of 10.1 bar at t ≈ 171 ms, which is higher than the Joukowsky pressure. The end of this pressure pulse corresponds to the arrival of the second water hammer wave at t ≈ 177 ms. This very short duration of 6 ms pressure pulse demonstrates the accuracy of the FVM. The three last pressure peaks were overestimated due to the inclusion of only steady-state friction. However, the timing of the pressure waves was well retrieved by the FVM, even for the later water hammer events. In a further study, a validation exercise was undertaken for first-order and secondorder FVMs as per the methodology presented in Pal et al. [53]. A frictionless pipe system with a large reservoir at the upstream and a valve at the downstream was considered. The internal diameter and length of the pipe were taken as 0.04 m and 1000 m, respectively. Initial steady velocity (i.e., prior to valve closure) was taken as 1.0 m/s, while constant upstream reservoir head was taken as 10 m. The wave celerity was taken as 1000 m/s, with a mass density of water = 1000 kg/m 3 . The domain was spatially discretized into 100 segments. Total simulation time was considered to be 20 s after the valve closure, with a constant time step ∆ of 0.01 s. The same space and time discretization was adopted for the MOC. The simulation results are presented in Figure 14. Both the first-order and second-order FVM (MUSCL-Hancock method) were considered, and their results were the same. Therefore, both are represented as FVM in Figure 14. The computational times of the first-order FVM and the standard MOC are comparable. However, the computational time of the second-order FVM is an order of magnitude more than the MOC. The comparisons in Figure 14 demonstrate that the results of the MOC and the FVM are the same for = 1. However, for values less than unity, such as values of 0.9 and 0.8, the FVM maintains excellent phase accuracy; however, as the Courant number decreases, the accuracy decreases, as shown by the rounding of the sharp corners of the exact solution, although the error is small. This particular characteristic of the FVM is highly advantageous in cases where wave celerity changes along the pipeline due to either air pockets or, for instance, a change in the pipe wall thickness or pipe diameter.
These results also demonstrate that the first-order FVM is as good as the MOC, which was also demonstrated by Pal et al. [53]. In addition, the FVM method [53] was also validated with experimental data, as demonstrated by USC-2 in the work of Reddy et al. [88].  In a further study, a validation exercise was undertaken for first-order and secondorder FVMs as per the methodology presented in Pal et al. [53]. A frictionless pipe system with a large reservoir at the upstream and a valve at the downstream was considered. The internal diameter D and length L of the pipe were taken as 0.04 m and 1000 m, respectively. Initial steady velocity V 0 (i.e., prior to valve closure) was taken as 1.0 m/s, while constant upstream reservoir head H 0 was taken as 10 m. The wave celerity was taken as 1000 m/s, with a mass density of water ρ = 1000 kg/m 3 . The domain was spatially discretized into 100 segments. Total simulation time was considered to be 20 s after the valve closure, with a constant time step ∆t of 0.01 s. The same space and time discretization was adopted for the MOC. The simulation results are presented in Figure 14. Both the first-order and second-order FVM (MUSCL-Hancock method) were considered, and their results were the same. Therefore, both are represented as FVM in Figure 14. The computational times of the first-order FVM and the standard MOC are comparable. However, the computational time of the second-order FVM is an order of magnitude more than the MOC. The comparisons in Figure 14 demonstrate that the results of the MOC and the FVM are the same for C r = 1. However, for C r values less than unity, such as C r values of 0.9 and 0.8, the FVM maintains excellent phase accuracy; however, as the Courant number decreases, the accuracy decreases, as shown by the rounding of the sharp corners of the exact solution, although the error is small. This particular characteristic of the FVM is highly advantageous in cases where wave celerity changes along the pipeline due to either air pockets or, for instance, a change in the pipe wall thickness or pipe diameter. presented in Figure 15. These results demonstrate that the dissipation in exper pressure is quite high, whereas numerical modelling with steady-state friction un dicted the decay. This can be addressed by implementing either IAB (instantaneou eration-based) models or by weighting function-based unsteady friction models. the similarities between the FVM data, MOC data, and experimental data suggest t ability of the FVM for modelling water hammer flows.  These results also demonstrate that the first-order FVM is as good as the MOC, which was also demonstrated by Pal et al. [53]. In addition, the FVM method [53] was also validated with experimental data, as demonstrated by USC-2 in the work of Reddy et al. [88]. The pipe system consisted of 156-m-long copper pipe of 0.025 m diameter. The estimated values of wave velocity and the Darcy-Weisbach friction factor were 1050 m/s and 0.026, respectively. The experimental flow conditions were initial steady-state flowrate and inlet pressure head equal to 4.8 × 10 −5 m 3 /s and 7.15 m, respectively. The valve was rapidly closed in 0.14 s, and it was modelled by the MOC and the first-and second-order FVMs. Since the MOC and FVM results are identical, only the FVM and experimental data are presented in Figure 15. These results demonstrate that the dissipation in experimental pressure is quite high, whereas numerical modelling with steady-state friction underpredicted the decay. This can be addressed by implementing either IAB (instantaneous acceleration-based) models or by weighting function-based unsteady friction models. Finally, the similarities between the FVM data, MOC data, and experimental data suggest the suitability of the FVM for modelling water hammer flows.

Two-Phase Flow
Multiphase flow in pipe systems is commonly observed in various industrial fi whether requiring of different types of fluids or different phases of the same fluid. In cases, both steady-state and transient predictions of flow and hydraulic shocks are cru A critical problem associated with the application of multiphase analysis is its high putational expense [89]. Taitel [90] presented a model that treats only the liquid as transient, whereas the gas is assumed to be in a quasi-steady state. However, such a m is invalid for large shocks. The early model was later modified by Minami and Sho [91], who advocated for an implicit computational scheme, and suggested an impr flow regime transition model. Li [92] considered the continuity equation to be trans but employed the momentum equation as a quasi-steady-state representation. Choi [93]  Pressure head (m) Figure 15. Comparison between the FVM and experimental pressure histories [88].

Two-Phase Flow
Multiphase flow in pipe systems is commonly observed in various industrial fields, whether requiring of different types of fluids or different phases of the same fluid. In such cases, both steady-state and transient predictions of flow and hydraulic shocks are crucial. A critical problem associated with the application of multiphase analysis is its high computational expense [89]. Taitel [90] presented a model that treats only the liquid as fully transient, whereas the gas is assumed to be in a quasi-steady state. However, such a model is invalid for large shocks. The early model was later modified by Minami and Shoham [91], who advocated for an implicit computational scheme, and suggested an improved flow regime transition model. Li [92] considered the continuity equation to be transient, but employed the momentum equation as a quasi-steady-state representation. Choi et al. [93] suggested a two-phase transient model with a modified continuity equation, and treated momentum as an extended drift-flux equation. Oloruntoba and Kara [94] presented a transient model based on the Navier-Stokes equations, without neglecting the convective term.
One of the most important applications of two-phase transient analysis is the safety analysis of nuclear reactors. Water hammer and two-phase transients can induce failure and accidents at such plants, making transient analysis and prediction of utmost importance. For the analysis of the safety of nuclear reactors, a number of computational techniques were introduced in the early 1980s. The two-phase flow equations represent a comprehensive presentation of the flow, but lead to a non-hyperbolic governing system. Nevertheless, the hyperbolicity can be established with a special treatment of the phase interaction terms, if one assumes equal velocities or equal temperatures and the number of equations drops. With equal velocities, the system of governing equations becomes hyperbolic. A few notable examples among those codes used for thermohydraulic analysis are THERMIT, TREMOLO, RELAP5, RETRAN, CATHARE, TRACE, etc. [95][96][97][98][99][100]. The basic equation of multiphase flow used in reactor safety analysis is in the following form [101]. The continuity equation: ∂ρ a ∂t M a (71) where mass conservation is represented as ∑ . M = 0 and a = l, g, where l and g are liquid and gas phases, respectively.
The momentum equation with an assumption of equal interfacial pressures is as follows: The study of nuclear reactors requires consideration of the energy conservation, as follows: ∂(ρ a S a ) ∂t where ρ a is the partial density of phase a, α a is the volume fraction of phase a, v a is the phase velocity, M a is the rate of formation of phase a, p is the pressure, A ab is the interphase area per unit volume, D ab g is the drag coefficient multiplied by the density and the absolute value of the velocity difference,v a is intrinsic velocity associated with the source term, F a x is the wall friction for phase a, S a is the entropy per unit mass, T a is the phase temperature, and ω a is the rate of dissipation of energy due to friction and mass transfer. The right-hand sides of the total energy equations for both of the phases are complicated, and contain terms rising from Equation (72) and heat transfer. All of these are grouped to the RHS term of Equation (73). Since the original variable is ρ a e + (v a ) 2 /2 , and the pressure work is usually included in the convection term to form a phase enthalpy, the actual equations that are solved may differ from Equation (73). The equation of state for phase a is given as ρ a = ρ a (S a , p).
RELAP (1)(2)(3)(4)(5), developed by the Nuclear Regulatory Commission, is widely considered to be one of the best estimating codes among those in use for light water reactors. It is essentially a one-dimensional, two-phase, non-homogeneous, and non-equilibrium model based on a partial implicit scheme that provides prime importance to mass and energy conservation over momentum conservation. RELAP (1)(2)(3)(4)(5) is popular due to its specific adaptation to this application.
Tiselj and Petelin [102] used a shock-capturing two-phase model based on the RELAP5 code, consisting of non-hyperbolic partial differential equations of mass momentum and energy balance. They used a basic averaging technique for the solution of the Jacobian, using conservative and non-conservative basic variables, and they used the method for the solution of a two-phase shock tube problem.
Based on the method developed by Tiselj and Petelin [102], Melikhov et al. [103] and Gale et al. [104] proposed a state-of-the-art WAHA (water hammer) code-which is a onedimensional water hammer model-to solve two phases, and they found good agreement with experimental measurements. The WAHA code is a six-equation model based on mass, momentum, and energy conservation. The model is mainly used for calculating hydraulic forces on the piping system during a water hammer event.

Development of Other Approaches
Various numerical schemes have also been used to analyze water hammer problems. For example, an MOC-like scheme inspired by the graphical method that tracks wavefronts is often called the wave plan or wave characteristic method [105]. Since the accuracy of this approach is first-order both in space and in time, a finer discretization is required in order to achieve high accuracy, but internal sections are largely avoided. The key assumption associated with this method is that friction is concentrated in a virtual orifice located at the midpoint of each pipe. Verwey and Yu [106] used a fourth-order scheme to describe steep pressure fronts accurately. Using three time levels and two spatial points better represented the boundary conditions. Chen et al. [107] modelled water hammer by using optimal boundary control for the nonlinear hyperbolic PDEs describing closed-conduit flow. They used a method of lines to solve the converted ordinary differential equations (ODEs), thus reducing numerical dissipation and associated fluctuations.

CFD Schemes
Zhang et al. [47] improved a Lax-Wendroff method (LWM) [1]-based approach for a gravitational reservoir pipe system with air entrainment, without neglecting the convective acceleration terms. Comparison of their experimental data with numerical data showed that although the LWM scheme accurately predicted the maximum pressures similarly to the MOC, it underpredicted the attenuation of the water hammer waves-something that they attempted to control through an additional friction function. Szymkiewicz and Mitosek [108] used the Galerkin-based finite element method (FEM) for simulating water hammer problems. The FEM is not favorable for modelling water hammer problems, for reasons such as (1) stability criteria still needing to be satisfied, (2) accuracy being no better than that of the MOC, (3) requiring excessive computational resources, and (4) difficulty to program. However, FEM is the most sought-after method for modelling water hammer problems with fluid-structure interaction (FSI).
Modelling of the water hammer phenomenon for dynamic characteristics within a turbine is impossible using the classical 1D MOC or FVM schemes, and such applications require more extensive 3D grids and turbulence models. Several commercial pieces of software for turbulence modelling available today can be effectively used for this type of study. Some well-known and well-applied turbulence models currently in use are FLUENT and CFX (https://www.ansys.com (accessed on 1 June 2021)). Commonly, user-defined function and 1D MOC coupled with a 3D CFD model are used for providing the complex input and boundary conditions for turbine transient problems [109]. Yuan et al. [110] presented another approach for coupling 1D MOC with a 3D CFD model, and validated surge tank oscillation, which was further modified by Fan et al. [111], Wu et al. [112], and Zhang et al. [113]. Yan et al. [114] and Yin et al. [115] used CFX and FLUENT (https://www.ansys.com (accessed on 1 June 2021)) for 3D simulation of compressible water hammer problems, in which a hydroelectric generator suddenly rejected its load in order to safeguard itself from an accident, and created a hydraulic transient in the penstock. Zhang et al. [116] modelled a transient associated with Francis turbines during load rejection by using a novel 1D-3D coupled method for a compressible liquid model, and compared their data effectively with standard classical MOC models.

Lagrangian-Based Numerical Methods
Recently, Lagrangian-based numerical methods such as the lattice Boltzmann method (LBM) and smoothed-particle hydrodynamics (SPH) are gaining popularity for modelling unsteady flows. The LBM and SPH models are meshfree, in which governing equations are solved at any point in the space-time plane to find local solutions without explicitly using previously stored data. In particular, the LBM is a modern approach in CFD, which belongs to a class of kinetic theory. The LBM is suitable for modelling (1) multiphase flows including tiny droplets and bubbles, (2) flow through porous media, (3) flow through complex geometries, (4) coupled flow with heat transfer, and (5) cardiovascular hemodynamics. The LBM was considered by Cheng et al. [117] for unsteady flow, and later a practical and simple implementation was carried out by Wu et al. [118]. Budinski [119] modelled a one-dimensional water hammer using the LBM, in which the grid was independent of the Courant number. To achieve this, the governing equations were transformed using an alternative coordinate system suitable for an adaptive grid. The LBM simulation was validated using the MOC for water hammer resulting from pump tripping, with excellent agreement. Furthermore, Louati et al. [120] modelled a two-dimensional water hammer using the LBM, and compared the results with a fifth-order FV scheme. According to Louati et al. [120], the higher order FDM and FVM are challenging to implement with complex boundary conditions such as pumps, turbines, and valves.
Overall, LBM schemes are a viable alternative. Louati et al. [120] observed that the LBM is computationally fast and equally accurate to the MOC for low-frequency water hammer waves. However, the LBM becomes unstable and inaccurate for high-frequency water hammer waves. The LBM is faster and more accurate, whereas the SPH method is more stable and has broader applications, such as one-dimensional wave propagation, slug flow, sediment flushing from sluice gates, and three-dimensional dam breaking.
Hou et al. [121] were the first to use SPH for simulating the rapid filling of a pipeline. The 1D water hammer equations were used along with a dynamic SPH pressure boundary condition. This method reduced the problem of boundary deficiency associated with conventional SPH methods. In this approach, a set of particles is introduced at the upstream boundary, depending on the smoothing length, with predefined velocity and pressure. As those particles are flowing through the system, the next set of particles is generated. The same procedure is also adopted for the downstream boundary condition. A cubic spline kernel function approximates the spatial derivatives. Good agreement is observed between SPH and rigid column theory, especially in the deceleration phase of the filling. Hou et al. [122] utilized a corrective smoothed-particle method (CSPM) to simulate water hammer in a simple reservoir, pipe, and valve system. A corrective kernel estimate of SPH approximated the spatial derivatives in the water hammer conditions. For the transient derivatives, the Euler-forward time joining calculation was utilized. An artificial viscosity term was added to the momentum equation in order to prevent numerical dispersion. Hou et al. [122] found that the CSPM outcomes are in good agreement with the results acquired by the MOC for three different cases. A parametric examination provides knowledge about the impacts of particle conveyance, smoothing length, and kernel work. Although this method is inefficient for traditional water hammer systems, it is a key step representing the free-surface water hammer problem. Furthermore, Hou et al. [123] used a second-order accurate CSPM to model dynamic interaction between air bubbles and water in a dead-end pipe. In this method, a new water particle is introduced at the upstream when the most upstream particle distance is more than one particle size from the inlet. Similarly, as soon as the most upstream particle leaves the inlet, it is deleted in backflows. This method can introduce an error due to the mispositioning of the upstream particle, but the authors emphasize that this effect will be negligible for smaller particles. Hou et al. [123] argued that this method is accurate, efficient, and straightforward for a water hammer with a nonlinear downstream boundary condition, such as a compressible air pocket at the dead end. Although a few researchers have obtained good results with SPH for modelling the water hammer, the method is still awkward for the implementation of complex boundary conditions. Furthermore, there are a few advanced schemes-namely, FEM, LBM, etc.-that can be used for transient analysis, but those are not of much relevance to practicing engineers.

Two-Dimensional Numerical Schemes
Although one-dimensional schemes are easy to understand and implement, and when used appropriately achieve excellent accuracy, stability, and monotonicity, their use in multidimensional problems with anisotropic turbulence is difficult and error-prone. In addition, 1D schemes for complex devices like pumps and turbines, and even unsteady friction, have their limitations. Therefore, in recent times, considerable effort has been invested into developing two-dimensional transient pipe flow models for calculating Reynolds shear stress.
Shear stress τ in Equation (17) is a summation of viscous and Reynolds stresses. The Reynold stress −ρu v is dominant in turbulent flows, and it can be computed using an appropriate turbulence model based on RANS equations. For water hammer applications, the most popular method for turbulence modelling is an algebraic model, where the eddy viscosity is expressed as a function of the flow field. Details of various turbulence models based on the same principle for water hammer phenomena, their theories, and their constraints can be found in the works of Vardy and Hwang [12], Eichinger and Lein [124], Silva-Araya and Chaudhry [17], Silva-Araya and Chaudhry [125], Pezzinga [14], and Ghidaoui et al. [126]. Vardy and Brown [127] followed Zielke's methodology by assuming an idealized turbulent flow with eddy viscosity varying linearly in the outer layer, which enabled an approximate analytical relation between the decay of wall turbulence and the sudden change in the velocity. Vardy and Brown [127] used the MOC and a five-layer turbulence model in the radial direction. Silva-Araya and Chaudhry [17] proposed a two-layer turbulence model in which Prandtl's mixing-length approximation was used for estimating eddy viscosity at the inner and outer layers. Pezzinga [14] implemented the mixing-length hypothesis in a quasi-two-dimensional (2D) model in order to simulate transient flows in simple pipe networks. A comparison between two popular algebraic models-two-and five-layer turbulence models-was presented by Ghidaoui et al. [126,128], showing that the two-layer model is preferable for unsteady flow. Vardy and Brown further improved their method by assuming the bilinear approach of eddy viscosity in the outer region, i.e., steady increase in the outer region close to the wall and an approximately uniform eddy viscosity in the core region. Moreover, Zhao and Ghidaoui [80] improved the speed and stability of the Vardy and Hwang [12] model by decoupling the overall solution matrix into two tridiagonal matrices consisting of pressures and velocities, respectively.
Furthermore, Zhao and Ghidaoui [122] developed a k − ε turbulence model based on the work of Fan et al. [105] for water hammer flows. Zhao and Ghidaoui [129] used the MOC to solve continuity and momentum equations for water hammer equations. Recently, numerical models based on a combination of the Runge-Kutta methods and two-equation turbulence models are receiving attention. Wahba [130,131] modelled twodimensional unsteady pipe flow using mixing-length Baldwin-Lomax and Cebeci-Smith turbulence models, along with the Runge-Kutta method, to integrate mass, momentum, and central differencing for spatial derivatives. Riasi et al. [132] used the fourth-order Runge-Kutta scheme along with k − ω and k − ε turbulence models for studying the behaviour of the unsteady turbulent flow. Both of these models reported good consistency with the experimental results. Later, some effort was also made in using LES (large eddy simulation), but this 3D-based model requires significantly smaller cells [133]. Direct numerical simulation (DNS) resolves all scales, including Kolmogorov scale eddies, without using any empirical model. DNS solves full Navier-Stokes equations, and hence, is computationally very intensive, although recent developments in computational power and parallelization of processing are encouraging for the use of DNS. However, DNS remains applicable only to small Reynolds numbers (Reynolds numbers less than 5000), small pipe lengths, and small diameters. Therefore, the use of the DNS technique has not gained sufficient popularity in application to transient analysis. A few studies that deal with the techniques include those by Eggels et al. [134] and Ghosh et al. [135], but most are not directly related to transient analysis. In order to obtain a representative ensemble, a number of simulations are required, which is often computationally infeasible. In recent years, some researchers [136][137][138] have used 3D CFD models for solving valveinduced water hammer based on MOC code. Mahdizadeh [139] used a modified flux wave approach and two-dimensional finite volume method with a ν 2 − f turbulence model, and got identical results for experimental damping. Finally, it is clear from the literature presented here that two-and three-dimensional models can accurately predict frictional damping. Mandair et al. [140] developed a tool called integrated total energy (ITE) for CFD modelling to analyze energy dissipation during water hammer. Figure 16 illustrates a brief review of all of the turbulence models explained previously for a better understanding. The available turbulence models for the Reynolds shear stress are based on either steady-state or unsteady-state eddy viscosity distribution. The research approaches [130][131][132][133][134][135]141,142] based on time-varying or time-invariant eddy viscosity distributions are given in Figure 16.
all of the turbulence models explained previously for a better understanding. The available turbulence models for the Reynolds shear stress are based on either steady-state or unsteady-state eddy viscosity distribution. The research approaches [130][131][132][133][134][135]141,142] based on time-varying or time-invariant eddy viscosity distributions are given in Figure  16.

Modelling of Unsteady Flows in Pipe Networks
The presence of complex configurations of pumps and pipe connections in networks gives rise to different types of flow disturbances. Thus, it is often vital to consider pipe network design along with the study of the unsteady flow itself. Karney

Modelling of Unsteady Flows in Pipe Networks
The presence of complex configurations of pumps and pipe connections in networks gives rise to different types of flow disturbances. Thus, it is often vital to consider pipe network design along with the study of the unsteady flow itself. Karney and McInnis [143,144] presented a simple algebraic framework for a simple pipeline network system, and proposed an explicit algorithm known as EED (external energy dissipator). Cabrera et al. [145] suggested that there are two types of basic models available for network analysis: (a) dynamic models (inertial and non-inertial), and (b) static or steady-state models. Among those models, only the elastic inertial model is the best for modelling water hammer in pipe networks. Ani and Khayatzadeh [146] introduced an interpolation technique (based on the finite difference method) coupled with a fixed-grid MOC to capture the exact solution where the Courant's condition does not satisfy the simple MOC at any specific location of the pipe network system. Gad and Mohammed [147] proposed an efficient technique for simplifying the water supply network using two different model networks based on numerical simulations carried with Water Hammer and Mass Oscillation (WHAMO) software (www.cecer.army.mil/usmt/whamo/whamo.htm (accessed on 1 June 2021)). They suggested hydraulic equivalence and demand concentration for the simplification.
Rigid water column (RWC) models are sometimes useful for modelling unsteady, incompressible flow in pipe networks. Their general applicability falls between the inertia model and the quasi-steady model [148,149]. Todini [150] provided an effective quasisteady model GGGA (generalized global gradient algorithm), but it was not helpful for inertia-enabled networks. Similar works were also carried out by various researchers including Shimada [151], Onizuka [152], Holloway [153], and Ahmed [154], and work on the application of RWC to a pressurized pipe was conducted by Islam and Chaudhry [155] and Nault and Karney [148,149]. Todini and Rossman [156] produced systems of nonlinear formulations for looped water distribution networks, and they used methods such as the NR (Newton-Raphson) method and LT (linear theory) to covert those into linear formulations, which show faster convergence rates than the nonlinear formulations. Nault and Karney [148] presented a modification of the global gradient algorithm (GGA) to produce the RWC GGA algorithm. The RWC CGA considers the inertial effects of pressured flow and improves the model applicability for simulating water hammer in the pipe networks. However, this model showed a higher computational cost than the conventional GGA method. Nault and Karney [149] proposed a combined GGGA and RWC GGA, known as HGGA, along with a variable time step approach to include the inertial effect of pressurized flow. This method is computationally costly, but it is highly accurate for dynamic hydraulics.
Jung and Karney [157] compared different water hammer models-including rigid water column analysis, quasi-steady analysis, and the Joukowski approach-based on the physical properties of the models, predictions of actual data, and stability and accuracy. Finally, they proposed some guidelines for the suitability of using any particular model according to the need. Nault et al. [158] proposed a flexible generalized characteristics method (GCM) by coupling a high-resolution MOC and a low-resolution WCM (wave characteristics method), which can generate accuracy similar to the MOC at a better computational cost. Later, Nault and Karney [159] proposed a novel adaptive hybrid transient model (AHTM) based on comprehensive GCA (CGCA) for accurate and efficient modelling of water distribution networks.

Future Research Directions for Water Hammer Modelling
Water hammer responses are now being used for leak detection and state assessment of pipe conditions in water pipeline systems. With the invention of efficient numerical schemes and CFD, a variety of commercial software such as RELAP5-3D (relap53d.inl. gov/ (accessed on 1 June 2021)), WANDA (deltares.nl/en/software/wanda/ (accessed on 1 June 2021)), HYTRAN (accutech2000.com.au/hytran (accessed on 1 June 2021)), PIPENET (sunrise-sys.com/pipenet (accessed on 1 June 2021)), KYPIPE (kypipe.com/ (accessed on 1 June 2021)), TRANSAM (hydratek.com/expertise (accessed on 1 June 2021)), and HAMMER (bentley.com/ (accessed on 1 June 2021)), etc., are presently available for practical use. Most are based on MOC-type numerical schemes for simplicity. Nevertheless, the MOC approaches have the notable disadvantages already introduced, such as nonconservation of mass momentum for non-unity and non-uniform courant numbers, and difficulty in including the convective acceleration term. Besides these concerns, complex device characteristics (e.g., for pumps and turbines), and even unsteady friction-not to mention certain air-water interactions-can be thought of as hints that simplified models themselves have limitations that can confound 1D analysis. By contrast, the FVM approach has key advantages, including an intrinsic satisfaction of the conservation of mass and momentum equations, as well as improved shock-capturing. However, the application of FVM methods has been limited due to their complexity and comparative newness when applied to water hammer problems.
The recent advances in genetic algorithms (GAs) and other evolutionary codes have created a huge opportunity for hydraulicians in their application to various unsteady flow problems. The GA method can achieve high efficiency in modelling the optimization problems due to its ease of application and its controlled computational effort. Tang and Karney [160] successfully applied GA techniques using data from two pump trip tests, and estimated wave speed and roughness factors. In addition, extensive analysis of pressure histories at various locations of pipe systems can be useful for the rehabilitation of pipes, calibration of pipe networks, leakage, and blockage detection by applying inverse analysis, which is also gaining interest.
Parallel computing techniques are gaining increased attention due to their potential to reduce computational time and costs. With the increasing availability of supercomputers, parallel computing has become progressively easier, cheaper, and more widely applied. Graphical processing units (GPUs) have attractive capabilities for large-scale computation. Wu et al. [161] successfully used GPUs for FSI applications, while Bonelli et al. [162] found GPUs to be up to 156 times faster than traditional CPU approaches. Zhang et al. [163] proposed and implemented a multiple-relaxation-time lattice Boltzmann model (MRT-LBM) on a single GPU. Griebel et al. [164] simulated three-dimensional multiphase flows on multiple GPUs. Meng et al. [165] successfully completed parallel integration of the MOC with GPUs, thus confirming its attractiveness for this application.
The published FVM schemes are mostly either first-order or second-order accurate. Third-order accuracy can be achieved by applying the piecewise parabolic method (PPM) technique, which may allow a more precise parabolic reconstruction for solving the water hammer problem using Godunov's scheme. This model requires further investigation in terms of computational speed, reliability, and applicability. A flexible error-detection technique is needed in order for higher order schemes to achieve excellent results. In addition, the feasibility and the usefulness of FVM schemes for their applicability in currently available commercial software need to be extensively studied for the advancement of software capabilities.
Future development in numerical techniques has the potential to reduce extensive and expensive testing. Efficient numerical models may enable leakage and blockage detection, and knowledge of whether to assess deposit thickness, rehabilitation needs in pipes, and loss of non-revenue water (NRW) or pipe wall condition. Two-and three-dimensional finite volume modelling with the ability to capture turbulence-driven vortices in transient flow might pave the way for groundbreaking studies of fluid-structure interaction. Possible problems such as pipe incrustation, wall erosion, and water contamination in water networks should be better analyzed, so that suitable preventive and rehabilitation measures can be taken. Thus, the further advancement of numerically more accurate models continues to have the potential to save costs and avoid trouble in the pressing pipe design, management, and operational domains.

Summary and Conclusions
Transient analysis is crucial for the design, operation, and risk analysis of pipeline systems. This work reviews historical developments, with a special emphasis on the more recent finite volume method due to its high potential of providing more accurate results related to the numerical prediction of transient flow. This paper briefly recaps the advancement of the MOC scheme, from fixed-grid MOC schemes to various interpolation techniques. Detailed analysis of water pipe networks might require even more sophisticated 1D, higher dimensional, and CFD schemes in order to face present and future challenges. In this direction, the application of the FVM to transient flow is concisely explained, considering Godunov's first-and second-order schemes for the transient flow problem. Numerical and experimental validation demonstrates that the first-order Godunov's scheme is as good as the MOC in terms of accuracy and computational speed.
In addition, recent developments related to the modelling of unsteady flows in complex pipeline networks are briefly summarized. Overall, this paper provides a systematic representation of past developments and directions for future research related to water hammer modelling in pipe systems.
Author Contributions: S.P.: writing and original draft preparation; P.R.H.: writing-review and editing, project administration, and funding acquisition; B.W.K.: writing-review and editing, and funding acquisition. All authors contributed to the work. All authors have read and agreed to the published version of the manuscript.