- Email: [email protected]
Considerations on soil response to the rotation of principal stress directions
Computers and Geotechnics 8 (1989) 89-110
C O N S I D E R A T I O N S ON SOIL R E S P O N S E TO THE ROTATION OF P R I N C I P A l . STRESS DIRECTIONS
S Pietruszczak and S. Krucinski(*) Department of Civil Engineering and Engineering Mechanics McMaster University Hamilton, Ontario, Canada L8S 4L7
ABSTRACT The paper evaluates various plasticityformulations in the context of their ability to describe the influence of the rotation of principal stress axes. The discussion extends to a new theoretical framework which incorporates an implicit measure of material microstructure, derived from the directional distribution of porosity
INTRODUCTION In many practical geotechnical situations, the principal stress directions, in the vicinity of an engineering structure, deviate from those imposed by the depositional conditions. This m a y have a pronounced effecton the resulting displacement field.In recent years a number of researchers (e.g. [1],[2],[3])performed material tests which invariably indicate the sensitivity of soil response to the rotation of principal stress axes. Although the experimental information is stillfragmentary, the evidence gathered so far is very convincing (at least in a qualitative sense) and should not be ignored when formulating appropriate constitutive relations. Over the last two decades a large number of phenomenological models for geomaterials have been proposed [4].These models are mostly cast in the framework of theory of plasticity and they do not, in general, account for the influence of the rotation of principal stress axes. [n order to demonstrate this, let us examine closely the mathematical structure of some of the well established formulations. Consider an arbitrary loading history confined to pure rotation of the principal stress system. In other words, assume that a sample, under initial stress oijo,is subjected to the stress increment Ao..=T. T. o ° - o °. q ip Jq pq q
(1)
(*)On leave of absence from Technical University of Lodz, Poland 89
Computers and Geotechnics 0266-352X/90/S03.50 © 1990 Elsevier Science Publishers Ltd, England. Printed in Great Britain
90 Here, Aoij and o°ij are referred to the same Cartesian coordinate system fixed within the material and Tij is the usual transformation tensor. The stress path (1), when imposed on a sample, results in the deviation of principal stress axes from those corresponding to o°ij, while the principal values of o°ij are preserved.
Suppose at first, that the material is isotropic,
linearly elastic. In this case, the Hooke's law takes the form h~ii= hOii/3K;
heij= Asij/2G
(2)
where sijand eijrepresent the deviatoric parts of the stress and strain tensors and K,G are the elasticconstants (i.e.bulk and shear modulus, respectively). The orthogonality conditions on the directioncosines Tij,i.e.
Tip Tiq = 8pq
(3)
constrain the stress path (1) to satisfy Aoii = 0;
ASij = Aoij
(4)
so that, according to eqs. (2) A~ii = 0;
1 Aeij= ~-~(TipTjq-SipSjq) O°pq
(5)
Thus, the imposed path results in generation of deviatoric strain eij, while the volume change remains zero. It should be noted now, that the existing experimental evidence ([1],[2]} clearly indicates that the stress trajectories satisfying eq. (1) induce not only reversible but also permanent deformations of both distortional as well as volumetric nature. Let us examine now, in a similar context, the formulations derived from the theory of plasticity. The simplest of all, are the descriptions admitting isotropic hardening. In the most common approach, the evolution of the yield surface, f = 0, is linked to the generation of irreversible volumetric strain, {:Pii, so that f(oij, cp) = O;
0
o~o.. U
(6)
In such a ease, the consistency condition (7) aoij
~ii
results in
(8) Since f - 0 is assumed to be an i s o t r o p i c function of oij, then for all histories consisting of pure rotation of principal stress axes, eq. (1),
91 af ao.. u
6
ij
=o
~
~,=o
(9)
i.e.the loading process is a 'neutral'one and does not produce any plasticdeformations. The same conclusion can also be reached when employing multi-yield loci theories of soildeformation. In this case, the associated flow rule reads
--
(10)
n aq..
1] n
l]
where the summation includes only those terms which correspond to f n : 0 and f n : 0 . Let fm ( f, (n= 1, 2, ...., N) be expressed in a form fm(Oij ' ¢'P'n' cP} = 0
(II)
where ~P = (elfeiip)I/2.The consistency condition for any f m =
}m = __~m6.. + ~ ~ aoij q
__ +
_n
n [ ~ii /k;kk
0 takes the form
n_
= 0
(12)
00{ijl °~o{ij)
where afn/ao{ij)designates the deviatoric part of the gradient tensor. Denoting (afmDoij)oij by tkrn and the term contained in square brackets by l.lmn, a set of simultaneous equations (12) can be rewritten in the form
{~n} = [[.Imn]-I {~krn}
(13)
It is obvious that if {/~m}= {0} then, for arbitrary [Pmn], {~n}= {0}. The effects of deformation-induced anisotropy are often described by employing the framework of kinematic hardening. In a typical formulation the yield surface is assumed to undergo translation guided by a tensorial argument aij, i.e. frn(Oij, Uij, ePi) = 0
(14)
where aij=aij(eP~), i.e. the evolution of this surface is implied by the generation of plastic deformations. According to eq. (14) the consistency condition becomes o~f , o~f o~aij 0]"'+ - - - - ~ 1 aOiJ 0aij ~ 1
--
o~/ + --~k
=0
~ii
and results in a similar functional form to that of eq.(8),i.e.
(15)
92 1 ( of
)
~f
Of
of °°ij
of
(16)
In this case again, the pure rotation of principal stress directions describes a 'neutral' process, for which ~ = 0 and thus ~P ---0. Finally, another common description is that of ' b o u n d i n g s u r f a c e ' plasticity [5],[6] Here, for stress histories penetrating the interior of the bounding surface, the deformation process is described in terms of evolution of an inner loading surface
fl(oij, oij) = 0
(17)
The kinematics of this surface (i.e. direction of translation) is guided by a 'conjugate' stress tensor o¢ij located on the bounding surface. A typical form of the evolution law for the components of all is
aij ---- ~lq/ij ;
tl/ij
q/ij(Okl,Okl)
(18)
where l~ is a scalar parameter, determined via consistency condition
f"l =
Ofl ° Ofl dO.'--"~OiJ + 0ct..-- Vii ~ ---- 0 ~j
According to eq.(19), the case of
(19)
~j
(dfl/OOij)oi j -- 0 results in la = 0 and the surface fl = 0 remains
stationary. Therefore, once again, the formulation does not respond to the rotation of principal stress directions. It is evident from the above considerations, that none of the existing frameworks, based on classical concepts of isotropic/kinematic hardening, is capable of taking into account the sensitivity of plastic flow to the rotation of principal stress axes. Indeed, if the evolution of the yield surface (regardless of its form) results exclusively from the generation of plastic deformations, then these effects can never be properly embraced. There are several ways of enhancing the plasticity framework so that it responds to a change in orientation of the principal stress system. This has been accomplished, for example, in the formulation proposed in ref.[7] (see also ref.[8]), whereby the macroscopic behaviour of a sample has been deduced from deformation pattern along an infinite number of randomly oriented planes, by using an appropriate averaging procedure. An a l t e r n a t i v e approach, initiated in ref.[9], may be derived by enriching the set of classical functions defining the state of material by new tensorial function(s) which is/are capable of representing the soil fabric and its evolution. In such a case, the resulting mathematical structure remains phenomenological; the material response however, is an implicit function of the evolution of the fabric.
93 In the next section, the framework outlined in ref.[9] is reviewed in the context of its ability to describe the influence of the rotation of principal stress axes. Subsequently, the details concerning the i d e n t i f i c a t i o n of the t e n s o r i a l m e a s u r e of i n t e r n a l s t r u c t u r e a r e provided, followed by some numerical simulations p e r t a i n i n g to drained and undrained tests involving a continous rotation of principal stress directions.
PLASTICITY FORMULATION FOR ANISOTROPIC SOIL In ref.[9] an a t t e m p t was made to describe the anisotropy of soil structure in terms of a n i s o t r o p y in its p h a s e d i s t r i b u t i o n . As a n i m p l i c i t m e a s u r e of i n t e r n a l s t r u c t u r e , a 'directional porosity' has been employed, which could be identified as a generalized, direction dependent counterpart of 'porosity', i.e. the scalar valued q u a n t i t y defining the void space fraction. In order to construct a n average, non-singular measure of voids distribution, let us isolate in the vicinity of a m a t e r i a l point (P) a unit sphere (S), which encloses a representative volume of the material (V). Consider now a test line of length [ , = 2R with the orientation v with respect to a C a r t e s i a n coordinate system fixed at P. Assume t h a t l(v) represents the total length of interceptions of this line with soil pores. Thus, the fraction of L occupied by voids, L(v), can be defined as
L(v) = l(v)~ ;
l(v) = ~
Ii (v)
(20)
The m e a n value of the quantity L (v) averaged over the domain S, is 1 Lay = 4n
f
S
L(v) f(v)dS ;
II
4n
S
f(v)dS = 1
(21)
where tlv) is a scalar valued function describing the spatial distribution of test lines. It can been shown t h a t for uniformly distributed test lines, the first integral in eq.(21) is the measure
of average porosity no in the vicinity of P, whereas the lineal fraction occupied by pores is an unbiased estimator of the volume fraction of voids in the direction v, i.e. n o = Lay ,
n(v) - L(v)
(22)
The scalar valued function n(v), defined over the unit sphere S, can be represented by the generalized double Fourier series. The desired best fit approximation can be established by the 'least square' method leading to the r e p r e s e n t a t i o n in t e r m s of symmetric traceless tensors ~ij, flijkl, ...., cf.[10] n(vi) ~ no(1 + ~ij vi vj + ~ijkl vi vj v k v I + . . . . )
(23)
94 The higher rank tenors ~ijkl ..... relate to higher order fluctuations in void space distribution. Thus, in order to describe a smooth orthogonal anisotropy it is sufficient to employ an approximation based on the first two terms of the expansion (23). Incorporation of the concept of 'directional porosity', eq.(23), into plasticity framework requires the formulation of an appropriate evolution law for the components of fill. The latter should refer the rate of change of f~ij to the deformation history. In general, the components of flij may be related to the strain rate deviator eij, through an isotropic tensor valued function • • ° ~ij = ~ij(~kl ' ekl' no) (24) The experimental evidence suggests that for a certain class of geomaterials (e.g. clays, silts) the principal axes of ~ij may be considered as coaxial with those of the strain rate deviator. The above assumption results in significant simplification of representation (24), namely
~ij ---- ~ eij ;
~ = al(no) + a2(no) flik ~ki
(25)
where a's are scalar valued functions of the average porosity. Based on existing experimental evidence, it is now postulated that the deformation process alone, is not capable of inducing changes in the microstructure which would affect the functional form of the failure criterion F = 0 or result in anisotropy in the elastic range. The latter effects can only be attributed to the inherent anisotropy in the fabric of the material, which manifests itself in the size, shape and geometrical arrangement of the particles. In order to describe this anisotropy the concept of the fabric tensor Oil has been introduced, with aij defined as a tensor-valued function of ~ij, i.e. aij -- aij(flkl). Thus, in materials with large grain size (e.g. gravel with flat grains, sand with elongated grains, etc.),
Oil : Oij(~kl , akl) F = F(oij,
in elastic range
(26)
aij)= 0
whereas in materials with relatively small grain size the fabric remains isotropic, i.e. aij = 8ij. The proper mathematical representation of both the elasticity tensor and the failure criterion in terms of fabric tensor is again provided in ref.[9]. The anisotropy in the spatial distribution of porosity, as described by flij 2 0 , will
inevitably affect the plastic flow. The procedure for incorporation of the tensor ~ij into plasticity framework is not unique. One way of formulating the problem is to assume that f/ij affects directly the functional form of the yield criterion. Then, Of ~P" = ~ - ~J clo.. Ij •
f = f(°ij' ~'~ij K); '
(27)
95
i.e.,f= 0 depends on the deformation history (recorded by K) and ten (in general) functionally independent invariants of both tensors oij and f/ij.Now, it is interesting to note that the yield surface, defined by eq.(27), will undergo progressive evolution not only during an active loading process, but for all histories (including those for which f<0) associated with the rebuilding of material microstructure. Consequently, such a formulation will account for the influence of the rotation of principal stress axes. In order to demonstrate this, let us assume again z = ePi,and write the consistency condition in the form =
~Oij
+
dOij
-dfl.. -~ ij + - ~# de~ii
~k =
0
(28)
Utilizing the evolution law (25) for the components of flij, i.e. fiij = l]eiej %13~'p" = ~[3 ~ +l]~t Of 'J
2 G ij
(29)
o~(ij}
and substituting the above representation back in eq.(28),one obtains =
l(~f
+2Gt3 0t" )
H =
'
af
of
~
0f
(30)
It is evident from eq.(30) that for (df/aoij)6ij = 0 gP = Cijkl °kl ;
13 of of Cijkl -- 2GH doij aflkl
(31)
which implies that for histories experiencing pure rotation of principal stress axes, the progressive generation of plastic strain takes place, provided that (of/af~ij)fij > 0. A simpler description than that based on eq.(23) may be obtained by resolving the plastic strain rate into components governing a pseudo-isotropic response and a diversion from it. The flow rule can then be expressed as =
ao.. + 1~Gij ," q
f = f(oij , K) = 0
(32)
where the first term represents the strain rates generated in a fictitious isotropic medium, whereas Gij is a symmetric second order tensor whose components are function of f/ij and the deformation history, i.e. Gij = Gij(flkl, K). In the first approximation, one may simply take Gij = h f~ij;
~ = ~t
(33)
where h is a scalar valued function of the deformation history. Substituting eq.(33) into eq.(32) yields
96
•
=
+h~ij
- -
(34)
t~
The above flow rule is analogous to a non-associated law, i.e. according to eq.(34), the induced/inherent anisotropy results in a progressive deviation of the direction of plastic flow from that specified by the gradient tensor Although the formulation (34) is very attractive, in the sense of its numerical simplicity and performance, it does not respond to the rotation of principal stress axes. Indeed, when K is identified with ~P ~, i.e. plastic volumetric strain generated in a fictitious isotropic medium, then the consistency condition (7) results in a similar functional form to that of eq.(3) H\0o.. ij/; tj
H - -
a ~-~.. C~kk 11
(35)
It appears that the only rational way of enhancing the representation (34}, is to define the condition of the neutral state relative to the actual direction of plastic flow, i.e. --
- -
a(i.. l]
+
h
(36)
~ij
In such a case, it is evident that the stress trajectories resulting from pure rotation of principal stress axes will generate the plastic strain rates
.
q = Cijkl Okl ;
)
Cijkl =
+ h ~ij ~kl
(37)
U
provided h flij (~ij > 0. Finally, an alternative mathematical formulation may be obtained by admitting a coupling between both terms in the flow equation (32). Such a flow rule is then analogous to that derived in multi-yield loci theories, implying that the direction of plastic flow is sensitive to the direction of stress rate. Assume, for this purpose, that the components of Gij, eq.(33), are derived from a potential function Q = Q(~ij) = const. ;
~ij -- no(~ij + 8ij)
(38)
where ~ij is defined in such a manner that the equation (23) becomes n(v i) ~ ~ij vi vj. Let
Gij = h C~jkl 0~kl;
l~: O~ij q
(39)
97
where C"jkl is the elastic compliance tensor and fl= fl(~ij)=const, represents a family of loading surfaces defining the neutral states for structural rearrangement. Then, according to eq.(32),for la > 0 ~p •
=
~ - af JchC e ( aQ ~Ofl ) ~ doij ijkl 8~kl a~pcl Pq
(40)
It should be noted that in the above flow rule, the ptastie strain rates e'ij = ~ Gii, eq.(32), can be interpreted as being derived from residual stress rates Sail,
o' __ C e • R Cij ijkl Okl ;
* R
Okl =
h
- -
--
O~ Pq] O~ pq kl
(41)
so that such an approach is analogous to strain-space plasticity formulation. By inspecting eq.(40), it is again obvious that this framework is eapoble of responding to pure rotation of the principal stress system.
Under (df/doij)oij = 0, one obtains the
contribution from the second term in eq.(40), which does not vanish provided that l~ > 0. In conclusion, the reviewed constitutive relations are capable of describing the influence of the inherent as well as induced anisotropy on the deformation process. Their mathematical structure permits, in general, the modelling of the sensitivity of soil response to the rotation of principal stress axes. ON S T E R E O L O G I C A L E S T I M A T I O N OF P O R O S I T Y D I S T R I B U T I O N If the approximation of voids distribution, eq.(23), is confined to the representation n(v) -= no(1 + f~ij vi vj)
(42)
i.e. spherical harmonies of degree higher than two are disregarded, and if, according to equation (22) n(v) = L(v) then, the components of tensor f~ij can be estimated from the observation of L(v) on a limited number of cross-sections. In what follows, the procedure for the stereological estimation of the directional porosity distribution, given by equation (42), is outlined (after general recommendations provided in ref.[ll]). The presented sample calculations are completed using the results of an experimental study of the structure of Champlain clay and its evolution during a standard one-dimensional consolidation, as reported in ref.[12]. Let us consider a cross-section of soil sample in a plane n(u) with unit normal u, with respect to an orthogonal Cartesian coordinate system fixed in the material. Place a bunch of
98 sampling lines in the plane n(u), having directions v and emanating from a common point P, being the geometric center of the cross-section. Now, consider the 'zero' moment and the 'second' skew-symmetric moment of quantity L(v) in the plane rl(u) with regard to a pole P, i.e., M(u) = I
L(v)ds(v)
(43)
v.v. L(v) de(v)
(44)
J Clu)
Mij(u) = J
C(u)
t
j
Here, C(u) is the unit circle rounding the vector u perpendicularly and de(v) is the arc length measure of C(u) normalized with respect to 2n, as sketch schematically in Fig.1. Substituting the approximations given by equations (42) and (22) in equations (43) and (44) one can express the moments of L(v) in terms offlij(cf.ref.[10]),
(1
M(u) = 2nn ° 1 - - 2 ~ij v i vj
)
(45)
1 Mij(u) = nno[ 12~ij - 21flk(i 'vj)Vk + ( 1 - ~ flkl Vk Vl) vi vj]
(46}
where round brackets ( ) designate the symmetrization of indices. Assume that the experimental observations are made on three perpendicular planes with unit normals u parallel to the base vectors el = {1, 0, 0}, e2 = {0, 1, 0}, e 3 = {0, 0, 1}. In this case, the components of tensor flij are given by the following expressions, M(e 1) f/ll = 2 - ~ ; nn °
M 12(e3) f~12= 2 - ~ Hn° ;
M(e 2) ~22 = 2 - - ; Fin
M(e 3) f$33 = 2 -
o
M23(e 1) fi23 : 2 -
nn o
nn
o
M la(e 2) ;
['/13=2 -
(47)
Hno
According to equations (47) and the definition of 'zero' moment, eq. (43), the diagonal terms of ~ij are proportional to average porosities measured on planes perpendicular to the base vectors, whereas the off-diagonal terms describe the deviation of principal directions of Qij relative to reference coordinate system. The integrals M(ek) and Mij(e k) can be calculated using the trapezoidal integration scheme N-1
M(e k) = 2n ~
L(a m)/N
(48)
m=0 N-1
Mij(e k) = n ' ~ m=0
L(a m) sin(2nm/N)/N
(49)
99 Here, L(am) is the fraction of the test line (occupied by voids), which orientation v is described in terms of an angle a, i.e.am = mn/N (m = 0, 2,...,N-l). The above estimates for the components of f~ij,in terms of the moments of L(v) were employed to evaluate the porosity distribution and its evolution induced by Ko-consolidation process. In particular, the experimental data reported in ref.[12] were used. According to ref.[12], the undisturbed samples of the Champlain clay were consolidated in a standard oedometer with a 1.5 load increment ratio applied every 24 h to 23 kPa, 124 kPa, 421 kPa and 1452 kPa, respectively. Load was then released in several stages and samples were left24 h under the load of 4 kPa. Specimens for microscopic structure identification were trimmed in a form of small sticks with a thin wire, quick frozen in Freon 22, cooled by liquid nitrogen, freeze-fractured and sublimated. The observations of the clay cross-sections were performed using J E O L 25 scanning electron microscope. A graphical representations of porosity, in a horizontal and a vertical plane, at two different levels of the consolidation process, are presented in Fig.2, after ref. [12]. For the purpose of lineal analysis, a discus p a r t of the cross-sections presented in Fig.2 was considered.
The set of eighteen probe lines was drawn (N = 18), with o r i e n t a t i o n s
am = (m-1)n/N (m -- 1, ...,18),as sketched in Fig.3. Then, the total length of the interceptions of each test line with pores was measured and, subsequently, the length of the interceptions per unit length of the test line L(am) was evaluated. The results of measurements and calculations are summarized in Table 1 (intact sample) and in Table 2 (the sample consolidated under the pressure 421 kPa).
Figure 1:
U n i t vector v of a circle and unit normal vector u to the plane of the circle.
1 O0
vertical plane
1 I nt ac t
421 kPa horizc~t~l plane
Intact Figure 2:
421 kPa Graphical representation of porosity at different stages of consolidation of Champlain c|ay (after ref. [12]); voidsare black and particles are white• 1
2
~
3
7
11
i2 13
Figure 3:
A discus part of the cross-sectionfrom Fig. 2, subtracted for the purpose of lineal analysis.
101 T a b l e 1. C a l c u l a t i o n s o f p o r o s i t y d i s t r i b u t i o n f o r t h e i n t a c t s a m p l e
Horizontal Plane
Vertical Plane
Direction C[m
1 (am)
L (am)
E 1 (a m)
L
(am)
(1)
0°
95.5 m m
0.530
48.0 m m
0.267
(2)
I0 °
109.5 m m
0.608
42.5 m m
0.236
(3)
20°
99.0 m m
0.550
24.0 m m
0.133
(4)
30 °
122.5 m m
0.680
29.0 m m
0.161
(5)
40 °
92.0 m m
0.511
20.0 m m
0.111
(6)
50 °
97.0 m m
0.539
32.0 m m
0.178
(7)
60 °
100.5 m m
0.558
47.0 m m
0.261
(8)
70 °
92.0 m m
0.511
30.0 m m
0.167
(9)
80 °
108.0 m m
0.600
27.0 m m
0.150
(10)
90 °
111.0 m m
0.616
13.5 m m
0.0750
(11)
100 °
110.5 m m
0.614
22.0 m m
0.122
(12)
110 °
110.0 m m
0.611
50.5 m m
0.280
(13)
120 °
86.0 m m
0.478
47.5 r a m
0.264
(14)
130 °
86.0 m m
0.478
58.0 m m
O.322
(15)
140 °
96.0 m m
0.533
31.0 m m
0.172
(16)
150 °
67.5 m m
0.375
26.5 m m
0.147
(17)
160 °
67.0 m m
0.372
21.0 m m
0.117
(18)
170 °
101.5 m m
0.564
53.0 m m
0.294
M (e3)
3.39
MI2 (e 3) = 0.0240
M (e l ) = 1.21
M23 (e l ) = - 0 . 0 0 2 5 8
102
Table 2. Calculations of porosity distribution for the consolidated sample (421 kPa) Horizontal Plane
Vertical Plane
Direction •m Z 1 (am)
L (am)
Z 1 (am)
L (am)
(1)
0°
47.2 mm
O,262
50.0 mm
0.278
(2)
10 °
55.0 ram
0 305
33.0 mm
0.183
(3)
20 °
50.0 mm
0.228
43,5 mm
0,242
(4)
30 °
58.5 mm
0.325
36.0 mm
0,200
(5)
40 °
33.0 mm
0.183
28.5 mm
0,158
(6)
50 °
40.0 mm
0.222
30.0 mm
0.167
(7)
60 °
55.5 mm
0.308
22.0 mm
0.122
(8)
70 °
94.0 mm
0.522
15.0 mm
0.0833
(9)
80 °
53.0 mm
0.294
12.5 mm
0.0694
(I0)
90 °
33.0 mm
0.183
4.0 mm
0.0222
(11)
100 °
0.0 mm
0.0
7.0 mm
0.0389
(12)
110 °
12.0 mm
0.0663
5.0 mm
0.278
(13)
120 °
14.0 mm
0.0777
29.0 mm
0.161
(14)
130 °
29.5 mm
0.164
23.0 mm
0.128
(15)
140 °
27.0 mm
0.150
23.5 mm
0.130
(16)
150 °
57.5 mm
0.319
15.5 mm
0.0861
(17)
160 °
35.0 mm
0.194
30.0 mm
0.167
(18)
170 °
51.5 mm
0.286
29.5 mm
0.164
M (e3) -- 1 . 4 5
M (el) = 0.847
M12 (e3) = 0 . 0 4 1 6
M23 ( e l ) = - 0 . 0 0 2 2 4
103 During the one dimensional consolidation process, both u n d e r in-situ conditions and in oedometer, a n axisymmetric deformation mode is imposed. Therefore, the evaluation of the porosity distribution was restricted to one vertical plane only. In calculating the components of tensor flij it was assumed that, M(e 2) = M(e 1)
and
M13(e2) =
M23(el)
(50)
The results, using expressions (47) and values given in Table I a n d Table 2, are presented below.
Intact sample: n o -- 0.308,
eo = 0.445
~11 -- 0.745,
~22 = 0.745,
~33 = -1.508
fl12 = 0.0,
fl23 = -0.00532,
fl13 = -0.00532
Consolidated sample: no = 0.167,
eo = 0.200
•]1 = 0.381,
fl22 = 0.381,
fl33 = - 0 . 7 6 2
~12 = 0.0,
~23 = 0.00858,
~13 = 0.00858
| .0
o
int~t ~ m m a eonsollda~ specimen
.......
!
O.5
0
0.5
l,O
I
i
1.0
Figure 4:
I
i
i
I
|
0.5 Radial
I
I
I
0.0 void
I
i
I
I
0.5 ratio
I
i
i
!
1.0
Spatial d i s t r i b u t i o n of void r a t i o for i n t a c t a n d c o n s o l i d a t e d s a m p l e s of C h a m p l a i n clay.
I O~
The porosity profiles, for both intact and consolidated sample, are presented in graphical form in Fig.4. It should be noted that the average porosities no, as estimated by the described stereological procedure, are in agreement with values reported in ref.[ 12].
NUMERICAL SIMULATIONS The performance of the mathematical framework incorporating the above described measure of material microstructure, has been examined in ref.[9]. In this section, some additional numerical examples are provided involving tests with a continuous rotation of principal stress directions• The main objective here is to illustrate the ability of the proposed formulation to respond to stress trajectories described by eq.(1). The quantitative aspects are, at this stage, rather speculative, as no experimental information is available. All numerical simulations pertain to a normally consolidated (under Ko-conditions) clay. The material is assumed to be initially isotropic and is described by the set of parameters specified in the Appendix. The simulations are based on the representation (34) with the condition of neutral state defined by eq.(36). The results shown in Fig.5 correspond to a specimen which has been consolidated, under Ko-conditions , to an axial stress of 1040 kPa and subsequently subjected to a progressive rotation of principal stress directions, from 0 ° to 180 °. Figures 5a and 5b show the history of the components of total and elastic strain induced by rotation. In the elastic range the response is governed by eq.(5), i.e. the imposed stress trajectory results in generation of reversible deviatoric strain only (Fig. 5b). The plastic flow contributes to both distortional as well as volumetric strain (Fig.5a). It should be mentioned that the only experimental information available, for a similar loading configuration, is that pertaining to sand (ref. [3]). The predicted deformation history seems to be in a qualitative agreement with results reported in [3]. Fig.5 is complemented by Fig.7 which shows the evolution of the spatial distribution of void ratio during the entire loading history. The Ko-consolidation results in an induced crossanisotropy which is displayed through a significant bias in the directional distribution of void ratio. Subsequent history, consisting of rotation of the principal stress directions, results in the deviation of principal axes of tensor ~ij from the coordinate system fixed in the sample. Fig.6 shows a hypothetical response of the same Ko-consolidated sample, subjected to an undrained uniaxial compression combined with a continous rotation of the principal stress axes. It is evident, from Fig.6a, that the rate of rotation influences the effective stress paths and thus, the predicted undrained shear strength. The deviatoric characteristics, as shown in Fig.6b, indicate a progressive reduction in ultimate strength with increasing amount of rotation. Finally, Fig.8 shows the evolution of the spatial distribution of void ratio during the considered loading programme•
105 "-
3
(a) t,otml mr,a-rain 0
2
e--Ill- IE • AEI%
o -1
-30
15 30 45 60 75 90 i 0 5 120 135 150 165 180 Rotation of principal stress axes (degree)
3 ~,
0 VJ
(b) elaJtJe 8b"tJn o--o--O Et -m-mE2 --,o,-- A E02
2
aJ
o
- ~--~-~..~..e
O e...) -3
0
ntliLlalloilltlollt 15 3 0 45
60
~ ~ e ~ - -
75
Rotation of principal Figure 5:
90
Jl n tlialiil' nl'J i 0 5 120 135 150 165 l O 0
stress
axes
(degree)
Strain history due to pure rotation of the principal stress system following Koconsolidation; (a) total strain, (b) elastic strain.
O
(~I~l)
\
lttt
,art a¢
o!.zo:~.~(I
~++i
~tt, _~_/~ s s ~ s
I
0
0 0
II
W
o o
&,;
o o
I::4 i I I
0 o
o
c
o
o
oo
=
C o
2
¢=k
B
107 1.00
o
0.50
0.00 t~
0.50
1.00 1.00
Figure 7:
0.50 0.00 Radial void
0.50 ratio
1.00
Evolution of spatial distribution of void ratio for loading history considered in Fig. 5.
l.O
rotatiom o
0.5
o
0.0
°~
F-4 °~
0.5
i.O
i
I.O
Figure 8:
i
i
i
|
I
i
i
0.5 Radial
I
|
0.0 void
i
i
!
]
!
0.5 ratio
l
I
I
1.0
Evolution of spatial distribution of void ratio for loading history shown in Fig. 6 (the rate ofrotationofprincipal stress axes: 0.2°/kPa).
108 F I N A L REMARKS The discussion presented in this paper focused on the sensitivity of plastic flow to the change in orientation of the principal stress system. It has been d e m o n s t r a t e d t h a t constitutive concepts built within the classical framework of theory of plasticity do not respond to trajectories described by eq.(1). Indeed, these trajectories result in a neutral process which does not generate any irreversible deformations. A mathematical framework has been introduced which employs a continuum measure of internal fabric. The proposed approach results in a phenomenological law with a 'structure' which can describe, at least in a qualitative manner, various manifestations of clay anisotropy [9]. In addition, the mathematical form of the governing equations permits, in general, the modelling of sensitivity of soil response to the rotation of principal stress axes.
REFERENCES Miura, K., Miura, S. and Shosuke T., Deformation behaviour of anisotropic dense sand under principal axes rotation, Soils and Foundations, 26 (1986), 36-52. Matsuoka, H., Suzuki, Y. and Murata, T., A constitutive model for soils accounting for principal stress r o t a t i o n and its application to finite element analysis, In: N u m e r i c a l Models in Geomechanics, Ed. by S. Pietruszczak and G.N. Pande, Elsevier Applied Science Publ., London (1989), 155-162. Ishihara, K. and Towhata, I. Sand response to cyclic rotation of principal stress direction as induced by wave loads, Soils and Foundations, 23 (1983), 11-26. 4.
Saada, A. and Bianchini, G. (Editors)
Constitutive equations for granular non-cohesive soils, Balkema Publ., Rotterdam (1988). 5.
Dafalias, Y.F. and Herrmann, L.R., Bounding surface formulation of soil plasticity, In: Soil M e c h a n i c s ; Transient and Cyclic Loads, Ed. by G.N. Pande and O.C. Zienkiewicz, John Wiley & Sons Inc., New York (1982).
6.
Pietruszczak S. and Mroz Z., On hardening anisotropy of K -consolidated clays, Int. J o u r n . Num. Anal. Meth. Geomeeh., 7 (1983), 19-38.
7.
Pande G.N. and Sharma, K.G., Multilaminate model of clays; a numerical evaluation of the influence of rotation of principal stress axes, Int. J o u r n . Num. Anal. Meth. Geomeeh., 7 (1983), 397-418.
109 8.
Pietruszczak S. and Pande, G.N., Multilaminate framework of soilmodels; plasticityformulation, Int. Journ. Num. Anal. Meth. Geomech., 11 (1987), 651-658.
9.
Pietruszczak, S. and Krucinski, S., Description of anisotropic response of clays using a tensorial measure of structural disorder, M e c h a n i c s o f M a t e r i a l s , (1989), in press.
10.
Kanatani, Ken-lchi, Distribution of directional data and fabric tensor, Int. Journ. Eng. Sci.,22 (1984), 149-164.
11.
Kanatani, Ken-lchi, Procedures for stereological estimation of structural anisotropy, Int. Journ. Eng. Sci.,23 (1985), 587-598.
12.
Delage, P. and Lefebrve, G., Study of the structure of a sensitive Champlain clay and its evolution during consolidation, Can. Geotechn. Journ., 21 (1984), 21-35.
A P P E N D I X : Selection of material functions and parameters The numerical simulations presented in this paper were completed using representation (34) in which the isotropicyield surface f = f(oijeiip) = 0 was assumed to consist of: (i)
an irregular paraboloid (forI ~ Ix): _ q g ( 0 ) ( I o - I x) f' = o + [(i_ix)2_i.]gx = 0 ,2x
(ii)
a n ellipsoidal cap (for I ~ Ix): f , = ~2 + rl2g2(e) [ ( i _ i ) 2 _ ( i o _ i x ) 2 1 = 0
T h e s t r e s s i n v a r i a n t s a p p e a r i n g in t h e s e e q u a t i o n s a r e defined as
I = --
O=
SijSij
;
0-----3 sin-I
where J3 = ~sijSjk Ski. The only material function employed above is that of Io = Io(eP). This function has been chosen, according to Critical State formulation, as / \ p e-- e p • Io = I ~ ° ' e x p ~ ~ _ ' ~ - ) ;
e P = (l+e)eiP i
w h e r e Jt a n d ~ a r e m a t e r i a l c o n s t a n t s identified from h y d r o s t a t i c c o m p r e s s i o n test.
110 The requirement of a continous coalescence of both segments of the yield surface results in the following equation for the position of the merging point q Io Ix = ( q + qf)'
6sin ~b qf
3-sinqb
where q = const, and ¢ is the angle of internal frietion. The numerical simulations carried out pertain to 0 = n]6, i.e. g(0)= 1. The material parameters defining isotropie response correspond to Weald clay and have been selected as: v = 0.27, X = 0.091, ~ = 0.03, ¢ = 26 °, q = 0.594, eo = 0.95 (at initial hydrostatic pressure of 50 kPa). The material functions entering the evolution law, eq.(25), have been chosen as e° - e s a 1 = A;
a 2 = B{1-exp(-In)} ;
I
n
-
eL _ es
where eL and es are void ratios corresponding to liquid and shrinkage limits, respectively. Finally, the function h, eq.(33), has been defined as
where ho and ¥ are constants. The simulations were completed assuming: ec=l.25,
e s = 0.10,
A=0.90,
B=425,
ho= -3.75×105 ,
¥=0.25
A more detailed discussion p e r t a i n i n g to the selection of m a t e r i a l f u n c t i o n s and t h e identification of parameters is provided in ref.[9].
Received 16 October 1989; revised version received 11 December 1989; accepted 11 December 1989
C O N S I D E R A T I O N S ON SOIL R E S P O N S E TO THE ROTATION OF P R I N C I P A l . STRESS DIRECTIONS
S Pietruszczak and S. Krucinski(*) Department of Civil Engineering and Engineering Mechanics McMaster University Hamilton, Ontario, Canada L8S 4L7
ABSTRACT The paper evaluates various plasticityformulations in the context of their ability to describe the influence of the rotation of principal stress axes. The discussion extends to a new theoretical framework which incorporates an implicit measure of material microstructure, derived from the directional distribution of porosity
INTRODUCTION In many practical geotechnical situations, the principal stress directions, in the vicinity of an engineering structure, deviate from those imposed by the depositional conditions. This m a y have a pronounced effecton the resulting displacement field.In recent years a number of researchers (e.g. [1],[2],[3])performed material tests which invariably indicate the sensitivity of soil response to the rotation of principal stress axes. Although the experimental information is stillfragmentary, the evidence gathered so far is very convincing (at least in a qualitative sense) and should not be ignored when formulating appropriate constitutive relations. Over the last two decades a large number of phenomenological models for geomaterials have been proposed [4].These models are mostly cast in the framework of theory of plasticity and they do not, in general, account for the influence of the rotation of principal stress axes. [n order to demonstrate this, let us examine closely the mathematical structure of some of the well established formulations. Consider an arbitrary loading history confined to pure rotation of the principal stress system. In other words, assume that a sample, under initial stress oijo,is subjected to the stress increment Ao..=T. T. o ° - o °. q ip Jq pq q
(1)
(*)On leave of absence from Technical University of Lodz, Poland 89
Computers and Geotechnics 0266-352X/90/S03.50 © 1990 Elsevier Science Publishers Ltd, England. Printed in Great Britain
90 Here, Aoij and o°ij are referred to the same Cartesian coordinate system fixed within the material and Tij is the usual transformation tensor. The stress path (1), when imposed on a sample, results in the deviation of principal stress axes from those corresponding to o°ij, while the principal values of o°ij are preserved.
Suppose at first, that the material is isotropic,
linearly elastic. In this case, the Hooke's law takes the form h~ii= hOii/3K;
heij= Asij/2G
(2)
where sijand eijrepresent the deviatoric parts of the stress and strain tensors and K,G are the elasticconstants (i.e.bulk and shear modulus, respectively). The orthogonality conditions on the directioncosines Tij,i.e.
Tip Tiq = 8pq
(3)
constrain the stress path (1) to satisfy Aoii = 0;
ASij = Aoij
(4)
so that, according to eqs. (2) A~ii = 0;
1 Aeij= ~-~(TipTjq-SipSjq) O°pq
(5)
Thus, the imposed path results in generation of deviatoric strain eij, while the volume change remains zero. It should be noted now, that the existing experimental evidence ([1],[2]} clearly indicates that the stress trajectories satisfying eq. (1) induce not only reversible but also permanent deformations of both distortional as well as volumetric nature. Let us examine now, in a similar context, the formulations derived from the theory of plasticity. The simplest of all, are the descriptions admitting isotropic hardening. In the most common approach, the evolution of the yield surface, f = 0, is linked to the generation of irreversible volumetric strain, {:Pii, so that f(oij, cp) = O;
0
o~o.. U
(6)
In such a ease, the consistency condition (7) aoij
~ii
results in
(8) Since f - 0 is assumed to be an i s o t r o p i c function of oij, then for all histories consisting of pure rotation of principal stress axes, eq. (1),
91 af ao.. u
6
ij
=o
~
~,=o
(9)
i.e.the loading process is a 'neutral'one and does not produce any plasticdeformations. The same conclusion can also be reached when employing multi-yield loci theories of soildeformation. In this case, the associated flow rule reads
--
(10)
n aq..
1] n
l]
where the summation includes only those terms which correspond to f n : 0 and f n : 0 . Let fm ( f, (n= 1, 2, ...., N) be expressed in a form fm(Oij ' ¢'P'n' cP} = 0
(II)
where ~P = (elfeiip)I/2.The consistency condition for any f m =
}m = __~m6.. + ~ ~ aoij q
__ +
_n
n [ ~ii /k;kk
0 takes the form
n_
= 0
(12)
00{ijl °~o{ij)
where afn/ao{ij)designates the deviatoric part of the gradient tensor. Denoting (afmDoij)oij by tkrn and the term contained in square brackets by l.lmn, a set of simultaneous equations (12) can be rewritten in the form
{~n} = [[.Imn]-I {~krn}
(13)
It is obvious that if {/~m}= {0} then, for arbitrary [Pmn], {~n}= {0}. The effects of deformation-induced anisotropy are often described by employing the framework of kinematic hardening. In a typical formulation the yield surface is assumed to undergo translation guided by a tensorial argument aij, i.e. frn(Oij, Uij, ePi) = 0
(14)
where aij=aij(eP~), i.e. the evolution of this surface is implied by the generation of plastic deformations. According to eq. (14) the consistency condition becomes o~f , o~f o~aij 0]"'+ - - - - ~ 1 aOiJ 0aij ~ 1
--
o~/ + --~k
=0
~ii
and results in a similar functional form to that of eq.(8),i.e.
(15)
92 1 ( of
)
~f
Of
of °°ij
of
(16)
In this case again, the pure rotation of principal stress directions describes a 'neutral' process, for which ~ = 0 and thus ~P ---0. Finally, another common description is that of ' b o u n d i n g s u r f a c e ' plasticity [5],[6] Here, for stress histories penetrating the interior of the bounding surface, the deformation process is described in terms of evolution of an inner loading surface
fl(oij, oij) = 0
(17)
The kinematics of this surface (i.e. direction of translation) is guided by a 'conjugate' stress tensor o¢ij located on the bounding surface. A typical form of the evolution law for the components of all is
aij ---- ~lq/ij ;
tl/ij
q/ij(Okl,Okl)
(18)
where l~ is a scalar parameter, determined via consistency condition
f"l =
Ofl ° Ofl dO.'--"~OiJ + 0ct..-- Vii ~ ---- 0 ~j
According to eq.(19), the case of
(19)
~j
(dfl/OOij)oi j -- 0 results in la = 0 and the surface fl = 0 remains
stationary. Therefore, once again, the formulation does not respond to the rotation of principal stress directions. It is evident from the above considerations, that none of the existing frameworks, based on classical concepts of isotropic/kinematic hardening, is capable of taking into account the sensitivity of plastic flow to the rotation of principal stress axes. Indeed, if the evolution of the yield surface (regardless of its form) results exclusively from the generation of plastic deformations, then these effects can never be properly embraced. There are several ways of enhancing the plasticity framework so that it responds to a change in orientation of the principal stress system. This has been accomplished, for example, in the formulation proposed in ref.[7] (see also ref.[8]), whereby the macroscopic behaviour of a sample has been deduced from deformation pattern along an infinite number of randomly oriented planes, by using an appropriate averaging procedure. An a l t e r n a t i v e approach, initiated in ref.[9], may be derived by enriching the set of classical functions defining the state of material by new tensorial function(s) which is/are capable of representing the soil fabric and its evolution. In such a case, the resulting mathematical structure remains phenomenological; the material response however, is an implicit function of the evolution of the fabric.
93 In the next section, the framework outlined in ref.[9] is reviewed in the context of its ability to describe the influence of the rotation of principal stress axes. Subsequently, the details concerning the i d e n t i f i c a t i o n of the t e n s o r i a l m e a s u r e of i n t e r n a l s t r u c t u r e a r e provided, followed by some numerical simulations p e r t a i n i n g to drained and undrained tests involving a continous rotation of principal stress directions.
PLASTICITY FORMULATION FOR ANISOTROPIC SOIL In ref.[9] an a t t e m p t was made to describe the anisotropy of soil structure in terms of a n i s o t r o p y in its p h a s e d i s t r i b u t i o n . As a n i m p l i c i t m e a s u r e of i n t e r n a l s t r u c t u r e , a 'directional porosity' has been employed, which could be identified as a generalized, direction dependent counterpart of 'porosity', i.e. the scalar valued q u a n t i t y defining the void space fraction. In order to construct a n average, non-singular measure of voids distribution, let us isolate in the vicinity of a m a t e r i a l point (P) a unit sphere (S), which encloses a representative volume of the material (V). Consider now a test line of length [ , = 2R with the orientation v with respect to a C a r t e s i a n coordinate system fixed at P. Assume t h a t l(v) represents the total length of interceptions of this line with soil pores. Thus, the fraction of L occupied by voids, L(v), can be defined as
L(v) = l(v)~ ;
l(v) = ~
Ii (v)
(20)
The m e a n value of the quantity L (v) averaged over the domain S, is 1 Lay = 4n
f
S
L(v) f(v)dS ;
II
4n
S
f(v)dS = 1
(21)
where tlv) is a scalar valued function describing the spatial distribution of test lines. It can been shown t h a t for uniformly distributed test lines, the first integral in eq.(21) is the measure
of average porosity no in the vicinity of P, whereas the lineal fraction occupied by pores is an unbiased estimator of the volume fraction of voids in the direction v, i.e. n o = Lay ,
n(v) - L(v)
(22)
The scalar valued function n(v), defined over the unit sphere S, can be represented by the generalized double Fourier series. The desired best fit approximation can be established by the 'least square' method leading to the r e p r e s e n t a t i o n in t e r m s of symmetric traceless tensors ~ij, flijkl, ...., cf.[10] n(vi) ~ no(1 + ~ij vi vj + ~ijkl vi vj v k v I + . . . . )
(23)
94 The higher rank tenors ~ijkl ..... relate to higher order fluctuations in void space distribution. Thus, in order to describe a smooth orthogonal anisotropy it is sufficient to employ an approximation based on the first two terms of the expansion (23). Incorporation of the concept of 'directional porosity', eq.(23), into plasticity framework requires the formulation of an appropriate evolution law for the components of fill. The latter should refer the rate of change of f~ij to the deformation history. In general, the components of flij may be related to the strain rate deviator eij, through an isotropic tensor valued function • • ° ~ij = ~ij(~kl ' ekl' no) (24) The experimental evidence suggests that for a certain class of geomaterials (e.g. clays, silts) the principal axes of ~ij may be considered as coaxial with those of the strain rate deviator. The above assumption results in significant simplification of representation (24), namely
~ij ---- ~ eij ;
~ = al(no) + a2(no) flik ~ki
(25)
where a's are scalar valued functions of the average porosity. Based on existing experimental evidence, it is now postulated that the deformation process alone, is not capable of inducing changes in the microstructure which would affect the functional form of the failure criterion F = 0 or result in anisotropy in the elastic range. The latter effects can only be attributed to the inherent anisotropy in the fabric of the material, which manifests itself in the size, shape and geometrical arrangement of the particles. In order to describe this anisotropy the concept of the fabric tensor Oil has been introduced, with aij defined as a tensor-valued function of ~ij, i.e. aij -- aij(flkl). Thus, in materials with large grain size (e.g. gravel with flat grains, sand with elongated grains, etc.),
Oil : Oij(~kl , akl) F = F(oij,
in elastic range
(26)
aij)= 0
whereas in materials with relatively small grain size the fabric remains isotropic, i.e. aij = 8ij. The proper mathematical representation of both the elasticity tensor and the failure criterion in terms of fabric tensor is again provided in ref.[9]. The anisotropy in the spatial distribution of porosity, as described by flij 2 0 , will
inevitably affect the plastic flow. The procedure for incorporation of the tensor ~ij into plasticity framework is not unique. One way of formulating the problem is to assume that f/ij affects directly the functional form of the yield criterion. Then, Of ~P" = ~ - ~J clo.. Ij •
f = f(°ij' ~'~ij K); '
(27)
95
i.e.,f= 0 depends on the deformation history (recorded by K) and ten (in general) functionally independent invariants of both tensors oij and f/ij.Now, it is interesting to note that the yield surface, defined by eq.(27), will undergo progressive evolution not only during an active loading process, but for all histories (including those for which f<0) associated with the rebuilding of material microstructure. Consequently, such a formulation will account for the influence of the rotation of principal stress axes. In order to demonstrate this, let us assume again z = ePi,and write the consistency condition in the form =
~Oij
+
dOij
-dfl.. -~ ij + - ~# de~ii
~k =
0
(28)
Utilizing the evolution law (25) for the components of flij, i.e. fiij = l]eiej %13~'p" = ~[3 ~ +l]~t Of 'J
2 G ij
(29)
o~(ij}
and substituting the above representation back in eq.(28),one obtains =
l(~f
+2Gt3 0t" )
H =
'
af
of
~
0f
(30)
It is evident from eq.(30) that for (df/aoij)6ij = 0 gP = Cijkl °kl ;
13 of of Cijkl -- 2GH doij aflkl
(31)
which implies that for histories experiencing pure rotation of principal stress axes, the progressive generation of plastic strain takes place, provided that (of/af~ij)fij > 0. A simpler description than that based on eq.(23) may be obtained by resolving the plastic strain rate into components governing a pseudo-isotropic response and a diversion from it. The flow rule can then be expressed as =
ao.. + 1~Gij ," q
f = f(oij , K) = 0
(32)
where the first term represents the strain rates generated in a fictitious isotropic medium, whereas Gij is a symmetric second order tensor whose components are function of f/ij and the deformation history, i.e. Gij = Gij(flkl, K). In the first approximation, one may simply take Gij = h f~ij;
~ = ~t
(33)
where h is a scalar valued function of the deformation history. Substituting eq.(33) into eq.(32) yields
96
•
=
+h~ij
- -
(34)
t~
The above flow rule is analogous to a non-associated law, i.e. according to eq.(34), the induced/inherent anisotropy results in a progressive deviation of the direction of plastic flow from that specified by the gradient tensor Although the formulation (34) is very attractive, in the sense of its numerical simplicity and performance, it does not respond to the rotation of principal stress axes. Indeed, when K is identified with ~P ~, i.e. plastic volumetric strain generated in a fictitious isotropic medium, then the consistency condition (7) results in a similar functional form to that of eq.(3) H\0o.. ij/; tj
H - -
a ~-~.. C~kk 11
(35)
It appears that the only rational way of enhancing the representation (34}, is to define the condition of the neutral state relative to the actual direction of plastic flow, i.e. --
- -
a(i.. l]
+
h
(36)
~ij
In such a case, it is evident that the stress trajectories resulting from pure rotation of principal stress axes will generate the plastic strain rates
.
q = Cijkl Okl ;
)
Cijkl =
+ h ~ij ~kl
(37)
U
provided h flij (~ij > 0. Finally, an alternative mathematical formulation may be obtained by admitting a coupling between both terms in the flow equation (32). Such a flow rule is then analogous to that derived in multi-yield loci theories, implying that the direction of plastic flow is sensitive to the direction of stress rate. Assume, for this purpose, that the components of Gij, eq.(33), are derived from a potential function Q = Q(~ij) = const. ;
~ij -- no(~ij + 8ij)
(38)
where ~ij is defined in such a manner that the equation (23) becomes n(v i) ~ ~ij vi vj. Let
Gij = h C~jkl 0~kl;
l~: O~ij q
(39)
97
where C"jkl is the elastic compliance tensor and fl= fl(~ij)=const, represents a family of loading surfaces defining the neutral states for structural rearrangement. Then, according to eq.(32),for la > 0 ~p •
=
~ - af JchC e ( aQ ~Ofl ) ~ doij ijkl 8~kl a~pcl Pq
(40)
It should be noted that in the above flow rule, the ptastie strain rates e'ij = ~ Gii, eq.(32), can be interpreted as being derived from residual stress rates Sail,
o' __ C e • R Cij ijkl Okl ;
* R
Okl =
h
- -
--
O~ Pq] O~ pq kl
(41)
so that such an approach is analogous to strain-space plasticity formulation. By inspecting eq.(40), it is again obvious that this framework is eapoble of responding to pure rotation of the principal stress system.
Under (df/doij)oij = 0, one obtains the
contribution from the second term in eq.(40), which does not vanish provided that l~ > 0. In conclusion, the reviewed constitutive relations are capable of describing the influence of the inherent as well as induced anisotropy on the deformation process. Their mathematical structure permits, in general, the modelling of the sensitivity of soil response to the rotation of principal stress axes. ON S T E R E O L O G I C A L E S T I M A T I O N OF P O R O S I T Y D I S T R I B U T I O N If the approximation of voids distribution, eq.(23), is confined to the representation n(v) -= no(1 + f~ij vi vj)
(42)
i.e. spherical harmonies of degree higher than two are disregarded, and if, according to equation (22) n(v) = L(v) then, the components of tensor f~ij can be estimated from the observation of L(v) on a limited number of cross-sections. In what follows, the procedure for the stereological estimation of the directional porosity distribution, given by equation (42), is outlined (after general recommendations provided in ref.[ll]). The presented sample calculations are completed using the results of an experimental study of the structure of Champlain clay and its evolution during a standard one-dimensional consolidation, as reported in ref.[12]. Let us consider a cross-section of soil sample in a plane n(u) with unit normal u, with respect to an orthogonal Cartesian coordinate system fixed in the material. Place a bunch of
98 sampling lines in the plane n(u), having directions v and emanating from a common point P, being the geometric center of the cross-section. Now, consider the 'zero' moment and the 'second' skew-symmetric moment of quantity L(v) in the plane rl(u) with regard to a pole P, i.e., M(u) = I
L(v)ds(v)
(43)
v.v. L(v) de(v)
(44)
J Clu)
Mij(u) = J
C(u)
t
j
Here, C(u) is the unit circle rounding the vector u perpendicularly and de(v) is the arc length measure of C(u) normalized with respect to 2n, as sketch schematically in Fig.1. Substituting the approximations given by equations (42) and (22) in equations (43) and (44) one can express the moments of L(v) in terms offlij(cf.ref.[10]),
(1
M(u) = 2nn ° 1 - - 2 ~ij v i vj
)
(45)
1 Mij(u) = nno[ 12~ij - 21flk(i 'vj)Vk + ( 1 - ~ flkl Vk Vl) vi vj]
(46}
where round brackets ( ) designate the symmetrization of indices. Assume that the experimental observations are made on three perpendicular planes with unit normals u parallel to the base vectors el = {1, 0, 0}, e2 = {0, 1, 0}, e 3 = {0, 0, 1}. In this case, the components of tensor flij are given by the following expressions, M(e 1) f/ll = 2 - ~ ; nn °
M 12(e3) f~12= 2 - ~ Hn° ;
M(e 2) ~22 = 2 - - ; Fin
M(e 3) f$33 = 2 -
o
M23(e 1) fi23 : 2 -
nn o
nn
o
M la(e 2) ;
['/13=2 -
(47)
Hno
According to equations (47) and the definition of 'zero' moment, eq. (43), the diagonal terms of ~ij are proportional to average porosities measured on planes perpendicular to the base vectors, whereas the off-diagonal terms describe the deviation of principal directions of Qij relative to reference coordinate system. The integrals M(ek) and Mij(e k) can be calculated using the trapezoidal integration scheme N-1
M(e k) = 2n ~
L(a m)/N
(48)
m=0 N-1
Mij(e k) = n ' ~ m=0
L(a m) sin(2nm/N)/N
(49)
99 Here, L(am) is the fraction of the test line (occupied by voids), which orientation v is described in terms of an angle a, i.e.am = mn/N (m = 0, 2,...,N-l). The above estimates for the components of f~ij,in terms of the moments of L(v) were employed to evaluate the porosity distribution and its evolution induced by Ko-consolidation process. In particular, the experimental data reported in ref.[12] were used. According to ref.[12], the undisturbed samples of the Champlain clay were consolidated in a standard oedometer with a 1.5 load increment ratio applied every 24 h to 23 kPa, 124 kPa, 421 kPa and 1452 kPa, respectively. Load was then released in several stages and samples were left24 h under the load of 4 kPa. Specimens for microscopic structure identification were trimmed in a form of small sticks with a thin wire, quick frozen in Freon 22, cooled by liquid nitrogen, freeze-fractured and sublimated. The observations of the clay cross-sections were performed using J E O L 25 scanning electron microscope. A graphical representations of porosity, in a horizontal and a vertical plane, at two different levels of the consolidation process, are presented in Fig.2, after ref. [12]. For the purpose of lineal analysis, a discus p a r t of the cross-sections presented in Fig.2 was considered.
The set of eighteen probe lines was drawn (N = 18), with o r i e n t a t i o n s
am = (m-1)n/N (m -- 1, ...,18),as sketched in Fig.3. Then, the total length of the interceptions of each test line with pores was measured and, subsequently, the length of the interceptions per unit length of the test line L(am) was evaluated. The results of measurements and calculations are summarized in Table 1 (intact sample) and in Table 2 (the sample consolidated under the pressure 421 kPa).
Figure 1:
U n i t vector v of a circle and unit normal vector u to the plane of the circle.
1 O0
vertical plane
1 I nt ac t
421 kPa horizc~t~l plane
Intact Figure 2:
421 kPa Graphical representation of porosity at different stages of consolidation of Champlain c|ay (after ref. [12]); voidsare black and particles are white• 1
2
~
3
7
11
i2 13
Figure 3:
A discus part of the cross-sectionfrom Fig. 2, subtracted for the purpose of lineal analysis.
101 T a b l e 1. C a l c u l a t i o n s o f p o r o s i t y d i s t r i b u t i o n f o r t h e i n t a c t s a m p l e
Horizontal Plane
Vertical Plane
Direction C[m
1 (am)
L (am)
E 1 (a m)
L
(am)
(1)
0°
95.5 m m
0.530
48.0 m m
0.267
(2)
I0 °
109.5 m m
0.608
42.5 m m
0.236
(3)
20°
99.0 m m
0.550
24.0 m m
0.133
(4)
30 °
122.5 m m
0.680
29.0 m m
0.161
(5)
40 °
92.0 m m
0.511
20.0 m m
0.111
(6)
50 °
97.0 m m
0.539
32.0 m m
0.178
(7)
60 °
100.5 m m
0.558
47.0 m m
0.261
(8)
70 °
92.0 m m
0.511
30.0 m m
0.167
(9)
80 °
108.0 m m
0.600
27.0 m m
0.150
(10)
90 °
111.0 m m
0.616
13.5 m m
0.0750
(11)
100 °
110.5 m m
0.614
22.0 m m
0.122
(12)
110 °
110.0 m m
0.611
50.5 m m
0.280
(13)
120 °
86.0 m m
0.478
47.5 r a m
0.264
(14)
130 °
86.0 m m
0.478
58.0 m m
O.322
(15)
140 °
96.0 m m
0.533
31.0 m m
0.172
(16)
150 °
67.5 m m
0.375
26.5 m m
0.147
(17)
160 °
67.0 m m
0.372
21.0 m m
0.117
(18)
170 °
101.5 m m
0.564
53.0 m m
0.294
M (e3)
3.39
MI2 (e 3) = 0.0240
M (e l ) = 1.21
M23 (e l ) = - 0 . 0 0 2 5 8
102
Table 2. Calculations of porosity distribution for the consolidated sample (421 kPa) Horizontal Plane
Vertical Plane
Direction •m Z 1 (am)
L (am)
Z 1 (am)
L (am)
(1)
0°
47.2 mm
O,262
50.0 mm
0.278
(2)
10 °
55.0 ram
0 305
33.0 mm
0.183
(3)
20 °
50.0 mm
0.228
43,5 mm
0,242
(4)
30 °
58.5 mm
0.325
36.0 mm
0,200
(5)
40 °
33.0 mm
0.183
28.5 mm
0,158
(6)
50 °
40.0 mm
0.222
30.0 mm
0.167
(7)
60 °
55.5 mm
0.308
22.0 mm
0.122
(8)
70 °
94.0 mm
0.522
15.0 mm
0.0833
(9)
80 °
53.0 mm
0.294
12.5 mm
0.0694
(I0)
90 °
33.0 mm
0.183
4.0 mm
0.0222
(11)
100 °
0.0 mm
0.0
7.0 mm
0.0389
(12)
110 °
12.0 mm
0.0663
5.0 mm
0.278
(13)
120 °
14.0 mm
0.0777
29.0 mm
0.161
(14)
130 °
29.5 mm
0.164
23.0 mm
0.128
(15)
140 °
27.0 mm
0.150
23.5 mm
0.130
(16)
150 °
57.5 mm
0.319
15.5 mm
0.0861
(17)
160 °
35.0 mm
0.194
30.0 mm
0.167
(18)
170 °
51.5 mm
0.286
29.5 mm
0.164
M (e3) -- 1 . 4 5
M (el) = 0.847
M12 (e3) = 0 . 0 4 1 6
M23 ( e l ) = - 0 . 0 0 2 2 4
103 During the one dimensional consolidation process, both u n d e r in-situ conditions and in oedometer, a n axisymmetric deformation mode is imposed. Therefore, the evaluation of the porosity distribution was restricted to one vertical plane only. In calculating the components of tensor flij it was assumed that, M(e 2) = M(e 1)
and
M13(e2) =
M23(el)
(50)
The results, using expressions (47) and values given in Table I a n d Table 2, are presented below.
Intact sample: n o -- 0.308,
eo = 0.445
~11 -- 0.745,
~22 = 0.745,
~33 = -1.508
fl12 = 0.0,
fl23 = -0.00532,
fl13 = -0.00532
Consolidated sample: no = 0.167,
eo = 0.200
•]1 = 0.381,
fl22 = 0.381,
fl33 = - 0 . 7 6 2
~12 = 0.0,
~23 = 0.00858,
~13 = 0.00858
| .0
o
int~t ~ m m a eonsollda~ specimen
.......
!
O.5
0
0.5
l,O
I
i
1.0
Figure 4:
I
i
i
I
|
0.5 Radial
I
I
I
0.0 void
I
i
I
I
0.5 ratio
I
i
i
!
1.0
Spatial d i s t r i b u t i o n of void r a t i o for i n t a c t a n d c o n s o l i d a t e d s a m p l e s of C h a m p l a i n clay.
I O~
The porosity profiles, for both intact and consolidated sample, are presented in graphical form in Fig.4. It should be noted that the average porosities no, as estimated by the described stereological procedure, are in agreement with values reported in ref.[ 12].
NUMERICAL SIMULATIONS The performance of the mathematical framework incorporating the above described measure of material microstructure, has been examined in ref.[9]. In this section, some additional numerical examples are provided involving tests with a continuous rotation of principal stress directions• The main objective here is to illustrate the ability of the proposed formulation to respond to stress trajectories described by eq.(1). The quantitative aspects are, at this stage, rather speculative, as no experimental information is available. All numerical simulations pertain to a normally consolidated (under Ko-conditions) clay. The material is assumed to be initially isotropic and is described by the set of parameters specified in the Appendix. The simulations are based on the representation (34) with the condition of neutral state defined by eq.(36). The results shown in Fig.5 correspond to a specimen which has been consolidated, under Ko-conditions , to an axial stress of 1040 kPa and subsequently subjected to a progressive rotation of principal stress directions, from 0 ° to 180 °. Figures 5a and 5b show the history of the components of total and elastic strain induced by rotation. In the elastic range the response is governed by eq.(5), i.e. the imposed stress trajectory results in generation of reversible deviatoric strain only (Fig. 5b). The plastic flow contributes to both distortional as well as volumetric strain (Fig.5a). It should be mentioned that the only experimental information available, for a similar loading configuration, is that pertaining to sand (ref. [3]). The predicted deformation history seems to be in a qualitative agreement with results reported in [3]. Fig.5 is complemented by Fig.7 which shows the evolution of the spatial distribution of void ratio during the entire loading history. The Ko-consolidation results in an induced crossanisotropy which is displayed through a significant bias in the directional distribution of void ratio. Subsequent history, consisting of rotation of the principal stress directions, results in the deviation of principal axes of tensor ~ij from the coordinate system fixed in the sample. Fig.6 shows a hypothetical response of the same Ko-consolidated sample, subjected to an undrained uniaxial compression combined with a continous rotation of the principal stress axes. It is evident, from Fig.6a, that the rate of rotation influences the effective stress paths and thus, the predicted undrained shear strength. The deviatoric characteristics, as shown in Fig.6b, indicate a progressive reduction in ultimate strength with increasing amount of rotation. Finally, Fig.8 shows the evolution of the spatial distribution of void ratio during the considered loading programme•
105 "-
3
(a) t,otml mr,a-rain 0
2
e--Ill- IE • AEI%
o -1
-30
15 30 45 60 75 90 i 0 5 120 135 150 165 180 Rotation of principal stress axes (degree)
3 ~,
0 VJ
(b) elaJtJe 8b"tJn o--o--O Et -m-mE2 --,o,-- A E02
2
aJ
o
- ~--~-~..~..e
O e...) -3
0
ntliLlalloilltlollt 15 3 0 45
60
~ ~ e ~ - -
75
Rotation of principal Figure 5:
90
Jl n tlialiil' nl'J i 0 5 120 135 150 165 l O 0
stress
axes
(degree)
Strain history due to pure rotation of the principal stress system following Koconsolidation; (a) total strain, (b) elastic strain.
O
(~I~l)
\
lttt
,art a¢
o!.zo:~.~(I
~++i
~tt, _~_/~ s s ~ s
I
0
0 0
II
W
o o
&,;
o o
I::4 i I I
0 o
o
c
o
o
oo
=
C o
2
¢=k
B
107 1.00
o
0.50
0.00 t~
0.50
1.00 1.00
Figure 7:
0.50 0.00 Radial void
0.50 ratio
1.00
Evolution of spatial distribution of void ratio for loading history considered in Fig. 5.
l.O
rotatiom o
0.5
o
0.0
°~
F-4 °~
0.5
i.O
i
I.O
Figure 8:
i
i
i
|
I
i
i
0.5 Radial
I
|
0.0 void
i
i
!
]
!
0.5 ratio
l
I
I
1.0
Evolution of spatial distribution of void ratio for loading history shown in Fig. 6 (the rate ofrotationofprincipal stress axes: 0.2°/kPa).
108 F I N A L REMARKS The discussion presented in this paper focused on the sensitivity of plastic flow to the change in orientation of the principal stress system. It has been d e m o n s t r a t e d t h a t constitutive concepts built within the classical framework of theory of plasticity do not respond to trajectories described by eq.(1). Indeed, these trajectories result in a neutral process which does not generate any irreversible deformations. A mathematical framework has been introduced which employs a continuum measure of internal fabric. The proposed approach results in a phenomenological law with a 'structure' which can describe, at least in a qualitative manner, various manifestations of clay anisotropy [9]. In addition, the mathematical form of the governing equations permits, in general, the modelling of sensitivity of soil response to the rotation of principal stress axes.
REFERENCES Miura, K., Miura, S. and Shosuke T., Deformation behaviour of anisotropic dense sand under principal axes rotation, Soils and Foundations, 26 (1986), 36-52. Matsuoka, H., Suzuki, Y. and Murata, T., A constitutive model for soils accounting for principal stress r o t a t i o n and its application to finite element analysis, In: N u m e r i c a l Models in Geomechanics, Ed. by S. Pietruszczak and G.N. Pande, Elsevier Applied Science Publ., London (1989), 155-162. Ishihara, K. and Towhata, I. Sand response to cyclic rotation of principal stress direction as induced by wave loads, Soils and Foundations, 23 (1983), 11-26. 4.
Saada, A. and Bianchini, G. (Editors)
Constitutive equations for granular non-cohesive soils, Balkema Publ., Rotterdam (1988). 5.
Dafalias, Y.F. and Herrmann, L.R., Bounding surface formulation of soil plasticity, In: Soil M e c h a n i c s ; Transient and Cyclic Loads, Ed. by G.N. Pande and O.C. Zienkiewicz, John Wiley & Sons Inc., New York (1982).
6.
Pietruszczak S. and Mroz Z., On hardening anisotropy of K -consolidated clays, Int. J o u r n . Num. Anal. Meth. Geomeeh., 7 (1983), 19-38.
7.
Pande G.N. and Sharma, K.G., Multilaminate model of clays; a numerical evaluation of the influence of rotation of principal stress axes, Int. J o u r n . Num. Anal. Meth. Geomeeh., 7 (1983), 397-418.
109 8.
Pietruszczak S. and Pande, G.N., Multilaminate framework of soilmodels; plasticityformulation, Int. Journ. Num. Anal. Meth. Geomech., 11 (1987), 651-658.
9.
Pietruszczak, S. and Krucinski, S., Description of anisotropic response of clays using a tensorial measure of structural disorder, M e c h a n i c s o f M a t e r i a l s , (1989), in press.
10.
Kanatani, Ken-lchi, Distribution of directional data and fabric tensor, Int. Journ. Eng. Sci.,22 (1984), 149-164.
11.
Kanatani, Ken-lchi, Procedures for stereological estimation of structural anisotropy, Int. Journ. Eng. Sci.,23 (1985), 587-598.
12.
Delage, P. and Lefebrve, G., Study of the structure of a sensitive Champlain clay and its evolution during consolidation, Can. Geotechn. Journ., 21 (1984), 21-35.
A P P E N D I X : Selection of material functions and parameters The numerical simulations presented in this paper were completed using representation (34) in which the isotropicyield surface f = f(oijeiip) = 0 was assumed to consist of: (i)
an irregular paraboloid (forI ~ Ix): _ q g ( 0 ) ( I o - I x) f' = o + [(i_ix)2_i.]gx = 0 ,2x
(ii)
a n ellipsoidal cap (for I ~ Ix): f , = ~2 + rl2g2(e) [ ( i _ i ) 2 _ ( i o _ i x ) 2 1 = 0
T h e s t r e s s i n v a r i a n t s a p p e a r i n g in t h e s e e q u a t i o n s a r e defined as
I = --
O=
SijSij
;
0-----3 sin-I
where J3 = ~sijSjk Ski. The only material function employed above is that of Io = Io(eP). This function has been chosen, according to Critical State formulation, as / \ p e-- e p • Io = I ~ ° ' e x p ~ ~ _ ' ~ - ) ;
e P = (l+e)eiP i
w h e r e Jt a n d ~ a r e m a t e r i a l c o n s t a n t s identified from h y d r o s t a t i c c o m p r e s s i o n test.
110 The requirement of a continous coalescence of both segments of the yield surface results in the following equation for the position of the merging point q Io Ix = ( q + qf)'
6sin ~b qf
3-sinqb
where q = const, and ¢ is the angle of internal frietion. The numerical simulations carried out pertain to 0 = n]6, i.e. g(0)= 1. The material parameters defining isotropie response correspond to Weald clay and have been selected as: v = 0.27, X = 0.091, ~ = 0.03, ¢ = 26 °, q = 0.594, eo = 0.95 (at initial hydrostatic pressure of 50 kPa). The material functions entering the evolution law, eq.(25), have been chosen as e° - e s a 1 = A;
a 2 = B{1-exp(-In)} ;
I
n
-
eL _ es
where eL and es are void ratios corresponding to liquid and shrinkage limits, respectively. Finally, the function h, eq.(33), has been defined as
where ho and ¥ are constants. The simulations were completed assuming: ec=l.25,
e s = 0.10,
A=0.90,
B=425,
ho= -3.75×105 ,
¥=0.25
A more detailed discussion p e r t a i n i n g to the selection of m a t e r i a l f u n c t i o n s and t h e identification of parameters is provided in ref.[9].
Received 16 October 1989; revised version received 11 December 1989; accepted 11 December 1989