Axisymmetric Consolidation of Unsaturated Soils by Differential Quadrature Method

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.

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 Figure 1, the unsaturated soil layer is described in a schematic diagram and thickness, 𝐻. 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 𝑟_𝑒 and 𝑟_𝑤. A surcharge 𝑞 is applied on the top surface of the soil layer. 𝑘_𝑤 and 𝑘_𝑎 are the permeability coefficients of water and air, respectively. The coordinated origin is selected at the center of the top surface and the 𝑧 coordinate is positive downward. In this paper, the vertical prismatic blocks of soil surrounding the sand drains are simulated by cylindrical blocks, of radius 𝑆. 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 𝑆=0.564𝑟_𝑒. 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)

(b)

(c)

Figure 1 

(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).

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𝑟𝜕𝑢_𝑤,𝜕𝑟 (1) where 𝑉₀ is the initial total volume of the soil, 𝜕(𝑉_𝑤/𝑉₀)/𝜕𝑡 is the net flux of water per unit volume of the soil, 𝑢_𝑤 is the excess pore-water pressure, 𝛾_𝑤 is the water unit weight, and 𝑡 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: 𝜕𝑉_𝑤/𝑉₀𝜕𝑡=𝑚^𝑤_1𝑘𝜕𝑞−𝑢_𝑎𝜕𝑡+𝑚^𝑤₂𝜕𝑢_𝑎−𝑢_𝑤,𝜕𝑡 (2) where 𝑚^𝑤_1𝑘 is the coefficient of water volume change with respect to a change in the net normal stress (𝑞−𝑢_𝑎), 𝑚^𝑤₂ is the coefficient of water volume change with respect to a change in matric suction (𝑢_𝑎−𝑢_𝑤), and 𝑢_𝑎 is the excess pore-air pressure.

For constant loading, 𝑑𝑞/𝑑𝑡=0, substituting (1) into (2), the governing equation for the water phase can be written as 𝜕𝑢^𝑤𝜕𝑡=−𝐶_𝑤𝜕𝑢^𝑎𝜕𝑡−𝐶^𝑤_𝑣𝜕²𝑢^𝑤𝜕𝑟²+1𝑟𝜕𝑢^𝑤,𝜕𝑟 (3) where 𝐶_𝑤=(𝑚^𝑤_1𝑘−𝑚^𝑤₂)/𝑚^𝑤₂ and 𝐶^𝑤_𝑣=𝑘_𝑤/𝛾_𝑤𝑚^𝑤₂.

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: 𝜕𝜌𝑉_𝑎/𝑉₀=𝑘𝜕𝑡𝑔^𝑎𝜕²𝑢_𝑎𝜕𝑟²+1𝑟𝜕𝑢_𝑎,𝜕𝑟 (4) where 𝜌 is the density of air phase, 𝜕(𝜌𝑉_𝑎/𝑉₀)/𝜕𝑡 is the net mass rate of air flow per unit volume of the soil, and 𝑔 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 air 𝑉_𝑎 can be related to the volume-mass properties of the soil: 𝑉_𝑎=(1−𝑆₀)𝑛₀𝑉₀, where 𝑛₀ and 𝑆_𝑟0 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] 𝜕𝑉_𝑎/𝑉₀=𝑘𝜕𝑡^𝑎𝑅𝑇𝑔𝑢⁰_𝑎𝑀𝜕²𝑢_𝑎𝜕𝑟²+1𝑟𝜕𝑢_𝑎−𝑢𝜕𝑟_atm𝑛₀1−𝑆₀𝑢⁰_𝑎²𝜕𝑢_𝑎,𝜕𝑡 (5) where 𝑅 is the universal air constant (8.314J/molK), 𝑇=293.15K is the absolute temperature, 𝑀=29kg/kmol is the molecular mass of air, 𝑢^𝑎₀ is absolute pore-air pressure (i.e., 𝑢^𝑎₀=𝑢^𝑎₀+𝑢^atm), and 𝑢^atm=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: 𝜕𝑉_𝑎/𝑉₀𝜕𝑡=𝑚^𝑎_1𝑘𝜕𝑞−𝑢_𝑎𝜕𝑡+𝑚^𝑎₂𝜕𝑢_𝑎−𝑢_𝑤,𝜕𝑡 (6) where 𝑚^𝑎_1𝑘 is the coefficient of air volume change with respect to a change in the net normal stress (𝑞−𝑢^𝑎) and 𝑚^𝑎₂ is the coefficient of air volume change with respect to a change in matric suction (𝑢^𝑎−𝑢^𝑤).

Substituting (5) into (6), the governing equation for the air phase under constant loading 𝑑𝑞/𝑑𝑡=0 can be written as 𝜕𝑢^𝑎𝜕𝑡=−𝐶_𝑎𝜕𝑢^𝑤𝜕𝑡−𝐶^𝑎_𝑣𝜕²𝑢^𝑎𝜕𝑟²+1𝑟𝜕𝑢^𝑎,𝜕𝑟 (7) where 𝐶_𝑎=𝑚^𝑎₂/(𝑚^𝑎_1𝑘−𝑚^𝑎₂−𝑢^atm𝑛₀(1−𝑆_𝑟0)/(𝑢^𝑎₀)²) and 𝐶^𝑤_𝑣=𝑘_𝑎𝑅𝑇/𝑔𝑢^𝑎₀𝑀(𝑚^𝑎_1𝑘−𝑚^𝑎₂−𝑢^atm𝑛₀(1−𝑆_𝑟0)/(𝑢^𝑎₀)²).

Further, two transformed governing equations of water and air phases from (3) and (7) can be obtained: 𝜕𝑢^𝑤𝜕𝑡=𝑊_𝑎𝜕²𝑢^𝑎𝜕𝑟²+1𝑟𝜕𝑢^𝑎𝜕𝑟+𝑊_𝑤𝜕²𝑢^𝑤𝜕𝑟²+1𝑟𝜕𝑢^𝑤,𝜕𝑟𝜕𝑢^𝑎𝜕𝑡=𝐴_𝑎𝜕²𝑢^𝑎𝜕𝑟²+1𝑟𝜕𝑢^𝑎𝜕𝑟+𝐴_𝑤𝜕²𝑢^𝑤𝜕𝑟²+1𝑟𝜕𝑢^𝑤,𝜕𝑟 (8) where 𝑊_𝑎=𝐶_𝑤𝐶^𝑎_𝑣/(1−𝐶_𝑎𝐶_𝑤), 𝑊_𝑤=−𝐶^𝑤_𝑣/(1−𝐶_𝑎𝐶_𝑤), 𝐴_𝑎=−𝐶^𝑎_𝑣/(1−𝐶_𝑎𝐶_𝑤) and 𝐴_𝑤=−𝐶_𝑎𝐶^𝑤_𝑣/(1−𝐶_𝑎𝐶_𝑤).

2.3. Boundary and Initial Conditions

As shown in Figure 1, the wall of sand drain is permeable and lateral boundary surface is impermeable or impeded; the boundary conditions for radial consolidation are 𝑢^𝑤𝑟_𝑤,𝑡=𝑢^𝑎𝑟_𝑤,𝑡=0,𝜕𝑢^𝑤𝑟_𝑒,𝑡𝜕𝑟+𝑎_𝑤𝑢^𝑎𝑟_𝑒,𝑡=𝜕𝑢^𝑎𝑟_𝑒,𝑡𝜕𝑟+𝑎_𝑎𝑢^𝑎𝑟_𝑒,𝑡=0, (9)(10) when 𝑎_𝑤=𝑎_𝑎=0 or 𝑎_𝑤=𝑘⁰_𝑤(𝑟_𝑒−𝑟_𝑤)/𝑘_𝑤𝑟₀, 𝑎_𝑎=𝑘⁰_𝑎(𝑟_𝑒−𝑟_𝑤)/𝑘_𝑎𝑟₀; (9) reflects that the lateral boundary surface is impermeable or impeded, respectively. 𝑘⁰_𝑤 and 𝑘⁰_𝑎 are the coefficients of permeability for water and air at lateral boundary, respectively. 𝑟₀ is the thickness of lateral boundary.

The initial condition can be written as 𝑢^𝑎(𝑟,0)=𝑢^𝑎₀,𝑢^𝑤(𝑟,0)=𝑢^𝑤₀, (11) where 𝑢^𝑤₀ and 𝑢^𝑎₀ 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 𝜑=𝜑(𝑥) on the domain 0≤𝑥≤𝑎 and the domain is dispersed as 𝑁 points. Then the general differential quadrature approximation of the function at the 𝑖th discrete point is given by 𝜕^𝑟𝜑𝜕𝑥^𝑟_{𝑥=𝑥𝑖}≅^𝑁_𝑘=1𝐷^(𝑟)_𝑖𝑘𝜑_𝑘,𝑖=1,2,…,𝑁, (12) where 𝐷^(𝑟)_𝑖𝑘 are the weighting coefficient of 𝑟th derivative, 𝑟<𝑁.

This paper adopts a method derived by Quan and Chang [21] which uses Lagrange polynomial to determine the weighting coefficients 𝐷⁽¹⁾_𝑖𝑘=∏^𝑁_{𝑣=1,𝑣≠𝑖}𝑥_𝑖−𝑥_𝑣𝑥_𝑖−𝑥_𝑘∏^𝑁_{𝑣=1,𝑘≠𝑖}𝑥_𝑘−𝑥_𝑣,(),𝐷𝑖,𝑘=1,2,…,𝑁𝑘≠𝑖_{(𝑟)𝑖𝑖}=−^𝑁_{𝑣=1,𝑣≠𝑖}𝐷^(𝑟)_𝑖𝑣𝐷,𝑖=1,2…,𝑁,1≤𝑟≤𝑁−1,_{(𝑟)𝑖𝑘}=𝑟𝐷^(𝑟−1)_𝑖𝑖𝐷⁽¹⁾_𝑖𝑘−𝐷^(𝑟−1)_𝑖𝑘𝑥_𝑖−𝑥_𝑘,𝑖,𝑘=1,2…,𝑁(𝑘≠𝑖),2≤𝑟≤𝑁−1. (13)

The soil layer is dispersed in the radial direction and the number of discrete points is 𝑁. Then, in order to solve the equations conveniently, the local coordinate 𝜀 is introduced. The relationship between local coordinate and integral coordinate is 𝑟=(0.5−𝜀)𝑟₁+(0.5+𝜀)𝑟_𝑁,(−0.5≤𝜀≤0.5), (14) where 𝑟 is the integral coordinate of soil layer and 𝑟₁ and 𝑟_𝑁 are the radius values of 1th and 𝑁th point, respectively.

The differential of (14) can be expressed as 𝑑𝑟=𝑟_𝑁−𝑟₁𝑑𝜀=𝑟_𝑒−𝑟_𝑤𝑑𝜀. (15)

The relationship between partial differential of local coordinate and integral coordinate is shown as follows: 𝜕𝑢^𝑤=1𝜕𝑟𝑟_𝑒−𝑟_𝑤𝜕𝑢^𝑤,𝜕𝜕𝜀₂𝑢^𝑤𝜕𝑟²=1𝑟_𝑒−𝑟_𝑤²𝜕²𝑢^𝑤𝜕𝜀²,𝜕𝑢^𝑎=1𝜕𝑟𝑟_𝑒−𝑟_𝑤𝜕𝑢^𝑎,𝜕𝜕𝜀₂𝑢^𝑎𝜕𝑟²=1𝑟_𝑒−𝑟_𝑤²𝜕²𝑢^𝑎𝜕𝜀². (16)

Hence, (16) can be approximated by DQM into 𝜕𝑢^𝑤𝜕𝑟_{𝑟=𝑟𝛼}=1𝑟_𝑒−𝑟_𝑤^𝑁_𝛽=1𝐷⁽¹⁾_𝛼𝛽𝑢^𝑤_𝛽,𝜕²𝑢^𝑤𝜕𝑟²_{𝑟=𝑟𝛼}=1𝑟_𝑒−𝑟_𝑤^2𝑁_𝛽=1𝐷⁽²⁾_𝛼𝛽𝑢^𝑤_𝛽,𝜕𝑢^𝑎𝜕𝑟_{𝑟=𝑟𝛼}=1𝑟_𝑒−𝑟_𝑤^𝑁_𝛽=1𝐷⁽¹⁾_𝛼𝛽𝑢^𝑎_𝛽,𝜕²𝑢^𝑎𝜕𝑟²_{𝑟=𝑟𝛼}=1𝑟_𝑒−𝑟_𝑤^2𝑁_𝛽=1𝐷⁽²⁾_𝛼𝛽𝑢^𝑎_𝛽(𝛼=1,2,…,𝑁), (17) where 𝐷⁽¹⁾_𝛼𝛽 and 𝐷⁽²⁾_𝛼𝛽 are the weighting coefficient matrices of the first order and second order derivatives, respectively.

Let 𝑢^𝑎=𝑢^𝑎/𝑞, 𝑢^𝑤=𝑢^𝑤/𝑞, and 𝑇=(−𝑘_𝑤/𝛾_𝑤𝑚^𝑠_1𝑘ℎ²)𝑡, where 𝑚^𝑠_1𝑘=𝑚^𝑤_1𝑘+𝑚^𝑎_1𝑘; substituting (16) into the governing equations (8), one obtains 𝜕𝑢^𝑤_𝛼=𝑊𝜕𝑇𝑘_𝑤/−𝛾^𝑎_𝑤𝑚^𝑠_1𝑘×^𝑁−1_𝛽=2𝐷⁽²⁾_𝛼𝛽𝑢^𝑎_𝛽+𝐷⁽²⁾_𝛼1𝑢^𝑎₁+𝐷⁽²⁾_𝛼𝑁𝑢^𝑎_𝑁+1𝑟_𝛼^𝑁−1_𝛽=2𝐷⁽¹⁾_𝛼𝛽𝑢^𝑎_𝛽+𝐷⁽¹⁾_𝛼1𝑢^𝑎₁+𝐷⁽¹⁾_𝛼𝑁𝑢^𝑎_𝑁+𝑊_𝑤𝑘_𝑤/−𝛾_𝑤𝑚^𝑠_1𝑘×^𝑁−1_𝛽=2𝐷⁽²⁾_𝛼𝛽𝑢^𝑤_𝛽+𝐷⁽²⁾_𝛼1𝑢^𝑤₁+𝐷⁽²⁾_𝛼𝑁𝑢^𝑤_𝑁+1𝑟_𝛼^𝑁−1_𝛽=2𝐷⁽¹⁾_𝛼𝛽𝑢^𝑤_𝛽+𝐷⁽¹⁾_𝛼1𝑢^𝑤₁+𝐷⁽¹⁾_𝛼𝑁𝑢^𝑤_𝑁,𝜕𝑢^𝑎_𝛼=𝐴𝜕𝑇𝑘_𝑤/−𝛾^𝑎_𝑤𝑚^𝑠_1𝑘×^𝑁−1_𝛽=2𝐷⁽²⁾_𝛼𝛽𝑢^𝑎_𝛽+𝐷⁽²⁾_𝛼1𝑢^𝑎₁+𝐷⁽²⁾_𝛼𝑁𝑢^𝑎_𝑁+1𝑟_𝛼^𝑁−1_𝛽=2𝐷⁽¹⁾_𝛼𝛽𝑢^𝑎_𝛽+𝐷⁽¹⁾_𝛼1𝑢^𝑎₁+𝐷⁽¹⁾_𝛼𝑁𝑢^𝑎_𝑁+𝐴_𝑤𝑘_𝑤/−𝛾_𝑤𝑚^𝑠_1𝑘×^𝑁−1_𝛽=2𝐷⁽²⁾_𝛼𝛽𝑢^𝑤_𝛽+𝐷⁽²⁾_𝛼1𝑢^𝑤₁+𝐷⁽²⁾_𝛼𝑁𝑢^𝑤_𝑁+1𝑟_𝛼^𝑁−1_𝛽=2𝐷⁽¹⁾_𝛼𝛽𝑢^𝑤_𝛽+𝐷⁽¹⁾_𝛼1𝑢^𝑤₁+𝐷⁽¹⁾_𝛼𝑁𝑢^𝑤_𝑁, (18) where 𝑟_𝛼=(𝑟_𝛼/(𝑟_𝑒−𝑟_𝑤))(𝛼=1,2,…,𝑁), 𝑟_𝛼 is the radius value of 𝛼th point, and 𝑢^𝑤_𝛼 and 𝑢^𝑎_𝛼 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 𝑢^𝑤₁=𝑢^𝑎₁𝜕=0,𝑢^𝑤_𝑁𝜕𝑟+𝑎_𝑤𝑢^𝑤_𝑁=𝜕𝑢^𝑎_𝑁𝜕𝑟+𝑎_𝑎𝑢^𝑎_𝑁=0,𝑢^𝑤_𝑡=0=𝑢^𝑤₀𝑞,𝑢^𝑎_𝑡=0=𝑢^𝑎₀𝑞. (19)(20)(21)(22)

According to (19) and (20), (18) can be rewritten as 𝜕𝐕^𝑤=𝑊𝜕𝑇𝐺^𝑎𝐀𝐕^𝑎+1𝑟𝐁𝐕^𝑎+𝑊_𝑤𝐺𝐀𝐕^𝑤+1𝑟𝐁𝐕^𝑤,𝜕𝐕^𝑎=𝐴𝜕𝑇𝐺^𝑎𝐀𝐕^𝑎+1𝑟𝐁𝐕^𝑎+𝐴_𝑤𝐺𝐀𝐕^𝑤+1𝑟𝐁𝐕^𝑤, (23) where 𝐕^𝑤=[𝑢^𝑤₂,𝑢^𝑤₃,…,𝑢^𝑤_𝑁−1]^′, 𝐕^𝑎=[𝑢^𝑎₂,𝑢^𝑎₃,…,𝑢^𝑎_𝑁−1]^′, 𝐀=𝐀_𝟏+𝐀_𝟐, and 𝐁=𝐁_𝟏+𝐁_𝟐; 𝐀_𝟏=𝐷⁽²⁾₂₂⋯𝐷⁽²⁾_2(𝑁−1)𝐷⋮⋱⋮_{(2)(𝑁−1)2}⋯𝐷⁽²⁾_{(𝑁−1)(𝑁−1)},𝐁_𝟏=𝐷⁽¹⁾₂₂⋯𝐷⁽¹⁾_2(𝑁−1)𝐷⋮⋱⋮_{(1)(𝑁−1)2}⋯𝐷⁽¹⁾_{(𝑁−1)(𝑁−1)},𝐀_𝟐=𝐌⁻¹𝐷×𝐐×⁽²⁾²¹𝐷⁽²⁾_2𝑁𝐷⋮⋮_{(2)(𝑁−1)1}𝐷⁽²⁾_{(𝑁−1)𝑁},𝐁_𝟐=𝐌⁻¹𝐷×𝐐×⁽¹⁾²¹𝐷⁽¹⁾_2𝑁𝐷⋮⋮_{(1)(𝑁−1)1}𝐷⁽¹⁾_{(𝑁−1)𝑁},𝑘𝐺=^𝑤−𝛾_𝑤𝑚^𝑠_1𝑘,𝐷𝐌=10_(1)𝑁1𝐷⁽¹⁾_𝑁𝑁,𝐷𝐐=−0⋯0_(1)𝑁2⋯𝐷_𝑁()⁽¹⁾_𝑁−1. (24)

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 phase 𝑈^𝑤 and the average degree of consolidation with respect to air phase 𝑈^𝑎. To obtain the average degree of consolidation, two formulations are given by Fredlund and Rahardjo [7]. Consider 𝑈^𝑤∫=1−^{(𝑟𝑒−𝑟𝑤)0}𝑢^𝑤𝑑𝑟∫^{(𝑟𝑒−𝑟𝑤)}₀𝑢^𝑤₀𝑈𝑑𝑟_𝑎,∫=1−^{(𝑟𝑒−𝑟𝑤)0}𝑢^𝑎𝑑𝑟∫^{(𝑟𝑒−𝑟𝑤)}₀𝑢^𝑎₀𝑑𝑟. (25)(26)

According to the two stress-state variable approaches [6, 7], volume strain is represented by the following constitutive equation for soil layer: 𝜕𝜀_𝑣=𝑚𝜕𝑡^𝑤1𝑘+𝑚^𝑎_1𝑘𝜕𝑞−𝑢^𝑎+𝑚𝜕𝑡^𝑤2+𝑚^𝑎₂𝜕𝑢^𝑎−𝑢^𝑤.𝜕𝑡 (27)

By integrating (26) with respect to time 𝑡 from 0 to 𝑡, we get the expression of volumetric strain 𝜀_𝑣: 𝜀_𝑣=𝑚^𝑤_1𝑘+𝑚^𝑎_1𝑘𝑞−𝑢^𝑎−𝑞−𝑢^𝑎₀+𝑚^𝑤₂+𝑚^𝑎₂𝑢^𝑎−𝑢^𝑤−𝑢^𝑎₀−𝑢^𝑤₀. (28)

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 𝑆_ℎ=23𝑚^𝑤_1𝑘+𝑚^𝑎_1𝑘^𝑟𝑒₀𝑞−𝑢^𝑎−𝑞−𝑢^𝑎₀+2𝑑𝑟3𝑚^𝑤₂+𝑚^𝑎₂^𝑟𝑒₀𝑢^𝑎−𝑢^𝑤−𝑢^𝑎₀−𝑢^𝑤₀𝑆𝑑𝑟,_𝑣=13𝑚^𝑤_1𝑘+𝑚^𝑎_1𝑘^𝐻₀𝑞−𝑢^𝑎−𝑞−𝑢^𝑎₀+1𝑑𝑧3𝑚^𝑤₂+𝑚^𝑎₂^𝐻₀𝑢^𝑎−𝑢^𝑤−𝑢^𝑎₀−𝑢^𝑤₀𝑑𝑧. (29)

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 Figure 1, a vertical loading 𝑞=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: 𝑢^𝑎₀=5kPa, 𝑢^𝑤₀=40kPa. Other parameters are 𝑛₀=0.5, 𝑆_𝑟0=0.8, 𝑘_𝑤=10⁻¹⁰m/s, 𝑚^𝑤_1𝑘=−0.15×10⁻⁴kPa⁻¹, 𝑚^𝑤₂=−0.61×10⁻⁴kPa⁻¹, 𝑚^𝑎_1𝑘=0.06×10⁻⁴kPa⁻¹, 𝑚^𝑎₂=0.26×10⁻⁴kPa⁻¹, 𝑟_𝑤=0.2m, 𝑟_𝑒=1.8m, and 𝑢_atm=101.3kPa. By applying the DQM solution presented above, the results of axisymmetric consolidation of unsaturated soils in a sand drain foundation were obtained.

Tables 1 and 2 list the solutions of excess pore-water and pore-air pressures at radius 𝑟=1m at different time factors and different equally spaced discrete points, respectively. The ratio of permeability coefficient is 𝑘_𝑎/𝑘_𝑤=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 𝑁=9 is adopted in most case studies.

Table 1 

Different excess pore-water pressures 𝑢_𝑤/𝑞 under time factor 𝑇 with different discrete points 𝑁 (𝑘_𝑎/𝑘_𝑤=0.1 and 𝑟=1 m).

Table 2 

Different excess pore-air pressures 𝑢_𝑎/𝑞 under time factor 𝑇 with different discrete points 𝑁 (𝑘_𝑎/𝑘_𝑤=0.1 and 𝑟=1 m).

Figure 2 shows the curves of average degree of consolidation with respect to water and air phases when 𝑁=9. The ratio of permeability coefficient is 𝑘_𝑎/𝑘_𝑤=1. From Figure 2, it can be observed that the consolidation of air phase and water phase begins when 𝑇=1.0×10⁻³. When 𝑇=1.0×10⁰ and 𝑇=1.0×10^1.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 factor 𝑇 (𝑘𝑎/𝑘𝑤=1).

Figure 3 shows the radial displacement development with time and Figure 4 describes the settlement of top surface (𝐻=10m) at different radii with time factor 𝑇. The number of discrete points is 𝑁=9. The ratio of permeability coefficient is 𝑘_𝑎/𝑘_𝑤=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 𝑟=1m with time factor 𝑇. When the height of soil layer 𝐻 is equal to 4m, the settlement is almost zero.

Figure 3 

The radial displacement 𝑆ℎ with time factor 𝑇 (𝑘𝑎/𝑘𝑤=0.1).

Figure 4 

Variations of soil layer settlements 𝑆𝑣 at different radii with time factor 𝑇 (𝑘𝑎/𝑘𝑤=0.1).

Figure 5 

Variations of soil layer settlements 𝑆𝑣 at different soil layer 𝐻 with time factor 𝑇 (𝑘𝑎/𝑘𝑤=0.1 and 𝑟=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 loading 𝑞=100kPa is applied on the top surface of soil layer. Other parameters are 𝑛₀=0.5, 𝑆_𝑟0=0.8, 𝑘_𝑤=10⁻¹⁰m/s, 𝑟_𝑤=0.2m, 𝑟_𝑒=1.8m, 𝑚^𝑤_1𝑘=−0.15×10⁻⁴kPa⁻¹, 𝑚^𝑤₂=−0.61×10⁻⁴kPa⁻¹, 𝑚^𝑎_1𝑘=0.06×10⁻⁴kPa⁻¹, 𝑚^𝑎₂=0.26×10⁻⁴kPa⁻¹, and 𝑢_atm=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 be 𝑁=17. The lateral boundary is considered as impermeable and the ratio of permeability coefficient with respect to air and water is 𝑘_𝑎/𝑘_𝑤=0.1.

Four different initial pore-water/air distributions are taken into account. In Figure 6(a), the case of uniform initial pore-air and pore-water pressures assumes that 𝑢^𝑎₀=5kPa and 𝑢^𝑤₀=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: 𝑢^𝑤𝑞=𝑟/𝑟_𝑒^0.25𝑏2−𝑟/𝑟_𝑒^𝑏10^0.105𝑏,𝑢^𝑎𝑞=𝑟/𝑟_𝑒^0.25𝑏2−𝑟/𝑟_𝑒^𝑏8×10^0.105𝑏, (30) where 𝑏 is the variable which can control the spread or skewness of the curves in Figures 6(c) and 6(d).

(a)

(b)

(c)

(d)

Figure 6 

Initial excess pore pressure distributions: (a) uniform, (b) linearly increasing, (c) lateral skewed (𝑏=2), and (d) lateral skewed (𝑏=6).

Under the condition shown in Figure 6(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 𝑟=0.6m to 𝑟=1.2m, it can be seen that the effect of the distance on the dissipation becomes less obvious.

Figure 7 

Variation of pore-water pressures 𝑢𝑤/𝑞 at different radii with time factor 𝑇 under Figure 6(a) distribution (𝑘𝑎/𝑘𝑤=0.1).

Figure 8 

Variation of pore-air pressures 𝑢𝑎/𝑞 at different radii with time factor 𝑇 under Figure 6(a) distribution (𝑘𝑎/𝑘𝑤=0.1).

In Figure 6(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 𝑟=0.4m to 𝑟=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 pressures 𝑢𝑤/𝑞 at different radii with time factor 𝑇 under Figure 6(b) distribution (𝑘𝑎/𝑘𝑤=0.1).

Figure 10 

Variation of pore-air pressures 𝑢𝑎/𝑞 at different radii with time factor 𝑇 under Figure 6(b) distribution (𝑘𝑎/𝑘𝑤=0.1).

Comparing Figures 7 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 Figures 6(c) and 6(d), the dissipation curves of the excess pore-water and pore-air pressures for initial distributions where 𝑏=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 pressures 𝑢𝑤/𝑞 at different radii with time factor 𝑇 under Figure 6(c) distribution (𝑘𝑎/𝑘𝑤=0.1).

Figure 12 

Variation of pore-air pressures 𝑢𝑎/𝑞 at different radii with time factor 𝑇 under Figure 6(c) distribution (𝑘𝑎/𝑘𝑤=0.1).

Figure 13 

Variation of pore-water pressures 𝑢𝑤/𝑞 at different radii with time factor 𝑇 under Figure 6(d) distribution (𝑘𝑎/𝑘𝑤=0.1).

Figure 14 

Variation of pore-air pressures 𝑢𝑎/𝑞 at different radii with time factor 𝑇 under Figure 6(d) distribution (𝑘𝑎/𝑘𝑤=0.1).

Here, excess pore pressure isochrones for the above four initial conditions are presented in Figures 15 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 𝑏=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.

(a)

(b)

(c)

(d)

Figure 15 

Normalized pore-water pressure isochrones for different initial excess pore-water pressure distributions (𝑘𝑎/𝑘𝑤=0.1): (a) uniform, (b) linearly increasing, (c) lateral skewed with 𝑏=2, and (d) lateral skewed with 𝑏=6.

(a)

(b)

(c)

(d)

Figure 16 

Normalized pore-air pressure isochrones for different initial excess pore-air pressure distributions (𝑘𝑎/𝑘𝑤=0.1): (a) uniform, (b) linearly increasing, (c) lateral-skewed with 𝑏=2, and (d) lateral skewed with 𝑏=6.

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 distributions 𝑢^𝑎₀=5kPa and 𝑢^𝑤₀=40kPa are adopted. Four different ratios 𝑘_𝑎/𝑘_𝑤 are considered, that is, 0.1, 1, 10, and 100 with 𝑘_𝑤=10⁻¹⁰m/s. Two different drainage conditions are analyzed: lateral surface is impervious and lateral surface is impeded. For the impeded surface, 𝑎_𝑤 and 𝑎_𝑎 in (9) are chosen to be 0.5.

Figures 17 and 18 shows the dissipation curves for the water and air pressures at radius 𝑟=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, 𝑘_𝑤 is the same in the cases of four 𝑘_𝑎/𝑘_𝑤 values. It is observed from the results that the bigger 𝑘_𝑎/𝑘_𝑤 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 𝑘_𝑎/𝑘_𝑤. 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 𝑘_𝑎/𝑘_𝑤 is, the earlier the excess pore-air pressure dissipates completely. Therefore, it is shown that the value of 𝑘_𝑎 affects the dissipation of pore-water pressure during the consolidation process.

Figure 17 

Variation of excess pore-water pressures 𝑢𝑤/𝑞 under different 𝑘𝑎/𝑘𝑤 with time factor 𝑇 (𝑟=1.2m).

Figure 18 

Variation of excess pore-air pressures 𝑢𝑎/𝑞 under different 𝑘𝑎/𝑘𝑤 with time factor 𝑇 (𝑟=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.

Notation

𝐶𝑤𝑣:	The consolidation coefficient of water
𝐶𝑎𝑣:	The consolidation coefficient of air
𝐷(𝑟)𝑖𝑘:	The weighting coefficient
𝜀𝑣:	Volumetric strain
𝑔:	The gravitational acceleration
𝐻:	The thickness of soil layer
𝑘𝑎:	The permeability coefficients of air
𝑘0𝑎:	The coefficients of permeability for water at lateral boundary
𝑘𝑤:	The permeability coefficients of water
𝑘0𝑤:	The coefficients of permeability for air at lateral boundary
𝑚𝑎1𝑘:	The coefficient of air volume change with respect to a change in the net normal stress
𝑚𝑤1𝑘:	The coefficient of water volume change with respect to a change in the net normal stress
𝑚𝑎2:	The coefficient of air volume change with respect to a change in matric suction
𝑚𝑤2:	The coefficient of water volume change with respect to a change in matric suction
𝑛0:	The initial porosity
𝑁:	The discrete point
𝑃𝑣:	A vacuum pressure
𝜌:	The density of air phase
𝑞:	The vertical loading
𝑟0:	The thickness of lateral boundary
𝑟𝑒:	The radius of soil layer
𝑟𝑤:	The radius of sand drain
𝑟𝛼:	The radius value of 𝛼th point
𝛾𝑤:	The unit weight of water phase
𝑅:	The universal air constant
𝑆:	The distance between two sand drains
𝑆𝑟0:	The initial degree of saturation
𝑆𝑣:	Settlement
𝑆ℎ:	The radial displacement
𝑡:	Time
𝑇:	The time factor
𝑢𝑎:	The excess pore-air pressure
𝑢atm:	Atmospheric pressure
𝑢𝑎0:	The initial excess pore-air pressure
𝑢𝑤:	The excess pore-water pressure
𝑢𝑤0:	The initial excess pore-water pressure
𝑈𝑎:	The average degree of consolidation with respect to air
𝑈𝑤:	The average degree of consolidation with respect to water
𝑉𝑎:	The air volume of soil layer
𝑉𝑤:	The water volume of soil layer
𝑉0:	Initial total volume of the soil.

Acknowledgments

The author gratefully acknowledges the financial support from the Macau Science and Technology Development Fund (Grant no. FDCT/011/2013/A1) and the University of Macau Research Fund (Grants nos. MYRG189(Y2-L3)-FST11-ZWH and MYRG067(Y2-L2)-FST12-ZWH).