Transp. Infrastruct. Geotech. https://doi.org/10.1007/s40515-018-0060-3 TECHNICAL PAPER
Modeling the Time- and Temperature-Dependent Response of Cohesive Soils in a Generalized Bounding Surface Framework Victor N. Kaliakin1 · Andr´es Nieto-Leal2 · Meysam Mashayekhi1
Accepted: 15 May 2018 © Springer Science+Business Media, LLC, part of Springer Nature 2018
Abstract The Generalized Bounding Surface Model (GBSM) synthesizes many previous bounding surface models for cohesive soils and improves upon their predictive capabilities. Following an overview of the GBSM and its extension to account for non-isothermal conditions, the predictive capabilities of the model are assessed. Keywords Cohesive soils · Bounding surface · Time-dependent · Temperature · Constitutive model
Introduction The results of experimental studies indicate that as a saturated cohesive soil is loaded, it continuously develops both elastic and inelastic strains without a distinct yield state that delineates elastic from inelastic response [69, 79]. The absence of a clearly defined yield state has been attributed to the highly redundant clay microstructure that contains a wide variety of stress magnitudes at points of interparticle contact and a wide range of bond strengths [8, 11].
Victor N. Kaliakin
[email protected] 1
Department of Civil and Environmental Engineering, University of Delaware, Newark, DE 19716, USA
2
Department of Civil Engineering, Universidad Militar Nueva Granada, Cajic´a, 250247, Colombia
Transp. Infrastruct. Geotech.
The absence of a sharp delineation between elastic and inelastic states has led to the development of constitutive formulations for cohesive soils based on the concept of a subloading surface [35], an extended subloading surface [34], and a bounding surface in stress space [20]. The essence of such “unconventional” plasticity models is the hypothesis that inelastic deformations can occur for stress states either within or on a properly defined outer loading or bounding surface [18]. Thus, unlike classical elastoplasticity models, inelastic states are not restricted only to those states lying on the aforementioned outer surface. Initially, constitutive models based on the bounding surface concept were used to successfully simulate the rate-independent response of cohesive soils [20, 21, 59–61]. These subsequently evolved into more formal two- [27, 28] and three-invariant [29, 31] rate-independent formulations that employed a critical state failure criterion [70]. The latter formulation was next synthesized [30] and subsequently simplified [44]. A bounding surface formulation able to account for inherent and stress-induced anisotropy was first discussed by Dafalias [24]. This led to the subsequent development of anisotropic, rate-independent bounding surface models [2, 54]. The theoretical basis for a time-dependent version of the bounding surface model for isotropic cohesive soils was presented by Dafalias [23, 24]. This model was subsequently refined and implemented numerically [41] and then formally presented [45] and verified [46]. The logical union of the anisotropic time-independent bounding surface model with the isotropic time-dependent one led to somewhat complex formulations that involved rather large numbers of model parameters [1, 76]. Realizing the overly complex nature of these earlier formulations, Ling et al. [57] developed a more rational model that reduced the total number of model parameters without detrimentally affecting its predictive capabilities [56, 80]. Jiang et al. [39] extended this model to include a non-associative flow rule, thus facilitating the simulation of cohesive soils that exhibit softening response. Kaliakin [43] first outlined the formulation and numerical implementation of a time-dependent bounding surface model extended for thermo-hydro-mechanical (THM) analyses. However, the lack of suitable experimental results precluded realistic application of the model. Subsequently, this limitation was alleviated, and several bounding surface-based THM models were developed [17, 37, 51–53, 58, 81, 82]. Although all of the aforementioned models for cohesive soils were based on the concept of a bounding surface in stress space, they improved upon earlier versions of such models, enhanced the predictive capabilities of earlier models by expanding their features, or accomplished both of these tasks. Missing from the earlier development of bounding surface models for cohesive soils was any attempt to synthesize the many previous forms of these models. Such a synthesis has now been realized through the Generalized Bounding Surface Model (GBSM) for cohesive soils [47, 48]. The GBSM also improves upon many aspects of previous models for cohesive soils and expands the scope of the model to THM analyses. This paper presents a detailed overview of the GBSM and discusses its extension to account for non-isothermal conditions. The predictive capabilities of the model are also assessed.
Transp. Infrastruct. Geotech.
Model Overview In its most general form, the GBSM for cohesive soils is a fully three-dimensional, temperature and time-dependent model that accounts for both inherent and stressinduced anisotropy. To better simulate the behavior of cohesive soils exhibiting softening, the model employs a non-associative flow rule. To more accurately simulate the behavior of cohesive soils subjected to cyclic loading, the GBSM uses a robust general methodology for locating the projection center at any point in stress invariant space. The microfabric-inspired rotational hardening rule associated with the model was developed after a thorough review of past modeling practice. Finally, the shape hardening function used in the GBSM was developed simplifying earlier versions without compromising the model’s predictive capabilities. In the subsequent development tensorial quantities are presented in indicial form with the indices obeying the summation convention over repeated indices. A single subscript for the internal variables is not a tensorial index but merely identifies the plurality of these quantities. Compressive stresses and strains are taken as positive, and infinitesimal strains are assumed.
Overview of the Bounding Surface Concept The bounding surface in stress space is smooth and convex [26, 63]. It always encloses the origin and is origin-convex, i.e., any radius emanating from the origin intersects the surface at only one point (Fig. 1). The bounding surface is defined analytically by F σ¯ ij , qn = 0 (1)
Fig. 1 Schematic illustration of the bounding surface, a loading surface, and radial mapping rule in multiaxial stress space
Transp. Infrastruct. Geotech.
where σij is the effective stress tensor and qn represents the inelastic internal variables. The bar over σij indicates an “image” point on the bounding surface. The essence of the bounding surface concept is the hypothesis that inelastic deformations can occur for stress states either within or on the surface. Thus, unlike classical yield surface elastoplasticity, the inelastic states are not restricted only to those lying on an outer surface. The bounding surface is instrumental in defining the direction of inelastic loadingunloading (Lij ) and the plastic modulus (Kp ) that enter into the formulation. The expression for Lij at the actual stress point (σij ) is defined as the gradient of F at the “image” point, viz., ∂F (2) Lij = ∂ σ¯ ij For any stress increment σ˙ ij causing inelastic loading, a corresponding “image” stress increment σ˙¯ ij occurs as a result of the hardening of the bounding surface by means of the internal variables qn . The superposed dot indicates a material time derivative or rate. Three additional expressions are required to complete the bounding surface formulation. First, a scalar loading index (L) must be defined in terms of Lij , σ˙ ij and σ˙¯ ij , the plastic modulus Kp (associated with σij ) and a “bounding” plastic modulus K¯ p (associated with σ¯ ij ). Second, an expression for K¯ p must be obtained from the consistency condition F˙ = 0. Finally, a state dependent relation between Kp and K¯ p , that is established as a function of the Euclidean distance (δ) between σij and σ¯ ij , must be determined, viz., δ Kp = K¯ p + Hˆ σij , qn (3) r −δ where r is a reference stress or distance from the projection center aij to the “image” stress on the bounding surface (Fig. 1) and Hˆ is a positive scalar “shape hardening” function of the state. All bounding surface formulations require a “mapping rule” that relates σij to σ¯ ij . The radial mapping rule [22], which combines simplicity and ease of numerical implementation, is used herein. It has been shown to accurately predict the rateindependent [2, 27–31, 39, 57] and time-dependent [1, 45] behavior of isotropic and anisotropic cohesive soils, as well as sands and clays in conjunction with “unified” forms of the model [13, 14, 16]. In the radial mapping form of the model, shown schematically in Fig. 1, the projection center tensor aij is assumed to always lie within the bounding surface and to never cross it. It evolves according to a proper rate equation, is one of the qn , but does not enter into Eq. 1. Using aij as the projection center, the “image” stress is obtained by the radial projection of σij onto the surface according to (4) σ¯ ij = b σij − aij + aij where the dimensionless parameter b (≥ 1) is determined in terms of the material state by substituting σ¯ ij from Eq. 4 into an explicit form of Eq. 1 and solving the resulting expression for b.
Transp. Infrastruct. Geotech.
From Eq. 4, b = σ¯ ij − aij / σij − aij . From Fig. 1, r = σ¯ ij − aij and r − δ = σij − aij . Comparing these equations, it follows that b = r/(r − δ), which leads to the useful relations r/δ = b/(b − 1) and δ/(r − δ) = b − 1. A consequence of assuming the radial mapping rule is that a loading surface homologous to the bounding surface with respect to the projection center aij and passing through the actual stress point (σij ) is implicitly defined (Fig. 1). This surface determines all the paths of neutral loading emanating from σij . The radial mapping rule also implicitly defines a quasi-elastic domain referred to as the “elastic nucleus” [30]. The presence of this domain in the formulation necessitates the replacement of the denominator (r − δ) in Eq. 3 by the quantity < r − sp δ >, where r and δ are as previously defined, and sp (≥ 1) is a dimensionless model parameter that defines the extent of the elastic nucleus. The symbols < > denote Macaulay brackets, which imply that, for some quantity d, < d >= d if d > 0 and < d >= 0 if d ≤ 0.
General Coupled Elastoplastic-Viscoplastic Formulation This section gives an overview of a general bounding surface elastoplasticviscoplastic formulation for cohesive soils. Additional details pertaining to this formulation are given elsewhere [48]. The following additive decomposition is assumed for the infinitesimal strain rate tensor ε˙ ij : p e i e v + ε˙ ij = ε˙ ij + ε˙ ij + ε˙ ij (5) ε˙ ij = ε˙ ij where the superscripts e, i, v and p denote its elastic, inelastic, viscoplastic (delayed) and plastic (instantaneous) components, respectively. The internal variables include, in general, proper measures of inelastic deformation that quantify the hardening of the bounding surface. Knowing the material state and the rate σ˙ ij , the constitutive relations are thus given by a proper analytical expression for the elastic strain rate and the q˙n . The following additive decomposition of p q˙n into a delayed or viscoplastic part (q˙nv ) and an instantaneous or plastic part (q˙n ) is assumed [45]: p q˙n = q˙nv + q˙n (6) Both parts are, in general, incrementally irreversible. The q˙nv , which accounts for the delayed deformation of the cohesive soil, depends only upon the state; it can, p under certain conditions, be zero. The q˙n , on the other hand, depends upon the state and the rates σ˙ ij and q˙nv are non-zero only if the scalar loading (L) index is positive. This is done in order to emphasize the general coupling that exists between plastic and viscoplastic processes, whereby in addition to the σ˙ ij , the q˙nv may affect the plastic loading process [23]. Finally, recalling that the projection center tensor aij is also one of the qn but distinguishing it for emphasis, it follows that aij evolves according to the following rate equation: p a˙ ij = a˙ ijv + a˙ ij (7)
Transp. Infrastruct. Geotech.
The Elastic Response The response associated with the elastic strain rate is assumed to be independent of the rate of loading. The elastic constitutive relations, in direct and inverse form, are thus e e ε˙ ij = Bij kl σ˙ kl or σ˙ ij = Dij kl ε˙ kl (8) where Bij kl and Dij kl are the fourth-order tensors of elastic compliance and moduli, respectively, being in general functions of the state. Historically, bounding surface models for cohesive soils have exclusively assumed elastic isotropy. The two material parameters associated with an isotropic elastic idealization are typically the bulk modulus (K) and either the shear modulus (G) or Poisson’s ratio (ν). The elastic parameters are assumed to be independent of the qn . In anisotropic forms of the model it is also possible to assume anisotropic elasticity. Although experimental findings [49] suggest that natural soils are elastically orthotropic, the difficulty associated with determining values for the nine elastic constants precludes the adoption of such an idealization. Instead, transversely isotropic elasticity is commonly assumed, thus reducing the number of elastic parameters to five. The Viscoplastic Response v and q v are based on a generalization of The viscoplastic rate equations for εij n Perzyna’s elastic/viscoplastic theory [66, 67]. The measure of “distance” in stress space from a “static” yield surface associated with this theory is taken to be the normalized overstress Δσˆ , where δˆ (9) Δσˆ = r 1 − s1v
where r is as defined in Fig. 1, sv (> 1) is a dimensionless model parameter, and δˆ is the Euclidean distance between σij and a second “image” stress σˆ ij on the boundary of a second elastic nucleus that is associated with the viscoplastic response and is homologous with the bounding surface (Fig. 2). This elastic nucleus assumes the role of the aforementioned “static” yield surface. Since the material is assumed to be v becomes zero when Δσ inviscid in the elastic region, ε˙ ij ˆ ≤ 0, i.e., when σij is either on or within the second elastic nucleus. The overstress Δσˆ enters the rate equations through the continuous scalar “overstress” function Φ, such that Φ > 0 when Δσˆ > 0 and Φ ≤ 0 when Δσˆ ≤ 0. The rate equations for the viscoplastic response are thus v ε˙ ij = Φ Rij
;
q˙nv = Φ rn
;
a˙ ijv = Φ rij
(10)
where Rij , rn and rij represent proper tensorial functions of the state. If a nonassociative flow rule is assumed, the tensor Rij is the gradient of a plastic potential Q, i.e., Rij = ∂Q/∂ σ¯ ij . If, on the other hand, an associative flow rule is assumed, then Rij = ∂F /∂ σ¯ ij = Lij (recall Eq. 2). A commonly used form of the overstress function is Φ = (Δσˆ )n /V , where n and V are positive model parameters.
Transp. Infrastruct. Geotech.
Fig. 2 Schematic illustration of the bounding surface, a loading surface, and radial mapping rule in multiaxial stress space
The Plastic Response As noted in the Overview of the Bounding Surface Concept section, one of the key quantities in formulating elastoplastic response is a properly defined scalar loading index L. Consistent with the framework of the radial mapping rule, L is defined as follows [23]: 1 1 ∂F 1 ∂F L= rn − Φ 1 − rij Lij σ˙ ij + Φ (11) Kp b ∂qn b ∂ σ¯ ij where Lij is given by Eq. 2 and b is as previously defined in conjunction with Eq. 3. When the stress point lies on the bounding surface, b = 1, Kp = K¯ p , and σ˙ ij = σ˙¯ ij ; Eq. 11 then simplifies accordingly. The presence of rn and rij in Eq. 11 shows that L couples the plastic-viscoplastic hardening response for states on and within the bounding surface. Kaliakin and Dafalias [45] give additional details pertaining to the three terms appearing in this equation. The rate equations for plastic response are obtained by assuming linear depenp dence of q˙n on the rates σ˙ ij and q˙nv , giving p
ε˙ ij = L Rij
;
p
q˙n = L rn
;
p
a˙ ij = L rij
(12)
In light of their dependence on the q˙nv through Eq. 11, the assumption of linearity p does not make the q˙n rate-independent. This means that as a result of an instantap p neous change of σ˙ ij to cσ˙ ij (where c is a positive scalar), q˙n does not become cq˙n v because the q˙n remains instantaneously unchanged. An instantaneous increase of the p q˙n occurs, however, because of its dependence upon σ˙ ij . The word “instantaneous”
Transp. Infrastruct. Geotech. p
attached to the q˙n must be viewed in a general manner, for due to the dependence p of the loading index on q˙nv , it is possible to have q˙n = 0 even though σ˙ ij = 0. The overall response is obviously time-dependent. In light of Eqs. 1 and 6, the consistency condition is written as follows: ∂F ∂F p (13) F˙ = Lij σ¯˙ ij + v q˙nv + p q˙n = 0 ∂qn ∂qn Evaluating (11) at the bounding surface (i.e., for b = 1), solving for Lij σ˙ ij , substituting this result into Eq. 13 and using the second of Eqs. 10 and 12 gives the following general expression for the “bounding” plastic modulus entering (3): ∂F K¯ p = − rn (14) ∂qn where the condition L > 0 (i.e., plastic loading) has been assumed. The Total Response Substituting (5) into the second of Eq. 8 gives a general expression for the inverse form of the constitutive relations, viz., p v σ˙ ij = Dij kl ε˙ kl − ε˙ kl − ε˙ kl = Dij kl (˙εkl − Φ Rkl − L Rkl ) (15) where the first of Eqs. 10 and 12 have been used. Using Eq. 15 to form the product Lij σ˙ ij , and substituting this expression into the first of Eqs. 11 gives, after some simple manipulation, the following expression for the loading index: 1 ∂F 1 Lij Dij kl ε˙ kl − Φ Lij Dij kl Rkl + Φ rn − Φ 1 − Lij rij b ∂qn b L= (16) Kp + Lab Dabcd Rcd where L > 0 has been assumed. Substituting (16) into (15) gives the final form of the constitutive relations, viz., σ˙ ij = Cij kl ε˙ kl − Vij ep
(17)
where the elastoplastic tangent stiffness tensor is given by H(L) (18) Dij st Rst Lmn Dmnkl D and the viscoplastic righthand side vector is
Vij = Φ Dij kl Rkl H(L) 1 ∂F 1 − rn + 1− Rmn rmn Rkl (19) Dij kl Lmn Dmnqr Rqr − D b ∂qn b ep
Cij kl = Dij kl −
In Eqs. 18 and 19
1 if L > 0 0 otherwise is the Heaviside step function and D = Kp + Lab Dabcd Rcd . H(L) =
(20)
Transp. Infrastruct. Geotech.
The extension to non-isothermal conditions requires the inclusion of two incremental strains that must be subtracted on the righthand side of Eq. 15. The first of these quantities, which represents the reversible part of the thermal strain, is βs (21) δkl T˙ 3 where β s is the coefficient of volumetric thermal expansion for the solid phase, δkl is the Kronecker delta, and T is the temperature. The second additional incremental strain is the irreversible (thermoplastic) part of the thermal strain, which is given by
Tp ε˙ kl = T˙ fT σij , qn , T Rkl (22) Te ε˙ kl =
where fT is a function of the effective stress state and certain internal variables qn . For a specific elastic idealization Dij kl , the complete definition of the GBSM constitutive relations given by Eqs. 18 and 19 thus requires explicit forms for the gradients Rij = ∂Q/∂ σ¯ ij and Lij = ∂F /∂ σ¯ ij , explicit expressions for K¯ p and Hˆ σij , qn entering (3), and suitable hardening rules for the plastic potential and bounding surfaces. The aforementioned gradients require explicit definition of the plastic potential surface (Q) and the bounding surface (F ). The definition of these surfaces requires the selection of suitable invariants, which are next discussed.
Definition of Invariants The GBSM accounts for inherent and stress-induced anisotropy. The former is attributed to preferred particle arrangement and interparticle forces developed during sedimentation; it is a physical characteristic of the soil that is independent of stresses that are subsequently applied to the soil [9]. Stress-induced anisotropy, on the other hand, is a physical characteristic that almost exclusively results from the deformation associated with changes in applied stress. Guided by the results of experimental studies, inherent anisotropy is accounted for by a suitable initial rotation of the bounding surface (BS). The conjecture underlying such an initial configuration is that the rotation corresponds to the predominant orientation of the particles or particle clusters in the soil. If a non-associative flow rule is used, the plastic potential surface (PPS) must likewise be rotated. The mathematical representation of stress-induced anisotropy typically requires the introduction of rotational hardening into a constitutive model [25]. Since the most general form of the GBSM employs a non-associative flow rule, the rotation of both the PPS and of the BS must be accounted for. The material state is thus not determined solely by the state of σij , but also by a suitable measure of anisotropy [3]. Typically, a symmetric anisotropic tensor is used to account for the anisotropy. In past anisotropic, non-associative bounding surface models [38, 40], the same symmetric second-order anisotropic tensor was used to represent the rotation of both of these surfaces. However, guided by recent findings [32], different anisotropic tensors are used to quantify the rotations of the PPS and BS [65] in the GBSM.
Transp. Infrastruct. Geotech.
In the anisotropic form of the GBSM employing a non-associative flow rule, the rotation of the PPS is described by the symmetric second-order anisotropic tensor αij [48]. The following invariants are thus used in the analytical description of the PPS: √ α α α 1 α α 3 sij sj k ski 1 −1 α α I = σij δij = σkk ; J = sij sij ; θ = sin (23) 2 3 2 (J α )3 where δij is again the Kronecker delta, I is the first invariant of σij , J α is the square root of the second reduced deviatoric stress invariant, and θ α (−π/6 ≤ θ α ≤ π/6) is the reduced Lode angle [15]. In Eq. 23 sijα is the reduced deviatoric stress tensor, which is related to the deviatoric stress tensor sij as follows: 1 sijα = sij − αij σkk 3
(24)
where 1 (25) sij = σij − δij σkk 3 In anisotropic forms of the GBSM, the rotation of the BS is described by the symmetric second-order anisotropic tensor βij . The analytical description of the BS is written in terms of I , as defined in Eq. 23, and the square root of the second reduced deviatoric stress invariant (J β ) and the reduced Lode angle θ β , viz., √ β β β 1 β β 3 sij sj k ski 1 −1 β β sij sij ; θ = sin (26) J = 3 2 3 2 Jβ where the reduced deviatoric stress tensor is now given by 1 1 β sij = sij − βij σkk = σij − (27) βij + δij σkk 3 3 where sij is the deviatoric stress tensor given by Eq. 25. In isotropic forms of the GBSM, both the PPS and BS are analytically defined by I , as given in Eq. 23, and by the following isotropic invariants: √ 3 sij sj k ski 1 1 −1 J = sij sij ; θ = sin (28) 2 3 2 (J )3 where sij is as defined by Eq. 25, and −π/6 ≤ θ ≤ π/6. In axisymmetric triaxial (p √ - q) stress space, where p = I /3 is the mean normal effective stress and q = ± 3 J is the deviator stress, the rotation of the PPS is quantified by the dimensionless scalar variable α, where 3 (29) αij αij α= 2 In a similar manner, the rotation of the BS is quantified by the dimensionless scalar variable β, where 3 (30) βij βij β= 2
Transp. Infrastruct. Geotech.
Lode Angle Dependence of Surfaces Experimental evidence indicates that, within the octahedral plane, the PPS, BS and the failure criterion should vary with Lode angle (in direct or reduced form). The specific functional form for this variation is by no means unique. In keeping with previous forms of the bounding surface model for cohesive soils, the variation with Lode angle of a generic model parameter P associated with one of the aforementioned surfaces is given by ˜ = g(θ, ˜ k)Pc P (θ) (31) where θ˜ is the direct or reduced Lode angle and k = Pe /Pc , with Pe = P (−π/6) and Pc = P (π/6) being the values of P associated with axisymmetric triaxial extension and compression, respectively. The dimensionless function g(θ˜ , k) must thus take on the values g(−π/6, k) = k and g(π/6, k) = 1. ˜ k , attributed to Gudehus [33] and The following general form of g θ, Argyris [83], has been used in conjunction with previous bounding surface models for clays [29–31, 44]: 2k ˜ k = (32) g θ, 1 + k − (1 − k) sin 3θ˜ The above form of g(θ˜ , k) gives convex failure surfaces only for k ≥ 0.778, which correspond to effective friction angles less than 22◦ [72]. A better choice for this function thus appears to be
˜ k = g θ,
2k 4 1 + k 4 − (1 − k 4 ) sin 3θ˜
14 (33)
which gives convex failure surfaces for k ≥ 0.60 (corresponding to effective friction angles less than or equal to 48.6◦ ) [72]. Due to this additional robustness, Eq. 33 now used exclusively in simulations obtained using the GBSM. Since derivatives of the PPS and BSwith respect to θ˜ are required, in light of Eq. 31, it follows that ∂P /∂ θ˜ = ∂g/∂ θ˜ Pc , where ∂g 3g(1 − k 4 ) cos 3θ˜ = ∂ θ˜ 4 1 + k 4 − 1 − k 4 sin 3θ˜
(34)
If Pc = Pe , then k = 1 and P does not vary with Lode angle. For this case, (34) gives ∂g/∂ θ˜ = 0, indicating that there is indeed no variation of P with Lode angle.
Definition of the Plastic Potential Surface The most general form of the PPS associated with the GBSM is given by ω2 2 R−2 α ¯α ¯ ¯ ¯ ¯ Q I , J , θ , Iα , α = I − Iα I + Iα + (R − 1)2 J¯α = 0 (35) 27 R
Transp. Infrastruct. Geotech.
where I¯ and J¯α are values of I and J α associated with the “image”√point on the PPS. The quantity Iα is the value of I at J¯α = α˜ I¯, where α˜ = α/(3 3). The value of Iα is adjusted so that Eq. 35 is satisfied for the current “image” stress point I¯, J¯β on the PPS. The dimensionless model parameter R (≥ 2.0) controls the shape of the PPS [44]. Finally, ω is the generalized expression for the critical state line given by
M −α 2 2 2 2 ω = (36) 2α (R − 1) + M − α + 4α (R − 1) M + (M − α) 2 where α is as defined in Eq. 29. The quantity M is the slope of the failure surface which, in a meridional section (i.e., for a specific value of the θ α ) in p - q space, is assumed to be a straight line that coincides with √ the critical state line [70]. The counterpart of M in I - J α space is M˜ = M/(3 3). Both M and M˜ vary with θ α according to (37) M θ α , kM = g θ α , kM Mc ; M˜ θ α , kM = g θ α , kM M˜ c where kM = Me /Mc = M˜ e /M˜ c , with Me = M(−π/6) and Mc = M(π/6) being the values of M associated with triaxial extension and compression, respectively. The dimensionless function g is given by Eq. 33, with θ α and kM replacing θ and k, respectively. Fig. 3 shows the variation of the PPS and failure surface with Lode angle within the octahedral plane for anisotropic forms of the GBSM. Solving (35) for Iα leads to the following quadratic equation: ¯ 2 2 I 27 − 1) R−2 (R 2 2 Iα − J¯α + I¯ =0 (38) (Iα )2 + R R ω2
Fig. 3 Plastic potential and failure surface for anisotropic form of the GBSM employing a non-associative flow rule
Transp. Infrastruct. Geotech.
For R > 2.0, the only meaningful root of Eq. 38 is 2 27 (R − 1)2 ¯α 2 ¯ 2 ¯ ¯ −I + + I I + R (R − 2) J ω2 Iα = (R − 2)
(39)
If R = 2.0, Eq. 38 reduces to a linear equation in Iα , the solution of which is 2 27 J¯α + I¯ (40) Iα = ω I¯ For stress states on the I -axis, J¯α = 0; Eqs. 39 and 40 then reduce to the desired result of Iα = I¯. The PPS given by Eq. 35 must satisfy two fundamental requirements. Fora given value of Lode angle (θ α ), it must pass through the “image” stress point I¯, J¯β . Secondly, it must give ∂Q/∂ I¯ = 0 at its intersection with the critical state (failure) line, as defined by ω; the latter requirement is in keeping with the concept of a critical state [70]. Experience [32, 38] shows that the simplest way in which to satisfy these requirements is to set R = 2.0 in Eq. 36. This results in ω2 = M 2 −α 2 , which holds for both the anisotropic (α = 0) and isotropic (α = 0) forms of the GBSM. In the former case, the formulation resembles the non-associative anisotropic bounding surface model developed by Jiang and co-workers [38, 40]. If R > 2.0, then it is only for the case of the anisotropic, associative flow version of the GBSM that ω, as given by Eq. 36 and with R > 2.0, is actually needed [15]. Finally, for the isotropic form of the GBSM employing a non-associative flow rule, α = 0 and Eq. 36 reduces to ω2 = M 2 . In light of the invariants defined in the Definition of Invariants section, the gradient to the PPS, required for Eqs. 15 and 16, is thus ∂Q ∂Q ∂ I¯ ∂Q ∂ J¯α ∂Q ∂ θ¯ α ≡ R = + + (41) ij ∂ σ¯ ij ∂ θ¯ α ∂ σ¯ ij ∂ I¯ ∂ σ¯ ij ∂ J¯α ∂ σ¯ ij where ∂Q/∂ θ¯ α = (∂Q/∂ω) (∂ω/∂M) ∂M/∂ θ¯ α . Kaliakin and Nieto-Leal [48] give explicit expressions for all of the partial derivatives appearing in Eq. 41.
Definition of the Bounding Surface The results of an extensive study into suitable analytical expressions for bounding surfaces [63] showed that an elliptical form avoids difficulties associated with higher-order expressions. In light of this finding, the following elliptical BS is used in conjunction with the GBSM: 2 N − β2 ¯ R−2 F I¯, J¯β , θ¯ β , Io , β = I − Io I¯ + Io 27 R 2 (42) + (R − 1)2 J¯β = 0
Transp. Infrastruct. Geotech.
where I¯ and J¯β are values of I and J β associated with the “image” point on the BS (recall Fig. 1), and R and β are as previously defined. The quantity Io represents the point of intersection of the bounding surface with the positive part of the hydrostatic axis in invariant stress space (Fig. 4). The parameter N defines √ the shape of the BS; ˜ where N˜ = N/(3 3). Both N and N˜ its counterpart in stress invariant space is N, vary with θ β according to (43) N θ α , kN = g θ α , kN Nc ; N˜ θ α , kN = g θ α , kN N˜ c where kN = Ne /Nc = N˜ e /N˜ c , with Ne = N(−π/6) and Nc = N(π/6) being the values of N associated with triaxial extension and compression, respectively. The dimensionless function g is given by Eq. 33, with θ β and kN replacing θ and k, respectively. Unlike Eqs. 35, 42 does not contain its own counterpart for ω; instead, ω is replaced by the quantity (N 2 − β 2 ). This simplification is explained by the fact that there is no need to ensure normality of Lij = ∂F /∂ σ¯ ij on the bounding surface when β = N. For the isotropic form of the GBSM employing a non-associative flow rule, β = 0. Eq. 42 then reduces to N2 2 R−2 ¯ ¯ ¯ ¯ ¯ F I, J, θ = (44) I − Io I + Io + (R − 1)2 J¯ = 0 27 R If an associative flow rule is used in conjunction with the anisotropic form of the GBSM, the bounding surface becomes identical to the plastic potential, i.e., F ≡ Q. Thus, in light of Eq. 35, and replacing βij by αij , β by α, and Iα by Io , the equation of the bounding surface becomes 2 ω2 ¯ R−2 ¯ (45) I − Io I + Io + (R − 1)2 J¯α = 0 F = 27 R where ω is given by Eq. 36 and R ≥ 2.0. The GBSM is then quite similar to the anisotropic bounding surface model developed by Ling et al. [57]. J N Elastic Nucleus
J1
Bounding Surface ( I, J ) ( I, J )
δ
r/sp
I =CI c
0
r
I
0
I
I1=I0/R
Fig. 4 Schematic illustration of elliptical bounding surface and radial mapping rule in meridional section in stress invariant space for isotropic forms of the GBSM
Transp. Infrastruct. Geotech.
The necessity of using ω in Eq. 45 is explained by the fact that if ω was replaced by M, there would be no way to guarantee that ∂F /∂I = 0 at the intersection of the bounding surface and the critical state line for values of R > 2.0. As a result, a fundamental requirement associated with Critical State soil mechanics would not be satisfied. If the GBSM is further specialized for isotropic cohesive soils, then α = 0, ω = M and J α = J β = J . Equation 42 then becomes M2 2 R−2 ¯ ¯ ¯ ¯ ¯ F I , J , θ, Io = (46) I − Io I + Io + (R − 1)2 J¯ = 0 27 R The GBSM is then quite similar to the isotropic bounding surface model developed by [44, 45]. For this form of the GBSM, differentiation of Eq. 46 with ω = M gives ∂F /∂I = 0 at the intersection of the bounding surface and the critical state line. In light of the invariants defined in the Definition of Invariants section, the gradient to the BS, required for Eqs. 15 and 16, is thus ∂F ∂F ∂ I¯ ∂F ∂ J¯β ∂F ∂ θ¯ β ≡ Lij = + + ∂ σ¯ ij ∂ θ¯ β ∂ σ¯ ij ∂ I¯ ∂ σ¯ ij ∂ J¯β ∂ σ¯ ij
(47)
where ∂F /∂ θ¯ β = (∂F /∂N )(∂N /∂ θ¯ β ). Kaliakin and Nieto-Leal [48] give explicit expressions for all of the partial derivatives appearing in Eq. 47.
Explicit Expression for Radial Mapping Rule Specializing the radial mapping rule given by Eq. 4 to stress invariant space gives β β (48) I¯ = b (I − IP C ) + IP C ; J¯β = b J β − JP C + JP C β
where (IP C , JP C ) are the coordinates of the projection center. For monotonic loading, the projection center must be an isotropic tensor with a principal value on the I -axis in invariant stress space, i.e., aij = δij IP C /3 [30]. β PC β = 0, then JP C = 0 (recall Eq. 26). The second of Eq. 48 thus reduces Since sij β β β β β β β β to J¯β = bJ β . From Eq. 26, it follows that s¯ = bs , s¯ s¯ s¯ = b3 s s s and β β β (¯sij s¯j k s¯ki )/(J¯β )3
ij
β β β (sij sj k ski )/(J β )3 .
ij
ij j k ki gives θ¯ β = θ β .
ij j k ki
= Thus, Eq. 26 To relate it to I¯, the aforementioned expression for aij is substituted into Eq. 4. Recalling the definition of I given by Eq. 23 then gives 1 1 I¯ = σ¯ ij δij = b σij δij − δij δij IP C + δij δij IP C 3 3 = b (I − IP C ) + IP C (49) where, from the definition of the Kronecker delta tensor, δij δij = 3. Commonly, IP C = CIo , where C is the dimensionless projection center parameter (0 ≤ C < 1). It is timely to note that realistic simulation of the response to cyclic loading necessitates a more general relocation of the projection center [62]. This is, however, beyond the scope of the present paper.
Transp. Infrastruct. Geotech.
The dimensionless quantity b (≥ 1) appearing in Eq. 49 is determined in terms of the material state by substituting Eqs. 48 into Eq. 42. The solution for the resulting quadratic equation in b is −B ∗ + (B ∗ )2 − A∗ C ∗ b= (50) A∗ where
2 β A∗ = D ∗ (I − IP C )2 + (R − 1)2 J β − JP C Io β β B ∗ = D ∗ (I − IP C ) IP C − + (R − 1)2 JP C J β − JP C R
2 2I R−2 o β C ∗ = D ∗ IP C IP C − − (Io )2 + (R − 1)2 JP C R R 2 2 N −β D∗ = 27
(51) (52) (53) (54)
Hardening of the Plastic Potential and Bounding Surfaces As in the case of previous isotropic and anisotropic bounding surface models, the BS undergoes isotropic hardening along the hydrostatic (I ) axis. The proper simulation of stress-induced anisotropy, necessary for anisotropic forms of the GBSM, requires the PPS and BS to also undergo rotational hardening (RH). Isotropic Hardening The BS is assumed to undergo isotropic hardening along the hydrostatic (I ) axis. The hardening is controlled by a single scalar internal variable that measures the increment in inelastic volumetric strain, viz., p
i v ε˙ kk = ε˙ kk + ε˙ kk
(55)
Denoting the total void ratio by e and by ee , ei , ev and ep its elastic, inelastic, viscoplastic and plastic parts, respectively, leads to the following relations between i and the ei are related these quantities: e = ee + ei = ee + ev + ep . The quantities εkk as follows: v p e˙i = − (1 + ein ) ε˙ kk + ε˙ kk = − (1 + ein ) (φ Rkk + L Rkk ) (56) where ein represents the initial total void ratio corresponding to the reference configuration with respect to which engineering strains are measured. For natural strains ein represents the current total void ratio. From Eq. 56, it follows that the ei may then be chosen as the single scalar internal variable that controls the hardening of the bounding surface. It is convenient to relate the evolution of the BS to the inelastic void ratio through the value of Io , which is the point of intersection of the BS with the positive part of
Transp. Infrastruct. Geotech.
the I -axis in invariant stress space. It follows that Io thus measures the amount of preconsolidation of the soil. For an isotropic soil consolidated isotropically (i.e., along the hydrostatic axis), only volumetric strains are generated. This fact, in conjunction with Eq. 56, implies that Io must depend only on the ei . An expression for dIo /dei , necessary for the analytical description of the isotropic hardening behavior, is given by [30, 48] Io − IL + IL dIo =− (57) i de λ−κ where the critical state parameters λ and κ are equal to the negative of the slopes of the virgin consolidation and swell/recompression lines, respectively, in a plot of e versus the natural logarithm of I [70]. In addition, λ is assumed to be the negative of the slope of the straight line trace of the failure surface in e versus ln I space. The nonzero limiting “transitional” stress IL is assumed such that for I < IL the relation between I and the elastic part of the void ratio (ee ) changes continuously from logarithmic to linear [30]. In this manner, the singularity of the elastic stiffness near I = 0 (resulting from excessive material softening) is thus removed. IL is not a model parameter; its value is typically taken equal to one-third of the atmospheric pressure (Pa ). Underlying (57) is the important assumption that this relation remains valid for ei changing plastically, viscoplastically, or as a combination of the two [10]. This assumption constitutes the basic link of hardening interaction between plastic and viscoplastic mechanisms. Rotational Hardening To properly account for the evolution of stress-induced anisotropy during shearing, a suitable RH rule must be provided so as to define α˙ ij . Recently, Nieto-Leal, et al. [65] proposed robust RH rules for the PPS and the BS that were inspired by experimental observations of microfabric evolution of cohesive soils and guided by earlier formulations [2, 7, 12, 19, 25, 32, 39, 54, 57, 68, 77, 78]. The following RH rule applies to the PPS: 2 1 + ein I i i α˙ ij = |χη η˜ − α| ˜ ˙εv + ψ2 fη fI0 M˜ − η ˜ |˙εs | ψ1 λ−κ I0 ×
sijα I0
sign(M − α)
(58)
where ε˙ vi and ε˙ si denote volumetric and deviatoric inelastic strain increments, respectively, and 10 0 if Δη = 0 and fI0 = (59) fη = 1 if Δη > 0 1 + 100 exp(−I0 /I ) In Eq. 58, λ, κ, ein , I0 , and sijα are as previously defined. The term (1 + ein )/(λ √ − κ) is included in order to maintain similarity with Eq. 57. The quantity α˜ = 2/3 α, and the stress ratio in this space is given by η˜ = J /I , where I and J is as previously defined. The corresponding stress ratio in axisymmetric triaxial
Transp. Infrastruct. Geotech.
stress space is η = q/p , where q and p are as defined in the Definition of Invariants section. The quantities χη , ψ1 , and ψ2 are positive dimensionless model parameters. The parameter χη ensures that the correct value of K = σ3 /σ1 will indeed be reached in the simulation of constant η˜ loading [2]. Values of the model parameter ψ1 primarily control the pace of RH in the simulation of constant η˜ loading; by contrast, the value of ψ2 primarily controls the pace of RH in the simulation of shearing. The Heaviside step function fη is equal to zero during constant η˜ loading (i.e., when Δη = 0). It is p introduced so as to eliminate the dependency of α˙ on ε˙ s for constant η˜ loading. The function fI0 is the so-called logistic function [74, 75]; it predicts almost exponential growth with increasing values of the ratio c = I0 /I , which is approximately equal to the overconsolidation ratio (OCR). For higher values of this ratio, the growth will significantly slow, as the function “saturates”. Such behavior is especially needed for the accurate simulation of K0 -consolidated cohesive soils when sheared in extension under undrained conditions (i.e., for CK0 U E tests). In a similar manner, the following RH rule is proposed for the BS [65]: 2 χη 1 + ein I β˙ij = ˜ |˙εsi | ψ1 η˜ − β˜ ˙εvi + ψ2 fη fI0 M˜ − η λ−κ I0 A β
×
sij I0
sign(N − β)
(60)
˜ I , I0 , η, ˜ α, fη , fI0 , ψ1 , ψ2 , and N are as previously defined where ein , λ, κ, M, The parameter A is a constant that defines the initial or inherent anisotropy [65]. The sign(N − β) term is the only one in Eq. 60 that carries sign. It thus prevents excessive rotation of the BS in a manner very similar to that provided to the PPS by the term sign(M − α) in Eq. 58. Nieto-Leal et al. [65] give additional details pertaining to the above RH rules.
The Bounding Plastic Modulus Having defined the isotropic and rotational hardening rules associated with the GBSM, it is possible to specialize the definition of the “bounding” plastic modulus (K¯ p ) defined in Eq. 14 and required in Eq. 3 for the determination of Kp . Since the most general form of the GBSM employs both isotropic and rotational hardening, Eq. 14 becomes ∂F ∂F r1 + r2 (61) K¯ p = − ∂q1 ∂q2 For versions of the GBSM employing only isotropic hardening, only the first term in Eq. 61 is required. Recalling the discussion of isotropic hardening given in the Isotropic Hardening section, noting that q1 = ei , and using Eq. 6 gives r1 = − (1 + ein ) Rkk
(62)
Transp. Infrastruct. Geotech.
The contribution of isotropic hardening to K¯ p is thus ∂F ∂F ∂Io ∂F (1 + ein ) r1 = (1 + ein ) Rkk = − Rkk (63) [Io − IL + IL ] ∂ei ∂Io ∂ei λ−κ ∂Io where Eq. 57 has been used. The specific form of ∂F /∂Io and Rkk = ∂Q/∂ I¯ depends on the version of the GBSM that is used. Kaliakin and Nieto-Leal [48] give explicit expressions for these derivatives. For anisotropic forms of the GBSM, the term (∂F /∂q2 )r2 accounts for the RH of the BS. For a non-associative flow rule, q2 = βij so ∂F /∂q2 = ∂F /∂βij . For this case, β sij 1 + ein r2 = sign(N − β) (64) β∗ λ−κ I0 −
where, recalling (60), 2 χ I η β ∗ = ψ1 ˜ |˙εsi | η˜ − β˜ ˙εvi + ψ2 fη fI0 M˜ − η I0 A
(65)
For the anisotropic form of the GBSM employing an associative flow rule, F ≡ Q and αij ≡ βij . Now, q2 = αij so ∂F /∂q2 = ∂F /∂αij . For this case, α sij 1 + ein r2 = sign(M − α) (66) α∗ λ−κ I0 where, recalling (58), I 2 i ∗ i ˜ |χη η˜ − α| ˜ ˙εv + ψ2 fη fI0 M − η ˜ |˙εs | α = ψ1 I0
(67)
Finally, Kaliakin and Nieto-Leal [48] give explicit expressions for the derivatives ∂F /∂βij and ∂F /∂αij .
The Shape Hardening Function The positive hardening function Hˆ defines the shape of the response curves during inelastic hardening (or softening) for points within the BS [26]. It relates the plastic modulus Kp to its “bounding” value K¯ p in the manner given by Eq. 3. The majority of past functional forms used for Hˆ have the following general form [48]: Hˆ = CH fF fh fn
(68)
where CH is a constant, fF is a function of the analytical expression for the bounding surface or its derivatives, fh is a function that involves the shape hardening parameters, and fn is a function of suitable unit normal vectors. For the GBSM, (1 + ein ) (69) CH = Pa λ−κ with ein , λ, κ, and Pa are as previously defined.
Transp. Infrastruct. Geotech.
For the anisotropic form of the GBSM employing a non-associative flow rule,
2 1 2 fF = 9 F,I¯ + (70) F,J¯β 3 where, for brevity, the partial derivative with respect to an invariant is denoted with a comma followed by the symbol of the invariant as a subscript. The function fh , appearing in Eq. 68, is given by fh = h θ β z0.02 + ho 1 − z0.02 (71) where the definition of the dimensionless variable z varies depending on the form of the GBSM used. For, example, for the anisotropic form√of the GBSM employing a β non-associative flow rule, z = J β /J 1 = J β /(N˜ I1 ) = 3 3 J β R/NIo (see Fig. 4). The dimensionless quantity h θ β in Eq. 71 has the following general form: (72) h θ β = g θ β , kh hc where kh = he / hc , with he = h(−π/6) and hc = h(π/6) being the values of h associated with axisymmetric triaxial extension and compression, respectively. The dimensionless function g θ β , kh is given by Eq. 33, with θ β and kh replacing θ and k, respectively. Finally, the expression for fn appearing in Eq. 68 is a generalization of the form given by Nieto-Leal and Kaliakin [64], viz., I 1 1 fn = (73) a + sign (nI ) (|nI |) 5 2 Io where a(> 1.0) is a dimensionless model parameter [44], and nI is the component, in the I¯ direction, of the unit outward normal to the bounding surface in stress invariant space. For the anisotropic form of the GBSM with a non-associative flow rule, F,I¯ 2 2 F,I¯ + F,J¯β
nI =
(74)
If an associative flow rule is used instead, F,J¯β is replaced by F,J¯α in Eqs. 71 and 74 and in the definition for z. Finally, if an isotropic form of the GBSM is used, F,J¯β is replaced by F,J¯ in the aforementioned equations.
Quantities Describing the Initial Material State The definition of the initial state requires knowledge of the initial stress state and the maximum past consolidation pressure to which the soil has been subjected. The quantities describing the initial stress state thus consist of the total confining stresses σ11c , σ22c , σ33c , σ12c , σ13c , σ23c , the total stresses σ11p , σ22p , σ33p , σ12p , σ13p , σ23p associated with the maximum past preconsolidation pressure, and the initial excess pore fluid pressure ϑin .
Transp. Infrastruct. Geotech.
Knowing the confining and preconsolidation total stresses, the OCR is defined by /σ , where the effective preconsolidation pressure to which the soil has OCR = pin c been subjected and the initial effective confining stress are given by 1 pin = (75) σ11p + σ22p + σ33p − ϑin 3 1 (76) σ11c + σ22c + σ33c − ϑin σc = 3 The definition of the initial state for all forms of the GBSM also requires knowledge of the initial void ratio ein . Anisotropic forms of the GBSM require the specification of the inherent anisotropy. This is achieved by specifying the initial values of αij , commonly determined from the results of laboratory tests on anisotropically consolidated samples. For the special case of a transversely anisotropic (or “cross-anisotropic”) soil, the inherent (initial) values of αij are given by 2A(1 − Kin ) −A(1 − Kin ) in in ; α22 = α33 = 1 + 2Kin 1 + 2Kin in in in in in = α21 = α13 = α31 = α23 = α32 = 0
in α11 = in α12
(77) (78) (σ3 /σ1 )
where A ≤ 1.0 and Kin denotes the constant initial effective stress ratio to which a soil has been subjected in a drained axisymmetric triaxial (isotropic) or oedometer (anisotropic) test [65]. In the anisotropic form of the GBSM employing a non-associative flow rule, the rotation of the bounding surface is quantified through the symmetric anisotropic tensor βij . The initial values of βij , corresponding to a transversely anisotropic stress, will be 2(1 − Kin ) −(1 − Kin ) in in ; β22 = β33 = 1 + 2Kin 1 + 2Kin in in in in in = β21 = β13 = β31 = β23 = β32 =0
in β11 =
(79)
in β12
(80)
where Kin is as previously defined. The initial values of αij will again be given by Eqs. 77 and 78. Nieto-Leal et al. [65] give additional details pertaining to the proper definition of the initial material state in conjunction with the GBSM.
Identification of the Model Parameters The model parameters associated with all versions of the GBSM are grouped by their type as follows: the critical state parameters λ, κ, Mc , and Me , the elastic parameters G or ν, the surface configuration parameters R, Nc , and Ne , the projection center parameter C, the elastic nucleus parameter sp , the shape hardening parameters hc , he , and a, the rotational hardening parameters χη , ψ1 , and ψ2 , and the viscoplastic parameters sv , V , and n (recall Eq. 10). All of the above parameters are positive. The aforementioned parameters are associated with the most general form of the GBSM. If less general forms of the model are appropriate, then the number of active model parameters reduces. For example, if only experimental results in axisymmetric triaxial compression are to be simulated, the values of Me , Ne , and he are immaterial.
Transp. Infrastruct. Geotech.
If only time-independent results are to be simulated, the values of sv , V , and n are immaterial. If only isotropically consolidated samples are to be simulated, then values of χη , ψ1 , and ψ2 are immaterial. Finally, if an associative flow is assumed, then values of Nc and Ne are immaterial.
Numerical Implementation The constitutive relations given by Eqs. 15 and 16, along with such enhancements as local iteration, sub-stepping and radial return [36] have been incorporated into a modular system of FORTAN-77 subroutines. This system consists of a master subroutine clay UD that is supported by numerous lower-level subroutines. Implementation of the GBSM in such a manner facilitates simple and inexpensive incorporation into existing and new computer programs for the analysis of geo-structures. For purposes of assessing the characteristics of the GBSM and for simulating experimental results, a means for predicting the results of simple homogeneous tests must be made available. The Fortran-90 computer program CALBR8 [42] has been developed for this purpose. CALBR8 evaluates the constitutive relations, as applied to a homogeneous specimen, subjected to arbitrary homogeneous stress, strain and temperature histories, under either drained or undrained conditions. Although CALBR8 has been written to evaluate the time- and temperature-dependent constitutive relations, the formulation can also be used to perform time-dependent and isothermal analyses.
Assessment of Model Capabilities The analyses presented in this section serve as an assessment of the simulative and predictive capabilities of the GBSM. Undrained and Drained Shearing of Bangkok Clay Balasubramaniam and Chaudhry [4] presented the results of an extensive laboratorytesting program on soft Bangkok (Nong Ngoo Hao) clay. This program consisted of isotropically consolidated undrained (CIUC) and drained (CIDC) shearing of normally consolidated specimens subjected to various stress paths in compression. This data set is used to assess the predictive capabilities of the simplest form of the GBSM, i.e., the isotropic version, employing an associative flow rule (referred to herein as the GBSMia form of the model). Since Balasubramaniam and Chaudhry [4] did not provide sufficient time-dependent experimental data, only time-independent simulations will be generated (with no loss in computational efficiency in the clay UD subroutine system). The natural water content of Bangkok clay is 112 to 130%. The specific gravity (Gs ) for the soil has been reported to be 2.72 [6] and 2.73 [5]. The Atterberg limits for the clay are a liquid limit (wL ) of 118.0 ± 2%, a plastic limit (wP ) of 43.0 ± 2%, and a plasticity index (Ip ) of 75 ± 4%.
Transp. Infrastruct. Geotech. Table 1 Initial conditions of Bangkok clay specimens as reported by Balasubramaniam and Chaudhry [4]
pin
CIUC test
CIDC test
(kPa)
ein
with constant cell pressure ein
138
2.184
2.345
207
2.053
2.124
276
1.925
1.968
345
1.780
1.862
414
1.720
1.769
All tests performed by Balasubramaniam and Chaudhry [4] were carried out under stress-controlled conditions. Table 1 summarizes the initial void ratios associated with the various drained and undrained tests performed. Since the experimental program was limited to normally consolidated specimens, the only model parameters active are λ, κ, Mc , ν, and R. From the results of the isotropic consolidation and rebound tests reported by Balasubramaniam and Chaudhry [4], values of 0.510 and 0.090 were determined for λ and κ, respectively. Based on the reported value of 26.0◦ for the effective friction angle (φ ) in axisymmetric triaxial compression [4], a value of Mc = 1.050 for the slope of the critical state line (CSL) was computed from Mc = 6 sin φ /(3 − sin φ ). Based on the variation of Poisson’s ratio (ν) with Ip [50], a value of 0.30 was assumed. To determine a suitable value for the bounding surface configuration parameter R, the results of CIUC tests at initial confining pressures ranging from 138 to 414 kPa were used. A best fit was obtained for a value of R equal to 2.00. Using this value in conjunction with the other model parameter values given above, the time-independent undrained simulations shown in Figs. 5, 6, and 7 were generated using the GBSMia form of
CSL (M = 1.05) 138.0 kPa GBSMia 207.0 kPa GBSMia 276.0 kPa GBSMia 345.0 kPa GBSMia 414.0 kPa GBSMia
350.0
Deviator Stress (kPa)
300.0 250.0 200.0 150.0 100.0 50.0 0.0 0.0
50.0
100.0
150.0
200.0
250.0
300.0
350.0
400.0
450.0
Mean Normal Effective Stress (kPa)
Fig. 5 Comparison of simulated and experimental undrained stress paths in triaxial compression for Bangkok clay [4]
Transp. Infrastruct. Geotech.
Principal stress difference (kPa)
300.0
250.0 138.0 kPa GBSMia 207.0 kPa GBSMia 276.0 kPa GBSMia 345.0 kPa GBSMia 414.0 kPa GBSMia
200.0
150.0
100.0
50.0
0.0 0.0
2.0
4.0
6.0
8.0
10.0
12.0
14.0
16.0
18.0
20.0
Deviatoric strain (%)
Fig. 6 Comparison of simulated and experimental stress-strain response in undrained triaxial compression for Bangkok clay [4]
the model, in conjunction with the CALBR8 computer program [42]. The agreement with experimental results is seen to be quite good. With the exception of the value of Mc = 0.920 reported by Balasubramaniam and Chaudhry [4], the aforementioned model parameter values were next used to predict the drained response of normally consolidated specimens with constant cell pressure. The initial confining pressures again ranged from 138 to 414 kPa. The initial void ratio values (ein ) associated with the various confining pressures are given in the third column of Table 1.
Excess pore pressure (kPa)
300.0
250.0 138.0 kPa GBSMia 207.0 kPa GBSMia 276.0 kPa GBSMia 345.0 kPa GBSMia 414.0 kPa GBSMia
200.0
150.0
100.0
50.0
0.0 0.0
2.0
4.0
6.0
8.0
10.0
12.0
14.0
16.0
18.0
20.0
Deviatoric strain (%)
Fig. 7 Comparison of simulated and experimental excess pore pressure-strain response in undrained triaxial compression for Bangkok clay [4]
Transp. Infrastruct. Geotech.
Principal stress difference (kPa)
600.0
500.0 138.0 kPa GBSMia 207.0 kPa GBSMia 276.0 kPa GBSMia 345.0 kPa GBSMia 414.0 kPa GBSMia
400.0
300.0
200.0
100.0
0.0 0.0
4.0
8.0
12.0
16.0
20.0
Deviatoric strain (%)
Fig. 8 Comparison of simulated and experimental stress-strain response in drained compression, with constant cell pressure, for Bangkok clay [4]
Using the aforementioned model parameter values, the drained predictions shown in Figs. 8 and 9 were generated using the GBSMia. The agreement with experimental results is seen to be good, particularly in light of the fact that the model parameter values used to simulate the CIUC tests were not changed to improve the numeric results (they are thus true predictions). In general, the simulations are not “stiff” enough as compared to the experimental results. That is, for a given level of deviatoric
500.0
Principal stress difference (kPa)
450.0 400.0
138.0 kPa GBSMia 207.0 kPa GBSMia 276.0 kPa GBSMia 345.0 kPa GBSMia 414.0 kPa GBSMia
350.0 300.0 250.0 200.0 150.0 100.0 50.0 0.0 0.0
3.0
6.0
9.0
12.0
Volumetric strain (%)
Fig. 9 Comparison of simulated and experimental stress-strain response in drained compression, with constant cell pressure, for Bangkok clay [4]
Transp. Infrastruct. Geotech.
strain, the simulated principal stress difference values are lower than the experimental ones. Undrained Shearing of Cardiff Kaolin Stipho [73] conducted a series of axisymmetric undrained triaxial compression (CIUC) and extension (CIUE) tests on remolded laboratory prepared specimens of isotropically and anisotropically consolidated Cardiff Kaolin (Gs = 2.65, wL = 52.0%, wP = 26.0%). The kaolin specimens were consolidated anisotropically before shearing. During the anisotropic consolidation and swelling to a prescribed OCR, the principal stress ratio (Kin = σ3 /σ1 ) was kept constant. For the present model assessment, the specific case of Kin = 0.667 is considered. Table 2 summarizes the initial conditions of the CIUC and CIUE specimens immediately before shearing for this case. In the past, these test results have been used to verify the predictive capabilities of various isotropic and anisotropic constitutive models [2, 7, 55, 57]. This data set is used to assess the predictive capabilities of the anisotropic form of the GBSM employing an associative flow rule (referred to herein as the GBSMaa form of the model). Since Stipho [73] did not provide sufficient time-dependent experimental data, only time-independent simulations will agin be generated. Values for the critical state model parameters λ = 0.140 and κ = 0.050 were obtained directly from Stipho [73], and are consistent with values used by others for the same soil [7]. The effective friction angle (φ ) was found to be 29.5o . The values of Mc and Me were thus determined as follows: Mc =
6 sin φ = 1.178 ; 3 − sin φ
Me =
6 sin φ = 0.846 3 + sin φ
(81)
A value of ν = 0.26 was assumed, based on the variation of Poisson’s ratio (ν) with Ip [50]. A suitable value for the bounding surface configuration parameter R = 2.60 was determined by matching the undrained stress paths for normally consolidated samples from a set of isotropically consolidated tests (i.e., for Kin = 1.000) on the same soil. The value for the projection center parameter C = 0.28 was determined by matching the shapes of the undrained stress paths for lightly overconsolidated (OCR = 2.0 and 4.0) samples from the same isotropic series of tests. Values for the shape hardening parameters hc = 4.00 and he = 5.60 were obtained by matching the deviatoric stress (q)-axial strain response for the same lightly over-
Table 2 Initial conditions of anisotropically consolidated Cardiff Kaolin specimens with Kin = 0.667 [73]
pin
σc
σ1c
σ3c
(kPa)
(kPa)
(kPa)
(kPa)
1.0
206.9
206.9
266.0
177.4
2.0
206.9
103.5
133.0
88.7
1.057
4.0
206.9
51.7
66.5
44.4
1.068
8.0
551.8
69.0
88.7
59.1
0.988
OCR
ein
1.004
Transp. Infrastruct. Geotech.
consolidated samples. A value of a = 1.20 was then determined by matching the deviatoric stress-axial strain response for samples at OCR = 8.0. Finally, a value of sp = 1.50 for the elastic nucleus parameter was required so as to give slightly “stiffer” response for the initial stages of response for the overconsolidated samples. Although the above model parameter values were determined by matching the experimental data for isotropically consolidated samples, they likewise apply to the subsequent simulation of anisotropically consolidated samples. Finally, values for the rotational hardening parameters ψ1 , χη , and ψ2 were next determined following the procedure outlined by Nieto-Leal et al. [65], giving ψ1 = 300, χη = 2.74, and ψ2 = 10.0, for a value of A equal to 0.64. Using (77) gives in = 0.182, and α in = α in = −0.091. The time-independent model simulations α11 22 33 and predictions were obtained using the above parameter values in conjunction with the anisotropic form of the GBSM with an associative flow rule (i.e., the GBSMaa version of the model), in conjunction with the CALBR8 computer program [42]. Figures 10 and 11 compare the simulated normalized undrained stress paths (i.e., vs. q/p ) with experimental values for the CIUC and CIUE tests, respecp /pin in tively. Figures 12 and 13 compare the simulated variation of normalized deviator ) with axial strain with experimental values. Finally, Figures. 14 and 15 stress (q/pin compare the simulated variation of the normalized excess pore pressure with axial strain with experimental values. Overall, the simulations generated using the GBSMaa version of the model are in good agreement with the experimental results. To illustrate the effect that the inclusion of anisotropy and RH has on the above results, the response of Cardiff Kaolin was simulated using the GBSMia form of the model. Although the initial conditions for each test are unchanged form the values listed in Table 2, the orientation of the BS is now symmetric with respect to the I in = α in = α in = 0.0). In Figs. 10, 11, 12, 13, 14, and 15, the dashed axis (i.e., α11 22 33
Fig. 10 Comparison of simulated and experimental normalized undrained stress paths in triaxial compression for data of Stipho [73]. No experimental results available for OCR = 8.0
Transp. Infrastruct. Geotech.
Fig. 11 Comparison of simulated and experimental normalized undrained stress paths in triaxial extension for data of Stipho [73]. No experimental results available for OCR = 8.0
lines depict the results obtained using the GBSMia form of the model. The absence of anisotropy has a particularly pronounced effect on the CIUE simulations, especially the normalized undrained stress paths (Fig. 11) and excess pore pressure response (Fig. 15). Clearly, the inclusion of anisotropy significantly increases the accuracy of the simulations.
Fig. 12 Comparison of simulated and experimental normalized stress-strain response in triaxial compression for data of Stipho [73]
Transp. Infrastruct. Geotech.
Fig. 13 Comparison of simulated and experimental normalized stress-strain response in triaxial extension for data of Stipho [73]
Simulation of the Effect of Sample Age on Undrained Creep Shen et al. [71] performed extensive undrained axisymmetric triaxial tests on undisturbed specimens of San Francisco Bay Mud, a soft normally consolidated silty clay (Gs = 2.65, wL = 88%, wP = 43%). The testing program was divided into two
Fig. 14 Comparison of simulated and experimental normalized excess pore pressure-strain response in triaxial compression for data of Stipho [73]
Transp. Infrastruct. Geotech.
Fig. 15 Comparison of simulated and experimental normalized excess pore pressure-strain response in triaxial extension for data of Stipho [73]
parts. The first part consisted of nine undrained axisymmetric triaxial strength tests. The effective isotropic consolidation pressure (σc = 196.2, 294.3, and 392.4 kPa) and the period of consolidation (Tc = 1, 3, and 7 days) were varied in these tests. At the end of their respective periods of consolidation, undrained, stress-controlled triaxial tests with pore pressure measurements were performed. The second part of the testing program also consisted of nine specimens. The same σc and Tc values were employed as in the strength tests. Following their respective periods of consolidation, a strength test was not, however, performed. Instead, the drainage valve was closed and an isotropic state of total stress was maintained constant for a period of 8 days. During these undrained creep conditions, the increase in pore pressure was measured continuously and no axial strain was detected. One additional undrained creep test was performed with σc = 196.2 kPa and Tc = 14 days. Figure 16 shows experimental results for Tc = 1 and 7 days, which are depicted by discrete symbols. In all the undrained creep tests, following an initial period of increase, the rate of pore pressure build-up decreased with time. For any given Tc , the amount of pore pressure build-up increased with σc . Also, for any given σc , the pore pressure build-up decreased with increasing Tc , i.e., with the “age” of the specimen. This is explained by the fact that “older” specimens afford greater time for particle rearrangement under drained conditions and thus greater degrees of secondary compression. The resulting reduction in potential particle rearrangement gives “older” specimens a lower potential for subsequent pore pressure build-up under undrained conditions. Figure 17 shows that if, for a given period of consolidation the pore pressure is normalized by σc , the resulting experimental points fall in fairly narrow bands. Since all specimens were isotropically consolidated, this data set is used to assess the time-dependent predictive capabilities of the GBSMia form of the
Transp. Infrastruct. Geotech.
Fig. 16 Variation of pore pressure with time during undrained creep of San Francisco Bay mud for Tc = 1 and 7 days [71]
model. Corresponding to the continuously increasing pore pressure measured during undrained creep is a continuously decreasing state of effective stress. For normally consolidated specimens subjected to isotropic states of total stress, Kaliakin and Dafalias [45] showed that, for the GBSMia form of the model, the mean normal effective stress decreases according to p˙ = −3 Φ K
∂F ∂ I¯
(82)
Fig. 17 Variation of normalized pore pressure with time during undrained creep of San Francisco Bay mud for Tc = 1 and 7 days [71]
Transp. Infrastruct. Geotech.
where K is the elastic bulk modulus, and the overstress function is given by 1 Φ= V
Io I
J exp N˜ I
Δσˆ
n
(83)
The inclusion of the (Io /I ) term in Eq. 83, which represents a new development in the definition of Φ, better accounts for the “age” of a specimen. The normalized overstress (Δσˆ ) appearing in Eq. 83 is as defined by Eq. 9. In Eq. 83 V and n are positive model parameters. The former influences the initial portion of the time-dependent response, while the latter controls the magnitude of the ultimate (“steady-state”) values that are attained. Whereas n is dimensionless, V has units of stress-time (i.e., FL−2 t). The initial void ratios (ein ), as well as the values of the critical state parameters λ = 0.369, κ = 0.100, and Mc = 1.49 were taken directly from the work of Shen et al. [71]. Based on the variation of Poisson’s ratio (ν) with Ip [50], a value of 0.285 was assumed. The shape parameter R = 2.80 was determined by matching the experimental undrained stress paths for the strength test associated with Tc = 1 day; the loading rate employed in the simulations was identical to that used by Shen et al. [71]. Values for the model parameters C = 0.20 and sv = 2.40 were determined by matching the changes in p measured during undrained creep. Values for the final two viscoplastic parameters, V and n, were determined by matching numerical predictions with experimental pore pressure time histories. In the course of determining the aforementioned values for V and n, the value of V giving the best agreement with experimental results was found to be constant for a given value of Tc . In particular, for Tc = 1 day, V = 9.0 × 106 kPa-min; for Tc = 3 days, V = 1.0 × 107 kPa-min; for Tc = 7 days, V = 5.0 × 107 kPa-min; and, for Tc = 14 days, V = 9.0 × 107 kPa-min. Although the values of n appear to increase with σc , the results are inconclusive - additional data sets need to be analyzed to investigate this issue in greater detail. The simulations of undrained triaxial creep obtained using the aforementioned parameter values are shown in Figs. 16 and 17 in the form of continuous and dashed curves. The overall agreement between model simulations and experimental results is quite good and shows the model to be capable of simulating the effects of sample age (under drained conditions) and its effect on subsequent undrained creep response. Qualitative Simulation of Drained Heating The final example qualitatively simulates the response of hypothetical normally consolidated and overconsolidated cohesive soil samples subjected only to changes in temperature under drained conditions. Experimental results [82] show that as the OCR increases, the volumetric strains resulting from a temperature cycle progressively become less compressive, eventually changing to dilatational. This response is simulated through the reversible and thermoplastic strain increments given by Eqs. 21 and 22, respectively. Figure 18 shows the qualitative response for OCR = 1, 2, 4, and 8.
Transp. Infrastruct. Geotech.
Fig. 18 Qualitative simulation of a temperature cycle under drained conditions
Conclusion This paper presented details pertaining to the Generalized Bounding Surface Model (GBSM) for cohesive soils. This model synthesizes many previous bounding surface models for cohesive soils, improves many aspects of these models, and expands the formulation to non-isothermal conditions. The simulative and predictive capabilities of the GBSM were assessed. Based on comparisons with experimental results, several forms of the GBSM were shown to realistically simulate the time-independent and time-dependent response of different cohesive soils. Finally, a temperature cycle under drained conditions at different levels of overconsolidation was qualitatively simulated. Acknowledgements The first author wishes to acknowledge past interactions with the research group of Professor H. I. Ling in the Department of Civil Engineering & Engineering Mechanics at Columbia University. The second author’s graduate studies were supported by the Fulbright Colombia-Colciencias Scholarship and Universidad Militar Nueva Granada. The third author’s graduate studies were supported by the University of Delaware’s Department of Civil & Environmental Engineering. This support is gratefully acknowledged.
References 1. Al-Shamrani, M.A., Sture, S.: A time-dependent bounding surface model for anisotropic cohesive soils. Soils Found. 38(1), 61–76 (1998) 2. Anandarajah, A., Dafalias, Y.F.: Bounding surface plasticity III: Application to anisotropic cohesive soils. J. Eng. Mech. ASCE, 112(12), 1292–1318 (1986) 3. Baker, R., Desai, C.S.: Induced anisotropy during plastic straining. Int. J. Numer. Anal. Methods Geomech. 8(2), 167–185 (1984) 4. Balasubramaniam, A.S., Chaudhry, A.R.: Deformation and strength characteristics of soft bangkok clay. J. Geotech. Eng. Div. ASCE 104(GT9), 1153–1167 (1978) 5. Balasubramaniam, A.S., Li, Y.G.: Stress-strain behaviour of soft bangkok clay in triaxial extension. In: Proceedings of the International Symposium on Soft Clay, pp. 105–116. Bangkok (1977)
Transp. Infrastruct. Geotech. 6. Balasubramaniam, A.S.: Waheed-Uddin: Deformation characteristics of weathered bangkok clay in triaxial extension. G˙eotechnique 27(1), 75–92 (1977) 7. Banerjee, P.K., Yousif, N.B.: A plasticity model for the mechanical behavior of anisotropically consolidated clay. Int. J. Numer. Analytic. Methods Geomech. 10(5), 521–541 (1986). https://doi.org/10.1002/nag.1610100505 8. Barden, L.: Primary and secondary consolidation of clay and peat. G´eotechnique 18(1), 1–24 (1968) 9. Casagrande, A., Carrillo, N.: Shear failure of anisotropic materials. Proc. Boston Soc. Civil Eng. 31, 74–87 (1944) 10. Cheney, J.A., Arulanandan, K.: Stress relaxation in sand and clay. In: Constitutive Equations of Soils, 9th International Conference on Soil Mechanics and Foundation Engineering, pp. 279–285. Tokyo (1977) 11. Christensen, R.W., Wu, T.H.: Analysis of clay deformation as a rate process. J. Soil Mech. Found. Div. ASCE 90(SM6), 125–157 (1964) 12. Coombs, W., Crouch, R., Augarde, C.: A unique critical state two-surface hyperplasticity model for fine-grained particulate media. J. Mech. Phys. Solids 61(1), 175–189 (2013) 13. Crouch, R.S., Wolf, J.P.: A unified 3-d anisotropic modular elliptic bounding surface model for soil. In: Pande, G.N., Pietruszczak S. (eds.) Proceedings of the 4th International Conference on Numerical Methods in Geomechanics, pp. 137–147 (1992) 14. Crouch, R.S., Wolf, J.P.: Unified 3-d critical state bounding-surface plasticity model for soils incorporating continuos plastic loading under cyclic paths i: Constitutive relations. Int. J. Numer. Anal. Methods Geomech. 18(11), 735–758 (1994) 15. Crouch, R.S., Wolf, J.P.: On a three-dimensional anisotropic plastic model for soil. G˙eotechnique 45(2), 301–305 (1995) 16. Crouch, R.S., Wolf, J.P., Dafalias, Y.F.: Unified critical state bounding surface plasticity models for soil. J. Eng. Mech. ASCE 120(11), 2251–2270 (1994) 17. Cui, Y.J., Sultan, N., Delage, P.: A thermomechanical model for saturated clays. Canad. Geotech. J. 37(3), 607–620 (2000). https://doi.org/10.1139/cgj-37-3-607 18. Dafalias, Y.: Discussion: Constitutive equations of elastoplastic materials with elastic-plastic transition (Hashiguchi, K., 1980, J. Appl. Mech., ASME, 47: 266–272). J. Appl. Mech. 48(1), 211–212 (1981) 19. Dafalias, Y., Taiebat, M.: Anatomy of rotational hardening in clay plasticity. G´eotechnique 63(16), 1406–1418 (2013). https://doi.org/10.1680/geot.12.P.197 20. Dafalias, Y.F.: A bounding surface plasticity model. In: Proceedings of the Seventh Canadian Congress of Applied Mechanics, pp. 89-90. Sherbrooke (1979) 21. Dafalias, Y.F.: A model for soil behavior under monotonic and cyclic loading conditions. In: Transactions of the Fifth International Conference on SMIRT, pp. 1–9. Berlin (1979) 22. Dafalias, Y.F.: The concept and application of the bounding surface in plasticity theory. In: Hult, J., Lemaitre, J. (eds.) IUTAM Symposium on Physical Non-Linearities in Structural Analysis, pp. 56–63. Springer, Berlin (1981) 23. Dafalias, Y.F.: Bounding surface elastoplasticity-viscoplasticity for particulate cohesive media. In: Vermeer, P.A., Luger, H.J. (eds.) Deformation and Failure of Granular Materials, IUTAM Symposium on Deformation and Failure of Granular Materials, pp. 97–107. A. A. Balkema, Rotterdam (1982) 24. Dafalias, Y.F.: On rate dependence and anisotropy in soil constitutive modeling. In: Gudehus, G., Darve, F., Vardoulakis, I. (eds.) Results of the International Workshop on Constitutive Relations for Soils, pp. 457–462. A. A. Balkema, Rotterdam (1982) 25. Dafalias, Y.F.: An anisotropic critical state soil plasticity model. Mech. Res. Commun. 13(6), 341–347 (1986) 26. Dafalias, Y.F.: Bounding surface plasticity. I: Mathematical foundation and the concept of hypoplasticity. J. Eng. Mech. ASCE, 112(90), 966–987 (1986) 27. Dafalias, Y.F., Herrmann, L.R.: A bounding surface soil plasticity model. In: Pande, G.N., Zienkiewicz, O.C. (eds.) Proceedings of the International Symposium on Soils Under Cyclic and Transient Loading, pp. 335–345. A. A. Balkema, Rotterdam (1980) 28. Dafalias, Y.F., Herrmann, L.R.: Bounding surface formulation of soil plasticity. In: Pande, G.N., Zienkiewicz, O.C. (eds.) Soil Mechanics – Transient and Cyclic Loads, chap. 10, pp. 253–282. Wiley, Chichester (1982)
Transp. Infrastruct. Geotech. 29. Dafalias, Y.F., Herrmann, L.R.: A generalized bounding surface constitutive model for clays. In: Yong, R.N., Selig, E.T. (eds.) Application of Plasticity and Generalized Stress-Strain in Geotechnical Engineering, pp. 78–95. ASCE, New York (1982) 30. Dafalias, Y.F., Herrmann, L.R.: Bounding surface plasticity II: Application to isotropic cohesive soils. J. Eng. Mech. ASCE 112(12), 1263–1291 (1986) 31. Dafalias, Y.F., Herrmann, L.R., De Natale, J.S.: The bounding surface plasticity model for isotropic cohesive soils and its application at the grenoble workshop. In: Gudehus, G., Darve, F., Vardoulakis, I. (eds.) Results of the International Workshop on Constitutive Relations for Soils, pp. 273–287. A. A. Balkema, Rotterdam (1982) 32. Dafalias, Y.F., Manzari, M.T., Papadimitriou, A.G.: SANICLAY: Simple anisotropic clay plasticity model. Int. J. Numer. Anal. Methods Geomech. 30(12), 1231–1257 (2006) 33. Gudehus, G.: Elastoplastische stoffgleichungen f˙ur trockenen sand. Ingenieur-Archiv 42(3), 151–169 (1973) 34. Hashiguchi, K.: Subloading surface model in unconventional plasticity. Int. J. Solids Struct. 25(8), 917–945 (1989) 35. Hashiguchi, K., Ueno, M.: Elasto-plastic constitutive laws of granular materials. In: Murayama, S., Schofield, A.N. (eds.) Constitutive Equations of Soils, pp. 73–82. Tokyo, Japan (1977) 36. Herrmann, L.R., Kaliakin, V.N., Shen, C.K., Mish, K.D., Zhu, Z.Y.: Numerical implementation of a plasticity model for cohesive soils. J. Eng. Mech. ASCE 113(4), 500–519 (1987) 37. Hong, P.Y., Pereira, J.M., Cui, Y.J., Tang, A.M.: A two-surface thermomechanical model for saturated clays. Int. J. Numer. Anal. Methods Geomech. 40(7), 1059–1080 (2016) 38. Jiang, J., Ling, H.I., Kaliakin, V.N.: A time-dependent anisotropic bounding surface model for clays. In: Proceedings of the Engineering Mechanics 2011 (EMI 2011), pp. CD–ROM. ASCE (2011) 39. Jiang, J., Ling, H.I., Kaliakin, V.N.: An associative and non-associative anisotropic bounding surface model for clay. J. Appl. Mech. ASME, 79(3), https://doi.org/10.1115/1.4005,958. (Jim Rice Special Edition) (2012) 40. Jiang, J., Ling, H.I., Kaliakin, V.N.: Simulation of natural pisa clay using an enhanced anisotropic elastoplastic bounding surface model. J. Appl. Mech. ASME, 80(2), https://doi.org./10.1115/1.4007, 964 (2013) 41. Kaliakin, V.N.: Bounding surface elastoplasticity-viscoplasticity for clays. Ph.D. thesis, University of California Davis,, Davis, CA (1985) 42. Kaliakin, V.N.: Calbr8, a simple computer program for assessing the idiosyncrasies of various constitutive models used to characterize soils. Tech Rep. 92–1. University of Delaware, Department of Civil Engineering, Newark, DE (1992) 43. Kaliakin, V.N.: Numerical implementation and solution strategies for a thermo-elastoplasticviscoplastic model for cohesive soils. Comput. Syst. Eng. 5(2), 203–214 (1994) 44. Kaliakin, V.N., Dafalias, Y.F.: Simplifications to the bounding surface model for cohesive soils. Int. J. Numer. Anal. Methods Geomech. 13(1), 91–100 (1989) 45. Kaliakin, V.N., Dafalias, Y.F.: Theoretical aspects of the elastoplastic-viscoplastic bounding surface model for cohesive soils. Soils Found. 30(3), 11–24 (1990). Japanese Society of Soil Mechanics and Foundation Engineering 46. Kaliakin, V.N., Dafalias, Y.F.: Verification of the elastoplastic-viscoplastic bounding surface model for cohesive soils. Soils Found. 30(3), 25–36 (1990). Japanese Society of Soil Mechanics and Foundation Engineering 47. Kaliakin, V.N., Nieto-Leal, A.: Towards a generalized bounding surface model for cohesive soils. In: Hellmich, C., Pichler, C.B., Adam, C. (eds.) Proceedings of the 5th Biot Conference on Poromechanics (BIOT-5), p. CD–ROM, Vienna (2013) 48. Kaliakin, V.N., Nieto-Leal, A.: Details pertaining to the generalized bounding surface model for cohesive soils: Revised & expanded. Tech. rep., Department of Civil and Environmental Engineering, University of Delaware. Newark, DE. USA (2017) 49. Kirkgard, M.M., Lade, P.V.: Anisotropy of normally consolidated San Francisco, bay mud. ASTM Geotech. Testing J. 14(3), 231–246 (1991) 50. Lade, P.V.: Calibration of the single hardening constitutive model for clays. In: Proceedings of the 11th International Conference of the Association for Computer Methods and Advances in Geomechanics (IACMAG), pp. CD–ROM. Turin (2005)
Transp. Infrastruct. Geotech. 51. Laloui, L., Cekerevac, C.: Thermo-plasticity of clays: An isotropic yield mechanism. Comput. Geotech. 30(8), 649–660 (2003) 52. Laloui, L., Cekerevac, C.: Non-isothermal plasticity model for cyclic behaviour of soils. Int. J. Numer. Anal. Methods Geomech. 32(5), 437–460 (2008) 53. Laloui, L., Franc¸ois, B.: ACMEG-T: Soil thermoplasticity model. J. Eng. Mech. 135(9), 932–944 (2009) 54. Liang, R.Y., Ma, F.: Anisotropic plasticity model for undrained cyclic behavior of clays. I: Theory. J. Geotech. Eng. ASCE 118(2), 229–245 (1992) 55. Liang, R.Y., Ma, F.: Anisotropic plasticity model for undrained cyclic behavior of clays. ii: Verification. J. Geotech. Eng. ASCE 118(2), 246–265 (1992) 56. Ling, H.I., Yue, D., Kaliakin, V.N.: Geosynthetic-reinforced containment dike constructed over soft foundation: Numerical analysis. In: Ling, H.I., Leshchinsky, D., Tatsuoka, F. (eds.) Reinforced Soil Engineering: Advances in Research and Practice, pp. 283–295. Marcel Dekker Inc., New York (2003) 57. Ling, H.I., Yue, D., Kaliakin, V.N., Themelis, N.J.: An anisotropic elasto-plastic bounding surface model for cohesive soils. J. Eng. Mech. ASCE 128(7), 748–758 (2002) 58. Modaressi, H., Laloui, L.: A thermo-viscoplastic constitutive model for clays. Int. J. Numer. Anal. Methods Geomech. 21(5), 313–335 (1997) 59. Mrˇoz, Z., Norris, V.A., Zienkiewicz, O.C.: An anisotropic hardening model for soils and its application to cyclic loading. Int. J. Numer. Anal. Methods Geomech. 2(3), 203–221 (1978) 60. Mrˇoz, Z., Norris, V.A., Zienkiewicz, O.C.: Application of an anisotropic hardening model in the analysis of elasto-plastic deformation of soils. G˙eotechnique 29(1), 1–34 (1979) 61. Mrˇoz, Z., Norris, V.A., Zienkiewicz, O.C.: Elastoplastic and viscoplastic constitutive models for soils with application to cyclic loading. In: Pande, G.N., Zienkiewicz, O.C. (eds.) Soil Mechanics Transient and Cyclic Loads, chap. 8, pp. 173–217. Wiley, Chichester (1982) 62. Nieto-Leal, A.: Generalized bounding surface model for cohesive soils: A novel formulation for monotonic and cyclic loading. Ph.D. thesis, University of Delaware Department of Civil and Environmental Engineering (2016) 63. Nieto-Leal, A., Kaliakin, V.N.: On soil yielding and suitable choices for yield and bounding surfaces. Tech. rep., University of Delaware, DE. Department of Civil and Environmental Engineering (2013) 64. Nieto-Leal, A., Kaliakin, V.N.: Improved shape hardening function for bounding surface model for cohesive soils. J. Rock Mech. Geotech. Eng. 6(3), 327–337 (2014) 65. Nieto-Leal, A.N., Kaliakin, V.N., Mashayekhi, M.: Improved RH rule for cohesive soils and inherent anisotropy. Int. J. Numer. Anal. Methods Geomech. 42(3), 469–487 (2018) 66. Perzyna, P.: The constitutive equations for rate sensitive plastic materials. Q. Appl. Math. 20, 321–332 (1963) 67. Perzyna, P.: Fundamental problems in viscoplasticity. Adv. Appl. Mech. 9, 243–377 (1966) 68. Pestana, J.M., Whittle, A.J.: Formulation of a unified constitutive model for clays and sands. Int. J. Numer. Analytic. Methods Geomech. 23(12), 1215–1243 (1999). https://doi.org/10.1002/(SICI)10969853(199910)23:12<1215::AID-NAG29>3.0.CO;2-F 69. Roscoe, K.H., Burland, J.B.: On the generalized stress-strain behaviour of wet clay. In: Heyman, J., Leckie, F.A. (eds.) Engineering Plasticity, pp. 535–609. Cambridge University Press, London (1968) 70. Schofield, A.N., Wroth, C.P.: Critical State Soil Mechanics. McGraw-Hill Book Co. Inc., London (1968) 71. Shen, C.K., Arulanandan, K., Smith, W.S.: Secondary consolidation and strength of a clay. J. Soil Mech. Found. Div. ASCE 99(SM1), 95–110 (1973) 72. Sheng, D., Sloan, S.W., Yu, H.S.: Aspects of finite element implementation of critical state models. Comput. Mech. 26, 185–196 (2000) 73. Stipho, A.S.: Theoretical and experimental investigation of the behavior of anisotropically consolidated kaolin. Ph.D. thesis, University College. Cardiff, England (1978) 74. Verhulst, P.: Recherches math´ematiques sur la loi d’accroissement de la population. Nouv. m´em. de l’Academie Royale des Sci. et Belles-Lettres de Bruxelles 18, 1–41 (1845) 75. Verhulst, P.: Deuxi`eme m´emoire sur la loi d’accroissement de la population. M´em. de l’Academie Royale des Sci., des Lettres et des Beaux-Arts de Belgique 20, 1–32 (1847) 76. Voyiadjis, G.Z., Kim, D.: Finite element analysis of the piezocone test in cohesive soils using an elastoplastic-viscoplastic model and updated lagrangian formulation. Int. J. Plast. 19(2), 253–280 (2003)
Transp. Infrastruct. Geotech. 77. Wheeler, S., Naatanen, A., Karstunen, M., Lojander, M.: An anisotropic elastoplastic model for soft clays. Canadian Geotechnical J. 40(2), 403–418 (2003). https://doi.org/10.1139/T02-119 78. Whittle, A., Kavvadas, M.: Formulation of mit-e3 constitutive model for overconsolidated clays. J. Geotech. Eng. 120(1), 173–198 (1994) 79. Wong, P.K.K., Mitchell, R.J.: Yielding and plastic flow of sensitive cemented clay. G˙eotechnique 25(4), 763–782 (1975) 80. Yue, D., Ling, H.I., Kaliakin, V.N., Themelis, N.J.: Formulation and calibration of an anisotropic bounding surface model for clay. In: Ling, H.I., Anandarajah, A., Manzari, M.T., Kaliakin, V.N., Smyth, A. (eds.) Constitutive Modeling of Geomaterials: Selected Contributions from Frank L. DiMaggio Symposium, pp. 137–144. CRC Press, Boca Raton (2003) 81. Zhou, C., Fong, K., Ng, C.: A new bounding surface model for thermal cyclic behaviour. Int. J. Numer. Anal. Methods Geomech. 41(16), 1656–1666 (2017) 82. Zhou, C., Ng, C.: A thermomechanical model for saturated soil at small and large strains. Can. Geotech. J. 52(January), 1–10 (2015) 83. Zienkiewicz, O.C., Pande, G.N.: Some useful forms of isotropic yield surfaces for soil and rock mechanics. In: Gudehus, G. (ed.) Finite Elements in Geomechanics, chap. 5, pp. 179–190. Wiley, London (1977)
