MSU LIBRARIES .—::—. RETURNING MATERIALS: PIace in book drop to remove this checkout from your record. FINES will be charged if book is returned after the date stamped below. QUALITATIVE ANALYSIS OF A DAHPED LINEAR DYNAHIC SYSTEH BY Jane L. Rankine A THESIS Submitted to flichigan State University in partial fulfillment ot the requirements for the degree of MASTER OF SCIENCE Department of flechanical Engineering 1985 38% £2.29 ABSTRACT QUALITATIVE ANALYSIS OF A DAHPED LINEAR DYNAHIC SYSTEM BY Jane L. Hawkins Frequently in system analysis a qualitative description of the system response can lead to useful insight into the solution. In the case of a damped linear dynamic system. knowledge of whether the system is overdamped. underdamped or a combination may be adequate. thus making the more expensive calculation of the modes unnecessary. This thesis describes a method and presents software (QUALADS) for determining the qualitative damping nature of a system. given a system description in terms of mass. stiffness and damping matrices. Chapter I. INTRODUCTION . 2. GUALADS . . . 2.1 Example 1 2.2 Example 2 2.3 Example 3 2.4 Example 4 3. USERS GUIDE T0 3.1 Sample Run 3.2 General Comments and Instructions 4. QUALITATIVE VERSUS QUANTITATIVE ANALYSIS 4.1 One Degree of Freedom System 4.2 flultidegree of Freedom System TABLE OF CONTENTS Underdamped System Critically Damped System flixed Damped System Dverdamped System QUALADS . . 4.3 Qualitative Hethod by Inman and Andry S. OUALADS CORPUTATIONAL INFORNATION 5.1 flethods of Computation 5.2 Computation Time 5. SUHHARY . . . 6.1 Conclusion 6.2 Extensions List of References Appendix - QUALADS Source 111 Page 16 2O 2O 23 24 27 27 29 31 31 32 33 34 Figure 1. Figure 2. Figure 3. LIST OF FIGURES Two Degree of Freedom System One Degree of Freedom System CPU Time Plot . . . . . iv Chapter 1 INTRODUCTION A primary question in the analysis of a damped linear dynamic system is whether or not the system can exhibit oscillatory behavior in free motion. This question can be answered quantitatively by calculating the eigenvalues of the differential equations of motion. The process is generally complex and burdensome for large systems. The question can also be answered qualitatively by determining the system damping type. Normally a qualitative analysis requires less time and effort than a quantitative analysis. making it a more attractive method for answering the question as stated above or for gaining preliminary information about a system’s dynamics. There are several techniques for making a qualitative analysis of a damped system. The method derived by Inman and Andry [1.2] proves to be one of the more computationally efficient and comprehensive techniques. It defines a matrix. termed the critical damping matrix. which. when compared with the system’s damping matrix. generates the system’s damping type. This thesis describes the software QUALADS (QUALitativs Analysis of Damped Systems). which uses the technique of Inman and Andry to determine the damping type of 55-" a damped linear system described by a set of second order differential equations. Chapter 2 describes the use of QUALADS and gives some simple examples. Chapter 3 presents an example run and users guide for QUALADS. Chapter 4 is a general discussion of quantitative versus qualitative analyses with a more detailed description of the Inman and Andry algorithm. Chapter 5 explains some of the methods of computation used by QUALADS and computation time for various segments. The final chapter provides a summary and suggestions for additional applications and extensions. The appendix contains the source code for QUALADS. Chapter 2 QUALADS OVERVIEU The QUALADS program performs a qualitative analysis of a damped linear dynamic system. The system must be asymptotically stable and defined by a set of second order differential equations. The inputs for QUALADS are the three coefficient matrices of the differential equations. commonly termed the mass. damping and stiffness matrices. The mass and stiffness matrices must be positive definite and symmetric and the damping matrix must be positive semi-definite and symmetric. The mass and stiffness matrices are used to compute a matrix called the critical damping matrix. By examining the matrix formed by the difference of the damping and critical damping matrices. QUALADS determines whether the system is overdamped. critically damped. underdamped or mixed damped. These damping types are the qualitative description which is the output from QUALADS. In the process of determining the damping type. QUALADS computes the undamped eigenvalues and eigenvectors as well as the critical damping matrix. If the user desires quantitative results then the damped eigenvalues may be computed as an additional option. Consider the following 2 degree of freedom spring-mass-damper system. x1 x2 r—’ r———’ k2= 1.5 N/m k3= l N/m —/\/\/\/~1.. 4\/\/\/— , k. *vav— A __J_]__ 3 c1 c2 c3 // ///// \\\\\\\: Lil—J Figure 1. Two Degree of Freedom System The mass. damping and stiffness matrices for this system are: 1 0 c1+c2 -c2 2.5 -l.5 [fl] = [C] = [K] = O 1 -c2 c2+c3 -1.5 2.5 Depending on the values of c1. c2 and c3 the system will exhibit a particular damping type response. 2.1 xam e - Underdam ed s at m If c1=c2=c3=1 kg/sec. then the damping matrix would be: 2 -l [C] 3 -1 2 After inputting the mass. stiffness and damping matrices. OUALADS would respond with NIN‘INNNNNINNNN‘IN‘I‘INNNINI‘INN The system is underdamped nunuauuaananuuauiaaeauaaara As a check the damped eigenvalues can be computed and are as follows. A 1.2 -- -.s 1 .856i A3.4 = -1.5 1; 1.3231 Since the eigenvalues are complex. both modes will underdamped and thus the system will be underdamped. predicted by QUALADS. 2. a 2 - iti all dam ed a stem If c1=c3=2 kg/sec and c2=1 kg/sec.then 3 -l [C] = -l 3 and QUALADS outputs *‘INN‘INI‘I!NININNNNNIINNINNN‘INN‘INN! The system is critically damped urn»!ninafinunnseruuuineiinureter» The damped eigenvalues for this system. as shown below. indicate that the system is critically damped. A1v2 = -l.O lam a -2.0 they be also Since this system is critically damped. the damping matrix above must be equal to the critical damping matrix. In accord. QUALADS computes the critical damping matrix to be: 3 -1 [C] = -1 3 2.3 Examplg § - flixed damped system Letting cl=c3=2.5 kg/sec and c2=.5 kg/sec will form a damping matrix equal to: [C] = I (II 1...) Using these damping conditions. GUALADS predicts NNNNNNINNNI’N‘IINNIN‘INI*‘I‘I‘l‘lilflfl The system has mixed damping Iflflflfllflflflflfiflfllflflflfllflfiiflfllliiii The damped eigenvalues are A1 = -.5 . A2= -2.0 13.4 = —1.75 1 .9551 Since one pair of eigenvalues is real and one pair is complex conjugate. one mode will be overdamped and one underdamped. respectively. This is mixed damping. 2.4 Example 4 - Overdamped system Given cl=c3=3 kg/ssc and c2=1 kg/sec. the damping matrix will be [CJ 8 -1 4 These damping values when inputted into QUALADS will produce: N‘INNIII‘I‘I‘I‘INII‘IIN‘II‘I‘II‘ll}! The system is overdamped {Nihilliflfllilliflifl*flflliill The damped eigenvalues are: -2.618 xi =- -.382 . A2 A3=-i.o . it: -4.0 which agrees with QUALADS conclusion. Chapter 3 OUALADS SAHPLE RUN This chapter provides instruction on how to use QUALADS. First an illustration of a QUALADS run will be presented. flora detailed explanations of the subsections denoted by fi3.2.n* appear at the end of the chapter. as well as general comments on QUALADS (section 3.2.0). 3.1 - Sample Run Consider the two degree of freedom system with its four variants that is examined in chapter 2 (Figure 1). [II] = [C] = [x] = O 1 -c2 c2+c3 -l.5 2.5 l) c18c2=c3= 1 kg/sec 2) cl=c3= 2 kg/sec . c2= 1 kg/sec 3) c1=c38 2.5 kg/sec . c2= .5 kg/sec 4) c1=c3= 3 kg/sec . c2= 1 kg/sec The following QUALADS run determines the damping type and the damped eigenvalues for each of the four systems shown above. The underlined items are the information entered by the user. 0K. SE6 #QUALADS Enter the size of the system. 2 . ”3.2.1“ Enter mass matrix h by rows now 1 H( 19 I): H( 10 2): *3.2.2* KN“ ROW 2 N( 2. I)‘ H( 2. 2)8 F40 Enter stiffness matrix K by rows ROW 1 K( 10 1,8 2.5 K( 1. 2)8 -1.5 ROW 2 m 2. 1): 1,5 m 2. 21- Li Enter the damping matrix C by rows ROW 1 C‘ 1v 1). 2 c1 1. 2): :1 R09 2 C‘ 20 1’3 2L c< 2. 2): g INNNINNINNNNNNUNNINNUQNiii! The system is underdamped runuuanaainuanucuuauukauaxx Do you want to compute the damped eigenvalues? 1 13.2.3» Vould you like to open a file in which to place system matrices and output values? 1 *3.2.4* Enter filename EXQUI 10 The following information can be printed to your output file. 1)Nass hatrix 2)Stiffness hatrix 3)Damping Matrix 4)Critical Damping Natrix 5)Undamped Natural Frequencies 5)Uodal flatrix 7)Damped Eigenvalues Enter the option numbers desired in a single line separated by commas. 1. '4' .6. ’7 *3e2e5* Do you wish to run again with a different C matrix? 1 .3.2.5. Enter the damping matrix C by rows Row 1 c1 1. 1): 1 c1 1. 21: :L ROW 2 C‘ 20 I). :L c< 2. 2)- ; NINNINNINNNNNNNNNNNNNINNNNNNNNIII The system is critically damped arrr»nuaauunaaruauuunurnruxuuuuau Do you want to compute the damped eigenvalues? Y The following information can be printed to your output file. l)nase hatrix 2)Stiffness Natrix 3)Damping flatrix 4)Critical Damping flatrix 5)Undamped Natural Frequencies Slflodal flatrix 7)Damped Eigenvalues Enter the option numbers desired in a single line separated by commas. 3.7 11 Do you wish to run again with a different C matrix? 1 Enter the damping matrix C by R0" 1 (:1 1. 1): g C( 19 2): :fii ROW 2 C‘ 29 I): :fii C( 29 2), ; NINNINNNNNNNINNNNNN”Ullliiiiflfl The system has mixed damping NNNNINNINNNNNNNINNINNINNNNNNNN rows Do you want to compute the damped eigenvalues? Y _The following information can be printed to your output file. 1)flass hatrix 2)Stiffnees flatrix 3)Damping flatrix 4)Critical Damping hatrix 5)Undamped Natural Frequencies 6)flodal flatrix 7)Damped Eigenvalues Enter the option numbers desired in a single line separated by commas. 3.7 Do you wish to run again with a different C matrix? X. Enter the damping matrix C by rows ROW 1 c1 1. 1):- g C( 1. 2): :_1_ R09 2 c1 2. 1): ;1_ C( 20 2): fl l2 IflififlflflflfllilflfifliflflfliiNIH!!! The system is overdamped Niii!!!NNNINNNNNNNNNNNININ Do you want to compute the damped eigenvalues? Y The following information can be printed to your output file. l)flass Natrix 2)Stiffnees flatrix 3)Damping flatrix 4)Critical Damping flatrix 5)Undamped Natural Frequencies 5)flodal flatrix 7)Damped Eigenvalues Enter the option numbers desired in a single line separated by commas. 3.7 Do you wish to run again with a different C matrix? 1‘. ii! Output data file is EXDUT OK. In the QUALADS run above the user requested that an output file be opened and named EXDUT. EXDUT was supplied with the data requested by the user from the menu of options presented after each system solution. The file EXDUT would appear as follows. 13 NINNM‘IMMMI’NIIMI‘IIINM Underdamped system Iilmflfleiiiifleiiefllflfl MASS MATRIX Column Column I l) I 2) 0.100000D+01 0.0000000+00 0.000000D+00 0.100000D+01 STIFFNESS MATRIX\ Column Column ( l) ( 2) 0.2SOOOOD+01 -0.1500000+01 -0.1500000+01 0.2500000+01 CRITICAL DAMPING MATRIX Column Column ( 1) ( 2) 0.300000D+01 -0.100000D+01 -0.100000D+01 0.300000D+01 UNDAMPED NATURAL FREQUENCIES (rad/sec) 1 O.IOOOOD+OI 2 O.2000OD+OI EIGENVECTORS (MODAL MATRIX) Column Column ( 1) ( 2) 0.7071 -0.7071 0.7071 0.7071 DAMPING MATRIX Column Column ( 1) ( 2) 0.2000000+01 -0.1000000+01 -0.1000000+01 0.2000OOD+01 DANPED EIGENVALUES (Unpaired) -0.500000D+00+-0.86502SD+001 -0.5000000+00+ 0.855025D+00i -0.150000D+01+ 0.132288D+Oli -0.1500000+01+-0.1322880+011 bUMH 14 NNNNQNINMNMMQNN'Ifiiliflfiiifi Critically damped system NNNIINQNNNNNNNININNNNNNNNN DAMPING MATRIX Column Column I 1) I 2) 0.3000000+01 -0.1000000+01 -0.1000000+01 0.300000D+01 DAHPED EIGENVALUES (Unpaired) 1 -0.100000D+01+ 0.0000000+00i 2 -0.100000D+01+ 0.000000D+00i 3 -0.200000D+01+ 0.0000000+00i 4 -0.200000D+01+ 0.000000D+001 {MMNNMNNINMMINNNNMNNMNNNNII System with mixed damping ruaaruuaaaauuaaaurutuuuuaar DAMPING MATRIX Column Column I I) I 2) 0.3000000+01 -0.500000D+00 -0.5000000+00 0.300000D+01 DAMPED EIGENVALUES (Unpaired) 1 -0.175000D+01+-0.958246D+001 2 -0.175000D+01+ 0.96824SD+00i 3 -0.5000000+00+ 0.000000D+00i 4 -0.2000000+01+ 0.0000000+001 MNINNNQMNNNINNIMMNN Overdamped system iliihfllilieillilifle DAMPING MATRIX Column Column I 1) I 2) 0.400000D+01 -0.1000000+01 -0.100000D+01 0.4000000+01 15 DAMPED EIGENVALUES (Unpaired) -O.3819560+00+ 0.000000D+001 -O.IOOOOOD+OI+ 0.000000D+001 -O.261803D+01+ 0.000000D+OOi -O.4000OOD+01+ 0.000000D+001 bIJb3H The command file for the QUALADS run shown would be as follows. m 0 *QUALADS wow- ‘ m Ulm IBJM ll NI-CDOF-th p..- p.- o m- I U- Ccr Overdamping C = Ccr Critical Damping C < Ccr Underdamping Note that the qualitative method allows for simple and almost computation free determination of damping type if the damping constant C is changed. For the one degree of freedom system the amount of work required in the determination of the eigenvalues for the quantitative analysis is just slightly more than the work required in comparing the critical damping to the damping constant for the qualitative analysis. The difference in workload increases when a multidegree of freedom system is examined. 23 4.2 - Multidegree of Freedom System A multidegree of freedom system is described by a group of coupled second order differential equations which are expressed in matrix form below. (1111;) + ((211.21 + [1(1le = o The coefficient matrices are the mass matrix [M]. the damping matrix [C] and the stiffness matrix [K]. all of which have real elements. If the system has N degrees of freedom. then the square matrices [M]. [C] and [K] are size N x N. The unknowns are the terms of the Ix} vector where each x represents the relative motion of the system at some location. For an N degree of freedom system there will be N pairs of eigenvalues. Each pair forms a normal mode which looks like the solution to a one degree of freedom system. Thus each mode can be defined as overdamped. critically damped or underdamped. The system’s solution is a linear combination of these normal modes. where the magnitude of each mode’e contribution depends on the initial conditions. If all of the normal modes are overdamped or nonoscillatory. then the system must be nonoscillatory: it is called an overdamped system. If all of the modes are critically damped. then the system is critically damped. If all of the modes are 24 underdamped. then the system is underdamped (oscillatory). There is of course the possibility that the modes are a mixture of the three damping types. In this case the system is said to be mixed damped. As in the one degree of freedom system. 'the principal step in a quantitative analysis of a multidegree of freedom system is to calculate the system’s eigenvalues. Examining the eigenvalues will tell if the system can oscillate. The most comprehensive method to solve for the eigenvalues of a size N system is to convert the system into a size 2 x N system of first order state equations. The determination of the 2N eigenvalues requires lengthy matrix operations. Taking this into account. as well as the fact that system size has doubled. may make a quantitative analysis less attractive than a qualitative analysis in a design context. 4.; - ualigative Method by Inman and Andry A qualitative analysis of a multidegree of freedom system should determine whether a system can oscillate without computing the solution. i.e. eigenvalues. A method which accomplishes this was derived by Inman and Andry [1.2]. Their algorithm defines and uses a matrix called the critical damping matrix. [Ccr]. [Ccr] = 2[Ml[Ul[wl[UT][Ml 25 The matrix [U] is the undamped modal matrix normalized such that [UTUIMJIUI = [I] and the matrix [w] is the diagonal matrix containing the undamped natural frequencies. Both matrices come from solving the eigenvalue problem of the undamped system [lexl + [KJle = o. The algorithm is restricted to asymptotically stable systems with symmetric. positive definite mass and stiffness matrices and a symmetric. positive semi-definite damping matrix. The critical damping matrix serves a function similar to that of the critical damping constant in the single degree of freedom system. The damping type of a system is ascertained by determining the definiteness of the matrix equal to the difference between the damping matrix [C] and the critical damping matrix [Ccr]. The relationships are as follows. - or System damping gype positive definite overdamping zero critical damping negative definite underdamping indefinite mixed damping Note again that changing the damping of the system does not affect the critical damping value. Thus to reexamine a system with a different damping matrix only requires 26 determining the definiteness of the new [C - Ccr]. This makes the qualitative method more attractive than the quantitative. which would require that the analysis be reworked from the beginning. Chapter 5 QUALADS COMPUTATIONAL INFORMATION This chapter presents a brief description of the major methods of computation used in QUALADS. It then compares the computational load for the qualitative method of determining damping type with that for the quantitative procedure of determining the damped eigenvalues. 5.; - Methods of Compptatign nvalu r l m - The undamped natural frequencies. modal matrix. and damped eigenvalues are all computed using the IMSL eigenvalue package EIGZF (3]. Given the mass and stiffness matrices. EIGZF computes the undamped eigenvalues and eigenvectors. which produce the natural frequencies and normalized modal matrix. For the damped eigenvalues the system equations are converted to state space equations [4]. The ”mass" and ”stiffness” matrices for state space are as follows. I 0 I M I I -M I 0 I ”M” a I----I----I "K" = I----I----: I M I C I I 0 I K I Given the above matrices EIGZF computes the damped eigenvalues unpaired and unordered. 27 28 Qefiniteness - To determine the damping type QUALADS must determine if the [C - Ccr] matrix is positive definite. negative definite. indefinite or zero. The procedure used involves determining the sign of the determinants of the aubmatrices. Let [A] = [C -Cch. For the square matrix [A] that is size n. the square submatrix tAkl is a size k matrix formed by the first k rows and first k columns of the [A] matrix. The following conditions determine the corresponding definite type [5]. on 'ti n Iype DetItAkl) > 0 for k=1.2.....n positive definite DetII-Akl) > 0 for k-1.2.....n negative definite If [A] is neither positive definite nor negative definite nor zero. then it is indefinite. tot' ta ' it - If the damping matrix [C] of a system is positive semi-definite. then the system could be asymptotically unstable. QUALADS uses the criterion that a system is asymptotically unstable if any of the modal vectors of the undamped system lie in the null space of the damping matrix [C] [61. This check is achieved by examining the matrix that is the product of the damping matrix (Cl and the modal matrix [U]. If none of the columns of this product matrix are zero vectors then the system is stable. 29 5.2 - Computation Time The major part of the computational load in determining the damping type occurs during the computation of the critical damping matrix [Ccr] as opposed to determining the definiteness of the [C - Ccr] matrix. If a set of several damping matrices is examined during one run then the computational load is considerably less after the initial system has been solved. This occurs because the critical damping matrix. which is determined by the mass and stiffness matrix. has already been computed. The computational load incurred by the quantitative analysis. which consists of computing the damped eigenvalues. is much greater than that of the qualitative analysis. Figure 3 shows the CPU (Central Processor Unit) time required to determine [Ccr]. the definiteness (damping type) given [Cl and [Ccr]. and the damped eigenvalues. All the calculations were for overdamped systems. which produces the maximum computational effort in the definiteness routine. Although the data suggests that CPU time is a function of the system size squared for all three computations. the controlling coefficient of the definiteness algorithm is small relative to the coefficient of the eigenvalue algorithm. Therefore. as system size increases the computational savings afforded by the qualitative method become substantial. 60 40 C P U T 20 I M E (sec) 0 A B c o 30 SYSTEM SIZE - Computation of critical damping matrix [Ccr] - Computation of definiteness of [C] - [Ccr] (damping type) given [Ccr] - Computation of damped eigenvalues - Computation of damping type given [M]. [C]. and [K] (i.e. A + 8) Figure 3. CPU Time Plot These results were obtained on a Prime 750 virtual memory super mini-computer with a 6 Mb memory operating under PRIMOS Revision 19.2. Chapter 6 SUMMARY A qualitative analysis of damped linear dynamic systems was considered. This analysis consists of an algorithm derived by Inman and Andry that determines a matrix called the critical damping matrix. which is then used to generate the damping type of a system described by a set of linear second order differential equations. Software called QUALADS was written to implement this qualitative analysis method. In addition to finding the damping type. QUALADS can execute a quantitative analysis of a system by computing the damped eigenvalues. 6.1 - Conclusions In the examination of a damped system the question of the system's ability to exhibit oscillatory behavior in free motion is generally an important question. The qualitative analysis described in this thesis offers an effective method for answering this question. As compared to a quantitative solution to the problem. the qualitative analysis requires less computational effort. The compute time and effort difference is even greater when several damping mechanisms for the same model are to be examined. There the qualitative 31 32 solution determined by QUALADS may yield sufficient information about a system or provide sufficient preliminary data to guide a further investigation of the damping behavior. 6.2 - Extensions An interesting extension to the qualitative analyses examined in this thesis is the inverse problem of designing a system that will not oscillate when perturbed. For an overdamped or nonoscillatory system. the analysis requires that the matrix [ C - Ccr ] be positive definite. This requirement would give rise to a set of nonlinear inequalities. An algorithm to solve these inequalities. along with additional user imposed constraints. could prove to be a very useful design tool. L I ST OF REFERENCES 33 LIST OF REFERENCES Inman. D. J. and Andry. A. N.. ”Some Results on the Nature of Eigenvalues of Discrete Damped Linear Systems". Journal of Applied Mechanics. Dec. 1980. Inman. D. J.. ”Nature of the Solutions of Damped Linear Dynamic Systems. Ph.D. Dissertation. Michigan State University. 1980. ”IMSL Library - User’s Manual”. IMSL Inc.. Texas. 1964. Chapter E. Meirovitch. L.. "Analytical Methods in Vibrations". The Macmillan Co.. N.Y.. 1967. pg. 411.412. Phillips. 3.. "Lecture Notes on Matrix Theory”. Michigan State University. 1981. Moran. T. J.. "A Simple Alternative to the Routh-Hurwitz Criterion for Symmetric Systems”. Journal of Applied Mechanics. Dec. 1970. APPENDIX 34 The following appendix consists of the source code for QUALADS. The code is written in FORTRAN 77. The first page of code is the program MAIN which calls the main subroutines of the program. An explanation of the purpose of each subroutine is given before its call. A calling tree for the subroutines of QUALADS follows. MAIN (QUALADS) DIMSIZ MATINP CHKM CRKSYM PDCRK DETER DIAGPR MATINP CHKK CNKSYM PDCHK DETER DIAGPR MATINP CRKC CHKSYM PDCRK DETER DIAGPR PSDCRK DETER DIAGPR EIGEN EIGZF (IMSL) SORT NORM CMKSTB COMP DEFNT DETER DIAGPR EIGQUE DMPEIG EIGZF (IMSL) OPNFIL OUTFIL RERUN 35 The subroutines shown in the calling tree appear in the code in the following order. The subroutines TRMTYP and PAGE perform screen paging for those terminals that require it and are called throughout QUALADS in several of its subroutines. MAIN DIMSIZ CRKM CRKK CHKC CHKSTB EIGEN SORT NORM DEFNT EIGQUE DMPEIG OPNFIL OUTFIL RERUN CRKSYM PDCRK PSDCHK DETER DIAGPR TRMTYP PAGE 36 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CMNNNINNNINMNNNNNNNIMINNIIIINNINNNIRIMININIIIIINIININIININNNIIIIMNNC C! PROGRAM MAIN MC C! !C C! MAIN calls the major subroutines and creates and !C C! dimensions the arrays they will need. !C C! ' !C CIMIINIIIMIN“III.MIMMI’IMIIQINIIINIIIMIMMIIIIMII'MM'MIIIIIMIIIMIIIIIC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C! !C C! !Variable List!. !C C! C - Damping matrix !C C! CCR - Critical damping matrix !C C! K -Stiffness matrix !C C! LAMDA - Damped eigenvalues !C C! M - Mass matrix !C C! U - Modal matrix (undamped eigenvectors) !C C! V - Undamped natural frequencies !C C! UK!! - General work matrices !C C! Logicals - .TRUE. meaning !C C! ABORT - M or K are illegal I i.e. not symmetric !C C! or positive definite) !C C! CONT - User has opted to input another C matrix !C C! DEVAL - User asked to compute damped eigenvalues !C C! FIRST - Analysis Is of first C matrix !C C! ON - C matrix is legal (i.e. symmetric and at !C C! least positive semi-definite) !C C! OPEN - User has opened a file for output values !C C! PSD - C matrix is positive semi-definite !C C! UCALC - The modal matrix U and the undamped !C C! eigenvalues have been computed !C C! IC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C REAL!8 MI4000O).CI4000O).KI4000O).UI4000O).CCRI4000O).UI200) COMPLEX'IB LAMDAI400) CHARACTERIS TYPE REAL'B UK33I4000O).UK44I4000O) REAL'B VKII(160000).UK22I160000).UKLI32000O).UKII400) COMPLEX!16 UKICIZOO).UX11C(160000) CHARACTER!32 FILEN LOGICAL FIRST.OK.OPEN.DEVAL.UCALC.ABORT.CONT.PSD COMMON/STUFF/M.C.X.U.CCR.U.UK33.UK44.UK11.9K22.UK1.UK1C COMMON/STUFZ/UKL.UK11C C C UCALC=.FALSE. FIRST=.TRUE. ABORT=.FALSE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CIIIIIIINIlMINOR!NNNMMINNIINMNNNNMNINNINIIIINMINMINONINRNNMINIINNNNC C! DIMSIZ asks the user for the system size (N). i.e. the !C C! dimensions of the matrices. !C CRINIIIIIMIN!MINIMININNNNIINNNINNNNNINNIININNINNNNNNINNIMIIMNNNNNNNC C C CALL DIMSIZIN) 37 NNBNIN N2=2IN NN=2!N2!N2 C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Cullrrunrrurnnnrnrerunno:neurnnuurlnulmunelnurluuour!surruunuurururc C! MATINP asks the user to input the mass. stiffness. and !C C! damping matrices (M.N.C). !C C! CRKM and CHKK check the M and K matrices. respectively. for !C C! symmetry and positive definiteness. !C C! CRKC checkc the C matrix for symmetry and positive semi- !C C! definiteness. !C C! CRKSTB checks to see if the system is asymptotically stable. !C C! (Since this check requires the modal matrix U. EIGEN is !C C! is called to compute it. See EIGEN below.) !C curluonlureuinureerroneous:neuronalnolunrmurlluuuunslrnunravuunuannc C CALL MATINPIM.N.’M') CALL CHKMIM.N.UKL.NN.ABORT) IFIABORTlGO TO 999 CALL MATINPIK.N.'K') CALL CRKKIK.N.UKL.NN.ABORT) IFIABORTlcO TO 999 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 10 CONTINUE CALL MATINPIC.N.’C') C C CALL CHKCIC.N.UKL.NN.OK.PSD) IFIOKITREN IFIPSDlTHEN IF(.NOT.UCALC)TREN CALL EIGENIM.K.U.U.N.NN2.UK1.UKL.UK1C.UK11.9K22.UK11C) UCALCI.TRUE. ENDIF CALL CRKSTBIC.U.UK11.N.OK) ENDIF ENDIF C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IFIOxlTNEN C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Crlonurnnnrlunnu!marina!nurrwrnruuuunluaururunner»:nirurlunruruluuuc C! EIGEN computes the undamped eigenvalues (natural !C C! frequencies)(v) and eigen vectors (modal matrix)(U). !C C! NN2 is the dimension of a workspace needed by the !C C! eigenvalue package EIGZF (IMSLD). !C CNNIINIIINNINNIN!!!NflelllellulmmeemlelllIINIMIIINQNNNNNIllllmlllllic C C IF(.NOT.UCALC)THEN CALL EIGENIM.K.U.U.N.NN2.UK1.UKL.UK1C.UK11.UK22.UK11C) UCALC8.TRUE. 38 ENDIF C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Cnnmurrkurulnruulruin.oururn!un!!!!rrirmlruuuunrlrrrunnrunununrurnlc C! COMP computes the critical damping matrix ICCR). !C CMNINNIINllllliilNNNNNNIINIMINNMNINNNNMINNIIIIIINNMIIININNIIINNNIMIC C C IFIFIRSTlTHEN CALL COMPIU.U.CCR.M.N.UK11.9K22.VK33.UK44) ENDIF C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CNMlIIINMIlIIMIMIIMNNNIINNMNMNMININIMIRIININNNINIMIINNNMNININNMIIIMC C! DEFNT computes the definiteness (TYPE) of the matrix which !C C! is the difference of the damping matrix (C) and the !C C! critical damping matrix ICCR). !C C! NN is the dimension of a vector formed by the !C C! aforementioned difference (C - CCR). !C C! The system description is printed out by DEFNT at this time.!C C! i.e. it is stated which type of system (overdamped. !C C! critically damped. underdamped. or mixed) is present. !C curnruururranour:re!!!nu:runnunrnnnnnuuuuuunruunupure...nunnnurrruuc C C CALL DEFNTIC.CCR.TYPE.N.NN.UK11.UKL) C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CINNIINNNININIININNNINMRINNI!!!IINIININNIINNNNNIINNNINNNIIIIIINNNNNC C! EIGGUE Asks the user if he wnats to compute the damped !C C! eigenvalues. DEVAL is true if answer is yes. !C C! DMPEIG computes them. !C CIINNNNRNMilli!!!NIINNIIIINNNINMNNIIDNNNINNNNINMIIIIINNIMIN!!!ICNNIC C C C C CALL EIGQUEIDEVAL) IFIDEVAL)TREN CALL DMPEIGIM.C.K.LAMDA.N.N2.NN2.VK1.UKL.UK11.9K22.UK11C) ENDIF C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CINNIINNINNMININIIMNINMNMNINNNMINNNIMill!MillIIINNINNNMNNIMIIINNNINC C! OUTFIL creates an output file with the users choice of !C C! values. such as any of the system matrices. !C CINNNNIINMMIINNNIINNNNIINIIMNNIMNMNIIININIMINIINNNNNNINIINIMMMNNMNNC C C IFIFIRST1CALL OPNFILIOPEN.FILEN) IFIOPEN1CALL OUTFILIM.C.K.CCR.U.U.N.TYPE.LAMDA.N2.DEVAL) C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FIRST'.FALSE. ENDIF C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CIMNNNMIMNNNNNMIMNNNINIMMNMNINMNNIMIININIIINNNNNNNNNIIIIMNNNNIIIIIIC C! RERUN asks the user if he wishes to rerun the problem !C 39 C! using a new damping matrix. but maintaining the !C C! same mass and stiffness matrices. !C CNINNNNNIIIIINNNNMINNINNNNNIINIIIIIMNNNININIIIIIMIIIIMINNNINNONINIIC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CALL RERUNICONT) IFICONT)GO TO 10 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IFIOPEN)THEN CLOSEI12) VRITEI1.100)FILEN 100 FORMATII.’!!! Output data file is ’.A32./) ENDIF 999 CALL EXIT END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4O CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine DIMSIZ C C C C DIMSIZ asks the user to input the dimension size N of the C C system. i.e. matrix size or number of degrees of freedom. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !Variable List! N - System size LOGICALS - .TRUE. meaning ERROR - Unrecognized character REENT - Variable needs to be reenterd CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE DIMSIZIN) 0 0000000 0 0000000 CHARACTER!32 STRING LOGICAL ERROR. REENT CALL TRMTYPINP) IFINP.EQ.1)CALL PAGEINL) REENT=.TRUE. UNILEIREENT)DO ERROR-.TRUE. PRINT!.' Enter the size of the system.’ READ(!.’(A32)’)STRING READISTRING.100.ERR-IO)N REENT=.FALSE. ERROR=.FALSE. 10 IFIERROR)THEN PRINT!.'Unrecognized character' REENT-.TRUE. ENDIF IFIN.GT.200.AND..NOT.ERROR)TREN PRINT!.'Size limit is 200’ REENT'.TRUE. ENDIF ENDURILE 100 FORMATII3) C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 41 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine MATINP C C C C MATINP asks the user to input the values of the mass. C C stiffness and damping matrices. entry by entry and row by row. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C !Variable List! C C A - System matrix (M.K. or C) to be entered C C N - System size C C MAT - Name of matrix to be entered IM.K. or C) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE MATINPIA.N.MAT) C C REALIB AIN.N) CRARACTERl32 STRING CNARACTERII MAT C C CALL TRMTYPINP) C C IFINP.EQ.1)CALL PAGEINL) PRINT!.’ ’ PRINT!.' ’ IFIMAT.EQ.'M’)TREN PRINT!.'Enter mass matrix M by rows' ELSEIFIMAT.EQ.’K’)TREN PRINT!.’Enter stiffness matrix K by rows’ ELSEIFIMAT.EO.’C’)TNEN PRINT!.’Enter the damping matrix C by rows’ ENDIF NL=3 DO 20 I-1.N IFINL.GE.30.AND.NP.EQ.1.AND.J.NE.1)CALL PAGEINL) PRINT!.’ ’ URITEI1.60)I 80 FORMATI’ROU ’.I2) NL'NL+2 DO 30 J81.N IFINL.GE.30.AND.NP.EQ.llCALL PAGEINL) URITEI1.40)MAT.I.J 40 FORMATIAI.’('.I2.’.'.I2.’)= '.$) 50 READ(!.’IA32)')STRING READISTRING.!.ERR=100)AII.J) NLINL+2 GO TO 30 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Error statment. C C C 100 PRINT!.’Unrecognized entry. Rs-enter value.’ NLINL+2 GO TO 50 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 30 CONTINUE 20 CONTINUE 42 RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine CRKM C C C C CRKM determines if the M matrix is legal. i.e. symmetric C C and positive definite. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !Variable List! M - Mass matrix C UK - Vork matrix C LOGICALS - .TRUE. meaning C C C C 0 ABORT - M illegal OK - M is legal . CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE CNKMIM.N.UK.NN.ABORT) 0 (10I1CIOITCIO REAL!8 MIN.N).UKINN) LOGICAL OK.ABORT CALL CRKSYMIM.N.OK) IFIOK)TREN CALL PDCRKIM.N.UK.NN.OK) IFI.NOT.OK)THEN PRINT!.’ M is not positive definite. Program abort.’ ENDIF ELSE PRINT!.’ M is not symmetric. Program abort.’ ENDIF IFI.NOT.OK)ABORT=.TRUE. RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 43 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine CHKK C C C C CHKK determines if the K matrix is a legal sytem matrix. C C i.e. symmetric and positive definite. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C !Variable List! C K - Stiffness matrix C C UK - Uork matrix C C LOGICALS - .TRUE. meaning C C ABORT - K not legal C C C C C C C 0 OK - K legal CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE CHKKIK.N.UKL.NN.ABORT) REAL!8 KIN.N).UKL(NN) LOGICAL OK.ABORT CALL CRKSYMIK.N.OK) IFIOK)TREN CALL PDCHKIK.N.VKL.NN.OK) IFI.NOT.OK)THEN PRINT!.’ K is not positive definite. Program abort.’ ENDIF ELSE PRINT!.’ K is not symmetric. Program abort.’ ENDIF IFI.NOT.OK)ABORT-.TRUE. RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 44 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine CNKC C C C C CHKC determines if the C matrix is legal. meaning it C C must be symmetric and either positive definite or positive C C semi-definite. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C !Varibale List! C C C - Damping matrix C C UK - Dork matrix C C LOGICALS - .TRUE. meaning C C OK - C matrix is legal C C PD - C matrix is positive definite C C PSD - Cmatrix is positive semi definite C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE CRKCIC.N.VK.NN.OK.PSD) C C REALuB C(N.N).UK(NN) LOGICAL OK. PSD. PD C C OK=.TRUE. PSD-.FALSE. C C CALL CRKSYMIC.N.OK) IFIOKlTREN CALL PDCRKIC.N.UK.NN.PD) IF(.NOT.PD)TREN CALL PSDCRKIC.N.VK.NN.P5D) IF(.NOT.PSD)TREN OK=.FALSE. PRINT!.’ C is not positive semi-definite.’ PRINT!.’ Illegal system.’ ENDIF ENDIF ELSE PRINT!.’ C is not symmetric. Illegal system.’ ENDIF C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 45 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 00000 Subroutine CHKSTB CRKSTB determines if the system is asymptotically stable. This check is only necessary if the C matrix is positive semi-definite. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 0 00000000 0 C C - Damping matrix CU - Product of the C matrix and U matrix U - Modal matrix. undamped eigen values LOGICALS - .TRUE. meaning STBL - System is asymptotically stable C C C C C C !Variable List! C C C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE CRKSTBIC.U.CU.N.STBL) REAL!8 C(N.N).UIN.N).CU(N.N) LOGICAL STBL C STBL=.TRUE. C CALL DMMLTICU.C.U.N.N.N) J-l URILEISTBL.AND.J.LE.N)DO l-i URILEIDABSICUII.J)).LT.1.E-6.AND.I.LE.N)DO I-I+l ENDURILE IFII.EO.N+1)STBL-.FALSE. J=J+1 ENDURILE C IFI.NOT.STBL)TREN PRINT!.’ System is not asymptotically stable.’ PRINT!.’ Illegal system.’ ENDIF C RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 45 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C SUBROUTINE EIGEN C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C EIGEN computes the eigenvectors and natural frequencies of C C the system given the mass and stiffness matrices. C C All numbers are in double precision. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C SUBROUTINE EIGENIMMAIN.KMAIN.V.UREAL.N.NN.BETA.UK.ALFA &.MASS.STIF.UCOMP) C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C REAL!8 MMAININ.N).KMAININ.N).VIN) REAL!8 BETAIN).UKINN).MASSIN.N).STIFIN.N).UREALIN.N) COMPLEX!16 ALFAIN).UCOMPIN.N) C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Rename the mass and stiffness matrices for preservation. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C DO 110 I-I.N DO 100 J81.N MASSII.J)=MMAINII.J) STIFII.J)-KMAIN(I.J) 100 CONTINUE 110 CONTINUE C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Compute the eigenvectors and eigenvalues C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C CALL EIGZFISTIF.N.MASS.N.N.2.ALFA.BETA.UCOMP.N.UK.IER) C C DO 210 I-I.N DO 200 J-i.N UREALII.J)-REALIUCOMP(I.J)) 200 CONTINUE 210 CONTINUE C C DO 220 IBI.N ALFAII)=ALFAII)/BETA(I) UIIl‘ALFAII) IF (VII) .GT. 1.D-8) THEN 47 UII)=DSQRTIUII)) ELSEIFIUII).GT.-1.D-5.AND.V(I).LT.1.D-6)THEN VII) 8 0.0 ELSE PRINT !.’THIS CAN NEVER HAPPEN -JANE HAUNINS’ PRINT I. ’NEGATIVE ARG TO SQRT IN EIGEN' END IF 220 CONTINUE C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SORT orders the modal matrix so that the first column C C is associated with the first eigenvector. etc. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CALL SORTIV.UREAL.N) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C NORM normalizes the modal matrix U so that UT!M!U-I. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CALL NORMIUREAL.MMAIN.N.MASS.STIF.VK) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 48 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C SUBROUTINE SORTIV.E.M) REALIB UIM).E(M.M) MMSM-1 DO 40 I=1.MM DO 30 J=I.MM IFIVIJ).LE.U(J+1)) GO TO 25 N=J+1 UTEMP‘UIJ) UIJISUIN) UIlevTEMP DO 20 K=I.M ETEMP-EIK.J) EIK.J)IE(K.N) EIK.N)8ETEMP CONTINUE CONTINUE CONTINUE CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 20 10 C C SUBROUTINE NORMIU.M.N.UT.UTM.UTMU) C REALIB UIN.N).M(N.N) REAL!8 UTIN.N).UTMIN.N).UTMU(N.N) CALL DMTRNIUT.U.N) CALL DMMLTIUTM.UT.M.N.N.N) CALL DMMLTIUTMU.UTM.U.N.N.N) DO 10 J81.N DO 20 I=I.N UII.J)=U(I.J)/DSQRTIUTMUIJ.J)) CONTINUE CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 49 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine COMP C C C C COMP computes the critical damping matrix. using the C C formula: Ccr=2!M!UT!VM!U!M . C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C !Variable list! C C C C CCR - Critical damping matrix. C C N - Size of matrices. C C U - Vector of natural frequencies (ordered). C C M - Mass matrix. C C UT - Transpose of modal matrix U. C C VM - Diagonal matrix of natural frequencies. C C U - Modal matrix (UT!M!U-I). C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE COMPIU.U.CCR.M.N.UT.VM.C1.C2) C REAL!8 UIN.N).VIN).MIN.N) REAL!8 UTIN.N).UM(N.N).CCRIN.N) REALIB CIIN.N).C2(N.N) DO 200 J=I.N DO 210 I=1.N VMII.J)=O. UMII.I)=U(I) 210 CONTINUE 200 CONTINUE CALL DMTRNIUT.U.N) CALL DMMLTIC1.M.U.N.N.N) CALL DMMLTIC2.C1.UM.N.N.N) CALL DMMLTIC1.C2.UT.N.N.N) CALL DMMLTIC2oC1.M.N.N.N) CALL DMSCLICCR.C2.N.N..2DOI) RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 50 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Subroutine DEFNT 0 DEFNT determines the type of damping present in the system. 1)C-Ccr=0: critically damped 2)C-Ccr positive definite: overdamped 3)C-Ccr negative definite: underdamped 4lC-Ccr indefinite: mixed damping The definiteness of the matrix is determined by examining the determinant’s of the submatrices: sizes I X I through N X N. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C C C C C C C C C !Varimble List! C C - Damping matrix. C CCR - Critical damping matrix. C DIF - Matrix that is the difference of C - CCR C TYPE - Type of definiteness. C N - Size of matrices. C DET - Determinant of aubmatrices. C MS - Size of submatrix. C DPR - Product of diagonal terms of submatrix. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE DEFNTIC.CCR.TYPE.N.NN.DIF.TM) 0 00000000000000000000000000 REAL'B CIN.N).CCRIN.N).DIFIN.N).TMINN).DET.DPR CRARACTER'B TYPE LOGICAL QUITP.QUITN C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 0 PRINT!.’ ’ C ---Compute C - Ccr. CALL DMSUBIDIF.C.CCR.N.N) 0 000 TYPE8’INDEF’ QUITP=.FALSE. QUITN'.FALSE. ---Check to see if C - Ccr = 0.. thus critically damped. 000 NC'O DO 10 1.19" DO 20 J31.N IF(DABS(DIF(I.J)).LE..OOOI)NC=NC+1 2O CONTINUE 10 CONTINUE IFINC.EQ.N!!2)TMEN TYPE"ZERO’ QUITP'.TRUE. QUITN'.TRUE. ENDIF 51 C --- Determine if C - Ccr is positive definite. negative definite C --- or indefinite. C MS=1 UHILEIMS.LE.N.AND..NOT.QUITP)DO CALL DETERIDET.DIF.TM.MS.N.NN.1.DO) CALL DIAGPRIDIF.N.MS.DPR) IFIDET/DPR.LT.1.E-I2)QUITP=.TRUE. MS=MS+1 ENDUHILE IFI.NOT.QUITP)TYPE=’POSDEF’ MS'I UHILEIMS.LE.N.AND.QUITP.AND..NOT.QUITN)DO CALL DETERIDET.DIF.TM.MS.N.NN.-1.DO) CALL DIAGPRIDIF.N.MS.DPR) IFIDET/DPR.LT.1.E-12lQUITN=.TRUE. MSIMS+1 ENDURILE IFI.NOT.QUITN.AND.QUITP)TYPE=’NEGDEF’ C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Print out the determined type C C C CALL TRMTYPINP) IFINP.EQ.1)CALL PAGEINL) C C IFITYPE.EQ.’POSDEF’)VRITEI1.100) 100 FORMATI//.26(’!’)./.’ The system is overdamped’./.26I’!’)) IFITYPE.EQ.’NEGDEF’)VRITEII.110) 110 FORMAT(//.27(’!’)./.’ The system is underdamped’./.27I’!’)) IFITYPE.EO.’INDEF’)VRITE(1.I20) 120 FORMATI//.30I’!’)./.’ The system has mixed damping’./. a 30(’!’)) IFITYPE.EO.’ZERO’)URITEII.130) 130 FORMATI//.33(’!’)./.’ The system is critically damped’./. & 33(’!’)) C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C RETURN END C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 52 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine EIGQUE C C C C EIGQUE asks the user if he wishes to compute the damped C C eigenvalues. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE EIGQUEIDEVAL) C C C C CHARACTER!1 ANSI3) LOGICAL DEVAL. UNREC C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C UNREC=.TRUE. PRINT!.' ’ PRINT!.’Do you want to compute the damped eigenvalues?’ VNILEIUNREC)DO READII.’(3AI)’)ANS CALL UPCASEIANS) IFIANSII).EQ.’Y’)THEN DEVAL=.TRUE. UNREC=.FALSE. ELSEIFIANSII).EO.’N’)TREN DEVAL=.FALSE. UNREC=.FALSE. ELSE PRINT!.’Unrecognizmd character. Enter yes or no.’ UNRECI.TRUE. ENDIF ENDUHILE RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Hi)“ '2'! J ' I, 1‘ 53 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine DMPEIG C C C C The subroutine DMPEIG calculates the damped eigenvalues. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C !Variab1e List! C C N - System or matrix size C C N2 - Two times the size N C C NN - workspace size of 2!N2!N2 C C M - Mass matrix as inputted C C K - Stiffness matrix as inputted C C C - Damping matrix as inputted C C LAMDA - Damped eigenvalues vector C C BETA - Output of EIGZF used to determine eigenvalues C C VK - workspace needed by EIGZF C C DM - Damped mass matrix. see below C C DK - Damped stiffness matrix. see below C C C C _ _ _ _ C C I I I I I I C C I 0 I M I I -M I 0 I C C DM = I----I----I DK 8 I----I----I C C I M I C I I 0 I K I C C I_ I _I I_ I _I C C C C C C C C UCOMP - Damped eigenvectors (not used) C C EIGZF - Library program that computes eigenvalues C C and eigen vectors C C C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE DMPEIGIM.C.K.LAMDA.N.N2.NN.BETA.UK.DM.DK.UCOMP) C C C C REALIB MIN.N).C(N.N).KIN.N).BETAIN2).VKINN) REAL!8 DMIN2.N2).DKIN2.N2) COMPLEX!18 LAMDAIN2).UCOMPIN2.N2).LAMTMP C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C The following loops construct the 2N by 2N mass and C C stiffness matrices. DM and OK as shown above. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DO 20 1'1.N DO 10 J81.N II=N+I JJ=N+J DMII.J)=O. DMII.JJ)=MII.J) DMIII.J)=MII.J) DMIII.JJ)8C(I.J) 54 DKII.J)=-M(I.J) DKII.JJ)=O. DK(IIIJ)=O‘- DKIII.JJ)=K(I.J) 10 CONTINUE 20 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Send the DM and DK matrices to EIGZF which computes the C C values LAMDA and BETA which define the eigenvalues. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CALL EIGZFIDK.N2.DM.N2.N2.2.LAMDA.BETA.UCOMP.N2.VK.IER) C C DO 40 I =1.N2 IFIDABSIBETAII)).LT.1.D-10)THEN LAMDAIIl=O ELSE LAMDAII)=LAMDA(I)/BETAII) ENDIF 40 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Separate the real from the cdmplex eigenvalues. C C the real going to the work array UK and the C C complex reforming the vector LAMDA C C NR is the number of real eigenvalues C C NC is the number of complex eigenvalues C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C NR=0 NC=0 DO 100 I81.N2 IFIDABSIDIMAGILAMDAII))).LT.I1.E-30!DREALILAMDA(I))))TREN NR-NR+1 VKINR)-DREALILAMDA(I)) ELSE NC=NC+1 LAMDAINC)=LAMDAII) ENDIF 100 CONTINUE IFINR.EQ.0)GD TO 140 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Order the real values from smallest to largest C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DO 120 IISI.NR-1 DO 110 I=1.NR-1 IFIUKII).LE.VKII+1))GO TO 110 UKTMP=UKIII UKII1=UKII+II UKII+1>=UKTMP 55 110 CONTINUE 120 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C If all the values are real put them back into the C C LAMDA vector and quit C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IFINC.EQ.0)TREN DO 130 [81.N2 LAMDAII)=DCMPLXIUK(I)) 130 CONTINUE ENDIF C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Order the complex values from smallest magnitude to C C largest by complex pair C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 140 D0 160 II=1.NC-1 ' DO 150 I31.NC-1 IFICDABSILAMDAII)).LE.CDABS(LAMDA(I+1)l)GO TO 150 LAMTMPSLAMDAII) LAMDAII)=LAMDA(I+1) LAMDAII+118LAMTMP 150 CONTINUE 150 CONTINUE C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C If all the values are complex then quit C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IFINC.EQ.N2)GO TO 180 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Place the real values back into the LAMDA vector after C C the complex. but rewrite them as complex numbers C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DO 170 I=NC+1.N2 [Isl-NC LAMDAII)=DCMPLXIUK(II)) 170 CONTINUE 180 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Since EIGZF outputs the eigenvalues as positive values. C C change them to the more common negative. C C C 56 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DO 200 I81.N2 LAMDAIIII-LAMDAII) 200 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 57 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine OPNFIL C C C C OPNFIL open a file for output values. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C !Variable List! C C LOGICALS - .TRUE. meanings C C OPEN - file opened for output values C C EX - filename already exists C C REENT - need to reenter filename C C UNREC - unrecognized character (i.e. not yes or no) C C LEGAL - legal filename C C FILEN - filename C C ANSI.ANS2 - yes or no answers C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE OPNFILIOPEN.FILEN) C LOGICAL OPEN.EX.REENT.UNREC.LEGAL CRARACTER!1 ANSII3).ANS2I3) CRARACTER!32 FILEN C C REENT=.TRUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PRINT!.’ ’ PRINT!.’Vould you like to open a file in which to place’ _ PRINT!.’ system matrices and output values?’ 5 READII.’(3A1)’)ANSI CALL UPCASEIANSI) IFIANSIII).EQ.’Y’)THEN VHILEIREENT)DO PRINT!.’Enter filenams’ READII.’IA32)’)FILEN C C LEGAL=.TRUE. IFIFILENII:1).GT.’0’.AND.FILEN(1:1).LT.’9’)TREN PRINT!.’IIIsga1 filename’ LEGAL-.FALSE. ENDIF C C UNREC=.TRUE. IFILEGAL)TREN INQUIREIFILE-FILEN.EXISTBEX) IFIEX)TREN PRINT!.’File already exists. Do you want to overwrite it?’ VHILEIUNREC)DO READII.’(3A1)’)ANS2 CALL UPCASEIANS2) IFIANS2I1).EQ.’Y’)TNEN REENT=.FALSE. UNREC=.FALSE. ELSEIFIANS2II).EQ.’N’)TREN REENT=.TRUE. UNREC=.FALSE. ELSE erw (Ta-“4' 58 PRINT!.’Unrecognized character. Enter yes or no.’ UNREC=.TRUE. ENDIF ENDVRILE ELSE REENT=.FALSE. ENDIF ENDIF ENDVRILE OPENI12.FILE=FILEN) OPEN=.TRUE. ELSEIFIANSIII).EQ.’N’)THEN OPEN=.FALSE. ELSE PRINT!.’Unrecognized character. Enter yes or no.’ GO TO 5 ENDIF RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 59 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine OUTFIL C C C C OUTFIL presents the user with the menu of options that C C can be printed to the output file. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE OUTFILIM.C.K.CCR.U.V.N.TYPE.LAMDA.N2.DEVAL) C C C C REALIB MIN.N).C(N.N).KIN.N).CCR(N.N).U(N.N).UIN) COMPLEXIIB LAMDAINZ) CHARACTER!1 RESPIIS) CHARACTER56 TYPE INTEGER COLI8) CRARACTER'32 TEXT LOGICAL DEVAL.CONT C C CONT-.TRUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C VRITEII2.’I//)’) IFITYPE.EO.’ZERO’)TNEN URITEI12.10) URITEI12.!) ’Critically damped system’ VRITEI12.10) ELSEIFITYPE.EQ.’NEGDEF’)TREN VRITEI12.20) VRITEI12.!) ’Underdamped system’ VRITEII2.20) ELSEIFITYPE.EQ.’INDEF’ITREN URITEI12.30) URITEI12.!) ’System with mixed damping’ URITEI12.30) ELSEIFITYPE.EO.’POSDEF’)TREN VRITEI12.40) URITEI12.!) ’Overdamped system’ URITEI12.40) ENDIF 10 FORMATI26I’!’)) 20 FORMATIZOI’!’)) 30 FORMATI27I’!’)) 40 FORMATIISI’!’)) C C PRINT!.’The following information can be printed’ PRINT!.’ to your output file.’ PRINT!. 1)Mass Matrix’ 9 PRINT!.’ 2)Stiffnees Matrix’ PRINT!.’ 3)Damping Matrix’ ’ 4)Critical Damping Matrix’ 0 PRINT!. 5)Undamped Natural Frequencies’ PRINT!.’ 6)Modal Matrix’ IFIDEVAL)PRINT!.’ 7)Damped Eigenvalues’ PRINT!.’ ’ PRINT!.’Enter the option numbers desired in a single line’ PRINT!.’ separated by commas. ’ PRINT!.’ ’ 50 READII.’(18A1)’)(RESPIII).II=1.16) II=1 URILEICONTlDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Urite Mass matrix C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 110 120 100 IFIRESPIII).EQ.’1’)TREN VRITEIIZ.’I )’) URITEI12.!) ’MASS MATRIX’ NltIN/Sll5 IFINI.EQ.N)NBLX3N/5 IFINI.NE.N)NBLKIN/5+1 DO 100 IB'1.NBLX IFIN.LE.I5IIB))TREN NCOL=N NUM'S-I5!IB-N) ELSEIFIN.GT.I5!IB))THEN NCOLISIIB NUM-S ENDIF DO 110 J35!IB-4.NCOL JJ'J-SIIIB-I) COLIJJIIJ CONTINUE IFIIB.EQ.1)TREN VRITEITEXT.SOO)NUM URITEI12.TEXT) URITEITEXT.910)NUM URITEI12.TEXT)ICOLIJJ).JJ.1.NUM) ELSE URITEITEXT.9201NUM URITEI12.TEXT) URITEITEXT.9301NUM URITEI12.TEXT)ICOLIJJ).JJ=1.NUM) ENDIF DO 120 I-l.N IFIIB.EQ.1)TNEN URITEI12.940)IMII.J).J=1.NCOL) ELSE URITEIl2.950)IMII.J).J85!IB-4.NCOL) ENDIF CONTINUE URITEIIZ.’( )‘l CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Brits Stiffness matrix C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSEIFIRESPIII).EQ.’2’)TNEN URITEI12.’( )’) URITEI12.!) 'STIFFNESS MATRIX’ NI‘IN/S)!5 IFINI.EQ.N)NRLK'N/5 IFINI.NE.N)NBLXIN/5+1 DO 200 IB'I.NBLK IFIN.LE.I5!IB))THEN NCOL=N NUMBS-I5!IB-N) ELSEIFIN.GT.I5!IB))TREN NCOL-5IIB 61 NUM=5 ENDIF DO 210 J85!IB-4.NCOL JJ=J-5!(IB-I) COLIJJ)=J 21o CONTINUE IFIIB.EQ.I)TREN URITEITEXT.900)NUM VRITEII2.TEXT) VRITEITEXT.910)NUM URITEI12.TEXT)(COLIJJ).JJ=I.NUM) ELSE VRITE(TEXT.920)NUM VRITEI12.TEXT) ‘ URITEITEXT.930)NUM E URITEI12.TEXT)ICOLIJJ).JJ=I.NUM) H ENDIF ' no 220 I=1.N IFIIB.EQ.1)TREN URITEII2.940)IKII.J).J=l.NCOL) ELSE URITEI12.950)IKII.J).J-5!IB-4.NCOL) ENDIF 220 CONTINUE URITEII2.’I 1') I zoo CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Urite Damping matrix C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSEIFIRESPIII).EO.’3’)THEN URITEI12.’( 1'1 VRITEII2.!) ’DAMPING MATRIX’ NI-(N/S)!5 IFINI.ER.N)NBLK=N/5 IFINI.NE.N)NBLK=N/5+1 DO 300 IB-1.NBLK IFIN.LE.(5!IB))TREN NCOLaN NUM=5-I5!IB-N) ELSEIFIN.GT.I5!IB))THEN NCOL-5!IB NUM=5 ENDIF DO 310 J35!IB-4.NCOL JJ-J-5!(lB-1) COLIJJ)-J 310 CONTINUE IFIIB.EQ.1)TNEN URITEITEXT.900)NUM VRITEII2.TEXT) URITEITEXT.910)NUM VRITEI12.TEXT)(COLIJJ).JJ=I.NUM) ELSE VRITEITEXT.920)NUM VRITEII2.TEXT) VRITEITEXT.930)NUM URITEII2.TEXT)(COLIJJ).JJ=1.NUM) ENDIF DO 320 I-1.N 62 IFIIB.EQ.1)THEN VRITEI12.940)(C(I.J).J=1.NCOL) ELSE URITEII2.950)(C(I.J).J=5!IB-4.NCOL) ENDIF 320 CONTINUE VRITEI12.’I )’) 300 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Urite Critical damping matrix C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSEIFIRESPIII).EQ.’4’)THEN URITEI12.’( )’) URITEI12.!) ’CRITICAL DAMPING MATRIX’ NI-(N/S)!5 IFINI.EQ.N)NBLK-N/5 IFINI.NE.N)NBLK=N/5+1 DO 400 IB-i.NBLK IFIN.LE.I5!IB))TREN NCOL-N NUM=5-I5!IB-N) ELSEIFIN.GT.I5!IB))THEN NCOL=5!IB NUM=5 ENDIF DO 410 J-5!IB-4.NCOL JJ-J-5!(IB-i) COLIJJ)-J 410 CONTINUE IFIIB.EQ.1)TREN URITEITEXT.900)NUM URITEII2.TEXT) URITEITEXT.910)NUM URITEI12.TEXT)(COLIJJ).JJ=1.NUM) ELSE URITEITEXT.920)NUM VRITEI12.TEXT) URITEITEXT.930)NUM VRITEI12.TEXT)ICOLIJJ).JJ=1.NUM) ENDIF DO 420 I-1.N IFIIB.EQ.1)THEN URITEI12.940)(CCRII.J).J=1.NCOL) ELSE VRITEI12.950)(CCRII.J).J=5!IB-4.NCOL) ENDIF 420 CONTINUE URITEI12.’( )’) 400 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Urite undamped natural frequencies C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSEIFIRESPIII).EQ.’5’)TREN VRITEII2.’( )’) VRITEI12.!) ’UNDAMPED NATURAL FREQUENCIES (rad/sec)’ VRITEI12.!) ’ ’ DO 500 I-l.N VRITEI12.510)I.UII) 510 FORMATII2.1X.EI3.6) l‘n_fl.1n'flj ‘3 63 500 CONTINUE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Vrite Modal matrix C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSEIF(RESP(II).EQ.’6’)THEN URITEIIZ.’( )’) URITEI12.!) ’EIGENVECTORS (MODAL MATRIX)’ NN=IN/8)!8 IFINN.EQ.N)NBLK=N/8 IFINN.NE.N)NBLK=N/8+1 DO 600 IJ=1.NBLK IFIN.LE.I8!IJ))TMEN NCOL=N NUM‘O'IOIIJ-N) ELSEIFIN.GT.(8!IJ))TREN NCOL=8!IJ NUM'B ENDIF URITEITEXT.650)NUM URITEI12.TEXT) DO 610 J88!IJ-7.NCOL JJ8J-8!(IJ-1) COLIJJ)=J 610 CONTINUE URITEITEXT.660)NUM URITEI12.TEXT)(COLIJJ).JJ=1.NUM) DO 620 I31.N URITEI12.67O)IUII.J).JI8IIJ-7.NCOL) 620 CONTINUE 600 CONTINUE 650 FORMATI’(/.’.I1.’(3X.6HColumn))’) 660 FORMATI’(4X.’.Il.’(IN(.I2.1R).5X))’) 670 FORMATI8F9.4) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Urite Damped Eigenvalues C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSEIFIRESPIII).EQ.’7’)THEN IFIDEVAL)TREN URITEI12.’I )’) URITEI12.!) ’DAMPED EIGENVALUES (Unpaired) ’ DO 700 I=I.N2 URITEI12.710)I.LAMDA(I) 710 FORMATII3.lX.El3.6.’+’.E13.6.’i’) 700 CONTINUE ELSE PRINT!.’!!You did not chose to compute damped eigenvalues!!’ ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Read delineators C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSEIFIRESPIII).EQ.’ ’.OR.RESP(II).EQ.’.’)TNEN II-II-l CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Unrecognized character C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ELSE PRINT!.‘Unrecognized character’ CONT=.FALSE. p.710; 'I'u. .- em? '- 64 ENDIF C C IIBII+2 IFIII.GE.16)CONT=.FALSE. ENDUMILE C C URITEIIZ.’I )’) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 900 FURNATI’(lr’rIlo’(QXDBHCOIUDNOAX))’) 910 FORMATI’I’.I1.’(5X.1H(.I2.1H).5X))’) 920 FORHATI’(l03XO’rIIO’IQXIBHCOIUIHOQX))’) 930 FORMATI’I3X.’.I1.’(5X.1H(.I2.1R).5X))’) 940 FORMAT(5(E13.6.1X)) 950 FORMATI2X.5(1X.E13.5)) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 65 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine RERUN C C C C RERUN asks the user if he wishes to continue and look at C C a new problem using a different damping matrix. C C CONT is true if the answer is yes. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE RERUNICONT) C CHARACTER!1 ANSIB) LOGICAL CONT PRINT!.’ ’ PRINT!.’ ’ PRINT!.’Do you wish to run again with a different C matrix?’ READII.’(3A1)’)ANS CALL UPCASEIANS) UNILEIANSII).NE.’Y’.AND.ANSI1).NE.’N’)DO PRINT!.’Unrecognized character. Answer yes or no.’ READII.’(3AI)’)ANS CALL UPCASEIANS) ENDVRILE IFIANSII).EQ.’Y’)CONT-.TRUE. IFIANSII).EQ.’N’)CDNT-.FALSE. RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine CHKSYM C C C C CNKSYM determines if the matrix MAT is symmetric. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE CRKSYMIMAT.N.SYM) C C REALNO MATIN.N) LOGICAL SYM SYM-.TRUE. [=1 URILEISYM.AND.I.LE.N)DO J=I+1 UNILEISYM.AND.J.LE.N)DO IFIDABSIMATII.J)-MATIJ.I)).GT.1.E-30)SYM=.FALSE. J=J+1 ENDUHILE I=I+1 ENDUHILE RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 66 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine PDCHK C C C C PDCRK determines if a matrix MAT is positive definite. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE PDCNKIMAT.N.TM.NN.PD) C C REAL!8 MATIN.N).TMINN).DPR.DET LOGICAL PD C C PD=.TRUE. C C MS=1 URILEIMS.LE.N.AND.PD)DO CALL DETERIDET.MAT.TM.MS.N.NN.1.DO) CALL DIAGPRIMAT.N.MS.DPR) IFIDET/DPR.LT.1.E-12)PD'.FALSE. MS8M5+1 ENDUHILE C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine PSDCRK C C C C PSDCHK determines if a matrix MAT is positive semi-definite. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE PSDCRKIMAT.N.TM.NN.PSD) C C REALIB MATIN.N).TMINN).DPR.DET LOGICAL PSD PSD8.TRUE. M581 UHILEIMS.LE.N.AND.PSD)DO CALL DETERIDET.MAT.TM.MS.N.NN.1.DO) CALL DIAGPRIMAT.N.MS.DPR) IFIDET/DPR.LT.-1.E-12)PSD8.FALSE. MS=MS+1 ENDURILE C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 67 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine DETER C C C C DETER computes the determinant of a matrix MAT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE DETERIDET.MAT.TM.MS.N.NN.C) REALIO MATIN.N).TMINN).DET.C' KK=1 DO 10 J81.MS DO 20 I81.MS ‘ TMIKK)=C!MAT(I.J) KK'KK+1 2O CONTINUE 10 CONTINUE CALL DMDETIDET.TM.MS.I1.I2.I3.I4.IERR) C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine DIAGPR C C c C DIAGPR computes the product of the nonzero diagonal of C C the matrix MAT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE DIAGPRIMAT.N.MS.DPR) REALIB MATIN.N).DPR DPR81.DO DO 10 I-1.MS IFIDABSIMATII.I)).GT.1.E-12)TREN DPR=DPRIMAT(I.I) ENDIF 10 CONTINUE DPR=DABSIDPR) C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 68 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine TRMTYP C C C C TRMTYP determines if the terminal being used is one that C C requires paging. NP=I if it Needs Paging. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE TRMTYPINP) . C C INTEGER ITI5) DATA IT/4014.4052.4006.4010.4114/ CALL TERMCOIITERM.ITYPE.IBAUD.IER) NP‘O DO 10 381.5 IFIITYPE.EQ.IT(J))NP-1 10 CONTINUE RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine PAGE C C C C PAGE pages the screen and resets the line counter NL C C to zero. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE PAGEINL) C C CALL NEUPAG CALL SLEEPOIZOOO) NL=O C C RETURN END C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MIIIIIIIIIHIIIIIIII“WWII II IIIIIBIIIIIIES 3 1293 03015 0352