Development of Smoothed Particle Hydrodynamics based water hammer model for water distribution systems

ABSTRACT

Smoothed Particle Hydrodynamics (SPH) method is used to solve water hammer equations for pipeline systems due to its potential advantages of easily capturing column separation and slug impact. Currently, the SPH-based water hammer model has been only developed to simulate single pipe flow with simple boundary conditions. It is still a challenge to apply the SPH-based water hammer model to practical water distribution systems (WDSs). To address this issue, this study develops a complete SPH-based Water Hammer model for Water Distribution System (SPH-WHWDS). Within the proposed method, the complex internal and external boundary condition treatment models of the multi-pipe joint junction and different hydraulic components are developed. Buffer and mirror particles are designed for boundary treatment coupling with the method of characteristics (MOC). Two benchmark test cases, including an unsteady pipe flow experiment and a complex WDS, are used to validate the proposed model, with the data from the experimental test in the literature and the simulation results by the classical MOC. The results show the proposed SPH-WHWDS model is capable to simulate transient flows with accurate and robust results for pipeline systems, which may provide further insights and an alternative tool to study water hammer phenomena in complex WDSs.

KEYWORDS:

• Smoothed Particle Hydrodynamics (SPH)
• water distribution systems
• water hammer
• boundary condition treatment

1. Introduction

Transient flow is a common phenomenon in practical water distribution systems (WDSs) (Alawadhi & Tartakovsky, 2020; Kim, 2022), which is usually caused by the rapid change of one or more boundary conditions in the system, such as pump start-up or shut-off, instantaneous valve operation, etc. This phenomenon, commonly known as pressure surge or water hammer, is characterised by fast and violent fluctuations in pressure and velocity in the pipeline systems. If not managed properly, it may damage the pipelines and associated hydraulic machinery. Therefore, the prediction, analysis, and prevention of water hammer in engineering practice are of great significance to pipeline design, construction and operation. To this end, an accurate simulation of water hammer becomes crucial to such practical applications and theory development.

Several numerical approaches are available to simulate water hammer problems, among which the one-dimensional (1D) models and methods are widely applied because of their efficiency for simulations and convenience for implementation (Duan et al., 2020). Wylie and Streeter (1978) used a set of two hyperbolic partial differential equations (PDEs) to describe the transient flows. These equations are solved by using the method of characteristics (MOC), which transforms the governing pair of PDEs into four coupled ordinary differential equations (ODEs). Wiggert and Sundquist (1977) predicted pipeline transients with fixed-grid intervals and demonstrated the effects of interpolation and grid size on numerical attenuation and dispersion, while Goldberg and Wylie (1983) proposed the interpolations in time had more advantages for applying the MOC to wave problems in hydraulics. Sibetheros et al. (1991) utilised the MOC with spline polynomials for interpolations on water hammer analysis, which was more accurate than linear interpolations and second-order explicit finite difference techniques. Later, Karney and Ghidaoui (1997) introduced a hybrid interpolation approach as a pre-processing step for discretization in pipe networks. In addition to the MOC, the finite-difference method (FDM) and the finite-volume method (FVM) have also been utilised to solve transient equations. Chaudhry and Hussaini (1985) analysed a second-order accurate explicit finite-difference scheme for water hammer simulation, while Samani and Khayatzadeh (2002) used the FDM coupled with a fixed-grid MOC to handle the disproportionate pipes in networks. Besides, Zhao and Ghidaoui (2004) solved the water hammer problems with first and second-order explicit finite volume Godunov-type schemes. Unlike the scheme of Zhao and Ghidaoui (2004), Leon et al. (2007) proposed a second-order accurate FV scheme for modelling water hammer flows by using the conservative form of the water hammer equations. Later, Leon et al. (2008) introduced an FV shock-capturing scheme for simulating one-phase and two-phase water hammer flows. Daude et al. (2018) utilised the quasi-1D FVM model to simulate water hammer events with cavitation-induced column separation. Furthermore, Zeng et al. (2022) proposed an elastic water column model coupling with the graph-theoretic approach for simulating transients in water distribution networks. So far, these developed 1D numerical models and simulation methods have been widely used in practical applications such as system design and evaluation as well as theoretical development such as transient-based pipeline assessment (Duan et al., 2020).

Recently, some meshless numerical methods, such as the lattice Boltzmann method (LBM) and smoothed particle hydrodynamics (SPH), are also applied to water hammer analysis. The LBM method was first used for water hammer simulation by Cheng et al. (1998). Budinski (2016) utilised the LBM method with the adaptive grid to simulate 1D water hammer and gave detailed complex boundary conditions (such as pump, surge tank, reservoir, and valve, etc.) in WDSs, but didn’t consider the unsteady friction model for water hammer pipe flow. Later, Louati et al. (2019) modelled the two-dimensional (2D) water hammer flows using the multiple relaxation times LBM. Zeidan and Ostfeld (2022) developed a friction model for Lagrangian methods in transient flows and validated it in a looped water system, but the study didn’t consider some hydraulic components like pumps. The SPH approach was first proposed for astrophysical simulations by Lucy (1977) and Monaghan and Lattanzio (1985) and has become increasingly popular in simulating fluid flow due to its flexibility in tracing moving boundaries and adaptability to complicated geometries (Liu & Liu, 2003). Thus, the SPH method is widely used in complex flow simulations, such as flooding (Albano et al., 2016; Vacondio et al., 2013), multi-phase flow(De Padova et al., 2022; Dong et al., 2022; Meister & Rauch, 2016), and fluid-structure interactions (Basser et al., 2019; Lin et al., 2015; Ming et al., 2018). In addition, the SPH method has been applied in simulating channel and pipe flows with open boundaries (Chang & Chang, 2021; Nomeritae et al., 2018; Sun et al., 2022). The first application of the SPH on water hammer simulation was considered by Hou et al. (2011). Later, Pan et al. (2022) further considered the unsteady friction model in the SPH-based water hammer model, but it was only suitable for the simple pipeline. Compared with the traditional grid-based numerical methods, the SPH method showed an advantage in capturing the discrete process in transient pipe flow systems, such as column separation, rapid filling, and slug impact, because of its meshless feature (Hou et al., 2012). However, the majority of existing SPH-based water hammer models have been only developed for single or simple pipeline systems, and it is still a challenging task for complex pipeline systems such as WDSs, which has greatly limited its applications in practical engineering. And the main difficulties and challenges focus on the boundary treatments of various hydraulic components. Therefore, it is valuable and significant to build appropriate boundary treatment models and further investigate the application of the SPH method in WDSs.

Compared with the previous studies, this paper proposes an improved SPH model for Water Hammer simulation in WDSs (SPH-WHWDS), and the main contributions are as follows: (1) developing a more systematic and complete SPH model for water hammer problems; (2) extending the boundary conditions of the SPH model to complex pipeline systems, such as N-way nodal junctions, valves, reservoirs, surge tanks, and pumps, which greatly broadens the model’s applicability and is capable to simulate transient flows in practical WDSs; and (3) exploring the effects of various parameters in the SPH-WHWDS model, such as particle spacing and smoothing length, compared with experimental data.

The rest of this paper is structured as follows. In section 2, the governing equations of the 1D transient flow are presented, and the SPH formulation and computational domain are introduced. Besides, a detailed introduction to the solution of boundary parameters for various WDS components is provided in this section. In section 3, two benchmark cases, including the unsteady pipe flow experiment and a complex WDS, are investigated to verify and validate the proposed SPH-WHWDS model and the protection of the air chamber against water hammer is analysed. Finally, section 4 summarises the results and conclusions.

2. Models and methods

2.1. Governing equations

In general, 1D transient flow (water hammer) in pipes can be expressed by a set of two partial differential equations, namely the continuity and momentum equation (Wylie et al., 1993): (1) $\begin{array}{rl}\frac{\mathrm{\partial }H}{\mathrm{\partial }t}+v\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{a^2}{g}\frac{\mathrm{\partial }v}{\mathrm{\partial }x}& =0\end{array}$(1) (2) $\begin{array}{rl}\frac{\mathrm{\partial }v}{\mathrm{\partial }t}+v\frac{\mathrm{\partial }v}{\mathrm{\partial }x}+g\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+g\left(S_f-S_0\right)& =0\end{array}$(2) where $H$ is the pressure head, $v$ is the mean flow velocity, $\rho$ is the fluid density, $g$ is gravitational acceleration, $S_0$ is the slope of the pipe, and $S_f$ is the friction slope. Wave speed $a$ is calculated as: (3) $a=\sqrt{\frac{K/\rho }{1+\varphi \frac{K}{E}\frac{D}{e}}}$(3) where $K$ is the bulk modulus, $D$ is the inner pipe diameter, $E$ is Young’s modulus, $e$ is the pipe wall thickness, and the parameter $\varphi$ is determined by the pipe anchors.

Equations (1)–(2) can alternatively be written in the Lagrange form: (4) $\begin{array}{rl}\frac{DH}{Dt}& =-\frac{a^2}{g}\frac{\mathrm{\partial }v}{\mathrm{\partial }x}\end{array}$(4) (5) $\begin{array}{rl}\frac{Dv}{Dt}& =-g\frac{\mathrm{\partial }H}{\mathrm{\partial }x}-g\left(S_f-S_0\right)\end{array}$(5) where $D/Dt$ is the total time derivative $\left(\frac{D}{Dt}=\frac{\mathrm{\partial }}{\mathrm{\partial }t}+v\frac{\mathrm{\partial }}{\mathrm{\partial }x}\right)$.

2.2. The SPH method

The formulation and principle of applying the SPH method for simulating water hammer flows are elaborated in this section.

2.2.1. SPH formulation

In the SPH method, the fluid is discretized into a set of particles. An interpolation procedure is used to approximate particle properties (e.g. velocity, pressure) over neighbouring particles within a domain determined by a weighting function (Liu & Liu, 2003). Any physical quantity of particle $i$ ($f_i$) can be computed by the integral approximation: (6) $f_i=\sum _j\frac{m_j}{{\rho }_j}{f}_{j}W\left(|x_i-x_j|,h_i\right)$(6) where $m_j$ is the mass of particle $j$, ${\rho }_j$ is the density of particle $j$, $x_i$ is the position of particle $i$, $h_i$ is the smoothing length of particle $i$, $W\left(|x_i-x_j|,h_i\right)$ is the kernel function of particle $i$. Herein, the cubic spline function (Monaghan & Lattanzio, 1985) is selected as the kernel function: (7) $W\left(R,h\right)=\frac{1}{h}\{\begin{array}{cc}\frac{2}{3}-R^2+\frac{1}{2}{R}^3,& 0\le R<1\\ \frac{1}{6}{\left(2-R\right)}^3,& 1\le R<2\\ 0,& R\ge 2\end{array}$(7) where $R$ is the relative distance between two particles, $R=|x-x^{\mathrm{\prime }}|/h$.

The following are two commonly used SPH formulas for the spatial derivative of $f_i$: (8) $\begin{array}{rl}(\mathrm{\nabla }\cdot f{)}_i& =\frac{1}{{\rho }_i}\sum _j{m}_j(f_j-f_i)\cdot \mathrm{\nabla }{W}_i\left(|x_i-x_j|,h_i\right)\end{array}$(8) (9) $\begin{array}{rl}(\mathrm{\nabla }\cdot f{)}_i& ={\rho }_i\sum _j{m}_j\left(\frac{f_j}{{{\rho }_j}^2}+\frac{f_i}{{{\rho }_i}^2}\right)\cdot \mathrm{\nabla }{W}_i\left(|x_i-x_j|,h_i\right)\end{array}$(9)

2.2.2. Discretization of the governing equation of water hammer

The discretized SPH form of Equations (4) to (5) can be approximated by Equations (8) and (9) as follows: (10) $\begin{array}{rl}\frac{D{H}_i}{Dt}& =-\frac{a^2}{g}\sum _j{m}_j(v_j-v_i)\cdot \widetilde{\mathrm{\nabla }{W}_i}\end{array}$(10) (11) $\begin{array}{rl}\frac{D{v}_i}{Dt}& =-{\rho }_{i}g\sum _j{m}_j\left(\frac{H_j}{{{\rho }_j}^2}+\frac{H_i}{{{\rho }_i}^2}+{\mathrm{\Pi }}_{ij}\right)\\ & \cdot \widetilde{\mathrm{\nabla }{W}_i}-g\left(S_{f,i}-S_{0,i}\right)\end{array}$(11) with (12) $\begin{array}{rl}{\mathrm{\Pi }}_{ij}& =\{\begin{array}{cc}{\displaystyle \frac{-\alpha a{\phi }_{ij}+\text{β}{\phi }_{ij}^2}{\rho }},& (v_i-v_j)(x_i-x_j)<0\\ 0,& (v_i-v_j)(x_i-x_j)\ge 0\end{array}\end{array}$(12) (13) $\begin{array}{rl}{\phi }_{ij}& =\frac{\overline{h_{ij}}(v_i-v_j)(x_i-x_j)}{{(x_i-x_j)}^2+\eta {\overline{h_{ij}}}^2}\end{array}$(13) where ${\mathrm{\Pi }}_{ij}$ is the artificial viscosity required for the numerical method to be stable; $\overline{h_{ij}}$ = $0.5(h_i+h_j)$ is the average smoothing length of particles $i$ and $j$; $\alpha$, $\text{β}$ and $\eta$ are constant coefficients; $\widetilde{\mathrm{\nabla }{W}_i}$ is the correction of the spatial derivative of kernel function as (Chen et al., 1999): (14) $\widetilde{\mathrm{\nabla }{W}_i}=\frac{\mathrm{\nabla }{W}_i\left(|x_i-x_j|,h_i\right)}{{\sum }_j{m}_j/{\rho }_j(x_j-x_i)\cdot \mathrm{\nabla }{W}_i\left(|x_i-x_j|,h_i\right)}$(14) In this study, the computational domain of each pipe is divided into three zones, i.e. upstream buffer zone, in-pipe fluid zone (between the inflow boundary and the outflow boundary), and downstream buffer zone, as shown in Figure 1. Additionally, the computational domain has three types of particles: upstream-downstream buffer particles, fluid particles, and mirror particles. The parameters of buffer particles are determined by upstream and downstream boundary conditions, while mirror particles are used to simulate solid wall boundaries, such as closed valves or dead-end pipelines. The upstream and downstream zones of the pipe are usually nodes, valves, pumps, or other hydraulic components. To achieve the connection between the pipe and the hydraulic components, the study first calculates the upstream and downstream boundary parameters at the next moment, such as velocity and pressure, then assigns them to the buffer particles, and finally applies the SPH method to update the parameters of the fluid particles in the pipe.

Figure 1. Sketch of the computational domain of pipe.

2.2.3. Time-marching scheme

To integrate the particle pressure and velocity in time, the conditionally stable Euler forward method with first-order accuracy is used (Hou et al., 2012): (15) $\begin{array}{rl}{H}_i^{n+1}& =H_i^n+{\left(\frac{D{H}_i}{Dt}\right)}^n\cdot \mathrm{\Delta }t\end{array}$(15) (16) $\begin{array}{rl}{v}_i^{n+1}& =v_i^n+{\left(\frac{D{v}_i}{Dt}\right)}^n\cdot \mathrm{\Delta }t\end{array}$(16) (17) $\begin{array}{rl}{x}_i^{n+1}& =x_i^n+\frac{v_i^{n+1}+v_i^n}{2}\cdot \mathrm{\Delta }t\end{array}$(17)

The time step ($\mathrm{\Delta }t$) must satisfy the Courant–Friedrichs–Levy (CFL) condition as: (18) $\mathrm{\Delta }t\le CFL\cdot \text{min}\left(\frac{\mathrm{\Delta }{x}_0}{|v_i|+a}\right)$(18) where $CFL$ is the Courant number ($CFL\le 1$), $\mathrm{\Delta }{x}_0$ is the initial particle spacing.

2.3. Boundary treatments for typical WDSs

WDSs typically have complex boundary conditions which include a variety of hydraulic components (nodal junctions, valves, surge tanks, pumps, etc.). Therefore, one of the important tasks for transient pipeline flow modelling in complex looped WDS is the proper treatment of boundary conditions. To ensure the accuracy of the SPH-WHWDS model, the values of the unknown variables (pressure and velocity) of buffer particles at in/out-flow boundary must be defined at the start of the next time step. In this section, the solution of boundary parameters of different hydraulic components is developed and introduced in detail.

In order to specify the proper in/out-flow boundary conditions, the classic MOC method is used. Based on the directions of characteristic lines, characteristic equations (Chaudhry, 1993) are divided into positive and negative forms, as shown in Equations (19) to (22): (19) $C^{+}:\frac{dH}{dt}+\frac{a}{g}\frac{dv}{dt}+a\left(S_f-S_0\right)=0$(19) along (20) $C^{+}:\frac{dx}{dt}=v+a$(20) and (21) $C^{-}:\frac{dH}{dt}-\frac{a}{g}\frac{dv}{dt}-a\left(S_f-S_0\right)=0$(21) along (22) $C^{-}:\frac{dx}{dt}=v-a$(22) In addition, the specified time interval method (Chaudhry, 1993; Sturm, 2010) is used to calculate the unknown variables at in/out-flow boundaries along the positive and negative characteristic lines for each time step, as shown in Figure 2. The pressure and velocity at points L and R are determined by linear interpolation, one can refer to Federico et al. (2012), Chang and Chang (2013), and Aristodemo et al. (2015) for more details about this method.

Figure 2. Schematic of the specified time interval method.

2.3.1. N-way nodal junction

The nodal junction is the most common structure in WDSs. As shown in Figure 3, an N-way nodal junction boundary consists of m₁ inflowing pipes and m₂ outflowing pipes, with N = m₁+ m₂. 2N + 1 variables are unknown at the N-way nodal junction boundary, namely the pressure head at the node, the flow velocity and the pressure head at each pipe boundary, which are defined before the next time step. Based on the assumption of pressure head equality at the node and the mass conservation relationship, we can obtain Equations (23) to (24), and the equations of m₁ positive characteristic lines and m₂ negative characteristic lines are given by Equations (25) to (26). The node demands are based on the instantaneous pressure head and demand discharge coefficients, which allow the actual demands to change with the instantaneous local pressure, simulating more practical conditions (Gorev & Kodzhespirova, 2013; Xing & Sela, 2020). Therefore, a nonlinear system of equations is presented, involving Equations (23) to (27), which can be solved by the Newton-Raphson method (He, 2004). (23) $\begin{array}{rl}& H_1=H_2=\dots =H_N=H_d\end{array}$(23) (24) $\begin{array}{rl}& Q_1+Q_2+\dots +Q_N=Q_d\end{array}$(24) (25) $\begin{array}{rl}& Q_i=A_i{v}_i\end{array}$(25) (26) $\begin{array}{rl}& C^{+}:\{\begin{array}{rl}& H_1-H_{R,1}+\frac{a}{g}\left(v_1-v_{R,1}\right)\\ & +a\left(S_{f,1}-S_{0,1}\right)\mathrm{\Delta }t=0\\ & \cdots \\ & H_{m_1}-H_{R,m_1}+\frac{a}{g}\left(v_{m_1}-v_{R,m_1}\right)\\ & +a\left(S_{f,m_1}-S_{0,m_1}\right)\mathrm{\Delta }t=0\end{array}\end{array}$(26) (27) $\begin{array}{rl}& C^{-}:\{\begin{array}{rl}& H_{m_1+1}-H_{L,m_1+1}-\frac{a}{g}\left(v_{m_1+1}-v_{L,m_1+1}\right)\\ & -a\left(S_{f,m_1+1}-S_{0,m_1+1}\right)\mathrm{\Delta }t=0\\ & \cdots \\ & H_N-H_{L,N}-\frac{a}{g}\left(v_N-v_{L,N}\right)\\ & -a\left(S_{f,N}-S_{0,N}\right)\mathrm{\Delta }t=0\end{array}\end{array}$(27)

Figure 3. SPH Schematic of an N-way nodal junction boundary.

with (28) $Q_d=\{\begin{array}{cc}k\sqrt{H_d},& H_d>H_{demand}\\ 0,& H_d\le H_{demand}\end{array}$(28) where $A_i$ is the cross-sectional area of pipe $i$, $Q_d$ is the demand discharge of the nodal junction, $k$ is the demand discharge coefficient, $H_{demand}$ is the minimum service head, subscripts $L$ and $R$ denote the foot of the negative characteristic line and positive characteristic line, respectively.

2.3.2. Valve

A pipe valve (Figure 4) can regulate the flow of water while causing large energy loss. By applying the conservation equation of mass and energy, the following equation can be obtained: (29) $\begin{array}{rl}{Q}_{p1}& =Q_{p2}=Q_v\end{array}$(29) (30) $\begin{array}{rl}{H}_{p1}& =H_{p2}-{\xi }_v\frac{Q_v|Q_v|}{2g{A}_v^2}\end{array}$(30) where ${\xi }_v$ is the loss coefficient, and $A_v$ is the cross-sectional area of the valve.

Figure 4. SPH Schematic of a valve component: (a) inline valve; (b) end valve.

Two characteristic equations are also included: (31) $\begin{array}{rl}& H_{p1}-H_{R,p1}+\frac{a}{g}\left(v_{p1}-v_{R,p1}\right)\\ & +a\left(S_{f,p1}-S_{0,p1}\right)\mathrm{\Delta }t=0\end{array}$(31) (32) $\begin{array}{rl}& H_{p2}-H_{L,p1}-\frac{a}{g}\left(v_{p2}-v_{L,p2}\right)\\ & -a\left(S_{f,p2}-S_{0,p2}\right)\mathrm{\Delta }t=0\end{array}$(32) By solving Equations (29) to (32), we can obtain the unknown variables at the valve. However, when the valve is fully closed, it can be considered as a solid wall boundary, and the virtual boundary particle (VBP) method (Ferrari et al., 2009) is used to handle this condition. The mirror particles are placed outside of the end valve. The mirror particle and the corresponding actual particle are at the same distance from the valve. The pressure of the mirror particle is the same as the corresponding fluid particle, while the velocity is the opposite.

2.3.3. Reservoir

For a reservoir located at the beginning or end of the pipe (Figure 5), providing a constant pressure head, the parameters of buffer particles are determined by the reservoir and the particles closest to it. Using the mass and energy equation, the following equations are derived (Figure 5): (33) $\begin{array}{rl}{Q}_{buffer}& =Q_{p1}\end{array}$(33) (34) $\begin{array}{rl}{H}_R& =H_{buffer}\pm \left(1+{\xi }_{in}\right)\frac{Q_{p,1}|Q_{p,1}|}{2g{A}_p^2}\end{array}$(34) where ${\xi }_{in}$ is the head loss coefficient.

Figure 5. SPH Schematic of a reservoir.

2.3.4. Surge tank or air chamber

The surge tank is usually used to reduce and eliminate surges caused by the water hammer effect, including the open surge tank and the air chamber, as shown in Figure 6. During the transient simulation, the open surge tank connected to the pipeline stores or supplies water according to pressure changes in the system, thereby moderating pressure transients. To introduce the open surge tank as a boundary condition, the following set of equations is used: (35) $\begin{array}{rl}{Q}_{p1}-Q_{p2}& =Q_{tank}\end{array}$(35) (36) $\begin{array}{rl}{H}_{p1}& =H_{p2}=H_{tank}+{\xi }_t\frac{Q_{tank}|Q_{tank}|}{2g{A}_{pt}^2}\end{array}$(36) (37) $\begin{array}{rl}\frac{d{H}_{tank}}{dt}& =\frac{Q_{tank}}{A_{tank}}\end{array}$(37) where $Q_{tank}$ is the flow into the open surge tank, $H_{tank}$ is the water head in the tank, $A_{pt}$ is the cross-sectional area of the connection pipe, ${\xi }_t$ is the loss coefficient of the connection pipe, $A_{tank}$ is the cross-sectional area of the open surge tank. With two characteristic equations [Equations (31) to (32)], six unknowns ($Q_{p1},Q_{p2},Q_{tank},H_{p1},H_{p2},H_{tank}$) can be solved at each time step.

Figure 6. SPH Schematic of surge tank: (a) open surge tank; (b) air chamber.

For small WDSs, an air chamber is a more practical solution to relieve pressure transients. Compared to the open surge tank, the air chamber utilises air pressure to slow down the deceleration or the acceleration of water flow, thus reducing the size and cost of the device. To introduce the air chamber as a boundary condition, the following set of equations is used: (38) $\begin{array}{rl}{Q}_{p1}-Q_{p2}& =Q_{ac}\end{array}$(38) (39) $\begin{array}{rl}{H}_{p1}& =H_{p2}=H_{ac}+{\xi }_{ac}\frac{Q_{ac}|Q_{ac}|}{2g{A}_{pt}^2}\end{array}$(39) (40) $\begin{array}{rl}\frac{d{z}_{ac}}{dt}& =\frac{Q_{ac}}{A_{ac}}\end{array}$(40) (41) $\begin{array}{rl}\frac{d{V}_{ac}}{dt}& =Q_{ac}\end{array}$(41) (42) $\begin{array}{rl}{H}_{ac}& =z_{ac}+\frac{P_{ac}}{\rho g}\end{array}$(42) (43) $\begin{array}{rl}\left(P_{ac}+P_{atm}\right)\cdot V_{ac}^c& =constant\end{array}$(43) where $Q_{ac}$ is the flow into the air chamber, $H_{ac}$ is the total head, $z_{ac}$ is the water level, $A_{ac}$ is the cross-sectional area of the air chamber, $V_{ac}$ is the volume of the air in the air chamber, $P_{ac}$ is the air pressure, $P_{atm}$ is the standard atmospheric pressure, $c$ is the polytropic index. With two characteristic equations [Equations (31) to (32)], nine unknowns ($Q_{p1},Q_{p2},Q_{ac},H_{p1},H_{p2},H_{ac},z_{ac},V_{ac},P_{ac}$) can be solved at each time step.

2.3.5. Pump

Pump is an important power device in the WDSs, and pressure fluctuations occur when the pump is start-up or shut-off. In this study, the simulation of pump start-up and shut-off is achieved by changing the pump rotation speed over time. The boundary condition at the pump can be solved by using one characteristic equation, the affinity law [Equation (46)], and the pump characteristic curve at a given rotation speed [Equation (45)]. Additionally, the study assumes that the check valve prevents the pump from reversing, and the pump shutting down due to power failure has not been considered yet. To introduce the pump as a boundary condition, the following set of equations is utilised: (44) $\begin{array}{rl}{N}_{pump}& =f\left(t\right)\end{array}$(44) (45) $\begin{array}{rl}{H}_{pump}& =C_1+C_2{Q}_{pump}+C_3{Q}_{pump}^2\end{array}$(45) (46) $\begin{array}{rl}\frac{H_{pump,1}}{H_{pump,2}}& ={\left(\frac{N_{pump,1}}{N_{pump,2}}\right)}^2\end{array}$(46) where $N_{pump}$ is the pump rotation speed, and $C_1,C_2{C}_3$ are the coefficients of the pump characteristic curves at a given rotation speed. With one characteristic equation [Equations (31) or (32)], two unknowns ($H_{pump}$, $Q_{pump}$) can be solved at each time step.

It is noted that only the typical elements in practical WDSs are developed and listed in this study, while for other specific or special elements, a similar treatment way can be applied as long as their hydraulic characteristics are known for the SPH formulation.

3. Results and discussion

In this section, the performance of the proposed SPH-WHWDS model is evaluated by two test cases: (1) the unsteady single-pipe flow experiment by Bergant et al. (2001), with two reservoirs, one pipe, and one valve; and (2) a complex WDS, with two reservoirs, two pumps, nineteen pipes, a valve, and an air chamber. In the first case, the proposed model is validated with experimental data, and the Brunone friction model (Brunone et al., 1991) is utilised to calculate the unsteady friction. In the second case, we choose valve closure and pump shut-off, which are the most common reasons for inducing transient flows in WDSs, for assessing the proposed numerical model on more practical conditions and comparing the results calculated by the proposed model with the classic MOC scheme. Additionally, the effect of the air chamber on water hammer protection is demonstrated in the second test.

3.1. Experiment test

In this test case, a simple reservoir-pipe-valve experiment by Bergant et al. (2001) is chosen to evaluate the effectiveness of the proposed model in unsteady pipe flow, and the transient event is generated by the instantaneous closing of the valve.

The Brunone friction model (Brunone et al., 1991), which is dependent on instantaneous mean flow velocity $v$, the local acceleration $\mathrm{\partial }v/\mathrm{\partial }t$, and the convective acceleration $\mathrm{\partial }v/\mathrm{\partial }x$, is used to calculate the friction slope $S_f$ in Eq.(2), and can be expressed as follows: (47) $\begin{array}{rl}{S}_f& =\lambda \frac{v|v|}{2gD}\end{array}$(47) (48) $\begin{array}{rl}\lambda & ={\lambda }_q+\frac{k_{b}D}{v|v|}\left[\frac{\mathrm{\partial }v}{\mathrm{\partial }t}-a\cdot sign\left(v\right)\left|\frac{\mathrm{\partial }v}{\mathrm{\partial }x}\right|\right]\end{array}$(48) where ${\lambda }_q$ is the steady friction factor, $sign\left(v\right)$ is the sign of velocity ($+1forv\ge 0or-1forv<0$), $k_b$ is Brunone’s friction coefficient, which is empirically estimated by trial and error.

The local acceleration is solved by the explicit scheme between the $n-1$ and $n$ time step for particle $i$: (49) $\frac{\mathrm{\partial }{v}_i}{\mathrm{\partial }t}=\frac{v_i^n-v_i^{n-1}}{\mathrm{\Delta }t}$(49) The convective acceleration is calculated by the implicit scheme, which can be derived from Equation (8) as follows: (50) $\frac{\mathrm{\partial }{v}_i}{\mathrm{\partial }x}=\frac{1}{{\rho }_i}\sum _j{m}_j(v_j-v_i)\cdot \mathrm{\nabla }{W}_i\left(|x_i-x_j|,h_i\right)$(50) As shown in Figure 7, the experimental apparatus consists of two pressurised tanks connected by a straight copper pipe measuring 37.23 m in length, with an internal diameter of 22.1 mm and a wall thickness of 1.63 mm, and a downstream instantly closed valve. A computerised pressure control system regulates each of the tanks’ pressure to a constant value. Other parameters of the experimental apparatus are as follows: wave speed $a$ = 1319 m/s; valve closure time $t_v$ =  0.009 s; pipe slope $s_0$ =  5.45%. Based on Bergant’s experiment (Bergant et al., 2001), three water hammer cases with different initial velocities (0.1, 0.2, 0.3 m/s) are chosen to verify the performance of the proposed model. The pressure fluctuations at the valve ($H_v$) and the middle of the pipe ($H_m$) are recorded during the experiment.

Figure 7. Schematic diagram of Bergant’s experimental test system.

In the simulations, the spacing between the particles is 0.15 m ($\mathrm{\Delta }x=0.15m$), which means 249 particles are evenly distributed along the pipe, the time step is $9\times 1{\text{0}}^{-\text{5}}\text{s}$ (CFL = 0.8), and the smoothing length $h=1.3\mathrm{\Delta }x$ is adopted here. Figures 8–10 depict the results of experimental data and numerical simulation for the three selected water hammer cases at the midpoint and valve, respectively. The results demonstrate the proposed SPH-WHWDS model can obtain almost the same pressure oscillations as the experimental data, which implies the method is reliable in unsteady pipe flow simulation. With the evolution of transient events, the pressure oscillation period of the SPH-WHWDS model is slightly shorter than the experimental data, which is similar to the numerical simulation results of Bergant et al. (2001), and this difference becomes more significant as the initial velocity increases. The reason may be that a constant wave speed is used in the simulation, whereas the actual wave speed may slightly decrease in that test system.

Figure 8. Comparison of numerical and experimental measured pressure heads traces (0.1 m/s): (a) at the midpoint ($H_m$) and (b) the valve ($H_v$).

Figure 9. Comparison of numerical and experimental measured pressure heads traces (0.2 m/s): (a) at the midpoint ($H_m$) and (b) the valve ($H_v$).

Figure 10. Comparison of numerical and experimental measured pressure heads traces (0.3 m/s): (a) at the midpoint ($H_m$) and (b) the valve ($H_v$).

In order to test the influence of particle spacing in the SPH-WHWDS model, we take the case with initial velocity $v_0=0.1\text{m/s}$ as an example and adopt five different particle spacing configurations ($\mathrm{\Delta }x$ = 0.05, 0.1, 0.15, 0.25, 0.5 m). Hence, the corresponding numbers of particles are 745, 373, 249, 149, and 75, while for the corresponding time steps, $\mathrm{\Delta }t$  = $3\times 1{\text{0}}^{-\text{5}}\text{s}$, $6\times 1{\text{0}}^{-\text{5}}\text{s}$, $9\times 1{\text{0}}^{-\text{5}}\text{s}$, $\text{1}.5\times 1{\text{0}}^{-\text{4}}\text{s}$, and $3\times 1{\text{0}}^{-\text{4}}\text{s}$ are utilised. A comparison of the pressure heads at the midpoint with five different particle spacing configurations is shown in Figures 11 and 12. Except for $\mathrm{\Delta }x$ = 0.5 m, it is clear that the other four curves can reproduce the historical traces of pressure fluctuations and only subtle differences can be seen in the later periods of the simulation. The comparison of results also indicates that the amplitude of pressure fluctuation decreases with the increase of particle spacing, especially after 0.4 s. As shown in Figure 13, with the decrease of particle spacings, the highest and lowest pressure heads near the upstream tank are more accurate, while the results at the middle pipe and the downstream valve are similar. Additionally, although the computation time reduces with fewer particles involving the computational domain, the model causes more pressure dissipation and the accuracy decreases gradually, which can be seen in the magnified view at 0.24s∼0.29s in Figures 11 and 12. Therefore, the appropriate particle spacing should be selected during the simulation to balance the efficiency and accuracy of the developed SPH-WHWDS model. Specifically, in this benchmark test case, particle spacing with $\mathrm{\Delta }x$ = 0.1m∼0.15 m is proven to be more appropriate.

Figure 11. Comparison of pressure head traces at the midpoint with different particle spacings (0.1 m/s).

Figure 12. Comparison of pressure head traces at the valve with different particle spacings (0.1 m/s).

Figure 13. Pressure head envelope curves with different particle spacings.

According to the conclusion of Hou et al. (2012), when the smoothing length coefficient is less than 1, it will lead to evident dispersion errors, so only cases with smoothing length coefficients greater than 1 are considered in the study. Therefore, we also take the situation with initial velocity $v_0=0.1\text{m/s}$ and particle spacing $\mathrm{\Delta }x$= 0.15 m for demonstration, applying five different smoothing length configurations ($h$=$1.0\mathrm{\Delta }x$, $1.3\mathrm{\Delta }x$, $1.5\mathrm{\Delta }x$, $2.0\mathrm{\Delta }x$, $2.5\mathrm{\Delta }x$). As shown in Figures 14 and 15, except for $h=2.5\mathrm{\Delta }x$, the remaining four curves almost coincide with each other, and the differences can only be seen in the later simulation stages. Based on the magnified plot at 0.24s∼0.29 s in Figures 14 and 15, it can be observed that the model’s prediction of the pressure amplitude at the midpoint decreases gradually with the increase of the smoothing length $h$. While the smoothing length $h$ is $1.0\mathrm{\Delta }x$, the predicted peak value of pressure fluctuation is slightly smaller at the valve, and an oscillation error occurs around 0.1 s, and this error is noticeable in pressure head envelope curves, as shown in Figure 16. Hence, for the unsteady pipe flow modelling investigated herein, the selection of smoothing length coefficients between 1.3 and 1.5 is proven to be more appropriate.

Figure 14. Comparison of pressure head traces at the midpoint with different smoothing lengths (0.1 m/s).

Figure 15. Comparison of pressure head traces at the valve with different smoothing lengths (0.1 m/s).

Figure 16. Pressure head envelope curves with different smoothing lengths.

3.2. A complex water distribution system

In the second case, a more practical and complex WDS, as shown in Figure 17, is utilised to evaluate the proposed SPH-WHWDS model. The system is composed of two reservoirs, two pumps, nineteen pipes, seventeen nodes, a valve, and an air chamber. All the pipes in this test case have some similar characteristics: Darcy–Weisbach friction factor = 0.02; bed slope = 0; wave speed = 1000 m/s. Other relevant parameters for the system are depicted in Figure 17 and Table 1. The initial steady state conditions, such as the flow of the pipelines and the pressure heads of the nodes, are calculated by EPANET. After establishing the initial state of the WDS, the transient events are generated by valve closure and pump shut-off respectively. To validate the protection of the air chamber against the water hammer effect, the transient events with and without an air chamber are compared, resulting in four different cases in total.

Figure 17. Schematic layout of a practical water distribution system.

Table 1. Parameters of pipeline in Case 2.

Pipeline	Diameter(mm)
p11	80
p7, p8, p14, p15, p16, p17	100
p6, p9, p10, p12, p13, p18, p19	150
p1, p2, p5	200
p3, p4	250

Before the simulation, the spacing or number of particles in each pipe needs to be specified. Due to the different lengths of pipes in the looped multi-pipe system, the particle spacing of each pipe can be the same or adjusted adaptively, and the time step is determined by the smallest particle spacing. In this study, we choose the same particle spacing ($\mathrm{\Delta }x=$ 0.8 m) for each pipe, for the longest pipes p₁ and p₂, the number of particles is 188; while for the shortest pipe p₁₆, the number of particles is 45. For testing and verification of the developed SPH-WHWDS model, a publishing transient simulation package TSNet (Xing & Sela, 2020) is used in this case, which adopts the classic MOC scheme as the solution technique. The grid size set in TSNet is the same as the particle spacing in the SPH-WHWDS model for a fair comparison.

In addition, to investigate the numerical accuracy, the $L_2$ relative error norm based on the variable $\phi$ as shown in Equation (45) is therefore introduced. (51) $L_2\left(\phi \right)=\sqrt{\frac{\sum {\left({\phi }_i^{SPH-WHWDS}-{\phi }_i^{MOC}\right)}^2}{\sum {\left({\phi }_i^{MOC}\right)}^2}}$(51) where ${\phi }_i^{SPH-WHWDS}$ and ${\phi }_i^{MOC}$ are the simulated results by the SPH-WHWDS model and the MOC scheme, respectively.

3.2.1. Scenario 1: water hammer by valve closure

First, the transient regime of the system is induced by the instantaneous closing of the valve after node N₁₇ in Figure 17. This triggers the water hammer effect, which is followed by intense pressure and flow changes throughout the system. An air chamber is designed and installed at node N₄ to protect the upstream pumps from severe damage caused by drastic pressure fluctuations during transient events. The basic parameters of the air chamber are as follows: the cross-sectional area $A_{ac}$ = 1 m², height $H$ = 3 m, initial water level $z_{ac0}$ = 1.5 m, initial air volume $V_{ac0}$ = 1.5 m³, and the loss coefficient of the connection pipe is ignored. For the time step and smoothing length, value $\mathrm{\Delta }t$ = $\text{6}.4\times 1{\text{0}}^{-\text{4}}\text{s}$ (CFL = 0.8) and $h$ = $\text{1}\text{.5}\mathrm{\Delta }x$ are adopted respectively. To compare results from the SPH-WHWDS and MOC method, the upstream node N₃ and the downstream node N₅ near the air chamber are selected. Comparisons of the simulated pressure heads are shown in Figures 18 and 19. Furthermore, the highest and lowest pressure head along the main pipe (p₁, p₃, p₄, p₉, p₁₂, p₁₃, p₁₉) are recorded and compared in Figure 20.

Figure 18. Comparison between the SPH-WHWDS and MOC results of pressure heads at node N₃ in scenario 1: (a) without an air chamber; (b) with an air chamber.

Figure 19. Comparison between the SPH-WHWDS and MOC results of pressure heads at node N₅ in scenario 1: (a) without an air chamber; (b) with an air chamber.

Figure 20. The SPH-WHWDS and MOC results of pressure head envelope curves along the main pipes in scenario 1: (a) without an air chamber; (b) with an air chamber.

As illustrated in Figures 18(a) and 19(a), without an air chamber in the system, the pressure heads of nodes N₃ and N₅ reach the maximum value of about 100 m at t = 0.6 s, and the minimum value of 5 m at t = 0.8 s. Subsequently, the energy gradually dissipates and the pressure tends to be stable. While with an air chamber equipped in the system, as shown in Figures 18(b) and 19(b), the pressure head at node N₅ still fluctuates violently, but the pressure at node N₃ fluctuates very slightly, remaining around 48 m. In Figure 20, it is evident that the air chamber at node N₄ can moderate and prevent the downstream violent pressure oscillation from propagating upstream when the downstream valve is suddenly closed, thus ensuring the normal operation of upstream pumps. According to Equation (51), without an air chamber, the $L_2\left(H_{N3}\right)$ and $L_2\left(H_{N5}\right)$ are calculated as 5.21% and 4.62%, respectively; while with an air chamber, the $L_2\left(H_{N3}\right)$ and $L_2\left(H_{N5}\right)$ are calculated as 0.83% and 5.01% respectively. Through the comparison of the results for different cases, it is clear that the SPH-WHWDS results show good agreement with the MOC solutions in valve closure events. In addition, at the crest and trough, the SPH-WHWDS model captures pressure fluctuations more smoothly, while the MOC is more detailed. Although the simulation results of node pressure heads are smoother, the proposed SPH-WHWDS model is still reliable in predicting the cycles and peaks of pressure fluctuations in complex WDSs.

3.2.2. Scenario 2: water hammer by pump shut-off

The second water hammer event is induced by a controlled pump shut-off at Pump 2 in Figure 17. During the simulation, the pump rotational speed of Pump 2 reduced linearly to zero within 5 s, while Pump 1 operated normally. Due to the check valve behind the pump, the reverse flow to the pump was not considered in this case. The time step and smoothing length are the same as in scenario 1, i.e. $\mathrm{\Delta }t$ = $\text{6}.4\times 1{\text{0}}^{-\text{4}}\text{s}$ and $h$ = $\text{1}\text{.5}\mathrm{\Delta }x$. The upstream and downstream nodes (N₃ and N₅) of the air chamber are also selected for analysis. Comparisons of the obtained pressure heads for these two nodes are presented in Figures 21 and 22, respectively. And the highest and lowest pressure heads along the main pipe are shown in Figure 23.

Figure 21. Comparison between the SPH-WHWDS and MOC results of pressure heads at node N₃ in scenario 2: (a) without an air chamber; (b) with an air chamber.

Figure 22. Comparison between the SPH-WHWDS and MOC results of pressure heads at node N₅ in scenario 2: (a) without an air chamber; (b) with an air chamber.

Figure 23. The SPH-WHWDS and MOC results of pressure head envelope curves along the main pipes in scenario 2: (a) without an air chamber; (b) with an air chamber.

Without an air chamber in the WDS, as shown in Figures 21(a) and 22(a), the pressure heads of nodes N₃ and N₅ drop to about 27 m within 3s, and the amplitude of pressure drops to 20 m approximately. Subsequently, due to the operation of pump 1, the pressure heads of nodes N₃ and N₅ gradually increase, and finally become stabilised at about 40 m. The results show that when the pump is shutting down, a significant transient can be produced in the WDS. Therefore, it is important to assess the impact of pump operation on piping and propose appropriate procedures to guide and protect pump operation, such as the installation of water hammer protection devices. As a result, when an air chamber is installed at node N₄, as shown in Figure 21(b) and 22(b), the pressure fluctuations of nodes N₃ and N₅ are greatly changed compared to former situations without an air chamber. Specifically, node N₃ is located between the pump and the air chamber, and the pressure heads reduce from 47.43 m to approximately 40 m within 2s (i.e. about 7 m reduction only) and then begin to fluctuate, with a range between 35 and 50 m. While node N₅ is located downstream of the air chamber, its pressure drops gradually from 45.78 m to 37 m, but without obvious pressure fluctuation. As shown in Figure 23, with an air chamber in the system, the amplitude of pressure fluctuation along the main pipe reduces from 20 m to 10 m in the downstream pipeline. In conclusion, when the upstream pump is shutting off, the air chamber can prevent the severe pressure fluctuation from spreading downstream, effectively protecting the downstream WDS.

Besides, without an air chamber, the $L_2\left(H_{N3}\right)$ and $L_2\left(H_{N5}\right)$ are calculated as 2.73% and 3.72% through Eq.(51), respectively; while with an air chamber, the $L_2\left(H_{N3}\right)$ and $L_2\left(H_{N5}\right)$ are calculated as 1.79% and 1.06%, respectively. Similar to scenario 1, although the SPH-WHWDS model may slightly underestimate the pressure head results in comparison with the classical MOC scheme, the overall simulation results demonstrate that these two models are almost consistent in terms of both the phase shift and amplitude attenuation of pressure heads. From this perspective, it is evident that the SPH-WHWDS model is capable of simulating water hammer events in complex WDSs, which can provide an alternative tool for practical applications. Based on the developed SPH-WHWDS model in this study, it is worthy of further exploring and visualising the complex flow phenomena and dynamic processes such as column separation and slug flow in practical WDSs in future work.

4. Conclusions

This study aims to extend the SPH-based water hammer model to more complex WDSs. Compared with the previous SPH model, which could only simulate a single pipe with reservoirs or valves, the developed SPH-WHWDS model in this study combines the MOC to solve the complex boundary conditions for simulating water hammer problems in complex WDSs. A detailed method for solving various and typical boundary conditions of pipeline and hydraulic components is formulated for the developed SPH-WHWDS framework. Two benchmark test cases, including a simple reservoir-pipe-valve experiment and a looped complex WDS, are investigated to validate and verify the performance of the developed SPH-WHWDS model. The key results and findings can be summarised as follows:

1. The developed SPH-WHWDS model in this study overcomes the limitations of the existing SPH-based water hammer model, which further extends the applicability of the advantageous SPH method to complex WDSs, thereby improving the foundation for future research on complex water hammer phenomena in WDSs.

2. The proposed SPH-WHWDS model is capable to reproduce the experiment results in unsteady pipe flows. Besides, with the increase of the smoothing length $h$ and particle spacing $\mathrm{\Delta }x$, the amplitude of pressure fluctuation decreases gradually in later wave periods. For Bergant’s unsteady pipe flow experiment investigated herein, to ensure numerical accuracy and computational efficiency, the particle spacing is recommended to be between 0.1 and 0.15 m, and the smoothing length factor is suggested to be between 1.3 and 1.5.

3. The proposed SPH-WHWDS model can also accurately simulate the water hammer events caused by both valve closing and pump shut-off in WDS. The boundary conditions of hydraulic components in the system are solved accurately and reliably. In addition, the numerical simulation results reveal that the air chamber can moderate the pressure fluctuations in the downstream pipelines, indicating the proposed model can evaluate the performance of water hammer protection devices in the WDS, which is helpful to guide pipeline design and operation.

Finally, the application results in this study demonstrate the capability of the developed SPH-WHWDS model to reproduce both the results from the experimental data and the MOC scheme. Such extension of the SPH model can provide an alternative tool and approach to understanding and capturing the complex flow phenomena in practical WDSs such as column separation, slug flows, and rapid filling, which deserves further investigations in future work.

Notations

$a$	=	wave speed
$A$	=	cross-sectional area
$c$	=	the polytropic index
CFL	=	the Courant number
$D$	=	pipe inner diameter
$g$	=	gravitational acceleration
$h_i$	=	smoothing length of particle i
$H$	=	pressure head
$\overline{h_{ij}}$	=	average smoothing length of particles i and j
$H_{demand}$	=	minimum service head
$k$	=	demand discharge coefficient
$k_b$	=	Brunone’s friction coefficient
$m_i$	=	mass of particle i
$N_{pump}$	=	pump rotation speed
$P_{ac}$	=	air pressure
$P_{atm}$	=	standard atmospheric pressure
$Q$	=	flow discharge
$S_0$	=	pipe slope
$S_f$	=	friction slope
$t$	=	time
$v$	=	mean flow velocity
$V_{ac}$	=	volume of the air in air chamber
$W$	=	kernel function
$\mathrm{\nabla }W$	=	the spatial derivative of kernel function
$\widetilde{\mathrm{\nabla }W}$	=	the correction of the spatial derivative of kernel function
$z_{ac}$	=	water level in air chamber
${\lambda }_q$	=	steady friction factor
$\xi$	=	loss coefficient
$\rho$	=	density
${\mathrm{\Pi }}_{ij}$	=	artificial viscosity

Abbreviations

SPH	=	smoothed particle hydrodynamics
WDSs	=	water distribution systems
MOC	=	method of characteristics
FDM	=	finite-difference method
FVM	=	finite-volume method

Acknowledgments

This work was financially supported by National Natural Science Foundation of China (Grant Nos. 51978493 and 51808396), the Shanghai Pujiang Program (Grant Nos.20PJ1417500) and the Hong Kong Research Grants Council (RGC) under project No.15200719. We also thank the anonymous reviewers and editors for their comments and suggestions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China [grant number 51808396, 51978493]; the Shanghai Pujiang Program [grant number 20PJ1417500]; the Hong Kong Research Grants Council (RGC) [grant number 15200719].