Research Article  Open Access
Yuebao Deng, Ganbin Liu, Rongyue Zheng, Kanghe Xie, "Finite Element Analysis of Biot’s Consolidation with a Coupled Nonlinear Flow Model", Mathematical Problems in Engineering, vol. 2016, Article ID 3047213, 13 pages, 2016. https://doi.org/10.1155/2016/3047213
Finite Element Analysis of Biot’s Consolidation with a Coupled Nonlinear Flow Model
Abstract
A nonlinear flow relationship, which assumes that the fluid flow in the soil skeleton obeys the Hansbo nonDarcian flow and that the coefficient of permeability changes with void ratio, was incorporated into Biot’s general consolidation theory for a consolidation simulation of normally consolidated soft ground with or without vertical drains. The governing equations with the coupled nonlinear flow model were presented first for the force equilibrium condition and then for the continuity condition. Based on the weighted residual method, the finite element (FE) formulations were then derived, and an existing FE program was modified accordingly to take the nonlinear flow model into consideration. Comparative analyses using established theoretical solutions and numerical solutions were completed, and the results were satisfactory. On this basis, we investigated the effect of the coupled nonlinear flow on consolidation development.
1. Introduction
The theory of consolidation, which forecasts the timedependent deformation of a soil foundation, is one of basic issues in the field of soil mechanics [1]. Terzaghi [2] made use of the effective stress principle and first put forward the onedimensional consolidation theory in the early 1920s. Biot [3] further studied the interaction between the dissipation of excess pore pressures and the deformation of the soil skeleton and introduced a coupled consolidation theory. Since then, Biot’s consolidation theory and its modifications have been used extensively, which has made more precise analyses much more possible [4–7]. It can be seen from these studies that Darcy’s flow law, which assumes a linear correlation between the flow and hydraulic gradient, has been commonly used for its simplicity.
The validity of Darcy’s law has been widely recognized; however, it has also been noted that the flow deviates from Darcy’s law in the case of large plastic sticky soil with a small hydraulic gradient. Hansbo [8] reported that there are several experimental cases indicating that the flow may correlate to an exponential trend with the gradient at low gradient values. In addition, Qi et al. [9] also found the fit results with nonDarcy flow models to be relatively good for Xiaoshan clay. Thus, the influence of nonDarcian flow, which was utilized to depict the deviation of water flow from Darcy’s law, has been assessed by a series of studies in the simulation of soft ground consolidation behavior. Hansbo [10–13] conducted the pioneering and systematic work in the application of nonDarcian flow to the uncoupled consolidation problem. His studies have aroused the attention of several researchers in regard to this issue. Pascal et al. [14] first considered threshold hydraulic gradient in a 1D consolidation problem and obtained the solution using the finite different method (FDM). Dubin and Moulin [15] extended Terzaghi’s onedimensional consolidation to consider the Hansbo nonDarcian flow model. In recent years, Xie et al. [16] and Li et al. [17, 18] published a series of papers regarding the analyses of onedimensional consolidation with nonDarcian flow. In general, studies that consider nonDarcian flow for multidimensional consolidation analyses have been relatively rare. Teh and Nie [19] derived the finite element (FE) equations for Biot’s axisymmetric consolidation problem based on the principle of virtual work and investigated the effect of nonDarcian flow on soft ground consolidation behavior. However, their work was not suitable for the consolidation analysis of multidrain systems and had not taken the variability of soil permeability into consideration.
Based on these previous studies, this paper aims to extend Biot’s general threedimensional consolidation theory to incorporate nonDarcian flow, together with the consideration of the variation of the soil permeability during the consolidation process. Governing equations for the force equilibrium condition and the continuity condition will be presented first, and the finite element method (FEM) will be utilized here to solve these equations. Some specific techniques in the programming process for the newly derived formulas will be stated in detail. The validation of the program will be performed for two cases.
2. A Coupled Nonlinear Flow Model
The fluid flow in the pore media obeys the Hansbo nonDarcy flow law (see Figure 1), which can be described as follows [8]:where is flow velocity; is permeability coefficient of the linear flow at high gradient; is a permeability parameter that is applied in the exponential part; is the exponent of exponential flow at low gradient; is a threshold gradient; and represents the gradient required to overcome the maximum binding energy of mobile pore water [11]. Several studies reported that the value of was within 1.0~2.0 and the typical value was 1.5, whereas the value of ranged in 0~40 [8]. Note here that when and , (1) degenerates to Darcy’s linear flow law; that is, .
Combined with the conditions of the continuity and differentiability at , the following expressions can be derived from (1) [19]:
In addition, a nonconstant permeability is assumed for soft soil. The coefficient of permeability, , is one of the most important parameters for studying nonlinear properties of soil, and a linear correlation between the void ratio and is valid for the most common clays [20], which can be expressed as follows:where and are the initial void ratio and the initial permeability coefficient, respectively, and is seepage index.
There is a linear relationship in  for normally consolidated soft soils according to the oedometer test. That is,where is compressive index; is the mean effective stress; is initial mean effective stress; the mean effective stress ; , , and are the effective stress in the , , and directions, respectively.
Combining (3) and (4), it can be obtained thatwhere ; is in 0.5~2, and the common value is in 0.5~1.0 [20].
A coupled nonlinear flow model, which takes the nonconstancy of soil permeability and nonDarcian flow law into account, is then formed by a combination of (2) and (5). Here it should be pointed out that the settlement of soft soil foundations will be affected by a series of other factors, such as the creep (viscous) behavior of clayey soil (e.g., [21, 22]), the overconsolidated condition, and loading pattern. This paper mainly focuses on the study of the influence of varied soil permeability and nonDarcian flow on consolidation process.
3. Governing Equations
The basic assumptions on the coupled consolidation problem are the same as Biot’s general consolidation theory except for the above newly presented flow model. Equations for force equilibrium and flow continuity relations are presented here for comprehensive consideration.
3.1. Force Equilibrium Equation
The compact form of the force equilibrium equation, which is the same as it is for a linear flow model, can be expressed as follows [23–25]:where the displacement vector ; , , and are displacements in the , , and directions, respectively; is the excess pore pressure (EPP); vector ; the elastic matrix about the isotropic linear elastic material and the differential operator are the same as those for a homogeneous element.
3.2. Flow Continuity Equation
Figure 2 shows a typical element with widths and and thickness . To satisfy the continuity condition in saturated soil, the net flow rate should be equal to the rate of volume change [3]. Consider where , , and are flow velocities in the , , and directions, respectively; is volumetric strain. The matrix form of (7) can be rewritten as follows:where the flow velocity vector ; the strain vector .
Substituting the flow velocity of three directions into (8), the continuity equation can be obtained. Note here that it is necessary to consider the flow velocity components individually because the correlation between the flow velocity and hydraulic gradient in each direction may follow either a linear or an exponential relationship. There are nine possible conditions when comparing the magnitudes of hydraulic gradients in the , , and directions with . To simplify the following derivation process, symbol is introduced to denote the , , or direction.
For any direction , it can be deduced from (2) thatin which (the subscript represents , , or ) indicates the direction of the hydraulic gradient. When the pore water flows to the positive direction of axis (see Figure 2), then the direction of the hydraulic gradient is positive and is +1; otherwise, it is −1. Additionally, indicates the permeability coefficient in the , , or direction. According to (5),
The derivation of (9) can be written as follows:where is the unit weight of water.
According to the study of Teh and Nie [19], item can be computed approximately by using the value of at the previous time step. Combined with (10), (11) can be rewritten in a uniform expression as follows:in whichSubstituting (12) into (8), the matrix form of flow continuity equation can be obtained:where is flow control matrix; is permeability matrix; andin which, , , and can be obtained from (13), respectively; , , and are the initial permeability coefficients in the , , and directions, respectively.
Note here that if parameter is set as 0, then the matrix degenerates to a unit matrix; if the value of is set as 0, the matrix for nonconstant permeability degenerates to the matrix for the constant case. Further, the flow continuity equation, that is, (14), will be reduced to Biot’s general consolidation equation.
3.3. Boundary Conditions
As for the coupled consolidation problem of flow in saturated twophase porous media, there are four categories of boundary conditions: the displacement boundary, the force boundary, the EPP boundary, and the flow velocity boundary. The first two boundary conditions are the same as general mechanical problems, while the EPP boundary usually refers to the zero value of EPP on the drainage surface. In addition, for a typical soil unit (see Figure 2), the flow velocity, , in the normal direction of a random boundary surface can be expressed as follows:where , , and denote the direction cosine of angles between surface and the positive direction of , , and axis, respectively. A flow velocity of zero () is assumed for completely impervious surfaces.
4. Finite Element Approximations
The Galerkin weighted residual method (e.g., [23, 25]) is adopted to transform the above governing equations into a form suitable for obtaining a solution by FEM. In addition, the spatial 8node element is utilized for spatial discretization (see Figure 3). Each node in this element has four degrees of freedom, that is, , , , and .
The FE approximations of displacements and EPP are given bywhere , , and are approximate solutions of displacement in the , , and directions, respectively; is the approximate solution of EPP; the nodal displacement vector is defined by , , ; , , and are nodal displacements in the three directions, respectively; nodal EPP vector is defined by ; is the nodal EPP; ; denotes the cell matrix; ; is the shape function andwhere , , and are local coordinates of each node.
4.1. Equilibrium Equation in FE Formulation
According to the weighted residual method, equations eliminating the residuals of the equilibrium equation and the force boundary equation can be obtained when selecting as the weight function. After reorganization, the following equilibrium equation can be obtained, which is the same as it is for the linear flow condition:in whichwhere , , and are the stiffness matrix, the coupling matrix, and the equivalent incremental nodal load vector, respectively; , , and are incremental vectors of the displacement, EPP, and force, respectively; the straindisplacement matrix , is a matrix with 6 × 3 order, and , ; , , and are surface incremental forces in the , , and directions, respectively. In addition, and are the domain and boundary surface of the operating element, respectively.
4.2. Continuity Equation in FE Formulation
Residuals of the continuity equation and flow velocity boundary condition areEquations eliminating the residuals are defined by the following expressions, respectively:
First, substituting (18) into (23), we havewhere . Combining the integration method and the residual of the flow velocity boundary, it can be obtained thatwhere . Combining (23), (25), and (26), we haveEquation (27) is established from node 1 to node 8. Thus, there are eight equations when (27) is expanded to each node. These equations can be expressed in a compact form as follows:in whichwhere is the flow matrix; is the nodal equivalent flow vector; and specifically, for the impermeable boundary, .
Integrating both sides of (28) with respect to time from to , we havewith the following approximation expression:The incremental form of continuity equation can be obtained as follows:in whichwhere is the incremental nodal equivalent flow vector and is the time increment. Additionally, is an integration constant. A value of is required to maintain the stability of the integration schemes [26].
Equations (20) and (32) are the FE formulations of Biot’s consolidation theory considering the coupled nonlinear flow model in the current study. Because the mathematical derivation of each item for general consolidation theory has been described elsewhere (e.g., [23, 25]), just the evaluation of matrices for nonlinear flow conditions will be discussed in the following section.
4.3. Matrix Evaluation and Programming
Compared with the element consolidation matrix for the linear flow model, matrices and are the same, while the flow matrix is different, leading to a different value for . Matrix evaluations of and are the same as it is for general consolidation theory but with an extra process for calculating the permeability matrix and the flow control matrix .
According to (10) and (16), the permeability matrix at various time steps can be evaluated by iterative computation. First, the initial permeability matrix is formed with the initial condition of soil permeability. Simultaneously, the initial stress field can be derived with the geostatic stress field and the lateral pressure coefficient, and then can be calculated. Substituting matrix into the FE iteration, the stress field of this step is obtained, followed by ascertaining using a similar approach. Substituting into (16), the flow matrix is obtained. The rest can be deduced by analogy.
The evaluation of matrix , which can be considered an adjusting matrix for determining the flow obeying either Darcy’s law or nonDarcian law, depends on the magnitude comparison between and . It will take on different forms under different gradient conditions. As referenced in the previous section, matrix can be evaluated approximately by using the value of the previous time step. Stated in detail, under the local coordinate (see Figure 3), the hydraulic gradient of the next time step, (the superscript denotes the time step), can be calculated by the nodal EPP at current time step . That is,where , , and are half sizes of the element in , , and directions, respectively; see Figure 3. The ascertained values of , can be calculated by (13) after a comparison between and is completed. Then the flow control matrix of the next time step is obtained.
Note here that, for Darcy’s flow, which is a special case of nonDarcian flow, matrix degenerates into a unit matrix. Additionally, two approaches can be used to ascertain the initial flow control matrix . Considering that the hydraulic gradients decrease gradually from the maximum values at the initiation of consolidation to zero when consolidation has fully completed, the first approach is setting matrix as a unit matrix directly. This approach denotes that the flow at the initial condition follows the trend for Darcy flow [19]. Another approach is setting as a unit matrix at first and then evaluating it using an iterative operation. The iterative process is the same as the process described above, which uses the EPP of the previous step to calculate the matrix of the current step. The flow control matrix for various time steps is also computed by iteration in this approach.
Based on an existing numerical code [25, 27, 28], which had been verified and applied in a series of studies, a new program NDFEA (nonDarcy finite element analysis) was developed to analyze the consolidation process for both natural and vertical drained foundations (either singledrain or multidrain conditions). To check for the correct implementation of the algorithms and formulations, two comparative analyses for the analytical and numerical solutions will be performed.
5. Verification
5.1. Case 1
The first case is a fully penetrating singledrain foundation with one clay layer, which has also been used by Hird et al. [29] and Teh and Nie [19]. The cylindrical influence zone diameter is denoted by m; the smear zone diameter is m; and the drain diameter is m. According to the principle of equation area, the cylindrical zone is translated into a block zone [25]; see Figure 4. The corresponding sizes , , and for the equivalent system are calculated by a simple function; for example, . Foundation depth or drain length m. The diagram in Figure 4 illustrates the FE discretization and boundary conditions: the total number of mesh elements is 3179; displacements at the outer boundaries are fixed in the horizontal direction, with only vertical displacement allowed; displacement at the top surface is free, whereas the bottom surface is totally fixed; regarding the drainage boundary, the periphery of the unit cell is set as impermeable, except the drain unit on the top surface, which is assumed to have been drained. Additionally, the soildrain system is assumed to be elastic, with regard to and , and is set for 1D deformation analysis. A surcharge loading of is also simulated by instantaneous application to the upper boundary, and the permeability parameters , , , and are used to capture the permeability differences of the different parts, respectively. The values of these parameters are listed in Table 1, in which four calculation conditions were set to examine the newly developed program. The above settings are aiming to be accordant with the assumptions of the available analytical solution.
 
: the largest difference between the FE results and the Hansbo solution [11] for nonDarcy flow under different calculation conditions. 
The program will be verified against the theoretical solutions based on both Darcy flow and nonDarcy flow. According to the Hansbo consolidation theory, including the effects of both smear and well resistance [8], the average degree of consolidation at a given elevation can be written as follows. For Darcy flow, For nonDarcian flow, in which is consolidation time; for instantaneous applied load, the initial EPP is . Additionally,
FE analyses with different iteration times (e.g., iteration time = 1, 2, 3, 6) for the time step cycling were conducted to check the influence of the hydraulic gradient approximation on numerical results. As Figure 5(a) shows, the numerical result of the second approach with six iteration times gets close to the theoretical solution (nonDarcy condition) when compared with the result with the first approach (iteration time = 1), but this improvement is very limited. The following FE results for different conditions are based on the first approach, in which and are directly ascertained using previous time step values without iteration.
(a) Condition 1
(b) Condition 2
(c) Condition 3
(d) Condition 4
From the plots in Figure 5, it can be found that the results calculated by FEM are roughly accordant with the closedform solutions. The last column in Table 1 shows the largest differences for different calculation conditions. Combining Figure 5 and Table 1, the following conclusions can be obtained: (1) the largest differences for the three calculation conditions are between 2 and 10%; (2) FE solutions are mostly between the analytical results based on Darcy flow and those based on nonDarcy flow; (3) these differences increase with the increase of the nonlinear control parameter . The cause of the difference is stated as follows: with the increase of , the nonlinearity of water flow becomes more obvious (see Figure 1); then, there is the approximation that computing the hydraulic gradient of the next step with pore water pressure of the present step leads to a more apparent difference between the actual value and calculating result by (34). However, the errors are not significant and the comparison results here matched perfectly with the previous comparative analyses by Teh and Nie [19].
5.2. Case 2
Another validation case will be conducted for the analysis of the 1D consolidation problem, which is a saturated homogeneous soft clay layer with a thickness of ( m). The top boundary is free and drained, while the bottom boundary is fixed and impervious. A surcharge ( = 100 kPa) is instantaneously applied at the ground surface. The permeability coefficient and mechanical parameters are assumed to be constant. In addition, the equation governing the 1D consolidation of soil with nonDarcy flow is given as follows [8, 17]:
Equation (38) is solved numerically using both FEM and FDM, and the CrankNicolson difference scheme is adopted here for its relative stability. The detailed scheme with the FDM has been presented by Li et al. [17], and it has been validated with good reliability for calculation accuracy.
Two calculation conditions have been conducted for this case. The parameters and their values are listed in Table 2. The comparison of results between the two numerical methods is shown in Figure 6. Additionally, Table 2 presents the largest differences under different calculation conditions. It can be noted that the largest difference between the two methods is less than 3%. These contrasting analysis results are satisfactory.
 
: the largest difference between the results by FEM and FDM under different conditions. 
6. Effect of Nonlinear Flow on Consolidation
6.1. Description of the Problem
An equivalent block unit cell will be studied to investigate the effect of nonlinear flow on consolidation. Information about the geometries, boundary conditions, and FE discretization for the case can also be seen in Figure 4. The sizes of , , and are listed in Table 3, respectively. Boundary conditions are the same as for the first case for validation, except for the drainage condition for the top surface, which is assumed to be drained. The total number of 1690 spatial eightnode block elements is employed in the case after an initial mesh sensitivity study was performed. In addition, the soil and drain are still modeled with a linear elastic model, to highlight the influence of the nonlinear flow. Here note that (4), which depicted a nonlinear compression property of normally consolidation soil, was just used for the derivation of a nonlinear flow expression, that is, (5). The permeability property of the soil is assumed to generally change with effective stress, and the initial values , , , and are listed in Table 3, in which parameters corresponding to the initial stress fields and are also presented.

As shown in Table 4, nine calculation conditions will be conducted in which is set as 1.5 because this value is validated by several studies, whereas and are assumed in a range of 0~20 and 0–1.0, respectively. Note here that the above nonlinear flow parameters are just applicable for the soil, while the permeability property of the drain is assumed constant. In addition, loading size is related to the hydraulic gradient. The higher the value of , the more the significant change of hydraulic gradient in the dissipation of EPP, and thus the influence of nonDarcy flow will be more obvious. In this consideration, a value range of 10 to 100 kPa is adopted. Among the nine calculation conditions, the first case is the linear flow condition, which is used to make a comparison with the other conditions.

6.2. Parametric Studies
The consolidation data are plotted in terms of average degree of consolidation (by (38)) with time in Figure 7. The time taken to achieve a fixed degree of consolidation and the largest differences for different calculation results are listed in Tables 5 and 6, respectively.
 
and : the time taken to achieve % and %, respectively. 
 
: the largest difference between the FE results for nonlinear flow under different conditions with linear flow model. 
(a) Assess the influence of nonlinear flow parameter
(b) Assess the influence of surcharge intensity
(c) Assess the influence of nonDarcy flow parameter
The influence of the nonlinear flow parameter (), which reflects the reduction rate of permeability, has been estimated by result comparison between conditions 2, 3, and 4; see Figure 7(a). As expected, the consolidation rate slows down when the variability of is considered, and the retardation on consolidation becomes significant with an increasing . As shown in Table 5, under the nonDarcy flow conditions with , 0.5, and 1.0, consolidation times () are approximately 23 days (102), 29 days (131), and 39 days (190), respectively. The time for the linear flow condition is 22 days (92); in addition, the largest differences for the three conditions are in the range of 2 to 15%; see Table 6. Among them, the biggest value is 14.58% under a condition where .
Conditions 3, 5, and 6 are used to analyze the effect of loading size by fixing and varying . The comparisons are shown in Figure 7(b) and Table 5. Conclusions can be obtained that (1) with reducing the consolidation rate is retarded and the consolidation time is prolonged; (2) the influence of on consolidation is insignificant when changes from 100 kPa to 50 kPa, whereas sharp retardation happens when further reduces to 10 kPa. It is to be expected that the fluid flow follows nonDarcy flow when is reduced, which results in a low flow velocity. Thus, it can be foreseen that as the load further reduces, the consolidation rate will slow further, and the EPP gradients will be less than the initial hydraulic gradient. It can be observed from Table 6 that the largest difference between condition 6 ( = 10 kPa) and condition 1 is 13.5%.
Figure 7(c) shows the comparison of consolidation results with different nonDarcy flow parameters. For conditions 3, 7, and 8, varying from 5 to 20, consolidation times () are approximately 29 days (131), 31 days (148), and 33 days (176), respectively. These data lend support to the observation that as the values of and/or increase, consolidation development slows down [19]. However, this effect is not so obvious as to be attributed simply to the influence of the nonlinear flow parameter. The consolidation curve for condition 9 (; ; ) is also plotted in Figure 7(c). It can be obtained from the comparison of conditions 8 and 9 that as the nonlinear flow degree increases, the consolidation rate is further slowed down, and the effect of is more obvious when compared with the influence of . , for example, changes from 29 days to 33 days for the previous three conditions, 3, 7, and 8, whereas it is prolonged to 43 days, which is almost twice as long as for the results based on linear flow formulation. The largest difference in the results of the degree of consolidation under different conditions also validates the above conclusion. Table 6 shows that the largest difference between conditions 1 and 9 is 22%.
7. Conclusions
Biot’s general consolidation theory was extended to incorporate a coupled nonlinear flow model for the timedependent settlement analysis of normally consolidated soft soil foundations under static loading. The presented model assumes that fluid flow follows the Hansbo nonDarcian flow and that the permeability coefficient of soil decreases with an exponential change of effective stress. The governing differential equations and the FE formulations were derived, and a programming code was developed to solve the obtained equations. We implemented two calculation cases to check the rationality of the method, and we found both cases to demonstrate that the consolidation results for the two comparative methods were in good agreement. Finally, we investigated the effect of nonlinear flow on consolidation, and we drew the following conclusions:(1)The consolidation rate slows down when coefficient or permeability variability is considered. With an increase in the value of (, value range is 0.5–2.0), the deduction of soil permeability in the EPP dissipation process becomes obvious, and the consolidation rate slows down. The effect of on the consolidation development is somewhat significant.(2)The hydraulic gradient of EPP in the consolidation process becomes uniform due to the load size decrease, which results as an apparent effect of nonDarcian flow on the degree of consolidation. Thus, with the loading size decrease, consolidation rate slows down and the time taken to achieve a fixed degree of consolidation increases. This trend of decreasing consolidation rate with decreasing loading size has a tendency to accelerate.(3)Increasing the values of nonDarcian flow parameters ( and ) can cause the consolidation rate to slow down. In relative terms, this effect is not as obvious as the influence of parameter on consolidation behavior.
Nomenclature
:  Half sizes of the element in , , and directions 
:  Widths of the equivalent unit cell, the smear zone, and the drain 
:  Diameters of the equivalent unit cell, the smear zone, and the drain 
:  Compressive index 
:  Seepage index 
:  Initial void ratio 
:  Young’s modulus 
:  Depth of the soft soil foundation 
:  Hydraulic gradient 
:  Threshold hydraulic gradient in the Hansbo nonDarcian flow model 
:  Permeability coefficient 
:  Initial permeability coefficient 
:  Permeability coefficients in the horizontal and in the vertical directions 
:  Initial permeability coefficients in the horizontal and in the vertical directions 
:  Permeability coefficients of the smear zone and the drain 
:  Initial permeability coefficients of the smear zone and the drain 
:  Length of drain 
:  The exponent of exponential flow at low gradient 
:  Excess pore water pressure 
:  Initial excess pore water pressure 
:  Mean effective stress, 
:  Initial mean effective stress 
:  The instantaneous load applied at the top surface of the ground 
:  Displacements in the three directions , , and 
:  Ratio 
:  Volumetric strain 
:  Unit weight of water 
:  The effective weight and friction angle of soil 
:  Permeability parameter applied in the exponential part 
:  Integration constant 
:  Poisson’s ratio 
:  The effective stress in the , , and directions 
:  The three directions of the local coordinate system. 
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was sponsored by K. C. Wong Magna Fund in Ningbo University and by the National Natural Science Foundation of China (nos. 51308309 and 51278256), which are gratefully acknowledged.
References
 L. Ho, B. Fatahi, and H. Khabbaz, “Analytical solution for onedimensional consolidation of unsaturated soils using eigenfunction expansion method,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 10, pp. 1058–1077, 2014. View at: Publisher Site  Google Scholar
 K. Terzaghi, Theoretical Soil Mechanics, John Wiley & Sons, New York, NY, USA, 1943. View at: Publisher Site
 M. A. Biot, “General theory of threedimensional consolidation,” Journal of Applied Physics, vol. 12, no. 2, pp. 155–164, 1941. View at: Publisher Site  Google Scholar
 R. W. Lewis and B. A. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, John Wiley & Sons, New York, NY, USA, 2nd edition, 1998.
 J. Huang, D. V. Griffiths, and G. A. Fenton, “Probabilistic analysis of coupled soil consolidation,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 136, no. 3, pp. 417–430, 2010. View at: Publisher Site  Google Scholar
 C. Menéndez, P. J. Nieto, F. A. Ortega, and A. Bello, “Nonlinear analysis of the consolidation of an elastic saturated soil with incompressible fluid and variable permeability by FEM,” Applied Mathematics and Computation, vol. 216, no. 2, pp. 458–476, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 L. Ho and B. Fatahi, “Analytical solution for the twodimensional plane strain consolidation of an unsaturated soil stratum subjected to timedependent loading,” Computers and Geotechnics, vol. 67, pp. 1–16, 2015. View at: Publisher Site  Google Scholar
 S. Hansbo, “Experience of consolidation process from test areas with and without vertical drains,” in Ground Improvement—Case Histories, B. Indraratna and J. Chu, Eds., vol. 3, pp. 3–49, Elsevier, London, UK, 2005. View at: Google Scholar
 T. Qi, K.H. Xie, A.F. Hu, and Z.Q. Zhang, “Laboratorial study on nonDarcy seepage in Xiaoshan clay,” Journal of Zhejiang University: Engineering Science Edition, vol. 41, no. 6, pp. 1023–1028, 2007 (Chinese). View at: Google Scholar
 S. Hansbo, Consolidation of Clay with Special Reference to Influence of Vertical Drains, vol. 18 of Swedish Geotechnical Institute Proceeding, Swedish Geotechnical Institute, Stockholm, Sweden, 1960.
 S. Hansbo, “Aspects of vertical drain design: Darcian or nonDarcian flow,” Geotechnique, vol. 47, no. 5, pp. 983–992, 1997. View at: Publisher Site  Google Scholar
 S. Hansbo, “Consolidation equation valid for both Darcian and nonDarcian flow,” Geotechnique, vol. 51, no. 1, pp. 51–54, 2001. View at: Publisher Site  Google Scholar
 S. Hansbo, “Deviation from Darcy's law observed in onedimensional consolidation,” Geotechnique, vol. 53, no. 6, pp. 601–605, 2003. View at: Publisher Site  Google Scholar
 F. Pascal, H. Pascal, and D. W. Murray, “Consolidation with threshold gradients,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 5, no. 3, pp. 247–261, 1981. View at: Publisher Site  Google Scholar
 B. Dubin and G. Moulin, “Influence of a critical gradient on the consolidation of clays,” in Consolidation of Soils: Testing and Evaluation, ASTM STP 892, pp. 354–377, American Society for Testing and Materials, West Conshohocken, Pa, USA, 1986. View at: Google Scholar
 K.H. Xie, K. Wang, Y.L. Wang, and C.X. Li, “Analytical solution for onedimensional consolidation of clayey soils with a threshold gradient,” Computers and Geotechnics, vol. 37, no. 4, pp. 487–493, 2010. View at: Publisher Site  Google Scholar
 C.X. Li, K.H. Xie, and K. Wang, “Analysis of 1D consolidation with nonDarcian flow described by exponent and threshold gradient,” Journal of Zhejiang University—Science A, vol. 11, no. 9, pp. 656–667, 2010. View at: Publisher Site  Google Scholar
 C.X. Li, K.H. Xie, A.F. Hu, and B.X. Hu, “Onedimensional consolidation of doublelayered soil with nonDarcian flow described by exponent and threshold gradient,” Journal of Central South University of Technology, vol. 19, no. 2, pp. 562–571, 2012. View at: Publisher Site  Google Scholar
 C. I. Teh and X. Y. Nie, “Coupled consolidation theory with nonDarcian flow,” Computers and Geotechnics, vol. 29, no. 3, pp. 169–209, 2002. View at: Publisher Site  Google Scholar
 F. Tavenas, M. Tremblay, G. Larouche, and S. Leroueil, “Insitu measurement of permeability in soft clays,” in ASCE Special Conference on Use of Insitu Tests in Geotechnical Engineering, pp. 1034–1048, ASCE, 1986. View at: Google Scholar
 B. Fatahi, T. M. Le, M. Q. Le, and H. Khabbaz, “Soil creep effects on ground lateral deformation and pore water pressure under embankments,” Geomechanics and Geoengineering, vol. 8, no. 2, pp. 107–124, 2013. View at: Publisher Site  Google Scholar
 T. M. Le, B. Fatahi, and H. Khabbaz, “Numerical optimisation to obtain elastic viscoplastic model parameters for soft clay,” International Journal of Plasticity, vol. 65, pp. 1–21, 2015. View at: Publisher Site  Google Scholar
 O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structure Mechanics, Elsevier, Singapore, 6th edition, 2009.
 R. S. Sandhu and E. L. Wilson, “Finite element analysis of seepage in elastic media,” Proceedings of the American Society of Civil Engineers, vol. 95, no. 3, pp. 641–652, 1996. View at: Google Scholar
 K. H. Xie and J. Zhou, Theory and Application of Finite Element Method in Geotechnical Engineering, Science Press, Beijing, China, 2002.
 J. R. Booker and J. C. Small, “An investigation of the stability of numerical solutions of Biot’s equations of consolidation,” International Journal of Solids and Structures, vol. 11, no. 78, pp. 907–917, 1975. View at: Publisher Site  Google Scholar
 G. X. Zeng, K. H. Xie, and Z. Y. Shi, “Consolidation analysis of sanddrained ground by FEM,” in Proceedings of the 8th Asian Regional Conference of Soil Mechanic and Foundation Engineering, pp. 139–142, Tokyo, Japan, 1987. View at: Google Scholar
 K.H. Xie, Y.B. Deng, K. Wang, and D.Z. Huang, “A 3D composite finite element used to model a vertical drain and its surrounding remoulded soil,” Computers and Geotechnics, vol. 38, no. 8, pp. 1078–1088, 2011. View at: Publisher Site  Google Scholar
 C. C. Hird, I. C. Pyrah, and D. Russell, “Finite element modelling of vertical drains beneath embankments on soft ground,” Géotechnique, vol. 42, no. 3, pp. 499–511, 1992. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2016 Yuebao Deng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.