THE USE OF COMPUTERS IN ENGINEERING MECHANICS EDUCATION Department of Engineering Mechanics The University of Michigan June 30, 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.

ABSTRACT During the past two years, the faculty of the Engineering College of The University of Michigan has been exploring the integration of computers into the various undergraduate curricula. This report describes the Engineering Mechanics program at the University and discusses the role which the faculty of the Department feels computers should play in this curriculum. Digital computer solutions to engineering problems have been demonstrated in Departmental courses, although students have not been assigned digital computer problems as part of their homework. The Department which operates an analog laboratory of its own, uses analog computers extensively for both demonstration and laboratory purposes. This report contains a selected set of three example problems with complete computer solutions prepared by a Departmental faculty member. These may be considered as a supplement to the 96 example engineering problems, including some related to Engineering Mechanics subject areas, which have been published previously by the Project. -F1

Table of Contents Page I. Introduction F3 II. The Program in Engineering Mechanics F4 III. The Use of Computers in the Mechanics Curriculum F6 IV. Example Problems F10 97 Numerical Solution of the Harmonic and Biharmonic Equations Fll 98 Joukowski Airfoil F23 99 Principal Axes of a Second Order Tensor F27 Tables IFA Elementary Course Work F5 IFB Intermediate Course Work F5 IFC Advanced Course Work in Mechanics on the Undergraduate- F5 Beginning Graduate Level IFD Advanced Courses in Mechanics on the Graduate Level F5 IFE Advanced Mathematics Courses Frequently Elected F5 IIF List of Example Problems F10 -F2

USE OF COMPUTERS IN ENGINEERING MECHANICS EDUCATION* I. INTRODUCTION The Ford Foundation Project on the Use of Computers in Engineering Education at The University of Michigan has enabled the faculty and students of the College of Engineering to test on a fairly large scale the role of the computer in undergraduate engineering. The present report attempts to set forth the views and efforts of the Department of Engineering Mechanics in this regard, in the hope that others can benefit in some measure from our experiences. It is recognized that the makeup and aims of the departments of mechanics vary to a large extent from school to school, and hence it is thought worthwhile to set forth our program and course offerings first, so that a frame of reference is in the reader's mind from the start. Our conclusions and opinions can then be modified to suit the conditions at the reader's institution. This report should be read in conjunction with the final report of the project. It is not intended to be self-contained. It should be appreciated also that due to changes which can be expected in the expanding fields of computer technology and numerical analysis, and the varying needs of science and industry in general, any and all of the ideas and conclusions presented here can become obsolete very rapidly. Whatever impression the reader receives as to the use of computers in the mechanics curriculum, it behooves him to be aware of the possibilities and potentialities as they develop, and to revise his ideas accordingly. One special elaboration on these last two sentences is of particular note. Much of the present use made of computers by engineers is of the "slide rule" type. That is, the same methods are used on the high speed computer as were found effective on a desk calculator, and the only advantage gained is in speed and relief from tedium. Recently in several areas of science, users of the computer have broken away from utilizing just the arithmetical hardware of the computer and have made extensive use of the logical or decision-making side of it. In the long run, this is probably the most exciting side of the computer, and the one to watch most closely. The direct application of this to mechanics is not immediately clear, but one can perceive potentialities which could open doors now almost impenetrably closed (such fields as the study of turbulence, fatigue, crack propagation, to name but a few of the possibilities). Here, considerable research is still necessary before any tangible rewards are gained. In the following, when the word computer is used without any modifying adjective, it is understood to mean digital computer. * This report represents the group opinion of the Department. In particular, Professors S. K. Clark, W. R. Debler, R. A. Dodge, J. H. Enns, W. P. Graebel, R. M. Haythornthwaite, B. Herzog, and M. J. Kaldjian have made contributions to it. The report was prepared by W. P. Graebel with their assistance. -F3

II. THE PROGRAM IN ENGINEERING MECHANICS The undergraduate engineering student would normally elect the engineering mechanics program at the beginning of his sophomore year. He would start with the usual elementary sequence of courses (Table IFA). Following these courses, at the second semester of the junior or senior level he takes a series of intermediate courses (Table IFB) which lay the foundations and give an overall view of the main branches of mechanics. He is also required to select a minimum of seven credits of advanced subjects in engineering mechanics, in which he has considerable latitude (Table IFC). These courses are usually offered at least once a year, so that by proper planning, a student can elect a sequence of courses in one particular area of mechanics. In addition, the student can select seven to nine credits in any college of the University, and a group option of twelve or fourteen hours from any other program in the College of Engineering. (Frequent choices for the latter are aerodynamics, electronics, instrumentation, metallurgy, naval architecture, nuclear engineering, physics, structural engineering, etc.) It is not uncommon to find that a student obtaining a bachelor's degree in engineering mechanics will, with perhaps an additional semester's work, receive an additional B.S. degree in applied mathematics, engineering physics, or in various other curricula available to him in the College. As is evident from the course outlines, the Engineering Mechanics Department does not offer any design courses as such. The student is required to take design courses in other departments, and can elect more in his group option if he so desires. The Department has felt that its principal contribution is not that of producing designers, but rather the training of the engineer in depth in the fundamentals of his field so that he is aware of the foundations and realizes the possibilities and limitations of mathematical models in his work. The student with such training at the bachelor's level receives ready acceptance from industries dealing with diverse technical interests. As might be expected, a large percentage of the students decide to continue their studies in graduate schools throughout the country, most continuing in mechanics, although a significant number go into allied fields such as mathematics, systems engineering, communications, physics, nuclear engineering, space sciences, etc. The student going to these other fields usually finds that his early training allows him to compete with his fellow students in these other disciplines at their level. Since many of the students do enter graduate school upon completion of their bachelor's training, our course work is designed with this in mind, and a brief description of this program is perhaps necessary to give a complete overall picture. E.oM. 412, 422, 441 (Table IFB), while intended as terminal courses at the B.S. level are beginning courses for students entering our our graduate program from other schools. Beyond these, no other specific mechanics courses are required at the graduate level and the student is free to specialize as he and his advisor see fit -F4

Use of Computers in Engineering Mechanics Education TABLE IFA Elementary Course Work Credits E.M. 208 or 218 Statics 3 E.M. 210 or 219 Strength of Materials 4 E.M. 212 or 402 Laboratory in Strength of Materials 1 or 2 E.M. 343 or 345 Dynamics 3 E.M. 324 or 326 Fluid Mechanics 3 or 4 TABLE IFB Intermediate Course Work Credits E.M. 403 Experimental Mechanics 2 E.M. 412 Intermediate Mechanics of Materials 3 E.M. 422 Intermediate Mechanics of Fluids 3 E.M. 441 Intermediate Mechanics of Vibrations 3 TABLE IFC Advanced Course Work in Mechanics on the Undergraduate-Beginning Graduate Level Credits E.M. 411 Structural Mechanics 3 E.M. 413 Photoelasticity 2 E.M. 416 Stress Analysis 2 E.M. 514 Theory of Elasticity I 3 E.M. 515 Theory of Plates 3 E.M. 518 Theory of Elastic Stability I 3 E.M. 519 Theory of Plasticity I 3 Dynamics E.M. 542 Advanced Dynamics 3 E.M. 543 History of Dynamics 2 E.M. 544 Dynamics and Stability of Rotors 3 E.M. 545 Vibrations of Continuous Media 3 E.M. 547 Theory of Gyroscopes 2 Fluid Mechanics E.M. 522 Mechanics of Inviscid Fluids I 3 E.M. 523 Mechanics of Viscous Fluids I 3 E.M. 529 Advanced Laboratory in Mechanics 2 of Fluids Thermodynamics E.M. 527 Thermodynamics 2 TABLE IFD Advanced Courses in Mechanics on the Graduate Level Solid Mechanics Credits E.M. 714 Theory of Elasticity II 3 E.M. 715 Theory of Shells 3 E.M. 718 Theory of Elastic Stability II 3 E.M. 719 Theory of Plasticity II 3 Dynamics E.M. 714 Theory of Vibrations 2 E.M. 745 Wave Motion in Continuous Media 3 Fluid Mechanics E.M. 721 Mechanics of Inviscid Fluids II 3 E.M. 723 Mechanics of Viscous Fluids II 3 Other E.M. 707 Theory of Continuous Media 3 TABLE IFE Advanced Mathematics Courses Frequently Elected Credits Math. 749 Methods of Partial Differential Equations 3 Math. 750 Methods of Mathematical Physics I 3 Math. 751 Methods of Mathematical Physics II 3 Math. 757 Special Functions in Classical Analysis 3 Math. 777 Tensor Analysis 3 -F5

The Ph.D. candidate is, however, required to demonstrate proficiency in a written qualifying examination in the subjects covered in Tables IFA and IFB. He also must take at a later date, an oral examination in the area in which he majored plus two minor areas, one in the Engineering Mechanics Department, the other outside of the Department. The minor in the Department usually consists of two courses selected from one of the groups in Table IFC. The major consists of course selected from groups in Tables IFC and IFD, plus individual reading courses and advanced courses in other departments. Mathematics is heavily stressed, and the student usually finds it desirable to elect several graduate mathematics courses as well. Courses which have been elected most regularly in recent years are given in Table IFE, III. THE USE OF COMPUTERS IN THE MECHANICS CURRICULUM During the period when the Project on the Use of Computers in Engineering Education was in effect at The University of Michigan, several points involving the use of digital computers in undergraduate education became apparent. Since these apply to any student (and teacher) user of the computer, they will be set forth first, serving then as a basis of discussion when the role of computers in the mechanics curriculum is discussed. Discussions of other factors can be found in the various reports of the project. The Computer Requires a Precise Statement of the Problem Since communicating with the computer is somewhat similar to communicating with a threeyear-old child, the student finds at an early stage that he must tell the computer exactly what he wants to do in order to obtain the correct results,. This requires an understanding of the mathematics of the problem by the student and also a clear statement of the problem in mathematical language, desirable goals in the training of the engineer. Of course, this clarity of problem statement should be the goal in all teaching, whether a computer is present or not. The teacher, however, as a human being, will usually accept the work if small errors are made. The computer is a much more stern taskmaster, and will produce satisfactory results only if the problem is presented to it letter perfect. Like most good things in life, this need for preciseness is not a complete blessing. First of all, it tests only one facet of the student's knowledge, that is, his problem solving ability in a particular subject. It does not in general test his understanding of the fundamentals, since he may be led to a solution by the fact that the problem appeared at a particular part of the course. Also, after the first (usually relatively small) percentage of his time spent on this particular problem, the student's concern is not in furthering his understanding of the engineering involved, but rather with the details of computer programming. Even with a problem-oriented computer language such as MAD, considerable practice, experience, and care are needed to make sure that "punctuation" and similar details do not cause the computer to return undesirable results. -F6

Use of Computers in Engineering Mechanics Education The Student Should Become Familiar with the Computer as a Mathematical Tool Just as we expect the student to be familiar with techniques such as separation of variables, power series solutions of differential equations, etc., so he should also recognize that numerical methods are powerful tools in problem solving. Through experience he should learn when one method has advantages over another, and what factors play important roles in determining how accurate his answer is. This requires then,,that besides learning to program in a computer language, the student learn numerical analysis as well, including error analysis along with numerical techniques. In fact this knowledge of numerical analysis may be the most significant aspect in the long run. Knowledge of a computer language is necessary in that it provides laboratory experience, but it is the least permanent part of a student's education in that it is subject to rapid obsolescence. The Computer Serves as a Laboratory An "exact" solution to a complicated problem provides some inherent satisfaction to many people; however, it is not always possible to see just what role the various parameters play in governing the behavior of the system to which the solution applies. The computer, either through a numerical solution of the basic problem or through evaluation from the exact solution, can provide a graphical presentation of the solution which is frequently more easily understood. Thus a computer can serve much the same function as, say, an elementary laboratory in strength of materials does in the student's education. Use of the Computer by Students Requires Considerable Availability of Knowledgeable Staff An instructor considering the assignment of a problem for execution on the computer should first be sure that someone (preferably the instructor himself, or at least someone familiar with the engineering problems) is available to help the students interpret error returns, core dumps, and the like. If the degree-of-problem-difficulty times number-of-students is at all large, it should be anticipated that this help will be required for considerable periods of time over an interval of several weeks per problem. Instructor time, not computer time, is perhaps then the chief "cost" factor involved. Student Problems on a Computer Require Days if not Weeks from Assignment to Completion During the existence of the Computer Project, the approach to the computer at The University of Michigan was made relatively easy for the student. Carrier service from locations in the engineering buildings provided delivery of programs to and from the Computing Center twice a day. Keypunches were provided in the engineering buildings as well as at the Computing Center. Personnel were available at several locations for consultation in case of student or faculty difficulty. High speed printers and special services were provided at the Computing Center. -F7

In spite of all these, 24 hours seemed to be the minimum time from delivery to receipt of the program. Breakdown of equipment and computer load at times lengthened this to as much as one week or more, so that average delivery-receipt time at the middle or end of the semester was probably at least 48 hours. Since the average number of passes through the computer for students with some acquaintance with programming, running reasonably short programs, was estimated to be in excess of three, it is easy to see the effect of this on the student. By the time the majority of the students in a class ran a simple problem on the computer, the class would have progressed to another topic, and hence the impact of the computer's contribution would be somewhat dulled. If, in addition, some of the students are unfamiliar with the particular language available, an appreciable amount of class time may be lost in an introduction to this languag A Certain Number of Students Will Become Completely Engrossed in the Computer to the Point of Losing Interest in Engineering The computer can be an engrossing object in the student's life, so much so that he feels compelled to devote much of his time and energy to this fascinating device, usually to the detriment of his engineering studies. This is perhaps due to the fact that in a reasonably brief period of time he may gain a high degree of proficiency and find himself in great demand by his teachers as well as his fellow students. This success quite often proves hazardous to his other studies and delays his progress towards a degree. The Computer Solution Gives only the "Expected" Answer While this is true to some degree of any method, there is a further danger when the computer is used. When a numerical method is decided upon, this pretty well fixes the type, or behavior, of the solution which will be found. An "Answer" based on this method will usually then be obtained, which however, may have little bearing on the true physical problem. To illustrate, suppose a problem in the form of a differential equation is to be solved which in some region has either a boundary-layer type behavior (that is, the solution changes rapidly in almost a discontinuous manner) or perhaps the family of solutions become so close together in a region that the c'omputer can jump easily from one solution to another. If the student is unaware of this, he might start off with a mesh size of, say, 0.1, on the next trial change this to 0.05, and next to 0.01. The solution may turn out to be relatively insensitive to this size of change; ergo, the student concludes that his method has converged to the solution. If, however, he would continue to reduce his mesh size until it was of the order of the boundary layer thickness (or the "distance" between radically differing solutions), he would find that his previous efforts were deceiving, and gave him false confidence in his "solution." This seems to emphasize again the need for proper training in error analysis and an understanding of the limitations of numerical computation. The above could happen even without a computer, and perhaps best serves as an indictment of a cook book approach to numerical computation. The reason that this point is of particular danger -F8

Use of Computers in Engineering Mechanics Education when a computer is used, is because when one is performing a hand computation, the relative size of various terms in an equation quite often becomes apparent with little effort, and curiosity as to where various terms become important is easily aroused. In the computer, however, the inner intricacies of the computation are hidden, and unless the programmer is alert to the dangers beforehand, the lack of imagination in the machine gives the wrong solution. Based on the above considerations and others mentioned in the various project reports, the Department of Engineering Mechanics presently makes only limited use of the digital computer in its existing undergraduate and graduate courses. For the courses in Table IFA which perform a service role besides giving an introduction to mechanics for our own students, there is simply not sufficient time to introduce any computer work without reducing the amount of mechanics taught in these courses. It may very well be desirable to use the computer in these courses, and if additional time were made available to us we would be willing to do so. We do not, however, believe it desirable to reduce the mechanics content of these courses below the level now being given. For the other courses, instructors' computer solutions have been used to illustrate what various solutions look like when numerical computations are carried out, and also the digital computer has been used to replace some hand calculations which were previously used. So far it has been left optional to the student as to whether he would use a desk calculator or the IBM 709. Those familiar with programming choose the latter, and in the process, convince friends as to the desirability of learning programming. When more students come to these courses with the ability to program fairly involved problems, undoubtedly more use will be made of the digital computer. However, the content of the existing courses were set up before the Computer Project made the computer so accessible. The arrangement and content of these courses will soon come up for review and it is likely that new courses decided on will use the computer in a more integrated fashion. In summary, then, we believe that our duty as a department is to train the student in the fundmentals of mechanics. We have adopted a wait and see attitude as to whether the computer can further these aims. It is becoming more apparent that it can do so, but at this time the advantages to be gained by widespread use in the undergraduate program do not seem to be economical timewise. However, through the aid made possible by the Computer Project, our future courses can expect to incorporate in varying degrees the use of the computer as a means of furthering our aims. Analog Computers The electronic differential analyzer (EDA) has been used by engineers long enough so that it is a familiar teaching tool. These computers are commonplace enough so that little need be said, but perhaps a brief comparison with the digital computer is in order. -F9

The EDA is suited to a much smaller class of problems than the digital computer, namely, the solution of ordinary differential equations. This lack of versatility goes hand in hand with easier operation; hence the student can be introduced to the principles of EDA operations in a few minutes, and can learn to operate a particular machine in an afternoon. The nuances of how to handle infinities, two-point boundary value problems, and nonlinear terms takes more time, but available manuals make such learning relatively simple compared with the digital computer. The quickness and ease of solution and presentation makes the EDA a tool well adapted to classroom demonstration. The Department of Engineering Mechanics presently owns two 10-amplifier computers plus a 24-amplifier model (the latter was provided through the Computer Project and is available to other departments in the College of Engineering). These are presently used to demonstrate to appropriate classes the behavior of various multi-mass systems, uniform vibrating strings and beams, flow in the boundary layer, the response of nonlinear materials, and other problems with only one independent variable. The students are also encouraged to use the EDA on their own for appropriate problems, and have relatively free access to it. Solutions on the EDA can be made to exhibit some of the dangers of "numerical' solutions in a simpler and more visual manner than on the digital computer. Thus, some of the benefits of the digital computer can be obtained with relatively little financial or time expense. IV. EXAMPLE PROBLEMS Included in this section are three problems and their solutions programmed in the MAD (Michigan Algorithm Decoder) language. They are listed in Table IIF. TABLE IIF List of Example Problems Number* Title Author Page 97 Numerical Solution of the Harmonic and W. P. Graebel Fll Biharmonic Equations 98 Joukowski Airfoil W. P. Graebel F23 99 Principal Axes of a Second Order Tensor W. P. Graebel F27 * These problems may be considered as a supplement to problems 1 through 96 published in previous reports of the Project. -F10

Example Problem No. 97 NUMERICAL SOLUTION OF THE HARMONIC AND BIHARMONIC EQUATIONS by William P. Graebel Department of Engineering Mechanics The University of Michigan The harmonic and biharmonic equations appear many times in problems in mechanics. Potential flow, stress distribution, deflection of a membrane, steady-state temperature distribution are examples. The solution of these by relaxation methods is well known due to the efforts of Southwell and his followers, hence there seems to be little need in recapitulating the problem. The techniques used by these workers to speed the convergence now seem to be of incidental interest if machine computation is to be used, since the added complexity to the program would only be worthwhile if there was need for an extreme number of mesh points. The arithmetic operations here are trivial; the main problem is to identify whether a point is interior, on the boundary, or exterior to the region of interest, and also how close are the neighbors of the point if any of the neighbors are boundary points. The present programs approach the problem by taking the simplest case, and then modifying this to meet more complex boundaries. Case 1 (Program Relaxati) Here we consider the case where all mesh points for a square grid lie on the boundary, hence all boundaries must be a series of straight line segments intersecting at prescribed mesh points. The information necessary to the program is simply the location of the boundary (given by the BOUND vector) and also the value of the function (PSI) at the boundary. For convenience, we also read in the maximum number of times we will allow a point value to be relaxed (NMAX) (to avoid a runaway program if something is wrong) as well as a means of stopping the program when we are content with the answer, that is, when the maximum change made in going through a complete iteration is less than ERRMAX. To avoid programming complexities, a code is used on the BOUND vector, and data is read in for every point on our net. The code used was a simple one, BOUND=O for an interior point, 1 for points on the boundary, 2 for exterior points, and 3 for interior points which have at least one nearest net neighbor which is an exterior point. Actually, for Case 1, any integer other than zero will work instead of 1 or 2, and 3. The code s stated here for use in Case 3. PSI is read in for every point, whether or not it is a boundary point. For nonboundary points, any value will do as it is never used. -Fll

Example Problem No. 97 The program now goes through each point in a rectangle XMAX by YMAX in a straightforward manner, checking to see the nature of a point, then either going on if it is not an interior point or else taking the average value of its neighbors if it is an interior point. If NMAX relaxations are undergone without meeting the ERRMAX condition, the program prints'out the last value of PSI at each point as well as the difference between this value and the one computed on the previous iteration. The problem presented represents the irrotational flow in a channel which suddenly increases in width. Program variables for quantities not described in the text are listed below and are applicable to all four cases. N Number of iteration NMAX Maximum number of iterations allowed C Switch to stop iteration when change becomes small enough PSI Solution of harmonic or biharmonic equation PSIB New value of PSI XMAX, YMAX Size of rectangle enclosing the boundary PSR, PSL, PSU, PSD Value of PSI to right, left, above, and below (I, J) NUM Numerator of Equation (1) DEM Denominator of Equation (1) T Entry into random sequence AVE Average number of throws needed to get to boundary A flow diagram for the program (RELAXATI) to solve problems of Case 1 is shown below. XMAX, YMAX, PSI(OO),... BOUND OO),... N=O START - N=O NMAX, ERRMAX PSI(XMAX,YMnX) B YMAXl ST RT X=i1, Y=i1, A P BOUND(X, Y)/O - N=N+i, C=0 -F12

Numerical Solution of the Harmonic and Biharmonic Equations / MID PSIDXY) BETA -- NAX X=O,1 Y=O,1 RROR(X,Y) MID FNAL x > XMAX Y >YMAX 4 x X=X+1 ALPHA <YMAX TES=T,YALPH F LAST N FINAL The MAD program is shown below. RELAXAI'IUN MEIHOLj FOR 2D LAPLACE EQUATION INTEUER N, C, X, Y, XMAX, YMAX, NMAX, BOUND READ FORMAT ABLE, XMAX, YMAX, NMAX, ERRMAX VECIOR VALUES ABLE = $3I5,F10.0*$ UIMB(2) = YMAX + 3 DIMB(3) = YMAX + 1 READ FORMAr EASY, PSI (O,O)...PSI(XMAX,YMAX) VECTOR VALUES EASY=$(8F10.6)*$ I)IMENSION PSI(625, _DIMB 1(.) ) READ FORMAf FOX, BOUND(O,O)..B OUND(XMAX,YMAX) VEC'[iLR VALUES FOX=$(4012)*$ DIMENSION BOU.ID(625,DIMB(1)), DIMB(3) VECTOR VALUES DIMB(1) = 2,27,11 N'= 0 START x = 1 N= N + 1 C = O ALPHA W.IHENVE BOUUN.)(X,Y).NE. O, TRANSFER TO DELTA PSIK = (PSI(X - 1,Y) + PSI(X + 1,Y) + PSI(X,Y - 1) + PSI(X,Y I + 1) )/4. ERROR(X,Y) =.ABS.(PSIB - PSI(XY)) DIM'ENSION ERROR(625, DIMH(1)) WHENEVER ERROR(XY).G. ERRMAX, C = 1 PSI(XY) = PSI3 WHENEVER N.G.NMAX, TRANSFER TO BETA X = X + 1 TRANSFEK r T ALPHA BETA PRINT FORMAT INTER, NMAX VECTUR VALUES INTER = $18HO EXCEEDED NMAX = I5*$ _ THROUGH MID, FOR X =. o,1, X.G. XMAX THROUGH MID, FOR Y = O,1, Y-G.YMAX M_ID PRINI' FORMAT ER, X, Y, ERROR(X,Y), X, Y, PSI(X,Y) -F13

Numerical Solution of the Harmonic and Biharmonic Equations MAD Program, continued VECTOR VALUES ER = $S2, 6HERROR(I3,H,I3,4H) =,E11.4,S10O,4H 1PSI(,I3,H,,I 1'3,4H) =,F10.6*$ TRANSFER'r0 FINAL DEL[A WHENEVER X.L.XMAX X =X + I x(=X+l TRANSFER TO ALPHA END OF CONDITIONAL WHENEVER Y.L.YMAX X = Y =Y+ I TRANSFER 0TO ALPHA END OF CONDITIONAL TESI WHENEVER C.E., TRANSFER TO START LAS'[ THROUGH LAST, FOR..X ='I_0X.G.XM.A.X... THROUGH LAST, FOR Y = O,1,Y.G.YMAX LAST PRINT FORMAT OUT, XYPSI(XY) VECTOR VALUES OUT $iHOS2,4HPSI(I5,iH,5,2H)=F1lO.6*$. PRINT FORMAT U, N VECTOR VALUES = $1HOS2,24HNUMNBER OF RELAXATIONS = 15*$ FINAL END OF PROGRAM A typical set of computer output is shown below. -.- EXCEEDED -— N M AX ERROR( O, 0) =.1885E-36 PSI( O, 0) =.000000 ERROR( 0... I) =. 1885E-36 PSI ( -__.200000 ERROR( 0, 2) =.1885E-36 PSI( 0, 2) =.400000 ERROR( 0, 3) =.188'E-36 PSI( 0. 3) =.600000 ERROR( 0, 4) =.1885E-36 PSI( 0, 4) =.800000 ERROR( 0, 5) =. 1885E-36-___ _____Y( PSI( 0, 5) = 1.000000_ - ERROR( 0, 6) =.1885E-36 PSI( 0, 6) = 7.000000 ERROR( Q,.7_) =..18.85E-3.6. PSI.._P _ I____0 _. 7) 7 =:0.000000 _L ERROR( 0, 8) =.1865E-36 PSI( 0, 8) = 7.000000 ERROR(,j 9) =.1885E-36 PSI( O0 9) = 7.000000 ERROR( 0, 10) =.1885E-36 PSI( 0, 10) = 7.000000 ERROR( 1, 0) =.18_85E-_36. -— PSI( 1, 0)=.000000 ERROR( 1, 1) =.2906E-06 PSI( 1, 1) =.199967 ERROR( 1, 2) =.4917E-06 PSI 1___ 2) =......399946_.. ERROR( 1, 3) =.4694E-06 PSI( 1, 3) =.599946 ERROR( 1, 4) =.2831E-O6 PSI( 1_ 4) =.799967 ERROR(- 1, 5) =.1885E-36 PSI( 1, 5) = 1.000000.RRUR( 1, 6) =. 13a5E-36 PSI 1, ___ 6).= 7.000000._._, _ _ ERROR( 1, 7) =.1885E-36 PSI( t1, 7) = 7.000000 ERROR( 1, 8) =.1885E-36 PSI( 1, 8) = 7.000000 ERROR( 1, 9) =.1885E-36 PSI( 1, 9) = 7.000000 ERROR( l, 10) =.1885E-36 PSI( 1, 10) = 7.000000 ERROR( 2, 0) =.1885E-36 PSI( 2, 0) =.000000 ERROR( 2, 1) =.6743E-06 PSI( 2, 1) =.199922 ERROR( 2, 2) =.1043E-05 PSI( 2, 2) =.399873 ERROR( 2, 3) =.1006E-05 PSI( 2, 3) =.599874 ERROR( 2, 4) =.5960E-06 PSI( 2, 4) =.799922 ERROR( 2, 5) =.t1885E-36._.PSIC.( 2_, 5) = 1...__.o.000.000..... ERROR( 2, 6) =.1885E-36 PSI( 2, 6) = 7.000000 ERRORC 2, 7) =.1885E-36 PSI( 2, 7) = 71.000000 ERROR( 2, 8) =.1885E-36 PSI( 2, 8) = 7.000000 ERROR( 2, 9) =.1.885E36 PSI( 2,r 9) = 7.000000.. ERRORC 2, 10) =.1885E-36 PSI( 2, 10) = 7.000000 -F14

Example Problem No. 97 Case 2 (Program Biharmon) This program is essentially the same as the preceding, except it is used for solving the biharmonic equation with boundaries of the type given in Case 1. The main difference is in the arithmetic operation necessary in computing a new value of PSI, and also the boundary conditions, since the normal derivative of PSI is also given on the boundary. These data were incorporated by simply inserting data in the form of MAD statement cards (cards 19-38) which took care of these boundary conditions. All else is similar to the preceding case. The example here was taken from Timoshenko and Goodier, Theory of Elasticity, McGraw-Hill, 1951, pp. 483-489. The flow diagram for Case 1 is to be used with the change indicated below going from 1 to 2. CALC PSI FOR EXTERIOR OF RECTANGLE The MAD program for Case 2 is shown below. E L -i'-~R T' f,,, T F- Hi,;'-,'i.- E _'i U I' H Ti E * E,: * "''i:.'.:.i!*414:',-,1 i'ii M i:', /1~ri-~.E i Ti L E.F F -' i' i- ir- E.-:! i;i 0flf R p': D I M E:':.' 2 -:.:-:: + 1. _______..j.. _HP i-i X -i r. i "'. -.T'..'''.-: l'- Ei i,: 0 _,.:,, _ i —,'.,.F i-.; I,E C T 0 R-,'. UF,. _. - - - - - - - ~ - -. L.' FT' ".' E~I Ms E;.'::7! ":..- I 2 4: 2! - — TRR T. - i V P- i..................... L<' T':= — b: FI L:;z i,: - = F-:, --. F -, 4: I".-:, -' F D 1,;'. F " —.i —._ -:: —; —-—': —-- -----—.; —----- ------------------ iF!,:_ 2:.: i? F _; _1 I,:: 2:_:, 1 _ - F I i,:. _ -.: - F 5 i -i. -:, ------ ------- 1 -, —. -,'.:, -' F-:;- - -. -. - -.:, - ='i.'F'.:-; i',-:, 4.': "- F':B!,,, -.: - i:::"_..................'_-,!- F.':'....-F'-:; i: ~ ~ - ~

Numerical Solution of the Harmonic and Biharmonic Equations MAD Program, continued':t i i U....~1 -, I;. ~ I' s E'TiRE'F L!4.:' I F:',:':', i "'. +" I —--— c~ — rTi-iHUHUH NiL Fr I i~:-:,.:;' I'':.:',:. III'- I*Ti.4 1Hr T T-. — r7-.. ~..- -T.-T —..-.: —:-.....-.............. TF'4 E Ti Fit- iL,. OELTP LIIH~~~~E R tV:' -r, -L Cj%: - rD p!- I E0..: i..': I.WH EiiE:' E R' E R R'i C,:Y:'... -::G..... " i'..:..C.. F' II,,'.,- F'I E_:, ------— H — --------- ------- 1.~,9-E HE,r.,~ R r:- - jG-:],pi i:.'; - R]], FC:TFFE:R..T -7 E.'i:i....................... E Ii U i T is lB r~~~~~~~~~~~~~~~~~~~~:i'-:;=;-' + EE T A'RFI'F.N T F' 0 R" A T I, P" T P E:',,E'CT RF.:),,FL LJE!_:;. i,"'.IT E R "I:3H-FO E.',:.;C.CEDE D E IM:.: J'i~.:.:' L T H FI L 1 H, i H'.' I D FVP THiIIti LH Fi Y ~ - -- 1_,~' r iEi;ji I DT F;F I['T F O F"RI 0 ii - 1T FFH iF IT'' tr....- ZD7 —,,,.-" —-:;c — - --- -- 7 7 -T -f.. P. i ST E, Tq]F,-'F FILHL~ ~ ~ ~ Ft iF 5}'i... 1 E " E.L: ~PIA -1 I -- -,-. 1 --. I. _-4::, - TRi". SF ER TF FFLF: H -. 1II H E PT F.: L -'', 1:: EF:...., —:C.... —i 1, 1 E 1F —--.... —- --..,.-.,~ti iE _i,.., +I~~"l A EFFC1F~~~~~~~~dHmH ~F.' I 1,,~tl F:'T -E.,.'r ~ I ~ *- + 1i E I —---- - ---------- I- - 1 1.1.1E EF ~ ~ ~ ~~~~''"-'-F,'' -I' -~ - - "~T............ LE'-,Ti T H R 0 Li 1 H L. R -- F R F.: -.i..'; -- i~7 i~1HROU1GH LA'~., i.,.'..',:': LRPST PR I N T FF N T` U.':. 7' O I,: EF:R~~~~i~~iF;':. L. 1 - I LI'- ~ I'.. Z I. IiJ tk'i~~T PRINT" FORF1I__1T f', 1T'1 ------------------------— ~ ~ ~~~~~ I._! - F' Ri I'"T F F.:"i -T.T,!:7l EFF1Pe?~reii:7-'Z;L~ ":=-" -- ~-i T-,- r?-~ -~:.'"'"f"'T-:I- - Ti- T:iT~- "~" T - -'.'.. -,"E v ___. q' -~:-, ~'i i... i~i..,S ~.:-..T-.]-;,:'.' r i R::.':,' " T t "I:.. F I iA. L E' H, D F' F' R F.C3 i': A GR!' A typical set of computer output is shown below. E X' CE E D E D,NH01 - ":.'' 0 _-7 -R,, i.-rm; — ~ —- - - - -- - - - ----- -;~- - --- -.-.-........ —, ERR. 0 I',:-, 1:....' L - - - -- - - - - -- -- - - - - - - - - - ii i~~~~~~~~~~~';r r~~~~~~~~~~~~~c: i~~~~,: ERROR-( I E ERRO-R 0. I - 8 E-.3., i,. ". 3. 7,-, 7 E RR R. 0.,:O 5 ],=~1',- 8.5 E - S,:'S!,] f.,5],: —.:.4'-': 2. T ~~~~~~~~~~~~~~~~~~~~~~77f —----------—...0 --— ~::, -.-1 ----— F I'!:.',,..S. 31'I17::-":[~~~~~~~~~~~~~~~~~~~:'-. 18E..-.. PI,:!..,' 2:,.:. Of0 0 ERORC iI-, -IIDOF E R R Cl R i:: 1.~ ~~~~~I-:, - ~1 8 8 5 E- Z 6 E T,: i.']" —,-i'-I 0 {- {- Cl E.R.:0 R,'"1.. 1:.' = ~ -':: 5 E -.31 rF'S00:' 1, IF: =- 7000tIqi' E R. R 0 R (: V.I I -' ~{ 8 5:: E -.E -,F'-. T 7;. 2'C:'-) —- 0 "{0 0'-ER.RORT,:I I'-' E-~-~ - _- P'' i-. 0 EF..R 0ROR:. 1. 4.': - ~:3'.::5 E -'G F'S':... i::, = — 7'0'-E'.k",' UP 7.'-'1.T -- - -—....I"-'%: —....... =-' C E RRO0R C: 1 I "S, [-:' -'i-:::8E-.'6F":':I.,I':".:'-,E' 2 0 A-I -ii.. -E F:. K. 0 R C' I. 7.. I ~ 85 E-7S;, F''-I,. 7.:'- 7':C-i U C"1Cl ERR:O.:0R.:: 1.,'-::, I ~,~:..5 E -...,:,'' —. i,: I. "',::: J-UL n ~~~~~~~~~~~~ —, —~ ——'....-,,,c —--------- E F:R 0R,:. 2., l:: -~i 8. —:5 E -. ID F':;i:]2,'l )::.iI. I"1-.i!i t"-C' -- E/..: F': - — ~...........,. —---—.... -E F.: F.:0R,':.,.3:, -R ~C3.'.: -.3 F' S!, 2..J:,. 7':, _,.~:-.'EF.:F:O R C: 2.-, 4).% -..'"J1 E - I_-.. F"-:!:]2.'4"'''1.:'22:4j EF.:.'OR..,." C2`'"' -'r 2'':-!.-I..F'i,:":.,;", "':r7'.4 ~,.:.,... ~~ - - - - - - —'. — -,-' — -: -~16-~~~~~~!

Example Problem No. 97 Oase 3 (Program Orelax) This deals with the same problem as Case 1, except that arbitrary (but continuous) boundaries are now acceptable. The nmain problem now is getting the additional data on the curved boundaries read into the machine. This is done in the rather awkward but effective manner of inserting statement cards to take care of troublesome points indicated by 3 in the BOUND code. U, D, L, R data tell the distance (zero to one) of the neighboring points (mesh or boundary) above, below, left or right of the point where BOUND=3, and PSU, PSD, PSL, PSR give the values of PSI at these neighboring points. Statements must appear for all neighbors, to avoid further complications in input. Except for the other modification in the arithmetic computation of PSI at a point, everything else is as in Case 1. The problem here is the same as for Case 1, except that the area changes more gradually, in that a quarter of a circle is substituted for the right angle corner. The flow diagram for Case 1 is applicable with the indicated changes shown below in going from 1 to 2 and DELTA to ALPHA1. APl CALC DELTA BOUND(X, Y) < PSR, PSL, PSU., PSD PI The MAD program for Case 3 is shown below. W P GRAEBEL S146D 005 005 000 ORELAX $COMPILE t4AD, EXECUTE, DUMP RRELAXATION METHOD FOR 2D LAPLACE EQUATION R RDATA CARDS IN THE FOLLOWING ORDER R RXMAX(I5),YMAX(15),NMAX(IS)tERRMAX(FlO) R RPSI(0,0)...PSI(XMAX#YMAX) (8F10) R RBOUNDARIES(0,0)...BOUNDARIES(XMAX,YMAX)(8011) R(THIS LAST CARD SHOULD HAVE ZEROS FOR INTERIOR POINTS, 1 FOR RMESH POINTS WHICH LIE ON THE BOUNDARY, AND 2 FOR MESH POINTS ROUTSIDE OF THE BOUNDARY. IF THE NEIGHBOR OF AN INTERIOR RPOINT LIES OUTSIDE THE BOUNDARY, THAT POINT SHOULD HAVE A 3 RINSTEAD OF A ZERO.) R RU(XY), D(XY)g L(XY)o R(XY) * (DATA) R(L, R, U. D CARDS INDICATE DISTANCE FROM INTERIOR POINT TO RNEIGHBORS, WHEN AT LEAST ONE NEIGHBOR IS AN EXTERIOR POINT. RALL 4 CARDS MUST BE PRESENT FOR EACH SUCH POINT, EVEN IF THE RDISTANCE IS 1.) R RPSL(X,Y), PSR(XY), PSU(XY), PSD(X,Y) (MAD STATEMENTS) RPSL(X,Y), PSR(X#Y), PSU(XY), PSD(X,Y) INDICATE THE VALUES OF RPSI FOR THE POINTS DESCRIBED ON THE PREVIOUS DATA CARDS. THEY RBELONG IMMEDIATELY AFTER DELTA IN THE PROGRAM.) R DIMENSION PSL(10009 DIMB(1)) DIMENSION PSR(1000, DIMB(1)) DIMENSION PSU(1000, DIMB(1)) DIMENSION PSD(1000, DIMB(1)) DIMENSION R(1000, DIMB(1)) DIMENSION U(1000, DIMB(1)) DIMENSION D(1000, DIMB(1)) DIMENSION L(1000, DIMB(1)) DIMENSION PSI(1000, DiMB(1)) DIMENSION ERROR(1000, DIMB(1)) DIMENSION BOUND(1000. DIMB(1)), DIMB(3) -F17

Numerical Solution of the Harmonic and Biharmonic Equations MAD Program, continued INTEGER N. Co X, Y. XMAX* YMAX9 NMAX INTEGER BOUND READ FORMAT ABLE, XMAX, YMAX, NMAX, ERRMiAX VECTOR VALUES ABLE = $3159,F10.0*$ 5ECTOR VALUES DIMB(1) = 2,27,11 DIMB(2) = YMAX + 3 DIMB(3) = YMAX + 1 READ FORMAT EASY, PSI(OU)..,PSI(XMAXpYMAX) SECTOR VALUES EASY=$(8F10.6)*$ READ FORMAT FOX, BOUND(,90)...BOUND(XMAX.YMAX) VECTOR VALUES FOX=$(80I11)*$ READ DATA U. Do L. R N = 0 START X = 1 8=1 N = N + 1 C = O ALPHA WHENEVER BOUND(X,Y).NEO, TRANSFER TO DELTA PSIB = (PSI(X - 1,Y) + PSI(X + 1jy) + PSI.(XY - 1) + PSI(XY + 1))/4. ALPHA1 ERROR(XY)=.ABS.(PSIB-PSI(XY)).WHENEVER ERROPR(X.Y) G. ERRMAX, C = 1 PSI(XY) = PSIB WHENEVER N.GoNMAX, TRANSFER TO BETA X = X + 1 TRANSFER TO ALPHA BETA PRINT FORMAT INTER. NMAX VECTOR VALUES INTER = $18H0 EXCEEDED NMAX = I5,'$ 3HROUGH MID, FOR X = 01, X.G.XMAX THROUGH MID. FOR Y = O,l Y*.G.YMAX MID PRINT FORMAT ER. X, Y. ERROR(X.Y), X, Y. PSI(XY) VECTOR VALUES ER = $52, 6HERROR(I3i1HI3S,4H) = E1.4- S109 4 1HPSI(I39 1H9 I3. 4H) = F10.6*$ TRANSFER TO FINAL DELTA 6HENEVER BOUND(XY).L.39TRANSFER TO DELTA1 R RPSL, PSR, PSU, PSD CARDS HERE PSR(10,5)=PSI(11,5) PSL(109,5)=PSI(9,5) PSD( 1O,5)=PSI(10,4) PSU( 10,5)=1. PSR(11.5)=PSI(1295) PSL(11,5)=PSI(10,5) PSD(11l5)=PSI(11,4) PSU(11,5)=1. PSR(12,6)=PSI(13,6) PSL(12,6)=1. PSU(12,6)=PSI(12,7) PSD(12.6)=PSI(12.5) PSR(12,7)=PSI (13.7) PSL(12,7)=1. PSU(12.7)=PSI(12.8) PSD(12,7)=PSI(12.6) R PSIB=(PSL(XY)/L(X.Y)+PSR(X.Y)/R(XY)+PSD(X,Y)/D(X,Y)+PSU(XY 1)/U(XgY))/(l/L(XY)+I"./R(XY)+I./U(XrY)+l./D(XY)) 3RANSFER TO ALPHA1 DELTA1 WHENEVER X.L.XMAX X = X + 1 TRANSFER TO ALPHA END OF CONDITIONAL WHENEVER Y* L. YMAX 7=1 y = Y + 1 TRANSFER TO ALPHA END OF CONDITIONAL 6HENEVER C".E.1 TRANSFER TO START THROUGH LAST. FOR X = 091,X.G.XMAX 3HROUGH LAST. FOR Y = U,1,Y.G.YMAX LAST PRINT FORMAT OUT, XYPSI(XY) VECTOR VALUES OUT = $1HU.S2.4HPSI(I3,1H1IS34H) = F10.6*$ PRINT FORMAT Q. N VECTOR VALUES Q = $1HUQS224HNUMBER OF RELAXATIONS = I5'$ FINAL END OF PROGRAM -F188

Example Problem No. 97 A typical set of results obtained from the computer program is shown below. EXCEEDED Nt"I"A 70I E'Rf{..PfiR,': - O,' 0':' -., a=,=;;; —_-,. XlsO.. -. _lOll.ll -O. -"1 =8'.. f.-:.:<-I,-,.' —- 0'.....0 —..c,0000'0 ERROR 0.. }. 1 885E-36 F' I. 0..,.,1 -.OOOOO ERRO"F.!~.... i - 1 E - 400000 FRRORT 0., 3 -=. 1531 E -'I, 0,.600000 ERR- 4', 1C,::,.'.1-. E-.3 6 F'S I, 0. 4': 80000' EF.RORF C:. 1:885E-.36 FP I I, - 1. 000000 -ERRR-...i.-0_. -t —---...-6 —----.' ~.-:E-r............''................................ Il~~-.~,p~: a-~~ —.~-t-irlll 1,_~~~~~~~~~~~~~~-,,-E —- —'7 0-. —-~-,'.-l -0' ERROR 0 7,. 1 7...., 1 000"-1: ERT r l L 1 f:E 81::s F T -=. 3 P 5L E...'. 0oJ~ %5E-$6 PI'0.000000 -EI~P.OPsC1F.-('=0->l —Oi.....' —''-.....'] -—'E-' —— g...........',: —-,-j;..g —....-. -- ----— 1O 6b'r[ERR OR, C -L 1 E I - 36 I:: - 7 0 I00 II ERROR 1. -*. -. E 1 T: - 1.000000 ERRF 1 _ 1.. 4lEIu FlT 1 O, F<- fiE:, -.:. 1,:,,:,._E-.36t- P-.I:f -," =. (ll~l~~ ERRORF 1, 0': - 8'-" E-376 FIC T,E 1 rI F O: 1. -ER O'R. (~;8C........'.......-.2E; E- 0t r.........._-' ~(....... - — _............."l I9' 9 -i I ERRORC.. 2':) --- ~SE CiE-, FPSI - 1 ) -..3~91': 8. ERROR' I, 4..5'E F' Ps:i 1 4.).7 S-9.5 E~~~'UP' 1~~~~~~~~ - 1~~~~~~E-~~~~Me-I' 1. ~ Y0Y ERROFR' 1 1 o._ E a,E-': FP I' 1 3, i' 000000_0'E_ R.'R 0....1. = t....'-L'~.....F'.........E 7 1 000000ERROR( 1,:) -. 1 88E- F I,1, 7. l000000 ER RUR' -I- I..-.. E a-'1 t... —.. 7' OOOOOERROR. e, 1 ) -. 185E-E P,: 1. 10 - 7,.000000 E-F~~~~~~~~~,,-U_,~'-,C T 8`- E_: -, - - - ~.000 ERC' L —2 I- - -._ —F --- F"',,: " 7. 000000I -RR OR -' —'-'-'-,'_.:7's.SS-...6'' -' - -~' —'- - - -..- - -,- -r —-~_'~~{~C 31 ERRORC 2, 1. 1..9E.05 F'S 2, 1" 199754 -E' R.OR...'C-..2.:; — 2- I.-....... 96=E.._'c 960 ERROR 2'2 3..1RS[2E-05 PT,: 2,.)'.599598........ r- ------ - -- l"'~: —--.-r.........-. —.-.-..-.-... —,-sT"~'-....'_......i —.... —-....,'-:39 — c-'~ERPORC 4'R 1'" -. E I......... 1... ERROR( 2,' 5) - I 88tE-. - PS T, 2, 5) 1000000 ERROR' -.1 885E-.6 FPS I,: 2 7 7.0 OOOOO — FR'C- -— 1 - -. I) -''88='.... - 7-3TS 0' r IER~ROR::2ct 1...E PS!," -. 0 ERR.R. 2 9E. 2'8E:-. -,' 000000 ERROP-.. C TO-Y —-. T..-f,:, 5-E:6T- F, T-,T........T.00000U Case 4 (Program Random) An interesting variation to the relaxation type solution is the use of random walk techniques. By flipping a four-sided coin to decide which direction to move to next, the value of PSI at the starting point is given by MMAX P Z- ni PSIi PSI= MMAX Z n. (1) i=l where ni are the number of "flips" needed to reach the boundary, and PSIi is the value of PSI at the boundary point reached. The program given enters a random sequence by means of an entry point determined from the time clock. As a safeguard, if it takes more than NMAX "flips" to reach the boundary, the trial is restarted. The program is suitable for the same boundaries as Case 1. The flipping procedure is by no means as efficient as the relaxation approach, and probably is of merit only if information is needed in a small region of a very large domain. -Fl9

Numerical Solution of the Harmonic and Biharmonic Equations A flow diagram of the program for Case 4 is shown below. XNMAX, YMAX, PSI(O,O)... Boum(oo)... X=l, Y=.l, NMAX, MMAX PSI(XJAXeYMAX) BOUND(XMAXOYMAX) CALC T (T) ~~~~NUM=O T BETA GANJYIA BOUND ~Z 0 DELTA RI=RAM2B _ T T O ORI <.-25 1=1+1 3.5 RI<7 J=J+l3 F F F T T.25~RI~~.5 I=I-1 3 7~R 1 J=J-l T~~~~ N=N~~l Y=< lX., EN F~~~~ T m w g r-* S~~~UD NU = NUM +4' No PSI

Example Problem No. 97 The MAD program is shown below. RANDOM WALK FOR 20 LAPLACE EWAIION I\I.'IRMAL NUDE- IS I1NTEGLR FLCIA'IN:.G P[][N[ PSI~,.PSI, RI, R~AM2H., NUR1 DEM READ FOIPMAT Al.LE, NMAX, MMAXt XMAX, YMAX VE —CT['J. VALUES ABLE=$4I5*$ i.'I M3( = Y NiAX + 3 [.UIMB(3) Y AX + X RLA\U F EAI\1 \SY, PSi(O,U)...PSI(XMAXYMAX) VLU'foR VALUES EASY-.$(8F10.6)*$ G01.0 I';IMENSI~i', PSI ( 62.5 D I Mb (I READ FPI'A'f FOX., BOU3'JD(0 0)...BOUND(XMAXYMAX) VECFOR VALUES FOX=$(4012)*$ GO12 1..IMENS[ON hOQUN.XI(625,I;IMPM(1)), DIMB(3) V tETOR VALUES DI Pi;(I') 2,27,11 x 1 0O013' =i 6014 I I IE. (0) I =-'~435-73836/ F1(ECU1'F RA\R'2A. (T) ALPHA t1 G015'..,?M I= 00 G19 I!EL=? - G020 B3 E I'A I =IX G016 J. Y GO 17 GAMMA V,,HLNEVtI, HUU')(I IJ:).h.E. U, TRANSFER T- DELTA': I 1},A i2h. (UJ)'.,,1tENEVFR RI.bE. 0..AND. RI.L. 0.25 =I! +1 0 G024 C.'R W —E\'UVEiR RI.(,E. (;.25.AND. RI.L. 0.50 1 =1-I.... 1..iE.F.V.R RI.(i-'_NU RI-,E.-\D Ri.L. 0. 75 J=J+ 1 G028 Ur; WH:N:VER kI.u'. 0.75.AND. RI.LE. 1.00 Jd j- 1 I..:L) OF k.]NLUI I IJNAL G031 I', —-=-N + 1- G032 RO[JI I AN. U.NiNAX, TRANSFER TOD GAMMA P,[INI Fi:ORMAl i[Go,X,Y VLL['1`.. VALU-S'' — _.22'-f,0 XCLELIED NNAX FOR X=15, 711 AND Y=15*$ G036 I RANSFI'1[0 ",i [A G037 DEL (A VHL:\L 1FJLUNV)( I, Jj ). V.-2, I'RANS:F ER I [) MU.".j'=:,L" I+: *PS I (I, U.; I-i'"'i: J)[- ~.~+iJ F039 EI,C N L"V _,t V.L. [.'MAX................... r' -- - --- -- -- --....................................... -- -- --- -- -- -- --- -- -- -- --- -- -- -- iAl'Ra NFt.'A, I 0,7'1'A F.J WIF.'iF C[!xl'.I I' I{..;A L I-' I(X,Y = E['/IE'INI' U'.MA' 1 HIJCLK,, X, Y, PSI (,'.,Y'i,'PSI B VF(;I'{.:-l, VALUJES Cf!(K:.'K=,IO, S2, 4HPSI(IP,11t,1I5,2H)=FlO.6tS5,46H 1 av';cU, A': NLit.',"'"" F' T.R.WS _!U2 REA GHI UNDARY,F1.6*-.$ MU N L)(41 ~rhL~ ~ i.:~ OUN:'CXY).i'. 03 T'RANSFt}R 11) ALPH-A i —L: X.L. X'.X, IHAS!'R 0 MI,,it " 1' I Y -L. Y;'%;\X I iaf' F' AI LJ;';J L' (,...! F. _C [L I I I (-.L\1 L_ ---- ------------- ---- ---- -------- -----------—......................:1; IF P! f\ tRA, Rk,!x.; -F2Z

Numerical Solution of the Harmonic and Biharmonic Equations A typical set of computer output is shown below. PSI( 1, 1)=.334311 AVERAGE NUMBER OF THROWS 10 REACH BOUNDARY = 3.410000 PSI( 2_ 1.-:.315789 AVERAGE NUMBER OF THROWS 1O REACH BOUNDARY = 6.460000 PSI( 3, 1)=.367710 AVERAGE NUMBER OF THROWS TO REACH BOUNDARY = 6.070000 PSI( 4, 1)=.296160 AVERAGE NUMBER OF THROWS I0 REACH BOUNDARY = 5.990000 PSI( 5, 1)=.347952 AVERAGE NUMBER OF THROWS V0 REACH BOUNDARY = 8.300000 PSI( 6. 1}=.295181 AVERAG= NUMBER OF THROWS TO REACH BOUNDARY = 9.960000 PSI( 7, 1)=.300597 AVERAGE NUMBER OF THROWS TO REACH BOUNDARY = 6.700000 PSI( 6, 1)=.354658 AVERAbE NUMBER OF THROWS [O REACH BOUNDARY = 8.910000 PSI[ 9, 1)=.287466 AVERAGE NUMBER OF THROWS 10 REACH BOUNDARY: 7.340000 EXCEEDED NMAX FOR X= 10 AND Y= 1 PSI( 10, 1)=.308601 AVERAGE NUMBER OF THROWS 10 REACH BOUNDARY = 7.790000 PSI( 11 1)=.476510 AVERAGE NUMBER OF THROWS TO REACH BOUNDARY = 8.940000 PSI( 12, 1)=.1392925 AVERAGE NUMBER OF IHROWS 0I[ REACH BOUNDARY: 10.460000 PSI( 13, 1)=.425653 AVERAFE NUMBER OF THROWS 10 REACH BOUNDARY = 10.726000 PSHI 14, ])=.418621 AVERAGE NUMBER OF THROWS 10 REACH BOUNDARY = 11.60000U PS[, 15, 1)=.399716 AVERAGE NUMBER OF THROWS TO REACH hOUNDARY = 14.160000 PSI( 16, 1)=.275203 AVERAGE NUMBER OF THROWS IFJ REACH BOGUNDARY = 14.760000 PSI( 17, 1)=.361626 AVERAGE NUMBER OIF THROCWS I0 REACH BOUNDARY = 14.020000 PSI( 18, 1)=.4q05120 AVERAGE NUMBER OF THROWS TO REACH BOUNDARY = 12.896000 PSI( 19, ] )=.23818.3 AVERAGE NUMBER (OF THRUWS 10 REACH BOUNDARY 12.44U000 PSI1 20, 1)=.383511 AVERAGE NUMBER OIF THROWS TO REACH BOUNDARY = 10.310000 PSI( 21, 1)=.445896 AVERAGF NUMBER F)F THROWS TO REACH BOUNDARY = 14.010000 PSI( 22, 1)=.260944 AVERAGE NUMBER OF THRO(WS 1'O REACH BOUNDAPY: 6.99U000 PSI( 23, 1)=.265636 AVERAE NUMBER OF IHROWS T') RFACH BOUNDARY = 5.820000 Case 3 has been used successfully in the Photoelasticity course for solving f or the first stress invariant to enable separation of the stresses. The instructor's program was used by all students. None of the other programs have been class tested. -F22

Example Problem No. 98 JOUKOWSKI AIRFOIL by William P. Graebel Department of Engineering Mechanics The University of Michigan Course: Mechanics of Inviscid Fluids I Credit Hours: 3 Level: Senior This is a problem whioh appears in all of the books dealing with classical hydrodynamics, but has become of perhaps decreasing interest with the advent of high-speed aerodynamics. The approach used in such books as Milne-Thomson (Theoretical Hydrodynamics, MacMillan, 1960) and Streeter (Fluid Dynamics, McGraw-Hill, 1948)., emphasizes the graphical method. The advantage accrued by use of the computer is to emphasize the fundamentals of the problem. Basically., this is the problem of the mapping of a circle into a shape with one cusp point. For potential flow with circulation around a circular cylinder in the complex plane., the complex potential is given by =U(3 e' 1A+ a2 eia)+ iP In~ 1 where U and 0, are the free stream speed and direction, a is the radius of the cylinder, and r is the magnitude of the circulation. Transforming this into the complex z plane by means of z + am I e + a b /(5 +am e'~ (2) with b2 +2mb cos G +m2 - 1=0,. (.3) the airfoil is obtained. For the cylinder of radius a,. a e on the cylinder, hence z/a = (cos G + m os6 ~) (1 + b2/'R) + i(sin G + m sins) (1 b2 /R)., (4) R = (cos G + m cos ~)2 + (sin G + m sinS6)2 (5) The parameters governing the shape of the airfoil are then m and.using equations (1), (2)., and (3), the streamlines in the z plane could also be found. In the present program, we have contented ourselves with plotting.jL~t the airfoils. A given number (MAX) of G's are used to carry out the computations given by equation (4), computing the airfoil shape at points 2w/MAX radians apart on the circle. For given m, S. the values of these points are given numerically and then graphically., so that the computer can giv3- the shapes of many airfoils in less time than the best draftsman can construct a single airfoil. (The plots are scaled by a factor L = 2 -(b + m cos G), which estimates the chord length of the airfoil.) The effect of variation of m and is thus easily grasped.

Joukowski Airfoil A 1is t of variables appearing in the program but not described in the text is given below. MAX Number of points on airfoil to be computed DELTA(I) The angle delta in degrees DELTA The angle delta in radians DELMAX Number of DELTA's to be read in M(I) M in the transformation (0 O1 MMVAX Number of M's to be read in B Positive root of equation (3) XI, ETA Coordinates in zeta plane XV,~ Y'T Coordinates of points on airfoil in Z plane A flow diagram for the program is shown below. EXCUTE i MAX so ~~~~~~~~~1LMXPLOT 1 DE1LTA(l)...DELTA(DELMVAX) M(I)....MOYM )(0.,3,15,6,15) FIN IN ~~~~~DELTA=DELTA(P)*Tw/180 MAS (A) *C A:> MMVAX P>MV]XS =SIN DELTA PB=-A2 ~ 1 —_(-M AS)~2 THETA=2iTN/18o XI=CO:S (THETA)~MAC L=-2*(B DELTA(P) B=B2 FNIOSH ETA==SIN (THETA)3 + MAC) M(A),L BS B N=O~~~~~l + MAS R = (4~)2 + (ETA)2 XSO=BSXI XV(N)X XY FNS {EXECUTE PLOT 2 10 =XI(1BSR)/,L j_____________________ -151O. Y =ETA(1-BSOR)/L ~~~~~~~~~~-0.5) 4 EXERNCUTE PLOT 3 EXECUTE PLOT 4 LABEL=:Y/LFI ($ —$.,XV.,YV.,MAX) (3., LABEL)

Example Problem No. 98 The MAD program is shown below. JOUKOWSKI AIRFOIL READ FORMAlT DATA1,MAX READ FORMAT DATA2,DELMAX,DELTA(l)...DELTA(DELMAX) READ FORMAl DATA2,MMAXM(l)...M(MMAX) INTEGER MAX t )ELMAX, MMAXA, P,N VECTOR VALUES DATAI=$I2*$ VECTOR VALUES DArA2=$12, 15F5.5*$ ) I MENS OlN DEL A (_1 5)_ M(15 ) DIMENSION XV(100), YV(100) D I MENS ON_.. IMLA E_(2000.... VECT(OR VALUES PI=3.14159 VECTOR VALUES rTWOPI=6.28318.............__-___ EXECUTE PLOTI. (0,4,15,6,15) P 10180=P I / 180. X H JH FIN FOR A=1 1.tA.G.MMAX THROUGH FIN,FOR P=1,1,P.G.MMAX DELTA=I)ELTA( P )*PIO180 ____ C=COS. ( DELTA(P) ) __S = S. I S=SIN.( DEL-1'A (_P )') _- ______-_-____-__________-___-_ __-_-____,AC=M ( A)*C MAS=M (A )*S__ __ _ P=-MAC +SQRT. (1.- MAS*MAS) L=2.*(B + MAC) PRINT F3RMAT TITLEDELTA(P),M(A),L VECTOR VALUES TIlLE=$13H FOR DELLA = FO10.5,16H DEGREES M/A 1 F10.5,51H THE FOLLOL ING ARE OBTAINED FOR X/L AND Y/L. L/A 2F 10.7*$ BROUGH FINISHFOR N=1N.E.MAX ITHETA=N*TPOMAX XI=COS. (THETA)+MAC FFA=SIN. (THETA)+PsAS R=X I*XI +ETA*ETA ________________ B SOR = S, R X=XI*( l.+BSOR)/L Y=ETA*(l1.-BSOR)/L XV(N) X.YV(N) =Y FINISH PRINT FORMAT RESULTXY ___ _ _______ EXECUTE PLOT2. (IMAGE, 1.5, -1.5, 1.0, -1.0) PRINT COMMENF $1$ E-XECUTE PLOT3. ($+$, XV, YV, MAX) EXECUTE PLO'14. (3, LABEL) ~PRI1NTh-CO1MMEN Il'i VECTOR VALUES LAbEL = ___ __Y_______L_ _____ FIN CONTINUE VECTOR VALUES RESULT =$1H S55F10.6,5S5,F1O.6*$ END OF PROGRAM --- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - F 2 —- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Joukowski Airfoil A typical set of computer output is shown below. FOR DELTA = 15.00000 DEGREES, N/A =.20000 THE FOLLOWING ARE'OBTAINED FOR X/L AND Y/L. L/A = 1.9'830129 1. 1872.56 -. 050904 1.036732 -. 126794+.708810 -. 137762. 487903 -. 104112.1lb5991 -. 052068 -.153671 -.003221 -.452478.029131! -.714320.040291 -.924972.033028 -1.072219.016012 - 1. 146066.001!71! -1. 138655.004253 -1.043833.03'5555 -.856800.101641 -. 575282.. 196586........0~o507?6.2.......~''9~ 677.........................................227432.360709___ -.661245. 345316 1.006939.239799__ 1.188652.085704 Y LI I I I I I I I I I I I I I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~... I I I I I I I I I I 1 1 I I I I I I I I I I ~ ~ ~~~I! I I I I I ~ ~ ~~~I I I. I II I I I I I 1 I 0...... * 500....................... ~......+.................+.....................~.........~........I....~................... I I I I I I I I I I I I- -I.. I I I I I I I I I I I I 1I -~~ - - I-I - - +. 1~{ I I +I I I I I....I 9oI.........................I...._..........I._..............I.....t..............,_Ii I I I I I I............................. +........... I I +I I I I I I I I! i +I+ + 1+ 1! I I _ i I I I +1I + II i!I I I _ I I I...... _!...............I..... I I..........+ _ 1+ I!I + I I_ I I I I I i I I II ____ I ________ I I I ~~~~~~~~~~~~~~~~~~~~I, i I I I I I! I... I +.......I I I I +I I I I I I I I I I I I I I'I I ~~~I I I +I I +I I - I I -I I +I+ I I I I I _I I _I I --- I i I I I I I I _ _ _ _ _ I _ _ _ _ _ I _ _I I _ _ IIIII I I I I I I I I I I I I I I I I I I I I I I I I I I i I I I I *I I.. 1 I I I I I I I I I I -100 + ~ + ~ + ~ +~+~+ ~ -.500...... -.000.............. -.500............ 0.000............... 0.500...... 1.000.. i.5...........0..... I I I I I I I~~~~F26

Example Problem No.- 99 PRINCIPAL AXES OF A SECOND ORDER TENSOR by William P. Graebel Department of Engineering Mechanics The University of Michigan The problem of determining principal axes occurs frequently in mechanics in connection with the stress., strain, and moment of inertia tensors. Using stress as an example, the law of transf ormation when axes are rotated is IT ~~~k 1~ i'j = a~a F (1) k=ljf=l where the a's represent direction cosines. As a consequence of the orthogonality of axes., Z ai O~1if i~j (2) k=l ik lf~ The problem is then to find a set of a's for which f =0 when i/Lj. (This problem is the same as the problem of diagonalizing a matrix.) This is done for given Tk. by solving 3 3 k and any 3 of equation (2). In the attached program, besides finding the a's,. the T are also solved for as well as the three principal invariants 3 ONE =Z -r.(4) i=l 2 l1 3 TWO = 0.5 (ONE) 2 z >z'Tij.. i i=l j=l u i THREE= det T.. The input data is XX = IXY = =' XZ = T xx ~ xy yxxz IzxYY =,YZ ='T PZ Z = T yy yz zy zz read in that order on two data cards with the format 3E14.8. The invariants, Tir and a1)s are then found in a straightforward manner. r- is found as roots of the equation 3- (ONE) ~V2 (TWO)' (THREE) = 0 (5) As a check, the computed a~'s are then put into the remaining 6 of equation (2),, and should

Principal Axes of a Second Order Tensor The program can be easily adapted to show the need for carrying sufficient figures. In the THROUGH statement ending at ROOT, the value of one of the roots of equation (5) is found by Newton's method to a greater nmumber of significant places than is required for engineering work. If, however, the number of significant figures retained were only two or three, the check will show that the resulting principal axes are far from being orthogonal. This "extra" precision then is important for the computations which follow it. This problem is suitable for introduction at the 412-422-542 level. It has not been used in the classroom. A list of program variables for quantities not described in the text is given below. AXX, AXY, AXZ, AYX, AYY, AYZ, Direction cosines of principal axes. AZX, AZY, AZZ TWOHAF Second moment of the tensor. Il, I2, I3 Three roots of cubic equation (5). I, IN Temporary names for Il. CHECK 1,...,CHECK 6 Accuracy check for orthogonal axes. These should all be zero. TEST Switch to stop iteration on Il A flow diagram for the program is shown below. XXXX XY., 2Z T XX XZ. XZ, YY., CAL ON W=~ OE 1 CALC THREE 1 m31-THIS lT/L ST_ _ CALC Il by Method C12 = 0.5-I II3C12 CAL BXv CA B, BZGX, C12 0.5 C1 BY G~~~~~~~~j AXZ.,AYX., AYY., 0C2 = TWO-Ii*O1I(l2)_C GY_,GZ AXZ. L AYZ,AZX,AZY, CALC CHECK1 Il, I2, ONE. TWO, AXX300C IH13 THREEAZCHK6 -F28

Example Problem No. 99 The MAD program is shown below. PRINC(IPAL ALLIJES AND I)IRECTIUNS OF A SECO(ND ORDER SYM. METRIC I EL S SU I J[EtJLP N, ILSr It- ='r XX+YY+ZZ T'-')f'iAF i-J-. 5* ( XX*XX+YY YY+ZZ*ZZ ) +XY*XY+XZ*XZ+YZ*YZ I =A. 5 ~*ON4E*IE-r OHAF I HRE-=. < ( iYY*ZZ-YZ*YZ ) fXY*(-ZZ*XY+YZ*XZ)+XL*(XY*YZ-YY*XZ) 1EST=U 1 JHENF,1J U-' [ FUR N= 1 i TEST. E ~ 1 L= NI-( I) l F C ( IO+ I r)-I'HEE)/I(3. I-2.0O E) )N L WHENLVF.'. AS. ( 1.0-I\./I).L..00000003, TESU'=1 ROOI I=IN C2= TWO- I *C. 1 (? =S< T. ( C1 1 2- C 1c-2 )! = 12 0 Rl I 3=C12-R L:X X= (XZ. XY+YZL-(-XX11 1 1)/(XY*YZ-XZ.(YY-I i)) FY-(XZ*XY+YZ*. -X+I2) )/(XY*YZ-XZ*(YY-I2)) Z=(XZ XY+YZ*(-XX+I3) )/(XY*YL-XZ(YY-I3) ) GX= ((XX- 1() * YY-1 1)-XY*XY)/(XY*YL-XZ*(YY-I1 )) LY= ( ( X- 12 ) * (YY- I 2 ) -~Y*XY ) / ( X YY*YZ-X Z* ( YY- 12 ) ) (( X- I 3 ) * ( YY-13 )-XY*XY ) / ( XY*Y-XZ* (YY-13 ) AXX= 1./SO(RT. ( 1.+HX BX+GX*X X) AYX=I.O(/SQRT.(1.(';+BY*tY+GY*GY) AZX=1. /SUR T.(1 I.+Z-HZHZ+bZ*GZ) AXY = IX*AXX A YY = Y' AYX AZY = FZ*ALX AXZ = (X*AXX AYZ = (',Y*AYX AZZ=UZ*AZX CHECKI = AXX*AYX + AXY*AYY + AXZ*AYZ HLECK2 = AXX*/\ZX + AXY*AZY + AXZ*AZZ (HECK3 = AYX-*AZX + AYY*AZY + AYZ*AZZ LHECK4 = AXX*AXY + AYX*AYY t AZX*AZY U-C —--- = A XX AX - P X-* —-- -----— X-ZZ LH-CK6 = AXY*AXZ + AYY*AYZ + AZY*AZZ PR IN I F.ORMAT 0UTi11I,1t2,I3 IRIN I i \ ORNMAIUTL'F1r,0N:,TWO~ 1'HR E L' RIt I vMl A T t2J 2 R INI f:ORMAl J0UT3,AXX,AXY,AXZ PRIN F'RMAT CIUT3,AY,, AYY,AYZ PRiN[ FJRMA'T t]UT3,AZ,<jAZYAZZ PP INT FORMAT ou r iI N1 F nRMAT 0TUT3, CI- E-CK, CtHECK2, CHECK3 P-1 IN IFRMAT IJOUT3, CI-[CK4, CHECKS, CHFCK6 -.-uTUR VIALUOS IJ) A $'A=tIOH,2,5i-_ XX -= F14.8,6-I, XY = E14.8t6H, X 17 = E14. 28,,1Hi YY = Ei4.8,6-t, YZ = E14.8,6H, ZZ = 614.8*$ VLCTOR VALUEfS 0)U1~1=11H,',t42Ht1FE PRINCIPAL VALUES ARE, IN X, 1YZ [:RTE, FZ-)RlR F4.8,2 t —El14.-8-2H t, C14.8*$ ___ _ VE TLUR VALUES UlJ1R = I -----— 1t.i(,S2,'t5HtH-'iE TiREE INVARIANTS ARE, IN 1, 2, 2 3, R,)ER.J(.14. 8,t1, ) *$ VECT("'R VAL()':,S UJI 2=$LHO, 52, 0HTHlIE I.JlRECTION COSIiNES AXX, AXY, 1 AXL,.Y., AYY, AYZ, A' X,, 7Y, ALL ARE',

Principal Axes of a Second Order Tensor A typical set of computer output is shown below. XX =.10880952E 02, XY =.2571428u)E 01, XZ =.42857140E 00, YY =.57380951E 01, YZ = o.85714280E 00, ZZ =.13738095E 02 OE THE PRINCIPAL VALUES ARE, I10 X,Y,Z ORDER,.11696239E 02,.14031550E 02,.46293523E 01 THE THREE INVARIANTS ARE, I4 1t 2, 3, ORDER.30351143E 02,.28321937E 03.15975236E 039 THE DIRECTION CO(USINES AXX, AXY, AXZ, AYX, AYY, AYZ AZX, AZY, ALL A —RE. 88412448E _0o____t._3466290E OROt)_ -. 32607464E 00.27816949E 00,).18368832E 00,.9428044')E 00.37541855E 0() -. 92426053E 00,.69306377E-01 THE CHECK ON OUR'TH),O'U"4ALITY OF ITHE PRINCIPAL AXES IS -.14510006E-04,.19116099E-05, -. 35120174E-05 -.44070184E-otq, -.12225104E-04, -.60535961E-0( -F30