Fractional viscoelastic analytical solution for the ground displacement of a shallow tunnel based on a time-dependent unified displacement function

Computers and Geotechnics 117 (2020) 103284

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Fractional viscoelastic analytical solution for the ground displacement of a shallow tunnel based on a time-dependent unified displacement function

T

⁎

Dechun Lua, Fanchao Konga, Xiuli Dua, , Chenpeng Shena, Cancan Sua, Jun Wangb a b

Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, Beijing University of Technology, Beijing 100124, China Key Laboratory of Engineering and Technology for Soft Foundation and Tideland Reclamation, Wenzhou University, Wenzhou, Zhejiang 325035, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Fractional Generalized Kelvin constitutive model Time-dependent unified displacement function Complex variable method Fractional viscoelastic analytical solution Analyses of the time effect of ground displacement field

A time-dependent unified displacement function expressed by a Fourier series is used as the displacement boundary condition of the tunnel cross-section (BCTC). Abel element expressed by the fractional differential is adopted to reflect the rheological behaviours of the surrounding rock. Using Laplace transform and the inverse transformation, the relaxation modulus of the Fractional Generalized Kelvin constitutive model for one-dimensional (1D) is obtained. The effects of the fractional order and viscous material parameters on the relaxation modulus are analysed. 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. Two factors are considered to capture the time effect of ground deformation: a time-dependent displacement BCTC and time-dependent mechanical parameters. A parameter study is conducted to investigate the time effect of the ground displacement. The deformation parameters of the tunnel cross-section are first calibrated by a part of the field data for the Longdong tunnel and Gaoqiao tunnel. Based on the calibrated model, the proposed method can well predict ground displacement due to the time effect.

1. Introduction Most types of geomaterials exhibit a significant time-dependent response [1,2]. The time-dependent rheological property of the ground deformation is an important problem in the tunnel projects. There are generally three methods to investigate the displacement field of the surrounding rock caused by the time effect: empirical methods [3], numerical methods [4], and analytical methods [5,6]. The analytical methods based on the rigorous mathematical deduction can quantitatively consider the effects of geological parameters and geometrical parameters. The ground time-dependent behaviours for a shallow tunnel can be reflected by the time-dependent boundary conditions of the analytical models and time-dependent material parameters. The motivation for the presented work is to analytically obtain the ground displacement field of shallow tunnels due to the time effect. The viscoelastic analytical solution is usually presented by following two steps at present [7]: (1) to obtain the elastic analytical solution corresponding to the mechanical analysis model; (2) to perform the Laplace computation for the functions of the obtained elastic solution only including the material parameters, or simultaneously including material parameters and the time-dependent boundary conditions.

The reasonable boundary conditions are the key factors for determining the elastic and viscoelastic solutions of a shallow tunnel. For the mechanical analysis model of the shallow tunnel, there are two kinds of the boundary conditions for the tunnel cross-section, i.e., stress BCTC [8] and displacement BCTC [9,10]. In recent years, the reasonable stress boundary functions [11] and displacement boundary functions [12,13] of the tunnel cross-section are constantly proposed. The corresponding elastic solutions are obtained by the complex variable method [14,15]. Although some works are done on the elastic solution of a shallow tunnel, few works have been carried out on the viscoelastic analytical solutions. The viscoelastic analytical solutions for a shallow tunnel considering the surcharge loads of the ground surface are determined in the Wang et al., works, where 0 stress [16] and time-dependent uniform support force [17] are used as the stress BCTC, respectively. In the tunnel engineering, the displacement monitoring is easier than the stress monitoring, and the accuracy of displacement monitoring is higher. However, there are no works on the viscoelastic analytical solutions of the shallow tunnel for the displacement BCTC. At present, the time-dependent material properties can be accounted for by the classical viscoelastic constitutive models, which are made up of all types of series and parallel connections of the elastic and

⁎

Corresponding author. E-mail addresses: [email protected] (D. Lu), [email protected] (F. Kong), [email protected] (X. Du), [email protected] (C. Shen), [email protected] (C. Su). https://doi.org/10.1016/j.compgeo.2019.103284 Received 1 August 2019; Received in revised form 26 September 2019; Accepted 27 September 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Newton viscous elements. In particular, the function of the Newton viscous element is 1 order differential between the strain and the time. Wang et al., [16] considered the classical Generalized Kelvin constitutive model to obtain the viscoelasticity analytical solution of the shallow tunnel by the correspondence principle. A fractional differential relation between the strain and the time can be established by a fractional element, namely Abel viscous element [18,19]. The models consisting of elastic and Abel viscous elements are called as fractional viscoelastic constitutive model. The relaxation and creep functions are modelled in a type of the power functions on the basis of the fractional viscoelastic constitutive model, which can reasonably capture the rheological behaviours of materials [20]. The fractional order reflects the rheological response of different materials. In addition, the classical viscoelastic constitutive model is only a special case of the fractional viscoelastic constitutive model when the fractional order is equal to 1. Unfortunately, the fractional viscoelastic constitutive models are not considered in the analytical solution of shallow tunnels. A fractional viscoelastic constitutive model for soil is considered to more reasonably describe the rheological behaviour of the ground. A time-dependent unified displacement function of the tunnel cross-section is proposed to reflect the non-symmetrical deformations on the horizontal and vertical centrelines. A fractional viscoelastic analytical solution of the ground displacement is presented by combining the complex variable method and the correspondence principle. The time effect of the ground displacement is reflected by two factors: the timedependent material parameters and time-dependent displacement BCTC. The rationality of the proposed method is verified by the field data of the Longdong tunnel and the Gaoqiao tunnel. 2. Fractional viscoelastic constitutive model for the surrounding rock In recent years, fractional calculus has been widely used in mechanical analyses [21,22]. In this work, the relaxation modulus of 1D based on the Fractional Generalized Kelvin constitutive model is obtained. The effects of the fractional order and viscous parameters on the relaxation modulus are analysed. In the viscoelastic constitutive models, the spring is used as an elastic element as shown in Fig. 1(a). The linear relationship between stress and strain is expressed as follows,

σ = Eε

(1)

where σ is the stress; ε is the strain; E is Young elastic modulus. When the Newton element is used as the viscous element as shown in Fig. 1(b), the relationship between stress and strain is expressed as follows,

σ=η

dε dt

Fig. 2. Creep and relaxation curves of Abel viscous model.

In Eq. (3), when σ remains constant, the creep behaviour of materials can be described. According to the definition of Riemann-Liouville fractional integral, the fractional integral is conducted for Eq. (3),

(2)

where η is the viscous coefficient of the Newton element; t is the time. When Abel element is used as the viscous element as shown in Fig. 1(c), the relationship between stress and strain can be expressed by a fractional differential as follows,

d με σ = ηA μ dt

ε (t ) =

σ0 tμ ηA Γ(1 + μ)

(4)

where 0⩽ μ ⩽ 1. In particular, the time fractional derivatives are taken over the total time of the analytical process, and the “full” memory [23,24] is adopted in this work; thus, the left terminal in the fractional operator is equal to 0. Γ(x ) is the gamma function and is expressed as follows,

(3)

where μ is the fractional order; ηA , which is the viscous coefficient of Abel element, has an anomalous dimension and depends on the value of μ [20]. In particular, when μ is equal to 0, Eq. (3) reduces to Eq. (1). When μ is equal to 1, Eq. (3) reduces to Eq. (2).

Γ(x ) =

∫0

∞

e−ξ ξ x − 1dξ

(5) Fig. 1. Viscous and elastic elements.

2

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

The recursive nature can be observed for the gamma function, i.e., Γ(x + 1) = x Γ(x ) . For example, ηA = 100 MPa Dayμ and σ0 = 1 MPa; for different fractional orders, the relationship between ε (t ) and t is shown in Fig. 2(a). When μ is equal to 0, ε remains constant, which is determined by σ0 and E. When μ is equal to 1, ε linearly increases with the increase of t. When μ is greater than 0 and less than 1, the nonlinear rheological behaviours of the materials are shown in Fig. 2(a). ε nonlinearly increases, and the degree of the increase becomes small with the increase of t. When ε remains constant in Eq. (3), based on the fractional integral of Riemann-Liouville, the relaxation behaviour of the materials is expressed as follows:

σ (t ) = ηA ε0

t −μ Γ(1 − μ)

(6) μ

For example, ηA = 100 MPa Day and ε0 = 0.05; for different fractional orders, the relationship between σ (t ) and t is shown in Fig. 2(b). When μ is equal to 1, dε d t is equal to 0, so σ is equal to 0. When μ is equal to 0, σ remains constant. When μ is greater than 0 and less than 1, for the same μ , σ first rapidly decreases and subsequently tends to a stable value with the increase of t. For the same t, a smaller value of μ corresponds to a larger value of σ , as shown in Fig. 2(b). The classical viscoelastic constitutive models are established by the series or parallel connection of spring elements and Newton viscous elements. Similarly, the fractional viscoelastic constitutive models are obtained by the series or parallel connection of the spring elements and Abel viscous elements. The Fractional Maxwell constitutive model, Fractional Kelvin constitutive model, Fractional Generalized Kelvin constitutive model and Fractional Burgers constitutive model are shown in Fig. 3. Taking the Fractional Generalized Kelvin constitutive model as an example, a method is provided to obtain the relaxation modulus. The 1D constitutive relation of the Fractional Generalized Kelvin constitutive model is expressed as,

d με 1 d μσ E E + EK + Kε= + M σ dt μ η′K EM dt μ η′K EM

(7)

where EM and EK are Young elastic modulus of two spring elements as shown in Fig. 3(c), and ηK′ is the viscous parameter. The method of determining parameters, EM , EK , ηK′ and μ are provided as follows. The stress σ is applied on the specimen at t = 0. At this point, the viscous element will not deform in time. The instantaneous strain ε of a spring element in series is produced. EM is determined by

EM =

σ ε

Fig. 4. Effects of μ and ηκ on Es(t).

transformation of the fractional differential operator of RiemannLiouville for Eq. (7),

(8)

EM can be considered as elastic modulus in the standard geotechnical investigation. The strain ε0 is generated under the external force at t = 0. Thereafter, the external force is suddenly released. The unloading curve can be determined to obtain the relationship between the strain and the time as follow, m

ε = ε0

∑ k=1

t kμ − 1 (−EK η) k − 1 Γ(kμ)

EM ⎧ Es (t ) = EM 1 − ⎨ EK + EM ⎩

Laplace

n=0

(−1)n [(EK + EM ) η′K ]n + 1t (n + 1) μ ⎫ ⎬ Γ[(n + 1) μ + 1] ⎭

(10)

For example, EM = 10 MPa, EK = 5 MPa, η'κ = 100 MPa Dayμ, and the effect of μ on Es(t) is shown in Fig. 4(a). For the same μ, Es(t) gradually decreases and the degree of the decrease becomes small with the increase of t. When μ = 0.5, the effect of η'κ on Es(t) is shown in Fig. 4(b). For the same η'κ, a smaller t corresponds to a faster decrease of Es(t).

(9)

EK, η'κ and μ can be determined by the test results. The relaxation modulus is obtained by

∞

∑

inverse

Fig. 3. Fractional viscoelastic constitutive models. 3

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

3. Time-dependent unified displacement function The consolidation and creep of the ground causes time-dependent deformation of the surrounding ground mass. The time span of ground deformation may be from the moment that the tunnel begins to be excavated to forever. For each moment, the asymmetrical displacement on the horizontal and vertical centre lines of the tunnel cross-section is always formed. In particular, when the tunnel lining contacts with the surrounding ground mass, the tunnel cross-section in this work refers to the interface between the ground and the outer of the lining. When the lining is not in contact with the surrounding ground mass, the tunnel cross-section refers to the excavation section of the soil. To quantitatively express this asymmetrical deformation of the tunnel cross-section, Eq. (11) is proposed to express the deformation of the tunnel crosssection over time, ∞

u (θ′, t ) = −u 0 (t ) +

∑ [a′n (t ) sin(nθ′) + b′n (t ) cos(nθ′)] n=1

(11)

After the construction of the whole tunnel has been completed, the time effect of the ground deformation is focused on in this work. The moment for the completion of the whole tunnel construction is defined as t = 0. The influence of the whole construction process is comprehensively reflected by the displacement BCTC at t = 0. When n is equal to 2, Eq. (11) is simplified as,

u (θ′, t ) = −u 0 (t ) + a1′ (t ) sin θ′ + b1′ (t ) cos θ′ + a2′ (t ) sin 2θ′ + b2′ (t ) cos 2θ′ (12) At any moment, the effects of five deformation modes on the tunnel contour are shown in Fig. 5. The basis deformation modes a1′ (t ) sin θ′ or a2′ (t ) sin 2θ′ are involved in the displacement BCTC to reflect the asymmetrical deformation on the horizontal centre line. The basis deformation modes b1′ (t ) cos θ′ or a2′ (t ) sin 2θ′ are involved in the displacement BCTC to reflect the asymmetrical deformation on the vertical centre line owing to the non-symmetrical geological contours on the vertical centrelines and the bias load on the ground surface or underground. The time-dependent unified displacement function can reflect the dynamic change of the tunnel cross-section deformation over time. 4. Viscoelastic analytical solution of the ground displacement The displacement BCTC is the result of the interaction between the ground and the lining. When the displacement BCTC is used as input of an analytical solution, the time-dependent ground displacement field is paid attention to in this work by combining the correspondence principle and the complex variable method. In particular, the assumptions in this work are as follows: (1) on the basis of the plane strain condition; (2) the surrounding rock is considered as isotropic homogeneous viscoelastic material. 4.1. Elastic analytical solution The mechanical model of a shallow tunnel is shown in Fig. 6. γ and ν are the uniform volumetric weight, and Poisson’s ratio of the surrounding rock, respectively; r0 and h are the radius and depth of the tunnel, respectively. G is the shear modulus. The displacement at any point in the ground is expressed as [8],

2G (u x + iu y ) = κφ (z ) − zφ′ (z ) − ψ (z )

(13)

Fig. 5. Contour of tunnel cross-section at any moment when series expansion n is equal to 2.

where κ = 3 − 4ν. In the z plane, zero stress is used as the boundary conditions of the surface. Eq. (12) is considered not to change over time and is used as BCTC. The solution process of the elastic solution for a shallow tunnel is as follows.

z = z : φ (z ) + zφ′ (z ) + ψ (z )=0

4

(14)

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 6. Mechanical model of the shallow tunnel excavation.

|z + ih| = r0: 2G (u x + iu y ) = κφ (z ) − zφ′ (z ) − ψ (z )

viscoelastic analytical solution with displacement BCTC be obtained. In particular, that Poisson’s ratio is equal to 0.5 corresponds to the undrained conditions. A similar consideration also appears in Wang et al., work. [16] When k in Eqs. (A5) and (A6) is very large, Laurent series satisfies the convergence, and ak + 1 and bk + 1 are equal to 0. The homogeneous form of the system of Eqs. (A5) and (A6) admits a solution [25], where,

(15)

When the “Buoyancy effect” is considered, φ (z ) and ψ (z ) are written as,

φ (z ) = −

ψ (z ) =

Fx + iFy 2π (1 + κ ) Fx − iFy

2π (1 + κ )

[κlog (z − z c ) + log (z − z c )] + φ0 (z )

[log (z − z c ) + κlog (z − z c )] + ψ0 (z )

(16)

(17)

where Fx and Fy are the concentration forces in the directions of x and y axes, respectively. Fx is equal to 0. Fy is equal to (1 − wt)γπr02. wt is the ratio of the average unit weight of the installed tunnel to the average weight of the soil [13,14]. Using the conformal transformation function, the physical plane is mapped to the annulus domain plane as shown in Fig. 6,

z = w (ζ ) = −ia

1+ζ 1−ζ

1 − α2 h 1 + α2

(18) 2α 1 + α2

r0 . h

∞

∞

∑ ak ζ k + ∑ bk ζ −k k=1

k=1

∞

ψ0 (z ) = c0 +

∑ ck k=1

+

∑ dk

bk = G·bk′ + fx , y, bk

(24)

c0 = G·f1 (a0′, a1′, b1′) + fx , y, c0

(25)

ck = G·f2 (ak′, bk′) + fx , y, ck

(26)

dk = G·f3 (ak′, bk′) + fx , y, dk

(27)

u x + iu y =

k=1

κφ (z ) − zφ′ (z ) − ψ (z ) 2G

(28)

When the “Buoyancy effect” is not considered, i.e., fx , y is equal to 0 in Eqs. (22)–(27); combining Eqs. (22)–(28), it can be seen that the viscoelastic analytical solution of the ground displacement is only related to the time-dependent deformation parameters of the tunnel crosssection, and irrelevant to G. When the “Buoyancy effect” is considered, the time effect of the ground displacement is only related to G by Eqs. (22)–(28). The Laplace transform and inverse transform need to be conducted for the items including the “Buoyancy effect”. The constitutive relation of the linear rheological models in the complex stress state is obtained by generalizing the constitutive relation of the uniaxial stress state. The Fractional Generalized Kelvin constitutive model of 3D can be expressed as,

(20)

4.2. Fractional viscoelastic analytical solution Based on the obtained elastic analytical solution in the Section 4.1, the Laplace transform and inverse transformation need to be conducted for the parameters of φ (z ) and ψ (z ) only including G and ν, or simultaneously including G, ν and the time-dependent unified displacement function. The expression of the viscoelastic Poisson's ratio ν(t) is,

3K (t ) − 2G (t ) 6K (t ) + 2G (t )

(23)

ζ −k

The method to obtain a0, ak, bk, c0, ck, and dk is provided in the Appendix A.

ν (t ) =

ak = G·ak′ + fx , y, ak

(19)

∞

ζk

(22)

where a0′, ak′ and bk′ are functions of α and the deformation parameters of the tunnel cross-section; f1, f2 and f3 are functions of a0′, ak′ and bk′; fx , y is the item including the “Buoyancy effect”. Eq. (13) is also expressed as,

= and where a = φ0 (ζ ) and ψ0 (ζ ) are represented by the Laurent series expansions, φ0 (ζ ) = a0 +

a0 = G·a0′ + fx , y, a0

P·sij = 2Q·eij (21)

(29)

where sij and eij are the deviatoric stress and deviatoric strain, respectively, and i, j = 1, 2, 3.

where K(t) is volumetric modulus; G(t) is shear modulus. Because the functions that consist of Poisson's ratio cannot be separated from the elastic solution of the shallow tunnel in the Section 4.1. Therefore, only when K(t) is considered the infinity at every moment, the result of the Laplace calculations for ν(t) remains constant with time and is equal to 0.5 whatever the expression of G(t) is. Only under this condition can the

P=1+

Q= 5

ηK sμ GM + GK

s μGM ηK GM GK + GM + GK GM + GK

(30)

(31)

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 7. Effect of μ on L−1 {1 [s·(Q P )]} .

where s is Laplace operator. The shear modulus G(t) is obtained according to Eqs. (29)–(31),

G (t ) =

sij 2eij

=

Q P

(32)

By Eqs. (30)–(32), the viscoelastic analytical solution of ground displacement is

u x + iu y =

[κφ (z ) − zφ′ (z ) − ψ (z )] −1 L {1 [s·(Q P )]} 2

(33)

n

where L−1 {1 [s·(Q P )]} =

∞ 1 1 + η ∑ GM K n=0 −1

When μ is equal to 1, L

L−1 {1 [s·(Q P )]} =

⎛− GK ⎞ t (n + 1) μ ⎝ ηK ⎠ . Γ[(n + 1) μ + 1]

⎜

⎟

{1/[s⋅(Q/P)]} can be simplified as follows,

1 1 1 G + − exp ⎜⎛− K t ⎞⎟ GM GK GK ⎝ ηK ⎠

(34)

t 0

which is identical to ∫ H (τ )dτ in this work [16]. It can justify the rationality of the proposed method. For example, GM = 10 MPa, GK = 5 MPa, ηK = 100 MPa·Day μ , and the effect of μ on L−1 {1 [s·(Q P )]} is shown in Fig. 7. For the same μ, L−1 {1 [s·(Q P )]} increases, and the degree of the increase decreases with the increase of t. When t is equal to 0, L−1 {1 [s·(Q P )]} is identical for different values of μ. Using the above method, the fractional viscoelastic analytical solution is presented when the time-dependent unified displacement function is set as the BCTC. 5. Analyses of the ground displacement field The common mechanical parameters of the surrounding rock and common geometrical parameters of the tunnel are used as an example to analyse the ground displacement. The geometrical parameters of the tunnel are: r0 = 3 m; h = 10 m; h′ = h − r0 ; the mechanical parameters of the surrounding rock are: γ = 17.64 kN/m3; wt = 0.3; ν =0.5; GM = 10 MPa; GK = 5 MPa; ηK = 100 MPa·Day μ . The time effects of the ground displacement field are analysed for two conditions: constant and time-dependent displacement BCTC over time. 5.1. Constant displacement BCTC Fig. 8. Shape of the surface settlement for different time.

When the stiffness of the supporting structure is notably large, the deformation of the tunnel cross-section has completely developed when the entire tunnel cross-section is completely excavated. The displacement BCTC is considered a rigid displacement boundary and does not change over time, the deformation parameters of Eq. (12) (u 0 (t ) , a1′ (t ) , 6

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 9. Ground settlement trough curves for different time.

caused by the “Buoyancy effect” are coupled to form a time-dependent displacement field. The settlement of ground surface centre point o gradually decreases as shown in Fig. 10(b), which depends on the degree of the increase of L−1 {1 [s·(Q P )]} over time.

b1′ (t ) , a2′ (t ) , b2′ (t ) ) are constant as follows to analyse the ground displacement, u (θ′, t ) = −0.075 − 0.025 sin θ′ − 0.05 cos θ′ + 0.025 sin 2θ′ + 0.05 cos 2θ′ (35) For different moments, the effects of the fractional order on the shape of the surface settlement are shown in Fig. 8. For each moment, the ground displacement field consists of two parts: downward settlement induced by the displacement BCTC and upper displacement due to the “Buoyancy effect”. The change in ground displacement is only affected by the material parameters due to the constant displacement boundary condition over time in this section. As expressed in Eq. (33), L−1 {1 [s·(Q P )]} is identical for different fractional order μ when t is equal to 0, so the shape of the surface settlement is identical in Fig. 8(a). When t increases, a larger μ corresponds to a larger L−1 {1 [s·(Q P )]} . Therefore, the upper displacement induced by the “Buoyancy effect” increases over time. As shown in Fig. 8(b) and (c), a larger μ corresponds to a more obvious uplift of the ground surface. For different moments, the effects of the fractional order on the ground settlement trough curves are shown in Fig. 9. At the same depth, the fractional order does not affect the ground settlement value when t is equal to 0, as shown in Fig. 9(a). For identical t and μ, with the increase of the depth, the maximum settlement value gradually increases, and the width of the curve of the ground settlement decreases. For the same t and depth, a larger μ corresponds to a smaller maximum settlement and a narrower curve of the ground settlement. When μ increases, the phenomenon of the uplift of the settlement curve becomes obvious. The deeper ground has more obvious uplift, as shown in Fig. 9(b) and (c). For example, with μ = 0.5, and the shape of the surface settlement and the settlement curve of ground surface centre point o over time are shown in Fig. 10(a) and (b), respectively. The maximum settlement decreases, and the degree of the decrease becomes small with the increase of t as shown in Fig. 10(a). The ground settlement caused by the disturbance of the tunnel excavation and the rebound deformation

5.2. Time-dependent displacement BCTC The deformation of the tunnel cross-section continuously develops when the stiffness of the supporting structure is not large. In the parameter analysis of this section, the time-dependent displacement BCTC is considered a simple mode as expressed in Eq. (36),

u (θ′, t ) = (−u 0 + a′1 sin θ′ + b′1 cos θ′ + a′2 sin 2θ′ + b′2 cos 2θ′) λ (t ) (36) In the analysis of this section, u 0 = 0.075, a1′ = −0.025, b1′ = −0.05, a2′ = 0.025, b2′ = 0.05. The change of each deformation parameter over time is consistent with the hyperbolic laws λ (t ) ,

λ (t ) = 1 +

va t (vb + t )

(37)

where va and vb are the parameters of function λ (t ) . The effects of va and vb on λ (t ) are shown in Fig. 11. va controls the limit value of the curve, and vb reflects the initial slope of the curve. For example, with va = 1.0 and vb = 3.0, the deformation of the tunnel cross-section, contour line of the ground displacement, shape of the surface settlement and settlement curve of the surface centre point over time are shown in Figs. 12–14. There is always asymmetrical deformation of the tunnel on the vertical and horizontal centre lines in different moments as shown in Fig. 12. a1′ sin θ′ and a2′ sin 2θ′ cause the asymmetrical deformation on the horizontal centre line; b1′ cos θ′ and a2′ sin 2θ′ cause the asymmetrical deformation on the vertical centre line. The deformation of the tunnel increases, and the degree of the increase becomes small with the increase of t, which is determined by λ (t ) . The ground displacement first rapidly increases, subsequently 7

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 10. Shape of the surface settlement and settlement curve of centre point o over time when displacement BCTC is constant. Fig. 11. Effects of va and vb on λ (t ) .

slowly increases and finally stabilizes with the increases of t, as shown in Figs. 13 and 14, which is determined by the coupling effects of the time-dependent unified displacement function and the “Buoyancy effect”. The ground displacement mainly concentrates on the upper left side of the tunnel as shown in Fig. 13. For each moment, the position corresponding to the maximum settlement of ground surface always lies on the left side as shown in Fig. 14(a). The settlement of ground surface central point o rapidly increases before approximately Day 20. At approximately Day 50, the settlement of point o reaches 0.1183 m. From Day 50 to Day 100, the settlement of point o hardly changes as shown in Fig. 14(b).

6.1. Longdong tunnel The mechanical parameters are considered as follows to analyse the ground displacement of the Longdong tunnel, i.e., γ = 18 kN/m3, ν = 0.5, GM = 8 MPa, GK = 15 MPa, μ = 0.5, ηκ = 150 MPa Day0.5 and wt = 0.7. The radius of the tunnel is 3.1 m. The depths of the two sections are 8.37 m and 10.37 m, respectively. For each section, the field monitoring data of surface settlement for four different moments, namely Day 30, Day 90, Day 150 and Day 360, are provided. With the section of h = 8.37 m as an example, the method of calibrating parameters are provided as follows. The deformation parameters of the tunnel cross-section at Day 30 are calibrated by the inversion method when the surface settlement values obtained by the analytical solution are close to the 13 field monitoring data, namely labelled points 1–13 in red circles shown in the Fig. 15(a). The deformation function of the tunnel cross-section at Day 30 is,

6. Analysis of engineering cases To access the potential of the proposed method in predicting the ground displacement owing to the time effect, Gaoqiao tunnel [26] and two sections of the Longdong tunnel [27] are used as the case study. The three-bench seven-step excavation method and the earth pressure balanced shield method are adopted in the Gaoqiao tunnel and the Longdong tunnel, respectively. The deformation law of the tunnel crosssection with time is considered to satisfy the Eq. (36). The monitoring points labelled in red circles are used to calibrate the deformation parameters of the cross-tunnel section. Other monitoring points are used to verify the rationality of the proposed method.

h = 8.37m: u (θ′, 30) = −0.03 − 0.021 sin θ′ + 0.033 cos 2θ′

(38)

Due to the symmetrical deformation of the ground surface, the asymmetrical deformation parameters on the y axis (b1′ cos θ′ and a2′ sin 2θ′) are equal 0. Considering Eqs. (36) and (38), the deformation function of the tunnel cross-section can also be expressed as at Day 30, 8

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 12. Deformation of the tunnel for different moments.

of the centre point by the proposed method for Day 90 and Day 150, namely labelled points 14 and 15, are equal to the field data. By combining the Eqs. (39)–(41), the following relations are obtained,

− 0.03 − 0.021 sin θ′ + 0.033 cos 2θ′ 30va ⎤ = (−u 0 + a′1 sin θ′ + b′2 cos 2θ′) ⎡1 + ⎢ ( v + 30) ⎥ b ⎦ ⎣

(39)

For Day 90 and Day 150, the boundary function of the ground surface is expressed as,

90va ⎤ 30va ⎤ ⎡1 + = 1.34 ⎡1 + ⎢ ⎢ (vb + 90) ⎥ (vb + 30) ⎥ ⎣ ⎦ ⎦ ⎣

(42)

(−0.03 − 0.021 sin θ′ + 0.033 cos 2θ′) x 90

150va ⎤ 30va ⎤ ⎡1 + = 1.48 ⎡1 + ⎢ ⎢ (vb + 150) ⎥ (vb + 30) ⎥ ⎣ ⎦ ⎦ ⎣

(43)

90va ⎤ = (−u 0 + a′1 sin θ′ + b′2 cos 2θ′) ⎡1 + ⎢ (vb + 90) ⎥ ⎦ ⎣

(40)

By Eqs. (42) and (43), va = 2.32 and vb = 54. Considering Eq. (39), u 0 , a1′ and b2′ can be determined.

(−0.03 − 0.021 sin θ′ + 0.033 cos 2θ′) x150 150va ⎤ = (−u 0 + a′1 sin θ′ + b′2 cos 2θ′) ⎡1 + ⎢ (vb + 150) ⎥ ⎦ ⎣

h = 8.37m: u (θ′, t ) 2.32t ⎤ = (−0.0166 − 0.0116 sin θ′ + 0.0182 cos 2θ′) ⎡1 + ⎢ (54 + t) ⎥ ⎦ ⎣

(41)

when x 90 and x150 are equal 1.34 and 1.48, respectively, the settlement 9

(44)

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 13. Contour line of the ground displacement in different moments.

ν = 0.5, GM = 15 MPa, GK = 10 MPa, μ = 0.5, ηκ = 150 MPa Day0.5, and wt = 0.5. The depth of Gaoqiao tunnel is 40 m. The tunnel crosssection is approximately circular, and its radius is 8 m. At Day 0, when the analytical results are close to the field data of the surface settlement of the labelled points 1–9, the deformation function of the tunnel cross-section are calibrated as expressed in Eq. (46).

When t are equal to Day 90 and Day 150, respectively, Eq. (44) is taken as the BCTC, the values of the surface settlement by the analytical solution are also in close agreement with the field data of other positions besides points 14 and 15. To further verify the rationality of the proposed method, the surface settlement curve by the calibration model and the field data are quantitatively consistent at Day 360 as shown in Fig. 16, which indicate that the analytical solution in this work can well predict the ground displacement induced by the time effect. For the section of h = 10.37 m, the time-dependent deformation function of the tunnel cross-section is also determined by the above same method,

u (θ′, 0) = −0.20 − 0.05 sin θ′ + 0.20 cos 2θ′

(46)

The deformation function of the tunnel cross-section at Day 5 can be expressed as follows,

u (θ′, 5) = (−0.20 − 0.05 sin θ′ + 0.20 cos 2θ′)·λ (5)

h = 10.37m: u (θ′, t )

(47)

For Day 90, Day 150 and Day 360, the analytical results are close agreement with the field data as shown in the Figs. 15(b) and 16.

When λ (5) is equal to 1.46, the settlement of the centre point by the analytical solution is equal to the field data of the labelled point 10. When Eq. (47) is used the displacement BCTC at Day 5, the field data of other positions approaches to the analytical results as shown in Fig. 17, which can verify the rationality of the proposed method.

6.2. Gaoqiao tunnel

7. Conclusions

The mechanical parameters are considered as follows to analyse the ground displacement of the Gaoqiao tunnel, i.e., γ = 17.64 kN/m3,

The fractional viscous element is adopted in the viscoelastic analytical solution of this work, which can be used to reasonably reflect the

1.72t ⎤ = (−0.0189 − 0.0135 sin θ′ + 0.0216 cos 2θ′) ⎡1 + ⎢ (75 + t) ⎥ ⎦ ⎣

(45)

10

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 14. Shape of the surface settlement and settlement curve of centre point o over time when displacement BCTC is time-dependent.

rheological behaviours of the surrounding rock by the fractional order. Based on the spring element and fractional viscous element, the Fractional Generalized Kelvin constitutive model is established. The Laplace transform and inverse transformation are adopted to obtain the relaxation modulus of the surrounding rock. A unified displacement function is adopted to capture the deformation response of the tunnel cross-section. Considering the complex variable method and unified displacement function of the tunnel crosssection, the elastic analytical solution of the ground displacement is obtained. Based on the correspondence principle, the fractional viscoelastic analytical solution is obtained by developing the obtained elastic solution. The viscoelastic analytical solution with the displacement BCTC is obtained for the first time. The displacement of the tunnel cross-section is easier to monitor than its stress. The analytical solution of this work uses this characteristic to take the displacement as the BCTC. Therefore, this work is suitable for engineering applications. Each deformation parameter of the tunnel cross-section is considered to be consistent with the same hyperbolic laws in the parameter analysis. The rheology behaviour of the ground is described by the fractional viscoelastic constitutive equation. For different ground conditions, the ground displacement can be well reflected by the fractional order. Considering a time-dependent unified displacement function and the fractional viscoelastic constitutive equation, the entire ground settlement increases,

Fig. 15. Calibration of parameters and prediction of surface settlement for two sections.

the degree of the increase becomes small, and the settlement finally stabilizes with the increase of time. For the field data of Longdong tunnel and Gaoqiao tunnel, the deformation parameters of the tunnel cross-section are first calibrated. The ground displacement field can be well predicted based on the calibrated model. The proposed method can be adopted to well evaluate the ground displacement field caused by the time effect. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement The authors gratefully acknowledge the financial support by the National Basic Research Program of China (Grant No. 2018YFC1504301), the National Natural Science Foundation of China (Grant NOs. 51421005; 51538001; 51778026) and the Natural Science 11

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

Fig. 16. Prediction of surface settlement at Day 360.

Fig. 17. Calibration of parameters and prediction of surface settlement for Gaoqiao tunnel.

Foundation of Beijing, China (Grant NO. 8161001).

Appendix A By Eqs. (14) and (16)–(20), coefficients c0 , ck and dk of two complex potential functions are expressed as,

c 0 = B0 − a 0 −

ck = B−k − bk +

1 1 a1 − b1 2 2 1 (k 2

− 1) ak − 1 −

1

(A1) 1 (k 2

+ 1) ak + 1 k = 1, 2, 3. ..

(A2)

1

dk = Bk − ak + 2 (k − 1) bk − 1 − 2 (k + 1) bk + 1 k = 1, 2, 3. ..

(A3)

where

k ⩽ −2 ⎧0 ⎪ κ (Fx − iFy) k = −1 ⎪ 4π (1 + κ ) ⎪ F − iF Bk = x y k=0 ⎨ 4π ⎪ Fx − iFy k=1 ⎪ 4π (1 + κ ) ⎪0 k⩾2 ⎩

(A4)

By the BCTC, coefficients a 0 , ak and bk of two complex potential functions are expressed as, 12

Computers and Geotechnics 117 (2020) 103284

D. Lu, et al.

α 2 (1 + κα 2k ) ak + (1 − α 2) kbk ⎫ (1 + κα 2k + 2) ak + 1 + (1 − α 2)(k + 1) bk + 1 = ⎧ k = 1, 2, 3. .. k+1 + B 2 ⎨ k + 1 − Bk α ⎬ ⎩ + Ak + 1 α ⎭

(A5)

(1 − α 2) α 2kkak − (κ + α 2k ) bk ⎫ (1 − α 2) α 2k (k + 1) ak + 1 − (κ + α 2k + 2) bk + 1 = ⎧ k = 1, 2, 3. .. ⎨+A−k α k + B−k α 2k − B−(k + 1) α 2k + 2 ⎬ ⎭ ⎩

(A6)

(1 −

α 2) a1

− (κ +

α 2) b1

= A0 − (κ + 1) a 0 + B0 −

α 2B−1

(A7)

(1 + κα 2) a1 + (1 − α 2) b1 = A1 α + (κ + 1) α 2a0 + B1 − α 2B0

(A8)

where 2

Ak =

Cg =

⎧ A′k + Fx + iFy ⎡ 1 + (α − 1) k ⎤ k ⩽ −2 2π ⎪ ⎣ k (k − 1) α k ⎦ ⎪ α (Fx + iFy ) ακ (Fx − iFy ) (2 − α 2) − 4π (1 + κ ) k = −1 ⎪ A′−1 + 4π ⎪ 2 ⎪ A′0 + Cg − α (Fx + iFy) − Fx − iFy + κ (Fx − iFy ) (1 − α 2) k = 0 ⎪ 2π 4π (1 + κ ) 4π (1 + κ ) ⎨ A′ − αC + ακ (Fx + iFy ) + Fx − iFy + g 2π 4πα (1 + κ ) ⎪ 1 ⎪ α2κ (Fx + iFy ) Fx − iFy + 4π (1 + κ ) ⎪ A′2 − 4π ⎪ ⎪ A′ − κ (Fx + iFy) ⎡ α k ⎤ ⎪ k 2π ⎣ k (k − 1) ⎦ ⎩

Fx + iFy 2π

κ log(−2ia) +

Fx + iFy 2π

α (κ − 1)(Fx − iFy )

log(2ia) +

4π (1 + κ )

k=1 k=2 k⩾3

κ (Fx + iFy ) π (1 + κ )

(A9)

log(α )

(A10)

G (α 2 − 1)2α −k − 1(ib′2 − a′2) k ⩽ −1 −G i[2u 0 α − a′1 (α 2 + 1) + 2b′2 α ] + G (1 − α 2)(b′1 + 2αa′2) k=0 Ak′ = k=1 ⎨ G i[2u 0 − a′1 (3α − α3) + 2b′2 α 2 (2 − α 2)] + G (1 − α 2)[b′1 α + 2α 2a′2] ⎪ ′1 (α 2 − 1)2α k − 2 + Gb′1 (1 − α 2)(α k − α k − 2) − G [3 + (k + 1)(α 2 − 1)](α 2 − 1)2α k − 3 [ib′2 + a′2] k ⩾ 2 G i a ⎩ ⎧ ⎪

(A11)

By the convergence of Laurent series [9], the coefficients of φ (z ) and ψ (z ) , a0, ak, bk, c0, ck, and dk are obtained.

2019;68:422–42. [14] Lu DC, Kong FC, Du XL, Shen CP, Gong QM, Li PF. A unified displacement function to analytically predict ground deformation of shallow tunnel. Tunn Undergr Space Technol 2019;88:129–43. [15] Zhang ZG, Huang MS, Xi XG, Yuan X. Complex variable solutions for soil and liner deformation due to tunneling in clays. Int J Geomech 2018;18(7). 04018074-1–19. [16] Wang HN, Chen XP, Jiang MJ, Song F, Wu L. The analytical predictions on displacement and stress around shallow tunnels subjected to surcharge loadings. Tunn Undergr Space Technol 2018;71:403–27. [17] Wang HN, Wu L, Jiang MJ. Viscoelastic ground responses around shallow tunnels considering surcharge loadings and effect of supporting. Eur J Environ Civ Eng 2018:1–23. [18] Koeller RC. Applications of fractional calculus to the theory of viscoelasticity. J Appl Mech 1984;51(2):299–307. [19] Kempfle S, Schäfer I, Beyer H. Fractional calculus via functional calculus: theory and applications. Nonlin Dyn 2002;29(1–4):99–127. [20] Alotta G, Barrera O, Cocks ACF, Di Paola M. On the behavior of a three dimensional fractional viscoelastic constitutive model. Meccanica 2017;52(9):2127–42. [21] Sumelka W, Nowak M. On a general numerical scheme for the fractional plastic flow rule. Mech Mater 2018;116:120–9. [22] Lu DC, Liang JY, Du XL, Ma C, Gao ZW. Fractional elastoplastic constitutive model for soils based on a novel 3D fractional plastic flow rule. Comput Geotech 2019;105:277–90. [23] Voyiadjis GZ, Sumelka W. Brain modelling in the framework of anisotropic hyperelasticity with time fractional damage evolution governed by the CaputoAlmeida fractional derivative. J Mech Behav Biomed Mater 2019;89:209–16. [24] Sumelka W, Voyiadjis GZ. A hyperelastic fractional damage material model with memory. Int J Solids Struct 2017;124:151–60. [25] Strack OE, Verruijt A. A complex variable solution for a deforming buoyant tunnel in a heavy elastic half-plane. Int J Num Anal Meth Geomech 2002;26(12):1235–52. [26] Li PF, Zhao Y, Zhou XJ. Displacement characteristics of high-speed railway tunnel construction in loess ground by using multi-step excavation method. Tunn Undergr Space Technol 2016;51:41–55. [27] Qu JL. Study of ground long-term settlement induced by shield construction. PHD Thesis Tongji University; 2002. [in Chinese].

References [1] Sulem J, Panet M, Guenot A. An analytical solution for time-dependent displacements in a circular tunnel. Int J Rock Mech Min Sci Geomech Abstr 1987;24(3):155–64. [2] Malan DF. Manuel Rocha medal recipient simulating the time-dependent behaviour of excavations in hard rock. Rock Mech Rock Eng 2002;35(4):225–54. [3] Sakurai S. Approximate time-dependent analysis of tunnel support structure considering progress of tunnel face. Int J Num Anal Meth Geomech 1978;2(2):159–75. [4] Hasanpour R, Chakeri H, Ozcelik Y, Denek H. Evaluation of surface settlements in the Istanbul metro in terms of analytical, numerical and direct measurements. Bull Eng Geol Environ 2012;71(3):499–510. [5] Zhang ZG, Huang MS, Xu C, Jiang YJ, Wang WD. Simplified solution for tunnel-soilpile interaction in Pasternak’s foundation model. Tunn Undergr Space Technol 2018;78(8):146–58. [6] Zhang ZG, Huang MS, Zhang CP, Jiang KM, Lu MH. Time-domain analyses for pile deformation induced by adjacent excavation considering influences of viscoelastic mechanism. Tunn Undergr Space Technol 2019;85(3):392–405. [7] Wang HN, Li Y, Ni Q, Utili S, Jiang MJ, Liu F. Analytical solutions for the construction of deeply buried circular tunnels with two liners in rheological rock. Rock Mech Rock Eng 2013;46(6):1481–98. [8] Verruijt A. Deformations of an elastic half plane with a circular cavity. Int J Solids Struct 1998;35(21):2795–804. [9] Verruijt A. A complex variable solution for a deforming circular tunnel in an elastic half-plane. Int J Num Anal Meth Geomech 1997;21(2):77–89. [10] Fu JY, Yang JS, Yan L, Abbas SM. An analytical solution for deforming twin-parallel tunnels in an elastic half plane. Int J Num Anal Meth Geomech 2015;39(5):524–38. [11] Fang Q, Song HR, Zhang DL. Complex variable analysis for stress distribution of an underwater tunnel in an elastic half plane. Int J Num Anal Meth Geomech 2015;39(16):1821–35. [12] Kong FC, Lu DC, Du XL, Shen CP. Displacement analytical prediction of shallow tunnel based on unified displacement function under slope boundary. Int J Num Anal Meth Geomech 2019;43(1):183–211. [13] Kong FC, Lu DC, Du XL, Shen CP. Elastic analytical solution of shallow tunnel owing to twin tunnelling based on a unified displacement function. Appl Math Model

13