# Analytical Study of Soil Displacement Induced by Twin Shield Tunneling in Semi-Infinite Viscoelastic Ground

Authors: Xin Shi; Chuanxin Rong; Houliang Wang; Linzhao Cui; Haibing Cai; Bin Wang
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Effective measures are needed to strictly control soil displacement caused during the process of shield construction excavation for urban subway tunnels. When calculating the displacement of soil caused by loading or unloading, many previous analytical studies have assumed that the soil was a linear elastic body and ignored the viscosity of the soil. In this study, the Boltzmann viscoelastic model and the Mindlin basic solution were combined to consider the effects of the additional support pressures, the shield shell frictions, the grouting pressures, and the ground loss, and a three-dimensional viscoelastic solution for soil displacement caused by shield tunneling was derived. According to the calculation results of an example, the analytical solution was able to consider the asynchronous construction of the left and right tunnels and the mutual influence of the double shield tunnel. The rationality of the approach proposed in this study was verified by comparing the theoretical solution with the measured settlement values. In addition, the influence of differences in the viscoelastic parameters (the viscosity coefficient, the shear modulus of the elastic element, and the shear modulus of the viscous element) and the geometric parameters (the distance from the excavation surface, the calculated depth, and tunnel spacing) on soil displacement is discussed. The calculation method in this study provides a theoretical basis for predicting the three-dimensional soil deformation caused by shield tunneling, especially in soft clays.


---

## Body

## 1. Introduction

Due to the continuous deterioration of surface traffic conditions, many large and medium-sized cities have built underground transportation networks to ease pressure on surface traffic. Tunnel construction is the first consideration when establishing an underground traffic system. In densely populated urban areas, soil displacement and soil-structure interaction caused by shield tunneling are issues that must be considered. Excavation will inevitably affect the displacement field and stress field of the soil. Therefore, the ground displacement pattern caused by construction disturbance has always been the focus of discussion.In terms of theoretical analysis, Zhang et al. [1] proposed a new analytical method to efficiently assess the ground surface settlements induced by twin tunneling in clay, where the nonuniform contraction deformation pattern at the tunnel opening boundary was intensively examined. Lu et al. [2] analyzed the effects of the fractional order and viscous material parameters on the relaxation modulus, and the correspondence principle and complex variable method are adopted to develop the elastic analytical solution into the fractional viscoelastic analytical solution of the ground displacement. Zhang et al. [3] presented a closed-form analytical solution for predicting long- and short-term ground deformations and liner internal forces induced by tunneling in saturated soils in which shield excavation effects with and without air pressure are both considered. Jin et al. [4] proposed an approach to estimate three-dimensional ground displacements induced by shield tunneling using the superposition method, and they also considered the effects of ground loss, additional support pressures, shield shell frictions, and cutter head rotation during shield tunneling. Hesham and Dipanjan [5] established the governing differential equations for beam and soil displacements using the extended Hamilton’s principle, and they proposed a new method for the dynamic analysis of the Euler–Bernoulli beams resting on multilayered viscoelastic soils.In terms of numerical simulations and field measurements, Zhang et al. [6] described the key inﬂuences of yaw excavation loadings on ground displacement and segmental stress for a curved shield tunnel. The inﬂuences are investigated through finite element models, the reliabilities of which are validated through comparisons to field data and analytical solutions. Lin et al. [7] established a series of simplified three-dimensional FEM models on the background of Changsha Metro Line 2 in China and investigated the stress redistribution and the soil arching evolution induced by the earth pressure balance shield (EPBS) tunneling. In addition, many scholars [8–10] have further studied the structural response caused by tunnel excavation based on the soil displacement caused by tunnel excavation, including the impact on adjacent piles [11–14], existing tunnels [15–19], and underground pipelines [20–23].In terms of experimental studies, Song and Marshall [24] believed that the choice of the model tunnel has an effect on the imparted tunnel boundary displacements and resulting ground deformations, and contrasting plane-strain centrifuge test results from experiments using a flexible membrane model tunnel with those from a newly developed eccentric rigid boundary mechanical model tunnel to quantitatively evaluate this effect. Loganathan et al. [25] studied the deformation of the clay foundation caused by tunnel excavation and its effect on the adjacent pile foundation using three centrifugal model tests. Long and Tan [26] investigated the geo-hazards and explored the associated failure mechanism due to tunnel leaking, and conducted experimental studies and extensive numerical simulations using the validated finite-difference-method and discrete-element-method (FDM-DEM) coupling method.When calculating the displacement of the soil caused by loading or unloading, it has primarily been assumed that the soil was a linear elastic body. This is not consistent with the characteristics of soil, especially in soft clays, and this assumption may even produce erroneous calculation results. Some scholars had deduced a two-dimensional viscoelastic solution that considers the effect of soil viscoelasticity. However, tunnel excavation is a three-dimensional dynamic process, and the disturbance caused by shield tunneling and the displacement of soil caused by additional forces change with time. Therefore, the research regarding the three-dimensional viscoelastic analytical solution of shield tunneling is of great significance to accurately understand and predict the soil displacement and soil-structure interaction caused by construction.This study combines the Boltzmann viscoelastic model and the Mindlin basic solution and deduces a three-dimensional viscoelastic solution of soil displacement caused by shield tunneling that considers the effects of the additional support pressures, the shield shell frictions, the grouting pressures, and the ground loss. The rationality of the solution is then verified using calculation examples and on-site measurements. In addition, the influence of differences in the viscoelastic parameters and geometric parameters on soil displacement is discussed. The calculation method proposed in this study provides a theoretical basis for predicting the three-dimensional soil deformation caused by shield tunneling.


## 2. Basic Assumptions and Theories

### 2.1. Basic Assumptions

This study makes the following assumptions about the soil: (1) Soil is a linear viscoelastic material. (2) It is a uniform and isotropic continuous deformation body that extends infinitely in the depth direction and is regarded as a semi-infinite body. (3) It is in a three-dimensional stress state under the action of internal forces, and the relationship between the deviatoric stress tensor and the deviatoric strain tensor is a viscoelastic stress-strain relationship. (4) The additional support pressures, the shield shell frictions, and the grouting pressures produced by twin shield tunneling are uniform distributed.


### 2.2. Viscoelastic Model

It is assumed that the viscoelastic stress-strain relationship of the soil conforms to the Boltzmann viscoelastic model, as shown in Figure1. This model is represented by a serial arrangement of the Kevin–Voigt model and an elastic relation. It introduces an elastic element, H, and a viscous element, K, connected in series. It can reflect the instantaneous elastic deformation of the soil under load and the viscoelastic deformation that gradually increases with time and approaches stability.Figure 1
The Boltzmann viscoelastic model.The constitutive relationship of the model is as follows:(1)G1+Gkσ+ηkdσdt=G1Gkε+G1ηkdεdt,where G1, Gk, and ηk are the shear modulus of the elastic element, shear modulus of the viscous element, and viscosity coefficient of the Boltzmann viscoelastic model, respectively.The rheological properties of the soil are represented by the viscoelastic relationship between the deviatoric stress tensor and the deviatoric strain tensor, and the volume change in the soil is represented by the elastic relationship between the spherical stress tensor and the spherical strain tensor. Therefore, from equation (1), the viscoelastic relationship of the deviatoric tensor is as follows:(2)G1+Gk⋅Sijt+ηkdSijdt=2G1Gk⋅eijt+2G1ηkdeijdt,where Sijt is the deviatoric stress tensor, and eijt is the deviatoric strain tensor.The elastic relationship of the spherical tensor is as follows:(3)σkkt=3K1εkkt,where σkkt is the spherical stress tensor; εkkt is the spherical strain tensor; K1 is the bulk modulus.According to the classical three-dimensional viscoelastic constitutive relationship, the unified expression is as follows:(4)P′Sijt=Q′eijt,P″σkkt=Q″εkkt,where P′, Q′, P″, and Q″ are linear differential operators of time variables, which can be expressed as follows:(5)P′=∑k=0mpk′dkdkt,Q′=∑k=0mqk′dkdkt,P″=∑k=0mpk″dkdkt,Q″=∑k=0mqk″dkdkt.From equations (2) and (3), the linear differential operators in equation (4) are as follows:(6)P′t=G1+Gk+ηkddt,Q′t=2G1Gk+2G1ηkddt,P″t=1,Q″t=3K1,The Laplace transform of equation (6) results in the following:(7)P′¯t=G1+Gk+ηks,Q′¯t=2G1Gk+2G1ηks,P″¯t=1,Q″¯t=3K1.According to the elastic-viscoelastic correspondence principle, the Laplace transform expressions of the elastic modulus, Poisson’s ratio, and the shear modulus are, respectively,(8)E¯s=3Q″¯sQ′¯s2Q″¯sP′¯s+Q′¯sP″¯s,v¯s=Q″¯sP′¯s−Q′¯sP″¯s2Q″¯sP′¯s+Q′¯sP″¯s,G¯s=E¯s21+v¯s=12⋅Q′¯sP′¯s=G1Gk+G1ηksG1+Gk+ηks.


### 2.3. The Mindlin Basic Solution of Semi-Infinite Space

Wei et al. [27] showed that the primary factors that caused surface deformation during shield tunneling were the additional support pressures, the shield shell frictions, grouting pressures, and ground loss. In addition, the calculation formula for the three-dimensional ground deformation caused by various factors was derived. Based on the Mindlin solution of elasticity, when a unit concentrated force is applied along the x-axis at an internal point (0, 0, h) of a semi-infinite space, the displacements along the x-, y-, and z-axes can be expressed as follows:(9)uxx,y,h=116πG1−v3−4vM+1N+x2M3+3−4vN3+2zhN31−3x2N2+41−v1−2vN+z+h1−x2NN+z+h,(10)vxx,y,h=xy16πG1−v1M3+3−4vN3−6zhN5−41−v1−2vNN+z+h2,(11)wxx,y,h=x16πG1−vz−hM3+3−4vz−hN3−6zhz+hN5+41−v1−2vNN+z+h,where M=x2+y2+z−h2, N=x2+y2+z+h2, G is the shear modulus of the soil, and v is Poisson’s ratio.


## 2.1. Basic Assumptions

This study makes the following assumptions about the soil: (1) Soil is a linear viscoelastic material. (2) It is a uniform and isotropic continuous deformation body that extends infinitely in the depth direction and is regarded as a semi-infinite body. (3) It is in a three-dimensional stress state under the action of internal forces, and the relationship between the deviatoric stress tensor and the deviatoric strain tensor is a viscoelastic stress-strain relationship. (4) The additional support pressures, the shield shell frictions, and the grouting pressures produced by twin shield tunneling are uniform distributed.


## 2.2. Viscoelastic Model

It is assumed that the viscoelastic stress-strain relationship of the soil conforms to the Boltzmann viscoelastic model, as shown in Figure1. This model is represented by a serial arrangement of the Kevin–Voigt model and an elastic relation. It introduces an elastic element, H, and a viscous element, K, connected in series. It can reflect the instantaneous elastic deformation of the soil under load and the viscoelastic deformation that gradually increases with time and approaches stability.Figure 1
The Boltzmann viscoelastic model.The constitutive relationship of the model is as follows:(1)G1+Gkσ+ηkdσdt=G1Gkε+G1ηkdεdt,where G1, Gk, and ηk are the shear modulus of the elastic element, shear modulus of the viscous element, and viscosity coefficient of the Boltzmann viscoelastic model, respectively.The rheological properties of the soil are represented by the viscoelastic relationship between the deviatoric stress tensor and the deviatoric strain tensor, and the volume change in the soil is represented by the elastic relationship between the spherical stress tensor and the spherical strain tensor. Therefore, from equation (1), the viscoelastic relationship of the deviatoric tensor is as follows:(2)G1+Gk⋅Sijt+ηkdSijdt=2G1Gk⋅eijt+2G1ηkdeijdt,where Sijt is the deviatoric stress tensor, and eijt is the deviatoric strain tensor.The elastic relationship of the spherical tensor is as follows:(3)σkkt=3K1εkkt,where σkkt is the spherical stress tensor; εkkt is the spherical strain tensor; K1 is the bulk modulus.According to the classical three-dimensional viscoelastic constitutive relationship, the unified expression is as follows:(4)P′Sijt=Q′eijt,P″σkkt=Q″εkkt,where P′, Q′, P″, and Q″ are linear differential operators of time variables, which can be expressed as follows:(5)P′=∑k=0mpk′dkdkt,Q′=∑k=0mqk′dkdkt,P″=∑k=0mpk″dkdkt,Q″=∑k=0mqk″dkdkt.From equations (2) and (3), the linear differential operators in equation (4) are as follows:(6)P′t=G1+Gk+ηkddt,Q′t=2G1Gk+2G1ηkddt,P″t=1,Q″t=3K1,The Laplace transform of equation (6) results in the following:(7)P′¯t=G1+Gk+ηks,Q′¯t=2G1Gk+2G1ηks,P″¯t=1,Q″¯t=3K1.According to the elastic-viscoelastic correspondence principle, the Laplace transform expressions of the elastic modulus, Poisson’s ratio, and the shear modulus are, respectively,(8)E¯s=3Q″¯sQ′¯s2Q″¯sP′¯s+Q′¯sP″¯s,v¯s=Q″¯sP′¯s−Q′¯sP″¯s2Q″¯sP′¯s+Q′¯sP″¯s,G¯s=E¯s21+v¯s=12⋅Q′¯sP′¯s=G1Gk+G1ηksG1+Gk+ηks.


## 2.3. The Mindlin Basic Solution of Semi-Infinite Space

Wei et al. [27] showed that the primary factors that caused surface deformation during shield tunneling were the additional support pressures, the shield shell frictions, grouting pressures, and ground loss. In addition, the calculation formula for the three-dimensional ground deformation caused by various factors was derived. Based on the Mindlin solution of elasticity, when a unit concentrated force is applied along the x-axis at an internal point (0, 0, h) of a semi-infinite space, the displacements along the x-, y-, and z-axes can be expressed as follows:(9)uxx,y,h=116πG1−v3−4vM+1N+x2M3+3−4vN3+2zhN31−3x2N2+41−v1−2vN+z+h1−x2NN+z+h,(10)vxx,y,h=xy16πG1−v1M3+3−4vN3−6zhN5−41−v1−2vNN+z+h2,(11)wxx,y,h=x16πG1−vz−hM3+3−4vz−hN3−6zhz+hN5+41−v1−2vNN+z+h,where M=x2+y2+z−h2, N=x2+y2+z+h2, G is the shear modulus of the soil, and v is Poisson’s ratio.


## 3. Viscoelastic Solution of Three-Dimensional Soil Displacement Caused by Shield Tunneling

### 3.1. Mechanical Calculation Model

The establishment of the coordinate system and the mechanical calculation model are shown in Figure2. It is assumed that the shield on the right is constructed first. The parameters involved in the calculation included the distance between the first and second excavation surface, K, the distance between the tunnel centerline, L, the buried depth of the tunnel, h, and the outer radius of the tunnel, R, the length of the shield machine, J, the additional support pressures on the excavation surface, P1, the shield shell frictions, P2, and the grouting pressure, P3.Figure 2
Mechanical calculation model of soil deformation caused by shield tunneling.


### 3.2. Soil Displacements Induced by Additional Support Pressures

The additional support pressure,P1, is applied to the excavation surface of the tunnel, and the differential area, rdrdθ, is within the excavation surface. Hence, the support pressures acting on this area is dp1=P1rdrdθ. The coordinates of any point in the left and right excavation surfaces are 0,−L/2+r⋅cosθ,h−r⋅sinθ and K,L/2+r⋅cosθ,h−r⋅sinθ, respectively, and the equivalent coordinates that can be brought into the Mindlin solution obtained by the coordinate transformation are (subscript l stands for the left tunnel, r stands for the right tunnel) as follows:(12)x1l=x,y1l=y+L2−r⋅cosθ,h1l=h−r⋅sinθ,(13)x1l=x−K,y1l=y+L2−r⋅cosθ,h1l=h−r⋅sinθ.By incorporating equations (12) and (13) into equations (10) and (11), and by integrating within the range of the excavation surface, the displacement of any point (x, y, z) in the Y and Z directions caused by the additional support pressures can be obtained as follows:(14)v1=v1l+v1r,w1=w1l+w1r,where for left tunnel, v1l=∫02π∫0RP1l⋅vxx1l,y1l,h1l⋅rdrdθ and w1l=∫02π∫0RP1l⋅wxx1l,y1l,h1l⋅rdrdθ, and for right tunnel, v1r=∫02π∫0RP1r⋅vxx1r,y1r,h1r⋅rdrdθ, and w1r=∫02π∫0RP1r⋅wxx1r,y1r,h1r⋅rdrdθ.According to the elastic-viscoelastic correspondence principle, under the condition of the same force on the viscoelastic body, the Laplace transformation ofv1 and w1 with respect to time, can be obtained as follows:(15)v1¯s=v1l¯s+v1r¯s=∫02π∫0RP1ls⋅wx¯x1l,y1l,h1l⋅rdrdθ+∫02π∫0RP1rs⋅vx¯x1r,y1r,h1r⋅rdrdθ,(16)w1¯s=w1l¯s+w1r¯s=∫02π∫0RP1ls⋅wx¯x1l,y1l,h1l⋅rdrdθ+∫02π∫0RP1rs⋅wx¯x1r,y1r,h1r⋅rdrdθ.The Laplace inverse transformation of equations (15) and (16) with respect to time, t, are then performed, and the viscoelastic solutions of the displacement components along the Y and Z directions generated at any point (x, y, z) caused by the additional support pressures in the semi-infinite space are v1¯x,y,z,t and w1¯x,y,z,t, respectively. The detailed formula derivation results are provided in the Appendix. v1¯x,y,z,t and w1¯x,y,z,t can be expressed as follows:(17)v1¯x,y,z,t=v1l¯x,y,z,t+v1r¯x,y,z,t,(18)w1¯x,y,z,t=w1l¯x,y,z,t+w1r¯x,y,z,t.


### 3.3. Soil Displacements Induced by the Shield Shell Frictions

For the frictional force,P2, between the shield shell and the soil, take the differential area of the shield machine surface, Rdjdθ, and the frictional force on the differential area is dp2=P2Rdjdθ. The length of the shield machine is J. The coordinate of any point on the interface between the shield shell on the left and the soil is −j,−L/2+R⋅cosθ,h−R⋅sinθ. In the same way, the right side is −j+K,L/2+R⋅cosθ,h−R⋅sinθ. The equivalent coordinates that can be brought into the Mindlin solution and after coordinate transformations are as follows:(19)x2l=x+j,y2l=y+L2−R⋅cosθ,h2l=h−R⋅sinθ,(20)x2r=x+j−K,y2r=y−L2−R⋅cosθ,h2r=h−R⋅sinθ.By incorporating equations (19) and (20) into equations (10) and (11), and by integrating within the range of the shield machine length J, the displacement of any point (x, y, z) in the Y and Z directions caused by the shield shell frictions can be obtained as follows:(21)v2=v2l+v2r,w2=w2l+w2r,where for the left tunnel, v2l=∫02π∫0JP2l⋅vxx2l,y2l,h2l⋅Rdjdθ and w2l=∫02π∫0JP2l⋅wxx2l,y2l,h2l⋅Rdjdθ, and for the right tunnel, v2r=∫02π∫0JP2r⋅vxx2r,y2r,h2r⋅Rdjdθ, and w2r=∫02π∫0JP2r⋅wxx2r,y2r,h2r⋅Rdjdθ.According to the elastic-viscoelastic correspondence principle, under the condition of the same force on the viscoelastic body, the Laplace transformation ofv2 and w2 with respect to time, t, can be obtained as follows:(22)v2¯s=v2l¯s+v2r¯s=∫02π∫0JP2ls⋅vx¯x2l,y2l,h2l⋅Rdjdθ+∫02π∫0JP2rs⋅vx¯x2r,y2r,h2r⋅Rdjdθ,(23)w2¯s=w2l¯s+w2r¯s=∫02π∫0JP2ls⋅wx¯x2l,y2l,h2l⋅Rdjdθ+∫02π∫0JP2rs⋅wx¯x2r,y2r,h2r⋅Rdjdθ.The Laplace inverse transformation of equations (22) and (23) with respect to time, t, is then performed, and the viscoelastic solutions of the displacement components along the Y and Z directions generated at any point (x, y, z) caused by the shield shell frictions in the semi-infinite space are v2¯x,y,z,t and w2¯x,y,z,t, respectively. The detailed formula derivation results are provided in the Appendix. It can be expressed as follows:(24)v2¯x,y,z,t=v2l¯x,y,z,t+v2r¯x,y,z,t,(25)w2¯x,y,z,t=w2l¯x,y,z,t+w2r¯x,y,z,t.


### 3.4. Soil Displacements Induced by Grouting Pressures

For the grouting pressure,P3, the derivation process of its viscoelastic solution was similar to the additional support pressures and the shield shell frictions. Taking the differential area of the shield tail, Rdjdθ, the concentrated force in this area was dp3=P3Rdjdθ, and decompose dp3 into horizontal force dp3h=−P3Rcosθdθdj and vertical force dp3v=−P3Rsinθdθdj.By rotating the coordinate system, it can be obtained that the displacement of the soil in they-axis direction when the unit concentrated force is applied along the y-axis at a point (0, 0, h) inside the semi-infinite space is as follows:(26)vYx,y,h=116πG1−μ3−4μM+1N+y2M3+3−4μy2N3+2zhN31−3y2N2+41−μ1−2μN+z+h1−y2NN+z+h.The displacement of the soil in thez-axis direction is as follows:(27)wYx,y,h=y16πG1−μz−hM3+3−4μz−hN3−6zhz+hN5+41−μ1−2μNN+z+h.In the same way, when the unit concentrated force along thez-axis is applied at a point (0, 0, h) inside the semi-infinite space, the displacement of the soil in the y-axis direction is as follows:(28)vZx,y,h=y16πG1−μz−hM3+3−4μz−hN3+6zhz+hN5−41−μ1−2μNN+z+h.The displacement of the soil in thez-axis direction is as follows:(29)wZx,y,h=116πG1−μ3−4μM+81−μ2−3−4μN+z−h2M3+3−4μz+h2−2hzN3+6zhz+h2N5.The coordinate of any point on the area of the grouting pressure on the left is−J−j,−L/2+R⋅cosθ,h−R⋅sinθ, and the right side is −J−j+K,L/2+R⋅cosθ,h−R⋅sinθ. The equivalent coordinates that can be brought into the Mindlin solution and after coordinate transformations are as follows:(30)x3l=x+J+j,y3l=y+L2−R⋅cosθ,h3l=h−R⋅sinθ,(31)x3r=x+J+j−K,y3r=y−L2−R⋅cosθ,h3r=h−R⋅sinθ.Then the solution method of the grouting pressures is referred to, and equations (30) and (31) are brought into equations (26), (27), (28), and (29), and then the Laplace transform and the inverse transform are performed. Therefore, the time-domain solutions of the displacement components along the Y and Z directions generated at any point (x, y, z) in the semi-infinite viscoelastic space caused by the grouting pressure can be expressed as follows:(32)v3¯x,y,z,t=v3l¯x,y,z,t+v3r¯x,y,z,t,(33)w3¯x,y,z,t=w3l¯x,y,z,t+w3r¯x,y,z,t.


### 3.5. Soil Displacements Induced by Ground Loss

For soil displacement caused by ground loss, since no viscoelastic parameters are involved, this study used the uniform ground movement model proposed by Wei [28] to consider different soil conditions for calculation. By revising the calculation formulas of Verruijt and Booker [29], a two-dimensional solution of soil deformation caused by ground loss during shield tunnel construction was derived [30]. It was assumed that the ground loss mainly produced vertical soil displacement along the excavation direction while ignoring the horizontal displacement. On the basis of the two-dimensional solution, Wei [31] defined the ground loss rate η=Vloss/πR2, Vloss as the amount of ground loss per unit length and proposed the expression of the ground loss caused by shield excavation along the X direction as follows:(34)ηx=η21−xx2+h2.From equation (34), it can be concluded that when x⟶+∞, ηx=0; when x⟶0, ηx=η/2; when x⟶−∞, ηx=η/2. At the same time, when x=−3h, ηx=0.974η. This indicates that the ground loss basically reaches the maximum value beyond three times the buried depth of the tunnel behind the shield excavation surface.Similarly, assuming that the tunnel on the right is excavated first, combined with the two-dimensional solution and equation (34), the horizontal and vertical displacement of any point, (x, y, z), in the soil caused by ground loss can be expressed as follows:(35)v4x,y,z=v4lx,y,z+v4rx,y,z=−R2y+L/2+b2⋅hh+dl1y+L/2+b2+h−z2+1y+L/2+b2+h+z24zh+zy+L/2+b2+h−z22⋅4Rg1x−gl2x4R2Blx⋅expy+L/2+b2Inλlxh+R2+z2Inλlx−Inδlxh+dl2+−R2y−L/22⋅hh+dr⋅1y−L/22+h−z2+1y−L/22+h+z2−4zh+zy−L/22+h+z22⋅4Rgrx−gr2x4R2Brx⋅expy−L/22lnλrxh+R2+z2lnλrx−lnδrxh+dr2,(36)w4x,y,z=w4lx,y,z+w4rx,y,z=R22h−zy+L/2+b2+h−z2+h+zy+L/2+b2+h+z2−2zy+L/2+b2−h+z2y+L/2+b2+h+z22⋅4Rglx−gl2x4R2Blx⋅expy+L/2+b2lnλlxh+R2+z2lnλlx−lnδlxh+dl2+R22h−zy−L/22+h−z2+h+zy−L/22+h+z2−2zy−L/22−h+z2y−L/22+h+z22⋅4Rgrx−gr2x4R2Brx⋅expy−L/22lnλrxh+R2+z2lnλrx−lnδrxh+dr2,where gx=2R1−1−ηx, Bx=4hh+d−h+d2−ηxR+d2/RηxR+d,(37)λx=14−gxπRηxarcsindR−gx/2+1−dR−gx/22−1,δx=12−gxπR2ηxR−gx4arcsindR−gx/4,where gx and ηx are the equivalent ground loss parameter and ground loss rate, respectively, at a distance of x from the excavation surface along the tunneling direction, and Bx, λx, and δx are the calculated parameters.The three-dimensional solution introduces the offset,b, in the horizontal direction, and at the same time differentiates the calculation parameters of the first and the second shield construction to consider the mutual influence of the first and the second tunnel construction. Refer to references [32, 33] for the values of the parameters η and d. Affected by the prior construction of the right tunnel, the ground loss rate caused by the second shield construction in this study can be expressed as follows [27]:(38)η1=0.0017h+DL−0.1454ηr.


### 3.6. The Three-Dimensional Viscoelastic Solution of the Total Soil Displacement

The numerical calculation program was written in MATLAB R2017. The horizontal displacement of the soil in theY direction was obtained by superposing equations (17), (24), (32), and (35), and the vertical displacement of the soil in the Z direction was obtained by superposing equations (18), (25), (33), and (36). By combining equation (9) with the calculation method in this study, the horizontal displacement in the X direction can also be derived. However, it is generally ignored in actual construction, so this study only discusses the changes in the soil displacement values in the Y and Z directions under different X coordinates. Then the time-domain solutions of the horizontal and vertical displacement of any point, (x, y, z), in the ground caused by the additional support pressures, shield shell frictions, grouting pressures, and ground loss can be expressed as follows:(39)vx,y,z,t=v1¯x,y,z,t+v2¯x,y,z,t+v3¯x,y,z,t+v4x,y,z,wx,y,z,t=w1¯x,y,z,t+w2¯x,y,z,t+w3¯x,y,z,t+w4x,y,z.


## 3.1. Mechanical Calculation Model

The establishment of the coordinate system and the mechanical calculation model are shown in Figure2. It is assumed that the shield on the right is constructed first. The parameters involved in the calculation included the distance between the first and second excavation surface, K, the distance between the tunnel centerline, L, the buried depth of the tunnel, h, and the outer radius of the tunnel, R, the length of the shield machine, J, the additional support pressures on the excavation surface, P1, the shield shell frictions, P2, and the grouting pressure, P3.Figure 2
Mechanical calculation model of soil deformation caused by shield tunneling.


## 3.2. Soil Displacements Induced by Additional Support Pressures

The additional support pressure,P1, is applied to the excavation surface of the tunnel, and the differential area, rdrdθ, is within the excavation surface. Hence, the support pressures acting on this area is dp1=P1rdrdθ. The coordinates of any point in the left and right excavation surfaces are 0,−L/2+r⋅cosθ,h−r⋅sinθ and K,L/2+r⋅cosθ,h−r⋅sinθ, respectively, and the equivalent coordinates that can be brought into the Mindlin solution obtained by the coordinate transformation are (subscript l stands for the left tunnel, r stands for the right tunnel) as follows:(12)x1l=x,y1l=y+L2−r⋅cosθ,h1l=h−r⋅sinθ,(13)x1l=x−K,y1l=y+L2−r⋅cosθ,h1l=h−r⋅sinθ.By incorporating equations (12) and (13) into equations (10) and (11), and by integrating within the range of the excavation surface, the displacement of any point (x, y, z) in the Y and Z directions caused by the additional support pressures can be obtained as follows:(14)v1=v1l+v1r,w1=w1l+w1r,where for left tunnel, v1l=∫02π∫0RP1l⋅vxx1l,y1l,h1l⋅rdrdθ and w1l=∫02π∫0RP1l⋅wxx1l,y1l,h1l⋅rdrdθ, and for right tunnel, v1r=∫02π∫0RP1r⋅vxx1r,y1r,h1r⋅rdrdθ, and w1r=∫02π∫0RP1r⋅wxx1r,y1r,h1r⋅rdrdθ.According to the elastic-viscoelastic correspondence principle, under the condition of the same force on the viscoelastic body, the Laplace transformation ofv1 and w1 with respect to time, can be obtained as follows:(15)v1¯s=v1l¯s+v1r¯s=∫02π∫0RP1ls⋅wx¯x1l,y1l,h1l⋅rdrdθ+∫02π∫0RP1rs⋅vx¯x1r,y1r,h1r⋅rdrdθ,(16)w1¯s=w1l¯s+w1r¯s=∫02π∫0RP1ls⋅wx¯x1l,y1l,h1l⋅rdrdθ+∫02π∫0RP1rs⋅wx¯x1r,y1r,h1r⋅rdrdθ.The Laplace inverse transformation of equations (15) and (16) with respect to time, t, are then performed, and the viscoelastic solutions of the displacement components along the Y and Z directions generated at any point (x, y, z) caused by the additional support pressures in the semi-infinite space are v1¯x,y,z,t and w1¯x,y,z,t, respectively. The detailed formula derivation results are provided in the Appendix. v1¯x,y,z,t and w1¯x,y,z,t can be expressed as follows:(17)v1¯x,y,z,t=v1l¯x,y,z,t+v1r¯x,y,z,t,(18)w1¯x,y,z,t=w1l¯x,y,z,t+w1r¯x,y,z,t.


## 3.3. Soil Displacements Induced by the Shield Shell Frictions

For the frictional force,P2, between the shield shell and the soil, take the differential area of the shield machine surface, Rdjdθ, and the frictional force on the differential area is dp2=P2Rdjdθ. The length of the shield machine is J. The coordinate of any point on the interface between the shield shell on the left and the soil is −j,−L/2+R⋅cosθ,h−R⋅sinθ. In the same way, the right side is −j+K,L/2+R⋅cosθ,h−R⋅sinθ. The equivalent coordinates that can be brought into the Mindlin solution and after coordinate transformations are as follows:(19)x2l=x+j,y2l=y+L2−R⋅cosθ,h2l=h−R⋅sinθ,(20)x2r=x+j−K,y2r=y−L2−R⋅cosθ,h2r=h−R⋅sinθ.By incorporating equations (19) and (20) into equations (10) and (11), and by integrating within the range of the shield machine length J, the displacement of any point (x, y, z) in the Y and Z directions caused by the shield shell frictions can be obtained as follows:(21)v2=v2l+v2r,w2=w2l+w2r,where for the left tunnel, v2l=∫02π∫0JP2l⋅vxx2l,y2l,h2l⋅Rdjdθ and w2l=∫02π∫0JP2l⋅wxx2l,y2l,h2l⋅Rdjdθ, and for the right tunnel, v2r=∫02π∫0JP2r⋅vxx2r,y2r,h2r⋅Rdjdθ, and w2r=∫02π∫0JP2r⋅wxx2r,y2r,h2r⋅Rdjdθ.According to the elastic-viscoelastic correspondence principle, under the condition of the same force on the viscoelastic body, the Laplace transformation ofv2 and w2 with respect to time, t, can be obtained as follows:(22)v2¯s=v2l¯s+v2r¯s=∫02π∫0JP2ls⋅vx¯x2l,y2l,h2l⋅Rdjdθ+∫02π∫0JP2rs⋅vx¯x2r,y2r,h2r⋅Rdjdθ,(23)w2¯s=w2l¯s+w2r¯s=∫02π∫0JP2ls⋅wx¯x2l,y2l,h2l⋅Rdjdθ+∫02π∫0JP2rs⋅wx¯x2r,y2r,h2r⋅Rdjdθ.The Laplace inverse transformation of equations (22) and (23) with respect to time, t, is then performed, and the viscoelastic solutions of the displacement components along the Y and Z directions generated at any point (x, y, z) caused by the shield shell frictions in the semi-infinite space are v2¯x,y,z,t and w2¯x,y,z,t, respectively. The detailed formula derivation results are provided in the Appendix. It can be expressed as follows:(24)v2¯x,y,z,t=v2l¯x,y,z,t+v2r¯x,y,z,t,(25)w2¯x,y,z,t=w2l¯x,y,z,t+w2r¯x,y,z,t.


## 3.4. Soil Displacements Induced by Grouting Pressures

For the grouting pressure,P3, the derivation process of its viscoelastic solution was similar to the additional support pressures and the shield shell frictions. Taking the differential area of the shield tail, Rdjdθ, the concentrated force in this area was dp3=P3Rdjdθ, and decompose dp3 into horizontal force dp3h=−P3Rcosθdθdj and vertical force dp3v=−P3Rsinθdθdj.By rotating the coordinate system, it can be obtained that the displacement of the soil in they-axis direction when the unit concentrated force is applied along the y-axis at a point (0, 0, h) inside the semi-infinite space is as follows:(26)vYx,y,h=116πG1−μ3−4μM+1N+y2M3+3−4μy2N3+2zhN31−3y2N2+41−μ1−2μN+z+h1−y2NN+z+h.The displacement of the soil in thez-axis direction is as follows:(27)wYx,y,h=y16πG1−μz−hM3+3−4μz−hN3−6zhz+hN5+41−μ1−2μNN+z+h.In the same way, when the unit concentrated force along thez-axis is applied at a point (0, 0, h) inside the semi-infinite space, the displacement of the soil in the y-axis direction is as follows:(28)vZx,y,h=y16πG1−μz−hM3+3−4μz−hN3+6zhz+hN5−41−μ1−2μNN+z+h.The displacement of the soil in thez-axis direction is as follows:(29)wZx,y,h=116πG1−μ3−4μM+81−μ2−3−4μN+z−h2M3+3−4μz+h2−2hzN3+6zhz+h2N5.The coordinate of any point on the area of the grouting pressure on the left is−J−j,−L/2+R⋅cosθ,h−R⋅sinθ, and the right side is −J−j+K,L/2+R⋅cosθ,h−R⋅sinθ. The equivalent coordinates that can be brought into the Mindlin solution and after coordinate transformations are as follows:(30)x3l=x+J+j,y3l=y+L2−R⋅cosθ,h3l=h−R⋅sinθ,(31)x3r=x+J+j−K,y3r=y−L2−R⋅cosθ,h3r=h−R⋅sinθ.Then the solution method of the grouting pressures is referred to, and equations (30) and (31) are brought into equations (26), (27), (28), and (29), and then the Laplace transform and the inverse transform are performed. Therefore, the time-domain solutions of the displacement components along the Y and Z directions generated at any point (x, y, z) in the semi-infinite viscoelastic space caused by the grouting pressure can be expressed as follows:(32)v3¯x,y,z,t=v3l¯x,y,z,t+v3r¯x,y,z,t,(33)w3¯x,y,z,t=w3l¯x,y,z,t+w3r¯x,y,z,t.


## 3.5. Soil Displacements Induced by Ground Loss

For soil displacement caused by ground loss, since no viscoelastic parameters are involved, this study used the uniform ground movement model proposed by Wei [28] to consider different soil conditions for calculation. By revising the calculation formulas of Verruijt and Booker [29], a two-dimensional solution of soil deformation caused by ground loss during shield tunnel construction was derived [30]. It was assumed that the ground loss mainly produced vertical soil displacement along the excavation direction while ignoring the horizontal displacement. On the basis of the two-dimensional solution, Wei [31] defined the ground loss rate η=Vloss/πR2, Vloss as the amount of ground loss per unit length and proposed the expression of the ground loss caused by shield excavation along the X direction as follows:(34)ηx=η21−xx2+h2.From equation (34), it can be concluded that when x⟶+∞, ηx=0; when x⟶0, ηx=η/2; when x⟶−∞, ηx=η/2. At the same time, when x=−3h, ηx=0.974η. This indicates that the ground loss basically reaches the maximum value beyond three times the buried depth of the tunnel behind the shield excavation surface.Similarly, assuming that the tunnel on the right is excavated first, combined with the two-dimensional solution and equation (34), the horizontal and vertical displacement of any point, (x, y, z), in the soil caused by ground loss can be expressed as follows:(35)v4x,y,z=v4lx,y,z+v4rx,y,z=−R2y+L/2+b2⋅hh+dl1y+L/2+b2+h−z2+1y+L/2+b2+h+z24zh+zy+L/2+b2+h−z22⋅4Rg1x−gl2x4R2Blx⋅expy+L/2+b2Inλlxh+R2+z2Inλlx−Inδlxh+dl2+−R2y−L/22⋅hh+dr⋅1y−L/22+h−z2+1y−L/22+h+z2−4zh+zy−L/22+h+z22⋅4Rgrx−gr2x4R2Brx⋅expy−L/22lnλrxh+R2+z2lnλrx−lnδrxh+dr2,(36)w4x,y,z=w4lx,y,z+w4rx,y,z=R22h−zy+L/2+b2+h−z2+h+zy+L/2+b2+h+z2−2zy+L/2+b2−h+z2y+L/2+b2+h+z22⋅4Rglx−gl2x4R2Blx⋅expy+L/2+b2lnλlxh+R2+z2lnλlx−lnδlxh+dl2+R22h−zy−L/22+h−z2+h+zy−L/22+h+z2−2zy−L/22−h+z2y−L/22+h+z22⋅4Rgrx−gr2x4R2Brx⋅expy−L/22lnλrxh+R2+z2lnλrx−lnδrxh+dr2,where gx=2R1−1−ηx, Bx=4hh+d−h+d2−ηxR+d2/RηxR+d,(37)λx=14−gxπRηxarcsindR−gx/2+1−dR−gx/22−1,δx=12−gxπR2ηxR−gx4arcsindR−gx/4,where gx and ηx are the equivalent ground loss parameter and ground loss rate, respectively, at a distance of x from the excavation surface along the tunneling direction, and Bx, λx, and δx are the calculated parameters.The three-dimensional solution introduces the offset,b, in the horizontal direction, and at the same time differentiates the calculation parameters of the first and the second shield construction to consider the mutual influence of the first and the second tunnel construction. Refer to references [32, 33] for the values of the parameters η and d. Affected by the prior construction of the right tunnel, the ground loss rate caused by the second shield construction in this study can be expressed as follows [27]:(38)η1=0.0017h+DL−0.1454ηr.


## 3.6. The Three-Dimensional Viscoelastic Solution of the Total Soil Displacement

The numerical calculation program was written in MATLAB R2017. The horizontal displacement of the soil in theY direction was obtained by superposing equations (17), (24), (32), and (35), and the vertical displacement of the soil in the Z direction was obtained by superposing equations (18), (25), (33), and (36). By combining equation (9) with the calculation method in this study, the horizontal displacement in the X direction can also be derived. However, it is generally ignored in actual construction, so this study only discusses the changes in the soil displacement values in the Y and Z directions under different X coordinates. Then the time-domain solutions of the horizontal and vertical displacement of any point, (x, y, z), in the ground caused by the additional support pressures, shield shell frictions, grouting pressures, and ground loss can be expressed as follows:(39)vx,y,z,t=v1¯x,y,z,t+v2¯x,y,z,t+v3¯x,y,z,t+v4x,y,z,wx,y,z,t=w1¯x,y,z,t+w2¯x,y,z,t+w3¯x,y,z,t+w4x,y,z.


## 4. Verification of the Time-Domain Solution

The accuracy of the viscoelastic parameters directly affects the calculation results of the viscoelastic solution proposed in this study. For a tunnel project, the viscoelastic parameter values should be determined through an indoor creep test first, and then the calculated results in this study should be compared with the field measurements to verify the rationality of the analytical solution. Based on the shield tunnel parameters and soil viscoelastic parameters proposed by Wei et al. [27, 34], the rationality of the method in this study is verified.The calculation parameters are as follows: the outer radius of the tunnelR = 3.17 m; the buried depth of the tunnel h = 19 m; the length of the shield machine J = 8.4 m; the distance between the first and second excavation face K = 10 m; the center of the two tunnels line spacing L = 12 m; the additional support pressures P1l = P1r = 20 kPa; the shield shell frictions P2l = P2r = 10 kPa; the grouting pressures P3l = P3r = 120 kPa. The viscoelastic parameters are G1 = 3 MPa, Gk = 2 MPa, K1 = 20 MPa, and ηk = 0.2 GPa·d. The distances from the soils moving the focus of the left and right tunnels to the center of the tunnel are dl = 0.24R and dr = 0.79R, respectively, the maximum settlement offset of the second tunnel b = –2.8 m, and the final ground loss rates in the left and right tunnel are ηl=0.301% and ηr=0.82%, respectively.According to equations (28) and (29), the time-domain curves of the horizontal and vertical displacement of the ground surface directly above the central axis of the two tunnels under the influence of additional support pressures, shield shell frictions, grouting pressures, and ground loss were then calculated. The results are shown in Figure 3. It can be concluded that the influence of various factors on soil displacement from high to low is a ground loss, shield shell frictions, grouting pressures, and support pressures.Figure 3
The time-domain curves of the horizontal and vertical displacement of the ground surface on the central axis of a double-line tunnel induced by various factors: horizontal displacement caused by support pressures (a), shield shell frictions (c), grouting pressures (e), and ground loss (g); vertical displacement caused by support pressures (b), shield shell frictions (d), grouting pressures (f), and ground loss (h).
(a)(b)(c)(d)(e)(f)(g)(h)From the time-domain curves shown in Figure3, it can be concluded that the three-dimensional viscoelastic analytical solution can consider the change in soil displacement over time. Under the action of various factors, the instantaneous displacement of the soil occurs when t = 0 d, which gradually increases and stabilizes with time. When t = 300 d, the soil displacement is basically stable, and the displacement value in the stable stage is approximately 2.5 times the instantaneous displacement. This calculation result is related to the value of the viscoelastic parameters.The soil displacement caused by the support pressures is distributed antisymmetrically to the position of the excavation surface (left linex = 0 m, right line x = 10 m). The soil in front of the excavation produces a slight horizontal displacement close to the central axis of the two tunnels, while the soil in the rear is the opposite. The vertical displacement of the soil induced by the support pressures is primarily manifested as the uplift at the front of the excavation and the settlement in the rear. The peaks of the horizontal and vertical displacement both appear at a position 14 m away from the excavation surface. The influence of the shield shell frictions was similar to that of the support pressures, but the soil displacement caused by the frictions was larger than that by support pressure, and the horizontal and vertical displacements in the stable stage reached 0.6 mm and 1.8 mm, respectively.The grouting pressures primarily occurred in the grouting area of the shield tail, so the horizontal and vertical displacement peaks of the soil appeared at the shield tail (left linex = –10 m, right line x = 0 m). The horizontal displacements of the left and right lines were distributed in an antisymmetrical distribution. The grouting pressures caused the soil in the middle of the two tunnels to move close to the central axis, and the peak value in the stable phase was 1.8 mm. The vertical displacement is represented by the settlement of the soil directly above the grouting area of the shield tail, with a peak value of 5.8 mm.Soil displacement caused by ground loss is the primary reason for ground deformation induced by shield tunneling. The advance shield construction on the right line caused a large horizontal displacement of the soil behind the excavation face, which was biased toward the axis of the right tunnel. The secondary disturbance caused by the tunnel construction on the left made the surface soil on the central axis move to the left again on the basis of the right displacement, but the displacement was reduced by approximately 80% compared to the right. Therefore, the horizontal displacement of the soil was inclined to the side excavated first. The vertical displacement caused by the ground loss on the left and right lines was asymmetrical. The secondary excavation caused further soil settlement on the basis of the primary excavation, which further increased the total displacement which was consistent with the research conclusions of many scholars. This phenomenon was consistent with numerous reports. Along the opposite direction of tunneling, the vertical displacement of the surface caused by ground loss gradually increased and finally reached a stable state.Many previous studies have used the equivalent ground loss to predict surface settlement, but the surface settlement obtained by this prediction method was often smaller than the measured value. This cannot explain the uplift of the soil before the excavation face. This study obtained the time-domain solution of the total displacement by superimposing the soil displacement under the influence of four factors and compared the final settlement with the measured data, as shown in Figure4. The total surface deformation curve caused by various factors was displayed as an uplift in front of the excavation face and settlement behind the face. The calculation results showed that the uplift in the soil in front of the excavation face was primarily caused by the shield shell frictions. At the back of the excavation surface, along the opposite direction of the excavation, as the distance increased, the surface settlement first increased and then gradually decreased, and finally reached a stable state. This trend was affected by the effect of shield tail grouting. The better the grouting effect, the more significant the inhibition in the soil settlement, and the more the maximum settlement value tended to be a stable value.Figure 4
Surface vertical displacement and total displacement curves generated by various factors in the stable phase (a) and comparison with the field measured data (b).
(a)(b)As the soil displacement caused by shield tunneling consisted of primarily vertical and horizontal displacement that was rarely monitored during construction, this study verified the rationality of the analytical solution by comparing the vertical displacement. Figure4(b) shows a comparison of the settlement value in the stable phase calculated by the viscoelastic solution on the y = 0 m and y = 6 m axes of the ground surface with the measured results. By comparing the results, it can be concluded that the theoretical solution was consistent with the measured value. The calculation method in this study clearly explained the phenomenon of the surface uplift in front of the shield excavation and the settlement of the rear surface under the action of multiple factors. However, the three-dimensional viscoelastic analytical solution involved many parameters, and the accuracy of the calculation result was closely related to the accuracy of the value of each parameter. Therefore, it was necessary to further analyze the influence of the parameter values on soil displacement.


## 5. Analysis of the Influence of Viscoelastic Parameters on the Vertical Soil Displacement

### 5.1. The Influence of Different Viscosity Coefficientsηk

To analyze the influence of the viscosity coefficient on the displacement of the soil, the vertical displacement was used as an example. Points (20, 0, 10) and (–10, 0, 10) in front of and behind the excavation surface, respectively, were used to calculate the settlement under the influence of different factors. The calculation results are shown in Figure5. ηk was set at 0.1 GPa · d, 0.2 GPa · d, 0.3 GPa · d, and 0.4 GPa · d. A positive value indicated that the soil had settled, and a negative value indicated that the soil had uplifted.Figure 5
Different viscosity coefficients,ηk, cause a change in the settlement value at point (20, 0, 10) and point (–10, 0, 10) with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)The time required for the settlement to stabilize increased with an increase in the viscosity coefficient, and the change of the settlement volume primarily occurred within the first 300 days. Under the action of various factors, the difference in the viscosity coefficient will not cause a change in the final settlement, and the settlement value gradually stabilizes with an increase of time. In addition, it can be seen from Figure5(d) that under the action of multiple factors, the uplift value in front of the excavation surface will gradually increase and become stable with an increase in time. With the gradual advancement of the shield, the calculation points in front of the excavation turned to the rear. When the shield tail lining segments are released, the soil behind the excavation face will undergo instant settlement. As time increases, the settlement value gradually decreases and tends to stabilize.


### 5.2. The Influence of Different Shear ModulusG1

Figure6 shows the change curve of the settlement values with the shear modulus, G1, at the point (–10, 0, 10) under different factors. G1 was set at 3 MPa, 6 MPa, 9 MPa, and 12 MPa. It was found that the shear modulus, G1, had a significant influence on the soil displacement. With an increase in the shear modulus, the settlement and uplift values of the soil caused by various factors decreased. As the selected calculation point was 10 m behind the excavation surface, this position is greatly affected by the grouting pressure. Therefore, as the uplift value of the soil caused by the grouting pressures decreased, the total settlement showed an increasing trend. In addition, it can be seen from the calculation results that the difference in the shear modulus will not cause changes in the time required for the settlement to stabilize, and the soil displacement basically reached a stable state at 400 d.Figure 6
Different shear moduli,G1, caused a change in the settlement value at point (–10, 0, 10) with time under the action of the support pressures (a), the shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


### 5.3. The Influence of Different Shear ModuliGk

Figure7 shows the change curve of the settlement value with the shear modulus, Gk, at the point (–10, 0, 10) under different factors. Gk was set at 2 MPa, 4 MPa, 6 MPa, and 8 MPa. From the calculation results, it was found that a change in the shear modulus, Gk, will not cause a change in the initial settlement value of the soil, but it will affect the soil displacement during the stable stage. As Gk increases, the settlement value or uplift value of the soil gradually decreases during the stable phase. In addition, the change in the shear modulus, Gk, has a significant effect on the time required for the soil to reach the stable stage. The larger the shear modulus, the shorter the time required. When Gk=2MPa, it was approximately 400 d, and when Gk=8MPa, it was approximately 100 d.Figure 7
Different shear moduli,Gk, cause a change in the settlement value at point (−10, 0, 10) with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


## 5.1. The Influence of Different Viscosity Coefficientsηk

To analyze the influence of the viscosity coefficient on the displacement of the soil, the vertical displacement was used as an example. Points (20, 0, 10) and (–10, 0, 10) in front of and behind the excavation surface, respectively, were used to calculate the settlement under the influence of different factors. The calculation results are shown in Figure5. ηk was set at 0.1 GPa · d, 0.2 GPa · d, 0.3 GPa · d, and 0.4 GPa · d. A positive value indicated that the soil had settled, and a negative value indicated that the soil had uplifted.Figure 5
Different viscosity coefficients,ηk, cause a change in the settlement value at point (20, 0, 10) and point (–10, 0, 10) with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)The time required for the settlement to stabilize increased with an increase in the viscosity coefficient, and the change of the settlement volume primarily occurred within the first 300 days. Under the action of various factors, the difference in the viscosity coefficient will not cause a change in the final settlement, and the settlement value gradually stabilizes with an increase of time. In addition, it can be seen from Figure5(d) that under the action of multiple factors, the uplift value in front of the excavation surface will gradually increase and become stable with an increase in time. With the gradual advancement of the shield, the calculation points in front of the excavation turned to the rear. When the shield tail lining segments are released, the soil behind the excavation face will undergo instant settlement. As time increases, the settlement value gradually decreases and tends to stabilize.


## 5.2. The Influence of Different Shear ModulusG1

Figure6 shows the change curve of the settlement values with the shear modulus, G1, at the point (–10, 0, 10) under different factors. G1 was set at 3 MPa, 6 MPa, 9 MPa, and 12 MPa. It was found that the shear modulus, G1, had a significant influence on the soil displacement. With an increase in the shear modulus, the settlement and uplift values of the soil caused by various factors decreased. As the selected calculation point was 10 m behind the excavation surface, this position is greatly affected by the grouting pressure. Therefore, as the uplift value of the soil caused by the grouting pressures decreased, the total settlement showed an increasing trend. In addition, it can be seen from the calculation results that the difference in the shear modulus will not cause changes in the time required for the settlement to stabilize, and the soil displacement basically reached a stable state at 400 d.Figure 6
Different shear moduli,G1, caused a change in the settlement value at point (–10, 0, 10) with time under the action of the support pressures (a), the shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


## 5.3. The Influence of Different Shear ModuliGk

Figure7 shows the change curve of the settlement value with the shear modulus, Gk, at the point (–10, 0, 10) under different factors. Gk was set at 2 MPa, 4 MPa, 6 MPa, and 8 MPa. From the calculation results, it was found that a change in the shear modulus, Gk, will not cause a change in the initial settlement value of the soil, but it will affect the soil displacement during the stable stage. As Gk increases, the settlement value or uplift value of the soil gradually decreases during the stable phase. In addition, the change in the shear modulus, Gk, has a significant effect on the time required for the soil to reach the stable stage. The larger the shear modulus, the shorter the time required. When Gk=2MPa, it was approximately 400 d, and when Gk=8MPa, it was approximately 100 d.Figure 7
Different shear moduli,Gk, cause a change in the settlement value at point (−10, 0, 10) with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


## 6. Analysis of the Influence of the Geometric Parameters on the Vertical Soil Displacement

### 6.1. The Influence of Different Distances from the Excavation Surface

Figure8 shows the relationship between the changes in the surface settlement caused by the various factors at different distances from the excavation face over time, where y = 6 m and z = 0 m. From the calculation results, it can be concluded that in front and behind the excavation face, as the distance from the face increased, the settlement gradually decreased. For different calculation points, the settlement value or the uplift value gradually increased and tended to be stable with an increase in time, which was consistent with the previous calculation results. The soil displacement caused by the additional support pressures and the shield shell frictions reduced by order of magnitude at a distance of 80 m (approximately four times the buried depth) from the excavation surface. The influence of grouting pressure was primarily located near the excavation face and more biased toward the back of the excavation face. This directly caused the total settlement at x = –20 m to be lower than the total displacement farther behind the excavation surface, and the uplift value at x = 20 m in front of the excavation surface is higher than that at the far front. The total displacement at a distance of 60 m from the excavation surface (approximately three times the buried depth) was primarily caused by ground loss, and this range can be considered as the area affected by the various factors.Figure 8
Different distances from the excavation face caused a change in the ground settlement with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


### 6.2. The Influence of Different Calculation Depths

Figure9 shows the relationship between the soil settlement caused by various factors at different depths with time. In this figure, the position of the excavation face on the left is x = 0 m, y = 6 m is the position of the central axis of the two tunnels, and z has four different calculation depths of 0 m, 5 m, 10 m, and 15 m. It can be seen from the calculation results that as the calculated depth gradually approaches the buried depth of the tunnel, the soil displacement caused by the support pressures first increases and then decreases, and the displacement caused by the shield shell frictions and the grouting pressures gradually increases. This phenomenon is related to the relative position of the force and the calculation points. From the calculation result of the total displacement, it can be concluded that the closer to the tunnel axis, the greater the increment of soil deformation, and the more obvious the influence of the time effect on soil deformation. As the calculated depth increases, the settlement of the surface soil gradually turns into the uplift of the deep soil, indicating that close to the tunnel axis, the grouting pressure will cause the soil to undergo a radial displacement centered on the tunnel axis. The magnitude of this displacement is related to the actual grouting pressure value.Figure 9
Different calculation depths cause a change in the ground settlement with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


### 6.3. The Influence of Different Tunnel SpacingsL

Figure10 shows the relationship between the settlement of points (–10, L/2, 0) caused by various factors under different tunnel spacings with time. With an increase in time, the settlement value or uplift value produced by various factors gradually stabilized. The soil displacement caused by the support pressures, the shield shell frictions, and the grouting pressures decreases with an increase in the distance between the tunnels. It shows that an increase in the spacing reduces the influence of the shield construction of the two tunnels on the soil displacement of the central axis.Figure 10
Different tunnel spacing causes a change in the ground settlement at points (–10,L/2, 0) with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)Figure10(d) shows the change curve of the total displacement. It can be seen from the figure that with the gradual increase in the tunnel spacing, the surface settlement on the central axis first rose slightly and then fell sharply. Furthermore, the settlement change caused by the viscoelastic properties of the soil decreases with an increase in the tunnel spacing. As mentioned above, the tendency of the total settlement to change is due to an increase in the tunnel spacing, which gradually reduces the influence of various factors on the central axis. This leads to a gradual decrease in the uplift value and settlement value.


## 6.1. The Influence of Different Distances from the Excavation Surface

Figure8 shows the relationship between the changes in the surface settlement caused by the various factors at different distances from the excavation face over time, where y = 6 m and z = 0 m. From the calculation results, it can be concluded that in front and behind the excavation face, as the distance from the face increased, the settlement gradually decreased. For different calculation points, the settlement value or the uplift value gradually increased and tended to be stable with an increase in time, which was consistent with the previous calculation results. The soil displacement caused by the additional support pressures and the shield shell frictions reduced by order of magnitude at a distance of 80 m (approximately four times the buried depth) from the excavation surface. The influence of grouting pressure was primarily located near the excavation face and more biased toward the back of the excavation face. This directly caused the total settlement at x = –20 m to be lower than the total displacement farther behind the excavation surface, and the uplift value at x = 20 m in front of the excavation surface is higher than that at the far front. The total displacement at a distance of 60 m from the excavation surface (approximately three times the buried depth) was primarily caused by ground loss, and this range can be considered as the area affected by the various factors.Figure 8
Different distances from the excavation face caused a change in the ground settlement with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


## 6.2. The Influence of Different Calculation Depths

Figure9 shows the relationship between the soil settlement caused by various factors at different depths with time. In this figure, the position of the excavation face on the left is x = 0 m, y = 6 m is the position of the central axis of the two tunnels, and z has four different calculation depths of 0 m, 5 m, 10 m, and 15 m. It can be seen from the calculation results that as the calculated depth gradually approaches the buried depth of the tunnel, the soil displacement caused by the support pressures first increases and then decreases, and the displacement caused by the shield shell frictions and the grouting pressures gradually increases. This phenomenon is related to the relative position of the force and the calculation points. From the calculation result of the total displacement, it can be concluded that the closer to the tunnel axis, the greater the increment of soil deformation, and the more obvious the influence of the time effect on soil deformation. As the calculated depth increases, the settlement of the surface soil gradually turns into the uplift of the deep soil, indicating that close to the tunnel axis, the grouting pressure will cause the soil to undergo a radial displacement centered on the tunnel axis. The magnitude of this displacement is related to the actual grouting pressure value.Figure 9
Different calculation depths cause a change in the ground settlement with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)


## 6.3. The Influence of Different Tunnel SpacingsL

Figure10 shows the relationship between the settlement of points (–10, L/2, 0) caused by various factors under different tunnel spacings with time. With an increase in time, the settlement value or uplift value produced by various factors gradually stabilized. The soil displacement caused by the support pressures, the shield shell frictions, and the grouting pressures decreases with an increase in the distance between the tunnels. It shows that an increase in the spacing reduces the influence of the shield construction of the two tunnels on the soil displacement of the central axis.Figure 10
Different tunnel spacing causes a change in the ground settlement at points (–10,L/2, 0) with time under the action of the support pressures (a), shield shell frictions (b), grouting pressures (c), and total displacement over time (d).
(a)(b)(c)(d)Figure10(d) shows the change curve of the total displacement. It can be seen from the figure that with the gradual increase in the tunnel spacing, the surface settlement on the central axis first rose slightly and then fell sharply. Furthermore, the settlement change caused by the viscoelastic properties of the soil decreases with an increase in the tunnel spacing. As mentioned above, the tendency of the total settlement to change is due to an increase in the tunnel spacing, which gradually reduces the influence of various factors on the central axis. This leads to a gradual decrease in the uplift value and settlement value.


## 7. Conclusion

Based on the Mindlin solution and the Boltzmann viscoelastic model, this study comprehensively considered the influence of the additional support pressures, the shield shell frictions, the grouting pressures, and the ground loss, and proposed a three-dimensional viscoelastic solution of the soil displacement caused by shield tunneling. The following conclusions were drawn:(1)
The calculation model considered the nonsynchronous construction of a double-line tunnel and the mutual influence between the two tunnels. It can be concluded from the calculation examples that the three-dimensional viscoelastic analytical solution can reflect the change in soil displacement caused by various factors with time. Under the action of various factors, whenT = 0 d, the soil had an instantaneous displacement that gradually increased with time and tended to be stable. The displacement in the stable phase was approximately 2.5 times the instantaneous displacement. The time required to reach a stable stage was related to the value of the viscoelastic parameters. The results obtained by the calculation method in this study were consistent with the measured results. In addition, it clearly explained the phenomenon of the surface uplift in front of the excavation surface and the settlement of the rear surface under the action of multiple factors.(2)
This study discussed the influence of changes in viscoelastic parameters on soil displacement. It can be concluded that the time required for the settlement to stabilize increased with an increase in the viscosity coefficient, and the change of the settlement volume primarily occurred within the first 300 days. With an increase in the shear modulus,G1, the settlement value and uplift value of the soil caused by various factors decreased, whereas the time required for the soil displacement to reach the stable stage was basically the same, approximately 400 d. The change in the shear modulus, Gk, will not cause a change of the initial settlement value of the soil, but it will affect the soil displacement value during the stable phase and the time it takes to reach the stable phase. The larger the shear modulus, Gk, the smaller the settlement value or uplift value of the soil in the stable stage, and the shorter the time required to reach stability.(3)
This study further discussed the influence of a change in the geometric parameters on soil displacement. In front of and behind the excavation surface, as the distance from the excavation surface increased, the settlement gradually decreased. The total displacement at a distance of 60 m from the excavation surface (approximately three times the buried depth) was primarily caused by ground loss, and this range can be considered as the area affected by various factors. The closer the calculated depth was to the buried depth of the tunnel, the larger the increment of soil deformation. Furthermore, the more obvious the influence of the time effect on soil deformation was. An increase in the tunnel spacing gradually reduced the influence of various factors on the central axis, which led to a gradual decrease in the uplift value and settlement value.


## 8. Discussion

(1)
This study aimed to propose a displacement calculation method considering the creep characteristics of soil, which was derived based on the Boltzmann viscoelastic model. For different geological conditions, it was necessary to conduct indoor tests and judge the applicability of the Boltzmann viscoelastic model. It could be considered to use different models combined with the solution method proposed in this article to obtain analytical solutions suitable for specific site conditions to meet engineering needs.(2)
The focus of this study was to analyze the viscoelastic effect and the time-domain deformation of the soil caused by shield tunneling. It should be noted that the analytical solution found in this study was limited to the soil displacement caused by construction factors. This study did not consider the influence of external factors, such as the discontinuity of the construction, fluid-solid coupling, and the influence of secondary grouting. A follow-up study should be performed that is combined with more detailed field measurements of soil displacement throughout the entire construction process for a more accurate assessment.



---

*Source: 8839010-2020-12-07.xml*