- Email: [email protected]
Coupled theory of mixtures for clayey soils
PII:
Computers and Geotechnics, Vol. 20, No. 3/4, pp. 195--222. 1997 Q 1997 Elsevier Science Ltd. All rights reserved Printed in Great Britain SO266-352X(97)00003-7 0266-352X/97 S 15.00 + 0.00
ELSEVIER
Coupled Theory of Mixtures for Clayey Soils G. Z. Voyiadjis Department
of Civil and Environmental
& M. Y. Abu-Farsakh
Engineering, Louisiana LA 70803, U.S.A.
State University,
Baton
Rouge.
ABSTRACT In this work, elasto-plastic coupled equations arc formulated in order to describe the time-dependent deformation of saturated cohesive soils (twophase state). Formulation of these equations is based on the principle of virtual work and the theory of mixtures for inelastic porous media. The theory of mixtures for a linear elastic porous skeleton was$rst developed by Biot (Theory of elasticity and consolidation for a porous anisotropic solid, Journal of Applied Physics, 1955, 26, 188-185). An extension of Biot’s theory into a nonlinear inelastic media was performed by Prevost (Mechanics of continuous porous media, International Journal of Engineering Science, 1980, 18, 787-800). The saturated soil is considered as a mixture of two deformable media, the solid grains and the water. Each medium is regarded as a continuum andfollows its own motion. Thejow, of pore-water through the voids is assumed to follow Darcy ‘s law’. The coupled equations are developed for large deformations with finite strains in an updated Lagrangian reference frame. The coupled behavior of the tlz’ophase materials (soil-water state) is implemented in a jnite element program. A modtfied Cam-clay model is adopted and implemented in the$nite element program in order to describe the plastic behavior of clayey9 soils. Penetration of a piezocone penetrometer in soil is numerically simulated and implemented into a finite element program. The piezocone penetrometer is assumed to be infinitely st@ The continuous penetration of the cone is simulated by applying an incremental vertical movement of the cone tip boundary. Results of the finite element numerical simulation are compared with experimental measurements conducted at Louisiana State University using the calibration chamber. The numerical simulation is carried out for two cases. In the jrst case, the interface friction between the soil and the piezocone penetrometer is neglected. In the second case, interface friction is assumed between the soil and the piezocone. The results of the numerical simulations are compared with experimental laboratory measurements. 0 1997 Elsevier Science Ltd.
195
196
G. Z. Voyiadjis and M. Y. Abu-Farsakh
INTRODUCTION The soil skeleton consists of solid grain particles and voids that are filled with water and/or air. Soil can be considered as a mixture of multiphase deformable medium of solid grains, water, and air, each of which is regarded as a continuum and follows its own motion. The general theoretical framework that forms the basis of the theory of mixtures was first developed by Truesdell and Toupin [l]. Further development of the modern theory of mixtures was presented by Green and Naghdi [2] and Eringen and Ingram [3]. The theory of mixtures for a small strain linear elastic porous media was first developed by Biot [4] assuming that fluid flow through solid phase obeys Darcy’s law. An extension of Biot’s theory for large strain nonlinear inelastic porous media was presented by Prevost [5]. The theory of mixtures was developed in a Eulerian reference frame. The formulation by Prevost [5%3] was also developed in a Eulerian reference frame. It was shown that the use of Jaumann stress rate in a Eulerian formulation may cause inaccuracies in large strain problems [9,10]. Kiousis and Voyiadjis [ 1 I] developed a theory for two-phase material based on the concept of the theory of mixtures in a Lagrangian reference frame using the material stress rate. In this paper, an elasto-plastic analytical model is formulated in order to describe the nonlinear elasto-plastic, time-dependent deformation of a saturated cohesive soil skeleton (two-phase state). Formulation of the coupled equations for the model is based on the principle of virtual work and the theory of mixtures for inelastic porous media. The theory for large strain nonlinear inelastic porous media used here is based on the work of Prevost [5] and Kiousis and Voyiadjis [ 111. However, the coupled equations are formulated for large deformations using finite strains in an updated Lagrangian framework. Recently, there has been an increased concern toward the evaluation of different engineering soil parameters using in situ testing. Piezocone penetrometer tests (PCPT) are one of the most widely used for in situ testing. Many interpretation procedures based on bearing capacity model [ 121, cavity expansion theories [13-l 51, strain path method [ 166181, semi-empirical methods [19], dislocation methods [20], and the finite element techniques [21], have been proposed by several investigators to analyze the piezocone penetration and to evaluate some engineering soil parameters. Most of the methods that are available assume certain deformation patterns and small strains that do not match with the experimental results. Strains are reported to exceed 100% [I 6,17,21]. To accurately analyze the PCPT, it is necessary to have a numerical model that takes into consideration the large deformations, nonlinear soil behavior, and changing of the boundary condition with penetration. The piezocone penetration test results depend on many factors,
Theory
of mixtures
for clayey soils
197
including soil type and compressibility characteristics, stress history, and the presence of fissures and cracks. To validate the numerical model, it is necessary to compare it with controlled laboratory calibration tests of known soil types, characteristics, and stress history. The proposed analytical model is used for the analysis of the piezocone penetration problem. A modified Cam-clay model is adopted to describe the plastic behavior of the clayey soil. The penetration of the piezocone penetrometer in cohesive soils is numerically simulated and implemented in the finite element program. The continuous penetration of the piezocone is simulated by applying an incremental vertical movement of the cone tip boundaries. The problem is treated as a boundary problem that includes changing of the boundary condition as penetration proceeds. The numerical simulation is performed for two cases. In the first case, no interface friction is assumed between the soil and the piezocone penetrometer while, in the second case, interface friction is assumed between the soil and the piezocone penetrometer. A coefficient of interface friction equal to 0.25 is used in this work, which corresponds to the angle of friction of 6= 14” between the soil and the piezocone. The results of the numerical simulation are compared with experimental laboratory measurements of the miniature PCPT performed on a clay calibration chamber conducted at Louisiana State University (LSU) [22,23].
COUPLED
FIELD
EQUATIONS
The coupled field equations proposed herein are formulated in an updated Lagrangian framework. The principle of virtual work is expressed here using the updated Lagrangian formulation given by Bathe [24] as follows:
s
;+‘&sS( ;+‘E,@‘%’ = “+‘R
(1)
“V
Ef16is the total strain and is decomposed into the linear strain, #As, and the nonlinear strain nqAs: n+l
,,
CAB =
&AB
+
nT]AB
nUA,B
+
db3.A)
+
nuK,B)
(2)
where neAB
=
;1(
1 nl]AB= 2 (nUK,A
(4)
G. Z. Voyiadjis and M. Y. Abu-Farsakh
198
where u is the displacement vector. In Eqn (1) :+~s.Q is the second PiolaKirchhoff stress tensor at n + 1 referred to n configuration as indicated in Fig. 1. It is expressed in terms of the Cauchy stress tensor “0~s such that
:+‘SAB= naAB+
(5)
naiad
where
naAB= In Eqn (5) “ASAB is the y1+ 1 as indicated in Fig. pore water pressure and external virtual work due as follows:
azB+ Pd.4~
(6)
incremental
stress tensor from configuration n to Cauchy stress tensor P, is the GABis the Kronecker delta. In Eqn (l), "+'Ris the to the applied loads and tractions and is evaluated 1. aiB is the effective
The relation between the coordinates solid medium is given by
at n + 1 and n configurations
I
for the
Yo ,Y. ,Yn+l
/0
Vo+l
0 V.
n+l configuration
n configuration
/ 0 configuration vo D
/
zo,Z”
Fig. I. Configuration
,Z”hl
of a body in motion
x0
,xn,xn+1
with a fixed Cartesian
coordinate
system
Theory of mixtures for clayey soils
199
Zi = n+lq = “Qy+ S&
(8)
The relation between the Cauchy stress tensor and the second Piola-Kirchhoff stress tensor is given by
where
and Js is the corresponding Elasto-plastic
Jacobian for solids.
constitutive equation
For finite deformations, assumed to be
the constitutive o’s (3
ab
=
equation for the solid skeleton is
%cdds,d
where Ds is the elasto-plastic stiffness matrix. The effective co-rotational stress rate tensor 0’s is given by [25] as
where a’” is the effective Cauchy stress rate tensor, 0’” is the effective Cauchy stress tensor, and
is the modified spin tensor. In Eqn (12a)
is the total spin tensor, and J+$Lis the plastic spin tensor. When I%$ = 0, the co-rotational stress rate tensor reduces to the Jaumann stress rate tensor. The deformation gradient ds is given as follows:
200
G. Z. Voyiadjis and M. Y. Abu-Farsakh (13)
The relation co-rotational
between the total co-rotational stress rate tensor is given by
where P, is the material
derivative
stress
rate and the effective
of P, given as V”.gradP,
By differentiating both sides of Eqn (9) with respect gian reference frame, one obtains
(15) to time in the Lagran-
Making use of expression (13) for the spatial deformation gradient (16) together with Eqns (6) and (1 l)-(14) one obtains the following sion:
Substituting
in Eqn expres-
Eqn (10) into Eqn (17) and using the relation
(18) one reduces
Eqn (17) into the following
If one neglects the plastic spin tensor relation: ,AB
=
DlZ,,cDic~
relation:
@,
+
Eqn (19) reduces to the following
px”,,,~B,b~w~ab
(20)
Theory of mixtures for clayey soils
201
where 1 is the strain rate tensor. In Eqn (19) D* is given as follows:
The left-hand side of Eqn (1) may be expressed as follows: ( naAB +
&AB)J(
neAB f
%AB)dnV
=
“V
“OAB~( @AB
+
“VAB)~~~
“V (22)
&SAB6(
+
neAB +
nl]AB)dnV
1 “V
However, Ans~a may be expressed as follows: [+A/ AnsAB
t+At iABdt
=
=
D&D
ic-dt
s
J
t
t
t+At +
A ~SAB
Substituting
=D&D
(23) &dt
‘?&%,a
AECD
+ f~A,a%,#%i
Eqn (23) into Eqn (22), one obtains the following expression:
DiBcD(A &CD
+
A nWD)&eABdp
-I
“V
DiBCDA
&‘CD&VABdnV
“V
+
J
(“&j
+
“V
=
%~AB)klABdnV+
s
JSXSA,,XB,bsabAPw(6neAB
“V
( nc7;B + nPw8A&,eABdnV
"+'RI
“V
Equation (24) is the first derived coupled equation.
+ &‘lAB)dnV
202
G. Z. Voyiadjis and M. Y. Abu-Farsakh
Theory of mixtures The saturated soil is able media, the solid uum and follows its (assuming chemically
considered here as a mixture consisting of two deformgrains and water, each of which is regarded as a continown motion. The balance of mass for each constituent inert) can be written as [5] (a = s or w)
ba + p”divv” = 0
(25)
where p” is the macroscopic mass density of the a! constituent. If one assumes both solid grains and water to be incompressible materials, then the balance of mass can be written in terms of the soil porosity n’+‘:,solid velocity vS and water velocity v”‘ as follows [5]:
div(v’) - div(v”‘) = $ [div(u’) + (v”’ - v”)grad(n”‘)] The general form of Darcy’s ium is given by
law for the flow of water through
(vl - 9) = - -&K”“(grad(P,)
- p,b)
YW
porous
med-
(27)
where K”‘” is the permeability tensor, b is the body force vector, and yW is the unit weight of water. Taking the divergence of both sides of Eqn (27) and substituting it in Eqn (26) together with Eqn (18) and making use of the following relations:
apw
x”-
a -
=
a.22
one then obtains the second formulation as follows:
apw
(284
Csbaxc
azb -=
XD,, g
coupled
D
equation
in an updated
Lagrangian
Theory of mixtures for clayey soils
NUMERICAL
203
FORMULATION
The finite element discretization is used here for the displacement pore water pressure P, as follows: u=M.U
u and the
(30a)
P,=N.W
(30b)
where M is the displacement shape function, N is the pore water pressure shape functions, U is the nodal displacement, and W is the nodal pore water pressure. The incremental strains can be expressed as follows: 6e=
BL.GU
(3la)
Jr] = BNL. 6U
and the pore water pressure
gradient
ap,
-=
axB
are given such that N,BW
where BL and BNL are the linear and nonlinear strain-displacement Substituting Eqn (31) into Eqn (24) one obtains ,,KAU +” S-JAW= ,,a where .K =
is the elasto-plastic
.KL
+
stiffness matrix,
&a_
(3lb)
+
nKzL + “KS
(3lc)
matrices.
(32)
204
G. Z. Voviadjis and M.
nKL
s
=
Y. Abu-Farsakh
&D*
,,BLd”V
&I*
,,BNLdnV
”v is the linear stiffness matrix,
nKNL
=
s
”v is the nonlinear
stiffness matrix.
Jc =
.I
,,B~,%,,B~,d”V
“V
is the geometric
is
the coupling
stiffness matrix,
matrix,
and
The condition that the continuity Eqn (29) applies throughout nuum using Galerkin’s weighted residual method requires that
the conti-
(33) The weak form of Eqn (33) is obtained this equation:
s
J’C$ i,P,d”V
I
”v
by applying
g J”C;; C&l l$d”B(s+p,Ba)
Green’s
theory
[26] to
$d”V
“V (34) -
c J “S
q,,P,d
“A = 0
Theory of mixtures for clayey soils
205
where &,=NG
(35a)
is the weighted residual (virtual pore pressure)
aL
z=
w
N
(35b)
,A
qn
is the seepage velocity normal to the boundary. Substituting Eqn (3 1) and Eqn (35) into equation Eqn (34), one obtains -,,fl=AU
+ ,,\k&AW = J-I
(36)
where
is the flow matrix, and “II = -6tG
- Gt\kW” +
J
NTq,dnA
“S
From Eqn (32) and Eqn (36), the complete coupled behavior of the twophase skeleton-fluid state can be expressed as follows: (37)
where AU is incremental nodal displacement, nodal excess porewater pressure.
NUMERICAL
and A W is the incremental
EXAMPLES
A number of examples are analyzed to demonstrate the accuracy of the analytical model described. In the analysis, the eight node quadratic isoparametric finite element (Q8P4) is used.
G. Z. Voyiadjis and M. Y. Abu-Farsakh
206
One-dimensional
elastic consolidation
The one-dimensional elastic consolidation problem is studied here. The finite element mesh used in the analysis of this problem and the corresponding boundary conditions are shown in Fig. 2. It consists of 20 Q8P4 finite elements with an H/B ratio of 5. The traction q is applied instantaneously to the top surface at time t = 0. Thereafter, the traction is held constant with drainage allowed only at the top surface. A closed form analytical solution for the small deformation case as a function of time can be found in Lambe and Whitman [27]. The average degree of consolidation as a function of the dimensionless time factor (T = Cvt/H2) is obtained from the finite element analysis for a load ratio of q/E = 0.1. This is compared with the analytical solution as shown in Fig. 3; q is the load intensity, E is the Young’s modulus, and C, is the coefficient of consolidation. The finite element results show good agreement with the analytical solution. Figure 4 represents the settlement of the top surface 6 as a function of the dimensionless time factor (T), for various q/E ratios (0.1 and 0.5). Figure 4 shows that the settlements obtained from the updated Lagrangian finite strain deformation analysis are lower than those obtained from small deformation analysis. It is noted that, as the q/E load ratio increases, the
Drainage
Fig. 2. One-dimensional
elastic consolidation:
finite element
mesh
Theory of mixtures for clayey soils
207
difference between the small and the finite deformation analysis also increases. This is expected since the large deformation solution is stiffer than the small deformation case.
0
,
0.2 -
e
.3 s ‘i; 9
0.4 -
a 8 % :
0.6 -
rQ B 0.8 -
Analytical 0
Finite Elements
0.1
0.01
0.001
Time
Fig. 3. One-dimensional
1
Factor (T)
elastic consolidation:
i
\
degree of consolidation.
\
WE - 0.6 \ ,
0.3 -
\ -
__
Small Deformation Updated
Fig. 4. One-dimensional
Lagrangian
_-__-_-_
Lqe
Deformation
elastic consolidation:
settlement.
G. Z. Voyiadjis and M. Y. Abu-Farsakh
208
Two-dimensional
elastic consolidation
The finite element mesh and the corresponding boundary conditions for the problem are shown in Fig. 5. A uniform distributed strip load of intensity q and width 2B is applied instantaneously to the top surface at time t = 0. Thereafter, the load is held constant with drainage allowed only at the top surface. Figure 6 shows the ratio of the computed excess pore pressure at 0.5B beneath the center of the footing to the initial excess pore pressure (t = 0) as a function of the dimensionless time factor (0. The finite element analysis results are compared with the closed form solution of the excess pore pressure for an infinite elastic half-space [28]. During the early times of consolidation, there is a slight increase in the pore pressure as a result of the particular distribution of the pore pressure. This phenomenon is known as the Mandel-Cryer effect. A good agreement is obtained between the numerical solution and the analytical solution, as shown in Fig. 6. The computed settlement of the point located at the center of the strip load 6 as a function of the dimensionless time factor T for various q/E ratios (0.1 and 0.3) are shown in Fig. 7. Similar behavior to the one-dimensional elastic consolidation can be seen here for the two-dimensional elastic consolidation. It is noted that, as the q/E load ratio increases, the difference between small and finite deformation analysis also increases.
+tti
Fig. 5. Two-dimensional
Drainage
elastic consolidation:
finite element
mesh.
Theory of mixtures for clayey soils
0.4
-
209
Analytical
0.01
0.1
Time Factor Q
Fig. 6. Two-dimensional
elastic consolidation:
pore water pressure at depth z =0.5B.
0.4
0.8
- - --
\
Small Deformation Updated Lagrangian Large Deformation
1.2
I
’ ’ “““I
0.001
0.01
’ “““‘I 0.1
’ “‘m’l
Time Factor
Fig. 7. Two-dimensional
elastic consolidation:
1
\
\
\
’ ’ “““~ 10
\
\
\
\
’ ’ “In 100
(T)
settlement of the mid-point.
210
G. Z. Voyiadjis and M. Y. Abu-Farsakh
Piezocone penetrometer analysis In this work, the penetration of a piezocone penetrometer in cohesive soils is numerically simulated assuming that the piezocone penetrometer is infinitely stiff and no tensile stresses are allowed. The projected piezocone penetrometer area is 1 cm2 with a cone apex angle of 60” which is the same size as the cone used in the experiments by Kurup et al. [23] and Voyiadjis et al. [22]. Simulation is performed for two cases. The first case assumes that the soilpenetrometer interface friction is negligible. In case 2, it is assumed that the coefficient of the soil-penetrometer interface friction is 0.25, which corresponds to an angle of friction 6 =14” between the soil and the piezocone surface. The continuous penetration of the piezocone is simulated by applying an incremental vertical movement of the cone tip boundary, as shown in Fig. 8. For case one, the nodes along the stress boundary of the cone tip and along the piezocone shaft are allowed to slide from the beginning of the incremental penetration, while, for case two, the nodes along the inclined tip surface and along the piezocone shaft are first prevented from sliding along the surface until sliding potential occurs. The sliding potential is reached along the surface when friction forces of these nodes exceed the allowable friction forces. After that, sliding along the stress tip surface and along the shaft surface is allowed. During the incremental penetration for cases one and two (after sliding potential is reached), the nodes are checked for possible changes in
Position (n)
Fig. 8. Incremental
Position (n+l)
penetration
of the piezocone
penetrometer.
Theory of mixtures for c1aye.y soils
211
the boundary conditions. When tensile stresses occur in the node immediately below the lower end of the cone at location (l), the node is released as indicated by location (2) and allowed to move freely. Once the free node at location (2) reaches the cone tip boundary, its movement for case one will be restricted along the skew surface of the cone tip as indicated by (3). In case two, the nodes are first fixed and vertically incrementally displaced until the sliding potential is reached. Then they are allowed to slide along the surface. When the horizontal co-ordinate of the nodes at location (3) exceeds the penetrometer radius, these nodes are released at location (4) until they return to the boundary of the shaft. They are then restricted from further movement in the horizontal direction, as shown in location (5). The numerical simulation is compared with the PCPT conducted on cohesive soil specimen prepared in a calibration chamber at LSU. The soil specimen used is a K-33 prepared by consolidating a soil slurry of 33% kaolinite and 67% Edgar fine sand (D6,,/Di0= 1.4) at a water content of twice the liquid limit. The soil specimen has a liquid limit, WL = 20% and a plastic limit, W, = 14%. The specimen is consolidated to an isotropic effective stress of 207 kPa against a back pressure of 138 kPa. Full details about specimen preparation, test procedures and type of penetrometer used are found in Kurup et al. [23] and Voyiadjis et al. [22] The filter used to measure the porewater pressure is located either in the lowest quarter of the cone at the very tip (Ul configuration) or starting 0.5 mm above the base of the cone with 2 mm vertical height (U2 configuration). The modified Camclay model is adapted in this work to describe the plastic behavior of soil. The yield locus for modified Cam-clay is given by f(P,
q,
P&))
z M2P2 -
M2PcP+ cl*= 0
(38)
It represents an ellipse in the p - q plane. In Eqn (38) P and q are the mean effective and deviatoric stresses. PC is the strain hardening parameter representing the apex of the yield locus ellipse in the P axis. The soil parameters used for the modified Cam-clay model are the slope of the virgin consolidation line in the e-1nP space, h=O.ll; the slope of the unloading-reloading line in the e-1nP space, K =0.024; the slope of the critical state line, M =1.2; r = 1.162; the coefficient of permeability for the soil is 0.5 x 1O-9 m/s; and the Poisson’s ratio u =0.3. A number of trial finite element meshes are studied first with different degrees of refinement in order to obtain the appropriate mesh for the analysis of the piezocone penetration problem. Figure 9 shows the computer dimensionless quantities qc/qc5 and Ul /Ul 5 as a function of the number of finite elements used along the inclined surface (as a measure of degree of
212
G. Z. Voyiadjis and M. Y. Abu-Farsakh
refinement) for a penetration depth of 1Omm. qC/qC5 is the ratio of the cone tip resistance to the cone tip resistance computed using five elements along the inclined surface. Ul /Ul s is the ratio of the excess pore pressure at Ul configuration to the one computed using five elements along the inclined surface. From Fig. 9, it is clear that, with the increase of refinement, the numerical results converge to a unique solution. Figure 10 shows the finite element mesh used in the analysis of the penetration of the piezocone penetrometer in cohesive soils. The calculated cone tip resistance profiles during piezocone penetration and those obtained experimentally are compared in Fig. 11. In Fig. 12, the developed excess pore pressure profiles during piezocone penetration obtained numerically are compared with those obtained from experimental testing. It can be seen that the cone resistance obtained from the numerical simulation increases rapidly until it reaches a steady-state condition at a depth of 25-30mm, as shown in Fig. 11. The same behavior is also obtained for the development of excess pore pressure. The computed cone tip resistance and excess pore pressures develop at faster rates than those obtained experimentally. This difference may be due to the errors encountered at the initial experimental stage of testing and also to the lack of confining stresses at the top surface of the soil specimen. This can be seen if one compares the
2 Number Fig. 9. Convergence
3 of Elements
4 5 along the Inclined
6 Surface
of the cone time resistance and pore water pressure with mesh refinement.
of Ul configuration
Theory of mixtures for clayey soils
213
two cone tip resistance test profiles. The only difference between them is the filter location. The final values obtained from the computational model after the steady-state condition is reached are quite close to those obtained experimentally. In case 1, in which no friction is assumed between the soil and the piezocone surface, the cone tip resistance value at the steady state is
Fig. 10. Finite element mesh of the piezocone
P
r C = :
-60.00
penetration
problem.
--._--_._ 1.---
Finka Elemmt IFrktionl
I
i 0.00
0.40
0.84
Cone Resistance
Fig. 11. Cone resistance
1.20
1 .BO
(MPa)
profiles during piezocone
penetration.
214
G. Z. Voyiadjis and M. Y. Abu-Farsakh
0.60 0.40 0.20 Excess Pore Pressure (MPa)
0.80
(a)
1
.200.00 1
-
Fin::
I 0.00
0.60 0.40 0.20 Excess Pore Pressure (MPa)
I
1
0.80
(b)
Fig. 12. Excess pore pressure
profiles during piezocone penetration; configuration U2.
(a) configuration
Ul, (b)
1.131 MPa (as shown in Fig. 11). In case 2, where the soil-piezocone coefficient of interface friction p is taken as 0.25, the cone tip resistance at the steady state is 1.294 MPa. The cone resistance obtained experimentally is 1.255 MPa. The computed excess pore water pressures at the steady state for configuration Ul (at the lowest quarter of the cone tip) are 0.565 and 0.774 MPa for cases 1 and 2, respectively, as shown in Fig. 12(a). However, the
Theory of mixtures for clayey soils
215
measured excess pore water pressure at the steady state is 0.608 MPa. The computed excess pore water pressures for the steady state in configuration U2 (above the cone base) are 0.532 and 0.542 MPa for cases 1 and 2, respectively, as shown in Fig. 12(b). The corresponding experimental measured excess pore pressure at the steady state is 0.593 MPa. Figure 13 shows the dissipation of excess pore water pressure with time. Figure 13(a) compares the dissipation of excess pore pressure for configuration
h
2 i+
P
0.6 -
t!
2 0 8
ti
0.4 -
0 t i
2c
0.2 -
-
---.0.0
,
, ,,,,,,,
,.O
,
,
,,,,,I,
,o.o
/
100.0
I
,,,111/
I
I
In-n
10000.
1000.0
Time (set) (a)
_ -
1.0
~
PcPTul
---
Finit. Element (NO Friction)
-.
-
(Kurup
dill.
lSs4,
Finita Ekment(FMtion) . .
10.0
100.0
Time (set) (b)
Fig. 13. Dissipation
of excess pore pressures;
(a) configuration
U1, (b) configuration
U2.
G. Z. Voyiadjis and M. Y. Abu-Farsakh
216
Ul obtained both numerically and experimentally. Figure 13(b) compares the dissipation of the excess pore water pressures for configuration U2. It is obvious that the finite element analysis results are close to those obtained experimentally with the exception that there is no sudden drop of the pore pressure when the penetration stops. The distribution of the excess pore water pressures developed around the piezocone penetrometer for a penetration depth of 25mm are shown in Fig. 14. For the case with no interface friction (Fig. 14(a)), the maximum concentration excess pore water pressure is obtained in the upper third of the cone, while for case 2 with the coefficient of interface friction p of 0.25, the maximum concentration of excess pore water pressure is around the lower of the cone as shown in Fig. 14(b). In both cases, negative pore water pressure is developed at about 2-3 diameters above the cone base. Figures 15-18 present the contours of the effective stresses around the piezocone penetrometer for a penetration of 25 mm. The initial stress state is isotropic with ai = 0: = ok = 207KPa, and trz = 0. The radial effective stress ai around the piezocone penetrometer for case 1, for which friction is neglected, increases in the region above the mid-point of the cone up to 424.1 KPa immediately above the cone base. It then drops to 48.3 KPa below the cone tip, as shown in Fig. 15(a). For case 2, in which interface friction is -
Normalized Radial Distance 14. Contours
1
rho
Normalized
of excess pore pressures around the piezocone of 25 mm; (a) no friction, (b) coefficient of interface
Radial Distance
rho
penetrometer for penetration friction p = 0.25.
Theory of mixtures for clayey soils
20.0
217
J
L P
20.0
24.0
0.00
¶.oo
4.00
0.00
1
24.0
(
I
0.0
Normalized Radial Distance rho
2.0
I
60
0.0
1
c.10
I
Normalized Radial Distance rho
Fig. 15. Contours of radial stresses 0: around the piezocone penetrometer for penetration 25 mm; (a) no friction, (b) coefficient of interface friction p = 0.25.
Normalized Radial Distance
rho
Normalized Radial Distance
of
rho
Fig. 16. Contours of axial stresses ai around the piezocone penetrometer for penetration 25 mm; (a) no friction, (b) coefficient of interface friction w = 0.25.
of
G. Z. Voyiadjis and M. Y. Abu-Farsakh
218
considered, the stress increases up to 405.7 KPa and drops down to 54.3 KPa, as shown in Fig. 15(b). The radial stress bulb region of the zone influenced by penetration extends to approximately four times the piezocone diameter. The axial effective stress a: arising from penetration is presented in Fig. 16. The stress increases in the region below the cone face and decreases in the region above it. The maximum stress at the end of a 25mm penetration is located below the cone tip, while the minimum stress is located above the base of the piezocone. The stress ranges from 305.5 to 62.1 KPa and from 347.3 to 45.6 KPa for cases 1 and 2, respectively. Figure 17 shows the distribution of the tangential stresses aa around the cone. The stresses decrease from 207 KPa (initial stress) to 43.4 and 47.6 KPa, respectively, for cases 1 and 2. The minimum tangential stress (maximum decrease) is located at the lower half of the cone. In Fig. 18, the developed shear stresses t ,.Z are presented. The maximum shear stresses induced are -166 KPa for case 1 and -201.6 KPa for case 2. The shape of the shear stress distribution around the piezocone is in the shape of a bulb with the maximum value located at the cone face and oriented at an angle of approximately 45” to the cone face. The shear stress bulb is also extended
0.00
2.00
Normalized
4.00
0.00
Radial Distance
0.00
rho
1-r 2.0
0.0
Normalized
Fig. 17. Contours of tangential stresses CJ~around the piezocone of 25 mm; (a) no friction, (b) coefficient of interface
Radial Distance
rho
penetrometer for penetration friction p= 0.25.
Theory of mixtures for clayey soils
c =
219
32.
p” P z
30.
6 z 23.
24.0 I 0.0
Normalized Radial Distance rho
, 2.0
I 4.0
I 6.0
Normalized Radial Distance rho
Fig. 18. Contours of shear stresses r,, around the piezocone penetrometer for penetration 25 mm; (a) no friction, (b) coefficient of interface friction p = 0.25.
of
behind the cone face. This extension is bigger for case 2 than it is for case 1. For case 2, there is a shear stress concentration around two points. One is at the lower face of the cone, and the other is at the upper portion of the cone face near the cone base.
CONCLUSIONS Coupled elasto-plastic equations are formulated in order to describe the time-dependent deformation of saturated cohesive soils. The formulation of these equations is based on the principle of virtual work and the theory of mixtures for inelastic porous media. The coupled equations are implemented in a finite element program. The modified Camclay constitutive model is used to describe the plastic behavior of the clayey soils. Penetration of the piezocone penetrometer is numerically simulated and implemented into the finite element program. Simulation is carried out for two cases. In the first case, the interface friction between the soil and the piezocone penetrometer is neglected. In the second case, the soil-piezocone penetrometer interface
220
G. Z. Voyiadjis and M. Y. Abu-Farsakh
friction is considered. Results obtained from the numerical simulations are compared with the experimental laboratory measurements conducted at the LSU calibration chamber. Good correlation is obtained between the experimental results and those simulated numerically using the proposed model. The following conclusions can be made from the analysis of the piezocone penetration. The steady-state condition is reached in the numerical analysis at a depth of 25 to 30mm, which is at a faster rate than those obtained experimentally. This may be attributed to some errors encountered at the initial experimental stage of testing and also due to the lack of confining stresses near the top surface of the soil specimen. The pore pressure distribution around the cone tip is concentrated around the upper third of the cone if friction is neglected and around the lower third of the cone if interface friction is taken into consideration. In addition, a negative pore pressure may be generated above the cone base due to soilpiezocone separation in this area, especially at the early stages of testing. This raises the question as to where is the best location to place the transducer in order to avoid possible misleading calculations. The authors believe that the middle third is the best location to place the transducer because it is the most representative of the average pore pressure which is more uniform and less sensitive. The bulb of stress changes surrounding the cone tip is extended radially outward to about four times the piezocone diameter. The stresses change faster in the vicinity of the cone tip. This makes the problem very localized around the cone tip, which may lead to incorrect calculation due to the high stress redistribution in this area. The effect of considering the interface friction between the soil and the piezocone penetrometer is paramount for the shear stress distribution around the cone tip and the generated excess pore pressures close to the cone surface, which has an affect on the calculated cone tip resistance. The effect is minimal for other stress distributions.
ACKNOWLEDGEMENTS This research was supported by grant MSS-9312707. The authors encouragement of Dr Priscilla P. Directorate for Engineering, NSF. portation Research Center under edged.
the National Science Foundation under wish to acknowledge the support and Nelson, Program Director of G3S/CMS, Partial support of the Louisiana TransTask Order No. 736-99-0293 is acknowl-
Theory of mixtures for clayey soils
221
REFERENCES 1. Truesdell, C. and Toupin, R., The classical field theories. In Hun&& der Physik, Vol. 3, No. 1, ed. S. Flugge. Springer-Verlag, Berlin, 1960. 2. Green, A. E. and Naghdi, P. M., A dynamical theory of interacting continua. International Journal of Engineering Science, 1965, 3, 23 1-241. 3. Eringen, A. C. and Ingram, J. D., A continuum theory of chemically reacting media-I. International Journal of Engineering Science, 1967, 3, 197-2 12. 4. Biot, M. A., Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 1955, 26, 182-185. 5. Prevost, J. H., Mechanics of continuous porous media. International Journal of Engineering Science, 1980, 18, 787-800. 6. Prevost, J. H., Consolidation of anisotropic porous media. Journal of Engineering Mechanics Division, AXE, 198 1, 107(EMl), 169-l 86. 7. Prevost, J. H., Nonlinear transient phenomena in saturated porous media. Computer Methods in Applied Mechanics and Engineering, 1982, 20, 3-l 8. 8. Prevost, J. H., Implicit-explicit schemes for nonlinear consolidation. Computer Methods in Applied Mechanics and Engineering, 1983, 39, 2255239. 9. Dienes, J. K., On the analysis of rotation and stress rate in deforming bodies. Acta Mechanica, 1979, 32, 217-232. 10. Dafalias, Y. F., Corotational rates for kinematic hardening at large plastic deformations. Journal of Applied Mechanics, ASME, 1983, 56 l-565. 11. Kiousis, P. D. and Voyiadjis, G. Z., Lagrangean continuum theory for saturated porous media. Journal of Engineering Mechanics, 1985, lll(lO), 1277-1288. 12. Meyerhoff, G. G., The ultimate bearing capacity of wedge-shaped foundation. Proceedings 5th International Conference on Soil Mechanics and Foundations, Vol. 2, 1961, pp. 103-109. 13. Torstensson, B. A., Pore pressure sounding instrument, Proceedings AXE Speciality Conference on In situ Measurement of Soil Properties, Vol. 2, ASCE, New York, 1975, pp. 48-54. 14. Torstensson, B. A., The pore pressure probe, Paper 34, Geotechnical meeting, Norwegian Geotechnical Society, Oslo, Norway, 1977, pp. 34.1-34.15. 15. Vesic, A. S., Expansion of cavities in infinite soil mass. Journal of Soil Mechanics, ASCE, 1972, 98(3), 265-290. 16. Levadoux, J. N. and Baligh, M. M., Pore Pressures during Cone Penetration in Clays, Report No. R.80- 15, Dept. of Civil Engineering, MIT, 1980, 3 10 pp. 17. Levadoux, J. N. and Baligh, M. M., Consolidation after piezocone penetration. I: prediction. Journal of Geotechnical Engineering, ASCE, 1986, 112(7), 707-726. 18. Houlsby, G. and Teh, C. I., Analysis of the piezocone in clay, Proceeding of ISOPT I, Vol. 2, A. A. Balkema, Rotterdam, The Netherlands, 1988, pp. 777783. 19. Gupta, R. C. and Davidson, J. L., Piezoprobe determination coefficient of consolidation. Soils and Foundation, 1986, 26(3), 12-22. 20. Elsworth, D., Analysis of piezocone dissipation data using dislocation methods. Journal of Geotechnical Engineering, ASCE, 1993, 119(10), 1601-1623. 21. Kiousis, P. D., Voyiadjis, G. Z. and Tumay, M. T., A large strain theory and its application in the analysis of the cone penetration mechanism. International Journalfor Numerical and Analytical Methods in Geomechanics, 1988, 12, 45-60.
222
G. Z. Voyiadjis and M. Y. Abu-Farsakh
22. Voyiadjis, G. Z., Kurup, P. U. and Tumay, M. T., Preparation of large size cohesive specimens for calibration chamber testing. Geotechnical Testing Journal, ASTM, 1993, 16(3), 339-349. 23. Kurup, P. U., Voyiadjis, G. Z. and Tumay, M. T., Calibration chamber studies of piezocone tm:st in cohesive soils. Journal of Geotechnical Engineering, 1994, 120(l), 81-107. 24. Bathe, K. J., Finite Element Procedures in Engineering Analysis. Prentice-Hall, Inc., New Delhi, India, 1990. 25. Voyiadjis, G. Z. and Kattan, P. I., Eulerian constitutive model for finite deformation plasticity with anisotropic hardening. Mechanics of Materials, 1989, 7, 279-293. 26. Zienkiewicz, 0. C., The Finite Element Method. McGraw-Hill Publishing Co., Inc., New York, 1977. 27. Lambe, T. W. and Whitman, R. V., Soil Mechanics. John Wiley, Inc., New York, 1969. 28. Lee, I. K., Ingles, 0. G. and White, W., Geotechnical Engineering. Pitman Publishing, Inc., Marshfield, MA, 1983.
Computers and Geotechnics, Vol. 20, No. 3/4, pp. 195--222. 1997 Q 1997 Elsevier Science Ltd. All rights reserved Printed in Great Britain SO266-352X(97)00003-7 0266-352X/97 S 15.00 + 0.00
ELSEVIER
Coupled Theory of Mixtures for Clayey Soils G. Z. Voyiadjis Department
of Civil and Environmental
& M. Y. Abu-Farsakh
Engineering, Louisiana LA 70803, U.S.A.
State University,
Baton
Rouge.
ABSTRACT In this work, elasto-plastic coupled equations arc formulated in order to describe the time-dependent deformation of saturated cohesive soils (twophase state). Formulation of these equations is based on the principle of virtual work and the theory of mixtures for inelastic porous media. The theory of mixtures for a linear elastic porous skeleton was$rst developed by Biot (Theory of elasticity and consolidation for a porous anisotropic solid, Journal of Applied Physics, 1955, 26, 188-185). An extension of Biot’s theory into a nonlinear inelastic media was performed by Prevost (Mechanics of continuous porous media, International Journal of Engineering Science, 1980, 18, 787-800). The saturated soil is considered as a mixture of two deformable media, the solid grains and the water. Each medium is regarded as a continuum andfollows its own motion. Thejow, of pore-water through the voids is assumed to follow Darcy ‘s law’. The coupled equations are developed for large deformations with finite strains in an updated Lagrangian reference frame. The coupled behavior of the tlz’ophase materials (soil-water state) is implemented in a jnite element program. A modtfied Cam-clay model is adopted and implemented in the$nite element program in order to describe the plastic behavior of clayey9 soils. Penetration of a piezocone penetrometer in soil is numerically simulated and implemented into a finite element program. The piezocone penetrometer is assumed to be infinitely st@ The continuous penetration of the cone is simulated by applying an incremental vertical movement of the cone tip boundary. Results of the finite element numerical simulation are compared with experimental measurements conducted at Louisiana State University using the calibration chamber. The numerical simulation is carried out for two cases. In the jrst case, the interface friction between the soil and the piezocone penetrometer is neglected. In the second case, interface friction is assumed between the soil and the piezocone. The results of the numerical simulations are compared with experimental laboratory measurements. 0 1997 Elsevier Science Ltd.
195
196
G. Z. Voyiadjis and M. Y. Abu-Farsakh
INTRODUCTION The soil skeleton consists of solid grain particles and voids that are filled with water and/or air. Soil can be considered as a mixture of multiphase deformable medium of solid grains, water, and air, each of which is regarded as a continuum and follows its own motion. The general theoretical framework that forms the basis of the theory of mixtures was first developed by Truesdell and Toupin [l]. Further development of the modern theory of mixtures was presented by Green and Naghdi [2] and Eringen and Ingram [3]. The theory of mixtures for a small strain linear elastic porous media was first developed by Biot [4] assuming that fluid flow through solid phase obeys Darcy’s law. An extension of Biot’s theory for large strain nonlinear inelastic porous media was presented by Prevost [5]. The theory of mixtures was developed in a Eulerian reference frame. The formulation by Prevost [5%3] was also developed in a Eulerian reference frame. It was shown that the use of Jaumann stress rate in a Eulerian formulation may cause inaccuracies in large strain problems [9,10]. Kiousis and Voyiadjis [ 1 I] developed a theory for two-phase material based on the concept of the theory of mixtures in a Lagrangian reference frame using the material stress rate. In this paper, an elasto-plastic analytical model is formulated in order to describe the nonlinear elasto-plastic, time-dependent deformation of a saturated cohesive soil skeleton (two-phase state). Formulation of the coupled equations for the model is based on the principle of virtual work and the theory of mixtures for inelastic porous media. The theory for large strain nonlinear inelastic porous media used here is based on the work of Prevost [5] and Kiousis and Voyiadjis [ 111. However, the coupled equations are formulated for large deformations using finite strains in an updated Lagrangian framework. Recently, there has been an increased concern toward the evaluation of different engineering soil parameters using in situ testing. Piezocone penetrometer tests (PCPT) are one of the most widely used for in situ testing. Many interpretation procedures based on bearing capacity model [ 121, cavity expansion theories [13-l 51, strain path method [ 166181, semi-empirical methods [19], dislocation methods [20], and the finite element techniques [21], have been proposed by several investigators to analyze the piezocone penetration and to evaluate some engineering soil parameters. Most of the methods that are available assume certain deformation patterns and small strains that do not match with the experimental results. Strains are reported to exceed 100% [I 6,17,21]. To accurately analyze the PCPT, it is necessary to have a numerical model that takes into consideration the large deformations, nonlinear soil behavior, and changing of the boundary condition with penetration. The piezocone penetration test results depend on many factors,
Theory
of mixtures
for clayey soils
197
including soil type and compressibility characteristics, stress history, and the presence of fissures and cracks. To validate the numerical model, it is necessary to compare it with controlled laboratory calibration tests of known soil types, characteristics, and stress history. The proposed analytical model is used for the analysis of the piezocone penetration problem. A modified Cam-clay model is adopted to describe the plastic behavior of the clayey soil. The penetration of the piezocone penetrometer in cohesive soils is numerically simulated and implemented in the finite element program. The continuous penetration of the piezocone is simulated by applying an incremental vertical movement of the cone tip boundaries. The problem is treated as a boundary problem that includes changing of the boundary condition as penetration proceeds. The numerical simulation is performed for two cases. In the first case, no interface friction is assumed between the soil and the piezocone penetrometer while, in the second case, interface friction is assumed between the soil and the piezocone penetrometer. A coefficient of interface friction equal to 0.25 is used in this work, which corresponds to the angle of friction of 6= 14” between the soil and the piezocone. The results of the numerical simulation are compared with experimental laboratory measurements of the miniature PCPT performed on a clay calibration chamber conducted at Louisiana State University (LSU) [22,23].
COUPLED
FIELD
EQUATIONS
The coupled field equations proposed herein are formulated in an updated Lagrangian framework. The principle of virtual work is expressed here using the updated Lagrangian formulation given by Bathe [24] as follows:
s
;+‘&sS( ;+‘E,@‘%’ = “+‘R
(1)
“V
Ef16is the total strain and is decomposed into the linear strain, #As, and the nonlinear strain nqAs: n+l
,,
CAB =
&AB
+
nT]AB
nUA,B
+
db3.A)
+
nuK,B)
(2)
where neAB
=
;1(
1 nl]AB= 2 (nUK,A
(4)
G. Z. Voyiadjis and M. Y. Abu-Farsakh
198
where u is the displacement vector. In Eqn (1) :+~s.Q is the second PiolaKirchhoff stress tensor at n + 1 referred to n configuration as indicated in Fig. 1. It is expressed in terms of the Cauchy stress tensor “0~s such that
:+‘SAB= naAB+
(5)
naiad
where
naAB= In Eqn (5) “ASAB is the y1+ 1 as indicated in Fig. pore water pressure and external virtual work due as follows:
azB+ Pd.4~
(6)
incremental
stress tensor from configuration n to Cauchy stress tensor P, is the GABis the Kronecker delta. In Eqn (l), "+'Ris the to the applied loads and tractions and is evaluated 1. aiB is the effective
The relation between the coordinates solid medium is given by
at n + 1 and n configurations
I
for the
Yo ,Y. ,Yn+l
/0
Vo+l
0 V.
n+l configuration
n configuration
/ 0 configuration vo D
/
zo,Z”
Fig. I. Configuration
,Z”hl
of a body in motion
x0
,xn,xn+1
with a fixed Cartesian
coordinate
system
Theory of mixtures for clayey soils
199
Zi = n+lq = “Qy+ S&
(8)
The relation between the Cauchy stress tensor and the second Piola-Kirchhoff stress tensor is given by
where
and Js is the corresponding Elasto-plastic
Jacobian for solids.
constitutive equation
For finite deformations, assumed to be
the constitutive o’s (3
ab
=
equation for the solid skeleton is
%cdds,d
where Ds is the elasto-plastic stiffness matrix. The effective co-rotational stress rate tensor 0’s is given by [25] as
where a’” is the effective Cauchy stress rate tensor, 0’” is the effective Cauchy stress tensor, and
is the modified spin tensor. In Eqn (12a)
is the total spin tensor, and J+$Lis the plastic spin tensor. When I%$ = 0, the co-rotational stress rate tensor reduces to the Jaumann stress rate tensor. The deformation gradient ds is given as follows:
200
G. Z. Voyiadjis and M. Y. Abu-Farsakh (13)
The relation co-rotational
between the total co-rotational stress rate tensor is given by
where P, is the material
derivative
stress
rate and the effective
of P, given as V”.gradP,
By differentiating both sides of Eqn (9) with respect gian reference frame, one obtains
(15) to time in the Lagran-
Making use of expression (13) for the spatial deformation gradient (16) together with Eqns (6) and (1 l)-(14) one obtains the following sion:
Substituting
in Eqn expres-
Eqn (10) into Eqn (17) and using the relation
(18) one reduces
Eqn (17) into the following
If one neglects the plastic spin tensor relation: ,AB
=
DlZ,,cDic~
relation:
@,
+
Eqn (19) reduces to the following
px”,,,~B,b~w~ab
(20)
Theory of mixtures for clayey soils
201
where 1 is the strain rate tensor. In Eqn (19) D* is given as follows:
The left-hand side of Eqn (1) may be expressed as follows: ( naAB +
&AB)J(
neAB f
%AB)dnV
=
“V
“OAB~( @AB
+
“VAB)~~~
“V (22)
&SAB6(
+
neAB +
nl]AB)dnV
1 “V
However, Ans~a may be expressed as follows: [+A/ AnsAB
t+At iABdt
=
=
D&D
ic-dt
s
J
t
t
t+At +
A ~SAB
Substituting
=D&D
(23) &dt
‘?&%,a
AECD
+ f~A,a%,#%i
Eqn (23) into Eqn (22), one obtains the following expression:
DiBcD(A &CD
+
A nWD)&eABdp
-I
“V
DiBCDA
&‘CD&VABdnV
“V
+
J
(“&j
+
“V
=
%~AB)klABdnV+
s
JSXSA,,XB,bsabAPw(6neAB
“V
( nc7;B + nPw8A&,eABdnV
"+'RI
“V
Equation (24) is the first derived coupled equation.
+ &‘lAB)dnV
202
G. Z. Voyiadjis and M. Y. Abu-Farsakh
Theory of mixtures The saturated soil is able media, the solid uum and follows its (assuming chemically
considered here as a mixture consisting of two deformgrains and water, each of which is regarded as a continown motion. The balance of mass for each constituent inert) can be written as [5] (a = s or w)
ba + p”divv” = 0
(25)
where p” is the macroscopic mass density of the a! constituent. If one assumes both solid grains and water to be incompressible materials, then the balance of mass can be written in terms of the soil porosity n’+‘:,solid velocity vS and water velocity v”‘ as follows [5]:
div(v’) - div(v”‘) = $ [div(u’) + (v”’ - v”)grad(n”‘)] The general form of Darcy’s ium is given by
law for the flow of water through
(vl - 9) = - -&K”“(grad(P,)
- p,b)
YW
porous
med-
(27)
where K”‘” is the permeability tensor, b is the body force vector, and yW is the unit weight of water. Taking the divergence of both sides of Eqn (27) and substituting it in Eqn (26) together with Eqn (18) and making use of the following relations:
apw
x”-
a -
=
a.22
one then obtains the second formulation as follows:
apw
(284
Csbaxc
azb -=
XD,, g
coupled
D
equation
in an updated
Lagrangian
Theory of mixtures for clayey soils
NUMERICAL
203
FORMULATION
The finite element discretization is used here for the displacement pore water pressure P, as follows: u=M.U
u and the
(30a)
P,=N.W
(30b)
where M is the displacement shape function, N is the pore water pressure shape functions, U is the nodal displacement, and W is the nodal pore water pressure. The incremental strains can be expressed as follows: 6e=
BL.GU
(3la)
Jr] = BNL. 6U
and the pore water pressure
gradient
ap,
-=
axB
are given such that N,BW
where BL and BNL are the linear and nonlinear strain-displacement Substituting Eqn (31) into Eqn (24) one obtains ,,KAU +” S-JAW= ,,a where .K =
is the elasto-plastic
.KL
+
stiffness matrix,
&a_
(3lb)
+
nKzL + “KS
(3lc)
matrices.
(32)
204
G. Z. Voviadjis and M.
nKL
s
=
Y. Abu-Farsakh
&D*
,,BLd”V
&I*
,,BNLdnV
”v is the linear stiffness matrix,
nKNL
=
s
”v is the nonlinear
stiffness matrix.
Jc =
.I
,,B~,%,,B~,d”V
“V
is the geometric
is
the coupling
stiffness matrix,
matrix,
and
The condition that the continuity Eqn (29) applies throughout nuum using Galerkin’s weighted residual method requires that
the conti-
(33) The weak form of Eqn (33) is obtained this equation:
s
J’C$ i,P,d”V
I
”v
by applying
g J”C;; C&l l$d”B(s+p,Ba)
Green’s
theory
[26] to
$d”V
“V (34) -
c J “S
q,,P,d
“A = 0
Theory of mixtures for clayey soils
205
where &,=NG
(35a)
is the weighted residual (virtual pore pressure)
aL
z=
w
N
(35b)
,A
qn
is the seepage velocity normal to the boundary. Substituting Eqn (3 1) and Eqn (35) into equation Eqn (34), one obtains -,,fl=AU
+ ,,\k&AW = J-I
(36)
where
is the flow matrix, and “II = -6tG
- Gt\kW” +
J
NTq,dnA
“S
From Eqn (32) and Eqn (36), the complete coupled behavior of the twophase skeleton-fluid state can be expressed as follows: (37)
where AU is incremental nodal displacement, nodal excess porewater pressure.
NUMERICAL
and A W is the incremental
EXAMPLES
A number of examples are analyzed to demonstrate the accuracy of the analytical model described. In the analysis, the eight node quadratic isoparametric finite element (Q8P4) is used.
G. Z. Voyiadjis and M. Y. Abu-Farsakh
206
One-dimensional
elastic consolidation
The one-dimensional elastic consolidation problem is studied here. The finite element mesh used in the analysis of this problem and the corresponding boundary conditions are shown in Fig. 2. It consists of 20 Q8P4 finite elements with an H/B ratio of 5. The traction q is applied instantaneously to the top surface at time t = 0. Thereafter, the traction is held constant with drainage allowed only at the top surface. A closed form analytical solution for the small deformation case as a function of time can be found in Lambe and Whitman [27]. The average degree of consolidation as a function of the dimensionless time factor (T = Cvt/H2) is obtained from the finite element analysis for a load ratio of q/E = 0.1. This is compared with the analytical solution as shown in Fig. 3; q is the load intensity, E is the Young’s modulus, and C, is the coefficient of consolidation. The finite element results show good agreement with the analytical solution. Figure 4 represents the settlement of the top surface 6 as a function of the dimensionless time factor (T), for various q/E ratios (0.1 and 0.5). Figure 4 shows that the settlements obtained from the updated Lagrangian finite strain deformation analysis are lower than those obtained from small deformation analysis. It is noted that, as the q/E load ratio increases, the
Drainage
Fig. 2. One-dimensional
elastic consolidation:
finite element
mesh
Theory of mixtures for clayey soils
207
difference between the small and the finite deformation analysis also increases. This is expected since the large deformation solution is stiffer than the small deformation case.
0
,
0.2 -
e
.3 s ‘i; 9
0.4 -
a 8 % :
0.6 -
rQ B 0.8 -
Analytical 0
Finite Elements
0.1
0.01
0.001
Time
Fig. 3. One-dimensional
1
Factor (T)
elastic consolidation:
i
\
degree of consolidation.
\
WE - 0.6 \ ,
0.3 -
\ -
__
Small Deformation Updated
Fig. 4. One-dimensional
Lagrangian
_-__-_-_
Lqe
Deformation
elastic consolidation:
settlement.
G. Z. Voyiadjis and M. Y. Abu-Farsakh
208
Two-dimensional
elastic consolidation
The finite element mesh and the corresponding boundary conditions for the problem are shown in Fig. 5. A uniform distributed strip load of intensity q and width 2B is applied instantaneously to the top surface at time t = 0. Thereafter, the load is held constant with drainage allowed only at the top surface. Figure 6 shows the ratio of the computed excess pore pressure at 0.5B beneath the center of the footing to the initial excess pore pressure (t = 0) as a function of the dimensionless time factor (0. The finite element analysis results are compared with the closed form solution of the excess pore pressure for an infinite elastic half-space [28]. During the early times of consolidation, there is a slight increase in the pore pressure as a result of the particular distribution of the pore pressure. This phenomenon is known as the Mandel-Cryer effect. A good agreement is obtained between the numerical solution and the analytical solution, as shown in Fig. 6. The computed settlement of the point located at the center of the strip load 6 as a function of the dimensionless time factor T for various q/E ratios (0.1 and 0.3) are shown in Fig. 7. Similar behavior to the one-dimensional elastic consolidation can be seen here for the two-dimensional elastic consolidation. It is noted that, as the q/E load ratio increases, the difference between small and finite deformation analysis also increases.
+tti
Fig. 5. Two-dimensional
Drainage
elastic consolidation:
finite element
mesh.
Theory of mixtures for clayey soils
0.4
-
209
Analytical
0.01
0.1
Time Factor Q
Fig. 6. Two-dimensional
elastic consolidation:
pore water pressure at depth z =0.5B.
0.4
0.8
- - --
\
Small Deformation Updated Lagrangian Large Deformation
1.2
I
’ ’ “““I
0.001
0.01
’ “““‘I 0.1
’ “‘m’l
Time Factor
Fig. 7. Two-dimensional
elastic consolidation:
1
\
\
\
’ ’ “““~ 10
\
\
\
\
’ ’ “In 100
(T)
settlement of the mid-point.
210
G. Z. Voyiadjis and M. Y. Abu-Farsakh
Piezocone penetrometer analysis In this work, the penetration of a piezocone penetrometer in cohesive soils is numerically simulated assuming that the piezocone penetrometer is infinitely stiff and no tensile stresses are allowed. The projected piezocone penetrometer area is 1 cm2 with a cone apex angle of 60” which is the same size as the cone used in the experiments by Kurup et al. [23] and Voyiadjis et al. [22]. Simulation is performed for two cases. The first case assumes that the soilpenetrometer interface friction is negligible. In case 2, it is assumed that the coefficient of the soil-penetrometer interface friction is 0.25, which corresponds to an angle of friction 6 =14” between the soil and the piezocone surface. The continuous penetration of the piezocone is simulated by applying an incremental vertical movement of the cone tip boundary, as shown in Fig. 8. For case one, the nodes along the stress boundary of the cone tip and along the piezocone shaft are allowed to slide from the beginning of the incremental penetration, while, for case two, the nodes along the inclined tip surface and along the piezocone shaft are first prevented from sliding along the surface until sliding potential occurs. The sliding potential is reached along the surface when friction forces of these nodes exceed the allowable friction forces. After that, sliding along the stress tip surface and along the shaft surface is allowed. During the incremental penetration for cases one and two (after sliding potential is reached), the nodes are checked for possible changes in
Position (n)
Fig. 8. Incremental
Position (n+l)
penetration
of the piezocone
penetrometer.
Theory of mixtures for c1aye.y soils
211
the boundary conditions. When tensile stresses occur in the node immediately below the lower end of the cone at location (l), the node is released as indicated by location (2) and allowed to move freely. Once the free node at location (2) reaches the cone tip boundary, its movement for case one will be restricted along the skew surface of the cone tip as indicated by (3). In case two, the nodes are first fixed and vertically incrementally displaced until the sliding potential is reached. Then they are allowed to slide along the surface. When the horizontal co-ordinate of the nodes at location (3) exceeds the penetrometer radius, these nodes are released at location (4) until they return to the boundary of the shaft. They are then restricted from further movement in the horizontal direction, as shown in location (5). The numerical simulation is compared with the PCPT conducted on cohesive soil specimen prepared in a calibration chamber at LSU. The soil specimen used is a K-33 prepared by consolidating a soil slurry of 33% kaolinite and 67% Edgar fine sand (D6,,/Di0= 1.4) at a water content of twice the liquid limit. The soil specimen has a liquid limit, WL = 20% and a plastic limit, W, = 14%. The specimen is consolidated to an isotropic effective stress of 207 kPa against a back pressure of 138 kPa. Full details about specimen preparation, test procedures and type of penetrometer used are found in Kurup et al. [23] and Voyiadjis et al. [22] The filter used to measure the porewater pressure is located either in the lowest quarter of the cone at the very tip (Ul configuration) or starting 0.5 mm above the base of the cone with 2 mm vertical height (U2 configuration). The modified Camclay model is adapted in this work to describe the plastic behavior of soil. The yield locus for modified Cam-clay is given by f(P,
q,
P&))
z M2P2 -
M2PcP+ cl*= 0
(38)
It represents an ellipse in the p - q plane. In Eqn (38) P and q are the mean effective and deviatoric stresses. PC is the strain hardening parameter representing the apex of the yield locus ellipse in the P axis. The soil parameters used for the modified Cam-clay model are the slope of the virgin consolidation line in the e-1nP space, h=O.ll; the slope of the unloading-reloading line in the e-1nP space, K =0.024; the slope of the critical state line, M =1.2; r = 1.162; the coefficient of permeability for the soil is 0.5 x 1O-9 m/s; and the Poisson’s ratio u =0.3. A number of trial finite element meshes are studied first with different degrees of refinement in order to obtain the appropriate mesh for the analysis of the piezocone penetration problem. Figure 9 shows the computer dimensionless quantities qc/qc5 and Ul /Ul 5 as a function of the number of finite elements used along the inclined surface (as a measure of degree of
212
G. Z. Voyiadjis and M. Y. Abu-Farsakh
refinement) for a penetration depth of 1Omm. qC/qC5 is the ratio of the cone tip resistance to the cone tip resistance computed using five elements along the inclined surface. Ul /Ul s is the ratio of the excess pore pressure at Ul configuration to the one computed using five elements along the inclined surface. From Fig. 9, it is clear that, with the increase of refinement, the numerical results converge to a unique solution. Figure 10 shows the finite element mesh used in the analysis of the penetration of the piezocone penetrometer in cohesive soils. The calculated cone tip resistance profiles during piezocone penetration and those obtained experimentally are compared in Fig. 11. In Fig. 12, the developed excess pore pressure profiles during piezocone penetration obtained numerically are compared with those obtained from experimental testing. It can be seen that the cone resistance obtained from the numerical simulation increases rapidly until it reaches a steady-state condition at a depth of 25-30mm, as shown in Fig. 11. The same behavior is also obtained for the development of excess pore pressure. The computed cone tip resistance and excess pore pressures develop at faster rates than those obtained experimentally. This difference may be due to the errors encountered at the initial experimental stage of testing and also to the lack of confining stresses at the top surface of the soil specimen. This can be seen if one compares the
2 Number Fig. 9. Convergence
3 of Elements
4 5 along the Inclined
6 Surface
of the cone time resistance and pore water pressure with mesh refinement.
of Ul configuration
Theory of mixtures for clayey soils
213
two cone tip resistance test profiles. The only difference between them is the filter location. The final values obtained from the computational model after the steady-state condition is reached are quite close to those obtained experimentally. In case 1, in which no friction is assumed between the soil and the piezocone surface, the cone tip resistance value at the steady state is
Fig. 10. Finite element mesh of the piezocone
P
r C = :
-60.00
penetration
problem.
--._--_._ 1.---
Finka Elemmt IFrktionl
I
i 0.00
0.40
0.84
Cone Resistance
Fig. 11. Cone resistance
1.20
1 .BO
(MPa)
profiles during piezocone
penetration.
214
G. Z. Voyiadjis and M. Y. Abu-Farsakh
0.60 0.40 0.20 Excess Pore Pressure (MPa)
0.80
(a)
1
.200.00 1
-
Fin::
I 0.00
0.60 0.40 0.20 Excess Pore Pressure (MPa)
I
1
0.80
(b)
Fig. 12. Excess pore pressure
profiles during piezocone penetration; configuration U2.
(a) configuration
Ul, (b)
1.131 MPa (as shown in Fig. 11). In case 2, where the soil-piezocone coefficient of interface friction p is taken as 0.25, the cone tip resistance at the steady state is 1.294 MPa. The cone resistance obtained experimentally is 1.255 MPa. The computed excess pore water pressures at the steady state for configuration Ul (at the lowest quarter of the cone tip) are 0.565 and 0.774 MPa for cases 1 and 2, respectively, as shown in Fig. 12(a). However, the
Theory of mixtures for clayey soils
215
measured excess pore water pressure at the steady state is 0.608 MPa. The computed excess pore water pressures for the steady state in configuration U2 (above the cone base) are 0.532 and 0.542 MPa for cases 1 and 2, respectively, as shown in Fig. 12(b). The corresponding experimental measured excess pore pressure at the steady state is 0.593 MPa. Figure 13 shows the dissipation of excess pore water pressure with time. Figure 13(a) compares the dissipation of excess pore pressure for configuration
h
2 i+
P
0.6 -
t!
2 0 8
ti
0.4 -
0 t i
2c
0.2 -
-
---.0.0
,
, ,,,,,,,
,.O
,
,
,,,,,I,
,o.o
/
100.0
I
,,,111/
I
I
In-n
10000.
1000.0
Time (set) (a)
_ -
1.0
~
PcPTul
---
Finit. Element (NO Friction)
-.
-
(Kurup
dill.
lSs4,
Finita Ekment(FMtion) . .
10.0
100.0
Time (set) (b)
Fig. 13. Dissipation
of excess pore pressures;
(a) configuration
U1, (b) configuration
U2.
G. Z. Voyiadjis and M. Y. Abu-Farsakh
216
Ul obtained both numerically and experimentally. Figure 13(b) compares the dissipation of the excess pore water pressures for configuration U2. It is obvious that the finite element analysis results are close to those obtained experimentally with the exception that there is no sudden drop of the pore pressure when the penetration stops. The distribution of the excess pore water pressures developed around the piezocone penetrometer for a penetration depth of 25mm are shown in Fig. 14. For the case with no interface friction (Fig. 14(a)), the maximum concentration excess pore water pressure is obtained in the upper third of the cone, while for case 2 with the coefficient of interface friction p of 0.25, the maximum concentration of excess pore water pressure is around the lower of the cone as shown in Fig. 14(b). In both cases, negative pore water pressure is developed at about 2-3 diameters above the cone base. Figures 15-18 present the contours of the effective stresses around the piezocone penetrometer for a penetration of 25 mm. The initial stress state is isotropic with ai = 0: = ok = 207KPa, and trz = 0. The radial effective stress ai around the piezocone penetrometer for case 1, for which friction is neglected, increases in the region above the mid-point of the cone up to 424.1 KPa immediately above the cone base. It then drops to 48.3 KPa below the cone tip, as shown in Fig. 15(a). For case 2, in which interface friction is -
Normalized Radial Distance 14. Contours
1
rho
Normalized
of excess pore pressures around the piezocone of 25 mm; (a) no friction, (b) coefficient of interface
Radial Distance
rho
penetrometer for penetration friction p = 0.25.
Theory of mixtures for clayey soils
20.0
217
J
L P
20.0
24.0
0.00
¶.oo
4.00
0.00
1
24.0
(
I
0.0
Normalized Radial Distance rho
2.0
I
60
0.0
1
c.10
I
Normalized Radial Distance rho
Fig. 15. Contours of radial stresses 0: around the piezocone penetrometer for penetration 25 mm; (a) no friction, (b) coefficient of interface friction p = 0.25.
Normalized Radial Distance
rho
Normalized Radial Distance
of
rho
Fig. 16. Contours of axial stresses ai around the piezocone penetrometer for penetration 25 mm; (a) no friction, (b) coefficient of interface friction w = 0.25.
of
G. Z. Voyiadjis and M. Y. Abu-Farsakh
218
considered, the stress increases up to 405.7 KPa and drops down to 54.3 KPa, as shown in Fig. 15(b). The radial stress bulb region of the zone influenced by penetration extends to approximately four times the piezocone diameter. The axial effective stress a: arising from penetration is presented in Fig. 16. The stress increases in the region below the cone face and decreases in the region above it. The maximum stress at the end of a 25mm penetration is located below the cone tip, while the minimum stress is located above the base of the piezocone. The stress ranges from 305.5 to 62.1 KPa and from 347.3 to 45.6 KPa for cases 1 and 2, respectively. Figure 17 shows the distribution of the tangential stresses aa around the cone. The stresses decrease from 207 KPa (initial stress) to 43.4 and 47.6 KPa, respectively, for cases 1 and 2. The minimum tangential stress (maximum decrease) is located at the lower half of the cone. In Fig. 18, the developed shear stresses t ,.Z are presented. The maximum shear stresses induced are -166 KPa for case 1 and -201.6 KPa for case 2. The shape of the shear stress distribution around the piezocone is in the shape of a bulb with the maximum value located at the cone face and oriented at an angle of approximately 45” to the cone face. The shear stress bulb is also extended
0.00
2.00
Normalized
4.00
0.00
Radial Distance
0.00
rho
1-r 2.0
0.0
Normalized
Fig. 17. Contours of tangential stresses CJ~around the piezocone of 25 mm; (a) no friction, (b) coefficient of interface
Radial Distance
rho
penetrometer for penetration friction p= 0.25.
Theory of mixtures for clayey soils
c =
219
32.
p” P z
30.
6 z 23.
24.0 I 0.0
Normalized Radial Distance rho
, 2.0
I 4.0
I 6.0
Normalized Radial Distance rho
Fig. 18. Contours of shear stresses r,, around the piezocone penetrometer for penetration 25 mm; (a) no friction, (b) coefficient of interface friction p = 0.25.
of
behind the cone face. This extension is bigger for case 2 than it is for case 1. For case 2, there is a shear stress concentration around two points. One is at the lower face of the cone, and the other is at the upper portion of the cone face near the cone base.
CONCLUSIONS Coupled elasto-plastic equations are formulated in order to describe the time-dependent deformation of saturated cohesive soils. The formulation of these equations is based on the principle of virtual work and the theory of mixtures for inelastic porous media. The coupled equations are implemented in a finite element program. The modified Camclay constitutive model is used to describe the plastic behavior of the clayey soils. Penetration of the piezocone penetrometer is numerically simulated and implemented into the finite element program. Simulation is carried out for two cases. In the first case, the interface friction between the soil and the piezocone penetrometer is neglected. In the second case, the soil-piezocone penetrometer interface
220
G. Z. Voyiadjis and M. Y. Abu-Farsakh
friction is considered. Results obtained from the numerical simulations are compared with the experimental laboratory measurements conducted at the LSU calibration chamber. Good correlation is obtained between the experimental results and those simulated numerically using the proposed model. The following conclusions can be made from the analysis of the piezocone penetration. The steady-state condition is reached in the numerical analysis at a depth of 25 to 30mm, which is at a faster rate than those obtained experimentally. This may be attributed to some errors encountered at the initial experimental stage of testing and also due to the lack of confining stresses near the top surface of the soil specimen. The pore pressure distribution around the cone tip is concentrated around the upper third of the cone if friction is neglected and around the lower third of the cone if interface friction is taken into consideration. In addition, a negative pore pressure may be generated above the cone base due to soilpiezocone separation in this area, especially at the early stages of testing. This raises the question as to where is the best location to place the transducer in order to avoid possible misleading calculations. The authors believe that the middle third is the best location to place the transducer because it is the most representative of the average pore pressure which is more uniform and less sensitive. The bulb of stress changes surrounding the cone tip is extended radially outward to about four times the piezocone diameter. The stresses change faster in the vicinity of the cone tip. This makes the problem very localized around the cone tip, which may lead to incorrect calculation due to the high stress redistribution in this area. The effect of considering the interface friction between the soil and the piezocone penetrometer is paramount for the shear stress distribution around the cone tip and the generated excess pore pressures close to the cone surface, which has an affect on the calculated cone tip resistance. The effect is minimal for other stress distributions.
ACKNOWLEDGEMENTS This research was supported by grant MSS-9312707. The authors encouragement of Dr Priscilla P. Directorate for Engineering, NSF. portation Research Center under edged.
the National Science Foundation under wish to acknowledge the support and Nelson, Program Director of G3S/CMS, Partial support of the Louisiana TransTask Order No. 736-99-0293 is acknowl-
Theory of mixtures for clayey soils
221
REFERENCES 1. Truesdell, C. and Toupin, R., The classical field theories. In Hun&& der Physik, Vol. 3, No. 1, ed. S. Flugge. Springer-Verlag, Berlin, 1960. 2. Green, A. E. and Naghdi, P. M., A dynamical theory of interacting continua. International Journal of Engineering Science, 1965, 3, 23 1-241. 3. Eringen, A. C. and Ingram, J. D., A continuum theory of chemically reacting media-I. International Journal of Engineering Science, 1967, 3, 197-2 12. 4. Biot, M. A., Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 1955, 26, 182-185. 5. Prevost, J. H., Mechanics of continuous porous media. International Journal of Engineering Science, 1980, 18, 787-800. 6. Prevost, J. H., Consolidation of anisotropic porous media. Journal of Engineering Mechanics Division, AXE, 198 1, 107(EMl), 169-l 86. 7. Prevost, J. H., Nonlinear transient phenomena in saturated porous media. Computer Methods in Applied Mechanics and Engineering, 1982, 20, 3-l 8. 8. Prevost, J. H., Implicit-explicit schemes for nonlinear consolidation. Computer Methods in Applied Mechanics and Engineering, 1983, 39, 2255239. 9. Dienes, J. K., On the analysis of rotation and stress rate in deforming bodies. Acta Mechanica, 1979, 32, 217-232. 10. Dafalias, Y. F., Corotational rates for kinematic hardening at large plastic deformations. Journal of Applied Mechanics, ASME, 1983, 56 l-565. 11. Kiousis, P. D. and Voyiadjis, G. Z., Lagrangean continuum theory for saturated porous media. Journal of Engineering Mechanics, 1985, lll(lO), 1277-1288. 12. Meyerhoff, G. G., The ultimate bearing capacity of wedge-shaped foundation. Proceedings 5th International Conference on Soil Mechanics and Foundations, Vol. 2, 1961, pp. 103-109. 13. Torstensson, B. A., Pore pressure sounding instrument, Proceedings AXE Speciality Conference on In situ Measurement of Soil Properties, Vol. 2, ASCE, New York, 1975, pp. 48-54. 14. Torstensson, B. A., The pore pressure probe, Paper 34, Geotechnical meeting, Norwegian Geotechnical Society, Oslo, Norway, 1977, pp. 34.1-34.15. 15. Vesic, A. S., Expansion of cavities in infinite soil mass. Journal of Soil Mechanics, ASCE, 1972, 98(3), 265-290. 16. Levadoux, J. N. and Baligh, M. M., Pore Pressures during Cone Penetration in Clays, Report No. R.80- 15, Dept. of Civil Engineering, MIT, 1980, 3 10 pp. 17. Levadoux, J. N. and Baligh, M. M., Consolidation after piezocone penetration. I: prediction. Journal of Geotechnical Engineering, ASCE, 1986, 112(7), 707-726. 18. Houlsby, G. and Teh, C. I., Analysis of the piezocone in clay, Proceeding of ISOPT I, Vol. 2, A. A. Balkema, Rotterdam, The Netherlands, 1988, pp. 777783. 19. Gupta, R. C. and Davidson, J. L., Piezoprobe determination coefficient of consolidation. Soils and Foundation, 1986, 26(3), 12-22. 20. Elsworth, D., Analysis of piezocone dissipation data using dislocation methods. Journal of Geotechnical Engineering, ASCE, 1993, 119(10), 1601-1623. 21. Kiousis, P. D., Voyiadjis, G. Z. and Tumay, M. T., A large strain theory and its application in the analysis of the cone penetration mechanism. International Journalfor Numerical and Analytical Methods in Geomechanics, 1988, 12, 45-60.
222
G. Z. Voyiadjis and M. Y. Abu-Farsakh
22. Voyiadjis, G. Z., Kurup, P. U. and Tumay, M. T., Preparation of large size cohesive specimens for calibration chamber testing. Geotechnical Testing Journal, ASTM, 1993, 16(3), 339-349. 23. Kurup, P. U., Voyiadjis, G. Z. and Tumay, M. T., Calibration chamber studies of piezocone tm:st in cohesive soils. Journal of Geotechnical Engineering, 1994, 120(l), 81-107. 24. Bathe, K. J., Finite Element Procedures in Engineering Analysis. Prentice-Hall, Inc., New Delhi, India, 1990. 25. Voyiadjis, G. Z. and Kattan, P. I., Eulerian constitutive model for finite deformation plasticity with anisotropic hardening. Mechanics of Materials, 1989, 7, 279-293. 26. Zienkiewicz, 0. C., The Finite Element Method. McGraw-Hill Publishing Co., Inc., New York, 1977. 27. Lambe, T. W. and Whitman, R. V., Soil Mechanics. John Wiley, Inc., New York, 1969. 28. Lee, I. K., Ingles, 0. G. and White, W., Geotechnical Engineering. Pitman Publishing, Inc., Marshfield, MA, 1983.