Elasto-plastic Constitutive Model for Overconsolidated Clays

Constitutive model for saturated cohesive soils based on the bounding surface plasticity notion with anisotropic hardening law is presented in the paper. The model predicts inelastic behaviour of overconsolidated cohesive soils. The projection centre is the only point in the stress space which represents elastic soil behaviour. Approximation of the plastic modulus within the preconsolidation domain is made using the radial mapping rule between a projection centre and a reflecting point on the bounding surface. The projection centre changes its position each time when stress path turns rapidly of more than 90°. The configuration of the elliptic bounding surface is governed by preconsolidation effective pressure p′c which depends on change of plastic both volumetric and deviatoric strain. Associated flow rule has been assumed in the formulation. Integration of constitutive relations is done according to forward Euler scheme with error control proposed by Sloan. The effectiveness of the proposed model is illustrated in both monotonic and cyclic loading in the homogeneous triaxial drained and undrained conditions.


Introduction
Critical State Soil Mechanics (CSSM) provides theoretical framework for description of soil behaviour [1]. Within this concept, reaction of soils to loading depends on overconsolidation ratio and varies for normally consolidated soils (OCR = 1) and for overconsolidated soils (OCR > 1). One of the main features of CSSM is the assumption of entirely elastic response of overconsolidated soil and sudden appearance of irreversible deformations when soil reaches normally consolidated state.
Experimental results have proved that irreversible deformations occur in each phase of loading. This fact stimulates researchers to seek theoretical solutions, which account for plasticity also within the range of overconsolidation. Amongst many formulations, two main concepts gained particular popularity in modelling of cohesive soils behaviour, i.e., Bounding Surface Soil Plasticity [2,3] and Subloading Surface models [4,5].
Technical development of laboratory measurement techniques enabled investigation of changes of stiffness during a loading process within the range of small strain. Strong decrease of soil stiffness during straining within the range up to 1% can be illustrated by reduction of the shear and bulk moduli even 15 times with respect to their initial val-1 3 soil mechanics. The concept of the model named NAHOS originates from CSSM and uses the notion of a bounding surface. Its main feature is reduction of the domain of elastic behaviour to a point.
Verification of model capability to predict soil behaviour under monotonic and cyclic loading has been done by simulating loading programmes typical for tests in the conventional triaxial apparatus.

General Remarks
Constitutive relations between Cauchy stress tensor and strain tensor for the sake of small deformation assumption and their symmetry are written in Voigt's vector notation.
Total stress vector can be separated into effective part σ′ and pore pressure u, according to the well-known Terzaghi's principle: where m is the unit vector (equivalent to Kronecker's delta): The superscript T stands for transposition. The stress vector decomposition into deviatoric and hydrostatic part takes the form: where is the stress deviator and (1) 22 , 33 , 12 , 23 , 31 T = x , y , z , xy , yz , zx T , is the mean stress (hydrostatic stress), i.e., first stress invariant.
The second stress invariant in geomechanics is usually defined as The strain vector can be decomposed into deviatoric and volumetric parts: where the volumetric strain (first strain invariant): The second strain invariant, i.e., shear strain, can be defined: where = diag 1, 1, 1, 1 2 , 1 2 , 1 2 . Incremental strain decomposition into elastic and plastic part reads: where the elastic part is given as  Fig. 1 Dependence of clay stiffness on stress and strain level: a undrained shearing of Pentre clay [4]; b torsional shearing of kaolin clay [20] (a) (b) with D isotropic elastic constitutive matrix: K t and G t are elastic tangent bulk and shear moduli, respectively. Plastic strain increment is calculated assuming the associated flow rule: where the plastic multiplier is K P is the plastic modulus associated with σ′ and n F is the unit vector normal to the bounding surface F = 0 (the normalised gradient of the bounding surface): where ‖⋅‖ stands for Euclidean norm.
Since it is assumed that the bounding surface undergoes isotropic changes due to changes of isotropic hardening parameter κ, K P is given in a general form: In case of the associated flow rule, constitutive relations for hardening/softening material are given as follows:

Bounding Surface Concept and the Radial Mapping Rule
Bounding surface concept was introduced independently by Krieg [3] and Dafalias and Herrmann [2]. Contrary to classical plasticity it assumes that unloading/reloading is inelastic which remains in agreement with experimental findings. Unloaded soils are overconsolidated. Stresses in overconsolidated soils are represented in the stress subspace enclosed by a bounding surface.
Another feature of the bounding surface theory is the reduction of the elastic domain to a point, so that the yield surface does not exist. Thus, soil is assumed to be elastoplastic material at the onset of loading.
The behaviour of a soil before failure is governed by variations of the plastic hardening modulus K P . For hardening material, K P > 0, and for softening material, K P < 0. Thus, the basic ingredient of constitutive equations is formulation of the rule governing changes of the plastic modulus inside a bounding surface. This rule may be written in a general form: where K p is the plastic modulus at the point representing current stress state σ′; K P is the plastic modulus at the socalled reflecting point ̄′ (Fig. 2); r, r 0 are distances illustrated in Fig. 2

and defined as
H is a hardening function that scales the value of K P between K P = ∞ at the back stress representing elastic state (elastic centre) and K P =K P .
There are no unique rules for defining the hardening function H. The main assumption concerns equality of unit vectors associated with the reflecting point and the current stress point (Fig. 2). Thus, the plastic multiplier is expressed as follows: Dafalias and Herrmann [2] assumed that the back stress coincides with the origin of the coordinate system. [12] made the back stress moving along p′ axis which is an element of isotropic hardening. In Hashiguchi's concept [5], the back stress moves in the same direction as the current stress. Allowing the back stress to move aside the hydrostatic axis made the concept belonging to models with anisotropic hardening. Dashed line in Fig. 2 represents a loading surface. It is homothetic to the bounding surface and passes through the current stress point. The elastic centre is the centre of homothety with the scale β: By virtue of (25), the reflecting point can readily be found:

Description of the Model NAHOS
The presented model assumes the bounding surface of an ellipsoidal shape in the principal stress space. It is given by equation proposed by Dafalias and Herrmann [2]- Fig. 3: Since the origin of the coordinate system must always lie within or on the surface the value of parameter R must be greater or equal to 2.
The shape of the bounding surface is similar to the wellknown ellipsoid of Modified Cam-clay model for which R = 2. The configuration of the ellipsoid is established by three parameters: R, the slope of the critical state line M (Fig. 2b) and pressure p′ c which is dependent on volumetric plastic variations according to the isotopic exponential hardening law: where p′ c0 and e 0 are initial values of the mean effective pressure for isotropic compression and the void ratio, respectively, e p is plastic part of the void ratio, Δε v p is plastic increment of volumetric strain, λ is the slope of isotropic compression line, and κ is the slope of elastic part during isotropic rebound (Fig. 4).
Elastic tangent bulk and shear moduli are assumed in the usual manner for CSSM: where e is the void ratio and ν is Poisson's ratio.
As pointed out by Borja et al. [13], the bulk modulus in (29) contains assumption that the hysteretic behaviour exhibited by soils during unloading and reloading is small. Assumption of the shear modulus is known to be nonconservative [14].
As mentioned earlier, after each abrupt turning of a stress path, soil stiffness increases significantly. This fact is modelled by assuming the change of a position of the elastic centre when a stress path changes its direction of at least 90°.
In Fig. 5, σ′ Sn−1 and σ′ Sn stand for the positions of the elastic centre before and after the turn of a stress path, respectively; similarly, ̄� n−1 and ̄′ n are the Since in numerical analyses of geotechnical BVP the direction of a stress path is not controlled, abrupt turning can be identified analytically by means of the condition: After detecting, the abrupt turning response of a soil becomes instantaneously elastic and the hardening modulus is interpolated between the new elastic centre and the reflecting point.
To ensure that the elastic centre remains inside or on the bounding surface, i.e., F ′ S ,e p ≤ 0, an additional assumption should be made. Let ′ S0 denotes the elastic centre immediately after changing of its position. Since then, it remains unchanged with respect to the current configuration of a bounding surface (expressed by p' c ). If p ′ c specifies the size of the bounding surface at the abrupt turning of a stress path, then the position of the elastic centre for the current size of a bounding surface must be updated as follows: For the described formulation of the model NAHOS, the normalised gradient of the bounding surface is given as follows: where The formula for the plastic modulus at the reflecting point ̄′ is obtained from the consistency condition dF = 0: The final element that completes the specification of the model is the hardening function H. In the present formulation, the following form of H is proposed: where C and µ are model constants. Function f was proposed by Kaliakin and Dafalias [15] in their model to adjust predictions of shear curves to experimental results for heavily overconsolidated soils (Fig. 6): The quantity n p represents the component in the p direction of the unit gradient of the bounding surface in the space of stress invariants: (35) Although in recent years, various return algorithms became very popular, Potts and Ganendra [23] have proved that substepping algorithm proposed by Sloan [19] can give results that are closer to existing closed solutions. As the particularly attractive approach for the sake of its simplicity and robustness, they recommend algorithm, which is based on the explicit Euler scheme with error control. In this algorithm, an increment of strain is divided into a number of substeps. The substeps vary proportionally over the increment. The lengths of the substeps are determined by setting an error tolerance on the numerical integration.
The explicit modified Euler scheme with error control has been used to integrate numerically equations of (38) Simulation of triaxial undrained cyclic loading at different sizes of applied axial strain increments dε 1 has been performed. The influence of the size of strain increments in the case of model NAHOS has proved to be insignificant (Fig. 7).

General Predictive Capabilities of the Model
The capability of the model to predict soil reaction to loading has been tested for various overconsolidation ratios OCR = 1; 1.2; 1.5; 2; 2.5; 4; 6; 10. In numerical tests, the following set of parameters has been used: Undrained and drained behaviour of the model is presented in Figs. 8 and 9, respectively.

Comparison with Experimental Results
Kaliakin and Dafalias [15] have quoted experimental results obtained on samples of kaolin clay at various values of OCR = 1; 1.3; 2; 6. Samples were tested in a standard triaxial apparatus. Initial state of stress was isotropic.
In each test, the initial effective mean pressure prior to shearing was p′ c0 = 392.2 kPa. Triaxial shearing was carried out in undrained conditions. Comparison of theoretical and experimental results is presented in Fig. 10.

Reduction of Stiffness
As mentioned in the first section, along with the progress of loading between two abrupt turnings of a stress path, there is a significant reduction of soil stiffness. Jardine et al. [7,8] have found that values of moduli may decrease even 15 times. The form of a hardening law in NAHOS aims to reproduce this mechanical feature. For the present elastic constitutive relation, variability of a shear modulus in undrained conditions for three different OCR values is presented in Fig. 11. The superscript "ep" at the modulus G ep has been introduced to distinguish between the elastoplastic (global) shear stiffness and the elastic one declared in (29). The relation for G ep in triaxial conditions may be obtained from the following form of constitutive equations: Then where:

Conclusions
The model described in the paper aims at description of the mechanical behaviour of overconsolidated clays. Experimental results obtained from laboratory and field investigations show that nonlinear irreversible deformations are observed from the beginning of loading. NAHOS makes use of the bounding surface concept to simulate the reaction of clay to monotonic and cyclic loading.
The important feature of the model is the hardening law. The value of the hardening modulus is interpolated between the infinitely large value at the elastic centre and the value at a reflecting point on a bounding surface. The point representing purely elastic behaviour of a soil may occupy different positions within a bounding surface. This enables proper hysteretic behaviour (closed hysteresis loops in Fig. 7) and smooth stress-strain curves at the transition from overconsolidation (stress states inside a bounding surface) to normal consolidation (stress states on a bounding surface).
Clay behaviour predicted in simulations of triaxial tests compared to experimental results exhibits capability of the model to properly reproduce the most important mechanical features of clays in monotonic and cyclic loading programmes. Degradation of soil stiffness, which has been presented in the paper, is strongly dependent on the elastic properties. There are many propositions for describing elastic part of deformations which can be found, e.g., in [11] and they should be implemented and tested within the general framework of NAHOS.
Presented version of the model still suffers from several drawbacks. It assumes circular shape of a bounding surface in deviatoric plane, which is equivalent to the assumption that the strength is independent of the direction of loading. Associated flow rule assumed in the model does not allow to predict material instabilities which may occur before reaching the failure condition [24,25]. These traits must be taken into account in the forthcoming versions of the model.