THE USE OF COMPUTERS IN MATERIALS AND METALLURGICAL ENGINEERING EDUCATION M. J. Sinnott and L. H. Van Vlack Department of Chemical and Metallurgical Engineering The University of Michigan September 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.

ABSTRACT During the past two years faculty members of the Materials and Metallurgical Engineering curricula have been taking an active part in computing work at The University of Michigan. All undergraduate students in these curricula are now required to take an introductory sophomore-level computer course taught by Computing Center and Department of Mathematics personnel. Digital and analog computer work has been assigned in seven departmental courses during the past year, giving students an opportunity to gain practice in the application of computers to the solutions of their engineering problems. A description of the curricula and a discussion of computer use in the classroom are included. This report also contains a selected set of nine example problems appropriate for classroom use prepared by departmental faculty. These may be considered as a supplement to the 87 example engineering problems, including several related to the materials and metallurgical subject area, which have been published previously by the Project. -El -

Table of Contents Page I. Introduction E3 II. The Student and Curricula E3 III. Selecting Engineering Problems for Computer Solution E5 IV. Discussion of Computer Use E6 V. Conclusions E7 VI. Example Problems E7 88 Furnace Efficiency as a Function of Stack Gas Temperature E9 89 Mass or Heat Transfer by Diffusion E15 90 Predicting the Scrap Requirement for the Oxygen Steel Converting Process E19 91 Precision Lattice Parameter Determination for a Cubic Material E30 92 Ionic Crystal Structure E40 93 Unsteady-State Heat Conduction in Solidifying Alloy E51 94 Cooling of Pig Iron in a Transfer Ladle E58 95 Digital Computer Analysis of Heat Flow and Temperature Distribution Around a Copper Converter Tuyere E64 96 Digital Computer Analysis of Galvanic Cell Data E71 VII. References Tables IE Curricula Requirements IIE List of Example Problems Figures IE Sequences for Stem Courses Utilizing Computer Problems E2

USE OF COMPUTERS IN MATERIALS AND METALLURGICAL ENGINEERING EDUCATION I. INTRODUCTION This report presents a set of problems which are suitable as assignments in undergraduate curricula in metallurgical and/or materials engineering. The curricula are outlined briefly to indicate bhe student's background and course sequences. A discussion of the computer experience in these undergraduate engineering curricula is included, citing both the advantages and limitations. The work covered in this report is a direct result of studies sponsored by The Ford Foundation at The University of Michigan on "The Use of Computers in Undergraduate Engineering Education." (See References 1, 2, 3). II. THE STUDENT AND CURRICULA During the past several years a course has been developed in the Mathematics Department at The University of Michigan which introduces the sophomore student to the principles of computer programming. This one credit hour course (Math 373), which is described in more detail in Reference 4, presents the introductory elements of computer organization and acquaints the student with programming procedures through the use of a formal computer language, MAD (Michigan Algorithm Decoder) (See Reference 5). This language is particularly suited for use at a university computing installation. Its use is relatively widespread and it is sufficiently similar to other available languages (e.g., FORTRAN, ALGOL, etc.), that the student has no major difficulty reorienting himself to other commonly used languages and systems which he may encounter in industry. In the introductory computer course the student solves four problems on the computer. These are normally related to numerical topics such as 1.) finding the roots of polynomial and transcendental equations, 2.) interpolation using Lagrange's or Newton's formula, 3.) GaussJordan reduction of a series of simultaneous linear equations, and 4.) least squares polynomial fitting. The above problems obviously do not make the student an expert. However, they do provide him with the basic computer training required to allow him to handle the computer problems in stoichiometry, thermodynamics, and rate processes which will be described later. Thus, the student is quickly in a position to use the computer on problems he would be required to solve in his normal course of studies, and a minimum of time is required in computer instruction for the computer's sake. The requirements for the two curricula, metallurgical engineering and materials engineering are given in Table IE. As might be expected, these two curricula are very similar, particularly in their main stems, or 'principles" courses. Figure IEa presents the processing stem of the common part of the two curricula, with an indication of the approximate semester for the scheduling of each course. Although the course label may not be identical with offerings in E3

other schools, the general goals are comparable. The "Introduction to Engineering Calculations," is primarily material and energy balances. "Thermodynamics" provides an introduction to the several laws with emphasis on property relationships. "Rate Processes" is a unified study of heat, mass, and energy transfer, and introduces chemical kinetics. The latter two courses have integrated laboratories. The materials stem parallels the processing stem as shown in Figure IEb. The "Structure of Solids course follows physical chemistry and is basically a solid state physics course designed for engineers. The balance of the materials courses are described sufficiently well by the indicated course titles. These are taken during the final year after the student has studied the full background of fundamentals. TABLE IE Curricula Requirements B.S.E. (Materials Engineer) Credit B.S.E. (Metallurgical Engineer) Credit Courses Common to All Engineers: 43 Courses Common to All Engineers: 43 Computing Techniques* 1 Computing Techniques* 1 Advanced Chemistry 17 or 18 Advanced Chemistry 15 Advanced Science 4 or 5 Advanced Science 7 Engineering Cognates 14 Engineering Cognates 14 (Including Mechanics and Circuits) Humanities, Social Sciences, Humanities, Arts, etc. 14 Arts, etc. 16 Economics 6 Thermodynamics and Rate Processes Thermodynamics and Rate Processes Introduction to Engr. Calculations* 3 Introduction to Engr. Calculations* 3 Thermodynamics* 5 Thermodynamics* 5 Rate Processes* 5 Rate Processes* 5 Electives 2 Separations 3 Materials Materials Structure of Solids* 3 Structure of Solids* 3 X-Ray Diffraction* 3 X-Ray Diffraction* 3 Physical Metallurgy* 7 Physical Metallurgy* 7 Physical Ceramics* 4 Polymers 4 Metals Processing Engineering Design 6 Processing of Cast Metal 2 Process Design* 4 Electives 3 Courses in which computer problems have been introduced. E4

Use of Computers in Materials and Metallurgical Engineering Education Figure IE Sequences for Stem Courses Utilizing Computer Problems (* indicates computer problems have been introduced) (a) (b) Processing Stem Approximate Semester Materials Stem Introduction to Engin- Physical Chemistry I eering Calculations* Computer Techniques* Computer Techniques* 3 Thermodynamics* 4 Physical Chemistry II Rate Processes* 5 or 6 Structure of Solids* 7 X-Ray Diffraction* Physical Metallurgy I* Metallurgical Process Physical Ceramics* (not required Design* (Not required for for Metallurgical Engineering) Materials Engineering) 8 Physical Metallurgy II* III. SELECTING ENGINEERING PROBLEMS FOR COMPUTER SOLUTION It seems appropriate to select computer oriented problems for materials and metallurgical engineering courses which enhance training in two areas primarily: 1.) understanding engineering principles, and 2.) implementing design optimization. The example problems included in section VI represent such a selection with the above goals in mind and result in the development of two major categories of problems. The first category includes relatively simple problems requiring for their solutions either manual or computer procedures which are repetitive in character. Admittedly, in some cases, these problems may require more 'time to program and completely check out on the computer than to solve manually by traditional procedures. However, the computer requires a sufficiently greater accuracy of thinking on the part of the student so that only one problem may need to be assigned on a topic rather than two or more. The calculation of furnace efficiency (Problem 88) is an example of this type of problem. The second category of computer usage does not involve the programming process but does require the student's knowledge of programming. This is illustrated in the diffusion problem (Problem 89) which is preprogrammed for the class, but requires the student to provide appropriate data for the correct solution. Thus, if the student is to determine carbon diffusion into steel he must define suitable starting and terminal conditions. The computer then makes the calculation E5

and supplies a more complete answer than could be demanded otherwise, with the consequence of better teaching illustrations. This particular problem has the further advantage that the dual use emphasizes the importance of the essential gradients and mobilities upon the unsteady-state transfer, whether it be energetic or atomic in character. This second category of preprogrammed computer usage in class may be incorporated into a course for another purpose as illustrated by the problem on ionic crystal structure (Problem 92). The purpose in assigning this problem is to have the student gain a better understanding of the several simpler ionic structures. More specifically, the student must understand these structures if he is to provide the appropriate data for the calculations. Following this, the computer calculates the Madelung Constant which provides a positive check on the student's work. This particular problem is also programmed so that this student may analyze his error if he receives an incorrect answer. IV. DISCUSSION OF COMPUTER USE This discussion assumes that the instructor has taken time to develop a familiarity with the computer so that he may be able to effectively use the computer himself. The use of computers at the undergraduate level is not recommended unless this familiarity on the part of the instructor exists. Beyond the above limitation there is one major disadvantage and there are four advantages which can be cited for the use of computers at the undergraduate level. The disadvantage of computer usage is the tendency on the part of the instructor to let the computer monopolize the time of the student in the class. Since the computer often becomes fascinating to both instructor and student, and since the number of problems suitable for solution using the computer is great, the instructor must "lean over backwards" to keep a balance in his course. The advantages of computer usage in the areas of materials and metallurgical engineering include the following: 1. There is an opportunity to acquaint the future engineer with a tool which certainly will receive wide usage during the course of his professional life. It is difficult to imagine that this basic knowledge will not be desirable in either an engineering or management position of tomorrow. 2. Computer programming requires the student to synthesize a problem solution as-well as to analyze the situation. Therefore, there are advantages in the incorporation of computeroriented problems so that the student may experience more training in synthesis skills which are unique to engineering curricula. The student soon realizes that his solution must be a rigorous one. An attempt to "second guess" the professor usually fails at this point. E6

Use of Computers in Materials and Metallurgical Engineering Education 3. Simple but time consuming optimization problems can be incorporated into design courses. This particular advantage is most apparent in the processing courses within the metallurgical and materials engineering curricula. However, it is anticipated that this aspect will become more significant in the future in other courses as well. 4. Through the use of problems whose solutions have been preprogrammed, the student may be required to consider problems which are impossible in normal assignments. V. CONCLUSIONS The advantages of introducing the computer into the undergraduate engineering curricula outweigh the disadvantages, particularly when consideration is given to the fact that today's student will be serving his professional life in a technical and management situation where the computer will be even more important than it is today. A major advantage is the increase in synthesis aspects of the engineering training. Finally, time consuming and repetitive problems which have been bypassed in the past may be used. The use of the computer is not recommended unless the instructor familiarizes himself with the computer's potentialities and limitations. This familiarity implies that the instructor can program his own problems, analyze his own errors, etc. Care must be taken by the instructor so that his enthusiasm for computers does not permit an imbalance in the courses where computers are used. VI. EXAMPLE PROBLEMS Since the beginning of the Project on the Use of Computers in Engineering Education, a number of participants have prepared solutions for example problems of their own choosing which they felt would be suitable for use in a classroom at the graduate or undergraduate level. Nine complete solutions of problems which have been used in undergraduate classes are included in this section. They are listed in Table IIE. All of the digital computer programs which were required in the solution of the problems, were programmed in the MAD (Michigan Algorithm Decoder) language, which is described in a number of places, including Reference 5. E7

TABLE IIE List of Example Problems Number* Title Author Page 88 Furnace Efficiency as a Function of Stack B. Carnahan E9 Gas Temperature 89 Mass or Heat Transfer by Diffusion M. J. Sinnott E15 90 Predicting the Scrap Requirement for the R. D. Pehlke E19 Oxygen-Steel Converting Process 91 Precision Lattice Parameter Determination J. V. Gluck E30 for a Cubic Material 92 Ionic Crystal Structure L. H. Van Vlack E40 93 Unsteady State Heat Conduction in a J. R. Street and E51 Solidifying Alloy J. 0. Wilkes 94 Cooling of Pig Iron in a Transfer Ladle R. D. Pehlke E58 95 Digital Computer Analysis of Heat Flow R. D. Pehlke E64 and Temperature Distribution Around a Copper Converter Tuyere 96 Digital Computer Analysis of Galvanic R. D. Pehlke and E71 Cell Data K. J. Guion These problems may be considered as a supplement to problems 1 through 87 published in previous reports of the Project. In addition to the problems presented here, several others related to materials or metallurgical engineering subject areas are included among the 87 problems previously published by the Project. In particular, see Problems 1, 2, 3, 4, 5, 13, 14, and 17 in Reference 1, and Problem 47 in Reference 2. The authors wish to acknowledge the contributions of several engineering faculty members in preparation of these problems. Several classes of students are to be thanked as well, particularly for their tolerance and interest during this early phase of computer integration into engineering course work. Finally the authors wish to thank the Project for its support and assistance in providing the necessary computer experience and also for help in preparing and reproducing this manuscript. E8

Example Problem No. 88 FURNACE EFFICIENCY AS A FUNCTION OF STACK GAS TEMPERATURE by Brice Carnahan Department of Chemical and Metallurgical Engineering The University of Michigan Course: Introduction to Engineering Calculations Credit hours: 3 Level: Sophomore Statement of Problem The burner of a domestic gas furnace is adjusted to use 20% excess air over that theoretically required for complete combustion. Write a MAD program which will prepare a table relating the efficiency of the furnace to the stack gas temperature at 200 degree intervals beginning at 700F. It may be assumed that the natural gas consists of pure methane (CH4), that the inlet air of 50% relative humidity has a temperature of 700F, and that reaction is complete at all temperatures. Solution Combustion proceeds according to the reaction CH4 + 202 -- C C02 + 2H20 (g) + QR (1) or that 1 mole of methane reacts with 2 moles of oxygen to yield 1 mole of carbon dioxide and 2 moles of water vapor with a release of heat QR. The percent efficiency can be calculated as eff. (%) = (QF/QR)X(1OO) where QF is the amount of heat absorbed by the furnace per mole of methane burned and is equal to the total amount available,QR,minus the heat content of all product gases relative to the inlet temperature. If we denote the heat carried away by the product gases as Qp, then the efficiency can be written as %R - % eff.(%) = " Q X(100) where (2) QR = z r 1 p i= ( At) ] mi where - (3) is the average specific heat of the ith component gas in the product mixture of N gases, mi is the number of moles of component i, and At=(T-TZ) where T is the stack gas temperature and TZ is the inlet feed gas temperature. Since the specific heat for the product gases vary as functions of temperature, Qp may be written more accurately as t=T Qp = z X cp (t) dt m (4) t=TZ E9

Furnace Efficiency as a Function of Stack Gas Temperature A mass balance for the system based on one lb.-mole of methane is tabulated below where PCXS is the percent excess air in the entering mixture, and Wa is the lbs. H20/lb. dry air. Component Into System Out of System CH4 1.0 moles 0.0 moles 02 2.0(100+PCXS)/100 2.0(PCXS/100) N2 2.0(100+PCXS) (79/21)/100 2.0(100+PCXS)(79/21)/100 H20 Wa(100+PCXS) 2.0*29.0/( 18.*100) 2.0+H20 into system C02 0.0 1.0 Data required for the solution of the problem is as follows: 1.) QR' the heat of combustion at 70~F. 2.) c, specific heat as a function of temperature for the four components in the stack gases. 3.) Wa, the pounds of water per pound of dry air for the indicated relative humidity of the inlet air. The standard heat of combustion for methane is found to be 1 51,571.4 Btu./lb. or 344,032 Btu./lb.- mole of methane at 25~C (77~F). Assume that the heat of combustion at 70~F is the same as at 770F. The specific heats of the four product gases as functions of temperature are as follows:2 0,2 cp = 8.27 + 0.000258t - 187700/t2 N2, Cp = 6.50 + 0.00100t H20, c = 8.22 + 0.00015t + 0.00000134t2 C02, c = 10.34 + 0.00274t - 195500/t2 where t is in ~K, and c is in cal./ (eK gm-mole.) p W for an air temperature of 70 F and relative humidity of 50% is given as 0.017 lb. water/l dry air.3 Calculation of Q In order to permit iterative calculation of Qp,a general cp equation of the form c = ai + bt + cit2 + di/t2 p. i 1 1 is used. Thus the coefficients for 4 product gases are as follows: Comp. i a b c d 02 1 8.27 0.000258 0.0 -187700 N2 2 6.50 0.00100 0.0 0.0 H20 3 8.22 0.00015 0.00000134 0.0 C02 4 10.34 0.00274 0.0 -195500 1. Perry, John H., Chemical Engineers Handbook, McGraw Hill, Inc., New York; 3rd edition, 1950, p. 244. 2. Ibid.; pp. 210-212. 3. Brown, G.G., Unit Operations, Wiley and Sons, New York, 1953, P. 545. E10

Example Problem No. 88 The integration required in Equation (4) can then be determined as t=T t=T cP(t) dt = (ai+bit + cit2 + di/t2)dt (5) 1 <y t=TZ t=TZ = ai (T-TZ)+bi (T2_TZ2)/2+ci (T-TZ3)/3+di (T - ) In MAD notation this would be written as A(I)*(T-TZ)+B(I)*(T.P.2-TZ.P.2)/2.+C(I)*(T.P.3-TZ.P.3)/3+D(I)*(l./TZ-1./t) (6) Since the constants given above are for temperature in ~K, T and TZ must be in 0K. The result of the integration will then be in units of cal./(gm7mol.~K) which can be converted to English units (BTU/ lb.- mol. OR) by multiplying by 1.8. If we then let mi in Equation (4) be represented by MOLOUT(I) in MAD notation, the total heat (in BTU's) carried away by component i is given by (A(I) *(T-TZ)+B(I)....etc. ) 1.8*MOLOUT(I) (7) Evaluating this expression for all values of 1(1=1,2,3,4) yields the contribution of 02, N2, H20, and CO2 respectively to the sum in Equation (4). These contributions can be added into a variable location (HTCAP) as they are calculated, so that after the summing operation is completed HTCAP will contain Qp, the total heat carried away by all gases leaving the stack. A flow diagram for the solution of the problem and a listing of the MAD program are given on the following pages. The MAD program listing is discussed below. Variable names in the program are as follows: A Array of 4 constants for the specific heat equation, i.e., al,a2,a3,a4 B Array of 4 constants for the specific heat equation, i.e., bl,b2,b3,b4 C Array of 4 constants for the specific heat equation, i.e., c-,c2,c3,c4 D Array of 4 constants for the specific heat equation, i.e., dl,,d,d3,d4 MOLOUT Array containing values of mi (Equation (4)) PCXS Percent excess air DELT Increment in ~F. for preparing the table of stack gas temperature vs. efficiency TZ Inlet feed gas temperature in ~K. TF Stack gas temperature in OF. T Stack gas temperature in ~K. HTCAP Variable which is used to keep the partial sum in Equation (4) The values of A(l), A(2), A(O), and A(4) are preset in memory by means of the VECTOR VALUES A(1) = statement. Ell

Furnace Efficiency as a Function of Stack Gas Temperature The executable portion of the program begins with the statement labeled START. The first statement causes one data card to be read containing the values of PCXS and DELT which is the increment in temperature to be used for the tabulation of efficiencies. The four statements which follow calculate the moles of the various components per mole of methane in the product mixture. MOLOUT(l), (2), (3), (4) are the moles of the various components O2, N2, H20 and C02. The PRINT statement which follows merely prints out the value of PCXS and the number of moles of the four components in the stack gases. TZ is the inlet temperature, 70~F, in OK. The iteration loop starting with THROUGH LAST and ending with the statement labeled LAST is the loop which prepares the table. TF is the temperature (in OF) which begins at 70 and is incremented by an amount DELT each time until TF >5000. The first statement inside the iteration loop determines T which is TF in ~K. The statement HTCAP=O simply clears the storage location where the partial sum to determine Qp will be saved. Qp is calculated in the 2 statements starting with THROUGH SUM and ending with the statement labeled SUM. If Qp (in location HTCAP after the THROUGH loop is completed) is greater than QR (344032), the table is finished and the program transfers back to the statement labeled START to read new values of PCXS and DELT. If HTCAP is not greater than 344032, a one line table entry is made consisting of TF and the efficiency. This is accomplished by the PRINT statement labeled LAST. In this case the program will not return to START but will automatically return to the THROUGH LAST statement where TF will be incremented by an amount DELT. The statements which follow the statement labeled LAST are merely the VECTOR VALUES statements which preset the formats for input and output statements, i.e. READ and PRINT statements, in the executable part of the program. The data cards for the program follow the $DATA card. In this case, two tables will be prepared, one with PCXS = 20% and DELT=2000F, and the other with PCXS = 10 and DELT=1000F. Flow Diagram CALCULATE MOLES OF EACH GAS IN STACK GAS PER r \ BREAD IMOLE OF METHANE START ) PCXS DELTCD 0: MOLOUT(l) = 2.0*FXCS/100. N2: MOLOUT(2) = 2.0*(100+PCXS)*(79. /21.)/100. H20: MOLOUT(3) = O.017*2.0*(100+PCXS)*29./(18.+100. )+2. C02: MOLOUT(4) = 1.0 E12

Example Problem No. 88 Flow Diagram (continued) PRINT PCXRINT Mo CONVERT INLET /HHIGH IAST\ CQNVERT TF TO ZERO PCXS & MOLES / \ O ACCUMULATING ( IOF EACH L TEMP TO ~K _TF=7-0 DELT, K VARIABLE ICOMPONENT IN TZ=(70+460YL.8 TF5000 / = (TF+460.) HTCAP=O THROUGH SIM\ EVALUATE QF Table Entry =< Ii= 1, 1 HTCAP = HTCAP+1.8*MOLOUT(I), ST >4 +f Cp(I)dt Efficiency TZ START MAD Program and Data BRICE CARNAHAN S225C 001 010 000 FURNACE BRICE CARNAHAN S225C 001 010 000 FURNACE $COMPILE MAD, PRINT OBJECT, EXECUTE, DUMP FURN 001 R R FURNACE EFFICIENCY AS A FUNCTION OF STACK GAS TEMPERATURE R DIMENSION MOLOUT(4) INTEGER PCXS, I VECTOR VALUES A(1) = 8.27, 6.50, 8.22, 10.34 VECTOR VALUES B(1) =.000258,.001, *J0015,.000274 VECTOR VALUES C(1) = 0.0, 0U0,.0000 134, 0. VECTOR VALUES D(1) = -18770U., U.0,.0, -195500. R R THE EXECUTABLE PART OF THE PROGRAM BEGINS HERE R START READ FORMAT INPUT, PCXS, DELT MOLOUT(1) = 2.*(PCXS/100.) MOLOUT(2) = 2.*(100 + PCXS)*(79./21.)/100. MOLOUT(3) = 0.017*2.0*(100 + PCXS)*29./(18.*100.) + 2. MOLOUT(4) = 1.0 PRINT FORMAT HEDINGPCXS, 0.0, MOLOUT(1)*..MOLOUT(4) TZ = 530./1.8 THROUGH LAST, FOR TF= 70.,DELT,TF.G.5000. T = (TF+460.)/1.8 HTCAP = 0. THROUGH SUM, FOR I = 1, 1, I.G.4 SUM HTCAP = HTCAP + ( A(I)*(T-TZ ) + B(I)*(T.P.2-TZ.P.2)/2. + 1C(I)*(T.P.3- TZ.P.3)/3. +D(I)*(1./TZ - 1./T))*1.8*MOLOUT(I) WHENEVER HTCAP.G. 344032., TRANSFER TO START LAST PRINT FORMAT RESULT, TF, (344032. - HTCAP)/344032.*100. R R FORMATS FOR THE READ AND PRINT STATEMENTS R VECTOR VALUES INPUT = $I5,F10.4*$ VECTOR VALUES HEDING = $22H1PERCENT EXCESS AIR = I4/45HOSTAC 1K GAS COMP. FOR ONE MOLE METHANE IN FEED /30HOMOLES OF METHA 2NE = F6.4/3UH MOLES OF OXYGEN = F6.4/ 330H MOLES OF NITROGEN = F6.4/ 30H MOLES OF WATER VA 4POR = F6.4/30H MOLES OF CARBON DIOXIDE = F6.4/ 5 9H4TEMP.(F) S5, 15HPCT. EFFICIENCY //*$ VECTOR VALUES RESULT = $1H FU1, F16.1*$ END END OF PROGRAM $DATA 20 200. 10 100. E13

Furnace Efficiency as a Function of Stack Gas Temperature Computer Output Output for Second Data Card The computer output has been modiPERCENT EXCESS AIR = 10 fled slightly to conserve space. STACK GAS COMP. FOR ONE MOLE METHANE IN FEED MOLES OF METHANE = 0,0000 MOLES OF OXYGEN = 0.2000 MOLES OF NITROGEN = 8.2762 MOLES OF WATER VAPOR = 2.0603 MOLES OF CARBON DIOXIDE = 1.Uo00 Output for First Data Card TEMP.(F) PCT. EFFICIENCY 70,0 100.0 170,0 97,6 PERCENT EXCESS AIR = 20 270,0 95.1 370,0 92.6 STACK GAS COMP. FOR ONE MOLE METHANE IN FEED 470,0 90.1 470,0 9011 570,0 87.5 MOLES OF METHANE = 0.0000 670.0 85.0 MOLES OF OXYGEN = 0.4000 70,0 82.4 MOLES OF NITROGEN = 9.0286 870 0 79.8 MOLES OF WATER VAPOR = 2.0657 970,0 7 1170,0 71.8 1270.0 69.1 1370.0 66.3 1470,0 63.6 TEMP.(F) PCT. EFFICIENCY 1570,0 60.8 1670,0 58.0 70,0 100.0 1770,0 55.2 270,0 94 7 1870,0 52.3 470.0 89.3 1970.0 49.4 670.0 83 8 2070,0 46.5 870,0 78.2 2170,0 43.6 1070,0 72.5 2270,0 40.6 1270,0 66 6 237000 37.6 1470,0 60.7 2470.0 34 6 1670.0 54.7 2570.0 31.5 1870,0 48.6 2670.0 28.5 2070,0 42.3 2770,0 25.3 2270,0 36,0 2870.0 22.2 2470,0 29 5 2970.0 19.0 2670,0 22.9 3070,0 15.8 2870,0 16.2 3170.0 12o6 3070,0 9o3 3270.0 9.3 3270,0 2.4 3370.0 6.0 3470,0 2.7 Discussion of Results The printed output from the program is shown for two different sets of input parameters. The first solution is for the case of 20% excess air and a 200F~ temperature interval in the table, i.e. the solution for the originally stated problem. The second set is for 10% excess air and a 100F~ table increment, and illustrates that with an almost trivial increase incomplexity (reading PCXS and DELT from data cards rather than incorporating them as constants into the body of the program) a computer program can be generalized to solve a more general class of problems of which the one at hand is simply a specific case. Note also that since the adiabatic flame temperature of the combustion gases is simply the temperature at which the furnace efficiency equals zero, this simple program with three or four additional statements could calculate the flame temperature of any mixture of air and methane as long as the original built-in assumptions (complete combustion, 50% relative humidity etc.) were satisfied. These restrictions could of course be removed also if additional generalization were introduced. E14

Example Problem No. 89 MASS OR HEAT TRANSFER BY DIFFUSION by M. J. Sinnott Department of Chemical and Metallurgical Engineering The University of Michigan Course: Physical Metallurgy Credit hours: 4 Level: Junior or Senior Statement of Problem Iron at elevated temperatures and in a CO atmosphere will dissolve carbon interstitially to produce a carbon gradient in the iron. The gradient produced is a function of time, and the diffusion equation relating concentrations, time and the diffusion constant is as follows: (C - CO)/(C1-CO) = 1 - erf(X/2(Dt)0 5) where C = Concentration of carbon at depth X CO = Initial carbon concentration in the iron C1 = Carbon concentration at the surface of the iron D = Diffusion constant t = Diffusion time erf = error function A library subroutine is available in The University of Michigan Executive System to calculate erf(x) for a given value of x. For problems where the error function is not used directly as such, but which requires the "inverse error function," (i.e., the inverse error function, x = erf- (y), is defined herein to mean the value, x, for which erf(x) = y) a more involved system of programming is required in order to obtain solutions. The problem given here is of this type. The following are experimental data of carbon concentration versus depth below the surface of an iron sample that was carburized for 10 hours at a temperature of 17000F. The initial carbon concentration in the iron was 0.0%, while the surface carbon potential is 1.50%. Distance Below Carbon Conc. Surface, Cm Per cent 0.025 1.00 0.050 0.80 0.075 0.59 0.100 0.42 0.125 0.30 0.150 0.20 0.175 0.14 0.200 0.10 0.225 0.08 0.250 0.03 0.275 0.01 Evaluate the diffusion coefficient, D, from the above data. El5

Mass or Heat Transfer by Diffusion Solution The method of solution is self explanatory from the flow diagram, MAD program, and computer output which follow. Since there is no subroutine to generate the inverse error function the program is written in such a fashion that a value is assumed for the error function, the argument calculated, then compared with the value being sought, until the assumed value yields an argument which agrees with the argument being sought to as close an approximation as is desired. The flow diagram, program and computer output is as follows: Flow Diagram NT PROUGH \INT ST.qR -- READ N, T N C END, FOR \ READ CO, C1 J=,1, /CONCDIST ETA AER3./DIST (C-CO)/(C1EPS-CO)-\ ERR. OR..ABS.(ERF. (3K)| E. \RFC=LE. 0.0001 DT.ABS. ERF.(3K) T — /ERR.LE. 0.0001 |. F=3./K ---- K=l^ (E 1CONC, --- U=3./K-1 2.0 TEST=ERF. ( -.-U) L-ERTEST.Z.ERRE) F _ ABS. (TEST- ___ TA ERR).L.O0.0001 F=~ U= AMMA ( BETA )^ AERF=3./K --- EPS ( DEL TA)^ AERF= - - EPS /-~^\ I --- —------- PRINT ^EPS D=- (4.0*T*AERF*AERF) -F.( ),) El6

Example Problem No. 89 MAD Program STAR.T t AD FLt,'RMAT I NPUT, N, T I CC( C1 VECTCR VALUES INPUT: I5, El.3, 2F10.5-$ _,,,____ _._ __'.. LC^ R_- k IJNPT _- A_ I5. ClC.._3__ _2F1lQ_._5 _$ _ ___________________________. __ F'RINT FUlRMAT CUT, ~, T, CC, C! '__CT(-R VALUES JOUTJ = tHiO, 13tlO)ATA PCOINtS = I5, S5, 6,TIMh =- E 11C.3, S5, 2C0-hjASL CCJ\;CENTRATIC4 -= F1u.b, S5, 23HSURFACE CUNCL 2FN [tAT I[!'i = F1C.5:$ ThRCOUGH ENDj, FCR J =1,1, J.G.N REA.) F RNIAT UATA, C, X -.J0-J ^.-. --- —-- -..-.-..-... --- —-- --- -------------------------- VLCTCR VALUES DAT A = $ 2F C.5$ PRINI FRKMAT CRIl., Ci X VECTOR VALUES CRIG = HO, S1, 1 HCONCENT RATIN = F1. 5, S10, 31CLIS1ACNCE = F1,t.5*4. --- —— ~_ --- —..-__ — ------- 3 1 C F t) I S I A N C I C, 5 I TE G E E J,N,K ERFC = (C-CC)/(C1-CC) E, R = I.C - E FC LCCP Tl-RGUJ- L LOCP, FOR K= 1,, ERF.(3.0/K).L. ERRR.CR..ABS. (ERF. (3. 1C/K)- LRR).LE. (,.CCG YHE'EVF.AbS.(ERF. [3.C/K) -ORR).LE./.COC1, TRANSFER TO ReTA L 3.C/K i,-' ENE V- R K.E....... --- —-—.-......... --- —. —. —.-. --- —---;=3. 5 ('t T - F K a I SE l = 3.('/(K-1) FND GF CONLC II NAL ~ ~ -,..... -._-.... ---- - - - - - - - - - - - - - - - - - - - - - - ------- - - - GAMMA IFSTr = E F. ( ( L +U)/2.(i ) _ iE'EVIR TESI.L. ER _ _ L = (L+U)/2.C i, = (L+U)/2.0 TRANSFIR' TC uAMMA LNL OF C}NUDITIGNAL _ BEIa A ERF = 3.0/K TRA,'SF.R FC EVS CELFA iE AREF = (L+U)/2.0......-_. ---_ J-I- J ----'-~-~ —.-. ^ =~ + --- —-----------------------— _ --- —-—. --- —- --- —..-. -. EPS 1. = X*(/(4.C*T*AER'F*AEQFF) EN.) PR INi[ RR MAT CUTPUT, O, C, X, AERF, ERF.(AERF), E RR --— _-_ --- — -— I - -------- p_ OK-MA-T —GUJ-pjJJJ —c-' -- R --- — VECTOR VALUES CUIPUT =$1H, 3hU = F10.3, S3, 6HCt:NC = FlC.5., 4S3, (6t)IST = F10.5, S3, 6HAERF = F1C.5, S3, 12HEiF.(AERF) = F 51C.', S3, 5HiTRR = Fl1:.5*t _ TRANSF-R TC STARl 'NDU OF POoR0, Ai' E17l

Mass or Heat Transfer by Diffusion Computer Output DATA _POINTS = 11 TIME j.360E 05 BASE CONCENRAl ION =.OOoc SSRFACE CONCtNTRAfION = 1.300 --- ----- CONCENTRArIGN = 1.0000 O I STANCE =__.02500_ D =.101E-06 CONC = 1.00CO0C DIST =.02500 AERF =.20748 LRF.(AEIF) =.23080 ERR =1.23077 CONCENTRATION =.8000 CISTANCE =.OCOO O =.138E-06 CONC =.8UCOO D1ST =__500__ 0 AERF 3=5531 ERF. (AERF) =.38467 ERR =.38462 COjNCENTRATICN =.590C0 _ISTANCE=.0 7500.. - =.139E-06 CONC =.59300 T =.075'0( AERF.=529t9 ERF.(AERF) =.54620 ERR 46 CONCENIRATION =.420CO DISTANCE =.1000 D =.142E-06 CONC.425000 DIST ='10000 AERF =.29 8 ERF.(AEfF) =.7697 ERK =.7692 CONCENTRATICN -=.30uC DISISANCE. 12_500 D =.151E-06 CONC =.30CCO 0 IS rT =.1250G AERF =.:84741 ERF.(AERF) =.76925 ERR =.16923 CONCENTRATICN =.20C00 DISTANCE =.15000 _ _.154E-06 CONC =.C42C000 IST =.15000 AERF =.30804 _ERF.(AEAF) =.'4622 ERR ='.84615 O:. 51'- — CONCENTRATION =.IOCOO) DIS0ANCE 2=:.20000 'E - 0 —.178E-06 CON =1.C10000 DIST =.20000 AEF = 1.25098 EFRR =:.92308 CONCENTRATICN =.000C DISTANCE =.22500 D =.214E-06 CONC =.C08000 OIST =.12500 AERF = 1.1322270 FEF.(AERF) =9,d6232 ERR =.923846 CONCENTRATICN =.0300 DISlANCL =.25000 - 0 =.168E-06 CONGC. 3CC000 01S T.215OU AERk = 1.61093 -F".(/AEPF).97695 ER =.9792 CO INCENTRATI OC =.10000 DISIANCFC =.27500 D =.148E-06 CONC =.CCOO0 DIST =.205000 A7ER2 F = 1.E237 E F.(AEF =.933 ERq =.99231 CtNCENrRATIrN = ~0300C DISIA:4CL =.2500b.168E- 06 COC _. 0_0 _S _=_ __.25C0U E:RF:1.b063 EkF. _...... 7E:: =..._.. 7 2 rCONCENTRATIO~ =_.OlO CI *Cv _ 1S!1 Arx:E -.27l)0 6 D =.140E-06 (:;U4C:.C1CO Fl)IST =1. 27500- AERF = 1.S 72 f~F. —(-~E-i<-tF) =:.. L9"7 ERtP =,9923i Discussion of Results The results obtained above required 0.70 minutes of IBM 704 computer time for compilation of the MAD program and 0.40 minutes for execution. By contrast, the identical program required 0.32 minutes for compilation on the IBM 709 and 0.58 minutes for execution. This same program can be used to solve unsteady-state heat transfer problems. In place of C, CO, and Cl where these are concentrations in the diffusion problem, one can use temperatures. For the diffusion constant, one can use the thermal diffusivity. E18

Example Problem No. 90 PREDICTING THE SCRAP REQUIREMENT FOR THE OXYGEN STEEL CONVERTING PROCESS by Robert D. Pehlke Department of Chemical and Metallurgical Engineering The University of Michigan Course: Metallurgical Process Design Credit hours: 4 Level: Senior Statement of Problem The computer analysis of oxygen steel converting has recently become operational at a Pennsylvania plant of one of the larger steel companies. The program considers 13 variables which include hot metal analysis, temperature, oxygen flow rate, process geometries, and several others which influence the control of the converting process. The commercial computer program is highly empirical, and although of limited accuracy, it has produced marked improvement in control of the operation. In an effort to illustrate the magnitude of the process, and some of the thermochemical interactions which play a role in design and operation of oxygen converters, a problem was selected for computer analysis which requires a mass and thermal balance for the oxygen steelmaking process. The relationship selected for investigation is the influence of silicon content of the hot metal on the required scrap addition. Computation of the mass balance involves a trial-and-error calculation, and the thermal balance requires lengthy calculations of the heat effect of each reaction, a task suitable for the digital computer, Write and test a MAD Program which will permit a calculation of the scrap requirement for an oxygen-steelmaking converter, and use the program to determine the influence of silicon content of the hot metal on scrap consumption. The following assumptions may be made. Heat losses are 2 x 106 BTU/hr Production rate is 80 tons/hr Flue dust losses are 3% of hot metal charged and exit at hot metal charging temperature. Use the following data to test your program. %Si %Mn %C %P Temperature Hot metal 0.5-2.0 1.30 4.40 0.13 2500~F Blown metal 0.02 0.20 0.05 0.02 3000~F Scrap Same as Blown metal 77~F Oxygen 99.5% oxygen 77~F Burnt lime Available CaO = 92% 77~F Mill Scale (0.15 lb/lb-burnt-lime) 50% FeO 77~F Slag Basicity ratio (%CaO/%Si02) = 3.0 3000~F Combustion gas 100% CO 5000~F El9

Scrap Requirement for the Oxygen Steel Converting Process The following thermodynamic data may be used. 1. Heats of Reaction at 77~F AH,BTU/lb-mole P2(g) + 3 CaO + 5/2 02(g)- >Ca3P208 -1,016,000 2 Fe + 3/2 02(g). — >Fe203 -354,000 Fe + 1/2 02 - FeO -114,800 Si + 02(g) -- Si2 -369,700 Mn + 1/2 02(g) -- MnO -165,600 2CaO + SiO --- 2CaOSi02 -53,640 FeO + Si02 > FeSiO3 -8,100 MnO + Si02 -- > MnSi0 -10,600 C(gr) + 1/2 02(g) - CO(g) -26,416 C(gr) + 02(g) -- C 02(g) -94,052 2. Heats of Solution, BTU/lb-mole T in ~F. Si77oF - SiT =H = 6.7 T - 34,150 Mn77oF _ >MnT AH = 11.0 T + 2,900 1/2 P2(g)o77F -- PT AH = 4.45 T - 92,500 C(gr)77oF > _ -T AH = 5.23 T + 10,020 3. Sensible Heats - Reference Temperature 77~F, T in ~F, Iron-Base Metal (liquid) AH = 0.184 T + 63 BTU/lb. Slag (liquid) AH = 875 + (T - 2820) x 0.3 BTU/lb. CO LH = 23,200 + (T - 3000) x 8.5 BTU/lb-mole CO2 AH = 37,400 + (T - 3000) x 14.4 BTU/lb-mole N2 AH = 23,000 + (T - 3000) x 8.6 BTU/lb-mole Solution The diagramming of a sequence to solve this problem requires a relatively complete understanding of the steelmaking process. The oxygen converting process for producing steel had been discussed in detail in class prior to presenting this problem, which is a detailed examination of the process in terms of its mass and energy requirements, and involves a complete examination of one of the principal process variables. The solution requires that a mass balance be completed which satisfies the production rate to which the heat losses are tied. Upon completion of the mass balance, the excess heat energy is calculated to determine the amount of scrap which is required to maintain the process temperature below its thermal limit. This scrap plus the refined steel produced by oxidizing the impurities from the blast furnace hot metal charged to the converter must satisfy ihe production specified. E20

Example Problem No. 90 The mass balance involves a trial and error computation since the burnt lime added depends upon the silicon which is oxidized from the hot metal, and the basicity ratio of the slag. Mill scale, added to supply additional oxygen to the process, is defined in terms of the lime added. The thread which runs through the several input and output streams of the converting process is the requirement that the charged iron, a constituent of all of the condensed phase streams (i.e., all except the oxygen lanced into the converter and the combustion gases exiting from the process ) must be equal to the iron removed. The mass balance is then defined in terms of an iron balance. The method of calculation is as follows. 1. Calculate the moles of silicon, manganese, carbon, and phosphorous in the hot metal which is assumed to be charged, correcting its value for any flue-dust losses. 2. Initialize the weight of scale added, and the moles of iron in the slag stream at zero. 5. Correct the value for the weight of iron in the hot metal by adding the weight of iron in the scale and subtracting the weight of iron in the slag. 4. Calculate the weight of refined steel produced from the hot metal charged on the basis of an iron balance. 5. Calculate the moles of silicon, manganese, carbon, and phosphorous which remain in the steel, subtracting them from the resultscf step 1 to determine the moles which enter the slag stream. 6. Calculate the mole fraction of iron oxide in the slag using a relationship involving the carbon content of the finished steel. 7. Calculate the weight of burnt lime added to provide the correct basicity ratio in the slag, assuming silica to be the only acidic slag component. 8. Calculate the weight of scale added on the basis of the result of step 7, and calculate the weight of slag produced. 9. Compare the weight of scale calculated with that used in step 3. If they agree within a limit of error continue on; otherwise, return to step 3, substituting the result of step 8. 10. Calculate the composition of the slag in terms of the compounds which form, assuming a preferred order to the reactions. 11. Compute the enthalpy of each of the streams, the hot metal charged, the slag, combustion gases, and steel produced. The inlet enthalpy of the hot metal (including the oxidizing reactions which occur) minus the enthalpy of the slag, combustion gases, steel, and heat lost from the vessel will yield the excess energy which is available to melt scrap. E21

Scrap Requirement for the Oxygen Steel Converting Process 12. Calculate the scrap which can be melted, and the metal produced which is the scrap plus the refined steel. 13. Compare the metal produced with that associated with the reported heat loss. If they do not agree, correct the assumed hot metal charge proportionately and return to step 1. Flow Diagram Calc. Print Convert NM' SCALE = 0 Calc. conc. in steel Calc. lime weight _ NSFEO = 0 Calc. moles in slag > Calc. scale weight F ^ ^<~~Calc. moles gas Calc. slag weight SCALE-WSCALEI >DELT WTH R ead Data HCAeG a- dTTP AML PHA /2 IOTA ^J i~e., cump- "~ NSCAO>(2*NSSIO + 3*NSP205) Calc. NMSI NMNSI <NSMNO 3 Calc. RNFEadSI. 0 H —) - p WTHM NFE I 3 N AE2SI = Calc. enthalpy Calc. WTSCP, WTH1VT Print CIOTA molSCALE = WSCALE - ALPHA BETA Ot -pDELTt PSC E22

Example Problem No. 90 List of Variables HM(1)...HM(3) Linear array for storage of the words, "hot metal" HM(4) Weight percent silicon in the hot metal HM(5) Weight percent manganese in the hot metal HM(6) Weight percent carbon in the hot metal HM(7) Weight percent phosphorus in the hot metal HM(8) Temperature of the hot metal in OF SCRAP(1)...SCRAP(3) Linear array for storage of the word, "scrap" SCRAP(4) Weight percent silicon in the scrap SCRAP(5) Weight percent manganese in the scrap SCRAP(6) Weight percent carbon in the scrap SCRAP(7) Weight percent phosphorus in the scrap SCRAP(8) Temperature of the scrap in OF STEEL(1)...STEEL(3) Linear array for storage of the word, "steel" STEEL(4) Weight percent silicon in the steel STEEL(5) Weight percent manganese in the steel STEEL(6) Weight percent carbon in the steel STEEL(7) Weight percent phosphorus in the steel STEEL(8) Finishing temperature of the steel in ~F POXY02 Percent oxygen in the oxygen stream A Available base in the lime flux PSCFEO Weight percent FeO in the scale RC02 C02/CO ratio in the stack gas T02 Temperature of the oxygen in OF TSC Temperature of the scrap in OF TS Temperature of the slag at tap in OF TCG Average temperature of the combustion gas in OF HL Heat loss in BTU's per hour PROD Production in lbs. per heat FDLOSS Flue dust lost as percent of hot metal weight charged WTHM Estimated weight of hot metal to be charged in tons DELT Allowable error in calculated scale weight EPS Allowable error in calculated production BR Basicity ratio, moles CaO/moles Si02 WTHM Weight of hot metal, converted to lbs. in computer program NHMSI Moles of silicon in hot metal NHMMN Moles of manganese in hot metal NHMC Moles of carbon in hot metal NHMP Moles of phosphorus in hot metal WHMFE Lbs. of iron in hot metal charged SCALE Weight of scale NSFEO Moles of FeO in slag WST Weight of steel in lbs. NSTSI Moles of silicon in steel NSTMN Moles of manganese in steel NSTC Moles of carbon in steel NSTP Moles of phosphorus in steel E23

Scrap Requirement for the Oxygen Steel Converting Process List of Variables, Continued NSSIO Moles of silica in slag NSMNO Moles of MnO in slag NSP205 Moles of P205 in slag NCCG Moles of carbon in combustion gas XSFEO Mole fraction of FeO in slag NSCAO Moles of lime in slag WTLIME Weight of lime charged in lbs. WSCALE Weight of scale charged in lbs. WTSLAG Weight of slag in lbs. NOXYO Moles of 0 supplied by oxygen stream SCFOX Standard cubic feet of oxygen required for heat N2CS Moles of 2CaO.Si02 in slag NMNSI Moles of MnO.Si02 in slag NFE2SI Moles of 2FeO.Si02 in slag HHM Enthalpy of charged hot metal HS Enthalpy of finishing slag HCG Enthalpy of combustion gases HST Enthalpy of finishing steel HSCRAP Heat available for melting scrap WTSCP Weight of scrap to be charged WTOUT Lbs. of steel produced PSC Percent scrap to be charged TSCP Tons of scrap to be charged E24

Example Problem No. 90 MAD Program and Data R. Do PEHLKE C07-N 0 3 030 000 CHE 13 R. D. PEHLKE C07-N 0 3 030 000 CHE 13 $COMPILE MAD, PRINT OBJECT, EXECUTE, DUMP DIMENSION HM(10)9 STEEL(10)* SCRAP(10) READ FORMAT INPUT1 HM(1)**.HM(8), STEEL(1) **STEEL(8 ) 2SCRAP( 1 )..SCRAP(8) VECTOR VALUES INPUT1 s $(3C6,4F6.3,F6.0)*$ READ FORMAT INPUT2, POXY02, A, PSCFEO.BRRCO2T02*,TSCTSTCG VECTOR VALUES INPUT2 = $5F8.4o4F50*$ READ FORMAT INPUT3# HLPRODFDLOSStWTHM VECTOR VALUES INPUT3 = $4F15*4*$ READ FORMAT INPUT4,DELT EPS VECTOR VALUES INPUT4 = $2F10*4*$ PRINT FORMAT TITLE VECTOR VALUES TITLE = $1Ht165H INFLUENCE OF SILICON ON SCRAP 2ADDITION TO BASIC OXYGEN CONVERTER ////*$ PRINT FORMAT DATA1,HM( 1,*HM(8) STEEL( 1)*..STEEL(8), 2SCRAP(1) 1*SCRAP(8) VECTOR VALUES DATA1= $5H DATA//S14,31H COMPOSITIONS IN WEIGHT 2 PERCENT/S19,3H SIS4.SH MNS4,2H CeS4*2H P*S3S5H TEt4P//(3C6# 34F6.3,F6,0)//*$ PRINT FORMAT DATA2, POXY02, A, PSCFEO, RC02, T02, TSC, TSt 2TCG, HL, PROD, FDLOSS. WTHMs DELT. EPS, BR VECTOR VALUES DATA2 = $26H PERCENT OXYGEN IN BLAST =F6.2/33H 2PERCENT AVAILABLE BASE IN LIME =F6~2/23H PERCENT FEO IN SCALE 3 =F6.2/33H C02/CO RATIO IN COMBUSTION GAS =F6.2/24H TEMPERATU 4RE OF OXYGEN =F6.0/23H TEMPERATURE OF SCRAP =F6.0/22H TEMPERA 5TURE OF SLAG =F6.0/32H TEMPERATURE OF COMBUSTION GAS =F6.0/ 620H HEAT LOSS. BTU/HR =F8.0/21H PRODUCTION, LBS/AR =F7*0 /35H 7 FLUE DUST LOSS, PCT OF HOT METAL =F6.2/43H WEIGHT OF HOT MET 8AL TO START TRIAL, TONS =F6.0/44H ALLOWABLE ERROR IN SCALE WE 9IGHT, DELT, LB =F6.0/50H ALLOWABLE ERROR IN WEIGHT OF METAL O 1UT. EPS, LB =F8.0/25H BASICITY RATIO OF SLAG = F6.4////*$ PRINT FORMAT HEAD VECTOR VALUES HEAD = $S5,20H PCT SI IN HOT METAMS3,20H PCT S 2CRAP IN CHARGES3.17H TONS SCRAP ADDED//*$ R MASS BALANCE WTHM = WTHM*2000, DELTA READ FORMAT INPUTS, HM(4) VECTOR VALUES INPUT5 = $F6.3*S GAMMA NHMSI =(100 -FDLOSS)/10 *WTHM*HM(4)/(100.*28 06) NHMMN (100*-FDLOSS)/10 **WTHM*HM(5)/(10 **54*94) NHMC = (100~-FDLOSS)/10O*WTHM*HM(6)/(100O*12.01) NHMP = ( 100-FDLOSS )/100*WTHM*HM(7)/(100O*30.98) WHMFE a (100~-FDLOSS)/10 **WTHM*(100,-HM(4)-HM(5 -HM(6)-HM(7) 2)/10, SCALE = 0.0 NSFEO = 0.0 ALPHA WHMFE = WHMFE+SCALE*PSCFEO*55.85/( 100.*71 85)+SCALE*( 100.2PSCFEO)/100*-NSFEO*55#85 WST = 100.*WHMFE/(100*-STEEL(4)-STEEL(5)-STEEL(6)-STEEL(7)) NSTSI = WST*STEEL(4)/(10.*28.06) NSTMM = WST*STEEL(5)/(10 **54.94) NSTC = WST*STEEL(6)/(100o*12*01) NSTP = WST*STEEL(7)/(100s*30.98) NSSIO = NHMSI-NSTSI NSMNO = NHMMN-NSTMN NSP205 = (NHMP-NSTP)/2i NCCG = NHMC-NSTC E25

Scrap Requirement for the Oxygen Steel Converting Process MAD Program and Data (continued) XSFEO = O.01/STEEL(6) 81fJ8 s?lii?5O+SMNO+NSP205+NSCAO)*XSFEO/( 1-XSFEO) WTLIME ' NSCAO*56.08*100./A WSCALE 0.15*WTLIME WTSLAG = NSSIO*60.6+N:SMN:O*70.94+NSP205*141*96+NSFEO*71.85 2+WTLIME WHENEVER*ABS. (SCALE-WSCALE).G*DELT# TRANSFER TO IOTA NOXYO = NSSIO*2 +NSMNO+NSP205*5.+NSFEO+NCCG* ( 1 +2. RC02)/( 1.+ 2RC02)-WSCALE*PSCFEO 16* /471*85*100) SCFOX = (IOXYO/2. ) 379*1 00 /POXY02 N2CS a NSSIO WHENEVER NSCAO.* ( 2.^*NSSIO+3.*NSP205) NMNSI = 0.0 NFE2SI = 0.0 OTHERWISE N2CS = (NSCAO-3.SNSP205)/2. NNSI ' NSSIO-N2CS WHENEVER NMNSI.*4NSMNO NFE2SI a NSSIO-N2CS-NSMN0 NMNSI x NSMNO OTHERWISE NFE2SI a 0.0 END OF CONDITIONAL END OF CONDITIONAL R HEAT BALANCE HHM u0184*HM(8)+.63* )*WTHM*(1 00.-FDLOSS ) /100.+(67*HM( 8 ) 2-3l4150 ) *NSSIO+( 11 HM( 8 )+2900. ) *NSMO+( 4.45*HM( 8 )-92500. ) * 3NSP205*2.+( 523*HM(8)+10 20.)*NCCG HS = (875*+(TS-2820e)*03)*WTSLAG-1016000**NSP205-53640.* 2NSSIO- 114800**(NSFEO-WSCALE*PSCFEO/ 100**71*85 ) -369700.* 3NSSIO-165l600 **NSMfO-5364 **N2CS-10600*NMNSI-8100.*NFE2SI HCG = (23200*+(TCG-3000)*8,5-26416,)*NCCG*(1l/(1.+RC02)+ ( 23740 *:+TCG-3000,)*14&4.-94052 )*NCCG*(RC02/1+RC2 RC2) HST a (0,184*TS+63*)*WST HSCRAP HHM-HST-HS-HCG-HL WTSCP = HSCRAP/{(0184*(STEEL(8)-TSC)+63.+14.168) WTOUT WST+WTSCP PSC * 100**(WTSCP/(WTSCP WTHM)) TSCP = WTSCP/2000, WHENEVERABS. iWTOUT-PROD).LEPSITRANSFER TO BETA WTHM z WTHM*PROD/WTOUT TRANSFER TO GAMMA BETA PRINT FORMAT OUTPUT. HM(4)) PSC* TSCP VECTOR VALUES OUTPUT a $F162tF24*.2F22*2*$ TRANSFER TO DELTA IOTA SCALE = WSCALE TRANSFER TO ALPHA END OF PROGRAM $OATA HOT METAL 500 1300 440 130 2500 STEEL 20 200 05 20 3000 SCRAP 20 200 05 20 77 995 92000 500000 3 77 77 3 3000 2000000 160000 3 60 100 1000 000 200 400 600 800 1000 1200 1400 1'00 1800 2000 E26

Example Problem No. 90 Computer Output INFLUENCE OF SILICON ON SCRAP ADDITION TO BASIC OXYGEN CONVERTER DATA COMPOSITIONS IN WEIGHT PERCENT SI MN C P TEMP HOT METAL *500 1.300 4.40 *130 2500. STEEL,020 0.200 0*05 *020 3000. SCRAP *020 020.0 0.05.020 77* PERCENT OXYGEN IN BLAST a 99*50 PERCENT AVAILABLE BASE IN LIME = 92 00 PERCENT FEO IN SCALE a 50.00 C02/CO RATIO IN COMBUSTION GAS a 0.00 TEMPERATURE OF OXYGEN = 77. TEMPERATURE OF SCRAP = 77, TEMPERATURE OF SLAG = 3000. TEMPERATURE OF COMBUSTION GAS = 3000. HEAT LOSS, BTU/HR =2000000, PRODUCTION, LBS/HR =160000* FLUE DUST LOSS, PCT OF HOT METAL - 3.00 WEIGHT OF HOT METAL TO START TRIAL, TONS - 60. ALLOWABLE ERROR IN SCALE WEIGHTS DELT# LB = 10. ALLOWABLE ERROR IN WEIGHT OF METAL OUT, EPSI LB = 1000 BASICITY RATIO OF SLAG 3,.0000 PCT SI IN HOT METAL PCT SCRAP IN CHARGE TONS SCRAP ADDED 0.00 9,31 8.08 0.20 12*09 10050 0 40 14*55 12.66 0,60 1:6.89 14*71 0,80 19,10 16466 1.00 21.19 18.51 1,20 23.18 20.27 1.40 25*07 21*95 1.60 26687 23*55 1,80 28.59 25.08 2000 30.22 26.54 Discussion of Results The results printed out by the computer are the percent silicon in the charged hot metal, the percent scrap in the charge which is the calculated scrap divided by the sum of the scrap and the hot metal charged, and the actual number of tons of scrap added to the operation which produces 80 tons per hour in the case under consideration. As the silicon content of the hot metal increases, the excess heat available for melting scrap increases resulting in a higher scrap percentage in the charge as shown in the figure below. The structure of the computer program is such that other variables may be investigated, including the influences of changing the temperatures of the process. One commercial proposal for increasing the ability of the process to handle scrap involves preheating the charged scrap, E27

Scrap Requirement for the Oxygen Steel Converting Process thus decreasing the heat units required to accomplish melting. This program can be used to study the effect of preheating to various temperature levels. Other compositional variables can also be studied. Effect of Silicon Content of Hot Metal on Percent Scrap in Charge of 80 ton/hr Oxygen Converter 35. g4 I-i Percent Scrap slc H Mt in Charge::::::::: j:::::: 25. 20. i- ^:^^ ^:^^ 0-0 0.4 0.8 1.2 1.6 Weight percent silicon in Hot Metal

Example Problem No. 90 Critique Of the 23 solutions attempted by the class to which this problem was assigned, only three were successful. These three solutions employed essentially the same approach outlined above, and yielded results in reasonable agreement. Nonetheless, this particular problem is an excellent one for use in the terminal course in extractive metallurgy. The relatively unsuccessful performance of the students may be an indication that the problem is too difficult as a class assignment, and might be better used in its finished form as a supplement to the lecture material. The problem is pertinent to the material covered in the course and was attacked with enthusiasm by the students. The low success rate may serve as an indication that in long and relatively complex problems, additional assistance in programming, and problem outline may be absolutely essential. The bulk of the class attempted this problem on the IBM 704 system. Two groups of two attempted the problem on the LGP-30, but were not successful in the alloted time of five weeks. Although the class worked individually on the problem, in the case of the small computer the four students were all taking their first look at a new language, and the grouping was suggested. The interpretive language, ACT III, was used but all of the students experienced considerable difficulty in mastering the programming for the LGP-30. This particular experience is a re-enforcement of the view that the large machine, if available, with its user-oriented language and easy accessibility is by far the most desirable for the user. E29

Example Problem No. 91 PRECISION LATTICE PARAMETER DETERMINATION FOR A CUBIC MATERIAL by Jeremy V. Gluck Department of Chemical and Metallurgical Engineering The University of Michigan Course: X-Ray Studies of Engineering Materials Credit Hours: 3 Level: Senior Statement of the Problem Write a computer program for the precision calculation of the lattice parameter of a known cubic material. The following conditions apply: A. X-ray diffraction exposures are made in a 114.6 mm. diameter cylindrical DebyeScherrer camera. B. Plane indices can be determined from rough measurements made with a plastic scale; however, the lattice parameter calculations should be based on line position measurements made on the optical comparator to a precision of + 0.0025 cm. C. The film measurements should be corrected for shrinkage. D. The Nelson-Riley extrapolation technique should be used to correct for sample absorption errors. A least-squares fit should be employed in the extrapolation. E. An estimate should be made of the standard deviation of the parameter. F. A plot should be generated showing the relationship between the parameter values and the Nelson-Riley factor. G. Provision should be made for discarding diffraction angles of less than 30~. H. In the back-reflection region, the a,' and a2 doublets should be used where resolved. Solution For a detailed discussion of the X-ray relationships and procedures, refer to Klug and Alexander (Reference 1). The computer program for the IBM 709 was written in the MAD language. Variable names used in the MAD program are defined in the body of the discussion which follows. They are capitalized and enclosed in parentheses. The solution employs the Bragg diffraction equation: a XA = 2dL sin GL, where dL = 0 h2 + k2 + 2 and the Nelson-Riley factor (NELRI), / ___2 2 1 - L + -— LFL 2 \ sin L + L / E30

where: X = wave length of the radiation used (LAM1 or LAM2). dL = interplanar spacing between successive atomic planes. L = angle between the atomic plane and the incident and reflected X-ray beam (THETA(L) in radians; ANGLE(L) in degrees). ao = the lattice parameter (side length) of the cubic crystal unit cell (LGOOD). h,k,j = Miller indices of the diffracting plane; ((h2+k2+92) is often referred to as N, (NN)). L = the diffraction line in question (L). All calculations are made in terms of Angstrom Units, A, (10-8cm.). The Bragg equation can be rearranged as follows: a /h2 k2+ 2 aL = 2 sin QL For a given radiation, a (LAM1 or LAM2), a lattice parameter, aL (LATPA(L)), can be calculated for each diffraction line, L, which is produced at an angle GL (THETA(L)) by a particular set of planes, h, k, A. In a cylindrical Debye-Scherrer camera, systematic errors in G values are possible due to film shrinkage during processing, displacement of the sample, and absorption of X-rays by the sample. In addition, there can be random errors occurring during subsequent film measurement. Various methods proposed to correct for such errors are discussed fully in Chapter 8 of Reference 1. The premise of these methods is that functions of G exist for which systematic errors in G vary linearly. The precision of sin Q values increases as G increases and systematic errors in determination of G approach zero as G approaches 90~. Thus, it is possible to extrapolate parameters calculated at intermediate angles to the value of some function at G = 900, where the error is presumably zero (Reference 5). The function proposed by Nelson and Riley (see Reference 1, pp. 464-467) has been found to be highly linear down to very low values of G and "for all practical purposes" all diffraction angles can be used in employing it. However, for optimum results Klug and Alexander (Reference 1) do recommend that its use be confined, where possible, to angles between 9 = 30~ and 900. The extrapolation is made by the method of least-squares. An equation of the slopeintercept form L = m FL + ao relates values of the lattice parameter, aL (LATPA), and the corresponding Nelson-Riley factors, FL (NELRI). The intercept of this line at FL = 0 (i.e., 0 = 90~) is a (LGOOD), the "true" lattice parameter. E31

Precision Lattice Parameter Determination for a Cubic Material In addition to determining the best straight line through the points, an estimate of random errors can be made by calculating the standard deviation of the points from the line. Useful formulae for these computations are given by Moroney (Reference 2) and are as follows: (ZFL)(ZFLaL) - (aZ)(2FL2) a 2 2 oa ( )2 -( QFL - ZFL m (FLa) - (2FL)(ZaL)/Q ZFL2 (FL)2/Q whereZ refers to the summation of the function for all diffraction lines and Q is the total number of lines. The "standard deviation", O, is caloulated from the equation: ( J =] (mFL+ao-aL)2 Q-1 where the expression (mFL+ao-aL) is the vertical deviation of an individual parameter, aL, from the least-squares line at FL. The film for the X-ray exposure is mounted in the camera as shown below: Film Sample Camera X-Ray Source Path of X-Ray Beam Paths of Diffracted Beams Following film development and numbering of the diffraction lines, measurements are made from the strip of film which has the following appearance: Back Reflection Region Front Reflection Region Line No. 6 78 8 76 5 4 3 2 1 1 2 (((i ) ) I i i1} RAD(1) < RAD(2)!RAD(3),I = ----— RAD (7) CENT2CEENT1 29 = 180~ = 0~ (inlet hole) (outlet hole) E32

Example Problem No. 91 Since the line pairs are not necessarily symmetrical about the inlet and outlet holes, the centers are determined by averaging the distances between corresponding pairs of lines. The geometry of the camera is such that the distances RAD(1), etc. (in millimeters) are approximately equal to twice the diffraction angle (ANGLE). In the process of converting the radius measurements, a shrinkage correction factor is introduced to adjust the results so that the distance CENT1-CENT2 is 90~. The plane indices are obtained from the ASTM Powder Data File Card for the substance in question. This identification is made from "d" values (see the Bragg equation) which can be determined from readily available charts or tables. A preliminary rough determination of the line positions Is made for this purpose. For the computer run, the wave length of the radiation is known, and the resolution of the al(J=l) and a2(J=2) doublets at high diffraction angles is determined by inspection of the film. The input to the computer consists of the wave length of the radiation (in the example, Cobalt) and an identification card describing the sample and film exposure conditions, and defining the total number of diffraction lines obtained, Q, and the number of lines at Q>30~, R. This is followed by a series of cards, one for each line (L), indicating the line number, the sum of the squares of plane indices (NN), the wave length to be used in the computation (J=l or 2), and the comparator measurement (RAD(L)). The program then does the following: 1. Computes the shrinkage correction factor (SHRINK). 2. Determines the angle (ANGLE) for each diffraction line, corrects it for shrinkage and converts to radians (THETA). 3. Computes the lattice parameter (LATPA) for each line using the appropriate wave length. 4. Computes the Nelson-Riley factor (NELRI), for each line. 5. Plots the parameter versus the corresponding Nelson-Riley factor. 6. Tabulates the results. 7. Performs a least-squares analysis to determine the extrapolated parameter (LGOOD) and slope (SLOPE) of the least-squares line. 8. Computes the standard deviation (VAR). 9. Repeats steps 7 and 8 for only those angles greater than 30~. (Q is now set equal to R, the number of lines at 0>30~.) 10. Prints a summary of the results from steps 7, 8, and 9. 11. Repeats for the next set of data cards for a new film. The plot can be used for the elimination of obviously incorrect points by inspection and the program can then be rerun with the corresponding data card or cards eliminated. E33

Precision Lattice Parameter Determination for a Cubic Material Flow Diagram Beginning\ LAMi, FILMNO,SPECNO,RADIAT, L \ NN(L),J(L) of LAM2 EXPTIM,CENT1,CENT2, RIOTA D ( Program FILTER, Q,R L>Q RAD(L) BETA ANGLE(L)=(CENT1-RAD(L))* z = 0 SHRINK = L =, 1 SHRINK10 y = 9/(CENT1-CENT2) L>Q THETA(L)=ANGLE(L)*0.01745329 T LATPA(L) = T LATPA(L) = 2 ( J(L)=M1 * LAM1 * [ A T2T / J(L)=2 LAM2 * L / 0 (2*Sin(THETA(L))) J L2 (2*Sin(THETA(L))) 3 NELRI(L)=4os2(THETA(L) )/Sin(THETA(L) )+Cos2(THETA(L) )/ (THETA(L))] / BETA Heading ILMNO, ILTER, Initialize Plot LATPA(1)...(Y) Plot SPECNO, CENT1, Plot vs. NELRI(1). (Y) Abcissa IAT, ENT2, ^ Subroutine using an asterisk(*) - b EXPTIM SHRINK NU S(L ), L=,1 RAD(L) NUU 0 <57 L=, 1 U = 0 / GLE (L) DET - - S = S + NELRI(L) 'F T = T + LATPA(L) T (6 L = 1, 1 U=U+NELRI(L) *LATPA (L) DELTA = 1 R \ L>Q MV= V+NELRI(L)2 LGOOD= (S*U-T*V)/(S'-V T o, 'L GLE(L-< 30.. SLOPE= (U-ST/Q)/(V-S~ Q W O ~ =_ ___~ __ ________.__________, \ _y ^My EPS W=W+X2 F VAR = rE/(Q -l) L O Z=l 9L Z y( ETA ETZ~~~~~ZETA ^YK~SLOPE ^ SLOPE - ') L-1 ) ANGLE(L)< 30.) ETA Y Z=l ^ ~KAPP^~ f ~AL ~PH;~ATPA(L)=0. NELRI(L)=O. E34

Example Problem No. 91 MAD Program and Data J. V. GLUCK S164P 002 020 000 J. V. GLUCK S164P 002 020 000 $COMPILE MAD, EXECUTE, DUMP RCALCULATION OF LATTICE PARAMETER INTEGER L, Q, SPECNO, RADIAT, FILTER, NN, J, R, Z, FILMNO 1,Y DIMENSION NN(35), RAD(35), THETA(35), LATPA(35), NELRI(35), 1J(35), ANGLE(35), DUMMY(1539) READ DATA ALPHA READ FORMAT DATA2, FILMNO, SPECNO, RADIAT, FXPTIM, CENT1, 1CENT2, FILTER, Q, R THROUGH IOTA, FOR L=1,1,L.G.Q IOTA READ FORMAT DATA3, NN(L), J(L), RAD(L) Z=O Y=Q SHRINK=9/(CENT1-CENT2) THROUGH BETA, FOR L=1,1,L.G4Q ANGLE(L)=(CENT1-RAD(L))*SHRINK*10 THETA(L)=ANGLE(L)*0e01745329 WHENEVER J(L).E.1, LATPA(L)=(LAM1*3QRT.(NN(L)+O.O) /(2*SIN.( 1THETA(L))) WHENEVER J(L).E.2, LATPA(L)=(LAM2*SQRT.(NN(L)+0.O))/(2*SIN.( 1THETA(L))) BETA NELRI(L)=(COS.(THETA(L)).P.2/SIN.(THETA(L))+COS.(THETA(L)).P. 12/THETA(L))/2 EXECUTE PLOT1.(0,8,10,11,10) EXECUTE PLOT2.(DUMMY, 5.50600.,8*600,8.200) EXECUTE PLOT3.($*$, NELRI(1), LATPA(1), Y) PRINT FORMAT GTITLE, SPECNO EXECUTE PLOT4.(45,ORD) PRINT FORMAT ABS VECTOR VALUES GTITLE=$1H1S9,14H SPECIMEN NO C6*$ VECTOR VALUES ORD=$ LATTICE PARAME 1TER*$ VECTOR VALUES ABS=$SO,22H NELSON-RILEY FUNCTION*$ PRINT FORMAT HEAD, FILMNO, SPECNO, RADIAT, EXPTIM, FILTER, 1CENT1, CENT2, SHRINK THROUGH NU, FOR L=1i1,L.G.Q NU PRINT FORMAT ANS1, L, NN(L), J(L), RAD(L), ANGLE(L), LATPA(L) 1, NELRI(L) GAMMA S=O T=0 U=O V=O THROUGH DELTA, FOR L=1,1,L.G.Q S=S+NELRI(L) T=T+LATPA(L) U=U+NELRI(L)*LATPA(L) DELTA V=V+NELRI(L).P*2 WHENEVER Z.E.1, Q=R LGOOD=((S*U)-(T*V))/(S.P.2-(Q*V)) SLOPE=(U-(S*T/Q))/(V-(S.P.2/Q)) L=0 W=O EPS L=L+1 WHENEVER L.G.Y, TRANSFER TO MU WHENEVER ANGLE(L).Le30*0.AND. Z.E1, TRANSFER TO EPS X=(SLOPE*NELRI(L)+LGOOD-LATPA(L)) W=(W)+(X.P.2) TRANSFER TO EPS E35

Precision Lattice Parameter Determination for a Cubic Material MAD Program and Data (continued) MU VARtSQRT.(W/(Q-1)) WHENEVER Z.E.1, TRANSFER TO ZETA PRINT FORMAT ANS2, LGOOD, SLOPE, VAR WHENEVER Z.E.O TRANSFER TO KAPPA ZETA PRINT FORMAT ANS3, LGOOD9 SLOPE, VAR WHENEVER Z.E.1~ TRANSFER TO ALPHA KAPPA THROUGH ETA, FOR L=11L.*G.Q WHENEVER ANGLE(L).L.30.0 LATPA(L)=0 NELRI(L)=O ETA END OF CONDITIONAL Z=1 TRANSFER TO GAMMA VECTOR VALUES HEAD=$1H1S9,9H FILM NO I59S5513H SPECIMEN NO C6 1/S10,12H RADIATION C4,S3~15H EXPOSURE TIME F3.1,4H HRS/S10, 1 9H FILTER C3/S10,10H CENTER 1 F6.3,S3.10H CENTER 2 F6.3,3S3 116H SHRINK FACTOR= F6.5///S10,59H LINc NO NNJ RAD(L) 1 THETA LATPA NELRI*$ VECTOR VALUES DATA2=$I5 C6, C2, F3.1, 2F6.3, C3, 213*$ VECTOR VALUES DATA3=$52, I2i1H,,I1# F6.3*$ VECTOR VALUES ANS1=$S13,I2,S7,I2,1HI1,S5,F6.3,S49F6.39S4, 1F8.5,S3,F8.5*$ VECTOR VALUES ANS2=$1HOS5,54HFOR ALL LINES, THE EXTRAPOLATED 1LATTICE PARAMETER IS F8.5/S19. 41H THE SLOPE OF THE LEAST SQ 1UARES LINE IS F8.5/S79 53H THE STANDARD DEVIATION OF THE LA 1TTICE PARAMETER IS F8.5*$ VECTOR VALUES ANS3=$1HOS5,54HFOR THETA +30, THE EXTRAPOLATED 1LATTICE PARAMETER IS F8.5/S19, 41H THE SLOPE OF THE LEAST SQ 1UARES LINE IS F8.5/S7, 53t THE STANDARD DEVIATION OF THE LA 1TTICE PARAMETER IS F8.5*$ END OF PROGRAM $DATA LAM1=1.78890, LAM2=1.79279* 42662NZ15DFC06.033.57515.572FE 21 16 1 3,131.420 2 8,130.040 311,129.415 412,129.220 516,128*505 624,127.260 727,126.830 8329126.140 940,125.085 1043,124.695 1144,124.565 1248,124.045 1459,122.580 1564,121.875 1772,120.640 1872,220.585 1975,120.110 2075,220.050 2176119.935 2280s119.110 2380,219.035 E36

Example Problem No. 91 Computer Output The following is the plotted output produced with the PLOT. subroutine along with a title and label for the abscissa produced using PRINT FORMAT statements. The output has been modified slightly (in an obvious fashion) to simplify its reproduction in this report. S P EC I M, N J:W [JZ 1 Z 'L;F........... 1.M~ -- - - - -..................~- --------:M ---- --- ------- " "y " ~ - - ~- ~.' 0~+~~- - i -- -- r —l,7 — ----- - -- - --- ------------------------------- -------- I ' I I I I I I I I I I I I I I I I I I i' 3,4%- +^ --- —--— + ----^ ----v --- —--— + --- —------- ---- -- - ----- I I I I I I I I I -L...... I..I I I I I I I I TFI _ _ I I I _ _ __I _I ____I I__. I A I I I I I I I I I '1..I. 4, +.I --- I I I I _ I I _I I_ T I I I I " I I I I M 8. 4____ +-_ -___ _ - -------------------- ------ -- -------- --- ------ --------- -------------- A I _ I I I I I II A __________I I__________ I I I I I I I E _ I I I I I I I I I..... I ** I I I I I I I I r 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 8.3L>0 --- + -- -- - -------— ` --- --- -+ —~ — --- --------- -----— I ---+ --- — -------------- ------- I I I I I I I I I I I I I I I I I 8-.- 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 -------- 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 ------------ -- -— j --- ~~- - - ---— ~ — - ~~- - -— y --- —Y --- —------------------ -------- ---- ------- -- -------- I I I I I I I I _I I I I I I I I I 8.200 ++ + ---+ --- —--— + --- —---- ___ 0 - 0.5 0 1. D- ( I I 0 ----— ^ --- --- -50 --- —- 1 0 3 —5. 5 iNFLSCN-RILi-Y FUNCTIGN E37

Precision Lattice Parameter Determination for a Cubic Material Computer Output (continued) The following is the printed output from the computer other than that associated with the plot. FIL —i Nfli 4"26 2 - SfELIM:.. f3i NC N,71 - [)F RAiIA I 1I CN Cti -:XPt tJIt-: [IF 6.O hS F IL TE R FL CEN TER 1 33. '75 CNI Ei 2 15.1572 SHRINK FgCTOR=.49992 LINE NOT N',J, A L) [HER [A L Tr P NE L I 1 3,l 31.420 10.773 8.28813 5.14172 2;, 1 30.040 17.672 8.33383 2.90699 3 If,1 29.415 20.797 8.35531 2.43461 4 1,1.29.220 21. 771 8.35383 2.29) 44 5 1 28.05 25.346 8.3 57 1. t8I1 ___________6_ _ __ _ 2 7.260 31.57C0 8.369_,0 1.3 2C0C 7 2, 1 26.33C 33.719 8.31234 1.2l109 8 32,1 2 6.14C 37.169 8.37482 1.149 '9 4, 1 25. 85 42.443 8.38253.7 7104 _ _ _ _. _-_-_ _ _.-.-.-.-..- - - _-_-_-_-, _ _ ___ _ _ _ _ ___ __ _ _ _ _____ __ _ _ _ _ _ _ _ _ _ _ 0 _, 2 4.-6 9 _5 4 4. 3 3 8. 8 __4 14 4__. _ _4 4 4 11 44,1 24.~565 45.042 8.38447.7C031 12 4, 1 24.045 47.642 8. 8613.58I13 13 51,1 22.580 54.966 8 8.39012.3/30( 14 6,1 21.875 58.490 8,.39316.29399 15, 1 20.640 64.664 8.39 136. 624 16 7,2 2C. J585 64. 93 9 08.:3- 664. 1781t L~ 7...7......,1 20.110 - 67.314. 3954.14392 18 Zb,2 20.750 67.614,.3t / 3.139F 19 76, 1 19. 935 68. 1o9 3. 39c;9 1 2 3% ------------- --- --— -- -L --- —0 —76- - --- -- -^*^ --- —------ ------- 20 ut 1 19.110 72.313 8.3 714.8.,8C1 21 8., 2 19.035 2. 688 8.3905.5u212, FCR ALL LINFS, THi EXTRitR,''L[TFDF) LATFICE PARl4i,"FT-i~, IS C.39945 lhE SLOPFE O;F TlTH LLAS1 SLItJAE_'-; L5 IS -.02130 THE Si AN)AR:L brIAT;T IC; 'F I —!E LTATiTIC PiA/ AL.I IS.,00231 FOK THE A +3C, Tfl L X iKA PELA I E) LATIICF PARA4'FFR IS.40 01 1-F SLOPE c-F I-3r LEAST SQUJAiFS LINt IS -.023C9 ThE STAN\DAll); EvL fIA [, 1 F- THE L/ ' I 'IC E PAR.A.FT-' IS ~. 0' 9~ Discussion of Results A program of the type developed is mainly useful where a number of samples are to be examined of a known material. In effect, this spreads out the necessary set-up time, line identification and determination of plane indices over a number of specimens. As an example, this program has been used to follow solubility effects in materials containing an excess constituent as affected by annealing atmosphere, time, and temperature. Another application would be the study of minor phase precipitates in heat-resistant alloys subjected to elevatedtemperature creep-exposure. Accurate cell measurements can also be used for such things as determining thermal expansion coefficients, phase boundaries in equilibrium diagrams, density measurement, and measurement of internal stresses (Reference 4). The accuracy of the parameter determination depends (ultimately) on the resolution of the doublets in the back-reflection region of the Debye-Scherrer film. Various authorities E38

Example Problem No. 91 differ somewhat on the degree of accuracy attainable. The following tabulation indicates two opinions: Accuracy Attainable - Per Cent Doublet Resolution Reference 3 Reference 4 1. Well resolved, sharp 0.005 or better 0.001 2. Resolved, but not sharp 0.02 0.2 - 0.02 3. Not resolved 0.1 4. High-angle lines too blurred to measure 1.0 1.0 In the example given, assuming a 2r scatter-band (95% confidence) and using only angles greater than 30~ in the extrapolation, ao = 8.4000 + 0.0020 A~ (rounding off to four places). This corresponds to an accuracy of 0.024 per cent, a figure which appears to be reasonable for the resolved doublets in the film used in the example. (The material was a nickel-zinc ferrite.) The plot function in this program has the main purpose of indicating which points have an excessive deviation from the least-squares extrapolation line. Since the experimental procedures are constant from film to film, it can be assumed that such a deviant point resulted from either an error in measurement or in indexing. It is then a simple matter to identify the suspect data, remove the card from the set of data cards, revise the identification card, and rerun the calculation. Alternatively, the program could be revised to perform a third extrapolation calculation, eliminating from consideration all points falling more than, say, 20 from the line calculated in the second extrapolation. The program reproduced above required 0.44 minutes for compilation and then executed five sets of data in 0.99 minutes. It is estimated that hand computation using a calculating machine and hand plotting of the results would require about three hours for each set of data. Much more time would be saved when examining the more geometrically complex structures. For example, a comparable program could be written for the determination of the two parameters of hexagonal materials by the method of successive approximations, a process that might require as many as five trials before consistent results were obtained (Reference 4, page 482). References 1. Klug, H. P. and L. E. Alexander, X-Ray Diffraction Procedures, John Wiley and Sons, Inc., New York, 1954. 2. Moroney, M. J., Facts from Figures, Penguin Books, Baltimore, Maryland, 1956. 3. Azaroff, L. V. and M. J. Buerger, The Powder Method in X-Ray Crystallography, McGraw-Hill Book Company, New York, 1955, p. 215. 4. Henry, N. F. M., Lipson, H. and W. A. Wooster., The Interpretation of X-Ray Diffraction Photographs, Macmillan and Company, London, 1960, pp. 191, 195. 5. Piazza, J. R., High Temperature Phase Equilibria in the System Carbon-Oxygen-Uranium, Ph.D. Thesis, The University of Michigan, 1961. E39

Example Problem No. 92 IONIC CRYSTAL STRUCTURE by L. H. Van Vlack Department of Chemical and Metallurgical Engineering The University of Michigan Course: Physical Ceramics Credit hours: 4 Level: Senior The following problem involves calculation of the Madelung Constant for an AX compound. It may be assigned in one of two ways: 1. The student may be required to decide on a suitable computational technique, program the solution, and check the program using appropriate test data. 2. The student may be required only to analyze the crystal structure for one or more ionic crystals and prepare the data in a form suitable for submission to a program pre-written by the instructor. This description of the problem and its solution are presented in detail sufficient for the first method of assignment (i.e., the student writes the program). In practice, however, the problem has only been assigned in the second manner (i.e., the student Just prepares the data) The advantages of the second method of assignment are discussed in the critique. Statement of Problem The Madelung Constant, Am, is a ratio between the potential energy, Ec, of an ion within a crystal lattice as compared with the energy, Epr, between a single pair of ions; thus Am = Ec/Epr The calculation of the Madelung Constant is achieved through the summation of coulombic attractions and repulsions. In strict mathematical terms the calculation may be achieved as m 2s q+ q- rij '7 However, any methodical procedure for selecting the three-dimensional lattice gives a series which converges extremely slowly, if at all, unless appropriate mathematical manipulations are made. (1) Select a procedure and write a program for the solution of the Madelung Constant. (2) Provide appropriate data from the crystal structure to solve the Madelung Constant, checking the result against known answers. (3) Discuss any limitations in the program imposed by the method selected for the program. E40

Solution The mathematical manipulation followed in the technique which is outlined below involves establishing a sphere surrounding a reference ion and calculating the attraction and repulsion forces between the reference ion and each ion within the selected radius; more remote ions are discarded. This method of calculation is valid for any AX compound regardless of crystal structure. It may also be used to determine the Madelung constant for ions on or in the surface, edge, or corner of the crystals with ((100)) faces. As presently programmed, the calculation can not be made for non-AX compounds or for other crystal surfaces*. The program could be modified to do so; however, it is not recommended that such an assignment be given to the ordinary student because of time limitations. The present program is sufficient to emphasize the order and pattern of NaCl-, ZnS-, ZnO-, CsCl-, and FeS-type structures. If desired, it may also be used to show the origin of surface energies. The student must supply data as follows which pertain to the crystal structure: 1. The crystal system and type. 2. The number of ions per unit cell. 5. The lattice constant and the reference dimension which is the closest approach of ions. 4. The ion charges, and their location within the unit cell. In addition he must put limits on his calculation as follows: 5. The sphere size for calculation (in angstroms). 6. The crystal block (in unit cell dimensions) which is large enough to enclose the sphere, but which avoids excessive calculations outside the sphere. Additional information may be supplied, stating the amount of intermediate printing which is desired as follows: 7. Whether intermediate parameters calculated for each cell, for only the cell with X = 0, Y = 0, and Z = 0, or no intermediate parameters are to be printed. The above required data are illustrated in the table below where MgO is used as an example. The data are given in a form which may be read by the computer program using the simplifiedinput (READ DATA) statement in the MAD language. The definitions of the symbols are given in the next section. * However, other crystal surfaces may be obtained by selecting non-conventional unit cells of these AX crystal types, e.g., the primitive cell of NaCl will give the ((111)) faces of the regular cell. E41

Ionic Crystal Structure EXAMPLE DATA SUPPLIED TO THE PROGRAM FOR MgO Crystal Properties STRUCT = $CUBIC$, TYPE = $MGO$, NION = 8 A = 4.203 REFDIM = 2.1015 CHARGE(O) = -2, U(O) = 0.0, V(O) = 0.0, W(O) = 0.0, CHARGE(1) = +2, U(1) = 0.5, V(1) = 0.0, W(1) = 0.0, CHARGE(2) = -2, U(2) = 0.5, V(2) = 0.0, W(2) = 0.5, CHARGE(3) = +2, U(3) = 0.0, V(3) = 0.5, W(3) = 0.0, CHARGE(4) = -2, U(4) = 0.5, V(4) = 0.5, W(4) = 0.0, CHARGE(5) = +2, U(5) = 0.0, V(5) = 0.0, W(5) = 0.5, CHARGE(6) = -2, U(6) = 0.0, V(6) = 0.5, W(6) = 0.5, CHARGE(7) = +2, U(7) = 0.5, V(7) = 0.5, W(7) = 0.5, Limits RLIMIT = 8.405 UPPERX = +2, UPPERY = +2, UPPERZ = +2, LOWERX = -2, LOWERY = -2, LOWERZ = -2, Print-out Options for Intermediate Calculations TRIAL = $CENTER$* List of Symbols The following are the important symbols used in the MAD program with their definitions. ADDON CHARGE(ION)/R3DIM, angstroms A unit cell dimension, angstroms BETA axial angle between Y and Z axes, degrees B unit cell dimension, angstroms CELLCK sum from cell for X = 0, Y = 0, Z = 0, normalized to REFDIM = 1, also equals (CHECK2-CHECK1)/REFDIM CHARGE(ION) array containing charge carried by each ion (ION) CHECK1 sum prior to beginning calculation for any cell CHECK2 sum after completing calculation for any cell COSBET cosine of BETA COSGAM cosine of GAMMA C unit cell dimension, angstroms GAMMA axial angle between X and Y axes, degrees ION counter for each ion within the unit cell E42

Example Problem No. 92 LOWERX lower limit of cell scans in X-direction, Y-direction, and Z-direction, LOWERY LOWERZ respectively MADEL Madelung Constant NEGIN negative charges within sphere* NEGOUT negative charges discarded (outside sphere) NION number of ions per unit cell POSIN positive charges within sphere* POSOUT positive charges discarded (outside sphere) R2DIM X-Y plane distance between reference ion and particular ion being considered, angstroms R3DIM three-dimensional distance between reference ion and particular ion being considered, angstroms RADIM A*(U(ION)+X) RBDIM B*(V(ION)+Y) RCDIM C*(W(ION)+Z) REFCHG charge of reference ion, in electron units REFDIM closest approach of two ions, angstroms RLIMIT radius of calculation sphere (It is best if this radius equals the distance to the closest surface ion.) SUM variable for storing cumulative results of integration of ADDON STRUCT variable containing alphabetic code for crystal structure: STRUCT = $CUBIC$, $TETRAG$, $ORTHO$, $HEXAG$, $RHOMB$, $MONOCL$, or $TRICL$, depending on whether structure is cubic, tetragonal, orthorhombic, hexagonal, rhombic, monoclinic, or triclinic, respectively. TRIAL variable containing alphabetic code to control intermediate printing: TYPE = $EACH1$, $CENTER$, or $NONE$ depending on whether intermediate printing is desired for each cell, only for the center cell, or not at all, respectively. TYPE variable containing alphabetic code for type of crystal, examples are $CSCL$ for CsC1, $NACL$, for NaCl, and $FES$ for FeS UPPERX upper limit of cell scans in X-direction, Y-direction, and Z-direction, UPPERY UPPERZ respectively U(ION) fractional cell location in X-direction, Y-direction, and Z-direction, V(ION W(ION, respectively of each ion (ION), X indices in matrix of cells in X-direction, Y-direction, and Z-direction, Z respectively * if POSIN f NEGIN, the program places the balancing number of ions on the surface of the sphere. This necessary approximation may place some variation in the second decimal place of the Madelung constant if there is a big difference between POSIN and NEGIN. E43

Ionic Crystal Structure SUM = 0 "Madlung CHECKl=O. Constant eginning RIAT= $NONE$ STRUCT and necessary CHECK2=0. Calcula- > Eof Program -data depending on POSOUT=O tion0 F I STCT POSIN= NEGOUT=0 NEGIN=0O START (U $B STRUCT = $TETRCUBIC$ STRUCT = $ORRUCT = HEXAG T T r T BETA =9M0. GAM24A = 90. GAMMA = 90. GAMMA = 120. BETA BETA=90. BETA= 90. BETA = 90. C A BETA = 90. BETA = 90. BETA = 90.: AB=A B B=A "Stru"Str"Strutc- S ture = tre =ture = "Struc- Cubic" 3 Tetrag- 3 rtho- 3< ture = onal't ^ rhombict Hexagl_^ l _^~ 0 1 ~onal _^^ _ _ _ _ _ _ _[,_ ~_. _______F_ ___"C he ck M TR ICL3 Structur C $RHOMB$s.STRUCT = $MONOCL$CA STRUCT = $TRaICL$l C M START T T T GAMMA = GAMMA = 90. All ParameBETA ters were B = A read in for "Struc- ItStruc- "StrucRhombic - | Mono- - 3 = Tri- 3 clinic"t Clinicl TYPE 0 /OSGAM C^OSGAMI<0.0 1 COSGAM 0 C CosBE.T9 RCosGmbic5.29F UPE "Data useO BETA /LOCAT \ CHARGE (ION) F~ -- -— ere-" GAMMA / A \ I (N ) COSBET<O. COSBET 0 A, B, C, COSGAM ION =, 1 ION -— NION COSBET ION > W ION -- "^ /^~^ \(NION-1)-J l U- _XPPERX LOWERX>UPPERY L- LATS ZH RTZC UPER_ LOWERY ATTX LAT (6 Z=LOWERZ, 1 TRIAL = $EACH1$ vTRIAL = $CENTER$ \Z>UPPERZ / ^Ax=oAy=oAZ=-o_ J E44

Example Problem No. 92 C ELL=, RADIM=A* ('U(ION)+X R2DIM = IM+BDIM2+2*RBDIM*RADIM*COSGAM ION =0O,1 _3 RBDIM=B* kV IONR2Dj M =-8 ION- RCDIM=C*(W ION +Z r3DIM = CDIM2R2DI+2RCDIM*R2DIM*OSBET (la:ON-1) r IF F T 1- -NEGIN = R3DIM >RLIMIT ADDON 0 11 NEGIN - REFCHG ' ADDON = CHARGE(ION)/R3DIM 11\ ADDON=OACHARGE(ION F DN AcHARGE(N F ADDON 13 AR3DIM >RLIMIT AR3DIM >RLIMIT ^ CHARGE(ION) >0 T iT POSOUT = NEGOUT = POSIN = POSOUT + CHARGE(IO 1 NEGOUT-CHARGE(IO POSIN+CHARGE(ION) ADDON / F SUM = RIAL = $EACH1$ V(TRIAL = T 3 A CHARGE (ION)< 10 SUM + ADDON ' $CENTER$AX=OAY=OAZ=0) 14 NEGIN = /\ NEGIN-CHARGE(ION) - |R3DIM IR3DIM 1 C~-LX(R3DIA0EM>RK2MIT (C CK'C E =I 0'.0 m ADDON ADDON VR3DIMR E SUf ADDON CELLCK = |F CELL X = OA Y = O A Z = OHECK2=SUMX (CHECK2-CHECKOl)* I SUM /'^ '\ /'^ '\ /'^ '\ P~OSOUT MADEL = l GLATTZ UATTY ' LATTX POSIN CELLCK (SUM+( NEGIN-POSIN)/ MADEL STAR ' — /.........~NGouT ' /'LIMIT)* START NEGIN I G ~ REFDIM/(-REFCHG) E45

Ionic Crystal Structure MAD Program L. H. VAN VLACK S225B 002 050 000 MADLUNG __ L. H VAN VLACK......___. S225B 002 050 Q O0 MADLUNG L. H. VAN VLACK S225B 002 050 000 MADLUNG -_ $COMPILE MAD EXECUTE DUMPPRINT OBJECT, PUNCH OBJECT INTEGER X,YZ, IONtNION,UPPERXUPPERYUPPERZLOWERX,LOWERY. 1 LOWERZSTRUCTTYPEHOLDTRIALPOSOUT POSINNEGOUT.NEGIN DIMENSION CHARGE(100),U(100),V(100),W(100) STABZ I R. ----..IAL = $NQONES -. _ READ DATA SUM=O CHECK1=0. CHECK2=0. POSOUT = 0 POSIN = 0 NEGOUT = 0 NEGIN = 0 PRINT COMMENT $1 ***** *** 1 *************** $ PRINT COMM4ENT $ * MADELUNG CONSTAN __ _1____ LT-. CALCULAT ION * —. -- —.PHERE$ PRINT COMMENT $ ****** 1************** $ WHENEVER STRUCT.E.$CUBIC$ GAMMA=90. BETA=90. B=A..._ _. C=A PRINT COMMENTSO STRUCTURE = CUBIC$ OR WHENEVER STRUCT*E.$TETRAG$ GAMMA=90. BETA=90. ___..................B=A................ PRINT COMMENT$O STRUCTURE = TETRAGONAL$ OR WHENEVER STRUCT.E.$ORTHO$ GAMMA=90. BETA=90. PRINT COMMENTSO STRUCTURE = ORTHORHOMBIC$..... OR WHENEVER STRUCT.E.$HEXAG$.___ GAMMA=120, BETA=90. B=A PRINT COMMENTSO STRUCTURE = HEXAGONAL$ OR WHENEVER STRUCT.E.$RHOMB$ GAMMA=BETA B=A C=A PRINT COMMENT$O STRUCTURE = RHOMBIC$ OR WHENEVER STRUCT.,E.$MONOCL$ GAMMA=90. PRINT COMMENTSO STRUCTURE = MONOCLINIC$ OR WHENEVER STRUCT.E.$TRICL$ PRINT COMMENTSO STRUCTURE = TRICLINIC$ OTHERWISE PRINT COMMENTSOCHECK STRUCTURE LABELS TRANSFER TO START END OF CONDITIONAL PRINT FORMAT XTAL, TYPE VECTOR VALUES XTAL = $1HO,S10, 7HTYPE = C6*$ PRINT COMMENTS..-..-. -$ $ COSGAM = COS,(GAMMA/57,2958) WHENEVER *ABS. COSGAM,Le 0.001, COSGAM = 0.0 COSBET = COS.(BETA /57,2958) WHENEVER *ABS. COSBET Le 0.U001, COSBET = 0.0 PRINT COMMENT$ODATA USED WERE-$ PRINT RESULTS A,BC,NION PRINT RESULTS BETA, GAMMA,COSGAM, COSBET THROUGH LOCATEr FOR ION=0,1,IONG.(NION-1) LOCATE PRINT RESULTS CHARGE(ION),U(ION),V(ION),W(ION) REFCHG=CHARGE (0) PRINT COMMENTS-REFERENCE CHARGE IS CHARGE OF ION AT U(O), V(O 1), W(O)$ PRINT RESULTS REFCHG E46

Example Problem No. 92 PRINT COMMENTS-REFERENCE DIMENSION IS THE CLOSEST APPROACH OF MAD Program (continued) 1 POSITIVE AND NEGATIVE IONS$ MAD Program (continued) PRINT FORMAT RFORM, $REFDIM$, REFDIM VECTOR VALUES RFORM = $1HuS8,C6,3H =,F12.6,10H ANGSTROMS *$ PRINT COMMENT S-LIMITS USED WERE-$ PRINT RESULTS UPPERXUPPERYUPPERZ PRINT RESULTS LOWERX#LOWERYLOWERZ PRINT COMMENT $ THE UPPER AND LOWER 1LIMITS ARE IN UNIT CELL DIMENSIONS$ PRINT FORMAT RFORM, $RLIMIT$, RLIMIT PRINT COMMENT $- $ THROUGH LATTX, FOR X=LOWERX,1,X.CUPPERX THROUGH LATTY, FOR Y=LOWERY,.1Y.G.UPPERY THROUGH LATTZ, FOR Z=LOWERZl.1Z.G*UPPERZ WHENEVER TRIAL.Eo$EACH1$.OR*( TRIAL.E.$CENTER$.AND. 1X.E.O.ANDY.E*O.AND*Z.E.O) PRINT RESULTS X, Y, Z PRINT COMMENT $ $ END OF CONDITIONAL THROUGH CELL. FOR ION=O,1,ION.G.(NION-1) RADIM=A*(U(ION)+X) RBDIM=B*(V(ION)+Y) RCDIM=C*(W(ION)+Z)' R2DIM=SQRT (RADIM.P.2+RBDIM* P2+2*RADIM*RBDIM*COSGAM) R3DIM=SQRT~(RCDIM.P. 2+R2DIM2+2*RCDIM*R2DIM*COSBET) WHENEVER R3DIM.E0.0 ADDON = (p.0 WHENEVER REFCHG.G.0,POSIN=POSIN+REFCHG WHENEVER REFCHG*L*0,NEGIN=NEGIN-REFCHG OR WHENEVER *ABSs R3DIM *G. RLIMIT ADDON = 0.0 OTHERWISE ADDON = CHARGE(ION)/R3DIM END OF CONDITIONAL WHENEVER ADDON.E.O.AND.CHARGE(ION)*G.O.AND.R3DIM.G.RLIMIT POSOUT=POSOUT+CHARGE(ION) OR WHENEVER ADDON.E.O.AND.CHARGE(ION).L.0.AND,R3DIM.G.RLIMIT NEGOUT=NEGOUT-CHARGE(ION) OR WHENEVER ADDON.NE.QoAND.CHARGE(ION).*.O POSIN=POSIN+CHARGE(ION) OR WHENEVER ADDON.NEO.QAND.CHARGE(ION).L.0 NEGIN=NEGIN-CHARGE(ION) END OF CONDITIONAL SUM = SUM + ADDON WHENEVER TRIAL.E.$EACH1$*OR.( TRIAL.E.$CENTER$.AND. 1X*.EO.AND.Y.E.O*AND.Z.E.O) WHENEVER R3DIM.G.RLIMIT PRINT COMMENT $+ OUTSIDE$ OR WHENEVER R3DIM.E.O.O PRINT COMMENT $+ SKIPPED$ END OF CONDITIONAL PRINT FORMAT ANSALL, R3DIM, SUM, ADDON VECTOR VALUES ANSALL = $S6, 8HR3DIM = F12.6,S22,6HSUM = 1F12.6,S16,12HINCREMENT = F12.6*$ END OF CONDITIONAL CELL CONTINUE WHENEVER X.E.O.AND.Y.E.O.AND.Z.E.O CHECK2=S5M CELLCK = (CHECK2-CHECK1)*REFDIM END OF CONDITIONAL CHECK1=SUM LATTZ CONTINUE LATTY CONTINUE LATTX CONTINUE PRINT RESULTS POSOUT, POSING NEGOUTt NEGIN PRINT FORMAT CFORM, $CELLCK$, CELLCK VECTOR VALUES CFORM = $1HO, S8, C6, 3H = t F12.6, S9, 1 78HCELLCK IS THE SUM FROM CELL FOR X = 0, Y =, Z 0 NOR 2MALIZED TO REFDIM = 1 *$ MADEL=(SUM+(NEGIN-POSIN)/RLIMIT)*REFDIM/(-REFCHG) PRINT COMMENT $0 ******************** 1 *************$ PRINT FORMAT ANSWER, MADEL VECTOR VALUES ANSWER = $47H+ MADELUN 1G CONSTANT =, F93-*$ PRINT COMMENT $ * 1 *$ PRINT COMMENT $ ******************* 1-*************$ TRANSFER TO START END OF PROGRAM E Eh7

Ionic Crystal Structure Data $DATA STRUCT=$CUBIC$,TYPE=$MGO $,NION=8.... TRIAL=$NONE$ A=4.203 REFDIM=221015..U....... PPERX=+2,LOWERX=-2.UPPERY=+2tLOWERY-=2-. --- UPPERZ=+2,LOWERZ=-2 RLIMIT = 8.405 CHARGE(0)=-2, U(0)=0.0, V(O)=O.Ot W(O)=0.0 CHARGE(1)=+2, U(1)=0.5, V(1)=0.0, W(1)=O.O CHARGE (2)=-2, U(2)=0.5, V(2)=0.0, W(2)=0.5 CHARGE (3)=+2, U(3)=00, V(3)=0.5, W(3)=OQ —.-0.-... CHARGE (4)=-2, U(4)=0.5, V(4)=0.5, W(4)=00 CHARGE (5)=+2, U(5)=0.0, V(5)=0.,0 W(5)=0.5 CHARGE (6)=-2, U(6)=00, V(6)=0.5, W(6)=0.5 CHARGE (7)=+2, U(7)=0.5, V(7)=0.5. W(7)=0.5 * STRUCT=$CUBIC$,TYPE=$BLENDE$, NION=8.......... A-= —A=5.41.2 RLIMIT = 10.82 REFDIM = 2.35 UPPERX=+2,LOWERX=-2,UPPERY=+2,LOWERY=-2..__...__.... UPPERZ=+2,LOWERZ=-2 CHARGE(O)=-2, U(O)=O.O, V(O)=OO, W(O)=OO CHARGE(1)=+2, U(1)=0.75tV(1)=0.75,W(1)=0.75 CHARGE (2)=-2, U(2)=0.5, V(2)=0.0 W(2)=0.5 CHARGE (3)=+2, U(3)=0.25, V(3)=0.25, W(3)=0.75 CHARGE (4)=-2. U(4)=0.5, V(4)=0.5, W(4)=0.0 CHARGE (5)=+2, U(5)0=.75,V(5)=0.25W( 5)=Q0.25 _._.... _ CHARGE (6)=-29 U(6)=0.0, V(6)=0.5, W(6)=0.5 CHARGE (7)=+2, U(7)=0.25tV(7)=0.75,W(7)=0.25 STRUCT=$HEXAG$o TYPE=$WURTZI$, NION=4 A=3.811 I C=6.234 RLIMIT = 11.25 REFDIM = 2.34 UPPERX=+3,LOWERX=-3,UPPERY=+3,LOWERY=-3 UPPERZ=+2PLOWERZ=-3 CHARGE(0)=-2, U(0)=0.0, V(0)=0.0, W(0)=0O0 CHARGE (1)=+2, U(1)=0.Q *V(1)=00, W(1)=0.375 CHARGE (2)=-2, U(2)=0*333333,V(2)=0.66667,W(2)=0.5 CHARGE (3)=+2, U(3)=0.333333,V(3)=0.6666667,W(3)=0.875.* STRUCT=$HEXAG$, TYPE=$FES$, NION=4 A = 3.535 C = 5.774 RLIMIT = 12.00 REFDIM=2.5 UPPERX=+3 LOWERX=-4, UPPERY=+3, LOWERY-.-4 UPPERZ=+2,LOWERZ=-3 CHARGE(0)=+2, U(0)=0.,V(0)=0. W(0O)=0. CHARGE(1)=-2# U(1)=0.666667,V(1)=0,333333,W(1)=0.25 CHARGE(2)=+2, U(2)=0333333,V(2)=0.666667,W(2)=0.50 CHARGE(3)=-2, U(3)=0.666667,V(3)=0.333333,W(3)=0.75 * _. STRUCT=$CUBIC$,TYPE=$CSCL$,NION=2....__.. A=4.110 REFDIM =3*559 RLIMIT = 8.21 UPPERX=+2,LOWERX=-2,UPPERY=+2,LOWERY=-2 UPPERZ=+2,LOWERZ=-2 CHARGE(O)=-1l U(O)=0,0 V(O)=OO W(O)=00.O. CHARGE(1)=+1i U(1)=0.5, V(1)=0.5# W(1)=05 * Computer Output A typical set of computer output is shown on the next page. This corresponds to the first set of data shown with the MAD program above. An example was selected with no intermediate printing to conserve space in this report. E48

Example Problem No. 92 Computer Output ********************** ************ * MADELUNG CONSTANT CALCULATION * SPHERE ---- * ** ** * * * ** * * * *** ** * **** * ** * ***** * * * * STRUCTURE = CUBIC TYPE = MGO DATA USED WEREA = 4.2030 B = 4.203U, C = 4 2030, NION = 8 - -- BETA = 90.0000, GAMMA = 90.0000, COSGAM =.J000, COSBET =.0000 CHARGE(O) = -2.0000, U(O) =,OO000 V(O) = G00OO, W(O) =.0000 CHARGE(1) = 2.000, U(1) =.5000, V(1) =.0000 W(1) =.0000 CHARGE(2) = -2.0000, U(2) = *5000, V(2) =.0000 W(2) =.5000 CHARGE(3) = 2.0000, U(3) =.0000 V(3) =.5000, W(3) =.0000 CHARGE(4) - -2000, U(4) =.50U0 "1(4) =.5000, W(4) =.0000 ----. CHARGE(5) = 2.0000, U(5) =.000, V(5) = *000, W(5) =.5000 CHARGE(6) = -2.0000, U(6) =.0000, V(6) = *5000, 4(6) =.5000 CHARGE(7) = 2.0000, U(7) =.5000, V(7) =.5000, W(7) =.5000 REFERENCE CHARGE IS THE CHARGE OF ION AT U(O), V(),t W(0) REFCHG = -2.000 REFERENCE DIMENSION IS THE CLOSEST APPROACH OF POSITIVE AND NEGATIVE IONS REFDIM = 2.1015 ANGSTROMS LIMITS USED ----E-RE — UPPERX = 2, UPPERY = 2, UPPERZ = 2 LOWERX -2 LW = -2, LOWRY = -2 LOWERZ = -2 THE UPPER AND LOWER LIMITS ARE IN UNIT CELL DIMENSIONS RLIMIT = 8.4050 ANGSTROMS POSOUT = 768, POSIN = 232. NEGOUT = 730 CELLCK = 2.912, NEGIN = 270 CELLCK IS THE SUM FROM CELL FOR X = U, Y = 0, Z = NORMALIZED TO REFOI --- —------ * MADELUNG. CONSTANT = 1.742 - - - E 9- -X- - -X - - - - i E49

Ionic Crystal Structure Discussion of Results Results have been obtained for several different crystal structures. The numerical values of the Madelung Constants which were calculated by the program are listed below along with the generally accepted value where available. Madelung Constant Generally Accepted Crystal Structure Calculated by Program Value* NaCl 1.744 1.747 CaCl 1.77 1.1762 ZnS (Blends) 1.633 1.658 ZnS (Wurtzite) 1.567 1.641 FeS 1.745 MgO 1.742 1.747 It can be seen that the computational technique utilized by the program give values of the Madelung Constant which is generally accurate to about two decimal places. Critique Experience in several classes has indicated that the calculation of the Madelung Constant serves as a very useful teaching procedure. Successful results have been obtained by requiring the student only to prepare appropriate data for submission to a program pre-written by the instructor. This avoids having the students spend an undue amount of time on the problem. In order to supply the correct data the student must understand the crystal structure of his problem more thoroughly than he would from simply reading descriptions in a textbook. v The problem solution then serves to check him and, in fact, provides him with a means of analyzing his structure for errors if he receives an incorrect answer. * Sherman, Chemical Review, 11, 93 (1952) E50

Example Problem No. 93 UNSTEADY-STATE HEAT CONDUCTION IN SOLIDIFYING ALLOY by J. R. Street and J. 0. Wilkes Department of Chemical and Metallurgical Engineering The University of Michigan Course: Rate Operations Credit hours: 5 Level: Junior Statement of Problem A cylindrical mold is to be charged with a molten alloy at an initial temperature of 4000F., the liquidus temperature of the alloy. The mold is well insulated except at one end, which is maintained at the solidus temperature of the alloy, 150~F., by a stream of oil. The geometry is indicated - - Oil at 150~ F. - - in the figure at the right. - By (a) setting up the partial differential equation and appropriate boundary and initial conditions to des- Aloy Mold /Mold cribe the temperature e(x,t) of the solidifying ingot as a function of distance x and time t, and (b) reduc- / = L 2 ft. x = L t=. It ing the partial differential equations and associated // Mold conditions to appropriate finite difference forms, calculate the temperature profiles in the ingot as a function of time. Also, determine the time required for the face x = 2 ft. to cool to 200~F. The following data are available for the alloy: Temperature Pliquid = Psolid = 540 lb./cu.ft. t kliquid = ksolid= 10.1 B.Th.U./hr.ft.~F. (c ) = (c ) = 0.038 B.Th.U./lb.~F. Liquid P liquid P solid 4000 - 400~F Liquid Liquid AH = heat of fusion of the alloy = 120.0 B.Th.U./lb. + Solid 150~F The quantity a( ), given in the table on the Solid following page as a function of temperature G, has Composition Alloy been determined from the phase equilibrium diagram Composition shown in the figure at the right. It gives the amount of solid formed corresponding to a small drop in temperature of the liquid/solid mixture. E51

Unsteady-State Heat Conduction in Solidifying Alloy Values of a(G) = lbs. solid formed per ~F. per lb. mixture of solid and liquid. (OF) -a(G) G(OF) -a(9) 9 (~F) -a( ) ((OF) -a( ) 150 0.0119 220 0.00454 290 0.00252 560 0.00173 160 0.0102 250 0.00412 300 0.00258 370 0.00167 170 0.00874 240 0.00375 310 0.00223 380 0.00160 180 0.00756 250 0.00341 320 0.00211 390 0.00153 190 0.00657 260 0.00315 330 0.00200 400 0.00147 200 0.00575 270 0.00291 340 0.00190 210 0.00509 280 0.00270 350 0.00181 Solution By performing a heat balance on a differential element, and taking into account a generation term arising from latent heat effects, thedifferential equation giving the alloy temperature Q(x, t) as a function of distance x and time t is found to be = P (c p-aH) (1) - k p < where a = a(G) is a known function of temperature, there being no heat conduction in the y direction. Equation (1) is subject to the boundary and initial conditions, G(O,t) = 150 ~F for t 7 0 (la) G-(Lt) = 0 for t o 0, at L = 2.0 ft. (lb) Q(x,0) = 400 ~F for 0 x L (1c) In the finite difference method of solution, the interval 0<x~ L is divided into a grid of points spaced Ax apart, and temperatures at all grid points are computed at timeintervals t. Replacement of /At.ax2 and 3G/3t by their finite difference approximations yields the following finite difference representation of equation (1); G(x+Ax,tiAt) - 20(x,t+At) + 9(x-Ax,t+At) = k (cp - aAH ) [G(x,t+At)-G(x,t)] (2) The choice of approximation of 92Q/9x2 from temperatures not at time t, but at time t+At, has led to an "implicit" finite difference form. In general, the G(x,t) will be known and the G(x,t-At) will have to be calculated from equation (2) applied at each grid point. This requires the solution of a set of N simultaneous equations (where N = number of grid points), but does not suffer from the restriction on the size of At which attends the alternative "explicit" form. Let I designate the grid point with x coordinate = IAx; I may have any value from 0,1,2...to N Then, in order to simplify programming, the number of subscripts on Q may be reduced from two to one by introducing the variables: S(I) = value of 9 at grid point I and time level t T(I) = value of 9 at grid point I and time level t+at E52

Example Problem No. 93 With this new notation, and substituting T(I) = (Cp -c a LH) t (3) Equation (2) becomes T(I) T(I+l) + T(I-1) +Tr(I)S(I) (4) =T(T(I) + (4) In equation (4), the "old" temperature S(I) will be known at all grid points I = 0,1,...,N. From a table look-up procedure the corresponding values of a and hence the values of P(I) can be found. Successively better approximations for the "new" temperatures T(I) are obtained by repeated application of equation (4). As a starting guess it is convenient to put T(I)=S(I). The iteration is repeated until the T(I) have converged sufficiently; a suitable criterion for stopping is when the:um n of the absolute values of the deviations of the T(I) from their last computed value is less than some arbitrarily small quantity max. At this stage the T(I) are then called "old" temperatures and are substituted into the S(I) and the process is carried on for as many time-steps as is desired. In the present case, computation was stopped when T(N) had dropped to 200 ~F. Boundary condition (la) dictates that S(O) = T(O) = 150 for all values of t. The initial condition (1c) sets S(I) = 400 for I=1,2,...,N at t = 0. Boundary condition (lb) requires that, instead of using equation (4), the temperature T(N) at the last grid point N should be computed fromequation (5): 2T(N -1) + P (N)S(N) (5) T(N) =N +2 List of Symbols MAD Symbol Definition T(I) "New" temperature at point I S(I) "Old" temperature at point I J iteration parameter for table look-up Y(J) table look-up values of temperature X(J) table look-up values of -a A -a GAM(I) value of F(I) = P (c - amH) corresponding to S(I) B P (Ax)2 k At DELT At (time increment) HEURE total time from start EPS sum of absolute values of deviations of the T(I) from their last computed values E53

Unsteady-State Heat Conduction in Solidifying Alloy EPSMAX x = largest tolerable value of EPS max K counter on the number of iterations required to compute T to required accuracy at any one time level. N index of last grid point (number of grid points = N+1) I index for grid points (0 ---IN) Flow Diagram (START |Beginnin ad Table o Print Initialize oR Values of H ead N,DELT Heading __S(O), T(O), \ ^ I _ e_ and _a(G) EPSMAX URE, and()..T(N) Put Find Value of Calculate new T S(1).. S(N) GAM(I for temps. T(N) Print n T HEURE S (l)NTN)+ eELT TB(N) -EPSKMSM? temp. (I) BE and uGIN table-lokup errors. MAD Program and Data Jo 0O WILKES X137N OC5 040 000 JOW08 —. J*- 0 WILKES X137N 005 040 000 JOW08 S COMPILE MAD, EXECUTE. _.. __. __._..-.. SQL DIEICATION QF AN INGOT.... INTEGER Is Je Ks N JOW08010 DIMENSION S(30)# T(30)# GAM(30), X(30)t Y(30) JOW08020 READ FORMAT TABLEt X(Q.O)*X(26), y(O)e..Y(26) JOW08025 VECTOR VALUES TABLE = S(8F10.5)*S JOW08026 START READ FORMAT IN, DELT. EPSMAX N JOW08030 VECTOR VALUES IN = $2F10.5I10*S JOW08040 PRINT COMMENT $1 UNSTEADY STATE HEAT CONDUCTIO N NSOLIDJOW08041 1IFYING INGOT WITH THE PARAMETERS$ JOW08042 PRINT RESULTS DELT, EPSMAX# N JOW08045 B = (540.0/101~)*(2*0/)(2/N)(2 )/DELT JOWO8050 HEURE = 0.0 JOW08055 __ O_) 150 _ _ _ JOW08060 S<O) = 150.0 JOW08070 THROUGH ONE, FOR I 1.i 1. IeGeN JOW08080 ONE T(I) = 400.0 JOW08090 BEGIN. THRQ.UGH THREE, FOR I = 1, 1, I.GoN JOW08100 S(I) ~ T(I) JOW08110 TWO._ THROUGH TWO. FOR J = 0, 1t Y(J).G.S(I) JOW08120 E54

Example Problem No. 93 MAD Program and Data (continued) A = X(J-1) + (X(J)-X(J-1))*(S(I)-y(J-1))/IYJ)-Y(J -1) JOW08130 THREE_ GAM I) _ B*(_0_038 + A*12 0 ) __ ___JOW08O 40____ EPS = 10.0 JOW081a50 THROUGH FIVE, _FOR K = 1, 1 EPS.L*EPSMAX JOW08160 EPS = 0 0 JOW08t7O THROUGH FOUR. FOR I = 1 1, I*G*N-1 _ JOW08180 DUM = T(l) JOW0&190 _- -T___ T} li{i1tl) +_ T (I-1) +_ GAM(I)*s(I/(rAM(Ij + 2.0) JOWOB200 FOUR EPS = EPS +.ABS.(T(1) - DUM) JW0O8210 __ DUM T(N) JO_ V.08220 T(N) (2.0*T(N-1) + GAM(N)*S(N))/(GAM(N) + 2*0) JOW08230 FIVE EPS = EPS + *ABS.(TLN - DUM W_ JO0Q82O HEURE = HEURE + DELT JOW08245 PRINT FORMAT OUT. HEURE, K, T(0).T(N) I__ JOW08250 VECTOR VALUES OUT $SS10F10*5 I1O/(551lF8.3)*S JOW08260 WHENEVER T(N).*G200*0. TRANSFER TO BEGIN _ JOWO&270 TRANSFER TO START JOW08280 END OF PROGRAM _ __.L __JOWOf290 S DATA. 0.0119_ _. O 102 0.00874 0. 756 0.00657.00575 0.00509 0.00454 0.00412 O 375 0*00341 0. 315 0.00291.00270 0.00252 0.00238 *__0_223_ O.0_._ 211_ 0.00200 0. 190 0 00181,00173 0 00167 0 o00160 0.00153 0. 147 0.00140 150. 160.0 170.0 180.0 1900 2.0j 0 a210. 220.0 2230. 240.0 2500 260 2600 270.0 280.0 290. 300.0 310. 320.0 330 0 340.0 350 0..360_.0 370. _ 380.0 390. 400 0 4100* 4.0 0.2 20 Computer Output A typical set of results obtained using the data shown with the MAD Program above is presented below. UNSTLADY SAITEl HFAI (L; i0iUCI [fN I 'lt SQLIUIIFYIN; Il, g I tF TioE P-: RA\MLAfRS iEL = 4.0.0 0 FPSFAX =.200C l 6 = 4. 00( C CL 13 -15C.LOO 188.633 221.,J11 249.273 272.611 292.291 3u8.8 1' 7 322.843 334.. 8 34. 44 449 352. 707 359 5)8 3b5 319 370. 04 373. 77 376.959 379.367 361.170 382.421 3~3 1:'~4 3~3 392.6 *. -.0, 1 0 1 I 1sc.00 1(74. 4" 1,)/.o b 1.72 237.230 254.4.3 27.028 263.t 296.435 37.4 41 317.17 325.6 1 333 3.C 33-.434 344. 86 349.210 352.757 355.470 357.382 358.517 358 ob9 215~C I-tH 1-c6, l-(.9;7 74. 36 20 2720.010 3 4, -2 748.448 260.9 tC '72.97 292.876. 251 3u0.4L4 30.- 71.5 -14.163 319.657 324.261, 327.9045 130.878 332.9,5 _34.i46 37A43 ^ 7 1 50. 0 0 165. 0 16 ".:~ 4.)2 I'. O9 942 15 4 9.213 2 126.1 1 34.116 2450.254 2 5. 5 1 6 2, 9 3. 44 1 281.122 2. 9 fl 59 2 3. 92 2)9.143 303. 5 9 317.C6 3U9.827 311. 1791 31. 96)5 313 lg.01 7.163 5.2 l2- j9 41 6 219.461 2 1.528 1 2.18 223.6 0 233.. 2 242. 66 21 '36 ^5o 9o5 265.948 2 2. 190 7 7.6 1 2. 452 286. 73 29. 757 292.37 2t4. 12-5 29-5 12 2 5.69,. - - _ _ _ _ _ _4 _ _ _ _ _ _ _ __. 3_ _ _ _ _ _ _ __ _ _ 3 __ _ _ _ _ - - - - - - -- 1 (.0CC 161.943 173.66 1 4.O 195 642 25_.87 215._48 224,365 232. 240, 42 247. 34J 2 — 3.6)4 259.3/3 264.3,5 2 68.129 2/2.40t 275.411 2I. 74 9 279.418 2'LC.416 28C6 144 210000 1.60.674 171 12 I 14} Iu 2 li5 190 937 20. 1 9. 2u.797 216.915 224.467 231 438 237.t;l 32..0O;, 4 u 24 3.60 6 24.7 2 245..73 c39 3. 7 L _ 7. 341 26o77W4 -3.44.-65.2 2&7.17 2I68. 9 26.18 0 L, & 15 0.0 u 1539.636 16.1091 1(6. 263 l 1.54 195.415 tL3.31.3 210. I(1 21.L2 22 4.iuC 2.6 44 235.147 239.901 244.105 247 153 20 6.42 2o3.371 213.337 256.739 9 )77.07 2 _7 7 27.4f[ 1 3 '-.? 0] 0 ' — ' -- 10 54, C 1 ) 8 *(I69 1 6 7. *3 8 175.1 3 1 3 *7)-0 1.445 198.65 2i5.43 2 11 02 2 17. 7. 37l. 221.922 32 335 36.t29 239.72 2 2.421 244.7 54 246.569 247. 64 24 8.cj 248 1 4 i. 3 Ci IJ 52C 1,0 CoOiO 1 8.i30 lb)29. )32 13.611 loh VO. 94 168i3.045 194.716 200.41 206.640 212 )59,17 j4 -221- 54 ~22.82 229.4LZ 232.l18 235. 16 237.318 238.998 240.17 6 242-821 2, 41'1".1 44 Ju_ b1 -C 0.u Z GI __l 3 r) 4 4. 67 76 19 178. 3 1 5. 100 1 91.2 7 197. 69 202 17 207 46 1 7 64 21 6.34 220.137 223.476,.26..77 2;8.3. 2. 471 32.412 233.528 2 4. 15 2 34,2 15.oG C 15 6.38 I..7 1707. 1 41 165. 4(1 1.*3 1O.3.5 13.o l3 1 7 25 203.410 7 (14 211.629215.148 218.2u4 220. 373 223.269 225. 150 22o.6i3 227.656 22AH.O 8:4o3 E55

Unsteady-State Heat Conduction in Solidifying Alloy Computer Output (continued) 5. O 0 C 4 ' 15C.GOC 16. 34c' 1'62.Slu 168. 713 1 /4.t6i& 16i,.245 i05.oiK) 190.61 1,).366 1)9. (45 2u3_ (69 207.431 210U-723 213.6 38 216.173 218.323 22.3 5 2 1.4 'o 222.433 i 223. L7 - -223.2 5- O 0. 'O 4.. 130.0CC 135. 1f6 1 2 1 61. 1 167.446 172. )49 1 8.2'1 lo3.231 17. - 2 b 16.s4 2 0.2.9 _203.695 2J6.64 2C9=.521 2 11.39_9 213.916 215 569 26lSi) 211.713 t3 1_____ 6..`J 0i. ', J: 150*.0`0 15.1 I)o 986. 114 I_1.4_0 1 -.1C i. 110 18. 41 15) 9. t i 3. 6' 1 Uj1 6-.. J,.>; 3 c. 200.342 203.249 2''.825 2(,(. j64 2,)9.'64 211.519 212.719 213.592:141. i1 214.275 64 * t; 00 ~ l.0 ~ 0`.o.ll 11 JII.29)b 1T665 97 ITO, 1 3T) l 7' t. ' 1(9 _20L 103. J7 o 87.2 9 0 10.: 92 4 194.I 2-2 191.320 230.003 2U2.494 4 204.6C7 2;',6.399, 2C/.48 209.01 209.625 219,j12 210.4/i5 6. C (i.(. CO 34 15U.C0O 154.83U 159 61/ 1,)4.378 168.932 1/3.307 1(7.414 11.414 151;8 l8.,j_41 1 ---._-4.. 194.5035 197.1/6 1)-).47/ 201.477 2C3.172 2G04.562 20(3.642 206.412 2o6.672 '21?..71 12.00(_00 _ _ _ __33 _ _ ____ _ ----_ ---— _ - 150.COO 154.' 83 1 l.12 13.' 43 1t67.8a38 1(1.966 175.'9(2 179. 625 163.118 1-;6.367 lc9.j361 192.090 194. 546 196.724 198.619) 2cUu.221 2(1.544 2*2,b60 23. 3C:, 2 3.-736 2,;, l7.77 76.30ii, 31 150.000 1 143 " j 3IJ. 13 h ' 70 V (4 1,6.q42 1(C.46 1(.4.411 17. 9) 1u1.3'6 1o4.&6 167.'''6 -ll~o~uj).? ------------- -------------- * -8^ -! -_-_________ -~ZJ_______________l~___ _ ___1 4_ (-_- -__ * — L -- __ 189.815 192.147 1)4.2i5 196.14 1~/7.540 I1 ).192 l19.765 200.459 2..br3 21 1.L;7 80. CUCb 3, 15-.O0 154.u-J8 1;8.13C 1b2.3;9 lo5.:) 30 169.628 173.158 176.t02 179.6,43 1i32.567 1,t,.2o5 187.125 l189.942 l)1.9:)8 1w3.519 1:95.071 19'6.260 197.165 197.345 198.239 19)8.366__ Discussion of Results In order to investigate the effect of choice of grid spacing, time increment and allowable error on the computed temperatures, the program was run for several different selections of N, At and E ma. The output printed above corresponds to one such choice, viz N=20, at = 4.0 hours, and ~max = 0.20, and the resulting temperatures at x = 0,1 and 2 feet are plotted against time in the figure below. The results for the remaining runs are summarized in the table below. 400 Variation in Temperature at Three Points in the Ingot: x=0, x=1.0 and x=2.0 feet Results are for N=20, at=4 hours, EPSMAX=0.2 300 Temperature of Ingot 200 0 20 40 60 Time, Hours E56

Example Problem No. 93 Effect of choice of N, at, and max on the time t (hours) max for the face x=2.0 feet to cool to 200~F Effect of time increment Effect of allowable Error Effect of grid spacing (For N=20, Emax=0.2) (N=20, at=l.0) (at=l.0, Emax=0.2) At t emax t N t 4.0 77.52 1.0 84.73 5 75.21 2.0 76.97 0.5 79.49 10 75.50 1.0 76.66 0.2 76.66 20 76.66 0.5 76.42 0.1 75.80 30 77.97 0.25 76.07 0.05 75.30 0.1 75.90 0.02 75.06 0.01 74.98 It is seen from the table that there is, perhaps, not so much need to make At small as to ensure, by reducing 6max that the sequence of iterations has converged sufficiently at each time step. The average numbers of iterations required for max values of 1.0, 0.1 and 0.01 max were 16, 28 and 41 respectively. The IBM 709 computer performed approximately 20 such iterations (over 20 grid points) per second. The slight divergence of t as N is increased has not been explained. Critique Although the present example, together with some notes on finite difference methods, was assigned at the start of the semester, many students waited until the 11th week, when heat transfer was first discussed in lectures, to start work on the problem. Half the class succeeded in obtaining correct answers to what amounted to a fairly difficult problem at the Junior level. The problem well illustrated the ability of the computer to handle an otherwise formidable task. E57

Example Problem No. 94 COOLING OF PIG IRON IN A TRANSFER LADLE by Robert D. Pehlke Department of Chemical and Metallurgical Engineering The University of Michigan Course: Rate Operations Credit Hours: 5 Level: Junior The pig iron produced in an iron blast furnace is transferred to steelmaking facilities in a ladle mounted on a railroad car. Several types of ladles are used, the most prominent one being the Pugh-type which is a cigar-shaped tilting ladle. The pig iron is tapped from the blast furnace into the ladle and then transferred to the steelmaking facilities which may be some distance from the blast furnace plant. It is of interest then to know the time which one can allow the metal to remain in the ladle without any solidification taking place. Statement of Problem Write and test a MAD program which will permit a calculation of the temperature drop of liquid pig iron held in a closed transport ladle of the Pugh-type (See Figure 1). Using the test data, carry the computations to either the initial solidification temperature or to an elapsed time of twenty-four hours. 150 Ton BOTTLE CAR Figure 1. Suggested Approach: Assume the following: 1. The temperature of 'the liquid metal is uniform and equal to the internal interface temperature of the refractory lining of the ladle car. 2. The heat losses may be estimated by determining the rate of loss over short intervals of time and approximating the heat loss as the rate-time product. 3. Neglect the thermal resistance of the thin steel shell which surrounds the ladle. Problem Preparation: Do the following in the given order: 1. Prepare a flow diagram for the problem in the form of concise algorithmic statements. 2. Convert the flow diagram to a set of MAD statements. 3. Keypunch the statements and required data cards. 4. Flow diagram, MAD statements, and keypunched cards should be submitted. E58

Data and Nomenclature: The program should be tested on a ladle car used for pig iron transport under the following conditions: TM Temperature of charged metal, 2700~F TA Temperature of surroundings, 80~F AO Area of external surface, 800 ft2 AI Area of internal surface, 530 ft2 K Thermal conductivity of refractory, 2.35 BTU-ft/hr-ft2-~R L Thickness of refractory, 1.5 ft HC Convective heat transfer coefficient, 0.28 BTU/~R0'25 E Emissivity of steel shell, 0.8 W Weight of metal charged, 150 tons F Factor to account for specific heat of ladle, 1.4 EPS Tolerable error in surfact temperature, 0.1 OR I Time interval, 0.5 hour J Number of time intervals IMAX Maximum time of holding, 24 hours NMAX Maximum number of iterations to determine surface temperature, 100 TLOW Temperature of initial solidification, 2100~F CP Specific heat of metal, 0.2 BTU/lb OR R Radiation from surroundings to the ladle AM Log mean area for heat transfer through the refractory CPM Gross heat capacity for ladle and metal, BTU/lb OF TO Temperature of metal at time, P, OR P Elapsed time from charging of ladle, hours TS Temperature of outside surface, OR N Number of iterations required for temperature convergence QC Conduction heat transfer, BTU/hr QR Radiation and convection heat transfer, BTU/hr TI Most recent value of TS during iteration, OR TN Metal temperature after iteration, i.e., latest TO, OR TX TN in OF TY TI in OF Solution Heat is transferred to the surroundings by conduction through the refractory and then by radiation and convection to the ambient surroundings at temperature TA. The conduction heat transfer during a unit time is QC = K.AM (TTI) where the cylindrical wall is being considered; the area for heat conduction is given by the relationship AO -AI AM AOAO AI The radiation loss during unit time is Q= 0.173 X 10-8 (E) (TI4-TA4)AO The convection loss is Qc = h AO (TI-TA) E59

Cooling of Pig Iron in a Transfer Ladle Here the convection coefficient h is a function of the temperature difference to the one-quarter power, (TI-TA)0'25. Therefore, one may write Q =HC AO (TI-TA)125 The combined convection and radiation loss is QR = r + The heat loss from the metal results in a decrease in its temperature. This decrease in energy is AE = 2000 W.CP-F (TO-TN) By energy balance QC = QR If either of these is multiplied by the time interval, I (it must be short enough so that the changes in the temperature during the interval are small), a quantity of heat is obtained which determines the decrease in metal temperature during the interval. Thus, I-QC = 2000 W-CP-F (TO-TN) = CPM (TO-TN) where CPM is the combination of a number of constants. At the start of any time interval the metal temperature TO and the ambient temperature TA are known. The surface temperature TI is notknown, so a first guess is made. This allows a QC and QR to be calculated. If TI has been guessed correctly, QC comes out equal to QR, and then the temperature of the metal at the end of the short time interval may be calculated as TN = TO - I CPM If TI has not been determined within 0.10R, it is necessary to correct the assumed temperature. The magnitude of the change to be made in TI depends upon how different QC and QR turn out to be for a given TI. In this solution it was decided that if ITS - TII| EPS where EPS is 0.1~R, the surface temperature would be changed to TI = TI + (100) (QQRQR) It is clear that a higher order approximation might be made in the heat transfer equations by putting in the average metal or average surface temperatures during the time interval. Thus, for conduction, one might write QC = KAM (TO+TN - TI ) L 2 Tave) The anticipated change in temperature, however, did not merit this refinement, so it was simply assumed for the heat transfer calculations that TO and TI remained constant during the short time interval. E60

Example Problem No. 94 Flow Diagram Cooling of Liquid-Metal Transport Ladle \ PH]a (IPRINT oalc:R,Al|, P=O Start) READ INPUT - TITLE OPM - N=0 --- |TM=TM+460- __EHROUGH I\ AP PRINT TA=TA+460 QC=O /FOR J=l1, (PRd FRESULT TS=TA QR=1 TO.LE.TLOW 4) ) TM,P,N TO=TM \.OR.P.GE, / TI=(TM+TA )/ 2MAX / IHROUGH G. OR N=,1,N\ QR TI=TI+ f h G.NMAX.OR. - TS=TI -calc: Q ( 100)*(Q- (GAMMA \ABS.(TS-TI)/ IL-QR)/QR. L.EPS / P INT FORMA RESULT, TPRINT PN TM TN —TO- TX=TN-460. RESULT, TS TA P=I *J (QC*I)/ TY=TI-460. TX1PN — TO2TN T CPM TYQ QR (TRANS ( ~ COH I T- TO J \ END v APHA / MAD Program and Data $COMPILE MAD, EXECUTEDUMP ALPHA READ FORMAT INPUT, TM, TA, AI, AO, Ko L, HCP E, W9 CP, 1 2F, EPS9 It IMAX, NMAX, TLOW 2 PRINT FORMAT TITLE1i TM, TA, AI, AO K, L, HC, Et Wt CP 3 PRINT FORMAT TITLE2 F, EPSN I, IMAXR NMAX, TLOW 4 R = (0.173E-8)*E*(TA*P.4) 5 AM = (AO-AI)/(+ELOG(AO/AI)) 6 CPM = W*2000.*CP*F 7 P = 0 8 N = 0 9 PRINT FORMAT RESULT, TMX P, N 10 TM = TM + 460. 11 TA = TA+ 460. 12 TS = TA 13 TO = TM 14 TI = (TM+TA)/2 15 QC = 000 16 QR = 10 17 THROUGH CHI# FOR J=l,i> TOeLE.(TLOW+460*.)ORP.eGE*IMAX 18 THROUGH GAMMA, FOR N=1v1v NeG.NMAX.OReeABSe(TS-TI)eLeEPS 19 TS = TI 20 QR= AO*(HC*((TI-TA-)*P.(125) }+(0173E-8)*E*(TI.P*4)-R) 21 QC = K*AM*(TO-TI)/L 22 GAMMA TI = TI + (100)*(QC -QR)/QR 23 P = I*J 24 TN = TO -(QC*I)/CPM 25 TX - TN - 460. 26 TY = TI - 460. 27 PRINT FORMAT RESULT, TX, P, N~ TY, QCt QR 28 TS = TA 29 CHI TO = TN 30 INTEGER J. N, NMAX 31 TRANSFER TO ALPHA 32 VECTOR VALUES INPUT = $,45,0,4F5.2,F4.0,4F4,2,F4,0I3,F50*$ 33 E6l

MAD Program and Data, Continued VECTOR VALUES TITLE1= $48H1TEMPERATURE DROP OF LIQUID METAL I 34 2N A LADLE CAR //5H DATA/34H INITIAL METAL TEMPERATURE Ft TM = 35 3F7.1/ 25H AIR TEMPERATURE Ft TA = F6.1/25H INSIDE AREA SQ FT 36 4, AI = F51l/ 26H OUTSIDE AREA SQ FT, AO = F5.1/ 27H THERMAL C 37 50NDUCTIVITYt K = F4l1/ 30H REFRACTORY THICKNESS FT, L = F4.1/ 38 630H CONVECTION COEFFICIENT, HC = F4.1/ 23H SHELL EMISSIVI 39 7TYt E = F4.1/ 24H METAL WEIGHT TONS, W = F5.1/ 34H SPE 40 8CIFIC HEAT METAL BTU/LB, CP = F4.1*$ 41 VECTOR VALUES TITLE2 = $18H HEAT FACTOR, F = F4.1/ 32H TcMPER 42 2ATURE DEVIATION Ft EPS = F4.2/24H TIME INCREMENT HR, I = F4.1 43 3/ 37H MAXIMUM NUMBER OF INCREMENTS, IMAX = F5,1/ 34H MAXIMUM 44 4NUMBER OF TRIALS. NMAX = I4/ 26H FREEZING POINT F, TLOW = F6. 45 51//// 30H METAL TEMPERATURE VERSUS TIME ///7H METALS5t4HTIM 46 6ES4t10HITERATIONStS4,7HSURFACE5S4,10HCONDUCTIONS4t8HRAD-CON 47 7V/8H TEMP FSS552HHRS19t6HTEMP FS4,11HLOSS BTU/HRtS3,11HLOS 48 8S BTU/HR//*$ 49 VECTOR VALUES RESULT = $S1,F8.2,F6.1,IllF14.2tF14.2,F13.2*$ 50 END OF PROGRAM 51 $DATA 2700 80 530 800 240 150 30 80 150 18 140 10 50 48100 2100 D-1 Computer Output TEMPERATURE DROP OF LIQUID METAL IN A LADLE CAR DATA INITIAL METAL TEMPERATURE Ft TM = 2700.0 AIR TEMPERATURE F. TA = 80.0 INSIDE AREA SQ FT. AI = 530,0 OUTSIDE AREA SQ FT. AO = 800*0 THERMAL CONDUCTIVITY, K = 2.4 REFRACTORY THICKNESS FTT L = 1*5 CONVECTION COEFFICIENT# HC = *3 SHELL EMISSIVITYt E = *8 METAL WEIGHT TONS. W = 150,0 SPECIFIC HEAT METAL BTU/LB, CP =.2 HEAT FACTOR, F = lf4 TEMPERATURE DEVIATION Ft EPS = *10 TIME INCREMENT HR. I =,5 MAXIMUM NUMBER OF INCREMENTS, IMAX = 48.0 MAXIMUM NUMBER OF TRIALS. NMAX = 100 FREEZING POINT F, TLOW = 2100.0 METAL TEMPERATURE VERSUS TIME METAL TIME ITERATIONS SURFACE CONDUCTION RAD-CONV TEMP F HR TEMP F LOSS BTU/HR LOSS BTU/HR 2700.00,0 0 2685.61.5 25 626.47 2175488.53 2177499.22 2671.31 1.0 7 624.59 2162391,78 2163937.19 2657.10 1.5 7 622.74 2149327.75 2150826.00 2642*97 2.0 7 620.90 2136346.56 2137822.84 2628*92 2.5 7 619.06 2123450.69 2124906.56 2614*96 3.0 7 617.22 2110639,69 2112075.47 2601.09 3.5 7 615.40 2097913*03 2099328*97 2587*30 4.0 7 613.57 2085270*08 2086666.47 2573.59 4.5 7 611.75 2072710.31 2074087.20 2559.96 5.0 7 609.93 2060233 08 2061590.89 2546.42 5,5 7 608.12 2047837.87 2049176.80 E62

Example Problem No. 94 Computer Output, Continued 2532.96 6.0 7 606.32 2035524.12 2036844 19 2519*57 6*5 7 604.51 2023291.25 2024592.89 2506.27 7.0 7 602.71 2011138.69 2012422.09 2493.05 7e5 7 600.92 1999065*89 2000331.17 2479.91 8,0 7 599.13 1987072.33 1988319.81 2466*85 8.5 7 597,34 1975157.45 1976387.30 2453.86 9.0 7 595.56 1963320.70 1964533.14 2440695 9.5 7 593*79 1951561.56 1952756.69 2428.12 10.0 7 592.01 1939879.47 1941057.56 2415.37 10.5 6 590.31 1928169.75 1930086.48 2402.69 11.0 7 588.48 1916739.84 1917921.75 2390.09 11.5 6 586.78 1905187.78 1907058.81 2377.57 12.0 7 584.97 1893907*50 1895055.41 2365.12 12.5 6 583,27 1882506.14 1884329.44 2352.74 13*0 6 581.53 1871269.09 1873123.59 2340*44 13.5 6 579.79 1860112*70 1861947.16 2328.21 14.0 6 578.05 1849029.78 1850840.55 2316.05 14.5 6 576.32 1838019.33 1839806.17 2303.97 15*0 6 574.59 1827080.78 1828844.14 2291.96 15.5 6 572.86 1816213.64 1817953.58 2280*02 16.0 6 571.14 1805417*45 1807134.25 2268*15 16.5 6 569*42 1794691.66 1796385.69 2256.35 17.0 6 567.71 1784035.87 1785707,27 2244.62 17.5 6 566.00 1773449.53 1775098.58 2232.96 18.0 6 564.29 1762932.20 1764559.14 2221.37 18.5 6 562.59 1752483.39 1754088.47 2209.85 19.0 6 560.90 1742102.62 1743686.16 2198*39 19.5 6 559.21 1731789.44 1733351.67 2187.01 20.0 6 557.52 1721543.39 1723084.47 2175.69 20.5 6 555.84 1711364.00 1712884.09 2164.44 21.0 6 554.16 1701250.78 1702750.19 2153.25 21.5 6 552.48 1691203.34 1692682.25 2142,13 22.0 6 550,81 1681221.19 1682679.97 2131.08 22.5 6 549.14 1671303*89 1672742.67 2120.09 23.0 6 547.48 1661450.97 1662869.97 2109.17 23.5 6 545.82 1651662.03 1653061.48 2098.31 24.0 6 544.17 1641936*62 1643316e69 Discussion of Results The result of the calculation is in keeping with observations from steel mill-operations, i.e., pig iron can be allowed to remain in a bottle car for periods ranging up to twenty-four hours. Critique This problem has not been used in teaching but appears to be well suited for illustrating the use of the computer in solving heat transfer problems of the type which commonly arise in equipment design. An extension of the problem is possible in terms of the cost of refractories which could be used as lining and insulating materials for the ladle. E63

Example Problem No. 95 DIGITAL COMPUTER ANALYSIS OF HEAT FLOW AND TEMPERATURE DISTRIBUTION AROUND A COPPER CONVERTER TUYERE1 by Robert D. Pehlke Department of Chemical and Metallurgical Engineering The University of Michigan Course: Rate Operations Credit Hours: 5 Level: Junior Statement of -Problem Write and test a MAD program which will perform the relaxation calculations outlined by Krivsky and Schuhmannl for determining the heat flow and temperature distribution around a copper converter tuyere. The cross-section of a copper converter through a tuyere center line is shown schematically in Figure 1. A relaxation technique was employed to obtain a solution and involved dividing the solid into equally-spaced grid points, specifying the boundary conditions, and assuming a steady state temperature distribution which dictates that the heat flow into each grid point shall be zero. The heat flow is idealized into a two-dimensional matrix with radial symmetry about the center line of the tuyere through which the blast air is supplied. The reference paper derives the equations for heat flow into each type of grid point which includes the interior points, interior or exterior surface points, tuyere surface points, and the pointsat the tuyere source and tuyere mouth. These heat flow equations are then solved in terms of the steady state temperature of the grid point itself and set into a series of conditional statements which permits the calculation of the steady state temperature at each grid point as a function of the temperature of the surrounding grid points or boundary temperatures. (It should be noted that in Figure 6 of the reference, the equation which reads Q4o0 = h4w(r + a ) (T4-T0) should read Q40 = h 4 Arr(r+ a) (T4-TO) 1 Krivsky, W. A., and R. Schuhmann, Jr. Heat Flow and Temperature Distribution Around a Copper Converter Tuyere," AIME Transactions, 215, February 1959, p. 82. E64

Figure 1....... Blast __ _ -- 1 1/2" tuyere (100~F) Ambient Refractory Charge air (2100~F) (1000F) I Idealized cross section of single 1 1/2 in. tuyere in 15 in. wall, showing grid of points for relaxation calculation. The Basic Approach: a. Initialize the temperature distribution by a programmed estimate of the temperature at each grid point. b. Make a point by point calculation of the steady state temperature by systematically moving through the grid, utilizing the surrounding temperatures and heat flow equations of each grid point. c. After calculating the temperature at each grid point, make a comparison of the calculated temperature with the temperature which existed on the previous pass through the grid. When there are no temperature adjustments which are greater than a specified temperature increment C for any of the points in the grid from one pass to the next, the calculation is considered complete and the temperature distribution is then printed out. Flow Diagram E65~^\ PRINT DIM T- START)-? READ input dat~ — DIM(2)= -^ title THICK & data A. / TH- ALPHA \ APHA \ nitial \ (o-( A=1,1 1(U=1,1,> rid ALPHA 9)+ ) / \ A.G.THICK / \U.G.DIST / T(A, U)= E65

Heat Flow and Temperature Distribution Around a Copper Converter Tuyere Flow Diagram, Continued Calculate / / TH- GAM!1A BBLAST, - N=O > BETA A=1,1 BMELT =O 0 A.G.THIOK / ABAIR H1'1 /-1 7- I ~.L.THICK A1-~U 'G' TN= TN=.V_ G71 --- U.E.1 F ~E,1 E /U.E:DIST T.1. THICK...THICK E.1 A...L.TICK r., - - N = TN= TN= TN= TN= \' r-', 'u.:ERR C=C+0 1 7,U)= 0+1 MADrg r ad DaPRINT t TRaAN 1 COMPILEMAD EXECE PanCH ECT D TAR MAD Program and Data $COMPILE MAD, EXECUTE, PUNCH OBJECT, PUNCH LIBRARY, DUMP TUYER DIMENSION T(38CODIM) 0 VECTOR VALUES DIM = 2, 1, 0 00 START READ FORMAT INPUT, DIST, THICK, NMAX, SCALE, RADIUS, TMELT, 1 2TBLAST, TAIR. HBLASTt HMELT# HAIR, KCOND, ERR 2 VECTOR VALUES INPUT = $3I 5F5S,3,F4F647F62,F4.2*$ 3 DIM(2) = THICK 00 PRINT FORMAT TITLE 4 VECTOR VALUES TITLE - $58H1TEMPERATURE DISTRIBUTION AROUND A 5 2COPPER CONVERTER TUYERE///*$ 6 PRINT FORMAT LIST, DIST. THICK, NMAX, SCALE, RADIUS. TMELT, 7 2TBLAST, TAIR, HBLASTe HMELTt HAIR, KCOND, ERR 8 VECTOR VALUES LIST =$5H DATA//33H DISTANCE FROM TUYERE IN SPA 9 ICES = I4/ 33H REFRACTORY THICKNESS IN SPACES = I4/ 31H MAXIMU 10 2M NUMBER OF ITERATIONS = I4/ 27H SCALE IN SPACES PER FOOT =F6 11 3,3/26H TUYERE RADIUS IN FEET = F5.4/ 41H MELT TEMPERATURE I 12 4N DEGREES FAHRENHEIT = F5.0/ 42H BLAST TEMPERATURE IN DEGREES 13 5 FAHRENHEIT = F5.0/40H AIR TEMPERATURE IN DEGREES FAHRENHEIT 14 6= F5.0/ 29H HBLAST. BTU PER HR-SQ FT-F = F5.1/ 28H HMELT. BT 15 7U PER HR-SQ FT-F: F5.1/ 27H HAIR, BTU PER HR-SQ FT-F = F5.1/ 16 E66

Example Problem No. 95 MAD Program and Data, Continued 860H THERMAL CONDUCTIVITY OF REFRACTORY, BTU-FT PER HR-SQ FT-F 17 9 = F5.1/ 21H ERROR IN DEGREES F = F4.2///*$ 18 THROUGH ALPHA, FOR A=1ltl A.G.THICK 19 THROUGH ALPHA, FOR U=1,1, U.G.DIST 20 ALPHA T(A*U) = TAIR + (TMELT - TAIR)*(A-1)/(THICK-1) 21 BBLAST = (HBLAST)/(KCOND*SCALE) 22 BMELT = (HMELT)/(KCOND*SCALE) 23 BAIR = (HAIR)/(KCOND*SCALE) 24 N = 0 25 BETA C = 0 26 THROUGH GAMMA, FOR A = 1,1. A.G.THICK 27 THROUGH GAMMA, FOR U = 1,1, U.G.DIST 28 R = RADIUS + (U-1)/SCALE 29 WHENEVER U.E.1lAND&A*G.1.AND.A.L.THICK 30 TN = ((2*BBLAST*TBLAST)+1++(l/(4*R*SCALE)))*(T(A-1,U)+T(A+1IU 31 2)+(2*T(AU+1))))/(2*BBLAST+4+(1/(R*SCALE))) 32 OR WHENEVER U.G.1.AND.A.E.1 33 TN= (T(AU-1)+T(A,U+1)+(2*T(A+1,U))+(2*BAIR*TAIR)+U(T(AU+1)- 34 2T(AtU-1))*(1/(2*R*SCALE))))/(4+(2*BAIR)) 35 OR WHENEVER U.Ge1*AND.A.EETHICK 36 TN = (T(AU-1)+T(AU+1)+(2*T(A-1,U))+(2*BMELT*TMELT)+(T(AU+1 37 2)-T(AU-1))*(1/(2*R AL/2SCALE) (4+(2*BMELT)) 38 OR WHENEVER Us.E1.AND*A.E*THICK 39 TN =((BBLAST*TBLAST)+(T(A-1,U)*(1+(1/(4*R*SCALE))))+(T(A,U+1) 40 2*(1+(1/(2*R*SCALE))))+(TMELT*BMELT*(1+(1/(4*R*SCALE)))))/(2+B 41 3BLAST+(1/(4*R*SCALE))+(1/(2*R*SCALE))+((BMELT)*(1+(1/(4*R*SCA 42 4LE))))) 43 OR WHENEVER U.E#1,AND.A.E.1 44 TN =((BBLAST*TBLAST)+(T(A+1,U)*(1+(1/(4*R*SCALE))))+(T(AtU+1) 45 2*(1+(1/(2*R*SCALE))))+(TAIR*BAR*(SCALE+(1/(4*R)))))/(2+BBLAS 46 3T+(1/(4*R*SCLE)) /(2CAL)+(1/(2*RSCALE)+(BAIR*(SCALE+(/(4*R))))) 47 OR WHENEVER U.E.DISTANDA.*G.1.AND.A.L.THICK 48 TN = T(1,DIST) + ((T(THICKDIST)-T(1,DIST))*(A-l)/(THICK-1)) 49 OTHERWISE 50 TN = (T(A-1,U)+T(A+1,U)+T(AU+1)+T(AU-1)+((T(AtU+1)-T(A~U-1) 51 2)/(2*R*SCALE)))/4 52 END OF CONDITIONAL 53 WHENEVER.ABS.(TN-T(A#U))*G.ERR C = C+1 54 GAMMA T(AU) = TN 55 N = N+1 56 WHENEVER *NOT. (N*.GNMAX.ORC.E.O)* TRANSFER TO BETA 57 PRINT FORMAT RESULT, N, C, T(1,1)..sT(THICKDIST) 58 TRANSFER TO START 59 INTEGER At U, THICK. DIST. NMAX, C. N 60 VECTOR VALUES RESULT = $4H N =I4, 4H C =I4/(10F10.2)*$ 61 END OF PROGRAM 62 $DATA 16 16 10012000 0625210000 10000 10000 20000 50000 600 200 50 D-9 16 16 10012000 0625210000 10000 10000 1000 50000 600 200 50 D-10 Computer Output TEMPERATURE DISTRIBUTION AROUND A COPPER CONVERTER TUYERE I) ATA DISTANCE FROM TUYERE IN SPACES = 16 REFRACTORY THICKNESS IN SPACES = 16 MAXIMUM NUMBER OF ITERATIONS = 10C SCALE IN SPACES PER FCOT =12.COO TUYERE RADIUS IN FEEl =.0625 MELT TEMPERATURE IN I)EGREES FAHRENHEIT = 2100 BLAST TEMPERATURE IN DEGREES FAHRENHEIT = 100 AIR FEMPERATURE IN DEGREES FAHRENHEII = 100 HBLAST, MTU PER HR-SQ FT-F =200.C HMELT, BTU PER HR-SQ FT-F_ =500.0 HAIR, BtiU PER HR-SQ FT-F = 6.0 THERMAL CCNDUCTIVIIY CF REFRACTORY, BTU-FT PER HR-SQ FT-F = 2.0 ERROR IN DEGREES F =.50 E67

Heat Flow and Temperature Distribution Around a Copper Converter Tuyere Computer Output, Continued N.= 58 C =0....__________________ 115.30 223.80 282.93 322.67 351.66 373.55 390.20 402.67 411.57 417.18 419.45 417.96 411.64 397.92 370.01 306.09 121.37 255.83 329.81 379.66 416.03 443.49 464.38 480.06 491.30 498.49 501.66 500.45 493.96 480.47 457.21 422.C01 126.09 289.32 379.02 439.43 483.47 516.70 542.01 561.08 574.92 584.10 588.86 589.13 584.51 574.43 558.48 537.93 130.98 324.66 430.98 502.46 554.50 593.72 623.59 646.19 662.80 674.27 681.08 683.48 681.58 675.51 665.77 653.85 136.17 362.21 486.15 569.30 629.66 675.03 709.55 735.72 755.17 769.02 778.02 782.73 783.60 781.14 776.11 769.77 141.72 402.42 545.12 640.55 709.53 761.15 800.29 829.95 852.14 868.26 879.34 886.19 889.47 889.90 888.30 885.69 147.72 445.84 608.63 717.00 794.84 852.68 896.26 929.17 953.81 971.91 984.75 993.32 998.48 1001.02 1001.77 1001.61 154.28 493.28 677.75 799.71 886.51 950.35 997.99 1033.68 1060.31 1079.91 1094.02 1103.78 1110.18 1114.05 1116.23 1117.53 161.57 545.98 754.05 890.20 985.81 1055.12 1106.12 1143.88 1171.80 1192.28 1207.03 1217.37 1224.32 1228.78 1231.56 1233.45 169.88 605.96 839.96 990.68 1094.46 1168.18 1221.43 1260.21 1288.50 1309.04 1323.73 1333.97 1340.80 1345.10 1347.69 1349.37 179.71 676.70 939.42 1104.41 1214.79 1291.04 1344.77 1383.10 1410.58 1430.24 1444.09 1453.54 1459.59 1463.05 1464.67 1465.29 192.01 764.72 1058.97 1236.20 1349.81 1425.39 1477.00 1512.90 1538.12 1555.83 1568.05 1576.08 1580.78 1582.76 1582.62 1581.21_ 208.98 883.56 1209.91 1392.94 1503.12 1572.83 1618.67 1649.66 1670.95 1685.61 1695.46 1701.59 1704.56 1704.63 1701.92 1697.13 239.07 1064.67 1411.-59 1583.52_____ 1678.04 1734.23 1769.64 1792.89 1808.53 1819.10 1826.00 1829.98 1831.23 1829.42 1823.59 1813.05 356.41 1389.48 1693.95 1816.03 1875.50 1908.51 1928.50 1941.31 1949.79 1955.44 1959.04 1960.94 1961.06 1958,53 1950.54 1928.98 1612.75 2060.73 2081.18 2086.95 2089.70 2091.22 2092.14 2092.73 2093.12 2093.38 2093.54 2093.63 2093.63 2093.49 2092.07 2044.90 TEMPERATURE DISTRIBUTION AROUND A COPPER CONVERTER TUYERE DATA -DISTANCE FROM -TUYERE- IN SPACES = 16 REFRACTORY THICKNESS IN SPACES = 16 MAXIMUM NUMBER OF ITERATIONS = ICC SCALE IN SPACES PER FCOT =12.CO0 TUYERE RADIUS IN FEET =.0625 MELT TEMPERATURE IN DEGREES FAHRENHEIT = 2100 BLAST TEMPERATURE IN DEGREES FAHRENHEIT = 100 AIR TEMPERATURE IN DEGREES FAHRENHEIT = 100 HBLAST, BTU PER HR-SQ FT-F = 10.C HMELT, BTU PER HR-SQ FT-F =500.0 HAIR, BTU PER hR-SQ FT-F = 6.0 THERMAL CONDUCTIVITY CF REFRACTORY, BTU-FT PER HR-SQ FT-F = 2.0 ERROR IN DEGREES F =.50 N =:. ___=.._........Q_......_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 187.54 324.10 374.13 401.84 420.29 433.60 443.47 450.79 456.00 459.35 460.86 _ 460.35 __._457.29 4_50.3_8 436.11 40_._13 3_._2 44,28 4 79,__9 501.81 518.22 530.49 539.60 546.14 550.38 552.42 552.11 549.02 542.28 530.50 512.60 377.84 470.96 524.26 560.18 586.05 605.31 619.83 630.71 638.60 643.90 646.79 647.26 645.17 640.23 632.30 622.07 436.75 543.59 603.32 643.72 673.08 695.09 711.78 724.37 733.64 740.09 744.00 745,56 744.86 742.00 737.29 731.54 493.56 616.72 684.48 730.00 763.02 787.77 806.57 820.81 831.40 838.96 843.93 846.63 847.30 846.26 843.93 841.02 551.17 691.70 768.31 819.36 856.15 883.61 904.40 920.15 931.94 940.51 946.43 950.10 951.91 952.26 9_ 51.59 _950.49 610.87 769.59 855.52 912.29 952.86 982.90 1005.50 1C22.56 1035.33 1044.70 1051.33 1055.73 1058.35 1059.61....10. ]0__59_. 98 1 59 6.__ 6.673.57 851.36 946.94 1009.47 1053.65 1086.01 1110, 12 112817 1141.63 1151.51 1158.56 1163.35 1166.38 1168.10 1168.96 1169.43 740.27 938.17 1043.61 1111.71 1159.13 1193.36 1218.52 1237.14 1250.90 _1260. 9'3 1268_Q5 7.l_? __ 8_... 1275.93 1277.64 1278.48 1278.90 812.36 1031.62 1146.95 1220.11 1270.06 1305.43 1330.99 1349.62 1363.19 1372.96 1379.79 1384.29 1386.96 1388.22 1388.53 1388.37 892.07 1134.17 1258.96 1336.09 1387.34 1422.76 1447.80 1465.71 1478.55 1487.61 -_1.- 1493,77.. 1497,60 _ L499.53 1499.92 1499.17 1497.84 983.39 1249.85 1382,59 1461.46 1511.99 1545.82 1569.13 1585.45 1596.91 1604.82 1609.98 1612.88 1613.80 1612.94 1610.57 1607.31 1Q94.5_5 _1385_...._5__2_.._.2 15, 45 96_ 16_74_92 1694,97_.. 1_70_6_-_ _9 1718.12 1724.46 1728.39 1730.23 1730.05 1727.73 1723.15 1716.79 1245.46 1554.54..83.46 _16781__4b _,_____ a __-._O_9._9____._4_,_____Q3__.___g_...___.__.. L...1A4_m 1_9,108....... 1848.68 1845.13 1837.94 1826.26 1494.42 1778.63 1872.76 1914.73 1936.93 1950.02.9__ 5S...1__ ____0163.74 1967.36 1969.68 1970.95 1971.19 1970.09 1966.63 1957.91 1935.73 2046.59 2084.80 2089.56 2091.50 2092.53 2093.13 2093.51 2093.76 2093.92 2094.03. 2094.09 2094, 10 2094.Q5 _______93_ _._8_6_ __0__2_._4_L_.... _?..__ ___.__________ ___.. _ *** ALL. INPUT DATA HAVE.... EEN_ PROCESSED. _ _ __ ___ AT LOC 24C32 E68

Example Problem No. 95 Discussion of Results The purpose of making the temperature distribution calculation is to determine the corner temperature at the tuyere mouth. This temperature is particularly critical in copper converter operation since a low temperature at this point permits the precipitation of magnetite from the slag which blocks the tuyere and requires that the operation be stopped and the tuyere punched out. Tuyere punching is presently one of the limiting steps in copper production rates and adds to the expense and difficulty of the operation. The results of the calculation show that a more highly insulating refractory at the tuyere wall, which in effect decreases the rate of heat transfer to the incoming air, would increase the theoretical temperature at the tuyere mouth and prevent the precipitation of magnetite from the slagthus eliminating the necessity of punching the tuyeres. The accuracy of the result which is obtained depends upon the number of grid points selected. The solution of Krivsky and Schuhmann was carried out by hand calculation and involved first selecting a five inch grid spacing (4 X 4 matrix), utilizing this calculated temperature distribution to estimate a temperature distribution for making the temperature calculation at the tuyere mouth based on a finer grid spacing in the lower right hand corner. The calculated result is markedly influenced by the grid spacing as shown by the data in Table I. The calculations carried out for presentation of the example problem were not extended to extremely fine grid sizes because of the time requirement on the computer. One of the major factors in the time requirement, in addition to the number of matrix points, is the accuracy of the initial temperature distribution. In the present problem a linear distribution was assumed between the melt and the outside or ambient air temperature which would be approximate in the case of a perfectly insulated tuyere, but is far from precise in the case considered in the example problem. TABLE I Computed Corner Temperatures for Various Grid Spacings (Allowable Temperature Change at Grid Point, 0.5~F) Tuyere Heat Corner Transfer Coefficient Grid Spacing Temperature BTU/hr-ftf-OF Matrix Size Inches Iterations OF 10 4x4 5 12 2075 10 6x6 3 23 2067 10 llxll 1.5 58 2054 10 16x16 1 93 2047 200 4x4 5 8 1829 200 6x6 3 17 1753 200 llxll 1.5 40 1659 200 16x16 1 58 1613 The program compiled in 0.;46 minutes and 7.01 minutes of execution time on the IBM 709 were required to generate Table I. E69

Heat Flow and Temperature Distribution Around a Copper Converter Tuyere Critique The problem has not been used in the classroom. It appears, however, to be an excellent application of the computer in the teaching of extractive metallurgy since it is a problem which is well adapted to the computer. As discussed by the authors of the reference, the calculation is an extremely laborious and time consuming process when carried out by hand. The presentation of the problem and its importance in copper processing as well as the details of the solution are more than adequately discussed in the reference paper which would serve as an excellent guide to the student in preparing a computer solution to the problem. E70

Example Problem No. 96 DIGITAL COMPUTER ANALYSIS OF GALVANIC CELL DATA by Robert D. Pehlke and Kenneth J. Guion Department of Chemical and Metallurgical Engineering The University of Michigan Course: Thermodynamics Credit Hours: 4 Level: Junior This course presents the basic thermodynamic principles as applied to engineering processes and includes a laboratory which stresses the experimental techniques for measuring thermodynamic parameters. Statement of Problem The problem selected for computer analysis is the processing of data taken in the laboratory on high temperature galvanic cells. The purpose of the experiment is to measure the thermodynamic properties of metal alloys. The experiment consists of equilibrating electrodes of binary metallic systems immersed in a molten salt bridge and determining the electromotive force between the electrodes and a pure metal standard electrode. These cell potentials are measured at several temperatures and for several compositions across the binary alloy system. The data is utilized to determine the activity and activity coefficient of one alloying element at various temperatures and composition levels. The problem is amenable to measurements made in chemical systems, and could be modified for use with activity measurements determined by other experimental techniques. This program could also be modified to include calculations of the slope of the EMF-temperature curve by a least-squares fit through the data. Utilizing these slopes, the heat of solution and entropy of solution could be computed. Write and test a MAD program which will permit the calculation of activities and activity coefficients from a series of electromotive force measurements on several binary alloys at discrete temperatures. Experiment: The data as measured in the laboratory are presented in the form of weights of the two elements in each electrode and a series of cell potentials measured at several temperatures. Test Data: Weighings gms Cadmium gms Tin Cell 1 1.1785 11.4400 Cell 2 2.2578 10.8039 Cell 3 2.4115 9.9616 Cell 4 2.6463 9.8081 Cell 5) 4.1233 7.7739 Cell 6 5.1567 4.6476 Cell 7I 9.3444 1.3422 E71

Galvanic Cell Data EMF Recordings Tempera- Cell Cell Cell Cell Cell Cell Cell ture, oC (1) (2) (3) (4) (5) (6) (7) 445 57.33 39.76 36.76 34.62 23.82 15.21 5.13 446 57.29 39.83 36.72 34.56 23.73 14.97 5.18 444 57.35 39.76 36.76 34.56 23.87 15.07 5.16 492 62.18 44.24 40.93 38.58 26.14 16.83 5.67 492 62.20 44.26 41.03 38.64 26.09 16.87 5.71 492 62.17 44.18 40.94 38.58 26.11 16.85 5.71 525 64.60 46.59 43.16 40.70 27.28 17.88 6.04 525 64.59 46.63 43.16 40.75 27.26 17.91 6.04 525 64.59 46.71 43.34 40.84 27.26 17.91 6.08 553 70.19 49.14 45.58 42.98 29.00 19.01 6.41 553 69.79 49.12 45.57 42.99 29.33 19.03 6.37 553 69.91 49.29 45.72 43.13 29.05 19.08 6.53 At. Wt. Cadmium 112.41 Tin 118.70 Solution In order to assure equilibrium conditions during the experiment, the cell potentials are measured at periodic intervals and at as constant a temperature as possible. The data are then in the form of several temperatures and the EMF associated with each temperature for a given series of alloy compositions. The program will consist of the following steps: 1. Utilizing the weights of alloying elements and the molecular weights of each element, the weight of each element, the weight percentage and the mole fraction of transferred component will be computed for each cell. 2. The average temperature and EMF associated with each set of cell measurements will be computed. 3. The activity for each cell at each temperature will be calculated. 4. The activity coefficient for each cell at each temperature will be computed. The program as outlined above consists of a series of loops which involve calculations of averages and evaluation of thermodynamic formuli. The relationships which are pertinent to the solution of this problem are summarized as: 1. The thermodynamic activity of an alloying element in solution, given by the equation: -T FE in a = RT 2. The activity coefficient of an element in solution, given by the equation: ai Ni E72

Example Problem No. 96 Flow Diagram -STAR T _ READ name & _ PRINT _ READ mole wts,: B(2)= element title m, cellstemps, readings THROUGH Al, ac { READ weights _ READ temps, FORN=1,1, S(N) (N emf's MFN") H /THROUGH A2,\ /THROUGH A2 TOTAL _ MFP-(N)AI — FOR AA=1,1, - FOR P=1,1, = ^^/:= V/ \AA.G.TEMPS/ \P.G.CELLS+1 THROUGH B, \ INT I PRINT ( t y- PFOR N=1,1, - TOTAL = - B1 VG(AP) -heading- N,G,Z / f = _ ^ A s THROUGH B2,9 PRINT 1 PRINT >( <- ~( FOR N=1,1, ) avgs CONST = heading ~-L \N, G. TEMPS / THROUGH A4, THROUGH A3,\PRINT (> V-~FOR N=1,1, FOR P —=2,1, ACT(NP) temps \N.G.TEMPS / \PG. (CELLS & act \ _____ / \+ 1 ) /' ' --- — PRINT THROUGH A9\ /THROUGH A8,\ ( heading — FOR N=1,1, 9 FOR P —2,1, -GAM(NP) \N.G.TEMPS / \G(CELS / PRINT ^) — temps ~A9 )-PRANSFER - END & gan's 0-/ PO START E73

Galvanic Cell Data MAD Program DIMENSION (' {1 OC,A), X 9 C, b),MF C ) W'[S 5 ),AVG(40C, ), 1 ACT(400,b), GAM(400,8) VECTOR VALUES A= 2,1,2 VECTOR VALUES B= 2,1,4 INIEGER TEtMPS,CELLS, P, N, P, AA,Z START READ FORMAT NAME,NA1,NA2,NA3,NANA, NBI PRINI FORMAT TIILE, NA1, NA2, NA3, NA4, NA5 READ FORMAT UAfA,MOLEB,MCLEC,M,CELLS, rEMPS, Z B(2)= CELLS+1 -.... READ FORMA1 'AtIIGhI, -rT(1,1)...IWTCELLS,2) READ FORMAl1 NUMBElR, X(1,1)...X(Z*TEMPS CELLS*1) _________ ikROUGH A1FCGR N=1,1,N.G.CELLS W fS(N)= 1./(1.+ iT(1N,2) /WT(Ni)) Al- MF(N)=l./(l.+('T((N,2)/MOLEC)*~CLEB/iTI(N,l)) IhROUGH A2, FGR AA=1,1, AA.G.TEMPS ThROUGH A2, FCRP=-1,1, P.G.CELLS+1 TCIAL =.C 1THROUGH B1, FCR N=I,1, N.G.I 81 ICI'AL = X(N+Z*(AA-1),P ) + TOTAL A2 AVG(AA,P) = IOTAL/Z PRIN1 FORMAT SUHL. PRINT FORMAT HEAD1,Nbl, MF 1)...t F(CELLS) PRINT FORMAT HELAU2, hTS(1 )... TSCELLS) PRINI FORMA IJEMPER -t ROUGH b2, FCR N=1,1,N.G.TEMPS 2 -- PRINT FORMAT THREE, AVG(N,)....AVG(N,CELLS+l) LCNST=-23. 66*M/1.981 PRINT FORMAT SUlB PRINT FORMAT HEAD1,NB1, MF(1 )... F (CELLS ) PRINf FORMAl HEAU2, TS(1)... WTS(CELLS) PRINI FORMAT [EMPER THROUGH A4, FCR N=,,N. G.'TEMPS THROUGH A3,FCR P=2,1,P.G.(CELLS+1)A3 ACT(N,P)= ELXP. (CONST*AV(N,P)/(AVG( N,)+2 13.16)) A4 PRINT FORMAT IWO, AVG(N,1), ACT(N,2)...AC'I(N,P-1) PRINI FORMAT SUB4 PRINT FORMAT HEAU1,NbI, MF(1)... F F(CELLS) PRINT FORMAT HEAU2, i^TS(1)...WiS(CELLS) PRINT FORMAT IEMPER lFROUGH A9,FLR N-I,1,N.G.TEMPS Tr kROUGH A8,F L R P=2,1,P.G. (CELLS+1) Ad GAM(N,P) = ACI(N,P)/ PF(P-1) A9 PRINT FORMAI TWO, AVu(i\,1), GAM (,2)... GAM( N,P-1) TRANSFER TC SIARI VECTOR VALUES NAME=$ 5C6, S8, C2 2 $ VEGCTuR VALUES TIILE=', H1, S40, 32H ELLCTRCLYTIC CELL MEASURE 1MENTS // S,6 H OF // S47,.- _5 /I///C.._ _ _____ ____ VECrfR VALUES Ul)AIA 2F1C.5, 4110 *$ VECTUR VALUES WEIGHT=$ (6F10.5) *$ VECTOR VALuES NUMbER=$(7F10.5) *, ECifOR VALUES SUtsb=$ 1Hu,S45,31H AVERA;E CF EXPERIMEltlAL DATA -—. ---- I /1 - $ _VECTCR VALUES SUlI=$ IHl, S52, 1H ACTIVITY ///*$ VECTOR VALUES SUt.4=$ 1H1, S45, 2LH AC IVITY CCEFF IENIL ///*$ VECTOR VALUES HIEAG1=$ 10H ELEMENI, C2, / 16H MOLE FRACIIONJ - -"-...... -------— S -—. ----.. - VECTOR VALUES HEAD2=$ 18H WEIGHI FRACTION, SI, 9Fll.6*i... --- —..-."" —V'TOECi ~-VA1-tUES~- R- 15 1 EMPkRATURE C / * VECTOR VALUES ANS=$ / S19, 9F11.4 *$ VECTOR VALUES TW=S$ /S3, F9.2, Si, 'iE11.5 *$ VECTOR VALUES rHRcE=_= / S3, F9.2, S7, 9F11.3 *$.... O- - ' PROGRAM-' - E74

Example Problem No. 96 Data The format used for putting the data into the computer program was as follows: First Data Card Columns 1-30 contain the title of the experiment. Columns 39-40 contain the symbol of the element under investigation. Second Data Card Columns 1-10 contain the atomic weight of the element under investigation. Columns 11-20 contain the atomic weight of the alloying element. Columns 21-30 contain the number of electrons involved in the transfer. Columns 31-40 contain the number of different electrodes investigated. Columns 41-50 contain the number of different temperatures investigated. Columns 51-60 contain the number of readings taken per temperature setting. Weight Information Cards Columns 1-60 are used on each card and are divided into six groups of ten. Three complete electrode weights can be placed on one card. The weight of the major element and then the alloying element are listed for each electrode. Temperature and EMF Cards Columns 1-70 are used on each card and are divided into seven groups of ten. The temperature and then the EMF readings are listed across the card. The list continues until the data is complete. Below is a sample printout of the Cadmium-Tin data used as the example problem. $DATA -.-_.. —.... C.CADMIUM_ t CD-S.ALLYS_. — C-. —..C —. —. —. g 112.41 118.7 2 7 4 3 _ ..178 5 11. *44 2i2 2578- --- -0 9 —...115 9. --- --- 2.6463 9.8081 4.1233 7.7739 5.1567 4.6476.-....__.......3__.... 9.3 44-.- 3...-. 22 --- ——....-..........4... _..._.............______ __..._ __.__.._.___ _... _ 445. 57.33 39.76 36.76 34.62 23.82 15.21. —.-.....-... - 5-........3. --- -- - 45- - -. --- —-- 572-9- 39.4-839a --- —-- -36-72...- 34-. -6~-5 --- —-.235-73 ---- 14,97 5.18 445. 57.35 39.76 36.76 34.56..-.......-............ 2-38- 1 ---.-.-..5..7 - 5, 6...- -.- 1. 492. +2.1. 4... --- 4 24 - 40. 9-3 38,58 26.14 16.83 5.67 492. 62*20 44*26 _ 4.1-4_ 4. 3 --- - 26.-Q9- -.....16.8-. -----— 7 ----.- -- — 4.2. - +-17 — 44.18 40.94 38.58 26.11 16.85 5.71 525..-_._...6....__- ___6 ------- 4.596_ __ ------— 40-7 ----- -------- 7. -. ----6 —+04 --- 525. 64,59 46*63 43.22 40.75 27.26 17.91.. _ -6..,O,Z5...=.6.: 525 4. ---.6. 4..5 9_-4 43L_.474_" N- *4O.: c& - ------ 27-267 17.91 6.08 553, 70,19 49,14 45.58 42.98.....290-.. -l...l._. 4-5___ _53.. 9__1..____.6979 —_I..45-157. 42.99 42.33 19.03 6.37 553. 69.91 49.29...-. ---. —..-45.- -712 --- —43 3- -S9*3.9.. -653- - ----.. ---- -----—. ---_ E75

Galvanic Cell Data Computer Output ELECTROLYTIC CELL MEASUREMENTS OF CADMIUM IN CU-SN ALLOYS AVERAGE OF EXPERIMENTAL DAiA ELEMENT CD MOLE FRACTION.098108.186780.203584.221732.359008.539515.880262 WEIGHT FRACTION.0933%5.172857.194899.212479.346577.525963.874403 IEMPERATURE C 445.00 57.323 39.783 36.747 34.580 23.807 15.083 5.157 492.00 62.183 44.227 40.967 38.600 26.113 16.850 5.697 525.00 64.593 46.643 43.363 40.763 27.267 17.900 6.053 553.00 69.963 49.183 45.623 43.033 33.487 19.040 6.347 ACTIVI Y ELEMENT CD MOLE FRACTION.0981C8.180780.203584.221732.359008.539515.880262 -WE-TTr-F~A-T-ON.3 --- —------.-593 --- —-..T7S7.T89 -. 1247..... -. 34 65-7-..- 5275 3.* 87440-^..3.. TEMPERATURE C 445.00.15674E 00.27634E 00.30484E 00.32696E 00.46318E CO.61409E 00.84645E 00 492.00.15156E C0.26134E 00.28851E 00.30999E 00.45278E 00.59973E 00.84126E 00 525.CO.15276E 00.25749E 00.28327E 00.30552E 00.45242E 00.59412E 00.83855E 00 553.00.14000E CC.25104E 00.27745E 00.29840E 00.39022E 00.58563E 00.83665E 00 ACTIVITY COEFFICIENI E-L-EM E N — CD --- ----- --------------------------------------- ACT --- ELEMENIr CU MOLE FRACTION.098108.180780.203584.221732.359008.539515.880262 ~-wrr-FwTr-tKt-CTTO-N '. --- —-— 9S3-5 --- -.T7 — 87. 194899.9.-2 12 479- ".-34-65-77- — 525)963..814403 FEMPERATURE C 445.00.15976E 01.15286E 01.14974E 01.14746E 01.12902E 01.11382E 01.96159E GO 492.00.15448E 01.14456L C1.14171E 01.1398QE 01.12612C 01.11116E 01.955708 00 525.00.15571E 01.14243E 01.13914E 01.13779E 01.12602E 01.11012E 01.95262E 0C 553.CO.14270E 01.13886E 01.13628E 01.13458E 01.10869E 01.10855E 01.95045E 0O,**** LL INPUT DATA HAVE BEEN PRLCESSED._... _ __.__. E76

Example Problem No. 96 Discussion of Results The computer program prints out three tables which list the averages of the data, the activities, and the activity coefficients for the various values of temperature and composition investigated. The MAD program compiled in 0.52 minutes and performed all the required operations for the Cadmium-Tin example in 0.34 minutes of execution time. Concerning the results, it is important to realize that the number of significant figures presented by the computer are not necessarily correct, but instead depend on the number of significant figures in the original data. Critique This particular problem is especially well suited to a laboratory course in which a series of data are taken and require considerable time for analysis. This particular problem has not been used in instruction of the course in thermodynamics, but would be especially well suited to the laboratory supplement to the course. One important use of the computer in laboratories and engineering sections of industrial organizations is for routine data analysis. In many cases, small computers are especially wellsuited for this application. To emphasize this aspect of computer use would be of particular benefit to the engineering student. The program for analyzing electromotive force measurements has been used in research work at The University of Michigan, and although it has not been prepared for the small computer, the program would be especially well suited to it. The program contains several computing steps involving data averaging, and formula evaluation. The problem, thus, appears to be well suited for the purpose of illustrating data analysis using the digital computer. E77

VII. REFERENCES 1. First Annual Report, Project on Use of Computers in Engineering Education, The University of Michigan, College of Engineering, Ann Arbor, Michigan, August, 1960. 2. Second Annual Report, Project on Use of Computers in Engineering Education, The University of Michigan, College of Engineering, Ann Arbor, Michigan, December, 1961. 3. Sinnott, M. J. and R. D. Pehlke, "Computers in Undergraduate Teaching of Metallurgical Engineering," Journal of Engineering Education, Vol. 52, No. 9, May, 1962. 4. Arden, B. W., An Introduction to Digital Computing, Addison-Wesley Publishing Co., Reading, Massachusetts, 1962. 5. Organick, E. I., A Computer Primer for the MAD Language, Cushing-Malloy, Inc., Ann Arbor, Michigan, 1961. E78