THE USE OF COMPUTERS IN CIVIL ENGINEERING EDUCATION Harold J. Welch Department of Civil Engineering The University of Michigan June 1, 1962 This material is distributed by The Ford Foundation Project on the Use of Computers in Engineering Education at The University of Michigan. This report appears in the library edition of the Final Report of the Project and is also issued as a separate booklet. Similar "Curriculum Reports" for other engineering disciplines are also available on request.

Table of Contents Page I. Introduction C3 II. The Civil Engineering Department of The University of Michigan: Faculty and Curricula c4 III. Introduction of Civil Engineering Students to the Use of the Computer Cll IV. Specific Uses of Computers in Various Fields of Study C16 V. Conclusions C20 VI. Example Problems Prepared by Project Participants C21 76 Use of the Computer as a Teaching Machine for Plane Analytic Geometry C23 77 Optimum Regulation of a Reservoir for Power Production C35 78 Economical Size of Canal-Feeding Hydro-Power Plant C41 79 Most Probable Number Method for Determining Bacterial Populations C47 80 Iteration of Sum of Principle Stresses at Interior Points of Plane Rings C56 81 Elastic Buckling Load for Columns of Non-Uniform Cross Section C67 82 A MAD Program for Truss Analysis C82 83 Stiffness Factors for a Flat Slab Column C97 84 Vibration of Beams on Spring Supports C104 VII. References C113 -C2 -

USE OF COMPUTERS IN CIVIL ENGINEERING EDUCATION ABSTRACT The Civil Engineering Department of The University of Michigan has taken an active part in computing work at the University with a total of seven faculty members serving as participants in the Project on the Use of Computers in Engineering Education during the past two-year period. In addition, twenty-six civil engineering professors from other universities have participated in Project activities for periods varying from one week to a full semester. Most undergraduate engineering students at the University are introduced to computer work through a required introductory course taught by Computing Center and Mathematics Department personnel. The Civil Engineering Department, however, has elected to introduce its students to computers via a departmental course, the second of the three required surveying bourses. Details of the organization of this course are included as part of this report. Digital and analog computer work has been assigned in more than 25 departmental courses during the past year, giving students an opportunity to gain practice in the application of computers in the solution of their engineering problems. The Civil Engineering curriculum, the areas where computer activity has been strongest, and a sampling of opinion of faculty members as to the effectiveness of computer usage in engineering instruction are included. A selected set of example problems prepared by faculy Project participants is also included. These may be considered as a supplement to the seventy-five example engineering problems, including several related to Civil Engineering subject areas, which have been published previously by the Project. I. INTRODUCTION Civil Engineering, originally named to distinguish it from military engineering, bhas always covered a wide field of engineering practice. Civil engineers plan, design and supervise the construction of roads, railroads, harbors, buildings, tunnels, waterways, bridges, dams, airfields, canals, water supply and sewerage systems, and the many other facilities necessary for public works and industrial development. They plan the conservation, utilization, and control of water resources. They operate in the field of surveying and mapping. In recent years, digital and analog computers have become important tools in many of these areas. With this in mind, the Civil Engineering Department of The University of Michigan has eagerly cooperated with the Ford Foundation Project on the Use of Computers in Engineering Education. Since the Project was initiated, seven members of the Civil Engineering Department at Michigan and twenty-three visiting Civil Engineering faculty members have participated in the Project's program. In addition, other members of the faculty have become familiar with computer -C3 -

philosophy and techniques as a result of weekly luncheon discussions and special lectures sponsored by the Project. It has been found that problem solving by computer methods adds to the depth, scope, and logical understanding of many problems. It is believed that civil engineers of the future will profit from an introduction to computers as part of their undergraduate training. II. THE CIVIL ENGINEERING DEPARTMENT OF THE UNIVERSITY OF MICHIGAN: FACULTY AND CURRICULA The courses which comprise civil engineering may be divided into seven general areas: Construction, Geodetic, Hydraulic, Sanitary, Soil-Mechanics and Foundation, Structural, and Transportation Engineering. These areas are referred to as technical options and represent fields of specialization for the civil engineer. In addition, there are a few basic courses dealing with legal problems, professional conduct, ethics, the scientific method, etc., which are common to all fields of civil engineering. Each of the general areas is described below. Members of the faculty primarily concerned with each field and specific courses are listed. Courses which are required for all civil engineering students are indicated with plus signs; courses in which the presentation is related to or supplemented by computers are indicated with asterisks. A more detailed description is given for those courses involving computer applications. Course numbers are given with each course to provide a convenient means of referring to them later. The number of semester credit hours is given in parentheses after each course. A. Courses Common to All Fields of Civil Engineering Lawrence C. Maugh, Professor of Civil Engineering,and Acting Chairman of the Department of Civil Engineering + CE400 Contracts, Specifications, Professional Conduct, and EngineeringLegal Relationships (2). CE500 Fundamentals of Experimental Research (2). CE501 Legal Aspects of Engineering (3). CE990 Civil Engineering Research (arranged). B. Construction Engineering The methods and techniques of modern construction; fundamental principles of construction applicable to all types of engineering structures; business and legal principles of contracting as applied in the field of construction. Glen L. Alt, Professor of Civil Engineering (Advisor) Frank E. Legg, Jr., Associate Professor of Construction Materials + CE350 Concrete Mixtures (1). CE531 Cost Analysis and Estimating (2). CE532 Construction Methods and Equipment (3). CE533 Estimating Practice (1). + Course required for all Civil Engineering students.

Use of Computers in Civil Engineering Education C. Geodetic Engineering Theory and practical applications of surveying and surveying techniques; theoretical and applied geodesy, figure of the earth, local and extended control surveys; precise measurements and the adjustment of observations; design and execution of municipal surveys, topographic mapping projects, boundary surveys and land subdivision; industrial applications of surveying techniques; research in fields of instrumentation, computation, mapping and photogrammetry, and the problems of land surveying and route location. Ralph Moore Berry, Professor of Geodetic Engineering (Advisor) Harold J. McFarlan, Associate Professor of Geodetic Engineering Harold J. Welch, Assistant Professor of Geodetic Engineering +* CE260 Basic Surveying (3). Use, care, and adjustment of basic surveying instruments; horizontal angle measurement, leveling, taping; circular curves, grades and vertical parabolic curves, profiles, earthwork; computation, significant figures, errors, desk calculators, U. S. Public Land System. +* CE261 Surveying Computations (3). Principles of horizontal control, unified system of geometric computations based on rectangular coordinates and the "Intersection Solution"; logical synthesis of computations for complicated geometrical engineering applications; introduction, applications, and use of electronic computers; principles of surveying astronomy, use of ephemeris and star catalog; elements of photogrammetry. +* CE362 Advanced Surveying Measurements (4). Precise observational methods for triangulation and control levels; use of theodolite and geodetic level, base-line measurements; electronic distance measurements; field assessment of errors and adjustments; astronomical observations for precise azimuth determinations; application of precise control to layout of engineering projects. CE505 Boundary Surveys (3). CE560 Photogrammetry (2). CE561 Geodesy (3). * CE562 Geodetic Field Methods (2). Reconnaissance for geodetic triangulation; special observing methods for first-order horizontal and vertical control; Laplace stations and deflection of the vertical; actual observations, reduction and adjustment of results in an actual field situation. Term paper or project report required of each student. * CE563 Adjustment of Geodetic Measurements (2). Theory of least squares, applications to the adjustment of geodetic observations; arrangement for solution of complex adjustments on electronic computers; actual solution of selected problems. * CE564 Special Problems in Advanced Surveying. Special advanced work can be provided for those who have received credit in Advanced Surveying Measurements (CE362). * CE565 Municipal Surveying (2). Control surveys, methods and adjustments for use in municipal mapping and administration, surveys for streets, utilities, property lines, tax maps, subdivision control and development. * CE960 Geodesy and Surveying Research. Assigned work in geodesy, or other special field in surveying of interest to the student and approved by the professor of geodesy and surveying. * Course involves computer applications. -C5 -

Hydraulic Engineering T'Sr application of the fundamental principles of hydraulics and hydrology to the deve.ir pment of water power, flood control, drainage, and improvement of rivers and harbors, and other hydraulic structures. Laboratory facilities and instruction are offered for students who wish to engage in research work in hydrology and hydraulics that will lead to advanced degrees. Ernest F. Brater, Professor of Hydraulic Engineering (Advisor) Victor L. Streeter, Professor of Hydraulics - CE420 Hydrology (3). +-* CE421 Hydraulics (2). Hydrostatic stability; orifices and weirs; Venturi meters; cavitation; pump characteristics; flow in pipes and fittings; unsteady-uniform flow; steadynonuniform flow. Lecture, laboratory, and computation. * CE522 Hydraulic Transients (2). Introduction to water power development, including selection of type of turbine; storage and pondage; surge in pipe lines; water hammer analysis; digital programming of unsteady flow situations. CE523 Flow in Open Channels (3). * CE524 Advanced Hydraulics (3). Two-dimensional potential flow; the flow net; percolation and hydrostatic uplift; side-channel spillways; boundary-layer; hydraulic similitude; hydraulic models; stilling pools. C E526 Hydraulic Engineering Design (3). Design of hydraulic structures such as diversion dams, head gates, control works, silt traps, siphon spillways, side-channel spillways, earth canals, and other structures involving accelerated flow, backwater, hydraulic jump, sedimentation, and erosion. CE623 Applied Hydromechanics (2). CE825 Seminar in Hydraulic Engineering. CE920 Hydrological Research. * CE921 Hydraulic Engineering Research. Assigned work in hydraulic research; a wide range of matter and method permissible. E. Sanitary Engineering The planning, construction, and operation of water works, sewerage and drainage systems, water-purification plants, and works for the treatment and disposal of city sewage and industrial wastes; improvement and regulation of natural waters for purposes of sanitation; air sanitation; principles and standards of ventilation. Jack A. Borchardt, Professor of Civil Engineering (Advisor) Lloyd L. Kempe, Professor of Civil Engineering and of Chemical Engineering Eugene A. Glysson, Associate Professor of Civil Engineering +* CE480 Water Supply and Treatment (3). Sources of water supply, quality and quantity requirements, design fundamentals of works for development, collection, purification, and distribution of water. +* CE481 Sewerage and Sewage Treatment(2). Requirements of residential and municipal sewerage systems, procedures for the design and construction of sewerage and sewage treatment works. -C6 -

Use of Computers in Civil Engineering Education CE580 Microbiology (3). CE581 Sanitary Chemistry (2-3). CE582 Sanitary Engineering Design (3). CE583 Water Purification and Treatment (3). CE584 Waste Water Treatment and Disposal (3). CE585 Municipal and Industrial Sanitation (3). CE586 Industrial Waste Treatment (2). * CE587 Industrial Bacteriology (3). Lectures and demonstrations to illustrate the application of microbiological principles and techniques to industrial processes. CE680 Microbiology II (2). CE681 Advanced Sanitary Chemistry (3). * CE682 Advanced Sanitary Engineering Design (3). Functional design of sanitary engineering structures and typical plant layouts; drafting room and field studies; preparation of design reports. CE685 Special Problems in Sanitary Engineering I. CE686 Special Problems in Sanitary Engineering II. CE780 Public Water Supply (3). CE880 Sanitary Engineering Seminar (1). CE980 Sanitary Engineering Research. F. Soil Mechanics and Foundation Engineering Soil mechanics; evaluation of soil properties and environmental conditions in foundations of earth-supported structures; mass stability in excavations and subsurface construction. Soil engineering; use of soil characteristics and properties and soil classification in design and construction of highways, railways, airports, and other surface facilities. William S. Housel, Professor of Civil Engineering (Advisor) Robert 0. Goetz, Instructor in Civil Engineering Ulrich W. Stoll, Instructor in Civil Engineering + CE445 Engineering Properties of Soil (3). CE543 Soils in Highway Engineering (2). CE544 Airport Design and Construction (3). CE545 Foundation and Underground Construction (3). CE546 Soil Mechanics Laboratory (1). CE946 Soil Mechanics Research. G. Structural Engineering The theory, design, and construction of structures, such as bridges, buildings, dams, retaining walls, and reservoirs, involving the use of steel, reinforced concrete, and lumber; the testing and utilization of soils in foundations and subsurface construction. Lawrence C. Maugh, Professor of Civil Engineering, and Acting Chairman of the Department of Civil Engineering (Advisor) Glen V. Berg, Professor of Civil Engineering Bruce G. Johnston, Professor of Structural Engineering Leo M. Legatski, Professor of Civil Engineering Robert B. Harris, Associate Professor of Civil Engineering Wadi S. Rumman, Assistant Professor of Civil Engineering -C7 -

+* CE312 Theory of Structures (3). Calculations of reactions, shears, and bending moments in simple, restrained, and continuous beams due to fixed and moving loads; simple trusses with fixed and moving loads; determinate frames; columns; tension members; girders; introduction to design. +* CE313 Elementary Design of Structures (3). Design and details of simple beams, girders, columns, and trusses. Computations, drawing work, and laboratory experiments. +* CE415 Reinforced Concrete (3). Properties of materials; stress analysis and design of reinforced concrete structures; introduction to prestressed concrete and ultimate strength analysis. Lectures, problems, and laboratory. CE310 Structural Drafting (2). CE511 Timber Construction (1). * CE512 Advanced Theory of Structures (3). Stresses in subdivided panel trusses; principle of virtual displacements and virtual work; energy theorems; graphical methods; analysis of statically indeterminate trusses and frames. * CE513 Design of Structures (3). Design of reinforced concrete and steel structures. Computations and drawing. * CE514 Rigid Frame Structures (3). Analysis of rigid frames by methods of successive approximations and slope deflections; special problems in the design of continuous frames. CE515 Prestressed Reinforced Concrete (2). * CE516 Advanced Design of Structures (3). Functional design of buildings; selection and analysis of structural elements; reinforced concrete flat slab design; digital computer applications. Lectures and computation laboratory. CE517 Bridge Engineering and Design (3). * CE611 Structural Dynamics (3). Structural vibrations. Transient and steady state response to dynamic forces. Response beyond the elastic range. Response to nuclear explosions. Earthquake forces. Structural response to earthquake. The response spectrum. Seismic building codes and their relation to structural dynamics. CE612 Structural Members (3). * CE613 Structural Plate Analysis (2). Stress analysis of flat plates loaded either in their plane or in bending. Numerical analysis. Applications to special problems in flat slab construction. * CE614 Advanced Problems in Statically Indeterminate Structures (3). Continuous truss bents; hinged and fixed arches; rings; frames with curved members; flexible members including suspension bridges; frames with semirigid connections. CE615 Analysis and Design of Folded Plates, Domes and Shells (3). CE616 Plastic Analysis and Design of Frames (2). CE617 Mechanical Methods of Stress Analysis (1). * CE618 Computer Analysis of Structures (3). The analysis of beams, frames, trusses, and arches by high speed digital computers. The general method of influence coefficients; matrix methods, algorithm development, flow charts, programming. Students will solve a sequence of problems on the high speed computer in the University Computing Center. CE810 Structural Engineering Seminar (1). CE910 Structural Engineering Research. -C8 -

Use of Computers in Civil Engineering Education H. Transportation Engineering Transportation Engineering is further subdivided into Highway, Railroad and Traffic Engineering. Highway Engineering involves the location, design, construction and maintenance of various types of roads and streets, including materials, surveys, plans, specifications, economics, and financing. Railroad Engineering involves the design, construction, and operation of railroad properties, including metropolitan terminals, statistical analysis of operating data, freight and passenger traffic, economics, financing, administration, and regulation. Traffic Engineering involves methods of increasing the efficiency and safety of traffic movement; street and off-street parking; traffic surveys, geometrical design of urban and rural highways, traffic control devices, and other means of regulating and controlling the use of highways. John C. Kohl, Professor of Civil Engineering and Director of the Transportation Institute (Advisor for Traffic Engineering) Ward K. Parr, Associate Professor of Highway Engineering Bruce D. Greenshields, Lecturer in Transportation Engineering Donald N. Cortright, Associate Professor of Civil Engineering (Advisor for Highway Engineering) Clinton L. Heimbach, Lecturer in Civil Engineering (Advisor for Railroad Engineering) + CE370 Transportation Engineering (3). CE550 Highway Materials (3). CE551 Physical Properties of Concrete Masonry (2). CE552 Bituminous Materials and Pavements (2). CE570 Highway Traffic (2). CE571 Traffic Engineering (3). CE572 Highway Economics (2). CE573 Highway Design (3). CE574 Railroad Engineering (3). * CE575 Terminal Design (3). Design of railroad, highway, waterways, and airport terminals, joint terminals, layout of the various types of yards and traffic facilities. CE576 Economics of Railroad Construction and Operation (2). CE652 Advanced Bituminous Materials and Flexible Pavement Design. * CE670 Transportation Planning (3). Analysis of supply and demand for transportation services, transport relationships to land use and other elements of regional and urban planning, and planning techniques applied to transportation problems. CE671 Advanced Highway Engineering (2). CE672 Transportation (3). CE673 Highway Transport (2). CE674 Industrial Transport Management (4). CE970 Highway Engineering Research. CE971 Transportation Engineering Research. -C9 -

A list of courses taken by Civil Engineering students is given in the table below. The table is arranged to indicate the semester in which the course is normally taken. Courses in which the presentation is related to or supplemented by computers are indicated by an asterisk (*). LIST OF COURSES TAKEN BY CIVIL ENGINEERING STUDENTS Total Hours Semester Course Name and Number Hours for Semester Analytic Geometry and Calculus I (Math 233) 4 General and Inorganic Chemistry (Chem 104) 4 1 Engineering Drawing (Engr. Graphics 101) 3 15 Freshman Composition I (English 111) 3 Introductory Speech (English 121) 1 Analytic Geometry and Calculus II (Math 234) 4 Descriptive Geometry (Engr. Graphics 102) 3 2 Freshman Composition II (English 112) 2 15 English Elective 2 General and Inorganic Chemistry (Chem 106) 4 Analytic Geometry and Calculus III (Math 371) 4 Mechanics, Sound, and Heat (Physics 145)5 3 Statics (Engr. Mech. 208) 3 16 Basic Surveying (CE260) 3 Concrete Mixtures (CE350) 1 Analytic Geometry and Calculus IV (Math 372) 4 Electricity and Light (Physics 146) 5 4 Mechanics of Material (Engr. Mech. 210) 4 17 Laboratory in Mechanics of Material (Engr. Mech. 212) 1 * Surveying Computations (CE261) 3 Summer * Advanced Surveying Me<- -urements (CE362) 4 Geology for Engineert.- Gieology 218) 4 Fluid Mechanics (Eng. i.-ich. 324) 3 * Theory of Structures::13) 3 5 Transportation Engirt - (CE370) 3 17 Hydrology (CE420) 3 (*) Advanced Mathematics:ive 3 Elective in Humanities...: Social Science 2 Dynamics (Engr. Mech. 343) 3 * Reinforced Concrete (CE415) 3 6 * Water Supply and Treatment (CE480) 3 17 * Hydraulics (CE421) 2 Engineering Properties of Soil (CE445) 3 Elective in Humanities or Social Science 3 Specifications and Contracts (CE400) 2 * Elementary Design of Structures (CE313) 3 7 Electrical Apparatus and Circuits (EE315)4 16 * Sewerage and Sewage Treatment (CE481) 2 English Elective 2 (*) Technical Elective 3 Thermodynamics and Heat Transfer (Mech. Engr. 333) 4 (*) Technical Electives 7 8 Modern Economic Society (Econ. 401) 3 17 Elective in Humanities or Social Science 3 * Course Related to or supplemented by computers. (*) Student may elect course related to or supplemented by computers. -C1C

Use of Computers in Civil Engineering Education III. INTRODUCTION OF CIVIL ENGINEERING STUDENTS TO THE USE OF THE COMPUTER Surveying Computations (CE261) has been selected as the medium for introducing Civil Engineering students to the use of the digital computer. This course is the second in a series of three required courses in surveying and is normally taken in the second semester of the sophomore year. Civil Engineering students are not required to take the one-hour course in introductory computing techniques (Math 373) which is required of students in most of the other fields of engineering. Surveying Computations has, therefore, been designed as an effective substitute. A. Preparation of Input for the Computer The first half of the semester is devoted to applications of analytic geometry for plane, rectangular coordinate computations. Emphasis is placed on five basic problem types. These types, described in detail in Example Problem No. 76, have the following names: Forward Computation, Inverse Computation, Line-Line Intersection, Line-Circle Intersection, and Circle-Circle Intersection. The students are shown the recommended computing methods, supported by a mathematical derivation, for each problem type. After studying the assignment, listening to the lecture, and assimilating the methods taught, each student is assigned a set of fifteen individual problems of each type. The student is instructed initially to perform the necessary calculations for only the first five problems in each set and to prepare both the given information and the solution on punched IBM cards. These cards are then submitted as data to a computer program prepared specifically for the course,called the "teaching machine program." This program, described in more detail in Example Problem No. 76, has access to all of the given data and computes the correct results for each of the students' problems. It compares all of the items on each of the students' cards with the known, correct values. If any single item in any one of the five solutions is found to be in error, the whole problem is considered incorrect. For each incorrect problem the student is required to correct all of the errors, and also solve an additional problem in the assignment. The original five problems (with errors corrected) plus the additional problems are then submitted in a second approach to the computer. Thus, for the second approach, the student may be required to enter from six to ten problems. The same rules apply again, for a third and final approach to the computer. The worst possible penalty would be that on the final approach to the computer the student would be required to present solutions to all fifteen of the assigned problems. If all of the problems submitted in the final approach to the computer are not correct, the student is penalized in his grade. On the other hand, if a correct solution is obtained -Cll

to all of the problems submitted on any given approach to the computer, the student is given full credit regardless of whether success is achieved on the first, second or third trial. It may be seen that the students quickly realize that it is to their advantage to obtain the correct results for submission on the first pass through the computer. In addition to the savings in their calculating time (on the desk calculator), by solving the first five problems of each type correctly and finishing the assignment, it has been noticed that a friendly spirit of rivalry generally develops between the students in which those who do not complete the assignment on the first pass receive frendly, but reslightly barbed comments on their proficiency. The "teaching machine" theory seems to be working quite well for these problems. The students quickly overcome their difficulties and by the time they have finished the five problem types, they show a marked proficiency in the use of trigonometric tables and the desk calculator. In order to meet the requirements for tolerance in the teaching machine program, the students learn by experience the necessity for using significant figures properly in their computations. For example, they are shown, by actually working problems both ways, that in an intersection of a line and a circle, it is necessary to interpolate fractions of seconds in order to obtain enough accuracy when the intersection angle is very small and the radius of the circle is very large. These are things which cannot be easily taught in classroom lectures but which can be illustrated very pointedly to the student during his laboratory work. This three-semester-hour course is handled with a one-hour lecture per week and two three-hour laboratory periods per week. The first hour of each period is devoted to a lecture, followed by two hours of actual computation. During this time, the students solve problems in the laboratory room using desk calculators. They are encouraged to ask the instructor for help at any time on problems involving the theory of the computations, the mechanics of the desk calculators, the manipulation of trigonometric tables, or an explanation of lack of accurate checking of their completed problems. The "teaching machine" idea serves also to introduce the students to the digital computer. They learn to use the keypunch to punch their own cards before they are actually assigned to write a program for the digital computer. They are strongly impressed with the capabilities of the computer when they see how it is able to solve their problems and check them at a rate of approximately two of their problems per second (on the IBM 709). They realize quite painfully that these problems cost them sometimes as much as an hour or more of desk calculator work. Because of this, when they are later permitted to begin programming for the computer, they do not have the reaction so often expressed by students who are given elementary problems, "I could have worked the problem faster on a desk calculator." -C12 -

Use of Computers in Civil Engineering Education B. Geometric Analysis Using the Interpretive Language, AGILE The writer has developed a computer program over the last several years which has recently been named "AGILE" (Analytic Geometry Interpretive Language for Engineers). A number of small programs to solve problems such as the five basic types described in Example Problem No. 76, have been combined with a reading, storing, and scheduling master program. In this way, any geometry problem (which may consist of a number of basic problems) can be solved with one pass through the computer. After learning to compute the basic problems already described, the students are taught to recognize them when they occur in such applications as the design of a subdivision or a highway interchange. The students are also taught to describe these problems in a logical manner, using the AGILE language. A complete description of the AGILE language is being prepared. If sufficient interest develops, a manual describing the language and its use will be made available. C. Computer Programming Beginning with the first Wednesday night of the semester, a series of three weekly lectures of two hours each, on computer programming, is presented for the entire College of Engineering. These lectures were presented in 1961-1962 by Mr. Brice Carnahan, Assistant Director of the Project on Computers in Engineering Education. The lectures introduce the students to the philosophy of computers, and specifically, to the MAD language. The students are not required to attend these lectures, but are strongly urged to do so, being told that no substitute for them will be presented in classroom lectures. Beginning immediately after the series of night lectures, the students are assigned their first problem which requires them to write a computer program for the course. This initial program is for a forward computation (the first problem type listed in Example Problem No. 76 ), and is handed to the students completely checked out and ready to run on the computer. They are required simply to prepare cards for communicating with the computer. They are also asked to punch data cards for a problem which is given to them. The form of the data cards is slightly different than the form in which the data are presented. In addition to the punched cards, the students are required to submit a flow diagram for the problem and a short discussion of the problem, its solution, and the method which a user must follow in preparing data for a run on the computer. By this means the students are taken over a second, short step towards computer programming by being introduced as painlessly as possible to the computer itself. The submission of a program which has already been written for them may sound like a wasted step but it has been found that this is a step which, when combined with other complications, can become a very frustrating experience to the neophyte. -C13 -

This first program is explained thoroughly to the students so that they understand what is done by the computer as a result of each instruction in the program. Two very elementary loops are included so that the student becomes familiar very quickly and painlessly with the concept of iterative techniques on the computer. The students are permitted to submit this program to the computer only once. It is felt that since everything has been handed to them they should be graded on their ability to follow instructions. The second computer program is for the inverse computation. In this case, they are given a clear concept of the theory and the method of solution of the problem but are required to make up their own instructions in MAD for the solution of the problem. Again, they are required to submit a flow diagram, and a written report which includes an explanation of the method necessary to use the program successfully on the computer. In the solution of this problem the students are required to write their own Newton-Raphson iteration for the square root solution of the distance. They are also asked to solve for the distance by using either the sine or the cosine of the angle which has previously been computed from the arctangent. This gives them their first introduction to the use of library subroutines, the formal part of MAD programming including the loop concept, and also gives them a chance to compare the relative accuracy by which they can obtain the distance by two completely independent methods. In the third program which is assigned in this course, the students are instructed to write a computer program which will calculate the volume of a borrow pit. In this problem, the students are asked to bring in the "before" and "after" elevation readings over a grid pattern which is basically rectangular. They are asked to work with these two matrices, subtracting the "after" readings from the "before" readings to obtain a new matrix called the "excavation." This subtraction is done by the use of a double loop. Following this, the computations are performed by means of some single and double loops. This helps the student become familiar with methods of handling large sets of subscripted variables. The fourth and last program is one which may vary considerably from semester to semester. It is designed to offer a challenge to the students. It is intentionally changed radically from one semester to the next in order to discourage any ideas of propagating information from one class to another. Examples of the problems which have been given previously are as follows: (1) Intersection of a straight line with a spiral envelope. (2) Approximate adjustment of a triangulation quadrilateral. (3) Progressive adjustment of interlocking traverses by the Bowditch (compass) Method. (4) Crandall method of adjustment of a traverse. At the time of this writing it is strongly believed that the inclusion of computer programming in this course has not detracted from the course work in any way. In fact, it is felt that the blending image of the desk calculator, the digital computer, and surveying computation methods, has resulted in a course which has materially strengthened the students' understanding in each of the three fields mentioned. -C14 -

Use of Computers in Civil Engineering Education Some of the leading exponents of computer philosophy and education feel that geometry is one of the best vehicles for teaching computer ph-il-sophy. The student's ability to follow graphically and pictorially every step in his solutionf makes it possible for him to concentrate on the programming necessary to solve the problem. By the same token, his detailed computer solution of the problem brings him into a very intimate relationship with the problem itself. By the time he has successfully programmed a solution, he is completely familiar with every aspect of the problem. It has been found that many of the students are completely apathetic to computers and the philosophy of computation. These students manage to pass the course by literally hanging on for dear life. Other students in the middle category passively accept the work and perform the assignments satisfactorily but with little imagination. Finally, we have the group who become quite interested in computers and their possibilities and who go beyond the call of duty, polishing up their programs to achieve more than the assignment requires. Certainly these students are in the minority, but it is encouraging to notice that their number seems to be increasing. It is believed that this enthusiasm is growing because the work is passed down to succeeding generations by those who have run into advanced problems which required computer applications. Also, it is quite apparent to the students that many interviewers from industrial firms who come to the campus are now asking about the amount of computer work students have been exposed to. These students, especially those who are doing well in their other subjects (primarily mathematics), are encouraged to sign up for Mathematics 473, Introduction to Digital Computers, a comprehensive, senior-graduate level com coter course. It is not our intention to force the majority of ->r students into computer technology, but it is felt that all students should have some kno;wle-ge of computers and that some few students should have a specialized knowledge of computer programming techniques. In Surveying Computations, it is possible for the students to decide which of these paths they wish to follow. Upon completion of this course, no student will feel that a computer is a mysterious black box which possesses magical properties. On the contrary, he will realize that it is most simple for him to solve elementary problems and that it is not impossible for him to solve even the most difficult problems by using it.

IV. SPECIFIC USES OF COMPUTERS IN VARIOUS FIELDS OF STUDY A. Geodetic Engineering The required courses in surveying are designed to give the student a realistic appreciation for the theory of measurements, application of mathematics to geometric figures, logical analysis, and systematic computation procedures. These courses provide six of the first seven semester hours of contact between the Civil Engineering Department and the students. The first course, Basic Surveying (CE260), does not lend itself to student use of computers. However, the instructor assigns a unique traverse computation problem to each student and compares his answers to those provided by the computer. Other similar uses are being explorec The second required course, Surveying Computations (CE261), has been described in detail in Part III of this report. The terminal required course, Advanced Surveying Measurements (CE362), is taught at Camp Davis, Jackson, Wyoming. This course is primarily one engaging in precise control surveys using theodolites, electronic distance measuring devices and invar tapes. Although no computers are available in the vicinity of the camp which can be used for the course purposes, programs have been written which will make it possible to send data from camp to the University Computing Center and receive the results back in time to use them for the purposes of the course. At the present time these programs are primarily concerned with reducing observation data, such as Tellurometer observations, to yield the actual measured values from the observations. These reductions, when performed by hand, are quite laborious. Even more important to the user, however, such reductions are filled with potential errors, which, because they deal with measurements, cannot be checked. The use of the computer for this purpose virtually eliminates all chance of such blunders in reductions. Geodetic Field Methods (CE562) consists of unique problems which vary with the interests of the students and the program of the department. Usually there are several applications for computer solutions of problems which occur in this course. Adjustment of Geodetic Measurements (CE563) depends very heavily on digital computer applications. Problems are programmed and solved by means of the Gauss-Seidel iterative technique. Matrix inversion and other methods of solving large systems of linear equations are also applied to the computer. Special Problems in Advanced Surveying (CE564). Problems in this course again are quite variable and may or may not require computer applications. Municipal Surveying (CE565) consists in part of several conducted tours to various private and municipal or governmental agencies. Oneprivate computing firm, the Bureau of Surveys of the City of Detroit, and the Michigan State Highway Department were three locations visited where computers were being used. Students observed the computers in operation and were told by the engineers in charge exactly what applications were being made in the specific fields of their needs. -C16 -

Use of Computers in Civil Engineering Education Geodesy and Surveying Research (CE960) is another example of a course which may have varying interests or needs. Certainly many students in the future who take this course will be carrying on extensive computer research in the field of geodetic surveying. B. Hydraulic Engineering During the last three years a majority of the courses offered in hydraulics and hydraulic engineering have been completely re-examined, and now require use of the digital computer as a part of the course work. The required course, Hydraulics (CE421), has in the past offered the student an opportunity to program one or two problems on a voluntary basis. Three or four students, out of about 40, usually wrote programs. Beginning with the spring semester (1962), however, two programs are required of all students: the first the determination of critical depth, normal depth, water surface profiles, and location of a hydraulic jump in a trapezoidal channel; and the second, the determination of the best formula for a v-notch weir, from student data, by the method of least squares. The course in Hydraulic Transients (CE522), has evolved from a hydropower course, primarily because of the great improvement in means of solving problems in unsteady flow. Problems in surge and water hammer have been programmed, and the instructor has set up the differential surge tank problem for class inspection. New methods for handling water hammer have been developed by graduate students as a result of this course and the availability of the digital computer. The course in Advanced Hydraulics (CE524), relies heavily on the computer, and has used it for such problems as relaxation methods for computing flow under a dam and for applications of the momentum theory. One significant problem, the design of a side channel spillway, has been given very thorough study, and a new method has been developed to optimize the design. n Two unknown parameters are needed: a and n in V = axn, for velocity at a cross section x distance along the spillway for minimum total excavation. A "gradient method" is used, to determine how a and n should be changed to move in the direction of minimum excavation. The Hydraulic Design course (CE526), is built entirely on the use of the computer. At present it consists of several programs relating to a project involving construction of a gravity dam, mass diagram study of the inflow hydrograph, flood routing study, reservoir operation study for flood control, power and irrigation supply, and finally an economic study which determines the height of the dam to be constructed. It is believed that the students in the hydraulic option are much better prepared to meet and adjust to the future demands of their profession as a result of their ability to think through the logical processes required in a program. It also gives them a definite competitive edge in securing interesting job assignments with greater growth potential. -C17 -

C. Sanitary Engineering The digital computer is being used for network analysis in Water Supply and Treatment (CE480), as an optional problem which has been worked out successfully by several students. In Sewerage and Sewage Treatment (CE481), the design of a storm sewer system is at present being worked on as an optional assignment. In Industrial Bacteriology (CE587), a problem has been developed requiring the solution on the computer of the MPN (Most Probable Number) equation for determining'bacterial populations. In Advanced Sanitary Engineering Design (CE682), the programming of the sewage network analysis is to be described and discussed but no actual programming is anticipated. However, the students may elect to do so if they desire. D. Soil Mechanics and Foundation Engineering The traditional aim of the courses in soil mechanics and foundation engineering at The University of Michigan has been to train the student to view problems not only in terms of mathematical formulations but also in relation to environmental associations. Particular emphasis has been placed on framing problems in terms of simplified basic relationships which can be conveniently handled by slide rule or desk calculator. One observes that the conscientious student develops a "feel" for the subject in this approach which is more readily adapted to variations in field conditions. It is desired to maintain this emphasis on the practical application of fundamental principles and, at the same time, to take full advantage of the spectacular developments in computational technique with present day electronic computers. In this respect, the computer makes possible more exacting solutions of certain soil mechanics problems and eliminates compromising assumptions and approximations which have heretofore been necessary. Consequently, the soil mechanics teaching and research staff have re-examined a number of the more important problems that require extensive computations and have developed computer programs for the following: 1. Analysis of the plate loading tests to determine the bearing capacity of soil. 2. Determination of pressure distribution from large loaded areas. 3. Static equilibrium of sheet pile retaining walls under a variety of nonhomogeneous soil conditions. Additional problems under investigation include computation of settlement in granular materials from volume-density relationships, a general solution for critical stresses under surface loading and the reproduction of pavement profiles for a variable length of base line reference. The solution of these more complex problems has not advanced to the stage where they can be utilized in undergraduate instruction. They are, however, being developed and used by graduate students working on these specialized problems. The objective in presenting the Buse of computers to students of soil mechanics is to demonstrate the value of these more exacting solutions and familiarize the student with their application to soil mechanics problems. It is felt that the real value of computer techniques

Use of Computers in Civil Engineering Education to soil mechanics lies in the rigorous organization of the problem in a sequence of logical steps leading from fundamental principles and realistic assumptions to a useful practical solution. It is contemplated that the use of computers in soil mechanics will be included in the advanced design course. E. Structural Engineering Computers are used in the structures sequence both at undergraduate and graduate levels. In Elementary Structural Design (CE313), some experiments have been carried on in the use of the computer in simple design problems. This has not been entirely successful, largely because the students have too little experience in the use of the computer in earlier work. There is not sufficient time to include programming in the course lectures. In Theory of Structures (CE512), and Rigid Frame Structures (CE514), the computer has been used to solve problems in truss and beam deflections, using conventional procedures. The use of computers in these courses is still on an experimental basis. In Structural Design (CE513), the computer is used regularly each semester and one computer problem is assigned. In a typical semester, the students are asked to prepare section modulus and shear connector tables for composite beams in accordance with the AISC 1962 Specifications. The problem is relatively simple but involves quite a number of logical decisions that require a more complete knowledge of composite beam behavior than would the comparable hand solution for only a few special cases. There is nothing artificial about the problem. It is the type of problem that one might expect to encounter not infrequently in practice. In Structural Dynamics (CE611), the computer is regularly used for one or two problems each semester. Ordinarily the problems involve numerical solution of the differential equations of motion and constructing response spectrum curves for some given dynamic load. In Structural Plate Analysis (CE613), one problem is assigned each semester in finding the deflections of a plate by finite difference methods. The problem is taken up in class at two levels; first the solution is obtained by hand for a coarse mesh, and later a finer mesh solution is prepared using the computer. Computer Analysis of Structures (CE618), is a course serving the specific purpose of providing the background needed for using the computer efficiently in structural analysis. Matrix methods are of special value in computer analysis of structures, so the course is built on matrix techniques. Most of the course is devoted to formulating the theory of structures in the language of matrices. Assigned problems include small, near-trivial problems where all of the matrix operations are carried out by hand, and larger, non-trivial problems for which the students program the matrix operations for the machine. The smaller problems serve the purpose of familiarizing the student with the details of the procedures and also retaining a physical interpretation for all of the various matrix operations. -C19 -

F. Transportation Engineering The staff use of computers in the area of traffic and transportation is somewhat limited at the present time, but will undoubtedly increase in the future. In CE575, Transportation Terminal Design (highway, rail, water, and air), and CE670, Transportation Planning, the problems, insofar as possible, are placed in a local Ann Arbor context. Both the City of Ann Arbor and Washtenaw County have acquired voluminous data on local transportation, land use, and population. Much of this data is available on punched cards so that it can be sorted in any manner desired. Thus in design problems the computer is used to summarize and compile data on specific geographic areas interconnected by transportation facilities. There are some iterative type solutions involved in layout and design of automobile parking lots and in the delineation of retail trade areas for shopping centers which lend themselves quite nicely to computer application. At the present, a study is being made of the extent to which it is desirable to require the student to use the computer in solving these problems. Transportation planning problems quite often are examined on a systems analysis basis. While the computer has not been used to any great extent in this field, there is a potential for widespread usage. To summarize, the various areas in traffic and transportation courses and the extent to which it is desirable for the student to become involved with the computer are currently being examined. Use at the present time, however, is largely experimental. V. CONCLUSIONS As in the case of every important innovation, it will take several years to evaluate the ultimate benefits resulting from the introduction of computers into Civil Engineering education. At present, opinions of the faculty are varied, as evidenced by the following sampling of statements. "The Civil Engineering Department at The University of Michigan has successfully initiated the use of computers into the curriculum and the faculty is gradually understanding the value of adapting their work to accept and include these new techniques." "The real benefits of using the computer come not from teaching the student how to program, but rather from forcing him to make all of the logical decisions that are necessary to arrive at a suitable algorithm for solving an engineering problem in a more general formulation than he usually encounters. It is often' said that a good way to learn a subject is to teach it. In a sense, the student is being put in a position of teaching the computer how to solve the problem at hand. In order to accomplish this he must necessarily understand the problem more completely than would be necessary for conventional hand solution, since all of the steps involving judgement that are at his disposal in hand methods must be analyzed and programmed when the machine is doing the work." -C20 -

Use of Computers in Civil Engineering Education "While it is agreed that computers are here to stay and are a valuable tool in the civil engineering profession, we must be very careful not to sacrifice the teaching of fundamental engineering principles in order to teach the use of computers." "The method of introducing Civil Engineering students at Michigan to computer techniques by integrating it with the teaching of surveying computations seems to be an excellent approach." "In introducing the use of computers, it is necessary to maintain the emphasis on the practical application of fundamental principles and, at the same time, to take full advantage of the spectacular developments in computational technique with present-day electronic computers. The real value of computer techniques lies in the rigorous organization of the problem in a sequence of logical steps leading from fundamental principles and realistic assumptions to a useful, practical solution." "We need to make fuller use of computers in all fields of Civil Engineering. In this way our graduates will be much better prepared to meet and adjust to the future demands of their profession as a result of their ability to think through the logical processes required to develop algorithms for the computer." VI. EXAMPLE PROBLEMS PREPARED BY PROJECT PARTICIPANTS Each participant prepared solutions for example problems which he felt would be useful in courses at the undergraduate or graduate level. Nine of these problems have been selected for inclusion in this report and are listed below. Number* 76 77 78 79 LIST OF PROBLEMS Title Use of the Computer as a Teaching Machine for Plane Analytic Geometry Optimum Regulation of a Reservoir for Power Production Economical Size of Canal-Feeding HydroPower Plant Most Probable Number Method for Determining Bacterial Populations Iteration of Sum of Principal Stresses at Interior Points of Plane Rings Elastic Buckling Load for Columns of NonUniform Cross Section A MAD Program for Truss Analysis Stiffness Factors for a Flat Slab Column Vibration of Beams on Spring Supports Author H. J. Welch J. W. Howe J. W. Howe L. L. Kempe W. Weaver, Jr. K. H. Chu K. H. Chu R. B. Harris S. S. Kuo Page C23 C35 C41 C47 C56 C67 C82 C97 C113 80 81 82 83 84 * These problems may be considered as a supplement to problems presented in tions of the Project (3, 4, 5, 6) and are numbered accordingly. previous publica -C21 -

All of these problems were solved on the digital computer using either the MAD or FORTRAN language. These languages are described in References 1 and 2. It should be noted that there are several other sources of example computer problems dealing with Civil Engineering subjects. See, for example, Example Problem Numbers 8, 9, 10, 21, 22, 23, and 27 in Reference 3; Example Problem Numbers 48, 49, and 50 in Reference 4; and Example Problem Numbers 63 and 64 in Reference 6. A number of problems prepared by Project participants are abstracted in the Second Annual Report of the Project (4). Some of these (listed above) have been included in this report. Most of the remaining problems are available on microfilm from the Project. -C22 -

Example Problem No. 76 USE OF THE COMPUTER AS A TEACHING MACHINE FOR PLANE ANALYTIC GEOMETRY by Harold J. Welch Department of Civil Engineering The University of Michigan Course: Surveying Computations Credit Hours: 2 Level: Sophomore Statement of Problem In learning surveying computations, students are required to solve a number of problems in plane, analytic geometry using a desk calculator and tables of trigonometric functions. The problems consist of several examples of each of five basic problem types illustrated below. 1. Forward Computation (X2,Y2) Given: X1, Y1, D1, 01 (X1,Y1) i D1 Find: X2' Y2 01 2. Inverse Computation (X2,Y2) Given: X1, Y1, X2, Y2 (Xl,) D1D Find: Di, Q 1 3. Line-Line Intersection Given: X1, Y1, G91 X2, Y2, 02 D D2 Find: X3, Y3, D1, D2 ( 3 1(X' Y2) 4. Line-Circle Intersection Given: X1, Y1, 91' X2, Y2, D2 x, Y3) Find: X3, Y3 82, D1 2 ( Y2) 3-9 3 2D -C23 -

The Computer as a Teaching Machine 5. Circle-Circle Intersection D \ Given: X1, Y1, D1, X2, Y2, D2 (X1 Y ), Find: X3, Y3, 01, G2 The only difference between two problems of the same type is that the numerical data are different. A computer program has been written to check the students' solutions to problems assigned in class and point out any errors. The program has been written to "penalize" students who have incorrect solutions by assigning them additional problems. On the other hand, students with correct solutions are "rewarded" by not being assigned additional problems. Because of this feature of penalties and rewards, the program which grades the student problems has been named the "Teaching Machine Program." Use of the Teaching Machine Program Each student is given an identification number and a set of basic data for each problem type. He then uses his identification number to modify the basic data. For example, if the student's number were 1234, selected distances and point coordinates would be changed by 4253.897 1.234 4255.131 Azimuths (line directions) would be changed according to 290~ 47' 18" 12 34 290~ 59' 52" In this manner, all the students work essentially the same problem, but with slightly different data and results. When the student has solved a problem to his satisfaction, he punches both his modified data and his computed results on IBM cards. On the first card, the station number and the coordinates of the first point are punched. The length and azimuth of the line between the first two points are punched on the second card. This sequence is repeated until all of the problem data havebeen given. Hence, for the Forward and Inverse Computations, each problem requires three cards. For each Line-Line Intersection problem, five cards are required. Examples of input cards are shown on the next page. -C24 -

Example Problem No. 76 A}, n Bra Col 1 3 M id E2 655 3730.751 3636.576 Station Coordinates Number Col 26 38 42 45 1415.702 332 36 13 Length Azimuth On the first run of each type of problem, the student must work 5 problems (3 for the Circle-Circle intersection). In front of the data cards, he places a card bearing his name, class number, computing center number, and the type number of the problems which follow. For example, a name card might be Col 1 21 32 39 /J Q STUDENT 1234 X106Z 4 Student Identifica- Computing Problem Name tion Number Center Type Number Number The data for each type of problem must be submitted separately, with an individual name card for the set of data for each type. At execution, the program first reads the current date to be used in output headings. Following this is the class roll. Each card in this deck contains a student's name and class number, the number of problems due and the number of the run to be made for each of the five types of problems. Next, the five sets of basic data are read and then the first student's data is considered. The class roll, the sets of basic data, and each set of student data is followed by a card having an equal sign punched in column 1. Whenever this equal sign is read, control is transferred to the portion of the program which handles the type of data expected to follow. In processing student data, the information on the student's name card is read and printed in a heading. Then the student's class number is checked to see if it is a legal number. The corresponding run number is checked to see if the student is attempting an extra run. Finally, the student's name is checked to see if it corresponds with the given class -C25 -

The Computer as a Teaching Machine number. If any of these tests fail, an appropriate comment is printed and the program searches for the name card of the next student. If the information is all legal, control is transferred to the section of the program which handles the type of problem given on the name card. Each of these five sections modifies the basic data in the same manner as the student and computes the correct solution to the problem. These results are then used by the internal function CHECK. CHECK reads the student's data and compares it with the computed data, with a tolerance of.003 feet on distances and coordinates and 3 seconds on azimuths. The student's data are then printed, with asterisks flagging any incorrect values. Whenever an error occurs, a boolean flag is set, and at the end of the problem, the value of the flag is returned to the main program. If the flag has been set, the main program increments an error counter. Then the data for the next problem (of the same type) are examined, the correct solutions are computed, and CHECK is again called. This process continues until the required number of problems for this sequence have been considered or until an equal-sign card is read. Extra problems are ignored and simply scanned for an equal-sign card. At the end of the set of data, if the error counter is not equal to zero, or if too few problems have been submitted, the number due is increased. On the next run, the student must resubmit all the correct problems, the problems which were incorrect or missing (hopefully corrected by this time), and an extra problem for each incorrect or missing problem. This new due number is printed following the data listing. If all the problems are correct, the due number is set to zero and a comment is printed. The run number is now incremented and the program goes on to the next student's data. After all the problem sets have been processed, the updated class roll is punched for use in the next run and printed for the information of the instructor. Flow Diagram A general flow diagram of the Teaching Machine Program is shown below. C 1 - Forward Computation C 2 - Inverse Computation C 3 - Line-Line Intersection C 4 - Line-Circle Intersection C ) - Circle-Circle Intersection -C26 -

Example Problem No. 76 MAD Program The Teaching Machine Program (which is written in the MAD language) is listed below. H J WELCH 217 W ENG X106W 002 050 050 9C8 H J WELCH 217 W ENG X106W 002 050 050 9C8 $COMPILE MADEXECUTEPRINT OBJECTPUNCH OBJECTDUMP 629C8001 EXECUTE SETEOF.(LAST) EXECUTE SETERR.(SRCH) READ FORMAT RDATEDATEoATE1DATE2 THROUGH ROLL. FOR N=11tN*G1000 READ)- FORMAT RD1,EQUALSNAME(N),NAME1(-N) NAME2(N —CN-(-N-T-DUEIITN 1),DUE2(N) DUE3(N),DUE4(N),DUE5(N),RUN1(N)RUN2(N),RUN3(N) RUN 24(N) RUN5(N ) ROLL WHENEVER EQUALS.E.$=$,TRANSFER TO OUT PRINT CO MMENT$4T00 MANY STUDENTS$ EXECUTE SYSTEM. OUT THROUGH LFOR- Ji,,J- G25 READ FORMAT RDEQUALS,FROM(J),TO(J),FRX(J),FRY(J),DTO(J),D(J) 1.M(J) S(J) L WHENEVER EQUALS*Ee$=$W1RANSFER TO NEXT3 RTINTCOfMNENT -'$4TO0 MUCH DATA FOR PROBLEM 1$ EXECUTE SYSTEM* NEXT3 THROUGH L2,FOR J=l,1,J1G425 READ FORMAT RD2,EQUALSFROM2(J),FRX2(J)FR(J) (J)T02(J)*TOX2(J 1ltTOY2(J) L2 WHENEVER EQUALS*E.$=$,TRANSFER TO NEXT4 PR I NT COMMENT $4TOO MUCH DATA O'R -PROBLEM —2$' EXECUTE SYSTEM. NEXT4 HTHROI -L3 FOR ~J1-i J G25 READ FORMATRD3,EQUALSFROM3(J),FRX3(J),FRY3(J),FD3(J),FM3(J), 1FS3(J) INT3(J) TD3(J),TM3 (J),TS3(J) T03T(J,TOX3(J) TOY3(J) L3 WHENEVER EQUALS.Ea$=$,TRANSFER TO NEXT1 'P-INT —CO-MMENT$4TO'-1o-MUH-' DATA 'F'R- PR-OBLEM 3$ —' EXECUTE SYSTEM. NEXT 1 THiROU-GFI- 4 ' FOR-J=-1I -i"-J.G,25 READ FORMAT RD4,EQUALSFROM4(J),FRX4(J),FRY4(J),FD4(J),FM4(J) -C27 -

The Computer as a Teaching Machine MAD Program, Continued IFS4(J),INT4(J),DTO4(J),TO4(J),TOX4(J),TOY4(J) L4 WHENEVER EQUALS.E.$=$,TRANSFER TO NEXT2 PRINT COMMENT$4TOO MUCH DATA FOR PROBLEM 4$ EXECUTE SYSTEM. NEXT2 THROUGH L5, FOR J=1~1,J.G.25 READ FORMAT RDSEQUALStFROM5(J),FRX5(J),FRY5(J),DFM5(J),INT5( 1J),DTO5(J),T05(J),TOX5(J),TOY5(J) L5 WHENEVER EQUALSE.$=$,TRANSFER TO NEWNAM PRINT COMMENT$4TOO MUCH DATA FOR PROBLEM 5$ EXECUTE SYSTEM. NEWNAM ERR=O TER=06B READ fORMAT RDNAMEEQUALSSNAME..,SNAME(2) SCN,PROB WHENEVER EQUALS.E.$=$.TRANSFER TO NEWNAM PRINT FORMAT TITLEBSNAME*.-SNAME(2),SCNPROB,DATE,DATE1 1,DATE2 THROUGH L1i FOR K=1,,K.E.NN L1 WHENEVER SCN.E.CN(K),TRANSFER TO TSTRUN(PROB) PRINT COMMENT $OINCORRECT CLASS NUMBER$ TRANSFER TO SRCH TSTRUN(1) WHENEVER RUN1(K).L.4,TRANSFER TO TSTNAM PRINT COMMENT $ORUN LIMIT EXCEEDED ON PROBLEM 1$ TRANSFER TO SRCH TSTRUN(2) WHENEVER RUN2(K).L.4oTRANSFER TO TSTNAM PRINT COMMENT $ORUN LIMIT EXCEEDED ON PROBLEM 2$ TRANSFER TO SRCH TSTRUN(3) WHENEVER RUN3(K).L*4#TRANSFER TO TSTNAM PRINT COMMLNT $ORUN LIMIT EXCEEDED ON PROBLEM 3$ TRANSFER TO SRCH TSTRUN(4) WHENEVER RUN4(K).L.4,TRANSFER TO TSTNAM PRINT COMMENT$ORUN LIMIT EXCEEDED ON PROBLEM 4$ TRANSFER TO SRCH TSTRUN(5) WHENEVER RUNS(K).L.4,TRANSFER TO TSTNAM PRINT COMMENT$ORUN LIMIT EXCEEDED ON PROBLEM 5$ TRANSFER TO SRCH TSTNAM WHENEVER NAME(K)*NE.SNAME.OR.NAME1(K).NE.SNAME(1).OR.NAME2(K 1)oNE.SNAME(2) PRINT COMMENT $OINCORRECT NAME$ SRCH READ FORMAT SEARCHEQUALS WHENEVER EQUALS.E.$=$*TRANSFER TO NEWNAM TRANSFER TO SRCH END OF CONDITIONAL TRANSFER TO C(PROB) C(1) PRINT FORMAT TITLECPROBDUE1(K),RUN1(K) THROUGH BIGLP1,FOR I=l,l,I.G.DUE1(K) STA=FROM(I) STA(1)=TO(I) X=FRX(I) Y=FRY(I) DIST=DTO(I)+SCN/1000. D1=D(I) M1=M(I)+SCN/100 S1=S(I)+SCN-100*(SCN/100) SEC=3600*Dl+60*M1+S1 RA=,48481368E-5*SEC X(1)=FRX(I)-SIN.(RA)*DIST Y(1)=FRY(I)-COS(RA)*DIST TER=CHECK.(THRU1) WHENEVER TER ERR=ERR+1 TER=OB BIGLP1 END OF CONDITIONAL WHENEVER ERR*E.O DUE1(K)=O PRINT FORMAT OKPROB TRANSFER TO RK1 END OF CONDITIONAL DUE1(K)=DUE1(K)+ERR TEST1 WHENEVER DUEl(K).G.15,DUE1(K)=15 PRINT FORMAT RESMITDUE1(K) RK1 RUN1(K)rRUN1(K)+l TRANSFER TO SRCH THRU1 DUE1(K)=2*DUE1(K)+l+ERR-I TRANSFER TO TEST1 C(2) PRINT FORMAT TITLECPROBDUE2(K),RUN2(K) THROUGH BIGLP2,FOR I=1,,1I. G.DUE2(K) STA=FROM2(I) -C28 -

Example Problem No. 76 MAD Program, Continued X=FRX2(I)+SCN/1000. Y=FRY2(I) DX=TOX2(I)-SCN/1000.-X DY=TOY2(I)-FRY2(I) DIST=SQRT.(DX*DX+DY*DY) SEC=ATN1.(-D —DY)/0.48481368E-5+0.5 STA(1)=T02(I) X(1)=X+DX Y(1)=TOY2(I) TER=CHECK (THRU2) WHENEVER TER ERR=ERR+1 TER=OB BIGLP2 END OF CONDITIONAL WHENEVER ERR.E.0 DUE2(K)=0 PRINT FORMAT OKPROB TRANSFER TO RK2 END OF CONDITIONAL DUE2(K)=DUE2(K)+ERR TEST2 WHENEVER DUE2(K)*G.15*DUE2(K)=15 PRINT FORMAT RESMITDUE2(K) RK2 RUN2(K)=RUN2(K)+1 TRANSFER TO SRCH THRU2 DUE2(K)=DUE2(K)*2+1+ERR-I TRANSFER TO TEST2 C(3) PRINT FORMAT TITLECPROBDUE3(K),RUN3(K) THROUGH BIGLP3# FOR I=1,lI.G.*DUE3(K) STA=FROM3(I) STA(1)=INT3(I) STA(2)=T03(I) X=FRX3(I) Y=FRY3(I) X(2)=TOX3-(I) Y(2)=TOY3(I) D1=FD3( I) D2=TD3(I) M1=FM3(I)+SCN/100 M2=TM3(I)+SCN/100 S1=FS3(I)+SCN-100*(SCN/100) S2=TS3(I)+SCN-100*(SCN/100) SEC=3600*D1+60*Ml+S1 SEC(1)=3600*D2+60*M2+S2 RA=.48481368E-5*SEC RB=.48481368E-5*SEC(1) SA=-SIN.(RA) SB=-SIN.(RB) CA=-COS.(RA) CB=-COS.(RB) DX=X(2)-X DY=Y(2)-Y DA=(DX*CB-DY*SB)/SIN.(RA-RB) X(1)=X+DA*SA Y(1)=Y+DA*CA DB=(DY-DA*CA)/CB DIST=DA DIST(1)=DB TER=CHECK.(THRU3) WHENEVER TER ERR=ERR+1 TER=OB BIGLP3 END OF CONDITIONAL WHENEVER ERR.E.O DUE3(K)=O PRINT FORMAT OK,PROB TRANSFER TO RK3 END OF CONDITIONAL DUE3(K)=DUE3(K)+ERR -C29 -

The Computer as a Teaching Machine MAD Program, Continued TEST3 WHENEVER DUE3(K)J.G15#DUE3(K)=15 PRINT FORMAT RESMITDUE3(K) RK3 RUN3(K)=RUN3(K)+1 TRANSFER TO SRCH THRU3 DUE3(K)=2*DUE3(K)+1+ERR-I TRANSFER TO TEST3 C(4) PRINT FORMAT TITLECPROBDUE4(K),RUN4(K) THROUGH BIGLP4~ FOR I1,1,I.G.DUE4(K) STA=FROM4(I) STA(1)=INT4(I) STA(2)=T04(I) X=FRX4(I) Y=FRY4(I) X(2)=TOX4(I) Y(2)=TOY4(I) D1=FD4(I) M1=FM4(I)+SCN/100 S1=FS4(I)+SCN-100*(SCN/100) SEC=3600*D1+60*Ml+S1 WHENEVER DT04(I).L.O. DIST(1)=DT04(I)-SCN/1Q000 OTHERWISE DIST(1)=DT04(I)+SCN/1000. END OF CONDITIONAL DX=X(2)-X DY=Y(2)-Y RA=.48481368E-5*SEC SA=-SIN.(RA) CA=-COS.(RA) Q=DX*SA+DY*CA H=DX*CA-DY*SA F=SQRT.(DIST(1)*DIST(1)-H*H) WHENEVER (DX*DX+DY*DY-DIST(1)*DIST(1) ).L*.0OR.DIST(1).L.*O DIST=Q+F OTHERWISE DIST=Q-F END OF CONDITIONAL X(1)=X+SA*DIST Y(1)=Y+CA*DIST RB=ATN1l(X(1)-X(2).Y(1)-Y(2)) SEC(1)=RB/*8448i368L-5+.5 TER=CHECK*(THRU4) WHENEVER TER ERR=ERR+1 TER=OB BIGLP4 END OF CONDITIONAL WHENEVER ERR.E.O DUE4(K)=O PRINT FORMAT OK,PROB TRANSFER TO RK4 END OF CONDITIONAL DUE4(K)=DUE4(K)+ERR TEST4 WHENEVER DUE4(K).G.15#DUE4(K)=15 PRINT FORMAT RESMITPDUE4(K) RK4 RUN4(K)=RUN4(K)+1 TRANSFER TO SRCH THRU4 DUE4(K)=2*DUE4(K)+1+ERR-I TRANSFER TO TEST4 C(5) PRINT FORMAT TITLECPROBDUE5(K),RUN5(K) THROUGH BIGLP5, FOR I=1,1,I*G.DUE5(K) STA=FROM5(I) STA(1)=INT5(I) STA(2)=TO5(I) X=FRX5(I) Y=FRY5(I) X(2)=TOX5(I) Y(2)=TOY5(I) DIST=DFM5(I)+SCN/1000. DIST(1)=DTO5( )+SCN/1000. DX=X(2)-X DY=Y(2)-Y DC=SQRT.(DX*DX+DY*DY) F=(DC*DC+DIST*DIST-DIST(1)*DIST(1))/(2.*DC) H=SQRT.(DIST*DIST-F*F) S3=DX/DC C3=DY/DC X(1)=X+F*S3-H*C3 Y(1 ) Y+F*C3+H*53 RA=ATN1.(X-X(1) Y-Y(1)) -C30 -

Example Problem No. 76 MAD Program, Continued RB=ATN1.(X(1)-X(2),Y(1)-Y(2)) SA=-SIN.(RA) SB=-SIN.(RB) CA=-COS.(RA) CB=-COS.(RB) SEC=RA/ 48481368E-5+.5 SEC(1)=RB/.48481368E-5+.5 TER=CHECK.(THRU5) WHENEVER TER ERR=ERR+1 TEK=OB BIGLP5 END OF CONDITIONAL wHENEVER ERR.E.O DUE5(K}=O PRINT FORMAT OK,PROB TRANSFER TO RK5 END OF CONDITIONAL DUE5(K)=DUE5(K)+ERR TEST5 WHENEVER DUE5(K).G.15,DUE5(K)=15 PRINT FORMAT RESMITtDUES(K) RK5 RUN5(K)=RUN5(K)+1 TRANSFER TO SRCH THRU5 DUE5(K)=2*DUE5(K)+l+ERR-I TRANSFER TO TEST5 LAST PRINT FORMAT PDATEDATE,DATE1,DATE2 THROUGH DONE. FOR I=1,1,I.E.N PRINT FORMAT PR1,NAME(I),NAME1(I),NAME2(I),CN(I),DUE1(I).DUE2 1(I),DUE3(I),DUE4(I),DUE5(I),RUN1(I),RUN2(I),RUN3(I).RUN4(I). 2RUN5(I) DONE PUNCH FORMAT PCHNAME(I)tNAME1(I),NAME2(I CNCN(I)DUE1(I),DUE2 1(I),DUE3(I),DUE4(I)~DUE5(I)tRUN1(I)tRUN2(I)gRUN3(I)~RUN4(I), 2RUN5(I) EXECUTE SYSTEM. INTERNAL FUNCTION (HERE) STATEMENT LABEL HERE BOOLEAN FLAG ENTRY TO CHECK. CK=O FLAG=OB NWCD READ FORMAT NEWCADEQUALS,ST(1),ST(3),ST(5) WHENEVER EQUALS.E.$=$,TRANSFER TO HERE WHENEVER ST(1).E.STA(CK) ST(2)=$ $ OTHERWISE ST(2)=$***$ FLAG=lB END OF CONDITIONAL WHENEVER.ABS.(FST(3)-X(CK)).G..003 ST(4)=$***$ FLAG=1B OTHERWISE ST(4)=$ $ END OF CONDITIONAL WHENEVER.ABS.(FST(5)-Y(CK)).G..003 ST(6)=$***$ FLAG=1B OTHERWISE ST(6)=$ $ END OF CONDITIONAL PRINT FORMAT OUT2,ST(1)...ST(6) WHENEVER CK.El.1AND.PROBI.L.3TRANSFER TO FINI WHENEVER CK.LE.1 READ FORMAT NEXCAD,EQUALS,ST(7),ST(9)...ST(11) WHENEVER EQUALStE.$=$,TRANSFER TO HERE WHENEVER.ABS.(FST(7)-DIST(CK)).G..003 ST(8)=$***$ FLAG=1B OTHERWISE ST(8)=$ $ END OF CONDITIONAL BIGSEC=3600*ST(9)+60*ST(10)+ST(11) -C31 -

The Computer as a Teaching Machine MAD Program, Continued WHENEVER.ABS*(SEC(CK)-BIGSEC) G.3 ST(12)=$***$ FLAG=lB OTHERWISF ST(12)=$ $ END OF CONDITIONAL PRINT FORMAT OUT3,ST(7).**ST(12) CK=CK+1 TRANSFER TO NWCD END OF CONDITIONAL FINI FUNCTION RETURN FLAG END OF FUNCTION EQUIVALENCE (FSTST) BOOLEAN TER BOOLEAN CHECK. INTEGER STDATE.DATE1tNEQUALSNAMENAME1,NAME2,CNDUE3~DUE4, 1DUE5,RUN3,RUN4,RUN54JFROM3,FD3,FM3,FS3,INT3,TD4,TM3*TS3,T03, 2FROM4tFD4,FM4tFS4,INT4,T04,FROM5,INT5 T05 ERRSNAMESCNPROB, 3KSTAIDl1D2oMl1M21SlS2,SECoCKSTtBIGSEC INTEGER DUElDUE2,RUN1,RUN2oFROMTODMS$FROM2,T02TDATE2 DIMENSION NAME(100),NAME1(100),NAME2(100),CN(100),DUE3(1 ),D 1UE4(100)tDUE5(100)~RUN3(100),RUN4(100),RUN5(100)tFROM3(25),FR 2X3(25) FRY3(25),FD3(25),FM3(25),FS3(25),INT3(25),TD3(25) 3TM3(25),TS3(25)*TO3(25) TOX3(25) TOY3(25)DFROM4(25),FRX4(25)~ 4FRY4(25),FD4(25),FS4(25)tFM4(25),INT4(25),DT04(25),T04(25), 5TOX4(25),TOY4(25),FROM5(25),FRX5(25),FRY5(25),DFM5(25),INT5(2 65),DT05(25), T05(25)tTOX5(25)~TOY5(25)S5NAME(2)~STA(2),X(2), 7Y(2),DIST(1),SEC(1)ST(12),FST(12) 8,DUE1(100),DUE2(100) RUN1(100)RUN2(100) DIMENSION FROM(25),TO(25),FRX(25),FRY(25),DTO(25),D(25),M(25) 1S(25),FROM2(25),FRX2(25),FRY2(25)~T02(25)~TOX2(25)>TOY2(25) VECTOR VALUES RDATE=$3C6*$ VECTOR VALUES TITLEC=$13H-INTERSECTIONI3,I6,26H PROBLEMS DUE 1 FOR RUN NO-I3//*$ VECTOR VALUES PRI=$1H S,20t3C6,S6 tI4,S4,10I6 *$ VECTOR VALUES RDNAME=$C1,3C6,I,5S14l1*$ VECTOR VALUES TITLEB=$1H1S5,3C6,S5*I4,S5,7HPROBLEMI3,S52C6*$ VECTOR VALUES SEARCH=$C1*$ VECTOR VALUES OK=$1HOS5,45HYOU HAVE DONE ALL THE PROBLEMS FOR 1 ASSIGNMENT I4,13H CORRECTLY./45HOSEE YOUR INSTRUCTOR FOR 2FURTHER ASSIGNMENTS.*$ VECTOR VALUES RESMIT=$1HO5555HYOU HAVE AT LEAST ONE ERROR (I 1NDICATED BY ASTERISKS). /34HOCORRECT THE DATA AND RESUBMIT T 2HEI4,20H PROBLEMS NOW DUE.*$ VECTOR VALUES PCH=$S1,3C69151I2,9I3*$ VECTOR VALUES RD=$(iA2I4~3F10.3,S55313*$ VECTOR VALUES RD1=$C1*3C6t151I2.9I3*$ VECTOR VALUES RD2=$C1,Z(15v2F10.3)*$ VECTOR VALUES RD3=$Cl13.2F,.3I5,213,2I5t213,15,2F9.3*$ VECTOR VALUES RD4=$Ct113t2F9.3,I5,213,A,FlO.3>I6t2F9.3*$ VECTOR VALUES RD5=$Cl1I3t3F9.316,FlO.316,2F9.3*$ VECTOR VALUES NEXCAD=$C1,S25tF9.3#16,2I3*$ VECTOR VALUES NEWCAD=$C1,I4,2F10.3*$ VECTOR VALUES OUT2=$1H S55I4,C3,F10.3C3,F10.3sC3*$ VECTOR VALUES OUT3=$1H S37,F10.3,C3,I5,213IC3 *$ VECTOR VALUES PDATE=$1HltS19,65H PLANE GEOMETRY TEACHING 1MACHINE MODEL 3-62 PROGRAM 9C8 3C6//520,12HSTUDENT NAME, 258912HCLASS NUMBERS3,58HDUE1 DUE2 DUE3 DUE4 DUE5 RUN1 3RUN2 RUN3 RUN4 RUN5//*$ END OF PROGRAM S DATA -C32 -

Example Problem No. 76 Data The following is an abbreviated set of data which might be submitted to the Teaching Machine Program. The class roll has been shortened to four names. Basic data for only one problem type are shown and only two sets of student's results are shown. $DATA MAY 17, 1962 D J PIKE J A FOX H W ANDRES S H HAMMOND 4116 3 1949 3 1144 3 27U3 3 3 3 3 3 1 1 1 1 1 3 3 3 3 1 1 1 1 1 3 3 33 1 1 1 1 1 3 3 33 1 1 1 1 1 Class Roll 147 148 161 233 168 167 234 99 85 18 19 21 45 36 46 89 254 162 234 227 228 161 26 29 15 20 16 41 44 1 5000.000 4569.374 3228.603 3325.440 3911.262 3803.7~6 3234.698 4896.323 4511.209 4303.953 3784.575 5110.979 4482.951 3539.884 4542.456 5000.000 4430.280 4458.202 4388.835 4533.392 4505.241 4344. 77 4957.408 4552.191 4984.691 4539.076 4 32.670 4623.0752 4589. 095 4538o.92 253.897 15'7 * 16 93.000 101. 16O 108. 887 103.256 114.2 8 1875.690 2165.140 908.394 7804284 Z89.450 099.568 543.067 65U. 158 341 1- 59 345 50 01 266 56 35 63 44 43 00 28 34 342 49 0U 176 56 34 157 39 59 171 17 46 324 45 23 272 27 17 157 27 17 145 12 50 272 54 44 224 43 41 Data for First Problem Type / H W ANDRES 147 5000.000 1144 5000.000 TOO3A 1 255. 41 341 29 43 158.160 346 01 45 89 5080.946 4758.145 148 4569.374 4430.280 254 4607.958 161 3228.603 4276.799 4458.202 } Data for Other Four Problem Types Would be Inserted Here First Student's Solution Seoond Student's Solution 94.144 267 08 19 162 3322.630 4462.996 D J PIKE 147 5000.000 89 5079.782 148 4569.374 254 4606.930 161 3228.603 162 3325.638 4118 5000.000 4574.630 4430.280 4273.564 a458. u2 4462.*16 T0c3P 1 258. 15 341 59 17 161.134 346 31 1i 97.113 267 37 53 -C33 -

The Computer as a Teaching Machine Computer Output The following is the computer output which corresponds to the data shown in the previous section. H W ANDkES 1144 PRO1bLEM 1 MAY 17, 1962 INTERSECTION 1 3 PROirLEMS DUE FOR RUN NO. 1 147 5000.0O 500o0.uO0 255.041 341 29 43 89 50g0.946 47)8.145 148 4569.374 4430.280 158.160 346 1 45 Th 254 4607.5U 8 42 6. 99 tc 161 3228.603 4458.202 94.144 267 8 19 162 3322.630 4462.996*** YOU HAVE AT LEAST ONE ERROR (INOICATED BY ASTERISKS). CORRECT THE DATA AND RESUBMIT THE 4 PROBLEMS NOW DUE. D J PIKE 4118 PROBLEM 1 _ __ MAY 17, 1962 INTERSECTION 1 3 PROBLEMS DUE FOR RUJN NO. 1 147 5000.000 0.0o —o - 258.015 341 59 17 - 89 5079.782 4574.630*** 148 4569.374 4430.280 161.134 346 31 19 254 46C6.930 4273.584 161 3228.603 4458.202 97.118 267 37 53 162 3325.638 4462.216. YOU HAVE AF LEAST ONE ERROR (INDICATED BY ASTERISKS). CORRECT [HE DATA AND RLSUHtMIT THE 4 PROBLEMS NOW DUE. This summary is printed for use by the instructor: INVERSE CUMPUTATION TEACHING MACHINE MODEL 3-62 PROGRAM 9C8 MAY 17, 196 STUDEN NAME CLASS NUMBER DUE1 DUE2 DUE3 DUE4 DUE5 RUN1 RUiN2 D J PIKE 4118 4 3 3 3 3 2 1 J A F-'X 1949.3 3 3 3 3 1 1 H W ANDRES 1144 4 3 3 3 3 2 l S H HAMMOND 2703 3 3 3 3 3 1 1 Number due on next run Run numb for each problem --- ype —. f-or each lis is returned ) first student. iis is returned D second student. 2 RUN3 RUN4 RUN 5 1 1 1 __ 1 1 1 1 1 _ 1 1 1 1 er of next run problem type. - This is the updated class roll which is used as input on the next run: UJ PIKe 411 4 3 3 3 3 2 1 1 1 1 J A FOX 1949 3 3 3 3 3 1 1 1 1 1 H W ANDRES 1144 4 3 3 3 3 2 1 1 1 1 S H HAMMOND 270 3 3 3 3 3 1 1 1 1 1 Critique A complete discussion of the effectiveness of the use of the Teaching Machine Program is given in the section "Introduction of Civil Engineering Students to the Use of the Computer," in the first part of this report. -C34 -

Example Problem No. 77 OPTIMUM REGULATION OF A RESERVOIR FOR POWER PRODUCTION by J. W. Howe Department of Mechanics and Hydraulics State University of Iowa Course: Water Power Engineering II Credit Hours: 1-3 Level: Graduate Statement of the Problem It is desired to determine the best way of operating the turbines at a reservoir in order to produce the most power in the period studied. The following physical relations are assumed: Reservoir volume = (Y)3'5/4 (acre feet) Spillway discharge = 120,000 h3/2 (acre feet/month) Turbine discharge = 12,000 H1/2 (acre feet/month) Tailwater discharge = 21,000 Z3/2 (acre feet/month) Reservoir spillway level at Y = 63 feet Reservoir inflows: 12, 16, 90, 210, 47, 44, 27, 19, 45, 64, 31, 30, 10, 50, 199, 137, 59, 231, 307, 239, 75, 50, 34, 34 (thousand acre feet/month) Determine the average power output per year if the turbines are shut down for various percentages of time when the upstream water level is below spillway level. Solution A cross-sectional diagram of the reservoir is shown below to illustrate the nomenclature. i |zH i Q TW QTURB Z -C35 -

Optimum Regulation of a Reservoir for Power Production Initial levels Yl and Z1 at beginning of month. Final levels YAV and Z2 at end of month. VOLYAV = YAV3 5 VOLO = 633'5 VOLYI = YI3'5 QSTOR = (VOLYAV-VOLY1)/4 YAV = (YMAX + YMIN)/2 = level at end of month AVHD = (YI + YAV)/2 STORAGE EQUATION: INFLOW - QSTOR = OUTFLOW = QSPLWY + QTURB = QTW INFLOW = Q(I) = QTW + QSTOR = QTOT (1) Initial values of YMIN and YMAX of 20 and 68 feet are assumed and YAV determined. Using YAV in determining terms of Eq. (1), if 1000 Q(I) > QTOT, Y should be increased. This is accomplished by letting YAV replace YMIN. In the opposite situation, YAV becomes YMAX. In 15 iterations YAV is within a few thousandths of the value it must have to satisfy Eq. (1). It is assumed that the monthly discharge through the turbines or over the spillway can be secured by using the average of the initial and final water levels. In the case of a month having some spillway discharge in the month but not during the entire month, due to rising or falling water levels, the spillway discharge for a full month is multiplied by: (VOLYAV - VOLO)/(VOLYAV - VOLYI) for a rising stage and by (VOLYI - VOLO)/(VOLY1 - VOLYAV) for a falling stage QZ+Z2 0.746X0.80 43560 Zi+Z2 Power production/month = QTURB (AVHD - 2 )x 7468 80 X 456 =.82 QTURB (AVHD Z- List of Symbols Q(I) Monthly flow, thousands of acre feet/month Yl Upstream depth at beginning of month, feet Y2 Upstream depth at end of month, feet Z1 Downstream depth at beginning of month, feet Z2 Downstream depth at end of month, feet YMIN Minimum upstream depth at end of month, feet YMAX Maximum upstream depth at end of month, feet YAV Average upstream depth at end of month, feet AVHD Average upstream depth during month, feet -C36 -

Example Problem No. 77 Lst of Symbols, Continued PURB Monthly flow through turbine, acre feet/month 3TOR Volume change in storage in one month, acre feet 3PLWY Monthly spillway discharge, acre feet/month DW Monthly outflow, acre feet/month (I) Monthly power output, KWH/month \POWN Annual power output, KWH/year JLT Fraction of time turbines are operated when upstream water level is below spillway level, dimensionless Low Diagram /^ ^ ( -- -- - / -. EXECUTE ZE1o o ODIN( Q... (.Q(24) ^ UT tMV T< z 2 SUKC)-S<. - L I, \,p L;,>u(I zz, s TOTPOW=-O U(C,)...~Cu(4) \^^ = 4 ______3__:T:, CUNT. ( = 1 "' 1. rO: V6Y.Y- 3$T. VOLYAV-VOLY1( IK CO\NTrI,,CGUNT),, 2. yNT 'AOOOVV2 D W R3. * -/ YMAW ^;D Y0 + ^ YAYV ________~~v,i:: Y1 +V 2. -C37 -

Optimum Regulation of a Reservoir for Power Production Flow Diagram, Continued KE \ 1 K>) Fv t AWM2K-TR(V AWMIN( o P(l-O. 20-QTaR(W-: YZ =2 SUM(H)- SU~~(trl tUM().-CVU<+ MULT (7)-^ <?f ) ---^~NWANPTO T WN VOW - sum(.IH ^ MAD Program and Data J. We HOWE S164A 005 010 000 STOREG J. W. HOWE S164A 005 010 000 STOREG SCOMPILE MAD EXECUTEDUMP R ROPTIMUM REGULATION OF A POWER DEVELOPMENT AT A RESERVOIR RASSUMED RELATIONSHIPS, VOLUMNE OF RESERVOIR~Y.Pe3S5/4. RQTW=21000*ZP* 1 5 OQTURB=12000*SORT (H) RQSPLWYz120000*(Y-63) P* 1 5 POWER*0*H*e820 KWH/MO RQUANTITIES IN ACRE FEET/MO. SPILLWAY CREST AT Y~63FT. R DIMENSION Q(24). SUM(4),CUM(4),P(24) INTEGER Q#COUNT#ItJ4K#SUMMCUM PRINT COMMENT S15 READ FORMAT DATA Q(1) * Q(24) VECTOR VALUES DATA-$(1216)*S PRINT FORMAT ECHO*Q(1)t l Q(24) VECTOR VALUES ECHOIS1H 818*S THROUGH OMEGA. FOR MULT'l*-lt1 MULT.LT.e.2 Y1l63. Z1l2.64 Z2'3, TOTPOW O EXECUTE ZERO.(SUM( 1 ) * SUM(4 ).CUM( 1...CUM(4) THROUGH ALPHAJ FOR I"11pIe. G24 YMIN"20 YMAXt68 THROUGH BETA*FOR COUNT"1, 1 COUNToG 15 YAV'tYMIN+YMAX)/2 AVHD (Y1+YAV)/2. QTURB=12000**SQRT (AVHD-( Z1+Z2 ) /2 ) VOLYAV=YAV~P~3*5 VOLO63. P*3*5 VOLYV1Yl1Pe3e5 QSTOR (VOLYAV-VOLY1) /4 WHENEVER YAV GE.634 WHENEVER AVHD.Ga63 QSPLWYu120000*(AVHD-63..P.1.5 OTHERWISE QSPLWY 120000 * YAV-63 )/2 ) P 1. 5* ( VOLYAV-VOLO )tVOLYAV1VOLY1) END OF CONDITIONAL QTW-QSPLWY+QTURB OTHERWISE WHENEVER AVHD.G~63o QSPLWYr120000* ((Y1-63 ) /2 ) P 15*(VOLY1-VOLO) / (VOL1-VOLYAV 1) QTW-QSPLWY+QTURB OTHERWISE OSPLWY-O QTUR BQTURB*MULT QTW-QTURB END OF CONDITIONAL END OF CONDITIONAL QTOT-QTW+QSTOR WHENEVER 1000*Q(I).G.QTOT YMIN-YAV OTHERWISE -C38 -

Example Problem No. 77 MAD Program and Data, Continued YMAX=YAV END OF CONDITIONAL BETA Z2 (QTW/21000. ) P.667 PRINT FORMAT CHECK YAVtQTURBQSPLWYtQTOT VECTOR VALUES CHECK$1HOFl0*2*3F100*$ VECTOR VALUES HWMAX(1)=65.64.564.,63. VECTOR VALUES HWMIN(1)=56*,58*,60,61* THROUGH ADD, FOR KltlvK*G.4 WHENEVER YAVEGE*HWMAX(K) SUM( K)SUM ( K )+ OR WHENEVER YAV.LE*HWMIN(K) CUM(K)=CUM(K)+1 ADD END OF CONDITIONAL P(I) lQTURB*(AVHD(Z1 +Z2)/2 **820 Y1=YAV Zl=Z2 ALPHA TOTPOW=TOTPOW+P (I) ANPOWN TOTPOW*12 /24 OMEGA PRINT FORMAT POWRtMULTtANPOWNtSUM( 1)...SUM 4) CUM(1) **CUM(4) VECTOR VALUES POWR$1lHO06H MULTFS5.2,12H ANNUAL KWH,Flll/ 117H MONTHS ABOVE 65,I3/19H MONTHS ABOVE 6-45,I3/ 217H MONTHS ABOVE 64,I3/17H MONTHS ABOVE 63I13/ 317H MONTHS BELOW 56,I3/17H MONTHS BELOW 58,I3/ 417H MONTHS BELOW 6013/17H MONTHS BELOW 61,I3*$ END OF PROGRAM $ DATA 012 016 090 210 047 044 027 019 045 064 031 030 010 050 199 137 059 231 307 239 075 050 034 034 Computer 12 45 59 YAV 59.91 56.68 56.76 61,81 60*05 58*08 55,18 51.53 49.12 47.90 4-407 39.49 30.28 27.21 45 20 Output 16 64 231 92015 89510 88262 90324 91600 90126 88189 85512 82963 81411 79182 75318 68527 61909 70035 90 210 31 30 307 239 0 0 0 0 0 47 10 75 11997 15981 90001 209970 47029 44001 26972 18985 45002 64014 31008 30002 9993 49998 199002 44 27 19 50 199 137 50 34 34 I % -C39 -

Optimum Regulation of a Reservoir for Power Production Computer Output, Continued 49040 80462 47T*4 81526 46s9s 84506 6444 9:1118 14' 6C9S 93039 1 9; 46t42 92167 26' 60*61 91934 58646 90575 55.94 a643 NrLT* 1.00 AWIIIAL KWH 424 MONTHS ABOVE 65 0 MONTHS ABOVE 64*5 0 QMNTHS ABOVE 64 1 MOSTHS ABOVE 63 2 MOTHS BELOW 56 12 MOTHS BELOW 58 15 O!HTHS BELOW 60 18 MONTHS BELOW 61 20 0 3 0 0; 785 229 0 0 0 15296 ~ 0 137016 59005 231023 107008 239057 74969 49989 34027 34014 Discussion of Results Only two years (24 months) of data were used which required 2.2 and 0.4 minutes, respectively, of execution and compiling time. Because of space limitations, the printed output from the computer is shown for only one percentage of turbine operation (i.e., the case of full-time turbine operation). The results indicate the desirability of decreasing the range of turbine operation examined and increasing the number of months studied. Eight years of record could be studied in five minutes of machine time with this modification. The results revealed an unexpected advantage to an operation involving considerable drafting of the reservoir even at substantial reductions of head since such operation produced the greatest amount of energy conversion. Critique This problem has been used in a second course in Water Power Engineering but, because of the complexity of numerical work,,full solutions have not been required previously. Since the computer makes a trial and error solution for the upstream and downstream water levels feasible, the problem can now be fully solved and will be assigned to students in the future. -C40 -

Example Problem No. 78 ECONOMICAL SIZE OF CANAL FEEDING HYDRO-POWER PLANT by J. W. Howe Department of Mechanics and Hydraulics State University of Iowa Course: Water Power Engineering Credit Hours: 3 Level: Senior and Graduate Statement of the Problem A canal having 1:2 side slopes is to feed a hydro-power plant. It is to have best hydraulic efficiency (R = D/2) when carrying the mean discharge of 1500 cubic feet per second, but capacity for the maximum discharge of 2000 cubic feet per second with a 2-feet freeboard. Excavation is limited to 60 per cent of the maximum depth with the excavated material forming berms on either side to provide the upper portion of the cross-section. Transverse slope of the ground is assumed to be zero but any longitudinal slope of the canal is possible by slight divergence from a contour. The following data are given: Canal roughness coefficient, n = 0.025; excavation cost, $2.00/cubic yard; annual charge for maintenance, interest, and taxes, 10 per cent; cost of power, $0.006/KWH; operation schedule, 16 hours/day, 300 days/year. Slopes varying from 0.00020 to 0.00046 are to be used. Determine the dimensions of the crosssection for which the annual cost is least. Solution A diagram of a typical canal cross-section is shown below to illustrate the important dimensions. A = B X Y + 2Y2 __^ ^^\. R R A/P BxY + 2Y2 Y R = A/P = B + 4.4722 x Y = 1^.z~ | VY B = 0.4722 X Y For a series of slopes varying from 0.00020 to 0.00046 by 0.00002 increments, the Manning formula is solved by a trial and error process. Since the channel is to have best hydraulic efficiency at mean discharge (R = Y/2): By the Manning formula = 15 A R2/3 S1/2 - 1.5 (BXY + 2Y2)5/3 /2 N (B + 4.4722 X Y)23 -C41 -

Economical Size of Canal Feeding Hydro-Power Plant Substituting B = 0.4722Y Q 1.5 (0.4722 y2 + 2y2)5/3 1/2 2.472 2)/ _(o.472 S (2.42 N (0.4722Y + 4.4722Y)2/3 (4.944 y)23 y = Q N)3/8 B Y 0.4722 (Q N)3/8 1.375 S3 1 1.375 S3/16 Hence, there is a particular bottom width for each slope which satisfies the condition R = Y/2 for the mean discharge. The required depth for maximum discharge at each slope and corresponding bottom width are next found by trial and error starting with given minimum and maximum values of Y and using their mean as the depth. After 11 iterations the mean should have the correct value. QMAX 1.5 [B X YAV + 2 x (YAV)2] 5/3 1/2 N (B + 4.4722 YAV)2/3 L = (QMAX X N)/1.5 X S1/2 R = [B X YAV + 2(YAV)2] 5/3 / (B + 4.4722 YAV)2/3 Whenever L > R Whenever L< R YMIN = YAV YMAX = YAV Depth = 0.6 x (Final value of YMAX + 2) Annual excavation cost/foot length = C X INT X (Vol. of Exc. C X INT X [B X 0.6 DEPTH + 2(0.6 DEPTH)2] Ann. Cost of Power loss/foot length = QMEAN X S X H X D X 1 X ~74 = MEAN X S X H X D X = QMEAN X S X H X D X ~475 Total Annual Cost = Annual chargefbr Excavation plus Annual Cost of Power Loss List of Symbols QMEAN average discharge rate, ft3/sec QMAX maximum discharge rate, ft3/sec A cross-sectional area, ft2 R hydraulic radius, ft N roughness coefficient n in Manning formula S slope of bottom in direction of flow, dimensionless B bottom width of channel, ft YMIN lowest depth of water, ft YAV mean depth of water, ft -C42 -

Example Problem No. 78 List of Symbols, Continued YMAX highest depth of water, ft DEPTH total depth of canal with 2-ft freeboard C cost of excavation, $/yd3 INT annual fixed charges (interest, taxes, maintenance), dimensionless P value of power, cents/KWH H hours of operation/day, hrs D days of operation/yr, days ANCSXC annual cost of canal, $/yr ANCSPR annual cost of power lost, $/yr ANCOST total annual cost, $/yr Flow Diagram V(PH De P r H (L) - lvl\\X2 - C^2.N r )AC. ^'^.^ -"^ ^ DePT*(x)EPT) A)bCSiR G ) DN 0 A P7.) PPN(SA- 2 7 475 0 - ~ "' H() I\ ft ~ _c Irrri).bDFTH().2 D~r~r)) rI i () I XDr 0( c- (f) ( A h.5Csr ANCO'.)" =T q5& &(s) S -1 5(14) E', rI) - ( DP(' 4) ( SLO PES A >ANC>A(J) W(..AL14) 1 RT \LU i/~S ~ANC& P-)" -P!F-itC'\ 4) (ANCS Pi().. ~ NcPRL(, ) T,, — /''' -C43 -

Economical Size of Canal Feeding Hydro-Power Plant MAD Program and Data Jo We HOWE S164A 001 005 000 CANAL J. We HOWE S164A 001 005 000 CANAL S COMPILE MAD, EXECUTE, DUMPt PRINT OBJECT R RCANAL IN EARTH WITH 1/2SIDE SLOPES. MANNING ROUGHNESS N, BEST RHYDRAYLIC RADIUS (0.4722*Y FOR QMEAN),OMEAN AND QMAX IN CFS. REXCAVATED DEPTH 0.6(YMAX+2). COST OF EXCAVATION C $/CUYDe RFIXED CHARGE INT, POWER COST P CENTS/KWH, OPERATION H HOURS/ RDAY, D DAYS/YR~ ANY SLOPE AVAILABLE BY DIVERGENCE FROM RCONTOUR. TRANSVERSE SLOPE ASSUMED ZERO* R DIMENSION S(14),B(14) ANCSXC(14),ANCSPR(14) ANCOST(14), 1DEPTH(14) INTEGER IJ PRINT COMMENT $1$ READ FORMAT DATAQMEANOMAXNCINTPHD*S(1).*.S(14) VECTOR VALUES DATA=$2F1O.OFlO~.3,3F8*22F8.0/( 7F105)*S PRINT FORMAT ECHOQMEANQMAX.NCINTPHDS( 1)**.4S( 14) VECTOR VALUES ECHO = $1H 2F10s0PF1O.3s3F82t22F8.0/ 1 (1H07OFIO5)*$ THROUGH SLOPESFOR Isl, l IG,14 B(I) =.4722*( MEAN*N )P 375/( 1375*S(! ) *P 187'5) YMIN=11~ YMAXuS'18 THROUGH TRIALSFFOR J="11.tJ.G.11 YAV=tYM+YN+YMAX)/2 Ls(QMAX*N)/( 15*SQRT (S( I ) Ra(B I)*YAV+2*YAV*YAV) P ~ 1667/(B(I )+4 4722*YAV ),P *667 WHENEVER Le G.R YMIN=YAV OTHERWISE YMAX=YAV TRIALS END OF CONDITIONAL DEPTH (I) -YMAX+2 ANCSXC(I)~C*INT*(B(I)*6+*DEPTH( I)+*72*DEPTH(I)*DEPTH( ))/27 ANCSPR(I)=QMEAN*S(I)*H*D*P/1475 SLOPES ANCOST(I) mANCSXC(I)+ANCSPR (I) PRINT COMMENT $4$ PRINT FORMAT RESLTSS(1)*..S(14),B(1)*.*.B(14),DEPTH(1)*.4DEPT 1H(14) ANCSXC(1) *l ANCSXC(14) ANCSPR(1) **ANCSPR(14) ANCOST(l) 2 ~ ANCOST(14) VECTOR VALUES RESLTS = S1H 7F10.5/1H07F10.5/(1HO7F10O2)*S END OF PROGRAM $ DATA 1500. 2000. *025 2.00 ~10 *60 16. 300. 00020 *00022 *00024.00026 *00028 @00030.00032 *00034 *00036 *00038.00040 *00042.00044.00046 -c44 -

Example Problem No. 78 Computer Output 1500 000 20 *00034.00020.00034 6.60 5.98 17 74 16.25 2.20 1.84.59 1.00 2,79 2.84 2000.00022.00036.00022.00036 6.48 5.91 17.46 16.10 2.13 1.81 *64 1.05 2.77 2.86 *025.00024.00038.00024,00038 6.38 5.85 17.21 15.96 2.07 1,77 *70 1.11 2*77 2,89 2.00,10 *60 16 300,00026,00028 *00030,00032.00040,00042 *00044 00046 *00026 *00040 6.28 5.80 16.99 15.83 2.01 1074 *76 1.17 2.78 2,92. 00028.00042 6.20 5.74 16.78 15.70 1,96 1.72 *82 1.23 2.78 2.95.00030.00044 6.12 5.69 16.59 15,58 1.92 1.69.88 1.29 2.80 2.98.00032. 00046 6.04 5,65 16.42 15.47 1X88 1o66.94 1,35 2.82 3.01 } } ) Input Data Slopes Bottom Widths Depths Annual Cost of Canal Annual Cost of Power Lost Total Annual Cost * Discussion of Results A plot is shown below of the total annual cost as a function of the canal slope. It is seen that a minimum cost of $2.77 occurs for a slope of 0.00023. E-4 Co 0 0 g 14I $3.00 2.90 2.80 2.70 0.0002 0.0003 0.0004 CANAL SLOPE The problem involves 0.3 minutes of compiling time and 0.7 minutes total time. Critique This is a problem regularly given to students in a Water Power course. Without the use of the computer, the students can obtain a satisfactory but laborious solution using special tables. Plotting of computed results is usually necessary for interpolation between points. -C45 -

Economical Size of Canal Feeding Hydro-Power Plant The use of the computer for obtaining a trial and error solution replaces the tables, makes it feasible to compute many more points and affords an easy way to secure solutions for different original data with little additional work. The problem, though short, illustrates the advantage of the computer approach to an otherwise laborious problem and provides a good exercise for students who are learning to program.

Example Problem No. 79 MOST PROBABLE NUMBER METHOD FOR DETERMINING BACTERIAL POPULATIONS by L. L. Kempe Sanitary Engineering The University of Michigan Course: Industrial Bacteriology Credit hours: 3 Level: Senior or Graduate Statement of Problem Write a MAD program that will solve, using the Newton-Raphson method, the MPN (Most Probable Number) equation. The input of the program will be N, ai, ni, Pi. The output should be X, the most probable number of bacteria, and P, the probability associated with X. Solution The MPN equation is based on the extinction dilution method for finding the number of viable bacteria in a solution. The development of the MPN equation is well documented ( ) and will only be reviewed here. Assuming a negative exponential probability of growth, the probability of bacteria growing in n tubes of culture medium innoculated with "a" milliliters of a solution containing X bacteria per mililiter is nj -aX)p /-aX) n-p p:(n-p)! (1 - ea ) (ea)n If several different dilutions are used, the probability of having Pi tubes show growth out of a total of n. tubes innoculated with a. ml. each is n n ni.p i=l Pi!n ei-p i) (1) This gives an expression for the probability of ni out of pi tubes showing growth when the solution contains X bacteria per ml. Now an expression will be found for the most probable X given a set of ai, ni, Pi. To do this, it is noted that if P (the probability) is a maximum, then In P will also be a maximum. Taking the natural logarithm of (1) ln p = i l [ln (p i-Pi - (ni-Pi)aiX + Pi ln (l-e-aiX. (2) Thus ln p= ai (nia -ni) =0, (3) dX i=l i a)X i

Method for Determining Bacterial Populations which is an implicit equation in X and gives the most probable number (MPN) of bacteria per ml. of a solution used in an extinction dilution. The approach to the problem solution is to use the Newton-Raphson method to find the X that satisfies the MPN equation. Then, given this X, the probability associated with it is calculated using equation (2). The probability will give a measure of the reliability of the experiment. A low probability would say that the experimental results were not reliable. Several extra features were incorporated in the solution. These include the ability to do a single calculation or a table of calculations. An automatic estimation of the initial value of X is also included. Any number of dilutions from 3 to 9 may be used for input data. The output format is automatically adjusted to the number of dilutions. This program is much more elegant than required, but illustrates the power of a reasonably short computer program. Further information is included as remarks in the program. List of Symbols Problem Notation Program N number of different dilutions used NBR ai number of ml. used for innoculation A(I) ni total number of tubesinnoculated N(I) Pi number of tubes innoculated which show growth P(I) X most probable number of bacteria X P(%) probability associated with X PPER Important additional variables: G = MPN function GP dG dX DELX = Ax Flow Diagram -c48 -

Example Problem No. 79 Flow Diagram (continued) P R I "T X \PPEK= (e- oo; EX P i L L L (,'LI T ' - p, oe ) ) / - pp re MAD Program and Data SCOMPILE MADePRINT OBJECTPUNCH OBJECT*EXECTEDUW4P ' PN 001 R PROGRAM FOR DETERMINING BACTERIAL POPULATION RIJING THE MPN EQUATION R R RCODED BY W.J. SANDERS 4/3/62 R R R R THIS PROGRAM USES A METHOD DEVELPED BY HALVORSOI AND RZIEGLER (Jo BACT.,25*101-121(1933)) AND PROGRAMMED FOR THE RIBM 650 BY NORMAN AND KEMPE (BIOTECHNOLOGY AND BIOENGINEERING RII-2,l57-1-3(1960))e* IT INVOLVES THE COMPUTATION OF THE RiOST PROBABLE NUMBER OF VIABLE BACTERIA IN A SAMPLE OF RMATERIAL. INPUT DATA ARE THE RESULTS OF AN EXTINCTION RDILUTION EXPERIMENT AND ARE ARRANGED AS FOLLOWS R R RCARD TYPE I R R COL R R 1-4 R 5-7 R 8-10 R 11-14 R 15-22 R 23.24 R 41 42 R R RCARD TYPE II R VAR CHECK HOW NBR RAT Al P(i) N(1) All) 25.26 P(2) *~. 39,40 P(9) 43,44 N(2) *.. 57,58 N(9) R R 1-8 9-16 A(2) ** 65-72 A(9) -C49 -

Method for Determining Bacterial Populations MAD Program and Data (Continued) R R.Rt~ST NOTAT ION FOLLOWS THAT OF HALVORSON AND ZIEGLER. R NBR - NUMBER OF DILUTIONS R gN(! NUMBER OF SAMPLES OF VOLUME A(I) ML. (.l Is.t1,*. NBR) R P(I NUMBER OF FAILURES OUT OF N(I) TRIALS* FAILURE R MEANS G6ROWTH IN TUBES R R R RTHERE ARE TWO METHODS OF IN PUT OF THE A (I) THEY DEPEND RUPON THE SETTING OF THE VARIABLE HOW* R RHOW a REG fR THIS OPTION IS USED WHEN A REGULAR DILUTION RATIO R (E.G* DECIMAL) S USED I IS QE IT I LY NECESSARY TO PROVIDE R THE FOLLOWING ADDITIONAL DATA R R Al Al ) A R RAT ~ DILUTION RATIO R R THE FORtULA A lI) a Al*RAT*P.(Il-) IS THEN USED TO R COMPUTE THE A(l), CARD TYPE II IS NOT NEEDED. R RHOW ~ VAR R THIS OPTION IS USED WHEN THERE IS NO REGULAR DILUTION R RATfO* VARIABLES Al AND RAT ARE IGNORED AND AN R ADOITIONAL CARD OF TYPE II IS REQUIRE)D R R R RTHERE IS ALSO THE OPTION OF DOING A SINGLE CALCULATION OR THE RPRINTING OF A TABLE OF CALCULATIONS* THIS DEPENDS UPON THE RSETTING OF THE VARIABLE CHECK. R RCOECK B (BLANKS) R A SINGLE CALCULATION IS MADE USING THE VALUES OF N(il R AND P(1 ) ON THE CARD. R RCHECK 3 TABL R A TABLE IS PRINTED4 ONLY P(1) IS READ. THIS IS TAKEN R AS THE UPPER LIMIT FOR ALL OF THE N(I: AND P(CI. ALL R FOLLOWING NUMBERS ON THE CARD ARE IGNORED. R R R RRESTRICTIONS ON VARIABLES R NBR nUST BE AN INTEGER BETWEEN 3 AND 9 R N(1) AND P(I) ARE INTEGERS R Alt RAT. AND THE All) MUST BE FLOATING POINT. R (IF THERE IS DIFFICULTY IN OBTAINING A SOLUTION, R THE All) SHOULD BE RESCALED TO A VALUE CLOSER TO.0O) R R R RSAMPLE INPUT CARDS ARE R R VAR 3 8 87 888 Rl. l.01 R RTABLREG 3 *. 1 l 8 -C50 -

Example Problem No. 79 MAD Program and Data (continued) R R R RTHE OUTPUT OF THE PROGRAM CONSISTS OF A LISTING OF THE INPUT RVARIABLES PLUS THE VARIABLES X AND P. THE VARIABLE X IS THE RMOST PROBABLE NUMBER OF BACTERIA PER ML* OF THE RSOLUTIONF P IS THE PERCENTAGE PROBABILITY THAT THE RCOMBINATION OF THE N(I) AND P(1) WOULD HAVE RESULTED IF X RWERE THE TRUE NUMBER OF BACTERIA. R R RIT SHOULD BE NOTED THAT ALL CALCULATIONS THAT YIELD A P RLESS THAN *01 ARE SUPPRESSED WHEN A TABLE IS PRINTED. R R R DIMENSION A( 10)P(10 )PI(10) *N(10) *NI(10) INTEGER CHECK, HOW, NBR.PI *INICOUNTtRET.KEY BOOLEAN ZEROSeTABL VECTOR VALUES BCD( 3 )S 36 S $4S8SS 60S5$72S S84S,$965.$108$ EXECUTE SETERRe START) R RREAD IN CARD TYPE I R START READ FORMAT IN. CHECK HOW,* NIBR.RAT.:Al*PI( 1})..**PI(9)NI( 1) 1. e*,N (9) VECTOR VALUES INMSC4,C3t13,F4.3.F8*6t181I2*S R RCHECK NBR R WHENEVER 1BR.L 3 OR NBR.G4 9 PRINT CO4MENTS4NBR NOT IN RANGES TRANSFER TO START R RREGULAR DILUTI ON RATIO R OR WHENEVER HOW*E*SREGS THROUGH OMEGA FOR Il1 t I*G.NBR OMEGA A(I)"Al*RATP ( I-1) PRINT FORMAT REG6NBRtAlRAT VECTOR VALUES REGaSx20H4NUMBER OF GROUPS IS.13t914 A(1) IS. 1F8.4*2H */18HODILUTION RATIO ISF6.4.2H **S R RVARIABLE DILUTIONS R OR WHENEVER HOW E.:,SVAR$ R RREAD IN CARD TYPE Ii R READ FORMAT IN1l A( 1):.:..A(NBR) VECTOR VALUES INl$=9F8.2*S PRINT FORMAT VAR*A:( l)+,eA(NBR) VECTOR VALUES VAR=S2H4A*S3,9E12 3*S OTHERWISE PRINT COMMENTS4INPUT ERRORS EXECUTE SYSTEM. END OF CONDITI ONAL R RPOSITION OUTPUT FORMAT R SLIDE(2)aBCD(NBR) OUT(2 )BCD(NBR) PRINT FORMAT SLIDE VECTOR VALUES SLIDE=$1H /1H+tS9,S,1HX.S695HP 0/O*S R RSINGLE CALCULATION R WHENEVER CHECK~E.S S TABL-OB RET-O -C51 -

Method for Determining Bacterial Populations MAD Program and Data (continued) R RFLOAT N(i- AND P(I R THROUGH PSI.FOR I-s1-.l: Ie6*BR PSI p(-)PI(I) PRINT FORMAT NOUT.*NI ( 1* )*NI(NBR) VECTOR VALUES NOUTS$2H NS399t I17'5S5)*$ TRANSFER TO 6RIND R R TABLE R OR WHENEVER CHECK*,E.4STALS TABL1 B RET-1 P a PI1l) EXECUTE SPRAY I(P,*N(l) 1 *N(NBR. tP:( 1 ).is,.*P( NBRi ) PRINT FORMAT NOUT.PI (1) PRINT COMNENTSOS R RCALCULATION OF P( I R START( 1 ) THROUGH CHIFOR KEY:NBR'-1 KEY.L. 1 P('KEY) wP(KEY)-1. R RALL P(I) O I.E~!* FINISHED R WHENEVER ZEROS. (P*NBR)AMD.-P ( 2 ). *E 0*.AND P ( 1 )1.E.0 * ITRANSPER TO START WHENEVER P KEY) *GE 0. R RFIX P( I FOR PRINTING R THROUGH PHI IFOR I 'I * I4.&.NBR PHt PI(I)-P(I) TRANSFER TO GRIND END OF CONDITIONAL CHI P (KEY) P OTHERWISE PRINT CO#IENTS4INPUT ERRORS EXECUTE SYSTEM. END OF CONDITIONAL R RCALCULATION OF INITIAL GUESS FOR X R GRIND THROUGH UPSe FOR I=l1.1IE.GNBR WHENEVER P( I)**0&.AND*P( I *L N( I) X.f( At NBR/2+1 ) * SP ( 1l.1/t;(N ( I I.-*A ( I ) TRANSFER TO TAU UPS END OF CONDITIONAL TAU DELX=1. R RNEWTON-RAPHSON SOLUTION FOR X R THROUGH ALPHAt FOR COUNT'l lCOUNT*e.T&50,OR..ABS*DELX.L...005 GP-O0 THROUGH BETA. FOR I1.l. I.'.NBR EXPON11*-EXP* (-A (1)*X) sG& A 1P 1 I /EXPN-N( I )) BETA GPGP-A ( I *A tA)*P(1 t ( l-EXPON)/ EXPONEXPON ) DELXG/GP X-X-DELX R R(SOLUTIO UNSTABLE FOR X LESS THAN 0) -C52 -

Example Problem No. 79 MAD Program and Data (Continued) R ALPHA WHENEVER X*L*0*.X=.000000001 R RITERATION LIMIT EXCEEDED WITHOUT REDUCING DELX TO LESS THAN R 60005 R WHENEVER COUNT wGe 50 PRINT FORMAT CMPTBLtPI (1)-...ePI (N R) PRINT FORMAT XCED VECTOR VALUES XCED=S1H,S113913HITER EXCEEDED*$ TRANSFER TO START(RET) END OF COND.ITIONAL R RCtOPUTATION OF P R LOGEPtO THROUGH GAMIMA FOR I l-.l ItG:NBR GAMMA LOGEPaLOGEP+ELOG ( B I NO t 14(N (Z t P ( I I i -A (.) * (N I(:.)-P (I Xt. X 1+P( I )*ELOG. ( 1,-EXPe (-A( I *X) PPER= OO *EXP (LOGEP) WHENEVER PPERsL. O01.AND*TABLtTRANSFER TO START(RET) PRINT FORMAT CMPTBLtPI ( 1):*PI(NBR) VECTOR VALUES CMPTBL=$1H /2H+P,S3,9(I7tS5)*S R RPRINT ANSWERS R PRINT FORMAT OUT:*X#PPER VECTOR VALUES OUTSIH *S4.S,E10*3tF7.2*$ TRANSFER TO STARTIRET) END OF PROGRAM $COIPILE MAD-PRINT OBJECT PUNCH OBJECT ZEROS001 R RPROGRAM FOR CHECKING TO SEE IF P(2),.**..*P(NBR) ARE ZERO R R R EXTERNAL FUNCTION (PNBR) INTEG4ER I NBR ENTRY TO ZEROS* THROUGH ALPHA FOR 1=21.-1I. eGNBR ALPHA WHENEVER P( I,.G.0.tFUNCTION RETURN OB FUNCTION RETURN lB END OF FUNCTION END OF PROGRAM $COMPILE MADtPRINT OBJECTTPUNCH OBJECT R RPROGRAM FOR COMPUTING THE BINOMIAL COEFFICIENT R R R BINOMOOl ALPHA BETA $DATA VAR 1. TABLREG EXTERNAL FUNCTION (UP.DWN) ENTRY TO BINOM. FACT=1l THROUGH ALPHAt FOR TOP=DWN+1 t,1,TOP.G4UP FACT=FACT*TOP THROUGH BETA. FOR BOT=2*l18BOT*G (UP-DWN) FACT=FACT/BOT FUNCTION RETURN FACT END OF FUNCTION END OF PROGRAM 3 1 *01 3 *1 1. 887 888 3 -C53 -

Method for Determining Bacterial Populations Computer Output A N P 100E 01 100OE-01 1OOOE-02 8 8 8 8 8 7 X P 0/0 *208E 03 39.27 NUMBER OF GROUPS IS 3* Al ) IS DILUTION RATIO IS.1000 N 3 P 3 3 P 3 3 P 3 3 P 3 2 P 3 2 P 3 2 P 3 2 P 3 1 P 31 P 3 1 P 3 1 P 3 0 P 3 0 P 2 3 P 2 3 P 2 2 P 2 2 P 2 1 P 2 1 P 2 1 P 2 0 P 2 0 P 2 0 P 3 P 1 2 1,0000. x X 2 *110E 03 1 *462E 02 0 *240E 02 3.292E 02 2.215E 02 1,149E 02 0 *933E 01 3 *159E 02 2 *115E 02 1 *749E 01 0.427E 01 2 *636E 01 1 *3-5E 01 1 *360E 01 0 *286E 01 1 *276E 01 0 *211E 01 2 *268E 01 1 *205E 01 0 147E 01 2 *199E 01 1 *143E 01 0 9-18E 00 0 *157E 01 1 *154E 01 P 0/0 44.44 42 77 36 59 *24 2.47 12. 51 32 82.03 *65 6.56 37 43 *16 3*10 *02 *22 *17 2 32 *02 *63 11 96 *02 1*12 31.93 *03 *03 -C54 -

Example Problem No. 79 Computer Output (continued) P 1 2 0 *114E 01 63 P 1 1 1 *112E 01 18 P 1 1 0 *7i6E 00 6.45 P 1 0 1 *72E 00 *62 P 1 0 0 *35E 00 39*20 P 0 2 0 *619E 00 *16 P 0 1 1 *11E 00 *05 p 0 1 0 *305E 00 3.36 Discussion of Results The output shown above was designed to be illustrative rather than useful. In practice more than three tubes would be used in each dilution group. Compiling time for the program was 1.1 minutes. Computing the results shown above took 0.5 minutes. Critique This problem is very appropriate for student solution. It gives the student an opportunity to become familiar with numerical solution of complex equations as well as an opportunity to become familiar with the use of the MPN equation. Solution of the MPN equation by hand is extremely tedious, so that the tables prepared by the instructor's program will be useful in later work. References 1. Halvorson, H. 0. and Ziegler, N. R. J. Bact., 25, 101-121 (1933). 2. Norman, R. L. and Kempe, L. L., Biotechnology and Bioengineering, II, No. 2, 157-163 (1960). -C55 -

Example Problem No. 80 ITERATION OF SUM OF PRINCIPAL STRESSES AT INTERIOR POINTS OF PLANE RINGS by William Weaver, Jr. Department of Civil Engineering Stanford University Course: Experimental Stress Analysis Credit Hours: 3 Level: Graduate Statement of the Problem Write a computer program which will evaluate the sum of principal stresses at interior points of plane rings for which the boundary stresses are known from photoelastic analysis. The recommended procedure is to iterate the Laplace difference equation applied to an x-y grid network of approximately square elements obtained by the conformal transformation of the R-theta polar network using the relationships: x = In (R) y = theta The program should be flexible enough to accommodate a whole ring or any segment of a ring for which the stresses on the radial boundaries are known to be equal by symmetry. Divide the central angle of the polar network into 5~ increments. Make the elements of the transformed rectangular grid as nearly square as possible. The number of radial increments will depend upon the ratio of the outer and the inner radii of the ring. Limit the number of storage locations required by the grid to about 3,000. Iterate the sum of principal stresses to within 3 decimal point accuracy. Limit the number of iteration cycles to 50. Solution The Laplace difference equation for the sum of principal stresses at a lattice point 0 and its four neighboring points 1, 2, 3, 4 (see Fig. 1) may be written: S3-S0 SO-S1 S2-S0 S0-S4 Ax3 Axl AY2 AY4 + = o (1) Ax1 + Ax3 AY2 + AY4 2 2 -C56 -

For the present problem let Ax1 = Ax3 = PAY2 = PAY4 Solving for SO in equation (1), 1 2 5\ P S - (s + + + S4) (2) 0 2(p2+1) (2 3 2 (S2 4) +1) +2+ / + s(2 in which, S = Sum of principal stresses at a point p = Ratio of x-intervals to y-intervals. Equation (2) may be used for the iterative evaluation of S at the lattice points of the transformed grid. In the photoelastic analysis of plane rings the principal stresses at the inner and outer edges are determined experimentally at even intervals of the central angle. A practical minimum interval is 5~, or r/36 radians. By conformal transformation: x = in (R); y = 9 Let Ay = AG = T /36. Define M as the number of angular increments, i.e., increments in the y-direction. M is given by M total angle of ring AG Thus, for a full ring and AG = 5~, M = 3600/5~ = 72. For a quarter ring M = 90~/5~ = 18. In order to obtain an exactly square grid, the number of subdivisions in the radial direction must be: xN-x0 ln(RN)-ln(RO) 36 - = - ^ - -A in(RN/RO) N1 = /36 = /36 - in which, R = Inner radius of ring. RN = Outer radius of ring. x = Transformed inner radius. xN = Transformed outer radius. In general, N1 will not be an integer, but an approximately square grid may be obtained by computing the actual number of intervals N used in the radial direction as the next higher integer value of N1. -C57 -

Iteration of Sum of Principal Stresses N = Xln(RN R)] Rounded up The dimensionless ratio, p, may be computed from: Ax N1 36 R P = y - = n(RN/R) (4) Radii of interior network circles are determined by the relationship: (x ) (xo+jAx) ln(RN/R) j/N= R(RN/R) Rj = e = e R(e )N = Ro(RN/RO) (5) in which, Rj = Radius of jth network circle. j = Radial index. Assuming that a storage space of about 3,000 locations is available for the M x (N+1) matrix of lattice points required for the solution, M X (N+l) _ 3000. (6) For a full ring the maximum value of M is 3600/50 = 72, and the value of N is given by substitution of equation (3) into equation (6). 72{ [ l(N Rounded u + I1 3000 (7) Equation (7) yields the maximum values of N and the corresponding ratio of radii. (N)max = 41; (RN/R)max 35.0 The minimum non-trivial values for these parameters correspond to a network of only two intervals in the radial direction. (N)min = 2; (RN/RO)min 1 19 Figures 2 and 3 show the polar and transformed lattic networks for a quarter ring having RN/R0 = 4, M = 18, and N = 16. This configuration, accompanied by a set of hypothetical boundary stresses, constitutes the first of two trial problems submitted to the computer. The second problem, for which no figure is given, consists of a quarter ring having RN/R0 = 4/3, M = 18, and N = 4. The first row of the transformed grid may be omitted because the stresses must be the same as in the last row. -C58 -

Example Problem No. 80 Flow Diagram M,)lINNER (i) P/REST I NDEX( R R-)R) a TITLE A INN^ER(,(M) vINE(I) ReRNeUTER(1) \I>M E I>M / iT~u(r'.- *UTER(e ) N= Nf + f F t F,2 P PZ F I zZO) R T)= R e\ Ri -b 5 (Il,J ) = INNFR(I) " I TER' z= If1= J = (FJ/IF.NNER(I) N 5WESWEPO 5EP=5WEEP+t ' \,,>l - eUTER(Z)) J COUNT '= 0 TEMP= FI(S(1,J-I)+S(i,J-+I)) +FZ(S(M, J) +S (,)) --- FE~R Tre) --- D Special Output Form when N < 10: INNER(I) = Boundary stress at inner edge. OUTER(I) = Boundary stress at outer edge. Fl = Weighting factor from first term in equation (2). F2 = Weighting factor from second term in equation (2). FJ = J expressed as a floating point number. FN = N expressed as a floating point number. -C59 -

Iteration of Sum of Principal Stresses MAD Program W. WEAVER D046N 005 005 000 W. WEAVER D046N 005 005 000 $COMPILE MAD, EXECUTE R RTHIS PROGRAM IS INTENDED FOR USE IN EVALUATING THE RSUM OF PRINCIPAL STRESSES- AT INTERIOR POINTS OF RPLANE RINGS ANALYZED EXPERIMENTALLY BY THE RMETHOD OF PHOTOELASTICITY* THE PROCEDURE CONSISTS OF RTHE ITERATION OF THE LAPLACE DIFFERENCE EQUATION RAPPLIED TO AN X-Y GRID NETWORK OF APPROXIMATELY RSQUARE ELEMENTS OBTAINED BY THE CONFORMAL RTRANSFORMATION OF THE R-THETA NETWORK USING RTHE RELATIONSHIPS R X=LN(R) R Y=THETA RTHE PROGRAM IS APPLICABLE FOR A WHOLE RING OR RANY SEGMENT THEREOF,THE RADIAL BOUNDARIES OF RWHICH LIE MIDWAY BETWEEN ADJACENT RADII OF SYMMETRY. RTHE CENTRAL ANGLE IS DIVIDED INTO FIVE DEGREE RINCREMENTStRESULTING IN M=72 INCREMENTS FOR A RFULL RINGM=36 INCREMENTS FOR A HALF RINGM=18 RFOR A THIRD RINGETC. THE NUMBER OF RADIAL INCRERMENTS DEPENDS UPON THE RATIO OF THE OUTER AND THE RINNER RADII OF THE RING. IN ORDER TO LIMIT THE NUMBER ROF STORAGE LOCATIONS REQUIRED BY THE GRID TO ABOUT R3,000, THE MAXIMUM RATIO OF RADII PERMITTED IS 35. RTHE MINIMUM NON-TRIVIAL RATIO IS 1.19. RTHE CALCULATIONS ARE CARRIED TO 3 DECIMAL PLACE ACCURACY ROR 50 ITERATION CYCLES, WHICHEVER OCCURS FIRST. R PRINT FORMAT TITLE VECTOR VALUES TITLE=$1H1,29HITERATION OF SUM OF PRINCIPAL/ 11HO.42HSTRESSES AT INTERIOR POINTS OF PLANE RINGS*$ R RREAD AND ECHO CHECK DATA R START READ FORMAT DATA1,M,INNER(1)...INNER(M) READ FORMAT DATA2,RO.RN,OUTER(1)...OUTER(M) VECTOR VALUES DATA1=$I2,S18S5F10.3/(7F10.3)*$ VECTOR VALUES DATA2=$2F10.45F10,3/(7F10.3)*$ THROUGH PRESET, FOR I=0,1,I.G.M PRESET INDEX(I)=I PRINT FORMAT DDATA1 PUNCH FORMAT DDATA1 VECTOR VALUES DDATA1-$HO5,4HDATA515t5HINDEXtS555HINNERt55 1,5HOUTER*$ THROUGH ALPHA, FOR I=1,lI.*G.M ALPHA PRINT FORMAT DDATA2,INDEX(I),INNER(I),OUTER(I) THROUGH ALPHA1i FOR I=1l1,t.G.M ALPHA1 PUNCH FORMAT DDATA2,INDEX(I),INNER(I),OUTER(I) VECTOR VALUES DDATA2=$1H,S27,I2,2F10.3*$ R RCOMPUTE ITERATION FACTORS R Nl(36./3.14159)*ELOG.(RN/RO) NN1l+1 Pz(36.*ELOG*(RN/RO))/(3.14159*N) F11.o/(2.*(P*P+1.)) F2zP*P*F1 R RCOMPUTE INTERMEDIATE RADII R(USING FLOATING POINT NF=N AND JF=J) R NF=N THROUGH BETA, FOR J=0,1,J.G N JF=J BETA R(J)xRO*EXP.((JF/NF)*ELOG.(RN/RO)) R RSET STARTING VALUES FOR S USING LINEAR APPROXIMATION RTO ESTIMATE INTERIOR VALUES BASED ON BOUNDARY VALUES, R(USING FLOATING POINT NF=N AND JF=J) R DIMS(2)}N+1 THROUGH GAMMA, FOR Il.tltI.G.M THROUGH GAMMA, FOR JOltJ*G.N JFzJ -c6o

Example Problem No. 80 MAD Program, Continued GAMMA S( IJ)=INNER( I )-(JF/NF)*( INNER( I )-OUTER( I)) R RINITIALIZE COUNTERS FOR NO. OF CYCLES AND ACCURACY R SWEEP=O R RITERATE FIRST ROW*USING VALUES FROM FIRSTLASTAND RSECOND ROWS. R ITER SWEEP=SWEEP+1 COUNT-O THROUGH DELTA, FOR J=l1lJ.G.(N-1) TEMP=Fl*(S(1,J-1 )+S(1 J+ )+F2*(S(MJ)+S(2,J)) WHENEVER.ABS (S(1,J)-TEMP )L.O. OO1COUNTCOUNT+l DELTA S(1,J)=TEMP R RITERATE INTERMEDIATE ROWS R THROUGH EPSIL, FOR I=2,1.I.G*(M-1) THROUGH EPSIL, FOR J=11,iJ.G.(N-1) TEMP=Fl*(S( IJ-1)+S( IJ+1))+F2*(S( I-1J)+S(i+1,J)) WHENEVER.ABS.(S( I,J)-TEMP).L.0*001COUNT=COUNT+l EPSIL S(IJ)=TEMP R RITERATE LAST ROWIUSING VALUES FROM FIRSTLASTAND NEXTRTO-LAST ROWS. R THROUGH ZETA. FOR J=ltlJ.G,(N-1) TEMP=F1*(S(MJ-1)+S(MJ+1))+F2*(S(M-1,J)+S( 1 )) WHENEVER.ABS.(S( MJ)-TEMP ) L.O0OO1COUNTzCOUNT+l ZETA S(MtJ)=TEMP R RTEST FOR NO.OF CYCLES AND ACCURACY R WHENEVER SWEEP.GE.50,TRANSFER TO QUIT WHENEVER COUNT.L*(M*(N-1)),TRANSFER TO ITER QUIT WHENEVER N.L.10 TRANSFER TO ANS2 OTHERWISE R ROUTPUT FORM FOR N GREATER THAN OR EQUAL TO 10. R PRINT FORMAT FORM1,R(Q)..*R(N) VECTOR VALUES FORM1=$1HO,11HVALUES OF R/(7F10.4)*$ PRINT FORMAT FORM1A.M,N VECTOR VALUES FORM1A=-1HO,15HITERATED VALUESS10,2HMI,12,S10O 12HN=*I2//*$ PRINT FORMAT FORM1BS(10,)..*S(MtN) VECTOR VALUES FORM1B=$(7F10.3)*$ PRINT FORMAT FORMlCSWEEPCOUNTM*(N-l) VECTOR VALUES FORM1C=$1HO~14HNO, OF CYCLES=,SllI2/1HO023HN0. 1 OF ACCURATE VALUES=,I4/1HO,23HNO. OF ITERATED VALUES=,I4*$ END OF CONDITIONAL TRANSFER TO START R RSPECIAL OUTPUT FORM FOR N LESS THAN 10. R ANS2 PRINT FORMAT FORM2,MNR(O)...R(N) VECTOR VALUES FORM2=$IHO,15HITERATED VALUES,S102HM=,I2,S1l0 12HN=:I2/1HOS3'6HVALUES /1HOQS522HOFS7,2HR=,S3,lOF10.4*$ PRINT FORMAT FORM2AINDEX(O),.*INDEX(N) VECTOR VALUES FORM2A=$1HOS4,5HTHETAS3,7HINDICES, 1S6,I1,9I10*$ THROUGH ETA. FOR I=1lltI.G.M ETA PRINT FORMAT FORM2B,5*INDEX(I),INDEX(I),S(IO)...S(I*N) VECTOR VALUES FORM2B=$55tI3~S7,I2,S3,1OF1O.3*$ PRINT FORMAT FORM2C9SWEEPCOUNTM*(N-1) VECTOR VALUES FORM2C=$1HO,14HNO. OF CYCLES=S11,I2/1HO, 123HNO. OF ACCURATE VALUES=vI4/lHO,23HNO. OF ITERATED VALUES= 2I4*$ TRANSFER TO START INTEGER M,NIJ9SWEEPCOUNTINDEX DIMENSION INNER(72),OUTER(72),INDEX(72),R(42) DIMENSION S(3000,DIMS) VECTOR VALUES DIMS=2,1,42 END OF PROGRAM -C61 -

Iteration of Sum of Principal Stresses Data $DATA 18 -5000 -3222 -1800 -457 613 1711 2641 3416 3800 4000 3800 3416 2641 1711 613 -457 -1800 -3222 15000 60000 3000 7 -2697 -5002 -6896 -8005 -8985 -9625 -9900 -10000 -9900 -9625 -8985 -8005 -6896 -5002 -2697 7 18 -5000 -3222 -1800 -457 613 1711 2641 3416 3800 4000 3800 3416 2641 1711 613 -457 -1800 -3222 45000 60000 3000 7 -2697 -5002 -6896 -8005 -8985 -9625 -9900 -10000 -9900 -9625 -8985 -8005 -6896 -5002 -2697 7 Computer Output DATA INDEX INNER OUTER 1 -5 000 3.000 2 -3.222 0.007 3 -1.800 -2.697 4 -0O457 -5.002 5 0.613 -6.896 6 1.711 -8.005 7 2.641 -8.985 8 3.416 -9.625 9 3.800 -9.900 10 4*000 -10.000 11 3.800 -9.900 12 3 416 -9.625 13 2.641 -8.985 14 1.711 -8.005 15 0.613 -6.896 16 -0.457 -5.002 17 -1.800 -2.697 18 -3.222 0.007 VALUES OF R 1.5000 1.6358 1.7838 1.9453 2.1213 2.3133 2.5227 2.7510 3.0000 3.2715 3.5676 3.8905 4.2426 4.6266 5.0454 5.5020 6.0000 ITERATED VALUES M=18 N=16 -5.000 -3.075 -2.290 -1.982 -1.920 -1.999 -2.160 -2.365 -2.584 -2.790 -2 953 -3.034 -2.977 -2.688 -1.988 -0.471 3.000 -3.222 -2.498 -2.049 -1.860 -1.852 -1.961 -2.142 -2.362 -2.595 -2.816 -2.998 -3 108 -3.098 -2.899 -2*404 -1.464 0.007 -1.800 -1.642 -1.545 -1 556 -1.666 -1.851 -2,085 -2.348 -2.620 -2.883 -3.119 -3*304 -3.412 -3.411 -3.271 -2.994 -2.697 -0.457 -0.725 -0.933 -1.153 -1.406 -1.691 -2.001 -2.325 -2.654 -2.979 -3.290 -3.578 -3.836 -4.063 -4.278 -4.543 -5.002 0.613 0.132 -0.307 -0.717 -1.114 -1.507 -1.901 -2.297 -2.693 -3.089 -3.484 -3 882 -4.290 -4.727 -5.234 -5 897 -6.896 1.711 0.945 0*288 -0.293 -0.824 -1.321 -1 799 -2.266 -2.730 -3.197 -3.674 -4.174 -4.713 -5.318 -6.031 -6.909 -8.005 2.641 1.649 0.807 0.079 -0.566 -1.154 -1.706 -2.238 -2,762 -3 291 -3.840 -4.425 -5.068 -5.799 -6.659 -7.699 -8.985 3.416 2.202 1.210 0.370 -0.363 -1.022 -1.632 -2.213 -2.785 -3.363 -3 966 -4.614 -5.332 -6.149 -7.103 -8.241 -9*625 3.800 2.529 1.461 0.554 -0.233 -0.936 -1.583 -2.196 -2.798 -3.407 -4.044 -4.730 -5.491 -6.357 -7 360 -8.533 -9.900 4.000 2.653 1.549 0.618 -0.188 -0*905 -1.564 -2.188 -2.800 -3.420 -4.068 -4 767 -5.544 -6.425 -7 443 -8.628 -10.000 3.800 2 530 1.463 0.557 -0.230 -0 931 -1.577 -2.190 -2.791 -3.401 -4. 038 -4.724 -5.487 -6.354 -7.358 -8.532 -9.900 3.416 2.204 1.214 0.376 -0.355 -1.013 -1.621 -2.202 -2.773 -3.351 -3.955 -4.604 -5.323 -6.142 -7.098 -8.239 -9.625 2.641 1.652 0.812 0.087 -0.556 -1.142 -1.692 -2.222 -2.746 -3 275 -3.825 -4.411 -5.056 -5.790 -6.653 -7.696 -8.985 1.711 0.948 0.294 -0.285 -0.812 -1.308 -1.783 -2.249 -2.712 -3.178 -062 -

Example Problem No. 80 Computer Output, Continued -3.657 0.613 -2.280 -5.227 -1.397 -3.564 -1.640 -2*609 -2 991 -1.957 -3.094 -4.158 0.135 -2.675 -5.893 -1.680 -3.824 -1.542 -2.872 -2.697 -2.137 -2.896 -4.700 -0.302 -3.071 -6.896 -1.987 -4.054 -1.551 -3.108 -3.222 -2.357 -2.402 -5.308 -0.708 -3 467 -0.457 -2.310 -4.272 -1.660 -3.294 -2.497 -2.590 -1 462 -6*024 -1.102 -3.866 -0.722 -2.639 -4.540 -1.843 -3.403 -2.047 -2.811 0.007 -6.905 -1.494 -4.276 -0.928 -2.964 -5.002 -2. 076 -3.404 -1.858 -2.993 -8.005 -1.886 -4 716 -1.146 -3.275 -1.800 -2.337 -3.267 -1.849 -3.103 NO. OF CYCLES= 50 NO. OF ACCURATE VALUES= 211 NO. OF ITERATED VALUES= 270 DATA INDEX 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 INNER -5.000 -3.222 -1.800 -0.457 0.613 1.711 2.641 3.416 3.800 4.000 3.800 3.416 2.641 1.711 0.613 -0.457 -1.800 -3.222 OUTER 3.000 0.007 -2.697 -5.002 -6.896 -8.005 -8.985 -9.625 -9.900 -10.000 -9.900 -9.625 -8.985 -8.005 -6,896 -5.002 -2.697 0.007 ITERATED VALUES M=18 N= 4 VALUES OF R 4.5000 4.8356 THETA INDICES 0 1 5 1 -5.000 -2.999 10 2 -3.222 -2.524 15 3 -1.800 -2 023 20 4 -0.457 -1 551 25 5 0.613 -1 132 30 6 1.711 -0.675 35 7 2.641 -0.248 40 8 3.416 0.114 45 9 3.800 0.326 50 10 4.000 0.414 55 11 3.800 0.326 60 12 3.416 0.113 65 13 2.641 -0.248 70 14 1.711 -0.674 75 15 0.613 -1.130 80 16 -0.457 -1 549 85 17 -1.800 -2 022 90 18 -3.222 -2 524 5.1962 2 -1.644 -1.845 -2.226 -2.607 -2.901 -3.039 -3.092 -3.086 -3.064 -3.050 -3.064 -3.086 -3.092 -3.039 -2.899 -2.606 -2.225 -1.844 5.5836 6.000 3 -0.018 -1.042 -2.430 -3.723 -4.775 -5.463 -5.974 -6.297 -6*447 -6.493 -6.447 -6.297 -5.974 -5.462 -4.774 -3.722 -2.429 -1.042 4 3.000 0.007 -2.697 -5.002 -6.896 -8.005 -8.985 -9.625 -9.900 -10.000 -9.900 -9.625 -8.985 -8.005 -6.896 -5.002 -2.697 0.007 11 NO. OF CYCLES= NO. OF ACCURATE VALUES= 54 NO. OF ITERATED VALUES= 54 -C63 -

Iteration of Sum of Principal Stresses Discussion of Results The first problem submitted to the computer involved 270 iterations in each cycle. After 50 cycles, 211, (or about 78%)of the values were accurate to three decimal places. Because N was large, values of R and S were printed sequentially according to the general format labeled QUIT. The second problem involved 54 iterations in each cycle, and all values were accurate to three decimal places after 11 cycles. Because N was, in this case, less than 10, values of R and S were printed in rows and columns corresponding to Figure 3, according to the special format labeled ANS2. Critique This problem is very well suited for digital computation. The analysis of plane rings is a recurring problem in photoelastic experiments, and the iteration of the sum of principal stresses at interior points is an essential, though tedious, step in the total analysis. The fineness of the lattice network is not significantly limited by the storage space of a computer of reasonable size. A more important factor in limiting fineness is the practical difficulty of measuring boundary stresses accurately at the inner edges of plane rings for small increments of the central angle. Increments of 5~ are therefore used as a practical minimum in establishing the lattice network. x yi 4 - Y4 ) 1 Axl 0 Ax3 t )3 AY2 2 ---------- ~~~~ -------- - - 2- - -- Figure 1. Lattice Point 0 and its Four Neighboring Points 1,2,3,4 -C64 -

Example Problem No. 80 o o O ') /c~~~~~ 0 /3 /0 7 S Figure 2. 4 3 2 /o 0 +R~0 0 2 3, 5' G 7 8 9 / /3 2 13 /6 I /- /."6,0 7 Figure 2. Quarter Ring: RN/R = 4 M = 18 N = 16 -C65 -

Iteration of Sum of Principal Stresses 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 N 0 0 1 2 3 4 5 6 7 I i 8 9 10 11 12 13 14 15 16 17 M 18 0 00 50 10~ 15~ 200 25~ 300 350 400 450 50o 550 60~ 65~ 70~ 75~ 800 850 90~ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Figure 3. Transformed Quarter Ring RN/R = 4 M = 18 N = 16 -C66 -

Example Problem No. 81 ELASTIC BUCKLING LOAD FOR COLUMNS OF NON-UNIFORM CROSS SFCTION by K. H. Chu Civil Engineering Department Illinois Institute of Technology Course: Advanced Structural Theory and Design Credit hours: 3 Level: Seniors or First Year Graduates Statement of the Problem The differential equation for elastic buckling of a column as shown below is as follows: P IY 2 d Y M = P*Y P*_y dX E*IA E*IA E*IO*I (1] where X is the coordinate along the longitudinal axis of the column, Y is the deflection of the column, M is the bending moment, P is the axial load, E is the modulus of elasticity, IA is the moment of inertia of the cross-section about the Z axis normal to the XY plane at the point under consideration, IO is the same at a certain reference point and I = IA/IO. The buckling load, PCR, is defined as the smallest value of P which will satisfy the above equation with appropriate boundary conditions. The value of PCR is to be determined by an iterative procedure using a digital computer. Solution Equation (1) is analogous to the following equation =W (2) dX where W is the distributed load on a beam whose longitudinal axis coincides with the X axis and M' is the bending moment due to W. W is positive when it -C67 -

Elastic Buckling Load for Columns of Non-Uniform Cross Section is acting upward. M' is positive when it causes compression in the upper fiber Qf a beam. If YQ is a solution of Eq. (1), then using - as elastic weight W, the bending moment obtained will be equal to M' = YR - dX*dX + A*X + B (3) where M' is renamed as YR and A and B are arbitrary constants determined by the boundary conditions. But, from Eq. (1) PCR YQ -EI * YR (4a) Therefore YQ _ PCR YR -E*IO (4b) or PC-R= - *E*IO (4c) The procedure of obtaining the buckling load then is as follows: The column length L is divided into N equal spacings* of length A. Assume an initial set of deflectionsYO(J) for each point at X=J*A. Let I(J) be the I at the same point. Then take - O(J) as the distributed elastic weight W at the point X=J*^. Obtain the bending moment Y(J) at each of the equally spaced points. If the ratio YO(J) =C for 0 < J < N (5) Y(J) where C is a constant (within a tolerable limit) then PCR = C*E*IO (6) If not, then take the Y(J)'S as the new values for YO(J)'S and repeat the procedure. Note that the real deflection Y(J) due to P=l is given by the following equation: Y(J) =P* Y = (J) (for P = 1) (7) E*IO E*I0 * Spacings do not have to be equal, but computations will be simpler if equal spacing is used. -C68 -

Example Problem No. 81 The problem is essentially an eigenvalue problem and it can be proved mathematically that the procedure will always converge to the eigenvalue which corresponds to the smallest value of P, i.e. PCR. Newmark (reference 1) has shown that the numerical integration of equations (1) and (2) can be facilitated by the following method. If W is a distributed load, then we may pass a parabola through three points of W, say W(1), W(2) and W(3), which are equally spaced at a distance A and convert the distributed load to concentrated loads by the formulas as shown below. CiCW I C2,W(3) CW() CW(2) W(3) CW() = [7*W(1) + 6*W(2) - W(5) * (8a) CW(2) = [(1) + 10*W(2) +W(5)] * (8b) CW(3) = [7*W(3) + 6*W(2) -W(l)1 (8c) The shear and moment due to a system of concentrated loads at equally spaced points can be easily obtained without further comment. The solution to the problem consists of writing a program using the above procedure which will determine the elastic buckling load of a column with simply supported ends and with either uniform or non-uniform cross-section. For columns of nonuniform cross-sections, the cross-sections may be varying continuously or discontinuously with steps occurring at any of the equally spaced points along the column. If the discontinuity occurs between two of the neighboring equally spaced points, the column will be treated as if the cross-section varies continuously between these two points. It should be noted that the solution technique, with slight modification, can be applied to a cantilever column with one end free and the other end fixed. For example, the cantilever column and the simply supported column as shown below will each have the same buckling load. -c69 -

Elastic Buckling Load for Columns of Non-Uniform Cross Section (b) 4 LI i LI ] I L r-... The following suggestions may be helpful in programming the solution. Let us define the following quantities: KY(J) Y(J) *E*I Y( (9) 2 N) 2 2 "' 2 N L2*YO() L2*YO( L2*YO() N2*YO(() E*IO*Y() Y(N) 2*Ky ( N) KY () T N2O(J) (11) (CPCR) *KY(N) Let "EPS" be a preassigned tolerance value. Then the condition for constancy of the ratio 0) (Eq. (5)) may be restated as follows: Y(J) IT - 11 < EPS and T >0 (for all N-l points) (with 0 <c J<c N ) (12) The formula for the critical load (Eq. (6)) may be rewritten into the following form PCR = (CPCR)* E*IO (13) L2 If the number of points satisfying the above condition (Eq. (12)) is Q, then if Q < N - 1, we shall take the KY's as the new values of YO's and repeat the process. In order to prevent the machine from doing an endless number of iterations, a maximum number of cycles should be specified. Moreover, in order to prevent the value of YO's from becoming very large after h/ ^ many iterations, we may take the new YO's in terms of the maximum value of YO * Note that N is the number of divisions and "integer division" is assumed for evaluating the quantity N 2' -C70 -

Example Problem No. 81 or some value of YO which is close to the maximum, say the YO at a point near the midpoint (J = ) of the column. In other words, we shall define a new variable YO, such that YO(J) = ^ N for initial values (14a) YO( 7a YO(J) = KY(J) for subsequent iterations (14b) KY () In Eq. (14b) N KY(j) KY(J ) KY( ) = N (15) YO(7 2 The KY's are obtained by taking YOJ = W(J) (16) EW(J) (17) A and KY(J) = - where M' is the bending moment (18) caused by EW(J) From the definitions of YO and KY, we have yo(N) TO(N) / To(N) _ o(N) KY(O ) / (N KY( 2 KY( N K'(N)/ O() N KY (N9 Therefore, the formulas for CPCR and T may be rewritten as follows N2*YO (N) CPCR= (20) KY(2) T N2*YO(J) (21) (CPCR)*KY(N) -C71 -

Elastic Buckling Load for Columns of Non-Uniform Cross Section If the moment of inertia varies discontinuously at any of the equally spaced points, then one way of solving the problem is as follows: Let the number of points with a step change in moment of inertia be M. For M >0, we will have two sets of index numbers for the division points as shown below. M= 2 I i - I I I I i: I I j I I.Dumm;n.ex.S 0, / - I 1 / Dumry Ind4ex S o / 2{ 3 I I I I I I, ' 2 Indx No. Set ( O0 2 3 4 6? 8 9 /0 - 1-119. 9- -. - I - 3 srpo(s)o I 3 I I 7 I I 4 Yo Yo(o) Y O() YY 3) YoW4) Yoro) Yo() o(7)J Y YOS) YO) Y.I.. t t ' I _ 5 Index No. Set (Z) 2 3 4 5 6 7 8 9 /O // /2 6 STPP(S) -/ 8 1 I 83 7JI (o) 1(O ) 1(2) 1i)(4)) 11(6) 1(7) (I() 1(0) I 8 Y1 Y(o6) Y/(O/ YI(29 YtI) Y o()YI)Y ) Y/( Y/2) q Ini-tirl Yalu/S (OW YOve) Y) YbO) Y(q) Y)Y)Y (1/0)| 0 VaO/us of YJ YO(o) Yo(l) YOZ() ) YYt ) YOf(4) YO(5) YO(6) Y)71YiK8) Yo{q) YO(IO) after frst mo) IVa/es of mv Y() YaO) Yoi O(2) "YO() YO(4) YO(5) W)( Vt7 t7) Y Y(9) YO(IO) Iaf t e r — d, I v -__ _ _ _ _ _ _ _ _I_ I The index number set (1) will have one number for each division point while the index number set (2) will have two numbers at each step. The I's are stored according to the index number set (2). (See line 7 of the figure.) At each step in I, the index number in set (1) of the point with the step is stored in the vector* "STPO". (The two ends are also considered as steps.) For example, the index number 3 will be in the location for STPO(l). (See line 3 of the figure.) The values of YO(J) are first stored in the vector Yl(J) with J corresponding to the index number set (1). (See line 9 of the figure.) Then the step numbers, each corresponding to the index number in set (2) on the left side of a step, are computed and stored in the vector "STPP". For example, the index number 8 will be in the location STPP(2). (See line 6 of the figure.) Then starting from the right end, the value of Yl(J) is put into the location of Yl(J+l) until the index number of STPP(l) is reached. (See line 10 of the figure.) The process is then repeated until STPP(2) is reached and so on until STPP(M). (See line 11 of the figure.) * The "vector" is a linear array of numbers as defined in modern algebra and used in computer language.

Example Problem No. 81 Next, the concentrated elastic weights (EW) at the ends are computed. Then the elastic weights at two neighboring points at each step and the elastic weights between the steps are computed. Then the two concentrated elastic weights at the two neighboring points of each step are combined into one value and the values EW(J) are moved back to the locations with index number set (1). The procedure is similar to the one for Yl(J) described above but in reverse. With the elastic weights known, the procedure will be the same as that for columns of continuously varying moment of inertia. Table of Symbols In the following, the definitions of the variables which have been introduced previously and will be used in the program are recapitulated. Some new variables needed for the programming are also defined. N number of divisions M number of division points with step change in moment of inertia EPS an assigned tolerance value STPO(O)...STPO(M + 1) Index numbers of points with steps corresponding to index number Set (1). The two end points are also considered as steps. I(O)...I(M + N) The ratio of the actual moment of inertia (IA) about an axis normal to the plane of bending to the moment of inertia at a reference point (IO). The I's carry the index numbers of Set (1) for continuous variation and they carry the index number of Set (2) for discontinuous variation at the division points. YO(O)...YO(N) The ratio of the initial deflection at a point YO(J) to the same at a reference point near the center of the span YO(N). The YO's initially read in are assumed. For subsequent iterations, YO(J) = The YO's carry the index numbers of Set (1). K() STPP(O)...STPP(M + 1) Index numbers of the steps each corresponding to the index number in Set (21 on the left side of a step. Y1(O)...Yl(M + N) YO's carrying the index numbers of Set (2). EW(O)...EW(M + N) Concentrated elastic weights due to the distributed elastic load - Y(J) with division length A as I(J) a unit. These elastic weights carry the index numbers of Set (2). EWNN refers to the set of elastic weights EW(O)...EW(M + N). -C73 -

Elastic Buckling Load for Columns of Non-Uniform Cross Section EW(0)...EW(N) Similar to the weights in the set EWNN except that the weights carry index numbers of Set (1). The two weights at each stepping point of the Set EWNN have been combined into one. EW refers to the set of elastic weight EW(0)...EW(N). RL, RR respectively left and right reactions due to EW SW shear due to EW A2*Yo(N) KY(0)...KY(N) deflections with E*I as a unit. They are equal to the bending moments (with division length A as a unit) due to the concentrated elastic weights EW. N2*YO(N) CPCR = KY(2) 2 T N2*YO(J) (CPCR) *KY(N) CYL = number of cycles of iterations The data to be read in are as follows: N, M, EPS STPO(O)...STPO(M + 1) I(O)...I(M + N YO(O)...YO(N) (assumed) The data read in should be printed out immediately as an echo check. The output should be as follows: CYL the final number of iteration when convergence is obtained. CPCR the coefficient of the critical load. The critical load is understood as given by Eq. 13 PCR = (CPCR)* E*I L Where E is the modulus of elasticity, IO is the moment of inertia at a reference point and L is the length of the column as defined previously. In order to check the results by manual computations, the following intermediate steps shall be printed out for each iteration: Yl(0)...Yl(N + M) EWNN tW(O)...EW(N + M)] EW CEW(o...EW(N KY(O)...KY(N) -C74 -

Example Problem No. 81 Flow Diagram 0= O,. G.(+l =5TP(S) TPO(5) -, 1, $. S, I 152 \. _______, \ PRIN T — E______, rw~o - r 7*Y( 6 Y/ 7 VW — > ~V J I G + * a') 24 ~B J= N+S-L zI ) (O) ~-) -J \\ EW(NN I(NN+) ~N+N 6*YI-(N)-I) Y.rNN-2 L. 571 P(S) TP(N+M) L3W(NN)= -G+ IZ4 -- -N- Zi-(u f) Z (NN-Z)) PP/T / L \ UWUS-PS)Tu, + I (U+I) + |(U) F | 3r PP( —Ir I, m<r7oo v Y(2)]on / /-= — S r*pp(s) - ~EW(-=- (71; I(J)G;~ - I(J) EW(U) — =* WvU) + 6 /J53 pP(s) L =0, I,$.G. U +(J- /OY/(.) yI(J+I) (?/ Ul OPP(S+/) 'T. &. =-STPPTS+I)4 LI(J-)^ ^^ )70 i Tg)-/ PRINT ----- / r^ \ KYW = —L P (=, I,) 7S.G. -+ Eu - = U+z, I EW(-N) JE G.(NAW-$+I)/ - EW 0co) 0c,) zO. E w (0).. r 7 my!2,'o N + (oYOoNmI).YO.(N-'2)1, T G.N- E-W( EW(N)- N) +r xON - N-Z). Ew()-. ) 0 Yo) - +f., J. C7. F W(J. L7 L8ysw(YI — R RL= -99 T= 0, I, XC. N RL=RL- EW(1) 5W(o)= -TL = o, Z 5L8 I L 9 KY JS o.- PRIN T C, ='/? )tlm(6XWO -T 0, X, It YO') + s wC Ti)J-~j.N- YO (NKZ) pR/iNr ( <f)A4t c J CYL- START CPcCR -C75 -

Elastic Buckling Load for Columns of Non-Uniform Cross Section MAD Program KUANG-HAN CHU D028N 001 010 000 KUANG-HAN CHU D028N 001 010 000 SCOMPILE MAD, EXECUTE R ELASTIC BUCKLING LOAD FOR COLUMNS R OF NON-UNIFORM CROSS SECTION DIMENSION YO(50),I(100),EW(100) SW(50),KY(50), 1 STPO(50),STPP(100),Y1(100) INTEGER N,JCYLQMNN SU U1 STPO,STPP START READ FORMAT DATA1,NMEPS VECTOR VALUES DATA1=$2I4,F7.5*$ READ FORMAT DATA2,STPO(0')**STPO(M+1) VECTOR VALUES DATA2=$(18I4)*$ READ FORMAT DATA3,I(0).**I(N+M) VECTOR VALUES DATA3=$(8F9.3)*$ READ FORMAT DATA4,YO(0) **YO(N) VECTOR VALUES DATA4=$(8F9.3)*$ PRINT FORMAT ECHO1,N,M.EPS VECTOR VALUES ECHO1=$S1H1S993HN =,I4,S3~3HM =I3,5S4~ 1 4HEPSzF6.4*$ PRINT FORMAT ECH02,STPO(0)* *STPO(M+1) VECTOR VALUES ECHO2=$1HO,59,8HSTEPS AT/(S10,1215)*$ PRINT FORMAT ECH03,I(0)***I(N+M) VECTOR VALUES ECH03=$1HOQS9o10HRELATIVE I/(S10, 6F10.3)$* PRINT FORMAT ECH04,YO(0)*o.YO(N) VECTOR VALUES ECH04=$1H0S9,2HYO/(S10, 6F10.3)*$ CYLz1 Al WHENEVER M.E*0,TRANSFER TO A2 THROUGH L1,FOR J-=11,J.G*N L1 Y1(J)=YO(J) THROUGH L2#FOR S=01,[S.G.M+1 L2 STPP(S)=STPO(S)+S-1 THROUGH LS2,FOR S=1,lS.G.M THROUGH LS2,FOR J=N+S-1,-1,J.L.STPP(S) LS2 Y1(J+1)=Yl(J) PRINT FORMAT CHECK1,Y1(0).**Yl(N+M) VECTOR VALUES CHECK1=$1HO,59,2HY1/(S10. 6F10.3)*$ NN~N+M EW(O)=-(7*Y1(O)/I(O)+6*Yl(1l)/I(1)-Yl(2)/I(2))/24 EW(NN)=-(7*Y1(NN)/I(NN) 1 +6*Y1(NN-1)/I(NN-1)-Yl(NN-2)/I(NN-2))/24 STPP(M+1)=NN THROUGH L3~FOR S=1a1,S.G.M UzSTPP(S) EW<U+1)=-(7*YI(U+1)/I(U+l) 1 +6*Y1(U+2)/I(U+2)-Y1(U+3)/I(U+3))/24 L3 EW(U)=-(7*YI(U)/I(U) 1 +6*Y1(U-1)/I(U-1)-Y1(U-2)/I(U-2))/24 THROUGH LS3,FOR S=0,1,S.G.M U=STPP (S) U1lSTPP(S+1) THROUGH LS3,FOR J=U+2t1,J.G.(U-1)) LS3 EW(J)z-(Y1(J-1)/I(J-1)+10*Y1(J)/I(J)+Y1(J+1)/I(J+1))/12 PRINT FORMAT CHECK2,EW(O).**EW(NN) VECTOR VELUES CHECK2=$1HOS9,4HEWNN/(S10, 6F10.3)*$ THROUGH L4tFOR S=1lS.G.MM U=STPP(S)-S+1 EW(U)=EW(U)+EW(U+1) THROUGH L4,FOR J=U+2,19J.G.(NN-S+1) L4 EW(J-1)=EW(J) PRINT FORMAT CHECK3,EW(0).**EW(N) VECTOR VALUES CHECK3=$1HOS9,2HEW/(S10, 6F10.3)*$ TRANSFER TO A3 A2 EW(O)m-(7*YO(0)/I(0)+6*YO(1)/I (1)-YO(2)/I(2))/24 EW(N)=-(7*YO(N)/I(N)+6*YO(N-1)/I(N-1)-YO(N-2)/I(N-2))/24 THROUGH L5 FOR J1,1,J.G~.(N-1) L5 EW(J)m-(YO(J-1)/I(J-1)+lYO*O(J)/I(J)+YO(J+1)/I(J+l))/12 A3 RR=O THROUGH L6~FOR J*0,1JeGeN L6 RR~RR-EW(J)*J/N RL'-RR THROUGH L7,FOR Jx01[J.G.NN -C76 -

Example Problem No. 81 MAD Program, Continued L7 RL=RL-EW(J) SW(O)=RL THROUGH L89FOR J=01,J.G.(N-1) L8 SW(J+1)=SW(J)+EW(J) KY(O)=O THROUGH L9,FOR J=Oi,JJG. (N-1) L9 KY(J+1)=KY(J)+SW(J+1) PRINT FORMAT CHECK4,KY(0). *KY(N) VECTOR VALUES CHECK4=$1HOS9 2HKY/(S10, 6F10.3)*S Q=O CPCR = N*N*YO(N/2)/KY(N/2) THROUGH LlOFOR J=l11,J.G.(N-1) T=N*N*YO(J)/(KY(J)*CPCR) L10 WHENEVER (oABS (T-1) L.EPS)IAND (TGG*0 ),Q=Q+1 WHENEVER Q.E.(N-1) TRANSFER TO A4 OTHERWISE CYL=CYL+l END OF CONDITIONAL WHENEVER CYL.G.20 PRINT FORMAT NGD VECTOR VALUES NGD=$1HQOS9.10HYO NO GOOD*$ TRANSFER TO START OTHERWISE THROUGH LllFOR J=1,1~J.G.N L11 YO{J)IKY(J)/KYIN/2) TRANSFER TO Al END OF CONDITIONAL A4 PRINT FORMAT CYNCYL VECTOR VALUES CYN=$iHOS9.15HNO. OF CYCLES =tI4*$ PRINT FORMAT CRLtCPCR VECTOR VALUES CRL=$lHO S9,15HCOEFF. OF PCR =tF7.2*$ TRANSFER TO START END OF PROGRAM Data $DATA 10 2 *005 0 3 7 10 *2 2 2. 2 *2 1. 1. 1 1 1. *2 *2 *2 *2 0. *4 *6 *8 *9 1. *9.8 *6 *4 0. 10 0 *005 O 10 *25 *36.49 *64.81 1..81 *64 *49.36.25 0. *4 *6.8.9 1..9.8.6.4 0. 12 0 *005 0 12.2.2 *2.2 1. 1. 1. 1. 1 *.2 *2 *2,2 0..4.6.8.9.95 1..95.9.8 *6.4 0. -C77 -

Elastic Buckling Load for Columns of Non-Uniform Cross Section Computer Output N a 10 M a STEPS AT 0 3 RELATIVE I 0.200 1.000 0.200 2 EPS=0.0050 7 10 0.200 0.200 0.200 1.000 1.000 1.000 1.000 0.200 0.200 0.200 YO Y1 EWNI EW KY Yl m 00 0.000 0.600 06 0.,800 0.900 0.900 0.800 0.600 0.400 0.000 1.000 0.000 1*000 0.000 -0.375 -0.983 -0.375 0.400 0.600 0.800 0.800 0.900 0.900 0.800 0.800 0.600 0.400 -1.917 -3.000 -1.833 -0.417 -0.900 -0.900 -0.417 -1.833 -3.000 -1.917 -0.375 -1.917 -3.000 -2.250 -0.900 -0.900 -2.250 -3.000 -1.917 -0.375 0.000 8.558 15.200 18.842 20.233 20.233 18.842 15.200 8,558 0.000 -0.983 20.725 0.000 1.000 0.000 EWNN -0.363 -0.996 -0.363 0.413 0.733 0.909 0.909 0.976 0.976 0.909 0.909 0.733 0.413 -2.026 -3.607 -2.157 -0,468 -0.973 -0.973 -0.468 -2.157 -3.607 -2.026 EW KY Yl -0.363 -2.026 -3.607 -2.624 -0.973 -0.973 -2.624 -3.607 -2.026 -0.363 0.000 9.728 17.429 21.524 22.995 22.995 21.524 17.429 9.728 0.000 -0.996 23.493 0.000 1 000 0.000 EWNN -0.363 -0.996 -0.363 0.414 0.742 0.916 0.916 0.979 0.979 0.916 0.916 0.742 0.414 -2.034 -3.646 -2*177 -0.470 -0.975 -0.975 -0.470 -2.177 -3.646 -2.034 EW -0.363 -2.034 -3.646 -2.647 -0.975 -0.975 -2.647 -3.646 -2.034 -0.363 -0.996 KY 0.000 23 162 NO. OF CYCLES 9.801 21.689 a 3 17.568 21.689 23.162 17.568 9.801 0.000 23.661 COEFF. OF PCR a 4.23

Example Problem No. 81 Computer Output, Continued N 10 M 0 EPS=0 0050 STEPS AT 0 10 RELATIVE I 0.250 0.810 0.360 0.490 0.640 0.810 0.640 0,490 0.360 0.250 YO KY KY KY 0.000 0.400 0.600 0.800 0,900 0.900 0.800 0.600 0.400 0,000 0.000 5.104 9,180 12.039 13,662 13.662 12.039 9.180 5.104 0.000 0,000 5.257 9.570 12*587 14e289 14.289 12.587 9.570 5.257 0,000 0.000 5.243 9.555 12.574 14,277 14.277 12,574 9.555 5.243 0,000 1.000 1.000 14.171 14.804 14.793 NO. OF CYCLES = 3 COEFF. OF PCR = 6.76 N = 12 M 0 EPS=0.0050 STEPS AT 0 12 RELATIVE I 0,200 1.000 0,200 YO 0.000 1.000 0.000 KY 0,000 31.267 0,000 KY 0.000 32.735 0.000 KY 0,000 32.778 0,000 NO. OF CYCLES COEFF. OF PCR 0.200 0.200 0.200 1,000 1.000 1.000 1.000 0.200 0.200 0.200 0.400 0.600 0,800 0.900 0.950 0.950 0.900 0.800 0,600 0.400 11.183 20,450 26.717 29.325 30.771 30.771 29.325 26.717 20,450 11.183 11.604 21.445 28.056 30,756 32.236 32.236 30.756 28.056 21.445 11.604 11.610 21*470 28.095 30,799 32.280 32.280 30.799 28.095 21.470 11.610 3 = 4.39 -C79 -

Elastic Buckling Load for Columns of Non-Uniform Cross Section Discussion of Results Three example problems are shown. The first example is a column shown below which is the same as shown in the previous figure. With 10 division points and values of EPS and YO as shown in the program, the procedure converges in 3 cycles and gives the answer of (CPCR) = 4.23. The coefficient given by exact solution of the differential equation is 4.22 (see reference 2). IO.2 I I= 0.I ~I= O'.IA -I I= /O.Z I i I 0 / 3 4 5 6 7 8 9 / The second example is for a column as shown below. The coefficient determined by the numerical method is 6.76, as compared with the value 6.72 which is interpolated from a table given in reference 3. o *^ ~ - -~ t ^ 6 4 ^^_, 0.2 0.'6 o. 4 64 0 0 0/ -4f 049 0 The third example for the column shown below is the same as the first, except that it uses 12 divisions and any discontinuity of moment of inertia between two neighboring division points is ignored and the column is treated as if it were a column with moment of inertia varying continuously. The value of (CPCR) obtained is 4.39 as compared with the true value 4.22. L, = o. < { g sL - I r I 0.2 ~.* **= ( ~ ~ ~ 0. 1X.0.0 0 2 3 4 5 6 7 8 9 / / 12 It should be noted that although the examples shown are for columns with moment of inertia varying symmetrically with respect to the mid-point of the column, the program is perfectly general and is applicable to columns with moment of inertia varying unsymmetrically. -c80 -

Example Problem No. 81 Critique The program is a natural one for the computer since the computer is most suitable for any iterative procedure which requires us to do the same computation process over and over again. References 1. N. M. Newmark,"Numerical Procedure for Computing Deflections Moments and Buckling Loads, ASCE Transactions,1943. 2. S. Timoshenko, Theory of Elastic Stability. 3. S. Timoshenko & J. M. Lessels, Applied Elasticity.

Example Problem No. 82 A MAD PROGRAM FOR TRUSS ANALYSIS by K. H. Chu Civil Engineering Department Illinois Institute of Technology Course: Structural Theory Credit hours: 3 Level: Junior Statement of the Problem The problem is to write a program which will determine stresses in members of trusses where the method of joints is applicable. The problem is of repetitive type which is well adapted for computers. Also it is of fundamental importance because if it is used as a subroutine, then a series of problems in truss analysis can be solved easily; for example, influence lines for trusses, deflection of truss by virtual work, stresses in indeterminate trusses, etc. The method of joints is chosen instead of setting up a matrix and solving a system of simultaneous equations because there will be too many zero elements in the matrix. (The number of simultaneous equations the computer can solve by elimination is usually limited, say to 30.) Therefore, there is a need to develop a program which will not only conserve storage space but also will do the problem as we do it manually, i.e., the machine would select the joint to be solved and proceed to a joint which can be solved next. Solution Let us consider the truss shown on the next page. We will number the joints and members in any fashion and set up an arbitrary but convenient coordinate system. -C82 -

yI The data to be read in are as follows: N the total number of joints M the total number of members EPS1 a small number to guard against the case IDI| 0 (explained later) JT(1)...JT(N) index number of the joints JX 1...JX N X and Y coordinates of the joints JY 1... JY NI FXO(l)...FXO(N) the X and Y components of the force acting on each joint. Positive FYO(l)...FYO(N) if the component is in the positive direction of the coordinate axis, zero if there is none. UNMB(1)...UNMB(N) the number of members meeting at each joint. Their stresses are unknowns to be determined. For example, in the figure, the value of UNMB(3) will be equal to 5. MB(1...MB M JN( 1...JN(M JF 1..JF(M the index number of the members and the corresponding index number of the joints at the ends of each member. It does not matter which end is referred to as the near joint (JN) and which end is referred to as the far joint (JF). For example in the figure the value of MB(15) will be 15 and the values of JN(15i and JF(15) may be 5 and 7 respectively or 7 and 5 respectively. -C83 -

A MAD Program for Truss Analysis So far as reactions are concerned, we need to read in the following data: JRO The joint number where the reaction of unknown magnitude and direction is applied. For example in the figure the value of JRO is 5. R1 This value is equal to 1. The purpose of reading in this number is for extending the program to the case of finding reactions along 3 given lines of applications. The reactions will be identified as FR(Rl), FR(R2), FR(R3) with R1, R2, R3 being equal to 1, 2, 3 respectively. For the sample problem under consideration, it serves no useful purpose except to identify the reaction with given line of application as FR(R1) or FR(1). JR1 The joint number where the reaction of given line of application is applied. For example, in the figure the value of JR1 is 6. LXR1, LYR1, LR1 These are the numbers which specify the given line of application in the following way. The line of action is assumed to be a vector directed away from the joint and length of the vector (LR1) may be taken as unity or any other value. The X and Y components of this vector are denoted as LXR1 and LYR1 respectively. In order that we understand the correlation of the various data which have been read in, a concrete example is shown below. This is the same truss as shown earlier except that the dimensions and loadings are specified and the origin of the coordinate axes is at joint 1 and the X axis coincides with the lower chord. The data for JT, JX, JY, FXO, FYO and UNMB are given in Table 1, and the data for MB, JN and JF are given in Table 2. The data for the reactions are given in Table 3. In the tables, there are many items which have not been defined. They will be explained later., (s) ~(~) FR( FXRI O/WXX FYRO 4()30' Having read in all the data, the first step is to compute the following quantities: LX(l)...LX(M) where L's are the length of the members and LX and LY are LY(l)..LY(M) respectively the components of the vector of length L which is L(1)..L(M) directed from the joint (JN) toward the joint (JF). The computed values of LX, LY, and L for the members of the truss shown are given in Table 2. -c84 -

Example Problem No. 82 Next we shall compute the reactions with the definition of terms as follows: SFX, SFY, SMO FR(R1) FXR1, FYR1 are respectively the sum of the X components of the external forces, the sum of the Y components of the external forces and the sum of moments of the external forces about the joint (JRO). is the force of the reaction of known line of action at the joint (JR1). It will be positive if it is directed away from the joint. This quantity is denoted as a subscripted variable with the intention of writing a subroutine (or external function) for computing 3 reactions with given lines of action. The method of section is a possible application of this subroutine. are respectively the X and Y components of FR(R1). The signs of FXR1 and FYR1 are positive if they are in the positive directions of the coordinate axes. FXO(JR1), FYO(JR1) FXRO, FYRO FXO(JRO), FYO(JRO) The external forces acting on the joint (JR1) must be modified to include the reaction components just computed. are values at the joint (JRO) corresponding to those as defined for the joint (JR1). The answers for the values of FR(R1), FXR1, FYR1, FXRO, FYRO for the example shown are given in Table 3. The values of FXO, FYO at joints 5 and 6 after the reactions have been computed are as shown by the circled values in of Table 1. lines 7 and 8 Table 1 lJne Dummy Index I Z 3 4 5 6 7 8 9 / JT(I) /_ 2 3 4 5 6 7 8 9 2 IX(I) o. 30. 60. 90. /2 /20. 90. 60. 3 0. 3 JY() o. o.. o. o. 60. 40. 40. 40. 4 FXO() 0....0.. 0. 0000. 0. 5 FYO(I) -/0. -10. -/0.. 0. O. 0. 0o. O 6 UNMB1) 2 3 5 3 3 2 5 3 4 7 Af rSoAy/ng FXO 0. 0. 0. 0. '-5050.. 0 0. 0. 8 +hereactions FYO -/0. -/0. -/0. O4 0. ~ o. o. o. 9 After I UNMAB @. 5 3 3 5 3 (3 /0 solv/ng the FXO0, A750 0. 0 -5'C1. o0.. O. 0, 756 // | first Jointj FYr -/o. /0.00 -/o. - o/. t40. O. 0. 0.aO Table 3 Dqtq Answers JRo 5 FR(R/) 50o. TR / 6 FXR 50. ZXRI.00 FY I o. LYRI 0.00 FXRO -50. LRI. oo0 FYIqO 40. Table 2 LinelDummY Index J I 2 3 4 5 6 7 8 /O /I 1/2 3 4 /5 / M/B(J) 2_ 3 4. 5 6 7 /o // /Z./3 /4 /5 2 TJN() / 2 3 4 5 6 7 9 2 3 3 45 3 JF() 2 3 478 5 6 9 7 78 4 LX(J) 30. 30.. 30.3. 0. -30. -30. -30. -30. 0. -30. o. 30, 0. -30. 5 LLy(J) o. 0. 0. 60. -Zo. 0. 0. -40a +4040. +4-40. 40. 40. 40. 6 L(7) 30. 130. 30. 60. 06.o6 30. 30. 50. 40. 50. 40. 50. 40. 50. 7 j (J ) I / 8 Va1ves SeLX(J) 30. u 30. for the joint'. 9 Es.-forYi LY(J), ______0 -. 40. 10 After LY)9 M I_____ I// so/ving SF(J) 9999 __999 _Z the first FX C) 0.0 _ 7__ 50 __1 /3 jointf F-YJ) -7.5 f/0O0 /4 __ F (J) -750o ___ _ /2.5_ -C85 -

A MAD Program for Truss Analysis Next we shall compute the stresses in the members. We shall use the following symbols: JT(I) is the index number of the joint to be solved as identified by the dummy index I. MS(l), MS(2) are the index numbers of the two members meeting at JT(I). F(MS(l)), F(MS(2)) are the two forces in the members to be solved. They are positive if they are directed away from the joint JT(I). JTS(MS(1)), JTS(MS(2)) both store the value of JT(I) in order that in the output, we know the forces in the members MS(1) and MS(2) are the results of the solution of that particular joint. FX(MS(1)), FY(MS(l)) are the components of the forces acting on JT(I) FX(MS(2)), FYMS 2 in the members MS() and MS(2). They are positive if they are along the positive direction of the coordinate axes. IS(1), IS(2) are respectively the joint numbers of the other ends of the members MS(l) and MS(2) which meet at the joint JT(I). FXO(IS(1)), FYO(IS(1)) The external forces FXO, FYO acting on the joints FXO IS(2), FYO IS(2)) IS() and IS(2) shall be revised in order to include the components FX(MS(1)), FY(MS(1)), etc. as known forces. It should be noted that the signs of FX(MS(1)), FY(MS(1)), etc. should be reversed since they are originally considered as applied at the other ends of the members (i.e. at the joint JT(I)). The joint equations can be written very easily. In writing the joint equations, it should be noted that the joint under consideration JT(I) is taken as the near joint (JN). The forces are acting on the joint JT(I) and shall be considered as positive when they are directed away from the joint JT(I). We should also guard against the case where the denominator of the joint equation is nearly equal to zero. This is really the case when the joint is unstable. The remaining problem is how the machine can tell which joint is to be solved first and which to be solved next. The following suggestions may be helpful. Have the machine examine a list of UNMB. If there is a joint with UNMB=2, the index number I of the joint JT(I) is picked up. First look into the list for JN. If at the index number J, the value of JN(J) is equal to JT(I), then the index number J is picked up and stored in the location MS(l). Notice that this number J identifies the member MB(J) (See Table 1) which is renamed as MS(l), since this is the first member to be solved for at the joint under consideration. Then,in order to prevent the number in JN(J) and JF(J) from being picked up in future examination of the list, we change the values of JN(J) and JF(J) into some unreasonable number, say 9999. -C86 -

Example Problem No. 82 Then look further down the list of JN. If we find that at another value of J (say j'), JN(J') is equal to JT(I), then we store J' in MS(2) and the value of JF(J') of the number MB(J') in IS(2). Also we change the values of JN(J') and JF(J') into 9999. If we examine the list of JN and find either none of the values or only one of the values is equal to JT(I), then we must examine the list of JF to find either both of the members MS(l) and MS(2) or only the member MS(2). The procedure is similar to the one described above with the following exceptions. Suppose that we find the value of JT(I) in the location of JF(J). We should store the value of JN(J) in the location of IS and we should reverse the signs of the values LX(MB(J)) and LY(MB(J)). The sign change is necessary because the joint equation is written for the joint JT(I) which is considered as the near joint (JN). Since we find the value of JT(I) in the list of far joints (JF), the vector which represents the positive direction of the member MB(J) must be reversed. (Note that the value of JT(I) is identical with the index number I and the value of MB(J) is identical with the index number J.) Having found the two members meeting at the joint JT(I) ready to be solved, we end the search. We may set up a counter P to see whether P reaches the value of 2. Also we would like to update the number in the list of UNMB for the two members meeting at JT(I). The value of UNMiB(I) corresponding to the joint solved JT(I) should be reduced by 2. Each of the two values of UNMB corresponding to the far ends of the 2 members solved as identified by the joint numbers corresponding to IS(1) and IS(2) shall be reduced by 1. After going through the above procedure of searching and updating, we can go through the following steps: (1) Applying the joint equations to find the forces in the two members MS(l) and MS(2), (2) Obtaining the components of these forces, and (3) Modifying the forces applied at the joints at the far end of the member to include the components obtained as known forces. As a concrete example of the above procedure, we examine the list of UNMB in Table 1 and find that UNMB=2 at I=1. Thus JT(1) is the joint to be solved. Next we examine the list for JN in Table 2 and find that JN(J)=JT(I)=I=l at J=l. Thus we have MS(l)=l, noting that MS(l) really represents the member MB(1). We set IS(l)=JF(1)=2 and change the values of JN(1) and JF(1) to 9999. -C87 -

A MAD Program for Truss Analysis Looking further down the list of JN, we do not find any value equal to JT(I)=I=l. When we look into the list of JF, we find JF(J)=JT(I)=I=1 at J=9. Thus MS(2)=MB(9)=9. We set IS(2)=JN(9)=9 and change the values JN(9) and JF(9) to 9999 (See Table 2, lines 10 and 11). Since the number JT(l) is found in the list of JF, we change the signs of LX(9) and LY(9) (See lines 8 and 9 in Table 2). We will stop the search since we have already found the two members to be solved. We then reduce the value of UNMB(JT(l)) by 2 and the values of UNMB(IS(1)) and UNMB(IS(2)) each by 1 (See circled numbers in line 9 of Table 1). Since the joint solved is JT(l), the values of JTS(l) and JTS(9) for the members MB(1) and MB(9) respectively are equal to JT(l)=l. (See line 7 of Table 2). The forces F(MS(1)) and F(MS(2)) are then determined by applying the joint equations to the joint JT(1). (See line 14 of Table 2.) Their components with respect to the joint JT(1) are then determined. (See lines 12 and 13 of Table 2.) The forces FXO and FYO applied at joints IS(1) and IS(2) are then modified by adding the force components of MS(1) and MS(2) with sign reversed to the external forces acting at these joints. (See lines 10 and 11 in Table 1). After the first joint has been solved, we re-examine the list of UNMB (line 9 of Table 1) and repeat the whole procedure. We shall stop the process when all stresses in all the M members are solved. Since after each application of the joint equation, the total number of unknowns is reduced by 2, we may set up a counter, say Q, which will be increased by 2 after each cycle. Since M may be either even or odd, finally we may have (M-Q) equal to 0 or 1. If it is equal to zero, the stresses in all the members have been determined. If it is equal to 1, we shall find the remaining member using a procedure similar to the one described previously. The force in the remaining member can easily be determined. It should be noticed, however, that in order to maintain accuracy, the force should be determined from the component with absolute value larger than the other. We would like to print out the data as originally read in. They are: JT, JX, JY, FXO, FYO, UNMB, MB, JN, JF, JRO, R1, JR1, LXR1, LYR1 and LR1. The answers to be printed out are: FR(R1), FXR1, FYR1, FXRO, FYRO, FXO(JRO), FYO(JRO), FXO(JR1), FYO(JR1), LX, LY, L, JTS, FX, FY and F. -c88 -

Example Problem No. 82 With the detailed instructions as given above, the programming becomes more or less straight forward. From statement labeled "Start" to the statement just preceding the through loop L5, the data are read in and printed out in thespecific form as stated in the instructions. The values of LX, LY and L are computed in the loop L5. The fifteen statements following the statement SFX=O are essentially for determining the reactions. The following quantities are computed: SFX, SFY, SMO, FR(R1), FXR1, FYR1, FXO(JR1), FYO(JR1), FXRO, FYRO, FXO(JRO) and FYO(JRO). The formula for SMO may be written in the following form: SMO = {[JX(I)-JX(JRO] *FYO(I)- EJY(I)-JY(JRO)] *FXO(I) (1) With computer language, in a THROUGH loop for I, the equation for SMO becomes SMO = SMO + [JX(I)-JX(JRO)] *FYO(I)- [JY(I)-JY(JRO)] *FXO(I) (la) The formula for FR(R1) is written in the following form: FR(R1) = SMO*LR1 (2) [JY(JR1)-JY(JRO)] *LXR1- LJX(JR1)-JX(JRO)J *LYR1 The searching for the members to be solved and updating various items by the procedure as described in the instructions are performed from twenty-seven statements after and including the statement labeled C3. The joint equations are solved and the quantities F(MS(1)), FX(MS(1)), FY(MS(1)), FXO(IS(1)), FYO(IS(1)), F(MS(2)), FX(MS(2)), FY(MS(2)), FXO(IS(2)), and FYO(IS(2)) are computed from the statement (D =...) to the statement (FYO(IS(2)) =...). The joint equations are written in the following form: -FXO(JT(I)) * LY(MS(2) + FYO(JT(I)) * LX(MS(2)) L(MS(2)) L(MS(2)) (3) F(MS(1)) D LY(MStl)) LX(MStl/t FXO(JT(I)) * LMS - FO(JT(I)) * LMS 1( F(MS(2)) (4) D where LX(MS(1)) * LY(MS(2)) LY(MS(1)) * LX(MS(2)) (5) L(MS(1)) * L(MS(2)) L(MS(1)) * L(MS(2)) -c89 -

A MAD Program for Truss Analysis The counter Q registers the number of members already solved. If only one member is left to be solved, (M-Q) = 1, the joint with UNMB=1 is searched out by the loop LL4. The member to be solved is searched out by the loops LL5 and LL6. The force in the member is found by the statement starting from the statement labeled C6 to the statement END OF CONDITIONAL just before the statement labeled C4. The rest of the statements are for printing out the answers. Flow Diagram START READ N, M,EP51 JTO) * ~ ~ 3JTN) JX(O) ~~* JX(N) JY('). Y(N) FXO(1). FXO(N) FYO()-. FYO (N) UNMB(I).. UN(MB(N) MB (0).. * M6 (M) J.TN (1).... JN(M) JF(l) ~ ~. rF (M) Print N,, A/, EPS1 lRead T'RO 0Prin -. YR 0 L__ LX)= SYX(JF(7))- TJX(.TN()) (s _Rea d TR/,Pri ntL J/R o R LR LY(7)=YY(7F(J))- Y(XN(J)) - / RI R.TR, JLXRI,LYR LRI R, LYRY LR JRI RILYRI, LR.M L,(. =J LX(J)-x L*YX(J) t LY(J) LY(J),~~~~~~~Z 5SFX=u LNI SFX= S SFX+ f -XuI () ( SFY= 0 L -,LN SFY=SFY+ FYo(ir) ASMo=o See E. (/a)in textf or SMO 5ee Eq.(Z) in tex- for FR(RI) FXR I FA (R 1) * LXR I/LR / FYKI = FR(RI) *LYRI/LRI FXO ( JR) -= FXO(J-RI) -- FXRI F YO (R I) = FYO (IR/) J FYR I I GH FXRO = - (FX- FXRI) - FYR O = - (SFY - FYfI F) Pr'n t R I, FR(RI), FXF,/, Fi /, FX'Rc, FYRO FXo (JRo) = FAX (J7Ro) FXRO JRO, FXO(FZRo) cJRo, FYO( ), /JRI, FXO(JRI), TRI, FY (JRI) E FYO(JRO) = FYO (JRo) t FYR OQ, -C90 -

Example Problem No. 82 Flow Diagram, Continued =-O!LL1 LL. = QOI~. GN. 'NP48E P= LO =. IT U NMB (JT(I)) = See E'.(3) in fext for F(MS(i)) UN ) MB (Jrr))-Z ee2 Eq.(5) i I< EF E X(t))= FS(x ) 7-UNW - tr SnFY(Ms I))= F(MSW()*Ly(LfM)) /NVM (I5(0->/ -Pr/intf L(MS(i)) UNMB (IS(2)) I UNSfrA8LE rAr) FO( ) IUT rsM(Z) ())T) 5 fxo(L5()) - FX(.MS()) JTS(MS( = ))=(I) / =FYO(I ()= JTS(M())= Jr(l)FYO(IS())- FYS()) See Eq. (4) in textfor F(MS(Z)) * FX( 2 = F(45(2))* LX(MS(Z)) L(MS(Z)) FY(S(z2)) = F(Pf.52)) LY(MS z)) FXO(IS(z))= L(MS(z)) FXO (IS())- FX(M"(z)) rFO( Y(IS()= ) - FYo (ISz)) -Fy(s5(z)) MAD Program KUANG-HAN CHU D028N 001 003 000 KUANG-HAN CHU D028N 001 003 000 $COMPILE MAD, EXECUTE R A MAD PROGRAM FOR TRUSS ANALYSIS DIMENSION JT(50),FXO(50),FYO(50),UNMB(50)tJX(50) JY(50) DIMENSION MB(100) JN(100) JF(100),LX(lOO) LY(100) L(100) DIMENSION JTS(lOO) FX( 100)F Y( 100),F(lUO IS(2),MS(2) FR(3) INTEGER JT,UNMB,JNJFMBJTS,IS,MSPQ,I,J,M,NJRO,JR1,R1 START READ FORMAT DATA1,N,M,EPS1 VECTOR VALUES DATA1i$2I4,F8.5*$ READ FORMAT DATA2,JT(1)*..JT(N) VECTOR VALUES DATA2z$(184)*$ READ FORMAT DATA3,JX(1)...JX(N) VECTOR VALUES DATA3=$(8F9.3)*$ READ FORMAT DATA4,JY(1)...JY(N) VECTOR VALUES DATA4=$(8F9.3)*$ READ FORMAT DATA5,FXO(1)..*FXO(N) VECTOR VALUES DATA5=$(8F9.2)*$ -C91 -

A MAD Program for Truss Analysis MAD Program, Continued READ FORMAT DATA6,FYO( 1).FYO(N) VECTOR VALUES DATA6S$(8F9-2)*$ READ FORMAT DATA7tUNMB(1)***UNMB(N) VECTOR VALUES DATA7-$(1814)*$ READ FORMAT DATAllMB(l) 1.MB(M) VECTOR VALUES DATA1l=$(1814)*$ READ FORMAT DATA12,JN(1),**JN(M) VECTOR VALUES DATA12=$(1814)*S READ FORMAT DATA13tJF(1),**JF(M) VECTOR VALUES DATA13=$(1814)*$ PRINT FORMAT ECHOlN#MtEPSl VECTOR VALUES ECHOls=$H1S7t16HNO. OF JOINTS, I2.54, 1 17HNO. OF MEMBERS =- I3vS57HEPS1 =,F8.5/*S WHENEVER N.G.5 THROUGH L1,FOR I5=1 5*I.G*N WHENEVER (I+4)*GE.N TRANSFER TO Al OTHERWISE PRINT FORMAT ECH02 JT(I)**.*JT(1+4) VECTOR VALUES ECH02=$lHOS7,2HJTtS2, 5I10/*$ PRINT FORMAT ECH03#JX(I)..JX(I+4) VECTOR VALUES ECH03=$S892HJXS5 5F10.3/*$ PRINT FORMAT ECH04,JY(I)**.JY(I+4) VECTOR VALUES ECH04=$S8,2HJYS5 5F10.3/*$ PRINT FORMAT ECH05,FXO(I).**FXO(I+4) VECTOR VALUES ECH05=$S8,3HFXO,549 5F10*2/*$ PRINT FORMAT ECH06#FYO(I)***FYO(I+4) VECTOR VALUES ECH06=$S8,3HFYOS4. 5F10.2/*$ PRINT FORMAT ECH07,UNMB(I)**.UNMB(1+4) VECTOR VALUES ECH07=$S8,4HUNMB, 51l0/*$ L1 END OF CONDITIONAL OTHERWISE I1l END OF CONDITIONAL Al PRINT FORMAT EECH02. JT(I)***JT(N) VECTOR VALUES EECH02=$S892HJTS2, 5I10/*$ PRINT FORMAT EECH03,JX(I)*..JX(N) VECTOR VALUES EECH03=$S8.2HJX,55 5F103/*$ PRINT FORMAT EECH04,JY(I)**.JY(N) VECTOR VALUES EECH04=$S8,2HJY55, 5F10*3/*$ PRINT FORMAT EECH05,FXO(I)**.FXO(N) VECTOR VALUES EECH05O$S8,3HFXOS4, 5F10.2/*$ PRINT FORMAT EECH06,FYO(I)e**FYO(N) VECTOR VALUES EECHO06-$S83HFYOS4, 5F10.2/*$ PRINT FORMAT EECH07,UNMB(I)**.UNMB(N) VECTOR VALUES EECH07=$S8t4HUNMB, 5I10/*$ WHENEVER M.G.5 THROUGH L2,FOR JIl, 5,J.G.M WHENEVER (J+4)GE.M TRANSFER TO A2 OTHERWISE PRINT FORMAT ECHOllMB(J)...MB(J+4) VECTOR VALUES ECHOll$lHOS722HMBS2, 5I10/*$ PRINT FORMAT ECH012,JN(J).*.JN(J+4) VECTOR VALUES ECH012=$S8,2HJN,52, 5I10/*$ PRINT FORMAT ECH013,JF(J)**.JF(J+4) VECTOR VALUES ECH013=$S8,2HJFtS2. 5I10/*$ L2 END OF CONDITIONAL OTHERWISE J1l END OF CONDITIONAL A2 PRINT FORMAT EECOlIMB(J)***MB(M) VECTOR VALUES EECOll=$lHO,57,2HMBS2o 5I10/*$ PRINT FORMAT EEC012,JN(J)*..JN(M) VECTOR VALUES EEC012=$S8,2HJNS2, 5I10/*$ PRINT FORMAT EEC013#JF(J)**.JF(M) VECTOR VALUES EEC013=$S8,2HJFS2, 5I10/*$ READ FORMAT ROoJRO VECTOR VALUES RO$I4*$ READ FORMAT RARlJR1,LXR1~LYR1,LR1 VECTOR VALUES RA=$I4,S29I4,3F10*.5* PRINT FORMAT RROJRO VECTOR VALUES RRO=$S8,4HJRO,I3/*$ -C92 -

Example Problem No. 82 MAD Program, Continued PRINT FORMAT RR1,JR1,LXR1,LYR1,LR1 VECTOR VALUES RR1=$S8,4HJR1 I 3/S8 4HLXR1 S3 F10.3 1 S5,5HLYR1,F10.3#S5,5HLR1,F10.3/*$ THROUGH L3,FOR J=-lltJ.G*M LX(J)=JX(JF(J))-JX(JN(J)) LY(J)=JY(JF(J))-JY(JN(J)) L3 L(J)=SQRT.(LX(J)*LX(J)+LY(J)*LY(J)) SFX=O SFY=O SMO=O THROUGH LN1,FOR I=11PI.G..N SFX=SFX+FXO(I) SFY=SFY+FYO(I) LN1 SMO=SMO+(JX(I)-JX(JRO))*FYO(I)-(Jy(I)-JY(JRO))*FXO(I) FR(R1)=SMO*LR1/{(JY(JR1)-JY(JRO))*LXR1 1 -(JX(JR1)-JX(JRO))*LYR1) FXR1=FR(R1)*LXR1/LR1 FYR1=FR(Rl)*LYR1/LR1 FXO(JR1)=FXO(JR1)+FXR1 FYO(JR1)=FYO(JR1)+FYR1 FXRO=-(SFX+FXR1) FYRO=-(SFY+FYR1) FXO(JRO)=FXO(JRO)+FXRO FYO(JRO)=FYO(JRO)+FYRO PRINT FORMAT ANS4R1i>FR(R1)tFXRltFYR1.FXROFYRO VECTOR VALUES ANS4=$1HOS7t3HFR(9I291H),S1lF1O.2,S5S4HFXR1, 1 S1tF102,S5F4HFYR1S1,lF1Oe2/S8,4HFXROS3tF1lO2,55$4HFYRO~S1l 2F10.2/*$ PRINT FORMAT ANS5,JROFXO(JRO),JROFYO(JRO), 1 JR1,FXO(JR1) JR1,FYO(JR1) VECTOR VALUES ANS5=$S8,4HFXO(,I2t1H),F10.2,S3,4HFYOC(I291H)t 1F10.2/S8,4HFXO(.I21H) F1O.2,S3~4HFYO(,12,1H) F10.2*$ Q=O C3 THROUGH LL1,FORI=,1 1I.GGN LL1 WHENEVER UNMB(I)*E.2,TRANSFER TO C1 Cl P=0 THROUGH LL2,FOR J=11J.*G.M WHENEVER JN(J)*E.JT(I) P=P+1 IS(P)=JF(J) MS(P)=J JN(J)=999 JF(J)=9999 WHENEVER P.E.2.TRANSFER TO C2 LL2 END OF CONDITIONAL THROUGH LL3OFOR J=11t~J.G.M WHENEVER JF(J).E.JT(I) P=P+1 IS(P)=JN(J) MS(P)=J JN(J)=9999 JF(J)=9999 LX(J)=-LX(J) LY(J)=-LY(J) WHENEVER P.E.2tTRANSFER TO C2 LL3 END OF CONDITIONAL C2 UNMB(JT(I))=UNMB(JT(I))-2 UNMB(IS() )=UNMB(IS(1))-1 UNMB(IS(2))=UNMB(IS(2))-1 JTS(MS(1) )JT(I) JTS(MS(2))=JT(I) D=LX(MS(1))*LY(MS(2))/(L(MS{(1)*L(MS(2))) 1 -LY(MS(1))*LX(MS(2))/(L(MS(l))*L(MS(2))) WHENEVER.ABS.(D).L.EPS1 PRINT FORMAT UNSTBL VECTOR VALUES UNSTBL=$S8,22HTHIS TRUSS IS UNSTABLE*S TRANSFER TO START END OF CONDITIONAL -C93 -

A MAD Program for Truss Analysis MAD Program, Continued F(MS(1))=(-FXO(JT(I)L*LY(MS(2))/L(MS(2)) 1 +FYO(T(I))*LX(MS(2) )/L(MS(2)))/D FX(MS(1))=F(MS(1))*LX(MS(1))/L(MS(1)) FY(MS(1))=F(MS(1))*LY(MS(1))/L(MS( 1)) FXO(IS(1 )=FXO(IS(1))-FX(MS(1)) FYO(IS(1) )FYO(IS1) )-FY(MS(1)) F(MS(2))-(FXO(JT(I))*LY(MS(l))/L(MS(l)) 1 -FYO(JT(I))*LX(MS(1))/L(MS(1)))/D FX(MS(2) )F(MS(2))*LX(MS(2))(2/L(MS)) FYMS(2)) =F(MS(2))*LY(MS(2) /L(MS(2)) FXO(IS(2) )FXO(IS(2))-FX(MS(2)) FYO(IS(2) )FYO(IS(2))-FY(MS(2)) Q=Q+2 WHENEVER (M-0).G61,TRANSFER TO C3 WHENEVER (M-Q).E.OTRANSFER TO C4 THROUGH LL4,FOR IltlI.G.N LL4 WHENEVER UNMB(I)*E.1,TRANSFER TO C5 C5 THROUGH LL5,FOR J=1,1,J.G.M LL5 WHENEVER JN(J).E.JT(I)JTRANSFER TO C6 THROUGH LL6,FOR J=1,1,J.G.M WHENEVER JF(J).*EJT(I) LX(J)=-LX(J) LY{J)=-LY(J) TRANSFER TO C6 LL6 END OF CONDITIONAL C6 FX(J)=-FXO(I) FY(J)=-FYO( I) JTS(J)=JT(I) WHENEVER *ABS(F Y(J))*GE.*ABS^(FX(J)) F(J)=FY(J)*L(J)/LY(J) OTHERWISE F(J)=FX(J)*L(J)/LX(J) END OF CONDITIONAL C4 WHENEVER M.G*5 THROUGH L4OFOR J=l1 5tJ*.GM WHENEVER (J+4) GE.M TRANSFER TO A3 OTHERWISE PRINT FOPMAT ECHO11,MB(J).**MB(J+4) PRINT FORMAT ANS1,LX(J)*..LX(J+4) VECTOR VALUES ANS1=$S892HLX#S5, 5F10.3/*$ PRINT FORMAT ANS2,LY(J)**.LY(J+4) VECTOR VALUES ANS2=$S8,2HLY#55t 5F10.3/*$ PRINT FORMAT ANS3.L(J)...L(J+4) VECTOR VALUES ANS3=$S8,1HL,56, 5F10.3/*$ PRINT FORMAT ANS11,JTS(J)***JTS{J+4) VECTOR VALUES ANS11=$S8,3HJTSS1s 5110/*$ PRINT FORMAT ANS129FX(J).**FX(J+4) VECTOR VALUES ANS12=$S8,2HFXS.5, 5F10.2/*$ PRINT FORMAT ANS13,FY(J).**FY(J+4) VECTOR VALUES ANS13=$S8,2HFYS55S 5F10.2/*$ PRINT FORMAT ANS14tF(J)**.F(J+4) VECTOR VALUES ANS14=$S8,1HF9S6, 5F10.2/*$ L4 END OF CONDITIONAL OTHERWISE Jul END OF CONDITIONAL A3 PRINT FORMAT EECO11lMB(J)***MB(M) PRINT FORMAT AANS1*LX(J)*.,LX(M) VECTOR VALUES AANS1=$S8,2HLX~S5 5F10.3/*$ PRINT FORMAT AANS2,LY(J)*..LY(M) VECTOR VALUES AANS2=$S8,2HLY,55, 5F10.3/*$ PRINT FORMAT AANS3tL(J)*..L(M) VECTOR VALUES AANS3=$58,1HLS6, 5F10.3*$ PRINT FORMAT AANS119JTS(J).**JTS(M) VECTOR VALUES AANS11=$S893HJTSS1* 5I10/*$ PRINT FORMAT AANS12,FX(J).**FX(M) VECTOR VALUES AANS12=$S8,2HFX9S5t 5F10.2/*$ PRINT FORMAT AANS13,FY(J)...FY(M) VECTOR VALUES AANS13S$S8,2HFYS59 5F10.2/*$ PRINT FORMAT AANS14$F(J).**F(M) VECTOR VALUES AANS14$S8,1HFS6, 5F10.2*$ TRANSFER TO START END OF PROGRAM -C94 -

Example Problem No. 82 Data $DATA 9 1 2 1 1 2 5 1 15 2 0O 30. 0. O. 40. 0O O. -10, 3 3 2 2 3 3 4001 4 5 30. 0. 0O -10. 3 3 3 4 5 4 5 5 6 6 7 8 9 60. 90* 120. 120* 90. 60. 40. 40. 0. 0. 0. 60. 0. 0. 00. 0. 0. 5 3 3 4 -10. -10. 2 5 3 4 6 7 8 9 6 7 8 9 7 8 9 1 0.0 O. 0. 0O 0O 10 2 9 1.0 11 3 9 12 3 8 13 3 7 14 4 7 15 5 7 6 1.0 Computer Output NO. OF JOINTS = 9 JT 1 JX 0*000 JY 0*000 FXO 0.00 FYO -10.00 UNMB 2 JT 6 JX 120.000 JY 60.000 FXO 0.00 FYO 0.00 UNMB 2 MB 1 JN 1 JF 2 MB 6 JN 6 JF 7 MB 11 JN 3 JF 9 JRO 5 JR1 6 LXR1 1.000 FR( 1) 50.00 FXRO -50.00 FXO( 5) -50.00 FX3( 6) 50.00 MB 1 LX 30.000 LY 0.000 L 30.000 JTS 1 FX -7.50 FY -0.00 F -7.50 MS 6 LX -30.000 LY -20.000 L 36.056 JTS 6 FX -50.00 FY -33.33 F 60.09 MB 11 LX -30.000 LY 40.000 L 50.000 JTS 3 FX 15.00 FY -20.00 F -25.00 NO. OF 2 30.000 0.000 0.00 -10.00 3 7 90.000 40.000 0.00 0.00 5 2 2 3 7 7 8 12 3 8 LYR1 FXR1 FYRO FYO( 5) FYO( 6) 2 30.000 0.000 30.000 2 -7.50 -0.00 -7.50 7 -30.000 0.000 30.000 7 -22.50 0.00 22.50 12 0.000 40.000 40.000 3 0.00 0.00 0.00 MEMBERS = 3 60.000 0.000 0.00 -10.00 5 8 60.000 40.000 0.00 0.00 3 3 3 4 8 8 9 13 3 7 0.000 50.00 40.00 40.00 0.00 3 -30.000 -0.000 30 000 4 45.00 0*00 -45.00 8 -30.000 0.000 30.000 8 -22.50 0.00 22,50 13 -30.000 -40.000 50.000 7 -22.50 -30.00 37.50 15 EPS1 = 0*00100 4 5 90.000 120.000 0.000 0.000 0.00 0.00 -10.00 0.00 3 3 9 30.000 40.000 0.00 0.00 4 4 5 4 5 5 6 9 10 9 2 1 9 14 15 4 5 7 7 LR1 1.000 FYR1 0,00 4 5 -30.000 -0.000 -0.000 -60.000 30.000 60.000 5 6 45.00 0e00 0.00 33.33 -45.00 -33.33 9 10 30.000 0.000 40.000 40.000 50.000 40,000 1 2 7.50 0.00 10.00 10.00 12.50 10.00 14 15 0.000 -30.000 40.000 40.000 40.000 50.000 4 5 0.00 5.00 10o00 -6.67 10.00 -8.33 -C95 -

A MAD Program for Truss Analysis Critique In the writer's opinion, the program is well suited for instructive purposes. The problem can be separated into three parts. In the first part, the students can be asked to write a program (or draw a flow diagram) for computing the lengths and components of the length vector of the members. In the second part, the students can be asked to write a program for computing the reactions and modifying the external forces at the joints where the reactions are applied. There are also possibilities of writing the program for computing length or reactions or both as an external function.* The students may be asked to write a program for computing three reactions along three given lines of action and applied at three different joints of the truss.* This program, if it is written as an external function, may be applicable to trusses which can be solved only by method of sections.* Then the students can be asked to write a program for computing two unknown forces in members meeting at a joint using the method of joints and modifying the external forces applying at the joint at the other ends of the members. The process of searching for the joint to be solved and updating the various items involved is a rather complicated procedure. Detailed instructions as given in the text seem to be necessary. However, this complicated procedure is very instructive since it dramatically illustrates the computer's ability to handle problems which are essentially logical rather than computational. Often a human being is scarcely aware of the logical detail which is involved where he solves a problem manually. * Due to lack of time, these problems have not been worked out by the writer. -C96 -

Example Problem No. 83 STIFFNESS FACTORS FOR A FLAT SLAB COLUMN by Robert B. Harris Department of Civil EngnQring The University of Michigan Course: Advanced Design of Structures Credit hours: 3 Level: Graduate Statement of Problem Estimate the slab thickness, drop panel thickness, and the column diameter for a typical panel in your flat slab structure. Use the ACI methods for this approximation. Using these results along with the story height, study the effect of column capital size upon the stiffness coefficients of the column and the carry-over factors. Program your solution for the digital computer using the MAD language. Solution Stiffness coefficients for members with variable cross sections may be found by expressing the relationship between unit end moments applied to a member and the rotations of the end tangents. If the ends are defined by a and b respectively, and any end moments Mab and Mba are applied simultaneously, the rotation of tjie end tangents a and Gb must be,according to reference 1: Ga = Mabal - Mbaa2 b -Mabbl + MbaSb2 Dal and Pbl are the angle changes due to a unit moment applied to end a and Pa2 and,b2 are the angle changes due to a unit moment applied at end b. By the reciprocal theorem,angle Pbl is equal to Pa2' Solving the above equations for Mab and Mba: b2 + bl Nab k= a -+ b Mba -. a A b -C97 -

Stiffness Factors for a Flat Slab Column in which -A.= alb2 - (bl)2 EI If I is the minimum moment of inertia then K = - and: 11 Pal 1 KA C bl c2 - K~-; c b2 3 = KA_ The following sketch defines lengths used. The following sketch shows the conjugate beams loaded with the appropriate diagrams from which the values of al bl and b2 EI diagrams from which the values of p B P an< Bb? can 1be 2 oun< B -C98 -

Example Problem No. 83 Assuming circular elements in the capital, the moment of inertia varies as the fourth power of the diameter, therefore I(H-Q) ^ D and (H-T) C (H-Q), D ) (DY)4 I(y) (Dy)4 In order to control the accuracy of the solution a variable was introduced called DEL which represents the size of increment AY. This was assigned the value of 0.2 in the sample program. This was particularly necessary since rectilinear integration was used. The following summation equations are the result of the application of the conjugate beam method and the letters appearing in brackets below the various terms are the letters used in the program for that quantity. n = H-Q iEI6A1 = n = A 2 n = H-T 2 Yn Y AY + D4 I -- bAY n =/ l (Q (DY) n = (H-Q A ) -" - Al A 2EIfBi = (H - Yr n^n = 2 B n = H-T 4 -- ) Yn AY + D T n = (H-Q4S ) 2 (H - Yn) Y 4(DY (DY) I I I L Bl I I G G I K n = H-Q H2IB2 =, b y aY n "-2 n = H-T (H - Yn)2 AY + D4 n = (H-Q+y-) 2 (H - Y ) (DY) I I L B2 I L F E -099 -

Stiffness Factors for a Flat Slab Column Letting X = H2EI, then A1 = Al DB1 Bl PB2 = B2 PB2 = i 1 r A =PAlB2 - Pl2 A = 5A1lB2 - 5B1 (Al) (B2) X2 - B1l' where U = (A1)(B2) - B12 U X2 X Al K U (Al)X2 XKU Coefficient at bottom, C1 = X HEI = H K EI H Hence, Coefficient at top, H3Al 1 U H3B2 C3 U For the carry-over factors the coefficient C2 is needed. C2 Carry over - top to bottom Carry over - bottom to top HBl1 U B1 B2 B1 A1 Flow Diagram (Start) —* Print IRead - HC = 0 /Through Title Data BACK D,T,H For C D,1.,HC.E.H/2 Calculate A = 0 /Through ) HC, Q E =E = 0 Y SUMA G = 0 For Y DEL = 0.2.l,01YG.H -C100 -

Example Problem No. 83 MAD Program, and Data R. B. HARRIS S164E 003 020 000 R. B. HARRIS S164E 003 020 000 $COMPILE MAD, EXECUTE R R PROGRAM FOR STIFFNESS COEFFICIENTS FOR A COLUMN WITH 90 DE 1GREE CAPITALS R START PRINT FORMAT TITLE R R HEAD INGS R VECTOR VALUES TITLE=$1H1,S23 26HCOLUMN STIFFNESS STUDY 1 //S26,19H90 DEGREE CAPITAL ///S5,3HCAPS3,3HCOL~S4,3HCAP, 2S3,4HSLABS5,3HCOLS4t3HTOPS6,3HBOTS5,5HCARRY*S4t5HCARRY/ 3S5 4HOIAM,S2 4HDIAM S2,5HDEPTHS2,5HTHICKS2 6HLENGTHS2, 45HCOEFF,S4,5HCOEFFS3,7HTOP-7BOT S2,7H0T-TOP*$ R R INPUT R READ FORMAT DATA D, T. H VECTOR VALUES DATA = $3F10.3*$ HC = 0 THROUGH BACK. FOR C a Dl HCHE.H/2 HC = C/2 - D/2 Q - T + HC R R INTEGRATING STRAIGHT LINE PORTION R A = 0.0 E - 0.0 G = 0.0 DEL = 0.2 THROUGH SUMA, FOR Y = *1~*2,Y*G.(H-Q) A = A +(Y.P.2)*DEL E = E +((H-Y).P.2)*DEL SUMA G = G +((H-Y)*Y)*DEL R R INTEGRATING CURVED PORTION R DY = 0.0 B 0.0 F a 0.0 K = 0*0 DEL = 0.1 THROUGH SUMB~ FOR Y = (H-Q+.05).S1,Y*G.(H-T) DY = DY + D + 2*(C-D)*(Q-H+Y)/HC B = B +(D.P.4)*DEL*Y*P.2/DY*P.4 F = F +(D.P.4)*DEL*(H-Y).P.2/DY.Ps4 SUMB K = K +(D.P.4)*DEL*(H-Y)*Y/DY#P.4 R R STIFFNESS FACTORS AND OUTPUT R Al = A+B B1 = G+K B2 = E+F U = A1*B2-B1.P*2 CBOT = A1*H*P.3/U CTOP = B2*H*P.3/U COBT = B1/A1 COTB = B1/B2 BACK PRINT FORMAT COEFFt C~ Dt HC. THt, CTOPCBOT.COTB.COBT VECTOR VALUES COEFF = $F8*1,4F7,1,4F9F4*$ TRANSFER TO START END OF PROGRAM SDATA 1*0 0,000 10.0 1.0 1,000 10.0 -C101 -

Stiffness Factors for a Flat Slab Column Computer Output COLUMN STIFFNESS STUDY ______..__. _.___ 90. DEGREE CAPITAL CAP COL CAP SLAB COL TOP OT CARRY CARRY DIAM DIAM DEPTH THICK LENGTH COEFF COEFF TOP-BOT BOT-TOP 1.0. ~~1.0.U.. 160.0 4.0012 4.0012.5002.5002 2.0 1.0..5.0 10.0 4.6159 4.1473_.4985 _.5548 3.0 1.0 1.0.0 10.0 5.9587 4.4210.4878.6575 4.0 1.0 1.5.0 10.0 7.1344 4.6274.4765.7346 5.0 1.0 2.0.0 10.0 9.4565 4.9705.4538.8634 6.0 1.0 2.5.0 10.0 11.5429 5.2332.4359.9616 7.0 1.0 3.0.0 10.0 15. 1811 5.6761.4056 1.1278 8.0 1.0 3.5.0 10.0 19.7202 6.0218.3838 1.2568 9.0 1.0 4.0.0 10.0 28.0289 6.6153.3492 1.4796 10.0 1.0 4.5.0 10.0 36.1126 7.0907.3253 1.6568 11.0 1.0 5.0.0 10.0 54.0806 7.9271.2889 1.970 -CULUMN STIFFNS S STUDY 90 DEGREE CAPITAL CAP COL CAP DIAM1 DIAM I-EP r 1.0 1.0.0 2.0 1.0.5 3.0 1.0 1.0 4.0 1.0 1.5 5.0 1.0 2.0 6.0 1.0 2.5 7.0 1.0 3.0 8.0 1.0 3.5 9.O 1.0 4.0 10.0 1.0 4.5 11.0 1.0 5.0 SLAB COLL TOP BOT CARRY CARRY THICK LENGTH COEFF CUEFF TOP-BUT bOT-TOP 1.0 10.0 6.C930 4.4461.4866.6669 1.0 10.0 7.1344 4.6274._ 4765.7346 1.0 lu.0 9.4565 4.9705.4538.8634 1.0 10.0 11.5429 5.2332.4359.9616 1.0 10.0 15.7811 5.6761.4056 1.1278 1.0 10.0 19.7202 6.0218.3838 1.2568 1.0 10.0 28.0289 6.6153.3492 1.4796 1.0 10.0 36.1126 7.0907.3253 1.6568 1.0 10.0 54.0806 7.9271.2889 1.9707 1.0 10.0 72.7273 8.6220.2643 2.2293 1.0 10.0 117.4918 9.8889.2276 2.7045 Discussion of Results Capital diameter, column diameter, capital depth, slab thickness, and column length must all be put in using the same units. In the sample these are all in feet and are printed out with the results for checking. The program limits the capital depth to one-half the column length and increments by one-half a length unit. Coefficients obtained give an accuracy of four significant figures. This accuracy may be improved by making the integration increments smaller in the SUMA and SUMB iterations. -C102 -

Example Problem No. 83 Critique This problem was selected as representative of the use of the computer in making repetitive computations which would take considerable time to do longhand. In a sense, the problem is synthetic. If, however, such information were available with the detail obtained, it would materially help in preliminary design and checking of flat slab construction. The problem was designed to serve as an introduction to the use of a more complicated program prepared by the instructor, the use of which would allow the student to make adjustments in columns, slabs, drop panels, and capitals, and exercise the broadest judgement in design. As an introduction, the student reaction was quite good and little difficulty was experienced in programming. About one week and a half (3 laboratory periods) were used in this activity, but actual results were not turned in for some time after this since mechanical delays with the Computing Center prolonged the submission date. This problem served the purpose for which it was intended and in addition firmed up in the students' minds the effectiveness of the column capital in flat slab construction. Reference 1. Maugh, L. C. Statically Indeterminate Structures, page 133. -C103 -

Example Problem No. 84 VIBRATION OF BEAMS ON SPRING SUPPORTS by Shan S. Kuo Department of Civil Engineering Yale University Course: Structural Dynamics Credit Hours: 2 Level: Senior, Graduate Statement of the Problem Write a computer program to determine the natural shapes for a beam on spring supports as shown in Figure 5 4 3 2 k frequencies and the modal 1. The beam is idealized 1 kII L h h h h X = 4h <,~~~~ Figure 1. Beam on spring supports as a multi-mass system. Each of the two end masses is taken as one half of an interior mass, and the sum of all masses is equal to the total mass of beam m. In order to study the effect of spring constant k, we take c R/ M H as a parameter where R = 2k/m = "rigid body"frequency 6H =7 2 EI/m = fundamental frequency of simply supported beam The program should be so written that any number of lumped masses can be handled. Solution We are concerned with the following matrix equation: X = FMX 60 2 where F = flexibility matrix due to the spring action and the beam action M = mass matrix, a diagonal matrix -C104 -

X = displacement of lumped mass O = frequency It should be noted that the matrix FM is not symmetrical; therefore, Jacobi's method cannot be applied to this case directly. The following is a general procedure for solving this problem using a digital computer. (a) Formulate the flexibility matrix F. It is considered here as the sum of three action matrices, one from the beam action and the other two from the spring action. In connection with the beam action matrix, we take kh T 4 R 2 EI 2(n+l)3 RH where n = number of interior masses, or number of total masses -2. (b) Formulate the mass matrix M. (c) Formulate the matrix M1/2 FM1/2 (d) Make use of the Jacobi's method to find the eigenvalues and the eigenvectors of the matrix found in the step c. (e) The natural frequencies are then equal to (eigenvalues)1/2 (f) The true eigenvectors are equal to M1/2 times the computed eigenvectors. The steps (a) and (b) are incorporated in the subroutine "INPUT", whereas the steps (c) through (f) are carried out by the subroutine "RITE'. Flow Diagram 44/ 7/ \- Lp^ /^7/ \-5TA 42 97 CALL C<L Pox [NIN -F/,- - wr < J / )- (T'wwHe-))CO C/ I )= O7 STp N~y Npz# XJ7 N 21/ 141 AA(X)~& / K a<ct) o. \/ -C105 -

Vibration of Beams on Spring Supports Flow Diagram, Continued hi; 2o/ \ /;l\ V CTOR /7 ps 2V E R02 /3 y r,:^ <:l r (xTJl + + /&CT, J P- ) ^ L53\ ^j / - 3 X4,J NpZT rL(x Cz,J) =JI1A P= /, 2 j/, NPZ T CrA, - -t) Tzp/, t f)/> 4- V C-C7T^ (I), -) -c106 -

Example Problem No. 84 Flow Diagram, Continued FORTRAN Programs and Data SHAN S. KUO (YALE UNIV) SHAN S. KUO (YALE UNIV) $COMPILE FORTRAN. EXECUTE D003N 9 DOO3N 9 004 027 200 004 027 200 SPRING 6 SPRING 6 C C C C C C VIBRATION OF BEAMS ON SPRING SUPPORTS SHAN S. KUO 7/29/61 DIMENSION FLEX(5,5),ZMASS(5,5) TEMP(5,5),VECTOR(5t5), 1 XX(5),AA(5) WRWH(5) COMMON NSPANBETA NFIRST IS THE FIRST N AND NLAST IS THE LAST N READ INPUT TAPE 7,12, NFIRSTNLASTNRWRWH 12 FORMAT(312) READ INPUT TAPE 7,14, (WRWH(I1),I1=1,NRWRWRH) 14 FORMAT(6F6.3) DO 471 N=NFIRSTNLAST SPAN=N+1 DO 471 J=1iNRWRWH BETA=48 *70455*WRWH(J)*WRWH( J ) / (SPAN*SPAN*SPAN) WRITE OUTPUT TAPE 6,223tWRWH(J) WRITE OUTPUT TAPE 5,223,WRWH(J) 223 FORMAT(67H1VIBRATION OF BEAM ON SPRING SUPPORTS. A PROBLEM OF SHOC 1K ISOLATION/64HORIGID BODY FREQUENCY / FUNDAMENTAL FREQ. OF HINGED 2-HINGED BEAM=F5.2) CALL INPUT(FLEXZMASSVECTORTEMP XXAA) CALL RITE (FLEX9ZMASSVECTORTEMP) 471 CONTINUE -C107 -

Vibration of Beams on Spring Supports FORTRAN Programs and Data, Continued SCOMPILE FORTRAN C C INPUT MATRICES SHAN S. KUO 7/29/61 8/5/61 C SUBROUTINE INPUT (FLEX#ZMASSVECTOR,TEMPXXAA) DIMENSION FLEX(595) ZMASS(5,5) VECTOR(5,5),XX(5) AA(5),TEMP(5.5) COMMON N SPANBETA EN=N NM1=N-1 NP1=N+1 NP2=N+2 N2'N*2 C C FORM THE INFLUENCE MATRIX* II AND JJ FOR MASS NUMBERING C DO 16 I=-1N EI=I IPl=I+1 XX(IP1)=EI 16 AA(I)=EI XX(NP2 )SPAN XX( 1 ) 0 DO 15 II=lN DO 15 JJ=1iN IP1=II+1 IF(XX(IP1)-AA(JJ) ) 17t17 18 17 BB=SPAN-AA(JJ) FLEX(II JJ)= (BB/(6*SPAN))*(XX(IP1)*XX(IP1)*XX(IP1) 1 -(SPAN*SPAN-BB*B8)*XX(IP1))*(- 1) GO TO 15 18 BBASPAN-AA(JJ) FLEX(IIJJ)=-(BB/(6**SPAN))*(XX(IP1)*XX(IP1)*XX(IP1) 1 -(SPAN*SPAN-BB*BB)*XX( IP) 2 -(SPAN*(XX( IP1)-AA(JJ) )*(XX( IP )-AA(JJ) )*(XX( IP1)-AA(JJ) )/BB)) 15 CONTINUE C C SHIFT FLEX MATRIX DUE TO BEAM ACTION. BORROW VECTOR C DO 201 I=1,NP2 DO 201 J=1iNP2 201 VECTOR(ItJ)=O0 DO 202 I=1tN IP1=I+1 DO 202 J=lN JP1=J+1 202 VECTOR(IP1,JP1)=(FLEX(I J))*BETA WRITE OUTPUT TAPE 6, 777, BETA WRITE OUTPUT TAPE 5. 777, BETA 777 FORMAT(7HOBETA =F6.3) C C FLEX IS NOW ONE OF THE 3 MATRICES TO FORM THE FINAL MATRIX FLEX C DO 203 I=1,NP2 EI=I DO 203 J=1,NP2 203 FLEX( I J)=(SPAN+1.-EI)/SPAN C C BORROW TEMP, ONE OF THE 3 MATRICES TO FORM THE FINAL FLEX MATRIX C DO 204 I=1,NP2 EI=I DO 204 J=1NP2 EJ=J 204 TEMP (IJ)=(EJ-1.)*(SPAN+2.*(1.-EI))/(SPAN*SPAN) C C HAPPY FINAL FLEX C DO 205 I=1,NP2 DO 205 J31,NP2 205 FLEX(I J)=FLEX( IJ)-TEMP (IJ)+VECTOR(I J) -clo8 -

Example Problem No. 84 FORTRAN Programs and Data, Continued C C FORM MASS MATRIX C DO 114 II=1iNP2 DO 114 JJ=1.NP2 IF(II-JJ) 112,111,112 111 ZMASS(IIJJ)=l./(EN+1.) GO TO 114 112 ZMASS(IIJJ)sO. 114 CONTINUE ZMASS( 11)=O05*ZMASS(1,1) ZMASS(NP2,NP2) 0.5*ZMASS(NP2,NP2) RETURN $COMPILE FORTRAN C C WRITE OUTPUT SHAN S. KUO 7/29/61 C SUBROUTINE RITE(FLEXZMASSVECTORTEMP) DIMENSION FLEX(5,5),ZMASS(5,5),TEMP(5,5),VECTOR(5,5) COMMON N,SPANBETA K=N+2 DO 157 I=1,K WRITE OUTPUT TAPE 5,158tI,(FLEX(IJ), J1lK) 157 WRITE OUTPUT TAPE 6,158.I,(FLEX(IJ), J=lK) 158 FORMAT(23H FLEXIBILITY MATRIX ROWI3/(1H 5F10.5)) DO 701 I=1,K WRITE OUTPUT TAPE 5711,I*(ZMASS(I9J),J=lzK) 701 WRITE OUTPUT TAPE 6,711,I,(ZMASS(IJ)*J=1~K) 711 FORMAT(16H MASS MATRIX ROWI3/(1H 5F9.5)) C C FORM TEMP MATRIX = SQ. RT. OF MASS MATRIX C DO 401 I=1iK DO 401 J=1,K IF(I-J) 412,411,412 C IN THE UNIV. OF MICHIGAN USE SORT NOT SQRTF 411 TEMP(I,J)=SQRT (ZMASS(I,I)) GO TO 401 412 TEMP(I,J)=0O 401 CONTINUE C C FORM H MATRIX = TEMP TIMES FLEX TIMES TEMP C STORE FLEX TIMES TEMP INTO VECTOR ( TEMPORARY USE) C CALL MATMPY(FLEXKKTEMPKVECTOR) CALL MATMPY(TEMPKK,VECTOR,KFLEX) DO 402 I=1,K WRITE OUTPUT TAPE 5,403,I,(FLEX(IJ),J=1,K) 402 WRITE OUTPUT TAPE 6,403,I,(FLEX(IJ),J=1,K) 403 FORMAT(13H H MATRIX ROWI3/(1H 5F10.5)) C C EIGN IS A SUBROUTINE AVAILABLE IN THE UNIVERSITY OF MICHIGAN C.CO C S=EIGN(FLEXK,1,VECTOR(1,1),F) IF (S-1.5) 13,13,14 13 WRITE OUTPUT TAPE 6,12,5 12 FORMAT(15HOBAD RETURN S =F10.6) GO TO 777 14 IF(S-2.5) 15.15,16 15 WRITE OUTPUT TAPE 6,11,F 11 FORMAT(15HOSCALE FACTOR =F10.6) C C GET TRUE VECTORS C 16 DO 404 I=1,K 404 TEMP(I,I)=1./(TEMP(II)) C C USE MASS AS A TEMPORARY STORAGE C CALL MATMPY(TEMPKKVECTORKZMASS) DO 405 I=1tK WRITE OUTPUT TAPE 5,406*I,(ZMASS(IJ),J=1,K) 405 WRITE OUTPUT TAPE 6,406,I,(ZMASS(IJ),Ju1lK) 406 FORMAT(16H TRUE VECTOR ROWI3/(1H 5F10.5)) -C109 -

Vibration of Beams on Spring Supports FORTRAN Programs and Data, Continued C C C TRANSFORM THE EIGENVALUES TO FREQUENCIES DO 150 I=1,K 150 FLEX(II)=SQRT (l./FLEX4I,Il) DO 177 I=1,K WRITE OUTPUT TAPE 5,178,IFLEX(I.I) 177 WRITE OUTPUT TAPE 6,178,I,FLEX(II) 178 FORMAT(28H FREQUENCY IN SQ. RT. OF K/MI6/(1H F30.6)) 777 RETURN C C MATRIX MULTIPLICATION C SUBROUTINE MATMPY(ANoMBLC) DIMENSION A(5,5) B(5,5) C(5,5) DO 75 I=1,N DO 75 J=1tL C(I,J)=0. DO 75 K=1,M 75 C(IJ)=C(IoJ)+A(I K)*B(K#J) RETURN $ DATA 030302 0.5 4, Computer Output VIBRE TIO t OF BERM i-OH SFRIHG SUPPORTS. FPROE:LEt OF SHOCK ISOCLATIOH RIGID BODV FREQUEHC'.' BETT = 0. 190 FLE>',IBILITYV MFPTRIX ROWd 1.00000 0.75000 FLEXIBILITY M1RTRIX FRO 0.75000 0. 76769 FLEXIBILITY MATRIX OWl '.50000 0.6744 FLEXIBILITY MATRIX ROWil 2. 25000. 485 9: FLE IBILIITY MATRIX ROW Q. 00. 25C00' MRSS MATRIX ROW 1 0.12500 0. 0. M SS MATRIX ROW 2 0. 0.25000 0. MASS MATRIX ROW 3 0!. 0. 0. MiSS MATRIX ROW 4 O. O. 0. MASS MRTRIX ROW 5 0. 0. 0. H MATRIX ROW 1 0.12500 0. 1325 H MPTRI ROW 2 O. 13 258 0. 1 9192 H MPTRIX ROW 3 0.08839 0.160 H M'TRIX ROW 4 0.04419 0.12150 H MRTRIX ROW 5 0i.. 0.0441 TRUE VECTOR ROW 1 1.60933 -1.18746 TRUE VECTOR ROW 2 0.83966 1.1 3797 TRUE VECTOR ROW 3 -0. 00000 -0. 00000 TRUE VECTOR ROW 4 -0. 83966 -1.1 3797 TRUE VECTOR ROW 5 -1.60933 1.18746 FUDARMENTPL FREQ. OF HINHED-HINGED EEAEl= 0.50 1 0.67440 0.48598 3 0.75367 0.67440 4 0.67440 0.76769 5 0.50000 0. '5000 O. 0. 0. 0. 25000 0...,..-._, 0.25000. 0. 0. 0.25000 0.50000 0.75000 1.00000 0. 0.12500 0.08839 0. 1 6860 0. 18842 0.16860 0.08839 -1.71950 0.16622 0.99402 0.16622. 04419 0. 0. 2150 0.0441.I6860 O.08839 -L. 191 92 -... J-L1c 325 -,:.i 1325 0. 12500 0. 33 96 0. 58978 1.02312 -0. 9208 1.10045 1.34 199 1.02312 -0. 96208 -1.71950 0. 3396 0.58973 -C110 -

Example Problem No. 84 Computer Output, Continued FREQUENCY IN SQ. RT. OF K/M 1 2.29S4......_.._________ ----------. --- —------------- FREQUENCY IN SQ. RT. OF K/M 2 -_ ~_ -.. -13- 8.55_16D-_______ --- —---------— _ FREQUENCY IN SQ. RT. OF K/M 3 5_ 91 5994 FREQUENCY IN SQ. RT. OF K/M 4 1. 294614 __ —___ _ -------------------------------- FREQUENCY IN SQ. RT. OF K/M 5 — ___ ___-_ _ _ ___ _24. 901JiS4 --- —------------— _______- — __________ _ _v IBRIT.I.Ot _OE.._BERMJ_ S._ PRI NG S_. PTS PBLEM SUP RTS_ P B LEE HO CKLSO _LLO N _ RIGIn BODn FPFRUFNCY / FlINDAMFNTrt FPFO- nF HTNffn-TMTNGFn P.FQM- 4. 0 BET..I2._76 __ --- —---- _.7 --- —------ -------------- FLEXIBILITY MATRIX ROW 1. —._ Il....D._,-_7QSt-__ — O5.. Q0o ---_ Qno --- _-______________ --- — -------------- FLEXIBILITY MATRIX ROW 2. 75noo 9. 7S71 n 11. 6, 1 4i. 7 47775 nrt. ^nnn FLEXIBILITY MATRIX ROW 3..-__-____Q5Q-__O.. LL._6_J_4. - __ J73_45 _1~6J- ___._-__.LQ ------------------— ____ — ___ FLEXIBILITY MATRIX ROW 4. —_ _ __2..50O_..-.- 7.-_47775___1 --- L46____O675n — __ -- ---------— 5LO- FLEXIBILITY MATRIX ROW 5 0c 0. 250-Inn n. snnn n- 70innn 1 nnnn MASS M.TRIX ROW 1 _....... 0 s- _10.. -.....0 ---— a --------------------— ______________ MASS MATRIX ROW 2 -_-....... O....25Q....-__-___._.Q._._.-.._____,____ _______ ------- MASS MATRIX ROW 3 0. n..o25n00 n n, MASS MATRIX ROW 4 -. -. _..............-____Q.__________ _____________ --- —----------— __ ------ MASS MATRIX ROW 5............. _......... _.. J2.. 25Q- ---— ______________ _ ----- H MATRIX ROW 1 0. 12500 0.13258 0.08839 0.0441 - 0. H MATRIX ROW 2.13258.6 -- -...2.4...92_2. 536_. ------- _,_ --- —-__1_-_ t-9 H MATRIX ROW 3 _. 0.08839 2. 91 536.4_.!i1 Z -S6._ --- —------- ____________ H MRTRIX ROW 4 __. 04419 1 86944 2.9153.?-4q98 II 1 9.q H MATRIX ROW 5 0. 0.0 449....... 0.889....... _.1-325.. -L-L2.l-= — ------------- TRUE VECTOR ROW 1 1.96434 0.37601 0.....0517_-. 7-.L7 r79.J4.5..__0_____ --- —-- TRUE VECTOR ROW 2 -0. 26588 _ -. 8- i 0. 07.,o 2 -? -7nin -n. 7A ---TRUE VECTOR ROW 3.2. 0.000000 -0. 000000. 1. 283. — _ -l -D6.36__ t.l.............. - ------- TRUE VECTOR ROW 4 0. 26588 -1. 38-00 1.00732..7.. _ 0 10 -.._..3-,i _..._........... TRUE VECTOR ROW 5 o -.96434....-_O_.360_ n-_ n.....5.l_7EL 1_ - 1.- 3.7979, 44 r:,9~ FREQUENCY IN SQ. RT. OF K/'M 1 3, 041763.. _.-._...... _..... FREQUENCY IN SQ. RT. OF K/M 2 1,305481............... —. --- —--- FREQUENCY IN SQ. RT. OF K-NM 3 0.. 345544........ FREQUENCY IN SQ. RT. OF K/.M 4 2. 530636 FREQUENCY IN SQ. RT. OF K M 5 3.407844 -Clll

Vibration of Beams on Spring Supports Discussion of Results The natural frequencies are equal to the computed values times the square root of k/m. The modal shapes are normalized in the subroutine BIGN according to the following relation: n+2. m. i 2= m j = 1, 2...,n+2 Figure shows the computed normal modes for two different alues. Figure 2 shows the computed normal modes for two different &6R/aH values. 1 44 RH i1 = 1.295 k M2 =2.293 2 V 1 = 0.346 - 1 m.37 () = 1.305 2 m is 1.377 t.ct~ ooc - ' 3 = 5.916 k 3 ~~m M3 = 2.531 3 Vm / \ I / - 1. oS I 9sll 4 = 13.855 a5 = 24.902 m 5 4) = 3.042 0 5=g 3.408 i Figure 2. Normal modes and frequencies __________ c R H = 0.5 ----- - 4.0 -C112 -

Example Problem No. 84 Critique A problem of this type is of great importance in the design of underground structures to isolate shock loads. Many textbooks on vibration do not mention Jacobi's method. With the advent of high speed digital computers, this method becomes very important. Its algorithm is simple and all the eigenvalues and eigenvectors are computed at the same time. Further extensions of this problem include the computation of bending moment matrix, shear matrix, reaction coefficients, and the moment-time relation for a given shock function. VII. REFERENCES 1. Organick, E. I., A Computer Primer for the MAD Language, University of Houston, Computing and Data Processing Center, 1961. 2. Organick, E. I., A Computer Primer for the FORTRAN Language, University of Houston, Computing and Data Processing Center, 1961. 3. Project on the Use of Computers in Engineering Education, First Annual Report, The University of Michigan College of Engineering, Ann Arbor, August 2b, 1960. 4. Project on the Use of Computers in Engineering Education, Second Annual Report, The University of Michigan College of Engineering, Ann Arbor, December 15, 1961. 5. Westervelt, F. H., Use of Computers in Mechanical Engineering Education, Project on the Use of Computers in Engineering Education, The University of Michigan College of Engineering, Ann Arbor, June 1, 1962. 6. Wilson, R. C., Use of Computers in Industrial Engineering Education, Project on the Use of Computers in Engineering Education, The University of Michigan College of Engineering, Ann Arbor, January 27, 1962. -C113 -

)r^ teo N 4 _~s~sY