EARTHQUAKEINDUCED SEPARATION
IN A COMPOSITE DAM
by
Fatili 6nci.il
B.Sc., Middle East Technical University, 1992
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment of
the requirements for the degree of
Master of Science
Civil Engineering
1995
This thesis for the Master of Science
degree by
Fatih Onciil
has been approved
by
Date
Dedicated
to my parents
Onciil, Fatih (M.S., Civil Engineering)
EarthquakeInduced Separation in a Composite Dam
Thesis directed by Professor N. Y. Chang
ABSTRACT
This thesis discusses the evaluation of possible seperation at the interface of a com
posite dam under earthquakes. Analysis is performed using a nonlinear dynamic
analysis code NIKE2D with RambergOsgood elastoplastic material model and
interface algorithm that allows separation and frictional sliding.
The Left Wing Dam of the Folsom Dam is analyzed because it is very well doc
umented. The concrete gravity main dam is partially embedded in the Left Wing
envelopment fill. The Folsom Darn is located on the American River 20 miles north
east of Sacremento, California.. Although it is not located near an active causative
fault capable of producing a strong motion earthquake, the chance exists that it
may be subjected strong shaking from an active distant fault. The earthquake
induced separation along the earthconcrete interface is a main potential problem
in the Wing Dams.
A crosssection in the Left Wing Dam was chosen in this study. NIKE2D was
used to analyze its behavior under Koyna Earthquake of India. The NIKE2D
iv
results are compared with the results of dynamic analysis using FLUSH. The
NIKE2D result confirms the suspicion of the interface separation in the FLUSH
analysis
The analysis result indicates: 1.) the soilconcrete interface of the Left Wing
Dam can separate under a strong motion earthquake, 2.) the separation can be
repetitive, 3.) the separation can reach a significant depth.
Both separation depth and maximum separation are related to maximum base
acceleration of the input ground motions used in the analyses.
This abstract accurately represents the content of the candidates thesis. I rec
ommend its publication.
Signed
NienYin CHANG.
CONTENTS
1. Introduction 1
1.1 Problem Statement................................................ 1
1.2 Significance of Research ........................................ 2
1.3 Research Objectives.............................................. 3
2 . Engineering Seismicity and Design Earthquake Selection 4
2.1 Introduction..................................................... 4
2.2 Determining Ground Motion........................................ 5
2.2.1 Geologic and Seismologic Considerations................. 5
2.2.2 Determining Design Ground Motion Parameters............. 9
2.2.3 Duration of the Earthquake............................... 10
2.2.4 Predominant period of rock acceleration................. 12
2.2.5 Maximum amplitude of acceleration....................... 13
2.3 Conclusion...................................................... 13
3 . Literature Review On Dynamic Analysis Methods of Earth Dams 15
3.1 Introduction.................................................... 15
3.2 Pseudostatic Stability Analysis ................................ 15
3.3 Newmarks Method................................................ 17
3.3.1 Determination of Yield Points............................ 19
3.3.2 Calculation of Displacement.............................. 23
3.4 Seed and Makdisi Simplified Method.............................. 26
3.4.1 Determination of the Yield Acceleration.................. 27
vi
3.4.2 Determination of Earthquake Induced Acceleration.......... 29
3.4.3 Calculation of Permanent Deformation...................... 32
3.5 Strain Potential Method.......................................... 36
3.6 Summary.......................................................... 37
4 . Folsom Dam 39
4.1 Introduction................................................... 39
4.1.1 General................................................... 39
4.1.2 Description of the Left Wing Dam, Interface Zone and Con
crete Gravity Dam................................................. 40
4.2 Geologic Conditions, Seismicity and Geotechnical Properties of Soils 41
4.2.1 Site Geology.............................................. 41
4.2.2 Seismological and Geological Investigations............... 41
4.2.3 Selection of Design Ground Motions........................ 42
t
4.3 Engineering Properties of Soil at the Left Wing Dam.............. 44
5 . NIKE2D A Nonlinear Dynamic Analysis Code 50
5.1 Introduction.................................................... 50
5.2 General Review................................................... 51
5.3 Material Models................................................. 52
5.4 Finite Element Equations......................................... 54
5.5 Quasistatic Anafysis ............................................ 55
5.6 Dynamic Analysis................................................ 57
5.7 Slidelines....................................................... 58
5.7.1 Penalty Formulation ...................................... 59
6. Previous Dynamic Analyses 63
6.1 Introduction.................................................... 63
vii
6.1.1 Description of Dynamic Analysis Code FLUSH .......... 63
6.2 Analysis Procedure to Estimate Possible Separation........... 64
6.2.1 Interface with HINGE Support.......................... 67
6.2.2 Interface with ROLLER Support......................... 67
6.3 Conclusion................................................... 73
7 . Dynamic Analysis with NIKE2D 74
7.1 Introduction................................................. 74
7.2 Input to NIKE2D............................................ 74
7.2.1 Finite Element Mesh................................... 75
7.2.2 Determination of Material Properties................. 76
7.2.3 Interface Contact Conditions.......................... 82
7.2.4 Ground Motion ..................................... 83
7.3 Results of Dynamic Analysis.................................. 84
7.4 Conclusion................................................... 87
8 . Summary, Conclusions and Recommendations 92
8.1 Brief Summary . ...................................... 92
8.2 Conclusion................................................... 93
8.3 Recommendations for Future Study .......................... 94
Appendix 95
A Plans and Details of Folsom Dam 95
Bibliography
via
FIGURES
Figure
1.1 Illustration of the separation problem in a composite dam........... 2
2.1 Procedure to determine ground motion at a site. (Boggs et. al. 1972) 6
2.2 Duration of strong shaking.(Boggs et. al, 1972) ................... 12
2.3 Predominant period for maximum acceleration.(Seed et. al, 1969) . 13
2.4 Variation of maximum acceleration with earthquake magnitude and
distance from causative fault.(Seed et. al, 1969).................. 14
3.1 Conventional Method for computing elTect of earthquake on stability
of a slope (after Terzaghi, 1950).................................. 16
3.2 Different sliding surfaces (after Newmark, 1965)................... 18
t
3.3 Rigid block on a moving support (after Newmark, 1965) ............. 23
3.4 Rectangular block acceleration pulse (after Newmark, 1965) .... 24
3.5 Velocity response to rectangular block acceleration (after Newmark,
1965).............................................................. 24
3.6 Determination of dynamic yield strength (Seed et. al,1978)......... 28
3.7 Calculation of avarage acceleration from finite element response
analysis (Seed at. a.1, 1978)...................................... 29
3.S Time Histories of Avarage Accelerations For Various Depth of Po
tential Sliding Mass (after Makdisi and Seed, 1978)..................... 30
3.9 Variation of effective peak acceleration, A~lnax, with depth of base of
potential slide mass (after Makdisi and Seed, 1978)................ 31
IX
3.10 Variation of permanent displacement with yield acceleration Sum
mary of all data (after Makdisi and Seed, 1978)................. 34
3.11 Variation of avarage normalized displacement with yield acceleration
(after Makdisi and Seed, 197S).................................... 35
4.1 Laboratoiy tests of variation of shear modulus with shear strain. . 45
4.2 Shear wave velocity distribution at the upstream site of Left Wing
Dam.................................................................. 46
4.3 Variation of shear wave velocity with confining pressure and degree
of saturation, (after Hardin, 1963).................................. 47
4.4 Material properties of the Left Wing Dam at ST 299+35. (Chang
et. al 1992)......................................................... 48
5.1 Typical Interface.(Hallquist, 1978).................................. 58
5.2 A method of imposing the prescribed displacement. (Cook et. al,
1989).............................................................. 59
5.3 Contact of node m with segment of jk. (Engelmann, 1991).............. 60
6.1 Schematic peocedure followed in the FLUSH analysis................... 65
6.2 Koyna Darn.Earthquake Record, 1967................................... 66
6.3 The distribution of horizontal normal stress along the boundary, at
the end of the static analysis. (Chang et. al, 1992) ............. 68
6.4 Maximum and minumum horizontal normal stresses at the hinge
boundary elements at ST. 299+35. Chang et. al, 1992. (Tension
negative)............................................................ 69
6.5 Superimposed horizontal normal stress in the hinge boundary ele
ments. Chang et. al, 1992. (Tension negative)........................... 70
x
6.6 Maximum and minumum horizontal normal stresses in the roller
boundary elements. Chang et. al, 1992. (Tension negative)........ 71
6.7 Superimposed horizontal normal sresses in the roller boundary ele
ments. Chang et. al, 1992. (Tension negative)................... 72
7.1 Finite Element Mesh of Left Wing Dam at ST. 299+35............... 75
7.2 Material properties of Left Wing Dam at ST. 299+35............... 77
7.3 Variation of Damping Ratio with Shear Strain for Sands (Seed at
al, 1984)...................>.................................... 78
7.4 Typical loading and unloading curve for RambergOsgood Model. . 79
7.5 Original and scaled 1967 Koyna Dam Earthquake motion.............83
7.6 Horizontal displacement time histories of nodes on the soil side at
the interface, (with original Koyna Dam Earthquake Record) .... 84
7.7 Horizontal displacement time histories of nodes on the soil side at
the interface, (with scaled Koyna Dam Earthquake Record by 0.75) 85
7.8 Horizontal displacement time histories of nodes on the soil side at
the interface, (with scaled Koyna Dam Earthquake Record by 0.5) . 86
7.9 Maximum displacement versus depth................................ 88
7.10 Deformed shape of Left Wing Dam at ST. 299+35 at the time of
maximum displacement (displacements are scaled 50 times).........89
7.11 Maximum separation versus maximum horizontal base acceleration. 89
7.12 Depth of separation versus maximum horizontal base acceleration. . 90
A.l Location of Folsom Dam and Reservoir.............................. 96
A.2 Layout of Folsom Dam and Reservoir amid surface geology and lo
cations and borrow sources...................................... 97
xi
A.3 Plan view of Right and Left Wing Dams and detailed plan view of
interface areas.................................................... 98
A.4 Left Wing Dam Typical Section....................................... 99
A.5 Plan of Concrete Gravity Dam........................................100
A.6 Geologic Map, parts of the Folsom Dam...............................101
A.7 Identification of study area (from USGS 7.5 and 15 ft topographic
maps after Tierra Engineering Consultants (1983)................... 102
A.8 Regional lineament map (after Tierra Engineering Consultants (1983)103
A.9 Acceleration histories recommended by Bolt and Seed (1983) .... 104
xn
TABLES
Table
2.1 Comparision of the Richter Scale Magnitude with the Modified Mer
calli Scale, (after Das, 1993)............................... 11
2.2 Duration of fault break for different magnitude of earthquakes, (af
ter Das, 1993) ............................................. 11
7.1 Comparison of NIKE2D and FLUSH analyses...................... 90
XIII
ACKNOWLEDGMENT
I would like first of all to express my gratitude to those who have contributed to
the development of this thesis. I am greatly indebted to Dr. N. Y. Chang for his
guidence and valuable comments on each part of the thesis.
I also wish to thank the members of the committee, Dr. John Mays, and Dr.
Dunja Peric for their helpful comments and suggestions.
The financial support of the Dumlupinar University in Kutahya, Turkey is
gratefully. acknowledged.
Finally, I am thankful to my parents and friends for their love and patience.
xiv
1 Introduction
1.1 Problem Statement
Observations of the behavior of earth dams in past earthquakes have showed us
the importance of seismic stability analysis for evaluating their performance during
strong motion earthquakes. Though many dynamic analysis methods have been
developed in the past two and a half decades, the complicated soilconcrete interface
behavior along the composite section of a soilconcrete dam has never received
significant attention.
In evaluation of the seismic stability of a composite dam, the main problem is
the dynamic interaction between concrete gravity dam and soil embankment. What
could happen to the interface during an earthquake ? Could debonding/separation
take place ? How much is the permanent deformation along the interface or will
there be a permanent gap after the earthquake ? What is the effect of excess
pore pressure on the seismic behavior of the interface ? Could erosion take place
during the pumping action ? Unfortunately, all these questions have been waiting
for answers. Due to the lack of necessary analysis tools it has been really difficult
to perform the rational analysis of the interface performance. Fig. 1.1 helps to
visualize the potential separation in a composite dam.
Nonlinear finite element computer code, NIKE2D, developed at the Lawrence
Livermore National Laboratory was made available to the University of Colorado
at Denver with a signed agreement. In this thesis, NIKE2D was used to study
the potential separation of the interface at the Left Wing Dam of Californias
Folsom Dam which was chosen for the interface study. The NIKE2D code incorpo
1
Possible seperation zone !1
Figure 1.1: Illustration of the separation problem in a composite dam.
rates nonlinear and implicit finite element analysis with multiple nonlinear material
models and interface algorithm that allows the interface debonding, rebonding, and
frictional sliding.
1.2 Significance of Research
Lysmer in 197S commented on the problems associated the dynamic analyses of
earth structures and soilstructure interaction :
Having taken a sober look at the possibilities for nonlinear analysis and with
this the dream of all soil engineers, dynamic effective stress analysis with evaluation
of dynamic pore pressures, you might conclude with the writer that thestateof
theart of dynamic analysis is still in its infancy and that for a long time to come
most soil dynamic problems will have to be solved by total stress analysis and
essentially linear methods of analysis.
Has the soil engineers dream come true ? The answer is not quite. There is
2
no computer code which can realistically simulate the nonlinear dynamic behavior
of an earth darn with the pore pressure generation and appropriate interface algo
rithm. Many people are still using linear and/or equivalent linear elastic models.
The present research can be considered as an important step forward in re
alizing the dream of the realistic simulation of the soilstructure behavior under
earthquakes. In this thesis, NIKE2D, a nonlinear, implicit, twodimensional fi
nite element code, is used. Of the many nonlinear models available in the code,
RambergOsgood elastoplasticity is used in this research to model soil behavior.
For interface simulation, the slideline algorithm which allows separation with fric
tional sliding is used. Even though RambergOsgood is not a perfect model for
soils and though is not possible to evaluate the pore pressure generation with
NIKE2D, this research is on the cutting edge of the todays technology and has a
big importance in soilstructure engineering.
1.3 Research Objectives
Since NIKE2D is a powerful dynamic analysis code, it is possible to perform many
different research with or without interface. Due to the high nonlinearity of the
problem investigated in this thesis, the analysis is quite time consuming. Thus,
the scope and objectives of this study is limited to :
evaluation of separation potential at the interface,
depth of separation along the interface,
effect of different input ground motions on the behavior of the interface.
The analysis procedure adopted in this study can be expanded into many dif
ferent types of soilstructure interaction problems.
3
2 Engineering Seismicity and Design Earthquake Selection
2.1 Introduction
Before starting any dynamic analysis of an engineering structure, it is necessary to
determine the seismicity of the area and the probable, site specific ground motion.
One should know what kind of ground motions used in the analysis so that results
can be interpreted correctly. Krinitzsky (19S3) describes the general principals of
selecting design ground motion as: The selection of earthquake motion is depen
dent on the engineering analysis to be performed. Each type of analysis require
different input motion. For example for pseudostatic analysis it is enough to deter
mine the coefficient of earthquake from maps prepared for this purpose. Whereas
for dynamic analysis it is required to find the acceleration time history specific
to a site. Determination of an acceleration time history of a site is not straight
forward and needs empirical knowledge, theoretical interpretation and professional
judgment.
Each local area has different special characteristics in its geology, seismic his
tory, attenuation, recurrence, etc. Therefore, a different decision level is required
in each part of the analysis. Since there is no unique and standard procedure it
is useful to use different approaches and compare them as much as possible.
In this chapter, the methods to determine the design earthquake will be ex
plained briefly. The importance of the methods will be discussed in terms of
determining the earthquake design parameters for an earth dam.
4
2.2 Determining Ground Motion
In order to determine the design ground motion for a dam site, a methodical
approach should be followed. A sample procedure is shown in Fig. 2.1. It is really
expensive and time consuming to follow a procedure like the one in Fig. 2.1, but for
huge structures like dams extra attention should be given to the selection of design
ground motion. As a first step, the magnitude and epicenter of the earthquake
that is expected to affect the dam site should be determined. The ground motion
at the site resulting from the estimated earthquake can then be determined using
various methods. The following sections explain how to determine the design
ground motion of a specific site.
2.2.1 Geologic and Seismologic Considerations
Motion that are specified should based on the following relationships:
The presence of active faults,
The magnitude and tjfe recurrence time oFthe earthquake,
The distance, type, orientation and the character of the surface displacement
of potential earthquake sources,
The peak motions (particle velocity, acceleration, displacement), as well as
duration and.predominant period that are associated with these events,
The attenuation characteristic of these earthquakes depending on the geology
of the site,
The overall effects of the site characteristics such as soil, rock, topography,
?
field conditions, etc... on the resultant motions.
5
Figure 2.1: Procedure to determine ground motion at a site. (Boggs et. al. 1972)
6
All of the items listed above can be determined from ;
1. the historical seismicity record,
2. the seismographic record,
3. the geologic study of the site.
When there is a large hazard to life, i.e. a dam above an urbanized area, care
should be taken in using the information sources given above, and the values used
should approximate the worst case.
After obtaining all the information above, dynamic analysis of an earth dam is
performed to evaluate:
maximum stresses induced by the ground motion,
the material strength to resist these maximum stresses,
pore pressure generation,
liquefaction potential,
post earthquake stability,
possible cracks or gaps at the interfaces for the composite dams,
excessive settlements which can cause overtopping of the water,
permanent deformations
possible movement at the foundation level.
7
This thesis doesnt address all the items above. As already mentioned, in this
study only the possible separation between concrete and soil, which is the case for
composite dams, is investigated. In the following sections, the earthquake infor
mation sources and methods determining peak motions are discussed.
1. ) Historical Seismicity:
An important first step is to gather all the information about the previous earth
quakes that have occurred at or near the site from newspapers, diaries, early sci
entific and historical works, and manuscripts written about a specific earthquake.
The information obtained should be examined critically, tabulated, and plotted
geographically along with the geology. It should be noted that: no earthquake
should be related to a fault unless there is an evidence that the fault actually
moved during that earthquake [14]. Since fault movement is not the only source
of earthquake the correct source should be determined very carefully.
This information may be very brief or may be extended over the past years. One
should be very carefull when extrapolating the magnitude of earthquakes beyond
their data base.
2. ) Seismograph Record:
The instrumental record of an earthquake is a very useful tool in a seismic hazard
evaluation. These records provide information about size, location, depth, time of
occurrence of earthquakes. It is not possible to determine this information from
the felt ground motions and damage patterns.
The instruments to record the ground motions have been available since the
1900 but have not been placed at every susceptible location around the world yet.
As more instruments are placed in more areas of the world, the usefulness of the
seismographic record in assessing the seismic hazard will increase.
8
3.) Geologic Study:
Geologic study includes review of available literature, especially in terms of struc
tural and tectonic history, the use of trenching, borehole, agedating, geophysical
techniques, the interpretations of satellite imagery, and plate tectonics. The pur
pose of air photo interpretation is to judge whether the faults are active or inactive.
Once a fault is said to be active, it is necessary to determine the whether it is ca
pable of generating earthquakes.
The larger the capable fault, the greater the potential earthquake. Thus, re
lationships have been developed between dimensions of faults and magnitude of
earthquakes. Dimensions include length of fault, displacement during the move
ment of the fault, whether the movement is on a primary fault or a branch fault or
an accessory fault. A useful summary relating faults to earthquake magnitude is
provided by Slemmons (1977). Use of the data is definitely a matter of responsible
judgment on the part of the investigator [14].
The seismic history, in combination with the geology and the evidence of fault
activity, can be combinedto identify earthquake zones.
2.2.2 Determining Design Ground Motion Parameters
Ground motion characteristics show significant differences depending on whether
the site is soil or rock. Because the area of interest of this thesis is a rock site,
soil site ground motions are not discussed here. The important characteristics of
earthquake induced ground motion of rock or rocklike materials are:
1. duration of the earthquake,
2. predominant period of acceleration,
3. maximum amplitude of acceleration.
9
In order to determine the above parameters, the magnitude of the earthquake
at the causative fault should first be estimated. Then, knowing the distance to
the center of the earthquake and the geology of the site or in other words the
attenuation characteristics, basic parameters of ground motion can be determined.
There are many methods to determine the magnitude of an earthquake.. Factors
affecting the magnitude can be listed as follows:
rupture length,
rupture area,
total displacement,
slip rate,
seismic moment,
historical seismicity.
Many different methods to evaluate each factor listed above. All available
methods and formulas are well documented in reference [15]. The use of methods
is totally dependant on the availability of the data and engineering experience as
well as the professional judgment.
On the other hand, besides the magnitude, strong motion parameters can also
be determined with regard to the intensity of the earthquake not the magnitude.
Table 2.1 can be used to convert intensity to magnitude (according to Richter
scale).
2.2.3 Duration of the Earthquake
Shaking is likely to continue at least as long as fault rupture is occurring. Housner
(1965) estimated the rate propagation of fault rupture as 2 mi/s. Based on this,
10
Richter Scale magnitude M Maximum intensity Modified Mercalli Scale
1
2 Ml
3 III
4 IV
5 VI. VII
6 VIII
7 IX, X
8 XI
Table 2.1: Comparision of the Richter Scale Magnitude with the Modified Mercalli
Scale, (after Das, 1993)
Magnitude of earthquake Duration of fault break
(Richter Scale) (s)
5 5
6 15
7 2530
Table 2.2: Duration of fault break for different magnitude of earthquakes, (after
Das, 1993)
11
50
40
 30
CO
o
 20 O
UJ
tn
 10
 0
I 23456789
MAGNITUDE
Figure 2.2: Duration of strong shaking.(Boggs et. al, 1972)
Table 2.2 can be used to estimate the duration. Also the graph proposed by Boggs
(1972) is in good agreement with Housners results and is shown in Fig. 2.2.
The other method for determining duration is called bracketed duration. It is
the inclusive time in which the acceleration level equals or exceeds the amplitude
of 0.05<7, or 0.10
2.2.4 Predominant period of rock acceleration
Predominant period is the period for which the spectral acceleration is at its max
imum. Based on the studies done by Gutenberg and Richter for earthquakes with
magnitudes ranging 5.5 and 6.5 and Figueroa for magnitudes greater than 7, it
was noted that the predominant period, for any epicentral distance, increases with
the magnitude of the earthquake [IS]. By making interpolation and extrapola
tion of the available data, predominant period is plotted versus distance from the
causative fault for different magnitudes of earthquake. This plot is shown in Fig.
2.3.
7^
/
o o
oo
0.0
/
/,
X
/
7^
/
/
12
1.2
0.2
0 I
0 80 160 240 320
Distance from causative fault (km)
Figure 2.3: Predominant period for maximum acceleration.(Seed et. al, 1969)
2.2.5 Maximum amplitude of acceleration
Maximum acceleration depends on the magnitude and the distance from the causative
fault. In order to determine the maximum acceleration for a specific site Fig. 2.4
by Seed et. al (1969) can be used.
2.3 Conclusion
Determining a design ground motion for a specific site entails a lot of research. As
already mentioned, there is no unique method. For critical structures like earth
dams, tall buildings, nuclear power stations, etc... special attention should be
given to every step of ground motion determining procedure. Moreover, different
approaches should be used to be able to make comparison.
13
\
o 40 80 120 160
Dittante from ctumlve faull (km)
Figure 2.4: Variation of maximum acceleration with earthquake magnitude and
distance from causative fault.(Seed et. al, 1969)
14
3 Literature Review On Dynamic Analysis Methods of
Earth Dams
3.1 Introduction
Dynamic analysis methods of earth dams have been improved since many years.
Due to high nonlinearity and complexity of the problem it is very difficult to come
up with a unique method. There are several methods and approaches. There
is no method for composite dams to calculate the separation at the interface in
between concrete and soil. Methods for determining the permanent deformations
are the only available methods which gives an idea about the complex behavior
of earth dams during earthquakes. In this chapter, literature survey of currently
available dynamic analysis methods for earth dams will be presented. Basically
these methods can be classified as; Pseudostatic Stability Analysis and the Methods
to Determine Permanent Deformations.
3.2 Pseiidostatic Stability Analysis
In this method the effects of earthquake on a potential slide mass are represented
by an equivalent static horizontal force determined as the product of a seismic
coefficient, k, or ng and the weight of the potential slide mass, as illustrated in Fig.
3,1. 
Terzaghi in 1950 described the method in the following words;
An earthquake with an acceleration equivalent ng produces a mass force acting
in a horizontal direction of intensity ng per unit of weight of the earth. The
resultant of this mass force, ngW, passes like the weight W, through the center
15
Figure 3.1: Conventional Method for computing effect of earthquake on stability
of a slope (after Terzaghi,1950)
of gravity 0\ of the sliding body abc. It acts at a lever arm with length F and
increases the moment which tends to produce a rotation of the slice abc about the
axis O by naFW. Hence the earthquake reduces the factor of safety of the slope
with respect to sliding from Gs, equation (1) to
' EW + naFW
(3.1)
The numerical value of ng depends on the intensity of the earthquake. Inde
pendent estimates (Freeman, 1932) have led to the following approximate values:
Severe earthquakes, RossiForel scale IX ng = 0.1
Violent, destructive, RossiForel scale X ng = 0.25
Catastrophic ng = 0.5
Equation 3.1 based on the simplifying assumptions that the horizontal accel
eration ngg acts permanently on the slope material and in one direction only.
16
i
Therefore the concept it conveys of earthquake effects on slopes is very inaccurate,
to say the least. Theoretically a value of G's 1 would mean a slide but in reality
a slope may remain stable in spite of G's being smaller than unity and it may fail
at a value of G's > 1, depending on the character of the slopeforming material.
In the above explanation Terzaghi states the inaccuracy of the method and
importance of the slope forming material to evaluate the failure. According to him
failure may occur even if Gs is bigger than unity which is the case in engineering
practice accepted as safe. Terzaghi strongly emphasizes that material property
plays an important role in evaluation of safety. Then it is necessary to make it
clear when and how to use that method in practice. What are the restrictions and
what is the reliability of this method ? According to Idriss et. al (1988) the answer
of this question is: The decision about whether to use the pseudostatic method
for any given dam is dependant on whether the embankment or the foundation
soils are likely to be vulnerable to major strength loss due to postulated level of
shaking. Therefore pseudostatic method can provide reasonable results for clay
or very dense cohesionless soils as a preliminary design. However it shouldnt be
considerate as a final design. It is better to couple pseudostatic method with an
evaluation of permanent deformations.
3.3 Newmarks Method
In 1965 Newmark proposed a method to determine the earthquake induced per
manent deformations of a dam or embankment [8]. In this method three different
possible sliding surfaces are assumed as shown in Fig. 3.2. For each case the dis
placement is considered to be rigidplastic and can be defined as no displacement
until the yield point is reached, after which the displacement may have any value.
In the next sections the determination of the yield point for each of these three
17
Figure 3.2: Different sliding surfaces (after Newmark, 1965)
18
cases and the calculation of the permanent displacement will be discussed.
3.3.1 Determination of Yield Points
1.) Circular Sliding Surface:
A sliding surface is considered in the embankment as shown in Fig. 3.2 where a
circular arc of radius, R, defines the sliding surface. A force, NW, where N and
W correspond to a coefficient of yield acceleration and weight of the sliding mass
respectively, is considered to act along the line shown making an angle a with the
horizantal. According to rigidplastic deformation assumption for accelerations less
than Ng no sliding will occur. Then the dynamic factor of safety is unity when
Ng = Ag, where Ag can be any value of induced acceleration.
In case when Ag = 0, static factor of safety can be defined as the ratio of
resisting moments to the driving moments. Equating the driving forces to resisting
forces:
Wb=RYrds
(3.2)
where
b is the lever arm of the weight.
and the factor of safety can be expressed as:
(3.3)
F ftpSgds Eyls
RY'rds Erds
where
r is shear force acting along the length ds,
sq is the undrained shear strength.
19
An approximate value of N which will cause sliding is calculated by equating
driving and resisting moments as follows:
Wb + NWh = Â£ Â£s
where
h is the lever arm of the force NW
By subtracting equation 3.2 from equation 3.4, one gets:
NWh = R'%2 s9ds .nyvis
(3.5)
Dividing the equation 3.5 by equation 3.2 and multiplying through b/h :
which can be written as
ESqds
Erds
(3.6)
N
E Â£q
 1
h LEr
substituting equation 3.3 in to equation 3.6:
(3.7)
A' = \ [F 1] (3.8)
The minumum value of N occurs when h is maximum or equal to d, where d is
equal to the distance from the center of arc to the center of gravity of the sliding
mass. For this case:
Nmm = 2[F1] = {F1)sW
(3.9)
20
where
/? is the angle between d and the vertical
when N acts horizontal:
Nmin = (F(3.10)
2.) Plane Sliding Surface:
For plane sliding case it is found that the most critical plane is the upper slope,
making an angle 0 with the horizantal. For this case with a material having internal
angle of friction , the factor of safety against sliding becomes:
p = tanM
tan(0)
The minumum value of N can be written as :
(3.11)
Amin = [F l)sin(0)
(3.12)
3.) Block Sliding:
As shown in Fig. 3.2, sliding can take place only on a horizontal plane and between
vertical fissures or embankment surfaces. For static case it is obvious that the shear
stress on the horizontal plane is zero. The only resisting force under earthquake
loading is the undrained shear strength sq.
Equating the driving forces and resisting forces per unit width:
NW = '^2sqds (3.13)
21
where
ds is the length of the element resistances act.
From equation 3.13 it can be concluded that N is ratio of total horizontal
resistance to weight. The pore pressure can be expressed as:
p = 7/1 up (314)
where
up is the pore pressure
7 is the bulk density
The weight per unit width is:
W = ^ 7/ids (3.15)
Undrained shear strength is a function of effective overburden pressure and
for the case of normally consolidated soil the ratio of sq to p is constant. Using
equations 3.13 and 3.14, for normally consolidated soil N can be determined as:
or
Sq E 7^ds XT Upds
P H 7/ids
Â£i T1 H Upds \
p V EtWs;
(3.16)
99
x=y+u
y(t)
Figure 3.3: Rigid block on a moving support (after Newmark, 1965)
N = 7 (1 r) (3.17)
P
where
_ D upds
7 U H 7/ids
(3.18)
The value ?u is not constant and should be calculated as a conservative average
value to take the pore pressure increase during the earthquake.
3.3.2 Calculation of Displacement
The movement of the embankment dam is based on the assumption that the whole
moving mass as a single rigid bodywith resistance mobilized along the sliding
surface. The relative motion of the mass, M, with respect to the ground can be
expressed as (see Fig. 3.3):
U = xy
(3.19)
23
i
Time', t
Figure 3.4: Rectangular block acceleration pulse (after Newmark, 1965)
Time, t
Figure 3.5: Velocity response to rectangular block acceleration (after Newmark,
1965)
24
where
x is the displacement of the mass,M
y is the displacement of the ground
The resistance of the motion is accounted for by a shearing resistance of mag
nitude, NW, which corresponds to an acceleration of the ground of the magnitude
Ng, that would cause the mass relative to the ground.
In Fig. 3.4, the accelerating forces on the mass are shown. The mass, M, is
assumed to be subjected a single pulse of magnitude, Ag, up to time t0. The
resisting accelerating force, Ng is shown by a dashed line. The accelerating force
acts for only a short time but resisting force lasts until the direction of motion
changes. In Fig. 3.5 the velocities are shown for both accelerating and resisting
forces. The maximum velocity for the accelerating force has a magnitude, V, given
by the expression:
V = Agt0 (3.20)
After the time t0 is reached the velocity due to accelerating force remains con
stant, while the velocity due to resisting acceleration has the magnitude Ngt. At
time tm the two velocities become equal and the net velocity becomes zero or in
other words the block comes to rest with respect to the ground. The time tm can
be expressed as:
V
*m~ Ng
(3.21)
The maximum displacement of the mass relative to the ground U is obtained
by computing the shaded triangular area in Fig. 3.5. The calculation is as follows:
25
1 T 1_,
q Vtm ~ ^ Wo
IZ1 IY1
2iVg 2 Ag
(3.22)
The accelerationtime pulse shown in Fig. 3.4 corresponds to an infinite ground
displacement which is not the case in real earthquake motion in which there are
many random positive and negative pulses. The result given in equation 3.22
generally overestimates the relative displacement for an earthquake because it does
not take into account the negative pulses. However, it renders a reasonable order
of magnitude for the relative displacement. Also it indicates that the displacement
is proportional to the square of the maximum ground velocity. In this method the
accelerations in the direction of downward slope causes the mass to move whereas
the accelerations in the direction of upward slope doesnt unless it is extremely
large.
3.4 Seed and Makdisi Simplified Method
In 1978, Makdisi and Seed [6] developed a method to determine the permanent dis
placement in embankments subjected to earthquake ground motion. This method
basically uses the principles defined in the Newmarks block sliding method [8], but
instead of rigid body motion it uses dynamic response of embankment as stated by
Seed and Martin (1966) [9]. It assumes that failure occurs on a well defined slip
surface and that the material behaves elastically at stress levels below failure but
develops a perfectly plastic behavior above yield. General outline of this method
can be stated as follows:
26
I
1. A yield acceleration i.e. an acceleration at which a potential yield surface
would develop a factor of safety of unity is determined. Yield acceleration is a
function of the geometry of embankment, undrained strength of the material
and the location of potential sliding mass.
2. Using finite element techniques together with the straindependant soil prop
erties, earthquake induced accelerations in the embankment are determined.
From these analyses timehistories of average accelerations for various sliding
masses can be determined.
3. For a given potential sliding mass, when the induced acceleration exceeds
the calculated yield acceleration, movements are assumed to occur along
the direction of the failure plane and the magnitude of the displacement is
evaluated by a simple double integration procedure.
The above steps will be discussed in detail together with the design charts.
3.4.1 Determination of the Yield Acceleration
The yield acceleration, k.y can be defined as the acceleration which causes horizontal
force on the sliding mass and produces factor of safety of unity. If the induced
acceleration exceeds ky value sliding mass experiences permanent displacement.
For soils which do not show large cyclic strains or pore pressures and maintain
most of their original strength after earthquake shaking, the value of ky can be
calculated by stability analysis using limiting equilibrium methods. However cyclic
loading condition is not the same as in static loading condition. The soil can be
subjected to stress pulses higher than its static strength and that might cause some
permanent deformations without a complete failure. Therefore, when the material
subjected to cyclic loading induced by earthquake shaking it is necessary to define
27
Figure 3.6: Determination of dynamic yield strength (Seed et. al,1978)
a yield stress below which the material undergoes nearly elastic deformations and
above which the material exhibits plastic deformation of magnitudes dependant on
the number and the frequency of the pulses applied.
Fig. 3.6 illustrates the concept of cyclic yield strength [6]. The material in this
case has a yield strength of about 90 % of its static undrained strength. In Fig.
3.6 , it can be seen that under 100 cycles of stress, which amounts to 80 % of static
undrained strength, the material behaves in a near elastic manner. However when
10 cycles of stress, which amounts to 95 % of static undrained strength, is applied,
substantially large permanent deformation is observed. Therefore it is reasonable
to say that this material has the cyclic yield strength which is 90 % of its static
undrained strength.
After determining the cyclic yield stress of the slope forming material from the
lab tests, the next step is to find the yield acceleration for a specific slip surface
by using one of the available methods of stability analysis.
Having determined the yield acceleration the next step is to find the time history
of earthquake induced average accelerations of that potential slip surface.
28
Element
n
F (t) = Jfrhv. (t )l4 + C^( t) di)
il
ne number of elements along the sliding surface
k =F(t)/W
av
Figure 3.7: Calculation of avarage acceleration from finite element response anal
ysis (Seed at. al, 197S)
3.4.2 Determination of Earthquake Induced Acceleration
In order to find the permanent deformation it is necessaiy to determine when the
average acceleration of the slip surface exceeds the yield acceleration. The steps
to find the average acceleration history can be summarized as follows:
1. Find the stress histories of each element on the slip surface by using finite
element computer codes using elastoplastic material model or equivalent
linear straindependent properties (e.g. Nike2d, QUAD4).
2. At each time step find the resultant of forces along the slip surface. This
procedure is illustrated in Fig. 3.7.
I
\
29
Figure 3.8: Time Histories of Avarage Accelerations For Various Depth of Potential
Sliding Mass (after Makdisi and Seed, 1978)
3. To find the average acceleration kav, divide the resultant force by the mass
of the sliding block.
Repeat the above procedure for every time step to get the time history of
average acceleration for a particular sliding mass.
As an example, for a 150fthigh dam subjected to 30 seconds of the Taft earth
quake record scaled to produce a maximum base acceleration 0.2g, the variation
of the time history of kav with the depth of sliding mass within the embankment
together with the time history of crest accelerations is shown in Fig. 3.8, [6].
Comparing the acceleration timehistory of the crest and the average accelera
tion time histories for sliding masses at different depths, the similarity in frequency
content can be seen easily. The amplitudes decrease as the depth of sliding mass
30
0
Figure 3.9: Variation of effective peak acceleration, A'max, with depth of base of
potential slide mass (after Makdisi and Seed, 197S)
increases towards the base of the embankment.Lets call the maximum crest ac
celeration as ii,nax, aiid the maximum average acceleration for a potential sliding
mass at a specified depth, y as A:lliax.
If a relationship can be constructed between the depth and fcmax/iimax for a
range of embankments and earthquake loading conditions, then it would be enough
to find umax to determine kmax at different depths of sliding masses. A simplified
procedure to find u,ax is explained in [7].
31
Makdisi and Seed [6] collected the previously published data for different em
bankment heights and different earthquake loading conditions including the finite
element results [6], and established a range of curves for design purpose. Those
curves are illustrated in Fig. 3.9 and considered accurate enough for practical pur
poses. If a conservative estimate of the acceleration is desired the upper bound
curve shown in Fig. 3.9 may be used leading to values which are 10 to 30 % higher
than those estimated using the average relationships.
3.4.3 Calculation of Permanent Deformation
Once the yield acceleration and the timehistory of average induced acceleration
for a potential sliding mass have been determined, the permanent displacements
can be calculated.
By assuming a direction of sliding surface and writing the equation of motion
for the sliding mass along that direction the displacements can be calculated using
numerical techniques. In this method the sliding surface is assumed to be hori
zontal which is not the case in real life. However studies made for other directions
of the sliding surface showed that this factor has little effect on the computed
displacements [6].
Makdisi and Seed (1978) analyzed many cases to determine an order of mag
nitude of the deformations induced in embankments due to strong shaking. They
used 3 different earthquake magnitudes : 61/2, 71/2, 81/4. For each earthquake,
a few embankments were used having different heights (ranging from 75150 ft),
different slopes and different material properties. For each embankment the fol
lowings were determined:
acceleration time history of the crest
average acceleration time history for a potential sliding mass extending through
32
almost the full height of the embankment
first natural period of the embankment, T0
ky / ^max for the crest and the full height of the embankment
permanent displacements at the crest and for the full height of the embank
ment (Makdisi and Seed (1978) showed that displacements change linearly
from top to the bottom of the embankment. Therefore it is enough to calcu
late the maximum (at the top), and the minumum (at the bottom) displace
ments).
The calculated displacements are plotted versus A'y/&max ratio for each different
magnitude of earthquakes in Fig. 3.10. According to Fig. 3.10 for low values of
ky/hnaN (< 0.1) the displacement values are extremely high. This doesnt show the
real field behavior of embankments. Therefore the portion of the graph which is
ky/km!ix smaller than 0.1 is not used for the design purposes. Moreover a large scat
ter in the computed results of displacements is observed. Makdisi and Seed (1978)
ended up with a conclusion that for a specific potential sliding mass with a yield
acceleration, the magnitude of permanent displacement induced by an earthquake
is controlled by:
1. the amplitude of induced average accelerations, which is a function of base
motion, the amplifying characteristics of the embankment and the location
of the sliding mass within the embankment
2. the frequency content of the average acceleration time history, which is gov
erned by the embankment height and stiffness characteristics, and is usually
dominated by first natural frequency of the embankment
33
ail11________i________l
C o 0.2 0.4 0.6 0.8
k / k
y/ max
Figure 3.10: Variation of permanent displacement with yield acceleration Sum
mary of all data (after Makdisi and Seed, 197S)
34
0 0.2 0.4 0.6 0.8 1.0
k A
y/ max
Figure 3.11: Variation of avarage normalized displacement with yield acceleration
(after Makdisi and Seed, 1978)
3. the duration of significant shaking, which is a function of the magnitude of
the specified earthquake
Therefore to reduce the large scatter in the data in Fig. 3.10 the permanent
displacements were normalized with respect to its calculated first natural period,
T0, and with respect to the maximum value of the average acceleration time history,
kmax used in the computation. Because of the scattered data of the new normalized
curves, average curve for each data range is calculated. Average curves for those
normalized curves are shown in Fig. 3.11. Fig. 3.11 is considered adequate for
design purposes to find the permanent displacements. It should be noted that for
yield acceleration values less than 0.2 the displacement values may be unrealistic.
Finally the steps finding the permanent displacements can be summarized as
follows:
1. Determine the maximum crest acceleration, itmax, and the natural period, T0
35
2. Using Fig. 3.9 find the maximum value of average acceleration history, &njax
for any level of specified sliding mass
3. Use Fig. 3.11 together with fcmax and T0 to find the permanent displacement
for any value of yield acceleration associated with that sliding surface.
3.5 Strain Potential Method
The basic steps and principles of this method are outlined by Seed [4]. This
method was developed upon the recognition of the limitations of the pesudostatic
analysis approach and the difficulties of evaluating a yield stress criterion for many
saturated cohesionless soils.
Basic principles of the method can be summarized as follows:
1. Determine the cross section of the dam to be used in the analysis.
2. Determine the design ground motion.
3. Determine the stresses in the dam before the earthquake.
4. Determine the dynamic properties of the soils.
5. Evaluate the response of the embankment to the base acceleration selected
in step 2.
6. Subject representative samples of the embankment materials to the com
bined effects of the initial stresses and superimposed dynamic stresses and
determine their effects in terms of the generation of pore water pressures
and development of strains. Perform a sufficient number of these tests to
permit similar evaluations to be made, b}' interpolation, for all elements in
the embankment.
36
7. From the estimated pore pressure generation and soil deformation charac
teristics an the strength characteristics, calculate the factor of safety against
failure of the embankment during and after the earthquake.
8. If the embankment is safe against failure, calculate the strains induced by
the ground motion to determine the overall deformation of the embankment.
The steps summarized above are the main principles of dynamic analysis of
dams. The most important thing is to use proper judgment and to best tools for
the analysis. Maybe those steps will remain the same but the improvements about
how to perform these steps will go on.
3.6 Summary
History of dynamic analysis of embankments started with Terzaghi (1950) with a
very simplified procedure called Pseudostatic Stability Analysis, later finite element
analysis introduced more sophisticated techniques in the analyses.
Pseudostatic stability anatysis should be performed for small embankments
having materials, like clay or dense sand,' that are not vulnerable to excessive
deformations due to ground shaking. In fact in any case pseudostatic stability
analysis should be considered as the preliminary design not the final design. A
finite element technique must be used in the final design.
Since Seed and Makdisi and strain potential methods use finite element tech
niques, they can be considered as better analysis tools than Newmarks Method.
The accuracy of Seed and Makdisi and Strain Potential method depends on the
computer codes used in the analysis. A good computer code has to have nonlin
ear computing techniques, nonlinear soil models which can simulate the dynamic
behavior of soil together with the pore pressure generation.
37
None of the methods explained in this chapter addresses to the separation prob
lem in Composite Dams, but at least they give some idea about the general seismic
behavior of embankments before going into more sophisticated computational tech
niques to evaluate the separation potential in a soilconcrete dam.
3S
4 Folsom Dam
4.1 Introduction
In this chapter some general information will be presented about Folsom Dam and
especially Left Wing Dam including location, geometry, geology and seismicity
of the area, materials used in the construction and dynamic properties of those
materials. All those information is obtained from references [1, 2].
4.1.1 General
Folsom Dam and Reservoir Project located on the American River in Sacramento,
Placer and El Dorado Counties, California, about 20 airline miles northeast of the
City of Sacramento.
The manmade water retaining structures at the Folsom Dam Project consist
of the concrete gravity dam constructed in the American River channel, the Right
and Left Wing Dams which flank the concrete structure, Mormon Island Auxiliary
Dam which was constructed in the Blue Ravine (an ancient channel of the American
River), and 8 saddle dikes which complete the reservoir rim. A location map and
plan of the project area shown in Figures A.l and A.2.
The Folsom Dam Project was designed and built by the Corps of Engineers in
the period 1948 to 1956, as authorized by the Flood Control Act of 1944 and the
American river Basin Development Act of 1949. Upon completion of the project in
May 1956, ownership of the Folsom Dam and Reservoir was transferred to the US
Bureau of Reclamation for operation and maintenance. As an integral part of the
central Valley project, the Folsom Project provides water supplies for irrigation,
39
domestic, municipal, industrial and power production purposes.
Folsom lake impounds the runoff from 1,875 square miles of rugged mountainous
terrain. The reservoir has a storage capacity of 1 million acreft at gross pool and
is contained by approximately 4.8 miles of manmade water retaining structures
that have a crest elevation of 4S0.5 ft. above the sea level. At gross pool,.elevation
466 ft, there are 14.5 ft of freeboard.
4.1.2 Description of the Left Wing Dam, Interface Zone and Con
crete Gravity Dam
The Wing Dams are zoned embankment dams founded on weathered quartz diorite
granite. A plan of the Wing dams is shown in Figure A.3. The Left Wing Dam
is approximately 2100 ft long and 167 ft high, the core consist of well compacted
decomposed granite and is flanked upstream and downstream by 12ft wide filters.
The Upstream and downstream shells are constructed of gravel, which come from
dredged tailings in the Blue Ravine The filters are the 2 in. fraction of the Blue
Ravine gravel. A plan and a typical section in the vicinit}' of the interface area of
the Left Wing Dam is shown in Figure A.4. In this thesis only the section at ST.
299+35, where concrete monoliths end and earth fill starts, is analyzed.
The Right and Left Wing Dams flank the Concrete Gravity Dam. The Concrete
Gravity Dam consists of 28 50ftwide monoliths founded on hard granodiorite rock.
The overall length of the concrete structure is 1,400 ft, the maximum height is 340
ft measured from the foundation to the crown of the roadway, elevation 480.5
ft, and the crest width is about 32 ft. Monoliths are numbered consecutively (1
through 28) beginning at the right abutment. The plan of Concrete Gravity Dam
is shown in Figure A.5.
Concrete Dam Monoliths 1 through 6 interface with the Right Wing Dam and
are fully to partially embedded in the Right Wing envelopment fill. Monoliths 22
40
through 28 interface with the Left Wing Dam and are partially to fully embedded
in the Left Wing envelopment fill.
4.2 Geologic Conditions, Seismicity and Geotechnical Properties
of Soils
4.2.1 Site Geology
The Folsom Dam and Reservoir Project is located in the low, westernmost foothills
of the Sierra Nevada in central California, at the confluence of the North and South
Forks of the American River.
The important formations at the dam site are: a. quartz diorite granite which
forms the foundation at the Concrete Gravity Dam, Wing Dams, and saddle dikes
1 through 7. Weathered granitic or metamorphic rock is present throughout the
area. Figure A.6 shows a geologic map of the project area
4.2.2 Seismological and Geological Investigations
Detailed geological and seismological investigations in the immediate vicinity of
Folsom Reservoir were performed by Tierra Engineering, incorporated to assess
the potential for earthquakes in the vicinity to estimate the magnitudes these
earthquakes might have, and to assess the potential for ground rupture at any
of the waterretaining structures (see Tierra Engineering Consultants, Inc. 1983,
for comprehensive report). The 12mile wide by 35mile long study area centered
on the Folsom Reservoir was extensively investigated using different techniques.
Figure A.7 shows a closeup of the study area as it surrounds the Folsom Project.
Figure A.8 also shows the major faults in the area. In the investigation of the
faults, shears, and lineaments, it was determined that the closest capable fault
is the East Branch of the Bear Mountains fault zone which has been found to
41
be capable of generating a maximum magnitude M=6.5 earthquake. The return
period for this maximum earthquake is estimated to exceed 400 years (Tierra
Engineering, Inc. 1983). Faults that underlie the water retaining structures at
the Folsom Project were found to be noncapable, so seismic fault displacement in
the foundations of the water retaining structures is judged to be highly unlikely.
The minimum distance between the East Branch of Bear Mountains fault zone
and the Concrete Gravity Dam is 9.5 miles. The focal depth of the earthquake
is estimated to be 6 miles. This hypothetical maximum earthquake would cause
more severe shaking at the project area than earthquakes originating from other
known potential sources.
4.2.3 Selection of Design Ground Motions
As mentioned in the preceding section the fault zone of concern is East Branch
of the Bear Mountains fault zone located at a distance about 15 kilometers from
the site. This fault zone has an extensional tectonic setting and a seismic source
mechanism that is normal dipslip. The slip rate from historic geomorphic and
geological evidence is very small, less than 103 cni/year.
After studying the horizontal ground accelerations recorded on an accelerome
ters normal to the Imperial Valley fault during the ImparialValley earthquake of
1979 and also many additional strong ground motion recordings, Bolt and Seed [3]
recommended the following design ground motions:
Peak horizontal ground acceleration = 0.35 g
Peak horizontal ground velocity = 20 cm/sec
Bracketed duration (>0.05(7) 16 sec
42
Because of the presence of granitic plutons at the site, it is expected that the
earthquake accelerations might be relatively rich in high frequencies. Bolt and Seed
[3] provided 2 accelerograms representative of the design ground motions expected
at the site as a result of a maximum magnitude M equal to 6.5 occurring at the
East Branch of the Bear Mountains fault zone. The accelerograms are designated
as follows:
M6.515K83A
This accelerogram is representative of the 84percentile level of ground mo
tions that could be expected to occur at a rock outcrop as a result of a
Magnitude 61/2 earthquake occurring 15 kms from the site. It has the
following characteristics:
Peak acceleration = 0.35 g
Peak velocity = 25 cm/sec
Duration =16 sec
M6.515K83B ~
This accelerogram is representative of the S4percentile level of ground mo
tions that could be expected to occur at a rock outcrop as a result of a
Magnitude 61/2 earthquake occurring 15 kms from the site. It has the
following characteristics:
Peak acceleration 0.35 g
Peak velocity = 19.5 cm/sec
Duration = 15 sec
43
Figure A.9 shows plots of acceleration as a function for the two design accelero
grams.
4.3 Engineering Properties of Soil at the Left Wing Dam
In this section materials properties of soil used in the Left Wing Dam will be
discussed. Especially the properties which are needed to perform the dynamic
analysis will be explained in detail.
The representative section was chosen at the interface of Left Wing Dam at
station 299+35. The cross section at the interface is shown in Fig. A.4. According
to that figure there are basically 3 different material types: shell, core and transition
zones. Transition zone is very small with respect to other parts. Thats why
transition zone is not considered as a separate material in the dynamic analysis.
Core material is classified as clayey sand and the. shell material is classified as
siltyclayeygravel [1]. It has been shown that clay type of materials can show
very good resistance during the earthquakes. On the other hand sandy soils can
have problems in terms of liquefaction and overall stability. The resistance of
Folsom Dam materials against liquefaction is very well reported in references [1]
and [2] and it is not included in this thesis. In order to reflect the effects of
overburden pressure on the material property, it is required to divide the core and
shell materials into layers. This classification of materials will be discussed later.
For dynamic analysis, the most important material parameters are Gmax values
and the variation of shear modulus G, and damping ratio with respect to shear
strain. Fig. 4.1 shows the laboratory tests results of variation of shear modulus
with shear strain. These results are in good agreement with Seeds (1984) standard
curves for gravelly soil.
From field tests shear wave velocities are determined as shown in Fig. 4.2. This
44
Figure 4.1: Laboratory tests of variation of shear modulus with shear strain.
DEPTH, ft
0
 20
 40
 60
 BO
!00
120
 140
 160
BORING
SETS
CREST EL. 480.5 FT
GRANITE FOUNDATION
480 
460 
440
420 
400 
380
360
340
320
DOWNSTREAM

NOTE: ALL VELOCITIES ARE IN fp9
SCALE IN FEET
30 0 30 60
' nLEj UKBFM
Figure 4.2: Shear wave velocities at the upstream site of Left Wing Dam.
46
ELEVATION, ft
Figure 4.3: Variation of shear wave velocity with confining pressure and degree of
saturation, (after Hardin, 1963)
figure shows the variation of shear wave velocity with depth. Shear wave velocity
can be related to maximum shear modulus as:
Therefore the effect of overburden pressure on the material property can be
reflected with varying shear modulus along the depth of the embankment. Fig.
4.2 shows the velocities at the downstream site of the Left Wing Dam. The only
material property difference between upstream and downstream site of the em
bankment is the degree of saturation. According to Hardin and Richard (1963)
there is not a very big difference in the shear wave velocities for dry and saturated
sands at a specific confining pressure. Therefore, it reasonable to reduce the shear
wave velocity by 50 100ft/s for saturated soil of the same type.
Finally considering all the effects of overburden pressure and saturation and
the test results reported in references [1, 2] material classification can be done
as shown in Fig. 4.4. Actually this classification was done by Chang (1992) to
analyze the separation problem with the computer code FLUSH. In Fig. 4.4, v is
47
Figure 4.4: Material properties of the Left Wing Dam at ST 299+35.
MATERIAL PROPERTIES USED IN OYNAMIC ANALYSIS
MAT'L TYPES Y (ocf) 'Is (ft/s) 900 G (XsO 3145 Gi (Xsf) 2400 DESCRIPTIONS
l Q. 30 137.0 0.065 Shell : coarse dredgetailings, unsaturated
2 0.33 137.0 1100 4700 3000 0.075
3 0.3S 137.0 1250 6065 4000 0.064
4 0.3S 138.0 1000 4190 3000 0.075 Core: decomposed granite, unsaturated
s 0.40 138.0 1300 7190 4000 0.064
5 0.45 142.0 1300 7450 4000 0.064 Core: decomposed granite, saturated
7 0.45 142.0 1700 12745 12000 0.050
8 0.35 140.0 300 2590 2400 0.065 Shell: coarse dredgetailings, saturated
9 0.35 140.0 1000 4040 3000 0.075
10 0.35 140.0 1150 5340 4000 0.064
ST. 299 + 35
poissons ratio, 7 is the unit weight, is the shear wave velocity, G is maximum
shear modulus, G1 is shear modulus for first iteration (special input for FLUSH),
Ai is damping ratio for the first iteration (special input for FLUSH).
49
5 NIKE2D A Nonlinear Dynamic Analysis Code
5.1 Introduction
The following paragraph, taken from Nike2ds users manual, gives a very good
brief description of NIKE2D:
Nike2d is an implicit finite element code for analyzing the finite deformation,
quasistatic and dynamic response of two dimensional, axisymmetric, plain stress
and plain strain solids. The finite element formulation accounts for both material
and geometric nonlinearities. A number of material models are incorporated to
simulate a wide range of material behavior including, elastoplasticity, anisotropy,
creep thermal effects, and rate dependence. Arbitrary contact between indepen
dent bodies is handled by a variety of slideline algorithms. These algorithms model
gaps and sliding along material interfaces, including interface friction and single
surface contact [10].
This program was first developed by Dr. John 0. Hallquist of the Methods
Development Group. Since then it has been extensively used at the Lawrance
Livermore National Laboratory. Nike2d has been improved over the years. It is a
public domain code and open to anybody who wants to study and improve it.
In terms of Geotechnical Engineering it seems to be a promising code to eval
uate nonlinear behavior of soils. Because it already has very powerfull solution
techniques and the capability of handling complicated finite element problems. In
terms of finite element programming it is at the cutting edge of the technology.
It just needs to be improved in the direction of meeting Geotechnical Engineer
ing needs like more realistic soil constitutive models and simulating pore pressure
50
generation.
In the next sections capabilities and some features of Nike2d will be discussed.
5.2 General Review
Nike2d has a wide variety of capabilities that allow the users to solve most of the
twodimensional problems. Nike2ds material models include:
elasticity and plasticity,
anisotropy ,
creep,
thermal effects,
rate dependence.
Element formulation that:
handle geometric nonlinearities,
do not lock for isochoric materials.
A wide variety of boundary conditions and loads including:
prescribed displacements,
body force loads,
concentrated nodal forces and follower forces,
pressure and shear loads,
loads due to thermal expansion.
51
Analysis types are:
quasistatic
dynamic
eigenvalue.
A general interface contact capability including:
sliding with or without separation,
friction,
single surface contact,
tied slidelines for mesh grading,
tie breaking slidelines for simulation of failure,
continuous surface,
one way slide with void
mergedemerge slideline.
5.3 Material Models
The following is a list of currently available material models for plain strain con
dition:
elasticity,
orthotropic elasticity,
elastoplasticity,
52
thermalelostoplasticity
soil and crushable foam
viscoelasticity,
thermalorthotropic elasticity
thermoelasticcreep
BlatzKo rubber
power law plasticity with failure
unified creep plasticity
power law thermoplasticity
strain rate sensitive elastoplasticity
extended two invariant cap
RambergOsgood elastoplasticity
thermal elastoplasticity with 8 curves
three phase thermal elastoplasticity,
strain rate sensitive power law plasticity,
power law thermoplasticity with failure,
nonlinear elastoplasticity
polynomial hyperelastic rubber,
primary, secondary, tertiary creep,
53
deformation mechanism,
GursonTvergaard void growth plasticity,
MooneyRivlin rubber
5.4 Finite Element Equations
The governing equation which is called equation of motion is:
Tij + bi = pui (5.1)
where r is the Cauchy total stress tensor, 6, is the body force per unit volume,
p is the density, U{ are the relative displacements, Cl represents the continuum
domain. The boundary of continuum can be divided into two parts as the boundary
ru where displacements are described and the boundary FT where stresses are
described and the conditions on the boundary can be defined as:
= U{
(5.2)
and
Hj T{
(5.3)
respectively.
The initial conditions are:
/
i(0) = U{0 (5.4)
54
(5.5)
^i(O) Uio
The rate deformation tensor dij is defined in terms of velocity u as
dij + Uj,i) (56)
The Cauchy stress is, in general, a function of rate of deformation dij a set of
history variables H, and the temperature T.
Tij = Tij(dijx H,T) (5.7)
In Nike2d, quadrilateral elements are used for the spatial discritization, yielding
a system of ordinary differential equations,
Mii + Fint(u,u,T) = P(u,b,f,T) (5.8)
Where M is the mass matrix, F is the internal nodal force vector, and P is the
external nodal force vector. P can be a function of nodal displacement u, body
force per unit volume b, time f, and nodal temperature T. In dynamic analysis
the equation 5.8 is solved by Newmark/? method. For quasistatic analysis the
equation 5.8 becomes,
Fint(u, u, T) = P(u, b, f, T) (5.9)
5.5 Quasistatic Analysis
In quasistatic analysis the equation 5.9 is solved incrementally by an iterative
strategy. Without thermal and viscous effects the equation 5.9 becomes:
55
Fint(un+i) = Pn+l (5.10)
where un+i and Pn+i are the nodal displacements and external load evaluation
at time tn+1 At each time step, quantities are known at tn, and the solution
involves finding the displacement un+i that satisfies the equation below:
K, Au = Pn+i F*nt(un+i) (5.11)
This new equilibrium solution is found by iteration. To obtain the solution at
tn+1, the finite element equations 5.9 are first linearized about the configuration at
tn and the displacements at tn+1 are determined. Later the equation 5.11 solved
iteratively.
Convergence is determined by examining the both the displacement and energy
norms. For displacement norm
JIM
^max
< ed
(5.12)
where
Umax is the maximum displacement norm obtained overall of the n steps
is used, and for energy norm;
(Atd)rQ;;+1 <
(Au)<3;+1
(5.13)
is used.
Where ed and ee are tolerances that
are typically 10 2 to 10 3 or
smaller.
56
5.6 Dynamic Analysis
For dynamic analysis the Newmark/? method is used to integrate the semidis
cretized finite element equations 5.1 in time. The Newmark/? family methods is
given by
un+1 = un + Atun + (1/2 /3)At2un + /?At2iin+i (5.14)
un+i = Un + (1 7)+ 7Atun+i (5.15)
In NIKE2D the default values for (3 and 7 are 1/4 and 1/2 respectively. This
choice of parameters represents the trapezoidal rule which is second order accurate
in time, is energy conserving for linear problems and does not introduce numerical
dissipation.
To obtain the dynamic solution at tn+1, the finite element equations 5.1 are first
linearized about the configuration at tn and using the expressions 5.14 and 5.15
the nodal acceleration, velocity, and displacements at tn+1 are obtained. Finally
the equilibrium iterations performed with
K*Au* = Pn+1F*+1)
where
(5.16)
F* = F(U.) + Mui+1+Dui,+1
(5.17)
K' = K + ^M + ^D
(5.18)
where M is the damping matrix. For dynamic analysis, damping can be incor
porated into the model at the element level in two ways: through the specification
57
ktti slave zone
Figure 5.1: Typical Interface.(Hallquist, 1978)
of dissipative material behavior or by a nonlinear adaptation of Rayleigh damping.
Dissipative mechanisms in material behavior are pointwise in nature and include
viscoelasticity and include hysteretic damping caused by cyclic plasticity. Gener
alized Rayleigh damping is more global in nature. However, it is implemented at
the element level to maximize computational efficiency.
5.7 Sidelines . . 
In NIKE2D sliding with no separation or sliding with separation, closure and fric
tion is available. The implementation technique of the slideline algorithms are well
explained in references [11] and [12]. In this section the logic of the slidelines will
be presented briefly.
A typical interface illustration is shown in Fig. 5.1 The nodal points along the
slave line will be called slave nodes, and, likewise, the nodal points on the master
line will be called master nodes. Slave and master segments are straight line
segments joining two adjacent slave and master nodes, respectively. Slave nodes
58
[K] {0} = {R} for active d.o.f.:
Kn Ki
Ki. K + k. Ki* v, = k,v,
Kji Kn K .K l*J
Figure 5.2: A method of imposing the prescribed displacement. (Cook et. al, 1989)
are constrained to slide on the master segments after impact and must remain on
a master segment until a tensile interface force develops. Any separation between
a slave and master line is called void [11].
5.7.1 Penalty Formulation
All of the contact algorithms in Nike2d are based on penalty formulation which
will be described briefly, in penalty formulation, any node that penetrates through
its respective contact surface causes a linear interface spring to be inserted into
the stiffness matrix.
In order to explain the penalty method a simple example taken from the refer
ence [13] will be given. Consider the truss in Fig. 5.2. Let the nodal forces be R\,
R2 and 1?3 and imagine that the vertical displacement at node 1, is to be forced
to have a value v\. This condition can be treated as follows. Let ka be a large
positive stiffness, say 106 times K22. Add a spring of stiffness as shown in Fig.
5.2, and apply the large force kavi in the direction of uj. Load R2 is discarded.
Loads i?! and R3 can continue to act. In this case whatever the values of R\ and
i?3 are, the vertical displacement at node 1 will be constant due to the inserted
59
3
Figure 5.3: Contact of node m with segment of jk. (Engelmann, 1991)
spring having a large stiffness constant. Now solve for all three degree of freedom
(i, ui, U3) in'the usual way. The equilibrium equation with the inserted spring con
stant in the stiffness matrix is shown in Fig. 5.2. The value >i = 0 is permissible.
Mathematically this procedure is called penalty method and ks is called penalty
number. As ks approaches infinity, the constrained iq = v\ is exactly forced. Of
course in the calculations a finite value is given to ks but it shouldnt be very large
to cause numerical problems. Therefore by trial and error method an appropriate
value for penalty number can be found and the corresponding displacements can
be controlled in the desired manner [13].
Now consider a slideline with separation and closure. The following illustrates
the Nike2ds argumentation of stiffness matrix Ks and the internal nodal force
F3 when penetration is detected. Fig. 5.3 shows an isolated portion of the inter
face where node m is penetrating through segment jk [10]. A local equilibrium
60
relationship can be written as
KsAu* = Ps Fs (5.19)
where Au is the incremental displacement vector containing the penalty spring
degreesoffreedom, Ks is the spring stiffness, Fs is the spring internal force and
Ps is the external force arising from internal stress in the interface elements. The
spring degrees of freedom are ordered
Aus = [At)m, Atum., Avj, Awj, Avk, A to*]
The spring stiffness matrix Ks is defined as
(5.20)
Ks = K
force Fs is defined by
s2 sc (1 a)s2
sc c2 (1 a)sc
(1 a)s2 (1 a)sc (1 a)2s2
(1 a)sc (1 a)c2 (1 a)2sc
0 as* asc (1 a)as2
asc ac2 (1 a)asc
(1 a)sc as2 asc
(1 a)c2 asc ac2
(1 a)2sc (1 a)as2 (1 a)asc
(1 a)2c2 (1 a)asc (1 a)ac2
(1 a)asc a2s2 a2sc
(1 a)ac2 a2sc a2c2
sin(0), and k is the penalty stiffness. Th
Fs = k6 s c (1 a)s (1 a)c as ac
(5.21)
(5.22)
61
where 8 is the amount of penetration of node m through segment jk. The
spring stiffness Ks and force F5 are computed for all active slideline nodes and
segments, and are assembled into the global finite element equations. Thus the
stiffness profile changes as analysis with slidelines evolve.
The penalty stiffness k is unique for each segment, and is based on the contact
area and bulk modulus of the penetrated material. If noticeable penetration is
observed, the penalty number should be increased. However high penalty numbers
can be detrimental to the convergence of global iterations. In Nike2d the default
value of k has been chosen to balance global convergence rate with slideline con
straint enforcement on a wide variety of problems. Therefore penalty number k
should not be changed unless it is extremely necassary. Mesh size and bulk mod
ulus, I\, of the penetrated material should be changed first to see if there is any
improvement in the penetration.
62
6 Previous Dynamic Analyses
6.1 Introduction
The dynamic analysis of Folsom Dam was first performed by US Army Corps of
Engineers in 1989, using the dynamic analysis code FLUSH. In this analysis it
was targeted to evaluate the overall seismic stability. The separation problem at
the interfaces was not addressed and recommended as a further study. Detailed
information about this analysis can be found in reference [1]. Investigation of
possible separation at the interface was first studied by Dr. N. Y. Chang and Dr.
H. H. Chiang at the University of Colorado at Denver in 1992, using the computer
code FLUSH. In the next sections only the analysis which addresses the separation
problem will be explained.
6.1.1 Description of Dynamic Analysis Code FLUSH
The computer code FLUSH [20] was developed for analyses of soilstructure inter
action effects during earthquakes. Nonlinear soil behavior is approximated with
an equivalent linear constitutive model which relates shear modulus and damping
ratio to the dynamic strain level developed in the material. FLUSH solves the
equation in the frequency domain and assumes plain strain conditions. It does not
take into account the effect of possible pore water pressure dissipation during the
earthquake. Therefore FLUSH can make only total stress analysis. It is necessary
to input unit weight, shear modulus, and strain dependant modulus degradation
and damping ratio curves for each element in the mesh as the material properties.
FLUSH doesnt have interface properties allowing separation and frictional sliding.
63
6.2 Analysis Procedure to Estimate Possible Separation
Due to the absence of interface properties allowing separation and frictional sliding,
the interface condition can be assumed as hinge or roller connection. Therefore
analyses can be grouped under two different category because of the interface
condition. The main idea is tracing the horizontal normal stress history of elements
at the interface to see if tensile stresses are developed during the earthquake.
Therefore when tensile stress occurs at the interface that indicates the presence
of possible separation. It is called possible separation, because it is really hard
to conclude a certain amount of separation due to very high complexity of the
problem and the lack of appropriate computing tools.
The results of FLUSH analysis will be presented for the cases as listed below:
Case 1: assuming hinge support at the interface
Case 2: assuming roller support at the interface
The classification of material types and representative section is shown in Fig.
4.4. As the strain dependant shear modulus degradation and damping ratio curves
Seeds (19S4) typical curves for sands were used.
The outline of the FLUSH analysis is shown schematically in Fig. 6.1 For
both cases first static analysis is performed to obtain the static stress distribution
at the interface. Then dynamic analysis is done. Using horizontal normal stress
history outputs maximum and minimum stress values are obtained in the boundary
elements. The results from both static and dynamic analyses are then superim
posed to get the resultant maximum and minimum horizontal stress distribution
at the boundary.
In all FLUSH analyses, 1967 Ivoyna Dam Earthquake Record was used as the
64
ure 6.1: Schematic procedure followed in the FLUSH ana
65
Koyna Oam Earthquake Record, M=6.5 at 1.9 miles
Max. Acc. > 0.87 (g)
Figure 6.2: Koyna Dam Earthquake Record, 1967.
66
input ground motion. Fig. 6.2 shows the acceleration vs time plot of the ground
motion.
In the next sections the results of FLUSH analysis will be. presented for both
cases.
6.2.1 Interface with HINGE Support
Fig. 6.3 shows the distribution of horizontal normal stresses in the elements along
the boundary at the end of the static analysis. This figure contains the results of
both cases, hinged and roller, and also for the full section without the concrete
gravity dam at the same crosssection.
Fig. 6.4 shows maximum and minimum horizontal stress distribution at the
boundary elements. Finally, the results of static and dynamic analyses are super
imposed to get the final maximum and minimum stress distribution. Superimposed
stress distribution is shown in Fig. 6.5 The numbers next to data labels on the
minimum stress line tells that how many times that stress level can occur during
shaking, and since tensile stress means a possible separation that number tells how
many times the interface can separate during the earthquake. According to Fig.
6.5 it is possible to have separation along the interface down to approximately 140
ft. from the top.
6.2.2 Interface with ROLLER Support
Fig. 6.6 shows earthquake induced maximum and minimum horizontal stresses in
the roller boundary elements and finally Fig. 6.7 shows the superimposed hor
izontal normal stress distribution. According to Fig. 6.7 it is possible to have
separation along the interface almost down to 120 ft. from the top for the roller
support.
67
\ (psl)
* Upstream soilconcrete section
Perfectly bonded Interface
b Upstream soilconcrete section
Roller interface
Figure 6.3: The distribution of horizontal normal stress along the boundary, at the
end of the static analysis. (Chang et. al, 1992)
6S
r, {psO
Elevitlon (feet)
Figure 6.4: Maximum and minumum horizontal normal stresses at the hinge
boundary elements at ST. 299f35. Chang et. al, 1992. (Tension negative)
69
(pit)
Elevation (feet)
Figure G.5: Superimposed horizontal normal stress in the hinge boundary elements.
Chang et. al, 1992. (Tension negative)
70
ff, (pU
Elevation (feet)
Figure 6.6: Maximum and minumum horizontal normal stresses in the roller bound
ary elements. Chang et. al, 1992. (Tension negative)
71
(piO
Elevation (feet)
Figure 6.7: Superimposed horizontal normal sresses in the roller boundary ele
ments. Chang et. al, 1992. (Tension negative)
72
6.3 Conclusion
In this FLUSH analysis, possible separation is investigated by tracing the normal
horizontal tensile stresses in the elements at the roller and hinge boundary. From
the outputs it is easy to see that in the upper part of the interface tensile stresses
can occur more than once which is likely to cause separation. Hinge boundary
case shows a 140 ft. separation depth and roller boundary case shows around 120
ft.. Maybe the results of roller boundary case can be accepted more reasonable,
because roller support can simulate the interface contact condition better by al
lowing movement in the vertical direction. Of course it is impossible to say that
these results are completely right, at least they give some idea to understand the
interface behavior.
It was necessary to analyze this separation problem with more powerful analysis
tools such as nonlinear dynamic analysis code, elastoplastic material models, and
a good interface modeling. Thats why same problem is analyzed with NIKE2D
which has better analysis capability than FLUSH.
73
7 Dynamic Analysis with NIKE2D
7.1 Introduction
As already recommended in the previous analysis with FLUSH, the same separa
tion problem is analyzed with a better computer code NIKE2D. Capabilities of
NIKE2D are explained in Chapter 5. The main reason of choosing that computer
code is the already implemented slideline algorithm that allows separation and
frictional sliding and nonlinear material models allowing to determine permanent
deformations. In the next sections the procedure followed in the analysis will be
explained in detail.
7.2 Input to NIKE2D
The steps of preparing input parameters can be listed as below.
preparing the finite element mesh,
determining the material properties,
specifying interface contact conditions.
All of the items above can be done by using the preprocessor MAZE which
prepares the input file for NIKE2D. In the following sections those items will be
explained.
74
Figure 7.1: Finite Element Mesh of Left Wing Dam at ST. 299+35.
7.2.1 Finite Element Mesh
Finite element mesh is prepared by using MAZE. In order to get meshes first several
lines are defined. These lines form the outline of different parts. Then parts are
defined with surrounding lines by dividing into meshes in the desired manner.
Fig. 7.1 shows the Finite Element of Left Wing Dam generated by MAZE. As a
boundary condition the bottom of the dam is assumed as a fixed connection which
doesnt allow movement in any direction. In real life it is not the case, because base
can slide in the horizontal direction too. Sliding base condition can introduce more,
complexity and difficulties in the analysis. Thats why this case is not included in
the present analysis and left as a further study.
Since the Left Wing Dam is located on a rock formation, foundation part is not
included in the mesh and assumed as a rigid foundation by assuming fixed support
condition at the base.
Finer mesh is used for concrete section in order to prevent the penetration of
soil into the concrete. Penetration problem is discussed in Chapter 5.
75
7.2.2 Determination of Material Properties
In the present analysis two different material models are used: elasticity and
RambergOsgood elastoplasticity. Elastic model is used for concrete and Ramberg
Osgood model is used for soil. Fig. 7.2 shows all the material properties of Left
Wing Dam at ST. 299135 including RambergOsgood parameters.
Elastic properties of concrete are obtained from reference [21] For soil part
Gmax values are obtained based on the shear wave velocity distribution showed in
Fig. 4.2. These are the same values used in the previous FLUSH analysis.
RambergOsgood parameters are obtained by the computational procedure
which is proposed by Ueng and Chen (1992). This procedure calculates the
RambergOsgood parameters based on the test results or typical relations in terms
of modulus and damping ratios versus strain. So first step to obtain material pa
rameters is to choose G/Gmax and damping ratio versus shear strain curves. For
G/Gmax versus shear strain, Seeds (1984) typical curve for gravelly soil, is used.
This curve shows a good agreement with the test result of Folsom Gravel as shown
in Fig. 4.1 As damping ratio versus shear strain the curve from Seed et al.
(1984) is used. Fig. 7.3 shows Seeds typical curve of damping ratio versus shear
strain used in the analysis. Ueng and Chens (1992) procedure to obtain Ramberg
Osgood material parameters is explained later in this section.
Determination of RambergOsgood Material Parameters: The backbone (mono
tonic loading) strainstress relation of the RambergOsgood elastoplastic model
can be expressed by:
76
Figure 7.2: Material properties of Left Wing Dam at ST.
MATERIAL PROPERTIES USED IN THE DYNAMIC ANALYSIS
RamberqOsqood Material Parameters
Description Mail. No Unit Weight (pcf) Shear Wave velocity (ft/s) Max. Shear modulus (Gnu) (KsO reference shear strain r,( 10) reference shear stress V (KsO R alia (a)
Stiell: Coarse Dredgetailings, unsaturated 137 900 3145 0.696 0.219 2.218 0.9013
2 137 1100 4700 0.696 0.327 2.218 0.9013
3 137 1250 6065 0.696 0.422 2.218 0.9013
Core: Decomposed granite, unsaturated 4 138 1300 7190 1.87 1.349 2.246 0.7961
5 138 1400 7450 1.87 1.398 2.246 0.7961
6 138 1700 12745 1.87 2.39 2.246 0.7961
Core: Decomposed granite, saturated 7 142 1700 12745 1.09 1.39 2.164 0.8606
8 142 1300 7450 1.09 0.811 2.164 0.8606
9 140 1000 4190 1.08 0.456 2.164 0.8606
Shell: Coaise dredgetailings, saturated 10 140 900 3145 0.491 0.154 2.276 0.9423
11 140 800 2590 0.491 0.127 2.276 0.9423
12 140 1000 4040 0.491 0.198 2.276 0.9423
13 140 1150 5340 0.491 0.262 2.276 0.9423
Concrete 14 E 8.49*104 (Ksfl v = 0.2
to
CO
CO
+
CO
. to
s
tf 12
, A WflumoH ond Mori (19611 llordln (19651 a o s 
Q MoliuiNlo, Klshldo ond KroC967l Sflvf and S*d (1969) A Oonovon 1(9691
Hardin ond v Klthldo end Drnevleh (19701 Tckono (19701 T / / / / / / / r /. /
/ / / 0 /
?yt / r A O
* 9* *'
SlltU STRAIN 7CRCCHI
Figure 7.3: Variation of Damping Ratio with Shear Strain for Sands (Seed at al,
1984)
78
Figure 7.4: Typical loading and unloading curve for RambergOsgood Model.
lv T<
r1'
(7.1)
where 7=shear strain, r=shear stress, 7y=reference shear strain, ry=reference
shear stress, a=constant > 0, and r=constant >1.
Fig. 7.4 shows a. typical loading and unloading curve for RambergOsgood
model. For unloading and reloading, according to Masings rule the relation be
comes:
7 7o T T0
27y
2 tv,
1 + a
T ~ Ty
2r
r~lN
(7.2)
where 70=shear strain at point of stress reversal, and r0=shear stress at point
of stress reversal.
79
The values of 7,,, Ty, a, and r are to be determined according to the material
properties. By rearranging Eqn. 7.1 the secant modulus for the backbone curve
can be expressed as:
(
G0 = =
7 7y
1
\
1 f a
r1
/
For very small strain, i.e., 7 .>0 and t 0, since r >1,
(73)
(G 0)7=0 = Gmax = ~
ly
Then the backbone relation can be rewritten as:
(7.4)
7
rr
y Q t r i M
ly ^ma x'y
Therefore, besides Gmax 5 other parameters 7y, q, r should be determined for
the RambergOsgood model.
Substituting r = Go'y and rearranging Egn. 7.5 we get:
Gt>
Go
 1
a
GoJ
Gmax7y
r1
!) = lga + (r l)log
G7
(7.6)
(7.7)
Gmax7y /
Eqn. 7.7 can be plotted using only G/GTOar data and then the values of r and
a can be determined from slope and intercept respectively. Next step is plotting a
similar curve using only damping ratio curve and obtaining the same parameters
one more time. Final values of r and a are determined by taking the average, a
and t are determined using only damping ratio curve as follows:
The equivalent critical damping .ratio, /3, for a hysteresis loop with the tip at
(7o? To) can be expressed as:
SO
AE 2a(rl). G0 >r,7oyi
27rr070 7r (r + 1) Gmax 7,j
(7.8)
where AE= energy dissipation in one loading cycle. Substituting Eqn. 7.6 in
Eqn. 7.8 we get;
or
/? =
2(r 1)
7r(r + 1)
(1
Go_
Jmax
)
G0 _ (r + 1)
Gmax 2(r 1)
Substitute Eqn. 7.10 in Eqn. 7.7 we obtain;
(7.9)
(7.10)
log
/?7r(r + 1)
2 (r 1) /?7r(r + 1)
= log(a) + (r l)log
1 
j3ir(r + 1)
2(r 1) J 7y
(7H)
Finally Eqn. 7.11 can be plotted like Eqn. 7.7. Thus, a best fit straight line,
and values of a and r can be found for data including both modulus and damping
data.
If we examine Eqn.s 7.7 and 7.11 we see that both have 7y. Therefore to
determine yy and the other parameters an iterative procedure is followed as below;
1. Assume a value for 7y and obtain the values of a and r by plotting the data
according to Eqn.s 7.7 and 7.11.
2. Compute 7y according to Eqn. 7.8 from the given modulus and obtain an
average value of 7y.
3. Compare the new value of 7y with the previous value. Repeat steps 1 and 2
if the difference is too great.
81
Finally a, 7, and r are all obtained. Then ry can be calculated by using Eqn.
7.4 and this finishes the procedure for obtaining the RambergOsgood material
properties.
7.2.3 Interface Contact Conditions
NIKE2D uses penalty formulation for the slidelines with separation and frictional
sliding. Penalty formulation is already explained in Chapter 5. The most important
point in penalty formulation is to control the penetration. Basically penetration is
dependent on the following criterions;
penalty number,
mesh size of penetrated material, and
bulk modulus of the penetrated material.
High penalty numbers decreases the penetration but also introduces huge num
bers in the calculation which can cause numerical errors. NIKE2D uses a default
penalty number unless another number is specified by the user. The default num
ber is chosen by authors of NIKE2D so that it can fit many problems without
causing numerical trouble. Depending on the problem penalty number may have
to be changed. According to the experience in preparing this thesis it is strongly
recommended to play with the mesh size of the penetrated material first. If there
is still significant amount of penetration then bulk modulus of the penetrated ma
terial can be increased. Changing penalty number should be the last thing to
try. .
In this study increasing the number of meshes in concrete eliminated the pene
tration of soil into the concrete. Fig. 7.1 shows relative mesh size of concrete with
respect to other parts.
82
Figure 7.5: Original and scaled 1967 Koyna Dam Earthquake motion.
7.2.4 Ground Motion
In the analysis with NIKE2D, Koyna Dam Earthquake Record (see Fig. 6.2 )
is used. This record is a rock base motion and has a maximum acceleration of
.87g and duration around 10 sec. By scaling Original Koyna Dam Earthquake
record by .75 and .5 three different rock base motions are obtained with maximum
accelerations .87g, .65g and .44g. These three input ground motions are shown
in Fig. 7.5 The purpose of scaling the original record is to observe the effect
of ground motions having different maximum accelerations on the separation of
the interface. Results of the three different analyses will be presented in the next
section.
S3
Figure 7.6: Horizontal displacement time histories of nodes on the soil side at the
interface, (with original Koyna Dam Earthquake Record)
7.3 Results of Dynamic Analysis
Each analysis is performed in two parts: static an dynamic. Simulation of gravi
tational force is done by applying gravitational acceleration in time domain; from
O.Og at time 0.0 sec. to l.Og at time 100.0 sec. in 1 sec. intervals. After time
100.0 sec. dynamic analysis is initiated with horizontal input acceleration data.
Gravitational acceleration is applied gradually because subjecting the structure to
gravity in a very short time ca.uses a kind of shock effect on the dynamic analysis.
84
Figure 7.7: Horizontal displacement time histories of nodes on the soil side at the
interface, (with scaled Koyna Dam Earthquake Record by 0.75)
S5
node 444 d = 0.0 ft
node 441 d s 54 ft
g
Q.
M
a
0.2
0.0
0.2
0.4
0.6
0.8
1.0
. i.i.
i i 1 1 1 1 1 1 1 1 I I ( j t l i 1 1 1 1 1 1 1 1 1 1 1 t 1 1 i : i t 1 1 1 1 1 1 1 1 1 1 1
fHHKiHhH'ji
Hi IHI HlH'4H i i 1 j H H il H 14
e
1A
a
0.2
0.0
0.2
0.4
0.6
0.8
1.0
i.i..i.i. U44.U..
1 1 I 1 1 .<.4.4.4>l>4.4>l. 1 1 t 1 I 1 1 1 1 i*!Vf *** ittii
. .t.4.444 444444 TiTHTf'!!* 4444444>4> ii i i ii i4 i 4i **!*!T!T! 4.4* Hiii4 rh*?' i i i i
oooaooooooooooao
d^Nri^irioNooidrnri'f'm
OOOOOOOOOOOOOOO
OTrinviri(DrcooidTnr)T
time (sec)
node 443 d = 18 ft
node 440
time (sec)
d = 72 ft
0.2 ,
0.2 t i i i\f < i i i i i i i i i i i i i i i i i i i I
0.4 44444l+44i44l4i44i444+l44i44
0.6 .4.1.1.l.l..1.2.4.1.i.i..i.t.l.i.j..l.j. 1.1. 4.1.1.11. 1.1.1.
0.8 
1.0 i i i i i vi f i > i i i i i i i* i ii i i i i i ii *
ooooooooooooooo
drrin>r)(i)Nioidrnn7
0.2 ri ititi
^ 0.0 Â£ 0.2 5. 0.4 2 0.6 Q 0.8 1.0 i i i i i i I i i < i i i i i t i i i i i
 i i i rrl'i1' .i.l.j.l. <.44.1 J j I j ;! T  ! : YT*:*Y f* .i.il.i.j.1. 1 1 I i i i ! > i *T*i,i*f* .1.1.i.l. .J. .1.2.1. i i i < i i i.j..Lj.i.i. ! ! ! rhrj till* S..I.4..S.I
OOOOOOOOOOOOOOO
oVtNrtvmiDNOoidtNrjn
time (sec)
time (sec)
scaled by 0.5 node 442
node 439 d90ft
e
a
*
a
0.2
0.0
0.2
0.4
0.6
0.8
1.0
ItrHvr .i.i.i.i.i.i.
rrrrvrt I'yrrn frnTfVITTTV'f
iiittit t > t i i i* i i i i i i i*i i i 7 i i i i i i i i i i > i i i j j* i rrrv ttv;
.
w
0.2
0.0
0.2
0.4
0.6
0.8
1.0
i i iiiitiitiii ittii
rrr U.i. T7T ; t t ? ? ; t j ; t i r i i i i r i i j Tt'i i'i"t'7'iT7i T .i 'iT'ji'f u~m. rrrr   j i.UU "V7*iii
ooooooooooooooo
d r ri ri v a d s d n o' r n ri v
ooooooooooooooo
drnnTiiiidNideidrHni
time (sec)
time (sec)
Figure 7.S: Horizontal displacement time histories of nodes on the soil side at the
interface, (with scaled Koyna Dam Earthquake Record by 0.5)
86
