An undrained expansion solution of cylindrical cavity in SANICLAY for K0-consolidated clays

Haohua Chen, Pin-Qiang Mo

a China University of Mining and Technology, Xuzhou 221116, China

b The University of Arizona, Tucson, AZ 85716, USA

Keywords:Undrained cavity expansion SANICLAY K0-consolidation Anisotropy

A B S T R A C T

1. Introduction

The in situ stress state of natural soil deposits is usually inherently anisotropic,and it is common to define the at-rest coefficient as the ratio of horizontal(σ′h0)to vertical(σ′v0)effective stress at the geostatic condition,i.e.K0= σ′h0/σ′v0.Considering the sedimentary process without lateral deformation, the initial stress-induced anisotropy is particularly crucial to the mechanical behaviours of saturated clays (e.g. Nakase and Kamei, 1983; Ghaboussi and Momen, 1984; Mayne, 1985; Wang et al., 2008). Following the early developments on plasticity theory for geomaterials, many anisotropic constitutive models for soft clays have been proposed considering anisotropicK0-consolidation. Originated from the anisotropic extension of work by Dafalias (1986), the anisotropic critical state models were extensively developed with various types of yield surfaces and hardening rules (e.g. Whittle and Kavvadas,1994; Wheeler et al., 2003; Dafalias et al., 2006; Yin et al., 2010;Rezania et al., 2016). The SANICLAY (an abbreviation of ‘Simple ANIsotropic CLAY’)model,extended by Dafalias et al.(2006),is one of the milestones for simulation of both undrained and drained rate-independent behaviours of sensitive clays, using a nonassociated flow rule, a combined rotational and distortional hardening rule and the evolution laws for hardening variables. Extensions of the SANICLAY model have also been developed in the recent years for further considerations of destructuration, plastic deformations within yield surface, and cyclic behaviour, etc. (e.g.Taiebat et al., 2010; Jiang et al., 2012; Dafalias and Taiebat, 2013;Seidalinov and Taiebat,2014; Yang et al., 2019).

The cavity expansion theory has been greatly developed and widely applied in geotechnical engineering since 1940s(e.g.Bishop et al.,1945; Hill,1950; Chadwick,1959). Considering more sophisticated constitutive models for geomaterials, cavity expansion solutions have been progressively improved during the past half century(e.g.Vesic,1972;Carter et al.,1986;Yu and Houlsby,1991;Cao et al., 2001; Chen and Abousleiman, 2012; 2018; Mo and Yu,2017, 2018; Chen et al., 2020; Zhao et al., 2020). The derivations for stress and displacement fields around cavities are then employed for analyses of complex geotechnical problems,including in situ soil testing,foundations,underground excavations,wellbore instability. However, most of the early solutions are developed for cavities embedded in an infinite medium and isotropic stress field,leading to inaccurate results for solving boundary value problems in numerous applications. Correction factors based on numerical methods or empirical relationships were then introduced for the compromised use of analytical results (e.g. Naggar and Naggar,2012).

Over the last decade, many attempts have been made to tackle this issue providing semi-analytical solutions of cavity expansion or contraction under anisotropic initial stress conditions. Following the framework of Chen and Abousleiman (2012, 2013), semianalytical solutions of undrained and drained cylindrical cavity expansion were proposed by Li et al.(2016a,2017) to consider the effects of initial stress anisotropy by employing an associativeK0-consolidated anisotropic modified Cam-clay (AMCC) model(Sekiguchi and Ohta,1977).It was reported that the degree of initial anisotropy has a notable influence on the stress distributions around the cylindrical cavity, though the rotation of yield surface representing the stress-induced anisotropy was neglected for simplicity. Chen and Liu (2019) reported a rigorous undrained solution in the AMCC model (Dafalias, 1987), which considers both rotation and distortion of yield surfaces. The effects ofK0-consolidation and the stress-induced anisotropy were clearly shown on the distributions of effective stresses and excess pore pressure around the cylindrical cavity.This work was then extended for the drained solution (Liu and Chen, 2018). Similarly, an undrained solution in the S-CLAY1 model(Wheeler et al.,2003) was developed by Sivasithamparam and Castro (2018), considering the fabric anisotropy of soft soils and its evolution with plastic strains.It was later updated using the S-CLAY1S model (Karstunen et al., 2005),accounting for both fabric anisotropy and destructuration(Sivasithamparam and Castro, 2020).

Since the strength of clay is typically overestimated by the Mises criterion-based models, the spatially mobilized plane (SMP) criterion is suggested for the three-dimensional (3D) strength of the cohesive-frictional soil (Li et al., 2016b). Following the advanced cavity expansion solutions of Li et al. (2016b) and Chen and Abousleiman (2012), Chen et al. (2019) presented a generic stress transform approach for both undrained and drained cylindrical cavity expansion solutions in an SMP criterion revised AMCC model(Yao and Wang, 2014). The solution was then modified for the undrained contraction problems by Zhang et al. (2020). Additionally, the effects of rotational hardening, ignored in Chen et al.(2019), have been recently included in Yang et al. (2020) under drained loading conditions. Some recent solutions for cavity expansion theory are summarized and compared in Table 1.

Table 1 Summary of recent solutions to cavity expansion problem.

In this paper,an undrained expansion solution in the SANICLAY model,incorporating 3D strength of soil,a non-associated flow rule and a combined rotational and distortional hardening rule, is proposed for cylindrical cavities consideringK0-consolidation and stress-induced anisotropy. The solution with attached corresponding MATLAB source code (see Appendix A) is validated against that based on the recovered isotropic model by Chen and Abousleiman (2012), and the results ofK0-consolidated clays are presented to investigate the effects of overconsolidation ratio and model constants on the cavity expansion responses and the stress paths, along with the evolutions of anisotropic parameters. The highlight of this context also lies in the provided generic framework for the derivation of undrained cavity expansion in an anisotropic critical state model, and the solution using the SANICLAY model could serve as a useful benchmark for further developments and numerical calculations, along with the source code and a workedout example.

2. Definition of cavity expansion problem

The problem of this study concerns with the evolutions of stress and displacement fields during the undrained expansion of a vertical and cylindrical cavity inK0-consolidated clays.Assuming the planestrain condition along the vertical direction,the in-plane schematic of the cylindrical cavity expansion problem is illustrated in Fig.1.The initial stress condition of the infinite soil mass is considered asK0-consolidated with in situ horizontal and vertical effective stresses σ′h0and σ′v0.The case σ′r0≠σ′θ0will not be considered in the present study,in which σ′r0and σ′θ0are the initial radial and circumferential stresses, respectively. The circular cavity expands from its initial cavity radiusa0to its current cavity radiusaunder an internal cavity pressure σa. During the process of expansion, soil mass near the cavity wall is squeezed outwards under the increase of internal cavity pressure,and a plastic region with a radius ofrpis developed around the cavity wall,which is embedded within the outer elastic region. In the framework of critical state soil mechanics, further expansion may create a critical state region adjoined to the cavity wall.The quasi-static stress state of an arbitrary soil element around the cylindrical cavity at any time is governed by the equilibrium equation in the radial direction alone:

Fig.1. Schematic illustration of cylindrical cavity expansion problem.

where σ′rand σ′θare the effective radial and circumferential stresses,respectively;uis the pore water pressure;ris the radius of a soil element;d()is the Eulerian derivative for every material particle at a specific moment.

To accommodate large strains in the plastic and critical state regions during the cavity expansion process,the logarithmic strains are assumed as (Chen and Abousleiman, 2012; Chen et al.,2021):

where εrand εθ are the radial and circumferential strains, respectively;r0is the initial position of a soil element before expansion.The undrained expansion condition used to formulate the problem later indicates that the specific volumevremains constant during expansion, and the volumetric strain εvof soil around the cavity vanishes everywhere,which gives

where εzis the vertical strain, which is equal to zero under the plane-strain condition.

3. SANICLAY model

The fundamental formulation of the SANICLAY model developed by Dafalias et al.(2006)will be adopted in the present solution.The tensorial form of the plastic potential in the general stress space is given as

wheresis the deviatoric stress tensor, defined ass=σ -p′I(σ is the effective stress tensor,Iis the identity tensor, andp′is the effective mean stress withp′=(1/3)trσ);α is the non-dimensional anisotropic variable tensor, which serves for the rotational hardening of the plastic potential surface; the symbol ‘:’ implies the trace of the product of two tensors;Mis the critical stress ratio;p′αis the value ofp′atq=p′α, whereqis the deviatoric stress and α = [(3/2)α:α]1/2. It should be noted thatp′αis determined by substituting the tensorssand α into Eq.(4),and no hardening law is adopted to updatedp′α. In the general stress space, the 3D soil strength is considered by means of the Lode’s angle θL, and the Lode’s angle dependent formulation ofMis defined as follows:

whereMcandMeare the critical stress ratios in compression and extension,respectively.According to Jiang and Pietruszczak(1988),the shape of Eq. (5) in the π plane might be concave whenMe/Mc< 7/9, which is a shortcoming of Eq. (5) compared to other methods for introducing 3D strength such as the stress transformation method (Matsuoka and Sun, 2006). The potential function is a curved triangle as shown in Fig.2,and bothMcandMeare subtly considered in SANICLAY.Hence,the present solution can be seen as a more general solution to the cavity expansion problem than the exiting circular critical state line-based solutions, such as Chen and Abousleiman(2012),Sivasithamparam and Castro(2018).

The non-associated flow rule was employed in the SANICLAY model, which means the yield surface does not overlap the plastic potential surface.The tensorial form of the yield surface function in the general stress space is defined as

Fig.2. Yield surface and plastic potential of the SANICLAY model in(a)p′-q space and(b) π plane.

where β is the rotational hardening variable tensor of yield surface;p′cis the value ofp′at η = β0, where η is the stress ratio defined asq/p′and β0= [(3/2)β:β]1/2, which represents the isotropic hardening variable;Nis a soil constant that serves as the bound for β and is assumed independent on the Lode’s angle.Again, the proposed solution with a non-associated flow rule can be regarded as a more general case than the existing associated constitutive model-based solutions, which is another predominant feature of this study.

Both yield surface and plastic potential surface are schematically illustrated in Fig. 2 in thep′-qspace and the deviatoric plane, respectively. The critical stress state at failure is unnecessarily located at the peak of the yield surface in thep′-qspace(Fig. 2a), contributing to the undrained softening afterK0-consolidation. The shape of yield function in the π plane is a circle, whereas the plastic potential surface is anisotropic owing to dependence on the Lode’s angle. To describe the evolution of the internal variables serving as three hardening laws, the isotropic hardening parameterp′cand the anisotropic rotational hardening parameter tensors α (for plastic potential) and β (for yield surface) are defined in terms of their rate forms as

whereLis the loading index,which will be specified later from the plastic consistency condition;〈〉 is the Macaulay bracket;e0is the initial void ratio, v0is the initial specific volume; λ and κ are the slopes of normal compression and swelling lines, respectively, in thee-lnp′space(eis the void ratio);Cis a model constant,which represents the rate of evolution of anisotropy; the subscript ‘b’ is used for the bounding ‘image’ of the variable;xis another model parameter to define the attractor tensorr/xfor α, indicating the saturation limit of anisotropy.

4. Undrained cavity expansion solution

4.1. Constitutive relationship in matrix form

In the classical plasticity theory,the incremental total strains in the plastic phase can be divided into elastic and plastic components, with superscripts ‘e’ and ‘p‘, respectively. In terms of the cylindrical coordinates in this problem, the rates of the elastic strains can be represented by the effective stresses in the matrix form, following the small-strain assumption and the Hooke’s law for isotropic mass:

where μ is the Poisson’s ratio;Eis the elastic modulus, which depends on the current effective mean stressp′,specific volume v and basic soil parameters as

Based on the plastic flow rule(see Eq.(4)),the rate of the plastic strain tensor εpcan be written as

where

For determination of the loading indexL,the plastic consistency condition needs to be fulfilled as follows:

Substituting Eq.(7)into Eq.(12)yields the expression of loading indexL:

where

The rates of the three plastic strain components can thus be rewritten in the matrix form based on Eq. (10):

where the matrix elements are given as

Combining the stress-strain relationships in Eqs. (8) and (15),the rates of stress components can be inversely expressed in terms of the incremental total strain components as follows:

where the matrix elements are given as

4.2. Solution in the elastic region

At aK0-consolidation condition before expansion, the initial stress relation gives σ′h0=K0σ′v0. Considering the rotational hardening of the yield surface and plastic potential, the initial values of the dimensionless anisotropic variable tensors α0and β0are related to the coefficient of lateral earth pressure at restK0in terms of the cylindrical coordinates,as follows(Dafalias et al.,2006):

The elastic region is located outside of the elastoplastic boundary withr>rpafter cavity expansion. The elastic solution can then be expressed following the small-strain elasticity (Yu,2020), as

where σ′rpdenotes the effective radial stress at the elastoplastic boundary,which will be determined in the following section;Uris the radial displacement;G0is the shear modulus in the elastic region, takingG=E/[2(1+μ)] with Eq. (9). Note the excess pore pressure in the elastic region remains to be zero, owing to the unchanged mean stresses.

4.3. Solution in the plastic region

4.3.1. Conditions at elastic-plastic boundary

Soil states at the elastoplastic boundary are taken as the initial conditions for the solution within the plastic region.Combining the yield surface function,initial anisotropic tensor,and the stress state at the elastic-plastic boundary (Eqs. (6), (19b), and (20a)-(c)), the value of σ′rpcan be obtained as

wherep′0is the in situ value of the effective mean stress(i.e.p′0=(1 + 2K0)σ′v0/3);OCRdenotes the overconsolidation ratio and is defined asp′c0/p′0. Note that the value ofOCRin this paper represents the isotropic overconsolidation ratio, in terms of the mean effective stress rather than the conventional vertical effective stress, following Mo and Yu (2017).

The ratio of the current position of one soil particle to its original position is introduced herein to indicate the relative deformation as a time variable,i.e.ρ =r/r0.From the displacement distribution in the elastic region(see Eq.(20d)),the value of ρ at the elastic-plastic boundary is derived as

The range of ρ within the plastic zone is given as

Considering the undrained condition with Eqs. (2) and (3), the current radial coordinate can be determined from ρ:

4.3.2. Governing differential equations

In conjunction with the introduced definition of ρ and the undrained condition for a cylindrical scenario,the incremental strains can be re-written as

where D()is the material derivative along the particle motion path using the Lagrangian description.

Considering the undrained condition and the definition of ρ,the Eulerian derivative for every material particle at a specific moment can also be related to ρ as

It should be noted that Eqs.(25)and(26)can be easily extended to their drained versions by taking the specific volume as an extra unknown. Then the proposed solving approach can be easily extended to the drained case, which shows the advantages of the utilized time variable ρ-based solving technique over the existing solutions.

By substituting the strain definitions of Eq. (25) into the constitutive matrix of Eq.(17),the components can be expressed as

It should be noted that bothp′candp′αare expressed by stresses,α and β according to Eqs. (4) and (6) during the updating of the constitutive matrix. Therefore,p′α andp′care not involved in the constitutive matrix and not updated by integration.

Combining Eq.(1)with Eq.(26),the excess pore water pressure can be obtained as

It should be noted that Eqs. (27a-c) and (28) are first-order differential equations with respect to ρ. Substituting Eq. (25) back into Eqs. (7) and (13a), one can also express the rates of isotropic hardening parameterp′cand anisotropic parameters α and β as functions regarding ρ, which are given in Appendix A. Therefore,the governing equations for the cavity expansion problem can now be simplified to a system of first-order differential equations,which can be solved by the Runge-Kutta method (Atkinson, 1989) via a computational software as shown in Fig. 3. The corresponding MATLAB source code can be found in Appendix A to facilitate the use of the presented solution, as well as a worked-out example.Since the cavity in the infinite medium expands in a self-similar way, all the soil particles experience the same stress path.Considering the relative deformation along with the position to the cavity wall,the stress and displacement distributions in both elastic and plastic regions are obtained,as well as their evolutions during the cavity expansion process.

5. Results and discussion

5.1. Validation against MCC solution

The proposed solution is firstly validated against the MCC solution proposed by Chen and Abousleiman (2012), excluding the anisotropy and non-associativity by settingN=Mc=Me,x= 1 andC= 0.The soil parameters are set equivalent to those for the normally consolidated Boston Blue clay withK0= 0:625, after Chen and Abousleiman (2012), as detailed in Table 2.

The normalized distributions of stress components after a certain expansion witha/a0= 2 are compared in Fig. 4. The effective stresses(σ′r,σ′θ,σ′z)and excess pore pressure(Δu)based on the recovered SANICLAY model are identical to those of the MCC solution, for both critical state and plastic regions, validating the developed formulation in the previous section. Additionally,comparing to the MCC solution, the advantages of the current solution lie on the rotation and distortion of both yield surface and plastic potential with respective evolution law consideringK0-consolidation and stress-induced anisotropy, as follows.

5.2. Influence of overconsolidation ratio

A parametric study is then conducted to present the results of cylindrical cavity expansion using the foregoing solution, examining the effects of overconsolidation ratio.All the SANICLAY model parameters are set as those from Dafalias et al.(2006)for the Lower Cromer Till(LCT)clay,and four tests with differentOCRin the range of 1-10 are considered with identical initial effective mean stress(p′0= 120 kPa). The corresponding initial state parameters are summarized with details in Table 3.

Table 2 Parameters of special case for comparison of present solution and MCC solution.

Table 3 Stress conditions and material properties of the LCT clay involved in parametric analysis. (Extended and revised from Chen and Abousleiman, 2012).

The variations of the normalized cavity pressure (σa/p′0) and excess pore pressure at the cavity wall (Δu/p′0) during expansion froma/a0= 1 toa/a0= 10 are shown in Fig. 5, for tests of soil with various values ofOCR.The cavity pressure increases rapidly at the initial expansion stage witha/a0<3,which tends to approach an asymptotic limit pressure at a larger expansion period. The normalized cavity pressure appears to increase withOCR, and similar trends were also reported by Sivasithamparam and Castro(2018). The excess pore pressure at the cavity wall, as depicted in Fig. 5b, shows generally a gradual increase ata/a0<4, and also approaches the limiting value at notable expansion.The soil with a larger value ofOCRis found to have a higher limiting value of normalized excess pore pressure. However, slight negative excess pore pressure appears during the very early expansion witha/a0<1:1 for heavily overconsolidated soil withOCR= 10, which is consistent with those in Chen and Abousleiman(2012)and Li et al.(2016b).

Fig. 3. Solving procedures for soil particle in plastic zone.

Fig. 4. Distributions of stress components of special case for comparison of present solution and MCC solution.

The distributions of effective stresses and excess pore pressure at an expansion instant ofa/a0= 2 are provided in Fig. 6, for all cases with various values ofOCRin the range 1-10. All stress components are normalized by the initial effective mean stressp′0,and the radial distance of soil element to the cavity centerris normalized by the current cavity radiusa. For normally consolidated clay (Fig. 6a), any degree of expansion causes the plastic yielding, leading to an infinite plastic region. For the isotropically consolidated clay withK0= 1 in Fig. 6b, the effective stresses converge to an identical constant at far field of the elastic region for undisturbed soil element after expansion. Stress changes at the elastoplastic boundary are shown unsmooth, especially for the circumferential stress. The radial stress in the elastic region decreases with radial distance, whereas the mean stress, vertical stress and pore pressure remain constants during the elastic stage,as reflected in Eq. (20). The effective stresses in the vicinity of the cavity wall are uniform within the critical state regions, while the excess pore pressure decreases almost linearly with the logarithmic scale of the radial distance and the size of the critical state region decreases with the value ofOCR.The size of the plastic region is also found to decrease with the overconsolidation ratio, and negative excess pore pressure is shown in the distribution curve for heavily overconsolidated clay (Fig. 6d). Note that the effective vertical stress is not always the intermediate principal stress in the critical state region, for example as shown in Fig. 6d, which was also reported by Chen and Liu (2019) but is against that of Chen and Abousleiman (2012), where the ultimate vertical stress equalled to the average of radial and circumferential stresses.

Fig. 5. Variations of (a) normalized cavity pressure and (b) normalized excess pore pressure at cavity wall during expansion in soil with different values of OCR.

Fig. 6. Distributions of normalized stresses and excess pore pressure after undrained expansion of a/a0 = 2.(a)OCR = 1,K0 = 0:625;(b)OCR = 3,K0 = 1;(c)OCR = 5,K0 = 1:3; (d) OCR = 10, K0 = 2.

Correspondingly, Fig. 7 presents the distributions of anisotropic parameters(α and β) after undrained expansion ofa/a0=2 for all tests with various values ofOCR. All parameters in both critical state and elastic regions remain constant, while the monotonical transitions appear in the plastic regions, indicating the rotation and distortion of yield surface and plastic potential with plastic yielding. It should be noted that the magnitudes of

radial and circumferential parameters are identical at the elastic stage (i.e. αr0= αθ0and βr0= βθ0), due to the initial vertical cross anisotropy. They start to deviate from each other at the elastoplastic boundary, and the amounts of changes of parameters represent the degree of induced anisotropy. Following the definitions of Eq. (19), the anisotropic parameters converge to zeros in the elastic region for isotropically consolidated clay(Fig. 7b). The radial parameters (αrand βr) tend to increase with expansion while the circumferential parameters show opposite trends for all tests with theOCRin range of 1-10. It is interesting to find that the vertical parameters (αzand βz) decrease with expansion for normally consolidated and lightly overconsolidated soils, whereas they increase slightly with expansion for heavily overconsolidated soil (Fig. 7d). However, the restrictive relations of αr+αθ+αz= 0 and βr+βθ+βz= 0 hold true during the whole stages of cavity expansion, as also reported by Chen and Liu (2019). This is because the following equations maintain zero throughout the expansion process:

The effective stress paths (ESP) of cavity expansion in both normally consolidated and heavily overconsolidated soils are depicted in thep′-qplane with normalization of the initial effective mean stress (Fig. 8). The initial stress states, represented by the notation‘O’in the figure,are located at theK0-line(q/p′=3|1 -K0|/(1 +2K0))forK0-consolidated soil,and both initial yield and plastic potential surfaces plotted in dash lines intersect at point‘O‘. As for critical-state theory, the hardening/softening behaviour of soil is associated with stress ratio η(q/p′)rather than shear stressqaccording to Wood (1990). The expansion for normally consolidated soil (Fig. 8a) causes immediately plastic yielding with rotation and distortion of yield and plastic potential surfaces. The ESP describes the evolutions of effective stresses during undrained expansion,and clearly reflects the hardening behaviour(increase of stress ratio η) with the loading history of cavity expansion. The stress state approaches to the final point (’F’ in the figures) at the critical state line (CSL). It should be noted that the ESP does not intersect with the final potential surface because the extension of the potential surface fromp′axis is much smaller during loading than its final state (Fig. 8b). For heavily overconsolidated soil in Fig. 8b, the vertical elastic trajectory is firstly shown to reach its initial yield condition at the‘dry’side of the CSL(q/p′>M).Plastic softening behaviour (decrease of stress ratio η) is then observed during the plastic stage, and the stress path eventually terminates at a critical state as well.

In terms of the anisotropic stress condition in conjunction with the distortion and rotation of yield and plastic potential surfaces,the ESPs are depicted in the deviatoric (π) plane (Fig. 9), for both normally consolidated and heavily overconsolidated soils. The initial‘O’points along the σ′zaxis indicate theK0stress anisotropy,and the axisymmetric location of the initial plastic potential implies the initial cross-anisotropy and depends on the magnitude ofOCRvalue.The shape of yield surface is a circle in the π plane,while the plastic potential function gives a smoothed-triangular surface owing to the non-associated flow rule and the 3D definition of critical stress ratio in Eq.(5).The shape of the potential function is concave as the ratio ofMetoMcin this study is only 0.729. This value is less than the convexity criterion value 7/9 for Eq. (5) according to the research on the convexness of Eq. (5) by Jiang and Pietruszczak (1988). It should be noted that the initial yield surface returns to a single point ‘O’ for normally consolidated soil(Fig. 9a). The evolutions of yield and plastic potential surfaces are clearly shown, along with the ESPs during undrained expansion.Note that the horizontal elastic trajectory in Fig. 9b represents the unchanged vertical stress during the elastic expansion,and the final critical states are unnecessarily located at a specific Lode’s angle(θL= π/2), as reported by Sivasithamparam and Castro (2018,2020). It should be noted that the elastic volumetric strain rate plasticevis null at the critical state.Then considering the undrained condition, the volumetric strain ratepvshould also vanish (i.e.tr (?g/?σ) = 0) at the critical state, thus no unique relation between α andrcan be found as indicated by Eq.(7b).Following Chen et al. (2019), combing Eq. (27a-c) gives

Due to the existence of ?g/3?θL(?θL/?σ′r+?θL/?σ′θ-2?θL/?σ′z)and the lack of determined relation between α andr,Eq.(30)does not guarantee σ′z= (σ′r+σ′θ)/2 and αz= (αr+αθ)/2 at the critical state and this is the reason that the ESPs would not always terminate at θL= π/2.

5.3. Effects of model constants

The test of normally consolidated clay in Table 3 (OCR= 1) is taken as a reference in this section, for investigating effects of the three new model constants of SANICLAY(i.e.N,xandC),compared with the MCC model. The parameter ‘N’ was used to define the shape of the yield surface in Eq. (6) that introduces the nonassociativity compared with the plastic potential defined by Eq.(4) with anisotropic value ofM. A restrained condition was also noted by Dafalias et al.(2006)withMe≤N

Fig.7. Distributions of anisotropy parameters after undrained expansion of a/a0 = 2.(a)OCR = 1,K0 = 0:625;(b)OCR = 3,K0 = 1;(c)OCR = 5,K0 = 1:3;(d)OCR =10, K0 = 2.

Fig. 8. Effective stress paths in p′ -q plane during undrained expansion of a/a0 = 2.(a) OCR = 1, K0 = 0:625; (b) OCR = 10, K0 = 2.

The parameter‘x’was introduced to define the‘attractor’tensor,for α, that accounts for the anisotropy in the plastic potential. In Fig.11, it shows the influences ofxon the stress distributions and the stress paths of cavity expansion.Whenxvaries between 1 and 2, the influence on the distributions of effective radial and circumferential stresses and the excess pore pressure is relatively limited, while the decrease of effective vertical stress with the magnitude ofxis more distinct in the plastic and critical state regions.The stress paths in thep′-qplane(Fig.11b)show that both effective mean and deviatoric stresses at failure decrease marginally withx, whereas the stress ratio at the critical state appears to increase slightly withx.It is worth noting that the CSL differs from each other,despite the identical stress ratios at the critical state in compression and extension (i.e.Mc= 1:18,Me= 0:86). It is attributed to the different rotational hardening of the plastic potential with changes ofx,and the critical stress ratioMdepends on Lode’s angle that can be observed from the π plane (Fig.11c). The ultimate Lode’s angle decreases nonlinearly withx,and the size of the critical state surface at failure in the π plane reduces withx,indicating less plastic hardening.

The parameter ‘C’ was adopted in Eqs. (7b) and (7d) for determination of the rates of evolutions for the anisotropy variables(i.e.and),and Dafalias et al.(2006)suggested that the magnitude ofCis usually between 3 and 20 for various clays.The effects ofCare presented in Fig. 12, with the variation ofCbetween 11 and 21.Higher effective stresses are noticed for tests with largerCvalue,and the excess pore pressure in the critical state region increases slightly withC. Fig. 12b shows that the stress paths in thep′-qplane for cavity expansion rise withCvalue,and the stress ratios at failure are quite close with a tiny growth againstC. It is also illustrated in Fig. 12c with the stress paths in the π plane, providing larger critical state surfaces and lower Lode’s angles at failure for tests with higherCvalue.

Fig.9. Effective stress paths in π plane during undrained expansion of a/a0 = 2.(a)OCR = 1, K0 = 0:625; (b) OCR = 10, K0 = 2.

Fig.10. Effects of N on(a)stress distributions after cavity expansion;(b)stress paths in p′-q plane; (c) stress paths and critical state surfaces in π plane.

It is noticed that the model constants have significant effects on the stress distributions and paths during the undrained cavity expansion,and the uniqueness of the critical state is affected by the anisotropy and the evolution laws. The developed solution can provide more possibilities for accurate predictions and backcalculations, with considerations ofK0-consolidation and the stress-induced anisotropy. The proposed solution in SANICLAY is rigorously formulated,and the governing differential equations are numerically solved for the stress and strain fields around the expanded cavities. Further considering the rate-dependent behaviour and destructuration upon loading,the presented framework is thus suitable for extensions using other SANICLAY based soil models (e.g. Taiebat et al., 2010; Seidalinov and Taiebat, 2014;Rezania et al.,2016).The current solution without introducing any assumptions, along with the provided source code, can also potentially serve as a benchmark for the extended analytical solutions and for the simplified numerical simulations of cavity expansion problems.

Fig.11. Effects of x on(a)stress distributions after cavity expansion;(b)stress paths in p′ -q plane; (c) stress paths and critical state surfaces in π plane.

Fig.12. Effects of C on(a)stress distributions after cavity expansion;(b)stress paths in p′-q plane; (c) stress paths and critical state surfaces in π plane.

6. Conclusions

A rigorous semi-analytical solution of undrained cylindrical cavity expansion in the SANICLAY model is proposed in this paper,considering theK0-consolidation and the stress-induced anisotropy. The corresponding MATLAB source code is attached to facilitate the use of the solution and to serve as a benchmark for validations and extensions. Comparing to the MCC model, three model constants were introduced to include the 3D strength,nonassociativity, rotational and distortional hardening, and the evolution laws in the multiaxial stress space.By utilizing a time variable ρ, both the position derivative and particle derivative can be converted to the same framework,thus the time variable-based solving technique can be used for both undrained and drained cavity problems. The exact derivation yields to the governing first-order differential equations, and the stress and strain distributions around the cavity after a certain expansion can be obtained effectively for all elastic, plastic and critical state regions, by adopting the Runge-Kutta method.The developed solution is firstly validated against the well accepted MCC solution, by setting constants for model recovering.Four tests in Lower Cromer Till clay with various overconsolidation ratio or alternatively the at-rest earth pressure coefficientK0are then conducted to examine the influences on the cavity pressure curves, stress distributions, evolutions of anisotropic parameters, and stress paths. The 3D evolutions of yield surface and plastic potential during cavity expansion inK0-consolidated clays are accurately captured by the solution in a non-associated anisotropic model without introducing assumptions and complexities.The effects of new model constants are also investigated, indicating the abilities with various stress paths,critical state surfaces and Lode’s angles at failure. The proposed procedure is also capable for further extensions,and the solution is potentially useful for verifying numerical results and for backanalyses of geotechnical problems, including pressuremeter tests and pile installation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We acknowledge the funding support from National Natural Science Foundation of China (Grant Nos. 51908546 and 52178374)and China Postdoctoral Science Foundation (Grant No.2020T130699).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.10.016.