Pipeline Network Modeling and Simulation for Intelligent Monitoring

Ind. Eng. Chem. Res. 1998, 37, 1033-1044

1033

Pipeline Network Modeling and Simulation for Intelligent Monitoring and Control: A Case Study of a Municipal Water Supply System Weerapong Kritpiphat,† Paitoon Tontiwachwuthikul,*,† and Christine W. Chan‡ Process Systems Laboratory, Industrial Systems Engineering, Faculty of Engineering, and Department of Computer Science, Faculty of Science, University of Regina, Regina, Saskatchewan, Canada S4S 0A2

This paper reports on an extended-period simulation for predictions of the flows and pressures of a municipal water supply system. The program has been developed within the scope of the Real-Time Intelligent Systems Project at the University of Regina, Saskatchewan, Canada. The aim of this project is to develop a real-time decision support system for monitoring and control of the water pipeline network operations. The system supports the operator’s decision-making in maintaining the distribution pressures of water at the pumping stations at all times. The system consists of three modules: an expert system module, a water demand predictions module, and a pipeline network simulation module. In this paper, we describe the simulation module which has been implemented using the Hardy-Cross and Newton-Raphson methods. The programs were implemented in MathCad and FORTRAN. A comparison of the two implementations is presented, and some results from the simulations are described. 1. Introduction Pipeline network simulation is an essential tool for control and operations of municipal water supply and distribution systems because it can be used to simulate and analyze network behaviors under different operating conditions. This analysis helps to ensure that the pressure and flow conditions are satisfactory to consumers. One of the challenges of pipeline network operations is how the operational procedures can be adjusted to meet the dynamic and future demands of customers (Brindle et al., 1993). Walski et al. (1990) suggested that analyses of pressures and flows in pipeline networks are needed whenever significant changes in patterns and magnitudes of demands or supplies occur. In the absence of such analyses, the operational procedures may not be optimal, resulting in unnecessarily high operating costs. To overcome these problems, pipeline network simulations can be used to predict the behaviors of pressures and flows throughout the systems. Quevedo et al. (1987) developed an interactive simulation system using linear theory and Hardy-Cross methods. The system has been employed as a decision support aid for the control of water distribution operations in Barcelona, Spain. Even though the system consists of an on-line advisory message system and an anomalous-input detection mechanism, a user is expected to know how to analyze results from the simulation. This is a weakness in the system that can be eradicated if it can automatically translate the numerical results to a recommendation to the user by means of an expert (or intelligent) system. An interactive hybrid system consisting of a pipeline network simulation and an expert system can be used for simulation * Author to whom correspondence should be addressed. Tel: (306) 585-4726. Fax: (306) 585-4855. E-mail: wee@cs. uregina.ca, [email protected], and [email protected]. † Process Systems Laboratory, Industrial Systems Engineering. ‡ Department of Computer Science.

and analysis of the network operations. Such a system would support the control operator in his or her decisions on switching on or off any of the elements in the network, determining the set points for all the valves, or modifying the network data interactively before deciding a control action for the real network. The simulation enables the operator to know in advance the results of his/her action. The system thus supports the operator’s decisionmaking in maintaining the distribution pressures of water at the pumping stations at all times. The system to be described here works in the following fashion. The intelligent system for monitoring and control of the water supply system consists of three modules which implement different artificial intelligence (AI) techniques and mathematical problemsolving methods: (1) expert system module, (2) water demand predictions module, and (3) pipeline network simulation module. Human expertise and heuristics on operations formed the basis of the conceptual or mental model which are represented as IF-THEN rules in the expert system (Chan et al., 1994; Kritpiphat et al., 1995). Modules which employ the artificial neural network approach, fuzzy set, case-based reasoning and machine-learning techniques are developed for water demand predictions (Lertpalangsunti et al., 1995; An et al., 1995). The expert system is responsible for detecting the faults and recommending a series of adjustments on equipment at the current state of operations, while the modules for water demand predictions provide an estimation of water usage in the future. The recommendations from the expert system and the water demand predictions are then input into an extended-period simulation program for verification purposes. An operational procedure of pumps and valves suggested by the expert system will be satisfactory if and only if (1) the pressure at the outlet of each pumping station is within the set-point limit, (2) the water level remains within the range between the minimum and maximum allowance of the reservoirs,

S0888-5885(97)00424-7 CCC: $15.00 © 1998 American Chemical Society Published on Web 01/24/1998

1034 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998

and (3) the combination of pumps and valves does not accelerate any reservoir to overflow or deplete within a short period of time. The simulation module will simulate the dynamic behavior of flows, pressures, and water levels; the results generated are then compared with the maximum and minimum allowance provided by the expert system. In this paper, we present in detail the simulation component of the intelligent system.

3. Problem-Solving Methods of Steady-State and Extended-Period Simulations

2. Pipeline Network Models

where hf represents a drop of pressure in the pipe and K denotes a resistance constant of the pipe. The pressure drop is usually related to flow in the pipe by the characteristic power n. Darcy-Weisbach and Hazen-Williams methods are two methods commonly used to define the values of n and K (Bhave, 1991). For the Darcy-Weisbach method, n is set at 2 and

The primary objective of the simulation study is to build a pipeline network model that can be used to simulate the dynamic behavior of the water supply system based on different pumps and valves combinations. The water level in each reservoir and pressure at the outlet of each pumping station are the key parameters of interest for the extended-period simulations. Before any pipeline network system is simulated, a model of the system must be formulated. Depending on the objective of the study, a water pipeline system can be modeled as a steady-state model, an extendedperiod model, a pipe-sizing model, a water quality model, and a fire-flow model (Walski et al., 1990). Since our focus is on timely simulation of the pipelines, the extended-period model is adopted and will be the subject of this paper. We also describe the fundamental concept of a steady-state model which is used as the basis of the problem-solving technique of the extended-period simulation. 2.1. Steady-State Model. A steady-state model is a typical model used in many fluid mechanics problems. It is assumed that we take a snapshot at one point in time during which conditions of the network are stable (Walski, 1994). When given pressure and flow points from the model, we can use continuity (mass balance) and energy balance equations to estimate flow and pressure at any point in the system. Several numerical problem-solving methods have been developed to solve the steady-state model, for example, Hardy-Cross, linear theory, and Newton-Raphson methods. Algorithms of these methods are available in Streeter and Wylie (1979), Hodge (1990), and Bhave (1991). Wood and Rayes (1981) examined several different techniques for obtaining a solution to pipeline networks and also provided a good comparison of the different methods. 2.2. Extended-Period Model. In reality, water demand is always changing with time. The transient model is used to forecast behaviors in terms of flow, pressure, and inventory supply of the system over an extended period of time. The equations in the transient model are functions of time and flowing distance, expressed as a complicated set of differential equations. To solve the set of equations representing complex pipeline network systems, it is necessary to adopt some assumptions. Generally, in water supply and distribution systems, it is assumed that the transient model can be represented by connecting a series of snapshots of the steady state model. The snapshots are developed based on a steady-state that has been achieved for a particular set of input conditions. Results from the current snapshot will then be fed as the new input conditions for the next steady-state simulation. When the input conditions are varied, a new steady state is determined. This process is repeated until an extended period is constructed. The model is referred to as an extended-period model (Walski et al., 1990).

3.1. Hardy-Cross Method for Steady-State Simulations. Let us start from the characteristic equation of a pipeline as follows:

hf ) KQn

K ) 8f(L/gcπ2D5)

(1)

(2)

where f, L, D, and gc denote a Darcy-Weisbach friction factor, a length of pipe, a diameter of pipe, and a forcemass conversion factor, respectively. For the Hazen-Williams method, the value of n in eq 1 depends on the type of liquid in the system. For water, n is set at 1.852 and

K)

k1L 1.852 4.8704

C

D

(3)

where k1 is a constant (10.675 for SI units) and C is referred to as the Hazen-Williams coefficient. The Hazen-Williams coefficient can be obtained from field experiments or from the literature (Williams and Hazen, 1933; Giles, 1962; Brater, 1976; Kevin and Mays, 1989). The basic idea of the Hardy-Cross method is that conservation of mass at each node can be established initially (Cross, 1936). This means we must first assume an initial guess of flows in every pipe element before starting the pressure drop calculation. For any pipe in which Q0 is assumed to be the initial flow rate, eq 1 can be estimated using a Taylor series expansion as follows: 2 n-2 hf ) K(Qn0 + nQn-1 0 ∆Q + n(n - 1)Q0 ∆Q + ...) (4)

where

∆Q ) Q - Q0

(5)

Q is the correct flow rate and ∆Q is the flow correction factor. If ∆Q is small compared with Q0, all terms of the series after the second may be omitted so that

hf ) K(Qn0 + nQn-1 0 ∆Q)

(6)

Now for any loop,

∑hf ) ∑KQn ) ∑KQ|Q|n-1 ) ∑KQ0|Q0|n-1 + ∆Q∑KnK|Q0|n-1 ) 0

(7)

in which ∆Q has been taken out of the summation since it is the same for all pipe sections in a loop. The absolute-value notations have been added to account for the direction of summation around a given loop. The sign conventions of flows used here are as follows: “If

Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 1035

∑KQ0|Q0|n-1 - pumps ∑ (A0 + A1Q0 + A2Q0|Q0|) + [8(Ce + Cv)/gcπ2D4]Q0|Q0|)/ ∑ min or losses (∑nK|Q0|n-1 - ∑ (A1 + 2A2Q0) + pumps 2 4 [16(C ∑ e + Cv)/gcπ D ]|Q0|) (11) min or losses

∆Q ) -(∆Z +

Figure 1. Sample network with imaginary links.

when passing around the loop in a clockwise fashion, we are moving in the same direction as the flow in a pipe, then Q (or Q0) is positive. If we are moving against the flow in a pipe, then Q (or Q0) is negative” (Watters, 1979). Equation 7 can then be solved for ∆Q.

∆Q ) -

∑KQ0|Q0|n-1 ∑nK|Q0|n-1

(8)

Next, ∆Q is applied to update the flow rate in eq 9 for each pipe section in a loop.

Q0,new = Q ) Q0 + ∆Q

(9)

For complex pipeline networks that contain multiple reservoirs, a special artifice must be introduced. Imaginary links are required to form the closed loops between reservoirs (see Figure 1). The links are defined as imaginary pipelines that interconnect each pair of fixed pressure levels; each carries no flow but maintains a constant pressure drop that is usually equal to the difference in geodetic elevation of water in the reservoirs (∆Z). After accounting for ∆Z and entrance and exit effects of reservoirs in the system, eq 8 becomes:

∆Q ) ∆Z +

∑KQ0|Q0|n-1 + min or∑losses(8Ce/gcπ2D4)Q0|Q0|

∑nK|Q0|n-1 + min or∑losses(16Ce/gcπ2D4)|Q0|

where A0 to A2 are pump coefficients and Cv is a loss coefficient of the valve. The pump coefficients can be determined from a pump-head-discharge curve using a multiple linear least-squares regression, while the loss coefficient of the valve can be found from field experiments or from the literature (Hodge, 1990). The additional calculations of reservoirs, pumps, and valves are easily included in a computer program. An arithmetic procedure for solving the pipeline network model was taken from Watters (1979), Streeter and Wylie (1979), and Hodge (1990). The procedure consists of fourteen steps, as shown in Table 1. The first step is for loop identification, while steps 2 and 3 acquire the input data and initial guesses of flows in the pipes. Next, steps 4-9 perform the iterative calculations to find the correct flow rates in the system. After the value of ∆Q satisfies the convergence criterion in step 10, the pressures are determined in steps 11-13 and the procedure ends. 3.2. Newton-Raphson Method for Steady-State Simulations. The Newton-Raphson method is a numerical method that can solve a set of equations simultaneously. The method is particularly convenient for solving differentiable equations when the value of the desired unknown parameters is known approximately (Pipes, 1958). The Newton-Raphson method for multiple variables is developed based on the Taylor’s series by linearizing a set of nonlinear equations to be a set of linear equations so that the model can be solved using a rootfinding method such as Gaussian elimination with a back-substitution technique. Let

f1(Q1, Q2, ..., Qn) ) 0 f2(Q1, Q2, ..., Qn) ) 0

(10) where Ce is a minor loss coefficient caused by a change in cross-sectional area of flow between pipe and reservoir. For large pipeline network systems, the entrance and exit effects may be negligible when their pressure losses are insignificant compared to the pipe pressure drops. Let us extend eq 10 to include pumps and valves in the pipeline network model. Pumps usually produce an increase in pressure, while flows in the pipes always create a decrease in pressure in the network. Thus, the sign of pumps in a new ∆Q equation must be opposite to the sign of the ΣKQ0|Q0|n-1 term, as shown in eq 11. Details of the derivation can be found in Kritpiphat (1996).

f3(Q1, Q2, ..., Qn) ) 0

(12)

l fm (Q1, Q2, ..., Qn) ) 0 be a set of nonlinear equations with an initial guess of the unknown variables Q1 ) Q1k, Q2 ) Q2k, ..., Qn ) Qnk. Using the Taylor’s series expansion, a first-order approximation can be written as

fi(Q1, Q2, ..., Qn) ) n

fi(Q1k, Q2k, ..., Qnk)

( )| ∂fi

∑ j)1 ∂Q

k j Qj

(QjQkj ) ) 0

for i ) 1, 2, ..., m (13a)

1036 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 Table 1. Arithmetic Procedure of the Hardy-Cross Method Used in This Study 1. Identify all loops in the network. Use only as many loops as are needed to ensure that each pipe is in at least one loop. 2. Input data such as fixed pressure nodes, reservoir levels, system input flows, demand flow rates, equipment status, and all network attributes. 3. Assume initial values of flows (Q0s) that satisfy mass balance equation at each node of the system. Each Q0 value should be accompanied by a positive or negative sign indicating the direction of flow. 4. Compute KQ0|Q0|n-1 and nK|Q0|n-1 for each pipe section. 5. Determine ΣKQ0|Q0|n-1 and ΣnK|Q0|n-1 for each loop. 6. Use eq 11 to determine ∆Q of each loop. 7. Apply the ∆Q correction algebraically to each member of the loop for all the loops. Be very careful of the sign during this correction and be sure to include both ∆Q corrections for the common member of loops. 8. Check the mass balance at each node to confirm a right correction of the sign during ∆Q correction. If a sum of flow at each node is not zero, examine the sign of each ∆Q and repeat step 7. Otherwise, go to next step. 9. Set the value of Q0 to be equal to Q0,new for each pipe section. 10. Repeat steps 4 and 9 until accuracy requirements are met in all loops (i.e., all flow correction factors ( ’s) are arbitrarily small or less than an acceptable level of error). 11. Set the correct flow Q equal to Q0 for each pipe section. 12. Compute hf of each pipe section. 13. Calculate pressure at each node of the system from hf starting from the fixed pressure node (e.g., the distribution outlet). 14. End the algorithm.

[ ][ ]

[ ] [] [ ] Q1k+1 Q2k+1 ‚ ‚ ‚ Qnk+1

This can be represented as a matrix:

∂f1k ∂Q1 ∂f2k ∂Q1 ‚ ‚ ‚ ∂fmk ∂Q1

∂f1k ∂Q2 ∂f2k ∂Q2 ‚ ‚ ‚ ∂fmk ∂Q2

∂f1k ∂Qn ∂f k ‚‚‚ 2 ∂Qn ‚ ‚ ‚ ‚ ‚ ‚ ∂f k ‚‚‚ m ∂Qn ‚‚‚

Q1k Q2k

Q1 Q2 ‚ ‚ ‚ Qn - Qnk

)

n×1

[ ]

m×n

f1(Q1k, Q2k, ..., Qnk) f2(Q1k, Q2k, ..., Qnk) - ‚ ‚ ‚ fm(Q1k, Q2k, ..., Qnk)

so that

[ ] Q1 - Q1k Q2 - Q2k ‚ ‚ ‚ Qn - Qnk

[ ]

∂f1k ∂Q1 ∂f2k ∂Q1 )- ‚ ‚ ‚ n×1 ∂fmk ∂Q1

∂f1k ∂Q2 ∂f2k ∂Q2 ‚ ‚ ‚ ∂fmk ∂Q2

∂f1k ∂Qn ∂f k ‚‚‚ 2 ∂Qn ‚ ‚‚‚ ‚ ‚ ∂f k ‚‚‚ m ∂Qn

(13b)

m×1

‚‚‚

[ ] m×n

f1(Q1k, Q2k, ..., Qnk) f2(Q1k, Q2k, ..., Qnk) ‚ ‚ ‚ fm(Q1k, Q2k, ..., Qnk)

(14)

m×1

and the new corrected values for the next iterative calculation become

Q1k Q2k ) ‚ ‚ ‚ Qnk n×1

Q1 - Q1k Q2 - Q2k + ‚ ‚ ‚ Qn - Qnk n×1

(15)

n×1

A matrix of fi/∂Qj is a Jacobian of all the fi functions. The Jacobian can be supplied either by the user or by a finite difference subroutine. When all the values in the matrix of Qj - Qjk are sufficiently close to zero, the procedure has converged and the iteration is terminated. The most recent value of Qj on the left-hand matrix constitutes the solution to the system. 3.3. Algorithm Used for Extended-Period Simulations. The extended-period model is based on a “steady-state-time-windows” concept, and consists of a series of the steady-state model (Walski et al., 1990). Although operating conditions such as water demand are changing with time, they are assumed to be constant within a time interval. With an initial estimated set of values for the water levels in the reservoirs, the steadystate model can be used to determine flows and pressures of individual pipelines at the first time interval. Then the results of flow can be used to calculate the water levels of reservoirs for the second time interval and so on. Using this procedure, temporal results of an extended-period model can be calculated over a period of time using a series of the steady-state simulation and sequential updating of the operating conditions. The main algorithm of the extended-period simulation is summarized in Table 2. It consists of two parts. The first part operates as a snapshot of steady-state simulations, while the second part is responsible for updating water reservoir levels and reassigning the temporal input to the first part. The second part is therefore an outer loop of the first part and includes five main tasks: (1) detect and update the time variable, (2) assign the temporal water demand flow rates into the first part, (3) feed data specification of that time interval to the first part, (4) receive results from the first part and save them into the output files, and (5) calculate and update the water reservoir levels. 4. Case Study The problem domain is the pipeline network of a municipal water supply system of an actual moderated-


Figure 2. Overall diagram of main pipelines, pumping stations, and reservoirs. Table 2. Algorithm of an Extended-Period Simulation

2nd part

[

1st part

[

1. Input the system attributes, such as the connection index of nodes and links, pipe resistance constants, characteristic constants for each equipment, and geodetic elevation of reservoirs. 2. Give a set of initial conditions for the fluid level in each reservoir. 3. Read a starting time, time interval, and ending time. 4. Set t ) starting time. 5. Assign the fluid level of reservoirs at time t to be equal to the initial value given in step 2. 6. Assign a set of system specifications at time t (e.g., pressure set points, equipment status, system input flows, and demand flow rates). 7. Calculate flows and pressures throughout the network using a Hardy-Cross or Newton-Raphson method. 8. Record flows, pressures, and reservoir levels at time t. 9. Update the time. t ) t + time interval. 10. Compute the fluid level of the reservoirs at new time t. 11. Go to steps 6-10 until t ) ending time. 12. End the algorithm.

sized prairie city in North America. The major sources of water include a lake and a number of underground wells (Karius, 1993). The lake water is treated at the water treatment plant to an acceptable standard of drinking water before being delivered to the system, while underground water is directly drawn from wells located in and adjacent to the city. Water is pumped from its sources to fill in reservoirs located at different elevations in the city and is pumped from the reservoirs to the distribution pipeline, or to another reservoir when it is necessary to adjust water levels (Brindle et al., 1993). Distribution pressures and rates of water transfer can be controlled by means of pumps and valves housed in remote hubs and pumping stations. Figure 2 illustrates an overall diagram of the main supply pipelines, pumping stations, and reservoirs in the problem domain. There are three pumping stations and five reservoirs. The three pumping stations have been named pumping station A, B, and C; and the reservoirs named reservoir no. 1, 2, 3, 4, and 5, respectively. Reservoirs no. 1 and 2 are connected to each other by an equalized pressure line so that the water levels in these reservoirs are always identical. There are three control valves next to reservoir no. 4, whose main function is to control the flow rate and flow direction of water between pumping station B and reservoirs no. 4 and 5. Figure 3 shows a flow diagram of pumps at pumping station A. This pumping station is the main station satisfying about 70% of the water demand of the entire

Figure 3. Flow diagram of pumping station A.

city. There are six pumps in this station, named NP1NP6. NP1-NP5 are fixed-speed pumps, while NP6 is a variable-speed pump. For all fixed-speed pumps, a discharged flow is regulated by a control valve installed in front of the pump. For a variable-speed pump, the flow is controlled by a change of speed of the motor. Pumping station B is a pumping station that supplies water to customers in the downtown and surrounding areas. The amount of water coming out from this station is about 30% of the daily water demand of the city. Figure 4 shows a flow diagram of four pumps at


Figure 4. Flow diagram of pumping station B.

Figure 5. Flow diagram of pumping station C.

pumping station B. Pump FP1 is a two-speed pump that can be operated at either low (79%) or high (100%) speed. FP2 and FP3 are fixed-speed pumps, and FP4 is a variable-speed pump. Pumping station C is a supplementary station. It is used when (1) a large volume of water is required in case of fires, for example, and (2) water demand is extremely high during a summer drought. There are three identical pumps at pumping station C as shown in Figure 5. They are fixed-speed pumps driven by electrical motors which can be run arbitrarily in single mode, in parallel of two, or in parallel of three pumps. 5. Implementation Procedure The implementation procedure consists of the following five steps: (1) collect appropriate data, (2) build a model, (3) code the model into a simulation program, (4) test the pipeline network simulation, and (5) verify results and employ the model. In this study, we have completed the first four steps of system implementation and partially verified the results from the model. 5.1. Collection of Appropriate Data. The data items that describe the components of a pipeline network are represented as objects and attributes of the pipeline network system; they are listed as follows: (1) system layout (a) maps showing the connections between individual reservoirs, pipes, pumps, and valves. (b) geodetic elevation at each node of the system

(2) pipes (a) diameter (b) length (c) roughness, friction factor, or Hazen-Williams coefficient (C factor) (3) pumps (a) pump type (i.e., fixed-speed or variable-speed pump) (b) pump characteristics of head vs flow (pump curve) (c) minimum and maximum flow delivered (d) minimum and maximum speed of motor for a variable speed pump (4) control valves (a) Cv factor (b) calibration curve of flow vs percentage of opening valve (if any) (5) reservoirs (a) shape (e.g., cylindrical, conical, or rectangular) (b) diameter or cross section area of reservoirs (c) geodetic elevation of the reservoir (d) minimum and maximum water level that must be maintained in the reservoir Additional data required to run the extended-period model of the pipeline network include the equipment status, an initial guess of the water level in the reservoirs, the control valve set point, pressure set point, and temporal water demand at each pumping station. Water demands at individual pumping stations vary as a function of time. The future water demands can be obtained using one of the three following methods. First, they can be entered directly by an operator as a temporal table of demand vs time. Second, the operator can estimate water demands from historical operation records using artificial intelligence techniques such as fuzzy logic, neural network, machine learning, and/or case-based reasoning. Third, water demands can be obtained based on numerical estimation using nonlinear equations built into the simulation program. 5.2. Building the Pipeline Network Model. A pipeline network model was developed on the basis of the data collected as outlined in the previous section. First, unnecessary components were eliminated in order to reduce the computational complexity of the model. The remaining components form the skeleton structure that represents the real system in the model. Then, all components were converted into either a link or node, as commonly used in graph theory (Mah, 1990). A link represents a pipe element, such as a pipe, pump, or valve. A node represents a junction between links and can also denote a fixed-pressure reservoir and an outlet of a pumping station where water flows in and out from the system (Streeter and Wylie, 1979). The model was developed based on the following assumptions: (1) All manual valves were omitted from the model. (2) The pipes between any two nodes were represented by a set of pipe length, roughness, and diameter. If two or more sections of pipe exist between any two nodes, all pipe sections will be standardized to be only one section of pipe using an equivalent length method. (3) The initial water level is known for all reservoirs. (4) Reservoirs are assumed to be a fixed-pressure node at each time interval. (5) The output nodes of pumping stations are assumed to be a fixed-pressure node at every time interval. The


Figure 6. A 28-nodes-34-links model.

distribution pressures at these nodes are equal to the pressure set points. (6) The water demand flow rate is assumed to hold constant at each time interval. (7) The flows of water coming from wells and the water treatment plant are assumed to be known and remain constant at every time interval of the simulation if no overriding information is provided by the user. The simulation was to be connected with an expert system that controls and monitors the water supply system. When the expert system recommends changes in flows from wells, changes in pump and valve combinations, or changes in pressure set points, the user will have an option to terminate its current calculations and reinitiate a new extended-period simulation of the model. The new simulation run will use the recommended changes from the expert system and the most updated results from the previous simulation as its starting point, thereby making the model responsive to variations in the water supply operations. Figure 6 illustrates a pipeline network model of the water supply system. The model consists of 28 nodes and 34 links. Any pumps that are identical were grouped together to simplify the model. For example, the three pumps at pumping station C are represented as one equivalent pump which can be run either in parallel or single mode. Although reservoirs no. 1 and 2 differ in size, their water level is identical because of an automatic pressure adjustment of the equalized pressure line connected between the reservoirs. Thus, reservoirs no. 1 and 2 were also grouped together in the model. The model can be further simplified to involve only 14 nodes and 14 links, as shown in Figure 7. While the simplified model will lose some flexibility such as pump option and the control of water flow direction, it still gives a reasonable estimate of water level in the reservoirs vs time. Since the water level is the output

information needed for feeding into the expert system, the simplified model is sufficient to satisfy the objective of the study. After completing a schematic diagram of the pipeline network, a set of mathematical equations representing the pipeline network model can be formed using the node-and-loop method (Stephenson, 1976; Edalat et al., 1987). A loop is defined as a series of pipes forming a closed path. Within the cyclic loops of pipe networks, the equations governing flow rates and pressure drops can be derived by analogy to Kirchoff’s laws for electrical circuits described as follows: Kirchoff’s First Law. The algebraic sum of flows (Q) at each node must be zero. This is identical to the continuity law of mass balance.

∑Q ) 0

(16)

∑Qin - ∑Qout ) 0

(17)

or

For example, Kirchoff’s first law yields the following equation when it is applied to node 10 in Figure 6.

Q26to10 + Q22to10 - Q10to19 - Q10to24 ) 0 (18) Kirchoff’s Second Law. The algebraic sum of pressure differences around each cyclic path must be zero. Kirchoff’s second law is also called the loop law. This law when applied to a sample loop of nodes 1019-20-21-22 in Figure 6 yields:

(P10 - P19) + (P19 - P20) + (P20 - P21) + (P21 - P22) + (P22 - P10) ) 0 (19) Each term of pressure difference (Pj - Pk) forms a relationship with the fluid flow, as given in the energy


Figure 7. Simplified pipeline network model.

balance equation between point j and point k.

(Pj - Pk) ) (static energy losses + kinetic velocity losses + energy losses in pipes + minor energy losses from the entrance and exit effects, valves, and fittings energy gain from pumps) × (conversion factor from energy to pressure) (20a)

(

( ) ∑ ( )

Pj - Pk ) (Zk - Zj) +

8

2

4

Q2 +

gcπ D 8(Ce + Cv)

KQn + ∑ pipes min or losses

2

4

gcπ D

Q2 -

)

(A0 + A1Q + A2Q2) ∑ pumps

γgc g

(20b)

where Zj and Zk are the geodetic elevations at points j and k. γ denotes the specific gravity of water and g represents a gravity constant (9.81 m/s2 or 32.2 ft/s2). After applying Kirchoff’s laws, a set of equations are formed which can then be solved by any systematic procedure, such as Hardy-Cross and Newton-Raphson methods. 5.3. Coding the Model. 5.3.1. MathCad Program. MathCad 4.0 (trademark of MathSoft Inc., Cambridge, MA), which runs on the PC computer platform under the Windows environment, was initially used to implement the model. MathCad was selected because it is a self-documenting software that allows inputs to the model to be programmed in a readable format of mathematical-like equations. It is also simple to use and consists of a number of built-in functions, thereby eliminating tedious coding of subroutines. The model was implemented using a multivariable Newton-Raphson method to solve a simplified network of 14 nodes and 14 links. The program uses a finite difference approximation to estimate the Jacobian of the

Figure 8. Part of the program in MathCad.

functions. Both the Newton-Raphson method and finite difference are included in the built-in function of “Minerr,” which is based on the algorithm given by Levenberg and Marquardt (MathSoft, 1993). The model was conveniently implemented in MathCad because its mathematical-like programming, builtin functions, and documentation utilities are simple to use. Figure 8 presents part of the simulation program in MathCad and illustrates the user-friendly and selfdocumenting nature of the language. The program took about 12 min to run a 10-h simulation with a 1-h time interval on a 486 66-MHz computer. The simulation time is reduced to approximately 5 min when a Pentium 90-MHz computer is used.


Figure 9. Example of the Visual Basic interface menu used for easy communication between the user and the FORTRAN program.

5.3.2. FORTRAN Program. A FORTRAN program was also developed for solving the extended-period pipeline network model. The program is portable to any operating system provided that a FORTRAN 77 (or higher) compiler is available on the system. The program was modified and extended from the “steadystate-model” source code, given in Streeter and Wylie (1979). The program can handle pipe friction factors either in terms of the Hazen-Williams coefficient (C factor) or pipe roughness. Either SI, English, or custom (mL/day and meter) units can be used by proper specification of input data. The program can analyze problems with multiple reservoirs and multiple booster or supply pumps. Not only can the program solve horizontal pipeline problems, but it can also handle different elevation pipeline and reservoir problems. Input of equipment status (on/off status) can be read into the program to prune impossible loops and accelerate the calculations. The program can handle pumps as fixed-speed vs variable-speed pumps, single vs parallel pumps, and electrical-motor-driven vs diesel-engine-driven pumps. It accepts raw data on pump characteristics as well as cubic-polynomial pump coefficients. The FORTRAN program was developed using the algorithm for extended-period simulations given in section 3.3. The program usually converges and yields solutions within a few seconds, but it is not userfriendly. A graphical user interface was implemented which allows the user to view and alter inputs and generates graphs automatically from the output data file. The interface was developed using Microsoft Visual Basic 4.0 on the Windows and Windows/NT platforms. Figure 9 shows an example of the interface menu used in the simulation program. 5.4. Testing and Verification of the Pipeline Network Simulation. The computational results from the MathCad and FORTRAN programs are comparable, but their efficiencies differ. Although it is simple to program using MathCad, the program took approximately 5 min to run a 10-h simulation with a 1-h time

Figure 10. Change of water level in the reservoirs over a period of time.

interval on a Pentium 90-MHz computer. In contrast, the Hardy-Cross and Newton-Raphson methods are both working effectively in the FORTRAN program which takes 3-15 s to run the 28-nodes-and-34-links model and 1-10 s for the smaller 14-nodes-and-14-links model. When the inputs are the normal pumps and valves combinations used in day-to-day operations, the outputs of the simulations are similar regardless of the model and method used. However, the larger model has greater flexibility in terms of control option which reflects the reality of the water supply operations. Figures 10 and 11 are sample results from the pipeline network simulations showing a change of water levels and flows over a period of time. The x-axis in Figures 10 and 11 denotes the future time (t) in hours when the simulation program starts predictions at t ) 0. The y-axis in Figure 10 represents the water level in the reservoirs in meters, while the y-axis in Figure 11 represents the flow rate of water in the pipe element (m3/s). The following discussion illustrates how the simulation program can complement the expert system.


Figure 11. Flow of water in the system over a period of time.

From Figure 10, it can be seen that reservoir no. 4 is filled with water from the current level of 5.7 m at t ) 0 to the level of approximately 7.5 m at t ) 2 h. However, the maximum level of this reservoir is set at 7.4 m in the knowledge base of the expert system. This means that the current operating conditions of pumps and valves would cause reservoir no. 4 to overflow within 2 h. Thus, an operator should consider an alternative pump and valve combinations to transfer the water from reservoir no. 4 to another reservoir. Figure 10 shows that reservoirs no. 1 and 2 will run dry in 8-9 h since the water level decreases gradually from 8.8 m at t ) 0 to approximately zero at t ) 8.7 h. The water levels at reservoir no. 3 and 5 are also being depleted steadily, but reservoirs no. 1 and 2 are the first two reservoirs to run out of water. This prediction can be fed to an expert system which provides suggestions to the operator to reduce the amount of water being drawn from reservoirs no. 1 and 2, and increase the amount of water being drawn from reservoir no. 4. From Figure 11, it is obvious that the flow from node 4 to node 9 (Q4to9) is always negative. This variable is defined in the model as an outlet flow of reservoir no. 4 (see Figure 7). The negative value reverses flow in the opposite direction; which in turn confirms the previous results of an increase in water level of reservoir no. 4. Other variables in Figure 11 are Q5to6, Q3to2, and Q1to13. These are examples of flows in the system, as shown in Figure 7. The stable condition happens at reservoir no. 5 (Q5to6) where water is being drawn constantly at the rate of approximately 0.32 m3/s across the time horizon. Flow from reservoir no. 3 (Q3to2) also has a relatively stable flow of water out from the reservoir, except the transition period from t ) 5 to t ) 8. The flow increases sharply from 0.13 m3/s at t ) 5 to 0.51 m3/s at t ) 6 and then deviates down to a new stable state of 0.38 m3/s at t ) 8. Obviously, the sharp increase in flow occurs at the same time with a sudden increase in demand at pumping station B. From the knowledge acquired for the expert system, it is known that reservoir no. 3 is the main reservoir that supplies 50-80% of water to pumping station B. Thus, any increase in demand at pumping station B indirectly causes an increase in flow of water being drawn from reservoir no. 3. Moreover, Figure 11 illustrates a large fluctuation of Q1to13 across the time horizon. This can be explained

by the fact that the pump in this pipe element is a variable-speed pump which responds quickly to the change of water demand in the system. In this particular example, the discharged flow of the pump varies greatly depending on the water demand at pumping station A. However, the pump cannot keep up with the demand after t ) 5 when the demand reaches its peak. The results from the simulation program will be fed to the expert system for data interpretation and the expert system will inform the user the future consequences of some control actions. In this case, the expert system will alert the user that the current pump combination at pumping station A can be used only for the next 5 h, and then another pump combination whose capacity is higher than that of the current pump must be used to handle the future water demand. The expert system will search through the knowledge base to determine the best alternative of pump combinations that can be used after t ) 5 and recommend it to the user. It is clear from Figures 10 and 11 that water levels and flows change with time. These changes are influenced by the change in water demand over a period of time. They also depend on differences in elevation among reservoirs and pumping stations, and operating conditions of the pump and valve equipment. For example, starting a pump may force a change of flow direction and cause one reservoir to deplete faster than the other. Because of the intricate relationships among the components in the system, pipeline network simulations are useful for predicting system behaviors. Setting the time interval or time step is one of the most important considerations in extended-period simulations. A small time step would yield better precision and better pump control, but increases computational time and requires more disk space to collect input and output data. Since the operator expects a time interval of at least 30 min after issuing a control-action command to visually see a change of water level in the reservoirs, a time step of 0.5-2 h is adopted in our simulation program. These time steps are typical for most pipeline network simulations in the water supply and distribution systems. From the experiments, it was found that reservoirs no. 1 and 2 are most sensitive to a change in time intervals. This is explainable by the fact that both reservoirs are located next to the main distribution outlet (pumping station A) which has a high customer demand that varies greatly with time. The more frequently the user provides values for the water demand to the simulation program, the higher the accuracy of water level predictions obtained from the program. Figure 12 compares three individual predictions of water level in reservoirs no. 1 and 2; all the predictions have identical input conditions but differ in the time step used. Figure 12 shows that the 2-h time step is too large, resulting in inaccurate predictions of the water level. The results from the simulations using a 0.5-h time steps closely approximate the actual data points, while those of a 1-h time step are also acceptable for the first 5 h, with a maximum error of less than 5%. Since the number of loop calculations in the 0.5-h time step simulation is double that of the 1-h time step simulation, a time step of 1 h is preferred for short-term predictions of 5 h or less. However, the smaller time step of 0.5 h is necessary for predictions lasting longer than 5 h.


Nomenclature

Figure 12. Effect of time steps on accuracy of simulation results.

6. Conclusions A simulation model of the municipal water supply network system of a moderated-sized prairie city in North America has been described in this study. The model can be used to simulate dynamic behavior of the system based on different pump and valve operation modes. The developed model can predict temporal flow rates and pressures at any location in the system and determine a change of water level in the reservoirs over an extended period of time. We found that the programming tools have a significant impact on the methods used to solve the problem. The FORTRAN program was found to be superior to MathCad in terms of computational and convergence speed. But the FORTRAN program requires a graphical user interface for enhancing communication between the user and the program. Since the computational time requirements for running the 28-nodes-and-34-links vs the 14-nodesand-14-links model in the FORTRAN program do not differ significantly, the former model is used as the working model. This generates results on the dynamic change of water level in the reservoirs, which are then used as input in the expert system developed for the monitoring and control of the water supply system. The simulation may be modified for gas transmission and distribution networks by introducing compressibility into the equations and substituting compressors for pumps. Furthermore, by substituting appropriate electrical terms such as voltage, resistance, and current for pressure, pressure drops, and flow rates, respectively, a similar package can be developed for the transmission and distribution of electricity. Acknowledgment We are grateful for the generous support of Telecommunications Research Laboratories (TRLabs) and the Natural Sciences and Engineering Research Council of Canada (NSERC). The cooperation of personnel at the Municipal Water Engineering and Water Supply Departments of the City of Regina is acknowledged. We also thank Shane Davidson, a former Computer Science undergraduate student at the University of Regina, for completing the graphical user interface of the FORTRAN program.

A0, A1, A2 ) coefficients of pump curve C ) Hazen-Williams coefficient Ce ) entrance or exit loss coefficient of flow between pipe and reservoir (dimensionless) Cv ) loss coefficient of valve and fitting (dimensionless) D ) diameter of pipe (m) f ) Darcy-Weisbach friction factor (dimensionless) fi ) function of pipe section i used in the Newton-Raphson method g ) gravity constant (9.81 m/s2 or 32.2 ft/s2) gc ) force-mass conversion constant (1 for SI units) hf ) a drop of pressure in the pipe (m) k1 ) constant of Hazen-Williams equation (10.675 for SI units) K ) resistance constant of the pipe L ) pipe length (m) n ) characteristic constant of the pipe Pj, Pk ) pressure at nodes j and k (kPa) Q ) liquid flow rate (m3/s) Qitoj ) flow of liquid from node i to node j (m3/s) Q0 ) initial guess of liquid flow rate used in the HardyCross method (m3/s) Qi ) liquid flow rate in pipe section i (m3/s) Qik ) initial guess of Qi used in the Newton-Raphson method (m3/s) Qin ) inlet flow (m3/s) Qout ) outlet flow (m3/s) ∆Q ) flow correction factor used in the Hardy-Cross method (m3/s) t ) time (h) Zj, Zk ) geodetic elevation at nodes j and k (m) ∆Z ) difference in geodetic elevation of water in the reservoirs (m) γ ) specific gravity of a fluid (dimensionless, 1 for water)

Literature Cited An, A.; Shan, N.; Chan, C. W.; Cercone, N. Discovering Rules from Data for Water Demand Prediction. Proceedings of the 1995 International Joint Conference on Artificial Intelligence, Workshop on Machine Learning in Engineering, Montreal, Canada August, 1995; p 187. Bhave, P. R. Analysis of Flow in Water Distribution Networks; Technomic Publishing: Lancaster, PA, 1991. Brater, E. F. Handbook of Hydraulics; McGraw-Hill: New York, 1976. Brindle, D.; Tunison, W.; Martin L. R. Water Distribution with TIRS. Proceedings of SHARE-80/Winter Conference; IBM Corp.: San Francisco, CA, 1993. Chan, C. W.; Cercone, N.; An, A.; Tontiwachwuthikul, P. ObjectOriented Modeling and Simulation of an Expert System for Monitoring and Control of a Water Distribution System. Proceedings of the 1994 International Conference on ObjectOriented Information Systems; London, U.K., December, 1994; British Computer Society: London, 1994; p 130. Cross, H. Analysis of Flow in Networks of Conduits or Conductors. Univ. Ill. Bull. 1936, November, 286. Edalat, M.; Zaman, A. A.; Tooba, H. A New Model for Designing Gas Distribution Networks. Iran. J. Chem. Chem. Eng. 1987, 9, 38. Giles, R. V. Fluid Mechanics and Hydraulics; McGraw-Hill: New York, 1962. Hodge, B. K. Analysis and Design of Energy Systems; Prentice Hall: Englewood Cliffs, NJ, 1990. Karius, R. H. Long-Term Water Utility Study: Report No. 6: Operational Considerations and Report No. 7: Major Facilities; Associated Engineering (Sask.): Regina, Canada, 1993. Kevin, E. L.; Mays, L. W. Network Simulation Models. In Reliability Analysis of Water Distribution Systems; Mays, L. W., Ed.; American Society of Civil Engineers: New York, 1989.

1044 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 Kritpiphat, W. Pipeline Network Simulations of Gas, Oil and Water Distribution Systems; Comprehensive Examination Report (Internal Document); Faculty of Engineering, University of Regina: Regina, Canada, 1996. Kritpiphat, W.; An, A.; Chan, C. W., Tontiwachwuthikul, P.; Cercone, N. Supervisory and Decision-Support System for Intelligent Monitoring and Control of a Pipeline Network. Proceedings of the 1995 International Conference on Intelligent Systems in Process Engineering, Snowmass Village, CO, July, 1995; paper no. 24. Lertpalangsunti, N.; Kritpiphat, W.; Chan, C. W.; Tontiwachwuthikul, P. Applications of Artificial Neural Networks for Demand Predictions of a Pipeline System. Proceedings of the Sixth Saskatchewan Petroleum Conference, the Petroleum Society of CIM, Regina, Canada, October, 1995; Paper CIM 95155. Mah, R. S. H. Chemical Process Structures and Information Flows; Butterworth-Heinemann: Stoneham, 1990. MathSoft, Inc. MathCad 4.0 User’s Guide Windows Version; MathSoft: Cambridge, 1993. Pipes, L. A. Applied Mathematics for Engineers and Physicists; McGraw-Hill: New York, 1958. Quevedo, J.; Cembrano, G.; Casanova, F.; Ros, E. A Contribution to the Interactive Dispatching of Water Distribution Systems. Artificial Intelligence, Expert Systems and Languages in Model-

ing and Simulation; Kulikowski, C. A., Huber, r. M., Ferrate, G. A., Eds.; North-Holland: Amsterdam, The Netherlands, 1987; p 349. Stephenson, D. Pipeline Design for Water Engineers; Elsevier Scientific: Amsterdam, The Netherlands, 1976. Streeter, V. L.; Wylie, E. B. Fluid Mechanics; McGraw-Hill: New York, 1979. Walski, T. M. Pipe Network Modeling: Video Learning System for CYBERNET Software; Haestad Methods: Waterbury, 1994. Walski, T. M.; Gessler, J.; Sjostrom, J. W. Water Distribution Systems: Simulation and Sizing; Lewis Publishers: Chelsea, MI, 1990. Watters, G. Z. Modern Analysis and Control of Unsteady Flow in Pipelines; Ann Arbor Science: Ann Arbor, MI, 1979. Williams, G. S.; Hazen, A. Hydraulic Tables; John Wiley and Sons: New York, 1933. Wood, D. J.; Rayes, A. G. Reliability of Algorithms for Pipe Network Analysis. J. Hydraul. Div., Am. Soc. Civ. Eng. 1981, 107, 1145.

Received for review June 13, 1997 Revised manuscript received December 9, 1997 Accepted December 10, 1997 IE970424A

Pipeline Network Modeling and Simulation for Intelligent Monitoring

Recommend Documents