H O S T E D  B Y
A stress–strain description of saturated sand under undrained cyclic torsional
increments in volumetric strain due to dilatancy, a normalized stress–plastic shear strain relationship is employed in combination with a novel empirical
Earlier experimental attempts to study the liquefaction
behavior of soils date back to the 1960s when Seed and Lee
(1966) conducted a series of undrained cyclic triaxial tests on
The Japanese Geotechnical Society
www.s
journal homepage: w
Soils an
Soils and Foundations 2015;55(3):559–574saturated sand and reported that the onset of liquefaction was
primarily governed by the relative density of the sand, the
confining pressure, the stress or strain amplitude and the
http://dx.doi.org/10.1016/j.sandf.2015.04.008
0038-0806/& 2015 The Japanese Geotechnical Society. Production and hosting by Elsevier B.V. All rights reserved.
E-mail addresses: nalinds@uom.lk (L.I.N. De Silva),
koseki@iis.u-tokyo.ac.jp (J. Koseki), chiaroga@iis.u-tokyo.ac.jp (G. Chiaro),
tsato@iis.u-tokyo.ac.jp (T. Sato).
Peer review under responsibility of The Japanese Geotechnical Society.1. IntroductionnCorresponding author.stress–dilatancy relationship derived for torsional shear loading. The proposed stress–dilatancy relationship accounts for the effects of over-consolidation
during cyclic loading. Numerical simulations show that the proposed model can satisfactorily simulate the generation of excess pore water pressure and
the stress–strain relationship of saturated Toyoura sand specimens subjected to undrained cyclic torsional shear loading. It is found that the liquefaction
resistance of loose Toyoura sand specimens can be accurately predicted by the model, while the liquefaction resistance of dense Toyoura sand
specimens may be slightly underestimated. (i.e., the liquefaction potential is higher). Yet, the model predictions are conservative.
& 2015 The Japanese Geotechnical Society. Production and hosting by Elsevier B.V. All rights reserved.
Keywords: Liquefaction behavior; Stress–dilatancy relationship; Over-consolidation; Sand; Constitutive model
JEL classification: IGC; A10; D07; E08; T06shear loading
Laddu Indika Nalin De Silvaa,b,n, Junichi Kosekic, Gabriele Chiaroc, Takeshi Satod
aDepartment of Civil Engineering, University of Moratuwa, Sri Lanka
bDepartment of Civil Engineering, University of Tokyo, Japan
cInstitute of Industrial Science, University of Tokyo, Japan
dIntegrated Geotechnology Institute Ltd., Japan
Received 21 January 2014; received in revised form 11 December 2014; accepted 2 February 2015
Available online 16 May 2015
Abstract
A constitutive model to describe the cyclic undrained behavior of saturated sand is presented. The increments in volumetric strain during undrained
loading, which are equal to zero, are assumed to consist of increments due to dilatancy and increments due to consolidation/swelling. This assumption
enables the proposed model to evaluate increments in volumetric strain due to dilatancy as mirror images of increments in volumetric strain due to
consolidation/swelling, thus simulating the generation of excess pore water pressure (i.e., reduction in mean effective principal stress) during undrained
cyclic shear loading. Based on the results of drained tests, the increments in volumetric strain due to consolidation/swelling are evaluated by assuming
that the quasi-elastic bulk modulus can be expressed as a unique function of the mean effective principal stress. On the other hand, in evaluating theciencedirect.com
ww.elsevier.com/locate/sandf
d Foundations
d Fnumber of loading cycles. Since then, extensive studies have
been conducted on soil liquefaction throughout the world
(Vaid and Thomas, 1995, among others) and a number of
attempts have been made to define proper constitutive models
to describe it (Liou et al., 1977; Liyanapathirana and Poulos,
2002, among others).
Based on the results of several series of experiments on saturated
Nomenclature
τzθ shear stress
σ
0
z; σ
0
r and σ
0
θ axial, radial and circumferential effective
stress, respectively
p0 mean effective stress
Drini relative density measured at confining pressure of
30 kPa
(τzθ=p0)max maximum shear stress ratio
τzθ max peak shear stress
Gzθ0 initial quasi-elastic shear modulus (¼dτzθ=dγezθ)
γzθ; γ
e
zθ; γ
p
zθ total, elastic and plastic shear strain, respec-
tively (engineering strain)
εpvol plastic volumetric strain
γref reference shear strain (¼ (τzθ=p0)/(Gzθ0=p0))
m, n, k material parameters that accounts for the stress
induced anisotropy of Young's moduli, shear
moduli and Poisson's ratio, respectively
CE, CG factors that account for the degradation of quasi-
elastic Young's and shear moduli, respectively
(assumed as zero in the present study)
L.I.N. De Silva et al. / Soils an560hollow cylindrical sand specimens, Towhata and Ishihara (1985a)
proposed a unique correlation between the shear work and the
generation of pore water pressure (PWP). Furthermore, the effects
of the rotation of the principal stress axes on sand liquefaction were
investigated by Towhata and Ishihara (1985b) using hollow
cylindrical specimens subjected to cyclic torsional shear loading.
However, compared to the large amount of experimental data
existing on liquefaction and the undrained behavior of soils, very
few models are available to successfully simulate the soil
performance under cyclic undrained loading. Ishihara et al.
(1975) proposed a model based on five postulates to trace the
generation of the excess PWP of sand subjected to undrained
irregular cyclic loading. This model qualitatively simulates the
stress–strain relationships and the shear stress versus mean effective
stress relationships.
A constitutive model to simulate the cyclic undrained behavior
of sand, based on the multi-spring concept, was developed by Iai
et al. (1992). In this model, commonly known as the “Towhata–Iai
model”, shear deformation is modeled by employing the multi-
spring concept, and the generation of excess PWP is modeled
using a unique correlation between the increments in excess PWP
and shear work, as proposed by Towhata and Ishihara (1985a).
Nishimura (2002) and Nishimura and Towhata (2004) modified
the above model by expanding the multi-spring concept from two
dimensions to three dimensions, while using an empirical stress–
dilatancy relationship to model the generation of excess PWP bycorrelating the stress–dilatancy relationship to consolidation. Never-
theless, these models do not consider the inherent anisotropy of
soils. Furthermore, the steady state during liquefaction and the
continuous increase in shear strain with cyclic loading cannot be
properly simulated.
An elasto-plastic constitutive model for sand, based on
the non-linear kinematic hardening rule, was employed to inves-
A Ez0/Eθ0, i.e. ratio of vertical to circumferential
quasi elastic Young's moduli at isotropic stress
state
Y normalized shear stress (¼ (τzθ=p0)/(τzθ=p0)max)
X normalized shear strain (¼γpzθ=γref )
D1 and D2 drag parameters
D plastic shear moduli immediately after reversal of
stress/initial plastic shear moduli (i.e., damage
parameter)
Dult minimum value for D
S amount of hardening
Sult maximum value for S
OC over-consolidation ratio
dεpvol=dγpzθ dilatancy ratio
Rk gradient of the empirical stress–dilatancy
relationship
Rm the maximum value for Rk
C intercept of the empirical stress–dilatancy
relationship
Cmin minimum value for C
oundations 55 (2015) 559–574tigate the effectiveness of the cement-mixing column method and
the gravel drain method as countermeasures against liquefaction by
means of a two-dimensional liquefaction analysis (Oka et al.,
1992). Later, Oka et al. (1999) further modified this model by
introducing a stress–dilatancy relationship that accounts for the
damage to plastic stiffness at large levels of shear strain. In
addition, several other constitutive models, based on the critical
state framework, are proposed in the literature. Jefferies (1993)
proposed a strain-hardening model, which utilizes the state
parameter, to explain the behavior of very loose to very dense
sand. A unified generalized plasticity model, based on the non-
linear critical state line, was proposed by Ling and Yang (2006).
It should be noted that all the above-described models are based
on either the critical state soil mechanics approach (e.g., Oka et al.,
1992) or the energetic approach (Iai et al., 1992; Nishimura,
2002). In the current study a different and original approach is
attempted by extending empirical relationships that are found to be
reasonably consistent with the experimental observations. The
undrained cyclic behavior of sand is simulated based on the
response whereby the same sand is shown during drained cyclic
loading. In fact, after appropriate normalization, the stress–strain
relationship is found to be unique for drained and undrained
conditions. Moreover, the generation of PWP during undrained
loading can be described based on the volumetric strain response
of sand during drained loading. This is done by improving the
model proposed by De Silva and Koseki (2012) that can
accurately simulate the drained cyclic behavior of sand, i.e., the
stress–strain relationship and the volumetric strain response.
However, no attempt has been made so far to utilize the above-
mentioned approach in simulating the cyclic undrained behavior (i.
e., liquefaction) of soil. The attempt is made in this paper, where a
cyclic constitutive model is presented to describe the undrained
cyclic behavior of sand. In the model, a simulation of the plastic
volumetric strain due to dilatancy (dεdvol) is combined with the
consolidation/swelling behavior of sand to simulate the generation
of excess PWP (i.e., a reduction in mean effective principal stress)
and stress–shear strain relationships.
2. Test material and procedures
In order to support the modeling work, a series of drained and
undrained cyclic torsional shear loading tests were conducted on
saturated Toyoura sand specimens (D50¼0.162 mm, emax¼0.966,
emin¼0.600, coefficient of uniformity Uc¼1.46). Hollow cylind-
rical specimens, having dimensions of 20 cm in outer diameter,
12 cm in inner diameter and 30 cm in height, were prepared at
initial relative densities (Drini) of 21%, 56% and 75%, as measured
at a confining pressure of 30 kPa. A modified air pluviation
technique, in which the sand pluviation was completed in a radial
Tokyo, was used for the testing program. A recently developed
local deformation measurement technique was employed in the
evaluation of quasi-elastic deformation properties and volumetric
strain during isotropic loading for drained specimens in the current
study. Refer to De Silva et al. (2005) for details on the torsional
shear apparatus and local deformation measurement system used.
In order to investigate the stress–dilatancy relationship of Toyoura
sand, a series of drained cyclic torsional shear tests was conducted
on loose and dense specimens, with Drini¼56% and 75%,
respectively, while keeping the mean effective principal stress
(p0) constant at 100 kPa. Details of the stress paths employed in the
drained cyclic tests are presented in De Silva and Koseki (2012). In
addition to the above, a series of constant stress amplitude
undrained cyclic torsional shear tests was conducted on Toyoura
sand specimens, while keeping the specimen height constant, for a
comparison with the model predictions. The stress paths of the
cyclic undrained tests are shown in Table 1.
3. Framework for modeling of liquefaction behavior
Changes in p0 during undrained loading cause the consoli-
dation/swelling of a specimen, while changes in shear stress (τ)
(1) IC (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 50-100 kPa)
(2) UTS (τzθ ¼ 0-22-22-0 kPa, until liquefaction at p΄0¼100 kPa)
SAT 31 79.4
L.I.N. De Silva et al. / Soils and Foundations 55 (2015) 559–574 561SAT 32 75.2
SAT 34 50.3
SAT 33 21.3
΄direction while slowly moving the nozzle of the pluviator in
alternate clockwise/anticlockwise directions, was employed in the
current study to minimize the degree of anisotropy of the horizontal
bedding plane of the hollow cylindrical specimens (see details in
De Silva et al. (2006)).
A high-capacity medium-sized hollow cylinder apparatus,
developed at the Institute of Industrial Science, University of
Table 1
Stress paths and test conditions of liquefaction tests.
Test Dr (%)
SAT 38 75.7
SAT 28 74.6Dr: relative density (%) measured at an isotropic stress state of σc¼30 kPa, IC: iso
mean effective stress at start of cyclic undrained loading.(1) IC (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 50-100 kPa)
(2) UTS (τzθ ¼ 0-30-30-0 kPa, until liquefaction at p΄0¼100 kPa)
(1) IC (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 50-100 kPa)
(2) UTS (τzθ ¼ 0-40-40-0 kPa, until liquefaction at p΄0¼100 kPa)
(1) IC (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 50-100 kPa)
(2) UTS (τzθ ¼ 0-60-60-0 kPa, until liquefaction at p΄0¼100 kPa)
(1) IC (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 50-100 kPa)
(2) UTS (τzθ ¼ 0-30-30-0 kPa, until liquefaction at p΄0¼100 kPa)
(1) IC (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 50-100 kPa)
(2) UTS (τzθ ¼ 0-30 kPa, flow failure at p΄0¼100 kPa)
0cause dilation. Therefore, the increments in volumetric strain
(dεvol) during undrained loading, which are equal to zero, are
assumed to consist of volumetric strain components due to
both dilatancy (dεdvol) and consolidation/swelling (dε
c
vol), as
expressed in Eq. (1).
dεvol ¼ dεcvolþdεdvol ¼ 0: ð1Þ
Stress pathstropic consolidation. UTS: undrained cyclic torsional shear loading, p0¼Initial
d FL.I.N. De Silva et al. / Soils an5624. Evaluation of dεcvol
Fig. 1(a) shows volumetric strain (εcvol) versus p
0 during the
swelling of isotropically consolidated Toyoura sand specimens. It
can be clearly seen that specimens of similar density have similar
swelling curves. The εcvol values reported in Fig. 1(a) are evaluated
by employing the local deformation measurement (De Silva et al.,
2005) and assuming the isotropy of the horizontal bedding plane
Fig. 1. Comparison of swelling curves for (a) loose and dense sand Toyoura
sand specimens, (b) 10 isotropic cycles between p0 ¼100 and 400 kPa of a
typical loose specimen and (c) 1st and 10th cycles of a typical loose specimen.(i.e., radial and circumferential strains, εr and εθ, respectively are
equal). Local deformation measurement transducers (LDTs) are
mounted on metal hinges, which are glued onto the membrane.
However, excessive hinge deformation may take place when the
confining stress becomes less than 50 kPa, causing an error when
evaluating εcvol, as shown in the upper left corner of Fig. 1(a).
In Fig. 1(a), it can also be observed that the swelling
curves of Toyoura sand, subjected to isotropic unloading and
reloading cycles in the range of p0 from 100 to 400 kPa, can be
expressed by Eq. (2).
dεcvol ¼
dp0
Ko
p0
p
0
0

 mk ð2Þ
where Ko is the bulk modulus at the reference mean effective stress
(p
0
0) and mk is a material parameter. As shown in Fig. 1(a), it is
found that the value for Ko (evaluated at p
0
0¼100 kPa) is 58 MPa
for dense Toyoura sand specimens (Dr¼75–80%) and 50 MPa for
loose Toyoura sand specimens (Dr¼54–57%). The value for mk,
for both dense and loose specimens, is taken as 0.9.
The effect of large isotropic loading/unloading cycles
(IC cycles) on the swelling curve is shown in Fig. 1(b).
Swelling curves for the first and tenth cycles are compared in
Fig. 1(c). No significant effect on the swelling curve due to the
application of isotropic cycles could be observed after applying
10 IC cycles (disregarding the effects of creep during the
application of 10 cycles). Therefore, a unique swelling curve is
employed in the current study to model εcvol.
5. Modeling of monotonic stress–shear strain relationship
It is a well-known fact that τzθ=p0 versus plastic shear strain
(γpzθ) during drained or undrained monotonic shear is of a non-
linear shape (e.g., Koseki et al., 1998; among others). A typical
relationship for a cyclic undrained test, conducted on a Toyoura
sand specimen of Drini¼75.7%, is shown in Fig. 2(a). γpzθ is
evaluated by deducting the elastic strain (γezθ) components from the
total shear strain (γzθ). The γ
e
zθ component is evaluated by
employing the quasi-elastic constitutive model (IIS model) pro-
posed by HongNam and Koseki (2005, 2008). The parameters
employed in the IIS model for a dense specimen with Drini of
about 75% (defined at a confining stress of 30 kPa) are listed in
Table 2.
It should be mentioned that, when the soil reaches the full
liquefaction state (i.e., p0E0), the τzθ=p0 values become extremely
sensitive to very small changes in p0. Hence, highly scattered data
can be observed. Koseki et al. (2005) investigated the liquefaction
properties of sand under low confining stress levels and proposed a
simplified procedure to estimate the liquefaction resistance by
introducing the concept of the apparent increase in effective mean
principal stress (Δp0), due to particle interlocking, as well as
parameter Δτ to correct τ due to possible measurement errors.
Therefore, a modified stress ratio was proposed by Koseki et al.
(2005): ðτzθΔτzθÞ=ðp0 þΔp0Þ. Chiaro et al. (2013) showed that
the values for Δτzθ and Δp0 may slightly vary during cyclic
oundations 55 (2015) 559–574loading. However, for simplicity, the values for Δτzθ and Δp0 were
assumed to be constant, as originally suggested by Koseki et al.
(2005) in the present study. The modified stress ratio is employed
in Fig. 2(a) with Δτzθ¼1.3 kPa and Δp0 ¼0.2 kPa. Typical
evaluations of Δτzθ and Δp0 are shown in Fig. 2(b). It can be
clearly seen that ðτzθΔτzθÞ=ðp0 þΔp0Þ versus γpzθ is of a non-
linear shape with hysteresis, which is similar to τzθ=p0 versus γ
p
zθ of
a drained cyclic test.
In order to quantitatively investigate the above, comparisons
of τzθ=p0 versus γ
p
zθ for the virgin loading of drained and
undrained specimens of similar densities are shown in Fig. 3. It
can be clearly seen that, for similar densities, τzθ=p0 versus γ
p
zθ
for the undrained test is very similar to that for the drained test.
This observation suggests that it is possible to evaluate dεdvol
(increment in volumetric strain due to dilatancy) during cyclic
undrained loading by combining the simulation of the stress–
shear strain relationship of a drained test with an appropriate
stress–dilatancy relationship.
Normalized stress–strain relationships during the virgin
loading (backbone curves) of undrained tests with different
densities are compared in Fig. 4. The peak shear stress ratio
((τzθ=p0)max) and the initial quasi-elastic shear modulus (Gzθ0)
at p
0
0¼100 kPa for a specimen of Drini¼75% were taken as
-1.6 -1.2 -0.8 -0.4 0.0 0.4 0.8 1.2
-0.8
-0.6
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
Δτ
Δ
(τ
zθ
-Δ
τ z
θ)/
(p
'+
Δ
p'
)
γpzθ (%)
-0.8 -0.4 0.0 0.4 0.8 1.2 1.6 2.0
-0.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Δτ = 1.3 kPa
τ z
θ
(k
Pa
)
p' (kPa)
Δp' = 0.2 kPa
Fig. 2. (a) ðτzθΔτzθÞ=ðp0 þΔp0Þ versus γpzθ relationship of undrained test and
(b) evaluation of Δτzθ and Δp0.
Table 2
Model parameters.
Test Quasi-elastic model parametersa
ram
L.I.N. De Silva et al. / Soils and Foundations 55 (2015) 559–574 563SAT 38 Ezo¼215 MPa, σ΄o¼100 kPa, νzθ0¼0.18, m¼0.5, n¼0.5,
k¼0.3, CE¼CG¼0.0, a¼0.7SAT 28
SAT 31
SAT 32
SAT 34 Ezo¼190 MPa, σ΄o¼100 kPa, νzθ0¼0.18, m¼0.5, n¼0.5, k¼0.3,
CE¼CG¼0.0, a¼0.7
SAT 33 Ezo¼145 MPa, σ΄o¼100 kPa, νzθ0¼0.18, m¼0.5, n¼0.5, k¼0.3,
CE¼CG¼0.0, a¼0.7
aDetails of the formulation of the quasi-elastic model (IIS Model) and its pa
bDrag and hardening parameters are the same as those employed in the simulati
loading in De Silva and Koseki (2012).
cDamage parameter Dult, for drained cyclic loading, is taken as 0.20.85 and 100 MPa, respectively, based on the results of a
drained test on a specimen of similar density. For Drini¼50%,
the (τzθ=p0)max and Gzθ0 values are 0.78 and 60 MPa, respec-
tively, and for Drini¼21%, the (τzθ=p0)max and Gzθ0 values are
0.60 and 42.9 MPa, respectively. Reference shear strain γref is
taken as the ratio of (τzθ=p0)max to Gzθ0=p
0
0.
Then, the Generalized Hyperbolic Equation (GHE), as
proposed by Tatsuoka and Shibuya (1991) (Eq. (3)), was
employed to simulate the backbone curve shown in Fig. 4.
Y ¼ X
1
C1ðXÞ þ
jXj
C2ðXÞ
ð3Þ
where
C1ðXÞ ¼
C1ð0ÞþC1ð1Þ
2
þ C1ð0ÞC1ð1Þ
2
cos
π
α0
X
 mt þ1
 !
ð4aÞ
C2ðXÞ ¼
C2ð0ÞþC2ð1Þ
2
þ C2ð0ÞC2ð1Þ
2
cos
π
β0
X

 nt þ1
0
B@
1
CA
ð4bÞ
and X and Y are normalized plastic shear strain and shear stress
parameters, respectively. Normalized stress and strain para-
meters, as defined below, are selected by following the
Drag parametersb Hardening parameterb, Sult Damage parameter
c, Dult
D1¼0.15 1.15 0.6
D2¼12
D1¼0.01
D2¼3.13
not required (specimen fails by flow failure)
eters are presented in HongNam and Koseki (2005, 2008)
on of the stress–shear strain relationship during drained cyclic torsional shear
Tatsuoka et al. (1997) reported that the stress–strain relationships
d F0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Drained 
backbone 
(OC = 1)
γpzθ (%) 
Drained loading
 SAT35(Drini = 79.3 %, OC = 1)
Undrained loading
 SAT32(Drini = 75.2 %, OC = 1)
Toyoura sand (saturated)
(σ'z = σ
'
r = σ
'
θ)initial = 100 kPa
τ z
θ/p
'
Undrained backbone
Fig. 3. Drained and undrained τzθ=p0 versus γ
p relationships during virgin
L.I.N. De Silva et al. / Soils an564procedure proposed by De Silva and Koseki (2012):
X ¼ γ
p
zθ
γref
and Y ¼ τzθ=p
0
ðτzθ=p0Þmax
ð5Þ
where γpzθ ¼ γzθγezθ and γref ¼ ðτzθ=p0Þmax=ðGzθ0=p0Þ.
The GHE has 8 parameters, i.e., C1(0), C1(1), C2(0),
C2(1), α0, β0, mt and nt, which can be determined from a
single monotonic drained torsional shear test. Refer to
Tatsuoka and Shibuya (1991) for the procedure for evaluating
these parameters.
It can be seen from Fig. 4 that the normalized stress–shear
strain relationships of undrained tests with different densities
are similar. Hence, they can be modeled by employing a single
set of GHE parameters obtained either by a normally con-
solidated drained test or an undrained shearing test. In order to
obtain a better fitting to the experiment data, GHE parameters
of undrained tests are selected by slightly modifying those of
drained tests (Fig. 4).
-5 0 5 10 15 20 25 30 35 40
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Simulation using GHE parameters  
that are slightly modified from SAT30
 SAT32 (Drini = 75.2 %)
 SAT28 (Drini = 74.6 %)
 SAT34 (Drini = 50.3 %)
 SAT33 (Drini = 21.3 %)
Toyoura sand(saturated)
p'o= 100 kPa
γref = (τz /p')max/ (Gz /p' )Y=
 ( τ
zθ
/p
')/
( τ
z θ
/p
') m
ax
X = γpzθ/γref
Simulation using GHE parameters 
of a drained test (SAT30, OC =1)
GHE parameters (modified from SAT 30)
C1(0) = 35, C1(∞) = 0.3, α = 0.000146, m = 0.25
C2(0) = 0.078, C2(∞) = 1, β = 0.178, n = 0.25
Fig. 4. Comparison of the backbone curves of undrained tests with different
densities.
zθ
loading.of soils are significantly influenced by cyclic strain hardening,
damage due to straining etc., which are caused by the rearrange-
ment of particles during cyclic loading, and proposed additional
rules to model these features to account for behaviors under more
general stress conditions. In this regard, a conceptual approach was
implemented for dense Toyoura sand under a plane strain condition
by applying a horizontal shift to the basic skeleton curves, i.e.,
dragging the basic skeleton curve along the X-axis (strain parameter
axis) in the opposite direction to its loading direction, while
applying the (extended) Masing's rule (refer to Masuda et al.
(1999) and Tatsuoka et al. (2003) for further details). It was
assumed that the amount of drag β, applied to one basic skeleton
curve in one loading direction, is a function of the plastic shear
strain accumulated in the opposite loading direction (Masuda et al.,
1999; Balakrishnaiyer and Koseki, 2000; Tatsuoka et al., 2003;
HongNam and Koseki, 2008). The same approach was employed
to model the cyclic stress–strain relationships of Toyoura sand
under cyclic torsional shear loading (HongNam and Koseki, 2008;
De Silva and Koseki, 2012). The dragged backbone curve can be
written as follows:
Y ¼ ðXβÞ
1
C1ðXβÞ þ
jXβj
C2ðXβÞ
ð6Þ
where β denotes the amount of drag, which can be evaluated by a
drag function as shown below.
β ¼ X
0
1
D1
þ X0D2
ð7Þ
where D1 is a fitting parameter (i.e., initial gradient of the drag
function) and D2 is the maximum amount of drag. D1 and D2 can
be experimentally determined. X0 is the accumulated normalized
plastic shear strain in one direction (positive or negative direction).
HongNam and Koseki (2008) employed D1¼0.45 and
D2¼3.13 for dense Toyoura sand (Drini¼71%) subjected to
cyclic torsional shear loading starting from an isotropic stress
state (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 100 kPa). However, it was observed by
De Silva and Koseki (2012) that the application of drag alone
is not sufficient for simulating the cyclic stress–shear strain
relationship close to the peak stress of the material.
6.1. Modification of extended Masing's rule
In view of the limitation of Masing's rule with drag in
simulating the stress–shear strain relationship close to the peak
stress of the material, De Silva and Koseki (2012) proposed
two conceptual modification factors, which take into account6. Modeling of cyclic stress–shear strain relationship by
using extended Masing's rule
Subsequent unloading/reloading cycles are modeled by
employing the procedure proposed by De Silva and Koseki
(2012), as briefed below.
oundations 55 (2015) 559–574the hardening behavior during cyclic loading (reduction in
damping ratio with constant stress amplitude cyclic loading)
to dilative at the phase transformation state (usually when
τzθ=p'E0.5 in the case of Toyoura sand).
In addition, De Silva and Koseki (2012) proposed a
conceptual equation for hardening parameter S by assuming
that S can be expressed as a hyperbolic function of the total
normalized plastic strain up to the current turning point, as
follows:
S¼ 1þ
PΔX
 
Upto current turing point
D2
D1
þ
PΔX
 
Upto current turing point
Sult1ð Þ
ð11Þ
where Sult is the maximum value for S after applying an infinite
number of cycles, and D1 and D2 are the same drag parameters
Toyoura sand with Drini¼50% (De Silva and Koseki, 2012).
d Foundations 55 (2015) 559–574 565and damage to plastic stiffness (dτzθ=dγ
p
zθ) at large stress
levels, while maintaining continuity in the simulation.
The two parameters in the GHE, C1(X¼0) and C2(X¼1),
represent the initial plastic stiffness and the peak strength,
respectively. Therefore, the damage occurring to the plastic
stiffness can be obtained by multiplying C1(X) by damage
factor D, and the hardening can be obtained by multiplying
C2(X) by hardening factor S. Note that in the approach used by
Masuda et al. (1999), Balakrishnaiyer and Koseki (2000), and
Tatsuoka et al. (2003), unique backbone curves were used to
model subsequent cyclic branches by employing extended
Masing's rules.
Therefore, the hysteresis curve, starting from an arbitrary
point A, can be obtained by employing the extended Masing's
rules with damage and hardening by using Eq. (8).
Y ¼ YAþ
XXA
1
C1
X  XA
np

 
D
þ jXXAj
npC2 X  XAnp

 
S
ð8Þ
In order to maintain continuity, the dragged backbone curve,
which is in the same direction as the curve to be modeled,
should be modified, as shown in Eq. (9), to determine the
proportional parameter np by using the extended Masing's
rules. Refer to Masuda et al. (1999) and Tatsuoka et al. (2003)
for details on the extended Masing's rules and sub-rules.
Y ¼ ðXβÞ
1
C1ðXβÞD þ
jXβj
C2ðXβÞS
ð9Þ
Note that parameters D and S were employed in Eqs. (8)
and (9) in such a way that D and S are constant for a given
hysteresis curve, but they change from curve to curve. An
evaluation of D and S for a particular hysteresis curve is
made based on a few empirical equations, as follows. The
plastic shear modulus (D) can be expressed as an “S curve”,
as proposed by De Silva and Koseki (2012) and shown in
Eq. (10).
D¼Dultþ
1þexp γnð Þ
1þexp
Δγpzθ
p
γn
  1Dultð Þ ð10Þ
where jΔγpzθjp is the total plastic shear strain (%) accumu-
lated between the current and the previous turning points,
and Dult is the minimum value for D, which would be
applied to evaluate the minimum value for the plastic shear
modulus. γn corresponds to the value for γpzθ (in %) at which
the volumetric behavior of the material changes from
contractive to dilative, and is taken as 0.8.
De Silva and Koseki (2012) proposed Dult¼0.2 for drained
cyclic torsional shear loading, although experimental evidence
suggests that Dult may lie between 0.2 and 0.6. In the current
study, for an improved simulation of the experimental results,
Dult¼0.6 is used for the case of cyclic undrained loading. Yet,
it should be noted that the D value is assumed to be equal to
L.I.N. De Silva et al. / Soils an1.0 (i.e., no damage to the plastic shear modulus) until the
volumetric behavior of the material changes from contractiveThe same values are employed for the respective densities in the
simulation of undrained behavior. The concepts of drag, hard-
ening and damage during drained cyclic torsional shear loading
are illustrated in Fig. 5.
7. Evaluation of dεdvol
In order to evaluate dεdvol, it is necessary to combine dγ
p
zθ
with an appropriate stress–dilatancy relationship. For the above
purpose, the stress–dilatancy relationship, proposed by De
Silva and Koseki (2012), is modified as follows.
During undrained cyclic loading, the mean effective stress
(p0) mainly decreases with the number of cycles. It is assumed
in this study that the above reduction in p0 is associated with
-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0
-80
-60
-40
-20
0
20
40
60
80
τ
γ  (%) used in Eq. (7). As suggested by De Silva and Koseki (2012),
Sult¼1.15 was employed in the current study.
Note that the values for D1 and D2 differ for the cases with
and without damage and hardening. HongNam and Koseki
(2008) proposed D1¼0.45 and D2¼3.13 for Toyoura sand
(Drini¼71%) subjected to cyclic torsional shear loading starting
from an isotropic stress state (σ΄z ¼ σ΄r ¼ σ΄θ ¼ 100 kPa). After
introducing the damage and hardening factors, D1¼0.15 and
D2¼12 were found to be appropriate for dense Toyoura sand
with Drini¼75%, and D1¼0.01 and D2¼3.13 was suggested forFig. 5. Illustration of the concepts of drag, hardening and damage using a
typical drained cyclic torsional shear test (De Silva and Koseki, 2012).
over-consolidation and cyclic mobility. Firstly, the soil under-
going a decrease in p0 is subjected to over-consolidation
until the stress state exceeds the phase transformation stress
state (Ishihara and Li, 1972) for the first time (i.e., the first
instance where the volumetric strain increment changes from
contractive to dilative behavior, i.e., dεdvolo0). Then, the soil
will enter the stage of cyclic mobility.
7.1. Stress–dilatancy relationship during virgin loading and
before exceeding the phase transformation stress state
It can be observed in Fig. 6(a) and (b) that the stress–dilatancy
relationship during cyclic loading before exceeding the phase
transformation stress state is different from that after exceeding the
phase transformation stress state (refer to the explanations given
within Fig. 6(a) and (b)). In addition, it can be observed from
Fig. 6(a) that the effects of over-consolidation significantly alter the
stress–dilatancy relationship during virgin loading (refer to the line
shown as “start of loading” in Fig. 6(a)). Since sand will be
subjected to over-consolidation/swelling during undrained cyclic
loading, the stress–dilatancy relationships during different stages
are addressed separately in the current study.
The stress–dilatancy relationships, employed in the evalua-
tion of dεdvol, consist of four equations, namely, Eqs. (A), (B),
(C) and (D) in Fig. 7(a). The stress–dilatancy model is
summarized in Table 3. According to the proposed model,
the stress–dilatancy relationship during the virgin loading of
normally consolidated Toyoura sand is given by Eq. (12) (refer
to Eq. (A) of Fig. 7(a)) with Rk¼1.3 and C¼0.6.
τzθ
p0
¼ Rk 
dεdvol
dγpzθ
 
7C for dγpzθ40 and dγ
p
zθo0; respectively
ð12Þ
Since Toyoura sand was subjected to undrained cyclic
torsional shear loading from a normally consolidated stress
state, dεdvol during virgin loading was evaluated in the current
study by applying the stress–dilatancy relationship for virgin
loading, as given by Eq. (A) in Fig. 7(a).
If the loading direction after virgin loading is reversed before
exceeding the phase transformation stress state, a different linear
relationship is employed, by referring to Fig. 6, to account for the
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
τ z
θ/p
'
-0.4
Rk= 1.5, C= -0.50C= -0.18
Cyclic loadings after exceeding 
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
D
C
BA
L.I.N. De Silva et al. / Soils and Foundations 55 (2015) 559–574566-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
-0.6 phase 
transformation 
stress state
-dεdvol/dγ
p
zθ
the phase transformation stress 
state for the first time
(-0.50 ≥ τz /p' ≥ 0.50)
Fig. 6. Stress–dilatancy relationships during cyclic torsional shear loading (a)
before exceeding the phase transformation stress state for the first time and (b)
bilinear stress–dilatancy relationship after exceeding the phase transformation-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
-0.8
-0.6
-dεdvol/dγ
p
zθ
-0.2
0.0
0.2
0.4
0.6
0.8
phase transformation 
stress state
Subsequent unloading
Subsequent 
unloading
Rk= 0.33, 
τ z
θ/p
'
Subsequent loading
Rk= 1.5, C= 0.50
Subsequent unloading
Rk= 0.33, C= 0.18
Start of loadingstress state.
Adopted from De Silva et al. (2014).-0.8 -0.4 0.0 0.4 0.8 1.2
-0.8
-0.6
-0.8 -0.4 0.0 0.4 0.8 1.2
-0.8
-0.6
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
Fig. 7. Proposed stress–dilatancy relationships for (a) different stress states
during undrained cyclic torsional shear loading of dense sand and (b) the stress
state within the over-consolidation boundary surface.
Adopted from De Silva et al. (2014).
.. (A
den
ely

Rk
) and
cu
d Fohardening behavior until the stress state exceeds the phase
transformation stress state for the first time (refer to steep dashed
lines in Fig. 6(a) and (b)). Refer to the stress–dilatancy relationship
denoted by Eq. (B) in Fig. 7(a) for further details. Eq. (B), as
shown in Fig. 7(a), is obtained by employing Rk¼2.2 and C¼0.5
in Eq. (12) for both dense and loose Toyoura sand (Drini¼75% or
Drini¼50%).
7.2. Effects of over-consolidation on stress–dilatancy
relationship
The over-consolidation ratio (OC) continuously changes (i.
e., increases) when the cyclic loading continues within the
Table 3
Stress–dilatancy relationships to evaluate dεdvol during undrained cyclic loading
Simulation case Stress–dilatancy relationships
Virgin loadinga Subsequent loading
Case 1 (Bilinear) Eq. (A) with
Rk ¼ 1:3; C¼ 0:6
Eq. (A) with Rk ¼ 1:5;C¼ 0:5 and Eq
Case 2 (Modified
bilinear)
Eq. (A) with
Rk ¼ 1:3; C¼ 0:6
If
τzθ=p0maxr C
Eq. (B) with Rk ¼ 2:2;C¼ 0:5for both
(Drini 50%) Toyoura sand, respectiv
Case 3 (Modified
bilinearþEffects
of OC)
Eq. (A) with
Rk ¼ 1:3; C¼ 0:6
If
τzθ=p0maxr Cτzθ=p04 C  lnðOCÞ τzθ=p0rC
Same as Case 2.1 Eq. (C) with
(Drini¼75%
respectively
OC: over-consolidation ratio.
aVirgin loading always starts from a normally consolidated stress state in the
specimen.
L.I.N. De Silva et al. / Soils anlimits of the phase transformation stress state (i.e., |τzθ=p0|
oC). Hence, the effects of the change in OC on the stress–
dilatancy relationship, within the above stress range, are taken
into account by applying Eq. (13) (Oka et al., 1999) (refer to
Eq. (C) in Fig. 7(a)). Note that the model by Oka et al. (1999)
was formulated using the stress and strain invariants, while
maintaining the objectivity. For simplicity, on the other hand,
the stress–dilatancy relationship formulated in this study is
described one-dimensionally using specific stress components.
As will be described at the end of Section 9, an attempt is
under way to extend the current stress–strain description into a
generalized three-dimensional modeling.
For dτzθ40 and dτzθo0:
 dε
p
vol
dγpzθ
 
¼Dk
τzθ
p0  τzθp0 =lnðOCÞ

 
Rk
0
@
1
A ð13Þ
where Dk ¼
 τzθp0 =ðC  lnðOCÞÞ
 1:5
The Dk value changes with OC. It is assumed that Eq. (13)
continues until Dk¼1 with Rk¼2.2 and C¼0.50 for both
dense and loose specimens and then follows Eq. (13) with
Dk¼1 (i.e., changes in Dk due to changes in the stress state arenot considered after Dk41). Refer to the stress–dilatancy
relationship denoted by Eq. (C) in Fig. 7(a) for details. The
stress state at which Dk¼1 is defined as the over-consolidation
boundary surface (Oka et al., 1999).
When Dk¼1
C¼ τzθ
p0
=lnðOCÞ ð14Þ
Then, Eq. (13) can be rewritten as follows:
 dε
p
vol
dγpzθ
 
¼
τzθ
p0 C
Rk
 !
ð15Þ
Eq. (15) corresponds to Eq. (B) in Fig. 7(a). After that, the
) with Rk ¼ 0:33; C¼ 0:18 for immediately after stress reversal
If
τzθ=p0max4C
se (Drini 75%) and loose
(Case 2.1)
Eq. (D) with Rmax ¼ 1:5; Cmin ¼ 0:36and
Eq. (A) with Rk ¼ 0:33;C¼ 0:18 for
immediately after stress reversal
Eq. (D) with Rmax ¼ 1:5; Cmin ¼ 0:36and
Eq. (A) with Rk ¼ 0:33;C¼ 0:18for
immediately after stress reversallnðOCÞ

¼ 2:2; C¼ 0:5 for both dense
loose (Drini¼50%) specimens,
rrent study
τzθ=p0max: the maximum value of τzθ=p0 currently applied to the
undations 55 (2015) 559–574 567stress–dilatancy relationship follows Eq. (15) until the stress
state exceeds the phase transformation stress state for the first
time and then follows the modified bilinear stress–dilatancy
relationship (Eq. (D) shown in Fig. 7(a)).
In short, if stress reversal occurs before exceeding the phase
transformation stress state, the stress–dilatancy relationship
follows a combination of Eqs. (B) and (C) until the stress state
exceeds the phase transformation stress state for the first time.
The combination of Eqs. (B) and (C), for dτzθ40 with
different values for OC, is illustrated in Fig. 7(b).
Rearranging the terms in Eq. (14), we get the over-
consolidation boundary surfaces for positive and negative
shear stress increments, as shown in Eqs. (16a) and (16b).
For dτzθ40:
τzθ ¼ C  p0  lnðOCÞ ð16aÞ
Similarly, for dτzθo0:
τzθ ¼ C  p0  lnðOCÞ ð16bÞ
Fig. 8 shows a typical over-consolidation boundary surface
for a normally consolidated specimen subjected to cyclic
undrained torsional shear loading starting from p
0
0¼100 kPa.
minimum value for C after application of a large number of
d Fconstant stress amplitude cyclic loadings (i.e., Cmin¼0.36
adopted from De Silva et al. (2014)) and D is the sameThe phase transformation stress state for the above specimen is
also indicated in Fig. 8.
7.3. Stress–dilatancy relationship after exceeding the phase
transformation stress state
If the loading direction during virgin loading is reversed
after exceeding the phase transformation stress state, the
stress–dilatancy relationship will follow the modified bilinear
stress–dilatancy relationship, as given below, for subsequent
cyclic loadings (refer to Eq. (D) in Fig. 7(a)).
τzθ
p0
¼ ðRmax  DÞ 
dεdvol
dγpzθ
 
7
Cmin
D
ð17Þ
where Rmax is the maximum value for Rk in Eq. (12) (i.e.,
Rmax¼1.5 adopted from De Silva et al. (2014)), Cmin is the
0 20 40 60 80 100 120
-40
-30
-20
-10
0
10
20
30
40
τ z
θ
(k
Pa
)
p' (kPa)
OC boundary surface
τz = × ×
Phase transformation 
line (τzθ/p' = C )
p'0
Fig. 8. Typical over-consolidation boundary surface and phase transformation
stress state for a normally consolidated specimen with no initial shear.
L.I.N. De Silva et al. / Soils an568damage parameter as in Eq. (10).
The following boundary conditions were specified for the
RmaxD and Cmin/D values by referring to the experimental
data (based on De Silva et al. (2014)). Note that, RmaxD
should not be less than 1.0. The maximum value for Cmin/D
should not be larger than 0.60. Therefore, the RmaxD value
in Eq. (2) varies between 1.5 and 1.0, while the Cmin/D value
varies between 0.36 and 0.60 depending on the accumulated
plastic strain between the current and the previous turning
points (i.e., damage parameter D).
The application of different stress–dilatancy equations for
different stress states during a typical cyclic undrained test is
illustrated in Fig. 9. The four-phase stress–dilatancy model,
employed in the current study, is summarized in Table 3.
By combining Eq. (1) with (2), the following equation can
be derived to evaluate the change in mean effective stress (dp0)
during undrained cyclic loading:
dp0 ¼ Ko
p0
p0o
 mk
ðdεdvolÞ ð18Þ8.1. Dense sand behavior
The experimentally obtained effective stress path (τzθ versus p0)
and stress–shear strain relationships (τzθ versus γzθ) of a dense
Toyoura sand specimen, subjected to undrained cyclic torsional
shear loading with stress amplitude (τzθ=p
0
o) equal to 0.22, are
shown in Fig. 10(a) and (a1), respectively. In order to show the
improvement in the simulation results with the introduction of the
modifications to the stress–dilatancy relationship, a simulation ofwhere dεdvol is the volumetric strain increment due to dilatancy,
which is evaluated by combining the simulation of stress–shear
strain relationship during drained cyclic torsional loading with
the proposed stress–dilatancy model. Then, by the numerical
integration of Eq. (18), the generation of excess PWP with the
shear stress level (or p0 versus τzθ) can be established. In
addition, the stress–strain relationship (i.e., τzθ versus γ
p
zθ) can
also be obtained.
8. Simulation of liquefaction behavior
0 20 40 60 80 100 120
-40
-30
-20
-10
0
10
20
30
40
D
C
B
B
τ z
θ
(k
Pa
)
p' (kPa)
A
Stress state exceeds 
the phase transformation 
state for the first time
application of Eqs. A, B, C, D for 
different stress-dilatancy regions
Fig. 9. Application of different stress–dilatancy relationships for different
stress states during undrained cyclic torsional shear loading.
oundations 55 (2015) 559–574the stress paths was carried out in three steps, as denoted by (b),
(c) and (d) in Fig. 10. The corresponding stress–strain relation-
ships are denoted by (b1), (c1) and (d1), respectively.
First, a unique bilinear stress–dilatancy relationship (refer to
De Silva et al. (2014) for details), without considering any
effects of over-consolidation or change in the stress–dilatancy
relationship with cyclic loading, is employed (refer to Case 1
in Table 3). The simulated stress path and the stress–strain
relationship that corresponds to each of the above stress–
dilatancy relationships are shown in Fig. 10(b) and (b1),
respectively. It can be seen that the reduction in the change
in p0, due to the effects of over-consolidation (see Fig. 10(a)),
cannot be accurately simulated by employing the above stress–
dilatancy relationship, as shown in Fig. 10(b). In addition,
dense Toyoura sand shows a continuous increase in shear
strain after the onset of cyclic mobility, as shown in Fig. 10
(a1), which is similar to that of a loose sand specimen, and
ends up yielding the same curve (closed loop) for subsequent
cycles, as shown in Fig. 10(b1).
-20 0 20 40 60 80 100 120
-30
-20
-10
0
10
20
30
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
-30
-20
-10
0
10
20
30
-20 0 20 40 60 80 100 120
-30
-20
-10
0
10
20
30
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
-30
-20
-10
0
10
20
30
-20 0 20 40 60 80 100 120
-30
-20
-10
0
10
20
30
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
-30
-20
-10
0
10
20
30
-20 0 20 40 60 80 100 120
-30
-20
-10
0
10
20
30
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
-30
-20
-10
0
10
20
30
Fig. 10. Simulations of stress paths and stress–strain relationships using different stress–dilatancy relationships (τzθ=p0 ¼0.22, Drini¼75.7%).
L.I.N. De Silva et al. / Soils and Foundations 55 (2015) 559–574 569
Fig. 11. Simulations of stress paths and stress–strain relationships using different stress–dilatancy relationships (τzθ=p0 ¼0.40, Drini¼79.4%).
L.I.N. De Silva et al. / Soils and Foundations 55 (2015) 559–574570
Consequently, the variations in the stress–dilatancy relation-
ship during cyclic loading are taken into account in the
simulation by employing the modified bilinear stress–dilatancy
relationship (refer to De Silva et al. (2014)) without consider-
ing any effects due to over-consolidation (refer to Case 2 in
Table 3). It can be seen from Fig. 10(c) that the stress path
after the onset of cyclic mobility is improved and, in addition,
the stress–strain relationship is improved to some extent
showing a continuous increase in shear strain with cyclic
loading after the onset of cyclic mobility, as shown in Fig. 10
(c1). However, when the stress path enters the steady state, the
stress–strain relationship becomes a closed loop. Furthermore,
since the effects of over-consolidation are not taken into
account in the above simulation, the reduction in the change
in p0 due to the effects of over-consolidation, as shown in
Fig. 10(a), cannot be accurately simulated.
Finally, the further modified stress–dilatancy model, sum-
marized in Case 3 (Table 3), which considers the effects of
over-consolidation and changes in the stress–dilatancy rela-
tionship during cyclic loading, was employed in the simulation
of the stress path and stress–strain relationships during
undrained cyclic loading, as shown in Fig. 10(d) and (d1),
respectively. It can be observed that the simulations for both
the stress path and the stress–strain relationship are certainly
-20 0 20 40 60 80 100 120
-40
-30
-20
-10
0
10
20
30
40 Toyoura sand (saturated)
Test SAT34
Drini =  50.3 %
τ z
θ
(k
Pa
)
p' (kPa)
Exper ime nt
OC boundary surface
phase transformation 
line
-10 -8 -6 -4 -2 0 2 4 6 8 10
-40
-30
-20
-10
0
10
20
30
40 Toyoura sand (saturated)
Test SAT34
Drini = 50.3 %
τ z
θ
(k
Pa
)
γzθ (%)
Experiment
-20 0 20 40 60 80 100 120
-40
-30
-20
-10
0
10
20
30
40
Test SAT34
D rini = 50.3 %
Bilinear stress-dilatancy 
model 
τ z
θ
(k
P
a)
p' (kPa)
Simulation
Phase 
transformation 
l ine
-10 -8 -6 -4 -2 0 2 4 6 8 10
-40
-30
-20
-10
0
10
20
30
40 Toyoura sand (saturated)
Test SAT34
Drini = 50.3 %
Bilinear stress-dilatancy 
model 
τ z
θ
(k
P
a)
γpzθ (%) 
Simulation
Toyoura sand (saturated)
L.I.N. De Silva et al. / Soils and Foundations 55 (2015) 559–574 571-20 0 20 40 60 80 100 120
-40
-30
-20
-10
0
10
20
30
40
Test SAT34
Drini =  50.3 %
Modified bilinear 
stress-dilatancy model 
τ z
θ
(k
P
a)
p' (kPa)
Simulation
Toyoura sand(saturated)Fig. 12. Simulations of stress paths and stress–strain relationships using-10 -8 -6 -4 -2 0 2 4 6 8 10
-40
-30
-20
-10
0
10
20
30
40 Toyoura sand (saturated)
Test SAT34
Drini = 50.3 %
Modified bilinear 
stress-dilatancy model 
γp (%) 
τ z
θ
(k
P
a)
Simulationzθ
different stress–dilatancy relationships (τzθ=p0 ¼0.30, Drini¼50.3%).
obtained from the experimental data is compared with that
obtained from the simulation results by employing different
stress–dilatancy relationships. However, it can be seen that the
simulation was significantly improved when the modified
bilinear stress–dilatancy relationship was employed while
considering the effects of over-consolidation. It can be seen
that the liquefaction resistance of dense Toyoura sand speci-
mens is underestimated (i.e., in the simulation liquefaction
occurs faster) by the simulated results. On the other hand, the
experimentally evaluated liquefaction resistance curves of the
loose Toyoura sand specimens (obtained using data from this
study and Chiaro et al., 2012) are similar to those obtained
from the simulation results after employing the modified
bilinear stress–dilatancy relationship with the effects of over-
consolidation (refer to Fig. 14(b)).
9. Conclusions
The following main conclusions can be derived from the
above study:
15a)
-1 0 1 2 3 4 5 6 7 8 9 10
0
5
10
γpzθ (%)
τ z
θ
(k
P
Experiment
Fig. 13. Comparison of (a) stress paths and (b) stress–strain relationships
(Drini¼21.3%).
d Fimproved after introducing the effects of over-consolidation
into the modified bilinear stress–dilatancy relationship. A
continuous increase in strain can be observed until the stress
path enters the steady state.
It should be noted that the simulation started producing the
same curve when jγpzθj was close to about 4%, while the
experimental data shows a continuous increase in shear strain
with cyclic loading. This strain level (jγpzθj¼4%) is close to the
strain level at which the peak stress state of dense Toyoura
sand is mobilized in drained shearing. Further modification of
the stress–strain relationship, which considers the strain-
softening behavior observed during cyclic torsional shear
loading (Kiyota et al., 2008), would be necessary to address
the above issue for dense sand. However, this was not
attempted in the current study. Since the simulation of the
stress–strain relationship of dense Toyoura sand subjected to
cyclic undrained loading using the bilinear stress–dilatancy
model (refer to Fig. 10(b1)) gives the same curve whenγpzθ 3%, liquefaction resistance is defined in the current
study as the number of cycles required to yield a double
amplitude shear strain of 6%.
A comparison of the experimental stress paths and the
stress–strain relationships of dense Toyoura sand specimens,
subjected to undrained cyclic torsional shear loading, with
their simulation results for stress amplitudes (τzθ=p
0
o) equal to
0.40, is shown in Fig. 11. The same stress–dilatancy relation-
ships as employed in Fig. 10(b)–(d) are utilized in Fig. 11.
8.2. Loose sand behavior
The experimental effective stress paths and stress–strain
relationships of a loose Toyoura sand specimen (Drini¼50%)
with stress amplitudes (τzθ=p
0
o) equal to 0.30 are shown in
Fig. 12(a) and (a1), respectively. It can be seen in Fig. 12(c)
and (c1) that the simulations of the stress path and the stress–
strain relationship, respectively, become consistent with the
corresponding experimental data when the modified bilinear
stress–dilatancy relationship is employed (compare Fig. 12(a)
and (a1) with Fig. 12(c) and (c1), respectively). The simulation
continues until the double amplitude of shear strain becomes
about 18% and then gives the same curve for further cycles.
Fig. 13(a) compares the experimentally obtained stress paths
of a very loose specimen (Drini¼21.3%) with its simulation,
while Fig. 13(b) compares the experimentally and numerically
obtained stress–strain relationships. Note that the specimen
shows flow failure and the simulation shows a similar
tendency. Therefore, only the stress–dilatancy relationship
during virgin loading (refer to Table 3) is required for the
simulation.
8.3. Liquefaction resistance curve
Fig. 14(a) and (b) shows the liquefaction resistance curves
for dense and loose Toyoura sand specimens, respectively. In
the current study, liquefaction resistance is defined as the
L.I.N. De Silva et al. / Soils an572number of cycles required to yield a double amplitude shear
strain of 6%. In each figure, the liquefaction resistance curve0 20 40 60 80 100 120
0
5
10
Bilinear stress-dilatancy 
model 
τ z
θ (
kP
p' (kPa)
Experiment
15
20
25
a)
Toyoura sand (saturated)
Test SAT33
Drini = 21.3 %
Bilinear stress-dilatancy 
model 
Simulation20
25
Toyoura sand (saturated)
Test SAT33
Drini = 21.3 %
Simulationounda(1)tions 55 (2015) 559–574A unique swelling curve, that does not change with the
number of cycles, has been proposed to evaluate the
(2)
(3)
(4)
Fig.
and
d Fo0.20
τ z0.25
0.30 Simulations
 Bilinear
 Modified bilinear
 Modified bilinear + OC
 Experiment
θ/p
' o0.1
0.2
0.3
2 4 6 8 10 20 40 60 80
Toyoura sand
Drini = 75 %
no.of cycles
τ
no.of cycles required to have DA = 6 %
0.350.4
0.5
0.6
Simulations
 Bilinear
 Modified bilinear
 Modified bilinear + OC
 Experiment
zθ
/p
' o
L.I.N. De Silva et al. / Soils anincrements in volume change due to consolidation/swelling
during undrained cyclic loading.
The normalized stress–shear strain relationships of undrained
tests with different densities were found to be similar. Hence,
they can be modeled by employing a single set of GHE
parameters obtained either by normally consolidated drained
tests or undrained shearing tests.
τzθ=p0 versus γ
p
zθ of an undrained specimen is very
similar to that of a normally consolidated drained
specimen of similar density. Hence, similar GHE para-
meters can be used for a drained normally consolidated
specimen to evaluate the increments in volume change
due to the dilatancy of an undrained specimen. In
addition, drag parameters, damage parameters and hard-
ening parameters were also similar for both drained and
undrained loadings.
The stress path during undrained cyclic loading is divided
into four sections, namely (1) virgin loading, (2) stress path
within the limits of phase transformation stress state, (3)
stress path within the limits of over-consolidation bound-
ary surfaces and (4) stress path after exceeding the phase
transformation stress state for the first time. Different
stress–dilatancy relationships are proposed for each section
of the stress path to evaluate the increments in volume
change due to dilatancy.
De
c
C
De
G
De S
c
S
0.10
0.15
0.6 0.81 2 4 6 8 10 20 40
Toyoura sand
Drini = 50 %
no.of cycles
no.of cycles required to have DA = 6 %
After Chiaro et al.(2012)
14. Liquefaction resistance curves: (a) dense Toyoura sand, Drini¼75%
(b) Loose Toyoura sand, Drini¼50%.1
De Sr
troceedings of the International Symposium on Geomechanics and
eotechnics of Particulate Media, Ube, Yamaguchi, Japan, pp. 29–33.
ilva, L.I.N., Koseki, J., Sato, T., Wang, L., 2005. High capacity hollow
ylinder apparatus with local strain measurements. Proceedings of the
econd Japan–U.S. Workshop on Testing, Modeling and Simulation, vol.
56. Geotechnical Special Publication, ASCE16–28.
ilva, L.I.N., Koseki, J., Wahyudi, S., Sato, T., 2014. Stress–dilatancyte
PSilva, L.I.N., Koseki, J., 2012. Modelling of sand behavior in drained
yclic shear. In: Miura (Ed.), Advances in Transportation Geotechnics II.
RC Press, pp. 686–691.
Silva, L.I.N., Koseki, J., Sato, T., 2006. Effects of different pluviation
chniques on deformation property of hollow cylinder sand specimens. In:References
Balakrishnaiyer, K., Koseki, J., 2000. Modelling of stress–strain relationships
of a reconstituted gravel subjected to large cyclic loading. In: Elnashai, A.
S., Antoniou, S. (Eds.), Implications of Recent Earthquakes on Seismic
Risk: Series on Innovation in Structures and Construction, vol. 2; 2000, pp.
105–114.
Chiaro, G., Kiyota, T., Koseki, J., 2013. Strain localization characteristics of
loose saturated Toyoura sand in undrained cyclic torsional shear tests with
initial static shear. Soils Found. 53 (1), 23–34.
Chiaro, G., Koseki, J., Sato, T., 2012. Effects of initial static shear on
liquefaction and large deformation properties of loose saturated Toyoura
sand in undrained cyclic torsional shear tests. Soils Found. 52 (3),
498–510.Acknowledgments
Our special appreciation goes to Prof. R. Uzuoka of the
University of Tokushima, Japan, for his valuable suggestions
in developing the cyclic elasto-plastic model. In addition, the
authors sincerely acknowledge the Ministry of Education,
Culture, Sports, Science, and Technology (Grants-in-Aid for
Scientific Research #21360221), Japan for providing the
financial assistance for the current research.(5) The simulations of stress paths and stress–strain relation-
ships of Toyoura sand subjected to undrained cyclic
torsional shear loadings were significantly improved after
employing the proposed four-phased stress–dilatancy
relationship.
(6) Complete liquefaction and steady state can be accurately
simulated by the proposed model; however, further modifica-
tions are required to address the strain-softening behavior at
large strain levels of dense sand.
In order to develop a generalized three-dimensional model based
on the proposed stress–strain description, extended research has
been initiated as preliminarily described by Namikawa et al. (2011)
on the results from a numerical simulation of drained cyclic
loading behavior. To maintain the objectivity, in terms of the
independence from the coordinate systems, the three-dimensional
model adopts the concept of an infinite number of nesting surfaces
(Mroz et al., 1978). Its application to the simulation of undrained
cyclic loading behavior, as is the case with the current study, is
under way. Details will be reported elsewhere.
undations 55 (2015) 559–574 573elationships of sand in the simulation of volumetric behavior during cyclic
orsional shear loadings. Soils Found. 54 (4), 845–858.
HongNam, N., Koseki, J., 2005. Quasi-elastic deformation properties of
Toyoura sand in cyclic triaxial and torsional loadings. Soils Found. 45 (5),
19–38.
HongNam, N., Koseki, J., 2008. Deformation characteristics of dry Toyoura
sand in large cyclic torsional loading and their modeling. In: Burns, Mayne,
Santamarina (Eds.), Proceedings of the Fourth International Symposium on
Deformation Characteristics of Geomaterials, vol. 2, IS Atlanta, USA,
pp. 801–807.
Iai, S., Matsunaga, Y., Kameoka, T., 1992. Strain space plasticity model for
cyclic mobility. Soils Found. 32 (2), 1–15.
Ishihara, K., Li, S., 1972. Liquefaction of saturated sand in triaxial torsion
shear test. Soils Found. 12 (2), 19–39.
Ishihara, K., Tatsuoka, F., Yasuda, S., 1975. Undrained deformation and
liquefaction of sand under cyclic stresses. Soils Found. 15 (1), 29–44.
Jefferie, M.G., 1993. Nor-Sand: a simple critical state model for sand.
Géotechnique 43 (1), 91–103.
Kiyota, T., Sato, T., Koseki, J., Abadimarand, M., 2008. Behavior of liquefied
sand under extremely large strain levels in cyclic torsional shear tests. Soils
Found. 48 (5), 727–739.
Koseki, J., Hamaya, S., Tatsuoka, F., Maeshiro, N., 1998. Elastoplastic
deformation characteristics of Toyoura sand during liquefaction. Geotech-
nical Engineering and Soil Dynamics III, vol. 75. Geotechnical Special
Namikawa, T., Koseki, J., De Silva, L.I.N., 2011. Three dimensional modeling
of stress-strain relationship of sand subject to large cyclic loading. In:
Proceedings of the 5th International Conference on Earthquake Geotechni-
cal Engineering. Santiago, Chile, paper-ID: TMSNA, CD-ROM, pp. 1–10.
Nishimura, S., 2002. Development of Three Dimensional Stress–Strain Model
of Sand Undergoing Cyclic Undrained Loading and Stress–Axes Rotation
(M.Eng. thesis). Dept. of Civil Engineering, The University of Tokyo,
Japan.
Nishimura, S., Towhata, I., 2004. A three-dimensional stress–strain model of
sand undergoing cyclic rotation of principal stress axes. Soils Found. 44
(2), 103–116.
Oka, F., Yashima, A., Kato, M., Sekiguchi, K., 1992. A constitutive model for
sand based on the non-linear kinematic hardening rule and its application.
In: Proceedings of 10th World Conference on Earthquake Engineering,
pp. 2529–2534.
Oka, F., Yashima, A., Tateishi, Y., Taguchi, Y., Yamashita, S., 1999. A cyclic
elasto-plastic constitutive model for sand considering a plastic-strain
dependence of the shear modulus. Geotechnique 49 (5), 661–680.
Seed, H.B., Lee, K.L., 1966. Liquefaction of saturated sands during cyclic
loading. Proc. ASCE SM6, 105–134.
Tatsuoka, F., Shibuya, S., 1991. Modelling of non-linear stress–train relations
of soils and rocks – Part 2: new equation, Seisan-kenkyu. J. IIS Univ.
Tokyo 43 (10), 435–437.
Tatsuoka, F., Jardine, R.J., Lo Presti, D., Di Benedetto, H., Kodaka, T., 1997.
Characterising the pre-failure deformation properties of geomaterials,
Theme Lecture for the Plenary Session No. 1, Proceedings of XIV IC on
stress–strain relations of sand in cyclic plane strain loading. J. Geotech.
L.I.N. De Silva et al. / Soils and Foundations 55 (2015) 559–574574generalized plasticity. J. Eng. Mech. 132 (12), 1380–1391.
Liou, C.P., Richart, F.E., Streeter, V.L., 1977. Numerical model for liquefac-
tion. J. Geotech. Eng. Div. ASCE, 103; 589–606.
Liyanapathirana, D.S., Poulos, H.G., 2002. A numerical model for dynamic
soil liquefaction analysis. Soil Dyn. Earthq. Eng. 22 (9–12), 1007–1015.
Masuda, T., Tatsuoka, F., Yamada, S., Sato, T., 1999. Stress-strain behavior of
sand in plane strain compression, extension and cyclic loading tests. Soils
Found. 39 (5), 31–45.
Mroz, Z., Norris, V.A., Zienkiewicz, O.C., 1978. An anisotropic hardening
model for soils and its application to cyclic loading. Int. J. Numer. Anal.
Methods Geomech. 2, 203–221.Geoenviron. Eng. ASCE, 129; 450–467.
Towhata, I., Ishihara, K., 1985a. Shear work and pore water pressure in
undrained shear. Soils Found. 25 (3), 73–85.
Towhata, I., Ishihara, K., 1985b. Undrained strength of sand undergoing cyclic
rotation of principal stress axes. Soils Found. 25 (2), 135–147.
Vaid, Y., Thomas, J., 1995. Liquefaction and postliquefaction behavior of
sand. J. Geotech. Eng. 121 (2), 163–173.Ling, H.I., Yang, S., 2006. Unified sand model based on the critical state and SMFE, vol. 4, Hamburg, pp. 2129–2164.
Tatsuoka, F., Masuda, T., Siddiquee, M.S.A., Koseki, J., 2003. Modeling thePublication, ASCE385–397.
Koseki, J., Yoshida, T., Sato, T., 2005. Liquefaction properties of Toyoura
sand in cyclic torsional shear tests under low confining stress. Soils Found.
45 (5), 103–113.