Loading ...
Global Do...
News & Politics
60
0
Try Now
Log In
Pricing
<p>CHAPTER 4 _________________________________________________________________________ PIPE NETWORK ANALYSIS 4.1 INTRODUCTION This chapter describes the analysis of steady flow in pipe systems. In any analysis problem all of the physical features of the network are known, and the solution process endeavors to determine the discharge in every pipe and the pressure, etc. at every node of the network. Therefore in this chapter the diameters of all pipes, their lengths and their roughnesses are known, as well as where reservoirs, pumps, pressure reduction valves, and other fittings are located. The ways in which these devices influence the hydraulics of the system will be specified. Design problems, on the other hand, try to select (wisely!) the diameters of pipes, the capacities of pumps, the water surfaces elevations of reservoirs, and so on. Thus, a design problem is distinguished from an analysis problem by the choice of the variables that are regarded as unknown. At some risk we dare to generalize by saying that design problems are usually more challenging to solve than are analysis problems, and design problems usually require the simultaneous solution of a larger system of equations than do analysis problems. A thorough understanding of the techniques of analysis for large networks that are composed of known physical features is a prerequisite to the understanding of the design of networks. The design of pipe networks is the focus of Chapter 5 and is not discussed directly in this chapter. The analysis of a pipe network can be one of the more complex mathematical problems that engineers are called upon to solve, particularly if the network is large, as occurs in the water distribution systems of even quite small cities. A significant fraction of the entire set of equations consists of nonlinear equations, and a large number of these equations must be solved simultaneously. Before digital computers were widely used in engineering practice, it simply was not practical to solve such network problems, and consequently many existing water distribution systems have "grown" with time, based primarily on the best professional judgment of engineers, without any thorough or detailed analysis of the pressures and discharges that could exist in the pipes of the network in response to various combinations of demands on the system. The computer has made it possible to solve such large network problems with ease, and as a result many municipalities and water districts have benefited from the results of relatively detailed computer analyses of their systems in recent years. We believe it is important for an engineer to understand what is being accomplished in these computer solutions. To aid engineers in gaining this knowledge, we begin with the basic principles, and the equations that embody them, that interrelate the discharge in each pipe and the pressure at each node of the network. The same few basic principles of fluid mechanics are the foundation of our work on pipe network analyses. These basic principles are (1) conservation of mass, or the continuity principle, (2) the work-energy principle, and (3) the relation between fluid friction and energy dissipation. Chapter 2 has already introduced these principles. The task here is to employ these ideas effectively in describing a large hydraulic system accurately by means of equations, and then to solve these simultaneous equations efficiently. The oldest systematic method for solving the problem of steady flow in a pipe network is the Hardy Cross method, which is itself an early adaptation of the method of moment dis-tribution from structural engineering in 1936. Before the ready availability of digital computers in the late 1960s, this method was prized because it is so well suited for hand computations. Then it became the basis of most early computer software, but because of convergence problems for large systems containing pumps and other appurtenances, it will not be discussed herein. Over the past quarter century the Newton method has proven to be superior in solving the nonlinear equations, and now networks of 2500 pipes or more can be analyzed successfully with a desktop computer. 4.1.1. DEFINING AN APPROPRIATE PIPE SYSTEM The first step in studying pipe systems is to decide what features are important and to retain them in defining the network of pipes. For large water distribution systems some "skeletonization" usually occurs in this process. In other words, not all pipes in the system are included in the analysis. This skeletonization can take on many forms, such as the following: 1. Not all connections to houses are considered as separate nodes or junctions, and all of the distributed demands along one block of a street, or even a small subdivision, may instead be aggregated or lumped at a single node; 2. Only those pipes that carry the water from the supply sources to the areas of demand are included, i.e., only the main transmission system is considered; 3. Only a few pipes and their associated appurtenances are considered; these compo- nents are regarded as vital to the proper operation of the system. Any study of a pipe system may include one or even all of these levels of skeletoniza- tion; the first preliminary study may start with a model of type (3), and subsequent, more complete analyses may proceed back to type (1) as the adequacy of each is verified, or as components are adjusted. After these analyses have been obtained and studied, it may then be desirable to study intensively the network of pipes within a city block, or the pipes within the area of a major water user, such as a large structure or a manufacturing facility. Thus analyses can treat an entire delivery system, which is generally skeletonized, or a more detailed analysis of the piping system or plumbing within a large building, or a golf course, etc. When an analysis of a building's piping system is conducted, the exterior pressures that are supplied by the larger system can be specified with some degree of confidence since the analyses of the larger "delivery" networks provide this information. There are no hard rules that dictate which pipes should be omitted. Such decisions are often left to the professional judgment of the supervising engineer, and sometimes these decisions are called "art," but the insight gained from analyses at different levels of skeletonization often indicate which pipes should be included in the next level of analysis. Computers can now analyze a problem consisting of many more pipes (e.g., several thousand) than the human mind can visualize in detail when deciding which features should be changed to improve the performance and reduce the costs of the system. Another vital part of defining the network problem is to determine which demands should be specified. The demands on an existing system can be obtained from water usage or billing records. Even for existing systems the data are seldom complete in describing how these demands vary during a day, or from day to day. Analyses are usually needed for a range of system demands, from peak hourly demands down to minimal demand periods (e.g., 2-3 a.m.). During above-average demand periods tanks will have their storage volumes partially depleted, but this loss of volume should be recovered when demands are small. Since a water system may be designed for a 50-year life, the specified demands must appropriately reflect future growth and increases (or possibly decreases) in per capita consumption. In the design of a new system, the demands may have to be based on comparisons with similar cities etc. However, if a system is to be designed to deliver known quantities at specified times, then the problem of determining appropriate demands does not exist. So we see that engineering experience, based on sound judgment, is often required in defining the most appropriate piping system problems to analyze. After the analyst has obtained one or several apparently reasonable solutions, the next step is to verify by measurements in the actual system that reasonable agreement exists be- tween the solution to the mathematical problem and the real system. This process is called network verification. If significant disagreements occur, their causes must be identified. Are some valves in the real system unknowingly closed or partly closed; do some major leaks exist in the real system; has the skeletonization process inappropriately excluded some pipes that carry large flows? These and other possibilities should be explored until reasonable agreement does exist. Firms specialize in field flow measurements to verify that the actual pipe system is modeled properly. After analysis has provided solutions to the network problem for various levels of de- mands, non-ideal or simply inadequate performance parameters can be identified. Some indicators of inadequate or poor performance consist of the following (many other possibilities that are peculiar to an individual system do exist): 1. Pressures at some nodes are too low; 2. Pressures are too high at some nodes (If water is pumped, excessively high pressures cost money, owing to larger power consumption than is needed, more frequent pipe ruptures and the premature replacement of facilities.); 3. Discharges are inadequate and/or pressures are too low to meet emergency demands, such as fire fighting; 4. Pumps are not operating near their peak efficiencies; 5. Some water storage facilities are always nearly empty, while others are nearly full or overtopping (Are the tanks under- or over-sized and located at the best elevations? Unless storage facilities perform near their capacities, they represent investments with cost/benefit ratios that are too large.); 6. Pressure reduction valves, or back pressure valves, are inactive or open (Perhaps they are not needed, or pipes should be removed or added.); 7. Too much of the supply is coming from expensive sources, etc. Upon identifying deficiencies, the engineer's next task is to determine the best, most economical means of overcoming these deficiencies and improving the performance of the system. How best to accomplish this will again require sound professional judgment, but sound judgment seldom occurs in the absence of relevant information, i.e., the engineer must understand the system. Section 5.7 of the next chapter discusses sensitivity analyses, which could materially aid this evaluation process. In the following work we will express the head loss in each pipe in a network by an exponential formula hf = KQn, Eq. 2.17, so one presentation of the theory covers all cases, regardless of whether the Darcy-Weisbach equation, the Hazen-Williams equation or the Manning equation is used to express the head loss as a function of discharge. Only the values for K and n change, as we saw in Chapter 2. 4.1.2. BASIC RELATIONS BETWEEN NETWORK ELEMENTS The two basic principles, upon which all network analysis is developed, are (1) the con- servation of mass, or continuity, principle, and (2) the work-energy principle, including the Darcy-Weisbach or Hazen-Williams equation to define the relation between the head loss and the discharge in a pipe. The equations that are developed from the continuity principle will be called Junction Continuity Equations, and those that are based on the work-energy principle will be called Energy Loop Equations. The number of these equations that constitutes a non-redundant system of equations is related directly to fundamental relations between the number of pipes, number of nodes and number of independent loops that occur in branched and looped pipe networks. In defining these relations NP will denote the number of pipes in the network, NJ will denote the number of junctions in the network, and NL will denote the number of loops around which independent equations can be written. In defining junctions, a supply source will not be numbered as a junction. A supply source is a point where the elevation of the energy line, or hydraulic grade line, is established; a junction, or node, is a point where two or more pipes join. A node can exist at each end of a "dead end" pipe; this instance is an exception to the usual rule, where only one pipe is connected to a node. In a branched system there are by definition no loops, and thus NL = 0 for any branched system. In branched systems the number of nodes is always one larger than the number of pipes, or NP = NJ - 1, unless a reservoir is shown at the end of one pipe and this is not considered to be a junction. Then NP = NJ. (This situation actually occurs.) Figures 4.1a and 4.1b depict a small branched network and also a small looped network. [1] (1) [4] [5] [6] [7] [3] [2] (5) (6) (4) (2) (3) [1] (1) [4] [5] [6] [7] [3] [2] (5) (6) (4) (2) (3) [8] [9] (10) (12) (11) (8) (9) (7) Figure 4.1 (a) A small branched system. (b) A small looped system. 6 pipes, 7 nodes 12 pipes, 9 nodes In the branched system the number of nodes is 7 and the number of pipes is 6 (one less than the number of nodes), whereas in the looped system there are 12 pipes and 9 nodes, i.e., the number of nodes is less than the number of pipes. For a looped network the number of loops (around which independent energy equations can be written) is given by NL = NP − NJ (4.1) if the network contains two or more supply sources, or NL = NP − (NJ −1) = NP − NJ +1 (4.2) If the network contains fewer than two supply sources and the flow from the single source is determined by adding all of the other demands, then this source is shown as a negative demand and the source is called a node. We note that this is the case in the small looped network in 4.1.b, so we have NP = 12, NJ = 9 and NL = 12 - (9 - 1) = 4. Equation 4.2 also applies to a branched system with NL = NP - NJ + 1 = 0, since a branched system can have at most one supply source. Actually, every pipe system must have at least one supply source, but sometimes the source is not shown since the discharge from this supply source is known, and the source is replaced by a negative demand, which is a flow coming into this junction, equal to the sum of the other demands. When this is done, the elevation of the energy line (or HGL or pressure) must be specified at a node so the other HGL elevations can be determined. Energy loops that begin at one supply source and end at another are called pseudo loops, i.e., these loops do not close on themselves. The number of pseudo loops, which are numbered as part of NL, equal the number of supply sources minus one. In forming pseudo loops all supply sources must be located at the end of a pseudo loop. It is generally possible to form more loops than are needed to produce a set of independent equations. As each new loop is formed, see that at least one pipe in the new loop is not a part of any prior loop; in this way the formation of redundant loops can usually be avoided. For special devices, such as pressure reduction valves, this rule of experience must be modified slightly, as will be described later. 4.2 EQUATION SYSTEMS FOR STEADY FLOW IN NETWORKS Three different systems of equations can be developed for the solution of network analysis problems. These systems of equations are named after the variables that are regarded as the principal unknowns in that solution method. These systems of equations are called the Q-equations (when the discharges in the pipes of the network are the principal unknowns), the H-equations (when the HGL-elevations, also simply called the heads H, at the nodes are the principal unknowns), and the ∆Q-equations (when corrective discharges, ∆Q, are the principal unknowns). Each of these three systems of equations will be studied separately. 4.2.1. SYSTEM OF Q -EQUATIONS The analysis of flow in pipe networks is based on the continuity and work-energy principles. To satisfy continuity, the volumetric discharge into a junction must equal the volumetric discharge from the junction. Thus at each of the NJ (or NJ - 1) junctions an equation of the form of Eq. 4.3 is obtained: QJ j − ΣQi = 0 (4.3) In this equation QJj is the demand at the junction j, and each Qi is the discharge in one of the pipes that join at junction j. These junction continuity equations are the first portion of the Q-equations. The work-energy principle provides additional equations which must be satisfied. These equations are obtained by summing head losses along both real and pseudo loops to produce independent equations. There are NL of these equations, and they are of the form of Eq. 4.4 or 4.5, depending upon whether the loop is a real loop or a pseudo loop, respectively, and they are the second portion of the Q-equations: Σh fi = 0 (4.4a) Σh fi = ∆WS (4.5a) When the head losses are expressed in terms of the exponential formula, then these equations take the forms Ki ∑ Qi n = 0 (4.4b) Ki ∑ Qi n = ∆WS (4.5b) in which the summation includes the pipes that form the loop. If the direction of the flow should oppose the direction that was assumed when the energy loop equations were written, such that Qi becomes negative, then there are two alternatives: One is to reverse the sign in front of this term, i.e., correct the direction of the flow. The second, which is generally preferred when writing a program to solve these equations, is to rewrite the equations as follows: Ki Qi Qi ∑ n−1= 0 (4.4c) Ki Qi Qi ∑ n−1= ∆WS (4.5c) To illustrate the system of Q-equations, consider the small 5-pipe network shown in Fig. 4.2. Since no supply sources are shown for this network, only NJ - 1 junction continuity equations are available. Writing these continuity equations for nodes 1, 2, and 3 leads to [1] (1) [4] [3] [2] (5) (4) (2) (3) CHW = 120 (all pipes) All elev. = 0' 10" - 35 00' 12" - 1500' 4.45 ft3/s 1.11 ft3/s 3.34 ft3/s 12" - 3000' 6" - 10008" - 1000' Figure 4.2 Small network. F1 = Q1 + Q3 − 4.45 = 0 F2 = − Q1 + Q2 + Q4 +1.11 = 0 F3 = − Q4 − Q5 + 3.34 = 0 (4.6) In these equations and throughout the text, Fi for any number i is any equation which has been arranged into the form Fi = 0; this format is useful for identification purposes and also for subsequent mathematical and numerical developments. The continuity equation at node 4 is - Q3 - Q2 + Q5 = 0. However, this equation is not independent of the other three nodal equations since it is, except for sign, the sum of these three equations. Now let us use the Hazen-Williams equation to define the head loss in each pipe. In expressing these head losses, the exponential equation will be used. From Eq. 2.18 the coefficients are the following: K1 = 2.018, K2 = 5.722, K3 = 19.674, K4 = 4.847, K5 = 1.009 (4.7) The energy equations around the two loops may be written as F1 = K1Q1 1.852 + K2Q2 1.852 − K3Q3 1.852 = 0 F2 = K4Q4 1.852 − K5Q5 1.852 − K2Q2 1.852 = 0 (4.8) or F1 = 2.108Q1 1.852 + 5.722Q2 1.852 −19.674Q3 1.852 = 0 F2 = 4.847Q4 1.852 −1.009Q5 1.852 − 5.722Q2 1.852 = 0 (4.9) which might alternatively be written as follows if the directions of the flows are uncertain: F1 = 2.108Q1 Q1 0.852 + 5.722Q2 Q2 0.852 −19.674Q3 Q3 0.852 = 0 F2 = 4.847Q4 Q4 0.852 −1.009Q5 Q5 0.852 − 5.722Q2 Q2 0.852 = 0 (4.10) These two work-energy equations are obtained by starting at nodes 1 and 2, respectively, and traversing the respective loops I and II in the clockwise direction. If the assumed direction of flow opposes this traverse, a minus sign precedes the head loss term for that pipe. The simultaneous equations, such as those appearing as Eqs. 4.6 and 4.10, are called Q-equations because it is the Q's, the discharges in the pipes, that are the set of primary unknown variables. After the Q's are found (and the head loss in each pipe is therefore also known) for each pipe, the HGL-elevations at the nodes can be found by starting at a known HGL-elevation and repeatedly applying the exponential formula for head loss to each pipe. If the network is a branched system, then the Q-equations consist of only the junction continuity equations. These can be solved, giving the discharge in every pipe, with a linear algebra solver, i.e. a pocket calculator that implements matrix algebra, a spreadsheet, TK- solver, MathCAD or a solver such as SOLVEQ. Thereafter the individual heads are com- puted from the head loss equation for each pipe. Methods for solving looped systems are described later. Example Problem 4.1 The coefficients K and n for the exponential formula are given in the table for each of the three pipes in this branched system. Find the discharge in each pipe and the pressure at each node. The elevation of the HGL at node 1 is H1 = 100 ft. [1] (1) [4] [3] [2] (2) (3) 17' 14' 15' 0.8 ft3/s 0.5 ft3/s 1.2 ft3/s 2.5 ft3/s HGL1 = 100 20' Pipe K n 1 3.772 1.944 2 5.730 1.926 3 16.29 1.889 A formal method for solving the Q-equations for this network is to use matrix algebra to write the coefficient matrix, an unknown vector, and a known vector in the following way: 1 −1 0 0 1 −1 0 0 1 Q1 Q2 Q3 = 0.8 1.2 0.5 The solution is Q1 Q2 Q3 = 2.5 1.7 0.5 We should note that it is easy to obtain this solution by inspection. Starting at the down- stream end of a branch, the analyst can progressively satisfy each junction continuity equa- tion while working upstream. After finding the discharges, the elevations H of the HGL are determined by starting where the HGL is known, in this case at node 1, and computing the head losses in the pipes that join this node; then the frictional head losses hf are sub- tracted from the known values of H, etc. until all of the nodal heads have been determined. The pressures are then determined by subtracting the nodal elevations z from the heads H and multiplying this by the specific weight, i.e., p = γ(H - z). The tables which follow present the computed values for this network: Pipe Q , ft3/s hf = KQ n , ft 1 2.5 22.395 2 1.7 15.922 3 0.5 4.398 Node H down = Hup - hf, ft. Pressure , lb/in2 1 100.0 Given 34.67 2 100.0 - 22.395 = 77.61 27.13 3 77.61 - 15.922 = 61.69 19.37 4 61.69 - 4.398 = 57.29 18.76 * * * Example Problem 4.2 Write the system of Q-equations for this network. In these equations use the parameters Ki and ni, in which i is the pipe number. [1] (1) [4] [3] [2] (5) (4) (2) (3) WS1 ∆Q2 QJ1 QJ3 WS2 QJ4 (6) ∆Q1 I II QJ2 Since two supply sources are present, four junction continuity equations are available. They are the following: F1 = Q1 − Q2 − Q4 − QJ1 = 0 F2 = Q2 − Q3 − QJ2 = 0 F3 = Q3 + Q4 + Q5 − Q6 − QJ3 = 0 F4 = Q6 − QJ4 = 0 The number of energy loop equations is NL = NP - NJ = 6 - 4 = 2 (one pseudo and one real loop). These equations follow: F5 = K2Q2 n2 + K3Q3 n3 − K4Q4 n4 = 0 F6 = K1Q1 n1 + K4Q4 n4 − K5Q5 n5 − WS1 + WS2 = 0 Since F4 requires Q6 = QJ4, this dead end pipe could be removed from the network, and the demand at node 3 would then be changed to QJ3 + QJ4. These steps would reduce NP to 5 and NJ to 3, and they would eliminate any need for equation F4. After the HGL elevation, H3, at node 3 has been determined by solving this equation set, then H4 can be found by computing hf6 and subtracting it from H3. * * * 4.2.2. SYSTEM OF H -EQUATIONS If the elevation of the energy line or hydraulic grade line throughout a network is initially regarded as the primary set of unknown variables, then we develop and solve a system of H-equations. One H-equation is written at each junction (or at NJ - 1 junctions if fewer than two supply sources exist). Since looped pipe networks have fewer junctions than pipes, there will be fewer H-equations than Q-equations. Every equation in this smaller set is nonlinear, however, whereas the junction continuity equations are linear in the system of Q-equations. To develop the system of H-equations, we begin by solving the exponential equation for the discharge in the form Qij = (h f ij / Kij ) 1/nij = [(Hi − H j ) / Kij ] 1/nij (4.11a) Here the frictional head loss has been replaced by the difference in HGL values between the upstream and downstream nodes. In addition, in this equation a double subscript notation has been introduced; the first subscript defines the upstream node of the pipe, and the second subscript defines the downstream node. Thus Qij and Kij denote the discharge and loss coefficient for the pipe from node i to node j. An alternative way of writing Eq. 4.11a is Qk = (h f k / Kk ) 1/nk = [(Hi − H j ) / Kk ] 1/nk (4.11b) in which k is the pipe number. Substituting Eq. 4.11 into the junction continuity equations, Eq. 4.3, yields QJ j − {[(Hi − H j ) / Kij ] 1/nij }in ∑ + {[(H j − Hi ) / Kij ] 1/nij }out ∑ = 0 (4.12a) in which the summations are over all pipes that flow to and from junction j, respectively. If it is desired to automate the choice of sign, then Eq. 4.12a can be written as QJ j − {[(Hi − H j ) / Kij ](Hi − H j ) / Kij 1/nij −1 }in ∑ + {[(H j − Hi ) / Kij ](H j − Hi ) / Kij 1/nij −1 }out ∑ = 0 (4.12b) As an application of the H-equations with HGL-elevations at the nodes as the un- knowns, consider the one-loop network in Fig. 4.3 which consists of three junctions and three pipes. In this network two independent continuity equations are available, and conse- quently the head at one of the junctions must be specified. In this case at node [1] the head loss in the pipe that connects the reservoir to the network can first be determined, and then this value can be subtracted from the reservoir water surface elevation to determine H1. [1] [3] [2] Figure 4.3 A network with three pipes and three junctions. The two continuity equations are Q12 + Q13 = QJ1 = QJ2 + QJ3 Q21 + Q23 = − QJ2 (or − Q12 + Q23 = − QJ2 ) (4.13) Although in the second equation the flow in pipe 1-2 is toward the junction, the discharge Q21 is not preceded by a minus sign since the notation 2-1 takes care of this. Alternatively the equations could have been written at junctions 2 and 3 instead of 1 and 2. Substituting Eq. 4.11 into these continuity equations gives the following two equations to determine H2 and H3: H1 − H2 K12 1/n12 + H1 − H3 K13 1/n13 = QJ2 + QJ3 − H1 − H2 K12 1/n12 + H2 − H3 K23 1/n23 = − QJ2 (4.14) Since a negative value cannot be raised to a power, a minus sign must precede any term in which the subscript notation opposes the direction of flow. Systems of these equations will be called H-equations, since the HGL-elevations are the primary unknowns. After the heads H are found, then the discharge in each pipe can be obtained from Eq. 4.11. Example Problem 4.3 Write the system of H-equations for the network in Example Problem 4.2. Refer to the figure in Example Problem 4.2. Only the junction continuity equations are used in forming the H-equations, and each Qi is replaced by [(Hui − Hdi ) / Ki ] 1/ni , in which subscript u is the upstream node and subscript d is the downstream node. The system is F1 = WS1 − H1 K1 1/n1 − H1 − H2 K2 1/n2 − H1 − H3 K4 1/n4 − QJ1 = 0 F2 = H1 − H2 K2 1/n2 − H2 − H3 K3 1/n3 − QJ2 = 0 F3 = H2 − H3 K3 1/n3 + H1 − H3 K4 1/n4 + WS2 − H3 K5 1/n5 − H3 − H4 K6 1/n6 − QJ3 = 0 F4 = H3 − H4 K6 1/n6 − QJ4 = 0 * * * 4.2.3. SYSTEM OF ∆Q -EQUATIONS The number of ∆Q-equations is normally about half the number of H-equations for a network. This reduction in number is not necessarily an advantage, since all of the equations are nonlinear and may contain many terms. These equations consider the loop corrective discharges or ∆Q's as the primary unknowns. These corrective discharges or ∆Q's will be determined from the energy equations that are written for NL loops in the network, and thus NL of these corrective discharge equations must be developed. To obtain these equations, we replace the discharge in each pipe of the network by an initial discharge, denoted by Qoi, plus the sum of all of the initially unknown corrective discharges that circulate through pipe i, or Qi = Qoi + ∆Qk ∑ (4.15) in which the summation includes all of the corrective discharges passing through pipe i. The initial discharges Qoi must satisfy all of the junction continuity equations. It is not difficult to establish the initial discharge in each pipe so that the junction continuity equa- tions are satisfied. However, these initial discharges usually will not satisfy the energy equations that are written around the loops of the network. Equation 4.15 is based on the fact that any adjustment can be added (accounting for sign) to the initially assumed flow in each pipe in a loop of the network without violating continuity at the junctions. It is very important to understand the validity of this decompo-sition; it may help to note that any ∆Q entering a junction as it flows around a loop must also leave that junction, and vice versa (See Fig. 4.4). Because of this fact, we decide I II Figure 4.4 A two-loop portion of a network. to establish NL energy loop equations around the NL loops of the network, in which each initial discharge plus the sum of corrective loop discharges Σ∆Qk is used as the discharge. The junction continuity equations are satisfied by the initial discharges Qoi and are not a part of the system of equations. The corrective discharges can be chosen as positive if they circulate around a loop in either the clockwise or counterclockwise direction. It is necessary to be consistent within any one loop, but the sign convention may change from loop to loop, if desired. A corrective discharge adds to the flow Qoi in pipe i if it is in the same direction as the pipe flow, and it subtracts from the initial discharge if it is in the opposite direction. To summarize how the ∆Q-equations are obtained, replace the Q's in the energy loop equations, Eqs. 4.4 and 4.5, by Qi = Qoi ± ∆Qk ∑ (4.16) Here the summation includes all corrective discharges which pass through pipe i, and the plus sign is used if the net corrective discharge and pipe flow are in the same direction; otherwise the minus sign is used before the summation. Thus Eqs. 4.4 and 4.5 become Ki ∑ Qoi ± ∆Qk ∑ { }ni = 0 for real loops (4.17a) and Ki ∑ Qoi ± ∆Qk ∑ { }ni = ∆WS for pseudo loops (4.18a) To automate the choice of sign, these equations can be rewritten as Ki ∑ Qoi ± ∆Qk ∑ { }Qoi ± ∆Qk ∑ ni −1 = 0 for real loops (4.17b) and Ki ∑ Qoi ± ∆Qk ∑ { }Qoi ± ∆Qk ∑ ni −1 = ∆WS for pseudo loops (4.18b) To illustrate the development of the system of ∆Q-equations, consider the network in Fig. 4.5. If the Q-equations were used, there would be five junction continuity equations and two loop equations, a total of seven equations. If the H-equations were used, there would be an equation for the HGL-elevation at each of the five nodes where the head is unknown (The head at one node must be known.). But there will be only two [1] (1) [4] [5] [6] [3] [2] (5) (6) (4) (2) (3) (7) 5 ft3 /s 1.1 ft3 /s 3.3 ft3 /s 3.5 ft3/s 0.7 ft3/s K=2.717n=1.945K=0.497n=1.938K=2.722n=1.942K = 4.108 n = 1.921 n = 1.878 K = 1.628 n = 1.917 K = 0.755 n = 1.929 K = 1.793 Figure 4.5 A seven-pipe network. ∆Q-equations, one for each real loop in this network. These two equations are F1 = K1 Qo1 + ∆Q1 ( )n1 + K2 Qo2 + ∆Q1 − ∆Q2 ( )n2 − K3 Qo3 − ∆Q1 ( )n3 − K4 Qo4 − ∆Q1 ( )n4 = 0 F2 = − K5 Qo5 − ∆Q2 ( )n5 + K6 Qo6 + ∆Q2 ( )n6 + K7 Qo7 + ∆Q2 ( )n7 − K2 Qo2 + ∆Q1 − ∆Q2 ( )n2 = 0 (4.19) The next step is to find an initial estimate for the discharge in each pipe that will satisfy all of the junction continuity equations. One possible set of initial discharges is Qo1 = 1.75 ft3/s, Qo2 = 3.55 ft3/s, Qo3 = 1.05 ft3/s, Qo4 = 1.75 ft3/s, Qo5 = 1.8 ft3/s, Qo6 = 1.5 ft3/s, and Qo7 = 0.4 ft3/s. When these initial discharges and the parameters that are listed on the network sketch are substituted into Eqs. 4.19, the following two equations result; their solution will yield values for the two unknowns ∆Q1 and ∆Q2: F1 = 1.793 1.75 + ∆Q1 ( )1.929 + 0.497 3.55 + ∆Q1 − ∆Q2 ( )1.938 − 4.108 1.05 − ∆Q1 ( )1.921 − 2.717 1.75 − ∆Q1 ( )1.945 = 0 F2 = − 0.755 1.8 − ∆Q2 ( )1.917 + 2.722 1.5 + ∆Q2 ( )1.942 +1.628 0.4 + ∆Q2 ( )1.878 − 0.497 3.55 + ∆Q1 − ∆Q2 ( )1.938 = 0 (4.20) Upon obtaining the solution to these two equations for the two unknowns, ∆Q1 and ∆Q2, the discharge in each pipe can easily be determined by adding these loop corrective discharges to the initial discharges. From these discharges the head loss in each pipe can be determined, and from these values the head at each node can be found. The nonlinearities in these systems of equations create difficulties when we seek the solution. Later in the chapter we apply the Newton method to overcome this problem. Example Problem 4.4 Write the system of ∆Q-equations for the network depicted in Example Problem 4.2. The ∆Q-equations are based on the energy loop equations alone. Therefore these equa- tions can be obtained by taking the equations for F5 and F6 directly from Example Problem 4.3 and replacing each discharge Qi by Qoi ± ∆Qk ∑ . The ∆Q-equations are F5 = K2 Qo2 + ∆Q1 ( )n2 + K3 Qo3 + ∆Q1 ( )n3 − K4 Qo4 − ∆Q1 + ∆Q2 ( )n4 = 0 F6 = K1 Qo1 + ∆Q2 ( )n1 + K4 Qo4 − ∆Q1 + ∆Q2 ( )n4 − K5 Qo5 − ∆Q2 ( )n5 − WS1 + WS2 = 0 With the writing of the ∆Q-equations we must also provide values for Qoi that satisfy all of the junction continuity equations. For this purpose we assume that all four demands are equal to 1.0. Then the following values could serve as an acceptable initialization of the discharges: Qo1 = 3, Qo2 = 1, Qo3 = 0, Qo4 = 1, Qo5 = 1, and Qo6 = 1. * * * 4.3 PRESSURE REDUCTION AND BACK PRESSURE VALVES A pressure-reducing valve (PRV) is designed to maintain a constant pressure at its downstream side, independent of the value of the upstream pressure at the device. Only two exceptions to the maintenance of this downstream pressure exist: (1) if the upstream pressure becomes less than the valve setting; or (2) if the downstream pressure exceeds the pressure setting of the valve so that flow would occur in the upstream direction if the PRV were not present. If the first exception occurs, the valve has no effect on flow conditions except to create a local loss; generally its effect is then like a globe valve in dissipating additional head beyond the friction loss in that line. If the second condition occurs, then the PRV acts as a check valve, preventing a reverse flow in the line. Then the PRV allows the pressure immediately downstream from the valve to exceed its pressure setting. In this way PRV's reduce pressures in portions of a pipe distribution system if the pressure would otherwise be excessive, and they may also be used to control the choice of a supply source in response to various demand levels. In the latter applications the PRV acts as a check valve until the pressure is reduced to a critical level by large demands, at which time additional supply sources become available. A back-pressure valve (BPV) is designed to maintain a constant pressure upstream from it, independent of the value of the downstream pressure. Like a PRV there are exceptions to this normal mode of operation. Should the upstream pressure become less than the pressure setting, the valve can not maintain the pressure setting since it is not an energy- creating device, and the most it can do is shut down the flow in its line. Should the flow want to reverse direction from the positive flow direction through the valve, the valve opens completely and acts as a local loss device. A BPV is used in situations where the pressure would otherwise become too low in elevated portions of the network. Such a situation arises, for example, where a pump is needed to sustain adequate pressures in a higher part of a network but is not needed in the lower portions of the network; without a BPV, or possibly several BPV's, the flow pattern might then lead to discharges through pressure relief valves in the lower portions of the network and possibly create excessively large pressures in the lower region. The equations that describe the behavior of a pipe network that contains PRV's or BPV's must include new features to account properly for the effects of these valves on the discharges and pressures throughout the network. Furthermore, the analysis of a pipe network with such devices must be able to determine the correct operational conditions, i.e., determine whether the PRV's and BPV's are operating in their normal modes or in one of their exceptional modes. Methods for adding such devices into network analyses are described in these sections. The discussion begins with pressure reduction valves. Underlying the writing of the three systems of equations described in Section 4.2 is the basic assumption that a relation exists between the magnitude of the discharge in a pipe and the amount of the head loss, or head difference, between the ends of this pipe. Such a relation does not exist if a PRV (or a BPV) is present in the pipe. Therefore a pipe with a PRV in it should not appear in a normal energy loop equation. However, in the usual mode of operation for a PRV a constant head is maintained at its downstream end; in this way it behaves like a reservoir. Furthermore, regardless of its mode of operation the discharge at the upstream node of a pipe containing a PRV will be the same as the discharge at the downstream node of this pipe. The details of developing a proper system of equations to describe a network containing one or more PRV's are different, depending upon whether a system of Q-equations, H-equations, or ∆Q-equations are desired. Therefore, each of these will be described in a separate section. 4.3.1. Q-EQUATIONS FOR NETWORKS WITH PRV'S/BPV'S The procedure for developing the Q-equation system for a network containing PRV's is as follows: (1) write the junction continuity equations in the usual manner, ignoring the PRV's; (2) replace each PRV with an artificial reservoir which has a water surface elevation equal to the HGL-elevation that is the sum of the pressure head set on the PRV and its elevation in the pipeline; finally (3) write the energy equations around the loops of this modified network. The resulting equations describe the normal mode of operation. Let's try this procedure on the seven-pipe network shown in Fig. 4.6, in which a PRV exists in pipe 6, located 500 ft. downstream from node 1, the upstream end of this pipe. Since a PRV is a directional device, we must always identify the upstream and downstream P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) All e = 0.02" 1.0 ft3/s 6" - 200'100' 50' 50' 50' 20' 90' 6" - 1000' HGL1 = 55' PRV1 1" - 1500 ' 6"-1000'6" - 800' 6"-2000'6"-1000'500'1.0 ft3/s HV1 Q ft3/s Head ft. 1.0 60 1.5 55 2.0 48 Figure 4.6 A seven-pipe network. ends of the pipe containing it. The system of Q-equations for this network consists of four junction continuity equations and three energy loop equations. According to the usual rules, an independent junction continuity equation can be written for each of the four junctions since there are two supply sources for this network. These junction continuity equations are F1 = − Q1 + Q2 + Q6 + Q7 = 0 F2 = 1.0 − Q2 − Q3 = 0 F3 = Q3 − Q4 + Q5 − Q7 = 0 F4 = 1.0 − Q5 − Q6 = 0 (4.21) These continuity equations are unaffected by the presence or absence of a PRV in the net- work. We next modify the network so the upstream portion of the pipe containing the PRV is removed and the PRV is replaced by a reservoir with a water surface elevation equal to the HGL of the pressure setting of the PRV (see Fig. 4.7). Of the three loops that exist in this modified network, only one is a real loop which traverses pipes 2, 3, and 7. Two pseudo loops also exist. One pseudo loop connects the two original supply sources. This loop can start at the reservoir and end at the source pump so it includes pipes 4, 7, and 1. The second loop must extend from the artificial reservoir created by the PRV to one of the other supply sources (or another artificial reservoir, if two or more PRV's exist in the network). The shortest path for this second pseudo loop will traverse pipes 4, 5, and 6. In writing the head loss in pipe 6, only that portion of the pipe downstream from the PRV is used. A modified loss coefficient K' will be used to denote this change in the exponential formula. The new coefficient K' equals the K for the pipe containing the PRV, multiplied by the ratio of the pipe length from the PRV to the pipe's downstream end divided by the total pipe length, or K' = K(Ld / L) (4.22) or in this example K6 ' = K6 (500 /1000) = 0.5K6 . P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) All e = 0.02" 1.0 ft3/s 6" - 200'100' 50' 50' 50' 20' 90' 6" - 1000' 1" - 1500 ' 6"-1000'6" - 800' 6"-2000'6"-500'1.0 ft3/s 55' Figure 4.7 The modified seven-pipe network. The energy equations are F5 = K2Q2 n2 − K3Q3 n3 − K7Q7 n7 = 0 (real loop) F6 = K4Q4 n4 − K7Q7 n7 − K1Q1 n1 + hp1 −100 + 90 = 0 ( pseudo loop) F7 = K4Q4 n4 + K5Q5 n5 − K6 ' Q6 n6 + 55 −100 = 0 ( pseudo loop) (4.23) The head produced by the pump hp1 can be defined by a second-order polynomial passing through three points of the pump curve, or hp1 = AQ1 2 + BQ1 + C (4.24) We have now formed seven independent equations that contain the seven unknown discharges Q1, Q2, . . . , Q7. In this example the real loop that was lost by having a PRV in pipe 6 is replaced by an additional pseudo loop. We see that the number of equations again equals the number of unknown discharges. To obtain a solution for this network by using the computer program NETWK, the in- put data can be prepared (see the NETWK user manual for input data requirements or the condensed description of this input on the CD) as listed in Fig. 4.8. The solution tables from NETWK are reproduced in Fig. 4.9. A study of this output will show that the PRV is operating in its normal mode of operation. Example of a network containing a PRV /* 6 1 4 1000 RESER PIPES 7 1 3 1500 1 4 100 1 0 1 1000 6 0.02 NODES PUMPS 2 1 2 1 0 50 1 1 60 1.5 55 2 48 90 3 3 2 800 2 1 VALVE 4 0 3 200 3 0 6 500 55 5 3 4 2000 4 1 20 RUN Figure 4.8 Input data for the network shown in Figs. 4.6 and 4.7. This solution indicates that the PRV is operating in its normal mode of maintaining the set pressure at its downstream end because the reported downstream HGL-elevation equals the value specified in the input data. If this had not been the case, the solution from NETWK would have indicated either that the PRV had shut off the flow in pipe 6 or that it was completely open and replaced by a minor loss. In solving the network equations, if the discharge in pipe 6 had been negative, then the program would have noted that the PRV would act as a check valve, preventing a reverse flow. If this situation should occur, then the network problem would be altered so it would only have six pipes instead of seven (pipe 6 would not exist in this modified network). The equations describing the flows in LOSSES DUE TO FLUID FRICTION IN ALL PIPES POWER LOSS = 11.51 H.P. = 8.585 KWATTS. ENERGY LOSS = 206.0 KWHRS/DAY PUMPS: PIPE HEAD Q HORSEPOWER KILOWATT KWATT-HRS/DAY 1 59.1 1.11 7.43 5.54 133.0 ELEVATION OF HGL UPSTREAM AND DOWNSTREAM OF PRVS: Figure 4.9 Output tables from NETWK. PIPE UPSTREAM DOWNSTREAM HGL HGL 6 121.79 55.00 UNITS OF SOLUTION ARE: DIAMETER, inch LENGTH, feet HEAD, feet ELEVATION, feet PRESSURE, lb/in2 DISCHARGE, ft3/s DARCY-WEISBACH FORMULA USED FOR COMPUTING HEAD LOSSES PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 3 Q VEL . HEAD LOSS HLOSS/ 1 0 0 0 ft. in in ft3/s ft/s ft. 1 0 1 1000 6.0 20.0 1.11 5.65 27.28 27.28 2 1 2 1000 6.0 20.0 1.07 5.43 25.26 25.26 *3 2 3 800 6.0 20.0 0.07 0.34 0.10 0.12 4 0 3 200 6.0 20.0 0.89 4.54 3.55 17.74 5 3 4 2000 6.0 20.0 0.96 4.91 41.47 20.74 6 1 4 1000 6.0 20.0 0.04 0.18 0.04 0.04 7 1 3 1500 1.0 20.0 0.01 1.31 25.56 17.04 AVE. VEL. = 3.19 ft/s, AVE. HL/1000 = 15.46, MAX. VEL. = 5.65 ft/s, MIN. VEL. = 0.18 ft/s *Flow direction is reversed from input data. NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. ft3/s gal/min ft. ft. lb/in2 ft. 1 0.00 0.0 50. 71.81 31.1 121.81 2 1.00 448.8 50. 46.55 20.2 96.55 3 0.00 0.0 50. 46.45 20.1 96.45 4 1.00 448.8 20. 34.98 15.2 54.98 AVE. HEAD = 49.95 ft., AVE. HGL = 92.45 ft. MAX. HEAD = 71.81 ft., MIN HEAD = 34.98 ft. Figure 4.9, concluded. Output tables from NETWK. this modified network would consist of the original equations with the last one omitted. If the HGL elevation at node 1 had been lower than the HGL setting of the PRV, then it would be known that the PRV would not be able to sustain its pressure setting, and the network problem must then be solved by using equations that replace the PRV with a minor loss device. For this last mode of operation the last energy equation would be replaced by a real loop equation traversing pipes 5, 6, and 7. Pipe 6 would contain a minor loss device to represent the PRV as being fully open. The procedure for writing the system of Q-equations should now be apparent for back- pressure valves (BPV's) in networks. As with PRV's, the junction continuity equations are written ignoring the presence of BPV's. The junction continuity equations are unaffected by the existence of a BPV. In writing the energy equations, the upstream side of each BPV is replaced by an artificial reservoir; in each case the pipe segment from the downstream end of the BPV to its downstream node is then removed, and the energy equations are written for this revised network. The writing of a system of Q-equations will be illustrated with the network in Fig. 4.10, which has 9 pipes and 6 nodes, is supplied by a source pump and has two tanks P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) 200 m (8) (9) [6] [5] 150 - 2000 Diameters in mm Lengths in m 0.015 m3/s 0.015 m3/s HGL = 195 m 200 - 2000 800 m BPV 200 - 1500 150m0.02 m3/s 300-1000140 m 250 - 1200 150-100060 m 135 m 70 m 0.03 m3/s 100m150 - 1500 180 m 150 - 1200 80 m 0.02 m3/s 0.02 m3/s 150 - 1000 All e = 0.02 mm Figure 4.10 A 9-pipe, 6-node network. (reservoirs) connected to it. Without a BPV (or some other device) this network would cause the lower reservoir at the end of pipe 9 to overflow. There are six junctions in this network. The corresponding six junction continuity equations are F1 = 0.015 − Q1 − Q2 + Q3 = 0 F2 = 0.020 − Q3 + Q4 + Q8 = 0 F3 = 0.015 − Q4 + Q5 = 0 F4 = 0.020 − Q5 + Q6 − Q9 = 0 F5 = 0.020 − Q7 − Q6 = 0 F6 = 0.030 − Q8 + Q7 = 0 (4.25) Before forming the loops around which the energy equations are written, an artificial reser- voir is placed on the upstream side of the BPV with a water surface elevation equal to the HGL resulting from the pressure setting of the valve. The pipe downstream from the BPV is removed. When these changes are completed, the network appears as in Fig. 4.11, and energy equations can next be written around the loops of this modified network. Three loops are needed, since NL = NP - NJ = 9 - 6 = 3. These are all pseudo loops and may be P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) 200 m (8) (9) [6] [5] 150 - 2000 Diameters in mm Lengths in m 0.015 m3/s 0.015 m3/s 800 m 200 - 1500 150m0.02 m3/s 300-1000140 m 250 - 1200 150-1000135 m 70 m 0.03 m3/s 100m150 - 1500 180 m 150 - 1200 80 m 0.02 m3/s 0.02 m3/s 150 - 1000 All e = 0.02 mm 195 m Figure 4.11 The modified 9-pipe, 6-node network. composed in the following way: the pipes in loop 1 are 1 and 2; the pipes in loop 2 are 1, 3, and 4 (upstream portion); and the pipes in loop 3 are 9, 6, 7, 8, and 4 (upstream portion). It is incorrect to write a loop through pipes 9, 5, and 4 (the upstream portion) because a BPV sets the pressure on its upstream side. Hence the energy equations in the Q-equation system are the following: F7 = K1Q1 n1 − hp1 − K2Q2 n2 −180 + 200 = 0 F8 = K1Q1 n1 − hp1 + K3Q3 n3 + K4 ' Q4 n4 −180 +195 = 0 F9 = K9Q9 n9 + K6Q6 n6 − K7Q7 n7 − K8Q8 n8 + K4 ' Q4 n4 −135 +195 = 0 (4.26) One possible input file to NETWK for the solution of this problem is presented in Fig. 4.12, and the resulting solution tables are presented in Fig. 4.13. Network Containing BPV 1 .015 140 /* 2 .02 140 $SPECIF NFLOW=3,NPGPM=3,NUNIT=4 $END 3 .015 70 PIPES 4 .02 60 1 0 1 1200 250 .02 5 .02 80 2 0 1 2000 150 6 .03 100 3 1 2 1000 300 RESER 4 2 3 2000 2 200 5 3 4 1000 150 9 135 6 4 5 1200 PUMPS 7 6 5 1500 1 .1 35 .15 32 .2 28 180 8 2 6 1500 200 BPVALVE 9 0 4 1000 150 4 1200 195 NODES RUN Figure 4.12 The input data file to NETWK for the 9-pipe, 6-node network. From this solution we see that the BPV dissipates 65.88 m of head to sustain the upstream HGL setting of 195 m. This value is obtained by subtracting the downstream HGL from the BPV setting. It is a worthwhile exercise to begin with the head losses in the PIPE DATA table and verify the HGL elevations reported in the NODE DATA table; it will lead to a better understanding of the BPV and its effect on pressures and discharges in this network as the BPV operates in its normal mode. If the solution had shown a negative flow through pipe 4, then the downstream pressure would actually be larger than the BPV setting, and the valve would open up completely. For this occurrence the BPV must be re-placed by a minor loss device, and then this modified network problem could be studied. If the HGL at node 2 (the node immediately upstream from the BPV) were less than the HGL established by the BPV setting, then the BPV would close completely. The pipe LOSSES DUE TO FLUID FRICTION IN ALL PIPES POWER LOSS = 65.18 H.P. = 48.63 KWATTS. ENERGY LOSS = 1167.0 KWHRS/DAY PUMPS: PIPE HEAD Q HORSEPOWER KILOWATT KWATT-HRS/DAY 1 34.88 10.0 46.9 35.0 839.8 HGL DOWNSTREAM AND UPSTREAM FROM BPV 4 129.12 195.00 Figure 4.13 The output tables from NETWK for the 9-pipe, 6-node network. PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 3 Q VEL. HEAD LOSS HLOSS/ 1 0 0 0 m mm mm m3/s m/s m 1 0 1 1200 250 20.0 0.102 2.09 15.58 12.98 2 0 1 2000 150 20.0 0.004 0.21 0.75 0.37 3 1 2 1000 300. 20.0 0.091 1.29 4.28 4.28 4 2 3 2000 300. 20.0 0.006 0.08 0.06 0.03 5 4 3 1000 150. 20.0 0.009 0.52 1.90 1.89 6 5 4 1200 150. 20.0 0.015 0.86 5.69 4.74 7 6 5 1500 150. 20.0 0.035 2.00 33.11 22.07 8 2 6 1500 200. 20.0 0.065 2.08 25.25 16.83 9 0 4 1000 150. 20.0 0.014 0.79 4.03 4.03 AVE. VEL. = 1.10 m/s, AVE. HL/1000 = 7.47, MAX. VEL. = 2.09 m/s, MIN. VEL. = 0.08 m/s NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. m3/s ft3/s m m kPa m 1 0.015 0.53 140.0 59.25 580.7 199.25 2 0.020 0.71 140.0 55.02 539.2 195.02 3 0.015 0.53 70.0 59.08 579.0 129.08 4 0.020 0.71 60.0 70.97 695.6 130.97 5 0.020 0.71 80.0 56.66 555.3 136.66 6 0.030 1.06 100.0 69.78 683.8 169.78 AVE. HEAD = 61.79 m, AVE. HGL = 160.13 m MAX. HEAD = 70.97 m, MIN. HEAD = 55.02 m Figure 4.13 (Concluded) The output tables from NETWK for the 9-pipe, 6-node network. containing the BPV should then be removed from the network, and the problem could then be solved by using the equations for this modified network; then the BPV could not maintain the pressure setting, and it would simply prevent any flow from passing through the pipe in which it is installed. 4.3.2. H -EQUATIONS FOR NETWORKS WITH PRV'S/BPV'S The procedure for writing the H-equations for a network that contains PRV's and/or BPV's is described here. First, view the HGL resulting from the pressure setting of the de- vice as a reservoir, since under normal operation the HGL is fixed by the device. Second, place an additional unknown variable on the other side of the device to represent the eleva- tion of the HGL there. We will denote this variable by Hvi, in which i is the number of the device. For the first PRV or BPV i = 1, for the second i = 2, etc. For a PRV the value of Hvi is the HGL-elevation immediately upstream from the valve, whereas Hvi is the HGL-elevation immediately downstream from the valve for a BPV. Third, the junction continuity equations are written in the usual way, with the difference between the upstream and downstream HGL-elevations, divided by K for this pipe, all raised to the reciprocal of the discharge exponent n, i.e., Qk = {(Hi − H j ) / Kk } 1/nk . Finally, since an additional unknown is introduced for each PRV or BPV, one additional equation must be added to the system of continuity equations for each device. These additional equations are obtained by noting that the head losses in the upstream and downstream portions of the pipe containing the device are proportional to these two lengths. For a PRV this equation is (HGL − Hd )Lu − (Hu − Hvi )Ld = 0 (4.27) in which Lu and Ld are the lengths upstream and downstream from the device, respec- tively, and Hd and Hu are the HGL-elevations at the downstream and upstream ends of pipe i containing the device. Using again the network depicted previously in Fig. 4.6 to illustrate the formation of the H-equation system, we would first modify the network as shown in Fig. 4.14: P1 [1] (1) [4] [3] [2] (5) (4) (2) (3) (7) All e = 0.02" 1.0 ft3/s 6" - 200' 100' 50' 50' 50' 20' 90' 6" - 1000' 55' 1" - 1500 ' 6"-1000'6" - 800' 6"-2000'500'500'1.0 ft3/s Hv1 (6) Figure 4.14 A 9-pipe, 6-node network containing a BVP, as modified. The final H-equations for this network are the following: − {(90 + hp1 − H1) / K1 } 1/n1 + {(H1 − H2 ) / K2 } 1/n2 + {(H1 − Hv1) / K6 " }1/n6 + {(H1 − H3 ) / K7 } 1/n7 = 0 1.0 − {(H1 − H2 ) / K2 } 1/n2 − {(H3 − H2 ) / K3 } 1/n3 = 0 {(H3 − H2 ) / K3 } 1/n3 − {(100 − H3 ) / K4 } 1/n4 + {(H3 − H4 ) / K5 } 1/n5 − {(H1 − H3 ) / K7 } 1/n7 = 0 1.0 − {(H3 − H4 ) / K5 } 1/n5 − {(HGL1 − H4 ) / K6 ' }1/n6 = 0 (HGL1 − H4 )Lu − (H1 − Hv1)Ld = (HGL1 − H4 ) − (H1 − Hv1) = 0 (4.28) in which K6 " and K6 ' are the coefficients for the upstream and downstream portions of pipe 6, respectively. The network in Fig. 4.10 that contains the BPV should now be viewed as shown in Fig. 4.15, and the H-equation system for this network is presented as Eqs. 4.29. F1 = 0.015 − {(180 + hp1 − H1) / K1 } 1/n1 − {(200 − H1) / K2 } 1/n2 + {(H1 − H2 ) / K3 } 1/n3 = 0 F2 = 0.020 − {(H1 − H2 ) / K3 } 1/n3 + {(H2 − HGL1) / K4 '' }1/n4 + {(H2 − H6 ) / K8 } 1/n8 = 0 F3 = 0.015 − {(Hv1 − H3 ) / K4 ' }1/n4 + {(H3 − H4 ) / K5 } 1/n5 = 0 F4 = 0.020 − {(H3 − H4 ) / K5 } 1/n5 + {(H4 − H5 ) / K6 } 1/n6 − {(135 − H4 ) / K9 } 1/n9 = 0 F5 = 0.020 − {(H4 − H5 ) / K6 } 1/n6 − {(H6 − H5 ) / K7 } 1/n7 = 0 F6 = 0.020 − {(H2 − H6 ) / K8 } 1/n8 + {(H6 − H5 ) / K7 } 1/n7 = 0 F7 = 1200(H2 − HGL1) − 800(Hv1 − H3 ) = 0 (4.29) P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) 200 m (8) (9) [6] [5] 150 - 2000 Diameters in mm Lengths in m 0.015 m3/s 0.015 m3/s 800 m 200 - 1500 15 0 m 0.02 m3/s 300-1000140 m 250 - 1200 150-1000135 m 70 m 0.03 m3/s 100m150 - 1500 180 m 150 - 1200 80 m 0.02 m3/s 0.02 m3/s 150 - 1000 All e = 0.02 mm 195 m 1200 m Hv1 Figure 4.15 A seven-pipe network, modified. 4.3.3. ∆Q-EQUATIONS FOR NETWORKS WITH PRV'S/BPV'S Let us begin this section by reviewing the underlying concept that is used in writing the ∆Q-equations: if the junction continuity equations are satisfied by the initial discharges Qoi, then a corrective loop discharge, ∆Q, can flow around a loop without violating the principle that the discharge into all junctions will still equal the discharge out of these junctions, regardless of the magnitude of ∆Q. These corrective loop discharges can be regarded as the primary unknowns, and the resulting solution to the system of equations will produce discharges that also satisfy the energy equations around the loops. Therefore the discharges Qi in the Q-equation loops were replaced by Qoi ± ∆Qk ∑ . For the junction continuity equations to remain valid for any values of ∆Qk, these corrective loop discharges must circulate around loops that are formed before any PRV's or BPV's are converted into artificial reservoirs. Thus it is necessary to consider two sets of loops with this method. The first set is the set of loops around which the ∆Q's circulate, and the second set is the set of loops that is used in writing the energy equations. These two sets of loops will be called the corrective discharge or ∆Q loops, and the energy loops. The ∆Q loops are formed while ignoring the existence of PRV's or BPV's. These devices are later replaced by artificial reservoirs, and the energy e quations are written for this modified network. Thus the energy set of loops will always contain more pseudo loops than does the ∆Q set of loops by the number of PRV's and/or BPV's that exist in the network. To track these two separate sets of loops in figures, the ∆Q loops will list ∆Qi by the arc denoting the loop, and the energy loops will be numbered by roman numerals I, II, etc. To illustrate the writing of the ∆Q-equations, we examine again the network with a PRV that is in Fig. 4.6. This network is redrawn below in Fig. 4.16 to display both the corrective discharge loops and the energy loops. To emphasize that ∆Q loops are different than energy loops, ∆Q3 is chosen to pass through pipes 4, 3, 2, and 1. A more effi-cient route for this corrective discharge loop would traverse pipes 4, 7, and 1, coinciding with energy loop II, because one less pipe is in this loop. To obtain the ∆Q-equations, we replace each Qi in the energy equation portion of the Q-equations by Qi = Qoi ± ∆Qk ∑ , in which the ∆Q's must be those circulating through pipe i, as defined by the ∆Q loops, and the sign before each term in the sum is determined by whether ∆Q agrees with or opposes the direction of the assumed discharge Qoi. If the directions agree, the sign is positive; otherwise the sign is negative. The resulting ∆Q-equations for this network are listed as Eqs. 4.30. P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) All e = 0.02" 1.0 ft3/s 6" - 200' 100' 50' 50' 50' 20' 90' 6" - 1000' HGL1 = 55' PRV1 1" - 1500 ' 6"-1000'6" - 800' 6"-2000'6"-1000'500'1.0 ft3/s III ∆Q 2 ∆Q3∆Q 1 II I Figure 4.16 The network in Fig. 4.6, modified for solution with the ∆Q-equations. F1 = K2 Qo2 + ∆Q1 − ∆Q3 ( )n2 − K3 Qo3 − ∆Q1 + ∆Q3 ( )n3 − K7 Qo7 − ∆Q1 + ∆Q2 ( )n7 = 0 F2 = K4 Qo4 + ∆Q3 ( )n4 − K7 Qo7 − ∆Q1 + ∆Q2 ( )n7 − K1 Qo1 − ∆Q3 ( )n1 + hp1 −10 = 0 F3 = K4 Qo4 + ∆Q3 ( )n4 + K5 Qo5 + ∆Q2 ( )n5 − K6 ' Qo6 − ∆Q2 ( )n6 − 45 = 0 (4.30) These equations are in a sense incomplete until each Qoi is replaced by a value. When this step is completed, then these three equations contain only three unknowns, ∆Q1, ∆Q2 and ∆Q3. In principle each reader could produce a different set of acceptable values for the initial discharges, so long as they do indeed satisfy each and every junction continuity equation. One valid set of Qoi's is Qo1 = 1.0, Qo2 = 1.0, Qo3 = 0.0, Qo4 = 1.0, Qo5 = 1.0, Qo6 = 0.0, and Qo7 = 0.0. Finally, the pump head now becomes hp1 = A Qo1 + ∆Q1 − ∆Q3 ( )2 + B Qo1 + ∆Q1 − ∆Q3 ( ) + C when the pump curve is fitted with a second-order polynomial. If desired, as an alternative either a linear or a higher-order polynomial could be chosen to describe the operating characteristics of this pump. Now let us revisit the network in Fig. 4.10 that contains a BPV as a second illustration of forming the ∆Q-equations. In this analysis we can visualize the two sets of loops as shown in Fig. 4.17. The ∆Q loops ignore the presence of the BPV in this network, but the energy loops will be written for the modified network with the BPV converted into an artificial reservoir. The resulting ∆Q-equations for this network appear as Eqs. 4.31. In F1 = K1 Qo1 + ∆Q2 ( )n1 − hp1 − K2 Qo2 − ∆Q2 − ∆Q3 ( )n2 + 20 = 0 F2 = K1 Qo1 + ∆Q2 ( )n1 − hp1 + K3 Qo3 − ∆Q3 ( )n3 + K4 ' Qo4 − ∆Q1 − ∆Q3 ( )n4 +15 = 0 F3 = K9 Qo9 + ∆Q3 ( )n9 + K6 Qo6 − ∆Q1 ( )n6 − K7 Qo7 + ∆Q1 ( )n7 − K8 Qo8 + ∆Q1 ( )n8 + K4 ' Qo4 − ∆Q1 − ∆Q3 ( )n4 + 60 = 0 (4.31) P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) 200 m ∆Q1 (8) (9) [6] [5] 150 - 2000 ∆Q2 ∆Q3 0.015 m3/s 0.015 m3/s HGL = 195 m 200 - 2000 800 m BPV 200 - 1500 150mII 0.02 m3/s 300-1000140 m 250 - 1200 150-10060 m 135 m 70 m 0.03 m3/s 100 m 150 - 1500 180 m III 150 - 1200 80 m 0.02 m3/s 0.02 m3/s 250 - 1000 I Diameters in mm Lengths in m Figure 4.17 The network of Fig. 4.10, modified for the ∆Q-equation system. these equations hp1 = A Qo1 + ∆Q2 ( )2 + B Qo1 + ∆Q2 ( ) + C . The initial flows that satisfy the junction continuity equations are chosen as Qo1 = 0.1, Qo2 = 0.0, Qo3 = 0.085, Qo4 = 0.015, Qo5 = 0.0, Qo6 = 0.0, Qo7 = 0.02, Qo8 = 0.05, and Qo9 = 0.02. The substitution of these values into Eqs. 4.31 yields the final set of ∆Q-equations. If large differences in ground elevation occur in a network, PRV's are often installed in a sequence of pipes to prevent excessively large pressures in the lower part of the network. Such a series of PRV's may cause pressures in one subregion to be completely independent of the remainder of the network. Such isolation creates what are commonly called separate pressure zones. When separate pressure zones are created, it is normally better to form sub- networks and analyze each one separately, starting with the subnetwork at the lower eleva- tion. The solution from the isolated lower subnetwork can then be used to determine the demands at the nodes of the next higher network, and so on. The 10-pipe, 6-node network in Fig. 4.18 contains three PRV's in pipes 4, 5, and 7, respectively; it typifies such a situation. In this network the three PRV's cause the pressures at nodes 4 and 6 to be independent of pressures in the remainder of the network. The best analysis, therefore, would begin by studying separately the subnetwork that is composed of pipes 5, 4, 7, and 8 downstream from the PRV's. In this subnetwork the PRV's are modeled as three constant-head reservoirs. The values of Q4, Q5, and Q7 from the solution of the subnetwork are next added to the other demands to determine the demands at nodes 3, 2, and 5, respectively, in organizing the remainder of the network for analysis. P1 400' P2 500' [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] (10) [6] 3.5 ft3/s 400' 6" - 1000' 6" - 2000' 10"-1000'400' 350' 1.0 ft3/s 0.5 ft3/s 0.5 ft3/s 60' 6" - 3300' 400' 6" - 800' (PRV)2 (HGL)2 = 150' 60' 6"-2000'8"-1700'(HGL)1 = 150' (PRV)1 500'400' (PRV)3 (HGL)3 = 150' 700'8"-3000'10"-1000'e = 0.002" All pipes 6"-2500' Pump Characteristics Pump 1 Pump 2 Q ft3/s H ft. Q ft3/s H ft. 1.5 110 1.5 15.0 2.5 104 2.0 12.0 3.5 92 3.0 7.5 Figure 4.18 A network with two pressure zones. While it is generally not difficult to determine by visual examination of a map of the piping system whether PRV's isolate a portion of a network into a separate pressure zone, in computer programs a simple test is needed to identify this situation. Such a test can be based on the fact that no series of connected pipes exists between any of the artificial reservoirs created by the PRV's and any of the other reservoirs and source pumps. That no connection exists in the network example can be seen by resketching the network, as shown in Fig. 4.19. As a consequence, if pseudo loops between artificial reservoirs or source pumps cannot be found by a computer program that uses its own internal loop- finding algorithm, then the PRV's isolate a subnetwork into a pressure zone that is separate from the remainder of the network. One difficulty with this kind of test, which relies on the inability to find paths which connect all supply sources, is that errors in the network input data or an ill-defined network itself can also cause this test to be satisfied; network computer programs are supposed to identify such input errors and terminate if any such errors are found. Thus it is desirable to have an independent verification, i.e., a separate test, that can indicate that separate pressure zones exist. This alternative or verification test could take the form that is described next. The goal in this "test" is to determine whether the sum NJ + NL, determined in the usual way, is equal to, or exceeds, the number of pipes in the network. P1 400' P2 500' [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] (10) [6] 3.5 ft3/s 1.0 ft3/s 0.5 ft3/s 0.5 ft3/s 150' 150' 150' Q4 Q5 Q7 Figure 4.19 PRV's create two separate pressure zones. To highlight the problem, let us examine closely the network in Fig. 4.19 which is known to have two separate pressure zones. Overall there are six junctions, so NJ = 6. There is one real loop, and the usual rule indicates there is Nres - 1 = 5 - 1 = 4 pseudo loops, or NL = 5. Using these values, we obtain NJ + NL = 6 + 5 = 11, which is larger than the actual number of pipes, which is NP = 10. In this instance the inequality occurs because only four independent loops exist, one real loop and three pseudo loops. These numbers will be found to be correct when we view the overall network as two separate networks. The higher network in Fig. 4.19 has NP = 6, one real loop and one pseudo loop, Ns = Nres - 1 and NL = 2. Since there are four junctions in the network with the higher pressure zone, NJ = 4, and NP = NJ + NL = 6. For the network with the lower pressure zone NP = 4, NJ = 2, there are no real loops, and the expected number of pseudo loops is Nres - 1 = 2, giving NL = 2. Again NP = NJ + NL. The verification test to determine whether PRV's isolate a portion of a network into a separate pressure zone might therefore be as follows: 1. Find the real loops which exist after pipes containing PRV's have been disconnected from their upstream junctions. 2. Compute NLs from NLs = Nres + Npump - 1. 3. Add the number of loops that were found in steps 1 and 2 to determine NL, and then determine NP from NP = NJ + NL. 4. If this computed NP exceeds the number of pipes in the network, then the PRV's isolate a portion of the network. The amount of the difference between the newly computed NP and the actual number of pipes in the network is the number of additional pressure zones that exist in the network; the total number of zones is one more than this number of additional zones. 4.4 SOLVING THE NETWORK EQUATIONS 4.4.1. NEWTON METHOD FOR LARGE SYSTEMS OF EQUATIONS In Sections 4.2 and 4.3 we explored the writing of systems of algebraic equations to de- scribe the relations between the discharges, pressures, and other variables and parameters in a pipe network. Many of the equations in these systems of equations are nonlinear. A good method for solving nonlinear equations is therefore needed. Numerous methods exist, but the Newton Method is the method of choice here. Its application to the solution of the Q-equations, the H-equations and the ∆Q-equations will be discussed in this section. To treat the unknown discharges (when using the Q-equations), the unknown heads (when using the H-equations), and the unknown corrective loop discharges (when using the ∆Q- equations) in a uniform way, the primary unknown variable in this section will be called the vector {x}. The Newton iterative formula for solving a system of equations can be written as {x }(m+1) = {x }(m) − [ D ]−1{F}(m) (4.32a) Here x is an entire column vector {x} of unknowns, {F} is an entire column vector of equations, and [D]-1 is the inverse of a matrix [D] which is the Jacobian. The Jacobian occurs in several applications in mathematics, and it represents the following matrix of derivatives: D [ ] = ∂F1 ∂x1 ∂F1 ∂x2 ⋅ ⋅ ∂F1 ∂xn ∂F2 ∂x1 ∂F2 ∂x2 ⋅ ⋅ ∂F2 ∂xn ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ∂Fn ∂x1 ∂Fn ∂x2 ⋅ ⋅ ∂Fn ∂xn (4.33) Likewise {x} and {F} are actually x { } = x1 x2 ⋅ ⋅ xn F { } = F1 F2 ⋅ ⋅ Fn (4.34) Equation 4.32a indicates that the Newton method solves a system of nonlinear equations by iteratively solving a system of linear equations because [D]-1{F} represents the solution of the linear system of equations [ D ]{ z } = {F} (4.32b) That is, the vector that is subtracted from the current estimate of the unknown vector {x} in Eq. 4.32a is the solution {z} to the linear system of equations that is Eq. 4.32b. In practice we therefore see that the Newton method solves a system of equations by the iterative formula {x }(m+1) = {x }(m) − {z } (4.32c) where {z} is the solution vector that is obtained by solving [D]{z} = {F}. If the system should actually contain only linear equations, then the first iteration will produce the exact solution. The development of Eq. 4.32 follows. We begin by using a multi-dimensional Taylor series expansion to evaluate the individual equations Fi in the neighborhood of an initial solution estimate that we call {x} which is presumed to be near the actual solution: F1 (m+1) = F1 (m) + ∂F1 ∂x1 ∆x1 + ∂F1 ∂x2 ∆x2 +⋅ ⋅ + ∂F1 ∂xn ∆xn + Ο ∆x 2 ( ) = 0 F2 (m+1) = F2 (m) + ∂F2 ∂x1 ∆x1 + ∂F2 ∂x2 ∆x2 +⋅ ⋅ + ∂F2 ∂xn ∆xn + Ο ∆x 2 ( ) = 0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ Fn (m+1) = Fn (m) + ∂Fn ∂x1 ∆x1 + ∂Fn ∂x2 ∆x2 +⋅ ⋅ + ∂Fn ∂xn ∆xn + Ο ∆x 2 ( ) = 0 (4.35) When we use matrix notation and make the substitution ∆xi = xi (m+1) − xi (m) , this sys- tem of equations becomes F1 F2 ⋅ ⋅ Fn (m) + ∂F1 ∂x1 ∂F1 ∂x2 ⋅ ⋅ ∂F1 ∂xn ∂F2 ∂x1 ∂F2 ∂x2 ⋅ ⋅ ∂F2 ∂xn ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ∂Fn ∂x1 ∂Fn ∂x2 ⋅ ⋅ ∂Fn ∂xn x1 (m+1) − x1 (m) x2 (m+1) − x2 (m) ⋅ ⋅ xn (m+1) − xn (m) = 0 (4.36) which can be written compactly as {F}(m) + [D](m)({x}(m+1) - {x}(m)) = {0} and solved for {x}(m+1) to produce Eq. 4.32a. Now let us see in some detail how the Newton method works in practice. First we shall examine the three-reservoir problem by forming and solving manually the appropriate systems of Q-equations, H-equations and ∆Q-equations. Then we will look at computer programs that could be used to find the solution to the Q-equations; the first program is simpler and more specialized, and the second program is longer but more versatile. Finally we examine a third program that will solve any equation system that is supplied to it in a subroutine. 100 m 60 m 85 m II I 0.3 - 2000 [1] (1) (2) (3) e = 0.0005 m v = 1.31 x 10-6 m2/s 0.25 - 3000 0.25 - 15 00 QJ1 = 0.06 m 3/s Diameters in m Lengths in m Pipe K n 1 1469 1.974 2 2432 1.927 3 5646 1.971 Figure 4.20 The three-reservoir problem. Figure 4.20 shows three reservoirs that are connected by three pipes with an external de- mand at the common junction of the pipes. The highest reservoir has a water surface elevation of 100 m; the middle reservoir water surface elevation is 85 m, and the lowest reservoir has a water surface elevation of 60 m. We will use the data in the figure and table to form and solve the three systems of equations. The Q-equations are F1 = Q1 + Q2 − Q3 − QJ1 = 0 F1 = Q1 + Q2 − Q3 − 0.06 = 0 F2 = K1Q1 n1 − K2Q2 n2 − WS1 + WS2 = 0 F2 = 1469Q1 1.974 − 2432Q2 1.927 −15 = 0 F3 = K1Q1 n1 + K3Q3 n3 − WS1 + WS3 = 0 F3 = 1469Q1 1.974 + 5646Q3 1.971 − 40 = 0 (4.37) To satisfy the junction continuity equation, equation F1, and also determine initial values for the Newton method, we can select Q1 (0) = Qo1 = 0.10 m3/s, Q2 (0) = Qo2 = 0.05 m3/s, and Q3 (0) = Qo3 = 0.09 m3/s. The superscript (0) denotes initial values for use by the Newton method in solving the Q-equations, and the subscripts denote initial discharge values for the ∆Q-equations. The initial values for use with the Q-equations are not required to satisfy the junction continuity equations, although this set of values does. The H-equation is F1 = WS1 − H1 K1 1/n1 + WS2 − H1 K2 1/n2 − H1 − WS3 K3 1/n3 − 0.06 = 0 F1 = 100 − H1 1469 0.507 + 85 − H1 2432 0.519 − H1 − 60 5646 0.507 − 0.06 = 0 (4.38) The ∆Q-equations are F1 = K1 Qo1 + ∆Q1 + ∆Q2 ( )n1 − K2 Qo2 − ∆Q1 ( )n2 − WS1 + WS2 = 0 F2 = K1 Qo1 + ∆Q1 + ∆Q2 ( )n1 + K3 Qo3 + ∆Q2 ( )n3 − WS1 + WS3 = 0 (4.39a) With the initial discharges that we have chosen, these equations become F1 = 1469 0.10 + ∆Q1 + ∆Q2 ( )1.974 − 2432 0.05 − ∆Q1 ( )1.927 −15 = 0 F2 = 1469 0.10 + ∆Q1 + ∆Q2 ( )1.974 + 5646 0.09 + ∆Q2 ( )1.971 − 40 = 0 (4.39b) We now begin the solution of the Q-equations by the Newton method using the equa- tions [D]{z} = {F}, {Q}(m+1) = {Q}(m) - {z}. According to Eq. 4.33, the Jacobian is D [ ] = 1 1 −1 n1K1Q1 n1−1 −n2K2Q2 n2 −1 0 n1K1Q1 n1−1 0 n3K3Q3 n3−1 (4.40a) or D [ ] = 1 1 −1 2900Q1 0.974 −4686Q2 0.927 0 2900Q1 0.974 0 11128Q3 0.971 (4.40b) For the first computational cycle we therefore solve the equation set 1.00 1.00 −1.00 307.89 − 291.57 0.00 307.89 0.00 1073.96 z1 z2 z3 = 0.00 − 6.97 24.64 (4.41a) and obtain the results z { } = − 0.0004 0.0235 0.0231 Q { } = 0.1004 0.0265 0.0669 (4.42a) We now iterate to obtain the following equation set and updated solution: 1.00 1.00 −1.00 309.09 − 161.86 0.00 309.09 0.00 805.20 z1 z2 z3 = 0.000 − 1.506 3.051 (4.41b) z { } = − 0.0017 0.0061 0.0044 Q { } = 0.1021 0.0204 0.0625 (4.42b) One more iteration leads to 1.00 1.00 −1.00 314.19 − 127.01 0.00 314.19 0.00 753.42 z1 z2 z3 = 0.000 − 0.095 0.151 (4.41c) z { } = − 0.0001 0.0004 0.0003 Q { } = 0.1022 0.0200 0.0622 m3/s (4.42c) This solution is now sufficiently accurate! The solution of the H-equations by the Newton method uses basically the same equa- tions [D]{z} = {F} and {H}(m+1) = {H}(m) - {z}. These relations lead to a single update equation H1 (m+1) = H1 (m) − F1 / ( dF1 dH1 ) (4.43) The derivative is dF1 dH1 = − 1 n1K1 100 − H1 K1 1/n1−1 − 1 n2K2 85 − H1 K2 1/n2 −1 − 1 n3K3 H1 − 60 K3 1/n3−1 (4.44a) or dF1 dH1 = − 1 2900 100 − H1 1469 −0.493 − 1 4686 85 − H1 2432 −0.481 − 1 11128 H1 − 60 5646 −0.493 (4.44b) If we initiate the solution procedure with the initial estimate of H1 (0) = 84 , then the first two iterative cycles produce F1 = − 0.00415, dF1 dH1 = − 0.0136, H1 (1) = 84 − 0.00415 / 0.0136 = 83.70 (4.45a) and F1 = − 0.000247, dF1 dH1 = − 0.01251, H1 (2) = 83.70 − 0.000247 / 0.01251 = 83.68m (4.45b) which will be regarded as adequate. Finally, we now solve the ∆Q-equations by the Newton method using again the equa- tions [D]{z} = {F} and {∆Q(m+1)} = {∆Q(m)} - {z}. In this case we solve repeatedly the two-equation system for updated correction vectors {z} until it is declared to be sufficiently small. Three cycles of computation yield these results: 599 307 307 1382 z1 z2 = − 6.97 24.6 z1 z2 = − 0.0234 0.0230 ∆Q1 ∆Q2 = 0.0234 − 0.0230 (4.46a) 472 309 309 1115 z1 z2 = −1.523 3.13 z1 z2 = − 0.0062 0.0045 ∆Q1 ∆Q2 = 0.0296 − 0.0275 (4.46b) 441 314 314 1068 z1 z2 = −0.0952 0.1507 z1 z2 = − 0.0004 0.0003 ∆Q1 ∆Q2 = 0.0300 − 0.0278 (4.46c) Now the discharges can be computed as Q1 = 0.1 + ∆Q1 + ∆Q2 = 0.1022 m3/s, Q2 = 0.05 - ∆Q1 = 0.0200 m3/s, and Q3 = 0.09 + ∆Q2 = 0.0622 m3/s. Computer programs of differing complexity and generality can also be developed for the solution of these equation systems by application of the Newton method. We will now look at two programs. The first program is relatively simple but must be recoded in part for each application; it will be applied to the solution of the Q-equations for the three- reservoir problem. The second program is more versatile. Program 4.2, the FORTRAN program listed in Fig. 4.21, is designed to solve three si- multaneous equations with the Newton method. It calls a matrix solver that has the coef- ficient matrix expanded by one column to contain the known vector, and it places the inverse in additional columns beyond the location of the known vector. The first part of the main program is currently written specifically to solve the Q-equations for the three- reservoir problem. However, the portion that numerically evaluates the derivatives in the Jacobian matrix is written more generally, with N giving the size of the matrix problem to be solved. Careful study of this listing will clarify considerably how the various tasks are performed. The subroutine INVERM employs a common method in linear algebra problems by using an expanded matrix. The coefficient matrix is square, here 3 rows by 3 columns. The known vector is placed in the next column, in this case column 4. The subroutine solves the system of equations and provides the inverse matrix. The solution is returned in the same column that initially contained the known vector, here column 4. ************************************************************************* * PROGRAM NO. 4.2, NEWTON, FORTRAN * THIS PROGRAM HAS BEEN INCLUDED FOR THE CONVENIENCE OF THE READER. * THE AUTHOR ACCEPTS NO RESPONSIBILITY FOR ITS CORRECTNESS. * USERS OF THIS PROGRAM DO SO AT THEIR OWN RISK. ************************************************************************* * IMPLEMENTS THE NEWTON METHOD IN SOLVING THREE EQUATIONS REAL X(3),F(3),D(3,7),RK(3),RN(3),K1,K2,K3,N1,N2,N3 DATA N,N1/3,4/,DX/.001/,MAX/15/,ERR/.0001/ WRITE(5,*)' GIVE: K1,K2,K3,N1,N2,N3,Q1,Q2,Q3' READ(5,*) RK,RN,X M=0 1 NT=0 5 F(1)=X(1)+X(2)-X(3)-0.06 F(2)=RK(1)*X(1)**RN(1)-RK(2)*X(2)**RN(2)-15.0 F(3)=RK(1)*X(1)**RN(1)+RK(3)*X(3)**RN(3)-40.0 IF(NT.NE.0) GO TO 15 DO 10 I=1,N 10 D(I,N1)=F(I) X(1)=X(1)-DX NT=1 GO TO 5 15 X(NT)=X(NT)+DX DO 20 I=1,N 20 D(I,NT)=(D(I,N1)-F(I))/DX NT=NT+1 IF(NT.GT.N) GO TO 30 X(NT)=X(NT)-DX GO TO 5 30 CALL INVERM(D,N) DIF=0. DO 40 I=1,N DIF=DIF+ABS(D(I,N1)) 40 X(I)=X(I)-D(I,N1) M=M+1 IF(DIF.GT. ERR .AND. M.LT.MAX) GO TO 1 WRITE(5,*)' THE SOLUTION IS ',X END Figure 4.21 Program 4.2 to use the Newton method to solve three equations. The listing of the inverse starts in the next column. In solving three equations the array D(3,7) therefore has 7 as its second subscript, and the last three columns contain the inverse. Currently the input data to this program includes the coefficients K and n for each pipe and an initial estimate of the discharge in each pipe: 1469 2432 5646 1.974 1.927 1.971 0.10 0.05 0.09. The program produces the same solution as we obtained manually. The next program, listed in Fig. 4.22, is essentially the same as the previous program, except that it calls the linear equation solver GAUSEL, described in Appendix A, rather than INVERM. GAUSEL is a more versatile subroutine that interchanges rows to minimize truncation error, applies one iterative correction to the solution vector, and returns an estimate of the relative error for each unknown in the array ERRNOR, so the user has parameters to determine the accuracy of the solution. However, the relative error is not printed in this program. This subroutine also illustrates the Microsoft FORTRAN ability to allocate the array sizes that are needed. ************************************************************************* * PROGRAM NO. 4.3, NEWTON, FORTRAN * THIS PROGRAM HAS BEEN INCLUDED FOR THE CONVENIENCE OF THE READER. * THE AUTHOR ACCEPTS NO RESPONSIBILITY FOR ITS CORRECTNESS. * USERS OF THIS PROGRAM DO SO AT THEIR OWN RISK. ************************************************************************* * IMPLEMENTS THE NEWTON METHOD IN SOLVING THREE EQUATIONS REAL X(3),F(3),D(3,3),F1(3),ERRNOR(3),RK(3),RN(3),K1,K2,K3,N1,N2,N3 DATA N/3/,DX/.001/,MAX/15/,ERR/.0001/ WRITE(5,*)' GIVE: K1,K2,K3,N1,N2,N3,Q1,Q2,Q3' READ(5,*) RK,RN,X M=0 1 NT=0 5 F(1)=X(1)+X(2)-X(3)-0.06 F(2)=RK(1)*X(1)**RN(1)-RK(2)*X(2)**RN(2)-15.0 F(3)=RK(1)*X(1)**RN(1)+RK(3)*X(3)**RN(3)-40.0 IF(NT.NE.0) GO TO 15 DO 10 I=1,N 10 F1(I)=F(I) X(1)=X(1)-DX NT=1 GO TO 5 15 X(NT)=X(NT)+DX DO 20 I=1,N 20 D(I,NT)=(F1(I)-F(I))/DX NT=NT+1 IF(NT.GT.N) GO TO 30 X(NT)=X(NT)-DX GO TO 5 30 CALL GAUSEL(3,3,D,F1,DET,ERRNOR) DIF=0. DO 40 I=1,N DIF=DIF+ABS(F1(I)) 40 X(I)=X(I)-F1(I) M=M+1 IF(DIF.GT. ERR .AND. M.LT.MAX) GO TO 1 WRITE(5,*)' THE SOLUTION IS ',X END Figure 4.22 Program 4.3, a more versatile implementation of the Newton method. With this introduction to the Newton method, let us look further at the structure of a computer program that would solve any system of equations. This program should consist of two primary elements: First, a main or driver program that accomplishes the following tasks: (a) it allows the user to assign values to known variables and initial estimates for unknown variables; (b) it creates the Jacobian matrix and the known vector and supplies values to these arrays; (c) it calls a linear algebra solver; (d) it implements the Newton iteration; and (e) it prints the solution. Second, it must contain a subroutine (or function subprogram) that defines the system of equations to be solved. The equations in this subroutine will change, depending upon the nature of the problem that is being solved, and therefore the statements in this subroutine would be changed as different types of problems are to be solved. A listing of such a general purpose program EQUSOL1.FOR will be found in Fig. 4.23. The subroutine FUNCT provides the equations for this program. The main program calls on the linear algebra solver SOLVEQ that is described in Appendix B. ************************************************************************* * PROGRAM NO. 4.4, EQUSOL1, FORTRAN * THIS PROGRAM HAS BEEN INCLUDED FOR THE CONVENIENCE OF THE READER. * THE AUTHOR ACCEPTS NO RESPONSIBILITY FOR ITS CORRECTNESS. * USERS OF THIS PROGRAM DO SO AT THEIR OWN RISK. ************************************************************************* * THIS EQUATION SOLVER IMPLEMENTS THE NEWTON METHOD LOGICAL IV[ALLOCATABLE](:) INTEGER*2 INDX[ALLOCATABLE](:) CHARACTER*3 SYMB[ALLOCATABLE](:),CH REAL X[ALLOCATABLE](:),F[ALLOCATABLE](:),F1[ALLOCATABLE](:), &D[ALLOCATABLE](:,:) * WRITE(*,*)' GIVE (1) NO. OF EQS., (2) NO. OF VARIABLES,', &' (3) INPUT UNIT, (4) OUTPUT UNIT' READ(*,*) N,NV,IN,IOUT ALLOCATE(X(NV),F(N),F1(N),D(N,N),SYMB(NV),IV(NV),INDX(N)) IF(IN.EQ.0 .OR. IN.EQ.5) WRITE(IN,100) NV 100 FORMAT(' GIVE',I3,' LINES WITH:',/,3X,'(1) SYMBOL FOR VAR.(3 CH)', &/,3X,'(2) K OR U FOR KNOWN OR UNKNOWN AND',/,3X,'(3) VALUE') J=0 DO 10 I=1,NV READ(IN,110) SYMB(I),CH,X(I) 110 FORMAT(A3,1X,A1,1X,F10.0) IF(CH.EQ.'U' .OR. CH.EQ.'u') THEN IV(I)=.TRUE. J=J+1 ELSE IF(CH.EQ.'K' .OR. CH.EQ.'k') THEN IV(I)=.FALSE. ELSE WRITE(*,*)' ERROR IN INPUT FOR VARIABLE', I STOP ENDIF 10 CONTINUE IF(J.EQ.N) GO TO 12 WRITE(*,*)' YOU GAVE',N,' EQS. BUT',J,' UNKNOWNS' STOP 12 NCT=0 15 CALL FUNCT(X,F) J=0 DO 30 JJ=1,NV IF(IV(JJ)) THEN XX=X(JJ) X(JJ)=1.005*X(JJ) J=J+1 CALL FUNCT(X,F1) DO 20 I=1,N 20 D(I,J)=(F1(I)-F(I))/(X(JJ)-XX) X(JJ)=XX ENDIF 30 CONTINUE CALL SOLVEQ(N,1,N,D,F,1,DD,INDX) NCT=NCT+1 SUM=0. Figure 4.23 A listing of the program EQUSOL1.FOR. J=0 DO 40 I=1,NV IF(IV(I)) THEN J=J+1 X(I)=X(I)-F(J) SUM=SUM+ABS(F(J)) ENDIF 40 CONTINUE WRITE(*,*)' NCT =',NCT,' SUM =',SUM IF(NCT.LT.20 .AND. SUM.GT. 0.0001) GO TO 15 WRITE(IOUT,120)(I,SYMB(I),X(I),I=1,NV) 120 FORMAT(I5,1X,A3,' =',F10.3) END * SUBROUTINE FUNCT(X,F) REAL F(22),X(41) DATA E,G2,P,AP/0.005,64.4,8.6714174E-6,5.4541539E-3/ F(1)=X(41)-X(1)-X(3)-X(35) !Unknowns F(2)=X(1)-X(2)-X(36) ! 1=Q2 15=H4 29=L1 F(3)=X(3)-X(4)-X(37) ! 2=Q3 16=H5 30=L2 F(4)=X(2)+X(4)-X(5)-X(38) ! 3=Q4 17=f1 31=L3 F(5)=X(12)-X(40)+X(6) ! 4=Q5 18=f2 32=L4 F(6)=X(13)-X(12)+X(7) ! 5=Q6 19=f3 33=L5 F(7)=X(14)-X(12)+X(9) ! 6=hf1 20=f4 34=L6 F(8)=X(15)-X(13)+X(8) ! 7=hf2 21=f5 35=QJ1 F(9)=X(16)-X(15)+X(11) ! 8=hf3 22=f6 36=QJ2 F(10)=X(7)+X(8)-X(10)-X(9) ! 9=hf4 Knowns 37=QJ3 RF1=1./SQRT(X(17)) ! 10=hf5 23=D1 38=QJ4 RF2=1./SQRT(X(18)) ! 11=hf6 24=D2 39=QJ5 RF3=1./SQRT(X(19)) ! 12=H1 25=D3 40=WS1 RF4=1./SQRT(X(20)) ! 13=H2 26=D4 41=Q1 RF5=1./SQRT(X(21)) ! 14=H3 27=D5 RF6=1./SQRT(X(22)) ! 15=H4 28=D6 F(11)=X(6)-X(17)*X(29)*12./X(23)*(X(41)/(AP*X(23)**2))**2/G2 F(12)=X(7)-X(18)*X(30)*12./X(24)*(X(1)/(AP*X(24)**2))**2/G2 F(13)=X(8)-X(19)*X(31)*12./X(25)*(X(2)/(AP*X(25)**2))**2/G2 F(14)=X(9)-X(20)*X(32)*12./X(26)*(X(3)/(AP*X(26)**2))**2/G2 F(15)=X(10)-X(21)*X(33)*12./X(27)*(X(4)/(AP*X(27)**2))**2/G2 F(16)=X(11)-X(22)*X(34)*12./X(28)*(X(5)/(AP*X(28)**2))**2/G2 F(17)=RF1-1.14+2.*ALOG10(E/X(23)+P*X(23)*RF1/X(41)) F(18)=RF2-1.14+2.*ALOG10(E/X(24)+P*X(24)*RF2/X(1)) F(19)=RF3-1.14+2.*ALOG10(E/X(25)+P*X(25)*RF3/X(2)) F(20)=RF4-1.14+2.*ALOG10(E/X(26)+P*X(26)*RF4/X(3)) F(21)=RF5-1.14+2.*ALOG10(E/X(27)+P*X(27)*RF5/X(4)) F(22)=RF6-1.14+2.*ALOG10(E/X(28)+P*X(28)*RF6/X(5)) RETURN END Figure 4.23 (Concluded) A listing of the program EQUSOL1.FOR. Let's examine how the main program does its tasks of providing values to the Jacobian matrix [D] and equation vector {F} and then carrying out a Newton solution. The key portion of the main program that implements the Newton method appears in bold characters in Fig. 4.23. Three tasks are accomplished by these statements: (1) defining the equation vector; (2) numerically evaluating the elements of the Jacobian matrix; and (3) solving the resulting linear system of equations and subtracting this solution from the current vector of unknowns, as described by Eq. 4.32a. The FORTRAN integer NCT is the iteration counter; it is set to 0 before beginning the Newton iteration. Statement 15 CALL FUNCT(X,F) has two arguments, an array X for the variables and an array F for the equations. The array X includes both the known and unknown variables of the problem. Upon returning from CALL FUNCT, the array F contains a set of equations that have been evaluated by using the initial estimates of the unknowns. Since the initial estimates are incorrect, the individual elements of {F} will not be zero, but subsequent Newton iterations will drive these elements progressively closer to zero. Statement DO 30 JJ=1,NV, in which NV is the total of all variables, evaluates individual columns of the Jacobian matrix D(I,J) by using a first-order numerical evaluation of the derivatives. Since IV(JJ) is .FALSE. for known variables and .TRUE. for unknown variables, we note that nothing happens in loop 30 if IV is .FALSE. Hence J, which identifies the column in which the Jacobian derivatives are entered, is incremented only for unknown variables. When an unknown is encountered, xj, which is X(JJ), is incremented by multiplying its current value by 1.005 before the equation is evaluated again by calling FUNCT. Upon returning from FUNCT, the array F1 now contains equation values based on 1.005xj, and then the statement D(I,J) = (F1(I) - F(I))/(X(JJ) - XX) numerically evaluates the derivatives of the equations by using a first-order approximation. The statement DO 20 I=1,N fills all row entries for column J of the Jacobian matrix [D]. Upon completing the DO 30 loop, the equation vector {F} and the Jacobian matrix [D] have been fully evaluated. The next statement CALL SOLVEQ calls a linear equation solver, which upon return has replaced the elements of the array F with the solution vector {z} found in Eq. 4.32b. The statement DO 40 I=1,NV implements Eq. 4.32c with SUM accumulating the absolute sum of the corrections applied to the unknown vector {x}. If this SUM is larger than the allowable error and fewer than 20 iterations have been completed, then GO TO 15 at the end of this code segment will begin another Newton iteration. In our example it would be relatively easy to derive the actual partial derivatives of each equation with respect to each unknown, and the elements of the Jacobian could be evaluated by using these derivatives. The length of the program would be longer if these derivative expressions were used. The numerical approximation of the derivatives requires extra arith- metic, particularly since many derivatives are zero, but the advantage of a shorter code makes the numerical approximation of the derivatives attractive. Example Problem 4.5 Use program EQUSOL1 to solve the 6-pipe, 5-node network shown below. Obtain this solution in four ways: (1) use the program as it now exists with subroutine FUNCT; (2) use the Q-equations; (3) use the H-equations; (4) use the ∆Q-equations. [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) [5] All e = 0.005" QJ5 = 0.25 ft 3/s QJ3 = 0.5 ft 3/s QJ2 = 0.35 ft 3/s QJ1 = 0.5 ft 3/s QJ4 = 0.5 ft 3/s All elev. = 350' 8" - 1500' 4" - 1000' WS1 = 500' 6"-1000'6" - 1500' 6" - 1 20 0' 6" - 1500' 1. The existing subroutine FUNCT explicitly defines the equations that we want to solve; there are 22 equations, F(1) through F(22). There are 41 variables associated with the solution; therefore 41 - 22 = 19 of these variables are known. These equations are as follows: Junction continuity equations: Q1 − Q2 − Q4 − QJ1 = 0 (1) Q2 − Q3 − QJ2 = 0 (2) Q4 − Q5 − QJ3 = 0 (3) Q3 + Q5 − Q6 − QJ4 = 0 (4) (The junction continuity equation at node 5 is not included here, but this simple equation Q6 - QJ5 = 0 establishes the discharge in pipe 6 as 0.25 ft3/s.) Head loss equations giving the HGL at a downstream node relative to the upstream node: H1 = WS1 − h f 1 (5) H2 = H1 − h f 2 (6) H3 = H1 − h f 4 (7) H4 = H2 − h f 3 (or H4 = H3 − h f 5) (8) H5 = H4 − h f 6 (9) Energy equation around a loop: h f 2 + h f 3 − h f 5 − h f 4 = 0 (10) Darcy-Weisbach equations to define the frictional head losses (pipe numbers i = 1, 6): h fi = f i Li Di Qi 2 2gAi 2 (11-16) Colebrook-White equations (pipe numbers i = 1, 6): 1 f i = 1.14 − 2 log10 ei Di + 9.35νDi 4/π ( )Qi f i (17-22) In the program listing the integer within X( ) identifies the variable in the array, as is seen by the comments following the exclamation points there. We note there are 22 equations: a continuity equation for NJ - 1 = 4 junctions, a separate head difference equation for each pipe, a Darcy-Weisbach equation for each pipe, a companion Colebrook- White equation for each pipe, and finally an energy loop equation, for a total of 3NP + NJ = 18 + 5 = 22 equations. Since the entire system demand must come from pipe 1, its discharge must be Q1 = 2.1 ft3/s, and the unknowns are 5 unknown discharges Q2 .. Q6, 6 unknown head losses hf1 .. hf6, 5 unknown heads H1 .. H5, and 6 unknown friction factors f1 .. f6, for a total of 22. The input and solution files for the program now follow: Input File Output File Q2 U 0.8 f6 U 0.020 1 Q2 = 0.820 22 f6 = 0.024 Q3 U 0.5 D1 K 8.0 2 Q3 = 0.470 23 D1 = 8.000 Q4 U 0.8 D2 K 6.0 3 Q4 = 0.780 24 D2 = 6.000 Q5 U 0.3 D3 K 6.0 4 Q5 = 0.280 25 D3 = 6.000 Q6 U 0.3 D4 K 6.0 5 Q6 = 0.250 26 D4 = 6.000 hf1 U 24.0 D5 K 6.0 6 hf1 = 23.952 27 D5 = 6.000 hf2 U 11.0 D6 K 4.0 7 hf2 = 11.279 28 D6 = 4.000 hf3 U 0.2 L1 K 1500.0 8 hf3 = 5.872 29 L1 = 1500.0 hf4 U 0.15 L2 K 1000.0 9 hf4 = 15.368 30 L2 = 1000.0 hf5 U 2.0 L3 K 1500.0 10 hf5 = 1.783 31 L3 = 1500.0 hf6 U 6.0 L4 K 1500.0 11 hf6 = 9.117 32 L4 = 1500.0 H1 U 476.0 L5 K 1200.0 12 H1 = 476.048 33 L5 = 1200.0 H2 U 465.0 L6 K 1000.0 13 H2 = 464.769 34 L6 = 1000.0 H3 U 460.0 QJ1 K 0.50 14 H3 = 460.680 35 QJ1 = 0.50 H4 U 458.0 QJ2 K 0.35 15 H4 = 458.897 36 QJ2 = 0.35 H5 U 450.0 QJ3 K 0.50 16 H5 = 449.780 37 QJ3 = 0.50 f1 U 0.020 QJ4 K 0.50 17 f1 = 0.019 38 QJ4 = 0.50 f2 U 0.020 QJ5 K 0.25 18 f2 = 0.021 39 QJ5 = 0.25 f3 U 0.020 WS1 K 500.0 19 f3 = 0.022 40 WS1 = 500.0 f4 U 0.020 Q1 K 2.1 20 f4 = 0.021 41 Q1 = 2.100 f5 U 0.020 21 f5 = 0.024 2. The Q-equations are F1 = Q1 − Q2 − Q4 − QJ1 = 0 F2 = Q2 − Q3 − QJ2 = 0 F3 = Q4 − Q5 − QJ3 = 0 F4 = Q3 + Q5 − Q6 − QJ4 = 0 F5 = K2Q2 n2 + K3Q3 n3 − K5Q5 n5 − K4Q4 n4 = 0 The K and n for each pipe must now be determined. Program 2.1, PIPK_N, or some other means will provide these values: Pipe K n 1 5.6845 1.9381 2 16.4967 1.9185 3 24.3685 1.8858 4 24.7450 1.9185 5 19.0411 1.8611 6 126.3843 1.8970 To compute the five unknown discharges Qi (i = 2, 6) (with Q1 = 2.1 ft3/s known), the subroutine FUNCT must be modified as follows: SUBROUTINE FUNCT(X,F) REAL F(5),X(11) REAL K2/16.4967/,K3/24.3685/,K4/24.745/,K5/19.0411/ REAL N2/1.9185/,N3/1.8858/,N4/1.9185/,N5/1.8611/ F(1)=X(6)-X(1)-X(3)-X(7) ! Unknowns Knowns F(2)=X(1)-X(2)-X(8) ! 1 = Q2, 4 = Q5, 6 = Q1, 9 = QJ3 F(3)=X(3)-X(4)-X(9) ! 2 = Q3, 5 = Q6, 7 = QJ1, 10 = QJ4 F(4)=X(2)+X(4)-X(5)-X(10) ! 3 = Q4, 8 = QJ2, 11 = QJ5 F(5)=K2*X(1)**N2+K3*X(2)**N3-K5*X(4)**N5-K4*X(3)**N4 RETURN END The input data (all in ft3/s) that were used to solve this problem (with 5 and 11 plus 2 and 3 for I/O units from the keyboard) and the solution are listed now: Input Data Solution Variable Type Initial value Index Value Q2 U 0.80 1 Q2 0.82 Q3 U 0.50 2 Q3 0.47 Q4 U 0.80 3 Q4 0.78 Q5 U 0.30 4 Q5 0.28 Q6 U 0.25 5 Q6 0.25 Q1 K 2.10 6 Q1 2.10 QJ1 K 0.50 7 QJ1 0.50 QJ2 K 0.35 8 QJ2 0.35 QJ3 K 0.50 9 QJ3 0.50 QJ4 K 0.50 10 QJ4 0.50 QJ5 K 0.25 11 QJ5 0.25 Since X(11) = QJ5 is not used in the equations, the keyboard input could have been changed to 5 10 2 3, and the last line of input could then be deleted. 3. The number of H-equations could be reduced below five, but we will use five head equations to determine the head at the five nodes. These equations are F1 = 500 − H1 K1 1/n1 − H1 − H2 K2 1/n2 − H1 − H3 K4 1/n4 − QJ1 = 0 F2 = H1 − H2 K2 1/n2 − H2 − H4 K3 1/n3 − QJ2 = 0 F3 = H1 − H3 K4 1/n4 − H3 − H4 K5 1/n5 − QJ3 = 0 F4 = H2 − H4 K3 1/n3 + H3 − H4 K5 1/n5 − H4 − H5 K6 1/n6 − QJ4 = 0 F5 = H4 − H5 K6 1/n6 − QJ5 = 0 These equations are arranged to allow us to find Hi (i = 1, 5) with a demand at each of the five nodes as additional variables, so there are 5 unknowns and 10 variables. The appropriate modifications of subroutine FUNCT are as follows: SUBROUTINE FUNCT(X,F) REAL F(5),X(10) REAL K1/5.6845/,K2/16.4967/,K3/24.3685/,K4/24.745/,K5/19.0411/ &,K6/126.3843/,R1/.515969/,R2/.52124/,R3/.53028/,R4/.52124/ &,R5/.53732/,R6/.52715/ C UNKNOWNS: 1 = H1, 2 = H2, 3 = H3, 4 = H4, 5 = H5; C KNOWNS: 6 = QJ1, 7 = QJ2, 8 = QJ3, 9 = QJ4, 10 = QJ5 F(1)=(ABS(500-X(1))/K1)**R1-(ABS(X(1)-X(2))/K2)**R2 &-(ABS(X(1)-X(3))/K4)**R4-X(6) F(2)=(ABS(X(1)-X(2))/K2)**R2-(ABS(X(2)-X(4))/K3)**R3-X(7) F(3)=(ABS(X(1)-X(3))/K4)**R4-(ABS(X(3)-X(4))/K5)**R5-X(8) F(4)=(ABS(X(2)-X(4))/K3)**R3+(ABS(X(3)-X(4))/K5)**R5 &-(ABS(X(4)-X(5))/K6)**R6-X(9) F(5)=(ABS(X(4)-X(5))/K6)**R6-X(10) RETURN END The input data (all in ft or ft3/s) for this problem (with an additional 5 and 10 plus 2 and 3 for I/O units from the keyboard) and the solution follow: Input Data Solution Variable Type Initial value Index Value H1 U 476.0 1 H1 476.1 H2 U 465.0 2 H2 464.8 H3 U 460.0 3 H3 460.7 H4 U 458.0 4 H4 458.9 H5 U 450.0 5 H5 449.8 QJ1 K 0.50 6 QJ1 0.50 QJ2 K 0.35 7 QJ2 0.35 QJ3 K 0.50 8 QJ3 0.50 QJ4 K 0.50 9 QJ4 0.50 QJ5 K 0.25 10 QJ5 0.25 4. For this problem there is only one ∆Q-equation, which is F1 = K2 Qo2 + ∆Q1 ( )n1 + K3 Qo3 + ∆Q1 ( )n3 − K5 Qo5 − ∆Q1 ( )n5 − K4 Qo4 − ∆Q1 ( )n4 = 0 The input, since it is only two lines, can be given directly from the keyboard as 1 1 5 6 DQ1 U 0.1 The estimate of 0.1 ft3/s is used because the main program uses 1.005 times the current value. This might be changed with an IF statement that adds 0.001 to the variable if its value is zero. The solution is DQ1 = 0.020 ft3/s. The subroutine FUNCT can be modified, with the initial discharge in each pipe chosen to be Qo2 = 0.8 ft3/s, Qo3 = 0.45 ft3/s, Qo4 = 0.80 ft3/s, and Qo5 = 0.30 ft3/s to satisfy continuity, as shown: SUBROUTINE FUNCT(X,F) REAL F(1),X(1) REAL K2/16.4967/,K3/24.3685/,K4/24.745/,K5/19.0411/,N2/1.9185/ &,N3/1.8858/,N4/1.9185/,N5/1.8611/ REAL QO2/0.80/,QO3/0.45/,QO4/0.80/,QO5/0.30/ C UNKNOWN: DQ1 F(1)=K2*(QO2+X(1))**N2+K3*(QO3+X(1))**N3- &K5*(QO5-X(1))**N5-K4*(QO4-X(1))**N4 RETURN END Writing every equation, as was done in the listing of EQUSOL1, makes it easy to fol- low the computational sequence in subroutine FUNCT. However, since there are as many Darcy-Weisbach equations as there are pipes and the Gauss-Seidel iteration could be used to solve the Colebrook-White equations internally within the system of equations, FUNCT can be simplified. A separate function can be written to evaluate the Colebrook-White equation, and the Darcy-Weisbach equations are in a DO loop. Now the equations to be solved are the four junction continuity equations and the six head loss equations that indicate the difference in head along a pipe is equal to the frictional head loss between the pipe ends. For variety, Q1 is now unknown, and in its place Q6 is assumed to be known. The listing of the modified subroutine FUNCT follows: SUBROUTINE FUNCT(X,F) INTEGER*2 ID(6)/7,8,10,9,10,11/ ! 1 = Q1, 10 = H4, 19 = D3 &,IU(6)/16,7,8,7,9,10/ ! 2 = Q2, 11 = H5, 20 = D4 REAL F(10),X(28) ! 3 = Q3, 12 = QJ1, 21 = D5 DATA G2/64.4/,P4/0.7853982/ ! 4 = Q4, 13 = QJ2, 22 = D6 F(1)=X(1)-X(2)-X(4)-X(12) ! 5 = Q5, 14 = QJ3, 23 = L1 F(2)=X(2)-X(3)-X(13) ! 6 = Q6, 15 = QJ4, 24 = L2 F(3)=X(4)-X(5)-X(14) ! 7 = H1, 16 = WS1, 25 = L3 F(4)=X(3)+X(5)-X(6)-X(15) ! 8 = H2, 17 = D1, 26 = L4 DO 10 I=1,6 ! 9 = H3, 18 = D2, 27 = L5 J=I+16 ! 28 = L6 10 F(I+4)=X(ID(I))-X(IU(I))+FR(I,J,X)*X(I+22) & /X(J)*(X(I)/(P4*X(J)**2))**2/G2 RETURN END * FUNCTION FR(I,J,X) REAL X(28) REAL FI(6)/6*.02/ DATA E/0.0004166667/,CCVISC/1.03543E-4/ F1=1./SQRT(FI(I)) 10 F2=F1 F1=1.14-2.*ALOG10(E/X(J)+CCVISC*X(J)*F2/X(I)) IF(ABS(F1-F2).GT. 1.E-6) GO TO 10 FR=1./F1/F1 FI(I)=FR RETURN END The input data file (all in ft or ft3/s) and the solution are given next: Input Data Solution Variable Type Initial value Index Value Q1 U 2.10 1 Q1 2.000 Q2 U 0.80 2 Q2 0.770 Q3 U 0.50 3 Q3 0.420 Q4 U 0.80 4 Q4 0.730 Q5 U 0.30 5 Q5 0.230 Q6 K 0.25 6 Q6 0.250 H1 U 476.0 7 H1 478.2 H2 U 465.0 8 H2 468.2 H3 U 460.0 9 H3 464.7 Continued: Input Data Solution Variable Type Initial value Index Value H4 U 458.0 10 H4 463.5 H5 U 450.0 11 H5 454.3 QJ1 K 0.50 12 QJ1 0.50 QJ2 K 0.35 13 QJ2 0.35 QJ3 K 0.50 14 QJ3 0.50 QJ4 K 0.50 15 QJ4 0.50 WS1 K 500.0 16 WS1 500.0 D1 K 0.667 17 D1 0.667 D2 K 0.50 18 D2 0.50 D3 K 0.50 19 D3 0.50 D4 K 0.50 20 D4 0.50 D5 K 0.50 21 D5 0.50 D6 K 0.333 22 D6 0.333 L1 K 1500 23 L1 1500 L2 K 1000 24 L2 1000 L3 K 1500 25 L3 1500 L4 K 1500 26 L4 1500 L5 K 1200 27 L5 1200 L6 K 1000 28 L6 1000 * * * 4.4.2. SOLVING THE THREE EQUATION SYSTEMS VIA NEWTON The Newton method will now be applied in turn to the solution of the Q-equations, the H-equations and the ∆Q-equations for network shown in Fig. 4.24. Considerable detail will be presented in these solutions so the details of applying the Newton method can be [1] (1) [3] [2] (5) (4) (2) (3) 0.8 ft3/s 100' 90' 1.5 ft3/s 1.0 ft3/s I II ∆Q1 ∆Q2 Pipe K n 1 7.59 1.936 2 9.63 1.901 3 48.6 1.882 4 39.7 1.768 5 16.5 1.935 Figure 4.24 A 5-pipe, 3-node network. examined. The reader is encouraged to check numerically some of these steps. In the Q- equations the elements of the Jacobian will either be ∂Fi / ∂Qj = ±1 or zero in row i for a junction continuity equation row. The Jacobian terms for the energy loop equation rows will either be ∂Fi / ∂Qj = ± njK jQj n j −1 or zero. The non-zero elements of the Jacobian for the H-equations are ∂Fi / ∂H j = ± {1 / (nmKm )}{(H j − Hk ) / Km } 1/nm −1 in which the sign is determined by the sign in front of this term in the equation and the sign before Hj within the parentheses. Non-zero terms in the Jacobian for the ∆Q-equations will be of the form ∂Fi / ∂∆Qj = ± nk Kk (Qok ± ∆Qm ∑ )nk −1. The Q-equations are F1 = Q1 − Q2 − Q4 −1.0 = 0 F2 = Q2 + Q3 −1.5 = 0 F3 = Q4 − Q3 + Q5 − 0.8 = 0 F4 = K1Q1 n1 + K4Q4 n4 − K5Q5 n5 −10 = 0 F5 = K2Q2 n2 − K3Q3 n3 − K4Q4 n4 = 0 (4.47) The Newton method is described by [D]{z}={F} and {Q}(m+1) = {Q}(m) - {z} with D [ ] = 1 −1 0 −1 0 0 1 1 0 0 0 0 −1 1 1 n1K1Q1 n1−1 0 0 n4K4Q4 n4 −1 − n5K5Q5 n5−1 0 n2K2Q2 n2 −1 − n3K3Q3 n3−1 − n4K4Q4 n4 −1 0 (4.48) If we choose the initial estimate of the solution vector to be Q { }(0) = 2.0 0.9 0.6 0.1 1.3 (4.49) with which we have been careful to satisfy the junction continuity equations, so these discharges can be used in the ∆Q-equations, the first evaluation of the Jacobian matrix and right-hand side leads to 1.0000 − 1.0000 0.0000 − 1.0000 0.0000 0.0000 1.0000 1.0000 0.0000 0.0000 0.0000 0.0000 −1.0000 1.0000 1.0000 28.1133 0.0000 0.0000 11.9749 − 40.8039 0.0000 16.6487 − 58.2888 − 11.9749 0.0000 z1 z2 z3 z4 z5 = 0.0000 0.0000 0.0000 − 7.6935 − 11.3783 (4.50a) with the solution z { } = − 0.1169 − 0.1470 0.1470 0.0301 0.1169 Q { }(1) = 2.1169 1.0470 0.4530 0.0699 1.1831 (4.51a) The Newton equations for the next cycle are 1.0000 − 1.0000 0.0000 − 1.0000 0.0000 0.0000 1.0000 1.0000 0.0000 0.0000 0.0000 0.0000 −1.0000 1.0000 1.0000 29.6481 0.0000 0.0000 9.0910 − 37.3636 0.0000 19.0804 − 45.4902 − 9.0910 0.0000 z1 z2 z3 z4 z5 = 0.0000 0.0000 0.0000 − 0.0682 − 0.7993 (4.50b) with the solution z { } = − 0.0022 − 0.0111 0.0111 0.0089 0.0022 Q { }(2) = 2.1191 1.0581 0.4419 0.0610 1.1809 (4.51b) One more cycle would yield the final solution Q { }(3) = 2.1191 1.0583 0.4417 0.0608 1.1809 (4.51c) Referring again to Fig. 4.24, since we have only three nodes, we must construct three H-equations. They are F1 = 100 − H1 K1 1/n1 − H1 − H2 K2 1/n2 − H1 − H3 K4 1/n4 −1.0 = 0 F2 = H1 − H2 K2 1/n2 + H3 − H2 K3 1/n3 −1.5 = 0 F3 = H1 − H3 K4 1/n4 − H3 − H2 K3 1/n3 + 90 − H3 K5 1/n5 − 0.8 = 0 (4.52) Using [D]{z} = {F} and {H}(m+1) = {H}(m) - {z} to implement the Newton method with an initial estimate of the nodal heads as H { }(0) = 93 85 88 (4.53) successive computational cycles produce − 0.166 0.060 0.035 0.060 − 0.100 0.040 0.035 0.040 − 0.162 z1 z2 z3 = −1.258 − 0.365 − 0.382 z { } = 15.953 17.241 10.088 H { }(1) = 77.047 67.759 77.912 (4.54a) − 0.171 0.056 0.075 0.056 − 0.078 0.023 0.075 0.023 − 0.134 z1 z2 z3 = − 0.095 − 0.084 − 0.499 z { } = 8.019 9.614 9.830 H { }(2) = 69.028 58.145 68.083 (4.54b) − 0.158 0.052 0.072 0.052 − 0.075 0.023 0.072 0.023 − 0.123 z1 z2 z3 = − 0.120 − 0.003 − 0.049 z { } = 1.547 1.353 0.770 H { }(3) = 67.480 56.793 67.312 (4.54c) − 0.239 0.052 0.153 0.052 − 0.075 0.022 0.153 0.022 − 0.202 z1 z2 z3 = 0.019 0.000 − 0.019 z { } = −0.032 0.002 0.070 H { }(4) = 67.513 56.791 67.242 (4.54d) − 0.210 0.052 0.124 0.052 − 0.075 0.023 0.124 0.023 − 0.174 z1 z2 z3 = 0.002 − 0.001 0.002 z { } = − − 0.005 0.001 0.006 H { }(5) = 67.517 56.793 67.236 (4.54e) Depending on the desired accuracy of this solution, the process might have been terminated one or two cycles sooner. One of the best features of the Newton method is a quadratic convergence rate as the solution is approached; here we can easily see the rapid reduction in the size of the corrections {z} in successive cycles. When this problem is solved by using the ∆Q-equations, we need only one energy loop and one pseudo loop, as is indicated on Fig. 4.24. The ∆Q-equations are F1 = K1 Qo1 + ∆Q1 ( )n1 + K4 Qo4 − ∆Q2 + ∆Q1 ( )n4 − K5 Qo5 − ∆Q1 ( )n5 −10 = 0 F2 = K2 Qo2 + ∆Q2 ( )n2 − K3 Qo3 − ∆Q2 ( )n3 − K4 Qo4 − ∆Q2 + ∆Q1 ( )n4 = 0 (4.55) The equations for the Newton method are [D]{z}={F} and {∆Q}(m+1) = {∆Q}(m) - {z} in which the Jacobian and vector of initial discharges are D [ ] = ∂F1 ∂∆Q1 ∂F1 ∂∆Q2 ∂F2 ∂∆Q1 ∂F2 ∂∆Q2 Qo { }(0) = 2.0 0.9 0.6 0.1 1.3 (4.56a,b) Three successive solution cycles produce 80.892 −11.975 −11.975 86.913 z1 z2 = − 7.694 −11.378 z { } = − 0.117 − 0.147 ∆Q { }(1) = 0.117 0.147 (4.57a) 76.103 − 9.091 − 9.091 73.662 z1 z2 = − 0.068 − 0.799 z { } = − 0.002 − 0.011 ∆Q { }(2) = 0.119 0.158 (4.57b) 75.163 −8.188 −8.188 71.954 z1 z2 = 0.004 − 0.009 z { } = 0.000 0.000 ∆Q { }(3) = 0.119 0.158 (4.57c) From these results we can compute the discharges themselves in ft3/s as Q1 = Qo1 + ∆Q1 = 2.119 Q2 = Qo2 + ∆Q2 = 1.058 Q3 = Qo3 − ∆Q2 = 0.442 Q4 = Qo4 + ∆Q1 − ∆Q2 = 0.061 Q5 = Qo5 − ∆Q1 = 1.181 (4.58) Readers will find it instructive to solve this same problem by modifying subroutine FUNCT in program EQUSOL1. For still more experience the reader should consider the use of an equation-solving software program such as MathCAD, TK-Solver, or MathLab. 4.4.3. COMPUTER SOLUTIONS TO NETWORKS In this section we concentrate on the implementation of solutions to networks using computers, and how pumps, local losses and/or PRV's are readily included. To begin this process consider first the eight-pipe network in Fig. 4.25 that includes a source pump that supplies some of the system demand and a booster pump in pipe 1. There are also local loss devices in pipes 7, 8, and 3, the first two of which have a loss coefficient of 10, and the third has a loss coefficient of 2. (The Roman numeral loop notation will be used in Example Problem 4.6.) After developing and solving the equations for this network, we will place a PRV in pipe 5. P1 P2 [1] (1) [3] [2] (5) (6) (4) (2) (3) (7) (8) 0.05 m3/s 200 m 170 m [5] 0.25 - 300 III I II 0.25 - 300 Globe valve Globe valve Meter KL3 = 2 KL1 = 10 KL2 = 10 0.03 m3/s 0.08 m3/s 0.08 m3/s 0.2 - 500 0.2 - 500 0.2-3000.2-3000.2-6000.2 - 500 [4] Pump Characteristics (Q in m3/s and h p in meters) Pump No. Point 1 Point 2 Point 3 Q h p Q h p Q h p 1 0.025 12.0 0.040 10.5 0.055 8.0 2 0.060 4.0 0.090 3.8 0.120 3.5 Figure 4.25 An eight-pipe network with pumps and local losses. For this network there are five junction continuity equations and three energy loop equa- tions. The Q-equations are F1 = −Q1 + Q4 + Q7 − 0.03 = 0 F2 = Q1 + Q2 − Q5 − 0.08 = 0 F3 = −Q2 + Q3 − Q6 − 0.05 = 0 F4 = −Q3 − Q4 + Q8 − 0.00 = 0 F5 = Q5 + Q6 − 0.08 = 0 F6 = K1Q1 n1 − hp2 − K2Q2 n2 − K3Q3 n3 − 2Q3 2 / (2gA3 2 ) + K4Q4 n4 = 0 F7 = K5Q5 n5 − K6Q6 n6 + K2Q2 n2 = 0 F8 = K7Q7 n7 − hp1 +10Q7 2 / (2gA7 2 ) − K4Q4 n4 − K8Q8 n8 −10Q8 2 / (2gA8 2 ) + 30 = 0 (4.59) In these equations the local loss coefficients 10, 10, and 2 have been inserted in the equations, but the pump heads are written as hp1 and hp2. These pump heads can be expressed as a function of discharge by fitting a second-order polynomial through three points on the pump characteristic curve over the range of expected operation. Alternatively, if the power supplied by the pump to the fluid is assumed to be constant, then these pump heads can be defined by hp = Power/(γQ). Example Problem 4.6 Solve the 8-equation system, Eqs. 4.59, for the network shown in Fig. 4.25 by using hand methods. Then verify this solution by using program EQUSOL1 and/or an equation- solving software package such as MathCAD or TK-Solver. This and the next few Example Problems will be solved by rewriting the subroutine FUNCT for use with EQUSOL1 in each case. MathCAD and TK-Solver models of these problems will be found on the CD that accompanies this book. The first step is to estimate the pipe discharges; based on these values, we then compute K and n for the 8 pipes. The listed discharges produce the K and n values in the table: Pipe No. 1 2 3 4 5 6 7 8 Q (m3/s) 0.100 0.015 0.100 0.080 0.030 0.050 0.070 0.170 K 1160 613 1160 690 1292 1115 322 239 n 1.827 1.788 1.827 1.824 1.801 1.812 1.772 1.832 f 0.0134 0.0314 0.0134 0.0212 0.0168 0.0152 0.0223 0.0127 The next step is to fit the three pairs of points for the two pumps to obtain the coefficients for the polynomials: A1 = - 2220, B1 = 44.4, C1 = 12.28 and A2 = - 55.6, B2 = 1.667, C2 = 4.10. These values can now be substituted into the equations, the equations can be differentiated to produce the Jacobian matrix and equation vector, and all of the terms can be evaluated by using the data in the table and figure. (The reader should evaluate some terms to verify that the process is fully understood.) The results are D [ ] = −1.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 1.0 1.0 0.0 0.0 −1.0 0.0 0.0 0.0 0.0 −1.0 1.0 0.0 0.0 −1.0 0.0 0.0 0.0 0.0 −1.0 −1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 325 −40.2 −337 158 0.0 0.0 0.0 0.0 0.0 40.2 0.0 0.0 140 −178 0.0 0.0 0.0 0.0 0.0 −158 0.0 0.0 370 −173 and {F}T = {0.020 0.005 - 0.015 -0.010 0.000 1.560 -2.226 4.903} with the superscript T indicating the transpose of the right-side equation vector. The solution to this linear system of equations produces {z}T = {- 0.0029 0.0008 - 0.0071 0.0021 - 0.0071 0.0071 0.0150 - 0.0150} so that the discharges after the first iteration are {Q}T = {0.1029 0.0142 0.1071 0.0779 0.0371 0.0429 0.0550 0.1850} After two additional iterations, the following solution, in m3/s, is obtained: {Q}T = {0.1028 0.0142 0.1072 0.0780 0.0370 0.0430 0.0548 0.1852} The subroutine FUNCT is on the CD under EPRB4_6Q.FOR for use with EQUSOL1, and the TK-Solver model is listed on the CD under EPRB4_6Q.TK. Upon supplying the input file in column one below to EQUSOL1, the solution in the second column is obtained: From the keyboard: 8 8 2 3 Solution (m3/s) Q1 U 0.100 1 Q1 = 0.103 Q2 U 0.015 2 Q2 = 0.014 Q3 U 0.100 3 Q3 = 0.107 Q4 U 0.080 4 Q4 = 0.078 Q5 U 0.030 5 Q5 = 0.037 Q6 U 0.050 6 Q6 = 0.043 Q7 U 0.070 7 Q7 = 0.055 Q8 U 0.170 8 Q8 = 0.185 * * * Next let us examine the formulation and the solution of the H-equations for the network in Example Problem 4.6. The H-equations are presented as Eqs. 4.60. The pump heads have been added to the upstream HGL-elevations by indicating these heads as hp. In a similar way the local losses, denoted simply as hL, have been subtracted from the HGL-elevations in pipes 3, 7, and 8. By using second-order polynomials for the pump characteristics and noting that each Q in these equations can again be replaced by similar F1 = H1 + hp2 − H2 K1 1/n1 − H4 − H1 K4 1/n4 − 170 + hp1 − hL1 − H1 K7 1/n7 + 0.03 = 0 F2 = H2 − H5 K5 1/n5 − H3 − H2 K2 1/n2 − H1 + hp2 − H2 K1 1/n1 + 0.08 = 0 F3 = H3 − H5 K6 1/n6 + H3 − H2 K2 1/n2 − H4 − hL3 − H3 K3 1/n3 + 0.05 = 0 F4 = H4 − hL3 − H3 K3 1/n3 + H4 − H1 K4 1/n4 − 200 − hL2 − H4 K8 1/n8 + 0.00 = 0 F5 = − H2 − H5 K5 1/n5 − H3 − H5 K6 1/n6 + 0.08 = 0 (4.60) head-difference terms, we find that we cannot create an equation that does not contain the pump discharge in the pipes that contain the pumps. The same is true for pipes that contain local losses because again the magnitude of the loss is a function of the discharge through the pipe. If an iterative approach were used to approximate the discharge in terms of the upstream and downstream nodal heads, the quadratic convergence rate of the Newton method would be sacrificed. This dilemma highlights a deficiency in using the H- equations when pumps, whose heads are strongly dependent upon the discharges passing through them, are present. If many pumps exist in a network and the H-equations are to be used, especially if equation-solving software such as MathCAD or TK-Solver is used, it might be advantageous to write the continuity equations in terms of the discharges, then add additional equations to describe these discharges in terms of nodal heads, and finally solve simultaneously for the heads and discharges. This approach will be taken in Example Problem 4.7. A successful technique that solves the H-equations in a computer program will be described in a subsequent section. This technique obtains the current value of the discharge in every pipe from the heads that exist during an iteration by calculating Q = [(Hi - Hj)/K]1/n, and when a pump exists in a pipe, the Newton method can find the discharge Q from the single equation for that pipe Q = [(Hi + hp - Hj)/K]1/n in which the pump head is hp = AQ2 + BQ + C. Example Problem 4.7 Solve the H-equations for the 8-pipe network shown in Fig. 4.25 that was the subject of study in Example Problem 4.6. The form of the subroutine FUNCT that is needed in EQUSOL1 to solve the H- equations is on the CD in EPRB4_7H.FOR with input data in EPRB4_7H.DAT; the corresponding TK-Solver model will be found there as EPRB4_7H.TK. * * * Finally, we now turn our attention to the ∆Q-equations, which are given in Eq. 4.61 for this same pipe network. In addition, we must select an appropriate set of initial discharges. F1 = K1 Qo1 + ∆Q1 ( )n1− hp2 − K2 Qo2 − ∆Q1 + ∆Q2 ( )n2 − K3 Qo3 − ∆Q1 ( )n3 − 2 Qo3 − ∆Q1 ( )2 / 2gA3 2 ( ) + K4 Qo4 + ∆Q1 − ∆Q3 ( )n4 = 0 F2 = K5 Qo5 + ∆Q2 ( )n5 − K6 Qo6 − ∆Q2 ( )n6 + K2 Qo2 − ∆Q1 + ∆Q2 ( )n2 = 0 F3 = K7 Qo7 + ∆Q3 ( )n7 +10 Qo7 + ∆Q3 ( )2/ 2gA7 2 ( ) − K4 Qo4 + ∆Q1 − ∆Q3 ( )n4 − hp1 − K8 Qo8 − ∆Q3 ( )n8 −10 Qo8 − ∆Q3 ( )2/ 2gA8 2 ( ) + 30 = 0 (4.61) In these ∆Q-equations the pump heads will be replaced by second-order polynomial equations in the forms hp1 = A1(Qo7 + ∆Q3 ) 2 + B1(Qo7 + ∆Q3 ) + C1 hp2 = A2 (Qo1 + ∆Q1) 2 + B2 (Qo1 + ∆Q1) + C2 (4.62) The local losses are replaced by hL = KLQ2/(2gA2), in which each discharge is written as the algebraic sum of Qoi and the corrective discharges in that pipe. The derivatives of these terms that contribute to elements of the Jacobian are then easily evaluated. Example Problem 4.8 Solve the ∆Q-equations for the 8-pipe network depicted in Fig. 4.25. To solve the ∆Q-equations using EQUSOL1, download EPRB4_8D.FOR from the CD; the appropriate input is found in EPRB4_8D.DAT. A TK-Solver model will be found as EPRB4_8D.TK. A set of initial discharges might be selected as follows: Qo1 = 0.100 m3/s, Qo2 = 0.015 m3/s, Qo3 = 0.110 m3/s, Qo4 = 0.060 m3/s, Qo5 = 0.035 m3/s, Qo6 = 0.045 m3/s, Qo7 = 0.070 m3/s, and Qo8 = 0.170 m3/s. * * * 4.4.4. INCLUDING PRESSURE REDUCING VALVES To acquire experience in analyzing networks containing pressure reducing valves (and similar appurtenances such as back pressure valves), let us assume that pipe 5 in our 8- pipe network contains a PRV that is 200 m downstream from junction 2, and the pres- sure setting of this valve is equivalent to a reservoir water surface elevation of 149 m. The five junction continuity equations are unchanged for this network. The energy equations now consist of one real loop and two pseudo loops, as the revised diagram of this network shows in Fig. 4.26. The pump data in Fig. 4.25 are unchanged. In the P1 P2 [1] (1) [3] [2] (5) (6) (4) (2) (3) (7) (8) 0.05 m3/s 200 m 170 m [5] III I 0.03 m3/s 0.08 m3/s 0.08 m3/s II [4] (HGL)1 (PRV)1 Figure 4.26 An eight-pipe network with pumps and local losses, now including a PRV. Q-equations, Eqs. 4.72, equations F6 and F8 are unchanged, but equation F7 must now be written around the new pseudo loop; it becomes F7 = K5 ' Q5 n5 − K6Q6 n6 − K3Q3 n3 − 2Q3 2 / (2gA3 2 ) − K8Q8 n8 −10Q8 2 / (2gA8 2 ) + 200 − HGL1 = 0 (4.63) in which K5 ' represents the portion of pipe 5 downstream from the PRV. Eight inde- pendent equations therefore exist, which may be solved for the discharges Qi, i = 1,8. If the PRV does not close completely, the solution is obtained with HGL1 in F7 equal to the head that is set at the valve, i.e., 149 m in this example. If the Newton solution process that is based on this assumption should produce a negative value for Q5, the PRV will close completely. Then the discharge in the pipe containing the PRV is no longer an unknown but is zero, i.e., Q5 = 0, and the value of the HGL immediately downstream from the PRV becomes the unknown. In other words, when this PRV closes, the same system of equations still describes the network operation, but the set of unknowns changes to {Q1, Q2, Q3, Q4, HGL1, Q6, Q7, Q8}. In a computer program this change can be accommodated by dropping the pipe number containing the PRV from the integer arrays that define the junction continuity equations and the energy loop equations. Also, a flag is set in the program to indicate in the solution array that the HGL of the PRV is stored in place of the discharge in that pipe. Subroutine FUNCT for use with EQUSOL1 for this problem is on the CD under the name EPRB4_VQ.FOR. The reader should study a listing of this file to understand the logic that will model the PRV when it closes. The input from the keyboard for this problem is 8 8 2 3 (Thus there are 8 unknowns and 8 variables, and the 2 and 3 are the input and output units.). The input file that was used in Example Problem 4.6 still applies. There are only two basic changes in the program that was used to solve this network problem in the absence of the PRV: (1) element X(5) is now either Q5 or HGL, depending upon whether it is negative, and (2) now F(7) has become a pseudo loop from the artificial reservoir at the PRV to the reservoir at the end of pipe 5. The computer provides the following solution (units are m3/s or m): Q1 = 0.102, Q2 = - 0.022, Q3 = 0.108, Q4 = 0.077, Q5 = 149.227, Q6 = 0.080, Q7 = 0.055, Q8 = 0.185. From this solution we see that the flow has tried to reverse direction in pipe 5, indicating that the PRV must then close. Thus Q5 = 0, and the reported value of 149.227 is the HGL on the downstream side of the PRV, which is slightly above its pressure setting. In- stead of operating in its normal mode, the PRV has acted as a check valve, not permitting the flow to reverse its direction. Let's increase the demand at node 5 to 0.100 m3/s. To obtain a solution for this de- mand, we must change the line that defines the continuity equation at node 5 in subroutine FUNCT to F(5)=Q5+X(6)-0.100. Now the execution of the program produces the following solution (units are m3/s): Q1 = 0.113, Q2 = 0.002, Q3 = 0.117, Q4 = 0.079 Q5 = 0.034, Q6 = 0.066, Q7 = 0.063, Q8 = 0.197 With this slightly larger demand at node 5 the PRV operates normally, maintaining an HGL = 149 m on its downstream side. From these discharges the pipe head losses, the pump heads, and other quantities can be evaluated, and from these the head at each node can be determined. Upon carrying out such computations, we find that the heads are H1 = 173.6 m, H2 = 155.7 m, H3 = 155.7 m, H4 = 180.3 m, and H5 = 147.2 m. The head on the upstream side of the PRV is 154.8 m, so the PRV loss is 154.8 - 149 = 5.8 m. If the demand at node 5 is QJ5 = 0.16 m3/s, the continuity statement in FUNCT for node 5 is changed to F(5)=Q5+X(6)-0.16, then the solution (units are m3/s) becomes Q1 = 0.145, Q2 = 0.050, Q3 = 0.145, Q4 = 0.088 Q5 = 0.115, Q6 = 0.045, Q7 = 0.087, Q8 = 0.233 Without a more complete examination of these solutions, we might be inclined to accept all of them as valid. However, upon computing some head losses and pump heads, the following are found: hf7 = 6.98 m, hp1 = - 28.0 m, hp2 = 4.04 m, and hf1 = 4.92 m, so the HGL at node 2 is H2 = 134.2 m. Obviously the negative head for pump 1 is unrealistic; it is caused by operating the pump with a discharge that is outside the range of the three pairs of points that were used to define this pump characteristic curve. Also H2 = 134.2 m is much smaller than the HGL setting of the PRV, and since this device can not act as a pump to increase the head across it, the most it can do is open completely. If a PRV opens completely, it typically still acts as a local loss device in a way that is similar to a globe valve with a loss coefficient of about 10. Thus the last solution is not valid, and the problem must be solved again with another local loss device in place of the PRV. Currently there exists no simple a priori test to learn that the PRV should open fully and that it is unable to maintain its pressure setting; a solution must first be obtained when we use the Q-equations, because the nodal heads are determined as a secondary step after the discharges are found. The same statements apply to the use of the ∆Q-equations. The three modes in which a PRV may operate are treated most conveniently with the H-equations, since it is then possible to check directly, as the solution is obtained, whether the head at the upstream side of the PRV is less than Hv (the HGL upstream of the PRV) and whether the head at the downstream pipe node is greater than the set HGL. If the head H at the downstream node becomes larger than the HGL setting of the PRV, then it should shut off the flow in the pipe, and if Hv becomes smaller than HGL, then the PRV should open fully. When the ∆Q-equations are used to analyze networks that contain PRV's, we must work with two different sets of loops, one around which the ∆Q's circulate, and one around which the energy equations are written. For our 8-pipe example, Fig. 4.27, the first and P1 P2 [1] (1) [3] [2] (5) (6) (4) (2) (3) (7) (8) 0.05 m3/s 200 m 170 m [5] III I 0.03 m3/s 0.08 m3/s 0.08 m3/s II ∆Q1 ∆Q2 ∆Q3 [4] (HGL)1 (PRV)1 Figure 4.27 The eight-pipe network with a PRV, with loop notation shown. third equations match the corresponding equations, Eqs. 4.61, for this network without a PRV in pipe 5. The second equation is replaced by F2 = K5 ' (Qo5 + ∆Q2 ) n5 − K6 (Qo6 − ∆Q2 ) n6 − K3(Qo3 − ∆Q1) n3 − 2(Qo3 − ∆Q1) 2 / (2gA3 2 ) − K8 (Qo8 − ∆Q3 ) n8 −10(Qo8 − ∆Q3 ) 2 / (2gA8 2 ) − HGL + 200 = 0 (4.64) It is notable that this equation does not contain ∆Q2 in every term and that the system of equations does not produce a symmetric Jacobian. To determine the correct operational mode for a PRV when using the ∆Q-equations is much the same as with the Q-equations. Should the flow in a pipe reverse direction, then the PRV should close, and if the HGL at the upstream end of the PRV is less than its setting, then the PRV should open fully; otherwise the PRV is operating in its normal mode. Logic can easily be included in the computer program to check whether the flow is negative in pipes containing PRV's and then change the nature of the problem being solved. The fully-open mode of operation can not be determined until the nodal heads are computed, as with the Q-equations. Should a PRV close, then the discharge in that pipe becomes zero, and the HGL becomes unknown and larger than the setting. If a pipe containing a PRV has only one ∆Q flowing through it, then that corrective discharge becomes known and is ∆Qj = ± Qoi, in which the minus sign applies if the assumed directions for Qoi and ∆Qj coincide, and the plus sign applies if these directions are opposed. In place of ∆Qj as the unknown, the HGL is unknown, and the number of unknowns still equals the number of equations. Should a PRV close that is internal, with two or more corrective discharges circulating through it, then one of these corrective discharges must be expressed in terms of the others, and the HGL of the PRV replaces this ∆Q in the list of unknowns. In our example, if the PRV were in pipe 2 instead of pipe 5, as shown in Fig. 4.28, then ∆Q2 = Qo2 + ∆Q1, and it is replaced by this quantity wherever else ∆Q2 appears, such as in the discharges for pipes 6 and 7. (1) (5) (6) (4) (2) (3) I II ∆Q1 (HGL)1 ∆Q2 (7) (1) (5) (6) (4) (2) (3) (7) H1 (HGL)1 H1 Corrective discharge loops Energy equation loops Figure 4.28 The modeling of a PRV in pipe 2 of the 8-pipe network. To study this problem further, the reader should obtain a listing of FUNCT under the name EPRB4_9.FOR. It can be used to solve this problem. One additional integer variable IDOO has been added to the list of arguments in FUNCT; it is given a value of 0 in the calling program when the equation vector is evaluated and 1 when the Jacobian is evaluated. This variable is needed because we do not want to close the PRV when we evaluate the derivatives. It is instructive to trace the logic that sets Q5 to zero when the PRV is closed, fixes the value of ∆Q2 = Qo5 and replaces ∆Q2 by the HGL as the unknown represented by X(2). These changes are made when Q5 becomes negative. Subsequent checks might determine whether the HGL becomes less than the PRV setting; if this occurs, the PRV should be reopened. Another modification of this subroutine allows the initial discharges Qoi, i = 1-8, to enter FUNCT through X(i), i = 4-11, thus making it possible to change the demands without changing Qoi within the subroutine. Example Problem 4.9 Solve the eight-pipe network shown in Fig. 4.27 by using the ∆Q-equations. Obtain this solution first for a demand at node 5 of QJ5 = 0.100 m3/s and then for QJ5 = 0.080 m3/s. The input data (EPRB4_9.DAT) to solve this problem with a demand of 0.100 m3/s at node 5 is listed below with the solution: Input Data Solution DQ1 U 0.00 1 DQ1 = - 0.00749 DQ2 U 0.00 2 DQ2 = - 0.00570 DQ3 U 0.00 3 DQ3 = - 0.01668 Qo1 K 0.12 4 Qo1 = 0.12 Qo2 K 0.00 5 Qo2 = 0.00 Qo3 K 0.11 6 Qo3 = 0.11 Qo4 K 0.07 7 Qo4 = 0.07 Qo5 K 0.04 8 Qo5 = 0.04 Qo6 K 0.06 9 Qo6 = 0.06 Qo7 K 0.08 10 Qo7 = 0.08 Qo8 K 0.18 11 Qo8 = 0.18 Applying these solution values for the ∆Qi, as appropriate, to the initial discharges gives the final discharges, and then the pipe head losses can be computed, using the proper K and n for each pipe, as listed in the table: Pipe 1 2 3 4 5 6 7 8 Q, m3/s 0.1125 - 0.0018 0.1175 0.0792 0.0343 0.0657 0.0633 0.1967 hL, m 21.40 0.0075 23.20 6.77 2.97 8.03 2.43 12.14 From these discharges the pump heads and local losses are hp1 = 6.18 m, hp2 = 3.58 m, hL1 = 0.848 m, hL2 = 8.18 m, and hL3 = 1.426 m. From these the nodal heads can be found as H1 = 175.3 m, H2 = 157.5 m, H3 = 155.1 m, H4 = 179.7 m, H5 = 147.0 m, and Hv1 = 156.5 m. We see that the head upstream from the PRV is 156.5 m which is less than H2 = 157.5 m, so the PRV has not opened fully. The head at node 5 down- stream, H5 = 147.0 m, is less than the HGL setting of the PRV (149 m) so the PRV has not closed but operates normally. When the demand at node 5 is QJ5 = 0.080 m3/s, then the input data and solution are Input Data Solution DQ1 U 0.00 1 DQ1 = - 0.01827 DQ2 U 0.00 2 DQ2 = 149.2 DQ3 U 0.00 3 DQ3 = - 0.02508 Qo1 K 0.12 4 Qo1 = 0.12 Qo2 K 0.00 5 Qo2 = 0.00 Qo3 K 0.09 6 Qo3 = 0.09 Qo4 K 0.07 7 Qo4 = 0.07 Qo5 K 0.04 8 Qo5 = 0.04 Qo6 K 0.04 9 Qo6 = 0.04 Qo7 K 0.08 10 Qo7 = 0.08 Qo8 K 0.16 11 Qo8 = 0.16 In this input file the initial discharge estimates Qoi have been altered from previous values so that all continuity equations remain satisfied with QJ5 = 0.080 m3/s. The solution values remind us that ∆Q2 is actually the HGL at the downstream end of the PRV since it has closed, and FUNCT has set ∆Q2 = - Qo5 and then used X(2) as the position for HGL. The table lists the discharges and head losses for this situation: Pipe 1 2 3 4 5 6 7 8 Q, m3/s 0.1017 0.0217 0.1083 0.0768 0.0 0.0800 0.0549 0.1851 hL, m 17.83 0.652 19.98 6.40 0.0 11.47 1.887 10.86 The pump heads and local losses are hp1 = 8.02 m, hp2 = 3.69 m, hL1 = 0.638 m, hL2 = 7.25 m and hL3 = 1.211 m, with nodal heads of H1 = 177.4 m, H2 = 163.2 m, H3 = 160.7 m, H4 = 181.9 m, and H5 = 149.2 m. Now the PRV has closed entirely so the flow in pipe 5 is zero, and the HGL at its downstream end is above its set point. * * * 4.4.5. SYSTEMATIC SOLUTION OF THE Q-EQUATIONS In earlier sections we have developed the three systems of equations that can be used to analyze pipe networks. We have written these equation systems for several small networks, seen how the Newton method can be applied to any system of nonlinear equations and how to solve a problem by using a general purpose program that implements the Newton method for all three equation systems, and finally we have carried out detailed computations by hand to obtain some solutions by iteration. In this section we will see how this knowledge can be used in developing computer programs that will analyze any pipe network by using the Q-equations, and the programs will require only enough information from the user to describe adequately the network and its connectivity. In the next two sections similar programs will be developed for the solution of the H-equations and the ∆Q-equations. Let us begin by assuming that there are no local losses. If they exist, they can be mod- eled simply by assigning a larger equivalent sand roughness, or Hazen-Williams CHW, to the pipes containing minor losses. Here we ignore the Manning equation. In describing any network of pipes, we need two types of information: (1) Pipe infor- mation consisting of the diameter, length, roughness coefficient, and end nodes for each pipe. This information can be called pipe-oriented data, since we assemble it by going though a list of the pipes in the network; (2) Junction information, including the demand at the junction, its elevation, and possibly the pipes that join at the junction. This information is called node-oriented data, since it is assembled by moving through the nodes of the network. Actually the connectivity of the network can be defined either by giving the nodes at the ends of each pipe, or by giving the pipe numbers that join at each node. We shall use this duplicative information to check that the user has not erred in defining the network. Now we shall describe the input data that are required. Details on the form of this input will be provided subsequently. Prior to study of this section the reader should obtain a listing of the program SOLQEQS.FOR (or C if you prefer) from the text CD. The information that is required from the user is the following: 1. A line that gives (a) the number of pipes, (b) the number of nodes, (c) the number of reservoirs that supply the network, (d) the number of pumps, and (e) the options which you wish to change from the default values. (The default options and how these are changed will be described later.) 2. For each pipe, list its (a) number, (b) upstream node, (c) downstream node, (d) length, (e) diameter, and (f) roughness coefficient. 3. For each node, list its (a) number, (b) demand, and (c) elevation, and (d) a list of pipes that join this node. 4. For each reservoir, list (a) the pipe number that connects this reservoir to the net- work, and (b) the water surface elevation of the reservoir. 5. For each pump, list (a) the number of the pipe that contains the pump, and (b) three (Q, hp) pairs of discharge vs. pump head that will allow its operating characteristics to be defined. 6. Finally, because the algorithms that could be used to determine the minimum set of independent loops for the energy equations are relatively complex, we require a list of pipe numbers around each loop (with a minus sign before a pipe number if the movement around the loop opposes the assumed direction of flow in that pipe). We require that pseudo loops be provided first, and then the real loops. The information in item one is used to dimension the arrays that will store the remaining information about the pipe network and to determine how many lines of each information type will be read from the input data file. The information must be provided in the sequence that is listed above. The program must perform five major tasks: 1. Read the input data that defines the network. 2. Develop from this information the system of Q-equations, i.e., the junction contin- uity equations and the energy equations around pseudo and real loops of the network. This task defines the equations and also forms each element of the Jacobian matrix. 3. Solve the system of equations. Here we will simply call a standard linear algebra solver. However, a program for larger network problems should have a special linear algebra solver that takes advantage of the special properties of a sparse Jacobian matrix. 4. Obtain the head H at each node after the pipe discharges have been found. 5. Write the solution results in tables that can be readily understood. We choose to provide these results in two tables: a pipe data table and a node data table. In reading the pipe numbers that connect at a node and the pipe numbers that define a loop, a pointer is used to separate data for consecutive nodes, and a second pointer separates data for consecutive loops. The pipes that join at nodes are placed consecutively in a one- dimensional array JN, with NN pointing to the position in this array that separates data for consecutive nodes. A similar one-dimensional array IK contains the pipe numbers that form the loops, with LP pointing to the first pipe number in each loop. The use of one-dimensional arrays with pointers is a more efficient use of storage than the use of two- dimensional arrays, because the second subscript of a two-dimensional array must then be at least as large as the maximum number of pipes that may exist in a loop. When solving the Q-equations (or ∆Q-equations), we compute the nodal heads after obtaining the solution for the discharges. These heads are found by starting at all reservoirs and computing each H at the node at the other end of a pipe from (to) the reservoir by subtracting (adding) the pipe head loss from (to) the reservoir water surface elevation. After these heads have been determined, the nodes one pipe away from these can be determined next, and so on. This process continues until the head at every node has been determined. In program SOLQEQS the computation of heads begins after the PIPE DATA table is written by the DO 130 loop. This loop finds each head at the other end of a pipe that is connected to a reservoir, and upon computing H the integer array INDX, with its argument equal to this node number, is set to 1 to show that nodal head has been computed. Now loop DO 160 passes through the nodes, but nothing is done if INDX(I) for node I is zero. Otherwise INDX(I) = 1, and then the pipes that join this node are accessed through array JN; if H at the other end of a pipe is not known, it is computed. Since not all nodal heads will be found from the first pass through the nodes, the integer IJ also accumulates the number of nodes whose head has been found. One way to learn if another pass is needed is to check whether IJ is less than NJ, the total number of nodes. Actually we see whether IJ has increased from the previous pass. If so, we pass through the nodes again. This method may result in passing through the nodes one extra time, but it prevents the creation of an infinite loop if there is an error in the network description so that fewer than NJ heads can be found. After finding every head, the NODE DATA table is written. The program then allows the user to solve another problem whose data is in a different file, or to change the peaking factor for the same network. Detailed instructions on the preparation of input data to SOLQEQS follows: Line 1: No. of Pipes (NP), No. of Nodes (NJ), No. of Reservoirs (NRES), No. of Pumps (NPUMP), No. of Options (NOPT), Option Pairs. The options consist of a letter in quotes followed by a value, as follows: Letter Controls What Value Default U or u ES or SI units 0 = ES, 1 = SI 0 Q or q Discharge units 0 = ft3/s, 1 = gal/min, 2 = mgd, 3 = m3/s, 4 = l/s 0 or 3 D or d Pipe diameter and roughness units 0 = in, 1 = ft, 2 = m, 3 = cm, 4 = mm 0 if ES 4 if SI F or f Peaking factor Multiplier of demands 1.0 V or v Kinematic viscosity ν = value ES, 1.217E-5 SI, 1.310E-6 G or g Specific weight γ = value ES, 62.4 SI, 9806.0 C or c Network check 1 = yes, 0 = no 1 Here is an example of specifying options: 2 'U' 1 'F' 1.5. The 2 indicates two options are being changed, to specify SI units and to specify a peaking factor of 1.5. In giving the options, the units (ES or SI) option should appear first if it is to be elected, but otherwise the options can be given in any order. Next group: Pipe data consisting of NP lines, each giving pipe number, node 1, node 2, length, diameter, and roughness. Pipes must be numbered consecutively, starting with 1, but they need not be entered consecutively. The roughness may be either the equivalent sand roughness e (in the same units as the diameter) for use in the Colebrook-White and Darcy-Weisbach equations, or a Hazen-Williams CHW, and these may be mixed. The program decides which equation to use, based on the roughness size. Next group: Node data consisting of NJ lines, each giving node number, demand, elevation, number of pipes at the node, and a list of these pipe numbers with a minus if the flow is from the junction. This information is used to define the junction continuity equations. Next group: Reservoir data consisting of NRES lines, each giving the number of the pipe connecting the reservoir to the network and the water surface elevation. Next group: Pump data consisting of NPUMP lines, each with the number of the pipe containing the pump, followed by three (Q, hp) pairs to define the pump curve. Next group: Loop data consisting of NL = NP - NJ lines, one loop on each line with the number of pipes in the loop and a list of these pipes. A negative sign must precede the pipe number if the direction around the loop opposes the assumed direction of flow in this pipe. Pseudo loops connecting reservoirs must appear first in this list, and the real loops follow. Example Problem 4.10 Use program SOLQEQS to solve the network of Example Problem 4.5. Obtain two solutions: (1) for the given demands and (2) with these demands multiplied by 2.0. The input data takes the form 6 5 1 0 1 'D' 1 1 0.50 350. 3 1 - 2 - 4 1 0 1 1500 0.667 0.000417 2 0.35 350. 2 2 - 3 2 1 2 1000 0.5 0.000417 3 0.50 350. 2 4 - 5 3 2 4 1500 0.5 0.000417 4 0.50 350. 3 3 5 - 6 4 1 3 1500 0.5 0.000417 5 0.25 350. 1 6 5 3 4 1200 0.5 0.000417 1 500 6 4 5 1000 0.333 0.000417 4 2 3 - 5 - 4 or, if the diameters and e are given in inches (Inches is the default; either giving 0 options as the last 0 in the first line below, or giving 1 'D' 0, tells SOLQEQS to use inches), then the input would be 6 5 1 0 0 1 0.50 350. 3 1 - 2 - 4 1 0 1 1500 8.0 0.005 2 0.35 350. 2 2 - 3 2 1 2 1000 6.0 0.005 3 0.50 350. 2 4 - 5 3 2 4 1500 6.0 0.005 4 0.50 350. 3 3 5 - 6 4 1 3 1500 6.0 0.005 5 0.25 350. 1 6 5 3 4 1200 6.0 0.005 1 500 6 4 5 1000 4.0 0.005 4 2 3 - 5 - 4 When prompted after the first solution, we give the peaking factor 2.0. The solution tables follow. In the last NODE DATA table we see that some heads are negative, so the network is unable to supply demands that are double the initial values. PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 4 Q VEL. HEAD LOSS HLOSS/ 1 0 0 0 ft. ft ft ft3/s ft/s ft. 1 0 1 1500 0.667 4.17 2.10 6.02 23.50 15.67 2 1 2 1000 0.500 4.17 0.82 4.18 11.00 11.00 3 2 4 1500 0.500 4.17 0.47 2.39 5.67 3.78 4 1 3 1500 0.500 4.17 0.78 3.97 14.97 9.98 5 3 4 1200 0.500 4.17 0.28 1.43 1.70 1.42 6 4 5 1000 0.333 4.17 0.25 2.87 8.83 8.83 NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. Estimate ft3/s ft. ft. lb/in2 ft. 1 0.5 0.500 350.0 126.50 54.82 476.50 2 0.3 0.350 350.0 115.50 50.05 465.50 3 0.5 0.500 350.0 111.53 48.33 461.53 4 0.5 0.500 350.0 109.82 47.59 459.82 5 0.3 0.250 350.0 101.00 43.77 451.00 For peaking factor = 2.0: PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 4 Q VEL. HEAD LOSS HLOSS/ 1 0 0 0 ft. ft ft ft3/s ft/s ft. 1 0 1 1500 0.667 4.17 4.20 12.03 91.53 61.01 2 1 2 1000 0.500 4.17 1.64 8.36 42.48 42.48 3 2 4 1500 0.500 4.17 0.94 4.79 21.53 14.35 4 1 3 1500 0.500 4.17 1.56 7.94 57.69 38.46 5 3 4 1200 0.500 4.17 0.56 2.85 6.32 5.27 6 4 5 1000 0.333 4.17 0.50 5.73 33.66 33.66 NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. Estimate ft3/s ft. ft. lb/in2 ft. 1 1.0 1.000 350.0 58.48 25.34 408.48 2 0.7 0.700 350.0 15.99 6.93 365.99 3 1.0 1.000 350.0 0.79 0.34 350.79 4 1.0 1.000 350.0 - 5.53 - 2.40 344.47 5 0.5 0.500 350.0 - 39.20 - 16.99 310.80 * * * Example Problem 4.11 Use program SOLQEQS to analyze the 5-pipe, 3-node network in the figure. In pipe 1 is a pump, with the characteristics given in the table, which is connected to a reservoir. Let ν = 1.417x10-5 ft2/sec. The elevation of all nodes is zero. P [1] (1) [3] [2] (5) (4) (2) (3) 1.2 ft3/s All e = 0.002" 8" - 4000' 6" - 3000 ' 8" - 6000' 100' 12" - 4000' 1.5 ft3/s 1.0 ft3/s 6" - 2000' 90' The pump curve is described by data in the following table: Q , ft3/s H, ft 4.5 54 4.0 50 3.5 44 The input data this problem are listed first in two columns, followed by the solution tables. 5 3 2 1 1 'V' 1.417E-5 2 1.2 0 2 3 - 4 1 0 1 4000 12 0.002 3 1.0 0 3 2 4 - 5 2 1 3 6000 8 0.002 1 100 3 1 2 4000 8 0.002 5 90 4 2 3 3000 6 0.002 1 4.5 54 4.0 50 3.5 44 5 3 0 2000 6 0.002 3 1 2 5 1 1.5 0 3 1 - 2 - 3 3 2 - 4 - 3 PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 3 Q VEL. HEAD LOSS HLOSS/ 1 0 0 0 ft. in in ft3/s ft/s ft. 1 0 1 4000 12.0 2.0 4.13 5.26 26.24 6.56 2 1 3 6000 8.0 2.0 1.21 3.45 29.13 4.85 3 1 2 4000 8.0 2.0 1.42 4.08 26.48 6.62 4 2 3 3000 6.0 2.0 0.22 1.14 2.64 0.88 5 3 0 2000 6.0 2.0 0.43 2.18 5.85 2.93 Pump 1 in pipe 1: Head = 52.21 ft, Q = 4.13 ft3/s NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. Estimate ft3/s ft. ft. lb/in2 ft. 1 1.5 1.500 0.0 124.98 54.16 124.98 2 1.2 1.200 0.0 98.50 42.68 98.50 3 1.0 1.000 0.0 95.85 41.54 95.85 * * * Example Problem 4.12 Solve the 7-pipe, 4-node network shown in Fig. 4.6, which contains a PRV in pipe 6, by using program SOLQEQS. The input data for this problem are listed after this paragraph. Then the two solution tables follow. Several observations should be made here: In the lines of nodal data the information after the nodal demand that lists the pipes that join at a node is used to define the junction continuity equations; therefore the list of pipes that join at node 1 must include pipe 6 with the PRV in it. The input that lists the pipes that define the loops of the network are used to define the energy loop equations; this group should therefore define a loop that starts (or ends) at the artificial reservoir created by the PRV, so for this network there will be two pseudo loops and one real loop. Also, since the downstream part of pipe 6 defines K', its length is 500 ft, and its upstream node is given as 0 (a reservoir). 7 4 3 1 1 'C' 0 3 0 50 4 4 7 -3 -5 1 0 1 1000 6 0.02 4 20 1 2 5 6 2 1 2 1000 6 0.02 1 90 3 3 2 800 6 0.02 4 100 4 0 3 200 6 0.02 6 55 5 3 4 2000 6 0.02 1 1.0 60 1.5 55 2.0 48 6 0 4 500 6 0.02 3 1 7 - 4 7 1 3 1500 1 0.02 3 6 - 5 - 4 1 0 50 4 1 - 2 - 6 - 7 3 2 - 3 - 7 2 1 50 2 2 3 PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 3 Q VEL. HEAD LOSS HLOSS/ 1 0 0 0 ft. in in ft3/s ft/s ft. 1 0 1 1000 6.0 20.0 1.11 5.65 27.28 27.28 2 1 2 1000 6.0 20.0 1.07 5.43 25.26 25.26 3 3 2 800 6.0 20.0 - 0.07 - 0.34 0.10 0.12 4 0 3 200 6.0 20.0 0.89 4.54 3.53 17.74 5 3 4 2000 6.0 20.0 0.96 4.91 41.47 20.74 6 0 4 500 6.0 20.0 0.04 0.18 0.04 0.04 7 1 3 1500 1.0 20.0 0.01 1.31 25.56 17.04 Pump 1 in pipe 1: Head = 59.09 ft, Q = 1.11 ft3/s NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. Estimate ft3/s ft. ft. lb/in2 ft. 1 0.0 0.000 50.0 71.81 31.1 121.81 2 1.0 1.000 50.0 46.55 20.2 96.55 3 0.0 0.000 50.0 46.45 20.1 96.45 4 1.0 1.000 20.0 34.98 15.2 54.98 * * * 4.4.6. SYSTEMATIC SOLUTION OF THE H -EQUATIONS This section is similar to Section 4.4.5, but now the objective is to describe a program that analyzes a network by solving the H-equations. This program will be restricted to the solution of the H-equations for networks that do not contain a PRV or a BPV and in which minor losses can be neglected. (Exercises to include these devices can be found in the end-of-chapter problems.) With these restrictions the Jacobian matrix of the equation system is symmetric. Symmetry occurs because the partial derivatives of terms which describe the discharge in pipe k between nodes i and j, such as {(Hi − H j ) /Kk } 1/nk , will be the same in the equations where this discharge occurs, so long as neither i nor j are the node for which this junction continuity equation is written. With the sign convention that flow to a junction is positive and flow from a junction is negative, this term will be preceded by a plus sign when j is the junction for which the equation is written. The derivative with respect to Hi will be positive. The derivative with respect to Hj will be negative. If the term describes a pipe whose flow leaves the junction, a negative sign will precede the term and i will be the junction for which the equation is written, and the derivative for the other node j will be positive. Thus the off-diagonal elements of the Jacobian matrix are positive, and the diagonal elements are negative, as we have already seen in an example. We will utilize this symmetry property in developing an algorithm to generate the Jacobian. However, we first examine alternative means for evaluating the derivatives. The direct way to differentiate {(Hi − H j ) /K } 1/n , in which the pipe number subscript k has been omitted, is to use the power rule of calculus to obtain ±[{(Hi − H j ) / K } (1−n)/n ] / (nK) (4.65) in which the minus sign applies when differentiating with respect to Hj. When a pump is present in the pipe, however, it is no longer a straightforward process to differentiate this term, as it now is {(Hi + hp - Hj)/K}1/n, in which hp = hp(Q) is normally expressed as hp = AQ2 + BQ + C. Another way to obtain these derivatives is to start with Q = Hi − H j K 1/n (4.66) and compute the differential of this formula as dQ = {Hi − H j } 1/n−1dH / (nK1/n ) = Q1−ndH / (nK) (4.67) The partial derivative with respect to Hi is then ∂Q /∂Hi = Q 1−n/(nK) (4.68) and the partial derivative with respect to Hj is identical, except for a minus sign. So Jacobian matrix elements can be obtained quickly via Eq. 4.68. With this approach we can compute the Jacobian terms for a pipe with a pump in it. Also write {(Hi - Hj)/K}1/n as |(Hi - Hj)/K|1/n-1(Hi - Hj)/K to allow a sign change for flows that oppose the assumed direction, which may occur during an intermediate iteration even if the assumed direction is correct for the final solution. When a pump is present in a pipe, then we can write Hi − H j + AQ 2 + BQ + C − KQn = 0 (4.69) Following the procedure of computing the differential of this equation, we find ∂Q / ∂H = ±1 / (nKQn−1 − 2AQ − B) (4.70) If the derivative is with respect to Hi, choose the plus sign; otherwise choose the minus sign for Hj. Thus, for a pipe containing a pump, Eq. 4.69 is first solved for Q, and this Q is then used in Eq. 4.70 to evaluate the derivatives for the Jacobian matrix. Now we can modify SOLQEQS to solve a system of H-equations. Now please obtain a listing from the CD of SOLHEQS and refer to it as you read this section. Here a NAMELIST (actually an extension of standard Fortran 77, but implemented in many Fortran compilers) will be used to handle options. The NAMELIST will also be used in programs in later chapters. The options that may be in the &OPTIONS list are the following: IU (to set ES or SI units), IQ (to set the discharge units), ID (to give the units for diameters and roughnesses), IC ( = 0 to omit checking the dual network con- nectivity description), VIS (kinematic viscosity), PF (peaking factor), GAMMA (specific weight) and ERR (Newton error criterion). With the H-equations there are no loop energy equations, so the input for loops is eliminated, as is the program segment that generates the loop equations. The section that creates the system of equations will include the junction continuity equations, but this section is modified to implement the new method of obtaining the system Jacobian and the H-equations. In SOLHEQS the length of array H, which stores the nodal heads, has been augmented to include the reservoir heads, so that Hi and Hj now provide the nodal heads at the ends of each pipe, including those that supply the network from a reservoir. So we can easily detect a pump in a pipe, its upstream node number is changed to a negative value. The function subprogram COMPK_N now supplies n , but (1 - 1/n ) is stored in array N for later use. In this program we must have access to discharge values during any iteration for any pipe containing a pump. We do this by computing the discharge from the heads that exist during any iteration by letting ARG = [(K/(Hi - Hj)] and DD = ARG(1 - 1/n); then we find that Q = ARG/DD. Statements following label 146 in the program listing carry out this step. When a pump exists in the pipe, then the Newton method is used to solve Eq. 4.69 by statements found in the DO 145 loop. Example Problem 4.13 Prepare suitable input data to analyze the network of Example Problem 4.11 by using pro- gram SOLHEQS. Only minor modifications to the input data in Example Problem 4.11 are needed. First, because the options are entered via a NAMELIST in program SOLHEQS, the first line of the input data now should contain only four values: the numbers of pipes, nodes, reservoirs, and pumps. The second line of input data begins with &OPTIONS, and the next entries contain the namelist variables that differ from the default values, each followed by an equals sign and the value of that variable. This list is ended with a /. Since no loop data are needed, the two lines of loop data are deleted from the input data for the solution to Example Problem 4.11. Since IQ = 0 is the default value, it need not be included in the list after &OPTIONS. With these changes, the input data are now as follows: 5 3 2 1 1 1.5 0 3 1 - 2 - 3 &OPTIONS IQ=0,VIS=1.417E-5/ 2 1.2 0 2 3 - 4 1 0 1 4000 12 0.002 3 1.0 0 3 2 4 - 5 2 1 3 6000 8 0.002 1 100 3 1 2 4000 8 0.002 5 90 4 2 3 3000 6 0.002 1 4.5 54 4.0 50 3.5 44 5 3 0 2000 6 0.002 * * * Example Problem 4.14 The network below is supplied by the source pump in pipe 1, and a booster pump is needed to get the water over the hill below nodes 2 and 3. A turbine is placed in pipe 6 to extract the extra head after the water is moved over the hill crest. Analyze this network using program SOLHEQS. Diameters are in mm, and lengths in m. The kinematic viscosity is ν = 1.31x10- 6 m2/s. P [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) [5] [6] 0.08 m3/s All e = 0.4 mm T P 80 m 400 - 500 75 m 150 - 2500 350 - 1500 0.01 m3/s 150 - 2000 5 m 0.04 m3/s 0.03 m3/s 0.05 m3/s 0.025 m3/s 200 - 1500 20 m 200-50020 m 300-1000350 - 1500 100 m Diameters in mm Lengths in m 90 m Source Pump Booster Pump Turbine Q m3/s hp m Q m3/s hp m Q m3/s ht m 0.20 50 0.20 15 0.15 - 30 0.30 47 0.25 14 0.25 - 25 0.50 43 0.30 12 0.35 - 18 The turbine can be modeled as a pump; the heads are recorded as negative values in preparing its operating characteristics. Since this network is described in SI units, the options for units, discharges and diameters must all be specified. The input data file for this problem, listed in two columns, is therefore 8 6 2 3 2 100 0.04 2 2 - 3 &OPTIONS IU=1,IQ=3,ID=4/ 3 90 0.03 2 3 - 6 1 0 1 500 400 0.4 4 10 0.05 3 4 5 - 8 2 1 2 1500 350 0.4 5 20 0.025 3 6 - 5 - 7 3 2 3 1500 350 0.4 6 5 0.01 1 8 4 1 4 2500 150 0.4 1 80 5 5 4 1500 200 0.4 7 20 6 3 5 1000 300 0.4 1 0.2 50 0.3 47 0.5 43 7 5 0 500 200 0.4 2 0.20 15 0.25 14 0.30 12 8 4 6 2000 150 0.4 6 0.15 - 30 0.25 - 25 0.35 - 18 1 75 0.08 3 1 - 2 - 4 The solution tables from SOLHEQS are PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 2 Q VEL. HEAD LOSS HLOSS/ 1 0 0 0 m mm mm m3/s m/s m 1 0 1 500 400 40.0 0.330 2.63 8.78 17.55 2 - 1 2 1500 350 40.0 0.217 2.26 23.02 15.35 3 2 3 1500 350 40.0 0.177 1.84 15.39 10.26 4 1 4 2500 150 40.0 0.033 1.87 76.55 30.62 5 5 4 1500 200 40.0 0.027 0.86 6.93 4.62 6 - 3 5 1000 300 40.0 0.147 2.08 15.87 15.87 7 5 8 500 200 40.0 0.095 3.03 27.83 55.66 8 4 6 2000 150 40.0 0.010 0.57 5.89 2.95 Pump 1 in pipe 1: Head = 46.22 m, Q = 0.330 m3/s Pump 2 in pipe 2: Head = 14.77 m, Q = 0.217 m3/s Pump 3 in pipe 6: Head = - 30.11 m, Q = 0.147 m3/s NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. Estimate m3/s m M kPa m 1 0.1 0.080 75.0 42.45 416.2 117.45 2 0.0 0.040 100.0 9.19 90.1 109.19 3 0.0 0.030 90.0 3.80 37.3 93.80 4 0.1 0.050 10.0 30.90 303.0 40.90 5 0.0 0.025 20.0 27.83 272.9 47.83 6 0.0 0.010 5.0 30.01 294.3 35.01 * * * 4.4.7. SYSTEMATIC SOLUTION OF THE ∆Q-EQUATIONS In this section we describe the development of the computer program SOLDQEQS that is based on the ∆Q-equations and analyzes pipe networks. This program requires the user to specify the initial discharges, Qoi, so they satisfy all of the junction continuity equations, because algorithms that do this automatically involve considerable logic. We will also omit the input that provides the dual description of the network connectivity; instead we will generate the pipe numbers that interconnect the network nodes from the data on the nodes at the ends of the pipes. This generated data will be used to verify that the input Qoi do satisfy the junction continuity equations. Finally, this program will not allow a PRV or any similar device in the network. With this restriction the Jacobian matrix will be symmetric and positive definite, thereby allowing a special linear algebra solver that requires only the upper (or lower) triangular and diagonal elements of the Jacobian to be available during the solution process. This approach provides us a solution variant that could also be used in solving the H-equations by the Newton method. To describe the computer program logic that forms the ∆Q-equations and the derivatives that form the Jacobian elements, it will be convenient to be able to refer to the equations and the nonzero derivatives with respect to ∆Q from an example. At this time obtain a listing of SOLDQEQS from the CD so it can be studied while you read the rest of this section. The network in Example Problem 4.14 will serve the purpose of illustrating the logic of this program. Since this network contains several pumps, one of which produces a negative head as a turbine, this example will help us incorporate pumps correctly into the code. The two ∆Q-equations for this network are F1 = K1(Qo1 + ∆Q1) n1 − hp1 + K4 (Qo4 + ∆Q1 − ∆Q2 ) n4 − K5 (Qo5 − ∆Q1 + ∆Q2 ) n5 + K7 (Qo7 + ∆Q1) n7 − 80 + 20 = 0 (4.71) and F2 = K2 (Qo2 + ∆Q2 ) n2 − hp2 + K3(Qo3 + ∆Q2 ) n3 + K6 (Qo6 + ∆Q2 ) n6 − hp3 + K5 (Qo5 − ∆Q1 + ∆Q2 ) n5 − K4 (Qo4 + ∆Q1 − ∆Q2 ) n4 = 0 (4.72) In these equations the head hpj of pump j is described by hpj = Aj (Qoi ± ∆Qk ∑ )2 +Bj (Qoi ± ∆Qk ∑ ) + Cj , and the coefficients A, B, and C are chosen to fit three pairs of points along the pump curve, as before. These energy equations are written around the same loops in which the corrective discharges ∆Q1 and ∆Q2 circulate. Therefore, every term in Eq. 4.71 contains a ∆Q1, and every term in Eq. 4.72 contains a ∆Q2. The Jacobian [D] will have two rows, one for each of the two equations, and two columns corresponding to the two unknowns ∆Q1 and ∆Q2, or D [ ] = ∂F1 ∂∆Q1 ∂F1 ∂∆Q2 ∂F2 ∂∆Q1 ∂F2 ∂∆Q2 (4.73) in which the individual elements are ∂F1 ∂∆Q1 = n4K4 (Qo4 + ∆Q1 − ∆Q2 ) n4 −1 + n5K5 (Qo5 − ∆Q1 + ∆Q2 ) n5−1 + n7K7 (Qo7 + ∆Q1) n7 −1 + n1K1(Qo1 + ∆Q1) n1−1 − 2A1(Qo1 + ∆Q1) − B1 (4.74) ∂F1 ∂∆Q2 = ∂F2 ∂∆Q1 = − n4K4 (Qo4 + ∆Q1 − ∆Q2 ) n4 −1− n5K5 (Qo5 − ∆Q1 + ∆Q2 ) n5−1 (4.75) ∂F2 ∂∆Q2 = n2K2 (Qo2 + ∆Q2 ) n2 −1 + n3K3(Qo3 + ∆Q2 ) n3−1 + n6K6 (Qo6 + ∆Q2 ) n6 −1 + n5K5 (Qo5 − ∆Q1 + ∆Q2 ) n5−1 + n4K4 (Qo4 + ∆Q1 − ∆Q2 ) n4 −1 − 2A2 (Qo2 + ∆Q2 ) − B2 − 2A3(Qo6 + ∆Q2 ) − B3 (4.76) To allow for the possibility that one of more flows might change direction and Qoi ± ∆Qk ∑ would become negative, the quantities Ki (Qoi ± ∆Qk ∑ )ni will be rewrit- ten as Ki |(Qoi ± ∆Qk ∑ )|ni −1(Qoi ± ∆Qk ∑ ) . Doing this will be convenient since all factors but the last are also needed to evaluate terms in the derivatives. For this program we must define the loops around which (1) the energy equations will be written, and (2) each ∆Q circulates. Thus the user must supply the pipe numbers which define each energy loop, with a negative pipe number whenever the direction around the loop opposes the assumed direction of flow in that pipe. This information was also required as input to SOLQEQS. These loop data determine the terms in each equation and the sign of each term. As in SOLQEQS, this data resides in a one-dimensional integer array IK, with a pointer LP to separate individual loops. The corrective loop discharge data for each pipe is in a similar array IK1, with a pointer LP1 to separate entries be- tween individual pipes. Thus the positions in array IK1 that will contain information on a corrective loop discharge through pipe I will start with subscript (argument of the array) LP1(I) + 1 and end with subscript LP1(I+1). Thus LP1 must have dimension NP + 1. In a similar way LP must have dimension NL + 1 = NP - NJ + 1. Let us now examine an algorithm to obtain the corrective loop discharges in each pipe from the loop information. The pipes around the two loops in the example network are Loop 1: 1 4 - 5 7 Loop 2: 2 3 6 5 - 4 and this data will be stored in IK as follows: IK(1) = 1, IK(2) = 4, IK(3) = - 5, IK(4) = 7, IK(5) = 2, IK(6) = 3, IK(7) = 6, IK(8) = 5 IK(9) = - 4, with LP(1) = 0, LP(2) = 4, and LP(3) = 9. Since ∆Q1 circulates through loop 1 and ∆Q2 circulates through loop 2, we see that the loop number (the argument of LP) gives the corrective loop discharge through a pipe when the pipe number occurs in the list of IK's for that loop. For example, since pipe 4 is a pipe number in loops 1 and 2, the corrective loop discharges ∆Q1 and ∆Q2 both circulate through it, and also ∆Q1 is in the same direction as the assumed flow in pipe 4 since it is positive in loop 1, whereas ∆Q2 opposes the assumed flow since it is nega- tive in the list of pipes in loop 2. The number of corrective loop discharges through a pipe is not known in advance, so it is simpler to use a two-dimensional array initially, with the pipe number as the first subscript and the number of corrective loop discharges through that pipe as the second subscript. Hence the second subscript of this array must equal or exceed the maximum number of ∆Q's passing through any pipe, so most of this array space will be unused; once these numbers are known, the information can be trans- ferred into the one-dimensional array IK1. Then the two-dimensional array can be deallo- cated and the memory used for other purposes. An alternative for this array is to EQUIVALENCE it to another array that is not used until later, such as the array for the Jacobian matrix. Figure 4.29 lists Fortran statements that could be used to generate these arrays, with the array LP1 zeroed before beginning this algorithm. A very similar algorithm can be used to generate the pipe numbers that join at each node. The essential difference is that the upstream and downstream nodes in the arrays L1 and L2 identify the node to which the pipes attach. In program SOLDQEQS this started in the DO 24 loop. Since we want to verify that the user-supplied initial discharges Qoi DO 50 I=1,NL DO 50 J=LP(I)+1,LP(I+1) II=IABS(IK(J)) NI=LP1(II)+1 IK2(II,NI)=ISIGN(I,IK(J)) 50 LP1(II)=NI NI=0 NCT=NI DO 54 I=1,NP DO 53 J=1,LP1(I) NI=NI+1 53 IK1(NI)=IK2(I,J) LP1(I)=NCT 54 NCT=NI LP1(NP+1)=NI Figure 4.29 Listing of Fortran code to generate arrays IK1 and LP1. do satisfy all of the junction continuity equations, this check is performed immediately after the pipes that join at each node are identified. This information makes it easy to obtain the heads H at the nodes after the discharges and head losses in the pipes are computed by using essentially the algorithm that is in SOLQEQS for this purpose. Now let's see how to obtain the system of ∆Q-equations and the Jacobian that are needed to implement the Newton method. The symmetry that occurs in the Jacobian, if devices such as PRV's do not exist, will be advantageously used, and a one-dimensional array will store the banded portion of the Jacobian. In SOLDQEQS these tasks are accomplished within the outer DO 90 loop. The index I in this loop tracks the NL loop equations, and the equation values are generated and stored in the array F. The process begins with F(I)=F(I)+FI* ... The columns of the Jacobian matrix are each related to a ∆Q, and these values are placed in the one-dimensional array IK1. The pipe numbers in each loop, which identify the terms that are needed to evaluate the equations and the Jacobian elements, are stored in the one-dimensional array IK. The array LP is a pointer that separates consecutive equations, e.g. loops, in array IK. In a banded matrix all elements which are displaced more than the band width from the diagonal are zero. If i is the row number and j is the column number, then the band width NBAND is the maximum difference between a nonzero element in any column and its row number, plus one, or NBAND = |j - i|max + 1 (4.77) In some literature this definition is the half band width since, if the matrix is not symmet- ric, as many elements must to be stored to the left of the diagonal as to its right. In any symmetric matrix [A] each element Aij = Aji. If a two-dimensional array is used in a computer program to store the elements of a banded matrix, the first subscript (for rows) must be at least as large as the number of equations to be solved, and the second subscript must be as least as large as 2NBAND - 1. A special algorithm that properly accounts for the matrix properties is needed to solve a banded matrix problem. If the banded matrix is symmetric, it is not necessary to store all of the elements above and below the diagonal if the solution algorithm accounts for this symmetry. Either the elements above and on the diagonal, or those below and on the diagonal, are all that must be stored. Program SOLDQEQS uses a one-dimensional array to store the banded elements of the Jacobian and calls a linear algebra subroutine SYMMAT to return the solution to the linear equation system in the array F. Before calling SYMMAT, the upper triangular portion of a banded symmetric Jacobian matrix is stored in a one-dimensional array DJ. In the declaration portion of SOLDQEQS we will find that DJ is a one-dimensional allocatable array with DJ[ALLOCATABLE](:) and that the number of real positions to store values in DJ is allocated with ALLOCATE(DJ(NL*NBAND-MM)), in which NBAND is the band width and MM = NBAND - 1. Thus a preliminary task is to determine the band width before allocating DJ and storing the nonzero derivative values in it. The listing in Fig. 4.30 determines the required band width. C FINDS BAND WIDTH MM=0 DO 56 I=1,NL DO 56 J=LP(I)+1,LP(I+1) IJ=IABS(IK(J)) DO 56 JJ=LP1(IJ)+1,LP1(IJ+1) II=IABS(IK1(JJ))-I IF(II.GT.MM) MM=II 56 CONTINUE NBAND=MM+1 Figure 4.30 Band width algorithm. The first position in array DJ is the diagonal element in the first row. The diagonal element of the second row is in position (2 - 1)NBAND + 1, the position of the diagonal element in the third row is (3 - 1)NBAND + 1, and in general the diagonal element in the ith row is in position id = (i - 1)NBAND + 1. An alternative formula for locating the diagonal position is id = i NBAND - MM in which MM = NBAND - 1, the number of elements beyond the diagonal. Thus we see that the storage that is needed by DJ is NL*NBAND - MM (NL is the number of equations), as given in the ALLOCATE state- ment. The position of off-diagonal elements in this one-dimensional array will be the diagonal position id plus the difference between the column number and the row number for the element. In any equation this position is iu = id + (j - i) = id - i + j. Thus in SOLDQEQS the statement after DO 90 I=1,NL that is used to define the NL equations is ID=NBAND*I-MM, which locates the position of the diagonal element for each row, and the statement NI = ID - I is an integer which locates the nonzero off-diagonal positions in DJ when the column position is added. Thus the statements that store the values in the proper locations of DJ are DJ(NI+JJ1)=DJ(NI+JJ1)+FI*FLOAT(IK1(JJ)/JJ1)*DD1 87 CONTINUE 90 DJ(ID)=DJ(ID)+DD1 The portion of the code, within the DO 90 loop in program SOLDQEQS, that generates the system of equations and the values for the Jacobian and then obtains the solution that is used as the Newton correction, consists of the lines listed in Fig. 4.31: DO 90 I=1,NL IB=NBAND*I-MM NI=IB-I II=LP(I)+1 II1=LP(I+1) DO 90 J=II,II1 IJ=IABS(IK(J)) IF(I.GE.NRES.OR.J.GT.II) GO TO 83 IJ1=IABS(IK(II1)) DO 80 JJ=1,NRES IF(IRES(JJ).EQ.IJ) F(I)=F(I)-ELE(JJ) IF(IRES(JJ).EQ.IJ1) F(I)=F(I)+ELE(JJ) 80 CONTINUE 83 FI=IK(J)/IJ QQ=Q(IJ) DO 85 JJ=LP1(IJ)+1,LP1(IJ+1) JJ1=IABS(IK1(JJ)) 85 QQ=QQ+FLOAT(IK1(JJ)/JJ1)*DQ(JJ1) DD=K(IJ)*ABS(QQ)**N(IJ) DD1=DD*(N(IJ)+1.) IF(L1(IJ).LT.0) THEN JJ=1 DO 86 WHILE (IPUMP(JJ).NE.IJ) 86 JJ=JJ+1 DD1=DD1-2.*AP(JJ)*QQ-BP(JJ) F(I)=F(I)+FI*(DD*QQ-(AP(JJ)*QQ+BP(JJ))*QQ-CP(JJ)) ELSE F(I)=F(I)+FI*DD*QQ ENDIF DO 87 JJ=LP1(IJ)+1,LP1(IJ+1) JJ1=IABS(IK1(JJ)) IF(JJ1.LE.I) GO TO 87 DJ(NI+JJ1)=DJ(NI+JJ1)+FI*FLOAT(IK1(JJ)/JJ1)*DD1 87 CONTINUE 90 DJ(IB)=DJ(IB)+DD1 C SOLVES LINEAR EQUATIONS CALL SYMMAT(NL,NBAND,DJ,F) Figure 4.31 The solution algorithm. To enhance solution efficiency we might try to arrange the equations to reduce the band width as much as possible. Not only will a smaller band width reduce the required amount of computer memory for a solution, but it also reduces the computational effort in solving the linear equation system. As the loop data are developed, the user can reduce the band width of the Jacobian matrix by trying to arrange the ∆Q numbering to be as close as possible to the equation numbering . The band width will equal the maximum difference in any equation between the equation number and the ∆Q number, plus 1. However, placing this burden on the user is not desirable, especially since a banding algorithm can readily be implemented in computer code that will probably achieve a tighter banding than the user could arrange, even after some attention is given to the order in which equations should be listed and loops formed. One approach to minimizing the band width is described by Jeppson and Davis (1976). This approach is implemented in SOLDQBAN.FOR, which is on the CD. Also on the CD is SOLDQEQ1 that does not use the band width of the Jacobian but instead uses the standard linear algebra solver SOLVEQ, as do SOLQEQS and SOLHEQS, as it solves the ∆Q-equations. Example Problem 4.15 In the sketch is a network with 10 pipes and 6 nodes which contains three pumps and one turbine. Prepare input data files for SOLQEQS, SOLHEQS and SOLDQEQS so these programs can be used to analyze this network. Use the pairs of (Q, hp) data in the table to define the pump curves. Then replace the pump curve for pump 1 with the new pump data listed later in the solution, and resolve the problem with all three programs. P1 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] [6] 0.04 m3/s All e = 0.0001 m P2 P3 T1 0.10 m3/s 0.07 m3/s 0.05 m3/s 0.04 m3/s 0.06 m3/s 0.05 m3/s 0.25 - 2000 0.25 - 2000 0.20 - 2000 245 m 0.45 - 1000 220 m 228 m 200 m 180 m 170 m 160 m 160 m 200 m 0.2-9000.2 - 600 0.2 - 80 0 (10) 0.2-8000.25-900v = 1.31 x 10-6 m2/s 170 m 0.35-800Diameters in mm Lengths in m Pump 1 Pump 2 Pump 3 Turbine Q m3/s H m Q m3/s H m Q m3/s H m Q m3/s H m 0.40 20.0 0.12 16.0 0.06 8.0 0.09 - 8.0 0.42 18.0 0.15 15.0 0.08 7.5 0.10 - 7.5 0.44 15.0 0.18 13.6 0.10 6.8 0.11 - 6.8 Since SI units are used, options must be changed from the default values. The input data file for each of these three programs are listed next, using two columns for each set: Input Data For SOLQEQS 10 7 2 4 2 'U' 1 'D' 2 4 0.06 180 2 6 -7 1 0 1 1000 0.45 0.0001 5 0.04 170 2 7 -8 2 1 2 800 0.35 0.0001 6 0.07 160 4 4 5 8 -9 3 2 3 2000 0.25 0.0001 7 0.04 160 2 9 -10 4 1 6 2000 0.25 0.0001 1 245 5 3 6 800 0.20 0.0001 10 200 6 1 4 900 0.25 0.0001 1 0.40 20 0.42 18 0.44 15 7 4 5 2000 0.20 0.0001 2 0.12 16 0.15 15 0.18 13.6 8 5 6 900 0.20 0.0001 4 0.06 8 0.08 7.5 0.1 6.8 9 6 7 600 0.20 0.0001 5 0.09 -8 0.10 -7.5 0.11 -6.8 10 7 0 800 0.20 0.0001 4 1 4 9 10 1 0.05 200 4 1 -2 -4 -6 4 2 3 5 -4 2 0.05 228 2 2 -3 4 4 -8 -7 -6 3 0.10 220 2 3 -5 Input Data For SOLHEQS 10 7 2 4 2 0.05 228 2 2 -3 &OPTIONS IU=1,IQ=3,ID=2/ 3 0.10 220 2 3 -5 1 0 1 1000 0.45 0.0001 4 0.06 180 2 6 -7 2 1 2 800 0.35 0.0001 5 0.04 170 2 7 -8 3 2 3 2000 0.25 0.0001 6 0.07 160 4 4 5 8 -9 4 1 6 2000 0.25 0.0001 7 0.04 160 2 9 -10 5 3 6 800 0.20 0.0001 1 245 6 1 4 900 0.25 0.0001 10 200 7 4 5 2000 0.20 0.0001 1 0.40 20 0.42 18 0.44 15 8 5 6 900 0.20 0.0001 2 0.12 16 0.15 15 0.18 13.6 9 6 7 600 0.20 0.0001 4 0.06 8 0.08 7.5 0.10 6.8 10 7 0 800 0.20 0.0001 5 0.09 -8 0.10 -7.5 0.11 -6.8 1 0.05 200 4 1 -2 -4 -6 Input Data For SOLDQEQS 10 7 2 4 3 0.10 220 &OPTIONS IU=1,IQ=3,ID=2/ 4 0.06 180 1 0 1 1000 0.45 0.0001 0.44 5 0.04 170 2 1 2 0800 0.35 0.0001 0.20 6 0.07 160 3 2 3 2000 0.25 0.0001 0.15 7 0.04 160 4 1 6 2000 0.25 0.0001 0.12 1 245 5 3 6 0800 0.20 0.0001 0.05 10 200 6 1 4 0900 0.25 0.0001 0.07 1 0.4 20 0.42 18 0.44 15 7 4 5 2000 0.20 0.0001 0.01 2 0.12 16 0.15 15 0.18 13.6 8 5 6 0900 0.20 0.0001 -0.03 4 0.06 8.0 0.08 7.5 0.10 6.8 9 6 7 0600 0.20 0.0001 0.07 5 0.09 -8.0 0.10 -7.5 0.11 -6.8 10 7 0 800 0.20 0.0001 0.03 4 1 4 9 10 1 0.05 200 4 2 3 5 -4 2 0.05 228 4 4 -8 -7 -6 The solution tables from SOLQEQS and SOLDQEQS are identical, as shown below. SOLHEQS failed to converge. The failure was caused by the relative inaccuracy of the initial values of the heads that were provided to the Newton method by the automated estimator in the code; the values were too crude in relation to the sensitivity of the code to the way that the three pairs of points for pump 1 define its operating characteristics. If PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 4 Q VEL. HEAD LOSS HLOSS/ 1 0 0 0 m m m m3/s m/s m 1 0 1 1000 0.450 1.0 0.436 2.74 12.61 12.61 2 1 2 800 0.350 1.0 0.163 1.70 5.39 6.73 3 2 3 2000 0.250 1.0 0.113 2.31 36.76 18.38 4 1 6 2000 0.250 1.0 0.118 2.40 39.64 19.82 5 3 6 800 0.200 1.0 0.013 0.42 0.74 0.92 6 1 4 900 0.250 1.0 0.105 2.14 14.29 15.88 7 4 5 2000 0.200 1.0 0.045 1.43 19.20 9.60 8 5 6 900 0.200 1.0 0.005 0.16 0.14 0.15 9 6 7 600 0.200 1.0 0.066 2.10 11.92 19.87 10 7 0 800 0.200 1.0 0.026 0.82 2.56 3.20 Pump 1 in pipe 1: Head = 15.71 m, Q = 0.436 m3/s Pump 2 in pipe 2: Head = 14.44 m, Q = 0.163 m3/s Pump 3 in pipe 4: Head = 6.02 m, Q = 0.118 m3/s Pump 4 in pipe 5: Head = - 5.17 m, Q = 0.013 m3/s NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. Estimate m3/s m m kPa M 1 0.1 0.050 200.0 48.10 471.7 248.10 2 0.1 0.050 228.0 29.15 285.8 257.15 3 0.1 0.100 220.0 0.39 3.8 220.39 4 0.1 0.060 180.0 53.81 527.7 233.81 5 0.0 0.040 170.0 44.61 437.5 214.61 6 0.1 0.070 160.0 48.46 475.2 208.46 7 0.0 0.040 160.0 42.56 417.3 202.56 this pump curve is plotted, we see immediately how the curve turns steeply downward outside each end of the given data. Using pump curves of this nature should be avoided. To obtain a solution from SOLHEQS, either the points that define the pump curve must be adjusted, or the code must be modified so the user can supply the initial estimates of the heads for the Newton method. If the pump curve for pump 1 is modified so the three discharge-head data pairs are (0.40, 16.0), (0.43, 15.8), and (0.46, 15.5), then SOLHEQS can solve the modified problem. * * * 4.5 CONCLUDING REMARKS This chapter concentrated on the analysis of pipeline networks. The first area of empha- sis was on the development of the three kinds of systems of equations to describe mathe- matically the flow in a pipe network, first for simpler networks and then for networks which contain pumps or turbines and loss-producing devices such as a pressure reducing valve or a back pressure valve. The Newton method for the solution of these equation systems was introduced and later included in computer programs to solve the equation systems. Later sections of the chapter developed solution routines and implemented them for each of the three types of equation systems. There are features that a production network program would usually include that these programs do not have. For example, rather than requiring the user to provide a set of esti- mated initial discharges Qoi that satisfy all of the junction continuity equations, the pro- gram should develop these values. One way to create these values is to reduce the network to a branched system by deleting some pipes with smaller diameters and then using the methods in this chapter to obtain a solution for the branched network. Another burden that would not be placed on the user of a production program is the need to supply the pipe numbers that define the loops around which the energy equations are written and the corrective loop discharges circulate. An algorithm for this task must satisfy two criteria: (a) the pipes that define any loop should be minimum in number, and also (b) these loops must lead to the creation of energy equations that are independent so that none of the equations are a combination of any group of the other equations. The first criterion can be achieved by using a "minimum path algorithm," and the second criterion requires each new loop to contain at least one pipe that does not exist in any of the previous loops. Production programs will also take full advantage of the sparsity of the Jacobian matrix in computing network solutions in an efficient manner. Network solvers can also allow the user to obtain time-dependent solutions. Such solu- tions, which do not account for the forces that are required to accelerate or decelerate the fluid columns, have become known as "extended time simulations." To develop an extended time simulation, additional information of several kinds is needed, such as demand functions which describe how one or several external nodal demands QJ vary with time, rules based on pressures at nodes or on water surface elevations in tanks or reservoirs can determine how many pumps should operate in series or parallel, and storage-elevation- capacity curves can be used to describe the behavior of tanks, and so on. The use of programs for network analysis can also allow designers to obtain answers for the many questions that naturally occur during the design process. For example, what head and capacity should a pump produce to maintain a prescribed pressure and/or discharge at the far end of the network? What is the discharge from a junction if the pressure is known from a measurement there? How much head should a PRV dissipate so the pressure does not exceed a set value? How much flow can be obtained from a fire hydrant if its discharge characteristics are known? What are the flows from sprinkler heads if their sizes are known? How does a contaminant spread through a pipe network if it is accidentally introduced at a point? Chapter 5 will explore the design of these pipe networks, and Chapter 6 will examine further several topics, including extended-time simulations. 4.6 PROBLEMS 4.1 For the two pipe networks shown below, write the system of Q-equations. In writ- ing these equations, use K and n with subscripts that correspond to the pipe number. [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) 0.5 m3/s 0.18 m3/s 0.15 m3/s 0.3 - 500 P [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) [5] [6] All e's for both networks = 0.00002 m 0.25 m3/s 0.1 m3/s 0.1 m3/s 0.5 - 500 0.5 m3/s 0.25 m3/s 0.2 - 500 0.4 - 500 HGL = 300 m 0.2-6000.5-6000.4-600100 m 0.45 - 1600 0.5 - 1500 0.6 - 900 0.25-14000.5 - 1800 0.4-90080 m 0.4-500Diameters in mm Lengths in m Pump Q m3/s h p m 0.2 30 0.4 27 0.7 21 4.2 Write the system of Q-equations for the network shown. It is not necessary to substitute the values of K and n from the table into the equations; instead use Ki and ni where i is the pipe number. If the discharge in pipe 1 is Q1 = 3.1 ft3/s, then what is the friction factor f for this pipe? Pipe K n 1 1.841 1.928 2 11.47 1.871 3 7.47 1.839 4 1.615 1.914 5 11.08 1.828 6 7.69 1.884 [4] [2] (5) (6) (2) (3) [1] (1) [3] (4) 1.1 ft3/s 120' 8" - 3000' 1.4 ft3/s 100' 8" - 3000' 8" - 2000'8" - 2000'10" - 1300' 10" - 1500' 1.0 ft3/s 1.3 ft3/s 4.3 A 5-pipe, 3-node network appears below. On this diagram the first number along each line is the pipe diameter in inches, and the second number is the pipe length in feet. All pipes have an equivalent sand roughness e = 0.001 ft = 0.012 inches. Do the follow- ing: (a) compute the values of K and n in hf = KQn for pipe 1, based on the Darcy- Weisbach equation, and (b) write the system of Q-equations for this network. (Use sub- scripts on K, n and Q corresponding to the pipe number.) [2] (5) (2) (3) [1] (1) [3] (4) 1.1 ft3/s 140' 1.2 ft3/s 100' 6" - 2500' 0.9 ft3/s All e = 0.001' = 0.012" v = 1.217 x 105 ft2/s 8" - 600' 8" - 1000' 8" - 800 ' 6" - 400 ' Pipe K n 1 2 3.53 1.961 3 4.44 1.929 4 48.6 1.934 5 6.40 1.817 4.4 In the sketch the network consists of 6 pipes and 3 nodes. A source pump and one reservoir supply the network, and the lower reservoir receives water. Do the following tasks: (a) write the system of Q-equations; (b) write the system of H-equations; (c) write the system of ∆Q-equations; (d) if the discharge in pipe 5 is Q5 = 0.026 m3/s into the reservoir, what is the elevation of the HGL at node 3; (e) if the discharge in pipe 6 is Q6 = 0.112 m3/s, what are the HGL and pressure at node 2? [2] (5) (6) (2) (3) 0.02 m3/s P [1] (1) [3] (4) 0.05 m3/s 300 m 280 m 170 m 180 m 190 m 280 m 0.03 m3/s Pipe Dia. m Length m K n 1 0.30 1000 543 1.886 2 0.20 2500 13700 1.946 3 0.20 1000 3270 1.839 P t . Q hp 4 0.30 1500 1077 1.965 m3/s m 5 0.15 1000 27400 1.974 1 0.05 35 6 0.35 800 260 1.968 2 0.10 31 3 0.15 24 4.5 For the network below: (a) write the Q-equations; (b) write the H-equations; and (c) write the ∆Q-equations. (Use the symbols K and n with correct subscripts for the pipes in the equations.) Pipe K n 1 3.59 1.922 2 7.97 1.917 3 7.94 1.821 4 28.80 1.809 [2] (2) (3) [1] (1) (4) 1.0 ft3/s 70' 100' 6" - 2000' 1.5 ft3/s 10" - 3000' 6" - 500 ' 8" - 2200' 4.6 For the network shown, write (a) the Q-equations, (b) the H-equations, and (c) the ∆Q-equations. The K and n values for the pipes in this network are given in the table which follows. (Your equations should contain only numerical values and unknowns.) Pipe D in. L ft. K n 1 8 1500 5.72 1.930 2 6 2000 33.00 1.931 3 6 1000 16.30 1.889 4 8 1700 6.53 1.913 5 6 2500 40.70 1.890 (2) (3) [1] (1) (4) 0.5 ft3/s 90' 120' 6" - 2500' 0.6 ft3/s 140' [2] (5) All e = 0.005" 8" - 1700' 6" - 200 0' 8" - 1500' 6" - 100 0' 4.7 Prepare the input data for, and obtain the solution from, NETWK for the network described in Problem 4.6. 4.8 A pipe branches into a 6-in diameter, 1500-ft long pipe and a 8-in diameter, 1400-ft long pipe and then rejoins, so the two pipes are in parallel. Pipe 1 contains an open globe valve with a local loss coefficient K = 10. If the total discharge is Q = 3 ft3/s, de- termine the discharges Q1 and Q2 in the individual pipes. For simplicity, we shall as- sume f1 = 0.018 and f2 = 0.015. (2) (1) Q = 3 ft3/s e = 0.005" 6" - 1500' f1 = 0.018 8" - 1400' Open globe valve K = 10 f2 = 0.015 4.9 This sketch of a small water system shows two reservoirs, with a pipe connected to the center node with an inflow of 1.0 ft3/s at the other end. Set up the three equation systems that could be used to solve this problem, and then obtain a solution by using one of them. (2) (3) [1] (1) 100' 8" - 1800' 1.0 ft3/s 80' [2] e = 0.012" 8" - 2000' 6"-1000' 4.10 In the diagram three pipes that form a triangle are supplied by a reservoir at one ver- tex of the triangle, and demands of 1.0 and 3.0 ft3/s are found at the other two vertices. A booster pump exists in pipe 3. What head should the pump in pipe 3 produce so the pressure at node 2 causes an HGL of 60 ft there. The ground elevation is everywhere 0 ft. You may solve this problem by using any, or all, of the equation systems that are available. Assume e = 0.012 in and ν = 1.41x10- 5 ft2/s. [2] (2) (3) 3.0 ft3/s P [1] (1) 1.0 ft3/s 100' Head = 60' 8" - 6000' 6" - 8000' Pump 10"-2000'hp = ? 4.11 In the network of Problem 4.10 the diameters of pipes 2 and 3 have been enlarged to 12 inches. In place of the pump a pressure reduction valve is now needed in pipe 3 to create a pressure head at node 2 of 60 ft. Determine the required setting, i.e. HGL, of the pressure reduction valve. You can use any of the equation systems. 4.12 A network is shown in the diagram. Write the system of ∆Q-equations for this network and complete one Newton iteration toward a solution of the problem. Verify this result by using an application software package such as MathCAD or TK-Solver. [4] [2] (5) (2) (3) 0.5 ft3/s 0.5 ft3/s P [1] (1) [3] (4) e = 0.005' 2.0 ft3/s 1.0 ft3/s 8" - 3000' 4" - 3000' 8"-1500'4"-1500'6" - 3500 ' 40' 30' 50' 160' Pump Characteristics Pipe K n Q h p 1 375.0 1.860 ft3/s ft 2 57.4 1.902 0.3 40 3 190.0 1.898 0.5 35 4 11.5 1.880 0.8 28 5 5.72 1.930 4.13 Write the system of H-equations for the two networks in Problem 4.1. 4.14 Write the system of ∆Q-equations for the two networks in Problem 4.1. 4.15 The following network contains a pressure reducing valve (PRV) that is set so it will produce a HGL of 145 m on its downstream side. This valve is 800 m downstream from node 1. Do the following: (a) write the system of Q-equations; (b) write the system of H-equations; (c) write the system of ∆Q-equations; (d) using the Newton method, solve the system of ∆Q-equations; (e) and what is the HGL on the upstream side of the PRV? [1] (1) [3] [2] (5) (4) (2) (3) 0.05 m3/s All e = 0.15 mm 0.08 m3/s 0.06 m3/s 200 m 140 m v = 1.31 x 10-6 m2/s HGL = 145 m 400 - 1000 250 - 2 500 250 - 2000 250 - 3500 300 - 500 800 m All elev. = 100 m Pipe 1 2 3 4 5 K 196 3520 2380 4130 192 n 1.819 1.955 1.895 1.892 1.834 4.16 The reservoir water surface elevation at the beginning of pipe 1 in Problem 4.15 is lowered by 50 m so it is WS1 = 150 m, and a pump with the characteristics given below is installed in pipe 1. Write the three equation systems and solve one of them, also finding the HGL elevation on the upstream side of the PRV. Q h p m3/s m 0.18 55 0.22 51 0.26 44 4.17 For the network shown below write (a) the Q-equations, (b) the H-equations, and (c) the ∆Q-equations. Pipe 3 contains a pressure reduction value 200 ft down- stream from node 1 that is set to maintain an HGL = 430 ft on its discharge side. Use the notation Ki and ni in these equations. Pipe K n 1 1.93 1.935 2 4.44 1.940 3 3.50 1.840 4 47.90 1.866 5 7.67 1.917 [1] (1) [3] [2] (5) (4) (2) (3) All e = 0.005' 520' 420' v = 1.27 x 10-5 ft2/s HGL = 430' 0.9 ft3/s 1.2 ft3/s 1.0 ft3/s 12" - 4000' 8" - 2000' 10" - 3500' 10"-3000'350' 280' 320' PRV 200' 6"-3000'4.18 For the network below: (a) write the Q-equations; (b) write the H-equations; (c) write the ∆Q-equations; and (d) solve the system of ∆Q-equations. The water surface elevation of the right reservoir is 300 ft, and the following three (Q, hp) pairs define the pump characteristic curve: (1.0, 26), (1.5, 24), (2.2, 20). [2] (5) (6) (2) (3) P [1] (1) [3] (4) 2.0 ft3/s 80' 2.0 ft3/s 4" - 1000' 4" - 3000' 6" - 5000' 6" - 6000' 8" - 30 00' All elev. = 0' All e = 0.005' 10" - 1000' 4.19 Solution tables from NETWK follow, with four values omitted. Fill in the missing values. What head drop occurs across the PRV? What horsepower does this loss represent? PIPE DATA PIPE NO. N O D E S FROM TO L DIA. e x10 3 Q VEL . HEAD LOSS HLOS S / 1 0 0 0 ft. in in ft3/s ft/s ft. 1 0 1 4000 12.0 5.0 4.72 6.01 9.69 2 1 2 3500 10.0 5.0 3.26 5.97 41.87 11.96 3 1 3 3000 10.0 5.0 0.46 0.85 0.91 0.30 4 2 3 3000 6.0 5.0 0.44 2.22 10.21 3.40 5 2 0 2000 8.0 5.0 1.62 4.64 19.36 9.68 AVE. VEL. = 3.94 ft/s, AVE. HL/1000 = 7.01, MAX. VEL. = 6.01 ft/s, MIN. VEL. = 0.85 ft/s NODE DATA NODE D E M A N D ELEV. HEAD PRESSURE HGL ELEV. ft3/s gal/min ft. ft. lb/in2 ft. 1 1.0 449 350 131.2 56.87 481.2 2 1.2 539 320 119.4 51.72 439.4 3 0.9 404 280 AVE. HEAD = 133.3 ft, AVE. HGL = 450.0 ft MAX. HEAD = 149.2 ft, MIN. HEAD = 119.4 ft 4.20 For the network shown: (a) write the Q-equations; (b) write the H-equations; (c) write the ∆Q-equations; and (d) solve the ∆Q-equation system. P [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) [5] [6] All e = 0.0005 m 0.025 m3/s 0.005 m3/s 0.02 m3/s 0.30 - 1000 0.015 m3/s 0.35 - 600 0.045 m3/s 0.20-7500.25-7500.35-50060 m 40 m 38 m 35 m 35 m 36 m 33 m 0.25 - 1000 0.025 m3/s 0.20-1500Diameters in m Lengths in m Q h p m3/s m 0.120 40 0.140 38 0.165 35 4.21 For the two networks in Problem 4.1, solve the Q-equation system using the New- ton method. 4.22 For the two networks in Problem 4.1, solve the H-equation system using the New- ton method. 4.23 For the two networks in Problem 4.1, solve the ∆Q-equation system using the Newton method. 4.24 Determine the pressures in lb/in2 at the six nodes of Problem 4.1a. 4.25 For the network below, write the ∆Q-equation system and solve them, and verify your solution by obtaining a computer solution by using NETWK. 90' 100' [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] [6] 1.0 ft3/s 0.6 ft3/s 0.5 ft3/s 6" - 1800' 8" - 2500' HGL = 80' P All e = 0.005" 0.5 ft3/s 10" - 1000' 8" - 1500' All elev. = 20' 1000' 6" - 2600' 6" - 2800' 0.6 ft3/s 6" - 1000' 6"-2000'1.5 ft3/s 8"-2000' Q h p ft3/s ft 1.0 50 2.0 48 3.0 45 4.26 Determine the discharge and head loss in each pipe of the networks shown on the following pages by first determining a set of values for the coefficients K and n; then setting up and solving the equations without using a computer, except perhaps to solve each linear system of equations that is formed as a part of the Newton method. (a) Analyze this network with the Hazen-Williams equation and CHW = 120 for all pipes. [1] (1) [4] [3] [2] (5) (4) (2) (3) 500 gal/min 6"-1000'12" - 3000' 12" - 1500' 8"-1000'10" - 3 50 0' 1500 gal/min 2000 gal/min I II (b) Analyze this network by using the Hazen-Williams equation; for pipes 1 through 5 use CHW = 120, and for pipes 6 through 11 use CHW = 100. [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] [6] 6.5 ft3/s 2.0 ft3/s 4.5 ft3/s 12" - 3000' 12" - 3000' 1.0 ft3/s 6" - 2000' 10"-2000'[7] (11) (10) 10"-2000'12" - 4250' 10" -2500' II I III IV 6" - 3000' 14"-5000'10"-1000'10"-1000'V 2.0 ft3/s 1.0 ft3/s (c) Use the Hazen-Williams equation to analyze this network; all pipes are cast iron. 0.15 m3/s 0.1 m3/s 0.1 m3/s 0.2 m3/s 0.1 m3/s 0.05 m3/s 0.2 m3/s [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] [6] 1500 m - 0.3 m [7] (11) (10) 600m-0.3mII I III IV 500 m - 0.2 m 700 m - 0.2 m 700 m - 0.4 m 1000m-0.5m1000m-0.5m900 m - 0.2 m 600m-0.3m700 m - 0.2 m 500m-0.2mV (d) Analyze this network with the Darcy-Weisbach equation, assuming e = 0.001 ft for the 8-in and 10-in pipes and e = 0.0005 ft for all other pipes. [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] [6] [7] (11) (10) II I III IV (13) (14) (12) [9] [8] 4 ft3/s 5 ft3/s 4 ft3/s 15" - 3000' 5 ft3/s 15" - 3000' 15"-2500'8"-1000'V 10 ft3/s 10 ft3/s 15" - 800' VI 10"-1500'10" - 1500' 12" - 1500' 12" - 25 00' 10"-1000'10"-3350'21"-2000'10 " - 21 20 ' 2 ft3/s 15" - 800' (e) Analyze the network of part (a) by using the Darcy-Weisbach equation with a rough- ness of e = 0.0005 ft for all pipes. (f) Analyze the network of part (b) by using the Darcy-Weisbach equation; for pipes 1 through 5 use e = 0.005 ft, and for pipes 6 through 11 use e = 0.006 ft. (g) Use the Hazen-Williams equation to analyze this network; all pipes are made of cast iron. The diameters are given in centimeters; lengths and elevations are in meters. Pump performance data are listed in the table. Pump 1 Pump 2 Q h p Q h p m3/s m m3/s m 0.10 50 0.05 10 0.15 48 0.10 8 0.20 44 0.15 5 20 - 1000 P2 [1] (1) [3] [2] (5) (4) (2) (3) 100 m 0.03 m3/s 0.09 m3/s 0.06 m3/s P1 15 - 10 00 40 m 30 m 60 m 0 m 15-150015 - 1500 20" - 1500' (h) Use the Hazen-Williams equation to analyze this network that is diagrammed atop the next page; all pipes are made of cast iron. The diameters are given in centimeters; lengths and elevations are in meters. Pump performance data are listed in the table. 120 m 150 m [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] (10) [6] 0.057 m3/s 120 15 - 300 15 - 600 20-300120 106 0.01425 m3/s 110 15 - 400 6 - 800 (PRV)2 HGL = 45' 15-60015-600HGL = 45 m (PRV)1 200 120 (PRV)3 HGL = 45 m 21015-55020-3000.01425 m3/s 0.0285 m3/s 115 P1 P2 15-760120 Diameters in cm Lengths in m Pump 1 Pump 2 Q h p Q h p m3/s m m3/s m 0.0425 36.6 0.0283 12.2 0.0708 30.5 0.0425 10.7 0.0991 22.9 0.0566 8.5 4.27 Only the Colebrook-White equation is used in subroutine COMPK_N that deter- mines the values of K and n in the exponential formula for programs SOLQEQS, SOLQEQS and SOLDQEQS. Modify this subroutine so it will allow laminar flow in the pipe. Also modify the subroutine so the discharge can become zero; e.g., this might commonly occur for initial discharges Qoi for use with the ∆Q-equations. 4.28 Assume that the flow in all pipes will always be turbulent; however, a user might select initial values for Qoi that are zero when solving the ∆Q-equations. Then the sub- routine COMPK_N would fail, as it is now written in SOLQEQS, SOLHEQS and SOLDQEQS. Modify COMPK_N so it can accept a value of zero for the discharge. 4.29 Modify SOLQEQS so an option allows the user to supply starting values for Q that will be used in the Newton method rather than generating these values internally. 4.30 Modify SOLHEQS so it has an option that allows the user to supply initial values for H that will be used in the Newton method rather than generating these values internal- ly. The additional input could be supplied from another read statement, or these heads could be listed after the nodal elevations in the node data. 4.31 Modify SOLDQEQS so it will allow PRV's to exist in the network. This change will require two sets of loop data to be read as input data (unless you wish to obtain these loops internally), one around which the ∆Q's circulate, and one around which the energy equations are written. Since the two sets of loops will not be identical, the Jacobian will not be symmetric. 4.32 Modify SOLHEQS and/or SOLDQEQS so it calls a symmetric matrix solver such as SYMMAT. 4.33 SOLQEQS, SOLHEQS and SOLDQEQS can all analyze networks that contain local losses if the user will provide the actual length of the pipe and the additional length of pipe that would cause a frictional pipe loss that is equivalent to what the local loss device would cause. Modify one, or all, of these programs so each equivalent pipe length is computed internally within the program and then added to the actual length before the problem is solved. 4.34 Rather than compute an equivalent length of pipe for a local loss, as in Problem 4.33, modify SOLQEQS so that local losses, where they occur, are treated by adding a head loss term of the form hL = KQ2/(2gA2) to the energy loop equation. 4.35 Repeat Problem 4.34, but modify SOLDQEQS. 4.36 Use SOLQEQS, SOLHEQS, and SOLDQEQS to analyze the network depicted in Problem 4.20. 4.37 Use SOLQEQS to analyze the network in Problem 4.25. This network contains a PRV in pipe 3 that is located 1000 ft downstream from the beginning of this pipe. 4.38 SOLQEQS, SOLHEQS and SOLDQEQS all represent pump performance by fit- ting three (Q, hp) pairs of pump characteristic curve data with a second-order polynomial. Modify one or all of these programs so they accept the normal capacity (discharge at best efficiency) and head at this discharge as input, and then the relation between hp and Q is obtained from the power equation P = γQhp under the assumption that the power P re- mains constant. 4.39 Modify the program that was developed in Problem 4.38 so the efficiency of the pump is a linear function of the difference of the discharge from its normal capacity. 4.40 Place a PRV in pipe 2 of the network in Example Problem 4.5 with a pressure setting of HGL = 445 ft. Obtain a solution for this network using SOLQEQS. Verify this solution using NETWK. 4.41 SOLQEQS contains a code segment that cross checks the connectivity of the net- work by looking at the two node numbers at the ends of a pipe and at the pipe numbers that join at a junction. It also checks that upstream node numbers are negative and that downstream node numbers are positive. But the algorithm currently can not determine whether an extra pipe might be connected to a node. Modify the code so a check can iden- tify any extra pipe(s) that might be specified in the data that lists the pipes that are con- nected to nodes. 4.42 Modify SOLDQEQS so PRV's can exist in the network. Now two separate kinds of loops will exist, those around which the corrective loop discharges circulate and those around which the energy equations are written. Therefore two sets of loop data must be in- cluded in the input data file. 4.43 Use the resulting computer program from Problem 4.42 to obtain a solution to the network in Example Problem 4.5 with a PRV in pipe 2 having a pressure setting that causes the downstream head to be HGL = 445 ft. Verify this solution by (1) using NETWK and by (2) changing subroutine FUNCT in program EQUSOL1. 4.44 The network diagram below lists average demands on it. The storage tank that is connected to the network by pipes 14 and 16 has a 20-m diameter; at 6 a.m. its water surface elevation should be 200 m. The demands at all nodes change according to the peaking factors in the table. The pump characteristics represent two pumps in parallel at each location. Obtain a series of solutions for the times at which the peaking factors are given. For each solution of this series determine the new water level in the tank and the electrical energy consumed by each pump during the latest time increment. Suggest when one pump at each station should be shut off. What might be done to improve the design and thereby the operation of the system? Time 6 a.m. 9 a.m. 12 Noon 4 p.m. 7 p.m. 10 p.m. 12 Mid. 3 a.m. PF 1.0 1.8 1.3 1.3 1.7 1.5 0.6 0.2 Pump 1 Pump 2 Q h p Q h p m3/s m m3/s m 0.15 50 0.20 30 0.25 47 0.25 28 0.35 42 0.30 25 [1] (1) [4] [3] [2] (5) (6) (4) (2) (3) (7) (8) (9) [5] [7] (11) (10) (13) (14) (12) [9] [8] (16) (17) (18) [10] (15) [11] [12] P1 P2 0.25 - 1500220 m 0.4 - 3500 0.3 - 2000 0.3 - 2000' 0.30 - 3000 0.25 - 1800 0.3 - 1800 0.20 - 2200 0.20 - 2200 210 m 200 m 210 m 200 m 200 m 135 m 130 m 185 m 180 m 175 m 200 m 165 m 170 m 175 m [6] 0.3 - 15000.3 - 15000.25 - 15000.25 - 15000.25 - 15000.20 - 15000.20 - 15000.04 m3/s All e = 0.0002 m 0.03 m3/s 0.06 m3/s 0.025 m3/s 0.05 m3/s 0.03 m3/s 0.02 m3/s 0.03 m3/s 0.03 m3/s 0.025 m3/s 0.03 m3/s 0.03 m3/s 0.35-5000.25-10000.25-1000Diameters in mm Lengths in m </p>