Analytical solutions for the circular stress transducer embedded in rheological rock mass

Applied Mathematical Modelling 81 (2020) 538–558

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Analytical solutions for the circular stress transducer embedded in rheological rock mass Yuanguang Zhu a,∗, Quansheng Liu b, Xuewei Liu a, Zhanbiao Yang c a State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China b The Key Laboratory of Safety for Geotechnical and Structural Engineering of Hubei Province, School of Civil Engineering, Wuhan University, Wuhan, Hubei 430072, China c State Key Laboratory of Coking Coal Resources Development and Comprehensive Utilization, China Pingmei Shenma Group, Pingdingshan, Henan 467099, China

a r t i c l e

i n f o

Article history: Received 15 October 2019 Revised 5 December 2019 Accepted 9 January 2020 Available online 15 January 2020 Keywords: In-situ stress Rheological stress recovery Transducer Soft rock Coalmine

a b s t r a c t The rheological stress recovery (RSR) method was proposed to measure the in-situ stress of rock mass with time-dependent property by drilling a hole and embedding transducers into it. To solve the stress distribution on the transducer, a viscoelastic axisymmetric plane model was firstly built considering an arbitrary stress boundary condition. The analytical solution was developed by dividing the stress boundary conditions into axisymmetric and anti-axisymmetric combining with Laplace transformation technique. The explicit stress expressions on interfaces of rock–grout and grout–transducer was obtained using Burgers and Boltzmann viscoelastic models, respectively. Furthermore, the variations of transducer surface final stress, which is related to rheological time, geometric and mechanical properties of rock mass, grout parameter, and transducer materials, was proposed for calculating the measured stress by RSR method. For both of Burgers and Boltzmann viscoelastic model, final stress increases as elastic modulus ratio increases when elastic modulus ratio under 20, and the final stress could be ignored when diameter ratio is over 1.4. The rheological time increases with increasing of viscosity coefficient and the modulus ratio, but decreases as the shear modulus increases. The results in here provide a simple method for stress analyzing and have great value for understanding the relationship between the initial stress of rock mass and the measured stress for the RSR method. © 2020 Elsevier Inc. All rights reserved.

1. Introduction Recently, more and more tunnels and roadways will be constructed in deep area, in which surrounding rock with characteristics of soft and broken [1–4]. Therefore, it is very important for In-situ stress measurement to obtain the mechanical properties and design the rock support system [5,6]. At present, two methods, namely hydraulic fracturing and borehole relief, are the most commonly used to measure the rock stress in field [7,8]. However, these two methods generally are used in the areas with characteristics of continuous, homogeneous, isotropic and linearly elastic [9]. In fact, most of rock in deep buried areas usually are soft, broken and show a significant time-dependent behavior [10,11], which causes these two commonly used methods are hardly applied. ∗

Corresponding author. E-mail address: [email protected] (Y. Zhu).

https://doi.org/10.1016/j.apm.2020.01.034 0307-904X/© 2020 Elsevier Inc. All rights reserved.

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

539

In order to achieve the stress measurement in rheological rock masses, the rheological stress recovery (RSR) method was proposed [12–14]. The RSR method indicates that rock stress can be measured by embedding stress transducers and includes the following three steps. Firstly, drilling a testing borehole in the surrounding rock mass and embedding a six-directional compressive stress transducer at the preset depth with grout filled around it. Then, the output stresses on the transducer gradually became stable and six normal stresses could be finally recorded with different directions because of the creep property of the rock mass. Finally, the initial stress of the rock mass can be analyzed and calculated with the obtained six normal stresses. During these three steps, one of the most important problems for the RSR method is to build the relationship between the measured normal stresses and the initial stress of the rock mass. For fluid materials, generally, the measured pressure is equal to the real pressure of the surrounding medium because the calibrating and application environments of the transducer are the same. However, for rock-like or soil-like materials, it is difficult to ensure the identity of the calibrating and application environments. As a result, many attempts have to be made to study the factors affecting the stress output and present corresponding formulas. There are two methods to solve this kind of problem. The first one is numerical simulation, including finite element method [13,15], discrete element method [16], boundary element method [17], and numerical manifold method [18]. Generally, the numerical methods can solve complex problems but sometimes they are time consuming and precision limitation. The second method is analytical solutions, which can obtain the solutions of the problem directly and analyze influences of parameters on problem widely. Recently, many researches focus on analytical solutions of stress distribution on rock mass. For instance, Chen et al. [19] studied the interaction between geomaterial and liners for twin tunneling; Song et al. [20,21] investigated the influence of interface condition on lined circular tunnels in viscoelastic rock; Lin et al. [22] developed an analytical solution to estimate the stress and displacement fields around a shallow tunnel; Wang et al. [23–25] systematically studied the stresses and displacements for deep buried, shallow, and non-circular tunnels by analytical method. Duc et al. [26–28] nonlinear buckling, post-buckling, static and dynamic response of sigmoid functionally graded circular cylindrical shells based on Reddy’s third-order shear deformation shell theory. From the above references, it can be found that most of work only considers the interaction between link and rock mass and the rock mass was assumed as linearly elastic material. Interaction among three different materials and the time-dependent behaviors of rock mass were less discussed. Weiler and Kulhawy [29] indicated that the factors affecting stress output including inclusion effect, cell-soil interaction, placement effect, and environmental effect. Further, Zhu and Liu [13] found that the output stresses of the transducer are related to the initial stress of the surrounding medium and the elastic ratio between the two materials by numerical simulation. Therefore, studies on influence factors and parametric analyses for stress distribution on transducer, which is embedded in rheological rock masses, are significant and necessary. To obtain this aim, a mathematical model, which considered the borehole in the rheological surrounding rock is axisymmetric and a ring-shape transducer is embedded in the borehole, was firstly established. The rheological properties of the rock mass were evaluated within a unified viscoelastic model. A time-dependent embedding process was considered while deriving the solutions. Then, the parametric analyses were conducted on basis of the model and the influence of the rheological properties, embedding time, and diameters of the borehole on the surface stress of the transducer was discussed. 2. Problem definition Drilling of an axisymmetric borehole in a rheological rock mass was considered. As introduced by Liu et al. [14], the shape of real transducer is round sphere, with a shell outside to fix the sensor element and the inside is empty. The shell of transducer is made by SUS304 stainless steel. Therefore, for the 2-dimentional condition, the transducer can be assumed as ring shape and the transducer was embedded in the center of the borehole with grouting material filled around it. The following assumptions are made: I. The rock mass is homogeneous, isotropic, and linearly viscoelastic. The transducer and the grouting material are homogeneous, isotropic, and linearly elastic. II. The transducer is embedded deep enough to ignore the axial deformation of the borehole. The issue could be simplified as a two-dimensional plane strain problem. III. The scale of the rock mass is much larger than that of the borehole. The stress distribution in X and Y directions are uniform, with the magnitude of principal stress P and Q, respectively (Fig. 1, the tensile stress was defined as positive in this study). IV. The time of drilling the borehole and embedding the transducer is ignored. The testing process is divided into two stages. The first stage spans from t = 0 to t = t0 , with t0 being the installing time of the pressure transducer. The second stage spans from the installation of the transducer t0 onward. In rectangular coordinates (x, y, z), the stress boundary conditions for the surrounding rock are σ x = –P and σ y = –Q. Defining

p0 =

P+Q P−Q , q0 = . 2 2

(2.1)

the stress boundary condition could be divided into

σr = −p0 , τrθ = 0

(2.2)

540

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

Fig. 1. Problem description of a ring-shaped transducer embedded in the surrounding rock.

and

σr = −q0 cos 2θ , τrθ = q0 sin 2θ

(2.3)

in the cylindrical coordinate (r, θ , z). For convenient analysis, the problem of determining the stress components with boundary conditions in Eqs. (2.2) and (2.3) was defined as Problem 1 and Problem 2, respectively (Fig. 1). Based on the superposition principle, the solution of the original problem could be determined by summing the solutions of the two problems. 3. Solution of the problem According to the elasticity viscoelasticity correspondence principle, if the solution of an elastic problem is known, the Laplace transform of the solution to the corresponding viscoelastic problem can be found by replacing the elastic constants and the actual loads by their Laplace transforms. The solutions of Problems 1 and 2 were first analyzed by considering the rock mass as a linearly elastic material. Then, the solutions for the viscoelastic surrounding rock were derived using the elasticity viscoelasticity correspondence principle. The following parameters were defined: the elastic modulus and Poisson’s ratio of the surrounding rock are Er and ν r , the elastic modulus and Poisson’s ratio of the grouting material are Eg and ν g , and the elastic modulus and Poisson’s ratio of the transducer are Et and ν t . The initial radius of the borehole is r1 , the outside and inside radii of the transducer are r2 and r3 (Fig. 1). The small deformation of the borehole generated from time 0 to t1 was not included. In elastic mechanics, the two-dimensional plane strain problem included four nonzero stress components (σ r, σ θ , τ rθ , and σ z ), three nonzero strain components (ε r, ε θ , and ε rθ ), and two nonzero stress components (ur and uθ ). If volume force was ignored, the equilibrium equations can be written as follows:

∂ σ r 1 ∂ τr θ σ r − σ θ 1 ∂ σθ ∂ τr θ 2 τr θ + + = 0, + + = 0. ∂r r ∂θ r r ∂θ ∂r r

(3.1)

The kinematic equations are expressed as follows:

  ∂ ur ur 1 ∂ uθ 1 1 ∂ ur ∂ uθ uθ εr = ,ε = + , εr θ = + − . ∂r θ r r ∂θ 2 r ∂θ ∂r r

(3.2)

The constitutive equations are expressed as follows:

εr =

1+ν 1+ν 1+ν τr θ . [(1 − ν )σr − νσθ ], εθ = [(1 − ν )σθ − νσr ], εrθ = E E E

(3.3)

where E and ν are the elastic modulus and Poisson’s ratio of linearly elastic materials, respectively. Problem 1 is completely axisymmetric, and the general solution of stress components for elastic materials could be written as follows:

σr = A/r2 + B, σθ = −A/r2 + B.

(3.4)

where A and B are two functions determined by stress boundary conditions. Substituting Eq. (3.4) into Eq. (3.3) and considering Eq. (3.2), the general displacement could be obtained as follows:





A 1+ν ur = − + (1 − 2ν )Br . E r

(3.5)

Problem 2 is axisymmetric in geometry, while the boundary conditions are related to angle θ . For this problem, the general solutions of stress components for elastic materials are listed as follows:

      2C 6D 4C 6D 6D σr = − cos 2θ 2B + 2 + 4 , σθ = cos 2θ 12Ar2 + 2B + 4 , τrθ = sin 2θ 6Ar2 + 2B − 2 − 4 . r

r

r

r

r

(3.6)

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

541

where parameters A, B, C, and D are four functions determined by stress boundary conditions. Similarly, using Eqs. (3.2) and (3.3), the general displacements are obtained as follows:









1+ν 4C (1 − ν ) 2D C ( 2 − 4ν ) 1+ν 2D ur = − cos 2θ 4Aν r 3 + 2Br − − 3 , uθ = sin 2θ A(6 − 4ν )r 3 + 2Br − + 3 . E r E r r r (3.7) 3.1. Solution of Problem 1 During the first stage (0 ≤ t ≤ t0 ), the stress boundary conditions of rock mass are listed as follows:

σr (r = ∞, t ) = −p0 , σr (r1 , t ) = 0.

(3.8)

After the transducer was embedded with grout filled around it (t > t0 ), the stress boundary conditions of the rock mass, grout, and transducer can be expressed as follows:

σr (r = ∞, t ) = −p0 , σr (r1 , t ) = −p1 (t ), σr (r2 , t ) = −p2 (t ), σr (r3 , t ) = 0.

(3.9)

where p1 (t) and p2 (t) are the undetermined functions and represent the pressure on rock–grout and grout–transducer interfaces, respectively. For 0 ≤ t ≤ t0 , we have p1 (t) = p2 (t) = 0. Substituting the Eqs. (3.8) and (3.9) into Eq. (3.4), it could be determined that the stress components are deduced as follows:

σrr = Ar /r2 + Br , σθr = −Ar /r2 + Br , (rock mass ) σrg = Ag /r2 + Bg , σθg = −Ag /r2 + Bg , (grout )

σrt = At /r2 + Bt , σθt = −At /r2 + Bt .(transducer )

(3.10)

where





Ar = r12 [ p0 − p1 (t )], Br = −p0 , Ag = m2 r22 [ p1 (t ) − p2 (t )]/(m2 − 1 ), Bg = −m2 p1 (t ) + p2 (t ) /(m2 − 1 ), At = n2 r32 p2 (t )/(n2 − 1 ), Bt = −n2 p2 (t )/(n2 − 1 ), m = r1 /r2 , n = r2 /r3 .

(3.11)

In view of Eq. (3.5), the displacements of the surrounding rock, grout, and transducer are listed as follows:

1 + νr [−Ar /r + (1 − 2ν )Br r], (rock mass ) Er 1 + νg  ugr = −Ag /r + (1 − 2ν )Bg r , (grout ) Eg urr =

utr =

1 + νt [−At /r + (1 − 2ν )Bt r].(transducer ) Et

(3.12)

Consider the rock mass as a viscoelastic material, the viscoelastic constitute equation could be expressed as follows:

P  Si j = Q  ei j , P  σii = Q  εii .

(3.13)

where

P =

l

p k

k=0

1 1



dk dk dk dk , Q = q k k , P  = p k k , Q  = q k k . k dt dt dt dt

r

l

r

k=0

k=0

k=0

(3.14)

and the parameters p k , q k , p  k , and q  k are determined using the material properties. The Laplace form of Eq. (3.13) is written as follows:

P  Si j = Q  ei j , P  σ ii = Q  εii .

(3.15)

Based on the correspondence principle, the elastic modulus and Poison’s ratio of surrounding rock in Laplace space are:

E r (s ) =

3Q  Q  2P  Q  + P  Q 

, ν r (s ) =

P  Q  − P  Q  2P  Q  + P  Q 

.

(3.16)

Then, replacing the elastic parameters Er ,ν r , Ar , and Br in Eq. (3.15) by viscoelastic parameters E r , ν r , Ar , and Br , the radial displacement of the rock mass in Laplace space can be obtained:

urr

=

P Q





A 3P  Q  − r + Br r .  r 2P Q  + P  Q 

(3.17)

According to the convolution theorem, we have

L[ f (t )∗ g(t )] = f (t ) · g(t ).

(3.18)

542

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

where (∗ ) was defined as follows:



f (t ) ∗ g(t ) =

t 0

f (t − τ )g(τ )dτ .

(3.19)

In view of Eq. (3.18), Eq. (3.17) could be rewritten as

1 urr = − L[Ar ∗ M (t )] + 3rL[Br ∗ N (t )]. r by defining



M (t ) = L

−1

P Q





, N (t ) = L

−1

P  P  2P  Q  + P  Q 

(3.20)

 .

(3.21)

The inverse transform of Eq. (3.20) is listed as follows:

1 urr = − Ar ∗ M (t ) + 3rBr ∗ N (t ) r

(3.22)

Substituting Eq. (3.11) into Eq. (3.22) and considering the definition of (∗ ) in Eq. (3.19), the radial displacement solution of the rock mass could be determined as follows:

urr = −

r12 p0 r



t 0

M (t − τ )dτ +

r12 r



t 0

p1 (τ )M (t − τ )dτ − 3r p0



0

t

N (t − τ )dτ .

(3.23)

While deriving stress solutions, the following two boundary compatibility conditions were applied:

urr (r1 , t ) − urr (r1 , t1 ) = ugr (r1 , t ) for t > t0 ugr (r2 , t ) = utr (r2 , t ) for t > t0

(3.24)

Introducing Eqs. (3.23) and (3.12) into Eq. (3.24), the following equation set can be obtained:

C11 p1 (t ) − C12 p2 (t ) = D1 , C21 p1 (t ) − C22 p2 (t ) = 0. where

(3.25)

 ( 1 + νg ) 1 + ( 1 − 2 νg ) m 2 2(1 − νg2 ) 2(1 − νg2 )m2 C11 = , C12 = , C21 = , 2 2 Eg ( m − 1 ) Eg ( m − 1 ) Eg ( m2 − 1 )   (1 + νt ) 1 + (1 − 2νt )n2 ( 1 + νg ) m 2 + 1 − 2 νg C22 = + , 2 Et (n − 1 ) Eg ( m2 − 1 )   t  t  t0  t  D1 = p0 M (t − τ )dτ − M (t0 − τ )dτ − p1 (τ )M (t − τ )dτ + 3 p0 N (t − τ )dτ − 0

0

t0

0

0

t0

 N (t0 − τ )dτ . (3.26)

By solving Eq. (3.25), the functions of p1 (t) and p2 (t) are determined as follows:

p1 (t ) =

C22 C D , p2 (t ) = 21 p1 (t ) C11C22 − C12C21 1 C22

(3.27)

3.2. Solution of Problem 2 During the first stage (0 ≤ t ≤ t0 ), the stress boundary conditions of rock mass are listed as follows:

σr (r = ∞, t ) = −q0 cos 2θ , τrθ (r = ∞, t ) = q0 sin 2θ , σr (r1 , t ) = τrθ (r1 , t ) = 0.

(3.28)

For t > t0 , the known stress boundary conditions of the rock mass and transducer are:

σr (r = ∞, t ) = −q0 cos 2θ , τrθ (r = ∞, t ) = q0 sin 2θ , σr (r3 , t ) = 0, τrθ (r3 , t ) = 0.

(3.29)

and the unknown stress boundary conditions on the interfaces were defined as follows:

σr (r1 , t ) = −qc1 (t ) cos 2θ , τrθ (r1 , t ) = qs1 (t ) sin 2θ σr (r2 , t ) = −qc2 (t ) cos 2θ , τrθ (r2 , t ) = qs2 (t ) sin 2θ qc1 (t )

qs1 (t )

qc2 (t )

qs2 (t )

(3.30)

where = = = = 0 when 0 ≤ t ≤ t0 . Substituting -(Eqs. (3.28)–(3.30) into Eq. (3.6), it could be determined that the stress components of the rock mass are the following equations:

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

 σrr = − cos 2θ 2B∗r +

 2

4Cr∗ 6D + 4r r2 r



, σθr = cos 2θ 12A∗r r 2 + 2B∗r +

  2C ∗ 6 D∗ τrrθ = sin 2θ 6A∗r r2 + 2B∗r − 2r − 4 r . r

543



6D∗r , r4 (3.31)

r

The stress components of the grout are listed as follows:

    4C ∗ 6 D∗ 6 D∗ σrg = − cos 2θ 2B∗g + 2g + 4 g , σθg = cos 2θ 12A∗g r2 + 2B∗g + 4 g , r r r   ∗ ∗ 2 6 C D τrgθ = sin 2θ 6A∗g r2 + 2B∗g − 2g − 4 g . r

(3.32)

r

The stress components of the transducer are expressed as follows:

    4C ∗ 6 D∗ 6 D∗ σrt = − cos 2θ 2Bt∗ + 2t + 4t , σθt = cos 2θ 12At∗ r2 + 2Bt∗ + 4t , r r r   2Ct∗ 6Dt∗ t ∗ 2 ∗ τrθ = sin 2θ 6At r + 2Bt − 2 − 4 . r

(3.33)

r

where

A∗r = 0, B∗r = A∗g = B∗g = Cg∗ = D∗g =

r4 r2 q0 ∗ , Cr = − 1 [2q0 − qc1 (t ) − qs1 (t )], D∗r = 1 [3q0 − qc1 (t ) − 2qs1 (t )], 2 2 6



1 6r22

(

m2

3

− 1)



1 2(

m2

3

− 1)

r22



3

2 ( m2 − 1 ) r24

3

6 ( m2 − 1 )





−(m4 + 3m2 )qc1 +(m4 − 3m2 )qs1 + (3m2 + 1 )qc2 + (3m2 − 1 )qs2 ,

(m6 + m4 + 2m2 )qc1 + 2m2 qs1 − (2m4 + m2 + 1 )qc2 − 2m4 qs2 ,

−(2m6 + m4 + m2 )qc1 − (m4 + m2 )qs1 + (m6 + m4 + 2m2 )qc2 + (m6 + m4 )qs2 ,

(3m6 + m4 )qc1 + 2m4 qs1 −(m6 + 3m4 )qc2 − 2m6 qs2 ,

(n4 + 3n2 )qc2 + (−n4 + 3n2 )qs2 1 ∗ (n6 + n4 + 2n2 )qc2 + 2n2 qs2 · 2 , Bt = , 3 3 r3 6 ( n2 − 1 ) 2 ( n2 − 1 ) (2n6 + n4 + n2 )qc2 + (n4 + n2 )qs2 2 ∗ (3n6 + n4 )qc2 + 2n4 qs2 4 Ct∗ = − · r3 , Dt = · r3 . 3 3 2 ( n2 − 1 ) 6 ( n2 − 1 ) At∗ = −

(3.34)

In view of Eq. (3.7), the radial and circumferential displacements of the rock mass are listed as follows:



urr



1 + νr 4C ∗ (1 − νr ) 2D∗r =− 4A∗r νr r 3 + 2B∗r r − r − 3 cos 2θ , Er r r





C ∗ ( 2 − 4 νr ) 1 + νr ∗ 2 D∗ uθ = Ar (6 − 4νr )r 3 + 2B∗r r − r + 3 r sin 2θ . Er r r r

(3.35)

The displacements of the grout are expressed as follows:



ugr = −



4Cg∗ (1 − νg ) 2D∗g 1 + νg 4A∗g νg r 3 + 2B∗g r − − 3 cos 2θ , Eg r r



ugθ =



Cg∗ (2 − 4νg ) 2D∗g 1 + νg ∗ Ag (6 − 4νg )r 3 + 2B∗g r − + 3 sin 2θ . Eg r r

(3.36)

The displacements of the transducer are expressed as follows:



utr



4C ∗ (1 − νt ) 2Dt∗ 1 + νt =− 4At∗ νt r 3 + 2Bt∗ r − t − 3 cos 2θ , Et r r





C ∗ (2 − 4νt ) 2 D∗ 1 + νt uθ = At∗ (6 − 4νt )r 3 + 2Bt∗ r − t + 3t sin 2θ . Et r r t

(3.37)

Similarly to Problem 1, replacing the elastic parameters Er , ν r , A∗r , B∗r , Cr∗ , and D∗r in Eq. (3.35) by viscoelastic parameters E r , ν r , A∗r , B∗r , Cr∗ , and D∗r , the displacement solutions of the rock mass in Laplace space are listed as follows:

544

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

 urr

= − cos 2θ

2B∗r r

 urθ = sin 2θ

2C ∗ 2 D∗ − r − 3r r r

2B∗r r +

2D∗r r3



P



P

6C ∗ P  P  − r · ,   r Q 2P Q  + P  Q 



6Cr∗ P  P  · .  r 2P Q  + P  Q 

−

Q



The inverse transform of Eq. (3.38) is shown as follows:





(3.38)

2Cr∗ 2 D∗ 6C ∗ − 3 r ∗ M (t ) − r ∗ N (t ) , r r r

 ∗ ∗ 6 2 D C urθ = sin 2θ 2B∗r r + 3 r ∗ M (t ) − r ∗ N (t ) . r r urr = − cos 2θ

2B∗r r −

(3.39)

Substituting Eq. (3.34) into Eq. (3.39) and considering the definition of (∗ ) in Eq. (3.19), the displacement solutions of the rock mass could be determined as follows:



urr

2r 2 r4 1 + 21 − 14 r r

= −r cos 2θ +



6r12 q0 r2

t

N (t − τ )dτ −

0



uθ = r sin 2θ r

r4 1 + 14 r



6r 2 + 21 q0 r

 

t

0

 

t

q0 0

t

q0 0

3r12 r2



M (t − τ )dτ +

t

t0

N (t − τ )dτ −

t

t0

t



t0

r14 r12 − 3r 4 r2





qc1

(τ ) +



2r14 r12 − 3r 4 r2



qs1

 (τ ) M (t − τ )dτ

[qc1 (τ ) + qs1 (τ )]N (t − τ )dτ ,

M (t − τ )dτ − 





t

t0

r14 c [q (τ ) + 2qs1 (τ )]M (t − τ )dτ 3r 4 1



3r12 c [q (τ ) + qs1 (τ )]N (t − τ )dτ . r2 1

(3.40)

Four boundary compatibility conditions should be considered to determine the stress solutions:

urr (r1 , t, θ ) − urr (r1 , t0 , θ ) = ugr (r1 , t, θ ), urθ (r1 , t, θ ) − urθ (r1 , t0 , θ ) = ugθ (r1 , t, θ ),

ugr (r2 , t, θ ) = utr (r2 , t, θ ), ugθ (r2 , t, θ ) = utθ (r2 , t, θ ).

(3.41)

Introducing Eqs. (3.36), (3.37), and (3.40) into Eq. (3.41), the following equation set was obtained:

⎡

∗ ∗ ∗ ∗ C11 C12C13C14

⎤⎡

⎤

qc1 (t )

⎡ ∗⎤ D1

∗ ∗ ∗ ∗ ⎥⎢ s ⎢C21 ⎥ ⎢ ∗⎥ ⎢ C22C23C24 ⎥⎢q1 (t )⎥ = ⎢D2 ⎥ ⎣C ∗ C ∗ C ∗ C ∗ ⎦⎣qc (t )⎦ ⎣0 ⎦ 31 32 33 34 ∗ ∗ ∗ ∗ C41 C42C43C44

2 qs2

where ∗ C11 = ∗ C12 = ∗ C13 = ∗ C21 = ∗ C22 = ∗ C23

=

∗ C31 = ∗ C33 =

1−

3

Eg ( m2 − 1 )

2

1 + νg 3

Eg ( m2 − 1 )



1 + νg 3

Eg ( m2 − 1 ) 1 + νg

3

− 1)

1−

3

Eg ( m2 − 1 )



1 + νg 3

Eg ( m2 − 1 )

3



1 + νt

3

1−

3

Et (n2 − 1 ) 1 + νg

4

3

3

Eg ( m2 − 1 )

−

3

νg −





1 + νg 8 ∗ m2 + 4νg − 4 , C14 = ( 4 νg − 4 ) m 4 + 3 3 Eg ( m2 − 1 )



4 − 2 νg , 3



− 2 νg





8



4 3

νg −



4 m2 , 3

1 + νg 4 ∗ m2 − 4 + 4νg , C24 = 3 3 Eg ( m2 − 1 )

( 4 − 4 νg ) m 6 + −

3

Eg ( m2 − 1 )

8

4 − 2 νg , 3

5 2 νg m 6 − ( 3 − 2 νg ) m 4 + ( 3 − 2 νg ) m 2 + − 2 νg , 3 3

νg −

5

1 + νg

νg m 6 − 2 νg m 4 + ( 4 − 2 νg ) m 2 +

 4

1 + νg



2 5 νg m 6 + ( 5 − 6 νg ) m 4 + ( 3 − 2 νg ) m 2 + − 2 νg ) , 3 3

( 4 νg − 4 ) m 4 + 3

3

0

νg m 6 − 2 νg m 4 + ( 4 − 2 νg ) m 2 +



1 + νg Eg (

3

2

Eg ( m2 − 1 )

− ∗ C34 =



1 + νg

m2

(t )

(3.42)

 8 3

νg −



8 m2 , 3

1 + νg 8 ∗ − νg m4 + (4 − 4νg )m2 , C32 = ≥ 3 3 3 Eg ( m2 − 1 )

 4 3

−



4 νg m 4 + ( 4 − 4 νg ) m 2 , 3

2 m 6 − ( 3 − 2 νg ) m 4 − ( 5 − 6 νg ) m 2 + νg − 1 3

2 νt n6 + (5 − 6νt )n4 + (3 − 2νt )n2 + 3



5 3

− 2νt

2 − 2 νg m 6 − ( 4 − 2 νg ) m 4 + 2 νg m 2 − νg ) − 3



,

2

1 + νt 3

Et (n2 − 1 )

3

νt n6 − 2νt n4 + (4 − 2νt )n2 +

4 − 2νt , 3

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

∗ C41 = ∗ C43

=

∗ C44 =



1 + νg Eg (

3

− 1)

m2

( 4 − 4 νg ) m 6 +

4

1 + νg

−

3

Eg ( m2 − 1 ) Eg ( m2 − 1 )

Et (n2 − 1 )

D∗1 = 2q0



t 0

+ 6q0 D∗2 = 2q0

t

0



t



0

+ 6q0

0

t

 8 3

−



8 νg m 4 , 3



2 1 + νt 2 4 νg − νt n6 − 2νt n4 + (4 − 2νt )n2 + − 2νt , 3 3 3 Et (n2 − 1 ) 3









t0

0

N (t − τ )dτ −

M (t − τ )dτ −



5 2 m 6 − ( 3 − 2 νg ) m 4 + ( 3 − 2 νg ) m 2 − 1 − νg 3 3

M (t − τ )dτ −





1 + νg 4 ∗ νg m4 , C42 = 3 3 Eg ( m2 − 1 )

2 5 νt n6 − (3 − 2νt )n4 + (3 − 2νt )n2 + − 2νt . 3 3

1−

3

and





1 + νt

−

3

− 2 νg m 6 − ( 4 − 2 νg ) m 4 + 2 νg m 2 −

2 νg −

3

−

3



1 + νg

4

545



t0

0



t0

0

N (t − τ )dτ −

 M (t0 − τ )dτ N (t0 − τ )dτ

0

 −

t

t0

1

t

3

t0

 N (t0 − τ )dτ

3 

−3

 t0

2

t

t0



M (t1 − τ )dτ



 −



−3

t

t0

qc1 (τ ) +

(3.43)

1 s q (τ ) M (t − τ )dτ 3 1

[qc1 (τ ) + qs1 (τ )]N (t − τ )dτ qc1 (τ ) +

2 s q (τ ) M (t − τ )dτ 3 1

[qc1 (τ ) + qs1 (τ )]N (t − τ )dτ

(3.44)

By solving Eq. (3.42), the functions of qc1 (t ), qs1 (t ), qc2 (t ), and qs2 (t ) are determined as follows:

qc1 (t ) = C1∗ D∗1 + C2∗ D∗2 , qs1 (t ) = C3∗ D∗1 + C4∗ D∗2 ,

(C4∗C5∗ − C3∗C6∗ )qc1 + (C1∗C6∗ − C2∗C5∗ )qs1

qc2 (t ) =

C1∗C4∗ − C2∗C3∗

(C4∗C7∗ − C3∗C8∗ )qc1 + (C1∗C8∗ − C2∗C7∗ )qs1

, qs2 (t ) =

C1∗C4∗ − C2∗C3∗

.

(3.45)

where

C1∗ = C5∗ =

A∗11

A∗21

A∗12

A∗22

det Ci j

det Ci j

det Ci j

det Ci∗j

A∗13

A∗23

A∗14

A∗24

det Ci j

det Ci j

det Ci j

det Ci∗j

 ∗  , C2∗ =  ∗  , C6∗ =

 ∗  , C3∗ =  ∗  , C7∗ =

 ∗  , C4∗ =  ∗  , C8∗ =

 ,  .

(3.46)

det is determinant and Aij is the algebraic complement of Ci∗j . The solutions of stress components on the interfaces for Problems 1 and 2 can be determined using Eqs. (3.27) and (3.45). 4. Solution of the rheological models Boltzmann viscoelastic model and Burgers viscoelastic model are frequently used to model a variety of rheological characteristics of rock masses. The Boltzmann viscoelastic model can be simplified from Burgers viscoelastic model (Fig. 2) by removing Newtonian viscous dashpot η2 . For Burgers viscoelastic model,

P  (s ) = 1 +

η2

G2

s+

η1 + η2 G1

s+

η1 η2

G1 G2

s2 , P  (s ) = 1, Q  (s ) = η2 s +

η1 η2 G1

s2 , Q  (s ) = 3K.

(4.1)

Substituting Eq. (4.1) into Eq. (3.21), the following equations can be obtained:



M (t ) = L

−1

N (t ) = L

−1





G 1 G 2 + G 1 η2 s + G 2 ( η1 + η2 ) s + η1 η2 s 2 , G 1 G 2 η2 s + G 2 η1 η2 s 2



G 1 G 2 + G 1 η2 s + G 2 ( η1 + η2 ) s + η1 η2 s 2 . [6 K ( G 1 G 2 + G 1 η 2 s + G 2 η 1 s + G 2 η 2 s + η 1 η 2 s 2 ) + G 1 G 2 η 2 s + G 2 η 1 η 2 s 2 ]

(4.2)

If the surrounding rock is incompressible, that is, Kb (t) = ∞, no displacement occurred before borehole drilling. The inverse transform of Eq. (4.2) is listed as follows:

M (t ) =

1 1 1 − Gη1 t δ (t ) + + e 1 , N (t ) = 0. G2 η2 η1

(4.3)

Inserting Eq. (4.3) into Eqs. (3.26) and (3.44), the following equations can be obtained:

D1 = p0

t − t η2

0

+



G G 1 − 1t − 1t e η1 0 − e η1 G1



−

p1 (t ) − G2



t

t0

p1 ( τ )

1

1 − Gη1 (t−τ ) + e 1 dτ , η2 η1

(4.4)

546

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

Fig. 2. Burgers viscoelastic physical model.

t − t

0

t − t

0

G1 2 qc1 (t ) 1 qs1 (t ) 1 − Gη1 t0 ( e 1 − e − η1 t ) − − η2 G1 3 G2 3 G2  t  t

1

G1 2 1 1 1 − η (t−τ ) 1 − Gη1 (t−τ ) − qc1 (τ ) + e 1 dτ − qs1 (τ ) + e 1 dτ , 3 t0 η2 η1 3 t0 η2 η1

D∗1 = 2q0

+





1 qc1 (t ) 2 qs1 (t ) − η2 3 G2 3 G2   1

1

1 t c 2 t s 1 − Gη1 (t−τ ) 1 − Gη1 (t−τ ) − q1 ( τ ) + e 1 dτ − q1 ( τ ) + e 1 dτ . 3 t0 η2 η1 3 t0 η2 η1

D∗2 = 2q0

G G 1 − 1t − 1t e η1 0 − e η1 G1

+

−

(4.5)

Considering the expressions of stresses on the interfaces and substituting Eq. (4.4) into Eq. (3.27), the following expressions can be obtained:





G G C G2 t − t0 1 − 1t − 1t p1 (t ) = p0 + e η1 0 − e η1 C + G2 η2 G1

whereC = C22 /(C11C22 − C12C21 ). Defining

f (t ) =



t

t0

p1 (τ )dτ , g(t ) =



t

t0

G1

−



1

η2

t

t0

p1 ( τ )d τ −

1

η1

G

e

− η1 t



t

1

t0

p1 ( τ )e

G1

η1

τ



dτ .

τ

p 1 ( τ ) e η1 d τ .

Eq. (4.6) could be rewritten as follows:

f  (t ) =







G G C G2 t − t0 1 − 1t − 1t p0 + e η1 0 − e η1 C + G2 η2 G1

(4.6)

(4.7)



−

1

η2

f (t ) −

1

η1

G

e

− η1 t 1



g(t ) .

(4.8)

Taking the derivative of both sides of Eq. (4.8) yielded

f  (t ) =

 1

1  1 − Gη1 t  1 − Gη1 t G1 1 − Gη1 t + e 1 − f  (t ) − e 1 g (t ) + e 1 g(t ) . η2 η1 η2 η1 η1 η1

C G2 p0 C + G2

(4.9)

In view of Eqs. (4.7) and (4.8), the expressions of g(t) and g (t) could be obtained, and substituting them into Eq. (4.9) yielded

f  (t ) +

CG  1  G 1

G1 1 G1 C G2 C G2 G1 1 2 1 + + f  (t ) + f (t ) = p0 + (t − t0 ) + e− η1 t0 . C + G 2 η1 η2 η1 η1 η2 C + G 2 C + G 2 η2 η1 η2 η1

(4.10)

Eq. (4.10) is a second-order linear nonhomogeneous differential equation and its solution is following equation:

f (t ) = m1 ex1 (t−t0 ) + m2 ex2 (t−t0 ) + p0 (t − t0 ) − p0 η2 where



x1,2

1 C G2 =− 2 C + G2

G1

G1

( 1 − e − η1 t 0 ) +

1 1 + . G2 C

(4.11)

1  G 1  C G  1  G 2 1 1 G1 G2 C 1 2 1 + + ± + + −4 . η2 η1 η1 2 C + G 2 η2 η1 η1 η1 η2 C + G 2

and

m1 = −

1





(4.12)

G1 G1 p0 p 0 η2 x 2 1 1 1 p0 p 0 η2 x 1 1 1 1 − ( 1 − e − η1 t 0 ) + + , m2 = + ( 1 − e − η1 t 0 ) + + . ( x1 − x2 ) ( x1 − x2 ) G1 G2 C ( x1 − x2 ) ( x1 − x2 ) G1 G2 C (4.13)

In view of Eq. (4.7), the expression of p1 (t) could be obtained as follows:

p1 (t ) = p0 + x1 m1 ex1 (t−t0 ) + x2 m2 ex2 (t−t0 ) .

(4.14)

For Boltzmann model, no η2 term existed, and Eq. (4.10) could be rewritten as follows:

p1 (t ) +

 CG 2

C + G2

·

1

η1

+

G1

η1



p1 (t ) = p0

C G2 1 − Gη1 t0 e 1 . C + G 2 η1

(4.15)

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

547

The solution of Eq. (4.15) is listed as follows: G

p1 (t ) =

p0 e

− η 1 t0 1



1+

1 − ex3 (t−t0 ) G1 G2

+

.

G1 C

(4.16)

where

x3 = −

C G2

(C + G2 )η1

−

G1

η1

.

(4.17)

Further, substituting Eq. (4.5) into Eq. (3.45) yielded



2C ∗ + C2∗ 3+ 1 G2



f  c (t ) +





G G C1∗ + 2C2∗  t − t0 1 − 1t − 1t f s (t ) = 6q0 (C1∗ + C2∗ ) + e η1 0 − e η1 G2 η2 G1



1

1

fc (t ) + K1 − (C1∗ + 2C2∗ ) fs (t ) + K2 , η2 η2   t − t  G1  G 2C3∗ + C4∗  C3∗ + 2C4∗ 1 − t − 1t 0 f c (t ) + 3 + f  s (t ) = 6q0 (C3∗ + C4∗ ) + e η1 0 − e η1 G2 G2 η2 G1 1

1

− (2C3∗ + C4∗ ) f (t ) + K1 − (C3∗ + 2C4∗ ) f (t ) + K2 . η2 c η2 s − (2C1∗ + C2∗ )

where fc (t), fs (t), K1 , and K2 were defined as follows:

fc (t ) =



t

t0

qc1 (τ )dτ , fs (t ) =



t

t0

qs1 (τ )dτ , K1 =

1

η1

G

e

− η1 t



t

1

t0

G1

τ

qc1 (τ )e η1 dτ , K2 =

According to Eq. (4.18), the expressions of K1 and K2 are obtained as follows

K1 = 2q0 K2 = 2q0

t − t

0

t − t

0

η2 η2



G G 1 − 1t − 1t + e η1 0 − e η1 G1

+



G G 1 − 1t − 1t e η1 0 − e η1 G1

 

− −

1

η2 1

η2



1

η1

G

e

− η1 t



fc (t ) −

C3∗ + 2C4∗ 1 + C1∗C4∗ − C2∗C3∗ G2

fs (t ) +

2C3∗ + C4∗  f (t ) − C1∗C4∗ − C2∗C3∗ c

f  c (t ) +



(4.18)



t

1

t0

G1

τ

qs1 (τ )e η1 dτ .

(4.19)

C1∗ + 2C2∗  f (t ), C1∗C4∗ − C2∗C3∗ s

2C1∗ + C2∗ 1 + C1∗C4∗ − C2∗C3∗ G2



f  s (t ).

(4.20)

Taking the derivative of both sides of Eq. (4.18) and inserting the expression of K1 and K2 into Eq. (4.19) yielded

(a1 D2 + a2 D + a3 ) fc (t ) + (a4 D2 + a5 D + a6 ) fs (t ) = 6q0 (C1∗ + C2∗ )



1

1 η2 + η1 e

(a7 D2 + a8 D + a9 ) fc (t ) + (a10 D2 + a11 D + a12 ) fs (t ) = 6q0 (C3∗ + C4∗ ) where D =

d dt



G

− η 1 t0 1

0 + Gη11 t−t η2 ,

G

1 1 1 − η1 t 0 0 + Gη11 t−t η2 + η1 e η2 .

(4.21)

and







2C1∗ + C2∗ 1 1 G1 G1 , a2 = (2C1∗ + C2∗ ) + + +3 , G2 η1 η2 η1 G 2 η1 C ∗ + 2C2∗ G1 a3 = (2C1∗ + C2∗ ), a4 = 1 , η1 η2 G2 1  2C ∗ + C4∗ 1 G1 G1 a5 = (C1∗ + 2C2∗ ) + + , a6 = (C1∗ + 2C2∗ ), a7 = 3 , η1 η2 η1 G 2 η1 η2 G2 1  ∗ C + 2C4∗ 1 G1 G1 a8 = (2C3∗ + C4∗ ) + + , a9 = (2C3∗ + C4∗ ), a10 = 3 + 3 , η1 η2 η1 G 2 η1 η2 G2 1 

1 G1 G1 G1 a11 = (C3∗ + 2C4∗ ) + + +3 , a12 = (C ∗ + 2C4∗ ). η1 η2 η1 G 2 η1 η1 η2 3 a1 = 3 +

(4.22)

Eq. (4.21) is a second-order linear nonhomogeneous differential equation set denoted as follows:

a = a1 a10 − a4 a7 , b = a1 a11 + a2 a10 − a4 a8 − a5 a7 , c = a1 a12 + a2 a11 + a3 a10 − a4 a9 − a5 a8 − a6 a7 , d = a2 a12 + a3 a11 − a5 a9 − a6 a8 , e = a3 a12 − a6 a9 .

(4.23)

The general solutions of Eq. (4.21) could be determined using Cramer’s rule as follows

fc (t ) =

4

i=1

Mi eXi (t−t0 ) + A1 t + A2 , fs (t ) =

4

i=1

αi Mi eXi (t−t0 ) + A3t + A4 .

(4.24)

548

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

where

X1,2,3,4

b 1 =− ± 4a 2

=



 



√ 3 21



3 a 3 2 +

3

8d − ab3 + 4abc 2c 1  b2 4c b2 2 − a − +  ± − −  − ,   3a 2 2a2 3a 4a2 2 4 4ba2 − 32ac + 



3

−431 + 22

+

2 +



−431 + 22 , 1 = c2 − 3bd + 12ae, √ 3 3 2a a1 Xi2 + a2 Xi + a3

2 = 2c3 − 9bcd + 27ad2 + 27b2 e − 72ace, αi = −

a4 Xi2 + a5 Xi + a6

( i = 1, 2, 3, 4 ).

(4.25)

and

A1 = A3 = 2q0 , G

A 2 = 2 q 0 η2

A 4 = 2 q 0 η2

e

− η 1 t0 1

−1

G1 e

G1

− η t0 1

G1

−1

−

−

!

t0

C ∗ + 2C ∗ − C3∗ − 2C4∗ 1 − + 1 ∗ 2∗ G2 C1 C4 − C2∗C3∗

t0

2C3 ∗ + C4∗ − 2C1∗ − C2∗ 1 − + G2 C1∗C4∗ − C2∗C3∗

η2 η2

! (4.26)

Substituting Eq. (4.24) into Eq. (4.19) yielded 4

Mi Xi eXi (t−t0 ) , qs1 (t ) = 2q0 +

4

αi Mi Xi eXi (t−t0 ) .

(4.27)

fc (t1 ) = 0, fc (t1 ) = qc1 (t0 ) = 0, fs (t1 ) = 0, fs (t1 ) = qs1 (t0 ) = 0.

(4.28)

qc1 (t ) = 2q0 +

i=1

i=1

Considering the initial conditions,

Inserting Eqs. (4.24) and (4.27) into initial conditions of Eq. (4.28) yielded the expressions of Mi are following:

⎡

⎤⎡

⎤

⎡

⎤

M1 −A1 t1 − A2 ⎢α1 α2 α3 α4 ⎥⎢M2 ⎥ ⎢−A3t1 − A4 ⎥ ⎣X X X X ⎦⎣M ⎦ = ⎣ − A ⎦. 1 2 3 4 3 1 α1 X1 α2 X2 α3 X3 α4 X4 M4 − A3 1111

(4.29)

The solution of Eq. (4.29) is expressed as follows:

Mi = − where

A 1i ( A1 t1 + A2 ) + A 2i ( A3 t1 + A4 ) + A 3i A1 + A 4i A3 . det(Hi j )

(4.30)

" " "1111 " " " "α1 α2 α3 α4 " Hi j = " ". "X1 X2 X3 X4 " "α1 X1 α2 X2 α3 X3 α4 X4 "

(4.31)

For Boltzmann viscoelastic model, a3 = a6 = a9 = a12 = 0. Then, Eq. (4.21) could be rewritten as follows:

(a1 D + a2 )qc1 (t ) + (a4 D + a5 )qs1 (t ) = 6q0 (C1∗ + C2∗ )

1

G

η1

(a7 D + a8 )qc1 (t ) + (a10 D + a11 )qs1 (t ) = 6q0 (C3∗ + C4∗ )

e

− η 1 t0

1

η1

1

G

e

,

− η 1 t0 1

.

(4.32)

The solutions of Eq. (4.32) are listed as follows:

qc1 (t ) = M5 eX5 (t−t0 ) + M6 eX6 (t−t0 ) + 2q0 A5 , qs1 (t ) = M7 eX5 (t−t0 ) + M8 eX6 (t−t0 ) + 2q0 A6 . where

√ √ b2 − 4ac −b − b2 − 4ac X5 = , X6 = , 2a 2a  

G1 3 − η t0 1 G1 G1 A5 = e 1 (C1∗C4∗ − C2∗C3∗ ) + + 3(C1∗ + C2∗ ) , c η1 η1 η1 G 2 η1 1 

3 − Gη1 t0 G1 G1 A6 = e 1 (C1∗C4∗ − C2∗C3∗ ) + + 3(C3∗ + C4∗ ) , c η1 η1 η1 G 2 η1 −b +

(4.33)

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558



549



(a1 X6 + a2 )(a4 X5 + a5 ) (a4 X6 + a5 )(a4 X5 + a5 ) + A6 , (a2 a4 − a1 a5 )(X6 − X5 ) (a2 a4 − a1 a5 )(X6 − X5 )

M5 = 2 q0 A5

M6 = −2q0 A5 − M5 , M7 = −

a1 x1 + a2 a1 x2 + a2 M5 , M8 = − M6 . a4 x1 + a5 a4 x2 + a5

(4.34)

According to Eqs. (3.9) and (3.30), the stress components on the rock–grout and grout–transducer contact surfaces are listed as follows:

σr (r1 , t ) = −p1 (t ) − qc1 cos 2θ , τrθ (r1 , t ) = qs1 sin 2θ , σr (r2 , t ) = −p2 (t ) − qc2 cos 2θ , τrθ (r2 , t ) = qs2 sin 2θ .

(4.35)

Substituting Eqs. (4.14), (4.27), (3.27), and (3.45) into Eq. (4.35), the stress components on the rock–grout and grout– transducer contact surfaces for Burgers model could be finally determined as follows:

σr ( r 1 , t ) = − p 0 +

2

!

mi xi e

xi (t−t0 )

i=1

τr θ ( r 1 , t ) = 2 q 0 +

4

− 2q0 +

4

αi Mi Xi e

Mi Xi e

cos 2θ ,

i=1

! Xi (t−t0 )

!

Xi (t−t0 )

sin 2θ .

(4.36)

i=1

# C σr (r2 , t ) = − 21 C22

p0 +

2

! mi xi exi (t−t0 )

i=1

4

$

!

i=1

C1∗C6∗ − C2∗C5∗ 2q0 + αi Mi Xi eXi (t−t0 ) cos 2θ , ∗ ∗ ∗ ∗ C1 C4 − C2 C3 4

+

!

C ∗C ∗ − C3∗C6∗ + 4∗ 5∗ 2q0 + Mi Xi eXi (t−t0 ) cos 2θ C1 C4 − C2∗C3∗

# τr θ ( r 2 , t ) =

i=1

C4∗C7∗ C1∗C4∗

− C3∗C8∗ − C2∗C3∗

2q0 +

4

!

Mi Xi e

Xi (t−t0 )

i=1

C ∗C ∗ − C2∗C7∗ + 1∗ 8∗ 2q0 + αi Mi Xi eXi (t−t0 ) C1 C4 − C2∗C3∗ 4

!$ sin 2θ .

(4.37)

i=1

Similarly, the stress components on the rock–grout and grout–transducer contact surfaces for Boltzmann model could be determined as follows:



G

σr ( r 1 , t ) = − 

p0 e 1+

τr θ ( r 1 , t ) = M 7 e

− η 1 t0

G1 G2

1

+

X5 (t−t0 )

G1 C

1−e

−



C G2 C+G2

−

τr θ ( r 2 , t ) =

t−t







− M5 eX5 (t−t0 ) + M6 eX6 (t−t0 ) + A5 cos 2θ ,



+ M8 eX6 (t−t0 ) + A6 sin 2θ . G1

t



C p e η1 0 − σr (r2 , t ) = − 21 · 0 G1 G1 1 − e C22 1 + + C G2 −



+G1 · η 0 1



C G2 C+G2



t−t

+G1 · η 0 1

(4.38)

 −

C4∗C5∗ − C3∗C6∗  M5 eX5 (t−t0 ) + M6 eX6 (t−t0 ) + A5 cos 2θ C1∗C4∗ − C2∗C3∗

C1∗C6∗ − C2∗C5∗  M7 eX5 (t−t0 ) + M8 eX6 (t−t0 ) + A6 cos 2θ , C1∗C4∗ − C2∗C3∗

C4∗C7∗ − C3∗C8∗  C ∗C ∗ − C2∗C7∗  M5 eX5 (t−t0 ) + M6 eX6 (t−t0 ) + A5 sin 2θ + 1∗ 8∗ M7 eX5 (t−t0 ) + M8 eX6 (t−t0 ) + A6 sin 2θ . ∗ ∗ ∗ ∗ ∗ ∗ C1 C4 − C2 C3 C1 C4 − C2 C3 (4.39)

5. Results and discussion When the RSR method is used for measuring rock stress, the rheological time and the final stress on the transducer surface are two of the most important factors that need attention. In theory, the parameters xi and Xi should be negative and −xi −1 and −Xi −1 actually represent the rheological time when the stresses on the transducer surface become stable. For convenience, −xi −1 and −Xi −1 are expressed as −x−1 = τi and −Xi−1 = ξi , respectively. The expression of xi and Xi shown in i Eqs. (4.12), (4.17), (4.25), and (4.34) shows that the rheological time parameters τ i and ξ i are not only dependent on the mechanical parameters of rock masses, grout, and transducer but also related to the geometrical parameters of the borehole and transducer (Table 1). Based on Eqs. (4.37) and (4.39) and considering the observation time t → ∞, the final stress on the transducer surface could be simplified as follows:

σr (r2 , ∞ ) = −p0Cp − 2q0Cqc cos 2θ , τrθ (r2 , ∞ ) = 2q0Cqs sin 2θ

(5.1)

550

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558 Table 1 Factors affecting the rheological time and final stress on the transducer surface. Rock mass model

Burgers model

Typical parameters

C p , Cqc , Cqs

Influence factors

Eg , νg , Et , νt , m, n.

Boltzmann model

τ 1 ,τ 2 Eg , νg , Et , νt , G1 , G2 , η1 , η2 , m, n.

ξ 1 ,ξ 2 ,ξ 3 ,ξ 4 Eg , νg , Et , νt , G1 , G2 , η1 , η2 ,

C p∗ , Cqc∗ , Cqs∗ Eg , νg , Et , νt , G1 , G2 , η1 , m, n, t0 .

m, n.

τ3

Eg , νg , Et , νt , G1 , G2 , η1 , m, n.

ξ 5 ,ξ 6 Eg , νg , Et , νt , G1 , G2 , η1 , m, n.

Table 2 Magnitude of parameters for Burgers model. Parameters

Cp

Cqc

Cqs

τ i / hour

ξ i / hour

Magnitude

1.062

0.595

1.590

τ 1 = 349.1, τ 2 = 32.01

ξ 1 = 33.01, ξ 2 = 35.28, ξ 3 = 352.9, ξ 4 = 365.8

Table 3 Magnitude of parameters for Boltzmann model. Parameters

C p∗

Cqc∗

Cqs∗

τ 3 / hour

ξ i / hour

Magnitude

0.597

0.333

0.877

τ 3 = 41.28

ξ 5 = 39.21, ξ 6 = 38.22

σr (r2 , ∞ ) = −p0Cp∗ − 2q0Cqc∗ cos 2θ , τrθ (r2 , ∞ ) = 2q0Cqs∗ sin 2θ

(5.2)

where −

G1

t

Cp =

C21 ∗ C21 C G 2 e η1 0 , Cp = · , C22 C G2 + G1 (C + G2 ) C22

Cqc =

C4∗C5∗ − C3∗C6∗ + C1∗C6∗ − C2∗C5∗ s C ∗C ∗ − C3∗C8∗ + C1∗C8∗ − C2∗C7∗ , Cq = 4 7 , ∗ ∗ ∗ ∗ C1 C4 − C2 C3 C1∗C4∗ − C2∗C3∗

 Cqc∗ =



C4∗C5∗ − C3∗C6∗ C ∗C ∗ − C2∗C5∗ A5 + 1∗ 6∗ A6 , Cqs∗ = ∗ ∗ ∗ ∗ C1 C4 − C2 C3 C1 C4 − C2∗C3∗





C4∗C7∗ − C3∗C8∗ C ∗C ∗ − C2∗C7∗ A5 + 1∗ 8∗ A6 . ∗ ∗ ∗ ∗ C1 C4 − C2 C3 C1 C4 − C2∗C3∗

(5.3)

Eqs. (5.1) and (5.2) show the final stress on the transducer surface, which depends on parameters Cp , Cqc , and Cqs for Burgers model and parameters C p∗ , Cqc∗ , and Cqs∗ for Boltzmann model. According to Eq. (5.3) parameters Cp , Cqc , Cqs , C p∗ , Cqc∗ , and Cqs∗ are defined as the ratio of several parameters and therefore, these parameters have no units. Besides, From the Eqs. (5.1) and (5.2), it also can be found the final normal stress increases as the values of parameters Cp , Cqc , C p∗ , and Cqc∗ . The final shear stress increases as the values of parameters Cqs and Cqs∗ . Considering the expression of the parameters shown in Eq. (5.3), the final stress on the transducer surface for Burgers model is further found to be related to the elastic parameters of the grout and transducer and the diameters of the borehole and transducer but independent of the rheological parameters of the rock masses. In comparison, the final stress on the transducer for Boltzmann model is more complicated because it is also related to the rheological parameters and the embedding time of the transducer. Using Eq. (5.3), the elastic parameters of the grout and transducer materials are easily obtained. The diameter parameters of the borehole and transducer are usually predetermined before rock stress measurement. Thus, the most difficult theoretical problem for the RSR method is to determine the rheological parameters of the rock masses. However, the results showed that the final stress on the transducer surface was independent of rheological parameters for Burges model, which is usually used to represent the rock masses with significant rheological behavior. The next problem is how to determine the parameters for the rock masses, which satisfied the Boltzmann model. An example was presented herein to illustrate the influence factors on the resulting stress distribution. Considering the experimental results, the following values were adopted for the Burgers model: G1 = 2.58 GPa, G2 = 1.90 GPa, η1 = 317 GPa·h, η2 = 169 GPa·h, and ν r = 0.35. For Boltzmann model, η2 =∞ while other parameters remained the same as for Burgers model. For the elastic parameters of the grout and transducer, Gg = 20 GPa, ν g = 0.35, Gt = 120 GPa, and ν t = 0.25 were assumed. The in situ stress at infinity of P = 15 MPa and Q = 10 MPa were assumed. The radius of the borehole was r1 = 55 mm, and the outside and inside radii of the transducer were r2 = 40 mm and r3 = 20 mm, respectively. The parameters defined in Eq. (5.3) were calculated, as shown in Table 2 (Burgers model) and Table 3 (Boltzmann model). The parameters presented in here were used to investigate the influence of material mechanical properties and transducer and borehole diameters on the final stress and rheological time.

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

551

Fig. 3. Variations in rheological time with rheological parameters for Burgers model.

5.1. Influence of the rheological properties Four rheological parameters, including two shear moduli G1 and G2 and two viscosity coefficients η1 and η2 , were evaluated for Burgers model. The variations in rheological time ξ i (i = 1, 2, 3, and 4) and τ i (i = 1 and 2) with these four parameters are shown in Fig. 3. In Fig. 3, the letters L and R indicate the left and right y-direction coordinate axis, respectively. For each figure, the three lines with letter L correspond to the coordinate axis on the left and the other three lines with letter R correspond to the coordinate axis on the right. As shown in Fig. 3(a) and (b), when the shear moduli G1 increases from 1.0 GPa to 2.4 GPa, rheological time ξ 1 , ξ 2 , and τ 1 decrease from about 40 h to about 32 h and rheological time ξ 3 , ξ 4 , and τ 2 decrease from about 550 h to about 300 h. Similarly, when the shear moduli G2 increases from 2.0 GPa to 3.4 GPa, rheological time ξ 1 , ξ 2 , and τ 1 decrease from about 40 h to about 27 h and rheological time ξ 3 , ξ 4 , and τ 2 decrease from about 380 h to about 330 h. Meanwhile, as shown in Fig. 3(c) and (d), when viscosity coefficient η1 increases from 100GPa·h to 240GPa·h, rheological time ξ 1 , ξ 2 , and τ 1 increase from about 22 h to about 45 h and rheological time ξ 3 , ξ 4 , and τ 2 increase from about 330 h to about 380 h. Besides, when viscosity coefficient η2 increases from 200GPa·h to 340GPa·h, rheological time ξ 1 , ξ 2 , and τ 1 increase from about 30 h to about 33 h and rheological time ξ 3 , ξ 4 , and τ 2 increase from about 250 h to about 370 h. Therefore, it is clearly from Fig. 3 that the rheological time decreases with increasing of shear modulus and increases as the viscosity coefficient increases. From the point of applications of the RSR method, the observation time of the transducer is about 1 month, which is acceptable for engineering requirements. Three rheological parameters, including two shear moduli G1 and G2 and a viscosity coefficient η1 , were evaluated for Boltzmann model. According to Eq. (5.2), the final stress increases with increasing of values of parameters C p∗ , Cqc∗ , and Cqs∗ . Therefore, as shown in Fig. 4(a), all the parameters decrease as G1 increases, which indicates that the final stress on the transducer surface decreases with increasing of G1. Similarly, the final stress increases as G2 increases (Fig. 4b).

552

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

Fig. 4. Influence of rheological parameters on final stress on the transducer surface.

Compared with G1 and G2 , the values of parameters C p∗ , Cqc∗ , and Cqs∗ nearly keep constant when η1 increases from 100GPa·h to 240GPa·h. Therefore, the influence of η1 on final stress was negligible (Fig. 4c). Similar to Burgers model, the increasing in shear moduli G1 and G2 reduced the rheological time while the viscosity coefficient η1 had an opposite effect (Fig. 5). Compared with Fig. 3 and Fig. 5, it can be found that the rheological time for Burgers model is about 300 h for different parameters, while this value decreased to about 40 h for Boltzmann model. The result in here indicates that the rheological time for Boltzmann model is less about 10 days than Burgers model. 5.2. Influence of the mechanical properties and solidification time of the grout Eqs. (3.26) and (3.46) show the elastic parameters of the transducer and grout included in Cij and Ci∗ , respectively. The influence of mechanical properties of the transducer and grout on final stress and rheological time could be reflected via the ratio of Et to Eg , that is, Et /Eg , as shown in Eq. (5.3) With increasing of Et /Eg , for Burgers model, the absolute value of normal stress σ r increases and the shear stress τ rθ decreases (Fig. 6). Besides, the changing of σ r and τ rθ became extremely small when the value of Et /Eg over 20. Meanwhile, the rheological time increases as the Et /Eg increases (Fig. 7). In contrast, the performance of Boltzmann model had the same characteristics (Figs. 8 and 9). The installing time of the transducer affected the final stress on the transducer surface for Boltzmann model (Table 1). From Fig. 10, it can be found that with the grout solidification time increases from 0 to 300 h, the values of parameters C p∗ ,

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

553

Fig. 5. Variations in rheological time with rheological parameters for Boltzmann model.

Cqc∗ , and Cqs∗ decreases quickly, which means final normal and shear stresses decrease with the increasing of grout solidification time.

5.3. Influence of borehole diameter Parameter m is the ratio of the borehole diameter and the external diameter of the transducer. The influences of parameter m on final stress and rheological time for Burgers model and Boltzmann model are given in Figs. 11–14, respectively. From Fig. 11, for Burgers model, it can be found that the values of parameter Cqc , and Cqs first increase from 0.5 and 1.45 to 0.6 and 1.6, respectively as the m increases from 1.25 to 1.4, and then these two values keep constant when m overs 1.4. Similarly, as shown in Fig. 13, when the value of parameter m is larger than 1.4, the values of parameter Cqc∗ , and Cqs∗ stop increasing and remain unchanged for Boltzmann model. According to Eqs. (5.1) and (5.2), the final stress increases as values of parameter Cqc , Cqs , Cqc∗ , and Cqs∗ increase, therefore, the result means that for different rheological models, the increasing in final stress is very small when the magnitude of m exceeded 1.40. However, from Figs. 12 and 14, it is easily found there is not evident influence of m on the rheological time for both Burgers and Boltzmann models. As a result, when ratio of the borehole diameter and the external diameter of the transducer is over 1.4, the final stress on the transducer has no change for different rheological rock.

554

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

Fig. 6. Influence of elastic ratio on final stress on the transducer surface for Burgers model.

Fig. 7. Variations in rheological time with elastic ratio for Boltzmann model.

Fig. 8. Influence of elastic ratio on final stress on the transducer surface for Boltzmann model.

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

Fig. 9. Variations in rheological time with elastic ratio for Boltzmann model.

Fig. 10. Influence solidification time on final stress on the transducer surface for Boltzmann model.

Fig. 11. Influence of borehole diameters on final stress on the transducer surface for Burgers model.

555

556

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

Fig. 12. Variations in rheological time with the borehole diameters for Burgers model.

Fig. 13. Influence of borehole diameters on final stress on the transducer surface for Boltzmann model.

Fig. 14. Variations in rheological time with the borehole diameters for Boltzmann model.

5.4. Validation and application The proposed viscoelastic axisymmetric plane model was validated in here by comparing the results from this model and analytical solutions (Song et al. [20]) and experimental results (Liu et al. [14]). Song et al. [20] proposed analytical stress and displacement solutions for lined circular tunnels in viscoelastic rock. In this model, there are two elastic liners, namely the first liner and second liner, which is corresponding to the grout and ring-shaped transducer in here. The influence of

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

557

diameter ratio of the first liner to the second liner on final stress was discussed and the results shown that normalized equivalent stresses of the second liner increases nonlinearly with increasing of diameter ratio. The results obtained by Song et al. [20] is consistent with the results proposed by the Figs. 11 and 13. Besides, Liu et al. [14] investigated the relationship between grout setting time and rate of stress recovery for the transducer. Seven different grout setting time, including 0 h, 10 h, 20 h, 50 h, 100 h, 200 h, and 500 h, were set and stress recovery value at different time were obtained. The results showed that stress recovery decreases as grout setting time increases. Generally, a larger stress recovery means a larger final measure stress on the transducer. Therefore, the normal and shear stresses decrease with the increasing of grout solidification or setting time, which agrees well with the analytical stress obtained by the proposed model (showed in Fig. 10). When using RSR method measuring in-situ stress in practice, the first step is drilling a testing borehole in the surrounding rock mass and then put transducer in the borehole with grout filled around it. Therefore, the borehole diameter and grout setting time is very important for the final measured stress on transducer. If the borehole diameter is too large, it will be time consuming and raise the construction cost. However, the transducer will be hard to put into the borehole if the diameter is too small. Therefore, a suitable and reasonable borehole diameter can improve measure accuracy and reduce cost. Theorical study by the proposed model shows that when ratio of the borehole diameter and the external diameter of the transducer is over 1.4, the final stress on the transducer has no change for different rheological rock. Generally, the transducer diameter is constant and then the borehole diameter is suggested as 1.4 times of transducer diameter for different rock in practice. Further, analytical stress solution indicates that the longer grout setting time, the smaller the final normal and shear stress of the transducer eventually will be. Therefore, after the transducer has been put into the borehole, the quick-setting grout materials should be grouted around the transducer immediately considering that the stress increases quickly at the first several hours. 6. Conclusions The stress distribution of the transducer embedded in the rheological rock mass was considered as a geometric axisymmetric plane problem based on the procedures of the RSR method. An arbitrary initial stress field was assumed, with the rock masses as linearly viscoelastic and the grout and transducer materials as purely elastic. The analytical solutions for the plane problems were put forward by dividing the initial stress into two symmetric stress boundaries and using the Laplace transform. The explicit expressions were given for the stress distribution on the rock–grout and grout–transducer interfaces. The solutions of the equations were obtained by assuming either Burgers or Boltzmann viscoelastic model for the rock mass so that the explicit analytical expressions for the time-dependent stress distributions in the interfaces could be derived. The expressions of final stress on the transducer surface and the rheological time for the stress to be stable were two of the most important concerns for the application of the RSR method. The factors affecting the final stress and rheological time were discussed, including viscoelastic parameters of the rock mass, elastic parameters of the transducer and grout, embedding time of the transducer, and diameters of boreholes. For Burgers viscoelastic model, which is commonly used to describe a rock mass with strong rheological behavior, the final stress on the transducer surface was independent of the viscoelastic parameters and the embedding time of the transducer. The absolute value of final stress increases positively with increasing of elastic modulus ratio Et /Eg and gradually becomes stable when Et /Eg exceeded 20. The influence of borehole diameter on final stress could be ignored when borehole diameter is 1.4 times larger than that of the outer diameter of the transducer. The rheological time increases positively with increasing of viscosity coefficient and the modulus ratio Et /Eg , but the shear modulus has an opposite effect. For the Boltzmann viscoelastic model, usually employed for the rock mass, exhibited limited viscosity. It showed similar influences on the final stress and rheological time as Burgers model except that the final stress was affected by the viscoelastic parameters and the embedding time of the transducer. The findings in the present study might be used to build the relationship between measured stress and initial stress for measuring rock stress using the RSR method. Acknowledgments The work was supported by the National Natural Science Foundation of China [grant numbers 51874275] and Natural Science Foundation of Hubei Province (grant numbers 2019CFB535). References [1] H.T. Yu, Z.W. Zhang, J.T. Chen, A. Bobet, M. Zhao, Y. Yuan, Analytical solution for longitudinal seismic response of tunnel liners with sharp stiffness transition, Tunn. Undergr. Space Technol. 77 (2018) 103–114. [2] X. Liu, Q. Fang, D.L. Zhang, Mechanical responses of existing tunnel due to new tunnelling below without clearance, Tunn. Undergr. Space Technol. 80 (2018) 44–52. [3] H. Zhou, C.Q. Zhang, Z. Li, D.W. Hu, J. Hou, Analysis of mechanical behavior of soft rocks and stability control in deep tunnels, J. Rock Mech. Geotechnol. Eng. 6 (2014) 219–226. [4] S.Q. Yang, M. Chen, H.W. Jing, K.F. Chen, B. Meng, A case study on large deformation failure mechanism of deep soft rock roadway in Xin’an coal mine, China. Eng. Geol. 217 (2017) 89–101. [5] X. Liu, Q. Liu, B. Liu, Y. Zhu, P. Zhang, Failure behavior for rocklike material with cross crack under biaxial compression, J. Mater. Civ. Eng. 31 (2019) 06018025.

558

Y. Zhu, Q. Liu and X. Liu et al. / Applied Mathematical Modelling 81 (2020) 538–558

[6] X. Liu, Q. Liu, L. Wei, X. Huang, Improved strength criterion and numerical manifold method for fracture initiation and propagation, Int. J. Geomech. 17 (2017) E4016007. [7] B.C. Haimson, F.H. Cornet, ISRM suggested methods for rock stress estimation—Part 3: hydraulic fracturing (HF) and/or hydraulic testing of pre-existing fractures (HTPF), Int. J. Rock Mech. Min. Sci. 40 (2003) 1011–1020. [8] X.R. GE, M.X. HOU, A new 3D in situ rock stress measuring method: borehole wall stress relief method (BWSRM) and development of geostress measuring instrument based on BWSRM and its primary applications to engineering, Chin. J. Rock Mech. Eng. 30 (2011) 2161–2180. [9] C. Ljunggren, Y. Chang, T. Janson, R. Christiansson, An overview of rock stress measurement methods, Int. J. Rock Mech. Min. Sci. 40 (2003) 975–989. [10] M.C. He, H.P. Xie, S.P. Peng, Study on rock mechanics in deep mining engineering, Chin. J. Rock Mech. Eng. 24 (2005) 2803–2813. [11] J. Sun, S. Wang, Rock mechanics and rock engineering in China: developments and current state-of-the-art, Int. J. Rock Mech. Min. Sci. 37 (20 0 0) 447–465. [12] F. Zhang, Q. Liu, C. Zhang, Measurement of geostress and sensor about rheological stress recovery method, Rock Soil Mech. 35 (2014) 3273–3279. [13] Y. Zhu, Q. Liu, A new three-dimensional pressure transducer for measuring soft rock stress, Geotech. Test. J. 39 (2016) 703–711. [14] B. Liu, Y. Zhu, Q. Liu, X. Liu, A novel in situ stress monitoring technique for fracture rock mass and its application in deep coal mines, Appl. Sci. 9 (2019) 3742. [15] Y.Y. Wang, J. Ji, J.D. Jiang, Simulation of stress recovery rule of three-dimensional pressure transducer embedded in rheological rock, Sci. Technol. Eng. 16 (2016) 82–87. [16] M.J. Jiang, Z.F. Shen, S. Utili, DEM modeling of cantilever retaining excavations: implications for lunar constructions, Eng. Comput. 33 (2016) 366–394. [17] M.Y. Fattah, Boundary element analysis of a lined tunnel problem, Int. J. Eng. 25 (2012) 89–97. [18] W. Wei, Q. Jiang, A modified numerical manifold method for simulation of finite deformation problem, Appl. Math. Model. 48 (2017) 673–687. [19] F. Chen, L. Lin, D. Li, Analytic solutions for twin tunneling at great depth considering liner installation and mutual interaction between geomaterial and liners, Appl. Math. Model. 73 (2019) 412–441. [20] F. Song, H. Wang, M. Jiang, Analytical solutions for lined circular tunnels in viscoelastic rock considering various interface conditions, Appl. Math. Model. 55 (2018) 109–130. [21] F. Song, H. Wang, M. Jiang, Analytically-based simplified formulas for circular tunnels with two liners in viscoelastic rock under anisotropic initial stresses, Constr. Build. Mater. 175 (2018) 746–767. [22] L. Lin, F. Chen, Z. Huang, An analytical solution for sectional estimation of stress and displacement fields around a shallow tunnel, Appl. Math. Model. 69 (2019) 181–200. [23] H.N. Wang, G.S. Zeng, S. Utili, M.J. Jiang, L. Wu, Analytical solutions of stresses and displacements for deeply buried twin tunnels in viscoelastic rock, Int. J. Rock Mech. Min. Sci. 93 (2017) 13–29. [24] H.N. Wang, X.P. Chen, M.J. Jiang, F. Song, L. Wu, The analytical predictions on displacement and stress around shallow tunnels subjected to surcharge loadings, Tunn. Undergr. Space Technol. 71 (2018) 403–427. [25] H.N. Wang, G.S. Zeng, M.J. Jiang, Analytical stress and displacement around non-circular tunnels in semi-infinite ground, Appl. Math. Model. 63 (2018) 303–328. [26] N.D. Duc, S.E. Kim, D.T. Manh, P.D. Nguyen, Effect of eccentrically oblique stiffeners and temperature on the nonlinear static and dynamic response of S-FGM cylindrical panels, Thin Wall. Struct. 146 (2020) 106438. [27] N.D. Khoa, H.T. Thiem, N.D. Duc, Nonlinear buckling and postbuckling of imperfect piezoelectric S-FGM circular cylindrical shells with metal-ceramic-metal layers in thermal environment using Reddy’s third-order shear deformation shell theory, Mech. Adv. Mater. Struct. 26 (3) (2019) 248–259. [28] N.D. Duc, H. Hadavinia, T.Q. Quan, N.D. Khoa, Free vibration and nonlinear dynamic response of imperfect nanocomposite fg-cntrc double curved shallow shells in thermal environment, Eur. J. Mech. A/Solid 75 (2019) 355–366. [29] W.A. Weiler Jr, F.H. Kulhawy, Factors affecting stress cell measurements in soil, J. Geo. Eng. Div. 108 (1982) 1529–1548.