Elastoplastic Model for Clay with Microstructural Consideration

Clay material can be considered as a collection of clusters, which interact with each other mainly through mechanical forces. From this point of view, clay is modeled by analogy to granular material in this paper. An elastoplastic stress-strain relationship for clay is derived by using the granular mechanics approach developed in previous studies for sand. However, unlike sand, clay deformation is generated not only by the mobilizing but also by compressing clusters. Thus, in addition to the Mohr-Coulomb’s plastic shear sliding and a dilatancy type flow rule, a plastic normal deformation has been modeled for two clusters in compression. The overall stress-strain relationship can then be obtained from the mobilization and compressing of clusters through a static hypothesis of the macro-micro relations. The predictions are compared with the experimental results for clay under both drained and undrained triaxial loading conditions. Three different types of clay, including remolded and natural clay, have been selected to evaluate the model’s performance. The comparisons verify that this model is capable of accurately reproducing the overall behavior of clay, which accounts for the influence of key parameters such as void ratio and mean stress. A section of this paper is devoted to show the model’s capability of considering the influence of inherent anisotropy on the stress-strain response under undrained triaxial loading conditions. keywords: Clays; Microstructures; Stress strain relations; Anisotropy; Elastoplasticity.


Introduction
Under a microstructural approach for stress-strain modeling, a set of microsystems ͑e.g., interparticle contacts͒ is used to represent the material.Material properties are defined for each microsystem and the overall stress-strain relationship for the material is obtained from averaging the behaviors of this set of microsystems.
If the microsystems are considered as a set of mobilized planes in the material, the approach used to estimate the overall behavior from this set of planes can be linked to G. I. Taylor's concept, developed long ago, as the slip theory of plasticity for polycrystalline materials by Batdorf and Budianski ͑1949͒.These ideas were applied by Pande and Sharma ͑1982͒ to rocks and soils in what they called the overlay model, and to concrete by Bazant et al. ͑1996͒ in the so-called microplane model.
For granular material, a set of particle pairs in contact are considered as the microsystems, and material constants are defined for the interparticle contacts.The approach used to estimate the overall stress-strain behavior of granular material can be found in the work done over the last two decades.For elastic behavior, work has been done by Jenkins ͑1988͒, Walton ͑1987͒, Rothenburg and Selvadurai ͑1981͒, Chang ͑1988͒, Emeriault and Cambou ͑1996͒, Liao et al. ͑2000͒, and Kruyt and Rothenburg ͑2002͒, among others.For inelastic stress-strain behavior, models can be found in Jenkins and Strack ͑1993͒, Matsuoka and Takeda ͑1980͒, Chang et al. ͑1989a,b͒, etc. Chang et al. ͑1992a,b͒, and Suiker and Chang ͑2004͒.These inelastic models have encountered some difficulties in predicting the correct shear strength in so far as the predicted strength is often overestimated for the models with a kinematic hypothesis ͑Chang and Misra 1990͒.These models also have difficulties in predicting correctly the behavior under different stress paths.To overcome these difficulties, we have adopted a static hypothesis proposed by Liao et al. ͑1997͒ and incorporated the critical state concept to formulate a microstructural model suitable for sand ͑Chang and Hicher 2005͒.
In this paper, we treat clay as a collection of clusters formed by groups of platy clay particles.At the scale of cluster sizes, long range forces are negligible, and the clusters interact with each other mainly through mechanical forces.Thus, clay material can be modeled by analogy to granular material.Some studies on clay fabric supporting this notion are discussed in the "Simple Model for Clay Fabric" section.A stress-strain model is then presented in the "Stress-Strain Model" section, in which we consider a clay cluster as a deformable grain.We then extend the granular mechanics approach ͑i.e., Chang and Liao 1990;Chang and Gao 1995;Chang and Hicher 2005͒ to derive the elastoplastic stressstrain relationship for clay.
The model's performance is then evaluated by comparing the predicted with the measured triaxial loading results for clay specimens of various overconsolidation ratios ͑OCRs͒, under various confining stresses, and in both drained and undrained conditions.The experimental results are obtained from two types of remolded clay with very different compressibilities.A natural clay with inherent anisotropy has also been studied to evaluate the model's ability to model the effects of inherent anisotropy on the stressstrain response under undrained triaxial loading condition.

Simple Model for Clay Fabric
The constituents of clay can be generally viewed across three scales 1. Clay particle: a clay particle is usually platy in shape.Its thickness and size can vary 100 times according to specific clay types, such as montmorillonite, illite, or kaolinite.The size for a platy particle generally ranges from 0.01-1 m. 2. Clay cluster as an aggregate of clay particles: clay particles attract each other due to surface forces among particles such as chemical, electrostatic, van der Waals forces, etc.These forces pull together the particles to form particle clusters, which have either a flocculate or dispersed type structure, as shown in Fig. 1.The size of the clusters continues to grow until the clusters are large enough so that the cluster weight, due to gravitation, becomes significantly larger than the interparticle surface forces.At this stage, the cluster loses its potential to attract further clay particles, and the size of clusters stops to grow.The ultimate cluster-size depends on the clay particle type, the liquid inside the pores, and its sedimentation history.3. Clay material as an assembly of clusters.
Since the clay-particles are strongly attracted to each other by surface forces, a cluster of particles does not deform very much under usual external stresses and can, therefore, be considered as an intact unit.On the other hand, clusters interact mainly through mechanical forces.Thus, as schematically shown on the left side of Fig. 2, clay can be regarded as an assembly of "grains," in which each grain is also a cluster.The main difference between sand grains and clay clusters is that the clay clusters, compared to sand grains, are more deformable.This deformability depends on the way the clay particles are assembled, which is a function of the mineralogy, adsorbed ions, etc.
A study undertaken by Hicher et al. ͑2000͒, by means of scanning and transmission electron microscopes, related the mechanical behavior of two saturated remolded clays, a kaolinite and a bentonite, to the evolution of their microstructural parameters such as particle shape, particle size, and particle orientation.The results of this study show that the main role played by the clay clusters is similar to the role played by the grains in the mechanical behavior of granular materials.This explains why sand and clay have similar qualitative behavior even though each material consists of different constituents ͑Biarez and Hicher 1994͒.The difference in nature between sand grains and clay clusters can nevertheless explain quantitative differences in their stress-strain relationship.In particular, the deformability of the clusters can have a significant effect on volume change behavior in drained condition and on the shape of stress path in undrained condition.The influence of the cluster's deformability on the elastic properties of clayey materials is particularly pronounced.Since the elastic domain is restricted to very small strains ͑Ͻ10 −5 ͒, for which the relative displacements of the clusters are negligible, the deformation of the material must be mainly due to the deformation of individual clusters.Thus, the elastic moduli measured in sands and gravels are much higher than those in clays ͑Biarez and Hicher 2008͒.
At higher strain level, the major deformation mechanism of the continuous medium is due to the relative displacements of the clay clusters.Under theses conditions, the main factor affecting the overall behavior is the friction between clusters, determined by the nature and the shape of these clusters.The photographs in Fig. 3͑a͒ show how kaolinite microstructure is made up of rigid small particles and aggregated in compact clusters.Fig. 3͑b͒ shows the clusters of both natural and remolded Saint Herblain clay.
An analysis of the pore size distribution in the kaolinite, by means of mercury intrusion porosimetry, confirmed the existence of two major groups of pore sizes: one centered around 1 m and the other around 10 nm.The first group corresponds to the intercluster pores, whereas the second one is related to the interparticle pores.The number and the size of large pores progressively decreased wit increasing consolidation stresses, while the small pore size remained relatively unchanged under moderate loading stresses.This result confirmed that the volume change during loading is mainly due to the rearrangement of the clayclusters ͑Hicher et al. 2000͒.

Stress-Strain Model
In this model, clay is regarded as an aggregate of clusters.The deformation of a representative volume of the material is generated by mobilizing and compressing the clusters.Thus, the stressstrain relationship can be derived as an average of the deformation behavior of local contact planes in all orientations.
For contact planes in the ␣th orientation, the local forces f j ␣ and the local movements ␦ i ␣ are denoted by , where the subscripts n, s, and t represent the components in the three directions of the local coordinate system, as shown in Fig. 4. The direction outward normal to the plane is denoted as n; the other two orthogonal directions, s and t, are tangential to the plane.

Density State of a Cluster Assembly
One of the important elements to consider in soil modeling is the critical state concept.At critical state, the clay material remains at a constant volume while it is subjected to a continuous distortion.The void ratio corresponding to this state is e c .
The critical void ratio e c is a function of the mean stress p = ͑ x + y + z ͒ / 3. The relationship has traditionally been written as follows: The two parameters ͑e cr0 , p cr0 ͒ represent a reference point on the critical state line.For convenience, the value of p cr0 is taken to be 0.01 MPa.The critical state line can be defined by two parameters e cr0 and .Using the critical state concept, the density state of an assembly is defined as the ratio e / e c , where e is the void ratio of the assembly.

Intercluster Behavior
In order to have a more apparent link between the micro and macro variables, we define a local stress i ␣ and a local strain ␥ i ␣ , which are directly related to the local forces f j ␣ and the local movements ␦ i ␣ at each contact by where l ␣ = length of the branch vector, which joins the centroids of two contact clusters, and V = volume of the representative element.It is to be noted that the local stress i ␣ is not the stress on the physical contact area between two clusters.It should be viewed rather as the average stress on the intercluster plane when the cluster and voids in the representative volume are homogenized into a continuum.For an isotropic medium, the local stress is identical to the tractions resolved on the plane due to global stress ͑i.e., i ␣ = ji n j ␣ ͒.The proof will be given later in Eq. ͑19͒.In the local coordinate system, the local stress and local strain are respectively denoted as For convenience, we use the notation ␣ = n ␣ for local normal stress and the notation ␣ = ␥ n ␣ for local normal strain in the following sections.

Elastic Part
The intercluster behavior can be characterized as the relationship between contact force and contact displacement, given by in which the stiffness tensor can be related to the contact normal stiffness, k n ␣ , and shear stiffness, k r In terms of local stress and local strain, an alternative intercluster stiffness is defined as where The intercluster stiffness can be expressed as the form adopted for sand grains by Chang et al. ͑1989a,b͒, given by where ␣ = local stress in a normal direction; p ref = standard reference pressure taken as 1 MPa.; and k rR = ratio of shear to normal stiffness.k ¯n0 ␣ , k rR and n are material constants.The value of n is found to be 0.33 for two elastic spheres according to Hertz-Mindlin's formulation ͑Mindlin 1969͒.Based on experimental measurements of the elastic modulus under different confining stress, the value of n has been found to be 0.5 for sand and 0.5-1.0 for clay.

Shear Sliding
The elastic part of the tangential movement between two clusters does not have a coupling effect ͑i.e., there is no shear induced normal movements͒.However, plastic sliding often occurs along the tangential direction of the contact plane with an upward or downward movement ͑i.e., dilation or contraction͒.Stress dilatancy is a well-known phenomenon in sand ͓see discussions in the work by Taylor ͑1948͒, Rowe ͑1962͒, Goddard ͑1990͒, etc.͔, and should be correctly modeled.The dilatancy equation used here is modified from the equation adopted for sand by Chang and Hicher ͑2005͒, given by The modified equation allows more flexibility in modeling performance for different behaviors.In this equation, a, b, and are intercluster property constants, e c is the critical void ratio for the clay, and ͑sgn͒ is the sign of ͑ / − tan ͒.When the void ratio e is equal to the critical void ratio, zero dilation holds.It is to be noted that the state variables e and e c of the clay are referred to the cluster assembly, which is used to regulate the dilation of individual intercluster contacts.It is rational to consider the micro variable as a function of the macrostate because the intercluster behavior is indeed influenced by the density state of the specimen.
In Eq. ͑7͒, is the intercluster friction angle, which, in value, is very close to the internal friction angle, measured at critical state.The values of a and b can be calibrated from experimental measurements of triaxial tests, which will be shown in the later section on numerical simulation.
Note that the shear stress and the rate of the plastic shear strain d␥ p in Eq. ͑7͒ are defined as The yield function is assumed to be of Mohr-Coulomb type, given by where 1 ͑␥ P ͒ = isotropic hardening/softening parameter.Plastic loading corresponds to dF 1 Ͼ 0. The hardening parameter is defined by a hyperbolic function in the 1 -␥ p plane, which involves two material constants: p and k ¯p through the following relationship: When plastic deformation increases, 1 asymptotically approaches tan p .For a given value of , the initial slope of the hyperbolic curve is k ¯p / .The flow rule is nonassociated.Under a loading condition, the shear plastic flow in the direction tangential to the contact plane is determined by a normality rule applied to the yield function.However, the plastic flow in the direction normal to the contact plane is governed by the stress-dilatancy equation in Eq. ͑7͒.
The value of k ¯p is found to be linearly proportional to k ¯n such that The ratio k pR is a material parameter.
The internal friction angle is a constant for a given material.However, the apparent friction angle, p , on a contact plane is dependent on the density state of neighboring clusters, which can be related to the void ratio e by tan p = ͩ e c e ͪ m tan ͑13͒ where m = material constant ͑Biarez and Hicher 1994͒.
In a loose structure, clusters can rotate more freely, preventing the intercluster shear force from fully mobilizing the sliding resistance.The apparent frictional angle p is smaller than .On the other hand, a dense structure has a higher degree of interlocking, which requires more effort to mobilize the contacts between clusters.In such a case, the apparent frictional angle p is greater than .When the loading stress reaches the apparent frictional angle p , the dense structure dilates and the degree of interlocking becomes relaxed.As a consequence, the apparent frictional angle ͑i.e., the peak angle in this case͒ is reduced, which results in a strain-softening phenomenon.

Normal Compression
To describe the compressible behavior between two clay clusters, a second yield surface is necessary.The second yield function is assumed to be as follows: where the local normal stress and local normal strain p are defined in Eq. ͑3͒.The hardening function 2 ͑ p ͒ is defined as where c p = compression index for the compression curve plotted on the p − log plane.When is less than p , the plastic strain produced by the second yield function is null.Thus, p in Eq. ͑12͒ corresponds to the preconsolidation stress in soil mechanics.

Elastoplastic Relationship
With the basic elements of intercluster behavior ͑elastic and plas-tic͒ discussed above, the final incremental local stress-strain relation of the intercluster contact can be derived Since it is a standard procedure to derive a detailed expression for the elastoplastic stiffness tensor k ¯ij ␣p by using the plastic yield functions given above, the derivation is not presented here.

Stress-Strain Relationship
The stress-strain equations at cluster scale represent the relationships between two vectors-intercluster stress i and intercluster strain ␥ j -whereas the stress-strain equations at the assembly scale represent the relationships between two tensors-stress ij and strain kl .Due to the two different levels of complexity, it is much easier to establish, from a phenomenological approach, the stress-strain equations at cluster scale than those at assembly scale.Thus, in the present model, the stress-strain equations at cluster scale are first established in a phenomenological manner, as shown in the last section.In this section, these equations are used to derive the stress-strain equations at assembly scale.This way is more rational and takes less effort than a straight phenomenological approach of establishing the stress-strain relationship at assembly scale.

Macro-Micro Relationship
The stress-strain relationship for an assembly of clay clusters can be determined from integrating the intercluster behavior at all contacts.During the integration process, a relationship is required to link the macro and micro variables.Using the static hypothesis proposed by Liao et al. ͑1997͒, we obtain the relation between the strain of assembly and intercluster strain where ␥ ˙j = local strain between two contact clusters; n k = unit vector of the branch joining the centers of two contact clusters; and N = total number of contacts, over which the summation is carried out.The tensor B ik ␣ in Eq. ͑16͒ is defined as Using the principle of energy balance, which states that the work done in a representative volume element is equal to the work done on all intercluster planes within the element, we have and using Eq.͑16͒, the local stress on the ␣th contact plane is derived as follows: For the case of isotropic fabric, it can be derived that B ik =3␦ ik / N, where ␦ ik is the Kronecker delta.Thus Eq. ͑19͒ is reduced to the usual form ˙j ␣ = ˙ij n j ␣ .The stress increment ˙ij can be obtained by adding the dyadic product of the contact force and the branch vectors for all contacts ͑Christofferson et al. 1981; Rothenburg and Selvadurai 1981͒.In terms of local stress, it is Applying the defined local stress in Eq. ͑19͒, Eq. ͑20͒ is unconditionally satisfied.By using Eqs.͑15͒, ͑16͒, and ͑19͒, the following relationship between stress and strain can be obtained: where The summation in Eq. ͑22͒ can be expressed by a closed-form solution for some limited conditions such as the elastic modulus of randomly packed equal-size particles ͑Chang et al. 1995͒.However, for elastic-plastic behavior, due to the nonlinear nature of the local constitutive equation, a numerical calculation with an iterative process is needed to carry out the summation in Eq. ͑22͒.

Computation Scheme
Initially, we know the global variables ͑ ij and ij ͒ for the assembly and the local variables ͑f j ␣ and ␦ j ␣ ͒ for each contact.For a given loading increment, which can be of a stress control, a strain control or a mixed mode, 6 out of the 12 variables ͑⌬ ij and ⌬ ij ͒ are unknown.The objective is to determine all global ͑ ij and ij ͒ and local variables ͑f j ␣ and ␦ j ␣ ͒ at the end of the load increment.For a system with N intercluster contacts, the number of unknowns is 3N for f j ␣ and 3N for ␦ j ␣ .The total number of unknowns is therefore 3N +3N +6.
The following constraints must be satisfied: 1.The local constitutive equation, i.e., Eq. ͑15͒.Since there are three equations for each contact, the total number of equations is 3N; N being the total number of intercluster contacts.2. The static hypothesis between global stress and local forces, i.e., Eq. ͑19͒: the number of equations is 3N. 3. The strain definition between global strain and local displacement, i.e., Eq. ͑16͒.The number of equations is 6 ͑strain tensor is symmetric͒.The total number of unknowns is the same as the total number of equations.Therefore, a solution can be determined.To facilitate the numerical calculation, the summation process in the above equations can be replaced by an integral process in a spherical coordinate system with an orientation distribution function E͑␣ , ␤͒ for the intercluster contacts, provided that the number of contact N is sufficiently large ͑Chang and Misra 1990͒.An example is the fabric tensor in Eq. ͑17͒ The surface integration on a sphere can be carried out numerically through Gaussian integration points with weight factors.Thus the integration process requires less effort than the summation process because it can be integrated over a small number of Gaussian integration points ͑i.e., selected number of contact ori-entations͒.Lebedev ͑1976, 1977͒ has pioneered an integration scheme using a set of integration points with octahedral symmetry.Bazant and Oh ͑1986͒ has devised an integration method similar to the ones proposed by Lebedev.Along this line, a more refined method has been proposed by Delley ͑1996͒ for applications in weather forecasting, quantum chemistry, wave scattering, and radiation studies.
For all simulations presented in this paper, we have performed using three numbers of integration points N = 56, N = 74, and N = 122.The predicted stress-strain curves show about 5-7% difference between the results obtained from integration numbers N = 56 and 74, whereas the curves show less than 1% difference between the results obtained from integration numbers N = 74 and 122.Thus N = 74 was found to be adequate.However, it is antici-pated that in the case of strong anisotropy, the number of orientations may require higher values to reach convergence.
For a strain-controlled test, Eq. ͑21͒ presents numerical difficulties at the after-peak range with strain softening.In this case, a method of "elastic predictor-plastic corrector" was adopted to obtain the solution.The most widely used method for the solution of nonlinear constitutive equations is in the category of elastic predictor-plastic corrector method ͑Ortiz et al. 1983͒, where a purely elastic trial state is followed by a plastic corrector stage ͑return mapping algorithm͒.The purpose of the return mapping is to enforce consistency, at the end of the load step, of the prescribed yield surface and flow rule.For simple classical plasticity models, the return path can be determined in closed form.However, for the present model that accounts for pressure sensitive nonlinear work hardening/softening, nonlinear elasticity, and two yield surfaces, it becomes necessary to compute the return path in an iterative fashion.An implicit integration scheme is more stable to obtain the solution of the system of nonlinear equations through iterative processes.The specific method adopted for computation is the single step backward Euler return method ͓see Simo and Hughes ͑1998͔͒.
It is to be noted that the implicit integration scheme with the single step backward Euler return method is applicable only for a strain controlled test.For a mix-mode loading condition ͑i.e., for some components, stresses are specified rather than strains such as that in triaxial compression test conditions͒, an additional iteration process is also needed to satisfy the condition of specified stresses.Since the procedure is straightforward, thus it is not included here.

Summary of Parameters
The material parameters are summarized as follows: 1. Microstructural descriptions ͑two parameters͒. a.
Contact number per unit volume, N / V and mean cluster size, d.Here we use this equation as a first-order approximation to estimate N / V for clay by treating d as the mean size of the clay clusters.It is to be noted that the value of contact number per unit volume changes with void ratio.The evolution during the deformation process is taken into account.
The intercluster parameters are not feasible to be determined from direct measurements on interclusters due to experimental difficulties.A possible way of parameter determination might resort to the numerically simulated cluster behavior by the discrete element method.However, this approach can only be applied after the discrete element simulation is fully verified by experiments.Thus, for convenience, the interclusters parameters in the present model are to be phenomenologically calibrated from the behavior of soil sample measured in conventional laboratory tests.
Among the interclusters parameters, the exponent n is generally between 0.5 and 1.0 for clay, and the typical value of exponent m is 1.The typical value is between 0.25 and 1 for the ratio k pR , and is about 0.5 for the ratio k rR .The other parameters can be easily obtained from standard laboratory experiments.
From an isotropic compression test, four parameters can be determined, namely, e 0 , , k ¯n0 ␣ , and c p .The void ratio e 0 and can be measured directly from the compression line.The contact stiffness k ¯n0 ␣ can be calculated by using the rebound index measured from the rebound curve where the factor is calculated given the pressure range p 1 to p 2 , from which the rebound index is estimated.
Using the compression index and the rebound index of an assembly of clay clusters, we can also determine the plastic intercluster compression index c p from the following equation: The other parameters , k pR , k rR , a , b, and e cr0 can be obtained from two drained triaxial tests.An example will be given in the following section treating parameter calibration for white kaolinite clay.

Results of Numerical Simulations on Clays
In order to evaluate the model's performance, the predicted results are compared with experimental measurements from triaxial compression tests under both drained and undrained conditions for white kaolinite clay and black kaolinite clay.The predicted intercluster behavior are illustrated to show the link between macro and micro phenomena.In addition, a prediction

White Kaolinite Clay
White kaolinite clay is a remolded clay prepared in the laboratory from a mixture of dry clay powder and water.The slurry is then progressively consolidated.The white clay has a plastic limit of about 30% and liquid limit of about 60%.According to scanning electron microscope photos, we assume that the mean cluster size d is 0.004 mm and the value of N / V is calculated from Eq. ͑23͒.

Calibration of Model Parameters
Calibration of the parameters for white kaolinite clay is illustrated in this section.The calibration requires one isotropic compression test and two drained triaxial compression tests under different confining pressures.The experimental results used for calibration are shown in Fig. 5.
1. Intercluster elastic constants: k ¯n0 , k rR and n-the exponent n can be determined from stress-strain curves considering very small strains.As reported by Hicher ͑2001͒, the value is 0.75 for white kaolinite clay.The value of k ¯n0 can be determined from Eq. ͑24͒, and its calibration is shown in Fig. 6͑a͒.k rR can be determined from the e-1 curve of a drained compression test at a small strain condition, as shown in Fig. 6͑b͒, in which the experimental data are taken from Fig. 5͑c͒.2. Intercluster friction angle: and m-the intercluster friction angle is the slope of the critical state line on the p-q plane, as shown in Fig. 5͑b͒.A typical value of 1 is used for m. 3. Intercluster normal hardening rule: c p and k pR -the value of c p can be determined from Eq. ͑26͒ based on the slopes of  The value of e 0 , e cr0 and can be determined from the isotropic compression line and the critical state line on the e-ln p plane.The model parameters determined for white kaolinite clay are listed in Table 1.

Drained Tests on Normally Consolidated Clay
Drained triaxial tests on normally consolidated white kaolinite clay have been reported and analyzed by Biarez and Hicher ͑1994͒.Two specimens isotropically consolidated up to respectively 0.4 and 0.8 MPa, were loaded to failure in drained condition.Using the parameters in Table 1, the predicted test results are plotted in Fig. 5.The stress-strain curves in Fig. 5͑a͒ and the void ratio change in Figs.5͑c and d͒ show good agreement between experimental and computed curves.The paths in the e-log pЈ plane ͓Fig.5͑d͔͒ show that the void ratio approaches the critical state line when the stress state approaches the failure line in the p-q plane.Since this set of experimental data are used for parameter calibration, a good comparison is expected.

Undrained Triaxial Tests on Normally and Overconsolidated Clay
Experimental data of undrained triaxial tests on normally and overconsolidated white kaolinite clay specimens were also reported by Biarez and Hicher ͑1994͒.Three samples were isotropically consolidated up to 0.8 MPa, and two of them were unloaded to 0.4 and 0.067 MPa, so that the OCRs of the three samples were equal to 1, 2, and 12, respectively.The same material parameters in Table 1 were used to predict the stress-strain behavior for the three undrained tests.
As shown in Fig. 7͑a͒, the computed and measured stressstrain curves are in good agreement.The stress paths shown in Fig. 7͑b͒ indicate that for the normally consolidated ͑OCR= 1͒ and slightly overconsolidated ͑OCR= 2͒ samples, the stress paths do not overpass the critical state line.While, for the strongly overconsolidated specimen ͑OCR= 12͒, the stress path goes above the critical state line, at which dilation occurs, leading to an increase of the mean effective stress.The pore pressure development and e-log pЈ curves are shown in Figs.7͑c and d͒.The comparisons show good agreement.
It is to be noted that parameters a and b have significant influence on the stress path of undrained test.Fig. 8 shows different effective stress paths as influenced by different values of a and b.Based on the comparison in Fig. 7, the model using parameters calibrated from drained tests is capable of predicting the stressstrain behavior of undrained tests for specimens under different overconsolidation ratios.

Black Kaolinite Clay
Black kaolinite clay is also a remolded clay mixed from clay powder with a darker color.The black kaolinite clay has a plastic limit p =30%, and liquid limit l =70%, which is much more compressible than the white kaolinite clay.The parameters presented in Table 1 are calibrated from triaxial compression tests on normally consolidated specimen, using the method described above.Again, the mean cluster size d is 0.004 mm and the value of N / V is calculated by Eq. ͑23͒.

Drained Tests on Overconsolidated Clay
Tests on black kaolinite clay samples were performed by Zervoyanis ͑1982͒ and analyzed by Biarez and Hicher ͑1994͒.The tests begin with an isotropic consolidation up to 0.8 MPa, then unloaded to 0.4, 0.2, and 0.1 MPa, respectively.The OCRs are 1, 2, 4, and 8, respectively.The predicted void ratio changes in Fig. 9͑c͒ show a contractive behavior for OCR= 1 and 2, and a dilative behavior for OCR= 4 and 8.The contractive and dilative behaviors can also be seen in the predicted paths on the e-log p curves in Fig. 9͑d͒.For OCR= 4 and 8, the stress-strain curves in Fig. 9͑a͒ show strain softening, which corresponds to the stress paths in Fig. 9͑b͒ above the critical state line.An overall good agreement is observed between experimental and predicted results for different OCRs.The two examples of black and white kaolinite clay demonstrate that the model is capable of reproducing with good precision the mechanical behavior of both stiff and soft remolded clays under drained and undrained triaxial test conditions.

Microplane Behavior
The model is capable of describing the intercluster strains for cluster contacts in all orientations.The orientation of a given contact plane is represented by an angle ␥ measured from the x axis to the branch vector, as shown in Fig. 4. The angles ␥ selected in this study are 18°, 45°, 55°, and 72°respectively ͑note that the vertical orientation, ␥ = 0, corresponds to a horizontal contact plane͒.The deformation behavior of the intercluster contact planes are discussed here for both drained and undrained conditions.

Drained Conditions
The drained test on black kaolinite clay with OCR= 1 and 3c = 0.8 MPa ͑see Fig. 9͒ is selected for showing intercluster behavior.Four steps are selected for this test ͓see circles marked on Fig. 9͑a͔͒.For each step, the intercluster stress and strain on contacts of various orientations are plotted to show their evolution.Fig. 10 shows the local stress-strain relationships for the contact planes in the selected orientations.From thecurve ͓Fig.10͑a͔͒, we can see that the local stress paths have different slopes for plane orientations from 18°to 72°.The four steps are also marked on each stress path.Under an increase of the vertical stress, the planes oriented near the horizontal direction ͑i.e., small values of ␥͒ are subjected mainly to a normal stress component ⌬.The shear component becomes more significant when the plane is inclined.
The local shear stress-strain curves ͓Fig.10͑b͔͒ show that every plane is mobilized to a different degree.The actively moving planes are located in a narrowly oriented band near the orientation of about 55°, which is responsible for the overall deformation of the specimen.Other planes are inactive with small movement.This clearly indicates that the local strains do not uniformly conform to the overall strain of the specimen.
The local stresses plotted in Figs.10͑c and d͒ indicate that the local stresses uniformly conform to the overall stress of the assembly, which proves that the static hypothesis ͑i.e., j ␣ = ij n i ␣ ͒ used in this model is satisfied in all load steps.The local shear strains plotted in Fig. 10͑e͒ indicate that the local shear strains are relatively uniform up to the load Step 2 and become highly nonuniform at Step 3 and Step 4. The largest shear strain occurs on the planes near the orientation at about 55°.Fig. 10͑f͒ shows a very small change of the local normal strains during the four load steps.However, at load Step 3 and Step 4, the plot illustrates nonuniform strains, where contraction occurs accompanying the The drained test on black kaolinite clay with OCR= 8 and 3c = 0.1 MPa ͑see Fig. 9͒ is also selected for showing the intercluster behavior.Four load steps are selected for this test ͓see circles marked on Fig. 9͑a͔͒ to plot intercluster stresses and strains to show their evolution.
The trends of our results are similar to that of OCR= 1.The major difference is that local stresses on some planes can exceed the line of tan due to the interlocking of aggregates in this overconsolidated specimen.Fig. 11͑b͒ shows that the plane of the orientation equal to 55°contributes mostly to the deformation of the assembly.The local stress plotted in Figs.11͑c and d͒ indicates that the local stresses uniformly conform to the overall stress of the assembly at all steps.The local shear and normal strains plotted in Figs.11͑e and f͒ indicate that the largest shear strain, occurring on the planes near the orientation equal to 55°, is accompanied by large dilation.Thus a dilation shear band may occur at failure for the overconsolidated clay.

Undrained Conditions
The undrained test on white kaolinite clay with OCR= 1 and 3c = 0.8 MPa ͑see Fig. 7͒ has been selected to demonstrate intercluster behavior.Four load steps are also selected for this test ͓see circles marked on Fig. 7͑a͔͒ to show the evolution of local stresses and strains.
From thecurve in Fig. 12͑a͒, we can see different local stress paths for undrained loading conditions.The local shear stress-strain curves ͓Fig.12͑b͔͒ show that every plane is mobilized to a different degree.The active planes are located in a narrowly oriented band near the orientation of 55°.
The local stresses plotted in Figs.12͑c and d͒ indicate that the local stresses uniformly conform to the overall stress of the assembly, which demonstrates that the static hypothesis ͑i.e., j ␣ = ij n i ␣ ͒ is satisfied in all load steps.The local shear strains plotted in Figs.12͑e and f͒ indicate that the local shear strains are relatively uniform up to the load Step 3, and become highly nonuniform at Step 4. The largest shear strain occurs on the planes near the orientation of 55°, but no specific volume change has been associated with it, because of the constraint of overall zero volume change due to undrained conditions.

Inherent Anisotropy
Natural clay often exhibits anisotropic behavior.Kirkgard ͑1991͒ performed undrained triaxial test on samples of San Francisco Bay mud to study the degree of anisotropy in this clay deposit.The soil samples were taken from a depth of 6.1 m with a water table depth of 6 m.Gradation analysis showed that the soil consists of 55% clay and 45% silt.The soil has the following properties: water content of 98.5%, unit weight of 14.4 kN/ m 3 , specific gravity of 2.55, liquid limit of 1.36, and initial void ratio of 1.35.The samples were tested by applying an axial load in the horizontal and vertical directions to investigate different behaviors due to material anisotropy.Material anisotropy is characterized by local material constants, which are orientational dependent, and can be expressed as a function of ␤ and ␥ ͑the two angles in the spherical coordinate system, as shown in Fig. 4͒.To describe such a parameter, a density function E͑␤ , ␥͒ has been introduced.Integration of E͑␤ , ␥͒ over all orientations should be equal to 1, i.e.
Thus, for an isotropic material, the function E͑␤ , ␥͒ =1/ 2. For an orthotropic material, the density function can be expanded to a series using the method of spherical harmonic expansion in three dimensions.The truncated form of the series consisting of only second-order terms is For a cross-anisotropic material, a 22 becomes equal to zero and Eq.͑28͒ is reduced to In three dimensions, the inherent anisotropy can be represented by a distribution whose major axis often coincides with vertical or horizontal directions.An example of distributions with different values of a 0 is shown in Fig. 13.The axes of anisotropy of the soil are identical to those of the axes of loading stresses.
It has been found ͓see Chang and Misra ͑1990͔͒ that the anisotropy can be characterized by a tensor that matches the coeffi- where ave = average value.In the case of ave = 31.2°anda o = −0.25, the friction angle in the minor axis is 1 = 23.4°and in the major axes 2 = 3 = 35.1°.For a soil layer with an inherent anisotropy due to the geological deposition process, the material properties are usually cross anisotropic, with a symmetry around its major axis that coincides often with the vertical direction.Two samples cored from vertical and horizontal directions are represented schematically as a cube in Fig. 14.The shaded area is perpendicular to the vertical direction ͑Direction 1͒.The properties in Directions 2 and 3 are the same, but different from the ones in Direction 1.Using the averaged behavior of vertical and horizontal samples, we can calibrate the parameters which are summarized in Table 1.
For the purpose of comparing the predictions with the experimental results, two triaxial compression tests were simulated by the model.The two soil samples were isotropically consolidated up to 0.125 and 0.175 MPa, respectively.Afterward, both samples were deviatorically loaded in undrained condition.To model the difference in shear strength for vertical and horizontal samples, we assign the following anisotropy of the friction angle: 3 = 35.1°and 1 = 25.4°.The predicted curves in Fig. 15͑a͒ show good agreement with the maximum strength of the two samples.To improve the predicted initial slopes of the undrained stress paths, an additional anisotropy is assigned to the elastic stiffness, ͑k ¯n0 ͒ 1 = 600 and ͑k ¯n0 ͒ 3 = 150.As a result, the predicted undrained stress paths in Fig. 15͑d͒ show more difference in the initial slopes between the vertical and horizontal samples, thus showing a better comparison with the experimental results.

Summary and Conclusion
The notion of treating clay clusters as grains makes it promising for extending the well tested methodology of modeling granular material to the modeling of clay.In the newly developed microstructural model for clay, the overall strain includes plastic sliding and plastic compression among clay clusters.Although the model involves the use of parameters at the scale of clay cluster, we do not need to obtain these parameters directly from experimental tests at this small scale.The parameters can be easily obtained from calibrating two or three conventional triaxial tests on regular size soil specimens.
The model was used to simulate the stress-strain behavior of two different remolded clay.For each type of clay, using the same set of parameters, we could simulate reasonably well both drained and undrained tests.The model was also used to simulate the stress-strain behavior of a natural clay, taking into account its anisotropic nature caused by the clay deposition during geological formation.Based on the model's performance on these three different clays, the microstructural approach seems to be applicable for capturing the effects of confining stress, OCRs, and inherent anisotropy.

Fig. 5 .
Fig. 5. Comparison of predicted results and experimental results for white kaolinite clay in drained tests: ͑a͒ stress-strain curves; ͑b͒ effective stress paths; ͑c͒ void ratio versus axial strain curves; and ͑d͒ void ratio versus mean stress curves

Fig. 6 .
Fig. 6.Calibration of model parameters for white kaolinite clay

Fig. 7 .
Fig. 7.Comparison of computed and experimental results for white kaolinite clay in undrained tests: ͑a͒ stress-strain curves; ͑b͒ effective stress paths; ͑c͒ pore pressure versus axial strain curves; and ͑d͒ void ratio versus mean stress curves

Fig. 8 .Fig. 9 .
Fig. 8. Influence of parameters a and b on undrained stress path

Fig. 10 .
Fig. 10.Local stresses and strains on planes of various orientations for drained test on black kaolinite clay ͑OCR= 1͒: ͑a͒ local stress paths; ͑b͒ local shear stress-strain curves; ͑c͒ local shear stresses; ͑d͒ local normal stresses; ͑e͒ local shear strains; and ͑f͒ local normal strains

Fig. 11 .
Fig. 11.Local stresses and strains on planes of various orientations for drained test on black kaolinite clay ͑OCR= 8͒: ͑a͒ local stress paths; ͑b͒ local shear stress-strain curves; ͑c͒ local shear stresses; ͑d͒ local normal stresses; ͑e͒ local shear strains; and ͑f͒ local normal strains

Fig. 15 .
Fig. 14.Schematic plot for vertical and horizontal samples

Table 1 .
Model Parameters for Clay