Axisymmetric Consolidation of Unsaturated Soils by Differential Quadrature Method

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


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

Mathematical Problems in Engineering
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][14][15][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.

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.

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  0 is the initial total volume of the soil, (  / 0 )/ 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   1 is the coefficient of water volume change with respect to a change in the net normal stress ( −   ),    2 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 where 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, (  / 0 )/ 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 −  0 ) 0  0 , where  0 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] where  is the universal air constant (8.314J/mol K),  = 293.15K is the absolute temperature,  = 29 kg/kmol is the molecular mass of air,   0 is absolute pore-air pressure (i.e.,   0 =   0 + 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: where   1 is the coefficient of air volume change with respect to a change in the net normal stress ( −   ) and   2 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 where . Further, two transformed governing equations of water and air phases from ( 3) and ( 7) can be obtained: where

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 9) reflects that the lateral boundary surface is impermeable or impeded, respectively. 0  and  0  are the coefficients of permeability for water and air at lateral boundary, respectively. 0 is the thickness of lateral boundary.
The initial condition can be written as where   0 and   0 are the initial pore-water and pore-air pressure distributions.

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 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 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 − )  1 + (0.5 + )   , (−0.5 ≤  ≤ 0.5) , (14) where  is the integral coordinate of soil layer and  1 and   are the radius values of 1th and th 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  (1)   and  (2)   are the weighting coefficient matrices of the first order and second order derivatives, respectively.
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 0 to , we get the expression of volumetric strain  V : 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 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  = 1 m 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.
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 −3 .When  = 1.0 × 10 0 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 3 shows the radial displacement development with time and Figure 4 describes the settlement of top surface ( = 10 m) 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  = 1 m with time factor .When the height of soil layer  is equal to 4 m, the settlement is almost zero.

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.  2 = −0.61× 10 −4 kPa −1 ,   1 = 0.06 × 10 −4 kPa −1 ,   2 = 0.26 × 10 −4 kPa −1 , and  atm = 101.3kPa.

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   0 = 5 kPa and   0 = 40 kPa. Figure 6(b) presents the linear increasing initial pore-air and pore-water distributions.The lateralskewed 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).
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 porewater 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.6 m to  = 1.2 m, 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  = 0.4 m to  = 1.6 m, the initial pore-water pressure will increase from 30 kPa to 110 kPa and the initial pore-air pressure will increase from 3.75 kPa to 13.75 kPa.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  = 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 porewater 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 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 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.

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   0 = 5 kPa and   0 = 40 kPa are adopted.Four different ratios   /  are considered, that is, 0.1, 1, 10, and 100 with   = 10 −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.2 m under different boundary conditions.First, the excess pore-water and poreair 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 porewater 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.

Conclusion
This paper obtains a general solution to the axisymmetric consolidation of unsaturated soils by using differential quadrature method (DQM) based on Fredlund's onedimensional 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 porewater 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
V : The consolidation coefficient of water   V : The consolidation coefficient of air  ()   : The weighting coefficient  V : Volumetricstrain : The gravitational acceleration : The thickness of soil layer   : The permeability coefficients of air  The density of air phase : The vertical loading  0 : The thickness of lateral boundary   : Theradiusofsoillayer   : 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  V : 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

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

Figure 2 :Figure 3 :
Figure 2: The average degree of consolidation for the water and air phases with time factor  (  /  = 1).

2 : 2 :
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   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  0 : The initial porosity : The discrete point  V : A vacuum pressure :