- Email: [email protected]
A two characteristic model for cohesionless soil and its calibration using optimization
Pergamon
International Journal of Plasticity, Vol. 10, No. 3, pp. 309-326, 1994 Copyright © 1994 Elsevier ScienceLtd Printed in the USA. All rights reserved 0749-6419/94 $6.00 + .00 0749-6419(93)E0006-M
A TWO CHARACTERISTIC MODEL FOR COHESIONLESS SOIL AND ITS CALIBRATION USING OPTIMIZATION
ABDULAZEEZ ABDULRAHEEM, M . OMAR FARUQUE, a n d MUSRARRAr ZAU_~',; The University of Oklahoma, Norman
Abstract--An elasto-plastic constitutive model is developed to characterize the stress-strain behavior of cohesionless soils. The model utilizes the concept of two characteristic states represented by two characteristic state lines in the stress space. The two characteristic state lines represent the state of the soil at failure and the state at which the soil passes from compressive to dilative mode of deformation during shearing. The proposed approach provides a better control in predicting volumetric behavior of cohesionless soils. To increase the accessibility and usefulness of the model and reduce the difficulty of the model calibration process, an optimization procedure is used to evaluate the material parameters associated with the model. The model is verified with respect to a test that is used for finding the parameters, and a number of other tests that are not used for finding the parameters. A good correlation between predictions and experimental data is observed.
!. INTRODUCTION
In the past, plasticity theory has been used extensively to formulate constitutive models for cohesionless soils (DRtrCKER et al., [1957]; DI~L~GGIO~ SASDLER, [1971]; LADE & DU~rCAr~, [1975]; Zms~:mw~cz et aL, [1975]; SASDLER, et aL, [1976]; DAr,~d~tAS& HAe,RUAr~N, [1982]; VAI.ASIS & READ, [1982]; D~NAT~E, [1983]; D~S~ ~ FARUQUE, [1984]; FARUQUE& DESAI, [1985]; DES~ etal., [1986]; FARUQtrE, [1987]; Hmxt, [1987]; Am)t~_JAUWAD et al., [1991]). Many of these models characterize stress-strain response in a reasonable manner. Quite frequently, the volumetric behavior is ignored or less emphasized by these models and sometimes it is not used in the model validation process. As evident from past experimental investigations, induced volume change during sheafing is an important characteristic of cohesionless soil. The nature of volume change is primarily controlled by the relative density of the soil and the confining pressure at which the shearing takes place. Depending upon the relative density and the applied confining pressure, a soil can experience either dilation or compression during shearing (LADE, [1977]; LOUNO, [1980]; GRAsoussI ~ M o ~ r , [1982]; FUKUSrm~A~, TArstror~, [1984]; POOROOSnASn~ PmXRUSZ~ZAK, [1985]; ASDULJAtrWADet aL, [1991]; H ~ I N ,
I1989]). Figure 1 shows the effect of relative density and the confining pressure on axial-volumetric response obtained from (stress controlled) triaxial compression tests. Evidently, a relatively dense soil at a low confining pressure experiences a very little shear-induced soil compression and in most part it dilates. On the other hand, at a very high confining pressure, cohesionless softs usually experience compressive volume change all through shearing up to failure. This behavior seems to be exhibited irrespective of the initial relative density of the soil. 309
A. ABDULRAI-IEEM et al.
310
High Confining Pressure
'
\el_
~.'"
1/ 2 /
Low Confining P,'~ssure
",\///'2
Radial Strain
Axial Strain
(a)
Dilatant
I /
Compressive
. -
"q ' -.
Strain
1', 2' Compressive
(b) Fig. 1. (a) Typical stress-strain response of cohesionless soil. (b) Typical volumetric-axial response of cohesionless soil.
This paper presents a general constitutive model for cohesionless soils based on the framework of plasticity theory. One of the unique features of this model is that it incorporates explicitly the two characteristic states of a cohesionless soil in the formulation. The proposed framework, in general, can be used to characterize shear-induced volume change behavior of these soils in a rational manner. In addition, to make the calibration process of the model easier, an optimization technique is also presented in detail. II. P R O P O S E D CONSTITUTIVE MODEL
Formulation of the proposed constitutive model is based on the nonassociated plasticity theory. Two distinct characteristic state lines are proposed to represent the two dif-
Cohesionless soil model
311
ferent characteristic states of cohesionless soils. The first characteristic state line (CSL-1) refers to the state of cohesionless soil at failure and is expressed as
J'~2o MJ~ ,
(1)
=
where M is a parameter related to the coefficient of internal frictional, which in turn depends on the relative density of the material. The second characteristic state line (CSL-2) refers to the state of cohesionless soil at which the rate of volumetric strain momentarily vanishes as the soil passes from compressive mode of deformation to dilative mode of deformation during shearing. Experimental evidence suggests that as the confining pressure increases, the dilation decreases and the volumetric response becomes predominantly compressive. At high confining pressure, dilation is not encountered at all. This phenomenon of decreasing dilatancy with increasing confining pressure can be modeled effectively by letting CSL-2 get closer to CSL-1 with increasing confining pressure. Consequently, the equation of this line is proposed as = M[ 1 - Nexp( -#J~/P.)] J~.
(2)
The exponential term in eqn (2) becomes negligible at higher values of Jl. As such eqn (2) approaches eqn (1) with increasing Jl. Typical plots of eqns (1) and (2) on the Jl - J~-)-22oplane are shown in Fig. 2. Efforts have been made in the past to model the onset of dilatancy and subsequent dilatant behavior by including the second characteristic state line, which essentially diverges from the first characteristic state line as the confining pressure increases (Gi-L~OUSSi~ Mor,m~, [1982]; FLrKuSrnMA~ TATSUOIOt, [1984];
os ~°"
.......
CSL-1
~-~ ......
CSL-2: LinearfolTn CSL-2: Nonlinear form
....
~
.....
Dilative yield surface
p
~
-'"
yield Surface
I"
s~
..
.oo
A
a
4. 4
o,,x-" .,,~"
4
|
/
~
\ C
X
J~ Axis Fig. 2. Schematic representation of the yield caps and the two CSLs.
312
A . ABDULRAHEEM
et al.
DESAI et al., [1986]; CRISTESCU, [1989]). The form of CSL-2, as suggested by eqn (2), is such that it yields a predominantly compressive response at high confining pressure and both the compressive and dilative responses at low and medium confining pressure ranges. Further, the effect of relative density can also be incorporated by making the parameters M and N functions of relative density. However, since experimental data used in this study are such that they fall approximately in the lower confining pressure range, eqn (2) is simplified to a linear form by assuming/, = 0. The simplified equation is given by (3)
~ 2 D = M[1 - N I J , .
Figure 2 shows the simplified linear form of the CSL-2 along with the nonlinear form as represented by eqn (2). Consistent with the definition of CSL-2, all states of stress below the second characteristic state line yield compressive volumetric strain only. A compressive yield surface of elliptical shape that satisfies the above requirements is proposed for the model and is shown schematically in Fig. 2. The tangent to the yield surface is vertical at the intersection with J~-axis and horizontal at the intersection with CSL-2. The normal to the yield surface at any point between the above two intersections has nonnegative slope. The equation of the compressive yield surface is proposed as
Fe = J~ZD -- \/a2 _ ()/ - - C )
2/~R2 = 0
,
(4)
where R is a material constant associated with the ratio of the axes of the ellipse and C is a response function that denotes the value of J~ at the point of intersection of a generic yield surface and the second characteristic state line. In eqn (4), a represents the axis of the elliptical compression yield surface parallel to the JX~2Daxis and is defined in terms of C as a = M[1 -- N e x p ( - ~ C / P a ) ] C.
(5)
The proposed model utilizes a nonassociated formulation to describe the volumetric response in the compressive regime. A plastic potential function, Qc, is defined in the form Qc =
J'~ZZD-- A ~ ~2 -- (J, - C) 2/R2,
(6)
where ,4 is a material constant. Taking the derivatives of Fc and Qc with respect to Jr, the following relationship can be obtained
(7)
aQ¢ _ `4 OF,,
As evident from eqn (7), the constant .4 signifies deviation from normality and should have a positive value in the range 0 ___,4 _< 1. The dilation yield surface refers to the stress states that cause dilating response at a material point. It is bounded by CSL-1 and CSL-2 at all times and identifies with CSL-1 when the confining pressure is very high. An equation of the dilation yield surface that satisfies the above requirements is proposed in the form Fa = J'~ZD -- a[1 -- (1
-
J 1 / C ) MC/a]
(8)
Cohesionless soil model
313
A schematic representation of a generic dilation yield surface is shown in Fig. 2. Evidently, the tangent to the dilation yield surface at the point of intersection with the second characteristic state line is horizontal, and thereby satisfies the slope compatibility with the compression yield surface, Fc = 0. At Jl = 0, the slope of the dilation yield surface is M, and therefore, tangential to the first characteristic state line at that point. Nonassociated formulation is used here to achieve zero rate of dilatancy when the stress state is on the characteristic state lines (both CSL-1 and CSL-2), and nonzero rate of dilatancy when the state of stress lies in between the characteristic state lines. To achieve this, a dilating potential function, Qd = 0, is defined by the following equations:
oo. =q(j,,s47Vo)af4 ,,Jl OQa 0J2D
-
OF,l 0J2D' - -
(9)
(10)
where q(Jl,J2D) = ct (/3 /3t37. + 7)~+"~ Y (1 - G)t~G v
(11)
M J, -G ( J1, J2D ) = M N exp( - M J l /P,,)Jl "
(12)
in which
In the above equations, o~,/3, and 7 are material constants. The function q(J1,J2D) assumes a value of zero at the boundaries (CSL-1 and CSL-2) and is nonnegative elsewhere. Explicit form of Qd is not required since the derivatives of Qa suffice to obtain the elasto-plastic constitutive matrix. The inclusion of a distinct dilating yield surface along with the nonassociated formulation, in terms of material constants a,/3, and 7, provides a much better control on the dilative response of the material. However, depending on the response of the soil, the three parameters, or,/3, and 7, can be reduced to one (i.e. or). Even an associated formulation can be adopted, if the nature of response dictates such simplifications, in both the compressive and the dilative regimes, in which case the total number of parameters will be reduced by four. It should be noted here that the simplification of the CSL-2 from nonlinear to linear form, as represented by eqn (3), does not affect any of the propositions mentioned above.
III. HARDENING BEHAVIOR
In this work, modeling of soil response under monotonic loading is attempted using the concept of isotropic hardening plasticity. Both the compression yield surface, Fc, and the dilation yield surface, Fa, are allowed to expand in the stress space to account for hardening during the elasto-plastic deformation. It is important to realize that the mechanism of inelastic volume change during hydrostatic compression is quite different from that in shear. As a result, a hardening function, with its constants determined
314
A. ABDULRAHEEMet al.
from hydrostatic compression test data, generally cannot describe elasto-plastic stressstrain and volumetric behavior during shearing in a rational manner. This problem can be resolved by expressing the hardening function in terms of two parameters, ~ and rd, which are defined in terms of the incremental plastic strain tensor de~ and the incremental deviatoric plastic strain tensor deft, as follows:
= f (de p deiPj) '/e
(13)
r d = (d/~,
(14)
~d = f (dei~ dei]) l/z.
(15)
where
In the above equations, the symbol f denotes the history. Note that ~d and rd have zero value throughout hydrostatic compression. Referring to Fig. 2, the compressive yield surface, Fc, originates from the Jl-axis (i.e. point X) and terminates at a point on the second characteristic state line (CSL-2) where Jl assumes a value of C. The dilating yield surface, Fd, on the other hand, starts from the origin 0 and terminates at CSL-2 where J~ = C. Since X = C + Ra, where R is the ratio of major to minor axis of the elliptical cap (Fig. 2), and a is the ordinate corresponding to J1 = C, expansion of both Fc and Fd can be described conveniently by making X a (hardening) function of the parameters ~ and r d as
X = 31~ ~,
(16)
where = ~l[1 + 32r321.
(17)
In eqns (16) and (17), ~ , ~2, ~11, and ~72 are material constants. An initially isotropic material, when subjected to a hydrostatic loading, experiences a spherical deformation only (i.e. rd = 0 and ~ = ~v, where ~v = f 1 / , f 3 d e p = 1 / x ~ p, e~ being the total volumetric plastic strain). The hardening function for this case reduces to
x= ~'=
~(c~/4~).,.
(18)
IV. ELASTO-PLASTIC CONSTITUTIVE RELATIONS
Using the concept of nonassociated plasticity theory, incremental plastic strain tensor, de p, can be expressed as de p = dA O_Q_Q OCTij '
(19)
Cohesionlesssoil model
315
where Q is the plastic potential function and dA is an unknown scalar to be determined from the consistency condition of Prager, D f = O, F being the yield function. Following the standard steps of the theory of plasticity (PRAo~R & HODOE [1951]), the elastoplastic constitutive relation tensor, CvPktcan be written as
aF
0Q
C~.,, aa.o aa,.. C~,.kt Cu.~, =
~F
e
a_QQ _
,
(20)
aapq C~rs c~crrs A(~,rd, a~s) where Cfjk~ is the elastic constitutive tensor, and A(~,ra,aij) is a measure of plastic modulus involving derivatives of the yield function F with respect to ~ and ra, and derivative of the potential function Q with respect to a;j. Equation (20) is valid in the compressive as well as in the dilative regimes of shear deformation provided Q is appropriately defined. In the compressive regime, Q = Q¢ (eqn 6), while in the dilative regime, Q = Qa, defined implicitly by eqns (9)-(12).
V. C A L I B R A T I O N OF T H E M O D E L
The proposed model involves different material parameters, which are broken down into directly measurable or "fixed" parameters and the remaining "free" parameters. Determination of the directly measurable or fixed parameters is quite straightforward, and it can be readily accomplished by using well-established laboratory experiments. Determination of the remaining free parameters, however, can make the model calibration procedure quite difficult. These free parameters are generally established by a trial and error method, with the objective being to obtain the best overall fit to a given set of experimental data (Dar~tAs ~ H,~aL~a,m, [1982]). As a result, the accuracy and efficiency of the calibration process can be strongly dependent on the subjectivity of the user as well as the user's knowledge about the material and expertise with the particular constitutive model. In view of this, a computer-aided procedure is developed by casting the model calibration problem into a form which makes it amenable to solution by conventional mathematical optimization techniques. Because this automated procedure greatly reduces the dependence of calibration success on user expertise, it significantly increases the accessibility and usefulness of this model to the general engineering community. V.1. Simplex optimization technique Simplex approach of function minimization is used here as a tool for the evaluation of the free parameters in an optimum manner (SPEI'ZOLEYet aL, [1962]; FLETCm~R, [1980]; FLETCHER, [1981]; D~.NxrALE, [1983]). The procedure is based on the evolutionary operation (EVOP) technique originally introduced by Box in 1957, in which the solution space is explored by means of geometrical configuration of points rather than a set of directions. A number of modifications to this basic approach have since been suggested, but the version attributed to Nelder & Mead, in which the simplex automatically rescales itself according to local geometry of the objective function, is commonly regarded as most effective (Nm~DER * M ~ , [1965]). The primary merits of the simplex method are related to its extreme simplicity and great generality. The procedure
316
A. ABDULRAHEEMet al.
requires only the ability to evaluate the objective function and, typically, makes no assumptions with respect to the behavior of the function and its derivatives. For this reason, it is particularly useful when the function is nondifferentiable, when the gradients are discontinuous, or when the value of the function depends upon some physical measurements and may therefore be subject to random error. As a result of its simplicity, the strategy can be readily implemented into a workable code. In addition, because the procedure normally employs no complicated numerical operations, it is traditionally immune from rounding o f f errors a n d / o r ill-conditioning. Furthermore, the routine is computationally compact and therefore requires only modest amounts of computer storage. According to SWANN [1972], the direct search methods, of which the simplex approach is most successful, have proved to be quite reliable in practice and only rarely fail to locate at least a local minimum. V.2. Numerical algorithm for simplex technique The simplex search, known as the Box algorithm (SPENDLEY et al., [1962]; RELAK~TIS et al., [1983]), which can handle problems with many constraints, is used here due to the fact that it does not require any derivative of the objective function and constraints. The following steps are intended to explain the application o f the Box algorithm:
Step 1: Determine an initial starting point X ° (where X represents the vector containing the individual values denoted by x~, of the ten parameters which need to be optimized), a reflection coefficient &, and termination parameters i and 6. Generate a set of ( P - 1) random feasible points, where P = 2N, N b e i n g the number of decision variables in the problem. One needs to sample Ntimes to generate one point. If the point X i is found to be infeasible, retract half the distance toward the centroid, X c, of the remaining feasible points: X i=X g+0.5 (X c-X~).
(21)
This is repeated until Xi is feasible. Step 2: To carry out a reflection step, select the worst point X R where f ( X R) = m a x l f ( X ~, f ( X 2) . . . . . f ( x P ) ] . Then calculate the centroid, X c, of the generated points without including the worst point. A new point can be generated using the following equation:
X m=X'+6~(X ~-XR).
(22)
Check the feasibility of X m. If f ( X m) is feasible and f ( X m ) ~_ f ( X R ) , retract half the distance toward the centroid. Continue this process until f ( X m) < f ( x R ) , and go to Step 4. If X m is infeasible, go to Step 3. Step 3: If x m < x/, where x/is the lower limit for the variable xi, set x m = x/. If x~ >_ xy, where xy is the upper limit for the variable xi, set x m = xj~. If the resulting x/~ is still infeasible, retract half the distance toward the centroid until it is feasible. Step 4: To check whether the desired optimization has been achieved, determine the values of the following parameters:
Cohesionless soil model
317
P
~_j f(xi) . ~ i=l
(23)
P P
~ xt i=1
(24)
2i-----~ P
The iterations within the algorithm are terminated when the following conditions are satisfied P
>- ~_j [f(xi) -f(x~)] 2
(25)
i=A P
~ ~ Ix~ - J?~I.
(26)
i=1
Since the Box algorithm is designed for a minimization type problem, the utility function optimization problem is modified to minimize the error subject to the parameters X within reasonable limits. The vector X mentioned above refers to the vector of ten parameters which need to be optimized. The error is defined as follows: ~-~j y i 2 -- ~
ERR = 1 . 0 -
/=1
i=1 m
(..Vie -- Yip) 2
,
(27)
Z i=1
where m is the number o f points on the strain axis at which vertical error is calculated and Yip, Yie are the corresponding deviatoric stress values for the predicted and experimental observations (Fig. 3). Since all practical optimization strategies are designed only to locate a local minimum in the vicinity o f the initial search position, it is generally necessary to perform several trials from different starting points to establish the true global minimum. In the present work also, several trials were made to achieve the global minimum. All of the computer trials were run in double precision to ascertain maximum possible accuracy.
VI. E V A L U A T I O N O F M A T E R I A L P A R A M E T E R S
The proposed constitutive model has a total of fourteen material constants: K, G, 31, ~l, 32, 72, M, N, #, R, A, a, 3, and 3', where K and G are, respectively, the elastic bulk and shear modulii. VI. 1. Determination o f K and G The bulk modulus K and the shear modulus G are determined using eqns (28) and (29), by plotting K and G as a function o f mean pressure, p (AnDt~_JAtrWAO et ai.,
318
A. AaDULRAm~EM et al.
f
J J
/
j/ ////
3
/ Y,
le
Predicted Experimental
I Strain Fig. 3. Error definition.
[1991]). Both K and G (in kPa) are found to be a linear function of p, within the range of values considered here.
K ( p ) = 27369 + 407.7 p
(28)
G(p) = 9598 + 317.6 p.
(29)
Note that p = oiJ3 in the above equations is the mean pressure in kPa. VI.2. Determination o f M Equation (1) can be written as
M =
J'f~-zo/lJ~g(O,J,)].
(30)
The value of M for each CTC, TC, or TE test can be determined by substituting the value of Jl and J2D at failure in eqn (30). An average of the M values is used in the model. VI.3. Determination o f N and Equation (7) can be rewritten as ln[1
"~z/~ ] = l n N - - - ~J~ P,, MJig(O,J1)
(31)
Cohesionless soil model
319
All the parameters in the above equation are known except ~ and N. The values of J1 and JzD corresponding to the points where dilation starts (as shown by points a and a' in Fig. 1) for different tests, are used in the above equation to yield different values for left hand side of eqn (31). The resulting data is then regressed to yield the values of /t a n d N . VI.4. Determination o f ~31 and ~l The parameters/~1 and 7/1 are evaluated by fitting eqn (18) to the corresponding values obtained from the hydrostatic compression tests. X =/31 ~ '
(32)
where ,% = l
(33)
VI.5. Determination of/32 and 72 The hardening function as given by eqns (16) and (17) can be rewritten as In [ In(X//31) [ ~1 ~-~
] 1 = In/32 + 7/2In rd.
(34)
Using a shearing path of a CTC, TC, or TE test, the values of rd, ~, and X are evaluated at different points along the shearing path. The data is then regressed to yield the values of/32 and 72.
VI.6. Determination o f R,.4, c~, /3, and "y These five parameters are determined using the optimization technique described earlier. The initial values and the range of each parameter have to be selected before using the optimization technique. This requires some knowledge about the material behavior and familiarity with the constitutive model. It may be difficult to reach an optimized set if the starting point and the range of parameters are not chosen well. The global minimum could not be established successfully from all starting points (DEN~AI~, [1983]).
VII. APPLICATION A computer procedure is developed to facilitate the prediction of stress-strain and volumetric response of cohesiouless soils. Laboratory test data for a fine dune sand is used here to evaluate the material constants and to verify the proposed constitutive model. A summary of the test results for this sand is given in Table 1 (AnDtrUAOWAD et al.,
[1991]).
320
A. ABDULRAHEEMet aL
Table 1. S u m m a r y of test results for the fine dune sand M a x i m u m density ( k N / m 3)
17.71
M i n i m u m density ( k N / m 3)
14.73
Measured density ( k N / m 3)
16.40
Average relative density (%)
81.0 0.15
Dlo (mm) D60 Uniformity coefficient Cu = - -
1.84
/gjo
D~
Coefficient of curvature C c = - Dio * D6o Angle of friction (~) Cohesion 'c' ( k N / m m 2)
0.998 33.0 ° 0.0
By using the parameter evaluation methods described in the previous section, the "fixed" parameters are obtained as M = 0.2, N = 0.363, ~ 1 = 1.61 * l06, */l -- 1.17, 32 -- 0.442 and ~'/2~ 2.48. Since a linear form of the second characteristic state line is used in this study, # is taken as zero. The "free" parameters are optimized using only one conventional triaxial compression (CTC) test data at 30 psi (207 kPa), and all other tests are back-predicted using the same set of optimized parameters. Figure 4 shows the CTC test data which is used for evaluation of material constants using the optimization technique discussed earlier. The parameters are found as R = 7.14, A = 1.66, c~ = 0.218, 3 = 0.645, and "y -- 0.639. It is evident from Fig. 4 that the response predicted by the proposed model is in excellent agreement with the experimental data. The above set of optimized parameters are used to predict the data corresponding to different stress paths and different confining pressures. A comparison o f the predicted stress-strain, (~l - aa) vs ~1 and (ol - 03) vs e3, and volumetric response, ev vs el with the experimental data for CTC tests performed at confining pressures of 20 psi (138 kPa), 35 psi (241 kPa), and 40 psi (276 kPa), are presented in Figs. 5-7. The predicted response is in good agreement with the experimental data for all the tests. At higher values o f axial strains (> 1°/0), some deviations between the predicted and experimental deviatoric stresses are observed. The observed stress-strain and volumetric-axial strain responses for a triaxial compression (TC) test at a confining pressure of 20 psi (138 kPa) are compared with the model predictions in Fig. 8. The volumetric behavior in the dilation zone shows some deviation from the experimental data. However, the prediction for the stress-strain behavior, the point of the onset of dilation, and the response in the compression region is excellent.
VIII. C O N C L U D I N G R E M A R K S
A constitutive model is developed to characterize the stress-strain and volumetric behavior of cohesiordess soils. The model utilizes the concept of two characteristic states. The two states are represented in the stress space as the two characteristic state lines rep-
Cohesionless soil model
800
.
.
.
.
.
.
.
.
i
.
.
.
.
!
321
.
.
.
.
.
.
.
.
Predicted (Axial Strain) o- . . . . e Experiment (A~dal Strain) ..... Predicted (Radial Strain) ~__~'-_............ ~...... B-----O F_xpenment (Radial Strain) .....-I ................. ~ ~ ' ~ i
t~
o. v
600
~)
400 0
-t-
>
t-t
200
................................ ! .......................... ~i.~
0 -0.02
....
i ....
......................... ~.................................... ~ ..............................
~ ....
-0.01
0
] .... 0.01
i .... 0.02
0.03
Radial and Axial Strain (a) 0.0050
.
.
.
.
!
0.0025
.
.
.
.
.
.
.
Experimentaldata
.=_
j.V"
i
j"
E _= 0 > -0.0025
-0.0050 0.005
0.010 Axial
0.015
0.020
0.025
Strain
(b)
Fig. 4. (a) Comparison of observed and predicted stress-strain response for CTC test at conf. pr. of 30 psi. (b) Comparison of observed and predicted volumetric-axial strain response for CTC test at conf. pr. of 30 psi.
resenting the state o f cohesionless soil at failure, and the state at which the soil passes f r o m the compressive m o d e to the dilative m o d e o f d e f o r m a t i o n during shearing. The suggested special f o r m o f the CSL-2 allows the model to predict dilative response at low and medium confining pressures and purely compressive response at high confining pres-
322
A. ABDULRAHEEMet al.
600 .... o- . . . . e ..... G------E~
ta
400 ~
t3--O
Prediciton (Axial strain) Experimental data (Axial strain) Prediction (Radial strain) Experimental data (Radial strain)
- ~.
[
.
.
.
~
i
¢/) O .tO o~ >
200
................
a
0 -0.02
).........Y
-0.01
. . . . . . . . .
0
0.01
0.02
Radial and Axial Strain (a)
0.004
......
0.002
............................!...............................
Predicted data E x p e ~ e n t a l data
•
::
:Um .................................................................................................
E
" ' . • m • gill •
o >
m
~ ; ....
mr-
...........7; .........................................
.-"
•
. . . . . . . ~. a ._=.__=.. . . . . . . . !-" -0.002
-0.004
. . . . . . . . . .
=
0
=
i
i
a
i
0.005
*
=
i
i
......
=
0.010
~
i
i
,
0.015
,
=
i
0.020
Axial Strain (b) Fig. 5. (a) C o m p a r i s o n o f observed and predicted stress-strain response for C T C test at conf. pr, of 20 psi. (b) Comparison o f observed and predicted volumetric-axial strain response for C T C test at eonL pr. of 20 psi.
sure. In addition, the inclusion of a distinct dilative surface with nonassociated formulation provides much better control in predicting volumetric behavior of cohesionless soils. A computer-aided automated procedure, which involves Simplex optimization technique, is used in this work to evaluate the "free" parameters. This approach reduces
Cohesionless soil model
.
800
.
.
.
.
.
.
.
i
.
.
.
i
.
.
.
.
i
----- Prediction (Axial strain) .................................. o- .... e Experimental data (Axial strain) ..... Prediction (Radial strain) ---_.~g--'-~
to
.
323
.
.
.
.
i I i................................... I!
Exp~rnental date (Radial strain) Ii
600
o
400
-r"
> ¢D t:2 ,
m
200
0
-0.02
-0.01
0
0.01
0.02
0.03
Radial and Axial Strain (a) 0.0050
0.0025
......................... ....... •
._.q
o
ExperimentaIPredi datadata cted ill............................................ i!..................................... • •
i ...................................... 0 ......................................
b
E
; .Era=t= •;
o ~>
,
"
•
i
~.'"
:, i; .........j::; ................... =i ...- .-
,.•
-0.0025
-0.0050
.
0
.
.
.
=
0.005
,
i
,
=
=
0.010
=
=
,
0.015
,
=
i
|
i
0.020
,
,
,
,
0.025
Axial Strain (b) Fig. 6. (a) Comparison o f observed and predicted stress-strain response for CTC test at conf. pr. o f 35 psi. (b) Comparison o f observed and predicted volumetric-axial strain response for CTC test at conf. pr. of 35 psi.
the complexity o f the model calibration process and increases the accessibility and usefulness o f the proposed model significantly. The model is verified with respect to laboratory test data for a fine dune sand used in the parameter evaluation. The model predictions are carried out using a number o f other data not used in evaluating the con-
324
A.
1000
.
.
.
.
.
.
.
. .
! .
.
" ~ "
n
u)
. . . .
i
. . . .
.
I o- . . . .
800 ~ . . . . .
et al.
ABDULRAHEEM
e Experimental data (Axial strain)
/
] ....... Prediction(Radial strain)
I
/ B------El Experimentaldata(Radialstrain)
/
i
~:
600
09 Q L..
o
. . . . . . .
400
%1
....
,"
,
,
t~
a 200
........................................................................~j': ........................ i ............................. i .............
0 -0.02
-0.01
0
0.01
0.02
Radial and Axial Strain (a) 0.010
T
................................
0.005
•. . . . . . •
Predicted data
[ .............................................................................
Experimental data
1
o .r-
~ i i l i l
E
•
and I
i~i
•
•
_. . . . . . .
!
o
>
:
i }
-0.005
i -0.010
.
0
.
.
.
i
.
0.005
.
.
.
i
i .
0.010
.
.
.
=
.
0.015
.
.
.
i
.
0.020
.
.
.
0.025
Axial Strain (b) Fig. ?. (a) Comparison of observed and predicted stress-strain response for TC test at conf. pr. of 40 psi. (b) Comparison of observed and predicted volumetric-axial strain response for CTC test at conf. pr. of 40 psi.
stants. G e n e r a l l y , a g o o d a g r e e m e n t is o b s e r v e d b e t w e e n the m o d e l p r e d i c t i o n s a n d the experimental data. T h e c o n c e p t i n t r o d u c e d here c a n p r o v i d e a b a s e l i n e m o d e l for i n c l u s i o n o f s o m e i m p o r t a n t f a c t o r s such as relative density, o r i e n t a t i o n o f stress p a t h , a n i s o t r o p y , a n d s o f t e n i n g b e h a v i o r in d e v e l o p i n g c o n s t i t u t i v e m o d e l s f o r cohesiordess soils.
Cohesionless soil model
.
.
.
.
.
.
.
.
!
o- . . . . e ~--- - - ~ ...... .......
300
........................
200
":........:"-=-:-":="~"~-
,.
•
•
•
325
i
.
.
.
.
Expedment~ data (Axial strain) Expenrnent~ data (Radial strain) Prediclion (Radial strain) Prediclion (Axial strain)
......................................
Q. v
i ........... ~.-~":~'~::::
::-:::"-":
- t0' - -
Q
100 (.
? 0 -0.02
-0.01
0
0.01
0.02
0.03
Radial and Axial Strain
(a) 0.0050
0.0025 j"
._¢
iI i i
1" t
"'~'"
I "
i !
J
~=~ ~:~:~
E
=
..................... ................................................................. i......................................
o >
[ ~--•
-0.0025
-0.0050
.
.
.
.
. . . 0.005
.
.
Prediction Experimental data
i i.
. . 0.010
.
.
.
. . 0.015
[ . 0.020
Axial Strain (b) Fig. 8. (a) Comparison of observed and predicted stress-strain response for TC test at conf. pr. of 20 psi. (b) Comparison of observed and predicted volumetric-axial strain response for TC test at conf. pr. of 20 psi.
Acknowledgements-The authors are grateful to Prof. N. Cristescu, University of Romania, for very helpful comments during his visit to the University of Oklahoma. Financial support for this work was provided in part by the Oklahoma Center for Advancement of Science and Technology (ORAL 125-8836) and the National Science Foundation through the Rock Mechanics Research Center (ORAL 125-5861). Technical assistance provided by Mr. Azecmuddin is gratefully acknowledged. Discussions with Dr. Abdelaziz Mohammed on the Simplex Algorithm provided by him also significantly benefitted this work.
326
A. ABDULRAHEEMet al.
REFERENCES 1951 1957 1962 1965 1971 1972 1975 1975 1976 1977 1980 1980 1981 t982 1982 1982 1983
1983 1984 1984 1985 1985 1986 1987 1987 1989 1989 1991
PRAGER,W. and HOlmE, P.G., JR., "Theory of Perfectly Plastic Solids," John Wiley and Sons, New York. DRUCKER,D.C., GIBSON, R.E., and HENKEL, D.J., "Soil Mechanics and Work Hardening Theories of Plasticity," Proc. ASCE, 122, 338. SPENDLEY,W., HEXT, G.R., and HIMSWORTH, F.R., "Sequential Application of Simplex Designs in Optimization and Evolutionary Design," Technometrics, 4, 441. NELDER,J.A. and MEAD, R., "A Simplex Method for Function Minimization," Computer Journal, 7, 308. DIMAC_~10,F.L. and SANDLER,l.S., "Material Model for Granular Soils," J. Eng. Mech. Div., ASCE, 97, No. EM3,935. SWANN,W.H., "Direct Search Methods in Numerical Methods for Unconstrained Optimization," in MURRAY, W. (ed.), Academic Press, London, pp. 13-28. LADE,P.V. and DUNCAN, J.M., "Elastoplastic Stress-Strain Theory for Cohesionless Soil," J. Geotech. Eng. Div., ASCE, 101, No. GTI0, 1037. ZIENKIEWlCZ,O.C., H ~ r m S O N , C., and LEwis, R.W., "Associated and Nonassociated Viscoplasticity and Plasticity in Soil Mechanics," Geotechnique, 25, No. 4, 671. SANDLER,I.S., DIMAC~IO, EL., and BALADI, G.Y., "Generalized Cap Model for Geological Materials," J. Geotech. Div., ASCE, 102, No. GT7, 683. LADE,P.V., "Elastic-Plastic Stress-Strain Theory for Cohesioniess Soil with Curved Yield Surfaces," International Journal of Solids and Structures, 13, 1019. FLETCHER,R., "Practical Methods of Optimization-Volume l: Unconstrained Optimization," John Wiley and Sons, Chichester, U.K. LOUNG,M.P., "Stress-Strain Aspects of Cohesionless Soils Under Cyclic and Transient Loading," in Proc. Int. Symposium on Soils Under Cyclic and Transient Loading, PANDE, G.N. and ZIENKIEWICZ, O.C. (eds.), A.A. Balkema, Rotterdam, Vol. 1, 315. FLETCHER,R., "Practical Methods of Optimization-Volume 2: Constrained Optimization," John Wiley and Sons, Chichester, U.K. GHA~OUSSl,J. and MOMEN, H., "Modelling and Analysis of Cyclic Behavior of Sand," in PANDE, G.N. and ZmNKmWlCZ, O.C. (eds.), Soil Mechanics-Transient and Cyclic Loads, Wiley, p. 313. DArALIAS,Y.E and HARRUANN, L.R., "Bounding Surface Formulation of Soil Plasticity," in PANDE, G.N. and ZIENKIEWlCZ, O.C. (eds.), Soil Mechanics-Transient and Cyclic Loads, p. 253. VALAmS,K.C. and READ, H.E., "A New Endochronic Plasticity Model for Soils," in PANDE, G.N. and ZmNKIEWlCZ, O.C. (eds.), Soil Mechanics--Transient and Cyclic Loads, p. 375. DENATALE,J.S., "On the Calibration of Constitutive Models by Multivariate Optimization. A Case Study: The Bounding Surface Plasticity Model," Ph.D. thesis, University of California, Davis. RELAKITIS,G.V., RAVINDRAN,A., and RAGSDELL,K.M., "Engineering Optimization: Methods and Applications," John Wiley and Sons, New York. DESAI,C.S. and FARUQUE, i . O . , "Constitutive Model for (Geological) Materials," Journal of Engineering Mechanics, ASCE, 110, 1391. FUKUSHIMA,S. and TATSUOKA,F., "Strength and Deformation Characteristics of Saturated Sand at Extremely Low Pressures," Soils and Foundations, 24, No. 4, 30. FARUQUE,M.O. and DESAI, C.S., "Implementation of a General Constitutive Model for Geological Materials," Int. J. Num. Analyt. Meth. Geomech., 9, 415. POOROOSHASB,H.B. and PmTRUSZCZAK,S., "On Yielding and Flow of Sand; A Generalized Two Surface Model," Computer and Geotechnics, 1, 33. DESAI,C.S., SOMASUNDARAM,S., and FRANTZISKONIS, G., "A Hierarchical Approach for Constitutive Modelling of Geologic Materials," Int. J. Num. Analyt. Meth. Geomech., 10, 3, 225. FARUQUE,M.O., "A Third Invariant Dependent Cap Model for Geological Materials," Soils and Foundations, 27, No. 2, 12. HIRAI,H., "A Combined Hardening Model for Plasticity for Sands," 2nd Int. Conf. Short Course on Constitutive Laws for Engineering Materials: Theory and Applications, DESAI, C.S. et al, (eds.), Elsevier, New York, 557. CmSTESCU,N., "Rock Rheology," Kluwer Academic Publishers, The Netherlands. HARDIN,B.O., "Low-Stress Dilation Test," J. Geotech. Eng., ASCE, 115, No. 6, 769. ABDLrLIA~AD,S.N., FARUQUE,M.O., and AZEEMUDDIN,M., "Experimental Verification and Numerical Evaluation of the Three-Invariant Dependent Cap Model for Cohesionless Soil," The Arabian Journal for Science and Engineering, 16, 107.
School of Civil Engineering and Environmental Science The University of Oklahoma, Norman Norman, OK 73017, USA (Received in final revised form 14 June 1993 )
International Journal of Plasticity, Vol. 10, No. 3, pp. 309-326, 1994 Copyright © 1994 Elsevier ScienceLtd Printed in the USA. All rights reserved 0749-6419/94 $6.00 + .00 0749-6419(93)E0006-M
A TWO CHARACTERISTIC MODEL FOR COHESIONLESS SOIL AND ITS CALIBRATION USING OPTIMIZATION
ABDULAZEEZ ABDULRAHEEM, M . OMAR FARUQUE, a n d MUSRARRAr ZAU_~',; The University of Oklahoma, Norman
Abstract--An elasto-plastic constitutive model is developed to characterize the stress-strain behavior of cohesionless soils. The model utilizes the concept of two characteristic states represented by two characteristic state lines in the stress space. The two characteristic state lines represent the state of the soil at failure and the state at which the soil passes from compressive to dilative mode of deformation during shearing. The proposed approach provides a better control in predicting volumetric behavior of cohesionless soils. To increase the accessibility and usefulness of the model and reduce the difficulty of the model calibration process, an optimization procedure is used to evaluate the material parameters associated with the model. The model is verified with respect to a test that is used for finding the parameters, and a number of other tests that are not used for finding the parameters. A good correlation between predictions and experimental data is observed.
!. INTRODUCTION
In the past, plasticity theory has been used extensively to formulate constitutive models for cohesionless soils (DRtrCKER et al., [1957]; DI~L~GGIO~ SASDLER, [1971]; LADE & DU~rCAr~, [1975]; Zms~:mw~cz et aL, [1975]; SASDLER, et aL, [1976]; DAr,~d~tAS& HAe,RUAr~N, [1982]; VAI.ASIS & READ, [1982]; D~NAT~E, [1983]; D~S~ ~ FARUQUE, [1984]; FARUQUE& DESAI, [1985]; DES~ etal., [1986]; FARUQtrE, [1987]; Hmxt, [1987]; Am)t~_JAUWAD et al., [1991]). Many of these models characterize stress-strain response in a reasonable manner. Quite frequently, the volumetric behavior is ignored or less emphasized by these models and sometimes it is not used in the model validation process. As evident from past experimental investigations, induced volume change during sheafing is an important characteristic of cohesionless soil. The nature of volume change is primarily controlled by the relative density of the soil and the confining pressure at which the shearing takes place. Depending upon the relative density and the applied confining pressure, a soil can experience either dilation or compression during shearing (LADE, [1977]; LOUNO, [1980]; GRAsoussI ~ M o ~ r , [1982]; FUKUSrm~A~, TArstror~, [1984]; POOROOSnASn~ PmXRUSZ~ZAK, [1985]; ASDULJAtrWADet aL, [1991]; H ~ I N ,
I1989]). Figure 1 shows the effect of relative density and the confining pressure on axial-volumetric response obtained from (stress controlled) triaxial compression tests. Evidently, a relatively dense soil at a low confining pressure experiences a very little shear-induced soil compression and in most part it dilates. On the other hand, at a very high confining pressure, cohesionless softs usually experience compressive volume change all through shearing up to failure. This behavior seems to be exhibited irrespective of the initial relative density of the soil. 309
A. ABDULRAI-IEEM et al.
310
High Confining Pressure
'
\el_
~.'"
1/ 2 /
Low Confining P,'~ssure
",\///'2
Radial Strain
Axial Strain
(a)
Dilatant
I /
Compressive
. -
"q ' -.
Strain
1', 2' Compressive
(b) Fig. 1. (a) Typical stress-strain response of cohesionless soil. (b) Typical volumetric-axial response of cohesionless soil.
This paper presents a general constitutive model for cohesionless soils based on the framework of plasticity theory. One of the unique features of this model is that it incorporates explicitly the two characteristic states of a cohesionless soil in the formulation. The proposed framework, in general, can be used to characterize shear-induced volume change behavior of these soils in a rational manner. In addition, to make the calibration process of the model easier, an optimization technique is also presented in detail. II. P R O P O S E D CONSTITUTIVE MODEL
Formulation of the proposed constitutive model is based on the nonassociated plasticity theory. Two distinct characteristic state lines are proposed to represent the two dif-
Cohesionless soil model
311
ferent characteristic states of cohesionless soils. The first characteristic state line (CSL-1) refers to the state of cohesionless soil at failure and is expressed as
J'~2o MJ~ ,
(1)
=
where M is a parameter related to the coefficient of internal frictional, which in turn depends on the relative density of the material. The second characteristic state line (CSL-2) refers to the state of cohesionless soil at which the rate of volumetric strain momentarily vanishes as the soil passes from compressive mode of deformation to dilative mode of deformation during shearing. Experimental evidence suggests that as the confining pressure increases, the dilation decreases and the volumetric response becomes predominantly compressive. At high confining pressure, dilation is not encountered at all. This phenomenon of decreasing dilatancy with increasing confining pressure can be modeled effectively by letting CSL-2 get closer to CSL-1 with increasing confining pressure. Consequently, the equation of this line is proposed as = M[ 1 - Nexp( -#J~/P.)] J~.
(2)
The exponential term in eqn (2) becomes negligible at higher values of Jl. As such eqn (2) approaches eqn (1) with increasing Jl. Typical plots of eqns (1) and (2) on the Jl - J~-)-22oplane are shown in Fig. 2. Efforts have been made in the past to model the onset of dilatancy and subsequent dilatant behavior by including the second characteristic state line, which essentially diverges from the first characteristic state line as the confining pressure increases (Gi-L~OUSSi~ Mor,m~, [1982]; FLrKuSrnMA~ TATSUOIOt, [1984];
os ~°"
.......
CSL-1
~-~ ......
CSL-2: LinearfolTn CSL-2: Nonlinear form
....
~
.....
Dilative yield surface
p
~
-'"
yield Surface
I"
s~
..
.oo
A
a
4. 4
o,,x-" .,,~"
4
|
/
~
\ C
X
J~ Axis Fig. 2. Schematic representation of the yield caps and the two CSLs.
312
A . ABDULRAHEEM
et al.
DESAI et al., [1986]; CRISTESCU, [1989]). The form of CSL-2, as suggested by eqn (2), is such that it yields a predominantly compressive response at high confining pressure and both the compressive and dilative responses at low and medium confining pressure ranges. Further, the effect of relative density can also be incorporated by making the parameters M and N functions of relative density. However, since experimental data used in this study are such that they fall approximately in the lower confining pressure range, eqn (2) is simplified to a linear form by assuming/, = 0. The simplified equation is given by (3)
~ 2 D = M[1 - N I J , .
Figure 2 shows the simplified linear form of the CSL-2 along with the nonlinear form as represented by eqn (2). Consistent with the definition of CSL-2, all states of stress below the second characteristic state line yield compressive volumetric strain only. A compressive yield surface of elliptical shape that satisfies the above requirements is proposed for the model and is shown schematically in Fig. 2. The tangent to the yield surface is vertical at the intersection with J~-axis and horizontal at the intersection with CSL-2. The normal to the yield surface at any point between the above two intersections has nonnegative slope. The equation of the compressive yield surface is proposed as
Fe = J~ZD -- \/a2 _ ()/ - - C )
2/~R2 = 0
,
(4)
where R is a material constant associated with the ratio of the axes of the ellipse and C is a response function that denotes the value of J~ at the point of intersection of a generic yield surface and the second characteristic state line. In eqn (4), a represents the axis of the elliptical compression yield surface parallel to the JX~2Daxis and is defined in terms of C as a = M[1 -- N e x p ( - ~ C / P a ) ] C.
(5)
The proposed model utilizes a nonassociated formulation to describe the volumetric response in the compressive regime. A plastic potential function, Qc, is defined in the form Qc =
J'~ZZD-- A ~ ~2 -- (J, - C) 2/R2,
(6)
where ,4 is a material constant. Taking the derivatives of Fc and Qc with respect to Jr, the following relationship can be obtained
(7)
aQ¢ _ `4 OF,,
As evident from eqn (7), the constant .4 signifies deviation from normality and should have a positive value in the range 0 ___,4 _< 1. The dilation yield surface refers to the stress states that cause dilating response at a material point. It is bounded by CSL-1 and CSL-2 at all times and identifies with CSL-1 when the confining pressure is very high. An equation of the dilation yield surface that satisfies the above requirements is proposed in the form Fa = J'~ZD -- a[1 -- (1
-
J 1 / C ) MC/a]
(8)
Cohesionless soil model
313
A schematic representation of a generic dilation yield surface is shown in Fig. 2. Evidently, the tangent to the dilation yield surface at the point of intersection with the second characteristic state line is horizontal, and thereby satisfies the slope compatibility with the compression yield surface, Fc = 0. At Jl = 0, the slope of the dilation yield surface is M, and therefore, tangential to the first characteristic state line at that point. Nonassociated formulation is used here to achieve zero rate of dilatancy when the stress state is on the characteristic state lines (both CSL-1 and CSL-2), and nonzero rate of dilatancy when the state of stress lies in between the characteristic state lines. To achieve this, a dilating potential function, Qd = 0, is defined by the following equations:
oo. =q(j,,s47Vo)af4 ,,Jl OQa 0J2D
-
OF,l 0J2D' - -
(9)
(10)
where q(Jl,J2D) = ct (/3 /3t37. + 7)~+"~ Y (1 - G)t~G v
(11)
M J, -G ( J1, J2D ) = M N exp( - M J l /P,,)Jl "
(12)
in which
In the above equations, o~,/3, and 7 are material constants. The function q(J1,J2D) assumes a value of zero at the boundaries (CSL-1 and CSL-2) and is nonnegative elsewhere. Explicit form of Qd is not required since the derivatives of Qa suffice to obtain the elasto-plastic constitutive matrix. The inclusion of a distinct dilating yield surface along with the nonassociated formulation, in terms of material constants a,/3, and 7, provides a much better control on the dilative response of the material. However, depending on the response of the soil, the three parameters, or,/3, and 7, can be reduced to one (i.e. or). Even an associated formulation can be adopted, if the nature of response dictates such simplifications, in both the compressive and the dilative regimes, in which case the total number of parameters will be reduced by four. It should be noted here that the simplification of the CSL-2 from nonlinear to linear form, as represented by eqn (3), does not affect any of the propositions mentioned above.
III. HARDENING BEHAVIOR
In this work, modeling of soil response under monotonic loading is attempted using the concept of isotropic hardening plasticity. Both the compression yield surface, Fc, and the dilation yield surface, Fa, are allowed to expand in the stress space to account for hardening during the elasto-plastic deformation. It is important to realize that the mechanism of inelastic volume change during hydrostatic compression is quite different from that in shear. As a result, a hardening function, with its constants determined
314
A. ABDULRAHEEMet al.
from hydrostatic compression test data, generally cannot describe elasto-plastic stressstrain and volumetric behavior during shearing in a rational manner. This problem can be resolved by expressing the hardening function in terms of two parameters, ~ and rd, which are defined in terms of the incremental plastic strain tensor de~ and the incremental deviatoric plastic strain tensor deft, as follows:
= f (de p deiPj) '/e
(13)
r d = (d/~,
(14)
~d = f (dei~ dei]) l/z.
(15)
where
In the above equations, the symbol f denotes the history. Note that ~d and rd have zero value throughout hydrostatic compression. Referring to Fig. 2, the compressive yield surface, Fc, originates from the Jl-axis (i.e. point X) and terminates at a point on the second characteristic state line (CSL-2) where Jl assumes a value of C. The dilating yield surface, Fd, on the other hand, starts from the origin 0 and terminates at CSL-2 where J~ = C. Since X = C + Ra, where R is the ratio of major to minor axis of the elliptical cap (Fig. 2), and a is the ordinate corresponding to J1 = C, expansion of both Fc and Fd can be described conveniently by making X a (hardening) function of the parameters ~ and r d as
X = 31~ ~,
(16)
where = ~l[1 + 32r321.
(17)
In eqns (16) and (17), ~ , ~2, ~11, and ~72 are material constants. An initially isotropic material, when subjected to a hydrostatic loading, experiences a spherical deformation only (i.e. rd = 0 and ~ = ~v, where ~v = f 1 / , f 3 d e p = 1 / x ~ p, e~ being the total volumetric plastic strain). The hardening function for this case reduces to
x= ~'=
~(c~/4~).,.
(18)
IV. ELASTO-PLASTIC CONSTITUTIVE RELATIONS
Using the concept of nonassociated plasticity theory, incremental plastic strain tensor, de p, can be expressed as de p = dA O_Q_Q OCTij '
(19)
Cohesionlesssoil model
315
where Q is the plastic potential function and dA is an unknown scalar to be determined from the consistency condition of Prager, D f = O, F being the yield function. Following the standard steps of the theory of plasticity (PRAo~R & HODOE [1951]), the elastoplastic constitutive relation tensor, CvPktcan be written as
aF
0Q
C~.,, aa.o aa,.. C~,.kt Cu.~, =
~F
e
a_QQ _
,
(20)
aapq C~rs c~crrs A(~,rd, a~s) where Cfjk~ is the elastic constitutive tensor, and A(~,ra,aij) is a measure of plastic modulus involving derivatives of the yield function F with respect to ~ and ra, and derivative of the potential function Q with respect to a;j. Equation (20) is valid in the compressive as well as in the dilative regimes of shear deformation provided Q is appropriately defined. In the compressive regime, Q = Q¢ (eqn 6), while in the dilative regime, Q = Qa, defined implicitly by eqns (9)-(12).
V. C A L I B R A T I O N OF T H E M O D E L
The proposed model involves different material parameters, which are broken down into directly measurable or "fixed" parameters and the remaining "free" parameters. Determination of the directly measurable or fixed parameters is quite straightforward, and it can be readily accomplished by using well-established laboratory experiments. Determination of the remaining free parameters, however, can make the model calibration procedure quite difficult. These free parameters are generally established by a trial and error method, with the objective being to obtain the best overall fit to a given set of experimental data (Dar~tAs ~ H,~aL~a,m, [1982]). As a result, the accuracy and efficiency of the calibration process can be strongly dependent on the subjectivity of the user as well as the user's knowledge about the material and expertise with the particular constitutive model. In view of this, a computer-aided procedure is developed by casting the model calibration problem into a form which makes it amenable to solution by conventional mathematical optimization techniques. Because this automated procedure greatly reduces the dependence of calibration success on user expertise, it significantly increases the accessibility and usefulness of this model to the general engineering community. V.1. Simplex optimization technique Simplex approach of function minimization is used here as a tool for the evaluation of the free parameters in an optimum manner (SPEI'ZOLEYet aL, [1962]; FLETCm~R, [1980]; FLETCHER, [1981]; D~.NxrALE, [1983]). The procedure is based on the evolutionary operation (EVOP) technique originally introduced by Box in 1957, in which the solution space is explored by means of geometrical configuration of points rather than a set of directions. A number of modifications to this basic approach have since been suggested, but the version attributed to Nelder & Mead, in which the simplex automatically rescales itself according to local geometry of the objective function, is commonly regarded as most effective (Nm~DER * M ~ , [1965]). The primary merits of the simplex method are related to its extreme simplicity and great generality. The procedure
316
A. ABDULRAHEEMet al.
requires only the ability to evaluate the objective function and, typically, makes no assumptions with respect to the behavior of the function and its derivatives. For this reason, it is particularly useful when the function is nondifferentiable, when the gradients are discontinuous, or when the value of the function depends upon some physical measurements and may therefore be subject to random error. As a result of its simplicity, the strategy can be readily implemented into a workable code. In addition, because the procedure normally employs no complicated numerical operations, it is traditionally immune from rounding o f f errors a n d / o r ill-conditioning. Furthermore, the routine is computationally compact and therefore requires only modest amounts of computer storage. According to SWANN [1972], the direct search methods, of which the simplex approach is most successful, have proved to be quite reliable in practice and only rarely fail to locate at least a local minimum. V.2. Numerical algorithm for simplex technique The simplex search, known as the Box algorithm (SPENDLEY et al., [1962]; RELAK~TIS et al., [1983]), which can handle problems with many constraints, is used here due to the fact that it does not require any derivative of the objective function and constraints. The following steps are intended to explain the application o f the Box algorithm:
Step 1: Determine an initial starting point X ° (where X represents the vector containing the individual values denoted by x~, of the ten parameters which need to be optimized), a reflection coefficient &, and termination parameters i and 6. Generate a set of ( P - 1) random feasible points, where P = 2N, N b e i n g the number of decision variables in the problem. One needs to sample Ntimes to generate one point. If the point X i is found to be infeasible, retract half the distance toward the centroid, X c, of the remaining feasible points: X i=X g+0.5 (X c-X~).
(21)
This is repeated until Xi is feasible. Step 2: To carry out a reflection step, select the worst point X R where f ( X R) = m a x l f ( X ~, f ( X 2) . . . . . f ( x P ) ] . Then calculate the centroid, X c, of the generated points without including the worst point. A new point can be generated using the following equation:
X m=X'+6~(X ~-XR).
(22)
Check the feasibility of X m. If f ( X m) is feasible and f ( X m ) ~_ f ( X R ) , retract half the distance toward the centroid. Continue this process until f ( X m) < f ( x R ) , and go to Step 4. If X m is infeasible, go to Step 3. Step 3: If x m < x/, where x/is the lower limit for the variable xi, set x m = x/. If x~ >_ xy, where xy is the upper limit for the variable xi, set x m = xj~. If the resulting x/~ is still infeasible, retract half the distance toward the centroid until it is feasible. Step 4: To check whether the desired optimization has been achieved, determine the values of the following parameters:
Cohesionless soil model
317
P
~_j f(xi) . ~ i=l
(23)
P P
~ xt i=1
(24)
2i-----~ P
The iterations within the algorithm are terminated when the following conditions are satisfied P
>- ~_j [f(xi) -f(x~)] 2
(25)
i=A P
~ ~ Ix~ - J?~I.
(26)
i=1
Since the Box algorithm is designed for a minimization type problem, the utility function optimization problem is modified to minimize the error subject to the parameters X within reasonable limits. The vector X mentioned above refers to the vector of ten parameters which need to be optimized. The error is defined as follows: ~-~j y i 2 -- ~
ERR = 1 . 0 -
/=1
i=1 m
(..Vie -- Yip) 2
,
(27)
Z i=1
where m is the number o f points on the strain axis at which vertical error is calculated and Yip, Yie are the corresponding deviatoric stress values for the predicted and experimental observations (Fig. 3). Since all practical optimization strategies are designed only to locate a local minimum in the vicinity o f the initial search position, it is generally necessary to perform several trials from different starting points to establish the true global minimum. In the present work also, several trials were made to achieve the global minimum. All of the computer trials were run in double precision to ascertain maximum possible accuracy.
VI. E V A L U A T I O N O F M A T E R I A L P A R A M E T E R S
The proposed constitutive model has a total of fourteen material constants: K, G, 31, ~l, 32, 72, M, N, #, R, A, a, 3, and 3', where K and G are, respectively, the elastic bulk and shear modulii. VI. 1. Determination o f K and G The bulk modulus K and the shear modulus G are determined using eqns (28) and (29), by plotting K and G as a function o f mean pressure, p (AnDt~_JAtrWAO et ai.,
318
A. AaDULRAm~EM et al.
f
J J
/
j/ ////
3
/ Y,
le
Predicted Experimental
I Strain Fig. 3. Error definition.
[1991]). Both K and G (in kPa) are found to be a linear function of p, within the range of values considered here.
K ( p ) = 27369 + 407.7 p
(28)
G(p) = 9598 + 317.6 p.
(29)
Note that p = oiJ3 in the above equations is the mean pressure in kPa. VI.2. Determination o f M Equation (1) can be written as
M =
J'f~-zo/lJ~g(O,J,)].
(30)
The value of M for each CTC, TC, or TE test can be determined by substituting the value of Jl and J2D at failure in eqn (30). An average of the M values is used in the model. VI.3. Determination o f N and Equation (7) can be rewritten as ln[1
"~z/~ ] = l n N - - - ~J~ P,, MJig(O,J1)
(31)
Cohesionless soil model
319
All the parameters in the above equation are known except ~ and N. The values of J1 and JzD corresponding to the points where dilation starts (as shown by points a and a' in Fig. 1) for different tests, are used in the above equation to yield different values for left hand side of eqn (31). The resulting data is then regressed to yield the values of /t a n d N . VI.4. Determination o f ~31 and ~l The parameters/~1 and 7/1 are evaluated by fitting eqn (18) to the corresponding values obtained from the hydrostatic compression tests. X =/31 ~ '
(32)
where ,% = l
(33)
VI.5. Determination of/32 and 72 The hardening function as given by eqns (16) and (17) can be rewritten as In [ In(X//31) [ ~1 ~-~
] 1 = In/32 + 7/2In rd.
(34)
Using a shearing path of a CTC, TC, or TE test, the values of rd, ~, and X are evaluated at different points along the shearing path. The data is then regressed to yield the values of/32 and 72.
VI.6. Determination o f R,.4, c~, /3, and "y These five parameters are determined using the optimization technique described earlier. The initial values and the range of each parameter have to be selected before using the optimization technique. This requires some knowledge about the material behavior and familiarity with the constitutive model. It may be difficult to reach an optimized set if the starting point and the range of parameters are not chosen well. The global minimum could not be established successfully from all starting points (DEN~AI~, [1983]).
VII. APPLICATION A computer procedure is developed to facilitate the prediction of stress-strain and volumetric response of cohesiouless soils. Laboratory test data for a fine dune sand is used here to evaluate the material constants and to verify the proposed constitutive model. A summary of the test results for this sand is given in Table 1 (AnDtrUAOWAD et al.,
[1991]).
320
A. ABDULRAHEEMet aL
Table 1. S u m m a r y of test results for the fine dune sand M a x i m u m density ( k N / m 3)
17.71
M i n i m u m density ( k N / m 3)
14.73
Measured density ( k N / m 3)
16.40
Average relative density (%)
81.0 0.15
Dlo (mm) D60 Uniformity coefficient Cu = - -
1.84
/gjo
D~
Coefficient of curvature C c = - Dio * D6o Angle of friction (~) Cohesion 'c' ( k N / m m 2)
0.998 33.0 ° 0.0
By using the parameter evaluation methods described in the previous section, the "fixed" parameters are obtained as M = 0.2, N = 0.363, ~ 1 = 1.61 * l06, */l -- 1.17, 32 -- 0.442 and ~'/2~ 2.48. Since a linear form of the second characteristic state line is used in this study, # is taken as zero. The "free" parameters are optimized using only one conventional triaxial compression (CTC) test data at 30 psi (207 kPa), and all other tests are back-predicted using the same set of optimized parameters. Figure 4 shows the CTC test data which is used for evaluation of material constants using the optimization technique discussed earlier. The parameters are found as R = 7.14, A = 1.66, c~ = 0.218, 3 = 0.645, and "y -- 0.639. It is evident from Fig. 4 that the response predicted by the proposed model is in excellent agreement with the experimental data. The above set of optimized parameters are used to predict the data corresponding to different stress paths and different confining pressures. A comparison o f the predicted stress-strain, (~l - aa) vs ~1 and (ol - 03) vs e3, and volumetric response, ev vs el with the experimental data for CTC tests performed at confining pressures of 20 psi (138 kPa), 35 psi (241 kPa), and 40 psi (276 kPa), are presented in Figs. 5-7. The predicted response is in good agreement with the experimental data for all the tests. At higher values o f axial strains (> 1°/0), some deviations between the predicted and experimental deviatoric stresses are observed. The observed stress-strain and volumetric-axial strain responses for a triaxial compression (TC) test at a confining pressure of 20 psi (138 kPa) are compared with the model predictions in Fig. 8. The volumetric behavior in the dilation zone shows some deviation from the experimental data. However, the prediction for the stress-strain behavior, the point of the onset of dilation, and the response in the compression region is excellent.
VIII. C O N C L U D I N G R E M A R K S
A constitutive model is developed to characterize the stress-strain and volumetric behavior of cohesiordess soils. The model utilizes the concept of two characteristic states. The two states are represented in the stress space as the two characteristic state lines rep-
Cohesionless soil model
800
.
.
.
.
.
.
.
.
i
.
.
.
.
!
321
.
.
.
.
.
.
.
.
Predicted (Axial Strain) o- . . . . e Experiment (A~dal Strain) ..... Predicted (Radial Strain) ~__~'-_............ ~...... B-----O F_xpenment (Radial Strain) .....-I ................. ~ ~ ' ~ i
t~
o. v
600
~)
400 0
-t-
>
t-t
200
................................ ! .......................... ~i.~
0 -0.02
....
i ....
......................... ~.................................... ~ ..............................
~ ....
-0.01
0
] .... 0.01
i .... 0.02
0.03
Radial and Axial Strain (a) 0.0050
.
.
.
.
!
0.0025
.
.
.
.
.
.
.
Experimentaldata
.=_
j.V"
i
j"
E _= 0 > -0.0025
-0.0050 0.005
0.010 Axial
0.015
0.020
0.025
Strain
(b)
Fig. 4. (a) Comparison of observed and predicted stress-strain response for CTC test at conf. pr. of 30 psi. (b) Comparison of observed and predicted volumetric-axial strain response for CTC test at conf. pr. of 30 psi.
resenting the state o f cohesionless soil at failure, and the state at which the soil passes f r o m the compressive m o d e to the dilative m o d e o f d e f o r m a t i o n during shearing. The suggested special f o r m o f the CSL-2 allows the model to predict dilative response at low and medium confining pressures and purely compressive response at high confining pres-
322
A. ABDULRAHEEMet al.
600 .... o- . . . . e ..... G------E~
ta
400 ~
t3--O
Prediciton (Axial strain) Experimental data (Axial strain) Prediction (Radial strain) Experimental data (Radial strain)
- ~.
[
.
.
.
~
i
¢/) O .tO o~ >
200
................
a
0 -0.02
).........Y
-0.01
. . . . . . . . .
0
0.01
0.02
Radial and Axial Strain (a)
0.004
......
0.002
............................!...............................
Predicted data E x p e ~ e n t a l data
•
::
:Um .................................................................................................
E
" ' . • m • gill •
o >
m
~ ; ....
mr-
...........7; .........................................
.-"
•
. . . . . . . ~. a ._=.__=.. . . . . . . . !-" -0.002
-0.004
. . . . . . . . . .
=
0
=
i
i
a
i
0.005
*
=
i
i
......
=
0.010
~
i
i
,
0.015
,
=
i
0.020
Axial Strain (b) Fig. 5. (a) C o m p a r i s o n o f observed and predicted stress-strain response for C T C test at conf. pr, of 20 psi. (b) Comparison o f observed and predicted volumetric-axial strain response for C T C test at eonL pr. of 20 psi.
sure. In addition, the inclusion of a distinct dilative surface with nonassociated formulation provides much better control in predicting volumetric behavior of cohesionless soils. A computer-aided automated procedure, which involves Simplex optimization technique, is used in this work to evaluate the "free" parameters. This approach reduces
Cohesionless soil model
.
800
.
.
.
.
.
.
.
i
.
.
.
i
.
.
.
.
i
----- Prediction (Axial strain) .................................. o- .... e Experimental data (Axial strain) ..... Prediction (Radial strain) ---_.~g--'-~
to
.
323
.
.
.
.
i I i................................... I!
Exp~rnental date (Radial strain) Ii
600
o
400
-r"
> ¢D t:2 ,
m
200
0
-0.02
-0.01
0
0.01
0.02
0.03
Radial and Axial Strain (a) 0.0050
0.0025
......................... ....... •
._.q
o
ExperimentaIPredi datadata cted ill............................................ i!..................................... • •
i ...................................... 0 ......................................
b
E
; .Era=t= •;
o ~>
,
"
•
i
~.'"
:, i; .........j::; ................... =i ...- .-
,.•
-0.0025
-0.0050
.
0
.
.
.
=
0.005
,
i
,
=
=
0.010
=
=
,
0.015
,
=
i
|
i
0.020
,
,
,
,
0.025
Axial Strain (b) Fig. 6. (a) Comparison o f observed and predicted stress-strain response for CTC test at conf. pr. o f 35 psi. (b) Comparison o f observed and predicted volumetric-axial strain response for CTC test at conf. pr. of 35 psi.
the complexity o f the model calibration process and increases the accessibility and usefulness o f the proposed model significantly. The model is verified with respect to laboratory test data for a fine dune sand used in the parameter evaluation. The model predictions are carried out using a number o f other data not used in evaluating the con-
324
A.
1000
.
.
.
.
.
.
.
. .
! .
.
" ~ "
n
u)
. . . .
i
. . . .
.
I o- . . . .
800 ~ . . . . .
et al.
ABDULRAHEEM
e Experimental data (Axial strain)
/
] ....... Prediction(Radial strain)
I
/ B------El Experimentaldata(Radialstrain)
/
i
~:
600
09 Q L..
o
. . . . . . .
400
%1
....
,"
,
,
t~
a 200
........................................................................~j': ........................ i ............................. i .............
0 -0.02
-0.01
0
0.01
0.02
Radial and Axial Strain (a) 0.010
T
................................
0.005
•. . . . . . •
Predicted data
[ .............................................................................
Experimental data
1
o .r-
~ i i l i l
E
•
and I
i~i
•
•
_. . . . . . .
!
o
>
:
i }
-0.005
i -0.010
.
0
.
.
.
i
.
0.005
.
.
.
i
i .
0.010
.
.
.
=
.
0.015
.
.
.
i
.
0.020
.
.
.
0.025
Axial Strain (b) Fig. ?. (a) Comparison of observed and predicted stress-strain response for TC test at conf. pr. of 40 psi. (b) Comparison of observed and predicted volumetric-axial strain response for CTC test at conf. pr. of 40 psi.
stants. G e n e r a l l y , a g o o d a g r e e m e n t is o b s e r v e d b e t w e e n the m o d e l p r e d i c t i o n s a n d the experimental data. T h e c o n c e p t i n t r o d u c e d here c a n p r o v i d e a b a s e l i n e m o d e l for i n c l u s i o n o f s o m e i m p o r t a n t f a c t o r s such as relative density, o r i e n t a t i o n o f stress p a t h , a n i s o t r o p y , a n d s o f t e n i n g b e h a v i o r in d e v e l o p i n g c o n s t i t u t i v e m o d e l s f o r cohesiordess soils.
Cohesionless soil model
.
.
.
.
.
.
.
.
!
o- . . . . e ~--- - - ~ ...... .......
300
........................
200
":........:"-=-:-":="~"~-
,.
•
•
•
325
i
.
.
.
.
Expedment~ data (Axial strain) Expenrnent~ data (Radial strain) Prediclion (Radial strain) Prediclion (Axial strain)
......................................
Q. v
i ........... ~.-~":~'~::::
::-:::"-":
- t0' - -
Q
100 (.
? 0 -0.02
-0.01
0
0.01
0.02
0.03
Radial and Axial Strain
(a) 0.0050
0.0025 j"
._¢
iI i i
1" t
"'~'"
I "
i !
J
~=~ ~:~:~
E
=
..................... ................................................................. i......................................
o >
[ ~--•
-0.0025
-0.0050
.
.
.
.
. . . 0.005
.
.
Prediction Experimental data
i i.
. . 0.010
.
.
.
. . 0.015
[ . 0.020
Axial Strain (b) Fig. 8. (a) Comparison of observed and predicted stress-strain response for TC test at conf. pr. of 20 psi. (b) Comparison of observed and predicted volumetric-axial strain response for TC test at conf. pr. of 20 psi.
Acknowledgements-The authors are grateful to Prof. N. Cristescu, University of Romania, for very helpful comments during his visit to the University of Oklahoma. Financial support for this work was provided in part by the Oklahoma Center for Advancement of Science and Technology (ORAL 125-8836) and the National Science Foundation through the Rock Mechanics Research Center (ORAL 125-5861). Technical assistance provided by Mr. Azecmuddin is gratefully acknowledged. Discussions with Dr. Abdelaziz Mohammed on the Simplex Algorithm provided by him also significantly benefitted this work.
326
A. ABDULRAHEEMet al.
REFERENCES 1951 1957 1962 1965 1971 1972 1975 1975 1976 1977 1980 1980 1981 t982 1982 1982 1983
1983 1984 1984 1985 1985 1986 1987 1987 1989 1989 1991
PRAGER,W. and HOlmE, P.G., JR., "Theory of Perfectly Plastic Solids," John Wiley and Sons, New York. DRUCKER,D.C., GIBSON, R.E., and HENKEL, D.J., "Soil Mechanics and Work Hardening Theories of Plasticity," Proc. ASCE, 122, 338. SPENDLEY,W., HEXT, G.R., and HIMSWORTH, F.R., "Sequential Application of Simplex Designs in Optimization and Evolutionary Design," Technometrics, 4, 441. NELDER,J.A. and MEAD, R., "A Simplex Method for Function Minimization," Computer Journal, 7, 308. DIMAC_~10,F.L. and SANDLER,l.S., "Material Model for Granular Soils," J. Eng. Mech. Div., ASCE, 97, No. EM3,935. SWANN,W.H., "Direct Search Methods in Numerical Methods for Unconstrained Optimization," in MURRAY, W. (ed.), Academic Press, London, pp. 13-28. LADE,P.V. and DUNCAN, J.M., "Elastoplastic Stress-Strain Theory for Cohesionless Soil," J. Geotech. Eng. Div., ASCE, 101, No. GTI0, 1037. ZIENKIEWlCZ,O.C., H ~ r m S O N , C., and LEwis, R.W., "Associated and Nonassociated Viscoplasticity and Plasticity in Soil Mechanics," Geotechnique, 25, No. 4, 671. SANDLER,I.S., DIMAC~IO, EL., and BALADI, G.Y., "Generalized Cap Model for Geological Materials," J. Geotech. Div., ASCE, 102, No. GT7, 683. LADE,P.V., "Elastic-Plastic Stress-Strain Theory for Cohesioniess Soil with Curved Yield Surfaces," International Journal of Solids and Structures, 13, 1019. FLETCHER,R., "Practical Methods of Optimization-Volume l: Unconstrained Optimization," John Wiley and Sons, Chichester, U.K. LOUNG,M.P., "Stress-Strain Aspects of Cohesionless Soils Under Cyclic and Transient Loading," in Proc. Int. Symposium on Soils Under Cyclic and Transient Loading, PANDE, G.N. and ZIENKIEWICZ, O.C. (eds.), A.A. Balkema, Rotterdam, Vol. 1, 315. FLETCHER,R., "Practical Methods of Optimization-Volume 2: Constrained Optimization," John Wiley and Sons, Chichester, U.K. GHA~OUSSl,J. and MOMEN, H., "Modelling and Analysis of Cyclic Behavior of Sand," in PANDE, G.N. and ZmNKmWlCZ, O.C. (eds.), Soil Mechanics-Transient and Cyclic Loads, Wiley, p. 313. DArALIAS,Y.E and HARRUANN, L.R., "Bounding Surface Formulation of Soil Plasticity," in PANDE, G.N. and ZIENKIEWlCZ, O.C. (eds.), Soil Mechanics-Transient and Cyclic Loads, p. 253. VALAmS,K.C. and READ, H.E., "A New Endochronic Plasticity Model for Soils," in PANDE, G.N. and ZmNKIEWlCZ, O.C. (eds.), Soil Mechanics--Transient and Cyclic Loads, p. 375. DENATALE,J.S., "On the Calibration of Constitutive Models by Multivariate Optimization. A Case Study: The Bounding Surface Plasticity Model," Ph.D. thesis, University of California, Davis. RELAKITIS,G.V., RAVINDRAN,A., and RAGSDELL,K.M., "Engineering Optimization: Methods and Applications," John Wiley and Sons, New York. DESAI,C.S. and FARUQUE, i . O . , "Constitutive Model for (Geological) Materials," Journal of Engineering Mechanics, ASCE, 110, 1391. FUKUSHIMA,S. and TATSUOKA,F., "Strength and Deformation Characteristics of Saturated Sand at Extremely Low Pressures," Soils and Foundations, 24, No. 4, 30. FARUQUE,M.O. and DESAI, C.S., "Implementation of a General Constitutive Model for Geological Materials," Int. J. Num. Analyt. Meth. Geomech., 9, 415. POOROOSHASB,H.B. and PmTRUSZCZAK,S., "On Yielding and Flow of Sand; A Generalized Two Surface Model," Computer and Geotechnics, 1, 33. DESAI,C.S., SOMASUNDARAM,S., and FRANTZISKONIS, G., "A Hierarchical Approach for Constitutive Modelling of Geologic Materials," Int. J. Num. Analyt. Meth. Geomech., 10, 3, 225. FARUQUE,M.O., "A Third Invariant Dependent Cap Model for Geological Materials," Soils and Foundations, 27, No. 2, 12. HIRAI,H., "A Combined Hardening Model for Plasticity for Sands," 2nd Int. Conf. Short Course on Constitutive Laws for Engineering Materials: Theory and Applications, DESAI, C.S. et al, (eds.), Elsevier, New York, 557. CmSTESCU,N., "Rock Rheology," Kluwer Academic Publishers, The Netherlands. HARDIN,B.O., "Low-Stress Dilation Test," J. Geotech. Eng., ASCE, 115, No. 6, 769. ABDLrLIA~AD,S.N., FARUQUE,M.O., and AZEEMUDDIN,M., "Experimental Verification and Numerical Evaluation of the Three-Invariant Dependent Cap Model for Cohesionless Soil," The Arabian Journal for Science and Engineering, 16, 107.
School of Civil Engineering and Environmental Science The University of Oklahoma, Norman Norman, OK 73017, USA (Received in final revised form 14 June 1993 )