Abstract

A nonlinear flow relationship, which assumes that the fluid flow in the soil skeleton obeys the Hansbo non-Darcian 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 time-dependent 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 one-dimensional 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 [47]. 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 non-Darcy flow models to be relatively good for Xiaoshan clay. Thus, the influence of non-Darcian 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 [1013] conducted the pioneering and systematic work in the application of non-Darcian 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 one-dimensional consolidation to consider the Hansbo non-Darcian flow model. In recent years, Xie et al. [16] and Li et al. [17, 18] published a series of papers regarding the analyses of one-dimensional consolidation with non-Darcian flow. In general, studies that consider non-Darcian 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 non-Darcian 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 three-dimensional consolidation theory to incorporate non-Darcian 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 non-Darcy 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 non-Darcian 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 non-Darcian 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 [2325]: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 two-phase 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 8-node 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 strain-displacement 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 non-Darcian 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 non-Darcian 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 (non-Darcy finite element analysis) was developed to analyze the consolidation process for both natural and vertical drained foundations (either single-drain 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 single-drain 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 soil-drain 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 program will be verified against the theoretical solutions based on both Darcy flow and non-Darcy 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 non-Darcian 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 (non-Darcy 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.

From the plots in Figure 5, it can be found that the results calculated by FEM are roughly accordant with the closed-form 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 non-Darcy 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 non-Darcy flow is given as follows [8, 17]:

Equation (38) is solved numerically using both FEM and FDM, and the Crank-Nicolson 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.

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 eight-node 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 non-Darcy 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.

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 non-Darcy 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 non-Darcy 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 non-Darcy 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 time-dependent settlement analysis of normally consolidated soft soil foundations under static loading. The presented model assumes that fluid flow follows the Hansbo non-Darcian 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 non-Darcian 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 non-Darcian 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 non-Darcian 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.