An Iterative Algorithm for Optimizing Pipe Diameter in Pressurized System

Davor STIPANIĆ, Vanja TRAVAŠ*, Lado KRANJČEVIĆ, Danko HOLJEVIĆ 14pt Abstract: In order to achieve optimal flow velocity in stationary flow condition in pressurized pipe systems, it is necessary to determine an optimal value of pipe diameter. It should be noted that even for a simple pipe network with pipes connected in only few numbers of loops, the determination of the optimal pipe diameter is an unintuitive problem and from the mathematical point of view a nonlinear problem that requires some iterative process with a relatively complicated update procedure between iterations. Usually, the optimization algorithms are based on stochastic methods. A deterministic approach is proposed that satisfies the prescribed constrains in terms of the specified flow velocity in the pressurized system of pipes. The presented iterative algorithm for optimizing the pipe diameter is implemented in computer code PIPENET3D (written in FORTRAN90) and used for hydraulic analysis of steady flow in an arbitrary pressured pipe system.


INTRODUCTION
For given steady flow conditions in an arbitrary pressurized pipe system with various pipes diameter and output flow rate of the pressurized system, it is practically impossible to achieve optimal values of flow velocity vj by manually selecting the pipe diameter D j , where the index j denotes a particular geometrical or kinematical property of a pipe associated to a count number j. Hence, the need for developing a computer algorithm which would independently find the optimal pipes diameter D is obvious. Namely, by an iterative process, the required algorithm should properly adjust the pipe diameter and consequently the flow velocity v which will lead to a reduction of operative cost by reducing the pressure head loss Δh in the pressurized system.
Currently, there are various optimization methods that can be used for solving complex engineering problems such as the problem of optimizing pipes diameter in pressurized pipe system. From the mathematical point of view, optimizing pipes diameter in interconnected pressurized pipe systems i.e. looped networks systems represents a so-called NP-hard problem as pipe flows velocities are unknown variables and the time to find a solution grows exponentially with problem size (number of pipes).
First method used for optimizing pressurized pipe system was linear programming (subgroup of convex programming). Alperovits and Shamir [1] suggested successive linear programming gradient method. Development of linear method approach continued but still had same basic problem. Namely, it uses linear functions for solving nonlinear problems and consequently does not provide accurate solution. On the other hand, dynamic programming approach requires too much computer storage capacity as it requires making huge amount of combinations. As a logical development, after the development of a linear approach, nonlinear methods were used for solving this problem. Accordingly, Lansey and Mays proposed [2] the generalized reduced gradient method combined with a network simulator while Varma et al. suggested successive quadratic programming [3].
Similar to a linear method approach, the development of this approach continued but still had its limits. Basically, in nonlinear programming finding the global optima hugely depends on the initial solution. In some cases, if initial solution is too far from the global optima, it cannot escape the local optima. All these attempts have led to the development of methods inspired by natural evolution phenomenon known as Darwin's principle of the survival of the fittest. These kinds of algorithms are classified as naturally inspired meta-heuristics. Among them many have been used to find optimal design system such as Genetic Algorithm [4], Tabu Search [5], Particle Swarm Optimization [6], Ant Colony Optimization [7], Simulated Annealing [8] and Harmony Search [9].
The newest optimization methods used for optimizing water network system are Genetic Heritage Evolution by Stochastic Transmission [10], NLP-Differential Evolution algorithm [11], Hybrid Particle Swarm Optimization, Differential Evolution [12] and Charged System Search algorithm [13]. It should be noted that the aim of all these studies has been to minimize the network investment cost. Basically, the goal is to find the minimum required cost to a certain network topology while satisfying boundary conditions which reduce a problem of constraints for velocity range, pipe size restrictions and pressure at demand nodes.
To minimize the pressure head loss Δh, an iterative optimization algorithm that will accordingly adjust the pipe diameter Dj is hereafter presented and proposed. The algorithm is implemented in a previously written computer code PIPENET3D written in FORTRAN90 and primarily developed for the steady flow analysis and afterward expanded to handle unsteady flow conditions [14]. As the main condition that a pipe diameter D has to satisfy, in order to be considered as optimal D opt , is to achieve optimal flow velocity. Furthermore, a low flow rate Q contributes to the process of depositions while high one is considered economically unjustified. Accordingly, higher flow rates Q lead to high energy losses which grow with mean square velocity. Optimal water velocity value should be chosen before starting computation, and for that condition algorithm will find the optimal pipe diameters. In this respect, this paper assumes that 1 m/s is the optimum water velocity value v within the circular pipes [15].

PRELIMINARY CONSIDERATIONS
Since the actual physical processes of the flow are a complex matter, some assumptions had to be made. In PIPENET3D major pressure head losses Δh are calculated by the Darcy-Weisbach equation [16] and the minor i.e. local losses are computed by common coefficients ξ related to geometry variation in the flow (valves, pipe diameter reductions, pipe diameter expansion, etc.). As usual, the dominant components of velocity vector are parallel to the pipe axis and therefore, the dimensional reduction of the flow field was implemented, and the original threedimensional flow field was reduced to one dimension. Consequently, the flow rate Q is calculated as average velocity v of the local velocity profile. The computational model is based on continuity equation written for incompressible fluid and on Bernoulli's equation for steady flow.

Kirchhoff's Laws
The governing equations for steady flow analysis are given by Kirchhoff. Generally, there are two Kirchhoff's laws. The First law can be interpreted as the principle of conservation of flow rate Q. The algebraic sum of all flow into and out of the junction must be equal to zero. Mathematically, this law can be written as: where Q j represents the volumetric flow rate in pipe j(L 3 /T) and q consumption in junction i(L 3 /T). The index j represents a pipe counter which goes from 1 to n, where n denotes the total number of pipes connected to junction i.
The Second law defines the algebraic sum of all pressure head loss Δh along pipes inside a closed loop. This sum should be equal to zero by performing the summation in respect to the assumed and obtained direction of the flow in pipes of the loop. In other words, this requirement is essential for avoiding different pressure values at the same junction. From the mathematical point of view, the requirement can be defined as: where n is total number of pipes which forms a loop. The pressure head loss Δh j for a pipe j can be expressed by using the Darcy-Weisbach equation: where L j denotes the pipe j length, D j denotes the pipe diameter and λ j denotes the Darcy friction factor which for turbulent pipe flow can be computed from the Haaland explicit approximation [17] of the implicit Colebrook-White equation developed for a full-flowing circular pipe (which is attractive for computational implementation).

Pipe Network Topology
In order to recognize the domain of application of the presented algorithm, it is appropriate to mention that the topology of pipe networks can be generally divided in two main categories: (i) branched or tree-like configurations and (ii) looped or closed configurations. Namely, as would be hereafter presented, the optimization of pipe diameter for the first type of pipe network is trivial. On the other hand, as previously mentioned, the second type needs some iterative procedure. However, it should be noted that water distribution networks are mainly designed as a combination of these two types which could be considered as a third category and in that case a selective procedure should be used. The selection of a particular type of pipe network depends on available funds as well as the required security of consumer supply.

Branched or Tree-Like Configuration
For a steady flow analysis in a branched network the flow direction in pipes is unambiguously identified. Accordingly, the flow occurs from the reservoirs or the water sources through the pipes to consumption in junctions. Thus, this case of configuration is significantly simpler than the looped one. A generic type of such pipe network is illustrated in Fig. 1. Similarly to flow direction, the flow amount is unambiguously defined by the water consumption. Therefore, the value of mean velocity is simply identified by the formula: v j = Q j / A j for each pipe j, where A j denotes the area of the flow of the same pipe. Consequently, the ideal pipe diameter, which follows from the same expression, is trivially determined by adjusting the diameter D j . However, it should be noted that this type of pipe network configuration is avoided for several reasons. Firstly, for this type of configuration the pressure loss is usually higher than the looped network defined for the same ground area. Furthermore, water staleness is more common and in the case of water network failure, flow to consumers behind the breakdown area is interrupted. Moreover, this configuration type is more susceptible to the adverse effects of water hammer. Therefore, notwithstanding its simplicity of pipe diameter optimization, it is more convenient to produce looped network configuration.

Looped or Closed Configuration
Looped or closed configuration type of pipe networks is a configuration where every junction I is connected to at least two neighboring junctions. A generic type of such pipe network is illustrated in Fig. 2. Consequently, if one or more pipelines are closed for repair, water can still reach its consumers by the alternate route [18]. Furthermore, looped network has a significantly better adjustment to the water consumption oscillations and to mitigation of the negative effects of water hammer [19].

Figure 2 Looped or closed configuration of pipes connections in a pressurized system
However, because the looped networks require more pipes installation, they consequently require higher initial financial investment. Furthermore, the flow distribution is much more undefined, because by changing consumption demand, not only amount of flow, but also its direction can be changed [18]. Besides, the flow distribution will not only depend on the network topology and consumption demand, but also on pipe diameter, which is an unintuitive attribute. Hence, changing just one pipe diameter affects flow distribution in the entire network. Consequently, the optimal pipe diameter determination is a very complex problem and impossible feat without the computer algorithms specifically developed for this task. As mentioned previously, even analytical analysis without optimization is impractical for this type of pipe network.
The basic equations for such an analysis are given by Kirchhoff. For a looped type of pipe network Eq. (2) and Eq. (3) form a system of nonlinear algebraic equations and to solve them different methods are used (e.g. Hardy Cross Method [20], Newton-Raphson Method [21] or Linear Theory Method [22]). The program code PIPENET3D is based on the Newton-Raphson Method.

PIPE DIAMETER OPTIMIZATION
The main criterion used for optimization is the value of flow velocity. Pipe diameter is considered as optimal if the value of velocity is equal to 1 m/s. However, developed algorithm, implemented in PIPENET3D, allows it to be manually selected, as well as the fluid properties. Likewise, it can be selected if steady flow analysis is carried out with or without diameter optimization. Pipe diameter optimization for branched or tree-like configuration, as well as looped configuration will be hereafter described.

Pipe Diameter Optimization for Branched or Tree-Like Configuration
For branched or tree-like configuration optimal pipe diameter can simply be found by using the formula for circular pipes: where D opt represents optimal pipe diameter, v opt represents the optimal velocity i.e. the velocity specified by the user and Q represents the flow rate in pipe. Since v opt is the same for every pipe, for water network distribution systems, it can be noticed that the optimal pipe diameter depends only on flow rate Q in each pipe. On the other hand, the flow rate Q depends only on consumption and is constant for any pipe diameter D. Under those circumstances, and for a generic case, the velocity dependence on pipe diameter can be seen in Fig. 3. As seen in the diagram, optimal velocity is uniquely determined and there can be only one diameter for which optimal velocity is achieved (Fig. 3). For looped or closed configuration, the solution is far from trivial.

Pipe Diameter Optimization for Looped or Closed Configuration
For networks with looped or closed configuration, it is not possible to solve the optimization problem analytically. Therefore, in PIPENET3D this nonlinear problem was solved by using the Newton-Raphson Method. The initial idea for optimizing the pipe diameter was to implement Eq. (4) after the network analysis was done for the given flow conditions. This computational procedure defines one iteration cycle. Those new values of pipe diameters would then be used in next analysis i.e. iteration. This iterative procedure would repeat itself until optimal values are found. However, there were two main problems in this simple idea. Firstly, if the pipe diameter in a looped system changes, so does the flow rate of each pipe inside a loop and consequently its velocities. Secondly, diametervelocity dependence diagram is different for looped configuration networks (Fig. 5).
To develop a control mechanism for the adoption of pipe diameters between the iteration, and overcome the mentioned problems, an extensive parametric analysis for a variety of flow situations was performed to study the behavior of the pipe diameter-velocity dependence between iterations. For all the analysis the hereafter conclusion was the same. Namely, for an arbitrary pressurized system, during the course of optimization i.e. iteration, it was noticed that for each pipe inside a loop the flow velocity decreased by reducing the pipe diameter. Depending on the initial kinematical and geometrical properties of the flow, this behavior was observed once the pipe diameter becomes smaller than a particular value (denoted with Dt) which numerically varies depending on the problem. The same behavior was detected regardless of the kinematical and geometrical properties of the flow. That was contrary to expectation and it was shown that the main reason for this unusual feature was flow variability observed during the optimization process conducted for any pipe in a loop. The observed behavior is illustrated in Fig. 4, where the flow rate Q remains constant within the iterations after increasing the pipe diameter above the value Dt. During iteration procedure it was also registered that reducing the pipe diameter below this value firstly gradually decreases the flow velocity and then starts to decrease it rapidly. As seen on a diameter-flow rate dependence diagram (Fig. 4) the flow variability is present between the iterations when the pipe diameter is changed but still smaller than the value D t . On the other hand, for relatively large pipe diameters, greater than D t , the flow rate remains constant during the iterations and adaptations of D, even by increasing the pipe diameter.
The tricky part of optimization was reduced to the question "how to use Eq. (4) in variable flow zone (Fig. 4) i.e. below the value D t ". Namely, if there was relatively small pipe velocity, the trivial application of Eq. (4) in the optimizing algorithm would reduce pipe diameter in order to increase its velocity which will result exactly in the opposite scenarios i.e. the flow velocity will decrease. In other word, the application of Eq. (4) in the in variable flow zone (Fig. 4) would case reverse effect from the expected.
To address this problem, a new general remark is drawn from extensive parametric analysis carried out for different pipe networks. Namely, by plotting the dependence between the velocity in a pipe in an arbitrary closed loop and the change in pipe diameters between iterations, a characteristic function always appears. For a generic case the function is illustrated on Fig. 5 where vopt denotes the prescribed i.e. desiderated flow velocity in a particular pipe. By increasing the pipe diameter above some value of D, the velocity decreases and finally obtains a constant value that does not change by further increasing the pipe diameter D. This is complementary to the observed behavior in Fig. 4, where the flow rate does not increase by increasing D above Dt in the next iteration. On the other hand, note that velocity rapidly decreases if in the next iteration pipe diameter decreases above some threshold value.

Figure 5 Optimal pipe diameter determination for looped or closed networks
The most important conclusion from Fig. 5 is dissimilarity to the branched network and illustrated by the fact that in the looped network two different pipe diameters exist for which optimal flow velocity v opt could be achieved. The question which arises is whether the optimal pipe diameter is smaller or larger. As shown in Fig. 5, the larger pipe diameter is considered as optimal (the nonoptimal diameter D nopt and the optimal one D opt ). There are two main reasons for such a conclusion. Firstly, the smaller pipe diameter is in a zone of highly changeable flow distribution, and thus more sensitive to variable consumer needs. In other words, by just partially changing a consumer demand, the pipe velocity would significantly increase or decrease. Secondly, for a larger pipe diameter there are smaller operating costs (lower energy losses).
It is important to note that the conclusion retrieved form Fig. 5 indicates that there is a maximal value of flow velocity that can be achieved. This is important to note in a context of an optimal velocity that is specified greater than the pick of the function. In that case, the optimal pipe diameter is one from which maximal pipe velocity vmax is achieved (Fig. 6) since it is closest to the optimal velocity. However, this value can be still smaller than the specified one for that pipe, therefore some other corrections were required.
Nevertheless, the conviction that the optimal pipe velocity could be achieved for each pipe brought us to a conclusion that some additional adjustments could have been made. In that sense, it was essential to understand why the maximal pipe velocity in the particular pipe is smaller than the specified optimal velocity of 1 m/s. It was noticed that the main cause for that was small flow rate through the considered pipe, and thus additional algorithm for increasing the flow rate between iterations was conducted.

Figure 6
Optimal pipe diameter determination for looped networks where optimal pipe velocity could not be achieved In order to obtain the desired flow velocity v opt , an additional algorithm is developed and implemented in PIPENET3D. In order to increase velocity in considered pipe from v max to v opt , it is necessary to increase the flow rate through that pipe. Therefore, the algorithm increases the hydraulic head gradient ΔH between the junctions that connects the considered pipe. The graphical explanation is given in Fig. 7 where the velocity in the pipe between hydraulic head H 1 and H 2 reached the value v max instead of the prescribed value v opt . In order to increase the flow rate in the next iteration, and consequently increase the velocity toward the value v opt , the adequate changes in the related pipe diameters should be done. For the given example, which can be analogically adapted to any situation, the arrows on the right side of the notations indicate if the value in considerations must be increased or decreased. To increase the difference of hydraulic head (H 1 − H 2 ), and consequently pipe flow Q 5 , it is necessary to increase the value H 1 and/or decrease the value H 2 . This is achieved by carefully adjusting the neighboring pipe diameters (D 1 , D 2 , D 3 and D 4 ). By increasing the pipe diameter D 1 , the hydraulic head loss is decreased and thus the hydraulic head value of H 1 is increased. By the same logic, decreasing the pipe diameter D 3 , increases the hydraulic head loss in that pipe and consequently decreases hydraulic head value of H 2 . How and why changing the pipe diameter D 2 and D 4 affects the pipe flow distribution inside pipe 5 is seemingly a more complex problem. This was tested and the conclusion is as follows. By decreasing the pipe diameter D 2 , the flow inside it decreases also. Since the flow rate inside the junction is preserved (which follows from the first Kirchhoff law or continuity equation), the flow inside the pipe 5 increases. Similarly, by increasing the pipe diameter D 4 , the flow inside pipe 4 increases also. In this way, by applying the small adjustment to the neighboring pipe diameters, significant increase in hydraulic head for a given pipe could be achieved. Therefore, the flow rate is increased, hence optimal pipe velocity can be reached.

THE ITERATIVE ALGORITHM
The algorithm process necessary to conduct the optimization, based on the previous conclusions is presented hereafter. It should be noted that each flow distribution calculation represents one repetition. Hence, this procedure is iterative (since the problem is nonlinear).

Computational Step #1
Before starting the proposed optimization algorithm, the pipe diameters should be initially assumed so the flow network distribution could be calculated. After that, rough optimization should be carried out which is done by implementing Eq. (4) i.e. according to the algorithm given in Fig. 8. Once new pipe diameters are calculated they become new inputs for flow rate distribution calculation.

Computational Step #2
After rough optimization is done, it is necessary to conduct more precise one. This should be done by optimizing the pipe diameter for which pipe velocity deviates the most from the optimal velocity (maximal |v − v opt |). This step refers to finding the worst optimized pipe diameter according to the algorithm given in Fig. 9.

Computational Step #3
In order to conduct the diameter optimization for the worst optimized pipe, it is necessary to know diametervelocity dependence diagram (Fig. 5) for that particular pipe. Therefore, quadratic equation for that pipe should be found. For this purpose, three points on diameter-velocity dependence diagram should be known of which two already are. Those two points are (0, 0) and (Dw, v w ). Index w like in algorithm 2 ( Fig. 9) represents the worst optimized pipe. Consequently, one additional point should be found in order to define the quadratic equation. In this case, the new diameter for the worst optimized pipe is assumed and once again flow network distribution and pipe velocities are calculated. Now, all three points are known and the quadratic equation v = aD 2 + bD + c can be defined by defining the coefficients a, b and c. In case optimal velocity from this quadratic equation can be obtained then two solutions are possible (since quadratic equation has two solutions), from which the larger one is considered as optimal (Fig. 5). On the other hand, if the optimal flow velocity v opt cannot be obtained from quadratic equation, optimal solution is the maximum of the function, as seen in Fig. 6 (highest possible velocity). The resulting algorithm is illustrated in Fig. 10.

Figure 11
Algorithm for a computational step #4

Computational Step #4
If the same pipe proves to be the worst optimized for more than three times in the row, neighboring pipes diameters should be corrected in order to increase the flow rate through that pipe (Fig. 6). This principle is previously described and shown in Fig. 7. Topological matrix TOP mentioned in the algorithm illustrated in Fig. 11 defines which two neighboring junctions form a pipe. Matrix rows represent the counting number of a pipe while columns represent first and second junction of the same pipe.

Computational Step #5
The penultimate algorithm step is to return to step 2 (after which step 3, 4 and 5 are repeated) until the difference |v − vopt| becomes less than required precision. If for some reason the required precision is a too small number (e.g. zero) this step would repeat itself endlessly. Hence, there should be a counter included which counts how many times this step is repeated. After some finite number of computational procedures is reached, the algorithm should skip the computation and go to the next step (this is needed only if the tolerance value is too small). For practical purposes it is recommended that the tolerance value should be greater that 10 −4 .

Computational Step #6
Final step is to compute the new adopted pipe diameters that will converge to the specified velocity in pipes. The resulting algorithm is illustrated in Fig. 12.  It should be noted that the algorithm presented in Fig.  12 is called in the computation procedure at the end of the algorithm presented in Fig. 11 and uses some array quantity that is defined according to the topological structure of the analyzed network. Namely, to adopt the presented algorithm to a generic pipe network, like every other algorithm with the same purpose, the interconnection data of junctions and pipes should be known in advance. Accordingly, the array idJUCTIONonJUNCTION defines for each junction the junctions that are attached at the considered one. Similarly, the array idPIPEonJUNCTION defines the counting number of pipes attached to a particular junction. Apart from the mentioned array in the algorithm 5 (Fig. 12), a correction factor denoted as cor is introduced to regulate the change in diameter of pipes connected to junctions of the pipe for which the optimization of the diameter is conducted. The performed numerical experiments indicate that the correction factor cor equal to 1/15 leads to a result with reasonable accuracy in a reasonable time of computations. However, if a more accurate solution is demanded, the correction factor can be reduced and consequently the computation time will increase.

NUMERICAL EXAMPLE
The numerical example will be conducted for a pipe system with geometrical and topological data illustrated in Fig. 13. For the considered numerical example the sections lengths are given as: Δx = 1000 m, Δy = 1000 m, where Δx and Δy define the orthogonal distance between junctions. The density of fluid is ρ = 1000 kg/m 3 and dynamic viscosity μ = 0.001307 Pa·s (usual value for water at temperature of 10ºC). It is interesting to examine the optimization for the same consumption load in all junctions, mainly to see how gradually pipe diameters decrease from water reservoir to the final junction. Hence, given consumptions in all junctions are q = 60 l/s. The only exception is junction 1, which represents the reservoir. Sections diameters are all 200 mm. The hydraulic head on junction 1 is set to H1 = 50 m. Absolute roughness is assumed to be ε = 0.01 mm for all pipes. As a result of numerical computation, the distribution of pipe flow rate Q, velocity distribution, hydraulic heads H in junctions and optimal pipe diameters D opt are obtained (so that the prescribed velocity v opt should be obtained). The results of the achieved optimal pipe diameters D opt and velocities for those diameters v are given in Fig. 14 and Fig.  15, respectively. In Fig. 14 the velocity values are given in m/s up to the third decimal. Pipe diameters are given in mm with the same precision. It is worth mentioning that the calculation precision goes up to eleventh decimal (if necessary, it could be carried out with greater precision).
As seen in Fig. 14, the achieved velocity distribution is in the range of 0.943 to 1.061 m/s and pipe diameters from 190 to 972 mm. The only condition for conducting the optimization analysis was achieving the optimal flow velocity, thus the values of achieved pipe diameters are not commercially available and are given as a result of the computational procedure. Usually commercial pipe diameters go from DN 100 mm, DN 150 mm, DN 200 mm and so on. If needed, those values could be implemented as additional term in algorithm. Moreover, some additional constrains can be incorporated in the presented procedure to force i.e. round obtained pipe diameters from iteration to iteration according to the commercial pipe diameters. However, in that case the convergence to a prescribed velocity vopt is not guaranteed.

Figure 15
Computed optimal pipe diameters.
The presented numerical example is characterized by 25 junctions, 40 pipes and 16 loops were given. The achieved range of velocity distribution was from 0.943 to 1.061 m/s. However, it should be noted that by increasing the topological complexity and the number of connected loops, that range is additionally increased, but at the same time it can be controlled by the correction factor cor specified in the algorithm 5 (Fig. 12). Accordingly, for a tested network with 45 junctions, 76 pipes and 32 loops range of velocity distribution was from 0.871 to 1.121 m/s.

CONCLUSIONS
For looped configuration systems it is more than enough to change the diameter of just one pipe to alter the flow distribution inside the pressurized pipe network.
Hence, optimization is a nonlinear problem and only solvable iteratively. For example, if the correction of pipe diameters is conducted based on the performed steady analysis, it has been shown that reduction of pipe diameter could reduce the water velocity through that pipe in the next steady flow analyses. It is an unintuitive feature which proved to be of great importance and should be considered. To resolve this problem in a deterministic manner, an iterative procedure was developed and here presented.
The iterative procedure is conducted by defining the pressure distribution for initial guess of pipe diameters and then by invoking the presented algorithm where the corrections of pipe diameters are performed according to the obtained result of the first iteration. The next iteration cycle is conducted for the new pipe diameters suggested by the presented algorithm. The iterative process is conducted until the variation of all pipe diameters in two neighboring iterations is smaller than the prescribed tolerance.
It is important to note that the presented algorithm is relatively simple and so it can be implemented in any computational procedure used for the computation of pressured distribution in pressurized pipe systems. Although the presented procedure is here elaborated in a context of water pipe system distribution, for which the velocity of 1 m/s can be interpreted as optimal, the presented iterative algorithm can easily handle different scenarios i.e. situations in which the demand for the flow velocity can be different. The iterative procedure will identify the pipe diameters needed to obtain the prescribed velocity in the pipes.