#### Abstract

Soil strength isotropy and coaxiality of the principal axes of stress and plastic strain rate tensors are often assumed in modelling the soils around tunnels, in the conventional plastic theory. However, the complexity of the soil response during the excavation of tunnelling cannot be well exhibited. In this paper, a recently developed two-dimensional (2D) elastoplastic constitutive model incorporating both soil strength anisotropy and non-coaxiality is firstly reviewed and then implemented in the finite element platform ABAQUS through the user-defined material subroutine (UMAT). This model is then numerically applied to analyze tunnel excavation problems. Two case studies, i.e., centrifuge testing and field testing, are carried out and compared with the numerical simulations to investigate the influences of soil anisotropy and non-coaxiality on the stress paths of soils around the tunnel crown and the subsurface settlement induced by tunnel excavations. The results show that the representative soil elements around tunnel crown experience severe principal stress orientations. The predictions of normalized subsurface settlement troughs can be improved by considering initial soil strength anisotropy compared to the conventional Mohr–Coulomb model. A larger value of the non-coaxial coefficient results in a larger magnitude of the maximum vertical displacement. The normalized subsurface settlement troughs are better improved in the case of centrifuge testing than that of the field testing.

#### 1. Introduction

In recent years, with the continuous development and utilization of underground space and the increasing construction of subway tunnels, a number of excavation problems that may occur attract the researchers to work on them [1, 2]. As first summarized by Peck [3] and developed by Ward and Pender [4], the three most significant requirements that satisfy the successful design and construction of a tunnel are stability, ground movements, and their effects and performance of linings. Over the decades, different design methods have been used in engineering practice; namely, simple empirical, analytical, and finite element analyses. Finite element analysis has widely been adopted over the recent years [5]. With the rise of computer capacity, complex geometries, excavation procedures, and support installation sequences of tunnelling can be simulated more realistically. These advantages made the FE method attractive in the realm of tunnelling design practice. However, a notable literature has shown that the numerical results of the subsurface settlement trough induced by tunnelling are too wide when compared with the field data [6]. Several studies have focused on the reasons that could account for this discrepancy. All the parameters that may lead to the discrepancy have to be thoroughly investigated as the major concern in constructing such tunnels in urban areas. The main purpose is to reduce and control the subsurface settlements. These studies can be classified into two categories of scenarios [5]:(1)The selection of proper constitutive models(2)The excavation and support sequence

Gunn [7] identified the reason why finite element calculations are much shallower and wider than those obtained from model tests or observations on site by investigating through a tunnelling research programme at Cambridge University in the late 1960s and early 1970s The reason was quickly identified as the elastic part of any constitutive model, which is used to present the stress-strain behaviour of soil, is assumed to be linear and isotropic. After that, many researchers have been trying to perform the numerical simulations of tunnels with various constitutive models, such as Addenbrooke and Potts [8] used the isotropic and anisotropic, linear, and nonlinear prefailure deformation behaviour as well as loading reversals and Dolezalova [9] used a linear elastic perfectly plastic and a nonlinear elastic perfectly plastic constitutive model. However, current assumptions of constitutive models cannot describe the response of soils with a desired accuracy.

A large number of laboratory tests [10–12] and micromechanics-based virtual tests [13, 14] have observed that when the principal axes of stress rotate in the process of shearing, the directions of the principal plastic strain rate and the principal stress are not coincident. These fundamental findings have been used to guide the development of realistic continuum soil models. In the past, analyzing the boundary value problems (e.g., bearing capacity of slopes, piles, and strip footings) normally assumed the stress and plastic strain rate to be coaxial [15–18], even for the seismic loads where severe principal stress rotations can observed [19, 20]. Although a couple of authors [21–25] have applied their non-coaxial model in analyzing strip footing problems, these models have not widely been applied to analyze tunnelling problems, in particular, in the frame work of soil strength anisotropy. Yang et al. [12] conducted a series of hollow cylinder apparatus (HCA) tests and suggested that non-coaxiality is a significant aspect of soil anisotropy. It is widely acknowledged that in the process of tunnel excavations, there exhibit obvious principal stress rotations, especially in the soil adjacent to the tunnel crown. Hence, it is necessary and attractive to incorporate both non-coaxiality and soil strength anisotropy to the conventional constitute model to describe the response of soils.

In this paper, a recent developed two-dimensional perfect elastic-plastic constitutive model, which incorporates both the soil strength anisotropy and non-coaxiality to the conventional Mohr–Coulomb model, is briefly reviewed. Then the non-coaxial soil model is implemented into finite element software ABAQUS via the user-defined material subroutine (UMAT). As case studies, a series of numerical simulations have been performed with the non-coaxial model employed to analyze the subsurface settlements in tunnel excavation problems. The numerical results have been compared with the centrifuge testing results and the field testing results of the Qinghuayuan Sub-tunnel in Beijing.

#### 2. Non-Coaxial Soil Model with an Anisotropic Yield Criterion and Its Implementation: A Review

The signs of the stress (rate) are chosen to be positive for compression. The total strain rate includes the elastic strain rate and the plastic strain rate and is shown as follows:where the superscripts e and p denote the elastic and plastic strains, respectively.

The general rate equation for an elastoplastic relationship iswhere denotes the elastic stiffness matrix. In our model, the elastic component is assumed to follow Hooke’s law under plane strain conditions. The isotropic elastic stiffness matrix is shown as

Following Yu and Yuan [21] and Yuan [25], the total plastic strain rate is divided into two parts: the conventional component, which is derived from the plastic potential theory, and the non-coaxial component, which is caused by stress rates that are tangential to the yield surface:

To account for the effects of soil strength anisotropy, the coaxial plastic strain rates are defined on the basis of a modified Mohr–Coulomb yield function [26], which can be written in a plain strain stress space as follows [24, 27]:wherewhere and Θ refers to the angle of deviation of the major principal stress direction to the *x*-axis, *c* represents the cohesion, *ϕ*_{max} is defined as the maximum peak internal friction angles, respectively, along all possible major principal stress directions, and *n* and *β* are the parameters newly introduced to represent the degree of anisotropy of soils. Validation of this anisotropic yield criterion can be found in Figure 1.

**(a)**

**(b)**

Following Yu and Yuan [21] and Harris [28], the non-coaxial plastic strain rate is given bywhere *k* is the non-coaxial coefficient and is the material derivative.

By combining the consistency condition equation for perfect plasticity, the general relationship between total strain rates and stress rates can be given as [24, 27]wherein which a modified elastic stiffness matrix is introduced aswhere *I* is the identity tensor and *N* is a tensor depending on the non-coaxial coefficient and stress states.

The above constitutive formulations are implemented into a finite element software ABAQUS as a user-defined subroutine. In the material subroutine, the explicit substepping scheme (an explicit forward Euler/modified Euler pair) with automatic error control is used to integrate the constitutive formulations [29, 30]. The modified regula falsi is used to solve the yield surface intersection problem. For details of the description of the formulations and their verifications and validations, refer [25, 27]. The schematic illustrations of the strain rates in the deviatoric plane are shown in Figure 2.

**(a)**

**(b)**

#### 3. Case Study: Centrifuge Testing of Tunnelling

A series of centrifuge tests for tunnelling are performed on sand to validate and compare with the numerical simulations. In the centrifuge tests, two degrees of volume loss have been investigated, namely, and , respectively. As shown in Figure 3(a), the model dimension set up in the ABAQUS platform is the same as that in the centrifuge tests, with 25.6 m in length and 30.8 m in height. The diameter of the tunnel is 7.2 m. The tunnel is assumed to have been excavated in sand, 14.4 m below the ground surface. The cohesion is assumed to be *c* = 1.0 kPa in order to avoid the singularity problem. The maximum friction angle is *ϕ*_{max} = 30° [5, 31]. Typical elastic constants are Young’s modulus (*E* = 40 MPa) and Poisson’s ratio (*ν* = 0.25). The material surrounding the excavation is discretized with first-order 4-node plane strain elements (element type CPE4). On the right-side boundary, the infinite extent of the soil is represented by a 25.6 m wide mesh that extends offset from the center line to a length of 51.2 m. On the bottom boundary, the infinite extent of the soil is represented by an 18 m wide mesh that extends from the surface to a depth of 48.8 m below the surface. The left-hand boundary represents a vertical symmetry axis. Far-field conditions on the bottom and right-hand side boundaries are modelled by infinite elements (element type CINPE4). The relationship between the shape parameter *β* and the deposition direction can be found in Figure 3(b). In all cases, an initial stress regime with *K*_{0} = 0.5 and an associated flow rule are used.

**(a)**

**(b)**

Figure 4 shows the stress paths of the typical soil elements around the excavation face of the tunnel crow. It is seen that the principal stresses of the elements rotate obviously. This observation in turn demonstrates that it is necessary to take non-coaxiality into consideration when modelling tunnel excavations.

Figure 5 presents numerical results of the subsurface tunnelling-induced vertical displacements with various values of anisotropic coefficients *n* and *β* and non-coaxial coefficient *k*. The conventional isotropic Mohr–Coulomb yield criterion is retrieved when . The settlement trough obtained by using the non-coaxial soil model is identical with that obtained using the Mohr–Coulomb model embed in ABAQUS. This testifies the correctness of the finite element implementation procedures of the newly proposed model. As shown in Figure 5, the maximum vertical displacement () is larger with a deduction in the anisotropic coefficient *n* and an increase in the non-coaxial coefficient *k*. The maximum vertical displacement is not linearly related to the magnitude of the anisotropic coefficient *β*. However, a smaller value of *β* results in a narrower and steeper shape of subsurface settlement trough. Based on the analysis, it is concluded that the shape of tunnelling-induced subsurface settlement trough is unlikely influenced by the magnitude of *n* and *k*. Hence, the anisotropic coefficient is chosen as 0.707 for comparing with the centrifuge testing for reducing the repeat work of the parametric study.

**(a)**

**(b)**

The normalized settlement, which is defined as the vertical displacement () normalized by the maximum vertical displacement (), is plotted in Figure 6. Comparing Figure 6(a) with Figure 5(b), when the volume loss is low (), the normalized settlement troughs obtained by the new model and the Mohr–Coulomb model in ABAQUS are very close and coincide with the centrifuge test results and the Gaussian settlement trough. However, the normalized settlement trough becomes wider than the Gaussian settlement trough and centrifuge test results when using the Mohr–Coulomb model with a high volume loss (). Adopting the anisotropic coefficient *β* improves the settlement curve further. The combination of anisotropic coefficients and *β* = 0° leads to a much narrower normalized subsurface settlement trough (the red line with triangular marker). This curve is almost the same as the Gaussian curve and fits to the experimental data.

**(a)**

**(b)**

#### 4. Case Study: Field Testing of Tunnelling

Qinghuayuan Tunnel of Jing-Zhang High-Speed Railway is located in the core prosperous area of Beijing, which is believed to be the largest shield tunnel in Beijing till now. The tunnel is divided into two shield sections, of which the length of the north section is 1700 m. The excavation diameter of the tunnel is 12.64 m. Due to the dense distribution of the construction (structure) and underground pipelines around the tunnel line, the subsurface settlement has been monitored while tunnelling in order to avoid risks. The subsurface settlement of the monitoring section is taken 1000 m away from the launching shaft of the tunnel (section mileage: DK17 + 200) for comparative analysis, where the tunnel axis is buried 23.32 m below the ground surface. The stress reduction method is applied to achieve the same volume loss of . The tunnel is assumed to have been excavated in pebble soils, and material properties for the pebble soil are Young’s Modulus *E* = 136 MPa, Poisson’s ratio , the maximum internal friction angle *f*_{max} = 45°, and cohesion *c* = 0 kPa. The material properties of lining, with a thickness of 0.55 m, are Young’s Modulus *E* = 35.5 GPa and Poisson’s ratio .

Considering the actual tunnel size and Saint-Venant’s principle, a two-dimensional plane strain finite element field tunnel model is established in ABAQUS, and the size of which is 70 m × 80 m. The model is symmetric along the middle axis, so only half of the model is simulated, as shown in Figure 7. The material surrounding the excavation is discretized with first-order 4-node plane strain elements (element type CPE4), and the plane strain incompressible complete integral elements (element type CPE4I) are used to discretize the lining. On the left and right sides of the boundary, the horizontal constraint is applied, while the horizontal and vertical constraints are applied on the bottom boundary. The lining is connected with the soil of the tunnel excavation face by using the MPC tie. Then the lining and the contact will be activated simultaneously. In all cases, an initial stress regime with *K*_{0} = 0.5 and an associated flow rule are adopted.

Figures 8(a) and 9(a) present the subsurface settlement troughs with respect to various values of the shape parameters *n* and *β* and non-coaxial coefficient *k*, respectively. The normalized settlement trough with respect to various values of the shape parameters *n* and *β*, and the non-coaxial coefficient *k* are presented in Figures 8(b) and 9(b). As shown in Figures 8(a) and 9(b), the subsurface settlement troughs obtained from the field data are deeper than the numerical simulations; however, the normalized troughs obtained by using the proposed model are much closer to the field data compared to the Mohr–Coulomb model. A larger value of the non-coaxial coefficient results in a larger maximum vertical subsurface displacement.

**(a)**

**(b)**

**(a)**

**(b)**

Although the conclusions are similar with those of centrifuge testing results, the normalized subsurface settlement troughs are better improved in the case of centrifuge testing than that of the field testing. This may be because that the real excavation and support sequence are simplified in the numerical simulations. In addition, we assume only one layer of pebble soils; however, the site situation may involve layers of different soil types, and the effects of the soil lay are not considered in the present study.

#### 5. Conclusion

The influence of the non-coaxiality and soil strength anisotropy on the tunnel excavation problems has been investigated in this paper. The modified two-dimensional Mohr–Coulomb yield functions with soil yield anisotropy were employed. Then the proposed model was implemented into the finite element platform ABAQUS via the user-defined material subroutine. A case study in terms of centrifuge tests and field data were performed to analyze tunnel excavation problems. The following conclusions can be drawn:(i)In the centrifuge testing case when (i.e., the conventional isotropic Mohr–Coulomb yield criterion was retrieved), the settlement trough obtained by using the non-coaxial soil model was identical with that obtained by using the Mohr–Coulomb model embedded in ABAQUS. Stress rotations can be observed around the tunnel crowns, which demonstrated that it is necessary to take non-coaxiality into consideration in tunnlling excavations.(ii)In both centrifuge testing and field testing, the maximum vertical subsurface displacement was increasing with an increase in the anisotropic coefficient *n* and non-coaxial coefficient *k*. The shape of the normalized subsurface settlement trough was inclined to the centrifuge experimental results and field data when soil anisotropy was considered. The non-coaxial coefficient *k* exhibited negligible effects on the normalized subsurface settlement.(iii)The normalized subsurface settlement troughs were better improved in the case of centrifuge testing than those of field testing.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Key Research and Development Program of China (Grant no. 2016YFE0205200), the Fundamental Research Funds for the Central Universities (Grant No. 2682016CX084), and the Subject of Scientific Research and Development Plan of China Railway Corporation (Grant No. 2017G007-B). The author would like to thank the thesis written by Yuan [25] for providing the valuable data.