UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analysis of fluid transients in large distribution networks Karney, Bryan William 1984

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1984_A1 K36.pdf [ 6.67MB ]
Metadata
JSON: 831-1.0062640.json
JSON-LD: 831-1.0062640-ld.json
RDF/XML (Pretty): 831-1.0062640-rdf.xml
RDF/JSON: 831-1.0062640-rdf.json
Turtle: 831-1.0062640-turtle.txt
N-Triples: 831-1.0062640-rdf-ntriples.txt
Original Record: 831-1.0062640-source.json
Full Text
831-1.0062640-fulltext.txt
Citation
831-1.0062640.ris

Full Text

A N A L Y S I S O F F L U I D T R A N S I E N T S I N L A R G E D I S T R I B U T I O N N E T W O R K S by B R Y A N W I L L I A M K . A R N E Y B.A.Sc.(Hons.) , Un ive rs i t y o f B r i t i sh Co l umb i a , 1980 M.Eng . , Un ive rs i t y o f B r i t i sh Co l umb i a , 1982 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S Depar tment o f C i v i l Eng ineer ing W e accept this thesis as conforming to the requ i red standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 1984 ® B ryan W i l l i a m Ka rney , 1984 In present ing this thesis i n part ia l fu l f i lment o f the requ i rements for an advanced degree at the The Un ive r s i t y o f Br i t i sh C o l umb i a , I agree that the L ib ra ry shal l make it freely ava i lab le for reference and study. I further agree that permiss ion for extensive copy ing o f this thesis for scholar ly purposes may be granted by the Head o f my Depar tment or by his or her representatives. It is understood that copy ing or pub l i ca t ion o f this thesis for financial gain shal l not be a l l owed wi thout my wr i t ten permiss ion. Department of C i v i l Engineering The Un ive r s i t y o f Br i t i sh C o l u m b i a 2075 Wesb rook Place Vancouver , Canada V 6 T 1W5 Date: September 1984 Abstract It is we l l k nown that the pressures generated dur ing transient condit ions shou ld be an impor tant cons iderat ion when p ipe l ine systems are be ing designed or constructed. I f the size and strength o f the requ i red p ipe is to be rat ional ly selected, i f the surge suppress ion equ ipment is to be log ica l ly s ized and i f system operat ing rules are to be inte l l igent ly speci f ied, re l iable transient analysis is essential. T h e a im o f this research has been to deve lop a general, eff ic ient and re l iable a lgor i thm for comput ing the response o f large d is t r ibut ion systems to rap id transient condit ions. A l t hough further re f inement and ver i f i cat ion work is st i l l required, it is cons idered that this goal has been ach ieved. S igni f icant ly , the a lgor i thm is general , since a wider class o f networks can be ana lyzed than was prev ious ly possible. The p rogram is also eff ic ient, since the so lut ion o f many o f the network boundary condi t ions is made expl ic i t for the first t ime and because computer storage is conserved i n the p rogram implementat ion . F i na l l y , the p rogram is re l iable, since a number o f numer i ca l exper iments have p roven its correctness. Spec i f ica l ly , the so lut ion procedure for obta in ing the transient response o f large networks inc ludes the fo l l ow ing new features: a coup l ing o f the steady and unsteady parts o f the analysis avoids dup l icat ion and improves p rogram eff ic iency; a general network and boundary cond i t ion c lassi f icat ion system wh i ch permi ts networks wi th any reasonable topology to be ana lyzed; so lut ion procedures for many common network boundary condi t ions are made expl ic i t or pu t i n a f o rm wh i ch a l lows eff ic ient numer i ca l ca lcu lat ion; and, f ina l l y an automat ic a lgor i thm for select ing the t ime step and d i v id ing the network p ipes into an integer number o f reaches. A l l o f these elements combine to produce eff ic ient and accurate transient analysis o f large d is t r ibut ion systems. In summary , the deve loped a lgor i thm a l lows comprehens ive transient analysis o f networks hav ing arbi trary geometry to be pe r fo rmed . T h e networks may inc lude any i i general combinat ions o f hydrau l i c devices and no restr ict ion is p laced on either the numbe r o f boundary condit ions or p ipes that connect to a g iven node or the number o f sections f ound i n any p ipe i n the network. Procedures have been ut i l i zed wh i ch keep data requirements to a m i n i m u m despite the general nature o f the p rob lem being so lved. F i na l l y , the program has been carefu l ly tested and compared wi th known solut ions for s imple p ipe l ines w i th very good agreement be ing observed. i i i Tahle of Contents Abstract » U s t o f Tab les v i i L i s t o f F igu res v i i i 1. I N T R O D U C T I O N 1 2. G O V E R N I N G E Q U A T I O N S 8 2.1 Introduct ion 8 2.2 Bas ic Equat ions 8 2.2.1 Assumpt ions 11 2.2.2 Wavespeed r 12 2.3 D i scuss ion o f So lut ion Me thods 17 2.4 M e t h o d o f Character ist ics 20 2.5 Conc lus ions 23 3. S T E A D Y - S T A T E N E T W O R K A N A L Y S I S 24 3.1 Introduct ion 24 3.2 Statement o f the P rob l em 26 3.3 Bas ic Equat ions : 28 3.4 So lu t ion Techn iques 31 3.4.1 The P ipe F l o w Equat ions 32 3.4.2 The N o d e Equat ions 33 3.4.3 The L oop Equat ions 34 3.5 Compa r i s on o f the So lu t ion Techn iques 36 3.6 Character ist ics o f the Adop ted So lut ion Techn ique 38 4. B O U N D A R Y C O N D I T I O N S 39 4.1 Int roduct ion 39 4.2 Boundary Cond i t i on C lass i f i cat ion System 41 4.3 Character i zat ion o f the Ne twork 46 iv 4:4 S imp le and O rd i na r y O n e - N o d e Boundary Cond i t i ons 50 . 4.4.1 S imp le Nodes and L i nea r Reservo i rs 52 4.4.2 Externa l F l ows w i th Energy D i ss ipa t ion 56 4.4.3 A i r Chambe r s 59 4.4.4 Power Fa i l u re o f a P u m p w i th Inert ia 62 4.4.5 O the r Boundary Cond i t i ons 65 4.5 O rd ina ry T w o - N o d e Boundary Cond i t i ons 65 4.5.1 Booster Pumps Unde rgo ing Powe r Fa i l u re 67 4.5.2 T w o - N o d e Boundary Cond i t i on w i th Energy D iss ipa t ion 68 4.5.3 Pressure Reduc i ng Va lves 69 4.5.4 P ipe Rep lacement E lements 71 4.6 Boundary Cond i t i on Extens ions 76 4.6.1 Gene ra l Junct ion Equa t i on w i th Losses ..V. 78 4.7 Conc lus ions 82 5. D E V E L O P M E N T O F T H E N E T W O R K M O D E L 85 5.1 Introduct ion 85 5.2 Ne two rk Representat ion 86 5.2.1 Internal Sect ions 91 5.3 Boundary Cond i t i ons Representat ion 94 5.3.1 Examp l e Ne two r k w i th L i n k e d Boundary Cond i t i ons 96 5.4 Ne two rk D iscret i za t ion: Select ing the T i m e Step 104 5.4.1 The D iscret i za t ion A l go r i t hm 108 5.5 Conc lus ions 112 6. V E R I F I C A T I O N O F T H E S O L U T I O N P R O C E D U R E 114 6.1 Ver i f i ca t ion C r i t e r i a 115 6.2 Boundary Cond i t i o n Fo rmu la t i on 116 6.2.1 Examp le : Va l v e I n - L i n e Bounda ry Cond i t i on 116 v 6.3 Ne tworks w i th Spec ia l Symmet ry 118 6.3.1 Insc r i p t i on o f the Test Ne two rk 119 6.3.2 E xamp l e 1: Sudden Va l ve C losure 122 6.3.3 E xamp l e 2: L i nea r Va l ve C losure w i th D iamete r Changes 126 6.3.4 E xamp l e 3: C omp l e x Va l v e Movemen t 128 6.4 Convergence to Steady-State 129 6.5 Conc lus ions 131 7. Conc lus ions 132 B I B L I O G R A P H Y 134 vi Tab le N o . Desc r ip t i on L is t o f Tab les Page 1 Da t a Representat ion o f Ne two rk i n F igu re 12 90 2 Da t a Representat ion for Character is t ic Sections o f Ne two rk i n F i gu re 13.93 3 Da t a Requ i r ed to Locate Boundary Cond i t i ons i n F i gu re 14 97 4 Da t a Structure for Ne two rk i n F igu re 14 98 v i i List of Figures F i gu re N o . Desc r ip t ion Page 1 S imp le P ipe l ine System 9 2 x - t P lane wi th Character ist ics G r i d 21 3 Ne two rk w i th T r e e - Shaped Structure 26 4 Looped Ne twork Topo logy 27 5 S imp le N o d e 42 6 C o m p l e x N o d e 43 7 L i n k e d Ne two rk Segment 47 8 O rd i na r y N o d e 50 9 Reservo i r w i th Head Losses and Inert ia 54 10 A i r C hambe r 59 11 Gene ra l T w o - N o d e Boundary Cond i t i o n 66 12 S imp le Ne twork Examp l e 90 13 S imp le Ne two rk Examp l e 93 14 Examp l e Ne twork w i th L i n k ed Boundary Cond i t i ons 96 15 x - t P lane wi th Character ist ics G r i d and Interpolat ion 105 16 Va l ve I n - L i n e 117 17 Sketch o f Test Ne two r k 120 18 N o d e o f Degree 4 a long P r ima ry A x i s o f Symmetry 121 19 Trans ient Response: Sudden Va l ve C losure 124 20 T i m e - H e a d Prof i les: Sudden Va l ve C losure 125 21 Test Ne twork wi th Banded Structure 127 22 Trans ient Response: L inea r Va l ve C losu re 128 23 Trans ient Response: Comp l e x Va l v e M o v e m e n t , 130 v i i i 1. INTRODUCTION The p rob l em o f ana lyz ing f lu id transients i n large d is t r ibut ion systems is one o f the most complex o f a l l unsteady f low prob lems. A large d ist r ibut ion system, such as a water supply network for a city, may inc lude thousands o f branched or looped pipes and a host o f comp lex operat ing devices. F o r example, the network may contain a numbe r o f pump ing stations, smal l and large storage reservoirs, various surge suppress ion devices as we l l as numerous k inds o f valves. U n d e r extreme operat ing condi t ions, local two -phase f low or f low wi th a free surface may occur. In a network, a l l o f these elements are comb ined together in to a topo log i ca l^ complex f ramework governed by non l inear part ia l d i f ferent ia l equations. N o t on ly is the arrangement o f the var ious hydrau l i c devices i n a network compl icated, but the f low condi t ions the network may be requ i red to meet are also ext remely var ied. Typ i ca l condi t ions may range f r om the smal l f l ow demands wh i ch occur i n a water d is t r ibut ion system in the winter, when energy may have to be diss ipated to prevent network pressures f r om becoming excessive, to the large discharges wh i c h occur i n the summer, when spr ink l ing demand can increase f l u i d velocit ies to the i r max imum leve l . In add i t ion to the da i ly and seasonal changes i n demand, most networks must be designed to supply emergency flow condit ions. The large, rap id ly occur r ing discharge requ i red to fight a major fire is an example. Th i s part icular s i tuat ion is further compl icated by uncertainty as to the size and locat ion o f the requ i red f low. Hyd r au l i c condi t ions i n water d i t r ibut ion systems are i n an almost cont inual state o f change. Industr ia l and domest ic users are constantly changing their f low requirements. Supp l y condi t ions undergo adjustment as the water levels i n reservoirs and storage tanks change and as var ious pumps are turned o f f and on . However , not a l l o f these changes requ i re that a fu l l transient analysis be per fo rmed. I f the changes 1 2 in flow condi t ions take place very s lowly, i t is possible to go gradual ly f r om one steady-state cond i t ion to another, w i thout any large f luctuat ions i n pressure head or f low velocity. In add i t ion, i f pressures are generated by surge waves ( in wh ich the effects o f f l u i d compress ib i l i ty are smal l ) , eff ic ient analysis is possible using mass analysis techniques. However , i f rap id changes occur, whether caused by standard operat ing procedures or by accidental events, large f luctuat ions in pressure head may be produced. The magni tude o f these pressure head changes depends on the details o f the system's conf igurat ion and the specif ic condit ions wh ich in i t iated the transients. I f the sizes and strengths o f the var ious p ipes in the network are to be rat ional ly selected, the magni tude o f these rap id pressure f luctuations must be accounted for in the p ipe l ine select ion process (Sharp, 1981). The analysis o f these rap id pressure changes is the subject o f the current research. F o r s imple p ipel ines, conta in ing on ly series p ipe connect ions or a few branch ing pipes, the analysis tools requi red to complete detai led transient calculat ions are now wel l deve loped. Howeve r , for very large and complex d is t r ibut ion systems, fu l l transient analysis has not yet been possible. A l t h ough attempts at this analysis p rob l em have been made, the procedures and approaches have on ly been successful on relat ively smal l networks. The lack o f a power fu l analysis tool for large networks has led to a number o f d i f f icul t ies and shortcomings i n des igning d istr ibut ion systems. A few typical examples are l isted be low. 1. Sharp and Cou l s on (1968) relate an inc ident i nvo l v ing a large thermal power station pump i ng out o f a l ow pressure supply ma in . W h e n this system was operated, "unexpected water hammer ' ' occured wh ich caused a major water supply fa i lure. 2. Koe l l e , in h is open ing address to the International Institute on Hydraulic Transients and Cavitation (Sao Paulo, 1982), relates how transient pressures 3 caused a number o f concrete p ipes to fa i l after on ly seven years o f operat ion. T h e p ipes were i n a u rban por t ion o f a B raz i l i an city's water supply network and fa i led as a direct result o f transient condi t ions i n the network. 3. G rea te r Vancouver Water D i s t r i c t has made a "substantial ef fort" to predict the surge pressures generated by sudden power fa i lure i n Vancouver ' s water supply system. Th i s attempt has been largely unsuccessful and the lack o f an adequate analys is tool has been a h indrance to the rat ional deve lopment and design o f this system. 1 As i de f r om speci f ic operat ional prob lems, such as those related above, there may be a large economic cost associated w i th over-conservat ive network design. Th i s aspect is re lated by Jeppson (1976): In the absence o f such analysis to determine the per formance o f the exist ing system under the new demands accompanied by studies for p rov i d i ng the most economica l alternatives for e l iminat ing def ic iencies, needlessly large investments are made for larger than necessary pipes, redundant l ines or dupl icate faci l i t ies. Par t ly because such errors are bur ied, the tax payer, or consumers, se ldom become aware o f these mistakes. A l t hough Jeppson was p r imar i l y concerned w i th steady-state analysis, the statement is even more appl icab le to transient design condit ions. F o r example, i t is st i l l current practice to design network p ipes w i th diameters i n the range o f four to sixteen inches w i th a factor o f safety o f "at least" four (see, for example, Ame r i c an Wate r Wo r k s Assoc ia t ion Standards C400 and C401 govern ing the design o f asbestos-cement p ipe) . It is reasoned that the damage potent ia l i n these p ipes is s ignif icant, but that "the analysis o f transient condit ions and the imp lementa t ion o f surge protect ion equ ipment are . bo th unnecessary and impract ica l for p ipes i n this size range" (Hodgson , 1983, p. 36). T h e current research is an attempt to overcome the prob lems associated wi th the earl ier attempts to mode l transient condit ions i n large networks. ' Re l a t ed i n a pr ivate commun ica t ion by A . Pu rdon , C h i e f Eng ineer and Manager o f Opera t ions for the Grea te r Vancouver Reg iona l Distr ic t , B.C., Canada . 4 %The primary objective of this study is to develop an efficient and reliable procedure for analyzing large distribution systems in the transient state by computer. Specifically, the procedure should satisfy five criteria: 1. it should be efficient in execution time, allowing large and complex systems to be economically modeled; 2. it should be efficient in the use of computer storage, allowing large networks to be studied when computer memory is limited; 3. the procedure should be easy to use, permitting input data to be simply prepared, corrected and modified; 4. the algorithm should be reliable, so that transient pressures can be predicted with reasonable confidence; 5. and, finally, the procedure should be sufficiently general to allow networks having a high degree of complexity to be accurately represented. Previous attempts to analyze large networks have all suffered in at least one of the points itemized above. Typically, the program deficiencies have been relatively unimportant when small or medium sized networks of relatively simple topology are being considered. In these applications, program inefficiencies may simply result in the computer executing for a slightly longer time in order to produce an answer or that a somewhat greater effort must be expended to prepare the input data. However, for very large networks of high complexity, the deficiencies are crucial. When these large systems are considered, the cost of the computer analysis is often significant, the memory requirements are large and the effort in preparing (and obtaining) input data can be substantial. In these cases, the usefulness and value of an algorithm or program strongly depends on resources being used in an efficient manner. Probably the key characteristic of the network problem is its generality. The pipes, pumps, valves, reservoirs and other hydraulic devices which form a large distribution system can be connected together into an almost infinite number of 5 combinat ions and arrangements. Desp i te this, a survey o f the avai lable l i terature re lat ing to transient analysis o f networks reveals an a lmost total lack o f any coherent conceptual f ramework for cons ider ing the large network p rob lem. Researchers have persisted i n cons ider ing the network p rob l em as noth ing more than a s imple extension o f procedures app l ied to s imple p ipe l ine systems. Such an approach, relegating as it does, network topology to a mere ly secondary role, ignores those aspects o f the network p r ob l em wh i ch make it both interest ing and chal lenging. In part icular, the techniques for so lv ing boundary condi t ions i n networks have features that are rarely encountered i n s imple p ipe l ine appl icat ions. The most important example occurs when boundary condi t ions interact direct ly w i th each other to fo rm, what is cal led i n this thesis, a linked group. In s imple p ipe l ine appl icat ions, a typica l example o f a l i nked group occurs i n m u l t i - u n i t pump ing stations. In order to solve the boundary condi t ions o f such stations, a l l o f the unknowns i n the p rob l em must be assembled into a set o f s imultaneous equations. Th i s system o f unknowns requires special attention i n the so lut ion procedure. However , i n networks, l i nked groups o f boundary condi t ions can occur much more frequent ly. In order to be able to recognize, ident i fy, solve and, at t imes, avo id l i nked boundary condit ions, a fresh look is needed at the who le quest ion o f network boundary condit ions. Chap te r 4 o f this thesis attempts to p rov ide the concepts and termino logy requ i red for a general approach to this p rob l em. Th i s dissertation is organized i n a topica l manner. Chapters 2 through 4 fo l low the classical o rder for def in ing a p r ob l em invo lv ing a t ime-dependen t part ia l d i f ferent ia l equat ion. Name l y , Chap te r 2 describes the d i f ferent ia l equations themselves, l is t ing the assumptions requ i red to der ive them and the techniques appl icab le to obta in their solut ion. Chap te r 3 then def ines the in i t ia l condi t ions requ i red to actual ly compute a so lut ion. De t e rm in i ng the in i t ia l condi t ions i n a network usual ly requires the steady-state network p rob l em to be solved. Th i s steady-state so lut ion is a d i f f i cu l t p rob l em in its own right; however, an accurate so lut ion is essential i f the transient 6 condi t ions are to be accurately computed. In essence, the mater ia l i n Chapters 2 and 3 funct ions as a l i terature rev iew by cons ider ing recent developments i n the field o f hydrau l i c transients and by discussing current techniques for steady-state network analysis. T he def in i t ion o f the part ia l d i f ferent ia l equat ion p rob l em is completed in Chap te r 4 when the network boundary condi t ions are cons idered. The boundary condi t ions i n a network are where the transient condi t ions are in i t iated and where the character o f the transient can be mod i f i ed by the mechanics o f wave ref lect ions and transmissions. Because o f this role, Chapte r 4 is the heart o f the network p rob l em as we l l as the central contr ibut ion o f this thesis. Chap te r 4 attempts to take a new look at the condit ions, arrangement and interactions possib le for boundary condi t ions i n a general network o f complex topology. The emphasis at this stage is to deve lop a conceptual f ramework wh i ch w i l l a l low the network p rob l em to be cons idered in a general manner. The restrict ions and assumptions on network topology that are no rma l l y made i n the l i terature w i l l be careful ly avo ided. Once the transient p rob lem has been def ined and the quest ion o f how solut ions can be obta ined has been discussed, it is possib le to consider the procedures and data structures requ i red to actual ly imp lement the so lut ion on a computer. Th i s is the subject o f Chap te r 5. Th i s topic o f imp lementat ion has, unt i l now, been largely ignored i n the field o f hydrau l i c transients, but is key i f the solut ion to very large systems is to be obta ined. Three quest ions are cons idered i n this chapter: ' how can the network be eff ic ient ly represented?'; ' how can the boundary condit ions be descr ibed and the solut ions deve loped i n Chap te r 4 appl ied? ' ; and 'how can the t ime step requ i red to pe r f o rm the transient analysis be selected to keep computer costs to a m in imum? ' . O f course, i f the procedures developed to analyze d is t r ibut ion systems are to be used w i th conf idence, they must be ver i f ied. Chap te r 6 thus considers a number o f numer i ca l exper iments wh i ch are designed to test the general va l id i ty o f the procedures 7 deve loped i n Chapters 2 to 5. The test procedure consists o f developing networks hav ing rather unusual symmetry propert ies. These networks can then be made equiva lent to s imple p ipe l ines where the numer ica l procedure for obta in ing the solut ion is we l l deve loped. I f the network so lut ion and the s imple p ipe l ine so lut ion agree, the general va l id i ty o f the network procedures is establ ished. F ina l l y , the thesis concludes w i th a summary o f what the author considers to be the most s igni f icant contr ibut ions presented i n this dissertation. The p rob l em o f ana lyz ing large d is t r ibut ion systems in the transient state is a complex and demand ing one. However , i t is hoped that one further step has been made to understanding and mode l i ng these systems wh i ch are o f such signif icant economic importance. 2. G O V E R N I N G E Q U A T I O N S 2.1 I N T R O D U C T I O N T h e purpose o f this chapter is to br ie f ly out l ine the equations descr ib ing the unsteady f low o f f lu ids i n c losed conduits. Standard references der ive the basic part ial d i f ferent ia l equat ions i n some detai l (e.g., Chaudhry , 1979; Fox , 1977; Watters, 1984; W y l i e and Streeter, 1983), and no purpose is served by reproduc ing these derivat ions here. Howeve r , where discrepancies or d i f f icu l t ies exist i n the l i terature, addi t ional comments are p rov ided i n order to c lar i fy the situat ion. F l u i d transients are generated whenever condi t ions are changed i n a hydrau l ic network, whether due to standard operat ing procedures or due to accidental events. Th i s chapter introduces the basic equations descr ib ing the phys ica l behav ior o f the f l u i d in the p ipe l ine and discusses the parameters that inf luence the system response. The assumptions, and the l imi tat ions these assumptions impose, are br ie f ly discussed. Standard so lut ion techniques appl icab le to the basic equations are next considered, fo l l owed by a b r i e f descr ipt ion o f the adopted so lut ion technique. F u l l descr ipt ion o f the transient network response addi t iona l ly requires deve lopment o f the boundary and in i t ia l condi t ions o f the d i f ferent ia l equations. These topics are discussed in subsequent chapters. 2.2 RASTC. E Q U A T I O N S T h e mathemat ica l descr ipt ion o f the flow o f a compress ib le fluid i n an elastic conduit general ly requires three equat ions to be wr i t ten: cons iderat ion o f the conservat ion o f mass results i n a k inemat i c or cont inuity re lat ion; the equations descr ib ing the conservat ion o f l inear momen t um lead to a dynamic equat ion; and, finally, considerat ion o f the conservat ion o f energy leads to a thermodynamic equat ion. The thermodynamic equat ion requires an equat ion o f state to be wr i t ten between the temperature, pressure 8 9 and density o f the f l u i d i n quest ion. The unsteady analysis o f l i qu id f low i n p ipe l ine systems is a lmost universal ly assumed to occur under i sothermal condit ions. I f this is the case, the equat ion o f state becomes part icu lar ly s imple and is usual ly expressed as a re lat ion between f lu id density and pressure known as the bu lk modu lus o f the f lu id , wh ich is usual ly assumed to be a constant over the range o f operat ing pressures o f interest i n hydrau l i c appl ications. The bu lk modu lus can be comb ined wi th the p ipe l ine elasticity re lat ions to form a wavespeed or celer ity equat ion. Th i s wavespeed equat ion effect ive ly replaces the thermodynamic re lat ion in l i qu i d transient appl icat ions, reduc ing the number o f pr imary equations that must be cons idered to two — the dynamic and cont inui ty equations alone. F igu re 1 shows a s imple p ipe carry ing l i qu id i n the transient state f rom a RESERVOIR H H +<jH GATE V *• v + d V •6 X dx L x = 0 x=L F igu re 1. S imp le P ipe System 10 reservoir to a valve at the downstream end o f the p ipe. In general, the f low area o f the p ipe l ine is a funct ion o f the centerl ine distance f r om the upstream reservoir (x). T h e dynamic equat ion is der ived by isolat ing a smal l e lement o f f l u i d o f length dx. Cons ide ra t i on o f the forces act ing on the mass o f this element, and the accelerat ion it undergoes as a result o f these forces, produces the dynamic equat ion. T h e result is v t + v v x + g H x + f 4 M = o (1) i n wh i ch V is the instantaneous average f lu id velocity i n the p ipe cross-sect ion, H is the hydrau l i c gradel ine elevation above an arb i t rary datum, g is the accelerat ion o f gravity, f is the D a r c y - W e i s b a c h fr ic t ion factor, D is the p ipe diameter and t is the t ime. In this and subsequent sections, part ia l der ivat ives w i l l be ident i f ied as lower case subscripts. F r i c t i on loss formulas other than the D a r c y - W e i s b a c h re lat ion can be accommodated w i th m ino r changes to the f inal term. Cons idera t ion o f mass conservat ion in the f lu id e lement in F igure 1 leads to the cont inu i ty equat ion. Dens i ty changes i n the f lu id and the stretching o f the p ipe wa l l in the c i rcumferent ia l and long i tud ina l direct ions are accounted for. T he cont inu i ty equat ion can be wr i t ten as H , + V H + — V = 0 (2) t x g x w i n wh i ch a is the pressure wave veloc i ty i n the p ipe l ine. M a n y textbooks inc lude an add i t iona l V s i n c term i n Equat ion (2), where a is the slope o f the p ipe. W y l i e (1984) has recently argued that this addi t iona l te rm shou ld not be p resent The convect ive V V ^ and V H ^ terms i n Equat ions (1) and (2) are f requent ly very smal l re lat ive to the other parts o f the equat ions and it is c ommon to e l iminate these terms (Wy l i e , 1981; W y l i e and Streeter, 1983). Re ta in ing the convect ive terms often introduces unwanted numer ica l attenuat ion and d ispers ion into the so lut ion wh i ch can lead to greater inaccuracy than that caused by neglect ing the convect ive terms 11 (Koe l l e ; 1982). Recent ly, more complex schemes have been deve loped to overcome the numer i ca l p rob lems associated wi th reta in ing the convective terms (Go l dbe rg and Wy l i e , 1983), but the improvements are on ly ach ieved at the expense o f greater computer t ime and storage requirements. F o r the purposes o f this thesis the equations w i l l usual ly be represented in their s imp l i f i ed fo rm. Thus the dynamic equat ion becomes V t + * H X + ^ = 0 (3) S imi l i a r l y , the cont inui ty equaiton is wr i t ten H f + — V = 0 (4) t g x w 2.2.1 ASSUMPTIONS The assumptions made in der iv ing the basic equat ions descr ib ing the unsteady f low o f f lu id in p ipe l ines are br ie f ly summar i zed i n order to exp la in the i r l imitat ions. • It is assumed that the l i qu i d f low is one -d imens iona l and homogeneous. Thus pressure head, density and veloc i ty are assumed to vary on ly in the ax ia l d irect ion, w i th values at any section equal to the average value o f the var iable at that section. Boundary layer and f lu id structure effects are on ly inc luded to the extent they inf luence these mean quantit ies. A l t hough changes i n the conduit c ross-sect iona l area are a l lowed, the p ipe walls cannot diverge too rap id ly o r f low separation w i l l occur and the one -d imens i ona l assumpt ion w i l l be v io lated. • Neg lec t o f the convect ive terms i n the govern ing equat ions is equivalent to V say ing that the rat io o f — is much less than one (Wy l i e , 1981). F o r Si almost a l l reasonable f lu id velocit ies and commerc ia l p ipe mater ia ls carry ing pure water, this rat io is less than 0.01. However , i f a ir is dispersed i n the 12 pipe, or is released from solution during the transient, the wavespeed will be greatly decreased and the fluid velocity may become a significant fraction of the wavespeed. • The pipe is assumed full at all times and the liquid in it does not boil. This means that the minimum pressure at any section must always be greater than the vapour pressure.2 • It is assumed that the pipe and fluid deform linearly as a perfect elastic continuum. The change in density of the fluid under isothermal conditions is described by its bulk modulus. • Friction is evaluated according to an empirical power law expression. In addition, it is assumed that the friction loss at any instant in unsteady flow is equal to the loss which would occur for the same flow under steady conditions. This assumption ignores the complex interaction between the developing boundary layers and the strong pressure gradients that exist in unsteady flow. Any additional assumptions will be listed when they are made. The above assumptions have not, in general, proved restrictive in the transient analysis of hydraulic systems and very good agreement has often been obtained between predicted and measured system response. Because of this, it is expected that the analysis of large networks will not be unduly impaired by these restrictions. 2.2.2 WAVESPEED The speed of propagation of a pressure impulse in a conduit is called the celerity or the water-hammer wave velocity of the pipe (designated as 'a' ' S ome account can be taken o f cavitat ion through the interna l boundary condit ions discussed i n Chapte r 4. 13 •above). The va lue ascr ibed to the celer ity is a funct ion o f the compress ib i l i ty o f the f lu id , the elasticity o f the p ipe l ine mater ia l and the detai ls o f how the p ipe is anchored. Gene ra l wavespeed equations must account for density changes in the f l u i d as we l l as changes i n the length and cross-sect iona l area o f the conduit . Nume rous texts and papers (e.g., Chaudh ry , 1979; Fox , 1977; Ha l l iwe l l , 1963; Watters, 1984; W y l i e and Streeter, 1983) develop expressions for the pressure wave velocity i n the p ipe , and these der ivat ions are not reproduced here. However , some comments are p rov ided where c lar i f icat ion o f the l i terature is i n order. In a recent technical note, W y l i e (1984) has argued that the troublesome and technical ly incorrect V s i n a te rm usual ly inc luded i n the cont inui ty equat ion probab ly arises f r om an incorrect fo rmu la t ion o f the celer ity equat ion. W y l i e notes that the bu lk modu lus and the wavespeed cou ld both be def ined i n terms o f the hydrau l i c gradel ine elevat ion H . Thus the bu lk modu lus o f the f lu id ( K ) is de f ined as K = (5) h i n wh i ch p is the f l u i d density and the dot is used to s igni fy the total der ivat ive o f the argument in quest ion. The wavespeed expression then becomes K a 2 = £ r- (6) i + A p g H in wh i ch A is the cross-sect ional area o f the c ondu i t In prev ious derivat ions, the total der ivat ive o f the pressure replaces the p g H expression wh i ch appears i n Equat ions (5) and (6). It is this pressure term wh i ch former ly led to the annoy ing p ipe l ine slope term in the cont inu i ty equat ion. A s written, the def in i t ions result i n the quoted expression for the cont inuity equat ion and invo lve 14 no inconsistency wi th the steady-state condit ions. 3 Before Equat ion (6) can be used as a pract ical method o f calculat ing celerit ies, the support condit ions o f the p ipe l ine must be cons idered. Usua l l y three cases o f long i tud ina l p ipe l ine support are ident i f ied: case a refers to a p ipe l ine anchored at the upstream end on ly; case b refers to a p ipe anchored against long i tud ina l mot ion throughout its length; and case c represents a p ipe l ine f itted wi th expansion jo ints throughout. Cons idera t ion o f the stresses i n the p ipe wal l leads to the expression K a 2 = ^r- (7) 1 + | * i n wh i ch E is Young ' s modu lus for the p ipe l ine mater ia l and ii is a factor depending on the support conditons. In t h i c k -wa l l e d p ipes the stresses are cons idered to vary across the thickness (t) o f the p ipe l ine mater ia l . In this case, the va lue o f $ can be wri t ten as fo l lows: case a, + = 2(1 + u) + ° ( - ^ - X l - f ) case b, * = 2(1 + n) + ^ ( - ^ X l - »*) (8) case c, * = 2(1 + M ) + r (^TT } i n wh i ch M is the Poisson rat io o f the wa l l mater ia l . 3 A n interest ing feature o f Equat ions (1) and (2) is that they inc lude the steady-state network equat ions as a special case. In steady-state conditons, the derivat ives w i th respect to t ime are zero and the veloc i ty is a lmost constant for conduits o f constant cross sect ion. Equa t i on (1) can then be integrated to recover the D a r c y - W e i s b a c h f r ic t ion loss equat ion and Equat ion (2) is sat isf ied i f density changes are very smal l under steady-state pressure changes. The V s i n a term former ly inva l idated these conclus ions and transient programs wou ld not " h o l d " steady-state values p r i o r to the in t roduct ion o f transients. Th i s annoy ing p rob l em disappears w i th the e l iminat ion o f the p ipe l ine slope term. ! 15 T w o l im i t ing values can be der ived f rom these t h i c k -wa l l ed equations. A s t becomes very large the second term in each expression becomes negl ig ib le and ^ = 2(1 +u) regardless o f the restraint condi t ion. Th i s value o f \p is often used to represent an un l ined tunnel . A l so , as t becomes very smal l , the fract ion j^_t approaches one and the first term becomes ins igni f icant relat ive to the second. The expressions that rema in are often refered to as the t h i n -wa l l e d equat ions and can be der ived independent ly by assuming hoop stresses i n the wa l l . 4 The term (1 - ^ ) m the expression for case a restraint quoted above frequent ly appears in the l i terature as ( ^ - u). Th i s discrepancy is caused by neglect ing the velocity o f the p ipe i n calculat ing the change i n velocity o f the f l u i d i n the alternative der ivat ion. Tha t is, when a p ipe is on ly anchored at one end, the condui t w i l l lengthen by a smal l amount dur ing the t ime it takes for a wave to travel its length. A l t hough the change in length is smal l , the t ime o f travel is also smal l wh i ch results i n a finite value o f the velocity o f the p ipe l ine wa l l itself. W h e n this effect is inc luded in the der ivat ion, the term re lat ing to the long i tud ina l strain cancels and the quoted expression results. Th i s d iscrepancy is not too s igni f icant f r om a practical po int o f view, since on ly a smal l change i n wavespeed results and since case a restraint is relat ively uncommon . Howeve r the discrepancy is annoy ing and its persistence i n the l i terature indicates the elusive nature o f the source o f the error. It shou ld be 4 W y l i e and Streeter (1983) recommend using the t h i n -wa l l e d expressions whenever the rat io o f D to t is greater than 25 wh i le Watters (1984) is more conservative and sets the l im i t at 40. Actua l ly , as Hodgson (1983) shows, the rat io o f D to t is an imperfect ind icator o f wh i ch expression is va l id since cons iderat ion must be given to rat io o f K to E i n the denominator o f Equa t i on (6) as we l l . In practice, Hodgson found that on ly for po l yv iny l ch lor ide and polyethy lene p ipe d i d the dif ference between the th in and thick wa l led formulat ions effect the wave speed by more than three pe r cen t In general, the t h i n -wa l l e d expressions were conservative in that they ove r -es t imated the value o f the wavespeed. 16 noted as we l l that for p ipes that change in length as the pressure in them changes (and therefore may have s igni f icant velocit ies), r igorous analysis requires in t roduct ion o f a th i rd govern ing equat ion. Th i s th i rd equat ion is needed to relate the changes i n p ipe length to f l u id pressure via the p ipe mater ia l and f l u i d propert ies. On l y for the s imple case o f instantaneous closure where fr ict ion is ignored can expl ic i t relat ions between the change in f l u id velocity and the result ing head change be writ ten. However , in practice, the change in p ipe length is a lmost a lways neglected when calculat ing head changes. O n e o f the most s ignif icant d i f f icu l t ies i n calculat ing the wavespeed relates to the poss ib le presence o f a ir and /o r other gases i n the p ipe l ine. It has been theoret ical ly shown ( Ko r ob i , et al., 1955) and exper imenta l ly ver i f ied (Pearsal l , 1965) that even very smal l quantit ies o f gases present i n the l i qu id can drast ical ly reduce celerity values. A i r and other gases can be entrapped in a p ipe l ine dur ing filling operations, they can be drawn into a l ine at the entrance or at a ir in let valves and can be released f r om solut ion when there is a temperature increase or a pressure reduct ion. The quantity o f gases invo lved i n this process and its behav ior under transient condit ions is extremely d i f f icu l t to predict. U sua l l y careful design el iminates most o f these effects and reasonable mode ls have been developed for common ly used devices such as air inlet valves. However , many di f f icu l t ies rema in and much research is current ly in progress to achieve a better understanding o f the processes invo lved. F o r the purpose o f this thesis, it w i l l be assumed that no a i r occurs i n the network except at isolated nodes where air pockets may ex i s t 5 5 Ex c l u d i n g considerat ion o f a i r released f r om solut ion is quite a reasonable assumption for water d is t r ibut ion systems. A s F o x (1977) points out, subatmospher ic pressures are usual ly not accepted i n bur ied p ipe systems because o f the risk o f po l lu t ion f rom groundwater drawn in through p ipe l ine joints. L i t t le or no air release cou ld take place w i thout these low pressure per iods. 17 Hodgson (1983) discusses a number o f other factors wh ich inf luence the celer ity calculat ion. These inc lude the effects o f soi l conf inement, p ipe l ine l iners, temperature, rate o f strain and var iab i l i ty i n mater ia l propert ies. F o r composi te materials, considerat ion must be g iven to wh i c h mater ia l actual ly carries the load. Var ia t ions i n condi t ions between upsurge and downsurge cou ld lead, i n rare instances, to di f ferent celerity values depend ing on the sign o f the pressure change. Usua l l y these effects are not s igni f icant and Hodgson 's tables p rov ide useful summar ies o f reasonable celer ity va lues for a l l common ly avai lable p ipe l ine mater ia ls and dimensions. F o r the purpose o f this thesis, it w i l l subsequent ly be assumed that p ipe l ine celerit ies are known , s i ng l e -va lued constants for each p ipe i n the network. 2.3 DISCUSSION OF SOLUTION METHODS It has been said that "one cannot i n general solve part ia l d i f ferent ia l equations, one can only • approx imate t h em" (Gustafson, 1980, p. 44). Th i s is true o f the equat ions descr ibing the transient f low o f f l u i d i n p ipe l ines. The dynamic and cont inu i ty equations are classif ied as quas i - l i nea r hyperbo l i c part ia l d i f ferent ia l equations, wi th nonl inear i t ies introducted in the f r ic t ion term and i n the wavespeed expression i f smal l quantit ies o f a ir o r other gases are p resen t In add i t ion, many o f the c ommon boundary condi t ions present i n hydrau l i c networks exh ib i t non l inear behavior. Because o f these nonl inear i t ies, no analyt ic solut ions o f the unsteady f low equations are avai lable. ' A large var iety o f graphica l and numer i ca l techniques have been deve loped i n an attempt to obta in approx imate solut ions to the transient govern ing equations. Ea r l y solut ions usual ly requ i red dropp ing the convect ive and f r ic t ion terms f rom the govern ing equations and so lv ing the resultant l inear system by graphica l or algebraic means. G raph i ca l techniques reached a h igh degree o f re f inement i n the m idd le o f this century, w i th the works o f Schnyder (1929, 1932), Bergeron (1961) and Pa rmak ian 18 (1954) " be ing representative. 6 However , in t roduct ion o f the digi ta l computer made ava i lab le a wide range o f improved techniques and these have a l lowed progressively-more comp lex systems to be accurately analyzed. The re are essentially three alternative numer ica l methods to consider for the analysis o f hydrau l i c transients i n c losed-condu i t systems. These inc lude the me thod o f characterist ics, f in i te dif ference techniques and f in i te element methods. O f these, by far the most popu lar is the method o f characteristics. Each o f these methods has an extensive l i terature associated wi th it and on ly the barest discussion w i l l be prov ided here. The method o f characteristics exploits the hyperbo l i c nature o f the governing equations to convert the two part ia l d i f ferent ia l equations into four ord inary di f ferent ia l equations. The result ing ord inary d i f ferent ia l equat ions are then integrated by using a finite di f ference technique to approx imate the f r ic t ion term. The great popu lar i ty o f the method is due to its accuracy, eff ic iency and the ease wi th wh i ch it can be p rog rammed. The method has the dist inct advantage o f be ing exp l ic i t — a property that usual ly a l lows each boundary condi ton to be solved independent ly o f the others, w i th large savings i n computer effort as a consequence. In add i t ion, the power and the accuracy o f the method has been adequately con f i rmed by exper imenta l results (eg., Chaudh ry , 1979; Chaudhry , 1981; W y l i e and Streeter, 1983). F i n i t e dif ference methods are obta ined by replac ing the derivat ives in the govern ing equations by dif ference quotients. These methods are h igh ly deve loped both f r om a mathemat ica l and an appl icat ions po in t o f v iew (e.g., Lamber t , 1973; Kre iss, 1978). Exp l i c i t f in i te di f ference schemes usual ly have s ignif icant restict ions on the m a x i m u m t ime step i f stable solut ions are to be achieved. However , imp l i c i t methods, ' T h e first computer techniques s imp ly emulated the prev ious ly avai lable graphical techniques. The commerc ia l software p rogram S U R N A L is an attempt to solve the transient network p rob l em using this approach. 19 wh i ch usual ly overcome the stabi l i ty l imitat ions, require a s imultaneous solut ion for every unknown i n the p r ob l em at each t ime step. The s imultaneous nature o f the so lut ion procedure greatly increases the computat ional effort at every t ime step, but i f much larger t ime steps are a l lowed, the increased effort may be jus t i f i ed . However , for the transient analysis o f hydrau l i c systems, this does not appear to be the case. Wy l i e and Streeter (1983) po in t out that the improvements gained by re lax ing the stabi l i ty cr i ter ia are somewhat i l lusory since the restr ict ive stabi l i ty cr i ter ia associated wi th the exp l ic i t procedures must st i l l be mainta ined i f a satisfactory leve l o f accuracy is to be ach ieved. 7 F o r these reasons, f in i te dif ference techniques can se ldom compete wi th the me thod o f characteristics for the transient analysis o f p ipe l ine systems. M a n y o f the same general conclus ions about f inite di f ference methods also app ly to the f in ite e lement approach. L i k e f in ite difference techniques, the f inite e lement method is h igh ly developed. However , the work o f Wat t et a l . (1980) shows the f in i te element method to have s igni f icant ly increased computer costs over the me thod o f characteristics, w i th the method o f characteristics sti l l setting the standard o f computat iona l accuracy. T h e power o f the f in ite element method comes i n its abi l i ty to accurately represent complex, h i ghe r -d imens i ona l domains. Thus, as long as p ipe l ine transients are represented as one -d imens iona l " waves, i t is h igh ly un l i ke ly that the f in i te e lement method w i l l be a serious compet i tor wi th the method o f characteristics 7 These conclusions may not app ly to parabo l i c part ia l d i f ferent ia l equations such as those descr ib ing transient f low in natural gas networks. In this app l i cat ion the inert ia l terms i n the govern ing equat ions are not important and the f in ite d i f ference techniques are of ten prefered over the characterist ic methods. 8 B y one -d imens i ona l i t is meant that the pos i t ion in any p ipe can be located by f i x ing a single length var iab le measured f r om one end. In such appl icat ions the character ist ic 'surfaces' are s imp ly l ines, and solut ions to these equations are easily obta ined. I f the doma in were even two d imens iona l , the characterist ic surfaces wou ld become planes and eff ic ient appl icat ion o f the characteristics method wou ld be imposs ib le . I f a two or h igher d imens iona l doma in were to be analyzed, it is l ike ly that the f inite e lement and f in i te di f ference techniques wou ld qu i ck l y prove their super ior i ty over the method o f characteristics. 20 in pipe network applications. For all of the above reasons, the method of characteristics will be adopted for the analysis of the hydraulic networks studied in this thesis. The basic equations and solution procedures used in this method will be briefly oudined in the next section. 2.4 METHOD OF CHARACTERISTICS The method of characteristics is presented clearly in a number of references (e.g., Evangelisti, 1969; Gray, 1953; Chaudhry, 1979; Wylie and Streeter, 1983) and only a sketch of the method will be outlined here. The method of characteristics proceeds by combining the dynamic and continuity equations together with an unknown multiplier. Suitably chosen values of this multiplier allows the partial differential terms to be combined together and replaced by total differentials. When using the simplified governing equations the result of this process is dH + _a_ dQ, a f Q | Q | _ Q ( 9 ) with the associated equation Equation (9) can be easily generalized for friction loss formulas other than the X O l O l n ~ * Darcy-Weisbach form shown here by replacing the final term with V | V | for 2 D m appropriate values of X, m and n. The two equations associated with the positive value are usually termed the C + equations and the remaining two relations associated with the negative value are called the C equations. Figure 2 shows how the solution develops. The head and flow values are known at time to and it is desired to know these values at time to + A t A typical such point is P, with unknowns Hp and Q p. The ratio of Ax to At is chosen so the associated relation — — a is satisfied between A and P, and 4 r = -a is dt dt 21 t 0 + 2 A t t 0 + A t / \ A x • A C B x F igu re 2. x - t P lane wi th Character ist ics G r i d satisf ied between B and P. Nex t , Equa t i on (9) is integrated between A and P for the C + equat ion and between B and P for the C ~ equat ion. The d i f f i cu l ty i n this integrat ion l ies solely i n the f r ic t ion term, since the other terms can be exactiy evaluated. The discharge i n the f r ic t ion te rm is usual ly s imply evaluated at A or B where its value is known . Th i s results i n a f irst o rder dif ference so lut ion to Equat ions (9) and (10). T h e first order so lut ion has been p roven suff ic ient ly accurate for a l l but the most severe appl icat ions hav ing large f r i c t iona l losses (quantitat ive l imi ts on the size o f the f r i c t ion term are suggested by Cha udh r y and Ho l l oway , 1984). I f such losses are present, the te rm can be approx imated by a h igher order integrat ion wh i ch overcomes the inaccuracies o f the first order scheme. ' ' I n a recent paper, W y l i e (1983) has proposed a new integrat ion o f the f r ic t ion te rm wh i ch preserves the l inear nature o f the first order technique but wh i ch achieves a h igher order o f accuracy. W y l i e c la ims that this approx imat ion is uncond i t iona l ly stable for a l l f r i c t ion values but the author 's exper ience shows the procedure can produce v io lent instabi l i t ies for even l ow f r ic t ion values. T h e reason for this b reakdown is uncerta in and the technique is avo ided for th is reason. 22 „ Following the notation of Wylie and Streeter the procedure can be written in the following form: C + , H p = C p - B Q p (11) C", H p = C M + BQp (12) in which C P = H A + B Q A - R Q A I Q A I -C M = H B " B Q B + R Q B I Q B I ' B = - | T gA' and f A x R = 2gDA2 ' For internal sections where two characteristics meet, as at point P in Figure 2, Equations (11) and (12) provide a simple solution for Hp. This can be written H p = ( C p + C M ) / 2 (13) Qp can then be found by back substitution into Equation (11) or (12). At each end of the pipe only one of Equations * (11) and (12) will be available. Solution for the unknown head and flow at this point can not then be obtained by the above procedure. Instead, consideration must be given to the nature of the boundary condition, which generally takes the form of a prescribed relation between head and flow. Since there are a large number of possible boundary condtions, this subject is treated in a subsequent chapter. 23 2.5 C O N C L U S I O N S 1. The p rob l em o f ana lyz ing hydrau l i c networks i n the transient state is governed by two quas i - l i nea r hyperbo l i c part ia l d i f ferent ia l equations for wh ich there exists no c l o s ed - f o rm solut ions. 2. The p r imary assumptions made in der iv ing the govern ing equations are that the f low is one -d imens i ona l and isothermal . The convect ive terms are usual ly assumed to be negl ig ib le re lat ive to the other in the govern ing equations. 3. The wavespeed o f the pressure impulse can be readi ly estimated f rom reported values o f f l u id and p ipe l ine propert ies. The discrepancy i n reported celer ity equations for type ' a ' restraint is due to neglect ing the longi tudinal velocity o f the p ipe l ine wa l l i n some derivations. 4. A wide range o f techniques are avai lable for the numer ica l so lut ion o f the govern ing equations. The method o f characteristics is accurate, eff icient, easily p rog rammed and, in general , we l l suited to one -d imens i ona l network analysis. 5. A wel l posed mathemat ica l p rob lem for a first order system o f hyperbo l i c part ia l differentia] equations requires not on ly descr ipt ion o f the governing equations, but descr ipt ion o f the in i t ia l and boundary condit ions as we l l . These topics are treated in the fo l l ow ing two chapters. 3. S T E A D Y - S T A T E N F T W O R K A N A L Y S T S 3.1 INTRODUCTION The p rob l em o f determin ing the steady-state d istr ibut ion o f pressure and f low in a complex f l u id network is o f great importance i n engineering. D i s t r i bu t ion systems are used to convey a var iety o f f lu ids for numerous end uses. Wa te r d istr ibut ion systems for large cit ies represent some o f the most sophist icated hydrau l i c networks in existence. T h e capital investment and operat ing expenses associated w i th these networks are considerable, mak ing carefu l p lann ing desirable before a new system is instal led or before s igni f icant changes to exist ing systems are undertaken. A prerequis i te to good design, operat ion or opt imizat ion o f hydrau l i c networks is the ab i l i ty to pe r fo rm accurate analysis o f exist ing or proposed distr ibut ion systems. H is tor i ca l l y a lmost a l l the effort app l ied to large network analysis has been exerted to accurately pred ic t the f lows and pressures throughout the network under steady-state condit ions. U sua l l y the goal o f such analysis is to insure that an adequate m i n i m u m pressure is ma in ta ined at a l l locat ions wh i l e supp ly ing speci f ied flows at prescr ibed locations i n the network. T h e extensive and growing l i terature that exists per ta in ing to steady-state network analysis is test imony to both the impor tance o f this p r ob l em and to the complex i ty o f its solut ion. The int roduct ion o f the computer has created a revo lut ion i n this area o f hydrau l i c engineer ing as it has i n many others, w i th a great many new techniques be ing introduced that d i d not exist twenty years ago. However , there sti l l exists much controversy i n the l i terature concerning wh ich analysis technique is 'best' f r om the v iewpo int o f accuracy, re l iab i l i ty , speed, storage requi rements and ease o f p rog ramming . Recent ly , many attempts have been made to look not only at the p rob l em o f analysis, but also at the p r ob l em o f design opt imizat ion. Nume rous references attempt to prov ide techniques for op t im i z ing network design under a var iety 24 25 o f constraints and assumptions. N o t on ly is steady-state network analysis an important p rob l em in its own right, but i t is also an essential first step in ana lyz ing a d is t r ibut ion system under transient condit ions. In order to have a we l l - po sed transient p rob lem, the hydrau l i c gradel ine elevations and f low velocit ies must be spec i f ied throughout the network at the in i t ia l t ime. These in i t ia l condi t ions are then mod i f i ed by the mechanics o f wave propagat ion w i th in the network and by the interact ion o f the transient waves w i th the network boundary conditons. Accurate speci f icat ion o f the transient condi t ions is c lear ly imposs ib le i f the in i t ia l condit ions are i n considerable error. Fu r the rmore , it is h igh l y desirable to direct ly couple the determinat ion o f the steady-state f lows and pressures in a network w i th the subsequent transient analysis. F r o m a user po in t o f v iew such a coup l ing s impl i f ies data input by requ i r ing the preparat ion o f on ly one data file — there is no need to labor ious ly feed the output o f a steady-state program back i n as input for the transient computat ion. F r o m a p rog ramming po in t o f v iew, many o f the requ i red data checks, in i t ia l i zat ion o f data structures and pre l im inary computat ion needs to take place on ly once, result ing i n a shorter, s imp ler and more eff ic ient p rogram. In addi t ion, the abi l i ty to accurately pred ic t the steady-state condi t ions f r om w i th in the transient calculat ions has the potent ia l for p rov id ing important ins ights in to the response and behav iour o f the network as a whole. Such ins ight cou ld have s ignif icant design and operat ional impl icat ions. Th i s aspect w i l l be discussed i n a subsequent chapter, fo l low ing the deve lopment o f the transient network mode l . T h e purpose o f this chapter is to out l ine the nature o f the steady-state analysis p rob l em, the solut ions p roposed i n the l iterature, and the characterist ics o f the adopted so lut ion technique. The chapter concludes wi th some imp lementat ion details necessary for the eff ic ient so lut ion o f the transient p rob lem. 26 3.2 S T A T E M E N T O F T H E P R O B L E M A d is t r ibut ion system must be capable o f supp ly ing peak demands wi th adequate res idual pressures wh i le convey ing water f rom sources o f supp ly to demand centers. The least cost so lut ion to this p r ob l em is to connect the demand centers to the supply sources by a path o f m i n i m u m length (a m in ima l spanning tree). A t ree -shaped structure, such as the one shown in F i gu re 3, is easy to analyze hydrau l ica l ly , since the f low i n any one p ipe is s imp ly the sum o f the downstream demands, and the pressure fo l lows immedia te ly once the f low is speci f ied. Howeve r , for the fo l low ing pract ica l reasons water d is t r ibut ion systems are not designed as m i n ima l spanning trees: • dead end pockets where the water may become stagnant are undesirable f r om a water qual i ty po in t o f v iew; • rout ine maintenance or emergency repairs cut o f f supply to a l l users downstream o f the interrupt ion; • the tree structure does not have the requ i red f lex ib i l i ty to meet large f lows at unspec i f ied locations (such as fire f lows). Wa te r d is t r ibut ion systems therefore have a complex topology w i th numerous redundant p ipes in t roduced to overcome the above objections. Th i s new structure is shown in F i gu re 4. However , the redundant p ipes greatly compl icate hydrau l i c analysis. I N ' J F igu re 3. Ne two rk w i th T r e e - Shaped Structure 27 F igu re 4. Looped Ne two rk Topo logy De te rm in i ng the p ipe f low rates and junct ion heads i n the d is t r ibut ion system now requires a large system o f non l inear equations to be solved by numer ica l means. The steady-state hydrau l i c network analysis p rob l em consists o f f i nd ing the f low in a l l p ipes, the hydrau l i c gradel ine elevat ion at a l l nodes and the external f lows at a l l nodes w i th a spec i f ied h e a d - f l o w re lat ionship (pumps, reservoirs, or i f ices, etc.). F o r the purposes o f this thesis the steady-state p r ob l em wi l l be assumed to be o f the fo l l ow ing f o rm: • the phys ica l descr ipt ion o f the network is known ; • the network may have a complex l ooped or t r ee - shaped topology w i th consumpt ion or supply at any node i n the network; • special node boundary condi t ions may exist where external f lows are related to the hydrau l i c gradel ine e levat ion at the node; • and special " i n - l i n e " boundary condi t ions may exist where the pressure loss or gain i n a l ine is a funct ion o f f low in the l ine. 28 The phys ica l descr ipt ion o f the network inc ludes the layout o f pipes as we l l as their diameter, length and loss coefficients. A l l p ipes not designated as 'specia l ' are assumed to behave accord ing to standard f r i c t ion loss formulas such as represented by the D a r c y - W e i s b a c h , M a n n i n g or H a z e n - W i l l i a m s formulas. A l s o assumed to be known is the phys ica l arrangement o f a l l pipes, pumps , valves and other appurtenances in the network. The on ly restr ict ions on network topology is that a p ipe is not a l lowed to j o i n a node w i th i tse l f (i.e., a single p ipe loop) . Mu l t i p l e p ipes between nodes and 3 -d imens iona l loops (where p ipes cross but do not j o i n ) are a l lowed. In addi t ion, i f it is cons idered desirable, a number o f unconnected networks, or port ions o f a network, can be analyzed s imultaneously. N o d e and p ipe number ing is arb i t rary except that the m a x i m u m assigned number cannot exceed the d imens ion o f the node or p ipe array. Specia l node boundary condit ions are in t roduced to represent reservoirs, pumps, or i f ices and other devices wh i ch al low external f low f r om or to the network to be computed as a funct ion o f junc t ion head. These nodes can be arbi trar i ly located throughout the network. " I n - l i n e " boundary condi t ions a l low booster pumps, pressure reducing valves ( PRV ' s ) , check valves, local losses as wel l as var ious other hydrau l i c controls to be represented. In these devices, there is an assumed un ique re lat ionship between f low in a l ine and the head change encountered by travers ing that l ine. 3.3 B A S I C E Q U A T I O N S In order to solve the steady-state network p rob l em, enough equat ions must be written to a l low determinat ion o f al l the unknowns. Essent ia l ly, al l networks are solved by app ly ing the two fundamenta l phys ica l laws o f the conservat ion o f mass and the conservat ion o f energy. W h e n app l ied to a hydrau l i c network, these two laws can be expressed i n the fo l l ow ing f o rm: 29 • conservat ion o f mass - the net sum o f a l l the f lows at a junc t ion must equal zero; • conservat ion o f energy - the net change in head around any closed loop must equal zero. These equat ions are sometimes refered to as K i r c h h o f f s first and second law in analogy w i th electr ical networks. In add i t ion to these equations, the head at at least one locat ion must be a known funct ion o f the external f low at that locat ion for the pressures to be un ique ly determined w i th in the network. Each special node boundary cond i t ion o f the type descr ibed above w i l l introduce an add i t iona l unknown (the external f low) and an addi t ional equat ion (the external f l ow/ junc t ion head relat ion). Thus , special nodes do not inf luence the balance between the number o f equations and the number o f unknowns. I f the head change across any l ine is known as a funct ion o f the f low in that l ine, the result ing system o f equations has been shown to have a so lut ion. T h e so lut ion has been shown to be unique for s imple systems conta in ing p ipes w i th f r ic t ion loss, pumps and check valves (Warga, 1954). However , when P R V ' s and other more complex devices are introducted into the system, the quest ion o f uniqueness has not been finally demonstrated, a l though no examples o f nonuniqueness have been reported (Gessler, 1981). The above general discussion can be given a somewhat more precise mathemat ica l fo rmulat ion. Le t the number o f nodes in a network be n, the number o f p ipes be m and the number o f special node boundary condi t ions be s. The fo l l ow ing var iables are unknown: • m interna l f lows i n the m p ipes; • s - 1 external f lows at the special node boundary condi t ions (one external f low is not independent since it can be expressed in terms o f the sum o f the known demands at the ord inary nodes and the sum o f the other s - 1 external f lows); • n - s heads at the ord inary network nodes. 30 In total there are m + n - 1 unknowns to be determined. T o solve for these unknowns there are m (un ique) re lat ionships between the f low and head change i n a l ine and n - 1 independent node balance equations. Thus the number o f independent basic equations equals the number o f independent unknowns. Th i s a l lows the network to be represented by a system o f m + n - 1 equations. Unfor tunate ly , this system o f equations is nonl inear . The re lat ionship between f low in a p ipe and the change i n head across that p ipe is a lmost a lways nonl inear. F r i c t i on losses i n p ipes are usual ly mode led i n the f o rm h . = K - C v l Q j l 1 1 - 1 : (1) i n wh i ch h . is the head loss in l ine i, Q . is the discharge i n l ine i and K . and n I i E I are parameters depending on the p ipe propert ies, the units used and the part icu lar f o rm to the loss equat ion adopted (e.g. Mann i ng , H a z e n - W i l l i a m s ' or D a r c y - W e i s b a c h formulas). In add i t ion the head changes across pumps, valves and other hydrau l i c contro l elements are usual ly nonl inear funct ions o f the f low, often approx imated by piecewise po lynomia l s over the operat ing range o f interest. W h e n the l ine nonl inear i t ies comb ine wi th the frequent ly non l inear re lat ion between the junc t ion head and external f l ow present at special nodes, i t is clear that the result ing system o f equations is non l inear for wh ich there is no direct so lut ion. The no rma l procedure for so lv ing a nonl inear system o f equations is to l inear ize the p rob lem, solve the result ing l inear system, and then improve on the l inear approx imat ion through iterat ion. A l l the techniques proposed for solv ing hydrau l i c networks i n the steady-state use essentially this technique, the differences coming i n the nature o f the l inear izat ion and the details o f how the resul t ing system is solved. These solut ion procedures w i l l be br ief ly rev iewed in the next section fo l lowed by a discussion o f the adopted so lut ion technique. 31 3.4 S O L U T I O N T E C H N I Q U E S It has been shown how the nonl inear system o f equations representing a hydrau l i c network i n steady-state is developed. The purpose o f this section is to out l ine br ief ly the techniques that have been proposed for so lv ing this system o f equations. The procedures out l ined are o f relevance to dig i ta l network solut ion and no ment ion w i l l be made o f techniques such as analogue computat ion wh ich have been al l but total ly replaced by digi ta l techniques (a good review o f these approaches is found in Raman ' s 1970 paper). A few miscel laneous techniques have been proposed but have not gained widespread favour. F o r example, Co l l i n s and Johnson (1975) proposed using f in ite elements to solve the network p rob l em and F u k u d a (1981) suggested using parametr ic analysis o f d i f ferent ia l equations der ived f r om the imp l i c i t funct ion theorem for nonl inear equations. These approaches are interest ing but w i l l not be discussed since they p rov ide no s ignif icant advantages over the mo r e standard techniques o f steady-state analysis. Cons ide r a network consist ing o f m p ipes and n nodes wh ich can be d iv ided into / c losed loops. A s imple result f r om graph theory (easily proved by induct ion) is that m = / + ( « - 1) (2) Whereas the mass conservat ion equations at the nodes and the h ead - f l ow relat ionships i n the p ipes suggest that it is necessary to solve a complete system o f m + n - 1 equations, app l i cat ion o f this graph re lat ion a l lows the number o f unknowns so lved for direct ly to be greatly reduced. B y var ious techniques the number o f equations to be so lved can be reduced to any o f the number o f p ipes (m), the number o f junct ions less one (n - 1) or to the number o f loops (/). A l t hough there are many variat ions i n the specif ics o f the algor i thms, there are essential ly f ive techniques wide ly advocated i n the l i terature for the steady-state so lut ion o f hydrau l i c networks. T w o o f these 32 procedures are based on so lv ing the / loop equations, two are based on so lv ing the n - 1 junc t i on equations and one procedure solves for the m unknown p ipe flows. A good elementary introduct ion to these techniques is p rov ided be Jeppson (1976). 3.4.1 T H E P I P E F L O W E Q U A T I O N S In 1972, W o o d and Char l es presented a new technique for p ipe network analysis based on lett ing the flow i n each p ipe be the system unknowns. Since there are m p ipes in the network, m independent equations need to be writ ten. The equations chosen are the n - 1 cont inu i ty equations at the nodes and the / loop equations, the total number o f equations equal l ing the requ i red number by the graph re lat ion shown above. Th i s system is so lved by l inear iz ing the h e a d - f l o w relat ionship in each p ipe and, since the junc t ion equations are l inear, p roduc ing a l inear system o f equat ions wh i ch can be easily so lved. F o r example, the f r i c t ion loss re lat ion i n a p ipe can be l inear ized by wr i t ing the equat ion as h . = k.Q (3) i n wh i ch k. = K . Q - . I f we have a estimate for Q . , the value o f k. can I l ' ^ i 1 i I be evaluated and a l inear equat ion results. I f a l l the non l inear elements in the network can be s imi la r ly represented, a l inear system results wh i ch can be solved for a l l the p ipe flows. These new p ipe flows can be used to update the value o f k. and the process repeated unt i l the flow rates converge to a prescr ibed tolerance. Once the flows are accurately known , the heads at the nodes can be easily determined f r om the known h e a d - f l o w re lat ionships i n each p ipe. Despi te appearances, in i t ia l values o f do not need to be p rov ided since the in i t ia l value o f ytj can be s imply est imated as K^. E f f i c ient and re l iab le convergence is reported for this procedure as we l l as good stabil ity propert ies. The l inear system to be solved is larger than that p roduced by other techniques in common use, 33 but qu ick so lut ion is possib le by taking advantage o f the sparcity and symmetry o f the result ing l inear system. 3.4.2 THE NODE EQUATIONS I f a head is in i t ia l ly assigned to each node i n the network, the loop equations w i l l be automat ica l ly satisfied and the number o f equations to be so lved can be reduced f r om m to n - 1. Once the heads are known or assumed, the f low in each p ipe can be easily calculated. However , the in i t ia l head assignment w i l l general ly mean that the cont inuity equat ion is not satisf ied at each node. Correct ions can be appl ied to the value o f head at each node so that the imbalance i n the mass conservat ion equat ion, over a number o f iterations, is made arb i t rar i ly smal l . Iterations are necessary since the head correct ions w i l l be based on a l inear ized vers ion o f the cont inui ty equat ion (expressed i n terms o f head losses or gains i n the p ipe members) . The two techniques that ut i l ize the f low balancing technique are essential ly ident ical i n formulat ion but d i f fer i n the way head correct ions are calculated and app l ied to each node i n the network. Ha rdy Cross (1936) was the first to suggest this "method o f ba lanc ing f l ows" as a pract ical so lut ion technique. Cross proposed that correct ions be calculated on a node by node basis. Tha t is, at each node we consider that we have a single non l inear equat ion (the cont inui ty equat ion) i n one unknown (the node head). The root o f this equat ion is approx imated us ing Newton ' s method wh i l e at the same t ime ignor ing the interact ion o f this node's head wi th the other nodes and elements i n the system. F o r the smal l networks to wh ich Cross was app ly ing his procedure, convergence was cons idered to be satisfactory. M u c h later, when digital computers had vastly increased our abi l i ty to pe r f o rm large scale computat ions, Shami r and H o w a r d (1968) developed an 34 improvement to the tradi t ional Ha rdy Cross method. They reasoned that the convergence rate cou ld be greatly increased i f the calculated correct ions to node heads took into account the interact ion between adjacent nodes. Th i s formulat ion resulted m a n - 1 by « — 1 system o f nonl inear equations wh ich they then solved s imultaneously by using the N e w t o n - Raphson approach. Th i s procedure represents a vast improvement in convergence rate over the Ha r d y Cross method, a l though the authors noted that convergence was sometimes s low i f the in i t ia l assignment o f heads was poor. Later wri ters have also noted convergence prob lems when using this technique. 3.4.3 T H E L O O P E Q U A T I O N S I f an in i t ia l f l ow is assigned to each p ipe in the network such that the junc t ion cont inui ty equations are ident ical ly satisf ied, the number o f equations to be so lved can be decreased st i l l futher. Once in i t ia l f lows are assigned, it w i l l be found that, i n general , the sum o f the (now known) head changes around each loop is not zero. Correct ions can then be app l ied to the f l ow . i n each p ipe i n the loop i n order to br ing this nonsatisfact ion o f the energy equat ion arb i t rar i ly close to zero. Th i s correct ion w i l l . i n general not be exact, be ing calculated f r om a l inear ized version o f the loop balance equat ion, and iterations w i l l thus be requ i red i n order to achieve convergence. It shou ld be noted that the f low correct ion does not mod i fy the mass balance at a node, leav ing the junc t ion cont inui ty equations ident ical ly sat isf ied throughout the iterat ion process. The two solut ion techniques that ut i l ize the loop balance approach are identical i n their fo rmu la t ion but di f fer s igni f icant ly i n the way the / loop f low correct ions are calculated and appl ied to each loop. The first so lut ion technique suggesting the use o f the loop equations was presented by Ha r d y Cross i n 1936. Essent ia l ly the procedure consists i n so lv ing a 35 single non l inear equat ion using Newton ' s method as was done wi th the junct ion equations. Howeve r , now an unknown loop flow rate correct ion is app l ied in order to correct the imbalance i n the loop energy equat ion. Th i s "method o f ba lanc ing heads," l i ke Cross 's other method, ignores the interact ion between adjacent system components. F o r many years this method was the on ly pract ical one ava i lab le for obta in ing the steady-state distr ibut ion o f flows in a looped network and it is st i l l synonymous i n many engineers' m inds w i th hydrau l i c network analysis i n general . However , the Ha r d y Cross technique has been found to have s low or non-ex i s ten t convergence when appl ied to very large systems, a l though convergence was found to be adequate on the or ig ina l smal l networks tested by Cross. A s w i th the node approach, the poor convergence o f the Ha r d y Cross procedure can be traced to ignor ing system interactions wh i l e calculat ing the s ingle head or f low corrections. E p p and Fow le r (1970) developed the s imul taneous loop so lut ion approach by app ly ing the N e w t o n - R a p h s o n technique to the system o f / non l inear equations obta ined by app ly ing the conservat ion o f energy equat ion to each loop i n the network. In this formulat ion, the f low rate correct ion i n each loop is considered as a vector o f unknowns. The f low rate correct ions are then calculated by so lv ing a system o f equations formulated to take in to account adjacent loop interactions. The equations are l inear ized by eva luat ing the derivat ives requ i red by the N e w t o n - R a p h s o n method wi th the k nown or assumed f low rates i n the p ipe members. Convergence has been found to be excel lent when good in i t ia l approx imat ions o f the flow rates are suppl ied. Th i s technique w i l l be discussed i n more detai l later. 36 3.5 C O M P A R I S O N O F T H F S O L U T I O N T R C H N T O U E S There has been considerable argument i n the l i terature as to wh i ch analysis technique is preferab le for the steady-state ca lculat ion o f f lows and pressures i n d istr ibut ion networks. Typ i ca l l y authors adopt one approach and al l alternative approachs are br ief ly d ismissed. O n l y a few writers have attempted a careful and comprehens ive compar ison wi th the goal o f determin ing the accuracy and re l iab i l i ty o f the var ious a lgor i thms. The most important o f these studies is that o f W o o d and Rayes (1981). W o o d and Rayes studied a l l f ive o f the solut ion techniques descr ibed above wi th the goal o f determin ing how re l iab ly the techniques predicted the f lows and the pressures i n a number o f test networks. The solutions were cons idered acceptable i f the average deviat ion between the exact and calculated solut ion i n either heads or f lows was less than 10%, p rov ided that the m a x i m u m deviat ion was less than 30%. T h e exact so lut ion was considered to be the calculated solut ion based on the l inear theory me thod when using a much more str ingent convergence cr i ter ia. A total o f 60 networks were ana lyzed for systems hav ing under 100 pipes and 31 networks hav ing more than 100 pipes. W o o d and Rayes conc luded that both o f the techniques based on the noda l fo rmu la t ion and the l oop -based Ha r d y Cross method are unsuitable as general analysis tools. These approaches "exh ib i ted s igni f icant convergence problems, and the frequency o f the p rob lems increased as larger systems were analyzed." In contrast, the s imultaneous loop formulat ion and the l inear theory method had "excel lent convergence characterist ics." Perhaps the most a la rm ing o f their conclusions is that, in the techniques that fa i led to converge, the fai lures were a lmost imposs ib le to detect. Tha t is, by a l l ind icat ions the system appeared to have converged but the calculated values d i d not meet the selected accuracy cr i ter ia. W o o d and Rayes suggest that this fa i lure is p robab ly due to smal l errors i n the pressure values result ing in much larger errors i n the f low values i n some o f the networks studied. 37 , In this thesis, the steady-state so lut ion is calculated p r imar i l y i n order to establ ish the in i t ia l condit ions requ i red by subsequent transient analysis. D u e to the •mult ipl icat ive effect o f the wavespeed, transient pressures are very sensit ive to assigned p ipe f low rates. F o r this reason it is essential that the approach adopted for the steady-state analysis be as accurate and re l iab le as possible. Gess le r and others have po in ted out that a h igh degree o f pressure accuracy does not guarantee a correspondingly h igh degree o f f low accuracy. Thu s convergence cr i ter ia must be based on f low as the pr imary unknown i f the transient pressures are to be re l iab ly predicted. W h i l e the two methods o f analysis proposed by Cross represented a great advance, numerous studies have demonstrated their unre l iab i l i ty when app l ied to large and complex networks. M o r e power fu l and re l iab le a lgor i thms have subsequently been developed and, wi th the trend to increasing computer power and decreasing computer costs, the object ions to these more sophist icated approaches are becoming progressively more indefens ib le. A l t hough a great deal o f effort has been expended to improve the basic single equat ion a lgor i thms, it is c lear that these methods w i l l at best on ly approx imate the s imultaneous correct ion techniques. It is past t ime that the Ha r d y Cross methods be f ina l ly la id aside. Thus , some o f the a lgor i thms avai lable for so lut ion o f the steady-state network p rob l em are by modern standards inef f ic ient or unre l iab le when app l ied to large networks. Other , more nove l , solut ions have not gained widespread acceptance. The two techniques most appl icable today are the s imultaneous loop approach and the l inear theory m e t h o d Bo th o f these techniques are re l iab le but the quoted graph re lat ion shows that the s imultaneous loop technique w i l l result i n a s igni f icant ly smal ler l inear system (since / is much less than m). F o r this reason, and since an eff ic ient ly written source l is t ing was readi ly avai lable, the s imultaneous loop equat ion approach was adopted to determine the steady-state network solut ion. The characterist ics o f the adopted so lut ion technique w i l l be br ie f ly discussed i n the next section. 38 3.6 C H A R A C T E R I S T I C S O F T H F . A D O P T E D S O L U T I O N T E C H N I Q U E The steady-state analysis technique adopted in this thesis is based on the work by E p p and Fow l e r (1970). The a lgor i thm deve loped by E p p and Fow l e r is eff icient, re l iable and overcomes almost a l l o f the object ions put forward concerning the l oop -based approach by those who prefer other analysis techniques. Recent ly , Oh tme r (1981) has descr ibed a so lut ion procedure a lmost ident ical to the approach proposed by E p p and Fow l e r twelve years earl ier. T he pr imary characteristics o f the so lut ion technique are l is ted br ief ly below. Some addit ional detai ls are p rov ided i n subsequent chapters where such in format ion relates to the imp lementat ion or p rogramming features o f the transient solut ion. 1. The a lgor i thm uses the N e w t o n - R a p h s o n method to solve a system o f s imultaneous nonl inear equations. Init ia l f lows satisfying the mass balance equations are automatical ly assigned using a m in ima l spanning tree a lgor i thm wh ich , i n practice, assures the convergence o f the iterative so lut ion technique. 2. T he solut ion is l oop -o r i en ted i n order to reduce the number o f equations to be solved and to take advantage o f the better convergence propert ies associated wi th this approach. The loops are ident i f ied automat ica l ly by the p rogram wh ich keeps the requ i red input data to a m i n imum . In addi t ion, branched appendages are 'purged ' p r io r to network solut ion further reduc ing the execution t ime. 3. A n automat ic loop renumber ing a lgor i thm is employed to produce a banded symmetr i c matr ix. Th i s a l lows computer storage to be conserved and also a l lows eff ic ient rout ines such as Cho lesky decompos i t ion to be used i n solv ing the system o f l inear ized equations. The a lgor i thm has been wide ly tested and been shown to be an eff ic ient and rel iable way o f obta in ing the steady-state so lut ion o f large and complex f l u i d networks. Thus, the s imultaneous loop approach o f E p p and Fow l e r represents an excel lent way o f obta in ing the in i t ia l condi t ions requ i red for subsequent transient analysis. 4. B O U N D A R Y CONDITIONS 4.1 I N T R O D U C T I O N In Chap te r 2 o f this thesis the basic part ia l di f ferent ia l equat ions govern ing the unsteady f low o f f l u i d i n c losed-condu i t s were developed. The method o f characteristics was app l ied i n order to calculate the so lut ion to these equations on the discret ized x - t plane. Howeve r , the complete so lut ion o f the governing equat ions also requires specif icat ion o f the condit ions at the beg inn ing o f the transient r un and calculat ion o f the boundary condi t ions present i n a network. A procedure for eff ic ient ly obta in ing the in i t ia l condi t ions has been descr ibed i n Chap te r 3; the rat ional so lut ion o f the boundary condi t ions is the subject o f the present chapter. The prev ious two chapters have emphas ized recent developments in the study o f hydrau l i c networks. A l t hough frequent reference is made to the avai lable l i terature, the pr imary purpose o f this chapter is to improve the current techniques o f so lv ing the boundary condi t ions common ly found i n hydrau l i c networks. In recent t ime, the field o f transient analysis has been dominated by s imp le series p ipe l ine appl icat ions, wh i le the field o f network analysis has been preoccup ied wi th obta in ing steady-state network solutions. Th i s chapter seeks to develop general techniques o f so lv ing the . network boundary condi t ions i n order to predict their transient behavior and response. The subject o f what constitutes a boundary cond i t ion is treated quite general ly i n this thesis. Chap te r 2 describes how the time domain is discret ized in to A t segments i n order to apply the method o f character ist ics. 1 0 A s a consequence o f the choice o f A t , most p ipes i n the network w i l l be d iv ided into a numbe r o f reaches o f length A x . T h e end o f each reach, where head and f low values must be calculated til or speci f ied, w i l l be cal led a section. I f the subscr ipt Y is used to designate the i " T h e quest ion o f how to select the A t va lue is discussed i n Chap te r 5 under the topic o f network discret izat ion. 39 40 pipe, then the number o f reaches i n p ipe i is N R j and the number o f sections is N S j . 1 1 T he sections present i n any p ipe can be d iv ided into two groups: boundary sections are located at each end o f the physical p ipe and exactly two wi l l be present i n every p ipe; the remain ing sections are termed internal sections and any number may be present The head and discharge at internal sections can be computed by stra ightforward appl icat ion o f Equa t i on (13) in Chap te r 2. However , values at the boundary sections can not be obta ined unt i l an aux i l i a ry re lat ion between head and discharge at the end o f each p ipe is specif ied. Such a head-d i scharge re lat ion at any boundary section w i l l be termed a boundary condition i n this thesis. T he number o f dist inct boundary condit ions that may exist in a general hydrau l i c network is very large. Networks may contain a host o f control and dist r ibut ion devices connected together i n complex arrangements wh i ch seem to defy general izat ion. E f f i c ient so lut ion o f these boundary condit ions requires a classif icat ion system wh i ch can exploi t the s imi lar i t ies between the var ious network elements wh i le st i l l ma in ta in ing their dist inct ive features. Thus, this chapter begins by descr ib ing a boundary cond i t ion classif icat ion system wh ich al lows the behav ior o f al l hydrau l i c devices common l y found in d ist r ibut ion systems to be treated i n a general manner. F o l l ow i ng this, solut ions o f several o f the common boundary condi t ions are either expl ic i t ly obta ined or put i n a f ramework wh i ch a l lows eff ic ient numer ica l calculat ion. T o the author 's knowledge, this general ized character izat ion o f the network boundary condi t ions is unique to this thesis. " W i t h i n a p ipe — where, by def in i t ion, a l l parameters are constant — a s imple compat ib i l i ty re lat ion can be app l ied . Th i s compat ib i l i ty re lat ion states that, for internal sections, the head and the f low are uniquely determined. W h e n this re lat ion is app l ied , the numbe r o f sections is equal to the number o f reaches p lus one. I f m i d - p i p e boundary condit ions are a l lowed, such as caused by water co l umn separation occurr ing w i th in a p ipe, this compat ib i l i ty condi t ion no longer appl ies, and the number o f sections is then equal to twice the number o f reaches. Thus , internal p ipe boundary condi t ions greatly increase both the storage and computat iona l requirements o f the solut ion procedure. 41 4.2 B O U N D A R Y C O N D I T I O N C L A S S I F I C A T I O N S Y S T E M L o n g lists o f hydrau l i c boundary condit ions occur frequent ly i n the l i terature perta in ing to transients, and each o f these cases must be ind iv idua l l y solved wh i le per fo rming transient analysis o f p ipe l ine systems. F o r s imple appl icat ions, few, i f any restrictions result f r om this ca se -by - case approach. However , for large networks o f arbitrary structure and arrangement, unnecessary and awkward restraints may be imposed on the ana lys t These restraints may conf ine analysis to s imple networks or they may force the analyst to undertake compl icated manipu lat ions i n order to contr ive the desired boundary condi t ions f r om the exist ing ones. In extreme cases, the network w i l l be altered i n order to apply the restricted set o f avai lable boundary condit ions, general ly w i th no clear idea o f the effect the changes i n the network w i l l have on the analysis. Thus , for example, F o x (1972) requires that a f i v e -p i pe junc t ion be spl i t into two th ree -p i pe junct ions separated by a A x length, wh i le the commerc ia l p rogram S U R N A L a l lows air chambers and booster pumps to occur at the junc t ion o f two pipes only. Unnecessary and troublesome restrict ions such as these are overcome by the general treatment presented i n this chapter. Before a c lassif icat ion system appl icable to a l l boundary condi t ions can be undertaken, a number o f terms need to be def ined. Prev ious ly in this thesis the term pipe has been used wi th the general ly accepted mean ing o f "a c losed-condu i t convey ing a f l u id . " However , there are dist inct advantages in d ist inguish ing very short p ipes f r om longer ones dur ing the analysis. F o r instance, the method o f characteristics general ly requires the time step to be equal throughout the network. I f very short conduits are present, this may make the analysis uneconomica l ly expensive, since even the short condui ts must contain at least one characterist ic reach. It is often desirable to replace these short p ipes w i th approx imate mathemat ica l representations and use a larger t ime step than wou ld otherwise be possible. These p ipe replacement elements (abbreviated P R E ) must be treated separately i n the analysis. T he te rm p ipe w i l l 42 henceforth be restr icted to only those conduits wh i ch contain at least one characterist ic reach. A numbe r o f terms re lat ing to nodes must also be def ined. F i rs t , a node is s imp ly a fixed locat ion where boundary sections mee t The degree o f a node is the number o f p ipes (i.e., characterist ic sections) connected to it Howeve r , i n a general network, not on ly p ipes may be connected to a node, but var ious other elements as we l l . F o r example , a node may represent the suct ion or discharge f lange o f a pump, the locat ion o f a valve d ischarg ing f r om the network or a po in t to wh i ch a pressure re l ie f va lve is connected. A l l such n o n - p i p e connect ions at a node w i l l be labe l led external and the number o f external connect ions w i l l be cal led the complexity o f the node. A node o f complex i ty zero w i l l be ca l led simple, a node o f complex i ty one ordinary, and a node o f complex i ty greater than one complex Since the above termino logy is new, an example w i l l be presented to c lar i fy the concepts invo lved and to present the shorthand notat ion used to describe nodes in this thesis. F i gu re 5 shows a junc t ion o f four p ipes wi th no external f lows or ig inat ing f rom the node. Thus, th is is a s imple node o f degree four and complex i ty zero and - 6 F i gu re 5. S imp le N o d e 43 0 2 w i l l be designated as a node. F i gu re 6 depicts a comp lex node o f degree f ive and complex i ty two. The node is complex since an external f l ow Q and a p ipe rep lacement e lement both connect to it No te that each f low i n the network must be g iven an or ientat ion but that the node termino logy is independent o f whether the assumed di rect ion is towards or away f r om the node i n quest ion. The termino logy related to nodes can be extended to networks as we l l . A network consist ing o f on ly s imple nodes is t r iv ia l and w i l l not be further considered (the on ly so lut ion is zero f low in a l l p ipes). Ne tworks conta in ing on ly nodes o f complex i ty zero or one w i l l be designated ordinary wh i le networks hav ing one or more nodes o f complex i ty two or greater w i l l be ca l led complex. The context w i l l make it c lear whether an ind iv idua l node or a who le network is be ing cons idered. It shou ld be noted that the above designations may be a property o f the network but may also be s imp ly a property o f the analysis tool adopted to study i t F o r example, p ipe rep lacement elements general ly increase the complex i ty o f the network w i thout changing the network topology. F i gu re 6. C o m p l e x Node 44 -Junct ions o f several pipes i n hydrau l i c networks are a lmost universal ly mode led as fr ict ionless in transient f low appl icat ions. The prob lems encountered by attempting to calculate junc t ion losses are very considerable and w i l l be discussed later in this chapter. I f the usual assumption is made that junct ion losses are negl ig ible, a number o f important results can be qu ick ly obta ined. Un less speci f ica l ly stated to the contrary, nodes w i l l be assumed to be fr ict ionless throughout the remainder o f this thesis. In general , the di f f icul ty o f so lv ing a network increases as the complex i ty o f the nodes i n the network increase but, w i th proper boundary cond i t ion formulat ion, is independent o f the degree o f any node in the network. Th i s s imple fact is either not recognized or, at least, not stated i n the hydrau l i c transient l i terature. 1 2 However , the effort that must be expended to solve a network does not depend on the noda l complex i ty alone, but also on the detai ls o f how the network is connected. Before an attempt can be made to characterize the network itself, a classif icat ion system must be deve loped for external f lows. Externa l f lows or connect ions can be classif ied on the basis o f how many noda l values must be known in order to def ine them complete ly . The s implest external f lows depend on ly on the condit ions at a single node. Examp les o f one -node boundary condit ions inc lude a valve or or i f i ce d ischarging to the atmosphere or a constant head reservoir located at the junct ion o f several pipes. In a s imi lar manner, external f lows can depend on condi t ions at two nodes. Th i s is a very c ommon boundary cond i t ion for represent ing booster pumps, i n - l i n e loss elements such as valves or m inor losses, special contro l valves such as pressure reduc ing valves and many other devices. In these two -node elements, the f low does not depend on head values at a single node, but rather on the difference i n hydrau l i c gradel ine elevat ion between the upstream and " A c t u a l l y , many authors incorrect ly assert that the d i f f i cu l ty in solv ing a network depends direct ly on the degree o f the nodes it contains and ignore the concept o f complex i ty as deve loped here. 45 downstream ends o f the e l emen t W i t h the increasing complex i ty o f hydrau l i c control systems an even h igher degree o f nodal dependence may be possib le. Thus , the head d i f ference between two or more nodes cou ld be used to change a performance parameter at a th i rd or even a fourth node. Obv ious l y such general izat ion can go on indef in i te ly i n this c lass i f icat ion system. F o r convenience, the numbe r o f nodes on wh i ch an external f low depends w i l l be ca l led its nodal dependence number and w i l l be abbrev iated N D N . The author is aware o f other ways o f c lassi fy ing boundary condit ions, but these are not as power fu l as the nodal c lassif icat ion system presented above. F o r example, some authors (e.g., F ox , 1977; Koe l l e , 1982; Watters, 1984; W y l i e and Streeter, 1983) and some programs (e.g., S U R N A L , 1977), classify many two -node elements as ' i n - l i n e ' boundary condit ions. Tha t is, devices such as booster pumps, contro l valves or air chambers must be p laced at the abutt ing ends o f two pipes or w i th in a single pipe. In such a classif icat ion system, a booster p u m p discharging to a man i fo ld , two pumps connected to the same node or a surge tank found at the junc t ion o f even three p ipes can not be direct ly mode led . These restr ict ions do not exist in the classif ication or deve lopment o f boundary condi t ions discussed here. The ' i n - l i n e ' boundary condi t ion is s imp ly a special case o f a two -node e lement where both o f the nodes have first degree. In the def in i t ions above and the developments that fo l low, no restrictions are p laced o n the degree o f the nodes def in ing a boundary condi t ion. The s imultaneous classif icat ion o f nodes and external f lows in a network is a power fu l general tool wh i ch a l lows networks to be qu ick ly character ized accord ing to the d i f f i cu l ty o f pred ic t ing the i r transient response. The character izat ion system for the network i tse l f is presented i n the next sect ion and w i l l he lp to c lar i fy the importance o f the above def in i t ions. 46 4.3 C H A R A C T E R I Z A T I O N O F T H E N E T W O R K The benefits that result f rom the boundary cond i t ion c lassi f icat ion system are not immed ia te ly obv ious f r om the above def in i t ions alone. Cons idera t ion must also be given to the network as a whole, especial ly to the quest ion o f how boundary condit ions wi th a N D N o f two or greater are l i nked together. A s has been ment ioned earl ier, the real power o f the me thod o f characteristics is a consequence o f its expl ic i t nature. The expl ic i t nature o f the so lut ion procedure means that the unknown values at a given t ime step can be computed sequential ly wi thout hav ing to assemble and then solve a large system o f equations. Th i s p o i n t - b y - p o i n t so lut ion procedure has the enormous advantage o f produc ing an approx imate ly l inear re lat ionship between computat ion t ime and the size o f the network being solved. B y contrast, the solut ion o f a general system o f equat ions increases approx imate ly as the cube o f the number o f unknowns i n the p rob l em. F r o m this fact alone, even discount ing the addi t iona l compl icat ions o f increased computer storage and the d i f f i cu l ty o f obta in ing a re l iable so lut ion procedure for large non l inear systems, it is clear that expl ic i t procedures offer dist inct advantages over the s imultaneous equat ion approach when large networks are be ing studied. Howeve r , the situation when using the method o f characterist ics is sometimes not as s imp le as the above discussion impl ies . Wheneve r a boundary condi t ion has a noda l dependence number o f two or greater, the solut ion o f the boundary equations at one node depends on the condi t ions at one or more other nodes. I f al l the nodes def in ing this boundary cond i t ion have a complex i ty o f exactly one, no di f f icul ty is presented, and the boundary cond i ton can be so lved independent ly o f a l l the others in the network. Such a boundary cond i t ion is ca l led isolated. Howeve r , i f any o f the nodes i n quest ion are i n turn connected to other nodes or boundary condit ions, it is apparent that a s imultaneous solut ion o f a l l the unknowns w i l l become necessary. The interact ion o f a l l boundary elements shar ing this common set o f nodes must now be 47 accounted for and the so lut ion procedure becomes imp l i c i t Boundary condi t ions that interact w i th each other i n this manner w i l l be termed linked. T h e ideas presented above can best be i l lustrated by an example. F igu re 7 depicts a p u m p connected on both the suct ion and discharge sides to a larger network. It is assumed that the conduit connect ing nodes 2 and 3 is too short relat ive to the other p ipes i n the network to be economica l ly mode led by the method o f characterist ics and so is approx imated by a p ipe replacement e l emen t T h e suct ion side o f the p u m p is connected to an ord inary node (node 1) wh i le the discharge side is 3 connected to a complex C node (node 2). N o d e 2 is complex since three external connect ions — a pump, a pressure re l ie f va lve and a p ipe re lacement e lement — j o i n w i th i t Node s 1 through 3 may have any number o f other p ipes connected to them, as l ong as they are true pipes as def ined i n this thesis. T h e d i f f i cu l ty i n this example arises when one attempts to solve the boundary cond i t ion equations. Bo th the p u m p and the P R E are two -node boundary condit ions, so there exists a relat ionship between a l l three nodes i n the d iagram. Th i s means, in effect, that four f lows (the f low through the pump , the P R E and the pressure re l ie f va lve as we l l as the unknown external f low, Q) and three noda l head values must be (PRESSURE RELIEF VALVE) F igu re 7. L i n k e d Ne two rk Segment 48 obta ined s imultaneously. In order to write a general rout ine for solv ing such a p rob lem, it is necessary to develop the equat ions representing each element ind iv idua l l y and then assemble the result ing equations together to fo rm a complete system. Fo r the network shown in F i gu re 7 , a system o f seven or eight equations wou ld result (depending o n whether or not the pump ' s speed is changing). S ince the boundary condi t ions i n a hydrau l i c network are, in general , nonl inear, the result ing system must be l inear i zed and then solved by an iterat ion procedure. Usua l l y the in i t ia l values and the part ia l der ivat ives requ i red by the N e w t o n - R a p h s o n approach are avai lable, mak i ng this method o f so lut ion an attractive alternative. It is useful to prov ide a more quant i tat ive f ramework for the character izat ion o f the network than the above descr ipt ion al lows. Le t W represent a given set o f l i nked boundary conditons and let N be the set o f nodes that are l inked by the elements o f W . The number o f elements i n N w i l l be denoted by | N | and w i l l be ca l led the link length wh i l e the cardinal i ty o f W w i l l be wri t ten | W | and referred to as the group number o f a given l inked se t I f E j is used to represent the number o f equations i n the result ing system, then E j = | W | + | N | + P (1) in wh i ch P is the number o f two-parameter boundary condi tons. 1 3 App l i c a t i on o f this procedure to the example i n F igure 7 is straightforward. In F i gu re 7 , the set W contains the pump, the pressure re l ie f valve, the p ipe replacement e lement and the external f low Q . The elements o f the set N are s imp ly nodes 1 , 2 1 3 A two -pa rame te r boundary cond i t ion is an element whose state can not be complete ly descr ibed by the value o f f low through it and the head values at its end nodes. A n example o f such a boundary cond i t ion is a pump exper iencing power fa i lure. I n th is case, determinat ion o f the boundary condi t ion requires the speed o f the p u m p to be calculated, i n addi t ion to the rate o f f low through it Th i s boundary cond i t ion w i l l be discussed i n more detai l later. Ano the r example o f a two-parameter boundary cond i t ion is a control valve i n wh i ch the valve's pos i t ion and the discharge through i t both depend on the pressure head values present at its def in ing nodes. 49 and -3.^  Thus, for this l i nked set, the l ink length is three and the group number is four. The va lue o f P i n this network cou ld be either zero or one depending on whether or not the speed o f the pump is changing. Assuming the pump fai ls (hence P = 1), Equa t i on (1) results i n a E j value o f eight. A system o f eight equations w i l l have to be developed and solved i n order to predict the response o f this group o f nodes. Equa t i on (1) represents a s ignif icant contr ibut ion to the question o f how to characterize hydrau l i c networks. A s condit ions i n the network develop, ind iv idua l terms in the equat ion may change (e.g., i f a check valve closes, an external f low is e l iminated). Howeve r , at any t ime Equat ion (1) w i l l y i e ld the number o f equat ions that must be solved, wh i le the elements o f W and N w i l l ident i fy the components that make up this system. Since the work required to solve a system o f equations increases very rap id ly as the number o f unknowns in the system increases, the magni tude o f E j is c losely re lated to the d i f f i cu l ty o f so lv ing a g iven hydrau l i c system. In add i t ion, the manner in wh i ch E j changes as P R E ' s are substituted for short pipes in the network has s igni f icant impl icat ions for the selection o f an appropr iate t ime step. Th i s aspect w i l l be discussed i n more detai l in Chapter 5. T o the author 's knowledge, no other researcher has a l lowed for the existence o f l i nked boundary condit ions i n a hydrau l i c network . 1 4 The abi l i ty to inc lude l i nked sets in the network descr ipt ion greatly increases the range o f hydrau l ic systems that can be real ist ical ly studied. However , before the response o f the network can be predicted, the boundary equat ions themselves must be presented and solved. The so lut ion procedures and mathemat ica l representations o f the common hydrau l i c boundary condi t ions is the 1 4 T h e one poss ib le except ion to this may be the p rogram o f Pearce and Evans (1984). In this paper the authors state, w i th no further explanat ion, that they used a technique wh i ch requires the s imultaneous solut ion o f a l l the unknowns in the p rob lem, a l though they were using an exp l ic i t method. Even a l low ing for the existence o f l i nked boundary conditons, this procedure is very unusual and must be extremely expensive in terms o f computer t ime. 50 subject^ o f the remainder o f this chapter. 4.4 S I M P L E A N D O R D I N A R Y O N E - N O D E B O U N D A R Y C O N D I T I O N S O n e o f the most important boundary condi t ions i n the deve lopment o f a general network p rog ram is the m u l t i - p i p e fr ict ionless junct ion. It has been stated ear l ier i n this chapter that the d i f f i cu l ty o f so lv ing a network boundary cond i t ion is independent o f the numbe r o f p ipes connected to the node i n quest ion. The truth o f this statement w i l l now be demonstrated. The general C ^ node is so lved by a number o f investigators, such as F o x (1977), Va r dy (1976) and Koe l l e (1982), but none o f these writers seem to use the deve lopment to fu l l advantage. The benef i t o f the formulat ion presented here is that i t is used to obta in expl ic i t solut ions o f many o f the common boundary condit ions. Ra the r than beg inn ing wi th the s imple C ^ node, it is more convenient to der ive the general C ^ boundary cond i t ion and obta in the equations descr ib ing the s imple node as a special case. F i gu re 8 depicts the general s i tuat ion. Le t N 1 be the set o f a l l p ipes whose assumed di rect ion is toward the" node i n quest ion; s imi lar ly , let F i gu re 8. O rd ina ry N o d e 51 N2 be the set o f a l l p ipes whose assumed di rect ion is away f r om the junct ion. Le t one f low be ident i f ied as external and governed by an aux i l ia ry relat ion. The s ign convent ion adopted i n this thesis is that posit ive f lows are directed away f rom the junct ion be ing considered. The assumpt ion that the junc t ion is fr ict ionless is equivalent to saying that the hydrau l i c gradel ine elevat ion at the node can be represented by a single number , designated H p . F o r a l l p ipes be longing to the set N 1 , the C + characterist ic equat ion holds, whi le for members o f N 2 the C characterist ic equat ion appl ies. Equat ions (11) and (12) o f Chap te r 2 can then be rearranged to obta in H P C P i 1 1 and H C J J J i n wh ich Q p . ( Q p j ) respresents the discharge at the boundary section o f p ipe i ( j ) . 1 5 The cont inuity equat ion for the junc t ion requires the sum o f the f lows enter ing the node to equal the sum o f the f lows leav ing the node. Thus , the cont inui ty equat ion can be wr i t ten ? Q K - I Q P j - Q e x t = 0 (4) 1 3 J Equat ions (2) and (3) can be substituted direct ly into Equat ion (4) and an expression for H p obta ined. The result is H D = C - B Q t (5) P c c v ext v 1 5 I n this thesis, membersh ip i n a set w i l l be denoted by the symbo l e ( lower case epsi lon). F o r example, i n Equa t i on (2) ' i e N ^ means that the equat ion holds for a l l p ipes i be longing to the set N 1 . 52 in wh i ch B . B . c 1 1 D J (6) and C = c (7) Equa t i on (5) represents the general m u l t i - p i p e junc t ion with one external f low. The fo rm o f this equat ion is equivalent to a single C + characteristic equat ion and shows, when wr i t ten i n this manner, that any o n e - n o d e boundary cond i t ion located at a general network junct ion can be evaluated i n exactly the same manner as i f the boundary cond i t ion occurred at the downstream end o f a single p ipe. There is l ittle ind icat ion that the current attempts to solve the network p rob l em explo i t this advantage. Equa t i on (5) represents a single re lat ionsh ip between junct ion head ( H p ) and external f low ( Q g x t ) - Once a funct ional re lat ionsh ip representing a part icu lar hydrau l i c device is substituted into this expression, a single equat ion in a single unknown results. I f this re lat ionsh ip is either l inear or quadrat ic, an expl ic i t f o rmu la for the unknown can be obta ined. Several examples wi l l be presented i n order to i l lustrate this process. ident ica l ly zero. Th i s boundary cond i t ion can be used to mode l a l l s imple nodes i n a hydrau l i c network. B y inspect ion, Equa t i on (5) reduces to 4.4.1 S T M P L E N O D E S A N D L I N E A R R E S E R V O I R S I f the relat ion between Q g x t and H p is l inear, a number o f important results can be obta ined. The s implest example is the case when Q is H D = C . P c (8) 53 Ano the r elementary re lat ion occurs when Qtxl is e i ther constant or a known funct ion o f t ime. In this case, the value o f Q . can be substituted ^ e x t direct ly into Equat ion (5) to obta in the value o f the junc t i on head. Th i s boundary condi t ion can be used to represent constant d isp lacement pumps or other known f lows at any node i n the network. Reservo i rs or storage tanks may also be connected to ord inary nodes. The cont inu i ty equat ion for the reservoir relates the change in the water surface elevat ion to the rate o f discharge into the reservoir: thus, H = H + B ( Q + Q J (9) r o o w e v e x r v ' i n wh i ch H R is the hydrau l i c gradel ine elevat ion o f the surface o f the reservoir, H Q is the reservoir head at the beg inn ing o f the t ime step, Q G is the in i t ia l external discharge and B q is a constant The constant B q is equal to in wh ich A is the surface area o f the resevoir. No te that i f B is equal to zero, r o * ~ i » a reservoir head independent o f the external discharge is obta ined. H Q may be either constant or a known funct ion o f t ime (e.g. waves on the reservoir may be represented by a t ime vary ing H Q ) . I f the connect ion between the node and the reservoir is such that a l l hydrau l i c losses are negl ig ib le, the head at the node can be approx imated by the surface elevat ion o f the reservo i r . " Th i s a l lows equat ions (5) and (9) to be solved s imultaneously and results i n the re lat ionship C - H - B Q O = - £ ° ° e (10) c o Equa t i on (10) is an expl ic i t expression for the external f low into or out o f a storage reservoir and can be used to represent the behav ior o f devices such as " T h i s is equiva lent to saying that the var iat ion in head between the reservoir surface and the p ipe junc t ion is hydrostat ic. 54 constant head reservoirs, storage tanks and s imple surge tanks. These devices w i l l be referred to col lect ively as simple reservoirs. W h e n the cross-sect ional area o f the reservoir is smal l , neither the head losses due to fr ic t ion w i th in the tank, nor the inert ia o f the water co lumn i n the reservoir, shou ld be neglected. The representat ion o f the reservoir must now account for the mass o f the water i n the storage tank as it accelerates under the in f luence o f the head di f ference between the base o f the tank and the water surface, as we l l as the head loss due to f r ic t ion w i th in the reservoir. F igu re 9 i l lustrates the general s i tuat ion. A convenient l inear representat ion o f the reservoir is g iven by in wh i ch H ^ is the head at the base o f the reservoir and C j and are constants depend ing on the characterist ics o f the reservoir and the condit ions at H r , H 0 -HK.H' t ° " ' F i gu r e 9. A Reservo i r w i th H e a d Losses and Inert ia 55 the beg inn ing o f the t ime step. 1 7 The values o f C j and C ^ can be expressed as C . = H - H ' + 2 R Q | Q I 1 1 " 1 - C . Q (11a) 1 o r ^ e ' ^ e 1 2 v e v ' and c 2 = IA^T <LLB> i n wh i ch H ' is the head at the base o f the reservoir at the beginning o f the t ime interva l , L r is the in i t ia l length o f the water c o l umn and R r is the f r ic t ion loss coeff ic ient o f the reservoir. B y def in i t ion, H = z + L (11c) o r v ' i n wh i ch z is the elevat ion o f the base o f the storage tank above an arbitrary datum. The boundary cond i t ion represented by Equat ion (11) is usual ly referred to as a lumped inertia element. Equa t i on (9) (representing the storage component o f the reservoir) and Equa t i on (11) (representing the inter ia and head losses o f the reservoir water co lumn) can be comb ined into a single equat ion. Th i s re lat ion can be expressed as H , = C . + B , Q t (12) b b b ^ e x t v ' i n wh i ch C , = H + B Q + C , (12a) b o o ^ e 1 and 1 7 T h e der ivat ion o f this equat ion is g iven i n W y l i e and Streeter (1983), page 87. 56 The s igni f icant po int about Equa t i on (12) is that the reservoir re lat ion is st i l l l inear and, therefore, is no more d i f f i cu l t to solve then either o f Equat ions (9) or (11) taken ind iv idua l ly . Equa t i on (9) can be recovered f r om (12) by s imply setting the values o f and to zero. In addit ion, i f B Q is zero the inf in i te area reservoir boundary cond i t ion is recovered. Equat ion (12), therefore, represents a general storage element and w i l l be referred to as a linear reservoir i n this thesis. T he l inear reservoir inc ludes the effects o f f l u id storage, inert ia and fr ic t ion i n its most general fo rm. I f the head losses between the base o f the l inear reservoir and the node to wh i ch it is attached are very smal l , the head at the node H p can be equated to the reservoir h e a d ' H ^ . In this case, Equat ions (5) and (12) can be solved s imultaneous ly to produce C - C , c b (13) c b Once the value o f Q £ x t has been calculated, a l l other unknowns can be easily obta ined. 4.4.2 E X T E R N A ! , F L O W S W I T H E N E R G Y D I S S I P A T I O N W h e n the external f low f r om an ord inary node passes through an or i f ice or other restr ic i ton, a drop o f pressure head occurs in the f low direct ion and a new boundary condi t ion is obta ined. I f care is taken i n the formulat ion, a very general expression can be produced wh i ch can be used to represent many d i f ferent hydrau l i c devices. F o r example, valves and ori f ices d ischarging to l inear reservoirs, restr icted or i f ice surge tanks, o n e - w a y surge tanks and pressure re l ie f valves are a l l special cases o f this boundary condi t ion. T he abi l i ty to write a general expl ic i t re lat ion for a l l o f these elements connected to any ord inary node 57 overcomes the need to solve each e lement ind iv idua l ly . Th i s general treatment s impl i f ies both the preparat ion and mod i f i ca t ion o f input data. The general case o f an energy dissipating boundary cond i t ion can be deve loped by cons ider ing a va lve w i th t ime-dependent open ing discharging to a l inear reservoir. The s i tuat ion cons idered is identical to that i l lustrated in F igure 9, w i th the except ion that a va lve may now occur between the junc t ion and the base o f the reservoir. T he external f l ow may be related to the pressure head at the junc t ion by the expression Q e x t = s r E s [ s ( H p - H b J ] 0 - 5 (14) i n wh i ch r and E g are va lve parameters and s takes the s ign o f the external f l ow. 1 8 The constants E^ and E _ j are determined by know ing the values o f H p , H ^ and Q E X L for a g iven T for one case o f both posi t ive and negative f low. I f either E^ or E_^ is zero, check valve act ion is imp l i ed , since the discharge into or out o f the reservoir w i l l be prevented. The value ascr ibed to r is a convenient way o f represent ing an opening or c los ing valve. The value o f T is often a funct ion o f t ime but a constant va lue can be used to represent an or i f ice. Usua l l y a T va lue o f one imp l ies the va lve is i n a fu l ly open pos i t ion, wh i le a value o f zero a lways impl ies a c losed valve. Frequent ly , the hydrau l i c gradel ine elevation H ^ is constant, as i n the case o f a valve d ischarging to the atmosphere, or when an or i f i ce occurs at the entrance to a constant head reservoir . However , few compl icat ions are encountered " T h a t is, s = s g n ( Q ). The express ion means that s = l when the external discharge is away f r om the node and s = - l when the discharge is toward the junct ion. A n equiva lent expression is s = s g n ( H p - H ^ ) , s ince the discharge must a lways be f rom a h igher pressure locat ion to a lower one. 58 by assuming that is governed by the l inear reservoir equat ion . 1 9 The solut ion to the p rob l em o f an external f low enter ing a l inear reservoir through a control valve is obta ined by squar ing Equat ion (14) and substitut ing Equat ions (5) and (12) in to the result. T he final equat ion is a quadrat ic in QQXl and can be wri t ten as Q e x t + 2 b ( < e x t + c = 0 < 1 5> in wh i ch b = | ( r E s ) 2 ( B b + B c ) (15a) and c = s( T E g F ( C b - (15b) Cons iderat ion o f the f low di rect ion in Equat ion (14) leads to the fo l low ing solut ion o f Equat ion (15): Q e x t = ~ b + W ~ 3 ° ' 5 ( 1 6 ) in wh i ch s can now be determined as s g n ( C ^ - C c ) . Equa t i on (16) is the general fo rmu la for the discharge through an energy diss ipat ing contract ion into a l inear reservoir f r om an ord inary network node. Equa t i on (16) is extremely general since the reservoir may represent the atmosphere, a constant head storage tank, a s imple reservoir or a l inear reservoir. S imi la r ly , the head loss may occur through a contro l valve or through an or i f ice, and check va lve act ion can be inc luded i f desired. A n y o f these opt ions can be selected by s imp ly choos ing values o f the " U s i n g the l inear reservoir equat ion invo lves no addi t ional l imi tat ions since, as has already been po inted out, the equat ion inc ludes both the s imple and constant head reservoirs as a special case. 59 constants def in ing Equat ions (15a) and (15b). 4.4.3 ATR C H A M B E R S Ano the r impor tant boundary cond i t ion is that o f an air chamber connected to an ord inary network node. A i r chambers are useful surge supression devices: they attenuate transient pressures by acting as temporary storage elements, exchanging water w i th the junc t ion p ipes as the pressure head at the node changes. A l t hough air chambers are discussed i n a numbe r o f standard references, the treatment is not as general as that presented here. The general arrangement o f the air chamber is i l lustrated i n F igure 10. In this d iagram a short connect ing p ipe jo ins an ord inary network node to an a i r chamber. The external f low passes through an or i f i ce or f low restr ict ion as it enters the chamber and the water surface i n the chamber is located at a distance o f z above the datum. Above the f lu id interface, a vo lume V- o f a ir or F i gu re 10. A n A i r Chambe r 60 other gases is present i n the tank. It is assumed that the connect ing p ipe can be represented as a l umped inert ia e lement and that the behav ior o f the gas is governed by a po ly t rop ic gas law expression. The gas law for the chamber can be written H V - r - C = 0 (17) a ac ' i n wh i ch H & is the absolute pressure head o f the gas in the tank, r is the po ly t rop ic exponent and C is a constant depending on the in i t ia l gas vo lume and the pressure present i n the chamber. The exponent r is equal to 1.0 for i sothermal processes and 1.4 for reversible adiabat ic processes. The actual response o f the air chamber depends on the rate o f heat transfer between the gas in the chamber and the surrounding env i ronment (Graze , 1968), but it is se ldom possib le to calculate the heat exchange in detai l . The usual procedure for most purposes is to assume a value for r o f 1.2 (Chaudhry , 1979; Wy l i e and Streeter, 1983). Account ing for the change i n air vo lume in the chamber due to the external f low, the gas law equat ion can be rewritten for the end o f the current t ime step. The result ing equat ion is < H P 2 + H b a r - ~ K>ext + " C a c = 0 < 1 8> in wh i ch H p i n the unknown hydrau l i c gradel ine elevat ion at the end o f the t ime step, is the barametr ic pressure head and z is the water surface elevat ion at the end o f the t ime step. Equa t i on (18) must be expressed i n terms o f a single unknown i f i t is to be eff ic ient ly solved. The unknowns H D and H D in F igure 10 can be related to the pressure head at the node by comb in ing the lumped inert ia e lement wi th the or i f i ce loss equat ion. Rewr i t i ng equat ion (14) i n terms o f these head values results i n 61 H p 2 = H p i - ^ Q ^ (19) s i n wh i ch a l l terms fo l low the i r prev ious def in i t ions. The heads H p and H p are re lated by the l umped inert ia re lat ion. Equat ion (11) is now wr i t ten as H p - H P i = C j + C 2 Q e x t (20) i n wh i ch the constants and C 2 can be calculated f r om extending Equat ions (11a) and ( l i b ) to fit the current situation. A single re lat ion between H D and Q £ x t can be obtained by subst itut ing Equat ion (5) into Equa t i on (20) and the result into Equat ion (19). The final expression for H p is H p j = C c - C l - ( B c + C l ) Q e x t - ( 7 ^ y 2 Q e 2 x t . ( 2 D Since, i n general, the water surface level i n the chamber changes as water enters or leaves it, z is also a funct ion o f the external f low. Usua l l y the s imple reservoir re lat ion is used: that is, i n wh i ch Z q is the water surface he ight at the beg inn ing o f the time step and B is a constant depend ing o f the d imensions o f the a i r chamber. Equat ions 3.C (21) and (22) can now be substituted into Equat ion (18) to produce a single n on - l i n e a r equat ion i n a single unknown ( Q ). The value o f Q can be C A L C A I obta ined f r om this final equat ion using Newton ' s method and the known values o f the in i t ia l condit ions. L i k e the case o f an external flow w i th energy dissipat ion, the air chamber boundary cond i t ion is ent irely general. T he network node to wh ich it is connected may have arbi trary degree, the connect ing p ipe may or may not be present, the or i f ice loss may be unequal for posit ive and negative flows and the value o f this loss may be t ime dependent Sett ing the value o f 62 . E to zero for either -positive or negative f low al lows a check valve to be inc luded. A n d , finally, the value o f r can be adjusted to produce the desired gas law express ion. 2 0 Desp i te its general ity, this fo rmulat ion o f the air chamber boundary cond i t ion is no more d i f f i cu l t to solve than the more restrict ive developments usual ly presented in the l i terature. 4.4.4 P O W E R F A I L U R E O E A P U M P W I T H I N E R T I A The last o ne - node boundary cond i t ion to be cons idered i n detai l w i l l be a pump subjected to power fai lure. The condi t ions engendered by loss o f power to a pump ing un i t represent the most compl icated o f the one -node boundary condit ions. The calculat ion o f the pump ' s response to power fai lure must account for three dist inct aspects o f the system's per formance: the characteristics o f the p u m p itself, the rotat ional inert ia and decelerat ion o f the pump motor and impe l l e r and, finally, the head discharge re lat ion on both the suct ion and discharge sides o f the pump . F o r the first t ime, the analysis o f a c | node must introduce a second var iable — the rotat ional speed o f the pump — into the boundary cond i t ion so lut ion. The analysis o f the power fa i lure cond i t ion is c lear ly deve loped i n a number o f references, such as Chaudh r y (1979) and Wy l i e and Streeter (1983). F o r this reason, the presentat ion given here w i l l be very brief. Deta i l s w i l l be p rov ided on ly where they a l low the t rad i t iona l treatments to be extended or general ized to cover a broader range o f appl icat ions. F o r example, the on ly 2 0 A n interest ing result is obta ined i f r is equal to zero. In this case, a constant pressure process is represented and the air chamber boundary cond i t ion reduces to the 'external f low wi th energy d iss ipat ion ' boundary condi t ion. If, i n addi t ion, the head losses at the entrance to the tank are smal l , then a l l the boundary condit ions relat ing to l inear reservoirs are also recovered. Thus , the a i r chamber boundary cond i t ion inc ludes the behav ior o f a l l other o f the boundary condi t ions discussed up to this p o i n t Th i s fact seems to be o f more academic than pract ical importance. 63 changes i n the pump equations appear i n the head balance expressions, and it is therefore not necessary to make any alterations to the inert ia balance relations. Thus , the to rque- ine r t i a equat ion w i l l n o t - be presented and the reader is re ferred to the standard references re lat ing to hydrau l ic transients. W h e n consider ing pumps undergo ing power fa i lure it is convenient to represent the characteristics o f the p u m p in dimensionless fo rm. The p r imary var iables o f head, discharge and rotat ional speed are made dimensionless by using the pump ' s per formance at the rated cond i t i on 2 1 as the reference quantit ies. Thus , the var iables become Q E X T K H D „ N « T , v = T V h = TTR' " = TTR- = <23a> i n wh i ch the subscript R refers to the var iable at the rated condi t ion, N is the rotat ional speed o f the pump, represents the dynamic head o f the pump, and T is the torque acting on the p u m p axle. I f the p u m p is assumed to obey the homo logous ratios der ived f r om d imens iona l analysis, then the p u m p characterist ics can be represented i n a s imple manner. In part icular, the d imensionless pressure head can be expressed as a funct ion o f the d imensionless speed o f the pump and the discharge though i t F o r computer solut ion, a usefu l f o rm o f this re lat ion is the p lo t o f — ^ — 2 vs. t a n - 1 ^ ) (23b) a1 + v The re lat ive torque can be plotted i n an analogous fashion. It is usual to represent the pump ' s characteristics by a series o f equal ly spaced t a n 1 ^ ) values and then to find the head values at intermediate points by l inear or parabo l i c in terpo lat ion. I f the l inear approach is adopted, the dynamic head can be " T h e rated cond i t ion o f a p u m p is its po in t o f peak or m a x i m u m eff ic iency. 64 expressed as H d = H R ( c 3 + v ' ) [ a i + . a j t a n - 1 ^ ) ] (23c) i n wh i ch aj and a 2 are known interpolat ion constants. The general head balance equat ion can be der ived by cons ider ing a pump undergo ing power fa i lure d ischarging to a l inear reservoir. The head balance re lat ionship can then be wri t ten as H b = H p ± H d - H y (24) i n wh i ch H ^ is the reservoir head and H y represents the head loss across any contro l valves present at the pump. The s ign o f the dynamic head term is used to dist inguish between a p u m p wh ich norma l l y discharges f r om the network (plus), f r om a pump wh i ch discharges to the network (minus). It is necessary to wr i te Equat ion (24) i n this manner i n order to be consistent wi th the sign convent ion o f Q , a l though the usual arrangement o f pumps is for them to CX.L discharge to the network. The loss across the p u m p valves can be wr i t ten H = ,—4rr 2Q 1 , (25) v (T E ) 2 V : e x t v ' s i n wh i ch a l l terms have been prev ious ly def ined. Th i s representat ion is useful since it a l lows check va lve act ion as we l l as opening or c los ing valves to be represented. A single equat ion invo lv ing v and o can now be obtained. U s i ng Equa t i on (5) to relate H p to QQXt and the l inear reservoir re lat ion (Equat ion (12)) to relate H ^ to Q g x t produces a funct ion ( F ) wh ich can be written F = C c " C b " QR<Bb + B c > v ± H d " QRVE7v2 = 0 ™ Equa t i on (26) can be used i n conjunct ion wi th Equat ion (23c) to obtain a single 65 n o n - l inear equat ion w i th unknowns v and c. Th i s equat ion can then be solved s imultaneous ly w i th the rotat ional inert ia! equat ion to complete ly def ine the pump fai lure boundary condi t ion. Aga i n , the N e w t o n - Raphson technique is a convenient way o f so lv ing this pa i r o f non - l i n ea r equations s imultaneously. 4.4.5 O T H E R B O U N D A R Y C O N D I T I O N S ' Th i s thesis has emphas ized a general treatment o f one -node boundary condi t ions connected to ord inary network nodes. T h e examples presented have been chosen to i l lustrate techniques wh i ch are more broadly appl icable than the usual approaches g iven i n the l i terature. Nume rou s other one -node boundary condi t ions exist, but these w i l l not be presented i n detai l . These boundary condi t ions are omi t ted since they can be so lved i n an exactly analogous manner to their s imple p ipe l ine counterparts. Tha t is, i t has been shown how the behav ior o f a general network node can be descr ibed by an equat ion wh ich is ident ical in f o rm to a C + equat ion. Thus , the so lut ion o f many o f the omit ted boundary condi t ions can be obta ined by s imp ly rep lac ing the usual C + o r C characterist ic equations by Equa t i on (5). Dev i ces such as pumps operat ing in the steady-state, a i r - i n l e t valves, l umped capacitance elements and sp r i ng -mass -dashpo t systems can a l l be analyzed i n this way. O f course, devices exist wh ich can not be def ined by the condi t ions occurr ing at a single node. These two - node boundary condi t ions w i l l be discussed i n the next section. 4.5 O R D I N A R Y T W O - N O D E B O U N D A R Y C O N D I T I O N S The hydrau l i c e lements discussed up to this po in t have requ i red the pressure head to be determined at on ly one network node. However , when the solut ion o f a boundary cond i t ion depends on the behav ior o f two nodes, a mod i f i ed procedure is ca l led for. E lements that demonstrate this two - node dependence are used i n essentially four 66 di f ferent ways: to supply energy to the f low, for energy diss ipat ion, for automatic f low control and, f ina l ly , as p ipe replacement elements. The standard device for increasing the energy i n the downstream direct ion is a booster pump; elements such as valves, or i f ices and other local losses dissipate f low energy, usual ly with the intent o f ach iev ing some measure o f passive f low contro l . Pressure regulat ing and pressure sustaining valves are examples o f automat ic f low control devices. These valves open and close, not accord ing to some predetermined pattern, but i n response to the changing condi t ions i n the network. The purpose o f this contro l is usual ly to mainta in the pressure head o f a g iven node at a more un i f o rm level than wou l d otherwise be possible. F ina l l y , as has been ment ioned ear l ier i n this thesis, p ipe replacement elements can be an economica l way o f represent ing the very short conduits present in the network, thus a l low ing a larger t ime step to be used in the analysis. A l t hough the boundary condi t ions exh ib i t ing '^two-node behav ior are quite diverse, a general f ramework can st i l l be developed. Usua l l y the behav ior o f a two -node boundary cond i t ion depends p r imar i l y on the di f ference in pressure head between the upstream and downstream nodes, and not on the actual head present at either node. F i gu r e 11 shows a general two - node boundary cond i t ion where, for convenience, the upstream node is ident i f ied by the subscr ipt 1 and the downstream TWO - NODE ELEMENT F igu re 11. Gene ra l T w o - N o d e Boundary Cond i t i on 67 node by the subscr ipt 2. Node s 1 and 2 may be o f arbitrary degree but must both have a complex i ty o f one. Equat ion (5) can be app l ied to each node i n turn, bear ing i n m i n d the change i n the sign o f the external f low between the upstream and downstream nodes. Th i s results i n H D = C - B Q t (27a) P i Ci Ci ^ ext, for node 1, and H P 2 = C c 2 + B c 2 < < e x t 2 ( 2 7 b ) for node 2. W h e n consider ing two -node boundary condit ions, i t is usual to assume that the f low between the upstream and downstream nodes is incompress ib le and can be represented by a single discharge value, Q e x t - I f this assumpt ion is made, Equa t i on (27b) can be substracted f r om (27a) to produce H p ] - = 6 - B Q e x t (28) A A i n wh i ch C = C - C and B = B + B . Equat ion (28) represents the Ci C2 Cj c 2 general expression for the head dif ference between the upstream and the downstream nodes o f a two - node boundary condi t ion. I f the boundary cond i t ion i tsel f prov ides a second re lat ion between this pressure head di f ference and the external f low, the value o f Q can be calculated. Th i s process w i l l be i l lustrated w i th a number o f c ommon two - node boundary condit ions. 4.5.1 B O O S T E R P U M P S U N D E R G O I N G P O W E R F A I L U R E Booster pumps are frequent ly used devices for increas ing the pressure head between two nodes. U n d e r steady-state condit ions, the increase i n head between the downstream and upstream flanges o f a p u m p is a funct ion o f the f low through the pump, the re lat ion often be ing expressed as a piecewise l inear 68 or quadrat ic po lynomia l over the operat ing range o f in terest W h e n the boundary cond i t ion is combined wi th Equa t i on (28), the solut ion fo l lows immediate ly . The more interest ing example o f the pump boundary condi t ion occurs when the booster p ump undergoes power fai lure. In this case, the head balance equat ion for the pump can be wri t ten F = H p i - H p j + H d - H y = 0 (29) i n wh i ch is again the dynamic head produced by the p u m p and H y is the head loss across any p u m p contro l valves that may be p resent I f the pump ' s characterist ics are represented in the same manner as prev ious ly developed, Equa t i on (23c) can be substituted for in this re lat ion. In add i t ion, Equat ion (25) can be used to represent the head loss across the p u m p valves and Equa t i on (28) can be substituted direct ly into Equat ion (29). The result is F = C - B Q v + H (o» + v ' ) [ a i + a 2 t a n - ' ( ° ) ] - _ * Q ^ v ' (30) v s ' T h i s equat ion replaces Equa t i on (26) in booster pump appl icat ions and can be so lved i n a s imi l iar fash ion to the corresponding one -node power fa i lure boundary condi t ion. 4.5.2 T W O - N O D E B O U N D A R Y C O N D I T I O N W I T H E N E R G Y D I S S I P A T I O N W h e n the f low between two nodes passes through an or i f ice or other flow restr ict ion, part o f the flow energy w i l l be dissipated as heat and sound. T h e representat ion o f this loss can be obta ined by a s l ight mod i f i ca t ion o f Equa t i on (14) to produce Q e x t = s ' E s r s ( H p f h P / ' 5 ( 3 1 ) i n wh i ch s = s g ^ H p - H p ). Equat ion (31) can be squared and Equat ion 69 (-28) substituted in to the resu l t Th i s produces i n wh ich b = | ( r E g ) 2 B (32a) and c = s( r E g ) 2 C (32b) Cons iderat ion o f the direct ion o f f low a l lows the sign o f the radical to be correct ly chosen when so lv ing this quadrat ic. Th i s results in an expl ic i t re lat ion for the external f low: Q e x t = " b + s [ b 2 + c ] ° ' 5 ( 3 3 ) ' i n wh i ch s now takes the sign o f C . A s was the case w i th the one -node boundary cond i t ion, a wide range o f loss elements can be accomodated by this formulat ion. 4.5.3 P R E S S U R E R E D U C I N G V A L V E S The pressure reduc ing valve ( P R V ) is an automatic f low control device wh i ch is f requent ly employed i n water supp ly networks. The purpose o f the P R V is to adjust the pressure head losses between an upstream and downstream node in order to ma inta in the hydrau l i c gradel ine elevat ion o f the downstream node at a predetermined leve l . A P R V that is ma inta in ing this predetermined head is said to be operating and the pressure head setting o f the valve is termed its operating head (designated H ). It is assumed that a P R V acts as a check valve and thus prevents reverse f low when the head at the downstream node exceeds 70 the head at the upstream node. The head loss across the pressure reduc ing valve w i l l be expressed as H 0 - H D = K Q 2 , (34) ?i P2 s v e x t v ' i n wh i ch K g is a loss coeff ic ient (wh ich , i n general, is not constant) and the subscr ipt s takes the sign o f the external f low. The assumpt ion o f check valve act ion means that K equals zero when H D is greater than H p . It is also S r 2 i i assumed that even a fu l ly open P R V w i l l have a smal l loss associated wi th it Th i s m i n i m u m value o f the loss coeff ic ient w i l l be designated K . Thus, i f the upstream head exceeds the head at the downstream node but is less than the operat ing head o f the valve, this m i n i m u m head loss coeff ic ient w i l l apply. The solut ion o f this pressure reduc ing valve boundary cond i t ion can now be developed. Since the P R V is essential ly an energy diss ipat ing device, the arguments i n the prev ious section apply. Th i s means that the direct ion o f the A A external f low is st i l l determined by the sign o f C . Thus, i f C is negative, check valve act ion is imp l i ed — the external flow w i l l be zero and both nodes 1 and 2 w i l l be s imple. However , i f C is posi t ive, two poss ib i l i t ies exist: either the P R V is operat ing, i n wh i ch case H p w i l l equal H , o r the P R V is fu l ly open, r 2 S i n wh i ch case K g w i l l equal K . In i t ia l ly , the assumpt ion can be made that the pressure reducing valve is operat ing. The value o f Q£Xl can then be obta ined f r om Equat ion (27a) wh i ch , i n turn, a l lows the head at the upstream node to be obtained f r om Equat ion (27b). F ina l l y , since H D , H p and Q are i i 12 ext a l l known , the value o f K g i n Equa t i on (34) can now be calculated. I f this computed value o f K g exceeds the m i n i m u m value, then the P R V is indeed operat ing and the solut ion is complete. However , i f K g is less than K m > the value o f K must be used instead o f K i n Equat ion (34). Wheneve r the value m s o f K g is known, Equat ions (28) and (34) can be equated to produce a quadrat ic 71 equat ion i n Q e x t - The result ing equat ion is K Q 2 f + B Q - C = 0 (35) m ext ^ e x t wh i ch can be easily solved. The determined value o f Qtxt can be used to compute the pressure head values at the end nodes using Equat ions (27a) and (27b). The so lut ion o f other boundary condit ions representing automatic control valves fo l lows a s imi la r pattern to that i l lustrated above for the pressure reducing valve. Usua l l y the behav ior o f these valves introduces no new equations into the solut ion procedure, a l though the so lut ion logic can be compl icated to some extent O f ten the state o f the valve w i l l have to be assumed in advance o f the computat ion, and the correctness o f this assumption either con f i rmed or contradicted by the subsequent analysis. I f a reasonable mode l can be developed to describe the valve 's performance, the solut ion procedure appl icab le to automatic valves can frequent ly be made exp l i c i t 4.5.4 P I P E R E P L A C E M E N T E L E M E N T S Large hydrau l i c networks w i l l f requent ly contain conduits hav ing a wide var iety o f lengths and wavespeeds. S ince the method o f characterist ics requires each p ipe to have at least one reach, the t ime step needed to accurately represent the shortest 2 2 o f these p ipes may be uneconomica l ly smal l . It is 2 2 B y 'short ' i t is meant p ipes in wh ich the rat io o f length to wavespeed is smal l re lat ive to the t ime step used i n the analysis. That is, i f even a single characteristic reach is imposed on the p ipe, the wavespeed imp l i ed by the length o f this reach w i l l be much greater than the actual wavespeed o f the condu i t Usua l l y some adjustment i n the wavespeed is a l lowable , since its va lue is se ldom known wi th great accuracy. Howeve r , most researchers seem to feel that adjustments o f more than ten o f f i f teen percent are not jus t i f i ed (e.g., Wy l i e and Streeter, 1983). It is interest ing to note, however , that K oe l l e (1982) suggests the analyst s imply impose at least one reach on each p ipe and accept the penalty o f the imposed wavespeed change. Koe l l e calls such a one reach p ipe ah inertia. 72 desirable, therefore, to replace these short p ipes by approx imate representations and thereby use a larger t ime step than wou ld otherwise be possible. Such a boundary cond i t ion is cal led a pipe replacement element i n this thesis and it is the last two - node boundary cond i t ion that wi l l be discussed. Two useful techniques w i l l be descr ibed for mode l i ng these short conduits: the first uses the lumped inert ia e lement discussed ear l ier and the second used a finite difference approx imat ion to the transient govern ing equations. One o f the most convenient p ipe replacement elements is the lumped inert ia element that was first int roduced a long wi th the l inear reservoir boundary condi t ion. In a l umped inert ia element, the f low between nodes 1 and 2 is assumed to be incompress ib le and can thus be represented by a single variable, Q e x t - I f L is the length o f the p ipe to be replaced, A its cross-sect ional area and R its resistance coeff ic ient, then Equat ions (11), (11a) and ( l i b ) can be rewri t ten to apply to the present s ituat ion. The result is H p ] - = C j + B j Q e x t (36) in wh i ch C j = H j - H i + 2 R Q e | Q e | n ~ 1 - B j Q e (36a) and B i = g-n-t • <36b> In these expressions, Q e > H , and H 2 represent the discharge va lue and the hydrau l i c gradel ine elevations at the beg inn ing o f the t ime step, wh i le Q D V f , H D ext r i and H p represent these same quant i t ies at the end o f the current t ime step. Equa t i on (36) can be substituted into Equa t i on (28) and the result ing expression solved for Q g x t - T he final equat ion is 73 O = ^ - C W e x t , • (37) B + Bj i n wh i ch a l l terms on the right hand side are known constants. The major weakness o f the l umped inert ia mode l is the assumpt ion o f incompress ib le f low. Th i s assumpt ion means that disturbances travel the length o f the p ipe instantaneously and that no f lu id storage is a l lowed in the p ipe. The second technique for p ipe replacement that w i l l be cons idered overcomes both o f these l imitat ions, but results i n a somewhat more compl icated so lut ion procedure. It has been prev ious ly ment ioned that f in ite d i f ference techniques are general ly not a pract ical method o f obta in ing transient solut ions to large hydrau l i c networks. However , smal l segments o f a network can somet imes be effect ively replaced by f in ite di f ference approx imat ions and treated as special boundary condi t ions w i th in the general f ramework o f a method o f characterist ics so lut ion (Chaudhry , 1979). W h e n this hyb r id approach is app l ied to the short conduits i n a network, a very power fu l p ipe replacement e lement is produced. Th i s new P R E accurately represents the storage and t ime delay characteristics o f the real p ipe, as we l l as the inert ia and fr ict ional resistance o f the f lu id i n the p ipe. In order to app ly the f in ite d i f ference approach to a p ipe l y ing between nodes 1 and 2 in F igu re 9, f in ite d i f ference quotients must be wr i t ten for the part ia l der ivat ive terms i n the govern ing equat ions . " F o r example, to replace the H t te rm i n the cont inui ty expression the fo l l ow ing quot ient can be used: H p + H p - H i - H 2 H = — ^ ^ (38) t 2 A t O the r fract ions representing the H x > Q t and Q x can be s imi la r ly writ ten. These f in i te d i f ference approx imat ions can then be substituted in to the govern ing " T h e representat ion o f the part ia l derivat ives presented here is g iven i n Streeter's 1972 paper, a l though a few m ino r changes have been made. t 74 equations (Equat ions (3) and (4) o f Chap te r 2) to produce a lgebraic expressions for a l l o f the unknowns i n the p rob lem. T h e one d i f f i cu l ty i n this process l ies i n the f r i c t ion term. The general expression for the head loss due to fr ict ion can be expressed as H f = RQ|Q| n _ 1 (39) Usua l l y , the discharge in this expression is evaluated as Q = 0.25(Q e ] + Q e j + Q e x t ] + Q ^ ) (40) but i f Equa t i on (40) is substituted in to Equa t i on (39), the result ing equat ion becomes non l inear (Streeter, 1972). In order to achieve an exp l ic i t so lut ion to this boundary condi t ion, the fo l l ow ing l inear izat ion o f the f r ic t ion term is suggested. Le t RQIQI""1 = 0.25(Qei + Qe2 + + Q ^ O S O ^ + Q j (41a) or, expressed more convenient ly, R Q | Q | n _ 1 = j | - x ( Q e i + Q e 2 + Q e x t ] + Q e x t 2 ) (41b) i n wh i ch Th i s l inear izat ion is s l ight ly more accurate than that proposed by Chaudh r y (1979). W h e n Equat ion (41b) is used to represent the head loss due to f r ic t ion i n a p ipe replacement element, and the appropr iate f in ite d i f ference quotients are substituted in to the govern ing equations, a set o f two l inear a lgebraic equations is produced. 75 T h e so lut ion o f the algebraic finite di f ference equations represent ing the cont inu i ty and dynamic equations is easily obtained. Equat ions (27a) and (27b) can be substituted into the f inite difference equations in order to e l iminate the head terms. The result ing two expressions can then be solved s imultaneously to calculate the external f low at either the upstream or downstream nodes. The f ina l result is Q = ^ + ^ (42) ^ex t ] E ] E S + E2E4 v ' i n wh i ch E i = B + R + (43a) Ci g A A t v • Ej = B + R + - A - £ r , (43b) ^ c 2 g A A t v ' E > = c C l " c c 2 + H - H ^ + (ini - R X Q e i + V - ( 4 3 0 E , = B + - ^ M - <43d) Ci g A A x v ' a* A t ' c 2 ' g A A x ' E s = B „ + (43e) Once Q has been found, al l other unknowns can be easily determined by ext i back subst itut ion. Equat ions (42a) to ( 420 are considerably more d i f f i cu l t to evaluate than the constants descr ib ing the l umped inert ia element, but they also represent a more accurate representation o f the phys ica l behav ior o f the c ondu i t T h e quest ion o f how both o f these elements can be effect ively employed i n the context o f a general network program wi l l be discussed in Chapter 5. 76 4.6 B O U N D A R Y C O N D I T I O N E X T E N S I O N S T h e basic elements and boundary condi t ions that compr ise a general hydrau l i c network have now been formulated. The solut ion to the equations descr ib ing the response o f the var ious hydrau l i c devices has been put into the context o f a general boundary cond i t ion classif icat ion system. The f ina l remain ing di f f icu l ty is to extend the current results to inc lude more compl icated arrangements and situations than those that have been considered so far. F o u r speci f ic extensions w i l l be considered: boundary condit ions invo lv ing more than two nodes, internal boundary condit ions wh i ch occur w i th in a single p ipe, the so lut ion to l i nked groups o f boundary condit ions and, finally, a proposa l for inc lud ing junc t ion losses at a general network node. The first three o f these extensions w i l l on ly be br ie f ly discussed here; the final extension is more compl icated and w i l l be cons idered i n more detai l . Boundary condi t ions wh i ch have a noda l dependence number o f three or greater are sti l l rare i n hydrau l i c networks. However , as automatic control systems become more c ommon and complex, i t is l i ke ly that their use wi l l increase. A typical example o f such a mul t ip le node boundary cond i t ion is a surge ant ic ipator valve. In these valves, a sensing mechan ism can be p laced at a locat ion remote f r om the valve itself. W h e n a prescr ibed state is sensed at this locat ion (e.g., closure o f a check valve on the discharge side o f a pump) , a preset va lve act ion is in it iated. The key characteristic o f such boundary condit ions is that their behav ior is determined by the condit ions at three or more nodes, but that the external f low is usual ly between two nodes only. Th i s means that the so lut ion o f ord inary three node boundary condit ions can be computed i n an analgous manner to the boundary condit ions already presented, w i th on ly m ino r compl icat ions i n the contro l logic. The procedure i l lustrated i n the descr ipt ion o f the pressure reduc ing va lve boundary condi t ion remains typical. The cr i t ical step i n so lv ing such a boundary cond i t ion is to develop a fo rma l mode l wh ich real ist ical ly describes the behav io r and response o f the hydrau l i c device be ing 77 considered. Ano the r compl icat ion o f the so lut ion procedure arises when m i d - p i p e boundary condi t ions are considered. U p to this point, it has been assumed that the head and discharge at a l l internal p ipe sections can be computed by stra ightforward appl icat ion o f Equa t i on (13) o f Chapter 2. However , i f an air pocket occurs w i t h i n a p ipe, or i f water c o l umn separation is a l l owed to take place, the s i tuat ion becomes more invo lved. In this case, an unique head and discharge must be assigned to both the upstream and downstream ends o f each reach i n the pipe, effect ively doub l ing both the storage and computat iona l time requ i red to solve a given hydrau l i c network. T o avo id these prob lems, a reasonable approach is to assume that internal boundary condit ions do not occur except i n special ly designated pipes. These special p ipes can be ana lyzed separately and the rest o f the network solved by the more standard techniques. The actual so lut ion procedure appl icab le to internal boundary .condit ions is we l l developed i n the standard reference l i terature and need not be cons idered here. A general network analysis tool shou ld be able to solve, not on ly isolated boundary condit ions, but l i nked groups o f boundary condi t ions as we l l . A s long as every node in a network has a complex i ty o f at most one, there is no poss ib i l i ty that boundary condi t ions can become l i nked , and each e lement can be so lved independent ly o f the others. In such cases, the analysis procedures presented above are sufficient, and no extensions are necessary. Howeve r as soon as more than one external flow connects to a single node, the junc t ion cont inui ty relat ion is a l tered and the final expressions for the external f low presented above no longer apply . T h e solut ion procedure must now work w i th more fundamenta l boundary cond i t ion equat ions and solve for each unknown i n the l inked set s imultaneous ly. The ma in di f ference i n this procedure is that when Equa t i on (5) is app l i ed to a complex node, the external f low is no longer a single value, but the sum o f a l l external flows that j o i n w i th the node. T h e technique for obta in ing the so lut ion must be able to recognize l inked boundary 78 cond i t ion elements and to correct ly assemble the ind iv idua l equations into to a total system. However , once this assembly is complete, standard equat ion so lv ing procedures can be app l ied i n order to determine the so lut ion to the total l i nked se t Because the ma i n d i f f icu l t ies is solv ing complex networks relate to the recogni t ion and assembly o f the proper equations, and not to the fo rm o f the equations themselves, further cons iderat ion o f this p rob lem is deferred to Chapte r 5. The f ina l subject to be considered i n this chapter is the p rob l em o f solv ing an arbi trary network junct ion when the head losses due to fr ic t ion at the junct ion are signif icant. Since the topic o f head losses at a mu l t i p l e - p i pe junc t ion is very-compl icated, this subject w i l l be considered separately. 4.6.1 GENERAL JUNCTION EQUATION WITH LOSSES T o the author's knowledge, no one has attempted a so lut ion o f a general j unc t i on boundary cond i t ion in the unsteady state when head losses due to f r ic t ion at the junct ion are not assumed to be negl igible. Th i s is not surpris ing, since a rev iew o f the avai lable l i terature reveals that much confus ion st i l l exists regard ing the losses at mu l t i p l e - p i p e junct ions under steady-state condit ions. F o r even a three p ipe junct ion, the head losses (and the f low distr ibut ion) has been shown to be a funct ion o f the Reyno lds number o f the f low i n each pipe, the di rect ion o f f low in each p ipe , ' 4 the geometry o f the p ipe junc t ion as wel l as the detai ls o f the junct ion construct ion. M i n o r factors, such as the presence o f mit re bends or whether i rregular i t ies have been left on the surface o f the p ipe wal ls, can have a surprs ingly large inf luence o f the head loss values. One o f the basic prob lems the analyst encounters is that the under ly ing assumption o f " T h a t is, d i f ferent values o f head loss are obta ined depending on whether the f low is converg ing (two f lows i n and one out) or d iverg ing (one f low i n and two out) at the node. 79 one -d imens i ona l flow breaks down at such a boundary condi t ion. Deta i l ed analysis wou l d require deve lopment a th ree -d imens iona l turbulence mode l wh ich accounts for viscous diss ipat ion, flow separation and other effects both in the junc t i on and in the adjacent pipes. Obv ious ly , such a complex mode l is not jus t i f i ed . Thus the analyst reaches an impass: impos ing a one -d imens iona l mode l o f the junc t ion flow is bound to give mis lead ing results under some circumstances, wh i le use o f a more real ist ic th ree -d imens iona l mode l is expensive, t ime - con sum ing and, perhaps, no more accurate than its one -d imens i ona l counte rpar t It is for al l o f the above reasons that junc t ion losses have general ly been ignored i n network appl icat ions and, i n many situations, such neglect is we l l jus t i f ied . However , dur ing some f low condit ions, h igh veloc i ty flows occur at junct ions o f cr it ical importance. A typical example occurs when an urgent ly needed f ire flow is pumped f rom a water d is t r ibut ion network. M a n y fai lures have occured as a result o f a greater f low be ing demanded o f the system than it is capable o f supp ly ing. Because o f the importance o f these h igh velocity f lows, the assumpt ion that junc t ion losses can be neglected (or that neglect ing them w i l l result i n conservative answers) is obv ious ly both unwarranted and dangerous. F o r this reason, a mode l was sought wh i ch wou ld prov ide a s imple but rat iona l ly defencib le representat ion o f junc t ion losses. The representation is desired to be s imple i n order to result i n on ly a smal l increase i n computer time and i n order to a l low rap id mod i f i ca t ion o f this p rogram segment as more accurate i n fo rmat ion on junct ion losses becomes avai lable. The representat ion is desired to be rat iona l i n order to a l l ow s imple conclusions to be drawn f rom the mode l . Spec i f i ca l ly , answers to the fo l low ing two questions are sought: "Unde r what condi t ions do junc t ion losses represent a s igni f icant part o f the network 80 response?" and " W h e n can junct ion losses be neglected w i th the reasonable expectation that their neglect w i l l be conservat ive?" The mode l suggested below meets these cr iter ia. The mode l has been tested on s imple series p ipe l ine systems and been found to give very good agreement wi th current ly used techniques. In addi t ion, appl icat ion o f the head loss mode l results i n on ly a sma l l increase in compuat ion t ime for the solut ion o f the network as a whole. W h e n losses are inc luded at a junc t ion , a number o f head loss relations must be spec i f ied in addi t ion to the usual junc t ion equations. F o r a ord inary node, there are now 2n + 2 unknowns to be calculated: the head and discharge at the boundary section o f each p ipe and the head and discharge govern ing the external f low. T o calculate these unknowns, there are n independent characterist ic equations, a junc t ion continuity equat ion and an auxi l iary re lat ion govern ing the external f low. Th i s means that n add i t iona l head loss relat ions must be writ ten. The fo rm o f these equations is Q P i 2 Q P i J P i 2 g A . J Pj v y y 2 8 A j in wh i ch the subscript i designates a p ipe or external f low wh i ch discharges towards the junc t ion and the subscr ipt j ident i f ies a p ipe wh i ch discharges away f rom the j unc t i on . 2 5 T he parameter ident i f ies the head loss coeff ic ient between in f l ow p ipe i and out f low p ipe j and w i l l be assumed to be a known constant i n this d eve l opmen t 2 4 In order to solve the boundary condit ion o f a arbi trary network junc t ion when losses are inc luded, each p ipe and external f low must appear at least once in a head loss equat ion; in o rder to develop n 2 5 WheTeas prev ious ly it was on ly necessary to work w i th assumed discharge directions, Equat ion (44) relates actual in f lows to actual outf lows. " T h e d i f f i cu l ty o f mode l ing the head losses at junct ions l ies i n the fact that the loss coeff ic ients are not easy to obtain. The K... va lues depend on many factors and on ly for rectangular branchings are the head loss 1 J coeff ic ients reasonably constant 81 independent equations o f this type, a number o f p ipes may have to appear more than once. A t f irst glance, it appears that a system o f 2n + 2 non l inear equations w i l l have to be assembled and so lved for each junc t ion o f this type i n the network. Such a procedure wou ld result i n a vastly increased execut ion t ime compared to an equivalent network where such losses are not inc luded. In add i t ion, the existence o f l i nked boundary condit ions cou ld qu ick ly make the p r ob l em a lmost insoluble. Fortunate ly , it does not appear to be necessary to solve this large system o f non l inear equations. I f it is assumed that head losses at the junc t ion affect the magni tude o f heads and discharges at the node, but not the direct ion o f f low in any one pipe, the so lut ion o f this boundary cond i t ion can be qu ick ly and accurately obtained. W i t h this assumption, the so lut ion o f this boundary cond i t ion can be obtained as fo l lows. 1. Assume in i t ia l ly that the junc t ion head losses are negl ig ib le and calculate the values o f head and discharge i n each p ipe and the value o f each external f low in the usual manner. Th i s a l lows the in f l ow and out f low p ipes to be ident i f ied. 2. Ass ign the junct ion head (i.e., H p in Equat ion (5)) to the boundary section o f the p ipe or external f low hav ing the largest discharge value o f a l l p ipes discharging toward the node i n quest ion. U s e this head and the calculated value o f discharge i n this p ipe to complete ly determine the left hand side o f Equat ion (44) for at least one associated f low leaving the junc t ion . F o r the r ight hand side o f this equat ion, the C + or C equat ion can be substituted i n as appropriate. Th i s procedure produces a quadrat ic i n Qp^ wh i ch can be expl ic i t ly solved. Once is obtained, H p j can be calculated by back substitut ion. In this manner , it is possib le to determine the head and discharge i n each p ipe or the external f low. 82 3. However , the cont inuity equat ion at the junc t ion w i l l no longer be complete ly satisf ied. Th i s equat ion can now be thought o f as a funct ion o f a single unknown (the -discharge i n the p ipe wi th the largest in f low). Newton ' s method can then be app l ied to the cont inuity equat ion to produce an improved value o f the largest in f low discharge. The head i n this p ipe can then be recalculated and, successively, the head and discharge in every other p ipe. 4. Step 3 can be repeated as many t imes as necessary to achieve the des ired accuracy for each f low and discharge. Exper ience has shown that, even for junct ions hav ing unusual ly large values o f head losses, on ly one or two iterations are necessary to obta in disharges and head values wh i ch agree wi th the exact numer ica l so lut ion to w i th in 0.1%. Thus , assuming that the K y valves are known, junct ions hav ing head losses due to fr ict ion can be accomodated w i th in the f ramework o f the boundary condi t ion solut ions presented in this chapter w i th on ly m ino r mod i f i ca t ion to the usual procedures. E ven i f the exact values o f the loss coeff ic ients are not known , incorporat ing reasonable K~ values in the analysis o f speci f ic networks shou ld a l low the analyst to begin to assess the importance o f these losses on the transient behav ior o f the network as a whole. 4.7 C O N C L U S I O N S The study o f hydrau l i c networks i n the transient state is essentially a study i n the boundary condi t ions that make up the network. I f the boundary condit ions are rat ional ly formulated and solved, then the response o f the network can be ef f ic ient ly predicted. Th i s chapter has emphas ized a general ized approach to the formulat ion and so lut ion o f the network boundary condit ions. A number o f specif ic terms and concepts 83 have been deve loped wh i ch facil itate the general ized approach. In add i t ion, some o f the common hydrau l i c boundary condit ions have been treated in a more comprehensive manner than is usual i n the l iterature. Th i s procedure al lows special cases to be extracted f r om the analysis by s imply selecting specif ic constants rather than by solv ing each case ind iv idua l l y . The p r imary developments descr ibed in this chapter are l isted below. 1. A new boundary cond i t ion classif icat ion system is developed: the complexity o f a node is the number o f external f lows connected to it wh i le the nodal dependence number o f a boundary cond i t ion is the number o f nodes required to def ine the element completely. Th i s c lassif icat ion ignores the concept o f the degree o f a node, since the number o f p ipes that connect to a node does not inf luence the solut ion procedure. 2. The classi f icat ion system for boundary condi t ions leads logica l ly to a classif icat ion system for networks as we l l . Boundary condi t ions that can be so lved indepenent ly o f one another are termed isolated and are the on ly type o f e lement possible in ordinary networks. Cu r ren t solut ions to the network p rob l em work exclus ively with networks o f this type. However , when complex nodes exist i n a network, boundary condi t ions can become linked and must be so lved s imultaneously. The terms link length and group number have been introduced to quant i fy the var iables i nvo l ved i n a set o f l i nked boundary condit ions. 3. S imp le and ord inary one -node boundary condi t ions are also discussed. The general C ^ node equat ion is deve loped and app l ied to a number examples o f c ommon external f low relations. T he concept o f a linear reservoir is developed in order to mode l a f in i te -a rea reservoir when the inert ia and fr ic t ion effects w i th in it are s ign i f i cant The general energy diss ipat ing boundary cond i t ion is formulated to a l low devices such as valves and or i f ices discharging to l inear reservoirs to be mode led . These general formulat ions a l low a wide range o f hydrau l i c devices and 84 boundary condi t ions to be real ist ical ly represented. In addi t ion, the equations descr ib ing a i r chambers and the power fa i lure condi t ion for pumps are wri t ten in a general manner. 4. O rd ina ry t w o - n o d e boundary condi t ions can often be solved by relat ively s imple extensions o f their one - node counterparts. Thus, energy diss ipat ing f lows and pumps undergo ing power fa i lure are qu ick ly developed. The pipe replacement element is solved using either a l umped inert ia or f inite di f ference representation. The so lut ion to the pressure reduc ing va lve boundary cond i t ion is inc luded to i l lustrate the procedure for so lv ing an automat ic valve boundary condi t ion. 5. F ina l l y , the chapter concludes w i th b r i e f discussion o f how the prev ious ly developed boundary condit ions can be extended to more general situations. Three extensions — internal boundary condit ions, complex nodes and boundary condit ions invo lv ing more than two nodes — are br ie f ly discussed. The final extension considers the quest ion o f how to inc lude junct ion losses at a general network node. A procedure is proposed wh i ch a l lows junct ion losses to be convenient ly and accurately calculated w i th on ly a smal l increase i n computat ional e f f o r t The quest ion o f how a l l . o f these boundary cond i t ion solut ions can be imp lemented wi th the context o f a general network p rogram is discussed in the fo l l ow ing chapter. 5. D H V F . I D P M F N T O F T H F N F T W O R K M O D E L 5.1 I N T R O D T JCTTTON T h e object o f this thesis is to deve lop an eff ic ient and re l iab le procedure for ana lyz ing large water d is t r ibut ion systems in the transient state by computer. In order to accompl i sh this a im , three dist inct stages must be completed. F i rst , a mathematical mode l has to be selected or der ived. Second, an a lgor i thm or procedure must be found for so lv ing the mathemat ica l mode l . T h i r d , a convenient representat ion has to be selected for the data on wh i ch the a lgor i thm is to operate. 2 7 Chapters 2, 3 and 4 have deve loped the equations govern ing the steady and unsteady f low o f f l u i d i n closed conduits as we l l as the equations descr ib ing the behav ior o f the var ious network boundary condit ions. These chapters have therefore completed the f irst stage o f the p rogram deve lopment by p rov id ing a fu l l descr ipt ion o f the mathemat ica l mode l used to character ize a hydrau l i c network. In addi t ion, the use o f the method o f characterist ics to solve the transient equations and the associated boundary condi t ions has been descr ibed in Chapters 2 and 4, wh i le Chap te r 3 has cons idered the so lut ion o f the steady-state network p rob lem. In effect, then, the second stage has also been completed. The th i rd stage i n the p rogram deve lopment has, unt i l now, received relat ively l i tt le attent ion i n hydrau l i c appl icat ions. The ma in concern at this stage is to formulate procedures and data structures wh i ch a l low the solut ions deve loped i n the first two stages to be eff ic ient ly imp lemented. I f a p rog ram is to have general appl icat ion to a w ide var iety o f systems, the data structures selected shou ld a l l ow the f ina l solut ion to be easily calculated wh i le keeping the input requ i red by the p rogram to a bare m i n i m u m . The entire quest ion o f how to imp lement the so lut ion procedures presented The s e three stages are ident i f ied i n Berzt iss (1971). 85 86 so far is the subject o f the current chapter. Th i s chapter w i l l be concerned wi th three speci f ic aspects o f the imp lementat ion quest ion. F i rs t , the p rob l em o f how to represent the topology o f the network w i l l be considered. The proposed procedure w i l l a l low the constants i n the boundary cond i t ion equations presented i n Chapte r 4 to be readi ly obta ined. Fo l l ow ing this, the prob lems ar is ing f rom l i nked groups o f boundary condit ions w i l l be discussed. L i n ked boundary condit ions occur i n complex networks and must be easily ident i f ied and solved i f the transient so lut ion for the network is to be calculated. F ina l l y , the di f f icu l ty o f choosing an appropr iate t ime step for the analysis w i l l be cons idered under the topic o f the network discret izat ion. A f lex ib le a lgor i thm for select ing the t ime step is presented wh ich al lows ind iv idua l networks to be eff ic ient ly analyzed and should a id i n the development o f general rules for the transient analysis o f a l l hydrau l i c systems. 5.2 N E T W O R K R E P R E S E N T A T I O N One o f the most important funct ions o f the computer representat ion o f the network is to accurately describe the hydrau l i c system wh i le requ i r ing a m i n i m u m o f input data . 2 8 The solut ion to the network boundary condit ions requires that every p ipe wh i ch jo ins to any node i n the network to be examined and processed i n some way at least once. Fo r eff ic iency, a data structure shou ld be chosen that on ly stores " in fo rmat ion that is real ly needed, and stores it in such a way that the operations required by the a lgor i thm can be done qu i ck l y " (Baase, 1978). " T h e concern at this stage is w i th the descr ipt ion o f the network topology and not wi th the basic p ipe and boundary condi t ion data pe r se. Thus, o f pr imary importance is the quest ion o f how the var ious p ipes and boundary condi t ions are jo ined together. It is assumed that the phys ica l characterist ics o f the network elements, such as p ipe lengths, f r i c t ion factors, condui t diameters, reservoir d imens ions etc. must be suppl ied to any p rogram i f the analysis p r ob l em is to be un ique ly speci f ied. 87 Cons ider , for example, the general o ne - node boundary cond i t ion wh ich is represented by Equa t i on (5) i n Chapte r 4. The constants B c and C c i n this expression are evaluated by Equat ions (6) and (7) wh i ch invo lve summations over every p ipe wh i ch connects to a g iven node. Thus , in order for B c and C c to be calculated, there must be a convenient way o f ident i fy ing a l l the pipes that j o i n w i th each node in the network. The re are a number o f ways o f ach iev ing this purpose. O n e o f the most popu lar is an " i ndex ing" system, first proposed by Streeter (1967) and later improved by W y l i e and Streeter (1983). Th i s technique uses an index vector wh i ch requires, "as a m i n i m u m , " the fo l lowing in format ion at each node: the node number, the number o f elements attached to the node, and for each e lement attached: the element number w i th an indicator as to whether it is the upstream or downstream end o f the element, and the number o f the e lement section adjacent to the n o d e . " However , this tradit ional approach has several drawbacks. First, i f the network being studied is mod i f i ed i n any way, such as by add ing an extra p ipe or node, the index vector must be careful ly altered i n a number o f places. Second, i f the level o f discret izat ion is var ied, say by selecting a smal ler t ime step for the analysis, almost the entire index vector must be changed — a t ime consuming and e r r o r - p rone task. 3 0 In add i t ion, the long lists o f integers requ i red by the index vector are tedious to prepare, d i f f i cu l t to check and in f lex ib le to app ly when large networks hav ing a wide var iety o f boundary condit ions are to be mode led . F ina l l y , a s igni f icant amount o f p re l im ina ry work and organizat ion must be done pr io r to the computer analysis in order to determine the elements that make up the index vector. F o r al l o f these reasons, a s imp le r and more flexible representat ion than the index vector approach was " Q u o t e d f r om Wy l i e ' s 1983 paper. 5 0 I f an a lgor i thm for choos ing the t ime step is to be developed, it is essential that the network descr ipt ion be independent o f the value o f the t ime step. Th i s is important since the program, and not the user, shou ld make the final selection o f the t ime step. 88 sough t Baase (1978) describes two other c ommon data structures that are useful for represent ing networks in computers. T he first uses a character izat ion ca l led an adjacency matrix In the adjacency matr ix , the matr ix entry i n the row and j ^ co lumn denotes a p ipe f r om node i to node j . I f such a p ipe does not exist, a zero or other suitable number is stored at this locat ion. A l t hough this is quite a power fu l approach in some appl icat ions, the adjacency matr ix requires a larger amount o f storage and processing when app l ied to hydrau l i c ne tworks 3 1 than some other data structures, and is therefore not we l l suited to the present purpose. The way the network w i l l be represented i n this thesis is s imi l ia r to the adjacency list structure descr ibed by Baase. T he actual fo rmu la t ion that w i l l be presented is deve loped i n the steady-state analysis p rogram o f E p p and Fow l e r (1970). Th i s character izat ion is s imple but power fu l and al lows eff ic ient programs to be wri t ten and developed. In add i t ion, the structure can be developed f r om a m i n i m u m o f input data. Since this l ist processing technique has been rarely app l i ed i n transient appl icat ions, it w i l l be descr ibed i n some detai l here. The on ly essential i n fo rmat ion that is requi red to complete ly describe the topo logy o f a network is the number o f the nodes at both ends o f each p ipe and the ident i ty o f the nodes wh i ch def ine each boundary cond i t ion i n the network. I f a p ipe connects node i to node j , then the nodes i and j are said to be adjacent to each other and incident w i th the p ipe j o i n i ng them. Thus, by on ly know ing the nodes inc ident w i th each p ipe i n the network, a l ist o f the p ipes wh i ch j o i n each node must be developed. Once this l ist has been assembled, the constants def ined in Equat ions (6) and (7) i n Chapte r 4 can be easily computed. Th i s process can be 3 1 I n a hydrau l i c network, usual ly on ly a few pipes connect w i th any one node. Th i s means that the major i ty o f the terms i n the adjacency matr ix w i l l be zero and are not strict ly needed to describe the network. 89 accompl i shed i n the fo l low ing three steps. 1. The degree o f each node must be first calculated. Th i s va lue is obtained by successively consider ing each p ipe i n the network; every t ime a given node attached to either end o f a p ipe is encountered, the current va lue o f the degree o f this node is increased by one. A zero value for the degree o f each node must in i t ia l ly be assigned to each node in the network. 2. A f te r the degree o f each node has been determined, a l ist o f pointers can be obta ined. The pointer to the first node i n the network is assigned a value o f one; the value o f the pointer at each successive node is de f ined by a recurs ion re lat ion. Tha t is, i f P . is the value o f the pointer for node i, and D . is the degree o f this node, then P . = P . , + D . , (1) I l - l l - l The pointers w i l l be used to locate the pipes connected to a g iven node in an inc ident p ipe l i s t 3. F ina l l y , a second pass is made through each pipe i n the network. F o r each p ipe, the nodes attached to this p ipe are found f rom the input data as i n Step 1. F o r each o f these end nodes, the current value o f the po inter is obta ined and the associated p ipe numbe r stored at this locat ion i n the inc ident p ipe l i s t T he va lue o f the pointer must then be incremented by one to po in t to the next (vacant) locat ion in the p ipe l i s t A f te r each p ipe i n the network has been considered, the value o f the pointer shou ld be reset to the beg inn ing o f the p ipe lists. F i gu re 12 and the associated Tab le 1 i l lustrate how this procedure is app l ied to a s imp le p ipe network. The p ipes and nodes are numbered consecut ively starting at one for s impl ic i ty , a l though the actual number ing app l ied to any speci f ic network is arbi trary. T h e table depicts the final appearance o f the data structure after each o f 90 F i gu re 12. S imp le Ne twork Examp le Node Degree Po in ter Incident P ipe L i s t (Step 1) (Step 2) (Step 3) 1 2 1 1 1 3 2 2 3 3 3 1 ' 3 5 4 3 4 6 6 2 7 4 8 5 9 7 4 3 10 1 0 3 1 1 5 1 2 6 5 2 13 1 3 6 1 4 7 Tab le 1. Da ta Representat ion o f Ne twork in F igure 12 the steps i temized above have been completed. The f ina l inc ident p ipe list can be stored in a s i ng l e -d imens ioned array by us ing the pointer for any g iven node to locate the first p ipe attached to this node. The smal l superscr ipt numbers in the inc ident p ipe l ist denote the actual array locat ion where a part icu lar p ipe number w i l l be stored. F u z y et a l . (1975) suggest us ing a two -d imens i ona l equiva lent o f the inc ident p ipe list but imp l y that such a l ist must be read as input data by the program. In add i t ion, such a two -d imens i ona l l ist wastes computer memory by storing many unnecessary zeros i n the matr ix. T he two -d imens i ona l system wi l l also have restr ict ions on the m a x i m u m degree o f any node whi le no such restr ict ions exist i n the data structure proposed here. 91 Once the p ipes that are inc ident w i th each node i n the network have been clear ly ident i f ied, the summations requ i red to calculate the node constants B c and C c can be easily per fo rmed. The value o f these constants w i l l , i n turn, a l low any isolated boundary cond i t ion i n the network to be solved by the methods out l ined in Chapter 4. The on ly d i f f i cu l ty arises when the boundary condit ions become l inked; the procedures app l i cab le to this case w i l l be discussed later i n this chapter. 5.2.1 INTERNAL SECTIONS In add i t ion to the equations wh i ch def ine the constants i n the boundary cond i t ion expressions, the method o f characteristics requires that the hydrau l i c var iables at the internal sections o f each p ipe be calculated dur ing the solut ion procedure. The usual way o f ident i fy ing the characteristic sections in each p ipe is to use a doub le -subscr ip t ing system: the first subscript is used to ident i fy the p ipe numbe r i n wh ich the section occurs and the second ident i f ies the locat ion of the ind iv idua l section i n the p ipe. Thus, the head and discharge values requ i red by the method o f characterist ics can be stored i n a two -d imens i ona l array. T h e number o f rows in this matr ix must be equal to the greatest number o f p ipes i t is desired to analyze, wh i le the number o f co lumns must be equal to the m a x i m u m number o f sections a l lowed in any p ipe. T h e obv ious p rob lem w i th this approach is the poor use it makes o f computer storage. The number o f characterist ic sections that must be used to mode l a large hydrau l i c network is very gTeat and it is therefore desirable to improve on the above procedure. The common ly used two-d imens iona l array requires a m i n i m u m storage for each basic hydrau l i c va r i ab l e 3 2 to be equal to the m a x i m u m number o f p ipes t imes the m a x i m u m number o f sections per p ipe. " T h e basic hydrau l i c variables inc lude the value o f head and discharge at both the current t ime step and the prev ious t ime step for each section in the network. 92 The data structure that w i l l be presented here on ly requires the max imum number o f pipes t imes the average number o f sections pe r p ipe as a m i n i m u m storage per var iable. F o r a network wi th a typical d is t r ibut ion o f p ipe lengths, it is est imated that a system two to three t imes larger can be represented i n the same computer memory by the proposed approach. In addit ion, the doub le -subsc r ip t ing system places an unnecessary restr ic i t ion on the max imum number o f sections that can occur i n any one p ipe; this restr ict ion does not exist in the data structure to be presented here. The representat ion o f the basic hydrau l i c variables pre fer red in this thesis is to use a one -d imens i ona l array w i th an associated l ist o f pointers. The pointers are used to locate the characterist ic sections associated w i th each p ipe i n a manner analagous to the nodal pointers used to describe the network itself i n the data structure presented above. The procedure used to develop this data structure can also be completed in three steps. 1. B y using the discret izat ion a lgor i thm to be discussed later i n this chapter, the number o f sections associated wi th each pipe i n the network can be determined. 2. A l ist o f pointers to the sections o f each pipe can then be computed. The po inter to the f irst section o f the f irst p ipe is assigned a value o f one; the pointers to subsequent p ipe sections are determined by the value o f the pointer to the preced ing p ipe p lus the numbe r o f sections in the preced ing p ipe. 3. F ina l l y , the hydrau l i c var iables for the ind iv idua l p ipe sections can be assigned to the vector wh i ch is to store this basic in format ion. The upstream boundary sect ion o f the p ipe is assigned to the f irst locat ion i n this array and the process cont inued through the interna l p ipe sections. The boundary section o f the downstream end o f the p ipe w i l l be assigned to 93 the f ina l storage locat ion reserved for that part icu lar p ipe. Once this process has been completed for each p ipe, the boundary and internal sections o f every p ipe i n the network can be easily ident i f ied. A n example o f a s imp le network is shown in F igure 13 and Tab le 2. F i gu re 13 shows the actual arrangement o f p ipes wi th their assumed f low di rect ion wh i le Tab le 2 i l lustrates P ipe Number Po inter Character i s t i c Sect ions of Sect ions (Step 1) (Step 2) (Step 3) F igure 13. S imp le Ne two rk Examp l e 94 the associated data structure. Th i s table ident i f ies the var ious steps i n constructing the computer representat ion for this network. It shou ld be noted that the number ing o f the nodes and pipes is arbi trary i n this system, a l though the example shows consecutive number ing beg inn ing wi th one. I f a p ipe number is not present, the number o f sections associated wi th this p ipe is zero and the characterist ics section array need not be scanned. P ipe 7 in Tab le 2 is inc luded to i l lustrate this s ituation. The important po in t about this data structure is that no in fo rmat ion about the network is lost by the one -d imens i ona l characterizat ion. A n y operat ions pe r fo rmed on the tradit ional two -d imens i ona l representation can also be per fo rmed on the one -d imens i ona l array. In part icular, calculat ing the so lut ion for the internal sections using Equa t i on (13) o f Chap te r 2 as wel l as determin ing the C + and C constants requ i red to determine the boundary condi t ions can be easily accompl i shed wi th either representat ion. 5.3 ROT JNDARY CONDITIONS REPRESENTATION It has been shown i n the prev ious section that the s imple knowledge o f wh i ch nodes are inc ident w i th each p ipe is suff ic ient i n order to complete ly determine the topology o f the network. If, i n addi t ion, the number o f sections requi red to app ly the method o f characterist ics to each p ipe is known, a one -d imens i ona l array represent ing a l l o f these sections can be developed. D i r e c t app l i cat ion o f the method o f characteristics us ing the in fo rmat ion i n both o f these structures a l lows the head and discharge to be computed at every inter ior section i n the network and the boundary cond i t ion constants B c and C c to be determined. The one rema in ing d i f f icu l ty is to use the in format ion the me thod o f characterist ics has p rov ided i n order to solve a l l o f the network boundary condit ions. 95 N o d i f f i cu l ty i n so lv ing the network boundary condi t ions arises when the hydrau l i c e lements are connected to ord inary nodes. Such devices are isolated by the expl ic i t nature o f the method o f characteristics and can be solved independendy o f al l the other boundary condit ions in the network by the procedures out l ined i n Chapte r 4. However , d i f f i cu l t ies arise i n complex networks. In these networks, an interact ion exists between the external f low equations at a l l nodes wh i ch have a complex i ty greater than one. The nature o f the interact ion between l i nked boundary condi t ions and the quest ion o f how these l i n ked elements can be ident i f ied is the subject o f this section. F o r s impl ic i ty in this discussion, several assumptions w i l l be made. F i rst , al l boundary condi t ions w i l l be assumed to occur on ly at (or between) specia l ly designated nodes. A n y node not specif ical ly referred to w i l l be assumed to be s imple . Second, a l l boundary condi t ions in this section w i l l be assumed to have a noda l dependence number o f one or two. N o serious d i f f icu l t ies aTe raised by t h ree -node boundary condit ions, but the notat ion and logic becomes more invo lved and the procedures are h igh ly device dependent Since boundary condi t ions hav ing a noda l dependence number o f three or more are rare, they can be convenient ly separated and treated as extensions o f the s impler boundary condit ions. F ina l l y , the usual assumpt ion is made that any node i n the network may have arb i t rary degree. The p r ob l em o f ident i fy ing a l inked group o f boundary condi t ions i n a network is s imi lar , in many ways, to the p rob l em o f descr ib ing the layout and arrangement o f the network pipes. The goal remains the same — to un ique ly relate each e lement to the other components in a manner that a l lows the solut ion to be easily constructed but wh i ch requires a m i n i m u m o f input data. Thus , i n e f fec t the p rob l em is to specify a network (the l inked boundary condit ions) w i th in a larger network (the pipes themselves). The elements i n c ommon w i th both o f these systems are the network nodes; these nodes must complete ly def ine both the p ipe and the boundary cond i t ion connect ions i f the network is to be un ique ly speci f ied. 96 I n order to develop the ideas re lat ing to l i nked boundary condit ions, an example w i l l be presented and discussed i n some detai l . T he data structure descr ibed i n this example is, i n many ways, s imi lar to the representations deve loped ear l ier i n this chapter. The example shou ld p rov ide suff ic ient details to show how this structure is deve loped and app l ied to a l low a formal descr ipt ion to be omit ted. F igu re 14 shows a smal l section o f a network i n the v ic in i ty o f a t w o - p u m p booster station. In add i t ion to the two pumps, this sect ion o f network contains a reservoir for temporary water storage and for the transient protect ion o f the pumps, two known external discharges, one va lve or energy diss ipat ing boundary cond i t ion and one p ipe rep lacement e l emen t The p ipe replacement e lement ( P R E ) is used i n the usual way to mode l a short condu i t that can not be easi ly represented by the method o f characteristics. It is assumed that the contro l valve can be used to divert a por t ion o f the f low into the part o f the network connected to node 3, thus avo id ing the head increase caused by P u m p A . F o r s impl ic i ty , i t w i l l be assumed that both pumps are wo rk i ng at a constant 5.3.1 E X A M P L E N E T W O R K W I T H L I N K E D B O U N D A R Y C O N D I T I O N S (5) PIPE 0(1) (4) (3) PIPE F igu re 14. E xamp l e Ne twork w i th L i n k e d Boundary Cond i t i ons 97 speed i n the i r fo rward direct ion, and that transients have been introduced through the nearby segments o f the network. F o r example, the known external discharges at nodes 6 and 7 may represent a rap id ly occurr ing fire flow demand; the purpose o f the analysis is to predict how the pump ing station responds to this suddenly imposed f low. I f both pumps are operat ing at constant speed, there are no two -pa ramete r boundary condi t ions present in this network segment Tab le 3 i l lustrates the input data requ i red to locate each o f the • boundary condi t ions for the network shown in F igu re 14. The ident i f icat ion number associated wi th each boundary cond i t ion can be assigned by the p rogram and does not need to be p rov ided as i n pu t F o r convenience, the boundary cond i t ion number ( in brackets) are also shown in F igu re 14. Each o f the boundary condi t ions w i l l have an associated nodal dependence number ( N D N ) and a list o f its de f in ing nodes. G i v e n the in fo rmat ion contained i n Tab le 3, the goal is to construct lists showing wh ich boundary condi t ions are l i nked , and wh ich are isolated. The B.C Number D i s c r i p t i on N D N Def in ing Nodes 1 K n o w n D e m a n d 1 7 2 K n o w n D e m a n d 1 6 3 P u m p A 2 2, 4 4 P u m p B 2 1, 2 5 Surge Tank 1 4 6 Con t r o l Va l ve 2 2, 3 7 - P R E 2 5, 7 Tab l e 3. Da ta Requ i r ed to Locate Boundary Cond i t i ons i n F igu re 14 98 Sode Complex i ty B .C . L i s t Adjacency L i s t | N | | W | E l (1) (2) (3) (4) (5) (6) (7) 1 1 4 2 8 2 3 3. 4, 6 4. 1. 3 4 4 8 3 1 6 2 8 4 2 3. 5 2 8 5 1 7 7 4 6 1 2 1 7 2 1. 7 5 2 2 4 Tab le 4. Da t a Structure for Ne twork Shown in F i gu re 14 procedure deve loped is considerably more compl icated than that used ear l ier to def ine the connectedness o f the network pipes, a l though the structure is deve loped a long analagous l ines. The stages o f the deve lopment are expla ined wi th reference to the network in F igu re 14 and are i l lustrated in Tab le 4. 1. The first step is to determine the complex i ty o f each node. Th i s is done by cons ider ing each boundary cond i t ion in turn and enumerat ing the references to each node in the ne two rk . " I f the m a x i m u m complex i ty o f any node i n the network is one, no boundary condi t ions are l inked. In such a case, the network is ord inary and al l boundary condi t ions are isolated. Howeve r , i f the m a x i m u m complex i ty is two or greater, as it is " T h i s procedure is s im i la r to determin ing the degree o f each node, except that now boundary condi t ions are counted instead o f pipes. The complex i ty o f each node can also be used to determine a l ist o f pointers so that the in fo rmat ion stored i n the boundary cond i t ion l ist wh i ch w i l l next be deveoped can be retained in a one -d imens i ona l array. Th i s step is not shown in Tab le 4 or referred to i n subsequent sections since the procedure for deve lop ing such a l ist has a l ready been descr ibed i n detai l . 99 i n F igu re 14, the network is complex and at least two boundary condit ions are l inked. 2. Once the complex i ty o f each node has- been determined, a second pass can be made through the l ist o f nodes attached to each boundary condi t ion. O n this pass, every t ime a g iven node is encountered the boundary cond i t ion associated wi th that node is added to a boundary condi t ion list ( 'B .C L i s t ' or co lumn (3) i n Tab le 4). Once this step has been completed, the part icu lar boundary condit ions associated wi th each node wi l l be determined. 3. The th i rd step is to develop the set o f nodes that are adjacent 3 4 to each node o f complex i ty one or greater. Th i s l ist is now ident ica l i n funct ion to the adjacency list structure deve loped by Baase. Th i s l ist is developed by cons ider ing each node i n turn. F o r each node, the boundary cond i t ion list is scanned. Eve ry t ime a boundary cond i t ion is encountered, i t def in ing nodes are checked. A n y def in ing node for a boundary cond i t ion that is d i f ferent f r om the node be ing considered, and has not already been inc luded, is added to the adjacency l i s t C o l u m n (4) o f Tab le 4 shows the adjaceny l ist for the boundary condi t ions o f the network in F igure 14. ( Fo r eff ic iency, this step wou l d norma l l y be completed s imultaneously with Step 2) 4. The next step is to ' ident i fy the l i n ked components o f the network. Th i s is done by complet ing a "b read th - f i r s t scan" o f the adjacency l ist structure i n Tab le 4. F i rst , the node hav ing the highest complex i ty value in the network is ident i f ied. In Tab le 4, this is node 2. The nodes that are 3 4 I n this context, adjacent means that a boundary cond i t ion connects o r l inks the two adjacent nodes. Fo rmer l y , this te rm was def ined i n terms o f p ipe connections but this def in i t ion has no mean ing here. N o confus ion shou ld arise, however, since the context shou ld make the part icu lar mean ing clear. 100 adjacent to this node are now ident i f ied (nodes 4, 1 and 3) and are added, a long wi th the first node, to a new array cal led a linked node list (this is not shown i n Tab le 4). Each node i n the l i nked node l ist w i l l now be considered in turn. Each t ime a node is considered, its own adjacency l ist is scanned. Eve ry t ime a new node is discovered in the adjacency lists o f these nodes that has not yet been considered, it is added to the end o f the l i nked node l i s t In Tab le 4, no new nodes w i l l be ident i f ied in this way since the adjaceny lists o f nodes 1, 3 and 4 contain on ly one e lement — node 2 — wh i ch has already been considered. Once a l l the entries i n the adjacency list for a g iven node have been considered, the node is marked so it w i l l not be considered again (this is also how new nodes can be ident i f ied). Once al l o f the nodes in the l i nked node list have been considered, this part icu lar l i nked boundary cond i t ion set has been ident i f ied. Thus, for the network shown i n F i gu re 14, the l i nked node l ist w i l l contain the nodes 1, 2, 3, and 4. Th i s process can be repeated for each complex node that has not yet been marked unt i l a l l o f the complex nodes i n the network have been added to a l i nked set It is often convenient to assign a reference number to each o f these sets for future reference. Once al l the l i nked nodes have been ident i f ied, the l i nked boundary condit ions associated w i th them can easily be obta ined f rom co lumn (3) o f Tab le 4. 5. N o w that the l i nked elements (both nodes and boundary condit ions) have been ident i f ied, the cardinal i ty o f these sets can be computed. W h e n this step is completed, co l umn (5) and (6) in Tab le 4 resu l t The entries i n these co lumns can then be summed and the number o f two-parameter boundary condit ions added to this total. T he final result is the number o f equations E . that must be so lved to complete ly def ine each l i nked se t 101 Th i s number , as we l l as the l inked set reference number, shou ld then be associated wi th each node in the set so that members o f each set can be easi ly ident i f ied. These five steps a l low the l inked set o f boundary condi t ions i n every network to be systematical ly ident i f ied. Desp i te the long verbal descr ipt ion, the steps i temized above can be qu ick ly pe r fo rmed on the computer. Usua l l y the number o f l i nked elements i n a hydrau l i c network is quite smal l , but the behav ior and response o f these elements is often cr i t ical f r om an analysis po int o f v iew. The abi l i ty to recognize, ident i fy and analyze these elements represents a power fu l add i t ion to the field o f transient network analysis. N o w that a procedure for ident i fy ing a set o f l i nked boundary condit ions has been developed, the remain ing d i f f i cu l ty is to assemble and solve the associated boundary cond i t ion equations. The procedure required to complete this process is stra ightforward and w i l l be i l lustrated on the example network in F igu re 14. It has just been shown that the boundary condi t ions invo lved in nodes 1, 2, 3 and 4 are l i nked ; the so lut ion procedure appl icable to this l i nked set w i l l now be developed. In order to solve l inked boundary condit ions, it is necessary to work w i th the fundamenta l boundary cond i t ion equations. Tha t is, the f o rm o f the boundary cond i t ion equat ions used can no longer assume that Equat ion (5) o f Chapte r 4 appl ies i n the usual way, since more than one external f low jo ins wi th every complex node i n the l i n ked set. T he compound nature o f the external f lows at these comp lex nodes must now be accounted for in the solut ion procedure. Th i s discussion can best be i l lustrated by example. The fundamenta l boundary cond i t ion equations for the elements in the example network w i l l now be writ ten. In the equations that fo l low, the subscript numbers on the head terms refer to the assigned node number, wh i le the 102 subscr ipt numbers on the external f lows refer to the boundary cond i t ion number associated wi th each external f low. F o r p ump A , i n wh i ch is the dynamic head produced by pump A and H y is the head loss due to this pump 's valves. S imi la r l y , H P 2 = H P l + H d 4 " H v 4 <3> can be wri t ten for pump B. In both o f these expressions, the dynamic head term and the valve loss te rm are on ly funct ions o f the external f low through the part icu lar pump in quest ion. The fundamenta l equat ion for the surge tank is obta ined by comb in ing Equat ions (12) and (14) o f Chapte r 4. The result is Q e x t , - s ^ E s > > P < - C b " B c < W = 0 W A s imi la r equat ion can be wr i t ten for the control valve: namely, Q e x u " S ' ( T E s > > P 2 " * V = ° : ( 5 ) In order to s impl i fy these equations, it is possible to app ly the new nodal equations. These equations are obta ined f r om Equat ion (5) o f Chap te r 4 by cons ider ing the mul t ip le nature o f the external f lows at the complex nodes. The resul t ing equations can be wr i t ten as H P 2 = C c 2 " V Q « U + <<ext3 " Qexu> <7> H P 3 = C c 3 + B c 3 Q e x U (8) 103 and H P 4 = C c « " V Q e x t 5 - % n ? • < 9 > Equat ions (2) through (9) represent the complete set o f equations for the boundary condit ions associated wi th nodes 1-4. In order to solve this l i nked group, these equations can be used direct ly as a set o f e ight equations and eight unknowns. However , it is possible to produce a more eff ic ient solut ion by e l iminat ing the head terms f r om this equations. Th i s is possib le because the re lat ion between head and external f low at a complex node remains l inear. I f a l l the head terms are e l iminated f r om equations (2) through (5), a system o f four equations and four unknows results. Since the part ia l der ivat ives o f al l o f the terms i n this re lat ion are easily determined, and since good in i t ia l estimates o f the external f low values are avai lable f r om the so lut ion at the prev ious t ime step, the N e w t o n - R a p h s o n approach is a useful technique for so lv ing this system. Thus, i t has been shown that the solut ion to comp lex networks can be obta ined in a s imi lar manner to the solut ion o f ord inary networks. The compl icat ions arise i n ident i fy ing and assembl ing the l i nked boundary cond i t ion equations; however, these di f f icul t ies can be overcome by careful data man ipu la t ion and a logica l s t ep -by - s t ep solut ion deve l opmen t The procedure for ident i fy ing and assembl ing the equations has been demonstrated by example. The rema in ing quest ion that must st i l l be answered is what t ime step should be used for the analysis o f complex networks. The answer relates to how the nature o f the l inked boundary condi t ions change as p ipe replacement is considered. Th i s topic is presented i n the fo l l ow ing section. 104 5.4 N E T W O R K D I S C R E T I Z A T I O N : S E L E C T I N G T H E TTMF. S T E P Perhaps the greatest disadvantage to the direct appl icat ion o f the method o f characterist ics is the need to ma inta in a t ime step that is common to a l l parts o f the network. In order to calculate the condi t ions at any node, it is necessary to account for the inf luence o f the p ipes inc ident wi th each node at the speci f ic inc ident being considered. It is on ly once these condi t ions are known that the appropr iate boundary condi t ions equat ions can be developed and their so lut ion determined. The constant t ime step requ i rement immedia te ly raises two questions for the analyst: " W h a t must be done i n order to achieve an integer number o f reaches i n each p ipe?" and " H o w can the size o f the t ime step best be chosen?" These two questions fo rm the heart o f the network discret izat ion p rob lem; a reasonable and f lex ib le a lgor i thm for so lv ing it w i l l now be presented. T h e p rob l em o f ach iev ing an integer number o f reaches in each p ipe has received considerable attention i n the l i terature relat ing to hydrau l i c transient analysis. The reason this p rob l em arises is due to the fixed A x to A t rat io that is requ i red to app ly the me thod o f characteristics. Since, for a general system, the wavespeed and the length o f the var ious conduits i n the network varies widely, no single A t value w i l l usual ly be ab le to be chosen that achieves an integer number o f reaches in each p ipe. M a n y analysts (e.g., W y l i e and Streeter, Chaudhry , Watters, Sharp) feel that it is permiss ib le to a l l ow some smal l adjustment i n the value o f the wavespeed (say a m a x i m u m change o f ten or fifteen percent), in order to produce an integer number o f reaches i n each p ipe. Th i s procedure is jus t i f ied by the uncertainty i n the actual wavespeed value and by the relat ive insensit iv i ty o f the transient calculat ions to the wavespeed value. Wavespeed adjustments are preferable to changing the length o f any condu i t , 3 5 since the steady-state so lut ion is also altered by this ad jus tment (Steady-state 3 5 Th i s procedure is used wi thout comment by the program S U R N A L 105 condi t ions are unaffected by changes i n wavespeed.) Ano the r technique for app ly ing the me thod o f characterist ics to a p ipe wh i ch wou l d not no rma l l y have an integer number o f reaches is to develop an interpolat ion procedure. Such a technique a l lows the re lat ion between A t and A x to be re laxed so that the me thod o f characterist ics can be app l ied to the condui t us ing an unadjusted va lue o f wavespeed. T h e technique for app ly ing this procedure is i l lustrate i n F i gu re 15. T h e diagonal l ines i n this figure no longer intersect the gr id points chosen for the me thod o f characterist ics; however , the slope o f these l ines is chosen to satisfy the requ i red re lat ion between A t and A x. In order to calculate the unknown condi t ions at po int P, the head and discharge values must be known at some po in t a long both o f the d iagonal l ines. These condit ions are determined by interpolat ing between ne ighbor ing gr id points where condi t ions have been determined dur ing prev ious t ime steps. F o r a numbe r o f years, the usual procedure was to use l inear interpolat ion i n the x d i rec t ion at the prev ious t ime step (i.e., calculate the head and discharge at the to + 3 A t to + 2 A t to + A t P A / B \ S C R'I / R A' B' \ C' A" B" c" < A x » X F i gu r e 15. x - t P lane wi th Character ist ics G r i d and Interpolat ion 106 points R and S i n F igu re 15 f rom the k n o w n condit ions at points A , B and C) . Th i s procedure is easy to app ly and requires no addi t iona l computer storage over the regular method o f characteristics in order to imp lement it on the computer . However , this technique o f interpolat ing at a g iven t ime instant led to unwanted numerica l attentuation and dispers ion prob lems — p rob lems that have greatly restricted the usefulness o f this approach. The numer i ca l p rob lems have led some researches to develop new and more compl icated interpo lat ion procedures. Examp les o f such developments inc lude the " f i x ed - g r i d character ist ics" approach o f W igger t and Sundquist (1977), the " t ime - l i n e in terpo lat ion" procedure o f Go l dbe r g and W y l i e (1983) 3 6 and the ingenious var iab le t ime step method o f T r i k h a (1977). T o date, no systematic compar isons have been completed on actual hydrau l i c systems to determine the relative advantages and disadvantages o f the var ious interpo lat ion techniques or to compare the interpolat ion procedures to the other ways o f ach iev ing an integer number o f reaches i n each p ipe length. The more recent in terpo la t ion procedures invar iab ly result in more compl icated computer programs and larger memo r y requirements: it is on ly i f these disadvantages are more than offset by the imp roved transient analysis that the use o f the more complex techniques can be jus t i f i ed . Since a comprehens ive benef i t/cost analysis has not been per formed, a f lex ib le a lgor i thm to discreuze the network that is not commit ted to a part icu lar approach w i l l be presented here. The second part o f the discret izat ion quest ion raises the quest ion o f what t ime step shou ld be used to pe r f o rm the transient analysis. T he p rob l em o f rat ional ly selecting an appropr iate t ime step has rece ived relat ively l i tt le attention i n hydrau l i c transient appl icat ions. Th i s neglect is surpr is ing, g iven that the execut ion t ime for a " T i m e - l i n e interpolat ion can be s imply i l lustrated on F igu re 15. N o w the diagonal l ines are extended backward in t ime unt i l they have completed a crossing o f a fu l l A x step. These intersections are labe l led as R ' and S'. The unknown condit ions are then determined f r om the ne ighbor ing gr id po ints at that x pos i t ion. The reason for the increased storage requ i rement is immed ia te l y obv ious f r om this figure. 107 g iven p rob l em w i l l increase by approx imate ly four t imes every t ime the t ime step is ha lved. Ideal ly, an automatic procedure shou ld be developed wh i ch wou ld a l low an analyst to select the desired accuracy for a g iven appl icat ion and let the procedure choose an appropr iate t ime step for each part o f the network i n order to achieve this accuracy. Th i s t ime step cou ld then be adjusted as the solut ion develops i n order to ma in ta in this accuracy throughout the calculat ion. However , such a procedure cannot, at present, be deve loped since the errors associated wi th the method o f characteristics can not yet be s imp ly related to the t ime step values. W h e n l i nked boundary condit ions are present, the relat ion between the choice o f the t ime step and the p rogram execut ion t ime becomes even more compl icated. I f a subscr ipt 1 is used to ident i fy l i nked groups o f boundary condit ions, and E j is the number o f equations i n this group, then the approx imate re lat ion between execut ion t ime and the t ime step used for the analysis can be written T - ' • Execut ion T i m e = [ E C . N S . + E C , E , ] (10) A t ^ i i j l l i n wh i ch C j and C j are constants (general ly unknown) and N S j is the number o f sections i n p ipe i. In this expression, the f irst sum is also a funct ion o f the t ime step, since the number o f reaches i n any p ipe is approx imate ly equal to the wave travel t ime i n that p ipe d iv ided by the t ime step used for the analysis. Th i s re lat ionship assumes a l inear re lat ionship between execut ion t ime and the number o f sections i n each p ipe for a g iven t ime step. Th i s is a reasonable assumption because o f the expl ic i t nature o f the method o f characterist ics so lut ion procedure. However , execut ion t ime w i l l have an approx imate ly cub ic dependence on the size o f the l i n ked boundary cond i t ion sets ow ing to the d i f f i cu l ty o f so lv ing a set o f equat ions s imultaneously. A l t hough the constants C . and C } i n Equa t i on (10) can not be evaluated w i th certainty, two qual i tat ive conclusions can st i l l be d rawn f rom this relat ion. F i rst , f r om 108 the po in t o f v iew o f execut ion t ime, large t ime steps are almost always desirable. O f course, us ing too large a t ime step w i l l misrepresent the transient i n some networks, so a m a x i m u m t ime step must usual ly be imposed by the ana lys t Second, the mot iva t ion for us ing p ipe replacement elements is der ived f rom the larger time step these elements make possible. However , p ipe replacement increases the complex i ty o f the network wh i ch tends to offset this advantage. Thus , p ipe replacement is only jus t i f i ed where the increase i n effort devoted to so lv ing the l i nked sets is more than offset by the decrease i n t ime devoted to comput ing the solut ion at the characterist ic sections i n the network pipes. N o w that some pre l im inary discussion re lat ing to the discret izat ion quest ion has been presented, the procedure recommended for choosing the time step can be descr ibed. 5.4.1 T H F D I S C R E T I Z A T I O N A L G O R I T H M The fo l l ow ing procedure for selecting the t ime step is in tended p r imar i l y as a useful tool for invest igat ing network response. M o r e work is needed i n the approx imat ion theory o f hydrau l i c networks i f a general and w ide ly appl icable procedure is to be developed. The a lgor i thm is based on a pr ior i ty system that considers some alternatives for ach iev ing an integer number o f reaches super ior to others. The procedure has p roved useful for analyz ing specif ic prob lems by a l low ing a good representat ion o f the network to be achieved by s imp ly a l ter ing a number o f input parameters. It is hoped that the a lgor i thm w i l l also a id in general network invest igat ion by p rov id ing in format ion on wh i ch adjustment or interpolat ion techniques work best under what circumstances. The a lgor i thm requires on ly four pieces o f input data: the m a x i m u m and m i n i m u m a l lowable time step, a wavespeed adjustment parameter and a p ipe replacement parameter. The m a x i m u m t ime step is selected by the analyst as the 109 largest t ime step value he feels can be used to accurately represent the transient condi t ions i n the network. The m i n i m u m t ime step is the smal lest t ime step the analyst feels cou ld possib ly be just i f ied on economic grounds for mode l ing the network. Th i s procedure di f fers f r om most other programs wh i ch require the analyst to p rov ide a single t ime step. The proposed approach is more realist ic in that the analyst must prov ide a range o f t ime step values, determined f rom experience, but a l lows the p rogram to make the final select ion depending on the detai ls o f the computer representat ion. The wavespeed adjustment and P R E parameters are used for two purposes: to contro l the a l lowable change in wavespeed ( • ) and to contro l wh i ch o f the two p ipe replacement elements w i l l be chosen by the p rogram ( N ). The details o f how these values w i l l be app l ied w i l l be discussed be low. The steps requ i red to choose the t ime step for the analysis and to achieve an integer number o f reaches i n each p ipe can now be i temized. 1. Set A t to be equal to the m a x i m u m avai lable t ime step. 2. U s i n g this t ime step, find the integer number o f reaches each p ipe. The number o f reaches, i n turn, determines the ' requ i red ' wavespeed adjustment (i.e., the adjustment i n wavespeed necessary to achieve the number o f reaches just selected). 3. F o r a l l p ipes i n wh i ch the requ i red wavespeed adjustment is less than the a l lowable adjustment ( * ), adjust the wavespeed by the requ i red amount in this p ipe. I f a l l p ipes in the network are inc luded i n this group, stop. 4. I f the requ i red wavespeed adjustment is too large, check the connect iv ity o f the o f fend ing pipes. Tha t is, for those pipes w i th a large required wavespeed adjustment, p ipe replacement w i l l be attemped. I f replac ing a p ipe causes a l i nk between two (otherwise unconnected) l i nked groups o f boundary condit ions, p ipe replacement w i l l greatly increase the work o f 110 so lv ing these l i nked sets (the second te rm i n Equat ion 10). Th i s increase is undesirable and p ipe replacement shou ld be ru led out i f tHis occurs. However , i f either o f the nodes connected to this p ipe do not belong to a l i nked group, replace this p ipe and update the l i nked node data structure. Rep lacement w i l l take place as fo l lows: i f the p ipe for wh i ch the P R E is be ing substituted wou ld have had less than N reaches, the l umped pre inert ia parameter w i l l be used. Howeve r , i f this p ipe wou l d have had the same or a greater number o f reaches than N , the Finite dif ference pre approx imat ion for this p ipe w i l l be used. I f no p ipe replacement is desired, N shou ld be set to zero, pre I f a l l p ipes i n the network now have an integer number o f reaches or have been replaced by PRE ' s , stop. However , i f some pipes st i l l have large requ i red wavespeed adjustments and these pipes connect sets o f l i nked boundary condit ions, a reduct ion o f the t ime step shou ld be investigated. Tha t is, decreasing the t ime step (whi le ma inta in ing A t to be greater than or equal to the m i n i m u m t ime step) may result i n the requ i red wavespeed adjustment i n these p ipes becoming 'a l lowable. ' I f a l l o f the p rob l em pipes can be treated i n this manner, adopt the new t ime step value. I f the t ime step is changed, a l l data shou ld be reset to its or ig ina l values and the p rogram shou ld return to Step 2. I f noth ing has been complete ly successful so far, use l inear interpolat ion in the x d i rect ion at the prev ious (known) t ime step. I f care is taken to satisfy the Couran t stabi l i ty requ i rement (i.e., i n F i gu re 15 po in t R shou ld l ie between A and C and po in t S between C and B) , the interpolat ion procedure w i l l a lways a l low a so lut ion to be obta ined at the m a x i m u m t ime step. I l l The above procedure a l lows a f lex ib le • range o f alternatives to be investigated. Obv ious ly , i f a single t ime step is requ i red, the analyst can s imply equate the m i n i m u m and max imum t ime step to this value. Choos i ng a large va lue for the wavespeed adjustment parameter w i l l force an integer number o f reaches on every p ipe and no pipe replacement w i l l be requi red. P ipe replacement can take place in two ways: by forc ing any p ipe to be a P R E in the input da ta 3 7 or by a l lowing the program to replace some pipes dur ing the operat ion o f the discret izat ion rout ine. I f no p ipe replacement is required, the N p r e value can be set to zero. F ina l l y , shou ld interpolat ion be the desired alternative for a l l pipes, the value o f the wavespeed adjustment parameter and the p ipe replacement parameter can both be set to zero. The f lex ib i l i ty o f this a lgor i thm shou ld greatly faci l itate a general study o f the advantages and disadvantages o f the various approaches to the discret izat ion p rob l em and shou ld be o f great benef i t to analyz ing ind iv idua l networks. The f ina l level o f discret izat ion produced by the above a lgor i thm wi l l be a compromise between desired p rogram eff ic iency and the requ i red accuracy o f the transient analysis. P ipe replacement and interpolat ion procedures both a l low a larger time step to be used than wou ld otherwise be possible; however, these techniques also introduce unwanted numer ica l attenuation and dispers ion into the ca lcu lat ion o f transient condit ions. W a v e speed adjustment, on the other hand, accurately mainta ins sharp wave fronts but alters both the magni tude o f the transient pressures and the wave travel t ime between var ious network locations. The relat ive meri ts o f these approaches w i l l depend on the nature o f the transient be ing analyzed, the characteristics o f the network and the purpose o f the transient analysis. The ul t imate decis ion o f wh i ch procedure (or combinat ion 3 7 I n the program, this opt ion is selected by s imp ly setting the wavespeed for this part icu lar p ipe to zero. The p rogram wi l l then recognize this p ipe as a P R E 112 o f procedures) is most appropr iate w i l l be an economic one. 5.5 CONCLUSIONS I f an eff ic ient p rog ram o f ana lyz ing networks in the transient state is to be developed, cons iderat ion must be g iven not on ly to der iv ing the solut ions themselves, but also to the quest ion o f how these solut ions can be effect ively imp lemented on the computer. In the past, very l i tt le attention has been devoted to this p rob lem as it relates to the study o f hydrau l i c transients. Th i s chapter has sought to look at the p rob l em o f computer imp lementa t ion f rom three d i f ferent points o f v iew: representing the network topology, ident i fy ing and so lv ing l i nked groups o f boundary condit ions and discret iz ing the network. 1. T h e first section o f this chapter has developed the data structures and procedures requ i red to describe how the var ious hydrau l i c elements found in the network are connected to each other. Da ta structures for ident i fy ing wh i ch pipes meet at each node are developed f r om a m i n i m u m o f input data. In add i t ion, a storage and retr ival system for the sections requ i red to apply the method o f characterist ics to the network is descr ibed. The technique uses a one -d imens i ona l array i n order to conserve computer memory . 2. Nex t , the p rob lems associated wi th l i nked boundary condit ions are cons idered. W h e n boundary condit ions become l i nked , the fundamenta l boundary cond i t ion equat ions do not change, a l though the f o rm o f the characterist ic equations app l ied at each node is a l t e red A procedure for indent i fy ing the l i nked elements and the data structure requ i red to imp lement this procedure is descr ibed. A n example is used to i l lustrate the appl icat ion o f these techniques to a smal l group o f l i nked boundary condit ions. 3. F i na l l y , the p r ob l em o f network discret izat ion is discussed. The quest ion o f what t ime step to use for the analysis depends o f the specif ics o f the network be ing 113 studied, the purpose o f the analysis and the accuracy o f the avai lable input data. T o a id i n the invest igat ion o f this important quest ion, a procedure is descr ibed wh i ch prov ides a pr ior i ty rank ing o f alternatives. The a lgor i thm first considers wavespeed adjustment fo l l owed by p ipe replacement, t ime step adjustment and, finally, l inear interpolat ion i n order to arr ive at the final level o f discretizat ion for the network. The so lut ion procedures and imp lementat ion details requ i red to analyze large d ist r ibut ion systems in the transient state are now in place. The method o f characterist ics has been app l ied to the network and careful cons iderat ion has been g iven to the network boundary condit ions. The quest ion that must now be answered is "does the procedure work?" The quest ion o f how to ver i fy the network mode l is discussed i n the next chapter. 6. V E R I F I C A T I O N O F T H E S O L U T I O N P R O C F D I T R F In general , the object o f a numer ica l technique for so lv ing a mathemat ica l p r ob l em is "to determine suff ic ient ly accurate approx imat ions w i th m in ima l effort" (Bu rden et al., 1981, p.194). In this thesis, the p rob l em o f pred ic t ing transient f lows and pressures i n large d ist r ibut ion systems has been considered. The method o f characterist ics has been app l ied i n order to obtain solut ions for the head and discharge at the inter ior sections o f each p ipe, and i n order to a l low the aux i l ia ry condit ions de f in ing the network boundary condi t ions to be determined. Procedures have been suggested for convenient ly represent ing the topology o f the network and for storing i n the computer the computat ional sections used to app ly the method o f characteristics. The subject that must now be cons idered is whether the a lgor i thms are correct and whether they are suff ic ient ly accurate to warrant appl icat ion to real networks. Essent ia l ly, three steps are requ i red i n order to establish the correctness o f a computer a l go r i t hm. 3 8 The first step is to def ine what "correct" means. Genera l l y , a procedure is sought wh i ch w i l l produce the right answer i n a reasonable amount o f computer time. However , since the method o f characterist ics is an approx imate solut ion to an approx imate set o f equations, def in ing what the right answer is requires some considerat ion. The computer a lgor i thm itsel f has two components: the solut ion procedure that i t is based on and the sequence o f computer instruct ions wh i ch are used to carry out this procedure. Obv ious ly , for an a lgor i thm to be correct, both the procedure and the instruct ions must themselves be cor rec t Steps 2 and 3, therefore, consist o f prov ing the correctness o f f irst the procedure and then the computer representat ion o f the procedure. Th i s chapter w i l l be loosely organized a long the l ines o f these three steps. 3 "These three steps are i temized i n Computer Algorithms by S. Baase (1978). 114 115 6.1 V E R I F I C A T I O N C R I T E R I A F o r s imple p ipe l ine systems consist ing o f s ingle pipes, or a number o f p ipes connected together i n series, the va l id i ty o f the method o f characterist ics has been f i rmly establ ished. The general . agreement between theoret ical and exper imenta l results has often been excel lent and is reported in many textbooks and i n numerous research papers. Indeed, the general procedures are so we l l establ ished that Va rdy (1981) recommends using the method o f characteristics as one o f the best ways o f establ ishing the general va l id i ty o f any new approach to hydrau l i c transients. T o quote, the va l id i ty o f any new theoret ical approach can be evaluated relat ively cheaply by compar ing its predict ions w i th those o f the M O C [method o f characterist ics]. Large differences indicate errors i n the new approach; smal l d i f ferences represent open questions that must be answered by precise experiments. Thus , for s imple p ipe l ines at least, i t can safely be assumed that the method o f characterist ics is a correct and va l id approach for pred ict ing the transient response o f hydrau l i c systems. The quest ion that now arises is whether the method o f characterist ics is a va l i d procedure for analyz ing transient f lows i n large d istr ibut ion systems and i n what way its correctness can be proved. The approach taken i n this thesis w i l l be to construct a series o f test networks wh i ch have unusual symmetry properites: the symmetry o f these networks w i l l be used to establ ish an equiva lence between the network and a s imple p ipe l ine system. I f the transient response o f the network is a lmost ident ica l w i th the response o f the s imple p ipe l ine when subjected to the same boundary condit ions, the val id i ty o f the network approach is ver i f i ed . Thus , i n summary, the p r imary assumpt ion that w i l l be made in order to ver i fy the correctness o f the a lgor i thms presented i n this thesis is that the method o f characterist ics is a va l id numer ica l procedure for calculat ing transient f lows in s imple pipel ines. There are two impl icat ions wh i ch immedia te ly arise f r om this assumption. F i rs t , the general boundary cond i t ion formulas and equations presented i n Chapter 4 116 must reduce to their s imple p ipe l ine equivalents when app l ied to s imple systems. Second, any network that can be made equiva lent to a s imp le p ipe l ine must respond in an a lmost ident ica l manne r 3 9 to that s imple p ipe l ine. H o w these two impl icat ions can be used to ver i fy the general network a lgor i thm and program developed in this thesis is the subject o f the remainder o f this chapter. 6.2 B O U N D A R Y C O N D I T I O N F O R M U L A T I O N The prev ious discussion has completed the first step o f the ver i f icat ion procedure by def in ing correctness in terms o f agreement wi th establ ished method o f characteristics procedures for s imple p ipel ines. The second step involves p rov ing the correctness o f the actual so lut ion procedures. S ince the network solut ion and the s imple p ipe l ine procedure both emp loy the method o f characteristics, the ver i f icat ion o f the ca lculat ion procedure for internal sections is automatic. The on ly d i f f i cu l ty arises i n ver i fy ing the boundary cond i t ion equat ions deve loped i n Chapter 4. Because the number o f boundary condi t ions that have been presented is qui te large, on ly one example w i l l be presented. However , the same approach can be used to ver i fy the correctness o f a l l o f the boundary cond i ton equat ions presented i n this thesis. 6.2.1 E X A M P L E : V A L V F T N - I J N F . B O U N D A R Y C O N D I T I O N A very wel l establ ished boundary cond i t ion occurs when a valve is insta l led between the abutt ing ends o f two pipes. The typical arragement is " E x a c t argeement between the two systems is not possible for two reasons: inacuracies i n the steady-state so lut ion and dif ferent effects o f computer r o u n d - o f f i n the two systems. Tha t is, since the steady-state so lut ion in the network must be establ ished by an i terat ive procedure, a final exact flow balance w i l l , i n general, not be possib le; however , the steady-state condi t ions in a s imple p ipe l ine can be exactly establ ished. In add i t ion , s ince the number o f error p roduc ing operat ions is greater in a network than i n a s imp le p ipe l ine , i t must be expected that the effect o f f i n i t e -d ig i t ar i thmet ic is more impor tant i n the network solut ion. I f the numer ica l procedure is stable, neither o f these two effects shou ld be impo r t an t 117 shown i n F igu re 16. The general network so lut ion to this boundary cond i t ion has been deve loped in Chapte r 4 and the so lut ion at that t ime was found to be Q e x t = ~ b + + 3°' 5 ( 1 ) . i n wh i ch b = | ( r E s ) 2 B (2) and c = s( r E s ) 2 C (3) In add i t ion , s was shown to take the sign o f C . In order to s imp l i f y these equations, i t is noted that both nodes 1 and 2 in F igure 16 have a degree o f A A one. The means that B = (Bi + B 2) and C = ( C p - C M ) . In addi t ion, i f we adopt the termino logy used i n W y l i e and Streeter, ( r E g ) 2 = 2 C y . In effect, then, b = s C y ( B 1 + B 2) (4) and F i gu re 16. Va l ve I n - L i n e 118 c = 2 s C y ( C p - C M ) : (5) F o r pos i t ive f low, ( C p - C j ^ ) is greater than zero and s = 1. Th i s produces an equat ion for the f low through the valve ident ical to Equat ion 3 -43 on page 49 o f fluid Transients by W y l i e and Streeter. O n the other hand, i f ( C p - C j ^ ) is negative, s = - 1 and Equat ion (1) above becomes ident ical to Equat ion 3 -44 i n Fluid Transients. Thus , Equat ion (1) is equivalent to the s imple p ipe l ine equations i n this special case. A l t hough tests o f this k i nd are not, i n themselves, conclusive, this k i nd o f procedure lends strong support to the val id i ty o f the general boundary cond i t ion expressions. A l l o f the equations i n Chapter 4 have been tested against the equiva lent s ingle p ipe l ine expressions, where such equations are avai lable, and have been found to agree i n every case. It is therefore considered that the general boundary cond i t ion formulat ions g iven i n Chapte r 4 are consistent w i th the s imple p ipe l ine appl icat ion o f the method o f characteristics. 6.3 NETWORKS WITH SPECIAL SYMMETRY The procedures for ver i fy ing ind iv idua l expressions, such as those presented above, are reassuring but, obvious ly , do not prov ide a rigorous test o f the network a lgor i thm. A more power fu l approach wou l d be to write a general network program and compare the results i t produces w i th some exact so lut ion for the same p rob lem. However , since no exact solut ions exist for even s imple p ipe l ines when fr ic t ion is inc luded, and since no so lut ion o f the network p rob l em is possible wi thout f r i c t i on , 4 0 a mod i f i ed procedure is ca l led for. " T h e steady-state d istr ibut ion o f f lows in a network is indeterminant when fr ic t ion is absent 119 •The assumpt ion o f the general va l id i ty o f the method o f characterist ics for s imple p ipe l ine systems al lows a more pract ical ver i f icat ion procedure to be suggested. I f networks can be contr ived wh i ch are 'equivalent ' to s imple p ipe l ine systems, the transient response o f the two systems can be compared. I f the network so lut ion and the s imple p ipe l ine so lut ion agree, the general val id i ty o f the procedure is establ ished. One part icu lar network wh ich is equiva lent to a s imple p ipe l ine system wi l l next be descr ibed, fo l l owed by a discussion o f a number o f s imple numer ica l exper iments pe r fo rmed on this system. The so lut ion to the equivalent s imple p ipe l ine system wi l l be obta ined by us ing a method o f characterist ics program taken f rom a wel l known reference text on hydrau l i c transients. 4 1 6.3.1 DTSCRIPTION OF THF. TEST NETWORK The ver i f icat ion procedure adopted i n this thesis requires a network to be found wh i ch is equivalent to a single p ipe l ine system. I f no such networks existed, ver i fy ing a general network p rogram wou ld be a much more d i f f i cu l t task. For tunate ly , however, a large number o f such systems can be developed. F i gu re 17 i l lustrates a network wh i ch is equivalent to a s imple p ipe l ine system. The network comprises 80 pipes, 52 nodes, two in f low reservoirs and two valves where outf low can occur. T h e length o f each p ipe, its diameter, wavespeed and f r i c t ion are a l l ident ica l . F o r c lar i ty i n the i l lustrat ion, the node and p ipe numbers are not shown. In practice, the p ipe and node number ing w i l l usual ly be sequent ia l , start ing at one and cont inu ing through each p ipe and node in the network. Howeve r , when this network was tested, a number o f node and p ipe numbers i n the sequence were intent ional ly omit ted i n order to ver i fy that the 4 1 T h e spec i f ic p rog ram used is based on the va lve -c losure p rogram found on pages 469 -473 o f Applied Hydraulic Transients by H . Chaudh r y (1979). T he p rogram was s l ight ly mod i f i ed to a l low th ree -d imens iona l representations o f the pred icted response to be p roduced . 120 O 6 0-^ — 9 — 9 — 1 o-o — Q 0 O Q-SUMMARY OF NETWORK 80 PIPES 52 NODES 29 LOOPS INITIALLY Q= 1.0 FT 3/SEC FOR ALL PIPES L = 600 FT D = 0.5 FT a = 2400 FT/SEC f = 0.020 F igu re 17. Sketch o f Test Ne twork procedure does not requ i re sequent ia l number ing . O f the 52 nodes i n the network, 24 are o f degree 2 and 28 are o f degree 4. The locations o f the in f lows and outf lows for this system are indicated i n F igu re 17. The in f lows represent two constant head reservoirs wh ich ma inta in the hydrau l i c gradel ine elevat ion o f these nodes at 500 fee t The outf lows represent valves d ischarging to the atmosphere. The specif ic values chosen for the in i t ia l discharges, the reservoir e levat ions and the values o f the var ious p ipe constants are arbitrary, a l though the wavespeed values are chosen to produce a integer numbe r o f reaches i n each p ipe. T he one essential restr ict ion is that the propert ies o f both the two reservoirs and the two valves must be ident ica l . That is, the valves must be operated i n exactly the same manner and the hydrau l i c gradel ine elevations o f both reservoirs must be ident ica l . 121 The p r imary characterist ic o f the test network is its symmetry. Because the propert ies o f a l l p ipes and devices i n the network are ident ica l , i f the network were to be spl i t between the two reservoirs and the two valves a long a l ine at 45 degrees to the hor izonta l , two ident ica l smal ler networks wou ld be fo rmed. Bo th o f these smal ler networks have an axis o f symmetry at 45 degrees and can themselves be spl i t into two s imple p ipe l ines hav ing ident ical propert ies. It is this bi lateral symmetry o f the network wh i ch establishes the equivalence between it and a s imple p ipe l ine system. The argument for prov ing the equiva lence o f the test network to a s imple p ipe l ine can be better exp la ined by referr ing to F igu re 18. In this f igure, one o f the nodes o f degree four a long the the pr inc ipa l axis o f symmetry for the text network is shown. The p ipes i n F i gu re 18 have been numbered for convenience and the axis o f symmetry has been shown as a dashed l ine. Because o f the symmetry in both the p ipes and the boundary condit ions, the f low in p ipes 1 and 4 must a lways be equal and the f low in p ipes 2 and 3 must also be equal . But, because o f the cont inu i ty equat ion at the node, the f l ow in a l l o f the p ipes must therefore be equal . Thus , the boundary cond i t ion at this node F i gu re 18. N o d e o f Degree 4 a long the P r ima ry A x i s o f Symmetry 122 is now indist inguishable f r om two junct ions o f degree two invo l v ing p ipes (1, 2) at the first junc t ion and (3, 4) at the second. The same arguement can be extended to app ly -to a l l o f the nodes i n the network. F o r the nodes o f degree 2, the argument is even s impler since cont inuity i tself proves that the in f low and out f low discharges must be equal. S ince the symmetry and cont inu i ty arguments ho l d at a l l nodes, the test network shown in F igu re 17 is ident ical to 4 s i d e - b y - s i d e s imple p ipe l ines hav ing no interact ion w i th one another. Thus , i n order to test the network i n F igure 17, it is on ly necessary to consider a s imple p ipe l ine consist ing o f twenty pipes connected i n series and d icharg ing one ha l f o f the steady-state discharge Q . It shou ld be noted that the argument does not depend i n any way on the number o f p ipes i n the network or on the number o f external discharges. A n y network hav ing s imi l iar symmetry propert ies wh ich a l lows it to be spl i t into ident ical single p ipes carry ing equal discharges at each node is also equiv lent to a single p ipe. F o r example, i f every p ipe i n the test network were to be doub led or t r ip led, and the steady-state discharges increased by an equiva lent amount, ident ical response wou ld be obta ined. The author has tr ied the test network w i th tr ip le the number o f pipes, and three times the value for the external flows, and found the single p ipe response was accurately repeated, a lbeit at a much greater computer cost. Thus , the degree o f the nodes i n the test •network is arb i t rary. A few typical examples o f actual transient situations w i l l now be related. 6.3.2 E X A M P L E 1: S I I D P F N V A L V E C L O S U R E The first transient condi t ion to be considered occurs when the two valves at the downstream end o f the test network are suddenly c losed. The discharge at the downstream end o f a l l four l ines i n the network is suddenly arrested by 123 this act ion, produc ing a shock wave. The shock wave propagates a long each p ipe i n the network at the prescr ibed value o f the water hammer wave velocity. D u r i n g the passage o f this wave through the network, no wave ref lect ions can occur at the p ipe junct ions, since the propert ies o f a l l p ipes in the network are ident ica l . However , when the wave arr ives at the reservoir, the constant head boundary cond i t ion is imposed and the p r imary shock wave w i l l be ref lected. The ref lected wave w i l l , i n turn, be i tsel f ref lected by the now closed discharge valve. Th i s process o f successive wave ref lect ions w i l l cont inue unt i l f r ic t ion establishes the new steady-state cond i t ion i n the network ( in this case, the new steady-state is zero f low in a l l p ipes and the reservoir head at a l l nodes). Exact ly the same sequence o f events occurs i n the equiva lent s imple p ipe system when subjected to the same boundary condit ions. In the s imple p ipe l ine case, the output o f the va lve closure program is used to generate a th ree -d imens iona l representat ion o f the transient pressure head response. The result is shown in F igu re 19. In the th ree -d imens iona l image, the valve is located nearest the v iewer at the 12,000 ft locat ion on the distance axis. The constant head reservoir is found at the opposite end o f the distance axis, furthest f r om the viewer. The th ree -d imens iona l image c lear ly i l lustrates the process o f wave propagat ion and ref lect ion i n the p ipe l ine . A t t imes, the th ree -d imens iona l transient t ime history p lo t is d i f f i cu l t to interpret i n detai l . F o r this reason a number o f pro f i les or sections through this surface are generated at the same t ime as the th ree -d imens i ona l p l o t F i gu re 20 i l lustrates the t ime - head history for a number o f key locations a long the p ipe l ine. T h e slope o f the head pro f i l e after the passage o f the wave f ront can c lear ly be seen i n F igure 20. T h e slope o f this l ine results f r om the head losses due to f r ic t ion i n the p ipe l ine. 124 F i gu re 19. Trans ient Response: Sudden Va lve C losure O n e o f the most st r ik ing features o f F igure 19 is the abruptness o f the shock wave p roduced by the sudden valve closure. The method o f characterist ics reproduces these sharp, wave fronts very accurately. Th i s is in contrast to f in i te d i f ference techniques wh i ch tend to attenuate or smooth the discont inuit ies w i th t ime (Bu rden et al. , 1981). In a real ist ic network situation, consist ing o f a large numbe r o f p ipes hav ing di f ferent propert ies, the progressive effect o f this smooth ing cou ld greatly impa i r the accuracy and value o f the f inite dif ference approach. F o r this reason, the method o f characteristics is prefer red to the f in ite 125 10.0 TIME 15.0 LEGEND • Q. POINT o M. POINT A 3Q. POINT v VALVE END 20.0 F igu re 20. T i m e - H e a d Prof i les: Sudden Va lve C losure di f ference technique o f Sevak (1980) when ana lyz ing rap id ly vary ing unsteady f lows i n networks. W h e n the output o f the network p rogram deve loped for this thesis and the s imp le p ipe l ine p rogram o f Cha udh r y were compared, very good agreement was observed. U s i n g E p p and Fow le r ' s steady-state p rogram as a subroutine, the network p rogram converged to the expected steady-state d ist r ibut ion o f f lows i n seven iterations. T he m a x i m u m error i n the steady-state flow was less than 0.000005 f t 3 /sec. Next , the transient so lut ion was compared. The transient analysis for both the network and the s imple p ipe l ine used four reaches i n each p ipe; this leve l o f disret izat ion is equiva lent to a time step o f 0.0625 seconds. The s imple p ipe l ine and the network so lut ion were compared at each node for the first 20 seconds o f s imulat ion on a one second interva l . In every case, the network so lut ion agreed w i th the s imple p ipe l ine so lut ion to w i th in 0.01 f t 126 The agreement between the s imple p ipe l ine pred ict ions o f the transient pressures and f lows and the output o f the network p rog ram demonstrates ..the correctness o f both the network p rog ram and the fo rmu la t ion o f the network a lgor i thms presented in this thesis. O f course, only two very s imple boundary condi t ions are ver i f ied by this test However , the va l id i ty o f the boundary cond i t ion equations can be ver i f i ed independent ly o f the network p rogram along the l ines suggested earl ier i n the chapter. In order to prove that the agreement between the s imple p ipe l ine and the network was not s imply fortuitous, a number o f addi t iona l numer ica l exper iments have been per fo rmed. The first and simplest o f these tests was to tr ip le al l o f the pipes i n the network as we l l as the va lue o f the steady-state discharge. The execut ion t ime was now approx imate ly three t imes greater, but the transient response was v i r tua l ly ident ica l . T he next tests i nvo l ved d iameter changes i n the network p ipes and more compl i ca ted valve manouvers. T w o o f these tests w i l l now be br ie f ly descr ibed. 6.3.3 E X A M P L E 7: L I N E A R V A L V E C L O S U R E W I T H D I A M E T E R C H A N G E S The symmetry arguments that were used to establ ish the equivalence between the test network and the s imple p ipe l ine do not, i n fact, require that every p ipe i n the network be ident ica l . T he essential e lement i n the argument is that the network must be symmetr i c about a bisecting l ine d rawn between the reservoirs and the valves. Thus , i f the network were to be b roken into bands, such as those shown i n F i gu re 21, and the propert ies o f the p ipes changed f r om one band to the next, the network wou ld remain equ iva lent to a s imple p ipe l ine system. F i gu re 22 shows the transient response o f the s imp le p ipe l ine system wh i ch is equiva lent to a banded network. In this example, the va lve is closed B o — 6 — — 6 — 6 o 6 6 6—-6 -O 6 9 9 9 (? Q 6 6 6 6 g 6—-6-—6 6 0 6 6 6 6 Q — - 6 6 6 6 Q 6 6 6 6 9 9 9 9 6 F igure 21. Test Ne two rk w i th Banded Structure 127 l inear ly i n 5 seconds. The only property a l lowed to vary between the bands is the p ipe diameter. D iamete r valves o f 0.75, 0.50, 0.50 and 0.40 ft are used in zones A , B, C and D respectively. The reduct ion in magni tude o f the p r imary pressure wave as it enters the por t ion o f the p ipe l ine hav ing the largest d iameter is apparent i n F igu re 22. Ano the r interesting feature o f this system is that the pressure head peak at the valve actual ly becomes steeper for the second posi t ive wave that it was for the first. Th i s unexpected response is a result o f ref lect ion f r om the d iameter changes a long the length o f the condu i t A l t h ough the transient response for this second example network is cons iderably more compl icated than that o f the prev ious network descr ibed, the transient pressure head pred icted by the network p rogram again agreed i n every impor tant detai l w i th the s imp le p ipe l ine system. The m a x i m u m observed 128 F i gu re 22. Trans ient Response: L inea r Va l ve C losu re discrepancy between the two programs was less than 0.02 ft. Thus, the correctness o f the network p rogram is again ver i f i ed by the s imp le p ipe l ine response. 6.3.4 E X A M P L E 3: C O M P L E X V A L V E M O V E M E N T Once the basic test structure has been developed, the number o f i nd i v i dua l tests that can be pe r fo rmed on the test network is un l im i ted . N o purpose is served by conduct ing endless var iat ions on this k i nd o f network test T h e author has tested e ight d i f ferent valve closure arrangements w i th var ious 129 d iameter changes i n the network bands. In a l l cases, agreement between the s imple p ipe l ine system and the network mode l has been excel lent. T he max imum recorded error occurred when a very comp lex valve movement was tested. Th i s case invo lved a smal l instantaneous valve closure at t ime zero, fo l lowed by a l inear open ing to the fu l ly open pos i t ion i n two seconds. The l inear open ing was then immediate ly fo l lowed by a l inear c losure i n three seconds. Th i s k i nd o f complex valve movement produces a very compl icated transient response. The output f r om the th ree-d imens iona l p lot t ing rout ine is i l lustrated in F igure 23. In this case, the ouput f rom the two mode ls d i f fe r red by as much as 0.10 ft for some o f the compared points. Th i s agreement is st i l l far better than that wh ich cou ld be ach ieved using real data, and so is not considered a serious shortcoming. In general, the correctness o f the network p rogram is strongly supported by the numer ica l tests per fo rmed on s imple p ipe l ine systems. 6.4 C O N V E R G E N C E T O S T E A D Y - S T A T E A l l o f the transient tests that have been pe r fo rmed to date have ver i f i ed the general va l id i ty o f the network p rogram developed for this thesis. However , the networks tested have a l l had unusual symmetry propert ies, very un l i ke networks found in actual d is t r ibut ion systems. The quest ion that arises is whether there is a way o f testing the network p rogram on a system hav ing arbi trary geometry. O f course, the best test wou ld be to predict the transient response and behav ior o f an actual hydrau l ic system, then measure that response and compare the results. However , for the purposes o f the present study, this approach was considered too expensive and t ime -consum ing . A much s imp ler and less costly procedure is to compute the transient response for a very long s imulat ion t ime. Since the steady-state equations are s imply special cases o f the transient equations, the transient heads and f lows should eventual ly stabi l ize at the correct new steady-state cond i t ion i n the network. In general, it is 130 poss ib le to consider transient phenomena in hydrau l i c systems as s imp ly the mechan ism for af fect ing changes between steady-state condit ions. A l t hough such an approach certa in ly does not verify the correctness o f the intervening states, convergence to steady-state is cons idered a strong test for the correctness o f any transient program (Wy l i e and Streeter, 1983). The one addi t ional advantage o f this steady-state test is that it w i l l c lear ly demonstrate whether or not the solut ion procedure is stable (an unstable procedure wou ld u l t imate ly diverge f r om the f ina l steady-state condi t ion) . F i gu re 23. Trans ient Response: Comp l e x Va l ve M o v e m e n t 131 The network p rog ram has been tested on a number o f occasions for convergence to steady-state. In each case, the new steady-state cond i t ion predicted by the transient network p rogram was ver i f ied by either logic (convergence to zero) or by the steady-state network p rogram o f E p p and Fow le r . Runn i ng a network p rogram to convergence requires a large amount o f computer t ime, but, i f the program does converge to the conec t steady-state condi t ion, the val id i ty o f the solut ion procedure is further ver i f ied . Th i s is another reason why it is so desirable to directly couple a power fu l steady-state p rogram wi th a general transient a lgor i thm. 6.5 C O N C L U S I O N S The goal o f this thesis is to produce eff ic ient and general procedures for ana lyz ing transient f lows i n hydrau l i c networks. Ear l i e r chapters have formulated the general so lut ion procedure, discussed the p rob lem o f obta in ing steady-state condit ions in a network and deve loped the network boundary cond i t ion equations. In add i t ion, the un ique prob lems associated wi th eff ic ient representat ion o f large networks have been discussed. The purpose o f the present chapter has been to ver i fy the correctness o f the prev ious developments. T o this end, a general network p rogram has been deve loped and tested. In order to numer ica l ly ver i fy this program, the assumpt ion has been made that the me thod o f characterist ics is a va l id and correct procedure when appl ied to s imp le p ipe l ine systems; an assumpt ion amp ly ver i f i ed by exper imenta l evidence. The ver i f i cat ion procedure consists o f deve lop ing a test network wh i ch is equivalent to a s imp le p ipe l ine system. The s imple p ipe l ine is then analyzed by establ ished method o f characterist ics procedures wh i l e the equivalent network is ana lyzed by the general network program. Since the transient response pred icted by the network program has i n every case tested agreed w i th the s imple p ipe l ine predict ions, the general va l id i ty o f the network p rogram has been establ ished. 7. C O N C L U S I O N S Th i s thesis has attempted to reformulate the procedures for pred ict ing the transient response o f large hydrau l i c systems wi th part icular reference to water d is t r ibut ion networks. Trans ient condit ions are both in it iated and mod i f i ed by the operat ing and contro l devices i n the network. F o r this reason, part icu lar attention has been pa id to the general deve lopment and fo rmula t ion o f the network boundary condit ions. Emphas i s has also been placed on deve lop ing data structures and procedures wh i ch a l low the response o f large networks to be convenient ly and eff ic ient ly obta ined. The a i m has been to develop a general , eff ic ient and re l iab le a lgor i thm for comput ing the transient response o f large d is t r ibut ion systems. A l t hough some further ref inement and ver i f i cat ion work is st i l l requ i red, it is considered that this goal has been achieved. S ign i f i candy, the a lgor i thm is general , since a wider class o f networks can be ana lyzed than was prev ious ly possib le; the p rogram is eff ic ient, s ince the solut ion o f many o f the network boundary condi t ions is made exp l ic i t for the first t ime and because computer storage is conserved i n the p rogram imp lementat ion; and, finally, the p rog ram is re l iable, since a number o f numer ica l exper iments have p roven its correctness. Spec i f i ca l ly , the so lut ion procedure deve loped i n this dissertation for obta in ing the transient response o f large networks inc ludes the fo l lowing features. 1. The unsteady or t ime-dependen t por t ion o f the p rogram has been successful ly l i nked w i th a very eff ic ient steady-state network p rogram (based on the work o f E p p and Fow le r ) . T he coup l ing o f the steady and unsteady parts o f the analysis improves p rog ram eff ic iency by avo id ing dup l i cat ion o f effort. 2. A network and boundary cond i t ion c lass i f icat ion system has been developed. Th i s c lassi f icat ion system, for the first t ime, a l lows more than one boundary cond i t ion to connect w i th any node i n the network. G r oup s o f connected boundary condi t ions, o r l i nked sets, must be so lved by special procedures invo lv ing 132 133 s imultaneous equations. The requ i red procedures for so lv ing these l i nked sets are log ica l ly and systematical ly deve loped. 3. Boundary cond i t ion equations for many c ommon hydrau l i c devices are developed i n a general and comprehens ive manner . These formulat ions permi t special devices and device appl icat ions to be extracted f r om more general formulas by selecting values o f input constants. Fu r the rmore , the solut ion to many boundary condit ions can be obta ined expl ic i t ly for the first t ime. O f part icular importance is a stable and h igh ly eff ic ient numer ica l procedure for calculat ing junc t ion losses in transient condit ions. 4. L i s t processing theory has been used throughout the analysis to insure inexpensive comput ing and op t imum use o f storage. P re -ana l ys i s steps, such as d i v id ing p ipes into an integer numbe r o f reaches and o f replac ing short pipes by p ipe replacement elements are now done automatical ly. 5. A systematic technique for numer i ca l l y testing the stabi l i ty and correctness o f the computer procedure has also been deve loped. The test procedure has conv inc ing ly ver i f i ed the correctness o f the network a lgor i thm and program. In summary, the deve loped a lgo r i thm a l lows comprehens ive transient analysis o f networks hav ing arbitrary geometry to be per fo rmed. The networks may inc lude any general combinat ions o f hydrau l i c devices and no restr ict ion is p laced on either the number o f boundary condit ions or p ipes that connect to a g iven node or the number o f sections found i n any p ipe i n the network. Procedures have been ut i l i zed wh ich keep data requirements to a r r i in imum despite the general nature o f the p rob l em be ing solved. F ina l l y , the p rogram has been carefu l ly tested and compared wi th known solut ions for s imple p ipe l ines w i th very good agreement be ing observed. The program is now ready for appl icat ion to actual d is t r ibut ion systems where it can be further tested and used for rat ional design and system development. B I B L I O G R A P H Y Ame r i c an Wate r W o r k s Assoc iat ion Standard C 4 0 0 - 8 0 , "Asbes tos -cement distr ibut ion p ipe, 4 inch through 16 inch , for water or other l iqu ids" . Ame r i c an Wate r W o r k s Assoc iat ion Standard C401 - 77 , "The selection o f asbestos-cement d is t r ibut ion p ipe , 4 inch through 16 inch , for water or other l i qu ids" . Baase, S., Computer Algorithms: Introduction to Design and Analysis, A d d i s o n - W e s l y , D o n M i l l s , Ontar io , 1978, p. 286. Bergeron, L , Water Hammer in Hydraulics and Wave Surges in Electricity, (translated under the sponsorsh ip o f the A S M E ) , John W i l e y and Sons, N e w Yo r k , 1961, p. 293. Berztiss, A T . , Data Structures: Theory and Practice, Academic Press, N e w Yo r k , 1971, p. 442. Burden , R.L., Fa i res, J .D. , and Reyno lds , A .C . , Numerical Analysis, 2nd Ed i t i on , Pr ind le , Webe r and Schmidt , Boston, 1981, p. 598. Chaudh ry , M .H . , Applied Hydraulic Transients, V a n Nos t rand R i enho l d , N e w Yo r k , N.Y. , 1979, p. 503. Chaudh ry , M . H . , "Nume r i c a l Me thods for So lut ion o f Unsteady R o w Equat ions," Closed-Conduit Flow, E d . Chaudh ry , M .H . , and Yev jev i ch , V., Wate r Resource Publ icat ions, 1981, pp . 167-191 . Chaudh ry , M . H . and Ho l l oway , M.B. , "Stabi l i ty o f M e t h o d o f CHaracter is t ics ," A S C E Proc. Water for Resource Development, Coeu r d 'A lene, Idaho, Aug . , 1984, pp. 216-220 . Co l l i n s , A . G . , and Johnson, R.L., "F in i t e E lement M e t h o d o f Wate r D i s t r ibu t ion Networks , " Journal of the American Water Works Association, V o l . 67, n4, July, 1975, pp. 385-389 . Courant , R., and Fr iedr i chs , K .O. , Supersonic Flow and Shock-Waves, Interscience Publ ishers , N e w Yo r k , 1948. Cross , H. , "Ana l y s i s o f F l o w i n Condu i t s or Conductors ," Bulletin No. 286, Un ivers i t y o f I l l ino is Eng inee r ing Exper imenta l Stat ion, U rbana , I l l inois, 1936. Epp , R., and Fow le r , A . G . , "E f f i c i en t Code for Steady-State F l ows i n Networks , " Journal of the Hydraulics Division, Proceedings o f the Ame r i c an Society o f C i v i l Engineers, V o l . 96, N o . H Y 1 , Jan., 1970, pp. 43 -56 . Evangel i s t i , G . , "Wa te rhammer Ana lys i s by the M e t h o d o f Character ist ics," L'Energ. Elec, V o l . X L V I , nos. 10,11,12, M i l a n , 1969. Fo x , J.A., "Pressure Trans ients i n P ipe Ne tworks — A Compu te r S imu la t ion , " 1st International Conference on Pressure Surges, B H R A F l u i d Eng ineer ing , C ran f i e ld , Bed fo rd , U .K . , 1972, pp. B l - 1 to B l - 9 . 134 135 Fox , L A . , Hydraulic Analysis of Unsteady Flow in Pipe Networks, M a c M i l l a n Press, L ondon , 1977, p. 216. ' Fox , P., " The So lut ion o f Hype rbo l i c Part ia l D i f f e ren t i a l Equat ions by D i f fe rence Methods , " Mathematical Methods for Digital Computers, edited by Ra ls ton, A . and W i l f , H.S., J ohn W i l e y and Sons, N e w Y o r k , 1960, pp. 180-188. Fukuda , T., "Parametr i c App r oa ch to Water Supp ly Ne two rk Ana lys is , " Electrical Engineering in Japan, V o l . 101, N o . 3, M a y - J u n e , 1981, pp. 117-122. Fuzy , O., Ha lasz , G . , and K u l l m a n , L., "S imula t ions o f Pressure Osc i l la t ions i n P ipe l ines ," Proc. of 5 Conf. on Fluid Machinery, Budapest, 1975, pp. 309-316. Gess ler , J . , "Ana lys i s o f P ipe Networks , " Closed-Conduit Flow, E d . Chaudhry , M .H . , and Yev jev i ch , V., Wa te r Resource Publ icat ions, 1981, pp . 61 -99 . Go l dbe rg , D.E., and Wy l i e , E B . , "Character ist ics M e t h o d Us i ng T i m e - L i n e Interpolat ions," Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 109, N o . H Y 5 , M a y , 1983, pp. 670-683 . G ray , C . A . M . , "The Ana lys i s o f D iss ipat ion in Wate r Hamme r , " Proc. ASCE, V o l . 119, Paper 274, 1953, pp . 1176-1194. G raze , H.R., " A Rat iona l The rmodynam i c Equat ion for A i r Chambe r Des ign , " Proc. Third Australasian Conference on Hydraulics and Fluid Mechanics, Sydney, Aust ra l ia , 1968, pp. 5 7 - 6 1 . Gusta fson, K . E , Introduction to Partial Differential Equations and Hilberi Space Methods, John W i l e y and Sons, Toronto, 1980, p. 270. Ha l l iwe l l , A .R . , "Ve loc i t y o f a W a t e r - H a m m e r Wave in an Elast ic P ipe , " Journal of the Hydraulics Division, Proceedings o f the Ame r i c an Society o f C i v i l Engineers, V o l . 89, N o . H Y 4 , Ju ly , 1963, pp. 1-21. Hodgson , J., "P ipe l ine Ce ler i t ies , " M.Eng. Thesis (Civil Engineering), the Un ive rs i t y o f A lbe r ta , Edmonton , A lbe r ta , 1983. Jeppson, R.W., Analysis of Flow in Pipe Networks, A n n A r b o r Science, A n n A rbo r , M i c h i g an , 1976, p. 164. Kap l an , M . , Belonogoff , G . , and Wentwor th , R.C., " E conom i c Me thods for M o d e l l i n g Hyd r au l i c Trans ient S imu la t ion , " 1st International Conference on Pressure Surges, B H R A F l u i d Eng ineer ing , C ran f i e l d , Bedford , U .K . , 1972, pp. A 4 - 3 3 to A 4 - 3 8 . Koe l l e , E , "Trans ient Ana lys i s o f Pressure Condu i t Hyd r au l i c Systems," Proceedings of the International Institute on Hydraulic Transients and Cavitation, Sao Paulo, Braz i l , 1982, pp . B1 .1 -B1 .38 . K o r ob i , T., Y o ko yama , S., and M iyash i r o , H., "Propagat ion Ve loc i ty o f Pressure wave i n P ipe l ine , " Hitachi Hyoron, V o l . 37, N o . 10, O c t , 1955. Kre iss , H . - 0 . , Numerical Methods for Solving Time-Dependent Problems for Partial Differential Equations, Un ive r s i t y o f Mon t rea l Press, Mon t r ea l , 1978, p. 114. 136 Lamber t , J .D. , Computational Methods in Ordinary Differential Equations, John W i l ey and Sons, Toronto , 1973, p. 278. L ister, M . , "The Numer i c a l So lut ion o f Hype rbo l i c Part ia l D i f fe rent ia l Equat ions by the M e t h o d o f Character ist ics," Mathematical Methods for Digital Computers, edited by Ra ls ton, A . and Wi l f , H.S., J ohn W i l e y and Sons, N e w Yo r k , 1960, pp. 165-179 . Oh tme r , O. , "Non l i nea r F l o w Ana lys i s o f P ipe Networks , " International Journal for Numerical Methods in Enineering, V o l . 19, Ma r c h , 1983, pp. 373-392. Pa rmak i an J . , Water- Hammer Analysis, Dove r Publ icat ions, N e w Yo r k , 1963, p. 161. Pearce, R., and Evans, E.P., "Surge protect ion o f Wate r Supp ly P ipe l ines for Juba i l Industr ia l C i t y , " Proc. Instn. Civ. Engrs., Part 1, V o l . 76, Feb., 1984, pp. 249-268 . Pearsa l l , I.S., "The Ve loc i ty o f Wate r H a m m e r Waves," Sym. on Surges in Pipelines, Proc. Ins t o f M e c h . Eng. , V o l . 180, part 3E, 1965-1966. Raman , V. , "Deve lopments in Water Ne two rk Des ign , " Journal of the Sanitary Engineering Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 96, N o . SA5 , O c t , 1970, pp. 1249-1263. Schnyder, O., "Druckstosse in Pumpen Steig le i t iungen," Schweiz. Bauzeitung, 1929. Schnyder, O., " Ube r Druckstosse i n Rohr le i tungen, " Wasserkraft und Wasserwirtschaft, 1932. Sevuk, S., "Uns teady F l ows i n Wate r D i s t r i bu t ion Systems," A S C E Proc. Computer and Physical Modeling in Hydraulic Engineering, Ch icago, Aug . , 1980, pp. 228-239. Shami r , V., and Howa rd , C .D .D . , "Wate r D i s t r i bu t ion System Ana lys i s , " Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 94, N o . H Y 1 , Jan., 1968, pp. 219-234. Sharp, B.B., Water Hammer: Problems and Solutions, Edwa rd A r no l d , L ondon , 1981, p. 144. Sharp, B.B., and Cou l son , G.T. , "Surge P rob lems i n the Wate r Supp ly Ma i n s to Ha ze lwood Power Stat ion," J. of I.E.A., Aust ra l ia , 1968, p. 83. Streeter, V .L . , "Wate r H a m m e r Ana lys i s o f D i s t r i bu t ion Systems," Journal of the Hydraulics Division, Proceedings o f the Ame r i c an Society o f C i v i l Engineers, Vo l . 93, N o . H Y 5 , Sep t , 1967, pp. 185-201 . Streeter, V .L . , "Unsteady F l o w Ca lcu la t ion by Numer i ca l Me thods , " / Basic Eng., Trans. ASME, V o l . 94, June, 1972, pp . 457-466 . S U R N A L User ' s M a n u a l , Compu te r Sciences Canada , July, 1977, p. 122. T r i k ha , A.K.., "Var iab le T ime Steps for S imu la t ing Trans ient L i q u i d F l o w by M e t h o d of Character ist ics," / Fluids Eng., Trans. ASME, V o l . 99, M a r c h , 1977, pp. 259 -261 . 137 Va rdy , A . E , " O n the Use o f the M e t h o d o f Character ist ics for the So lut ion o f Unsteady F l ows i n Networks , " 2nd International Conference on Pressure Surges, B H R A F l u i d Eng ineer ing, C ran f i e l d , Bed fo rd , U .K . , 1976, pp. H 2 - 1 5 to H 2 - 3 0 . Va rdy , A .E . , D i scuss ion o f "Hyd rau l i c Transients Fo l l ow i ng Va l ve C losure , " (by C.S. Watt , J . M . H o b b s and A . P . Bo ldy) , Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 107, N o . H Y 7 , July, 1981, pp. 956 -957 . Wa rga , J . , "De te rm ina t i on o f Steady-State F l ows and Cur rents i n a Ne twork , " Proceedings of the Instrument Society of America, V o l . 9, pL 5, Paper 54 -43 , 1954. Wat t , C.S., Hobbs , J .M . , and Bo ldy , A.P. , "Hyd rau l i c Transients Fo l l ow ing Va l ve C losure , " Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 106, N o . H Y 1 0 , Oct., 1980, pp . 1627-1640. Watters, G .Z . , Analysis and Control of Unsteady How in Pipelines, 2nd ed., But terworth Publ ishers, Stoneham, M A , 1984, p. 349. Wigger t , D.C., and Sundquist , M.J . , " O n the F i x e d G r i d Character ist ics for P ipe l ine Transients," Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 103, N o . H Y 1 2 , D e c , 1977, pp . 1403-1416. W o o d , D.J., and Char les , C O . A . , "Hyd rau l i c Ne twork Ana lys i s U s i ng L inear Theory , " Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 98, N o . H Y 7 , Ju ly , 1972, pp . 1157-1170. W o o d , D.J., and Rayes, A . G . , "Re l i ab i l i t y o f A l go r i t hms for P ipe Network Ana lys i s , " Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, V o l . 107, N o . H Y 1 0 , OcL , 1981, pp. 1145-1161. W y l i e , E.B. , "Uns teady F l ow : Introduct ion, Gove r n i ng Equat ions ," Closed-Conduit Flow, E d . Chaudh ry , M .H . , and Yev jev i ch , V., Wa te r Resource Publ icat ions, 1981, pp . 145-166. W y l i e , E.B., " The M i c rocompu te r and P ipe l ine Transients," Journal of the Hydraulics Division, Proceedings o f the Ame r i c an Society o f C i v i l Engineers, V o l . 109, N o . H Y 1 2 , D e c , 1983, pp. 1723-1739. Wy l i e , E B . , "Fundamenta l Equat ions o f Wate rhammer , " Journal of the Hydraulics Division, Proceedings o f the Amer i c an Society o f C i v i l Engineers, Vo l . 110, N o . H Y 4 , A p r i l , 1984, pp. 539-542 . W y l i e , E.B., and Streeter, V .L . , Fluid Transients, F E B Press, A n n Arbo r , M i ch i gan , 1983, p. 384. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062640/manifest

Comment

Related Items