.‘ - #5,...» ‘ ‘ )1" . ’ 1 ‘x - .A, ~- w, :.:.>.1‘r ‘3';$ p w» ~¢ .- ” V. » ‘Egfigf‘ ‘ uzwxm a " v ‘ V ; «1 I ' -¢ I ‘ wrnr‘. . 971133 B n .‘i‘wm-I _" ,,., “33$ MICHIGAN STATE UN I I! ll! lllllllll/Ull’I’Illllllll 3 129020 4919 ll This is to certify that the thesis entitled A Condensed Finite Element Analysis of Microirrigation Hydraulics Which Incorporates Pipe Components presented by Philip Gerrish has been accepted towards fulfillment of the requirementsA for Agricultural M.S. Engineering degree in Major professor Date December 222 1993 0-7639 MS U is an Affirmative Action/Equal Opportunity Institution LIBRARY Michigan State Unlverslty PLACE ll RETURN BOXtomnovoMchockoufmn want-cord. To AVOID FINES Mum on or baton date duo. DATE DUE DATE DUE DATE DUE MSU Is An Nflnmtlvo Action/Equal Opportunity Inltltwon Wanna-m A CONDENSED FINITE ELEMENT ANALYSIS OF MICROIRRIGATION HYDRAULICS WHICH INCORPORATES PIPE COMPONENTS By Philip Gerrish A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Agricultural Engineering 1993 ABSTRACT A CONDENSED FINITE ELEMENT ANALYSIS OF MICROIRRIGATION HYDRAULICS INCLUDING PIPE COMPONENTS By Philip Gerrish The hydraulic design of microirrigation systems is a tedious and time- consuming task. It was the goal of this research to simplify the design of microirrigation systems in order to conserve energy and water. Numerical solutions to microirrigation hydraulics are problematic because of the prohibitive number of network nodes which need to be analyzed. The erdsting finite element solutions for pipe networks represent pipe components as separate elements, thereby increasing an already formidable number of nodes. In this research, a partial differential equation was developed which incorporates the effect of pipe components at nodes rather than in separate elements. The equation further condenses the network matrices by making use of the virtual node concept in which laterals with evenly spaced emitters are considered single elements with derivative boundary conditions. Results are strongly correlated with existing finite element solutions which show strong correlation with empirical data. Due to the reduced number of nodes, the solution converges rapidly in few iterations. A large network was solved using the condensed finite element analysis developed in this research; where previous methods would require over 12,000 nodes to solve this network, the method developed here required only 80 nodes. Results correlate strongly to those obtained using the Backstep method. A development is proposed for a two-dimensional equation describing irrigated sub-plots with evenly spaced laterals as single, two-dimensional elements with derivative boundary conditions. Some preliminary results were obtained and found to be promising. Approved Major Professor Approved Department Chairman Copyright by PHILIP JOHN GERRISH 1993 ACIWOWLEDGEMENTS Special thanks to the members of my committee: thanks to Dr. Larry Segerlind and Dr. Raymond Kunze for their comments and involvement with the preparation of this thesis, thanks to Dr. Vincent Bralts and Dr. Walid Shayya who provided the support, encouragement, and intellectual guidance vital to the deve10pment of this thesis. And thanks to Dr. Harold W. Belcher for the initial support and guidance in my studies. Also deserving of thanks and acknowledgement are Jim Schaper, Misael Miranda, Neba Ambe, Alexander White, and Theophilus Okosun for providing an intellectually stimulating environment. My parents deserve special thanks. And 03 course, I thank my lovely wife, Lucy, for her patience and continual insistance that life should not be too serious. Table of Contents List of Tables ................................................................ ix List of Figures ............................................................... x I. Introduction .............................................................. 1 Scope and Objectives .................................................. 3 II. Review of Theory and Literature .............................................. 5 Hydraulics of Irrigation ................................................. 5 Equation Governing Emitter Flow .................................... 5 Equations Governing Pipe Flow ..................................... 7 Equations for Approximating Lateral Flow ............................ 10 Equation Governing Component Head Loss ........................... 12 Analysis of Hydraulic Networks .......................................... 13 Conservation of Mass ........................................... 13 Conservation of Energy .......................................... 14 Choice of Unkown ............................................. 14 Notation ..................................................... 15 Some Methods of Network Analysis ....................................... 16 The Backstep Method ........................................... 16 vi The Newton-Raphson Method ..................................... 20 The Linear Theory Method ....................................... 22 Finite Element Formulation of Linear Theory Method .................... 23 The Finite Element Method ............................................. 28 Ill. Methodology ........................................................... 32 Research Approach .................................................. 32 Eflect of Components on Network Solution ................................. 34 Algebraic Development ................................................ 38 Linearization of Flow Equations including Minor Losses .................. 38 Calculation of Component Head Loss. ............................... 43 Algebraic Incorporation of Component Head Loss ...................... 46 Development of Partial Differential Equation ................................. 51 Incorporation of pipe components .................................. 56 IV. Results and Discussion ................................................... 63 Evaluation of Solution without Virtual Nodes ................................. 63 Evaluation of Solutions with Virtual Nodes .................................. 69 Error introduced by virtual node concept ............................. 69 Reducing a large network to a small, manageable size ......................... 74 Stability ........................................................... 77 Methods 1 and 2 .............................................. 77 Collective emitter flow in virtual node concept ......................... 77 vii The Newton Raphson Method ..................................... 78 Continuity Requirements ............................................... 78 Discussion ......................................................... 79 Ideas for further investigation ........................................... 79 Some ideas for a two dimensional development ........................ 79 Adding the third dimension ....................................... 89 v. Conclusions and Recommendations .......................................... 93 Appendix A: Computer Code ................................................... 95 Appendix B: An explanation of the structure of input and output data files .................. 132 Appendix C: Calculation of average head along a lateral .............................. 141 Appendix D: A comparison of stability ............................................ 143 Appendix E: An alternative development for the two-dimensional flow problem ............... 1 47 List of References ........................................................... 152 viii Table List of Tables Degription $92 List of coefficients for four different approaches. 61 Element contributions to global system of equations. 62 Shape functions for nine-node Lagrangian element. 82 ix Figure 10. 11. 12. 13. 14. 15. 16. 17. 18. List of Figures Description A microirrigation system. Approximation of pressure head along a lateral. A generic element. A lateral with seven emitters. Flow chart for Backstep method. Submain unit with finite element labeling. Two ways of analyzing the same network. A comparison of solution with and without including the effect of pipe components. Graph of solution which includes the effect of pipe components against the solution which does not. A schematic diagram of a pipe component. Energy grade line of "downstream element". A explanation of references to component diameters. Selection of component coefficient. Pipe element with elbow and several emitters for development of differential continuity of mass equation. Energy grade line for element with component. Network labeled for analysis by ANALYZER. Network labeled for analysis by ALGNET. ALGNET vs. ANALYZER regression plotted. Page 11 15 17 19 24 35 36 37 39 41 44 49 52 54 64 65 66 Figure 19. 20. 21 . 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 1 b. 2b. 3b. Description ALGNET vs. KYPIPE regression plotted. ALGNET compared to ANALYZER and KYPIPE. Network labeled for analysis by ALGNET. Network labeled for analysis by DIFNET. DIFNET results compared to ALGNET. Effect of partitioning the laterals on total error. A large submain unit with 12,000 emitters. Backstep and DIFNET solutions plotted together for a single lateral. Sub-plot of irrigated field with nodes numbered for two- dimensional analysis using nine-node Lagrangian element. Describing differential conservation of mass as continuous function in two dimensions. Flow path for water to get from point x to point x+dx in irrigated sub-plot. A hydraulic topology produced by LAGRANGE. Example of 3-D finite element grid in situ. Figures In Appendix Network labeled for analysis by ALGNET. Structure of input data files for ALGNET. A network labeled for analysis by ALGNET demonstrating ALGNET’s ability to handle pipe components with more than three fittings. xi 71 72 73 75 76 80 83 87 90 92 132 133 134 Ejgpgg 4b. 1d. 2d. 3d. Description Network labeled for analysis by DIFNET. Convergence of ALGNET1. Demonstration of the effect of averaging results from the previous two iterations to calculated linearizing constants for the current iteration. Convergence of ANALYZER. Convergence of Newton-Raphson method. xii Page 138 143 144 145 146 I. Introduction The Preclassic Period, starting around 2500 B.C., marks the beginning of the development of farming communities on this continent. This date corresponds to the rudimentary origins of irrigation practice worldwide. As population densities grew, the development of more sophisticated techniques of water management became necessary. Today, with high-tech irrigation and sophisticated methods of analysis, we struggle more than ever to keep up with the growing number of mouths to feed. While doomists point out the trends in population increase and the seeming impossibility of feeding all these mouths, scientists and engineers collaborate in an attempt to outwit nature, pointing out that man’s wit is in fact a part of nature. Irrigation technology is one of the most dramatic examples of man’s intelligence affecting the course of nature. In affecting nature’s course, man must be aware of its limits. Because the earth’s total area of suitable cropland and freshwater resources are limited, there has been a push for more efficient use of existing cropland and water resources--in other words, a preference in the development of intensive as opposed to extensive agriculture. It has been estimated that by the year 2000, the area of cropland in use will be double the area in use in 1985 (Power, 1986; Holy, 1981), pushing to its limits the world’s supply of suitable cropland. Development of intensive agriculture, however, is not without some costs to the environment. The amount of top-soil decreases more rapidly under intensive Page 1 Page 2 cultivation. Increases in chemical application are often paralleled by an increase in groundwater and runoff contamination. The soil, physically and chemically, is fatigued more rapidly. Any way of reducing these costs to the environment must be investigated thoroughly if we are to look to our future, near and far. A fine case for the use of environment-friendly intensive agriculture can be made by pointing to the Mayan Civilization-~one of the greatest pre-columbian empires of this continent. The case of the Maya is a classic example of the rise and fall of an irrigation society. Their rise is due in great part to the extensive development of sophisticated irrigation systems and practices (Turner and Harrison, 1983). Their fall, although still a great mystery, is thought by many to be related to environmental decay as a result of their highly intensive agriculture (Wiseman, 1989). This history and others like it must affect our thinking today, lest we lose "the ability to understand recorded historical materials." (Okosun, 1993). The case for more efficient use of water is made simply by noting that, on a worldwide average, our present efficiency is around 37%. Add to that the fact that 80% of all freshwater resources in use today are being used for irrigation, and the case for water efficiency becomes urgent (Power, 1986). One requirement for the improvement of water efficiency is the improved control of water application. This is one of the many benefits of microirrigation. Other benefits include zero runoff, reduced labor costs, ease of chemical and fertilizer injection, and higher salinity tolerance due to stable moisture conditions. Microirrigation is highly efficient and environment-friendly; it is therefore a good candidate for future-oriented agriculture. The design of microirrigation systems is tedious and capital costs are high; Page 3 these are the two most prohibiting factors in most cases. While capital costs cannot be changed, it is the aim of this research to facilitate design by improving computer models of the hydraulics involved. A. Scam and Objectives Hydraulic networks have traditionally been solved using the backstep, Hardy- Cross, Newton-Raphson, and Linear Theory methods. Bralts and Segerlind (1985) first put linearized flow equations into Finite Element formulation. Wood (1981) and Finkel (1982) point out that pressure losses across pipe components such as elbows, tees and valves may significantly afi‘ect pressure heads in a network. These minor losses were put into Finite Element formulation by Haghighi et al. (1988). A drip irrigation system is well designed if the water is applied uniformly throughout. As a measure of uniformity, the Statistical Uniformity Coefficient was introduced for microirrigation by Bralts, et al.(1987); it is defined as one minus the coefficient of variation of emitter outputs. The objective of computer models developed for micro irrigation design, therefore, is to accurately predict the output of each emitter and give the uniformity coefficient as an indication of the design’s quality. A shortcoming of existing models is the awkwardness with which pipe components are handled; inclusion of pipe components in the network analysis makes solution cumbersome and unstable because new nodes are added to the system. The overall goal of this research is to conserve water, chemicals and energy used for plant growth through improved hydraulic design of micro irrigation systems. More specifically, the focus of this research will be to develop an improved finite Page 4 element formulation of pipe network systems in order to facilitate design. The specific objectives are as follows: 1. Assess the effect of pipe components such as tees, elbows, crosses, expansions, and contractions on the solution of a hydrach network 2. Develop a condensed finite element formulation for the incorporation of pipe components which would not increase the number of nodes. 3. Include the virtual node concept in the finite element analysis, thereby further condensing the network to be analyzed. 4. Apply the condensed finite element formulation to the design of microirrigation systems, and compare the results with those of other methods. II. Review of Theory and Literature An irrigation system is characteristically a "tree" hydraulic network system (no closed loops) with a main as the "trunk", submains as primary "branches" and laterals as secondary "branches" (see Figure 1). Along the length of each lateral are the emitters which are water outlets. The emitters are where the water is applied directly to the plant in a drip irrigation system. In a sprinkle irrigation system, the water outlets are sprinklers. A. Hydraulics of Irfiggtion 1. Equation Governing Emitter Flow For the purposes of this study, all water outlets will be called emitters. The equation describing flow in an emitter is: qe = kh" (1) where q, = emitter discharge k = emitter discharge coefficient h = pressure head x = emitter discharge exponent Equation (1). in general, describes orifice flow to the atmosphere (Wu, et al., 1979). Page 5 Page 6 I A I u—u—oo—u ,--_\/ !' ‘r '7 ELECTRKEQL‘ , , CONTRO ’/ I r---"‘ . l o ’- " ————— ' L \ I ’ I l \ \ J b / / \\ \ \\ I o’ I \ / / renrruzrn \\ ‘ I I Imessunt / c‘m‘ neouuuon \ \\ \\ I I - ’ GAUGE . \ \ ,. mun own ‘0 - . ‘l ‘ hon—o- 00—00 m. m .—-————--—.—..‘—_.Q\‘\“————_—— 2“" w , are: 1' 1 .~. WATER now MAW OAT! I 0-. FILTER METER VALVE l s / I ,1 / | ‘NELL I g I ‘ / / 5' , 1‘ ,t / Awnrsn TABLE ‘ 0 l 3 , o EMITTERS 1' m: , \ . 9 , “Id :. l g I 8> LATERAL\ / s3 / / if 02 U 83' A microirrigation system. Figure 1 Page 7 The coefficient and exponent characterize the emitter and are different for different emitters. For a pressure-compensating emitter, for example, the value of x is close to or equal to zero, as flow does not depend on pressure head. The value of x, in other cases, is indicative of the flow regime in the emitter. A value of 0.5 indicates fully turbulent flow, whereas a value of 1 indicates laminar flow. 2. Equations Governing Pipe Flow A generally accepted form of the differential equations governing fluid flow is the Navier-Stokes equation: Pg—Z=-Vp+pe+uV2V (2) where; v = Laplacian operator, p = pressure, p = density, g = acceleration due to gravity, u = viscosity, V = velocity, Page 8 The particular solution of the Navier-Stokes equation for pipe flow is the Darcy- Weisbach equation: = E .23 3 AH, f D 29 l l where AH, = pressure head loss due to pipe friction f = friction factor L = length of pipe D = diameter of pipe V = velocity of fluid g = acceleration due to gravity The fi'iction factor, f, in the above equation is a dimensionless wall shear, f= t° 1 . _ V3 89 It is found to be f =%: (where Re = Reynold’s number) for laminar flow. This value, however, is not so nicely defined in the transition and turbulent flow regimes. To determine f in these flow regimes, it is common practice to refer to a chart called the Moody diagram. Some empirical equations have been derived and are successful in approximating the friction factor, f. Another equation which may be used to calculate head loss due to pipe friction is the Blasius equation: = 8(4lba 4:224: ht. W—gv DS+bL (4) Page 9 where; v = kinematic viscosity q = flow rate L = length of pipe D = pipe diameter a = constant b = constant g = acceleration due to gravity For PVC pipe, the values of a and b were determined by von Bernuth and Wilson (1989) to be: a = 0.316 and b = -0.25. Hence the equation, hr. = KLvo.2sq1.7SD-4.7s (5) An empirical equation often prefered for its simplicity is the Hazen-Williams equation: - ksysL 1.952 AHD- Cfiasthmq (6) where k,,, = 4.73 for English units (D and L in feet) km = 10.7 for International System of units q = flow Cm, = Hazen-Williams roughness coefficient This equation approximates head-loss over a limited range of Reynolds numbers, a drawback which should be considered especially when working in laminar or Page 10 transition flow regimes. While any of the above pipe flow equations may be used, the algebraic development which follows uses the Hazen-Williams equation. This equation will be used for the sake of comparison with other methods which use the same equation. 3. Equations for Approximating Lateral Flow Flow in a microirrigation lateral line (a pipe with emitters) can be approximated by constructing a dimensionless energy gradient curve (Wu and Gitlin, 1975). Basic assumptions are (a) flow from all emitters along the lateral is the same, and (b) emitters along the lateral are evenly spaced. Dimensionless head drop (AH/AH) is plotted against dimensionless length (x/L) and the curve derived is an exponential decay function: R1=i—P:;-=1-(1-i)m*1 (7) where, AH, = head drop at point i AH = total head drop in lateral i = x/L, where x = distance from origin, and L = total length of lateral. m = 1.852 from Hazen-Williams equation Page 11 Pressure Head Along a Lateral 120 100 ._ so 0 60 fi \\ l- 40 - 20 O s : : t : : : : : 0 0.2 0.4 0.6 0.8 1 Length Ratio, x/L Pressure Head Figure 2. Approximation of pressure head along a lateral. Page 12 Head at any point i along the lateral can then be described as: h, = Ho — R1AH a: Ri'AH' (8) where, H0 = head at i = 0 RAH = head loss due to pipe friction at point i R1’AH’ = head loss or gain due to elevation at pointi The curve described by equation (8) is shown graphically in Figure 2. With emitter flow described by q = kh", and constant emitter flow along the lateral, flow at point i along the lateral is: a, = k(h,)" = km, - 12,1111 :1: Ri'AH')‘ (9) Or, substituting the equality, q, = kHo" for flow from the first emitter on the lateral gives an equation for flow fiom the lateral at any point i along the lateral: (11 = qoll - R1(AH/Ho) 2t R1'(AH'/Ho)]“ (10) 4. Equation Governing Component Head Loss An abrupt change in pipe geometry causes turbulence. Where turbulence occurs, energy is lost. Pipe components (tee’s, elbows, valves, etc.) present abrupt changes in pipe geometry. Their presence therefore implies loss of energy. A pump, the exception of course, would increase energy. The energy lost1 at a pipe component 1Note here that the terms energy and head are used interchangeably. This is because in hydraulics, total energy is commonly expressed in potential form. 'Refer to next section under the heading, Conservation of Energy. Page 13 is equal to some fraction of the fluid’s kinetic energy at that point: 2 AHc=kc2Lg (11) where AH, = pressure head loss due to pipe component V = the velocity of the fluid g = acceleration due to gravity ke = a constant peculiar to the component used 2 The velocity head of water, 2V3 , is the water’s kinetic energy expressed as potential energy (head). B. Analysis of Hydraulic Networks As with any physical system, both conservation of mass and conservation of energy must be satisfied. These dictate the continuity equations throughout the system. 1. Conservation of Mass The conservation of mass of a fluid implies that flow in equals flow out: 2 q1n=2 gout: (12) When this continuity is met at several points in a system, the result is a system of Page 14 simultaneous equations. 2. Conservation of Energy The conservation of energy of a fluid is expressed by Bernoulli’s equation: 2 {—g + h + z + h1, = constant (13) where V = velocity g = acceleration due to gravity h = pressure head 2 = elevation h, = head loss due to fiiction Simply stated, this equation says that the sum of the kinetic energy (in potential 2 form), 2V3 , plus the potential energy, h+z, plus the heat energy lost due to fiiction (in potential form), 11 1., is equal to the total energy and must therefore be constant throughout the system. 3. Choice of Unkown The set of simultaneous equations describing a hydraulic network can be expressed in terms of either hydraulic head (potential energy) or flow as the unknown. Choosing flow as the unkown has the advantage that many of the equations in the set of simultaneous equations will be linear (Jeppson, 1976). In fact, if no closed Page 15 loops exist (i.e. a tree network, as commonly encountered in irrigation systems), all of the equations will be linear. While this is clearly a great advantage, it is countered by the disadvantage that the boundary conditions are expressed in terms of potential (head). A pump, for example, will deliver a certain hydraulic head, the outflow of which depends on the solution of the entire network. Also, the derivative boundary condition described later as the Virtual Node concept is a potential (head) gradient. The great disadvantage of choosing potential as the unkown is that the resulting equations are non-linear. The attractiveness of this choice, however, is due to (a) the systemmatic facility with which the set of equations is formed, and (b) the ability to apply boundary conditions essential for solution. Note that the unkown chosen in this research is potential energy which, as stated earlier, is the sum of elevation and pressure head, 2 + h. This sum is also refered to in this thesis as hydraulic head, and is denoted simply by h for hydraulic head within an element and H for hydraulic head at a node. 4. Notation The network analysis technique used in this thesis is called the Finite Element Method. The use of Finite Element notation throughout this section will facilitate presentation of the different network analysis i j techniques and will assure consistency of . , Figure 3. A generic element notation throughout this thesis. The reader, therefore, will be acquainted with Finite Element notation at this point. Page 16 Figure 3 shows a generic element--a one-dimensional element as defined by Segerlind (1984). The element, (e), lies between two nodes, i and j. In the case of hydraulic network systems, an element refers to a length of pipe. A node must be assigned to every point at which (a) head is known, (b) head is to be calculated, or (c) there is an abrupt change ("discontinuity") in head. For example, to calculate emitter flow, it is necessary to know hydraulic head at that point; so a node is assigned to that point. A variable with a superscript in parenthesis pertains to an element, whereas a variable with a subscript pertains to a node. The variable, q"), for example refers to flow through pipe element 7, while q, refers to flow through emitter node 5. Equations developed here are non-linear and their solution is numerical. Throughout this thesis, the subscript, n, denotes the number of the current iteration. Likewise, n-l refers to the previous iteration. C. Some Methods of Network Analysis 1. The Backstep Method For a string of emitters (a lateral), where one end represents a source with known pressure, the Backstep Method can be used to calculate the pressure heads at each of the emitter nodes along the lateral. In Figure 4, for example, node 1 represents a known head while nodes 2 through 8 represent emitters at which head is to be calculated. To begin the solution, an initial guess is made for head at node 8. Next, the Page 17 Reservoir (known head) l"" Emitter' J--E|enent 4 rigure 4. A lateral with seven emitters. Page 18 head at node 7 is calculated by satisfying the conservation of mass at node 8. That is, flow out of node 8 (emitter flow) is set equal to flow in (pipe flow). Pipe flow is calculated by rearranging the Hazen-Williams equation, ql7>=( i:§a)T;—53 (14) where, km: ksysL (15) 1.852 CH" 04.87 L = length of pipe element 7 CW = Hazen-Williams roughness coefi'icient for element 7 D = diameter of pipe element 7 For example, the equation expressing conservation of mass at node 8 is: q(7) _ qa : o (16) Substituting equations ? and (14) into equation (16) gives: _ 1 I—i— where, H - z = hydraulic head minus elevation = pressure head. Page 19 Initial guess forhatnodeB Increnent h at node B - l Satisfy Cons. of Mass at each previous node consecutively 00 atmm1:kmmhmfl Figure 5. Flow chart for Backstep Method. Page 20 With the initial guess for H8, H7 is calculated. From the expression for conservation of mass at node 7, H, is calculated, and so on until Hl is calculated. This value for H1 is then compared with the known value for head at that point. If these values are within a specified error interval, then the procedure is finished; otherwise, the initial guess for h, is incremented, and calculations are repeated. The solution is thereby arrived at iteratively as shown in Figure 5. The advantage of the backstep method is its straight-forward nature and hence its facility of formulation. The great disadvantage it has is the long convergence time required due to the large number of iterations required. Of the three methods presented in this section, this is by far the slowest. 2. The Newton-Raphson Method This method is based on a truncation of the Taylor series which takes the form, R (x -1) x12 arm1 R’( "-1) ( l where x = the unknown R(x) = the residual equation 11 = the number of current iteration This is an iterative solution to the equation R(x) = 0. An equation of this form is known as a residual equation. In the example of a string of emitters (Figure 4), a residual equation and its first derivative must be written for each node. The residual equation for any node in a pipe network system is simply, mm = qin - q“, = 0, where h Page 21 (head) is now the unknown. Written for all nodes as a system of equations, the matrix formulation takes the following form: {ha} ={h,,-1} - [D] '1{R(h,,-1) } (19) where [D] = the J acobian matrix of derivative elements {R(h)} = the residual vector {h} = the vector of unknowns (head) '8121 are, 8121‘ a—h'; 6—112 0 e e "a—II'; BE 2‘12 332 [D] = ah, ah, an, (20) 812,, 812,, 812,, _'d’IT, an, 517:2. Determining the J acobian matrix makes this method rather cumbersome when seeking a generic formulation for networks. Also, a disadvantage of this method is its instability, especially when the unknown is head. (Greater stability is achieved when flow is the unknown.) Because of its instability, this method depends greatly on the accuracy of the initial guess. Otherwise, this method has the great advantage of quadratic convergence, which means rapid solution. Page 22 3. The Linear Theory Method This method sets up a system of linear equations by "linearizing" the flow equations (Jeppson, 1976). This is achieved by observing that, -o.« (H1 -HJ ) "'1 and k (s) = kSYUL . (k(0))0.56 ’ Cfiésssz,” q‘.) =C‘C’ (Hi-Hj)n’ where C(‘lg The term, 0‘", is treated as a constant and its value can be determined because (Hi-H1) rm is known. (Note that n-l refers to the previous iteration.) Linearization of the emitter flow equation in terms of hydraulic head is achieved by noting that, q, = k,(H-z)" = (k,‘—H"—?II>H (21) sothat q, = Cfin , where C,, = krfl'fiz—E n-l The conservation of mass at node 6, for example, is now expressed in the following form: c151 (HS-H6) -c‘°’ (Hg-H7) -c,H,=o (22) Again, notation is important. 0‘" refers to the linearizing constant for pipe element flow. C, refers to the linearizing constant for emitter node flow. This method in general is quite stable and is not highly dependent on the accuracy of the initial guess. The solution converges in relatively few iterations. Page 23 4. Finite Element Formulation of Linear Theory Method One-dimensional Finite Element notation (Segerlind, 1984) is used to number nodes and elements of a hydraulic network system (Bralts et al., 1987). Figure 4 and Figure 6 show examples of this numbering system. A network numbered in this fashion is readily put into matrix form by adding the contribution of each element to the global system of equations. The global system of equations takes the form, [K] {H} ={F} (23) which can be written in residual form as, {R}=[K] {HI-{F}={0} (24) where {R} = global residual vector. [K] = global stiffness matrix. {H} = global vector of unknown heads. {F} = global force vector. The element stifl'ness matrix represents an element’s contribution to the global stiffness matrix. The element stiffness matrix, in this case, consists of the linearizing constant pertaining to an element: (25) (e) .. (e) [klel]=[c C ] -Cle) C(e) Page 24 I n) (a) 3 0 Figure 6 Submain unit with Finite Element labeling of nodes and elements Page 25 where; C") = linearizing constant, k“) = element stiffness matrix. The element force vector is an element’s contribution to the global force vector. When minor losses are not accounted for, it equals zero. The procedure by which element vectors and matrices are added to the global matrices is called the Direct Stifi‘hess Method (Segerlind, 1984). The generic form of an element stiffness matrix is, 2.1 k2.2 m = ' ' (26) [k l [k 1 This matrix is then added to the global stiffness matrix as follows: k1.1 is added to K1,: k1.2 is added to Kid ls,’1 is added to 1%, k.” is added to K”. where, k = value in element matrix K = value in global matrix Emitters may be incorporated into the global stiffness matrix by considering them to be separate elements (Bralts et al., 1987). Page 26 q=Cg(H1-Hacm) (27) where Hm is equal to zero (atmosphere) and C, is the linearizing constant for the emitter flow equation. The result is that, C. is added to K“; emitter contributions and up on the diagonal of the global stifl‘ness matrix. The global stiffness matrix which results is a banded symmetric matrix. The iterative solution of the resulting system of equations was written in computer code by Bralts and Segerlind (1985). An advantage of the Finite Element formulation for the solution of hydraulic network systems is its systematic simplicity. Because of this, it is well suited to the development of a universal network solver. Bralts and Segerlind (1985) list several other advantages. Accomodating Pipe Components. The cummulative efi‘ect of several pipe components in a hydraulic network is often substantial. While losses due to pipe wall friction often dominate, the "minor losses" due to turbulence in pipe components is seldom negligible. As pointed out by Villemonte (1977), the term "minor losses" is often a misnomer as their effect is frequently "major". Haghighi et al. (1988 and 1989) proposed that pipe components may be added to the system of equations as separate elements. This is done by adding an element matrix of linearized component coefficients to the global stiffness matrix for each pipe component. The component head loss in a component element is calculated by: Page 27 Va 8 3 AH: 3 H1 ' H: ‘ I‘d—2'3 ' kch-ga This can be linearized and rearranged to form, q 3 Cc(H1 " H1) fladl where C, = 8kg 1 c n- The element matrix for a pipe component is then added to the global stiffiiess matrix by the Direct Stifiness Method. The component element matrix is similar in form to those for pipe elements: Cc 'ce (cl . [kn ] 'Cc Cc The above development works only for certain two-node components such as elbows and valves. A tee, however, is a component with three "fittings"; it is represented by an element with three nodes (Haghighi et al., 1989). This approach has the advantage that the equations are assured global continuity of zeroth order, C°, meaning that discontinuities do not exist at nodes (see discussion of continuity in Results and Discussion section). The additional nodes added by pipe components, however, mean additional equations, which means increased computer time and memory requirements. Page 28 D. The Finite Element Method One way of approximating the solution to an equation is to multiply its residual form by a weighting function and set the integral of the product equal to zero. This procedure is especially useful in the solution of differential equations. While there are several weighting functions to choose from, the weighting function employed in Galerkin’s method is perhaps the most widely used. Consider the equation, 0%‘94’4'0'0 (28) Employing the product rule for derivatives together with Greene’s theorem, the second-derivative term can be broken down into two first-derivative terms, one of which represents the intra-element residue which should go to zero (refer to Segerlind, 1984 or Dhatt and Touzot, 1984). The function is multiplied by its weighting function (or shape function) and integrated over an element at each node of the element. A linear element has two nodes, i and j; the integration at node 1 yields: 1’: 31634; 11 3: fD-a—x- axd" -fc1v,¢dx +IQN1dx =- 0 X1 11 ’1 where Ni = shape function at node i. Page 29 Integration at nodej yields: x, SN ’9 "1 £03353}!th - £aNjodx + imam: -- o where Nj = shape function at nodej An element’s contribution, therefore, to the global system of equations is written in matrix form: x, a [N] r X: x, 29 - r r = Inwadx learn ¢dx+foINI dx 0 r ) 1 1 I where [N] = [N, Ni] = shape function matrix for element (e) Shape functions are derived so that, ¢lel = [N] {Qlel} Linear shape functions, therefore, satisfy the following: ¢(°) = N101 + N101 where potential, ¢‘°’, is a straight line connecting nodes i and j. Hence they take the form: 3 ll '3 l >< “2 l Page 30 where Xi=xatnodei Xj=xatnodej L = length of element While this research employs linear shape functions only, the same concepts developed here can be used in conjunction with shape functions of higher-degree polynomials to achieve more accurate results. Continuity Requirements. From equation (29), note that a first derivative term is integrated. The first derivative of the function must therefore be defined, requiring continuity of first order, 0‘, throughout an element. Also, interelement continuity of zeroth order, C°, is required, meaning that the value for nodej of an element must equal the value for node i of the next element. In the case of linear elements, the first derivative is undefined at the nodes. The Virtual Node Technique. A way to reduce the size of the global stiffness matrix is by merging the effect of a string of several emitters into a single element. This is achieved by considering the efi‘ect of these emitters to be continuous throughout a single element, thus having the effect of a derivative boundary condition along the element. The nodes of such an element are called virtual nodes or virtual emitters (Kelly, 1989; Bralts et al., 1993). This technique considerably reduces the number of nodes in a system and hence the size of the global stifl'ness matrix. Page 31 The implementation of this technique requires that the flow equations be rearranged to describe head as a residual equation in differential form. This equation is then solved simultaneously for all nodes in the system using the Finite Element Method. III. Methodology When applying the Finite Element formulation of the Linear Theory Method to the analysis of a medium-size microirrigation system, a problem which arises is the prohibiting size of the global stiffness matrix. A drip~irrigated plot of one hectare, for example, with emitters spaced one meter apart and laterals spaced 1.5 meters apart would contain 6,600 emitters and would require 6,667 nodes for head calculations. The size of the global stiffness matrix would then be 6,667 x 6,667. Depending on the arrangement, this number could be increased up to two fold due to the presence of pipe components. Such a matrix would prohibit the use of conventional personal computers. Also, when pipe components are taken into account, an elaborate "first guess" procedure is necessary to initialize iterations. In this research, a formulation will be developed to eliminate the extra nodes added by pipe components and thereby eliminate the need to closely approximate nodal values for the solution to converge. A. Research Approach The objectives of this research are restated below, each followed by the research approach employed in its development. 1. Assess the cumulative effect of pipe components such as tees, elbows, crosses, expansions and contractions on the solution of a hydraulic Page 32 Page 33 network. This objective will be addressed by using an existing pipe-network program called ANALYZER (Haghighi et al., 1988, 1992; Mohtar et al., 1991; Shayya et al., 1988) to compare the solution of a hydraulic network not including pipe components to the solution which includes the effect of pipe components. 2. Develop a condensed finite element formulation for the incorporation of pipe components without increasing the number of elements. This will be accomplished by developing a formulation in which a pipe component is represented at a node rather than as a separate element. Because a node at a junction is needed anyway, no new node is added as a result of a component at that junction. The efl'ect of a pipe component will be accounted for in pipe elements downstream of that component. 3. Include the virtual node concept in the finite element analysis, thereby further condensing the network to be analyzed. A string of evenly spaced emitters (a lateral) will be considered a single element with a derivative boundary condition describing out-flow as a continuous function along the length of the element. This way, where previous methods required a node at each emitter along a lateral, this method requires only two nodes (one element) to represent the entire lateral or lateral segment. Page 34 4. Apply the condensed Finite Element formulation to the design of microirrigation systems and compare the results with those of other methods. A pipe network will be analyzed and the results will be compared to those of ANALYZER and KYPIPE (Wood, 1980; Wood and Charles, 1972, 1973; Wood and Rayes, 1981). B. Effect of Commnents on Network Solution The effect of pipe components on the solution of a hydraulic network is often substantial. The following scenario emphasizes the need to include the effect of pipe components if the solution is to be meaningful. Figure 7 shows a hydraulic network; for now, consider only the first diagram in the figure. The solution to this network is achieved using the ANALYZER program. In Figure 8 and Figure 9, two solutions are compared: the solution not including pipe components is compared with the solution which includes the effect of pipe components. Note that if the effect of pipe components is not included, the head calculated at node 1 is not significantly difl‘erent than zero, because the error produced by not including pipe components is greater that the head at node 1. It is apparent here that the solution which does not include the effect of pipe components has limited meaning. Now consider the second diagram in Figure 7. This shows the same network but with pipe components represented by nodes instead of elements. The total number of nodes is reduced from 40 to 27--about two thirds! Page 35 Network In lthh 0 (mt ls mpresmted by m elemt Nelmkl Ihld1 i representedn Wow 5 Figure 7 Two ways of analyzing the same network. Page 36 Ind “it H did it I II anal-b - mid- III! and and I 3a 1.3 m w: 1 :0 I! m 0 3 31 LII I at i Y E WI! 4 W IJI II III-II I w M him 37 I M 2.5 M. d Huh 3 7 4.14 2.! I u 2.3 x 1.1m I u an fl fir d we II III M II I! 2. I! u: 1.15 I: III 1.! 14 I'D II? II I: 73 II I37 III I1 I37 III II nu ma II III: 37.30 to In I! 1| III I! 22 III In 23 III I'lI II II III II (II III a In I. 17 7.04 III a 78 w 2! 7.! u: m 7.! w 31 III III] I III 191 53 14.23 IIII II III! IIII 3 III: DA 3 IIII 17.! 37 III 47.! I 109 IN Figure 8 Comparison of solution including pipe components and solution which does not include the effect of pipe components. Page 37 50 .p. O Mflflisduiiml'lhcamuds m (.4 o o I —L O 1 J 0 51015 20 25 3O 35 4O ANALYZER solution without components 0 [- output data—y = x line I Figure 9 Comparison of network solution with and without pipe components. Page 38 The above observations demonstrate the benefits of developing a simple and efficient method for including the effect of pipe components in hydraulic network analysis without increasing the number of nodes. C. Alggbraic Development The development which follows will result in an algebraic solution for pipe network analysis. Because the residual equations must be linearized, the solution to the resulting system of equations is therefore iterative. This formulation will include the minor losses2 due to components such as tees, elbows, and crosses. Component head loss3 will be calculated as a drop in hydraulic head in the elements immediately downstream of a component. An advantage of this formulation is noted in the fact that no new elements are added to the network, thereby adding no new equations to the system. 1. Linearization of Flow Equations including Minor Losses Node i in Figure 10 represents a tee at a pipe junction. Any node such as node i which represents a pipe component will from here on be called a component node. Flow from node i to node j will be called positive flow. If flow is positive, then the 2 Note here that the term minor losses refers to the cumulative effect of head loss due to components in an entire network. 3 Note here that the term component head 1088 refers to the hydraulic head loss in a particular pipe element due to the presence of a component attached to the upstream side of the element. Page 39 Tee---’e" :5 a IL... .sr»-- Upstream Element i l ‘§L\‘EE§\-- Downstream Elements Figure 10. A schematic diagram of a pipe component. Page 40 total head difi'erence between i andj is: AHT = H1 - Hj _ AHC (30) where AHc = head loss due to the pipe component, AH... = total head loss due to pipe friction H, = head at node i H). = head at nodej Head loss at a component node is presented in Figure 11 as a sharp drop, or discontinuity, in the energy grade line; here, the logic leading to equation (30) is visually apparent. The energy grade line is a graphic representation of total energy-- kinetic plus potential--and is given in terms of hydraulic head. Thus, the general equation describing the relationship between flow and head difi‘erence in an element is: H1 .. H1 - AHc = kle)q1.852 (31) or, rearranging slightly, H1 _. Hj = klelq1.852 + AHc (32) While equations (31) and (32) appear trivially different, they are quite different in terms of global formulation and stability. The two methods shall be discussed separately as Method 1 (equation (31)) and Method 2 (equation (32)). Page 41 :: I I\_k(e)ul.852 :1: 1 Energy (lend) -" Figure 11. Energy Grade Line of downstream element. Page 42 Method 1. Rearranging equation (31) gives an expression for flow: 0.54 quI=[H1’HJ'AHc) (33) k (a) This flow equation can then be linearized by observing that q" = (H1'HJ‘AHC) ”C(e) (34) where CM = (Hr‘Hj‘AHclfagiflikm)'0'“ Method 2. Equation (32) is rewritten as follows: 2 H1 _ H] = klelql.852 + k _V_ Substituting Q/A for V, H1 .. H: = klelq1.852 + kc H1 _ Hj = (kIquo.852 + T‘9’q)q where k‘°’= 4'73L 0&352D6J‘7 Page 43 (e) = 8 and T kc “204 Or, in linearized form, qn = (H1_Hj)ncle) (35) where C ‘°’ = 1 k (C) q:::52 +T (G) qn-1 2. Calculation of Component Head Loss. The component head loss term requires that velocity be known. The velocity head will therefore be based on hydraulic heads calculated in the previous iteration. Note here that at least one iteration must be performed before component head loss is taken into account. This has the advantage of automatically approximating nodal values before including the effects of components. Because continuity is expressed in terms of flow, equation (11) will be rearranged by replacing velocity with flow over area: = i. = aqz AHc kc ngz kc 91:sz (36) where A is the cross-sectional area of the component Di is the inside diameter of the component The diameter, D,, in the above equation refers more specifically to the inside Page 44 Component heod loss colculoted using Kan W Dle) “To. 3. ‘m\‘g] VL /— Component heod loss colculoted using KC180 Figure 12. An explanation of references to component diameters. Page 45 diameter of the component (see Figure 12) where the pipe element and component are joined. In the case of positive flow, this occurs at node i of the element, hence the subscript. It is important to note that this diameter is specific to the pipe element and not to the component node. The component diameter data is therefore stored as part of the element data file for the computer program discussed in a later section. Method 1. Substituting Q from equation (33) into equation (36), we get: 2 He - Cr 2 I kle) 9'" D1 3 H -H -AH 1'“ AHC = kc: 2 I[ 1 kjle) C) (37) 97‘ D: This equation must be solved numerically for AH... but a first approximation can be made by noting that the exponent is approximately equal to one (1.08 s 1). Setting the exponent equal to one and solving for AHc yields: c (38) (1((9) +T) where T, = k 91:le Page 46 Equation (38) can then be used as a first approximation to solve for AH, in equation (37). The numerical method employed here is the Newton-Raphson method, chosen for its rapid convergence. Method 2. While component head loss is not incorporated as a separate head loss term, it is accounted for by the T‘"q term. Again, the resulting equation must be solved numerically for flow: __ (Hi-H1) n-1 k (e) gist-1952 +Tle) gm1 11-1 (39) A first approximation is made by setting qua“ = q‘. q can then be solved for: q = ML. k (9) +T(O) This first approximation is then used as a seed value for the Newton-Raphson numerical solution of equation (39). 3. Algebraic Incorporation of Component Head Loss The global system of equations is made up of several element contributions. An element is affected by a component only if it is located immediately downstream of a component node. For every node which represents a component, therefore, the computer algorithm must first locate all elements touching the node and then determine which of these elements are downstream of the node. Page 47 Elements downstream of a component node. The elements found to be immediately downstream of a component node will from here on be called simply downstream elements (see Figure 10). These elements are treated differently than other elements. First, for each downstream element, the component head loss, AHc, is calculated as indicated in equation (11). Component head loss is based on the heads calculated as a result of the previous iteration and is treated as a known in the current iteration. This head loss, AH” is incorporated in the linearized coefficients (C‘°”s) for downstream elements as seen in equation (7). Method 1. Component head loss contributes to both the linearized coefficients as in equation (34) and to the right-hand side of the equations as follows. The equation describing conservation of mass at node i of a downstream element has the constant, C ‘" T1, added to the right-hand side. Added to the right side at nodej is -C "’ T1. In matrix formulation, this is expressed as a contribution to the forcing vector: 6 - c(elT1 4 {f()}-{_C(.)T1} ( 0) The above equation applies only to the condition of positive flow as previously defined. When flow is negative (i.e. when flow is from node j to node i) and nodej represents a component, the element contribution is then: Page 48 _C(e) T {f(e)}={C(°) T11} (41) Method 2. Component head loss is incorporated only in the linearized coefficient, 0‘", as defined in equation (35). Selection of Component Coefficient, 12,. A computer program called ALGNET was written based on the preceding algebraic development. The algorithm used for selecting a component’s loss coefficient is described here. For a component with only two fittings, e.g. a coupling or an elbow, only one coefficient is needed, whereas a component with more than two fittings, e.g. a tee or a cross, requires at least two coefficients for the calculation of minor losses. The nodal data file for the computer program, ALGNET, contains two coeficient values for each component node. The first value, kcm, is used if an element exists upstream of the component node which lies in a straight line (180 degrees) with the element under scrutiny. The second value, km, is the default value. If no upstream element exists at 180 degrees, then the coefficient, k,, is assigned the value of kw. In Figure 13, for example, the pipe component coefficient chosen for downstream element (2) is km, because an upstream element exists, element (1), which lies between 160° and 200° from element (2). The coefficient chosen for element (3) is k,90 because no upstream Page 49 (3) ‘-}§;5 (2) 2 o : 20 degrees Figure 13. Selection of component coefficient. Page 50 element exists between 160° and 200° from element (3). See also Figure 12. Page 51 D. Development of Partial Differential Quation Consider the pipe element depicted in Figure 14. Conservation of mass dictates that flow at x equals flow at x+dx plus the out-flow from the section, dx. As discussed in the literature review, a virtual node approach requires that the out-flow be considered a continuous function of x. The flow lost per length along a section of pipe can be formulated as a constant flow gradient: Zia. .. ”#557” (42) 6:: L where n. = number of emitters on lateral k = emitter constant (same for all emitters on lateral) h = average head along lateral (see Appendix C for calculation) 2 = average elevation along lateral L = length of lateral Conservation of mass for the pipe element in Figure 14 requires that: qx 3 qxedx + aaqx‘dx (‘3) Equations relating head loss to flow can be given the general form, Ah 2 aq‘x Page 52 Hi E>idfi-a Hj | l x x+dx Component Emitter Figure 14. Pipe element with elbow and several emitters. Page 53 where, ab = head loss due to pipe friction x = length of pipe m = exponent: 1.852 for Hazen-Williams or 2 for Darcy-Weisbach. a = coefficient: A for Hazen-Williams or f 16 for Darcy- C A106 . 87 1; 3 D5 Weisbach (see lit. review) Over a distance, dx, this equation takes the form, 612 5?: dx = aqmdx Solving for flow, . = e a) '3 and linearizing the partial derivative, q = 0%}; where, Substituting into equation (43) gives, ah ah D a. aq. xx + = D x x+dx 8:: X dx (45) Page 54 I I l I \-_- Figure 15. Energy grade line for element with component. Page 55 where, db = _8_I3 62h Dxa—x x+dx 0" B): x + Dx axzdx (46) Substituting equation (46) into equation (45) gives, all: 3a. = 4., x5; + ax 0 ( l where 232'- can be treated as a constant gradient, Q: x _ aq. _ n,k(H-§)‘ (43) o " “a? ' L or as a linearized function of h, Gh: G h = aq. = (flak (17:31,, (49) 6x I. h Results of both approaches are compared in the Results and Discussion section. Equation (47) is the governing equation for lateral flow which is was the equation presented by Bralts et. al (1993). This equation, however, does not account for the presence of pipe components. Page 56 1. Incorporation of pipe components The fundamental problem of pipe components is the energy discontinuity which they present. This becomes especially problematic when a component is treated as a part of its downstream element. As a result, the energy grade line (Figure 11) contains a discontinuity within an element as shown in Figure 15‘. Intra-element discontinuities are often prohibitive because their derivatives are undefined; however, this one can be dealt with because it occurs at a node. Again, two approaches shall be considered. Method 1. The first approach accommodates the discontinuity in the head equation, using a linear shape function to weight the residuals and employing Galerkin’s method. Essentially, the function shown in Figure 15 is a linear function with the following boundary conditions: The equation for head is then derived as a linear function of s: h = (Hi—AHC) +( Hj-H2+AHC)S ‘Note here that the coordinate system has been changed from the global system, x, to the local system, s. Because scale remains the same, derivatives are not changed. Page 57 The linear shape functions are: Integration at node i for the second-partial-derivative term (or D-term) then yields: dN f D d 51‘ %h-ds= €(H1-HJ-AHC) and at nodej: dN dh [Dial d—sds= (—-H1+HJ+AHC) Equation (50) and equation (51) combine to take on the matrix form, dINIT db D 1 '1 1 _ D 1 fD— ds dads: L[-—1 1E3} 'L'AH°{-1} or, in matrix notation, L a [N] T ah = (e) (e) _ (e) (e) {a as Eds [kn 1w } If, lAHc (50) (51) (52) (53) Page 58 Likewise, integration of the G~term yields: folNl Th ds =-%—L[: :Jgj} - Est-A343 (54) = [kéfll {Hm} - {fé°)lAHé°’ And integration of the Q-term yields: ‘ 1 [om Tds = 9f:{ } = {11"} (55) o 2 1 If emitter flow is considered a linearized function of head (in other words, if emitter flow is expressed in the G-term), the following residual equation results for element, (e): {R W} = ( mg") + [ké‘H ) {Hm} - mi" } AHé" (56) Otherwise, if emitter flow is treated as a constant (in other words, if emitter flow is expressed in the Q-term), the residual equation is: {em} = [kg°’1 {HM} - affixing" - {153"} (57) Method 2. The second approach accommodates the energy discontinuity again by adding Page 59 pipe head loss, AI'IP, to component head loss, chz AHé" = mg” + AHé" Total head loss over length, L, (Figure 14) takes the form: 92. = m 2 58 ex aqL+Tq (I Site de‘ where T = as in the algebraic development. L = length of element Equation (58) can be solved for linearized flow as follows: gfiL = (agfi’ffiL + Tgn_1)qn L .33 aqg‘IfL + an,1 3" q”: Similar logic to that depicted in equations (45) through (47) results in the following general equation: 62h D +Gh-Q=o (59) where, D: X L aq‘HL + TqL.1 Page 60 Equation (59) is the general form of a second-order partial differential equation. The coefficients, D,, G, and Q are listed in Table I for both methods. Page 61 Table 1: List of coefficients for four different approaches. D G Method 1: 1 db (71(1) n 0'1 midi-115) " n-1 n-J. L aq‘HL + Talus midi-a " n-1 L aq‘HL + TqL.1 n.k (IT-ax L E n-1 It is worth noting that the results of the algebraic development are in agreement with the equations developed in this section. Each method’s contribution to the global stiffness matrix and forcing vector is listed in table II. Page 62 Table 11: Element contributions to global system. Stiffness Forcing Matrix Vector Meth°d1= mg") {133"} + {$th 3c. _ 0 3:: Method 1: [kg-I] + [kéell {3.)}AH‘5-I + {falel}AHéei aq, - o h 6:: Method 22 [1:30)] {fold} ac. ,_ '3; 9 MethOdz [kD(.)] + [ka‘Ol] 0 3g, - G h 3:? The vectors and matrices are listed below: (.) I 2 1 '1 [kn 1 Ll-l 1] (0) .fl2 1 [kg] 6L 2] (e) = CL 1 m .2 1 if. l L{_1} (e) 3 GL 2 u... .41} IV. Results and Discussion A. Evaluation of Solution without Virtual Nodes The effect of pipe components in a hydraulic network system is incorporated in the ALGNET and DIFNET computer programs (see Appendix A for code). ALGNETI encodes method 1 as described in the Algebraic Development section; ALGNETZ encodes method 2. The results of both ALGNET programs were compared with the ANALYZER and KYPIPE programs. ANALYZER makes use of the Finite Element Method and incorporates pipe components as separate elements. KYPIPE uses the Linear Theory method and includes pipe components by adding their effect to the pipe elements. Both ANALYZER and KYPIPE programs have been empirically tested and found to be accurate. Several different hydraulic networks were analyzed by these three programs, and the results were found to be very strongly correlated in all cases. As an example, a hydraulic network system is shown in Figure 16 and Figure 17 . Figure 16 depicts how the network would be labelled for analysis by the ANALYZER program, whereas Figure 17 shows the same network labelled for analysis by the ALGNET program. Note here the reduction in the number of nodes to be analyzed, from 12 to 6. A full explanation of how ALGNET data files are set up for this same network is given in Appendix B. The hydraulic head values calculated by ALGNET are compared to the values Page 63 Page 64 'h 31 [“3 12 Figure 16 Network labeled for analysis by ANALYZER Page 65 Ten (:3 9.... H11. LE] Q Figure 17 Network labeled for analysis by ALGNET Page 66 240 220 / I: 200 Egg 0 2‘ 180 / 160 K 14{) I I I I I I I I f’ 140 160 180 200 220 240 ANALYZER Figure 18 Regression plotted. Page 67 K 160 /////( 140 160 180 200 220 240 KWWPE Figure 19 Regression plotted. Page 68 demonstrated here. hedcdcdalicm ANALYZE ALGNET DIFNETI DIFNETZ nodes nodes man mmsrz ANALYZE KYPIPE 12 0 231.00 231.00 231.00 231.00 10 1 212.51 212.51 212.34 212.50 8 2 200.81 200.81 200.57 200.82 6 3 188.80 188.80 188.60 188.83 4 4 176.79 176.79 176.64 176.84 1 5 157.73 157.73 157.74 157.74 ALGNETI-anidale ALGNETz-dependstnvaiahle ANALYZER-WM ANALYZER-WM mm mm Conant 0 Germ 0 StiEncIYEst 0.105699 SUErroIYEsI 0.105625 NW 0.999964 RSquaed 0.999984 new 6 No.0IObsetvsticns 6 DeucesolFteetbm 5 Deg'eesolFteedom 5 X Confident“) Im0632 X Coefficians) 1.0(XJ637 86 End Cost. 003022 St! Err oICost. 0.00022 ALGNETI adspermnwuiabls ALGNETZ-depsmbntvaiable KYPIPE-WW6 KYPIPEaindspendsmvsndale Regesdcnomut ngeubnOquut Comm 0 Consist 0 StIEIroIYEsI 0.023208 suEnoIYEst 0.02227 Rsmered 0.999999 RSquaed 0.999999 Nectweewes'om 6 No.0!0beervat‘cns 6 Dagmadfieedom 5 DegeesoIFnaedom 5 X Coefia'enls) 0.999925 X Coefficienqs) 0.99993 St! Err of Coal. 4.83E-05 SI! Err of Coal. 4646-05 Figure 20 The strong correlation between methods is Page 69 calculated by ANALYZER and KYPIPE in Figure 20. Comparison of the two methods is somewhat complicated by the fact that ANALYZER takes pipe components as separate elements while ALGNET does not. The effect of a component in ALGNET is incorporated into the downstream element; as a result, the ANALYZER nodes chosen for comparison with ALGNET should be nodes situated "upstream" of components. The strong correlation (Figure 20) between a method which accommodates pipe components as separate elements and a method which accommodates pipe components as energy discontinuities within elements indicates that the error introduced as a result of the discontinuities is negligible for practical purposes. B. Evaluation of Solutions with Virtual Nodes The computer code, DIFNETl, incorporates the virtual node concept using Method 1 as defined in the Theoretical Development; DIFNET2 incorporates this concept using Method 2. The network shown in Figure 21, for example, is reduced to the network shown in Figure 22 by incorporation of the virtual node concept. DIFNET takes each lateral to be a single linear element, thus greatly reducing the number of nodes--in this case, from 21 nodes to 9 nodes. A complete explanation of how DIFNET data files are set up for the same network is found in Appendix B. 1. Error introduced by virtual node concept While the number of nodes in a network is greatly decreased by using virtual nodes, error is slightly increased. Results, however, seem still quite acceptable. Figure 23 shows a strong correlation between DIFNET and ALGNET. If greater Page 70 "2 "7 "12 " 20 " 19 " 18 " l7 Emitters Figure 21. Network labeled for analysis by ALGNET. Page 71 Emitters Figure 22. Network labeled for analysis by DIFNET. Page 72 node ALGNET DIFNET 2 residuals 0 100.00 100.00 0.00 Regressoon Output: 1 73.11 71.92 4.20 Constant 0 2 54.10 60.27 6.17 Std Err of Y Est 2.649425 3 37.56 35.51 -2.05 R Squared 0.991964 4 27.64 29.61 . 1.96 No. of Observations 9 5 24.12 21.99 -2.13 Degrees at Freedom 8 6 17.68 18.27 0.59 7 20.87 18.75 -2.12 X Coefficiento) 1.0065 8 15,28 15.55 038 Sod Err at Cool. 0.017858 average at residuals . 0.168767 mm elevation of residuals - 2.512838 an: nemnmmmsnhumnu 0 100.00 100.00 0.00 Regression Output: 1 73.11 71.52 -l.59 Cancun: O 2 54.10 59.99 5.89 Std Err 01 Y Est 2.509788 3 37.56 36.04 - l .52 R Squared 0.992619 4 27.64 30.08 2.44 No. 01 Observations 9 5 24.12 22.65 — l .47 Doomed 01 Freedom 8 6 [7.68 18.84 I .16 7 20.87 19.40 - l .47 X Coefficient“) 11178183 8 15% 16.12 0.84 Std Err of Coot. 0.016917 average of residuals - 0.476142 mndard deviation of residuals - 2.352907 Figure 23 Assessment of Virtual Node technique by comparing DIFNET against ALGNET. Note that results are slightly different for DIFNETl and DIFNET2. Page 73 3:4; 5-— E E c 'C c: g 0.5 O : 1 : : : 1 : i : I. 1 2 .3 4 5 6 Number of Partitions Figure 24. Effect of partitioning the laterals on total error. accu line 1018 isn DUI fur C01 16' of Page 74 accuracy is desired, error can be reduced by subdividing the laterals into two or more linear elements. Figure 24 shows the effect which partitioning the laterals has on the total error. Note that the error in this graph seems to approach an asymptote which is not zero. The fact that total error does not approach zero with an increasing number of partitions is largely due to the error inherent in apprordmating a discrete function (emitters on a lateral) with a continuous function (derivative boundary condition)--refered to in this thesis as continuization error. Another strategy for error reduction, perhaps a little closer to the true nature of flow in pipe networks, is the use of quadratic elements in place of linear elements. This strategy was used by Kelly (1989) and Bralts et al. (1993). C. Reducing a large network to a 9mg, manaflble size What follows is a demonstration of the benefits of the concepts developed in this thesis. Figure 25 shows the system to be analyzed, a medium sized submain unit consisting of 20 laterals and 600 emitters per lateral for a total of 12,000 emitters in an area of 1.37 acres. A one-percent slope goes downhill from the submain.‘ Previous methods would require over 12,000 nodes for analysis of this network. The methods developed in this research require only 80 nodes for analysis, about 0.7 percent of 12,000. (This is achieved by partitioning each lateral into three virtual elements as seen in Figure 25.) A great reduction in computer time and memory requirements is evident as a result of this drastic reduction of nodes. Figure 26 compares the solution to the Backstep solution on the last lateral. Values between virtual nodes are calculated by interpolation. Correlation between the two solutions was calculated Page 75 ?°°9°?°?°9?T9°°°f°T° ll— ’ ooiooeidTooooofiiodood g- Wit. oooéoooof$o¢$¢oo¢9¢d W '-— if wait. -I rigure 25. A large submain unit consisting of 20 laterals spaced 5 ft. apart and 600 emitters per lateral spaced 1 ft. apart. Previous methods would require 12,060 nodes for analysis. Here, only 80 nodes are used. Page 76 30 E20 2: 3.10 0 ::::¥r:::::: O 100 200 300 400 500 600 Emitter Number —Bockstep solution—DIFNET solution I Figure 26. The head calculations along a single lateral are compared; The Backstep solution and DIFNET solution are plotted together. Page 77 using all 600 data points; R-squared = 0.998. D. Stability 1. Methods 1 and 2 Stability is a problem with the solution derived by Method 1 (see Appendix D). This is most likely due to the fact that the energy discontinuities at component nodes are incorporated into the forcing vector. The system of equations is more sensitive to values in the forcing vector than to values in the stiffness matrix. Although stability is a problem with Method 1, the solution has eventually converged in every case tested. Method 2 proves to be the more practical of the two methods. 2. Collective emitter flow in virtual node concept When the flow of several emitters is incorporated into a single element using the virtual node concept, their treatment as a constant, Q, in the difl'erential equation caused instability. Again, this is because, as a constant, their effect ends up all in the forcing vector. If treated as a function of head Gh, however, instability ceases to be a problem because their effect is incorporated into the stiffness matrix. The instability caused by treating collective emitter flow as a constant, Q, was great enough to cause divergence in some cases. Page 78 3. The Newton Raphson Method The computer program, NEWTON, was written to explore the possibilities of using the Newton-Raphson method so as to achieve faster convergence. This proved to be impractical, however, due to instability (see Appendix D). E. Continuity mun—n ments As seen in Figure 15, the intra-element function is discontinuous at node i. The discontinuity, however, does not prohibit solution for the following reasons. The discontinuity exists at a node. The derivative at that node is therefore undefined. The derivative of a linear element without this discontinuity, however, is discontinuous at both nodes i and j. The only difference, therefore, lies in the size of the jump; the jump in both cases is finitely defined as the difi‘erence in slope between the two adjacent elements. In the case of the discontinuity, however, the jump goes to negative infinity and back in the process. Still, the final result in both cases is a finite jump. Continuity of first order, 0‘, is conserved throughout an element because 4H,, is treated as a constant and therefore does not appear in the derivative. Also, continuity of zeroth order, C°, is conserved because head at node j of one element is equal to head at node i of the next element. So in practice, the continuity requirements for solution of a second-order partial differential equation by the Finite Element Method are met. Page 79 F. Discussion For most practical purposes, the accuracy of the DIFNET programs should be sufficient. Where this level of accuracy is not sumcient, the laterals can be broken down into two or more segments, increasing the number of segments to increase accuracy; or, another way to increase accuracy would be to implement higher order shape functions. G. Ideas for further investigation The equations developed in this thesis describe flow through a pipe element in one dimension. The virtual node concept converts a series of discrete flow losses (a series of emitters) to a continuous derivative boundary condition, thereby reducing a lateral with several emitter nodes to one element. Perhaps the virtual node concept could also be applied in two dimensions where a derivative boundary condition is applied to a surface. Or for that matter, why not include the third dimension to incorporate infiltration into the soil, making the model complete. 1. Some ideas for a two dimensional development Consider, for example, a rectangular field with a submain and several laterals as shown in Figure 27 . The thicker lines represent pipes whereas the thin lines represent the boundaries of the rectangular field. The larger dots represent the nodes Page 80 NUDE EMITTERS Figure 27 Sub-plot of irrigated field with nodes numbered for two-dimensional analysis using rectangular Lagrangian element. Page 81 of a two-dimensional rectangular Lagrangian element. The Lagrangian element has nine nodes; its polynomial is second degree. The shape functions for a rectangular Lagrangian element are listed in table IV. These shape functions are given in the Natural Coordinate System (see Segerlind, 1984). Page 82 Table III: Shape functions for nine-node Lagrangian element. Node Shape Function 1 551(5-1) (11-1) 2 «21 (11-1) (9-1) 3 551 (5+1) (11-1) 4 --§ (5+1) (1.2-1) 5 5} (5+1) (n+1) 6 --'21(n+1) (53-1) 7 551 (5-1) mm 8 -—§ (e-l) (113-1) 9 (53-1) (113-1) Page 83 q y+dy/2 i_ J; dy qx-dx/2 -> ~5— qx+dx/2 T l X qy—dy/2 point (x,y) | l "ldxi‘" Figure 28. Differential conservation of mass as continuous function in two dimensions. Page 84 It may be possible to approximate a hydraulic head topology of the surface shown in Figure 27 by considering its discrete function to be continuous. As an extension to this research, a partial differential equation was developed to describe in two dimensions this topology created by an irrigation network. This was achieved by formulating conservation of mass as a continuous function in two dimensions. Initial observations are presented in Figure 28. Conservation of mass then dictates that flow into the control "volume" equal flow out: ' aq aq qx_d7x + qy- £2! a”? qy. 32; —ade —aydy Or, shifting the point (my) to the lower left corner of the control ”volume", a, + q, - am, - cm, - g-gdx - ggdy = o (60) where, q. = ngfi x (61) and as before, qmdx = Dig; x + D‘s—fgdx (52) The same applies in the y-direction. Substituting equations (61) and (62) into equation (60) gives the desired partial differential equation: Page 85 5311 D _ x 6::2 6y 6x dx+D,—a“—’:dy-iqu-—g—3dy=o (53) where '3ng + ggdy can be considered either a constant, Q, or a function of head, Gh: axdx + ggdy " 0 ”wish: or, am 2..., - ab = (”fir where n, = number of emitters per lateral n1: number of laterals L = length of element W = width of element The general form of equation (63) is: 063114.196“): x-gt-S yW-Gh-Q=O (64) (55) (56) In this form, the operators D, and D, resemble conductivity of heat transfer problems. Here they represent the conductivity of pipe flow in a two-dimensional grid. Their Page 86 determination is not straight-forward and is likely the key to solving this problem. Some observations are (1) the conductivity, Dy, looks like it should be the sum of the lateral conductivities divided by the length of the element‘: D = :11 1 (db)‘%.'“ Y dx2—_f d—y (67) 5 a1 and (2) the conductivity, D” depends on the position in y. The second observation becomes obvious with a glance at Figure 29. Water at point (x,y) to get to point (x+dx,y) must first go back to the main through a lateral (distance, y), then through a piece of the main (distance, dx), and then back through a lateral (distance, y). The conductivity of this route is then calculated by adding resistances. The total resistance is equal to the sum of the resistances in each segment: rx=rzy+rm+rzy where r, = resistivity“ in x-direction at any point (x,y) r,0 = resistivity in x-direction at any point (x,0) (in other words, resistivity of main) R, = resistance in y-direction (over length y) Or, in terms of conductivity, 5Note that one over dx is squared. The extra dx comes from the division of dxdy in equation (63). Likewise, the equation for I% has one over dy in it. ‘Resistivity is the inverse of conductivity; resistance is resistivity times length. Here, resistivity in the x-direction is equal to the resistivity of the main plus the resistance in the y- direction (a function of y). Page 87 I Loterols 1 Figure 29 Flow path for water to get from x to C1: in irrigated sub-plot. Page 88 ._1_=_:.L_+_2_.y Dx D by D = 1 i " _1_ .21 dy <58) Dxa Dy where, :(a—xh” (—-1) where, 3,0 = constant, a, for main and D, is as previously defined. Integration of equation (66) times the shape functions is facilitated by recalling that g = [N] (M. So the integration, Zamrggdx xdx becomes, 5812‘ Q) E. Q) T [N] dx{¢‘°’} x 5? m Page 89 The integration of equation (66) looks like: _a___[N1 aw] . aIMT aw] . [De ———dA{O"} +fn an an dAiO"} - [GEMTIMdAWml - [lemma = o -1 -1 (69) The fact that Dll is a function of y makes integration messy. The resulting equations, while quite messy indeed, were found to have four basic forms. This made computer coding of the 9 x 9 stifl‘ness matrix feasible. The results of these integrations are encoded in the computer program LAGRANGE, found in Appendix A. A typical output of this program is plotted in Figure 30. While results are known to be "in the ball park" of the correct values, time did not allow for a thorough development and evaluation of these ideas. Hence, the preliminary results achieved at this point are inconclusive. If continuization error is significant, perhaps similitude modeling techniques could be used to correct for this error. 2. Adding the third dimension A topology of hydraulic head is easily converted to a topology of flow. Described in terms of flow, the two-dimensional analysis could be expanded to include the third dimension, modeling infiltration as well as distribution. Such a model would be based on Darcy’s equation for flow through porous media: Page 90 Hydraulic Topology 200 1 50 1 OO 50 head Figure 30 A hydraulic topology produced by LAGRANGE program Page 91 g= —KVH where, q = flow K = hydraulic conductivity .6}! +6}! +6}! AH—E—{I .3333 82‘ The hydraulic topology created by the irrigation system would then be applied as a boundary condition on the soil surface. The concepts used in the development of the FINDIT computer program (Kunze and Shayya, 1990) would be essential to the development of such a model. A complete three dimensional model would have many obvious benefits in the design and evaluation of an irrigation system, not to mention applications to environmental studies. Figure 31 shows the geometry of a possible three-dimensional finite element grid for modeling distribution and infiltration. Page 92 1 \fl .'. in‘ijg ! L '0 4 1 x I 5 l Infiltration Figure 31 Example of 3-D Finite Element grid in situ. V. Conclusions and Recommendations Conclusions of this research are listed below: The cumulative efi'ect of pipe components in a hydraulic network is significant. Neglecting the efi‘ect of pipe components may introduce error large enough to misguide the design of a microirrigation system. Existing pipe-network programs which include components were found to significantly increase the number of nodes necessary for analysis. A partial differential equation was developed which incorporates pipe components at nodes rather than as separate elements. Galerkin’s weighted residual method was employed in the Finite Element solution of this equation. This solution was found to agree with an algebraic development of the Linear Theory method. The virtual node concept was successfully incorporated in the partial difi'erential equation to further reduce the number of nodes required for analysis. The results of these developments were compared with those of existing pipe- network analysis programs, namely ANALYZER and KYPIPE. Correlation among all methods was found to be very strong. As expected, some error was Page 93 Page 94 introduced by the virtual node concept. While this error was small (certainly negligible), it was found to be further reduced by partitioning the laterals (the "virtual elements") into two or more segments. Some recommendations for further investigation follow: The Ideas for Further Investigation section in Chapter N outlines a possible development of a two-dimensional flow equation which describes a microirrigation system as having a continuous topology of hydraulic head. This development is supplemented by an alternative approach outlined in Appendix E. A thorough evaluation of these ideas was not performed as a part of this research and would be a logical next step. The reader is encouraged to explore and improve upon the possibilities opened up by these developments. Also included in the Ideas for Further Investigation section is a suggestion that one might apply the two-dimensional topology of hydraulic head as a boundary condition to a three-dimensional porous-media flow equation in order to model flow of irrigated water through soil. The benefits of the equations derived in this thesis will only be realized when they are incorporated into a presentable, user-fi'iendly computer program. Appendix Page 95 Appendix A: Computer Code The computer language used is Quick Basic version 4.5. ALGNETl 0 'T T T Program to solve for pipe network (low T T T I DECLARE SUE 8Ccalc (Natrixlt). RHSNatrixl(). NunRCt, BCnodett). knownvall1), NunNodeO) DECLARE SUD Nelton (T1. Nil, Lot. al. kl) DECLARE SUB ConstSortZ (iit(). jjtl). alt). y!(). it. Cq90!(). Cq160!(). qu. NunElent) DECLARE SUB ConstSortl (11.1). 3391). Eli). yll). it. Cq90lt). quOOlil. qu. NunEle-t) DECLARE SUB Cachength (x!(). y!(). alt). 1191). jj‘l). ElLenqthl, it) DECLARE SUB Stats (NunArrayl1). countt. Neanl, Radiant. StanDevl. Ninl. Naxl) DECLARE SUB NatSave (NatriaAlt). Natriasiret) DECLARE SUB AddSqusre (Nit. Lot. valuet. Nattialt). NunNodet) DECLARE SUB NatShou (NstrixAlt). Natriasizet) DECLARE SUB NatAdd (NattixAll). Matrixatt). NatrixCl(). Matrixsiret) DECLARE SUB AddDiagonal (positiont, valuel. Nstrialtl. NunNodet) DECLARE FUNCTION Deternl (NstriaAltl. Nstriasiretl DECLARE SUB DoDst (Ternl. Ni. kt. NatriaAl1), Donet1). ValDeternl. Natriasiret) DECLARE SUB GetReply (Firsts. Lasts, Replys) DECLARE SUB KeyZArr (NunZArrayl(). RceCountt) DECLARE SUD NatInv (HatriRAll). MatrixBli). Matrixsirei, OK A3 INTEGER) DECLARE SUB Natuult (HatrixAli). Matrix81(). HetrixCl(). Matrixsizet) DECLARE SUB Shoe2Arr (NunZArrayl(). RoeCountt. N9. N9. NumColunnst) DECLARE SUE NaitNey (1 6 CL! PRINT ”Program to calculate heads of hydraulic network system" PRINT INPUT "Enter name of element data file: ', elements INPUT "Enter name of nodal coordinate data file: ", nodes INPUT "Enter name of boundary conditions tile ( " THEN INPUT ”Enter number of boundary conditions: '. NumBCt DIN knownvall1NunBC6), acnodettNunact) END IF INPUT ”Enter value for initial head. 810) (m): ". H0! 'T T T start timer T T T Stertrine! - TIMER vsseeeeseseee I 'T T T read data files T T T 0 OPEN elementt FOR INPUT A8 .1 INPUT 01. NunElent it - NumElemt Dru elenttit). 11.11.). 330(10). dial(i§), st(1e), idial(i§). jdial1it) FOR ii - 0 TO NumElemi INPUT .1. elen8118). 119119), jj‘(i§). dial(i§). HN‘(i§), idial(ib). jdiel(i§) NEXT it CLOSE (1) OPEN node8 FOR INPUT A5 .1 INPUT .1. NunNode‘ it - NumNodst DIN nodettit). xzris). ysrit), 2111‘). ctype9(i§). eq1s0111e1. Cq901(i§) FOR it - 0 TO NunNodet 1 INPUT 01. nodettit). xllit). yltit). 21111). ctypot111). Cq1s011111. Cq90l(i§) t CLOSE”) I! be: <> "” THEN OPEN bc$ FOR INPUT AS 01 FOR it - 1 TO NumBC‘ INPUT .1. BCnodettit). knownval!(i§) NEXT 18 CIDSEU) END IF 'T T T dimension arrays T T T DIN kl(NumElem| + l), kel(NunNode9), kstarl(NumElemt + I). H!(NunNodet) DIN NatrixltNunNodet. NumNodet). RHSHatrixltNumNodet. NumcheO) DIN InvertedNatrixl(NumNodet. NunNodet), AnswerNatrixltNumNodet. NunNodet) DIN NunArrayltNunNode‘) Page 96 DIN OltNumEIemt + 1) DIN flouttNumElemt). CompTermltNumElemt) DIN OK AS INTEGER 0 counter\ - 1 'T T T constants T T T HItO) - H01 gl - 32.2 'acc. due to gravity ft/s‘2 pil - 3.1415926540 NumItert - 2 'Number of iterations before including effect of components HNconstl - 4.73 'T T T calculate constant, k T T T 0 FOR it - 0 TO NumElemt CALL Cachengthtalt). yltl. :lt). ii‘t). jj\(). ElLengthl. it) kill.) - HHconeti ' ElLengthl I (HN‘tl‘) “ 1.852 T dial(l§) ‘ 4.87) NEXT i. 'T T T initialize head values T T T 0 FOR it - 1 TO NumNodet Hltit) - NltO) - it NEXT 1‘ I LOCATE 10. 301 COIDR 31. 0: PRINT "-- Converging --": COLOR 7, O :T T T Start iterative procedure T T T Again: :T T T First few iterations do not include effect of components T T T IF counter. <- NumItert THEN FOR it - 0 TO NumElemO IF Hltiittitl) - Hltjjttitl) - 0 THEN kstarltit) - 0 ELSE kstsrltii) - (Asstflltiittitl) - Hltjjttitllll T -(1 - l / 1.652) T kltit) T -(1 / 1.852) END IF NEXT it GOTO skipl END IF 'T T T reset values T T T FOR it - 0 To NumElemt flovttit) - 0 CompTermltit) - 0 NEXT it 'T T T Calculate new values for component head-loss terms T T T FOR it - 0 TO NunNode. IF ctype‘ti‘) - 2 THEN FOR 3! - 0 TO NunElen‘ IF 11‘111) - it AND ultiitrjsl) - H!(jj\tjt)) > 0 runs flouttjt) - l ' 1 - positive flow CALL ConstSortitiitt). jj‘t). alt). ylt). j:. Cg90|(). Cq180!t). qu. NumElemi) r! - qu - s / (g . pi! . 2 - idialtjt) * CompTermltjtl - T! T (H!(iit(jt)) - Hltjjttjtlll I (kltj‘) + T!) CALL NewtontTl, H!tiittjt)). Hltjjttjtl). CompTermltjt). k11j91) 'PRINT "Left: ”: CompTermltjt) 'PRINT "Right: ”: T! T ((Hltiittjtll - Hltjjttjtl) - CompTermlt191) / kl(i\)) T 1. 06 ELSEIF jj‘tjt) - it AND H!(jj§tj§)) ’ Hltiittjt)) > 0 THEN flowttjt) - -l ' -1 - negative flow CALL ConstSortztiittl. jjtl). alt). yltl. 3:. Cq90lt). qulOlt). qu. NumElemO) r1 - qu - s / (gl - pi! ~ 2 - jdialtjt) * CompTermltjt) - (T! T (Hltjj‘tj 9)) - H!(ii§tjt))) / kltji)) / (1 t T! I kltjt)) cart NevtontTl, Hl(jjt(j§)). Hittittjtll. CompTermltjt). k!(j9)) END IF NEXT jt NEXT 19 'T T T calculate new kstar's T T T FOR it - 0 TO NumElemt IF Hltiittitll - Hltjjttitl) - 0 THEN kstarltit) - 0 ELSEIF flouttit) - 1 THEN kstarltit) - (H!(iitti§)) - Hltjjttitl) - CompTerm!(it)) T -(1 - 1 / 1.852) T k!(i\) T -(1 / 1.852) ELSEIF flouttit) - -1 THE EN kstarltit) ' (Hltjj\tit)) - Hltiittitl) - CompTerm!ti!)) T -(1 - 1 / 1.852) T kltit) T -11 / 1.852) ELSE kstarlti‘) ' (ABS(Hl(ii‘(i§)) - H!(jj§(i§)))) ‘ -(1 - l / 1.852) T kltit) T -(1 / 1.852) NEXT 1‘ pkiplr :T T T Calculate new values for emitter constants T T T FOR it - 1 TO NumNodet IF ctypetti‘) - 1 THEN x1 - eq1s0111t) x! - Cq90ltit) Netti.) - k! T (Nitit) - rl(i\)) T x! / H!(i§) END IF NEXT ii 0 'T T T Set Up Global Stiffness Matrix T T T 0 Page 97 FOR it - 1 TO NumNodet FOR 3‘ - 1 TO NumNodet Matrixltit. j!) - 0 NEXT 19 NEXT 18 TECOO!IIIIIOICOIIIIIIIIIOIOOOIODIIIOIICC... 'T Add element matrices T 'OOIOIIIOOIOOIIIOOOI-CIIIOIIIOOIOOOOIOOIOO. FOR it - 1 TO NumElemt valuel - kstarlti‘) Hi9 - jjttia) Lot - iittit) CALL AddSquaretHit. Lot. valuel, Natrixl(), NumNodet) NEXT 1. l 'T Add emitter constants to diagonal of matrix T 'IIOOIIIIIIOIIIOOIIO.ICCIIIDIIOIIIOIIIIIIIIII... 0 FOR it - i To NumNodet IF ctypetti!) - 1 THEN valuel - -l T keltit) CALL AddDiagonaltit. valuel. Matrixlt), NumNodet) ND IF NEXT 18 l 'TTT Add hmainltO) to position 11.1) TTT I value! - -l T kstarltO) CALL AddDiagonaltl, valuel. Nstrialt). NumNodet) F 'T T T set up forcing vector as matrix T T T 0 FOR 3. - 1 TO NumNodet FOR kt - 1 TO NumNode6 RHSNatrialtjt. ht) - 0! NEXT kt _ NEXT 3! F 'T T T include boundary conditions T T T I If ch <> " THEN CALL BCcalctNatrixlt). RHSNatrialt). NumBC9. BCnodett). knovnvallt), NumNodet) END IF I 'T T T first time thru. don't include effect of components T T T I IF countsrt < NumItert THEN COTO skipz 0 'T T T effect of components in forcing vector T T T 0 FOR )9 - 0 TO NumElemt IF flovttjt) - 1 THEN RHSHatrixltiittji). I) - RHSNatrixltiiitjt). l RHSHatrixltjjttjt). l) - RHSNatrixltjjtht). l ELSEIr flow‘tj‘) - '1 THEN RHSHatrixltjjttjt). 1) - RHSNatrixltjjttj‘). l) - CompTermltj\) T kstarltjt) nusnatrixltiietjt), 1) - RHsNatrixltiie(jt). 1) + CompTermltjt) T kstarltjt) v - CompTermltjt) T kstarltjt) + CompTermltjt) T kstarltjt) ‘0 an IF urxr je skipZ: l 'T T T add initial head Boundary Condition to forcing vector T T T 0 RHSNAtrixltl. 1) - RHSNatrlxltl. 1) - kstar!(0) ‘ Hl(0) I ' OIIIIOIOIIIIIIIIOIICO...01...... 'T solve simultaneous equations T ' IIIOIIOOIIIIICOO-ICIIOIOCOIICCOI I CALL HatInvtNatrixlt). InvertedNatrialtl. NumNodet. OR) IF NOT ON THEN BEEP PRINT "Bad input. no solution is possible." END END IF CALL NatNulttInvertedNatrixlt). RHSNatrialt). AnswerHatrix!(). NumNodet) 'OOIOII. cnte - 0 FOR it - 1 TO NumNodet IF ABStAnsverNatrixlti‘. l) - Hlti§)) > .01 THEN GOTO OtraVer END IF NEXT 1. 'T T T stop timer T T T TotalTime! - TIMER - StartTimel 00.090.900.90 FOR it - 1 TO NumNodet Hltit) - AnseerNatrialti‘. 1) NumArrsyltit) - Hltitl NEXT 1‘ I 'T T T calculate coefficient of uniformity T T T CALL StatstNumArraylt), NumNodet. Neanl. Hedianl, StanDevl. Hinl. Max!) Ul - 100 T (l - StanDev! / Mean!) 'T T T print results T T T F CLS : LOCATE I. 1 it - 0 PRINT "Head at Node 1': it: "): "; Hltit) FOR it - 1 TO NumNodet Page 98 PRINT "Read at Node ("r 1‘: ")2 ": H!(i.) NEXT 19 PRINT PRINT "Coefficient of Uniformity - ": U1: ".” PRINT "Converged after ”; counter.: ” iterations." PRINT ”Total time of convergence - ”: TotalTimel: " seconds" PRINT 0 'T T T save data in DIFNET format for comparison T T T 'T T T (save only those points which correspond to DIFNET output T T T 0 biggestyl - 0 FOR 1. - 0 TO NumNode. IF ylti.) > biggesty! THEN biggesty! - ylti.) outfilel - "a:anl_' + NID$(node8. 7. l) + ”.out" OPEN outfile! FOR OUTPUT A3 01 FOR 1. - 0 TO NumNode. IF ylti.) - 0 THEN PRINT .1. Hitl.) ELSEIF ylti.) - biggestyl THEN PRINT .1. Hill.) END IF NEXT 1‘ CLOSE (1) I 'T T T save all data points 7? T T T I INPUT "Save results 7 (y/n) ”. ans: IF UCASEstansS) - ”Y" THEN INPUT ”Enter name of output file: ", filenames OPEN filename’ FOR OUTPUT AS 01 i. - 0 PRINT 01. ”Read at Node 1”: 1.: ”): ": Hlti.) PRINT 01. FOR 1. - 1 TO NumNode. PRINT .1. "Head at Node 1": 1.: "): "; Hlti.) NEXT 1. PRINT 01. PRINT .1. "Coefficient of Uniformity - ”: U1: "." PRINT .1. "Converged after ”: counter.: " iterations." PRINT .1. "Total time of convergence - ": TotalTimel: ” seconds" END IF END OtraVeax FOR i. - 1 TO NumNode. Nlti.) - Answerflatrialti.. 1) NEXT 1. LOCATE 15, 20: PRINT "Number of Iterations: ": counter. counter. - counter. + l GOTO Again SUD AddDiagonal tposition.. valuel. Natrixlt). NumNode.) Natrixltposition., position.) - Hatrialtposition., position.) + valuel END SUE SUD AddSquare (Hi.. Lo.. valuel, Natriattl. NumNode.) 'T T T adds square element matrix to global stiffness matrix T T T 0 NatrixltHi.. Lo.) Natria!(Lo.. Hi.) NatrixltLo.. Lo.) Nattixltflib. Ni.) END SUE NatrixltHi.. Lo.) + valuel Natrix!1Lo.. Hi.) + value! Nattix!(Lo.. Lo.) - valuel NatrialtHi.. Hi.) - valuel SUE DCcalc (Natrixlt). RHSNatrixltl, NumBC.. BCnode.t). knownvallt), NumNode.) I 'T T T include known values in matrices T T T 0 FOR 1. - 1 TO NumBC. son J. - 1 TO NumNode. IF 3. <> BCnode.ti.) THEN RHSHatrixltj.. 1) - RHSNatrixltj.. 1) - Natrix!(j., BCnode.ti.)) T knownvallti.) Hatria!tj.. BCnode.ti.)) - 0 Natrix!tBCnode.(i.), j.) - 0 ELSE RHSHatria!(j.. 1) - Nattix!tj., BCnode.(i.)) T knownvallti.) END IF NEXT 3 NEXT 1. END SUD SUE Cachength (alt). yltl. 211). ii.(), jj.(), ElLengthl, i.) I 'T T T calculate length of element T T T F xtemp! - x1tjj.(i.)) - x!(ii.ti.)) ytemp! - yltjj.(i.)) - y!(ii.(i.)) rtempl - rllji.ti.)) - r!(ii.(i.)) ElLengthl - (atemp! T 2 + ytemp! T 2 + ztempl T 2) T .5 END SUD SUB ConstSortl (ii.(). jj.(). xlt). ylt). 1.. Cq90lt). Cq160!t). qu, NumElem.) O 'T T T sort out which constant applies T T T roe 3s - 0 re NumElem. IF 11.11.) - jj.tj.) rues IF ta!(ii.(j.)) - x!(jj.(i.))) on (yltii.tj.)) - yltjj.ti.))) rssn Page 99 qu - quSOltii.(i.)) EXIT SUB ' 180 - degree constants have preference ELSEIF x!tii.(j.)) - xlt j.(j.)) on x!(ii.(i.)) - xltjj.(i.)) THEN Cq! - Cq90ltii.( .)) ELSEIF ABSltyltii.tj.)) - yltjj.tj\))) I txltii.tj‘)) - xltiJthtll) - 1y1<111¢1Tll - Y‘liiTli‘)ll / (xltii.(i.)) - xltjj.(i.)))) < .5 THEN qu - CgllOltii.ti.)) ELSE qu - CgSOltii.(i.)) END 1? ELSEIF ii.(i.) - 11.111) AND 33.11.) <> jj.(j.) THEN 1r txltjjitjtll - xltjjttitlll on (ylljjtliTll - yiljittit))) THEN qu - Cq180)(ii.(i.)) EXIT SUB ' 180 - degree constants have preference ELSEIF x!(ii.(j.)) - xltjj.tj.)) on xl(ii.(i.)) - x1133111.)) THEN Cq! - Cq90!tii.(i.)) ELSEIF Aasttyltii.tj.)) - y!(jj.(j.))) I (x. (11.13.)) - xlliitljtlll - (ylliiiliill ' Y'lfii‘l1T’l’ / txltii.ti.)) - xl(jj.(i.)))) < 5 rue 3N qu - Cgl.0|tii.(i.)) ens: qu - CqSOltii.(i.)) END It can xr urxr 3. END sue SUD ConstSortz (ii.(). jj.t). alt). ylt). i.. Cq90lll. Cq160!(). qu. NumElem.) l 'T T T sort out which constant applies T T T FOR 3. - 0 TO NumElem. IF 13.11.) - 11.13.) THEN IF xltii.(i.)) - xltjj.( 1.)), OR y!tii.ti.)) - yltjj.(j.)) THEN qu - quSOltjj. ti.) EXIT SUB rtsrxr xltii.tj.)) - alt tjj.(j.)) on xltii.ti.)) - x!(jj.ti.)) THEN Cq! - Cq901tjj.(i.)) rtsrxr A85ttyltii.tj.)) - y!(jj.(j.))) / (x!tii.(j.)) — xl(jj.(j.))) - (Yltii.(i.)) - y!(jj.ti.))) / (x!(ii.ti.)) - a!(jj.(i.)))) < 5 TH ”8 qu - Cq1801tjj.(i.)) ELSE qu - CqSOltjj.(i.)) ELSEIF jj.(i.) - 32.tj.) AND ii.(i.) <> ii.(j.) THEN I: xltii.t s1) - xl(ii.(j.)) on yltii.(i.)) - yltii.(j.)) THEN qu - qu.01(jj.ti.)) EXIT SUD ELSEIF x!(ii.(j.)) - xltjj.tj.)) on xltii.ti.)) - xltjj.ti.)) THEN Cq! - Cq90!tjj.(i.)) stsrxr A881ty!tii.tj.)) - yltjj.(j.))) / (x!tii.tj.)) - xltjj.(j.))) - (yltii.(i.)) - yltjj.ti.))) / (Xltii.(i§)) - X!(jj.ti|)))) < .5 THEN qu - Cq160!(jj.ti.)) ELSE qu - CqSOItjj.(i.)) END IF END IF NEXT 1 END SUD FUNCTION Determl (MatrixAt). Matrixsize.) l 'T T T evaluate determinant of matrix T T T F CONST False - 0 CONST True - NOT False DIN Done.(HstrixSize.) FOR 3. - 1 To NatrixSize. Done.(j.) - False urxr j. Vleeterm! - 0! CALL DoDettll. 0, 1, NatrixAlt), Done.t), ValDeterml. HatrixSize.) Determl - Vleeterml END FUNCTION SUD DoDet (Terml. N.. k., HatrixAlt), Done.t), ValDeterml, MatrixSire.) I 'T T T evaluate determinant of matrix T T T I CONST False - 0 CONST True - NOT False IF k. > HatrixSire. THEN Sign. - -1 IF (N. NOD 2) - 0 THEN Sign. - l ValDeterml - ValDeterml + Sign. T Term! E IF Term! <> 0! THEN N. - 0 FOR j. - HatriXSize. TO 1 STEP -1 IF Done.tj.) THEN N. - N. + 1 ELSE Done.tj.) - True a1! - Term! - NatrixA!(k.. j.) A2. - H. + N. A3. - k. + 1 A4. - Hatrixsise. CALL DoDettAll. A20, A39. HetfiXAlt). Done.(). ValDeterml. A4.) Done.(j.) - False END IF NEXT 3. END IF END IF END SUD Page 100 SUB GetReply (Firsts. Lasts. Reply:) La: - LEFT$(First8. 1) His - LEFT$(Last$. 1) IF L03 > H18 THEN SNAP L03. H13 PRINT USING "Enter reply from a to I": Los: His 00 Reply! - INKEY! LOOP UNTIL (Reply: >- L03) AND (Reply: <- H15) END SUB SUB KeyzArr (NumZArrayll). RouCount.) CONST EndNumber! - 10101 ’Nod. 01 Rowsize. - UBOUND(Num2Array!.1) ColSize. - UBOUND!Nun2Array!. 2) ' PRINT '--- Enter data for 2-dimensional array. ---' 'Nod. 02 ' PRINT '--- Naximun array size -": RouSize.: ” rows. --" 'Hod. .2 ' PRINT '--- There are ": Colsize.: ' columns per row. ---" 'Hod. 02 ' PRINT '--- Enter ": EndNunberl: " to end data entry. ---' 'Mod. 02 IF RovCount. >- RovSiae. THEN PRINT CHR$(7) PRINT "--- RovCount. too large. ---” EXIT SUD END IF DO NHILE RovCount. < RovSize. RowCount. - RovCount. + 1 PRINT "<<< Row number": RovCount.: ">>>' IF RovCount. - Rousize. THEN PRINT ”TTT Last row TTT“ ColCount. - 0 DO WHILE ColCount. < ColSize. ColCount. - ColCount. + 1 PRINT ” Entry 1”: RouCount.: ".": ColCount.: ”):": INPUT ' ', Entry! IF Entry! - EndNunber! THEN IF ColCount. - 1 THEN EXIT DO ELSE PRINT CHR$(7): PRINT "TTT Cannot end now. Please reenter number. TTT” ColCount. - ColCount. - 1 END IF ELSE Nun2Arrayl1RovCount.. ColCount.) - Entry! END IF IDOP IF Entry! - EndNumber! THEN RovCount. - RovCount. - 1 EXIT DO END IF LOOP PRINT '---": RovCount.: 'rovs entered. ---' PRINT "--- Data entry complete —~-” END :03 SUB NatAdd (MatrixAlt). Natrix8!(). NatrixC!(). Matrixsize.) l 'T T T matrix addition subroutine T T T I FOR 1. - I TO Natrixsize. FOR 3. - I TO NatrixSire. NatrixC!(i.. j.) - NatrixAl1i.. j.) + NatrixB!(i.. j.) err 3. NEXT 1. END SUB SUD NatInv (HALFIXAI‘). NatrixD!(). HALFixSiEOQ. OK AS INTEGER) I 'T T T matrix inversion subroutine T T T l CONST ErrorBound! - .0000000010 'Hod. 01 CONST False - 0 CONST True - NOT False DIN NatrixC!(NatrixSize.. Matrixsize.) FOR 1. - 1 TO NatrixSize. FOR 3. - I TO NatrixSize. NatrixC!(i.. j.) - HatrixA!(i.. j.) 1r 1. - 3. THEN Natrixal(i.. j.) - 1! LEE NatrixBl(i.. j.) - on END 1r FORij. -jl TO Natrixsire. . - . NHILE ABS(NatrixC!(i.. 3.)) < ErrorBoundl IF 1. - HatrixSize. THEN OR - False EXIT SUD END I? i. - i. + l NEND FOR k. - I To NatrixSIze. SNAP NatrixC!(i.. x.). NatrixC!(j.. k.) SNAP NatrixB!(i.. k.). NatrixE!(j.. k.) NEXT k. Factor! - l! I HatrixC!(j., j.) FOR k. - I TO Natrixsize. NatrixC!(j.. k.) - Factor! T MatrixC!(j.. k.) NattixB!(j.. k.) - Factor! T NatrixB!(j., k.) NEXT k. FOR N. - I TO Matrixsize. IF H. <> j. THEN Factor! - -HatrixC!(N., j.) FOR I. - I TO NatrixSize. Page 101 HatrixC!(H.. k.) - MatrixC!(H.. k.) 0 Factor! T HatrixC!!j.. k.) HatrixB!(H.. k.) - HatrixB!(H.. k.) 9 Factor! T Hatrix8!(j.. k.) NEXT k. END IF NEXT N. NEXT 1. OK - True END SUB SUE NatNult (NatrixAl!). NatrixB!(). MatrixC!(). Matrixsize.) I ‘T T T matrix multiplication subroutine T T T 0 FOR 1. - 1 TO NatrixSize. FOR 1. - 1 TO NatrixSize. TempSum! - 0! FOR X. - 1 TO Matrixsize. TempSum! - TempSum! + NatrixAl1i.. k.) T Hatrix8!(k.. 3.) NEXT x. HatrixC!(i.. j.) - TempSum! NEXT 3. NEXT 1. END SUB SUB NatSave (MatrixA11). NatrixSize.) OPEN ”a:natrix.dat” FOR OUTPUT AS 01 HetFormatS - ” 00.0““” 'Hod. 01 FOR 1. - 1 TO NatrixSize. FOR 3. - 1 TO HatrixSire. PRINT .1. USING NatFormats: NatrixAl(i.. 3.): NEXT 3. PRINT .1. NEXT 1. PRINT .1. CLOSE (1) END SUB SUE NatShov (NatrixAlt), NatrixSize.) NatFormatS - " 00.0TTTT" FOR 1. - 1 TO NatrixSize. FOR 1. - 1 TO Matrixsize. PRINT USING NatFormats: MatrixA!(i.. j.): NEXT 3. PRINT NEXT 1. PRINT END SUD 'Mod. 01 SUE Newton (T!. Hil. Lol. XI. k!) 0 'T T T solve for deltah 1x! here) using Newton-Raphson method T T T DO F! - T! T ((Hi! - Lo! - x!) I k!) ‘ (2 I 1.052) - x! dF! - 1.0. T T! T ((Hi! - Lo! - x!) / k!) T (2 / 1.052 - 1) T (-1 / k!) - 1 nevxl - x! - Fl / dF! x! - nevxl LOOP UNTIL Fl / dF! < .01 T x! END SUB SUB ShowZArr 1Num2Array!(). RowCount., N., H., NumColumns.) NumRovs. - 20 'Hod. .1 IF (RovCount. < 1) OR (N. < 0) OR (N. < 0) OR (NumColumns. < 1) THEN BEEP PRINT ”--- Parameter error in ShowZArr ---" EXIT SUB END IF CoISize. - UEOUND(Num2Array!. 2) LastElement. - RouCount. T ColSize. PageSize. - NumRovs. T NumColumns. ElementCount. - 0 NumOnPage. - 0 IndexFormat. - "(6.E)" IF N. - 0 THEN NumFormatS - ' 00." T STRINGS(H.. "0") T "““" ELSE NumFormatS - STRING$(N.. "0”) + ".” T STRING$!N.. "0") END IF IF N. - 0 THEN ColNidth. - 27 L8 E Collidth. - N. + N. + 10 ’Hod. 02 ND IF CLS FOR j. - 1 TO RowCount. StrJO - RIGHT$(STR$(j.). LEN!8TR$(1.)) - 1) FOR X. - I TO Colsize. StrKS - RIGHT$(STR$!N.). LEN(STR$(k.)) - 1) RovLoc. - (NumOnPage. \ NumColumns.) o 1 ColLoc. - (NumOnPage. HOD NumColumns.) T ColNidth. + 1 LOCATE RovLoc.. ColLoc. PRINT USING IndexFormatS: StrJS: StrKS: PRINT USING NumFormats: NumZArray!(j.. k.) ’Nod. 02 ElementCount. - ElementCount. + 1 NumOnPaqe. - ElementCount. NOD Pagesize. IF (NumOnPaqe. - 0) OR (ElementCount. - LastElement.) THEN PRINT "-- Press a key to continue --" DO LOOP UNTIL INKEYS <> "' IF ElementCount. <> LastElement. THEN CLS LSE PRINT END IF END IF Page 102 NEXT k. NEXT 3. END sun SUB Stats (NumArray!(), count.. Mean!. Median!. StanDev!. Min!, Max!) IF count. < 1 THEN EXIT SUB FOR 3. - 2 TO count. Temp! - NumArray!(j.) k. - j. - 1 DO NHILE ((Temp! < NumArray!(k.)) AND (N. > 0)) NumArray11k. + 1) - NumArray!(k.) k. - x. - 1 10°F NumArrayl(X. + 1) - Temp! NEXT 3. FOR j. - 1 TO count. ValueSum! - ValueSum! + NumArrayl(j.) SquareSum! - SguareSum! + NumArrayl(j.) T 2 NEXT 3. Min! - NumArray!(1) Max! - NumArray!(count.) IF ((count. + l) \ 2) - count. \ 2 THEN Mid. - count. \ 2 Median! - (NumArray!(Mid.) + NumArray!(Mid. + 1)) / 2! ELSE Median! - NumArray!((count. + 1) \ 2) END IF Mean! - ValueSum! / count. IF count. - 1 THEN StanDev! - 0! ELSE StanDev! - SOR((SguareSum! - count. T Mean! T Mean!) / (count. - 1)) END IF END SUD SUB NaitRey PRINT 'Mod. 01 PRINT ”--- PRESS ANY KEY TO CONTINUE ---" DO WP WHILE INKEYS ' '"' END SUB Page 103 ALGNETZ I ""Programtoaolveforflawinpipenetworks'" DKILARB sun BCcalc Mauro, Memo, Numncs, acnodeso, kmnvaflo. NumNodes) DmstmmmI-w,w, 11,111) DECLARE sun mason: amo, fixo, x10, yto, 115, 01900, Cqmmo, Cq1, NumElerni) DECLARE sun Con-sum (no, uso, x10, yio, 1s, Cq9010, quso1o,Cq!, NumElemS) DECLARE sun CakLength 1x10. ylo, :10, 11110, i750. Eflmgthl, 1s) DECLARE sun 9m (NumAmle, cams, Meanl. Medias-ll, 5mm Min), Max!) om SUB Mafiave NITHXMO. Karim) DELARE SUB AddSquare (HIS, L05, valuel, Man-1x10, NumNodaN) DKZLARE SUB MatShow (MEMO, MatrlehES) DKZLARE SUB W1! W0. Man-1:810, MatrhCIO, MattixSIuN) 0mm; SUB Adleagonal (poduonfi, valuel, MatTiXIO. NumNodai) DKZLARB FUNCTION Deccan! Wan-M10, W5) DKZLARB SUB DoDet (Tamil, MN, 11%, W10, Doneflo, ValDeunnl, MatrixSIni) Om SUB Cahply (Fm. DIS. Reply!) DKZLARE SUB KeyZAn NumZAmyIO, W) DKZLARE SUB Madnv (Man-MO, Man-1x310, MatuxSInfi, 0K AS INTEGER) om SUB MatMult Natl-(NAM. MatrixBIO, MattiaCIO, WW) om SUB ShowZArr (NumZAI-rayio, W, NS, MS, NunColuTnnsS) DKZLARE SUB WW 0 (18 PENT'PTopImtooalmlata heads 1! hydraulic natural: aymm" PRINT INPUI'Enmnameofelemdataflle: ”,aletnen!‘ INPUT ”Enter name at nodal modulate data file: ', nodeS INPUI'EnMnanIofbwndaryoDndiumfllakENIEIanone): ”.136 11' b6 0 "' THEN INPUT 'Enter number :1 baandary audition: ", NumBC$ DIM kmnvallmumxfl. BCnodeMNumKfi) END IF MUI'EnmnhebrWMHw) (m): CHI! neemmeae “nail-m "OIOOOIIOOOI l "“Teaddatafllea'T' OPEN elem-us TOR INPUT AS SI INPUT II. NumElcmi 1% - NumEleTni DIM elemfiufi), 111.0%), ififlfi). dialai), HW‘OE), MINCE). )dhKiS) IOR 1% - 0 TO NumElemi INPUT .1, elemfiufi). UNIS), 1550‘), dlaw‘b), HWNUN), 1111-1115), idialfls) NEXT 1% cwsa (1) OPEN node. TOR INPUT AS .1 INPUT .1. NumNod¢$ 1% - NumNodei DIM mum), x1015), was), 21111.), was), quaaos), 019011115) 10R 1% - 0 TO NumNodeS INPUT II, node‘OS), 110%), ”(1%), 210%). WON), (21180105), C490!(1%) NEXT 1% Page 104 CLOSE (1) IF b6 <> "' THEN OPEN ch RJR INPUT AS '1 ICE i% - 1 TO NumBC% INPUT ll, K1no:k%(i%), kmnvaii(i%) NEXT i% CLOSE (1) END IF I 'T'Tdinmnsionarraya'" DIM ki(NumElem% + I), hiWumNodefl, kstariNunEiem% + 1), Hi(NumNode%) DIM Mauixi(NunNoda%, NumNode%), W!!NumNoda%. NumNod¢%) DIM InvertedMatdxiONumNodd. NumNode%), AnswethbixiWumNodeE NumNode%) DIM NumAmyimumNode” DIM QiNumBJam% + 1). “NW + 1) DIN flow%Nuanlam\%). CompTarnuaNumElemfl DIM OK AS INTEGER mam-I 0 I. O O m 0 O O HIM-HG g! - 32.2 ’am. due to gravity In/sTI pi! -3.14IS92650 Numlm - 5 ’Nunber of iterations baton including effect of components NumIter% - Numlm + 1 Wound-.13 OeeeeeeeeeeemMm'keee ”klfi-OTOW CALLCMMWOI ”0: d0: “‘0: ”‘0; mafia“. 1‘) HONMI mei 'w I (HW‘O‘) O 1.852. dUU‘) " (.37) NEXT“ "'Tinitiaiiaheadvahea'" K)Rl%-1TONumNode% main-Hm.“ NEXT1% mm 10. I): COIDR 31. 0: PRINT '- Converging -": COLOR 7, 0 '"Taaniteaadvepmaduu'" Again: l "T'Puetiawmsdonotmmeihaofmponents'" IF counters < Numltefi THEN 10R 1% - 0 TO NumElem% IF mamas» - Hi(i%(i%)) - 0 THEN was) - 0 ELSE Maria” - (ABSO-Iiai%(i%)) - I-I!(jj%0$)))) “ -(1 - l / 1.852) ' ”(1%) " -(1 / 1.852) END 11' NEXT 1% 0010 skipl END IF MOOmmOOO IORi%-0TONuTnElem% MOM-0 NEXT!% ""Caioalaunewvaheaioroomponent head-bastanns'" RDR 1% - 0 TO NumNode% m was) - 2 THEN ma )1. - 0 To NumElem% n: 1111115) - 1% AND mamas» - waxes» > 0 THEN mom-1 '1 -poailivei|o\v CALI. ConsiSou-t1(ii%0, fi%0, x10, yio, i%, Cq90|0, quBGO, qu, NumElem%) non-qu-a/ (si'pil TZTHwQ'NMQ 010%) - «HURON» - Pumas)» / (”0%) + 116%)» T 5 IP 010%) > 0 TI-iEN Page 105 CALI. Newtonmofl, 111(11%(1%)). H1(]j%(j%)). 030%). 1110‘”) END IF EISEIF ij%(j%) - 1% AND H!(jj%(j%)) - H101%(j%)) > 0 THEN flow%(j%) - -l ' .1 - negattve flow CALL 01113150112050, jj%0, x10, y10, j%, Cq9010, CqISOEO, 01!. Nuchm%) ‘1‘1(j%)- qu ’ 8 / (31’ p11 A 2 ' jdh1(j%) " 4) 0105) - «11101505» . H101%(j%))) / (1110” + 110%)» " 5 IF 010%) > 0 THEN CALL Nmnfl!(j%), I'I1(jj%(j%)). H1050”). 010%). 111115)) END IF END IF NEXT 1% END IF NW 1% '- - - cum m w. - - - 10R 1% - 0 ‘10 Nam 11: mama» - missus» - 0 THEN was» - o EISEIF flow%(1%) - 1 THEN Inc-rum - 1 I was) ' was) A 352 + 110%) ' 010%» m now%(1%) - -1 ‘IHEN WU” - l / (”(1%) ' 010%) A .852 + 110%) ' 0105)) 3153 was) - 01350110111011» . mason») A -(1 - 1 / 13521' was) A 41 I135» END 11? NEXT 17. “cpl: “"Calmhumnhufamm"' ED: 115 - 1 '10 NumNodc% m was) - 1 THEN 1a - Cq1w1(1%) 11 - quaas) h10%)- u - was) - z1(1%)) A :1 / was) END 1? NEXT 1% M-saupcmlsumum-H IO! 1% - I TO NumNodo% R)! 1% - 1 IO Nuu'Noda WK“. 1‘) - 0 NEXT 1% NEXT 1% " Add claim mafia: ’ NR 1% - I TO NumEIcm% vahul - macs) I-I1% - fi%(1%) 1.0% - 11%(1%) CALI. AddSqunnG-Il%, Lo%, valuel, Mun-11110, NumNode%) NEXT 1% "AddmmbdhgomldM' 0 m3 1% . I TO NumNodc% IF dypc%(1%) - 1 THEN valud - -I ' halt“) CALL AddDhgomKfl. valuel, Man-1:110, NumNode%) END IF NEXT 1% "" Add W0) to 96111011 (1,1) "° valuel - ol ‘ WM) CALL AddDhgomKl, valuel, Macaw, NumNodc%) ""utupfcrdngvmuautux'“ 903511-1me¢“ Page 106 NE k% - 1 IO NumNode% Wm10%, 11%) - 01 NEXT 11% NEXT 1% " ' ' 111:de bwnduy common ' ‘ ' IF b6 <> "' THEN CALL Kammxw, W10, NumBC%, BCnode%O, knawnvallo, NumNodc%) END IF ""addInmnlhud BmxduyCdetIontofwdngvectot'" W0. 1) - W10. I) - knuKO) ' 111(0) "tolwlmuhmoqum' CALL WNMIO, Wmlo, NumNod¢%, 010 IF NOT OK THEN BEEP PEINT '30:! Input, no Iduuon b potable.“ END END IF CALL WWmIO, Wan-NO, Wuhan-1:10, NumNode%) I. O O O I 0 0 cat% - 0 NE 1% - I TO NumNodc% IF WHIKI‘, 1) - H1119) > .01 THEN GOIO Ctr-Va END IF NEXT 1% '0 O 0 M m D O O Tot-maul - TIMER - Stuffing! "0.0.0.0.... NE 1% - I TO Nam I-Il(1%) - W10%, 1) NumAmyl0%) - mass) mm as M-mummotuwmcy-H CALL WWyIO, NumNod¢%, Metal, Medhnl, StanDevl, M1111, Max!) UI-IW'fl-W/Munl) "DOW”... CLS:IDCATE8,I 1%-0 PEM'HnddNoch';1%;'):';I-I1(1%) IOE1%-ITONumNode% PEINT'Had ItNodc C'; 1%:"):";I-11(1%) NEXT“ PRINT PEWTWfldcmdUnflomRy-‘;U1;‘%' PEINT‘Conm'pd 31h: ';countcr%;' mun: PEINT'I'oqumofconvmu - ";'1'oul'l'1nu1;' seconds" PEINT l "“uwdmInDIFNEl’mfam‘" ""hwoflypdnbmpofikghmww0"' -o Poms-010W IPle%)>b1WTI—IEN WWW” ENDIF Nmm outfits - "13:12: + MIDSCnodeS, 7, I) + "an" OPEN and!!! ICE OUTPUT AS 01 NE 1% - 0 TO NumNodd IF ”(1%) - 0 m PEINT II, I-I1(1%) EISEIF yKM) - 13W THEN PRINT II, ”0%) END IF Page 107 NET“ m0) ”"uwandahpdnhflu' INPUT “Save vault: 7 (y/n) “, nus IPIXIASBQM) - ”Y'TIEN INPUT “Enact nun: {or wtput file: '; 111cm OPEN 11km RDE OUTPUT AS II 1% - 0 PRINT II, 'Hud at Node 1"; 1%; "): ‘; Hum) PEINTDI, RJE 1% - ITO NumNodc% PRINT II, 'I-Iud at Node (';1%;'):";HI(1%) NEX'T1% PRINT”. PENTII,WdUnlmy-';U1;'%' PEINT II, ‘Convupd after '; mm; " Random! PEINTOL'TDnIUmdmm -';ToulTluu1;'oecondo" c1053 (1) ENDIF END OmVa: TOE 1% - I TO NumNodc% I-I!(1%) - Mouth-NOE I) NEXT 1% LCXIATE 15, mzPEINT'NunMoflmn: "; mum mm - mm +1 0010 A3111! SUB AddDhgoul (podtba%, valuel, Mal-1:110, NumNode%) WME m5) - mmzms. W” + value! END SUB sun Addsqum aim. IA%, valuel, mo, NuaNM) wens. 1.0%) - Matt-111015. 1.0%) + «111.1 muons. Hm - 11.111.111.95, 111%) + mm was, 1.0%) - umzaos, Lo%) - valuel mamas, Hm - mums, Hm - mm END sun SUB BCaIc Man-1:110, ”W10. NumBC%, Knodc%0, knownvaflo, NumNode%) NE 1% - I TO NumK:% ICE 1% - 1 'IO NumNode% IF j% o BCnode%(1%) THEN W10%, 1) - W105 I) - Matux1(j%, BCnodc%(1%)) ' knownval!(1%) mm, BCnod¢%(1%)) - O WWO”, is) - o m W10%, 1) - mm, BCnode%(1%)) ' knownvnl!(1%) END 11' NEXT 1% NEXT 1% END SUB SUB Calling“: 0:10, y10, 210, 11%0, fi%0, Ella-3:111, 1%) W - x1(jj%a%)) - 111(11%(1%)) ytunpl - y1(fi%(1%)) . y1(11%(1%)) W - d(jj%0%)) - 21(11%(1%)) w-W*2+ml*2+aanpl*2)*5 END M SUB ComtSDttI (11%0, i%0, x10, y10, 1%, Cq9010, CqIBOIO, Cq1, NumElem%) TOE 1% - 0 IO NumEch IF 11%(1%) - 11%(1S) THEN IF (x101%(j%)) - x1(i%(1%))) OR 01050”) - y1(ji%0%))) THIN Cq! - qu ”(1509) EXITSUB 'Im-degrummnuhnveprcfenna m x1(11%(j%)) - x1(ij%(1%)) OE x1(11%(1%)) - x1(jj%(1%)) THEN Cq! - Cq901(u%0%)) W ABS((y!(11%(j%)) - y1(ij%(j%))) / (Id(11%(j%)) - x1(jj%(j%))) - 01(11%(1%))- y1(jj%(1%))) I (xl(11%(1%)) - x1(jj%(1%)))) < .5 ITEN qu - cqmmamas» ELSE c4: - cmmax» END 11: EISEIF uses) - 1mm AND fi%(1%) o 1150*) THEN 112 moms» - x1(ij%(1%)))OE (y!(fi%(i%)) - y1(i%(1%)))T1-IEN \ 5 \ § \ \ \ ~ Page 108 Cq! - Cq1801(11%(1%)) DOT SUB ’ 180 - degree comm: have peeferenoe ELSE]1 x1(11%(1%)) - x1(jj%(j%)) OR x1(11%(1%)) - x1(jj%(1%)) THEN qu - Cq901(11%(1%)) ELSEIP A86((y1(11%(1%)) . y1(11%(1%))) / (x1(11%(j%)) - x1(1j%(j%))) - (y!(11%(1%)) - y1(jj%(1%))) / (x1(11%(1%)) - x1(ij%(1%)))) < 5 THEN Cq! - OqIK)!(11%(1%)) qu - 01190101505» END IF END 111 NEXT1% ENDsun SUB Con-190112 (11%0, 1%0, x10, y10,1%, Cq9010, quBflO, Cq1, NumEIem%) 10E 1% - 0 TO NumElem% Q Q ‘ \ h \ \ \ m was) - 11%(1%)'1'HEN 111 110111015» - 11(i%(1%))OE muses» - y1(11%(1%)) m qu - 0113011115015» EXIT sun m 11mm» - «was» on 111011.05» - 11135051111131 qu - Cq901(i%(1%)) m ABSWK11%(1%)) - y1(fl%(1%)))/ mums» - noises») - mamas» - ytwsw» / (x1111%(1%)) . xiqisosn» < .5 “mm c4: - (11110111511111) Cq1 - qumisas» m m 3505) - 11%(1%)AND mas) 0 mos) mm 111 muses» - muses» on y1(11%(1%)) - y1(11%(1%)) THEN Cq! - quw1(11%0%)) EXIT sun m 310505)) - x1(11%(1%)) OE x1(11%(1%)) - x1(jj%(1%)) TIM Cq1 - Cq90|(i%(1%)) m ABS((y1(11%(1‘l» - ”01‘0”” / WONG”) - XKfl‘G‘») - M050”) - y1(jj%(1%))) / (x1(11%(1%)) - x1(11%(1%)))) < .5 THE! qu . Cqu101%0%)) qu - 01190105011» END IF END 11! NEXT1% ENDSUB mm Datum! ammo. Metrbflze%) CONST Pelee . 0 CONST True - NOT Pelee DIN Done%0htfllsm%) ICE 1% - I '10 W% Doae%(1%) - Fake NEXT 1% VelDeterml - (I CALL DoDeKIL 0, 1, Mlh'IxA10, Done%0, VaIDetermL MIMXSlu%) Dew - VelDetu-ml END FUNCTKDN SUB DoDet (Tami. 34%, 11%, MEMO, Done%0, VelDeterml, MetflxSln%) (IJNSTPnlne-O (DNST'True-NOI'Fahe mk%>Methh%II-EN Slum-d IFN%MODZ)-0THENSIgn%-I ValDetemi-VelDetez-znl+515n%"reml ELSE IFTeunIOOITI-EN N% - 0 ICE 1% - MAW TO I STEP -1 IF Done%(1%) THEN N%-N%+I EISE Done%(1%) - True A11 - Term! ' Man-11AM“. 1%) A2% - 51% + N% A3% - 11% + I M% - New CALL DoDet(A11, A2%, A3%, MMXAIO, Done%0, ValDeterml, A4%) Done%(1%) - False END IF NEXT 1% ENDIF Page 109 ENDIF ENDSUE sun Geflleply (Fm has. EeplyS) Ills-WA) HIS-ma,” IFIAS>HBTIENSWAPIA£HB PENTUSIM'Enmanyfm5tok';L¢-fisffl5 IX) EeylyO-W IDOPUN'TILCEeplys>-I.OS)ANDCEep)y$<-I*fl$) ENDSUB MWWle.W%) CONSTEndNuubefl-IOIOI ’ModJI W—UMUNDNWJ) CoISUe%-UBOUND(NumZAmyL2) 'PENT'-EnterdmfoTZ-d1membm)emy.—' 'Mosz 'PEINT'mMquumemyehe-flmfl‘mJ ’ModJZ ’ '—Tbenue ';Co)Sln%;'wh1mperm.—' 'ModJZ ' '-Enter';EndNuabet1;‘toeMdatnemy.—' ’ModJZ IFanCwnt%>-Enwflze%m PEINTCHESU) PEINT'-— Mama too Inge. —" EXITSUB DIDIF wWanCmm%>>" UWnfi-WDENPENT“L¢1DW“ W%-0 WWI-TEECDIOount%40. 1%, 019010, Cq18010, C41. NumElem%) DECLARE SUB Confionl (11%0. jj%0, x10, y10, 1%, Cq9010, Cq18010, qu, NumElem%) om SUB Caldmflh 0‘10, y10, 210. 11%0, i%0, Ellangthl, 1%) COMMON mm m1, 5:, p11 common 51mm mo N...amm... I CLS INPUT “Enact m 01 elem data 111:: ", clemflle‘ INPUT 'Entu mm of nodal data 111:: ', 11011111195 PRINT 'Entu m pal-am " INPUT ' k - '; b1 INPUT " x - '; tel INPUT "Enter 1n111al head, H10). (3.): ", HOI OPEN 01¢!!fo m1! INPUT AS II OPEN 11011211105 1011 INPUT AS '2 INPUT II, NumElcm% DIM clem%(NumElem%). 11%(Nuchm%). fi%(NumE1¢m%), d1a1(NumBem%), HW%(NumElem%), 1d1a1(NumElem%), flhlcNumEkmfl, NumEm1t%(NumE1¢m%) K311 1% - 0 TO Numflenfi INPUT '1, elem%(1%), 11%(1%), jj%(1%), d1al(1%), HW%(1%), 1d1a1(1%), )d1a1(1%), NumEm11%(1%) Nm 1% INPUT '2. NumNodd DIM noddmuwNode%), xlwumNode%). y1(NumNod¢%). 21(NumNode%), dype%(NumNode%), Cq1801(NumNode%). Cq901(NumNode%) FDR 1% - 0 TO NumNodc% INPUT .2. nod¢%(1%). 110%). MG”. 210%). dyp¢%(1%). Cq18010%). Cq901(1%) NEXT 1% CIDSE (1) CLOSE a) mud-I DIM H!(NumNodc%), DI(NumElan%), WINumEIeM) DIM DMatrlxlmumNodd, NumNod¢%), WMleumNode%, NumNod¢%), KMatflxKNumNode%, NumNode%) DIM HnWWumNodc%, NumNodfi), FMatrleWuuNodek NumNode%) DIM anaKNumNode%, NumNode%) DIM al(NumF.1¢a\%), dleumElcm%), d111¢NumElem%) DIM flw%(NumElcm%), Q11NumEl¢m%), deltammuuflem‘b) ""comm‘u 0 111(0) - 1'10! 51 - 32.2 p11 - 11415926540 HWconau - 4.73 IN - 1.852 Page 114 I ""calculataelemfllcngth,dx.andmnatant,a"' 10R 11. - 0 TO NumElcm% CALL Calcbangthwo, M0, 210, 11150, jj%0, ElLengfld, m was) - mum: a10%) - 11me / (stas) A ml - d1a10%) A 437) "1.166 NEXT 1s M-mmmuhudnhmu- NR 1% - 1 TO NumNodc% H10%) - HMO) - 1% NEXT 1% “"atanwerattvepmdun'" Nextltentlon: '0 O O M “1 mm D O I no: 11. - 0 '10 mm an is - 0 1o NumNodev. DMatdx10%. 1%)- o mamas, is) - 0 mm, is) . 0 names, 11;) - 0 NEXT 191 m as "“calculaudh‘" non m - 0 1o NumElem% was) - mums» - racism» NEXT as “"Calmlaudeltah"' 10R 1% - 0 TO NumElem% 1P dype%01%(j%)) - 2 AND dh1(j%) > 0 THEN MW) - I ’ I - poatdve flow CALL Conan-11050, fi%0. x10. y10. 1%. 01900. (3418010. 011. Nummcmfl T1-qu'8/(d'p11‘2'1d1a10fl‘0 CALI. Newtonfl'l, a16%), dh1(j%), dx1Q%). ml. deltahlfi”) m ctypc%(jj%(i%)) - 2 AND dth%) < 0 THEN M09 - -I ' -1 - negaflva flow CALL Conusonzwso, jj%0, x10, y10, is, Cq9010, quauo, cqt, NumElem%) ‘1'1-Cq1'8/(g1'p11‘2'1d1a10%)*1) CALL Nmnm, a1(j%), ABS(<11\10%)). dx1(j%), ml, deltah!(j%)) use ddtah!(j%) - 0 END 117 NEXT j% NOOWQhuD.IOI IO! 1% - 0 TO NumEleufl 010%) - l / a10%) " (I / ml) ' A351(dh10%) - deltahl0%)) / dx10%)) " (1 l m! - 1) NEXT 1% ""CalmlahM‘aandQ'a"' 103 1% - 0 TO Nam 11' NumEmK%0%) > 0 THEN Havd - (I / (m1 + 1)) ' (H!01%0%)) - deltah10%)) + (l - I / (ml 4» 1)) ' H!(jj%0%)) Zavd - (21(11%0%)) + 21(ij%0%))) / 2 MMK1%) - NumEtMt%0%) ' 1:21 ' (Havel - Zave1) " xel / Have! / dx1(1%) 010%) - NumEmlt%0%) ‘ b1 ‘ (Have! - lave!) " u! / dx!(1%) ELSE WK“) - 0 010%) - 0 END IF NEXT 1% "' ' ' Construct Global Stiffness Matt-1x ‘ ' ' l Page 115 '° ’ - Add 012mm conciliation to 0mm: ' ° ° FOR 1s - o m NumElem% DMau-1x101%0%). uses» - DMm-uuuses), uses» + mes) / dx1(1%) DMau-lx1(ij%0%), jjses» - DMatux1(jj%0%). uses» + D10%)/ dx10%) Duauurelses), jj%(1%)) - Dummuses), 11%0%)) . mes) / dx1(1%) omnmses), uses» - Dummses), uses» - 010%) / dx10%) NEXT 1s ""AddclmflmtflbuflmtoM-Matflx"' mk1% - 0 TO NumElem% W1(11%0%), 11%0%)) - MMatrlxKll%0%). 11%0%)) + W10%) ' dxl0%) / 3 563181111631“). ”(1%” - MMatrlx1(i%(1%), i%(1%)) + MM!(1%) ' dx!0%) / 3 menses), fi%(1%)) - MMItrlxl01%0%), jj%0%)) + MM1(1%) ' dx!(1%) I 6 Wl(i%0%), 11%0%)) I MMJWl(jj%0%). 11%0%)) 4* MM1(1%) ' dx1(1%) / 6 NEXT 1% ""adddlammlnultyammmdontofmvedor'" FOR 1s - 0 TO NumElcm% muses). a) - mewses), 0) + D1(1%)' deltah1(1%) / dx10%) Wmses), 0) - FMatrIx1(fi%0%), 0) - mes) ' deltah!0%) / dx1(1%) NEXT 1s ""adddhman-WfloatOM-Maut'" 1011 as - o 10 News mmses), a) - maximises), 0) + mes) - dx10%)' deltah10%) / 3 W(fl%(1%), 0) - mmxmses), 0) + MMIes) ' dx10%) ' deltah10%) / 6 NEXT 1s '0 I I m h a.” O D 0 CALL MathODMau-wo, MMaMxIO, man-1x10, NumNodeM " ‘ ' add boundary conduct: Oman value at node I) ' ' ' NumKZ% - I BCnOdc%0) - 0 knownvalta) - 1110’) CALL MWIO, Nah-1:10. NumBC%, BCnOde%0, knownvallo, NumNOd¢%) '00. I CALL Mauveoumo, Khwio, NumNOde%, 01(%) m NOI‘ OK% THEN BEEP PRINT ”No ”mun pea-able.- END END 11' CALL Mummvzo, memo, analo, NumNodc%) "OOMl‘ymb... CLS R31! 1% - 0 TO NumNod¢% PRINT 1%, anal0%, 0) NEXT 1% :'"d1eckagamatpmlous1tcraum"' 1011 1% - I TO NumNod¢% IF ABSO'H0%) - anal0%, 0)) > .01 “II-LN 1011 is - 1 '10 NumNod¢% 11101:) - IMKfi. 0) NEXT 1% CD10 thmeon END IF Nm 1% PRINT "done“ I "00mm”... 'INPIJT'Savcnmluwln)‘;m-p$ ’IFIXZASEKMPQO'N‘THEN ’INPUI'Entatnamofomputflle: “,Outflld outflleI-‘aadnl_”+ WMJ, I) +‘.Out" Page 116 OPEN outfits IUR OUTPUT AS 83 NR 1% - 0 T0 NumNOde% PRINT '3, an|10%. 0) NEXT 1% CLOSE 9) 'END IF END SUB BCcalc M10, mm, NumBC%, BCnOde%0, knownvallo, NumNOde%) "“mhmuflutohdudaknownvahua'" ICE 1% - I TO NumBC% 103 1% - 0 IO NnmNod0% IF 1% o BCnodc%0%) THEN WK“, 0) - W101», 0) - Mmmqs, BCnOde%0%)) ' knmva110%) Matdx1(1%, Knodc%0%)) - 0 Matflx1m0%). 1‘) - 0 mm, 0) - Matdx1(1%, BCnOde%(1%)) ‘ knownvalKII’.) sun aw etc, we #0. 11%0.1j%0. £11:th 1%) ""ubmflutocalmlau lengthotclam'" W1 - x1(11%0%)) - x101%(1%)) m1 - y1(jj%0%)) - y101%0%)) Mp1 - 11%0%)) - 2101%0%)) Engagihl-(xuu'lp152+yumpl"2+amp1*2)* 5 END SUB SUB Comm (11%0, i%0. x10, y10,1%, Cq9010, Cq18010, qu, NumElem%) ""aubmflutoaouwhldt Wantappuu'“ 10R 1% - 0 TO NumElem% IF 11%(1%) - 11%(1%) THEN 11' (x!01%(i%)) - x1(fi%(1%))) OR (yi01%(1%)) - y1(1j%0%))) “II-{EN th - Cq180101%0%)) EXITSUB 'Iw-depucomnuhavapufamca ELSEIF x101%(i%)) - x!(jj%(1%)) OR x101%0%)) - x1(1j%0%)) THEN Cq! - Cq901111%0%)) EISEIF ABS((yl01%(i%)) - y111i%0%))) / (x!01%(1%)) - x1(1j%(1%))) - 0101%0%)) - y101%0%))) / (x!(11%0%)) - x1(jj%0%)))) < .5 THE! Cq! . Cq1w101%0%)) Cq! - 0190050”) END IF ELSEIF 11%(1%) - 11%(1%) AND 11%0%) <> 11%0%1 THEN 11' 01101509) - 31(11%(1%))) OR (lej%(1%)) - y1(11%(1%))) THEN qu - Cqu101%0%)) EXITSUB ’Iw-degreeeomlm have pnfcnnce ELSEF x101%(j%)) - x1(jj%(1%)) OR x1(11%(1%)) - ”(11%0%)) THEN qu - Cq91l01%0%)) ELSEIF ABS((yl01%(1%)) - y1(ij%(1%))) I (xl(11%(1%)) - x1(11%(j%))) - M01%0%)) - y1(jj%0%))) / (1101%0%)) - xlfij%0%)))) < 5 THEN qu - Cq1w101%0%)) m Cq! - Cq901(fl%0%)) END IF END 11' NEXT 1% END SUB \ ‘ \ \ SJB Cameos-12 (11%0, i%0, x10, y10, 1%, Cq9010, CqIBOlO, qu, NumElem%) ""ubmflutoaonwhlchmantappllu"' 9011 is - 0 TO Nunflems m 11%0%) - uses) THEN I? x101%0%)) - x1(i%(1%)) on y101%(1%)) - y1(1j%(j%)) THEN Cq! - Cq1801(jj%0%)) am sun 5155117 x101%(1%)) - x1(1j%(j%)) OR muses» - muses» m 0;: . Cq901(i%(1%)) £15m ABS((yl01%(1%)) - y111j%(1%))) / memes» - x!(ij%(js))) - ezetses» - y1(ij%0%))) / euuses» - xze'iwwv < 5 “EN Page 117 0;: - quoijses» qu - Cq901(i%(1%)) IF EISEIF uses) - 11%U%) AND uses) 0 uses) THEN 1F muses» - muses» OR y101%0%)) - y1(11%(1%)) THEN Cq1 - Cq1801(1j%(1%)) EXIT sun ELSEIF xxelses» - x1(1j%(1%)) 0R muses» - x1(jj%0%)) THEN qu - quouises» ELSEIF ABS((y101%(i%)) - News») / muses» - noises)» - Moises» - y!(jj‘7o(i%))) / creases» - 11(jj%0%)))) < 5 THEN qu - quw1(11%0%)) Cq1 - Cq901(i%(1%)) END 112 ENDIF NEXT1% \ Q ~ \ END SUB SUB W11 W0, Man-1x310, Man-b100, MatrszIuM ""1111Ub1addkbnauhmt1u"' 3311 as - 0 TO mums EOE 1s - 0 TD MatrIxSlu% W10%, 1s) - MatrIxA!(1%, 1s) + Man-1113115, 1s) m 1s Nm 1s END sun SUB Mann (Memo, MatrixBlO. Mutant 0K AS INTEGER) ""mmvmubmtm"' CONST ErmrBautdl - WI 'Mod. 81 CONST Pale - 0 CONST Tmc - NOT Fall! DIM WWWE W%) EOE 1s - o ‘10 W% EOE 1s - 0 TD Tums muncxes, is) - Names, is) 1F1% - 1% THEN annexes. 1s) - u Man-111310%. 1%) - (1 END IF NEXT 1% NEXT 1% 1011 1% - 0 ID ”613% 1s - is WHILE ABS(Ma1:111C10%, 1%)) < EnorBoundl IF 1% - Man-ban% THEN OK - Falu EXIT SUB END IF 1% - 1% + I WEND 10R 11% - 0 '10 Mata-IxSIn% SWAP memes, ks), luau-1200*, ks) SWAP Man-111810%, 11%). Man-1,1810%, 11%) NEXT 11% Factofl - 11 / MatrtsC!(j%, 1%) 1011 11% - 0 TO MatrIxSIn% MatrGC1(j%, 11%) - Faded ' MaUhC!(i%, 11%) MatruBKfi, 11%) - Factod ' Matrle!(j%, 11%) NEXT 11% RR 111% - 0 'IO MatrleLn% IF 111% <> 1% THEN Faaofl - MatflxC!(m%, 1%) R311 11% - 0 “PO Matrixsm W10“, 11%) - Man-1:001“, 11%) + Factofl ' MatrIxC!(j%, 11%) Man-111811115, 11%) - Manxmens, 11%) + Emu-1 ‘ Man-1.111310%, 11%) NEXT 11% END IF NW 111% NDCT 1% Page 118 OK- True ENDSUB SUB MatMult «guano, Mauxazo, MatrIxCIO, MatrixSlzc%) " ' ' matrlx nultIpUcauon subroutine ' ' ' FOR 1% - 0 IO MatrixSlu% 10R 1% - 0 TO MatrIxSIu% TempSuml - 01 R311 11% - 0 TO Man-1.116122% TempSuml - TempSuI-nl + MatthA10%, 11%) ' Matr1x8101%, 1%) NEXT 11% MacuC10%. 1%) - TempSuml NEXT 1% NET 1% END SUB SUB MatShow NatflxAlo, MatflxSIu%) ""auhtouunnodbphy 111-1111““ him-“NM“ - 'MOdJI poms-omrmmsms mkfi-oma'wmxsms PEINI'USINGMatFm-Ims; MattIxAI(1%.1%); Nm1% PRINT NExrts PRINT ENDSUB SUB Newton (11. Il. 11111, 11111, 1111, deltahl) ""auluoudno toalmlau dehh unlng Newton-Raphson method'" «ham-NI DO Fl-TI'O/a1‘(d111-de11ahl)ldx1)*a/ml)-d¢11a111 dnc-I'TI/ml/al/dxd'a/a1'(dhl-dcltal\1)/dx1)"Q/ml-I)-I mxl-dduhl-Flldm dame-news! ImPUNlnFl/dfi<.01'd¢11a111 ENDSUB SUBWaIdCcy PRINT 'MOdJI PRM‘—PE1~33ANYI 0 THEN mason-1 '1 -po.mv¢aow CALL ComtSonl(11%0, ij%0, x10, y10, is, Cq9oxo, cqmmo, Cq1, NumElzm%) Tram-qua/ (31'p11 Az-mmqvu) CALL Newton(n(j%). 110%), amass), dx1(i%),ml, qe1(j%)) 51.5511: dype%(jj%(j%)) - 2 AND was) < 0 mm nowws) - -1 ' -1 - negative flow CALI. Con-1501121150. ij%0. x10, yIO. is. wow, qusao, qu, NumI-‘Jcm%) non-cw“ (gi'pfl *2'jdh10%)“0 CALL Nmnmqs), mm, dh1(j%), was), m1, qetos» ELSE «1%)-(dh109/ dos) / duos» A (1 I ml) ms) - 0 END 11: NEXT :11 “magnum"- ma 1% - 0 TO NumElem% 010%) - dxl(1%) / was) ° qe1(1%) A (m1 - 1) ° dx10%) + ms) ' qc!(1%)) NEXT 1% ""CalmhuM'o'" NR 1% - 0 TO NumEIcM IF NumEnut%(1%) > 0 mm Havel - 1'1101%(1%)) ’ 011(11%(1%)) + 3 ' HKjfiO'n)» / ‘ WK“) - NumEm11%(1%) ' 1:21 ‘ O-Iavcl - lave!) " u! I Hive! I dx!(1%) 010%) - Mama” ' 1:21 ‘ (Have! - lave!) " xel / dx1(1%) ELSE WK“) - 0 Qlas) - 0 END IF NEXT 1% “"Conmaclobalsuffnunuamx"' Page 121 I " ' ’ Add elemnt contflbutiom to D-Matnx ‘ ' ' RDR 1s - 0 T0 MmEIem% Dumwwses), uses» - DMau'ix1(11%(1%), uses» + D1(1%) / dx1(1%) DMmmejses). uses» - DMunxujjses), uses» + mes) / dx1(1%) Duauwcuses), fi%(1%)) - DMatrlx1(u%(1%), jj%(1%)) - mes) / dx!(1%) DMmuquses), uses» - Dummses), uses» - 010%) / dues) Nm 1s ""AddochmtflbudmtoM—Mnflx'" FDR 1s - 0 T0 NumEIcm% mauses). uses» - mwswses), uses» + Limes) - dx1(1%) / 3 menses), ises» - MMmuises), ises» + MM1(1%) - dx!(1%) / 3 mmuuses). ijses» - Wretses). jj%(1%)) + mes) ' dx1(1%) / 6 Wmuises), uses» - mwmjses), uses» + mes) ' dx1(1%) / 6 NEXT 1s '0 0 O m {a W O I 0 CALL Mdmmlo, 1101111111110, KMm-lxlo, NumNode%) ""hdudcbwnduymndfiommwnnhuaMO'" NumK% - 1 80101150) . 0 knownvwfl) - 111(0) CALL EQKKMIWIO, Eula-1110, NumBC%, BCnode%0, knownvauo, NumNode%) I... I CALL MWO, IGIMO. NumNodd, 0107.) IF NOT OK% THEN BEEP PRINT ”No uoluthn pou1bk.‘ END END IF CALI. MaMultGOIMO. FMmhdo, undo, NumNode%) 'OOOdwhymb... (:18 RM! 1% - 0T0 NumNodd PRINTuIIKM, 0) NEXT1% ""checkagunnpnvbuuwenum'" 10R 1% - 1 TO NumNodd IF ABSGfl0%) . alt-10%, 0)) > .01 THEN RJR 1% - 1 TO NumNode% I-I!(j%) - |m10%, 0) NW j% GOTO Mather-don END IF NEXT 1% PRINT “done“ I 'ODQWMII. outfits - "uan_" + MIDSOIodeflles, 8, I) + ".out“ OPEN 0136116 Kill OUTPUT AS ‘3 MR 1% - 0 TO NumNod¢% PRINT 03, 111.10%, 0) NE“ 1% CLOSE 9) END SUB BCalc Man-1x10, RHSMItrde, NumBC%, Knode%0, knownvallo, NumNode%) ""mchudekmnvfluulnmtrm'" ICE 1% - I TO NumBC% R38 1% - 0 TO NumNode% IF 1% o BCnode%(1%) THEN Wigs, a) - Wmujs, a) - Matr1x1(j%, BCnode%(1%)) ' knownval!(1%) Mumlfifl Knode%0%)) - 0 Page 122 Matt-1:1(BCmde%(1%), 1%) - 0 mees, 0) - Mdrlij%, BCnode%(1%)) ' knownval1(1%) SUB Caldength (x10, y10, :10, 11%0, 11%0, Bungthl,1%) ""eubmuunetocalmhte lemhofelemeut'" new» - x1(jj%(1%)) - x101%(1%)) yumpl - y1(11%(1%)) . y1(11%0%)) m1 - z1(jj%(1%)) - z!(11%(1%)) EMQN-(xhutpl‘2+yump1*2+aanp1*2)“ .5 END SUB SUB Con-150m (11%0, i%0, x10, y10, 1%, Cq9010, qusolo, Cq1, NumElem%) ""munwhkhmuuuapplb'" ICE 1% . 0 TO NumElem% IF 11%(1%) - 11%0%)‘1'HEN IF (x101%Q%)) - 11(i%(1%))) OR 01(11%(1%)) - y1(jj%(1%))) THEN Cq! - qufi)!(11%(1%)) EXIT SUB ‘ 1w - degree content: have preference ELSEIF x1111%(1%)) - x!(11%(1%)) 0R x1(11%(1%)) - x1(1j%(1%)) THEN Cq1 . Cq901(11%(1%)) Elm ABS((yl(l1%(1%)) - y1(ij%(i%)» / (x1(11%(i%)) - x1(1j%(1%))) - 04(11%(1%)) - y1(11%(1%))) / (x!(11%(1%)) - x1(1j%0%)))) < 5 THEN Cq1 - Cq1w1c11%(1%)) Cq1 - anouuses» 112 m uses) - uses) AND uses) 0 WW THEN m clauses» - muses») on muses» - 140115015)» WEN Cq1 - anmwses» EXITSUB 'lw-degnemu have preference ELSEIF menses» - muses» on muses» - uejses» m 0.1 - cqoouuses» m ABS((yi(11%(1%)) - y1(11%(1%))) / (xKusew - muses») - M01950“) - Ylfiiwv» / WWW") ' "(WWW ‘ 5 ““5" C111 - equenses» qu - eqoamses» END IF END 11: NM 1s END sun SUB Confined (11%0, 5%0, x10, y10, 1%, Cq9010, Cq18010, Cq1, NumElem%) ""mmWMMepplb"‘ XOR 1% - 0 TO NumElem% IF 11%0%) - 11%(1%) THEN IF x!(11%(1%)) - x1(i%(1%)) OR y1(11%(1%)) - y1(1j%(1%)) THEN qu - Cq1w1(1j%(1%)) EXIT SUB EISEIP x1(11%(1%)) - x!(1j%(1%)) OR muses» - x!(1j%(1%)) THI-N qu - cq90uises» m ABS((y1(11%(j%)) - y1(jj%(1%))) / muses» . muses)» - (y!(11%(1%)) - y1(jj%(1%))) / muses» - x1fij%(1%)))) < 5 THEN Cq1 - Cq1801(11%(1%)) qu - mousse”) END 1? 81.881? 11%(1%) - ij%(1%) AND uses) 0 uses) mm D: x1(11%(1%)) - muses» on yuuses» - y1(11%(j%))THEN Cq! - Cquejses» Exrr SUB ELSEIP menses» - x1(jj%(j%))OR x1111%(1%)) - x1(11%(1%)) THEN Cq! - Cq901(1i%(1%)) ELSEIF ABS((y!(11%(j%)) - y1(1i%(1%))) / (x!(fl%(j%)) - xuijses») - etetses» - yun'ses») I muses» - x!(jj%(1%)))) < 5 THEN th - cqwoxqjses» ELSE Cq! - Cq901(ji%(1%)) END IF END IF Page 123 END SUB SUB MltAdd (MltflxAIO, MIMXBIO, MatrixCIO, MmSM‘I») " ‘ ' matrix edd1t1on subroutme ' ' ' FOR 1% - 0 TO Matr1x512e% 10R 1% - 0 '10 Metr1xS1ze% MetrhC!(1%, 1%) - Metr1xA1(1%, 1%) + MetrixBl(1%, 1%) ND‘T 1% NWT 1% END SUB SUB Methw W0, Max:310, Wk OK AS INTEGER) ""mtrtxlnvenknmbmtlne'" CONST ErrorBound! - nooooooou 'Mod. 01 CONST Pelee - 0 CONST True - NOT Pike DIM MIMIW, Mettquem Fox 1s - o “no MetrkSlu% RJR 1s - 0 TO MatrixSlze% Machetes, 1s) - MetrtxA!(1%, is) 1? 1s - 1s THEN hummus. 1s) - e ELSE mates. 1%) - 0: END IF NEXT 1s NEXT 1s Eon 1s - o ‘10 Mums 1s - is WHILE Wres, is» < EmrBoundl 117 1s - humans THEN 101! k% - 0 TO Metr1x$12e% SNAP 51111111015, 11%), MetflxCI(1%, 1:1.) SWAP Mltfbt3w%, 11%), MetrIxB1(1%, 11%) NEXT 1t% Pm“ - 11 / MetrbCKfi, 1%) 1m k% - 0 TO Mam MIMNC10%, 11%) - Factor! ' MetrhC!(1%, 11%) Mate-1x811”, 11%) - Fmfl ' MmB:(1%, 1t%) NEXT 1t% R311 121% - 0 TD Metr1x512e% IF m% <> 1% THEN Factor! - -Metr1xC1(m%, 1%) ICE k% - 0 TD Metflxsm‘l MetrGC!(m%. 11%) - Mun-1110015, 11%) + Peder! ‘ MetflxC!(1%, k%) WNW, 11%) - Metrlx310n%, 11%) + Faded ‘ MmBKj‘l, 11%) NEXT k% END 1? NEXT m% NEXT 1% OK - True END SUB SUB MetMult W0, Mauxmo, Men-Milo, Men-12151216) ""uutrixuulttpllaumeubmune'" 10R 1% - 0 TO Metr1xStze% 103 1% - 0 TO Mate-1.612% Tempamtl - 1! 10R 11% - 0 TO Metr1xSlze% Tempfiam! - TempSuml + Matrle!(1%, 1t%) ’ MatflxB1(k%, 1%) NE“ k% MetflxCKlS i5) - TempSuml NEXT 1% NEXT 1% END SUB Page 124 SUB MatShow Men-11AM, MatrhrStze%) "OOdWhymu-u... mm - "uMM - ’Mod. an R)RI%-0108’Metrbz$lu% R311 15 - 0 To a ' Metrthbe% PRINT USING Moms; Mamms, 1%); Nm 115 PRINT NEXT 15 PRINT END sun SUBNewton m, e1, dhl, dxl, m1, qel) ""edveiotqeuetngNewm-Repheoumhod'" gel-m DO F1-e1'qe1"ml'dx1+'fl'qe1*2-dh1 dH-ml'el'dxl'qel“(,ml-l)+2'T!'qe1 newxl-qel-H/dfl (pl-nan! WPUNFEABSG’l/dFD<.OI'ABS(qu ENDSUB suawmxcy PRINT ’ModJl PRM'—PRBSSANYKEY'IOCONTINUB—" Do wormmmm-* BNDSUB Page 125 LAGRANGE I " ° ' program {at common 0! 2-D hydraulic topography mg Lagranghn element 0m SUB MatXCoM M10, Nubian-1:10, kl, MatSlu%) Om SUB Methane (Matteo, Mambo, Matrthtu%) Om SUB Matsub (Matqu, Matrlxbo, W10, Manama” Om SUB Matuult (Memo. Mambo. W10, mans) Dm SUB Madnv (Mata-Inc, Mambo, We, 0K AS MEET!) Om SUB MatAdd (bun-mo, Macaw. W10, Martian” 0m SUB Mr 0. 1:. F10) DKIIARB SUB Ky (I. 13. KDle) DELARE SUB K: (a, b, KDKIO, Dad, Dy1) Dm SUB K3 (8. b. (310) Dm SUB MatShow (Memo, MatrixSlu%) Om SUB BCalc W10, RHSMnuxlO, NumBC%, BCnode%0, knownvallo, NumNodc%) W - 9 DIM 01W, W), NewGIMaQSln%, WW) DIM KDRKMIW, W), NuKDxIOdatShe%, MatSlu%) DIM KDleMatSlze%, MM), NewKDyKMMSluE MatSluS) DIM Humans, MatShe%), NWRW, W) DIM mamas-s. Ma”). KDIannSM, W%). aneKMatSlufi MatSlu%) ”Oiwuau... CLS W'Enfierlengthdelemt: ”,Ll INPUT'BMetwidthoCeleumn: “,Wl 'GO‘IOaklp INPUT'Entetmateedlatex-ala: ',NumLat% W'Enutmntetdemepalaml: “,NumEmM PRM'Encermpat-anuum ' INPUT"k-';ld INPUT'x-";xl INPUT'Enmdlamdmaln: “,MalnDla INPUT 'Enter 11W. coefficient Iceman: ',HWma1n% INPUT 'Entetdlamdlaterab: “,LatDla INPUT "Enter HAN. coefficient for laterals: ", HWlat% lNPUT'Enmmlheadlnfieet: “,Hla) "" mutant: "" p1! - 3.14159 Ween-t1 - 3.027 mm-mmmnz-pu/a Ala!!-LatDla*2'p11/6 FOR 1% - 2 TO MatSlu% ”(1%) - HM) - 1% NEXT 1% ex! - HWanltl / Wailn% " 1.852 ' MalnDla " 1.166) Iyl - HWcoM / 0M“ 5 1.852 ' LatDla “ 1.166) “CalculauhgrangiandeMMue” W b-Ll/Z a-W1/2 CALLKguhcm) CALLKy(a.b.KDy10) CALLPvedoflab.F!0) ""3mb pmdunnmhen'" Page 126 Nextlteranon: Dxl-l /axl*54'(ABSO'Il(1)-1'fl(3))/L!)*~.(6 Dy1-1/ay1*54'(ABSG'fl(1)+wo)-1~11(S)-H1(7))/2/W1)“-4.6 IFH1(9)>0'IHEN Geonn1-Nurdat%'NumEn'ut%'kl'H1(9)"(xl-1)/L1/W1 Qe1--1'NuM'NuM'U'H(9)*x1/Ll/W1 Qd-ol'Gconed ELSE Guam-0 Qel-O Qd-O ENDIP '“CalmlatemvahxeeforKDxmu-lx“ I CALL Kan. b. KDxlo, 0x1, Dy1) "“Calmlate mvahwaforothernutrue“ Kd-Geonetl’a'b/ZZS CALL MatXConeflCHO, NewGlO, Kd, Mafln%) Kd-l /(2‘a) CALL MnXCoMtKDxlO, NewKDxlO. Kd.MatSlu%) Kd-Dyl'b/Oo'dla'b) CALL WWO. NewKDy10, Rd, mans) Kd-O 'Qe1'a'b/9 "'Kd-OwhenuelngCOI) CALLMnXCoMG-TO. NewFlo. Kd. Mafia” "' ' (Q - 0) ’“SolveMau-leqIL‘“ CALL MdeKDalO, NewKDyIO, K010, MatSlu%) CALL MathKDlo, NewGlO, K010, MatSlu%) " ' ' add boundary condltlon Omawn value at node 1) ' ' ' NumBC% - 1 801011611) - 1 knmvwa) - H10) CALL BCalcO .01 THEN 0010 TryAgaln 1? NEW 1% PRINT 'Done“ END TryAgIln: NR 1% - 1 To MISIZQ% 110%) - lluKl%. 1) NEXT 1% 6010 Nexfltenrlon SUB BCalc 04mm, Wan-1:10, NumBC%, BCnode%O, knownvallo, NumNode%) " ' ' um euhrouune evaluatee man-lea to helude known valuee ‘ ’ ' KJR 1% - 1 TO NumBC% H38 1% - 1 '10 NurnNode% 1P 1% o BCnode%(l%)TH'B\1 ELSE Page 127 RI-xsumms, 1) - ”15141me 1) - Matr1x1(1%, BCnode%(1%)) ° knownva11(1%) Matrlx1(1%, BCnode%0%)) - 0 MatrixKBCnode%(1%), 1%) - 0 W10%, 1) - Matrile%, BCnodc%(l%)) ' knownval!(l%) ENDIP NEGfi NEG“ ENDSUB SUB Pvector (a. b, F10) ”"thlamhrwmaetamnetamamfordngvector'" I P10.1)- 1 F10. 1) - 6 p10,”- 1 P1“. 1) - 6 P16.1)- 1 Fm. 1) - 6 F117.1)- 1 P10. 1) - 6 131(9, 1) - 16 10Rl%-1TO9 Roars-2109 mam-o: Name NEXTl% ENDSUB sun K; a. b. 610) ""thhmhwflneeeueonetantalnG-mm‘" Cl“, 1) u 16 CKL 2) - 6 Cl“: 3) - 4 GIG, 6) - '2 G10, 5) - 1 C1“. 6) - -2 GIG, 7) - -6 610. 5) - 8 Cl“, 9) - 6 01(2. 2) - 66 010. 3) - 8 61(2, 6) - 6 010. 5) - -2 (:10. 6) - -16 61(2. 7) - -2 GIG. 8) II 6 GIG, 9) - 32 (310. 3) - l6 GIG. 6) - 3 01(3. 5) - -6 61(3, 6) I- -2 010, 7) - 1 GIG. 8) - -2 G10, 9) - 6 61(6, 6) - 66 Cl“, 5) - 8 Cl“. 6) - 6 01(6. 7) - -2 Cl“, 8) - -16 61(6. 9) - 32 CNS, 5) - 16 CNS. 6) - 8 GIG, 7) - -6 61(5, 3) - -2 616, 9) - 6 (“(6. 6) - 66 01(6. 7) - 8 01(6, 8) - 6 (51(6. 9) - 32 01(7, 7) - 16 61(7, 3) - 6 6117. 9) - 6 (3118. 3) - 66 Page 128 01(8, 9) - 32 61(9, 9) - 256 '“mpyvaluatobottomhalfofnutfixm FOR 195 - 1 To 9 R)R1% - 1% TO 9 cross. m - cm. 1%) NEXT 11; NEXT 1% ENDSUB SUBKxub,KDx10.DxLDY') :"'thbeuhrwt1neevahramer-mau'lx"' 1304MB) Z". calculate c1 m a for Weflmm) “' c1-1/Dx1 c2-2/Dy1 '“Mquaflonewerefaumitocomelnburbaakfom” ""Theaeforuu are calculated below. "° lnIM-6'a'c1"3'd+30'a"2'c1"2'd‘2+50'a"3'c1'c2“3 MQ-w'a’u'd’HaJ’cl*6'1m1e1)+18’a'c1*3'd'bOG(cl) lnlO)-39'a*2'clA2'c2"2'Lm(cl)+36‘a"3’c1'd"3'106(c1) 1110(0-12'a‘6'c2"6'Lm(cl)-3'cl*6'IDC(c1+2'a'c2) lan)--18'a'cl"3'c2'Lm(c1+2'a'c2) lnl(6)--39'a"2’c1*2'd‘2’1061c1+2'a'c2) wm-fi-ua-a-aaa-mcmnan-a) lnl“)--12‘a*6'c2*6'100(c1+2'a'c2) Abram-U 10115-1108 Ame-anmas) ND‘T1% M0)-6'a'c1"2'd+12’a*2'¢1‘d*2+2'a*3'c2"3 MQ-S'cl*3'LOG(c1)+9'a'c1*2'd'wC(cl) 1MG)-6'a"2'c1'd‘Z‘lDle)-3'cl*3'1m(c1+2'a'c2) le-o9'a'c152'd'1mm+2'a'a) lMQ--6'a*2'c1'd‘Z’lDGk1+2'a'c2) Marni-(l 10115-1105 Bforrrl-Bforrrl+1ni(1%) NEXT1% 1M0)-6'a'c1“3'a+18'a*2'c1"2'c262+8'a"3‘c1'c2"3 lMQ-4'a"6'd"6+3'c1“6'Lm(c1)+12'a'c1"3‘c2‘IDC(c1) 1M0)-12'a"2'c1"2'c2‘2'LGKcD-3'c1*6'100(c1+2'a'd) 1M(6)-~12'a'c1*3'c2'1.m(c1+2'a‘c2) lMG)--12'a*2°c1"2‘c2"2’1.06(c1+2‘a’c2) Chm-OI EDEN-1105 Cforui qurrIi+ woe) NEXT1% 1M0)-6'a'c1"3'd+6'a“2'c1*2'd2“2+2'a"3'c1'c2"3-2'a‘6'c2"6 MQ-3‘c1*6'IDC(c1)+6'a'c1“S’Q'IDCKCD WG)-3'a"2'cl*2'd"2'106(c1)-3'c1"6'IDG(c1+2'a'c1) 1M(0--6'a'c1"3'c2'IDG(c1+2‘a'c2) MQ-J'a’W'cl"2'c2‘2'L(X?(c14-2’a'c2) Diem-N R3111%-1T05 DfouM-Dfomi+1r\6(1%) NEXT1% "" Valuo are now calculated for the upper trhngle matrix "‘ XDxla,”--7'Aforui/(72'a*6'b'c2“5) KDR10,2)-Afomi/(9'a“6'b'c2"5) KDI!(1,3)--1'Aforrrl/(72'a"6'b'c2"5) KDx1(1,6)-(cl+2'a'cD'BfornI/G6'a06'b'c2‘5) KDxlO.5)--(cl+I'C2)'BforTM/(72'a"6'b'c2"5) KDxlO,6)-(c1+a'Q'anM/(9'a‘6'b'd‘5) KDX!(1,7)--7'(c1+a'd)'Bforrri/(72‘a“6'b‘c2"5) KDx10,8)-7'(CI+2‘a'd)’BforrM/(36'1“6'6'405) KDxKl,9)--2'(c1+2'a'Q'Bfonrl/(9'a‘6'b'c2‘5) Page 129 KDx!Q.D--2'A10rml/(9'a“6'b'c2“5) KDx!O.3)-Morrri/(9'a"6’b'd"5) KDx1C2,0--2‘(cl+2'a'cb'BfonM/G'a‘6‘b'c2‘5) KDx1Q5)-(c1+a'cfl'BfonM/(9'a‘6'b'd‘5) KDx1a.6)--2'(c1+a'Q'BfoM/(9’a‘6'b'd‘5) KDx1O.7)-(cl+a'cm‘BfonM/(9'a06'b'd05) KDlefl-d'kl+2'a‘c2)'Bfoml/(9'a*6'b'c2"5) KDlefl-2'(c1+2'a'a)'Bfor-ul/(9'a06'b'c255) KDle.3)--7'Abrrri/02'a“6'b'c2“5) KDle,6)-7'(e1+2'a'c2)'Bfor-ui/06'a66'b'c265) KDx1G,5)--7'(c1+a'cD'Bbmi/(fl'a‘6‘b'd‘5) KDxlafl-(cl+a'Q'BW/(9'a‘6'b'd‘5) KDle,7)--1'(c1+a'w‘3bml/(72'a66'b'c265) KDx1G,8)-(e1+2'a'c2'BW/G6'a‘6'b‘d65) KDR!G,9)--2'(c1+2'a'cb'Bbrui/O‘a‘6'b'd05) KDx1(66)--7'Cfom0/(18'a‘6'b'd‘5) KDx1(6,5)-7'c1'Bfomi/(36'a‘6'b’c205) KDx1(6.6)--2'c1'Bforrri/(9'a‘6'b'c2‘5) KDx1(6.7)-c1'Bfarri/66'a*6‘b'd"5) KDxl(L8)--1'CforrrI/(18‘a*6'b'c2“5) KDxl(6,9)-6'Cbrrrl/(9'a"6'b'c2"5) KDle,S)--7'Dfomi/(72'a"6'b'd"5) KDle,6)-Dbmfl/G’a*6'b’d"5) KDle,7)--1'Dforui/(72'a*6'b'd"$) KDx!G,8)-c1'Bfor-ui/Gb'a‘l'b‘d‘S) KDx!G,9)--2'c1'BfomilG'a66'b'd05) KDxl(6.6)--2'Dforui/(9'a06'b'd05) KDxl(6,7)-Dfomi/(9'a“6'b'd“5) KD:1(6.8)--2'c1'Bfor-ui/(9‘a‘6‘b‘c2‘5) KDxl(6,9)-6‘¢1'Bluui/O'a’fl'b'dAS) KDx10,7)--7'Dforu*/02'a*6‘b'd‘5) KD:1(7,8)-7'c1'BfauI/(36'a*6'b'a*5) KDI10,9)--2'e1'Bfoml/(9'a56'b'c265) KDx1(B.B)--7'C£onri/(18'a“6'b'c2*5) KDx1mfl-6'Cbml/O'a66'b'd65) KDx10,9)--8'CforrrI/(9'a*6'b‘c2*5) "“eopyvahaeecolwertrungleofmamm Roma-nos: mus-195109 mos, 1%)-1(Dx1(1%, 11.) NEXT]; NEX'I'1% mason SUBKy(a.h.KDy10) ""thbeubrwflneeetaconeranumDy-M'" KDYKL 1,-2.8 KDy1(1.2)-14 tummy--7 mam-a KDy10,5)--1 KDyta.6)-z KDy10.7)-6 Rowan-.32 KDy!(1.9)--16 KDy10.2)-112 KDy10.3)-16 KDy1(2.6)--16 KDMQ.5)-2 wee-16 KDy1C2.7)-2 KDy1C2.8)--16 mam-.123 mam-23 KDy16.6)--32 wee-4 mama-2 mom-.1 mom-s KDylG,9)--16 KDy116.0-66 KDy1(6.5)--32 KDyl(6.6)--16 Page 130 KDyM. 7) - a 10th s) - .15 Roma 9) - 32 KDy1c5, 5) - 23 KDyKS. 6) - u KDyIG, 7) - .7 KDyKS, s) - a KDyKS, 9) - -16 1(Dy1(6. 6) - 112 Rome, 7) - u KDyKG. a) - -16 W16. 9) - 423 amen-u we. a) - -32 W0. 9) - 46 mm. s) - a nmmw-u we. 9) - 256 '“mpyvahanbhottomhaflofmmx’” Roms-ITO9 10R1%-1%'109 W115.1%)-1o¢>04>o Except for the last column, this data file has the same format as that for ALGNET. Nodal data file This data file has the same format as that for ALGNET Page 140 — Output from DIFNm Head at node (0): Head at node (1): Head at node (2): Head at node (3): Head at node (4): Head at node (5): Head at node (6): Head at node (7): Head at node (8): Head at node (0): Head at node (1): Head at node (2): Head at node (3): Head at node (4): Head at node (5): Head at node (6): Head at node (7): Head at node (8): 100 71.52046 59.9921 1 36.04306 30.07906 22.65207 18.83755 19.40233 16.1 1636 Output from DIFNET2 100 71 .91549 60.27168 35.51284 29.61253 21 .98885 18.26993 18.74747 15.55783 Page 141 Appendix C Calculation of average head along a lateral The shape of the energy grade line along a lateral element may be approximated by the equation developed by Wu and Gitlin (1975): 11,—). = 1 - (1 - 1)” 111-HI L where, Hi = head at node i (upstream node) I-lLj = head at nodej (downstream node) h = head at any point along the lateral element 3 = position on lateral element (local coordinate) L = length of lateral element m = velocity exponent (1.852 for Hazen-Williams, 2 for Darcy-Weisbach) Solving for h gives: 1: = H, - (11,-1.1) [1 -(1 «Er-:1] Page 142 Average head is then solved for by letting u = 1 - i- and integrating from u=0 to u=1: 3‘ 1 ll 9‘“ hdu =11 H,+[1- 11H, 0 n+1 n+1 Page 143 Appendix D A comparison of stability In this section, the stability of each method is evaluated by inspection of its convergence curve. These curves which rapidly approach an asymptote represent methods with rapid convergence. Those which oscillate before converging represent methods which are unstable. A convergence curve for each method follows. CONNIOICO Of laet MC. .L 1 1V .1 1.: ‘ 1 \ ‘ Overehoot IJ I -HWWWWWWW —*Neaseanasenseseeeasaeeeeeaeasa Ider of Iteration Figure 1d Convergence of last node in an irrigation network of 36 nodes, using ALGNET1. Note slight oscillation before convergence; this method is somewhat unstable for larger networks. Page 144 COIV3.XLC culparlaan at calveroence 7 a q l k 11(1)) I [Kn-1) 9 "(n-2)] I 2 Page 1 Figure 2d Effect of averaging values from previous two iterations to calculate linearizing constants for current iteration. Oscillation is not damped much and convergence is slower. Page 145 ANACONV.XLC Analyzer Convergence llnd u. a U 0 t 1 t t : i : : : t t t t t t t L l 1 2 3 4 5 8 7 I 9 10 11 12 13 14 15 18 17 10 19 lueber or Iteratlene Page 1 Figure 3d Convergence of ANALYZER is rapid. Page 146 WNW-WW 1 w ‘ .1 _150 lllllLlllllllllLlllLlJLllLllllllllllllllllllllLLJl lUIUUTUIITIUUIUUUUUUIIIIIIUIIIUIIITIIUIUUIUIIerV Iterations Figure 4d. Convergence of Newton-Raphson method for the last node in a 20-node network. While this method is highly unstable, it did eventually converge in this case. Page 147 Appendix E An alternative development for the two-dimensional flow problem Perhaps it makes sense to look at the two-dimensional flow problem in terms of velocity. Continuity then takes the form, E+ivz+am X By 32 I'0 where velocity in the x-direction is defined as, and a ' 1%). k C's-”D (6 . 87 -2al) The head difference between point (x,y) and (x+dx,y) may then be thought of as dependent on the flow path joining these points: Ah, = ‘3ng = axv:dx + 2ayvy'y or, rearranging and solving for V," so that, Page 148 1 V . i 3}; ZaIEV;y T- x 61,752? axdx de g 1 1 _d£1_ 2axv;y 1"“ 331; .55? max a, 6x axdx 8x3 Substituting the equality, yields, 2 6h (%-1) av: . 1 6h_ 3h '3; "_m 1/‘n' ax dx Derivation of the y-term is straight forward: so that, _[_a_h., <--n 63h mag/In Page 149 The z-term from the continuity equation will represent the "velocity" with which water is leaving the system through the emitters (total flow out of the system divided by the system’s surface area). This term will be taken to he the total flow from the element divided by the element’s area: av: a qu a nonlkh'x '3? A dxdy Again, the resulting equation takes the general form, 6'!) 3'1) D— ‘o-Dy — + G h + = ’82:“ ’dy’ 0 where, 23 ‘71:") D 1 db __ d X was" 6* dx 6h (--1) =a1/n[ay and, G = n‘nlkihrll dxdy n-1 Page 150 or, 12,111ka dxdy n-1 These results were encoded in computer program LAGRANGZ which utilizes the Lagrangian element. Maclaurin expansion of Dx Integration of the Dll term from Chapter IV may be facilitated by replacing D, by its Maclaurin expansion. Consider the expansion: 1 a §xk l-x o D, can be rearranged to take the appropriate form: . on 1 a _on D, dyl 11,) where p Dyy so that, D N D = _59. k x 32:?" This series would then be truncated at an appropriate k so as to approximate Dll with reasonable accuracy. Because this expansion has the constraint, |p| < 1 , the Natural Page 151 Coordinate System must be chosen and DN must be less than or equal to D,. If Dx0 is greater than D,, which is unfortunately the case for conventional irrigation systems, the integration limits must be changed. Page 152 LIST OF REFERENCES Bralts, V.F., D.M. Edwards, and I. Wu. 1987. Drip irrigation design and evaluation based on the statistical uniformity concept. Advances in Irrigation, Vol. 4, Hillel. Academic Press Inc. Bralts, V.F. and L.J. Segerlind. 1985. Finite element analysis of drip irrigation submain units. Transactions of the ASAE 28(3):809-814. Bralts, V.F., S. Kelly, W.H. Shayya and L.J. Segerlind. 1993. Finite element analysis of microirrigation hydraulics using a virtual emitter system. Accepted for publication in the TRANSACTIONS of the ASAE 36(3):?17-725. Dhatt, G. and G. Touzot. 1984. The Finite Element Method Displayed. John Wiley and Sons, New York. Finkel, H. 1982. CRC Handbook of Irrigation Technology, CRC Press Inc., Boca Raton, Florida. Vol. 1, Ch. 8, pp. 171- 193. Haghighi, K., V.F. Bralts and L.J. Segerlind. 1988. Finite element formulation of tee and bend components in hydraulic pipe network analysis. Transactions of the ASAE 31(6): 1750-1758. Haghighi, K., R. Mohtar, V.F. Bralts and L.J. Segerlind. 1992. A linear model for pipe network components. Published in Computers and Electronics in Agriculture, Elsevier Scientific Publishing Co., Amsterdam, The Netherlands 7:301-321. Holy, M. 1981. Irrigation systems and their role in the food crisis. International Commission on Irrigation and Drainage. International Conference for Cooperation in Irrigation, N.D. Gulhati Memorial Lecture. Sept. 1981. Jeppson, R.W. 1976. Analysis of Flgw in Pipe Networks. Ann Arbor Science Publishers, Inc. Ann Arbor, Michigan. Kelly, S.F. 1989. Finite Element Analysis of Drip Irrigation Hydraulics Using Quadratic Elements and A Virtual Emitter System. MS Thesis, Dept. of Agricultural Engineering, Michigan State University. Page 153 Kunze, R.J. and W.H. Shayya. 1990. FINDIT: a new 3- directional infiltration model with graphics. ASAE paper presented at the "1990 International Winter Meeting sponsored by the American Society of Agricultural Engineers," Dec. 18-21, 1990, Chicago, IL. Mohtar, R.H., V.F. Bralts and W.H. Shayya. 1991. A finite element model for the analysis and optimization of pipe networks. TRANSACTIONS of the ASAE 34(2):393-401. Okosun, T.Y. 1993. How Much Longer? All Out Publishing. Holland, Michigan. Power, J.W. 1986. Sharing irrigation know-how with developing countries. .Agricultural Engineering, Sept. 1986, pp. 15-18. Segerlind, L.J. 1984. Applied Finite Element Analysis. John Wiley and Sons. Shayya, W.H., R.H. Mohtar, and V.F. Bralts. 1988. ANALYZER: A computer model for the hydraulic analysis of pipe- networks. Department of Agricultural Engineering, Michigan State University, East Lansing, Michigan. Turner II, B.L. and Peter D. Harrison. 1983. Pulltrouser Swamp. University of Texas Press. Austin, Texas. Villemonte, J.R. 1977. Some basic concepts on flow in branching conduits. Journal of the Hydraulics Division, ASCE 103(HY7):685-697. von Bernuth, R.D. and T. Wilson. 1989. Friction factors for small diameter plastic pipes. J. Hydr. Engrg., ASCE 115(2):183-192. Wiseman, Frederick M. 1989. Agriculture and vegetation dynamics of the Maya collapse in Central Peten. Peabody Library Press, Harvard. Wood, D.J. and C.O. Charles. 1973. Closure of hydraulics network analysis using linear theory. Journal of the Hydraulics Division, ASCE 99(HY11):2129. Wood, D.J. and A.G. Rayes. 1981. Reliability of Algorithms for Pipe Network Analysis. Journal of the Hydraulics Division, ASCE 107(HY10):1145-1161. Wood, D.J. and C.O. Charles. 1972. Hydraulic network analysis using linear theory. Journal of the Hydraulics Division, ASCE 98(HY7):1157—1170. Page 154 Wood, D.J. 1980. Users Manual for computer analysis of flow in pipe networks including extended period simulations. Office of Engineering Continuing Education, University of Kentucky, Lexington, KY. Wu, I.P., and Gitlin, H.M. 1975. Energy gradient line for drip irrigation laterals. Journal of the Irrigation and Drainage Division of the American Society of Civil Engineers. 10(IR4), 321-326. Proceedings Paper no. 11750. Wu, I.P., and Gitlin, H.M. 1974. Design of drip irrigation lines. HAES Technical Bulletin. University of Hawaii, Honolulu, (96) 29. Wu, I., T.A. Howell, and E.A. Hiler. 1979. Hydraulic design of drip irrigation systems. HAES Technical Bulletin 105, University of Hawaii, Honolulu, Hawaii. .ul MICHIGAN smTE UNIV. LIBRARIES IIHI”NIWIIWHWl‘llWllHl HIIHIHH‘NIWW 31293010204919