UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analytical and experimental studies of the behaviour of equipment vibration isolators under seismic conditions Lam, Frank C. F. 1985

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1985_A7 L35.pdf [ 9.14MB ]
Metadata
JSON: 831-1.0062953.json
JSON-LD: 831-1.0062953-ld.json
RDF/XML (Pretty): 831-1.0062953-rdf.xml
RDF/JSON: 831-1.0062953-rdf.json
Turtle: 831-1.0062953-turtle.txt
N-Triples: 831-1.0062953-rdf-ntriples.txt
Original Record: 831-1.0062953-source.json
Full Text
831-1.0062953-fulltext.txt
Citation
831-1.0062953.ris

Full Text

A N A L Y T I C A L AND EXPERIMENTAL STUDIES OF THE BEHAVIOUR OF EQUIPMENT VIBRATION ISOLATORS UNDER SEISMIC CONDITIONS by F R A N K C. F. L A M B.A.Sc, The University of British Columbia, 1982 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF MASTER OF APPLIED SCIENCE in THE F A C U L T Y O F G R A D U A T E STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA April H85 e Frank C. F. Lam,/985 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, 1 agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It is understood that the copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Civil Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: April 85 ABSTRACT Analytical and experimental studies of the behaviour of equipment vibration isolators under seismic conditions are presented. A preliminary parametric study of the effect of equipment-structure interaction on the ultimate equipment response of general equipment-structure systems is considered first The results of this study indicate the conditions under which a non-interactive approach can yield adequate ultimate equipment response estimates. A model of a prototype air handling unit mounted on vibration isolators was constructed for use in the experimental studies. Two types of vibration isolators -elastomeric isolators and open spring isolators with uni-directional restraint - were tested under static and dynamic conditions. The frequency and seismic response characteristics of these vibration isolated systems were obtained. The experimental results indicate that the vibration isolators have nonlinear stiffness characteristics and high damping values. The results also show that the elastomeric isolators can survive a substantially higher level of base excitation than the open spring isolators with uni- directional restraint Analytical models of the vibration isolated systems, based on the model identification test results, have been formulated. A numerical procedure, ultilizing time series analysis, was used to solve the equations of motion of the systems. Good agreement between the experimental results and the analytical results was observed. This study indicates that the analytical procedure can be used to accurately predict the response characteristics of vibration isolated equipment systems subjected to known base excitation inputs. ii Table of Contents ABSTRACT ii LIST OF FIGURES vi LIST OF TABLES xiii ACKNOWLEDGMENTS xiv 1. INTRODUCTION 1 1.1 GENERAL : 1 1.2 LITERATURE REVIEW 3 1.2.1 The Seismic Response of Linear Equipment Systems 3 1.2.2 The Seismic Response of Nonlinear Equipment Systems 6 1.2.3 Dynamic Testing Techniques and Seismic Qualification Practice 7 1.3 OBJECTIVES AND SCOPE 7 2. ANALYTICAL STUDY ON THE EFFECT OF EQUIPMENT-STRUCTURE INTERACTION 9 2.1 INTRODUCTION 9 2.2 STEADY STATE RESPONSE '..9 2.2.1 Coupled Differential Equations of Motion 10 2.2.2 Uncoupled Differential Equations of Mouon 12 2.2.3 Comparsion of the Coupled and Uncoupled Steady State Response ....14 2.3 SEISMIC RESPONSE OF EQUIPMENT-STRUCTURE SYSTEM 15 2.3.1 System Parameters Studied 16 2.3.2 Solution Method 21 2.3.3 Response of Linear Equipment 22 2.3.4 Method of Shake Table Testing 23 3. TESTING FACILITIES AND TEST MODEL PARAMETERS 42 3.1 INTRODUCTION 42 3.2 TESTING FACILITIES 42 3.3 TEST MODEL CHARACTERISTICS 44 iii 4. EXPERIMENTAL STUDY 53 4.1 INTRODUCTION 53 4.2 STATIC TESTING PROGRAM 53 4.3 DYNAMIC TESTING PROGRAM 55 4.3.1 Sinusoidal Testing 55 4.3.2 Sinusoidal Decay Testing 56 4.3.3 Random White Noise Testing 57 4.3.4 Seismic Motion Testing 58 5. ANALYTICAL STUDY 62 5.1 INTRODUCTION 62 5.2 M O D E L IDEALIZATIONS 62 5.3 DEVELOPMENT O F T H E MATHEMATICAL M O D E L 65 5.3.1. The Structural Model 65 5.3.2 The Equipment Model 67 5.3.3 The Combined Equipment Structure Model 77 5.4 M E T H O D OF SOLUTION 80 6. EXPERIMENTAL A N D ANALYTICAL RESULTS 85 6.1 INTRODUCTION 85 6.2 M O D E L INDENTIFICATION TEST RESULTS (STATIC) 85 6.3 M O D E L INDENTIFICATION TEST RESULTS (DYNAMIC) 89 6.4 R A N D O M MOTION TEST RESULTS 91 7. CONCLUSIONS 125 7.1 SUMMARY A N D CONCLUSIONS 125 7.2 F U T U R E STUDIES 127 BIBLIOGRAPHY 128 APPENDIX A 131 APPENDIX B 134 iv A P P E N D I X C ! 139 v LIST OF FIGURES Figure Page 2.1 Two degree of freedom equipment-structure system 26 2.2 Uncoupled and coupled equipment response transfer function. e=0.01, $. = 0.01 27 2.3 Uncoupled and coupled equipment response transfer function. e=0.001, $ . = 0.01 28 i 2.4 Uncoupled and coupled equipment response transfer function. e =0.0001, $ . = 0.01 '. 29 2.5 Uncoupled and coupled equipment response transfer function. 6=0.01, £ . = 0.02 30 2.6 Uncoupled and coupled equipment response transfer function. 6=0.001, $. = 0.02 31 2.7 Uncoupled and coupled equipment response transfer function. e =0.0001, $. = 0.02 32 2.8 Uncoupled and coupled equipment response transfer function. 6=0.01, $. = 0.05 33 2.9 Uncoupled and coupled equipment response transfer function. e=0.001, $j = 0.05 34 2.10 Uncoupled and coupled equipment response transfer function. e=0.0001, $ . = 0.05 35 2.11 Uncoupled and coupled equipment response transfer function. e=0.01, $ . = 0.10 36 ° I vi Figure Page 2.12 Uncoupled and coupled equipment response transfer function. e =0.001, $. = 0.10 37 o 1 2.13 Uncoupled and coupled equipment response transfer function. e =0.0001. 5 . = 0.10 38 2.14 Three degree of freedom structural system 39 2.15 Single degree of freedom equipment system 39 2.16 Combined equipment-structure system 39 2.17 Coupled and uncoupled peak response of tuned secondary system to El Centra earthquake N-S component. 2% damping in structure, 5% damping in equipment. 40 2.18 Coupled and uncoupled peak response of tuned secondary system to El Centro earthquake N-S component 2% damping in structure, 2% damping in equipment 41 3.1 Overall function of U.B.C. shaking table 48 3.2 Experimental model showing instrument locations 49 3.3 Prototype air handling system 50 3.4 Experimental model 51 3.5 Open spring isolator with uni- directional restraint 52 3.6 Elastomeric isolator 52 4.1 Schematic respresentation of bi-directional static test 60 4.2 Static test set up . 60 4.3 Close up view of bi-directional static test 61 4.4 Typical measured load-deflection curve 61 5.1 Idealized model of the vibration isolated equipment system 82 vii Figure Page 5.2 Nonlinear load-deflection curve (static test) and piecewise linear load-deflection curve (analytical model) 83 5.3 Schematic representation of Coulomb friction between isolator and snubber device 83 5.4 Equilibrium forces in a displaced vibration isolated equipment system 84 6.1 Axial load-deflection curves of an elastomeric vibration isolator. ..94 6.2 Lateral load-deflection curves of an elastomeric vibration isolator. 95 6.3 Axial load-deflection curves of an open spring isolator with uni-direcuonal restraint. 96 6.4 Lateral load-deflection curves of.an open spring isolator with uni-directional restraint. 97 st 6.5 Decaying sinusoidal response time history at 1 natural frequency (1.8 Hz): open spring isolator with uni-direcuonal restraint 98 6.6 Decaying sinusoidal response u'me history at 3r<^ natural frequency (4.9 Hz): open spring isolator with uni-directional restraint .' 99 6.7 Horizontal acceleration frequency response function: open spring isolator with uni-directional restraint 100 6.8 Rotational acceleration frequency response function: open spring isolator with uni-directional restraint 101 6.9 Predicted and measured mode shapes: open spring isolator with uni-directional restraint 102 viii Figure Page st 6.10 Decaying sinusoidal response time history at 1 natural frequency (4.6 Hz): elastomeric isolator 103 6.11 Decaying sinusoidal response time history at 3TC^ natural frequency (15 Hz): elastomeric isolator 104 6.12 Horizontal acceleration frequency response function: elastomeric isolator 105 6.13 Rotational acceleration frequency response function: elastomeric isolator 106 6.14 Predicted and measured mode shapes: elastomeric isolator 107 6.15 Predicted and measured horizontal acceleration response time history: open spring isolator with uni- directional restraint Base motion: El Centro 1940 N-S. Peak acceleration scaled to O.lg ;..108 6.16 Predicted and measured rotational acceleration response time history: open spring isolator with uni-directional restraint Base motion: El Centro 1940 N-S. Peak acceleration scaled to O.lg 109 6.17 Predicted and measured horizontal acceleration response time history: open spring isolator with uni-directional restraint Base motion: San Fernando 1971 N21E. Peak acceleration scaled to 0.18g 110 6.18 Predicted and measured rotational acceleration response time history: open spring isolator with uni-directional restraint Base motion: San Fernando 1971 N21E. Peak acceleration scaled to 0.18g I l l ix Figure Page 6.19 Predicted and measured horizontal acceleration response time history: open spring isolator with uni-directional restraint. Base motion: Taft 1950. Peak acceleration scaled to O.llg 112 6.20 Predicted and measured rotational acceleration response time history: open spring isolator with uni-directional restraint. Base motion: Taft 1950. Peak acceleration scaled to O.llg 113 6.21 Predicted and measured horizontal acceleration response time history: open spring isolator with uni-directional restraint Base motion: Band limited (0 - 25 Hz) white noise signal. Peak acceleration scaled to O.lg 114 6.22 Predicted and measured rotational acceleration response time history: open spring isolator with uni-directional restraint Base motion: Band limited (0 - 25 Hz) white noise signal. Peak acceleration scaled to O.lg 115 6.23 Predicted and measured horizontal acceleration response time history: elastomeric isolator. Base motion: El Centro 1940 N-S. Peak acceleration scaled to O.lg 116 6.24 Predicted and measured rotational acceleration response time history: elastomeric isolator. Base motion: El Centro 1940 N-S. Peak acceleration scaled to O.lg 117 x Figure Page 6.25 Predicted and measured horizontal acceleration response time history: elastomeric isolator. Base motion: San Fernando 1971 N21E. Peak acceleration scaled to 0.18g 118 6.26 Predicted and measured rotational acceleration response time history: elastomeric isolator. Base motion: San Fernando 1971 N21E. Peak acceleration scaled to 0.18g 119 6.27 Predicted and measured horizontal acceleration response time history: elastomeric isolator. Base motion: Taft 1950. Peak acceleration scaled to 0.1 lg 120 6.28 Predicted and measured rotational acceleration response time history: elastomeric isolator. Base motion: Taft 1950. Peak acceleration scaled to 0.1 lg 121 6.29 Predicted and measured horizontal acceleration response time history: elastomeric isolator. Base motion: Band limited (0 - 25 Hz) white noise signal. Peak acceleration scaled to O.lg 122 6.30 Predicted and measured rotational acceleration response time history: elastomeric isolator. Base motion: Band limited (0 - 25 Hz) white noise signal. Peak acceleration scaled to O.lg 123 xi Figure Page 6.31 Failed open spring isolator with uni-directional restraint. 124 6.32 Failed open spring isolator with uni-directional restraint. 124 A.1 Power spectral density Function of a band limited white noise signal 133 A. 2 Power spectral density function 133 B. l Single input / mutiple output system 135 xii LIST OF TABLES Table Page 6.1 Model and prototype inertia properties. 86 xiii A C K N O W L E D G M E N T S 1 wish to express my sincere gratitude and appreciation to my three advisors, Dr. S. Cherry, Dr. N. D. Nathan, and Dr. D. L. Anderson, for their encouragement and guidance throughout the research and preparation of this thesis. The cooperation received from Mr. R. A. Strachan of Brown Strachan Associates was very helpful in this investigation and it is sincerely appreciated. My thanks are also extended to fellow graduate students, Mr. Bill Lipsett, Mr. Ken Tarn, and Mr. Bill Barwig, for their valuable advice during the course of this thesis. This thesis was made possible by the financial assistance of the National Research Council of Canada in the form of an research assistantship. The donation of the vibration isolators by Airmaster Sales Ltd. and Air System supplies Ltd. is also appreciated. 1 also wish to thank Genstar Structures Ltd. for partially donating the experimental model. Finally 1 wish to express my thanks to the faculty and technical staff of the Civil Engineering Department at U.B.C, fellow graduate students, and my parents for their encouragement, advice, and assistance. xiv Chapter 1 INTRODUCTION 1.1 GENERAL The protection of equipment from seismic hazards is of wide engineering interest. Essential equipment such as pressure vessels, pumps, power generators, fire sprinkling systems, control equipment, and piping are vital to certain lifeline and emergency systems. The operability of these systems during and after an earthquake is of extreme importance. Within these lifeline and emergency systems there are some items of equipment, such as air handling units and chillers, which are not vital to the operation of these systems. Compared with essential equipment, it is not as important that these non-essential items remain operative during or after an earthquake. However, it is dangerous to suggest that their seismic restraint criteria can be relaxed; during an earthquake such units could be displaced from their foundations, causing damage to neighboring essential equipment Past earthquakes have shown that the failure of equipment systems which lead to the breakdown of emergency services can produce loss of human life and extensive property damage. Also, the failure of lifeline systems such as transportation and energy transmission facilities, and the disruption of communication and water supply networks can cause extreme disorder in a city, which may in turn increase the potential for various secondary disasters. In addition, the economic consequences of equipment failure can be significant Specialized equipment systems may be much more expensive than the value of their housing structures because of the considerable investment in their engineering development, manufacturing, and construction costs. Furthermore, the time required to replace damaged specialized equipment may involve a revenue loss which further compounds the initial financial loss caused by an earthquake. In 1973, the California Division of Mines and Geology estimated that, in California alone, 1 2 twenty-one billion dollar worth of damage could be expected between the years 1970 and 2000 as a result of equipment failure due to earthquakes if the then current equipment design principles were not improved [1]. Thus the main motivations for developing criteria for the seismic protection of equipment are improvement in safety and protection of capital investment. Although the consequences of equipment failure can be severe, most equipment systems, with the exception of nuclear power plants, power transmission stations, and defence installations, are rarely designed with the degree of care actually warranted. The difficulties in providing for the proper seismic restraint of equipment systems arise from three sources: 1. Most equipment systems and their housing structures are analysed and designed as separate units. This practice neglects the effect of interaction between the primary structural system and the secondary equipment system; therefore, the design may be inadequate. 2. The seismic response of an equipment system depends on the characteristics of the equipment, of the housing structure, and of the ground motion, and also on the location of the equipment within a building. Therefore, the designer faces a very difficult task when attempting to represent the seismic environment accurately. 3. Most equipment systems have not been fully qualified, either analytically or experimentally, to withstand the seismic environment The risk of catastrophic failure of these systems during an earthquake may therefore be high. Vibration isolated equipment is particularly vulnerable to seismic forces. Typically, this type of equipment fs mounted on special spring units, or other resilient supports, in order to significantly reduce any tranmission of shock or periodic forces into the primary system as a result of normal equipment operation. Motion restraining devices are sometimes used to prevent the development of excessive response in the 3 resiliently mounted equipment Experience from past earthquakes has demonstrated that resilienfjy mounted equipment without motion restraining devices may be totally destroyed even though their housing structures may suffer relatively minor damage. In such cases, failure may result due to the dislodgement of the isolator spring units, which can damage anything situated near the equipment Moreover, some resilient mountings with motion restraining devices have not been extensively tested either in actual earthquakes or in a simulated seismic environment A reliable qualification program, including a consistent analysis and design procedure, is needed if the resiliently mounted equipment is to be protected from earthquake damage. 1.2 LITERATURE REVIEW During the past decade a great deal of research has led to a considerable advance in our understanding of the response of equipment systems under earthquake excitation. These studies can be broadly categorized into three major groups: 1. The analytical study of seismic response of linear equipment systems. 2. The analytical study of seismic response of nonlinear equipment systems. 3. Techniques of dynamic testing and seismic qualification. 1.2.1 T H E SEISMIC RESPONSE OF LINEAR EQUIPMENT SYSTEMS The earliest methods of seismic analysis of equipment systems adopt a non-interactive approach, in which the primary structural system is analysed separately from the secondary equipment system. The time history response of the primary system at the point of attachment of the secondary system is first obtained from a time step analysis. A response spectrum, generated from the computed time history record, is then used as the input excitation to the secondary system. Since the generation of a time history record is usually a lengthy and costly procedure, various authors have developed simpler approaches. 4 Bigg and Roessei [2], Amin et al [3], and Kapur and Shao [4] proposed empirical rules for generating floor response spectra directly from ground response spectra and the modal properties of the primary structure. Singh [5], Chakravarti and Vanmarke [6], Singh [7], and Vanmarke [8] developed alternate approaches based on random vibration theory; Scanlan and Sachs [9] suggested a technique involving the Fourier transform method. All of these conventional floor response spectrum methods have several important drawbacks. They consistently experience problems when the mass of the secondary system is not negligible in comparision with the mass of the primary system and when a natural frequency of the secondary system coincides with a natural frequency of the primary system. Crandall and Mark [10], Singh [5], Kapur and Shao [4] point out that under these conditions the effect of interaction between the primary system and the secondary system is important and significant errors are introduced if an uncoupled analysis is used. Der Kiureghian and Asfura • [11] suggest that proper modal combination rules cannot be developed if the response is artificially separated into a dynamic and a pseudo static part They also state that the effect of non-classical damping in the combined system, and the cross correlation between closely coupled modes in the primary system and the secondary system, is often neglected or improperly treated using a conventional floor response spectrum method of analysis. Various researchers have investigated different analytical methods designed to take it into account the important effect of dynamic interaction between the primary and secondary systems. In these methods, the combined equations of motion of both the primary and the secondary systems must be considered. Crandall and Mark [10] have found the root mean square response to stationary excitation of a two degree of freedom system using the combined equations of motion. Penzien and Chopra [12] proposed an innovative method which reduces 5 an N degree of freedom primary system and a single degree of freedom secondary system into N, two degree of freedom systems. However, the input excitation for the above method requires a two degree of freedom ground response spectrum, which is not usually available. In another development, Newmark [13] first introduced the concept of effective mass ratio to arrive at approximate solutions to simple equipment-structure interaction problems. Nakahata et al [14] later improved upon this concept to approximate modal properties and modal responses with various rules for modal combination. Sackman and Kelly [15] recognized that the perturbation method can be used to formulate equations of motion for a general primary system and a light single degree of freedom secondary system. They arrived at closed form solutions of the perturbed equations of motion using Laplace transform techniques. Perturbation methods are mathematical tools which can be used to solve differential equations containing small parameters. In this study the small parameters are the mass, stiffness, and damping values of the secondary system. Robinson and Ruzicka [16], Villaverde and Newmark [17], and Der Kiureghian et al [18] also used perturbation methods and modal combination rules to arrive at approximate modal properties of the combined system. Robinson and Ruzicka [16] used a Fourier transform approach to solve the perturbed differential equations. Villaverde and Newmark [17] considered secondary systems with multiple attachment points and non-classical damping. Der Kiureghian, Sackman, and Nour-Omid [18] studied a single degree of freedom secondary system and derived a floor response spectrum which can include the effect of equipment-structure interaction by using random vibration techniques. Der Kiureghian and Igusa [19] extended this method to multiply tuned and arbitrarily supported secondary systems with non-classical damping. Finally, Der Kiureghian and Asfura [11] develped a new floor response spectrum method which includes 6 the effect of equipment-structure interaction for multiply supported secondary systems. Good analytical solutions to the problem of linear equipment-structure interaction were obtained by these researchers. 1.2.2 T H E SEISMIC RESPONSE OF NONLINEAR EQUIPMENT SYSTEMS The analysis of the seismic response of a nonlinear equipment system is much more difficult than that of a linear equipment system. A popular analytical approach to the nonlinear problem makes use of numerical time step analysis which accounts for a proper modelling of the equipment nonlinearities. A good example of a nonlinear equipment system is illustrated by the case of resiliency supported equipment with motion limiting constraints. The nonlinear characteristics arise when the equipment system comes into contact with the stiff motion limiting device. Iwan [20, 21] studied this problem analytically by replacing the set of nonlinear springs by a set of equilvalent linear springs whose stiffness depends on the displacement amplitude. From a conventional floor response spectrum, the maximum kinetic energy of the equilvalent linear system can then be obtained as a function of natural frequency. Rayleigh's principle states that, in any one mode of a conservative system, the maximum potential energy, which depends on the displacement amplitude, must equal the maximum kinetic energy, which depends on the natural frequency. Using this principle, Iwan arrived at a function which relates the displacement amplitude to the natural frequency of the equipment. He obtained another function relating the displacement amplitude to the natural frequency of the equilvalent linear support stiffness, which depends on displacement amplitude. The transcendental equation obtained by combining these relations is solved by an iterative procedure which yields an estimate of the seismic response of the equipment system. This method offers an equilvalent linear, non-interacting seismic analysis for isolated equipment 7 systems with motion restraint devices. 1.2.3 DYNAMIC TESTING TECHNIQUES AND SEISMIC QUALIFICATION PRACTICE Rigorous guidelines for seismic qualification tests are given by the IEEE standards for seismic qualification of Class IE equipment for nuclear power generating stations [22], and by Silva [23], McGavin [24], and Wilson et al [25]. Rainer [26] and Nielson [27] have suggested procedures for the dynamic testing of structures. Bendat and Piersol [28, 29] and Newland [30] improved the techniques of spectral analysis for dealing with random data. Sabnis et al [31] provided insights into structural modelling and experimental techniques. Some of these experimental techniques are used in the testing program in the present study. 1.3 OBJECTIVES AND SCOPE The main objectives of this research are: 1. to investigate the seismic response of vibration isolated equipment systems 2. to study the effect of equipment-structure interaction 3. to evaluate the seismic performance of two typical equipment vibration isolators. The detailed investigation presented in this thesis can be divided into three parts. Chapter two involves a parametric study of the effect of equipment-structure interaction on the ultimate equipment response of a general equipment-structure system. It is found that the effect of equipment-stucture interaction can significantly influence the dynamic response of some equipment-structure systems. Next, consideration is given to the experimental investigation of the structural dynamic properties (natural frequencies, mode shapes, and damping values) and the seismic response of a model of an equipment system mounted on commercially available vibration isolator units. The 8 modelling techniques and experimental methods are presented in Chapter three and four. The final stage of the study pertains to the development of an analytical model of vibration isolated equipment units and the solution of the theoretical problem by means of a time series analysis method. The analytical study is discussed in Chapter five and the analytical and experimental results are presented in Chapter six. Chapter seven contains the summary and conclusions of the results of the investigation. Chapter 2 ANALYTICAL STUDY ON T H E EFFECT OI EQUIPMENT-STRUCTURE INTERACTION 2.1 INTRODUCTION As stated in the literature review, the effect of equipment-structure interaction, if improperly' treated, can introduce significant errors in the dynamic analysis of the equipment system. The effect of equipment-structure interaction is most significant when the mass of the secondary system is not negligible in comparision with the mass of the primary system or when a natural frequency of the secondary system is tuned to a natural frequency of the primary system. The purpose of this chapter is to study the effect of equipment-structure interaction by analysing the response of tuned equipment-structure systems and by reviewing the important parameters which can affect their response. The first section of this chapter deals with the steady state response of an item of equipment having a single degree of freedom which is attached and tuned to a single degree of freedom structure. The second section examines the seismic response of such an unit when attached to a three degree of freedom structure. Finally, the last section discusses shake table testing methods which can include the effect of equipment-structure interaction. 2.2 STEADY STATE RESPONSE The study of equipment-structure interaction is best introduced with the simplest and most fundamental two degree of freedom assemblage, such as shown in Fig. 2.1. This simple system contains all the essential properties which characterize the more general multi- degree of freedom equipment-structure system; therefore, the results from this study can easily be extended to more complicated equipment-structure systems. 9 10 2.2.1 COUPLED DIFFERENTIAL EQUATIONS OF MOTION The coupled equations of motion of the simple two degree of freedom system are formulated in terms of the parameters of the individual fixed based single degree of freedom subsystems. The physical properties of the two subsystems are mass influence coefficients m., damping influence coefficents c., and stiffness influence coefficents k ,^ where i = l refers to the primary system and i = 2 refers to the secondary system. Let M , C, and K, be the mass, damping, and stiffness matrix of the combined system respectively. Then: m, 0 Ci + C2 -c 2 kj + k2 - k : M = C = K = _0 m2 _ _-C2 C2 _-k2 k2 _ The equation of motion of the two degree of freedom system (2.1) horizontal base motion x (t) is g Mx + Cx + Kx=-MI x (t) g (2.2) T T T where x = [xj, x2] , x = [x,, x2] , X =[X J ( X2] , are the horizontal displacement, velocity, and acceleration vectors measured relative to the base motion and I =[1, 1] is the influence vector coupling the input base motion to the individual degrees of freedom of the system. Now consider harmonic base motion. The base motion is given in complex form as X (t)=X e g v ' g ILOt (2.3) 11 where X is the base acceleration amplitude and CJ is the base motion driving g frequency. The steady state response of the two subsystems can be expressed in complex form as xI(t) = BI e i ( w t " * ) = A l e i w t x,(t) = B, e i ( u t ^ ) = A ! eiu)t (2.4) (2.5) where A, and A 2 are complex; they contain information on the amplitudes as well as the phases of the steady state damped response of the subsystems due to a harmonic base excitation. Substituting equations (2.1), (2.3), and the derivatives of equations (2.4) and (2.5) into equation (2.2) yields the following matrix equation (k 1+k 2 -m 1a> 2)+ia>(c l + c2) - ( k 2 + iwc 2) - ( k 2 + ia>c2) (k 2 -m 2cj 2)+icL)c 2 A = -ni iX 6 v -m 2 x (2.6) T Equation (2.6) can be solved to find the response vector, A = [A l t A j , and subsequent substitution of A into equation (2.4) and equation (2.5) gives the subsystem responses. The absolute acceleration of the secondary system is given by x2 a(t)=x I(t) + X (t) (2.7) The complex transfer function of the secondary system from a coupled analysis, H 2 C ( w ) , represents the amplitude and the phase of the steady state damped response of the secondary system due to a unit harmonic base excitation. 12 Therefore, by setting x (t) = l e 1 £ J ^ , the amplitude of the complex equipment c transfer function is: 1 n | H 2 c ( « ) | = ((A2 + B 2 ) i / Z ) / C (2.8) where A = D 2 + E 2 - (DF+EG)u 2 B = ( E F - D G > J ; C = D 2 + E 2 D=k 3 2-(k, + k 2-m 1oj3Xk2-m ;cj 2) + ((ci + c2)c2-2c22>j2 E=-aj((c1 + c2)k2 + c2k1 + k 2- 2k2c2) + (m2(ci + c2) + m!C2Xj3 F=m2(k, + k2) + mik2-mimjO)2 G = (m2Ci + c2 + miC2)cj From equation (2.8), | H 2 C (<^ ) | can be obtained as a function of driving frequency for various systems, having damping ratios, $ j, and mass ratios, e, where $ . = C./(2CL> .m.) (i = l , 2) and e =m I/m 1 . 2.2.2 UNCOUPLED DIFFERENTIAL EQUATIONS OF MOTION For comparative purposes, now consider an uncoupled analysis of the same two degree of freedom system. Assuming that the response of the uncoupled subsystems is harmonic, and using the same approach illustrated in the previous section, the absolute acceleration response of the uncoupled single degree of freedom structure, X ^a(0. is: X l a (t)=(l + ((Q-Ri)m1a;J)/(Q2 + R 2 )X g e l u t (2.9) 13 where Q = k 1 -w J m 1 R=wc, Using x 2 a (0 a s ^ D a s e moiion applied to the uncoupled equipment, and repeating the above procedure, the absolute acceleration response of the uncoupled single degree of freedom equipment system, ^ a ^ ' ' s : x 2 a(t)=(l + S)X g e l c J t (2.10) where S = (T-Ui)mju l /V T = ( Q Y - R Z ) u ' m , + ( Q 2 + R 2 ) Y U = - ( ( Q 2 + R 2 ) Z + ( Q Z + YR)m,ca 2) Y = ( Q 2 + R J ) ( Y 2 + Z 2 ) Y = kj-w2m2 Z = C J C 2 The amplitude of the complex equipment transfer function , l ^ u c ^ ^ ' ^ r o m the uncoupled analysis then becomes: H 2 u c ( C J ) | =((V + Tm2c;2)2 + (Um 3 a; 2 ) 2 ) 1 / 2 )/V (2.11) Again,| H (^) | as a function of driving frequency for systems with different damping ratios, $ ., and mass ratios, e, is easily obtained. 14 2.2.3 COMPARSION OF T H E COUPLED AND UNCOUPLED STEADY STATE RESPONSE Fig. 2.2 - Fig. 2.13 show equipment transfer functions as a function of driving frequency obtained from the coupled and the uncoupled analysis for systems with different parameters. The results from a steady state uncoupled analysis, in some cases, are very conservative as compared with the results from a steady state coupled analysis; but in other cases, the uncoupled analysis underestimates the equipment response. The accuracy of. the uncoupled analysis significantly depends on' the mass ratio and the damping ratio of the equipment-structure system, and the driving frequency of the harmonic base motion. In general, when the driving frequency coincides with or is close to the tuned natural frequencies of the subsystems, the uncoupled analysis grossly overestimates the equipment response. However, there exists a range of driving frequencies on either side of the tuning frequency for which the uncoupled analysis yields nonconservative results. Moreover, as the level of damping of the subsystems increases, or as the mass of the secondary system becomes negligible as compared to the mass of the primary system, the accuracy of the uncoupled analysis increases. For a tuned two degree of freedom equipment-structure system, the combined coupled system has two natural frequencies on either side of the tuning frequency, i.e. there is a frequency shift in the combined system away from the tuning frequency of the uncoupled system. This frequency shift depends on the various system parameters and it is around these shifted tuning poles that a band of high amplification occurs. Therefore, if the forcing frequency happens to lie in this region, depending on the system parameters, the uncoupled analysis may be unsafe. The response of the tuned coupled combined system involves energy exchanges between the equipment and its housing structure at a 15 beat frequency which is much lower than the tuned frequency of the subsystems. It is this classic beating phenomenon which the uncoupled analysis fails to take into consideration. The coupling between the subsystems depends on the magnitude of the damping and mass ratios of the subsystems. If the mass of the secondary system is negligible compared to the primary system, the coupling between the subsystems is weak and the location of the tuning poles obtained from a coupled analysis does not differ very much from the location of the tuning pole obtained from an uncoupled analysis. Therefore, in this situation, using an uncoupled analysis may be justified. Futhermore, the energy transfer between the subsystems takes time and requires many cycles of oscillation. In heavily damped systems, some of the energy may be dissipated before it is transmitted between the subsystems. As a result, in systems with high damping the error of an uncoupled analysis is usually not significant. 2.3 SEISMIC RESPONSE OF EQUIPMENT-STRUCTURE SYSTEM The seismic response characteristics and the steady state response characteristics of the tuned equipment-structure system can be quite different Earthquake motions are not purely harmonic in nature, since they contain more than one predominant forcing frequency. Moreover, the duration of the ground motion is such that the disturbance can be considered as a transient disturbance. In the following sections, the seismic response of a tuned, linear single degree of freedom equipment system attached to the top floor of a three degree of freedom structure is obtained from a time step analysis. An uncoupled analysis of the equipment-structure system is considered first Then a coupled analysis of the equipment-structure system is" presented. Various system parameters corresponding to different tuning conditions and different damping and mass ratios are also studied. 16 2.3.1 SYSTEM PARAMETERS STUDIED The structural system under consideration is a simple three story shear frame building with masses lumped at each floor as shown in Fig. 2.14. The mass of each floor is 175,000 kg and the interstory stiffness is 350,000,000 N/m. The damping of the structure is assumed to be classical and the damping ratio is taken as 2% for each mode ($ g . = 2% i = l , 2, 3). The equation of motion of the structural system subjected to base motion X (t) is: M X +C x +K x = - M I X (t) (2.12) s s s s s s s g T T where x g = [ x s j . XS2> ' s ^ e relative displacement, x g = t x s ^ XS2> *sj\ T is the relative velocity, and X g = [ x s | . xs2> * s the relative acceleration of each floor measured with respect to the base. The influence vector 1 =[1, 1, T 1] couples the base input motion to the individual floor degrees of freedom. M , C , K are the mass, damping, and stiffness matrix of the structure. M s s s  v ° s and K are as follows: s 175 175 x l O 3 175 (2.13) 7.0 -3.5 K = s -3.5 7.0 -3.5 x l O 8 N/m -3.5 3.5 The structural damping ..matrix can be found from the structural mass matrix M 17 and the structural stiffness matrix K s - By considering the structure without damping, the structural mode shape matrix $ and the structural natural period vector T can be found from the eigenvalue solution to the following undamped free vibration problem: M X +K x =0 (2.14) s s s s where 0 is a zero vector. Assuming the free vibration response is harmonic, the equation of motion can be expressed as [K - u 2 M ] $ r=0 * (2.15) 1 s r s s where $ s r is the mode shape for the r * mode of the structure's vibration and cJ r is the corresponding natural frequency. A nontrival solution is possible only when the determinant of equation (2.15) vanishes, yielding the frequency equation of the system: | K s - w 2 M s | | = 0 (2.16) Expanding the determinant yields an degree algebraic equation in u | for a N degree of freedom system. The N roots of this equation are the square of the N natural frequencies of the system. Assuming an arbitrary value for the th r first element in the r mode shape, $ s can be found by substituting in turn the natural frequencies for each mode into equation (2.15). The resulting structural modal matrix $ and the structural natural periods T are as follows: s s 18 * = s L s s .0780 -.1762 -.1413 .1413 -.0784 0.1762 .1762 0.1413 -.0784 x l O - 2 m T S = [T,. T,, T J = [.316. .113, .078] sec (2.17) We assume the modal matrix $ s can uncouple the damping matrix, as occurs in the classical approach. Thus, pre- and post- multiplying the damping matrix by the modal matrix yields a diagonal matrix of generalized damping coefficients: C *= * T C * = 2 s s s s 5 2 u>2m2 $3 cj3m3 (2.18) Then the structural damping matrix C s is: - T * -1 s s s s (2.19) -1 T Instead of inverting <i> one can make use of the relationship $ = $ M s s s s to yield the symmetrical structural damping matrix as: C = M * C $ M = s s s s s s 424 -124 - 270 397 -151 SYM. 273 xKT NS/M (2.20) 19 Now consider the characteristic of the linear single degree of freedom equipment item shown in Fig. 2.15. Assuming the secondary system's natural period is tuned to one of the natural periods of the structure, and defining the ratio, e, as the mass ratio between the mass of the secondary system, M f , and the mass of a floor of the building, m^ = 175,000 kg, one can select the equipment stiffness characteristics for a given mass ratio and a given tuning condition. For an equipment system having a mass ratio, e and tuned to the til i mode of the structure, the equipment stiffness is K e = ((27r)/Tsj)2 enij. The damping of the equipment system is assumed to be classical and the modal damping ratio is varied from 2% to 5% ($ e = C £ / (2M e w ) = 2% or 5%). The equation of motion of the uncoupled equipment system subjected to a base motion (x, + x ) is: M e X e + 2$ e ( K e M £) 1 / 2 x e + K £ x e = - M e(Xg(t)+x,(t)) (2.21) where X3, which is obtained from an uncoupled analysis of the building alone, is the acceleration of the third floor relative to the ground which is undergoing an acceleration X . g Consider next the equations of the coupled system. The combined system of a single degree of freedom equipment with a given tuning condition and mass ratio attached to the top floor of the three story structure is shown in Fig. 2.16. The equation of motion of the coupled system is: MX + Cx + K x =-MI X (t) (2.22) T T where x =[x i , x 2, x 3, x j is the displacement, x =[% u x 2, x 3, x j is the T velocity, and X =[Xi, X2, X3, x j is the acceleration of each floor relative to 20 j the base. Here I =[1, 1, 1, 1] is the influence vector coupling the base input motion to the individual degrees of freedom. The mass matrix and the stiffness matrix of the combined system are as follows: M = 175 175 175 el75 xlO 3 kg (2.23) K = 7.0x10' -3.5x10* •3.5x10° 7.0x10 8 -3.5x10' -3.5x10 7.0x10 8 + K ( - K •K K N/m (2.24) where K £ depends on the choice of the mass ratio e and the tuning condition. When dealing with tuned subsystems with large differences in system parameters, such as large differences in mass or stiffness, the damping of the combined system needs to be treated with extra care. In such cases, as pointed out by Hurty and Rubinstein [32], if the damping ratios between the subsystems are different, the combined system can no longer be considered as classically damped even though the individual subsystems have classical damping. This is because the mode shapes of the combined system are complex and there is a phase relationship between the movement of each degree of freedom when the system undergoes free vibration in a given mode. Failure to include the effect of non-classical damping in the analysis can introduce significant errors. The combined system damping matrix C is built from the damping matrix of the individual subsystems C s and C^. 21 424x10 -124x103 C = 3 -270xl0 J These system parameters will be used in the solution to the equations of motion; the method of solution is examined in the next section. 2.3.2 SOLUTION M E T H O D The method of solution used here for both the uncoupled and the coupled analysis is the well known Newmark Beta method, which assumes a specified variation in the change of acceleration response over a small time interval. The differential equation of motion in incremental form for a given time history base acceleration record is given by Clough and Penzien [33] as: MAX(t) + CAx(t) + KAx(t) =-MIAX (t) (2.26) % The application of the numerical step by step integration technique to this differential equation of motion yields the following results for each time step: -124x10' 397x10' -151x10" -270x10 3 -151x10 J 273xl03 + C - C -c Ns/m (2.25) X(t+At) = X(t)+Ax(t)/(/3At2 )-X(t)/(0At)-X(t)/(2/3) x(t+At) = x(t)+X(t)At(l-a/(20)hx(t)a/0 +Ax(t)a/(0At) x(t+At) = x(t) + K * _ 1 F * (2.27) K * = M/( |3At 2 ) + aAx(t)/(/3At) + K F *= -M(x (t+At)-X (t)) + M(X(t)/2+x(t)/At)//3 +C(-X(t)(l-a /(20 ))At+ x(t)a /0) 22 The values of a and 0 depend on the assumption of the variation of the response acceleration within each time step. If a constant variation in acceleration, taken as the average value between the acceleration at the start and at the end of a time interval, within each time step is assumed, a = .5 and /3 = .25. If a linear variation in acceleration within each time step is assumed, a = .5 and 0 = .16667. The average acceleration method is unconditionally stable, i.e. the solution is stable for any chosen length of time step, and the accuracy will, in general, increase as the chosen time step interval is decreased. The linear acceleration method, however, is not unconditionally stable, i.e. the solution may oscillate with increasing amplitude if too large a time step is chosen. The stablility criterion for the linear acceleration method is that the chosen time step must be less than 0.551 of the smallest natural period of the system. Moreover, both methods generally yield good results if the chosen time step interval is less the one-tenth of the smallest natural period of the system. 2.3.3 RESPONSE OF LINEAR EQUIPMENT A comparison of the results from coupled time step analyses and uncoupled time step analyses for the linear single degree of freedom equipment, with various mass ratios and tuning conditions, attached to the three story building is shown in Fig. 2.17 - Fig. 2.18. These figures indicate that the major difference between the seismic response of equipment systems and steady state response of equipment systems is that the uncoupled seismic analyses generally overestimate the peak equipment response, whereas the uncoupled steady state analyses may sometimes underestimate the peak equipment response. However, there are some similarities between the seismic response analyses and the steady-state response analyses. 23 1. The errors from the uncoupled analyses are most severe when the mass ratios are high. 2. Tuning conditions also affect the accuracy of the uncoupled analyses. In general, the error from an uncoupled analysis is most significant when the equipment is tuned to the fundamental natural frequency of the structure. 3. Damping also influences the accuracy of the uncoupled analyses. As the level of damping of the subsystems increases, the uncoupled seismic analyses become more accurate. These similarities arise because the mass ratios, the damping ratios, and the tuning conditions determine the strength of coupling between the subsystems. 2.3.4 M E T H O D OF SHAKE T A B L E TESTING In conventional methods of shake table testing.the equipment system is subjected to a prescribed simulated seismic motion which closely resembles the predicted building floor motion corresponding to the floor response spectrum supplied by the design requirement. This response spectrum is called the required response spectrum (RRS) and it is often obtained from a non-interactive study of the housing structure alone; therefore, the effect of equipment-structure interaction is usually neglected in conventional practice. Shake table testing can diTectly include the effect of equipment-structure interaction. This is achieved by building and testing a scaled model of the combined equipment-structure system. The scaled parameters are the dimensions, mass and stiffness of the actual system. However, this approach has several major drawbacks. First, the physical size of the shake table dictates the allowable physical parameters of the testing model; this limits applicable equipment-structure systems. Second, the use of scaled models presents added difficulties when fragility testing is of interest In such tests the maximum capacity of the 24 equipment system is determined for both single and multiple frequency waveforms that comply with future requirements. Fragility tests are conducted until the specimen fails and the failure mechanism may be operational or structural. But operational type failure is not obtainable from dynamic testing of scaled models. Moreover, the monitored structural failure may not be valid because the scaled model and the actual equipment system may have different ultimate strengths since such models cannot accurately represent all aspects of the prototype. Finally, perhaps the most important drawback of scaled model testing originates from the difficultes in accurately modelling the nonlinearities and the damping of the actual system. Without an accurate model, the results of shake table testing are not valid. A more practical method uses a combined analytical and experimental approach. A mathematical model of the combined system is needed for the analysis. This involves the determination of the various system parameters, such as the mass, damping, and stiffness of the individual uncoupled subsystems. The system parameters can be obtained from static testing or from the known physical characteristics of the subsystems. These parameters can be checked and upgraded from the performance of dynamic tesung. This process is known as system identification. Knowing the system parameters, an accurate mathematical model of the combined system can be built. For a given ground motion, which closely resembles the predicted ground motion producing the response spectrum provided by the design requirements, a floor motion response record at the point of equipment attachment can be generated from a time step analysis of the combined equipment-structure system. This floor motion record, which includes the effect of equipment-structure interaction, can then be used as the input to the shake table to test the actual equipment system. 25 Another similar approach entails the determination of a new floor response spectrum as defined by some of the authors mentioned in section 1.2.1 such as Sackman and Kelly [15], Villaverde and Newmark [17], Der Kiureghian, Sackman and Nour-Omid [18], and Der Kiureghian and Asfura [11]. This new floor response spectrum, which also includes the effect of equipment-structure interaction, can be generated directly from the ground response spectrum rather than a time step analysis. The system parameters and a^  prescribed ground motion spectrum are again needed in order to obtain the new floor response spectrum. One can then use simulated seismic motions, which closely resemble the motions corresponding to the new response spectrum, as a shake table input to the actual equipment system. 26 - U n c o u p l e d A n a l y s i s - Coupled A n a l y s i s o •—to z-o; (X /A 1 1 1 r— r 0.0 0.2 0.4 0.6 0.6 1.0 FREQUENCY (HZ) — 1.2 1.4 1.6 1.6 2.0 Figure 2.2 Plot or uncoupled and coupled equipment response transfer function, e = 0.01 5. = 0.01 t U n c o u p l e d A n a l y s i s C o u p l e d A n a l y s i s 1 1 1 1 = * V J T — 0.0 0.3 0.4 0.6 0.8 1.0 1.2 FREQUENCY (HZ) 1.4 —I 1.6 1.6 2.0 Figure 2.3 Plot of uncoupled and coupled equipment response transfer function. e = 0.001 5j = 0.01 T i X " a: £si-a 0.0 Peak of coupled analysis U n c o u p l e d A n a l y s i s — — C o u p l e d A n a l y s i s —i 1 1 f» • P 0 2 0.4 0.6 0.8 1.0 |'2 FREQUENCY (HZ) ' ~ l — 1.4 —I— 1.6 ~ l 1.8 2.0 Figure 2.4 Plot of uncoupled and coupled equipment response transfer function. e = 0.0001 5. = o.oi Uncoupled Analysis Coupled Analysis V 1 1 — f - ^ — I 1 1 ^ f— 0.0 0.2 0.4 0.S 0.8 1.0 1.2 1.4 FREQUENCY (HZ) ~1— t.e 1.6 2.0 Figure 2.5 Plot of uncoupled and coupled equipment response transfer function. e = 0.01 $ . = 0.02 T 3 o a " " UJ fM R R P-—|— 0.2 U n c o u p l e d A n a l y s i s — C o u p l e d A n a l y s i s i i !. i ! it ii i \ V i •• • i - i 1 1 ~-r~- \-0.4 0.6 0.8 1.0 1.2 1.4 1.6 FREQUENCY IHZ) 0.0 1.8 2.0 Figure 2.6 Plot of uncoupled and coupled equipment response transfer function, e = 0.001 5 j "= 0.02 o I— u tl_«l a " (jj u. tn P e a k o f c o u p l e d a n a l y s i s i t "1 j . . . ._ U n c o u p l e d A n a l y s i s I! — _ Coupled A n a l y s i s H M l! !! i 1 11 1 1 1 I 1 j 1 1 / / / / 1 I •v. 1 • I 0 0.2 0.4 0.6 1 1 0.6 1.0 F R E Q U E N C Y 1 — • 1.2 ( H Z ) •-•-»t— r - — 1 4 1.6 1.6 2.0 Figure 2.7 Plot of uncoupled and coupled equipment response transfer function, e = 0.0001 $. = 0.02 t-J _ U n c o u p l e d A n a l y s i s — — C o u p l e d A n a l y s i s <_> LU u. <n z a " i \ a v. /»' I! / I I 1 • T — 1 1 1 I ~ 1 0.0 0.2 0.4 0.6 0.6 1.0 1.2 t .4 1.6 1.6 FREQUENCY IHZ) 2.0 Figure 2.8 Plot of uncoupled and coupled equipment response transfer function. e = 0.01 5j = 0.05 p-(_> LJ 01 I-w-i I ; ! U n c o u p l e d A n a l y s i s — C o u p l e d A n a l y s i s w V 1 1 1 1 1 1 1 i"~ ——r— 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 FREQUENCY (HZ) 2.0 Figure 2.9 Plot of uncoupled and coupled equipment response transfer function. f = 0.001 5 = 0.05 4i * 3 -o (_) z =3o a: u. 01 I*? 0.0 Peak of coupled analysis 11 / t i \ i l • / w U n c o u p l e d A n a l y s i s C o u p l e d A n a l y s i s — i — 0.2 —1 0.4 i 1 1 r— 0.8 0.8 1.0 I 2 FREQUENCY (HZ) ' - i f——• 1 , 1.8 1.8 2.0 Figure 2.10 Plot or uncoupled and coupled equipment response transfer function, e = 0.0001 "J = o.05 n \ I; V H Ii I! n i > i • U n c o u p l e d A n a l y s i s C o u p l e d A n a l y s i s I i w \ \ \ N 1 1 1 1 1 1 1 1— 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 FREQUENCY (HZ) i— 1.8 2.0 Figure 2.11 Plot of uncoupled and coupled equipment response transfer function. e = 0.01 5j = 0.10 i! \ Ii \ /' * i; I; Uncoupled AnalyBiB Coupled A n a l y s i s /•' i I, Ii I! \-t W I, Ii / / V W V I — 1.6 I — 1.8 0.0 I — 0.2 -1 1 1 1 1— 0.4 0.6 0.8 1.0 1.2 FREQUENCY [HZ) i 1.4 2.0 Figure 2.12 Plot of uncoupled and coupled equipment response transfer function. e = 0.001 $j = 0.10 'A i, /; U n c o u p l e d A n a l y s i s C o u p l e d A n a l y s i s I. I i /• I I, l'< I! i V. V. \ \ w \ - 1 — 0.2 -1 1 1 1 1— 0.4 0.6 0.8 1.0 1.2 FREQUENCY (HZ) — i — 1.6 - 1 — 1.8 0.0 1.4 2.0 Figure 2.13 Plot of uncoupled and coupled equipment response transfer function, c = 0.0001 $i = 0.10 39 s3 s3 ms2 ks2 m s l k s l M K C e e e 7 7 ^ Figure 2.14 Three D.O.F. structural system. Figure 2.15 S.D.O.F. equipment system. e e e s3 s3 ms2 ks2 m . k . si si Figure 2.16 Combined equipment structure system. Analysis type Structure tuning mode Uncoupled Coupled First — — Uncoupled Coupled Second - — - Uncoupled Coupled Third 2 % Damping in structure 5 % Damping in equipment r— , , , — 1 0.00001 0.0001 0.001 0.01 0.1 Mass Ratio Figure 2.17 Coupled and uncoupled peak response of tuned secondary system to El Centro earthquake N-S component G 15 10 \ Analysis type Structure tuning mode Uncoupled Coupled First — — Uncoupled Coupled Second Uncoupled Coupled Third \ 2 % Damping in all modes \ \ . . . . \ . ... \ \ \ \ 000001 0.0001 0.001 0.01 0.1 Figure 2.18 Coupled and uncoupled peak response of tuned secondary system to El Centro earthquake N-S component Mass Ratio Chapter 3 TESTING FACILITIES AND TEST M O D E L PARAMETERS 3.1 INTRODUCTION One of the main objectives in this study requires the testing and modelling of vibration isolated equipment systems subjected to single axis horizontal base excitation. The tests were divided into two categories a) static testing and b) dynamic testing. This chapter discusses the facilities at the University of British Columbia which were used to undertake these experiments. The properties of the model chosen to represent an actual air handling equipment system, and the properties of the chosen vibration isolation spring units, are also discussed. 3.2 TESTING FACILITIES The static tesung program was performed in the Structural Engineering Laboratory of the Department of Civil Engineering at the University of British Columbia. The static testing program employed a hydraulically actuated piston system which is capable of applying loads simultaneously in two perpendicular directions. This biaxial hydraulic system consists of a 11 kip capacity (MTS model 204.61) and a 50 kip capacity (MTS model 202.01) hydraulic actuator with corresponding loading rigs mounted perpendicular to each other. Each actuator has a load cell (MTS model 661.21A-02 and MTS model 661.22 respectively) incorporated into the arm of the hydraulic ram. The input to each hydraulic ram can either be load controlled or displacement controlled by an MTS control console (MTS model 483.02). This console can also monitor the ram displacements and the applied loads. The results can be displayed on X - Y plotters. The dynamic testing program was performed in the Earthquake Engineering Laboratory of the Department of Civil Engineering at the University of British 42 43 Columbia. This test facility consists of a 3.3 m x 3.3 m earthquake simulation shaking table. It has a capacity to drive a payload of 16000 kg with a maximum acceleration of 2.5 g and a maximum displacement of 130 mm peak to peak. The table is a cellular aluminum structure weighing 2090 kg. It has a grid system of threaded steel inserts into which the test specimens are anchored. The shaking table is driven by an MTS 16000 kg hydraulic actuator and the feedback is controlled by an MTS Systems Corporation electronic servo valve system. Only one-directional horizontal motion is attainable by the table; all other movements are restrained by hydrostatic slip bearings and by supporting columns. A PDP 11/04 minicomputer commands the table by mean of digitized prescribed motion records stored on floppy disks. The records can be actual recorded seismic motions or synthesized motions to comply with the relevant requirements. The minicomputer system also provides data acquisition capability. Analog/digital signal conversion and real time clock functions are controlled by an AR-11 real time module. The analog signals received from the data acquisition instruments are preconditioned by amplifiers. The signals are reconditioned by variable gain buffers and variable cut-off filters before conversion into digital form by the AR-11 real time module. Finally, the digitized records are stored on floppy disks. Up to 16 channels of input signal can be handled by the present system. The sampling rate depends on the total number of channels sampled. The maximun sampling rate is 2 kHz in continuous mode and 35 kHz in burst mode. A communication link between the UBC Amdahl 470 V8 main frame computer and the PDP 11/04 is possible through a H D L C (High-level Data Link Control). This connection allows the user to make use of the powerful Amdahl computer for data storage and analysis. Fig. 3.1 shows the overall function of the shaking table system. The test instrumentation consisted of accelerometers and linear voltage differential transducers (LVDT). Five accelerometers were used to monitor the 44 acceleration of the specimen and the table. A KisfJer servo accelerometer with a range of ± 5 0 g was mounted on the table. Equipment mounted accelerometers were of the strain-gauge type with a range of ± 2 0 g. They were located on the specimen in order to monitor the horizontal, the vertical, and the rotational acceleration of its centre of mass. Three LVDTs were connected to the specimen at the level of the isolator springs. They measured the horizontal and the vertical displacements of the vibration isolator springs relative to the table. The table displacement was measured with a fourth LVDT. Fig. 3.2 shows the mounting position of the instruments. The instrumentation was calibrated periodically to ensure accurate and consistent results. 3.3 TEST M O D E L CHARACTERISTICS An air handling system, Chicago Vane Axial Fan, Design 34 size 54 1/4 arrangement 9 with a 405 US 100 bhp motor frame and inertia base, mounted on vibration isolators was modelled for the purpose of the expermental study. In order to accurately represent the prototype system, the model characteristics must be chosen carefully. In the present study there were several idealizations adopted. 1. The actual air handling fan unit was assumed to be very rigid as compared to the relatively soft vibration isolators, i.e., the fan unit was assumed to undergo rigid body motion under seismic excitation. 2. The mass of the fan unit was assumed to be evenly distributed over its volume. 3. The response of the air handling system under its normal operation was considered to be negligible compared with its seismic response. 4. The vibration isolators were assumed to meet standard specifications for noise and vibration control. Based on the above assumptions, characteristics such as mass, mass moment of inertia, dimensions, location of centre of mass, location of vibration isolators, and types of vibration isolators were chosen to model the actual system. 45 From the specification of the Chicago vane axial fan, the mass of the actual system without an inertia base was estimated to be 1123 kg and the location of the centre of mass was estimated to be 1.367 m above the vibration isolators. Fig. 3.3 shows the location of the basic components of the air handling system. The isolator spacing along the horizontal axes xx and zz was 1.55 m and 0.9 m respectively. Standard specifications for noise and vibration control specify that the minimum distance between the adjacent corner isolators should at least be equal to the height of the centre of gravity of the equipment Since this specification cannot be met with the normal arrangement an inertia base with a mass of 730 kg was constructed for the air handling system. The addition of an inertia base lowers the centre of gravity of the system which reduces the rocking motion and eliminates excessive movement during equipment start up. However, with the additional mass, a different set of vibration isolators are needed, which results in a different dynamic system. Therefore, there are major differences between the seismic response of equipment systems with and without inertia bases. This point will not be examined further in the present study. The total mass of the air handling system with an inertia base was estimated to be 1853 kg and the location of the centre of mass was estimated to be 0.95 m above the vibration isolators. The isolator spacing along the horizontal axes xx and zz was the same as before, 1.55 m and 0.9 m respectively. The mass moment of inertia about the horizontal axes xx and zz was estimated to be 1772.4 kg.m.m and 1507.7 kg.m.m respectively. The above characteristics of the actual system had to be modelled in the study. An inverted concrete box mounted on an inertia base was constructed to model the actual system. This model was designed to withstand an acceleration of 2 g's in all directions. The concrete box, shown in Fig. 3.4, has five components: a top piece and four side pieces. The perimeter of each component is encased by steel angles. 46 The top panel is 0.152 m x 0.953 m x 0.889 m. Two of the side panels are 0.076 m x 1.054 m x 1.499 m and the other two side panels are 0.051 m x 0.737 m x 1.499 m. The five pieces were assembled by welding along the perimeters of the panels. Two A325 M10 threaded rods are implanted longitudinally into each 0.076 m thick side panel at a distance of 0.08 m from the base and 0.06 m apart. They are used to connect the box to its supporting angles and inertia base. The two 1.1 m long supporting angles are 200 x 200 x 200 mm. The inertia base consists of two 2 m long W410 x 39 supporting I beams and a 0.11 m x 0.88 m x 1.91 m concrete base. The I beam was needed to raise the centre of mass of the system to the correct height These components are bolted together to form the final model. The characteristics of the model were estimated to be close to the actual system; this was confirmed later in the system identification part of the experimental study. The selection of the vibration isolators was based on the estimated fan speed of 886 rpm under normal operation. Assuming the damping to be negligible and defining the transmissibilty T^ as the ratio between the maximum force transmitted to the floor and the maximum applied force due to the rotating fan, the vertical stiffness required for the vibration isolators was: K = (27rf(l/T r +l) 1 / 2 /60) 2 m/4 ' (3.1) where f= driving frequency =886 rpm m= total mass of the system Assuming T r can vary between 5% and 20%, the vertical stiffness of the vibration isolators can vary between 209 kN/m and 730 kN/m. Based on the above results, two 47 types of isolators were selected for testing: an open spring isolator with a uni-directional restraint, as shown in Fig. 3.5 (Airmaster ASL-2-208) and an elastomeric isolator as shown in Fig. 3.6 (Korfund series F F D D 1800). These isolators are commerically available and their respective rated vertical stiffness were 630 kN/m and 350 kN/m. I n p u t O u t p u t C A L C O M P P l o t t e r T e l e t y p e F l o p p y D i s k l< 1 , 1 ! i i i if I H D L C A M D A H L 1 7 0 VB Computer A c t u a t o r a n d F o u n d a t i o n R e s t r a i n t Data Acquisit ion Do to Process ing _ Generat ion of Command Signal Figure 3.1 Overall function of U.B.C. shaking table. End view Side view Accelerometers • O 0.89 m B Accelerometer 5 o LVDTs 0 0.78 m 1.03 m Accelerometer ± 77IT „ 0.78 m LVDTs 2.00 m Figure 3.2 Experimental model showing instrument locations. Figure 3.3 Prototype air handling system - Chicago vane axial Tan (Design H sl;e 54 1/4. arrangement 9). o Figure 3.4 Experimental model. 52 Figure 3.6 FJastomeric isolator Chapter 4 EXPERIMENTAL STUDY 4.1 INTRODUCTION The experimental program was divided into two parts: i) static testing and ii) dynamic testing. The static testing program was conducted prior to the dynamic testing program to confirm and explore experimentally the characteristics of the experimental model. The model characteristics of interest were the axial and lateral stiffness of the chosen isolators, the mass of the model and the location of centre of mass of the model. Without the above information, an analytical model of the vibration isolated system, based solely on the commercially rated isolator stiffness and the calculated inertia properties of the system, may be erroneously constructed. The dynamic testing program was subsequently performed to further identify experimentally the characteristics of the vibration isolated equipment model and to investigate experimentally the behaviour of the model under seismic excitation. Based on the above results, an accurate analytical model of the vibration isolated equipment system can be developed to predict its response in a seismic environment The present chapter describes the testing methods and procedures adopted in the experimental program. 4.2 STATIC TESTING PROGRAM Two types of isolators were tested in the static testing program: i) elastomeric isolators and ii) open spring isolators with uni-directional motion restraints. Both types of isolators are commercially available. Only the isolator stiffness in the axial direction is available from the manufacturers; however, the stiffness characteristics in both the axial and the lateral directions are needed to construct the analytical model. The bi-directional hydraulic system described in Chapter three was used in the static testing program. Each type of isolator tested was mounted horizontally in pairs in an inverted 53 54 series as shown schematically in Fig. 4.1; the actual test arrangement is illustrated in Figs. 4.2 and 4.3. The isolator assemblage was initially loaded to a specific axial load using the hydraulic load control system, and then laterally displaced using the hydraulic displacement control system. This procedure was repeated for a range of axial preloads so that the lateral load-deflecton relationships for each pair of isolators, at various levels of axial preloads, were obtained. The lateral stiffness of the isolators was determined from these lateral load-deflection curves. Also, the effect of axial load on the lateral stiffness of these isolators was found. With the same experimental set up, the isolator assemblage was brought to an initial lateral preload using the load controlled hydraulic system. Maintaining the lateral preload, the isolators were displaced axially by the displacement controlled hydraulic ram, which yielded an axial load-deflection curve having the characteristics indicated in Fig. 4.4. Repeating the above procedure for a range of lateral preloads resulted in axial load-deflection curves for the isolator assemblage at various levels of lateral preloads. Again, from the load-deflection curves, the axial stiffness characteristics of the isolators were obtained for comparison with the commercially rated values. Moreover, the influence of lateral load on the axial stiffness characteristics was obtained. The mass of the model was measured by suspending its individual components from a 10 kip capacity MTS load cell. The location of the centre of mass of the model was determined by suspending the model along one of its top edges from a crane. A ' plumb line was dropped from the point of suspension. This procedure was repeated by suspending the model along the opposite top edge and the intersection of the plumb lines yielded the location of centre of mass. of the model. Finally, based on these results, the mass moment of inertia of the model about the centroidal axis was calculated. 55 4.3 DYNAMIC TESTING PROGRAM The dynamic testing program was conducted after the completion of the static tesung program. The first series of dynamic tests were exploratory in nature. They consisted of low-level sinusoidal testing, low-level sinusoidal decay testing, and random white noise testing. From these tests, the system characteristics such as the natural frequencies, the mode shapes, and the modal damping ratios were' found. The next series of dynamic tests were seismic motion tests. The seismic response and the failure modes of the model were obtained from these tests. 4.3.1 SINUSOIDAL TESTING The sinusoidal tests were performed at relatively small fixed base displacement amplitudes. The magnitude of the amplitude was limited to avoid the possibility of damage to the model at system resonance. As a result, the base motion displacement amplitudes for the elastomeric isolators and the open spring isolators with uni-directional restraints were 1.8 mm and 1.3 mm respectively. Based on a preliminary estimation of the natural frequencies of the model from the static testing, a range of forcing frequencies (0 - 15 Hz) was chosen to excite the model. Approximately 25 frequencies were used within this range. Low pass filters with a cut off frequency of 10 Hz were used at test frequencies below 10 Hz, and 40 Hz low pass filters were used at test frequencies above 10 Hz. The low pass filters eliminated extraneous instrumentation noise which might otherwise have affected the accuracy of the results. A sampling duration of 10 seconds was chosen to ensure that the model responded in a steady state vibration. Digitized time history records of the model response and of the table motion were recorded and stored on floppy disks as discussed in Chapter three. 56 The ratio between the steady response amplitude for each degree of freedom of the model and the forcing motion amplitude of the table can be obtained by examining the response records at each forcing frequency. A plot of this response amplitude ratio at each forcing frequency then yielded the frequency response curve for the system. The natural frequencies and mode shapes of the system were obtained from these results. 4.3.2 SINUSOIDAL DECAY TESTING Damping is perhaps the most important and very often the most difficult dynamic characteristic to measure in a vibrating system. The purpose of the sinusoidal decay tests was to determine the modal damping ratios of the system. From the results of the static and sinusoidal testing, a good estimate of the natural frequencies of the system was obtained. In the sinusoidal decay lest, the system was excited at each of its natural frequencies so that it responded by vibrating in the corresponding natural mode. When the table motion was suddenly cut off, the system responded at its natural frequency in a mode of free damped vibration. The resulting decaying response was recorded and the modal damping ratios were determined by the well known logarithmic decrement method. Several assumptions are made in using the logarithmic decrement method: i) the damping of the system is viscous, ii) the system is underdamped, and iii) the different modes of vibration in a multi-degree of freedom system are well separated. The percentage of critical damping in each mode of vibration is given as: D ^ ln(x,/x ) / ( J 2 T T ) (4.1) 57 where D = percentage of critical damping in a mode th Xj= acceleration amplitude of the j cycle of decay This approximate method of determining the damping is valid for systems having a percentage of critical damping in each mode which is less than 20%, which is valid in our case. 4.3.3 R A N D O M WHITE NOISE TESTING Another method of system identification involves the spectral analysis of the system response due to white noise excitation. Random, band limited (0 -25 Hz) white noise signals can be generated by the technique of summation of sinusoids (see Appendix A). A band limited (0 - 25 Hz) white noise signal is defined as a random variable having a constant power spectral density function for frequencies from 0 to 25 Hz. The frequency band of the white noise signal was chosen to span the estimated natural frequency range of the model. When such a record is used as an input to the system, all the frequencies of the system within the range of 0 to 25 Hz are excited with equal intensity. The resulting input/output relationships can be analysed by the techniques of spectral analysis to yield the frequency response functions for the system of interest The theory of system identification using spectral analysis techniques is presented in Appendix B. A computer program, FOUR.S, based on the above principles has been written to analyse the time series records. The computer listing for this program is provided in Appendix C. In the random white noise tests, the model was excited by a randomly-generated band limited white noise signal. The sampling duration of each test was 26 seconds, with a recording interval of 0.006 second. The recording interval 58 was kept as short as possible to ensure good resolution and to avoid aliasing. Ideally, the sampling duration should be kept as long as possible. A practical choice of sampling duration was chosen to minimize possible bias errors or truncation errors in the spectral analysis. The input base motion was scaled to a peak aceleration of 0.1 g on the MTS system, which generated strong model responses without causing damage to the model. Low pass filters with a 40 Hz cut off frequency were used. The model responses and the table motions were recorded and stored on floppy disks. The recorded data were transferred to the main frame computer for analysis to yield the frequency response curves and the phase relationships for the system. The experiments were repeated three times and the resulting frequency response curves and phase angle relationships were averaged to minimize possible random errors. 4.3.4 SEISMIC MOTION TESTING The acceleration records of actual earthquake ground motions were used as input motions to the shaking table to generate seismic responses of the model. The choice of these input motion records was based on their frequency content in comparison to the natural frequencies of the equipment model. The chosen records were: i) El Centro N-S (1940), ii) San Fernando N21E (1971), and iii) Taft (1950). The input acceleration interval was 0.02 second and the sampling interval was 0.006 second. The sampling duration was the same as the input duration of the recorded earthquakes. The low pass filters were set at 10 Hz. The peak acceleration of the input motion was scaled to appropriate levels ranging from 0.1 g to 0.2 g, so that strong model responses were recorded without damaging the model. The model responses were, recorded on floppy disks and later transferred onto the main frame computer for further analysis. 59 After the completion of the seismic ground motion tests, the resulting response time histories were compared with the predicted model response time histories from an analytical study. The method of analysis will be discussed in the next chapter. Once satisfactory agreement between the predicted and measured responses was obtained an actual floor motion record was used as the input motion. The chosen floor motion record was the San Fernando 2500 Wilshire Blvd. roof response N61W (1971). Since the natural frequencies of the structure and the equipment were well separated, and the mass ratio between the equipment and the structure was assumed small, it follows that equipment-structure interaction was not significant in this case. The peak input floor acceleration was gradually increased until the isolators failed and the failure mode was observed and recorded. This represented the end point of all the experimental work in this thesis. The results of the experimental studies are presented in Chapter six. The next chapter describes the analytical method used to predict the seismic response of vibration isolated equipment-structure systems. 60 Cyclic lateral load Fixed axial load -0-Cydic axial load Fixed lateral load Figure 4.1 Schematic representation of bi-directional static test 61 Figure 4.4 Typical measured load - deflection curve Chapter 5 ANALYTICAL STUDY 5.1 INTRODUCTION The dynamic response of vibration isolated equipment housed in a structure subjected to strong motion earthquakes may also be investigated by an analytical approach. This approach requires the formulation of an accurate mathematical model of the dynamic system of interest, for which a set of equations that approximately represent the dynamic characteristics of the system can be written. The solution to this set of equations will yield a prediction of the dynamic response of the system to prescribed ground motions. In- this chapter, the assumptions and idealizations in the mathematical modelling of vibration isolated equipment systems and their housing structures will be presented. A fairly general mathematical model of equipment-structure systems, which can incorporate typical nonlinear properties of vibration isolators, will also be developed. Finally, a time step numerical integration scheme will be employed to solve the set of equations, yielding the dynamic response of vibration isolated equipment systems. 5.2 M O D E L IDEALIZATIONS For the purpose of developing a general mathematical model of vibration isolated equipment-structure systems, it is necessary to make simplifying assumptions and idealizations regarding the dynamic properties of the system. The dynamic properties of interest are the mass, stiffness, and damping parameters of the system. There are two broad categories of system modelling: i) lumped-parameter models and ii) continuous-parameter models. Both models are idealized; in general, lumped-parameter models are more common than continuous-parameter models but continuous parameters models are more accurate. In the present study, the simpler 62 63 lumped-parameter modelling method is adopted. In lumped-parameter models, various system characteristics are lumped into representative elements which either store or dissipate the seismic energy input of the system. Examples of energy storage elements are the mass and spring elements which store the kinetic energy and the potential energy respectively. In contrast, dampers are typical energy dissipating elements. These idealized lumped-parameter elements are assembled mathematically to form the final analytical model. Let us now examine the various idealized lumped-parameter elements in the mathematical modelling of the housing structure. It is commonly assumed that the housing structure of interest is a multistory simple shear building frame. Each floor of the building is assumed to act as a rigid diaphragm in its own plane. In addition, seismic motion is assumed to act in a single horizontal direction normal to a vertical face of the shear building frame. This idealization allows each rigid floor diaphragm to translate in a single horizontal direction only. The second idealization is that the mass of each floor is entirely lumped into the rigid diaphragm. The only inertia forces experienced by the rigid floor diaphragm are the inertia forces in the horizontal direction. Interstory stiffness of the building is provided by the elastic lateral stiffness of the interstory columns. The damping characteristics of the structure are assumed to be classical such that a completely diagonal matrix of generalized damping coefficients may be obtained by pre- and post- multiplying the damping matrix by the modal matrix of the structure. In general, this assumption is valid for structures which have critical damping ratios less than about twenty percent However, as stated in Chapter two, the assumption of classical damping fails for systems which have tuned frequencies and widely varying material properties, such as large differences in the interstory stiffness or mass. With the above assumptions, the idealized lumped parameters of the mathematical model of the housing structure are the mass of each floor, the interstory stiffness, and the damping of the building. 64 The formulation of the mathematical model of the vibration isolated equipment also involves many assumptions. First, the inertia properties, such as the mass or the mass moment of inertia of the vibration isolated machinery are assumed to be lumped entirely at the centre of mass of the system. This assumption dictates that all the inertia forces experienced by the equipment system act through the centre of mass of the system. Moreover, the machinery is assumed to be infinitely rigid as compared to the vibration isolators; therefore, the machinery is modelled as a rigid unit mounted on vibration isolators with the proper inertia properties at the centre of gravity. The forcing base motion is assumed to act in a single horizontal direction normal to the vertical plane of a pair of isolators. The response of the system to such base motion can be completely described by the three components of rigid body motion: the lateral movement, the rotation, and the vertical motion at the centre of gravity of the system. This system has three degrees of freedom as shown in Fig. 5.1. The vibration isolators provide the equipment system with its lateral, vertical, and rotational spring stiffness. Typically, vibration isolators exhibit some form of nonlinear stiffness characteristics. These nonlinearities can arise from i) materially nonlinearity of the isolators and ii) the geometric nonlinearity of the system. For a material nonlineaT system, the displacements and strains are considered small and the nonlinear effect lies in the nonlinear load deformation relationship of the system. In geometrically nonlinear systems, the nonlinearities are the results of large deformations in the system. The analytical model used in this study attempts to model the material nonlinearity of the vibration isolators only. The nonlinear isolator is idealized as a group of discontinuous linear elastic springs as shown in Fig. 5.1. Depending on the nonlinear stiffness characteristic of the vibration isolators, each of the discontinous linear elastic springs in the model has a specific spring constant Therefore, a piecewise linear stiffness characteristic is assumed to model the nonlinear isolator stiffness as shown in Fig. 5.2. The damping characteristic of the equipment system is assumed to be 65 provided solely by the damping in the vibration isolators. Classical damping is again assumed. Modal damping ratios of the equipment system are needed to formulate the mathematical model of the equipment system. These system parameters were obtained from the experimental tests. 5.3 DEVELOPMENT OF T H E MATHEMATICAL M O D E L The idealizations in the mathematical modelling of vibration isolated equipment-structure systems have been presented in the last section. Now, based on those idealizations, a mathematical model of the system will be developed. First, a mathematical model of the housing structure will be presented. Secondly, a mathematical model of the equipment system, which can incorporate the nonlinearity of the vibration isolators, will be studied. Finally, a mathematical model combining the equipment and the structure will be examined. 5.3.1 T H E STRUCTURAL M O D E L The structural model is very similar to the model of the simple three story shear building frame studied in Chapter two. In the present model, however, the system is expanded to a maximum of ten floors. The mass of each floor is denoted bv m . (i = l,2,..,r) where r ^ l O denotes the number of si floors in the structure. The interstory stiffness is given by k s . (i = l,2,...,r). This structural system can have a maximum of ten modes of vibration resulting in a maximum of ten modal damping ratios denoted by $ . (i = l,2,...,r). The equation of motion of the structural system when excited by base motion x (t) is: g M x +C x +K x = - M I 1 (t) s s s s s s s g w (5.1) where x s is the relative displacement vector, x s 66 = [k sj,x ^ x s r l is the relative velocity vector, and x g = t x s j ' X s 2 - — X S P ' s ^ relative acceleration vector of the floors with respect to the moving base. I T = [1,1, ,1] is the influence vector coupling the base motion input to the individual floor. The structural mass matrix M and the structural stiffness matrix s K are: s M m si m s2 m sr (5.2) K k s l + ks2 ~ ks2 -k s2 k -+k , . s2 s3 -k sr -k k sr sr Finally, the structural damping matrix C s is obtained from the structural mass matrix M , the structural stiffness matrix K , and the structural generalized s s damping matrix C g as illustrated in Chapter two. The structural generalized damping matrix and the structural damping matrix are as follows: 67 C = 2 C = M * C s s s structural mode shape matrix th the r natural frequency of the structure the percentage of critical damping of the r^1 mode of the structure Knowing the mass, stiffness, and damping matrix, the equations of motion of the structural system can be formulated. 5.3.2 T H E EQUIPMENT M O D E L The formulation of the mathematical model of the vibration isolated equipment will now be presented. First consider the idealized vibration isolated equipment system shown in Fig. 5.1. The system consists of a ; rigid mass with its inertia properties lumped at the centre of gravity of the model. The rigid mass is supported by four sets of springs which provide the system with its lateral, rotational, and vertical stiffness. Each set of springs consists of a linear elastic spring which provides a restoring force to the system when the rigid mass is displaced from its equilibrium position. The development of excessive equipment seismic motion is 5 si u s l m s l * T * M 5 S S K u> m J sr sr sr (5.3) where <i> = u> = sr *sr = 68 prevented by the provision of snubbers which act as motion restraint devices. Such devices provide an abrupt increase in restoring force to the equipment system when the displacement of the unit exceeds a prescribed value. In order to model the characteristics of these devices, a discontinuous linear spring is added to the assemblage, which provides the required addtional stiffness (see Fig. 5.1). This is a simple spring assemblage modelling a bilinear hardening system. In order to model a system which exhibits degrading stiffness characteristics, the discontinuous linear spring can be assigned a negative stiffness, such that the system softens when the displaced equipment system exceeds a prescribed displacement as shown in Fig. 5.2. From the static testing program, it was discovered that when the horizontal snubber units in the open spring isolators with uni-directional motion restraints were engaged at some Fixed lateral load, the axial suffness of the springs increased dramatically. However, under a continuously increasing axial load, the vibration isolators eventually regained their original axial stiffness as shown in Fig. 5.3. This phenomenon is the result of coulomb damping. With the horizontal snubbers engaged, the springs were held axially by friction between the snubber and the spring assemblage. An increased axial load was needed to overcome this friction force before slippage between the snubber and the spring assemblge could occur in the axial direction; at this point, the spring regained its normal axial stiffness. To model this effect, a friction spring is included in each set of vertically mounted springs in the mathematical model and the friction spring becomes active only when the horizontal snubber units are engaged. Finally, besides the effect of coulomb damping, the effect of material damping of the vibration isolators needs to be considered. Material damping is modelled by the provision of linear viscous dashpot dampers in the mathematical model. 69 The equations of motion of the model to base excitation xb(t) at the point of equipment support is formulated about the centre of mass of the model and with reference to the static equilibrium position of the dynamic system. The equations of motion are given as: M x +C x +K'(x )X = - M IX. (t) (5.4) e e e e ev c' e e b w v ' where x = [x .,x _,x ,] is the relative displacement vector, x = [x ,,x _,x ,] e el e2 e3 v e el e2 e3 is the relative velocitv vector, X = [x ,,X ~,x ,] is the relative acceleration e el e2 e3 vector of the equipment with reference to the base motion X^(t). The three degrees of freedom x ^, x^ , and x ^ are respectively the lateral motion, the rotation, and the vertical movement of the centre of gravity of the equipment system. The sign convention for positive movements is defined as horizontal motion to the right, clockwise rotation, and downward motion. The vector I 7 = [1,1,1] is the influence vector coupling the base input motion to the individual degrees of freedom of the equipment The stiffness matrix is nonlinear because it depends on the equipment displacement response, x . The equipment mass, damping, and stiffness matrices are formulated in the following manner. First, x £ , the relative displacements of the equipment are assumed' to be small. Based on this assumption, second order terms in the equations of motion can be neglected. Also, the linearization assumption sin(xe2) = x£2 and cos(xe.j) = l are assumed to apply. Summing all the forces in the direction of each degree of freedom of the system and setting the result to zero results in the desired equations of motion. So, assuming the system is undamped and excited by base motion X^(t), the inertia and spring forces of the displaced system are shown in Fig. 5.4. A compressive spring force is defined as positive by our sign convention. Equilibruim in the direction of the horizontal, 70 rotational, and vertical degrees of freedom of the equipment results in the following equations of motion: Mx . + C + D=-Mx.(t ) el b v ' Jx £ 2-A(arm a) + b(arm b)-C(arm c)-D(arm d) = 0 (5.5) Mx ,+A+B=0 e3 where M = mass of the equipment system J = mass moment of inerua of the equipment system about the centroidal axis A = total vertical spring force on the left side B = total vertical spring force on the right side C = total horizontal spring force on the right side D = total horizontal spring force on the left side arm a = horizontal distance between the eg. and the vertical spring on the left side arm b = horizontal distance between the e.g. and the vertical spring on the right side arm c = vertical distance between the eg. and the horizontal spring on the right side arm d = vertical distance between the eg. and the horizontal spring on the left side The various distances between the springs and the centre of mass of the displaced system is determined from the geometry and the displacements of the system. Assuming small rotations and displacements, the various distances are as 71 follows: arm a = b. + ax . 1 e2 arm b = b -ax 0 r e2 arm c = a + b r x e2 arm d = a-b,x e 2 where a = vertical distance between the eg. and the undisplaced vibration isolators bj = horizontal distance between the e.g. and the undisplaced vibration isolators on the left b r = horizontal distance between the eg. and the undisplaced vibration isolators on the right At this point, let us examine the restoring spring forces in the equations of motion. These spring forces may be nonlinear because of the discontinuous springs. First, let us assume that the system remains in its linear state. And let us denote K - ^ ^ , as the initial lateral stiffness of the vibration isolator on the left and the right side of the system respectively, and K - j y L » ^ J V R A S the initial vertical stiffness of the vibration isolator on the left and right side of the system respectively. From the geometry of the system, the various displacements of the vibration isolators due to the movements at the centre of gravity of the system are given by equation (5.7). x c = horizontal displacement of the vibration isolators = x .-ax . el e2 72 Vj = vertical displacement of the left vibration isolator = " b l x e 2 + xe3 ( 5 J ) V = vertical displacement of the right vibration isolator = b r x e 2 + x e 3 From the sign convention, positive displacements mean shortening of the springs. The linear restoring spring forces as a result of the various displacements are: A = K 1 V L V 1 = K l V L ( - b I X e 2 + X e3 ) B = K l V R V r = K l V R ( b r x e 2 + x e 3 ) C = K l H L X c = K l H L ( x e r a x e 2 ) D = K l H R X c = K l H R ( x e r a X e 2 ) (5.8) Substituting equations (5.6) and (5.8) into equation (5.5) yields the equation of motion for the linear three degree of freedom system. Let us now examine the situation in which the horizontal discontinuous spring units are engaged. Let the clearance between the mass and the horizontal discontinuous spring units be denoted by c and the stiffness of the horizontal discontinuous springs by K ^ H ' " ^ n e s e discontinuous springs are engaged when the absolute value of the horizontal displacement of the vibration isolators x £ exceeds the clearance c. The Final horizontal stiffness of the left and the right spring units are ^2HL ~ K 1 H L + K 2 H a n d K 2 H R = K 1 H R + K 2 H r e s P e c u v e l > - T n e horizontal spring forces of the system with the horizontal discontinuous spring units engaged are: C = K 2 H R X c - c K 2 H ( s g n ( x c ) ) ( 5 9 ) D = K 2 H L X c - c K 2 H ( s g n ( X c ) ) 73 where sgn(xc) = the algebraic sign of x c The expressions for the vertical restoring spring forces remain unchanged since they are still linear. Again, the expressions for the spring forces and the spring distances can be substituted into equation (5.5) to yield .the equations of motion for the case where only the horizontal discontinuous springs are engaged. The situation in which only the vertical discontinuous spring on the left side of the model is engaged will now be studied. In this situation, there are two possible cases: i) the discontinuous spring is engaged in a tensile mode when the left side of the model has lifted beyond the static equilibrium position, and ii) the discontinuous spring is engaged in a compressive mode when the left side of the model has dropped beyond the static equilibrium position. Some vibration isolators have very different tensile and compressive characteristics. Since the equations of motion are written about the static equilibrium position, and the isolators are initially in compression due to the dead weight of the machinery, such vibration isolators may experience an abrupt change in stiffness when the lifting vibration isolator overcomes the initial compressive dead load of the machinery. Moreover, some open spring isolators may not provide any motion restraining devices against uplift, whereas their compressive stiffness may increase abruptly when the spacing between the spring coils is closed in a compressive mode. In order to model the above situations, the vertical clearance, measured from the static equilibrium position, between the vertical spring and the rigid mass during an uplift is allowed to take on a different value from that during a downward movement. Also the tensile and compressive stiffness of the discontinuous vertical springs may be assigned different values. 74 Let us denote the uplift . clearance between the discontinuous spring and the rigid mass by D j and the downward clearance by D^- Keeping in mind that a positive vertical displacement is defined as downward, the discontinuous vertical spring on the leftside is engaged in tension when Vj is less than - D j and it is engaged in compression when Vj is greater than D^. Let ^2\2 denote the compressive suffness of the discontinuous vertical springs and Kjyi denote the tensile stiffness of these springs. The total vertical stiffness of the springs on the left side of the model is denoted by K^vx, w m c n equals * M V L + ^ 2 V 1 m a n uP^'ft mode and ^ 2 Y L + ^ 2 V 2 m a downward mode. The new vertical restoring spring forces for the case where the vertical discontinuous spring units on the left side are engaged in a downward mode and in an uplift mode are given by equation (5.10) and (5.11) respectively. A = V , K 2 V L - D 2 K 2 V 2 = ( - b j X ^ + x ^ X K j ^ + K ^ h D ^ ^ (5.10) A = V j K ^ - f D j K ^ = (-b 1 x e 2 + x e 3 X K 1 V L + K 2 V i ) + D 1 K 2 v l ( 5 . 1 1 ) Repeating the above procedure for the situation when only the vertical discontinuous spring on the right side is engaged, yields the following results. When is less than ~^>y the vertical discontinuous spring on the right side is engaged in tension. When V r is greater than D 2 , the vertical discontinuous spring on the right side is engaged in compression. The total vertical stiffness of the springs on the right side of the model is denoted by K- 2 v-^, which equals F£ 1 V R + ^ 2 V 1 in an uplifting mode and ^ • J Y R + ^ 2 V 2 m a n downward m ° d e . The new restoring spring forces for the case where the vertical discontinuous spring unit on the right side is engaged in a downward mode and in an uplift mode are given by equation (5.12) and (5.13) respectively. 75 B = V r K 2 V R - D 2 K 2 V 2 = <brXe2 + X e 3 ) ( K l V R + K 2 V 2 > " D 2 K 2 V 2 ( 5 J 2 ) B = V K r K 2 V R + D l K 2 V l = <brXe2 + X e 3 ) ( K l V R + K 2 V 2 > + D l K 2 V l ^ 1 3 ) Depending on which nonlinear springs are activated, equation (5.9) -equation (5.13) and equation (5.6) can be substituted into equation (5.5) to yield the equations of motion for the nonlinear three degree of freedom equipment system. Finally equation (5.5) can be rearranged into matrix form to yield the mass and stiffness matrices in equation (5.4). Let us now formulate the equipment damping matrix C g from the equipment mass matrix M e > the equipment stiffness matrix K £ , and the equipment generalized damping matrix C g as illustrated in Chapter two. The equipment natural frequencies u . (i = 1,2,3) and the equipment mode shape matrix $ are obtained from an eigenvalue solution of the frequency equation for the undamped free vibrating equipment system. The equipment generalized damping matrix and the equipment damping matrix are as follows: C = 2 e C = M * C *<i> T M e e e e e e * e2 w e2 M e2 S e 3 w e 3 M e 3 (5.14) where th $ £T=percentage of critical damping of the r mode of the equipment Knowing the equipment mass, stiffness and damping matrix, the equation of motion of the equipment system can be formulated. 76 So far, the treatment of the friction spring has been omitted in the analysis. Let us focus our attention on the mathematical formulation of these springs. The laws of Coulomb are assumed to govern and the friction force is assumed to be proportional to the normal force acting between the sliding surfaces. In the case of the open spring isolators with uni-directional restraints, the normal force acting between the sliding surfaces is the lateral load acting on the snubber units when the snubber units are engaged. The friction force always acts in a direction opposite to the velocity of the system, as in the case of the dashpot damper. However, the magnitude of the frictional resistance is independent of the velocity. The frictional resistance is formulated as an extra vertical spring load, which is proportional to the lateral spring forces C and D. The constant of proportionality u is determined from the experimental results. The vertical spring forces on the left and right side of the model, when the horizontal discontinuous springs are engaged, are given by equation (5.15) and equation (5.16) respectively A = K ^ V j + s g n K M A x ^ ^ V j ) (5.15) B = K 1 V R V r + s g n ( ( u A x c K 2 H ) , A V r ) (5.16) where A x £ = xc(t+At)-xc(t) = change in the lateral spring displacement ^ x c ^ 2 H = s n u ' 3 ^ e r s P n n § f ° r c e A V j = Vjd+AO-Vjd) = change in the vertical left spring displacement A V r = V r(t+At)-V r(t) = change in the vertical right spring displacement 77 Substituting equation (5.9) - equation (5.16) and equation (5.6) into equation (5.5) for the appropriate cases yields the complete system of nonlinear equations of motion for the vibration isolated equipment system. 5.3.3 T H E COMBINED EQUIPMENT STRUCTURE M O D E L The mathematical modelling of the combined equipment-structure system will now be presented. The equations of motion for such a system under a ground excitation Xg(t)is as follows: Mx + Cx+K(x)x=-MIx (t) (5.17) Again x, x, X are respectively the relative displacement, velocity, and acceleration of each degree of freedom to the ground and I is the influence vector coupling the input ground motion to each degree of freedom. The structure of interest will be considered to have r stories (r ^ 10); the overall system matrices depend upon the location of the equipment system in the building. Assume that the vibration isolated equipment is attached to the i^1 floor of an T story building. The various system matrices can be found in the usual fashion, with the mass matrix of the system given as: 78 m, m i + 1 M = m. i + 2 (5.18) m. i + 3 m where .th m. = mass of the i floor 1 m. + j — ni. + 2 = mass of the equipment system m i + 2 = m a s s m o m e n t 0 I " i n e r u a ° f the equipment system The suffness matrix of the combined system is formulated by the same method and is given by: K = k, + k2 . k. + k„ -kn -k 1 2 — - k n kn kn k -k 2 1 k 2 a k„ k k3. k 3 2 k -k -k k r r (5.19) 79 where k. = i ^ 1 interstory suffness of the structure k.j = the equipment restoring force corresponding to the i 1 * 1 equipment degree of freedom due to unit displacement of the j ^ 1 equipment degree of freedom (ij = 1,2,3) The various equipment stiffness entries kjj in the combined stiffness matrix can be obtained directly from the suffness of the vibration isolated equipment system formulated in the previous section. Finally, the combined equipment-structure damping matrix can be formulated in a similar fashion and is given by: C = c +c n - C n "Cn -c - C n Cn Cn Cn - C 2 1 Cn Cn c23 -c 3 1 CJI c32 c33 rr (5.20) M _ entries from the structural damping matrix Cy = entries from the equipment damping matrix The technique used to solve the resulting nonlinear system of equations of motion will be discussed in the next section. 80 5.4 M E T H O D OF SOLUTION The method of solution of the nonlinear system of equations of motion is the well known Newmark-Beta method. The solution scheme is similar to that discussed in Chapter two. The major difference between the present and the former scheme is the treatment of the nonlinear stiffness terms. In order to solve the nonlinear system of equations of motion by this step-by-step integration method, the equations must be expressed in incremental form: MAX(t)+CAx(t) + K(x)Ax(t) =-MIAx (t) (5.21) s where A represents an incremental change in value within a time step At. The nonlinearity in equation (5.21) exists because the stiffness matrix is dependent upon the displacement solution vector x. However, within a time step At each nonlinear spring can be assumed to have a constant stiffness. The stiffness of each nonlinear spring is assumed to be the tangent stiffness of the spring at the beginning of a time increment Therefore, the stiffness matrix only needs to be examined for modification at the conclusion of each time increment Then the incremental equations of motion within a time increment can be considered to be linear. Either a constant or a linear variation in the change of acceleration response over a small time interval can be adopted in the Newmark-Beta method. The resulting expressions of the response from the application of the numerical integration scheme are given by equation (2.27). The above time step integration scheme may lead to erroneous results if the chosen time step is too laTge. In addition, the nonlinear springs are modelled by a piecewise linear scheme which means that a change in stiffness can occur abruptly. When an abrupt stiffness change occurs, the calculated responses are incorrect because they have been computed on the basis of linear stiffness within a time step. Moreover, dynamic equilibrium will not generally be satisfied at the conclusion of the 81 lime step. The solution scheme adopts two strategies to reduce the possible error caused by these conditions: i) equilibrium correction and ii) event to event integration scheme. In the equilibrium correction scheme, the dynamic equilibrium is checked at the conclusion of each time interval. Any unbalanced loads are treated as external loads and are added to the load increment for the next time step. Very small time increments must be used for this method, otherwise substantial departure from the real load-displacement path is possible. The alternative scheme subdivides the time increments into sub-increments when an abrupt change in stiffness, an event, occurs. By using an iterative technique, the solution can proceed from event to event. There is never an unbalanced load in this solution scheme and the solution follows the exact load deformation relationship. However, the event to event method has the disadvantage of being computationally much more expensive than the equilibrium correction scheme. In the present study the equilibrium correction scheme is used during every time increment The event to event method is only employed when the horizontal discontinuous springs are engaged. This concludes the discussion on the analytical study of the seismic response of the vibration isolated equipment-structure system. A computer program, EVIS.S, based on the above principles has been written to analyse the seismic response of such systems. The computer listing for this program is provided in Appendix C. Discontinuous spring \ '2H 1HL \ -KNH -mm Discontinuous spring 2V1 (tens.) K.JY2 (comp.) I 1VL 1VR 1HR T a L - r D 2 Friction spring Discontinuous springs H 4 Figure 5.1 Idealized model of the vibration isolated equipment system. 83 Experimental M o d e l Figure 5.2 Nonlinear load- deflection curve (static test) and piecewise linear load- deflection curve (analytical model) Figure 5.3 Schematic representation of coulomb friction between isolator and snubber device. Figure 5.4 Equilibruim forces in a displaced vibration isolated equipment system. Chapter 6 EXPERIMENTAL AND ANALYTICAL RESULTS 6.1 INTRODUCTION The experimental and analytical results obtained in this study are presented in this chapter. The results are divided into two categories: i) model identification test results and ii) random excitation test results. As described in Chapter five, a mathematical model of the physical system of interest has been developed. The physical parameters of the mathematical model are based on the model identification test results obtained from the static and the dynamic testing program. Also, a computer program utilizing time step analysis has been developed to generate numerical solutions of the system of dynamic equations of motion. In this chapter, the model identification test results and the physical parameters adopted in the mathematical model are presented first. Then the random excitation test results are compared with the numerical solutions of the model equations of motion subjected to the same random base motion inputs. Finally the performance of the two vibration isolated equipment systems are discussed. 6.2 M O D E L INDENTIFICATION TEST RESULTS (STATIC) The static test program yielded the inertia properties of the experimental model and the stiffness characteristics of the vibration isolators. Table 6.1 shows the measured inertia properties of the model and the calculated inertia properties of the prototype equipment system. The prototype equipment system is the Chicago Vane Axial Fan Design 34 size 54 1/4 arrangement 9 with a 405 US 100 bhp motor frame and inertia base. The idealizations and the characteristics of this system have been discussed in Chapter four. 85 86 Mass Mass moment of inertia Location of eg. w.r.t z-z axis above base Model 2136.6 kg 1510.8 kg.m.m .965 m Prototype 1852.9 kg 1772.4 kg.m.m .950 m Table 6.1 Model and prototype inertia properties. A fairly close match can be seen between the model and prototype inertia properties. Exact matching of the model inertia properties with the prototype inertia properties proved to be a difficult task. However, this difference is not significant as the aim of the present study is to observe the seismic response of typical vibration isolated equipment systems which include a wide range of inertia characteristics. The prototype equipment system serves to provide the inertia characteristics of one of these systems. As long as the model inertia properties are reasonably close to the prototype inertia properties, and they do not exceed the dead load capacity of the selected vibration isolators, the model is considered acceptable. Instead of modelling exactly the particular prototype equipment system, our experiment models a similar vibration isolated equipment system which possesses the measured inertia properties summarized in Table 6.1. The lateral and axial stiffness characteristics of the elastomeric isolators and the open spring isolators with motion restraint were obtained from the static testing program. Fig. 6.1 shows the axial load-deflection curves of the elastomeric isolators under three different constant lateral preloads (0, 1.1, 2.2 kN). The following conclusions can be drawn from these results: 1. The elastomeric isolators have nonlinear axial stiffness characteristics. The isolator tends to stiffen under an increasing axial load and the loading path differs 87 substantially from the unloading path. This indicates that high hysteretic damping characteristics can be expected from the elastomeric isolators. 2. The effect of lateral preload on the axial stiffness characteristics is not significant, at least in the preload range examined. The axial load-deflection curves are similar and consistent under the various lateral preloads imposed. 3. The measured axial stiffness value differs substantially from the commercially rated axial stiffness value. Under an axial load of 5240 N (1200 lb), the approximate dead load of the equipment model supported by a vibration isolator, the measured tangent stiffness was 1350 kN/m (7700 lb/in); the commercially rated axial stiffness is 635 kN/m (3600 lb/in). The lateral stiffness characteristics of an elastomeric isolator under various constant axial preloads are shown in Fig. 6.2 For sake of clarity the complete unloading branch of only one test is shown in this figure. The results indicate the following: 1. The loading path differs from the unloading path; therefore, the lateral stiffness characteristics are nonlinear. Hysteretic damping characteristics can be expected . 2. The influence of axial preloads on the lateral stiffness characteristics of the elastomeric isolator is significant The lateral stiffness tends to increase with increasing axial preload. The lateral stiffness during loading of one elastomeric isolator varied from 232 kN/m (1333 lb/in) to 800 kN/m (4570 lb/in) when the axial preloads were 2.2 kN (500 lb) in tension to 18 kN (4000 lb) in compression respectively. The lateral stiffness during unloading of one elastomeric isolator varied from 455 kN/m (2600 lb/in) to 1750 kN/m (10000 lb/in) under axial preloads of 2.2 kN (500 lb) in tension to 18 kN (4000 lb) in compression. Fig. 6.3 shows the axial load-deflection characteristics of an open spring isolator with uni-directional restraint under two lateral preload conditions. (0 and 4.4 kN). Two conclusions are drawn from these results: 1. UndeT low lateral preloads, the axial stiffness characteristics of the vibration 88 isolator are trilinear. The isolator cannot carry any tensile load; therefore, the tensile stiffness is zero. Under low compressive load, a linear axial stiffness of 394 kN/m (2251 lb/in) was measured. The commercially rated axial stiffness of these springs is 350 kN/m (2000 lb/in). Under high compressive loads, the spacings between the spring coils are closed; thus, the isolator axial stiffness increases abruptly to 3500 kN/m (20000 lb/in). 2. As explained in Chapter five, under a high lateral load, when the horizontal snubber unit is engaged, the effect of Coulomb damping can be observed to act between the spring unit and the snubber unit. The coefficient of friction u which relates the extra axial frictional load to the lateral normal load was measured as 0.08. The lateral load-deflection characteristics of the open spring isolators with uni-directional restraint under various axial preloads are shown in Fig. 6.4. The results indicate the following: 1. The lateral stiffness characteristics of the isolator are bilinear. The initial lateral suffness is 98 kN/m (558 lb/in) and the second lateral stiffness is 839 kN/m (4784 lb/in). The abrupt increase in suffness occurs when the gap between the snubber unit and the spring unit is closed. The measured gap is 50 mm (.2 in). 2. The influence of axial preloads on the lateral stiffness characteristics is not significant. Based on the results from the static testing program, the inertia and the suffness characteristics to be used in the analytical models of these vibration isolated equipment systems can now be defined. For the system containing an open spring isolator with uni-directional restraint, the measured piecewise linear axial and lateral stiffness characteristics noted above are incorporated into the analytical model, which also involves friction springs with a friction coefficient of 0.08. In the analytical model of the elastomeric isolator, a linear axial stiffness of 2700 kN/m (7700 lb/in) and a 89 linear lateral stiffness of 3500 kN/m (9990 lb/in) was used. The use of linear stiffnesses to represent the nonlinear suffness characteristics of the elastomeric isolators is strictly incorrect; however, the dynamic tesung results reveal that this analytical model provides reasonable response predictions for various base motion inputs. Also, the lateral stiffness of the elastomeric isolator adopted in the analytical model seems to be high. The choice is justified however by the dynamic model identification test results presented in the next section. 6.3 M O D E L INDENTIFICATION TEST RESULTS (DYNAMIC) The dynamic model identification tests yielded decaying sinusoidal curves and frequency response curves for the two vibration isolated equipment systems. Fig. 6.5 - Fig. 6.9 show the dynamic model identification test results of the system incorporating an open spring isolator with uni-directional restraint By exciting the system at each of its natural frequencies (1.8 Hz and 4.9 Hz) until a steady state response was reached, and then suddenly stopping the table motion, the decaying sinusoidal response of the centre of gravity of the system at each natural frequency was obtained as shown in Fig. 6.5 and Fig. 6.6. Based on the logarithmatic decrement method, the modal damping ratios were found to be 0.103 and 0.105 for these two tests. These values seemed to be quite high; however, a close examination of the vibration isolators revealed thai the springs were not uniform in height, which caused the model to lean slightly toward one side. In the experiment when the system was excited, this misalignment caused the vibration springs to rub against their casings and created extra damping load. Based on the experimental results, 10% damping in all modes was used in the analytical model. Fig. 6.7 and Fig. 6.8 show the horizontal and rotational acceleration frequency response functions of this system. These responses were measured at the centre of gravity of the system. The results were obtained from experimental sinusoidal tests, 90 experimental white noise tests, analytical sinusoidal tests, and analytical white noise tests. The experimental sinusoidal and white noise tests have been described in Chapter four. The analytical sinusoidal and white noise tests were similar to their experimental counterparts; however, these results were generated analytically from the computer solution of the mathematical model for the system. In general, the results from these four types of tests show good correlation in the natural frequencies and response amplitudes. However, the experimental white noise test results do not clearly show the response peaks in the frequency range of 4.3 Hz - 4.9 Hz. The nonlinearity in the system caused the spectral analysis approach to perform poorly since the technique of spectral analysis is strictly valid for linear systems only. Based on these results, the natural frequencies of this system are 1.8 Hz, 4.3 Hz and 4.9 Hz. The predicted and measured mode shapes are shown in Fig. 6.9. The second mode predicted by the analytical approach is a pure vertical translational mode, which cannot be duplicated experimentally by exciting the system with a lateral base motion. As predicted by the analytical model, the first and the third modes of vibration are out of phase by 180 degrees. Fig. 6.10 - Fig. 6.14 show the dynamic model identification test results for the elastomeric isolator system. By exciting the system at each of its natural frequencies (4.6 Hz and 15.0 Hz) until a steady state response was reached, and then suddenly stopping the table motion, the decaying sinusoidal response of the centre of gravity of the elastomeric isolator system at each of the natural frequencies was obtained, as shown in Fig. 6.10 and Fig. 6.11. The modal damping ratios as calculated' by the logarithmatic decrement method were 0.17 for the first mode and 0.18 for the third mode. These damping values are fairly high, but the static test results have demonstrated that high hysteretic damping can be expected from an elastomeric vibration isolator system. These values were adopted as the modal damping ratios in the analytical model. 91 Fig. 6.12 and Fig. 6.13 show the horizontal and rotational acceleration frequency-response functions of the elastomeric isolator. Again these responses were monitored at the centre of gravity of the system. Based on the static test results, which yielded the nonlinear stiffness characteristics of the elastomeric isolators, it was decided not to test this system with the white noise since the spectral analysis on the response is valid for linear systems only. The experimental and analytical results show good agreement between the measured and predicted natural frequencies and response amplitudes. The natural frequencies of the system were 4.6 Hz, 8 Hz and 15.0 Hz. The three mode shapes of this system are shown in Fig. 6.14. The predicted mode shapes agree well with the corresponding measured results. 6.4 RANDOM MOTION TEST RESULTS The response of the analytical models of the two vibration isolated equipment systems were compared with the results obtained experimentally under random motion conditions. The experimental models were excited by four different random base motions: 1. El Centro 1940 N-S (maximum acceleration scaled to 0.1 g), 2. San Fernando 1971 N21E (maximum acceleration scaled to 0.18 g), 3. Taft 1950 (maximum acceleration scaled to 0.11 g), 4. Band limited (0 - 25 Hz) white noise record (maximum acceleration scaled to 0.11 g). The peak base motion acceleration was kept within 0.1 g and 0.18 g so that strong model response could be generated without damaging the experimental model. These same base motion records were used as the input excitations for the analytical models of the two vibration isolated equipment systems, thereby allowing a comparision to be made between the numerical solutions and the measured responses. 92 Fig. 6.15 - Fig. 6.22 show the predicted and measured horizontal and rotational acceleration responses at the centre of gravity of the open spring system with uni-direcuonal restraint For this system, good agreement is observed between the responses of the analytical and the experimental models. The results also show that, in some cases, the analytical model results overestimate the actual measured response values; this may be the result of assuming damping values for the analytical model which were, in fact slightly lower than those existing in the actual model. Fig. 6.23 - Fig. 6.30 show the predicted and measured horizontal and rotational acceleration responses at the centre of gravity of the elastomeric isolator equipment system. Good agreement is again observed between the predicted and measured responses. The analytical model predictions are slightly more conservative than the experimental results. The frequency content of the analytical response is higher than that of the measured response, which indicates that this analytical model is stiffer than the experimental model. This may be > due to modelling the nonlinear equipment by a linear analytical model, which is strictly incorrect The nonlinearity of the actual isolator may have softened the system when subjected to dynamic excitation. The final stage of the experimental program involved the determination of the failure modes of the two vibration isolated equipment systems. The open spring system was excited with an actual floor motion record: San Fernando 2500 Wilshire Blvd roof response N61W (1971). The base motion acceleration was increased until the vibration isolator failed. Fig. 6.31 and Fig. 6.32 show the failed vibration isolators; the open springs were about to dislodge from their casing. The recorded peak base acceleration was 1.2 g; however, the exact base motion acceleration at failure could not be determined since the isolators failed at some acceleration level during a spike in the input base motion. The elastomeric isolator equipment system was excited by the same input base motion and it survived with no apparent damage; however, the table did not experience an acceleration spike in this case. In order to determine the failure 93 mode of the elastomeric isolators, the system was excited with a sinusoidal input at the fundamental natural frequency of the system (4.6 Hz). The base motion acceleration amplitude was then increased to 2.5 g. The elastomeric isolators remained intact and suffered only minor damage: there was some tensile tearing along the base of the vibration isolators. These results indicate that the elastomeric isolator, with its high damping characteristic and nonlinearity, survived the extreme simulated seismic environment Axial load (lb) Lateral deflection (in) Lateral load-deflection curves of an elastomeric vibration isolator. T I M E ( S E C I Ct Figure 6.5 Decaying sinusoidal response time history at 1 natural frequency (1.8 Hz). oo a T I M E ( S E C ) Figure 6.6 Decaying sinusoidal response time history at 3 r d natural frequency (4.9 Hz). "i 1 r 4.0 5 0 6.0 FREQUENCY IHZI Figure 6.7 Horizontal acceleration frequency response function. p p 00. -e- Experimental sinusoidal results Analytical sinusoidal results Experimental White Noise results Analytical White Noise results Open spring isolator with uni-direcuonal restraint. ~ ~ i — 8.0 — i — 9.0 1 1— 4.0 5 0 6.0 FREQUENCY (HZ) " i — 7 0 10.0 Figure 6.8 Rotational acceleration frequency response function. Open spring isolator with uni-directional restraint. Figure 6.9 Predicted and measured mode shapes. •© Experimental sinusoidal results - — Analytical sinusoidal results FREQUENCY (HZ) Figure 6.13 Rotational acceleration frequency response function First mode Elastomeric isolator. Predicted Measured 4 Second mode Predicted and measured mode shapes. Third mode o Experimental response - Analytical response ( , Open spring isolator with uni-directional restraint. Base motion: El Centro 1940 N-S. i j 3.0 4 0 5 0 6.0 7.0 8.0 9.0 10.0 T I M E ( S E C ) Predicted and measured horizontal acceleration response time history. o o in m Q (X z o o <M_ Experimental response Analytical response Open spring isolator with uni-directional restraint. Base motion: El Centro 1940 N-S. Peak acceleration scaled to O.lg. i— 9.0 0.0 "I— 1.0 - 1 — 2.0 Figure 6.16 ~ i — 3.0 ~i 1 4 0 5 0 T I M E ( S E C ) i— 6.0 7.0 8 0 Predicted and measured rotational acceleration response time history 10. Experimental response Analytical response Open spring isolator with uni-directional restraint Base motion: San Fernando 1971 N21E Peak acceleration scaled to 0.18g. - 1 — i.o i— 2.0 3.0 4.0 SO 6.0 7.0 8.0 9.0 10. T I M E ( S E C ) Figure 6.17 Predicted and measured horizontal acceleration response time history. o.o — i — ~I— 7.0 I— 8.0 Experimental response Analytical response Open spring isolator with uni-directional restraint. Base motion: San Fernando 1971 N21E Peak acceleration scaled to 0.18g. i— 2.0 I 1 1 1— 3.0 4.0 5.0 6.0 T I M E ( S E C ) I— 8.0 0.0 Figure t o 6.18 —|— 7.0 9.0 10.0 Predicted and measured rotational acceleration response time history Experimental response m ° ~ - Analytical response ?- Base motion: Taft 1950. Peak acceleration scaled to O.llg. INI ? H 1 1 1 ! , 1 , , , 0.0 1.5 3.0 4 5 6.0 7.5 9.0 10.S 12.0 13.S T I M E ( S E C ) Figure 6.19 Predicted and measured horizontal acceleration response lime history 15 o » V C o' 1 tnrM a ? -P Experimental response Analytical response Open spring isolator with uni-directional restraint. Base motion: Taft 1950. Peak acceleration scaled lo 0.1 lg. I i i 1 r— 6.0 7.5 9.0 T I M E ( S E C I 0.0 1 1.5 3.0 ~I " 5  10.5 12.0 13.5 Predicted and measured rotational acceleration response ume history. Figure 6.20 15.0 CM Experimental response Analytical response Open spring isolator with uni-directional restraint. Base motion: Band limited (0 - 25 Hz) white noise signal. Peak acceleration scaled to O.lg. -i 1 1 1 1 7 ~ T 10.0 12.5 15.0 17.5 20.0 22 5 T I M E I S E C I 0.0 I— 2.5 - 1 — s.o 7.5 25 Figure 6.21 Predicted and measured horizontal acceleration response time history. o Experimental response Analytical response Open spring isolator with uni-directional restraint. Base motion: Band limited (0 - 25 Hi) white noise signal. Peak acceleration scaled to O.lg. 22.S 1 1 1 1 1 1 0.0 2.S 5.0 7.5 10.0 12.5 15.0 T I M E C S E C ) 17.5 20.0 25.0 Figure 6.22 Predicted and measured rotational acceleration response lime history. Experimental response Analytical response Elastomeric isolator. Base motion: El Centro 1940 N-S Peak acceleration scaled to O.lg. ~ ~ i — 6 3 o.o ~ ~ i — 0.7 1.4 i— 2.1 I 1 1— 2.0 3.5 4.2 TIME ISEC) ~I— 4.9 —1— 5.6 Figure 6.23 Predicted and measured horizontal acceleration response time history. 7.0 Experimental response Analytical response Elastomeric isolator. Base motion: El Centro 1940 N-S. Peak acceleration scaled to O.lg. o.o i — 0.7 ~T— 1.4 -1 1 1 1 1 I 2.1 2.8 3.5 4.2 4.9 5.6 T I M E I S E C ) Figure 6.24 Predicted and measured rotational acceleration response time history. Experimental response Analytical response Elastomeric isolator. Base motion: San Fernando 1971 N21E Peak acceleration scaled to 0.18g. o.o —j— 0.7 - 1 — 4.9 i 6.3 1.4 2.1 2.a 3 5 4.2 4.9 5.6 T I M E I S E C ) Figure 6.25 Predicted and measured horizontal acceleration response lime history. 7.0 Experimental response Analytical response Elastomeric isolator. Base motion: San Fernando 1971 N21E. Peak acceleration scaled lo 0.18g. - l — 6 3 — i — 0.0 0.7 1.4 2.1 Figure 6.26 Predicted and measured rotational acceleration response lime history ~ i — .4 i 1 1— 2.a 3.5 4.2 T I M E I S E C I - 1 — 4.9 I— 5.6 7.0 Experimental response Analytical response Flasiomcric isolator. Base motion: Taft 1950. Peak acceleration scaled to 0.1 lg. — i — I S 3.0 4.5 6.0 7.5 9.0 10.S 12.0 13.5 T I M E ( S E C ) Figure 6.27 Predicted and measured horizontal acceleration response lime history. ~ i — 12.0 o.o 15. Experimental response 1 1 1 1 1 1 1 i i i i 0.0 I.S 3.0 4.5 6.0 7.5 9.0 10.5 12.0 13.3 ISO T I M E ( S E C I Figure 6.28 Predicted and measured rotational acceleration response lime history. Experimental response Analytical response Base motion: Band limited (0 - 23 Hz) white noise1 signal. Peak acceleration scaled to O.lg. o.o i 0.5 I I.S I 1 1— 2.0 2.5 3.0 T I M E ( S E C ) 7 — l — 3.5 I— 4.0 1.0 1 5 2.0 2.5 3.0 3.5 4.0 4.5 Figure 6.29 Predicted and measured horizontal acceleration response lime history. 5.0 « o" s i-Experimental response Analytical response Elastomeric isolator. Base motion: Band limited (0 - 25 H i ) while noise Peak acceleration scaled to 0.1g> signal. I o.o - 1 — o.s -T— t.o - i 1 1 1 I 5 2.0 2.5 3.0 TIME (SEC) I 3.5 I— 4.0 - 1 — 4.5 Figure 6.30 Predicted and measured rotational acceleration response time history s.o to 124 Figure 6.32 Failed open spring isolator with uni-directional restraint Chapter 7 CONCLUSIONS 7.1 SUMMARY AND CONCLUSIONS Analytical and experimental studies on the behaviour of equipment isolators under seismic conditions have been investigated in three stages as follows: 1. A parametric study was conducted to investigate the effect of equipment-structure interaction on the maximum equipment response of general equipment-structure systems. The steady state response of a general two degree of freedom equipment-structure system subjected to harmonic base excitation was considered first Then the study was extended to examine the transient response characteristics of a general four degree of freedom equipment-structure system subjected to seismic excitation. From the results of this study, one can ascertain the validity of the conventional practice of non-interactive dynamic analysis and testing of equipment The effect of equipment-structure interaction depends on the parameters of the system of interest and the characteristic of the base excitation. In heavily damped equipment-structure systems, the non-interactive approach yields satisfactory ultimate equipment response predictions. The non-interactive approach is also . adequate for ultimate equipment response predictions when the equipment is very light in comparision with its supporting floor in the housing structure. For harmonic base motion inputs, the non-interactive results overestimate the equipment response when the driving frequency coincides with the tuning frequency of the subsystems and underestimate the equipment response when the driving frequency is close to the tuning frequency of the subsystems. The error in a non-interactive analysis is greatest when the driving frequency coincides with the tuning frequency of the subsystems. In the study with seismic base motion 125 126 inputs, the non-interactive results are generally conservative. 2. An experimental study of the seismic response characteristics of two different vibration isolated equipment units was conducted. A model of a prototype air handling unit mounted on vibration isolators was constructed for testing purposes. The vibration isolators investigated in this study were of the elastomeric type and the open spring type with uni-directional restraint. The two isolator models were tested both statically and dynamically. The stiffness and inertia properties of the systems were obtained in the static testing program. The dynamic testing program yielded the damping and dynamic response characteristics of the vibration isolated systems. Finally, the failure modes of the vibration isolated equipment systems were obtained by subjecting the systems to simulated seismic conditions. The static and dynamic tests provided information on the nonlinear axial and lateral stiffnesses and on the high damping characteristics of these vibration isolators. The elastomeric isolator, with its nonlinearity and high damping characteristics, survived the extreme simulated seismic conditions. The open spring isolators with uni-directional restraint failed in a brittle manner under a substantially lower level of base excitation. 3. Based on the system identification test results, analytical models of the two vibration isolated equipment systems were formulated. A computer program, utilizing time step analysis and incorporating the nonlinearity of the systems, was developed to solve the equations of motion of the systems. The analytical results, such as the frequency response curves, mode shapes, and seismic responses, were compared with the experimental results to verify the analytical model. The predicted response using the analytical model agreed well with the response values obtained experimentally. It can be concluded that the analytical model is adequate for the seismic analysis of general vibration isolated equipment systems. 127 7.2 FUTURE STUDIES The work presented in this thesis has shown that the current design rules for vibration isolated equipment systems may be inadequate and therefore there is a need to develop a rational seismic design method for vibration isolated equipment systems. In developing such a procedure the following items should be addressed: 1. The effect of equipment-structure interaction on the ultimate equipment response. 2. The nonlinearities and the damping characteristics of the vibration isolators. Studies should also be undertaken to examine the seismic response of vibration isolated equipment systems subjected to oblique base excitation. And shake table tests which can include the effect of equipment-structure interaction should be implemented in future experimental programs. Since there exists a large variety of commericial vibration isolators, it is not feasible to qualify each type of vibration isolator; however, vibration isolators which are capable of resisting expected seismic forces can be developed through the experimental and analytical approach. Finally, the effect on the seismic response characteristics of vibration isolated equipment without an inertia base needs futher investigation. Such studies are important in order to adequately protect the vibration isolated equipment systems from seismic hazards. BIBLIOGRAPHY 1) California Division of Mines and Geology Bulletin 198, 1973. 2) J. M. Biggs and J. M. Roesset, "Seismic Analysis of Equipment Mounted on a Massive Structure," in Seismic Design of Nuclear Power Plants, Ed. R. J. Hansen, pp 319-343, M.I.T. Press, 1970. 3) M. Amin, W. J. Hall, N. M. Newmark and R. P. Kasawara, "Earthquake Response of Multiple Connected Light Secondary System by Spectrum Methods," Proc. AS ME First National Congress on Pressue Vessel and Piping Technology, San Francisco, California, May 1971. 4) K. K. Kapur and L. C. Shao, "Generation of Seismic Floor Response Spectra for Equipment Design," Proc. Specialty Conference on Structural Design of Nuclear Power Plant Facilities, ASCE, Chicago, Illinois, December 1973. 5) M. P. Singh, "Generation of Seismic Floor Response Spectra," Journal of the Engineering Mechanics Division, ASCE, Vol. 101, No. EM5, PP 594-607, October 1975. 6) M. K. Chakravard and E. H. Vanmarke, "Probabilistic Seismic Analysis of Light Equipment Withih Building," Proc. Fifth World Conference on Earthquake Engineering, Vol. II, Rome, Italy, 1973. 7) A. K. Singh, "A Stochastic Model for Predicting Maximum Ressponse of Light Secondary Systems," thesis presented to the University of Illinois at Urbana, Illinois, in 1972, in partial fullfillment of the requirements for the degree of Doctor of Philosophy. 8) E. H. Vanmarke, "A Simplified Procedure for Predicting Amplitude Response Spectra and Equipmnent Response," Proc. Sixth World Conference on Earthquake Engineering, Vol. Ill, New Delhi, India, 1977. 9) R. H. Scanlan and K. Sachs, "Earthquake Time Histories and Response Spectra," Journal of the Engineering Mechanics Division, ASCE, Vol. 100, No. EM4, PP 635-655, August 1975. 10) S. H. Crandall and W. D. Mark, Random Vibration of Mechanical Systems, Academic Press, New York, N.Y., 1963. 11) A. Der Kiureghian and A. Asfura, "A New Floor Response Spectrum Method for Seismic Analysis of Multiply Supported Secondary Systems," Report No. UCB/EERC-84-04, Earthquake Engineering Research Center, University of California, Berkeley, California, June 1984. 12) J. Penzien and A. K. Chopra, "Earthquake Response of Appendage on a Multi-story Building," Proc. Third World Conference on Earthquake Engineering, New Zealand, 1965. 13) N. M. Newmark, "Earthquake Response Analysis of Reactor Structures," Nuclear Engineering and Design, Vol. 20, No. 2, PP 303-322, 1972. 128 129 14) T. Nakahata, N. M . Newmark, and W. J, Hall, "Approximate Dynamic Response of Light Secondary Systems," Structural Research Series Report No. 396, Civil Engineering Studies, Unversity of Illinois, Urbana, Illinois, 1973. 15) J. L. Sackman and J. M . Kelly, "Rational Design Methods for Light Equipment in Structures Subjected to Ground Motion," Report No. UCB/EERC-78-19, Earthquake Engineering Research Center, University of California, Berkeley, California, August 1978. 16) A. R. Robinson and C. G. Ruzicka, "Dynamic Response of Tuned Secondary Systems," Report Number UILU-Eng-80-2020, Department of Civil Engineering, Unversity of Illinois, Urbana, Illinois, November 1980. 17) R. Villaverde, N. M. Newmark, "Seismic Response of Light Attachments to Buildings," Structural Research Series Report No. 469, Civil Engineering Studies, Unversity of Illinois, Urbana, Illinois, February 1980. 18) T. Igusa and A. Der . Kiureghian, "Dynamic Analysis of Multiply Tuned and Arbitary Supported Secondary Systems," Report No. UCB/EERC-83-07, Earthquake Engineering Research Center, University of California, Berkeley, California, June 1978. 19) A. D. Kiureghian, J. L. Sackman, and B. Nour Omid, "Dynamic analysis of Light Equipment in Structures: Reponse to Stochastic Input," Journal of the Engineering Mechanics Division, ASCE, Vol. 109, No. EMI, Feburary 1975. 20) W. D. Iwan, "Predicting the Earthquake Response of Resiliently Mounted Equipment with Motion Limiting Constraints," Proc. Sixth World Conference on Earthquake Engineering, New Delhi, India, 1977. 21) W. D. Iwan, "The Earthquake Design and Analysis of Equipment Isolaton Systems," Journal of Earthquake Engineering and Structural Dynamics, ASCE, Vol. 6, PP 523-534, 1978. 22) IEEE Standard 344-1975, "IEEE Recomanded Practice for Seismic Qualification of Class IE Equipment for Nuclear Power Generating Stations," New York, 1975. 23) C. W. de Silva, Dynamic Testing and Seismic Qualification Practice, Lexington Books, D. C. Heath and Company, Lexington, Massachusetts, Toronto, 1983. 24) G. L. McGavin, Earthquake Protection of Essential Building Equipment, Wiley-Interscience, New York, 1981. 25) J. C. Wilson, W. K. Tso, and A. C. Heidebrecht, "Seismic Qualification by Shake Table Testing," Proc. Third Canadian Conference on Earthquake Engineering, Montreal, June 1979. 26) J. H. Rainer, "Dynamic Testing of Civil Engineering Structures," Proc. Third Canadian Conference on Earthquake Engineering, PP 551-574, Montreal, June 1979. 27) N. N. Nielson, "Dynamic Response of Multistory Buildings," thesis presented to the California Institute of Technology, Pasadena, California, June 1964, in partial fullfillment of the requirements for the degree of Doctor of Philosophy. 130 28) J. S. Bendat and A. G. Piersol, Random Data: Analysis and Measurement Procedures, Wiley- Interscience, New York, 1971. 29) J. S. Bendat and A. G. Piersol, Engineering Application of Correlations and Spectral Analysis, Wiley-Interscience, New York, 1980. 30) D. E. Newland An Introduction to Random Vibrations and Spectral Analysis, Longman, London, 1983. 31) G. M. Sabnis, H. G . Harris, R. N. White, and M. S. Mirzo, Structural Modelling and Experimental Techniques, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1983. 32) W. C. Hurty and M . F. Rubinstein, Dynamics of Structures, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1964. 33) R. W. Clough and J. Penzien, Dynamics of Structures, McGraw Hill, Inc., New York, 1975. 34) D. Gasparini and E H. Vanmarcke, "Evaluation of Seismic Safety of Buildings Report Number 2. Simulated Earthquake Motions Compatible with Prescribed Response Spectra," Research Report No. R76-4, Department of Civil Engineering, M.I.T., Cambridge, Massachusetts, January 1976. APPENDIX A Generation of Band Limited White Noise Record The theory of generation of band-limited white noise record presented in the following discussion is quoted without reference from Gasparini and Vanmarcke [34]. A band-limited white noise record is defined as a process having a uniform power spectral density over a specified frequency range, cjj to CJ2, as shown in Fig. A . l . By this definition, it is clear that such a process contains frequency components of equal intensity over the specified frequency range (based on squared amplitudes as a measure of intensity or power). Random records with specified power spectral density functions, such as a band-limited white noise record, can be generated by the method of summation of sinusoids. Any periodic function can be expanded into a series of N sinusoidal waves: N x(t)=Z A.sin(u> .t+0.) (A.l) i = l 1 1 1 where A j= the amplitude of the i ^ sinusoid. cjj= the frequency of the i ^ sinusoid. 0 j= the phase angle of the i ^ 1 sinusoid. Based on the assumption of squared amplitudes as a measure of intensity or N power, the steady state motion x(t) contains a total power of Z (A2./2). When the i = l 1 chosen number of sinusoids is large, the total power of x(t) is equivalent to the area under the power spectral density function G(co) shown in Fig. A.2. A sinusoid with a 131 132 frequency OK and a frequency interval ACJ =CJ . + ^- OK , contains A3./2 as power. This is approximately equivalent to the area under the power spectral density function at CJ . with a frequency interval of Aw so that (AJ./2) = G(o) .)Aw. Knowing a target power spectral density function G(a>), an array of sinusoid amplitudes A^ can be chosen to simulate a record which contains the target frequency content By keeping the array of amplitudes A . fixed, different motion records which aree similar in frequency content but different in form can be generated from different arrays of phase angles </» j. A random band-limited white noise record can be generated by randomly choosing an array of phase angles <p.. The present study uses a computer random number generator routine to obtain a random array of phase angles <j> ^  which has a uniform distribution between 0 and 27r. The band-limited white noise record of interest has a constant power spectral density between the interval 0 to 25 Hz. Five hundred sinusoids with a frequency-interval of 0.5 Hz were used to obtain the target band-limited white noise record. Figure A.1 Power spectral density funciion of a band-limited white noise signal 1 Figure A.2 Power spectral density function. APPENDIX B Spectral Analysis The theory of spectral analysis is well documented and detailed discussions on the subject can be found in Bendat and Piersol [28, 29] and Newland [30]. The following discussion is quoted without reference from the above mentioned authors. In working with the response characteristics of an ideal physical system, it is necessary to develop an approximate input/output relationship for the system. An ideal physical system process exhibits the following characteristics: 1. It is physically realizable which means that the system cannot respond to an input until the input has been applied to the system. 2. It is stable which means that a bounded input produces a bounded response. 3. It is linear which implies that the principle of superposition is valid such that the system is additive and homogeneous. 4. It has constant parameters which means that the physical characteristics of the system are time invariant The response of an ideal physical system which contains m degrees of freedom subjected to a single input x(t) can be described mathematically by a linear differential equation of the following form: a .(d ny./dt n) + .„. + a n .y. = b ( d r x / d t r ) + „ . . + b n x i = l m (B.l) ni l Or l r 0 where the coefficients a and b are the constants for a time invariant system and denotes the response of the i 1 * 1 degree of freedom of the system. 134 135 For a given excitation x(t) and given initial conditions, equation (B.l) can be solved completely to yield the response y. if the coefficients a and b were known. However, simple experimental techniques for determining these coefficients are not available. Therefore, it is more convenient to represent the relationship between the input and output signals by an alternate frequency response approach. For any arbitrary input x(t), the system output >'j(t) is given, in the time domain, by the superposition or convolution intergal yj(t)=J h.(T)x(t-r) dr i = l,....m (B.2) where h.(r) = y.(r) when x(t)=6(t) = unit impulse excitation. x(t): Figure B.1 Single input / multiple output system. Taking a finite Fourier transform of equation (B.2) over a long record length T yields the following results in the frequency domain: Y.(f ,T)=H . (0X(f ,T) i = l,...,m (B.3) 136 where Y i ( f - T ) = J y i ( i - ) e _ j ( 2 7 r f T ) d T X(f,T)=J x ( r ) e - j ( 2 , r f T ) dr H i (0=J h . ( r ) e - j ( 2 7 r f T ) dr The transformed functions in the frequency domain are complex such that each function contains a real part and an imaginary part It follows that Y*(f,T) = H*(f)X*(f.T) i i |Y.(f,T)| J = |H i(f)|'|X(f.T)|' (B.4) X *(f,T)Y .(f,T)=H j(f)|X(f,T)|2 where * denotes the complex conjugate. In terms of the one-sided spectral density functions, the cross-spectra and the auto-spectra are defined as G y v ( f ) = lim (2/T)X*(f,T)Y.(r.T) } i T - * - » (B.5) G (f) =lim (2/T)|X(f,T)|2 Subsututing equation (B.4) into equation (B.5) yields the following spectral density function relations: G (f)=H(0G (f) x y i 1 x x (B.6> G (0=|H.(0|2G (f) y. y . v ' 1 i W | xxv ' 3 v 1 Knowing the recorded input and output signals, the input/output cross-spectral and auto-spectral density functions can be determined by the techniques of fast Fourier 137 transform. And the system transfer function H.(f) can be obtained as H.(f) = G x v ( O / G m i xy. xx ( B 7 ) |H.(0| 2 = G (0/G (0 I ,v >\ V . V . XXV ' " 1" 1 The theory of discrete Fourier transform is well documented in the literature and it will not be presented further in this study. However, some important results from the sampling theorems and smoothing techniques used in this study will be presented. The most important aspect of the sampling theorems lies in the choice of the sampling rate. Suppose a band-limited signal extending from t = - ° ° to t = ° ° is sampled at a constant rate, one sample per interval A t At is chosen such that at least two samples per cycle occur in the highest frequency component of the band-limited function. Defining the Nyquist frequency as the frequency which corresponds to l/2At, if the band-limited signal is sampled at a rate lower than twice the Nyquist, it will cause aliasing problems. Aliasing means that higher frequencies present in the signal fold back toward the lower frequencies, therefore distorting the true frequency content of the signal. If the information at the Nyquist frequency is important, one has to sample at an even higher rate because the sine function in the mathematical expressions becomes identically zero at the Nyquist frequency. Therefore, information cannot be transmitted unless a higher sampling rate is chosen. Smoothing techniques are needed to reduce the errors which occur in the use of the discrete Fourier transform of the recorded signal. Ideally, the recorded signal is a continuous infinite time series; however, realistically one can only obtain a discrete finite time series in an experiment. One can regard the recorded signal as the product of a continuous infinite time series and a rectangular function of the form 0,0,0 1,1,1 1,1,1,...,0,0,0. Performing a discrete Fourier transform on such a signal 138 will yield erroneous results. The real spectrum may be distorted because the Dirichlet kernel function does not focus properly and the secondary maxima of the kernal function contributes improperly to the estimation of the Fourier coefficient in the resulting spectrum, causing leakage problems. To reduce such errors, smoothing techniques are introduced. A cosine taper window is used in place of a rectangular window to increase the focusing power of the main lobe of the kernel and to reduce leakage. APPENDIX C Program listing 139 140 L i s t i n g of EVIS.S 1 o cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 4. 3 c M A I N P R O G R A M E V I S . S C 4 c C 5 a cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 7 \* c THIS IS A TIME STEP PROGRAM WHICH CALCULATES THE RESPONSE OF 8 c A NH STORY STRUCTURE AND/OR A 3 DOF EQUIPMENT SYSTEM UNDER A 9 c BASE EXCITATION. IT ASSUMES EITHER A LINEAR VARIATION IN 10 c ACCELERATION OR AN AVERAGE ACCELERATION BETWEEN A TIME STEP. 11 c AND IT USES AN ITERATION AND EOUILIBRUIM CORRECTION PROCEDURE 12 c TO REDUCE THE ERROR IN RESPONSE INTRODUCED BY THE ABRUPT 13 c CHANGE IN SPRING STIFFNESS. THE MAXIMUM NUMBER OF STORIES IS 14 c 10. THE TOTAL DOF OF THE SYSTEM IS NB=NH+2. THE TIME STEP. 15 c DETAT, IS CHOSEN SUCH THAT XXX=BTI/DETAT IS A WHOLE NUMBER. 16 c BTI IS THE TIME INCREMENT OF THE BASE MOTION. (BTI=0.02 SEC) 17 c INPUT BASE ACCELERATION IN MM/SEC/SEC. 1 o 19 c DEFINE ALL VARIABLES 20 c 21 c NB THE TOTAL NUMBER OF DOF OF THE SYSTEM (MAX NB=13) 22 c NC NB- 1 23 c NBB NB-2 24 c AF ACCEL. RECORD OF THE GROUND AT 0 .02 SEC INCREMENT 25 c MM/SEC/SEC 26 c AFF ACCEL. RECORD OF THE GROUND M/SEC/SEC 27 c SK VECTOR OF STIFFNESSES OF SYSTEM (NB) N/M 28 c SKH1R = HORIZONTAL STIFFNESS OF RIGHT SPRING OF EOUIPIMENT N/M 29 c SKH1L = HORIZONTAL STIFFNESS OF LEFT SPRING OF EOUIPIMENT N/M 30 c SKV1R = VERTICAL STIFFNESS OF RIGHT SPRING OF EOUIPIMENT N/M 31 c SKV1L = VERTICAL STIFFNESS OF LEFT SPRING OF EOUIPIMENT N/M 32 c SKH2 » HORIZONTAL SNUBBER STIFFNESS OF EQUIPMENT N/M 33 c SKV2T = VERTICAL SNUBBER STIFFNESS OF EQUIPMENT AGAINST LIFT N/M 34 c SKV2C = VERTICAL SNUBBER STIFFNESS OF EQUIPMENT AGAINST DROP N/M 35 c AM VECTOR OF MASSES AND MOMENT OF INERTIA OF SYSTEM (NB) 36 c KG AND KG*M*M RESPECTIVELY 37 c ZETA = VECTOR OF PERCENTAGE OF CRITICAL DAMPING (NB) 38 c A VERTICAL DISTANCE BETWEEN CENTER OF MASS AND SPRINGS M 39 c B = HORIZONTAL DISTANCE BETWEEN SPRINGS M 40 c BL HORIZONTAL DISTANCE BETWEEN C . G . AND LEFT SPRING M 41 c BR HORIZONTAL DISTANCE BETWEEN C . G . AND RIGHT SPRING M 42 c c GAP BETWEEN SPRING AND HORIZONTAL SNUBBER M 43 c D1 GAP BETWEEN SPRING AND VERTICAL SNUBBER AGAINST LIFT M 44 c D2 GAP BETWEEN SPRING AND VERTICAL SNUBBER AGAINST DROP M 45 c E PERCENTAGE ECCENTRICITY BETWEEN C . G . AND SPRINGS LOCATION 46 c DETAT = TIME STEP INCREMWNT SEC 47 c DURA = DURATION OF GROUND MOTION SEC 48 c STIFF = STIFFNESS MATRIX OF THE SYSTEM (NB.NB ) 49 c SSK STIFFNESS MATRIX OF THE STRUCTURE (NBO.NBD) 50 c ESK STIFFNESS MATRIX OF THE EQUIPMENT (3.3) 51 c DC DAMPING MATRIX OF THE SYSTEM (NB.NB) 52 c SDC DAMPING MATRIX OF THE STRUCTURE (NBD,NBD ) 53 c EDC DAMPING MATRIX OF THE EQUIPMENT (3.3) 54 c AMASS = DIAGONAL MASS MATRIX OF THE SYSTEM (NB.NB) 55 c SMAS = DIAGONAL MASS MATRIX OF THE STRUCTURE (NBD.NBD) 56 c EMAS = DIAGONAL MASS MATRIX OF THE EQUIPMENT (3.3) 57 c AMASIN= INVERSE OF THE MASS MATRIX (NB.NB) 58 c DETAF = FORCE INCREMENT VECTOR (NB) 141 5 9 C P S E U F 0 = P S E U D 0 F O R C E I N C R E M E N T V E C T O R ( N B ) 6 0 C P S E U S T = P S E U D O S T I F F N E S S M A T R I X ( N B . N B ) 6 1 C D E T A X = D I S P L A C E M E N T I N C R E M E N T V E C T O R ( N B ) 6 2 C D E T A V = V E L O C I T Y I N C R E M E N T V E C T O R ( N B ) 6 3 C D E T A A = A C C E L E R A T I O N I N C R E M E N T V E C T O R ( N B ) 6 4 C F I I N E R T I A F O R C E V E C T O R ( N B ) 6 5 C S P F H = A C T U A L H O R I Z O N T A L S P R I N G F O R C E O F E Q U I P M E N T S Y S T E M 6 6 C F D D A M P I N G F O R C E V E C T O R ( N B ) 6 7 C F A F O R C E V E C T O R ( N B ) 6 8 C F U U N B A L A N C E D F O R C E V E C T O R ( N B ) 6 9 c X R E L A T I V E D I S P L A C M E N T V E C T O R ( N B ) 7 0 c V R E L A T I V E V E L O C I T Y V E C T O R ( N B ) 7 1 c A C R E L A T I V E A C C E L A T I O N V E C T O R ( N B ) 7 2 c A B C A B S O L U T E A C C E L A T I O N V E C T O R ( N B ) 7 3 c V E R T L = V E R T I C A L D I S P L A C E M E N T O F L E F T S P R I N G ( C O M P R E S S I O N + V E ) 7 4 c V E R T R = V E R T I C A L D I S P L A C E M E N T O F R I G H T S P R I N G ( C O M P R E S S I O N + V E ) 7 5 c S P F V R = V E R T I C A L S P R I N G F O R C E O F R I G H T S P R I N G ( C O M P R E S S I O N + V E ) 7 6 c S P F V L = V E R T I C A L S P R I N G F O R C E OF L E F T S P R I N G ( C O M P R E S S I O N + V E ) 7 7 c X C H O R I Z O N T A L D I S P L A C E M E N T A T S P R I N G 7 8 c XMAX = M A X . R E L A T I V E D I S P L . V E C T O R ( N B ) 7 9 c V M A X = M A X . R E L A T I V E V E L . V E C T O R ( N B ) 8 0 c A C M A X = M A X . R E L A T I V E A C C E L . V E C T O R ( N B ) 8 1 c A B C M A X = M A X . A B S O L U T E A C C E L . V E C T O R ( N B ) 8 2 c F U M A X = M A X . U N B A L A N C E D F O R C E V E C T O R ( N B ) 8 3 c S P M A X H = M A X . H O R I Z O N T A L S P R I N G F O R C E O F E Q U I P M E N T S Y S T E M 8 4 c S P M X V L = M A X . V E R T I C A L L E F T S P R I N G F O R C E O F E Q U I P M E N T S Y S T E M 8 5 c S P M X V R = M A X . V E R T I C A L R I G H T S P R I N G F O R C E O F E Q U I P M E N T S Y S T E M 8 6 c X C M A X = M A X . H O R I Z O N T A L D I S P L A C M E N T O F E Q U I P M E N T A T S P R I N G L E V E L 8 7 c V L M A X = M A X . V E R T I C A L D I S P L A C M E N T O F L E F T S P R I N G 8 8 c V R M A X = M A X . V E R T I C A L D I S P L A C M E N T O F R I G H T S P R I N G 8 9 c D F S M A X = M A X . E R R O R I N S P R I N G F O R C E D U E T O A B R U P T C H A N G E 9 0 c I N S T I F F N E S S . 9 1 c A C F A B S O L U T E A C C E L A T I O N R E C O R D DF T H E S Y S T E M 9 2 c DA V E C T O R W H I C H C O N T A I N S T H E LOWER H A L F B A N D O F S T I F F 9 3 c D B V E C T O R W H I C H C O N T A I N S T H E LOWER H A L F B A N D O F A M A S S 9 4 c D E V E C T O R W H I C H C O N T A I N S T H E E I G E N V A L U E S OF S Y S T E M 9 5 c T A C V E C T O R O F T I M E T O M A X I M U M R E L A T I V E A C C E L . 9 6 c T A B C = V E C T O R O F T I M E T O M A X I M U M A B S O L U T E A C C E L . 9 7 c T X V E C T O R O F T I M E T O M A X I M U M R E L A T I V E D I S P L . 9 8 c T V V E C T O R O F T I M E T O M A X I M U M R E L A T I V E V E L . 9 9 c T F U V E C T O R O F T I M E T O M A X I M U M U N B A L A N C E D F O R C E . 1 O 0 c T S P T I M E T O M A X I M U M H O R I Z O N T A L S P R I N G F O R C E 1 0 1 c T S P V R = T I M E T O M A X I M U M V E R T I C A L R I G H T S P R I N G F O R C E 1 0 2 c T S P V L = T I M E T O M A X I M U M V E R T I C A L L E F T S P R I N G F O R C E 1 0 3 c T D F T I M E T O M A X . E R R O R I N S P R I N G F O R C E 1 0 4 c T X C T I M E T O M A X I M U M H O R I Z O N T A L S P R I N G D I S P L A C M E N T 1 0 5 c T V L T I M E T O M A X I M U M V E R T I C A L L E F T S P R I N G D I S P L A C M E N T 1 0 6 c T V R T I M E T O M A X I M U M V E R T I C A L R I G H T S P R I N G D I S P L A C M E N T 1 0 7 c D V M A T R I X W H I C H C O N T A I N S T H E E I G E N V E C T O R O F S Y S T E M 1 0 8 c GAMMA = V E C T O R O F M O D A L P A R T I C I P A T I O N F A C T O R S 1 0 9 c A L P H A = M U L T I P L I C A T I O N F A C T O R D E P E N D I N G O N T H E T H E T I M E S T E P 1 1 0 c M E T H O D C H O S E N . 1 11 c B E T A « M U L T I P L I C A T I O N F A C T O R D E P E N D I N G O N T H E T H E T I M E S T E P 1 1 2 c M E T H O D C H O S E N . 1 1 3 c 1 1 4 c M I S C E L F L A G S A N D C H E C K S : 1 1 9 1 1 6 c I F L . I F A , I A , I Q , J Q . K Q . I E R , I F L A G . I F L A A . I F L A A A , C H E C K , I C H E C K . 1 1 C H E C K 142 1 17 C 1 18 IMPLICIT REAL*8(A-H.0-Z) 1 19 DIMENSION STIFF(13.13),AMASS(13.13).AMASIN(13.13),X(13).V(13) 120 DIMENSION ESK(3.3).EMAS(3.3),SMAS(10.10),SDA(100),EDA(9) 121 DIMENSION SDB(55).EDB(6).SDE( 10).EDE(3).SDV( 10. 10).EDVI3.3) 122 OIMENSION DETAF(13) ,PSEUFO(13) ,PSEUST(13,13),DETAX(13) 123 DIMENSION F I ( 1 3 ) . F S ( 1 3 ) , F A ( 1 3 ) . A C ( 1 3 ) ,DETAVl 13),ACF(9G012) 124 DIMENSION DB(91),DA(169 ).DE(13),AFF(60000),AF(7600) 125 DIMENSION DC(13. 13),FD< 13),DV< 13. 13) .IQ( 18 ) .IA< 3) ,dQ< 18 ) ,KQ(18) 126 DIMENSION XMAX(13).VMAX(13),ACMAX(13),ABC(13),ABCMAX( 13) 127 DIMENSION TAC(13),TABC(13).TX( 13 ) .TV(13),DETX(13).SDC<10,10) 128 DIMENSION AM(13),SK(13),ZETA( 13).IPERM(24 ) ,T(13,13) ,EDC(3,3) 129 DIMENSION EZETA(3,3),SZETA(10.10).YACL(100O0).YACR(10000) 130 DIMENSION DETAA(13).FU(13),FUMAX(13),TFU(13) 131 INTEGER STAR(1)/'*'/ 132 C 133 C READ IN PROGRAM OPTIONS 134 C 135 WRITE(5,199) 136 READ(5,STAR)MET 137 WRITE(5,200) 138 READ(5,STAR ) NB 139 WRITE(5,201) 140 READ(5,STAR ) IFL 141 IF(IFL.NE.O) IFA = 7 142 IF(IFL.NE.O) GOTO 40 143 WRITE(5,202) 144 READ(5.STAR)IFA 145 40 WRITE(5,203) 146 READ(5,STAR ) ICHECK 147 IF(ICHECK.EO.O) GOTO 25 148 WRITE(5,204) 149 READ(5,STAR ) IICHEC 150 C 151 C CALL SUBROUTINE TO READ THE INPUT STRUCTURAL AND 152 C EQUIPMENT DATA 153 C 154 25 CALL READ(SK.ZETA,AM,AFF,DURA.N.SKH2,A,BL,BR,C,NB,DETAT,IFL.IT, 155 1T0L,IFA,SKH1L.SKH1R,SKV1L.SKV1R.SKV2T.SKV2C.D 1 ,D2,EMAS.SMAS 156 1,EZETA,SZETA,B,E,FMU) 157 c 158 c SET INITIAL CONDINTIONS 159 c 160 DO 100 1=1,NB 161 X(I ) =0.DO 162 V(I)=0.D0 163 AC(I)=0.D0 164 ABC(I )=0.D0 165 ACMAX(I)=O.DO 166 ABCMAX(I)=0.D0 167 XMAX(I)=0.D0 168 DETAX(I)=0.DO 169 DETAA(I)=0.D0 170 DETAV(I)=0.DO 171 FU(I)=0.DO 172 FUMAX(I)=0.00 173 TAC(I)=0.00 174 TABC(I)=0.DO 143 175 T X ( I ) = 0 . D 0 176 TFU(I)=O.DO 177 TV(I)=O.DO 178 VMAX(I)=O.DO 179 100 CONTINUE 180 DO 101 1=1.18 181 IO(I)=0 182 J 0 ( I ) = O 183 K0(I)=O 184 101 CONTINUE 185 I FLAA=0 186 XCMAX=0.D0 187 SPMAXH'O.DO 188 DFSMAX=O.DO 189 SPMXVL=O.DO 190 SPMXVR=O.DO 191 VRMAX=O.DO 192 VLMAX=0.D0 193 VBL=O.DO 194 VBR=0.D0 195 TSPVL=0.D0 196 TSPVR=O.DO 197 TVL=O.DO 198 TVR=0,. DO 199 TDF=o\DO 200 TSP=0.D0 201 TXC=0.D0 • 202 ICOUNT=0 203 DT1=0.D0 204 NC=NB-1 205 NBB=NB-2 206 NBD=NB-3 207 NBV=NB 208 BT=DETAT 209 DETVL=0.D0 210 DETVR=0.D0 211 I F ( I F A . E Q . I ) NBV=NB-3 212 I F ( I F A . E O . O ) NBV = 3 213 I F ( M E T . E O . O ) BETA=1.DO/6.DO 214 I F ( M E T . E O . O ) ALPHA=0.5D0 215 I F ( M E T . E Q . I ) BETA=0.25DO 216 I F ( M E T . E O . I ) ALPHA=0.5D0 217 ALBE *ALPHA/BETA 218 ALBR=1.DO-(0 .5D0*ALBE) 219 C 220 C SET TIME STEP P INTERVAL < (SMALLEST NATURAL PERIOD /10) 221 C 222 XX=DURA/DETAT+1.5DO 223 L=IDINT(XX) 224 WRITE(6 .205)L 225 WRITE(6 ,206) 226 TIME=-DETAT 227 C 228 C SET UP MASS MATRIX 229 C 230 CALL MASS(AM,AMASS.AMASIN,NBV, IFL) 231 C 232 C BEGIN TIME STEP 144 233 C 234 DO 30 I - 1 . L 235 IERR=0 236 IER=0 237 NUM=0 238 11=1+1 239 JE=1 240 IFLAG=0 241 IFLAAA=0 242 IDIREC=0 243 IXX=0 244 A S = ( A F F ( I I ) - A F F ( I ) ) / D E T A T 245 CALL D R I V E ( X . S K . A M A S S . N B . I F L . A , B L , B R , C , S K H 2 . D E T A T . AM 246 , A F F , I . N . N B B . N C . I T . V . A C , I ERR,ABC.ZETA,DETAX.TIME,NBV 247 . J E , I F A . I C O U N T . X C C , I E R , A M A S I N , A S . I F L A G . X C . D X 1 . D X 2 . T 0 L , B T 248 . D E T A F S . I F L A A , I X X , I D IREC .ALPHA.BETA .ALBE ,ALBR.SDC .EDC 249 ,VERTL .VERTR,SKH1L .SKH1R,SKV1L .SKV1R.SKV2T .SKV2C ,D1 ,D2 250 , V B L , V B R . S P F H . S P F V R . S P F V L . I A .10 . J O . K Q . I F L A A A , E M A S . S M A S 251 . E Z E T A . S Z E T A , N B D , E , B , Y A L , Y A R , F U . D E T V L , D E T V R , F M U ) 252 C 253 C PERFORM ITERATION PROCEDURE TO MINIMIZE SPRING FORCE ERROR. 254 C 255 I F ( I E R . N E . 1 ) GOTO 310 256 DT1=DETAT 257 IFLAG= 1 258 10 DT1=DT1/DX2*DX1 259 IFLAA= 1 260 CALL D R I V E ( X , S K , A M A S S , N B . I F L . A , B L , B R . C , S K H 2 . D T 1 , A M 261 , A F F , I . N . N B B . N C , I T , V . A C . I ERR.ABC,ZETA.DETAX.TIME.NBV 262 , J E , I F A , I C O U N T . X C C , I E R . A M A S I N . A S . I F L A G , X C . D X 1 . D X 2 . T 0 L . B T 263 . D E T A F S . I F L A A . I X X , I D I R E C . A L P H A . B E T A . A L B E . A L B R . S D C , E D C 264 .VERTL .VERTR.SKH1L ,SKH1R.SKV1L ,SKV1R,SKV2T ,SKV2C .D1 .D2 265 . V B L , V B R . S P F H . S P F V R . S P F V L . I A .10. J Q . K O . I F L A A A , E M A S . S M A S 266 . E Z E T A . S Z E T A , N B D , E . B . Y A L . Y A R . F U . D E T V L . D E T V R . F M U ) 267 I F ( I F L A G . E 0 . 2 ) GOTO 177 268 NUM=NUM+1 269 IF(NUM.GE .100) GOTO 12 270 GOTO 10 271 12 WRITE(6.207)TIME 272 GOTO 400 273 C 274 c COMPLETE THE TIME STEP AFTER THE SUCCESSFUL ITERATION. 275 c 276 177 DT2=DETAT-DT1 277 IFLAA= 1 278 CALL D R I V E ( X , S K , A M A S S , N B , I F L . A . B L . B R , C , S K H 2 . D E T A T . AM 279 , A F F , I . N . N B B . N C , I T , V . A C , I ERR.ABC.ZETA,DETAX.TIME,NBV 280 , J E , I F A . I C O U N T . X C C . I E R , A M A S I N , A S . I F L A G , X C , D X 1 . D X 2 . T 0 L . B T 281 . D E T A F S , I F L A A , I X X , I D I R E C . A L P H A . B E T A . A L B E . A L B R . S D C , E D C 282 ,VERTL ,VERTR.SKH1L ,SKH1R,SKV1L .SKV1R.SKV2T ,SKV2C ,D1 ,D2 283 , V B L , V B R . S P F H . S P F V R . S P F V L . I A .10, J O . K Q . I F L A A A , E M A S . S M A S 284 . E Z E T A . S Z E T A , N B D , E . B . Y A L , Y A R . F U . D E T V L . D E T V R , FMU) 285 c 286 c FIND MAX. X, V. AC, ABC, XC. SPFH SPFVR. SPFVL, VERTR, VERTL 287 c AND DETAFS AT THEIR RESPECTIVE TIME OF OCCURENCE 288 c 289 310 TIME=TIME+DETAT 290 DO 37 K=1,NBV 145 291 IF(DABS(AC(K)) .GT .DABS(ACMAX(K) ) ) TAC(K)=TIME 292 IF(DABS(ABC(K) ) .GT.DABS(ABCMAX(K))) TABC(K)=TIME 293 IF(DABS(V(K)) .GT .DABS(VMAX(K ) ) ) TV(K)=TIME 294 IF(DABS(X(K)) .GT .DABS(XMAX(K ) ) ) TX(K) = TIME 295 IF(DABS(FU(K)) .GT .DABS(FUMAX(K) ) ) TFU(K)=TIME 296 IF(DABS(X(K)) .GT .DABS(XMAX(K ) ) ) XMAX(K)=X(K) 297 IF(DABS(V(K)) .GT .DABS(VMAX(K ) ) ) VMAX(K) = V(K) 298 I F ( D A B S ( A C ( K ) ) . G T . D A B S ( A C M A X ( K ) ) ) ACMAX(K )=AC(K ) 299 I F ( D A B S ( A B C ( K ) ) . G T . D A B S ( A B C M A X ( K ) ) ) ABCMAX(K)=ABC(K) 300 IF(DABS(FU(K)) .GT .DABS(FUMAX(K ) ) ) FUMAX(K)=FU(K ) 301 37 CONTINUE 302 IF(DABSCXC).GT.DABS(XCMAX ) ) TXC = TIME 303 IF(DABS(VERTR).GT.DABS(VRMAX ) ) TVR = TIME 304 IF(DABS(VERTL) .GT .DABS(VLMAX) ) TVL = TIME 305 IF(DABS(SPFH).GT.DABS(SPMAXH ) ) TSP = TIME 306 IF(DABS(SPFVL) .GT .DABS(SPMXVL)) TSPVL=TIME 307 IF(DABS(SPFVR) .GT .DABS(SPMXVR)) TSPVR=TIME 308 IF(DABS(DETAFS) .GT .DABS(DFSMAX)) TDF=TIME 309 IF(DABS(XC) .GT.DABS(XCMAX)) XCMAX = XC 310 IF(DABS(SPFVL ) .GT.DABS(SPMXVL)) SPMXVL = SPFVL 311 IF(DABS(SPFVR).GT.DABS(SPMXVR ) ) SPMXVR = SPFVR 312 IF(DABS(SPFH).GT.DABS(SPMAXH) ) SPMAXH=SPFH 313 IF(DABS(VERTR).GT.DABS(VRMAX ) ) VRMAX = VERTR 314 IF(DABS(VERTL) .GT.DABS(VLMAX ) ) VLMAX = VERTL 315 DFSMAX=DMAX1(DFSMAX.DABS(DETAFS)) 316 C 317 C STORE ABSOLUTE ACCEL. RECORD OF EACH FLOOR IN ACF 318 C 319 ID=0 320 I F ( I C H E C K . E O . O ) GOTO 30 321 DO 78 IP=1,NBV 322 ACF(ID+I)=ABC(IP) 323 ID=ID+8001 324 78 CONTINUE 325 Y A C L ( I ) = Y A L * 1000.DO 326 YACR(I)=YAR*1000.DO 327 30 CONTINUE 328 C 329 C OUTPUT THE ACCEL. RECORD OF EACH DOF. (MM/SEC/SEC) IN UNIT 330 C AND THE MAXIMUM VALUES IN UNIT 6. 331 C 332 I F ( I C H E C K . E O . O ) GOTO 197 333 ID=0 334 DO 170 I Z M . N B V 335 I F ( I I C H E C . E O . O ) GOTO 172 336 I F ( I I C H E C . N E . I Z ) GOTO 173 337 172 IDD=ID+1 338 IDL=ID+L 339 WRITE(8 ,208)L ,DETAT 340 DO 169 I ' I D D . I D L 341 A C F ( I ) = A C F ( I ) * 1 0 0 0 . D 0 342 169 CONTINUE 343 W R I T E ( 8 , 2 0 9 ) ( A C F ( I ) . I = I D D . I D L ) 344 173 ID=ID+8001 345 170 CONTINUE 346 I F ( I I C H E C . N E . O ) GOTO 197 347 C WRITE(10 .208)L .DETAT 348 C W R I T E ( 1 0 . 2 0 9 ) ( Y A C R ( I ) . I = 1 . L ) 146 349 C WRITE(11 ,208)L ,DETAT 350 C W R I T E ( 1 1 , 2 0 9 ) ( Y A C L ( I ) , I = 1 . L ) 351 C 352 197 WRITE(6 .210) 353 WRITE(6.21 1 ) 354 WRITE(6 .212) 355 WRITE(6 .213)(XMAX(K) ,K=1,NBV ) 356 W R I T E ( 6 . 2 1 4 ) ( T X ( K ) . K = 1 . N B V ) 357 WRITE !6 .215)(VMAX(K) .K=1 .NBV) 358 W R I T E ( 6 . 2 1 4 ) ( T V ( K ) , K = 1 , N B V ) 359 WRITE(6,216 )(ACMAX(K) .K=1 .NBV) 360 W R I T E ( 6 . 2 1 4 ) ( T A C ( K ) , K = 1 . N B V ) 361 WRITE(6 ,217)(ABCMAX(K) ,K=1,NBV) 362 W R I T E ( 6 . 2 1 4 ) ( T A B C ( K ) , K = 1 , N B V ) 363 WRITE(6 .218)(FUMAX(K) ,K=1,NBV) 364 W R I T E ( 6 . 2 1 4 ) ( T F U ( K ) , K = 1 , N B V ) 365 WRITE(6.219)SPMAXH,TSP 366 WRITE(6.220)XCMAX,TXC 367 WRITE(6.221)IC0UNT 368 WRITE(6.222)DFSMAX,TDF 369 WRITE(6.223)SPMXVL,TSPVL 370 WRITE(6.224)VLMAX.TVL 371 WRITE(6,225)SPMXVR,TSPVR 372 WRITE(6,226)VRMAX,TVR 373 GOTO 400 374 26 WRITE(5 .227) 375 WRITE(6 ,227) 376 C 377 199 FORMAT('WHICH TIME INTEGRATION METHOD DO YOU WANT?', 378 1 '(0=LINEAR ACCEL . , 1=AVERAGE A C C E L . ) ' ) 379 200 FORMAT('HOW MANY DOF ARE THERE IN THE SYSTEM? (MAX NB=13) ' ) 380 201 FORMAT('TO WHICH FLOOR IS THE EOUIPMENT ATTACHED?' , 381 1 '(0=NOT ATTACHED.1=FIRST,2 = SEC0ND ) ' ) 382 202 FORMAT('WHICH RESPONSE DO YOU WANT?(0=EOUIPMENT.1^STRUCTURE)') 383 203 FORMAT('DO YOU WANT THE ACCEL. RECORD?(YES=1,N0 = 0 ) ' ) 384 204 FORMAT('WHICH DOF RESPONSE RECORD DO YOU WANT?', 385 1 ' (0=ALL . 1=FIRST.2 = SEC0ND . . . ) ' ) 386 205 FORMAT(I 10) 387 206 FORMAT(/, '******** RESULTS **•****• ' ) 388 207 FORMAT('ITERATION FAILED @ TIME = ' . F 1 5 . 7 ) 389 208 F O R M A T ( I 1 0 . F 1 5 . 6 ) 390 209 FORMAT(8F10.0) 391 210 FORMAT(//, 'UNITS OF MAXIMUM DISPLACEMENTS = M OR RAD') 392 211 FORMAT('UNITS OF MAXIMUM VELOCITIES = M/SEC OR R A D / S E C ) 393 212 FORMAT('UNITS OF MAXIMUM ACCELERATIONS = M/SEC/SEC OR ' . 394 1 'RAD/SEC/SEC') 395 213 FORMAT(//, 'MAXIMUM REL. DISPL . ' . 5 F 1 6 . 5 ) 396 214 FORMAT( ' «» RESP. TIME = ' . 5 F 1 6 . 4 ) 397 215 FORMAT(/,'MAXIMUM REL. VEL . ' . 5 F 1 6 . 5 ) 398 216 FORMAT(/,'MAXIMUM REL. A C C E L . ' , 5 F 1 6 . 5 ) 399 217 FORMAT(/, 'MAXIMUM ABS. ACCEL . ' . 5 F 1 6 . 5 ) 400 2 18 FORMAT( / , 'MAX. UNBALANCED FORCE' , 5 F 1 6 . 5 ) 401 219 FORMAT(/,'MAX HORIZONTAL SPRING FORCE = ' . F 1 5 . 5 . 402 1 ' P TIME = ' . F 1 0 . 5 ) 403 220 FORMAT(/, 'MAX HORIZONTAL SPRING DISPL = ' . F 1 5 . 5 . 404 1 ' (S TIME = ' . F 1 0 . 5 ) 405 221 FORMAT(/, 'THE H OF TIMES THAT THE EQUIP. HITS HORIZONTAL SNUBBER 406 1 ,110) 147 407 222 FORMAT(/,'MAX SPRING FORCE ERROR = ' . F 1 5 . 5 , ' © TIME = ' . F 1 0 . 5 ) 408 223 FORMAT(/,'MAX VERTICAL LEFT SPRING FORCE = ' . F 1 5 . 5 . 409 1 ' & TIME = ' . F 1 0 . 5 ) 410 224 FORMAT(/,'MAX VERTICAL LEFT SPRING DISPL = ' . F 1 5 . 5 , 411 1 * TIME = ' . F 1 0 . 5 ) 412 225 FORMAT(/,'MAX VERTICAL RIGHT SPRING FORCE = ' . F 1 5 . 5 , 413 1 9 TIME = ' . F 1 0 . 5 ) 414 226 FORMAT(/,'MAX VERTICAL RIGHT SPRING DISPL = ' . F 1 5 . 5 . 415 1 0 TIME = ' , F 1 0 . 5 ) 416 227 FORMAT('SOLUTION F A I L E D ) 417 400 STOP 4 18 END 4 19 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 420 C C 421 C S U B R O U T I N E R E A D C 422 C C 423 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 424 C 425 C THIS SUBROUTINE READS FROM UNIT 4 THE STRUCTURAL DATA FLOOR 426 C BY FLOOR AND THEN THE EOUIPMENT DATA. IT ALSO READS THE 427 C GROUND MOTION DATA FROM UNIT 7. THE GROUND ACCELERATION TIME 428 C INCREMENT IS B T I . 429 C 430 SUBROUTINE R E A D ( S K , Z E T A , A M , A F F , D U R A , N , S K H 2 , A . B L , B R , C . N B . D E T A T , 431 1 I F L , I T . T O L , I F A . S K H 1 L , S K H 1 R , S K V 1 L , S K V 1 R , S K V 2 T . S K V 2 C , 0 1 , D 2 , E M A S , 432 1SMAS ,EZETA ,SZETA ,B ,E ,FMU) 433 IMPLICIT R E A L * 8 ( A - H , 0 - Z ) 434 DIMENSION T I T L E ( 8 0 ) 435 DIMENSION A F ( 7 6 0 0 ) , A F F ( 6 0 0 0 0 ) , A M ( 13 ) . S K ( 1 3 ) , Z E T A ( 1 3 ) 436 DIMENSION EMAS ( 3 , 3 ) , SMAS ( 10, 10 ) . E ZET A ( 3 . 3 ) ,-SZET A ( 10, 10) 437 SI=1000.DO 438 NBB=NB-2 439 NBD=NB-3 440 NC=NB-1 441 DO 8 1=1.NBD 442 R E A D ( 4 , 2 0 ) S K ( I ) , Z E T A ( I ) , A M ( I ) 443 8 CONTINUE 444 READ(4 ,20)SKH1L,SKH1R,ZETA(NBB) ,AM( NBB) 445 READ(4 ,20)SKV1L ,SKV1R ,ZETA(NC) .AM(NC) 446 READ!4 ,20)SKH2 ,SKV2T ,SKV2C ,ZETA(NB) ,AM(NB) 447 R E A D ( 4 , 2 0 ) A , B , C . D 1 , D 2 . D E T A T . E 448 READ(4 ,21)T0L .FMU 449 R E A D ( 7 , 2 2 ) T I T L E 450 R E A D ( 7 , 2 3 ) N , B T I 451 R E A D ( 7 , 2 4 ) ( A F ( L + 1 ) , L = 1 , N ) 452 C 453 C CALCULATE LOCATION OF C . G . KNOWING THE ECCENTRICITY 454 C 455 BL=B/2.D0-E*B/100.D0 456 BR=B/2.D0+E*B/100.D0 457 C 458 C ECHO INPUT STRUCTURAL AND EOUIPMENT DATA 459 C 460 10=5 461 7 W R I T E ( I 0 . 2 2 ) T I T L E 462 I F ( I F L . N E . O ) W R I T E ( 1 0 , 2 5 ) I F L 463 I F ( I F A . E Q . O ) W R I T E U 0 . 2 6 ) 464 I F ( I F A . E O . I ) W R I T E U 0 . 2 7 ) 148 465 W R I T E ( 1 0 . 2 8 ) 466 DO 10 1=1,NBD 467 W R I T E ( 1 0 , 2 9 ) 1 , 5 K ( I ) , ZETA(I ) ,AMI I) 468 10 CONTINUE 469 I F ( I 0 . E 0 . 5 ) W R I T E ( 5 , 5 1 ) 470 I F ( I 0 . E Q . 5 ) R E A D ! 5 . 5 2 ) N X X 47 1 I F ( ( I 0 . E 0 . 5 ) . A N D . ( N X X . E Q . 1) ) STOP 472 W R I T E ( 1 0 , 3 0 ) 473 WRITE(I0 ,31 ) 474 W R I T E ( 1 0 . 3 2 ) A M ( N B ) . A M ( N C ) 475 W R I T E ( I 0 . 3 3 ) Z E T A ( N B B ) . Z E T A ( N O . Z E T A ( N B ) 476 W R I T E ( 1 0 . 3 4 )SKH1L 477 WRITE(I0.35)SKH1R 478 WRITE(I0 .36)SKV1L 479 WRITE(I0 ,37)SKV1R 480 WRITE(I0 .38)SKH2 481 W R I T E ( 1 0 . 3 9 ) S K V 2 T 482 WRITE(I0 .40)SKV2C 483 WRITE(10,41 )A 484 W R I T E ( 1 0 . 4 2 )BL 485 WRITE(IO,43)BR 486 W R I T E ( I 0 . 4 4 ) C 487 WRITE(I0 .45)D1 488 WRITE(I0 .46)D2 489 W R I T E ( I 0 . 4 7 ) E 490 I F ( I 0 . E Q . 5 ) W R I T E ( 5 , 5 1 ) 491 I F ( I 0 . E Q . 5 ) R E A D ( 5 . 5 2 ) N X X 492 I F ( ( 1 0 . E Q . 5 ) .AND.(NXX.EO. 1)) STOP 493 WRITE(I0 ,48)DETAT 494 W R I T E ( I 0 . 4 9 ) T 0 L 495 W R I T E ( I 0 . 5 0 I N . B T I 496 10=10+1 497 I F ( I 0 . E 0 . 6 ) G0T07 498 C 499 20 F0RMAT(9F20.4) 500 2 1 FORMAT(2F19.15) 501 22 FORMAT(80A1) 502 23 F O R M A T ( I 1 0 . F 1 5 . 6 ) 503 24 F0RMAT(8F10.0) 504 25 FORMAT(/. 'COUPLED ANALYSIS--EQUIPMENT ATTACHED TO FLOOR /»'.I3) 505 26 FORMAT(/, 'UNCOUPLED ANALYSIS--EQUIPMENT RESPONSE ONLY' I 506 27 FORMAT!/, 'EQUIPMENT NOT ATTACHED - STRUCTURAL RESPONSE ONLY') 507 28 FORMAT(/ , ' FLOOR' ,3X , ' INTERSTORY STIFFNESS ( N / M ) ' . 3 X . 508 1 'DAMPING R A T I O ' , 6 X . ' M A S S ( K G ) ' ) 509 29 FORMAT!15 ,3F20 .4) 510 30 FORMAT!//. 'EQUIPMENT DATA ' ) 51 1 31 FORMAT( '*****»•»***»*•• ) 512 32 FORMAT!'MASS !KG) = ' , F 1 1 . 4 . 4 X . 513 1 'MOMENT OF INERTIA !KG*M*M) = ' . F 1 1 . 4 ) 514 33 FORMAT!'EQUIPMENT DAMPING RATIOS = ' . 3 F 1 0 . 4 ) 515 34 FORMAT!'HORIZONTAL STIFFNESS OF LEFT SPRING (N/M) = ' , F 1 3 . 4 ) 516 35 FORMAT!'HORIZONTAL STIFFNESS OF RIGHT SPRING (N/M) = ' . F 1 3 . 4 ) 517 36 FORMAT! 'VERTICAL STIFFNESS OF LEFT SPRING (N/M) = ' . F 1 3 . 4 ) 518 37 FORMAT! 'VERTICAL STIFFNESS OF RIGHT SPRING (N/M) = ' , F 1 3 . 4 ) 519 38 FORMAT!'HORIZONTAL SNUBBER STIFFNESS (N/M) = ' . F 1 6 . 4 ) 520 39 FORMAT! 'VERTICAL SNUBBER STIFFNESS AGAINST LIFT (N/M) = ' . F 1 G 4 ) 521 40 FORMAT*'VERTICAL SNUBBER STIFFNESS AGAINST DROP !N/M) = ' , F 1 6 . 4 ) 522 4 1 FORMAT!'VERTICAL DISTANCE BETWEEN CENTRE OF MASS AND SPRINGS' 149 523 1 , ' (M) = ' , F 7 . 4 ) 524 42 F0RMAT('HORIZONTAL DISTANCE BETWEEN C . G . AND LEFT SPRINGS (M) =' 525 1 , F 7 . 4 ) 526 43 FORMAT!'HORIZONTAL DISTANCE BETWEEN C . G . AND RIGHT SPRINGS (M) =' 527 1 . F 7 . 4 ) 528 44 FORMAT('GAP BETWEEN SPRING AND HORIZONTAL SNUBBER (M) = ' , F 7 . 5 ) 529 45 FORMAT('GAP BETWEEN SPRING AND VERTICAL SNUBBER AGAINST L I F T ' 530 1 . ' (M) = ' . F 7 . 5 ) 531 46 FORMAT('GAP BETWEEN SPRING AND VERTICAL SNUBBER AGAINST DROP' 532 1 . ' (M) = ' . F 7 . 5 ) 533 47 FORMAT('PERCENTAGE ECCENTRICITY BETWEEN C . G . AND SPRINGS = ' . F 7 . 4 ) 534 48 FORMAT(//. 'TIME INCREMENT (SEC) = ' . F 7 . 4 ) 535 49 FORMAT('TOLERENCE (M) = ' , F 1 5 . 1 3 ) 536 50 FORMAT(15. ' INPUT ACCEL. PT . * ' . F 8 . 4 . ' SECS. ' ) 537 51 FORMAT(//, 'DO YOU WISH TO CONTINUE ? 0 => YES : 1 = > NO') 538 52 FORMAT(11 ) 539 C 540 C ASSUMING A LINEAR VARIATION BETWEEN TIME STEP SUBDIVIDE 541 C THE INPUT ACELERATION ACCORDINGLY. 542 C 543 AF(1)=0 .DO 544 DURA=BTI*N 545 XXX=BTI/DETAT+0.5D0 546 IST=IDINT(XXX) 547 I F ( I S T . L E . 1) GOTO 11 54B ISS=IST-1 549 M=N+1 550 DO 1 d=1.M 551 dd=J*IST- ISS 552 A F F ( d d ) = A F ( d ) / S I 553 1 CONTINUE 554 DO 2 K=1,N 555 DO 3 d=1,ISS 556 d«Jd=K*IST+1 557 dd=IST*K+d-ISS 558 dddd=K*IST-ISS 559 A F F ( J d ) = ( A F F ( J d J ) - A F F ( d d d d ) ) / I ST*J+AFF(ddJJ) 560 3 CONTINUE 561 2 CONTINUE 562 AFF(ddd+1)=O.DO 563 GOTO 12 564 1 1 00 9 1=1,N 565 A F F ( I ) = A F ( I ) / S I 566 9 CONTINUE 567 AFF(N+1)=0.D0 568 C 569 C ZERO EOUIPMENT AND STRUCTURE MASS MATRIX AND DAMPING VECTOR 570 C 571 12 CALL D G S E T ( E M A S , 3 , 3 . 3 . 0 . D 0 ) 572 CALL D G S E T ( S M A S , 1 0 . 1 0 . 1 0 , 0 . D O ) 573 CALL D G S E T ( E Z E T A , 3 . 3 , 3 . 0 . D 0 ) 574 CALL D G S E T ( S Z E T A , 1 0 , 1 0 , 1 0 , 0 . D O ) 575 C 576 C SET UP EOUIPMENT AND STRUCTURE MASS MATRIX 577 C 578 DO 13 1=1,NBD 579 S M A S ( I . I ) = A M ( I ) 580 13 CONTINUE 150 581 DO 14 1=1,3 582 EMAS(I,I)=AM(I+NBD) 583 14 CONTINUE 584 00 15 1=1.NBD 585 S Z E T A ( I , I ) = Z E T A ( I ) 586 15 CONTINUE 587 00 16 1=1,3 588 EZETA(I , I )=ZETA(I+NBD) 589 16 CONTINUE 590 I F ( I F L . E 0 . N B - 3 ) GOTO 5 591 I F ( I F L . E O . O ) GOTO 6 592 C 593 C REARRANGE THE MASSES AND DAMPING VECTORS TO THE EOUIPMENT 594 C LOCATION. 595 C 596 EMASSX=AM(NBB) 597 EMASSY=AM(NC) 598 E A<J = AM( NB) 599 EOZETA=ZETA(NBB) 600 EOZETB=ZETA(NC) 601 EOZETC=ZETA(NB) 602 IFJ=IFL+4 603 IFK=IFL+3 604 IFD=IFL+2 605 IFE=IFL+1 606 DO 4 I=IFd,NB 607 AM(I)=AM(1-3) 608 ZETA(I ) = Z E T A ( 1 - 3 ) 609 4 CONTINUE 610 AM(IFE)=EMASSX 611 AM(IFD)=EAd 612 AM(IFK) = EMASSY 613 ZETA(IFE)=EOZETA 614 ZETA(IFD)=EOZETB 615 ZETA(IFK)=EQZETC 616 GOTO 5 617 6 I F ( I F A . E O . 1 ) GOTO 5 618 AM(1)=AM(NBB) 619 AM(2)=AM(NC) 620 AM(3)=AM(NB) 621 ZETA(1 ) = ZETA(NBB) 622 ZETA(2)=ZETA(NC) 623 ZETA(3)=ZETA(NB) 624 5 CONTINUE 625 RETURN 626 END 627 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 628 C C 629 C S U B R O U T I N E M A S S C 630 C C 631 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 632 C 633 C THIS SUBROUTINE SETS UP THE MASS MATRIX AND ITS INVERSE 634 C 635 SUBROUTINE MASS(AM,AMASS,AMASIN,NBV,IFL) 636 IMPLICIT R E A L * 8 ( A - H . O - Z ) 637 DIMENSION AMASS!13, 13),AMASIN! 13, 13),AMI 13) 638 CALL D G S E T ( A M A S S . 1 3 . 1 3 . 1 3 . 0 . D O ) 151 639 CALL D G S E T ( A M A S I N , 1 3 , 1 3 , 1 3 , 0 . D O ) 640 DO 1 1=1,NBV 641 A M A S S f I , I ) = A M ( I ) 642 AMASIN(I , I )=1 .DO/AMASS(1 ,1) 643 1 CONTINUE 644 WRITE(6 ,74) 645 74 FORMAT (// , ' »»»»* MASS MATRIX ***<•*<) 646 CALL D P R M A T ( A M A S S , 1 3 . 1 3 . N B V . N B V , 1 , 1 , 1 3 . 1 ) 647 RETURN 648 END 649 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 650 C C 651 C S U B R O U T I N E S T I F F S C 652 C C 653 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 654 C 655 C THIS SUBROUTINE SETS UP THE STIFFNESS MATRIX 656 C 657 SUBROUTINE S T I F F S ( S K . S T I F F . X C , A . B L . B R . C . S K H 2 . N B . I F L . I F A . I A . 1 0 658 1 ,VERTL ,VERTR,SKH1L .SKH1R.SKV1L .SKV1R,SKV2T,SKV2C,D1,D2 659 1 .SSK.ESK) 660 IMPLICIT R E A L * 8 ( A - H , 0 - Z ) 661 DIMENSION STIFF( 13,13 ) .SK< 1 3 ) , S S K ( 1 0 , 1 0 ) , I A ( 3 ) , I 0 ( 1 8 ) , E S K ( 3 . 3 ) 662 CALL D G S E T ( S T I F F , 1 3 , 1 3 , 1 3 , 0 . D O ) 663 CALL D G S E T ( S S K , 1 0 , 1 0 , 1 0 , 0 . D O ) 664 CALL D G S E T ( E S K . 3 . 3 , 3 . 0 . D 0 ) 665 NC=NB-1 666 NBB=NB-2 667 NBD=NB-3 668 NBS=NB-4 669 IFJ=IFL+4 670 IFK=IFL+3 671 IFD=IFL+2 672 IFE=IFL+1 673 U1=D1*SKV2T*A 674 U2=D2*SKV2C*A 675 SK(NBB)=SKH1L+SKH1R 676 SKDUMY=SKV1L*BL*BL+SKV1R*BR*BR 677 SK(NC)=SKV1L+SKV1R 678 SK(NB)=SKV1R*BR-SKV1L*BL 679 I F ( I F A . N E . O ) GOTO 77 680 C 681 C SET UP A STIFFNESS MATRIX FOR A UNCOUPLED EOUIPMENT ANALYSIS 682 C 683 IF ( D A B S ( X C ) . G E . C ) SK(NBB ) = SKH1L + SKH2+SKH1R + SKH2 684 IF ( V E R T R . L E . ( - D 1 ) ) SKDUMY = (SKV1R+SKV2T)»BR*BR+SKV1L*BL*BL-U1 685 IF ( V E R T R . L E . ( - D 1 ) ) SK(NC)=SKV1R+SKV2T+SKV1L 686 IF ( V E R T R . L E . ( - D 1 ) ) SK(NB) = (SKV1R+SKV2T)*BR-SKV1L*BL 687 IF ( V E R T R . G E . ( D 2 ) ) SKDUMY«(SKV 1R+SKV2C)*BR*BR+5KV1L*BL*BL + U2 688 IF ( V E R T R . G E . ( D 2 ) ) SK(NC ) = SKV1R+SKV2C+SKV1L 689 IF ( V E R T R . G E . ( D 2 ) ) SK(NB ) = (SKV1R+SKV2C)»BR-SKV1L*BL 690 IF ( V E R T L . L E . ( - D 1 ) ) SKDUMY = (SKV1L + SKV2T)*BL *BL + SKV1R*BR*BR-U1 691 IF ( V E R T L . L E . ( - D 1 ) ) SK(NC)=SKV1L+SKV2T+SKV1R 692 IF ( V E R T L . L E . ( - D 1 ) ) SK(NB ) = SKV1R*BR-(SKV1L + SKV2T)*BL 693 IF ( V E R T L . G E . ( D 2 ) ) SKDUMY = (SKV1L + SKV2C)*BL*BL + SKV1R*BR*BR+U2 694 IF ( V E R T L . G E . ( D 2 ) ) SK(NC ) = SKV1L + SKV2C+SKV1R 695 IF ( V E R T L . G E . ( D 2 ) ) SK(NB)=SKV1R*BR-(SKV1L+SKV2C)*BL 696 I F ( ( V E R T R . L E . ( - D 1 ) ) . A N D . ( V E R T L . G E . ( D 2 ) ) ) SKDUMY = (SKV1L+SKV2C)* 152 6 9 7 1 B L * B L + ( S K V 1 R + S K V 2 T ) * B R * B R - U 1 + U 2 6 9 8 I F ( ( V E R T R . L E . ( - D 1 ) ) . A N D . ( V E R T L . G E . ( D 2 ) ) ) S K ( N C ) = S K V 1 L + S K V 2 C + 6 9 9 1 S K V 1 R + S K V 2 T 7 0 0 I F ( ( V E R T R . L E . ( - D 1 ) ) . A N D . ( V E R T L . G E . ( D 2 ) ) ) S K ( N B ) = ( S K V 1 R + S K V 2 T ) * 7 0 1 1 B R - ( S K V 1 L + S K V 2 C ) * B L 7 0 2 I F ( ( V E R T R . G E . ( D 2 ) ) . A N D . ( V E R T L . L E . ( - D 1 ) ) ) S K D U M Y = ( S K V 1 L + S K V 2 T ) * 7 0 3 1 B L * B L + ( S K V 1 R + S K V 2 C ) * B R * B R - U 1 + U 2 7 0 4 I F ( ( V E R T R . G E . ( D 2 ) ) . A N D . ( V E R T L . L E . ( - D 1 ) ) ) S K ( N C ) = S K V 1 L + S K V 2 T + 7 0 5 1 S K V 1 R + S K V 2 C 7 0 6 I F ( ( V E R T R . L E . ( D 2 ) ) . A N D . ( V E R T L . L E . ( - D 1 ) ) ) S K ( N B ) = ( S K V 1 R + S K V 2 C ) * 7 0 7 1 B R - ( S K V 1 L + S K V 2 T ) * B L 7 0 8 E S K ( 1 , 1 ) = S K ( N B B ) 7 0 9 E S K ( 1 , 2 ) = - S K ( N B B ) * A 7 1 0 E S K ( 2 , 1 ) = - S K ( N B B ) * A 7 1 1 E S K ( 2 . 2 ) = S K ( N B B ) * A * A + S K D U M Y 7 1 2 E S K ( 2 . 3 ) = S K ( N B ) 7 1 3 E S K ( 3 . 2 ) = S K ( N B ) 7 1 4 E S K ( 3 . 3 ) = S K ( N C ) 7 1 5 C 7 1 6 C P U T E Q U I P M E N T S T I F F N E S S M A T R I X I N T O S Y S T E M S T I F F N E S S M A T R I X 7 1 7 C 7 1 8 D O 9 1 = 1 . 3 7 1 9 D O 1 0 J = 1 , 3 7 2 0 S T I F F ( I , J ) = E S K ( 1 , 0 ) 7 2 1 1 0 C O N T I N U E 7 2 2 9 C O N T I N U E 7 2 3 G O T O 7 9 7 2 4 C 7 2 5 C P U T S T R U C T U R E S T I F F N E S S E S I N T O A D U M M Y M A T R I X 7 2 6 C 7 2 7 7 7 D O 1 1 I = 1 , N B S 7 2 8 1 1 = 1 + 1 7 2 9 S S K ( I . I ) = S K ( I ) + S K ( I I ) 7 3 0 S S K ( I . I I ) = - S K ( 1 1 ) 7 3 1 S S K ( I I , I ) = - S K ( 1 1 ) 7 3 2 1 1 C O N T I N U E 7 3 3 S S K ( N B D , N B D ) = S K ( N B D ) 7 3 4 I F Z = I F L 7 3 5 I F ( I F L . E Q . O ) I F Z = N B D 7 3 6 C 7 3 7 C I F T H E E Q U I P M E N T I S N O T A T T A C H E D T O T H E S T R U C T U R E . F I L L 7 3 8 C I N T H E S T I F F N E S S M A T R I X W I T H T H E S T R U C T U R A L S T I F F N E S S . 7 3 9 C 7 4 0 D O 1 2 1 = 1 , I F Z 7 4 1 D O 1 3 J = 1 , I F Z 7 4 2 S T I F F ! I . J ) = S S K ( I , J ) 7 4 3 1 3 C O N T I N U E 7 4 4 1 2 C O N T I N U E 7 4 5 I F ( I F L . E Q . O ) G O T O 7 9 7 4 6 I F ( I F L . E Q . N B D ) G O T O 7 8 7 4 7 C 7 4 8 C I F T H E E Q U I P M E N T I S N O T C O N N E C T T O T H E T O P F L O O R . 7 4 9 C R E A R R A N G E T H E S T I F F N E S S A C C O R D I N G L Y . 7 5 0 c 7 5 1 D O 1 4 I = I F L , N B S 7 5 2 D O 1 5 K = 1 , N B D 7 5 3 S T I F F ( K . I + 4 ) = S S K ( K , 1 + 1 ) 7 5 4 1 5 C O N T I N U E 153 755 14 CONTINUE 756 DO 16 I=IFL.NBS 757 DO 17 KK=NBB,NB 758 S T I F F ( I + 4 . K K )=STIFF(1+1,KK) 759 STIFF(1+1,KK)=O.DO 760 17 CONTINUE 761 16 CONTINUE 762 DO 18 1=1,NB 763 DO 19 J=1,NB 764 S T I F F ( J , I ) = S T I F F ( I . J ) 765 19 CONTINUE 766 18 CONTINUE 767 C 768 C PUT EOUIPMENT STIFFNESSES INTO THE STIFFNESS MATRIX. 769 C 770 78 IF (OABS(XC) ,GE .C ) SK(NBB ) =SKH1L + SKH2 + SKH1R + SKH2 771 IF ( V E R T R . L E . ( - D 1 ) ) SKDUMY = (SKV1R+SKV2T)*BR*BR+SKV1L*BL*BL-U1 772 IF (VERTR . LE . ( -D1 ) ) SK(NC ) = SKV1R+SKV2T + SKV1L 773 IF ( V E R T R . L E . ( - D 1 )) SK(NB ) = <SKV 1R+SKV2T)*BR-SKV1L>BL 774 IF (VERTR.GE . (D2 ) ) SKDUMY=(SKV1R+SKV2C)*BR*BR+SKV1L*BL*BL+U2 775 IF ( V E R T R . G E . ( D 2 ) ) SK(NC ) = SKV1R + SKV2C + SKV1L 776 IF ( V E R T R . G E . ( D 2 ) ) SK(NB ) = (SKV1R + SKV2C)*BR-SKV1L*BL 777 IF ( V E R T L . L E . ( - D 1 ) ) SKDUMY = (SKV1L + SKV2T ) *BL*BL+SKV1R*BR*BR-U1 778 IF (VERTL. LE.(-DD) SK( NC ) = SKV 1 L + SKV2T + SKV1R 779 IF (VERTL . LE . (-D1 ) ) SK(NB ) =SKV 1R*BR-(SKV1L+SKV2T)*BL 780 IF ( V E R T L . G E . (D2)) SKDUMY = (SKV1L + SKV2C)*BL*BL+SKV1R*BR*BR + U2 781 IF ( V E R T L . G E . ( D 2 ) ) SK(NC)=SKV1L+SKV2C+SKV1R 782 IF ( V E R T L . G E . ( D 2 ) ) SK(NB ) =SKV1R*BR-(SKV1L + SKV2C)*BL 783 I F ( ( V E R T R . L E . ( - D 1 ) ) - AND . ( V E R T L . G E . ( D 2 ) ) ) SKDUMY = (SKV1L + SKV2C) 784 1BL*BL+(SKV1R+SKV2T)*BR*BR-U1+U2 785 I F ( ( V E R T R . L E . ( - D 1 ) ) . A N D . ( V E R T L . G E . ( D 2 ) ) ) SK(NC)=SKV1L+SKV2C+ 786 1SKV1R+SKV2T 787 I F ( ( V E R T R . L E . ( - D 1 ) ) . A N D . ( V E R T L . G E . ( D 2 ) ) ) SK(NB ) = (SKV1R+SKV2T) 788 1BR-(SKV1L+SKV2C )*BL 789 I F ( ( V E R T R . G E . ( D 2 ) ) . A N D . ( V E R T L . L E . ( - D 1 ) ) ) SKDUMY=(SKV1L+SKV2T) 790 1BL»BL+(SKV1R+SKV2C)*BR*BR-U1+U2 791 IF( ( VERTR.GE. ( D 2 ) ) . A N D . ( V E R T L . L E . ( - D 1 ) ) ) SK(NC) = SKV1L + SKV2T + 792 1SKV1R+SKV2C 793 I F ( ( V E R T R . L E . ( D 2 ) ) . A N D . ( V E R T L . L E . ( - D 1 )) ) SK(NB ) = (SKV1R + SKV2C) 794 1BR»BR-(SKV1L+SKV2T)*BL *BL 795 S T I F F ( I F E , I F E ) = S K ( N B B ) 796 STIFF(IFE,IFD)=-SK(NBB)«A 797 S T I F F ( I F D , I F E ) = - S K ( N B B ) * A 798 S T I F F ( I F D , I F D ) = SK(NBB ) *A*A+SKDUMY 799 S T I F F ( I F D , I F K ) = S K ( N B ) 800 S T I F F ( I F K . I F D ) = S K ( N B ) 801 S T I F F ( I F K . I F K ) = S K ( N C ) 802 S T I F F ( I F L . I F L ) = S T I F F ( I F L . I F L ) + S K ( N B B ) 803 S T I F F ( I F L . I F E ) = - S K ( N B B ) 804 S T I F F ( I F E , I F L ) = - S K ( N B B ) 805 S T I F F ( I F L , I F D ) r SK(NBB)* A 806 S T I F F ( I F D , I F L) = SK(NBB)* A 807 ESK( 1 , 1 ) = SK(NBB) 808 ESK( 1 , 2 ) = -SK(NBB)*A 809 ESK(2 ,1)=-SK(NBB)*A 810 ESK(2,2)=SK(NBB)*A*A+SKDUMY 811 ESK(2 ,3)=SK(NB) 812 ESK(3 ,2)=SK(NB) 154 813 E5K(3.3)=SK(NC) 814 C 815 C OUTPUT STIFFNESS MATRIX IF NEEDED 816 C 817 79 IF(DABS(XC) LT.C) IA(1 )=0 618 IF(DABS(XC).GE.C) IA(1)=1 819 IF (( VERTR. GT.(-DD). AND. (VERTR.LT.(D2))) IA<2)=0 820 IF (VERTR.LE.(-DD) IA(2)=1 821 IF (VERTR.GE.(D2) ) IA(2)=2 822 IF((VERTL GT.(-D1) ) .AND.(VERTL.LT.(D2))) IA(3 )=0 823 IF (VERTL.LE.(-DO) IA(3) = 1 824 IF (VERTL GE.(D2)) IA(3)=2 825 I JK= 1 826 IF(IA(1 ) .EO.1) IJK=10 827 IF(IA(2).EO.O IJK=IJK+3 828 IF(IA(2).E0.2) IJK=IJK+6 829 IF(IA(3).EO.O) IQ(IJK)=IQ(IJK)+1 830 IF((IA(3).E0.0).AND.(I0(IJK).E0.O) GOTO 81 831 IF(IA(3) .EO.1) I0(IJK+1) = I0(IJK+1) + 1 832 IF((IA(3).E0.O.AND.(10(1 JK+O.E0.1)) GOTO 81 833 IF(IA(3) EQ.2) I0(IJK+2)=IQ(IJK+2)+1 834 IF((IA(3).EO.2).AND.(I0(IdK+2).EO.1)) GOTO 81 835 GOTO 82 836 81 WRITE(6.84) 837 84 FORMAT(//, ' STIFFNESS MATRIX ***•*') 838 NBV=NB 839 IF(IFA.EO. 1 ) NBV=NB-3 840 IF(IFA.EO.O) NBV = 3 84 1 CALL DPRMAT(STIFF, 13, 13,NBV.NBV. 1,1.13,1) 842 82 CONTINUE 843 RETURN 844 END 845 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 846 c C 847 c SUBROUTINE FORCE C 848 c C 849 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 850 c 851 c THIS SUBROUTINE SETS UP THE INCREMENTAL FORCE VECTOR 852 c 853 SUBROUTINE FORCE(AFF,DETAT,DETAF.AMASS.I,N,FA,XC.A,SKH2,C, 854 1NBV,IFL,I FLAG,JE.IT,A 1,AS.IFA,SKV2T,SKV2C,D1,D2,VERTL.VERTR, 855 1,DETVL,DETVR,ROR,DETAX,FMU,DETAXC) 856 IMPLICIT REAL*8(A-H.O-Z) 857 DIMENSION AFF(60000).DETAF( 13) ,FA( 13) .AMASS( 13, 13 ),DETAX( 13 ) 858 11=1+1 859 A2 = AFF(11 ) 860 A 1 =AFF( I ) 861 IF(IFLAG.E0.2) GOTO 8 862 IF(IFLAG.EO.O) GOTO 5 863 c 864 c IF ITERATION IS NEEDED. USE THE APPROPRIATE PORTION OF 865 c THE INPUT ACCELERATON BASED ON THE ITERATION TIME STEP 866 c AND THE ASSUMPTION OF A LINEAR VARIATION IN INPUT 867 c ACCELERATION. 868 c 869 A2=AS*DETAT+A1 870 GOTO 5 155 871 C 872 C AFTER THE ITERATION PROCEDURE IS COMPLETED USE UP THE 873 C REMAINING INPUT ACCELERATION IN ORDER TO COMPLETE THE 874 C TIME STEP. 875 C 876 8 A1=A2 877 A 2 = A F F ( I I ) 878 5 DO 3 J=1.NBV 879 D E T A F ( d ) = - ( A 2 - A 1 ) * A M A S S ) J , J ) 880 F A ( J ) = - A M A S S ( J , J ) * A 1 881 3 CONTINUE 882 I F ( I FA.EO•1 ) GOTO 4 883 DUMMY=0.D0 884 DUMMY 1=0.DO 885 DUMMY2=0.D0 886 DUMMY3=0.DO 887 DUMMY4=0.D0 888 IF(DABS(XC) . G E . C ) DUMMY = DSIGN(2.D0*SKH2*C,XC) 889 I F ( D A B S ( X C ) . G E . C ) FFR=DSIGN(FMU'DETAXC*SKH2,DETVR) 890 I F ( D A B S ( X C ) . G E . C ) F F L = DSIGN(FMU*DETAXC*SKH2,DETVL) 891 I F ( D A B S ( X C ) . G E . C ) DUMMY3=FFL*BL-FFR*BR+(FFL+FFR)»A*ROR 892 I F ( D A B S ( X C ) . G E . C ) DUMMY4=FFL+FFR 893 IF ( V E R T R . L E . ( - D 1 ) ) DUMMY 1=-BR*D1*SKV2T 894 IF ( V E R T R . G E . ( D 2 ) ) DUMMY1=BR*D2*SKV2C 895 IF ( V E R T L . L E . ( - D I ) ) DUMMY2 = BL*D1*SKV2T 896 IF (VERTL G E . ( D 2 ) ) DUMMY2=-BL*D2*SKV2C 897 I F ( I F A . E O . O ) GOTO 6 898 F A ( I F L ) = - A M A S S ( I F L , I F L ) * A 1-DUMMY 899 6 FA(IFL+1) = -AMASS(IFL+1.IFL+1)*A 1+DUMMY 900 F A ( I F L + 2) = - DUMMY * A+DUMMY1+DUMMY 2+DUMMY 3 901 FA(IFL+3)=DUMMY1/BR-DUMMY2/BL-DUMMY4 902 DETAF(IFL+2)=0.D0 903 DETAF(IFL+3)=0.DO 904 4 CONTINUE 905 RETURN 906 END 907 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 908 c C 909 c S U B R O U T I N E P F O R C E C 910 c c 911 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 912 c 913 c THIS SUBROUTINE SETS UP THE PSEUDO FORCE VECTOR 914 c 915 SUBROUTINE PFORCE(DETAF,DETAT,AMASS,AC,V ,PSEUFO.DC,NBV , BET A 916 1,ALBE,ALBR ) 917 IMPLICIT REAL*8(A-H .O-Z ) 918 DIMENSION D E T A F ( 1 3 ) . A M A S S ( 1 3 , 1 3 ) , A C ( 1 3 ) . V ( 1 3 ) . P S E U F O ( 1 3 ) 919 DIMENSION DUM1(13),DUM( 13) .DC( 13, 13) ,DUM2( 13 ) ,DUM3( 13) 920 DO 1 1=1,NBV 921 DUM(I)=0.5D0*AC(I)/BETA+V(I )/DETAT/BETA 922 1 CONTINUE 923 CALL DGMATV(AMASS,DUM.DUM1,NBV.NBV.13) 924 DO 2 1=1,NBV 925 DUM2(I)"ALBE*V(I ) -AC(I ) *DETAT * ALBR 926 2 CONTINUE 927 CALL DGMATV(DC,DUM2,DUM3,NBV,NBV,13) 928 DO 3 1=1.NBV 156 929 P S E U F O U )=DETAF(I)+DUM1(I)+DUM3<I) 930 3 CONTINUE 931 RETURN 932 END 933 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 934 C C 935 C S U B R O U T I N E P S T I F F C 936 C C 937 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 938 C 939 C THIS SUBROUTINE SETS UP THE PSEUDO STIFFNESS MATRIX 940 C 941 SUBROUTINE P S T I F F ( S T I F F , D E T A T , A M A S S , P S E U S T , D C . N B V , B E T A . A L B E ) 942 IMPLICIT R E A L * 8 ( A - H . O - Z ) 943 DIMENSION S T I F F ( 1 3 . 1 3 ) , A M A S S ( 1 3 . 1 3 ) . P S E U S T ( 1 3 . 1 3 ) . D C ( 1 3 . 1 3 ) 944 DO 1 1=1.NBV 945 DO 2 d =1,NBV 946 P S E U S T ( I . J ) = A M A S S ( I . J ) / ( D E T A T * D E T A T * B E T A ) + S T I F F ( I . J ) + 947 1 A L B E * D C ( I , d l/DETAT 948 2 CONTINUE 949 1 CONTINUE 950 RETURN 951 END 952 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 953 C C 954 C S U B R O U T I N E S T E P C 955 C C 956 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 957 C 958 C THIS SUBROUTINE CALCULATES THE EIGENVECTORS. EIGENVALUES AND 959 C MODAL PATICIPATION FACTORS OF THE SYSTEM OF EQUATIONS TO BE 960 C SOLVED SO THAT A CHECK ON THE TIME STEP CAN BE MADE. IT USES 961 C A UBC LIBARY SUBROUTINE TO SOLVE THE EIGENVALUE PROBLEM. 962 C 963 SUBROUTINE S T E P ( S T I F F , A M A S S . C H E C K . D E . D V . X C . C . N B V , I F L . I A . J Q 964 1 .VERTR ,VERTL ,D1 .D2 .SDE ,SDV .EDE .EDV .SMAS .EMAS .SSK .ESK .NBD) 965 IMPLICIT R E A L * 8 ( A - H , 0 - Z ) 966 DIMENSION S T I F F ( 1 3 , 1 3 ) , A M A S S ( 13, 13),DB(91 ),DA( 169),DE( 13) 967 DIMENSION DV( 13. 13) ,GAMMA(13) ,VI(13) ,DUM( 13).DUM1( 13) .VV( 13) 968 DIMENSION E S K ( 3 . 3 ) , S D A ( 1 0 0 ) . E D A ( 9 ) , S D B ( 5 5 ) . E D B ( 6 ) . S D E (10) 969 DIMENSION EDE(3 ) ,SDV( 10. 10 ) , E D V ( 3 , 3 ) . S S K ( 1 0 , 1 0 ) 970 DIMENSION SMAS( 1 0 . 1 0 ) . E M A S ( 3 , 3 ) . I A ( 3 ) . J Q ( 18) 971 CALL D G S E T ( D V . 1 3 , 1 3 , 1 3 . O . D 0 ) 972 CALL D G S E T ( S D V , 1 0 , 1 0 . 1 0 . 0 . D O ) 973 CALL D G S E T ( E D V , 3 . 3 . 3 , O . D O ) 974 KK=0 975 DO 1 1=1.NBV 976 DO 2 K= 1 , I 977 KK=KK+1 978 DA(KK ) = S T I F F ( I , K ) 979 DB(KK )=AMASS(I,K) 980 2 CONTINUE 981 1 CONTINUE 982 CALL DRSGAL(DA.DB.NBV,DE.I ERROR.1) 983 DO 3 1 = 1 ,NBV 984 IOI = (1-1 )*NBV+1 985 DO 4 J= 1,NBV 986 D V ( J . I )=DA(IOI ) 157 937 101=101+1 988 4 CONTINUE 989 3 CONTINUE 990 KK=0 991 DO 11 1=1,NBD 992 DO 12 K=1,I 993 KK=KK+1 994 SDA(KK)=SSK(I ,K) 995 SDB(KK)=SMAS(I ,K) 996 12 CONTINUE 997 11 CONTINUE 998 CALL DRSGAL(SDA,SDB.NBD.SDE.IERROR,1) 999 DO 13 1=1,NBD 1000 IOI=(1-1)*NBD+1 1001 DO 14 J=1,NBD 1002 SDV(J,I)=SDA(101) 1003 101=101+1 1004 14 CONTINUE 1005 13 CONTINUE 1006 KK=0 1007 DO 21 1=1,3 ' 1008 DO 22 K=1,I 1009 KK=KK+1 1010 EDA(KK)=ESK(I ,K) 1011 EDB(KK ) = EMAS(I ,K) 1012 22 CONTINUE 1013 21 CONTINUE 1014 CALL DRSGAL(EDA,EDB,3 ,EDE,IERROR,1 ) 1015 DO 23 1=1,3 1016 I0I=(I -1)*3+1 1017 DO 24 J=1 ,3 1018 EDV(J ,I)=EDA(101) 1019 101=101+1 1020 24 CONTINUE 1021 23 CONTINUE 1022 C 1023 C OUTPUT NATURAL PERIODS, EIGENVALUES AND EIGENVECTOR IF NEEDED 1024 C 1025 I F ( D A B S ( X C ) . L T . C ) IA(1)=0 1026 I F ( D A B S ( X C ) . G E . C ) IA(1)=1 1027 I F ( ( V E R T R . G T . ( - D 1 ) ) . A N D . ( V E R T R . L T . ( D 2 ) ) ) IA(2)=0 1028 IF ( V E R T R . L E . ( - D 1 ) ) IA(2)=1 1029 IF ( V E R T R . G E . ( D 2 ) ) IA(2)=2 1030 I F ( ( V E R T L . G T . ( - D 1 ) ) . A N D . ( V E R T L . L T . ( D 2 ) ) ) IA(3)=0 1031 IF ( V E R T L . L E . ( - D 1 ) ) IA(3)=1 1032 IF ( V E R T L . G E . ( D 2 ) ) IA(3)=2 1033 IJK=1 1034 I F ( I A ( 1 ) . E O . 1 ) IJK=10 1035 I F ( I A ( 2 ) EO .1) IJK=IJK+3 1036 I F ( I A ( 2 ) . E 0 . 2 ) IJK=IJK+6 1037 I F ( I A ( 3 ) . E O . O ) J O ( I J K ) = JO(IJK)+1 1038 I F ( ( I A ( 3 ) . E O . O ) . A N D . ( J 0 ( I J K ) . E O . 1)) GOTO 81 1039 I F ( I A ( 3 ) .EO. 1) J0(IJK+1)=J0(IJK+1)+1 1040 I F ( ( I A ( 3 ) . E O . 1 ) . A N D . ( J 0 ( I J K + 1 ) .EO. 1)) GOTO 81 1041 I F ( I A ( 3 ) . E O . 2 ) J 0 ( I J K + 2 ) = J0(IJK+2)+1 1042 I F ( ( I A ( 3 ) . E Q . 2 ) . A N D . ( J O ( I J K + 2 ) . E O . O ) GOTO 81 1043 GOTO 82 1044 C 158 1045 C COMPUTE THE MODAL PARTICIPATION FACTORS 1046 C 1047 81 DO 5 1=1,13 1048 GAMMA(I)=0.D0 1049 VI(I)=O.DO 1050 VV(I)=0.D0. 1051 5 CONTINUE 1052 DO 6 1=1,NBV 1053 VI(I ) = AMASS(I,I) 1054 6 CONTINUE 1055 DO 7 1=1,NBV 105G DO 8 J=1,NBV 1057 VV(J)=DV(J.I ) 1058 8 CONTINUE 1059 VB=DGVV(VV,VI.NBV) 1060 CALL DGMATV(AMASS,VV,DUM1,NBV.NBV, 13) 1061 VD=DGVV(VV.DUM1.NBV) 1062 GAMMA(I)=VB/VD 1063 7 CONTINUE 1064 WRITE(6.77) 1065 77 F0RMAT(//,5X, 'EIGENVALUES' ,7X, 'NATURAL PERIODS' ,6X. 1066 1 'MODAL PARTICIPATION FACTORS) 1067 WRITE(6,78) 1068 78 FORMAT(/ ) 1069 DO 79 1=1,NBV 1070 PERI0D=(2.D0*3.141592654D0)/DSORT(DE(I)) 1071 WRITE(6,80)DE(I),PERIOD,GAMMA(I ) 1072 80 F0RMAT(3F20.5) 1073 79 CONTINUE 1074 82 CHECK=(0.2D0*3.14 15926540O)/DSORT(DE(NBV)) 1075 RETURN 1076 END 1077 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1078 c C 1079 c SUBROUTINE DAMP C 1080 c C 108 1 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1082 c 1083 c THIS SUBROUTINE SETS UP THE DAMPING MATRIX ASSUMING 1084 c CLASSICAL DAMPING. 1085 c 1086 SUBROUTINE DAMP(ZETA,AMASS.DE.DC.DV,XC.C.NBV.IFL.ST IFF,IA.KQ 1087 1.VERTR,VERTL.D1.D2,SDC,EDC,SDV.EDV,NB.EZETA,SZETA.IFA,SDE,EDE 1088 1,EMAS,SMAS ) 1089 IMPLICIT REAL*8(A-H.O-Z) 1090 DIMENSION DC( 13. 13).AMASS(13.13 ) ,DE( 13 ) ,DV(13. 13 ) ,SZ(10. 101 1091 DIMENSION VDUM(13.13),VDUM1(13,13),VDUM2(13,13),EZ(3.3) 1092 DIMENSION EVDUM(3.3).EVDUM1(3,3) ,EVDUM2(3.3).ETDV(3,3) 1093 DIMENSION SVDUM(10.10),SVDUM1(10,10),SVDUM2(10,10),EZETA(3.3) 1094 DIMENSION STDV(10.10),TDV(13,13),ZETA(13).STIFF(13.13).IA(3) 1095 DIMENSION K0(18),SDC( 10. 10).EDC(3.3),SDV( 10, 10) .EDV(3.3) 1096 DIMENSION SDE(10),EDE(3),EMAS(3.3).SMASf10.10).SZETA(10.10) 1097 CALL DGSET(DC. 13. 13. 13.0.DO) 1098 CALL DGSET(SDC.10,10.10,0.DO) 1099 CALL DGSET(EDC.3.3,3.0.DO) 1 100 CALL DGSETCSZ.10.10.10.0.D0) 1 101 CALL DGSET(EZ.3,3.3.0.DO) 1 102 NC=NB-1 159 1103 NBB=NB-2 1104 NBD=NB-3 1105 NB5=NB-4 1106 IFd=IFL+4 1107 IFK=IFL+3 1108 IFD=IFL+2 1109 IFE=IFL+1 11 10 C 1111 C FIND STRUCTURE STIFFNESS MATRIX ASSUMING CLASSICAL DAMPING 1112 C 1113 DO 1 1=1,NBD 1114 S Z ( I , I ) = 2 . D 0 * S Z E T A ( I . I ) * ( D S Q R T ( S D E ( I ) ) ) 1115 1 CONTINUE 1116 CALL DGTRAN(SDV,STDV,NBD.NBD,10,10) 1117 CALL DGMULT(STDV,SMAS,SVDUM,NBD,NBD,NBD,10,10,10) 1118 CALL DGMULT(SZ,SVDUM,SVDUM1.NBD,NBD,NBD,10,10,10) 1119 CALL DGMULT(SDV,SVDUM1.SVDUM2.NBD,NBD,NBD,10,10,10) 1120 CALL DGMULT(SMAS,SVDUM2.SDC,NBD,NBD,NBD,10,10,10) 1121 C 1122 C FIND EOUIPMENT STIFFNESS MATRIX ASSUMING CLASSICAL DAMPING 1 123 C 1 124 DO 2 1 = 1 .3 1125 E Z ( I , I ) = 2 . D 0 * E Z E T A ( I , I ) * ( D S O R T ( E D E ( I ) ) ) 1126 2 CONTINUE 1127 CALL D G T R A N ( E D V , E T D V , 3 , 3 , 3 . 3 ) 1128 CALL DGMULT(ETDV .EMAS ,EVDUM,3 ,3 ,3 ,3 ,3 ,3) 1129 CALL DGMULT(EZ ,EVDUM,EVDUM1,3 ,3 .3 ,3 ,3 ,3) 1130 CALL DGMULT(EDV,EVDUM1,EVDUM2,3 .3 ,3 ,3 ,3 ,3) 1131 CALL DGMULT(EMAS ,EVDUM2,EDC ,3 .3 .3 ,3 ,3 ,3) 1132 I F ( I F A . E O . O ) GOTO 3 1133 GOTO 4 1 134 C 1135 C BUILT DAMPING MATRIX 1136 C 1137 3 DO 6 1 = 1 , 3 1138 DO 7 d=1 ,3 1139 O C ( I . d ) = E D C ( I , d ) 1140 7 CONTINUE 1141 6 CONTINUE 1 142 GOTO 80 1143 4 IFZ=IFL 1144 I F ( I F A . E O . I ) IFZ=NBD 1 145 C 1146 C IF THE EOUIPMENT IS NOT ATTACHED TO THE STRUCTURE. FILL 1147 C IN THE DAMPING MATRIX WITH THE STRUCTURAL DAMPING. 1148 C 1149 DO 8 1 = 1 , IFZ 1150 DO 9 d=1 , IFZ 1151 D C ( I , d ) = S D C ( I , d ) 1152 9* CONTINUE 1153 8 CONTINUE 1154 I F ( I F A . E O . I ) GOTO 80 1155 IF ( I F L . E O . N B D ) GOTO 79 1 156 C 1157 C IF THE EOUIPMENT IS NOT CONNECT TO THE TOP FLOOR, 1158 C REARRANGE THE DAMPING ACCORDINGLY. 1 159 C 1160 DO 10 I=IFL,NBS 160 1 161 DO 11 K=1,NBD 1 162 DC(K,I+4)=SDC(K.1+1) 1 163 1 1 CONTINUE 1 164 10 CONTINUE 1 165 DO 12 I=IFL.NBS 1 166 DO 13 KK=NBB,NB 1 167 DC(I+4,KK)=DC(1+1.KK) 1 168 DC(I+1,KK)=0.D0 1169 13 CONTINUE 1170 12 CONTINUE 1 171 DO 14 I=1.NB 1 172 DO 15 d = 1 ,NB 1 173 DC(d,I)=DC(I.J) 1 174 15 CONTINUE 1 175 14 CONTINUE 1 176 79 DC(IFE.IFE)=EDC(1.1) 1 177 DC(IFE,IFD)=EDC(1,2) 1 178 DC(IFE,IFK)=EDC(1,3) 1 179 DC(IFD,IFE)=EDC(2, 1 ) 1 180 DC(IFD,IFD) = EDC(2.2 ) 1 181 DC(IFD,IFK)=EDC(2,3) 1 182 DC(IFK,IFE)=EDC(3. 1 ) 1 183 DC(IFK,IFD)=EDC(3.2) 1 184 DC(IFK,IFK)=EDC(3,3) 1 185 DC(IFL,IFL)=DC(IFL,IFL)+EDC(1,1) 1 186 DC(IFL,IFE)=-EDC(1.1) 1 187 DC(IFL,IFD)=-EDC(1,2) 1 188 DC(IFL,IFK)=-EDC(1,3) 1 189 DC(IFE,IFL)*-EDC(1.1) 1 190 DC(IFD,IFL)=-EDC(2.1) 1191 DC(IFK,IFL)=-EDC(3.1) 1 192 C 1 193 C OUTPUT DAMPING MATRIX IF NEEDED 1 194 c 1 195 80 IF(DABS(XC) LT.C) IA(1)=0 1 196 IF(DABS(XC).GE.C) IA(1)=1 1 197 IF((VERTR.GT.(-D1)).AND.(VERTR.LT.(D2))) IA(2)= 0 1 198 IF (VERTR.LE.(-D1)) IA(2)=1 1 199 IF (VERTR.GE.(D2)) IA(2)=2 120O IF((VERTL.GT.(-D1)).AND.(VERTL.LT.(02))) IA(3)= 0 1201 IF (VERTL.LE.(-DD) IA(3) = 1 1202 IF (VERTL.GE.(D2)) IA(3)=2 1203 IJK=1 1204 IF(IA(1).EO.1) IJK=10 1205 IF(IA(2).EO.1) IJK=IdK+3 1206 IF(IA(2) .E0.2) IJK = IdK+G 1207 IF(IA(3).EO.O) K0(IdK)=KQ(IdK)+1 1208 IF((IA(3).EO.O).AND.(K0(IdK).EO.1)) GOTO 81 1209 IF(IA(3) .EO.1) KQ(IJK+1)=K0(IdK+1)+1 1210 IF((IA(3 ) .EO. 1 ) .AND.(K0(IdK+1).EO. 1)) GOTO 81 121 1 IF(IA(3) .EO.2) K0(IdK+2)=K0(IdK+2)+1 1212 IF((IA(3).EO.2).AND.(K0(IJK+2).EO.1)) GOTO 81 1213 GOTO 82 1214 81 WRITE(6,83) 12 15 83 FORMAT(//,' * * * * * EIGEN VECTOR •»***' ) 1216 CALL DPRMAT(DV. 13. 13,NBV,NBV. 1.1.13,1) 1217 WRITE(6.84) 1218 84 FORMAT(//.' * * * * * DAMPING MATRIX * * * * + 161 1219 CALL DPRMAT(DC . 13, 13,NBV,NBV. 1 . 1 , 1 3 . 1 ) 1220 82 CONTINUE 1221 RETURN 1222 END 1223 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1224 c C 1225 c S U B R O U T I N E D R I V E C 1226 c C 1227 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1228 c 1229 c THIS SUBROUTINE IS THE DRIVER OF THE PROGRAM. IT CALCULATES 1230 c THE DISPLACEMENTS, VELOCITIES , ACCELERATIONS. AND FORCES. 1231 c IT USES A UBC LIBARY SUBROUTINE (SLE) TO SOLVE THE SYSTEM 1232 c OF EQUATIONS. 1233 c 1234 SUBROUTINE D R I V E ( X , S K , A M A S S , N B . I F L , A , B L , B R , C . S K H 2 , D E T A T , 1235 1 A M , A F F , I , N . N B B . N C , I T , V . A C , I ERR ,ABC ,ZETA ,DETAX,TIME ,NBV.UE . 1236 11 F A , I C O U N T . X C C , I E R . A M A S I N , A S , I F L A G , X C , D X 1 , D X 2 , T 0 L , B T , 1237 1 D E T A F S , I F L A A , I X X , I D I R E C , A L P H A , B E T A , A L B E . A L B R . S D C , E D C 1238 1 ,VERTL ,VERTR,SKH1L ,SKH1R,SKV1L ,SKV1R,SKV2T,SKV2C,D1,D2 1239 1 , V B L , V B R . S P F H . S P F V R . S P F V L , I A , 1 0 . J Q . K Q , I F L A A A , E M A S , S M A S 1240 1 , E Z E T A , S Z E T A , N B D , E , B , Y A L , Y A R , F U , D E T V L , D E T V R . F M U ) 1241 IMPLICIT R E A L * 8 ( A - H , 0 - Z ) 1242 DIMENSION STIFF( 13, 13 ) .AMASS( 13, 13),AMASIN( 1 3 , 1 3 ) . X ( 1 3 ) . V ( 1 3 ) 1243 DIMENSION DETAF(13) .PSEUF0( 13),PSEUST( 1 3 , 1 3 ) , D E T A X ( 1 3 ) 1244 DIMENSION F I ( 1 3 ) . F S ( 1 3 ) . F A ( 1 3 ) , A C ( 13) ,DETAV(13) .ACF(96013) 1245 DIMENSION D B ( 9 1 ) , D A ( 1 6 9 ) , D E ( 1 3 ) , A F F ( 6 0 0 0 0 ) , A F ( 7 6 0 0 ) , E D C ( 3 . 3 ) 1246 DIMENSION DC(13 . 13 ) ,FD( 13) ,DV(13 , 13 ) .DETX( 13 ) ,ABC(13 ) 1247 DIMENSION A M ( 1 3 ) , S K ( 1 3 ) ,ZETAI 13) . I P E R M ( 2 4 ) , T ( 1 3 , 1 3 ) , S D C ( 1 0 . 1 0 ) 1248 DIMENSION I A ( 3 ) . I Q ( 1 8 ) . J Q ( 1 8 ) , K Q ( 1 8 ) , S M A S ( 1 0 . 1 0 ) , E M A S ( 3 . 3 ) 1249 DIMENSION E S K ( 3 , 3 ) , S D A ( 1 0 0 ) , E D A ( 9 ) . S D B l 5 5 ) , E D B ( 6 ) . S D E ( 1 0 ) 1250 DIMENSION EDE(3) ,SDV( 10, 1 0 ) , E D V ( 3 . 3 ) , S S K ( 1 0 , 10) .SZETA( 10. 10) 1251 DIMENSION E Z E T A ( 3 , 3 ) , D E T A A ( 1 3 ) , FU(13) 1252 c 1253 c FIND SPRING DISPLACEMENT 1254 c 1255 IC=IFL+1 1256 ID=IFL+2 1257 I F ( I F A . E Q . O ) X C = X ( I C ) - A * ( X ( I D ) ) 1258 I F ( I F A . N E . 1 ) VERTL=-BL*X(ID)+X(ID+1) 1259 I F ( I F A . N E . 1 ) VERTR=BR*X(ID)+X(ID+1) 1260 I F ( I F A . N E . 1) YAL = -( .5842D0-E*B/100.DO)*ABC(ID) + ABC ( ID+1) 1261 I F ( I F A . N E . 1 ) YAR=(.5842D0+E*B/100.DO)*ABC(ID)+ABC(ID+1) 1262 I F ( I F A . E O . 1 ) XC=O.DO 1263 I F ( I F A . E Q . 1 ) VERTL=O.DO 1264 I F ( I F A . E Q . 1 ) VERTR=O.DO 1265 I F ( I F L . N E . O ) X C = X ( I C ) - X ( I F L ) - A * ( X ( I D > ) 1266 c 1267 c CHECK IF NEW STIFFNESS AND DAMPING MATRIX ARE NEEDED 1268 c 1269 I F ( ( V B L . L T . D 2 ) . A N D . ( V E R T L . G E . D 2 ) ) IF LAAA=1 1270 I F ( ( V B L . G E . D 2 ) . A N D . ( V E R T L . L T . D 2 ) ) IF LAAA-1 1271 I F ( ( V B L . G T . ( - D 1 ) ) . A N D . ( V E R T L . L T . ( - D 1 ) ) ) IFLAAA= 1 1272 I F ( ( V B L . L E . ( - D 1 ) ) . A N D . ( V E R T L . G T . ( - D 1 ) ) ) IF LAAA= 1 1273 I F ( ( V B R . L T . D 2 ) . A N D . ( V E R T R . G E D2) ) I FLAAA= 1 1274 I F ( ( V B R . G E , D 2 ) . A N D . ( V E R T R . L T .D2 ) ) IF LAAA= 1 1275 IF( ( V B R . G T . ( - D 1 ) ) . A N D . ( V E R T R . L T . ( - D I ) ) ) IFLAAA = 1 1276 I F ( ( V B R . L E . ( - D 1 ) ) . A N D . ( V E R T R . G T . ( - D 1 ) ) ) I FLAAA= 1 162 1277 I F ( ( I F L A A . N E . 2 ) . O R . ( I F L A A A . E Q . 1 ) ) GOTO 75 1278 GOTO 76 1279 75 IF ( I F L A A . N E . 2 ) I FLAA = IFLAA+1 1280 C 1281 C SET UP STIFFNESS MATRIX 1282 C 1283 CALL S T I F F S ( S K . S T I F F . X C . A . B L , B R . C . S K H 2 , N B , I F L . I F A . I A . 1 0 1284 1 .VERTL .VERTR.SKH1L .SKH1R.SKV1L .SKV1R.SKV2T ,SKV2C .D1 .02 1285 1 .SSK .ESK) 1286 C 1287 C CHECK IF THE TIME INCREMENT IS TOO LARGE 1288 C 1289 CALL STEP(ST I F F . A M A S S . C H E C K . D E . D V . X C . C . N B V , I F L . I A. JO 1290 1 .VERTR .VERTL .D1 .D2 .SDE .SDV .EDE .EDV .SMAS .EMAS .SSK .ESK .NBD) 1291 I F ( C H E C K . G E . B T ) GOTO 77 1292 WRITE(5 .61) 1293 WRITE(5 ,62)CHECK,BT 1294 WRITE(6 ,62)CHECK,BT 1295 WRITE(6 ,61) 1296 62 FORMAT('MAX TIME INCREMENT = ' . F 1 5 . 8 , ' DETAT = ' , F 1 5 . 8 ) 1297 61 FORMAT('***»*• TIME INCREMENT IS TOO LARGE ****»**') 1298 STOP 1299 C 1300 C SET UP DAMPING MATRIX 1301 C 1302 77 CALL D A M P ( Z E T A , A M A S S , D E , D C , D V . X C . C , N B V . I F L . S T I F F , I A . K Q 1303 1 . V E R T R , V E R T L . D 1 , D 2 , S D C , E D C , S D V . E D V , N B . E Z E T A , S Z E T A , I F A , S D E , E D E 1304 1.EMAS,SMAS) 1305 76 CONTINUE 1306 C 1307 C SET UP THE PSEUDO STIFFNESS MATRIX 1308 C 1309 CALL P S T I F F ( S T I F F . D E T A T , A M A S S . P S E U S T , D C . N B V , B E T A , A L B E ) 1310 C 1311 C SET UP FORCE INCREMENT VECTOR 1312 C 1313 CALL F O R C E ( A F F . D E T A T . D E T A F . A M A S S . I . N , F A . X C , A , S K H 2 , C , N B V , I F L 1314 1 , I F L A G . J E , I T . A 1 , A S , I F A , S K V 2 T , S K V 2 C . D 1 .D2 .VERTL .VERTR ,BL .BR 1315 1.DETVL,DETVR,ROR,DETAX,FMU,DETAXC) 1316 C 1317 C COMPUTE THE ERROR IN THE SPRING FORCE INTRODUCED BY THE 1318 C ABRUPT CHANGE IN SPRING STIFFNESS 1319 c 1320 I F ( I F A . E O . 1 ) GOTO 799 1321 VBL = - B L * ( X ( I D ) - D E T A X ( I D ) )+(X(ID+1)-DETAX(ID+1)) 1322 VBR=BR*(X(ID)-DETAX(ID))+(X(ID+1)-DETAX(ID+1)) 1323 IF( IFA . E0..0) XB = X ( I C ) - D E T A X ( I C ) - A * ( X ( I D ) - D E T A X ( I D ) ) 1324 I F ( I F L . N E . 0 ) XB = X ( I C ) - D E T A X ( I C ) - X ( I F L ) - D E T A X ( I F L ) - A * ( X ( I D 1325 1 ) -DETAX(ID)) 1326 I F ( ( D A B S ( X B ) . L T . C ) .AND.(DABS(XC) G E . C ) ) ICOUNT = ICOUNT+1 1327 I F ( I F L A G . N E . 2 ) GOTO 799 1328 GAP=C 1329 I F ( ( X B . L T . O . D O ) . A N D . ( X C . L T . 0 . D O ) ) GAP=-C 1330 IF((DABS(XB ) .GT C ) .AND. (DABS(XCI .LE C ) ) GOTO 797 1331 DETAFS=SKH2*(XC-GAP) 1332 GOTO 799 1333 797 DETAFS=-SKH2*(XC-GAP) 1334 c 163 1335 C APPLY EQUILBRUIM CORRECTION PROCEDURE 1336 C 1337 799 CALL D G M A T V ( S T I F F , X , F S , N B V , N B V , 1 3 ) 1338 CALL DGMATV(DC,V ,FD,NBV,NBV,13) 1339 CALL DGMATV(AMASS,AC,FI ,NBV,NBV,13) 1340 CALL D G A D D ( F I , F D , F I . N B V , 1 . 1 3 . 1 3 . 13) 1341 CALL D G A D D ( F I . F S . F I . N B V . 1 . 13, 13 .13) 1342 CALL D G S U B ( F A . F I . F U . N B V , 1 , 1 3 . 1 3 . 1 3 ) 1343 CALL D G A D D ( D E T A F , F U , D E T A F , N B V . 1 , 1 3 , 1 3 . 1 3 ) 1344 C 1345 C FIND SPRING FORCE 1346 C 1347 DUMX=O.DO 1348 IF(DABS(XC) . G E . C ) DUMX=DSIGN(SKH2*C,XC) 1349 SPFH=FS(IC)-DUMX 1350 SPFVR=VERTR*SKV1R 1351 SPFVL=VERTL*SKV1L 1352 IF ( V E R T R . L E . ( - D 1 ) ) SPFVR = (SKV2T + SKV1R)*VERTR+SKV2T*D1 1353 IF ( V E R T R . G E . ( D 2 ) ) SPFVR=(SKV2C+SKV1R)*VERTR-SKV2C*D2 1354 IF ( V E R T L . L E . ( - D 1 ) ) SPFVL=(SKV2T+SKV1L)*VERTL+SKV2T*D1 1355 IF ( V E R T L . G E . ( D 2 ) ) SPFVL=(SKV2C+SKV1L)*VERTL-SKV2C*D2 1356 I F ( I F A . E Q . I ) SPFH=O.DO 1357 C 1358 C SET UP THE PSEUDO FORCE VECTOR 1359 C 1360 CALL PFORCE(DETAF,DETAT,AMASS,AC,V ,PSEUFO,DC,NBV,BETA 1361 1 ,ALBE,ALBR) 1362 C 1363 C SOLVE FOR DISPLACMENT INCREMENT DETAX BY SUBROUTINE (SLE) 1364 C 1365 CALL SLE(NBV, 13.PSEUST, 1, 13,PSEUFO,DETX,IPERM, 1 3 , T , D E T , J E X P ) 1366 I F ( D E T ) 2 1 5 , 2 6 , 2 1 5 1367 215 I F ( I F A . E O . 1 ) GOTO 231 1368 I F ( I F A . E Q . O ) XCC = X(IC )+DETX(IC ) -A*(X(ID)+DETX(ID)) 1369 I F ( I F L . N E . O ) XCC = X(IC ) + D E T X ( I C ) - X ( I F L ) - D E T X ( I F L ) - A * 1370 1 (X(ID )+DETX(ID)) 1371 C 1372 C CHECK FOR ABRUPT CHANGE IN SPRING STIFFNESS, THE LOADING 1373 C DIRECTION AND PREPARE FOR ITERATION. 1374 C 1375 I F ( ( D A B S ( X C ) . L T . C ) . A N D . ( D A B S ( X C C ) . G T . C ) ) IER=1 1376 I F ( ( D A B S ( X C ) . G T . C ) . A N D . ( D A B S ( X C C ) . L T . C ) ) I E R M 1377 I F ( I X X . E O . I ) GOTO 269 1378 I F ( ( O A B S ( X C ) . L T . C ) . A N D . ( D A B S ( X C C ) . G T . C ) ) IDIREC=1 1379 I F ( ( D A B S ( X C ) .GT .C ) . A N D . ( D A B S ( X C C ) . L T . C ) ) IDIREC = 2 1380 IXX=1 1381 269 CONTINUE 1382 GAP=C 1383 I F ( X C C . L T . O . D O ) GAP*-C 1384 DX1=GAP-XC 1385 DX2=XCC-XC 1386 BT0L=DABS(DX1)-DABS(DX2) 1387 I F ( I F L A G . N E . 1 ) GOTO 271 1388 I F U D I R E C . E Q . 2 ) GOTO 273 1389 I F ( ( D A B S ( B T O L ) . L E . T O L ) . A N D . ( D A B S ( X C C ) . G T . C ) ) IFLAG = 2 1390 GOTO 271 1391 273 IF( (DABS(BTOL) . L E . T O L ) . A N D . ( D A B S ( X C C ) . L T . C ) ) IFLAG = 2 1392 271 I F ( I F L A G . E Q . 2 ) GOTO 231 164 1393 IF (I F L A G . E O . 1 ) GOTO 400 1394 I F ( I E R . E Q . I ) GOTO 400 1395 C 1396 C FIND DETAX. DETAV. X. AND V 1397 C 1398 231 DO 200 K=1,NBV 1399 DETAX(K)=DETX(K) 1400 DETAV(K)=ALBE*(DETAX(K)/DETAT-V(K))+DETAT•AC(K)*ALBR 1401 DETAA(K ) = (DETAX(K )/DETAT/DETAT-V(K)/DETAT-0.5D0*AC(K))/BETA 1402 X(K)=X(K)+DETAX(K) 1403 V(K)=V(K)+DETAV(K) 1404 AC(K)=AC(K)+DETAA(K) 1405 200 CONTINUE 1406 DETVL=-BL*DETAX(ID)+DETAX(ID+1) 1407 DETVR=BR*DETAX(ID)+DETAX(ID+1) 1408 ROR=DETAX(ID) 1409 I F ( I F A . E O . O ) DETAXC=DETAX(IC)-A*DETAX(ID) 1410 I F ( I F L . N E . O ) DETAXC = DETAX<IC )-DETAX(IFL )-A*DETAX(ID ) 1411 C 1412 C FIND ABSOLUTE ACCEL. OF EACH FLOOR 1413 C 1414 DO 87 JJ=1,NBV 1415 ABC(Jd)=AC<JJ)+A1 1416 87 CONTINUE 1417 I F O F A . N E . 1 ) ABC( ID ) = AC( ID ) 1418 I F ( I F A . N E . I ) ABC( ID+1 ) =AC(ID+1) 1419 GOTO 400 1420 26 WRITE(5 .60) 1421 WRITE(6 ,60) 1422 60 FORMAT('SOLUTION F A I L E D ) 1423 400 RETURN 1424 END 165 L i s t i n g o f FOUR.S 1 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 2 C C 3 C M A I N P R O G R A M F O U R . S C 4 C C 5 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 6 C 7 C THIS PROGRAM READS IN A DIGITIZED TIME SERIES RECORD 8 C AND COMPUTES THE POWER SPECTRAL DENSITY OF THE 9 C RECORD. 10 C INPUT RECORD X ( I ) IN I/O UNIT 7. 11 C OUTPUT RECORD Y d ) IN I/O UNIT 4. 12 C RESULTS IN I/O UNITS 1 , 2 . 3 . 6 . 8 13 C 14 C DEFINE VARIABLES 15 C 16 C N - NUMBER OF POINTS IN THE RECORD 17 C DT - TIME INCREMENT OF RECORD 18 C N*DT - LENGTH OF RECORD IN SECONDS 19 C 20 IMPLICIT R E A L * 8 ( A - H . O - Z ) 21 INTEGER STAR( 1 )/' * '/ 22 C0MPLEXM6 DATAX ( 4096 ), DATA Y ( 4096 ). CXY ( 4096 ) 23 DIMENSION X(4096 ) ,XSI4096) ,Y< 4096) , Y S ( 4 0 9 6 ) . X Y ( 4 0 9 6 ) . X Y S ( 4 0 9 6 ) 24 DIMENSION C H ( 4 0 9 6 ) . A H ( 4 0 9 6 ) . P H I ( 4 0 9 6 ) . C 0 ( 4 0 9 6 ) . P H I S ( 4 0 9 6 ) 25 C 26 C SET VALUES 27 C 28 PI=3. 1415926535897933 29 C 30 C READ IN DATA 31 C 32 WRITE(5 .87) 33 87 FORMAT('SPECTRUM TYPE . . . O => AUTOSPECTRUM 1 => CROSSSPECTRUM') 34 READ(5 ,STAR)ISPEC 35 READ(4 ,88)INY ,DTY 36 I F ( I S P E C . N E . 1 ) GOTO 90 37 READ(7 ,88) INX.DTX 38 88 F O R M A T ( I 1 0 . F 1 5 . 6 ) 39 I F ( D T X . N E . D T Y ) WRITE(5 .89) 40 89 FORMAT('ERROR :- INPUT/OUTPUT RECORDING INTERVAL ARE DIFFERENT ' ) 41 I F ( ( I N X . N E . I N Y ) . O R . ( D T X . N E . D T Y ) ) STOP 42 90 CONTINUE 43 WRITE(5 ,91) 44 91 FORMAT('INPUT FILE FORMAT . . . O => 8110 1 => 8 F 1 0 . 5 ' ) 45 READ(5 ,STAR)IFLAG 46 IN=INY 47 DT =DTY 48 W R I T E ( 5 , 9 3 ) I N . D T 49 93 FORMAT('# OF POINTS IN THE RECORD = ' . I 5 , 2 X , 50 1 '• A RECORDING INTERVAL (SEC) = ' . F 1 0 . 5 ) 51 WRITE(5 ,94) 52 94 FORMAT('N OF POINTS SAMPLED = 2**IP ' ) 53 WRITE(5 ,95) 54 95 F O R M A T ( ' I P . G E . O AND I P . L E . 1 2 IP=? ' ) 55 READ(5 ,STAR)IP 56 N=DINT(2.D0**IP) 57 WRITE(5 ,96)N,DT 58 96 FORMAT('M OF POINTS SAMPLED = ' . I 5 . 2 X . 166 59 1 '@ A SAMPLING INTERVAL (SEC) = ' , F 1 0 . 5 ) 60 DURA=DT*(2.D0**IP) 61 WRITE(5,97)DURA 62 97 F0RMAT( 'DURATION OF RECORD SAMPLED 1 ' , F 1 2 . 5 ) 63 I F ( I F L A G . E O . O ) GOTO 99 64 R E A D ( 4 , 9 8 ) ( Y ( I ) , I = 1 , N ) 65 I F ( I SPEC.EO. 1 ) READ(7,98 ) ( X ( I ) ,1 = 1,N) 66 98 FORMAT(8F10.5) 67 GOTO 101 68 99 R E A D ( 4 , 1 0 0 ) ( Y ( I ) . 1 = 1 , N ) 69 I F ( I SPEC .EO. 1) R E A D ( 7 , 1 0 0 ) ( X ( I ) , 1 = 1,N) 70 100 F0RMAT(8F10.0) 71 101 WRITE(5 .102) 72 102 FORMATCNFS = 0 => NO FREQUENCY SMOOTHING' ) 73 WRITE(5 .103) 74 103 FORMATCNFS = NFS => SMOOTH 2*NFS+1 ADJACENT SPECTRAL ORDINATES' 75 READ(5,STAR)NFS 76 WRITE(5 ,104) 77 104 FORMAT('NW = 0 => NO DATA TAPERING PRIOR TO FIND PSD ' ) 78 WRITE(5 ,105) 79 105 FORMAT('NW = 1 => DATA TAPERING PRIOR TO FIND PSD ' ) 80 READ(5,STAR)NW 81 WRITE(5 .106) 82 106 FORMAT('SCALE FACTOR = ? ' ) 83 READ(5,107)SCALE1 84 107 F0RMAT(F12.5) 85 AX=O.DO 86 AY=O.DO 87 DO 108 I=1,N 88 Y ( I ) = Y(I )*SCALE 1 89 AY = AY + Y( I ) 90 I F ( I S P E C . E Q . O ) X ( I ) = 0 . D 0 91 I F ( I S P E C . E O . 1) X(I ) =X(I ) *SCALE1 92 • I F ( I S P E C . E Q . 1 ) AX = AX + X ( I ) 93 108 CONTINUE 94 C 95 C TAKE OUT THE MEAN 96 C 97 AX=AX/N 98 AY=AY/N 99 DO 109 I=1.N 100 Y ( I ) = Y ( I ) - A Y 101 X ( I ) = X ( I ) - A X 102 109 CONTINUE 103 CALL PSD(X ,DATAX ,Y ,DATAY ,DT ,N ,NF S ,NW,PI , I S P E C . A H , C H , C O , P H I . X Y ) 104 C 105 C DT - SAMPLING INTERVAL 106 C N - NUMBER OF SAMPLE POINTS 107 C 108 DELF= 1 .DO/(N*DT) 109 N2=N/2 1 10 WRITE(6,STAR)N2 1 1 1 IF ( I S P E C . N E . 1 ) GOTO 111 112 WRITE(1,STAR)N2 113 WRITE(2,STAR)N2 1 14 WRITE(3.STAR)N2 115 WRITE(8,STAR)N2 1 16 WRITE(10,STAR)N2 167 1 17 111 CONTINUE 1 18 DO 112 I=1,N2 1 19 F=(I-1)« DELF 120 WRITE(6. 110) F . Y(I ) 121 I F ( I S P E C . N E . 1 ) GOTO 112 122 WRITE(1. 110) F . C H ( I ) 123 WRITE(2 ,110) F . P H I ( I ) 124 W R I T E O . 1 10) F , AH( I ) 125 WRITE(8. 110) F . C O d ) 126 WRITE(10 .110) F , XY(I ) 127 112 CONTINUE 128 1 10 FORMAT(F 12 . 5 , 1 X . F 3 0 . 2 0 ) 129 C 130 STOP 131 END 132 C 133 C SUBROUTINE PSD COMPUTES THE POWER SPECTAL DENSITY OF 134 C A DISCRETE OUTPUT TIME SERIES Y ( I ) OR COMPUTES THE CROSS 135 C SPECTRAL DESNSITY BETWEEN TWO DICRETE TIME SERIES INPUT X ( I ) 136 C AND OUTPUT TIME SERIES Y ( I ) USING FFT TECHNIQUES. 137 C 138 C Y(I ) - DISCRETE OUTPUT TIME SERIES 139 C - THE OUTPUT PSD IS RETURNED IN Y d ) 140 C X(I ) - DISCRETE INPUT TIME SERIES 141 C - THE INPUT PSD IS RETURNED IN X ( I ) 142 C X Y ( I ) - THE CROSS PSD IS RETURNED IN X Y ( I ) 143 c AH( I ) - THE AUTO-TRANSFER FUNCTION IS RETURNED IN AH(I) 144 c CH( I ) - THE CROSS-TRANSFER FUNCTION IS RETURNED IN C H d ) 145 c PHI(I ) - THE PHASE RELATIONSHIP BETWEEN INPUT AND OUTPUT 146 c - IS RETURNED IN PHI(I ) 147 c C0( I ) - THE COHERENCE BETWEEN INPUT AND OUTPUT IS RETURNED 148 c - IN C0( I ) 149 c 150 c - THE ABOVE RESULTS ARE RETURNED AT FREQUENCIES 151 c - O . D F , 2 * D F . . . .WHERE DF = 1.0/(N*DT ) 152 c 153 c DT - SAMPLING INTERVAL OF THE TIME SERIES 154 c N - 2**P,P INTEGER NUMBER OF SAMPLES OF Y d ) 155 c NFS - 0 NO FREQUENCY SMOOTHING 156 c - NFS SMOOTH 2*NFS+1 ADJACENT SPECTRAL ORDI NATES 157 c SEE NEWLANDS P.145 158 c NW - 0 NO DATA TAPERING PRIOR TO FINDING PSD 159 c - 1 COSINE DATA TAPER APPLIED TO THE TIME SERIES 160 c PRIOR TO FINDING THE PSD 161 c SEE NEWLANDS P. 146-147 162 c DATAX - COMPLEX ARRAY OF X(I ) 163 c DATAY - COMPLEX ARRAY OF Y ( I ) 164 c 165 SUBROUTINE PSD(X .DATAX,Y ,DATAY .DT .N ,NFS .NW.PI , I SPEC.AH.CH.CO,PHI 166 1 ,XY) 167 IMPLICIT R E A L * 8 ( A - H . O - Z ) 168 DIMENSION X ( 4 0 9 6 ) , X S ( 4 0 9 6 ) , Y ( 4 0 9 6 ) . Y S ( 4 0 9 6 ) . X Y ( 4 0 9 6 ) , X Y S ( 4 0 9 6 ) 169 DIMENSION CH( 4096 ) ,AH( 4096) .PHI (4096 ) . CO ( 4096 )'. PHI S( 4096 ) 170 COMPLEX* 16 D A T A X ( 1 ) , O A T A Y ( 1 ) . C X Y I 4 0 9 6 ) 171 PER=DT*N 172 DF=1.DO/PER 173 DO 1 I=1.N 174 DATAY(I ) =DCMPLX(Y(I) .0 .DO) 168 175 I F ( I SPEC.EO. 1 )DATAX(I )=DCMPLX(X(I) ,0 .DO) 176 I F ( I S P E C . E 0 . O ) D A T A X ( I ) = D C M P L X { Y ( I ) , 0 . D O ) 177 1 CONTINUE 178 IF(NW.EQ.O) SCALEX=1.DO 179 IF(NW.EO.O) SCALEY=1.DO 180 IF(NW.EQ.O) GO TO 8 181 C 182 C APPLY COSINE DATA TAPER 183 C 184 XJ=0.9D0*N 185 SUMY=0.DO 186 SUMX=0.D0 187 DO 2 1 = 1,N 188 J=I-1 189 T=J*DT 190 D=1.D0 191 I F ( J . L T . N / 1 0 . D 0 ) D=0.5D0*(1.DO-DCOS(10.DO*PI*T/PER)) 192 I F ( J . G T . X J ) D=0.5D0*(1.DO+DC0S(10.DO*PI* ( T -0 .9D0*PER)/PER)) 193 DATAY(I)=DATAY(I )*D 194 SUMY=SUMY+D**2 195 I F ( I S P E C . N E . 1 ) G 0 T 0 2 196 DATAX(I)=DATAX(I)*D 197 SUMX=SUMX+D**2 198 2 CONTINUE 199 SCALEY=SUMY/N 200 SCALEX=SUMX/N 201 C 202 C CALCULATE FOURIER AMPLITUDES 203 C 204 8 CALL DF0UR2(DATAY.N. 1 , - 1 , 1 ) 205 I F ( I S P E C . E Q . 1 ICALL DF0UR2(DATAX.N, 1 , -1 , 1 ) 206 DO 4 1=1,N 207 DATAY(I)=DATAY(I)/SCALEY 208 I F ( I S P E C . E Q . 1 ) DATAX(I)=DATAX(I)/SCALEX 209 4 CONTINUE 210 C 211 DO 5 1=1,N 212 REY=DREAL(DATAY(I)) 213 AIY =DAIMAG(DATAY(I)) 214 Y(I ) = (REY**2+AIY**2)*2.DO*DT*DT/PER 215 I F ( I S P E C . N E . 1 ) GOTO 5 216 REX=DREAL(DATAX(I)) 217 AIX=DAIMAG(DATAX(I)) 218 X(I ) = (REX**2+AIX**2)*2.D0*DT*DT/PER 219 CXY(I)=(DCONJG(DATAX(I))*DATAY(I))*2.DO*DT*DT/PER 220 REXY=DREAL(CXY(I)) 221 AIXY = DAIMAG(CXY(I)) 222 XY(I)*DSQRT(REXY*»2+AIXY**2) 223 I F ( R E X Y . E Q . O . O ) PHI(I)=O.DO 224 I F ( R E X Y . E O . O . O ) GOTO 5 225 PHI(I)=DATAN(-AIXY/REXY)* 180.DO/PI 226 5 CONTINUE 227 C 228 I F ( N F S . E Q . O ) GO TO 12 229 DO 6 1=1,N 230 Y S U ) - Y ( I ) 231 I F U S P E C . N E . 1 ) GOTO 6 232 X S ( I ) = X ( I ) 233 X Y S ( I ) = X Y ( I ) 234 P H I S ( I ) = P H I ( I ) 235 6 CONTINUE 236 C 237 C APPLY FREQUENCY SMOOTHING 238 C 239 NF2=2*NFS+1 240 DO 7 1 = 1 ,N 241 SUMY=0.DO 242 SUMX=O.DO 243 SUMXY=0.D0 244 SUMPHI=O.DO 245 MSUM=NF2 246 DO 9 J M . N F 2 247 K=I+J-NFS-2 248 I F 1 K . L T . 0 ) K=-K 249 K=K+1 250 I F ( K . G T . N ) MSUM=MSUM-1 251 I F ( K . G T . N ) GO TO 9 252 SUMY=SUMY+YS(K) 253 I F ( I S P E C . N E . 1 ) GOTO 9 254 SUMX=SUMX+XS(K) 255 SUMXY=SUMXY+XYS(K) 256 SUMPHI=SUMPHI+PHIS(K) 257 9 CONTINUE 258 Y(I)=SUMY/MSUM 259 I F ( I S P E C . N E . 1 ) GOTO 7 260 X(I)=SUMX/MSUM 261 XY(I)=SUMXY/MSUM 262 PHI(I)=SUMPHI/MSUM 263 7 CONTINUE 264 C 265 C COMPUTE AUTO GAIN FACTOR SQUARED AH . CROSS GAIN FACTOR 266 C SQUARED CH AND COHERENCE FUNCTION 267 C 268 12 CONTINUE 269 I F ( I S P E C , N E . 1 ) GOTO 10 270 DO 11 1 = 1.N 271 CH(I ) = XY(I )*XY( I )/X(I )/X( I I 272 AH(I )=Y(I )/X(I ) 273 C0( I ) = XY( I )«XY< I )/Y(I )/X( I ) 274 11 CONTINUE 275 10 CONTINUE 276 C 277 RETURN 278 END 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062953/manifest

Comment

Related Items