**3. Numerical Solution of Transient Flow Equations**

In order to numerically solve the set of Equations (2) and (10), the MacCormack predictor–corrector finite difference explicit scheme was used [20,21]. This method is widely used in computational fluid dynamics [22], but it was only recently adapted to simulate water hammer in elastic pipes by [23] and also successfully used in [24]. In this

model, the time step is no longer subjected to the spatial step, and thus, it is more convenient for meshing variable property series pipes than the classic method of characteristics is. The numerical grid was divided into sections that represent individual pipes with given properties. This allows for assigning different values of input data corresponding to each pipe (i.e., inner diameter, creep compliance, retardation time, number of Kelvin–Voigt elements). The pressure wave travels in the pipeline through the junction node connecting individual pipes. To illustrate the discretisation scheme and solving process, a numerical mesh is shown in Figure 2.

**Figure 2.** Numerical grid for the MacCormack explicit scheme.

Furthermore, subscripts denoting different types of numerical nodes are as in Figure 2.


For a pipe connected to a constant head reservoir, the boundary condition for the first node is as follows:

$$h\_1^{t + \Delta t} = h\_0 \tag{12}$$

$$v\_1^{t+\Delta t} = -\frac{c^2 v\_2^{t+\Delta t}}{\left(h\_2^{t+\Delta t} - h\_1^{t+\Delta t}\right)g - c^2} \tag{13}$$

where Δ*t* is the time step and *h*<sup>0</sup> is the initial head.


For a closed valve, the outlet condition is defined as follows:

$$v\_s^{t+\Delta t} = 0\tag{14}$$

Using the backward difference, the outlet boundary condition can be calculated:

$$h\_s^{t+\Delta t} = h\_s^t + \Delta t \left(\frac{c^2}{g} \frac{v\_{s-1}^{t+\Delta t}}{\Delta x}\right) \tag{15}$$


The derivation of the similar numerical discrete format for the internal nodes of the elastic pipeline system is in [12]. Here, Equations (1) and (10) are presented in basic time-marching format only:

$$\left(\frac{\partial h}{\partial t}\right)\_{n}^{t} = -\upsilon\_{n}\frac{h\_{n+1}^{t} - h\_{n}^{t}}{\Delta x} - \frac{c^{2}}{g}\frac{\upsilon\_{n+1}^{t} - \upsilon\_{n}^{t}}{\Delta x} - \frac{2c^{2}}{g}\frac{\partial \varepsilon\_{r}^{t-\Delta t}}{\partial t} \tag{16}$$

$$\left(\frac{\partial v}{\partial t}\right)\_{n}^{t} = -g\frac{h\_{n+1}^{t} - h\_{n}^{t}}{\Delta x} - v\_{n}^{t}\frac{v\_{n+1}^{t} - v\_{n}^{t}}{\Delta x} - \frac{f v\_{n}^{t}|v\_{n}^{t}|}{2g} \tag{17}$$

The retarded strain in Equation (10) is calculated as a sum of each Kelvin–Voigt element:

$$\frac{\partial \varepsilon\_r^{t-\Delta t}}{\partial t} = \sum\_{k=1}^{N} \frac{\partial \varepsilon\_{rk}^{t-\Delta t}}{\partial t} = \sum\_{k=1}^{N} \left[ \frac{aD}{\Delta s} \frac{f\_k}{\tau\_k} \rho g \left( h^{t-\Delta t} - h\_0 \right) - \frac{\varepsilon\_{rk}^{t-\Delta t}}{\tau\_k} \right] \tag{18}$$

Discrete approximation of each retarded strain is equal to the following:

$$\varepsilon\_{rk}^{t-\Delta t} = f\_k F^{t-\Delta t} - f\_k e^{-\frac{\Delta t}{\tau\_k}} F^{t-2\Delta t} - f\_k \tau\_k \left( 1 - e^{-\frac{\Delta t}{\tau\_k}} \right) \frac{F^{t-\Delta t} - F^{t-2\Delta t}}{\Delta t} + e^{-\frac{\Delta t}{\tau\_k}} \varepsilon\_{rk}^{t-2\Delta t} \tag{19}$$

with the function *F* at time *t*-Δ*t* described as follows:

$$F^{t-\Delta t} = \frac{aD}{2s} \frac{f\_k}{\tau\_k} \rho g \left( h^{t-\Delta t} - h\_0 \right) \tag{20}$$

It should be noted that in this numerical model, the constant value of the pressure wave velocity describes the velocity of the transient wave travelling through the entire pipeline system.
