Computers and Geotechnics 38 (2011) 227–232
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Analytical layer-element solution to axisymmetric consolidation of multilayered soils Zhi Yong Ai ⇑, Yi Chong Cheng, Wen Ze Zeng Department of Geotechnical Engineering, Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai, China
a r t i c l e
i n f o
Article history: Received 9 September 2010 Received in revised form 28 November 2010 Accepted 28 November 2010 Available online 24 December 2010 Keywords: Analytical layer-element solution Axisymmetric consolidation Multilayered soils Negative exponential function
a b s t r a c t After the application of a Laplace–Hankel transform, the governing equations of Biot’s consolidation are solved analytically by using the eigenvalue approach. Then the analytical layer-element of a single soil layer can be obtained in the transformed domain by synthesizing the generalized displacements and stresses, which are both expressed by six arbitrary constants. The elements of the analytical layer-element only contain negative exponential functions, which leads to a considerable improvement in computation efficiency and stability. The global stiffness matrix equation of multilayered soils is further obtained by assembling the interrelated layer-elements, and the actual solution is achieved by numerical inversion of the Laplace–Hankel transform after the solution in the transformed domain is obtained. Numerical examples are given to demonstrate the accuracy of this method and to study the influence of the layered soil properties and time history on the consolidation behavior. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction The problem of the time-dependent behavior of saturated soil is of great importance in foundation engineering. Up to now many theories of consolidation have been established to solve this problem, among which the theory of Biot’s consolidation [1] considering the coupling between the solid and the fluid in saturated soils is regarded as the most reasonable one to model the timedependent behavior of saturated soil. In the past several decades, numerous methods have been proposed to obtain the solutions of the equations of Biot’s consolidation, and one of the most successful methods is the displacement function method. McNamee and Gibson [2,3] proposed two displacement functions to solve plane strain and axisymmetric consolidation problems with a semi-infinite body, and their approach was further extended to get the solutions of non-axisymmetric consolidation problems by Schiffman and Fungaroli [4] through three displacement functions. Though the displacement function method is convenient in deriving the solution of the equations of Biot’s consolidation, it is difficult to find the displacement functions satisfying more complex consideration problems like as a consideration problem with anisotropy. In consideration of the fact that most soils are naturally layered in character, some researchers shifted focus to the study of the consolidation problems of multilayered soils. The finite layer ⇑ Corresponding author. Address: Department of Geotechnical Engineering, Tongji University, 1239 Siping Road, Shanghai 20092, China. Tel.: +86 21 65982201; fax: +86 21 65985210. E-mail address:
[email protected] (Z.Y. Ai). 0266-352X/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2010.11.011
method [5–7] and the transfer matrix method [8–12] were usually adopted to study the consolidation of multilayered soils. However, for the elements of the stiffness matrix and transfer matrix used in the conventional methods mentioned generally contain positive exponential functions, researchers have to make great efforts against ill-conditioned matrices and computation overflow. In this paper, the analytical layer-element method is presented for the purpose of eliminating the influence caused by the positive exponential functions. By employing the eigenvalue approach after the application of a Laplace–Hankel transform, the governing equations of Biot’s consolidation are solved analytically. Then the generalized displacements and stresses, which are both expressed by six arbitrary constants, are synthesized to acquire the analytical layer-element of a single soil layer in the transformed domain. In order to get the analytical layer-element matrix with its elements expressed by negative exponential functions, great efforts have been made in simplifications. Because of the elimination of positive exponential functions in the elements of the analytical layer-element, it is numerically efficient and stable. By combining the boundary conditions and the continuity conditions between adjacent layers, the global stiffness matrix equation of multilayered soils is gotten. It has to be noted that the analytical layer-element is different from the transfer matrix of a single soil layer, the former is a symmetric stiffness matrix describing the relationship between generalized displacements and stresses, while the later is about the relationship between the surface boundary vectors (including both displacements and stresses) and the bottom boundary vectors of a single soil layer. It also can be seen that, as to multilayered problems, the final transfer matrix is obtained by
228
Z.Y. Ai et al. / Computers and Geotechnics 38 (2011) 227–232
multiplying matrices, while the globe stiffness matrix is acquired by assembling the interrelated layer-elements in the analytical layer-element method, which reduces the operation of matrix. Once the solution in the transformed domain has been found, the actual solution is obtained by numerical inversion of the Laplace–Hankel transform. At last, numerical examples are given to demonstrate the accuracy of this method and study the influence of the layered soil properties and time history on the consolidation behavior.
2 r d u du nr r þ ng z ¼ n2 ð1 þ gÞu ; dz2 dz G 2 z d u 1 du 1 dr z ng r þ ; ¼ n2 u 2 1þg dz dz G dz 2 1 2 d r du þ scw z þ nscw u r : n kr ¼ k dz2 dz
ð5aÞ ð5bÞ ð5cÞ
The system of Eq. (5) may be written as:
d Wðn; z; sÞ ¼ Aðn; sÞWðn; z; sÞ; dz
2. Governing differential equations
ð6Þ
in which, In the axisymmetric consolidation problem, all variables are independent of the coordinate h, therefore the governing differential equations can be written as follows:
1 @e 1 @ r ur þ g ¼ 0; r2 @r G @r @e 1 @ r r2 uz þ g ¼ 0; @z G @z @e k ¼ r2 r; @t cw
r2 ur
2
ð1aÞ ð1bÞ ð1cÞ
2
@ where r2 ¼ @r@ 2 þ 1r @r@ þ @z 2 is the Laplacian operator for the axisymmetric consolidation; ur, uz are the displacements in the r, z direcz r tions, respectively; e ¼ @u þ urr þ @u denotes the dilatation; @r @z Eð1mÞ 1 g ¼ 12m ; M ¼ ðg þ 1ÞG ¼ ð1þmÞð12mÞ in which E, G and m are Young’s modulus, the shear modulus, and Poisson’s ratio of the soil, respectively; C denotes the coefficient of consolidation, which can be expressed as C ¼ ck M; k, cw and r are the coefficient of w permeability, unit weight of water and excess pore pressure, respectively. According to the constitutive equations of Biot’s consolidation, the following relationships can be obtained:
@ur @uz þ ; @z @r @ 1 @uz þ ur þ Gðg þ 1Þ ¼ Gðg 1Þ r: @r r @z
" W¼
U dU dz
#
;
I ; A1
0 A¼ A2
2
3 r u 6 7 U ¼ 4u z 5; r
2
3 0 ng 0 6 7 0 1=ðGð1 þ gÞÞ 5; A1 ¼ 4 ng=ð1 þ gÞ 0 scw =k 0 2 6 A2 ¼ 4 2
n2 ð1 þ gÞ
0
n=G
0
n2 =ð1 þ gÞ
0
nscw =k
0
n2
0 0 0
2
3
6 7 0 ¼ 4 0 0 0 5;
1 0
6 I ¼ 40 0
0 0 0
0
3 7 5;
3
7 1 0 5: 0 1
To solve Eq. (6), we assume that:
Wðn; z; sÞ ¼ Xðn; sÞehz :
ð7Þ
According to Eqs. (6) and (7), we have:
rrz ¼ G
ð2aÞ
hWðn; z; sÞ ¼ Aðn; sÞWðn; z; sÞ;
rz
ð2bÞ
which leads to the eigenvalue problem. The eigenvalue can be obtained by solving the characteristic equation as follows:
In accordance with Darcy’s law, the total flow Q in the z direction within the time from 0 to t is defined as:
Q ¼
Z
t
0
k @r dt: cw @z
ð3Þ
ð8Þ
detðA hIÞ ¼ 0;
ð9Þ
and the eigenvalues are:
h1;2 ¼ n;
h3;4 ¼ n;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h5;6 ¼ n2 þ s=C ¼ q:
ð10Þ
Thus the solution of Eq. (6) can be expressed as: 3. Analytical layer-element for Biot’s consolidation of a single soil layer In order to reduce partial differential equations of consolidation to ordinary differential equations, we take the Laplace–Hankel transform of the variables of the equations above. The Laplace– Hankel transform (mth order) of function f(r, z, t) with respect to the variables t and r respectively and its inversion are defined by Sneddon [13]:
f ðn; z; sÞ ¼
Z
1
0
f ðr; z; tÞ ¼
1 2pi
Z Z
1
f ðr; z; tÞJ m ðnrÞre
0 cþi1
ci1
st
Z 0
1
dr dt;
f ðn; z; sÞJ ðnrÞnest dn ds; m
ð4aÞ
Wðn; z; sÞ ¼ ðC 1 þ C 2 zÞenz þ C 3 eqz þ ðD1 þ D2 zÞenz þ D3 eqz ;
ð11Þ
where Ci and Di (i = 1–3) associated with n and q are not completely independent functions. The respective expressions of these six functions can be obtained by substituting Eq. (11) into Eq. (6), so we get:
C 1 ¼ m1 l1 þ m2 l2 ;
C 2 ¼ m1 l2 ;
C 3 ¼ m3 l3 ;
¼ m4 l4 þ m5 l5 ;
D2 ¼ m4 l5 ;
D3 ¼ m6 l6 ;
D1 ð12Þ
where li (i = 1–6) are arbitrary constants, and
m1 ¼ ½1=n; 1=n; 0; 1; 1; 0T ; m2 ¼ ½0; 1=n2 ; 2G=n; 1=n; 0; 2GT ;
ð4bÞ
where f ðn; z; sÞ is the corresponding variable of f(r, z, t) in the Laplace–Hankel transformed domain; s and n are the Laplace and Hankel transform parameters, respectively; Jm(nr) denotes the mth order Bessel function. Taking the Laplace and the Hankel transform of Eqs. (1a–c) with respect to t and r of order 1, 0, and 0, respectively, we have:
m3 ¼ ½kn=qscw ; k=scw ; 1=q; kn=scw ; kq=scw ; 1T ; m4 ¼ ½1=n; 1=n; 0; 1; 1; 0T ; m5 ¼ ½0; 1=n2 ; 2G=n; 1=n; 0; 2GT ; m6 ¼ ½kn=qscw ; k=scw ; 1=q; kn=scw ; kq=scw ; 1T : Substituting Eq. (12) into Eq. (11), W(n, z, s) can further be written as a matrix form:
229
Z.Y. Ai et al. / Computers and Geotechnics 38 (2011) 227–232
Wðn; z; sÞ ¼ Nðn; z; sÞL;
ð13Þ
U(ξ , 0, s)
0
V(ξ, 0, s)
where
a singer layer
L ¼ ½l1 ; l2 ; l3 ; l4 ; l5 ; l6 T ; Nðn; z; sÞ ¼ 2
N1
N2
N3
N4
kneqz =qscw
enz =n2 þ zenz =n
keqz =scw
2
2Ge =n
enz =n
kneqz =qscw
6 N 2 ¼ 4 enz =n enz =n2 zenz =n nz
N3 ¼
2Ge
dN 1 ; dz
N4 ¼
7 5;
Since the matrix X has the character of the analytical solution, in this paper it is called the analytical layer-element for axisymmetric consolidation problems, whose elements are provided in Appendix A after sufficient simplifications.
e =q
zenz =n
0
3
7 keqz =scw 5;
4. Analytical layer-element solution to Biot’s consolidation of multilayered soils
eqz =q
=n
dN 2 : dz
Similarly, the applications of the Laplace and the Hankel transform of Eqs. (2a), (2b) and (3) with respect to t and r of order 1, 0, and 0, respectively, result in:
r du z ; nu dz du r r z ¼ G ðg 1Þnur þ ðg þ 1Þ z ; dz G k dr : Q ¼ scw dz
r rz ¼ G
ð14aÞ ð14bÞ ð14cÞ
rz ; r z , Q may be represented by variables Obviously, variables r and their first derivatives with respect to z, which have alr ; u z ; r u ready be obtained in Eq. (13). Then, substituting Eq. (13) into Eq. (14), we get:
Vðn; z; sÞ ¼ Mðn; z; sÞL;
ð15Þ
2
Mðn; z; sÞ ¼ ½ M 1
2
2
M 2 ;
2zenz
2kneqz =scw
0
2kenz =scw
keqz =Gscw
2enz
2zenz
2enz 6 M 1 ¼ G4 2enz
6 M 2 ¼ G4 2enz 0
A soil system, which has been divided into n layers determined by the number of natural layers and calculation planes, is shown schematically in Fig. 2. Let the thickness of the ith layer be hi = Hi Hi1, in which Hi and Hi1 are the depths from the surface to the bottom and top of the ith layer, respectively. An axisymmetric normal load q(r, Hj, t) at the depth Hj is considered. Obviously, j = 0, is a special case, in which the load is acting on the surface of the soil system. Supposing that the surface of the soil system is permeable and the bottom of the soil system is fixed and impermeable, so we have:
rrz ðr; 0; tÞ ¼ rðr; 0; tÞ ¼ rz ðr; 0; tÞ ¼ 0; ur ðr; Hn ; tÞ ¼ uz ðr; Hn ; tÞ ¼ Q ðr; Hn ; tÞ ¼ 0; Applying Eq. (16) to each layer yields:
Vðn; Hi1 ; sÞ Vðn; Hi ; sÞ
¼ XðiÞ
Uðn; Hi1 ; sÞ Uðn; Hi ; sÞ
;
3
7 2ðzenz enz =nÞ 2kn2 eqz =qscw 5;
2kneqz =scw
3
7 2ðzenz þ enz =nÞ 2kn2 eqz =qscw 5: 2kenz =scw
keqz =Gscw
According to Eqs. (13) and (15), the relationship between U(n, z, s) and V(n, z, s) may be expressed as follows:
Vðn; 0; sÞ
Vðn; z; sÞ
¼X
Uðn; 0; sÞ Uðn; z; sÞ
;
ð16Þ
normal load q
where U(n, 0, s), V(n, 0, s) are the values when z = 0, respectively, and
X¼
M 1 ð0Þ M 2 ð0Þ M 1 ðzÞ
M 2 ðzÞ
N 1 ð0Þ N 2 ð0Þ N 1 ðzÞ
N 2 ðzÞ
ð18aÞ ð18bÞ
ð19Þ
where X(i) = X(n, hi, s) denotes the analytical layer-element of the ith layer. The global stiffness matrix of the multilayered soils is assembled through considering the continuity conditions of the interfaces between the adjacent layers:
where
3 r rz 6 7 Vðn; z; sÞ ¼ 4 rz 5; Q
V(ξ, z, s)
3
qz
nz
0
U(ξ , z, s)
Fig. 1. The generalized stresses and displacements for a single soil layer.
zenz =n
enz =n
6 N 1 ¼ 4 enz =n
z
;
1 ;
ð17Þ
is an exact symmetric matrix of order 6 6, which establishes the relationship between the generalized displacements and stresses in the transformed domain of a single soil layer, as shown in Fig. 1.
Fig. 2. Multilayered soils under circular uniform loadings.
230
Z.Y. Ai et al. / Computers and Geotechnics 38 (2011) 227–232
domain can be derived through solving Eq. (20), and the real solutions can be obtained by taking the inversion of the Laplace–Hankel transform to these variables.
.. .
..
.
5. Numerical results and discussion
.. .
In order to get the real solutions, it is necessary to take the inversion of the Laplace–Hankel transform to the solutions mentioned above. In this study, the method proposed by Talbot [14] for the numerical inversion of the Laplace transform is applied, whose feasibility and accuracy has been verified by Booker and Small [5–7]. As to the inversion of the Hankel transform, the technique suggested by Ai et al. [15] is adopted.
ð20Þ
.. . .. .
5.1. Verification In order to verify the accuracy of the foregoing theory, some numerical results are presented. In Fig. 3, the consolidation of a single soil layer under a uniform circular vertical loading is considered to validate the axisymmetric solution. At surface region r 6 b, a loading is applied at t = 0+ and is thereafter kept constant, while the boundary conditions are that the surface is permeable and the bottom of the soil layer is fixed and impermeable. Calculation results based on the proposed theory and the solutions given by
where [0, 0 , . . . 0, V(n, Hj, s) , 0 . . . 0, V(n, Hn, s)]T are the external forces of the soil system in the transformed domain; [U(n, H0, s), U(n, H1, s) . . . U(n, Hn, s)]T are the generalized displacements in the transformed domain at the interfaces. Combined with the boundary conditions Eq. (18), the unknown variables ½Uðn; H1 ; sÞ; Uðn; H2 ; sÞ . . . Uðn; Hn ; sÞT in the transformed
τ = ct / b 2 0.001 0
0.01
0.1
1
10
Gu z / bq
0.1
0.2 b
b
0.3
h / b = 1, v = 0
q
r / b = 0, z / h = 0 c = 2Gk / γ w
r h
0.4
z 0.5 Fig. 3. Consolidation of a single soil layer under axisymmetric loading.
τ = ct / b 2 0.0001 0.10
0.001
0.01
0.1
1 Approximate result Exact result
G1u z / bq
0.20
b
b
0.30
q
r h1, k1 , G1
b / h = 1 5, v = 0.3, r / b = 1, z / h = 0, c = 2G1k1 / γ w .
h2 , k2 , G2
0.40
h3 , k3 , G3 h4 , k4 , G4
z
0.50 Fig. 4. Displacement of the four-layered soil under axisymmetric loading.
10
231
Z.Y. Ai et al. / Computers and Geotechnics 38 (2011) 227–232
σ q 0
0
0.1
0.2
0.3
0.4
0.5
0.6
h / b = 1, c = 2G1k1 / γ w ,
0.2
τ = ct / b 2 . 0.4
b
z h
b
q
r h1, k1 , G1
0.6 τ=0.1 0.8
τ=0.05 τ=0.005
h2 , k2 , G2 h3 , k3 , G3 h4 , k4 , G4
z
1 Fig. 5. Excess pore pressure of the four-layered soil under axisymmetric loading.
Booker and Small [6] are both shown in Fig. 3, where the dimensionless settlement at the centre of the loaded area is plotted against time factor s. It may be clearly seen that the two results are in good agreement, which gives evidence of the accuracy of the foregoing theory. 5.2. Examples of multilayered soils In this section, we take several examples to study the consolidation of multilayered soils due to an axisymmetric loading. As shown in Figs. 4 and 5, a vertical loading of intensity q and radius b is uniformly applied at the surface of the four-layered soils, and the relationships of the parameters of the four soil layers are: h1:h2:h3:h4 = 1:2:1:1, G1:G2:G3:G4 = 1:3:4:2, k1:k2:k3:k4 = 2:1:1:4, where hi, Gi, ki are the thickness, shear modulus, and permeability of the ith layer, respectively, and the total thickness of four-layered soil is h = h1 + h2 + h3 + h4. Poisson’s ratio is assumed to be a constant value of m = 0.3 for each layer. The boundary conditions are the same as that of the previous section. In Fig. 4, the dimensionless vertical displacement at the center of the loaded area on the surface is plotted against the time factor s. A homogeneous singlePsoil layer with the weighted average shear modulus (i.e. G ¼ 4i¼1 Gi hi =h ¼ 2:6G1 Þ and the weighted P average permeability (k ¼ 4i¼1 ki hi =h ¼ 1:8k1 Þ are introduced to evaluate the influence of layered soil properties on the consolidation of multilayered soils. The difference between the actual results and approximate results is remarkable, which demonstrates layered soil properties have a great effect on the consolidation of multilayered soils. Another example is shown in Fig. 5 to study the variation of excess pore pressure along direction z with time growth. The dimensionless excess pore pressure at the points (r = b) along the direction z is plotted when s = 0.005, s = 0.05 and s = 0.1. We note that at the beginning there is a peak value at the depth of approximately 0.3 h, which disappears with the growth of the time. While both in the cases of s = 0.05 and s = 0.1, the excess pore pressure increases with the depth growth without a peak value.
analytically. Based on the obtained solutions the analytical layerelement for a single soil layer can be derived. Combined with the continuity conditions between adjacent layers and the boundary conditions of the layered soils, the global stiffness matrix is acquired by assembling the interrelated layer-elements. The unknown generalized displacements variables in the transformed domain can be derived in accordance with the synthesized global stiffness matrix, and the real solutions are obtained by taking the inversion of the Laplace–Hankel transform of these variables. One example is carried out to validate the accuracy of the theory and computer program in this study, and the others are performed to further study the consolidation of multilayered soils subjected to axisymmetric loadings, which prove that the differences of soil properties between layers have a great effect on the consolidation of multilayered soils. In comparison with conventional methods like as the finite layer method and the transfer matrix method, the analytical layer-element method possesses the advantage of no existence of positive exponential functions and reduction of matrix operation, which leads to a considerable improvement in computation efficiency and stability. Acknowledgements The work reported here was supported by the National Natural Science Foundation of China (Grant No. 50578121). The authors wish to express their gratitude for this financial assistance. Appendix A. The elements of analytical layer-element X (a symmetric matrix)
X11 ¼ 2GMsnð2CGnðd1 q d1 d2 d5 nÞ d1 Msd2 Þ=T; X12 ¼ 4Gn2 ð2a22 d1 M2 s2 z2 n þ 2C 2 d2 G2 nd3 þ CGMsðd1 d22 n qd4 ÞÞ=T;
X13 ¼ 2CGnð2d5 CGnd3 þ Msð4a2 znd5 þ d2 d6 ÞÞ=T; X14 ¼ 4GMsnða2 d1 Msðd7 2d2 Þ þ 2CGnðd2 d5 2qd5 d6 ÞÞ=T;
6. Conclusions By employing the eigenvalue approach after the application of a Laplace–Hankel transform, the governing equations are solved
X15 ¼ 4d2 GMsn2 ða2 d1 Msz 2d6 CGqÞ=T; X16 ¼ 4CGnð2a2 CGnd3 Msðd2 d5 þ a2 znd6 ÞÞ=T;
232
Z.Y. Ai et al. / Computers and Geotechnics 38 (2011) 227–232
X22 ¼ 2GMsnð2Cd2 Gnd6 d1 Msd8 Þ=T; 2
X23 ¼ 2CGnð2CGnd2 d3 Msðd1 d2 n þ qðd4 þ 4a2 d6 znÞÞÞ=T;
d1 d2 q2 2d3 qn þ d1 d2 n2 ; d4 ¼ d2 d3 8a2 d6 zn, d5 ¼ a2 d1 n a1 d2 q; d6 ¼ d2 d4 q d1 d5 n, d7 ¼ d2 þ d5 zn; d8 ¼ d2 d5 þ 4a22 zn, d9 ¼ 2 2 d2 4a22 z2 n2 ; T ¼ 4C 2 d2 G2 n2 d3 þ d1 M2 s2 d9 þ 4CGMsnðd1 d2 n qðd4 þ 4a2 d6 znÞÞ.
X24 ¼ X15 ; References
X25 ¼ 4GMsnð2CGnd2 d5 a2 d1 d7 MsÞ=T; 2
X26 ¼ 4CGMsnða1 qd8 þ a2 d1 d2 zn qa2 d7 d4 Þ=T; X33 ¼ k0 ð4C 2 d5 G2 n3 d3 2CGMsnð8a1 a2 d2 qzn2 þ d1 ðn2 d2 q2 d8 ÞÞ þ d4 M2 qs2 d9 Þ=sT;
X34 ¼ X16 ; X35 ¼ X26 ; 0
X36 ¼ 2k ð2a2 CGnð2CGn2 d3 þ Msðd1 q2 d7 d2 n2 ðd1 þ d4 qzÞ zn2 d6 ÞÞ þ a1 d9 M2 qs2 Þ=sT;
X44 ¼ X11 ; X45 ¼ X12 ; X46 ¼ X13 ; X55 ¼ X22 ; X56 ¼ X23 ; X66 ¼ X33 ; 0
where k ¼ k=cw ; q2 ¼ n2 þ s=C; a1 ¼ eqz ; a2 ¼ enz ; d1 ¼ a21 þ 1; d2 ¼ a22 þ 1, d3 ¼ ða1 a2 Þ2 þ ða1 a2 1Þ2 ; d4 ¼ a21 þ 1; d5 ¼ a22 þ 1, d6 ¼ ða1 a2 1Þða1 a2 Þ; d1 ¼ d3 d5 þ 4a2 d6 ; d2 ¼ d2 d5 4a22 zn, d3 ¼
[1] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155–64. [2] McNamee J, Gibson RE. Displacement functions and linear transforms applied to diffusion through porous elastic media. Quart J Mech Appl Math 1960;13:98–111. [3] McNamee J, Gibson RE. Plane strain and axially symmetric problem of the consolidation of a semi-infinite clay stratum. Quart J Mech Appl Math 1960;13:210–27. [4] Schiffman RL, Fungaroli AA. Consolidation due to tangential loads. In: Proceedings of the 6th international conference on soil mechanics and foundation engineering, Montreal, Canada, vol. 1; 1965. p. 188–92. [5] Booker JR, Small JC. Finite layer analysis of consolidation I. Int J Numer Anal Methods Geomech 1982;6:151–71. [6] Booker JR, Small JC. Finite layer analysis of consolidation II. Int J Numer Anal Methods Geomech 1982;6:173–94. [7] Booker JR, Small JC. A method of computing the consolidation behavior of layered soils using direct numerical inversion of Laplace Transforms. Int J Numer Anal Methods Geomech 1987;11:363–80. [8] Pan E. Green’s functions in layered poroelastic half-space. Int J Numer Anal Methods Geomech 1999;23:1631–53. [9] Wang JG, Fang SS. The state vector solution of axisymmetric Biot’s consolidation problems for multilayered poroelastic media. Mech Res Commun 2001;28:671–7. [10] Wang JG, Fang SS. State space solution of non-axisymmetric Biot consolidation problems for multilayered poroelastic media. Int J Eng Sci 2003;41:1799–813. [11] Ai ZY, Han J. A solution to plane strain consolidation of multi-layered soils. In: Luna R, Hong ZS, Ma GW, Huang MS, editors. Soil and rock behavior and modeling. ASCE Geotechnical Special Publication. Proceedings of the GeoShanghai international conference 2006, Shanghai, China, vol. 150, June 6–8, 2006. p. 276–83. [12] Ai ZY, Cheng ZY, Han J. State space solution to three-dimensional consolidation of multi-layered soils. Int J Eng Sci 2008;46:486–98. [13] Sneddon IN. The use of integral transform. New York: McGraw-Hill; 1972. [14] Talbot A. The accurate numerical inversion of Laplace transforms. IMA J Appl Math 1979;23:97–120. [15] Ai ZY, Yue ZQ, Tham LG, Yang M. Extended Sneddon and Muki solutions for multilayered elastic materials. Int J Eng Sci 2002;40:1453–83.