Consolidation analysis for layered transversely isotropic viscoelastic media with compressible constituents due to tangential circular loads

Computers and Geotechnics 117 (2020) 103257

Contents lists available at ScienceDirect

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

Research Paper

Consolidation analysis for layered transversely isotropic viscoelastic media with compressible constituents due to tangential circular loads

T

Zhi Yong Ai , Zi Ye, Yong Zhi Zhao ⁎

Department of Geotechnical Engineering, Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, College of Civil Engineering, Tongji University, 1239 Siping Road, Shanghai 200092, PR China

ARTICLE INFO

ABSTRACT

Keywords: Extended precise integration method Consolidation Transverse isotropy Viscoelasticity Compressible constituents Tangential circular loads

In soft soil area, the tangential load induced by the vehicle brake, tidal action and lateral loaded pile is inevitable for the consolidation problem. This article investigates the consolidation of a multilayered transversely isotropic viscoelastic media with compressible constituents subjected to tangential circular loads. A general solution scheme is established, in which the integration transform is firstly used to convert the governing equations into the ordinary differential equations, and then the extended precise integration method is employed to solve the aforementioned equations. The Merchant model is introduced to simulate the viscoelasticity of medium with the aid of the elastic-viscoelastic correspondence principle. After the verification implemented by existing analytical solutions and ABAQUS, numerical examples are given to further reveal the effects of viscoelastic parameters, compressibility of soil particles and the stratification on the consolidation behavior.

1. Introduction In geotechnical engineering, porous materials consisting of soil skeleton and voids, such as rocks and soils, are widely investigated in the fluid-solid coupling analysis [1], heat transfer analysis [2] and wave propagation analysis [3]. Particularly, consolidation problem is the most enduring and common topic due to its wide application. Since Biot [4–6] firstly established rigorous consolidation theory considering the coupling of pore water pressure and displacement, numerous researchers [7–15] have followed his work to investigate this complicated and interesting problem based on the assumption that geomaterial is regarded as a homogeneous half space or a finite soil layer. As is well-known, natural geotechnical materials show obvious stratification owing to long-time complex deposition, which significantly affects the mechanical behavior [16–19]. Therefore, it is more reasonable to take the horizontal stratification of the geotechnical materials into account. Booker and Small [20–22] utilized the finite layer method to investigate the consolidation problem of multilayered media. Meanwhile, there are also some researchers investigating this problem by the boundary element method [23,24] and the transfer matrix method [25]. Subsequently, Ai et al. [26] presented the analytical layer-element method that only possesses decaying exponential functions in the stiffness matrix of layer-element to deal with the unstable condition in transfer matrix.

⁎

Nevertheless, the majority of the aforementioned research aimed at the axisymmetric consolidation problem, which cannot accurately simulate the general three-dimensional condition that not only includes the vertical load, but also the tangential load. In soft soil engineering, the tangential load induced by the vehicle brake, tidal action and lateral loaded piles is inevitable in the analysis of consolidation problem. Hence, in order to simulate more generalized consolidation problems, it is necessary to enrich the research mentioned above to the general three-dimensional and non-axisymmetric conditions. Vardoulakis and Harnpattanapanich [27,28] investigated the distribution of displacement along depth caused by a vertical load in the Cartesian coordinate system. Senjuntichai and Rajapakse [29] derived the exact stiffness matrix for 3-D quasi-static response of poroelastic medium in the cylindrical coordinate system. Pan [30] derived the available Green’s functions when a dislocation source is acting in a layered and poroelastic system. Wang and Fang [31] employed the transfer matrix to obtain the state space solution of non-axisymmetric consolidation problem and discussed the influence of three boundary conditions. Mei et al. [32] as well as Chiou and Chi [33] investigated the consolidation problem in three-dimensional condition. Later, Ai and his cooperators [34–36] continued to make an effort on the non-axisymmetric consolidation problem. It is worth mentioning that most foregoing works assume that the deformation of soil particle is extremely minuscule compared with the

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

https://doi.org/10.1016/j.compgeo.2019.103257 Received 28 May 2019; Received in revised form 5 August 2019; Accepted 11 September 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

deformation of soils, in other words, the compressibility of soil particle is negligible and the settlement of consolidation is completely subjected to the compression of soil pore. In fact, in terms of soft soils and soft rocks, it is inevitable for the deformation of fluid and soil particle caused by their own compression [37]. Rice and Cleary [38] discussed the stress diffusion of soils subjected to a vertical load considering the compressibility of constituents. Chen [39], Singh et al. [40] as well as Ai and Hu [41] investigated the effect of the anisotropic permeability, the compressibility of fluid as well as soil particle. However, the compressibility of soil particle, in general, shows obvious transversely isotropy due to a long-time deposition, which has not been considered so far. This is one of the reasons that we conduct the present research. In the above works, the settlement of soils terminates when the excess pore pressure caused by load fully converts into the effective stress, and the secondary consolidation induced by the viscosity of soils is not considered. It should be pointed that for geotechnical material, especially soft clay, as a typical viscoelastic medium, its rheology determines that it will face long-term settlement problem and the secondary consolidation is non-ignorable. Therefore, it is meaningful to consider the viscosity of soft soil in the study of consolidation problem. The early meaningful work was done by Taylor and Merchant [42] who firstly provided typical Kelvin model to simulate the deformation of soil skeleton, and subsequently this work is extended to multi-element model by Gibson and Lo [43] as well as Lo [44]. Biot [6] also presented a three-dimensional viscoelastic medium model, nevertheless, the complexity of its mathematical solution limits its application. Based on the Voigt and Maxwell viscoelastic model, Singh and Rosenman [45] found that the displacements induced by vertical dip-slip fault for the viscoelastic case were consistent with that for the elastic case. Cai et al. [46] and Qin et al. [47] discussed the one-dimensional consolidation of saturated and unsaturated viscoelastic soils, respectively. Xie et al. [48] presented a four-element rheological model, and addressed the influence of dimensionless parameters subjected to different loading types. Recently, Ai et al. [49] analyzed viscoelastic consolidation behavior of transversely isotropic layered soils under 3-D and axisymmetric condition. It should be pointed that due to the involved complexity, the multi-dimensional consolidation solution of viscoelastic material is seldom found, and the consolidation problem for transversely isotropic viscoelastic layered soils due to tangential circular loads is not studied so far. Thus, this paper will fill this research gap and investigate on the response of viscoelastic multilayered media with a more generalized transversely isotropic characteristic under tangential load. To deal with the complicated governing equations, the integration transform and extended precise integration method [50] are well combined to solve them. Compared with traditional transfer matrix method [25,31,34], there is no any negative exponential function in the precise integration method, which is unconditionally stable during the calculation process. Meanwhile, compared with the analytical layer-element method [26,35,36,41], since the dimensions and values of relational matrices used in the calculation process are invariable, the whole theory is significantly easy to be implemented by numerical platforms (e.g. MATLAB), especially when the analytical layer-element is hard to obtain for the problems associated with material anisotropy. Thereafter, verification examples are provided and several numerical examples are conducted to study the influence of the viscoelastic parameters in the Merchant model on the consolidation behavior of the medium. In order to clearly address the breakthrough and novelty, Ref. [49] is presented as an example here to be compared with this work: (1) The first novelty of this work, as compared with Ref. [49], is that a more generalized transversely isotropy soil model is presented, and this model can take the transversely isotropy of the compressibility of soil particle into account, which ensures our research item to be more realistic and comprehensive. Meanwhile, the compressibility of fluid in Ref. [49] is reflected with a simplified parameter whilst it is strictly derived and has a clear expression in this work. (2) The non-axisymmetric consolidation problem is solved in this paper with the aid of Fourier expansion and

Hankel transform, whereas Ref. [49] mainly discussed the three-dimensional and axisymmetric conditions based on the double-Fourier transform and Hankel transform, respectively. (3) This paper mainly discusses the viscoelastic consolidation subjected to the tangential loads whilst Ref. [49] focuses on the vertical ones. (4) Moreover, more systematic and comprehensive parameters analysis is given in this paper, especially for the viscoelastic parameters and the compressibility of soil particles. 2. Ordinary differential equation matrix for the medium studied Ignoring the body force in multilayered medium, we can establish the following equilibriums equations in the cylindrical coordinate system: r

r r

r rz

r

1 r

+

+

1 r

+

1 r

r

+

+

z

+

rz

r

+

z z

z z

z

=0

r

2

+

r

=0

r rz

+

(1a) (1b)

=0

r

(1c)

where r , , z represent the normal stresses in the r , , z directions, respectively; while r , rz , z denote the shear stresses in the r - , r - z, - z planes, respectively. Combining with the principle of effective stress law, the constitutive equations are defined as follows:

s

c11 c12 = c 13

c12 c13 c11 c13 03× 3 c13 c33 03 × 3 c3 × 3

s

(2)

where s = [ r , , z , r , rz , z ]T is the total stress vector and = [ h , h , v , 0, 0, 0]T is the excess pore water pressure vector;

the strain vector is besides, c11 =

u ur ur , + , r r r

s=

2 vh ) ,

uz , z

ur + z

uz uz , + r r

u u , r + z r

u r

u r

T

;

2 c12 = ( h + vh ) , c13 = vh (1 + h ) , c33 = c44 0 0 2 0 = (1 c3 × 3 = 0 c44 c44 = Gv , , h ), 0 0 (c11 c12 ) 2 2 E v [(1 + h )(1 2 vh )], = Eh E v ; Gv , E v and Eh are the shear h modulus, vertical Young’s modulus and horizontal Young’s modulus, respectively; h and vh are the Poisson's ratio of the horizontal strain caused by the horizontal stress and the Poisson's ratio of the horizontal strain caused by the vertical stress, respectively. Furthermore, (c11 + c12 + c13) 3Ks and v =1 (2c13 + c33) 3Ks are the Bioth = 1 Willis coefficients in the horizontal and vertical directions, respectively. Here Ks is the bulk modulus of soil particle [51]. Considering the transversely isotropy of seepage and compressibility, we can establish the seepage constitutive equation as follow [52]:

(1

ur u 1u + r + r r r

h

=

1 w

kh

2

r2

+

+

1 r r

z

uz +M z t

+ kz

2

z2

(3)

M = [(1 n) Ks + n Kf (2c11 + c33 + 2c12 + 4c13) where represents the Biot’s modulus [51], here n is the porosity and Kf is the bulk modulus of fluid. When Kf = Ks = , the constituents of soil are assumed to be incompressible; v denotes the volume strain; kh and k z are the permeability coefficients in the horizontal and vertical directions, respectively. It is assumed that the seepage of pore fluid obeys the Darcy’s Law, then the total fluid flow in the z direction Qz with respect to the time from 0 to t can be obtained: 9Ks2] 1

2

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

kv

t

Qz =

0

z

w

dt

(6e), (6f). Then, in order to eliminate variable , the application of Fourier expansion with regard to coordinate leads to:

(4)

Given that these partial differential equations mentioned above are hard to solve, the Laplace and Hankel transform are applied to obtain variables s and in lieu of variables t and r, which converts the partial differential equations into ordinary differential equations. The expression of Laplace transform with respect to t and its inversion transform are as follows [53]:

f (r , z, s ) =

1 2 i

f (r , z , t ) =

st dt

f (r , z, t )e

0

c+i

f

( , z, s ) =

(5b)

f (r , z , s ) Jm ( r ) r dr

0

f (r , z, s ) =

m

f

0

u cos m , m = 0 hm

uz =

m=0

T=

u 1 = z c44

z

uz 1 = ( z c33

z

(6a)

1 uz r

(6b)

+

c13 c33

)

v

ur u 1 u + r + r r r

=

rz

z

=

c33

2 c13

+

c11 c33

r2

1 2 r r

c11 + c12 2

c33

(

2

1 r r

+

1 r2

3c11

+

)

c11 2 c13

c12 2

c12 1 2 r2 1 r2

c33

2 2

u +

c13 c33

ur

(

h

c13 c33

z

r

)

b3 = Fz =

z

2 c13

=

2 c13

+

c11 c33 c33

2

r2

1 2 r r

c11 + c12 2

c33

(

1 r r

+

1 r2

2 c13

+

)

c11 2

3c11

c33

c12

c12 1 r2

1 r2

2

2 2

u +

c13 1 c33 r

ur

(

h

c13 c33

z

v

=

cos m m = 0 zm m=0

Qzm cos m

(8)

(9)

0

b2 b2 0 0 0 0 b1 b1

2 0 0

b4 2 0 0 0

b4 2 0 0

b3 b3 0 0 0 0 b1

(b5 + b6 ) b7 2 0 0 0 0 b2 b3

2

b7 2 (b5 + b6 ) 0 0 0 0 b2 b3

c13 c33 v sw kv

,

h,

=

2 v

c33

2 c11 c33 c13 1 , b5 = , c33 c44 1 kh 2 + s . M w

b4 =

+

m+1 vm ,

m 1 hm ,

z

=

Qz z

=

r

(

c13 c33

h

2 v

c33

1 r

+

+

1 M

1 r

rz

v

)(

kh sw

(

ur r 2

r2

+

z

ur r

+

+

c11

c12 2

r

1 r

1 r r

(6g) 1u r r

+

1 r2

)+ 2 2

v

c33 z

+

)

(6h)

Subsequently, for the sake of simplification, we assume:

uv = v

=

ur + u u u , uh = r 2 2 rz

+

r

,

h

=

rz

(7a) r

m zm,

0 0 0 F 0 0 0 0

m T m]

33

b6 =

(6f)

z

0 0 0 0 2 2 0 0

2

, b7 =

2 c13

c33

33

c11 + c12 2

With the aid of the series and parallel connections of spring and dashpot, three typical rheological model, namely the Maxwell model, the Kelvin model and the Merchant model, are established to represent different rheology characteristics. The diagrams of the three models are shown as Fig. 1. The constitutive equations of the Maxwell model, the Kelvin model and the Merchant model are:

+

)

z

3. Elastic-viscoelastic correspondence principle

(6e) z

cos m m = 0 hm

Wm (z ) = [Vm, Um Vm = [ where and m T m +1 m 1 m Um = [u vm , uhm , uzm , Q zm ] ; the subscript and superscript m indicate the mth term of Fourier expansion and mth-order Hankel transforms, respectively. When m = 1, the aforementioned matrix is used to denote the tangential circular loads whilst the case of vertical circular loads c can be obtained via m = 0 Furthermore, b1 = c v , b2 = c13 ,

+ v

h

Qz =

vm cos m

(10)

(6d)

2 c13

m=0

=

]T ,

(6c)

s = w Qz z kv

0 0 2

(5d)

uz r

=

H1 H2 H3 H4 0 0

where f ( , z, s ) is the corresponding function of f (r , z, s ) in the Hankel transformed domain; is the Hankel transform parameter with respect to r; Jm ( r ) represents the mth –order Bessel function. Firstly, with the aid of the Laplace transform for Eqs. (1)–(4), we can derive the following equations: rz

,

v

d Wm (z ) = T ·Wm (z ) dz

m

ur 1 = z c44

u zm cos m ,

cos m m=0 m

(5c)

( , z , s ) Jm ( r ) d

u vm cos m ,

Substituting the Eq. (8) into Eqs. (6)–(7), mth term of Fourier expansion processed by mth-order Hankel transform with regard to r can be expressed as an ordinary differential equation matrix:

where f (r , z, s ) represents the corresponding function after f (r , z , t ) is applied by Laplace transform; s is the Laplace transform parameter with 1. respect to t; i = On the basis of Laplace transform, the expression of Hankel transform and its inversion transform are as follows [53]: m

m=0

=

(5a)

f (r , z , s )e st ds

c i

uv = uh =

(7b)

It is noted that the partial derivation of the four variables in Eq. (7) can be obtained by the combination of Eqs. (6a), (6b) as well as Eqs.

Fig. 1. The viscoelastic model. 3

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

Fig. 2. The combination of adjacent layer elements.

+

0 d = E 0 dt

= E1 + +

1

0

d dt

=

(11a)

d dt

d E0 E1 1 E0 d = + E0 + E1 dt E0 + E1 E0 + E1 dt

1 1 + E0 s 0

(11c)

1 E1 + s

By Fourier expansion, Laplace and Hankel transform, the ordinary differential equation matrix established in the Section 2 forms the relationship of eight variables and their partial derivatives with respect to z . In this section, the precise integration method is extended to solve this matrix, and then the real solution can be obtained by the Laplace and Hankel inversion. For a multilayered medium with total depth L , we divide the whole medium into 2N micro layers with length l = L 2 N , named as layer elements, as shown in Fig. 2. Here, in order to satisfy the accuracy of the calculation, N are supposed to be set as a sufficiently large value such as 15–30. Next, we will establish the relationship between the relational matrix of interlevel state vector and the ordinary differential equation matrix summarized above. Here for any arbitrary adjacent layer elements, the diagram of state vectors is shown as follows: In terms of the above two adjacent layer elements with the depth of [z a , z b] and [z b , z c ], we define them as layer element 1 and 2, respectively. The inherent relationship between the generalized state vectors (V and U ) at the top and bottom of a layer element, corresponding to Wm (z ) = [Vm, Um]T in Eq. (9), can be established by the relational matrix [56] as follows: For layer element 1:

G1Ub

(14a)

Ua = Q1Va + F1T Ub

(14b)

Vb = F1Va (12b)

1

(13b)

4.1. The establishment and combination of layer elements

For the Merchant model:

1 1 V (s ) = + E0 E1 + s

v

4. Extended precise integration solution for layered media

(12a)

1

1 1 + E0v E1v + s

where is the proportional coefficient between the horizontal and G vertical direction and = Ev . v Subsequently, substituting the reciprocal of these flexibility coefficients into the ordinary differential equation matrix in the transformed domain, we can obtain the viscoelastic solution for this problem.

For the Kelvin model:

V (s ) =

(13a)

v1

Vg (s ) = Vh (s ) = Vv (s ) =

where E0 is the Young’s modulus of spring; 0 represents the viscosity coefficient of the Maxwell model. E1 and 1 denote the Young’s modulus and the viscosity coefficient of the Kelvin model, respectively. Based on the series of a linear elastic spring described by Hooke's law and a linear viscous dashpot described by Newton's law, it is convenient for the Merchant model to simulate the creep and stress relaxation behaviors of soft soils synthetically. It is interesting to note that with t tending to 0, the stress reaches a stable value and strain reaches 0, which reflects that this model can clearly describe both the timelimited stress relaxation and creep. Thus, the Merchant model is applied in this study to describe the viscous behavior of soil skeleton. For the above equations, Lee [54] proposed the elastic-viscoelastic correspondence principle, in which he found that the governing equations of elastic theory is similar to that of viscoelastic theory in the Laplace domain, which means that the solution process of viscoelastic theory can be simplified by replacing the corresponding viscoelastic modulus with elastic modulus in the transformed domain. Kovarik [55] provided that this principle is applicable for the quasi-static linear elastic problem that the boundary conditions are independent of time, such as the consolidation. For the isotropic medium, under the assumption that the Poisson's ratio remains unchanged, the stress-strain relationship is completely determined by the Young’s modulus E . From this, we can use the reciprocal of flexibility coefficient V (s ) = (s ) (s ) of viscoelastic medium in the Laplace domain to establish constitutive equations instead of elastic modulus. Then we have: For the Maxwell model:

V (s ) =

h1

For the Merchant model, we have

(11b)

1

Eh0 E = h1 = Ev0 E v1

Similarly, for layer element 2:

G2 Uc

(15a)

Ub = Q2 Vb + F2T Uc

(15b)

Vc = F2Vb

(12c)

Nevertheless, the above derivation is suitable for the isotropic medium rather than the transversely isotropic medium. That is to say, for the transversely isotropic medium, there are three flexibility coefficients, Vh (s ) = h (s ) h (s ) , Vv (s ) = v (s ) v (s ) and Vg (s ) = g (s ) g (s ) , corresponding to Eh , E v and Gv , respectively. Here, with the assumption that parameters Eh and E v as well as h and v in horizontal and vertical directions keeps the same proportion relationship, a simplification is put forward and applied in our work, i.e.

where Fi , Ei and Qi (i = 1, 2) represent 4 × 4 dimensional relational matrices. Under the condition that the thickness of layer element l is small enough, the three relational matrices can be expressed as a Taylor expansion (16) with the matrices H1, H2 , H3 and H4 mentioned in Eq. (10) and the thickness of layer element l . Meanwhile, based on the premise of calculation accuracy, we cut off the high-order terms of the Taylor expansion to accelerate the calculation. 4

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

F (x ) = I + F (x ), 2

+

3l

2

+

3l

Q (x )

1l

+

2l

G (x )

1l

+

2l

3

+

3

+

(16a)

f1 l + f2 l2 + f3 l3 + f4 l 4

F (x )

4l 4l

(16b)

4

(16c)

4

The matrices fi , i , i (i = 1, 2, 3, 4) can be expressed by the matrices H1, H2 , H3 and H4 :

H1 f1 + 1 H3 , 2 H1 f3 + 3 H3 + 2 H3 f1 +

f1 = H1, =

1

=

H3,

=

2

H3 f2 +

H1 f2 +

2 H3

+

1 H3 f1

3

,

f4

1 H3 f2

(17a)

=

H3 f1 + f1T H3 2

f2T H3

,

3

f1T H3 f1

+ , 4 3 T T T H3 f3 + f3 H3 + f1 H3 f1 + f1 H3 f2

= =

f3 =

4

=

1

f2 =

(17b)

4 H2, H1

3

2

+

=

1 H3 2

1 H4

+ H1 1 , 3 = 2 + 2 H3 1 + 3 H4 4

2 H4

+ H1

2

3

+

1 H3 1

,

4

Fig. 3. The diagram of layer block division.

(17c)

Vc = F3Va

G3 Uc

Ua =Q3 Va +F3T Uc where,

G3 = G2 + F2 (G1 1 + Q2) 1F2T

(19b)

Q3 = Q1 + E1 (Q2 1 + G1) 1F1

(19c)

^ V +F ^ TU Ua = Q 1 a 1 b

(20b)

^ V + Uc = Q 3 c

^ TU F 3 d

(22a) (22b)

^ Q ^ 1 ^ Vc = (I + G 12 3) (F12 Va ^ V + Uc = Q 3 c

^ F ^T G 12 3 Ud

^P R 1 V

S^1PU )

^ TU F 3 d

U+b

(23a) (23b)

^ represents the where the subscript ‘12’ of relational matrices like F 12 merge of layer block 1 and layer block 2, and they can be expressed as:

For a multilayered medium to be solved, the calculation plane Hc and load plane interface Hl divide the whole medium into three zones, which are named as layer block 1, layer block 2 and layer block 3, separately (see Fig. 3). Through the method in the last section, the state vectors of each layer block surface can be obtained by the merge of layer elements. It is noteworthy that this problem is supposed to be considered as two conditions, i.e., Hc > Hl and Hc < Hl . Firstly, we discuss the condition Hc > Hl , and we can establish the relationship as follows. It can be seen that interface c is the calculation plane in this condition. For layer block [a , b]: (20a)

^ U G 3 d

U b and where the superscript “ ” and “+” in V b , represent the upper and lower surfaces of load plane, in which the inherent rePV and U+b = U b PU can be established with lationship V +b = V b the generalized load stress vector PV and generalized load displacement vector PU . By eliminating the state vector of interface b, the coalescence of Eqs. (20)–(22) leads to:

4.2. The merge of layer block with applied load

^ U G 1 b

^V Vd = F 3 c

V +b ,

There are 2N layer elements separated in the whole medium. We can continuously apply the same merge operation to every adjacent layer element; that is to say, the number of layer elements is halved after each operation step. Furthermore, similar state vector equations like Eq. (18) also exist in each new merging layer element. Repeating this coalescence operation N-1 times, the divided layer elements are recombined into the original medium and the expression of state vector can be derived as well.

^V Vb = F 1 a

(21b)

For layer block [c, d]:

(18b) (19a)

(21a)

^ V+ + F ^ TU U+b = Q 2 b 2 c

(18a)

F3 = F2 (I + G1Q2) 1F1

^ U G 2 c

^ V+ Vc = F 2 b

Subsequently, we will merge the layer element 1 and layer element 2 to build new layer element named layer element 3, as shown in Fig. 2. Accordingly, the relationship expression of the state vectors at the top and bottom of the new layer element can be obtained:

^ Q ^ 1^ ^ =F ^ (I + G F 12 2 1 2 ) F1

(24a)

^ =G ^ +F ^ Q ^ 1^ ^T ^ (I + G G 12 2 2 1 2 ) G1 F 2

(24b)

^ =Q ^ + Q 12 1

^ T (I F 1

^ G ^ 1^ ^ +Q 2 1) Q2 F1

^ 1 ^ =F ^ (I + G Q R 1 2 1 2) , ^ ^G S^1 = R 1 1,

S^1 =

^ = R 2

(24c)

^ S^2 Q 2

^ G ^ ^ T (I + Q F 2 1) 2

(24d) (24e)

1

From which we have discussed above, the state vector of calculation plane can be obtained as above when Hc > Hl . Furthermore, when Hc < Hl , namely, the load plane is located below the calculation plane. In this condition, interface b becomes the calculation plane and interface c becomes the load plane. Similarly, the state vector of the calculation plane can also be obtained:

^V Vb = F 1 a

^ U G 1 b

^ G ^ 1 ^ ^ ^T Ub = (I + Q 23 1) (Q23 F1Va + F 23Ud

For layer block [b , c]: 5

(25a)

^ P R 4 V

S^4 PU )

(25b)

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

^ are corresponding to that of Eq. where the relational matrices like F 23 (24) obviously, and it is convenient to obtain their expressions by replacing the subscript “1” and “2” with “3” and “4”. Particularly, when we face the case that Hc = Hl , Hl = 0, and H c = 0 , the thickness of special layer block is 0. According to the definition of Eq. (16), the calculation can still continue ^ , and Q ^ as I , 0, and 0 ^, G by setting the relevant relational matrices F respectively. In terms of the load term PV and PU , with the aid of Laplace and Hankel transform, we can obtain the tangential circular uniform load along the horizontal direction ( = 0 ):

PV =

0,

2pb J1 ( b), 0, 0 s

and PU = 0

(26a)

Similarly, for vertical circular uniform load and point sink induced by pumping, we can also get the corresponding load term separately:

PV =

pb J1 ( b), 0, 0, 0 s

PV = 0 and

PU =

and PU = 0

0, 0, 0,

Q 2 s2

Fig. 5. The comparison of results with Ref. [36].

(26b)

dimensionless factor is

=

2Gk z 2t. wb

Meanwhile, for a multi-layered

medium example in Ref. [36], as shown in Fig. 5, we have the diG1 uz 2Gk z mensionless factor = 2 t and uz = pb that are used to demon-

(26c)

These problems are not the focuses of our work; hence we no longer introduce that in the following sections and further work will focus on the tangential load. Based on the foregoing work, we have obtained the solution of this problem, nevertheless, it can be learned that this solution is located at the Laplace-Hankel transformed domain. Hence, after the precise integration method, the Laplace and Hankel inverse transforms [57,58] are applied to obtain the real solution.

wb

strate the accuracy of the vertical displacement in our work. It can be observed that the results of our work are both of good coincidence with those in Ref. [8] and Ref. [36]. 5.2. Verification with the FEM

In view of the limited credibility of degenerate solution, this section will employ the finite element method to further verify the correctness of our solution. Commercial software ABAQUS cannot provide the stress-strain relationship about the Merchant model to describe viscoelastic behavior of soils as well, but the user-defined material subroutine (UMAT) can achieve that. Based on the previous work [49], we will apply the UMAT subroutine to describe the stress increments in accordance with the strain increments for the Merchant model. According to that, we will build a viscoelastic model in ABAQUS to verify the correctness of our theory. As shown in Fig. 6, for simulating the infinite space and weakening the size effect, the soil model of ABAQUS is defined as a 1/4 cylindrical

5. Verification and numerical analysis 5.1. Verification with the analytical solution In order to demonstrate the convergence and accuracy of the employed method, two verification examples compared with the analytical solution are presented in this section. There are numerous researches for Biot’s consolidation of poroelastic medium with incompressible constituents. Nevertheless, to the best of authors’ knowledge, the analytical solution for the non-axisymmetric consolidation of transversely isotropic viscoelastic medium is not derived so far. Hence, we have degraded our problem to the non-axisymmetric consolidation of isotropic linear-elastic medium induced by tangential circular loads, and then verified the solution with the existing solutions [8,36]. For a single-layered medium example in Ref. [8], as shown in Fig. 4, the calculation is located at the r b = 0.5, z b = 0.1 and = 0 , which is used to verify the accuracy of the pore pressure in our work. The

Fig. 4. The comparison of results with Ref. [8].

Fig. 6. The ABAQUS model. 6

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

Fig. 8. The evolution of lateral displacement with time for different E v0 at point A.

Fig. 7. The comparison of results between ABAQUS and our method.

body with length 10b and depth 10b, where b represents the radius of load. The element type is C3D8P. As for boundary conditions, we assume that the top of soil is permeable whilst the bottom of soil is impermeable, and a uniform circular tangential load p = 1 MPa with radius b = 1 m is applied to the top of soil. The parameters of soils are 7 as follows: E v0 = 10 MPa , E v1 = 10 MPa , v1 = 1 × 10 MPa s , Kh = K z = 1 × 10 9 m s , = 1, Ks = 20 GPa , µh = µ vh = 0.25, Kf = 2 GPa . Then we calculate this model by ABAQUS and our method, respectively, and select the r = 0, z = 0 (A) and r = 0, z = b (B) as the calculation points. The results of horizontal displacement ur in Fig. 7 reveal that with the increase of time, the curves of displacement keep unchanged at first, then rise dramatically and eventually reach a steady value. A conclusion can be drawn that the results of our method have a great coincidence with those of ABAQUS, which proves the applicability and feasibility of our method. It is noted that the slight difference between the two results may stem from the size effect of the model. It can be seen distinctly that the horizontal displacement ur at r = 0, z = 0 is larger than that at r = 0, z = b , which is completely due to that the calculation A is close to the tangential circular load.

Fig. 9. The evolution of excess pore pressure with time f or different E v0 at point B.

with time factor = 86400t , namely, 1 = 1 day . It is noted from Fig. 8 that the elastic modulus E v0 has a pronounced effect on the evolutions of the lateral displacement and excess pore pressure. For different E v0 , with the decrease of the elastic modulus E v0 , the displacement increases gradually. It is worthy noted that the displacement difference between different E v0 remains unchanged, indicating that the displacement difference is mainly generated in the initial stage, i.e. elastic stage, and do not develop with time. Hence it can be deduced that the elastic modulus E v0 only has a significant effect on the instant deformation and no effect on the deformation caused by rheology. Fig. 9 shows that the curves of excess pore pressure stabilize for a period of time and then drop to zero smoothly. In order to describe the pore pressure at the initial stage of consolidation, the minimum time factor is as small as 10−7 day. The four cases reveal that the drop of elastic modulus E v0 decreases the initial excess pore pressure slightly, and postpones the start time of dissipation of pore pressure. The larger E v0 means the earlier dissipation of pore pressure, and then they lead to the relatively larger pore pressure.

5.3. Parameter analysis For the Merchant model, there are three parameters involving the time-dependent behavior of soils, i.e. the elastic modulus E v0 of spring, the elastic modulus E v1 of the Kelvin element and the viscosity coefficient of the Kelvin element v1. We will discuss the influence of the three parameters on consolidation. Meanwhile, there is no denying that the compressible constituents will also affect the consolidation of soils. In view of that the compressibility of fluid has been discussed fully in the previous studies and the compressibility coefficient of fluid does not change a lot under different cases, we will concentrate on the influence of the compressibility of soil particles. Finally, the stratification influence will be analyzed as well. 5.3.1. Elastic modulus of Hooke’s spring In order to analyze the influence of the elastic modulus E v0 of Hooke’s spring on consolidation, four cases, in which E v 0 E v1 = 0.5, 1.0, 2.0, 3.0 , are set up in this section. The other para7 E 1 = 10 MPa , meters are as follows: 1 = 1 × 10 MPa s , Gv E v1 = 0.4 , K z = 2 × 10 8 m s , = 1, Kh = 1 × 10 8 m s , Ks = 20 GPa , Kf = 2 GPa . The intensity of uniform tangential load is p = 1MPa and the radius is b = 1 m . Fig. 8 describes the variations of the lateral displacement at point A, namely, the central point of load (r b = 0, z b = 0 ). On the other hand, the dissipation of excess pore pressure is an essential signal of primary consolidation. Fig. 9 describes the dissipation of excess pore pressure at point B (r b = 0.5, z b = 0.1)

5.3.2. Elastic modulus of the Kelvin element The elastic modulus E v1 of Kelvin element is a key parameter for the secondary consolidation. Hence, it is of great importance to discuss the influence of E v1 on the consolidation of viscoelastic porous medium. In this section, four cases with E v1 E v0 = 0.5, 1.0, 2.0, 3.0 are investigated, respectively. Furthermore, the elastic modulus of Hooke’s spring is E v0 = 10 MPa , and other parameters are the same as those of previous section. Figs. 10 and 11 depict the variation of the lateral 7

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

Fig. 10. The evolution of lateral displacement with time for different E v1 at point A.

Fig. 12. The evolution of lateral displacement with time for different point A.

v1

at

Fig. 11. The evolution of excess pore pressure with time for different E v1 at point B.

Fig. 13. The evolution of excess pore pressure with time for different point B.

v1

at

displacement of calculation point A and excess pore pressure of point B against time, respectively. It is shown from Fig. 10 that, at the beginning stage, the lateral displacement does not have any changes for a long time with the fall of the elastic modulus E v1 of Kelvin element. Finally, the displacement of the case with the lowest E v1 develops quickly and reaches a highest stable value. Indeed, the reason for this is that the elastic modulus E v1 has an significant effect on the deformation subjected to rheology, as the spring E v is in parallel with a viscosity pot v1. Hence the instant displacement is identical due to the same E v0 and the difference occurs in the following consolidation. As for the excess pore pressure, the curves of four cases in Fig. 11 are identical completely. It is because that the pore pressure dissipates completely with the end of the primary consolidation, and the E v1 only alters the time-dependent behavior of rheology.

consolidation rate. Obviously, the decreasing v1 leads to the earlier beginning of secondary consolidation. Hence the curve of the lowest v1 reaches final displacement the earliest. The reason for the above phenomenon is due to that the larger viscosity coefficient, the stronger the ability of soil particles to adsorb bound water; and the strong bound water prevents the rheology of particles. Hence, the lower viscosity coefficient accelerates the secondary consolidation of soils. Nevertheless, this influence is only limited to the rate rather than the initial and final value. For excess pore pressure, the curves of four cases in Fig. 13 are coincident in relative terms, but just a little difference arises in the curve of lowest v1, which demonstrates that the viscosity coefficient v1 has a subtle effect on the dissipation of excess pore pressure when it is low enough. 5.3.4. The compressibility of soil particle The compressibility of soil particle in sand, soft rock such as halite, in fact, is non-negligible in the course of consolidation analysis. Hence, the aim of this section is to investigate the influence of the compressibility of soil particle. For the two different geomaterials such as clay and halite, six cases are divided into two group to discuss, respectively, i.e.Ks = , 4 × 10 4 MPa, 1 × 10 4 MPa and Ks = 4 × 10 4 MPa, 2 × 10 4 MPa, 1 × 10 4 MPa . Here Ks = is used to simulate the incompressible situation. The elastic modulus of Hooke’s E v 0 = E v1 = 10 MPa for elastomer is the clay while E v 0 = E v1 = 1200 MPa for the halite [59], and other parameters of both materials are identical with those of Section 5.3.1. Fig. 14 depicts the variations of lateral displacement and excess pore pressure with time

5.3.3. The viscosity coefficient of the Kelvin element The section provided four cases with 5 6 7 8 v1 = 1 × 10 , 1 × 10 , 1 × 10 , 1 × 10 MPa s to investigate the influence of the viscosity coefficient v1 in the Kelvin element on the timedependent behavior of soils. The elastic modulus of Hooke’s spring is E v0 = 10 MPa , and other parameters are identical with that of Section 5.3.1. Fig. 12 describes the time-lateral displacement curves under the four cases whilst corresponding time-excess pore pressure curves is shown in Fig. 13. As given in Fig. 12, it is obvious that the effect of the viscosity coefficient v1 of the Kelvin element on the initial and final lateral displacement is almost negligible whilst it significantly changes the 8

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

(a) Lateral displacement

(a) Lateral displacement

(b) Excess pore pressure

(b) Excess pore pressure

Fig. 14. The evolution of lateral displacement and excess pore pressure with time for different K s at point A and B (clay).

Fig. 15. The evolution of lateral displacement and excess pore pressure with time for different K s at point A and B (halite).

under various Ks for clay, and correspondingly Fig. 15 is for the results of halite. It is clearly shown in Fig. 14 that the bulk modulus of soil particle Ks does not have any influence on the lateral displacement and excess pore pressure, since the three curves overlap definitely. The reason for this is that comparing with the elastic modulus E v0 and E v1 of clay, the bulk modulus Ks is so large that the deformation of soil particle can be ignored. Nevertheless, for the halite with relatively larger E v0 and E v1, the initial lateral displacement has a slight increase due to the fall of the bulk modulus Ks as shown in Fig. 15(a); and the increase results from the compression of particles. Fig. 15(b) indicates clearly that the initial excess pore pressure drops with the decrease of Ks , which is attributed to that the compression of particles induced by the decreasing Ks allows the particles to bear more load, and the load bearded by the pore water is thus reduced.

coordinate r b = 0, z b = 1 is set as point C, and point D is the point r b = 0.1, z b = 0.1. The following three diagrams illustrate the timedependent behavior of consolidation considering the stratification and the displacement distribution along depth when the load is applied at the inner of soil, respectively. Fig. 16 shows the variations of lateral displacement with time at the point A (r = 0, z = 0 ) and point C (r = 0, z = b ), and Fig. 17 describes the variations of excess pore pressure with time at the point B (r = 0.1b, z = 0.1b ) and point D (r = 0.1b, z = 0.5b ). Obviously, the displacement at point A exceeds that at point C on account of the former near the position of load. It can be learned that the stratification has a significant influence on the time-dependent behavior of displacement and excess pore pressure. Both the displacement and excess pore pressure of multilayered soils always vary from those of single-layered soil, which is caused by the diversity of multilayered soil parameters. Turning to the reasons, it is conspicuous that there is a slight MandelCryer effect appearing at the curve of point D which proves that the Mandel-Cryer effect is supposed to be more and more obvious as the calculation depth increases. Fig. 18 depicts the distribution of lateral displacement, with = 0.01, 10, 10000 , subjected to a tangential load applied at the depth of 35 m. The three times is used to investigate the initial value, intermediate value and final value of lateral displacement, respectively. From the depth −0.5b to 0.5b around the load depth, the difference between multilayered medium and single-layered medium established by the weighted average method is remarkable. However, the difference between the two media from the depth −1.5b to −0.5b and 0.5b

5.3.5. The stratification of soil Stratification is a typical characteristic of geomaterials and widely appears in various engineering. To prove the application of the method to multilayered medium, a four-layered medium is established with the thickness ratio h1: h2 : h3: h4 = 4: 3: 6: 2 . For convenience, we assume that the relationship of several parameters remains unchanged as follows: K z Kh = 2 , µ vh µ v = 1.2 , and E v 0 = 1.25E v1. The corresponding parameters are illustrated in Table 1. Meanwhile, a single-layered medium that soil parameters are calculated by the weighted average 4 4 method, i.e., K z1 4 = i = 1 K zi hi i = 1 hi , is designed as a comparative case to discuss the influence of stratification. Here the point of 9

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

Table 1 The basic parameters of each layer. Layer number

E v1 (MPa)

1 2 3 4 Weighted average

8 9.6 6.4 12 10.27

v1 × 10

−7

1 2 1.5 1.8 1.51

µv

Kz × 108 (m/s)

0.25 0.27 0.22 0.24 0.245

2.0 4.0 1.6 2.0 2.0

(MPa·s)

Fig. 16. The evolution of lateral displacement with time for multilayered and single layered medium at point A and point C.

1.0 0.8 1.3 0.9 1.07

Ks (GPa)

Ks (GPa)

20 20 20 20 20

2 2 2 2 2

Fig. 18. The distribution of lateral displacement along the depth for multilayered and single layered medium ( = 0.01, 10, 10000 ).

Kelvin element v1, the compressibility of particles as well as the stratification. Several crucial conclusions can be drawn as follows: (1) The elastic modulus E v0 of spring will affect the initial displacement of soils dramatically due to a tangential loadings, in which increasing E v0 will result in the decrease of initial displacement and final displacement. Meanwhile, increasing E v0 will also alleviate the rate of excess pore pressure dissipation. (2) The elastic modulus E v1 of the Kelvin element and the viscosity coefficient of the Kelvin element v1 are both the key parameters for rheology. The increase of E v1 considerably reduces the final displacement, but remains the initial displacement unchanged. And increasing v1 will accelerate the rate of secondary consolidation remarkably. Nevertheless, the two parameters nearly have no effect on the dissipation of excess pore pressure. (3) For the clay, the influence of compressibility of soil particles on consolidation is almost negligible. However, for several soft rocks such as the halite, there is a slight growth of displacement resulting from particles compression. Furthermore, the excess pore pressure will soar initially due to the increase of bulk modulus of soil particle Ks . (4) The stratification of soils will impact the consolidation of soil significantly, and the weighted average method for soil parameters cannot accurately simulate the actual engineering although it can alleviate the amount of calculation to a certain extent.

Fig. 17. The evolution of excess pore pressure with time for multilayered and single layered medium at point B and point D.

to 1.5b is limited relatively, and the displacement outside this range is almost negligible. Generally speaking, the influence of stratification on consolidation is relatively inevitable in the engineering. 6. Conclusion With the aid of the integration transform and the extended precise integration method, this paper presents an exact solution for the consolidation problem of multilayered transversely isotropic medium with the compressible constituents. Two verification examples compared with the analytical solutions and a verification example implemented by ABAQUS with UMAT are conducted to demonstrate the correctness of this method. Then several numerical examples are calculated to investigate the influence of the elastic modulus E v0 of spring, the elastic modulus E v1 of the Kelvin element, the viscosity coefficient of the

Our work can be used to investigate the quasi-static response of soft soil induced by the vehicle brake and tidal action. Meanwhile, the presented consolidation solution can also be extended to further analyze the structure-soil interaction considering the effect of rheology, such as pile and caisson foundations. Nevertheless, this study also has several shortcomings facing the actual problem. For example, this proposed model can only handle the linear viscoelasticity problem, but

10

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al.

the geotechnical material is generally viscoelasto-plastic material. Also, this model can only accommodate small strain problem, and cannot deal with the large strain problem. We will try our best to solve these problems in the future.

Acknowledgements This research is supported by the National Natural Science Foundation of China (Grant No. 41672275).

Appendix A Nomenclature Greek letters , z r, r z, r , z

vh, v vh h, v 1 w

Jm ( r ) Symbol ur , uz Eh Ev E0 E1 Gv c11, c12, c13, c33, c44 K s , Kf

M n k h, k v Qz s t

and z directions [kPa] normal stress components in r, shear stress in r - z, r - and - z plane [kPa] excess pore pressure [kPa] Poisson’s ratios characterizing the lateral strain response in the plane of transverse isotropy to a stress acting parallel and normal to it [–] Biot-Willis coefficients in the horizontal and vertical directions [–] the viscosity coefficient of the Kelvin element [kPa s] the unit weigh of pore fluid [kN/m3] Hankel transform parameter with respect to r [–] the mth order Bessel function [–]

horizontal and vertical displacements [m] horizontal Young’s modulus [kPa] vertical Young’s modulus [kPa] Young’s modulus of spring [kPa] Young’s modulus of Kelvin element [kPa] vertical shear modulus [kPa] independent elastic parameters [kPa] the bulk modulus of soil particle and fluid [kPa]

Biot’s modulus [kPa] porosity of the material [–] horizontal and vertical coefficients of permeability [m/s] the total flow in the z direction [L] Laplace transform parameter with respect to time t [–] time [s]

[19] Yue ZQ, Yin JH. Backward transfer-matrix method for elastic analysis of layered solids with imperfect bonding. J Elasticity 1998;50(2):109–28. [20] Booker JR, Small JC. Finite layer analysis of consolidation I. Int J Numer Anal Meth Geomech 1982;6(2):151–71. [21] Booker JR, Small JC. Finite layer analysis of consolidation II. Int J Numer Anal Meth Geomech 1982;6(2):173–94. [22] Booker JR, Small JC. A method of computing the consolidation behaviour of layered soils using direct numerical inversion of Laplace transforms. Int J Numer Anal Meth Geomech 1987;11(4):363–80. [23] Aramaki G. Application of the boundary element method for axisymmetric Biot's consolidation. Eng Anal 1985;2(4):184–91. [24] Dargush GF, Banerjee PK. A boundary element method for axisymmetric soil consolidation. Int J Solids Struct 1991;28(7):897–915. [25] Wang JG, Fang SS. The state vector solution of axisymmetric Biot's consolidation problems for multilayered poroelastic media. Mech Res Commun 2001;28(6):671–7. [26] Ai ZY, Cheng YC, Zeng WZ. Analytical layer-element solution to axisymmetric consolidation of multilayered soils. Comput Geotech 2011;38(2):227–32. [27] Vardoulakis I, Harnpattanapanich T. Numerical Laplace-Fourier transform inversion technique for layered-soil consolidation problems: I. Fundamental solutions and validation. Int J Numer Anal Meth Geomech 1986;10(4):347–65. [28] Harnpattanapanich T, Vardoulakis I. Numerical Laplace-Fourier transform inversion technique for layered soil consolidation problems: II. Gibson soil layer. Int J Numer Anal Meth Geomech 1987;11(1):103–12. [29] Senjuntichai T, Rajapakse RKND. Exact stiffness method for quasi-statics of a multilayered poroelastic medium. Int J Solids Struct 1995;32(11):1535–53. [30] Pan E. Green's functions in layered poroelastic half-space. Int J Numer Anal Meth Geomech 1999;23(13):1631–53. [31] Wang JG, Fang SS. State space solution of non-axisymmetric Biot's consolidation problems for multilayered poroelastic media. Int J Eng Sci 2003;41(15):1799–813. [32] Mei GX, Yin JH, Zai JM, Yin ZZ, Ding XL, Zhu GF, et al. Consolidation analysis of a cross-anisotropic homogeneous elastic soil using a finite layer numerical method. Int J Numer Anal Meth Geomech 2004;28(2):111–29. [33] Chiou YJ, Chi SY. Boundary element analysis of Biot consolidation in layered elastic soils. Int J Numer Anal Meth Geomech 2010;18(6):377–96. [34] Ai ZY, Wang QS, Han J. Analytical solutions describing the consolidation of a multilayered soil under circular loading. J Eng Math 2010;66(4):381–93. [35] Ai ZY, Zeng WZ. Analytical layer-element method for non-axisymmetric consolidation of multilayered soils. Int J Numer Anal Meth Geomech 2012;36(5):533–45. [36] Ai ZY, Cang NR. Non-axisymmetric Biot consolidation analysis of multi-layered saturated poroelastic materials with anisotropic permeability. Soils Found

References [1] Barrat JL, Bocquet L. Large slip effect at a nonwetting fluid-solid interface. Phys Rev Lett 1999;82(23):4671–4. [2] Wang W, Kosakowski G, Kolditz O. A parallel finite element scheme for thermohydro-mechanical (THM) coupled problems in porous media. Comput Geosci 2009;35(8):1631–41. [3] Parra JO. The transversely isotropic poroelastic wave equation including the Biot and the squirt mechanisms: theory and application. Geophys 1997;62(1):309–18. [4] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12(2):155–64. [5] Biot MA. Theory of elasticity and consolidation for a porous anisotropic solid. J Appl Phys 1955;26(2):182–5. [6] Biot MA. Theory of deformation of a porous viscoelastic anisotropic solid. J Appl Phys 1956;27(5):459–67. [7] McNamee J, Gibson RE. Plane strain and axially symmetric problem of the consolidation of a semi-infinite clay stratum. Q J Mech Appl Math 1960;13(2):210–27. [8] Schiffman RL, Fungaroli AA. Consolidation due to tangential loads. Proc 6th Int Conf Soil Mech and Found Eng 1965;1:188–92. [9] Gibson RE, Schiffman RL, Pu SL. Plane strain and axially symmetric consolidation of a clay layer on a smooth impervious base. Q J Mech Appl Math 1970;23(4):505–20. [10] Booker JR, Randolph MF. Consolidation of a cross-anisotropic soil medium. Q J Mech Appl Math 1984;37(3):479–95. [11] Yue ZQ, Selvadurai APS, Law KT. Excess pore pressure in a poroelastic seabed saturated with a compressible fluid. Can Geotech J 1994;31(6):989–1003. [12] Chen SL, Zhang LM, Chen LZ. Consolidation of a finite transversely isotropic soil layer on a rough impervious base. J Eng Mech-ASCE 2005;131(12):1279–90. [13] Singh SJ, Rani S, Kumar R. Quasi-static deformation of a poroelastic half-space with anisotropic permeability by two-dimensional surface loads. Geophys J Int 2007;170(3):1311–27. [14] Ai ZY, Wang QS. A new analytical solution to axisymmetric Biot’s consolidation of a finite soil layer. Appl Math Mech-Engl 2008;29(12):1617–24. [15] Cai YQ, Geng XY. Consolidation analysis of a semi-infinite transversely isotropic saturated soil under general time-varying loadings. Comput Geotech 2009;36(3):484–92. [16] Pan E. The static response of multilayered foundations to general surface loading and body force. Acta Mech Sinica 1989;21(3):344–53. [17] Yue ZQ. On elastostatics of multilayered solids subjected to general surface traction. Quart J Mech Appl Math 1996;49(3):47–499. [18] Pan E. Static Green’s functions in multilayered half spaces. Appl. Math Model 1997;21(8):509–21.

11

Computers and Geotechnics 117 (2020) 103257

Z.Y. Ai, et al. 2013;53(3):408–16. [37] Skempton AW. Effective stress in soil, concrete and rocks. In: Proc Conf Pore Press and Suct Soils, Butterworths, London; 1960. p. 4–16. [38] Rice JR, Cleary MP. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents. Rev Geophys Space Phys 1976;14(2):227–41. [39] Chen GJ. Consolidation of multilayered half space with anisotropic permeability and compressible constituents. Int J Solids Struct 2004;41(16–17):4567–86. [40] Singh SJ, Kumar R, Rani S. Consolidation of a poroelastic half-space with anisotropic permeability and compressible constituents by axisymmetric surface loading. J Earth Syst Sci 2009;118(5):563–74. [41] Ai ZY, Hu YD. Multi-dimensional consolidation of layered poroelastic materials with anisotropic permeability and compressible fluid and solid constituents. Acta Geotech 2015;10(2):263–73. [42] Taylor DW, Merchant W. A theory of clay consolidation accounting for secondary compressions. J Math Phys 1940;19(3):167–85. [43] Gibson RE, Lo KY. A theory of consolidation for soils exhibiting secondary compression. Acta Polytech Scand 1961:296 Chapter 10. [44] Lo KY. Secondary compression of clays. J Soil Mech Found Div 1961;87(SM4):61–87. [45] Singh SJ, Rosenman M. Quasi-static deformation of viscoelastic half-space by a displacement dislocation. Phys Earth Planet In 1974;8(1):87–101. [46] Cai YQ, Xu CJ, Yuan HM. One-dimensional consolidation of layered and viscoelastic soils under arbitrary loading. Appl Math Mech-Engl 2001;22(3):353–60. [47] Qin AF, Sun DA, Zhang JL. Semi-analytical solution to one-dimensional consolidation for viscoelastic unsaturated soils. Comput Geotech 2014;62:110–7.

[48] Xie KH, Xie XY, Li XB. Analytical theory for one-dimensional consolidation of clayey soils exhibiting rheological characteristics under time-dependent loading. Int J Numer Anal Meth Geomechh Geomech 2008;32(14):1833–55. [49] Ai ZY, Zhao YZ, Song XY, Mu JJ. Multi-dimensional consolidation analysis of transversely isotropic viscoelastic saturated soils. Eng Geol 2019;253(10):1–13. [50] Ai ZY, Cheng YC. Extended precise integration method for consolidation of transversely isotropic poroelastic layered media. Comput Math Appl 2014;68(12):1806–18. [51] Cheng AHD. Material coefficients of anisotropic poroelasticity. Int J Rock Mech Min 1997;34(2):199–205. [52] Simon BR, Zienkiewicz OC, Paul DK. An analytical solution for the transient response of saturated porous elastic solids. Int J Numer Anal Meth Geomech 2010;8(4):381–98. [53] Sneddon IN. The Use of Integral Transforms. New York: McGraw-Hill; 1972. [54] Lee EH. Stress analysis in visco-elastic bodies. Q Appl Math 1955;13(2):183–90. [55] Kovarik V. Distributional concept of the elastic-viscoelastic correspondence principle. J Appl Mech 1995;62(4):847. [56] Zhong WX. Duality system in applied mechanics and optimal control. Boston: Kluwer Academic Publisher; 2004. [57] Abate J, Valkó PP. Multi‐precision Laplace transform inversion. Int J Numer Meth Eng 2010;60:979–993. [58] Ai ZY, Yue ZQ, Tham LG, Yang M. Extended Sneddon and Muki solutions for multilayered elastic materials. Int J Eng Sci 2002;40(13):1453–83. [59] Mctigue DF. Thermoelastic response of fluid-saturated porous rock. J Geophys ResSol Ea 1986;91(B9):9533–42.

12