Unique critical state single-surface anisotropic hyperplasticity

XII International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XII E. O˜nate, D.R.J. Owen, D. Peric and B. Su´arez (Eds) UNIQUE CRITICAL STATE SINGLE-SURFACE ANISOTROPIC HYPERPLASTICITY WILLIAM M. COOMBS∗ ∗School of Engineering and Computing Sciences Durham University Science Site, South Road, Durham, DH1 3LE, UK. e-mail: w.m.coom[email protected], web page: http://www.dur.ac.uk Key words: Computational Plasticity, Critical State Soil Mechanics, hyperplasticity Abstract. This paper presents the theoretical development and algorithmic implementation of a single surface anisotropic hyperplasticity model. The model extends the isotropic family of models developed by Coombs and Crouch (2011) through (i) introducing anisotropic shearing into the yield surface and (ii) using a more physically realistic pressure sensitive elastic free energy function. This model overcomes the difficulty of determining the constants of the isotropic two-parameter surface by analytically relating them to a single experimentally measurable physical quantity, namely the normalised hydrostatic position of the Critical State. This link results in a unique Critical State surface, invariant of the level of anisotropy inherent in the yield envelope. The model is compared with experimental data on Lower Cromer Till and contrasted against the SANIclay model. 1 INTRODUCTION A significant number of constitutive models intended to capture the anisotropic behaviour of fine grain particulate media (such as clays) have been proposed previously in the literature. The majority of these models have their roots within the classical framework of Critical State Soil Mechanics (CSSM) developed in Cambridge in the 1950s and 60s by Roscoe and co-workers[11], founded on the earlier work of Casagrande[1]. Many of the modifications to the original CSSM conceptual framework were motivated, not by deep insights into the underlying physics, but rather through a wish to improve curve-fits. This paper presents the theoretical development and algorithmic implementation of a single surface anisotropic hyperplasticity model constructed within a CSSM framework. The layout of the paper is as follows. Section 2 presents the theoretical development of the anisotropic single surface model including: (i) elasticity relationship, (ii) dissipation, (iii) Lode angle dependency (LAD), (iv) parameters controlling the shape of the yield surface, 1 Unique Critical State single-surface anisotropic hyperplasticity 1360 William M. Coombs (v) isotropic expansion or contraction and (vi) development of anisotropic shearing. The model’s stress integration is detailed in Section 3. Section 4 compares the proposed model with experimental data on Lower Cromer Till[9] and with the SANIclay model[8]. Brief conclusions are drawn in Section 5. 2 ANISOTROPIC CONSTITUTIVE FORMULATION This section presents the theoretical development of the single surface anisotropic model and draws together the key equations required for its implementation in Section 3. Section 2.1 and 2.2 describe the free-energy and dissipation functions, respectively, which are used to develop the stress versus elastic strain relationship, the yield function and the direction of inelastic straining. Section 2.3 discusses the implementation of a LAD in the model and Section 2.4 derives a relationship for the yield surface shape parameters based on the level of induced anisotropy. The model’s hardening laws are detailed in Sections 2.5 and 2.6. 2.1 Elastic free-energy function Here we use an elastic free energy function that provides pressure sensitive bulk and shear moduli[10] Ψ1=κprexpΩ+Gγe ijγe ij,where Ω = εe v−εe v0 κand G=G0+αeprexpΩ.(1) The elastic strain measures are given by εe v=εe ii and γe ij =εe ij −εe vδij/3, where δij is the Kronecker delta tensor. κis the bi-logarithmic elastic compressibility index (the gradient of the drained unloading line in the bi-logarithmic void ratio versus hydrostatic pressure plane), G0is the constant component of the shear modulus, pris the reference pressure, εe v0is the elastic volumetric strain at that reference pressure and αeis a dimensionless variable that controls relative sizes of the moduli. Taking the partial derivative of (1) with respect to the elastic strain, the Cauchy stress is given by σij =pδij +2Gγe ij where p=prexpΩ1+αe κγe ijγe ij. (2) Taking the second derivative of the free energy function with respect to elastic strain, the non-linear elastic stiffness matrix subsequently follows as De ijkl =p κ−2G 3δijδkl +2GIijkl+2prαe κexp(Ω)δijγe kl +γe ijδkl,(3) where Iijkl is a fourth order identity tensor. It is important to note that the form of elasticity presented here includes stress-induced anisotropy. 2 1361 William M. Coombs 2.2 Dissipation Starting from the following dissipation function ˙ Φ=(˙εp v+βij ˙γp ij)2A2+(˙εp γB)2,where ˙εp v=˙εp ii,and ˙γp ij =˙εp ij −˙εp vδij/3.(4) The stress like quantities are given by A= (1 −γ)p+γpc/2 and B=¯ρ(θ)M(1 −α)p+ αγpc/2, where p=σii/3, q=√sijsij and sij =σij −pδij.βij links the volumetric and deviatoric dissipation components, pcand Mcontrol the size and the axis-ratio of the yield surface, αand γcontrol the shape of the yield surface in the p-qplane and ¯ρ(θ) controls the deviatoric section. (4) was first introduced (in triaxial stress space) by Collins and Hilder[3] as an extension to the isotropic family of CS models. However, the model was only presented conceptually, and limited to the axi-symmetric triaxial case. Following the standard procedure (as given by Coombs and Crouch[4]), we obtain the dimensionless anisotropic yield surface in true stress space as f=γ(2 −γ)(¯p−1) ¯ B2+rβ ijrβ ij ¯p¯ A2=0,(5) where rβ ij =rij −βij,rij =sij/p,¯p=p/pcand ¯ A= (1 −γ)¯p+γ/2 and ¯ B=¯ρ(θ)M(1 −α)¯p+αγ/2.(6) The direction of plastic flow in true stress space similarly follows as (g,σ)ij =2/3¯ B2(¯p−γ/2) −¯ A2¯p(rβ klβkl)δij +2¯ A2¯pr β ij.(7) From (5) it is apparent that introducing a cross-coupling in the rate of dissipation function results in the yield surface being sheared off the hydrostatic axis, where βij is a second order, traceless (deviatoric), tensor measure of this inclination. If βij = 0 we recover an isotropic yield surface, with the ellipsoid’s major axis coincident with the hydrostatic axis. 2.3 Lode angle dependency The yield surface, (5), includes a dependency on the Lode angle, θ, through the normalised deviatoric yield radius, ¯ρ(θ), in ¯ B. Here, the model is presented with a Willam- Warnke (W-W)[13] LAD that can be expressed as[12] ¯ρ(θ)=a1C+√2a1C2+a2 2a1C2+1 ∈[¯ρe,1] where a1=2(1 −¯ρ2 e) (2¯ρe−1)2,a 2=5¯ρ2 e−4¯ρe (2¯ρe−1)2(8) and C= cos(π/6−θ). This W-W LAD is based on a local measure of the Lode angle, θ, from the major, βij, axis of the surface. This is achieved by measuring the second and third deviatoric stress invariants (J2and J3) based on the deviatoric distance from 3 1362 William M. Coombs the axis of anisotropy rather than the standard deviatoric measure sij. This definition ensures convexity of the yield surface and in this case the Lode angle is calculated from θ=1 3arcsin −3√3 2 J3 J3/2 2∈−π/6,π/6,(9) where J2=1 2rβ ijrβ jiand J3=1 3rβ ijrβ jkrβ ki. Note that βij corresponds to a shearing of the yield surface in the deviatoric direction, rather than a rotation away from the hydrostatic axis. This distinction is important, as an initially convex yield surface will remain convex for any degree of shearing. 2.4 Yield surface shape parameters The yield surface and direction of plastic flow can be tuned to simulate the behaviour of different fine-grain media using the two yield surface shape parameters parameters αand γ(in addition to the classical constants M,pcand ¯ρe). However, introducing anisotropy into the dissipation function results in the loss of uniqueness of the position of isochoric flow. Although this does not imply that the Critical State surface is no longer unique, it does remove one of the elegancies of the isotropic two-parameter Critical State model proposed by Collins and co-workers[2,3] and later developed further and implemented for finite-element analysis by Coombs and Crouch[7]. A constant Critical State surface in stress space requires the following: (i) the ratio of hydrostatic pressure to the size of the yield surface where ˙εp v= 0 is constant for any level of anisotropy, that is ¯pcs =(p/pc)cs is constant (where the subscript (·)cs denotes a quantity at the Critical State); and (ii) the stress ratio where ˙εp v= 0 is constant for any level of anisotropy, that is (q/p)cs = M¯ρ(θ). In order to recover this fundamental property of CSSM, first we equate the volumetric component of the direction plastic flow (7) to zero, giving ¯ B2=¯ A2¯p(η−β)β ¯p−γ/2, (10) where β=βijβij and η=q/p. Note that here for simplicity (but without loss of generality of the final result) the equations are presented in pversus qspace. The yield function (5), provides an alternative expression for ¯ B2 ¯ B2=¯ A2¯p(η−β)2 γ(2 −γ)(1 −¯p). (11) 4 1363 William M. Coombs Combining (10) and (11) provides an equation linking γwith the stress ratio at the Critical State, ηcs, and the current level of anisotropy in terms of a normalised parameter, ¯ β, γ2¯ β(¯pcs −1)+γ(1 −¯ β)/2+2¯ β(1 −¯pcs)−¯pcs(1 −¯ β)=0.(12) ¯ β=β/M ¯ρ(θ) is the ratio of the gradients of the current level of anisotropy and the position of the Critical State. γ∈[0,1] can be obtained by solving the quadratic (12) and selecting the positive root. The variation of γwith normalised anisotropy, ¯ β, for ¯pcs ∈[0.2,0.5] is shown in Figure 1 (i). In the limiting case of isotropy; ¯ β= 0 and γ= 2¯pcs, thereby recovering the isotropic formulation of[7]. Figure 1: Yield surface parameter variation with anisotropy for ¯pcs ∈[0.2,0.5] and M= 1: (i) γversus β,(ii) αversus βand (iii) αversus ¯pcs when β= 0. Rearranging (10) in terms of αallows the second shape parameter to be expressed in terms of the current level of anisotropy and the normalised pressure at the Critical State α=¯ AcsAβ−¯pcs γ/2−¯pcs ,where Aβ=¯pcs ¯ β1−¯ β ¯pcs −γ/2and ¯ Acs = (1 −γ)¯pcs +γ 2.(13) The variation of αwith normalised anisotropy for ¯pcs ∈[0.2,0.5] is shown in Figure 1 (ii). Specifying αand γthrough (13) and (12), respectively, maintains a constant Critical State stress ratio, ηcs. The limit of (13) when β→0 is not well defined due to the fact that when β=0a unique Critical State surface is obtained for any α∈[0,1]. The value of αwhen β≤tol was set to the value obtained when β= tol, where a tolerance (tol) of 1 ×10−5was found to give a stable result. The variation of the value of αwith ¯pcs ∈[0.2,0.5] when β≤tol is shown in Figure 1 (iii). Note that if ¯pcs =0.5 and β= 0, we recover the classical MCC ellipsoidal yield envelope, albeit with a non-circular deviatoric section. The variation in yield surface shape in normalised pversus qstress space is shown in Figure 2 for β= 0, 0.2, 0.4and0.6 with M= 1, ¯ρe=0.8 and ¯pcs =0.5and0.2. Note 5 1364 William M. Coombs Figure 2: Yield surface shape variation with anisotropy for β= 0, 0.2, 0.4and0.6 with M= 1, ¯ρe=0.8 and ¯pcs =0.5 and 0.2. that, due the relative evolutions of αand γ, increasing the level of anisotropy has the effect of reducing the deviatoric radius of the yield envelope. The yield surface of the proposed model has the following properties: (i) a unique Critical State surface for any degree of induced anisotropy; (ii) a constant ratio of deviatoric yield stress above and below the axis of anisotropy independent of βij; (iii) a narrowing of the deviatoric yield radius with increasing ¯ βdue to the reduction of γ, consistent with experimental findings on K0consolidated soils; and (iv) a requirement that the level of anisotropy must be restricted to β≤¯ρ(θ)M. The final point has important implications for the evolution of anisotropic shearing of the yield surface and is discussed in more detail in Section 2.6 below. The uniqueness of the position of the Critical State is demonstrated in Figure 3 (i), where the dilation angle, defined as φg= arctan˙εp v/˙εp γ, is plotted against the friction angle, arctan(q/p). The position of isochoric plastic flow remains at a friction angle coincident with the Critical State line (arctan(M)). On the compressive side of the Critical State line, for a given friction angle, increasing the level of anisotropy increases the level of plastic compaction. The level of anisotropy on the dilative side of the Critical State line (friction angles greater than π/4) has little influence on the dilation angle. 6 1365 William M. Coombs Figure 3: Direction of plastic flow: (i) dilation angle versus friction angle and (ii) associated yield surfaces for ¯pcs =0.5, M= 1 and ¯ρe=0.9. 2.5 Isotropic hardening/softening Following the same approach as[7,5,6], the rate of evolution of the size of the yield surface is given by ˙pc=˙εp vpc/(λ−κ).(14) This hardening law is equivalent to specifying a linear relationship between the bi-logarithm of the specific volume, v, and the pre-consolidation pressure, pc. 2.6 Anisotropic shearing Preserving the uniqueness of the Critical State by setting αand γas a function of ¯pcs limits the level of allowable anisotropy to β<¯ρ(θ)M. In order to maintain this condition, the evolution of anisotropy is specified through ˙ βij =||˙εp ij||Cβxβrij −βij,(15) 7 1366 William M. Coombs where ||(·)ij|| denotes the L2-norm of (·)ij. The parameter controlling the target level of anisotropy, xβ, is dependent on the current deviatoric stress ratio xβ=1−exp(−aβ)/aβwhere aβ=bβ||rij||/¯ρ(θ)M.(16) Cβand bβ≥1 are constants controlling the rate of evolution of anisotropy and the level of anisotropy under constant rij loading, respectively. For a constant stress ratio load path (that is, constant rij), increasing bβwill reduce the level of anisotropy generated in the model. The level of anisotropy at the Critical State is obtained by setting ||rij|| =¯ρ(θ)M βcs =¯ρ(θ)M1−exp(−bβ)/bβ.(17) 3 ALGORITHMIC IMPLEMENTATION The previous section has detailed the algorithmic development of the anisotropic constitutive model. However, these equations only provide a rate description of the model. In order for the model to be used in practical boundary value simulations (or even at a material point simulation level), these rate equations must be reformulated into an incremental relationship. Here a fully implicit backward Euler (bE) stress integration scheme is used. Note, that in this section we shift to matrix/vector notation to provide enhanced clarity for numerical implementation. First, implicit integration of (14) and (15) yields the following evolution laws ˜pc=pcn (1 −∆εp v/(λ−κ)) and {˜ βn+1}={˜ βnum} (˜ βdnm)={βn}+Cβxβ||{∆εp}||{r} 1+Cβ||{∆εp}|| . (18) The subscript ndenotes the previously converged solution associated with the last step (or the initial state at the start of the analysis). Here, we denote these evolution equations with a tilde to distinguish them from incremental updating through the bE method. Using the following fourteen dimensionless residuals {b}={εe}−{εe t}+∆γ{g,σ}1−˜pc/pc{β}−{˜ β}fT(19) and taking the derivative of the residuals respect to the following set of unknowns {x}= {εe}pc{β}∆γT, the (14×14) Hessian matrix is obtained as [A]=∂b ∂x=    [I]+∆γ[g,σσ ][De]∆γ{g,σpc}∆γ[g,σβ ]{g,σ} −(˜pc/pc),σT[De]−(˜pc/pc),pc−(˜pc/pc),βT−(˜pc/pc),∆γ −[˜ β,σ][De]−{˜ β,pc}[I]−[˜ β,β]−{˜ β,∆γ} {f,σ}T[De]f,pc{f,β}T0    . (20) The iterative increment in the unknowns is given by {δx}=−[A]−1{b}.(21) 8 1367 William M. Coombs The iterative increment of (·) is denoted by δ(·) using a lower-case delta to denote that this increment is the contribution to the unknowns for a given iteration. The total increment in the unknowns, {∆x}, is given by the summation of the iterative increments over the number of iterations required to converge within a specified tolerance. The iterative procedure starts from the following initial conditions {0εe}={εe t},0∆γ= 0, 0pc=pct and {0β}={βt}, the pre-superscript denotes the iteration number and (·)tdenotes a trial value. That is, ftis the value of the yield function at the trial stress state and pctand {βt} are the trial values of the size of the yield surface and the level of anisotropic shearing. For small strain analysis, these trial values are equal to the value of the parameter determined at the previously converged state. The Newton-Raphson iterative process continues until the residuals converge to within a specified tolerance on each of the four grouped residuals, typically 1×10−9. Throughout the stress return, all of the derivatives are evaluated at the current state. This requires the repeated evaluation of the derivatives at each iteration in addition to the inversion of the Hessian matrix (20). For the sake of brevity, the lengthy derivatives are omitted. 4 PHYSICAL COMPARISONS This section compares the ability of the proposed model to predict the experimental behaviour of Lower Cromer Till[9] (LCT). Figure 4 compares the proposed model1((i) to (iv)) with SANIclay ((v) to (viii))[8] under one-dimensional consolidation ((i) and (v)) and swelling ((ii) and (vi)), followed by undrained triaxial compression ((iii) and (vii)) or extension ((iv) and (viii)). The stress paths for SANIclay were obtained from[8], where the eight constants required for the model were calibrated on the same LCT experimental data as used in this paper. Unfortunately, the paper did no present the full one-dimensional loading and unloading behaviour of the model. The portions of the paths presented in that paper have bee reproduced in Figure 4 (v) and (vi). Dafalias et al.[8] allowed their model to start at a stress state in agreement with the experimental data for each of the individual triaxial simulations rather than simulating the material’s full stress history. The one-dimensional drained compaction followed by unloading of LCT (discrete points) is compared with the numerical prediction of the proposed model (solid line) in Figures 4 (i) and (ii). The proposed constitutive model stated from a hydrostatic stress and isotropic material state with a reference pressure and size of the yield surface equal to 75kPa. The model was then subjected to a one-dimensional compressive strain path in increments of ∆εz=1×10−4to a pressure of 233kPa followed by unloading to 62kPa. The model offers reasonable agreement with the experimental data under both (i) loading between Aand Band (ii) unloading, Bto D. Between Band Cthe model predicts elastic behaviour. The onset of yield occurs at Cand the model’s response has a notable change in gradient until arriving at D. 1The constants for the proposed single surface model are as follows: κ=0.007, G0= 2MPa, αe= 75, λ=0.044, M=0.96, ¯pcs =0.45, ¯ρe=0.73, Cβ= 80 and bβ=1.5. 9 1368