A Computational Package for Design and Process Evaluation in Bandsawing - A User's Manual by James E. Borchelt A. Galip Ulsoy and Panos Papalambros UM-MEAM-83-14 Department of Mechanical Engineering and Applied Mechanics and Integrated Design and Manufacturing Unit of the Center for Robotics and Integrated Manufacturinq The University of Michigan Ann Arbor, MI 48109

, v I t1 ^

Table of Contents Page Description of Program 2 Sample Terminal Session 11 Output from Sample Terminal Session 19 Program Listing 48 1

Description of Program The stress analysis program developed was combined with the previously developed vibration analysis program ("Two Computer Codes for Band Saw Vibration and Stability Analysis", A.G. Ulsoy, University of California, Berkley, 1978.) to form a computationally efficient interactive design program. The report "Development of an Efficient Computational Procedure for Evaluating Band Saw Blade Stresses" (J.E. Borchelt, A.G. Ulsoy, P. Papalambros, University of Michigan, 1983) presents the theory used to develop the stress analysis program. The report "Vibration and Stability of Band Saw Blades: A Theoretical and Experimental Study" (Ph.D. Dissertation, A.G. Ulsoy, University of California, Berkley, 1979.) presents the theory used to develop the vibration analysis program. Figure 1 shows the total desired automatic design computation package. The process parameters, are P, a are the in-plane stresses, w are the natural frequencies, and J is the performance criterion. The existing program consists of the stress and vibration analysis segments and can stand alone as a design package. The user can iteratively use the existing program as a computer-aided design package. The stresses and/or natural frequencies can be computed for various design parameters and process conditions, allowing the designer to make direct comparisons. Three possible paths exist that can be taken by the program. Figure 2 shows which path(s) the user can choose. The path taken by the program can be changed once a given problem has been solved without restarting the program. Figure 3 shows the subroutines called when only the in-plane stress analysis, Path 1, is chosen 2

3 L -H OPTIMIZATIONI PROGRAM R^*, STRESS ANALYS IS PROGRAM I - -L. -- PERFORMANCE CRITERION.R. CRITERION. J FIGURE 1 SCHEMATIC OF DESIGN APPROACH.

4 PATH 1: ANALYSIS OF IN-PLANE STRESSES ONLY PATH 2: VIBRATION ANALYSIS ONLY PATH 3: STRESS AND VIBRATION ANALYSIS TOGETHER NOTE: THE ORDER OF THE APPROXIMATING FUNCTIONS FOR THE DISPLACEMENT FIELDS MUST BE 4 WHEN PATH 3 IS CHOSEN. FIGURE 2 POSSIBLE PATHS THROUGH PROGRAM

1AIN PROG RAM1 1E RHS STROUT TEMP EXPAND I -— I.CM O D DLUD 2 MATDIM DBS CMOD FIGURE 3 SUBROUTINES CALLED FOR PLANE STRESS ANALYSIS

MAIN PROGRAM VIBRA an FIGURE 4 SUBROUTINES CALLED FOR VIBRATION ANALYSIS

7 by the user. The subroutines called when the vibration analysis, Path 2, alone is chosen by the user are shown in Figure 4. When the combined analysis, Path 3, is chosen by the user subroutines from both Figures 3 and 4 are used. A brief description of what each subroutine does follows: SUBROUTINE TDRIVE: SUBROUTINE SDRIVE: SUBROUTINE TCOFF: SUBROUTINE ATAM: SUBROUTINE TOUT: SUBROUTINE FORCE: SUBROUTINE TEMP: SUBROUTINE RHS: SUBROUTINE EXPAND: SUBROUTINE DLUD: SUBROUTINE DBS: SUBROUTINE MATDIM: SUBROUTINE CMOD: Reads temperature data for various points on the band saw blade. Calls subroutine TCOFF. Reads (Pre) stress data for various points on the band saw blade. Calls subroutine TCOFF. Calculates temperature Field coefficients using the data read in subroutine TDRIVE by solving the linear regression problem. Temperature field coefficients are used in subroutine TEMP. Calls subroutines ATAM, DLUD, and DBS. Calculates the matrix A-transposed times A for use in finding the temperature field coefficients. Produces a table of the blade temperatures at various points using temperature field coefficients calculated in subroutine TCOFF. Calculates the effect of the tensile, normal, and/or tangential loads for the righthand side vector. Also determines number of terms in approximating function for displacement field, (U). Calculates effect of a temperature gradient if present for the righthand side vector. Adds the load and temperature effects together to form the righthand side vector. Calculates submatrices as given in report. Determines the displacement field coefficients for U and V. Prints out the displacement field coefficients A and B, if desired. Calls subroutines MATDIM, DLUD, DBS, DOUT, and SCOFF. Decomposes a general matrix (A) into lower and upper triangular matrices. Solves a system of linear equations by forward elimination and back substitution. Calculates nondimensional integrals for use in subroutine expand to calculate submatrices K1,K2,K3, and K4. Calls subroutine CMOD. Calulates the powers of x and y used to determine the elements of K1,K2,K3, and K4. Bridge subroutine that calls AMOD 1, AMOD2, AMOD3, AMOD4, AMOD5, AMOD6, AMOD7.

8 SUBROUTINE AMOD 1: SUBROUTINE AMOD 2: SUBROUTINE AMOD 3: SUBROUTINE AMOD 4: SUBROUTINE AMOD 5: SUBROUTINE AMOD 6: SUBROUTINE SUBROUTINE AMOD 7: DOUT: Calculates the nondimensional integrals that are multiplied by (EH)/(1-NU**2)*(B/L) in Kl. Calculates the nondimensional integrals that are multiplied by (EH)/(1+NU)*(L/B) in Kl. Calculates the nondimensional integrals that are multiplied by (EH)/(1-NU) in K2. Calculates the nondimensional integrals that are multiplied by ((EH)NU)/(1-NU**2) in K2. Calculates the nondimensional integrals that are multiplied by (EH)/(1+NU) in K2. Calculates the nondimensional integrals that are multiplied by (EH/(1-NU**2)*(L/B) in K4. Calculates the nondimensional integrals, that are multiplied by (EH)/(l+NU)*(B/L) in K4. Produces displacement field tables for U and V at various points on the blade, using the displacement field coefficients determined in subroutine expand. Calculates the stress field coefficients for Sigma-X, Sigma-Y and Tau-XY, using the displacement field coefficients and Hooke's Law. Produces tables of stress at various points on the blade. Tables made are Sigma-X, Sigma-Y, shear stress, and the maximum principal normal stress. Reads dimensionless integrals for vibration analysis from file. Calculates mass, gyro, and stiffness matrices. Picks out the right dimensionless integral to be used in calculating the mass, gyro, and stiffness matrices. Changes the system of second order differential equations to a system"of first order differential equations. Calculates the natural frequencies of the band saw blade, and prints out the 5 lowest eigenvalues. SUBROUTINE SCOFF: SUBROUTINE STROUT: SUBROUTINE VIBRA: SUBROUTINE DECODE: SUBROUTINE REFORM: SUBROUTINE EIGENS: The temperature field and/or (pre) stress field coefficients are found by linear regression. The data points, blade coordinates x and y, and the temperature or (pre) stress values, are entered by the user or read from a data file. ATAM then calculates a NxN norm matrix, where N is the number of terms in the approximating function for the given temperature, or (pre) stress field. The norm matrix is formed by matrix multiplication of the transpose of a MxN matrix times the MxN matrix. M is the number of data points.

9 This method is inherently ill-conditioned so the method of singular-value decomposition and pseudo-inverse, (Numerical Methods, G. Dahlquist, and A. Bjorck, Prentice-Hall, Englewood Cliffs, N.J., 1974. Pp. 143-144), should be used to reduce the illconditioning problem. Important remarks conerning the use of this program are: 1) When a combined stress and vibration analysis are chosen the order of the displacement field approximating functions must be 4. If a combined analysis is chosen initially the program automatically set the order of the displacement field approximating functions equal to 4. However, this is not the case if the combined analysis is not chosen initially. An error message is printed out, and the run is terminated if the above condition is not met. 2) The highest order of the displacement field approximating functions is 7. The program automatically sets the order to 6 if it is greater than 7, and sets the order to 2 if it is less than 2. 3) A subroutine is needed to combine the stress and prestress field coefficients. The subroutine for calculation of prestress coefficients is SDRIVE subroutine SDRIVE is not working properly. 4) Input/Output file numbers: 5 - Terminal input 6 - Terminal output 7 - Input file that contains the nondimensional integrals used for the vibration analysis 8 - Output file 10 - Input file that contains the blade coordinates, and stress values that are used with vibration analysis alone, or prestress values. 11 - Input file that contains the blade coordinates, and temperature data for use in calcualting the tempera

ture field coefficients. 5) Maximum number of data points that can be entered in subroutines TDRIVE, or SDRIVE is 70. Further information on subroutines DLUD, and DBS can be obtained from the University of Michigan Computing Center CCMEMO 426. Information concerning subroutines BALANC, ELMHES, and HQR can be obtained from the University of Michigan Computing Center by using the Archival utility program (NAASA: AUSP) to list documentation on the above subroutines. Also see the program listings of the above suroutines in a subsequent section for further references. The next section gives a sample of a terminal session using the present program. A listing of the program is provided in the last section.

SAMPLE TERMINAL SESSION 11

12 %LF10: LA36 MTS Ann Arbor (LA36,LF10-GLAB-MP02-CD23,01142) * M.TS. will be down from 11 30 Pm Sat till 12 noon Sun (6-12-83) *SIGNON K921 *Enter Password. *Terminal NormalUrniv/Gov't 4Last sirnon was at 16t11:51, Fri Jun 10/83 *User K921 signed on at 16:48:57, Fri Jun 10/83 *Messases for mailbox K921: 2 new* Use $MESSAGE to retrieve then. *CON *SINK* UC=ON #RUN *FTN SCARDS=BANDSAW SPUNCH=-B Exeecution begins No errors in MAIN *Execution term inated iRER SCARDS=STRESS.CODE SPUNCH=-S *$Run *FTN SCARDS=STRESS.CODE SPUNCH=-S #Execution begins No errors in EXPAND No errors in ZEROC No errors in MATDIM No errors in CMOD No errors in AMOD1 No errors in AMOD2 No errors in AMOD3 No errors in AMOD4 No errors in AMOD5 No errors in AMOD6 No errors in AMOD7 No errors in SCOFF No errors in ZEROS No errors in DOUT No errors in STROUT No errors in TDRIVE No errors in TCOFF No errors in ATAM No errors in TOUT No errors in TEMP No errors in FORCE No errors i n RHS No errors in SDRIVE *Execution terlainrated #RER SCARDS=VIBRATION CO SPUNCH=-V *$Run *FTN SCARDS=VIBRATION CO SPUNCH=-V EExecutior besins No errors in VIBRA No errors in DECODE No errors in REFORM No errors in EIGENS fExecution terrminrted $

13 RUN -B+-S+-V+NAAS:NAL+NAAS:EISPACK 7=TAPE 8=-OUT 11=IN1 lExecution begins BEGIN PROBLEM DESCRIPTION ENTER THE TITLE OF THIS RUN (50 CHARACTERS MAX) TEMPERATURE GRADIENT WITH A TENSILE LOAD ENTER NUMBER OF PROBLEM TO BE SOLVED 1-PLATE STRESSES 2-FLATE VIBRATION 3FPLATE STRESSES AND VIBRATIONS 1 BEGIN DESCRIPTION OF APPROXIMATING FUNCTIONS FOR DISPLACEMENT FIELD COEFFICIENT CALCULATION ENTER DEGREE OF THE POLYNOMIAL FOR THE APPROXIMATING FUNCTIONS INTEGER 2 TO 7 4 BEGIN DESCRIPTION OF APPROXIMATING FUNCTION FOR TEMPERATURE FIELD COFFICIENT CALCULATION ENTER DEGREE OF THE POLYNOMIAL FOR THE APPROXIMATING FUNCTION IF ZERO NO TEMPERATURE FIELD IS PRESENT OTHERWISE ENTER INTEGER 1 TO 3 3 BEGIN ENTRY OF TEMPERATURE DATA IS THE DATA TO READ FROM A DATA FILE 0-NO 1-YES 1 ENTER THE NUMBER OF TEMPERATURE DATA POINTS (12 FORMAT) 66 BEGIN DESCR-IPTION OF GEOMETRY OF BANDSAW BLADE ENTER LENGTH OF BANDSAW BLADE 39.37 ENTER BANDSAW BLADE WIDTH 10.827 ENTER BANDSAW BLADE THICKNESS 0.059055 BEGIN DESCRIPTION OF MATERIAL PROPERTIES OF THE BANDSAW BLADE ENTER MODULUS OF ELASTICTY (YOUNG'S MODULUS) 29010000.0 ENTER FOISSON"S RATIO 0.3

14 ENTER THE MASS DEENSITY JF iF HMAITERAL 0.283 ENTER COEFFICIENT OF THERMAL EXPANSION 0,0000084 BEGIN DESCRIPTION OF OPERATING CONDITIONS OF THE BANDSAW BLADE ENTER THE VELOCITY OF THE BANDSAW BLADE 0.0 ENTER WHEEL SUPPORT COEFFICIENT 0.0 ENTER THE INITIAL TENSION ON THE BANDSAW BLADE 11240.0 ENTER NORMAL CUTTING FORCE ON THE BANDSAW BLADE 0.0 ENTER TANGENTIAL CUTTING FORCE ON THE BANDSAW BLADE 0.0 BEGIN OUTPUT CONTROL DESCRIFTION DO YOU WISH THE PRINTING OF THE DISPLACEMENT FIELD COEFFICIENTS (Y/N)? Y DO YOU WISH THE PRINTING OF THE DISPLACEMENT FIELD TABLES (Y/N)? Y DO YOU WISH THE PRINTING OF THE STRESS FIELD TABLES (Y/N)? Y DO YOU WISH THE PRINTING OF THE TEMPERATURE FIELD COEFFICIENTS/TABLE (Y/N)? Y ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APPROXIMATING FUNCTION FOR TEMPERATURE 4- NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OPERATING CONDITIONS 8- SOLVE NEW PROBLEM 9- STOP 1 ENTER NUMBER OF PROBLEM TO BE SOLVED 1-PLATE STRESSES 2-PLATE VIBRATION 3-PLATE STRESSES AND VIBRATIONS BEGIN DESCRIPTION OF EIGENFUNCTIONS FOR BANDSAW BLADE VIBRATIONS ENTER NUMBER OF SIMPLY SUPPORTED BEAM EIGENFUNCTIONS INTEGER 1 TO 6 4 ENTER NUMBER OF FREE-FR EE BEAM EIGENFUNCTIONS INTEGER 1 TO 6 2

15 ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APPROXIMATING FUNCTION FOR TEMPERATURE 4- NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OPERATING CONDITIONS 8- SOLVE NEW PROBLEM 9- STOP 8 ENTER THE TITLE OF THIS RUN (50 CHARACTERS MAX) PLATE VIBRATION WITH TEMPERATURE GRADIENT DO YOU WISH TO USE THE STRESS FIELD COEFFICIENTS CALCULATED PREVIOUSLY (Y/N)? Y ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APPROXIMATING FUNCTION FOR TEMPERATURE 4- NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OPERATING CONDITIONS 8- SOLVE NEW PROBLEM 9- STOP 1 ENTER NUMBER OF PROBLEM TO BE SOLVED I-PLATE STRESSES 2-PLATE VIBRATION 3-FLATE STRESSES AND VIBRATIONS BEGIN DESCRIPTION OF APPROXIMATING FUNCTIONS FOR DISPLACEMENT FIELD COEFFICIENT CALCULATION ENTER DEGREE OF THE POLYNOMIAL FOR THE APPROXIMATING FUNCTIONS INTEGER 2 TO 7 6 ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APPROXIMATING FUNCTION FOR TEMPERATURE 4- NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OPERATING CONDITIONS 8- SOLVE NEW PROBLEM 9- STOP 3 BEGIN DESCRIPTION OF AP'PROXIMATING FUNCTION FOR TEMPERATUJII: FEIrD rlnFTICIENT CALCULATION ENTER DEGREE OF THE POLYNOMIAL FOR THE APPROXIMATING FUNCTION IF ZERO NO TEMPERATURE FIELD IS PRESENT OTHERWISE ENTER INTEGER 1 TO 5 0

16 ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APPROXIMATING FUNCTION FOR TEMPERATURE 4- NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OPERATING CONDITIONS 8- SOLVE NEW PROBLEM 9- STOP 5 BEGIN DESCRIPTION OF GEOMETRY OF BANDSAW BLADE ENTER LENGTH OF BANDSAW BLADE 1.0 ENTER BANDSAW BLADE WIDTH 0.25 ENTER BANDSAW BLADE THICKNESS 0.00107 ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APF'FROXIMATING FUNCTION FOR TEMPERATURE 4 — NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OFERATING CONDITIONS 8- SOLVE NEW P'ROBLEM 9- STOP 6 BEGIN DESCRIPTION OF MATERIAL PROPERTIES OF THE BANDSAW BLADE ENTER MODULUS OF ELASTICTY (YOUNG'S MODULUS) 20700000000.0 ENTER POISSON'S RATIO 0.:3 ENTER THE MASS DENSITY OF THE MATERIAL 0.283 ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APPROXIMATING FUNCTION FOR TEMPERATURE 4 — NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OPERATING CONDITIONS 8- SOLVE NEW FROBLEM 9- STOP

17 7 BEGIN DESCRIPTION OF OPERATING CONDITIONS OF THE BANDSAW BLADE ENTER THE VELOCITY OF THE.BANDSAW BLADE 48.74 ENTER WHEEL SUPPORT COEFFICIENT 0.2 ENTER THE INITIAL TENSION ON THE BANDSAW BLADE 50000.0 ENTER NORMAL CUTTING FORCE ON THE BANDSAW BLADE 250.0 ENTER STARTING X-COORDINATE OF NORMAL CUTTING FORCE 0.487 ENTER ENDING X-COORDINATE OF NORMAL CUTTING FORCE 0.513 ENTER TANGENTIAL CUTTING FORCE ON THE BANDSAW BLADE 500.0 ENTER STARTING X-COORI:DINATE OF TANGENTIAL CUTTING FORCE 0.487 FNTFR ENDING X —CrORrITNATE CF TANGENTIAI. CUITTING FORCFE 0.513 ENTER TANGENTIAL CUTTING FORCE ON THE BANDSAW BLADE 500*0 ENTER STARTING X-COORDINATE OF TANGENTIAL CUTTING FORCE 0.487 ENTER ENDING X-COORDINATE OF TANGENTIAL CUTTING FORCE 0.513 ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TYPE 2- NEW APPROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW APPROXIMATING FUNCTION FOR TEMPERATURE 4- NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW OPERATING CONDITIONS 8- SOLVE NEW PROBLEM 9- STOP 8 ENTER THE TITLE OF THIS RUN (50 CHARACTERS MAX) TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 ENTER NUMBER OF NEXT INSTRUCTION 1- NEW PROBLEM TY'E 2- NEW APPF'ROXIMATING FUNCTION FOR DISPLACEMENTS 3- NEW AIP:'PROXIMATING FUNCTION FOR TEM:'ERATURE 4- NEW EIGENFUNCTIONS FOR VIBRATION ANALYSIS 5- NEW BLADE GEOMETRY 6- NEW MATERIAL PROPERTIES 7- NEW (OPERATING CONDITIONS 8- SOLVE NEW PROBLEM 9- STOP ENTER NUMBER OF PROBLEM TO BE SOLVED 1-PLATE STRESSES 2-PLATE VIBRATION 3-PLATE STRESSES AND VIBRATIONS 2 BEGIN DESCRIPTION OF EIGENFUNCTIONS FOR BANDSAW BLADE VIBRATIONS

18 ENTER NUMBER OF SIMPLY SUPPORTED BEAM EIGENFUNCTIONS INTEGER 1 TO 6 5 ENTER NUMBER OF FREE-FREE BEAM EIGENFUNCTIONS INTEGER 1 TO 6 **** ERROR THE STRESS FIELD POLYNOMIALS DO NOT MATCH THOSE REQUIRED FOR THE VIBRATION ANALYSIS! *** tExecution terminated fRUN *PAGEF'R SCARDS=-OU PAR=ONESIDED.Version Jun 02/83. $Copv *PAGENEWS(1,30) for description of changes* EEx)ecut ion begins *F'RINT* assigned receipt rumber 667585 *FRINT* 667585 held This *PAGEPR run benerated 634 linesy 28 Pagesv 28 images arnd 28 sheets* *FPRINT 667585 released, 30 Pages, route=CNTR, Printer=PAGE* lExecution te rrinated tSIG tReceipt SuJnmmaryr *.FPRINT* 667499 61 Pages, route=CNTR, Printer=PAGEo;*F'RINT* 667585 30 pames, route=CNTR, Printer=PAGE0 IK921 16:48:57 to 17:24:39, tTerminralNormalUniv/Gov' t I:Elapsed time 35.684 FCPU time used 5.553 tCFPU storage VMI 7.483 IWait storage VMI 67,475 1F'ame Printer lines 2357 IFPase Printer Pases 91 WPase Printer images 91 WPame printer sheets 91 WPage-ir s 577 -Dlisk. I/O 1639 \1 Approximate cost of th Fri Jun 10/83 minutes seconds Pa i e- m i r. page-hro $.81 $2.36 $*95 uis run,: $6,26 tDisk. storaEe charge 101 pase-hr. t Approximate remainirn balarnce: $.01 $43.43

OUTPUT FROM SAMPLE TERMINAL SESSION 19

TEMPERATURE GRADIENT WITH A TENSILE LOAD THE DEGREE OF THE POLYNOMIAL FOR THE DISPLACEMENT FIELD CALCULATION IS 4 THE DEGREE OF THE POLYNOMIAL FOR THE TEMPERATURE FIELD CALCULATION IS 3 THE LENGTH OF THE BANDSAW BLADE IS 0.393700E+02 THE WIDTH OF THE BANDSAW BLADE IS 0.108270E+02 THE THICKNESS OF THE BANDSAW BLADE IS 0.590550E-01 THE MODULUS OF ELASTICITY IS 0.290100E+08 POISON"S RATIO IS 0.300000E+00 THE MASS DENSITY IS 0.283000E+00 THE COEFFICIENT OF THERMAL EXPANSION IS 0.840000E-05 THE VELOCITY OF THE BANDSAW BLADE IS 0.0 THE WHEEL SUPPORT COEFFICIENT IS 0.0 THE INITIAL TENSION IS 0.112400E+05 THE NORMAL CUTTING FORCE IS 0.0 THE TANGENTIAL CUTTING FORCE IS 0.0 t) 0

TEMPERATURE GRADIENT WITH A TENSILE LOAD TEMPERATURE FIELD COEFFICIENTS T(O,O)= 0.13925E+03 T(1,0)=-0.31567E+01 T(2,0)= 0.19117E+00 T(3,0)=-0.30164E-02 T(1,1)= 0.52621E+00 T(2,1)=-0.14790E-01 T(1,2)= 0.84362E-02 T(0,1)=-0.51549E+01 T(0,2)= 0.17958E+00 T(0,3)=-0.33431E-03 I' —

TEMPERATURE GRADIENT WITH A TENSILE LOAD TEMPERATURE DISTRIBUTION Y I I COORDINATE I I I I 0.108E+02 I 0.104E+03 0.118E+03 0.132E+03 0.145E+03 0.156E+03 0.163E+03 0.165E+03 0.162E+03 0.152E+03 0.134E+03 0.107E+03 I I I 0.866E+01 I 0.108E+03 0.117E+03 0.126E+03 0.136E+03 0.144E+03 0.149E+03 0.151E+03 0.149E+03 0.140E+03 0.125E+03 0.102E+03 I I I 0.650E+01 I 0.13E03 0.117E+03 0.12E+03.129E+03 0.135E+03 0.13E+03 E+03 0. 0E +0 3 03 0.133E+03 0.120E+03 0.10E+03 I I I 0.433E+01 I 0.120E+03 0.119E+03 0.121E+03 0.125E+03 0.129E+03 0.132E+03 0.134E+03 0.134E+03 0.129E+03 0.120E+03 0.105E+03 I I I 0.217E+01 I O.129E+03 0.123E+03 0.122E+03 0.123E+03 0.125E+03 0.129E+03 0.131E+03 0.132E+03 0.130E+03 0.124E+03 0.114E+03 I I I 0.0 I 0.139E+03 0.130E+03 0.125E+03 O.i24E+03 0.125E+03 0.128E+03 0.132E+03 0.134E+03 0.135E+03 0.133E+03 0.127E+03 I I I I 0.0 0.394E+01 0.787E+01 0.118E+02 0.157E+02 0.197E+02 0.236E+02 0.276E+02 0.315E+02 0.354E+02 0.394E+02 I I X-COORDINATE I N)

TEMPERATURE GRADIENT WITH A TENSILE LOAD COEFFICIENTS A( 1, O)= A( 2, 0)= A( 3, 0)= A( 4, 0)= A( 1, 1)= A( 2, 1)= A( 3, 1)= A( 1, 2)= A( 2, 2)= A( 1, 3)= A( O, 1)= A( 0, 2)= A( 0. 3)= A( 0, 4)= OF APPROXIMATING FUNCTION FOR DISPLACEMENT FIELD U 0.54463E-03 -0.16677E-04 0.59267E-06 -0.62723E-08 -0.10565E-04 0.36739E-05 -0.67150E-07 -0.11390E-05 0.42909E-07 -0.41193E-07 -0.66577E-03 0.31867E-05 0.62317E-06 0.76584E-08,) WA

TEMPERATURE GRADIENT WITH A TENSILE LOAD COEFFICIENTS B( 1, O)= B( 2, 0)= B( 3, O)= B( 4, 0)= B( 1, 1)= B( 2, 1)= B( 3, 1)= B( 1, 2)= B( 2, 2)= B( 1, 3)= B( 0, 1)= B( O, 2)= B( 0. 3)= B( 0, 4)= OF APPROXIMATING FUNCTION FOR DISPLACEMENT FIELD V 0.67075E-03 0.56762E-05 -0.12378E-05 0.16788E-07 -0.59324E-05 0.10317E-05 -0.26166E-07 -0.20579E-05 0.66575E-07 -0.30633E-07 -0.31032E-04 -0.71422E-04 0.50149E-05 -0.11227E-06

TEMPERATURE GRADIENT WITH A TENSILE LOAD DISPLACEMENT FIELD U Y I I COORDINATE I I I I 0.108E+02 I -.594E-02 -.455E-02 -.236E-02 0.529E-03 0.398E-02 0.781E-02 0.118E-01 0.157E-01 0.193E-01 0.222E-01 0.240E-01 I I I 0.866E+01 I -.508E-02 -.345E-02 -.128E-02 0.140E-02 0.450E-02 0.789E-02 0.114E-01 0.149E-01 0.180E-01 0.207E-01 0.224E-01 I I I 0.650E+01 I -.401E-02 -.222E-02 -.115E-03 0.231E-02 0.501E-02 0.793E-02 0.110E-01 0.140E-01 0.168E-01 0.191E-01 0.208E-01 I I I 0.433E+01 I -.277E-02 -.886E-03 0.109E-02 0.321E-02 0.550E-02 0.794E-02 O.105E-01 O.130E-01 0.155E-01 O.176E-01 0.193E-01 I I I 0.217E+01 I -.142E-02 0.505E-03 0.231E-02 0.410E-02 0.595E-02 0.791E-02 0.996E-02 0.121E-01 0.142E-01 0.162E-01 0.179E-01 I I I 0.0 I 0.0 0.192E.352E02 0496E02 0637E-002 0.784E-02 0.942E-02 O.111E-1 0.130E-01 O.148E-01 0.167E-01 I I I I 0.0 0.394E+01 0.787E+01 0.118E+02 0.157E+02 0.197E+02 0.236E+02 0.276E+02 0.315E+02 0.354E+02 0.394E+02 I I X-COORDINATE I N) U,

TEMPERATURE GRADIENT WITH A TENSILE LOAD DISPLACEMENT FIELD V Y I I COORDINATE I I I I 0.108E+02 I -.389E-02 -.231E-02 -.466E-03 0.123E-02 0.246E-02 0.301E-02 0.276E-02 0.168E-02 -.159E-03 -.258E-02 -.532E-02 I I I 0.866E+01 I -.300E-02 -.103E-02 0.107E-02 0.291E-02 0.418E-02 0.471E-02 0.438E-02 0.319E-02 0.123E-02 -.130E-02 -.412E-02 I I I 0.650E+01 I -.204E-02 0.227E-03 0.251E-02 0.443E-02 0.572E-02 0.620E-02 0.579E-02 0.451E-02 0.247E-02 -.112E-03 -.294E-02 I I I 0.433E+01 I -.111E-02 0.137E-02 0.376E-02 0.572E-02 0.699E-02 0.741E-02 0.694E-02 0.559E-02 0.352E-02 0.952E-03 -.179E-02 I I I 0.217E+01 I -.354E-03 0.225E-02 0.469E-02 0.664E-02 0.787E-02 0.823E-02 0.770E-02 0.634E-02 0.429E-02.180E-02 -.762E-03 I I I I 0.0 I 0.0 0.266E-02 0.509E-02 0.700E-02 0.817E-02 0.848E-02 0.792E-02 0.657E-02 0.460E-02 0.229E-02 0.139E-16 I I 0.0 0.394E+01 0.787E+01 0.118E+02 0.157E+02 0.197E+02 0.236E+02 0.276E+02 0.315E+02 0.354E+02 0.394E+02 I I X-COORDINATE I

TEMPERATURE GRADIENT WITH A TENSILE LOAD COEFFICIENTS OF NORMAL AX(O,O)= 0.17066E+05 AX(1,O)=-0.11200E+04 AX(2,0)= 0.66548E+02 AX(3,0)=-0.10501E+01 AX(1,1)= 0.19488E+03 AX(2,1)=-0.51487E+01 AX(1,2)= 0.18569E+01 AX(O,1)=-0.17029E+04 AX(0,2)= 0.10757E+03 AX(0.3)=-0.56083E+01 STRESS IN X-DIRECTION to -i

TEMPERATURE GRADIENT WITH A TENSILE LO, COEFFICIENTS OF NORMAL AX(O,O)= 0.42194E+04 AX(1,0)=-0.50810E+03 AX(2,0)= 0.49895E+02 AX(3,0)=-0.10741E+Oi AX(1,1)=-0.60939E+02 AX(2,1)= 0.23181E+01 AX(1,2)=-0.21090E+01 AX(0,1)=-0.46548E+04 AX(0,2)= 0.46872E+03 AX(0,3)=-0.14711E+02 STRESS IN Y-DIRECTION 00

TEMPERATURE GRADIENT WITH A TENSILE LOAD COEFFICIENTS OF SHEAR STRESS AX(O,O)= 0.55516E+02 AX(1,O)= 0.87845E+01 AX(2,0)=-0.44214E+00 AX(3.0)= 0.0 AX(1,1)=-0.23936E+01 AX(2,1)= 0.81674E-01 AX(1,2)= 0.10680E+00 AX(O,1)= 0.49201E+01 AX(0, 2)=-. 21024E+01 AX(0.3)= 0.0

TEMPERATURE GRADIENT WITH A TENSILE LOAD SIGMA X Y I I COORDINATE I I I I 0.108E+02 I 0.412E+04 0.898E+04 0.138E+05 0.182E+05 0.217E+05 0.241E+05 0.248E+05 0.236E+05 0.201E+05 0.138E+05 0.433E+04 I I I 866E I 0.674E+04 0.980E+04 0.132E+05 0.164E+05 0.192E+05 0.212E+05 0.219E+05 0.209E+05 0.180E+05 0.127E+05 0.453E+04 I I I 0.650E+01 I 0.901E+04 0.103E+05 0.123E+05 0.145E+05 0.166E+05 0.182E+05 0.189E+05 0.183E+05 0.161E+05 0.118E+05 0.507E+04 I I I 0.433.3E0 9E05 16E05 2805 0143E05 56E+05 0.164E+05 0.162E+05 0.147E+05 0.115E+05 0.627E+04 I I I 0.217E+01 I 0.138E+05 0.119E+05 0.113E+05 0.117E+05 0.126E+05 0.137E+05 0.146E+05 0.149E+05 0.143E+05 0.122E+05 0.848E+04 I I I 0.0 I 0.171E+05 0.136E+05 0.119E+05 0.114E+05 0.118E+05 0.128E+05 0.139E+05 0.148E+05 0.150E+05 0.142E+05 0.120E+05 I I I I 0.0 0.394E+01 0.787E+01 0.118E+02 0.157E+02 0.197E+02 0.236E+02 0.276E+02 0.315E+02 0.354E+02 0.394E+02 I I X-COORDINATE I 0

TEMPERATURE GRADIENT WITH A TENSILE LOAD SIGMA Y II COORDINATE I I I I 0.108E+02 I -.990E+04 -.144E+05 -.69E+05 -.179E+05 -.178E+05 -.169E+05 -.156E+05 -.144E+05 -.136E+05 -.137E+05 -.149E+05 I I I 0.866E+01 I -.105E+05 -.142E+05 -.161E+05 -.166E+05 -.161E+05 -.151E+05 -.138E+05 -.127E+05 -.122E+05 -.127E+05 -.146E+05 I I I 0.650E+01 I -.103E+05 -.132E+05 -.146E+05 -.147E+05 -.140E+05 -.128E+05 -.116E+05 -.108E+05 -.107E+05 -.117E+05 -.142E+05 I I I 0.433E+01 I -.834E+04 -.107E+05 -.115E+05 -.113E+05 -.105E+05 -.929E+04 -.823E+04 -.767E+04 -.801E+04 -.963E+04 -.129E+05 I I 0.217E+01 I -.381E+04 -.558E+04 -.605E+04 -.560E+04 -.462E+04 -.352E+04 -.268E+04 -.250E+04 -.337E+04 -.568E+04 -.983E+04 I I I 0.0 I 0.422E+04 0.293E+04 0.279E+04 0.341E+04 0.440E+04 0.536E+04 0.590E+04 0.563E+04 0.415E+04 0.108E+04 -.399E+04 I I I I 0.0 0.394E+01 0.787E+01 0.118E+02 0.157E+02 0.197E+02 0.236E+02 0.276E+02 0.315E+02 0..354E+02 0.394E+02 I I X-COORDINATE I i-J

TEMPERATURE GRADIENT WITH A TENSILE LOAD SHEAR STRESS Y I COORDINATE I I I I 0.108E+02 I -.138E+03 -.149E+03 -.147E+03 -.130E+03 -.IO1E+03 -.571E+02 0.121E+00 0.711E+02 0.156E+03 0.254E+03 0.366E+03 I I I 0.866E+01 I -.596E+02 -.710E+02 -.741E+02 -.691E+02 -.558E+02 -.343E+02 -.453E+01 0.334E+02 0.796E+02 0.134E+03 0.197E+03 I I I 0.650E+01 I-. 125E+01 -.876E+01 -. 135E+02 -. 156E02 -. 149E+02 -.114E+02 -.523E+01 0.370E+01 0.154E+02 0.298E+02 0.469E02 I ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I I 0.433E+01 I 0.374E+02 0.377E+02 0.352E+02 0.300E+02 0.221E+02 0.114E+02 -.200E+01 -.182E+02 -.371E+02 -.587E+02 -.831E+02 I ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I I 0.217E+01 I 0.563E+02 0.684E+02 0.722E+02 0.678E+02 0.551E+02 0.343E+02 0.519E+01 -.321E+02 -.776E+02 -.131E+03 -.193E+03 I I 0.0 I 0.555E+02 0.832E+02 0.973E02 0.976E+02 0.842E+02 0.571E+02 O.163E+02 -.382E+02 -.106E+03 -.188E+03 -.284E+03 I I I I 0.0 0.394E+01 0.787E+01 0.118E+02 0.157E+02 0.197E+02 0.236E+02 0.276E+02 0.315E+02 0.354E+02 0.394E+02 I I X-COORDINATE I t(o Ix)

TEMPERATURE GRADIENT WITH A TENSILE LOAD MAXIMUM PRINCIPAL NORMAL STRESS Y I I COORDINATE I I I I I I I I 0.650E+01 I 0.901E+04 0.103E+05 0.123E+05 0.145E+05 0.166E+05 O.182E+05 0.189E+05 0.183E+05 O.161E+05 0.118E+05 0.507E+04 I 0.433E+01 I 0.113E+05 0.109E+05 0.116E+05 0.128E+05 0.143E+05 0.156E+05 0.164E+05 0.162E+05 0.147E+05 0.115E+05 0.627E+04 I I I 0.217 I 138E+05 19E05 13E7E+05 0.126E+05 0.137E+05 0.146E+05 0.149E+05 0.143E+05 0.122E+05 0.848E+04 I I I 0.0 I 0.171E+05 0.136E+05 0.119E+05 0.114E+05 0.118E+05 0.128E+05 0.139E+05 0.148E+05 0.150E+05 0.142E+05 0.120E+05 I I I I 0.0 0.394E+01 0.787E+01 0.118E+02 0.157E+02 0.197E+02 0.236E+02 0.276E+02 0.315E+02 0.354E+02 0.394E+02 I I X-COORDINATE I

PLATE VIBRATION WITH TEMPERATURE GRADIENT THE NUMBER OF SIMPLY SUPPORTED BEAM EIGENFUNCTIONS IS 4 THE NUMBER OF FREE-FREE BEAM EIGENFUNCTIONS IS 2 THE LENGTH OF THE BANDSAW BLADE IS 0.393700E+02 THE WIDTH OF THE BANDSAW BLADE IS O.108270E+02 THE THICKNESS OF THE BANDSAW BLADE IS 0.590550E-O1 THE MODULUS OF ELASTICITY IS 0.290100E+08 POISON"S RATIO IS 0.300000E+00 THE MASS DENSITY IS 0.283000E+OO THE COEFFICIENT OF THERMAL EXPANSION IS O 840000E-05 THE VELOCITY OF THE BANDSAW BLADE IS 0.0 THE WHEEL SUPPORT COEFFICIENT IS 0.0

PLATE VIBRATION WITH TEMPERATURE GRADIENT EIGEN REAL VALUE PART 1 -0.531418E+02 3 -0.299176E+02 5 -0.299176E+02 7 -0.299176E+02 9 -0.299176E+02 11 -0.299176E+02 13 -0.299176E+02 15 -0.299176E+02 LOWEST 5 EIGENVALUES IN IMAGINARY PART 0.0 0.844834E+02 0.233046E+02 0.318921E+02 0.714623E+02 0.468079E+02 0.489795E+02 0.648342E+02 INCREASING ORDER NON-DIMENSIONAL NATURAL FREQUENCY 0.455240E+03 0.767769E+03 0.324870E+03 0.374600E+03 0.663667E+03 0.475889E+03 0.491666E+03 0.611685E+03 5 7 1 11 13 LO U1

TENSILE. NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 THE DEGREE OF THE POLYNOMIAL FOR THE DISPLACEMENT FIELD CALCULATION IS 6 THE DEGREE OF THE POLYNOMIAL FOR THE TEMPERATURE FIELD CALCULATION IS 0 THE LENGTH OF THE BANDSAW BLADE IS 0.I00000E+01 THE WIDTH OF THE BANDSAW BLADE IS 0.250000E+00 THE THICKNESS OF THE BANDSAW BLADE IS. 107000E-02 THE MODULUS OF ELASTICITY IS 0.207000E+12 POISON"S RATIO IS 0.300000E+00 THE MASS DENSITY IS 0.283000E+00 THE VELOCITY OF THE BANDSAW BLADE IS 0.487400E+02 THE WHEEL SUPPORT COEFFICIENT IS 0.200000E+00 THE INITIAL TENSION IS 0.500000E+05 THE NORMAL CUTTING FORCE IS 0.250000E+03 THE TANGENTIAL CUTTING FORCE IS 0.500000E+03 THE STARTING X-COORDINATE FOR THE NORMAL CUTTING FORCE IS 0.487000E+00 THE ENDING X-COORDINATE FOR THE NORMAL CUTTING FORCE IS 0.513000E+00 THE STARTING X-COORDINATE FOR THE TANGENTIAL CUTTING FORCE IS 0.487000E+00 THE ENDINGX-COORDINATE FOR THE TANGENTIAL CUTTING FORCE IS 0.513000E+00 (JJ 0~n

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 COEFFICIENTS A( 1, O)= A( 2, 0)= A( 3, 0)= A( 4. 0)= A( 5, 0)= A( 6, 0)= A( 1, 1)= A( 2,1) A( 3, 1)= A( 4, 1)= A( 5, 1)= A( 1, 2)= A( 2, 2)= A( 3, 2)= A( 4, 2)= A( 1. 3)= A( 2, 3)= A( 3, 3)= A( 1, 4)= A( 2, 4)= A( 1, 5)= A( O. 1)= A( O, 2)= A( 0, 3)= A( 0, 4)= A( O, 5)= A( 0, 6)= OF APPROXIMATING FUNCTION FOR DISPLACEMENT FIELD U 0.83999E-03 0.16685E-03 -0.25155E-03 0.35664E-03 -0.35950E-03 0.14366E-03 0.10311E-02 -0.25467E-02 0.17339E-02 0.42432E-03 -0.55051E-03 -0.58239E-02 0.13222E-01 -0.10183E-01 0.20595E-02 0.13402E-01 -0.24878E-01 0.13451E-01 -0.91731E-02 0.66890E-02 0.32471E-02 -0.99129E-04 0.64076E-03 -0.17807E-02 0.24074E-02 -0.17646E-02 0.13378E-03 (AJ ',,

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 COEFFICIENTS B( 1, 0)= B( 2, Q)= B( 3. 0)= B( 4, 0)= B( 5, O)= B( 6, 0)= B( 1, 1)= B( 2, 1)= B( 3, 1)= B( 4, 1)= B( 5, 1)= B( 1, 2)= B( 2, 2)= B( 3, 2)= B( 4, 2)= B( 1, 3)= B( 2, 3)= B( 3, 3)= B( 1, 4)= B( 2, 4)= B( 1, 5)= B( O, 1)= B( 0. 2)= B( O, 3)= B( O, 4)= B( O, 5)= B( 0. 6)= OF APPROXIMATING FUNCTION FOR DISPLACEMENT FIELD V 0.97751E-05 -0.14137E-03 0.12220E-03 0.37973E-03 -0.63100E-03 0.26066E-03 -0.18996E-04 0.21425E-03 -0.84791E-03 0.11563E-02 -0.50844E-03 0.10928E-04 0.74170E-03 -0.10066E-02 0.28240E-03 -0.11924E-02 -0.29921E-03 0.12896E-02 0.38049E-02 -0.30997E-02 -0.80268E-03 -0.26631E-03 -0.71350E-04 0.44606E-03 -0.10174E-02 0.55038E-03 -0.16236E-03 (J 00

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 DISPLACEMENT FIELD U Y I I COORDINATE I I I I 0.250E+00 I -.484E-05 0.860E-04 0.176E-03 0.265E-03 0.353E-03 0.442E-03 0.531E-03 0.621E-03 0.712E-03 0.802E-03 0.893E-03 I I I 0.200E+00 I -.514E-05 0.854E-04 0.175E-03 0.264E-03 0.353E-03 0.443E-03 0.532E-03 0.622E-03 0.713E-03 0.803E-03 0.893E-03 I I I 0.150E+00 I -.538E-05 0.850E-04 0.175E-03 0.264E-03 0.353E-03 0.443E-03 0.533E-03 0.623E-03 0.714E-03 0.804E-03 0.894E-03 I I I 0.100E+00 I -.506E-05 0.847E-04 0.174E-03 0.264E-03 0.354E-03 0.443E-03 0.534E-03 0.624E-03 0.714E-03 0.804E-03 0.895E-03 I I I 0.500E-01 I -.356E-05 0.848E-04 0.174E-03 0.263E-03 0.353E-03 0.444E-03 0.534E-03 0.625E-03 0.715E-03 0.805E-03 0.896E-03 I I I 0.0 I 0.0 0.854E-04 0.173E-03 0.262E-03 0.353E-03 0.444E-03 0.535E-03 0.626E-03 0.716E-03 0.806E-03 0.896E-03 I I I I 0.0 0.100E+00 0.200E+00 0.300E+00 0.400E+00 0.500E+00 0.600E+00 0.700E+00 0.800E+00 O.900E+00 0.iOOE+01 I I X-COORDINATE I '4

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 DISPLACEMENT FIELD V Y I I COORDINATE I I I I 0.250E+00 I -.675E-04 -.681E-04 -.698E-04 -.721E-04 -.740E-04 -.748E-04 -.743E-04 -.727E-04 -.707E-04 -.690E-04 -.682E-04 I I I 0.200E+00 I -.540E-04 -.545E-04 -.564E-04 -.587E-04 -.60GE-04 -.613E-04 -.607E-04 -.590E-04 -.571E-04 -.555E-04 -.546E-04 I I I 0.150E+00 I -.405E-04 -.409E-04 -.428E-04 -.452E-04 -.471E-04 -.478E-04 -.471E-04 -.455E-04 -.436E-04 -.420E-04 -.410E-04 I I I 0.100E+00 I -.270E-04 -.273E-04 -.293E-04 -.317E-04 -.336E-04 -.342E-04 -.335E-04 -.319E-04 -.300E-04 -.286E-04 -.273E-04 I I I 0.500E-01 I -.134E-04 -.137E-04 -.157E-04 -.182E-04 -.200E-04 -.207E-04 -.200E-04 -.183E-04 -.165E-04 -.151E-04 -.136E-04 I I I 0.0 I 0.0 -.282E-06 -.230E-05 -.476E-05 -.656E-05 -.709E-05 -.632E-05 -.472E-05 -.299E-05 -.156E-05 -.542E-19 I I I I 0.0 0.-COORDINAT00E+00 0.20E00 0.300E+00 0.400E00 0.500E+00 0.600E+00.700E+0 0.800E+00.900E+00 0. E+01 I I X-COORDINATE I 0

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 COEFFICIENTS OF NORMAL STRESS IN X-DIRECTION AX(O,O)= 0.17290E+09 AX(1,0)= 0.74610E+08 AX(2,0)=-0.15704E+09 AX(3,0)= 0.26664E+09 AX(4,0)=-0.32997E+09 AX(5,0)= 0.16137E+09 AX(1,1)=-0.11571E+10 AX(2,1)= 0.12845E+10 AX(3,1)= 0.24870E+09 AX(4,1)=-0.58759E+09 AX(1,2)= 0.57713E+10 AX(2,2)=-0.70102E+10 AX(3,2)= 0.21379E+10 AX(1,3)=-0.10280E+11 AX(2,3)= 0.83331E+10 AX(1,4)= 0.27692E+10 AX(0.1)= 0.22482E+09 AX(0,2)=-0.12335E+10 AX(0,3)= 0.27709E+10 AX(0,4)=-0.18988E+10 AX(0,5)= 0.67216E+09 H

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 COEFFICIENTS OF NORMAL AX(0,0)=-0.32551E+07 AX(1,0)= 0.18451E+08 AX(2,0)=-0.27620E+07 AX(3,0)=-0.95527E+08 AX(4.0)= 0.14036E+09 AX(5,0)=-0.56837E+08 AX(1.1)=-0.34261E+09 AX(2,1)= 0.69241E+09 AX(3,1)=-0.34212E+09 AX(4,1)=-0.59365E+08 AX(1,2)= 0.99088E+09 AX(2,2)=-0.22889E+10 AX(3,2)= 0.14422E+10 AX(1,3)= 0.66624E+08 AX(2,3)=-0.66624E+08 AX(1,4)= 0.0 AX(0,1)= 0.37907E+08 AX(0,2)=-0.93037E+08 AX(0,3)=-0.11104E+08 AX(0,4)- 0.0 AX(0,5)= 0.0 STRESS IN Y-DIRECTION N)

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 COEFFICIENTS OF SHEAR STRESS AX(0O,)=-0.71140E+07 AX(1,O)= 0.59585E+08 AX(2,0)=-0.17357E+09 AX(3.0)= 0.25898E+09 AX(4,0)=-0.21740E+09 AX(5,0)= 0.80685E+08 AX(1,1)=-0.89323E+09 AX(2,1)= 0.19029E+10 AX(3,1)=-0.12532E+10 AX(4,1)= 0.12554E+09 AX(1,2)= 0.33192E+10 AX(2,2)=-0.61824E+10 AX(3.2)= 0.33027E+10 AX(1,3)=-0.29689E+10 AX(2.3)= 0.24382E+10 AX(1,4)= 0.79904E+09 AX(O,1)= 0.10052E+09 AX(0,2)=-0.42443E+09 AX(0,3)= 0.67173E+09 AX(0,4)=-0.39952E+09 AX(0,5)= 0.0 4-11

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 SIGMA X Y I I COORDINATE I I I I 0.250E+0O I 0.189E+09 0.187E+09 0.185E+09 0.183E+09 O.183E+09 0.184E+09 0.185E+09 0.187E+09 0.188E+09 0.188E+09 0.188E+09 I I I 0.200E+00 I 0.188E+09 0.187E+09 0.185E+09 0.184E+09 0.184E+09 0.185E+09 0.186E+09 0.187E+09 0.187E+09 0.187E+09 0.186E+09 I I I 0.150E+00 I 0.187E+09 0.186E+09 0.185E+09 0.185E+09 0.185E+09 0.186E+09 0.186E+09 0.187E+09 0.187E+09 0.187E+09 0.186E+09 I I I 0.100E+00 I 0.186E+09 0.186E+09 0.185E+09 0.185E+09 0.186E+09 0.186E+09 0.187E+09 0.187E+09 0.187E+09 0.187E+09 0.187E+09 I I I O.500E-01 I O.181E+09 0.184E+09 0.185E+09 0.186E+09 0.187E+09.8E+097 9 0..187E+09 87E09.188E7E+09 7E+09 E I I I 0.0 I 0.173E+09 0.179E+09 0.183E+09 0.186E+09 0.188E+09 0.189E+09 0.189E+09 0.188E+09 0.186E+09 0.186E+09 0.189E+09 I I I I 0.0 0.100E+00 0.200E+00 0.300E+00 0.400E+00 0.500E+00 0.600E+00 0.700E+00 0.800E+00 0.900E+00 0.100E+01 I I X-COORDINATE I

;ILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 SIGMA Y V I COORDINATE I I I 0.250E+00 I 0.233E+06 -.620E+04 -.108E+06 -.250E+06 -.448E+06 -.617E+06 -.649E+06 -.470E+06 -.119E+06 0.192E+06 0.526E+04 I I 0.200E+00 I 0.516E+06 -.132E+06 -.382E+06 -.495E+06 -.558E+06 -.556E+06 -.440E+06 -.189E+06 0.113E+06 0.215E+06 -.372E+06 I I I 0.150E+00 I 0.300E+06 -.350E+06 -.532E+06 -.544E+06 -.507E+06 -.427E+06 -.271E+06 -.299E+05 0.211E+06 0.204E+06 -.527E+06 I I 0.100E+00 I -.406E+06 -.657E+06 -.558E+06 -.402E+06 -.297E+06 -.232E+06 -.146E+06 0.401E+04 0.173E+06 0.162E+06 -.452E+06 I I I 0.500E-01 I -.159E+07 -.105E+07 -.458E+06 -.695E+05 0.664E+05 0.233E+05 -.683E+05 -.892E+05 0.446E+03 0.926E+05 -.136E+06 I 0.0 I -.326E+07 -.152E+07 -.233E+06 0.451E+06 0.581E+06 0.335E+06 -.422E+05 -.312E+06 -.306E+06 0.464E+03 0.426E+06 I I I I 0.0 0.100E+00 0.200E+00 0.300E+00 0.400E+00 0.500E+00 0.600E+00 0.700E+00 0.800E+00 0.900E+00 0.100E+01 I I X-COORDINATE I a(.J1

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 SHEAR STRESS Y I I COORDINATE I I I 0.250E+00 I 0.423E+06 0.142E+06 -.387E+06 -.798E+06 -.981E+06 -.981E+06 -.902E+06 -.809E+06 -.635E+06 -.804E+05 Q.148E+07 I I I 0.200E+00 I 0.747E+06 -.216E+06 -.997E+06 -.132E+07 -.119E+07 -.776E+06 -.309E+06 -.142E+05 0.992E+04 -.121E+06 0.241E+05 I I I 0.150E+00 I 0.479E+06 -.519E+06 -.128E+07 -.152E+07 -.126E+07 -.668E+06 -.421E+05 0.351E+06 0.342E+06 -.358E+05 -.452E+06 I I I 0.100E+00 I -.675E+06 -.898E+06 -.124E+07 -.131E+07 -.103E+07 -.506E+06 0.321E+05 0.361E+06 0.343E+06 0.250E+05 -.261E+06 I I I 0.500E-01 I -.307E+07 -.153E+07 -.915E+06 -.624E+06 -.384E+06 -.137E+06 0.573E+05 0.116E+06 0.294E+05 -.390E+05 0.341E+06 I I I 0.0 I -.711E+07 -.265E+07 -.390E+06 0.568E+06 0.784E+06 0.592E+06 0.189E+06 -.263E+06 -.545E+06 -.281E+06 0.116E+07 I I I I 0.0 0.100E+00 0.200E+00 0.300E+00 0.400E+00 0.500E+00 0.600E+00 0.700E+00 0.800E+00 0.900E+00 0.100E+01 I I X-COORDINATE I

TENSILE, NORMAL, TANGENTIAL LOAD DEGREE EQUAL 6 MAXIMUM PRINCIPAL NORMAL STRESS Y I I COORDINATE I I I I 0.250E+00 I 0.189E+09 0.187E+09 0.185E+09 0.183E+09 0.183E+09 O.184E+09 O.185E+09 0.187E+09 0.188E+09 0.188E+09 0.188E+09 I I I 0.200E+00 I 0.188E+09 0.187E+09 0.185E+09 0.184E+09 0.184E+09 0.185E+09 0.186E+09 0.187E+09 0.187E+09 0.187E+09 0.186E+09 I I I 0.1507E+09 0.186E+09 0.185E+09 0.185E+09 0.185E+09 0.186E+0 9 0.187E+09 0.187E+09 O.187E+09 O.186E+09 I I I 0.100E+00 I 0.186E+09 0.186E+09 0.185E+09 0.185E+09 0.186E+09 0.186E+09 0.187E+09 0.187E+09 0.187E+09 0.187E+09 0.187E+09 I I I 0.500E-01 I O.i81E+09 0.184E+09 0.185E+09 0.186E+09 0.187E+09 0.187E+09 0.187E+09 0.187E+09 0.187E+09 0.187E+09 0.188E+09 I I I 0.0 I 0.173E+09 0.179E+09 0.183E+09 0.186E+09 0.188E+09 0.189E+09 0.189E+09 0.188E+09 0.186E+09 0.186E+09 0.189E+09 I I I I 0.0 0.100E+00 0.200E+00 0.300E+00 0.400E+00 0.500E+00 0.600E+00 0.700E+00 0.800E+00 0.900E+00 0.100E+01 I I X-COORDINATE I

PROGRAM LISTING 48

MAIN PROGRAM 49

1 C 2 C 3 C 4 C 5 C 6 C 7 C 8 C 9 C 10 C 11 C 12 C 13 C 14 C 15 C 16 C 17 C 18 C 19 C 20 C 21 C 22 C 23 C 24 C 25 C 26 C 27 C 28 C 29 C 30 C 31 C 32 C 33 C 34 C 35 C 36 C 37 C 38 C 39 C 40 C 41 C 42 C 43 C 44 C 45 C 46 C 47 C 48 C 49 C 50 C 51 C 52 C 53 C 54 C 55 C 56 C 57 C 58 C 59 C 60 C 'BANDSAW' IS THE MAIN PROGRAM TO CALCULATE THE STRESS FIELD COEFFICIENTS, DISPLACEMENT FIELD COEFFICIENTS, TEMPERATURE FIELD COEFFICIENTS, AND NATURAL FREQUENCIES OF A GIVEN BANDSAW BLADE. THE INPUTS INTO THE PROGRAM ARE THE GEOMETRY OF THE BANDSAW BLADE, THE OPERATING CONDITIONS, TEMPERATURES AT VARIOUS POINTS ON THE BANDSAW BLADE, AND THE MATERIAL PROPERTIES OF THE BANDSAW BLADE. THE OPERATING CONDITIONS INCLUDES THE MILL SUPPORT COEFFICIENT, THE BLADE VELOCITY, THE INITIAL TENSILE LOAD, THE NORMAL AND TANGENTIAL CUTTING FORCES. A DISCUSSION OF THE THEORY LEADING TO THE EQUATIONS USED BY THE VARIOUS SUBROUTINES CALLED BY 'BANDSAW' SEE "DEVELOPMENT OF AN EFFICIENT COMPUTIONAL PROCEDURE FOR EVALUATING BAND SAW BLADE STRESS" BY J. E. BORCHELT, A. G. ULSOY, P. PAPALAMBROS (1983). EXAMPLE RUNS OF THIS PROGRAM ARE GIVEN IN "A COMPUTIONAL PACKAGE FOR DESIGN AND PROCESS EVALUATION IN BANDSAWING- A USER'S MANUAL". 'BANDSAW' CALLS SUBOUTINES TDRIVE, TOUT, FORCE, TEMP, EXPAND, STROUT, AND VIBRA. SUBROUTINE TDRIVE: READS TEMPERATURE DATA FOR VARIOUS POINTS ON THE BAND SAW BLADE. CALLS SUBROUTINE TCOFF SUBROUTINE TCOFF: CALCULATES TEMPERAURE FIELD COEFFICIENTS USING THE DATA READ IN SUBROUTINE TDRIVE BY SOLVING THE LINEAR REGRESSTION PROBLEM. TEMPERATURE FIELD COEFFICIENTS ARE USED IN SUBROUTINE TEMP. CALLS SUBROUTINES ATAM, DLUD, AND DBS. SUBROUTINE ATAM: CALCULATES THE MATRIX A-TRANSPOSED TIMES A FOR USE IN FINDING THE TEMPERATURE FIELD COEFFICIENTS. SUBROUTINE TOUT: PRODUCES A TABLE OF THE BLADE TEMPERATURES AT VARIOUS POINTS USING TEMPERATURE FIELD COEFFICIENTS CALCULATED IN SUBROUTINE TCOFF. SUBROUTINE FORCE: CALCULATES THE EFFECT OF THE TENSILE,NORMAL, AND/OR TANGENTIAL LOADS FOR THE RIGHTHAND SIDE VECTOR. ALSO DETERMINES NUMBER OF TERMS IN APPROXIMATION FUNCTION FOR DISPLACEMENT FIELD,(U). SUBROUTINE TEMP: CALCULATES EFFECT O(F A TEMPERATURE GRADIENT IF PRESENT FOR THE RIGHTHAND SIDE VECTOR. SUBROUTINE RHS: ADDS THE LOAD AND TEMPERATURE EFFECTS TOGETHER TO FORM THE RIGHTHAND SIDE VECTOR,WHICH IS STORED IN VECTOR BX. SUBROUTINE EXPAND: CALCULATES SUBMATRICES AS GIVEN IN REPORT. DETERMINES THE DISPLACEMENT FIELD COEFFICIENTS FOR U AND V. PRINTS OUT THE DISPLACEMENT FIELD COEFFICIENTS A AND B, IF DESIRED. CALLS SUBROUTINES MATDIM, DLUD, DBS, DOUT, AND SCOFF. SUBROUTINE DLUD: DECOMPOSES A GENERAL MATRIX(A) INTO LOWER AND UPPER TRIANGULAR MATRICES. THE DECOMPOSED MATRICES ARE STORED IN MATRIX(T). SUBROUTINE DBS: SOLVES A SYSTEM OF LINEAR EQUATIONS BY FORWARD ELIMINATION AND BACK SUBSTISUTION. SUBRUOTINE MATDIM: CALCULATES NONDIMENSIONAL INTEGRALS FOR USE IN SUBROUTINE EXPAND TO CALCULATE SUBMATRICES K1,K2,K3, AND K4. CALLS SUBROUTINE CMOD. SUBROUTINE CMOD: CALCULATES THE POWERS OF X AND Y USED TO DETERMINE THE ELEMENTS OF K1,K2,K3, AND K4. 0n 0

61 C 62 C 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 I 18 1 1A C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C SUBROUTINE SUBROUTINE SUBROUTINE SUBROUTINE SUBROUTINE SUBROUTINE SUBROUTINE SUBROUTINE SUBROUTINE SUBROUTINE BRIDGE SUBROUTINE THAT CALLS AMOD1, AMOD2, AMOD3,AMOD4,AMOD5,AMOD6,AMOD7. AMODI: CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1-NU**2)*(B/L) IN K1. AMOD2: CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1+NU)*(L/B) IN K1. AMOD3: CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1-NU) IN K2. AMOD4: CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY ((EH)NU)/(I-NU**2) IN K2. AMODS: CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1+NU) IN K2. AMOD6: CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1-NU**2)*(L/B) IN K4. AMOD7: CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1+NU)*(B/L) IN K4. DOUT: PRODUCES DISPLACEMENT FIELD TABLES FOR U AND V AT VARIOUS POINTS ON THE BLADE, USING THE DIPLACEMENT FIELD COEFFICIENTS DETERMINED IN SUBROUTINE EXPAND. SCOFF: CALCULATES THE STRESS FIELD COEFFICIENTS,FOR SIGMA-X, SIGMA-Y, AND TAU-XY, USING THE DISPLACEMENT FIELD COEFFICIENTS AND HOOKE'S LAW. STROUT: PRODUCES TABLES OF STRESSES AT VARIOUS POINTS ON THE BLADE. TABLES MADE ARE SIGMA-X, SIGMA-Y, SHEAR STRESS, AND THE MAXIMUM PRINCIPAL NORMAL STRESS. VIBRA: READS DIMENSIONLESS INTEGRALS FOR VIBRATION PROBLEM FROM FILE. CALCULATES MASS, GYRO, AND STIFFNESS MATRICES. CALLS SUBROUTINES DECODE, REFORM, AND EIGENS. DECODE: PICKS OUT THE RIGHT DIMENSIONLESS INTEGRAL TO BE USED IN CALCULATING THE MASS, GYRO, AND STIFFNESS MATRICES. REFORM: CHANGES THE SYSTEM OF SECOND ORDER DIFFERENTIAL EQUATIONS TO A SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS. EIGENS: CALCULATES THE NATURAL FREQUENCIES OF THE BAND SAW BLADE, -AND PRINTS OUT THE 5 LOWEST EIGENVALUES. SUBROUTINE Un SUBROUTINE SUBROUTINE SUBROUTINE ** ** ** IMPORTANT CONTROL VARIBLES USED ** IN PROGRAM ** ** ** ** IPROB: PROBLEM TYPE BEING SOLVED ** ** 1) PLANE STRESSES ONLY ** 2) VIBRATION PROBLEM ONLY ** 3) BOTH STRESS AND VIBRATION PROBLEM ** ** ** ** KOUT: INITIALLY ZERO SET EQUAL TO ONE IF ** ** SINGULAR MATRIX IS ECOUNTERED IN EITHER ** SUBROUTINE TCOFF, OR EXPAND. PROGRAM WILL ** ** STOP EXECUTION, IF KOUT=1! ** ** IS SET EAL ZERO IF A NEW ORDER OF * ** KON(1): IS SET EQUAL TO ZERO IF A NEW ORDER OF ** **

121 C ** SET EQUAL TO ONE UNLESS NEW ORDER IS** 122 C ENTERED. INITIALLY ZERO. USED IN ** 123 C **SUBROUTINE EXPAND. * 124 C ** ** 125 C * KON(2): INITIALLY ZERO. SET EQUAL TO ONE UNLESS ** 126 C NEW BLADE GEOMETRY, OR MATERIAL PROPERTIES ** 127 C ** ARE ENTERED. USED IN SUBROUTINE EXPAND. ** 128 C ** 129 C ** KON(3): INITIALLY ZERO. SET EQUAL ONE AFTER** 130 C PRINTING OF TEMPERATURE TABLE, UNLESS * 131 C ** NEW ORDER OF TEMPERATURE FIELD POLYNOMIAL ** 132 C IS ENTERED. PREVENTS REPEATED PRINTING ** 133 C OF TEMPERATURE FIELD TABLE DURING RUN WHEN ** 134 C ** THE TEMPERATURE DISTRIBUTION IS THE SAME. ** 135 C USED IN MAIN PROGRAM.** 136 C ** 137 C * KON(4): INITIALLY ZERO. SET EQUAL TO ONE IF ** 138 C ** PRINTING OF DISPLACEMENT FIELD ** 139 C **COEFFICIENTS IS NOT TO BE DONE. * 140 C **USED IN SUBROUTINE EXPAND.** 14 ** ** 142 C ** KON(5): INITIALLY ZERO. SET EQUAL TO ONE IF ** 143 C **PRINTING OF DISPLACEMENT TABLES IS NOT ** 1'44 C **TO BE DONE. USED IN SUBROUTINE EXPAND ** 145 C ** 146 C ** KON(6): INITIALLY ZERO. SET EQUAL TO ONE IF* 147 C PRINTING OF STRESS TABLES IS NOT TO ** 148 C ** BE DONE. USED IN MAIN PROGRAM.** 149 C ** ** 150 C ** KON(7): INITIALLY ZERO. SET EQUAL TO ONE IF** 151 C ** PRINTING OF TEMPERATURE TABLE IS NOT ** 152 C ** TO BE DONE. USED IN MAIN PROGRAM.** 153 C ** *+ 154 C ** KON(8): INITIALLY ZERO. SET EQUAL TO ONE AFTER * 155 C **READING OF DIMENSIONLESS INTEGRALS FOR ** 156 C VIBRATION ANALYSIS, SO DATA IS NOT* 157 C ** REREAD. USED IN SUBROUTINE VIBRA.* 158 C **** 159 C ** KON(9): INITIALLY ZERO. SET EQUAL TO ONE AFTER ** 160 C COMPLETION OF FIRST PROBLEM IN A RUN. ** 161 C ** USED IN MAIN PROGRAM.** 162 C ** ** 163 C * KON(10): INITIALLY ONE. SET EQUAL TO ZERO IF** 164 C **. NEW TITLE WAS ENTERED FOR NEW PROBLEM. * 165 C ** USED IN MAIN PROGRAM. * 166 C **** 167 C ** KON(11): USED TO CALL SUBROUTINE SDRIVE, WHICH ** 168 C **CALCULATES THE STRESS FIELD COEFFICIENTS ** 169 C ** NEEDED WHEN ONLY THE VIBRATION ANALYSIS ** 170 C ** DESIRED.** 171 C ** ** 172 C ** NEW: TELLS PROGRAM WHERE TO GO WHEN A CHANGE IN ** 173 C **THE PROBLEM HAS BEEN MADE.** 174 C 174 C ** ** 175 C*********** **************************** 176 C 176 C ***********************************+**********************+ 177 C 178 C 179 C 180 C U1 to

181 C 182 C 183 ***********C************C* ***+**+***************** 184 C ********************************+************+*** 185 ** ** 186 C ** NOTE: WHEN A STRESS AND VIBRATION PROBLEM ARE RUN ** 187 ** TOGETHER THE ORDER OF THE DISPLACEMENT FIELD ** 188 ** APPROXIMATING FUNCTIONS MUST BE 411 * 189 C ** THIS IS DONE AUTOMATICALLY FOR THE INITIAL * 190 C * COMBINED ANALYSIS. * 191 C ** * 192 C C ******************************************** 193 C ************************** ******************** * 194 C 195 C 196 IMPLICIT REAL*8(A-H,O-Z) 197 REAL*8 L.NU,KAPPA 198 DIMENSION PSAX(35,3),KON(11),ITITL(50) 199 DIMENSION BX(70),TX(35),AX(35,3),XI(2),X2(2) 200 DIMENSION TIX(35),TIY(35),FORCX(35),FORCY(35) 201 DATA NO /'N'/ 202 C 203 C THIS PART OF THE PROGRAM SETS THE CONTROL PARAMETERS 204 C KON. EQUAL TO ZERO. ESTABLISHES MAXN, MAXM, AND 205 C MAXNM. THE CONTROL FOR A SINGULAR MATRIX,KOUT, IS 206 C SET EQUAL TO ZERO, IF A SINGULAR MATRIX IS FOUND 207 C KOUT IS SET EQUAL TO 1. 208 MAXM=6 209 MAXN=6 210 MAXNM=36 211 KAPPA=O.OD+00 212 MAXT=O 213 C=O.OD+00 214 KOUT=O 215 DO 60 1=1,2 216 XI(I)=1.OD+00 217 60 X2(I)=2.0D+00 218 DO 65 1=1,9 219 65 KON(I)=O 220 KON(10)=1 221 C 222 C THIS PART OF THE PROGRAM BEGINS THE INTERACTIVE 223 C COMPUTER SESSION. FIRST THE TITLE IS READ INTO 224 C THE LABEL ITITL. 225 WRITE(6,610) 226 610 FORMAT(///,3X,'BEGIN PROBLEM DESCRIPTION') 227 WRITE(6,611) 228 611 FORMAT(/,6X,'ENTER THE TITLE OF THIS RUN (50 CHARACTERS MAX)') 229 READ(5,512) ITITL 230 512 FORMAT(50A1) 231 C 232 C THIS PART OF THE PROGRAM READS THE TYPE OF PROBLEM 233 C TO BE SOLVED. 234 2 WRITE(6,613) 235 613 FORMAT(/,6X,'ENTER NUMBER OF PROBLEM TO BE SOLVED',/, 236 $6X,'1-PLATE STRESSES'./, 237 $6X,'2-PLATE VIBRATION ',/, 238 $6X,'3-PLATE STRESSES AND VIBRATIONS') 239 READ(5.513) IPROB Ln wJ

241 IF(IPROB.EQ. 2.AND. KON(9).EQ. O) KON(11)=0 242 C 243 C THIS IF STATEMENT CHECKS TO MAKE SURE THAT A PROPER 244 C PROBLEM TYPE HAS BEEN CHOSEN. 245 IF(IPROB.GT. 3) GO TO 2 246 IF(IPROB.LE. O) GO TO 2 247 C 248 C THIS IF STATEMENT CHECKS TO SEE IF ONLY THE VIBRATION 249 C PROBLEM WAS CHOSEN. 250 IF(IPROB.EQ. 2) GO TO 20 251 C 252 C THIS PART OF THE PROGRAM ASKES THE USER THE DEGREE OF 253 C THE POLYNOMIAL TO BE USED IN THE STRESS AND TEMPERATURE 254 C FIELD COEFFICIENT CALCULATIONS. 255 5 WRITE(6,624) 256 624 FORMAT(//.3X,'BEGIN DESCRIPTION OF APPROXIMATING FUNCTIONS FOR 257 $,'DISPLACEMENT FIELD COEFFICIENT CALCULATION') 258 WRITE(6,614) 259 WRITE(6.638) 260 614 FORMAT(/,6X,'ENTER DEGREE OF THE POLYNOMIAL FOR THE APPROXIMATING 261 $,' FUNCTIONS') 262 638 FORMAT(6X,'INTEGER 2 TO 7') 263 READ(5,513) MAXS 264 513 FORMAT(II) 265 KON(1)=0 266 KON(2)=0 267 IF(KON(9).EQ. 1) GO TO 52 268 C 269 C THE IF STATEMENT CHECKS THAT MAXS IS SMALL ENOUGH SO 270 C ILL CONDITIONING OF COEFFICIENT MATRIX DOES NOT OCCUR. 271 IF(MAXS.GT. 7) MAXS=6 272 IF(MAXS.LT. 2) MAXS=2 273 615 FORMAT(/,6X,'ENTER DEGREE OF THE POLYNOMIAL FOR THE APPROXIMATING 274 $,' FUNCTION') 275 6 WRITE(6,636) 276 636 FORMAT(//.3X.'BEGIN DESCRIPTION OF APPROXIMATING FUNCTION FOR 277 $'TEMPERATURE FIELD COFFICIENT CALCULATION') 278 WRITE(6,615) 279 I2=MAXS-1 280 WRITE(6,637) 12 281 637 FORMAT(6X.'IF ZERO NO TEMPERATURE FIELD IS PRESENT',/, 282 $6X,'OTHERWISE ENTER INTEGER 1 TO',12) 283 READ(5,513) MAXT 284 IF(MAXT.EQ. O) KON(3)=1 285 IF(KON(9).EQ. 1.AND. MAXT.EQ. O) GO TO 52 286 C 287 C THE IF STATEMENT CHECKS TO SEE IF A TEMPERATURE FIELD 288 C IS PRESENT. 289 IF(MAXT.EQ. O) GO TO 30 290 IF(MAXT.GE. MAXS) MAXT=MAXS-1 291 CALL TDRIVE (MAXT,TX,KOUT) 292 KON(3)=0 293 C 294 C THIS IF STATEMENT AND OTHERS LIKE IT CHECK TO SEE IF 295 C A SINGULAR MATRIX WAS GENERATED IN ANY SUBROUTINE THAT 296 C USES THE MTS SYSTEM SUBROUTINES DLUD AND DBS. 297 IF(KOUT.EQ. 1) GO TO 57 298 IF(KON(9).EQ. 1) GO TO 52 299 C 300 C THIS IF STATEMENT CHECKS TO SEE IF ONLY THE STRESS FIELD Un VI

301 C COEFFICIENTS ARE TO BE CALCULATED. 302 30 IF(IPROB.EQ. 1) GO TO 10 303 C 304 C THIS PART OF THE PROGRAM ASKES THE USER THE NUMBER OF 305 C EIGENFUNCTIONS WHICH ARE TO BE USED IN THE NATURAL 306 C FREQUENCY CALCULATION. 307 20 WRITE(6,625) 308 625 FORMAT(//,3X,'BEGIN DESCRIPTION OF EIGENFUNCTIONS FOR 309 $'BANDSAW BLADE VIBRATIONS') 310 WRITE(6,634) 311 634 FORMAT(/,6X,'ENTER NUMBER OF SIMPLY SUPPORTED BEAM EIGENFUNCTIONS' 312* $,/,6X,'INTEGER 1 TO 6') 313 READ(5,513) N 314 WRITE(6,635) 315 635 FORMAT(/,GX,'ENTER NUMBER OF FREE-FREE BEAM EIGENFUNCTIONS',/, 316 $6X,'INTEGER 1 TO 6') 317 READ(5,513) M 318 C 319 C THIS IF STATEMENT CHECKS TO MAKE SURE THAT THE DEGREE OF THE 320 C CALCULATED STRESS FIELDS MATCHES THAT OF THE ONES REQUIRED IN 321 C THE VIBRATION ANALYSIS. 322 IF(KON(9).EQ. 0.AND. IPROB.EQ. 3) GO TO 825 323 IF(KON(9).EQ. 0.AND. IPROB.EQ. 2) GO TO 825 324 IF(KON(9).EQ. 1.AND. MAXS.EQ. 4) GO TO 825 325 WRITE(6,650) 326 GO TO 57 327 650 FORMAT(//,3X,'**** ERROR THE STRESS FIELD POLYNOMIALS DO NOT ', 328 1'MATCH THOSE REQUIRED FOR THE VIBRATION ANALYSIS! ***') 329 825 MAXS=4 330 IF(KON(9).EQ. 1) GO TO 27 331 C 332 C THIS PART OF THE PROGRAM ASKES THE USER FOR IFORMATION 333 C CONCERNING THE GEOMETRY OF THE PLATE (LENGTH, WIDTH, 334 C AND THICKNESS). 335 10 WRITE(6,'626) 336 626 FORMAT(//.3X,'BEGIN DESCRIPTION OF GEOMETRY OF BANDSAW BLADE') 337 WRITE(6,616) 338 616 FORMAT(/,6X,'ENTER LENGTH OF BANDSAW BLADE') 339 READ(5,514) L 340 514 FORMAT(D15.6) 341 WRITE(6,617) 342 617 FORMAT(/,6X,'ENTER BANDSAW BLADE WIDTH') 343 READ(5,514) B 344 WRITE(6,618) 345 618 FORMAT(/,6X,'ENTER BANDSAW BLADE THICKNESS') 346 READ(5,514) H 347 KON(2)=O 348 IF(KON(9).EQ. 1) GO TO 52 349 C 350 C THIS PART OF THE PROGRAM ASKES THE USER TO ENTER THE 351 C MATERIAL PROPERTIES OF THE BANDSAW BLADE (MODULUS OF 352 C ELASTICITY, POISSON'S RATIO, MASS DENSITY, AND 353 C COEFFICIENT OF THERMAL EXPANSION). 354 12 WRITE(6,627) 355 627 FORMAT(//,3X,'BEGIN DESCRIPTION OF MATERIAL PROPERTIES ', 356 $'OF THE BANDSAW BLADE') 357 WRITE(6,620) 358 620 FORMAT(/,6X,'ENTER MODULUS OF ELASTICTY (YOUNG"S MODULUS)') 359 READ(5,514) E tn U,

JV T I 10 1- r- u-~JtI I v / VA, Cl ciEcm rU133U i a IMa M I U I 362 READ(5,514) NU 363 WRITE(6,621) 364 621 FORMAT(/,6X,'ENTER THE MASS DENSITY OF THE MATERIAL') 365 READ(5,514) RHOH 366 RHI=RHOH 367 RHOH=RHOH*H 368 KON(2)=0 369 IF(KON(9).EQ. 1.AND. MAXT.EQ. O) GO TO 52 370 IF(MAXT.EQ. 0) GO TO 35 371 WRITE(6,640) 372 640 FORMAT(/,6X.'ENTER COEFFICIENT OF THERMAL EXPANSION') 373 READ(5,514) ALPHA 374 IF(KON(9).EQ. 1) GO TO 52 375 C 376 C THIS PART OF THE PROGRAM ASKES THE USER TO SUPPLY 377 C INFORMATION CONCERNING THE OPERATING CONDIONS OF 378 C THE BANDSAW. 379 35 WRITE(6,628) 380 628 FORMAT(//,3X,'BEGIN DESCRIPTION OF OPERATING CONDITIONS OF THE 381 $,'BANDSAW BLADE') 382 WRITE(6,622) 383 622 FORMAT(/,6X,'ENTER THE VELOCITY OF THE BANDSAW BLADE') 384 READ(5,514) C 385 WRITE(6.641) 386 641 FORMAT(/,6X,'ENTER WHEEL SUPPORT COEFFICIENT') 387 READ(5,514) KAPPA 388 IF(IPROB.EQ. 2) GO TO 27 389 WRITE(6.623) 390 623 FORMAT(/,6X,'ENTER THE INITIAL TENSION ON', 391 $' THE BANDSAW BLADE') 392 READ(5,514) RO 393 RIN=RO 394 C 395 C THIS STATEMENT CHANGES RO FROM A FORCE TO A FORCE PER UNIT LENGTH. 396 RO=RO/B 397 WRITE(6,629) 398 629 FORMAT(/,GX,'ENTER NORMAL CUTTING FORCE ON', 399 $' THE BANDSAW BLADE') 400 READ(5,514) FNOR 401 FNRI=FNOR 402 IF(FNOR.EO. 0.0) GO TO 43 403 WRITE(6,631) 404 631 FORMAT(/,9X,.'ENTER STARTING X-COORDINATE OF NORMAL CUTTING 405 $,'FORCE') 406 READ(5,514) X1(1) 407 WRITE(6,632) 408 632 FORMAT(9X,'ENTER ENDING X-COORDINATE OF NORMAL CUTTING FORCE') 409 READ(5,514) X2(1) 410 C 411 C THIS STATEMENT CHANGES FNOR FROM A FORCE TO A FORCE PER 412 C UNIT LENGTH. FNOR IS THE NORMAL CUTTING FORCE. 413 FNOR=FNOR/(X2(1)-X1(1)) 414 43 WRITE(6,630) 415 630 FORMAT(/,6X,'ENTER TANGENTIAL CUTTING FORCE ON', 416 $' THE BANDSAW BLADE') 417 READ(5,514) FTAN 418 FTNI=FTAN 419 IF(FTAN.EQ. 0.0) GO TO 27 420 WRITE(6,633) Ua 0%

421 633 FORMAT(/,9X,'ENTER STARTING X-COORDINATE OF TANGENTIAL CUTTING ', 422 $'FORCE') 423 READ(5,514) X1(2) 424 WRITE(6,647) 425 647 FORMAT(9X,'ENTER ENDING X-COORDINATE OF TANGENTIAL CUTTING ', 426 $'FORCE') 427 READ(5,514) X2(2) 428 C 429 C THIS STATEMENT CHANGES FTAN FROM A FORCE TO A FORCE PER 430 C UNIT LENGTH. FTAN IS THE TANGENTIAL CUTTING FORCE. 431 FTAN=FTAN/(X2(2)-XI(2)) 432 C 28 WRITE(6,827) 433 C 827 FORMAT(//,6X,'IS THERE A PRESTRESS ON THE BAND SAW BLADE (Y/N)?') 434 C READ(5,801) ICPS 435 C 436 C THIS IF STATEMENT CHECKS TO SEE IF A PRESTRESS FIELD IS PRESENT. 437 C IF(ICPS.EQ. NO) GO TO 27 438 C WRITE(6,828) 439 C 828 FORMAT(/,9X,'ENTER ORDER OF THE POLYNOMIAL FOR PRESTRESS FIELD 440 C 1,'CALCULATION',/,9X,'INYERGER 1 TO 6') 441 C READ(5,513) NPPS 442 C 443 C HERE SDRIVE CALCULATES THE PRESTRESS FIELDS. '444 C SDRIVE DOES NOT WORK PROPERLY AT THIS TIME. 445 C CALL SDRIVE(PSAX,NPPS,KOUT) 446 C IF(KOUT.EQ. 1) GO TO 57 44'7 27 IF(KON(9).EQ. 1) GO TO 52 448 C 449 C THIS PART OF THE PROGRAM ASKS THE USER WHAT OUTPUT IS TO 450 C BE WRITTEN TO THE OUTPUT FILE 451 WRITE(6,802) 452 802 FORMAT(///,3X,'BEGIN OUTPUT CONTROL DESCRIPTION') 453 IPRD=4 454 IF(MAXT.EQ. O) IPRD=3 455 DO 810 IPT=1,IPRD 456 IF(IPT.EQ. 1) WRITE(6,800) 457 IF(IPT.EQ. 2) WRITE(6,803) 458 IF(IPT. EQ. 3) WRITE(6,804) 459 IF(IPT.EQ. 4) WRITE(6,805) 460 IKON=IPT+3 461 READ(5,801) IDP 462 IF(IDP.EQ. NO) KON(IKON)=1 463 800 FORMAT(/,6X,'DO YOU WISH THE PRINTING OF THE DISPLACEMENT', 464 1' FIELD COEFFICIENTS (Y/N)?') 465 803 FORMAT(/,6X,'DO YOU WISH THE PRINTING OF THE DISPLACEMENT', 466 1' FIELD TABLES (Y/N)?') 467 804 FORMAT(/,6X,'DO YOU WISH THE PRINTING OF THE STRESS FIELD', 468 1' TABLES (Y/N)?') 469 805 FORMAT(/,6X,'DO YOU WISH THE PRINTING OF THE TEMPERATURE', 470 1' FIELD COEFFICIENTS/TABLE (Y/N)?') 471 801 FORMAT(IAI) 472 810 CONTINUE 473 C 474 C THIS PART OF THE PROGRAM WRITES TO THE OUTPUT FILE, DEVICE 475 C #8, THE USER INPUT FROM THE INTERACTIVE SESSION. 476 25 WRITE(8,710) ITITL 477 IF(IPROB.EQ. 2) GO TO 80 478 WRITE(8,720) MAXS,MAXT 479 IF(IPROB.EQ. 1) GO TO 81 U,

401 01 WK 1 It \, I 4U L, t,M 482 WRITE(8,750) E,NU,RHI 483 IF(MAXT.GT. O) WRITE(8,760) ALPHA 484 WRITE(8,770) C,KAPPA 485 IF(IPROB.EQ. 2) GO TO 55 486 WRITE(8,780) RIN,FNRI,FTNI 487 IF(FNRI.NE. 0.0) WRITE(8,790) Xl(1),X2(1) 488 IF(FTNI.NE. 0.0) WRITE(8,795) X1(2),X2(2) 489 710 FORMAT(1H1,2(/,IH-),40X,50A1) 490 720 FORMAT(1HO,5X,'THE DEGREE OF THE POLYNOMIAL FOR THE DISPLACEMENT' 491 1,' FIELD CALCULATION IS',I2,/,1HO,5X,'THE DEGREE OF THE ', 492 2'POLYNOMIAL FOR THE TEMPERATURE FIELD CALCULATION IS',12) 493 730 FORMAT(1H0,SX,'THE NUMBER OF SIMPLY SUPPORTED BEAM ', 494 1'EIGENFUNCTIONS IS',I2,/,IHO,5X,'THE NUMBER OF FREE-FREE ', 495 2'BEAM EIGENFUNCTIONS IS',I2) 496 740 FORMAT(1H0,5X,'THE LENGTH OF THE BANDSAW BLADE IS',D15.6,/, 497 11HO,5X,'THE WIDTH OF THE BANDSAW BLADE IS',D15.6,/,1HO,5X, 498 2'THE THICKNESS OF THE BANDSAW BLADE IS',D15.6) 499 750 FORMAT(1HO,5X,'THE MODULUS OF ELASTICITY IS',D15.6,/,1HO,5X, 500 1'POISON"S RATIO IS',D15.6,/1HO,5X,'THE MASS DENSITY IS',D15.6) 501 760 FORMAT(1HO,5X,'THE COEFFICIENT OF THERMAL EXPANSION IS',D15.6) 502 770 FORMAT(IHO,5X,'THE VELOCITY OF THE BANDSAW BLADE IS',D15.6, 503 1/,1HO,SX,'THE WHEEL SUPPORT COEFFICIENT IS',D15.6) 504 780 FORMAT(1HO,5X.'THE INITIAL TENSION IS',D15.6,/,1HO,5X, 505 1'THE NORMAL CUTTING FORCE IS',D15.6,/,1HO,SX,'THE TANGENTIAL', 506 2' CUTTING FORCE IS',D15.6) 507 790 FORMAT(IHO,SX,'THE STARTING X-COORDINATE FOR THE NORMAL' 508 1,' CUTTING FORCE IS',D15.6./1HO,5X,'THE ENDING X-COORDINATE', 509 2' FOR THE NORMAL CUTTING FORCE IS',D15.6) 510 795 FORMAT(1HO,5X,'THE STARTING X-COORDINATE FOR THE TANGENTIAL ' 511 1'CUTTING FORCE IS',D15.6,/,1HO,5X,'THE ENDING X-COORDINATE FOR', 512 2' THE TANGENTIAL CUTTING FORCE IS',D15.6) 513 C 514 C THIS EQUATION CALCULATES THE RUNNING TENSION ON THE BLADE. 515 R=RO-RHOH*KAPPA*C**2 516 IF(MAXT.EQ. O) GO TO 40 517 IF(KON(3).EQ. 1) GO TO 40 518 IF(KON(7).EQ. 1) GO TO 40 519 CALL TOUT(MAXT,L,B,TX,ITITL) 520 KON(3)=1 521 C 522 C SUBROUTINES FORCE,TEMP,RHS, AND EXPAND ARE USED 523 C TO CALCULATE THE STRESS FIELD COEFFIECIENTS. 524 40 CALL FORCE(MAXS,R.FTAN,FNOR,XI,X2,L,B,FORCXFORCYMSIZE) 525 CALL TEMP(MAXS,MAXT,TX,MSIZE,L,B,H,E,NU,ALPHA,TIX,TIY) 526 CALL RHS(FORCX,FORCY.TIX.TIY.MSIZE,MAXS,BX) 527 CALL EXPAND(MAXS,MSIZE,L,B,H,E,NU,BX,AX,KOUT,KON,ITITL) 528 IF(KOUT.EQ. 1) GO TO 57 529 IF(KON(6).EQ. 1) GO TO 50 530 CALL STROUT(MAXS,L,B,AX,ITITL) 531 IF(IPROB.EQ. 1) GO TO 50 532 55 CONTINUE 533 IF(N.GT. MAXN) N=MAXN 534 IF(M.GT. MAXM) M=MAXM 535 NM=N*M 536 C 537 C D IS THE FLEXURE MODULUS OF THE PLATE, WHICH IS USED IN 538 C THE VIBRATION ANALYSIS. 539 D=(E*H**3)/((1.2D+01)*((1.0D+00)-NU**2)) 540 IF(IPROB.NE. 2) GO TO 77 UL o0

541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 C 557 C 558 C 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 End of file IF(KON(11).EQ. O) GO TO 76 WRITE(6,850) 850 FORMAT(///,3X,'DO YOU WISH TO USE THE STRESS FIELD COEFFICIENTS' 1' CALCULATED PREVIOUSLY (Y/N)?') READ(5,801) ISDS IF(ISDS.EQ. NO) KON(11)=0 IF(KON(11).EQ. 1) GO TO 77 76 NCSF=3 IF(KON(11).EQ. O) CALL SDRIVE(AX,NCSF,KOUT) KON(11)=1 IF(KOUT.EQ. 1) GO TO 57 77 CONTINUE CALL VIBRA(N,M,NM,MAXNM,KON,L,B,H,D,RHOH,NU,CAXITITL) 50 CONTINUE KON(9)=1 THIS PART OF THE PROGRAM ASKS THE USER WHAT IS TO BE DONE NEXT. ANOTHER PROBLEM, OR STOP. 52 WRITE(6,700) 700 FORMAT(//,3X,'ENTER NUMBER OF NEXT INSTRUCTION',/, $6X,'1- NEW PROBLEM TYPE',/,6X,'2- NEW APPROXIMATING', $' FUNCTION FOR DISPLACEMENTS',/,6X,'3- NEW APPROXIMATING', $' FUNCTION FOR TEMPERATURE',/,6X.'4- NEW EIGENFUNCTIONS FOR $'VIBRATION ANALYSIS',/,6X,'5- NEW BLADE GEOMETRY', $/,6X,'6- NEW MATERIAL PROPERTIES',/,6X,'7- NEW', $' OPERATING CONDITIONS',/,6X,'8- SOLVE NEW PROBLEM' $,/,6X,'9- STOP') READ(5,513) NEW IF(NEW.EQ. 9) GO TO 57 IF(NEW.EQ. 8) KON(10)=0 IF(KON(10).NE. O) GO TO 62 WRITE(6,611) READ(5,512) ITITL 62 KON(1O)=l GO TO(2,5,6,20,10,12,35,25),NEW 57 STOP END Cn 0l k0

STRESS ANALYSIS SUBROUTINES 60

1 SUBROUTINE EXPAND(N,MSIZE,EL,B,H,E,UN.BX,AX,KOUT,KON,ITITL) 2 C 3 C EXPAND CALCULATES THE DIMENSIONIZED MATRIX (A), WHICH IS 4 C USED TO SOLVE FOR THE COEFFICIENTS OF THE DISPLACEMENT FIELD 5 C APPROXIMATING FUNCTIONS FOR U AND V. THE SUBMATRICES Ki, K2 6 C, K3, AND K4 ARE CALCULATED USING THE NONDIMENSIONAL INTEGRALS 7 C CALCULATED BY SUBROUTINE AMOD. THE MATRIX (A) IS THEN ASSEMBLED 8 C WITH THE ROW CORRESPONDING TO THE VARION OF B(N,O), AND THE 9 C COLUMN CORRESPONDING TO B(N,O) BEING ELIMINATED. 10 C 11 C EXPAND CALLS SUBROUTINES MATOIM, DLUD, DBS, DOUT, AND SCOFF. 12 C SUBROUTINE MATDIM: CALCULATES THE NONDIMENSIONAL INTEGRALS 13 C FOR USE IN CALCULATING THE SUBMATRICES 14 C K1, K2, K3, AND K4. 15 C SUBROUTINE DLUD: DECOMPOSES THE MATRIX (A) INTO LOWER AND 16 C UPPER TRIANGULAR MATRICES. 17 C SUBROUTINE DBS: CALCULATES THE UDETERMINED COEFFICIENTS A, 18 C AND B BY FORWARD ELIMATION AND BACK SUBSTITUTION. 19 C SUBROUTINE DOUT: PRODUCES DISPLACEMENT FIELD TABLES FOR 20 C U AND V AT VARIOUS POINTS ON THE BLADE. 21 C SUBROUTINE SCOFF: CALCULATES THE STRESS FIELD COEFFICIENTS 22 C USING THE DISPLACEMENT FIELD COEFFICIENTS 23 C AND HOOKE"S LAW. 24 C 25 C ******************************************************* 26 C ********************************************************** 27 C ** ** 28 C ** IMPORTANT VARIBLES IN SUBROUTINE EXPAND 29 C ** ** 30 C EL: BLADE LENGTH BETWEEN SUPPORTS. * 31 C ** VARIBLE (L) IN MAIN PROGRAM AND SUBROUTINES t 32 C ** VIBRA, AND EIGENS. ** 33 C ** ** 34 C ** B: BLADE WIDTH. ** 35 C ** ** 36 C ** H: BLADE THICKNESS. ** 37 C ** 38 C ** E: MODULUS OF ELASTICITY. 39 C ** ** 40 C ** UN: POISSON"S RATIO. VARIBLE (NU) IN MAIN ** 41 C ** PROGRAM AND SUBROUTINE VIBRA. ** 42 C ** 43 C ** N: ORDER OF APPROXIMATING FUNCTIONS FOR ** 44 C ** DISPLACEMENT FIELDS U AND V. ** 45 C ** 46 C ** MSIZE: NUMBER OF TERMS IN APPROXIMATING FUNCTION ** 47 C ** FOR DISPLACEMENT FIELD (U). 48 C ** * 49 C ** NSIZE: NUMBER OF TERMS IN APPROXIMATING FUNCTION 50 C ** FOR DISPLACEMENT FIELD (V). ** 51 C ** ** 52 C ** d2: SIZE OF THE MATRIX (A). MSIZE+NSIZE. 53 C 53 C ** ** 54 C ** NM: MAXIMUM SIZEOF THE MATRIX (A). ** 55 C ** ** 56 C * KON(1).: O- NEW ORDER OF APPROXIMATING FUNCTIONS. ** 57 C ** 1- ORDER OF APPROXIMATING FUNCTIONS DID ** 58 C ** NOT CHANGE. 59 C ** ** 60 C ** KON(2): 0- NEW BLADE GEOMETRY, MATERIAL PROPERTIES ** I

61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C c * ** ** * *, OR ORDER OF APPROXIMATING FUNCTIONS. 1- NO CHANGE IN BLADE GEOMETRY, MATERIAL PROPERTIES, OR ORDER OF APPROXIMATING FUNCTIONS. ** ** KON(4): O- PRINT OUT DISPLACEMENT FIELI ** COEFFICIENTS FOR U AND V TO ** 1- DO NOT PRINT OUT DISPLACEMEI ** COEFFICIENTS TO FILE. * ** KON(5): O- PRODUCE DISPLACEMENT FIELD * ~* 1- DO NOT PRODUCE DISPLACEMENT ** TABLES * * * COFF: DISPLACEMENT FIELD COEFFICIENTS. ** ** AX: STRESS FIELD COEFFICIENTS FOR SIGM ** SSIGMA, AND TAU-XY. ** D FILE. NT FIELD TABLES. FIELD 9* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** '9* ** ** A-X, ** ***IMPLICIT REA*8(A-H,O-Z)********** IMPLICIT REAL*8(A-HO-Z) REAL*8 K1,K2,K3,K4 COMMON /SUB2/ XP,YP,XP1,YPI COMMON /SUB4/ AMOD DIMENSION ITITL(50),AMOD(7,1230),XPi(1230),YPI(1230),XP(730) DIMENSION YP(730),A(70,70),KI(35,35),K2(35,35),K3(35,35),K4(35,35) DIMENSION KON(11),AX(35,3),T(70,70),IV(70),BX(70),COFF(75) INITIALIZE THE STRESS COEFFICIENTS TO ZERO DO 52 I=1.3 DO 53 d=1,35 53 AX(JI)=O.OD+00 52 CONTINUE NM IS THE MAXIMUM MATRIX SIZE FOR SUBROUTINES DLUD AND DBS NM=70 MSIZE IS THE DIMENSION OF MATRICES K1,K2,K3,K4 ALSO, MSIZE IS THE NUMBER OF TERMS IN THE APPROXIMATING FUNCTION FOR U. JJ=MSIZE NSIZE IS THE NUMBER OF TERMS IN THE APPROXIAMATING FUNCTION FOR V WHEN C(N,O) IS REMOVED NSIZE=MSIZE-1 J2 IS THE DIMENSION OF THE COFFECIENT MATRIX A. J2=dJ+NSIZE SUBROUTINE MATDIM FORMS THE NONDIMENSIONAL INTEGRATIONS OF THE VARIATIONAL EQUATION FOR U AND V. THIS SET OF IF STATEMENTS CHECKS TO SEE IF THE NONDIMENSIONAL INTEGRALS NEED TO BE RECALCULATED. IF(KON(1).EQ. 0 ) GO TO 73 IF(KON(1).EQ. 1) GO TO 68 N) to>

121 KON(1)-1 122 C 123 C THIS IF STATEMENT CHECKS TO SEE IF NEW BLADE GEOMETRY, OR NEW 124 C MATERIAL PROPERTIES HAVE BEEN ENTERED. 125 68 IF(KON(2).EQ. 1 ) GO TO 65 126 JK=1 127 DO 15 I=1,MSIZE 128 DO 15 J=I,MSIZE 129 C 130 C K1 IS THE COEFFICIENT OF A(N,M) MULTIPLIED BY THE 131 C VARIATION OF A(N,M) 132 K1(I,J)=(E*H/((1.0D+0)-UN**2)*AMOD(1iJK)*(B/EL) 133 1+E*H/((1.0D+0)+UN)*AMOK)AMO(2K)*(EL/B))*EL**XP(K)* 134 2B**YP(JK) 135 K1(d,I)=K1(I,J) 136 C 137 C K4 IS THE COEFFICIENT OF B(N,M) MULTIPLIED BY THE 138 C VARIATION OF B(N,M) 139 K4(I,d)=(E*H/((1.OD+00)-UN**2)*AMOD(6,K)*(EL/B)+E*H/((1.OD+00) 140 1+UN)*AMOD(7,dK)*(B/EL))EL EL**XP(JK)*B**YP(K) 141 K4(, I)=K4(I,J) 142 15 JK=JK+1 143 JK-1 1'44 DO 16 I=1,MSIZE 145 DO 16 J=1.MSIZE 146 C 147 C K2 IS THE COEFFICIENT OF B(N,M) MULTIPLIED BY THE 148 C VARIATION OF A(N,M) 149 K2(I, J)(E*H/((.00+00)-UN)*AMOD(3,dK)-E*H*UN/((1.OD+00)-UN**2)* 150 1AMOD(4,JK)-E*H/((1.OD+00)+UN)*AMOD(5.UK))*EL**XP1(JK)*B**YP1(JK) 151 C 152 C K3 IS THE COFFICIENT OF A(N,M) MULTIPLIED BY THE 153 C VARIATION OF B(N,M), K3 IS EQUAL TO K2 TRANSPOSED. 154 K3(d,I)=K2(I,J) 155 16 JK=JK+1 156 C 157 C THIS NESTED DO LOOP ASSEMBLES THE COEFFICIENT MATRIX,A, 158 C ELIMINATING THE COLUMN CORRESPONDING TO B(N,O) AND ROW 159 C CORRESPONDING TO THE VARIATION OF B(N,O) 160 DO 25 I1.,MSIZE 161 IA=I+dd 162 DO 26 J=1,MSIZE 163 d1=d 164 11=I 165 dA=J+Jd 166 A(I,d)=KI(I1,1 ) 167 IF(J1.EQ. MSIZE) GO TO 30 168 IF(d.GE. N) d1=J1+1 169 A(I,dA)=K2(I1,J1) 170 30 d1=d 171 IF( I1.EQ. MSIZE) GO TO 26 172 IF(I.GE. N) 11=11+1 173 A(IA,J)=K3(I1,d1) 174 IF(d1.EQ. MSIZE) GO TO 25 175 IF(d.GE. N) d1=d1+1 176 A(IA,JA)-K4(I1,J1) 177 26 CONTINUE 178 25 CONTINUE 179 C 180 C SUBROUTINE DLUD IS A SYSTEM SUBROUTINE THAT DECOMPOSES THE O) LA.

181 C MATRIX(A) INTO LOWER AND UPPER TRIANGLULAR MATRICES. 182 C THE DECOMPOSED MATRICES ARE STORED IN MATRIX(T). 183 CALL DLUD(02,NM,A,NMT.IV) 184 IF(IV(d2).EQ. O) GO TO 50 185 KON(2)=1 186 C 187 C SUBROUTINE DBS IS A SYSTEM SUBROUTINE THAT SOLVES 188 C THE SYSTEM OF EQUATIONS [A](X}={BX} BY FORWARD 189 C ELIMINATION AND BACK SUBSTITUTION 190 65 CONTINUE 191 CALL DBS(J2,NM,T,IV,BX) 192 C 193 C INITIALIZE DISPLACEMENT FIELD COEFFICIENTS FOR U AND V TO ZERO 194 DO 62 d=1,45 195 62 COFF(J)=O.OD+00 196 JM=d2+1 197 MINUS=O 198 DO 75 d=1,dM 199 J1=d 200 IF(MINUS.EQ. 1) d1=J1-1 201 IF(d.EQ. MSIZE+N) GO TO 60 202 C?03 C COFF ARE THE COEFFICIENTS A(N,M) FOR DISPLACEMENT FIELD U, 204 C AND B(N,M) FOR DISPLACEMENT FIELD V. 205 COFF(d)=BX(J1) 206 GO TO 75 207 C 208 C THIS DO LOOP CALCULATES B(N,O) KNOWING 209 C B(1,0),B(2,0),....B(N-1,0) 210 60 DO 85 I=1,N 211 IF(N-I.EQ. O) GO TO 75 212 J3=J1-N+I 213 F=-BX(d3)/EL**(N-I) 214 COFF(J)=F+COFF(d) 215 MINUS=1 216 85 CONTINUE 217 75 CONTINUE 218 C 219 C SUBROUTINE ZEROC IS USED TO, SET SMALL DISPLACEMENT FIELD 220 C COEFFICIENTS TO ZERO. 221 MD=d2+1 222 CALL ZEROC(COFF,MD,N,EL,B) 223 WRITE(8,500) ITITL 224 WRITE(8,101) 225 C 226 C THIS NESTED DO LOOP WRITES OUT THE COEFFICIENTS 227 C FOR THE APPROXIMATING FUNCTIONS U AND V 228 C 229 C THIS IF STATEMENT CHECKS TO SEE IF THE DISPLACEMENT COEFFICIENTS 230 C ARE TO BE PRINTED TO THE OUTFILE. 231 IF(KON(4).EQ. 1) GO TO 79 232 IJ=1 233 DO 95 J=1,2 234 IF(J.EO. 2) WRITE(8,500) ITITL 235 IF(d.EQ. 2) WRITE(8,102) 236 DO 105 I1=1,N 237 I=I1-1 238 DO 115 J1=1,N 9r~o r OC at 4^

241 C COEFFICIENT, AND I IS THE ORDER OF Y ASSOCIATED WITH EACH 242 C DISPLACEMENT FIELD COEFFICIENT. 243 IF(J1+I.GT. N) GO TO 105 244 IF(J.EQ. 1) WRITE(8,100) d1,I,COFF(Id) 245 IF(J.EQ. 2) WRITE(8,103) d1I,,COFF(IJ) 246 115 IJ=1J+1 247 105 CONTINUE 248 1=0 249 DO 120 d1=1,N 250 IF(J.EQ. 1) WRITE(8,100) I,d1,COFF(Id) 251 IF(J.EQ. 2) WRITE(8,103) I,d1,COFF(IJ) 252 120 Id=IJ+1 253 95 CONTINUE 254 100 FORMAT(8X,'A(',12,',',12,')=',D14.5) 255 101 FORMAT(1H-,4X,'COEFFICIENTS OF APPROXIMATING FUNCTION FOR' 256 $,' DISPLACEMENT FIELD U') 257 102 FORMAT(1H-,4X,'COEFFICIENTS OF APPROXIMATING FUNCTION FOR' 258 $,' DISPLACEMENT FIELD V') 259 103 FORMAT(8X,'B(',12,',',12,')=',D14.5) 260 500 FORMAT(1H1,2(/,1H-),40X,50A1) 261 C 262 C THIS IF STATEMENT CHECKS TO SEE IF THE DISPLACEMENT FIELD 263 C TABLES ARE TO BE PRODUCED. 264 79 IF(KON(5).EQ. 1) GO TO 97 265 C 266 C DOUT PRODUCES DISPLACEMENT FIELD TABLES OF U AND V. 267 CALL DOUT(N.MSIZE,EL,B,COFF,ITITL) 268 C 269 C SCOFF CALCULATES THE STRESS FIELD COEFFICIENTS, AX, FOR SIGMA-X, 270 C SIGMA-Y, AND TAU-XY USING THE DISPLACEMENT FIELD COEFFICIENTS 271 C AND HOOKE'S LAW. 272 97 CALL SCOFF(EL,B,N,MSIZE,COFF,E,UN,AXITITL) 273 GO TO 110 274 C 275 C SINGULAR MATRIX ENCOUNTERED. 276 50 WRITE(6,106) 277 KOUT=1 278 106 FORMAT(4X,'*** ERROR ** MATRIX IS SINGULAR ', 279 $'IN SUBROUTINE EXPAND') 280 110 RETURN 281 END 282 SUBROUTINE ZEROC(ZCOF,MD,N,EL,B) 283 C 284 C ZEROC SETS SMALL DISPLACEMENT COEFFICIENTS TO ZERO. EACH COEFFICIENT 285 C IS MULTIPLIED BY THE LENGHT(EL) RAISED TO THE APPROPIATE POWER 286 C AND THE WIDTH(B) RAISED TO THE APPROPIATE POWER AND STORED IN ZCF. 287 C THE MAXIMUM MAGNITUDE OF ZCF IS THEN FOUND. EACH VALUE OF ZCF IS 288 C THEN DIVIDED BY THIS MAXIMUM VALUE, IF THE RESULT IS LESS THAN 1.OE-6. 289 C THE COEFFICIENT ZCOF CORRESPONDING TO ZCF IS SET EQUAL ZERO. ZCOF IS 290 C TRANSFERED BACK TO SUBROUTINE EXPAND AS COFF. 291 IMPLICIT REAL*8(A-H,O-Z) 292 DIMENSION ZCF(75),ZCOF(75) 293 JBL=1 294 C 295 C THIS SET OF DO LOOPS MULTIPLIES THE DISPLACEMENT FIELD COEFFICIENTS 296 C ZCOF BY EL AND B RAISED TO THE APPROPIATE POWER. 297 30 DO 10 K=1,N 298 ZCF(JBL)=ZCOF(JBL)*EL**K 299 10 dBL=JBL+1 300 DO 15 L'1,N cn U'I

301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 End of file C DO 20 K1=,N IF(K+L.GT. N) GO TO 15 ZCF(JBL)=ZCOF(JBL)*EL**K*B**L 20 dBL=JBL+1 15 CONTINUE DO 25 L=I,N ZCF(JBL)=ZCOF(dBL)+B**L 25 JBL=JBL+1 IF(JBL.LT. MD) GO TO 30 THIS DO LOOP FINDS THE MAXIMUM VALUE OF ZCF. ZMAX=DABS(ZCF(1)) DO 35 I=2,MD ZC=DABS(ZCF(I)) IF(ZC.GT. ZMAX) ZMAX=ZC 35 CONTINUE THIS DO LOOP SETS SMALL VALUES OF THE DISPLACEMENT FIELD COEFFICIENTS ZCOF EQUAL TO ZERO. DO 40 I=1,MD ZC=DABS(ZCF(I))/ZMAX IF(ZC.LT. 1.OD-6) ZCOF(I)=O.OD+00 40 CONTINUE RETURN END C 0cn

326 SUBROUTINE MATDIM(N) 327 C 328 C ************+*************************** ********+*****+* 329 C ********************************************************** 330 C ** ** 331 C ** IMPORTANT VARIALES IN SUBROUTINE MATDIM 332 C ** * 333 C ** N: ORDER OF THE APPROXIMATING FUNCTIONS FOR * 334 C ** DISPLACEMENT FIELDS U AND V. ** 335 C ** ** 336 C ** NN: ORDER OF X IN THE APPROXIMATING FUNCTIONS. 337 C ** ** 338 C ** MM: ORDER OF Y IN THE APPROXIMATING FUNCTIONS. 339 C ** * 340 C KK: ORDER OF X IN THE VARIATION OF THE ** 341 C APPROXIMATING FUNCTIONS. ** 342 C ** ** 343 C * LL: ORDER OF Y IN THE VARIATION OF THE 344 C APPROXIMATING FUNCTIONS. 345 C ** * 346 C AMOD: NONDIMENSIONAL INTEGRALS. ** 347 C ** +* 348 C *********************** ************************** 349 C *********** ********************* ***************** 350 C 351 IMPLICIT REAL*8(A-H.O-Z) 352 COMMON /SUBI/ NN,MM,KK,LL 353 COMMON /SUB2/ XP,YP,XP1,YP1 354 COMMON /SUB3/ IGO 355 COMMON /SUB4/ AMOD 356 DIMENSION AMOD(7,1230),XPI(1230),YP1(1230),XP(730),YP(730) 357 C 358 C SUBROUTINE MATDIM CALCULATES THE NONDIMENSIONAL 359 C INTEGRATIONS OF THE VARIATIONAL EQUATION FOR U AND V. 360 C THE NONDIMENSIONAL INTEGRALS ARE FORMED IN THE FOLLOWING MANNER: 361 C STARTING WITH A(1,0) MULTIPLIED BY THE VARIATION OF A(1,0). THEN 362 C A(2,0) TIMES THE VARIATION OF A(1,0) TO A(N,O). THEN A(1,1) TO A(N-1,1) 363 C ARE MULTIPLIED BY THE VARIATION OF A(1,0). THE LAST INTEGRALS TO BE 364 C FORMED INVOLVE POWERS OF Y qNLYTIMES THE VARIATION OF A(1,0). THE 365 C ABOVE PROCESS IS REPEATED FOR THE VARIATION OF A(2,0) TO THE 366 C VARIATION OF A(O,N). THE ORDER OF THE VARIATIONS IS THE SAME AS 367 C THAT USED WHEN THE COEFFICIENTS ARE MULTIPLIED BY THE VARIATION 368 C OF A(1,0). THE NONDIMENSIONAL INTEGRALS OF B(NN,MM) TIMES 369 C THE VARITION'OF B(KK,LL) ARE FORMED IN THE SAME MANNER. 370 C 371 C NOTE: THE FIRST IN ( ) CORRESPONDS TO THE POWER OF X ASSOCIATED 372 C WITH EACH COEFFICIENT A AND B, AND THE SECOND NUMBER IN 373 C ( ) CORRESPONDS TO THE POWER OF Y ASSOCIATED WITH EACH 374 C COEFFICIENT A AND B. 375 C 376 C INITIALIZE NONDIMENSIONAL INTEGRALS TO ZERO. 377 DO 3 1=1,7 378 DO 2 IJ=1,730 379 2 AMOD(I.IJ)=O.OD+00 380 3 CONTINUE 381 J=O 382 JL=0 383 DO 20 LI=1,N 384 LL=LI-1 385 DO 30 KK=1,N a'

386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407.408 409 410 411 412 413 414 415 416 417 End of file IGO=O IF(KK+LL.GT. N) GO TO 20 DO 40 MI=1,N MM=M1-1 DO 50 NN=1,N IF(NN+MM.GT. N) GO TO 40 CALL CMOD(N,J,JL) 50 CONTINUE 40 CONTINUE NN=O DO 60 MM=1,N CALL CMOD(N,J,JL) 60 CONTINUE 30 CONTINUE 20 CONTINUE KK=O DO 70 LL=1,N IGO=O DO 80 M1=1,N MM=M1-1 DO 90 NN=1,N IF(NN+MM.GT. N) GO TO 80 CALL CMOD(N,J,JL) 90 CONTINUE 80 CONTINUE NN=O DO 95 MM=I,N CALL CMOD(N.J,JL) 95 CONTINUE 70 CONTINUE RETURN END oo

419 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 44,1 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 End of file I I C C C C C C C C C C C C C C C C C C C LU5KUUIINt LMUUMN,U,UL) CMOD IS USED AS A BRIDGE SUBROUTINE BETWEEN SUBROUTINES MATDIM AND AMOD1, AMOD2, AMOD3, AMOD4, AMODS, AMOD6, AND AMOD7. CMOD ALSO CALCULATES THE APPROPIATE POWERS OF X AND Y USED IN SUBROUTINE EXPAND TO DIMENSIONALIZE THE NONDIMENSIONAL INTEGRALS. **************** *************************************** ** ** ** IMPORTANT VARIABLES IN SUBROUTINE CMOD ** IGO: ** ** IGO: 0- DIAGONAL ELEMENT NOT ENCOUNTERED ** ** 1- DIAGONAL ELEMENT ENCOUNTERED ** ** ** AMOD: NONDIMENSIONAL INTEGRALS. ** * * ** IMPLICIT REAL*8(A-H.O-Z) COMMON /SUB1/ NN,MM,KK,LL COMMON /SUB2/ XP,YP,XP1,YP1 COMMON /SUB3/ IGO COMMON /SUB4/ AMOD DIMENSION XP(730),YP(730),XP1(1230),YPI(1230),AMOD(7,1230) C C AMOD(1,J),AMOD(2,d) IS USED TO CALCULATE K1 C AMOD(3,J), AMOD(4.J), AND AMOD(5,J) IS USED TO C CALCULATE K2, AMOD(6,J) AND AMOD(7,J) IS USED C TO CALCULATE K4, IN SUBROUTINE EXPAND. IF(IGO.EQ. 1) GO TO 10 IF(NN.EQ. KK.AND. MM.EQ. LL ) GO TO 10 GO TO 15 10 IGO=1 JL=JL+I XP(JL)=NN+KK YP(JL)=MM+LL CALL AMODI(dL) CALL AMOD2(JL) CALL AMOD6(JL) CALL AMOD7(N,JL) 15 J=J+1 XPl(d)=NN+KK YPI(d)=MM+LL CALL AMOD3(J) CALL AMOD4(J) CALL AMOD5(N,J) RETURN. END to

4b8 469 C 470 C 471 C 472 C 473 474 475 476 477 478 479 480 End of file bUBRUUIINE AMUU1()O AMOD1 CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1-NU**2)*(B/L) IN SUBMATRIX K1, WHICH IS CALCULATED IN SUBROUTINE EXPAND. IMPLICIT REAL*8(A-H,O-Z) COMMON /SUBI/ NN,MM.KK,LL COMMON /SUB4/ AMOD DIMENSION AMOD(7,1230) IF(KK+NN.EQ. 1) GO TO 20 AMOD(1,)=-DFLOAT(NN*KK)/DFLOAT((NN+KK-1)*(MM+LL+1)) 20 RETURN END 0

481 482 483 484 485 486 487 488 489 490 491 492 493 494 End of file SUBROUTINE AMOD2(J) C C AMOD2 CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED C BY (EH)/(1+NU)*(L/B) IN SUBMATRIX K1, WHICH IS CALCULATED IN C SUBROUTINE EXPAND. IMPLICIT REAL*8(A-H,O-Z) COMMON /SUB1/ NN,MM,KK,LL COMMON /SUB4/ AMOD DIMENSION AMOD(7,1230) IF(MM+LL.EQ. 1) GO TO 20 AMOD(2,J)=-DFLOAT(MM*LL)/DFLOAT((NN+KK+I)*(MM+LL-1)) AMOD(2,J)=AMOD(2.J)/(2.0D+00) 20 RETURN END H-I )I —

495 496 C 497 C 498 C 499 C 500 501 502 503 504 505 506 507 508 509 510 End of file SUBROUTINE AMOD3(J) AMOD3 CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1-NU) IN SUBMATRIX K2, WHICH IS CALCULATED IN SUBROUTINE EXPAND IMPLICIT REAL*8(A-H,O-Z) COMMON /SUB1/ NN,MM,KK,LL COMMON /SUB4/ AMOD DIMENSION AMOD(7.1230) IF(NN.EQ. O) GO TO 20 IF(MM.EQ. O) GO TO 20 D1=DFLOAT(NN)/DFLOAT(NN+KK) D2=DFLOAT(MM)/DFLOAT(MM4LL) AMOD(3,J)=(Dl*D2)/(2.0D+00) 20 RETURN END

511 512 513 514 515 516 517 518 519 520 521 522 523 524 End of file SUBROUTINE AMOD4(J) C C AMOD4 CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED C BY ((EH)NU)/(1-NU**2) IN SUBMATRIX K2, WHICH IS CALCULATED IN C SUBROUTINE EXPAND. IMPLICIT REAL*8(A-H,O-Z) COMMON /SUB1/ NN.MM,KK.LL COMMON /SUB4/ AMOD DIMENSION AMOD(7,1230) IF(NN+KK.EQ. O) GO TO 20 IF(MM.EO. O) GO TO 20 AMOD(4,.)=DFLOAT(MM)/DFLOAT(MM+LL) 20 RETURN END

525 526 C 527 C 528 C 529 C 530 531 532 533 534 535 536 537 538 539 540 541 542 543 End of file SUBROUTINE AMODS(N.J) AMOD5 CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1+NU) IN SUBMATRIX K2, WHICH IS CALCULATED IN SUBROUTINE EXPAND IMPLICIT REAL*8(A-H,O-Z) COMMON /SUB1/ NN,MM,KK,LL COMMON /SUB4/ AMOD DIMENSION AMOD(7.1230) IF(MM+LL.EO. 0) GO TO 20 IF(MM.EQ. O) GO TO 10 IF(NN.EQ. O) GO TO 20 AMOD(5,J)=(DFLOAT(NN)/DFLOAT(2*(NN+KK))) GO TO 20 10 D1=DFLOAT(NN)/DFLOAT(NN+KK) D2=FLAT(NDFLOAT(N)/DFLOAT(KK+N) AMOD(5,J)=(DI-D2)/(2.OD+00) 20 RETURN END

545 C 546 C 547 C 548 C 549 550 551 552 553 554 555 556 End of file ZUbKUUIINt AMUUtUl) AMOD6 CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1-NU**2)*(L/B) IN SUBMATRIX K4, WHICH IS CALCULATED IN SUBROUTINE EXPAND. IMPLICIT REAL*8(A-H,O-Z) COMMON /SUB1/ NN,MM,KK,LL COMMON /SUB4/ AMOD DIMENSION AMOD(7,1230) IF(MM+LL.EQ. 1) GO TO 20 AMOD(6,J)=-DFLOAT(MM+LL)/DFLOAT((NN+KK+I)*(MM+LL-i)) 20 RETURN END -1 U1

557 558 C 559 C 560 C 561 C 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 End of file SUBROUTINE AMOD7(N,J) AMOD7 CALCULATES THE NONDIMENSIONAL INTEGRALS THAT ARE MULTIPLIED BY (EH)/(1+NU)*(B/L) IN SUBMATRIX K4, WHICH IS CALCULATED IN SUBROUTINE EXPAND. IMPLICIT REAL*8(A-H,O-Z) COMMON /SUB1/ NN,MM,KK,LL COMMON /SUB4/ AMOD DIMENSION AMOD(7,1230) IF(LL.EQ. O) GO TO 10 IF(MM.E. 0) GO TO 15 IF(NN+KK.EQ. 1) GO TO 20 AMOD(7,d)=-DFLOAT(NN*KK)/DFLOAT((NN+KK-1)*(MM+LL+1)) AMOD(7,J)=AMOD(7, )/(2.OD+00) GO TO 20 10 IF(MM.EQ. O) GO TO 30 IF(NN+KK.EQ. 1) GO TO 20 D1=DFLOAT(NN*(NN-1))/DFLOAT((NN+KK-1)*(MM+1)) D2=DFLOAT(NN*(NN-1))/DFLOAT((NN+N-1)*(MM+I)) AMOD(7,J)=(D1-D2)/(2.0D+00) GO TO 20 15 D1=0.OD+00 IF(NN+KK.EQ. 1) GO TO 25 D1=DFLOAT(NN*(NN-1))/DFLOAT((NN+KK-1)*(LL+1)) 25 D2=DFLOAT(N*(N-1))/DFLOAT((KK+N-1)*(LL+1)) D3=DFLOAT(NN-N)/DFLOAT(LL+1) AMOD(7,d)=(D1-D2-D3)/(2.0D+00) GO TO 20 30 D1=O.OD+00 02=0.00+00 IF(NN+KK.EQ. 1) GO TO 40 D1=DFLOAT(NN*(NN-1))/DFLOAT(NN+KK-1) D2=DFLOAT(NN*(NN-1))/DFLOAT(NN+N-1) 40 D3=DFLOAT(N*(N-1))/DFLOAT(KK+N-1) D4=DFLOAT(N*(N-1))/DFLOAT(2*N-1) AMOD(7,d)=(DI-D2-D3+D4)/(2.OD+00) 20 RETURN END..j (n

502 C NAASA 2.1.001 DLUD FTN-A 10-29-75 THE UNIV OF MICH COMP CTR 503 SUBROUTINE DLUD (N,ADIM,A,TDIM,T,IV) 504 C 505 C COMPUTES THE LU-DECOMPOSITION OF THE N X N MATRIX A USING 506 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING. THIS FACTOR507 C IZATION MAY BE EXPRESSED IN THE FORM 508 C L(N-1)(N N-1)*....*L(1)*P(1)*A = U 509 C WHERE EACH L(J) IS THE IDENTITY MATRIX EXCEPT FOR THE SUB510 C DIAGONAL ELEMENTS IN COLUMN d, EACH P(d) IS A PERMUTATION 511 C MATRIX, AND U IS AN UPPER TRIANGULAR MATRIX. THIS IS THE 512 C PREPARATORY STEP IN SOLVING A SYSTEM OF LINEAR EQUATIONS, 513 C INVERTING A MATRIX, OR CALCULATING A DETERMINANT. A DISCUSSION 514 C OF GAUSSIAN ELIMINATION AND THE LU-DECOMPOSITION AND THEIR 515 C RELATIONSHIP TO THE NUMERICAL SOLUTION OF SYSTEMS OF LINEAR 516 C EQUATIONS MAY BE FOUND IN EITHER WILKINSON (1965,CHAPTER 4) 517 C OR FORSYTHE AND MOLER (1967). 518 C 519 INTEGER N,ADIM,TDIM,IV(1) 520 DOUBLE PRECISION A(ADIM,N),T(TDIM,N) 521 C 522 C N -> ORDER OF THE MATRIX A. 523 C ADIM -> ROW DIMENSION OF THE ARRAY A. BECAUSE A IS AN N X N 524 C MATRIX, ADIM SHOULD NOT BE LESS THAN N. IF ADIM IS LESS 525 C THAN N, THE CONTENTS OF A ARE IGNORED, AND THE MATRIX 526 C TO BE FACTORED IS ASSUMED TO BE STORED IN THE ARRAY T. 527 C SINCE ADIM MUST BE A POSITIVE INTEGER. IT IS RECOMMENDED 528 C THAT THE ACTUAL ARGUMENTS A AND T COINCIDE WHEN ADIM IS 529 C LESS THAN N TO AVOID THE INCONSISTENCY WHICH ARISES WHEN 530 C N EQUALS 1. 531 C A -> TWO-DIMENSIONAL ARRAY CONTAINING THE N X N MATRIX TO 532 C BE FACTORED, I.E., THE COEFFICIENT MATRIX OF THE SYSTEM 533 C OF LINEAR EQUATIONS OR THE MATRIX TO BE INVERTED. THE 534 C CONTENTS OF A ARE NOT ALTERED. 535 C TDIM -> ROW DIMENSION OF THE ARRAY T. 536 C T <- TWO-DIMENSIONAL ARRAY FOR RETURNING THE LU-DECOMPOSITION 537 C OF A. THE SUBDIAGONAL ELEMENTS OF THE d-TH COLUMN OF THE 538 C L(d) AND THE UPPER TRIANGULAR MATRIX U ARE RETURNED IN 539 C THE CORRESPONDING ELEMENTS OF T. IF ADIM IS LESS THAN N. 540 C T MUST CONTAIN THE MATRIX TO BE FACTORED WHEN THIS SUB541 C ROUTINE IS CALLED. 542 C IV <- VECTOR OF LENGTH N DEFINING THE PERMUTATION MATRICES 543 C P(J): MULTIPLICATION ON THE LEFT BY P(d) INTERCHANGES 544 C ROWS d AND IV(d). IF IV(d) IS NOT EQUAL TO d, THEN 545 C DET(A) = - DET(P(J)*A). AND TO AID IN THE COMPUTATION 546 C OF DET(A), IV(N) WILL CONTAIN +1 IF AN EVEN NUMBER OF 547 C INTERCHANGES ARE PERFORMED AND -1 IF AN ODD NUMBER. THUS 548 C DET(A) = IV(N)*T(1,1)*... *T(N,N). 549 C IV(N) WILL CONTAIN 0 IF A IS COMPUTATIONALLY SINGULAR. 550 C 551 INTEGER I,J,K,KP1,L 552 DOUBLE PRECISION TMP 553 DOUBLE PRECISION PIV 554 C 555 C IF ADIM IS GREATER THAN OR EQUAL TO N. THE CONTENTS OF A ARE MOVED 556 C TO T; WHILE IF ADIM IS LESS THAN N, THIS INITIAL DATA MOVEMENT IS 557 C SKIPPED. SINCE THE DATA MOVEMENT TIME IS PROPORTIONAL TO N**2 AND 558 C THE COMPUTATION TIME PROPORTIONAL TO N**3/3, SIGNIFICANT SAVINGS 559 C SHOULD NOT BE EXPECTED. 560 C 561 IF ( ADIM.LT. N ) GO TO 8110

562 DO 8100 J = 1, N 563 DO 8100 I = 1, N 564 8100 T(I,) = A(I,J) 565 8110 CONTINUE 566 C 567 C GAUSSIAN ELIMINATION CONSISTS OF N-I STAGES. DURING THE K-TH 568 C STAGE, THE PERMUTATION MATRIX P(K), THE LOWER TRIANGULAR MATRIX 569 C L(K), AND THE K-TH ROW OF U ARE COMPUTED BASED ON THE CURRENT 570 C ELEMENTS OF THE K-TH RESIDUAL MATRIX, I.E., THE ELEMENTS T(I.,), 571 C I,J=K...N. ONLY THE ELEMENTS OF THE K-TH RESIDUAL MATRIX ARE 572 C REFERENCED DURING THE K-TH STAGE OF GAUSSIAN ELIMINATION. 573 C 574 IV(N) = 1 575 DO 8260 K = 1, N 576 PIV = DABS(T(K,K)) 577 IF ( K.GE. N ) GO TO 8260 578 C 579 C SELECT THE PIVOT ROW FOR THE K-TH STAGE BY PARTIAL PIVOTING, 580 C I.E., THE MAXIMUM ELEMENT IN THE 1-ST COLUMN OF THE K-TH RESIDUAL 581 C MATRIX, AND SET IV(K) ACCORDINGLY. THE VARIABLE L HOLDS THE 582 C SUBSCRIPT OF THE PIVOT ROW, AND PIV THE ABSOLUTE VALUE OF THE 583 C PIVOT ELEMENT. 584 C 585 L = K 586 KPI = K + 1 587 DO 8200 I = KP1, N 588 IF ( PIV.GE. DABS(T(I,K)) ) GO TO 8200 589 PIV = DABS(T(I,K)) 590 L = I 591 8200 CONTINUE 592 IV(K) = L 593 C 594 C SAVE THE PIVOT ELEMENT IN TMP. IF P(K) IS NONTRIVIAL, I.E., IV(K) 595 C IS NOT EQUAL TO K, THE PIVOT ELEMENT IS ALWAYS NONZERO; OTHERWISE, 596 C THE PIVOT ELEMENT MUST BE CHECKED. IF THE PIVOT IS ZERO, I.E., 597 C THE MATRIX IS COMPUTATIONALLY SINGULAR, THEN T(I,K) IS ZERO FOR 598 C I=K...N, AND THE COMPUTATION MAY PROCEED TO THE NEXT STAGE. 599 C 600 TMP = T(L,K) 601 IF ( K.NE. L ) GO TO 8210 602 IF ( PIV.GT. O.DO ) GO TO 8220 603 IV(N) = 0 604 GO TO 8260 605 8210 CONTINUE 606 IV(N) = -IV(N) 607 T(L,K) = T(K,K) 608 T(K,K) = TMP 609 8220 CONTINUE 610 C 611 C COMPUTE THE NONTRIVIAL ELEMENTS OF L(K). BECAUSE OF THE PARTIAL 612 C PIVOTAL STRATEGY, THE ABSOLUTE VALUE OF L(I,K)=-T(I,K)/T(K,K) IS 613 C LESS THAN OR EQUAL TO 1. IF T(K) DENOTES THE K-TH RESIDUAL MATRIX, 614 C THEN THE SUBDIAGONAL ELEMENTS OF THE 1-ST COLUMN OF THE MATRIX 615 C L(K) P(K) T(K) ARE ZERO. THESE ELEMENTS ARE NOT ACTUALLY CALCU616 C LATED, THEY ARE REPLACED BY THE ELEMENTS OF L(K). 617 C 618 TMP = -TMP 619 DO 8230 I = KP1, N 620 8230 T(I,K) = T(I,K) / TMP 00 co

622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 End of ftll C APPLY P(K) AND L[K) TO THE K-1H RESIDUAL MAIRIX CULUMNWIbt, l.t., C FOR d=K+1...N, INTERCHANGE T(K,J) AND T(L,J), THE (K,d)-ELEMENT C OF U, AND THEN FOR I=K+1...N, REPLACE T(I,i) BY C T(I.J) + L(K)(I,K) * T(K,J). C ( DO 8250 d = KPI, N TMP = T(L,J) T(L,J) = T(K,J) T(K,d) = TMP DO 8240 I = KP1, N T(I,J) = T(I,J) + T(I.K) * TMP CONTINUE CONTINUE IF (PIV.EQ. 0.00) IV(N) = O RETURN 8240 8250 8260 C C FORSYTHEG.E. AND MOLERC.B. 1967. COMPUTER SOLUTION OF LINEAR C ALGEBRAIC SYSTEMS. ENGLEWOOD CLIFFS,N.J.: PRENTICE-HALL. C WILKINSON,J.H. 1965. THE ALGEBRAIC EIGENVALUE PROBLEM. OXFORD: C CLARENDON PRESS. C THE UNIVERSITY OF MICHIGAN COMPUTING CEN1 TER ( a C NUMERICAL ANALYSIS LIBRARY - JULY 1975 END

b4) NAAbA Z.I.UUJ Ubb I-IN —A 1U-y-/I IHt UNIV Ut Ml1n" UMv LIK 646 SUBROUTINE DBS (N.TDIM,T,IV,B) 647 C 648 C SOLVES THE SYSTEM OF LINEAR EQUATIONS AX=B, WHERE A DENOTES 649 C THE N X N COEFFICIENT MATRIX AND X AND B ARE N-VECTORS, BY 650 C BACK-SUBSTITUTION IN THE LU-DECOMPOSITION OF A. THIS 651 C DECOMPOSITION MUST BE PROVIDED IN THE ARRAY T AND VECTOR IV 652 C VIS-A-VIS THE SUBROUTINE DLUD. THE LU-DECOMPOSITION MAY BE 653 C EXPRESSED IN THE FORM 654 C L(N-1)*P(N-1)*...*L(1)*P()*A = U, 655 C WHERE EACH L(d) IS THE IDENTITY MATRIX EXCEPT FOR THE SUB656 C DIAGONAL ELEMENTS IN COLUMN,J EACH P(J) IS A PERMUTATION 657 C MATRIX, AND U IS AN UPPER TRIANGULAR MATRIX. USING THIS 658 C NOTATION, THE BACK-SUBSTITUTION CONSISTS OF FORMING 659 C Y = L(N-1)*P(N-1)*...*L(1)*P(1)*B 660 C AND SOLVING THE UPPER TRIANGULAR SYSTEM OF LINEAR EQUATIONS 661 C UX = Y. I.E., FOR I=N...1 662 C X(I)=(Y(I)-U(I,N)*X(N)-... -U(II+1)*X(I+1))/U(I,I). 663 C THIS BACK-SUBSTITUTION YIELDS A VECTOR X WHICH IS THE EXACT 664 C SOLUTION OF A SYSTEM OF LINEAR EQUATIONS (A+E)X=B, WHERE 665 C //E// IS GENERALLY ON THE ORDER OF N*//A//*MACHEPS. THIS METHOD 666 C OF SOLVING SYSTEMS OF LINEAR EQUATIONS IS DESCRIBED IN BOTH 667 C WILKINSON (1965,CHAPTER 4) AND FORSYTHE AND MOLER (1967). 668 C 669 INTEGER N,TDIM,IV(1) 670 DOUBLE PRECISION T(TDIM,N),B(1) 671 C 672 C N -> ORDER OF THE SYSTEM OF LINEAR EQUATIONS. 673 C TDIM -> ROW DIMENSION OF THE ARRAY T. 674 C T -> TWO-DIMENSIONAL ARRAY CONTAINING THE LU-DECOMPOSITION 675 C OF THE COEFFICIENT MATRIX VIS-A-VIS THE SUBROUTINE DLUD. 676 C IV -> VECTOR CONTAINING THE INTERCHANGE INFORMATION GENERATED 677 C BY THE SUBROUTINE DLUD DURING THE COMPUTATION OF THE 678 C LU-DECOMPOSITION OF THE COEFFICIENT MATRIX. 679 C B -- VECTOR CONTAINING THE RIGHT-HAND SIDE OF THE SYSTEM OF 680 C LINEAR EQUATIONS. THE CONTENTS OF B ARE REPLACED BY THE 681 C ELEMENTS OF THE SOLUTION. 682 C 683 INTEGER I,K,KP1,L 684 DOUBLE PRECISION TMP 685 C 686 C REPLACE THE CONTENTS OF THE VECTOR B BY THE VECTOR 687 C (L(N-1) (P(N-1)... ( L(1) (P() B)...). 688 C THIS COMPUTATION IS PERFORMED IN N-1 STAGES, WHERE DURING THE 689 C K-TH STAGE, THE CONTENTS OF B ARE REPLACED BY L(K)( P(K) B ). 690 C THE PERMUTATION MATRICES P(K) ARE DEFINED BY THE VECTOR IV, I.E., 691 C MULTIPLICATION BY P(K) SIMPLY INTERCHANGES THE K-TH AND IV(K)-TH 692 C ELEMENTS OF THE VECTOR. THE LOWER TRIANGULAR MATRIX L(K) IS THE 693 C IDENTITY MATRIX EXCEPT FOR THE SUBDIAGONAL ELEMENTS OF THE K-TH 694 C COLUMN, WHICH ARE STORED IN THE CORRESPONDING ELEMENTS OF THE 695 C ARRAY T. 696 C 697 DO 8110 K = 1, N 698 IF ( K.GE. N ) GO TO 8110 699 L = IV(K) 700 TMP = B(L) 701 B(L) = B(K) 702 B(K) = TMP 703 KP1 = K + i 704 DO 8100 I = KPi, N co 0

705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 ( 8100 B(I) = B(I) + T(I,K) * TMP 8110 CONTINUE C REPLACE THE CONTENTS OF THE VECTOR B BY THE SOLUTION TO THE SYSTEM C OF LINEAR EQUATIONS WITH UPPER TRIANGULAR COEFFICIENT MATRIX U C AND RIGHT-HAND SIDE VECTOR B. THE USUAL FORMULAS FOR THE BACKC SUBSTITUTION, WHICH ARE BASED ON THE SUCCESSIVE ROWS OF THE MATRIX C AND ARE SUITABLE WHEN INNER-PRODUCTS ARE ACCUMULATED, ARE NOT C EMPLOYED. THE COMPUTATION HAS INSTEAD BEEN ARRANGED TO REFERENCE C THE SUCCESSIVE COLUMNS OF U. THUS AFTER B(I) HAS BEEN COMPUTED, C IT IS REMOVED FROM THE SYSTEM BY SUBTRACTING B(I) TIMES THE I-TH C COLUMN OF U FROM THE RESIDUAL VECTOR B(1)...B(I-1). K = N 8200 B(K) = B(K) / T(K,K) IF ( K.LE. 1 ) RETURN TMP = -B(K) KP1 = K K K - DO 8210 I = 1, K 8210 B(I) = B(I) + T(I,KPI) * TMP GO TO 8200 FORSYTHE,G.E. AND MOLER,C.B. 1967. COMPUTER SOLUTION OF LINEAR ALGEBRAIC SYSTEMS. ENGLEWOOD CLIFFS,N.d.: PRENTICE-HALL. WILKINSON,d.H. 1965. THE ALGEBRAIC EIGENVALUE PROBLEM. OXFORD: CLARENDON PRESS. 727 ( 728 ( 729 ( 730 ( 731 C 732 C 733 C 734 End of file c c c c c C *~ THE UNIVERSITY OF MICHIGAN COMPUTING CENTER NUMERICAL ANALYSIS LIBRARY - JULY 1975 END 00 H —

595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 6'18 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 C C C C C C C C C C C C C C C C C C C C C C C SUBROUTINE SCOFF(EL,B,NMSIZE,COFF.EUNAX,ITITL) SCOFF CALCULATES THE STRESS FIELD COEFFICIENTS FOR SIGMA-X, AX(d,1), SIGMA-Y, AX(d,2), AND TAU-XY,AX(d,3) USING THE DISPLACEMENT FIELD COEFFICIENTS FOUND IN SUBROUTINE EXPAND AND HOOKE'S LAW. IMPLICIT REAL+8(A-H,O-Z) DIMENSION COFF(75),AX(35,3),DUX(35),DUY(35),DVX(35),DVY(35) DIMENSION ITITL(50) d=1. THIS SET OF DO LOOPS CALCULATES THE DERIVATIVES OF THE THE APPROXIMATING FUNCTIONS FOR THE DISPLACEMENT FIELDS U AND V WITH RESPECT TO X AND WITH RESPECT TO Y. DO 10 MI=1,N MM=MI-1.DO 20 NN=1,N IF(NN+MM.GT. N) GO TO 10 J1=MSIZE+d DUX IS THE DERIVATIVE OF THE A FIELD, U, WITH RESPECT TO X. DUY IS THE DERIVATIVE OF THE A FIELD, U, WITH RESPECT TO Y. DVX IS THE DERIVATIVE OF THE Al FIELD, V, WITH RESPECT TO X. DVY IS THE DERIVATIVE OF THE Al FIELD, V, WITH RESPECT TO Y. DUX(d)=DFLOAT(NN)*COFF(J) DUY(d)=DFLOAT(MM)*COFF(d) DVX(J)=DFLOAT(NN)*COFF(di) DVY(J)=DFLOAT(MM)*COFF(J1) 20 J=J+1 10 CONTINUE NN=O DO 30 MM=1,N d1=MSIZE+d DUX(d)=DFLOAT(NN)*COFF(d) DUY(J)=DFLOAT(MM)*COFF(d) DVX(d)=DFLOAT(NN)*COFF(d1) DVY(J)=DFLOAT(MM)*COFF(dl) 30 J=J+1 J1=MSIZE-N+1 E1=E/((1.0O+00)-UN**2) G=E/((2.0D+00)*((1.OD+00)+UN)) THESE STATEMENTS CALCULATE THE FIELD COEFFICIENTS. AX(,11)=E1*(DUX(1)+UN*DVY(d1)) AX(1,2)=E1*(UN*DUX(1)+DVY(d1)) AX(1,3)=G*(DUY(d1)+DVX(1)) NN=2 NI=N IF=O J=2 PPROXIMATING FUNCTION FOR DISPLACEMENT PPROXIMATING FUNCTION FOR DISPLACEMENT PPROXIMATING FUNCTION FOR DISPLACEMENT PPROXIMATING FUNCTION FOR DISPLACEMENT 00 fl) CONStANT TERMS FOR THE STRESS THIS PART OF SCOFF CALCULATES THE VARIOUS TERMS OF THE STRESS FIELDS THAT ARE ASSOCIATED WITH VARIOUS POWERS OF X AND Y. 40 IF(NN.GT. N1) GO TO 50 d1=d 2=N+-J-1

655 IF(IF.GE. 1) J1=d1+IF 656 AX(J,1)=E1*(DUX(J1)+UN*DVY(J2)) 657 AX(d,2)=E1*(UN*DUX(dJ)+DVY(d2)) 658 AX(d,3)=G*(DUY(J2)+DVX(J1)) 659 J=J+1 660 NN=NN+1 661 GO TO 40 662 50 NN=2 663 IF=IF+l 664 NI=N1-1 665 IF(NN.GT. N1) GO TO 60 666 GO TO 40 667 60 J1=N 668 IF=1 669 IF1=1 670 65 J1=J1+IF 671 IF(J1.GE. MSIZE-N+1) GO TO 70 672 J2=N+J 673 AX(d,1)=EI*(DUX(J1)+UN*DVY(d2)) 674 AX(,J2)=E1*(UN*DUX(J1)+DVY(d2)) 675 AX(d,3)=G*(DUY(J2)+DVX(d1)) 676 J=J+1 677 IF=N-IF1 678 IF1=IF1+1 679 GO TO 65 680 70 J=d-1 681 C 682 C ZEROS SETS STRESS FIELD COEFFICIENTS THAT ARE SMALL IN RELATIVE 683 C MAGNITUDE EQUAL TO ZERO. 684 CALL ZEROS(EL,B,NAX) 685 C 686 C THIS PART OF SCOFF PRINTS THE STRESS COEFFICIENTS TO THE OUTFILE. 687 C dd IS THE POWER OF X ASSOCIATED WITH EACH STRESS FIELD COEFFICIENT 688 C AND JM IS THE POWER OF Y ASSOCIATED WITH EACH STRESS FIELD 689 C COEFFICIENT. 690 WRITE(8,500) ITITL 691 WRITE(8,100) 692 DO 80 I=1,3 693 IF(I.NE. 1) WRITE(8,500) ITITL 694 IF(I.EO. 2) WRITE(8,101) 695 IF(I.EQ. 3) WRITE(8,102) 696 100 FORMAT(1H-,3X,'COEFFICIENTS OF NORMAL STRESS IN X-DIRECTION' 697 $) 69,8 101 FORMAT(1H-,3X,'COEFFICIENTS OF NORMAL STRESS IN Y-DIRECTION' 699 $) 700 102 FORMAT(1H-,3X,'COEFFICIENTS OF SHEAR STRESS') 701 500 FORMAT(1H1,2(/,IH-),40X,50A 1) 702 J=1 703 JM=0 704 NI=N-1 705 DO 90 J1=1,N 706 JJ=J1-1 707 WRITE(8,103) dJJ,M,AX(d,I) 708 90 d=J+1 709 DO 91 JM=1,N1 710 DO 92 JJ=1,N1 711 IF(dJ+dM.GT. N1) GO TO 91 712 WRITE(8,103) JJ,JM,AX(J,I) 713 92 J=J+1 00 L'J

I 15 716 717 718 719 720 721 722 End of file JdJ= DO 93 JM=1,N1 WRITE(8,103) JJ,JM,AX(J,I) 93 J=J+l 80 CONTINUE 103 FORMAT(6X.'AX(',I1, ',.I1,')=',D12.5) RETURN END 00 U:

' 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 End of file SUBROUTINE ZEROS(EL,B,N,AX) IMPLICIT REAL*8(A-H,O-Z) DIMENSION AXC(35,3),AX(35,3) NI=N-1 DO 10 K=1,3 d=1 DO 15 II=1,N I-=II-1 AXC(J.K)=AX(d,K)*EL**I 15 J=d+1 DO 20 L=1,N1 DO 25 I=1,N1 IF(I+L.GT. N1) GO TO 20 AXC(d,K)=AX(d,K)*EL**I*B**L 25 d=J+1 20 CONTINUE DO 30 L=1,N1 AXC(J,K)=AX(J,K)*B**L 30 J=J+1 J=J-1 AXMAX=DABS(AXC(1,K)) DO 35 I=2,J AC=DABS(AXC(I,K)) IF(AC.GT. AXMAX) AXMAX=AC 35 CONTINUE DO 40 I=1,d AC=DABS(AXC(I,K)) AB=AC/AXMAX IF(AB.LT. 1.OD-06) AX(I,K)=O.OD+00 40 CONTINUE 10 CONTINUE RETURN END 00 Vn

4I aoo UOKUU 1I 11C UUU I t1,|ILC,ll, DbUrr, I I I I L 757 C 758 C DOUT PRODUCES TABLES OF THE DISPLACEMENT FIELDS U AND V AT VARIOUS 759 C POINTS ON THE BAND SAW BLADE USING THE DISPLACEMENT FIELD COEFFICIENTS 760 C FOUND IN SUBROUTINE EXPAND. 761 IMPLICIT REAL*8(A-H,O-Z) 762 DIMENSION DISV(6,11),COFF(75),DISU(6. 11),XX111),Y1(6) 763 DIMENSION ITITL(50) 764 Y2=O.OD+00 765 DELX=EL/(1.OD+01) 766 DELY=B/(5.0D+00) 767 DO 20 I=1,6 768 X2=O.OD+00 769 DO 30 d=1,11 770 J1=1 771 d2=1+MSIZE 772 DISU(I,J)=O.OD+00 773 DISV(I,J)=O.OD+00 774 DO 40 K=1,N 775 DISU(I,J)=DISU(I.J)+COFF(dJ1)X2**K 776 DISV(I,J)=DISV(I,d)+COFF(d2)*X2**K 777 d2=J2+1 778 40 d1=J1+1 779 DO 50 L=1,N 780 DO 60 K=1,N 781 IF(K+L.GT. N) GO TO 50 782 DISU(I.J)=DISU(I,J)+COFF(d1)*X2**K*Y2**L 783 DISV(I,d)=DISV(I,J)+COFF(d2)*X2**K*Y2**L 784 J2=d2+1 785 60 d1=J1+1 786 50 CONTINUE 787 DO 70 L=1,N 788 DISU(I,d)=DISU(I,J)+COFF(dt)*Y2**L 789 DISV(I,d)=DISV(Id.)+COFF(d2)*Y2**L 790 J2=02+-1 791 70 01=d1+1 792 XI(d)=X2 793 30 X2=X2+DELX 794 Y1(I)=Y2 795 20 Y2=Y2+DELY 796 WRITE(8,500) ITITL 797 WRITE(8,100) 798 100 FORMAT(1H-,56X,'DISPLACEMENT FIELD U') 799 DO 90 IJ=1,2 800 IF(IJ.EQ. 2) WRITE(8,500) ITITL 801 IF(Id.EQ. 2) WRITE(8,120) 802 120 FORMAT(1H-,56X,'DISPLACEMENT FIELD V') 803 WRITE(8,107) 804 107 FORMAT(1H,124(' ')) 805 101 FORMAT(12X.'I'.50X,'X-COORDINATE',T125,'I') 806 102 FORMAT(12X,'I',T125,'I',/,12X,'I',11(1X,D9.3),T125,'I') 807 108 FORMAT(IX,124('-')) 808 109 FORMAT(IH.124(' ')) 809 500 FORMAT(IH1,2(/,1H-),40X,50A1) 810 WRITE(8,105) 811 105 FORMAT(5X,'Y',6X,'I',T125,'I') 812 WRITE(8,106) 813 106 FORMAT(1X,'COORDINATE',iX,'I',T125,'I') 814 DO 80 J=1,6 815 dJ17-d 00 o~

816 817 818 819 820 821 822 823 824 825 End of file IF(Id.EQ. 1) WRITE(8,104) Y1(d1).(DISU(d1.I),I=1,11) 80 IF(Id.EQ. 2) WRITE(8.104) Y1(d1),(DISV(d1,I),I=1,11) WRITE(8,108) WRITE(8,102) (X1(I).I=1,11) WRITE(8,101) WRITE(8,109) 104 FORMAT(12X,'I'.T125,'I',/,1X,D9.3.2X.''',11(1X,D9.3),T125,'I') 90 CONTINUE RETURN END 00 -,,!

826 SUBROUTINE STROUT(N.EL.BAX,ITITL) 827 C 828 C STROUT PRODUCES TABLES OF THE STRESS FIELDS SIGMA-X, SIGMA-Y, 829 C TAU-XY, AND THE MAXIMUM NPRINCIPAL NORMAL STRESS. 830 IMPLICIT REAL*8(A-H,O-Z) 831 DIMENSION AX(35,3),SIGX(6,11),SIGY(6,11),TAUXY(6,11) 832 $,SIGMAX(6,11),X1(1l),Y(6),ITITL(50) 833 N1=N-1 834 DELX=EL/(1.0D+01) 835 DELY=B/(5.0D+00) 836 Y2=0.OD+00 837 DO 20 1=1,6 838 X2=O.OD+00 839 DO 30 J=1,11 840 d1=2 841 SIGX(I.J)=AX(1,1) 842 SIGY(I,J)=AX(1,2) 843 TAUXY(I,d)=AX(1.3) 844 DO 40 K=1,NI 845 SIGX(I,J)=SIGX(I,J)+AX(d1,1)*X2**K 846 SIGY(I,J)=SIGY(I,J)+AX(J1,2)*X2**K 847 TAUXY(I,J)=TAUXY(I,d)+AX(dl.3)*X2**K 848 40 J1=J1+1 849 DO 50 L=1,NI 850 DO 60 K=1,N1 851 IF(L+K.GT. NI) GO TO 50 852 SIGX(I,d)=SIGX(I,J)+AX(J1,1)*X2**K*Y2**L 853 SIGY(I,J)=SIGY(I,J)+AX(d1,2)*X2**K*Y2**L 854 TAUXY(I,d)=TAUXY(I,d)+AX(d1,3)*X2**K*Y2**L 855 60 J1=J1+1 856 50 CONTINUE 857 DO 70 L=1,N1 858 SIGX(I,J)=SIGX(I,d)+AX(Jd,1)*Y2**L 859 SIGY(I,J)=SIGY(I,d)+AX(dJ,2)*Y2**L 860 TAUXY(I,J)=TAUXY(I,.J)+AX(d,3)*Y2**L 861 70 J1=J1+1 862 X1(J)=X2 863 CEN=(SIGX(I,d)+SIGY(I,d))/(2.00+00) 864 SRAD=((SIGX(I,d)-SIGY(I,J))/(2.OD+00))**2+TAUXY(I,d)**2 865 RAD=DSQRT(SRAD) 866 SIGMAX(I,d)=CEN+RAD 867 30 X2=X2+DELX 868 Y1(I)=Y2 869 20 Y2=Y2+DELY 870 DO 90 IJ=1,4 871 WRITE(8,500) ITITL 872 500 FORMAT(1HI,2(/,1H-),40X,50A1) 873 IF(Id.EQ. 1) WRITE(8,100) 874 IF(IJ.EQ. 2) WRITE(8,108) 875 IF(IJ.EQ. 3) WRITE(8,109) 876 IF(IJ.EQ. 4) WRITE(8,110) 877 100 FORMAT(1H-,64X,'SIGMA X') 878 108 FORMAT(1H-,64X,'SIGMA Y') 879 109 FORMAT(1H-,61X,'SHEAR STRESS') 880 110 FORMAT(1H-,48X,'MAXIMUM PRINCIPAL NORMAL STRESS') 881 WRITE(8,107) 882 107 FORMAT(IH,124(' ')) 883 101 FORMAT(12X,'I',50X,'X-COORDINATE',T125,'I') 884 102 FORMAT(12X,'I',T125,'I'./,12X,'I',I1(IX,D9.3),T125,'I') 885 103 FORMAT(IX,124('-')) 00 00

886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 End of file WRITE(8,104) 104 FORMAT(5X,'Y',6X.'I',T125,'I') WRITE(8,105) 105 FORMAT(1X,'COORDINATE',IXX,'I',T125,'I') DO 80 J=1,6 J1=7-J IF(IJ.EQ. 2) WRITE(8,106) Y1(J1),(SIGY(d1,I),I=,111) IF(Id.EQ. 3) WRITE(8,106) Y1(J1),(TAUXY(Jl,I),I=l,11) 80 IF(IJ.EQ. 4) WRITE(8.106) Y1(J1),(SIGMAX(d1.I),I=1,11) WRITE(8,103) WRITE(8,102) (X1(I),I=1,11) WRITE(8,101) WRITE(8,107) 106 FORMAT(12X,'I' T125,'I',/,IX,D9.3,2X,'I',11(IX,D9.3),T125.'I') 90 CONTINUE RETURN END 00 DC

904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 End of f le SUBROUTINE TDRIVE(NT,TX,KOUT) C C SUBROUTINE TDRIVE IS USED TO LET THE USER ENTER C X, Y COORDINATES AND TEMPERATURE DATA FOR VARIOUS C POINTS ON THE BANDSAW BLADE. IMPLICIT REAL*8(A-H,O-Z) DIMENSION TX(35),XT(70),YT(70).TXI(70) 10=5 WRITE(6,102) 102 FORMAT(/,9X,'BEGIN ENTRY OF TEMPERATURE DATA') WRITE(6,106) 106 FORMAT(/,I1X,'IS THE DATA TO READ FROM A DATA FILE',/, $11X,'0-NO',/,llX.'1-YES') READ(5.107) IR 107 FORMAT(I1) IF(IR.EQ. 1) IO=11 IF(IO.EQ. 11) REWIND 11 WRITE(6,104) 104 FORMAT(/,11X,'ENTER THE NUMBER OF TEMPERATURE DATA POINTS', $' (12 FORMAT)') READ(5,105) NPTS 105 FORMAT(I2) 103 FORMAT(/,12X,'ENTER X-COORDINATE.Y-COORDINATE,AND TEMPERATURE', $' FOR POINT',I3, ' ON THE BANDSAW BLADE') DO 60 I1=1,NPTS IF(IO.EQ. 5) WRITE(6,103) I1 60 READ(I0,201) XT(Ii),YT(II),TXI(I1) ( C C SUBROUTINE TCOFF CALCULATES THE TEMPERATURE FIELD C COEFFICIENTS USING LINER REGRESSION. CALL TCOFF(NT,NPTS,XT,YT,TXI,TX,KOUT) 201 FORMAT(3D15.6) RETURN END 0! b I

938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 98'1 982 983 984 985 986 987 988 989 End of file l SUBROUTINE TCOFF(NT.NPTS,XT,YT,TXI,TX,KOUT) C C SUBROUTINE TCOFF CALCULATES THE TEMPERATURE FIELD C COEFFICIENTS GIVEN TEMPERATURES AT VARIOUS POINTS C ON THE BANDSAW BLADE. THE TEMPERATURE DATA IS READ C IN SUBROUTINE TDRIVE. THE TEMPERATURE COEFFICIENTS C ARE USED IN SUBROUTINE TEMP TO CALCULATE THE EFFECTS C OF A TEMPERATURE FIELD FOR THE RIGHT HAND SIDE VECTOR C OF THE SYSTEM OF EQUATIONS. IMPLICIT REAL*8(A-H,O-Z) DIMENSION IV(35),TX(35),A(70,35),T(35,35) DIMENSION ATA(35.35).TXI(70),XT(70),YT(70) NM=35 MTZ=2*NT+1 DO 10 I=1.NT 10 MTZ=MTZ+(NT-I) 1 C C C C C c C THIS PART OF TCOFF CALCULATES THE COEFFICIENT MATRIX USED BY SUBROUTINES DLUD AND DBS. SUBROUTINE DLUD DECOMPOSES THE MATRIX(A) INTO LOWER AND UPPER TRIANGLUR MATRICES. SUBROUTINE DBS CALCULATES THE TEMPERATURE FIELD COEFFICIENTS USING THE DECOMPOSITION OF THE MATRIX(A) AND THE TEMPERATURE DATA TX. DO 45 I=1.NPTS 0=2 A(I,1)=1.OD+00 DO 25 K=1,NT A(I,J)=XT(I)**K 25 J=J+1 DO 30 L=1,NT DO 35 K=I,NT IF(K+L.GT. NT) GO TO 30 A(I,J)=XT(I)**K*YT(I)**L 35 J=J+1 30 CONTINUE DO 40 L=1,NT A(I,J)=YT(I)**L 40 J=d+1 45 CONTINUE SUBROUTINE ATA FORMS THE NORM MATRIX (ATRANSPOSE) TIMES A AND (ATRANSPOSE) TIMES TXI CALL ATAM(MTZ,NPTS,A,TXI,TX,ATA) CALL DLUD(MTZ,NM,ATA,NM,T,IV) IF(IV(MTZ).EQ. O) GO TO 75 CALL DBS(MTZ,NM,T,IV,TX) GO TO 80 75 WRITE(6,102) 102 FORMAT(4X,'***ERROR** MATRIX IS SINGULAR IN SUBROUTINE TCOFF') KOUT=1 80 RETURN END Ho H C C C

990 991 992 993 C 994 C 995 996 997 998 999 1000 1001 1002 1003 1004 1005 End of file SUBROUTINE ATAM(MTZ,NPTS,A,TXI,TX,ATA) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(70,35).TXI(70),ATA(35,35),TX(35) CALCULATE (ATRANSPOSE) TIMES A AND (ATRANSPOSE) TIMES TXI DO 25 I=1,MTZ DO 25 J=1,MTZ SUMA=O.OD+00 SUMB=O.OD+00 DO 30 K=1,NPTS SUMA=SUMA+A(K,J)*A(K,I) 30 SUMB=SUMB+TXI(K)*A(K,I) ATA(I, )=SUMA 25 TX(I)=SUMB RETURN END I',,J

1006 SUBROUTINE TOUT(N,EL,BTFLITITL) 1007 C 1008 C TOUT PRODUCES A TABLE OF TEMPERATURES AT VARIOUS POINTS ON THE 1009 C BAND SAW BLADE USING THE TEMPERATURE FIELD COEFFICIENTS FOUND IN 1010 C SUBROUTINE TCOFF. 1011 IMPLICIT REAL*8(A-H,O-Z) 1012 DIMENSION ITITL(50).TFL(35),TOU(6,11) X1(11).Y1(6) 1013 WRITE(8,500) ITITL 1014 WRITE(8.510) 1015 510 FORMAT(1H-,4X,'TEMPERATURE FIELD COEFFICIENTS') 1016 11=1 1017 N2=0 1018 M2=0 1019 WRITE(8,520) N2,M2,TFL(I1) 1020 11=2 1021 DO 55 I=1,N 1022 WRITE(8,520) I.M2,TFL(I1) 1023 55 11=I1+1 1024 DO 65 d=1.N 1025 DO 75 I=1,N 1026 IF(I+J.GT. N) GO TO 65 1027 WRITE(8,520) I, J TFL(I1) 1028 75 I1=I1+1 1029 65 CONTINUE 1030 DO 85 J=1,N 1031 WRITE(8,520) N2,d,TFL(I1) 1032 85 I1=11+1 1033 520 FORMAT(8X,'T(',I1,',',I1,')=',D12.5) 1034 Y2=O.OD+00 1035 DELX=EL/(1.0D+01) 1036 DELY=B/(5.0D+00) 1037 DO 20 1=1,6 1038 X2=0.OD+00 1039 DO 30 d=1,11 1040 d1=2 1041 TOU(I.J)=TFL(1) 1042 DO 40 K=1,N 1043 TOU(I,d)=TOU(I,.)+TFL(dJ)*X2**K 1044 40 d1=d1+1 1045 DO 50 L=1,N 1046 DO 60 K=1,N 1047 IF(K+L.GT. N) GO TO 50 1048 TOU(I,J)=TOU(I,d)+TFL(d1)*X2**K*Y2**L 1049 60 dJ1=1+1 1050 50 CONTINUE 1051 D0 70 L=1,N 1052 TOU(I,d)=TOU(I.,)+TFL(d1)*Y2**L 1053 70 d1=d1+1 1054 X1(d)=X2 1055 30 X2=X2+DELX 1056 YI(I)=Y2 1057 20 Y2=Y2+DELY 1058 WRITE(8,500) ITITL 1059 500 FORMAT(1H1,2(/,1H-),56X,50A1) 1060 WRITE(8,100) 1061 100 FORMAT(1H-,56X,'TEMPERATURE DISTRIBUTION') 1062 WRITE(8,107) 1063 107 FORMAT(IH,124(' ')) 1064 101 FORMAT(12X,'I',50X,'X-COORDINATE',T125,'I') 1065 102 FORMAT(12X,'I',T125,'I',/,12X,'I',11(IX,D9.3),T125.'I') (/3

1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 End of file 108 FORMAT(IX,124('-')) WRITE(8,105) 105 FORMAT(5X,'Y',6X,'I',T125,'I') WRITE(8,106) 106 FORMAT(1X,'COORDINATE',1X,'I',T125,'I') DO 80 J=1.6 1d=7-J 80 WRITE(8,104) Y1(J1),(TOU(J1,I).I=1,11) WRITE(8 108) WRITE(8,102) (X1(I),I=1, 11) WRITE(8,101) WRITE(8, 107) 104 FORMAT(12X,'I'.T125.'I',/,IX,D9.3,2X,'I',11(IX,D9.3),T125,'I') RETURN END %D J^N

1081 SUBROUTINE TEMP(N,NT,TX,MSIZE,EL,B,H,E,UN,ALPHA,TIX,TIY) 1082 C 1083 C TEMP CALCULATES THE EFFECTS OF A TEMPERATURE GRADIENT, IF PRESENT, 1084 C ON BLADE STRESS THROUGH ITS CONTRIBUTION TO THE RIGHT HAND SIDE 1085 C VECTOR, WHICH IS FORMED IN SUBROUTINE RHS. 1086 IMPLICIT REAL*8(A-H,O-Z) 1087 DIMENSION TX(35),TIX(35).TIY(35) 1088 DO 85 dJ=1,MSIZE 1089 TIX(JJ)=O.OD+00 1090 85 TIY(dJ)=O.OD+00 1091 IF(NT.EO. 0) GO TO 25 1092 KK=1 1093 LL=O 1094 DO 70 JJ=1,MSIZE 1095 DTX=O.OD+00 1096 DTY=O.OD+-00 1097 l1=2 1098 DO 60 I=1,NT 1099 IF(KK.EO. 0) GO TO 35 1100 DI=-DFLOAT(KK)/DFLOAT((I+KK)*(LL+1)) 1101 DTX=TX(d1)*EL**(I+KK)*B**(LL+1)*DI+DTX 1102 35 IF(LL.EQ. O) GO TO 60 1103 D2=-1.0/DFLOAT(I+KK+1) 1104 DTY=DTY+D2*TX(d1)*EL**(I+KK+1)*B**LL 1105 60 J1=d1+1 1106 DO 50 J=1,NT 1107 DO 55 I=1,NT 1108 IF(I+d.GT. NT) GO TO 50 1109 IF(KK.EO. O) GO TO 75 1110 D1=-DFLOAT(KK)/DFLOAT((I+KK)*(d+LL+1)) 1111 DTX=TX(d1)*EL**(I+KK)*B**(J+LL+1)*D1+DTX 1112 75 IF(LL.EO. O) GO TO 55 1113 D2=-DFLOAT(LL)/DFLOAT((I+KK+I)*(J+LL)) 1114 DTY=TX(J1)*EL**(I+KK+1)*B*(d+LL)*D2+DTY 1115 55 d1=d1+1 1116 50 CONTINUE 1117 DO 40 d=1,NT 1118 IF(KK.EQ. O) GO TO 65 1119 D1=-1.0/DFLOAT(d+LL+1) 1120 DTX=DTX+D1*TX(J1)*EL**KK*B**(dJ+LL+1) 1121 65 IF(LL.EQ. O) GO TO 40 1122 D2=-DFLOAT(LL)/DFLOAT((KK+1)*(J+LL)) 1123 DTY=TX(J1)*EL**(KK+1)*B**(d+LL)*D2+DTY 1124 40 d1=J1+1 1125 TIX(dJ)=E*ALPHA*H*DTX/((1.0D+OO)-UN) 1126 TIY(Jd)=E*ALPHA*H*DTY/((1.OD+00)-UN) 1127 IF(KK.EQ. 0) GO TO 20 1128 KK=KK+1 1129 IF(KK+LL.GT. N) GO TO 15 1130 GO TO 70 1131 15 LL=LL+1 1132 KK=1 1133 IF(KK+LL.GT. N) GO TO 10 1134 GO TO 70 1135 10 KK=O 1136 LL=1 1137 GO TO 70 1138 20 LL=LL+1 1139 IF(LL.GT. N) GO TO 2 1140 KK=O '10 Ul

1141 1142 1143 End of file 70 CONTINUE 25 RETURN END 0% tlo

1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 End of file C C C C C C C SUBROUTINE FORCE(N,R.FTANFNOR,XIX2.EL,BFORCX,FORCYMSIZE) FORCE CALCULATES THE BLADE STRESSES CAUSED BY A NORMAL. TANGENTIAL AND/OR TENSILE LOADS THROUGH ITS CONTRIBUTION TO THE RIGHT HAND SIDE VECTOR, WHICH IS FORMED IN SUBROUTINE RHS. IMPLICIT REAL*8(A-H,O-Z) DIMENSION X1(2),X2(2),FORCX(35),FORCY(35) MSIZE IS THE NUMBER OF TERMS IN THE APPROXIMATING FUNCTION FOR DISPLACEMENT FIELD, U. MSIZE=2*N DO 10 I=1,N 10 MSIZE=N-I+MSIZE J=1 DO 20 L1=1,N LL=L1-1 DO 30 KK=1,N IF(KK+LL.GT. N) GO TO 20 FORCX(J)=-R*EL**KK*B**(LL+1)/DFLOAT(LL+I) $+FTAN*(X2(2)**(KK+1)-X1(2)**(KK+1))*B**LL/DFLOAT(KK+1) FNI=FNOR*(X2(1)**(KK+1)-X1(1)**(KK+1))*B**LL/DFLOAT(KK+1) IF(LL.NE. O) GO TO 40 DI=DFLOAT(KK+1) D2=DFLOAT(N+1)*EL**(N-KK) D3=X2(1)**(KK+1)-X1(1)**(KK+1) D4=X2(1)**(N+1)-XI(1)**(N+1) FNI=FNOR*((D3/D1)-(D4/D2)) IF(KK.EQ. N) FNI=0.0D+00 40 FORCY(J)=FNI 30 J=J+1 20 CONTINUE DO 50 LL=1,N FORCX(J)=FTAN*(X2(2)-X1(2))*B**LL FORCY(J)=FNOR*(X2(1)-X1(1))*B**LL 50 J=J+1 RETURN END 0D.,,.j

1181 1182 C 1183 C 1184 C 1185 C 1186 1187 1188 1189 190 1191 1192 1193 1194 1195 1196 1197 End of file SUBROUTINE RHS(FORCX,FORCY,TIX.TIY,MSIZE,N,BX) RHS COMBINES THE EFFECTS OF THE EXTERNAL LOADS, STORED IN FORCX AND FORCY, WITH THE EFFECTS OF A TEMPERATURE GRADIENT, STORED IN TIX AND TIY, TO FORM THE RIGHT HAND SIDE VECTOR, WHICH IS STORED IN BX. IMPLICIT REAL*8(A-H,O-Z) DIMENSION FORCX(35).FORCY(35),TIX(35).TIY(35),BX(70) DO 10 J=1,MSIZE d1=J+MSIZE J2=d IF(J.GE. N) d2=J2+1 BX(J)=FORCX(d)+TIX(J) IF(J2.GT. MSIZE) GO TO 10 BX(l1)=FORCY(d2)+TIY(d2) 10 CONTINUE RETURN END oo 00

1198 SUBROUTINE SDRIVE(AX,NTS,KOUT) 1199 C SUBROUTINE SDRIVE DOES NOT WORK PROPERLY AT THIS TIME. 1200 C SDRIVE IS USED TO LET THE USER ENTER X, AND Y COORDINATES AND 1201 C PRESTRESS OR STRESS DATA FOR VARIOUS POINTS ON THE BAND SAW BLADE. 1202 C THUS, STRESS FIELD COEFFICIENTS CAN BE CALCULATED IF ONLY THE 1203 C VIBRATION PROBLEM IS BEING SOLVED. PRESTRESS FIELDS COEFFICIENTS 1204 C CAN BE CALCULATED TO BE COMBINED WITH THE STRESS FIELD COEFFICIENTS 1205 C LATER IN THE PROGRAM. 1206 IMPLICIT REAL*8(A-H,O-Z) 1207 DIMENSION AX(35,3),XT(70),YT(70),SXI(70),SX(35) 1208 IPI=O 1209 10=5 1210 DO 5 1=1,3 1211 DO 7 d=1,35 1212 7 AX(I,d)=O.OD+00 1213 5 CONTINUE 1214 WRITE(6,102) 1215 102 FORMAT(///,3X,'BEGIN ENTRY OF STRESS DATA') 1216 15 WRITE(6,103) 1217 103 FORMAT(//,6X,'ENTER THE NUMBER OF THE DESIRED STRESS FIELD',/,9X, 1218 1'1- SIGMA X',/,9X,'2- SIGMA Y',/,9X,'3- TAU XY',/,9X. 1219 2'4- ALL THREE STRESS FIELDS'./,9X,'5- STOP') 1220 READ(5,100) IS 1221 100 FORMAT(I1) 1222 IF(IS.NE. 4) GO TO 30 1223 IPI=1 1224 IS=1 1225 30 IF(IS.EQ. 5) GO TO 20 1226 WRITE(6,104) 1227 104 FORMAT(/,12X,'IS DATA TO BE READ FROM A DATA FILE',/,12X,'O- NO, 1227.2 1,'1- YES') 1228 READ(5,100) IR 1229 IF(IR.EQ. 1) 10=10 1230 IF(IO.EQ. 10) REWIND 10 1231 WRITE(5,'105) 1232 105 FORMAT(/,9X,'ENTER THE NUMBER OF DATA POINTS (12 FORMAT)') 1233 READ(5,101) NPTS 1234 101 FORMAT(I2) 1235 40 DO 45 I1=1,NPTS 1236 IF(IO.EQ. 5) WRITE(6,106) I1 1237 106 FORMAT(12X.'ENTER X-COORDINATE, Y-COORDINATE, AND STRESS FOR POINT 1238 1 ',12,' ON THE BAND SAW BLADE') 1239 45 READ(IO,201) XT(I1),YT(I1),SXI(I1) 1240 201 FORMAT(3D15.6) 1241 CALL TCOFF(NTS,NPTS,XT,YT,SXI,SX,KOUT) 1242 MT=(NTS+1)*(NTS+2)/2 1243 DO 50 IJ=1,MT 1244 50 AX(IJ,IS)=SX(IJ) 1245 IS=IS+1 1246 IF(IS.EQ. 4) GO TO 20 1247 IF(IPI.EQ. 1) GO TO 40 1248 GO TO 15 1249 20 RETURN 1250 END End of file

VIBRATION ANALYSIS SUBROUTINES 100

1 SUBROUTINE VIBRA(N,MNM,MAXNM,KON,L,B,H,D,RHOH.NU,C,AX ITITL) 2 C 3 C SUBROUTINE VIBRA SOLVES FOR THE FIVE LOWEST EIGENVALUES OF 4 C THE BAND SAW BLADE,BASED ON THE CLASSICAL RITZ METHOD CODE 5 C PRESENTED IN THE REPORT "TWO COMPUTER CODES FOR BAND SAW 6 C BLADE VIBRATION AND STABILITY ANALYSIS" BY A.G. ULSOY, 1978. 7 C THE THEORY IS PRESENTED IN THE DOCTORAL DISSERTATION BY 8 C A. G. ULSOY "VIBRATION AND STABILTY OF BAND SAW BLADES: A 9 C THEORETICAL AND EXPERIMENTAL STUDY". 10 C 11 C VIBRA USES NON DIMENSIONAL INTEGRALS STORED ON TAPE(FILE) 12 C TO COMPUTE THE MASS, GYROSCOPIC AND STIFFNESS MATRICES. 13 C VIBRA CALLS SUBROUTINES DECODE, REFORM, EIGENS. 14 C 15 C SUBROUTINE DECODE: PROVIDES A CONVERSION BETWEEN THE INDEXING 16 C NOTATION USED IN VIBRA AND THAT USED FOR 17 C THE NONDIMENSIONAL INTEGRALS ON TAPE(FILE). 18 C SUBROUTINE REFORM: CALCULATES A MATRIX (A) FROM THE MASS, 19 C GYRO, AND STIFF MATRICES TO CONVERT FROM 20 C A SYSTEM OF SECOND ORDER DIFFERENTIAL 21 C EQUATIONS TO A SYSTEM OF FIRST ORDER 22 C DIFFERENTIAL EQUATIONS. 23 C SUBROUTINE EIGENS: CALCULATES THE FIVE LOWEST EIGENVALUES OF 24 C THE BLADE, USING THE MATRIX (A) FORMED 25 C IN SUBROUTINE REFORM. 26 C 27 C ********** ********************** ***************** 28 C ************************************************ * 29 C ** ** 30 C ** IMPORTANT VARIBLES IN SUBROUTINE VIBRA ** 31 ** ** 32 C ** KON(8): PREVENTS REREADING OF NONDIMENSIONAL ** 33 C ** INTEGRALS AFTER INITIAL READING. ** 34 C ** ** 35 C ** N: NUMBER OF SIMPLE SUPPORTED BEAM EIGENFUNCTIONS. ** 36 C ** ** 37 C M: NUMBER OF FREE-FREE BEAM EIGENFUNCTIONS. 38 C ** ** 39 C ** NM: N TIMES M. ** 40 ** ** 41 C ** AX: STRESS FIELD COEFFICIENTS FOR SIGMA-X, ** 42 C ** SIGMA-Y, AND TAU-XY. ** 43 C ** ** 44 C ************** ************ ******************* 45 C *******t******************** ******* ************** 46 C 47 C 48 IMPLICIT REAL*8(A-HO-Z) 49 REAL*8 L.NU,MASS(36) 50 DIMENSION ITITL(50),KON(11),GYRO(36,36),A(72,72),STIF(36,36) 51 DIMENSION AX(35,3),SX01(6,6),SX02(6.6),SX03(6,6),SX04(6,6), 52 $SX05(6,6),SX06(6,6),SX07(6,6),SX08(6,6).SX09(6,6),SX10(6,6), 53 $SX11(6,6),SX12(6,6),SX13(6.6),SX14(6,6),SY01(6,6),SY02(6.6), 54 $SY03(6,6),SY04(6,6),SY05(6,6),SY06(6,6),SY07(6,6),SY08(6,6), 55 $SY09(6,6),SY10(6,6),SY11(6,6),SY12(6,6),SY13(6,6),SY14(6,6), 56 $SXT1(6),SXT2(6),SXT3(6),SXT4(6),SYT1(6),SYT2(6),SYT3(6),SYT4(6) 57 IF(KON(8).EQ. 1) GO TO 60 58 REWIND 7 59 C 60 C READ NONDIMENSIONAL INTEGRALS FROM TAPE (FILE). 0 H

61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 i00 101 102 103 104 105 106 107 108 109 110 112 112 113 End of fi1i READ(7) SXOl,SX02,SX03,SX04,SXO5,SXO6,SX07 READ(7) SX08,SX09,SXO1,SX11,X122,SX13,SX14 READ(7) SY01,SY02,SY03,SY04,SYO5,SY06,SY07 READ(7) SY08.SY09,SYIO,SY11,SY12,SY13,SY14 READ(7) SXTI,SXT2,SXT3,SXT4,SYT1,SYT2,SYT3,SYT4 C FORMULATE THE MASSGYRO,AND STIFFNESS MATRICES 60 CONTINUE DO 50 I=1,NM CALL DECODE(I,N.M,IX,IY) MASS(I)=RHOH*L*B*SXOi(IX,IX)*SYO1(IY,IY) DO 50 J=1,NM CALL DECODE(J,N,M,dX.JY) GYRO(I,J)=(2.OD+0)*RHOH*B*C*SX02(IX,JX)*SY01(IY,JY) $+L*B*SXO1(IX,dX)*SYOj(IYdY) STIF1=(B/L)*(AX(1,1)*SX04(IX,JX)*SY01(IY,JY)-AX(2,1)* 1 SX08(IX,JX)*SY01(IY,JY)-AX(8,1)*SX03(IX,JX)*SY06(IY,JY)2 AX(3,1)*SX11(IX,dX)*SY01(IY,JY)-AX(5,1)*SX08(IX,dX)* 3 SY06(IY,dY)-AX(9,1)*SX03(IXJX)*SY09(IY,JY)-AX(4,1)* 4 SX14(IX,dX)*SY01(IY,JY)-AX(6,1)*SX11(IX,dX)*SY06(IY,dY)5 AX(7,1)*SX08(IX.JX)*SY09(IY,JY)-AX(10.1)*SX03(IXdX.)* 6 SY12(IY,dY)) STIF2=(L/B)*(AX(1,2)*SX01(IX,JX)*SY04(IY,JY)-AX(2,2)*SX06(IX,dX) 1 *SY03(IY,JY)-AX(8,2)*SX01(IXJX)*SY08(IY,dY)-AX(3,2)* 2 SX09(IX,JX)*SY03(IY,dY)-AX(5,2)*SX06(IXdX)*SY08(IY,dY)3 AX(9,2)*SX01(IX,dX)*SY11(IY,JY)-AX(4,2)*SX12(IX,JX)* 4 SY03(IY,dY)-AX(6.2)*SX09(IXJX)*SY08(IY,dY)-AX(7,2)* 5 SX06(IX,dX)*SY11(IY,dY)-AX(10,2)*SX01(IX,JX)*SY14(IY,JY)) STIF3=(2.OD+00)*(AX(1,3)*SX02(JX,IX)*SY02(IYdY)*(5.0D-01)+ 1 AX(1,3)*SX02(IX,dX)*SY02(JY,IY)*(5.0D-01)-AX(2,3)* 2 SX07(IX,JX)*SY02(IY,dY)-AX(8.3)*SX02(IX,dX)*SY07(IY.JY)$ AX(3,3)*SX1O(IX,JX)*SY02(IYdY) 3 -AX(5,3)*SX07(IX,dX)*SY07(IY,dY)-AX(9,3)*SX02(IXdX)* 4 SY1O(IY,dY)-AX(4,3)*SX13(IXdX)*SY02(IYdY)-AX(6,3)* 5 SX10(IX,JX)*SY07(IY,dY)-AX(7,3)*SX07(IX,JX)*SY10(IY,dY)6 AX(10,3)*SX02(IX,dX)*SY13(IY,dY)) STIF(I,J)=D*((B/(L**3))*SX05(IXdX)*SY01(IYdY)+(L/(B**3))* 1 SX01(IX,dX)*SY05(IYdY)+(((2.0D+00)*((1.0D+00)-NU))/(L*B))* 2 SX04(IX,dX)*SY04(IY,dY)+(NU/(L*B))*(SX03(JX,IX)*SY03(IYdY)+ 3 SX03(IX,dX)*SY03(dY,IY)))+(STIF1+STIF2+STIF3)*H 0" t~i e C 3 50 CONTINUE KON(8)=1 NMT=2*NM MAXNMT=2*MAXNM SUBROUTINE REFORM REFORMULATES SET OF SECOND ORDER DIFFERENTIAL EQUATIONS AS A SET OF FIRST ORDER DIFFERENTIAL EQUATIONS. CALL REFORM(NM,NMT,MASS,GYRO,STIF,A) SOLVES THE EIGENVALUE PROBLEM FOR THE FIVE LOWEST EIGENVALUES OF THE BLADE. CALL EIGENS(NMT,MAXNMT,KON,L,RHOH,D,A,ITITL) RETURN END

<14 lUOKUUIINt UttUUt[1,N,M,IA,I~y 115 C 116 C USED TO PROVIDE A CONVERSION BETWEEN THE INDEXING NOTATION 117 C USED IN VIBRA AND THAT USED FOR THE NONDIMENSIONAL INTEGRALS 118 C STORED ON TAPE (FILE). 119 C 120 DO 10 K=1,N 121 KMI=K-1 122 IF(I.LE. (KMI*M).OR. I.GT. (K*M)) GO TO 10 123 IX=K 124 IY=I-KMI*M 125 RETURN 126 10 CONTINUE 127 RETURN 128 END End of file 0-' L,)

129 130 C 131 C 132 C 133 C 134 C 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 End of f le SUBROUTINE REFORM(NM,NMT,MASS,GYRO,STIF,A) CALCULATES A MATRIX (A) FROM THE MASS, GYRO, AND STIFF MATRICES TO CONVERT FROM A SYSTEM OF SECOND ORDER DIFFERENTIAL EQUATIONS TO A SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS. IMPLICIT REAL*8(A-H,O-Z) DIMENSION GYRO(36,36),STIF(36,36),A(72,72) REAL*8 MASS(36) DO 10 I=1,NM II=I+NM DO 10 d=1,NM dJ=d+NM A(I,J)=-GYRO(I,J)/MASS(I) A(I,JJ)=-STIF(I,d)/MASS(I) A(II,d)=0.OD+00 IF(II.EO. dJ) A(II,J)=1.0D+00 A(IIJJ)=O.OD+00 10 CONTINUE RETURN END 0-j 0 $pb

150 151 C 152 C 153 C 154 C 155 C 156 C 157 C 158 C 159 C 160 C 161 C 162 C 163 C 164 C 165 C 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 End of file SUBROUTINE EIGENS(N,NM,KON.L,RHOH,D,A,ITITL) EIGENVALUES ARE OBTAINED USING EISPACK ROUTINES BALANC, ELMHES. AND HQR. THE FIVE LOWEST EIGENVALUES AND NONDIMENSIONAL NATURAL FREQUENCIES ARE PRINTED TO THE OUTPUT FILE. THE EIGEN SOLUTION ROUTINES USED ARE FOR THE CASE OF A REAL MATRIX (A) ** ** ** NOTE: SUBROUTINE HQR DESTROYS THE ORIGINAL ** ** MATRIX (A). ** ** ** ********************************************************** IMPLICIT REAL*8(A-H,O-Z) REAL*8 L DIMENSION A(72,72).SCALE(72),WR(72),WI(72),OMEGA(36),INT(72), $ITITL(50).KON(11),LEV(5) CALL BALANC(NM,N,A,LOW,IGH,SCALE) CALL ELMHES(NM.N,LOW,IGH.A,INT) CALL HQR(NM,NLOW.IGHA,WR.WI,IERR) IF(IERR.NE. O) WRITE(6,680) IERR WRITE(8,610) ITITL WRITE(8,611) NM1=N-1 DO 20 I=1,NM1,2 II=(I/2)+1 FREQ=DSORT(RHOH/D)*DSQRT((WR(I)**2)+(WI(I)**2)) OMEGA(II)=L*L*FREO 20 WRITE(8,620) I,WR(I),WI(I),OMEGA(II) DO 30 1=1,5 LEV(I)=1 OMEGAO=OMEGA(1) DO 40 d=i,NM1,2 JJ=(d/2)+1 IF(OMEGA(dJ).LT. OMEGAO) LEV(I)=d IF(OMEGA(JJ).LT. OMEGAO) OMEGAO=OMEGA(JJ) 40 CONTINUE, J=(LEV(I)/2)+1 OMEGA(dJ)=1.0D+20 30 CONTINUE WRITE(8,625) (LEV(I),I=1,5) 610 FORMAT(1H1,/,IH-,50A1,/,1H-,4X,'EIGEN',6X,'REAL',9X,'IMAGINARY' $,9X,'NON-DIMENSIONAL') 611 FORMAT(SX,'VALUE',6X,'PART',12X,'PART',9X,'NATURAL FREQUENCY') 620 FORMAT(5X,I3,2X,2(1X,D14.6).5X,D14.6) 625 FORMAT(3X,'LOWEST 5 EIGENVALUES IN INCREASING ORDER ',5(2X,I2)) 680 FORMAT(//,'ERROR**** SUBROUTINE HQR HAS FAILED TO CONVERGE', $' TO EIGENVALUE NUMBER',I5) RETURN END H Ln

1 C NAASA 3.2.001 BALANC FTN 06-24-75 THE UNIV OF MICH COMP CTR 2 C 3 C -------------------- ------------ 4 C 5 SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) 6 C 7 INTEGER I,J,K,L,M,N,.d,NM,IGH,LOW,IEXC 8 REAL*8 A(NM,N),SCALE(N) 9 REAL*8 C,F,G,R,S,B2,RADIX 10 REAL*8 DABS 11 LOGICAL NOCONV 12 C 13 C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, 14 C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. 15 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). 16 C 17 C THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES 18 C EIGENVALUES WHENEVER POSSIBLE. 19 C 20 C ON INPUT: 21 C 22 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 23 C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 24 C DIMENSION STATEMENT; 25 C 26 C N IS THE ORDER OF THE MATRIX; 27 C 28 C A CONTAINS THE INPUT MATRIX TO BE BALANCED. 29 C 30 C ON OUTPUT: 31 C 32 C A CONTAINS THE BALANCED MATRIX; 33 C 34 C LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) 35 C IS EQUAL TO ZERO IF 36 C (1) I IS GREATER THAN J AND 37 C (2) J=,.....LOW-1 OR I=IGH+1,....N; 38 C 39 C SCALE CONTAINS INFORMATION DETERMINING THE 40 C PERMUTATIONS AND SCALING FACTORS USED. 41 C 42 C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH 43 C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED 44 C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS 45 C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN 46 C SCALE(J) = P(d), FOR J = I.....LOW-1 47 C = D(J,J), J = LOW,....IGH 48 C = P(d) d = IGH+1....N. 49 C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, 50 C THEN 1 TO LOW-1. 51 C 52 C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. 53 C 54 C THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN 55 C BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS 56 C K,L HAVE BEEN REVERSED.) 57- C 58 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, 59 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 60 C 0 <0~

61 C ----------------------------------------------- ---------- 62 C 63:::::::::: RADIX IS A MACHINE DEPENDENT PARAMETER SPECIFYING 64 C THE BASE OF THE MACHINE FLOATING POINT REPRESENTATION. 65 C RADIX = 16.0DO FOR LONG FORM ARITHMETIC 66 C ON S360:::::::::: 67 DATA RADIX/Z4210000000000000/ 68 C 69 B2 RADIX * RADIX 70 K = 1 71 L N 72 GO TO 100 73 C:::::::::: IN-LINE PROCEDURE FOR ROW AND 74 C COLUMN EXCHANGE:::::::::: 75 20 SCALE(M) = d 76 IF (d.EQ. M) GO TO 50 77 C 78 DO 30 I = 1, L 79 F = A(I,J) 80 A(I,J) = A(I,M) 81 A(I.M) = F 82 30 CONTINUE 83 C 84 DO 40 I = K, N 85 F = A(J,I) 86 A(J,I) = A(M,I) 87 A(M,I) = F 88 40 CONTINUE 89 C 90 50 GO TO (80,130), IEXC 91 C:::::::::: SEARCH FOR ROWS ISOLATING AN EIGENVALUE 92 C AND PUSH THEM DOWN:::::::::: 93 80 IF (L.EO. 1) GO TO 280 94 L = L - 1 95 C::::::::: FOR J=L STEP -1 UNTIL 1 DO --:::::::::: 96 100 DO 120 JJ = 1, L 97 d = L + 1 - dJ 98 C 99 DO 110 I = 1, L 100 IF (I.EQ. d) GO TO 110 101 IF (A(d,I).NE. O.ODO) GO TO 120 102 110 CONTINUE 103 C 104 M = L 105 IEXC = 1 106 GO TO 20 107 120 CONTINUE 108 C 109 GO TO 140 110 C:::::::::: SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE 111 C AND PUSH THEM LEFT:::::::::: 112 130 K = K + 1 113 C 114 140 DO 170 d = K, L 115 C 116 DO 150 I = K, L 117 IF (I.EQ. d) GO TO 150 118 IF (A(I,J).NE. O.ODO) GO TO 17C 119 150 CONTINUE 0 -j

121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 End of file 4 M = K IEXC = 2 GO TO 20 170 CONTINUE C:::::::::: NOW BALANCE THE SUBMATRIX IN ROWS K TO L:::::::::: DO 180 I = K, L 180 SCALE(I) = 1.000 C:::::::::: ITERATIVE LOOP FOR NORM REDUCTION:::::::::: 190 NOCONV =.FALSE. C DO 270 I = K, L C = O.ODO R = O.ODO C DO 200 J = K. L IF (J.EQ. I) GO TO 200 C = C + DABS(A(J,I)) R = R + DABS(A(I,J)) 200 CONTINUE C G = R / RADIX F = 1.0DO S =C + R 210 IF (C.GE. G) GO TO 220 F = F * RADIX C = C * B2 GO TO 210 220 G = R * RADIX 230 IF (C.LT. G) GO TO 240 F = F / RADIX C = C / B2 GO TO 230 C:::::::::: NOW BALANCE:::::::::: 240 IF ((C + R) / F.GE. 0.95DO * S) GO TO 270 G = 1.0DO / F SCALE(I) = SCALE(I) * F NOCONV =.TRUE. H0 00 ( ( ( ( ( C DO 250 J = K, N 250 A(I,J) = A(I,J) * G C 260 DO 260 J = 1, L A(J,I) = A(JI) * F C.~ 270 CONTINUE IF (NOCONV) GO TO 190 280 LOW = K IGH = L RETURN C:::::::::: LAST CARD OF END p BALANC::::::::::

1/4 U NAASA 3.2.006 ELMHES FTN 06-24-75 THE UNIV UF MICH CUMP CIR 175 C 176 C 177 C 178 SUBROUTINE ELMHES(NM,N,LOW,IGH.A,INT) 179 C 180 INTEGER I,,M,N,LA,NM,IGH,KPI,LOW,MM1,MP1 181 REAL*8 A(NM,N) 182 REAL*8 X,Y 183 REAL*8 DABS 184 INTEGER INT(IGH) 185 C 186 C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES, 187 C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. 188 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). 189 C 190 C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE 191 C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS 192 C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY 193 C STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. 194 C 195 C ON INPUT: 196 C 197 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 198 C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 199 C DIMENSION STATEMENT; 200 C 201 C N IS THE ORDER OF THE MATRIX; 202 C 203 C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING 204 C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, 205 C SET LOW=1, IGH=N; 206 C 207 C A CONTAINS THE INPUT MATRIX. 208 C 209 C ON OUTPUT: 210 C 211 C A CONTAINS THE HESSENBERG MATRIX. THE MULTIPLIERS 212 C WHICH WERE USED IN THE REDUCTION ARE STORED IN THE 213 C REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX; 214 C 215 C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS 216 C INTERCHANGED IN THE REDUCTION. 217 C ONLY ELEMENTS LOW THROUGH IGH ARE USED. 218 C 219 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, 220 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 221 C 222 C 223 C 224 LA = IGH - 1 225 KP1 = LOW + 1 226 IF (LA.LT. KP1) GO TO 200 227 C 228 DO 180 M = KP1, LA 229 MM1 = M - 1 230 X = O.ODO 231 I = M 232 C 233 DO 100 J = M, IGH 0 %:

234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 4 273 274 < 275 End of file IF (DABS(A(J,MM1)).LE. DABS(X)) GO TO 100 X = A(J,MM1) I = 100 CONTINUE C INT(M) = I IF (I.EQ. M) GO TO 130 C:::::::::: INTERCHANGE ROWS AND COLUMNS OF A:::::::::: DO 110 J = MMI, N Y = A(I,J) A(I,J) = A(M,J) A(M.J) = Y 110 CONTINUE C DO 120 d = 1, IGH Y = A(J,I) A(J,I) = A(J.M) A(J,M) = Y 120 CONTINUE C:::::::::: END INTERCHANGE:::::::::: 130 IF (X.EQ. O.ODO) GO TO 180 MP1 = M + 1 C DO 160 I = MP1, IGH Y = A(I,MM1) IF (Y.EQ. O.ODO) GO TO 160 Y = Y/ X A(I,MMI) = Y C DO 140 J = M. N A(I,J) = A(I,d) - Y * A(M,J) 140 o 0 — 0 C 150 160 DO 150 d = 1. IGH A(J.M) = A(d,M) + CONTINUE Y * A(J,I) C c c c 180 CONTINUE 200 RETURN::::::: LAST CARD OF ELMHES::::::: END

275 END 276 C NAASA 3.2.010 HQR FTN 06-24-75 THE UNIV OF MICH COMP CTR 277 C 278 C 279 C 280 SUBROUTINE HQR(NMN,LOW,IGH,H,WR,WI,IERR) 281 C 282 INTEGER I,J,K,L,M,N,EN,LL,MM,NA, NM,IGH, ITS,LOWMP2,ENM2,IERR 283 REAL*8 H(NM,N),WR(N),WI(N) 284 REAL*8 P,Q,R,S,T,W,X,Y,ZZ,NORM,MACHEP 285 REAL*8 DSQRT,DABS,DSIGN 286 INTEGER MINO 287 LOGICAL NOTLAS 288 C 289 C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, 290 C NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. 291 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). 292 C 293 C THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL 294 C UPPER HESSENBERG MATRIX BY THE QR METHOD. 295 C 296 C ON INPUT: 297 C 298 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 299 C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 300 C DIMENSION STATEMENT; 301 C 302 C N IS THE ORDER OF THE MATRIX; 303 C 304 C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING 305 C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, 306 C SET LOW=1, IGH=N; 307 C 308 C H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT 309 C THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG 310 C FORM BY ELMHES OR ORTHES, IF PERFORMED. IS STORED 311 C IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. 312 C 313 C ON OUTPUT: 314 C 315 C H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED 316 C BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND 317 C BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED; 318 C 319 C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, 320 C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES 321 C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS 322 C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE 323 C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN 324 C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT 325 C FOR INDICES IERR+1,....N; 326 C 327 C IERR IS SET TO 328 C ZERO FOR NORMAL RETURN, 329 C d IF THE J-TH EIGENVALUE HAS NOT BEEN 330 C DETERMINED AFTER 30 ITERATIONS. 331 C 332 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, 333 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 334 C i-. I-i H

335 C -- ----- ------------ ------------- 336 C 337 C:::::::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING 338 C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. 339 C MACHEP = 16.ODO**(-13) FOR LONG FORM ARITHMETIC 340 C ON 5360:::::::::: 341 DATA MACHEP/Z3410000000000000/ 342 C 343 IERR = 0 344 NORM = O.ODO 345 K = 1 346 C:::::::::: STORE ROOTS ISOLATED BY BALANC 347 C AND COMPUTE MATRIX NORM:::::::::: 348 DO 50 I = 1, N 349 C 350 DO 40 J = K, N 351 40 NORM = NORM + DABS(H(I,J)) 352 C 353 K = I 354 IF (I.GE. LOW.AND. I.LE. IGH) GO TO 50 355 WR(I) = H(I,I) 356 WI(I) = O.ODO 357 50 CONTINUE 358 C 359 EN = IGH 360 T = O.ODO 361 C:::::::::: SEARCH FOR NEXT EIGENVALUES:::::::::: 362 60 IF (EN.LT. LOW) GO TO 1001 363 ITS = 0 364 NA = EN - 1 365 ENM2 = NA - 1 366 C:::::::::: LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT 367 C FOR L=EN STEP -1 UNTIL LOW DO --:::::::::: 368 70 DO 80 LL = LOW, EN 369 L = EN + LOW - LL 370 IF (L.EQ. LOW) GO TO 100 371 S = DABS(H(L-1,L-1)) + DABS(H(L.L)) 372 IF (S.EQ. O.ODO) S = NORM 373 IF (DABS(H(L,L-1)).LE. MACHEP * S) GO TO 100 374 80 CONTINUE 375 C.:::::::::: FORM SHIFT:::::::::: 376 100 X = H(EN,EN) 377 IF (L.EQ. EN) GO TO 270 378 Y = H(NA,NA) 379 W = H(EN,NA) * H(NA,EN) 380 IF (L.EQ. NA) GO TO 280 381 IF (ITS.EQ. 30) GO TO 1000 382 IF (ITS.NE. 10.AND. ITS.NE. 20) GO TO 130 383 C:::::::::: FORM EXCEPTIONAL SHIFT:::::::::: 384 T = T + X 385 C 386 DO 120 I = LOW, EN 387 120 H(I,I) = H(I,I) - X 388 C 389 S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) 390 X = 0.75DO * S 391 Y = X 392 W = -0.4375DO * S * S 393 130 ITS = ITS + I )-J K)

395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 44 445 446 447 448 449 450 451 452 453 454 C C SUB-DIAGONAL ELEMENTS. FOR M=EN-2 STEP -1 UNTIL L DO —::::::::: DO 140 MM = L, ENM2 M = ENM2 + L - MM ZZ = H(M,M) R = X - ZZ S = Y - ZZ P = (R * S - W) / H(M+I,M) + H(M,M+1) Q = H(M+1,M+1) - ZZ - R - S R = H(M+2,M+1) S = DABS(P) + DABS(Q) + DABS(R) P = P/ S Q = /S R = R / S IF (M.EQ. L) GO TO 150 IF (DABS(H(M,M-1)) * (DABS(Q) + DABS(R)).LE. MACHEP * DABS(P) X * (DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1.M+1)))) GO TO 150 140 CONTINUE C 150 MP2 = M + 2 C DO 160 I - MP2, EN H(I,I-2) = O.ODO IF (I.EO. MP2) GO TO 160 H(I,I-3) = 0.000 160 CONTINUE C:::::::::: DOUBLE QR STEP INVOLVING ROWS L TO EN AND C COLUMNS M TO EN:::::::::: DO 260 K = M, NA NOTLAS = K.NE. NA IF (K.EQ. M) GO TO 170 P = H(K,K-1) O = H(K+1,K-1) R = O.ODO IF (NOTLAS) R = H(K+2,K-1) X = DABS(P) + DABS(Q) + DABS(R) IF (X.EQ. 0.000) GO TO 260 P =P / X Q = Q/ X R = R / X 170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P) IF (K.EQ. M) GO TO 180 H(K,K-1) = -S * X GO TO 190 180 IF (L.NE. M) H(K,K-1) = -H(K,K-I) 190 P = P + S X= P/ S Y =0 / S ZZ = R / S 0 = Q / P R = R / P C:::::::::: ROW MODIFICATION:::::::::: DO 210 d = K. EN P = H(K,J) + 0 * H(K+1,J) IF (.NOT. NOTLAS) GO TO 200 P = P +. R * H(K+2,d) H(K+2,d) = H(K+2,d) - P * ZZ 200 H(K+1,J) = H(K+1,.) - P * Y H(K,J) = H(K,d) - P * X 210 CONTTMlIF F-J

455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 End of file C d = MINO(EN,K+3) C:::::::::: COLUMN MODIFICATION:::::::::: DO 230 I = L, d P = X * H(I,K) + Y * H(I,K+1) IF (.NOT. NOTLAS) GO TO 220 P = P + ZZ * H(I,K+2) H(I,K+2) = H(I,K+2) - P * R 220 H(I,K+1) = H(I,K+1) - P * Q H(I,K) = H(I,K) - P 230 CONTINUE C 260 CONTINUE C GO TO 70 C:::::::::: ONE ROOT FOUND:::::::::: 270 WR(EN) = X + T WI(EN) = 0.ODO EN = NA GO TO 60 C:::::::::: TWO ROOTS FOUND:::::::::: 280 P = (Y - X) / 2.000 Q = * P + W ZZ = DSQRT(DABS(Q)) X = X + T IF (Q.LT. O.ODO) GO TO 320 C:::::::::: REAL PAIR:::::::::: ZZ = P + DSIGN(ZZ,P) WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ.NE. O.ODO) WR(EN) = X - W / ZZ WI(NA) = O.ODO WI(EN) = O.ODO GO TO 330 C:::::::::: COMPLEX PAIR::::::::: 320 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ 330 EN = ENM2 GO TO 60 C:::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS:::::::::: 1000 IERR = EN 1001 RETURN C:::::::::: LAST CARD OF HQR:::::::::: END ~ F-J HH^

UNIVERSITY OF MICHIGAN 3 9015 02527 8105