A NEW, NON-LINEAR ANALYSIS OF LAYERED SOIL-PILESUPERSTRUCTURE SEISMIC RESPONSE by ANUZ K . K H A N B. E.,(Hons.), Bharathiar University, India, 1991 Diploma in Civil Engineering, Government Polytechnic, India, 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF M A S T E R OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming to the requirecLstandard THE UNIVERSITY OF BRITISH C O L U M B I A AUGUST 1995 © Anuz K. Khan, 1995 In presenting degree freely at the available copying of department publication this of in partial fulfilment of the University of British Columbia, I agree for this or thesis reference thesis by this for his thesis and study. scholarly or her for I further purposes financial gain shall permission. Department of £ \ \ M L B^&\ The University of British Columbia Vancouver, Canada Date DE-6 (2/88) OCT 14 , W ECRW/n that agree may representatives. requirements be It not that the be Library an by understood allowed advanced shall permission for granted is for the that without make it extensive head of copying my my or written ABSTRACT A non-linearfiniteelement analysis for determining the seismic response of piles in multilayered soil deposit is presented. The analysis incorporates pile material and geometric non-linearities, P-A effects, non-linear load deformation behavior of soil-pile interaction, pile separation and the non-linear shear stress-strain relationship in the freefield soil. The soil shear stress-strain relation is expressed based on Hardin and Drenvich's model. The pile material stress-strain relation is based on a elastic-plastic model. Near field soil is replaced with a continuous foundation based on Yan-Byrne or American Petroleum Institute P-y curves. Time integration is performed utilizing Newmark-P integration scheme. Numerical integration of the virtual work equations has been carried out with Gauss quadrature. The resulting set of non-linear equations is solved by using the Newton-Raphson scheme. The method of analysis are incorporated in two programs : QUIVER and PILEUBC. The program QUIVER captures the seismic response of soil deposit and later is used as a module in the main program PILEUBC, in which the response of a single pile is predicted for an earthquake excitation. The program predicts time varying acceleration, velocity, displacement, bending/axial stress and strain, bending moment and shear force in the pile. The program is able to successfully simulate the seismic response of piles. Numerical investigations have been carried out to verify the results and test the capabilities of the program. A linearized closed form solution based on Laplace transformation forfree-fieldis developed for comparison with thefiniteelement analysis. ii TABLE OF CONTENTS ABSTRACT ii TABLE OF CONTENTS iii LIST OF TABLES v LIST OF FIGURES ACKNOWLEDGMENTS vi viii DEDICATION ix C H A P T E R 1: 1 INTRODUCTION 1.1 Introduction 1.2 Research Objectives 1.3 Organization O f Thesis 1 6 7 C H A P T E R 2: 8 REVIEW OF SINGLE PILE ANALYSIS 2.1 Introduction 2.2 Review O f Single Pile Analysis 8 8 CHAPTER 3:.... 17 FINITE ELEMENT FORMULATION OF SOIL-PILE-SUPERSTRUCTURE INTERACTION 3.1.0 Introduction 3.1.1 Problem Formulation And Finite Element Discretization 3.1.2 Material Stress-Strain Relationships 3.1.3 Near Field and Pile Relationships 3.1.4 Pile Separation (Gapping) 3.1.5 Kinematical Relations and Shape Functions 3.1.6 Virtual Work Equations 3.1.7 Consistent Mass Matrix for Pile Element 3.1.8 Numerical Integration 3.1.9 Discretization Error and Convergence 3.1.10 P I L E U B C - The Computer Program 17 19 22 22 27 28 31 35 35 36 38 3.2 FINITE ELEMENT FORMULATION OF FREE-FIELD SOIL DEPOSIT 3.2.0 Introduction 3.2.1 Finite Element Discretization of Soil Deposit 3.2.2 Soil Shear Stress-Strain Relationships 38 38 39 39 iii 3.2.3 Computation of Hyperbolic Model Parameters 3.2.4 Problem Formulation 3.2.5 Shape Functions 3.2.6 Virtual Work Equations for Free-Field 3.2.7 Consistent Mass Matrix for Soil Element 3.2.8 Time and Numerical Integration 3.2.9 Discretization Error and Convergence 3.2.10 Q U I V E R - The Computer Program C H A P T E R 4: 44 46 47 48 50 50 52 53 54 VERIFICATION AND N U M E R I C A L RESULTS 4.1 Introduction 4.2 Verification of Free-Field Movement Against Developed Closed Form Solution 4.3 Results for Varying Gmax and x in a Heterogeneous Deposit 4.4 Results From an Analysis For C A N L E X Site, Richmond 4.5 Results From Finite Element Analysis O f Soil-Pile-Superstructure Interaction 4.6 Results of an Example Problem at C A N L E X Site 55 56 63 C H A P T E R 5: 81 CONCLUSIONS AND SCOPE FOR FURTHER R E S E A R C H 81 REFERENCES 83 APPENDIX I 88 f A P P E N D I X II 54 68 74 93 iv LIST O F T A B L E S T A B L E 4.1 Soil Propreties for Non-Linear Analysis of Example 1 60 T A B L E 4.2 Soil Properties for Non-Linear Analysis of C A N L E X Site 68 LIST O F F I G U R E S F I G U R E 3.1 General nature of soil-pile-superstructure system during seismic excitation F I G U R E 2.1 Pile model used in S P A S M analysis (Finn and Gohl, 1990) F I G U R E 2.2 Analytical model for p~y evaluation and mathematical model for 2 12 soil-pile springs (Kagawa and Kraft, 1980) 14 F I G U R E 3.2 Finite element discretization for pile 20 F I G U R E 3.3 Transformation of typical pile sections 21 F I G U R E 3.4 Elastic, perfectly plastic pile material 23 F I G U R E 3.5 Soil-pile interaction 25 F I G U R E 3.6 Soil deposit and finite element representation 40 F I G U R E 3.7 Loading-unloading curve for shear stress-strain modeling 42 F I G U R E 3.8 Modified hyperbolic shear stress-strain formulation (Byrne and Mclntyre, 1994) 43 F I G U R E 3.9 Virtual work equation for soil-deposit 49 F I G U R E 4.1 Comparison of finite element and closed form solution at the surface for uniform stress field 57 Absolute acceleration along the profile using 10 elements for uniform stress field 58 Absolute acceleration time history from Q U I V E R for simulated sinusoidal excitation 59 F I G U R E 4.2 F I G U R E 4.3 F I G U R E 4.4 Absolute acceleration time history from Q U I V E R for earthquake excitation 61 F I G U R E 4.5 Comparison of cyclic stress ratio Vs depth between Q U I V E R and S H A K E F I G U R E 4.6 62 Typical shear stress-strain plot from the finite element program...64 vi F I G U R E 4.7 F I G U R E 4.8 F I G U R E 4.9 F I G U R E 4.10 F I G U R E 4.11 SPT data from test data at No. 4 road at River Road, Richmond 65 Absolute acceleration time history from Q U I V E R for Richmond site (earthquake base excitation) 66 Absolute acceleration time history from S H A K E for Richmond site (earthquake base excitation) 67 Comparison between lab test results, C Y C P I L E (Vazhinkhoo, 1995) 70 Comparison of closed form and finite element solutions for a cantilever beam 71 F I G U R E 4.12 Relative displacement time history of pile for example problem .. .73 F I G U R E 4.13 Relative displacement profile of pile at selected time steps for example problem 1 Bending moment profile of pile at selected time steps for example problem 1 F I G U R E 4.14 F I G U R E 4.15 F I G U R E 4.16 F I G U R E 4.17 F I G U R E A . 1.1 75 76 Relative displacement time history of pile for example problem 2 78 Relative displacement profile of pile at selected time steps for example problem 2 79 Bending moment profile of pile at selected time steps for example problem 2 80 Soil deposit as a one dimensional shear wave propagation medium 89 F I G U R E A.2.1 Generalized flow chart for finite element program, P I L E U B C 95 F I G U R E A.2.2 Generalized flow chart for finite element program, Q U I V E R 96 ACKNOWLEDGMENTS I am greatly indebted to my supervisor Dr. Ricardo O. Foschi, who suggested this problem and so willingly shared his time and ideas in countless discussions during different stages of this thesis. His versatile experience in numerical analysis helped in successful completion of this project. I thank Dr. Peter M. Byrne who gave me a definite environment to work with and the discussions were of particular benefit to me. I am grateful to Dr. Don L. Anderson who gave his time in reading this thesis. Last, but farfromleast, I thank my parents to whom this thesis is dedicated; for their continued encouragement and many personal sacrifices to help me learn, I shall always be grateful. Anuz K Khan viii Inspired by The Earthquake When the Earth is shaken to her utmost convulsion And the Earth throws up Her burdensfromwithin, And man cries distressed: 'What is the matter with her ?' On that Day will she declare her tidings: For that thy Lord will have given her inspiration On that Day will men proceed in companies sorted out To be shown the Deeds that they had done. Then shall anyone, who has done an atom's weight of good, see it!. And anyone who has done an atom's weight of evil shall see it THE HOLY QURAN (Surah X C K ) to- tny <utt*Hi cutd aMa, ix CHAPTER i INTRODUCTION 1.1 INTRODUCTION Pile-supported structures and foundations are man's oldest method of construction for overcoming the difficulties of founding on soft soils. The use of piles dates back to the prehistoric period. However, up until the late nineteenth century the design was entirely based on experience, or even divine providence (Poulos and Davis, 1980). Significant damage to pile supported structures during major earthquakes such as Niigata and Alaska earthquakes of 1964 and the San Francisco earthquake of 1906 led to an increase in demand to reliably predict the response of piles. Since then, extensive research have been carried out and several analytical and numerical procedures (Matlock and Reese, 1960; Novak and Aboul-EUa, 1978; Kuhlemeyer, 1979; Kagawa and Kraft, 1980; Velez et. al., 1983; Wu, 1994) have been developed to determine the static and dynamic response of piles subjected to horizontal or vertical loads. Also, full scale experimental observations on the pile's behavior (Abe et. al., 1984; Blaney and O'Neill, 1986; Sy and Siu, 1992) and numerous model testing (Finn and Gohl, 1987, Yan and Byrne, 1991) have been carried out. The problem with which this thesis is concerned is schematically shown in Figure 1.1. It shows a pile, embedded in a soil deposit, and supporting a superstructure above ground level. The source of the earthquake motion is assumed to be at the base or till and this triggers a vibration of the soil deposit above. In turn, this forces the pile to move under the action of the vibrating soil. Although the actual problem is three dimensional, l Structural Motion Acceleration at Till Fig. 3.1 General Nature of SoU-PHe-Superstructure System During Seismic Excitation this thesis would only consider a two-dimensional approximation representing the effect of the soil on the pile by p-y curves relating the applied pressure to the relative movement between pile and soil. The region of applicability of this model is called the "near field". The results from most of the numerical analysis give an insight to the problem and indicate that the key features of the problem are now better understood. However, complexity of analysis, attribution of failure, damage observation difficulties, and the problems in conducting in-situ tests have hampered the routine methods of analysis for pile foundations subjected to earthquake loading. Also, it has long been recognized that nonlinear effects, such as nonlinear soil behavior, slippage, eventual gapping, and material non-linearity play a fundamental role in the response of piles to seismic excitations. A number of approaches have been formulated for the analysis of dynamic soil-pile interaction in the past years. The research work carried out in the area of seismic soil-pilefoundation-structure interaction could be most generally classified into the determination of kinematic seismic response, determination of pile-head impedance and determination of superstructure seismic response (Gazetas et. al., 1992). A coupled 3-D analysis of pile foundation and superstructure is generally considered not to be a feasible engineering option. As mentioned earlier, a simpler procedure based on substructuring is usually adopted. Soil-pile interaction analysis is conducted on the foundation alone to determine the pile foundation impedances. These impedances are then incorporated into the model of the superstructure and the dynamic response is calculated using ground motions at the ground surface as input (Luco, 1982). It is assumed that these motions are not affected by the pile itself. But this may not be a tenable assumption for earthquake loading, although it may prove sufficient for pile head loading. 3 There is a significant difference between seismic excitation of ground and foundation, and cyclic pile head excitation. Pile head excitation creates disturbances in the soil adjacent to the pile only, whereas the entire near-field and its interaction plays a vital role in a pile subjected to earthquake loading. The distinguishing characteristic of strong shaking is the far reaching effects of non-linearity. This type of problem, typical of pile response subjected to strong shaking under maximum probable earthquakes requires a nonlinear analysis. In addition, in most of the procedures for evaluation of pile response a linear or an equivalent linear analysis is usually done. The whole idea of evaluation of impedance functions exemplifies the methodology of analysis as either linear or equivalent linear to handle non-linearity in the problem, if any. For realistic computations of stressdeformation and ultimate behavior of pile-founded structures it becomes necessary to consider factors like non-linear behavior of soil, interaction effects between soil and structure and nature of loading. If a long pile and its foundation experience large displacements it may become necessary to consider geometric non-linearity in addition to material non-linearity. The first step in any soil-structure interaction seismic analysis is the calculation of the free-field response of the site, that is, the spatial and temporal variation of motion before excavating or rigging the soil and superimposing the structure. Several methods for evaluating the effect of local soil conditions on the ground response during earthquakes are presently available. One such is the widely accepted procedure incorporated in the program S H A K E (Scnabel etal., 1972). This method is based on the assumption that the main response in the soil deposit is caused by the upward propagation of shear waves from the underlying rock formation. The seismic response analysis taking into account the 4 effects of nonlinear nature of soil deformation can be made either by equivalent linear method (Seed and Idriss, 1969) or the step-by-step integration scheme. If the earthquake induced shear strains are large, say to the order of 1 percent, the accuracy of the equivalent linear solution fails to yield an acceptably good solution (Ishihara, 1982). For the problem of inducing larger shear strains in the soil deposits a step-by-step time integration procedure using a variable tangent stiffness must be used, even i f it requires a greater computational effort to obtain the solution in the time domain. At present, the state of practice in analyzing nonlinear lateral pile response represents the pile by beam-column elements and the soil by compliance springs. The nonlinear response of the soil is captured by using non-linear springs termed as p~y curves (API Code, 1987). These curves are specified for soft clay and sandy materials and are based on a hyperbolic tangent function. Yan and Byrne (1990) have developed parabolic p~y curves for sandy soils. Instead of using discrete spring system (lumped model) a nonlinear continuous foundation would provide an improved solution. It was observed during cyclic loading that pile deformation, vertical heave and soilpile separation (gaps) occur near ground surface. As a result of gapping, the lateral resistance of the soil-pile system is significantly reduced. When piles are closely spaced, gaps around each pile may overlap and the effect of group interaction may be very pronounced. It is difficult to incorporate all the major factors in a solution procedure, even in a numerical procedure. Numerical procedures such as the finite element method, can be used to handle many complexities. The resulting solution strategies are difficult from the view point of mathematical derivation and involve large computational effort. 5 . In general the solution procedure for soil-pile-superstructure interaction should consider the three dimensions of the problem. A two dimensional solution with interface springs connecting the pile to the free-field is considered a feasible engineering option. These solution techniques could be in the frequency domain utilizing an elastic or equivalent elastic approach. It could also be a time domain inelastic solution. At this point, a nonlinear solution using a continuum approach to the problem of lateral response of piles subjected to earthquake loading is deemed to be needed. In this thesis an attempt is made to obtain such a solution, based on a finite element formulation of the nonlinear soil-pilesuperstructure interaction, when subjected to seismic excitations. The effect of gapping is also considered. 1.2 R E S E A R C H O B J E C T I V E S This study aimed at the following objectives : 1. Development of a finite element analysis model for a reliable prediction of the seismic response of piles, based on beam on an inelastic continuous foundation approach using a step-by-step time integration procedure. This analysis will include the nonlinear behavior of the near-field soil. The finite element model will take into account the soil-pile interaction using p~y response curves to depict the nonlinear continuous foundation. Also, the analysis will introduce nonlinear material and geometric properties for the pile as well as account for the development of gaps in the pile-soil interface. It will be shown that the continuum based computational nonlinear approach for soil-pile-superstructure interaction is a viable method for estimating the seismic response. 2. Implementation of the analysis into a computer program which would allow flexibility in placing the layers of soil in different configurations. 6 3. Validation of this computer program by comparing its predictions with developed analytical and numerical results. 1.3 ORGANIZATION O F THESIS Chapter 2 provides a critical review of the current methods for estimating dynamic response of pile foundations. Emphasis is given to the study of single piles in a layered media. Chapter 3 presents the finite element formulation of the free-field soil and soil-pilesuperstructure interaction. The method of analysis to determine the solution to the problem is discussed. The implementation of computer program is described. The assumptions are pointed out and the advantages and limitations are criticized. Chapter 4 presents typical analysis of the study and example problems. Chapter 5 presents the results and validity of the analysis on mathematical basis. Summary and conclusions from the research are given in Chapter 6. 7 CHAPTER 2 REVIEW OF SINGLE PILE ANALYSIS 2.1 INTRODUCTION Like other dynamic analysis of civil engineering systems, those for single piles may be classified into two major categories: continuum models coupling the soil and pile as a unified system, and lumped mass-spring-dashpot models. Soil-pile interaction for dynamic loading conditions have been represented analytically by discrete models such as the subgrade reaction theory, by continuous models such as the elastic solutions ( Tajimi, 1969; Nogami and Novak, 1977), or by the semi continuous models such as the finite element method (Blaney, Kausel, and Rosset, 1976; Kuhlmeyer, 1979; Rosset and Angelides, 1979). Also Wu (1994) has developed solution procedures for three dimensional idealizations of the coupled soil-pile system. However, a coupled dynamic continuous foundation based approach, incorporating all the important features of the soil-pilesuperstructure interaction has not been developed yet. 2.2 REVIEW O F S I N G L E PILE A N A L Y S I S Linear dynamic analysis in which the soil and pile are considered as an elastic continuum was first formulated by Tajimi (1966) for lateral response of an end bearing pile. Significant contributions followed from Novak, Nogami and their co-workers (Novak, 1974; Nogami and Novak, 1977; Novak and Aboul-EUa, 1978). Their formulation was based on an elastic beam vibrating in a homogeneous or non- homogenous elastic isotropic continuum subjected to dynamic pile head loading. The main 8 assumptions in their model were: the soil was composed of a horizontal layers that are homogenous, isotropic, and linearly viscoelastic with material damping of a frequencyindependent hysteretic type. The soil properties were constant within each layer but may be different in individual layers. The pile was vertical, linearly elastic, of circular cross section that may vary stepwise at the interfaces of the layers, and it is bonded to the soil. If pile head lies above the grade or the pile is assumed to be separated from soil, the adjacent layers were modeled as void. The soil reactions acting on a unit length of pile were described by a complex soil stiffness (Novak et al, 1978) associated with vertical, horizontal, anti-symmetrical (rocking), and torsional displacements of the pile. Closed form expressions were developed for this linear coupled response analysis. This solution procedure offered an insight to the soil-pile system. However, it does not contribute towards the development of a solution for the problem of earthquake excitation where the strains in the free and near-field soil and in the pile itself are subjected to significant nonlinearities. In addition, the reduction of soil stiffness and the increase in damping associated with strong shaking are sometimes modeled crudely by making arbitrary changes in complex stiffness of soil. Hence these studies have not proved very useful for evaluating the response of pile foundations to earthquake loading. Gazetas and his co-workers (Krishnan et al, 1983; Gazetas, 1984; Makris et al, 1992) compiled available literature and analysis procedures to study the influence of kinematic interaction on differences in pile accelerations at the ground surface relative to that of free field. Their work suggested (Gazetas et al, 1992) that soil-foundationstructure interaction analysis under seismic excitation could be organized as follows: (1) Obtain the motion of foundation in the absence of superstructure. (2) Determine the dynamic impedances associated with swaying, rocking and cross swaying-rocking 9 oscillations of the foundation. (3) Compute the seismic response of the superstructure supported on the foundation impedances and subjected at its base to the foundation input motion that was obtained earlier. In general, the above procedure is referred as uncoupled analysis as it does not solve the soil-pile-superstructure vibration in an amalgamated form. For each step of the analysis several alternative formulations have been developed and published in the literature, including finite element formulations (Blaney et al, 1976; Gazetas, 1984), boundary element (Kaynia et al, 1982), semi-analytical and analytical solutions (Banerjee et al, 1977) and a variety of simplified methods (Nogami, 1985). The influence of kinematic interaction in an uncoupled superstructure analysis appears to be invalid, provided the free-field surface motions are dominated by relatively low amplitude vibrations. This in result could over predict the dynamic response of the system. Whitman and Bielak (Rosenbleuth, 1980) define soil-structure interaction as " If the motion at any point on the soil-structure interface differsfromthe motion that would occur at this point in thefreefield if the structure were not present, there is soilstructure interaction ". This is a very general and widely accepted definition. It exemplifies the significance of the behavior of the coupled soil-pile system at the interface or contact point. Hence the solution procedure discussed by Gazetas et al (1992) may not prove adequate as it fails to address this contact or no contact situation in the soil-pile system. In addition the motion of a point in free-field definitely differs from motion of the same point in the near-field. However, they concluded that the effects of kinematic interaction are relatively minor for excitation frequencies of up about 1.5 times that of fundamental frequency of the freefield. These conclusions were again based on a elastic linear dynamic analysis of soil-pile system. 10 Kuhlmeyer (1979) developed a three dimensional linear finite element analysis to represent the beam bending aspect of the problem. The results were similar to those of Novak, suggesting that for a linear solution the plain strain models can provide a reasonably adequate prediction of pile head behavior. The analysis procedures discussed up until now are essentially linear and are only adequate to predict minuscule vibrations. This is usually not the case in a strong motion problem. During an earthquake, piles may behave in a materially nonlinear manner, and a long pile in the case of an off-shore foundation may exhibit geometric non-linearities. Also, additional effects of soil non-linearity at high strain levels, pile separation, slippage and friction are predominant. It is relatively difficult to develop a solution procedure incorporating all the above mentioned behavior. Lumped mass-spring-dashpot models seemed to be an obvious choice to be employed in a nonlinear analyses based on an approach in which pile foundation and the superstructure are analyzed as a combined system. The interaction between the pile and the near-field soil are replaced using a series of nonlinear springs derived from full scale (Matlock et al, 1978) or model test measurements (Yan and Byrne, 1990) or non-linear finite element solutions. The stiffness of these springs represent the combined stiffness of the soil in the near-field and exterior free-field whose properties are governed by the intensity of shaking. Difficulties exists in incorporating all the nonlinear features at high strain levels in a continuum model for the coupled soil-pile-superstructure system. It could be that these difficulties stem in the formulation of a mathematical solution procedure and the computational effort required. In a lumped mass-spring-dashpot model, interaction between the pile and the nearfield soil is usually accounted for by the use of linear or non-linear lateral compliances (springs and equivalent viscous dashpots) placed along the length of the pile. These 11 compliance springs model the forces acting on the pile due to relative movement between the pile and soil. One particularly attractive method of analysis, because of its ability to simulate non-linear soil-pile interaction, is embodied in the computer program S P A S M (Matlock et al., 1978a-b). In using S P A S M , the pile is modeled as a linearly elastic beam column incorporating the flexural rigidity of the pile, EI, and the effects of spatially varying axial load, P(z), on pile bending (P-A effects) during lateral seismic loading. The model used in S P A S M is shown in Fig. 2.1. RlGIO CL-NDR-CAL MASS STATION a £ - 5 0 * . SURFACE MMT FREE FiELO OISPLACEMCMTS AT EACH STATION CONNECTEO ro FREE F<U) SOIL DISPLACEMENT,,*. - ABSOLUTE REFERENCE AXtS INTERACTION ELEMENT SCA^f Fig. 2.1 Pile model used in SPASM analysis. (after Finn and Gohl, 1990) Interaction between the near-field soil and the pile during shaking is represented using non-linear lateral springs k placed along the length of the pile. The spring stiflhess s 12 ks for a particular pile deflection level is defined from tangent slope of the non-linear p~y curves. A back bone p~y curve f(y) is specified in the positive and negative quadrants of the p~y space. Unloading p~y response from a peak loading point (p,^ , ymax) is simulated by following a path in the p~y space that mimics the initial backbone response specified using a function f(y- ymax ). Equivalent viscous dashpots Cr are also placed parallel with the near-field springs to simulate the radiation of P and S waves away from the pile (Fig. 2.1). Time varying free-field displacements, which are varied along the length of the pile, are applied to the free-field end of the spring-dashpot assemblage and represent the input seismic excitation applied to the pile. The input ground displacements could be computed using free-field response analysis, S H A K E (Schnabel et al., 1972). In dynamic analysis of the above system radiation damping coefficients are assumed to remain constant during shaking. Verification studies of SPASM's method of analysis were done by Gohl (1992) using data from centrifuge tests. They showed that the dynamic response of the structure is sensitive to the time history of free-field displacements. The S P A S M program may underpredict pile flexural response. A key difficulty in using S P A S M is the accurate determination of free-field input motions to be used along the embedded length of the pile. In realistic terms the so called coupled analysis in S P A S M is actually a semi-coupled method. The method only couples the superstructure with piles, but it does not couple piles with the surrounding soil directly. In addition, the program S H A K E , in predicting the free-field ground motion, adopts a equivalent-linear dynamic analysis. If the earthquake induced shear strains are large, say to the order of 1 percent, the accuracy of the equivalent linear solution fails to yield an acceptably good solution (Ishihara, 1982). For 13 the problem of larger shear strains in the soil deposits a step-by-step time integration procedure must be used. Kagawa and Kraft (1980a, 1980b) carried out a comprehensive study on the soilpile stiffness using a closed-form solution for linear soil-pile-superstructure systems and a finite element analysis. They studied the p~y response for homogenous and linear soil conditions as well as for heterogeneous and nonlinear soil conditions. The free-field soil was adjusted to accompany hyperbolic models (Konder, 1963 Hardin and Drnevich, 1972; Ramberg-Osgood model(1943), and Martin-Davidenkov model (1975)) with backbone ^MASS-1, mi 1 [STORY-1 —•MASS-2. m 2 ISTQRY-2 h 2 h 4 FLEXURAL y* RIGIDITY. E!j G$ • Snear .vocuius ol Soil P • Mass Oensiiy ol Soil » • Poissor. sfiaiiool Sal (1 • Oampmg Patio oi Soil El • Fienurai B^.any oi Pile r - RaOnjs o' P.'ls . M • HeigM o! Serf Laver p - Mas* Derm ol Pile s IH 0 2 0 FREE TO ROTATE Fig. 2.2 Analytical model for p~y evaluation and mathematical model for soil-pile springs (after Kagawa and Kraft, 1980) 14 type curves using Masing rule (1926). Their analysis procedure is illustrated in Fig 2.2. Based on this analysis they concluded that 70 percent of the pile displacement is concentrated in the soil mass within 2 pile radius distance. And the increase in shear strain due to soil-pile interaction is concentrated in this zone. Although it appears to be an effective solution technique it does not include pile separation, nonlinear material and geometric characteristics of pile and the dynamic P-A effects. The "state of the art" procedure for three dimensional dynamic soil-structure interaction was developed by Wu (1994). It includes a quasi 3D finite element method of analysis to determine the dynamic response of pile foundations in a coupled form subjected to horizontal loading. Both elastic response of piles and non-linear soil-pilesuperstructure were performed for single piles and pile groups. In the soil domain, movement in vertical direction and horizontal cross shaking were neglected and hence a quasi 3D situation was established. The free and near-field soil were modeled using brick elements with horizontal displacement as a degree of freedom in each of the nodes. The pile was modeled with beam elements with horizontal translation and its slope as the degrees of freedom in each node. The linear analysis was carried out in the frequency domain and the accuracy of the model was compared with that of Kaynia and Kausel (1982), Novak et al., (1977) and Fan et al., (1991). The computed impedances and kinematic interaction were found to be similar. It may be concluded that for a linear solution the plain strain models can provide reasonably adequate prediction of pile head behavior. Extensive time domain non-linear finite element analysis were also carried out by Wu to study the soil-pile-superstructure problem and they are incorporated in the program 15 PILE3D (Wu, 1994). It adopts a equivalent-linear approach in which the stiffness was represented by a complex shear modulus G* = G(l + /'. 2X). In it, the stiffness is associated with the real part of the complex shear modulus and hysteretic damping with the imaginary part, using equivalent viscous damping for soil (Seed and Idriss, 1970). As equivalent viscous damping is frequency dependent in Raleigh damping (Humar, 1990), frequencies were computed for the degrading stiffness of the system at selected steps and updated. As mentioned earlier, the equivalent linear approach fails to predict the response of the free-field soil at high strain levels in the order of 1%. Secondly, in long piles the geometric non-linearity becomes significant since the beam-column effect may not be discounted. Thirdly, in an earthquake loading situation, nonlinearity due to yielding in the pile material at specific points where the state of stress is high should be taken into account. Wu's (1994) work does not consider these issues. However, it is one of the best available analysis procedures available to obtain the three dimensional response of a coupled soil-pile-superstructure system. The preceding section discussed some of the formulation and analysis procedures available for single pile response. On this bases, the study of a coupled approach for dynamic soil-pile-superstructure interaction analysis, accommodating the free and nearfield non-linearities, pile-gapping, the material non-linearity and P-A effects for earthquake loading was considered worthwhile. A n attempt is made in this thesis to study all the above mentioned effects. A finite element formulation is presented in next Chapter. In Chapter 4 analysis for sample problems are performed. Discussions of the results are presented in Chapter 5. 16 CHAPTER 3 FINITE ELEMENT FORMULATION OF SOIL-PILE-SUPERSTRUCTURE INTERACTION 3.1.0 INTRODUCTION Piles subjected to earthquake loading are of particular interest in connection with founding structures on soft soils, drilling platforms and other pile-supported industrial and defense installations. Lateral dynamic loads induced from wind, wave and earthquakes are frequently the most critical factor in the design of such structures. Solutions in a general form of this problem could be applicable in a variety of cases on shore, including power poles, pile-supports for earthquakes and pile-supported structures which may be subjected to lateral blast forces. The problem of earthquake loading in a pile is closely related to the familiar problem of a beam vibrating on a elastic foundation; however, in one respect, it represents a more specialized case. All time varying external forces and moments applied to the soilpile system are introduced through boundary conditions existing at the top of the pile, while time varying pressures are applied along the beam. Solutions of dynamic soil-pilesuperstructure interaction problems require generalization of the beam vibrating on a elastic foundation theory to account for non-linear soil and pile behavior. Fig 1.1 shows the general nature of soil-pile-superstructure subjected to seismic excitations. If a long pile (for example, as used in an offshore environment) and its foundation experience large displacements, it may become necessary to consider geometric non-linearity in addition to material non-linearity in the problem. During an earthquake the soil-pile-superstructure system is subjected to loading, unloading and reloading which in turn induces significant non-linearities in the near-field and free-field soil and in the pile itself. Utilizing a continuum approach the geometry of the pile at any instant, including gaps, if any, could modeled by cushioning or iterating to attain that shape. A finite element analysis procedure and a computer program, P I L E U B C , to obtain the soil-pile-superstructure interaction was developed based on a coupled and unified system. The basic concept underlying the finite element method is that a physical domain can be modeled by subdividing it into a finite number of elements. Within each finite element, a set of functions are assumed to approximate the stresses or displacements in that region. The set of interpolation functions contain unknown parameters and are chosen to ensure continuity throughout the domain. Application of the principle of virtual work yields a system of algebraic equations for the parameters in the shape functions. This Chapter sets out to explain the procedures of finite element formulation of soil-pilesuperstructure system. The conventional 'stiffness' approach or the formulation based on an assumed displacement field is followed here. The virtual work method is then used to set up the mass and stiffness matrices and to determine nodal forces equivalent to body forces. 18 3.1.1 P R O B L E M F O R M U L A T I O N A N D FINITE E L E M E N T DISCRETIZATION In this new formulation a single pile is allowed to experience large deformations due to ground excitation and could behave as an elastic-plastic material. Superstructure effect could be incorporated by introducing a concentrated mass at the top of the pile in the global domain. To model the soil-pile interaction, a continuous non-linear foundation is defined using p~y curves of either A P I , Yan-Byrne or defined by data points. The timevarying free-field displacements, which vary along the length of pile buried inside the soil are obtained from a continuum based one dimensional finite formulation (Section 3.2). The attenuated ground displacement time-histories are obtained from a finite element program called Q U I V E R exclusively developed for this purpose. The theory and formulation behind this program considers the response associated with vertical propagation of shear waves, which in turn induces transverse vibration throughout a heterogeneous soil deposit. A one dimensional element of length A with two end nodes, as shown in Fig. 3.2, p is used in the finite element formulation of the soil-pile-superstructure interaction. With the element shape functions expressed in local coordinates £(-1 < i; < 1) and r\(-l < TJ < 1) the coordinates of any point inside the element can be expressed as «* = * 2(z-z') y c A ri = — ' d ) and and A dz=-d% 2 * dy = — dr\ * 2 ' (3.1.1.1) (3.1.1.2) where z is the z-coordinate of the center of the element length, d is the diameter of the e c pile in a circular section. The width of the cross-section, b(rj), may be constant or varying with TJ, as shown in Fig 3.3. The consistent element mass matrix is evaluated in closed 19 w w W , II w u Degrees of freedom w u Fig. 3.2 Finite Element Discretization for File 20 b(r,) = b form as will be seen later. The pile is represented by an assemblage of one dimensional finite elements that are connected through nodal points. A finite element of length A between nodes i and j is shown in Fig. 3.2. The degrees of freedom at each node of the element are the bending (lateral) deflection, w, its slope and curvature, the axial translation, u, and its first derivative. Thus there are five degrees of freedom per node and are arranged in a vector {a} as follows a = {w w' w"i, u ui, w , w' , w" , u , u } (3.1.1.3) T h lt h 2 2 2 2 2 The interpolation function that describes the variation of the unknown lateral displacements within the element in terms of nodal values is quintic, and it produces a cubic variation in translational lateral strain within the element. The interpolating function for the axial displacements u is cubic. The different pile elements and the associated layers are numbered from bottom to top. 3.1.2 M A T E R I A L S T R E S S - S T R A I N RELATIONSHIPS The pile is assumed to behave as an elastic-perfectly plastic body. B y defining the yield stress of the material and the initial modulus of the stress-strain curve, the stressstrain behavior is fully defined. Unloading and reloading is assumed to be parallel to the initial loading curve resulting in hysteresis loop as shown in Fig 3.4. 3.1.3 N E A R FIELD A N D PILE RELATIONSHIPS A n earthquake basically has two effects on a pile foundation. It excites the superstructure supported by the piles so that the inertia forces (lateral shear, overturning moment and axial load) are applied to the pile at the ground surface. And it excites the 22 Fig. 3.4 Elastic, perfectly plastic pile material 23 ground surrounding the pile so that additional soil "interaction" forces are applied to the piles along their embedded depth due to relative movements between the pile and the moving ground. Extensive studies, mentioned earlier, have shown that non-linear soil-pile interaction dramatically influences pile response. In particular it is indicative that a strain softening is predominant whose stiffness and hysteretic damping are load level dependent. This is important from a practical point of view, since the majority of dynamic seismic analysis carried out use equivalent linear characterizations of pile head stiffness and damping (referred to as pile head compliance functions). This is indeed practical for low amplitude vibrations representative of machine vibrations. The near-field soil in this formulation is modeled as continuous non-linear foundation using p~y curves as shown in Fig 3.5. One of the p~y curves generally used in analysis are those specified in the American Petroleum Institute (API) Code (1982) for design of fixed offshore platforms. These curves are specified for soft clay and sandy materials. The API p~y curve was initially developed by Cox et al., and Reese et al.,(1974) and modified by Murchison and O'Neill (1984). The p~y curves are expressed by a hyperbolic tangent function: p = r\A p tanh (3.1.3.1) u in which p is the ultimate soil resistance taken as the lesser value from the following: u P = (CiZ + C D)yz (3.1.3.2) p = &Dyz (3.1.3.3) u 2 u in which 24 Fig. 3.5 Soil-Pile Interaction CI, C2 and C3 are a function of soil friction angle, ^ y is the unit weight of soil 7j is a factor used to describe the pile shape effect z is the soil depth at consideration D is the pile diameter and rihi is the coefficient of subgrade reaction modulus and is a function of soil relative density. Yan and Byrne (1992) conducted extensive experimental testing to verify the validity of the API p~y curves based on hydraulic gradient similitude technique (HGS) for sandy soils. They concluded that the A P I curves based on the hyperbolic tangent tend to over-estimate the pile stiffness and proposed a parabolic p~y curve that better resembles experimental curves. Their average curve is expressed as: P (3.1.3.4) where D = diameter of the pile. Emax = Initial modulus at that particular depth. a is a function of soil density, and P has a value of about 0.5. The computer program P I L E U B C written for non-linear seismic soil-pile-superstructure interaction is capable of handling the API curve , the Yan-Byrne curve or any other user defined curve through data points. In addition the program allows for different type of curves at different depths. The initial modulus in a heterogeneous soil deposit, 26 E max =2G Gmax, max ( l + v), where Gmax is the shear modulus of the soil and v is the Poison's ratio. in turn, is a function of the effective pressure at the particular depth. 3.1.4 PILE S E P A R A T I O N (GAPPING) Soil response to cyclic pile head loading in saturated clay has been summarized by Bea (1980). For example, it has been shown that lateral resistance p of a soft clay to lateral pile motions y, during a sufficient number of cycles of loading, soil-pile separation (gapping) and soil remolding occurs. The confining stresses in the soil near the pile head are not sufficient to close the gaps that may develop between soil and the pile during the cyclic loading, resulting in dog-boned hysteresis loops of the p~y response. The development of gaps can be expected to lead to significant increases in pile displacement, shifting into a longer period of the coupled system and reduction or no contact in pile to soil interaction. Soil gapping in this continuum dynamic formulation is modeled by assuming that the soil is incapable of carrying any tensile forces and monitoring whether there is soil-pile contact or not. Numerical studies on static response were made by Vazhinkhoo (1994) on single pile, incorporated in the program CYCPLLE, subjected to cyclic loading. In comparing field tests it has been observed that this assumption gives reliable results. Further work (Vazhinkhoo, 1994) is underway to develop a mathematical model for cyclic p~y behavior. 27 3.1.5 KINEMATICAL R E L A T I O N S A N D S H A P E F U N C T I O N S For the lateral response of pile the external work comprises the work done by the inertial forces, the soil reaction due to relative movement of soil and the pile and the damping forces. The internal work done corresponds to internal stresses which depend on the geometry of the pile and its constitutive relations. The stresses in the pile are a function of strain as the pile undergoes large geometric deformations and behaves inelastically. Let the stresses a obey the constitutive equation (3.1.5.1) a = where s is the strain at any point in the pile. This strain is the sum of axial strain plus bending strain. Axial strain, en 1 fdw\2 du e o = -dx~ 2{-dx~ + (3.1.5.2) Bending strain, ej, Sb = dwb dx = -z Idx ) 2 (3.1.5.3) Therefore total strain, e (3.1.5.4) 28 Continuity conditions require a displacement field with continuous derivatives of one order less than the highest derivative appearing in the strain-displacement relation. Thus w, w', and u must be continuos. Specifying w, w' as degrees of freedom ensure continuity of displacements and slopes, while w" ensures continuity of bending moment in the elastic case. In general the inclusion of w" and u' as degrees of freedom permits the use of fewer elements in obtaining the solution. Since it is necessary to choose the displacement function such that w' be continuous at the nodes, w must be assumed at least a cubic polynomial. However, a complete quintic interpolation is used to approximate w within each element. A complete quintic polynomial requires 6 parameters to define a function. For u, a linear interpolation would be sufficient, but a cubic polynomial is used instead, requiring 4 parameters. The lateral displacements, their first and second derivatives , and the axial displacement and its first derivative at the two nodes provide sufficient parameters to fully describe the specified polynomial functions. Hence we could express w, w', w ", u and u as w = Niwt + N2 dw} + N N6wj + Ni + NAU, dw\ \dx) j + N 5 \dx j 2 du; .dx) j t + N Uj 9 + N10 du {dxjj (3-1.5.5) where subscripts /' and j refer to the first and second nodes respectively of each element. Alternatively, Equation 3.1.6.1 could be written as w = MI{O],w' = MJ {a},w" = MJ {a},u = N$ {a}andw' = N? {a} (3.1.5.6) In which 29 M (2,4) = (-12 + 604 +124 - 60i| )(2 /16A) M o © 2 3 2 M (U) = ( 8-15^ + 10^.3^5)/16 0 M (2,4) = (5-7cV6;; +10ci +c; -34 )/(A/32) M (3,4) = ( -4 +124 + 124 - 204 ) /16 M (3,4) = ( 1^-2^2+2^+^ ^5 )/( 2 M (4,4) = 0 2 3 5 4 0 0 A /64) 2 3 2 2 M (4,4) = 0 M (5,4) = 0 M (5,4) = 0 M (6,4) = ( -604 + 604 ) (4 / 16A ) 2 0 3 2 2 0 M (6,4) = ( 8 + 154 -104 + 34 ) / 16 M (7,4) = (12 + 604 -124 - 604 ) (2 /16A) M (7,4) = ( -5-7cl+6c;2 i0c|3^4.3^5 /32) M (8,4) = (-4 -124 + 1242 + 204 ) /16 M (8,4) - ( l+£=-2^2_2^3+^4 +^5 )/( 2/ 4) M (9,4) = 0 M (9,4) = 0 M (10,4) = 0 3 2 5 0 + )/(A 0 A 3 2 0 6 3 2 2 2 0 N © M (10£) = 0 0 0 N (i,4) = o 0 N (2,4) = 0 0 1^(1,4) = ( - 15+ 304 - 15i; ) (2 / 16A) N (3,4) = 0 M!(2,§) = (-7-12^+30ct2 ^3_ ^4y N (4,4) = (2-34 + 4 )/4 2 4 +4 15 16 Mj(3,4) = (-l-44+64 +44 -54 )(A/32) 2 3 4 0 3 0 N (5,4) = ( l - 4 - ^ 0 2 +^ )(A/8) 3 N (6,4) = 0 M (4,§) = 0 0 1 N (7,4) = 0 0 M!(5,iD = 0 1^(6,4) = ( 15- 304 +154 ) (2 / 16A) 2 N (8,4) = 0 0 4 N (9,4)=(2 + 34-4 )/4 3 Mi(7,§) = (-7-T-124+304 - 44 -15^ )/16 2 3 4 0 N (10,4) = (-I-4 + 42 +4 )(A/8) 3 0 M!(8,ij) = ( 1-44-64 + 4$ + 54 )(A/32) 2 M (9£) L 4 Ni© =0 Mi(10,4) = 0 Ni(l,4) = 0 Ni(2,4) = 0 N,(3,4) = 0 M (l,4) = ( 60S " 604 ) (4 / 16A ) 3 2 2 N!(4,4) = (-3 + 342) /(2A) N!(5,§) = ( - l - 2 § + 3§2 ) / 4 Ni(8£) = 0 Ni(6,§) = 0 Ni(9,§) = ( 3 - 3 § ) /(2A) Ni(7,§) = 0 N (10^H-l+2^+3^)/ 2 1 4 where £ is the Gaussian coordinate (-1<£<1) between nodes i and j . 3.1.6 V I R T U A L W O R K E Q U A T I O N S The principle of virtual work is used to derive the equations of motion for the soilpile-superstructure system. This is illustrated in Fig. 3.4. Let Sa be a virtual perturbation of the vector a. This results, by Equation 3.1.5.6 in displacements and strains within the element equal to 6w = MT {8a}, &V = MT {da}, 6V'= M \ {da} ,SU = NT {da} and dil = NT {&*} (3.1.6.1) respectively. The internal virtual work per unit volume done by the stresses is (3.1.6.2) In this formulation the interaction between soil and pile during shaking could be represented using nonlinear continuous foundation along the length of the pile. The nearfield stiffnesses for a particular deflection level are defined from the tangent slope of nonlinear p-y curves defined earlier. Soil gapping could be heeded by keeping track of the hysteresis loop. In other words while unloading the previous displacement could be considered as gap. The soil deposit could be heterogeneous stratum with modulus varying non linearly. Let v(x,t) is the ground displacement relative to the till, the displacement of which is v (t). Then the relative displacement of the pile to which the soil would react would be g 31 (3.1.6.3) s = w-v(x,t) The virtual work done by the soil reaction forces is J P\ ~( > tVjSnvdx L " w v x (3.1.6.4) Similarly, the work done by the inertial forces in the pile, assuming no rotatory inertia and that there is no vertical accelerations at the till is ( •• • • \ pA j I w + v g pw + u Sudx (3.1.6.5) L Since it is a two dimensional domain, equating the external work with the total internal work obtained by integrating over the volume of the element, V, we have (3.1.6.6) pAl(w + vg)dw +dudx + j p{w-v(x,t))owdx + jo(e)dedV-0 cA ' o v u This must hold good for any da, substituting Equations 6.1.6.2 through 6.1.6.4 in Equation 6.1.6.6 to work in terms of local coordinates and introducing vector notations form, we obtain the equation of motion pA^- )^Mo(Z)Ml(tya + N fe)No($'a)dtl +M 0 s t r u c ,urefy + Ad 4 A - 1 - — P^j Mo dt,V + Mstructure { } e g V s (3 16 7) in which M a t u r e is the mass matrix for the superstructure and {e} is unit vector . 32 This equation of motion is integrated using the Newmark-P average acceleration method (Newmark, 1959) .Thus we write a function W(a) : A f - l - > / A ,+ ; M 4 A JHe, + i)JM© -(rt^MQ + M © ) +M©M r 1 ( ^ ) a ^ r ^ + | j/^) - 1 Vgr+i "I" A / structure {^} Vgj+i ^ -1 pA^ }(MO{Z>)MI(£>)+No(Z,)N% ($) +)d^ +structure M i (3.1.6.8) = 0 Solution of the problem then requires finding so that W(a = 0. However, it is not possible to achieve an exact solution in this nonlinear problem, and an iterative procedure is required. It is necessary to ensure that the total error W(a) is reduced to a given tolerance. The Newton-Raphson method is a widely used iterative technique to solve non-linear equations. It uses a truncated Taylor series expansion of the function *¥(a + Aa) = Y ( a ) + ^^-Aa da Aa = fflF(a) =0 (3.1.6.9) (3.1.6.10) da As *F(a+Aa) = 0, Aa = a -a, = M ~d¥(a)~ da (3.1.6.11) 33 or in matrix language ay(q) da ai+i - at (3.1.6.12) where {¥{} is the vector at previous iteration and is the solution vector at the next iteration. It should be realized that these iteration are in between time steps t and t i. The t matrix i+ dW(a)~ Equation 3.1.6.9 is da structure +^ Ha(8){Mft)Mi ft)}KT0*ftl + f r }^A/oft)A/Jft)* (3.1.6.13) The displacement solution vector will take the following form at the end of every iteration within a time step. w=w-[vvr{-^r , i (3.1.6.14) The displacement vector obtained through the solution procedure as in Equation 3.1.6.14 is checked for convergence of the error, y/, to a tolerance. The tangent stiffness matrix is evaluated numerically using Gauss quadrature and the consistent mass matrix is evaluated in closed form. 34 3.1.7 C O N S I S T E N T M A S S MATRIX F O R PILE E L E M E N T The consistent mass matrix of the pile element was evaluated analytically and is given as follows: 181pA 462 311p A 4620 52 A 3465 3 P 2 281p A 55440 23 A 18480 PA 9240 4 P 3 0 0 25pA 231 151pA 4620 181p A 55440 151p A 4620 19p A 1980 13 A 13860 0 0 0 0 0 0 181pA 462 311p A 4620 52pA 3465 2 0 0 0 0 13pA 35 Hp A 210 PA 105 3 5 2 2 3 181p A 55440 13 A 13860 PA 11088 3 4 P 4 0 0 0 0 0 0 9pA 70 13pA 420 13pA 420 PA 140 0 0 0 0 0 0 13pA 35 Hp A 210 PA 105 5 P 3 2 2 3 2 281p A 55440 23 A 18480 PA 9240 3 3 4 P 5 3 Symmetrical about the leading diagonal (3.1.7.1) 3.1.8 N U M E R I C A L INTEGRATION The near-field continuous foundation, pile material and geometric non-linearity leads to changing stiffness in the system when excited. Hence it is convenient to numerically integrate Equation 3.1.6.14 during every iteration in the Newton-Raphson procedure. The well proved Gauss quadrature is used. In general, an integral I / = / f(£)d£ -l J = SHifU) l then (3.1.8.1) (3.1.8.2) 35 2 where are fixed abscise and H corresponding weights for the n chosen sampling points. t If a polynomial expression is to be integrated it is easy to see that for n sampling points, there are 2n unknowns (£, and H,) and hence a polynomial of degree 2n-l can be constructed and exactly integrated. The error thus is of the order of 0(A ). 2n As indicated earlier, n Gauss points integrate a polynomial of order (2n-l) exactly. In the gradient of Equation 3.1.6.13 in the Newton-Raphson procedure, there is a polynomial of order 8. Hence a 5-point Gauss Quadrature would suffice to integrate exactly. However it is noteworthy to mention that in the Program P I L E U B C there is a flexibility in choosing the number of Gauss points from 1 to 9. It is understood that as numerical integration using Gaussian scheme is employed the pile stresses and strains could only be evaluated at the sampling points. A suitable interpolation should be adopted to evaluate them at other locations. 3.1.9 DISCRETIZATION E R R O R A N D C O N V E R G E N C E The assumed shape functions limit the infinite degrees of freedom of the continuous soil-pile-superstructure system, and the true minimum of the energy may never be reached, irrespective of the fineness of the subdivision (Zienkiewicz, 1977). However it is possible to minimize the total error to a tolerated value. A n Euclidean norm criterion has been used to check the convergence at the end of every iteration in the Newton-Raphson procedure. If \ar\ is the solution vector at the end 36 of r" iteration and is the solution at the end of (r+7/ 1 Wr+1 A iteration, then il/2 Aa= \ar+i~a \ r represents the difference between the two displacement vectors. In which, (3.1.9.1) a, NDOP \a i\ (3.1.9.2) 1=1 r+ where a (i) and a +i(i) are components of a and r r r a i r+ respectively. Then Aa could be written as NDOF -r r Aa = £ V(a-i(0-ar(0) (3.1.9.3) The convergence criterion for this norm is given as Aa < specified tolerance (3.1.9.4) Secondly, it is verified that the total error in the problem {¥*} is zero. The mathematical form are w-w* — -^r—<TOL\ 1 a and <TOL2 (3.1.9.5) (3.1.9.6) The convergence tolerance for Equation 3.1.6.14 needs special mention. As the unit for {T} is that of force, a tolerance IN might be appropriate. If the solution obtained has not converged then the whole procedure is repeated by updating the solution vector in Equation 3.1.6.8. If converged, then the corresponding 37 acceleration and velocity vectors are calculated by Newmark-P method. The entire procedure is repeated for all the time steps. 3.1.10 P I L E U B C - T H E C O M P U T E R P R O G R A M A computer program, P I L E U B C , based on the aforementioned nonlinear finite element procedure was developed. Its salient features include changing Emax for every single element as a function of effective or total stress, multi-materialing of the pile, multilayering in the soil deposit, energy dissipation (hysteretic damping) by following through the elastic-plastic pile material behavior, non-linear P~y curves for near field in every single Gauss sampling point throughout the depth for the entire time history, evaluation of cyclic stress ratio and flexibility in giving the input motion at anywhere in the soil profile. Results of a few runs are enclosed in the next chapter. 3.2 FINITE E L E M E N T F O R M U L A T I O N O F F R E E - F I E L D SOIL D E P O S I T 3.2.0 INTRODUCTION In the coupled problem of seismic soil-pile interaction, the pile displacement are sensitive to spatial variation of the free-field soil. Various idealized models and analytical techniques available are, in general, equivalent-linear, and their validity is questionable when the earthquake-induced shear strains are high. Whatever procedure is followed, it is necessary to determine the appropriate shear stress-strain and energy absorbing properties in the deposit. During an earthquake a soil deposit is subjected to an irregular loading pattern due to a vertically propagating wave which consists of intervals of loading, 38 unloading, and reloading. It has different behavior characteristics in each of the loading phases. The finite element method is a versatile tool to handle the complexities of the soil problem. 3.2.1 FINITE E L E M E N T DISCRETIZATION O F SOIL D E P O S I T The soil deposit is modeled as an assemblage of one dimensional finite elements interconnected at a finite number of nodal points. Details of the element arrangement and the finite element discretization is shown in Fig. 3.6. The domain could be subdivided into elements of length A. The degrees of freedom at each node of the element are the horizontal translation and its first derivative. The interpolation function that describes the variation of the unknown displacement within the element in terms of nodal displacements is cubic, and it produces a quadratic variation in strain within the element. The different soil layers are numbered from top to bottom. 3.2.2 SOIL S H E A R S T R E S S - S T R A I N RELATIONSHIPS Vertically propagating seismic waves at a site produces transverse movement throughout the depth of the soil stratum. This movement may induce large shear strains in the soil. The shear stress-strain modeling of soils for drained and undrained conditions is usually done by the hyperbolic model. In addition, seismic loading imposes irregular loading pulses which consist of loading, unloading and reloading. Modeling of this soil behavior in cyclic loading is usually made by first specifying the stress-strain relation in the virgin or initial loading, using a backbone or skeleton curve. Subsequent unloading and reloading characteristics are handled using the Masing criterion (1926) as illustrated in 39 G L node 1 e 6 M node 2 * =1 Finite Element • Fig. 3.6. Soil Deposit 40 Fig. 3.7. In the program QUIVER, the behavior of soil is treated to be non-linear and hysteretic, exhibiting Masing behavior during loading and unloading. The tangent shear modulus G for any one particular level of strain could be t represented by the three parameter hyperbolic model (Duncan and Barkeley) given by 5T f _ dy ~ ^ m a x 1 T (3.2.2.1) R f V - in which, G = the maximum shear modulus that occurs at zero strain at the beginning of m a x cycling, r = the shear stress, Tf = the failure shear stress, Rf = the failure ratio T/TUU in which T„it is the ultimate strength. Rf may also be considered a factor that defines the strain,Yf, at which the failure stress occurs, Yf =VG m a x (l-R ) f The parameter Rf is used to modify the hyperbola to fit the laboratory data. Rf = 0 specifies a linear elastic plastic material with yf = Tf/Gmax as shown in Fig. 3.8. Rf = 1 specifies a strain hardening material with yf equal to infinity. For most sands the , Rf lies between 0.5 and 0.9 (Byrne and Mclntyre). The modified shear stress-strain relationship under cyclic loading could be modeled as a hyperbola and the tangent stiffness at any stress level is expressed as f * Gt - ^ - Gmax I * \ 2 i-^rRf (3.2.2.2) V J the shear modulus immediately upon unloading, 41 Fig. 3.7 Loading-Unloading Curve for Shear Stress-Strain Modeling 42 X Fig. 3.8 Modified Hyperbolic Shear Stress-Strain Formulation (after Byrne and Mclntyre, 1994) 43 and Tf = (TA + Tf) TA = the shear stress at the reversal point as shown in Fig. 3.7. Rearranging and integrating Equation 3.2.2.2. gives equations 3.2.2.3 through 3.2.2.5, hence the generalized expression for rfor all types of loading could be obtained. 1 f • > \ Tf J aty = 0; z= * 0 ^G * Rf (3.2.2.3) Tf y (3.2.2.4) Umax I If 1+ GmaxY or, Gmax Y (3.2.2.5) 1 + ^-GmaxY 3.2.3 C O M P U T A T I O N O F H Y P E R B O L I C M O D E L P A R A M E T E R S The hyperbolic parameters in Equation 3.2.2.5 are influenced by a number of factors. In order to obtain a reasonable calculation it is necessary that the most important factors are accommodated appropriately. The maximum shear modulus G max obtained generally from in-situ measurements, , could be laboratory experiments or empirical equations (Byrne, 1994). There is no one procedure that is accepted by different 44 researchers to evaluate the maximum shear modulus. In Q U I V E R Gmax for sand is calculated using the following equation Gmax (3.2.3.1) kg Pa in which kg varies from 500 to 2000 for sands and is a function of normalized SPT blow count, a m = mean normal effective stress, P = atmospheric pressure, a m = constant depending on the type of soil. The equation for k suggested by Byrne (1994) takes the form g where, (Ni) = is the normalized SPT blow count. 60 The mean normal stress <r is calculated from the following expression, m G m = y"0- + °) 2fc (3.2.3.3) in which, CT = the effective normal stress, V £ = 0.5 0 The failure shear strength, Tf for soils is dependent on the stress system by means of which the soil element is brought to failure. Hardin and Drenvich (1972) suggested that the value of Tf calculated using Mohr-Coulomb failure envelope (defined by the static strength parameters such as c , effective cohesion, and <f> , internal angle of friction) is 45 adequate for dynamic loading. Modified Coulomb's model for shear strength of a soil might be expressed in the form (Scott, 1980) Xf = c' + a'v tan <J>' (3.2.3.4) where Tf cr is the shear stress in the soil at failure, v is the effective normal stress. c, <f> are approximately constant parameters expressing the shear strength in terms of effective stress. For sandy soils c is close to zero. 3.2.4 P R O B L E M F O R M U L A T I O N A one dimensional element of length A with two end nodes, as shown in Fig. 3.6 is used in the finite element formulation of the free-field. The displacement u is relative to the till, which itself displaces an amount Vg, with the acceleration v g . With the element shape functions expressed in local coordinate c; (-1 ^ \ ^ 1) the axial coordinate of any point inside the element can be expressed as e 2(z-z*) , ^ where z e c , A , * = 7« ' (3.2.4.1) is the z-coordinate of the center of the element length. Numerical integration has to be performed inside every element and summed over to obtain element stiffness matrix. The consistent element mass matrix is evaluated in closed form as will be seen 46 later. The displacement vector consists of nodal degrees of freedom u and its first derivative u. They are of the following form a = {u T h Hi, U2, (3.2.4.2) u} 2 This sequence gives the flexibility of having different G max for each element in a layer and more than one layer in the problem domain. It permits strain discontinuity between elements, a parabolic shear strain within an element and satisfaction of continuity in the shear stresses. The grouping of nodal degrees of freedom is worth mentioning. The nodes at ground surface and at the till have two degrees of freedom, and the other nodes have three degrees of freedom, u and w/, in which /' and j correspond to two adjacent h elements. Thus, except at the ends, there are three degrees of freedom per node. 3.2.5 S H A P E F U N C T I O N S For the one dimensional soil element the strain-displacement relationship contains first derivative in the lateral translations. Hence it is necessary to choose the displacement function such that u is continuous at the nodes. This can best achieved by adopting a linear displacement for u. However, complete cubic interpolations are used to approximate u within each element. A complete cubic polynomial requires 4 parameters to define a function. The lateral displacements and their first derivatives at the two nodes provide sufficient parameters to fully describe a cubic polynomial function. Hence we could express u and u u = N m + N2 (du\ + N3UJ Vdx); + N4 (du\ (3.2.5.1) 47 where subscripts / and j refer to the first and second nodes respectively of each element and Ni=(2-3£ +£ ) / 4 (3.2.5.2) 3 N2= ( 1 - -£ + £ ) (A / 8) 2 (3.2.5.3) 3 N =( 2 + 3 £ - £ 3 ) / 4 (3.2.5.4) 3 N = ( -1 - £ + £ + £3) (A / 8) 2 4 (3.2.5.5) Alternatively, u and u' could be written as u = Ml{a] u=Ml{a) (3.2.5.6) where Mo = { N i N N N } and M i = {N,' N ' N ' N ' } T T 2 3 4 2 3 4 (3.2.5.7) 3.2.6 V I R T U A L W O R K E Q U A T I O N S F O R F R E E - F I E L D Principle of virtual work is used in formulating the element-stiffness and mass matrices. This is illustrated in Fig. 3.9. Let such a virtual displacement be Set at the nodes. This results, by Equation 3.2.5.6 in virtual displacements and strains within the element equal to Su = M {5a} T 0 oy = M {Sa} T l (3.2.6.1) respectively. The internal work per unit volume done by the stresses t(y) is T(Y)5Y (3.2.6.2) Equating the external work done by the inertial forces with the total internal work obtained by integrating over the volume of the element, one obtains Jx(Y)5y& = - J p | i + v j8«ak s (3.2.6.3) 48 dx 0 Work Equations A' iiiiiiiiijjjjjjjijii &W =1p(u+" )dudz e 0 £ Till V g Fig. 3.9 Virtual work for soil deposit where v is the acceleration at the till. This must hold good for any du, substituting g Equation 6.2.6.1 in 6.2.6.3 to work in terms of global system, we obtain the equation of motion ^ J P Mo Ml a%\a\ + ^ j x(y) M i a\ = - ^ ) p Mo v s (3.2.6.4) 3.2.7 CONSISTENT MASS MATRIX FOR SOIL ELEMENT The consistent mass matrix for the soil element evaluated analytically is of the form, 13pA 35 Hp A 210 13pA 420 9pA 2 2 70 11/0 A 210 PA 105 PA 140 13/> A 420 13pA 420 PA 140 PA 105 Hp A 210 2 2 3 3 3 70 13pA 420 UpA 210 13pA 2 3 2 9pA 2 2 35 3.2.8 TIME AND NUMERICAL INTEGRATION Newmark-P integration is adopted for solving the time dependent soil vibration problem. Equation 3.2.6.4 is rewritten below Y W o M o ^ { a , i } + | }T(Y,. )MI^ = ~ + +1 Wo^jv,'J (3.2.8.1) After approximating with Newmark-P procedure (average acceleration method), 50 A 1 — JpMoMo^ 2 -i ^ JpMoM^ (3.2.8.2) > -i where At is the time integration step. Let the right hand side of the above equation be {R}. Thus, we must find {a 4 JpMoMo^ % , 2 { « « } + f 2-i =O A Using the Newton-Raphson approach, a i+1 such that (3.2.8.3) can be found starting from an initial a . Then i+1 at+i - aj + (3.2.8.4) where ^ is the vector Wat a* and the matrix i+l 8W _ i / 4 JpMoMj^ 5a ~ / a A r 2 -i 2 A \ dz 2-idy T „ (3.2.8.5) As indicated in section 3.1.10, n Gauss points integrate a polynomial of order (2n-l) exactly. In the soil stiffness expression Equation 3.2.8.5 in the Newton-Raphson procedure, there is a polynomial of order 4. Hence a 3-point Gauss quadrature would suffice to integrate exactly. However it is noteworthy to mention that in the program Q U I V E R there is a flexibility in choosing the number of Gauss points from 1 to 9. It is understood that as numerical integration using Gaussian scheme is employed the shear stresses and strains could only be evaluated at the sampling points. A suitable interpolation should be adopted to evaluate them at other locations. 51 3.2.9 DISCRETIZATION E R R O R A N D C O N V E R G E N C E An Euclidean norm criterion has been used to check the convergence at the end of every iteration in the Newton-Raphson procedure. If of/* iteration and is the solution vector at the end is the solution at the end of (r+l) iteration, then Aa= th \Clr+l\ OTr+l ~ Or] represents the difference between the two displacement vectors. In which, NDOP (3.2.9.1) a, NDOP (3.2.9.2) where a (i) and a i(i) are components of a, and r r+ a +i r respectively. Then Aa could be written as NDOP ~z r. Aa = 2 yj(ar (J)-ar(i)) (3.2.9.3) +l The convergence criterion for this norm is given as Aa < specified tolerance Secondly, it is verified that the total error in the problem (3.2.9.4) } is zero. The mathematical form are H-H1 < TOL1 and (3.2.9.5) M" TOL2 (3.2.9.6) 52 The convergence tolerance for Equation 3.2.9.6 needs special mention. As the unit for {*P} is that of force, a tolerance of say IN ox so might be appropriate. If the solution obtained has not converged then the whole procedure is repeated by updating the solution vector in Equation 3.2.8.2. If converged, then the corresponding acceleration and velocity vectors are calculated by Newmark-P method. The entire procedure is repeated for all the time steps. 3.2.10 QUIVER - T H E C O M P U T E R P R O G R A M A computer program, QUIVER, based on the aforementioned nonlinear finite element procedure was developed. Its salient features include changing G m a x and T U for U every single element as a function of effective or total stress, multilayering in the soil deposit, energy dissipation (hysteretic damping) by following through the Duncan and Barkaley shear stress-strain model in every single Gauss sampling point throughout the depth for the entire time history, evaluation of cyclic stress ratio and flexibility in giving the input motion at anywhere in the soil profile. Results of sample runs are enclosed in the next Chapter. 53 CHAPTER 4 VERIFICATION AND NUMERICAL RESULTS 4.1 INTRODUCTION Numerical results obtained through the proposed formulation described in Chapter 3, and incorporated in the finite element programs Q U I V E R and P I L E U B C are now reported. As in any other analytical procedure the validity of the formulation was verified against linearized closed form solutions to check the accuracy of predictions. The results obtained from the programs are presented in two stages. First, the cases for the finite element analysis of the free-field soil deposit are presented. It includes a comparative analysis with a program S H A K E (Scnabel et. al., 1972). Secondly, the analysis results from the program P I L E U B C are shown. It includes comparative solutions using assumed linear systems and two example problems are dealt. The results comprise of the pile's relative displacement time history at specific nodes, relative displacement profile and bending moment profile at selected time intervals. It is understood that a variety of data could also be obtained from the program. For example, the time varying bending moment, shear force, displacement, velocity, acceleration, bending stress, axial and bending strain and gap history. Explanation of input and sample input data file are attached in Appendix II. 54 4.2 VERIFICATION O F F R E E - F I E L D M O V E M E N T A G A I N S T D E V E L O P E D C L O S E D FORM SOLUTION In order to verify the proposed finite element formulation for free-field soil as illustrated in Section 3.2, a closed form solution for a linear uniform stress field system was developed. If for instance the soil deposit is subjected to a harmonic base excitation and it is idealized as a uniform stress field as described in Appendix A , it is possible to obtain a time dependent exact closed form solution. In this situation, whatever the element subdivision, the finite element solution will coincide with the exact one (Zienkiewicz, 1977). This is an obvious corollary of the formulation and it is useful as a first check of the computer program QUIVER. This solution involves both the transient response and the steady state response of the dynamic soil problem. With this as intention, utilizing Laplace transformation a closed form expression is developed for harmonic base excitation of an elastic soil deposit. The soil deposit is idealized as one layer with constant Gmax throughout the depth with no damping in the system. The shear stress-strain behavior is assumed to be linear in this boundary value problem. It is important that all the assumptions inherent in the closed form are precisely represented in the finite element program. The soil deposit used for the comparison between the closed form solution and the finite element approach (QUIVER) was : the depth of the deposit (H) was 30m, the maximum dynamic shear modulus (Gmax) was 275000 kN/m , the failure shear stress was 2 375 kN/m , and the simulated sinusoidal base acceleration at the till was 0.5g Sin(cot), 2 where, g is the acceleration due to gravity and 0) is the frequency of excitation. 55 To conform with the assumptions used in the closed form solution, the maximum dynamic shear modulus and the failure shear stress were assumed to be constant throughout the depth. In addition the free-field was assumed to behave as a undamped linear system. For comparison, the time history of relative displacements and absolute accelerations on the surface computed from closed form and QTJTVER are plotted in Fig. 4.1. The time history for relative displacements throughout the profile are plotted in Fig. 4.2. It can be seen from both the above mentioned figures that the predicted accelerations and displacements match exactly with the closed form. 4.3 R E S U L T S F O R V A R Y I N G G m a x A N D x, IN A H E T E R O G E N E O U S SOIL DEPOSIT The maximum dynamic shear modulus is a function of effective pressure, cr , and v normalized SPT blow count, (Nj)6o . And the failure shear is a function of effective pressure, effective friction angle <> / ,and the cohesion intercept, c. These properties vary with depth and type of soil layer. It may not be possible to develop a closed form solution for this continuous heterogeneous deposit. Results of such a nonlinear finite element analysis for a simulated sinusoidal base acceleration of 0.2g Sin(eot) is enclosed in Fig 4.3. The pre-earthquake soil properties of the deposit used in this example problem are shown in Table 4.1. 56 0.84 g Ground Level 0.831 g 0.769 G =275000kNfeqm 0.738 g T = 350kN/sq.n| uH 0.7023 g •y = 20 kN/m' 0.666 g H = 30m 0.626 0.587 g 0.546 g 10 10.00 0.5 g Till Soil Deposit —• 11 0.00 -10.00 FE Representation Base Acceleration = 0.5 g sin (wt) w= 6.28 tad/see Fig 4.2 Absolute acceleration along the profile using 10 elements for Uniform Stress Field Peak acceleration = 0.38 g 5.00 0 m. 0.00 -5.00 4.00 1 Pea k acceleration = 0.36 g 3.0 m . Peak acceleration = 0.32 g 4.00 12.0 m_ SAND Y=20kN/rr? 3.00 Peak acceleration = 0.29 g 1 0.00 18.0 m- -3.00 Simulated acceleration Induced at tilL 30.0 m. Till STANDARD LAYER 0.00 1.00 2.00 3.00 Time (Sec) 4.00 5.00 Fig. 4.3. Acceleration Time History From QUIVER Simulated Sinusoidal Excitation 59 Table 4.1 Soil Properties Used In The Example Problem 1 Soil Depth (m) (NOso kg m <> l Rf y (kN/m ) Crust 3.00* 10 440(Ni)6o 0.4 33 0.9 20 Sand 30.00 10 440(N ) o 0.4 27 0.7 20 G m a x = kgPa(o' /Pa) x 6 3 m ir Time history of absolute accelerations of the same soil deposit of Table 4.1, for an earthquake excitation (San Fernando, 1971) are shown in Fig. 4.4. A comparative cyclic stress ratio plot from Q U I V E R and S H A K E is shown in Fig. 4.5. Cyclic stress ratio is defined (Seed and Harder, 1992) by the following relation, CSR = 0.65Ov (4.1) In which Tmax = Maximum driving shear stress ' = Effective overburden pressure. CT The results obtained from the finite element analysis are in good agreement with S H A K E except near the water table. It could be observed from Equation 3.2.3.1 that the shear modulus varies parabolically with depth and is a function of effective pressure. If a water table is introduced, the rate of change of Gmax changes rapidly near the water table. As a result the shear strain in the soil at the location of water table would be more, resulting in swift changes in cyclic stress ratio profile. However, the accelerations predicted by 60 Peak acceleration = 0.52 g 5.00 0 m. 5.00 1 Peak acceleration = 0.46g 0.00 3.0 m . -5.00 Peak acceleration = 0.391 g 4.00 12.0 m_ S A N D 7=20kN/irr 18.0 m_ Peak acceleration = 0.30 g 3.00 — i 0.00 — -3.00 — 1 Earthquake acceleration induced at till (San Fernando, 1971) 30.0 m . Till STANDARD LAYER T 0 5 10 15 T 20 T 25 30 35 40 Time (Sec) Fig.4.4 Acceleration Time History From QUIVER for Earthquake Excitation 61 Cyclic Stress Ratio 0.10 0.20 0.30 0.40 0.50 Fig. 4.5. Cyclic Stress Ratio V s Depth between Q U I V E R and S H A K E S H A K E at relatively deeper space in the soil deposit is higher than the finite element solution. Fig. 4.6 shows the variation of shear stress-strain at one of the Gauss sampling points in the soil deposit obtained from Q U I V E R for the above mentioned analysis. It could be observed from Fig. 4.6 that irregular loading the soil deposit and hence the response could be captured appropriately. 4.4 R E S U L T S F R O M A N A N A L Y S I S F O R C A N L E X SITE, A T RICHMOND Another example was chosen at a site in Richmond, British Columbia, Canada from field data available through seismic cone penetration (SCPTU) tests performed by the In Situ Testing Group, The University of British Columbia. The analysis aimed at the study the free-field response subjected to a base acceleration of the San Fernando, 1971 earthquake record. The site reference indicated in the project report, C A N L E X , P H A S E II, (Campanella, 1995) was KD9033. This site was chosen for the reason that interpreted data on soil properties was made available upto the till at 48.83m. Fig 4.7 shows the variation of normalized SPT blow count (Ni) with depth for this site. The pre-earthquake 60 soil properties used in the analysis are tabulated in Table 4.2. The input data file for the program is included in Appendix n. The absolute acceleration time history obtained from the finite element analysis and from S H A K E analysis is shown in Fig. 4.8 and 4.9 respectively. The base acceleration for this site was the acceleration record from San Fernando earthquake record, 1971. It is observed that S H A K E analysis predicts relatively a higher acceleration at one layer above the till. 63 30.00 , Fig. 4.6 Typical Shear Stress-Strain Plot from the Finite Element Program Fig. 4.7 SPT Data for Test No. KD9303 No. 4 Road at River Road, Richmond, BC, Canada 65 5.00 -0 m C L A Y E Y SILT 0.00 -5.00 1.25 m- y =18.86 kN/m 3 No. of Finite Elements = 10 1 Peak acceleration = 0.34 g 3.25 m- SA N DY SILT Y =18.86 lcN/m 3 No. of Finite Elements = 10 6.25 m- FINE SAND (NiXr -4.00 —" 24 Y =19.65 kN/m 3 No. of Finite Elements = 20 23 m_ 3.00 Peak acceleration = 0.26 g 0.00 C L A Y -3.00 Y =18.86 kN/m 3 No. of Finite Elements = 20 arthquake acceleration induced at till (San Fernando, 1971) Peak acceleration = 0.192 g 48.2 m Till " ttjf*> CANLEX, RICHMOND, KIDD 2 SITE Site data from, In-Situ Testing Group, UBC, Canada -2.00 T T 10 T 15 20 T 25 Time (Sec) Fig. 4.8. Acceleration Time History F r o m Q U I V E R T 30 35 40 Peak acceleration = 0.46 g 5.00 0 m. 0.00 C L A Y E Y SILT -5.00 — 1 1.25 m- (N,X.= 6 y =18.86 kN/m 3 n Peak acceleration = 0.46 g 3.25 m_ S A N D Y SILT -5.00 (N,) „= 11 6 y =18.86 kN/m 3 Peak acceleration = 0.5 g 6.25 m_ FINE SAND y =19.65 kN/m 3 Peak acceleration 0.S4 g B3 6.00 23 m_ — i 0.00 C L A Y (NiX.= -6.00 — 1 6 Y =18.86 kN/m 3 Earthquake acceleration Induced at till (San Fernando, 1971) 48.2 m. Till CANLEX, KIDD 2 SITE Ate data from, In-Situ Testing Group, UBC, Canada 10 15 20 25 Time (Sec) Fig. 4.9. Acceleration Time History From SHAKE 30 35 40 Table 4.2 Soil Properties for Nonlinear Analysis of C A N L E X SITE Soil Depth c (m) (kN/m ) Clayey Silt 1.25* 70.0 0 Clayey Silt 3.25 55.0 Sandy Silt 6.25 Fine Sand Clay (Nl>60 No. of OCR Y Finite (kN/m ) (ratio) Elements 6 18.86 100.0 5 0 6 18.86 35.0 5 0.0 41 11 18.86 0.0 10 23.0 0.0 41 24 19.65 0.0 20 48.83 65.0 0 6 18.86 1.0 20 2 3 * Water table ocation It could be observed from Fig. 4.4 and Fig. 4.8 that the time varying relative displacement obtained by the finite element approach continues to vibrate even after the base acceleration has attained very small amplitudes. The reason is that as only hysteretic damping is considered in the soil stress-strain model, on account of minuscule vibrations, the strain levels will be in the linear range and hence energy dissipation is infinitesimal. As a result, the undamped free vibration response is observed. 4.5 R E S U L T S F R O M FINITE E L E M E N T A N A L Y S I S O F SOIL-PILES U P E R S T R U C T U R E INTERACTION The results for seismic soil-pile-superstructure is presented in this section. Comparison is made with a linear system. Later results from two examples problems are 68 illustrated. A similar formulation for soil-pile interaction for static cyclic loading was formulated and incorporated in the finite element program C Y C P I L E at the University of British Columbia (Vazhinkhoo, 1995). A comparison between analysis and experimental for a case study is shown in Fig. 4.10. The analysis was done (Vazhinkhoo, 1995) for a laboratory test on an instrumented aluminum pile embedded in a dense sand subjected to monotonic pile head loading. Stresses in the soil were increased to field stress levels using Hydraulic Gradient Similitude technique (Yan and Byrne, 1992). Comparison was made between C Y C P I L E predictions and test data. The predicted results were observed to be in excellent agreement with experimental data. As C Y C P I L E and P I L E U B C are based on the same theory, the former only considering static cyclic loads, the agreement shown in Fig. 4.10 can also be seen as a verification of P I L E U B C for static loads. In order to check the dynamic component of PILEUBC, the vibration of a simple cantilever was studied. If the pile is assumed to be a linear cantilever beam of length, /, with lumped mass, m, vibrating at the free end for a harmonic sinusoidal base excitation with no viscous damping (Fig. 4.11), closed form solution is available. The time varying relative displacement at the free end is u(x, t) = UQ COS at +—sin cor —^—r (sin Clt - 0 sin cor) co A: 1 — P ^ in which u(x,t) = time varying relative displacement u = initial displacement 0 Vo = initial velocity m = lumped mass 69 X J X LATPILE — CYCPILE prediction - load input CYCPILE prediction - displacement input _ g 40 pradiction —£— CD CO (0 O 20 H -#—» (0 0 I I I I I I I I I J I I I I I I I II 0.0 0.5 | I I I II I I I I j I I I I I I I I I | 1.0 1.5 2.0 Lateral Pile Deflection at Loading Point (mm) Pile Bending Moment - N-mm -3000 -100 —f -2000 -1000 0 iiitiniiliiiiiiniliiiiiinilii 1000 il o-q E E 100 H 200 H sz +—> Q. Q 300 -=j H Fig. 4.10 Comparison between Experimental Results and CYCPILE Predictions 70 k = stiffness of the system, k = 3EI L 3 I = moment of inertia a> = natural frequency of the system, co = Q= frequency of excitation /?= frequency ratio, P = — co The cantilever beam used in a comparative analysis between the closed form and finite element approach (PILEUBC) had a length (L) of 12m, the external and internal diameter of 2.0 and 1.8 m respectively. Young's modulus ( 2 w ) was 210000 MPa, the yield stress was 700 Mpa, and the simulated sinusoidal base acceleration at the till was 0.2g Sin(fit). This system will have a natural frequency of 13.56 rad/sec. The frequency of excitation, /?, was chosen to be was 11.51 rad/sec. A flag was introduced in the finite element program such that there was no near-field soil and the yield stress of the pile material was set to be a high number such that there was no hysteretic damping in pile. P-A effects were neglected. Results of this comparative analysis are shown in Fig. 4.11. The finite element solution coincides with the closed form solution. The proposed non-linear formulation was applied to predict the pilesuperstructure behavior for an example problem shown in Fig. 4.12. The pile was a steel pipe, 32.4 m long and the external and internal diameter were 0.9 m and 0.7 m respectively. The pile was assumed to be embedded upto the till in sandy soil for a depth of 30 m. Weight of the superstructure was assumed to be 500 kN. The average normalized 72 standard penetration value, ( N i ) , for sandy soil was assumed to be 10. The following 60 values were adopted for the soil: (a) Bulk unit weight = 20 kN/m 3 (b) Angle of internal friction = 3 3 0 A simulated sinusoidal base acceleration, 0.5 g Sin (6.28t), was imposed at the till. The number of finite elements for soil and pile were 30 and 15 respectively. Boundary condition for the pile was such that axial displacement at the pile tip was zero. The relative displacement time history of the structural mass and the pile at the surface is shown in Fig. 4.12. It also includes the relative displacement time history of the soil at ground level. Figs. 4.13 and 4.14 show the profiles of relative displacement and bending moment in the pile at selected intervals. 4.6 R E S U L T S O F A N E X A M P L E P R O B L E M A T C A N L E X SITE, RICHMOND To simulate an actual situation another example was considered at the site in Richmond, British Columbia. The soil properties of the site are tabulated in Table 4.2. The same cross section used in the previous example was tested herein. The embedded length of the pile was 48.83 m and length above the ground level was assumed to be 2.4 m. Structural weight was assumed to be 500 k N . The boundary condition for the pile was such that axial displacement at the pile tip was zero. The relative displacement of the pile is relative to the till and it is not constrained (as i f the tip of the pile were placed on rollers). A n earthquake acceleration (San Fernando, 1971) was applied at the till, and the boundary condition for the soil was such that the relative displacement at the till was zero. 74 time = 0.25 sec Base Acceleration = 0.5 g Sin(6.281) Example No. 1 -3.00 0.00 3.00 6.00 \ time = 0.5 sec time = 0.75 sec : \ / ; \ 9.00 ; \ / / \j - 1 f — —- 1200 - Q 15.00 - \ ; - \; 18.00 21.00 24.00 11 27.00 li li — II 30.00 l -0.10 -0.05 0.00 Relative Displacement (m) i 0.05 0.10 Fig. 4.13. Relative Displacement Profile at Selected Time Steps 75 -3.00 Base Acceleration = 0.5 g Sin(6.28 t) Example No. 1 time = 0.25 sec — — - time = 0.5 sec time = 0.75 sec — — - time = 1.25 sec 0.00 3.00 a t Q 6.00 9.00 12.00 -8000.00 -4000.00 0.00 Bending Moment (kN. m) 4000.00 Fig. 4.14 Bending Moment Profile at Selected Time Steps for Example Problem 1 76 The input data file is enclosed in Appendix TJ. Results of relative displacement time history of the structural mass and the pile at ground surface is shown in Fig. 4.15. Profiles of the pile relative displacement and bending moment are shown from Fig. 4.16 and Fig. 4.17. From Fig. 4.15 it is observed that, when the amplitude of base acceleration is small at the end of the quake, the soil-pile-superstructure system undergoes free vibration. The natural period of vibration may be obtained from this time history plot. It could be observed from Fig. 4.17 that during an earthquake a pile would have a maximum bending moment near the surface, approximately at about 1/7^-1/^ the embedded depth of pile below the ground level. The reason could be that as gapping forms the deformation of the pile could be large as compared to the situation with no gapping at all. This is a crude approximation as it depends on number of factors, such as the cross section and material of pile, embedded depth, type of the soil deposit and frequency and amplitude of earthquake excitation. 77 -5.00 time = 2.0 sec time = 4.0 sec time = 6.0 sec time = 12.0 sec Base Acceleration = San Fernando Earthquake, 1971 Example No. 2 \ 0.00 \\ s 5.00 / 10.00 y —— yyy t\ \ 15.00 \K !\ \ - 1 \ 1 \V 20.00 \, V 1 1 j i\y 1 25.00 30.00 / / s y y / f y y y t• 1 35.00 \ 40.00 \ \ 45.00 50.00 I 1 \ 1 \ \ i * 1 1 1 1 1 I 1 -0.1000 -0.0500 0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 Relative Displacement (m) Fig. 4.16. Relative Displacement Profile of Pile at Selected Time Steps for Example Problem 2 79 -5.00 Base Acceleration = San Fernando Earthquake, 1971 Example No. 2 time = 2.0 sec time = 4.0 sec time = 6.0 sec 0.00 5.00 10.00 15.00 20.00 +^4 25.00 30.00 35.00 40.00 45.00 50.00 -4000.0000 -2000.0000 0.0000 Bending Moment (kN. m) 2000.0000 Fig. 4.17 Bending Moment Profile of Pile at Selected Time Steps for Example Problem 2 80 CHAPTER 5 CONCLUSIONS AND SCOPE FOR FURTHER RESEARCH A finite element seismic response analysis procedure for soil-pile-superstructure interaction is developed incorporating the geometric and material non-linearities of the pile, non-linear near-field-pile relationship, P-A effects and the soil shear stress-strain nonlinearities. The formulations for free-field and the soil-pile interaction are based on the principle of virtual work. The formulation has been incorporated in two finite element programs : Q U I V E R and P I L E U B C . Finite element predictions were compared with linearized closed form solutions wherever possible. In particular, a linearized closed form solution was developed for shear wave propagation in a uniform stress field subjected to a base acceleration. This solution was in exact agreement with the finite element analysis. The significance of the procedure presented in this thesis is that the usual assumptions made in other procedures like, "sub-structuring" of the soil-pile- superstructure and hence a equivalent linear analysis has been avoided. Damping in this dynamic system is obtained purely through hysteretic energy dissipation, by following through the non-linear constitutive relations for the pile and soil. Soil-pile separation was considered, keeping track of whether there was contact or not between the soil and the pile throughout the depth. In this formulation an assumption was made that there is no rise in pore pressure when a deposit with sandy soil is subjected to an earthquake excitation. This assumption could be invalid, as has been observed that the maximum dynamic shear modulus of sandy 81 soil degrades in an undrained situation during excitation. However, it could be incorporated in the finite element program Q U I V E R using shear-volume coupling model developed by Byrne, 1990. B y incorporating this model in a subroutine called H Y S T E R in the Q U I V E R program, the seismic response of piles on liquefied soil could be studied. The effect of accounting for uncertainties in the geometry or material properties of the pile and soil could be studied using stochastic finite elements. Uncertainties may exist in, for example, the initial elastic modulus of the pile or the shear modulus in the soil which are usually spatially distributed over the region of the soil-pile system and could be modeled as a spatial random process. Since pile foundations areto be designed with a high level of reliability, it would be worthwhile to study effects of such spatial variability of material properties. The analysis developed in this thesis could thus be coupled with a general reliability analysis for earthquake loading, including the uncertainties in the seismic excitation itself. 82 REFERENCES [I] Abe, et. al., (1984), " Dynamic Behaviour of Pile Foundations During Earthquakes," Proceedings on 8th World Conference on Earthquake Engineering, Vol.3, pp. 585-592. [2] American Petroleum Institute, (1979), "Recommended Practice for Planning, Designing and Constructing Fixed Offshore Platforms," 10th Edition , Washington, D.C. [3] Banerjee, P. K., and Davies, T. G., (1977), "The Behaviour of axially and laterally loaded single piles embedded in non-homogeneous soils," Geotechnique 28, No. 3. [4] Bea, R. G., (1980), "Dynamic Response of Piles in Offshore Platforms," Proceedings - Dynamic Response of Pile Foundations: Analytical Aspects, A S C E National Convention, ed. M . W. O'Neill and R. Dobry. [5] Blaney, et. al., (1976), "Dynamic Stiffness of Piles," Proc. 2nd Int. Conference on Numerical Methods in Geomechanics, A S C E , Blacksburg, V A , 1001-1012. [6] Blaney, G. W., and O'Neill, M . W., (1986), "Analysis of Dynamic Laterally Loaded Pile in Clay," Journal of Geotechnical Engineering Division, Vol. 112, No. 9, pp. 827-841. [7] Byrne, P. M . , and Mclntyre, J., (1994), "Deformations in Granular Soils due to Cyclic Loading," Settlements '94, A S C E Specialty Conference, Texas, U S A . [8] Byrne, P. M . , (1994), "Lecture Notes - C I V L 581 Soil Dynamics," Department of Civil Engineering, The University of British Columbia, Vancouver, Canada. [9] Byrne, P. M . , (1990), " A Cyclic Shear-Volume Coupling and Pore Pressure Model for Sand," Soil Mechanics Series No. 145, Department of Civil Engineering, The University of British Columbia, Vancouver, Canada. [10] Duncan, J. M . , et. al., (1980), "Strength, Stress-Strain and Bulk Modulus Parameters for Finite Element Analyses of Stresses and Movements in soil Masses," Geotechnical Research Report No. UCB/ GT/ 80-01, University of California, Berkeley. [II] Finn, W. D . L and Gohl, W. B . , (1987), "Centrifuge Model Studies of Piles Under Simulated Earthquake Loading," Dynamic Response of Pile Foundations - Experiment, Analysis and Observation, A S C E Convention, Atlantic City, New Jersey, Geotechnical Special Publication No. 11, pp. 21-38. [12] Fan, K . , et. al., (1991), "Kinematic Seismic Response of Single Piles and Pile Groups," Journal of Geotechnical Engineering Division, A S C E , V o l . 117, No. 12, pp. 1860-1879. 83 [13] Foschi, R. O., (1988), "User's Manual: R E L A N (RELiability Analysis)," Civil Engineering Department, University of British Columbia, Vancouver, Canada. [14] Gazetas, G., (1984), "Seismic Response of End-Bearing Single Piles," Soil Dynamics and Earthquake Engineering, Vol. 3 (2), 82-93. [15] Gazetas, G., et. al., (1992), "Seismic Pile-Group-Structure Interaction," Piles Under Dynamic Loads," Geotechnical Special Publication No. 34, A S C E , New York, New York, pp. 56-93. [16] Gohl, W. B . , (1991), "Response of Pile Foundations to Simulated Earthquake Loading : Experimental and Analytical Results," Ph. D. Thesis, Department of Civil engineering, The University of British Columbia, Vancouver, B . C . , Canada. [17] Hardin, B . O., and Drenvich, V . P., (1972), "Shear Modulus and Damping in Soils, Design Equations and Curves," Journal of Soil Mechanics and Foundation Division, A S C E , Vol. 98, SM7, Proc. Paper 9006, July, pp. 667-692. [18] Humar, J. L . , (1990), "Dynamics of Structures," Prentice Hall, Englewood Cliffs, New Jersey, U S A . [19] Ishihara, K . , (1982), "Evaluation of Soil Properties for Use in Earthquake Response Analysis," International Symposium on Numerical Models in Geomechanics, Zurich, pp. 237-259. [20] Kuhlemeyer, R. L . , (1979), "Static and Dynamic Laterally Loaded Floating Piles," Journal of the Geotechnical Division, Vol. 105, GT2, pp. 289-304. [21] Kagawa, T., and Kraft, L . M . , (1980), "Lateral Load Deflection Relationships of Piles Subjected to Dynamic loadings," Soils and Foundations, Vol. 20, No. 4, pp. 19-36. [22] Kaynia, A . M . , and Kausel, E . , (1982), "Dynamic Behavior of Pile Groups," Proceedings of Conference on Numerical Methods in Offshore Piling, The University of Texas, Austin, Texas, pp. 509-532. [23] Konder, R. L . , (1963), "Hyperbolic Stress-Strain Response; Cohesive Soils," Journal of Soil Mechanics and Foundation Division, A S C E , Vol. 89, No. S M I . [24] Luco, J. E . , (1982), "Linear Soil-Structure Interaction : A Review," Earthquake GroundMotion and Effects on Structures, A S M E , A M D , Vol. 53, pp. 41-57. [25] Liou, D . D . , and Penzien, J., (1977), "Seismic Analysis of an Offshore Structure Supported on Pile Foundations," EERC Report 77-25, University of California, Berkeley. 84 [26] Makris, N . , and Gazetas, G., (1992), "Dynamic Pile-Soil-Pile Interaction. Part U : Lateral and Seismic Response," Journal. Earthquake Engineering and Structural Dynamics, Vol. 21, pp. 145-162. [27] Masing, G., (1926), "Eigenspannungen and Verfestigung Beim Messing," Proceedings, 2nd International Congress of Applied Mechanics, Zurich, Switzerland. [28] Matlock, H . , and Reese, L . C , (1960), "Generalized Solutions for Laterally loaded Piles," Journal of the Soil Mechanics and Foundations Division, A S C E , Vol. 86, No. SM5, pp. 63-91. [29] Murchison, J. M . , and O'Neill, M . W., (1984), "Evaluation of P - Y Relationships in Cohesionless Soils," from Analysis and Design of Pile Foundations, A S C E Geotechnical Engineering Division, A S C E National Convention, San Francisco, pp. 174-192. [30] Newmark, N . M.,(1959), " A Method of Computation for Structural Dynamics," Journal of Engineering Mechanics Division, A S C E , Vol. 85, pp. 67-94. [31] Nogami, T., (1985), "Non-Linear Dynamic Winkler Model for Lateral Cyclic Response Analysis of Single Piles," Proceedings of 2nd International Conference on Soil Dynamics and Earthquake Engineering, Vol. 4, pp. 51-61. [32] Nogami, T., and Novak, M . , (1977), "Resistance of Soil to a Horizontally Vibrating Pile," International Journal of Earthquake Engineering and Structural Dynamics, Vol. 5, pp. 249-261. [33] Novak, M . , (1974), "Dynamic Stiffness and Damping of Piles," Canadian Geotechnical Journal, National Research Council of Canada, Toronto, Ontario, Canada, Vol. 11, No. 4, pp. 574-598. [34] Novak, M . , and Aboul-Ella, F., (1978), "Impedance Functions of Piles in Layered Media," Journal of the Engineering Mechanics Division, A S C E , Vol. 104, No. E M 6 , pp. 643-661. [35] Ogata, N . , and Kotsubo, S., (1966), "Seismic Force Effect on Pile Foundation," Proceedings. Japan Earthquake Engineering Symposium, Tokyo, Japan. [36] Poulous, H . G., and Davis, E . H , (1980), "Pile Foundations Analysis and Design," John Wiley and Sons, Inc., I S B N 0-471-02084-2. [37] Powers, D . L . , (1989)," Boundary Value Problems," 3rd Edition, Harcourt Brace Jovanovich Publishers, Orlando, Florida, U S A . [38] Prakash, S., and Chandrasekaran, V . , (1973), "Pile Foundations under Dynamic loads," Symposium on Behavior of Earth and Earth Structures Subjected to Earthquakes and Other dynamic Loads, Roorkee, India, Vol. 1, pp. 165-173. 85 [39] Ramberg, W., and Osgood, W. R , (1943), "Description of Stress-Strain Curves by Three Parameters," Technical Note 902, N A C A . [40] Rosenbleuth, E . , (1980), "Design of Earthquake Resistant Structures," Pentech Press, Plymouth, London, pp. 223-260. [41] Rosset, J. M . , and Angelides, D., (1980), "Dynamic Stiffness of Piles," Numerical Methods in Offshore Piling, Institution of Civil Engineers, pp. 75-81. [42] Scnabel, P. B . , et. al., (1972), " S H A K E - A Computer Program for Earthquake Response Analysis of Horizontally Layered Sites," Report No. EERC 72 - 12, University of California, Berkeley. [43] Seed, R.B., and Harder, L.,(1992), "SPT-Based Analysis of Cyclic Pore Pressure Generation and Undrained Residual Strength", 'H. Bolton Seed', V o l II, Selected Papers, 1957-1987, BiTech Publishers, Vancouver, B . C , Canada [44] Scott, C. R , (1980), " A n Introduction to Soil Mechanics and Foundations," 3rd Edition, Applied Science Publishers Ltd., London. [45] Seed, H . B . , and Idriss, I. M . , (1970), "Soil - Moduli and Damping Factors for Dynamic Response Analysis," Report No, EERC 70 - JO, University of California, Berkeley. [46] Sy, A . , and Siu, D., (1992) "Forced Vibration Testing of an Expanded Base Concrete Pile," Piles Under Dynamic Loads, Geotechnical Special Publication No. 34, A S C E , New York, New York, pp. 170-186. [47] Tajimi, H . , (1969), "Dynamic Analysis of a Structure Embedded in an Elastic Stratum," Proceedings, Fourth World Conference on Earthquake Engineering, Santiago, Chile, Vol. 3, pp. 54-69. [48] Vazhinkhoo, S. (1996), "Modeling of Single Piles Subjected to Monotonic and Cyclic Lateral Loads Including Free-Field Displacements", Forthcoming, Master's Thesis, The University of British Columbia, Vancouver, Canada. [49] Velez, A., et. al., (1983), "Lateral Dynamic Response of Constrained Head Piles," Journal of Geotechnical Engineering, Vol. 109, No. 8, pp. 1063-1081. [50] Wu, G., (1994), "Dynamic Soil-Structure Interaction : Pile Foundations and Retaining Structures," Ph. D. Thesis, The University of British Columbia, Vancouver, Canada. [51] Yan, L . , and Byrne, P. M . , (1991), "Lateral Pile Response to Monotonic Pile Head Loading,",Canadian Journal of Civil Engineering. 86 [52] Yoshida, L . , and Yoshinaka, R., (1972), " A Method to Estimate Modulus of Horizontal Subgrade Reaction for Pile," Soils and Foundation, Vol. 12, No. 3, pp. 1-17. [53] Zienkiewicz, O. C , (1977), "The Finite Element Method," 3rd Edition, McGrawHill Book Company(U K) Limited, Berkshire, England. 87 APPENDIX I ANALYTICAL SOLUTION FOR A SINGLE LAYER SOIL SYSTEM USING LAPLACE TRANSFORMS The Laplace transform serves as a device for obtaining the solution of ordinary or partial differential equations (Powers, 1989). It associates a function /ft) with a function of another variable, F(s), Equation A. 1, from which the original function can be recovered. (A.1) Fig. 6 shows a soil particle of density p, undergoing a horizontal translation u(x,t). For an equilibrium to exist in this assumed linear undamped system with G as the maximun dynamic shear modulus, the following one dimensional wave propagation equation has to be satisfied + a(t) dx 1 (A.1) or a(0 + dt 1 =c dx 1 (A.2) In which, p = mass density of soil, G = shear modulus of the soil, c = shear wave velocity of the layer, aft) = harmonic base acceleration, 88 -x = 0 cfu p(-2 at + a(t)) O 1 x=G dx H T + 0 dr dx TiU a(t) Fig. A. 1.1 Soil deposit as a one dimensional shear wave propagation medium 89 u = relative displacement with respect to till. We want to find the Laplace transform of the one dimensional wave equation whose appropriate initial and boundary conditions are, u(x,0) = 0 u(x,0) = Y 0 =|>,t; = o (A.3) ox. u(H,t) = 0 Using the properties of Laplace transformation , we transform both sides of the equation, obtaining a(s) + u(s) = c 2 s 1 — u ( s ) ( A 4 ^ or * i(,)-4i -»a W c c (A.5) Solving Equation A.5 u ( s ) = C e" + C e - ' - ^ and x x l ^ 2 IT =G^'-GV** O X c (A.6) (A.7) c where C i and C2 are arbitrary constants. Using the boundary conditions that the relative displacement is zero at base and the shear strain at the surface is zero, the constants Ci and C2 are evaluated : 90 Therefore the transformed relative displacement response at any point can be evaluated, 2nr\ u(s)-—- e s e + (\+ -< y 2 a(s) \~ e (A.9) > H e or, u(s) = If z = e ~ then y 2sH/c B> (A. 10) -1.0 s 2 (1+e° ) H = ( - l ) " / i ! (1+z) expansion would be sufficient to numerically obtain the n derivative of y. Substituting this in Equation A . 10 as th 0«-7") =„? <- „ _2mK e 1) 0 a (s) sr e c l«=o ^0 v (All) c (A. 12) 1.0 +e As indicated earlier, it is now possible to recover the original function in the time domain «(x,0=H-l)" , ,lrH 1H - c c c) A<-\—+• A / +j t + \ c c c) +- \ c cJ b t-\ +- V c c) -fit) (A. 13) In which f(t) is the base displacement (excitation). For instance, say, we induce a sinusoidal acceleration #(f) = A sin(or), the corresponding displacement would be At A f(t) = — —sin(at). CO Where A is amplitude and co is the frequency of excitation. This will be the function f(t) to be used in the time domain solution. Equation A . 13, as pointed out earlier, would be the closed form expression for the relative displacement of a undamped single layer soil deposit with constant Gmax thoughout the depth. This involves both transient and steady state solutions. Also this expression holds 91 good for obtaining the response history at anywhere in the soil profile. If the surface relative displacement is of interest, then x = H in Equation A. 13, therefore (inH c if ... c ,2nH H c c -fit) (A.14) or, (A.15) This solution procedure is compared with the finite element program and it is observed that the finite element solution coincides with exact solution. 92 APPENDIX II PROGRAM SETUP AND INPUT DATA FILES A.2.1 INTRODUCTION The dynamic finite element computer codes Q U I V E R and P I L E U B C was written in F O R T R A N 77 and compiled with Microsoft® Fortran with Microsoft Visual Work Bench® Version 3.2. The machine used for developing the program was an IBM® compatible P C with 486DX33 micro processor. The minimum hardware required to run P I L E U B C are: • a 386 based or higher central processing unit. • a math co-processor • at least 12.0 Mb of free hard disk space (depends on the size of the problem) • 8.0 Mb of random access memory (RAM) Due to large number of arrays defined in the program, the compiler is capable of making an executable file such that memory allocation for all arrays is done at run time rather than including it in the memory of the executable file. The program has been tested with a wide range of cases and is believed to be free from serious defects. Troubles are usually found to be caused by user-oriented errors in the input files or misrepresentation of the physical system. This in return, could result in unexpected response of the soil-pile-superstructure system. The program has been designed to work with any consistent system of units by specifying the value of atmospheric pressure. On account of the integration scheme employed and the nodal 93 degrees of freedom relatively fewer elements may be used to represent the pile and the free-field soil for predicting the time varying response. The general flow of the coupled soil-pile-superstructure seismic response program (PILEUBC and QUIVER) are illustrated in Fig. A.2.1 and A.2.2 A.2.2 P R O G R A M STRUCTURE The program P I L E U B C consists of a main routine and the following fifteen subroutines : QUIVER, SHAPES, SHAPES1, SHAPES2, G A U S S , GAUSS1, SOMASS, E M A S S , STRESS, PSUP, H Y S T E R , S O A C C V E C , P I L A C C V E C , D E C O M P , and SOLV. The main routine is responsible for all input and output, initialization of arrays, construction of local and global stiffness matrices, construction of global mass matrix, inclusion of boundary conditions for the pile and near-field soil. Subroutine Q U I V E R computes the free-field soil profile as a function of space and time. Subroutine SHAPES obtains the values of shape functions at each Gaussian sampling point for each pile element in a layer. These are used in construction of local stiffness matrix for the pile element. SHAPES 1 is similar to SHAPES, but it calculates the values of shape functions for freefield soil elements. Subroutine GAUSS return appropriate Gaussian coordinates and weights for desired number of Gaussian integration points for pile elements. G A U S S 1 does it for free-field soil elements. Subroutine SOMASS computes the local consistent mass matrix for the free-field. E M A S S returns back the same for the pile elements. Subroutine STRESS looks up the appropriate axial and bending stresses for a given strain level at a Gaussian point in the pile. PSUP computes the tangent stiffness and pressure for 94 READ INPUT DATA FOR SOIL & PILE INITIALIZE ARRAYS AND OBTAIN GAUSSIAN COORDINATES AND WEIGHTS CALCULATE PILE ELEMENT MASS MATRIX AND ACCELERATION VECTOR AND GLOBALIZE IT CALL GAUSS & SHAPES CALL EMASS & PTLACCVEC REPEAT FOR EACH TIME STEP CALL QUIVER READ INPUT BASE ACCELERATION FROM FILE CONSTRUCT RHS f(t) INVOKE NEW MARK - Beta PROCEDURE CONSTRUCT DERIVATIVE MATRIX FOR THE NEWTON RAPHSON PROCEDURE CONSTRUCT PSI VECTOR IMPOSE BOUNDARY CONDITIONS DECOMPOSE AND SOLVE FOR INCREMENTAL DISPLACEMENT CALL DECOMP & SOLV CALL STRESS & PSUP YES EXTRACT FREE-FIELD RELATIVE DISPLACEMENTS -0- UPDATE STRESSES IN THE PILE AND SOIL PRESSURES ACTING ON THE PILE FOR RECENTLY OBTAINED DISPLACEMENTS CALCULATE FORCES, MOMENTS, VELOCITY, ACCELERATION VECTORS PRINT OUTPUT INTO PRECRIBED FILES ( END ) Fig. A.2.1 Generalized Flow Chart for Finite Element Program, PILEUBC 95 INITIALIZE ARRAYS AND OBTAIN GAUSSIAN COORDINATES AND WEIGHTS CALCULATE SOIL ELEMENT MASS MATRIX AND ACCELERATION VECTOR AND GLOBALIZE IT CALL GAUSS1 & SHAPES1 CALL SOMASS & SOACCVEC CONSTRUCT RHS f(t) INVOKE NEW MARK - Beta PROCEDURE CONSTRUCT DERIVATIVE MATRIX FOR THE NEWTON RAPHSON PROCEDURE CONSTRUCT PSI VECTOR IMPOSE BOUNDARY CONDITIONS DECOMPOSE AND SOLVE FOR INCREMENTAL DISPLACEMENT CALL DECOMP & SOLV UPDATE SHEAR STRESSES AND TANGENT MODULUS TN FREE-FIELD SOIL FOR RECENTLY OBTAINED DISPLACEMENTS CALCULATE, VELOCITY, ACCELERATION VECTORS PRINT OUTPUT INTO PRECRIBED FILES RETURN TO THE MAIN ROUTINE Fig. A.2.2 Generalized Flow Chart for Finite Element Program, Q U I V E R 96 a particular p~y curve depending on the amount of lateral movement of the pile in the near field. Subroutine H Y S T E R gives the tangent shear modulus and shear stress in the soil at a sampling point in the free-field element, based on the Hardin-Drenvich model. Subroutine S O A C C V E C calculates the acceleration vectors for the driving forces in free-field element, locally. P I L A C C V E C does for the pile element. Subroutine D E C O M P decomposes the one dimensional derivative of the matrix obtained from the truncated Taylor series for the Newton Raphson procedure. Subroutine S O L V solves the system of algebraic equations to get the incremental displacement of the system. Both D E C O M P and S O L V are used in free-field response and pile response. A.2.3 C U R R E N T P R O G R A M CAPABILITIES The program Q U I V E R and P I L E U B C are capable of analyzing any of the combination of the following : 1. Total stress free-field soil response subjected to base acceleration at the till or any point in the soil profile. It incorporates shear stress-strain nonlinearities in the soil. Also includes degradation or upgradation of the dynamic stiffness of free field depending on the loading, unloading or reloading conditions. 2. Different soil layers in the deposit. 3. Soil nonlinearities and yielding in the near-field through p~y curves. Also includes degradation or upgradation of the dynamic stiflhess of near field depending on the loading, unloading or reloading conditions. 4. Automatic computation of p~y curves based on either API or Yan-Byrne method. 97 5. Soil-pile gapping at the interface. 6. Earthquake loading. 7. Varying pile cross section along the length of the pile. 8. Different pile materials along the length of the pile. 9. Pile material nonlinearity, yielding, and geometric nonlinearity. 10. Direct loading at any node along the length of the pile. 11. P-A effects from applied axial load. 12. Non-linear vibration of a beam. A.2.4 DESCRIPTION O F INPUT D A T A FILES The input data files can be given in free format. The example files enclosed (Section A.2.4) contain the data using comma separated variables. As mentioned earlier any system of units may be employed by specifying the atmospheric pressure in the units of choice. However, other data such as specific unit weight, depth of soil layer, Young's modulus and yield stress of pile material must be in consistent units. Automatic node number generation is incorporated. It goes from bottom to top in the finite elements. The following describes the data file for soil deposit: Line No. Variable Name Format Description 1 TITLE A80 TITLE FOR THE SOIL DEPOSIT 2 NLAYERS 14 NUMBER OF LAYERS IN THE DEPOSIT 3 SDEPTH(NLAYERS) E15.6 DEPTH OF L A Y E R 98 3 NSOIELE(NLAYERS) 14 3 DENSITY(NLAYERS) E15.6 DENSITY OF THE L A Y E R 3 R N 1 (NL A Y E R S ) F15.6 NORMALIZED SPT BLOW COUNT 3 PHI(NLAYERS) F15.6 SOIL FRICTION ANGLE 3 SU(NLAYERS) F15.6 COHESION INTERCEPT NUMBER OF ELEMENTS IN T H A T LAYER REPEAT THE ABOVE FOR EACH L A Y E R NUMBER POINTS OF GAUSSIAN SAMPLING 4 NGS 14 5 STOL F15.6 TOLERANCE FOR INCREMENTAL RELATIVE DISPLACEMENT AND T O T A L ERROR. 6 NLIN 14 7 SDAMP F15.6 OPTION FOR LINEAR OR NONLINEAR ANALYSIS = 0 FOR LINEAR ANALYSIS = 1 NONLINEAR ANALYSIS MASS PROPORTIONAL RALEIGH DAMPING COEFFICIENT The data file for the earthquake acceleration should be made available in following format: Line No. Variable Name Format Description of Variable 1 TITLE2 A80 TITLE FOR T H E INPUT ACCELERATION RECORD 2 TSTEP F15.6 3 NPOINTS 18 NUMBER OF D A T A POINTS 4 GRACC(NPOINTS) E15.6 ACCELERATION IN CONSISTENT UNITS FOR NPOINTS TIME STEP FOR THE ACCELERATION RECORD BASE The data file P I L E U B C should be created as indicated below : Line No. Variable Name Format Description 1 TITLE A80 TITLE FOR THE PROBLEM 2 ELEMENTS 14 NUMBER OF ELEMENTS A L O N G T H E PILE 99 NUMBER OF DIFFERENT SOIL LAYERS 2 SOILLAYERS 14 2 PILEMATS 14 NUMBER OF MATERIALS 2 BCJSTODES 14 NUMBER OF DIFFERENT NODES WITH BOUNDARY CONDITIONS 2 INPUT PTS 14 NUMBER OF LOAD INPUT POINTS, IF ANY, ENTER 0 FOR SEISMIC LOADING 2 IS_CYC 14 2 FREEFLDNODES 14 2 PA F15.6 = 'Y' IF LOADING IS CYCLIC = 'N' OTHER WISE, USED IN CALCULATIONS FOR API CURVES NUMBER OD NODES ASSOCIATED WITH LATERAL STATIC LOADING, ENTER 0 FOR EARTHQUAKE LOADING ATMOSPHERIC PRESSURE IN DESIRED UNITS 2 STATICVERLOAD F15.6 WEIGHT OF THE STRUCTURAL MASS 2 TOLERANCE F15.6 TOLERANCE FOR INCREMENTAL DISPLACEMENT AND TOTAL ERROR 3 NODE 14 3 X_COORD(NODE) F15.6 NODE NUMBER. NODES ARE NUMBERED FROM BOTTOM TO TOP. THE TOTAL NUMBER OF NODES IS ALWAYS ONE PLUS TOTAL NUMBER OF ELEMENTS COORDINATE OF THE NODE IN THE VERTICAL DIRECTION 3 IS_FF(NODE) 12 = '0' FREE FIELD = '1'OTHERWISE ENTER '0' FOR EARTHQUAKE LOADING DIFFERENT PILE REPEAT THE ABOVE FOR EACH NODE 4 ELEM 14 ELEMENT NUMBER 4 XSEC_TYPE(ELEM) 12 4 OUT_DIA(ELEM) F15.6 THE SHAPE OF THE PILE CROSS SECTION, THE CHOICES ARE : = 1 - SOLID CIRCULAR = 2 - HOLLOW CIRCULAR =3 -RECTANGLE IFXSEC TYPE0=1: = DIAMETER OF PILE IFXSEC TYPE0 = 2: = OUTER DIAMETER OF PILE IFXSEC TYPE() = 3: = DEPTH OF PILE (DIRECTION OF BENDING) 4 LN_DIA(ELEM) F15.6 IFXSEC TYPE() = 1: = 0.0 IFXSEC TYPE0 = 2: = INNER DIAMETER OF PILE IFXSEC TYPE0 = 3: = WIDTH OF PILE 4 R P M A S S(ELEM) F15.6 MASS DENSITY OF THE PILE MATERIAL 100 4 MAT_NUM(ELEM) 14 PILE MATERIAL NUMBER FOR THIS ELEMENT 4 IN_LAYER(ELEM) 14 SOIL LAYER ASSOCIATED WITH THIS ELEMENT REPEAT THE ABOVE FOR EACH ELEMENT 5 LAYER 14 SOIL LAYER NUMBER 5 SOIL_TYPE(LAYER) A4 5 PY_TYPE(LAYER) A4 SOIL TYPE FOR THE LAYER: = CLAY = SAND = USER THIS OPTION DEFINES HOW TO CALCULATE P-Y CURVES FOR THE SOIL LAYER METHOD OF CALCULATING P-Y CURVES: = APIC, FOR APPROACH USED IN THE API CODE; = YANB, FOR YAN-BYRNE APPROACH = USER, FOR USER DEFINED P-Y CURVES IF PY_TYPEO = APIC, ENTER THE FOLLOWING: 6 GAMMA(LAYER) F15.6 THE UNIT WEIGHT OF SOIL FOR THIS LAYER 6 DR(LAYER) F15.6 THE RELATIVE DENSITY OF SOIL FOR THIS LAYER 6 ETA(LAYER) F15.6 THE FACTOR ETA FROM API CODE FOR THIS SOIL TYPE 6 N_HI(LAYER) F15.6 6 CI (LAYER) F15.6 THE FACTOR n hi FROM API CODE FOR THIS SOIL TYPE (COEFFICIENT OF SUBGRADE REACTION MODULUS) FACTOR CI FROM API CODE 6 C2(LAYER) F15.6 FACTOR C2 FROM API CODE 6 C3 ( L A Y E R ) F15.6 FACTOR C3 FROM API CODE IF PY_TYPE0= YANB, ENTER THE FOLLOWING 6 GAMMA(LAYER) F15.6 THE UNIT WEIGHT OF SOIL FOR THIS LAYER 6 DR(LAYER) F15.6 THE RELATIVE DENSITY OF SOIL FOR THIS LAYER 6 EMAX(LAYER) F15.6 OPTIONAL, ENTER 0.0 IF AUTOMATIC COMPUTATION IS REQUIRED AS A FUNCTION OF EFFECTIVE PRESSURE F15.6 THE UNIT WEIGHT OF SOIL FOR THIS LAYER IF PY_TYPE()= USER, ENTER THE FOLLOWING 6 GAMMA(LAYER) 101 6 NUMPY(LAYER) 14 NUMBER OF POINTS TO BE ENTERED FOR THE P-Y CURVE FOR THIS LAYER 6 NODE_PY(LAYER,NUM_PY) 14 NODE AT WHICH THIS PY-TYPE IS USED 6 EMAX(LAYER) F15.6 OPTIONAL, ENTER = 0.0 AUTOMATIC CALCULATION 6 Y_PY_IN(NODE,J) F15.6 DISPLACEMENT COORDINATE, Y 6 P_PY_IN(NODE,J) F15.6 PRESSURE COORDINATE. P0 7 PILEMAT 14 PILE MATERIAL NUMBER 7 MAT_TYPE(PILE MAT) A4 7 E(PTLE_MAT) F15.6 DEFINES THE STRESS-STRAIN BEHAVIOR FOR PILE ELEMENT: = ELPL, FOR ELASTIC PERFECTLY PLASTIC CURVE = WOOD, FOR TIMBER PILES YOUNG'S MODULUS FOR THE MATERIAL 7 YIELD_STRS(PILE_MAT) F15.6 YIELD STRESS OF THE MATERIAL FOR THE ELEMENT. FOR IF MATTYPEO = WOOD, ENTER FOLLOWING 7 E(PLLE_MAT) F15.6 YOUNG'S MODULUS 7 YTELD_STRS(PILE_MAT) F15.6 YIELD STRESS OF THE MATERIAL IN COMPRESSION IF BC NODES IS NOT ZERO, FOR EACH BCNODES ENTER THE FOLLOWING 8 BC_NODE(I) 14 NODE NUMBER WITH BOUNDARY CONDITION, I = 1, BCNODES 9 NUM_BCS(BC_NODE(I)) 14 NUMBER OF DIFFERENT BOUNDARY CONDITIONS TO BE SPECIFIED 9 B C ( B C NODE(I), J, J = l , NUM_BCS(BC_NODE(I))) 14 = 1,IFW= 0 = 2,IFW' = 0 = 3,IFW" = 0 = 4,IFU= 0 = 5,IFU'=0 N O T E : T H E N O D E N U M B E R S SPECIFIED H E R E M U S T E X A C T L Y COINCIDE W I T H T H A T G I V E N I N T H E INPUT D A T A F I L E 102 A.2.5 S A M P L E INPUT D A T A FILES The input data files for a few examples are given below : 4 LAYER 15 ELEMENTS INPUT 4 .325E-K)1,1,0.1886E+02,6.0D0,0.0D0,89.769 .325E+01,1,0.905E+01,11.0D0,0.0D0,180.00 .1625E+02,5,0.984E+Ol,24.0D0,41.0D0,0.00 .26E+028,0.984E+01,23.0D0,40.0D0,0.00 4 1.0E-04 16 150 1 0.00 ( CANLEX SITE RICHMOND 1995 15,4,1,1,0,0,^,0,101.325,-500.0,0.001 1,48.75,0 14,0.0,0 15, -1.20,0 16, -2.40,0 I, 2,0.900,0.7,7.4,1,4 2,2,0.900,0.7,7.4,1,4 3,2,0.900,0.7,7.4,1,4 4,2,0.900,0.7,7.4,1,4 5,2,0.900,0.7,7.4,1,4 6,2,0.900,0.7,7.4,1,3 7,2,0.900,0.7,7.4,1,3 8,2,0.900,0.7,7.4,1,3 9,2,0.900,0.7,7.4,1,3 10,2,0.900,0.7,7.4,1,3 II, 2,0.900,0.7,7.4,1,3 12,2,0.900,0.7,7.4,1,2 13,2,0.900,0.7,7.4,1,1 14,2,0.900,0.7,7.4,1,0 15,2,0.900,0.7,7.4,1,0 l.'SANDVAPIC' 0.1886E+02,60.0,0.0 2, SAND7APIC' 0.905E-K)1,60.0,0.0 3, SAND , YANB 0.984E-K)1,70.0,0.0 4, SAND , APIC 0. 984E+01,60.0,0.0 1. 'ELPL' 0.210E+O9,0.300E+O6 1,1,4 1,10 10,100,105,110,115,120,130,140,150,200,210 , , , , , , , , , 103 A.2.6 P R O G R A M O U T P U T There are 22 output files. Three set of five files consists of relative displacement, velocity and acceleration time histories at selected nodes. Five files consists of time varying bending moment at selected elements. A file contains the profile of the pile at selected time steps. The output consists of time step number, lateral deflection, shear force, bending moment and curvature. Also a file consists of absolute maximum bending moment in the pile. With very little modification in the program other desired data could also be obtained. The CONFIG.SYS file should be modified when run in batch mode to accomadate file management. A.2.7 P R O B L E M S E N C O U N T E R E D IN RUNNING P I L E U B C Although the program has been thoroughly checked and is believed to be free of errors, as it is a nonlinear finite element program, a few unforeseen errors might have been overlooked. In general, however, mostly problems occur in incorrect or incomplete input data files. Numerical inconsistencies depending on the sudden increase in successive earthquake acceleration (spikes), time step and very small tolerance may also create problems in attaining realistic computations. 104
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A new, non-linear analysis of layered soil-pile-superstructure...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A new, non-linear analysis of layered soil-pile-superstructure seismic response Khan, Anuz K. 1995
pdf
Page Metadata
Item Metadata
Title | A new, non-linear analysis of layered soil-pile-superstructure seismic response |
Creator |
Khan, Anuz K. |
Date Issued | 1995 |
Description | A non-linear finite element analysis for determining the seismic response of piles in multilayered soil deposit is presented. The analysis incorporates pile material and geometric non-linearities, P-A effects, non-linear load deformation behavior of soil-pile interaction, pile separation and the non-linear shear stress-strain relationship in the freefield soil. The soil shear stress-strain relation is expressed based on Hardin and Drenvich's model. The pile material stress-strain relation is based on a elastic-plastic model. Near field soil is replaced with a continuous foundation based on Yan-Byrne or American Petroleum Institute P-y curves. Time integration is performed utilizing Newmark-P integration scheme. Numerical integration of the virtual work equations has been carried out with Gauss quadrature. The resulting set of non-linear equations is solved by using the Newton-Raphson scheme. The method of analysis are incorporated in two programs : QUIVER and PILEUBC. The program QUIVER captures the seismic response of soil deposit and later is used as a module in the main program PILEUBC, in which the response of a single pile is predicted for an earthquake excitation. The program predicts time varying acceleration, velocity, displacement, bending/axial stress and strain, bending moment and shear force in the pile. The program is able to successfully simulate the seismic response of piles. Numerical investigations have been carried out to verify the results and test the capabilities of the program. A linearized closed form solution based on Laplace transformation for free-field is developed for comparison with the finite element analysis. |
Extent | 4232221 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-01-30 |
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. |
DOI | 10.14288/1.0050346 |
URI | http://hdl.handle.net/2429/4002 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1995-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1995-0613.pdf [ 4.04MB ]
- Metadata
- JSON: 831-1.0050346.json
- JSON-LD: 831-1.0050346-ld.json
- RDF/XML (Pretty): 831-1.0050346-rdf.xml
- RDF/JSON: 831-1.0050346-rdf.json
- Turtle: 831-1.0050346-turtle.txt
- N-Triples: 831-1.0050346-rdf-ntriples.txt
- Original Record: 831-1.0050346-source.json
- Full Text
- 831-1.0050346-fulltext.txt
- Citation
- 831-1.0050346.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-0050346/manifest