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

## 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.else...

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 ﬁrstly 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 ﬁnal 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, ﬁnal stress increases as elastic modulus ratio increases when elastic modulus ratio under 20, and the ﬁnal stress could be ignored when diameter ratio is over 1.4. The rheological time increases with increasing of viscosity coeﬃcient 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 ﬁeld [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 signiﬁcant time-dependent behavior [10,11], which causes these two commonly used methods are hardly applied. ∗

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

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 ﬁlled around it. Then, the output stresses on the transducer gradually became stable and six normal stresses could be ﬁnally 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 ﬂuid 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 diﬃcult 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 ﬁrst one is numerical simulation, including ﬁnite 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 inﬂuences 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 inﬂuence of interface condition on lined circular tunnels in viscoelastic rock; Lin et al. [22] developed an analytical solution to estimate the stress and displacement ﬁelds 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 inﬂuence factors and parametric analyses for stress distribution on transducer, which is embedded in rheological rock masses, are signiﬁcant 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 ﬁrstly established. The rheological properties of the rock mass were evaluated within a uniﬁed 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 inﬂuence of the rheological properties, embedding time, and diameters of the borehole on the surface stress of the transducer was discussed. 2. Problem deﬁnition 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 ﬁx 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 ﬁlled 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 simpliﬁed 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 deﬁned 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 ﬁrst 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. Deﬁning

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 deﬁned 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 ﬁrst 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 deﬁned: 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 ﬁrst 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 ﬁlled 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 deﬁned 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 deﬁning



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 deﬁnition 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 ﬁrst 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 deﬁned 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 deﬁnition 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 simpliﬁed 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 ). Deﬁning

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:

p1 (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 deﬁned 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 2 1

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

−4 31 + 22

+

2 +

−4 31 + 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 ﬁnally 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 ﬁnal 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 ﬁnal stress on the transducer surface could be simpliﬁed 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 ﬁnal stress on the transducer surface. Rock mass model

Burgers model

Typical parameters

C p , Cqc , Cqs

Inﬂuence 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 ﬁnal 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 deﬁned 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 ﬁnal normal stress increases as the values of parameters Cp , Cqc , C p∗ , and Cqc∗ . The ﬁnal shear stress increases as the values of parameters Cqs and Cqs∗ . Considering the expression of the parameters shown in Eq. (5.3), the ﬁnal 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 ﬁnal 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 diﬃcult theoretical problem for the RSR method is to determine the rheological parameters of the rock masses. However, the results showed that the ﬁnal stress on the transducer surface was independent of rheological parameters for Burges model, which is usually used to represent the rock masses with signiﬁcant rheological behavior. The next problem is how to determine the parameters for the rock masses, which satisﬁed the Boltzmann model. An example was presented herein to illustrate the inﬂuence 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 inﬁnity 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 deﬁned 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 inﬂuence of material mechanical properties and transducer and borehole diameters on the ﬁnal 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.

552

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

Fig. 4. Inﬂuence of rheological parameters on ﬁnal 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 inﬂuence of η1 on ﬁnal stress was negligible (Fig. 4c). Similar to Burgers model, the increasing in shear moduli G1 and G2 reduced the rheological time while the viscosity coeﬃcient η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. Inﬂuence of the mechanical properties and solidiﬁcation 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 inﬂuence of mechanical properties of the transducer and grout on ﬁnal stress and rheological time could be reﬂected 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 ﬁnal stress on the transducer surface for Boltzmann model (Table 1). From Fig. 10, it can be found that with the grout solidiﬁcation 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 ﬁnal normal and shear stresses decrease with the increasing of grout solidiﬁcation time.

5.3. Inﬂuence of borehole diameter Parameter m is the ratio of the borehole diameter and the external diameter of the transducer. The inﬂuences of parameter m on ﬁnal 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 ﬁrst 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 ﬁnal stress increases as values of parameter Cqc , Cqs , Cqc∗ , and Cqs∗ increase, therefore, the result means that for different rheological models, the increasing in ﬁnal 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 inﬂuence 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 ﬁnal 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. Inﬂuence of elastic ratio on ﬁnal stress on the transducer surface for Burgers model.

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

Fig. 8. Inﬂuence of elastic ratio on ﬁnal 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. Inﬂuence solidiﬁcation time on ﬁnal stress on the transducer surface for Boltzmann model.

Fig. 11. Inﬂuence of borehole diameters on ﬁnal 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. Inﬂuence of borehole diameters on ﬁnal 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 ﬁrst liner and second liner, which is corresponding to the grout and ring-shaped transducer in here. The inﬂuence of

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

557

diameter ratio of the ﬁrst liner to the second liner on ﬁnal 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 ﬁnal measure stress on the transducer. Therefore, the normal and shear stresses decrease with the increasing of grout solidiﬁcation 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 ﬁrst step is drilling a testing borehole in the surrounding rock mass and then put transducer in the borehole with grout ﬁlled around it. Therefore, the borehole diameter and grout setting time is very important for the ﬁnal 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 ﬁnal 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 ﬁnal 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 ﬁrst 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 ﬁeld 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 ﬁnal 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 ﬁnal 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 ﬁnal stress on the transducer surface was independent of the viscoelastic parameters and the embedding time of the transducer. The absolute value of ﬁnal stress increases positively with increasing of elastic modulus ratio Et /Eg and gradually becomes stable when Et /Eg exceeded 20. The inﬂuence of borehole diameter on ﬁnal 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 coeﬃcient 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 inﬂuences on the ﬁnal stress and rheological time as Burgers model except that the ﬁnal stress was affected by the viscoelastic parameters and the embedding time of the transducer. The ﬁndings 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 modiﬁed numerical manifold method for simulation of ﬁnite 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 simpliﬁed 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 ﬁelds 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-inﬁnite 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.