Undrained expansion of a cylindrical cavity in clays with fabric anisotropy: theoretical solution

This paper presents a novel, exact, semi-analytical solution for the quasi-static undrained expansion of a cylindrical cavity in soft soils with fabric anisotropy. This is the first theoretical solution of the undrained expansion of a cylindrical cavity under plane strain conditions for soft soils with anisotropic behaviour of plastic nature. The solution is rigorously developed in detail, introducing a new stress invariant to deal with the soil fabric. The semi-analytical solution requires numerical evaluation of a system of six first-order ordinary differential equations. The results agree with finite element analyses and show the influence of anisotropic plastic behaviour. The effective stresses at critical state are constant, and they may be analytically related to the undrained shear strength. The initial vertical cross-anisotropy caused by soil deposition changes towards a radial cross-anisotropy after cavity expansion. The analysis of the stress paths shows that proper modelling of anisotropic plastic behaviour involves modelling not only the initial fabric anisotropy but also its evolution with plastic straining.

Abstract This paper presents a novel, exact, semi-analytical solution for the quasi-static undrained expansion of a cylindrical cavity in soft soils with fabric anisotropy. This is the first theoretical solution of the undrained expansion of a cylindrical cavity under plane strain conditions for soft soils with anisotropic behaviour of plastic nature. The solution is rigorously developed in detail, introducing a new stress invariant to deal with the soil fabric. The semianalytical solution requires numerical evaluation of a system of six first-order ordinary differential equations. The results agree with finite element analyses and show the influence of anisotropic plastic behaviour. The effective stresses at critical state are constant, and they may be analytically related to the undrained shear strength. The initial vertical cross-anisotropy caused by soil deposition changes towards a radial cross-anisotropy after cavity expansion. The analysis of the stress paths shows that proper modelling of anisotropic plastic behaviour involves modelling not only the initial fabric anisotropy but also its evolution with plastic straining.
Stress ratio: g = q/p 0 or g ¼ r d =p 0 (tensor) h Lode's angle: Slope of swelling line from t À ln p 0 space k Slope of post-yield compression line from t À ln p 0 space m Poisson's ratio r,r 0 Total and effective stresses r a Internal cavity pressure r p Total radial stress at the elastic/plastic boundary r 0 d deviatoric stress tensor / Friction angle x, x d Absolute and relative effectiveness of rotational hardening

Introduction
There is a wide variety of practical problems that may be modelled as the expansion (or contraction) of a spherical or cylindrical cavity in a solid mass. The mathematical solutions to those problems are usually categorized within the cavity expansion theory. The first solid mechanics applications of the cavity expansion theory were for metal indentation problems (e.g. [1,11]). For geomaterials, the application came later [10] but has been widely investigated because of its utility in many practical situations. Some examples comprise the interpretation of in situ tests like pressuremeter (e.g. [19,20]) or penetrometer tests (e.g. [6,29]) and the study of the installation disturbances caused by foundation elements like driven piles (e.g. [22]) or stone columns (e.g. [5]). It is also useful for wellbore instability and deep tunnels because, although in these cases the cavity is contracted instead of being expanded, the problem is similar mathematically (e.g. [23,32,35]). Here, the analysis limits to the quasi-static expansion of a cylindrical cavity in plane strain conditions because the solid mass is assumed as infinite. In clays, the cavity is usually expanded in a short period of time, and then, no drainage is allowed (undrained conditions).
Cavity expansion solutions for soils under undrained conditions use mainly isotropic elastic-perfectly plastic models, using, for example, Tresca criterion (e.g. [10,12,24]), or isotropic hardening constitutive models, such as the modified Cam clay (MCC) model (e.g. [4,7,8,25]). There are some solutions for anisotropic materials within elasticity, e.g. cross-anisotropic material [14,15], but, to the authors' knowledge, there were not any theoretical solutions for anisotropic behaviour of plastic nature. Only very recently, Li et al. [17] have published an analytical solution that accounts for an initial stress-induced anisotropy using a rotated yield surface. The rotated yield surface is fixed, and consequently, stresses near the cavity are not usually at a zero Lode's angle [21], which corresponds to plane strain conditions r 0 Á . This leads to unrealistic results in materials whose anisotropy evolves with plastic strains (e.g. clays).
This paper presents a novel, exact, semi-analytical solution of the undrained expansion of a cylindrical cavity in natural clays, which exhibit fabric anisotropy. The presented solution uses a constitutive model that considers anisotropy of plastic nature that evolves with plastic strains, both volumetric and deviatoric strains. The solution goes a step further than Li et al. [17] by introducing the evolution of anisotropy with plastic strains. The solution is semi-analytical as it requires the numerical evaluation of a system of six first-order ordinary differential equations. The results show the influence of the anisotropic plastic behaviour. Soil initial fabric is generated by soil deposition and consolidation due to the vertical compression caused by the gravity acceleration (vertical axis), whereas during cylindrical cavity expansion, the soil is compressed radially and its fabric anisotropy changes accordingly (radial axis).
The paper reviews existing analytical solutions of the cylindrical cavity expansion problem (Sect. 2) and details the assumptions and constitutive model used to develop the present solution (Sect. 3). Next, the full analytical development and its solution procedure are presented (Sect. 4).
Boston blue clay, whose properties are depicted in Sect. 5, is the soft clay used for validation against finite element analyses, which are presented in Sect. 6. Finally, the results using the semi-analytical solution are discussed (Sect. 7) and some conclusions are derived.

Undrained cavity expansion theory in soft soils
By definition, the volumetric strain is null under undrained conditions. If additionally, the direction of the displacement vector is known at each point, only its value being unknown, then, the strain field is independent of the stresses and of the constitutive model and may be obtained using the boundary conditions. This happens, for example, for the present case of the expansion of a cylindrical cavity in plane strain conditions, in which the displacement vector at any point is horizontal and passes through the axis of symmetry. Initial stresses and material properties must also satisfy those symmetry conditions. So, the problem reduces to a one-dimensional boundary value problem. The strain field is first obtained from the incompressibility condition, and then, the constitutive law is used to derive the effective stresses. Finally, equilibrium conditions may be imposed to get the internal pressure of the cavity.
A quite comprehensive review of solutions for different constitutive models may be found in Yu [34]. For critical state models, Collins and Yu [8] developed a general approximate large-strain solution for both original and modified Cam clay models. For the original Cam clay model, they found a closed form solution, while for the modified Cam clay model numerical integration is needed. As pointed out by Silvestri and Abu-Samra [25], Collins and Yu [8] solution is approximate because it uses a simplified definition of the deviatoric stress, q, as will be explained in the next section. Chen and Abousleiman [7] were the first ones to obtain an exact analytical solution using the rigorous definition of the deviatoric stress, q, and a shear modulus, G, that varies with the mean pressure, p 0 . Vrakas [30] developed a general exact solution for different Cam clay models and presented a critical evaluation of the various simplifying expressions used for stress invariants.
3 Definition of the problem 3

.1 Geometry and assumptions
The quasi-static expansion of a cylindrical cavity of initial radius a 0 is studied. The axis of the cylindrical cavity is assumed as the vertical axis, and the initial stress state is homogeneous and consists of a horizontal effective stress and a vertical effective stress r 0 H ; r 0 V À Á . The initial stress state may be also formulated in terms of total stresses considering the initial pore water pressure (u 0 , r H , r V ). The initial horizontal stress on the cavity is also r H , and it increases up to r a , upon expanding the cavity to a final radius a (Fig. 1). The soil constitutive model will be detailed below but, despite being anisotropic, it will have initial cross-anisotropy (transversely isotropic material) with the main axis being the vertical one. So, the problem has axial symmetry and, due to the initial uniform stress state and the infinite extent of the soil and the cavity, plane strain conditions hold. In this way, the strain field is easily obtained as for isotropic incompressible materials. Cylindrical coordinates (r, h, z) are used throughout the paper because they are principal directions for this problem. The equilibrium equation in the radial direction for cylindrical coordinates that are principal directions may be written as or using effective stresses as Initially, the clay may be normally consolidated or overconsolidated, i.e. in a plastic state or inside the elastic region. For the latter case, the entire soil may remain elastic after expanding the cavity or there may be a plastic annulus around the cavity. Its external radius will be denoted as r p and defines the current elastic/plastic boundary ( Fig. 1). Obviously, if the cavity is further expanded, the value of r p gradually increases.
The cylindrical cavity expansion problem is an appropriate example of the importance of accounting for large displacements [31]. For large displacements, the internal cavity pressure approaches an asymptotic limit value, while for small displacements, the internal pressure continuously increases because it is considered to be applied on a cylinder with a smaller diameter (a 0 ) than the real one (a). In this paper, the material incompressibility makes it easier to account for large displacements. So, the current position of an arbitrary point, r x , is directly related to the initial position of the point, r x0 , and the initial and current radii of the cavity, a 0 and a, respectively or, using dimensionless terms Large-strain deformation is considered in the plastic region using natural (or logarithmic) strains, but smallstrain deformation is used in the elastic region. This simplifying assumption does not affect the results because in the elastic region, the strains are much smaller than in the plastic annulus. The solution by Vrakas [30] considers large-strain formulation also in the elastic zone, but the differences are negligible.

Constitutive model
Natural soft clays exhibit a significant degree of anisotropy in their fabric, which initially is derived from the shape of the clay platelets, deposition process and one-dimensional consolidation. Fabric anisotropy of natural clays is modified due to subsequent irrecoverable straining (e.g. [33]). Reorientation of particles and changes in particle contacts, i.e. changes in fabric anisotropy, cause changes in the mechanical response of the soil. Additionally, the more the subsequent loading path differs from the loading path that has created the current anisotropy, as happens for the cylindrical cavity expansion problem, the more anisotropy changes [13]. Therefore, it is important to use a constitutive model able to reproduce not only the initial anisotropic behaviour of the soft clay but also its evolution with plastic straining. In this study, the S-CLAY1 model [33] is used. S-CLAY1 is a Cam clay-type model with an inclined yield surface to model inherent anisotropy and a rotational component of hardening to model the development or erasure of fabric anisotropy during plastic straining. For the simplified stress space of triaxial compression (r 2 = r 3 ) and for an initial cross-anisotropy fabric with the main axis being the vertical one (e.g. a vertically cut sample), the yield curve is a sheared ellipse [9] f y ¼ q À ap where q is a deviatoric stress (q = r 1 -r 3 ), p 0 is the mean effective stress, M is the critical state value of the stress ratio (where g = q/p 0 ), and p 0 m and a define the size and inclination of the yield curve, respectively (Fig. 2).
S-CLAY1 incorporates two hardening laws. The first describes the change in size of the yield curve, which is assumed to be related solely to plastic volumetric strain (as in MCC) where t is the specific volume, k is the slope of the postyield compression curve in the t-lnp 0 plane for a constant g Fig. 2 Yield curve of the S-CLAY1 model using triaxial invariants and visualization of new invariant q stress path involving no change in anisotropy (e.g. isotropic loading of an isotropic sample), and j is the slope of the swelling line in the compression plane. The second hardening law (rotational hardening) describes the change in inclination of the yield curve produced by plastic straining, both volumetric and shear strains.
where x is a material constant that controls the absolute effectiveness of plastic strains in rotating the yield surface towards the target value. Similarly, x d controls the relative effectiveness of the deviatoric plastic strain, de p d ; and the volumetric plastic strain, de p v . This constitutive model is rate-independent and does not consider interparticle bonding, and the elastic behaviour within the yield surface is isotropic.

Definition of invariants
Current cavity expansion solutions for Cam clay models are formulated using p 0 and q because they do not consider material anisotropy. Some researchers, such as Collins and Yu [8] and Cao et al. [3], use simplified definitions of . For the cylindrical cavity expansion problem in a Cam clay material, Chen and Abousleiman [7] are apparently the first ones to use the full 3-D definition: For the sake of simplicity, S-CLAY1 has been presented in the previous section using the simplified triaxial compression stress space and p 0 and q. However, the cavity expansion problem in plane strain conditions produces stress paths out of the triaxial compression stress space and causes the yield surface to rotate out of it, and the use of q is no longer valid (as will be explained below).
The yield surface of the model (Fig. 2) can be expressed in generalized form as where and Although some Lode's angle dependency may be incorporated in the model, here a constant value of M is considered, and consequently, the critical state surface in the stress space coincides with the Drucker-Prager yield criterion. The polar coordinates are principal directions for this problem, and consequently, the shear components of both stress and fabric tensors are not considered (Eqs. 11 and 12).
In this paper, a new invariant q is proposed instead of q to derive a mathematical formulation for the cavity expansion problem in anisotropic plastic materials. This new invariant makes the derivation of the solution possible without modifying the constitutive model.
and s i are the following deviatoric stresses and a d i are deviatoric components of the fabric tensor Then, the yield surface may be expressed in a similar form as isotropic Cam clay models To understand the necessity of defining the new invariant q, it is useful to study the cross section of the yield surface with the p-plane (hydrostatic or constant p 0 plane). For isotropic Cam clay models, the yield curve in the p-plane is a circle centred in the origin (p 0 axis) and q is the radial distance to the origin, which coincides with the radius of the yield surface (Fig. 3). However, for S-CLAY1 the yield surface is centred in the a axis and the radius of the yield surface is no longer q (Eq. 9), but q (Eq. 13).

Analytical solution 4.1 Elastoplastic stiffness matrix
The usual decomposition of the strain increment vector de in an elastic, reversible contribution de e , and a plastic, irreversible contribution de p , is used.
The increment of elastic strain, de e , is expressed in terms of the effective stress using Hooke's elastic constitutive law where Young's modulus E is defined in terms of shear modulus G and Poisson's ratio m In S-CLAY1 model, G depends on the current stress state and is given by As the model considers an associated flow rule, the three components of the plastic strain increment, de p , are where de p r , de p h and de p z are plastic strain increments in r, h and z directions, respectively, and K is the plastic multiplier. Derivatives of the yield function in terms of stresses To derive the plastic multiplier, the consistency condition (df y ¼ 0) can be applied to the yield surface such that stresses cannot exist outside the yield surface or in terms of plastic strains where de p v and de p d are obtained under the associated flow rule: By substituting Eqs. (25) and (26) into Eq. (24), the plastic multiplier K can be derived as The plastic multiplier can be rewritten in a matrix form as where All required derivatives are given in Appendix 1. By substituting Eq. where notations in Eq. (29) are defined as Equations (19) and (29) Since the strains will be first obtained and then the constitutive equation will be used to get the effective stresses, inversion of Eq. (31) where 3 5 :

Rotational hardening rule
The rotational hardening rule of S-CLAY1 [33] gives the change in the fabric components (da d r , da d h and da d z ) By substituting Eqs. (25) and (26) into Eq. (33), the changes of the fabric components in terms of the plastic multiplier are obtained. where The partial derivative of Eq. (34) with the radial direction provides the changes in the fabric components with the radial direction

Governing equations
As the solution for the elastic zone is already known (e.g. [34]) and given in Appendix 2, here, the governing equations are derived just for the plastic zone. The deformation in the plastic zone should be considered as a large-strain problem (e.g. [8]); so, the radial and tangential strain increments can be defined in natural strain form as where r and dr are position of a material particle in the radial direction and change in the position of that particle, respectively.
Under undrained and plane strain conditions, the volumetric and vertical strains are zero, i.e. de v ¼ de z ¼ 0.
In addition, the change in anisotropy, which is given by Eq. (37), should also be solved simultaneously. By substituting Eq. (41) into Eq. (37), the following differential equations are developed.
So, Eqs. (41) and (42) provide the system of six firstorder ordinary differential equations (ODE) that governs the problem in the plastic region. To solve the system, initial values (i.e. boundary conditions) and numerical integration are required. Initial values are those corresponding to the elastic/plastic boundary (as provided in the following section) and numerical integration is performed in the radial direction from r xp to r x . Here, r x is the position of any particle located in the plastic zone, and r xp is the position of that particle when it was just entering into the plastic state.

Elastic/plastic boundary
The initial values at the elastic/plastic boundary that have to be determined are the position of that boundary, r xp , the corresponding stresses r 0 rp , r 0 hp , and r 0 zp , and the corresponding fabric tensor, which is the initial one (a r0 , a h0 , a z0 ) as it does not change in the elastic zone. The elastic/plastic boundary (i.e. initial yield surface) can be defined using an isotropic overconsolidation ratio, R, in terms of mean effective stresses as where p 0 m is a preconsolidation mean stress (see Fig. 2) and p 0 i is an initial mean stress that may be obtained using the initial stress state (p 0 0 and q 0 ) and the yield surface (Eq. 5). Note that R and the traditional overconsolidation ratio (OCR), which is expressed in terms of effective vertical stresses, are interrelated.
By using the coefficient of earth pressure at rest, K 0 , the initial effective stresses can be written as Deviatoric stress at initial yielding, q (Eq. 9), can be derived as shown in Fig. 3 as Using the initial stress state (K 0 and p 0 0 ), R and q (Eq. 45), the stress state at the elastic/plastic boundary (i.e. initial yielding) can be derived as From the radial displacement given by Eq. (71) in Appendix 2, the position of the material particle at the instant when the particle becomes plastic, r xp; can be obtained as where r x0 is the initial position of the particle and can be obtained using Eq. (4). So, r xp in terms of r x , a and a 0 can be determined substituting Eq. (4) into Eq. (49) The location of the current elastic/plastic interface r p can be obtained by equating both r xp and r x to r p in Eq. (50) because it corresponds to the particle that is just entering the plastic zone right now.

Excess pore pressures
The excess pore pressure, Du, at point r x can be calculated applying equilibrium of radial stresses, i.e. by integrating Eq. (2) from the elastic/plastic interface up to the point r x as

Solution procedure
Stresses, changes in fabric anisotropy, excess pore pressures around the cavity and other results presented below were obtained using a standard differential solver available in GNU Octave v4.0. An existing solver 'lsode' was utilized to solve the system of six first-order ordinary differential equations derived in Eqs. (41) and (42). Figure 4 summarizes the solution procedure used here to solve the cavity expansion problem.

Boston blue clay
To validate the analytical solution and to illustrate its application, Boston blue clay was chosen since it is a welldocumented clay that has already been used in previous cylindrical cavity expansion studies (e.g. [7,22]). The soil parameters are detailed in Table 1. For the sake of comparison, basic Cam clay parameters are taken from Chen and Abousleiman [7]. Additional anisotropy parameters are deduced following Wheeler et al. [33]; for example, the initial inclination of the yield surface, a 0 , is deduced from the stress history of the soil (one-dimensional compression) to fit the coefficient of lateral earth pressure at rest for normally consolidated conditions, K 0NC . OCR of Boston blue clay varies with depth; so, the corresponding state parameters are shown in Table 2. K 0NC is estimated using Jaky's expression: where / is the friction angle of the soil in triaxial compression. For higher OCR values, K 0 is estimated by numerical simulation of the corresponding loading and unloading process starting from the normally consolidated case. G 0 is calculated using Eq. (21), and the undrained shear strength for triaxial compression, c u,TX , is also analytically obtained from the previous values using the following expression, which was derived using a similar procedure as that used by Potts and Zdravkovic [21] for MCC In this case, failure is reached for plane strain conditions and, consequently, the undrained shear strength is defined, according to Tresca failure criterion, for a zero Lode's angle, i.e. plane strain conditions (e.g. [22]): However, for this anisotropic model, Eq. (55) only holds for normally consolidated conditions. When the soil is overconsolidated (OCR [ 1), the mean effective stress at critical state is different in triaxial compression (h = -308) and in plane strain conditions because the q values are different at yielding (horizontal path from the initial point) (see Fig. 3). Consequently, Eq. (55) is no longer valid and c u,PS values had to be obtained numerically ( Table 2).

Validation
To validate the semi-analytical solution, finite element simulations have been performed using the commercial code PLAXIS 2D 2015 [2]. The S-CLAY1 model has been implemented as user-defined soil model in PLAXIS, using Fig. 4 Solution procedure for solving ordinary differential equations of cylindrical cavity expansion in GNU Octave an automatic substepping in combination with a modified Newton-Raphson integration scheme [26,27].
The geometrical modelling of the cylindrical cavity expansion problem (Fig. 5) is based on that proposed by Burd and Houlsby [3] using a correcting layer. That is necessary because the semi-analytical solution assumes a material of infinite extent. In order to model that using a mesh of finite dimensions, a correcting layer is added to the perimeter of the mesh. Burd and Houlsby [3] show that the properties of the correcting layer should satisfy the following relationship where G is the shear modulus of the material, G c and m c are the shear modulus and the Poisson's ratio of the correcting layer, and r c and r e are the outer and inner radius of the correcting layer, respectively.
Here, common values of m c = 0.3, r e = 65a 0 and r c = 2r e were assumed and G c was calibrated using Eqs. (21) and (56). Besides, G c was varied within a limited range to confirm that the results were not affected.
In the analytical solution, the clay is perfectly incompressible, while in the numerical simulations a high but finite value of the bulk modulus is considered. Parametric analyses were performed to confirm the negligible influence. Only, when very high values of the bulk modulus are used, the quality of the stresses is poor (they oscillate).
A prescribed displacement (a-a 0 ) was imposed at the cavity, and initial homogeneous stresses and clay properties were input. To account for large displacements, the numerical code uses an updated Lagrangian formulation [18] and adopts the co-rotational rate of Kirchhoff stress (also known as Hill stress rate). The details of the implementation can be found in Van Langen [28].
Comparison between finite element simulations and the semi-analytical solution gives negligible differences. For illustrative purposes, finite element results are shown in Fig. 6a. To avoid duplication, finite element results have not been included in other figures.
7 Results and discussion 7.1 Internal cavity pressure The internal pressure (radial stress) necessary to expand the cavity, r a , is one of the important variables of the problem. Its value increases as the cavity is expanded. When the cavity has been notably expanded (a/a 0 [ 2), r a approaches an asymptotic limit value, sometimes called pressuremeter limit pressure. Figure 7 shows its variation with the normalized cavity radius for different OCR values. For the sake of comparison, ambient pore pressures are not included in r a . Excess pore pressures at the cavity wall are also depicted in Fig. 7. For high OCR values, slight negative excess pore pressures could be generated at the beginning of the cavity expansion (small a/a 0 values). The final limit values of r a and Du(a) decrease with the OCR in a roughly logarithmic way as proposed by Randolph et al. [22].
The difference between r a and Du(a), i.e. the radial effective stress r 0 r , quickly reaches a constant value that is independent of the OCR ðr Þ c u for a=a 0 [ 1:3Þ as explained in the next section.  7.2 Stresses around the cavity Figure 6 shows the stresses around the cavity for different OCR values, namely 1, 1.5 and 5, and when the cavity radius is twice the initial one (a/a 0 = 2). The stresses are normalized by the undrained shear strength of each case for plane strain conditions ( Table 2). Near the cavity, the stresses are constant with the radius because they are at critical state (CS). The extension of the CS region is smaller in this anisotropic solution (r/a \ 2), than in previous isotropic solutions (r/a [ 2), because important plastic strains are necessary to rotate the yield surface until its CS position. The CS region has been determined assuming a tolerance of 0.1% in any stress change and its extension is around r/a = 1.6 for any OCR. For OCR = 1, all the material points yield just when the cavity expansion starts, so there is no elastic region. The size of the elastic region increases with the OCR value. For high OCR values (for example, OCR = 5), there are negative values not only of the excess pore pressure but also of the hoop effective stresses. Attention should be paid to these negative effective stresses because they may be not realistic depending on the allowed tensile stress in the soil.
The effective stresses at CS may be analytically obtained as a function of the undrained shear strength for plane strain conditions. The vertical stress is the mean value of the radial and hoop stresses because plane strain conditions hold, and given that q f = H3c u = MÁp' f , the following values are obtained at CS:

Stress paths
To provide an understanding of the cavity expansion problem and the role of the fabric anisotropy of the clay is important to analyse the stress paths followed by a point at the cavity wall. Points at further distances follow the same path, but if they are outside the CS region, they stop earlier. Figure 8 shows the effective stress paths (ESP) in the p 0q diagram for OCR = 1, 1.5 and 5. As is common practice, the intersection of the yield surface (YS) with the triaxial plane is also depicted in Fig. 8. It is worth noting that for isotropic yield surfaces, there is a unique representation in the p 0 -q diagram, but for anisotropic yield surfaces, it depends on the intersection plane. As the stress paths go outside the triaxial plane, the intersection of the yield surface with the triaxial plane, as plotted in Fig. 8, is meaningless. For example, the stress paths in Fig. 8b, c go vertically upwards until they reach the yield surface and they bend when the yield surface is reached. As the yield surface is reached outside the triaxial plane, the q value is lower than expected from the intersection of the yield surface with the triaxial plane plotted in Fig. 8. Rotation of the yield surface also plays an important role. To perceive that, another type of stress space should be used, e.g. the pplane (Fig. 9). The cross section of the yield surface with the p-plane is a circle, whose centre is in the a-axis (ap 0 ) and its radius is q. Within the elastic region, the ESP goes straight to the right until it reaches the yield surface (Fig. 9b, c) because r 0 z does not change and the amount that r 0 r increases is the same as r 0 h decreases. Once the ESP touches the yield surface, it bends towards plane strain conditions (h = 08) r Á and the yield surface also rotates. The final point, i.e. at critical state, corresponds to q f = H3c u , is in plane strain (h = 08) and at the rightmost point of the yield surface in the p-plane. The centre of the yield surface, i.e. the a-axis, is also in the plane strain plane (h = 08), and the model (S-CLAY1) predicts a unique inclination of the yield surface at critical state, namely a = M/3. In Fig. 9a, the path followed by the aÁp 0 vector is also plotted to highlight how the yield surface rotates towards a zero Lode's angle (h = 08).
When the K 0 line is below the a-axis, i.e. high OCR values (e.g. Fig. 9c), the ESP goes upward a bit after reaching the yield surface because it goes towards the

Evolution of fabric anisotropy
One of the main capabilities of the used constitutive model (S-CLAY1) is the possibility to reproduce the evolution of fabric anisotropy with plastic straining. So, fabric anisotropy changes during the cavity expansion process as already shown in Fig. 9. If just the inclination of the yield surface changes, the analysis may be reduced to the scalar value a, but if the yield surface also rotates as in this case, it is necessary to analyse the fabric tensor a. Figure 10 shows the variation of the scalar value a (inclination of the a-axis in a p 0 -q diagram) and the different components of the fabric tensor, a i . The inclination of the yield surface at critical state is a = M/3, and since this inclination corresponds to plane strain conditions (a z = (a r ? a h )/2), the fabric tensor at CS, i.e. near the cavity, is a r a h a z ½ The results for different OCR only vary in the plastic region, depending on its extent.
To help to visualize the changes in the soil fabric, Fig. 11 shows the (a z -a h , a r -a h ) vector in arbitrary points. This vector tries to resemble the orientation of the clay platelets, and in this problem, it changes from horizontal direction for an initial vertical cross-anisotropy towards a nearly vertical one for radial cross-anisotropy after cavity expansion.

Influence of anisotropy
The initial inclination of the yield surface (a 0 ) has been varied to analyse its influence. Between the limit cases of For an initially inclined yield surface (a 0 = 0.46), the radial extent of the excess pore pressures is larger as the yield surface quickly starts to rotate and at the cavity boundary its normalized value is consequently higher. On the other hand, the rotational hardening rule causes that in both cases the final inclination of the yield surface is the same.
To really evaluate the influence of considering plastic anisotropy, the rotational hardening rule must be also deactivated (x = 0), yielding into the MCC model. For the MCC, the present solution exactly coincides with the solution by Chen and Abousleiman [7]. Comparison between MCC and S-CLAY1 is depicted in Fig. 13. When considering plastic anisotropy, the extend of the excess pore pressures is larger as the yield surface quickly starts to rotate, but, on the other hand, this rotation reduces pore pressure generation and the maximum excess pore pressure at the cavity wall is lower. Effective stresses at CS are exactly the same as already justified. In the plastic zone, there are small differences, for example, the hoop stress monotonically decreases for MCC, while for S-CLAY1 it decreases and slightly increases near the CS zone due to rotation of the yield surface and the shape of the rotational hardening law (Eq. 7) [33].

Case without rotational hardening
The constitutive model used to develop the analytical solution (S-CLAY1) is able to reproduce not only the initial anisotropic behaviour of the clay but also its evolution with plastic straining. The present problem (cavity expansion) is a very good example of the importance of considering the changes in fabric anisotropy when the loading path notably differs from the loading path that created the initial anisotropy (one-dimensional compression). To highlight the differences, the rotational hardening law is deactivated, but the soil is initially anisotropic (x = 0) ('no RH') and the results compared with the reference case where rotational hardening was activated ('with RH') in the anisotropic soil. Figure 14 shows the stress path at the cavity wall when rotational hardening is deactivated (OCR = 1, a 0 = 0.46). As for the case with RH, the stress path goes towards a zero Lode's angle (h = 08) (Fig. 14a). However, the stress path does not reach h = 08 because the yield surface is inclined in the vertical axis and it cannot rotate towards h = 08. So, the stress state must remain within the vertically inclined yield surface, which is not realistic. The final point of the stress path is at the rightmost point of the yield surface in the p-plane, which has the shape of a circle. In this case, the effective mean stress also decreases (Fig. 14b) to allow the final point to be closer to h = 08.
Neglecting fabric evolution with plastic strains in clays leads to unrealistic results (Fig. 15), particularly in the most important one, the radial stress, which is notably underpredicted if rotational hardening is not allowed. So, cavity expansion problems in natural clays should not be analysed with a plastic anisotropic model that does not consider rotational hardening, as done in Li et al. [17].
Since the formulation of the yield surface in Li et al. [17] is different from the S-CLAY1 model, the predictions are expected to be different from the present ones.

Conclusions
A semi-analytical exact solution for the undrained expansion of a cylindrical cavity in clays with fabric anisotropy is rigorously developed. The adopted constitutive model, namely S-CLAY1, is a Cam clay model that reproduces plastic anisotropy, both initial fabric and its evolution with plastic straining. To develop the solution, a new stress invariant q, which represents the radius of the yield surface in the p-plane, is introduced. The solution involves the numerical integration of a system of six first-order ordinary differential equations, three of them corresponding to the effective stresses in cylindrical coordinates and the other three to the components of the fabric tensor. The semi-analytical solution is validated against finite element analyses, using Boston blue clay as the reference clay. After cavity expansion, three zones may generally be S-CLAY1 predicts a unique fabric tensor at CS for plane strain conditions, which may be analytically obtained. The initial vertical cross-anisotropy caused by the soil deposition and consolidation changes towards a radial cross-anisotropy after cavity expansion. Normalized maximum excess pore pressures at cavity wall are slightly lower for S-CLAY1 than for MCC. Analytical solutions for soils with fabric anisotropy must consider fabric evolution with plastic straining because otherwise stresses near the cavity are not usually at a zero Lode's angle, which leads to unrealistic results.