#### 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 . 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)**

##### 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: 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: where 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, , substituting (1) into (2), the governing equation for the water phase can be written as where 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: 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: , where and 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] where is the universal air constant , is the absolute temperature, is the molecular mass of air, is absolute pore-air pressure (i.e., ), and 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: where 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 can be written as where and .

Further, two transformed governing equations of water and air phases from (3) and (7) can be obtained: where , , and .

##### 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 when 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 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 and the domain is dispersed as points. Then the general differential quadrature approximation of the function at the discrete point is given by where are the weighting coefficient of derivative, .

This paper adopts a method derived by Quan and Chang [21] which uses Lagrange polynomial to determine the weighting coefficients

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 where is the integral coordinate of soil layer and and are the radius values of and point, respectively.

The differential of (14) can be expressed as

The relationship between partial differential of local coordinate and integral coordinate is shown as follows:

Hence, (16) can be approximated by DQM into where and are the weighting coefficient matrices of the first order and second order derivatives, respectively.

Let , , and , where ; substituting (16) into the governing equations (8), one obtains where , is the radius value of point, and and are the pore-water and pore-air pressures at the point.

The boundary condition and initial condition are also approximated by DQM; (9)–(11) are transformed into

According to (19) and (20), (18) can be rewritten as where , , , and ;

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

According to the two stress-state variable approaches [6, 7], volume strain is represented by the following constitutive equation for soil layer:

By integrating (26) with respect to time from to , we get the expression of volumetric strain :

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

#### 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 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: , . Other parameters are , , , , , , , , , and . 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 at different time factors and different equally spaced discrete points, respectively. The ratio of permeability coefficient is . 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 is adopted in most case studies.

Figure 2 shows the curves of average degree of consolidation with respect to water and air phases when . The ratio of permeability coefficient is . From Figure 2, it can be observed that the consolidation of air phase and water phase begins when . When and , 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 3 shows the radial displacement development with time and Figure 4 describes the settlement of top surface () at different radii with time factor . The number of discrete points is . The ratio of permeability coefficient is . 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 with time factor . When the height of soil layer is equal to , the settlement is almost zero.

#### 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 is applied on the top surface of soil layer. Other parameters are , , , , , , , , , and .

##### 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 . The lateral boundary is considered as impermeable and the ratio of permeability coefficient with respect to air and water is .

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 and . 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: 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)**

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 to , it can be seen that the effect of the distance on the dissipation becomes less obvious.

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 to , the initial pore-water pressure will increase from to and the initial pore-air pressure will increase from to . 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.

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

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 and the peak of curves in Figure 13(a) approaches to . 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 . 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)**

**(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 distributions and are adopted. Four different ratios are considered, that is, 0.1, 1, 10, and 100 with . 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 .

Figures 17 and 18 shows the dissipation curves for the water and air pressures at radius 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.

#### 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 |

: | The coefficients of permeability for water at lateral boundary |

: | The permeability coefficients of water |

: | The coefficients of permeability for air at lateral boundary |

: | The coefficient of air volume change with respect to a change in the net normal stress |

: | The coefficient of water volume change with respect to a change in the net normal stress |

: | The coefficient of air volume change with respect to a change in matric suction |

: | The coefficient of water volume change with respect to a change in matric suction |

: | The initial porosity |

: | The discrete point |

: | A vacuum pressure |

: | The density of air phase |

: | The vertical loading |

: | 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 |

: | The initial degree of saturation |

: | Settlement |

: | The radial displacement |

: | Time |

: | The time factor |

: | The excess pore-air pressure |

: | Atmospheric pressure |

: | The initial excess pore-air pressure |

: | The excess pore-water pressure |

: | 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 |

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