## Abstract

Despite the applications of the Method of Characteristics (MOC) for analyzing the unsteady flows, using this method in networks with variable elevations still has many challenges. In this paper, by developing modified correlations as a computer code, the possibility of analyzing inclined pipelines has been evaluated. For validation and calibration, the results of MOC were compared with the results of EPANET software as well as experimental data. To extract experimental data, the water network of Babol Noshirvani University of Technology (NIT) with a constant head of 7 m three loops, and four inclined branches were employed. While evaluating the capabilities of the developed computer code, the results show that for all pipes, as the number of pressure fluctuations in a specific period increases, the intensity of the pressure fluctuations decreases, and the damping speed increases as well. Moreover, in inclined pipes, unlike noninclined pipes, the intensity of pressure fluctuations will increase as the elevation increases and the cross-sectional distance from the transient event increases as well. The evaluation of the effect of space steps on the accuracy of the solution to the MOC shows that in the study network, considering 20 segments for each pipe, the fastest response time with an error of less than 1% is obtained.

## HIGHLIGHTS

The modified correlations for use in a developed computer code based on the method of characteristics are presented.

The possibility of analyzing inclined water networks has been evaluated.

Transient signal properties and pressure fluctuations in the inclined experimental network were discussed.

Selecting the appropriate space step increases the solution speed and achieves acceptable solution accuracy.

## INTRODUCTION

The presentation and developments of models that can precisely simulate the flow in the pipes and water distribution networks (WDNs) are significantly important (Kapelan 2002). Nowadays, different methods have been presented for hydraulic analysis of WDNs with various purposes, such as estimating the amount of leakage (Negharchi & Shafaghat 2020). To select a method or an optimized tool for analyzing the flow in WDNs, the flow state in terms of being steady or unsteady is considerably determining. Different tools in software packages have been developed that present accurate and reliable results for analyzing steady flows. EPANET and Watergems are some of these software packages. It is noteworthy that usually, in WDNs, considering the unsteady pressure and quick changes in the flow rate is essential. Thus, considering the importance of the unsteady flows in the various evaluation of WDNs, during the past three decades, efforts are made to develop hydraulic analysis methods for unsteady flows (Ammar & Al-Zahrani 2015; Abdel-Gawad & Djebedjian 2020; Bohorquez *et al.* 2020a).

Various techniques have been introduced for analyzing the unsteady flows in pipelines and WDNs. MOC, implicit finite difference method (IDM), explicit first-order finite difference method (FDM), Lax scheme, second-order FDM, MacCormack scheme, finite volume method (FVM), spectral method (SM), and boundary-integral method (BIM) are some of these methods (Gao *et al.* 2018; Abdel-Gawad & Djebedjian 2020). Among the mentioned numerical methods, MOC is popular for handling hydraulic transients. Simplicity and reasonable accuracy in the solution, especially for one-dimensional unsteady flows with constant wave speeds, can be considered the advantages of this method. It converts the two partial differential equations (PDEs) of continuity and momentum into four ordinary differential equations that are solved numerically using finite difference techniques (Walski *et al.* 2003). The development of friction models, applying coefficients, and formulating in MOC, mixing MOC with other models, and employing different metaheuristic algorithms can be considered as some of the approaches for developing solution capability and applications of this method.

One of the complications of using MOC for solving unsteady flows is that the changes in the friction are affected by the alternations in the flow (Abdel-Gawad & Djebedjian 2020). In 1991, the development of the Brunone friction model significantly improved the unsteady friction model for transient flow inside the pipes (Brunone *et al.* 1991). In 1994, to find more reliable results in discharges (leaks and unauthorized use), Liggett and Chen added friction loss coefficients as extra unknown factors to the unknown properties (Liggett & Chen 1994). In 2004, the transient hydraulic idea for viscoelastic pipes was presented by Covas *et al.* (2004, 2005), and in the subsequent research studies, the constant coefficient *k*, which was a function of Reynolds number, was used in the Brunone friction formulation and led to higher solution speed (Vítkovský *et al.* 2007). Direct backward transient analysis MOC was introduced to reduce the error of modeling valves and other system parameters such as roughness factor (Haghighi *et al.* 2012).

The development of the MOC method using RMOC (Reconstructive MOC) added the ability to recognize local destructions such as thickness reduction due to corrosion in the water supply pipelines (Gong *et al.* 2014). In 2017, Capponi *et al.* studied leakage in a system with several branches, with the admittance matrix method and inverse transient analysis (Capponi *et al.* 2017). Rahmanshahi *et al.* showed that the developed MOC (which includes unsteady dynamic loss phenomena and viscoelastic characteristics of pipe walls) could adequately predict the transient pressure in polyethylene pipes (Rahmanshahi *et al.* 2018). In 2019, Sarkamaryan *et al.* developed a mathematical model for hydraulically analyzing unsteady flows in pipelines using the MOC. In this model, a modifying factor was considered for instantaneous and local accelerations and the wave speed. Then, a nonlinear optimizing was conducted to minimize the square difference of fluctuations in the calculated and measured pressures (Sarkamaryan *et al.* 2019).

Pressure fluctuations are a part of natural behavior in WDNs. The common perception is that transient events dissipate quickly in a network without significant consequences (Bohorquez *et al.* 2020b); however, previous studies show that analyzing these phenomena accurately can result in getting useful information about the performance of WDNs, such as leak detection, recognizing bursts and blockage of pipes in the network. Comprehensive investigations have been conducted on methods based on transient flow in pipes. However, the performed studies on the transient signal characteristics or pressure fluctuations are not significant for this approach (Lee *et al.* 2015). Thus, the development of MOC as a reliable method for different applications, including detecting leaks, bursts, and blockage of pipes, can be useful (Kapelan *et al.* 2003a; Covas *et al.* 2005; Shamloo & Haghighi 2009). Based on the research conducted by Covas and Ramos, by comparing the pressure fluctuations in different sections of a pipe and by considering the transient event, the number of leaks can be determined (Covas & Ramos 2010). It was indicated that by comparing the maximum fluctuations in pressure in both upstream and downstream waves of a blockage pipe, the damping level could be obtained (Lee *et al.* 2015). By considering Kelvin–Voigt viscoelastic parameters, it was shown that the calibration of the unsteady flow model in plastic pipes is independent of the elastic module, and it depends on pressure fluctuation parameters like the pipe length (Pezzinga *et al.* 2016).

Zhang *et al.* presented the HBMOC (head-based method of characteristics), and by analyzing fluctuations of the transient flow, they concluded that their proposed method for the transient flow analysis is more accurate, and in terms of calculation, it is four times more efficient than the conventional methods (Zhang *et al.* 2018). Moreover, in 2020, Abdel-Gawad and Djebedjian examined transient flow in viscoelastic pipes using the wave characteristic method (WCM) (Abdel-Gawad & Djebedjian 2020). For the first time, they discussed the influence of space steps and the number of pipe sections on the solution accuracy. In this research, the minimum space steps for the studied case were 0.36, 1.0, 2.0, and 2.72 m (Abdel-Gawad & Djebedjian 2020).

In another category of research studies, different metaheuristic methods were employed to improve the solution speed, reduce the computational cost of MOC, and solve unsteady flows. The hybrid optimization program (Kapelan 2002; Kapelan *et al.* 2003b), Shuffled Complex Evolution Metropolis algorithm (SCEM – UA) (Kapelan *et al.* 2007), GA and Particle Swarm Optimization (PSO) metaheuristic optimization methods (Jung & Karney 2008), Central Force Optimization (CFO) method (Haghighi & Ramos 2012), compound method of ITA and SA (known as Leak Detection Simulated Annealing or LDSA) (Huang *et al.* 2015), Gaussian function parameters (Sarkamaryan *et al.* 2018) and GA-Kriging alternative method (Sarkamaryan *et al.* 2020) can be considered as the used methods in these studies.

Most research in this field using the transient flow solution has been limited to numerical methods, and few studies have evaluated these methods using laboratory data (Brunone *et al.* 2000; Vítkovský *et al.* 2007) or field data (Covas *et al.* 2005; Ranginkaman & Haghighi 2017), especially in plastic pipes. A comparison of previous experimental research with the present study is provided in Table 1.

Tank– . | Imperial College. London, UK . | Robin Hydraulics Laboratory, University of Adelaide, Australia . | The experimental facility at Instituto Superior Técnico, Lisbon, . | Hydraulic Laboratory, Shahid Chamran University of Ahvaz . | Sea-Based Energy Research Laboratory, Babol Noshirvani University of Technology, . |
---|---|---|---|---|---|

pipe–valve cases Authors . | Covas et al. (2004, 2005)
. | Wang (2002); Wang et al. (2003)
. | Portugal Covas & Ramos (2001) . | Iran, Rahmanshahi et al. (2018)
. | Iran Current study . |

Pipe type | HDPE | Copper | HDPE | HDPE | HDPE |

Pipe length L (m) | 272 | 85.4 | 38 | 158 | 24 |

Pipe nominal Diameter (mm) | 63 | 3″ | 50 | 63 | 50 |

Pipe internal diameter D (mm) | 50.6 | 72.94 | 45 | 55.4 | 42.6 |

Pipe wall thickness e (mm) | 6.3 | 1.63 | 2.4 | 3.8 | 3.7 |

Pressure head of the tank H_{T} (m) | 30–65 | 3.86 | 26 | 44 | 6.11 |

Head downstream the valve (m) | 0 | 0 | 0 | 0 | 0 |

Steady discharge Q_{0} (l/s) | 0.05–2.0 | 0.125, 0.25 | 1.0 | 0.8–1.1 | 2.14 |

Number of loops | 1 | 3 | 6 | 1 | 3 |

The purpose of the investigation | Analyzing the ITA method and leakage detection and location | Leakage and blockage detection | Increase reliability of the ITA method | Evaluation of viscoelastic effects of the pipe material and water | Analyzing the inclined pipelines |

Tank– . | Imperial College. London, UK . | Robin Hydraulics Laboratory, University of Adelaide, Australia . | The experimental facility at Instituto Superior Técnico, Lisbon, . | Hydraulic Laboratory, Shahid Chamran University of Ahvaz . | Sea-Based Energy Research Laboratory, Babol Noshirvani University of Technology, . |
---|---|---|---|---|---|

pipe–valve cases Authors . | Covas et al. (2004, 2005)
. | Wang (2002); Wang et al. (2003)
. | Portugal Covas & Ramos (2001) . | Iran, Rahmanshahi et al. (2018)
. | Iran Current study . |

Pipe type | HDPE | Copper | HDPE | HDPE | HDPE |

Pipe length L (m) | 272 | 85.4 | 38 | 158 | 24 |

Pipe nominal Diameter (mm) | 63 | 3″ | 50 | 63 | 50 |

Pipe internal diameter D (mm) | 50.6 | 72.94 | 45 | 55.4 | 42.6 |

Pipe wall thickness e (mm) | 6.3 | 1.63 | 2.4 | 3.8 | 3.7 |

Pressure head of the tank H_{T} (m) | 30–65 | 3.86 | 26 | 44 | 6.11 |

Head downstream the valve (m) | 0 | 0 | 0 | 0 | 0 |

Steady discharge Q_{0} (l/s) | 0.05–2.0 | 0.125, 0.25 | 1.0 | 0.8–1.1 | 2.14 |

Number of loops | 1 | 3 | 6 | 1 | 3 |

The purpose of the investigation | Analyzing the ITA method and leakage detection and location | Leakage and blockage detection | Increase reliability of the ITA method | Evaluation of viscoelastic effects of the pipe material and water | Analyzing the inclined pipelines |

Also, in most previous studies (Creaco *et al.* 2017; Page & Creaco 2019; Galuppini *et al.* 2020), development objectives were achieved by solving the unsteady flow by MOCs method in flat networks. In research (Triki 2018), numerical computations used the MOCs to discretize unconventional transient models for two operating conditions associated with up-and-down surge frames. It was first validated using a non-slope model. The slope was then applied as a static pressure difference at the beginning, end and along the entire length of the pipe (Triki 2018; Bohorquez *et al.* 2019). Flow rates will also be constant along the pipe. So, in the development of MOC for unsteady flow analysis, one of the topics that seem to be neglected is slope coefficients in inclined pipelines.

Calibration coefficients of the water distribution system, number, and intensity of pressure fluctuations due to the transient event are among the physical concepts that change under the influence of slope. On the other hand, most of the methods for increasing the solution speed by the MOC method have focused on using metaheuristic algorithms, and the definition of the initial simulation framework has received less attention. In this research, because most of the WDNs have variable topographies, the capability of the MOC for the hydraulic analysis of these networks has been developed. In the MOC, the two PDEs of continuity and momentum that govern the unsteady flow in the pipes are simplified and solved as a one-dimensional system by applying some coefficients. These coefficients are obtained based on some correlations, and they have specific values for each place and time. In this research, after expressing the usual correlations in the MOC solution, some parameters regarding the slope were applied to these coefficients. They were modified in computer code to add the ability to solve MOC for every slope value of pipes. The Babol Noshirvani University experimental results, Iran (NIT) water network, were used for the validation and calibration. This network had a constant head of 7 m, three loops with 2 in 2 m dimensions, and four pipes with a slope of 0.25. The obtained results for the developed code were also compared with the results from the EPANET hydraulic analysis software with similar conditions. The alternations in the intensity and number of pressure fluctuations and the damping time are all influential in predicting the major events. Therefore, after the calibration was done based on the steady flow rate of 2.21 l/s, the important results regarding the pressure in different network points, the inclined pipes, and non-inclined pipes, were extracted and evaluated. Furthermore, by considering the importance of space steps in the speed and precision of the solution, its influence on obtaining an accurate solution was studied.

## THE EXPERIMENTAL NETWORK STRUCTURE

The experimental network of NIT includes three loops with 2 in 2 m dimensions. The number of pipes and nodes are 12 and 10, respectively. The pipelines are made of polyethylene with a nominal diameter of 50 mm and a wall thickness of 3.7 mm. This network has a slope of 0.25 along the horizontal direction of the network. To keep the inlet pressure of the system constant, an elevated tank is employed, which has an overflow pipe over the tank. Figure 1 shows the introduced laboratory network with the schematics of the equipment and fluid flow passage. Some of the network characteristics are presented as follows:

The water level in the tank is 7 m higher than the end of the network. In this situation, higher flow rates with fewer fluctuations can be obtained in transient conditions.

By adjusting the outlet valve, a transient event is formed in the network.

The pipe lengths are equal (with reasonable precision); thus, their Courant number is identical.

The intersection where the valves and sensors are installed (using a bifurcation belt) is placed on points that are a multiple of the obtained space step based on the Courant number.

A turbine flowmeter is installed at the end of the system.

As shown in Figure 1, when the pump is turned on, the water is pumped from the tub into the elevated tank at the height of 7 m. After the water flows in the network circuit, the outlet flowrate of the network is read and recorded using the turbine flowmeter. To control the water flow in the network, a conical valve was used at the end of the network, and by adjusting the valve at the end of the passage, the total outlet flowrate is recorded as 2.21, 1.85, 1.45, and 0.85 l/s. The calibration coefficient of the turbine flowmeter in the experiments was obtained to be 0.968. After adjusting the steady outlet flow on the value of 2.21 l/s, the passing flow rate through each of the pipes is measured using a portable ultrasonic flowmeter. These tests were done three times for each of the pipes. The used flowmeter was Flexim Fluxus ADM6725 (FLEXIM.Co 2020), and the evaluation of flow speed was done based on the time difference for the transit of the ultrasonic signals between the sensors (time difference method). The data transfer was done without a computer using the data storage in a device. Table 2 shows the technical specifications of this flowmeter. It is noteworthy that during the tests, for keeping the head constant, the inlet valve of the tank was set in a way that the minimum amount of overflow occurred in the elevated tank. The equipment and structure of the laboratory network have been illustrated in Figure 2.

Item . | Value and specifications . |
---|---|

Measurement method | Ultrasonic measurement of flow velocity |

Flow velocity | 0.01–25 (m/s) |

Repeatability of reading | 0.15% of reading ±0.01 (m/s) |

Configurable pulse and relay output | Digital RS 232, RS485 |

Accuracy of reading | 0.1% of reading |

Damping rate | 1 (s) |

Clamp and IP rating | On M2N (on encapsulated IP 68 sensors require no pipe cutting) |

Wide operating temperature range | −20 °C to 60 °C |

Item . | Value and specifications . |
---|---|

Measurement method | Ultrasonic measurement of flow velocity |

Flow velocity | 0.01–25 (m/s) |

Repeatability of reading | 0.15% of reading ±0.01 (m/s) |

Configurable pulse and relay output | Digital RS 232, RS485 |

Accuracy of reading | 0.1% of reading |

Damping rate | 1 (s) |

Clamp and IP rating | On M2N (on encapsulated IP 68 sensors require no pipe cutting) |

Wide operating temperature range | −20 °C to 60 °C |

## THE GOVERNING EQUATIONS ON THE UNSTEADY FLOW

*et al.*2015):where,

*x:*coordinate along the pipe axis*t:*time*H:*piezometric head*Q*: flow rate*a*: celerity or elastic wave speed*g*: gravity acceleration*A*: cross-sectional pipe area*ɛ*_{r}: retarded strain*hf*: head loss per unit length*d/dt*: total derivative described by*d*/*dt*= ∂/∂*t**+**v*∂/∂*x**v*: absolute fluid velocity of the flow

The MOC has solved the set of differential equations. Completing the space step calculations, the boundary conditions of upstream (the tank with constant pressure) and downstream (the valve that creates the transient event) are considered.

### Method of characteristics

*x*), the time steps of the problem analysis will be equal to . This means in the current numerical calculations the Courant number will be equal to 1. The algebraic equation for the linear elastic conduits and the linear viscoelastic conduits can be written in this simplified form (Soares

*et al.*2008; Chaudhry 2014):

*C*,

_{n}*C*, and

_{a}*C*are constants that depend on the used numerical approach for expressing steady-state friction, the selected equations for the unsteady friction, and the chosen mechanical model for explaining the viscoelastic behavior of the ducts. These constants are defined as follows (Soares

_{p}*et al.*2008; Chaudhry 2014):

*et al.*(2008). The

*B*parameter is a function of the fluid and pipe properties:

*R*(the pipe resistance factor) is defined as follows (Covas 2003):

The solution algorithm of the network using the MOC has been shown in Figure 3. In the examinations, the following assumptions were considered for the pipes:

The behavior of the pipe wall material is linear-elastic, and the pipe material is high-density polyethylene (HDPE).

The fluid is homogenous, single-phase, and compressible.

The pipe diameter is constant, completely uniform, and constrained against all axial or lateral movements.

The flow is under pressure.

The influence of cavitation (including the separation of the fluid column and the trapped air bubbles), leaks, and the interaction of fluid and the structure are ignored.

#### Estimating friction loss

*et al.*2000; Covas 2003).

In these equations, Re indicated the Reynolds dimensionless number (), *V* is the average flow speed, *v* is the kinematic viscosity, *f _{q}* is the steady friction factor, and

*Є*indicates the average wall roughness.

*et al.*2000). In this method, the influence of fluctuations and unsteady friction is simulated based on the momentary local acceleration and momentary transitional acceleration.where

*k′*is the Brunone friction factor,

*θ*is the Relaxation coefficient, and the sign operator is the signal function. The flow parameter in the

*i*section and for the

*j*time for all of the inner sections of the pipe is modeled as follows:

Also, for the two ends of each pipe, some extra correlations that define the boundary conditions should be specified (Covas 2003).

#### Developing the correlations for an inclined pipe under-pressure

### Validation of the computer code

As discussed, for calibrating the studied network, the local losses from the fittings and the correct roughness values of the pipes have been applied. The equivalent length can be used to calculate the pressure loss fittings (Covas 2003). For this purpose, at first, using the EPANET.2 hydraulic analysis software, the system condition was simulated without considering the losses. Based on the initial conditions, the network was modeled in the computer code using the MOC solution method. Based on Table 2, using both of these methods, the distributed flow in each branch was calculated as identical values. The flow rate through each of these branches was measured and recorded using portable ultrasonic flowmeters three times.

To perform a calibration process and make the model conditions more realistic, all of the fittings were considered pipes so that instead of a three-way tee, three pipes that are connected were used. Thus, the typical structure of the network, which included 10 nodes and 12 pipes, was modeled with 32 nodes and 34 branches (Figure 5); then the equivalent length of each of these added pipes was considered, and finally, the measured flowrate values for the pipes and the pressure in the main nodes of the network were modeled. The final calculated values using the MOC and their accuracy have been presented in Table 3. According to Table 3, the average error of the developed computer code was 0.26% compared with the measured values, which shows a good solution accuracy. Moreover, in the MOC method, the initial flow rate based on the experimental values was selected to be 2.21 l/s, and a flowrate change of 0.03 l/s was considered for the transient event.

Pipe number . | EPANET (non-calibrated) . | MOC (non-calibrated) . | Measuring . | MOC (calibrated) . | Error (%) . |
---|---|---|---|---|---|

1 | 2.21 | 2.2096 | 2.21 | 2.2099 | 0.00 |

2 | 1.30 | 1.2978 | 1.74 | 1.7316 | 0.48 |

3 | 1.11 | 1.1074 | 1.40 | 1.3922 | 0.56 |

4 | 0.91 | 0.9118 | 1.07 | 1.0734 | 0.32 |

5 | 0.91 | 0.9118 | 0.47 | 0.4782 | 1.74 |

6 | 1.11 | 1.1049 | 0.81 | 0.8176 | 0.94 |

7 | 1.30 | 1.2977 | 1.14 | 1.1364 | 0.32 |

8 | 0.91 | 0.9118 | 0.47 | 0.4782 | 1.74 |

9 | 0.20 | 0.1931 | 0.34 | 0.3394 | 0.18 |

10 | 0.20 | 0.1929 | 0.32 | 0.3188 | 0.38 |

11 | 0.91 | 0.9117 | 1.07 | 1.0734 | 0.32 |

12 | 2.21 | 2.2095 | 2.21 | 2.2099 | 0.00 |

Pipe number . | EPANET (non-calibrated) . | MOC (non-calibrated) . | Measuring . | MOC (calibrated) . | Error (%) . |
---|---|---|---|---|---|

1 | 2.21 | 2.2096 | 2.21 | 2.2099 | 0.00 |

2 | 1.30 | 1.2978 | 1.74 | 1.7316 | 0.48 |

3 | 1.11 | 1.1074 | 1.40 | 1.3922 | 0.56 |

4 | 0.91 | 0.9118 | 1.07 | 1.0734 | 0.32 |

5 | 0.91 | 0.9118 | 0.47 | 0.4782 | 1.74 |

6 | 1.11 | 1.1049 | 0.81 | 0.8176 | 0.94 |

7 | 1.30 | 1.2977 | 1.14 | 1.1364 | 0.32 |

8 | 0.91 | 0.9118 | 0.47 | 0.4782 | 1.74 |

9 | 0.20 | 0.1931 | 0.34 | 0.3394 | 0.18 |

10 | 0.20 | 0.1929 | 0.32 | 0.3188 | 0.38 |

11 | 0.91 | 0.9117 | 1.07 | 1.0734 | 0.32 |

12 | 2.21 | 2.2095 | 2.21 | 2.2099 | 0.00 |

## RESULTS AND DISCUSSIONS

In this research, the hydraulic conditions of steady and unsteady flows have been studied for a pipe network with various elevations. A network model has been simulated using the MOC in the MATLAB software and calibrated based on the measured values from field examinations.

The wave speed is assumed to be 314 m/s, and the Darcy–Weisbach factor *f* is assumed to be 0.0015 for all the pipes (based on the Reynolds number).

Recognizing the pressure wave parameters is very important in predicting significant events such as the location and amount of leakage and the place of blockage or burst. The intensity and number of fluctuations in a specific period and the damping duration of pressure waves are some of these parameters. Analyzing the changes in these parameters allows us to predict the status and intensity of some phenomena by having the data from a section or a pipe in the network and its position. Here, the influence of the distance between the section of a pipe and the location of the transient event will be analyzed on the pressure damping. Also, the effects of the pipe slope on the intensity and number of fluctuations will be examined.

Moreover, defining the framework and the comparison criteria in which the solution accuracy and the solution speed both matter is crucial. The size of space steps is one of the parameters that are decreasing leading to higher solution accuracy, but the solution speed will be reduced with an exponent of 2. Thus, the influence of space steps on obtaining accurate results will be evaluated. For the pipes with a length of 2 m, the space steps of 0.5, 0.25, 0.1, and 0.05 will be tested.

### The influence of the distance between the pipe and the transient event on the pressure damping

In this part, the intensity and the number of fluctuations in different sections of the pipe will be investigated. To compare other states, pipes 2 and 6 (without slope) and pipe 11 (inclined pipe) will be examined (Figure 2). Between the start and the end of pipe 11, there is a 0.5 m of elevation difference. In other words, a slope of −0.25 for the unit of length. Fluctuations and pressure differences in the middle section of these pipes have been presented in Figure 6. The unsteady flow simulations of the network were performed considering Δ*t* = 0.00125 s.

As shown in Figure 6, since pipe 2 is placed further from the transient event compared with other pipes, it reaches the damped condition faster than the other pipes. Figure 7 shows the magnification for a part of the time interval (from 1 to 1.6 s). Examining Figure 7 shows that the number of pressure fluctuations in each pipe negatively affects these fluctuations' intensity. Therefore, in the pipe networks, with a higher distance from the transient event, the number of pressure fluctuations in a specific period will be higher, although their intensity reduces. Also, as this distance increases, the damping speed rises.

To evaluate the damping speed of the pressure fluctuations in different parts of the network, the criteria fluctuation amplitude for the pressure fluctuations is considered to be 0.2 m (Figure 8).

### Influence of slope on the intensity and number of pressure fluctuations

Figure 9(a), 9(c) and 9(e) shows the pressure fluctuations in the initial moments of the transient event (2–2.3 s), and Figure 9(b), 9(d) and 9(f) shows the pressure fluctuations some moments after the formation of the transient event and before the damped state is reached (the time duration in the range of 7–7.3 s). Comparing the results of pipes 2 and 6 indicated that with a higher distance between the pipe location and the transient event, the intensity of the fluctuations increases. Also, examining the fluctuations in different sections of these pipes shows that the pressure fluctuation intensity will be lower in that section with a higher distance between the section and the transient event. This is clear in Figure 9(c) and 9(e). However, unlike the two previous cases, the pressure fluctuations in different sections of pipe 11 (Figure 9(a)) show that this parameter is higher at the start of the pipe (the section which is further compared with the transient event). This is due to the influence of the slope on the pipe. Also, the unsteady flow simulations of the network were performed considering Δ*t* = 0.00025 s.

After some time, the pressure difference in different pipe sections becomes closer to the final pressure difference. For the pipes with more intense pressure fluctuations, this happens faster. Based on this, in the time duration between 7 and 7.3 s, pipe 11 reached the final pressure difference in the pipe sections (based on Figure 9(b)), while pipes 2 and 6 did not get this condition at that time (based on Figure 9(f) and 9(d)). Also, evaluating the results in different periods shows that the converging period of the pressure towards the final pressure in each pipe has a negative relationship with the intensity of the pressure fluctuations in that pipe. Although, this relationship is not valid for different sections of a pipe.

### The influence of the space steps on the solution accuracy

One of the most influential parameters in the solution accuracy and speed is the space step size to solve the transient flow using the MOC. Considering the elastic pipe, the highest acceptable error for obtaining reliable results is selected to be 2% (Ramalingam *et al.* 2009; Abdel-Gawad & Djebedjian 2020). For this purpose, by reducing Δ*x*, the results of the MOC will be more precise. On the other hand, for the simulation of the transient flow in viscoelastic pipes with the MOC, the time interval (Δ*t*) should be less than half of the minimum retardation time (Covas 2003). To keep the solution stable and considering the Courant number equation (CN = *a*Δ*t*/Δ*x*), the reduction of space steps leads to decreasing the time step; thus, in a specified period, by increasing the space steps, the solution speed decreases.

This section aims to examine and specify a proper space step for reaching acceptable responses in less time. To evaluate the influence of the space steps on the results, the pressure fluctuations in three different points, including node 3 (the start of pipe 2), node 22 (the start of pipe 6), and node 15 (the start of pipe 11) have been evaluated (within the node at distances 12.5, 8.5, and 6.5 m from the transient event, respectively). The selected space steps include 0.5, 0.25, 0.1, and 0.05 m. Considering the equality in the length of the network pipes (2 m), the number of sections based on each of these space steps included 4, 8, 20, and 40 nodes, respectively.

Comparing the solution behavior with different space steps for nodes 3, 15, and 22 for the first 3 s has been presented in Figure 10. As it can be observed, Figure 10(a), 10(c) and 10(e) shows the pressure fluctuations for a 5 s interval. It is noteworthy that in all of the presented plots, the transient flow has been applied after 1 s. Also, in Figure 10(b), 10(d) and 10(f), by concentrating on the time interval between 1.4 and 1.5 s, the figures are magnified to investigate the changes in the plots based on different space steps more precisely. As can be seen, by reducing the space step, the solution accuracy increases, and more details are observed from the behavior and fluctuations of pressure. Furthermore, by reducing the space steps from 0.1 to 0.05, no significant change is observed in the plots, thus for continuing the solution, the space step of 0.1 can be used.

For evaluating the results mode precisely, the obtained results with the space step of 0.05 m were chosen as the reference, and the relative error of the solution with other space steps was compared with it. For this purpose, 4, 8, and 20 segments (*N*_{s}) are used, which are based on the space steps of 0.5, 0.25, and 0.1. As can be observed in Figure 11, the solution with *N*_{s} = 4 exhibits high error values, so that the error in node 15 reaches up to 9%. On the other hand, when *N*_{s} = 20 is used, the error is always lower than 1%, and that verifies the reliability of the solution.

## CONCLUSION AND FUTURE WORKS

The hydraulic analysis of steady and unsteady flows using MOC is a powerful solution for many practical applications, and the development of this method has always been an area of interest for many researchers. The development of MOC as a reliable method for different applications, including the detection of leaks, bursts, and blockage of pipes, can be useful. Also, hydraulic analysis of WDNs with various demand consumption patterns in different zones is another application of this method. Most of the developmental research studies on this method are related to improving the solution speed using metaheuristic algorithms, developing computer codes, modifying the MOC, and performing field and laboratory tests for calibration and enhancing the estimation methods of the solution, such as the unsteady friction. Thus, for this purpose, in this research, the data of the laboratory network that was constructed in NIT were numerically analyzed using the MOC method by developing the computer code. Since this network has various elevations, it can represent the conditions close to the real WDNs. The most important findings and results of this study can be summarized as follows:

The losses were specified as the equivalent length and the network was calibrated by transforming the fittings to a pipe structure.

The C

^{+}and C^{−}equations in the MOC were developed and introduced for inclined pipes. This method is used in the hydraulic analysis of WDNs to manage pressure, reduce leakage, etc.In the pipe networks with a higher distance between the pipe and the transient event, the number of pressure fluctuations in a specified period increases, but the intensity of these fluctuations decreases. Also, as this distance increases, the damping speed rises.

In the inclined pipes, lower pressure fluctuations are observed at the end of the pipe, meaning that it is closer than the transient event but has a lower elevation.

After some time, the pressure difference in different pipe sections becomes closer to the final pressure difference. In the pipe, with higher pressure fluctuation intensity, this happens earlier.

The number and duration of each pressure fluctuation are identical in different space steps.

The changes of

*N*_{s}in the range of 4–40 have been examined. Based on the examinations, even with one segment, acceptable results can be obtained. Also, it was observed that increasing*N*_{s}leads to numerical results closer to the real values.To increase the solution speed and get acceptable solution accuracy, the proper space step was selected. The results indicated that by choosing a space step (

*x*= 0.1 m or smaller), results with less than 1% error can be obtained.

Analysis of changes in the intensity and number of pressure fluctuations in a specific longitudinal range to detect leakage (as a transient event) can be considered as an important goal of this research in future work. Also, combining the two methods of using metaheuristic algorithms and modeling by specifying a proper space step can be considered to improve the solution speed of the MOC method.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.