A TWO-DIMENSIONAL NON-LINEAR STATIC AND DYNAMIC RESPONSE ANALYSIS OF SOIL STRUCTURES By RAJARATNAM SIDDHARTHAN B.Sc. The University of Sri-Lanka, Peradeniya Campus, 1977 M.A.Sc. The University of B r i t i s h Columbia, 1981 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN THE DEPARTMENT OF CIVIL ENGINEERING We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA MARCH 1984 (c) Rajaratnam Siddharthan, 1984 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission f o r extensive copying of t h i s t h e s i s f o r scholarly purposes may be granted by the head of my department or by h i s or her representatives. I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of C l y\/^ max Tmax* These values should be updated as r e s i d u a l porewater pressure develops. Hard in , et a l . (1972) assumed that G m a x i s independent of s tress h i s t o r y and suggested, * 1/2 G = K (a ' ) ' (3.7) max m i n which K i s a constant , depends on s o i l type and r e l a t i v e d e n s i t y . The i n i t i a l and current e f f ec t ive stress condi t ions i n a simple shear apparatus with zero i n i t i a l porewater pressure are shown i n F i g . 3 .3 . Here i t i s assumed that the r a t i o between h o r i z o n t a l and v e r t i c a l F i g . 3.3(a). Initial Effective Stress Condition in a Simple Shear Apparatus. m i«2K»«r- -u) 3 F i g . 3.3(b). Intermediate Effective Stress Condition i n a Simple Shear Apparatus. Co 24 e f f e c t i v e stresses i s a constant K Q, where K Q i s the c o e f f i c i e n t of l a t e r a l pressure at r e s t . For the i n i t i a l e f f e c t i v e stress condition, ^ max^ l s 8 i v e n b v> 1/2 1 +9 V ' (G ) = K* f V > l (Ko)1/2 max 3 v ; o For the current stress conditions, 1/2 * 1+2K (G ) = K (if*) r;D- u) i / 2 <3-9> max v 3 ; ; n on d i v i d i n g the equation (3.9) by (3.8) one obtains, ( Gmax ) - U . n _ r vo i l / 2 .) L a' J (G o-' J (3.10) max vo o T h e r e f o r e , knowing ( ^ m a x ^ o ' a ' v o a n c * ^> t n e maximum s h e a r modulus at the current e f f e c t i v e stress condition can be calculated using equation 3.10. The shear s t r e n g t h ( T m a x ^ o ^ o r t n e I n i t i a l e f f e c t i v e stress condition i s given by (Finn, et a l . 1977). f ^ 1 + K 9 o i 1 - K 1/2 ( W = {(— - —°) s i n 2 0 - ( ^ - ^ ) 2 } 0' = C a' (3.11) o 1K 2 J K 2 ' ' vo vo 25 * i n which 0' i s the angle of i n t e r n a l f r i c t i o n and C i s a constant which depends on s o i l properties. For the current stress condition, ( W ) = C*(o-' - U) (3.12) n vo Dividing equation (3.12) by (3.11) one obtains, ( t ) max * ^ • (3.13) ( T ) a* max vo o ^ T m a x ^ o ' "vo and U a r e known, t h e maximum s h e a r strength at current e f f e c t i v e stress condition can be calculated from equation (3.13). 3.1.4 Influence of S t r a i n Hardening During seismic loading of dry sand or saturated sand under drained conditions, the sand structure hardens due to grain s l i p . Finn, et a l . (1977), used f o l l o w i n g equations to modify G m a x and T • max e v d (G ) = (G ) {1 + — ~ ~ } (3.14a) max max 1 H, + H„E ,J nn n 1 2 vd and (3.14b) 26 i n w h i c h , ( G m a x ) n n a n d ( T _ a x ) n n a r e t h e m o d i f i e d maximum shear modulus and shear strength i n the n th cyc le and , » H3 and are hardening constants . The s t r e s s - s t r a i n behaviour for one-dimensional a n a l y s i s i s now completely defined by equations, 3 .1 , 3.2, 3.4, 3 .5, 3 .6, 3.10, 3.13 and 3.14. In laboratory c y c l i c simple shear tests most of the volume changes i n dry sands and the increases i n porewater pressure i n undrained saturated sands occur during the unloading por t ion of the load c y c l e . Therefore , F i n n , et a l . (1977), used modi f icat ions to the s t r e s s - s t r a i n curve to take account of s t r a i n hardening and porewater pressure only during the unloading phases of the l o a d i n g . 3.1.5 D i s s i p a t i o n of Porewater Pressure I f the saturated sand deposit can dra in during shaking there w i l l be simultaneous generation and d i s s i p a t i o n of porewater pressure . The rate of increase of porewater pressure w i l l be le s s than that i n completely undrained sand. The amount of drainage depends on the permea-b i l i t y and c o m p r e s s i b i l i t y of the sand, drainage path and durat ion of shaking. The d i s t r i b u t i o n of porewater pressure at time t i s given by, _ . _ E j _ <_z_ _J-\ E 8 £ v d ot 5z W 5 z J ot 'w (3.15) 27 i n which, k z i s the permeabi l i ty and y w t n e un i t weight of water. The term c o n t a i n i n g e y ( j r e p r e s e n t s the i n t e r n a l g e n e r a t i o n of porewater pressure ( F i n n , et a l . 1977). The s t r e s s - s t r a i n r e l a t i o n s h i p out l ined above can be very e a s i l y extended to non-uniform loading using an incremental ly e l a s t i c ana lys i s i n time domain. The dynamic response can be computed for each time step by numer ica l ly so lv ing equation (3.15) and the equation of motion, as explained by F i n n , et a l . (1977). 3.2 LABORATORY VERIFICATION OF EFFECTIVE STRESS RESPONSE ANALYSIS The bas ic assumptions made i n the formulat ion of the s t re s s -s t r a i n r e l a t i o n s h i p presented above can be broadly categorized in to two groups: Those made i n the formulat ion of porewater pressure model and those made i n model l ing load ing , unloading and re load ing . The fundamental assumption that was made i n the formulat ion of the porewater pressure model, was that the porewater pressures i n an un-drained test can be obtained from volumetr ic s t ra ins measured i n a drained tes t on a s i m i l a r sample with same h i s t o r y of shear s t r a i n l o a d i n g . This means that there should be a unique r e l a t i o n s h i p between volumetr ic s t r a i n s i n drained tests and porewater pressures i n undrained tests for a given sand at corresponding s t r a i n h i s t o r i e s . F inn (1981) reported re su l t s of an extensive laboratory program to inves t iga te th i s bas ic assumption. Volumetric s t r a in s were measured i n drained Ottawa sand samples at r e l a t i v e dens i t i e s D r = 45% and 60% when subjected to constant s t r a i n cycles i n a simple shear apparatus. Porewater pressures were a lso measured i n undrained c y c l i c tes t s at the 28 same r e l a t i v e d e n s i t i e s and i n i t i a l e f f e c t i v e confining pressures. Volumetric s t r a i n s E V ( J are shown plotted against porewater pressure r a t i o s U/o^. i n F i g . 3.4 for D„ = 45%. Each point on the curve represents a set of values of e y (j and U/a v o f o r a given number of cycles with equal c y c l i c s t r a i n amplitudes. I t can be noticed that there i s a s l i g h t difference i n the applied shear s t r a i n amplitudes. But these small deviations are not important. The data i n d i c a t e a unique r e l a t i o n s h i p between volumetric s t r a i n and porewater pressure r a t i o s . The slope of t h i s curve normalized with respect to confining pressure w i l l give the rebound modulus E_. Martin, et a l . (1975) suggested that E_ can be evaluated from the unloading curve i n an oedometer, under s t a t i c conditions. But Finn (1981) showed that the rebound modulus measured i n the oedometer i s higher than the modulus computed from the slope of the curve shown i n F i g . 3.4. He used the E„ values computed from the slope of the curve shown i n F i g . 3.4 to v e r i f y the porewater pressure model. Finn (1981) maintained that the the s t r a i n hardening e f f e c t (equation. 3.14) should not be included i n the s t r e s s - s t r a i n r e l a t i o n s h i p when pred i c t i n g the behaviour of sands under undrained conditions. This i s because net volumetric s t r a i n s do not occur during undrained conditions. If drainage i s allowed to occur, then the e f f e c t s of s t r a i n hardening should be included. The porewater pressure model coupled with the s t r e s s - s t r a i n r e l a t i o n s h i p can be employed to predict l i q u e f a c t i o n strength curves. The strength curve plots of the c y c l i c shear stress r a t i o t / o v o versus the number of cycles to cause i n i t i a l l i q u e f a c t i o n , for normally consolidated (0CR=1) and over consolidated sands, obtained a n a l y t i c a l l y o b D •> M IO 0) V-1.00 -•> o» c >nfin 0.75 -u o |c 0 5 0 -0 ^ _ 3 to • 025 -CL • fc. O 0. 0 L T T T T T Sand type : Ottawa sand (C-109) a ' s 200 kN/m 2 . Relative density * 45% vo Legend Shear strain amplitudes Drained Undrained 0.056 % 0.100 % 0.210 % 0 3 0 0 % o 0 0 5 6 % • 0.100% * 0.200% A 0 3 1 4 % ) 0 2 0.4 0.6 0 8 1.0 12 1.4 Volumetric strain in percent, € y ( j % Fig. 3 .4 . Relationship Between Volumetric Strains and Porewater Pressures in Constant Strain Cyclic Simple Shear Tests, D_ = 45%. 1.6 18 ho 30 and exper imenta l ly , are shown i n F i g . 3 .5 . The experimental curve was obtained from undrained constant volume c y c l i c simple shear t e s t s . The comparison between the computed and measured l i q u e f a c t i o n s trength curves i s very good. This means that the assumptions made i n the formulat ion of the non- l inear h y s t e r e t i c e f f ec t ive s t r e s s - s t r a i n r e l a t i o n -ship are v a l i d . But, l i q u e f a c t i o n res i s tance curve p r e d i c t i o n i s an extreme case. F inn (1981) used th i s e f f e c t i v e s tress model to pred ic t porewater pressure development during undrained tests when subjected to constant c y c l i c shear s tress i n a simple shear apparatus. This constant c y c l i c shear s tress loading re su l t s i n an i r r e g u l a r s t r a i n h i s t o r y as the porewater pressures develop. Further the model parameters ( i = 1,4) used were obtained from constant c y c l i c shear s t r a i n t e s t s . T y p i c a l r e s u l t s obtained from two undrained tests are shown i n F i g . 3 .6 . The agreement between measured and computed porewater pressures i s remarkably good, i n d i c a t i n g that a l l the assumptions made i n the porewater pressure model and s t r e s s - s t r a i n r e l a t i o n s h i p are reasonable. 3.3 FIELD VERIFICATION OF EFFECTIVE STRESS RESPONSE ANALYSIS A unique opportunity to inves t iga te the c a p a b i l i t y of the one-dimensional e f f e c t i v e s tress response ana lys i s was provided r e c e n t l y when data became a v a i l a b l e on the dynamic response of an a r t i f i c i a l i s l a n d i n Tokyo Bay to the Mid-Chiba earthquake of 1980. Owi Is land No.1 i s an a r t i f i c i a l i s l a n d located on the west s ide of Tokyo Bay. It was constructed with mater ia l s dredged from the nearby sea. A tes t s i t e at the south end of the i s l a n d i s instrumented to record porewater pressures and ground acce lera t ions during earthquakes. 0.40 o - > b 0.30 o o w S 0.20 w •> O *•= 0.10 o Sand type ; Ottawa sand (C-109) orv'0 «200 kN/m 2 , Relative density » 4 5 - 4 7 % Analytical curve O , « , A , A Experimental data 3 6 10 30 60 100 Number of cycles for initial liquefaction, N L F i g . 3 . 5 . C y c l i c S t r e s s Rat io v s . N u m b e r of C y c l e s for I n i t i a l L i q u e f a c t i o n for V a r i o u s S a n d s . Sand type : Ottawa sand (C-109) o-v'0 « 200 kN/m 2 , Relative density = 45% T / c r ' =0.074 T/crV*a065 vo / / o Experimental curve Analytical curve 6 10 30 Number of cycles , N 60 100 200 Fig. 3.6. Predicted and Measured Porewater Pressures in Constant Stress Cyclic Simple Shear Tests, D r = 45%. 33 Porewater pressures are recorded by piezometers i n s t a l l e d at depths of 6m and 14m. The Mid-Chiba earthquake, with magnitude M = 6 . 1 , shook the Tokyo Bay area on September 25, 1980. F i n n , et a l . (1982) computed the response of Owi Is land No.1 to the Mid-Chiba earthquake using a one-dimensional e f f e c t i v e s tress response a n a l y s i s . The f i r s t 10 sees. of the recorded ground acce lera t ions are shown i n F i g . 3 .7 ( a ) . During the f i r s t 4 sees, very low acce lera t ions occurred . S i g n i f i c a n t acce le ra t ions developed between 4 and 6 sees . , and thereaf ter only low l e v e l e x c i t a t i o n was recorded. The ground acce lera t ions computed by F i n n , et a l . (1982) are shown i n F i g . 3 .7 (b ) . Except for some minor d i f ferences between 8 -10 sees, range the computed recording was very s i m i l a r to the recorded motions. The porewater pressures recorded at the 6m depth on Owi Is land No. 1 are shown i n F i g . 3 .8 (a ) . The recorded porewater pressure has two components: t rans ient and r e s i d u a l . The t rans ient porewater pressures are instantaneous response of porewater to changes i n t o t a l app l ied stresses and r e s i d u a l porewater pressures occur due to p l a s t i c volume changes. The one-dimensional response ana lys i s used by F i n n , et a l . (1982) computes the r e s i d u a l porewater pressure component and i s shown i n F i g . 3 . 8 (b ) . Comparison between recorded and computed porewater pressures i s very good. 3.4. POREWATER PRESSURE MODEL IN PRACTICE To apply the porewater pressure model i n dynamic e f f e c t i v e analyses , 7 constants must be known; four C. ( i = 1,4) constants to CJ a.: s (a) o.o i.o i i — 4.0 1.0 TIME 0.0 10.0 (b) i i — 4.0 (.0 TIME 0.0 —I 10.0 Fig. 3.7. Measured (a) and Computed (b) Ground Accelerations (Acc. in ft/sec2. Time in Sees). 8-;HWv(M*w r-TIME ' • • o.o 10.0 Fig. 3.8. Measured (a) and Computed (b) Porewater Pressure at a Depth of 6m. (Porewater Pressures in lb/ft2, Time in Sees). LO 35 compute incremental volumetric s t r a i n and 3 constants K r, m and n to represent rebound c h a r a c t e r i s t i c s . C y c l i c simple shear apparatus has to be used to obtain these constants. A number of laboratories s t i l l do not have simple shear apparatus to do these t e s t s . Over the years a procedure has evolved from a number of p r a c t i c a l experiences by which the d i r e c t measurement of these constants can be avoided (Finn, et a l . 1982). This i s done by modifying the model parameters such that i t w i l l match with the experimental l i q u e f a c t i o n strength curve and give the r i g h t rate of porewater pressure generation. The l i q u e f a c t i o n strength curve and the rate of porewater pressure development can be experimentally obtained by doing c y c l i c t r i a x i a l tests or c y c l i c simple shear tests on f i e l d samples. A study of a number of t r i a l analyses to predict the undrained behaviour of samples i n simple shear has revealed the following: a) The shape of the l i q u e f a c t i o n resistance curve i s s e n s i t i v e to the constants C^, e s p e c i a l l y C^. b) The v a r i a t i o n of K r s h i f t s the l i q u e f a c t i o n p o t e n t i a l up or down without changing the shape appreciably. In r e a l i t y the shape of the l i q u e f a c t i o n resistance curve for a number of sands i s s i m i l a r and the values of C^ given i n the l i t e r a t u r e give the shape of t y p i c a l l i q u e f a c t i o n p o t e n t i a l curves. In pr a c t i c e a t r i a l and error procedure i s adopted to get values for the model constants. The procedure i s outlined below: 36 1) Performing a number of analyses by varying K„, and select the value for K_ such that the computed l i q u e f a c t i o n resistance curve matches the experimental l i q u e f a c t i o n resistance curve. 2) For t h i s selected K_ value, c a l c u l a t e the development of porewater pressure with number of cycles and compare with the laboratory porewater pressure curve. 3) If these porewater pressure curves are not s i m i l a r , a l t e r and repeat the a n a l y s i s . It should be noted that i s the only parameter that i s used i n the c a l c u l a t i o n of incremental volumetric s t r a i n i n the f i r s t c y c l e . Therefore, estimates of C^, can be interpreted from the r e s i d u a l porewater pressure recorded i n the f i r s t c y c l e . This type of t r i a l and error procedure can be employed to obtain relevant model constants such that the corresponding porewater pressure development and l i q u e f a c t i o n resistance curves are s u f f i c i e n t l y close to the ones observed i n the laboratory. 3.5 DISCUSSION In the response analysis of h o r i z o n t a l l y layered deposits subjected to h o r i z o n t a l accelerations, a shear beam type of deformation pattern i s assumed i n the deposit. Therefore, only the s t r e s s - s t r a i n r e l a t i o n s h i p i n shear i s required. The tangent shear modulus i s used as the e l a s t i c parameter i n the incrementally e l a s t i c response analysis proposed by F i n n , et a l . (1977). To extend t h e i r model to two-dimensions, two e l a s t i c parameters are required. A d e t a i l e d 37 d e s c r i p t i o n of the extension of the one-dimensional s t r e s s - s t r a i n r e l a t i o n s h i p to two-dimensions i s discussed i n Chapter 4. It has been observed i n the laboratory that the presence of s t a t i c shear s tress a f fects the porewater pressure response of samples subjected to c y c l i c loading ( F i n n , et a l . 1978; V a i d , et a l . 1979). The porewater pressure model of F i n n , et a l . (1977) i s s t r i c t l y app l i cab le to one-dimensional depos i t s , where s t a t i c shear s tress i s ze ro . Therefore , i n extending t h e i r model to two-dimensions, the inf luence of s t a t i c shear must be accounted f o r . 38 CHAPTER 4 TWO DIMENSIONAL STATIC ANALYSIS OF SOIL STRUCTURES A s t a t i c response analysis to compute i n - s i t u e f f e c t i v e stresses i s necessary because the dynamic s o i l properties, such as strength and s t i f f n e s s , depend on i n - s i t u e f f e c t i v e streses. A number of s a t i s f a c t o r y incremental e l a s t i c methods, which model the construction sequence of the s o i l structures, are already a v a i l a b l e (Ozawa, et a l , 1973; Duncan, et a l , 1978, Byrne, et a l , 1982). The s t a t i c analysis presented i n t h i s theis i s based mainly on the methods proposed by these authors. The method proposed i n t h i s thesis uses a consistent set of material parameters i n both the s t a t i c and dynamic analyses; procedures also have been incorporated to apply c o r r e c t i o n forces during the a p p l i c a t i o n of the load increments. 4.1 STRESS-STRAIN RELATIONSHIP A number of s t r e s s - s t r a i n r e l a t i o n s have been proposed i n the computation of i n - s i t u s t a t i c stresses i n s o i l deposits. They can be divided broadly into l i n e a r , b i l i n e a r , e l a s t o - p l a s t i c , v i s c o - p l a s t i c and non-linear models. Some of these models are very complex and even for simple monotonic types of loading are expensive to use i n computational schemes. 39 Two e l a s t i c constants are required for any two-dimensional i s o t r o p i c , e l a s t i c or incrementally e l a s t i c a n a l y s i s . Tangent shear modulus G t and tangent bulk modulus B t were selected as the e l a s t i c constants. In s e l e c t i n g these e l a s t i c constants, considerations were also given to the f a c t that the s t r e s s - s t r a i n formulation proposed here has to be extended to model dynamic loading conditions. It w i l l be shown i n Chapter 5, that the s e l e c t i o n of these parameters g r e a t l y reduces the amount of computation time i n the dynamic a n a l y s i s . The s t r e s s - s t r a i n model proposed here, l i k e almost a l l other s t a t i c s t r e s s - s t r a i n models for s o i l s , can model only saturated s o i l s under f u l l y drained or undrained conditions, and dry s o i l s . The parameters selected to model the s t r e s s - s t r a i n behaviour should be based on test r e s u l t s which represent as c l o s e l y as possible the loading conditions that e x i s t i n the f i e l d . For example, i n the analysis of long term s t a b i l i t y of earth dams, one should chose parameters from drained test r e s u l t s . A d e s c r i p t i o n of the s t r e s s - s t r a i n model, and the s e l e c t i o n of relevant parameters for the model are discussed i n d e t a i l , i n t h i s chapter. 4.1.1 Reasons for Selecting G t and B_t In general, s t r a i n i n an i s o t r o p i c , homogeneous, l i n e a r e l a s t i c medium can be divided into two components: volumetric s t r a i n and d e v i a t o r i c s t r a i n . The volumetric s t r a i n i s related to mean normal stress through the bulk modulus. The d e v i a t o r i c s t r a i n i s related to d e v i a t o r i c stress through the shear modulus. These two independent material moduli can be e v a l u a t e d i n d e p e n d e n t l y by a p p l y i n g u n i f o r m changes i n AO corresponding stresses. Therefore, by sel e c t i n g tangent bulk modulus and tangent shear modulus as two independent e l a s t i c constants, better controls on stresses and stra i n s can be imposed. 4.1.2 Hyperbolic Shear S t r e s s - S t r a i n Relationship A number of researchers have used a hyperbolic s t r e s s - s t r a i n r e l a t i o n s h i p to predict the behaviour of a s o i l deposit (Konder, et a l . 1963). The hyperbolic r e l a t i o n s h i p between i and y i s given i n terms of Gmax a n d xmax a s> G y max ' T = Gmax M l max ( l + W m ) i n which, T,y = are the shear stress and shear s t r a i n G_ a x = tangent shear modulus as y+0 Tmax = u l t l m a t e shear strength 4.1.2.1 Estimation of Gmax Experimental data have shown that f o r sands and s i l t s under drained conditions, Gmax " f K » e» 0CR) <4-2) i n which, 41 o"m = current e f f e c t i v e mean normal stress e = void r a t i o OCR = over consolidation r a t i o here OCR i s defined as: _ Maximum past ma.jor p r i n c i p a l stress Current major p r i n c i p a l stress The following non-dimensional equation i s widely used f o r Gmax' Gmax * K G P a t ^ j " 2 (4-3) i n which, K Q = a non-dimensional constant for a given s o i l . Pfl = atmospheric pressure. The value of K Q depends mainly on void r a t i o or r e l a t i v e density of the s o i l , grain contact c h a r a c t e r i s t i c s such as angularity and roughness of the s o i l p a r t i c l e s etc., and also on previous loading h i s t o r y . An equation s i m i l a r to (4.3) has been proposed by Hardin, et a l . (1972) and Seed, et a l . (1970) for the computation of G M O „ for sandy s o i l s for dynamic analyses. The equation given by Hardin, et a l , (1972) includes the e f f e c t of previous stress h i s t o r y . They proposed, G = K _ P fm/P ) (OCR) (4.4) max G a ^ a ; v ' 42 i n which the exponent n depends on the p l a s t i c index of the s o i l (Hardin, et a l . 1 9 7 2 ) . Values of n are given i n Table 4 . 1 . Table 4 . 1 V a r i a t i o n of Exponent, r\ with P l a s t i c Index, PI PI% 0 0 20 0 . 1 8 40 0 . 3 0 60 0 . 4 1 80 0 . 4 8 >100 0 . 5 For normally consolidated non-plastic s o i l s under drained condition t y p i c a l values for KQ varies between 200 and 8 0 0 (Byrne, 1 9 7 9 ) . For clayey s o i l s under undrained conditions, G M A X can be related to the undrained strength S U through an equation, G M A X = K S U ( 4 . 5 ) where K i s a constant for a given c l a y . 43 4.1.2.2 Estimation of Tmax For s o i l s under drained conditions, t m a x > i s the maximum i shear stress that can be applied by keeping a x i a l stresses o x, and i a v at the i r respective values, where the x-axis i s taken h o r i z o n t a l and the y-axis v e r t i c a l for convenience. In p r i n c i p l e , this i s s i m i l a r to the e s t i m a t i o n of maximum d e v i a t o r i c s t r e s s ( o " j ) m Q v i n the a n a l y s i s presented by Kulhawy, et a l . (1969). They e s t i m a t e ( a ( j ) m a x assuming tha t the minor p r i n c i p l e s t r e s s remains constant. Let us consider a case where the i n i t i a l v e r t i c a l , h o r i z o n t a l i i and s h e a r s t r e s s e s a r e av, av, and -r r e s p e c t i v e l y ( F i g . y x xy 4.1). F i g . 4.2, shows the corresponding Mohr c i r c l e diagram and the Mohr envelope. The Mohr envelope i s defined by c' and 0 ' . The points L and M represent the i n i t i a l stress state. The a p p l i c a t i o n of shear stress w i l l increase the size of the Mohr c i r c l e and the largest Mohr c i r c l e i s the one that i s tangent to the Mohr envelope. i OA = OP radius of the largest c i r c l e keeping a y and i a x constant and i t i s given by, i i OP = r c f g x + q y ) i s i n 0 ' (4.6) "•tan0' 2 ' From t r i a n g l e ABO, A B - S a x plane = maximum shear s t r e s s e s on the h o r i z o n t a l F i g . 4.2. M o h r C i r c l e Diagram . 45 = /(OA 2 - OB 2) and so, f - * } 2 s i n 2 0 ' - H^} 2 ] 1/2 (4.7) x + max This equation reduces to the equation presented by Hardin and D r n e v i c h ( 1 9 7 2 ) i f a v and a a r e r e p l a c e d by K Q a o u vo and a, vo* The estimation of T max for a s o i l element under undrained conditions can be made based on standard f i e l d t e s t s , laboratory tests or may be based on estimation of i n s i t u e f f e c t i v e stress conditions. 4.1.2.3 Influence of Over-Consolidation Compaction i s generally used to obtain a c e r t a i n density i n dam construction. So some parts of a dam are over consolidated due to the compaction pressure. The e f f e c t of over consolidation on G__„ i s Illcl A. already shown i n the equation (4.4). Over consolidation has an influence on the value of x m a x a l s o . A t y p i c a l Mohr envelope for p l a s t i c s o i l s would look l i k e that i n F i g . 4.3. Dif f e r e n t c 1 and 0' values should be used depending on whether the s o i l i s i n NC state or OC state. 4.1.2.4 E f f e c t s of Unloading In geotechnical problems which involve excavation or reduction 46 co CO Q> CO u cd u CO O-C Region N«C Region CO CO 0) tH * J CO M Rt CO Normal Stress Fig. 4.3. Mohr Envelope. S lope G 2-(Loadin Slope G 4 (Unloading^ G x > G 2 Shear Strain Fig. 4*4 . Effects of Loading and Unloading, 47 i n applied load, some s o i l elements w i l l experience unloading. The modulus, G^, during unloading from a s t r a i n , y> a s shown i n F i g . 4.4, i s higher than the modulus, G 2 , corresponding to loading from the same s t r a i n l e v e l i f the mean normal stress remains a constant. Unloading and reloading can be modelled using the procedures presented i n Chapter 3. However, i f the s t r a i n ranges of i n t e r e s t are small, the differ e n c e between G 2 and i s not large. Under these circumstances changes i n modulus need not be modelled. 4.1.3 Tangent Bulk Modulus B t It i s assumed that the tangent bulk modulus (B t) i s e l a s t i c for any s o i l under drained conditions and i t i s a function of the current t mean normal stress a m only. Duncan, et a l . (1978), suggested that, B = K, P A}n (4*8) t b a LP ' where, = Bulk Modulus constant n* = exponent Typical values of bulk modulus constant vary between 300-1000 and the exponent varies between 0.3 and 0.6. depends mainly on r e l a t i v e density of the s o i l , s o i l type and previous loading h i s t o r y . In e l a s t i c (or incrementally e l a s t i c ) analysis i n i s o t r o p i c and homogeneous materials, a change i n shear stress with constant mean normal stress w i l l not r e s u l t i n any volume change. But, s o i l s under 48 constant mean-normal stress exhibit volume change when subjected to shearing s t r e s s . In the analysis presented herein, allowance has been made to include the volumetric s t r a i n s that occur due to shear stresses. A d e t a i l e d d e s c r i p t i o n of how th i s i s done i s given i n Section 4.4. It should be noted that the bulk modulus constant defined here includes only the e f f e c t of mean normal st r e s s . Therefore, for the estimation of K^, i s o t r o p i c consolidation tests performed i n t r i a x i a l t est equipment must be used. Duncan, et a l . (1978) described a procedure for determining from conventional t r i a x i a l test data i n which the mean-normal stress i s not held constant. Values of determined i n thi s manner must be considered approximate. 4.2 PHYSICAL MODELLING The domain of i n t e r e s t i s assumed to be an assembly of a f i n i t e number of elements, connected at the nodal points. The formulation of the f i n i t e element equations including the e f f e c t of porewater pressure ( C h r i s t i a n , et a l , 1970) i s presented i n Appendix I. The equations are, {P} = [ K T ] {A} + [K*] {U} (4.9) where, {p} = global column vector of incremental applied loads [ K F C ] = global tangent s t i f f n e s s matrix {A} = global column vector of incremental displacements [ K * ] = global s t i f f n e s s matrix defined i n the Appendix I 49 {u} = column vector of incremental porewater pressures i n the elements. The tangent element s t i f f n e s s matrix [ k t ] depends on two fac t o r s ; the tangent moduli and the shape functions adopted i n the f i n i t e element formulation. The shape functions give the v a r i a t i o n of displacements within an element i n terms of the nodal displacements. The simplest shape function, which assumes a l i n e a r v a r i a t i o n of displacements, gives constant s t r a i n within a tria n g u l a r element. But experience has shown that the r e s u l t s obtained from such elements do not predict stresses and s t r a i n s accurately. Therefore, q u a d r i l a t e r a l elements which have a l i n e a r s t r a i n v a r i a t i o n within an element are used. For s o i l structures such as dams, layered deposits etc., elements of a r b i t r a r y q u a d r i l a t e r a l shape are very appropriate because they are f a i r l y simple and can be used to model the geometry of these s o i l structures q u i t e a c c u r a t e l y . The element s t i f f n e s s m a t r i x , [ k t ] f o r an isoparametric q u a d r i l a t e r a l element i s given i n Appendix 1. 4 . 3 . SIMULATION OF CONSTRUCTION SEQUENCE Dams are constructed sequentially. Since the behaviour of dam m a t e r i a l s are n o n - l i n e a r and s t r e s s path dependent, a r e a l i s t i c computation of stresses and s t r a i n s requires that the construction sequence be modelled. An analysis based on single stage construction or gravity switch on, w i l l give f i n a l stresses and st r a i n s d i f f e r e n t from those calculated by following the construction sequence. 50 4.3.1 Method of Ana ly s i s F i g . 4.5 shows a schematic representat ion of the sequent ia l procedure involved i n dam c o n s t r u c t i o n . There may be p r e - e x i s t i n g elements on which subsequent layers w i l l be p laced . The cons t ruc t ion l i f t s are compacted u n t i l the required density i s obta ined . Thi s type of l ayer by layer cons t ruc t ion procedure i s c a r r i e d out u n t i l the required dimensions of the dam are obta ined. In m o d e l l i n g the c o n s t r u c t i o n sequence , the i n c r e m e n t a l s t re s se s , s t r a i n s and deformations are computed for every new layer added. This i s done by so lv ing equation (4.9) for the incremental loads caused by p lac ing a f resh l a y e r . The f i n a l s t res ses , s t r a i n s and displacements of the dam are simply the a lgebra ic sum of a l l the incremental values computed for a l l the l a y e r s . 4 .3.2 Incremental Porewater Pressure S t a t i c ana lys i s can be c a r r i e d out i n a t o t a l or e f f e c t i v e s tress mode or a combination of both . In the combination mode, some elements may be i n an e f f e c t i v e s tress mode and some may be i n a t o t a l s t ress mode. When the e f f e c t i v e s tress p r i n c i p l e i s used i n the f i n i t e element formula t ion , the porewater pressure term i s introduced as shown i n the r i g h t hand s ide of the equation (4 .9 ) . For the elements for which the t o t a l s tress mode i s assumed to be a p p l i c a b l e , the element s t i f fne s s matrix [ k t ] i s based on a t o t a l s t r e s s - s t r a i n r e l a t i o n s h i p . Furthermore, the porewater pressure , u Q , i n these elements should be set to zero . However, i f the e f f e c t i v e s tress 52 mode i s s e l e c t e d f o r some e l e m e n t s , t h e element s t i f f n e s s f o r t h e s e e l e m e n t s s h o u l d be based on an e f f e c t i v e s t r e s s - s t r a i n r e l a t i o n s h i p and a v a l u e f o r u Q i s r e q u i r e d . I t s h o u l d be remembered t h a t the r e a s o n f o r t h e s t a t i c a n a l y s i s i s t o e s t i m a t e the i n - s i t u s t a t i c c o n d i t i o n w h i c h i s a " l o n g t e r m " c o n d i t i o n . I n v i e w of t h i s , the g l o b a l p o r e w a t e r p r e s s u r e v e c t o r {Ti ] used i n the e q u a t i o n (4.9) s h o u l d c o r r e s p o n d t o t h e l o n g term v a l u e . E s t i m a t e s of u Q f o r t h e e lements can be made u s i n g a number of methods s u c h as h y d r a u l i c model t e s t s , e l e c t r i c a l a n a l o g y e t c . Measured p o r e w a t e r p r e s s u r e s i n t h e f i e l d a l s o may be u s e d . The m a t r i x {IT} , f o r m u l a t e d u s i n g element p o r e w a t e r p r e s s u r e s u Q can now be used i n e q u a t i o n (4.9) t o compute s t r e s s e s and d i s p l a c e m e n t s . 4.3.3 C o m p u t a t i o n o f I n c r e m e n t a l S t r e s s e s and S t r a i n s S t r a i n s g i v e n by the f i n i t e element a n a l y s i s a r e a measure o f changes i n shape of the e l e m e n t s from some r e f e r e n c e s t a t e . I t i s assumed t h a t the c o n d i t i o n of newly p l a c e d e lements a f t e r t h e y have s e t t l e d under t h e i r own w e i g h t i s the r e f e r e n c e s t a t e (Ozawa, et a l . 1 9 7 3 ) . The t o t a l s t r a i n s a r e o b t a i n e d by a d d i n g i n c r e m e n t a l s t r a i n s c aused by t h e c o n s t r u c t i o n l a y e r s about t h i s r e f e r e n c e s t a t e . An i n c r e m e n t a l e l a s t i c a n a l y s i s can be c a r r i e d out i n a number o f ways ( D e s a i and A b e l , 1 9 7 2 ) . The a p p r o a c h a d o p t e d h e r e i s shown s c h e m a t i c a l l y i n F i g . 4.6. E s t i m a t e s of the i n c r e m e n t s i n s t r e s s e s and s t r a i n s due t o a l o a d i n c r e m e n t a r e d e t e r m i n e d u s i n g m o d u l i v a l u e s c o r r e s p o n d i n g t o s t r e s s - s t r a i n l e v e l b e f o r e t h e l o a d i n c r e m e n t was a p p l i e d . New m o d u l i c o r r e s p o n d i n g t o the average of t h e s t r a i n s b e f o r e I n i t i a l Estimates or Computed Values o^ 't t 1 — r Layer Loading (incremental) Estimated New Moduli Stresses and Strains Calculate K G ' T , B[ from W - W + M — 1 > Average Strains M = K I + M Moduli and Stresses V. i I n i t i a l Estimates or computed Values with modified G'T and B£ Layer Loading (incremental) F i n a l Stresses and Strains {e} = {e} + {AE} {a} = {a0} + {Aa} Next Layer Loading Fig. 4.6. Schematic Diagram Showing the Incremental Analytical Procedure 54 the load increment i s applied and the s t r a i n s computed a f t e r the increment are used to compute more correct incremental stresses and s t r a i n s for the same load increment. These incremental stresses and s t r a i n s are added to i n i t i a l stresses and s t r a i n s to obtain the i n i t i a l condition for the next load increment. Re c a l l from section 4.1.2, that the r e l a t i o n s h i p between shear stress and shear s t r a i n i s assumed to be hyperbolic. Therefore, the shear stresses are computed using the shear s t r a i n s at the end of the load increment. In doing t h i s , as pointed out by Desai and Abel (1972), equilibrium i s not n e c e s s a r i l y s a t i s f i e d . Under these circumstances equilibrium c o r r e c t i o n forces may be applied to s a t i s f y equilibrium condition. In the method adopted here, correction forces that correspond to the changes i n the shear stresses computed using the procedure outlined i n Appendix I, are applied to the next load increment. Before placing a fresh layer the stress condition i n previously placed elements (pre-existing) are known. Therefore, i n the i n i t i a l a nalysis for the load increment, moduli for the pre-existing elements are known. However, for the f r e s h l y placed elements, moduli must be based on estimated stresses i n these elements. Ozawa, et a l . (1973) suggested that the stresses can be estimated using the equations, and y ' s T = 0.5 v d s i n <* lxy is o (4.10) i n which d i s the depth of center of gravity of the element from the top surface, <=0 i s the slope of the top surface and y s i s the unit weight 55 of s o i l . For t o t a l l y submerged elements Ys should be replaced by the submerged unit weight y'« 4.4. SHEAR-VOLUME COUPLING The tangent bulk modulus Bfc defined i n section 4.1.3, re l a t e s an increment i n volumetric s t r a i n , A e ^ , to an increment i n e f f e c t i v e i mean normal s t r e s s , Aa m« But i n s o i l s volumetric s t r a i n s occur also due to changes i n shear stresses. This a d d i t i o n a l volumetric s t r a i n must be accounted for i n any r e a l i s t i c modelling of s o i l behaviour. The c h a r a c t e r i s t i c drained behaviour of i n i t i a l l y loose and dense sand samples i n a simple shear apparatus i s shown i n F i g . 4.7. I n i t i a l l y for small shear st r a i n s y, both the loose and the dense samples undergo volume reduction. But l a t e r , over a considerable range of s t r a i n , they e xhibit volume expansion ( d i l a t i o n ) . For both samples i n the d i l a t i o n range, the v a r i a t i o n of volumetric s t r a i n e v vs Y ^ s l i n e a r i n i t i a l l y and then e v approach fixed values at very high s t r a i n l e v e l s . The region of i n t e r e s t i n the e y vs Y p l o t i n t y p i c a l geotechnical problems would be the i n i t i a l compaction region and the l i n e a r d i l a t i o n region. The rate of volume change i n the l i n e a r d i l a t i o n region, i s larger for the dense sand than for the loose sand. Hansen, (1958) suggested using d i l a t i o n angle v Q to characterize the d i l a t i o n r a t e . He defined v as, sinv = o dy~ = t a n P o (4.11) 5 6 F i g . 4.7. (b). Typ ica l Plots for £ Vs T for Dense and Loose Sands . 57 where 8 Q i s the slope angle. The negative sign i s introduced since compressive volumetric s t r a i n i s considered to be p o s i t i v e . For a given type of sand the angle v Q was found to be a function of the r e l a t i v e density and confining pressure. The d i l a t i o n angle v Q increases with the r e l a t i v e density of the s o i l . This was c l e a r l y shown by Vaid, et a l . (1981). They performed drained simple shear tests with constant v e r t i c a l stress a y 0 = 200 kPa, ( F i g . 4.8) on Ottawa sand, (C-109), at various r e l a t i v e d e n s i t i e s . The d i l a t i o n angle, which i s the slope of the l i n e a r d i l a t i o n portion of the plo t e v and y i s found to increase with the r e l a t i v e density of the s o i l . The d i l a t i o n i s also a function of the mean normal stress l e v e l . This was shown by Lee (1965) who performed drained t r i a x i a l tests on dense Sacramento River sand samples of constant D r = 100%. F i g . 4.9 shows a series of tests with consolidation pressures varying from 0.1 MPa to 13.7 MPa. Several important features of the test data can be noted i n F i g . 4.9. F i r s t l y , dense samples at high consolidation pressures behave l i k e loose samples; secondly, f a i l u r e i n terms of maximum p r i n c i p a l stress r a t i o occurs at increasing s t r a i n l e v e l s as the consolidation pressure increases; and t h i r d l y , the d i l a t i o n angle decreases and becomes negative (compaction) with increasing consolidation pressure. F i g . 4.10 shows the v a r i a t i o n of the d i l a t i o n angle v Q with mean normal stress for a number of sands which were at an i n i t i a l r e l a t i v e density of 80 percent. It i s i n t e r e s t i n g to observe that the v a r i a t i o n of v Q l i e s within a narrow band (Robertson, 1982) and also the v a r i a t i o n i s li n e a r with logarithm of mean normal s t r e s s . Based on the experimental data presented above, the following approximation for a n a l y t i c a l purposes can be made for medium dense and 58 Fig. 4 . 8 . Stress-Strain Behaviour of Ottawa Sand in Drained Simple Shear. (After Vaid et. al., 1981) F i g . 4 .9 T y p i c a l Drained T r i a x i a l T e s t Results on Dense Sacramento Ri v e r Sand. (a) Pr i n c i p a l Stress Ratio V s A x i a l Strain (b) Volumetric Strain Vs A x i a l Strain ( A f t e r Lee, 1965) ~ 20 V « T J Ul -I o z < z o 10 ^'/(degree*) Dj.- 802 32 3 2 3 5 •37 3 3 3 5 x V Chattahoochee Sand Mol Sond Monterey Sond Glacial Sand SATAF Leigh ton Buzzard Sond Vesie and Clough 1968 OeBeer 1965 ViKet ond Mitchell 1981 Hirshfteld ond Pbuloe 1963 Boldi et ol. 1981 Cole 1967 0.1 0.5 I 5 10 50 100 .1 500 1000 MEAN NORMAL STRESS, crm Kg/cm3Fig. 4.10. Variation of Dilation Angle, with Mean Normal Stress for Various Sands. (After Robertson, 1982). ON o 61 dense sands: the volumetric s t r a i n due to shear stresses ( e v ) i s n e g l i g i b l y small (Varadarajan, et a l . 1980, Byrne, et a l . 1982) at low shear s t r a i n l e v e l s and above t h i s l e v e l i t varies l i n e a r l y with y. This means that the p l o t of e y vs y can be i d e a l i s e d as i n F i g . 4.11, where yQ i s the shear s t r a i n above which the volumetric s t r a i n due to shear stress i s important. It should be noted that the value of v Q should be modified for the changes i n mean normal stress according to some v a r i a t i o n such as shown i n F i g . 4.10. 4.4.1 A n a l y t i c a l Formulation There are a number of ways of modelling shear-volume coupling. One i s to modify the e l a s t i c i t y matrix I) (Appendix I) such that A e x and Ae y depend also on shear stress increment ( V e r r u i j t , 1977). But t h i s type of approach w i l l give r i s e to an unsymmetrical s t i f f n e s s matrix, which unduly complicates the computations. A simpler way i s to keep the I) matrix as i t i s and to incorporate the volume change the same way as the temperature v a r i a t i o n s are analysed i n s t r u c t u r a l mechanics (Zienkiewicz, et a l . 1967; Byrne, 1979, Byrne, et a l , 1982). This i s accomplished i n the following manner. a) The incremental shear st r a i n s i n a l l elements are computed for the increment i n load, neglecting shear-volume coupling. b) v Q can be estimated for the new mean normal s t r e s s , using F i g . 4.10, and then Ae y i s computed using equation 4.11. c) The v o l u m e t r i c s t r a i n A e v then i s s p l i t i n t o A e x and ^ f t . Here i t i s assumed that, 63 A e d = Aev = 0*5 Aev* L e t u s d e f i n e a s t r a i n x y v T v e c t o r Ae such that i t s components are the estimated "0 d i l a t i o n a l s t r a i n s . Then A e T i s g i v e n by; JAe^, ~ 0 Ae y, 0}. d) The nodal forces corresponding to this s t r a i n vector Ae^ can be computed as, / / / I*?. A lo d v (4.12) V (see Appendix I) Now these forces can be added to the applied incremental load i n a) and new s t r a i n s and stresses can be computed. For computing incremental stresses, the following equation should be used, Aa = D_(Ae - A_eQ) (4.13) where, A£ = s t r a i n vector computed for the modified applied load. e) Now steps b) ->• d) can be ca r r i e d out u n t i l convergence occurs i n stress and s t r a i n increments under the applied incremental load. 4.5 INTERFACE REPRESENTATION I t may be necessary to allow r e l a t i v e displacement to occur at the interface between two f i n i t e elements to model s l i p surfaces i n the 64 f i e l d . S l i p elements, which are sometimes referred to as elements of zero thickness, can be used to model t h i s r e l a t i v e displacement. S l i p elements can be assumed to be placed along the boundaries between the two-dimensional elements representing s o i l and s t r u c t u r a l elements or wherever i t i s anticipated that r e l a t i v e movements or separation between elements may occur. The s l i p i s assumed to occur only along t h i s d i r e c t i o n and thi s occurs when the shearing forces i n the s l i p element exceeds the shear strength at the i n t e r f a c e . Goodman, et a l . (1968) have developed a two-dimensional s l i p element with eight degrees of freedom to represent j o i n t and f a u l t behaviour i n rock mechanics problems. F i g . 4.12 shows a s l i p element with nodes I,J,K and L, i n global and element axes. The forces at any point i n a s l i p element are the shear force f„ and the normal force f„ expressed per unit area of the element. The force-displacement r e l a t i o n s h i p i s assumed to be: (4.14) n n n where K g, K n = j o i n t s t i f f n e s s per unit length i n shear and normal d i r e c t i o n s r e s p e c t i v e l y . wg, wn = Shear and normal displacement at the point of i n t e r e s t . The d e f i n i t i o n of unit j o i n t s t i f f n e s s needs c l a r i f i c a t i o n . Imagine a d i r e c t shear test being performed along an i n t e r f a c e element of unit thickness. At f i r s t , when a normal force i s applied, the element 65 X Shear/Normal Displacementa v^ ,*^ Fig. 4.13. Plot of Typical Shear/Normal Stress vs Shear/Normal Displacement. 66 shortens as the a s p e r i t i e s i n the j o i n t deform. A t y p i c a l plot of normal deformation at the j o i n t and the force applied per unit length i s shown i n F i g . 4.13. For a n a l y t i c a l purposes the r e l a t i o n s h i p can be approximated to a s t r a i g h t l i n e and the slope i s given by 1^. Similar tests can be performed i n the tangential d i r e c t i o n and a plot between f„ and w„ can be obtained. The slope of th i s curve w i l l give K„. Using the equation 4.14 and also assuming a l i n e a r v a r i a t i o n of displacement within the j o i n t element, a s t i f f n e s s matrix K _ can be — - o i l obtained i n l o c a l or element co-ordinates. This s t i f f n e s s matrix r e l a t e s the nodal forces and the nodal displacements. The displacement vector here i s simply, u T = {u I f V-p U j , v j , u K , v K , u L , v L} It has been shown i n Appendix II that K s n i s given by, 0 K s 0 ~ K s 0 -2KS 0 2K n 0 0 " K n 0 -2K, 2 K s 0 " 2 K s 0 " K s 0 2 K n 0 - 2 * n 0 ~ K n 2 K s 0 K s 0 sym 2 K n 0 2 K s Kn 0 2 K n n (4.15) 6 7 where, L i s the length of the element. To get the s t i f f n e s s matrix i n global co-ordinates a simple transformation i s used, where T= transformation matrix containing d i r e c t i o n cosines and i s given by; (4.16b) i n which, cosec sin<* ., . r N PQ = r o o-i (4.16c) '--sinoc cos« J o o and a = angle of i n c l i n a t i o n of the s l i p element with the h o r i z o n t a l . Even though t h i s type of formulation does not include r o t a t i o n of the element d i r e c t l y , this i s taken into account since a l l 8 d.o.f. have been considered. The displacement f i e l d v a r i a t i o n assumed within the s l i p element i s consistent with the displacement f i e l d i n an isoparametric q u a d r i l a t e r a l f i n i t e element. Furthermore, both elements have the same degrees of freedom. Thus e s t a b l i s h i n g a global s t i f f n e s s matrix can be evaluated simply t r e a t i n g the s l i p element l i k e any other element. 68 4.5.1 Method of Analysis of S l i p Elements I t was assumed tha t the t a n g e n t i a l s t r e s s - d i s p l a c e m e n t r e l a t i o n s h i p ( f g vs w g) i n the s l i p element i s e l a s t i c - p e r f e c t l y p l a s t i c and the p l a s t i c region where s l i p occurs i s defined by a simple Mohr-Coulomb type of y i e l d c r i t e r i o n . For incrementally e l a s t i c a n a l y s i s , the values of K„ and K_ should be kept constant u n t i l y i e l d occurs. A f t e r the y i e l d , K„ i s set to a very small value. This small value can be viewed as the r e s i d u a l shear s t i f f n e s s . The value of K^ i s kept at i t s o r i g i n a l value. But, i f the normal force on the element i s negative, meaning that the separation of the j o i n t occurs, then the v a l u e s K g and IC^ should be set to a small value. To investigate whether y i e l d i n g i s possible or not, the shear and normal stress i n the s l i p element should be determined. Since a l i n e a r v a r i a t i o n of displacement i s assumed within an element, the shear and normal stresses vary from point to point within the element. The average incremental values of shear stress A f g and normal stress A f n for a load increment can be estimated from equation (4.14) as, A f g = K s (Awg) and A f n = 1 ^ (Awn) (4.17) where Aw. and Aw are i n c r e m e n t a l average shear and normal 5 XI displacement i n the element r e s p e c t i v e l y . Now expressions for Aw and Awn are, 69 A w s = < A u t o p > a v e " ( A u b o t t o m > a v e = (Au K+Au L)/2 - (Au I+Au J)/2 (4.18) and Awn = ( A v t Q p ) a v e - ( A v b o t t o m ) a v e = (Av K+Av L)/2 - (A V l+Avj)/2 (4.19) From equations 4.17 to 4.19, Af„ and Af_ can now be w r i t t e n as, A f s = K s [(Au K+Au L)/2 - (Au I+Au J)/2] and Af~n = 1^ [(Av K+Av L)/2 - (Av I+Av J)/2] (4.20) Tot a l shear and normal stresses f g and f can be obtained by adding up the incremental stresses for a l l the load steps. Mohr-Coulomb f a i l u r e c r i t e r i o n g i v e s the shear s t r e n g t h f m a x i n the element at any time as, fmax = c s + f n t a n < (4.21) where c and 0 are the cohesion and f r i c t i o n angle required to define s s the f a i l u r e c r i t e r i o n . If f m a x i s greater than the absolute value of fg then the s l i p element nodes could separate and t h i s i s modelled as mentioned above, by reducing the K g to a small r e s i d u a l value. The separation of a s l i p element i s indicated by negative f n . I t should be noted that a l l c a l c u l a t i o n s for the computation of f g and f are performed i n the l o c a l axes. Since the displacements from f i n i t e element analysis are given with respect to global axes, they 70 must be transformed to get displacements i n the l o c a l axes by using the inverse of the transformation matrix T. 4.5.2 Factors that Influence J o i n t Parameters In the a n a l y t i c a l model presented above, three d i s t i n c t j o i n t parameters were introduced. (1) K„, the unit s t i f f n e s s along the element. (2) K R, the unit s t i f f n e s s across the element. (3) shear strength, f m a x defined by c g and 0 g. These parameters model the behaviour of s l i p elements adequately. The v a l u e s of K„, K_ and fmav w i l l depend on (1) the S I I 111 cL 2v surface roughness of the adjacent elements (2) shape and c h a r a c t e r i s t i c s of the a s p e r i t i e s , and (3) contact area r a t i o between the j o i n t w a l l s . D e t a i l s on how these parameters can be obtained i n the laboratory are given by Goodman, et a l . (1968). 4.6 SELECTION OF SOIL PARAMETERS The process of obtaining representative values for s o i l properties i s probably one of the d i f f i c u l t tasks i n stress a n a l y s i s . It should be emphasized that i n d i v i d u a l estimation of s o i l properties i s not important. But the s t r e s s - s t r a i n v a r i a t i o n given by the selected s o i l parameters should give the best f i t to the observed laboratory behaviour of s o i l samples i n the stress (or st r a i n ) range of i n t e r e s t . 71 4.6.1 Obtaining Shear Stress - S t r a i n Parameters Aft e r deciding what drainage condition i s l i k e l y to occur i n the f i e l d , a laboratory test can be performed simulating the drainage conditions. For example to obtain parameters for an e f f e c t i v e s t r e s s -s t r a i n r e l a t i o n s h i p , a serie s of drained simple shear tests can be performed to obtain plots % vs y for various constant mean normal stress l e v e l s . Simple shear tests are i d e a l since the mean normal stress during the test remains reasonably constant. A simple, t r i a l and error method can be employed to o b t a i n v a l u e s f o r G m o v and that f i t these curves i n the stress or s t r a i n ranges of i n t e r e s t . Then knowing the stress l e v e l s corresponding to a te s t , the best estimates f o r the parameters K Q , C' and 0 ' can be obtained e a s i l y . The e f f e c t i v e stresses i n simple shear at the beginning of the dr a i n e d l o a d i n g c o n d i t i o n s can be assumed to be o" v o and K Q t t o " v o . Then the mean n o r m a l s t r e s s a m i s g i v e n by (1+2K D) a V Q/3 and t h i s i s i n general assumed to remain constant during the t e s t . But i f conventional t r i a x i a l tests are performed on the samples, then Ojjj varies as the a x i a l load v a r i e s , and therefore, i t i s not easy to obtain these parameters. If t r i a x i a l test data only are a v a i l a b l e then again the above procedure can be c a r r i e d out by considering the shear i stress and the shear s t r a i n on the f a i l u r e plane. However, since o"m increases during shearing the shear s t r e s s - s t r a i n curve obtained by th i s i procedure may be interpreted for a constant average am- The average i i o"m can be assumed to be the mean of o m , c o r r e s p o n d i n g to the beginning and the end of the t e s t . 7 2 4.6.2 Obtaining Bulk Modulus Parameters The tangent bulk modulus B t was assumed to be (equation 4.8) given by, da' a' n B^ " ~A = K, P {—) (4.22) t de b a KV ' v ' vm a Integrating t h i s equation one gets, f K p (4.23) 1-n vm b a ' taking logarithms both sides, (l-n)log(o-' m) = log {K bP^ _ n(l-n)} + log (e f f l) (4.24) i . e . , l o g ( e v m ) = ( l - n ) l o g ( a ' m ) - log {K bp( 1 _ nh-n)} The slope of the p l o t l o g ( e v m ) vs l o g ( a m ) o b t a i n e d from drained i s o t r o p i c t r i a x i a l test r e s u l t s w i l l give (1-n), and from t h i s n can be determined. Using the value of n, and the inter c e p t on log (e.yjjj) axis, can be c a l c u l a t e d . Now knowing the bulk modulus parameters, and n, the tangent bulk modulus can be computed at any given mean normal s t r e s s . 73 CHAPTER 5 TWO-DIMENSIONAL DYNAMIC ANALYSIS 5.1 FORMULATION OF THE PROBLEM The general dynamic equilibrium equations for a l i n e a r system at any time are given by (Clough and Penzien; 1975), [M] {X} + [C] {X} + [K] {X} = {P} (5.1) i n which {x}, {£}> { x} = column v e c t o r s whose components X^, X^, and X^ g i v e the r e l a t i v e a c c e l e r a t i o n , v e l o c i t y and displacement with respect to the base motion re s p e c t i v e l y of a node, [M] = mass matrix [c] = damping matrix [K] = s t i f f n e s s matrix {p} = i n e r t i a force vector, which i s defined as, -[MJ{l}x b where, {i} i s a vector with a l l components unity and X^ i s the base a c c e l e r a t i o n . In two-dimensional problems the base a c c e l e r a t i o n may have two components: h o r i z o n t a l and v e r t i c a l . If the i t ' 1 equation i n (5.1) i s w r i t t e n f o r the h o r i z o n t a l d i r e c t i o n t h e n X^ i s the base acceleration i n the ho r i z o n t a l d i r e c t i o n , and i f the equation i s for the 74 v e r t i c a l d i r e c t i o n then X^ i s the base ac c e l e r a t i o n i n the v e r t i c a l d i r e c t i o n . 5.1.1 Incremental Eq u i l i b r i u m Equations f o r Non-Linear Systems In the analysis of non-linear systems, the material properties change with time. An incrementally e l a s t i c approach has been adopted to model the non-linear material behaviour. Incremental dynamic equilibrium equations for any time i n t e r v a l , At, can be obtained by replacing the variables i n equation (5.1) by incremental values. This leads to, [M] {AX} + [C] {AX} + [K] {AX} = {AP} (5.2) i n which, [M], [C] and [K] are the mass, damping and s t i f f n e s s matrices relevant to the time i n t e r v a l f o r which the above equations are written. It i s always assumed that the mass matrix i s constant. The mass matrix can be obtained by two ways: lumped mass method and consistent mass method. In the lumped mass method, the mass matrix i s obtained by lumping the mass of a f i n i t e element equally at i t s nodes. The consistent mass matrix, i s obtained using the same i n t e r p o l a t i o n functions used i n the f i n i t e element formulation. The lumped mass matrix i s very simple to compute and i t has only diagonal terms, whereas the consistent mass matrix i s somewhat harder to compute and has off-diagonal terms. Even though the consistent mass method i s more accurate, the presence of off-diagonal terms greatly increases the computation time required to solve the dynamic 75 equilibrium equations. For the accuracy l e v e l required i n t y p i c a l geotechnical problems, the lumped mass method i s considered appropriate. In general, the damping matrix [ c ] and the s t i f f n e s s matrix [K] which are introduced i n the incremental equilibrium equation (5.2) depend on the d i s t r i b u t i o n of v e l o c i t y and displacement i n the structure. Appropriate values for the time i n t e r v a l between any time t and t+At can be determined only by an i t e r a t i o n procedure, because the v e l o c i t y and displacement at the end of a time increment depend on the i n i t i a l s t i f f n e s s and damping values. This type of i t e r a t i v e s o l u t i o n scheme for every time step i s very expensive. Therefore, i n practice tangent damping and tangent s t i f f n e s s matrices which correspond to time t are used with appropriate corrections to the r e s u l t s . It w i l l be explained l a t e r how correction forces can be introduced into the so l u t i o n scheme. The s t r e s s - s t r a i n law used to determine the tangent s t i f f n e s s matrix at time t, [K F C] i s described i n Section 5.2. The [ c ] matrix w i l l be assumed to be a constant throughout the dynamic analysis and the procedure to evaluate t h i s i s presented i n Section 5.3. With, [K] = [ K t ] t and (5.3) [C] = [C] (5.4) the dynamic incremental equilibrium equations can be re-written as, 76 [M] {AX} + [C] {AX} + [ K t ] t {AX} = {AP} (5.5) When computing response for a random loading h i s t o r y , equations (5.5) have to be solved for every time step. During the procedure [ K t ] t and {AP} are updated. The step by step i n t e g r a t i o n procedure proposed by Newmark (1959) or Wilson's 0-method (Wilson, et a l . 1973) have been adopted to integrate the equations. These procedures are described i n Appendix I I I . 5.1.2. Correction Forces The n u m e r i c a l i n t e g r a t i o n procedure i s based on three s i g n i f i c a n t assumptions: (1) the v a r i a t i o n of acc e l e r a t i o n within a time step i s assumed to vary i n some known fashion e.g. l i n e a r or constant (2) the damping and s t i f f n e s s properties remain constant during a time step and (3) the response at time t + At can be evaluated from the known response at time t. In general neither of these assumptions i s e n t i r e l y c o r r e c t , even though the errors may be small when the time step i s short. If errors accumulate from step to step gross errors and even s o l u t i o n i n s t a b i l i t i e s may occur. These problems can be avoided by imposing the condition of global equilibrium at each step of the a n a l y s i s . The global equilibrium equations at time t i n terms of a l l force components are, { f l i t + { f D l t + { fslt = { P l t ( 5 . 6 ) i n w h i c h {^i}t» I^D^t' a n d representing i n e r t i a , damping, and mass system at any time t and {p}t t . Since the mass matrix and constant during dynamic a n a l y s i s , by, { f g } t a r e the column v e c t o r s spring forces acting on the d i s c r e t e i s the i n e r t i a force vector at time the damping matrix were assumed to be { f j } t and {felt a r e simply given { f l i t = [M] ( X } t and ( 5 . 7 ) { f D } t = [C] {X} t ( 5 . 8 ) The spring forces {^slt» c a n ^ e computed by expressing element stresses i n terms of nodal forces, applied at the nodes of the elements. The nodal forces that correspond to the dynamic stresses i n an element {f^elem* ^ r o m Appendix I are, i n which _B i s the matrix that r e l a t e s s t r a i n to nodal displacements. In t h i s manner nodal forces for a l l the f i n i t e elements can be computed and 78 the vector sum of a l l these nodal forces w i l l give the global vector { f S } f If the solutions obtained at time t are accurate then the ri g h t and l e f t hand sides of the equation (5.6) w i l l be i d e n t i c a l . But, i n general i t w i l l not be so. The corr e c t i o n force vector { p c o r r} i s given by, ( Pcorr} = l p } t " { f l l t " { f D } t " { f s k <5-10) From equations (5.6) to (5.9), equation (5.10) can be re-written as, { Pcorr> - ^ t " M W t - [C] {*} - J /// B C O j dv a l l elements (5.11) To impose equilibrium, the corr e c t i o n force vector { P C O r r ^ c a n ^ e added to the incremental equilibrium equations formulated at time t . This i s accomplished by adding { p c o r r} to the rig h t hand side of the equation (5.5). 5.2 DYNAMIC STRESS-STRAIN RELATIONSHIP In the proposed incrementally e l a s t i c analysis i n the time domain, the tangent shear modulus G t and tangent bulk modulus B t were selected as the two " e l a s t i c parameters". Some reasons for s e l e c t i n g G t and B t for s t a t i c analysis were presented i n Chapter 4 . There i s a further very important reason for adopting these parameters i n the dynamic 79 a n a l y s i s . In dynamic an a l y s i s , the moduli have to be changed for every time step. This means that the element s t i f f n e s s matrix has to be re-formulated each time. This time consuming procedure can be s i m p l i f i e d somewhat i f G t and B t are used as the " e l a s t i c " constants. The e l a s t i c i t y matrix D_ (Appendix I) under plane s t r a i n conditions, i s given by, D = B t + 3 G t B t " 3 G t 0 = B. 1 0 1 0 B t - 3 G t B t + 3 G t 0 0 + G. 0 G 4. 3 -2 3 0 -2 3 4. 3 0 (5.12) (5.13) B t % + Gt*2 (5.14) where and Q 2 are two constant matrices. From the formulation of s t i f f n e s s matrix presented i n Appendix I, the element s t i f f n e s s matrix i s , k = J / / B D B dv V (5.15) Now s u b s t i t u t i n g for D from (5.14) the equation (5.15) can be rewritten as, 80 k. = B J / J ^ £ B_ dv + G /// B 1 L B dv (5.16) V V i . e . kfc - B t Rx + G t R2 ( 5 > 1 7 ) where R != / / / I * £L B dv V x (5.18) (5.19) The constant matricies R^ and R^ have to be computed only once. The current l c t matrix may now be obtained by multi p l y i n g the constant matrices R, and R 9 by the current values of B t and G t. 5.2.1 Volume Change Behaviour With regards to volume change behaviour during dynamic loading, the s o i l deposits can be divided into two basic groups. The f i r s t group includes s o i l s which can undergo volume changes under the load increments induced by the base e x c i t a t i o n and the second group includes s o i l s which cannot. Saturated gravels and dry deposits belong to the f i r s t group of s o i l s . R e c a l l equation (4.8) which re l a t e s B,. to a', 81 nun 3 J (5.20) This equation may be used to compute B t. This means that B t has to be modified for every time step. However, i t i s known that the changes i n mean normal stresses i n the s o i l elements, due to seismic e x c i t a t i o n i s small and furthermore the volume change behaviour does not influence the response of s o i l structures s i g n i f i c a n t l y . Therefore, f o r s i m p l i c i t y , B t may be kept constant throughout the dynamic a n a l y s i s . An average B t for elements can be evaluated based on i n - s i t u mean normal stresses using equation (5.20). This i s because the load pulses during seismic loading induce stresses such that t h e i r mean values are i n i t i a l i n - s i t u i stresses. It should be noted here that even i f the changes i n a m are considered, f o r t y p i c a l values of n the changes i n B t w i l l be small. Laboratory r e s u l t s have revealed that, as p l a s t i c volume changes occur during c y c l i c loading, the s o i l samples harden leading to higher bulk modulus. This increase i n bulk modulus due to s t r a i n hardening e f f e c t can be modelled the same way as the increase i n maximum shear modulus was modelled by Finn, et a l . (1977). This was accomplished by introducing hardening constants (equation 3.14). Loose saturated, sandy s o i l s and saturated clays belong to the second group of s o i l s . In saturated s o i l s volume change can occur only by porewater drainage. Within the short duration of t y p i c a l seismic e x c i t a t i o n , the occurrence of a p p r e c i a b l e amount of drainage i s questionable i n s o i l s of low permeability. In view of t h i s , the dynamic analysis proposed here assumes that no drainage occurs during the dynamic 82 loading. In saturated gravels appreciable r e s i d u a l porewater pressure does not develop because of i t s high permeability. For the second type of s o i l s , to simulate the condition of no volume change, the bulk modulus i s set to a very high value during dynamic loading conditions. 5.2.2 Dynamic Shear S t r e s s - S t r a i n Behaviour In the formulation of a complete dynamic e f f e c t i v e s t r e s s - s t r a i n r e l a t i o n s h i p the following basic aspects should be considered; 1) s o i l behaviour under i n i t i a l loading, unloading and reloading phases. 2) r e s i d u a l porewater pressure generation and i t s e f f e c t s . 5.2.2.1 Skeleton Curve f o r Dynamic Loading Under dynamic loading conditions the r e l a t i o n s h i p between dynamic (or c y c l i c ) shear stress, % c > and dynamic shear s t r a i n , y c» I s assumed to be hyperbolic. The hyperbolic r e l a t i o n s h i p (equation 4.1) i s defined by G m a x and - c m a x . Seed et. a l , (1970) proposed that maximum shear modulus, ^max ^ o r s a n d y s o i l s i s a function of e f f e c t i v e mean normal stress only, and given by, (5.21) 83 i n which, ( K 2 ) m a x i s a constant which depends on the r e l a t i v e density for a given s o i l . Seed et. a l , suggested that for sands ( K 2 ) m a x varies between 20 and 100. Hardin et. a l , (1972) suggested that the ultimate shear strength Tmax c a n ^ e calculated using i n s i t u e f f e c t i v e stresses and s t a t i c shear strength parameters such as c' and 0' (equation 4.7). They pointed out that for the l e v e l of dynamic s t r a i n (y c °3c n ( 5 . 3 5 ) I t should be remembered that i n t y p i c a l s o i l structures only lower modes of v i b r a t i o n govern the response, and therefore, i t i s unnecessary to include higher mode components. The equation ( 5 . 35 ) implies that the damping r a t i o increases l i n e a r l y with the frequency. Therefore, the response due to high frequency components of the input motion w i l l be damped out s i g n i f i c a n t l y . This i s one of the advantages of using s t i f f n e s s proportional damping r a t i o s for s o i l s . From the equation (5 .35 ) the damping r a t i o at the fundamental frequency i s given by, \ L = 0 . 0 0 2 5 W l ( 5 . 36 ) Typical periods of v i b r a t i o n of s o i l structures may vary between 0 .5 to 1.5 sec. This means that the t y p i c a l damping r a t i o for the fundamental mode at the s t a r t of the dynamic loading varies between 1% - 3%. 102 The average s t i f f n e s s of a s o i l structure during dynamic loading i s less than the s t i f f n e s s at the s t a r t of the dynamic loading. Therefore, when u s i n g the damping m a t r i x [c] computed u s i n g [ K ] t = 0 , the e f f e c t i v e damping r a t i o may be higher than the range shown above. 5 . 4 BOUNDARY CONDITIONS A p p r o p r i a t e boundary c o n d i t i o n s i n terms of f o r c e s or displacements have to be s p e c i f i e d at a l l boundaries. In the dynamic analyses involving earthquake e x c i t a t i o n , two types of bottom boundary conditions are often used: 1. A f i x e d bottom boundary located at the top of a r i g i d l a y e r . In general, base rock or a s t i f f s o i l layer can be assumed to be r i g i d . 2. A bottom boundary located at the top of a s o i l layer or soft rock with constant e l a s t i c properties. This type of boundary i s generally known as transmitting boundary. For both the above boundary conditions the earthquake motion i s s p e c i f i e d at the bottom boundary. If the second type of bottom boundary conditions i s used, the domain of i n t e r e s t need not be extended down to a r i g i d l a y e r . This procedure reduces the number of degrees of freedom leading to a reduction i n computational costs. In the method of analysis presented here only the f i r s t type of boundary condition i s considered. Three types of l a t e r a l boundary conditions are commonly used i n two-dimensional dynamic problems, involving f i n i t e element procedures: 103 1) boundaries are located s u f f i c i e n t l y far away from a structure so that wave r e f l e c t i o n does not occur during analysis or i s minimized. On these boundaries, forces, displacements or a combination of forces and displacements can be s p e c i f i e d . 2) viscous boundaries are used which attempt to absorb the ra d i a t i n g waves by a series of dashpots and springs with constant or va r i a b l e properties (Lysmer and Kuhlemeyer; 1969) 3) consistent boundaries can be provided close to the foundation of structures. These boundaries attempt to reproduce the far f i e l d response i n a way consistent with the f i n i t e element formulation used to model the dynamic problem. This i s accomplished by formulating a frequency dependent boundary s t i f f n e s s matrix which can be obtained by solving the e l a s t i c wave propagation problem i n a layered half-space (Lysmer and Wass; 1972). Of the three types of l a t e r a l boundaries, the analysis with the consistent boundaries i s by far superior to the others with respect to accuracy. In the analysis with consistent boundaries, only a small region needs to be considered, thus reducing the number of degrees of freedom. But unfortunately, the formulation i s s t r i c t l y applicable only to l i n e a r (or i t e r a t i v e l i n e a r ) problems and for solutions i n the frequency domain only. The l a t e r a l boundaries i n an incrementally e l a s t i c approach i n the time domain therefore should be located as far away from a structure as p r a c t i c a b l e . The usual way to model the l a t e r a l boundaries i s to allow the nodes on these boundaries to move only i n the h o r i z o n t a l d i r e c t i o n . 104 5.5 ANALYSIS OF SOIL - STRUCTURE SYSTEMS The response of a structure founded on a s o i l deposit i s a f f e c t e d by the l o c a l s o i l conditions. The peak ac c e l e r a t i o n , frequency content and the s p a t i a l d i s t r i b u t i o n of the response c h a r a c t e r i s t i c s may a l l be a f f e c t e d . By including the structure i n the f i n i t e element domain for the response an a l y s i s , the coupled seismic response of the s o i l and structure may be determined. The presence of the structure has two major e f f e c t s on the s o i l deposit. I t increases the e f f e c t i v e stresses and i t also provides a d d i t i o n a l i n e r t i a forces. Therefore, for s o i l s which exhibit non-linear stress dependent behaviour, an uncoupled analysis i n which s o i l and s t r u c t u r a l systems are uncoupled may not be a p p l i c a b l e . Uncoupled analysis cannot predict the response of buried structures where strong s o i l - s t r u c t u r e i n t e r a c t i o n occur. 5.5.1 S l i p Elements i n Dynamic Analysis Relative displacement which may occur between s o i l and structure during strong shaking can be modelled using s l i p elements. However, the s l i p element model described i n Chapter 4 has to be modified for use under dynamic loading conditions. Ideally, once the slippage stops i n a s l i p element, the top and bottom nodes of the element should have the same acce l e r a t i o n and v e l o c i t i e s (Nadim and Whitman; 1982). However, when computing accelerations and v e l o c i t i e s by numerically i n t e g r a t i n g the equation (5.5), d i f f e r e n t values w i l l be unavoidable for the top and bottom nodes. In order to overcome t h i s d i f f i c u l t y , when no slippage 105 occurs In a s l i p element, the v e l o c i t i e s and accelerations of the bottom nodes were made equal to those of the top nodes. 5.6 SOLUTION SCHEME A step by step s o l u t i o n scheme i s ca r r i e d out to obtain the dynamic response i n the time domain. A b r i e f outline of the procedure i s given below: 1) based on the c u r r e n t v a l u e s of G , T and v. max' max ' t at time t the tangent shear modulus G t i s calculated for a l l the elements using the current s t r e s s - s t r a i n curve, for eit h e r i n i t i a l l o a d i n g , u n l o a d i n g or r e l o a d i n g whichever i s appropriate. 2) the global s t i f f n e s s matrix [ K t ] t i s determined. 3) knowing the base a c c e l e r a t i o n value at (t+At), new values f o r {x}, {x}, {x} at time ( t + At) and increments, Ay and {ACT} are computed by employing any of the d i r e c t i n t e g r a t i o n methods to solve the equations (5.5) as de t a i l e d i n Appendix III 4) i f stress or s t r a i n r e v e r s a l occurs i n any element during this time i n t e r v a l , At, the dynamic analysis i s repeated f o r a shorter time i n t e r v a l . 5) using the shear s t r a i n increment, increments i n volumetric s t r a i n and then porewater pressure are computed using equation (5.28) and (5.29). 106 6) using increments i n residual porewater pressure U, as v i r t u a l loads s t a t i c analysis i s performed to determine current e f f e c t i v e stresses. These e f f e c t i v e stresses are then used to update and X ] n a x values. A computer program TARA-2 has been developed incorporating a l l these basic steps. 5.7 COMPUTATION OF POST EARTHQUAKE DEFORMATION I t i s often required to predict the displacements at various points on the s o i l structure at the end of an earthquake. This i s referred to as dynamic r e s i d u a l displacements. To compute these, an earthquake record with enough t a i l i n g zeros should be used so that the free damped v i b r a t i o n response can be included i n the a n a l y s i s . The c y c l i c and permanent components of displacement response for saturated sands and s i l t s , are assumed to occur under undrained conditions. There w i l l be a d d i t i o n a l deformations i n these s o i l s when the resi d u a l porewater pressure d i s s i p a t e s . In undrained simple shear c y c l i c tests on saturated sands and s i l t s the p o t e n t i a l volumetric s t r a i n e y < j , which occurs due to grain s l i p , i s r e f l e c t e d as r e s i d u a l porewater pressure. When the re s i d u a l porewater pressure d i s s i p a t e s , volumetric s t r a i n , e V (j occur i n the sample r e s u l t i n g i n settlement. The pr e d i c t i o n of deformations due to the d i s s i p a t i o n of r e s i d u a l porewater pressure can be accomplished as follows. 107 One method i s to treat t h i s problem as a two-dimensional consolidation problem with known i n i t i a l porewater pressures, and to solve for the deformation at di s c r e t e time i n t e r v a l s as the drainage occurs. A second method i s to compute deformations using the volumetric s t r a i n is accumulated at the end of the dynamic loading, eV(j» This can be accomplished by tr e a t i n g t h i s volumetric s t r a i n the same way the v o l u m e t r i c s t r a i n e y i n shear-volume c o u p l i n g was modelled i n Chapter 4. In the program TARA-2 the second method i s used. The f i n a l or post earthquake deformation a f t e r the r e s i d u a l porewater pressure has dissipated i s the sum of the deformation calculated ic from the m o d i f i e d e y ( j and the r e s i d u a l dynamic deformation predicted at the end of the dynamic a n a l y s i s . 108 CHAPTER 6 VERIFICATION OF THE METHOD OF ANALYSIS The v a l i d a t i o n of computational techniques requires good prototype data for comparison. The common shaking table test cannot represent the range of i n - s i t u pressures experienced by the s o i l elements i n the f i e l d . For a p p l i c a t i o n to f u l l scale design problems we require data from c e n t r i f u g a l models where s i m i l a r i t y can be achieved with the self-weight of the prototype. Here i n - s i t u stresses are d i r e c t l y scaled to those of the f u l l scale event and therefore, the stress dependent s t r e s s - s t r a i n properties of the s o i l ( e s p e c i a l l y cohesionless s o i l s ) can be matched i n model and prototype. Centrifuge modelling laws are used to deduce prototype response from model response. The p r i n c i p l e s of c e n t r i f u g a l modelling have been discussed i n d e t a i l by Schofield (1981). 6.1 CAMBRIDGE CENTRIFUGE TESTS A number of centrifuge tests on submerged islands were conducted by Lee (1983) using the Cambridge University Centrifuge. F u l l d e t a i l s of Cambridge centrifuge equipment and test procedures have been given by Schofield (1981). F i g . 6.1 shows the model i s l a n d used i n the t e s t s . The i s l a n d was a 90mm high with side slopes at 3:1 and a crest width of 200mm. The centrifuge a c c e l e r a t i o n used i n the tests was 40g. This means that the corresponding prototype i s l a n d i s of height 3.6m, with side slopes 3:1 and has a crest width of 8m. The s t r u c t u r a l loading on the LVDT 82260 90mm I 720mm • I I LVDT 82273 A C C 734 A C C 1244 • PPT 2 338 OPPT 2332 Mild s t e e l p l a t e s 3 PPT 2331 Concrete base ACC 2 5 8 ^ Container Base 750 mm 776 mm • D i r e c t i o n of p o s i t i v e a c c e l e r a t i o n measured on accelerometers LVDT Pore p r e s s u r e transducer * • Accelerometer Fig. 6.1. Submerged Island Showing Transducer Locations. o KO 110 i s l a n d was simulated by using mild s t e e l plates of various thicknesses. The i s l a n d rests on a concrete base which i n turn i s bolted to the centrifuge container. The base shaking was generated when the rotary arm on wheels follows a track mounted on the wall of the centrifuge chamber. The i s l a n d was instrumented by 6 DJB A23 p i e z o e l e c t r i c accelerometers, 6 Druck PDCK81 pore pressure transducers and 2 LVDTs. The l o c a t i o n of these instruments are also shown i n F i g . 6.1. The test s o i l was f i n e Leighton-Buzzard sand, mostly of s i z e passing between B.S. sieve si z e No. 120 and No. 200. The r e l a t i v e density was 60 - 70%, with e m a x = 1.03 and e m i n = 0.63 (Lee, 1983). In c e n t r i f u g a l modelling, i f the model pore f l u i d i s the same as the prototype pore f l u i d , excess pore pressures would be able to d i s s i p a t e N 2 times faster i n the model than i n the prototype, whilst the earthquake would only occur N times f a s t e r . Here N i s the scale factor of the c e n t r i f u g a l a c c e l e r a t i o n given as a r a t i o of g r a v i t y . Therefore, to model the prototype porewater drainage condition, i t i s necessary to use a pore f l u i d of v i s c o s i t y N times that of water i n the model t e s t . A s p e c i a l s i l i c o n e o i l was used as the model pore f l u i d i n Lee's t e s t s . In the tests c a r r i e d out on the i s l a n d , the t h e o r e t i c a l input wave-form was intended to be 12 sinusoidal pulses with a constant period of 0.5secs. However, the actual input motion to the i s l a n d was more com p l i c a t e d due to resonances and mechanical l i n k a g e c l e a r a n c e s i n t e r f e r i n g with the input motion, e s p e c i a l l y during the i n i t i a t i o n of the base motion. The r e s u l t s of two centrifuge tests were made a v a i l a b l e to the w r i t e r . The average contact prototype pressure on the islands f or these two tests were 15kPa (Test 1) and 31kPa (Test 2) r e s p e c t i v e l y . The input I l l accelerations measured by the accelerometer mounted on top of the concrete base for Test 1 and Test 2 are shown i n F i g . 6.2(a), and ( b ) . The recorded maximum base accelerations i n Test 1 and Test 2 were O . l l g and 0.17g r e s p e c t i v e l y . The i s l a n d responded very d i f f e r e n t l y i n these two t e s t s . The comparative study c a r r i e d out to predict performance of both the tests are reported i n t h i s chapter. 6.2 COMPARATIVE STUDY The sand used i n the model test was at an average r e l a t i v e density D r = 65%. T y p i c a l properties which are consistent for medium dense sand of D r = 65% for s t a t i c and dynamic analysis are given i n Table 6.1. The l i q u e f a c t i o n resistance curve for the sand obtained by using UBC simple shear apparatus without any s t a t i c bias (t g=0) i s shown i n F i g . 6.3. This l i q u e f a c t i o n resistance curve matches the predicted l i q u e f a c t i o n resistance curve when porewater pressure model constant, K r = 0.012. As explained i n Chapter 5, the porewater pressure model constants can be selected appropriately to account for the behaviour of samples with i n i t i a l s t a t i c s t r e s s , T G . Since test data on samples with T G, are not a v a i l a b l e , i t was decided to ignore the changes i n the pore-water pressure generation rates, and to account only for the changes i n l i q u e f a c t i o n r e s i s t a n c e . The l i q u e f a c t i o n curves which correspond to non zero t s/°" v 0 were obtained by s c a l i n g the laboratory l i q u e f a c t i o n curve. The s c a l i n g was done by using a v a i l a b l e laboratory data on medium dense Ottawa sand (Vaid, et a l . 1979). The changes i n the l i q u e f a c t i o n resistance curve can be s p e c i f i e d by associating the appropriate K r ~i 1 1 r 3.0 4.6 5.0 Time in Sees Fig. 6.2(a). Input Base Motion in Test 1. • rc l_ tu _| —i max = 0.17g — i 1 1 1 1 1 1 1 1 i r -2.0 3.0 4.0 5.0 6.0 7.0 Time in Sees 8.0 Fig. 6.2(b). Input Base Motion in Test 2. 113 Table 6.1 S o i l Properties Properties Total unit weight kN/m3 Bulk modulus constant Bulk modulus exponent constant n Shear modulus parameter ( K 2 ) m a x Angle of i n t e r n a l f r i c t i o n E f f e c t i v e cohesion C o e f f i c i e n t K Q a, b values used to compute [C] C^-K^ Constants Rebound modulus constants m,n S t a t i c Dynamic 18.1 800 0.4 19.3 38.0 0.0 0.45 18.1 high 55.0 38.0 0.0 0.0,0.005 0.75,0.79,0.459,0.73 0.43,0.62 04\ Number of Cycles to Initial Liquefaction, Fig. 6 . 3 . Liquefaction Resistance Curve of Medium Dense Leighton-Buzzard Sand. 115 v a l u e w i t h each s t a t i c shear s t r e s s r a t i o , "^/^vo* L i n e a r i n t e r p o l a t i o n may be used to get the K r value corresponding to an" intermediate -u s/a v o. For each model test two analyses were performed. One with no s l i p elements between s o i l and structure and the other with s l i p elements. The following s l i p element properties were used, K g = 6.3 x 10 5 kN/m2/m, % » 6.3 x 10 5 kN/m2/m C g = 0.0 0 g = 35° 6.2.1 Results of Test 1 The recorded a c c e l e r a t i o n time h i s t o r i e s of the model i s l a n d have very high frequency components which contain n e g l i g i b l e energy. This type of high frequency e l e c t r i c a l noise i s unavoidable i n centrifuge te s t i n g as i t may orig i n a t e due to ambient sources such as the e l e c t r i c motor d r i v i n g the centrifuge, and also due to centrifuge v i b r a t i o n s . Dean (1983), suggested i t i s necessary to f i l t e r out very high frequency components from output q u a n t i t i e s . The computed and recorded a c c e l e r a t i o n time h i s t o r i e s shown here have been smoothed once using a three point average scheme, suggested by Dean (1983). In using t h i s scheme, the current value at any point i n time i s computed as the sum of 1/4 of the value the previous point, 1/2 the value of the current point and 1/4 of the value of the next point. F i g . 6.4 to F i g . 6.6 show the smoothed recorded and computed acceleration time h i s t o r i e s of accelerometers A1244, A1225 and A734. During the f i r s t 1.5 seconds of shaking low accelerations with very high 116 ~i 1 1 1 r aO 1.0 2.0 3.0 4.0 5.0 Tine in Sees 6.0 r "5.0 8.0 Fig. 6.4(a). Recorded Acceleration of ACC1244 in Test 1. -i 1 r Q.D 1.0 2-0 3.0 4.0 5.0 Time I n Sees r T "7.0 e.c Fig. 6.4(b). Computed Acceleration of ACC1244 in Test 1. (with and without Slip Elements). 117 !" c Tine in Sees 6.0 -r r — T \ 7.0 8.0 Fig. 6.5(a). Recorded Acceleration of ACC1225 in Test 1. Acceleration of ACC1225 * Teat 1. Fig. 6.5(b). Computed (with and without SUp Elements). Fig. 6.6(a). Recorded Acceleration of ACC734 in Test 1. Pig. 6.6(b). Computed Acceleration of ACC734 in Test 1. (with and without Slip Elements). 119 frequency were recorded i n a l l accelerometers. A f t e r that the amplitude of a c c e l e r a t i o n s t e a d i l y increased as i n the case of input motion, upto 5.5 seconds and then subsided. The a c c e l e r a t i o n time h i s t o r i e s computed by TARA-2 with and without s l i p elements are very s i m i l a r and therefore, only one of them i s shown i n F i g . 6.4.(b), through F i g . 6.6(b). The frequency and magnitude of the computed ac c e l e r a t i o n response are very s i m i l a r to corresponding recorded response values. Table 6.2 shows the computed and recorded maximum acc e l e r a t i o n of a l l three accelerometers. Table 6.2 Recorded and Computed Maximum Accelerations Maximum Acceleration % 8 Instrument Location Computed by TARA-2 Recorded Without S l i p Elements With S l i p Elements A1244 13.3 11.6 11.6 A1225 15.9 12.5 12.5 A734 13.9 12.7 12.7 The comparison between recorded and computed maximum acceleration values are very good. F o u r p o r e w a t e r p r e s s u r e d e v e l o p m e n t p l o t s o b t a i n e d experimentally and computed by TARA-2 are presented i n Figures 6.7 (a), (b), (c) and (d) . In t h i s test very low porewater pressures were developed. During the low l e v e l shaking of the f i r s t second, the response 120 Recorded Computed With and Without S l i p Elements I I 1 1 1 1 1 1 1 1 1 4.0 5.0 -6.0 "7.0 6.0 at Tine In Sees Fig. 6.7(a). Recorded and Computed Porewater Pressure of PPT2330 in Test 1. x x x— Recorded Computed Without S l i p Elements Computed With S l i p Elements i 1 r i.O 5.0 Time In Seas Fig. 6.7(b). Recorded and Computed Porewater Pressure of PPT68 in Test 1. 121 Fig. 6.7(c). Recorded and Computed Porewater Pressure of PPT2338 in Teat 1. Recorded Computed With and Without S l i p Elements i e.c Fig. 6.7(d). Recorded and Computed Porewater Pressure of PPT2342 in Test 1. 122 of the i s l a n d i s e s s e n t i a l l y e l a s t i c and porewater pressures recorded are the instantaneous response to changes i n t o t a l s t r e s s . Such porewater pressures develop from the e l a s t i c coupling of s o i l and water. Later during the period of more severe shaking, p l a s t i c volumetric s t r a i n s occur r e s u l t i n g i n the development of r e s i d u a l porewater pressures which are independent of the instantaneous states of s t r e s s . The recorded porewater p r e s s u r e s d u r i n g t h i s time have both r e s i d u a l and i n s t a n t a n e o u s components. A f t e r 6 seconds of shaking the input motion subsides over the next two seconds to zero. During t h i s time the magnitude of the instantaneous component of porewater pressure i s small. Drainage may occur during the e x c i t a t i o n depending on the drainage c h a r a c t e r i s t i c s of the sand. Since generation of porewater pressure a f t e r 6 seconds of e x c i t a t i o n i s very small, changes i n porewater pressure can be caused only by drainage during t h i s time. A close examination of recorded porewater pressures a f t e r 6 seconds of e x c i t a t i o n reveals that the porewater pressures i n a l l four locations except at the transducer P2342, which i s located at the middle of the i s l a n d are more or les s a constant. At t h i s l o c a t i o n porewater pressure increases s l i g h t l y i n d i c a t i n g movement of water towards the center of the i s l a n d . However, since these changes are small i t i s reasonable to assume that drainage i n the i s l a n d i s n e g l i g i b l e during the base e x c i t a t i o n . TARA-2 computes only r e s i d u a l porewater pressures, so there are no f l u c t u a t i o n s due to changes i n instantaneous stress l e v e l s i n the computed curves. Furthermore, since no drainage was assumed during the e x c i t a t i o n the computed r e s i d u a l porewater p r e s s u r e s i n c r e a s e co n s i s t e n t l y . However, the rate of increase i n porewater pressures during 123 low l e v e l e x c i t a t i o n which occur before 1.5 seconds and a f t e r 6 seconds i s r e l a t i v e l y small. When r i g i d j o i n t connection i s assumed between a heavy, s t i f f s t r u c t u r a l element and an adjacent s o i l element i n an a n a l y s i s , the dynamic str a i n s developed i n the adjacent s o i l element are very small due to compatibility requirements i n displacements. However, by introducing s l i p elements between the structure and s o i l , the r e l a t i v e movement which may occur between s o i l and structure during strong base e x c i t a t i o n can be accounted f o r . The r e s u l t s from TARA-2 analyses, with and without s l i p elements, ind i c a t e that computed porewater pressures are d i f f e r e n t only i n the transducers located just below the structures. The predicted porewater pressures just below the structure i n the analysis which incorporates s l i p elements, are s l i g h t l y higher than the analysis that assumes r i g i d connection between s o i l and structure. The comparison between recorded and computed porewater pressures i s good. Only four porewater pressure time h i s t o r i e s from the model tests are a v a i l a b l e . However, maximum re s i d u a l porewater pressures, which can be interpreted as the mean recorded values a f t e r the e x c i t a t i o n has subsided, are a v a i l a b l e for a l l s i x transducers. Table 6.3 shows the computed and recorded maximum re s i d u a l porewater pressures at a l l the t r a n s d u c e r l o c a t i o n s . 6.2.2 Results of Test 2 The smoothed recorded and computed a c c e l e r a t i o n time h i s t o r i e s of accelerometers A1244, A1225 and A734 are shown i n F i g s . 6.8 through 6.10. Table 6.3 Recorded and Computed Maximum Residual Porewater Pressures Residual Porewater Pressure, kPa Transducer Location Computed by TARA-2 Recorded Without S l i p Elements With S l i p Elements P2330 0.4 0.4 0.4 P2331 0.4 0.3 0.4 P2332 4.0 2.9 3.7 P68 3.0 1.5 4.4 P2338 2.4 0.6 1.8 P2342 4.0 4.2 4.2 125 T I 2.0 T" 3.D 4.0 5.0 Time In Sees Fig. 6.8(a). Recorded Acceleration of ACC1244 in Test 2. Fig. 6.8(b). Computed Acceleration of ACC1244 in Test 2. (with Slip Elements). + 25% 1 1 1 1 1 1 1 1 1 1 1 1 I 1 I I 13 LC 2.0 3.0 4.0 5.0 6.0 ">.0 8.0 Time In Sees Fig. 6.9(a). Recorded Acceleration of ACC1225 in Test 2 . + 2 6% "1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ' a0 1.0 2.0 3.0 4.0 5.0 6.C "7.0 8.TJ Time i n Sees Fig. 6.9(b). Computed Acceleration of ACC1225 in Test 2 . (with Slip Elements). 127 0.0 "i i 1 1 1 r 2.0 3.0 4.0 5.0 Time in Sees i 1 1 r 6.0 "7.0 s.c Fig. 6.10(a). Recorded Acceleration of ACC734 in Test 2. Time in Sees Fig. 6.10(b). Computed Acceleration of ACC734 in Test 2 (with Slip Elements). 128 The base e x c i t a t i o n for Test 2, shown i n F i g . 6.2(b), has the following c h a r a c t e r i s t i c s . Low l e v e l e x c i t a t i o n occurs i n the f i r s t second and then the amplitude of acceleration increases s t e a d i l y to a constant maximum value i n the next two seconds. This maximum amplitude i s maintained for 2 seconds and then i t subsides over the next 3 seconds. The acceleration records, except the record obtained i n accelerometer A734 which i s located on top of the structure, show the v a r i a t i o n of acceleration amplitude being s i m i l a r to that of the input motion. The acce l e r a t i o n h i s t o r i e s recorded on the top of the structure dropped to very low values a f t e r 4 seconds of e x c i t a t i o n . Unlike Test 1, the response computed with s l i p elements were found to be d i f f e r e n t from those computed without s l i p elements. The response analysis also showed that when s l i p elements are used s i g n i f i c a n t amount of s l i p occurs between s o i l and st r u c t u r e . The computed accelerations without s l i p elements are lower than the values computed when s l i p elements are used i n the an a l y s i s . In the computed ac c e l e r a t i o n h i s t o r i e s shown i n F i g s . 6.8 to F i g . 6.10, only accelerations computed using s l i p elements are presented. Lines of constant accelerations have been drawn i n Fig s . 6.8 to 6.10 to aid the i n t e r p r e t a t i o n of r e s u l t s . The v a r i a t i o n of computed ac c e l e r a t i o n h i s t o r i e s are very s i m i l a r to that of the input motion. The maximum recorded accelerations which were observed to be associated with high frequencies (indicated by sharp s p i k e s ) and maximum computed a c c e l e r a t i o n s f o r a l l three accelerometers i n the i s l a n d are shown i n Table 6.4. The comparison between recorded and computed accelerations are not good. However, a closer look at the other peak values of corresponding recorded and computed ac c e l e r a t i o n h i s t o r i e s suggests s a t i s f a c t o r y comparison. 129 The recorded a c c e l e r a t i o n at any time has two components: the acc e l e r a t i o n component transmitted through the s o i l from the base and the acce l e r a t i o n component transmitted to s o i l through side walls and top of the centrifuge container due to the container i t s e l f v i b r a t i n g . The computed ac c e l e r a t i o n h i s t o r y accounts only for the a c c e l e r a t i o n component transmitted through the s o i l . The accelerations transmitted through the container have high frequency components. The presence of these high frequency components may be responsible for the discrepancies between the recorded and computed ac c e l e r a t i o n h i s t o r i e s . Table 6.4 Recorded and Computed Maximum Accelerations Maximum Ac c l e r a t i o n , % g Accelerometer No. Computed by TARA-2 Recorded Without S l i p Elements With S l i p Elements A1244 24.0 15.1 18.2 A1225 42.5 15.5 23.1 A734 23.9 15.8 18.2 Four recorded and computed porewater pressure development plots are presented i n F i g s . 6.11 (a), (b), (c) and (d). In th i s test unlike Test 1, very high porewater pressures were developed. During the low l e v e l of e x c i t a t i o n of the f i r s t second very low porewater pressures were recorded. With the onset of more severe shaking, very high porewater 130 Fig. 6.11(a). Recorded and Computed Porewater Pressure of PPT2330 in Test 2. Fig. 6.11(b). Recorded and Computed Porewater Pressure of PPT68 in Test 2. Fig. 6.11(c). Recorded and Computed Porewater Pressure of PPT2338 in Test 2. Fig. 6.11(d). Recorded and Computed Porewater Pressure of PPT2342 in Test 2. 132 pressures were developed i n a l l transducers except i n the transducer P2330. A close examination of the recorded porewater pressure plots a f t e r 6 seconds of e x c i t a t i o n reveals that, except for the porewater pressure transducer P2342, which i s located at the middle of the i s l a n d , s i g n i f i c a n t d i s s i p a t i o n of porewater pressure occurred. This i s because i n t h i s test very high porewater pressures were developed leading to high pressure gradients. The tranducer P2342 behaved d i f f e r e n t l y because i t i s too far from the free draining boundaries and at t h i s l o c a t i o n inward flow of water occurs. The comparison between computed and recorded porewater pressures are very good for the two transducers (P2330, P2342) which are located well inside the i s l a n d . At these transducer locations the analysis with and without s l i p elements gave very s i m i l a r r e s u l t s . The porewater pressures predicted under the structure by the analysis without s l i p elements are very low. But, when s l i p elements were provided the comparison between predicted and computed porewater pressures improved. When s l i p elements were used i n the response analysis high shear s t r a i n s developed i n the elements below the structure r e s u l t i n g i n higher porewater pressures i n those elements. Table 6.5 shows the computed and recorded maximum r e s i d u a l porewater pressures at a l l the transducer l o c a t i o n s . 6.3 APPLICABILITY OF THE METHOD OF ANALYSIS The recorded r e s i d u a l porewater pressures are interpreted as the steady increase i n porewater pressures. Therefore, any high frequency noise from ambient sources on recorded values do not a f f e c t the comparison Table 6 .5 Recorded and Computed Maximum Residual Porewater Pressures Maximum Residual Porewater Pressure, kPa Transducer No. Computed by TARA-2 Recorded Without S l i p Elements With S l i p Elements P2330 1.0 1.5 1.5 P2331 0.9 1.1 1.1 P2332 10.5 12.0 12.1 P68 38.0 6.4 38.1 P2338 18.0 2.2 . 18.9 P2342 22.0 19.8 21.3 134 between recorded and computed porewater pressures. However, since the comparison i n a c c e l e r a t i o n response i s made point by point, any high frequency noise a f f e c t s the comparison. Therefore, caution should be exercised when comparing maximum accelera t i o n s . When a sand sample i s subjected to undrained loading, an abrupt change i n the d i r e c t i o n of stress path occurs, as the stress path approaches the f a i l u r e l i n e (Ishihara, et a l . 1975). The points at which various stress paths change d i r e c t i o n abruptly are assumed to l i e on a s t r a i g h t l i n e , c a l l e d the phase transformation l i n e . The slope of the phase transformation l i n e i s a few degrees less than the f a i l u r e l i n e . Any c y c l i c (or monotonic) loading beyond the phase transformation l i n e may r e s u l t i n very low e f f e c t i v e stresses i n very loose sands and increased e f f e c t i v e stresses due to d i l a t i o n i n medium dense or dense sands. The hyperbolic s t r e s s - s t r a i n r e l a t i o n s h i p assumed i n t h i s thesis i s s t r i c t l y a pplicable only for the region of stress space below the phase transformation l i n e . F i g s . 6.12 (a) and (b) show the stress paths i n a q, p' plot that were followed by four elements which correspond to the locations for which porewater pressure time h i s t o r i e s are a v a i l a b l e . The stress paths reported here are for the analysis i n which s l i p elements were included. The stress paths followed by the elements i n Test 1 are well below the f a i l u r e l i n e , where as i n Test 2, two elements are on the f a i l u r e l i n e for sometime during the dynamic loading. Under these circumstances the v a l i d i t y of response computed i n Test 2 a f t e r elements have reached the f a i l u r e l i n e i s questionable. Fig. 6.12(a). Effective Stress Paths in Test 1 Fig. 6.12(b). Effective Stress Paths in Test 2. O N 137 CHAPTER 7 APPLICATION OF THE METHOD OF ANALYSIS: TANKER ISLAND RESPONSE 7.1 INTRODUCTION Man made i s l a n d s of c o h e s i o n l e s s s o i l s have been used extensively as d r i l l i n g platforms for o i l and gas exploration i n the Beaufort Sea. Recently, as exploration has moved to deep waters, more complex forms of i s l a n d construction procedures have been introduced. The caisson-retained i s l a n d (De Jong, et a l . (1978), and s t e e l tanker i s l a n d are two t y p i c a l examples. These newer type of construction procedures greatly reduce the amount of f i l l material required and also reduce some of the hazards of wave loading on exposed i s l a n d beaches. The maximum set down water depth for the caisson-retained i s l a n d and the tanker islands depth i s f i x e d , generally around 6 to 9 metres. Therefore, i n the case of deep water a underwater sand berm i s constructed up to the set down water depth. Most of the sand berms are constructed by dumping sand excavated by suction dredges from an offshore and/or onshore borrow p i t and pumped as a s l u r r y through a p i p e l i n e d i r e c t l y onto the l o c a t i o n of the i s l a n d . Once the sand berm i s ready, a series of caissons or tanker i s brought to the l o c a t i o n and b a l l a s t e d onto the berm, and b a c k f i l l e d with sand, gravel or water. The d r i l l i n g i s then c a r r i e d out from the upper structure. Because of the nature of the i s l a n d construction, the dumped sand i s often loose and therefore the deformation, s t a b i l i t y and 138 l i q u e f a c t i o n p o t e n t i a l of the i s l a n d berm during earthquakes are of major concern. 7.2 ANALYSIS OF A TYPICAL TANKER ISLAND F i g . 7.1, shows, schematically a tanker i s l a n d . This i s l a n d i s provided with a cover of about 2m of rock f i l l . T y p i c a l properties of rock f i l l and sand for s t a t i c analysis are given i n Table 7.1. TABLE 7.1. S t a t i c S o i l Properties Properties Rock F i l l Sand Total unit weight kN/m3 18.7 18.1 Bulk modulus constant K^ 1000 800 * Bulk modulus exponent constant n 0.40 0.40 Shear modulus parameter ( K 2 ) m a x 24.0 16.0 Angle of i n t e r n a l f r i c t i o n 38.0 32.0 E f f e c t i v e cohesion 0.0 0.0 C o e f f i c i e n t K Q 0.45 0.45 The tanker i s assumed to weigh 200,000 metric tons when f u l l y ballasted with plan dimensions 170m and x 60m and 21m high. In the case presented here, i t i s assumed that the hyperbolic dynamic s t r e s s - s t r a i n r e l a t i o n s h i p has equal shear strength i n both Fig. 7.1. Schematic of Tanker Island 140 d i r e c t i o n s of shearing. T y p i c a l properties used for the dynamic analysis are give i n Table 7.2 TABLE 7.2. Dynamic S o i l Properties Properties Rock F i l l Sand 1300 V. High 0.4 70 45 0.8, 0.79, 0.459 and 0.730 Rebound modulus Constants m,n 0.43, 0.62 During shaking the rock f i l l i s assumed to be free draining and no drainage i s assumed i n dumped sand. A very high value was assigned to B t to simulate very low compressibility imparted to the saturated sand by the water i n the pores which i s not allowed to drain. Liquefaction resistance curves are required for d i f f e r e n t s t a t i c stress r a t i o s i n the i s l a n d . These are s p e c i f i e d by associating the appropriate K r value with each s t a t i c shear stress r a t i o t s / o " v o . The values used i n the example are presented i n Table 7.3. Bulk modulus constant ft Bulk modulus parameter n Shear modulus parameter (K 2) C^ •* CK Constants TABLE 7.3. T s / a v o and K r values T /a* K s vo r 0.0 0.004 0.10 0.015 0.20 0.05 The s t a t i c shear stress r a t i o s for the example problem considered here vary between 0.0 to 0.13 and therefore, the values of K r corresponding to r a t i o s T g / o y o , above 0.15 are not necessary. The input motion used for the analyses i s the S00E a c c e l e r a t i o n component of the Imperial V a l l e y Earthquake of May 18, 1940 scaled to O.lg. The input motion was applied at the bottom boundary of the i s l a n d . Three dynamic analyses were performed: i s l a n d alone, i s l a n d plus structure with and without s l i p elements. The properties of s l i p elements were selected so that some s l i p could occur between the structure and the i s l a n d . The s l i p element properties were assumed to be, K s = 6.3 x 10 5 kN/m2/m, K N = 6.3 x 10 5, kN/m2/m C s = 0 and 0 g = 30° A complete response study of the tanker i s l a n d could be c a r r i e d out by d i s c r e t i z i n g the enti r e domain into f i n i t e elements. However, since the s t i f f n e s s of tanker wall i s very much higher than s o i l , the tanker and i t s contents would respond l i k e a r i g i d box. In view of 142 t h i s , the tanker and i t s contents were modelled as a uniform r i g i d box. The s t i f f n e s s of s t r u c t u r a l elements were selected as 10 3 of the rock f i l l elements. 7.2.1. Results f o r Tanker Island Problem One of the factors which influence the development of r e s i d u a l porewater pressure i s c y c l i c shear stress (or c y c l i c s t r a i n ) . Since the generation of r e s i d u a l porewater pressure i s possible only i n the sand, the maximum dynamic shear stresses induced i n the dumped sand along section T-T which runs through the centre of the i s l a n d , are shown i n F i g . 7.2 for a l l three cases. This fi g u r e indicates that higher dynamic shear stresses are developed when the tanker i s i n place due to the i n e r t i a forces on the tanker. The induced shear stresses i n the dumped sand when s l i p i s allowed to occur between structure and adjacent s o i l are s l i g h t l y higher than when no s l i p elements were provided. When s l i p occurs, the magnitudes of shear stress that can be transmitted to the structure i s l i m i t e d , dictated by the shear strength of the s l i p elements. Therefore, i s l a n d responses with and without s l i p occurring between the structure and the i s l a n d may be expected to d i f f e r . Despite very high c y c l i c shear stresses generated i n the sand when tanker i s i n place, the greatest porewater pressure r a t i o s are developed i n the unloaded i s l a n d ( F i g . 7.4). This i s because the v e r t i c a l over-burden pressures are very much greater when the tanker i s present, so that, the c y c l i c shear stress r a t i o s , T c v / O y 0 , which i s the most important parameter c o n t r o l l i n g the development of porewater pressure i n a given sand, are a c t u a l l y smaller ( F i g . 7.3). I t can be r e a d i l y seen that Maximum Dynamic Shear Stress (kN/m2) 0 20 40 60 I 1 1 1 Fig. 7.2. Distribution of Maximum Dynamic Shear Stress in Sand. Fig. 7.3. Distribution of Maximum Dynamic Shear Stress Ratio in Sand. 144 the d i s t r i b u t i o n of maximum c y c l i c shear stress r a t i o s are proportional to the d i s t r i b u t i o n of r e s i d u a l porewater pressure r a t i o s . F i g . 7.4, shows the r e s i d u a l p o r e w a t e r r a t i o H/aLQ d i s t r i b u t i o n for a l l three analyses. The r e s i d u a l porewater pressure r a t i o s obtained for the unloaded i s l a n d are higher than for the loaded i s l a n d . The r e s u l t s obtained i n the analysis without the tanker can be viewed as the solutions at a section where the influence of the structure i s n e g l i g i b l e . Therefore, the d i s t r i b u t i o n of r e s i d u a l porewater pressure r a t i o s when the structure i s i n place w i l l vary from lower values at the middle to higher values as one moves away from the s t r u c t u r e . Same conclusions were drawn by Yoshimi and Tokimatsu (1977) who studied the response of a r i g i d structure subjected to base e x c i t a t i o n on a shaking ta b l e . The r e s i d u a l porewater pressure d i s t r i b u t i o n given i n the analysis with s l i p elements i s co n s i s t e n t l y higher than the analysis without s l i p elements. This i s because lower shear stresses are induced i n the l a t t e r case. The d i s t r i b u t i o n of maximum dynamic shear s t r a i n s f or section 1-1 i s shown i n F i g . 7.5. Even though the shear stresses induced i n the unloaded i s l a n d are smaller than the loaded i s l a n d , higher shear s t r a i n s have developed i n the unloaded i s l a n d . This i s because of two f a c t o r s . F i r s t l y , the i n - s i t u overburden stresses i n the unloaded i s l a n d are very much smaller and therefore the s t r e s s - s t r a i n curve for a given element i s s o f t e r . Secondly, when an e f f e c t i v e s t r e s s - s t r a i n r e l a t i o n s h i p i s used, any increase i n r e s i d u a l porewater pressure w i l l soften the s t r e s s - s t r a i n curve. The generation of higher r e s i d u a l porewater pressures and low over burden pressures have contributed to high shear s t r a i n s i n the unloaded i s l a n d . Fig. 7.4. Distribution of Residual Porewater Pressure Ratio in Sand. Fig. 7.5. Distribution of Maximum Dynamic Shear Strain in Sand. 146 The maximum hori z o n t a l displacements which occur during the earthquake for the section 1-1 are shown i n F i g . 7.6. Much smaller dynamic displacements are computed when the tanker i s i n place. When s l i p elements were provided, s l i p occurred between the s o i l and structure and the displacements are about twice the r e s u l t s obtained without s l i p elements. F i g . 7.7 a and b show the post earthquake deformations i n the X and Y d i r e c t i o n s . The post earthquake displacement i s the sum of the dynamic r e s i d u a l displacement and the displacement due to volumetric * s t r a i n component e y ^ • Two o b s e r v a t i o n s can be made from the r e s u l t s presented i n t h i s f i g u r e . F i r s t l y , the amount of the post earthquake deformations i n the X- d i r e c t i o n are proportional to the maximum dynamic deformations. Secondly, the X-component of the displacement i s of the same order of magnitude as of the Y-component of the displacement ( s e t t l e m e n t ) . The main c o n t r i b u t i o n to the X-component of the displacement comes from the dynamic r e s i d u a l displacement and for the settlement the main contribution i s from the volumetric s t r a i n component ft evd* In the dynamic response of structures the maximum induced acceleration i n the structure i s one of the main design concerns. The maximum induced accelerations given by TARA-2 i n the structure with and without s l i p elements are 0.15g and 0.17g. This means that, i f s l i p i s prevented, the acc e l e r a t i o n induced may be higher by as much as 15% of the acceleration when s l i p i s allowed. The maximum ac c e l e r a t i o n computed on top of the unloaded i s l a n d i s 0.15g. Fig. 7.6. Distribution of Maximum Dynamic Displacement. Post Earthquake X Disp.. (cm). Post Earthquake Y Disp., (cm) Fig. 7 . 7 . Post Earthquake X and Y Displacements. 00 1 4 9 One of the ways of presenting dynamic response i s to present i t i n terms of response spectra. Response spectra for displacement, v e l o c i t y and acceleration are often presented for the motions at the base of the stru c t u r e . These r e s u l t s are then used by the engineers to predict the behaviour of structures and also to compute design forces, such as base shear. F i g . 7.8 shows the response spectrum for a c c e l e r a t i o n of the motion at the berm surface for a l l three cases considered. The damping r a t i o used i n the computation was 3%. Inspection of t h i s f i g u r e suggests that for the example problem considered here, the a c c e l e r a t i o n response predicted using the response spectrum of the unloaded i s l a n d w i l l be higher for structures with a very low period. However, for the structures with a period greater than 0.5 s e c , the response predictions w i l l be s i m i l a r . The predominant motion of a tanker during e x c i t a t i o n are s l i d i n g and rocking. The r e l a t i v e importance of these two modes can be studied by comparing r e s u l t s obtained by a two-dimensional and one-dimensional response a n a l y s i s . F i g . 7.9, shows the computed d i s t r i b u t i o n of r e s i d u a l porewater pressure r a t i o u / o y o from three two-dimensional analyses which were reported e a r l i e r and also the d i s t r i b u t i o n from a one-dimensional response a n a l y s i s . The one-dimensional case considered was the i s l a n d with tanker, without any s l i p elements. The r e s u l t s c l e a r l y show that the maximum U/Oy0 of some elements may be predicted as low as 30% of those predicted by a two-dimensional response a n a l y s i s . This means a response analysis which neglects the rocking mode of v i b r a t i o n i s non-conservative. However, i t should be mentioned that the tanker considered i n t h i s example i s very t a l l (21m) and, r i g i d , and therefore the rocking mode of v i b r a t i o n may have been more important than usual. 600 1 r- r —i 500 Island Alone — 400 | t " J l II •*—. Island + Structure (Slip) -300 \f\ Island + Structure ?A (No Slip) -200 -100 0 10 2-0 3-0 4-0 5-0 Period, (sec) Fig. 7 . 8 . Acceleration Response Spectrum for the Motion at Berm Surface. 7.9. Distribution of Residual Porewater Pressure Ratio in Sand. 151 7.3 SOME PRACTICAL CONCLUSIONS In practice a number of s i m p l i f y i n g assumptions are made when computing the response of structures founded on s o i l deposits. The procedure outlined by the National Building Code of Canada has two basic steps. The f i r s t step i s to compute the response of the s o i l deposit alone for the given e x c i t a t i o n . The second step i s to compute the response of the structure, to the base accelerations obtained i n step one. In predicting the performance of the structure, the r e s u l t s such as porewater pressure, induced s t r a i n l e v e l etc., which are obtained from step one are also considered. This means that the code i n essence suggests the s o i l - structure systems be uncoupled and analysed independently. Figures 7.2 to 7.8 c l e a r l y show that the response of structures computed using the procedures outlined i n the National Building Code of Canada may be i n e r r o r . The presence of the structure has two basic influences on the s o i l deposit. It increases the e f f e c t i v e stresses and i t also provides a d d i t i o n a l i n e r t i a forces. Therefore, for s o i l s which exhibit non-linear stress dependent behaviour, the uncoupled analysis proposed by the code may not be a p p l i c a b l e . From t h i s t y p i c a l example, three basic conclusions can be drawn. F i r s t of a l l i t r a i s e s questions about the merit of any response analysis based on uncoupled s o i l - s t r u c t u r e systems. Secondly, one-dimensional representation of the domain which neglects the rocking degrees of freedom may not be applicable to t a l l , heavy and r i g i d s t r u ctures. T h i r d l y i t demonstrates the importance of incorporating s l i p elements i n the a n a l y s i s . Because of the great weight of the caissons or tankers, and 152 t h e i r large l a t e r a l dimensions, s o i l - s t r u c t u r e i n t e r a c t i o n e f f e c t s w i l l always be important. In these type of problems a coupled analysis of the i s l a n d and structure i s required. 153 CHAPTER 8 SUMMARY AND CONCLUSIONS 8.1 SUMMARY The main purpose o f t h i s r e s e a r c h was t o d e v e l o p a two-d i m e n s i o n a l s t a t i c and s e i s m i c r e s p o n s e a n a l y s i s o f s o i l d e p o s i t s i n c l u d i n g s o i l - s t r u c t u r e i n t e r a c t i o n . The new method f o r s t a t i c and dynamic a n a l y s e s can be pe r f o r m e d i n e i t h e r e f f e c t i v e o r t o t a l s t r e s s modes or a c o m b i n a t i o n o f b o t h . Non-l i n e a r s t r e s s - s t r a i n b e h a v i o u r of s o i l was m o d e l l e d by u s i n g an i n c r e m e n t a l l y e l a s t i c a p p r o a c h i n w h i c h t a n g e n t shear modulus and t a n g e n t b u l k modulus were t a k e n as the two " e l a s t i c " p a r a m e t e r s . The m a t e r i a l r e s p o n s e i n s h e a r was assumed t o be h y p e r b o l i c w i t h M a s i n g b e h a v i o u r d u r i n g u n l o a d i n g and r e l o a d i n g . Response t o changes i n mean normal s t r e s s was assumed t o be n o n - l i n e a r , e l a s t i c and s t r e s s dependent. When a s t a t i c a n a l y s i s i s per f o r m e d i n the t o t a l s t r e s s mode, t h e s h e a r s t r e n g t h , T „ , 0 „ , G m a v , a n d t a n g e n t b u l k m o d u l u s , B t , o f an element a r e k e p t c o n s t a n t t h r o u g h o u t t h e a n a l y s i s . The tan g e n t shear modulus, G t, i s m o d i f i e d f o r c o r r e s p o n d i n g s h e a r s t r a i n s d e v e l o p e d d u r i n g the a n a l y s i s . When e f f e c t i v e s t r e s s mode i s u s e d , the p a r a m e t e r s , t m a x > ^max a n c * B t a r e c o m p u t e d f r o m t h e e f f e c t i v e s t r e s s e s . The e f f e c t o f d i l a t i o n d u r i n g s h e a r on volume change i s t a k e n i n t o a c c o u n t . I n t he s t a t i c a n a l y s i s p r o p o s e d h e r e , g r a v i t y may be s w i t c h e d on 154 at once for the completed s o i l structure or the construction sequence can be modelled by layer a n a l y s i s . The s t r e s s - s t r a i n conditions determined by the s t a t i c analysis give the i n - s i t u stress conditions before the dynamic a n a l y s i s . S l i p or contact elements have been incorporated i n the analysis to represent the i n t e r f a c e c h a r a c t e r i s t i c s between s o i l and s t r u c t u r a l elements. The properties of the s l i p element were assumed to be e l a s t i c p e r f e c t l y p l a s t i c , with f a i l u r e at the interface given by the Mohr Coulomb f a i l u r e c r i t e r i o n . In the dynamic e f f e c t i v e stress response a n a l y s i s , r e s i d u a l porewater pressures are calculated using a modification of the model proposed by M a r t i n e t . a l , (1975). The parameters, G m a x> and Tmax» a r e m°dified for the e f f e c t s of r e s i d u a l porewater pressure. The dynamic response study includes the p r e d i c t i o n of post earthquake deformations. An extensive study c a r r i e d out to v e r i f y the proposed method of analysis using centrifuge test data suggests that the proposed method can be successfully used to predict seismic response of structures. Seismic response of a t y p i c a l tanker i s l a n d computed by t h i s method i s presented. 8.2 CONCLUSIONS The work that has been presented i n t h i s thesis leads to the following conclusions. 1. A consistent and r e l i a b l e method for computing transient and 155 permanent deformations i n two-dimensional s o i l structures i s needed. 2. A two-dimensional dynamic response analysis which takes into account the non-linear h y s t e r e t i c stress-dependent properties of s o i l s , has been developed i n terms of both t o t a l and e f f e c t i v e stresses. 3. The method has been v e r i f i e d by comparing data from centrifuged models with predictions of the method. Comparison between predicted and measured response parameters i s generally very good. 4. Allowing for s l i p to occur between s o i l and s t r u c t u r a l elements i s very important. Analyses which allow for s l i p have c o n s i s t e n t l y lead to higher displacements i n the structure and higher porewater pressures i n the s o i l deposit. 5. The method has been applied to compute seismic response of a t y p i c a l tanker i s l a n d . The re s u l t s of t h i s study suggests that i t i s important that the response of structures founded on s o i l deposits be analysed as a coupled s o i l - s t r u c t u r e systems. 6. The v a l i d i t y of a one-dimensional response analysis instead of a two-dimensional analysis f or tanker type of structures i s questionable. The porewater pressures predicted by using a one-dimensional response analysis model may be as low as 30% of those predicted by a two-dimensional response a n a l y s i s . 156 SUGGESTIONS FOR FURTHER STUDY 1 . A d d i t i o n a l comparative studies should be c a r r i e d out so that greater confidence could be placed on the v a l i d i t y of th i s method. Comparative studies may be performed with data from centrifuge tests or f i e l d studies. 2. Sandy materials exhibit p a r t i a l s t a b i l i z a t i o n at low confining pressures due to d i l a t i o n . Therefore, i n the response evaluation near l i q u e f a c t i o n , i t i s important that the method of analysis Include the d i l a t a n t behaviour of sands. 3. In the response evaluation of more permeable s o i l s , drainage during the seismic loading may be s i g n i f i c a n t and procedures should be developed to take t h i s into account. 157 REFERENCES 1. Byrne, P.M., (1979), Class Notes: Numerical Methods i n S o i l Mechanics (CE573). University of B r i t i s h Columbia, Vancouver, B.C., Canada. 2. Byrne, P.M., and Eldridge, T.L., (1982), "A Three Parameter D i l a t a n t E l a s t i c S t r e s s - S t r a i n Model for Sand", International Symposium on Numerical Models i n Geomechanics, Zurich, Sept., pp73-79. 3. Chang, C.S., (1982), "Residual Undrained Deformation from C y c l i c Loading", Journal of the Geotechnical Engineering, ASCE, V o l . 108, GT4, A p r i l , pp637-646. 4. Chern, J . C , (1981), " E f f e c t of S t a t i c Shear on Resistance to Liq u e f a c t i o n " , M.A.Sc. Thesis, The U n i v e r s i t y of B r i t i s h Columbia, Vancouver, May. 5. C h r i s t i a n , J.T. and Boehmer, J.W., (1970), "Plane S t r a i n Consolidation by F i n i t e Elements", Journal of the S o i l Mechanics and Foundation Engineering D i v i s i o n , ASCE, V o l . 96, No. SM4, July, ppl435-1455. 6. Clough, R.W., and Penzien, J . , (1975), "Dynamics of Structures", McGraw H i l l Book Co., New York. 7. De Jong, J.J.A., and Bruce, J . C , (1978), "Design and Construction of a Caisson-Retained-Island D r i l l i n g Platform for the Beaufort Sea", Proceedings, 10th Annual Offshore Technology Conference, Houston, Texas, May 8-11, OTC Paper No. 3294. 8. Dean, E.T.R., (1983), "FLY-14 Program Suite: i n F l i g h t Data Handling and Analysis Manual", Report, Cambridge Unversity Engineering Department, Cambridge, U.K. 158 9. Desai, C.S., and Abel, J.F., (1972), "Introduction to the F i n i t e Element Method", L i t t o n Educational Publishing, Inc. 10. Desai, C.S., and C h r i s t i a n , J.T. (1977), "Numerical Methods i n Geomechanics", McGraw H i l l Book Co., New York. 11. Duncan, J.M., Byrne, P.M., Wong, K.S., and Mabry, P., (1978), "Strength, S t r e s s - S t r a i n and Bulk Modulus Parameters f o r F i n i t e Element Analyses of Stresses and Movements i n S o i l Masses", Report No. WCB/GT/78-02, to the National Science Foundation, A p r i l . 12. Finn, W.D.L. and M i l l e r , R.I.S., (1973), "Dynamic Analyses of Plane Non-Linear Earth Structures", 5th World Conference on Earthquake Engineering, Rome, Session ID, Paper No. 42, pp360-367. 13. Finn, W.D.L., and Byrne, P.M., (1976), "Estimating Settlements i n Dry Sands During Earthquakes", Canadian Geotechnical Journal, V o l . 13, No. 4, pp355-363. 14. Finn, W.D.L., Lee, K.W., and Martin, G.R., (1977), "An E f f e c t i v e S t r e s s Model f o r L i q u e f a c t i o n " , J o u r n a l of the G e o t e c h n i c a l Engineering D i v i s i o n , ASCE, V o l . 103, No. GT6, Proc. Paper 13008, June, pp517-533. 15. Finn, W.D.L., Martin, G.R., and Lee, K.W., (1978a), "Comparison of Dynamic Analyses of Saturated Sands", Proc. ASCE Geotechnical Engineering D i v i s i o n , S p ecialty Conference on Earthquake Engineering and S o i l Dynamics, Pasadena, C a l i f o r n i a , June 19-21, pp472-491. 16. Finn, W.D.L., Lee, K.W., Maartman, C.H., and Lo, R., (1978b), " C y c l i c Pore Pressures Under Anisotropic Conditions", Proc. ASCE Geotechnical Engineering D i v i s i o n , S p ecialty Conference on Earthquake Engineering and S o i l Dynamics, Pasadena C a l i f o r n i a , June 19-21, pp457-468. 159 17. Finn, W.D.L., (1981), "Liquefaction P o t e n t i a l Development Since 1976", Proceedings, International Conference on Recent Advances i n Geotechnical Earthquake Engineering and S o i l Dynamics, St. Louis, Missouri, A p r i l 26-May 2, pp655-681. 18. Finn, W.D.L., I a i , S., and Ishihara, K., (1982), "Performance of A r t i f i c i a l Offshore Islands under Wave and Earthquake Loading", Offshore Technology Conference, V o l . 1, OTC, Paper 4220, Houston, Texas, May, pp661-672. 19. Finn, W.D.L., (1983), "Analyses of Cumulative Deformations Under C y c l i c Loading", Invited Challenge Contributed to National Science Foundation Workshop on Expected Research i n S o i l Engineering, Blackberg, V i r g i n i a , U.S.A., Aug. 22-24. 20. Goodman, R.E., and Seed, H.B., (1966), "Earthquake Induced Displacements i n Sand Embankments", Journal of the S o i l Mechanics and Foundations D i v i s i o n , ASCE, V o l . 92, SM2, March, ppl25-146. 21. Goodman, R.E., Taylor, R.L., and Brekke, TL.L, (1968), "A Model for the Mechanics of Jointed Rock", Journal of the S o i l Mechanics and Foundations D i v i s i o n , ASCE, Vol. 94, SM3, May, pp637-659. 22. Hansen, B., (1958), "Line Ruptures Regarded as Narrow Rupture Zones, Basic Equation Based on Kinematic Considerations', Proceedings, Brussels Conference on Earth Pressure Problems, Brussels, V o l . 1, pp. 39-48. 23. Hardin, B.O., and Drenevich, V.P., (1972), "Shear Modulus and Damping i n S o i l s , Design Equations and Curves", Journal of the S o i l Mechanics and Foundation D i v i s i o n , ASCE, V o l . 98, SM7, Proc. Paper 9006, July, pp667-692. 160 24. I a i , S., and Finn. W.D.L., (1982), "DONAL-2, A Computer Program for Dynamic One Dimensional Analyses of Slope Layers with Energy Transmitting Boundary", S o i l Dynamics Group, U n i v e r s i t y of B r i t i s h Columbia, Vancouver, B.C., Feb. 25. Ishihara, K., Tatsuoka, F. and Yasuda, S., (1975), "Undrained Deformation and Liquefaction of Sand Under C y c l i c Stresses", S o i l s and Foundations, Vol. 15, No. 1, pp29-44. 26. Konder, R.L., and Zelasko, J.S.,, (1963), "A Hyperbolic S t r e s s - S t r a i n Formulation of Sands", Proceedings of 2nd Pan American Conference on S o i l Mechanics and Foundation Engineering, V o l . 1, B r a z i l , 1963, pp289. 27. Kulhawy, F.H., Duncan, J.M., and Seed, H.B., (1969), " F i n i t e Element A n a l y s i s of S t r e s s e s and Movements i n Embankments During Construction", Geotechnical Engineering Research Report No. TE-69-4, Department of C i v i l Engineering, Unviersity of C a l i f o r n i a , Berkeley, Nov. 28. Lee, F.H., (1983), " P a r t i a l Liquefaction i n Centrifuge Model Embankment i n an Earthquake", M.Phil. Thesis, Engineering Department Cambridge Uni v e r s i t y , Cambridge, J u l y . 29. Lee, K.W., (1965), " T r i a x i a l Compressive Strength of Saturated Sands Under Seismic Loading Conditions", Ph.D Thesis, U n i v e r s i t y of C a l i f o r n i a , Berkeley, ppl-521. 30. Lee, K.W., (1975), "Mechanical Model for the Analysis of Li q u e f a c t i o n of Horizontal S o i l Deposits", Ph.D Thesis, Department of C i v i l Engineering, U n i v e r s i t y of B r i t i s h Columbia, Vancouver, B.C., Canada, Sept. 161 31. Lee, K.W., and Finn, W.D.L., (1975), "DESRA-1: Program for Dynamic E f f e c t i v e Stress Response Analysis of S o i l Deposits Including L i q u e f a c t i o n Evaluation", S o i l Mechanics Series, No.36, Dept. of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver, B.C., Canada. 32. Lysmer, J . , and Kuhlemeyer, R.L., (1969), " F i n i t e Dynamic Model for I n f i n i t e Media", Journal of the Engineering Mechanics, D i v i s i o n , ASCE, V o l . 95, EM4, Sept. pp859-877. 33. Lysmer, J . , and Waas, G., (1972), "Shear Waves i n Plane I n f i n i t e Structures", Journal of Engineering Mechanics D i v i s i o n , ASCE, V o l . 98, EMI, Feb., pp85-105. 34. Martin,. G.R., Finn, W.D.L., and Seed, H.B., (1975), "Fundementals of Liquefaction Under C y c l i c Loading", Journal of the Geotechnical Engineering D i v i s i o n , ASCE, V o l . 101, GT5, May, pp423-438. 35. Masing, G., (1926), "Eigenspannungen and Verfestigung beim Messing", Proceedings, 2nd International Congress of Applied Mechanics, Zurich, Switzerland. 36. Mroz, Z., Norris, U.A., and Zi e n k i e i c z , O.C, (1979), "Application of an Anisotropic Hardening Model i n the Analysis of E l a s t o - P l a s t i c Deformations of S o i l s " , Geotechnique, 29(1), ppl-34. 37. Nadim, F., and Whitman, R.V. (1982), " S l i p Elements for Earthquake Response", Proceedings, 4th International Conference on Numerical Methods i n Geomechanics, Edmonton, Alberta, Canada, pp61-67. 38. Newmark, N.M., (1959), "A Method of Computation f o r St r u c t u r a l Dynamics", Journal of the Engineering Mechanics D i v i s i o n , ASCE, V o l . 85, EM3, Ju l y . 162 39. Newmark, N.M., (1965), " E f f e c t s of Earthquakes on Dams and Embankments", 5th Rankine Lecture, Geotechnique 15, No.2, ppl39-160. 40. Newmark, N.M., and 'Rosenblueth, E., (1971), "Fundementals of Earthquake Engineering", Prentice-Hall Inc., Englewood, C l i f f , N.J., ppl62-163. 41. Ozawa, Y., and Duncan, J.M., (1973), "ISBILD: A Computer Program for A n a l y s i s of S t a t i c S t r e s s e s and Movements i n Embankments", Geotechnical Engineering Research Report No. TE-73-4, Department of C i v i l Engineering, University of C a l i f o r n i a , Berkeley, Dec. 42. Prevost, J.H., (1978), "Mathematical Modelling of S o i l S t r e s s - S t r a i n Strength Behaviour", Proceedings, 3rd International Conference on Numerical Methods i n Geomechanics, Aachen, Germany, A p r i l 6, pp347-361. 43. Pyke, R., Seed, H.B., and Chan, C.K., (1975), "Settlement of Sands Under M u l t i d i r e c t i o n a l Shaking", Journal of the Geotechnical Engineering ASCE, V o l . 101, No. GT4, A p r i l , pp379-398. 44. Rahman, M.S., Seed, H.B., and Booker, J.R., (1977), "Pore Pressure Development Under Offshore Gravity Structures", Journal of the Geotechnical Engineering D i v i s i o n , ASCE, V o l . 103, GT12, D e c , ppl419-1436. 45. Robertson, P.K., (1982), " I n - s i t u Testing of S o i l with Emphasis on i t s A p p l i c a t i o n to Liquefaction Assessment", Ph.D Thesis, submitted to Department of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver, B.C., Canada, Dec. 46. Schnabel, P.B., Lysmer, J . , and Seed, H.B., (1972), "SHAKE: A Computer Program for Earthquake Response Analysis of H o r i z o n t a l l y 163 Layered S i t e s " , Report No. EERC 72-12, Earthquake Engineering Reseach Center, University of C a l i f o r n i a , Berkeley, Dec. 47. S c h o f i e l d , A.N., (1981), "Dynamic and Earthquake Geotechnical C e n t r i f u g e M o d e l l i n g " , I n t e r n a t i o n a l Conference on Recent Development i n Geotechnical Earthquake Engineering and S o i l Dynamics, Missouri, U.S.A., A p r i l 28-May 2. 48. Seed, H.B., and Lee, K.L., (1969), "Porewater Pressures i n Earth Slopes Under C y c l i c Loading Conditions", Proceeding, 4th World Conference on Earthquake Engineering, Chile, V ol. I l l , p p l l l . 49. Seed, H.B., and I d r i s s , I.M., (1970), " S o i l Moduli and Damping Factors for Dynamic Response Analyses", Report No. EERC 70-10, Earthquake Engineering Research Center, Un i v e r s i t y of C a l i f o r n i a , Berkeley, December. 50. Seed, H.B., Lee, K.L., I d r i s s , I.M., and Makdisi, F.I., (1973), "Analysis of the Slides i n the San Fernando Dams During the Earthquake of Feb. 9, 1971", Report No. EERC 73-2, Earthquake Engineering Research Center, University of C a l i f o r n i a , Berkeley, June. 51. Seed, H.B., Pyke, R. and Martin, G.R. (1975), "Analysis of the E f f e c t of M u l t i - d i r e c t i o n a l Shaking on the Liquefaction C h a r a c t e r i s t i c s of Sands", Report NO. EERC 75-41, Earthquake Engineering Research Center, University of C a l i f o r n i a , Berkeley, C a l i f . Dec. 52. Seed, H.B., (1979), "Considerations i n the Earthquake-Resistant Design of E a r t h and R o c k f i l l Dams", 19th Rankine L e c t u r e , Geotechnique 29, No.3, pp215-263. 53. Seed, H.B., (1983), "Earthquake-Resistant Design of Earth Dams", Proceedings of Seismic Design of Embankments and Caverns, pp41-64. 164 54. S e l i g , E.T., and Chang, S.C. (1981), " S o i l F a i l u r e Modes i n Undrained C y c l i c Loading", Journal of the Geotechnical Engineering D i v i s i o n , ASCE, Vol. 107, GT5, May, pp539-551. 55. S e r f f , N., Seed, H.B., Makdisi, F.I., and Chang, C.Y., (1976), "Earthquake Induced Deformations of Earth Dams", Report No. EERC 76-4_, Earthquake Engineering Research Center, U n i v e r s i t y of C a l i f o r n i a , Berkeley, Sept. 56. U n i t e d S t a t e s N a t i o n a l Research C o u n c i l (1982), "Earthquake Engineering Research-1982", Report by Committee on Earthquake Engineering Research, National Acadamy Press, Washington, D.C. 57. Vaid, Y.P., and Finn, W.D.L., (1979), " E f f e c t of S t a t i c Shear on Liquefaction P o t e n t i a l " , Journal of the Geotechnical Engineering D i v i s i o n , ASCE, Vol. 105, GT10, OCT., ppl233-1246. 58. Vaid, Y.P., Byrne, P.M., and Hughes, J.M.O., (1981), " D i l a t i o n Angle and Liquefaction P o t e n t i a l " , Journal of the Geotechnical Engineering D i v i s i o n , ASCE, V o l . 107, GT7, J u l y . 59. Varadarajan, A., and Mishra, S.S., (1980), "Stress-Path Dependent S t r e s s - S t r a i n Volume Change Behaviour of G r a n u l a r S o i l " , I nternational Symposium on S o i l s under C y c l i c and Transient Loading, Swansea, ppl09-119. 60. V e r r u i j t , A., (1977), "Generation and D i s s i p a t i o n of Porewater Pressures", F i n i t e Elements i n Geomechanics, Edited by, Gudehus, G., John Wiley and Sons, pp293-319. 61. Wilson, E.L., Farhomand, I., and Bathe, K.J., (1973), "Non-Linear Dynamic Analysis of Complex Structures", International Journal of Earthquake Engineering and S t r u c t u r a l Dynamics, V o l . 1, pp241-252. 165 62. Wedge, N.E. (1977), "Problems In Non-Linear Analysis of Movements i n S o i l s " , M.A.Sc. Thesis, Submitted to the Dept. of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver, B.C., Canada. 63. Yoshimi, Y., and Tokimatsu, K. (1977), "Settlement of Buildings on Saturated Sand During Earthquakes", S o i l s and Foundations, Japanese Society of S o i l Mechanics and Foundation Engineering, V o l . 17, No. 1, March, pp23-38. 64. Zienkiewicz, O.C., and Cheung, Y.K., (1967), F i n i t e Element Method i n S t r u c t u r a l and Continuum Mechanics", McGraw H i l l Book Company. 166 APPENDIX I FINITE ELEMENT FORMULATION In the f i n i t e element analysis the entire domain of i n t e r e s t i s divided into a f i n i t e number of elements. The v a r i a t i o n of displacement within a f i n i t e element i s assumed to be given by, {tf} = [N] {6} where ( A l . l ) {tf} = displacement vector, giving x and y displacements at any point within an element, here TJ*" = {u,v} {6} = displacement vector, giving x and y displacements of the nodes, and [N] = i n t e r p o l a t i o n function. The type of element used i n the a n a l y s i s i s a 4 node isoparametric element. The term isoparametric implies common (i s o - ) parametric d e s c r i p t i o n of the unknown displacements and the geometry of the element. The same i n t e r p o l a t i o n functions N^ are used to express both the displacement and the geometry of the element. Isoparametric element formulation has a number of advantages; i t o f f e r s e f f i c i e n t i n t e g r a t i o n s , and d i f f e r e n t i a t i o n s and i t can handle curved and a r b i t r a r y geometrical shapes. The i n t e r p o l a t i o n function [N] can be selected such that i t can be expressed i n natural coordinates (s,t) which i s a system of co-ordinates i n t r i n s i c to an element ( F i g . A l . l ) . 168 In the isoparametric concept, the coordinates of a point i n the element i s given by, {X} = [ N ] {X±} (A1.2) If we consider a four-node q u a d r i l a t e r a l isoparametric element, the matrix [N] i s composed of, „ . d - s ( l - t ) _ ( l + s ) ( l - t ) 1 4 W2 4 N = (i+s)(i+t) _ q - s ^ (i+t) 3 4 4 4 Now equations ( A l . l ) and (A1.2) can be rewritten as, u l v l u 2 j U , r N , 0 N 2 0 N 3 0 N^ 0 -I v 2 M L 0 N x 0 N 2 0 N 3 0 N ^ u 3 (A1.3) v 3 u 4 v 4 and, ( x } = f N l 0 N 2 0 N 3 0 S l f 0 , y, V L 0 N x 0 N 2 0 N 3 0 N^-1 x 3 X l y i x 2 ?3 X 4 169 It can be seen that the transformation shown i n (A1.4), maps the q u a d r i l a t e r a l element into a square as shown i n F i g . A l . l . If plane s t r a i n conditions are assumed the s t r a i n vector e i s s i m p l y g i v e n by _et = { e x , e y , Y x y } where e x and e y a r e normal s t r a i n s and y x y 1 S t n e shear s t r a i n . These s t r a i n s are given i n terms of displacements as, x .e. e X _ M ox e y av ay Y ' xy 8u ay ax s t r a i n matrix e from A1.3 and A1.5, e x ax 0 ax 0 ax 0 wk ax 0 = 0 ay 0 a N 2 ay 0 3N 3 ay 0 ay aN 1 aN1 aN 2 8N 2 a \ at^ Y 'xy ey ax ay ax ay ax ay ax > {e} = [B] {6} But = f ( s , t ) and also x,y are functions of s,t. should be computed using the following r e l a t i o n s h i p , u. u„ (A1.5) 2 (A1.6) (A1.7) So any d e r i v a t i v e aN ± aN ± ax aN ± ay as ax as ay as (A1.8) In matrix form, the global and l o c a l d e r ivatives can be written as, 170 5N± ax ax aN ± 8N ± — — as as - — - — as ax r - i ax (AI.9) BN. 3x ay. 3N, L J J o N . l A i . y ; at a t 8 t ay ay where, [ j ] = Jacobian matrix The de r i v a t i v e s of with respect to x, y can now be obtained from A1.9, as, 8N± aN± 5± • [ J ] _ 1 3N, a y a t where, [ j ] - 1 = Inverse of Jacobian matrix J , which i s simply, a t a s Now from A1.4, the Jacobian matrix can be written as, ax By. 1 ^ 1 t^± R -I _ as as _ i as x i • i as y i ... L J J " . 8N u 8N (A1.12) a x a y . _ i fc f \ a t a t ^ 5 t x i A o t y± 171 The components of matrix [ j ] can be evaluated since oN^/os, oN^/Qt and (x^, y^) are known. So, knowing the matrix J , [ j ] - 1 also can be computed. Say, = [ h i Ti2) •21 (A1.13) 22 It should be noted that these I ^ j are f ( s , t ) . Therefore, s u b s t i t u t i n g this i n equation ( A L I O ) , y i e l d s , oN ± oNjL By" t l 1 1 L21 22 8N. as" aN. i at (A.14) The matrix [B] (equation A1.6) which relates the s t r a i n vector to nodal displacements, has derivatives of i n t e r p o l a t i o n function with respect to x and y. Now knowing these derivatives from equation (A1.4), [B] can be rewritten as, [ B ] = 3x8 aN, aN •11 0 S + X12 9 t | o i aN I aN 2i a s + 122 at | A n bi I aN, aN, + i . 2i as 22 at aN. aN, + i 3N 2 3N 2 X n + xi2 a~T aN. 9N„ 12 at | x2i a£ + i 22 at | 172 0 3N 3 aN 3 0 hi as + at aN 2 3N 3 3N 3 as + I22 at 0 X 2 1 as + I22 a t a N 2 aN 3 m3 aN 3 B N 3 as + 112 at I21 as + I 2 2 at i n as + h2 a t aN^ a i ^ | [ n o s + T i 2 a t | aN„ aN, o A 2 i a s + x22 a t aN 4 mh x n a s - + x i 2 IT I 2 1 as~ + *22 IT | (A1.15) The s t r e s s e s and s t r a i n s a r e c o n n e c t e d t h r o u g h e l a s t i c i t y m a t r i x g i v e n by, {a'} = [ D ] {e} (A1.16) where, [ a ' ] = e f f e c t i v e s t r e s s v e c t o r . F o r 2D p l a n e s t r a i n c o n d i t i o n s , i t i s g i v e n by, i o ' } * - = { 0 X , O y , T x y } [ D ] = e l a s t i c i t y m a t r i x f o r 2D p l a n e s t r a i n c o n d i t i o n , [ D ] -B + ^ G Sym B - f G 4 B + j G 0 0 (A1.17) 173 Using the v i r t u a l work p r i n c i p l e , the i n t e r n a l work (Wj^) done by a p p l y i n g i n f i n i t e s i m a l v i r t u a l nodal displacement {£} i s , WIN = /// M dv (A1.18) where {e} = v i r t u a l s t r a i n s due to v i r t u a l displacement {£} and {a} = t o t a l stress vector (A1.19) The t o t a l stress vector can be s p l i t into e f f e c t i v e stress and porepressure vectors. i . e . : {a} = {a'} + {uQ} E f f e c t i v e stress Porepressure (A1.19) Vector Vector Here {a'= {a x, a y , T x y } and {xxj*- = {u D, u D , o} i n which u Q i s the porewater pressure. Now, s u b s t i t u t i n g {a} from equations (Al.19) and (A1.18), one gets, WIN = /// {e}* [{o'} + {u }] dv (A1.20) V Substituting for {a'} from A1.16 th i s reduces to, WIN = /// M + K) D V (AI.21) v 174 using equation (A1.7), {e} can be replaced by [B] {£}, then, W I N - J / J [B]' [ D ] [B] {6} + J J J {6} f c [ B ] t { u J dv (A1.21) V V But, external work done by the load vector {p} r i d i n g through the v i r t u a l displacement {?>}, i s simply, WEX = t p J (A1.22) The v i r t u a l work p r i n c i p l e gives, W I N WEX i .e. {*}*{*} = {*}'/// [B] t[D][B] dv {6} v (A1.23) + {«}'/// [ B ] ' {uQ} dv V Noting that, dv = dxdydz and also = | j | dsdtdz where dz = thickness of the element. Aft e r s u b s t i t u t i n g t h i s i n (A1.23), and d i v i d i n g both sides by {?} t one gets, 175 {P} = [k]{6} + [k*] {uD} (A1.24) where, [k] = element s t i f f n e s s matrix l i t = [ / J [ B] [D][B] l Jl d s d t ] (A1.25) (assuming unit thickness) and, [k ] = porewater pressure matrix 1 1 t = [/ / [B] | j | dsdt] (A1.26) -1-1 The integrations shown above have to be evaluated numerically. The Gauss inte g r a t i o n technique has been employed and the number of points used are 2 x 2 . The formulation presented here i s for any l i n e a r e l a s t i c m a t e r i a l . For incrementally e l a s t i c a n a l y s i s , the displacements, stresses and moduli values should be simply replaced by incremental displacements, incremental stresses and tangent moduli. Af t e r evaluating the incremental load vector j), element tangent s t i f f n e s s m a t r i x [ k t ] , porewater p r e s s u r e m a t r i x [ k * ] , and a l s o estimating the incremental porewater pressure u Q , for a l l the elements the global incremental load-displacement r e l a t i o n s h i p can be formed. This w i l l lead to, {P} = [K t]{A} + [K*]{U} (A1.27) 176 i n which {P}, [ K T ] , {A}, [ K ] and {u} are relevant v a r i a b l e s i n glo b a l axes. By solving t h i s equation the displacement f i e l d {A} can be obtained, and i t can be used to c a l c u l a t e element s t r a i n s and stresses using equations (A1.6.) and (A1.16) r e s p e c t i v e l y . Since the shape function gives l i n e a r s t r a i n v a r i a t i o n within an element, the s t r a i n s and therefore, stresses vary within an element. For convenience, average stress and s t r a i n of an element are computed at the centre of gravity of the element. In Chapter 4, and Chapter 5, i t i s required to express s t r a i n s and stresses i n terms of nodal forces. R e c a l l from equation (A1.24) the nodal forces are given by, {P} = J / J !>]' [D] [B] dv {£} (A1.28) V But, str a i n s are connected to the matrix [B J i n equation (A1.7), {e} = [B] {6} (A1.29) Therefore, from (A1.28) and (A1.29), the nodal forces can be written i n terms of str a i n s as, {P} " J J J [ B ] 1 [D] {e} dv (A1.30) V Now, from equation (A1.16), k ' l = [D] { e } (A1.31) 177 and from equations, (A1.30), and (A1.31) the nodal forces can be written i n terms of stresses as, {P} = /// [ B f {a'} dv (A1.32) V The equation (A1.30) and (A1.32) can now be used to express s t r a i n s and stresses i n terms of element nodal forces. 178 APPENDIX I I STIFFNESS MATRIX FORMULATION FOR THE SLIP ELEMENT As o u t l i n e d i n Chapter 4, that the f o r c e - d i s p l a c e m e n t r e l a t i o n s h i p at any point within a s l i p element has been assumed to be given by, f K 0 w (fS} = U S K ] [ W S ] (A2.1) i . e . , f_ = ICQ W where, f g and f n = shear and normal stresses Kg, ^ = j o i n t s t i f f n e s s i n shear and normal d i r e c t i o n s wg, wn = shear and normal displacements The e l a s t i c stored energy, 0 E i n a s l i p element due ( F i g , (A2.1) to applied forces can be obtained by, 0 E = 2 t f dZ (A2* 2) i n which L i s the t o t a l length of the s l i p element. A factor half i s included because the r e l a t i o n s h i p between f_ and w i s assumed to be l i n e a r . From (A2.1), 0g now can be written as, Fig. A2.1. Slip Element. s = tangential direction U j = tangential displacement o f node i n = normal direction V j = normal displacement of node i 180 0 E = 2 /o W ' k o W d A ( A 2 * 3 ) Since the v a r i a t i o n of displacements (u, v) within an element i s l i n e a r the displacement at any point which i s at a distance, A, from node I on the bottom edge I J of the element i s , u U) = £ U j + (1 - L> U I ( A 2 - 4 > bottom In a s i m i l a r manner the following equations can be written f o r , u (A), v (1) and v U ) , top bottom top l .e. and u U) = L ° K + ( 1 ~ L ? " L (A2.5) top v(A) = L V J + ( 1 + L } V I (A2.6) bottom v (x) = L VK + ( 1 " L? V L ( A 2 ' 7 ) top where, u^, v^ r e f e r to displacement i n tangential and normal d i r e c t i o n of the nodes I, J , K and L . Shear and normal displacements at any point are, w = u (A) - u (A) (A2.8) top bottom and, 181 w n = v (A) - v U ) (A2.9) top bottom Now s u b s t i t u t i n g f o r u t o p > u b o t t o m , v t o p a n d v b o t t o m f r o m equations (A2.4) to (A2.7), i n equations (A2.8)-and (A2.9), w = H 1 -i> ~i i u - ft] (A.2.10) and, w_ = [-(1 - f) f f (1 ( A 2 . l l ) From equation (A2.1), w i s , U l v, w = [ S ] = — L w J -A 0 -B 0 B 0 A 0 [ ] 0 -A 0 -B 0 B 0 A J \ V K \ (A2.12) i n which, A = 1 " L A N D B - L In matrix form the equation (A2.12) i s , 182 w = C 6 2x1 2x§ 8X1 ( A 2 . 1 3 ) where, = the i n t e r p o l a t i o n matrix and 6^ = nodal displacement matrix Substitution of ( A 2 . 1 3 ) i n ( A 2 . 3 ) gives, K = 9 t s T cl K CJ d* E 2 •'o o o o ( A 2 . 1 4 ) Performing the matrix m u l t i p l i c a t i o n , T C k C o o o - A 0 0 - A - B 0 0 - B K 0 B 0 [ 8 0 B 0 K A 0 n 0 A r - A 0 - B 0 B 0 A On - A 0 - B 0 B 0 A - l K A 2 0 A B K 0 - A B K 0 - A 2 K 0 s s s s 0 A 2 K 0 A B K 0 - A B K 0 - A 2 K n B 2 K n n i A B K 0 0 - B 2 K 0 - A B K 0 s s s s 0 A B K 0 B 2 K 0 - B 2 K 0 - A B K n n n - A B K 0 - B 2 K 0 B 2 K 0 A B K 0 s s s B 2 K s 0 - A B K 0 - B 2 K 0 0 A B K n n n n - A 2 K 0 - A B K 0 A B K 0 A 2 K 0 s s s s 0 - A 2 K 0 - A B K 0 A B K 0 A 2 K ( A 2 . 1 5 ) To perform the in t e g r a t i o n shown i n equation A 2 . 1 4 , one should know integrations of, 183 / L A 2 dA, f L B 2 dA and f L AB dA J o J o J o These are simply, J L J o A 2 dA = 3 o (1 - ft2 dA L 3 / L J 0 B 2 dA = = J o ( f )2 dA L 3 (A2.16) and, Jo AB dA = • J L ; o ( 1 - f t f dA 6 Now the equation (A2.14) can be written as, 0 E " 2 ° T K s n 6 ( A 2' 1 7> where, sn 6 2K s 0 K s 0 -K s 0 -2K s 0 2K n 0 K n 0 -K n 0 -2K n 2K s 0 -2K s 0 -K s 0 2K n 0 -2K n 0 -K n Sym 2K s 0 2K n K s 0 2K s 0 K n 0 2K R e c a l l that e l a s t i c stored energy 0g i n the formulation of a l i n e a r e l a s t i c f i n i t e element i s , 0E = j 6 T K 6 (A2.18) 184 Where K i s the s t i f f n e s s matrix of the f i n i t e element. Now, comparing (A2.17) and (A2.18), the s t i f f n e s s matrix for s l i p element can be deduced as K„_. —sn 185 APPENDIX III STEP BY STEP INTEGRATION For the proposed incrementally e l a s t i c dynamic analysis i n the time domain the " e l a s t i c " properties have to be modified for every time step. The e l a s t i c parameters depend on the l e v e l of s t r a i n i n the deposit. Therefore, the displacement f i e l d i n the deposit should be evaluated at every time step. This requires that the incremental dynamic equilibrium equations (equation 5.5) have to be solved numerically for every time step. Newmark's method (1959) of step by step i n t e g r a t i o n i s very popular and extensively used i n dynamic analyses. This method b a s i c a l l y provides numerical s o l u t i o n i n time domain, where the s o l u t i o n i s advanced by one d i s c r e t e step at a time. In t h i s method, two parameters <* and p are used so that the v e l o c i t y and displacement at time t+At can be expressed i n terms of a c c e l e r a t i o n , v e l o c i t y and displacement at time t, and of the known acc e l e r a t i o n at time t+At. For convenience l e t us define that, T = t+At Then the r e l a t i o n s h i p i n terms of <* and 8 are, {X} T = {X}t + (1 - «) At {X} t + ocAt {X} T (A3.1) and, 186 {X} T = {X}t + At {X} t + (± - 6) ( A t ) 2 {X} t + B A t 2 {X} T (A3.2) Newmark (1959) proposed that « = 1/2 and B = 1/4 be used for an u n c o n d i t i o n a l l y s t a b l e i n t e g r a t i o n procedure, which i n c i d e n t a l l y corresponds to a constant average acceleration method of i n t e g r a t i o n . I f = = 1/2 and 8 = 1/6 are used then t h i s method gives a l i n e a r v a r i a t i o n of acc e l e r a t i o n within the time step. Re-writing the incremental equilibrium equations from Chapter 5, [M] {AX} + [C] {AX} + [ K t ] t {AX} = {AP} (A3.3) Now su b s t i t u t i n g f o r , {AX} {AX} {AX} (A3.4) and, i n equation (A3.3), one gets, [M] {^ - X j + [C] {Xj, - X t} + [ K t ] t {^ - X t} = {AP} (A3.5) From e q u a t i o n s (A3.1) and ( A 3 . 2 ) , {x}^, and {^}^ can be expressed i n terms of other variables as follows, 187 [ | £ 2 + [ c ] p £ + [ * t ] t ] (AX} = {AP} + [M] {E} t + [C] {F} t W T = m 2 [ { A X } " A t { x } t " (i ~ P> A t 2 W J ( A 3 * 6 ) and, W T - (X} t + (!--) At {X} (A3.7) + ^ [{AX} - At{x}t - B) At* {X}J Substituting for {x}-£ and {x},p i n equation (A3.5), W 2 [ {AX} - A t { X } t - (j - 6 ) At* {X} t - BAt 2 {x}J + [C] [ ( 1 — ) At {X}fc + ^ [{AX} - At {X} t (A3.8) - £ - B) At2 {X}J + [ K t ] t {AX} = {AP} Co l l e c t i n g terms, and defining following s i m p l i f y i n g symbols, and, . W t = f W t + ' ( 2 p " 1 } A t ^ t ( A 3 ' 1 0 ) the equation (A3.8) can be reduced to, (A3.11) 188 Recall from Chapter 5, Section 5.1.2, cor r e c t i o n forces should be applied to restore t o t a l equilibrium at time T. The { P c o r r } evaluated using equation (5.11) can be added to r i g h t hand side of the equation (A3.11). Then the equation (A3.11) i s , ^ + [ c ] p t r + [ K T ] T ] w = {AP} + [M] {E} t + [C] {F} t (A3.12) where, {AP} = {AP} + { P c o r r } (A3.13) The only unknown i n the above equation i s {AX} and therefore, {AX} can be obtained as, {AX} = [D]" 1 [{AP} + [M] {E} t + [C] {F}t] (A3.14) in which, [»]-[{S.]*WJE+[" tU (A3-15> Now knowing {AX}, the unknowns {x}^ and {x}T and {x}T can be evaluated. {x}T i s simply, {X} T = {AX} + {X} t (A3.16) 189 From equation (A3.2), an expression for {x}T i s , W T • "p^2 [{AX} - At {X} t " (| - 8) A t 2 {X} t] (A3.17) S u b s t i t u t i n g f o r {E} t from (A3.9), the eq u a t i o n (A3.17) can be s i m p l i f i e d as, W T = "pit"2 {AX} - {E} t + {X} t (A3.18) From equations (A3.1) and (A3.2) an e x p r e s s i o n f o r {x}^ , a f t e r rearranging terms i s , {X}T = {X}t + ( I - ) At {X} t + ccAt {X} T (A3.19) Knowing {AX} by solving equation (A3.14), the response at time T can be computed using equations, (A3.18) and (A3.19). In the numerical step by step i n t e g r a t i o n , the following sequence of c a l c u l a t i o n s have to be performed for every time step. 1. I n i t i a l v e l o c i t y {x} t and displacements {x} t are known eith e r from values at the end of the preceding increment or as i n i t i a l conditions of the problem. Based on these values, { E } t , { F } t , and { p C O r r } » a r e e v a l u a t e d u s i n g equations (A3.9), (A3.10) and (5.11). 2. With these values and the known non-linear properties of the s o i l deposit, the damping matrix [c] and [ K t ] t are evaluated according to appropriate equations i n Chapter 5. 190 3. The matrix [D] i s then calculated using equation (A3.15). 4. Using the increment i n base acceleration value at time t, i t i s possible to evaluate the r i g h t hand side of the equation (A3. 12). 5. The equation (A3.12) can now be solved for {AX} and then displacement, a c c e l e r a t i o n and v e l o c i t y vectors can be evaluated from equations (A3.16), (A3.18) and (A3.19) r e s p e c t i v e l y . When step 5 has been completed, the analysis f or t h i s time increment i s f i n i s h e d and the entire process may be repeated for the next time step. Obviously t h i s process can be c a r r i e d out consecutively for any desired number of time increments; thus the complete response h i s t o r y can be computed. Two important aspects have to be considered i n any numerical i n t e g r a t i o n procedure. They are the accuracy and s t a b i l i t y of the procedures. Accuracy ref e r s to how well the numerical s o l u t i o n matches the exact continuous s o l u t i o n . S t a b i l i t y r e f e r s to whether extraneous solutions are introduced i n such a way that they increase rather than decay, and thus come to dominate the r e s u l t s . Usually there i s an upper l i m i t to At that i s necessary to guarantee s t a b i l i t y , and the value of that l i m i t depends on the type of element s t i f f n e s s and mass matrix as well as on « and 8. With a l i n e a r a c c e l e r a t i o n assumption (<= = 1/2, 6 = 1/6) the analysis w i l l give good accuracy i f the shortest period of the deposit i s 5 to 10 times greater than At (Clough and Penzien 1975). The l i n e a r a c c e l e r a t i o n method i s only c o n d i t i o n a l l y stable, and i t w i l l blow up i f i t i s applied to s o i l structures with the shortest period le s s than about 1.8 times the i n t e g r a t i o n i n t e r v a l . Thus the time increments must be made 191 short r e l a t i v e to the least period of v i b r a t i o n contained i n the system regardless of whether the higher modes contribute s i g n i f i c a n t l y or not. In the analysis using f i n i t e element procedures the shortest periods of v i b r a t i o n i n general may be several orders of magnitude less than the periods associated with the s i g n i f i c a n t response. In these cases, the l i n e a r a c c e l e r a t i o n method cannot be used because of the very short time step required to avoid i n s t a b i l i t y ; instead, an unconditionally stable method i s required which w i l l not blow up regardless of the time step. Several unconditionally stable methods are a v a i l a b l e . The constant average a c c e l e r a t i o n method (« = 1/2 , 8 = 1/4) i s one of the simplest of these methods. But t h i s assumption has been reported not to give good r e s u l t s than the methods with l i n e a r a c c e l e r a t i o n assumption. The Wilson 9-method (Wilson, et a l . 1973), i s also an unconditionally stable method. This method i s a modification of l i n e a r a c c e l e r a t i o n method and i s reported to be the best of a l l unconditionally stable methods (Clough, et a l . 1975). The Wilson 8-method i s based on the assumption that the acceleration varies l i n e a r l y over an extended computation i n t e r v a l , %, such that, x = 9 At where 0 > 1.37 (A3.20) When 0 = 1 , t h i s method reverts to the standard l i n e a r a c c e l e r a t i o n method. The analysis procedure i s exactly the same as the procedure presented above except that i n the equations the time step At has to be r e p l a c e d by T and a l s o the equations to e v a l u a t e {x}-j,, {&} ,^ and { x } T have to be modified. 192 Since t h i s i s e s s e n t i a l l y a l i n e a r a c c e l e r a t i o n method, here = 1/2 and 8 = 1 / 6 . By inspection, the required equations can be rewritten as follows. The equation (A3.12) can be rewritten as, f^M] + 1[C] + [ K f c ] t ] { A X } m { ~ p } + [ M ] [ c ] { p } t ( A 3 > 2 1 ) i n which, { E*t = x ^ t + 3 M t ( A 3 ' 2 2 ) and, {F} t = 3 {X} t + \ {X} t (A3.23) After evaluating {AX}, which i s over an extended time increment the displacement, v e l o c i t y and a c c e l e r a t i o n values should be computed at time t = T. The a c c e l e r a t i o n at t = T can be evaluated using, W T = [*2 {AX} + { E } J {X} t (A3.24) {x},p, from equation (A3.19) with « = 1/2 and 8 = 1/6 i s , {X} T = W t + f [ W t + {X} T] (A3.25) 193 and also the displacement {x}^, from (A3.2) i s , {X} T = ( X } t + At ( X } t + *f - [ {x } T + 2 { X } J (A3.26) f I t must be remembered that s t a b i l i t y i n numerical i n t e g r a t i o n does not guarantee accuracy or vice versa. The Wilson 0-method imposes a r t i f i c i a l damping i n higher modes. But knowing that the response due to higher modes of v i b r a t i o n contributes very l i t t l e to the true response of structures, this method i n a way f i l t e r s out the high frequency response. Therefore, t h i s method has been found to y i e l d r e a l i s t i c r e s u l t s i n a number of dynamic analyses.