Analytical solution to axisymmetric consolidation in unsaturated soils with linearly depth-dependent initial conditions

Computers and Geotechnics 74 (2016) 102–121

Contents lists available at ScienceDirect

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

Research Paper

Analytical solution to axisymmetric consolidation in unsaturated soils with linearly depth-dependent initial conditions Liem Ho, Behzad Fatahi ⇑, Hadi Khabbaz School of Civil and Environmental Engineering, University of Technology Sydney (UTS), Sydney, Australia

a r t i c l e

i n f o

Article history: Received 11 August 2015 Received in revised form 25 December 2015 Accepted 31 December 2015

Keywords: Axisymmetric consolidation Unsaturated soils Analytical solution Vertical drain wells Excess pore pressures Degree of consolidation

a b s t r a c t This paper introduces an analytical solution for the axisymmetric consolidation of unsaturated soils subjected to constant external loading. The analytical procedure employs variables separation and Laplace transformation techniques while capturing the uniform and linear initial excess pore pressure distributions with depth. Excess pore-air and pore-water pressures as functions of time, radial and vertical flows are determined using Laplace transforms, Fourier Bessel and sine series, respectively. In this study, the consolidation behavior, in terms of changes in excess pore-air and pore-water pressures and the average degree of consolidation, are investigated against the air to water permeability ratio. The effects of radial distance from the drain well on the dissipation rate are likewise highlighted in worked examples. Excess pore pressure isochrones and the matric suction varying with time are also presented. Crown Copyright Ó 2016 Published by Elsevier Ltd. All rights reserved.

1. Introduction In major ground improvement projects, vertical drain assisted preloading has been a cost-effective method to accelerate drainage in soil deposits and shorten the compression process. Vertical drain consolidation can be modeled assuming axisymmetry around the drain. The basic theory of the axisymmetric consolidation was developed based on the traditional theory proposed by Terzaghi [1]. Barron [2] first introduced solutions for axisymmetric consolidation of saturated soils while considering radial flow only. Yoshikuni and Nakanodo [3] introduced a rigorous free strain equation adopting permeable top and bottom boundaries of a cylindrical soil mass while considering the well resistance. Based on the equal strain condition, Hansbo [4] developed a simplified closed-form solution for the degree of consolidation induced by drain wells capturing the well resistance and peripheral smear effects. Further studies have included the effects of the smear zone and time-dependent loading [5–11], varying radial permeability coefficients [12–14], new analysis for soil–drain system conceptually based on the double porosity model (DPM) [15], vertical drain consolidation with discharge capacity varying linearly with depth

⇑ Corresponding author at: School of Civil and Environmental Engineering, Faculty of Engineering and Information Technology, University of Technology Sydney (UTS), City Campus, PO Box 123, Broadway, NSW 2007, Australia. Tel.: +61 (2) 9514 7883, mobile: +61 413573481; fax: +61 (2) 9514 2633. E-mail address: [email protected] (B. Fatahi). http://dx.doi.org/10.1016/j.compgeo.2015.12.019 0266-352X/Crown Copyright Ó 2016 Published by Elsevier Ltd. All rights reserved.

and decreasing exponentially with time [16], or axisymmetric model including both radial and vertical flows for electroosmotic consolidation [17]. The above mentioned studies are only valid for saturated soils. In real practice, drain wells may be installed through multi-layered geomaterials, where groundwater table is found to be reasonably below the ground surface, and often unsaturated soil layers may be found near the soil surface. In other words, the capillary zones do not extend to the soil surface in many situations. This unsaturated soil portion could be a result of past earthworks (e.g. excavation, compaction) or climate changes (e.g. arid or semi-arid climates) [18]. It is believed that the presence of unsaturated soil deposits within the entire strata, where the drain wells are installed, may significantly influence the consolidation behavior and likewise characterize the settlement pattern of a soil. Thus, there has been an essential need for a generalized model that can be applicable for unsaturated soils. Consolidation studies for unsaturated soils have progressed rigorously since the inception of one-dimensional (1D) consolidation theory proposed by Fredlund and Hasan [19]. This theory introduces the nonlinear governing equations describing independent flows of air and water in unsaturated soil deposits. Later, Dakshanamurthy and Fredlund [20], and Dakshanamurthy et al. [21] expanded the existing equations to two-dimensional (2D) and three-dimensional (3D) models, respectively. This set of theories has inspired a large number of recent studies, of which analytical and numerical models have been conducted in the Cartesian

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

103

List of symbols A B Ca Cw cav r cav z cw vr cw vz g H i j kar kaz kwr kwz M ma1 ma2 ms1 ms2 mw 1 mw 2 n q0 R S Sr s T T a ðtÞ T w ðtÞ t U ua ua;t ua;r

constant of eigenfunction with respect to the radial domain constant of eigenfunction with respect to the radial domain interactive constant associated with the air phase interactive constant associated with the water phase coefficient of consolidation with respect to the air phase in the radial domain coefficient of consolidation with respect to the air phase in the vertical domain coefficient of consolidation with respect to the water phase in the radial domain coefficient of consolidation with respect to the water phase in the vertical domain gravitational constant thickness of the soil stratum integer for the Fourier Bessel series integer for the Fourier sine series coefficient of permeability for the air phase in the radial domain coefficient of permeability for the air phase in the vertical domain coefficient of permeability for the water phase in the radial domain coefficient of permeability for the water phase in the vertical domain molecular mass of the air phase coefficient of volume change of the air phase with respect to the change of net stress coefficient of volume change of the air phase with respect to the change of suction coefficient of volume change of the soil with respect to the change of net stress coefficient of volume change of the soil with respect to the change of suction coefficient of volume change of the water phase with respect to the change of net stress coefficient of volume change of the water phase with respect to the change of suction porosity external loading universal air constant center to center spacing in the vertical drain system degree of saturation matric suction time factor generalized Fourier coefficients varying with time with respect to the air phase generalized Fourier coefficients varying with time with respect to the water phase elapsed time average degree of consolidation excess pore-air pressure first order of partial differential equation (PDE) of air with respect to time first order of partial differential equation (PDE) of air with respect to radius

coordinate system. Further unsaturated soil studies include 1D consolidation settlement of a single layer unsaturated soil [22–27], consolidation of multi-layered soil [28], consolidation for viscoelastic soil [29], and 2D plane strain consolidation problems [30–32]. There have been, however, few attempts to model axisymmetric unsaturated consolidation particularly with analyti-

ua;rr

h h° kija

second order of partial differential equation (PDE) of air with respect to radius second order of partial differential equation (PDE) of air with respect to depth atmospheric pressure initial excess pore-air pressure at the soil surface initial excess pore-air pressure at depth z average excess pore-air pressure pore-water pressure first order of partial differential equation (PDE) of water with respect to time first order of partial differential equation (PDE) of water with respect to radius second order of partial differential equation (PDE) of water with respect to radius second order of partial differential equation (PDE) of water with respect to depth initial excess pore-water pressure at the soil surface initial excess pore-water pressure at depth z average excess pore-water pressure initial volume of the soil element volume of air within soil element volume of water within soil element eigenfunction with respect to the excess pore-air pressure for domain r eigenfunction with respect to the excess pore-water pressure for domain r investigated radius radius of zone of influence radius of the drain well eigenfunction with respect to the excess pore-air pressure for domain z eigenfunction with respect to the excess pore-water pressure for domain z investigated depth water unit weight suction change total volumetric strain average total volumetric strain parameter controlling the distribution of the initial excess pore-air pressure parameter controlling the distribution of the initial excess pore-water pressure absolute temperature in Kelvin polar angle temperature in degree celsius separation constant with respect to the air phase

kijw

separation constant with respect to the water phase

lj

eigenvalue for vertical boundary condition (i.e. PTIB and PTPB) eigenvalue for radial boundary condition total stress in r-domain total stress in z-domain total stress in h-domain

ua;zz uatm u0a uaz
a u uw uw;t uw;r uw;rr uw;zz u0w uwz
w u V0 Va Vw Ra ðrÞ Rw ðrÞ r re rw Z a ðzÞ Z w ðzÞ z

cw

Ds

ev


ev fa fw

H

ni

rr rz rh

cal approaches. Among pioneer studies, Conte [33] introduced the finite element technique to obtain a solution for the coupled consolidation under the plane strain and axisymmetric conditions. Qin et al. [34] dealt with the drain well consolidation problem in unsaturated soils using the modified Bessel functions and the Laplace transformation. On the other hand, Zhou and Tu [35],

104

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

and Zhou [36] presented the differential quadrature method (DQM) to estimate the axisymmetric consolidation behavior in the unsaturated soil stratum. Nevertheless, the aforementioned approaches are generally confined to complex numerical procedures. Moreover, some prediction methods may be impractical due to the lack of consideration of the vertical flow of air and water phases which occur in the field. There has hitherto been limited analytical research about the axisymmetric consolidation in unsaturated soils. Among existing studies, Qin et al. [34] introduced the method that comes closest to reproducing analytical solutions. However, the proposed method given by Qin et al. [34] adopted Crump’s Fourier series approximation [37] to numerically solve the complex Laplace transformed equations and subsequently obtain the solution. This paper proposes an analytical solution to predict the consolidation of unsaturated soils considering both radial and vertical flows. In addition, both uniformly and linearly distributed initial excess pore-air and pore-water pressures along the depth will be captured in this study. The uncoupled polar governing equations of flow with respect to air and water phases are first derived from the 3D Cartesian equations given by Fredlund et al. [38]. Based on an assumption of the free strain case, separation of variables and Laplace transformation techniques are used in this study to obtain the final solution. Graphical presentations are provided to demonstrate the effects of the air to water permeability ratio, different radii, and depth-dependent initial conditions on the changes in excess pore-air and pore-water pressures. Additionally, the matric suction change, the degree of consolidation and the excess pore pressure isochrones against time are graphically presented.

of the influence zone. When a constant loading is applied to the soil, the air and water phases will dissipate through the boundary of the drain well and through the permeable top boundary only or permeable top and base boundaries of the soil simultaneously. In the literature for consolidation theory of unsaturated soils, the flows of air and water phases within the soil are assumed to be independent and continuous [19–21]. Particularly, the air flow follows Fick’s law, whereas Darcy’s law is employed to describe the flow of water. Both air and water flows are usually simulated under the Cartesian governing equations for unsaturated soils [19–21]. Prior to the evaluation of the axisymmetric consolidation, the Cartesian equations should be transformed into the polar forms. The detailed polar transformation of governing equations is provided in Appendix A. Through mathematical derivations obtained from the net flux of the air and water phases per unit volume and constitutive equations for soil deformation, a set of uncoupled governing flow equations can be presented as follows: (a) Polar governing equations considering the radial flow only:


 1 ua;t þ C a uw;t þ cav r ua;rr þ ua;r ¼ 0 r

ð1aÞ


 1 uw;t þ C w ua;t þ cw v r uw;rr þ r uw;r ¼ 0

ð1bÞ

(b) Polar governing equations considering both radial and vertical flows:

2. Polar governing equations of flow The installation of vertical drain wells can be categorized in square and triangular patterns. According to Barron [2], most drain well systems are installed in triangular patterns due to economic reasons. Fig. 1a presents the plan of the drain well pattern and the influence zone of each well. Fig. 1b depicts the details of the typical well within an unsaturated soil stratum. Dimensions of the soil system include a depth H and a radius of the influence zone re . The vertical drain well with the radius rw is located at the center


 1 ua;t þ C a uw;t þ cav r ua;rr þ ua;r þ cav z ua;zz ¼ 0 r

ð2aÞ


 1 w uw;t þ C w ua;t þ cw v r uw;rr þ r uw;r þ cv z uw;zz ¼ 0

ð2bÞ

where ua and uw are the excess pore-air and pore-water pressures ðkPaÞ, respectively; ua;t and uw;t are the first order derivatives against time for excess pore-air and pore-water pressure functions, respectively; ua;r and uw;r are the first order derivatives against radius for excess pore-air and pore-water pressure functions, respectively; ua;rr and uw;rr are the second order derivatives against radius for excess

Vertical sand drain

Plan view

S re rw

S

r z

re H

Unsaturated soil

Zone of influence

(a)

(b)

Side view Fig. 1. Vertical drain system: (a) triangular drain well pattern and (b) details of a typical well.

105

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

pore-air and pore-water pressure functions, respectively; ua;zz and uw;zz are the second order derivatives against depth for excess pore-air and pore-water pressure functions, respectively. Additionally, the consolidation coefficients are:

C a ¼  a m

1 ma2

cav r ¼



1

rÞ  1  ma nuð1S 0 þu atm Þ 2ð a



ð3aÞ

kar RH 1 h   i gM ma u0 þ uatm maa1  1  nð1  Sr Þ a 2 m

ð3bÞ

kaz RH 1 h   i gM ma u0 þ uatm ma1a  1  nð1  Sr Þ a 2 m

ð3cÞ

2

cav z ¼

uniform initial condition may be only applicable when the external applied loading is indefinite and uniformly applied on the surface of the soil. In a case when footings built on the ground are small or the soil stratum is thick, the initial excess pore pressures generated by loading may decrease with depth. A simplified simulation for such initial condition can be presented adopting a linear distribution of excess pore pressures [40,42–44]. In this study, simulations for the linear distribution of initial excess pore pressures at r 2 ðr w ; re Þ and z 2 ð0; HÞ, are adopted as follows:

 z ua ðr; z; 0Þ ¼ u0a 1  fa H

ð6aÞ

 z uw ðr; z; 0Þ ¼ u0w 1  fw H

ð6bÞ

2

mw C w ¼ 1w  1 m2

ð3dÞ

cw vr ¼


 1 kwr mw cw 2

ð3eÞ

cw vz ¼


 1 kwz mw cw 2

ð3fÞ

where kar and kwr are the coefficients of air and water permeability in the radial domain ðm=sÞ, respectively; kaz and kwz are the coefficients of air and water permeability in the vertical domain ðm=sÞ, respectively; ma1 and mw 1 are the coefficients of air and water volume change with respect to the change of net stress (1/kPa), respectively; ma2 and mw 2 are the coefficients of air and water volume change with respect to the change of suction (1/kPa), respectively; g is the gravitational constant (i.e. 9.81 m/s2); u0a is the initial pore-air pressure at the soil surface (kPa); uatm is the atmospheric pressure (i.e. 100 kPa); R is the universal air constant (i.e. 8.314 J/mol K); H ¼ ðh þ 273Þ, is the absolute temperature (K); h° is the temperature in degree Celsius; M is the molecular mass of air (i.e. 0.029 kg/mol); n is the porosity during consolidation process; Sr is the degree of saturation during consolidation process; and cw is the water unit weight (9.81 kN/m3). 3. Analytical solution

where u0a and u0w are the initial pore-air and pore-water pressures at the soil surface, respectively; and fa and fw are the dimensionless parameters controlling the gradients of initial pore-air and porewater pressure distributions. Note that the proposed parameters fa and fw would range from 0 to 1. Fig. 2 shows the distribution of initial pore pressures in response of changes of fa and fw . When fa ¼ fw ¼ 0, the initial excess pore pressures are distributed uniformly along z-direction (Fig. 2a); when fa and fw parameters are between 0 and 1, the initial condition presents trapezoidal distribution (Fig. 2b); and when fa ¼ fw ¼ 1, the triangular distribution of initial excess pore pressures is produced (Fig. 2c). 3.2. Excess pore pressure dissipation A natural soil usually consists of nonlinear properties which may further require arduous numerical approaches to estimate the soil compression. In order to achieve an analytical solution, the following assumptions have been made in this study: (i) (ii) (iii) (iv)

The soil stratum is homogeneous. The air and water flows are continuous and independent. The soil grain and pore-water are incompressible. Environmental factors (e.g. air diffusion, temperature change) are neglected. (v) Well resistance and smear effects are discarded. (vi) Consolidation parameters for air phase (C a ; cav r and cav z ) and w water phase (C w ; cw v r and cv z ) are constant for the applied stress increment.

3.1. Boundary and initial conditions It is assumed that the boundary of the influence zone of the soil (i.e. r ¼ r e ) is impermeable whereas the surface of the drain cylinder (i.e. r ¼ r w ) is permeable, thus,

ua ðr w ; z; t Þ ¼ uw ðrw ; z; t Þ ¼ ua;r ðre ; z; t Þ ¼ uw;r ðre ; z; t Þ ¼ 0

ð4Þ

In addition, the pervious top–impervious base (PTIB) and pervious top–pervious base (PTPB) boundary conditions can be presented as below: (a) The PTIB condition:

ua ðr; 0; tÞ ¼ uw ðr; 0; t Þ ¼ ua;z ðr; H; t Þ ¼ uw;z ðr; H; tÞ ¼ 0

ð5aÞ

(b) The PTPB condition:

ua ðr; 0; tÞ ¼ uw ðr; 0; t Þ ¼ ua ðr; H; tÞ ¼ uw ðr; H; t Þ ¼ 0

ð5bÞ

Most literature and laboratory experiments have assumed that, immediately after loading, the initial excess pore pressures are distributed uniformly throughout depth [39–41]. However, this

It should be noted that the unit weight of water ðcw Þ is assumed to remain constant, implying that pore-water is incompressible, as mentioned in Assumption (iii). On the other hand, Assumption (vi) may not be completely accurate for some cases. In particular, the soil properties such as permeability (ka and kw ), degree of saturation ðSr Þ, and porosity ðnÞ in a soil stratum vary during the consolidation induced by a constant load, thus the consolidation parameters would eventually change. However, it may be acceptable to assume that these parameters are constant during the transient process for a particular stress increment. For the sake of generality, Eqs. (2a) and (2b) are considered for the analytical procedure. Firstly, general solutions for ua and uw can be written as products of functions with respect to dimensions r and z and time t. In addition, based on the homogeneous boundary conditions for the depth given in Eqs. (5a) and (5b), Fourier sine series can represent functions of the depth z. Hence, 1 X ua ðr; z; t Þ ¼ Ra ðrÞ  Z a ðzÞ  T a ðtÞ ¼ U a ðr; t Þ sin j¼0

pffiffiffiffiffiffi !

lj

H

z

ð7aÞ

106

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

Surcharge at t = 0 0

Unsaturated soil

Drain well

u 0a

u 0w

u 0a

u 0w

u 0a

u 0w

u 0a

u 0w

z

H

0

Elevation of soil (a) Uniformly distributed initial

0

0

(c) Triangularly distributed initial

(b) Linearly distributed initial

excess pore-air and pore-water pressures

0

excess pore-air and pore-water pressures

excess pore-air and pore-water pressures

Fig. 2. Initial conditions: (a) uniform, (b) trapezoidal and (c) triangular distributions of initial excess pore pressures.

1 X uw ðr; z; t Þ ¼ Rw ðrÞ  Z w ðzÞ  T w ðtÞ ¼ U w ðr; t Þ sin

pffiffiffiffiffiffi !

j¼0

lj

H

ð7bÞ

z

where Ra ðrÞ ¼ Rw ðrÞ ¼ RðrÞ, are the eigenfunctions for the domain  pffiffiffiffiffiffi  r; Z a ðzÞ ¼ Z w ðzÞ ¼ sin z l j =H , are the eigenfunctions for the domain z; T a ðtÞ and T w ðtÞ are the generalized Fourier coefficients with respect to air and water phases, respectively; U a ðr; t Þ ¼ RðrÞT a ðtÞ and U w ðr; t Þ ¼ RðrÞT w ðtÞ. Note that Ra ðrÞ resembles to Rw ðrÞ, and Z a ðzÞ resembles to Z w ðzÞ due to similar boundary conditions, as provided in Eqs. (4) and (5). Under the PTIB boundary condition, the term l j can be defined as:

l j ¼ ½ð2j þ 1Þp=22 ; ðj ¼ 0; 1; 2; . . .Þ

ð8aÞ

ð8bÞ

Then substituting Eqs. (7a) and (7b) into Eqs. (2a) and (2b), respectively, results in:

  1 U a;t ðr; t Þ þ C a U w;t ðr; t Þ þ cav r U a;rr ðr; t Þ þ U a;r ðr; tÞ r

lj H

2

U a ðr; t Þ ¼ 0

ð9aÞ

  1 U w;t ðr; tÞ þ C w U a;t ðr; t Þ þ cw v r U w;rr ðr; t Þ þ r U w;r ðr; t Þ  cw vz

lj H2

U w ðr; tÞ ¼ 0

ð9bÞ

Applying the separation of variables to Eqs. (9a) and (9b) gives:

 kar

   R;rr ðrÞ 1 R;r ðrÞ lj 1 T a;t ðtÞ T w;t ðtÞ þ  kaz 2 ¼  a þ Ca ¼ kija RðrÞ r RðrÞ cv T a ðtÞ T a ðtÞ H

ð10aÞ

 kwr

8 i < R;rr ðrÞ þ 1r R;r ðrÞ þ ðrn Þ2 RðrÞ ¼ 0 w   ; ði ¼ 0;1;2; .. .; j ¼ 0; 1; 2;. ..Þ : T ðtÞ þ C T ðtÞ  cw kij T ðtÞ ¼ 0 w;t w a;t v w w ð12bÞ owing to the fact that

8 ij < ka ¼ kar

ni ðr w Þ2

þ kaz Hl2 j

;

ði ¼ 0; 1; 2; . . . ; j ¼ 0; 1; 2; . . .Þ

ð13Þ

w

l j ¼ ðjpÞ2 ; ðj ¼ 0; 1; 2; . . .Þ

 cav z

ð12aÞ

lj ni : kij ¼ k wr ðr Þ2 þ kwz H2 w

l j for the PTPB boundary condition is given by:

or

8 i < R;rr ðrÞ þ 1r R;r ðrÞ þ ðrn Þ2 RðrÞ ¼ 0 w   ; ði ¼ 0; 1;2; .. .; j ¼ 0;1; 2; .. .Þ : T ðtÞ þ C T ðtÞ  ca kij T ðtÞ ¼ 0 a;t a w;t v a a

   R;rr ðrÞ 1 R;r ðrÞ lj 1 T w;t ðtÞ T a;t ðtÞ  k wz 2 ¼  w ¼ kijw þ þ Cw RðrÞ r RðrÞ cv T w ðtÞ T w ðtÞ H

Note that ni and l j are eigenvalues while kija and kijw are constant values ði ¼ 0; 1; 2; . . . ; j ¼ 0; 1; 2; . . .Þ. Based on Eq. (12), Fourier Bessel series are taken as the solution for the radial function Ri ðrÞ ði ¼ 0; 1; 2; . . .Þ, resulting in:

0qffiffiffiffi 1 0qffiffiffiffi 1 i n ni Ri ðrÞ ¼ AJ 0 @ rA þ BY 0 @ rA rw rw

ð14Þ


qffiffiffiffi  where A and B are the Bessel constants; and J0 r ni =r w and
qffiffiffiffi  Y 0 r ni =r w are the Bessel functions of the first kind and the second kind of order zero, respectively. Incorporating the radial boundary condition in Eq. (4) into Eq. (14) yields in:

8
qffiffiffiffi
qffiffiffiffi > i > > n þ BY ni ¼ 0 AJ 0 < 0
qffiffiffiffi
qffiffiffiffi > > > AJ 1 . ni þ BY 1 . ni ¼ 0 :

where . ¼ r e =rw ; and J 1

ð15Þ


qffiffiffiffi
qffiffiffiffi . ni and Y 1 . ni are the Bessel func-

where

tions of the first kind and the second kind of order one at r ¼ re . To assure A – 0 as well as B – 0, the following condition must be satisfied:

cav ¼

J0

ð10bÞ
 RH 1 1 1 h  ma  i ; and cw : ¼ w v m gM ma u0 þ uatm c 1 w 2  1  nð1  Sr Þ a 2 ma


qffiffiffiffi
qffiffiffiffi
qffiffiffiffi
qffiffiffiffi ni Y 1 . ni  J 1 . ni Y 0 ni ¼ 0

ð16Þ

2

ð11Þ Then, two sets of ordinary differential equations (ODEs) can be obtained as below:

Note that the eigenvalue ni ði ¼ 0; 1; 2; . . .Þ is the ith root of Eq. (16). Besides, the relationship between constants A and B in Eq. (15) can be obtained as follows:

107

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121


qffiffiffiffi J 1 . ni B ¼ A
qffiffiffiffi Y 1 . ni

ð17Þ

Combining Eqs. (14) and (17), the radial function Ri ðrÞ ði ¼ 0; 1; 2; . . .Þ becomes:

0qffiffiffiffi 1 ni A i@ R ðrÞ ¼
qffiffiffiffi D rA rw Y 1 . ni i

ð18Þ

0qffiffiffiffi 1 0qffiffiffiffi 1 0qffiffiffiffi 1
qffiffiffiffi
qffiffiffiffi ni ni ni i i D rA ¼ J0 @ r AY 1 . n  J 1 . n Y 0 @ r A: rw rw rw

To obtain the generalized Fourier coefficients T ija ðtÞ and T ijw ðtÞ ði ¼ 0; 1; 2; . . . ; j ¼ 0; 1; 2; . . .Þ, the Laplace transformation should be first applied to the function of time in Eqs. (12a) and (12b), thus,

h i  i ij ij ij ij ij s  cw v kw T w ðsÞ þ sC w T a ðsÞ  T w ð0Þ þ C w T a ð0Þ ¼ 0

Y1 T ija ðtÞ ¼

ð20bÞ

T ijw ðtÞ ¼

h i T ijw ð0ÞðC w C a  1Þs þ cav T ijw ð0Þ þ C w T ija ð0Þ kija   T ijw ðsÞ ¼ ij a w a ij ij ðC w C a  1Þs2 þ kijw cw v þ ka cv s  cv cv kw ka

ð21bÞ

Constant terms T ija ð0Þ and T ijw ð0Þ ði ¼ 0; 1; 2; . . . ; j ¼ 0; 1; 2; . . .Þ in Eq. (21) can be obtained using the orthogonalities of Bessel and sine functions, as shown below:

0qffiffiffiffi 1 pffiffiffiffiffi ! ni lj ua ðr; z; 0ÞD r A sin z rdrdz H rw 0 rw 2 0qffiffiffiffi 1 pffiffiffiffiffi !32 Z H Z re ni lj 5 i 4D @ r A sin rdrdz z rw H rw 0

A

Y1


qffiffiffiffi . ni

¼

A

Na

A

Y1 ¼ where


qffiffiffiffi . ni A

ð23Þ


qffiffiffiffi     . ni X0 eaij1 t  eaij2 t þ W0 eaij1 t þ eaij2 t 2gij

A

ð24aÞ

ð24bÞ

 12 2 ij ij ij a w cav kija  cw v kw þ 4cv cv C w C a ka kw ;

! ! ij ij a ij w ij ij 1 cav kija þ cw v kw þ g ; aij ¼ 1 cv ka þ cv kw  g ; 2 2 2 1  CwCa 1  CwCa   X ¼ cav kija  cwv kijw Na  2cwv C a kijw Nw ; W ¼ gij Na ;

aij1 ¼





X0 ¼ cwv kijw  cav kija Nw  2cav C w kija Na ;

and W0 ¼ gij Nw :

ð25Þ

Substituting Eqs. (18) and (24) into Eq. (7) gives the solutions as shown below:  ij   ij 3 0qffiffiffiffi 1 ij ij pffiffiffiffiffi !2 1 X ea1 t  ea2 t þ W ea1 t þ ea2 t j X ni l i@ 5 z 4 ua ðr;z;tÞ ¼ D rA sin H rw 2gij j¼0 ð26aÞ    ij 3 0qffiffiffiffi 1 ij ij pffiffiffiffiffi !2 0 aij1 t i 1 X e  ea2 t þ W0 ea1 t þ ea2 t j X n l 5 z 4 uw ðr; z;t Þ ¼ Di @ r A sin rw 2gij H j¼0 ð26bÞ


qffiffiffiffi Y 1 . ni T ijw ð0Þ ¼

i@

2gij

A Y1

ð21aÞ

T ija ð0Þ ¼


qffiffiffiffi     . ni X eaij1 t  eaij2 t þ W eaij1 t þ eaij2 t

ð20aÞ

h i T ija ð0ÞðC w C a  1Þs þ cw T ija ð0Þ þ C a T ijw ð0Þ kijw v   T ija ðsÞ ¼ ij a w a ij ij ðC w C a  1Þs2 þ kijw cw v þ ka cv s  cv cv kw ka

re

dz;

obtain the generalized Fourier coefficients T ija ðtÞ and T ijw ðtÞ ði ¼ 0; 1; 2; . . . ; j ¼ 0; 1; 2; . . .Þ:

gij ¼

Z

z


qffiffiffiffi
qffiffiffiffi
qffiffiffiffi
qffiffiffiffi2 2  i ¼ ðr e Þ J 0 . ni Y 1 . ni  J 1 . ni Y 0 . ni

neously for T ija ðsÞ and T ijw ðsÞ ði ¼ 0; 1; 2; . . . ; j ¼ 0; 1; 2; . . .Þ gives:

H

lj

H

pffiffiffiffiffi !# Z " Xi 4 H lj qffiffiffiffi z dz; uw ðr;z; 0Þsin H 0 H  i ni 
qffiffiffiffi
qffiffiffiffi
qffiffiffiffi
qffiffiffiffi 2 ni  J 1 ni Y 1 . ni ; and X i ¼ ðr w Þ J 1 . ni Y 1

where

Z

ua ðr; z;0Þsin 0

Nw ¼

where T ija ðsÞ and T ijw ðsÞ ði ¼ 0; 1; 2; . . . ; j ¼ 0; 1; 2; . . .Þ are the Laplace transformed equations with complex subjugate s. Solving simulta-


qffiffiffiffi Y 1 . ni

H

Combining Eqs. (21)–(23) and then taking the Laplace inverse to

ð19Þ

h

Z

2

i@

h i  i s  cav kija T ija ðsÞ þ sC a T ijw ðsÞ  T ija ð0Þ þ C a T ijw ð0Þ ¼ 0

Xi 4 qffiffiffiffi H  i ni

pffiffiffiffiffi !#

!2 
qffiffiffiffi
qffiffiffiffi
qffiffiffiffi
qffiffiffiffi2 Xi i i i i  ðr w Þ J 0  : n Y 1 . n  J1 . n Y 0 n rw

where

h

Na ¼

"

ð22aÞ

0qffiffiffiffi 1 pffiffiffiffiffi ! ni lj uw ðr;z; 0ÞD r A sin z rdrdz rw H rw 0 2 0qffiffiffiffi 1 pffiffiffiffiffi !32 Z H Z re ni lj 5 i 4D @ rA sin rdrdz z rw H rw 0

Z

H

Nw

Z

re

i@

ð22bÞ

Note that the Bessel constant A presented in Eqs. (18) and (24) is canceled out when the solutions are formed. The proposed equations provided in Eqs. (26a) and (26b) present the time-dependent changes in the excess pore-air and pore-water pressures, respectively, along dimensions r and z. Similar approach can be conducted to obtain solutions for Eqs. (1a) and (1b), in which the vertical permeability coefficient ðkz Þ is discarded and both top and base boundaries of the soil are considered to be undrained. Complete set of analytical solutions for axisymmetric consolidation with or without the vertical flow are presented in Appendix B. In addition, for the consolidation with the vertical flow, both cases of uniform and linear distributions of initial excess pore pressures along the soil depth are also captured in the solutions presented in Appendix B.

108

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

 Dimensions

3.3. Average degree of consolidation
a The average excess pore-air and pore-water pressures (i.e. u
w , respectively) can be obtained by averaging the pore presand u sures in both r- and z-directions as follows:


a ðtÞ ¼ u


w ðtÞ ¼ u

1

1

p½ðre Þ2  ðrw Þ2  H 1

1 p½ðre Þ  ðrw Þ  H 2

2

Z

H

Z

2pua ðr; z; tÞrdrdz

ð27aÞ

2puw ðr; z; tÞrdrdz

ð27bÞ

H

Z

re

rw

0

The constitutive model for the soil skeleton can be adopted from Fredlund et al. [38]. It is assumed that the volume change coefficients with respect to air and water phases are constant during the consolidation process at a particular stress increment. When @ rr =@t ¼ @ rh =@t ¼ @ rz =@t ¼ 0, the constitutive relationship under the 3D loading condition is presented as below:

@

  DV V0

@t

¼ ðms2  ms1 Þ


a
w @u @u  ms2 @t @t

ð28Þ

where ms1 and ms2 are the coefficients of volume change of the soil element with respect to the change in the net stress and change in suction (1/kPa), respectively. It is worth noting that s a w ms1 ¼ ma1 þ mw 1 and m2 ¼ m2 þ m2 . By integrating Eq. (28) against the time variable t, the average volumetric strain
ev is presented as follows:


ev ðtÞ ¼ ðms2  ms1 Þ½u
a ðtÞ  u
a ð0Þ  ms2 ½u
w ðtÞ  u
w ð0Þ

ð29Þ

The average degree of consolidation, U, based on the volumetric strain
ev , can be derived as follows:

  
ev ðtÞ    U¼
ev ð1Þ

ð30Þ

where
ev ð1Þ is the final average volumetric strain when time t approaches infinity. Eq. (30) describes the time-dependent response of consolidation settlement of the unsaturated soil due to a constant load at a particular time. Complete equations presenting the average degree of consolidation U are provided in Appendix B. 4. Examples In this study, the axisymmetric consolidation with and without the vertical flows of air and water phases are investigated in two examples. The uniformly and linearly distributed initial excess pore pressures, simulated by the parameters fa and fw , are highlighted in the axisymmetric consolidation considering both radial and vertical flows. Following Conte [30,33] and Fredlund et al. [38], this study adopts:  Volume change coefficients 1

ms2 ¼ 0:2ms1 ;

1

w mw 2 ¼ 1:8m1 ;

ms1 ¼ 5:64  104 kPa ; 4 kPa ; mw 1 ¼ 1:13  10

 Material properties

n ¼ 0:50;

 Other properties kw ¼ 1010 m=s; q0 ¼ 100 kPa; u0a ¼ 20 kPa; u0w ¼ 40 kPa:

ð31Þ

rw

0

Z

re

r w ¼ 0:2 m; r e ¼ 1:8 m; H ¼ 5 m;

Sr ¼ 80%;

 Physical properties

H ¼ ðh þ 273:16Þ K; h ¼ 20  C; R ¼ 8:314 J=mol K; M ¼ 0:029 kg=mol; uatm ¼ 100 kPa;

The above soil properties are assumed to be constant during the consolidation process. An instantaneous compression induced by the external applied load q0 generates initial excess pore-air and pore-water pressures (i.e., u0a and u0w , respectively). Considering that the soil is loaded three-dimensionally under isotropic conditions, changes in excess pore pressures can be determined using a method given by Fredlund and Hasan [19], and Fredlund et al. [38]. Detailed evaluations for the initial excess pore-air and porewater pressures are provided in Appendix C. On the other hand, when considering the triangular vertical drain system, both adopted values of radii of influence zone r e and sand drain r w are realistic and practical for evaluations of axisymmetric consolidation. In particular, the center to center spacing ðSÞ for the adopted re ¼ 1:8 m is determined to be 3.4 m (i.e. re ¼ 0:525S), which is in the desirable range of spacing ð1:2—3:6 mÞ recommended by Holtz [45]. Additionally, the adopted rw ¼ 0:2m is also in the design range ð0:075—0:3 mÞ given by Holtz et al. [46] and Smoltczyk [47]. For the sake of simplicity, the diagonal permeability krz is discarded and only isotropic permeability condition (i.e. kar ¼ kaz ¼ ka and kwr ¼ kwz ¼ kw ¼ 1010 m=s) is considered in this study. Thus, the coefficients of consolidation can be obtained based on given properties presented in Eq. (31):

C a ¼ 6:53  102 ; C w ¼ 44:32  102 ;

cav r ¼ cav z ¼ 52  103 ka m2 =s 8 w 2 cw v r ¼ cv z ¼ 5:03  10 m =s

ð32Þ

Example 1 demonstrates changes in excess pore pressure dissipation and settlement due to the only radial flows of air and water. These changes are examined against the permeability ratio ðka =kw Þ and radial distances from the well ðrÞ. Example 2, on the other hand, presents the consolidation process including both radial and vertical flows, and captures: (a) uniform initial condition (i.e. fa ¼ fw ¼ 0); and (b) linear initial condition (i.e. fa ; fw > 0). The effects of ka =kw on the axisymmetric consolidation behavior under the PTIB and PTPB boundary conditions are investigated in Example 2(a) whereas significant effects of fa and fw are highlighted in Example 2(b). Additionally, the excess pore pressure isochrones and the matric suction changes ðDsÞ are presented in both examples. 4.1. Example 1 – Axisymmetric consolidation with radial flow only Since only radial flow is considered in this example, the uniform initial excess pore pressures over the radial domain, u0a and u0w , may be adopted. The changes in the normalized pore-air ðua =u0a Þ and pore-water ðuw =u0w Þ pressures, and the degree of consolidation ðUÞ are analyzed in this section. Adopting Eq. (B.3) in Appendix B, Fig. 3 illustrates dissipation curves of the excess pore pressures at different radii, while adopting ka =kw ¼ 10. In general, it can be observed that the excess poreair (Fig. 3a) and pore-water (Fig. 3b) pressures dissipate more quickly for the points closer to the drain well. When the radius considerably increases (i.e. r P 1 m), the effects of radial distance on the dissipation rate become insignificant. While the excess pore-air pressure is fully dissipated, there is a delay period in the excess pore-water pressure curve, known as plateau (Fig. 3b). It

109

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

can be observed that the increase in the radius of the point of interest results in longer plateau and thus a clear double S-shaped pattern for the excess pore-water pressure dissipation. Fig. 4 demonstrates changes in excess pore pressures with different ka =kw values. It is worth noting that the air permeability ka is varying while the water permeability kw is kept constant at

analysis. The results in this study generally demonstrate expected dissipation characteristics comparable to the given boundary condition (i.e. Eq. (4)), in which excess pore pressures are zero at r ¼ r w and approach highest values at r ¼ re . It is also noted that the dissipation rate of the excess pore-air pressure is faster than that of the excess pore-water pressure.

1010 m=s. In addition, the radius r ¼ 0:5ðr w þ r e Þ is taken for the investigation. Fig. 4a illustrates parallel curves of excess pore-air pressure due to different ka =kw values. It can be noted that, when ka =kw increases, the excess pore-air pressure is prone to dissipate

4.2. Example 2 – Axisymmetric consolidation with both radial and vertical flows

faster. When ka =kw is very high, for instance ka =kw P 103 , the excess pore-air pressure may dissipate instantaneously. This observation was also reported by Conte [30]. Fig. 4b shows the typical excess pore-water pressure dissipation varying with ka =kw when time elapses. It can be observed that the dissipation at the early stages of consolidation proceeds more quickly as the result of increasing ka =kw . After the excess pore-air pressure completely diminishes, a plateau may occur in the excess pore-water pressure patterns when ka =kw > 1. It should be noted that the plateau gets longer as ka =kw increases. Afterwards, the curves of excess porewater pressure dissipation converge to a single curve and complete dissipation occurs at approximately the same time (i.e. after t ¼ 108 s). Fig. 5 shows the normalized matric suction change ðDs=q0 Þ varying with ka =kw values at r ¼ 0:5ðrw þ re Þ. Initially, the matric suction decreases and attains the lowest value due to the excess pore-air pressure dissipation. It can be observed that the decrease in the suction proceeds faster as the result of increasing ka =kw . After the plateau period, the suction gradually increases as the excess pore-water pressure dissipates while the excess pore-air pressure diminished completely. As observed, there is no further change in the matric suction after t ¼ 108 s. Based on Eq. (B.5) in Appendix B, Fig. 6 illustrates the variations of the average degree of consolidation ðUÞ against time factor T (i.e.   T ¼ kw t= ms1 cw H2 ) with different ka =kw values. It is obvious that U

1 r = 0.6 m

ka/kw = 10

r = 0.8 m r = 1.0 m

0.75

r = 1.2 m r = 1.4 m

0.5

0.25 (a)

0 1

102

106

104 Time (s)

108

idation, U, will be investigated in this section. Throughout the analysis, the point of investigation is located at r ¼ 0:5ðr w þ re Þ and z ¼ 0:5H. Considering Eqs. (B.9) and (B.10) in Appendix B, Fig. 8 depicts the dissipation rate of excess pore pressures varying with ka =kw under the PTIB and PTPB boundaries, along with the radial drainage boundary. As expected, variations of ka =kw result in single inverse S curves for the excess pore-air pressure (Fig. 8a) and typical double inverse S curves for the excess pore-water pressure (Fig. 8b). It can be also observed that the excess pore pressure dissipations assisted by the PTIB boundaries of the soil have comparable rates with those in PTPB boundaries. This result indicates that the radial flow, induced by drain wells, controls the dissipation while the effects of the vertical flow on the dissipation process are less significant. It should be noted that the dissipation rate under the PTIB boundary condition comparatively resembles to that presented in Fig. 4. Besides, the Mandel–Cryer effect during the compression [48,49], which has been widely known as the phenomenon being characterized by the Poisson’s ratio and leading to increases in the excess pore-water pressure prior to the dissipation, can be discarded in this analysis. Rigorous numerical studies conducted by Conte [30] and Wong et al. [50] elucidate that the mentioned effect plays a minor role in the consolidation of unsaturated soils. Moreover, Fredlund et al. [38] also confirms that the consideration of Mandel–Cryer effect in the computation may be unnecessary since the uncoupled theory of consolidation provides a decent approximation to the fully coupled approach. On the other hand, Fig. 9 shows the normalized matric suction change ðDs=q0 Þ varying with time while adopting the PTIB and PTPB boundary conditions. As observed, although the matric suction in the PTPB boundary system tends to initially decrease and then increase slightly faster than that in the PTIB boundary system, the matric suctions in both boundary conditions eventually become constant at about the same time (i.e. after t ¼ 108 s). It is worth noting that, while the total stresses remain unchanged due

Normalized Pore-Water Pressure (uw/u0w)

Normalized Pore-Air Pressure (ua/u0a )

patterns consist of double inverse S curves when ka =kw P 0:1, similar to the excess pore-water pressure dissipation process. The early stages of consolidation are governed by the simultaneous dissipation of both excess pore-air and pore-water pressures; and once the excess pore-air pressure is fully dissipated in the soil system, U patterns converge to only one curve and gradually approach 1 at the later stages. It is worth noting that the later stages of consolidation of unsaturated soils resembles to the classical consolidation introduced by Terzaghi [1] since the process during this time is excess pore-air pressure-free. On the other hand, the soil stratum tends to settle more quickly at the beginning of the consolidation as ka =kw increases. Fig. 7 represents groups of excess pore pressure isochrones along the radius varying with T while ka =kw ¼ 0:1 is adopted in this

4.2.1. Uniform initial condition ðfa ¼ fw ¼ 0Þ The condition fa ¼ fw ¼ 0 indicates that, under the undrained compression, the initial excess pore pressures u0a and u0w are ubiquitously distributed throughout the soil depth. Variations of normalized pore pressures, ua =u0a and uw =u0w , and degree of consol-

1010

1 ka/kw = 10

r = 0.6 m r = 0.8 m r = 1.0 m

0.75

r = 1.2 m r = 1.4 m

0.5

0.25 (b)

0 1

102

106

104

108

1010

Time (s)

Fig. 3. Dissipation of (a) excess pore-air and (b) excess pore-water pressures varying with different radii (for radial flow only).

110

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

Normalized Pore-Water Pressure (uw/u0w)

r = 0.5(rw + re)

0

Normalized Pore-Air Pressure (ua/ua)

1

0.75

102

ka/kw = 103

0.5

10

10-1

1

0.25 (a)

0 104

102

1

r = 0.5(rw + re) ka/kw = 103

102

10

1 10-1

0.75

0.5

0.25 (b)

0 106

104

102

1

1010

108

106

1

108

1010

Time (s)

Time (s)

0

to the constant external loading, the consolidation process may be
Þ, and subseattributed initially to the increase in the net stress ðr quently to the increase in the matric suction ðsÞ during the later stages of consolidation. Fig. 10 presents changes in the average degree of consolidation   ðUÞ against time factor T (i.e. T ¼ kw t= ms1 cw H2 ) under the PTIB

-0.1

-0.2

ka/kw = 103

102

and PTPB boundary conditions, as referred to Eq. (B.12) in Appendix B. Similar to Fig. 6, the results of U show the expected double inverse S curves, when ka =kw P 0:1. It can be observed that, although there are clear differences in U for the PTIB and PTPB boundary systems at the early stages, the consolidation process

10-1

1

10

-0.3

r = 0.5(rw + re) -0.4 10 6

10 4

10 2

1

10 10

10 8

Time (s)

Average Degree of Consolidation (Ū)

Fig. 5. Matric suction change varying with different ka =kw values (for radial flow only).

0

0.25 102

ka/kw = 103

10

1

10-1

0.5

0.75

1 2 10-9

2 10-5

2 10-7

2 10-3

2 10-1

20

Time Factor T

1

time T ¼ 103 at various radii while adopting ka =kw ¼ 0:1. It can be observed that the excess pore-air and pore-water pressures tend to increase when the point of interest is farther away from the drain well. In contrast, the excess pore pressures reduce as the radius decreases. This is due to the fact that a considerable portion of pressures preferably dissipate through the permeable boundary. When the radius considerably increases (i.e. r P 1 m), the effects of radial distance on the pore pressure isochrones will be less pronounced.

0

0

Normalized Pore-Air Pressure (ua/ua)

Fig. 6. Average degree of consolidation varying with different ka =kw values (for radial flow only).

ends almost at the same time (i.e. after T ¼ 0:2 or after t ¼ 108 s). Figs. 11 and 12 represent groups of excess pore pressure isochrones varying with T (along the r- and z-directions) under the PTIB and PTPB boundary systems, respectively. The ratio ka =kw ¼ 0:1 is taken for this investigation. As observed, the isochrones in Figs. 11a and 12a satisfy the radial boundary condition given in Eq. (4), and also those in Figs. 11b and 12b satisfy the vertical boundary condition, as provided in Eqs. (5a) and (5b), respectively. It can be found that the excess pore pressures in the PTPB boundary system dissipate just slightly faster than those in the PTIB boundary system. This suggests that the effects of the vertical flow may be less significant as the soil thickness increases. Fig. 13 depicts excess pore pressure isochrones at a particular

Normalized Pore-Water Pressure (uw/uw)

Normalized Matric Suction Change (Δs/q0)

Fig. 4. Dissipation of (a) excess pore-air and (b) excess pore-water pressures varying with different ka =kw values (for radial flow only).

(a)

0.75

ka/kw = 10-1

0.5

0.25

0 0

0.25

0.5

0.75

Radial Ratio (r - rw)/(re - rw)

1

1 (b)

0.75 ka/kw = 10-1

0.5

0.25

0 0

0.25

0.5

0.75

Radial Ratio (r - rw)/(re - rw)

Fig. 7. Distribution of (a) excess pore-air and (b) pore-water pressures along the radial domain (for radial flow only).

1

111

1 r = 0.5(rw + re) z = 0.5H

0.75 PTIB PTPB

0.5

ka/kw = 103

102

10

1

10-1

0.25

(a)

0

108

106

104

102

1

1010

Normalized Pore-Water Pressure (uw/u0w)

Normalized Pore-Air Pressure (ua/u0a)

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

1 r = 0.5(rw + re) ka/kw = 103

102

10

1 10-1

z = 0.5H

0.75

0.5

PTIB PTPB

0.25

(b) 0 1

102

104

106

108

1010

Time (s)

Time (s)

Normalized Matric Suction Change (Δs/q0)

Fig. 8. Dissipation of (a) excess pore-air and (b) pore-water pressures varying with different ka =kw values under the PTIB and PTPB boundary conditions (for radial and vertical flows, uniform initial condition).

demonstrated in Fig. 14a, the increasing fa values lead to a reduction in the average initial excess pore-air pressure at the beginning of consolidation. The excess pore-air pressures (with different fa

0

PTIB PTPB

-0.1

values) begin to dissipate significantly at about t ¼ 103 s and are

-0.2

ka/kw = 103

102

10

1

10-1

-0.3

r = 0.5(rw + re) z = 0.5H -0.4

1

106

104

102

108

1010

Time (s) Fig. 9. Matric suction change varying with different ka =kw values under the PTIB and PTPB boundary conditions (for radial and vertical flows, uniform initial condition).

Average Degree of Consolidation (Ū)

0 PTIB PTPB

0.25

0.5

ka/kw = 103

102

10

1

10-1

0.75

1 2 10-9

2 10-7

2 10-5

2 10-3

2 10-1

20

Time Factor T Fig. 10. Average degree of consolidation varying with different ka =kw values under the PTIB and PTPB boundary conditions (for radial and vertical flows, uniform initial condition).

4.2.2. Linear initial condition ðfa ; fw > 0Þ This section highlights the effects of the proposed parameters fa and fw on the axisymmetric consolidation, as presented in Eqs. (6a) and (6b), respectively. Since fa ; fw > 0, the initial excess pore pressures are linearly distributed along the soil depth. In this study, the normalized pore-air ðua =u0a Þ and pore-water ðuw =u0w Þ pressures are presented and discussed. Additionally, the point at r ¼ 0:5ðr w þ re Þ and z ¼ 0:5H is selected for further analysis. Based on Eq. (B.17) in Appendix B, Fig. 14 shows the effects of fa on the excess pore pressure dissipation rates when ka =kw ¼ 10 and PTIB boundary condition are adopted. In this analysis, the parameter fa varies from 0 to 1 while fw remains constant equal to 0.25. As

then fully dissipated at almost the same time (i.e. after t ¼ 105 s). In Fig. 14b, it is observed that the increasing fa values result in a slower rate of the excess pore-water pressure dissipation during the plateau period. However, at the later stages of consolidation, excess pore-water pressures (with different fa values) fully dissipate with approximately the same rate. Fig. 15, on the other hand, investigates the influence of fw on the dissipation rates while considering ka =kw ¼ 10 and PTIB boundary condition. The parameter fw is varying from 0 to 1 whereas fa is kept constant equal to 0.25 in this case. As observed, there is no variation in the excess pore-air pressure dissipation despite different fw values (Fig. 15a). When fw increases, the average initial excess pore-water pressure decreases at the beginning of consolidation (Fig. 15b). Nevertheless, at the later stages of consolidation, all excess pore-water pressures (with different fw values) diminish at almost the same time (i.e. after t ¼ 108 s). Fig. 16 demonstrates changes in the matric suction due to variations of fa and fw when adopting ka =kw ¼ 10 and PTIB boundary condition. As shown in Fig. 16a, increasing fa leads to smaller suction at the beginning of consolidation since the average initial excess pore-air pressure reduces while the initial pore-water pressure remains unchanged. This observation is contrary to the effects of fw , in which the increasing fw leads to higher suction at the beginning of consolidation due to a considerable reduction in the average initial excess pore-water pressure, as presented in Fig. 16b. At the later stages, in both cases, matric suctions will increase gradually and approach a constant value at approximately the same time (i.e. after t ¼ 108 s). Fig. 17, referred to Eq. (B.20) in Appendix B, illustrates the changes in the average degree of consolidation ðUÞ against time   factor T (i.e. T ¼ kw t= ms1 cw H2 ) due to the effects of fa and fw , while ka =kw ¼ 10 and PTIB condition are considered. As observed, U patterns initially increase with almost similar rates and then begin to vary during the later stages of consolidation, at which the increasing fa decelerates the consolidation (Fig. 17a) whereas the increasing fw accelerates the consolidation process (Fig. 17b). This is due to the fact that U curves are dependent on the dissipation rate of excess pore-water pressure, particularly a faster rate of excess pore-water pressure dissipation leads to faster consolidation rate and vice versa. It can be observed that consolidation completes almost at the same time regardless of values of fa and fw (i.e. after t ¼ 108 s). Figs. 18 and 19 demonstrate the excess pore pressure isochrones varying with T (along r- and z-directions) considering

0

1

0.75

0.5

z = 0.5H

0.25

0 0

(a)

0.25

0.5

0.75

1

Radial Ratio (r - rw)/(re - rw)

1

0.75 ka/kw = 10-1 z = 0.5H

0.5

0.25

0 0

0.25

0

0

0.25

0.25

0.5

0.75

0.75

1

0.5

0.75 ka/kw = 10-1

ka/kw = 10-1 r = 0.5(rw + re)

1

(b)

0.5

Radial Ratio (r - rw)/(re - rw)

Depth Ratio z/H

Depth Ratio z/H

ka/kw = 10-1

Normalized Pore-Water Pressure (uw/u0w)

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

Normalized Pore-Air Pressure (ua /u a)

112

0

0.25

0.5

r = 0.5(rw + re)

1 0.75

1

0

Normalized Pore-Air Pressure (ua/u0a)

0.25

0.5

0.75

1

Normalized Pore-Water Pressure (uw/uw0 )

1

0.75

ka/kw = 10-1 z = 0.5H

0.5

0.25 2 10-2

0 0

0.25

(a)

0.5

0.75

1

Normalized Pore-Water Pressure (uw/uw0 )

Normalized Pore-Air Pressure (ua/u0a)

Fig. 11. Distribution of excess pore-air and pore-water pressures along (a) radial and (b) vertical domains under the PTIB boundary condition (for radial and vertical flows, uniform initial condition).

ka/kw = 10-1

0.5

z = 0.5H

0.25 2 10-1

0 0

0.25

0.5

0.75

1

0.75

1

Radial Ratio (r - rw)/(re - rw)

0.75

ka/kw = 10-1 r = 0.5(rw + re)

0.25

0.5

10-1

0.5

0.75 2 10-1

Depth Ratio z/H

1

2 10-2

Depth Ratio z/H

0.75

Radial Ratio (r - rw)/(re - rw)

1

ka/kw = 10-1 r = 0.5(rw + re)

0.25

0

(b)

1

0 0

0.25

0.5

0.75

1

Normalized Pore-Air Pressure (ua/ua0)

0

0.25

0.5

Normalized Pore-Water Pressure (uw/uw0 )

Fig. 12. Distribution of excess pore-air and pore-water pressures along (a) radial and (b) vertical domains under the PTPB boundary condition (for radial and vertical flows, uniform initial condition).

the PTIB and PTPB boundary systems, respectively, while adopting fa ¼ fw ¼ 0:5 and the ratio ka =kw ¼ 0:1. The pressure isochrones in Figs. 18a and 19a satisfy the radial boundary condition presented in Eq. (4). Under an undrained compression, linearly distributed initial excess pore pressures generated by the applied constant load ðq0 Þ are clearly captured in Figs. 18b and 19b. It is observed that the excess pore pressures under the PTIB boundary conditions

undergo the pressure redistribution process, in which some pressures dissipate through the permeable top surface while the rest is transferred toward the base of the soil (i.e. at z ¼ H). This phenomenon results in notable increases in both excess pore-air and pore-water pressures with time towards the bottom of the soil deposit (Fig. 18b). Considering the PTPB boundary condition, the redistribution of excess pore pressures also occurs to achieve the

113

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121 0

0.25

Depth Ratio z/H

Depth Ratio z/H

0

1.8 m 1.4 m

0.5 r = 0.4 m

0.6 m 0.8 m 1 m

0.25

1.8 m 1.4 m

0.5

1m r = 0.4 m

0.8 m

0.75

0.75

T=

T = 10-3 ka/kw = 10-1

10-3

ka/kw = 10-1 1

1 0

0.2

0.4

0.6

0.8

1

0

1.2

0.2

0.4

0.6

0.8

1

1.2

Normalized Pore-Water Pressure (uw/u0w)

Normalized Pore-Air Pressure (ua/u0a)

(a) 0

Depth Ratio z/H

0

0.25

Depth Ratio z/H

0.6 m

1.8 m 1.4 m

r = 0.4 m

0.5

0.8 m 1 m

0.6 m

0.75

1.8 m

0.25

1.4 m 0.5

1m

0.6 m

r = 0.4 m

0.8 m 0.75

T = 10-3 ka/kw = 10-1

T = 10-3 ka/kw = 10-1

1

1

0

0.2

0.4

0.6

0.8

1

1.2

0

0.2

Normalized Pore-Air Pressure (ua/u0a)

(b)

0.4

0.6

0.8

1

1.2

Normalized Pore-Water Pressure (uw/uw0 )

0

)

0

)

Fig. 13. Distribution of excess pore pressures along the vertical domain at T ¼ 103 at different radii under (a) PTIB and (b) PTPB boundary conditions (for radial and vertical flows, uniform initial condition).

1

1 (b)

ζw = 0.25 ζa = 0

0.75

ka/kw = 10

0.25 0.5

0.5

0.75 PTIB condition

0.25

1

r = 0.5(rw + re) z = 0.5H

0

102

1

104

106

108

1010

Normalized Pore-Water Pressure

Normalized Pore-Air Pressure

(a)

ζw = 0.25 ka/kw = 10

0.75 ζa = 1 0.75 0.5

0.5 0.25 PTIB condition

0.25

0

r = 0.5(rw + re) z = 0.5H 0 102

1

104

Time (s)

106

108

1010

Time (s)

0

0

)

)

Fig. 14. Dissipation of (a) excess pore-air and (b) pore-water pressures varying with fa (while fw is constant) under the PTIB boundary condition (for radial and vertical flows, linear initial condition).

1

1 (b)

ζa = 0.25 ka/kw = 10

0.75 ζw = 1

0.5

ζw = 0

PTIB condition

0.25

r = 0.5(rw + re) z = 0.5H

0 1

102

104

106

Time (s)

108

1010

Normalized Pore-Water Pressure

Normalized Pore-Air Pressure

(a)

ζa = 0.25 ζw = 0

0.75

ka/kw = 10

0.25 0.5

0.5

0.75 1 PTIB condition

0.25

r = 0.5(rw + re) z = 0.5H

0 1

102

104

106

108

1010

Time (s)

Fig. 15. Dissipation of (a) excess pore-air and (b) pore-water pressures varying with fw (while fa is constant) under the PTIB boundary condition (for radial and vertical flows, linear initial condition).

0 (a)

ζw = 0.25 ζa = 0

-0.1

0.25

0.5

0.75

1

-0.2

PTIB condition

-0.3

r = 0.5(rw + re) ka/kw = 10

z = 0.5H

-0.4 10 2

1

10 10

10 8

10 6

10 4

Normalized Matric Suction Change (Δs/q0)

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

Normalized Matric Suction Change (Δs/q0)

114

0 ζa = 0.25

(b) -0.1

ζw = 1

-0.2

0.5

0.75

0.25 PTIB condition

-0.3

0

r = 0.5(rw + re)

ka/kw = 10

z = 0.5H -0.4

10 4

10 2

1

Time (s)

10 10

10 8

10 6

Time (s)

0 (a)

ζw = 0.25 ζa = 1

0.25

0.75 0.5

0.5

0.25 0

0.75 PTIB condition 1 2 10-9 2 10-7

ka/kw = 10 2 10-5

2 10-3

2 10-1

20

Average Degree of Consolidation (Ū)

Average Degree of Consolidation (Ū)

Fig. 16. Matric suction change due to variations of (a) fa and (b) fw under the PTIB boundary condition (for radial and vertical flows, linear initial condition).

0 (b)

ζa = 0.25

0.25 ζw = 0 0.25

0.5

0.5 0.75

0.75

1

PTIB condition 1 2 10-7 2 10-9

ka/kw = 10 2 10-5

2 10-3

2 10-1

20

Time (s)

Time (s)

Fig. 17. Average degree of consolidation due to variations of (a) fa and (b) fw under the PTIB boundary condition (for radial and vertical flows, linear initial condition).

) 0

0

)

1

ka/kw = 10-1 z = 0.5H

0.75

0.5

0.25 ζa = ζw = 0.5

2 10-2 0 0

(a)

0.25

0.5

0.75

0.75

0.5

0.25 ζa = ζw = 0.5

2 10-1 0

1

0

0.25

Radial Ratio (r - rw)/(re - rw)

0.5

0

0

0.25

0.25

0.5

ka/kw = 10-1

0.75

1

0.5

ka/kw = 10-1

0.75

r = 0.5(rw + re)

r = 0.5(rw + re)

ζa = ζw = 0.5

ζa = ζw = 0.5

1

(b)

0.75

Radial Ratio (r - rw)/(re - rw)

Depth Ratio z/H

Depth Ratio z/H

1

z = 0.5H

Normalized Pore-Water Pressure

Normalized Pore-Air Pressure

ka/kw = 10-1

1 0

0.25

0.5

Normalized Pore-Air Pressure

0.75 0

1

)

0

0.25

0.5

Normalized Pore-Water Pressure

0.75 0

1 )

Fig. 18. Distribution of excess pore pressures along (a) radial and (b) vertical domains under PTIB boundary condition while adopting fa ¼ fw ¼ 0:5 (for radial and vertical flows, linear initial condition).

pressure equilibrium. Evidently, the maximum values of excess pore pressures at a particular time no longer occur consistently at the mid-depth but move toward the base of the soil due to the pressure redistribution (Fig. 19b).

Fig. 20 shows excess pore pressure isochrones at T ¼ 103 at various radii while adopting fa ¼ fw ¼ 0:5 and ka =kw ¼ 0:1. Similar to the result presented in Fig. 13, the excess pore-air and porewater pressures increase when the point is farther away from

115

0

0

)

)

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

1 ka/kw = 10-1

ka/kw = 10-1

z = 0.5H

z = 0.5H

0.75

0.5

0.25 ζa = ζw = 0.5

2 10-2

0 0

0.25

(a)

0.5

0.75

1

Normalized Pore-Water Pressure

Normalized Pore-Air Pressure

1

0.5

0.25 ζa = ζw = 0.5

2 10-1

0 0

0.25

Radial Ratio (r - rw)/(re - rw)

0.5

0.75

0

ka/kw = 10-1

0.75

0.5

10-1

10-2

2 10-2

0.5

0.25 2 10-1

Depth Ratio z/H

0.25

ka/kw = 10-1

0.75

r = 0.5(rw + re)

r = 0.5(rw + re)

ζa = ζw = 0.5

ζa = ζw = 0.5

1

1 0

(b)

1

Radial Ratio (r - rw)/(re - rw)

0

Depth Ratio z/H

0.75

0.25

0.5

0.75

1 0

Normalized Pore-Air Pressure

0

0.25

)

0.5

0.75

1 0

Normalized Pore-Water Pressure

)

Fig. 19. Distribution of excess pore pressures along (a) radial and (b) vertical domains under the PTPB boundary condition while adopting fa ¼ fw ¼ 0:5 (for radial and vertical flows, linear initial condition).

0

Depth Ratio z/H

Depth Ratio z/H

0

0.25 r = 0.4 m

0.6 m

0.8 m 1 m

0.5 1.8 m

T = 10-3

0.75

1.8 m 1.4 m 1m T = 10-3 ζa = ζw = 0.5

1

1 0

0.25

0.5

0.75

1 0

Normalized Pore-Air Pressure

0

0.25

)

0.5

0.75

1 0

Normalized Pore-Water Pressure

0

)

0

Depth Ratio z/H

Depth Ratio z/H

0.8 m

ka/kw = 10-1

ζa = ζw = 0.5

(a)

0.6 m

0.5

0.75

1.4 m

ka/kw = 10-1

r = 0.4 m

0.25

0.25 r = 0.4 m

0.6 m

0.8 m 1 m 1.8 m

0.5

1.4 m T = 10-3 ka/kw =

0.75

10-1

r = 0.4 m 0.25

1.8 m 1.4 m

0.5

1m 0.8 m

T = 10-3

0.75

ζa = ζw = 0.5

0.6 m

ka/kw = 10-1 ζa = ζw = 0.5

1

1

(b)

0

0.25

0.5

Normalized Pore-Air Pressure

0.75

1 0

)

0

0.25

0.5

Normalized Pore-Water Pressure

0.75

1 0

)

Fig. 20. Distribution of excess pore pressures along the vertical domain at T ¼ 103 at different radii under (a) PTIB and (b) PTPB boundary conditions while adopting fa ¼ fw ¼ 0:5 (for radial and vertical flows, linear initial condition).

the drain well. The effects of radial distance on the pressure isochrones become less significant when the radius considerably increases (i.e. r P 1 m). 5. Conclusions This paper provides an analytical solution to predict the axisymmetric consolidation of unsaturated soils using variables separation and Laplace transformation techniques. The mathemat-

ical method captures the uniformly and linearly distributed excess pore pressures as initial conditions. The 3D Cartesian governing equations were first transformed into the polar equations. Subsequently, final solutions predicting excess pore pressure dissipation rates were obtained using Fourier Bessel and sine series, and Laplace transforms. In this study, two examples have been presented to illustrate the capabilities of the proposed analytical approach. The axisymmetric consolidation with only radial flow has been presented in

116

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

Example 1, whereas the case considering both radial and vertical flows has been investigated in Example 2. Overall, both examples show the significant effects of the air to water permeability ratio ðka =kw Þ on the dissipation process, in which changes in the ratio result in a single inverse S curves for the excess pore-air pressure, and double inverse S curves for the excess pore-water pressures and degree of consolidation. Besides, the excess pore pressures tend to dissipate faster for the points closer to the drain well. In Example 2, considering the uniform distribution of initial excess pore pressures, it can be concluded that the radial flow governs the axisymmetric consolidation while effects of vertical flow are less notable. It should also be noted that the compression process is ascribed to the increase in matric suction at the later stages of consolidation. On the other hand, the linear distribution of initial excess pore pressures with depth has considerable impacts on the dissipation process. In particular, the rapid reduction in both initial excess pore pressures with depth would lead to reduced average excess pore-air and pore-water pressures at the beginning of consolidation. Furthermore, the effects of initial pore pressures cause the pressure redistribution phenomenon in which, under the applied loads, some excess pore pressures are transferred toward the base of the soil, resulting in notable increases in the excess pore pressures. Appendix A. Polar transformation for x- and y-coordinates The net flux of air and water per unit volume of the element in Cartesian coordinates ðx; y; zÞ can be presented as follows: (a) Volume change in the 2D Cartesian system:

@



DV a V0



¼

@t

@



DV w V0

"

RH @ 2 ua kax 0 @x2 gM ua þ uatm
 nð1  Sr Þ @ua  0 @t ua þ uatm 

 ¼

@t

!

1

" kwx

cw

@ 2 uw @x2

! þ kwy

þ kay

@ 2 ua @y2

!#

Differentiating the equation of the radius r in Eq. (A.3) with respect to x and y results in:

@r x x ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ¼ cos h 2 2 @x r x þy

ðA:5aÞ

@r y y ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ¼ sin h @y x2 þ y 2 r

ðA:5bÞ

Then, differentiating the equation of the polar angle h in Eq. (A.4) with respect to x and y gives:

" !# @h @ x y y sin h 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 2 ¼ cos ¼ 2¼ 2 2 2 @x @x x r r þ y x þy

ðA:6aÞ

" !# @h @ y x x cos h 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ sin ¼ ¼ ¼ 2 @y @y x þ y2 r 2 r x2 þ y 2

ðA:6bÞ

Considering excess pore pressures ua and uw as functions of both r and h: (a) In the x-direction:

@ua @ua @r @ua @h @ua 1 @ua þ ¼ ¼ cos h  sin h r @h @x @r @x @h @x @r

ðA:7aÞ

@uw @uw @r @uw @h @uw 1 @uw þ ¼ ¼ cos h  sin h r @h @x @r @x @h @x @r

ðA:7bÞ

(b) In the y-direction:

@ua @ua @r @ua @h @ua 1 @ua þ ¼ ¼ sin h þ cos h r @h @y @r @y @h @y @r

ðA:8aÞ

@uw @uw @r @uw @h @uw 1 @uw þ ¼ ¼ sin h þ cos h r @h @y @r @y @h @y @r

ðA:8bÞ

ðA:1aÞ

@ 2 uw @y2

Taking the second derivatives of Eqs. (A.7) and (A.8) with respect to x and y, respectively, gives:

!# ðA:1bÞ

(b) Volume change in the 3D Cartesian system:   " ! ! !# @ DVV0a RH @ 2 ua @ 2 ua @ 2 ua  ¼ kax þ ka y þ ka z @t @x2 @y2 @z2 gM u0a þ uatm
 nð1  Sr Þ @ua  0 @t ua þ uatm

(a) In the x-direction:



 @ 2 ua @ 1 @ @ua 1 @ua  sin h ¼ cos h cos h  sin h @r r @h r @h @x2 @r ðA:9aÞ

 @ 2 uw @ 1 @ @uw 1 @uw ¼ cos h  sin h cos h  sin h 2 @r r @h r @h @x @r ðA:9bÞ

ðA:2aÞ

@



DV w V0



@t

¼

1

cw

"

2

kwx

@ uw @x2

!

2

þ kwy

@ uw @y2

!

2

þ kwz

@ uw @z2

(b) In the y-direction:

!#



 @ 2 ua @ 1 @ @ua 1 @ua þ cos h ¼ sin h sin h þ cos h @r r @h r @h @y2 @r ðA:10aÞ

ðA:2bÞ Cartesian coordinates ðx; yÞ can be transformed into the following set of polar coordinates ðr; hÞ:

(



pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ y2  h ¼ tan1 yx

r¼

ðA:3Þ

Also,

x ¼ r cos h y ¼ r sin h

ðA:4Þ



 @ 2 uw @ 1 @ @uw 1 @uw þ cos h ¼ sin h sin h þ cos h @r r @h r @h @y2 @r ðA:10bÞ Thus,

! @ 2 ua @ 2 ua @ 2 ua sin h cos h @ua sin h cos h 2  ¼ cos h  2 r r @x2 @r2 @h@r @h 2

þ

2

@ua sin h @ 2 ua sin h þ r @r @h2 r2

ðA:11aÞ

117

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

! @ 2 uw @ 2 uw @ 2 uw sin h cos h @uw sin h cos h 2  ¼ cos h  2 r r @x2 @r 2 @h@r @h 2

þ

Following Fredlund et al. [38], the 3D constitutive equations linking the soil deformation and stress state variables can be expressed in the polar coordinates ðr; h; zÞ, as presented below:

2

@uw sin h @ 2 uw sin h þ r r2 @r @h2

ðA:11bÞ

@

and

@ 2 ua @ 2 ua @ 2 ua sin h cos h @ua sin h cos h 2  ¼ sin h þ 2 2 2 r r @y @r @h@r @h þ

!

@

DV a V0

ðA:12aÞ

! @ 2 uw @ 2 uw @ 2 uw sin h cos h @uw sin h cos h 2  ¼ sin h þ 2 r r @y2 @r 2 @h@r @h

DV w V0

ðA:12bÞ

ðA:13aÞ

@ 2 uw @ 2 uw @ 2 uw 1 @uw 1 @ 2 uw þ ¼ þ þ 2 r @r r @h2 @x2 @y2 @r 2

ðA:13bÞ

@ 2 ua @ 2 ua @ 2 ua 1 @ua þ ¼ þ r @r @x2 @y2 @r 2

ðA:14aÞ

@ 2 uw @ 2 uw @ 2 uw 1 @uw þ ¼ þ r @r @x2 @y2 @r 2

ðA:14bÞ

Assuming kax ¼ kay ¼ kar and kwx ¼ kwy ¼ kwr , Eq. (A.14) can be substituted into Eqs. (A.1) and (A.2): (a) Volume change in radial direction only:

 k RH @ 2 ua 1 @ua  ar þ ¼ r @r @r 2 gM u0a þ uatm
 nð1  Sr Þ @ua  0 @t ua þ uatm

@t

@



DV w V0

 ¼

@t

@ 2 uw 1 @uw þ r @r @r 2

k wr

cw

!

ðA:15aÞ



DV a V0



"

ðA:15bÞ

RH @ ua 1 @ua  kar ¼ þ r @r @r 2 gM u0a þ uatm
 nð1  Sr Þ @ua  0 @t ua þ uatm

@t

2

!

2

þ k az

@ ua @z2

!#



DV w V0

@t

 ¼

1

cw

" kwr

ðA:17bÞ

!
 @uw @ua @ 2 uw 1 @uw w ¼0 þ þ Cw þ cv r r @r @t @t @r 2

ðA:18aÞ

ðA:18bÞ

(b) Polar governing equations considering both radial and vertical flows:

! !
 @ua @uw @ 2 ua 1 @ua @ 2 ua a a þ cv z ¼0 þ þ Ca þ cv r r @r @t @t @r2 @z2

ðA:19aÞ ! !
 @uw @ua @ 2 uw 1 @uw @ 2 uw w þ c ¼0 þ þ Cw þ cw v r @r 2 v z @z2 r @r @t @t ðA:19bÞ Eqs. (A.18) and (A.19) present the governing flow equations in the axisymmetric condition. Simplified forms of these equations can be found in Eqs. (1) and (2). Note that the consolidation coefw ficients for air (i.e. C a ; cav r and cav z ) and water (i.e. C w ; cw v r and cv z ) phases are defined in Eq. (3). Appendix B. Solutions for excess pore pressure dissipation and average degree of consolidation

@ 2 uw 1 @uw þ r @r @r 2

! þ kwz

@ 2 uw @z2

Example 1 (Axisymmetric Consolidation with Radial Flow only). In this example, during the consolidation process, pore-air and porewater dissipate in the radial direction only. Thus, initial conditions can be presented as follows:

ua ðr; 0Þ ¼ u0a

ðB:1aÞ

uw ðr; 0Þ ¼ u0w

ðB:1bÞ

Constant terms Na and Nw ði ¼ 0; 1; 2 . . .Þ then become:

ðA:16aÞ @

 @ rr þ rh þ rz @ ðua  uw Þ  ua þ mw 2 @t @t 3

¼ mw 1

!

(b) Volume change in radial and vertical directions:

@

ðA:17aÞ

!
 @ua @uw @ 2 ua 1 @ua þ cav r ¼0 þ þ Ca r @r @t @t @r2

For the case of axial symmetry, the pore pressures ua and uw are independent of h, thus,

DV a V0

 @ rr þ rh þ rz @ ðua  uw Þ  ua þ ma2 @t @t 3

(a) Polar governing equations considering the radial flow only:

@ 2 ua @ 2 ua @ 2 ua 1 @ua 1 @ 2 ua þ ¼ þ þ 2 r @r r @h2 @x2 @y2 @r 2



¼ ma1

where rr ; rh and rz are the total stress in r-, h- and z-directions, respectively. Under the constant loading condition, let set @ rr =@t ¼ @ rh =@t ¼ @ rz =@t ¼ 0. Eqs. (A.15) and (A.16) can be equated to Eq. (A.17), giving:

Adding Eqs. (A.11a) and (A.12a) as well as Eqs. (A.11b) and (A.12b) yields:

@



@t   @t

@ua cos2 h @ 2 ua cos2 h þ r @r @h2 r 2

@uw cos2 h @ 2 uw cos2 h þ þ r r2 @r @h2



Na ¼

2X i 0 qffiffiffiffi ua  i ni

ðB:2aÞ

Nw ¼

2X i 0 qffiffiffiffi uw  i ni

ðB:2bÞ

!# ðA:16bÞ

118

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

By repeating the analytical procedure in Eqs. (9)–(26) while discarding the permeability coefficients in the domain z (i.e. kaz and kwz ) and the Fourier sine series, solutions predicting excess poreair and pore-water pressure dissipation for Eq. (1) can be obtained and then combined with Eq. (B.2) become:   9 8  i i i i a1 n 2 t a2 n 2 t a1 n 2 t a2 n 2 t > 0qffiffiffiffi 1> ðr w Þ ðr w Þ ðr w Þ ðr w Þ > > X  e W þ e þ i 1 e 1 e i < = n X ua ðr; tÞ ¼ qffiffiffiffi Di @ rA > > g rw > > : ;  i ni

Na ¼

4X i Zj qffiffiffiffi pffiffiffiffiffiffi u0a j  i ni l

ðB:7aÞ

Nw ¼

4X i Zj qffiffiffiffi pffiffiffiffiffiffi u0w j  i ni l

ðB:7bÞ

where

Zj ¼ 1

ðB:3aÞ

¼ 1  ð1Þ j

  9 8  i i i i a2 n 2 t a1 n 2 t a2 n 2 t > 0qffiffiffiffi 1> 0 a1 ðr n Þ2 t 0 ðr w Þ ðr w Þ ðr w Þ w > > X e  e W e þ e þ i = < 1 1 n X uw ðr;t Þ ¼ qffiffiffiffi Di @ rA > > rw g > > ; :  i ni

for the PTIB boundary condition; or for the PTPB boundary condition: ðB:8Þ

i

Solutions predicting excess pore-air and pore-water pressure dissipation can be further expanded from Eq. (26) while combining with Eq. (B.7), resulting in: (i) The PTIB boundary condition:

ðB:3bÞ

where

h

g ¼ cav r  cwv r

i12

2

þ 4cav r cw vr CwCa ;
a 
 w 1 cv r þ cv r þ g 1 cav r þ cw vr  g ; a2 ¼ ; a1 ¼ 2 2 1  CwCa 1  CwCa  a X1 ¼ cv r  cwv r u0a  2cwv r C a u0w ; W1 ¼ gu0a ;  X01 ¼ cwv r  cav r u0w  2cav r C w u0a ; and W01 ¼ gu0w :



ms2  ms1

( R re 2p

rw

ua ðr;t Þrdr

p½ðre Þ2 ðrw Þ2 



¼

)  u0a

( R re  ms2

2p

rw

uw ðr;t Þrdr

p½ðre Þ2 ðrw Þ2 

ðB:4Þ

 u0w

rw

h i 0qffiffiffiffi 1 pffiffiffiffiffiffi ! j 1 X ni 2X i 1  ð1Þ lj qffiffiffiffi pffiffiffiffiffiffi Di @ ua ðr; z; t Þ ¼ r A sin z j rw H l j¼0  i ni   ij 3 2  ij ij ij X2 ea1 t  ea2 t þ W2 ea1 t þ ea2 t 5 ðB:10aÞ 4 ij

   9 ni ni ni ni  2 8 >X ea1 ðrw Þ2 t  ea2 ðrw Þ2 t þ W ea1 ðrw Þ2 t þ ea2 ðrw Þ2 t > > 1 = < 1 Xi > ; ua ðr; t Þrdr ¼ i i > g n > > > ; :

and Z

re

uw ðr; tÞrdr ¼

   9 ni ni ni ni  2 8 >X0 ea1 ðrw Þ2 t  ea2 ðrw Þ2 t þ W0 ea1 ðrw Þ2 t þ ea2 ðrw Þ2 t > > = < 1 1 Xi >

g

 i ni > > :

rw

> > ;

g

h i 0qffiffiffiffi 1 pffiffiffiffiffiffi ! j 1 X ni 2X i 1  ð1Þ lj qffiffiffiffi pffiffiffiffiffiffi Di @ r A sin uw ðr; z; tÞ ¼ z H rw lj j¼0  i ni   ij 3 2 0  ij ij ij X2 ea1 t  ea2 t þ W02 ea1 t þ ea2 t 5 ðB:10bÞ 4 ij

:

g

ðB:6Þ

where   X2 ¼ cav kija  cwv kijw u0a  2cwv C a kijw u0w ; W2 ¼ gij u0a ;   X02 ¼ cwv kijw  cav kija u0w  2cav C w kija u0a ; and W02 ¼ gij u0w :

Example 2 (Axisymmetric Consolidation with both Radial and Vertical Flows). (a) Uniform Initial Condition ðfa ¼ fw ¼ 0Þ When fa ¼ fw ¼ 0, the constant terms Na and Nw ði ¼ 0; 1; 2 . . . ; j ¼ 0; 1; 2 . . .Þ, as presented in Eq. (23), can be obtained as follows:

  
ev ðtÞ  ¼ U ¼  2
ev 2 ð1Þ



ms2



ms1



( R H R re 2p

0

p½

ua ðr;z;t Þrdrdz rw ðr e Þ2 ðr w Þ2 H

) ms2

2p

0

ðB:11Þ

The average degree of consolidation U capturing the uniform initial condition can be further expressed from Eq. (30):

( R H R re

  p½   s s s 0 0 m1  m2 ua þ m2 uw u0a

ðB:9bÞ

(ii) The PTPB boundary condition:

ms1  ms2 u0a þ ms2 u0w

where re

0qffiffiffiffi 1 pffiffiffiffiffiffi ! 1 X ni 2X i 1 lj i@ qffiffiffiffi pffiffiffiffiffiffi D r A sin uw ðr; z; tÞ ¼ z j H rw j¼0  i ni l   ij 3 2 0  ij ij ij X2 ea1 t  ea2 t þ W02 ea1 t þ ea2 t 5 4

gij

)

ðB:5Þ

Z

ðB:9aÞ

g

In addition, the average degree of consolidation U can be presented as below:

  
ev ðtÞ   U ¼  1
ev 1 ð1Þ

0qffiffiffiffi 1 pffiffiffiffiffiffi ! 1 X ni 2X i 1 lj i@ A q ffiffiffiffi p ffiffiffiffiffi ffi D ua ðr; z; t Þ ¼ r sin z j rw H j¼0  i ni l   ij 3 2  ij ij ij X2 ea1 t  ea2 t þ W2 ea1 t þ ea2 t 5 4 ij

uw ðr;z;tÞrdrdz rw ðr e Þ2  ðr w Þ2 H



) 

u0w ðB:12Þ

119

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

where (i) The PTIB boundary condition: Z

H

0

Z

re

ua ðr;z;t Þrdrdz ¼

 2 i 1 2 X X H

 i ni l j  ij   ij 3 ij ij X2 ea1 t  ea2 t þ W2 ea1 t þ ea2 t 5 ðB:13aÞ 4

rw

H

0

Z

re

uw ðr;z; tÞrdrdz ¼

 2 i 1 2 X X H

 i ni l j  ij   ij 3 ij ij X02 ea1 t  ea2 t þ W02 ea1 t þ ea2 t 4 5ðB:13bÞ  ij

rw

2

(ii) The PTPB boundary condition:

0qffiffiffiffi 1 pffiffiffiffiffiffi ! 1 X ni 2X i 1 lj i@ qffiffiffiffi pffiffiffiffiffiffi D ua ðr; z; t Þ ¼ r A sin z j H rw j¼0  i ni l   ij 3 2  ij ij ij X4 ea1 t  ea2 t þ W4 ea1 t þ ea2 t 5 ðB:18aÞ 4 ij

g

0qffiffiffiffi 1 pffiffiffiffiffiffi ! 1 X ni 2X i 1 lj i@ qffiffiffiffi pffiffiffiffiffiffi D uw ðr; z; t Þ ¼ r A sin z j rw H j¼0  i ni l   ij 3 2 0  ij ij ij X4 ea1 t  ea2 t þ W04 ea1 t þ ea2 t 5 ðB:18bÞ 4 ij

gij

0

Z

re

ðB:17bÞ

g

j¼0

(ii) The PTPB boundary condition:  2 h i2 j i Z H Z re 1 2 X H 1  ð1Þ X ua ðr;z;t Þrdrdz ¼ lj  i ni 0 rw j¼0   ij 3 2  ij ij ij X2 ea1 t  ea2 t þ W2 ea1 t þ ea2 t 5 ðB:14aÞ 4

H

0qffiffiffiffi 1 pffiffiffiffiffiffi ! i i 1 X 2X 1 i @ n A lj qffiffiffiffi j D r sin uw ðr; z; t Þ ¼ z r l H w j¼0  i ni   ij 3 2 0  ij ij ij X3 ea1 t  ea2 t þ W03 ea1 t þ ea2 t 5 4 ij

2

g

Z

ðB:17aÞ

gij

j¼0

gij

Z

0qffiffiffiffi 1 pffiffiffiffiffiffi ! i 1 X 2X i 1 i @ n A lj qffiffiffiffi j D ua ðr; z; t Þ ¼ r sin z l H r w j¼0  i ni   ij 3 2  ij ij ij X3 ea1 t  ea2 t þ W3 ea1 t þ ea2 t 5 4

uw ðr;z; tÞrdrdz ¼

 2 h i2 j i 1 2 X H 1  ð1Þ X

lj  i ni  ij   ij 3 ij ij X02 ea1 t  ea2 t þ W02 ea1 t þ ea2 t 5ðB:14bÞ 4 ij

rw

j¼0

2

g

g

where



(b) Linear Initial Condition ðfa ; fw > 0Þ When both fa > 0 and fw > 0, the constant terms Na and Nw ði ¼ 0; 1; 2 . . . ; j ¼ 0; 1; 2 . . .Þ, as presented in Eq. (23), are as follows: (i) The PTIB boundary condition:

qffiffiffiffiffiffi

Na ¼

4X i 1 qffiffiffiffi pffiffiffiffiffiffi u0a i i lj  n

Nw ¼

4X i 1 qffiffiffiffi pffiffiffiffiffiffi u0w i i lj  n

l j  ð1Þ j fa

qffiffiffiffiffiffi

X3 ¼ cav kija  cwv kijw

ðB:15aÞ  ðB:15bÞ

(ii) The PTPB boundary condition:

Na ¼

h i 4X i 1 qffiffiffiffi pffiffiffiffiffiffi u0a 1 þ ð1Þ j ðfa  1Þ j  i ni l

ðB:16aÞ

Nw ¼

h i 4X i 1 qffiffiffiffi pffiffiffiffiffiffi u0w 1 þ ð1Þ j ðfw  1Þ j  i ni l

ðB:16bÞ

Solutions predicting excess pore-air and pore-water pressure dissipation can be expanded from Eq. (26) while combining with Eqs. (B.15) and (B.16), giving: (i) The PTIB boundary condition:

  
ev ðtÞ  ¼ U ¼  3
ev 3 ð1Þ



ms2  ms1



( R H R re 2p

0

p½

ua ðr;z;t Þrdrdz rw ðre Þ2 ðr w Þ2 H





  u0a 1  f2a

)

2p

ðB:19Þ

The average degree of consolidation U capturing the linear initial condition is expressed as follows:

( R H R re  ms2



l j  ð1Þ j fa u0a

qffiffiffiffiffiffi  ij  2cw l j  ð1Þ j fw u0w ; v C a kw qffiffiffiffiffiffi  W3 ¼ gij l j  ð1Þ j fa u0a ;   qffiffiffiffiffiffi X03 ¼ cwv kijw  cav kija l j  ð1Þ j fw u0w qffiffiffiffiffiffi   2cav C w kija l j  ð1Þ j fa u0a ; qffiffiffiffiffiffi  W03 ¼ gij l j  ð1Þ j fw u0w ; i  h X4 ¼ cav kija  cwv kijw 1 þ ð1Þ j ðfa  1Þ h i j ij 0 u0a  2cw v C a kw 1 þ ð1Þ ðfw  1Þ uw ; h i W4 ¼ gij 1 þ ð1Þ j ðfa  1Þ u0a ; i  h X04 ¼ cwv kijw  cav kija 1 þ ð1Þ j ðfw  1Þ h i u0w  2cav C w kija 1 þ ð1Þ j ðfa  1Þ u0a ; and h i W04 ¼ gij 1 þ ð1Þ j ðfw  1Þ u0w :



l j  ð1Þ j fw

qffiffiffiffiffiffi

0

uw ðr;z;t Þrdrdz rw ðr e Þ2 ðrw Þ2 H

p½   f ms1  ms2 u0a 1  2a þ ms2 u0w 1  f2w



  u0w 1  f2w

) ðB:20Þ

120

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121

where h is the volumetric coefficient of solubility; and

where (i) The PTIB boundary condition:  2 i Z H Z re 1 2 X X H ua ðr;z;t Þrdrdz ¼ 3  i ni ðl j Þ2 0 rw j¼0    ij 3 2 ij ij ij X3 ea1 t  ea2 t þ W3 ea1 t þ ea2 t 5 4

bw ¼ 

1 dV w ; is the compressibility of water: V w duw

Combining Eqs. (C.2a) and (C.3a) and Eqs. (C.2b) and (C.3b), resulting in:

gij

ms1 dðrr  ua Þ þ ms2 dðua  uw Þ ¼

ð1  Sr þ hSr Þn dua þ Sr nbw duw ua þ uatm ðC:5aÞ

ma1 dðrr  ua Þ þ ma2 dðua  uw Þ ¼

ð1  Sr þ hSr Þn dua ua þ uatm

ðB:21aÞ Z

H

0

Z

re

uw ðr;z;t Þrdrdz ¼

 2 i 1 2 X X H

 i ni ðl j Þ2   ij 3 2 0  ij ij ij X3 ea1 t  ea2 t þ W03 ea1 t þ ea2 t 5 4

rw

3

j¼0

g

ij

ðB:21bÞ

(ii) The PTPB boundary condition:  2 h i j i Z H Z re 1 2 X H 1  ð1Þ X ua ðr;z;t Þrdrdz ¼ lj  i ni 0 rw j¼0   ij 3 2  ij ij ij X4 ea1 t  ea2 t þ W4 ea1 t þ ea2 t 4 5  ij

g

0

H

Z

re

uw ðr;z;t Þrdrdz ¼

rw

ðC:5bÞ

Simultaneously solving Eqs. (C.5a) and (C.5b) to determine the changes in excess pore-air and pore-water pressures in response to the change in total stress, thus,

Ba ¼

dua Rs2 Ra1  Ra2 ¼ drr 1  Rs1 Ra1

ðC:6aÞ

Bw ¼

duw Rs2  Rs1 Ra2 ¼ drr 1  Rs1 Ra1

ðC:6bÞ

where ðB:22aÞ

Z

ðC:4Þ

 2 h i j i 1 2 X H 1  ð1Þ X

lj  i ni   ij 3 2 0  ij ij ij X4 ea1 t  ea2 t þ W04 ea1 t þ ea2 t 5 4

Rs1 ¼ Ra1

j¼0

¼

r þhSr Þn ms2  ms1  ð1S ua þuatm

ms2 ma2



þ nSr bw ma2

ma1



ð1Sr þhSr Þn ua þuatm

ms1 ; þ nSr bw

;

Rs2 ¼

;

and Ra2 ¼

ms2

ma1 ma2



ma1

r þhSr Þn  ð1S ua þuatm

: ðC:7Þ

gij

ðB:22bÞ

Appendix C. Initial excess pore pressures in response to constant loading Following Fredlund et al. [38], the 3D constitutive equation for the soil structure and air phase in the polar coordinates ðr; h; zÞ can be expressed considering the total stress increments drr ; drh and drz as follows:

r þ r þ r  dV r h z ¼ ms1 d  ua þ ms2 dðua  uw Þ V0 3

ðC:1aÞ

r þ r þ r  dV a r h z ¼ ma1 d  ua þ ma2 dðua  uw Þ V0 3

ðC:1bÞ

Under an isotropic loading condition, total stress increments in the polar system would be equal (i.e., drr ¼ drh ¼ drz ), thus Eq. (C.1) becomes:

dV ¼ ms1 dðrr  ua Þ þ ms2 dðua  uw Þ V0

ðC:2aÞ

dV a ¼ ma1 dðrr  ua Þ þ ma2 dðua  uw Þ V0

ðC:2bÞ

The volume changes, dV and dV a , due to pore fluid compression and air compression, respectively, can be expressed as:

dV ð1  Sr þ hSr Þn ¼ dua þ Sr nbw duw V0 ua þ uatm

ðC:3aÞ

dV a ð1  Sr þ hSr Þn ¼ dua ua þ uatm V0

ðC:3bÞ

The terms Ba and Bw are the pore pressure parameters with respect to the air and water phases, respectively, under the isotropic loading condition. Since terms Rs1 ; Ra1 and Ra2 contain the absolute pore-air pressure, an iterative technique is required to solve for dua and duw . Referring to Fredlund et al. [38], it is worth noting that suction (i.e. ua  uw ) reduces with the increasing isotropic stress. As the stress increases significantly, the parameters Ba and Bw would approach 1 and subsequently suction becomes zero. References [1] Terzaghi K. Theoretical soil mechanics. New York: John Wiley and Sons; 1943. [2] Barron RA. Consolidation of fine-grained soils by drain wells. Trans Am Soc Civil Eng 1948;113:718–42. [3] Yoshikuni H, Nakanodo H. Consolidation of soils by vertical drain wells with finite permeability. Soils Found 1974;14:35–46. [4] Hansbo S. Consolidation of fine-grained soils by prefabricated drains. In: Proceedings of the 10th international conference on soil mechanics and foundation engineering, vol. 3; 1981. p. 677–82. [5] Tang XW, Onitsuka K. Consolidation by vertical drains under time-dependent loading. Int J Numer Analyt Meth Geomech 2000;24:739–51. [6] Leo CJ. Equal strain consolidation by vertical drains. J Geotech Geoenviron Eng 2004;130:316–27. [7] Zhu G, Yin JH. Consolidation analysis of soil with vertical and horizontal drainage under ramp loading considering smear effects. Geotext Geomemb 2004;22:63–74. [8] Conte E, Troncone A. Radial consolidation with vertical drains and general time-dependent loading. Can Geotech J 2009;46:25–36. [9] Geng X, Indraratna B, Rujikiatkamjorn C. Analytical solutions for a single vertical drain with vacuum and time-dependent surcharge preloading in membrane and membraneless systems. Int J Geomech 2010;12:27–42. [10] Lu M, Wang S, Sloan SW, Indraratna B, Xie K. Nonlinear radial consolidation of vertical drains under a general time-variable loading. Int J Numer Analyt Meth Geomech 2015;39:51–62. [11] Lei GH, Zheng Q, Ng CWW, Chiu ACF, Xu B. An analytical solution for consolidation with vertical drains under multi-ramp loading. Géotechnique 2015;65:531–47. [12] Indraratna B, Rujikiatkamjorn C, Sathananthan I. Radial consolidation of clay using compressibility indices and varying horizontal permeability. Can Geotech J 2005;42:1330–41.

L. Ho et al. / Computers and Geotechnics 74 (2016) 102–121 [13] Walker R, Indraratna B. Vertical drain consolidation with parabolic distribution of permeability in smear zone. J Geotech Geoenviron Eng 2006;132:937–41. [14] Walker R, Indraratna B. Vertical drain consolidation with overlapping smear zones. Géotechnique 2007;57:463–7. [15] Wang XS, Jiao JJ. Analysis of soil consolidation by vertical drains with double porosity model. Int J Numer Analyt Meth Geomech 2004;28:1385–400. [16] Deng YB, Xie KH, Lu MM. Consolidation by vertical drains when the discharge capacity varies with depth and time. Comp Geotech 2013;48:1–8. [17] Wu H, Hu L. Analytical solution for axisymmetric electro-osmotic consolidation. Géotechnique 2013;63:1074–9. [18] Fredlund DG. The scope of unsaturated soil mechanics: an overview. In: Proceedings of the 1st international conference on unsaturated soils, vol. 3; 1996. p. 1155–77. [19] Fredlund DG, Hasan JU. One-dimensional consolidation theory: unsaturated soils. Can Geotech J 1979;16:521–31. [20] Dakshanamurthy V, Fredlund DG. Moisture and air flow in an unsaturated soil. In: Proceedings of the 4th international conference on expansive soils, vol. 1; 1980. p. 514–32. [21] Darkshanamurthy V, Fredlund DG, Rahardjo H. Coupled three-dimensional consolidation theory of unsaturated porous media. In: Proceedings of the 5th international conference on expansive soils; 1984. p. 99–103. [22] Qin AF, Chen GJ, Tan YW, Sun DA. Analytical solution to one-dimensional consolidation in unsaturated soils. Appl Math Mech 2008;29:1329–40. [23] Qin AF, Sun DA, Tan YW. Analytical solution to one-dimensional consolidation in unsaturated soils under loading varying exponentially with time. Comp Geotech 2010;37:233–8. [24] Shan Z, Ling D, Ding H. Exact solutions for one-dimensional consolidation of single-layer unsaturated soil. Int J Numer Analyt Meth Geomech 2012;36:708–22. [25] Zhou WH, Zhao LS, Li XB. A simple analytical solution to one-dimensional consolidation for unsaturated soils. Int J Numer Analyt Meth Geomech 2014;38:794–810. [26] Zhou WH, Zhao LS. One-dimensional consolidation of unsaturated soil subjected to time-dependent loading with various initial and boundary conditions. Int J Geomech 2014;14:291–301. [27] Ho L, Fatahi B, Khabbaz H. Analytical solution for one-dimensional consolidation of unsaturated soils using eigenfunction expansion method. Int J Numer Analyt Meth Geomech 2014;38:1058–77. [28] Shan Z, Ling D, Ding H. Analytical solution for the 1D consolidation of unsaturated multi-layered soil. Comp Geotech 2014;57:17–23. [29] Qin A, Sun DA, Zhang J. Semi-analytical solution to one-dimensional consolidation for viscoelastic unsaturated soils. Comp Geotech 2014;62:110–7.

121

[30] Conte E. Consolidation analysis for unsaturated soils. Can Geotech J 2004;41:599–612. [31] Ho L, Fatahi B. Analytical solution for the two-dimensional plane strain consolidation of an unsaturated soil stratum subjected to time-dependent loading. Comp Geotech 2015;67:1–16. [32] Ho L, Fatahi B, Khabbaz H. A closed form analytical solution for twodimensional plane strain consolidation of unsaturated soil stratum. Int J Numer Analyt Meth Geomech 2015;39:1665–92. [33] Conte E. Plane strain and axially symmetric consolidation in unsaturated soils. Int J Geomech 2006;6:131–5. [34] Qin AF, Sun DA, Yang LP, Weng YF. A semi-analytical solution to consolidation of unsaturated soils with the free drainage well. Comp Geotech 2010;37:867–75. [35] Zhou WH, Tu S. Unsaturated consolidation in a sand drain foundation by differential quadrature method. Proc Earth Planet Sci 2012;5:52–7. [36] Zhou WH. Axisymmetric consolidation of unsaturated soils by differential quadrature method. Math Prob Eng 2013;2013:1–14. [37] Crump KS. Numerical inversion of Laplace transforms using a Fourier series approximation. J ACM 1976;23:89–96. [38] Fredlund DG, Rahardjo H, Fredlund MD. Unsaturated soil mechanics in engineering practice. New Jersey: John Wiley and Sons; 2012. [39] Craig RF. Craig’s soil mechanics. New York: Taylor and Francis; 2004. [40] Venkatramaiah C. Geotechnical engineering. New Delhi: New Age International Publishers; 2006. [41] Coduto DP, Yeung MR, Kitch WA. Geotechnical engineering: principles and practices. New Jersey: Pearson Higher Education; 2011. [42] Rowe RK. Geotechnical and geoenvironmental engineering handbook. Boston: Kluwer; 2001. [43] Aysen A. Soil mechanics: basic concepts and engineering applications. Amsterdam: A. A. Balkema; 2002. [44] Aysen A. Soil mechanics: basic concepts and engineering applications. Amsterdam: A.A. Balkema; 2002. [45] Holtz RD. Preloading with prefabricated vertical strip drains. Geotext Geomemb 1987;6:109–31. [46] Holtz RD, Jamiolkowski MB, Lancellotta R, Pedroni R. Prefabricated vertical drains: design and performance. Boston: Butterworth; 1991. [47] Smoltczyk U. Geotechnical engineering handbook. Berlin: Earnst and Sohn; 2003. [48] Mandel J. Consolidation Des Sols. Géotechnique 1953;3:287–99. [49] Cryer C. A comparison of the three-dimensional consolidation theories of Biot and Terzaghi. Quart J Mech Appl Math 1963;16:401–12. [50] Wong TT, Fredlund DG, Krahn J. A numerical study of coupled consolidation in unsaturated soils. Can Geotech J 1998;35:926–37.