Evaluation of service pressure regulation strategy on the performance of a rural water network based on pulse demand; using the method of characteristics

Symbols Roman

A

cross-sectional pipe area (m²)

a

celerity or elastic wave speed (m/s)

C

steady and unsteady constants (–)

D

pipe diameter (m)

g

gravitational acceleration (m/s²)

H

head loss (m)

I

pulse intensity (m³/s)

L

pipe length (m)

N

number of customers (–)

NDF

night-to-day factor (hr/day)

P

pressure head (m)

Q

flow rate, leakage (m³/s)

Sl

slope of the pipeline (–)

T

pulse duration (s)

t

time (s)

V

pulse volume (m³)

v

absolute velocity of the fluid (m/s)

x

coordinate along the pipe axis (m)

Greek symbols

β

model parameter (–)

Strain (–)

λ

frequency of the received pulse (1/s)

τ

pulse time arrival (s)

Superscripts

leaks

leakage output

α

model parameter

N

leakage exponent

Subscripts

d

day

f

loss friction

h

homogeneous

i

input or initial points

_j

final points

k

pipe number

loss

water loss

mean

average pipe pressure

MNF

minimum night flow

n

reverse characteristic line

p

forward characteristic line

r

retarded

Water distribution networks (WDNs) are complex systems of interconnected nodes and pipes. The design of WDNs should be such that water is transferred to the consumers at appropriate times and with acceptable quality (by the standard) (Creaco et al. 2019). Increasing the network life, reducing leakage, reducing bursts, reducing energy costs, increasing revenues and ultimately increasing customers' satisfaction are among the important goals to increase the efficiency of WDNs; therefore, in recent decades, WDN pressure management (PM) has been increasingly considered by researchers and managers of water supply companies (Bello et al. 2019).

Given the above priorities, leakage reduction has always been the most important goal of research: a phenomenon that is inevitable even in new networks (Taha et al. 2020). As a result, PM methods can be considered one of the best methods to achieve these goals (Creaco et al. 2019). In low pressure, the flow rate less than the demand is distributed in the nodes. In contrast, excessive pressure has various adverse effects such as increased leakage, a pipe burst and shortening of infrastructure life (Vicente et al. 2015). The optimal case occurs when the pressure in all nodes is proportional to the demand; this is called the optimum pressure profile or target pressure level.

Over the past decades, the effectiveness of the pressure control approach and its encouraging results have been repeatedly expressed in various studies. Studies suggest that this method can be an effective tool for decision making and determining the optimal solution to improve WDN resilience (Covelli et al. 2015; Samir et al. 2017). Conversely, knowing various methods for implementing PM is of utmost importance. Several methods can help reduce or eliminate the excessive pressure on the network. These include reducing the pressure head, creating pressure zones, using flow control valves or pressure-reducing valves (PRVs). The focus of initial research has been on determining the number, location, and setting of such PRVs while maintaining the minimum required operational pressure at every node and minimizing the leakage in the system as problem objectives (Javadinejad et al. 2019b). Replacing the pipe and installing the control valve (Creaco & Pezzinga 2014), providing a mixed-integer nonlinear program (MINLP) for minimizing the mean pressure through optimal placement and operation of the control valves (Pecci et al. 2015), finding the optimal location of the turbines (Corcoran et al. 2015), optimizing the water storage level in the tank (Gupta et al. 2018) and applying variable speed pumps (VSP) (Monsef et al. 2018) are other methods to implement pressure control programs. In conclusion, the type of pressure settings includes fixed outlet modulation (FO), time-based modulation (TM), flow-based modulation (FM) and remote node-based modulation (RNM) (Vicente et al. 2015). The advantages of remote control of the target node using PRVs in WDN were first presented (Campisano et al. 2010). In the following studies, the dimensionless procedure for the calibration of the algorithms (Campisano et al. 2011) or development and presentation of a new control algorithm (Creaco & Franchini 2013) and then numerical simulation of the application of target node control in a case study (Creaco et al. 2016, 2017) were performed. Real-time control (RTC) of pressure is a more effective method due to the increased ability to control the service pressure at the critical node because the pressure is controlled at shorter time intervals (Creaco et al. 2018). In the study (Galuppini et al. 2019), an advanced control algorithm with a 5-minute control time step was used as a criterion. Determining the appropriate time interval for pressure control using numerical simulation of the dynamic behavior of the WDN allows the simulation of various case studies.

Service pressure regulation is a complete approach that encompasses all the goals of other PM methods. Service pressure regulation in WDN should be high enough to provide customers' satisfaction (Galuppini et al. 2019). In this approach, a node with the lowest pressure, normally located at the end of the WDN, is selected as the controlled node or the target node. Pressure control conditions are then set by PRVs located at the first point of the WDN, and the pressure at the target node is continuously equal to the minimum head. In the physical implementation, it is assumed that using a function or programmable controller, the deviation of the actual pressure value from the set point is measured. Aperture settings of PRVs are applied by the actuators every few minutes. This approach is called Advanced Pressure Control (APC).

The steady-state hydraulic model limits the reliability and efficiency of WDNs. In these models, only a short-term sample of a subset of thousands of different and unknown behaviors of hydraulic data are approximated. Therefore, calibration results cannot accurately represent system conditions for a wide range of operational conditions that may occur. A real-time dynamic hydraulic model that continuously considers online hydraulic measurements will provide realistic predictions (Abu-Mahfouz et al. 2019). Transient pressure and water column theory provide useful information for describing transient modes in the pipes (Kaltenbacher et al. 2017).

Social parameters are the most important element in rural areas. Also, the indexes related to spatial changes in rural areas are affected by climate change (Javadinejad et al. 2019a). Therefore, research in this field can help managers in better strategy planning at this event. In some areas (even Mediterranean countries), customers try to deal with water scarcity by private tanks (De Marchis & Freni 2015; Negharchi & Shafaghat 2020). Private tanks often magnify the actual needs of customers. They cause the network to work in conditions far from its design. In this way, in the lower parts of the WDN, until the tanks are filled, the flow is allocated much more than the design. Thus, supplying to the low-pressure districts of the network is not done well. In such cases, the pressure on the network is generally less than the design and is controlled by the levels in the tanks (Freni et al. 2014). In addition, subsequent processes of filling and emptying private tanks change network pipes' pressure and flow rate and cause changes in energy production, leakage, and pressure fluctuations in the WDN (De Marchis & Freni 2015). This configuration reduces the application of conventional steady-state models since the private tank filling process causes a constant change in the hydraulic behavior of the WDN. Dynamic and pressure-driven models are required to analyze this continuous change in network state variables (Freni et al. 2014).

Conversely, highly variable demands in a WDN cause significant fluctuations in the PRV set point. Transient phenomena are then created and propagated through the network, which may lead to water quality problems, unequal distribution of resources among customers, and premature deterioration of the pipeline infrastructure (Freni et al. 2014). By determining the appropriate spatial step in the pipes and avoiding additional parameters such as pipe elasticity or fluid compressibility, computational effort and complexity are significantly reduced (Kaltenbacher et al. 2017; Negharchi & Shafaghat 2021).

The dynamic model uses real-time measured data for current network conditions and automatically sends control signals to various network components. This approach adjusts the network's performance and makes it more efficient (Abu-Mahfouz et al. 2019). A model was developed in previous studies and an additional pressure control module was implemented to analyze PRVs in a fully dynamic numerical framework. This model was robust and reliable in implementing PMAs in the network. This model was used for a district of the Palermo (Italy) network. Pressure and flow data were monitored and were available for model calibration (Freni et al. 2014). De Marchis & Freni (2015) developed a special dynamic module capable of simulating PATs for energy recovery in Palermo WDN. Numerical simulation by MOC method based on nodal demand model reproduced the effect of private tanks (De Marchis & Freni 2015). Rasekh and Brumbelow developed an evolutionary-computation-based dynamic optimization model to recognize the optimum time-dependent qualitative measures in an emergency (Rasekh & Brumbelow 2015).

One of the ideas studied in the unsteady flow analysis due to demand pulses is to examine the consumption pattern data in short-time intervals. In one study (Creaco et al. 2017), the potential of unsteady flow modeling for simulation of pressure control in WDNs was investigated. The created model was a simulator of unsteady flow for generating nodal demand pulses and dynamic adjustment of PRVs in the network. Different pressure control scenarios were solved using the MOC with the Manning formula in this study. The network consists of 26 nodes with unknown heads with a ground elevation equal to sea level and a reservoir with a ground level of 35 m above sea level. A critical node and a control valve were used in this network. The results of the skeletonized model of a real network show that the unsteady flow model provides a sounder description of the amplitude of the pressure head variations at the controlled node (Creaco et al. 2017). Due to the importance of unsteady flows in the evaluation of WDNs, various methods have been introduced to analyze unsteady flows in pipelines and WDNs (Abdel-Gawad & Djebedjian 2020). Among numerical methods, the most common method in calculating unsteady flow is the method of characteristics (MOC). The research results (Galuppini et al. 2019) proved that demand pulses are an efficient factor for modeling unsteady flows; therefore, the control error of all new control algorithms with the presence of demand pulses was reduced to about 40%.

Several software tools provide steady-state hydraulic analyzing. These tools include EPANET and WaterGems software. In steady-state analysis, network information such as tank level, consumption, pump and valve operation do not change over time; whereas, regarding the normal conditions in WDNs, it is essential to consider the rapid change of flow rate and unsteady pressures. It is necessary to simulate the effects of controlling node pressures in WDNs using a suitable model. Regarding the nature of demand pulses, unsteady flow modeling enables considering the hydraulic transients in WDNs as an appropriate tool for creating dynamic adjustments.

This study has developed a hydraulic model for sloping networks as a laboratory network as most WDNs have variable topography. Then, considering the case study network (Gavankola village), the solution's effective coefficients, including friction loss coefficients of steady-state and unsteady-state flows and minor loss coefficients, were calibrated. Due to the mentioned advantages in using the dynamic modeling of WDN by MOC, the capabilities of pressure-driven analysis and analysis by automatic adjustment of PRVs or VSP have been added to a computer code. This research aims to regulate the service pressure to achieve all the goals of PM. One of the innovations is combining three approaches to determine the demand pulses pattern: water demand time series, consumption accumulation pattern of each component (user), and leakage pattern. The measuring equipment included an ultrasonic flowmeter and a pressure sensor with the ability to send data remotely every 5 minutes. Most networks of the previous studies (Creaco et al. 2017, 2018) had a high average consumption flow rate. Therefore, the network under study must be a large one. However, in practice, the number of nodes and network pipes was much smaller than the normal number of real networks, which means that many details were overlooked in the network, which can sometimes affect the results. However, the selected network of this study, with a length of about 8 km and an average water demand of 4.2 L/s, has all the features of a complete network. Due to the topographic conditions, this network has excess pressure, while there is low pressure in another district.

Consumption curves are used to determine the nodal demand. Two different approaches can be adopted to present the consumption curve. The first approach, known as the top-down approach, is the most definitive and common. In this approach, the amount of water entering the network is measured at different time intervals. Then the consumption multiplier is created based on the average consumption. Finally, considering this correction factor, the nodal demand is obtained at any time. The implementation of this approach is simple. However, it has failed to succeed in large spatial-temporal changes in the water demand (Creaco et al. 2015). The second approach, known as the bottom-up approach, creates a consumption pattern or demand time series based on the aggregation of the nodal demand process of individual customers' demands. This approach is called demand pulses. Measuring the flow rate using an ultrasonic meter with small-time damping can clearly respond to demand changes in the network.

In recent years, hydraulic analysis of unsteady flows using the concept of demand pulses has received more attention (Creaco et al. 2017). In the demand pulses approach, the amount of common demand is called pulse intensity and the demand withdrawal time is called pulse duration. There is a positive correlation between these variables in some aspects and types of household water consumption. For example, when filling a bathtub, which can be classified as a long-term pulse, the customer typically opens the valve completely. Therefore, a high-intensity pulse is produced to fill the bathtub quickly. Similarly, washing hands is classified as a short-term pulse, and the customer tends to open the valve in a limited way. Therefore, the opening of the valves occurs in high or low intensities, which can be associated with high and low time pulses, respectively (Creaco et al. 2017). Therefore, input demand pulses into a network can be obtained through spatial–temporal accumulation of pulses (Creaco et al. 2015).

Pulse generation is used to define the demand in each node of the network. In this method, to create water demand, pulse properties are considered in terms of pulse time arrival (τ (s)), pulse durations (T (s)) and pulse intensity (I (m³/s)). Volume pulse (V) can be obtained by multiplying the pulse durations in pulse intensity (T.I). In pulse generation, the pulse durations (T) and the pulse intensity (I) are independent of the pulse time arrival (Figure 1(a)). In particular, the nodal demand at any time step is obtained as the sum of the intensities of all active pulses (Figure 1(b)). This feature is affected by the uniform parameters of T and I.

Most customers of each network have similar consumption behavior, which is called homogeneous customers. The different but homogeneous customers' demands can be grouped in the same upstream node in such circumstances. Therefore, for a given node with a homogeneous number of consumers N_h, if λ_i (s⁻¹) is the frequency of the received pulse corresponding to the customer n_i, then the frequency of the received pulse λ (s⁻¹) in the same node is obtained from the sum of the consumers as ⁠. Moreover, if there are different consumption categories, i.e., different groups of homogeneous consumers, in a particular node, the received frequencies should be considered separately for each type (Creaco et al. 2015) (Figure 2). Finally, for a perfectly homogeneous network or a perfectly homogeneous district metered areas (DMA), after generating a pulse in the generic network's node or DMA, the pulses can be collected on a specific time scale (for example, one second or one hour). Typically, the described model can be used to generate long pulses (for example, one month) in which model parameters including T and I can be assumed to be constant. Instead, the parameter λ changes according to a particular pattern to describe the daily changes in pulse input. Finally, the daily consumption pattern of each node is obtained (Creaco et al. 2017).

The lowest flow in the DMA over 24 hours is called the MNF, which generally occurs during the early morning due to reduced consumption. In this method, considering the customers' legitimate night-time consumption (LNC) of a DMA and the amount of leakage at that time (Q_loss) is obtained by subtracting them from the MNF value.

Creating WDNs as DMAs is a known leakage monitoring technique that was introduced in the UK in the early 1980s (Vicente et al. 2015). This technique involves dividing the network into locations with permanently defined boundaries and continuous measurement of input and output flow from DMA. Therefore, during low demand, such as overnight, leakage can be estimated, and bursts can be detected in any DMA. This helps water companies plan for pipeline repair and replacement solutions. In the UK, leakages have dropped by up to 30% over the past 25 years in areas where DMA has occurred (Wright et al. 2015). Water pressure control when entering DMA is one of the complete forms of PM. This approach is called PMAs (Vicente et al. 2015).

As DMAs or PMAs are being designed, several factors including size (geographic areas, length of pipes, and number of connections), topological considerations (number of inputs, areas supplying water to an adjacent area, and boundary features), changes in elevation, type of consumption specifications (uniform demand, change in characteristics over time with fixed patterns, and change in shape over time with random demand) and water quality considerations (Vicente et al. 2015) should be considered. This research tries to create PMAs based on the pressure behavior of pipes. This method, known as DMA performance with dynamic topology, combines several advances in hydraulic monitoring, modeling, and controlling, and can address the potential disadvantages of DMA while maintaining its advantages (Wright et al. 2015).

MOC discretizes and solves the hyperbolic partial differential equations (PDEs) governing the transient flow into four ordinary differential equations by the finite difference method. Simplicity and reasonable accuracy in solving equations, especially problem solving of one-dimensional unsteady flow with constant wave velocity, is one of the advantages of this method (Negharchi & Shafaghat 2021).

Two independent equations are needed to determine the flow parameters at any time and space to describe the unsteady flow. Equations of mass and momentum survival describe unsteady flows in closed conduits. Considering the two dynamic terms of unsteady loss and dynamic effect of pipe-wall viscoelasticity, these equations are PDEs.

Signs ′, ″, ′′′ refer to steady-state friction, unsteady-state friction, and rheological behavior of the pipe wall, respectively. The numerical definition of each coefficient has been provided previously (see table 1 in the reference Soares et al. 2008). The schematic of the mathematical model of Equations (9) and (10) is as shown in Figure 3. The boundary conditions and relationships of PRVs to solve the unsteady flow model are presented in Table 1 (as from Janus & Ulanicki 2018).

In this research, to estimate steady flow loss, Hagen–Poiseuille explicit correlations for slow flow regime (Re < 2,000) and Colebrook equation for turbulent flow regime (Re > 4,000) were employed. The effects of unsteady friction fluctuations are also simulated based on the momentary local acceleration and momentary transitional acceleration (Negharchi & Shafaghat 2021).

Moreover, additional relationships must be specified at both ends of each pipe to describe the boundary conditions. As mentioned earlier, in this study, three consumption patterns were used to generate pulsed node demands by measuring the flow rate (at the outlet of the reservoir) and network pressure (at the target point) and the results of the MNF method. The dynamic hydraulic model flowchart is presented in Figure 4.

For validation and calibration purposes, MOC results were compared with the experimental data. The experimental WDN of Babol Noshirvani University of Technology (NIT) was used to extract these data. As this network has a variable altitude structure, it can represent the conditions of the real WDNs. This network consists of three square loops with dimensions of 2 meters by 2 meters. The number of main pipes is 12, and the number of nodes is 10. Polyethylene pipelines are used, and each pipe has a nominal diameter of 50 mm with a 3.7 mm wall thickness. This network consists of four sloping branches with a transverse slope of 0.25. In order for the input pressure of the system to be constant, an elevated tank is employed about 7 meters above the end of the network. Setting the output valve, a transient flow is created in the network. Figure 5 shows the introduced experimental network with a schematic of the equipment and fluid flow passage.

After setting the steady output flow of the network at 2.21 l/s, the flow rate through each of the pipes was measured using a portable ultrasonic flowmeter. This experiment was performed with three replicates for each pipe. The flowmeter used was the Fluxus ADM6725 model made by Flexim company (FLEXIM 2020).

For validation and calibration, the results of MOC were compared with the results of EPANET hydraulic analysis software as well as experimental data. It has been presented previously (see table 2 from reference Negharchi & Shafaghat 2021). Then, based on these initial conditions, the Gavankola WDN was modeled using the MOC solution method in the computer code. In this regard, model structure concerning nodal demand generation, unsteady flow modeling, and APC were implemented in the MATLAB environment. Calculations were carried out with a laptop with an Intel Core i7-2630QM CPU with a frequency of 2.0 GHz and 8 GB of RAM. The results of the MOC showed that, in the studied network, with 20 segments for each pipe, the fastest response with an error of less than 0.5% is achieved. The turbine flow meter calibration was 0.968 after the experiments (Negharchi & Shafaghat 2021).

After determining the effective coefficients of unsteady flow, longitudinal friction loss coefficients and minor loss coefficients of steady flow are among the numerical solution unknowns that can be achieved using performing experimental experiments. To do so, pipes and fittings similar to the structure of the Gavankola rural network (made of polyethylene and the same density) were tested in the laboratory network of the NIT. The longitudinal loss coefficient was considered based on the Hazen–Williams equation coefficient, and minor losses were applied as equivalent lengths. To connect the pipe to the contraction, sudden, 90 deg. bend, and tee pipe, the minor loss coefficients are considered as k₁, k₂, k₃, k₄ and k₅, respectively, according to Figure 6. As the number of unknown variables is more than known variables, the test and numerical solution process was performed in several steps. First, pressure data were collected at four network points by performing three tests at eight different outflows.

By calibrating laboratory data with numerical solution results, minor loss coefficients were determined for different Hazen–Williams coefficients (from the value of 100 to 150). Then, by displacing one of the sensors and measuring the pressure changes at both ends of one 90 deg. bends in the test with four different flows, the value of the coefficients was determined according to Table 2. The calibration accuracy was very high in both steps of the experiment so that the error of the difference between the measured and modeled values was less than 1% in both tests. The pressure sensors used are the Atek BCT-22-6B-A-G1/4-C-S30-E153 model (Atek Co. 2021) and their specifications are presented in Table 3.

In this study, a real WDN has been selected to apply the above methodology in order to apply service pressure regulation. The Gavankola village WDN is located in the north of Iran in Babol city and serves up to 900 people. In the Gavankola WDN, the groundwater is chlorinated and pumped into the elevated water reservoir. Then, the water in the reservoir is distributed in the network by gravity. The height difference of the distribution system is about 91 meters and it has 129 pipes with a length of 8 kilometers, one reservoir and a pump station with one pump. This network is made up of nine commercially available pipe sizes (32 to 125 mm). Due to altitude difference there are two consumption patterns, and the presence of a household water tank and a pumping system in the house of this area. The pipe configuration and altitude of WDN are shown in Figure 7.

The selected area in this study has features that have been less addressed in previous research. These include being in a rural area, having tourist attractions, and various water distribution pressures. The importance of this issue is that the water distribution with variable pressure will affect the service demand and consumption pattern. It is worth noting that the probability of burst and leakage in the WDNs with high-pressure fluctuations is much greater.

This research is based on measuring the data over 3 months and with the frequency of data transmission every 10 minutes for the input flow to the network and frequency of data transmission every 5 minutes to measure the target point pressure of the network. The ultrasonic flowmeter used was the SonixMeter SL300 model made by the Soleyman company (Soleyman Co. 2020). The flow velocity measurement evaluates the time difference of ultrasonic signals between transit-time sensors (table 4 in Negharchi & Shafaghat 2020). In this regard, after installing the clamp sensors on the pipe, an accurate portable ultrasonic flowmeter (pre-calibrated) was used to calibrate the flow rate. In addition, the pressure sensor records the pressure while it is connected to a data logger used to record the remotee values.

The pressure sensor used is the Rayan logger (Rayan Electronic 2020), and its specifications have been presented previously (see table 5 in Negharchi & Shafaghat 2020). The data loggers, equipped with an internal antenna, send gathered data via the Short Message Service (SMS) at a user-defined period (Figure 8).

A combination of three approaches of water demand time series, consumption accumulation pattern of each component and leakage pattern has been used to create the demand pulse pattern in nodes. Based on the approach, first, by measuring the inlet flow rate of the network, the pattern of the demand factor of the whole network can be determined. Then, based on the required flow of each pipe, the demand of each node is assigned. Although this approach is simple to implement, three factors can challenge the accuracy of this model as demand pulses:

• 1.

  The inlet flow rate of the network due to leakage cannot indicate the consumption pattern. This challenge is solved by measuring the leakage.

• 2.

  The consumers' heterogeneity in a network makes this approach confront mismatched pulse intensity and volumetric pulse synchronization. Given that more than 97% of customers are household customers, this challenge does not exist in this network.

• 3.

  If there are various service pressures in the network, the first approach is problematic. This challenge can be addressed using the second approach.

The inlet flow rate of the network was collected over 3 months (Negharchi & Shafaghat 2020). In the WDN of Gavankola village, the altitude difference between the reservoir and the consumption points is very diverse, resulting in difficult operating conditions and low pressures for some customers during peak hours. Therefore, most customers use water storage tanks to store water at night within the specified limits.

Part of inlet water entering the network always leaks out of the network pipes. The leakage from the network varies at different times of the day. The leakage can be calculated using the MNF method (Negharchi & Shafaghat 2020). Figure 9 shows the network leakage pattern. The difference between the two inlet flow rates to the network and the amount of leakage indicates the flow rate consumed by customers (Figure 10(a)).

The pressure sensor was installed at the critical node of the network, with the least service pressure to customers (Figure 6). Since this pressure is service pressure or discharge pressure, the trend of the pressure pattern (unlike the pressure and flow behavior at the network input) will indicate the pattern of demand pulses. Considering the customers' demand specified in Zone B (Figure 6), the pattern of demand pulses of customers in Group B was obtained according to Figure 10(b). By differentiating the demand pattern of Groups A and B (Figure 10(a)) from the pattern of zone B, the pattern of demand pulses of Group A (Figure 10(c)) will be obtained.

Before implementing pressure control methods and strategies, it is first essential to analyze the hydraulic conditions of the network over 24 hours. To analyze the hydraulic behavior of the Gavankola village network, the coefficients of steady and unsteady flows were first obtained by analyzing the laboratory network using the MOC. PRVs and pumps based on the boundary conditions presented in Table 1 were added to the computer code. The way for the application of a skeletonized network structure includes the following:

• 1.

  Determining the specifications of network nodes

• • •

    The geographical location of nodes

  • •

    The altitude values of the nodes

  • •

    Based demand of each node

  • •

    Node consumption pattern

• 2.

  Determining the specifications of network pipes

• • •

    Nodes connected to the pipe

  • •

    The length of the pipes

  • •

    Inner diameter

  • •

    Slope of pipes

• 3.

  Determining the consumption pattern

• • •

    Applying pulse demand to network nodes was previously presented. For each node following its location in zones A and B, in addition to the leakage pattern, one of these two patterns were applied based on Equation (3). Demand pulses were measured at 10-minute intervals or 144 consumption coefficients in 24 hours.

The general skeletonized network structure in this village has 100 main pipes. Based on the developed code structure, we first needed to consider an estimated value as the initial pipes' pressure. Here the static pressure of each section of the pipe was assigned. One of the challenges of hydraulic analysis using the MOC is the impossibility of assigning demands to the nodes, after which there is another node. Because by applying the required demand to the hypothetical node J (Figure 11(a)), its unsteady flow analysis is performed according to the same demand. In other words, in this case, the unsteady flow fluctuations resulted from the total demand of pipes numbers 2 and 3 will not be considered. To keep the network structure and analyze steady and unsteady flows, this problem can be solved by creating a new pipe and node according to Figure 11(b). By applying the required demand to the hypothetical node J in pipe number 4, it is possible to analyze the total of unsteady flows fluctuations in these nodes.

Another critical point is that the length of pipe number 4, which also represents the branch pipe, should not be less than the length of the space step of the MOC solution. Therefore, according to the above, the network structure was composed of 166 pipes and 167 nodes. Input data from each pipe include start node number, end node number, pipe length, unsteady wave velocity, inner diameter, slope, and initial pressure and flow rate (in accordance with the MOC (Chaudhry 2014)).

In this research, the control conditions are considered as follows:

• 1.

  The pressure in the critical node and other nodes should not be less than 5 meters (Campisano et al. 2010, 2016).

• 2.

  The pressure in any node should not exceed 50 meters (Guideline 556 of Iran MOE 2017).

• 3.

  Service pressure to customers remains constant range while staying within the standard range (Galuppini et al. 2019).

Hydraulic analysis was done for the skeletonized network of Gavankola village with 166 pipes with a maximum length of 316 m (number of space steps 316). This analysis was performed in a simulation period of 72 min (allocating 30 s for each application of the demand) and considered the urge to create a steady-state condition in the pressure pattern. The solution time step was 0.003184 seconds. The solution matrix had dimensions of 167 × 316×678,240. Pressure values were achieved for all the segments at all times. In the computer solution using the MOC, pressure values were collected at the end of every 30 seconds and for the end section of every 100 main pipes. In a computer solution using the MOC, pressure values were collected in 30-second intervals and for the final segment of all the 100 main pipes.

By examining the pressure behavior of each pipe, five pressure patterns were observed. Figure 12 shows the pipe numbers and the divisions of the pressure isolation zones with letters A, B, C, D and E for the skeletonized network. The number of pipes in these zones was 13, 18, 18, 25, and 26, respectively. To understand the difference in pressure fluctuations in each of the PMAs, the pressure pattern for the pipes with the highest and lowest pressure head in each PMA is presented in Figure 13. For example, the similarity of the pressure fluctuation behavior of zone B is illustrated in Figure 14.

According to the presented topographic map (Figure 8) and the high values of the difference between the start and end of the pipes of this zone, zone A is highly affected by the slope, and the pressure of different periods of this zone is close to its hydrostatic pressure. As zone B is located at the end of the PMA, pressure fluctuations due to changes in consumption pattern completely affect the pipe pressure. The low diameter of the pipes in this zone and the high static pressure (which increases the initial velocity of the water in the pipe) cause such fluctuations. In zones C, D and E, the pipes' pressure fluctuations increase as they gets closer to the end of the PMA. According to the pressure patterns presented in Figure 13, the pressure fluctuations in zones A, B, C, D and E in a 24-hour period were 3.55, 16.3, 4.55, 8.34 and 14 meters, respectively.

The results of steady-state solution by EPANET and unsteady-state flow analysis by MOC were compared. According to Figure 15, the average errors of the results of the MOC with spatial steps of dx = 20 m, dx = 10 m and dx = 5 m are 9.2%, 1% and 0.5%, respectively. In comparison, the average error of EPANET results was 5.6%. This comparison means that the more the pipe is segmented, the more accurate the result. The results of Figure 15(a) present for the pipes with the highest and lowest pressure head in each PMA and show that, by choosing a spatial step dx = 10 m or less, more accurate results can be expected than using EPANET software. Therefore, dynamic hydraulic analysis of the network by MOC provides more accurate results.

Figure 15(b) shows the head deviation between the values ​​calculated using these two methods with a precise value (dx = 2 m) for the entire network pipes. The results ​​are relevant before implementing the PM approach. As can be seen, the amount of head deviation obtained using EPANET software at the initiative of the network (location of distribution reservoir) was less than the precise value and was higher than at the end of the network. However, the values ​​obtained by the MOC method and the spatial step of 5 m differed slightly. The reason for increasing accuracy of the MOC method is the addition of unsteady-state flow coefficients to steady-state coefficients. It is worth mentioning that the numbering of pipes in the skeletonized network structure is such that the number of pipes increases from the beginning to the end of the network. In other words, the place of the reservoir is pipe No. 1 and the end of the network in zone E, is pipe No. 100 (Figure 12).

As mentioned earlier, APC and constant service pressure delivery to customers eliminate excess pressure or low pressure in the network to prevent pressure fluctuations. That will increase customer satisfaction, reduce the number of bursts and reduce network leakage. Figure 16 illustrates a complete image of the hydraulic condition of the network pipes. Accordingly, the average pressure is 63.4 meters and the average pressure fluctuations are 9.7 meters. The highest pressure in the network is in zone B due to the considerable difference in altitude values, which in some places reaches 96 meters. Most pressure fluctuations occur in zones B and E. Variations in consumption patterns as well as being located at the bottom of the network are some of the reasons.

It is important to select the target node to select control strategies. In zone A, an excess pressure of up to about 76 meters is created in a small network district. As this district is separated, the installation of a PRV is not economically justified. Therefore, using other valves (such as flow control valve, fixed outlet PRV, reducing the pipe size, etc.) can solve the problem. As can be seen in zone B, excess pressure occurs in all cases. Due to the proximity of the pressure values, the pressure values of this zone will be corrected by installing a fixed outlet PRV. According to the pressure of other zones, two target nodes are defined:

• 1.

  Minimum pressures at the critical node of pipe No. 97 in zone E.

• 2.

  Maximum pressures at the critical node of pipe No. 56 in zone D.

Here, the end node of pipe No. 97 in zone E is selected as the target node due to having the lowest pressure head for 24 hours (Figure 13). It is worth mentioning that the pressure measuring data related to this node also played a role in determining the pattern of demand pulses (Figure 10(b)). According to the definition of control conditions, it is necessary to adjust the outlet pressure of the RTC-PRV so that the service pressure in this node is equal to 5 meters at any time. According to Figure 13, the minimum and maximum pressure head delivered to the customers in this node occurs in time steps 52 (4:30) and 120 (22:10), respectively. According to the computer code results, the outlet pressure of the RTC-PRV for both minimum and maximum pressure was set at 37.7 m and 46.5 m, respectively. However, it is worth noting that the whole network's minimum and maximum input flows were different from these hours.

The diagrams in Figure 17 show the pressure variations in the pipes with the highest and lowest heads in each PMA. The reason is that supplying service pressure in the desired range (5 to 50 meters) in these pipes will indicate the supply of pressure in the whole PMA. As can be seen, if the 5 m pressure is set for pipe No. 97, it is impossible to remove the excess pressure in many network districts. Therefore, it can be concluded that, in this network, in addition to using the RTC-PRV, it will be necessary to use a VSP.

Using APC strategy, the network was once again solved by considering the service pressure regulation at the end of pipe No. 56 in zone D at head equal to 50 meters. Then, in zones with service pressure less than 5 meters, a VSP installation strategy was adopted. Therefore, the control target node is the end of pipe No. 87 in zone E in this study.

Moreover, according to Figure 16, the minimum pressure in zone B in different states was 56.8, which is more than the defined range. Therefore, in the MOC model, one can expect a minimum pressure of 9.7 meters and a maximum pressure of 29 meters by adjusting the outlet of the PRV with a constant value of 10 meters. It is worth mentioning that PAT can be used instead of PRV here.

The location of the APC equipment was applied in Figure 12. Finally, by using the RTC-PRV strategy and one PRV in zone B, it will be essential to apply two VSPs in a district of the network of zone D and E by modifying a head pressure of 20 and 22 meters. As a result, we will be able to keep the whole network in the standard pressure range. According to Figure 18(a), with APC implementation, all districts of the network will reach standard pressure. The average pressure will drop by 54%, from 63.4 meters to 29.3 meters. That will reduce the number of bursts because of the high pressure of WDN and reduce leakage flow according to Equation (3) due to a reduction in the mean pressure of the network.

Furthermore, another advantage of using the APC method and delivering constant service pressure to the customers will be reducing the pressure fluctuations, reducing the number of bursts (Figure 18(b)). Even if the pressure is within the standard range, the pressure fluctuations can cause bursts. According to the results, the mean pressure fluctuations in the network have decreased by 64%, from 9.7 meters to 3.5 meters. In addition, in the pipes closer to pipe No. 87, the pressure fluctuations are reduced. It is worth mentioning that, in zone B, a fixed outlet valve pressure (i.e. not being affected by upstream pressure changes) will also reduce pressure fluctuations.

This study was aimed to determine the APC strategy and provide constant service pressure of the Gavankola rural network in northern Iran. Due to topographic conditions, this rural network has both excess pressure and low pressure. In addition, high water losses and a large number of bursts had caused customers' dissatisfaction. Dynamic hydraulic analysis of the network was performed by the MOC using various demand pulses. The effective coefficients in the simulation were obtained by performing a laboratory test in the experimental network of the NIT. Field measurements included measuring the target node pressure every 5 minutes and measuring the inlet flow rate of the network every 10 minutes. The most significant achievements of this research include the following:

1. Creating PMAs based on pressure patterns effectively takes the most appropriate action to apply pressure control methods.

2. Reducing pressure fluctuations while implementing PM methods should be one of the research objectives.

3. Here, the average pressure fluctuations in the network pipes decreased from 9.7 m to 3.5 m. That will increase customer satisfaction and will also reduce the number of bursts.

4. Adjusting the outlet of the PRVs by keeping the pressure constant in the control node will reduce the pressure fluctuations near the target node. Whereas in a control state with fixed outlet PRV, the least pressure fluctuations will occur near the valve.

5. The network studied in this research can be a suitable case for other investigations because this network has a variable topography and has unavoidable U-shaped passages.

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