See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/283503214

Application of SANISAND Dafalias-Manzari model in FLAC3D

Conference Paper · October 2013

CITATIONS READS

6 1,947

3 authors, including:

Zhao Cheng
Itasca Consulting Group, Inc
25 PUBLICATIONS   233 CITATIONS   

SEE PROFILE

All content following this page was uploaded by Zhao Cheng on 05 November 2015.

The user has requested enhancement of the downloaded file.

 Continuum and Distinct Element Numerical Modeling in Geomechanics -- 2013 -- Zhu, Detournay, Hart & Nelson (eds.) Paper: 09-03
©2013 Itasca International Inc., Minneapolis, ISBN 978-0-9767577-3-3

Application of SANISAND Dafalias-Manzari model in FLAC3D

Z. Cheng
Itasca Consulting Group, Inc., Minneapolis, MN USA
Y.F. Dafalias
University of California, Davis, CA, USA & National Technical University of Athens, Hellas, Athens
M.T. Manzari
The George Washington University, Washington, DC USA

ABSTRACT: The Dafalias-Manzari (2004) version of SANISAND model, an elastoplastic con-

stitutive model for granular liquefiable soils, is implemented as a user-defined model in the plat-
form of FLAC3D, a three-dimensional explicit finite difference program using fast Lagrangian
analysis of continua for engineering mechanics computation. This is a trial of fully nonlinear
analysis method, which contrasts with the commonly used equivalent-linear method, in ge-
otechnical earthquake engineering. The simulative capability is validated by comparing the re-
sults of simulation with drained and undrained triaxial compression laboratory data over various
mean effective pressures and void ratios. Simulations of cyclic stress ratio versus number of cy-
cles, back-calculation of strain-compatible shear moduli and damping ratios, free-field nonlinear
response, and slope dynamic analysis in terms of effective stresses are presented using the mod-
el in FLAC3D.

1 INTRODUCTION

Two essential soil properties, when soils are subjected to dynamic loading, are strain-dependent
shear modulus and damping ratio (Ishihara 1996). The equivalent-linear method uses an itera-
tion procedure to account for the changes of shear modulus and damping ratio with the shear
strain level, while during the solution process within each iteration, the soil is still linear-elastic
(Kramer 2005). The equivalent-linear method has been in use for many years to simulate wave
propagation in geotechnical media subjected to earthquake excitation. The equivalent-linear
method has the merit that it is user-friendly and accepts experimental cyclic test results directly.
It should be noted that the equivalent-linear model is only an approximation of the actual fully
nonlinear stress-strain relations of the soil. The method does not directly capture any nonlinear
constitutive effects and strain-dependent modulus and damping ratio are only taken into account
in an average sense. With equivalent-linear methods, the soil will be over-damped during quiet
periods in the shaking history and under-damped during strong excitation. The equivalent-linear
method cannot directly predict irreversible or permanent deformations and failures. Equivalent-
linear methods are based on total stresses and cannot track the build-up, redistribution, and dis-
sipation of excess pore pressure during and after a seismic loading. Fully nonlinear methods can
be formulated in terms of effective stresses and the fluid-mechanics interaction behaviors can be
simulated for stages of before, during and following dynamic shaking. Elastoplastic constitutive
models are shown to capture the essential features of soil's stress-strain-strength behavior. It is
noted however that despite their versatility and effectiveness in modeling soil behavior, the
adoption of these models for use in engineering practice has been rather slow (Anandarajah
2005). A key reason for this slow adoption in practice is the apparent need for many model pa-
rameters and a lengthy calibration process. Moreover, efficient calibration of these models usu-
ally requires strong background in elastoplastic theory, which is not an integral part of graduate
programs in geotechnical engineering. Hence to make these models more accessible to practic-
ing engineers, it is very important that they are implemented in the software packages that are
commonly used by geotechnical engineers. Moreover, these implementations must be supple-
mented by solution of practical problems that geotechnical engineers face in practice. The goal
of this paper is to provide an example of such implementations for an elastoplastic constitutive
model for granular liquefiable soils.
The plasticity sand model developed by Dafalias & Manzari (2004) is a stress-ratio con-
trolled, critical state compatible sand model, which takes into account the effect of fabric chang-
es in the multi-axial generalization. It is a version of a family of SANISAND (Simple ANIso-
tropic SAND) models developed by Dafalias and collaborators (Manzari & Dafalias 1997, Li &
Dafalias 2000, Dafalias & Manzari 2004, Dafalias et al 2004, Taiebat & Dafalias 2007). One of
the merits of this model is that one set of material parameters can be applied to dramatically dif-
ferent stresses and densities. The material parameters have been well calibrated for Toyoura
sand. Some material parameters are considered default values for most sands and need not to be
determined separately for other sands. The Dafalias-Manzari (2004) model involves a number of
highly nonlinear rate equations that are difficult to integrate using fully-implicit integration
techniques as the integration process involves second derivatives of yield and plastic potential
functions and leads to rather long expressions involving various operations on tensorial quanti-
ties. On the other hand, fully implicit techniques are shown to be the most robust, efficient, and
accurate methods when used in implicit nonlinear finite element analyses. Manzari and
Prachathananukit (2001) have presented a comprehensive analysis of explicit, semi-implicit, and
fully implicit integration techniques for an earlier version of the current model within the con-
text of finite element analysis. It is shown that for quasi-static problems, implicit methods are
clearly superior to the other techniques and can provide highly accurate and robust solution par-
ticularly when highly strained zones are localized within the soil body. However for seismic
analysis when a large number of small steps are needed, semi-implicit techniques (e.g., sub-
stepping technique) can provide very accurate and efficient solutions (Sloan 1987, Sloan et al
2001, Zhao et al 2005). Jeremić et al (2008) and Cheng & Jeremić (2009) have also implement-
ed and applied a semi-implicit integration technique within the finite element framework.
FLAC3D uses an “explicit” time integration scheme which can follow arbitrary nonlinearity in
stress-strain laws in almost the same computer time as linear laws. Explicit time integration
methods also have the advantage that they do not require assembly of global matrices. The
“mixed discretization” scheme (Marti & Cundall 1982) is used for accurate modeling of plastic
collapse loads and plastic flow. This scheme is believed to be physically more justifiable than
the “reduced integration” scheme commonly used with finite elements (Itasca 2012a). Users
may create their own constitutive model (user-defined model) for use in FLAC3D and take ad-
vantages of all plausible modules of FLAC3D, including fluid-mechanics coupling, multi-
threaded calculation, built-in free-field and quiet (viscous) boundary condition options for dy-
namic analysis (Itasca 2012b), and FISH, an embedded programming language. The above fea-
tures make FLAC3D a suitable platform to apply the Dafalias-Manzari model, particularly for
dynamics time history analysis.
In this paper, the Dafalias-Manzari model is implemented as a user-defined model for
FLAC3D. Validations of this model in FLAC3D are presented by comparison with a wide range of
triaxial laboratory data. Several examples including simulation of cyclic stress-strain behavior of
Toyoura sand, evolution of shear modulus and damping ratio with shear strain, free-field seis-
mic response of saturated sand deposit, and seismic response of a saturated soil slope are pre-
sented and discussed.

2 IMPLEMENTATION

The detailed description of this model is available in Dafalias & Manzari (2004), and the reader
can also refer to the earlier version of Manzari & Dafalias (1997) for more explanations of the
basis for its conception. This model is intended to simulate both monotonic and cyclic loadings.
In this paper, both stress and strain are assumed positive in tension, and pressure (e.g., mean
pressure, pore pressure) is assumed positive in compression, following the sign conversion in
FLAC3D. If not explicitly stated, the stresses are effective by default in this paper; the super-
scripts and denote the elastic and plastic parts of strain, respectively. The general stress
and strain tensors are expressed in the ( ) space by and , respectively. The mean

2
pressure is defined by ; the deviatoric stress component is defined by
, where for and for . The equivalent scalar-
valued deviatoric stress is expressed by .

2.1 Elastic Law

The elastic part of Dafalias-Manzari model is assumed to be hypoelastic. The expression of the
law in term of the rate (material time derivative, denoted by a superposed dot) has the form of:

(1.a)

(1.b)
Where is the elastic volumetric strain and is the elastic deviatoric shear strain, with the
definitions and . The hypoelastic bulk modulus and shear
modulus are not constants, but functions of the mean effective pressure and the current
void ratio :

(2.a)

(2.b)

where is a material constant, is the constant Poisson’s ratio and is a reference

pressure for normalization purpose (usually the value of the standard atmospheric pressure is
adopted).

2.2 Critical state line

The critical void ratio is related to the critical pressure , by the power relation (Li & Wang,
1998)

(3)

where , and ( for most sands) are material constants. The ultimate critical-
state stress ratio (q/p) is assumed a constant . The state parameter is defined as the differ-
ence between the current void ratio and the critical void ratio at the current pres-
sure: . The state parameter is an important measure of the material state.

2.3 Yield function

The yield function is expressed by

(4)
where is a material constant ( = 0.02-0.05) and is the deviatoric back-stress ratio ten-
sor characterizing the axis of the yield surface (Manzari & Dafalias 1997, Dafalias & Manzari
2004). Other two useful expressions are the deviatoric stress-ratio tensor
(5)
And the normalized tensor

(6)

where denotes the norm of a tensor. Note that and .

3
2.4 Evolution of plastic deviatoric strain
The rate of plastic deviatoric strain is following the law

(7)
where is the loading index (or plastic multiplier) to be defined later, the MacCauley brackets
is defined as if and if . is a scalar associated with
dilatancy, which will be defined in the following section 2.5. and account for the effect
of Lode angle on the direction of deviatoric plastic strain rate and are expressed by

(8)

(9)

where is a material constant denoting the ratio of the triaxial extensive strength to compres-
sive strength ( ) with θ the Lode angle; notice that when c=1, B=1 and C=0;
and the function are expressed by

(10)

(11)

2.5 Evolution of plastic volumetric strain

The rate of plastic volumetric strain is

(12)
Where the following quantities associated with dilatancy are defined by

(13)

(14)

(15)

and are material constants. The fabric-dilatancy tensor is defined in section 2.6.

2.6 Internal variables update

The rate of fabric-dilatancy tensor is

(16)
where and are material constants ( = 4.0-5.0 for most sands).
The rate of back-stress ratio tensor is

(17)
and

(18)

4
 (19)

where and are material constants. is the initial value of at initiation of a

new loading process and is updated not only when such process begins but also when the de-
nominator of Eq. (19) becomes negative.

2.7 Implementation procedure

The Dafalias-Manzari model is only applicable to sands or sand-like materials in which the
stresses correspond to a compressive mean effective stress. The model is not designed to predict
the behavior of materials in which this condition is not satisfied. In particular, before the appli-
cation of the model, the initial states (two important states are the initial mean effective pressure
and the initial void ratio) of the material must be specified. An error message will be issued if a
state is invalid (e.g., negative initial mean effective stress) or an invalid material parameter is
input (e.g., a Poisson’s ratio is greater than 0.5).
The elastic moduli are updated at the beginning of the time step according to the current mean
pressure and assumed constant during this time step. An elastic trial state is computed after add-
ing the old stress components to the stress increments obtained by assuming an elastic loading.
If the tried elastic stress violates the criterion for yielding, plastic deformation takes place and
the stresses, as well as the back-stress ratio tensor and fabric dilatancy tensor, must be updated
following the conventional elastoplastic theory. The stress tensor, the back-stress ratio tensor
and the fabric dilatancy tensor updates all depend on the loading index . The loading index
can be derived in terms of plastic modulus according to

(20)

(21)
Once the loading index is available, the rate of deviatoric and volumetric plastic strains and the
rates of internal variables (the back-stress ratio tensor and the fabric dilatancy tensor) can be
calculated according to the previous relations.
At the end of the time step, the rate of stresses are eventually calculated by the expression

(22)
The new values of the state variables and the new values of shear and bulk moduli are then
stored for use in the next time step. The material properties thus lag one time step behind the
corresponding calculation. In an explicit analysis, this error is small because the time steps are
small.
It is worth noting that the implementation has an option whereby the back-stress ratio tensor
can be initialized by inputting its six components as material parameters, or its initial values can
be automatically set to be components of the current deviatoric stress-ratio tensor based on the
initial stress state. The fabric dilatancy tensor can also be initialized by inputting its six compo-
nents as material parameters.

3 SIMULATIONS

Several simulation examples are presented here to demonstrate the simulative capability of
Dafalias-Manzari model in FLAC3D. In this paper, all examples use the same material parame-
ters for Toyoura sand which have been summarized and tabulated in Dafalias & Manzari (2004).

5
3.1 Validation to the triaxial test data
Triaxial laboratory data on Toyoura sand (Verdugo & Ishihara 1996), both drained and un-
drained, now serve as “standard” validation tests for the implementation of the Dafalias-Manzari
model in any computational platform. These validation tests consist of simulating material be-
havior in a wide range of initial mean pressures and void ratios.
The simulations on the triaxial compression undrained tests are for three initial void ratios of
0.735, 0.833 and 0.907 to stand for the dense, medium dense and loose sands, respectively. For
each void ratio case, there are several scenarios of the initial mean pressures (Figs. 1-3). From
Figure 2 for the medium dense sand, it is interesting to observe that the responses of contraction
together with softening for the higher values of mean pressures, and correspondingly dilatancy
together with hardening for the smaller mean pressure values.

p0 = 3000 kPa
p0 = 3000 kPa
p0 = 1000 kPa
p0 = 100 kPa

Figure 1. Simulative results (in blue) and experimental data (in red), undrained triaxial tests (e0 = 0.735).

p0 = 3000 kPa
p0 = 2000 kPa
p0 = 1000 kPa
p0 = 100 kPa

Figure 2. Simulative results (in blue) and experimental data (in red), undrained triaxial tests (e0 = 0.833).

p0 = 2000 kPa
p0 = 1000 kPa
p0 = 100 kPa

Figure 3. Simulative results (in blue) and experimental data (in red), undrained triaxial tests (e0 = 0.907).

6
 The simulations of the triaxial compression drained testes are for two initial mean pressures
of 500 and 100 kPa. For initial mean pressure of 500 kPa, there are three scenarios of the initial
void ratios: 0.960, 0.886 and 0.810; for initial mean pressure of 100 kPa, there are also three
scenarios of the initial void ratios: 0.996, 0.917 and 0.831 (Figs. 4 & 5).
It should be noted that in FLAC3D, the stress-strain relationship is computed within a tetrahe-
dron (subzone) and the zone (corresponding to element in finite element method) behavior is
averaged from its subzones by the weights of the subzone volume, e.g., a hexahedron zone has
two overlays and each overlay has five subzones (Itasca 2012a); while in finite element method,
the stress-strain relationship is computed at each integration point.

e0 = 0.810
e0 = 0.886
e0 = 0.960

Figure 4. Simulative results (in blue) and experimental data (in red), drained triaxial tests (p0 = 500 kPa).

e0 = 0.831
e0 = 0.917
e0 = 0.996

Figure 5. Simulative results (in blue) and experimental data (in red), drained triaxial tests (p0 = 100 kPa).

3.2 Simulation of cyclic shear

A one-zone example with an initialized mean effective pressure of 100 kPa and initial void ratio
of 0.833 is subjected to simple shear in the horizontal direction (X-direction) by velocity control
(the top surface is moving horizontally). The zone size is of unit length. This zone can be
coarsely considered as a sand sample under a depth of about 50 meters (assume the soil density
is 2×103 kg/m3) in an undrained condition (since the volume is kept constant) subjected to a cy-
clic horizontal shear.
Figure 6 shows the simulated shear stress versus vertical stress and shear stress versus shear
strain response for K0 = 0.5 & 1.0, where the coefficient K0 is defined as the ratio of horizontal
stress to vertical stress, . From these figures the vertical stress decreases to a value close
to zero after about three cycles and the shear modulus is reduced when the shear strain increases.
The degradation of the vertical effective stress, as well as the mean stress, implies that the soil
sample loses its vertical bearing capacity and any structure which is constructed above the posi-
tion may be damaged after this kind of cyclic disturbance. If the soil is saturated, the pore pres-
sure usually increases as the mean pressure is decreasing. This rapid decline of mean effective

7
stress due to the dynamic excess pore pressure build-up that reaches the initial vertical effective
stress leads to the loss of soil shear strength and is known as soil liquefaction.
It is interesting to plot the relationship of the cyclic shear ratio (CSR) and the required num-
ber of cycles to trigger liquefaction under different initial void ratios (e= 0.961, 0.833 & 0.735)
and under the different K0 values (K0 = 0.5 & 1.0) based on the simulations. The CSR is defined
as , where is the initial effective vertical stress, and is the horizontal shear
stress. The shear loading direction reverses when the absolute value of horizontal shear stress
reaches . In this example, when the vertical stress is less than 5% of its initial value (100 kPa
in this example) or the shear strain is greater than 1%, whichever comes first, the liquefaction is
considered to occur and the number of cycles is recorded at the current time step. The model
gives a reasonable representation of typical sand cyclic responses as shown in Figure 7. With the
same number of cycles, looser sand will require a lower CSR to trigger liquefaction. The coeffi-
cient K0 effect is small for loose sand, but bigger for denser sands.

Figure 6. Cyclic response of single zone subjected to simple shear.

Figure 7. Liquefaction resistance with an initial mean effective pressure of 100 kPa.

3.3 Back-calculation of strain-compatible shear moduli and damping ratios

The curves of modulus (usually represented by G/Gmax) and damping ratio (D) versus shear
strain are the key cyclic properties used for analyzing the dynamic response of soils and are
widely used in the equivalent-linear method. The important factors that affect the relations of
G/Gmax and damping ratio with shear strain include soil type, plasticity index, mean effective
pressure and geologic age. Other factors including frequency of loading, number of loading cy-

8
cles, over-consolidation ratio, void ratio, degree of saturation and grain characteristics appear to
have less-important effect on G/Gmax and damping ratio versus shear strain (Zhang et al., 2005).
Because the Dafalias-Manzari model is developed for clean sand with the plasticity index close
to zero, the only factor which can significantly affect strain-compatible G/Gmax and damping ra-
tios is the mean effective pressure if we neglect the geologic age.
It is noted that in the equivalent-linear method, the material constitutive model assumes that
the stress-strain hysteretic loop is in a shape of an ellipse. In nonlinear materials, a frequency-
independent hysteretic curve in the form of an ellipse is physically impossible. The damping ra-
tio in the nonlinear material is only approximately back-calculated from the loop area surround-
ed by the loading and unloading curve in one cycle and the triangle area representing elasticity
energy.
Numerical simulations of monotonic and cyclic triaxial tests at various constant mean pres-
sures (10, 20, 50, 100, 200 and 500 kPa) using Dafalias-Manzari model in FLAC3D are conduct-
ed to back-calculate the relations of G/Gmax and damping ratio with a shear-strain range of
0.0001% to 1%. From the results shown in Figures 8 & 9, it is observed that the model correctly
predicts that the cyclic threshold shear strain is greater at high mean pressures than at low mean
pressures and at the same shear strain, G/Gmax increases but the damping ratio decreases when
the mean pressure increases.

Figure 8. Back-calculated shear modulus versus shear strain.

Figure 9. Back-calculated damping ratios versus shear strain.

9
 The results are compared with the curves extended based on Seed & Idriss (1970), Idriss
(1990) & EPRI (1993), in which the curves at 0-20 feet represent a lower-bound (red dash lines
in Figs. 8 & 9) and the curves at 500-1000 feet for an upper-bound (blue dash lines in Figs. 8 &
9). The calculated G/Gmax and damping ratios versus shear strain at the different mean pressures
generally fit well between the upper- and lower-bound curves. From the back-calculated damp-
ing radio curves, the damping ratios are zero at shear strains below the cyclic threshold shear
strain because no yielding occurs, and thus no energy dissipation takes place. Experimental data
show that the damping ratio is never zero, even at very low strain levels (Kramer 1996). The
possible reason is that the back-calculated damping ratio is limited to the material internal ener-
gy dissipation by the constitutive law, and do not include other energy dissipations which occur
even at very small shear strains. In practice, this can be compensated by assigning small hysteric
damping or small stiffness-proportional Rayleigh damping which is automatically inhibited by
FLAC3D when material yielding occurs. The results in Figures 8 & 9 show that, the Dafalias-
Manzari model can characterize the essential features of sand cyclic behavior in terms of the two
cyclic material parameters of the equivalent-linear method: strain-compatible G/Gmax and damp-
ing ratios.

3.4 Free-field response

In this section, we simulate the problem of wave prorogation within a geological profile using
FLAC3D with Dafalias-Manzari constitutive model. Results are compared to those obtained from
SHAKE91 (Idriss & Sun 1992), which uses an equivalent linear method. The FLAC3D and
SHAKE91 models simulate the following problem: a horizontally layered saturated clean sand
deposit overlying rigid bedrock is subjected to a horizontal acceleration bedrock motion. The
model is 150 feet (45.72 m) deep. In FLAC3D, the evolution of the soil shear moduli and damp-
ing ratios are entirely specified by the constitutive law. The dynamic characteristics of these
soils in SHAKE91 are governed by the EPRI (1993) depth-dependent curves of strain-
compatible moduli (by G/Gmax) and damping ratios for sands. We increase the stiffness (shear
modulus) of the bedrock to a large value in the SHAKE91 data file in order to simulate rigid
bedrock for the purpose of comparison. The input bedrock acceleration motion is the ground
motion recorded in the Loma Prieta Earthquake (1989). The input acceleration history is shown
in Figure 10. The record has a peak acceleration of approximately 0.11 g (g is the acceleration of
gravity), duration of 40 seconds, and dominant frequency of 2.47 Hz.

Figure 10. Input Acceleration at the bedrock.

In FLAC3D, the model consists of a 1D column composed of 30 cubic zones, each with a
depth of 5 feet (1.524 m). Although uniform soil properties are assigned, the vertical stresses are
a function of depth, so the initial stresses are different in each zone. Model boundary conditions
are configured to simulate the free-field motion of a layered soil deposit with a rigid base. The
horizontal acceleration record is applied as an input boundary condition at the base of the model.

10
 In SHAKE91, the soil shear moduli are equal to Gmax used in FLAC3D. An iterative solution is
performed until the secant shear moduli and the damping ratios throughout the column converge
to a given tolerance corresponding to 50% of the peak shear strain. Calculations are made for
different input accelerations scaled to the peak base acceleration with different scaling factors.
Figure 11 shows the simulated horizontal stress versus vertical stress at Zone 25 (at depth of
5.03 m or 16.5 ft) for the peak acceleration scaled to 0.1 g and 1.0 g in the FLAC3D model. For
the peak acceleration scaled to 0.1 g, the zone vertical stress reduction is small; the zone loses
all vertical stress for the peak acceleration scaled to 1.0 g, for which liquefaction is likely to oc-
cur.
For validation purposes, acceleration amplifications near the ground surface, defined as the
ratio of the peak acceleration value at the top to the peak acceleration value at the base, are
compared between NL method in FLAC3D and EL method in SHAKE91. We measure the ratios
when applying scaled maximum input acceleration value of 0.0001 g, 0.0005 g, 0.001 g, 0.005
g, 0.01g, 0.05 g, 0.1 g, 0.5 g and 1.0 g. FLAC3D simply takes the scaled base motion from
SHAKE91 as the input motion. Two cases are considered here: one for K0 =1.0 and the other for
K0 = 0.5. SHAKE91 does not directly have the K0 effect. The K0 effect is indirectly included in
SHAKE91 by inputting corresponding shear moduli under the condition of K 0, because Gmax of
the sand is a function of the initial mean pressure and therefore a function of K 0 coefficient as
well. According to the calculation results shown in Figure 12, it is confirmed that FLAC3D and
SHAKE91 predict generally similar acceleration amplification.

Figure 11. Horizontal shear stress (Sxz) vs. vertical stress (Szz) at Zone 25 (at depth of 5.03 m or 16.5 ft).

Figure 12. Acceleration amplification comparison at the top of the model for K0 = 1.0 (left) and K0 = 0.5
(right) in SHAKE91 and FLAC3D.

11
3.5 Seismic response of a slope
This example is a soil slope subjected to seismic response. The slope angle is as small as 5.71
(V: H = 0.1) degree and the height is 36 m at one end and 30 m at the other end. The slope is as-
sumed lying on relatively hard bedrock. The model is divided into three layers in which a sand
layer is sandwiched by two clay layers. The slope is assumed to be a plane-strain problem for
simplicity although the grid is three dimensional. The finite differential grid is shown in Figure
13. There are four analysis stages for this simulation.
The first stage is pure fluid flow calculation to initialize the pore pressures (fluid on, mechan-
ics off) for the dynamic analysis. Hydrostatic pore pressures are assumed fixed at two slope ends
(X = 0 and 120 m). At the slope surface, the pore pressure is set to zero. The tensile limit for wa-
ter is set to zero so that a phreatic surface develops. After the model achieves equilibrium, the
pore pressures within the whole model with the corresponding phreatic surface are obtained.
The second stage is to set the initial stresses (fluid off, mechanics on). Roller conditions are
set up for the sides and the base. The model reaches equilibrium with zero shear stresses on the
boundary planes. All material groups are assumed Mohr-Coulomb materials with large cohe-
sions and large tension limits at this stage for stress initialization purpose. The water bulk modu-
lus is set to zero to avoid any pore pressure alteration due to the change of the soil matrix vol-
ume.

Figure 13. Finite differential grid in the XZ plane.

The next stage in generating the initial static condition consists of setting realistic strength
properties for the clay layers and assigning the Dafalias-Manzari material properties to the sand
layer (initial void ratio is set to 0.86). The stresses obtained from the previous stages are auto-
matically taken as the initial stresses for this stage. After the model reaches the equilibrium
again, all zone groups have the stresses and pore pressures with realistic material properties. The
deformations and strains are then set to zero again but stresses remain in equilibrium with the
gravity so that the dynamic analysis in the next stage begins from a non-deformed grid but with
realistic stresses, pore pressures and void ratios.
The final stage is the dynamic simulation with full mechanics-fluid interaction (fluid on, me-
chanics on). The boundaries at the lateral ends are changed to free-field boundaries (More de-
tails on free-field boundary can be found in Itasca 2012b). In this way, plane waves propagating
upward suffer no distortion at the lateral boundaries because the free-field grid supplies condi-
tions that are identical to those in a laterally infinite model. At the base, the boundary is as-
sumed rigid and an acceleration time history representing a SH-wave is input. The input acceler-
ation motion is shown as GP 2039 in Figure 14, which is scaled to make the peak acceleration to
be 0.1 g (red line of the left graph in Fig. 14) and is matched to make the peak response spec-
trum to be about 0.234 g (red line of the right graph in Fig. 14). The duration of the input accel-

12
eration motion is 44 seconds. At this stage, a realistic value of the water bulk modulus and real-
istic permeability coefficients are required. The water bulk modulus is set to be 2 GPa and the
permeability coefficients are 1×10-5 m/s for the sand layer and 1×10-7 m/s for the clay layers.
The plastic flow associated with the Mohr-Coulomb model in the clay layers dissipates consid-
erable energy at high excitation levels. The selection of additional damping is less critical when
yielding occurs. Here the default hysteretic damping with parameters calibrated for clays (Itasca
2012b) is used to simulate the cyclic behaviors of the clay. In the sand layer modeled by the
Dafalias-Manzari model, no hysteretic damping is added. In order to reduce the high frequency
noise, a small amount of stiffness-proportional Rayleigh damping is applied to the whole model.
Both hysteretic and Rayleigh damping are switched off automatically by FLAC3D when plastic
flow is occurring at any zone.
In Figure 14, the blue line is the calculated acceleration motion at GP 5581, which is the mid-
dle point on the slope. The peak acceleration at GP 5581 is approximately 0.2 g which is ap-
proximately double the input peak acceleration at the base. There are still some high frequency
noises after the excitation duration (after 44 seconds. The calculated response acceleration spec-
trum at GP 5581 is generally larger than that at the base at longer periods (above 0.5 seconds in
this example), which is consistent to the local site condition effect that soft soil deposits yield
greater proportions of high frequency motion (Seed et al 1976).
Figure 15 shows the displacement time histories at selected grid points. Irrecoverable perma-
nent displacements at the grid points on the slope are predicted due to the substantial plastic
flow in the fully nonlinear analysis. The left graph in Figure 15 plots X-displacement time histo-
ries at three selected grid points; where GP 2039 is at the base of the model, representing the in-
put displacement motion, GP 1921 and 5641 are at the top and the toe of the slope, respectively.
It is observed that even the input motion at the model base has been baseline-corrected, the out-
put displacement motions on the slope can still be drifted from its original position, which can-
not be functionally predicted by the equivalent-linear method. The right graph in Figure 15 plots
the Z-displacement time history at the middle grid point on the slope (GP 5581). A permanent
settlement of about 0.05 m at this grid point is predicted by the model.
Figure 16 presents the pore pressure time histories at two selected zones in the sand layer.
Figure 17 shows the excess pore pressure snapshot at the end of the excitation (at 44 seconds).
The excess pore pressure build-up in the sand layer which modeled by the Dafalias-Manzari
model is as expected because the model includes the effect of cyclic shear strain on volumetric
response of the soil matrix when the shear wave is propagating through the layer. At the quiet
periods of the excitation and even after the excitation is ceased (after 44 seconds), the pore pres-
sures are still dissipating slowly, which are aftereffects observed following numerous earth-
quakes. The generation and dissipation of the excess pore pressure cannot be simulated by the
equivalent-linear method because it is in terms of total stresses, not effective stresses.

Figure 14. Acceleration time histories (left) and acceleration response spectra (right) at GP 2039 (on the
base) and at GP 5581 (on the slope).

13
Figure 15. Displacement time histories (left: in X-direction; right: in Z-direction) at selected grid points.

Figure 16. Pore pressure time histories at selected zones.

Figure 17. The excess pore pressure distribution at the end of the excitation (44 seconds).

14
4 SUMMARY

The SANISAND Dafalias-Manzari model was implemented and applied in FLAC3D as a user-
defined model. The model validations were presented by comparing the simulation results with
the experimental data of triaxial compression tests from a wide range of initial mean pressures
and void ratios in both undrained and drained conditions. Cyclic shear loading of sand, compari-
son of the derived strain-compatible shear moduli and damping ratios with those used by the
equivalent-linear method, an example of response of saturated sand deposit, and dynamic re-
sponse analysis of a slope with fluid-mechanics interactions subjected to seismic excitation are
simulated using the Dafalias-Manzari model in FLAC3D. The model is appropriate to be used for
more complex static and dynamic problems to take advantages of all features of fully nonlinear
analysis.

ACKNOWLEDGEMENTS

Z. Cheng acknowledges constructive discussions with Dr. C. Detournay and Dr. Varun. Y.F.
Dafalias and M.T. Manzari acknowledge funding from the European Research Council under
the European Union's Seventh Framework Program (FP7/2007-2013) / ERC IDEAS Grant
Agreement No. 290963 (SOMEF). Y.F. Dafalias acknowledges partial support by NSF project
CMMI-1162096. M.T. Manzari acknowledges the support provided by the National Science
Foundation Geotechnical Engineering Program.

REFERENCES

Anandarajah, R. (organizer) 2005, Workshop on Nonlinear Modeling of Geotechnical Problems: From

Theory to Practice, Johns Hopkins University, Maryland, November 3&4, 2005.
Cheng, Z. & Jeremić, B. 2009. Numerical modeling and simulation of pile in liquefiable soil. Soil Dynam-
ics and Earthquake Engineering 29(11-12): 1405-1416.
Dafalias, Y.F. & Manzari, M.T. 2004. Simple plasticity sand model accounting for fabric change effects.
Journal of Engineering Mechanics 130(6): 622-634.
Dafalias, Y.F., Papadimitriou, A.G. & Li, X.S. 2004. Sand plasticity model accounting for inherent fabric
anisotropy. Journal of Engineering Mechanics 130(1): 1319-1333.
EPRI 1993. Guidelines for determining design basis ground motions. Final Report, No. TR-102293, Elec-
tric Power Research Institute, Palo Alto, California.
Idriss, I.M. 1990. Response of soft soil sites during earthquakes, Proceedings of Memorial Symposium to
Honor Professor H. Bolton Seed, Berkeley, California, Vol. II, May.
Idriss, I.M. & Sun, J.I. 1992. SHAKE91: A computer program for conducting equivalent linear seismic
response analyses of horizontally layered soil deposits. Center for Geotechnical Modelling, Depart-
ment of Civil and Environmental Engineering, University of California, Davis.
Ishihara, K. 1996. Soil Behaviour in Earthquake Geotechnics. Clarendon Press and Oxford University
Press.
Itasca Consulting Group, Inc. 2012a. FLAC3D – Fast Lagrangian Analysis of Continua in Three-
Dimensions, Ver. 5.0, User’s Guide Manual. Minneapolis: Itasca.
Itasca Consulting Group, Inc. 2012b. FLAC3D – Fast Lagrangian Analysis of Continua in Three-
Dimensions, Ver. 5.0, Dynamic Analysis Manual. Minneapolis: Itasca.
Jeremić, B., Cheng, Z., Taiebat, M. & Dafalias, Y. 2008. Numerical simulation of fully saturated porous
materials. International Journal for Numerical and Analytical Methods in Geomechanics 32(13): 1635-
1660.
Kramer, S.L. 1996. Geotechnical Earthquake Engineering. Upper Saddle River: Prentice Hall.
Li, X.S. & Dafalias, Y.F. 2000. Dilatancy for cohesionless soils. Géotechnique 54(4): 499-460.
Li, X.S. & Wang, Y. 1998. Linear representation of steady-state line for sand. Journal of Geotechnical
and Geoenvironmental Engineering 124(12): 1215-1217.
Manzari, M.T. & Dafalias, Y.F. 1997. A critical state two-surface plasticity model for sands.
Géotechnique 42(2): 255-272.
Manzari, M.T. & Prachathananukit, R. 2001. On integration of a cyclic plastic model. International Jour-
nal for Numerical and Analytical Methods in Geomechanics 25(6): 525-549.

15
 Marti, J. & Cundall, P. 1982. Mixed discretization procedure for accurate modelling of plastic collapse.
International Journal for Numerical and Analytical Methods in Geomechanics 6(1): 129-139.
Seed, H.B. & Idriss, I.M. 1970. Soil moduli and damping factors for dynamic response analysis. Report
No. EERC 70-10. Earthquake Engineering Research Center, Berkeley, California.
Seed, H.B., Murarka, R., Lysmer, J. & Idriss, I.M. 1976. Relationships of maximum acceleration, maxi-
mum velocity, distance from source and local site conditions for moderately strong earthquakes, Bulle-
tin of the Seismological Society of America 66(4): 221-243.
Sloan, S.W. 1987. Substepping schemens for the numerical integration of elastoplastic stress-strain rela-
tions. International Journal for Numerical methods in Emgineering 24: 893-911.
Sloan, S.W., Abbo, A.J., Sheng, D. 2001. Refined explicit integration of elastoplastic models with auto-
matic error control. Engineering Computations 18(1/2): 121-154.
Taiebat, M. & Dafalias, Y.F. 2007. SANISAND: simple anisotropic sand plasticity model. International
Journal for Numerical and Analytical Methods in Geomechanics 32(8): 915-948.
Verdugo, R. & Ishihara, K. 1996. The steady state of sand soils. Soils and Foundation 36(2): 81-92.
Zhang, J.F., Andrus, R.D. & Juang, C.H. 2005. Normalized shear modulus and material damping ratio re-
lationships. Journal of Geotechnical and Geoenvironmental Engineering 131(4): 453-464.
Zhao, J., Sheng, D., Sloan, S.W. 2005. Explicit stress integration of complex soil models. International
Journal for Numerical and Analytical Methods in geomechanics 29(12):1209-1229.

16

View publication stats