# Axisymmetric Consolidation of Unsaturated Soils by Differential Quadrature Method

Authors: Wan-Huan Zhou
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Axisymmetric consolidation in a sand drain foundation is a common problem in foundation engineering. In unsaturated soils, the excess pore-water and pore-air pressures simultaneously change during the consolidation procedure; and the solutions are not easy to obtain. The present paper uses the differential quadrature method (DQM) for axisymmetric consolidation of unsaturated soils in a sand drain foundation. The radial seepage of sand drain foundation is considered based on the framework of Fredlund’s one-dimensional consolidation theory in unsaturated soils. With the use of Darcy’s law and Fick’s law, the polar governing equations of excess pore-air and pore-water pressures of axisymmetric consolidation are derived. By using DQM, the two governing equations are transformed into two sets of ordinary differential equations. Then the solutions of excess pore-water and pore-air pressures can be obtained by Rong-Kutta method. The DQM solution can be used to deal with the case of nonuniform initial pore-air and pore-water distributions. Finally, case studies are presented to investigate the behavior of axisymmetric consolidation of unsaturated soils. The convergence analysis and average degree of consolidation, the settlements in radial and vertical direction, and the effects of different initial excess pore pressure distributions are presented, and discussed in this paper.


---

## Body

## 1. Introduction

As a common phenomenon of geotechnical engineering, consolidation is a process of decreasing soil volume when soil is subjected to increased stress. Terzaghi [1] established a classical theory for the analysis of one-dimensional (1-D) consolidation in saturated soils, which is still widely used in engineering practice. But in real cases, soils related to engineering are usually in a state of unsaturation and are defined by more characteristics than saturated soils, such as the coupled dissipation of pore-water and pore-air pressures during consolidation. In the case of unsaturated soil, the excess pore-water pressure and pore-air pressure change simultaneously during the consolidation process.For the consolidation of unsaturated soil, Biot [2] derived a general theory of consolidation for unsaturated soil which had occluded air bubbles. The theory provides a coupling between the magnitude and progress of settlement. Blight [3] presented a consolidation equation for the air phase, where the air was in a dry, rigid, and unsaturated soil. Based on Biot’s theory of consolidation, Scott [4] derived a consolidation equation for unsaturated soils with occluded air bubbles by changing void ratio and degree of saturation. Subsequently, Barden [5] presented an analysis for 1-D consolidation of compacted unsaturated clay. Fredlund and Hasan [6] proposed a 1-D consolidation theory for unsaturated soils. This formulation is based on two continuity partial differential equations, one for the water phase and the other for the air phase, which have to be solved simultaneously to give water and air pressures at any time and elevation throughout the soil. According to Fredlund’s consolidation theory [7], with the use of Darcy’s law and Fick’s law, Qin et al. [8] obtained a semianalytical solution for consolidation of unsaturated soils with free drainage well. They obtained the excess pore-air and pore-water pressures and the soil layer settlement in the Laplace transformed domains by utilizing the Bessel function. More recently, Zhou et al. [9] presented a simple analytical solution to 1-D consolidation for unsaturated soil, which can be degenerated into Terzaghi consolidation solution for fully saturated condition.As the consolidation equations of unsaturated soil are coupled, the analytical solutions cannot be easily obtained. Therefore, numerical solutions such as the finite element method (FEM) and the finite difference method (FDM) are traditionally employed. These two methods approximate the partial derivatives of a function at a grid point, only by using a limited number of function values in the vicinity of that grid point. The accuracy and stability of these methods depend on the size of the grid spacing. Compared with these two numerical methods, differential quadrature method (DQM) is a more efficient numerical method for the rapid solution of linear and nonlinear partial differential equations involving one dimension or multiple dimensions. The fundamental idea of DQM is that the derivatives of each node in a continuous function can be expressed as a weighted linear sum of function values at all grid points. Malik and Civan [10] have made a comprehensive comparison for linear and nonlinear convection-diffusion-reaction problems and have shown that the DQM is superior to FEM and FDM in numerical accuracy as well as computational efficiency. Many studies have successfully employed DQM for solving engineering problems [11, 12]. Examples for the applications of DQM in solving complex problems in geotechnical engineering are presented in these literatures [13–16] among others.In this paper, DQM is introduced into the analysis of axisymmetric consolidation of unsaturated soils in a sand drain foundation. Following Fredlund’s one-dimensional consolidation theory of unsaturated soils, with the use of Darcy’s law and Fick’s law, the polar governing equations of excess pore-air pressure and pore-water pressure of axisymmetric consolidation are obtained. The behaviors of axisymmetric consolidation of unsaturated soils in the sand drain foundation have been analyzed. The convergence analysis and average degree of consolidation, the effects of settlement in radial and vertical directions, the effects of different initial excess pore-air and pore-water pressure distributions, and the effects of different boundary conditions are presented and discussed.


## 2. Mathematical Model and Governing Equations

### 2.1. Mathematical Model for Axisymmetric Consolidation of Unsaturated Soils

As shown in Figure1, the unsaturated soil layer is described in a schematic diagram and thickness, H. Figures 1(a) and 1(b) are the sectional view and vertical view of the unsaturated soil layer, respectively. The radii of soil layer and sand drain are re and rw. A surcharge q is applied on the top surface of the soil layer. kw and ka are the permeability coefficients of water and air, respectively. The coordinated origin is selected at the center of the top surface and the z coordinate is positive downward. In this paper, the vertical prismatic blocks of soil surrounding the sand drains are simulated by cylindrical blocks, of radius S. Square pattern and triangular pattern are two common distributions of sand drains. For instance, if square pattern is adopted as shown in Figure 1(c), then S=0.564re. The basic assumptions are the same as those of Fredlund’s one-dimensional consolidation theory for unsaturated soils and the other assumptions are listed as follows.(1)
The seepage of water and air phases is considered only in radial direction.(2)
The top and bottom surfaces are impermeable and the right boundary surface is impermeable or impeded.(3)
The permeability coefficients of water and air phases are assumed to be constant.(4)
In the process of consolidation, the strain is small enough and the density of water is assumed to be constant.(a) The sectional view of sand drain foundation, (b) the vertical view of sand drain foundation, and (c) the distribution of sand drains (square pattern).
(a)
(b)
(c)


### 2.2. Governing Equations

The net flux of water through the soil layer is computed from the volume of water entering and leaving the soil layer within a period of time, with respect to Darcy’s law in the polar coordinate system:(1)∂(Vw/V0)∂t=kwγw(∂2uw∂r2+1r∂uw∂r),
where V0 is the initial total volume of the soil, ∂(Vw/V0)/∂t is the net flux of water per unit volume of the soil, uw is the excess pore-water pressure, γw is the water unit weight, and t is the time variable.The net flux of water per unit volume of the soil can be obtained by differentiating the water phase constitutive relation with respect to time:(2)∂(Vw/V0)∂t=m1kw∂(q-ua)∂t+m2w∂(ua-uw)∂t,
where m1kw is the coefficient of water volume change with respect to a change in the net normal stress (q-ua), m2w is the coefficient of water volume change with respect to a change in matric suction (ua-uw), and ua is the excess pore-air pressure.For constant loading,dq/dt=0, substituting (1) into (2), the governing equation for the water phase can be written as
(3)∂uw∂t=-Cw∂ua∂t-Cvw(∂2uw∂r2+1r∂uw∂r),
where Cw=(m1kw-m2w)/m2w and Cvw=kw/γwm2w.The net flux of air through the soil layer is computed from the volume of air entering and leaving the soil layer within a period of time, with respect to Fick’s law in the polar coordinate system:(4)∂(ρVa/V0)∂t=kag(∂2ua∂r2+1r∂ua∂r),
where ρ is the density of air phase, ∂(ρVa/V0)/∂t is the net mass rate of air flow per unit volume of the soil, and g is the gravitational acceleration.The total volume change of an unsaturated soil can be assumed to be small during the consolidation process. The volume of airVa can be related to the volume-mass properties of the soil: Va=(1-S0)n0V0, where n0 and Sr0 are the initial porosity and initial degree of saturation before loading. According to Boyle’s law, assuming there is no initial excess air pressure in the soil before loading, we have [17, 18]
(5)∂(Va/V0)∂t=kaRTgu-a0M(∂2ua∂r2+1r∂ua∂r)-uatmn0(1-S0)(u-a0)2∂ua∂t,
where R is the universal air constant (8.314J/molK), T=293.15K is the absolute temperature, M=29kg/kmol is the molecular mass of air, u¯0a is absolute pore-air pressure (i.e., u¯0a=u0a+uatm), and uatm=101.3kPa is atmospheric pressure.The derivative of the air phase constitutive relation with respect to time is equal to the net flux of air per unit volume of the soil:(6)∂(Va/V0)∂t=m1ka∂(q-ua)∂t+m2a∂(ua-uw)∂t,
where m1ka is the coefficient of air volume change with respect to a change in the net normal stress (q-ua) and m2a is the coefficient of air volume change with respect to a change in matric suction (ua-uw).Substituting (5) into (6), the governing equation for the air phase under constant loading dq/dt=0 can be written as
(7)∂ua∂t=-Ca∂uw∂t-Cva(∂2ua∂r2+1r∂ua∂r),
where Ca=m2a/(m1ka-m2a-uatmn0(1-Sr0)/(u¯0a)2) and Cvw=kaRT/gu¯0aM(m1ka-m2a-uatmn0(1-Sr0)/(u¯0a)2).Further, two transformed governing equations of water and air phases from (3) and (7) can be obtained:
(8)∂uw∂t=Wa(∂2ua∂r2+1r∂ua∂r)+Ww(∂2uw∂r2+1r∂uw∂r),∂ua∂t=Aa(∂2ua∂r2+1r∂ua∂r)+Aw(∂2uw∂r2+1r∂uw∂r),
where Wa=CwCva/(1-CaCw), Ww=-Cvw/(1-CaCw), Aa=-Cva/(1-CaCw) and Aw=-CaCvw/(1-CaCw).


### 2.3. Boundary and Initial Conditions

As shown in Figure1, the wall of sand drain is permeable and lateral boundary surface is impermeable or impeded; the boundary conditions for radial consolidation are
(9)uw(rw,t)=ua(rw,t)=0,(10)∂uw(re,t)∂r+awua(re,t)=∂ua(re,t)∂r+aaua(re,t)=0,
when aw=aa=0 or aw=kw0(re-rw)/kwr0, aa=ka0(re-rw)/kar0; (9) reflects that the lateral boundary surface is impermeable or impeded, respectively. kw0 and ka0 are the coefficients of permeability for water and air at lateral boundary, respectively. r0 is the thickness of lateral boundary.The initial condition can be written as(11)ua(r,0)=u0a,uw(r,0)=u0w,
where u0w and u0a are the initial pore-water and pore-air pressure distributions.


## 2.1. Mathematical Model for Axisymmetric Consolidation of Unsaturated Soils

As shown in Figure1, the unsaturated soil layer is described in a schematic diagram and thickness, H. Figures 1(a) and 1(b) are the sectional view and vertical view of the unsaturated soil layer, respectively. The radii of soil layer and sand drain are re and rw. A surcharge q is applied on the top surface of the soil layer. kw and ka are the permeability coefficients of water and air, respectively. The coordinated origin is selected at the center of the top surface and the z coordinate is positive downward. In this paper, the vertical prismatic blocks of soil surrounding the sand drains are simulated by cylindrical blocks, of radius S. Square pattern and triangular pattern are two common distributions of sand drains. For instance, if square pattern is adopted as shown in Figure 1(c), then S=0.564re. The basic assumptions are the same as those of Fredlund’s one-dimensional consolidation theory for unsaturated soils and the other assumptions are listed as follows.(1)
The seepage of water and air phases is considered only in radial direction.(2)
The top and bottom surfaces are impermeable and the right boundary surface is impermeable or impeded.(3)
The permeability coefficients of water and air phases are assumed to be constant.(4)
In the process of consolidation, the strain is small enough and the density of water is assumed to be constant.(a) The sectional view of sand drain foundation, (b) the vertical view of sand drain foundation, and (c) the distribution of sand drains (square pattern).
(a)
(b)
(c)


## 2.2. Governing Equations

The net flux of water through the soil layer is computed from the volume of water entering and leaving the soil layer within a period of time, with respect to Darcy’s law in the polar coordinate system:(1)∂(Vw/V0)∂t=kwγw(∂2uw∂r2+1r∂uw∂r),
where V0 is the initial total volume of the soil, ∂(Vw/V0)/∂t is the net flux of water per unit volume of the soil, uw is the excess pore-water pressure, γw is the water unit weight, and t is the time variable.The net flux of water per unit volume of the soil can be obtained by differentiating the water phase constitutive relation with respect to time:(2)∂(Vw/V0)∂t=m1kw∂(q-ua)∂t+m2w∂(ua-uw)∂t,
where m1kw is the coefficient of water volume change with respect to a change in the net normal stress (q-ua), m2w is the coefficient of water volume change with respect to a change in matric suction (ua-uw), and ua is the excess pore-air pressure.For constant loading,dq/dt=0, substituting (1) into (2), the governing equation for the water phase can be written as
(3)∂uw∂t=-Cw∂ua∂t-Cvw(∂2uw∂r2+1r∂uw∂r),
where Cw=(m1kw-m2w)/m2w and Cvw=kw/γwm2w.The net flux of air through the soil layer is computed from the volume of air entering and leaving the soil layer within a period of time, with respect to Fick’s law in the polar coordinate system:(4)∂(ρVa/V0)∂t=kag(∂2ua∂r2+1r∂ua∂r),
where ρ is the density of air phase, ∂(ρVa/V0)/∂t is the net mass rate of air flow per unit volume of the soil, and g is the gravitational acceleration.The total volume change of an unsaturated soil can be assumed to be small during the consolidation process. The volume of airVa can be related to the volume-mass properties of the soil: Va=(1-S0)n0V0, where n0 and Sr0 are the initial porosity and initial degree of saturation before loading. According to Boyle’s law, assuming there is no initial excess air pressure in the soil before loading, we have [17, 18]
(5)∂(Va/V0)∂t=kaRTgu-a0M(∂2ua∂r2+1r∂ua∂r)-uatmn0(1-S0)(u-a0)2∂ua∂t,
where R is the universal air constant (8.314J/molK), T=293.15K is the absolute temperature, M=29kg/kmol is the molecular mass of air, u¯0a is absolute pore-air pressure (i.e., u¯0a=u0a+uatm), and uatm=101.3kPa is atmospheric pressure.The derivative of the air phase constitutive relation with respect to time is equal to the net flux of air per unit volume of the soil:(6)∂(Va/V0)∂t=m1ka∂(q-ua)∂t+m2a∂(ua-uw)∂t,
where m1ka is the coefficient of air volume change with respect to a change in the net normal stress (q-ua) and m2a is the coefficient of air volume change with respect to a change in matric suction (ua-uw).Substituting (5) into (6), the governing equation for the air phase under constant loading dq/dt=0 can be written as
(7)∂ua∂t=-Ca∂uw∂t-Cva(∂2ua∂r2+1r∂ua∂r),
where Ca=m2a/(m1ka-m2a-uatmn0(1-Sr0)/(u¯0a)2) and Cvw=kaRT/gu¯0aM(m1ka-m2a-uatmn0(1-Sr0)/(u¯0a)2).Further, two transformed governing equations of water and air phases from (3) and (7) can be obtained:
(8)∂uw∂t=Wa(∂2ua∂r2+1r∂ua∂r)+Ww(∂2uw∂r2+1r∂uw∂r),∂ua∂t=Aa(∂2ua∂r2+1r∂ua∂r)+Aw(∂2uw∂r2+1r∂uw∂r),
where Wa=CwCva/(1-CaCw), Ww=-Cvw/(1-CaCw), Aa=-Cva/(1-CaCw) and Aw=-CaCvw/(1-CaCw).


## 2.3. Boundary and Initial Conditions

As shown in Figure1, the wall of sand drain is permeable and lateral boundary surface is impermeable or impeded; the boundary conditions for radial consolidation are
(9)uw(rw,t)=ua(rw,t)=0,(10)∂uw(re,t)∂r+awua(re,t)=∂ua(re,t)∂r+aaua(re,t)=0,
when aw=aa=0 or aw=kw0(re-rw)/kwr0, aa=ka0(re-rw)/kar0; (9) reflects that the lateral boundary surface is impermeable or impeded, respectively. kw0 and ka0 are the coefficients of permeability for water and air at lateral boundary, respectively. r0 is the thickness of lateral boundary.The initial condition can be written as(11)ua(r,0)=u0a,uw(r,0)=u0w,
where u0w and u0a are the initial pore-water and pore-air pressure distributions.


## 3. Differential Quadrature Formulation

DQM is a numerical solution technique for initial and/or boundary value problems, proposed by Bellman et al. [19] and Bellman and Casti [20]. According to DQM, a partial derivative of a function with respect to a variable can be approximated by a weighted linear sum of the function values at given discrete points. Chen et al. [13] employed DQM to solve one-dimensional consolidation problems in multilayered soils.To show the mathematic detail of DQM, consider a functionφ=φ(x) on the domain 0≤x≤a and the domain is dispersed as N points. Then the general differential quadrature approximation of the function at the ith discrete point is given by
(12)∂rφ∂xr|x=xi≅∑k=1NDik(r)φk,i=1,2,…,N,
where Dik(r) are the weighting coefficient of rth derivative, r<N.This paper adopts a method derived by Quan and Chang [21] which uses Lagrange polynomial to determine the weighting coefficients
(13)Dik(1)=∏v=1,v≠iN(xi-xv)(xi-xk)∏v=1,k≠iN(xk-xv),kkkkkkkkkki,k=1,2,…,N(k≠i),Dii(r)=-∑v=1,v≠iNDiv(r),i=1,2…,N,1≤r≤N-1,kkiDik(r)=r[Dii(r-1)Dik(1)-Dik(r-1)xi-xk],i,k=1,2…,N(k≠i),2≤r≤N-1.The soil layer is dispersed in the radial direction and the number of discrete points isN. Then, in order to solve the equations conveniently, the local coordinate ε is introduced. The relationship between local coordinate and integral coordinate is
(14)r=(0.5-ε)r1+(0.5+ε)rN,(-0.5≤ε≤0.5),
where r is the integral coordinate of soil layer and r1 and rN are the radius values of 1th and Nth point, respectively.The differential of (14) can be expressed as
(15)dr=(rN-r1)dε=(re-rw)dε.The relationship between partial differential of local coordinate and integral coordinate is shown as follows:(16)∂uw∂r=1re-rw∂uw∂ε,∂2uw∂r2=1(re-rw)2∂2uw∂ε2,∂ua∂r=1re-rw∂ua∂ε,∂2ua∂r2=1(re-rw)2∂2ua∂ε2.Hence, (16) can be approximated by DQM into
(17)∂uw∂r|r=rα=1re-rw∑β=1NDαβ(1)uβw,∂2uw∂r2|r=rα=1(re-rw)2∑β=1NDαβ(2)uβw,∂ua∂r|r=rα=1re-rw∑β=1NDαβ(1)uβa,∂2ua∂r2|r=rα=1(re-rw)2∑β=1NDαβ(2)uβa(α=1,2,…,N),
where Dαβ(1) and Dαβ(2) are the weighting coefficient matrices of the first order and second order derivatives, respectively.Letua¯=ua/q, uw¯=uw/q, and T=(-kw/γwm1ksh2)t, where m1ks=m1kw+m1ka; substituting (16) into the governing equations (8), one obtains
(18)∂uαw¯∂T=Wakw/-γwm1ks×[(∑β=2N-1Dαβ(2)uβa¯+Dα1(2)u1a¯+DαN(2)uNa¯)+1rα¯(∑β=2N-1Dαβ(1)uβa¯+Dα1(1)u1a¯+DαN(1)uNa¯)]+Wwkw/-γwm1ks×[(∑β=2N-1Dαβ(2)uβw¯+Dα1(2)u1w¯+DαN(2)uNw¯)+1rα¯(∑β=2N-1Dαβ(1)uβw¯+Dα1(1)u1w¯+DαN(1)uNw¯)],∂uαa¯∂T=Aakw/-γwm1ks×[(∑β=2N-1Dαβ(2)uβa¯+Dα1(2)u1a¯+DαN(2)uNa¯)+1rα¯(∑β=2N-1Dαβ(1)uβa¯+Dα1(1)u1a¯+DαN(1)uNa¯)]+Awkw/-γwm1ks×[(∑β=2N-1Dαβ(2)uβw¯+Dα1(2)u1w¯+DαN(2)uNw¯)+1rα¯(∑β=2N-1Dαβ(1)uβw¯+Dα1(1)u1w¯+DαN(1)uNw¯)],
where rα¯=(rα/(re-rw))(α=1,2,…,N), rα is the radius value of αth point, and uαw and uαa are the pore-water and pore-air pressures at the αth point.The boundary condition and initial condition are also approximated by DQM; (9)–(11) are transformed into
(19)u1w¯=u1a¯=0,(20)∂uNw¯∂r+awuNw¯=∂uNa¯∂r+aauNa¯=0,(21)uw¯|t=0=u0wq,(22)ua¯|t=0=u0aq.According to (19) and (20), (18) can be rewritten as
(23)∂Vw∂T=WaG(AVa+1r¯BVa)+WwG(AVw+1r¯BVw),∂Va∂T=AaG(AVa+1r¯BVa)+AwG(AVw+1r¯BVw),
where Vw=[u2w¯,u3w¯,…,uN-1w¯]′, Va=[u2a¯,u3a¯,…,uN-1a¯]′, A=A1+A2, and B=B1+B2;
(24)A1=[D22(2)⋯D2(N-1)(2)⋮⋱⋮D(N-1)2(2)⋯D(N-1)(N-1)(2)],B1=[D22(1)⋯D2(N-1)(1)⋮⋱⋮D(N-1)2(1)⋯D(N-1)(N-1)(1)],A2=M-1×Q×[D21(2)D2N(2)⋮⋮D(N-1)1(2)D(N-1)N(2)],B2=M-1×Q×[D21(1)D2N(1)⋮⋮D(N-1)1(1)D(N-1)N(1)],G=kw-γwm1ks,M=[10DN1(1)DNN(1)],Q=-[0⋯0DN2(1)⋯DN(N-1)(1)].Therefore, the governing equations of water and air are translated into two sets of ordinary differential equations (ODEs). The solutions of ODEs can be obtained by using Rong-Kutta method. Then we can apply the solutions into the following formulas: the average degree of consolidation, the radial displacement, and the vertical settlement.In unsaturated soils, the average degree of consolidation can be divided into two parts: the average degree of consolidation with respect to water phaseUw and the average degree of consolidation with respect to air phase Ua. To obtain the average degree of consolidation, two formulations are given by Fredlund and Rahardjo [7]. Consider
(25)Uw=1-∫0(re-rw)uwdr∫0(re-rw)u0wdr,(26)Ua=1-∫0(re-rw)uadr∫0(re-rw)u0adr.According to the two stress-state variable approaches [6, 7], volume strain is represented by the following constitutive equation for soil layer:
(27)∂εv∂t=(m1kw+m1ka)∂(q-ua)∂t+(m2w+m2a)∂(ua-uw)∂t.By integrating (26) with respect to time t from 0 to t, we get the expression of volumetric strain εv:
(28)εv=(m1kw+m1ka)[(q-ua)-(q-u0a)]+(m2w+m2a)[(ua-uw)-(u0a-u0w)].Volumetric strain consists of vertical strain and horizontal strain. For the axisymmetric consolidation, it is assumed in this paper that one-third of the volume strain is contributed from vertical strain and two-thirds are from horizontal strain. So the radial displacement and vertical settlement can be obtained, respectively. Consider(29)Sh=23(m1kw+m1ka)∫0re[(q-ua)-(q-u0a)]dr+23(m2w+m2a)∫0re[(ua-uw)-(u0a-u0w)]dr,Sv=13(m1kw+m1ka)∫0H[(q-ua)-(q-u0a)]dz+13(m2w+m2a)∫0H[(ua-uw)-(u0a-u0w)]dz.


## 4. An Example and Convergence Analysis

In this section, the convergence analysis, the average degree of consolidation, and the effects of displacements in radial and vertical directions are discussed by using a simple example.As shown in Figure1, a vertical loading q=100kPa is applied on the top surface of the soil layer. The lateral boundary is considered as impermeable. Considering a uniform initial excess pore-water and pore-air pressure distribution: u0a=5kPa, u0w=40kPa. Other parameters are n0=0.5, Sr0=0.8, kw=10-10m/s, m1kw=-0.15×10-4kPa-1, m2w=-0.61×10-4kPa-1, m1ka=0.06×10-4kPa-1, m2a=0.26×10-4kPa-1, rw=0.2m, re=1.8m, and uatm=101.3kPa. By applying the DQM solution presented above, the results of axisymmetric consolidation of unsaturated soils in a sand drain foundation were obtained.Tables1 and 2 list the solutions of excess pore-water and pore-air pressures at radius r=1m at different time factors and different equally spaced discrete points, respectively. The ratio of permeability coefficient is ka/kw=0.1. From Tables 1 and 2, the accuracy of solutions increases when the number of discrete point becomes big. It is obvious that 9 equally spaced grid points are sufficient to obtain the converged results of excess pore-water and pore-air pressures. Therefore, the number of discrete points N=9 is adopted in most case studies.Table 1
Different excess pore-water pressuresuw/q under time factor T with different discrete points N (ka/kw=0.1 and r=1 m).
T
N
=
5
N
=
7
N
=
9
N
=
11
N
=
13
N
=
15
N
=
17
1
×
10
-
6
0.400
0.400
0.400
0.400
0.400
0.400
0.400
1
×
10
-
4
0.400
0.400
0.400
0.400
0.400
0.400
0.400
1
×
10
-
1
0.395
0.395
0.395
0.395
0.395
0.395
0.395
2
×
10
-
1
0.386
0.386
0.385
0.385
0.385
0.385
0.385
5
×
10
-
1
0.359
0.354
0.350
0.350
0.349
0.349
0.349
1
×
1
0
0
0.326
0.316
0.310
0.310
0.309
0.309
0.308Table 2
Different excess pore-air pressuresua/q under time factor T with different discrete points N (ka/kw=0.1 and r=1 m).
T
N
=
5
N
=
7
N
=
9
N
=
11
N
=
13
N
=
15
N
=
17
1
×
10
-
6
0.050
0.050
0.050
0.050
0.050
0.050
0.050
1
×
10
-
4
0.050
0.050
0.050
0.050
0.050
0.050
0.050
5
×
10
-
2
0.048
0.048
0.048
0.048
0.048
0.048
0.048
1
×
10
-
1
0.046
0.045
0.045
0.045
0.045
0.045
0.045
2
×
10
-
1
0.042
0.041
0.040
0.040
0.040
0.040
0.010
5
×
10
-
1
0.035
0.033
0.032
0.032
0.032
0.031
0.031Figure2 shows the curves of average degree of consolidation with respect to water and air phases when N=9. The ratio of permeability coefficient is ka/kw=1. From Figure 2, it can be observed that the consolidation of air phase and water phase begins when T=1.0×10-3. When T=1.0×100 and T=1.0×101.5, the consolidation of air phase and water phase is almost finished, respectively. At the early stages, soil consolidation is mainly caused by the dissipation of excess pore-air pressure. But in the later stages, soil consolidation is mainly caused by the dissipation of excess pore-water pressure.Figure 2
The average degree of consolidation for the water and air phases with time factorT (ka/kw=1).Figure3 shows the radial displacement development with time and Figure 4 describes the settlement of top surface (H=10m) at different radii with time factor T. The number of discrete points is N=9. The ratio of permeability coefficient is ka/kw=0.1. From Figure 4, we can see that the soil settles earlier at the points with smaller radius, that is, closer to the sand drain. Figure 5 shows the settlement at different heights of soil layer at radius r=1m with time factor T. When the height of soil layer H is equal to 4m, the settlement is almost zero.Figure 3
The radial displacementSh with time factor T (ka/kw=0.1).Figure 4
Variations of soil layer settlementsSv at different radii with time factor T (ka/kw=0.1).Figure 5
Variations of soil layer settlementsSv at different soil layer H with time factor T (ka/kw=0.1 and r=1m).


## 5. Cases of Different Initial and Boundary Conditions

As the initial conditions do not need to be constant in the present DQM solution, it is easily used to analyze nonuniform initial pore-water and pore-air distribution problems. In this part, four different initial pore-water and pore-air distributions and some different boundary conditions are considered. A vertical loadingq=100kPa is applied on the top surface of soil layer. Other parameters are n0=0.5, Sr0=0.8, kw=10-10m/s, rw=0.2m, re=1.8m, m1kw=-0.15×10-4kPa-1, m2w=-0.61×10-4kPa-1, m1ka=0.06×10-4kPa-1, m2a=0.26×10-4kPa-1, and uatm=101.3kPa.


### 5.1. Effect of Different Initial Pore-Water/Air Pressure Distributions

In this section, in order to make the curves smooth, the number of discrete points is chosen to beN=17. The lateral boundary is considered as impermeable and the ratio of permeability coefficient with respect to air and water is ka/kw=0.1.Four different initial pore-water/air distributions are taken into account. In Figure6(a), the case of uniform initial pore-air and pore-water pressures assumes that u0a=5kPa and u0w=40kPa. Figure 6(b) presents the linear increasing initial pore-air and pore-water distributions. The lateral-skewed distributions initial pore-water/air pressures is shown in Figures 6(c) and 6(d) and they are represented by the following equations:
(30)uwq=(r/re)0.25b(2-r/re)b100.105b,uaq=(r/re)0.25b(2-r/re)b8×100.105b,
where b is the variable which can control the spread or skewness of the curves in Figures 6(c) and 6(d).Initial excess pore pressure distributions: (a) uniform, (b) linearly increasing, (c) lateral skewed (b=2), and (d) lateral skewed (b=6).
(a)
(b)
(c)
(d)Under the condition shown in Figure6(a), the dissipation curves of the excess pore-water and pore-air pressure at different distances from the sand drain are presented in Figures 7 and 8, respectively. The excess pore-air and pore-water pressures dissipate faster at the points with smaller radius, that is, closer to the sand drain. As the interval between the two curves becomes smaller, from r=0.6m to r=1.2m, it can be seen that the effect of the distance on the dissipation becomes less obvious.Figure 7
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(a) distribution (ka/kw=0.1).Figure 8
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(a) distribution (ka/kw=0.1).In Figure6(b), the initial pore-air and pore-water pressures are assumed not to be constants and to be linear, increasing when the radius increases. The initial pore-air pressure is assumed one-eighth of the initial pore-water pressure. In other words, from r=0.4m to r=1.6m, the initial pore-water pressure will increase from 30kPa to 110kPa and the initial pore-air pressure will increase from 3.75kPa to 13.75kPa. Figures 9 and 10 show the dissipation of excess pore-water and pore-air pressures at different distances from the sand drain under Figure 6(b) distributions.Figure 9
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(b) distribution (ka/kw=0.1).Figure 10
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(b) distribution (ka/kw=0.1).Comparing Figures7 and 9, it is found that the excess pore-water pressure dissipates completely at almost the same time. With the change of the initial pore-water pressure distribution, the time of dissipation with respect to water is not significantly affected. Comparing Figures 8 and 10, it can be seen that the excess pore-air pressure dissipates completely at almost the same time. The initial pore-air pressure distribution has little influence on the dissipation time of air phase.Under the conditions of Figures6(c) and 6(d), the dissipation curves of the excess pore-water and pore-air pressures for initial distributions where b=2,6 are shown in Figures 11, 12, 13, and 14. From these figures, whatever the initial pore-water and pore-air pressures are, the excess pore-water and pore-air pressures dissipate completely at almost the same time. The plateau period is longer at the points with larger radius.Figure 11
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(c) distribution (ka/kw=0.1).Figure 12
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(c) distribution (ka/kw=0.1).Figure 13
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(d) distribution (ka/kw=0.1).Figure 14
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(d) distribution (ka/kw=0.1).Here, excess pore pressure isochrones for the above four initial conditions are presented in Figures15 and 16. Figures 15(a) and 16(a) show the excess pore-water and pore-air pressure isochrones when the initial excess pore pressure is uniform at all radii. In Figures 15(b) and 16(b), the initial excess pore pressure distributions are linearly increasing. The pore pressure isochrones for lateral-skewed initial distributions are shown in Figures 15(c), 16(c), 15(d), and 16(d), respectively. From Figures 15(a) and 16(a), as the initial excess pore pressure is constant, the peak of curves in Figure 15(a) approaches to 0.4 and the peak of curves in Figure 13(a) approaches to 0.05. Comparing Figures 15(a) and 15(b), as the initial excess pore-water pressure is linearly increasing, the skewness is more evident during early stages of dissipation in Figure 15(b). But in the later stages, the influence of initial excess pore-water pressure becomes less significant. So are the isochrones of the excess pore-air pressure in Figures 16(a) and 16(b). Figures 15(c), 16(c), 15(d), and 16(d) show the excess pore-water and pore-air pressure isochrones for the lateral-skewed initial distribution when b=2,6. It is seen that some intersections appear in the isochrones during the early stages. It means that a redistribution of excess pore pressures occurs toward the region of minimal initial excess pore pressure during the early stages. With the time increasing, the isochrones of excess pore pressures become regular.Normalized pore-water pressure isochrones for different initial excess pore-water pressure distributions (ka/kw=0.1): (a) uniform, (b) linearly increasing, (c) lateral skewed with b=2, and (d) lateral skewed with b=6.
(a)
(b)
(c)
(d)Normalized pore-air pressure isochrones for different initial excess pore-air pressure distributions (ka/kw=0.1): (a) uniform, (b) linearly increasing, (c) lateral-skewed with b=2, and (d) lateral skewed with b=6.
(a)
(b)
(c)
(d)


### 5.2. Effect of Different Boundary Conditions

In this part, the effect of different lateral boundary conditions is studied. Uniform initial excess pore-water/air pressure distributionsu0a=5kPa and u0w=40kPa are adopted. Four different ratios ka/kw are considered, that is, 0.1, 1, 10, and 100 with kw=10-10m/s. Two different drainage conditions are analyzed: lateral surface is impervious and lateral surface is impeded. For the impeded surface, aw and aa in (9) are chosen to be 0.5.Figures17 and 18 shows the dissipation curves for the water and air pressures at radius r=1.2m under different boundary conditions. First, the excess pore-water and pore-air pressures dissipate faster in impeded surface than that in impermeable surface. Second, kw is the same in the cases of four ka/kw values. It is observed from the results that the bigger ka/kw is, the faster the excess pore-water and pore-air pressures dissipate. Third, the dissipation curves of the excess pore-air and pore-water pressures have almost the same shape under different values of ka/kw. Comparing Figure 17 with Figure 18, when the excess pore-air pressure dissipates almost completely, the dissipation of excess pore-water pressure enters into a plateau period. The bigger ka/kw is, the earlier the excess pore-air pressure dissipates completely. Therefore, it is shown that the value of ka affects the dissipation of pore-water pressure during the consolidation process.Figure 17
Variation of excess pore-water pressuresuw/q under different ka/kw with time factor T (r=1.2m).Figure 18
Variation of excess pore-air pressuresua/q under different ka/kw with time factor T (r=1.2m).


## 5.1. Effect of Different Initial Pore-Water/Air Pressure Distributions

In this section, in order to make the curves smooth, the number of discrete points is chosen to beN=17. The lateral boundary is considered as impermeable and the ratio of permeability coefficient with respect to air and water is ka/kw=0.1.Four different initial pore-water/air distributions are taken into account. In Figure6(a), the case of uniform initial pore-air and pore-water pressures assumes that u0a=5kPa and u0w=40kPa. Figure 6(b) presents the linear increasing initial pore-air and pore-water distributions. The lateral-skewed distributions initial pore-water/air pressures is shown in Figures 6(c) and 6(d) and they are represented by the following equations:
(30)uwq=(r/re)0.25b(2-r/re)b100.105b,uaq=(r/re)0.25b(2-r/re)b8×100.105b,
where b is the variable which can control the spread or skewness of the curves in Figures 6(c) and 6(d).Initial excess pore pressure distributions: (a) uniform, (b) linearly increasing, (c) lateral skewed (b=2), and (d) lateral skewed (b=6).
(a)
(b)
(c)
(d)Under the condition shown in Figure6(a), the dissipation curves of the excess pore-water and pore-air pressure at different distances from the sand drain are presented in Figures 7 and 8, respectively. The excess pore-air and pore-water pressures dissipate faster at the points with smaller radius, that is, closer to the sand drain. As the interval between the two curves becomes smaller, from r=0.6m to r=1.2m, it can be seen that the effect of the distance on the dissipation becomes less obvious.Figure 7
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(a) distribution (ka/kw=0.1).Figure 8
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(a) distribution (ka/kw=0.1).In Figure6(b), the initial pore-air and pore-water pressures are assumed not to be constants and to be linear, increasing when the radius increases. The initial pore-air pressure is assumed one-eighth of the initial pore-water pressure. In other words, from r=0.4m to r=1.6m, the initial pore-water pressure will increase from 30kPa to 110kPa and the initial pore-air pressure will increase from 3.75kPa to 13.75kPa. Figures 9 and 10 show the dissipation of excess pore-water and pore-air pressures at different distances from the sand drain under Figure 6(b) distributions.Figure 9
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(b) distribution (ka/kw=0.1).Figure 10
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(b) distribution (ka/kw=0.1).Comparing Figures7 and 9, it is found that the excess pore-water pressure dissipates completely at almost the same time. With the change of the initial pore-water pressure distribution, the time of dissipation with respect to water is not significantly affected. Comparing Figures 8 and 10, it can be seen that the excess pore-air pressure dissipates completely at almost the same time. The initial pore-air pressure distribution has little influence on the dissipation time of air phase.Under the conditions of Figures6(c) and 6(d), the dissipation curves of the excess pore-water and pore-air pressures for initial distributions where b=2,6 are shown in Figures 11, 12, 13, and 14. From these figures, whatever the initial pore-water and pore-air pressures are, the excess pore-water and pore-air pressures dissipate completely at almost the same time. The plateau period is longer at the points with larger radius.Figure 11
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(c) distribution (ka/kw=0.1).Figure 12
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(c) distribution (ka/kw=0.1).Figure 13
Variation of pore-water pressuresuw/q at different radii with time factor T under Figure 6(d) distribution (ka/kw=0.1).Figure 14
Variation of pore-air pressuresua/q at different radii with time factor T under Figure 6(d) distribution (ka/kw=0.1).Here, excess pore pressure isochrones for the above four initial conditions are presented in Figures15 and 16. Figures 15(a) and 16(a) show the excess pore-water and pore-air pressure isochrones when the initial excess pore pressure is uniform at all radii. In Figures 15(b) and 16(b), the initial excess pore pressure distributions are linearly increasing. The pore pressure isochrones for lateral-skewed initial distributions are shown in Figures 15(c), 16(c), 15(d), and 16(d), respectively. From Figures 15(a) and 16(a), as the initial excess pore pressure is constant, the peak of curves in Figure 15(a) approaches to 0.4 and the peak of curves in Figure 13(a) approaches to 0.05. Comparing Figures 15(a) and 15(b), as the initial excess pore-water pressure is linearly increasing, the skewness is more evident during early stages of dissipation in Figure 15(b). But in the later stages, the influence of initial excess pore-water pressure becomes less significant. So are the isochrones of the excess pore-air pressure in Figures 16(a) and 16(b). Figures 15(c), 16(c), 15(d), and 16(d) show the excess pore-water and pore-air pressure isochrones for the lateral-skewed initial distribution when b=2,6. It is seen that some intersections appear in the isochrones during the early stages. It means that a redistribution of excess pore pressures occurs toward the region of minimal initial excess pore pressure during the early stages. With the time increasing, the isochrones of excess pore pressures become regular.Normalized pore-water pressure isochrones for different initial excess pore-water pressure distributions (ka/kw=0.1): (a) uniform, (b) linearly increasing, (c) lateral skewed with b=2, and (d) lateral skewed with b=6.
(a)
(b)
(c)
(d)Normalized pore-air pressure isochrones for different initial excess pore-air pressure distributions (ka/kw=0.1): (a) uniform, (b) linearly increasing, (c) lateral-skewed with b=2, and (d) lateral skewed with b=6.
(a)
(b)
(c)
(d)


## 5.2. Effect of Different Boundary Conditions

In this part, the effect of different lateral boundary conditions is studied. Uniform initial excess pore-water/air pressure distributionsu0a=5kPa and u0w=40kPa are adopted. Four different ratios ka/kw are considered, that is, 0.1, 1, 10, and 100 with kw=10-10m/s. Two different drainage conditions are analyzed: lateral surface is impervious and lateral surface is impeded. For the impeded surface, aw and aa in (9) are chosen to be 0.5.Figures17 and 18 shows the dissipation curves for the water and air pressures at radius r=1.2m under different boundary conditions. First, the excess pore-water and pore-air pressures dissipate faster in impeded surface than that in impermeable surface. Second, kw is the same in the cases of four ka/kw values. It is observed from the results that the bigger ka/kw is, the faster the excess pore-water and pore-air pressures dissipate. Third, the dissipation curves of the excess pore-air and pore-water pressures have almost the same shape under different values of ka/kw. Comparing Figure 17 with Figure 18, when the excess pore-air pressure dissipates almost completely, the dissipation of excess pore-water pressure enters into a plateau period. The bigger ka/kw is, the earlier the excess pore-air pressure dissipates completely. Therefore, it is shown that the value of ka affects the dissipation of pore-water pressure during the consolidation process.Figure 17
Variation of excess pore-water pressuresuw/q under different ka/kw with time factor T (r=1.2m).Figure 18
Variation of excess pore-air pressuresua/q under different ka/kw with time factor T (r=1.2m).


## 6. Conclusion

This paper obtains a general solution to the axisymmetric consolidation of unsaturated soils by using differential quadrature method (DQM) based on Fredlund’s one-dimensional consolidation theory for unsaturated soil. With the use of DQM, the two governing equations are transformed into two sets of ordinary differential equations. The solutions to the transformed differential equations are obtained by Rong-Kutta method under different initial and boundary conditions. A case study has been presented for the analysis of unsaturated consolidation in a sand drain foundation. The convergence analysis, the average degree of consolidation for water and air phases, the settlement in radial and vertical directions, effects of different initial excess pore-air and pore-water pressure distributions, and effects of different boundary conditions have been presented. From the analyses, some conclusions can be drawn: (a) through compiling the programs, it is found that the DQ solution delivers accurate and stable results for unsaturated soil consolidation. Due to the uniform matrix structure of the DQ formulas, it is easy to compile computer programs when considering complicated conditions such as nonuniform pore-water/air distribution conditions; (b) soil consolidation is mainly caused by the dissipation of excess pore-air pressure in the early stages and by the dissipation of excess pore-water pressure in the later stages; (c) the bigger the ratio of the permeability coefficients for water and air is, the faster the excess pore-water and pore-air pressures dissipate; (d) The excess pore-air and pore-water pressures dissipate faster at the points closer to the sand drain; (e) the initial distribution has some effects on the early consolidation process of unsaturated soil and boundary condition has a significant effect on the whole consolidation process.



---

*Source: 497161-2013-12-11.xml*