Analysis of Soil-Compacting Effect Caused by Shield Tunneling Using Three-Dimensional Elastoplastic Solution of Cylindrical Cavity Expansion

1Department of Civil Engineering, Shanghai University, 99 Shangda Road, Shanghai 200444, China 2Department of Civil Engineering, Shanghai University, 99 Shangda Road, Shanghai 200072, China 3Department of Civil and Environmental Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong 4Department of Engineering, College of Engineering Mathematics and Physical Sciences, University of Exeter, Exeter, Devon EX4 4QF, UK


Introduction
Theoretical analysis of cylindrical cavity expansion theory in geotechnical engineering has been widely used in pile sinking, shield tunneling, static cone penetration test (CPT), pressuremeter test (PMT), and mixing pile construction disturbance problems [1][2][3][4][5][6].Hill (1950) first proposed the spherical cavity expansion method and derived a general solution of stresses and displacements in Tresca material [7].Vesic (1972) obtained the classical solutions considering the plastic zone volume change of ideal elastic-plastic expansion problem based on Mohr-Coulomb model [8].Cao et al. (2001) presented closed form solutions for both spherical and cylindrical cavities by assuming the ultimate deviator stress distributions in the plastic zone and taking the in situ stress conditions into consideration [9].Alonso  obtained the self-similar solution for the circular tunnel in strain-softening rock masses [10].Wheeler et al. (2003) and Nakano et al. (2005) adjusted yield surface equation on the basis of  0 -based modified Cam-clay, considering  0 consolidation caused by induced anisotropy and its evolution law in the loading process [11,12].Yin and Hicher (2008) identified parameters of modified Cam-clay model with viscosity of soil from the cavity expansion [13].The problem of cavity expansion and cavity contraction has attracted much attention in geotechnical problems with application to the bearing capacity of deep foundations, interpretation of pressuremeter tests, breakout resistance of anchors, pile driving, wellbore instability, underground excavation, and blasting fracturing by explosives.
However, based on the above-mentioned researches, the conventional cylinder expansion theory assumes that initial stress is isotropic (static earth pressure coefficient  0 = 1.0).Due to the sedimentation history, mineralogical characteristics, and other several factors, the initial stress of soil layer is usually anisotropic.For example, the vertical stress and horizontal stress of the cylinder bore are not the same in the horizontal tunneling construction of the underground pipe laying, tunnel engineering.As a result, the traditional theory of cylindrical expansion generally used in practical engineering fails in explaining some practical phenomenon.Also, the traditional cylindrical expansion theory assumes that the cylindrical boundary condition is controlled by displacement, which presents around the cylinder in the form of symmetrical distribution as a function of distance.The traditional control theory can be applied in geotechnical problems analysis such as cylindrical pile and static cone penetration test.But in the practical engineering of shield tunnel, boundary condition of the cylinder expansion is controlled by stress and displacement of the cylinder, which was unable to satisfy the condition of axial symmetry.Hence, it may not be reasonable to explain such geotechnical problems considering traditional theory of cylindrical expansion.
In this paper, a three-dimensional elastoplastic solution of cylindrical cavity expansion is obtained based on the theory of critical state soil mechanics and engineering characteristics of shield tunneling to calculate the undrained cylindrical cavity expansion considered anisotropic clay of  0 .This research work not only contributes to the theoretical values in field of geotechnical engineering but also has an important practical engineering significance.

Definition of Soil-Compacting Effect and Mechanical Model
2.1.Definition of Soil-Compacting Effect.Shield tunneling is a typical three-dimensional problem.Soil deformation is closely related to the relative position of the shield machine.
In actual shield construction, earth pressure and support pressure ahead of excavation face are not completely balanced [15].When support pressure is greater than passive earth pressure, the soil is extruded and then soil-compacting effect presents.And the free surface may collapse if support pressure is less than active earth pressure.This paper mainly aims to study the soil-compacting effect produced when the support pressure is larger than passive earth pressure.Soilcompacting effect of shield tunneling is shown in Figure 1.
(a) Along with the shield approach, the soil, which is squeezed by the excavation face, will be pushed around the shield.And then, the soil is output by the screw conveyor to form a void.The surrounding soil is excavated and unloaded, the stress is released to cause soil expansion, and the front soil moves to the void (Figure 1(a)).
(b) As the shield machine moves forward, the soil will be squeezed around the shield and excavation face.In the plane perpendicular to the axis of the tunnel, the soil is assume subjected to uniform radial extrusion force, which generates an equal amount of outward radial movement of the tunnel (Figure 1(b)).
(c) The extruded soil moves outwards of the shield machine and is also affected by the friction force arise from the driving process of shield machine.Radial deformation along the outer side of the tunnel finally formed under the influence of initial stress field (Figure 1(c)).It can be seen that the squeezing effect of shield construction makes the surface subsidence or uplift and layered soil movement, soil stress, water content, pore water pressure, elastic modulus, Poisson's ratio, strength, and bearing capacity of the physical and mechanical properties are likely to change.It is feasible to simulate the soil-compacting effect of shield tunneling by using cylindrical cavity expansion theory.

Mechanical Model.
Figure 2 shows a cylindrical cavity expansion in an infinite  0 consolidated saturated clay with an in-plane initial horizontal pressure stress  ℎ =  0  V and vertical pressure stress  V at infinity, as well as the hydrostatic pressure  0 .Note that cylindrical polar coordinate system is adopted for present analysis, and the occurring compressive stresses and strains are taken as positive.As average internal pressure gradually increases from  ℎ to   , the soil adjacent to the cavity wall yields first.Further increase in average internal cavity wall pressure will lead to current cavity radius to the elastic-plastic (EP) boundary   with a further increased internal cavity pressure with the continuation of the shield excavation.The symbol  0 represents the initial position of the soil when the soil becomes plastic state.
Under undrained condition, relationship between current radius   , initial radius  0 of a material particle, and current radius , and initial radius  0 of the cylindrical cavity is as follows: At any stage of the shield excavation, any soil element within the surrounding soil mass satisfies the following equilibrium equation in both elastic and plastic regions (effective stress form): or alternatively in the total stress form where    and    are effective radial and tangential stresses, respectively;   and   are total radial and tangential stresses, respectively;  is pore water pressure.

Constitutive Relations.
The  0 consolidated clay is assumed to be linearly elastic and infinitesimal deformation until the onset of yield.Thus, according to the Hook's Law, the elastic stress-strain relationship can be expressed in the increment form as follows: where  and  are the elastic modulus and the shear modulus, respectively, which are defined in the  0 -MCC model as follows [6]: where  is the slope of swelling line,  is the specific volume, and ] is the Poisson's ratio.After yielding, elastoplastic behavior of soil can be described by large deformation theory and  0 -MCC model, where a relative stress ratio  * (=  −  0 ) is adopted instead of the stress ratio (= /  ) as used in the MCC model to consider the effect of initially stress-induced anisotropy on mechanical behavior of  0 consolidated clay.Yield function () for associative plasticity also serves as a plastic potential function () of the  0 -MCC model can be expressed as follows [6,11,12]: where  0 = |3(1 −  0 )/(2 0 + 1)| is stress ratio under  0 consolidated condition,  = 6 sin   /(3 − sin   ) is slope of critical state line,   is effective friction angle,    is effective yield pressure, and stress parameters   and  are defined respectively as follows: where   is Kronecker's delta.After taking strain hardening rule into consideration, yield function or plastic potential function can be given as follows: where  and  are slopes of consolidation and swelling lines, respectively;  0 is initial void ratio of soil for   =   0 ;   V is plastic volumetric strain which is used as a hardening parameter.
As Figure 3 shows, plastic potential surface of  0 -MCC model in   − space is a rotated ellipse, and degree of rotation portrays the extent of anisotropy.When shearing starts from isotropic stress  * =  = /  because  0 = 0, the  0 -MCC model becomes the MCC model.The yield surface rotation can also be found in more detail in Yin et al. (2011Yin et al. ( , 2015) ) [16,17].
For elastoplastic constitutive model, the total strain increments   are decomposed into elastic and plastic strain increments; i.e., where elastic strain increments    can be still calculated using (3), and plastic strain increments    in the  0 -MCC model can be derived by taking account the associated plastic flow rule as follows: where the scalar multiplier Ω is given as follows: where Substituting ( 13)-( 16) into ( 18) yields the following equation: where   = ( − )/(1 +  0 ).Substituting ( 18) and ( 13)-( 18) into (11), the scalar multiplier Ω can be obtained as follows: Then, detailed relationship between the plastic strain increment    and the stress increment   can be expressed as follows: (21)

Elastic Zone. Based on the small-strain theory and
(2) and ( 3), stress and displacement in the spherical coordinates of the elastic region can be obtained as follows [9,14,18]: where   is the total radial stress in the phase plane at EP boundary.It is worth noticing that there is no excess pore water pressure in plastic region, because effective mean stress is constant under undrained condition ( = 3(1 − 2V)  / = 0) and total mean stress, , also keeps constant in elastic zone.

Elastoplastic Zone.
For the problem of shield excavation, soil position around the tunnel cavity can be described as , , and  in cylindrical polar coordinate system.Correspondingly, there are only variables of principal stress and principal strain without the deviator stress and deviator strain induced during shield excavation process.Therefore, elastoplastic constitutive relation of the  0 -MCC model can be simplified to a much simple matrix equation: where Taking the inverse calculation to (27) gives where For undrained cylindrical cavity expansion,  V =   ; thus, according to the large deformation theory, radial and tangential strains can be given in incremental form as follows: Substituting (32) into (29), constitutive matrix can be reduced to a set of first-order ordinary differential equations as follows: which can be solved by Lagrangian analysis method as an initial value problem with  starting at   , provided that initial values of    (  ),    (  ), and    (  ) are given.Here   denotes the position of a specific particle which comes into plastic state instantly.All of the initial values required for solving the differential equations can be derived from the EP boundary conditions, which are further analyzed as follows.
3.3.Elastic-Plastic Boundary Conditions.As shown in ( 22), ( 23), (24), and (26), the change of mean effective stress, Δ  , is equal to zero in elastic region, which is still valid at the EP boundary; i.e., where    denotes mean effective stress at the EP boundary.Substituting (36) into ( 6), the deviator stress at the EP boundary,   , can be obtained as follows: where OCR is overconsolidation ratio, defined as   0 /  0 , and pressure,   0 , is the maximum mean preconsolidation stress.According to (7), the deviator stress at the EP boundary can also be written as where    ,    , and    are effective radial, tangential, and vertical stresses, respectively, at the EP boundary.
For the  0 consolidated clay, relationship between in situ horizontal effective stress,   0 (or   0 ), and in situ vertical stress,   0 , can be given as follows: where symbol  0 denotes the anisotropic degree of clay.Substituting (39), (40), and (41) into (38), effective radial and tangential stresses at the EP boundary can be derived as Mathematical Problems in Engineering 7 The position,   , where the specific particle becomes plasticity, can be derived by taking the displacement at EP boundary,   , and the undrained condition into consideration.According to (25), one can obtain Substituting ( 1) into (43) gives [14] Furthermore, current location of the EP boundary,   , which is required for full analysis of stress distributions around the cavity, can be obtained by equating both   and   to   [18]: So far, (40), (42), and (44) provide initial conditions for solving the first-order ordinary differential equations (33), (34), and (35).
After solving these first-order ordinary differential equations by numerical method, excess pore water pressure at an arbitrary location in the plastic zone, Δ  , can be computed by integrating (2) from the EP boundary   to the point   :

Closed Form 3D Analytical Solution
In the front section (refer to Figure 3), the semianalytical solution of expansion of the cylindrical hole has been presented with the limitation of high computation cost.In order to reduce the computational cost and achieve the purpose of practical application, this paper deduces a closed 3D total stress solution for cylindrical cavity undrained expansion in  0 consolidated clay.Under the undrained condition, the total volume of the soil is equal to zero.That is, Taking EP boundary conditions into consideration in (47), the following equation can be found: According to (37) and (48), the relationship between average effective stress and the plastic zone can be presented as follows: where Λ is the plastic volumetric strain ratio, defined as (1 − /).
For the  0 -MCC model, critical state condition can be expressed as follows: Substituting (50) into (49), ultimate effective mean stress,    , and ultimate deviator stress,   , can be obtained as follows: On the other hand, the closed form 3D solutions in elastic region have been already derived, as shown in ( 22)-( 26).Only the solutions in plastic region are derived in this study.Total radial stresses at an arbitrary location in plastic region,   , can be obtained by integrating (4) from the EP boundary,   , to the point,   , where   is the total radial stress at the EP boundary.The deviator stress, , is related to the location   , so that (53) cannot be integrated directly.Therefore, in order to obtain the closed form 3D solution, the following two assumptions have to be made at first.

(a) Relationship between Deviator Stress and Radial Distance
in the Plastic Region.The deviator stress, , in the plastic zone is assumed to be equal to the ultimate deviator stress,   , which is similar to .
(b) Relationship between and among Total Radial, Tangential, and Vertical Stresses in the Plastic Region.The vertical total stress   is assumed to be the mean of radial stress   and tangential stress   after soil yields.Although these assumptions seem to be bold and strong, the approximate analytical results match the exact numerical results fairly well in the plastic region, which will be compared and discussed in detail later in the paper.
According to assumption (a), substituting (25) into (45) and ignoring higher order terms of √ 2  − (1 −  0 ) 2  2 V0 /, position of the EP boundary, which is related to the current cavity radius, can be written as Then, taking assumption (a) into consideration, radial stress in the plastic zone can be derived from (53) and can be written as follows: According to the previous two assumptions, substituting (55) into (8) and taking (52) into consideration, total tangential and vertical stresses in the plastic region can be obtained:

𝑥
) . (56) Furthermore, excess pore water pressure in plastic region can be derived by effective stress principle: Equating   to , pressure-expansion relationship and excess pore water pressure at cavity wall can be derived from ( 55) and (57), respectively, as follows: From (58), as the cavity continues to expand to infinity (/ 0 → ∞), the limit cavity wall pressure,   , and excess pore water pressure, Δ  , for undrained cavity expansion in  0 consolidated anisotropic clay can be obtained: ] .
(61)  [9].Especially for isotropic clay ( 0 = 1), the closed form 3D solutions presented in this paper reduce to the solutions obtained by Cao et al. [9], which demonstrates that the solution obtained by Cao et al. [9] is just a special case of the presented solution.Actually, for nature clay, the horizontal in situ stress is different from the vertical one, and the soil is initial anisotropy (i.e.,  0 ̸ = 1).Note that the term of √ 2  − (1 −  0 ) 2  2 V0 in the presented solution reflects the effect of initially stress-induced anisotropy on shield excavation.Therefore, presented solutions can yield more reasonable results for cylindrical cavity expansion in anisotropy clay.In addition, the term of ( 0 − 1)  0 /(1 + 2 0 ) in ( 57), (59), and (61) implies that initial anisotropy also has notable effects on excess pore water pressure during shield excavation process.

Calculation Result Analysis and Discussion
In order to verify the accuracy and suitability of this method considering the influencing factors of initial anisotropy and induced anisotropy, calculation results of present method and Chen and Abousleiman's solution [18] using MCC constitutive model results were, respectively, validated in this paper.
Physical and mechanical parameters of soil used in this paper are taken from Chen and Abousleiman [18].Table 1 shows the properties of three typical  0 consolidated clays.

Pressure-Expansion Relationship.
As shown in Figure 4, comparison of pressure-expansion relationship around shield tunnel is calculated by the proposed 3D closed form solution and Chen's solution.All calculated stresses have been normalized.Effective radial stress and vertical stress have been overestimated, while the tangential stress is noted to be underestimated in the plastic zone based on Chen's solution apart from  0 = 1.Therefore, the anisotropy induced by initial stress has a significant influence on the stress distribution in the plastic zone.If the tunnel is in the isotropic stress conditions, that is to say  0 = 1,  0 -MCC constitutive model will degenerate into a MCC constitutive model; in this case, the results of the two calculations are in complete agreement.
It is also interesting to note in Figure 4 that vertical stress is always equal to the average of radial and tangential stresses in critical state zone for all cases, which can be proved exactly as follows.Combining (33)-(35) gives Since the stress components keep constant in the critical region, therefore, the left side of (62) is equal to zero; thereby (62) can be reformulated as Equation ( 63) indicates that the second assumption for approximate closed form 3D solution holds in the critical region.
It can also be seen from Figure 4(d) that excess pore water pressure is underestimated near the cavity wall but overestimated as the radial distance increases.Besides, the excess pore water pressure keeps constant in critical region, and thus the second assumption for the approximate closed form 3D solution still holds well in the critical region.When the shield is compacted to a certain degree, excess pore water pressure in the soil squeezing effect zone is gradually stabilized, and the same level of pressure can be achieved, regardless of the value of  0 .

Shield Excavation Process Analysis.
According to cylindrical cavity expansion theory, effect of soil compaction caused by shield tunneling is mainly revealed in dynamic change of the normalized excess pore water pressure Δ  /  0 , the soil pressure inside the cavity   /  0 , and the cavity radius   /. Figure 5 shows variations of normalized cavity pressure, soil pressure, and excess pore water pressure with cavity radius during shield excavation process.Results computed by the three different typical values of  0 (0.625, 1, and 2) [14], as well as the results from the MCC model based solutions, are presented in the same figures for comparison.Similar to previous analysis, all calculation results have been normalized.
It can be seen from Figure 5 that excess pore water pressure and internal cavity pressure are affected by the initial stress-induced anisotropy significantly when / 0 > 2. In the process of gradual expansion of inner wall of the tunnel, the result of MCC constitutive model is smaller than that of the  0 -MCC constitutive model.Stress values at the inner wall of the tunnel and excess pore water pressure around the tunnel increase dramatically when / 0 < 2. Both parameters gradually achieve a stable level as the shield continues to advance.Compared to Figures 5(a), 5(b), and 5(c), it is clear that, as the  0 value increases, internal cavity pressure caused by shield tunneling is effectively suppressed, while the impact of excess pore water pressure and plastic radius is not very obvious.This is because, at the initial stage of pore expansion, negative pore pressure first occurs at the hole wall.Then, with the increase of pore size, the excess pore pressure will turn to a positive value from negative to a critical state.Furthermore, when the surrounding soil is in isotropic initial stress state condition ( 0 = 1), the  0 -MCC constitutive model is reduced to the MCC constitutive model, and the results are consistent with numerical results.However, effect of anisotropy degree on discrepancy of excess pore water pressure calculated from the two solutions is not obvious, which implies that the OCR also has a pronounced effect on excess pore water pressure at cavity wall during cavity expansion process.
Evolution of normalized plastic radius,   /, with normalized instant cavity radius, / 0 , for three different values of  0 and OCR, together with the MCC model based results is also plotted in Figure 5 for comparison.It is seen that, for normally consolidated clay with  0 = 0.625, plastic radius during whole cavity expansion process is significantly affected by initially stress-induced anisotropy.However, for heavily overconsolidated clay, discrepancy between presented solution and MCC model based solution is not apparent, because OCR also plays a major role in plastic radius for heavily consolidated clay during cavity expansion process.

Effective Stress
Path in   - Plane.Figure 6 shows effective stress path (ESP) in   - plane for a soil particle located at cavity wall in three typical  0 consolidated clays, where two stress parameters   and  have been normalized by   0 .Note that points O, P, and F in Figure 6 denote the in situ stress point, yield stress point, and failure stress point, respectively.It is also clear that initial yield loci are rotational ellipses for anisotropic clays ( 0 ̸ = 1), and degree of rotation reflects degree of anisotropic.For isotropic clay ( 0 = 1),  0 line overlaps with   axis, and yield locus is symmetric to   axis, and thus  0 -MCC model reduces to MCC model.For anisotropic clays, yield locus and CSL of the  0 -MCC model are higher than that of MCC model.Therefore, yield and failure stresses during cavity expansion process are underestimated by previous MCC model based solutions.
For normally (OCR = 1) consolidated clay with  0 = 0.625 in Figure 6, yielding occurs immediately after cavity expansion, so that in situ stress point O overlaps the yield stress point P. The corresponding ESP starts form point O which is located in  0 line and moves up to the left until ESP reaches critical state line (CSL) at point F. For lightly (OCR = 3) and heavily overconsolidated clay (OCR = 10), ESP firstly moves vertically and soil keeps elastic until the ESP approaches initial yield locus at point P. Then soil yields, and the ESP turns up to the right until reaching CSL at point F.

Conclusions
This paper aims to study the soil-compacting effect produced by shield tunneling by a precise semianalytical solution.Particular emphasis is given to develop an approximate closed 3D solution of cylindrical cavity undrained expansion in saturated anisotropic clay based on the  0 -based modified Cam-clay model.Furthermore, efforts have been made to examine (a) the effects of the initial anisotropic and initially stress-induced anisotropy on shield excavation and (b) stress distribution of shield tunneling under different  0 conditions calculated by the approximate closed form 3D solution.Some of key observation and findings from the study are as follows.
(a) Both numerical solution and closed from 3D solution presented in this study improve the understanding of conventional solutions by considering anisotropic properties of natural clay and reducing MCC model based solutions when  0 = 1.Therefore, presented solutions can yield a more realistic stress field induced by cylindrical cavity expansion in natural clay.
(b) Vertical total stress   , radial stress   , and tangential stress   become intermediate principal stress and major and minor principal stress, respectively, when soil came into critical state, and vertical total stress   is equal to the mean of radial stress   and tangential stress   in critical state region around the cavity for all cases.
(c) Initially stress-induced anisotropy has significant influence on stress field and plastic radius around the cavity, especially in normally consolidated clay.Excess pore water pressure is not only affected by the initially stress-induced anisotropy but also by the initial anisotropy of clay.The stress state of the soil around the borehole wall changes greatly due to the expansion of the column hole in the  0 consolidated clay, and the stress around the hole decays logarithmically with the increase of the radial distance and tends to the stress state at the time of departure.
(d) The initial stress anisotropy and stress history have obvious influence on the expansion process of cylindrical cavity, especially for the ultimate expansion pressure of cylindrical cavity expansion process.Stress field caused by shield excavation is overestimated, whereas plastic radius during cavity expansion process is underestimated by MCC constitutive model based on solutions in anisotropic clays.

Figure 1 :
Figure 1: Schematic diagram of compacting soil by shield.

Figure 2 :
Figure 2: Expansion of a cylindrical cavity in  0 consolidated anisotropic clays.

Figure 4 :
Figure 4: Comparison of pressure-expansion relationships around shield tunnel calculated by 3D closed form solution and Chen's solution.

2 Figure 5 :
Figure 5: Evolution of cavity pressure and excess pore water pressure with cavity radius during shield excavation process.

2 Figure 6 :
Figure 6: ESP in   - plane for a soil particle located at cavity wall during shield excavation process.

Figure 7 :
Figure 7: Stress distribution of shield tunneling under different  0 conditions calculated by the approximate closed form 3D solution.
Tunneling under Different  0 Conditions.Figure 7 shows stress distribution of shield tunneling under different  0 conditions calculated by the approximate closed form 3D solution, where normalized vertical total stress   /  0 , mean of radial stress   /  0 , tangential stress   /  0 , and excess pore water pressure Δ  /  0 are plotted against normalized radial distance / 0 .It is seen when / < 2, the rise of  0 leads to a significant change of stress distribution, which will reach a relatively stable state after EP boundary presents.Besides the value of  0 has no effect on the final stress except the vertical total stress   /  0 .Stress values under different  0 values finally reach the same magnitude near the cavity wall.It also indirectly verifies that the two basic assumptions of the closed form 3D solution is feasible and reasonable.Therefore, the closed form 3D solution is exact enough to interpret and solve general geotechnical problems.