MECHANICAL MODEL FOR THE ANALYSIS OF LIQUEFACTION OF HORIZONTAL SOIL DEPOSITS by KWOK WING LEE B.So. (Eng.), University of Hong Kong, I966 M.Sc. (Eng.), University of Hong Kong, I968 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of C i v i l Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September, 1975 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of B r i t i s h Columbia, I 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 representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of C i v i l E ngineering The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date Sept. 1975 i ABSTRACT During the development of liquefaction In a s o i l deposit subjected to vibration there are two processes which work in opposite directions. The volume compaction tendency under cyclic loading causes the pore water pressure to rise, and the dissipation of excess pore water pressure (consolidation) decreases i t . Recently, Martin, Finn and Seed(l975) studied the mechanics .of pore water pressure generation of a s o i l sample subjected to cyclic loading and a relation between shear strain cycles, volume compaction and pore water pressure increment was established. A material model based on this relationship is developed in this thesis for saturated granular s o i l under cyclic simple shear conditions. The model includes a hysteretlc stress-strain relationship, volume compaction, pore water pressure rise and dissipation. Using this proposed comprehensive material model, a global mechanical model is constructed to simulate the liquefaction (including consolidation) behavior of a thick horizontal deposit when subjected to horizontal base motion. In this way the coupled problems of dynamic response, pore water pressure rise and consolidation of the deposit under seismic loading can be analysed. The numerical techniques used to solve such problems are discussed in detail. The response of a typical saturated 11 sand deposit under earthquake loading i s determined using the proposed model and the re s u l t s show that the model can predict the various phenomena that saturated sand deposits exhibit during earthquakes. The global model also makes clear the influence of permeability on the li q u e f a c t i o n p o t e n t i a l of the s o i l deposit. i l l TABLE OF CONTENTS Page ABSTRACT 1 TABLE OF CONTENTS i i i LIST OF TABLES v l LIST OF FIGURES v i l NOTATION x i i ACKNOWLEDGEMENTS r v i i CHAPTER 1. INTRODUCTION. 1.1 General Characteristics of S o i l Liquefaction. 1 1.2 Description of Some Liquefaction Case Histories. 2 1.3 Spontaneous Liquefaction, Cyclic Mobility and Liquefaction. 7 1. k Purpose and Scope of Research. 10 CHAPTER 2. REVIEW OF PAST INVESTIGATIONS ON LIQUEFACTION OF SATURATED SANDS. 2.1 Liquefaction Potential of a Deposit Based on Field Experience. 12 2.2 Laboratory Investigation of Liquefaction Using Scaled Sand Models. 18 2.3 Laboratory Investigation of Liquefaction Using Small Sand Samples. 26 2. JJ- Some Methods of Determining Liquefaction Potential of a So i l Deposit. 38 CHAPTER 3. MODELLING OF SATURATED GRANULAR SOIL UNDER CYCLIC SIMPLE SHEAR CONDITIONS. 3.1 Volume Change Ch a r a c t e r i s t i c s During C y c l i c Drained Shear. 3.2 S t r e s s - S t r a i n Relationship Under C y c l i c Simple Shear Conditions. 3.3 Modelling Of Saturated Granular S o i l Under C y c l i c Simple Shear. 3.4 Computation of S o i l Behavior Using the Proposed S o i l Model. CHAPTER 4. MECHANICAL MODEL FOR THE ANALYSIS OF LIQUEFACTION. 4.1 Introduction. 4.2 Formulation of the Problem. 4.3 Solution Scheme. CHAPTER 5. APPLICATION OF THE MECHANICAL MODEL -A DYNAMIC ANALYSIS. 5.1 A Hypothetical S o i l Deposit. 5.2 Response Cha r a c t e r i s t i c s of the S o i l Deposit. 5*3 Magnitude of Pore Pressure Rise and i t s E f f e c t on the Dynamic Response of the Deposit. 5.4 E f f e c t of Pore Pressure D i s s i p a t i o n . CHAPTER 6 SUMMARY AND CONCLUSIONS. 6.1 Summary. 6.2 Conclusions. V 6.3 Suggestions For Further Studies. 147 LIST OF REFERENCES 148 APPENDIX I 155 APPENDIX II 162 APPENDIX III 169 APPENDIX IV 173 v i LIST OF TABLES Table gaffe 5-1 Data for Layer Approximation to Soil Deposit 118 5-2 Frequencies and Phase Angles for the Simulated Earthquake Accelerogram 130 5-3 S'l-Values for Different Layers 130 v i i LIST OF FIGURES Figure Page 1- 1 Modes of Foundation Failure 5 2- 1 C r i t i c a l N-Value versus Depth 15 2-2 Zone of Liquefaction 17 2-3 Probability of Liquefaction 19 2-4 C r i t i c a l Acceleration for Various Sands 21 2-5 Effect of Acceleration Level on Resistance to Liquefaction 2k 2-6 Effect of Surcharge Pressure on Resistance to Liquefaction 25 2-7 Idealised Field Loading. Conditions 28 2-8 Idealised T r i a x i a l Test Conditions 29 2-9 Cyclic Shear Stress Wave Forms 33 2-10 Typical Cyclic Simple Shear Test Data Jk 2-11 Typical Cyclic Tr i a x i a l Test Data 35 2-12 Resistance to Liquefaction in Simple Shear 37 e F Line for Sand B Method of Evaluating Liquefaction Potential Compaction of Dry Coarse Sand Effect of vertical Stress on Compaction Effect of Relative Density and Shear Strain on Compaction Incremental Volumetric Strain Curves Comparison of Measured and Computed Volumetric Strain Hyperbolic Stress-Strain Curve Equivalent Shear Modulus for I n i t i a l Loading Definition of Equivalent Shear Modulus And Damping Ratio Symmetrical Loop Hysteretic Characteristics One-Dimensional S o i l Model Force-Displacement Relationship of Slip Elements Volume Compaction Caused by Cyclic Loading 'Reduced' Incremental Volumetric Strain Increments of Volumetric Strain, €. Hyperbolic Stress-Strain Loops (Drained Test) Showing Effect of Hardening Shear Modulus as a Function of Cyclic Shear Strain and Volumetric Strain Shear Modulus as a Function of Cyclic Shear Strain and Number of Cycles Computed Damping Ratios Shear Modulus as a Function of Shear Strain and Total Compaction Damping Ratios for Monterey Sand One Dimensional Unloading Curves Hyperbolic Stress-Strain Loops (Undrained Test) Showing Softening due to Pore Pressure Rise Strain and Pore Pressure during Undrained Cyclic Shear Test X Figure Page J 4-1 Idealization of the Response Problem 97 4-2 Approximation by Lumped Mass System 100 4-3 Consolidation Model 107 4- 4 Flow Chart 114 5- 1 Properties of a S o i l Deposit 116 5-2 Input Base Acceleration 120 5-3 Surface Acceleration Response 122 5-4 Stress-Strain Response 123 5-5 Surface Acceleration Response Curve 124 5-6 Surface Acceleration from 'SHAKE1 Program 126 5-7 Simulated Earthquake Accelerogram 129 5-8 Surface Acceleration Response (No Pore Pressure) 131 5-9 Strain and Stress Responses for Layer 11 (No Pore Pressure) 133 5-10 Surface Acceleration and Displacement Responses (Pore Pressure Generated) 136 x i Figure Page 5-11 Stress, Strain and Pore Pressure Responses for Layer 11 135 5-12 Stress-Strain Response for Layer 11 (Pore Pressure Generated) 136 5-13 Pore Pressure Distribution at Various Times (No Dissipation) 139 5-14 Pore Pressure Distribution for Different k Values 141 x i i NOTATION A Amplitude of vibration A^, A2. A-j Constants for the shear modulus function a^ t &2 Constants for the increase in shear modulus due to compaction B^, B2, B 3 Constants for the shear modulus function bj_, b2 Constants for the increase in limiting shear stress due to compaction Cit C2, C3, C4 Constants for the volume change function [C ] Damping matrix C eq Equivalent viscous damping factor c± Damping coefficient for the i-th layer c Specific heat of a solid D Damping ratio D m a x Maximum damping ratio D r Relative density D 5 0 T n e maximum particle-size of the smallest f i f t y percent of a s o i l e Void ratio e"F C r i t i c a l void ratio line f i Stress-strain-rate function Stress-strain function Shear modulus Maximum shear modulus I n i t i a l maximum shear modulus Maximum shear modulus at the n-th cycle Thickness of a s o i l deposit Thickness of the I-th layer Coefficient of the earth pressure at rest Bulk modulus of water One-dimensional rebound modulus of the s o i l skeleton Stiffness matrix Stiffness value of the i-th layer Value of permeability Thermal conductivity Mass matrix Mass lumped at the top of the i-th layer Coefficient of volume compressibility C r i t i c a l standard penetration value Number of cycles to liquefaction Equivalent number of constant amplitude shear stress cycles Porosity Inertial force vector In s i t u vertical stress Position vector Shear strength Time Time Volume change function Uniformity coefficient Velocity Energy dissipated per cycle of loading Finite difference approximation to crw Displacement of the i - t h - layer relative to the base Displacement of the i-th layer Base displacement Height measured from the bedrock Height of the top of the i-th layer C r i t i c a l acceleration Constants for the numerical integration Damping coefficients Shear strain Hyperbolic shear strain Unit weight of water Amplitude of a shear strain cycle Increment in shear strain Volume change increment Recoverable component of the volume change increment Compaction compnent of the volume change Increment Strain Volumetric strain Recoverable component of the volumetric strain Compaction component of the volumetric strain Temperature Mass density-Stress Mean normal stress Octahedral normal stress Vertical Stress I n i t i a l v ertical stress Pore water pressure Minor principal stress at failure Effective stress Shear stress Amplitude of a shear stress cycle Cyclic horizontal shear stress Equivalent 10 cycle shear stress amplitude that i s developed Cyclic shear stress amplitude required to xvl cause liquefaction in 10 cycles T j n a x Maximum shear stress that can be applied to a s o i l element Maximum shear stress at the n-th cycle <j>' Effective angle of shearing resistance +1» +2* 3^* 4*4 Constants for the volumetric compaction function x v i l ACKNOWLEDGEMENTS The author wishes to express his gratitude to his principal adviser Dr. W.D. Liam Finn for his invaluable advice and guidance during the entire course of the research. Thanks are also due to Dr. P.M. Byrne, Dr. R.G. Campanella, Dr. G.R. Martin and Dr. Y. Vaid for their encouragement and constructive discussion at various stages of the research. The National Research Council of Canada provided the financial assistance which made this investigation possible. This assistance is gratefully acknowledged. The author deeply appreciated the support and consideration given to him by his wife during the period of his studies. 1 CHAPTER 1 INTRODUCTION 1.1 GENERAL CHARACTERISTICS OF SOIL LIQUEFACTION. It has been observed that when saturated and r e l a t i v e l y loose sandy s o i l deposits are subjected to r e p e t i t i v e or c y c l i c shear loading, such as those Imposed by an earthquake, they often exhibit the phenomena of being turned into a 'quick' or l i q u e f i e d c ondition. In the extreme case where great earthquake shocks occur In an a l l u v i a l region, the pore water i n the l i q u e f i e d s o i l which i s under very large pressure may be forced to the surface forming temporary geysers. In the less severe cases the pore pressure generated i s only a small f r a c t i o n of the o r i g i n a l e f f e c t i v e normal pressure and minor settlement may occur when excess pressure i s d i s s i p a t e d . It i s generally accepted that dry granular s o i l masses compact or contract In volume when subjected to v i b r a t i o n . I f the s o i l i s saturated and drainage prevented then i t i s not d i f f i c u l t to see that the pore water pressure w i l l b u i l d up. For a granular s o i l , the maximum shearing resistance that can be developed on a plane, or shear strength, s, i s given by the following expression, s = (<rv-<rw)tan<t>' ( i . - D 2 where s Is the shear strength, 0*v i s the t o t a l normal pressure, cr w i s the pore water pressure, and, $' i s the e f f e c t i v e angle of i n t e r n a l resistance. As the pore pressure r i s e s the e f f e c t i v e normal str e s s , 0^.-0^, decreases causing a decrease i n the shear strength, s. In the extreme case where the v i b r a t i o n a l disturbance i s of large magnitude and r e l a t i v e l y long duration the pore water pressure may r i s e to such an extent that the e f f e c t i v e normal stress may become zero. The condition of zero or near zero confining pressure makes the s o i l behave l i k e a l i q u i d , since l i t t l e or no shear resistance can be developed according to equation (1-1) and hence the term LIQUEFACTION i s used to describe the phenomenon. There have been instances of sewer tanks and other buried objects f l o a t i n g up to the ground surface. Foundation, slope and dam f a i l u r e s caused by t h i s type of strength reduction i n the supporting s o i l have been reported i n numerous investigations and they are described i n the next section. 1.2 DESCRIPTION OF SOME LIQUEFACTION CASE HISTORIES. In the San Francisco Earthquake of 1906 s o i l l i q u e f a c t i o n phenomena were observed i n numerous a l l u v i a l areas. The quake occurred on A p r i l 18, 1906 with a magnitude of about 8.2. The following i s quoted from the o f f i c i a l report of the earthquake: "One of the most common phenomena i n 3 such a situation was the expulsion of water in jets from apertures which suddenly appeared in the f l a t lying grounds. The water was usually thrown into the a i r for several feet; in some cases i t was reported to be as much as 20 feet, and the ejection continued for several minutes after the earthquake. The continuance of the ejection after the shock indicates that an elast i c stress has been generated in the saturated ground, which thus found r e l i e f In the expulsion of contained water . . . The water in Its passage to the surface brought up considerable quantities of fine sand, which from i t s prevailing light bluish gray colour was evidently derived from considerable depth. On the flood plain of the Salinas river, the sand was recognized by the people of the neighbourhood to be the same as that of a stratum of sand pierced by wells at a depth of 80 feet." Usually the sand brought up in this way formed funnel-shaped sand craters on the ground surface. The abundance of saturated a l l u v i a l f l a t s with such craters gave an Indication of the size of the area which liquefied as well as of the severity of the liquefaction problem. Fortunately, not many structures were built in such areas at the time of the earthquake and l i t t l e or no damage to structures due to s o i l liquefaction was reported. Numerous cases of earth-flows, which may have been due to liquefaction, were also recorded. Liquefaction of a saturated sandy deposit causing damage to buildings and structures occurred in the Nllgata 4 Earthquake of 1964. The quake took place at around 13 hours on June 16, 1964 in the Niigata and Yamagata area in Japan. The magnitude of the quake was less than 7.5 in the Hi enter Scale and the Intensity in the Niigata City was around 8 in the Modified Mercalli Scale. A description and analysis of the effects of the quake have been published in a detailed report compiled by the E d i t o r i a l Committee of "General Report on the Niigata Earthquake of 1964,"(1968). An interesting observation was that modern buildings, bridges and other structures(up to 1964) built on the a l l u v i a l formation around the Shinano river were damaged more readily than others. Osaki(1968) reckoned that two thirds of the damaged buildings In these areas Just settled, t i l t e d or overturned as a whole without much damage in the superstructure, indicating severe foundation f a i l u r e . Moreover, for this type of failure the usual upheaval of s o l i on the opposite side of t i l t i n g or overturning predicted by conventional foundation failure theory, as shown in figure 1-la, was not observed but a deformation pattern similar to figure 1-lb was obtained. In contrast to the settlement and t i l t i n g of heavy superstructures as described, light buried structures, such as concrete septic tanks, 'floated' during the quake and remained protruding above ground surface after the disturbance ceased. In the light of these data, Osaki reasoned that the s o i l at such sites behaved like a liquid and he concluded that the foundation s o i l liquefied under the earthquake vibration. Also an area near the same location, but on which no (a) ( b ) Fig. 1-1 Modes of Foundation Failure (after Osaki, 1968 ) 6 superstructures were built, mud volcanoes and temporary geysers erupted during the earthquake and these outbursts, which continued for a while ' after the tremor passed away, fi n a l l y died down giving rise to sand craters. Not only was loose sandy deposit susceptible to liquefaction under earthquake vibrations, medium dense sand could also liquefy when the Intensity and duration of the quake were both large. Seed(1972) reported that in the San Fernando Earthquake of February 9, 1971, the liquefaction of a medium dense sand foundation beneath the Jensen F i l t r a t i o n Plant induced large la t e r a l movements in the f i l l above the sand. As a result of this, extensive damage was done to the site as well as to structures which were under construction. Failures caused by s o i l liquefaction are not restricted to foundations of structures. Numerous earth-flows and slope failures have occurred during earthquakes due to liquefaction of granular s o i l masses or seams or lenses of similar material. Seed(1968) in his Terzaghi Lecture gave a large number of examples of slope failures or slides resulting from s o i l liquefaction and he examined the mechanism of failure of each case. Using a recently developed method of dynamic analysis, Seed et a l (1975) reconstructed the mechanics of sliding of the upstream part of the embankment of the Lower Fernando Dam, which moved some 70 feet or more into the reservoir. They found that an extensive zone of liquefaction developed near the base of the embankment and 7 their conclusion was confirmed by the observed phenomena of sand bolls, cracks f i l l e d with liquefied sand in the slide mass and spreading of the embankment s o i l 250 feet beyond the toe of the embankment. 1.3 SPONTANEOUS LIQUEFACTION, CYCLIC MOBILITY AND LIQUEFACTION. With regards to the term liquefaction there i s a controversy on i t s usage. There has been no definition for i t in terms of stress or strain nor has there been any in terms of other phenomenological parameters. Castro(l969) in his Ph.D. thesis pointed out that Allen Hazen(1920) was the f i r s t one to use the word "liquefies" to describe flow failures and he continued to use the term in association with flow slides in his thesis. He wrote : "Liquefaction or flow failure of sand Is caused by substantial reduction of i t s shear strength. The mass of sand actually flows spreading out u n t i l the shear stresses acting within the sand become so small that they are compatible with the reduced shear strength of the sand". However, Terzaghi and Peck(l968) used the term spontaneous llquefaction to describe this phenomenon of flow failure which occurred mostly in loose fine sand when subjected to mild shock. They suggested that the metastable structure formed by the sand grains was responsible for this kind of behaviour. It was because of such a metastable structure of the sand 8 aggregate that the failure process could take place in a very short time. Indeed, the term spontaneous liquefaction i s very appropriate in describing this type of flow failures. Several researchers in the f i e l d of s o i l dynamics, e.g.. Seed et a l at the University of California and Finn et a l at the University of British Columbia, carried out cyclic loading undrained tests on t r i a x i a l as well as simple shear specimens of granular material and observed In general that the pore pressure in the specimens continued to rise, leading to a decrease in the effective confining (or normal) pressure, and f i n a l l y reached a stage where large deformations occurred under constant cyclic stress amplitude. The term "liquefaction" was also used in connection with this phenomenon. This is quite suitable since the shear resistance of the s o i l decreased with increasing number of cycles of stress u n t i l i t f i n a l l y failed when i t offered l i t t l e or no resistance at a l l to large deformation. However, depending on the denseness of the sand specimen there are varying degree of dilative response at the peak shear so that for example very large strains cannot develop abruptly in f a i r l y dense sand, although for very loose and loose samples an additional cycle of shear stress near liquefaction can increase the maximum shear strain by ten to fifteen per cent. In view of the above described dilative response of some samples, Castro(1969)• following the suggestion of Dr. A. Casagrande, referred to the process of 9 the gradual development of pore pressure with cyclic shear straining as cyclic mobility and reserved the term liquefaction for the failure of s o i l under shear when the loss of strength due to rise in pore pressure was not regained (due to contractive response) even at large strains. It i s well-known that vibration causes compaction in a granular s o i l mass. In other words, It can be said -that cyclic stresses, especially shear stresses, acting on a granular s o i l mass always cause an increase in volumetric strain (contraction). Each additional cycle of stress or strain always causes an additional Increment in contraction though for constant amplitude stress or strain the increment is a decreasing function of the number of cycles. If the s o i l mass i s saturated and drainage (volume change) prevented, a. portion of the confining (or normal) pressure acting on the s o i l skeleton would be transferred to the pore water, causing an incremental rise in pore pressure. Since the volume of s o i l always decreases when drained, as described above, the pore pressure w i l l always rise in the undrained case. As a result there is always the possibility of making the effective confining pressure zero by continuing the process of cycling the stresses or strains. Though the term cyclic mobility describes the cyclic nature of the process quite adequately the word mobility does not give the notion of the gradual softening of the material like the word liquefaction does. It may be more suitable to use the term cyclic liquefaction J 10 Terzaghl and Peck(1968) In the recent edition of their book, "Soil Mechanics In Engineering Practice", added a new section on 'Liquefaction under Reversals of Stress or Strain* to describe the behavior of the rise in pore pressure in a s o i l mass due to cyclic loading which may lead to f a i l u r e . Thus the broad usage of the word 'liquefaction' to Include flow failures, with sharp rise in pore pressure and catastrophically large strains (spontaneous liquefaction), as well as limited failure, with incremental or continual rise in pore pressure and comparatively smaller strains (liquefaction under reversals of stress or strain), is not an over-generalization. As a matter of fact, the word liquefaction w i l l be used in this broad sense throughout this thesis. More discussion of the similarities of these behaviors w i l l follow in the coming chapters. 1.4 PURPOSE AND SCOPE OF RESEARCH In the liquefaction process of a s o i l deposit subjected to vibration there are two processes which work in opposite directions. The volume compaction tendency under cyclic loading causes the pore water pressure to rise, as described above, and the dissipation of excess pore water pressure (consolidation) decreases i t . As yet, no comprehensive material model has been presented to explain this coupled problem of liquefaction. Recently, Martin, Finn and Seed(1975) studied the mechanics of pore water pressure 11 generation of a s o i l sample subjected to cyclic loading and a relation between shear strain cycles, volume compaction and pore water pressure increment was established. In this research, a material model for granular s o i l under cyclic simple shear conditions has been constructed based on this relationship that includes a hysteretlc stress-strain relation, volume compaction, pore water pressure rise and dissipation. Using this proposed comprehensive material model, a global mechanical model is constructed to simulate the liquefaction (including consolidation) behavior of a thick horizontal deposit when subjected to base motion. Liquefaction may then be investigated as the cummulatlve consequences of the following interacting processes: 1. Dynamic response to base motion, 2. Increase in volumetric compaction strain due to cyclic straining, 3. Pore water pressure generation as a result of (2), and, 4. Dissipation of excess pore water pressure (consolidation). Numerical techniques for the solution of the liquefaction problem are developed and an example problem is solved showing that the model is capable of simulating the actual behaviour of a s o i l deposit as well as giving the influence of permeability or coefficient of consolidation on liquefaction potential. 12 CHAPTER 2 . REVIEW OF PAST INVESTIGATIONS ON LIQUEFACTION OF SATURATED SANDS 2.1 LIQUEFACTION POTENTIAL OF A DEPOSIT BASED ON FIELD EXPERIENCE. As early as 1906, after the San Francisco Earthquake, i t was observed that the places that were most susceptible to s o i l liquefaction were the a l l u v i a l f l a t s where loose to very loose saturated sand aggregates were the main constituent of the s o i l deposits. It was recognised at that time that these loose sands were responsible for the liquefaction of the s o i l deposits as well as for many earthflows. However, apart from these qualitative descriptions of the soils that were most susceptible to liquefaction, the values of the s o i l parameters, i f ever known at that time, were not reported. Florin and Ivanov(196l) published some of their f i e l d investigations into the s o i l liquefaction problem. They studied the liquefaction behavior of relatively loose sandy strata which were around eight to ten metres (around 30 feet) in depth using explosives as an energy source for vibration. Throughout the investigation they used a charge consisting of 5 kilograms of ammonite buried at a depth of 4.5 13 meters (about 15 feet). They observed that after a blast, the loose saturated sandy s o i l deposit liquefied with large surface settlements after re-consolidation whereas in dense sands no evidence of liquefaction could be found and l i t t l e or no surface settlement was observed. Also, in the case of a loose sandy deposit liquefaction always occurred within a few seconds of f i r i n g the explosives. As a result of these experiments they were able to establish a criterion for s o i l liquefaction using the magnitude of settlement due to blasting as a key variable, namely:-1. If the average settlements (due to the blast created by the explosion of the 5 kg of ammonite described above) in a radius of 5 meters (16 feet) was less than 8 to 10 cm (around 3.5 inches), no provision against liquefaction should be made, and, 2. If the ratio of settlements between successive blasts was greater than 1:0.6 liquefaction of the stratum would be very l i k e l y to occur. It was also reported that this method of standard blasting provided quite reliable results. However, i t suffers the following disadvantages:-1. Blasting has to be done on the particular site under investigation which may be quite objectionable in some cases, 2. Uncertainty involved when applying the criterion to soils with properties d i f f e r from those described in the 14 literature, and, 3. Magnitude and duration of vibrational disturbance cannot be taken into account. Koizumi(I966) examined the boring data for the s o i l conditions in the Niigata area before and after the 1964 Niigata earthquake in an effort to study the changes in s o i l density due to the earthquake. From standard penetration N-values for about 20 locations in the Niigata area he noticed that some loose sand strata were compacted by the earthquake while other not-so-loose sand strata were loosened. Based on these observations he reasoned that for the sandy s o i l at the Niigata area a c r i t i c a l void ratio as a function of overburden pressure should exist. A curve was proposed to relate c r i t i c a l standard penetration values, N c r, which were used as a measure for the c r i t i c a l void ratio at various depths, with depth of overburden as shown in Figure 2-1. Thus, according to Koizumi, a sand deposit having the same N values as given by N__ would not suffer any change in density when subjected to an earthquake similar to that of Niigata. However, for deposits with N-values above or below the N c r values, loosening or compaction would occur i f the sand deposit were subjected to the same quake again. Whenever compaction or volume contraction occurred in the saturated sandy deposit during an earthquake pore water pressure would rise causing the effective stress in the s o i l to decrease. Liquefaction would be the result when the effective stress approached a N-Value 0 10 20 30 UO 50 • - N c r ig.2-1 Critical N-Value versus Depth (after Koizumi, 1966) 16 zero or near zero value. Consequently, the Ncr curve could also be viewed as a c r i t e r i o n f o r s o i l l i q u e f a c t i o n . Using t h i s c r i t e r i o n , Osaki(1968) ' determined the zone of l i q u e f a c t i o n of" sand deposits i n his study i n order to re l a t e the depth and thickness of l i q u e f i e d s o i l s t r a t a to the amount of damage to super-structures. The act u a l N-value curve f o r a deposit and the Ncr curve were plotted i n the same graph, and the zone of l i q u e f a c t i o n was obtained from the points of in t e r s e c t i o n , a and b, as shown i n Figure 2-2. Kishida(1964,1969,l970) published several papers on the c h a r a c t e r i s t i c s of l i q u e f i e d sands i n an attempt to a r r i v e at a general c r i t e r i o n f o r accessing the l i q u e f a c t i o n p o t e n t i a l of a l e v e l deposit and he suggested that the following items determined the occurrence of l i q u e f a c t i o n of a deposit :-1. Earthquake i n t e n s i t y around V to VI according to the Japanese Meteorological Agency Intensity Scale. 2. E f f e c t i v e overburden pressure less than 2 kg/sq.cm. 3. Relative density less than 15%, 4. Saturated. 5. Coarse-grained s o i l with 0.074<D^ Q<2.0 mm. 6. S o i l grains are f a i r l y uniform with U.<10. From item number 3 a curve of N-values versus depth was obtained. The curve, which has a s i m i l a r usage as the Ncr curve, i s shown i n Figure 2-3. In the same figure the Nrty. curve which was obtained by Koizumi i s also plotted f o r N-Value 17 Fig.2-2 Zone of Liquefaction (after Osaki, 1968) 18 comparison. It i s quite obvious that these curves are applicable to the Niigata deposit subjected to the earthquake of 1964. Use of these curves to predict the liquefaction potential of deposit located elsewhere would be quite dubious since the dynamic stresses developed depend on many parameters relating to the s o i l properties, as well as those relating to the base rock motion. Nevertheless, these c r i t e r i a are quite useful in the preliminary stage of site investigation so that a direction for detailed site investigation may be decided. 2.2 LABORATORY INVESTIGATION OF LIQUEFACTION USING SCALED SAND MODELS. As early as 1936 Casagrande(l936) constructed model dam sections to demonstrate the s t a b i l i t y of sandy slopes. Two model dam sections, one loose sand and one of dense, were built of ordinary beach sand and placed in tanks mounted on casters. Water level on the upstream side of the dam model was kept at a small distance below the crest while on the downstream side the base was covered with a pervious f i l t e r blanket to keep the line of seepage away from the face of the slope. The tanks were shaken and i t was observed that the loose sand model liquefied (failure of flow slide type) very readily while the dense model remained intact. However, no directly useful numerical data were obtained from the model testing. N-Value 0 10 2.0 30 U0 High -^ \/ /\ N Y / / \ \ y / / \ \> / s \ \ N. / / . X / v / — - Low Possib ility Possibil Relative Density ty «—s\ ^ " \ / \ / ^ \ N ' X ^ \ P / X s dium D s s i b i l i t y ' Fig. 2-3 Probability of Liquefaction (after Kishida, 1969) 20 Maslov(1957) and Florin and Ivanov(l96l) performed model testing, as well as f i e l d testing, for the liquefaction of saturated sandy s o i l masses. Complete details of their models and experimental procedures have not been published so that accessment of boundary effects cannot be made. However, their efforts seemed to have been focused on the distribution and dissipation of excess pore water pressure that had been already generated by, seismic (vibrational) disturbances. Maslov suggested the existence of a c r i t i c a l acceleration, °crit» w n o s e value depended on the s o i l density and type, and that i f the maximum acceleration of the vibrational disturbance was less than the required for the particular s o i l mass no excess pore water pressure would be generated. A graph of c * c r i t for various soils is shown in Figure 2-4. The work of Florin and Ivanov was mainly concerned with the consolidation or (re-consolidation) process of a liquefied s o i l stratum, though, as discussed in the previous section, they established a criterion for liquefaction based on the amount of settlement due to blasting at the s i t e . They also made a comparison between the distribution of excess pore water pressure obtained theoretically and that obtained through model testing and seemingly good agreement was arrived at. However, neither of these investigators attempted to access the amount of excess water pressure generated as a function of the s o i l parameters and the acceleration level and duration of the vibrational disturbance. 21 POROSITY (Vo) Fig. 2-4 Critical Acceleration for Various Sands. (After Masloy. 1957) 22 Yoshlml(1967) placed loose saturated sand in a r i g i d box, 50cm x 25cm x 25cm in dimension, and applied horizontal vibration to the box in an effort to simulate the seismic effect on sand deposits. An impervious flexible membrane was put over the sand in the box so that a surcharge could be applied by means of a i r pressure on the membrane. Using such a set up, he observed that the liquefaction process consisted of two stages, an i n i t i a l compaction stage in which pore pressure rose in a near-hyperbolic fashion and a sudden liquefaction stage when failure occur in a short time. Although quantitative results were not directly applicable to f i e l d problems because of the difference in boundary conditions, qualitative extrapolations could be made to describe the behavior of a sand stratum located at depth and sandwiched between impervious layers. Recently, Finn et al(1970,1971) reported tests performed on saturated sands using relatively large sand samples. The size of sample used was approximately 72 inches by 18 Inches by 7 inches and the direction of base motion input was along the longest side of the sample so as to minimize, i f not eliminate, the boundary effects Imposed on the sample by the container. An electro-mechanical servo-controlled shake table was used as a source for imparting base motion to the sample. The advantage of using such a shake table is that various controlled input formats can be used. The input can be either acceleration, velocity or displacement 23 controlled and the wave form of the input can be either sinusoidal, triangular or square or irregular earthquake records. Test equipments and procedures were developed so that sample uniformity was within reasonable tolerance of +2% variation in void ratio, see e.g., Emery, Finn and Lee(1973). For sinusoidal base input the ratio of the response acceleration amplitude to that of the base input at various points within the sample was also measured in several samples and was shown to be close to unity when the frequency of input was less than 10 hz. Numerical results were obtained for constant amplitude sinusoidal base motion input relating the amplitude of acceleration to the number of cycles of vibration to cause liquefaction, N T, as shown in Figure 2-5. Effect of i n i t i a l effective normal (or surcharge) pressure, o-', on N V J j was also examined and the result is as shown in Figure 2-6. The general conclusion that can be obtained from these model tests can be summarized as follows :-1. High degree of saturation of the sand model is required. 2. The looser the sand the shorter the time (or smaller number of cycles) to liquefaction. 3. The higher the i n i t i a l normal pressure (or I n i t i a l effective confining stress) the longer the time (or larger number of cycles) to liquefaction. 4. The larger the acceleration input the shorter the time (or smaller number of cycles) to liquefaction. It should be pointed out that these conclusions agree 0.5 0.U h o I z o I— < DC LU _ l LU O O < 0.3 0.2 0.1 0.0 10 N, 100 1000 . I 10000 CYCLES TO LIQUEFACTION Fig. 2-5. Effect of Acceleration Level on Resistance to Liquefaction, (after Finn et al, 1971) 8.0 ~ 6.0 Q . LU O < o or 3 CO U.0 2D 10 N L 10000 CYCLES TO LIQUEFACTION Fig. 2-6. Effect of Surcharge Pressure on Resistance to Liquefaction. (After Finn et al, 1971) ro 26 extremely well with f i e l d observation. In addition to this, model testing provides a basis for evaluating the effect of a single variable on the liquefaction potential while from f i e l d observations the effect of a single variable i s seldom presented in a clear-cut manner. In passing, mention should be made that in most, i f not a l l , model testing the results obtained are useful as a stepping stone for predicting f i e l d behavior. If any theory of material (or soil) behavior i s proposed to describe or predict the behavior of an actual structure (including s o i l deposit), the theory should be able to predict the behaviour of the model with better accuracy since, 1. the loading on the model can be controlled and measured, 2. the boundary conditions of the model are known, and, 3. the material in the model can be made more homogeneous and isotropic than any f i e l d deposit i s l i k e l y to be. In order to arrive at some 'theory* for the liquefaction behavior of sandy material, testing of small samples under presumably uniform stress or strain conditions should not be overlooked. These are described in the next section. 2.3 LABORATORY INVESTIGATION OF LIQUEFACTION USING SMALL SAND SAMPLES During an earthquake a typical s o i l element in a 27 deposit Is deformed by a system of dynamic stresses as a result of earthquake waves passing through the s o i l medium. The most significant set of dynamic stresses are considered generally to be'caused by the propagation of shear waves from the bedrock to the ground surface. For a horizontal s o i l deposit on a horizontal bedrock formation this dynamic stress system can be idealized as shown in Figure 2-7. I n i t i a l l y , the s o i l element is subjected to the 'at rest* principal stresses, cr^0 and K0<r^0, where cr^Q i s the i n i t i a l effective normal stress and K Q i s the coefficient of earth pressure at rest. During the earthquake the dynamic simple shear stress system as shown in Figure 2-7b is imposed resulting in the stress system as shown in Figure 2-7c. Test procedures have been developed to apply a uniform stress system similar or equivalent to that of Figure 2-7c and three of the most useful and commonly employed ones are :-1. Cyclic Loading T r i a x i a l Test - This is quite similar to the conventional t r i a x i a l test though a cycling deviator load i s used instead of the unl-directional deviator load. The cylindrical saturated s o i l specimen is f i r s t consolidated under an ambient pressure of <T% then, with o further drainage from the sample prevented, a cycling axial load is applied so that the axial stress varies between <r0+/^a:D a n d °o~ A 0D P e r i o d i c a l l y • F o r f u l l y saturated s o i l specimens the effective stress applied can be represented as shown in Figure 2-8a and b. It should be noted that the actual f i e l d conditions are not a. b. c. Fig. 2-7. Idealised Field Loading Conditions. CO 29 Fig. 2 - 8 I dea l i sed T r i a x i a l Test Cond i t i o n s . 30 exactly duplicated since the actual principal stress system rotates approximately +40° periodically while that of the laboratory test conditions fluctuates by +90°. The i n i t i a l stress systems are also different; the s o i l element in the f i e l d is K -consolidated while the sample o in the t r i a x i a l apparatus is isotropically consolidated. 2. Cyclic Simple Shear Test - In this test a s o i l sample of rectangular cross-section i s confined in a 'shear box' with r i g i d sides. It i s then consolidated one-dimensionally under a vertical stress of <T^ 0. During testing drainage i s prevented and a cycling shear stress is applied to the top and bottom of the sample. Moreover, two sides of the 'shear box' rotates in such a way so that the uniform stress or strain condition as shown in Figure 2-7c i s duplicated closely. Since the stress conditions in the simple shear test closely resembles those acting on a s o i l element in the f i e l d test data obtained can be applied directly to f i e l d conditions without any correction. 3. Cyclic Torsion Test - The set up of equipment in this test is similar to that of cyclic t r i a x i a l test. But the cyclic shear stresses are applied through the use of torsional load applied to both ends of the cylindrical sample while the axial load is kept constant. Ishihara and Li(1972) used this method of applying shear deformation to the samples and they designed a t r i a x i a l torsion shear test apparatus that can be used to run 31 tests with constant confining pressure as well as tests with no lateral strain. It was claimed that the f i r s t type of tests corresponds tc the cyclic t r i a x i a l test while the * second type corresponds to the simple shear test. However, more correlation of this test with the other two is needed before the test data can be applied to f i e l d conditions since the usage of average stress and strain to represent the varying stress and strain fields Inside the cylindrical specimen under torsion needs more substantiation. Recently Ishibashi and Sheriff 1974) solved the problem of varying stress and strain f i e l d Inside a solid cylindrical specimen under torsion by using an ingenious sample cross-section. They used a 'donut-like 1 sample in order to achieve a uniform shear strain f i e l d and consequently they were able to compare their test results with those obtained from the t r i a x i a l and the simple shear apparatus based on the actual stress and strain values. Investigations of sand liquefaction using the above mentioned types of laboratory test have been made by Seed and Lee(1966), Peacock and Seed(1968), Finn et al(1971), Ishlhara and Ll(1972) and Ishibashi and Sherif(1974). Most of the tests performed involved the use of a constant amplitude cyclic shear stress, x^t as shown in Figure 2-9. Usually during the test the pore water pressure, cr , cyclic strain, 7 t and the cyclic stress, T , were continuously monitored by 32 using e l e c t r i c a l transducers and chart recroders. Typical records of simple shear and t r i a x i a l tests are as shown in Figures 2-10 and 2-11. It can be seen that during the i n i t i a l portion of the test where the pore pressure rises gradually, the amplitude of cyclic strain developed i s almost zero (of the order decimal of a percent) but becomes increasingly large magnitude (+15 to 20^ ) towards liquefaction failure when the pore pressure rises more rapidly. It was found that most sands behaved in this way and that their resistance to liquefaction could be described in a simplified manner by using the following variables :-N L: number of stress cycles to cause liquefaction. D rt relative density, defined as :-D r = ** a* " 6 x 100% (2-1) max min where e i s the current void ratio of the sand sample and emax a n d emin a r e t n e v o i d ratios at the loosest and densest state respectively. This variable is generally used to indicate the state of denseness of the sand specimen. t d ! cy°ll c s n e a r stress amplitude, and, C"yQ: the i n i t i a l effective ve r t i c a l (or confining) stress. It was found that the cyclic stress amplitude, T d» and the i n i t i a l effective vertical stress, (j-' . can be vo combined and (TA/crJrs) can be considered as a single 33 Fig. 2 - 9 Cyclic Shear Stress Wave Forms TEST 58 e = 0.617 cc UJ — < e w UJ * cc ~" 2 => 2.0n i.oH UJ ce 3 Z (O < UJ X cc o-1 ! PORE PRESSURE iff rfflllhl nflWi'rr SHEAR STRAIN «f nr f. r n f. r=Ml Hi h n I I i i111 II I;II :.| i l l Mi l l : i i l i I i . l .uhll • i.i'i t.111y11U11' ¥jr»r r rr i r i r i M i i i i i i i i i i i i i i i i i • i i i i i i i i i i 1 i i i i i i i I i i i i i 11 1 I M i l ni r 1 5 -10 0.4-cc e < o £ 2> 0.2-05 — 0 ^ X o p z < P I 0.2 £ S 0.4 b £ < c/) n fi i i II ii ! II II II n i SHEAR STRESS II M II M » M i M II II III II II l | II III II II I II III II II ii ii ii Ii ii ii ii ii iii ii ii H II II II II II II II n in II it II i II ii i II H ii ii H A A A A II n II II n n II II ii :i ii i! H II n ni .. Ml II III II A A ft II III II II I! Ii II II II II fl II II II II II II il II II v y v II n II II II n ii II II H II hi v H y i H n II ii II II II H i i i i II II II H n II n n ii i ,i M ii II H II ii II H II II n u u u u v u y n n J U U I it I I A: •10 L-15 x I z < — cc sS cc < UJ X w I 11 y v u ., II H H i H II ii V v y SHEAR FRICTION =0.03 Kg/cm' . 2-10 Typical Cyclic Simple Shear Test Data (after Finn et al, 1971) 36 Independent variable. Finn et al(197l) used a summary plot to present the test results relating N_ , D and (T./(J' ) in a L r d vo more efficient way as shown in' Figure 2-12, where each curve represents a constant (T,/o"f ) value. It can be seen from o. vo this figure that the resistance to liquefaction, as measured by N T, decreases as :-1. the relative density of the sample becomes smaller, 2. the cyclic shear amplitude becomes larger, and, 3 . the i n i t i a l effective vertical (or confining) stress becomes smaller. Although the cyclic t r i a x i a l test imposes a stress condition on the sand sample different from that in the f i e l d , Finn et al(l971) found that i f due consideration is given to the i n i t i a l effective normal pressure acting on the 45°-plane in the t r i a x i a l test specimen correlation between the cyclic simple shear test and t r i a x i a l test is possible. In addition, the recent development of the novel Torsional Simple Shear Device by Ishibashi and Sheriff 1974) makes i t possible to run liquefaction tests with varying degree of lateral confinement. Moreover, they found that the plot of AT^/<7^ct, the ratio of maximum change in the shear stress to the octahedral normal stress, versus the number of cycles to liquefaction bears close resemblance to that obtained by Finn et al(1971). This correlation has a non-trivial implication since i t indicates that the liquefaction of sand produced in the laboratory is not an a r t i f i c i a l , or apparatus-dependent. 37 (07')0 = 2.0 Kg/cm1 X- w Marked for Edch Line .75 a> I o < o > .70 .65 .60 .55 .50 .21 VI5 VI0 \ .34 V \ \ y.50 \ 10 100 1000 10000 CYCLES TO LIQUEFACTION - T 1 , Fig. 2-12 Resistance to Liquefaction in Simple Shear (After Finn et al, 1971) 38 phenomenon and the data obtained is good for application to the analysis of actual s o i l structures. The method of adapting liquefaction data obtained In the laboratory to f i e l d problem is discussed in the following section. 2.4 SOME METHODS OF DETERMINING LIQUEFACTION POTENTIAL OF A SOIL DEPOSIT. In sections 2.1 and 2.2 of this chapter, some of the methods of determining liquefaction potential were discussed, for example :-1. Florin and Ivanov(l96l) - using settlement due to blasting at the site as a key variable, 2. Osaki(l968) - using the c r i t i c a l standard penetration values, N c r, as a criterion, 3. Kishlda(1969) - using the 75% relative density as a criterion for the occurrence of liquefaction, and, 4. Maslov(l957) - using the concept of c r i t i c a l acceleration level, c^Q-pif to determine the probability of having a liquefaction, and they are not repeated here. The following is mainly concerned with the application of laboratory test data to fi e l d problems. Castro(I969) following the suggestion of Prof. A. Casagrande used the concept of c r i t i c a l void ratio to determine the liquefaction potential of a sandy s o i l . A series of stress-controlled t r i a x i a l undrained tests, which .3 39 were termed R tests, were performed on representative samples from the f i e l d for various void ratio or relative density values and at different i n i t i a l consolidation pressures, cr^. For samples that developed continuous pore pressure rise at constant axial load, (SPONTANEOUS LIQUEFACTION), the values of the effective minor principal stress, ^ f * w e r e plotted against the void ratio of the sample, e^, A curve was drawn through the points and a e~p line was obtained as shown In Figure 2-13. For a s o i l element with consolidation pressure a^ a and void ratio e^, represented by point A in Figure 2-13, would have a higher liquefaction potential than an element of the same s o i l with <r^ and e^, represented by point B, since the horizontal distance AC and BC indicate the amount of reduction in effective minor principal stress, or pore pressure r i s e , at liquefaction. A s o i l element with (r*^ and e^f represented by point D, would not develop continuous pore pressure rise at constant loading and is said to be d i l a t i v e . However, in this method of determining liquefaction potential there are a number of factors not being taken into account and two of the important ones are :-1. Magnitude and duration of the disturbance that cause liquefaction, and, 2. Effect of a further stress reversal on the sample that exhibit dilative behavior. Nevertheless, the method is a good guideline for recognizing deposits with metastable sand aggregates which would develop spontaneous liquefaction. 0*7 41 Seed and Idriss(1967.1971) proposed a method f o r determining the l i q u e f a c t i o n p o t e n t i a l of s o i l deposit subjected to earthquake forces. The time h i s t o r y response of ho r i z o n t a l shear stresses, T , at various depths was xy evaluated using dynamic a n a l y s i s . Based on the i r r e g u l a r time h i s t o r y stress response the amplitude of an equivalent 10 cycles of constant c y c l i c shear stress was then determined at each depth so that a plot of equivalent 10-cycle shear stress amplitude, T ^ Q , versus depth was obtained as shown i n Figure 2-14, curve AB. From the laboratory c y c l i c simple shear test data, or c y c l i c t r i a x i a l test data m u l t i p l i e d by a c o r r e l a t i o n f a c t o r , the shear stress amplitude that caused l i q u e f a c t i o n i n 10 cycles, T ^ Q » was determined at each depth using the p a r t i c u l a r values of D r and cr^ at the p a r t i c u l a r depth and the curve CD as shown i n Figure 2-14 was obtained. The zone of l i q u e f a c t i o n , i f there i s one, i s determined by the depth i n t e r v a l within which Td| 0 > Tno» e»S«» z a to z^ i n Figure 2-14. To sim p l i f y the procedure of determining l i q u e f a c t i o n p o t e n t i a l Seed and Idriss(197l) proposed that instead of approximating the stress response by an equivalent 10 cycles of constant amplitude c y c l i c shear s t r e s s , T,,rt, a d i r e c t dlO approximation was used to determine the equivalent cycles , N , of constant c y c l i c shear of amplitude, x . However, eq av the basic method of determining l i q u e f a c t i o n p o t e n t i a l was not changed. At f i r s t glance, i t would seem that conclusions 42 Stress Depth FIG.Z-14 Method of Evaluating Liquefaction Potential (After Seed et al, 1971) 43 about liquefaction potential based on the c r i t i c a l void ratio concept contradicts that using cyclic test data. According to the c r i t i c a l void ratio concept an increase in the i n i t i a l consolidation pressure would increase the liquefaction potential while according to the latter the potential is reduced. However, closer examination reveals that the term 'liquefaction potential' was used to convey different meanings in the two cases. In actual fact when the i n i t i a l consolidation pressure is increased, a. the c r i t i c a l void ratio concept Indicates a higher potential of failure in the spontaneous liquefaction mode but higher stresses are needed to i n i t i a t e such a failure, while, b. the cyclic test data indicates that for the same magnitude of disturbance more cycles of stresses are needed to cause the pore pressure to rise to the failure level (smaller liquefaction potential) though that the dilative response which has a stabilizing effect may be reduced is not mentioned. It can be seen that while the c r i t i c a l void ratio method emphasizes the classification of failure modes, dilative or contractive behaviour at failure, the cyclic stress method is mainly concerned with the magnitude and duration of the disturbing stresses to cause fai l u r e . In fact, the combination of the two methods would give the complete picture of liquefaction failure of s o i l structures when subjected to earthquake loading. Even though the method of determining liquefaction potential using cyclic test data is able to provide some approximate solutions which are useful in the design of actual of actual structures, there are s t i l l some questions l e f t unanswered by this method, e.g.:-1. pvnamlo Response - In arriving at values for N and T <' a dynamic analysis using equivalent shear moduli and damping coefficients are used. Seed and Idrlss(l969)• The non-linearity in s o i l stress-strain is approximated but the change of shear modulus due to changes in effective normal stress and compaction (hardening) is not taken into account, 2. Time to Liquefaction - Since the time history of pore pressure response cannot be determined by the method the time to liquefaction cannot be determined, 3. Maximum Strain at Liquefaction - As a result of approximating the actual stress-strain by a straight line the value of maximum strain obtained is not a r e a l i s t i c one and to obtain such strain value the actual stress-strain should be used, and, 4. The Effect of Consolidation during pore pressure rise on the liquefaction potential of a s o i l has not been considered. It is the purpose of this thesis to try to arrive at a method in which the above mentioned factors can be taken into account. 45 CHAPTER 3 MODELLING OF SATURATED GRANULAR SOIL UNDER CYCLIC SIMPLE SHEAR CONDITIONS 3.1 VOLUME CHANGE CHARACTERISTICS DURING CYCLIC DRAINED SHEAR. It is a commonly observed phenomenon that a dry granular s o i l compacts or reduces i t s volume when subjected to cyclic loading, e.g., machine foundation vibration and earthquake forces. The general observation from the performance of model as well as existing foundations is that the amount of volume reduction Increases with the magnitude and duration of vibration. Apart from this, the quantitative results obtained by many early investigators using laboratory testing are quite scattered and uncorrelated and sometimes contradictory, see e.g. Mogaml and Kubo(1953)t Barkan(l962) and Bazant and Dvorak(1965). 'Rigid' containers were used in their experiments to enclose the s o i l specimen during the course of vibration so that the state of stress and strain experienced by the specimen was not known. Invariably, the amount of compaction was related to the frequency and acceleration level, or velocity level in some cases, of the input vibration. It i s not d i f f i c u l t to realise that results obtained would be very much apparatus-dependent since the 46 actual stress and strain experienced by the s o i l specimen depends not only on the magnitude of input vibration but also on the dynamic response properties of the container-and-soil system. For example, i f one examines the data obtained by Schaffner(l962), as shown in Figure 3-1. one may observe that each of the f i n a l void ratio curves obtained for a constant acceleration level strongly suggests the presence of resonance of the container-and-soil system at around 80 to 100 Hz. Moreover, the effect of the superimposed vertical component of oscillation, which might be present due to the method of applying vibration, complicated the interpretation of the data. Consequently, i t is more r e a l i s t i c , in the sense of applying laboratory test data to the analysis of actual s o i l structures, to consider the volume change of a s o i l element as a function of the stress and strain experienced by i t . In the dynamic response of a horizontal s o i l deposit to earthquake base motion the cyclic simple shear stress system caused by the propagation of shear waves from the base to the surface of the deposit, as shown in Figure 2-7. is the most significant one. With this in mind many of the recent investigators, e.g., Silver and Seed(197D. Hardin and Dmevich(1972,1972a) and Youd(1971#1972), used the simple shear device to determine the relation between volume compaction and cyclic stress/strain conditions. More recently, Martin, Finn and Seed(1975) studied the fundamental mechanism of s o i l liquefaction and a relationship between the FIG. 3-1 Compaction of Dry Coarse Sand (After Schaffner, 1962) 48 volume change in the drained case and the pore pressure increment in the undrained case during cyclic shear was proposed. The construction of a material model to include both the volume change and pore pressure build-up phenomena is made possible by such a relationship. Before going into constructing the equations of the material model, the various important observations concerning the volume changes under drained cyclic shear are discussed f i r s t . In their investigations into the compaction of sand subjected to cyclic simple shear loading, Youd(1971,1972), Silver and Seed(197l) and Seed and Silver(1972) observed that the amount of volume change (compaction) caused by cyclic simple shear loading on dry granular s o i l depends on, 1. the relative density of the s o i l , 2. the magnitude of the cyclic shear strain, and, 3. the total number of strain cycles the s o i l specimen has been subjected to. Examples of experimental data obtained are as shown in Figures 3-2 and 3-3. It should be noted that the value of the vertical effective stress, which is constant during a particular test, does not influence the volume change characteristics significantly, as shown in Figure 3-2. There is growing evidence that the volume change associated with compaction is caused by the slight re-arrangement of the grains into a more 'stable 1 configuration. The amount of . 5 5 0 I 10 I O 2 I O 3 | 0 4 10 Number of Cycles ^ FIG.3-2 Effect of Vertical Stress on Compaction (After Youd, 1972) i — T " r i i i n r — i " ' i i i 111 D r = 4 5 % / • - 1 1 1 1 1 1 1 1 Cycle 10 - • -- / 6 0 % / -' / 8 0 % / i i i i i i t i i i t i i i i i 0.01 0.1 1.0 10.0 Cyclic Shear Strain - 7^y (%) 3 - 3 Effect of Relative Density and Shear Strain on Compaction (After Seed et al, 1972) 51 grain re-arrangement is governed by how much adjacent grains can move cycl i c a l l y relative to one another, and, since the cyclic shear strain amplitude is directly related to this type of relative displacement the strong dependence of compaction volume change on cyclic shear strain amplitude is not surprising. Furthermore, mention should be made that most of the non-recoverable volumetric strain (grain re-arrangement) caused by one-dimensional compression occurs in the i n i t i a l loading stage and that, 1. during drained cyclic shear, the vert i c a l effective stress is held constant, and, 2. during undrained cyclic shear, the vert i c a l effective stress decreases, so that the effect of vertical effective stress on the •compaction curves* shown in Figures 3-2 and 3-3 which are for constant cyclic strain amplitudes is of minor importance. On the other hand, the vertical effective stress affects the shear modulus and so has a significant effect on the shear strains which result from a given ground motion and therefore on the compaction which results from such an excitation. This topic is discussed in section 3.2. Seed and Silver(19?2) demonstrated the applicability of data obtained from drained cyclic simple shear tests to the determination of total settlement of a model dry sand deposit subjected to harmonic base motion. However, for non-uniform 52 cyclic loading, such as that Imposed by earthquakes, test data in their present form as shown in Figures 3-2 and 3-3 are not directly applicable. Martin, Finn and Seed(1975) proposed that for a s o i l element with given i n i t i a l void ratio (relative density) the total or accummulated volumetric compaction strain, Gyp» can be used as a strain history parameter so that a further increment in volume compaction, A£ yp, caused by an additional strain cycle, f& , can be determined - this i s quite similar to the use of strain hardening parameter used in the theory of p l a s t i c i t y . An expression based on a curve-fitting method was proposed to describe the experimental data:-A e v p = U ( 6 v P ' V"a> O" 1* 2 "' C l ( * a - C2 €vt> ) + T-2"^ ( 3 " l a ) where C^, C^. C^ and are constants determined experimentally. A method of graphically presenting the function given by equation (3-la) is as shown in Figure 3-4 and the values of the four constants were found to be 0.80, 0.79, 0.45 and 0.73 respectively for a Crystal S i l i c a Sand at a relative density, Dr, of 45$. It should be pointed out that the expression is an empirical one and in addition the incremental nature (in discrete steps) of the implied total volume compaction strain function, e v p , has to be accepted since the process of counting the number of shear strain FIG. 3 - 4 Incremental Volumetric Strain Curves (after Martin et al, 1974) V*) 5^ cycles completed is also of a discrete incremental nature. For drained cyclic simple shear tests in which the strain amplitude histories are monitored the volume compaction characteristics can be computed when the four constants associated with the particular s o i l are known. For example, i f is the total volumetric compaction at the end of the 1 i i-th strain cycle and and are the incremental vp s volumetric compaction and strain amplitude respectively for the i-th strain cycle then, for equation (3-1), 6vp = € ^ 1 + *4p (3-3) Accordingly the total volumetric compaction strain after the completion of the N-th cycle can be calculated. In fact, such a computation was carried out by Martin, Finn and Seed(l975) and a comparison between calculated and measured values was made as shown in Figure 3-5. 3.2 STRESS-STRAIN RELATIONSHIP UNDER CYCLIC SIMPLE SHEAR CONDITIONS. From a series of static t r i a x i a l drained tests on a sand Kondner and Zelasko(1963) concluded that the relationship between the axial stress, cr , and axial strain, £ , of granular s o i l under axial compression could be closely represented by a hyperbolic curve, 55 to I o I I I I I £ 0 5 10 15 20 Number of Cycles 1.2 1 c '5 0.8 GO o 'xl •f 0.4 E O > 0 1 1 1 Measured Computed i 1 i 5 10 Number of Cycles 15 20 FIG. 3-5 Comparison of Measured and Computed Volumetric Strain (after Martin et al, 1974) 56 c -e a + b (3-4) where a and b are constants determined from experimental data. The curve given by equation (3-4) i s shown in Figure 3-6a and i t can be seen that the reciprocals of a and b give the i n i t i a l tangent modulus and maximum stress level respectively. In terms of shear stress, T , and shear strain, T , the following equation can be written, see Figure 3-6b, T = 1 . + •* (3-4a) Gmax m^ax where G m a x i s the i n i t i a l tangent modulus, and, T i s the maximum shear stress that can be applied to the sample, Hardin and Drnevich(1972a) showed that for sands G „ and r max max can be obtained from the following expressions, max max = 1230.(2-973 -1 + e 1+K. - s in 4 - f f i T ' .cr; (3^ 5) (3-6) where G m o v and T n a v are i n p s i , max vmax r ' e i s the void r a t i o , (T0 i s the mean p r i n c i p a l e f f e c t i v e stress i n p s i , 0-^ . i s the v e r t i c a l e f f e c t i v e stress i n p s i , K Q i s the c o e f f i c i e n t of earth pressure at r e s t , and, 4>' i s the e f f e c t i v e angle of shearing resistance. FIG. 3-6 Hyperbolic Stress-Strain Curve. 58 Also, equation (3-^ ) can be rewritten into, < f c = i • mv °-7) Where G = T/ /" , and f = T m Q / G m o . and values of f were tabulated for various s o i l s . In addition, they used a hollow cylindrical sand specimen and applied cyclic torsion to both ends of the sample so that a typical s o i l element of the specimen was subjected to a cyclic simple shear stress system. From the test data on a large number of soils they were able to show that the hyperbolic stress-strain relationship can be used for s o i l elements subjected to cyclic simple shear loading and expressions for the shear modulus and damping ratio as a function of strain level were obtained. Furthermore, modifications to the idealized hyperbolic stress-strain proposed were made and are briefly l i s t e d as follows»-I, A hyperbolic shear strain, ^ , i s defined so that instead of equation (3-7) the following expression i s used, (3-8) G i + r. max h with f. . [ l + ae r J (3-9) h 'r 2. Assumptions were made concerning the damping ratio D so that at any cycle and strain level T , 59 D (3-10) max with defined similarly to (3-9) though the constants a and b involved are different from those appeared in equation (3-9). Numerical values for a's and b*s and D m Q V and f are given by Hardin and Drnevich(1972a). However, even with the modifications the expressions are only applicable to the determination of average shear modulus, G, and damping ratio, D, for a complete stress/strain cycle, as shown i n Figure 3-7 and 3-8. Neither the actual stress-strain curve i s traced nor the change in shear modulus due to volume compaction and change in effective stress are taken into account. the liquefaction process the change i n shear modulus due to volume compaction and change in effective normal stress could be taken into account by using the following expression, Martin, Finn and Seed(1975) suggested that during (3-U) hv a + hf where 1*hv i s the horizontal shear amplitude, f i s the shear strain amplitude, (Ty i s the current effective normal stress, and, a and b are functions of volume compaction strain, €vp* as follows t — Fig. 3-7 Equivalent Shear Modulus for Initial Loading ON O FIG. 3~8 Definition of Equivalent Shear Modulus and Damping Ratio for Symmetrical Loop ON 62 1 A2+A3€vp b = B , ^ 2 (3-13) 1 W v p with A 1 # A 2, A^, B l t B 2 and are constants determined from experimental data.. Comparing equations (3-11), (3-12) and (3-13) with (3-4), i t can be seen that i n effect the expressions in (3-12) and (3-13) give the "hardening' effect of volume compaction - i.e., increase in shear modulus and maximum shear stress due to volume compaction. In the next section, a general stress-strain relationship i s proposed based on the above described s o i l behavior. 3.3 MODELLING OF SATURATED GRANULAR SOIL UNDER CYCLIC SIMPLE SHEAR. 3.3.1 Hysteretic Hyperbolic Stress-strain Relationship. In order to arrive at expressions describing the stress-strain behavior of an idealized s o i l element under cyclic simple shear conditions i t i s assumed here that the Masing(1926) type of hysteretic relationship i s followed and that a hyperbolic stress-strain relationship similar to that proposed by Martin, Finn and Seed(1975) i s used. These are discussed in the following paragraphs. For the i n i t i a l loading phase of the granular s o i l the hyperbolic stress-strain relationship formulated by 63 Kondner and Zelasko(1963) and also adopted by Hardin and Drnevich(1972) i s used. Therefore, the i n i t i a l response of the granular s o i l to simple shear i s given by the following equation, Gmo^ T = 1 + mo' mo (3-14) where T i s the shear stress, "t i s the shear strain, G ^ i s the i n i t i a l maximum tangent modulus and TL« i s mo mo the maximum shear stress that can be applied without failure. Following the suggestion of Hardin x and Drnevich i t i s assumed that these quantities G and TL n are determined ^ mo mo by the in-sit u condition of the granular s o i l and are given by the following expressions, G m o = llt76o,(2.973 - e ) 2 J -mo 1 + e 1 u o mo 1+K„ 2 1-Kn 2 ( - y ^ i n * ' ) - ( ) > 1 * (3-15) (3-16) where and l " are in psf, mo mo r * e i s the void ratio, <TQ i s the mean principal effective stress in psf, is the vertical effective stress in psf, K Q i s the coefficient of earth pressure at rest, and, 64 4>' i s the effective angle of shearing resistance. Since we are considering the response of horizontal layers only, i t w i l l prove more convenient to express the stress conditions in terms of the vertical effective stress, a~y» Then equation (3-15) i s expressed in the form, Therefore, for a given in-situ condition, Gmo = c*/^ (3-18) = pjffi (3-19) ''•mo For unloading and reloading the Masing(1926) type of hysteretic behavior i s used. This assumes that i f the i n i t i a l loading curve, or skeleton curve, which i s described by equation (3-1*0 can be represented by, X s = f 2(D "(3-20) Then the unloading or reloading curve can be obtained from, where (f ,T_„) i s the last reversal point in the stress-strain plot. Consider a s o i l element being strained from i t s undeformed state to a shear strain level fa, with 81 65 .corresponding shear stress level I _. The element i s now unloaded and reloaded in the opposite direction so that the shear strain or stress undergoes a complete reversal when the shear strain reaches a f i n a l value of f-^, with corresponding shear stress, T^. The stress-strain diagram i s shown i n Figure 3-9a. Curve OA i s given by (3-1^ ) while curve AB i s given by, T g ~ = f 2 ( f~2 T&) (3-22) The reversal curve AB intersects the negative branch of the skeleton curve at f^'-fa. a n d when i s increased beyond fa i t i s assumed that the stress-strain path followed would be the continuation of the negative branch of the skeleton curve, BC, with (/, T Q) erased from the stress-strain history. In other words, when the last reversal point i s on the skeleton curve the absolute value of the stress or strain can be used as a limiting value for .switching back to the skeleton curve, thus 'erasing' (/"_, T ) cl S a from our calculation. A general loading history i s shown in Figure 3-9b. Assuming that the current loading would bring the stress-strain point, P to P', i t can be seen that two stress-strain branches, AB and BP'A' are required. The reversal points ( f a . T s a ) . (^D.^sb) and ( T c.X s c) have to be recorded for the determination of the current branch of the stress-strain curve, CPBP'A'. (/„,!) i s required for the C S C branch CPB, A* (a) first unloading FIG. 3-9 Hysteretic (b) general reloading Characteristics 67 2V 2 (3-23) (fo.X^o) i s required f o r the branch BP'A', a sa 2 = l2' — (3-24) While (^b'tgt)) i s u s e d a s a c r i t e r i o n f or switching from (3-23) to (3-24). Note that as soon as the stress s t r a i n point P passes the point B the reversal points B and C would be case, the stress or s t r a i n value of the l a s t but one rev e r s a l point can be used as a l i m i t i n g value f o r switching back to the l a s t but one unloading or reloading branch and also as a c r i t e r i o n f o r 'erasing' the l a s t and the l a s t but one reversal points from future c a l c u l a t i o n . I t should be pointed out that the above description constitutes a complete s t r e s s -s t r a i n r e l a t i o n s h i p under general loading and unloading sequence. Herrera(1964) has shown that combinations of Coulomb f r i c t i o n and l i n e a r e l a s t i c elements e x h i b i t t h i s type of behavior. Newmark and Rosenblueth(1971) suggest the Masing model to be a good representation of s o i l behavior and some experimental support for the suggestion i s provided by Marsal(1963). 3.3.2 Volume Compaction and Pore Pressure Increments. A t y p i c a l one dimensional load-deformation curve for a l a t e r a l l y confined sand specimen i s shown i n Figure 3-10a where ab represents loading, be unloading and ce reloading. eliminated from our c a l c u l a t i o n at the same time. In t h i s Log P (a) FIG. 3-10 One-Dimensional Soil Model ON 00 69 It i s generally observed that while ab i s quite different from be and cd the latter two are reasonably close together to allow their representation by one curve, ce. A conceptual model as shown in Figure 3-10b can be used to describe, to some degree of approximation, the behavior of sand under one dimensional compression. The spring element in the model can be assumed to be non-linear so that the idealised unloading and reloading curve ce can be simulated. The s l i p element in the model should be of a 'hardening' type with a load-deformation curve as shown in Figure 3-11&, i.e., the force required to cause s l i p increases with the magnitude of sl i p displacement. Using such a conceptual model i t can be seen that the volumetric strain, € v > i s given by the following equations;-€ v = < Klr + Ks>'P P*Pmax <3-25> A £ v = K l r.A P p < p m a x (3-26) where K l r i s the modulus (non-linear) of the spring element, K_ i s the pseudo-modulus of the s l i p element, and, s p „ i s the past maximum vertical effective stress, 'max When horizontal cyclic shear i s applied to a sand sample at equilibrium with an vertical load, e.g., cyclic simple shear test, volumetric strain occurs at constant vertical load so that a load-deformation curve as shown in Figure 3-12 i s obtained. Since the additional volume reduction caused by cyclic loading i s non-recoverable an additional feature has to FIG. 3-11 Force-Displacement Relationship of Slip Elements -s3 O 71 Fig. 3-12 Volume Compaction Caused by Cyclic Loading 72 be incorporated into the s l i p element of the conceptual model - namely, the addition of s l i p component caused by cyclic shear strain as described in previous section. Equations (3-25) and (3-26) have to be replaced by the following, *v = < Klr + Ks>'P + <Srp P*Pmax ^-27) A € y = K l r.AP p - c p J J ^ (3-28) where € i s defined as in equation (3-1). Since compaction V £J changes the properties of the s o i l element, e.g. void ratio and shear modulus, more or less in the same fashion as the maximum past pressure, P max' a n e c l u ^ v a i e n " t maximum past pressure, p!LV, i s used in equation (3-27) and (3-28) instead max °f P™<,v. The determination of p*~QV i s shown conceptually in max max Figure 3-12. To account for the undrained condition i n a s o i l element, i t may be assumed that the conceptual model shown in Figure 3-10b i s f i l l e d with water and drainage prevented. Under these circumstances, volume reduction cannot occur at the same time when an increment in volumetric compaction (grain slip) strain, AG , i s induced by a shear strain cycle. As a result, elastic strain in the spring element has to be relaxed causing a portion of i t s stress transferred to the pore water. Description of this phenomenon at the grain level can be obtained from Martin, Finn and Seed(1975)» and an expression for the pore water pressure increment, A(TW» was 73 found to be, . A C 7" W = ( x ' nK l r j » Klr' A €vp -^29) K w where n i s the porosity of the s o i l element, K w i s the bulk modulus of water. When water can be assumed to be incompressible or K w > > K i r » equation (3-29) can be simplified to, ACT w = K l r.A € v p (3-30) It should be pointed out that during the settlement or liquefaction process of a hroizontal deposit the past maximum vertical stress i s never exceeded and equation (3-27) i s not needed at a l l . This eliminates the determination of one-dimensional loading characteristics of the granular material under consideration, nevertheless, the one-dimensional rebound characteristics, K^r, has to be quantitatively known in order that equation (3-22) can be used. In the determination of the magnitude of volumetric compaction strain due to horizontal cyclic strain expressions (3-1) and (3-la) are used. Recently, Pyke(1973) obtained experimental data for a Monterey Sand at various normal pressures and relative densities. The independence of volumetric compaction strain on vertical stress under the conditions described in the previous paragraphs was again confirmed. A comparison of the 6 ^ function for sands of 74 different relative densities can be made by plotting data obtained by Pyke(1973) into one graph. A 'reduced' incremental volumetric compaction strain, ^C^p, i s obtained by the product of relative density in per cent, Dr, and incremental volumetric compaction strain, AC^, i.e., vp D .A€ r vp (3-31) The graphs are plotted in Figure 3-13. It may be seen that the curves f a l l more or less to an unique pattern suggesting the possibility of using one shape function, U s, for the incremental compaction strain for sands at different relative densities, namely, A € v p = ^i-W^ (3-32) vp (3-33) (3-34) where ^ , ^ 2, ^ and ^ are constants similar to those used in equation (3-la). In this usage i s inversely proportional to the relative density value, Dr, while i(/2, ^ and ij;^ are assumed to be constant for the same type of sand at different densities. Martin, Finn and Seed(1975) calculated the volume change increments by assuming that these occur at the completion of a strain cycle as shown in Figure 3-l4a. However, i t i s generally observed that in a drained cyclic 0.1 0.2 Shear Strain Amplitude — *f % FIG. 3-13 'Reduced' Incremental Volumetric Strain ->3 77 shear test on granular s o i l gradual volume compaction (or pore pressure build up in the undrained test) occurs during the stage of unloading from maximum shear strain. In view of this, i t i s not unreasonable to assume a linear variation of volume compaction increment with the reduction of shear strain from i t s maximum value as shown in Figure 3-l4b. The magnitude of volume change increment i s s t i l l assumed to be determined by an expression similar to that of equation (3-3*0 though the value of i f ^ i s halved in order that the volume change increment i s added every half cycle. For computations involved in this thesis the numerical values obtained by Martin, Finn and Seed(1975) were used to determine the values of q/2, y 3 and fy. 3.3.3 Hardening and Effect of Pore Pressure Rise In a study on the fundamental behavior of sands under cyclic loading in simple shear, Martin, Finn and Seed suggested that equation (3-11) might be used to represent the stress-strain behavior. In this equation the parameters a and b are constants for a given load cycle but in general are functions of the volumetric compaction strain, ^yp. and are given by equations (3-12) and (3-13). The maximum value of the shear modulus, G m n, in any cycle of loading n can be obtained in the following manner, Gmn = < d Thv / a t * = 0 (3-35) 78 From equation (3-H)t therefore, we obtain Gmn = JK'a (3-36) For the case of i n i t i a l loading when no hardening has yet occurred we obtain Gmo = ' J o T / A l (3-37) After n cycles, substitution for a from equation (3-12) gives Gmn " / A, - — € v p 1 A2 + A3 €vpJ (3-38) in which € i s the volumetric strain accumulated to the nth cycle of loading. With algebraic manipulation of equation (3-12) and substitution of G m Q for fcr^ / A 1 according to equation (3-37) we obtain in which a^ and a 2 are constants. We can obtain an expression for the maximum shear stress from equation (3-11) by considering very large strains. The maximum shear stress, T , in the nth cycle of loading as t—oo i s from equation (3-11), = ! ^ / o (3-40) During i n i t i a l loading when no hardening has taken place we 79 obtain and a comparison with equation (3-19) gives B l = 1 /PlK (3"42) Following a procedure similar to that used in deriving G m n we find that after n cycles of loading 6. tan " Tmo l1 + b. + IZt] P^3) L 1 2 VB -I where again i s the accumulated volumetric strain and b, vp 1 and b 2 are constants. It should be noted that Martin et al(1975) used six constants to define the stress strain relationship of sand. Two of the constants A. and B1 were used to determine G ^ and 1 1 mo X and the remaining four constants were used to define the strain hardening. Herein only four hardening constants are used because G m Q and T m 0 are determined independently using the Hardin-Drnevich equations. The four hardening constants, a 1, a 2, b^ and b 2, are obtained by f i t t i n g equations (3-39) and (3-^ 3) to the results of constant strain cyclic loading tests. In the case of undrained shear the increment in pore water pressure, A<7"w, may be computed using equations (3-30) and (3-34). In one load cycle the total change in effective 80 stress w i l l be Ac^, and the e f f e c t i v e stress changes from <r^ .0 to cr^ = <fy0 - Acr w. The new l e v e l of e f f e c t i v e stress w i l l also a f f e c t the maximum shear modulus, G , and the maximum mn shear stress, '^xm* applicable to the next cycle of loading. I t i s d i f f i c u l t to apportion the hardening between grain s l i p and volume compaction. Finn et al(1970) suggests that the major contribution to hardening comes from the creation of more stable grain to grain contacts due to c y c l i c s l i p at contact points rather than from increased density. Although the matter cannot be considered s e t t l e d , herein, we attri b u t e most of the hardening to grain s l i p and so we conclude that hardening occurs also i n undrained shear. In the case of undrained c y c l i c shear response the grain s l i p i s represented by the p o t e n t i a l volumetric s t r a i n , ^ -yp* For undrained sands, therefore, the maximum shear moduli and maximum allowable shear stresses f o r the nth cycle are related to the i n i t i a l values by the following equations: mn "mn = G. mo 'mo VP a T F 2fcvpJ [^voj 1 + K—Z" L b i + vp V v p u v 'vo (3-44) (3-45) i n which cr^Q i s the i n i t i a l v e r t i c a l e f f e c t i v e stress, and ^ i s the current v e r t i c a l e f f e c t i v e s t r e s s . The r e s u l t i n g s t r e s s - s t r a i n r e l a t i o n s h i p i s then given by the following: 81 G r m n . (3-46) G f 1 + m n 3.4 COMPUTATION OF SOIL BEHAVIOR USING THE PROPOSED SOIL MODEL. ~ . In this section the values of stress and strain and other related variables, e.g., shear modulus, damping, volume compaction and/or pore pressure build-up, etc., are computed based on the equations proposed in the previous section and comparison i s made, whenever possible, with experimental data. A hypothetical sandy s o i l i s assumed to have the following properties:-e = 0.83 ; * = 36.0° ; K Q = 0.50 ; (Ty0 = 1600 psf. With these and with the help of equations (3-5) and (3-6) the i n i t i a l maximum shear modulus, G , and i n i t i a l maximum shear mo stress, T m Q, can be calculated. To calculate the compaction characteristics values of 0.400, 0.790, 0.563 and O.730 are used for i j ^ , q>2 • ^3 a n d respectively. These values were obtained based on the discussion given in section 3.3. 3.4.1 Hypothetical Drained Cyclic Simple Shear Test. In the drained cyclic simple shear test the 82 effective stress i s kept constant. Consequently, equation (3-44) and (3-45) can be reduced to the following!-The values of a 1, a 2 > b^ and b 2 take the values of 0.754, 0,406, 0.550 and 0.500 respectively as discussed in section 3.3. By using a strain increment of 0.005$ a drained cyclic simple shear test with constant shear strain amplitude i s calculated. The procedures and equations described in section 3.3 are followed. The f i r s t four stress-strain hysteresis loops were computed and the results are shown in Figure 3-15. It can be observed that the effect of hardening i s well represented by the increasing level of shear stress required to bring the sample to the same peak strain at each cycle. The average shear modulus or secant modulus i s also computed for various shear strain levels and compaction strain values as shown in Figure 3-l6. Similar relationships in terms of loading cycles are shown in Figure 3-17. The data shown that at larger strains hardening has a proportionately larger effect on the average shear modulus. For example, at 0.5% strain the average shear modulus after 5 cycles i s twice that after 1 cycle while at 0.05% i t is only 10% greater. 600. Fig. 3-15 Hyperbolic Stress-Strain Loops (Drained Test) Showing Effect of Hardening .02 .05 0.1 0.2 0.5 Cyclic Shear Strain Amplitude (%) Fig. 3-16 Shear Modulus as a function of Cyclic Shear Strain and Volumetric Strain Fig. 3-17 Shear Modulus as a function of Cyclic Shear Strain and Number of Cycles CO • V n 86 These observations indicate the d i f f i c u l t i e s and approximation inherent in selecting a single shear modulus-strain curve for using in dynamic analysis of earth structures composed of granular materials. The damping ratio defined by (see Figure 3-B). D = 1 area of loop APQ (i LQ) ' area of AOAB v i s also calculated. For convenience of evaluating the area of the stress-strain hysteresis loop the volumetric compaction strain i s assumed to be constant so that the following expression for the damping ratio can be obtained D = | [l + i ] [ l - log e(l+Y )A] - -f- (3-50) G f where Y = -22- (3-51) Sin The variation of damping ratio, D, with shear strain amplitude, f , at different values of compaction strain, € , are calculated and the results shown in Figure 3-18. Pyke(l973) has presented data on the variation of average shear modulus and damping ratio with strain for Monterey #0 Sand obtained from the Geonor simple shear apparatus. These experimental data are shown in Figures 3-19 and 3-20, It can be seen that the behavior of the proposed s o i l model, namely shown in Figures 3-16 and 3-18 bear remarkable similarity to the experimental data. FIG. 3-18 Computed Damping Ratios FIG. 3*19 Shear Modulus as a Function of Shear Strain and Total Compaction (after Pyke, 1973) Damping Ratio 68 90 3.4,2 Hypothetical Undrained Cyclic Simple Shear Test. In this section the assumed s o i l data mentioned in the previous section i s used. In addition, the power law obtained by Martin, Finn and Seed(1975) i s used for the one-dimensional rebound stress-strain relationship, namely, € v r = ^ K o ^ K ^ (3-52) where i s the elastic (rebound) volumetric strain, <fVQ i s the past maximum effective vertical stress from which unloading i s taking place, <7-y i s the current effective normal stress, and, m, n and k 2 are constants taken to be 0.430, 0.620 and 0.0025 respectively. The set of unloading curves given by (3-52) i s shown in Figure 3-21. The one-dimensional rebound modulus can then be obtained by differentiating equation (3-52), K vr (cr;) 1"" 1 / mk 2(<r; 0) n- m (3-53) The procedures in computing an undrained cyclic test are similar to those used for drained test. There i s one additional step involved i n the calculation and that i s the calculation of pore water pressure increments using equation (3-30). Small strain increments of 0.005% a r e again used. 91 However, to il l u s t r a t e the softening behavior of the s o i l model during undrained cyclic shear a constant stress amplitude cyclic test i s chosen. The f i r s t five stress-strain loops are shown in Figure 3-22. The decrease in shear modulus due to rise in pore water pressure i s indicated by the increases in strain amplitude generated by successive cycles of constant amplitude shear stress. A schematic representation of the development of shear strain and pore water pressure during cyclic loading are shown in Figure 3-23. The progressive increase in the amplitude of cyclic shear strain and the double curvature in the pore water pressure development curve, which are very characteristic of actual granular s o i l under undrained cyclic shear, are also obtained from the model. Recoverable Volumetric Strain - e v r FIG. 3"2I One Dimensional Unloading Curves (after Martin et al, 1974) 93 FIG. 3 -22 Hyperbolic Stress-Strain Loops (Undrained Test) Showing Softening due to Pore Pressure Rise Time (sec) FIG. 3~23 Strain and Pore Pressure during Undrained Cyclic Shear Test 95 CHAPTER 4 MECHANICAL MODEL FOR THE ANALYSIS OF LIQUEFACTION 4.1 INTRODUCTION From the discussion of s o i l liquefaction in earlier chapters i t i s seen that the pore water pressure, <r , generated during shaking controls the occurrence and extent of s o i l failure. In view of this, i t i s necessary to set up an analysis of liquefaction potential to calculate the maximum rise of pore pressure, or maximum reduction in effective normal stress, in a saturated horizontal sand deposit when subjected to an oscillatory base motion. As an approximation, liquefaction failure i s assumed to take place when the effective normal stress, o~y = <r^.o-0^, i s reduced to zero, or, for practical purposes, to less than five per cent of the i n i t i a l effective normal stress, 0~yQ. A complete analysis of liquefaction should include the following processes, namely, 1. Dynamic response to base motion, 2. Increase in volumetric compaction strain due to cyclic straining, 3. Pore water pressure generation as a result of (2), and, 4. Dissipation of excess pore water pressure(consolidation). There are two main sources of non-linearity i n the 96 proposed analysis of liquefaction: the non-linearity of the s o i l stress-strain relationship and that which arises from the coupling of the above-mentioned four processes. As a result of the complex coupling of the processes, which w i l l be evident in the following sections when the whole problem i s formulated, numerical procedures have to be used for the solution of the equations obtained. For convenience, equations for numerical computations w i l l be derived immediately after each 'exact' equation i s obtained. 4.2 FORMULATION OF THE PROBLEM. A level sand deposit i s assumed to be of in f i n i t e extent so that propagation of shear waves through the deposit can be treated as a one-dimensional problem. Consequently, without any loss of generalization considerations can be confined to a s o i l column of unit cross-section through the deposit, as shown in figure 4-1. Soil under the water table i s assumed to be completely saturated so that when flow of water occurs as a result of s o i l compaction the direction of flow i s upward towards the ground surface. Governing equations and relevant s o i l parameters w i l l be introduced and discussed in the following sub-sections. 4.2.1 Dynamic Response to Base Motion Consider a s o i l element of thickness dz located at height z from the base of the deposit as shown in Figure 4-1. Let x be the horizontal displacement, 1 the horizontal shear FIG. 4-1 Idealization of the Response Problem 98 stress and y the t o t a l mass density of the s o i l element. Dynamic equilibrium requires that, 2 p . d z . - | | . d z . l 0 , 3t or, . || - .o (4-x) For v i s c o - e l a s t i c materials i n which the s t r a i n and s t r a i n rate e f f e c t s are additive and/ non-coupled the s t r e s s - s t r a i n r e l a t i o n s h i p can be represented by, T = f ^ i T ) + f 2 ( f ) (4-2) where f i s the shear s t r a i n , jr i s the shear s t r a i n rate, and f^ and f 2 are functions of f and f r e s p e c t i v e l y . In the case where the s o i l properties are constant with depth equation (4-2) can be substituted into (4-1) d i r e c t l y to give, Since if = i f , d Z 1£ = 2 2* .sir _ A and (4-4) 9z dz ^ , and wri t i n g t[Ct) = -$r> flCf) , and f 2 '(r) = ^ . f 2 ( / ) , (4-5) 99 Equation (4-3) can be rewritten in the following form by using the relations given by (4-4) and (4-5) »-2 2 px - f' l y i ) - f' 2_(x) = 0 (4-6) J dz 3z t t It can be seen that when f^=0 and ^2~G* w n e r e G represents the shear modulus, the familiar equation for the dynamic response of an undamped linear shear beam results, namely, j>x - G . ^ - J ( X ) = 0 (4-7) In general, the s o i l properties are functions of depth, (H-z), so that the stress-strain relationship given by equation (4-2) has to be replaced by a more general equation which contains s o i l parameters as functions of depth also. In this case an equation similar to (4-6) cannot be obtained in a straight forward manner since the partial differentiation carried out in equation (4-3) would involved more terms. It i s more advantageous to apply numerical approximations to equations (4-1), It i s shown in Appendix I that through the use of the fini t e difference method equation (4-1) can be discretized in the z-domain to give, rm. [ M ] ( X } + [C]{X) + [K]{X> = -L m n .'X b (4-8) so that a lumped mass system as shown in Figure 4-2 can be used for approximate analysis. The components of the matrices given in the Appendix I are rewritten here for 100 X-hi m, n "'I "•^j'a&5S.'^5'iS^~ m; — — — — — — — — c i »ki wmmmmmm. h 3 h 2 l-c 3,k 3 m-. 2 ^2* K 2 Layer I FIG. 4 - 2 Approximation by Lumped Mass System 101 convenience :-[ M ] = rm1 0 0 0 " 0 m2 0 * 0 0 0 • 0 • • • o nj [c] = ' C l + C 2 -c 0 0 • • 0 c2+c3 -c3 0 • • 0 0 c3+c4 • • 0 0 . -c c <* n nj [K] = V k 2 "k2 0 0 • • 0 ' "k2 k2+k3 0 • « 0 0 • • 0 0 • § • • -k • k n (4-9a) (4-9b) (4-9c) where = fyi-*bihi n c i + % + i b i + l h i + l i = 1 » 2 n - 1 = fA'ti A ) .b,/(h. .f, i) , and 1 1 - 2 1 - 2 1 1 1 - 2 = foU* i.z. i).b./(h..f. i) . C 1—2 1—2 1 1 1—2 It should be pointed out that h^ and b^ are the thickness and width of the i-th layer while p. i , f,(f- i,z- i ) and J 1 - 2 1 1 - 2 1 - 2 fo(fi j,) are values referred to the centre of the 1-th C- 1 - 2 1 -T2 102 layer. 4.2,2 Hyperbolic Stress-Strain Relationship, f^, and Material Damping Values, f 1 . It was suggested in chapter 3 that for a s o i l element under simple shear condition, such as in the present case of a horizontal s o i l deposit transmitting shear waves from the base to the surface, a hyperbolic hysteretic stress-strain representation i s a close approximation to the actual stress-strain. Consequently, the, fg function described in section 4.2.1 takes the form of equation (3-46) together with the Masing behavior given in section 3.3.1.1. When vibrational energy i s being transmitted through a material medium a portion of i t i s dissipated internally due to a number of mechanisms at the grain level, e.g. Coulomb fr i c t i o n , plastic deformation, viscosity, dislocation, etc. This dissipation^ of energy which causes a decrease in the amplitude of vibration i s generally called material damping or internal damping. Generally, material damping is classified into two main categories - a frequency dependent type and a frequency independent type. Incidentally, i f one integrates equation (4-2) over one strain cycle, assuming that the stress-strain loop i s similar to the loop AB as shown in Figure 3-3a and that the integrals involved are integrable, one should be able to obtain the energy dissipated during the strain cycle ABA, Wd , as 103 JX .df = y q c h . d f + jf2(f).df (4-10) or, wd(ir,f) = w'(ir) + vq(0 (4-n) Equation (4-2) can be used to describe a simple linear visco-elastic material, e.g., a Kelvin solid where, fAf) = c.'f , 1 (4-12) and f 2(jf) = K.r . where C and K are the visco-elastic constants. For harmonic vibration, f=Asin«/t, the energy dissipated per cycle, Wd, can be obtained from equation (4-10). Since, M j f^'f) .df = Jc.f.r. dt o 2 and j f 2(T) .df = 0 , giving Wd = C7TWA2 ^ ( 1^_ 1 3 ) It can be seen that the energy dissipated per cycle, or the damping depends on the frequency w . For this particular case of f^ where the damping stress i s linear with strain rate, the term viscous damping applies. For f r i c t i o n a l materials, the components of the constitutive equation (4-2) can be written as follows :-t^'f) = 0 , 104 fM) = +f for df>0 ^ y f 2 ( f ) = -f for df<0 where f i s a positive quantity. The energy dissipated per cycle can "be obtained similarly by (4-10) to give, wd = f f 2 ( n . d f = 4.f .A (4-i4) This type of damping which is called Coulomb Damping, i s frequency independent. It can be observed from -6he above two examples that the constitutive relation (4-2) can be used to include frequency dependent damping, for which the strain-rate dependent stress component i s responsible, and frequency independent damping, or structural damping, which i s associated with the non-linear hysteretic stress-strain under static loading. • >- • In the equivalent linear analysis where a non-linear hysteretic material i s represented by a linear visco-elastic model, damping has to be introduced a r t i f i c i a l l y through W as W" vanishes. In this case the equivalent viscous damping factor, C , has to be obtained by equating the energy eq dissipated per cycle in the non-linear material and i t s linear equivalent, so that, 105 In a truly non-linear analysis, as i s in the present analysis, where the hysteretic stress-strain relationship i s used, a r t i f i c i a l viscous damping coefficient given "by (4-15) i s no longer needed, instead, the function f^ can be used to describe the actual viscous damping effect present in the material. Or, in a fashion similar to that of the method using equivalent viscous damping, i t can be used to compensate for the insufficient damping caused by the approximation of the actual stress-strain by the function f^(f). Note that the uncoupling of Wd(f,T) into W and W" i s not a necessary one though i t brings about convenience in discussion as well as in the experimental determination of damping values. In this uncoupled form the f i r s t type of damping can be obtained through viscosity measurements while the latter type can be obtained by testing the material in question in a quasi-static manner. For the present purposes, i s assumed to be of the linear viscous type, i.e., f x ( f ) = Cf (4-16) to take into account of the damping effect caused by the flow of water occurring as a result of the deformation of the s o i l skeleton, see e.g. Biot(1956). Yen(1967) attempted to measure the viscosity of saturated sands near liquefaction. However, i t i s not clear whether the viscosity values obtained referred to the actual strain-rate dependent stress component 106 or to the ' a r t i f i c i a l ' equivalent viscous damping. 4.2.3 Volume Compaction and Pore Pressure Generation. Figure 4-3b shows the material model that was proposed in chapter 3. The spring-and-slip element represents the s o i l skeleton while the flu i d inside the container represents the pore water of the granular material. It i s assumed that the total volumetric strain, £ v , i s composed of two components, namely, the rebound (elastic) component, and the compaction (slip) component, ^vp» i.e., € = € + € (4-1?) v vr vp . " Incremental volumetric compaction strain, &€.y^> i s given by the following expression, Aevp = +1-e 2 ( f a " * 2 V + * 3 r a + ^ . (4vi8) where fQ i s the amplitude of shear strain, or the maximum shear strain level reached, in the last strain cycle (or half strain cycle) and the vj/s are s o i l constants determined from cyclic simple shear or t r i a x i a l tests. In the case of dry granular material an application of a shear strain cycle with amplitude f w i l l cause the total volumetric strain to change by an amount equal to that of the volumetric compaction strain, AE^, since there i s n_o change in the elastic component when the effective normal stays constant, i.e., Grd. Surface 02 0 - w + ^ d Z | oz l i t Base Rock , (a) dz / / / / MR ( b ) FIG. 4 - 3 Consolidation Model 108 A € v r * ' 0 • (4-19) and A£y = &vp. I f the s o i l i s saturated and drainage cannot take place instantaneously, relaxation of the elastic strain component w i l l take place. The result i s that a portion of the vertical stress i s transferred to the pore water so that an increment of pore water pressure, A(Tr, occurs, as explained in chapter 3 . A<TW - K i r ' A € v p ( 4 - 2 0 ) where K< i s the one-dimensional'rebound modulus of the s o i l l r skeleton. Detail discussion of the volumetric strain caused b,y cyclic loading i s given in section 3 . 3 . 2 and i s not repeated here. 4.2.4 Reconsolidation. The equation for reconsolidation with pore water pressure generation as a result of internal volumetric compaction can be derived in the same manner as that for the conventional consolidation equation. The assumptions used i n the derivation are i -1. Complete saturation. 2. Negligible compressibility of pore water. 3 . One dimensional compression. 4. One dimensional flow of water. 109 5. Validity of Darcy Law. Consider a typical s o i l element of unit cross-section as shown in figure 4-3a, where v i s the velocity of flow of water into the s o i l element. For the continuity of flow, the difference in volume of flow of water into and out of the element should equal the volume change of the element i t s e l f , so that we have, v + & z - vldt = f^.dz.dt £ - U <*-2L> I f the flow satisfies Darcy's Law, v = -k where k i s the coefficient of permeability, f i s the unit weight of water , and w tj i s the excess pore water pressure. W I f the expression for v i s substituted into equation (4-21) the following equation i s obtained, W From equation (4-17) i t can be seen that, Since, A e ^ = jp - .A<r v l r 110 = V - (A(T - ACT) (4-25) l r During the liquefaction of a horizontal s o i l deposit the total overburden pressure i s not changed, i.e., so that equations (4-23) , (4-24) , (4-26) and (4-27) can be combined to give, + K - ^ ' (4-28) When the coefficient of permeability, k, i s constant with depth and i f ( Ki r« k/^ w) i s denoted by c y, equation (4-28) t can be simplified to, •Wvt ^ v r . Equation (4-28) or (4-29) i s the desired equation governing the dissipation of excess pore water pressure with continuous pore water pressure generation due to internal volumetric compaction. It can be seen that by setting the volumetric compaction strain function identically equal to zero, i.e., €^ =0, the well known Terzaghi's Equation of Consolidation i s obtained. Inspite of the extra term in equation (4-29), i t i s s t i l l a typical parabolic partial d i f f e r e n t i a l equation in one I l l dimension and i s identical in form to that for the conduction of heat in a bar with heat sources distributed along the bar. For example, Sneddon(1957) gave the following equation for the conduction of heat in solids, = (cv 2 0 + H(r , e,t) (4-30) 3t j>c where Q i s the temperature at a point with position vector r, k i s the thermal conductivity, j). i s the mass density of the material, c i s the specific heat of the material, and, K = k / ( p c ) . The function H(f,0,t) gives the amount of heat generated per unit volume per unit time. Physically (H/j)c) gives us the rate of rise of temperature when there i s no conduction. This i s quite consistent with the physical meaning implied in the last term of equation (4-29), namely, -the rate of rise of pore pressure when there i s no drainage! Note that with the help of equation (4-18) the last term in equation (4-28) or (4-29) can be theoretically evaluated so that the internal volumetric compaction strain can be expressed as a function of the total number of shear strain cycles (or half cycles) to which the s o i l element has been subjected, Tf^ . The f i n i t e difference scheme for equation (4-29) i s derived in Appendix I I I . When a step by step integration method is used to determine the dynamic response of the s o i l deposit the shear strains can be computed for a l l 112 the layers for each time increment and the pore pressure generation and dissipation can he incorporated. Discussion of the solution procedures is vgiven in the next section, 4.3 SOLUTION SCHEME. A numerical scheme using a step by step method i s used to solve the governing equations. In equation (4-8) since [K] i s a function of the displacement vector, {X}, a direct integration procedure i s used so that the displacement, velocity and acceleration values at any time t can be related to values at (t-At), where At i s the time step used in the direct integration procedure. There are several methods for numerical integration of equations similar to equation (4-8) and a comparison of the methods based on their s t a b i l i t y and accuracy has been given by Bathe and Wilson(1973). For the present prupose Newmark's method(1959) of step by step integration was chosen for i t s unconditional s t a b i l i t y in the analysis of linear systems. The particular equations and procedures used in the actual computation are derived and discussed in Appendix II. Also, in equation (4-29), the volumetric compaction function, € , i s handled as a discontinuous function, which has an incremental increase in value every time the shear strain, f , associated with the same layer decreases from a maximum value, f„n„, as shown in J max Figure 3-14. In this case the f i n i t e difference method of solution becomes quite suitable. The relevant equations and 113 expressions arising from the use of the f i n i t e difference method of solution for equation (4-29)t or equation (4-28) for the more general case, are given in Appendix III. A brief outline of the step by step solution scheme can be given as follows :-1. Based on the current values of f, and <ry at time t, the shear modulus i s calculated as described in section 4.3.2. 2. The matrix [K] in equation (4-8) i s then updated. 3. With the base acceleration value at (t+At) new values for {X} , {X} , {X} and {f} can be calculated using a direct integration method. 4. When a layer has i t s shear strain decreasing from the maximum value i t s volumetric compaction strain i s incremented according to the equations described in section 4.2.3. 5. Current value of pore water pressure, j * w , are calculated using the fi n i t e difference method as discussed in Appendix III. 6. Current effective normal stresses, g-^, are then calculated "by 0~v- (7"v-(7^ , for use in the next time step. A computer program has been written based on this solution scheme and the overall flow chart i s shown in figure 4-5. A l i s t i n g of the computer program i s given in Appendix 114 READ IN DATA NUMERICAL INTEG-RATION FOR X, X, X & f ETC. PORE PRESSURE, VOLUME COMPAC-TION, IP REQD. RECONSC CALCUL/ NIDATION LTION CHANGE SHEAR MOD-ULUS & OTHER PARA-METERS. PRINT OUTPUT &/0R TAPE OUTPUT FOR PLOTTING. MORE [NCREMENTS ?. YES NO FIG. 4 -4 Flow Chart 115 CHAPTER 5 APPLICATION OF THE MECHANICAL MODEL - AN ANALYSIS OF LIQUEFACTION OF A SOIL DEPOSIT 5.1 A HYPOTHETICAL SOIL DEPOSIT The method of analysis developed in earlier chapters was applied to a hypothetical s o i l deposit described in this section. The same deposit was also analysed by using the equivalent linear analysis so that comparison can be made. Assumed values for the deposit are as shown in Figure 5-1. The variation of relative density with depth was made similar to that of the Niigata deposit analysed by Seed et al(l96l) so that the problem is as close to real situation as possible. However, mention should be made that i t is not the intention of this thesis to project the result of the analysis to the actual f i e l d condition at Niigata. The maximum and minimum void ratios were assumed to be 1.00 and 0.50 so that the void ratios at various relative density can be calculated. For convenience, the bulk density of the s o i l above water table, which was 5*0 feet beneath ground surface in this case, was assumed to be 122.0 pcf while the buoyant unit weight was assumed to be 60.0 pcf. The s o i l deposit was sub-divided into 14 layers in order that application of the previously described method of analysis "T^5J0 ft Ground Surface o ===== or Water Table SAND 200.0ft T = 1 2 2 0 PCf ' Bedrock ' ^ ^ ^ ^ ^ ^ ^ (a) Tbuoy. s 6 0 0 P c f H 150 50 100 a. UJ a 2001 FIG. 5-1 Properties of a Soil Deposit could be made. The properties of each layer are as shown in Table 5-1. Notice the numbering of the layers is from bottom up. Effective shearing a n g l e w a s assumed to be 30 and the coefficient of earth pressure at rest, KQ, was calculated from an equation given by Blshop(1958). K Q = 1. - sin d)' (5-1) The i n i t i a l mean effective stress <r£ was calculated from effective vertical stress, <rv , by, <T0 = ^.(l.+2.K0) .<4 (5-2) Through the use of equations (3-5) and (3-6) the i n i t i a l maximum shear modulus, G m o, and the i n i t i a l maximum shear s t r e s s , t m Q , as shown in columns 7 and 8 of Table 5-1 were obtained. Column 9 gives the masses of the lumped mass system calculated based on a unit cross-section s o i l column while column 10 gives the i n i t i a l spring stiffnesses. Choice of damping values is discussed in section 5*2. 5.2 RESPONSE CHARACTERISTICS OF THE SOIL DEPOSIT. In order to investigate the steady-state response of the s o i l deposit to harmonic base motion through the use of the proposed non-linear analysis the volume change parameter ij> used in equation (4-29) was set to zero. The reason for this i s obvious since transient response is included in the step-by-step integration technique, steady-state vibration TABLE 5-1 • VALUES FOR LATER APPROXIMATION TO SOIL DEPOSIT LATER THICK- DEPTH EFF VERT D G T _ m. K NO. NESS TO CTR STRESS r mo mo i 1 0 f t . f t . psf % kip/ft psf slug kip/ft 14 5.0 2.5 305.0 50. 594. 85. 9.5 118.8 13 5.0 7-5 760.0 50. 938. 212. 19.0 187.6 12 10.0 15.0 1210.0 50. 1184. 338. 28.4 118.4 11 10.0 25.0 1810.0 50. 1448. 506. 37.9 144.8 10 10.0 35.0 2410.0 55. 1733. 674. 37.9 173.3 9 10.0 45.0 3010.0 60. 2010. 841. 37.9 201.0 8 10.0 55.0 3610.0 65. 2283. 1009. 37.9 228.3 7 10.0 65.0 4210.0 70. 2557. 1177. 37.9 255.7 6 20.0 80.0 5110.0 75. 2923. 1428. 56.9 146.1 5 20.0 100.0 6310.0 80. 3369. 1764. 75.8 190.6 4 20.0 120.0 7510.0 85. 3813. 2099. 75.8 190.6 3 20.0 HO.O 8710.0 85. 4106. 2434. 75.8 205.3 2 20.0 160.0 9910.0 85. 4380. 2770. 75.8 219.0 1 30.0 185.0 11410.0 85. 4700. 3189. 94.8 156.7 119 occurs only after a few cycles of vibration and i f hardening or softening, i.e. change in maximum shear modulus, G m n» and maximum shear stress, T^, is allowed for every cycle the steady state would never be achieved. Moreover, to minimize as far as possible the duration of the transient vibration, the amplitude of base acceleration input was made to increase gradually in the f i r s t few cycles to the required level and kept constant, thereafter. After a few t r i a l s i t was found that i f the amplitude of base acceleration i s made to rise as a sine function within the f i r s t two and half cycles transient vibration has i t s least influence. The response of non-linear system depends not only on the frequency of the input loading but also on the magnitude. Consequently, complete response characteristics requires a lot of computer runs. In view of this, an a r b i t r a r i l y fixed base acceleration level of 0.065g was used. A time acceleration plot for the base input is as shown in Figure 5-2. Even though the hysteretic s o i l model contributes quite large damping, e.g. 20 to 30% in terms of damping ratio defined by equation (3-30), at strain levels around 0.10 per cent, viscous damping is s t i l l needed to take into account of the effect of the flow of water inside the s o i l structure. In addition, due to the stepwise-llnear nature of the integration procedure a small amount of viscous damping is needed to eliminate high frequency vibration at the reversal points. Small amount of viscous damping was added through f n Time (sees) FIG. 5 -2 Input Base Acceleration 121 defined in equation (4-2). Also i t was assumed that linear viscous damping can be used and that the value of the viscous damping coefficients, c i t are proportional to the i n i t i a l stiffness values, k ± t=0* namely, c i 8 8 jS-ki,t=0 (5-3) When a ft value of 0.0008 was used surface acceleration response as shown in Figure 5-3& was obtained while p =0.002 gave a response as in Figure 5-3b. It should be noted that the amplitude of response is slightly changed but the high frequency vibrations at reversal points are eliminated. By calculating the energy dissipated per cycle for the viscous force as given by equation (4-24) i t can be shown that a /$ value of 0.0008 gives approximately a damping ratio of 0.6% while 0.002 gives 2 %. For a l l analyses described in this chapter a 5^ value of 0.002 was used. Step-by-step integration procedures were carried out for at least 10 cycles of vibration and i t was found that in most cases steady-state response could be reached. Figure 5-3b shows a typical time history response of the surface acceleration. A typical stress-strain response of a layer is as shown in Figure 5-4. Twelve frequencies were used and the ratio of the amplitude of the steady-state surface acceleration to that of the base acceleration Input was computed for each frequency. The result obtained is plotted in Figure 5-5. For comparison purposes, the "SHAKE" program 122 c u u o> in c o Q> U U < Time (sees) u - 2. c o ^ 0 o v. 0) O - Z u -4. f\ VJ VI /I /I VI VI /I /I VI VJ /I _L 2. 3. Time (sees) 4. FIG. 5-3 Surface Acceleration Response 123 .RESPONSE OF HORIZONTAL SOIL DEPOSIT TO HARMONIC BASE ACC. STRESS-STRAIN FOR LAYER 12 0 . 0 4 0 . 0 2 0 0 . 0 2 0 . 0 4 Shear Strain (%) FIG. 5 - 4 Stress - Strain Response Steody State Surface Acceleration Base Acceleration 125 written by Schnabel et al(1972) was also used and the result is plotted in the same figure. It can be seen that while the equivalent linear analysis used in the "SHAKE" program gives a pronounced peak at the resonant frequency, the non-linear analysis gives a much less pronounced peak in the frequency-response curve. A plausible explanation is that when both the equivalent linear and the non-linear models are subjected to the same strain level, the former has a larger amount of strain energy which is responsible for the peaked response curve. Apart from this, the two methods give reasonably close acceleration ratio within the frequency range of 2 to 10 hertz. A typical surface acceleration response calculated by the "SHAKE" program is shown in Figure 5-6. Comparing with Figure 5-3, i t can be observed that the equivalent linear analysis gives the usual sinusoidal surface acceleration while in the non-linear case the acceleration response is no longer sinusoidal. Jennings(1964) used a numerical integration method to obtain the response of a single-degree-of-freedom system with non-linear hysteretic force-deflection relationship and observed that the acceleration response was not a sinusoidal one. This is in agreement with the result obtained by the proposed non-linear analysis. It is interesting to note from Figure 5-3 and Figure 5-6 that both the response computed by the 'SHAKE' program and that computed by the proposed non-linear analysis have more or less the same 127 phase shift with respect to the input acceleration wave form. 5.3 MAGNITUDE OF PORE PRESSURE RISE AND ITS EFFECT ON THE DYNAMIC RESPONSE OF THE DEPOSIT. In this section a simulated earthquake record was used as the base input to the deposit described in the previous section. The acceleration values were calculated from a random acceleration function suggested by Bogdanoff, Goldberg and Bernard(I96I), namely, J x g(t) = • .1. jjt.exp(-5jt).cos(wjt+ej) t>0 (5-4) where j^y ^ are real positive numbers, w l < w2 < w3 < w4 ' • • » W J » a n d » S^i 62t Ojt • • » t 0j are real random variables uniformly distributed over the interval 0 to 27T. A large number of terms can be used when best representation of an accelerogram is needed. In a study of the seismic response of structure-foundation systems Parmelee et al(l968) followed the suggestion of Bogdanoff et al(196l) and assumed and to be constant when generating a simulated accelerogram. Ten terms were used and the acceleration function used was given by, 10 x g(t) = O.50.t.exp(-0.333t).£1cos(wjt+elj) t>0 (5-5) in which the frequencies, Wj, and phase angles, 8-j, are shown in Table 5-2. In this section, the same accelerogram scaled 128 to various values of maximum acceleration was used for - a l l analyses. Figure 5-7 shows a plot of the accelerogram. It can be seen that after 15 seconds the amplitude of motion is relatively small so that the analyses can be stopped at this time without losing any significant data. For the I n i t i a l computer run, the volume change parameter, was s t i l l kept Identically zero so that a comparison could be made with results obtained otherwise. Figure 5-8a shows the surface acceleration obtained by the non-linear analysis while Figure 5-8b shows the same obtained by using the computer program "SHAKE". In both cases the acceleration values of the base motion were scaled so that the maximum acceleration was 0.065g. Though direct point to point comparison of the surface acceleration responses is not possible the same general characteristics seem to be present in both responses. Thus, so far, i t can be seen when there is no material hardening due to compaction nor material softening due to rise in pore water pressure, both the method of direct integration of the non-linear dynamical equations and the method of wave propagation with strain-compatible material constants give more or less the same response to harmonic as well as to random excitation. Figure 5-9a and b show the typical strain and stress response of a layer obtained by the non-linear analysis. In order that the pore pressure generated during shaking might be accounted for the respective «J/. values for FIG. 5~7 Simulated Earthquake Accelerogram TABLE 5-2 FREQUENCIES AND PHASE ANGLES FOR THE SIMULATED EARTHQUAKE ACCELEROGRAM 3 UK (rad/sec) 8. (radians) «J 1 6.00 3.7663 2 8.00 1.3422 3 10.00 4.8253 4 11.15 0.2528 5 12.30 4.5204 6> 13.25 1.8834 7 14.15 1.3320 8 16.20 1.7852 9 17.35 0.1517 10 19.15 2.4881 TABLE 5-3 %-VALUES FOR DIFFERENT LAYERS Layer No. Relative Density it) +1 14 50. 0.0000 13 50. 0.5200 12 50. 0.5200 11 50. 0.5200 10 55. 0.4727 9 60. 0.4333 8 65. 0,4000 7 70. 0.3714 6 75. 0.3467 5 80. 0.3250 4 85. 0.3059 3 85. 0.3059 2 85. 0.3059 1 85. 0.3059 131 Time (sec.) FIG. 5-8 Surface Acceleration Response (no Pore Pressure ) 132 each layer were calculated based on the Inverse proportional relationship between 4^ a n d relative density, Dr, as discussed in section 3.3.2 of chapter 3. Assuming that the present s o i l deposit has the volume change characteristics given by Figure 3-11, the vj/^ values listed in Table 5-3 are obtained. Since the top layer (layer 14) is above the water table i t s volume change parameter \\>^ was set to zero so that excess pore water pressure due to compaction could not be generated. As before, the acceleration values were scaled down so that a maximum of 0.065g was obtained. In this analysis the layers were assumed to be separated by impervious boundaries and re-distribution of pore pressure was not allowed. As previously discussed the integration procedures were carried out up to 15 seconds. Figure 5-10a and b show "the surface acceleration and displacement responses and Figure 5-Ha, b and c show the strain, stress and pore water pressure of the 11th layer which is situated between depths of 20.0 to 30.0 feet. Figure 5-12 shows a stress-strain plot for the 11th layer. The time at which liquefaction occurred was around 8.5 second and i t is marked in a l l the response plots. From Figure 5-8a and 5-10a i t can be seen that as the pore pressure rises during the liquefaction process surface acceleration tends to be attenuated. The relatively large acceleration response that occurs at a time of 9 seconds when pore pressure is not generated is not seen in Figure 5-10a where pore water pressure is taken into account. As 133 FIG. 5-9 Strain and Stress Responses for Layer II (no Pore Pressure) FIG. 5-10 Surface Acceleration and Displacement Response (Pore Pressure Generated) c to 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 — . T — [ ..., , I , , , — i 1 1 r i i •A J\ * ~ A \ A A——-\ / 1 \ / v ~ 1 1 I I 1 1 1 I — (a) : 1 1 1 1 1 1 400.h T 1 1 1 1 1 r ^2000. a. CO CO £ 1000. Q_ <P i— O 0. 1 1 1 1 r (c) 5. 10. Time (sec) T — | — i — i — i — r I I L T 1 1 1 1 T r 15. FIG. 5-11 Stress, Strain and Pore Pressure Response for Layer 11 FIG. 5-12 Stress - Strain Response for Layer II (Pore Pressure Generated ) 137 expected the acceleration response within the f i r s t three seconds was not much affected by the Increase in pore pressure which as shown in Figure 5-13 did not rise to a high value within this time interval. Figures 5-9a and 5-Ha shows the dramatic effect of the pore pressure rise on the strain response of the 11th layer. In the case where no pore pressure was generated the peak strain response was +0.05$ at around 2 seconds. The same peak was also obtained in the case where pore pressure was generated. In the latter case, however, the +0.05$ strain response at around 2 seconds was dwarfed by the +0.25$ strain at 8 seconds, a short time before liquefaction, as shown in Figure 5-lla. The softening behaviour of the s o i l ls well illustrated by the increasing strain response even when the input base acceleration as well as the stress response are decreasing. Moreover, the large difference in the maximum strain levels developed in the two cases casts some doubt on the validi t y of using an equivalent linear analysis for the study of s o i l liquefaction as was done by Seed and I d r l s s ( I 9 6 7 ) . The offset of the strain response of the 11th layer, and hence the surface displacement response, to one side of the zero axis after liquefaction is probably due to the slight non-symmetry of the input base motion. Figure 5-12 shows the stress-strain response of the 11th layer and i t can be seen that even after liquefaction the stress and strain were s t i l l 138 describing the kind of loops which they were supposed to describe. Close examination of the stress-strain plot reveals the slight one-sidedness of the stress response after 9 seconds. The reason is that at such low level of confining stress the failure stress is very small so that a small stress increment near the.failure stress can induce a very large strain. Figure 5-13 shows the distribution of pore pressure in the 200-foot deep s o i l deposit. I n i t i a l l y , the excess pore water pressure was an increasing function of depth. However, as time went on the s o i l layers situated between 10 and. .50 feet depth interval began to develop higher excess pressure. Consequently, these layers developed larger strains at later times which led to even higher excess pore water pressures. The interaction seemed to go on u n t i l liquefaction occurred in the 11th layer which is situated between 20 and 30 feet. 5.4 EFFECT OF PORE PRESSURE DISSIPATION. To investigate the effect of pore pressure dissipation on the liquefaction process the complete coupled analysis, i.e. including dyanmic response, pore water pressure generation and dissipation, was performed on the s o i l deposit described earlier. The simulated earthquake input described in section 5*3 of this chapter was again used. Analyses were carried out using three different values of the permeability, Pore Pressure (psf) (Note • values for I. sec. curve should be divided by I0<) FIG. 5-13 Pore Pressure Distribution at Various Times (no Dissipation) 140 k. In each analysis, the value of the permeability was assumed to be constant for the entire depth of the deposit. Even though the computer program was set up to accept permeability values that vary from layer to layer the capability was not used since the purpose was to obtain a general idea on the influence of pore pressure dissipation. As a matter of fact, such investigation is very l i k e l y the f i r s t of i t s kind carried out to date(1975). Terzaghi and Peck(l96?) l i s t e d a range of permeability values for various types of s o i l that are encountered in engineering practice. It was decided to analyse three different cases with permeability values of 0.000003, 0.0003 and 0.03 feet per second respectively and these approximately correspond to nearly impervious fine sands, medium permeable sands and very permeable clean sands. From here onwards these three cases w i l l be referred to as case Dl, D2 and D3 respectively. Pore pressure distribution for these cases are shown in Figure 5-l4a, b and c respectively. Without pore pressure dissipation, s o i l layer 11 liquefied at around 8 second as described in section 5«3 of this chapter. When dissipation of pore water pressure was allowed but with a very low value of permeability, as in case Dl, liquefaction occurred at the 12th layer, i.e. the zone of liquefaction was shifted upwards. This behaviour agrees with the boundary conditions of the problem since Pore Pressure (psf) 0 IOOO. 2000. 100 I50.h 200 (a) Pore Pressure (psf) 0 1000. 2000. 50. 100. I50.h 200. \ \ • • • I I / • • • . / / / / / / • • • 1/ / // • • • 3 CD 01 CO CD O CO (b) Pore Pressure (psf) .0 1000. 2000. 200.6 (c) FIG. 5-14 Pore Pressure Distribution for Different k Values 142 drainage can only take place at the top surface. When the permeability value of the s o i l was made larger, as in case D2, liquefaction did not occur though similar pore pressure distributions occurred at much later time, as shown by the 8-second and 15-second curves in Figures 5-l4a and b respectively. When the permeability value was made 0.03 feet per second the pore pressure for s o i l strata situated between depths 0 and 100 feet actually decreased during the interval of time 8 to 15 seconds. In general terms the following statements could be made in relation to the particular s o i l deposit analysed. 1. Relatively low permeability values (around 0.000003 feet per second) tends to shift the zone of liquefaction upwards and dissipation of pore pressure does not improve the resistance to liquefaction failure of the s o i l deposit to any significant extent. 2. Medium to large permeability value (0.0003 to 0.03 feet per second) Improves the resistance to liquefaction failure by dissipating the high pore pressure generated. Pore pressure may decrease instead increase after the peak disturbance is passed even though there are t r a i l i n g weak disturbances fed to the deposit. It may be visualized that i f the earthquake motion consists of strong pulses separated by weak ones the maximum pore pressure near the surface would be even smaller. 3. Within the time interval of 15 seconds the magnitude and 143 distribution of pore pressure below the 100-foot depth is l i t t l e affected, i f any, by the dissipation process. If the earthquake waves consist of a relatively low-valued and uniform amplitude acceleration waves with long duration, say 60 to 80 seconds, liquefaction of deeper zones is possible. Quantitative results obtained in this chapter are s t r i c t l y applicable to the particular deposit described in section 5*1. Qualitative extrapolation to other deposits in terms of the effects of some of the s o i l parameters on the liquefaction process can be made. In addition, the numerical solution of the coupled problem Is shown to be practicable as well as feasible. Computer time for numerical integration for any one case described in this section, i.e. integrating for a time history of 15 seconds at 0.01 second interval, is approximately 15 seconds CPU time. 144 CHAPTER 6 SUMMARY AND CONCLUSIONS / 6.1 SUMMARY The problem of s o i l liquefaction has been a main concern in the design of foundations for earthquake-resistant structures on saturated granular s o i l s . Previous investigations indicated that for horizontal s o i l deposits the important factor is the, amount of pore water pressure accumulated during the period of shaking. In this thesis i t is assumed that for the saturated granular s o i l comprising the deposit, liquefaction occurs when the effective normal stress reaches a value of 0.05 of i t s i n i t i a l value. A material model is proposed whereby such changes in pore water pressure can be calculated by following the strain history of a s o i l element subjected to shear deformation. The s o i l deposit is assumed to be of inf i n i t e extent and f i n i t e thichness resting on a firm bedrock in order that simplification of the problem to one of the shear-beam type is possible. Equations are derived to describe the coupled processes involved during the progress of liquefaction, namely, 1. Dynamic response to base motion, 2. Increase in volumetric compaction strain due to cyclic straining, 145 3. Pore water pressure generation as a result of (2), and, 4. Dissipation of excess pore water pressure(consolidation). Solution of the problem was achieved by using numerical integration. The f i n i t e difference method was used to arrive at an approximation of the continuous s o i l deposit by a lumped-mass-system in a rational way. A step-by-step integration method was then used to solve the resulting dynamical equation. A separate f i n i t e difference scheme was used to solve the consolidation equation. Several example problems were solved for an assumed deposit by varying some of the s o i l parameters so that factors affecting the liquefaction process can be studied. 6.2 CONCLUSIONS Conclusions reached in this research are as follows:-1. The hyperbolic hysteretic stress-strain relationship with volumetric compaction 'strain serving as a 'hardening' parameter, can be used effectively to describe the liquefaction behavior of s o i l under dynamic loading. 2. Using the proposed s o i l model for a deposit a relatively • f l a t ' frequency response is obtained. This indicates that the resonant capacity of the non-linear model is lower than i t s linear equivalent. 3. When pore water pressure accumulated during shaking reaches 146 A value of about 90.% of the i n i t i a l effective normal stress, very large strains develop In comparison with those computed for the same deposit but without taking into account of pore pressure r i s e . This throws some doubt on the validit y of using strain compatible equivalent linear analysis to obtain information concerning liquefaction potential. For relatively low permeability values (0.000003 to 0.0003 feet per second) in the deposit a localized type of liquefaction is produced so that stiffness and strength deterioration due to pore pressure rise for a particular layer occur at a much faster rate than others. In the usual case when the relative densities of the surface strata are lower than those of the deeper strata the zone of liquefaction occurs near the surface, i.e. at depth between 20.0 to 50.0 feet. Relatively large permeability values (around 0.03 feet per second) reduce the tendency to localize stiffness and strength deterioration to a single layer. For strata near the surface the maximum pore water pressure accumulated during the time of earthquake distrubance, is reduced and larger magnitude base motion is require to cause liquefaction. In general, the ratio of the pore pressure accumulated to the i n i t i a l effective normal stress decreases with depth after the 50.0 feet mark. For s o l i deposits without especially weak layers liquefaction at depth greater than 14? 100.0 feet is very unlikely. For s o i l stratra located 50.0 feet or more beneath the ground surface, the pore pressure accumulated is not much affected by the dissipation process even when the permeability value is as high as 0.03 feet per second. Therefore, when a permeable s o i l deposit is shaken by earthquake motions consisting of relatively low, uniform amplitude and long duration acceleration waves, liquefaction at larger depths ( i . e . 50.0 and 100.0 feet) is possible. 3 SUGGESTIONS FOR FURTHER RESEARCH In view of the parti a l stabilization effect produced by sandy material due to i t s dllatant behavior at low confining pressures, incorporation of dilatancy into the stress-strain relationship i s needed so that the analysis may be continued near and at liquefaction (zero or near zero effective stress). A thorough investigation of the volume compaction characteristics is necessary in order that the volume compaction function can be generalised to Include, 1. a continuous variation of volumetric compaction as a function of strain or stress, and, 2. effect of a general state of stress. A correlation study with f i e l d data should be carried out. 148 LIST OF REFERENCES Ambraseys, N., and Sarma, S., I969. "Liquefaction of Soils Induced by Earthquakes." Bull, of the Seism. Soc. of Amer., vol. 59» no. 2, April 1969, PP. 651-664. Barkan, D.D., I 9 6 2 , "Dynamics of Bases and Foundation." McGraw H i l l Book Co., Inc., I962. Bathe, K.J., and Wilson, E.L., I973t "Stability and Accuracy Analysis of Direct Integration Methods." Internl. J. Earthquake Engg. and Struct. Dyn., vol. 1, 1973. PP. 283-291. Bazant, Z., and Dvorak, A., 1965. "Effects of Vibrations on Sand and the Measurement of Dynamic Properties." Proc. 6th Internl. Conf. Soil Mech. & Found. En*rg., vol 1, 1965. PP. 161-164. Blot, M.A., 1956, "Theory of Propagation of Elastic Waves in Fluid-Saturated Porous Solid. I - Low-Frequency Range." J . Acoustical Soc. of Amer., vol. 28, no. 2, March 1956, pp. 168-178. -Bogdanoff, J.L., Goldberg, J.E., and Bernard, M.L., I 9 6 I , "Response of a Simple Structure to a Random Earthquake Type Disturbance." Bull. Seism. Soc. Amer., vol. 511 no. 2, April 1961, PP. 293-310. Casagrande, A., 1936, "Characteristics of Cohesionless Soils Affecting the Stability of Slopes and Earthf i l l s . 1 1 J . Boston Soc. C i v i l Engrs., Jan. 1936. 149 Castro, G., 19^9. "Liquefaction of Sands." Ph.D. Thesis, Harvard Univ., Cam., Mass., I969. E d i t o r i a l Committee of the "General Report.", I968, "General Report on the Niigata Earthquake of 1964." Finn, W.D.L., 1972, "Soil Dynamics - Liquefaction of Sands." Proc. Internl. Conf. Microzonation for Safer Construction Research arid Application, Seattle, Washington, USA, Oct., 1972, pp. 87-111. Finn, W.D.L., Emery, J.J., and Gupta, Y.P., 1970, "A Shake Table Study of the Liquefaction of Saturated Sands during Earthquakes.", Proc. 3rd European Sym. Earthquake Engg., Sept. 1970, PP. 253-262. Finn, W.D.L., Emery, J.J., and Gupta, Y.P., 1971, "Liquefaction of Large Samples of Saturated Sand on a Shake Table.", Proc. 1st Can. Conf. Earthquake Engg., Van., UBC, May 1971, PP. 97-110. Finn, W.D.L., Emery, J.J., and Gupta, Y.P., 1971a, "Soil Liquefaction Studies Using a Shake Table." Closed Loop, System Corporation, Fall/Winter, 1971. Finn, W.D.L., Pickering, D.J., and Bransby, P.L., 1971. "Sand Liquefaction in Tri a x i a l and Simple Shear Tests." J. S o i l Mech. & Found. Div. t ASCE, vol. 97, SM 4, April 1971. PP. 639-659. Florin, V.A., and Ivanov, P.L., 1961, "Liquefaction of Saturated Sandy Soils." Proc. 5th Internl. Conf. on Soil Mech. and Found. Engg., vol. 1, 1961, pp. 107-111. Forssblad, L., I967, "New Method for Laboratory Soil Compaction by Vibration." Highway Research Board, no. 177, 1967, pp. 219-224. 150 Hardin,. B.O., and Dmevich, V.P., 1972, "Shear Modulus and Damping in Soils : Measurement and Parameter Effects." J. Soi l Mech. & Found. Div., ASCE, vol. 98, SM 6, June 1972, pp. 603-624. Hardin, B.O., and Drnevich, V.P., 1972a, "Shear Modulus and Damping in Soils : Design Equations and Curve s." J. Soi l Mech. & Founds. Div., ASCE, vol. 98, SM 7, July 1972, pp. 667-692. Hazen, A., 1920, "Hydraulic F i l l e d Dams." ASCE Transaction, vol. 83, 1920, pp. I9I7-I945. Herrera, I., I 9 6 5 , "Modelo Dinamicos para Materiales y Estructuras del Tipo Masing." 'Boletin Sociedad Mexicana de Ingenieria Sismica, 3(1), 1965, PP. 1-8. Ishibashi, I., and Sherif, M.A., 1974, "Soil Liquefaction by Torshional Simple Shear Device." J. Geot. Engg. Div., ASCE, vol. 100, no. GT 7, July, 1974, pp. 871-904. Ishihara, K., and L i , S.L., 1972, "Liquefaction of Saturated Sand in Tria x i a l Torsion Shear Test." Soils and Foundation, vol. 12, no. 3. June 1972, pp. I939. Iwan, W.D., I 9 6 6 , "A Distributed-Element Model for Hysteresis and i t s Steady-State Dynamic Response." J. Appl. Mech., vol. 33, no. 4, I966, pp. 893-900. Jennings, P.L., 1964, "Periodic Response of a General Yielding Structure." J. Engg. Mech. Div., ASCE, vol. 90, EM~ 2, April 1964 pp. 131-195. 151 Kishlda, H., 1964, "Damage to Reinforced Concrete Buildings in Niigata City with Special Reference to Foundation Engineering." So i l and Foundation, vol. 7, no. 1, March 1964, pp.71-88. Kishida, H., 1969, "Characteristics of Liquefied Sands during Mino-Owari, Tohnankal and Fukui Earthquakes." Soils and Foundations, vol. 9, no. 1, March 1969t PP. 75-92. Kishlda, H., 1970, "Characteristics of Liquefaction of Level Sandy Ground during the.Tokachioki Earthquake." Soils and Foundations, vol. 10, no. 2, June 1970, pp. 103-111. Koizumi, Y., 1966, "Changes in Density of Sand Subsoil Caused by the Niigata Earthquake." Soil and Foundation, vol. 6, no. 2, June 1966, pp. 38-44. Kondner, R.L., and Zelasko, J.S., I963. "A Hyperbolic Stress-Strain Formulation for Sands." Proc. 2nd Pan Amer. Conf. Soil Mech. & Found. Engg., 1963, pp. 289-324. Marsal, R.J., 1963. Internal Report, Institute of Engineering, National University of Mexico, Mexico, 1963. Martin, G.R., Finn, W.D.L., and Seed, H.B., 1975, "Some Fundamental Aspects of Liquefaction under Cyclic Loading." J. Geot. Engg. Div., ASCE, vol. 101, no. GT 5, May 1975, PP. 423-438. Maslov, N.N., 1957, "Questions of Seismic Stability of Submerged Sandy Foundations and Structures." Proc. 4th Internl. Conf. Soil Mech. & Found. Engg., vol. 1, 1957, PP. 367-372. 152 Maslng, G., 1926, "Eigenspannlngen und Verfestigung beim Messing." Proc. 2nd Internl Cong. Applied Mech., 1926, PP. 332-335. Mogami, R., and Kubo, 1953, "The Behaviour of Soil during Vibration." Proc. 5th Internl. Conf. Soil Mech. and Found. Engg., vol. 1, Zurich, 1953, PP. 182-185. Newmark, N.M., 1959 "A Method of Computation for Structural Dynamics." J. Engg. Mech. Div., ASCE, vol. 85, no. EM3, July 1959. Newmark,.N.M., and Rosenblueth, E., 1971, "Fundamentals of Earthquake Engineering." Prentic Hall, Inc., Eaglewood C l i f f s , N.J., USA, 1971. PP. 162-163. Osaki, Y., 1968, "Building Damage and Soil Condition." General Report on the Niigata Earthquake of 1964. Committe of "General Report", I968, pp. 355-383. Parmelee, R.A., Perelman, D.S., Lee, S.L., and Keer, L.M., 1968, "Seismic Response of Structure Foundation System." J. Engg. Mech. Div., ASCE, vol. 94, EM 6, Dec. 1968, pp. 1295-1315. Peacock, W.H., and Seed, H.B., I968, "Sand Liquefaction under Cyclic Loading Shear Conditions." J. S o i l Mech. & Found. Div., ASCE, vol. 94, SM 3, May 1968, pp. 689-708. Pyke, R.M., 1973. "Settlement and Liquefaction of Sands under Multi-Directional Loading." Ph.D. Thesis, Univ. of Calif., Berkeley, Calif., 1973, P. 280. 153 Schnabel, P.B., Lysmer, J., and Seed, H.B., 1972, "A Computer Program for Earthquake Response Analysis of Horizontally Layered Sites." Earthquake Engineering Research Centre, Report Number EERC 72-12, Dec. 1972. Seed, H.B., I968, "Landslides during Earthquakes due to Liquefaction." J. So i l Mech. & Found. Div., ASCE, vol. 94, SM 5, Sept. 1968, pp. 1053-1122. Seed, H.B., 1972, "Dams and Soils." The San Fernando Earthquake of February 9, 1971, and Public Policy. Special Subcommittee of the Joint Committee on Seismic Safety, California Legislature, July 1972. Seed, H.B., and Idriss, I.M., I 9 6 7 , "Analysis of Soil Liquefaction : Niigata Earthquake." J. Soi l Mech. & Found. Div., ASCE, vol. 93, SM 3, May 1967, PP. 83-108. Seed, H.B., and Idriss, I.M., 19^9, "Influence of Soil Conditions on Ground Motions during Earthquakes." J. Soi l Mech. & Found. Div., ASCE, vol. 95, SM 1, Jan. I969, PP. 99-137. Seed, H.B., and Idriss, I.M., 1971, "Simplified Procedures for Evaluating S o i l Liquefaction Potential." J. So i l Mech. & Found. Div., ASCE, vol. 97, SM 9, Sept. 1971, pp. 1249-1273. Seed, H.B., and Lee, K.L., I 9 6 6 , "Liquefaction of Saturated Sands during Cyclic Loading." J. S o i l Mech. & Found. Div., ASCE, vol. 92, SM 6, Nov. 1966, pp. 105-134. Seed, H.B., and Silver, M.L., 1972, "Settlement of Dry Sands during Earthquakes." J. Soi l Mech. & Found. Div., ASCE, vol. 98, SM 4, April 1972, pp. 381-397. 154 Silver, M.L., and Seed, H.B., 1971. "Volume Changes in Sands during Cyclic Loading." J. Soil Mech. & Found. Div., ASCE, vol. 97. SM 9. Sept. 1971, PP. 1171-1182. Sneddon, I.N., 1957. "Elements of Partial Differential Equations." McGraw H i l l Book Company, Inc., 1957. State Earthquake Investigation Commission, 1908, "The California Earthquake of April 18, 1906." Report of the State Earthquake Investigation Commission, Published by the Carnegie Institute of Washington, 1908. Terzaghi, K., and Peck, R.B., I968, "Soil Mechanics in Engineering Practices." John Wiley Publishing Company, First Corrected Printing, 1968. Yen, B.C., I967, "Viscosity of Saturated Sand near Liquefaction." Proc. Internl. Sym. Wave Propag. & Dyn. Properties of Earth Materials, New Mexico, I967. PP. 677-888. Yoshimi, Y., 1967. "An Experimental Study of Liquefaction of Saturated Sands." Soli and Foundation, vol. 7. no. 2, June 1967, PP. 20-32. 155 APPENDIX. I LUMPED MASS SYSTEM APPROXIMATION FOR NON-LINEAR SHEAR BEAMS The equilibrium equation, equation (4-1), obtained in chapter 4 is re-written here for convenience, s$ - 8 = 0 When the s o i l parameters are functions of depth, (H-z), where H i s the total depth of the s o i l deposit and z i s the height from the surface of the bedrock, equation (4-2) has to be replaced by a more general equation, namely, T = f^JT.z) + f 2 ( f , z ) (1-2) and also, f = j>(z) (1-3) Equations ( 1 - 2 ) and (1-3) can be substituted into (1-1) to give, Consider, for example, the z-domain i s sub-divided into intervals, so that O ^ ^ z ^ z ^ . , . <zn=H , and l e t the length of each interval be denoted by h^, i.e., 156 Also, l e t x., f. and p. represent the values of x, t and 1 1 / 1 c j> at ,. and denote the particular function at z^ by f 2 , i ' i , e ' f2,i = f 2 ( fi' zi> ^- 6) The central difference quotient for Tz^2^'z^ a ^ z i ^ s given by, Af2(irH (f2.1+i " f 2 , i - i > / < i h i + l + * h i > and a similar expression can be written for f^(f,z) so that when these expressions are substituted into (1-3)• the following equation can be obtained, i(h i + 1+h i).^ 1.x i + ( f ^ . . ^ - f 1 § i + i ) + <f2,i-i " f2,i+*> = 0 tt-8) Now, write, tM,z) tz(f,z) = -L- . t (1-9) Since 3x dz fi-Jr ( x i * x i - l ) / h i ( I - 1 0 ) Consequently, equations (1-9) and (I-10) give us, f o ( f ? l ) 1 1 - 2 I f a new symbol i s defined as follows :-i h. . f. I 157 (1-12) equation ( I-11) can be written i n a simple form, f 2 , i - i *' k i ' ( x i " x i - i } (1-13) S i m i l a r l y , i t can be shown that, f 2.1+4 F ' f i . i + i ,=' k i + l * ( x i + l " x i ^ ^ i ' ( x i ~ x i - l ^ c i + l * ( x i + l x i ^ (1-14) (1-15) (1-16) where, c i A . Z . i ) 1 1 - 2 (1-17) Equations (1-13) through ( I - l 6 ) can be substituted into equation (1-8) and the following 'di s c r e t i z e d • dynamic equation i s obtained, V x i + (-^ c i + c i + 1 - c i + 1 ) . ^ X i - 1 X.. > X i + 1 + ( " k i k i + k i + l - k i + l ) * i - l X i x i + l * = 0 (1-18) Note that, i ( h i + 1 + h i ) . j > i 158 and since the average density within the interval i s usually known instead of the density value at z^, m^ i s "better calculated by, V = * \ f i - i - h i + ^ / i + i ^ i + l ( I " 1 9 ) It should be pointed out that the f i n i t e difference equation (1-18) holds in general for i=2,3,4,... ,n-l. For i=l and i=n, i.e. at both end points, the same equation can be used by observing that when i=l the variables x Q and x Q appearing in equation (1-18) are in fact the velocity and displacement of the base respectively so that the symbols x^ and x^ are used instead. And, when i=n the variables with the subscript of n+l are set identically equal to zero. A system of n equations thus obtained can be represented by a simple matrix equation!-[M].{'x} + [C].{x} + [K].{x> = {Pb> (1-20) where {x}, {xj- and {x} are column vectors whose components x i ' x i a n d x i Sive us the acceleration, velocity and displacement respectively of the point at z^, {x} = (x x x £ x^ . . . x n ) T {x} = (x x x 2 x^ . . . x n ) T {x} = (x x x 2 x 3 . . . x n ) T {Pb} = (k 1x b+c 1x b 0 0 . . . 0 ) T 159 where the superscript T r e f e r s to the transpose of the row vector into a column vector and x^ and x^ r e f e r to the base v e l o c i t y and displacement r e s p e c t i v e l y . The matrices [M], [C] and [K] as obtained from equation (1-20) are as follows / [M] = m1 0 0 m2 0 0 • • 0 ' 0 0 0 m^ , • 0 (I -21) • • 0 0 • * 0 • • m n j [C] = c 1 + c 2 ~C2 0 0 • • 0 " C 2 C 2 + C 3 " c3 0 • • 0 0 - c 3 c 3 + c 4 • • 0 (I -22) • 0 * 0 • 0 • 0 • • • - C n • c n [K] = " k l + k 2 0 0 • • -0 " k 2 k 2 + k 3 -k 3 0 t • 0 0 -k 3 k 3+k 4 • • 0 ( I -23) • 0 • 0 • 0 • 0 • • • k n For base acceleration type of loading the r e l a t i v e displacement vector i s preferred. Let {x} be the r e l a t i v e displacement vector, {X) and {X} be the corresponding 160 velocity and acceleration vectors and l e t x^ be the variable describing the displacement of the base as a function of time, i t can be seen that {x} (1-24) Similar expressions can be written for the velocity and acceleration. When equation (1-24) i s substituted into equation (1-20), i t can be shown that the following equation i s obtained, [Ml{X) + [ C ] . { X ) + [K].{X} m n k x b (1-25) It i s quite apparent that the f i n i t e difference equation obtained i s the same as that for the dynamic response of a system of n discrete masses connected by springs and dashpots. Thus, a lumped mass model can be used to approximate the continuous shear beam model for the dynamic response of the s o i l column discussed in chapter 4, The value of the masses are obtained by lumping half the mass of each layer to i t s corresponding nodes, z^ ^ and z^. The stiffnesses of the springs and the damping coefficients for the dashpots are obtained from f 1 and f 2 of the stress-strain relationship, equation (4-2), and through the use of (1-12) and (I - l ? ) , Also, the width of the s o i l column can be made 161 d i f f e r e n t from unity and i t s e f f e c t can be taken into account by multiplying the corresponding width of layer, b^, int o the appropriate terras. So that, i n general, the components m^ ,. c^ and k^ of the matrices [M], [C] and [K] are given by the following i -m i = * « / > i - * ' V h i + £-j>i+£' bi+l-hi+l c i = tilh-i'H-tf'*! / ih-h-tf ( I - 2 6 ) k i = f 2 ^ i - i ' Z i - ^ - b i '/ <Vi-*> where i re f e r s to the layer number of the lumped mass system. 162 APPENDIX II NUMERICAL INTEGRATION OF THE DYNAMIC RESPONSE EQUATION The dynamic response equation for. the lumped mass system i s re-written here for convenience, [M].{X} + [C].{X} + [KJ.{X} = {P} (II-l) where {p} i s the inertia force vector, {p} -•m. (II-2) II.1 Incremental Stress-strain. Let T t + A t , ( H - 3 ) where t i s an instant of time at which a l l the quantities in equation (II-l) are known and At i s a small time increment. Since equation (II-l) has to he satisfied at every instant, we have, [M]T{X}T + [CJT{X}T + [K]T{X}T = (P}T [M] t{x} t + [C]t{x}t + [K]t{X}t = {P}t (II-5) (II-6) where the subscript refers to the instant of time at which the particular quantity takes on i t s value. Since the mass 163 matrix i s a constant matrix, [M] T = [M] t = [ M], so that [M] T{X} T - [MJ t{x} t = [M] (X} T - {X}t (II-7) However, the [Cj and [Kj matrices depend on the shear strain of the s o i l , as a result of the non-linear stress-strain relation, and thus on the displacement of the masses so that equations similar to that of (II-7) cannot be obtained for these matrices in a straight forward manner. There are two ways to arrive at expressions similar to (II-7) using approximation methods. A crude method i s to replace [K]^ by [K]^, so that, [K] T{X} T - [K] t{X} t = [K] t A better method i s as follows »-{x}T - {x} t 1 (II-8) Define, {F} [K].{X} ( H - 9 ) so that, [K] T{X} T - [K] t{X} t W t where [K], 3F; 5X jj {X}, = [ K] t.|{X} T -- W t ^X>t} 3F^ ax. J -(11-10) (11-11) 164 and i s the component of the column vector defined in (II-9 ) . Therefore, in general, [K] T{X} T - [K] t{X> t = LK] t.|{X} T - {X} t| (11-12) where [K]^ can take the value of [K]^ for crude approximation or i t can be evaluated according to (11-11) for better approximation. Similar procedures can be used to obtain IF equation (II-6) i s subtracted from equation (II-5) and taking note of (II-7) and (11-12), an incremental form for the dynamic response equation can be obtained as, [M]|{X}T - { X } t j + [ C ] t | { x } T - { X } t J + [Kj t[{x> T - { x > t J = {P}T - {P}T (11-13) It can be seen that the only unknown in this equation are the acceleration, velocity and displacement vectors at time T. As such, the numerical procedures proposed by Newmark(1959) can be applied to the equation so that {X}^, {x} T and {x} T can be expressed in term of {*X}t, {X}^ and {x} t. II.2 Step by Step Integration In Newmark's method of step by step integration two parameters, o< and |3 are used so that the velocity and displacement at time T can be expressed in terms of the acceleration, velocity and displacement at time t, and of the 165 unknown acceleration at time T. The relations can be written as follows :-{X} T = {X} t + (l-cx).At.{X}t + *.At.{X}T (11-14) {X} T = {X} t + At.{X}t + (i-/5).(At)?{X} t '+ j8.(At)?{x} T (11-15) Newmark( 1959) proposed that and /3=i for unconditionally-stable integration procedure, which incidentally corresponds to a constant average acceleration method of integration. If «=(l/2) and /6=(l/6) are used the linear acceleration method as proposed by Wilson and Clough(1952) i s obtained. For the present purpose °<=i and /3=^ are used. If the acceleration increment during the time interval At=T-t i s represented by {Q}^ ,. i.e., {Q) T = ( x } T - {x} t (11-16) and i f the following simplifying symbols are used, {a} t = At.{x} t (H-17) {b} t = At.{X}t + 2-.(At)?{X}t (11-18) then equations (II-14) and (11-15) can be written into a simple form, namely, {X} T-{X} t = {a) t + c*.At.{Q}T (II-19) {X}T {X>t = {b}t + y5.(At)?{Q}T (11-20) 166 Equations (II-16), (11-19) and (11-20) can be substituted into (11-13) to give a matrix equation involving only one unknown, {Q} , namely, [D] t{Q} T = {P}T (11-21) where [D] t = [M] + ot.At. [C] t + p. (At) 2.[K] t (11-22) and {P} T = { P } T - {P} t - [C] t{a} t - [K] t{b) t (11-23) Consequently, {Q}T = [Dj^P}^, (11-24) and {x}^, {x} T and {X) T can be calculated using equations (11-16), (11-19) and (11-20). In the numerical step by step integration, a typical sequence of calculations that have to be performed in one time step, At, i s as follows »-1. Based on {x} t > {X} t and {X} t at time t, (11-17) and (II-18) are evaluated, {a} t = At.{x} t {b}t = At,{x} t + i.( A t ) 2 . { x } t 2. Based on {x}^ and {x} t, {f}t and {f} t are calculated by, m t = CaJ{x}t (11-25) {i}t = [^]{x} t (11-26) where 167 [a] = l /h 1 0 0 0 - l / h 2 l /h 2 0 0 0 - l/h^ l/h^ 0 - l A l/h n n 3. With the {i\t and {t}t values, [ C ] t and [ K ] t are evaluated according to the appropriate expressions in appendix I. 4. [D]^. i s then calculated using equation (11-22). 5 . Using the base acceleration value at time T, x b T, i t i s possible to evaluate, r -\ T m n bT and ( P } T = {P} T - {P} t - [CJ t{a} t - [K] t{b} t 6. {Q}T i s solved by (11-24) so that the acceleration, velocity and displacement vectors, {X)T, {x)T and fX}T can be found by, {X}n {X)T = {X}t + {a} t +o<.At.{Q}T (X}T = {X}t + {b} t + |3.(At)2.{Q}T. 7. The steps 1 through 6 are repeated for the next time increment. 168 An option i s built into the analysis to bring the stress-strain point closer to the actual stress-strain curve. I f the strains and strain rates determined in each integration are accepted as true strains and true strain rates, the stress-strain and stress-strain-rate relationships can be used to determined the 'true' restoring force, [K] T.{X) T, and the 'true' damping force, t C ] T . {X} T. However, these restoring and damping forces do not necessarily satisfy the equilibrium equation (II-5). An a r t i f i c i a l 'external' force, { P . ™ }, defined by the following can be applied to the system so that equilibrium i s restored, {P e r r} = {P}T - [M].{X}T - [C] T.{X} T - [K] T.{X} T . The {P_„_ } can be added to {P}m in step number 5 for for the next time increment. It should be pointed out that this method of correction forces i s only good where the magnitude of {P-„_ } i s small as in the present case C X X of step-by-step integration with small time increment. 169 APPENDIX III FINITE DIFFERENCE EQUATIONS FOR PORE PRESSURE DISSIPATION In the general case where the permeability and compressibility of the s o i l are functions of depth equation (4-19) is used, *crw 3 k d<r„ a e ^ I F = K l r - ^ ^ ' 1 7 ) + K l r - " i T Following the same approach as in appendix I, i t is assumed that the z-domain is divided into n intervals so that, 0<z,<zo<zQ . . ,<z =H, l c. } n where H is the total depth of the deposit. For convenience the following notations w i l l be used :-hj to denote the length of the i-th interval, i.e., h i = Z i " Z i - 1 (IH -2) W to denote the f i n i t e difference approximation 1 ,T to at z=z^ and t=T, or, W i , T = °w< zf T> (IH-3) 170 Moreover, attention should be given to the following symbols mv 1 # k^ and ^ which are used to denote the respective values at the centre of the i-th interval. To indicate their values at the 'nodes* the symbols m (z ), k(z,) and v 1 i coefficient of volume €vp^ zi^ are used. compressibility so that, mv is the m 1/K l r (III-4) For the f i n i t e difference approximation of equation (III-l) consider f i r s t the term on the l e f t hand side of the equation and use a forward difference quotient for approximation so that, acrw St t=T J z=z, ( W i . T + A T - W i . T ) / A T (III-5) A central difference quotient w i l l be used for the approximation of the f i r s t term on the right hand side of equation ( I I I - l ) , namely, r . A ( * 2£) l r 3z 2fw dz t=T J z=z W Li+1*' 1+1,T - W 111 l i + l W1,T " W1-1,T h. V ( * h l+1 + * h 1 ) ^ v ( z 1 ) (III-6) And for the second term on the right hand side of equation (III-l) a backward f i n i t e difference quotient w i l l be used, 17.1 l r 3t m v ( Z l ) T . A T - T v * (III-7) It can be shown that by combining equations ( I I I - l ) , ( I I I - 5 ) , (III-6) and (III-7) a f i n i t e difference equation of the following form can be obtained, W i . T + A T " M V l + l . I + "l-l.T * <Vl>-Wi.T] + W 1 . T + A M 1 . T ( I 1 I - 8 » where, R, = / — (III-9) 1 h i + l / h i 2 AT N = — — r . ~ (111-10) 1 *w' ( hi+l + hi)»n v(z 1) h 1 a n d W f , T = - ^ U l f ( 1 1 1 - U ) It can be seen that equation (III-8) resembles the usual f i n i t e difference solution of a one dimensional consolidation problem except for the last term which represents the •generation 1 of pore water pressure. In order that the values of mv and £ v p at the centre of the layer can be used instead of the respective values at the nodal points the following average values are used, / . mv,i+l' h H-l * m v , i ' h l ffiv<zl> = h 1 + 1 + h l (111-12) 172 cvp 1 T m (z, )_ v I T A£vp,l+1,T A G v p , i , ' m v.i+l m v . i (111-13) Expressions (III-10) and ( I I I - l l ) can be modified into the following, 2. AT * i r w- ( m v, i+l- h i+l + m v , i ' V ' h i (111-14) AW g l.T = 4. vp.i+l.T ASrp,i,T + m v,i+l m v . l (111-15) The f i n a l form of the equation is s t i l l (III-8) but expressions (III-9), (111-14) and (111-15) w i l l be used in the analysis. L I S T I N G O F L I Q U E F A C T I O N P R O G R A M - A U G , 7 5 L I Q P G 1 173 APPENDIX IV C**************************************** C* I N P U T C A R D S F O R L I Q U E F A C T I O N P R O G R A M * C**************************************** C C 1 - T I T L E ( 2 0 A U ) (1 C A R D ) C C 2 . N M A T , N P T Y P E , N V O L , N L A Y E R , N L I N E L , N D A M P , N B T , N B B , N C R V . . ( 2 0 I U ) ( 1 C C W H A T = N U M B E R O F M A T E R I A L T Y P E S . C N P T Y P E = P R O B K E M T Y P E N U M B E R , C 1 F O R D Y N A M I C R E S P O N S E O N L Y , C 2 F O R D Y N . R S S P . P L U S P O R E P R E S S . R I S E O R V O L . C H A N G C 3' F O R D Y N . R E S P . P L U S P O R E P R E S S . P L U S D I S S I P A T I O N . C N V O L = 0 F O R P O R E P R E S S U R E R I S E . C = 1 F O R V O L U M E C H A N G E ( S E T T L E M E N T ) . C N L A Y E R = N U M B E R O F L A Y E R S . C N L I N E L = 0 F O R N O N - L I N E A R M A T E R I A L ( C H A N G E M O D U L I ) C = 1 F O R L I N E A R M A T E R I A L ( C O N S T A N T M O D U L I ) C N D A M P = 1 F O R D A M P I N G C O E F F I C I E N T S I N P U T . C = 2 F O R D A M P I N G P R O P P O R T I O N A L T O M A S S A N D S T I F F N E S S . C N B T = 0 F O R Z E R O P O R E P R E S S U R E A T T O P B O U N D A R Y . C = 1 F O R Z E R O G R A D I E N T A T T O P B O U N D A R Y . C N B B = 0 F O R Z E R O P O R E P R E S S U R E A T B O T T O M B O U N D A R Y , C = ' 1 F O R Z E R O G R A D I E N T A T B O T T O M B O U N D A R Y . C N C R V = 0 F O R N O T C A L L I N G S U B R O U T I N E T O C H E C K S T R A I N R E V E R S A L C = 1 F O R C A L L I N G S U B R O U T I N E T O C H E C K F O R S T R A I N R E V E R S A L C C 3. M A T Y P ( I ) , N S U B D ( I ) , H ( I ) , W I D T H ( I ) ( 2 1 4 , 2 F 1 0 . U) , C ( N L A Y E R C A R D S ) C M A T Y P ( I ) = M A T E R I A L T Y P E N U M B E R O F I - T H L A Y E R . C N S U B D ( I ) = N U M B E R O F S U B D I V I S I O N O F I - T H L A Y E R F O R P P D I S S I P A T C H ( I ) = T H I C K N E S S O F I - T H L A Y E R . C W I D T H ( I ) = W I D T H O F I - T H L A Y E R . C U. A 1 ( I ) , A 2 ( I ) , A 3 ( I ) , B 1 ( I ) , B 2 ( I ) , B 3 ( I ) , C 1 ( I ) , C 2 ( I ) , C 3 ( I ) , C 4 ( I ) , C D 1 ( I ) , D 2 ( I ) , D E N - H ( I ) , D E N - V ( I ) , P E R M ( I ) , A L F A ( I ) , B E T A ( I ) . C C ( 6 F 1 2 . 4 ) C ( 3 * N M A T C A R D S C 5 . N E Q , I N T Y P , N C , N C P R , N C P R M , N P C O N , N C O N T R N P L D , D T — ( 5 I 4 , F 1 0 . 3 ) ( 1 C A C N E Q = 0 F O R C O N S T A N T A M P L I T U D E H A R M O N I C I N P U T . C " = 1 F O R E A R T H Q U A K E R E C O R D I N P U T . C = 2 F O R G R A D U A L L Y I N C R E A S I N G A M P L I T U D E T O T Y P E 1 . C I N T Y P = 0 F O R W I L S O N A N D C L O U G H ' S M E T H O D O F I N T E G R A T I O N . C = 1 F O R N E W M A R K ' S M E T H O D O F I N T E G R A T I O N . C N C = N U M B E R O F I N T E G R A T I O N S C N C P P = N U M B E R O F I N T E G R A T I O N S F O R E A C H P R I N T O U T C N C P R M = N U M B E R O F I N T E G R A T I O N S F O R P R I N T I N G M A X I M U M V A L . C N P C O N = 1 F O R P R I N T I N G O U T F I N A L V A L U E S F O R C O N T I N U I N G N E X T T C = 0 F O R N O P R I N T I N G O U T C N C O N T = 0 F O R Z E R O I N I T I A L C O N D I T I O N S . C = 1 F O R N O N - Z E R O I N I T I A L C O N D I T I O N S . C N P L D = 0 F O R N O O U T P U T O N U N I T 7 F O R P L O T . C = 1 F O R O U T P U T O N U N I T 7 F O R P L O T . C LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 2 174 C 6. FOR NEQ=1, NCARD,NREC,NFTS,RDFR C (2014) (1 CARD) C NCARD = NUMBER OF EQ REC CARDS. C NREC = NUMBER OF EQ REC PER CARD. C NFTS = 0 GOR NUMBER OF GEE'S INPUT. C . RDFR = FACTOR TO BE MULTIPLIED TO ACC RECORDS. C ' = 1 FOR FPS UNIT INPUT. C FMT C (20A4) C FORMAT FOR EQ REC INPUT. C NEQ=0 C NEQ=0, NC, NFTS,AMP,OMEGA C (2014) C NC NUMBER OF SINE CYCLES. C NCFTS C NFTS = 0 FOR NUMBER OF GEE'S EQ INPUT C = 1 FOR FPS UNIT INPUT C AMP = AMPLITUDE OF SINUSOIDAL INPUT C OMEGA - CRICULAR FREQ. OF SIN. INPUT. C 7. TITLE C (20A4) (1 CARD) C DESCRIPTION FOR BASE MOTION INPUT. C C 8. NEQ=0 S NCONT=0 THAT IS ALL NO MORE CARDS REQD. C NEQ=0 & NCONT=1 INITIAL VALUES INPUT ACCORDING TO TM1VAL. C NEQ=1 & NCONT=0 INPUT EARTHQUAKE RECORDS. C NEQ= 1 5 NCONT=1 INPUT INITIAL. VALUES AND EQ. REC. C** DYNAMIC ANALYSIS FOR LIQUEFACTION OF A HORIZONTAL SOIL DEPOSIT C SIX TYPES OF PROBLEMS CAN BE SOLVED, NAMELY: C 1. DYNAMIC RESPONSE OF SOIL DEPOSIT ONLY C 2. DYNAMIC RESPONSE WITH PORE PRESSURE GENERATION C 3 . DYNAMIC RESPONSE WITH PORE PRESSURE GENERATION AND DISSIPATION C AND FOR EACH OF THE ABOVE EITHER C I) LINEAR CASE (NLINEL= 1) , I.E. STIFFNESS ETC IS CONSTANT C II) NON-LINEAR CASE (NLINEL=0) C SO NPTYPE=1 CORRESPONDS TO TYPE ONE PROBLEM, ETC. C LCPPG IS USED TO CONTROL CALLING OF PORE PRESS GEN SUB-PROG. C LPPDIS IS USED TO CONTROL CALLING OF DISSIPATION SUB-PROG. C LPLOT IS USED TO CONTROL OUTPUT FOR PLOTTING PROGRAM. C INTYP IS USED TO INDICATE TYPE OF INTEGRATION SCHEME USED: C INTYP=0 WILSON AND CLOITGH'S METHOD IS USED C INTYP=1 NEWMARK'S METHOD OF INTEGRATION IS USED. C LDOF1 IS USED TO INDICATE ONE-DEGREE OF FREEDOM C LOGICAL* 1 LVDIR ( 2 0 ) ,LREV ( 2 0 ) ,LRG ( 2 0) ,DODEVP ( 2 0 ) , LLIQ ( 2 0 ) , LREAD,LCPPG,LPPDIS,LPLOT,LDOF1,LLIN,LSTOP,LMCLG, LJRV,LCKF3,LIT,LPMAX COMMON /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX . /LOGCOM/LVDIR,LR EV,LRG,DODEVP,LLIQ . /PROPTY/TITLE (20) ,PROP(20,17) , H ( 2 0 ) ,DEPTH ( 2 0 ) ,WIDTH ( 2 0 ) ,WT ( 2 0 ) , MATYP(20) ,SIGZ ( 2 0 ) ,CC ( 2 0 ) ,GG ( 2 0 ) , NSUBD ( 2 0 ) , GAMR ( 2 0 ) , LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 3 175 c 1 2 9001 1003 1004 7 C** C C c** 1001 8 LLIN=.TRUE. LPPDIS=.TRUE. LDOF1=. TRUE. LJRV=.TRUE, TAUR(20) ,CM(40) ,SK(40) ,GMN(20) ,TMN(20) /DYNEQT/DYM (40) /PEP (20) ,FORER(20) /CALV AL/DIRN (20) ,XD2 (20) , ABXD2 (20) ,XD1 (20) ,XD0 (20) ,BD(20) , BA (20) ,BB (20) , GD2 (20) , GD 1 (2 0) , G A MM A (20) , G A M M AO (20) , GAMMAX (2 0) ,TAXY(20) ,EVP (20) ,SK1R (20) , DPP (20) ,PP(20) /CPVALU/BAP(20) ,BBP(20) ,XD2P (20),XD1P (20) ,XD0P (20) ,GD2P (20) , GD1P (20) ,TAXYP (20) /REV PA R/NRPT (20) ,GGAMR(20,8) ,TTAUR(20,8) ,DIRNR(20,8) DIMENSION FMT (20) , AE (10) , DTJMIN (17) DATA LREAD/.TRUE./,LCPPG/.FALSE./, LPPDIS/. FALSE./,LPLOT/. FALSE./ DATA TNOT/0.0/,UG1/0.0/,NCPRC/0/,NRECC/0/,NPRMC/0/ READ(5,1) TITLE FORMAT(20 AU) WRITE (6,2) TITLE FORMAT(1H1,10X,20A4) READ (5,3) NMAT,NPTYPE,NVOL,NLAYER,NLINEL,NDAMP,NBT,NBB,ITERRV, ICKFB FORMAT (2014) I F (NLINEL.EQ. 1) IF (NPTYPE.GE.3) IF (NLAYER. EQ. 1) LJRV = .FALSE. LCKFB = .FALSE. I F (ITERRV. EQ. 0) IF (ICKFB. GE.1) LCKFB=. TRUE. IF (LLIN) WRITE (6, 9001) FORMAT ('0** LINEAR MATERIAL IS USED **•) N = NLAYER N2=N*2 N21=N2-1 IF(NLAYER.GE.1) GO TO 1003 WRITE (6,6) FORMAT(1H0,•** ERROR 2, ERROR IN INPUT **•) STOP DO 1004 1=1,N READ(5,7) MATYP (I) ,NSUBD (I) ,H (I) , WIDTH (I) FORMAT(2I4,2F10.4) IF (NMAT.GE.1) GO TO 1001 WRITE (6,4) FORMAT (1 HO,'** ERROR 1,ERROR IN INPUT **«) STOP ORDER OF INPUT CONSTANTS ARE: A 1,A2,A3,B1,B2,B3,C1,C2,C3,C4,D1,D2, DEN-H,DEN-V,PERM.,ALFA,BETA WRITE INPUT DATA WRITE (6,8) NMAT,NPTYPEfNVOLrNLAYER,NLINEL,NDAMP,NBT, NBB, ITERRV, ICKFB FORMAT ('0 NMAT =',13,' ; NPTYPE =',13,' ; NVOL = ' , 1 3 , ' ; NLAYER =«,I3,' ; NLINEL =',I3,' ;'/' NDAMP = ' , 1 3 , « ; NBT =',I3,« ; NBB =',I3, 1 ; ITERRV = ' , 1 3 , ' ; ICKFB =',13) WRITE(6,9) LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 4 176 9 FORMAT(1H0,'***********************'/* * MATERIAL PROPERTIES *'/ i * * * * * * * * * * * * * * * * * * * * * * * 1 C IF NDAMP=1, CC INPUT AS FORCE/UNIT VELOCITY C - 2, ALFA=PROP(1, 16) AND BETA=PROP (1,17) DO 1002 I=1,NMAT READ(5,5) (DUMIN (IJ) ,IJ=1, 17) 5 FORMAT (6F12.U) WRITE(6,10) I , (DUMIN (I J) , IJ= 1,17) 10 FORMAT (1X, 'MATERIAL TYPE', 13,' : A " S AND B«'S : ' , 1 P6 E1 4. 3/20 X , •C^S AND D"S :',6E14. 3/20 X, • DEN. - H = ,,0PF5.1,« PCF;',3X, •DEN.-V =',F5.1,« P C F ; PERM. =',1PE10.3,' ; ', * ALFA =',E11.3,3X,'BETA = ',E11.3) DO 1012 IJ=1,N IF (MATYP (IJ) . NE. I) GO TO 1013 DO 1014 IK=1,17 1014 PROP (IJ,IK) =DUMIN (IK) 1013 CONTINUE 1012 CONTINUE 1002 CONTINUE WRITE(6,11) 11 FORMAT(1H0,'*********************/1X,'* LAYER PROPERTIES *'/1X, m t * * * * * * * * * * * * * * * * * * * * i / j WRITE (6, 12) 12 FORMAT (2X f'LAY.NO. MATYP NSUBD THICK WIDTH') WRITE (6,13) (I, MATYP (I) , NSUBD (I) ,H(I) , WIDTH (I) ,I=1,N) 13 FORMAT(1X,I4,3X,2I8,2F10.2) C C ASSIGN LOWER MOST BOUNDARY AS PROBLEM' BOUNDARY C FORM DIAGONAL MASS MATRIX C FIND DEPTH AND SIGZ AT CENTRE OF EACH LAYER C DO 1007 1=1,N I J = N-I+1 IK = IJ+1 DT2 = 0. 5*H (IJ) UGI = DT2*PROP (I J , 13) CA = DT2*PROP ( I J , 14) IF (I.GT. 1) GO TO 1009 WT(N) = UGI DEPTH (N) = DT2 SIGZ(N) = CA GO TO 1016 1009 WT(IJ) = UGI+UG1 DEPTH(IJ) = DEPTH (IK) +DT1 + DT2 SIGZ(IJ) = SIGZ (IJ+1) +CA+BK 1016 DT1 = DT2 UG1 = UGI BK = CA 1007 CONTINUE C C CHANGE MASS UNIT TO LB-SEC.SQ./FT. UNIT CA = 1./32. 172 DO 1033 1=1,N LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 177 1033 WT(I) = WT(I) *CA*WTDTH(I) C C FIND DAMPING VALUE AND STIFFNESS VALUE FOP. EACH LAYER C** FORM STIFFNESS AND DAMPING MATRIX. NOTE HALF BAND WIDTH = 2 C CALL SYSCKM (NCONT) C C WRITE INITIAL PROPERTIES OF THE LUMPED MASS SYSTEM WRITE(6,15) 15 FORMAT(1 H1,• ** INITIAL PROPERTIES OF THE LUMPED MASS SYSTEM** /1H0,' NO. ALFA BETA»,10X,»G•,11X,•C»,11X,•M' 7X,'DEPTH',8X,•SIGZ'/) WRITE (6, 16) (I,PROP (1,16) ,PROP (1,17) ,-GG (I) ,CC(I) , WT (I) ,DEPTH (I) SIGZ(I) ,1=1 ,N) 16 FORMAT(1X,14,1P7E12.3) WRITE (6,24) 24 FORMAT('0***INITIAL STIFFNESS MATRIX***') J = 2*N-1 WRITE (6, 25) (SK (I) ,1=1, J , 2) 25 FORM AT(' DIAGONAL TERMS :•/(10X,1P10E11.3)) WRITE(6,26) (SK (I) ,1=2, J , 2) 26 FORMAT(' OFF-DIAGONAL TERMS :•/(10X,1P10E11.3) ) WRITE(6,27) 27 FORMAT ('0***INITIAL DAMPING MATRIX***') WRITE(6,25) (CM (I) ,1=1, J , 2) WRITE(6,26) (CM (I) ,1=2, J,2) C C** . START DYNAMIC ANALYSIS C READ (5,18) NEQ,INTYP,NC,NCPR,NCPRM,NPCON,NCONT,NPLD,DT 18 FORMAT(8I4,3F10.0) WRITE (6,28) NEQ,INTYP,NC,NCPR,NCPRM,NPCON,NCONT,NPLD 28 FORMAT ('6*** CONTROL NUMBERS FOR DYNAMIC RESPONSE ANALYSIS ***• 10X,'NEQ =',I4,» ; INTYP =',I4,« ; NC =',14 ' ; NCPR =',I4,» ;•/10 X,'NCPRM =',14,' ; NPCON = 14,' ; NCONT =',14,' ; NPLD =',I4,' ;•) LPLOT = .FALSE. RDFR = 1. IF (NPLD. EQ. 1) LPLOT=. TRUE. IF (NCPRM.EQ.0) NCPRM=NCPR IF (NEQ. NE. 1) GO TO 1056 READ(5,30) NCARD,NREC,NFTS,RDFR,FMT 30 FORMAT(3I4,F10. 0/20A4) WRITE (6,29) NCARD,NREC,NFTS,RDFR,FMT 29 FORMAT('0*** EARTHQUAKE RECORD USED AS INPUT ***'/10X,'NCARD =• 14,' ; NREC =',I4,» ; NFTS =',I2,' ; MULTIPLYING 'FACTOR =',F5.2/10X,'FORMAT OF INPUT : ',20A4) NCINP = NCARD*NREC GO TO 1057 1056 READ(5,22) NCINP,NFTS,AMP,OMEGA 22 FORMAT (214,2F10.4) LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 1?8 WRITE(6,31) NCINP,AMP,OMEGA,NFTS 31 FORMAT(* 0*** SINUSOIDAL BASE INPUT IS USED ***»/10X, , . * NUMBER OF BASE ACC. INPUT =',15,' ; AMPLITUDE = », 1PE11.3,' ; OMEGA =\,E11.3,' ; NFTS =',0PI2,« .») NREC = 1 0 1057 READ (5,3) NS.UBDT WRITE (6,3333) NSUBDT IF (NSUBDT.EQ.0) NSUBDT=1 3333 FORMAT(»0***** SUB-DIVISION OF TIME INTERVAL FOR », •CONSOLIDATION =',I3,» ****••) ACSC = RDFR I F (NFTS.EQ.0) ACSC=32.17*RDFR I F (DT.LE.0.1) GO TO 1050 WRITE(6,19) DT 19 FORMAT(»0*** TIME INTERVAL TOO LARGE, DT=«,F6.3,» ! » * * * ' ) STOP 1050 READ (5,1) TITLE WRITE(6,23) TITLE 23 FORMAT(«0*****»,20A4,•*****»/) WRITE (6,20) DT 20 FORMAT (* 0** TIME INTERVAL USED, DT= ,,F6.3,» SEC. ***») IF(NCONT.EQ.0) GO TO 1052 CALL TM1VAL(UG1,DT,TNOT) C C INITIAL CYCLE C 1052 WRITE(6,21) 21 FORMAT(1H1,* OUTPUT VALUES ARE IN THE FOLLOWING ORDER : **'/ » ** LAY.NO.,ACC., VEL., DISPL., GAMMA., TAU, EVP, PP. . • **«//) C C START REPETITIVE INTEGRATION, (PP GEN. AND PP DIS. IF REQUIRED) C UG1 = 0.0 100 IF (. NOT. LREAD) GO TO 1034 CALL GETACC(AE,FMT,DT,AMP,OMEGA,NEQ,NREC,NCINP,TNOT) DO 1051 1=1,NREC 1051 AE(I) = ACSC*AE(I) LREAD = .FALSE. 1034 NRECC = NRECC+1 NCPRC = NCPRC+1 NPRMC = NPRMC+1 INCR = INCR+1 ITERNO = 0 CA = INCR TIME = CA*DT+TNOT DDT = DT IF (INCR.GT.NC) GO TO 9999 DGI = AE (NRECC) -UG 1 UGII = UGI UG1 = AE (NRECC) IF (NRECC. NE. NREC) GO TO 1023 NRECC = 0 LISTING OF LIQUEFACTION PROGRAM -• AUG, 75 LIQPG 7 179 LREAD = .TRUE. 1023 CONTINUE C C************** 1054 CALL SYSCKM(NCONT) r j * * * * *********** C 1029 DO 1063 1=1,N BAP (I) = BA(I) BBP(I) = BB(I) XD2P(I) = XD2(I) XD1P (I) = XD1 (I) XDOP (I) = XDO (I) GD2P(I) = GD2(I) GD1P (I) = GD1 (I) GAMMAO (I) = GAMMA (I) 1063 TAXYP(I) = TAXY(I) C CALL NUMERICAL INTEGRATION SCHEME C c * * * * * * * * * * * * * * * * LIT = .FALSE. 1064 CALL INTSCH(DDT,UG1,UGII) C**************** C IF (LJRV) GO TO 1062 C ( 3 * * * * * * * * * * * * * * * * 1065 CALL CHEKRV (DT,DDT,UG1,UGT,UGII,ITERNO,NRETN) C**************** C GO TO (1062,1054,1064,1065),NRETN WRITE(6 ,93) NRETN 93 FORMAT (' 0 * * SOMETHING WRONG, NRETN = ,,I4,» ***•) STOP C 1062 CONTINUE IF (NPTYPE. LT. 2) GO TO 1035 LSTOP = .FALSE. C c * * * * * * * * * * * * * * * CALL PPGEN(TIME,LCPPG,LSTOP) IF(LSTOP) GOTO 9 9 9 9 C*************** c 1061 IF (.NOT.LPPDIS) GO TO 1058 C c * * * * * * * * * * * * * * * * CALL PPDISN(DT,NSUBDT,LCPPG) C**************** C GO TO 1035 1058 IF (.NOT.LCPPG) GO TO 1035 DO 1059 1=1,N LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 8 180 IF (.NOT.DODEVP(I)) GO TO 1059 PP (I) = PP (I) +DPP (I) 1059 CONTINUE C , 1035 DO 1666 1=1,N IF (PP(I) .GT,SIGZ(I)) PP (I) =SIGZ (I) 1666 CONTINUE C • . LPMAX = .FALSE. CALL MAXVAL (TIME,UG1) C IF(LPLOT) WRITE (7) TIME, UG1 , (I, ABXD2 (I) , XD1 (I) , XDO (I) , GAMMA (I) , TAXY (I) , EVP (I) ,PP(I) ,I=1,N) 1055 IF (NCPRC.LT.NCPR) GO TO 1667 NCPRC = 0 C DISPLACEMENT IS REFERRED TO THE TOP OF THE SOIL LAYER WRITE(6,90) TIME,UG1 , GG (1) 90 FORMAT (1H0, •TIME =,,F6.3,' SEC; ACC. =,,F8.t»,' FPSS', 1P2E13.4) IF (NPTYPE.LE. 1) WRITE (6,91) (I, ABXD2 (I) ,XD1 (I) ,XD0 (I) , GAMMA (I) , TAXY (I) ,1=1, N) 91 FORMAT (1X,IU, 1P5E1 2. 3) IF (NPTYPE. GT. 1) WRITE(6,92) (I,ABXD2(I) ,XD1 (I) ,XD0 (I) ,GAMMA (I) , TAXY (I) ,EVP (I) ,PP (I) ,1=1 ,N) 92 F0RMAT(1X,IU,1P7E12.3) 1667 IF (NPRMC.LT.NCPRM) GO TO 100. NPRMC = 0 LPMAX = .TRUE. CALL MAXVAL (TIME,UG1) . GO TO 100 9999 IF (NPCON.NE.1) GO TO 9998 WRITE(8) TIME,UG 1 , DT DO 1060 1=1,N 1060 WRITE (8) I,XD2 (I) ,XD1 (I) ,XD0 (I) ,TAXY (I) , EVP (I) , PP (I) ,GAMMA (I) , GAMMAO (I) , GAMR (I) ,TAUR (I) , DIRN (I) 9998 IF (LPLOT) END FILE 7 C LPMAX = .TRUE. CALL MAXVAL (TIME,UG1) C WRITE (6,150) INCR 150 FORMAT(* 0******************************************* */ • * NORMAL TERMINATION OF PROGRAM *'/ » * NUMBER OF INCREMENTAL CALCULATION =',14,' *•/ • * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * i ) STOP C *********************** C * END OF MAIN PROGRAM * C *********************** END BLOCK DATA LOGICAL*1 LVDIR (20) ,LREV (20) ,LRG (20) , DODEVP (20) ,LLIQ (20) , LDOF1,LLIN,LSTOP,LMCLG LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 9 181 COMMON . /NCONTL/N,NDAM P,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG . /LOGCOM/LVDIR,LREV,LRG fDODEVP,LLIQ „ /PROPTY/TITLE (20) ,PROP(20, 17) ,H (20) ,DEPTH(20) ,WIDTH(20) ,WT (20) , . MATYP(20),SIGZ(20) ,CC(20) ,GG (20) ,NSUBD(20) ,GAMR(20) , TAUR (20) ,CM (40) , SK (40) ,GMN (20) ,TMN (20) . /DYNEQT/DYM (40) , PEP (20) ,FOREP (20) . /CALVAL/DIRN(20) ,XD2(20) ,ABXD2(20) ,XD1 (20) ,XD0(20) ,BD(20) , BA (20) ,BB(20) ,GD2 (20) ,GD1 (20) ,GAMMA (20) , GAM MAO (20) , . GAMMAX (20) ,TAXY (20) , EVP (20) ,SK1R (20) ,DPP (20) ,PP (20) . /CPVALU/BAP (20) ,BBP(20) ,XD2P (20) ,XD1P(20) ,XD0P(20) ,GD2P(20) , GD1P(20) ,TAXYP(20) . , /REVPAR/NRPT (20) ,GGAMR (20, 8) ,TTAUR (2 0,8) , DIRNR (20,8) DIMENSION BLOC(360) ,BLOC2 (160) EQUIVALENCE (BLOC (1) , DIRN (1) ) , (BLOC2 (1) ,BAP (1) ) DATA INCP/0/, LDOF1/.FALSE./, LLIN/.FALSE./, LSTOP/.FALSE./, LMCLG/.FALSE./, LVDIR/20*.FALSE./, LEEV/2 0*.FALSE./, LRG/20*.FALSE./, DODEVP/20*.FALSE./, LLIQ/20*.FALSE./, DEPTH/20*0.0/, WT/20*0.0/, SIGZ/20*0.0/, CC/20*0.0/, GG/20*0.0/, BLOC/360*0.0/, BLOC2/16O*0.0/, FORER/20*0.0/, NRPT/20*0/, GGAMR/20*0.0/, TTAUR/20*0.0/ END SUBROUTINE TM1VAL(TNOT,UG1,DT) COMMON " . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /PROPTY/TITLE(20) ,PROP(20,17) ,H(20) ,DEPTH (20) ,WIDTH (20) ,WT(20) , HATYP(20) ,SIGZ (20) ,CC (2 0) ,GG (20) ,NSUBD (20) ,GAMP. (20) , .. TAUR (20) ,CM (40) , SK (40) ,GMN (20) ,TMN (20) . /CALVAL/DIRN (20) ,XD2 (20) ,ABXD2 (20) ,XD1 (20) ,XD0 (20) ,BD(20) , BA(20) ,BB(20) ,GD2(20) ,GD1 (20) ,GAMMA(20) ,GAMMAO(20) , GAMMAX (20) , TAXY (20) ,EVP (20) , SK1R (20) , DPP (2 0) ,PP (20) . /REVPAR/NRPT (20) ,GGAMR (20,8) ,TTAUR (20,8) , DIRNR (20,8) WRITE (6,21) 21 FORMAT ('0*** SUBROUTINE TM1VAL IS CALLED ***•/) READ (8) TN0T,UG1,DT1 WRTTE(6,22) TN0T,UG1,DT1 22 FORMAT(* TIME =',F7.3,' SEC.; A C C = « , F 7 . 3 , ' FPSS; TIME', » INTERVAL = ,,F5.3,' SEC.'/) TEM = ABS (DT-DT1) IF (TEM. GT. 0. 0005) GO TO 888 WRITE (6,23) 23 FORMAT('ONO. ACC',9X,'VEL•,8X,'DISP',7X,'TAUXY•,9X,'EVP•, 10X,*PP',7X,'GAMMA',6X,'GAMMAO',8X,'GAMR »,8X,'TAUR', 7X,'DIRN'/) C DO 101 1=1,N READ (8) a,XD2 (I) , XD1 (I) ,XD0 (I) , TAXY (I) ,EVP(I) ,PP(I) ,GAMMA (I) , GAMMAO (I) , GAMR (I) , TAUR (I) , DIRN (I) IF ( J . NE.I) GO TO 887 BA (I) = XD2 (I) *DT BB(I) = (XD1 (I)+0. 5*BA (I) ) *DT LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 10 182 101 WRITE(6,24) J,XD2(I) ,XD1 (I) ,XD0(I) ,TAXY(I) ,EVP(I) ,PP(I) , GAMMA (I) , GAMMAO(T) , GAMR (I) ,TAUR(I) ,DIRN(I) 24 FORMAT (1X,12,1X,1P10E12.4,0PF6.1) RETURN' 888 WRITE (6, 25) 25 FORMAT (»0**** SOMETHING WRONG, DT"S ARE NOT EQUAL !!! ****«) GO TO 999 887 WRITE(6,26) 26 FORMAT(* 0**** SOMETHING WRONG, I NOT EQUAL TO J I!! ****») 999 STOP END SUBROUTINE GETACC{AE,FMT,DT,AMP,OMEGA,NEQ,NREC,NCINP,TNOT) LOGICAL*1 LDOF1,LLIN,LSTOP,LMCLG,LO,L1,L2,L3 COMMON . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG DIMENSION AE (10) ,FMT(20) DATA LO/.FALSE./,L1/.FALSE./,L2/.FALSE./,L3/.FALSE./ , IF (L3) GO TO 999 IF (INCR. GT. NCINP) GO TO 99 > IF (LO) GO TO 11 IF (L1) GO TO 12 IF (L2) GO TO 13 IF (NEQ. NE. 0) GO TO 14 LO = .TRUE. GO TO 11 14 • IF (NEQ.NE.1) GO TO 15 L1 •= .TRUE. GO TO 12 15 . FNEQ = FLOAT (NEQ) RNPI = FNEQ*3.1415926 RNPI = 1./RNPI L2 = .TRUE. GO TO 18 11 FI = INCR T = FI*DT DO 17 1=1,10 T = T+DT TEM = T+TNOT TEM = TEM*OKEGA AE(I) = AMP*SIN (TEM) 17 CONTINUE RETURN 12 READ (5, FMT, END=99) (AE (IJ) ,IJ= 1 , NREC) RETURN 18 IF (TNOT.LT.DT) GO TO 13 WRITE(6,102) TNOT 102 FORMAT(•0**** GRADUALLY INCREASING AMPLITUDE REQUIRED, BUT TNOT *, • IS NON-ZERO, TNOT = ,,F6.2,« SECS. ****') GO TO 99 13 FI = INCR T = FI*DT DO 19 1=1,10 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 11 183 T = T+DT TEM = T*OMEGA ACCT =. AMP*SIN (TEM) TEM = TEM*RNPI C TEM IS DIVIDED BY N PI IF (TEM. GE. 1.0) GO TO 21 TEM = TEM*1.570796 TEM = SIN (TEM) ACCT = ACCT*TEM GO TO 20 21 LO = .TRUE. L2 = .FALSE. 20 AE(I) = ACCT 19 CONTINUE RETURN 99 DO 22 1=1,10 22 AE(I) = 0.0 L3 = .TRUE. 999 RETURN END SUBROUTINE SYSCKM(NCONT) LOGICAL*1 LVDIR (20) ,LREV (20) ,LRG (20) ,DODEVP (20) ,LLIQ (20) , LDOF1,LLIN,LSTOP,LMCLG,LJUMP1,LJUMP2,L3,L4 COMMON . /NCONTL/N,NDA MP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG . /LOGCOM/LVDIR ,LR EV,LRG,DODEVP, LLTQ. . /PROPTY/TITLE(20) ,PROP(20,17) ,H(20) ,DEPTH(20),WIDTH(20) ,WT(20) , MATYP{20) ,SIGZ(20) ,CC(20) ,GG(20) ,NSUBD(20) ,GAMR(20) , TAUR (20) ,CM (40) ,SK (40) ,GMN (20) ,TMN (20) . /CALVAL/DIRN(20) ,XD2(20) , ABXD2(20) , XD1 (20) ,XDO (20) , BD (20) , BA(20) , BB ( 20) ,GD2(20) ,GD1 (20) ,GAMMA (20) ,GAM M AO (20) , GAMMAX (2 0)",TAXY (20) ,EVP (20) , SK1R (20) , DPP (20) ,PP (20) . /REVPAR/NRPT (20) ,GGAMR (20,8) , TT AUR (20,8) , DIR NR (20,8) DIMENSION FK (20) DATA LJUMP1/.FALSE./,LJUMP2/.FALSE./,KPRT/0/ DATA RLIM1/0.0500/, RLIM2/0.05/ C C THIS PORTION FINDS SHEAR MODULUS FOR A LAYER c IF (LJUMP1) GO TO 14 C C HERE IS FOR LINEAR CASE OR FIRST INCREMENT 27 DO 15 1=1,N GMN (I) = PROP (1,1) TMN(I) = PROP(I,4) 15 GG (I) = GMN (I) IF (.NOT. LLIN) LJUMP1 = .TRUE. GO TO 50 C C HERE IS NON-LINEAR CASE, AND INCR > 1 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 12 184 14 DO 16 1=1,N IF (.NOT. LVDIR (I) ) GO TO 16 IF (.NOT.DODEVP(I)) GO TO 17 EVPI = EVP (I) CEP 1 = 1. +EVPI/ (PROP (1,2) + PROP (1,3) *EVPI) CEP2 = 1.+EVPI/(PROP (1,5)+PROP (1,6) *EVPI) RATS = 1.-PP (I)/SIGZ (I) SATT = RATS IF (RATT.LT. RLIM1) RATT=RLIM 1 TMN(I) = PROP (1,4) *CEP2*RATT RATT = SQRT (RATS) IF (RATT. LT.RLIM2) RATT=RLIM2 GMN (I) = PROP (I, 1) *CEP1*RATT 17 CONTINUE J = NRPT(I) GMNT = GMN (I) TMNT = TMN (I) DIRNI = DIRN (I) GAMI = GAMMA (I) TAUI = TAXY (I) GAMAB = ABS (GAMI) TAUAB = ABS (TAUI) IF (. NOT. LRG (I) ) GO TO 20 LRG(I) = -FALSE. C C HERE REVERSAL OCCURS, INCREASE NRPT BY 1 19 IF (J.LT.R) GO TO 29 - . IF (J.EQ.8) GO TO 29 WRITE (6,2005) INCR,I,J 2005 FOR MAT (' **** AT INCR',14," FOR LAY.»,I3,« J > 8 AND =',19, » ****!) STOP 29 J = 7 NRPT (I) = J GO TO 40 28 TEM1 = GAMR(I)*GMNT/TMNT TEM2 = TAUR (I) /TMNT IF (J.EQ.O) GO TO 30 TDIF = TEM2-TTAUR (I, J) TPRD = DIRNR (I, J) *TDIF IF (TPRD.GT.0.0) GO TO 31 J = J+1 NRPT (I) = J 31 DIRNR (I, J) = SIGN (1. ,TDIF) GGAMR (I, J) = TEM 1 TTAUR (I, J) = TEM2 GO TO 40 30 J = 1 NRPT (I) = 1 DIRNR (1,1) = SIGN (1. ,TEM2) GGAMR (1,1) = TEM 1 TTAUR(I,1) = TEM 2 GO TO 40 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 13 185 C 20 IF (J.EQ.O) GO TO 18 C CHECK FOR REDUCTION IN NRPT, IE. CHECK FOR ENVELOP CURVES :-C NOTE THAT THE LAST BUT ONE REVDRSAL POINT IS ALSO THE LIMIT C FOR SWITCHING TO ENVELOP CURVE IF (J.GT. 1) GO TO 21 TEMP = ABS (TTAUR (I, 1) ) *TMNT IF (TAUAB.LT.TEMP) GO TO 40 J = 0 NRPT (I) = 0 GO TO 18 C 21 TEMP = TTAUR (I,J-1) *TMNT C CHECK FOR CONSISTENCY BETWEEN DIRN AND TTAUR'S TDIF = TEMP-TTAUR (I, J) *TMNT TPRD = DIRNI*TDIF IF (TPRD) 24,25,25 24 KPRT = KPRT+1 J = J-1 NRPT (I) = J IF (KPRT. GT. 50) GO TO 20 IF (J.NE.O) GO TO 26 WRITE(6,2004) INCR,I,J GO TO 20 26 J1 = J+1 WRITE (6,2004) TNCR,I,J, (TTAUR(I,K) ,K=1,J1) ,TAUI 2004 FORMAT (1X,314,1P9E13.3) GO TO 20 25 IF (TDIF) 22,23, 23 23 . IF (TAUI.LT.TEMP) GO TO 40 J = J-2 NRPT (I) = J GO TO 20 22 IF (TAUI.GT.TEMP) GO TO 40 J = J-2 NRPT (I) = J GO TO 20 40 TEMP = GGAMR(I,J)*TMNT/GMNT GGT = 1. 0/(1.0 + 0.005*GMNT*ABS (GAMI-TEMP)/TMNT) GO TO 116 18 GGT = 1.0/(1.0+0.01*GMNT*GAMAB/TMNT) 116 GGT = 0.98*GGT*GGT*GMNT RLG = 1./PROP (1,1) IF (GGT*RLG.LT.0.005) GGT=0.005/RLG GG(I) = GGT C IF (TAUAB. LT. TMNT) GO TO 16 KPRT = KPRT+1 IF (KPRT.GT.50) GO TO 16 IF (KPRT. EQ. 50) GO TO 45 IF (TAUI) 42,43,43 42 IF (DIRNI.GT.0.0) GO TO 16 WRITE(6,2001) INCR,I LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 186 2001 FORMAT(' ** AT INCR ' ,14 , ' LAY.• , 1 3 , • TAXY < -TMN *) GO TO 16 43 IF (DIRNI.LT.0.0) GO TO 16 WRITE(6,2002) INCR,I . , 2002 FORMAT (' ** AT INCR ' ,14, ' LAY. ' ,13 ,• TAXY > TMN') GO TO 16 45 WRITE (6,2003) 2003 FORMAT(* *** NO MORE WRITING FOR |TAXY| > |TMN| ***') C 16 CONTINUE C 50 CONTINUE C THIS PORTION FORMS SYSTEM DAMPING AND STIFFNESS MATRICES C C** CC AND G VECTORS CONTAIN LAYER DAMPING AND STIFFNESS VALUES C C AND K ARE BANDED DAMPING AND STIFFNESS MATRICES C IF (LJUMP2) GO TO 58 LJUMP2 = .TRUE. AA1 = PROP(1,16) BB1 = PROP (1,17) C DO 56 1=1,N IF (NDAMP.EQ.1) GO TO 57 CC(I) = BB1*GG(I) GO TO 56 57 CC(I) = PROP (I, 16) 56 CONTINUE DO 59 1=1,N 12 = 2*1 I2M1 = 12-1 CCI = CC (I) *WIDTH (I)/H (I) CCI1 = CC (1 + 1) *WIDTH (1+1)/H (1+1) CM(I2M1) = CCI IF (NDAMP.EQ. 2) CM (I2M1) =CM (I2M1) +AA1*WT (I) IF (I.GE.N) GO TO 59 CM(I2M1) = CM (I2M1)+CCI1 CM (12) = - e c u 59 CONTINUE C 58 DO 52 1=1,N 52 FK(I) = GG (I) *WIDTH (I)/H (I) C DO 51 1=1,N 12 = 2*1 I2M1 = 12-1 FKI = FK(I) FKI1 = FK (1+1) SK (I2H1) = FKI . IF (I.GE.N) GO TO 51 SK(I2M1) = SK (12M1) + FK11 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 15 187 SK(I2) = -FKI1 51 CONTINUE C IF (INCR. LT. 113) GO TO 9999 C IF (INCR.GT. 120) GO TO 9999 ; C WRITE(6,2999) INCR, (GAMMA (K) , K= 1 , N) C WRITE (6,2999) INCR, (TAXY(K) ,K=1,N) C2999 FORMAT (1X,I4,1P10E12.3/(5X,10E12. 3) ) 9999 RETURN END SUBROUTINE SOLB2 COMMON /NCONTL/N,NDAMP,INCR,INTYP,N70L,NBT,NBB . /DYNEQT/DYM (40) ,PEP(20) ,FORER(20) DIMENSION T(20) C DYM BECOMES THE LOWER TRIANGULAR MATRIX C T IS THE UPPER TRIANGULAR MATRIX C NOTE L33 IS DYM (5) , ETC, T12 IS T(1), T23 IS T (2) ETC. IF (ABS(DYM(1)).GT.1.0E-75) GO TO 103 WRITE (6, 10) DYM (1) 10 FORMAT («0*** ZERO ENCOUNTERED IN ROW 1 !!! *** DYM (1)=•,E12.4) STOP 103 CONTINUE PEP(1) = PEP (1)/DYM (1) IF(N.EQ.1) RETURN DO 101 1=2,N IM1 = 1-1 MN = 2*IM1 NN = MN+1 MM = MN-1 AMN = DYM (MN) • T (IM1) = AMN/DYM (MM) DYM (NN) = DYM (NN)-T (IM1) *AMN IF (ABS (DYM (NN) ) .GT. 1. OE-75) GO TO 104 WRITE(6,11) I, NN , DYM (NN) 11 FORMAT (* 0*** ZERO ENCOUNTERED IN ROW',13,' !!! *** DYM(',I2, •)=',E12.4) STOP 104 CONTINUE PEP(I) = (PEP(I) -AMN*PEP(IM1) )/DYM(NN) 101 CONTINUE C BACK SUBSTITUTION DO 102 1=2,N MN = N-I + 1 NN = MN+1 102 PEP,(MN) = PEP(MN) -T(MN) *PEP(NN) RETURN END SUBROUTINE CHEKSV(DT,DDT,UG1,UGI,UGII,ITER NO,NRETN) LOGICAL*1 LDOF1,LLIN,LSTOP,LMCLG,LJUMP,LITE,LCPPG,LPRT, LVDIR (20) , LREV (20) , LRG (20) ,DODSVP (20) , LLIQ (20) , LRGL(20),LJRV,LCKF3,LIT,LPMAX COMMON /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 16 188 /LOGCOM/LVDIR,LREV,LRG,DODEVP,LLIQ /PROPTT/TITLS (20) , PROP (20, 17) , H (20) , DEPTH (20) , WIDTH (20) ,WT (20) , MATYP (2 0) ,SIGZ (20) ,CC (20) ,GG (20) , NSUBD (20) , G AMR (20) , TAUR (20) ,CM (40) ,SK{40) ,GMN (20) ,TMN (20) /DYNEQT/DYM (40) ,PEP (20) ,FORER (20) /CALVAL/DIRN (20) ,XD2(20) ,ABXD2(20) ,XD1 (20) , XDO (20) ,BD(20) , BA (20) , BB (20) ,GD2 (20) ,GD1 (20) ,GAMMA (20) ,GAMMAO (20) , GAMWAX(20) ,TAXY(20) ,EVP(20) ,SK1R(20) ,DPP(20) ,PP(?0) . /CPVALU/BAP (20) ,BBP (20) ,XD2P (20) ,XD1P (20) ,XD0P (20) ,GD2P (20) , GD1P (20) ,TAXYP (20) DIMENSION NREV (20) ,GAMRO (20) ,TAURO (20) C DATA JCALL/0/,LJUMP/.FALSE./,LPRT/.TRUE./,JPRINT/0/, LRGL/20*. FALSE. / C IF (INCR. EQ. 50) LPRT=. FALSE. IF (LJUMP) GO TO 100 IF (JCALL.GE.1) GO TO 101 DO 102 1=1,N 102 LVDIR (I) = . FALSE. 101 DO 103 1=1,N IF (GAMMA (I) .EQ.0.0) GO TO 103 IF (LVDIR (I) ) GO TO 103 DIRN(I) = SIGN (1. , GAMMA (I) ) LVDIR (I) = .TRUE. 103 CONTINUE JCALL = JCALL+1 IF (JCALL. EQ. 1) GO TO 131 DO 104 1=1,N • IF (. NOT. LVDIR (I) ) GO TO 100 104 CONTINUE LJUMP = .TRUE. 100 IF (ITERNO.NE.O) GO TO 111 NDT = 0 NDDT = • 10 NTDT = 0 111 ITERNO = ITERNO+1 IF (ITERNO.LT.11) GO TO 124 WRITE (6,9) ITERNO 9 FORMAT(' ***ITERNO =',14,' ***») STOP 124 CONTINUE C CHECK GAMMA FOR REVERSAL LITE = .FALSE. DO 105 1=1,N LREV (I) = .FALSE. IF (LRGL (I) ) GO TO 105 IF (.NOT.LVDIR(I) ) GO TO 105 TEMY = GAMMA (I)-GAMMAO (I) IF (TEMY. EQ. 0.0) GO TO 105 TEMY = DIRN (I) *TEMY IF (TEMY. GT. 0.0) GO TO 105 LREV (I) = .TRUE. LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 17 189 LRG (I) = .TRUE. GAMRO (I) = GAMR (I) TAURO (I) = TAUR (I) GAMR (I) = GAMMAO (I) TAUR (I) = TAXYP(I) GAMMAX (I) = GAMMAO (I) DIRN (I) = -DIRN (I) IF (ABS(GAMMAO(I)) .LT.1.E-8) GO TO 105 LITE = .TRUE. 105 CONTINUE IF ( (.NOT.LITE) .OR.NTDT. EQ. 10) GO TO 131 C SINCE RE-CALCULATION IS REQUIRED, CHANGE GAMR'S AND DIRN'S BACK C TO ITS VALUES BEFORE CHECKING FOR REVERSAL. DO 123 1=1,N IF (. NOT. LREV (I) ) GO TO 123 GAMR (I) = GAMRO (I) TAUR (I) = TAURO (I) DIRN(I) =-DIRN(I) 123 CONTINUE C NREV IS A NO. USED FOR INDICATING NEAREST TENTH OF (DT) FOR C (GAMMA-DOT)=ZERO. DO 106 1=1,N IF (. NOT. LREV (I) ) GO TO 106 IF (ABS (GAMMAO (I) ) .GT. 1. E-8) GO TO 120 NREVI =10 GO TO 121 120 ' TEMY = GD1P (I) TEMP = TEMY-GD1 (I) IF (ABS (TEMP) .GT. 1. E-20) GO TO 133 IF (LPRT) WRITE (6,3) I NCR , ITERNO, I 3 FORMAT(12X,'INCR =»,I4,», ITERNO =',13,', LAYER',13, ', (GD1-GD1P)=0«) NREVI =10 GO TO 121 133 TEMY = 10.0*TEMY/TEMP NREVI = IFIX (TEMY) IF (NREVI.LE.O) NREVI=0 IF (NREVI.GE.10) NREVI=10 121 CONTINUE NREV (I) = NREVI 106 CONTINUE C DO 107 KNT=1,11 K = KNT-1 DO 108 1=1,N IF (.NOT.LREV(I)) GO TO 103 IF (NREV (I) .EQ*. K) GO TO 109 108 CONTINUE 107 CONTINUE 109 NDT = K IF (NDT.GT.NDDT) NDT=NDDT C IF (NDT. NE. 0) GO TO 112 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 18 190 DO 113 1=1,N LRG (I) = - FALSE. GAMMA (I) = GAMMAO(I) IF (LRGL (I) ) GO TO 113 IF (. NOT. LREV (I) ) GO TO 11 3 IF (NREV (I) . NE. 0) GO TO 113 LRG (I) = .TRUE. LRGL (I) = .TRUE. DIRN(I) = -DIRN(I) IF (LPRT) WRITE(6,4) INCR,ITERNO,I,NDT 4 FORMAT(12X,«INCR =«,I4,', ITERNO =',13,»: ',13, ,» CHANGES MOD. , NDT =«,I2) GAMR (I) = GAMMAO (I) TAUR (I) = TAXYP (I) 113 CONTINUE CALL SYSCKM(O) NRETN = 3 LIT = .TRUE. RETURN 112 IF(ITERNO.EQ.1.AND.NDT.EQ.10) GO TO 134 FN = FLOAT (NDT) UG1 = UG1-UGII DDT = 0.1*FN*DT UGII = 0.1*FN*UGI UG1 = UG1+UGII LIT = .TRUE. CALL INTSCH (DDT,UG1,UGII) C 134 DO 116 1=1,N LRG (I) = . FALSE. IF (LRGL (I) ) GO TO 1 16 IF (. NOT. LREV (I) ) GO TO 116 IF (NPEV(I) .EQ.NDT) GO TO 117 TEMY = GAMMA (T)-GAMMAO (I) IF (TEMY. EQ. 0.0) GO TO 116 TEMY = DIRN (I) *TEMY IF (TEMY. GT. 0.0) GO TO 116 117 LRG (I) = .TRUE. LRGL (I) = . TRUE. IF (LPRT) WRITE (6,4) INCR,ITERNO,I,NDT IF (LPRT) WRITE (6,5) I, ABXD2 (I) , XD1 (I) ,XDO (I) , GAMMA (I) , TAXY (I) ,EVP (I) ,PP (I) 5 FORMAT(20X,I4,1P7E12.3) DIRN (I) = -DIRN (I) GAMR(I) = GAMMA (I) TAUR (I) = TAXY (I) 116 CONTINUE NTDT = NTDT+NDT NDDT = 10-NTDT FN = FLOAT(NDDT) DDT = 0.1*FN*DT UGII = 0.1*FN*UGI LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 19 191 DG1 = UG1+UGII NRETN = 2 IF (NTDT.GE.10) NRETN=4 RETURN 131 NRETN =1 DO 132 1 = 1 , N . 132 LRGL (I) = .FALSE. RETURN END SUBROUTINE INTSCH(DT,UG1,UGI) LOGICAL*1 LJUMP,L2,L3,L4,LDOF1,LLIN,LSTOP,LMCLG, LJRV,LCKFB,LIT,LPMAX COMMON . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX . /PROPTY/TITLE(20) ,PROP(20,17) ,H(20) ,DEPTH(20) ,WIDTH(20) , W T(20) , MATYP (20) ,SIGZ (20) ,CC (20) ,GG (20) ,NSUBD (20) ,GAMR ( 2 0 ) , TAUR(2 0) ,CM (40) , SK (40) , GMN (2 0) , TMN (20) . /DYNEQT/DYM (40) ,PEP (20) ,FORER(20) . /CALVAL/DIRN(20) , X D 2 ( 2 0 ) , A B X D 2 ( 2 0 ) , X D 1 (20) , X D 0 (20) ,BD (2 0 ) , BA (20) , B B (20) , G D 2 (20) , G D 1 (20) , G AMMA (20) , GAMMAO ( 2 0) , GAMMAX(20),TAXY(20) , EVP ( 2 0 ) , S K 1 R ( 2 0 ) , D P P (20) , P P ( 2 0) . /CPVALU/BAP (20) ,BBP(20) , X D 2 P ( 2 0 ) , X D 1 P ( 2 0 ) , X D 0 P ( 2 0 ) , G D 2 P ( 2 0 ) , G D 1 P (20) ,TAXYP (20) . /REVPAR/NRPT(2 0 ) ,GGAMR (20,3) ,TTAUR ( 2 0,8) , DIRNR (20,8) DIMENSION FCC E M (20) , FOCEC (20) , FOCEK (20) ,FORERO(20) DATA LJUMP/.FALSE./ IF (LJUMP) GO TO 100 LJUMP = .TRUE. DFINT =0.5 AFINT = 0.16666667 IF (INTYP. EQ. 1) AFINT=0.25 100 CONST 1 = DFINT*DT C O N S T 2 = AFINT*DT*DT C FORM RIGHT HAND SIDE OF DYNAMICAL EQUATION. DO 1 0 1 1 = 1 , N TEM1 = 0.0 TEM2 = 0.0 DO 102 IJ=1,3 II = I+IJ-2 IK = 2* (.1-1) +IJ-1 IF (IK.LE.O) GO TO 103 IF (IK.GE.2*N) GO TO 103 TEM 1 = T E M 1+CM (IK) *BAP (II) T E M 2 = T E M 2 + SK(IK) *BBP(II) 103 CONTINUE 102 CONTINUE PEP(I) = -WT (I) * U G I - T E M 1 - T E M 2 IF (.NOT. LCKFB) GO TO 1 0 1 IF (.NOT.LIT) GO TO 1 2 1 PEP (I) = PEP (I)-FORERO (I) GO TO 1 0 1 121 PEP (I) = PEP(I)-FORER (I) LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 20 192 101 CONTINUE C FORM LEFT HANDSIDE OF DYNAMICAL EQUATION DO 10S 1=1,N I J = 2*1-1 DYH (IJ) = WT (I) +CONST1 *CM (IJ) +CONST2*SK (IJ) I J = IJ+1 DYM(IJ) = CONST1*CM(IJ)+CONST2*SK (IJ) 105 CONTINUE C CALL SOLB2 C C CALCULATE ACC., VEL.,DISPL., GAMMA AND TAXY ETC. C DO 106 1=1,N TEH 1 = PEP (I) XD2 (I) = XD2P (I) + TEM1 ABXD2 (I) = XD2 (I) +UG1 XD1 (I) = XD1P (I) + BAP (I) +TEM1*C0NST1 XDO (I) = XDOP (I) + BBP (I) +TEM1 *CONST2 BA (I) = XD2 (I) *DT BB(I) = XD1 (I) *DT+BA (I) *CONST1 C IF (I.NE. 1) GO TO 107 TEM 1 = 1./H(1) GD2 (1) = XD2 (1) *TEM 1 GD1 (1) = XD1 (1) *TEM 1 GAMMA (1) = 100. 0*XD0 (1) *TEM1 GO TO 106 107 TEM 1 = 1./H(I) GD2(I) = (XD2 (I) -XD2 (1-1) ) *TEM1 GD1 (I) = (XD1 (I)-KD1 (1-1) ) *TEM1 GAMMA (I) = 100. 0* (XDO (I)-XDO (1-1) ) *TEM1 106 CONTINUE DO 114 1=1,N GMNT = GMN (I) TMNT = TMN (I) IF (LLIN) GO TO 116 J = NRPT (I) IF (J.EQ.O) GO TO 115 TEM 1 = 0. 01* (GAMMA (I) -GGAMR (I, J) *TMNT/GMNT) TAXY (I) = TTAUR (I, J) *TMNT +GMNT*TEM1/(1.0+0.5*GMNT*ABS(TEM1)/TMNT) GO TO 114 115 TEM 1 = 0.01*GAMMA(I) TAXY (I) = GMNT*TEM1/(1.0 + GMNT*ABS (TEM 1) /TMNT) GO TO 114 116 TAXY (I) = 0.01 *GMNT*GAMMA (I) 114 CONTINUE C DO 111 1=1,N FOCEM (I) = WT (I) *ABXD2 (I) FOCEC(I) = CC (I) *WIDTH (I) *GD1 (I) FOCEK(I) = TAXY (I) *WIDTH (I) LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 21 193 IF (I.EQ.N) GO TO 113 11 = 1+1 FOCEC(I) = FOCEC (I) -CC (11) *WIDTH (11) *GD1 (11) FOCEK (I) = FOCEK (I) -TAXY (11) *WIDTH (11) 113 CONTINUE TEM 1 = FOCEH(I) TEM2 = FOCEC(I) TEM3 = FOCEK (I) TEM4 = TEM1+TEM2+TEM3 FORERO(T) = FORER (I) FORER (I) = TEM4 TEM4 = ABS(TEM4) TEM5 = ABS (TEM1) +ABS (TEM2) +ABS (TEM3) IF (TEM4.LE.TEM5) GO TO 111 LSTOP = .TRUE. WRITE (6,2001) INCR,.I,GDI (I) ,GD2 (I) ,TEM 1, TEM2, TEM 3 , FORER (I) 2001 F0RMAT(1X,2I4,2X,1P6E14.4) 111 CONTINUE IF (LSTOP) STOP RETURN END SUBROUTINE PPGEN(TIME,LCPPG,LSTOP) LOGICAL* 1 LVDIR (20) ,LREV (20) , LRG (2 0) , DODEVP (20) , LLIQ (20) , LLPR(2 0) , LUPMP (20) ,LUPMC (20) ,LNUPK (20) ,LXZY (20) , LPRINT,LCALL,LCPPG,LSTOP COMMON . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB /LOGCOM/LVOIR,LREV,LRG,DODEVP,LLIQ . /PROPTY/TITLE(20) ,PROP (20,17),H (20),DEPTH (20),WIDTH (20) ,WT (20) , MATYP (20) ,SIGZ (20) ,CC (20) ,GG (20) ,NSUBD (20) ,GAMR (20) , TAUR (20) ,CM (40) , SK (40) ,GMN (20) ,TMN (20) . /CALVAL/DIRN (20) ,XD2 (20) ,ABXD2 (20) ,XD1 (20) ,XD0 (20) ,BD(20) , BA (20) ,BB (20) ,GD2 (20) ,GD1 (20) ,GAMMA (20) ,GAMMAO (20) , GAMMAX (20) ,TAXY (20) ,EVP(20) ,SK1R(20) ,DPP (20) ,PP (20) . /CPVALU/BAP (20) ,BBP(20) ,XD2P(20) ,XD1P(20) ,XD0P(2O) ,GD2P(20) , GD1P (20) ,TAXYP (20) DIMENSION YMXP (20) ,DEVP (20) ,DEVPP (20) ,DPPP (20) ,TFYMP (20) , DEVPT(20) ,YMXT(20) ,VCC1 (20) ,RDENK (20) DATA LCALL/.FALSE./ C , c ********** C INITIALIZATION FOR FIRST CALL C NOTE PP(I) AND EVP (I) HAS TO BE INITIALIZED IN THE MAIN PROGRAM C ********** IF (LCALL) GO TO 105 LCALL = .TRUE. NPRINT = 0 LPRINT = .TRUE. NWLIQ =0 VCC2 = 0.790 VCC3 = 0.5625 VCC4 = 0.730 POWTM = 0.43 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 22 194 POWTN = 0.620 SK2T = 0.0025 POWTMC = 1.-POWTM DO 106 1=1,N YMXT(I) = GAMMA (I) BD(I) = GAMMA (I) YMXP (I) = O.'O TFYMP(I) =0.0 DPPP(I) =0.0 LUPMP(I) = .TRUE. VCC1 (I) = 0. 5*PROP (1,7) TEMP = POWTN-POWTM TEMP = POWTM*SK2T*SIGZ (I)**TEMP RDENK(I) = 100.0/TEMP SK1R(I) = RDENK (I) *SIGZ (I) **POWTMC 106 LLPR(I) = • FALSE. RETURN C Q ********** C CHECK GAMMA FOR CROSSING ZERO VALUE AND DETERMINE GAMMAX C LXZY = .TRUE. FOR GAMMA CROSSING ZERO STRAIN AXIS C LUPMP = .TRUE. FOR GAMMA INCREASING IN MAG. DURING LAST INTERVAL C LUPMC = .TRUE. FOR GAMMA INCREASING IN MAG. DURING CUR. INTERVAL C LNUPK = .TRUE. FOR NEW PEAK OCCURING DURING CURRENT INTERVAL C ********** 105 CONTINUE LCPPG = .FALSE. DO 122 1=1,N DPP (I) =0.0 LUPMC (I) = .FALSE. DODEVP (I) = .FALSE. LXZY (I) = .FALSE. TEMP = BD (I) TEMC = GAMMA (I) TPRD = TEMP*TEMC IF (TPRD) 102,103, 101 C C GAMMAO AND GAMMA ARE ON THE SIDE, CHECK FOR GAMMAX 101 TEMPI = ABS(YMXT(I)) IF (ABS (TEMC) .LE.TEMPI) GO TO 131 YMXT(I) = TEMC IF (ABS (GAMMAX (I)) . LT. ABS (YMXT (I) ) ) GAMMAX (I) =YMXT (I) 131 IF (ABS (TEMC) .GT. ABS (TEMP) ) LUPMC (I) =. TRUE. GO TO 104 C 102 CONTINUE C GAMMA CHANGES SIGN HERE - I.E. CROSSES ZERO, AXIS, STOP PPGEN LXZY (I) = .TRUE. GO TO 104 C 103 CONTINUE C EITHER GAMMA OR GAMMAO OR BOTH EQUAL 0 IF (TEMP. NE. 0.0) GO TO 102 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 23 195 LOPHC(I) = .TRUE. GO TO 121 104 LNUPK(I) = .FALSE. IF (LUPMP (I) . AND. . NOT. LUPMC (I) ) LNUPK (I) = . TRUE. IF (LUPMC (I)) GO TO 121 LCPPG = .TRUE. DODEVP (I) = '.TRUE. 121 LUPMP(I) = LUPMC (I) 122 CONTINUE C c ********** C CHECK FOR PRINTING GAMMAX C ********** IF {.NOT. LCPPG) GO TO 888 DO 107 1=1,N IF (LLPR (I) . AND.LXZY (I) ) GO TO 109 107 CONTINUE GO TO 108 109 IF (. NOT- LPRINT) GO TO 108 NPRINT = NPRINT*1 IF (NPRINT.GT.20) LPRINT=.FALSE. WRITE (6,5) INCR 5 FORMAT('0 INCR =',I4,» GAM-MAX VOL CHANG PRES CHAN', " TIME •) DO 110 1=1, N IF (. NOT. LLPR (I) ) GO TO 110 WRITE (6,6) I,YHXP(I) ,DEVPP(I) ,DPPP(I) ,TFYMP(I) 6 FORMAT (1 OX, 14 ,1P3E 12. 3, 0PF8. 2) LLPR (I) = -FALSE. 110 . CONTINUE 108 CONTINUE C c ********** C CALCULATE VOLUME CHANGE OR PORE PRESSURE INCREMENT Q ********** DO 130 1=1,N DEVP(I) =0.0 DPP (I) =0.0 IF (.NOT.DODEVP(I)) GO TO 130 IF (. NOT. LNUPK (I) ) GO TO 116 TEMC = EVP (I) TEMP = ABS (YMXT (I) ) IF (TEMP.GT. 1.E-10) GO TO 115 114 DODEVP(I) = .FALSE. YMXT (I) =0.0 DEVPT(I) =0.0 GO TO 130 115 TPRD = VCC1 (I)*(TEMP-0.79*TEMC+0.5625*TEMC*TEMC/(TEMP+0.73*TEMC) ) IF (TPRD.LT.0.0) GO TO 114 DEVPT (I) = TPRD 116 IF (LUPMC(I)) GO TO 130 TEMP = ABS (YMXT (I) ) IF (TEMP.LT. 1.E-10) GO TO 130 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 196 TEMC = ABS (BD (I) -GAMMA (I) ) IF (LXZY (I)) TEMC=ABS (BD (I) ) DEVP(I) = TEMC*DEVPT (I)/TEMP EVP (I) = EVP (I) +DEVP (I) IF (NVOL.EQ.O) GO TO 117 ' . • , DPP (I) = DEVP(I) PP (I) = EVP (I) GO TO 118 C C CALCULATE PORE PRESSURE INCREMENT 117 SGVP = SIGZ (I)-PP (I) IF (SGVP.GT.0.05*SIGZ (I) ) GO TO 112 SGVP = 0.05*SIGZ (I) IF (LLIQ (I) ) GO TO 112 LLIQ(I) = -TRUE. WRITE ( 6 , 3 ) INCR, I 3 FORMAT(' **** AT INCR =',I4,» ; LAYER',13,' LIQUEFIED ! ****') NWLIQ = NWLIQ+1 N02 = N/2 IF (NWLIQ.LE.N02) GO TO 112 WRITE(6,200U) 2004 FORMAT('0** HALF OF DEPOSIT HAS LIQUEFIED ! **') LSTOP = .TRUE. GO TO 888 112 DEVPS = 0.01*DEVP(I) TEMC = RDENK(I)*SGVP**POWTHC DPP(I) = TEMC*DEVPS SK1R(I) = TEMC 118 CONTINUE IF (. NOT. LXZY (I) ) GO TO 130 LLPR(I) = .TRUE. YMXP(I) = GAMMAX (I) GAMMAX (I) = 0.0 YMXT(I) = GAMMA(I) DEVPP (I) = DEVPT (I) DPPP(I) = DPP (I) 130 CONTINUE 888 DO 889 1=1,N 889 BD(I) = GAMMA (I) RETURN END SUBROUTINE PPDISN(DT,NSUBDT,LCPPG) C SOLUTION OF ONE-DIMENSIONAL HEAT OR CONSOLIDATION PROBLEMS. C MAXIMUM 20 LAYERS C FINITE DIFFERENCE SOLUTION C PROP (I,J) = PROPERTIES OF SOIL LAYER FROM MAIN PROGRAM C SK1R(I) = ONE-D REBOUND MODULUS FOR THE LAYER C H(I) = THICKNESS OF SOIL LAYER C NSUBD (I) = NUMBER OF SUBDIVISIONS IN THE LAYER C PPO (I) = PORE PRESSURES AT THE START OF TIME INTERVAL C DPP (I) = INCREMENTS OF PORE PRESSURE FOR THE TIME INTERVAL C REFERRED TO MID-LAYER C DT = TIME INTERVAL FOR ITERATION LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 197 C N — NUMBER OF SOIL LAYERS C NTB,NBB = 0 FOR W=0 AT ALL TIMES AT TOP OR BOTTOM BOUNDARY, C 1 FOR (DW/DN) =0 AT ALL TIMES. C DZ(I) = DEPTH INTERVAL C CV(I) = COEFFICIENTS OF CONSOLIDATION C C AM (I) C.V*DT/(DZ (I) **2) C 101 130 131 132 133 C 15 102 103 49 C C LOGICAL*1 LVDIR (20) ,LREV(20) ,LRG (20) ,DODEVP(20) ,LLIQ(20) , LDOF1,LLIN,LSTOP,LMCLG,LCPPG,LJUMP,L3,L4 COMMON /LCONTL/LDOF1,LLIN,LSTOP,LMCLG /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB /LOGCOM/LVDIR,LREV,LRG,DODEVP,LLIQ /PROPTY/TITLE(20) ,PROP (20,17) ,H (20) ,DEPTH (20) ,WIDTH (20) ,WT (20) HATYP (20) ,SIGZ (20) ,CC (20) ,GG (20) ,NSUBD (20) ,GAMR (20) , TAUR (2 0) ,CM (40) ,SK (40) ,GMN (2 0) ,TMN (20) /CALVAL/DIRN (20) ,XD2(20) , ABXD2(2 0) ,XD1 (20) ,XD0 (20) ,BD(20) , BA(20) ,BB(20) ,GD2 (20) ,GD1 (20) , GAMMA (20) , GAMMAO (20) , GAMMAX (2 0) ,TAXY (20) ,EVP (20) , SK1R (20) , DPP (20) ,PP (20) DIMENSION W (100) ,TW (100) ,CV (20) , AM (20) ,DZ (20) DATA LJUMP/.FALSE./, W/100*0.0/, TW/100*0.0/ IF (LJUMP) GO TO 25 LJUMP = .TRUE. WRITE(6,101) FORMAT (•0*** DESCRIPTION OF CONSOLIDATION PROBLEM ***«/) IF (NBT. EQ. 0) WRITE(6,130) IF (NBT. EQ. 1) WRITE(6,131) IF (NBB. EQ. 0) WRITE(6,132) IF (NBB. EQ. 1) WRITE(6,133) FORMAT (' AT THE TOP BOUNDARY THE CONDITION IS ( W=0 )') FORMAT (' AT THE TOP BOUNDARY THE CONDITION IS ( DW/DN = 0 )') FORMAT (' AT THE BOTTOM BOUNDARY THE CONDITION IS ( W=0 )•) FORMAT (' AT THE BOTTOM 30UNDARY THE CONDITION IS ( DW/DN = 0 )' NTOT IS THE TOTAL NUMBER OF SOIL POINTS FOR PP CALCULATION NTOT = 0 DO 15 1=1,N IF (NSUBD (I) . EQ. 0) NSUBD(I)=2 FN = FLOAT(NSUBD(I)) DZ (I) = H (I) /FN NTOT = NTOT+NSUBD(I) CONTINUE NTOT = NTOT+1 WRITE (6, 102) NTOT FORMAT (» TOTAL NUMBER OF SOIL POINTS FOR PP CALCULATION =',I4) IF (NTOT.LE.100) GO TO 49 WRITE(6,103) FORMAT ('0*** NTOT IS GREATER THAN 100, PROCESSING STOPS.***') STOP CONTINUE INITIALIZING THE VARIABLES DO 16 IJ=1,NTOT LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 26 198 W (IJ) = 0. 0 TW(IJ) = 0.0 16 CONTIN.UE WRITE (6,105) 105 FORMAT('0*** INITIAL PROPERTIES OF SOIL LAYERS ***•/) C C CHECK THAT ALL AM'S ARE LESS 0.5 LSTOP = .FALSE. IF (NSUBDT.LE.10) GO TO 201 WRITE (6,1 35) NSUBDT 135 FORMAT(•0***** NSUBDT IS TOO LARGE !! =•,14,' ***•*•) STOP 201 RFNDT = FLOAT (NSUBDT) DDT = 0.01605*DT/RFNDT 50 DO 14 1=1,N TEMP = DZ (I) TEMP = TEMP*TEMP AH (I) = DDT*PROP (I, 15) *SK1R (I)/TEMP IF (AM (I).LT.0.5) GO TO 14 WRITE(6,109) I 109 FORMAT (•0*** AM VALUE IS TOO HIGH FOR LAYER',13,' ***•) LSTOP = .TRUE. 14 CONTINUE WRITE(6,107) 107 FORMAT ('0 LAY.NO. K-S K1R-S THICK. •AM-S NSUBD-S'/) WRITE (6, 108) (I, PROP (I, 15) ,SK1R(I) ,H(I) ,AM(I) , NSUBD (I) ,I=1,N) 108 FORMAT (1X,I6,4X,1P2E12.3,0PF12.2,1PE12.3,3X,0PI6) IF (LSTOP) STOP RETURN C C PORE PRESSURE REDISTRIBUTION STARTS HERE 25 CONTINUE DO 37 K=1,NSUBDT NP = 0 DO 26 1=1,N NN = NSUBD(I) TEMP = DZ (I) TEMP = TEMP*TEMP AMI = DDT*PROP (1,15) *SK1R (I) /TEMP IF (AMI.LT.0.5) GO TO 17 LSTOP = .TRUE. GO TO 50 17 CONTINUE DO 26 IJ=1,NN NP = NP+1 IF (I.NE. 1) GO TO 19 IF (IJ.NE. 1) GO TO 21 C CALCULATION FOR BOTTOM BOUNDARY IF (NBB. EQ. 0) GO TO 18 TW(1) = 2. *AMI* (W (2)-W (1) )+W (1) GO TO 18 19 IF (IJ. NE. 1) GO TO 21 LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 27 199 C CALCULATION FOR INTERFACE OF 2 LAYERS TEM 1' = DZ (I) /SK1F. (I) TEMP = DZ (1-1)/SK1R (1-1) TEMP = TEMP+TEM1 CONST = 2.*DDT/TEMP A = PROP(I-1.,15)/DZ(I-1) B = PROP (I, 15) /DZ (I) TW(NP) = CONST* (A*W (NP-1) + B*W (NP+1) - (A + B) *W (NP) )+W (NP) GO TO 18 C CALCULATION FOR POINTS IN A LAYER 21 TW(NP) = AMI*(W(NP-1) + W (NP+1)-2. *W (NP) ) + W (NP) 18 CONTINUE 26 CONTINUE C CALCULATION FOR TOP BOUNDARY IF (NBT. EQ. 0) GO TO 27 NP = NP+1 TW(NP) = 2. *AMI* (W (NP-1)-W (NP) )+W (NP) 27 CONTINUE C C ADD DPP TO HT(I) TO GIVE W (I) NP = 0 DO 30 1=1,N DPI = DPP (I) *RFNDT NN = NSUBD (I) .IF (T.EQ. N) NN=NN+1 DO 30 IJ=1,NN NP = NP+1 IF (NP.EQ.LAND.NBB.EQ.O) GO TO 33 IF (NP.EQ.NTOT.AND. NBT. EQ.O) GO TO 33 IF (I.NE. 1. AND.IJ. EQ. 1) GO TO 32 C NORMAL SOIL POINTS, BOTTOM BOUND. PT. AND TOP BOUND. PT. W(NP) = TW(NP)+DPI GO TO 33 C INTERFACE 32 W (NP) = TW (NP) +0 . 5 * (DPI + DPP (1-1) ) 33 CONTINUE 30 CONTINUE 37 CONTINUE C FIND PP(I) WHICH IS AVERAGE OF W (I) FOR A LAYER NP = 1 DO 35 1=1,N NP = NP-1 NN = NSUBD (I) +1 FN = FLOAT (NN) A = 0.0 DO 36 IJ=1,NN NP = NP+1 36 A = A+TW (NP) 35 PP(I) = A/FN C WRITE (6,9999) (TW (K) , K= 1 , NP) C999 FORMAT (' OWT (HP) : 1 , 1P5E14. 3, (/9X, 1P5E1 4. 3) ) C WRITE (6,9998) (DPP (K) , K= 1 , N) C998 FORMAT (' DPP (I) :»,1P5E14.3, (/9X,1P5E14.3)) LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 200 C WRITE (6, 9997) (W (K) , K=1 , NP) C997 FOR MAT(' W(NP) :»,1P5E14.3, (/9X,1P5E14.3)) C WRITE (.6,9996) (PP (I) ,T= 1, N) C996 FORMAT (' PP (I) :», 1P5E14.3, (/9X,1P5E14.3)) RETURN END SUBROUTINE MAXVAL(TIME,UG1) LOGICAL*1 LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX COMMON . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /CALVAL/DIRN (20) ,XD2 (20) ,ABXD2 (20) ,XD1 (2 0) ,XD0 (20) ,BD (20) , . BA(20) ,BB(20) ,GD2(20) ,GD1 (20) , GAMMA (20) ,GAMMAO(20) , GAMMAX (20) ,TAXY (20) , EVP (20) ,SK1R (20) ,DPP (20) ,PP (20) DIMENSION ACCX(20) ,GAMX(20) ,TMAX(20) ,EVPX(20) ,PPMX(20) , TFM (20,5) DATA ACCX/20*0.0/, GAMX/20*0.0/, TMAX/20*0.0/, EVPX/20*0.0/, PPMX/20*0.0/, TFM/100*0.0/, ACCINX/0.0/, T1/0.0/, TP/0.0/ IF (LPMAX) GO TO 200 TEM = ABS (UG1) IF (ACCINX.GE.TEM) GO TO 120 ACCINX = TEM T1 = TIME 120 CONTINUE C DO 150 1=1,N TEM = ABS (ABXD2(I)) IF (ACCX (I) .GE.TEM) GO TO 101 ACCX(I) = TEM TFM (1,1) = TIME 101 TEM = ABS (GAMMA (I) ) IF (GAMX (I) .GE.TEM) GO TO 102 GAMX (I) = TEM TFM (1,2) = TIME 102 TEM = ABS (TAXY (I) ) IF (TMAX (I) .GE.TEM) GO TO 103 TMAX(I) = TEM TFM (1,3) = TIME 103 TEM = ABS (EVP (I) ) IF (EVPX (I) .GE.TEM) GO TO 104 EVPX(I) = TEM TFM (1,4) = TIME 104 TEM = ABS (PP (I) ) IF (PPMX(I) .GE.TEM) GO TO 150 PPMX(T) = TEM TFM (1,5) = TIME 150 CONTINUE RETURN C C PRINT OUT MAXIMUM VALUES BEFORE TERMINATION 200 CONTINUE WRITE (6,2001) TP,TIME 2001 FORMAT ( LISTING OF LIQUEFACTION PROGRAM - AUG, 75 LIQPG 29 201 . »0********************************************* / • * MAXIMUM VALUES OCCURRED -',F8.3,' TO',F8.3,' SEC *'/ • ******************************************* ' LAY. TIME ABS ACC -TIME STRAIN TIME ', •STRESS TIME VOL STRAIN TIME PORE PRESS'/) DO 2000 1=1,;N 2000 WRITE (6,2002) I, TFM (1, 1) , ACCX (I) ,TFM (I , 2) , GAMX (I) ,TFM (I, 3) , TM AX (I) , TFM (1,4) ,EVPX(I) ,TFM (1,5) ,PPMX (I) 2002 FORMAT (14,1X,5 (0PF7. 2, 1PE11. 3) ) WRITE (6,2003) ACCINX,T1 2003 FORMAT('0** MAX BASE ACC =',1PE11.3,' : TIME =',0PF5.2, ' SECS **') C C RESET VARIABLES TO ZERO ACCINX = 0,0 DO 160 1=1,N ACCX(I) =0.0 GAMX(I) = 0.0 TMAX(I) = 0.0 EVPX(I) =0.0 PPMX(I) =0.0 160 CONTINUE TP = TIME RETURN END EXECUTION TERMINATED $SIG
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Mechanical model for the analysis of liquefaction of...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Mechanical model for the analysis of liquefaction of horizontal soil deposits Lee, Kwok Wing 1975
pdf
Page Metadata
Item Metadata
Title | Mechanical model for the analysis of liquefaction of horizontal soil deposits |
Creator |
Lee, Kwok Wing |
Date Issued | 1975 |
Description | During the development of liquefaction in a soil deposit subjected to vibration there are two processes which work in opposite directions. The volume compaction tendency under cyclic loading causes the pore water pressure to rise, and the dissipation of excess pore water pressure (consolidation) decreases it. Recently, Martin, Finn and Seed (l975) studied the mechanics of pore water pressure generation of a soil sample subjected to cyclic loading and a relation between shear strain cycles, volume compaction and pore water pressure increment was established. A material model based on this relationship is developed in this thesis for saturated granular soil under cyclic simple shear conditions. The model includes a hysteretic stress-strain relationship, volume compaction, pore water pressure rise and dissipation. Using this proposed comprehensive material model, a global mechanical model is constructed to simulate the liquefaction (including consolidation) behavior of a thick horizontal deposit when subjected to horizontal base motion. In this way the coupled problems of dynamic response, pore water pressure rise and consolidation of the deposit under seismic loading can be analysed. The numerical techniques used to solve such problems are discussed in detail. The response of a typical saturated sand deposit under earthquake loading is determined using the proposed model and the results show that the model can predict the various phenomena that saturated sand deposits exhibit during earthquakes. The global model also makes clear the influence of permeability on the liquefaction potential of the soil deposit. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0063017 |
URI | http://hdl.handle.net/2429/19649 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1975_A1 L43_5.pdf [ 8.54MB ]
- Metadata
- JSON: 831-1.0063017.json
- JSON-LD: 831-1.0063017-ld.json
- RDF/XML (Pretty): 831-1.0063017-rdf.xml
- RDF/JSON: 831-1.0063017-rdf.json
- Turtle: 831-1.0063017-turtle.txt
- N-Triples: 831-1.0063017-rdf-ntriples.txt
- Original Record: 831-1.0063017-source.json
- Full Text
- 831-1.0063017-fulltext.txt
- Citation
- 831-1.0063017.ris
Full Text
Cite
Citation Scheme:
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>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0063017/manifest