# Finite Element Analysis of Biot’s Consolidation with a Coupled Nonlinear Flow Model

Authors: Yue-bao Deng; Gan-bin Liu; Rong-yue Zheng; Kang-he Xie
Journal: Mathematical Problems in Engineering
 (2016)

---

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


---

## Body

## 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 [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 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 [10–13] 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 i0 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 Figure1), which can be described as follows [8]:(1)v=κimwheni<il,v=ki-i0wheni≥il,where v is flow velocity; k is permeability coefficient of the linear flow at high gradient; κ is a permeability parameter that is applied in the exponential part; m is the exponent of exponential flow at low gradient; i0 is a threshold gradient; and il represents the gradient required to overcome the maximum binding energy of mobile pore water [11]. Several studies reported that the value of m was within 1.0~2.0 and the typical value was 1.5, whereas the value of il ranged in 0~40 [8]. Note here that when m=1 and i0=0, (1) degenerates to Darcy’s linear flow law; that is, v=ki.Figure 1
The Hansbo non-Darcy flow model.Combined with the conditions of the continuity and differentiability ati=il, the following expressions can be derived from (1) [19]:(2)v=kimmilm-1wheni<il,v=ki-m-1milwheni≥il.In addition, a nonconstant permeability is assumed for soft soil. The coefficient of permeability,k, is one of the most important parameters for studying nonlinear properties of soil, and a linear correlation between the void ratio e and lg⁡k is valid for the most common clays [20], which can be expressed as follows:(3)e=e0+Cklg⁡kk0,where e0 and k0 are the initial void ratio and the initial permeability coefficient, respectively, and Ck is seepage index.There is a linear relationship ine-lg⁡p′ for normally consolidated soft soils according to the oedometer test. That is,(4)e=e0-Cclg⁡p′p0′,where Cc is compressive index; p′ is the mean effective stress; p0′ is initial mean effective stress; the mean effective stress p′=(σx′+σy′+σz′)/3; σx′, σy′, and σz′ are the effective stress in the x-, y-, and z-directions, respectively.Combining (3) and (4), it can be obtained that(5)k=k0p′p0′-α,where α=Cc/Ck; α 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 [23–25]:(6)∂TD∂f-∂TMp=0,where the displacement vector {f}=uvwT; u, v, and w are displacements in the x-, y-, and z-directions, respectively; p is the excess pore pressure (EPP); vector {M}=111000T; the elastic matrix [D] about the isotropic linear elastic material and the differential operator [∂] are the same as those for a homogeneous element.


### 3.2. Flow Continuity Equation

Figure2 shows a typical element with widths dx and dy and thickness dz. To satisfy the continuity condition in saturated soil, the net flow rate should be equal to the rate of volume change [3]. Consider (7)∂vx∂x+∂vy∂y+∂vz∂z=-∂εv∂t=-∂∂t∂u∂x+∂v∂y+∂w∂z,where vx, vy, and vz are flow velocities in the x-, y-, and z-directions, respectively; εv is volumetric strain. The matrix form of (7) can be rewritten as follows:(8)MT∂v=-∂∂tMTε,where the flow velocity vector v=vxvyvzT; the strain vector ε=εxεyεzγxyγyzγzxT.Figure 2
Schematic diagram of the soil unit.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 v and hydraulic gradient i 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 x-, y-, and z-directions with il. To simplify the following derivation process, symbol Ψ is introduced to denote the x-, y-, or z-direction.For any directionΨ, it can be deduced from (2) that(9)vΨ=SΨkΨiΨmmilm-1wheniΨ<il,vΨ=SΨkΨiΨ-m-1milwheniΨ≥ilin which SΨ (the subscript Ψ represents x, y, or z) 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 SΨ is +1; otherwise, it is −1. Additionally, kΨ indicates the permeability coefficient in the x-, y-, or z-direction. According to (5),(10)kΨ=kΨ0p′p0′-α.The derivation of (9) can be written as follows:(11)∂vΨ∂Ψ=kΨγwmilm-1∂p∂Ψm-1∂2p∂Ψ2wheniΨ<il,∂vΨ∂Ψ=kΨγw∂2p∂Ψ2wheniΨ≥il,where γw is the unit weight of water.According to the study of Teh and Nie [19], item ∂p/∂Ψ can be computed approximately by using the value of p at the previous time step. Combined with (10), (11) can be rewritten in a uniform expression as follows:(12)∂vΨ∂Ψ=HΨkΨ0γwp′p0′-α∂2p∂Ψ2in which(13)HΨ=1γwm-1ilm-1∂p∂Ψm-1wheniΨ<il,1wheniΨ≥il.Substituting (12) into (8), the matrix form of flow continuity equation can be obtained:(14)MT∂Hk∂TMp-∂∂tMT∂f=0,where [H] is flow control matrix; [k] is permeability matrix; and(15)H=Hx000Hy000Hz,(16)k=1γwp′p0′-αkx0000ky0000kz0in which, Hx, Hy, and Hz can be obtained from (13), respectively; kx0, ky0, and kz0 are the initial permeability coefficients in the x-, y-, and z-directions, respectively.Note here that if parameterm is set as 0, then the matrix [H] degenerates to a unit matrix; if the value of α is set as 0, the matrix [k] 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 Figure2), the flow velocity, vn, in the normal direction of a random boundary surface S can be expressed as follows:(17)l-vx+m-vy+n-vz=vn,where l-, m-, and n- denote the direction cosine of angles between surface S and the positive direction of x-, y-, and z-axis, respectively. A flow velocity of zero (vn=0) is assumed for completely impervious surfaces.


## 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]:(6)∂TD∂f-∂TMp=0,where the displacement vector {f}=uvwT; u, v, and w are displacements in the x-, y-, and z-directions, respectively; p is the excess pore pressure (EPP); vector {M}=111000T; the elastic matrix [D] about the isotropic linear elastic material and the differential operator [∂] are the same as those for a homogeneous element.


## 3.2. Flow Continuity Equation

Figure2 shows a typical element with widths dx and dy and thickness dz. To satisfy the continuity condition in saturated soil, the net flow rate should be equal to the rate of volume change [3]. Consider (7)∂vx∂x+∂vy∂y+∂vz∂z=-∂εv∂t=-∂∂t∂u∂x+∂v∂y+∂w∂z,where vx, vy, and vz are flow velocities in the x-, y-, and z-directions, respectively; εv is volumetric strain. The matrix form of (7) can be rewritten as follows:(8)MT∂v=-∂∂tMTε,where the flow velocity vector v=vxvyvzT; the strain vector ε=εxεyεzγxyγyzγzxT.Figure 2
Schematic diagram of the soil unit.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 v and hydraulic gradient i 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 x-, y-, and z-directions with il. To simplify the following derivation process, symbol Ψ is introduced to denote the x-, y-, or z-direction.For any directionΨ, it can be deduced from (2) that(9)vΨ=SΨkΨiΨmmilm-1wheniΨ<il,vΨ=SΨkΨiΨ-m-1milwheniΨ≥ilin which SΨ (the subscript Ψ represents x, y, or z) 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 SΨ is +1; otherwise, it is −1. Additionally, kΨ indicates the permeability coefficient in the x-, y-, or z-direction. According to (5),(10)kΨ=kΨ0p′p0′-α.The derivation of (9) can be written as follows:(11)∂vΨ∂Ψ=kΨγwmilm-1∂p∂Ψm-1∂2p∂Ψ2wheniΨ<il,∂vΨ∂Ψ=kΨγw∂2p∂Ψ2wheniΨ≥il,where γw is the unit weight of water.According to the study of Teh and Nie [19], item ∂p/∂Ψ can be computed approximately by using the value of p at the previous time step. Combined with (10), (11) can be rewritten in a uniform expression as follows:(12)∂vΨ∂Ψ=HΨkΨ0γwp′p0′-α∂2p∂Ψ2in which(13)HΨ=1γwm-1ilm-1∂p∂Ψm-1wheniΨ<il,1wheniΨ≥il.Substituting (12) into (8), the matrix form of flow continuity equation can be obtained:(14)MT∂Hk∂TMp-∂∂tMT∂f=0,where [H] is flow control matrix; [k] is permeability matrix; and(15)H=Hx000Hy000Hz,(16)k=1γwp′p0′-αkx0000ky0000kz0in which, Hx, Hy, and Hz can be obtained from (13), respectively; kx0, ky0, and kz0 are the initial permeability coefficients in the x-, y-, and z-directions, respectively.Note here that if parameterm is set as 0, then the matrix [H] degenerates to a unit matrix; if the value of α is set as 0, the matrix [k] 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 Figure2), the flow velocity, vn, in the normal direction of a random boundary surface S can be expressed as follows:(17)l-vx+m-vy+n-vz=vn,where l-, m-, and n- denote the direction cosine of angles between surface S and the positive direction of x-, y-, and z-axis, respectively. A flow velocity of zero (vn=0) 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, u, v, w, and p.Figure 3
Spatial 8-node element.The FE approximations of displacements and EPP are given by(18)f≈f~=u~,v~,w~T=Nδe;p≈p~=N-pe,where u~, v~, and w~ are approximate solutions of displacement in the x-, y-, and z-directions, respectively; p~ is the approximate solution of EPP; the nodal displacement vector is defined by δe=δ1Tδ2T⋯δ8TT, δi=uiviwiT, i=1,2,…,8; ui, vi, and wi are nodal displacements in the three directions, respectively; nodal EPP vector is defined by pe=p1p2⋯p8T; pi is the nodal EPP; N=N1IN2I⋯N8I; I denotes the cell matrix; N-=N1N2⋯N8; Ni is the shape function and(19)Ni=181+ξξi1+ηηi1+ζζi;i=1,2,…,8,where ξi, ηi, and ζi 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 selectingNi 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:(20)KeΔδe+KcΔpe=ΔRfein which(21)Ke=∭VeBTDBdV;Kc=∭VeBTMN-dV;ΔRfe=∬SeNTΔFedS,Δδe=Δδ1TΔδ2T⋯Δδ8TT;Δpe=Δp1Δp2⋯Δp8T;ΔFe=ΔFx,ΔFy,ΔFzT,where [Ke], [Kc], and ΔRf are the stiffness matrix, the coupling matrix, and the equivalent incremental nodal load vector, respectively; Δδ, Δp, and ΔF are incremental vectors of the displacement, EPP, and force, respectively; the strain-displacement matrix B=-∂N=Β1B2⋯B8, Bi is a matrix with 6 × 3 order, and Bi=-∂Ni, i=1,2,…,8; ΔFx, ΔFy, and ΔFz are surface incremental forces in the x-, y-, and z-directions, respectively. In addition, Ve and Se 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 are(22)RIqe=∂∂tMT∂f~-MT∂Hk∂TMp~,RBqe=l-vx+m-vy+n-vz-vn.Equations eliminating the residuals are defined by the following expressions, respectively:(23)∭VeNi∂∂tMT∂f~-MT∂Hk∂TMp~dV=0;i=1,2,…,8,(24)∬SeNil-vx+m-vy+n-vz-vndS=0;i=1,2,…,8.First, substituting (18) into (23), we have(25)∭VeNi∂∂tMT∂f~dV=∭VeNiMT∂Nδ˙dV,where δ˙e=∂δe/∂t. Combining the integration method and the residual of the flow velocity boundary, it can be obtained that(26)-∭VeNiMT∂Hk∂TMp~dV=∬SeNivndS+∭Ve∇TNiHk∇N-pedV,where ∇=∂/∂x∂/∂y∂/∂zT. Combining (23), (25), and (26), we have(27)∭VeNiMT∂Nδ˙dV+∭Ve∇TNiHk∇N-pedV=-∬SeNivndS;i=1,2,…,8.Equation (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:(28)KcTδ˙e-Kspe=Rqe,in which(29)Ks=∭VeBsTHkBsdV;Bs=∂∂x∂∂y∂∂zTN1N2⋯N8;Rqe=∬SeN-TvndS,where [Ks] is the flow matrix; Rqe is the nodal equivalent flow vector; and specifically, for the impermeable boundary, Rqe=0.Integrating both sides of (28) with respect to time from tn to tn+1, we have(30)∫tjtj+1KcTδ˙edt-∫tjtj+1Kspedt=∫tjtj+1Rqedtwith the following approximation expression:(31)∫tjtj+1pedt=Δtpne+θΔpe;∫tjtj+1Rqedt=ΔtRqne.The incremental form of continuity equation can be obtained as follows:(32)KcTΔδe-θΔtKsΔpe=ΔRpein which(33)ΔRpe=ΔtRqje+Kspje,where ΔRpe is the incremental nodal equivalent flow vector and Δt is the time increment. Additionally, θ is an integration constant. A value of θ≥0.5 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[Ke] and [Kc] are the same, while the flow matrix [Ks] is different, leading to a different value for ΔRpe. Matrix evaluations of [Ks] and ΔRpe are the same as it is for general consolidation theory but with an extra process for calculating the permeability matrix [k] and the flow control matrix [H].According to (10) and (16), the permeability matrix at various time steps ki can be evaluated by iterative computation. First, the initial permeability matrix k0 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 p0′ can be calculated. Substituting matrix k0 into the FE iteration, the stress field of this step is obtained, followed by ascertaining p1′ using a similar approach. Substituting p1′ into (16), the flow matrix k1 is obtained. The rest can be deduced by analogy.The evaluation of matrix[H], 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 iΨ and il. It will take on different forms under different gradient conditions. As referenced in the previous section, matrix [H] can be evaluated approximately by using the p value of the previous time step. Stated in detail, under the local coordinate (see Figure 3), the hydraulic gradient of the next time step, ij+1 (the superscript j denotes the time step), can be calculated by the nodal EPP at current time step pj. That is,(34)ixj+1=1γw∂pj∂x=p1j+p4j+p5j+p8j-p2j+p3j+p6j+p7j8γwa,iyj+1=1γw∂pj∂y=p1j+p2j+p5j+p6j-p3j+p4j+p7j+p8j8γwb,izj+1=1γw∂pj∂z=p1j+p2j+p3j+p4j-p5j+p6j+p7j+p8j8γwc,where a, b, and c are half sizes of the element in x-, y-, and z-directions, respectively; see Figure 3. The ascertained values of iΨ, HΨ can be calculated by (13) after a comparison between iΨ and il is completed. Then the flow control matrix [H] of the next time step is obtained.Note here that, for Darcy’s flow, which is a special case of non-Darcian flow, matrix[H] degenerates into a unit matrix. Additionally, two approaches can be used to ascertain the initial flow control matrix H0. 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 H0 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 H0 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 [H] of the current step. The flow control matrix Hi 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.


## 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 selectingNi 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:(20)KeΔδe+KcΔpe=ΔRfein which(21)Ke=∭VeBTDBdV;Kc=∭VeBTMN-dV;ΔRfe=∬SeNTΔFedS,Δδe=Δδ1TΔδ2T⋯Δδ8TT;Δpe=Δp1Δp2⋯Δp8T;ΔFe=ΔFx,ΔFy,ΔFzT,where [Ke], [Kc], and ΔRf are the stiffness matrix, the coupling matrix, and the equivalent incremental nodal load vector, respectively; Δδ, Δp, and ΔF are incremental vectors of the displacement, EPP, and force, respectively; the strain-displacement matrix B=-∂N=Β1B2⋯B8, Bi is a matrix with 6 × 3 order, and Bi=-∂Ni, i=1,2,…,8; ΔFx, ΔFy, and ΔFz are surface incremental forces in the x-, y-, and z-directions, respectively. In addition, Ve and Se 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 are(22)RIqe=∂∂tMT∂f~-MT∂Hk∂TMp~,RBqe=l-vx+m-vy+n-vz-vn.Equations eliminating the residuals are defined by the following expressions, respectively:(23)∭VeNi∂∂tMT∂f~-MT∂Hk∂TMp~dV=0;i=1,2,…,8,(24)∬SeNil-vx+m-vy+n-vz-vndS=0;i=1,2,…,8.First, substituting (18) into (23), we have(25)∭VeNi∂∂tMT∂f~dV=∭VeNiMT∂Nδ˙dV,where δ˙e=∂δe/∂t. Combining the integration method and the residual of the flow velocity boundary, it can be obtained that(26)-∭VeNiMT∂Hk∂TMp~dV=∬SeNivndS+∭Ve∇TNiHk∇N-pedV,where ∇=∂/∂x∂/∂y∂/∂zT. Combining (23), (25), and (26), we have(27)∭VeNiMT∂Nδ˙dV+∭Ve∇TNiHk∇N-pedV=-∬SeNivndS;i=1,2,…,8.Equation (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:(28)KcTδ˙e-Kspe=Rqe,in which(29)Ks=∭VeBsTHkBsdV;Bs=∂∂x∂∂y∂∂zTN1N2⋯N8;Rqe=∬SeN-TvndS,where [Ks] is the flow matrix; Rqe is the nodal equivalent flow vector; and specifically, for the impermeable boundary, Rqe=0.Integrating both sides of (28) with respect to time from tn to tn+1, we have(30)∫tjtj+1KcTδ˙edt-∫tjtj+1Kspedt=∫tjtj+1Rqedtwith the following approximation expression:(31)∫tjtj+1pedt=Δtpne+θΔpe;∫tjtj+1Rqedt=ΔtRqne.The incremental form of continuity equation can be obtained as follows:(32)KcTΔδe-θΔtKsΔpe=ΔRpein which(33)ΔRpe=ΔtRqje+Kspje,where ΔRpe is the incremental nodal equivalent flow vector and Δt is the time increment. Additionally, θ is an integration constant. A value of θ≥0.5 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[Ke] and [Kc] are the same, while the flow matrix [Ks] is different, leading to a different value for ΔRpe. Matrix evaluations of [Ks] and ΔRpe are the same as it is for general consolidation theory but with an extra process for calculating the permeability matrix [k] and the flow control matrix [H].According to (10) and (16), the permeability matrix at various time steps ki can be evaluated by iterative computation. First, the initial permeability matrix k0 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 p0′ can be calculated. Substituting matrix k0 into the FE iteration, the stress field of this step is obtained, followed by ascertaining p1′ using a similar approach. Substituting p1′ into (16), the flow matrix k1 is obtained. The rest can be deduced by analogy.The evaluation of matrix[H], 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 iΨ and il. It will take on different forms under different gradient conditions. As referenced in the previous section, matrix [H] can be evaluated approximately by using the p value of the previous time step. Stated in detail, under the local coordinate (see Figure 3), the hydraulic gradient of the next time step, ij+1 (the superscript j denotes the time step), can be calculated by the nodal EPP at current time step pj. That is,(34)ixj+1=1γw∂pj∂x=p1j+p4j+p5j+p8j-p2j+p3j+p6j+p7j8γwa,iyj+1=1γw∂pj∂y=p1j+p2j+p5j+p6j-p3j+p4j+p7j+p8j8γwb,izj+1=1γw∂pj∂z=p1j+p2j+p3j+p4j-p5j+p6j+p7j+p8j8γwc,where a, b, and c are half sizes of the element in x-, y-, and z-directions, respectively; see Figure 3. The ascertained values of iΨ, HΨ can be calculated by (13) after a comparison between iΨ and il is completed. Then the flow control matrix [H] of the next time step is obtained.Note here that, for Darcy’s flow, which is a special case of non-Darcian flow, matrix[H] degenerates into a unit matrix. Additionally, two approaches can be used to ascertain the initial flow control matrix H0. 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 H0 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 H0 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 [H] of the current step. The flow control matrix Hi 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 De=3 m; the smear zone diameter is ds=0.2 m; and the drain diameter is dw=0.1 m. According to the principle of equation area, the cylindrical zone is translated into a block zone [25]; see Figure 4. The corresponding sizes Be, bs, and bw for the equivalent system are calculated by a simple function; for example, Be=0.886De. Foundation depth or drain length H=5 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 E and μ, and μ=0.0 is set for 1D deformation analysis. A surcharge loading of q0 is also simulated by instantaneous application to the upper boundary, and the permeability parameters kw, ks, kh, and kv 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.Table 1
Calculation conditions and results for the first verification case.
Cases
k
h (m/s)
k
h/kv
k
h/ks
k
w/kh
E (MPa)
μ
m
i
l
α
Δ
Max
⁡
1 (%)
Condition 1
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.5
20
0.0
9.56
Condition 2
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.5
10
0.0
8.08
Condition 3
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.5
5
0.0
4.13
Condition 4
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.2
5
0.0
2.62
Δ
Max
⁡
1: the largest difference between the FE results and the Hansbo solution [11] for non-Darcy flow under different calculation conditions.Figure 4
Calculation example for the first verification case analyses.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,(35)U-h=1-exp⁡-8chtμ-De2,For non-Darcian flow,(36)U-h=1-1+λtα-De2Δu-0Deγwm-11/1-min which t is consolidation time; for instantaneous applied load, the initial EPP is Δu-0=q0. Additionally,(37)μ-=μ1+μ2+μ3;μ1=De2De2-dw2ln⁡Deds+khksln⁡dsdw-34;μ2=πkhqwz2H-z1-dw2De2;μ3=ds2De2-dw21-ds24De2+khksdw2De2-dw2ds4-dw44De2dw2-ds2dw2+1;ch=khmvγw;mv=1+μ1-2μE1-μ;qw=kwπdw24;κhκs=khks;λ=κhγwmv;α-=m2mβm4m-1m+1;β=β1+β2-β3+β4;β1=13m-1-m-1m3m-15m-1-m-122m25m-17m-1;β2=12mκhκs-1Deds1/m-1-κhκsDedw1/m-1;β3=12m-13m-1κhκs-1Deds1/m-3-κhκsDedw1/m-3;β4=κh2qwπz2l-z1-1mDedw1/m-11-dw2De21/m.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 Figure5(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 H0=I and Hi are directly ascertained using previous time step values without iteration.Figure 5
Comparison of FE results and analytical solutions under different conditions.
(a)
Condition 1
(b)
Condition 2
(c)
Condition 3
(d)
Condition 4From the plots in Figure5, 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 il. The cause of the difference is stated as follows: with the increase of il, 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 ofH (H=5 m). The top boundary is free and drained, while the bottom boundary is fixed and impervious. A surcharge q0 (q0 = 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]:(38)cvilm-11γw∂p∂zm-1∂2p∂z2=∂p∂twheni<il;cv∂2p∂z2=∂p∂twheni≥il.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 Table2. 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.Table 2
Calculation conditions and parameter values and calculation results for Case 2.
Conditions
m
i
l
k
v (m/s)
E (MPa)
μ
Δ
Max
⁡
2 (%)
Condition 1
1.5
20
1
E
-
8
3.0
0.35
2.96
Condition 2
1.2
5.0
1
E
-
8
3.0
0.35
2.47
Δ
Max
⁡
2: the largest difference between the results by FEM and FDM under different conditions.Figure 6
Comparative analysis of FEM and FDM for 1D consolidation with non-Darcy flow model.


## 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 De=3 m; the smear zone diameter is ds=0.2 m; and the drain diameter is dw=0.1 m. According to the principle of equation area, the cylindrical zone is translated into a block zone [25]; see Figure 4. The corresponding sizes Be, bs, and bw for the equivalent system are calculated by a simple function; for example, Be=0.886De. Foundation depth or drain length H=5 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 E and μ, and μ=0.0 is set for 1D deformation analysis. A surcharge loading of q0 is also simulated by instantaneous application to the upper boundary, and the permeability parameters kw, ks, kh, and kv 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.Table 1
Calculation conditions and results for the first verification case.
Cases
k
h (m/s)
k
h/kv
k
h/ks
k
w/kh
E (MPa)
μ
m
i
l
α
Δ
Max
⁡
1 (%)
Condition 1
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.5
20
0.0
9.56
Condition 2
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.5
10
0.0
8.08
Condition 3
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.5
5
0.0
4.13
Condition 4
1
E
-
8
4.0
4.0
1
E
4
10.0
0.0
1.2
5
0.0
2.62
Δ
Max
⁡
1: the largest difference between the FE results and the Hansbo solution [11] for non-Darcy flow under different calculation conditions.Figure 4
Calculation example for the first verification case analyses.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,(35)U-h=1-exp⁡-8chtμ-De2,For non-Darcian flow,(36)U-h=1-1+λtα-De2Δu-0Deγwm-11/1-min which t is consolidation time; for instantaneous applied load, the initial EPP is Δu-0=q0. Additionally,(37)μ-=μ1+μ2+μ3;μ1=De2De2-dw2ln⁡Deds+khksln⁡dsdw-34;μ2=πkhqwz2H-z1-dw2De2;μ3=ds2De2-dw21-ds24De2+khksdw2De2-dw2ds4-dw44De2dw2-ds2dw2+1;ch=khmvγw;mv=1+μ1-2μE1-μ;qw=kwπdw24;κhκs=khks;λ=κhγwmv;α-=m2mβm4m-1m+1;β=β1+β2-β3+β4;β1=13m-1-m-1m3m-15m-1-m-122m25m-17m-1;β2=12mκhκs-1Deds1/m-1-κhκsDedw1/m-1;β3=12m-13m-1κhκs-1Deds1/m-3-κhκsDedw1/m-3;β4=κh2qwπz2l-z1-1mDedw1/m-11-dw2De21/m.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 Figure5(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 H0=I and Hi are directly ascertained using previous time step values without iteration.Figure 5
Comparison of FE results and analytical solutions under different conditions.
(a)
Condition 1
(b)
Condition 2
(c)
Condition 3
(d)
Condition 4From the plots in Figure5, 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 il. The cause of the difference is stated as follows: with the increase of il, 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 ofH (H=5 m). The top boundary is free and drained, while the bottom boundary is fixed and impervious. A surcharge q0 (q0 = 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]:(38)cvilm-11γw∂p∂zm-1∂2p∂z2=∂p∂twheni<il;cv∂2p∂z2=∂p∂twheni≥il.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 Table2. 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.Table 2
Calculation conditions and parameter values and calculation results for Case 2.
Conditions
m
i
l
k
v (m/s)
E (MPa)
μ
Δ
Max
⁡
2 (%)
Condition 1
1.5
20
1
E
-
8
3.0
0.35
2.96
Condition 2
1.2
5.0
1
E
-
8
3.0
0.35
2.47
Δ
Max
⁡
2: the largest difference between the results by FEM and FDM under different conditions.Figure 6
Comparative analysis of FEM and FDM for 1D consolidation with non-Darcy flow model.


## 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 Figure4. The sizes of Be, bs, and bw 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 kh0, kv0, ks0, and kw0 are listed in Table 3, in which parameters corresponding to the initial stress fields γ′ and ϕ′ are also presented.Table 3
Parameters and their values.
B(m)
H (m)
d
s (cm)
d
w (cm)
E  (MPa)
μ (/)
γ
′ (kN/m3)
ϕ
′ (°)
k
h
0 (m/s)
k
h
0/kv0 (/)
k
h
0/ks0 (/)
k
w
0/kh0 (/)
1.15
10.0
15.0
7.0
3.0
0.35
10
30
10−9
1.0
4.0
104As shown in Table4, nine calculation conditions will be conducted in which m is set as 1.5 because this value is validated by several studies, whereas il and a 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 q0 is related to the hydraulic gradient. The higher the value of q0, 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.Table 4
Calculation conditions and input parameters.
Conditions
m
i
l
α
=
C
c/Ck
q
0 (kPa)
Condition 1
1.0
0
0.0
100
Condition 2
1.5
5
0.0
100
Condition 3
1.5
5
0.5
100
Condition 4
1.5
5
1.0
100
Condition 5
1.5
5
0.5
50
Condition 6
1.5
5
0.5
10
Condition 7
1.5
10
0.5
100
Condition 8
1.5
20
0.5
100
Condition 9
1.5
20
1.0
100


### 6.2. Parametric Studies

The consolidation data are plotted in terms of average degree of consolidationU (by (38)) with time t 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.Table 5
Time (days) taken to achieve specific degrees of consolidation under different calculation conditions.
Conditions
C1
C2
C3
C4
C5
C6
C7
C8
C9
U
50
22
23
29
39
30
35
31
33
43
U
90
92
102
131
190
133
167
148
176
259
U
50 and U90: the time taken to achieve U=50% and U=90%, respectively.Table 6
The largest difference of consolidation results for different calculation conditions.
Conditions
C1
C2
C3
C4
C5
C6
C7
C8
C9
Δ
Max
⁡
3 (%)
0.0
2.05
6.95
14.58
7.08
13.50
9.64
13.79
22.16
Δ
Max
⁡
3: the largest difference between the FE results for nonlinear flow under different conditions with linear flow model.Figure 7
Result comparisons under different calculation conditions.
(a)
Assess the influence of nonlinear flow parameterα
(b)
Assess the influence of surcharge intensityq0
(c)
Assess the influence of non-Darcy flow parameterilThe influence of the nonlinear flow parameter (α=Cc/Ck), 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 k 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.0, 0.5, and 1.0, consolidation times U50 (U90) 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 α=1.0.Conditions 3, 5, and 6 are used to analyze the effect of loading size by fixingα and varying q0. The comparisons are shown in Figure 7(b) and Table 5. Conclusions can be obtained that (1) with q0 reducing the consolidation rate is retarded and the consolidation time is prolonged; (2) the influence of q0 on consolidation is insignificant when q0 changes from 100 kPa to 50 kPa, whereas sharp retardation happens when q0 further reduces to 10 kPa. It is to be expected that the fluid flow follows non-Darcy flow when q0 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 (q0 = 10 kPa) and condition 1 is 13.5%.Figure7(c) shows the comparison of consolidation results with different non-Darcy flow parameters. For conditions 3, 7, and 8, il varying from 5 to 20, consolidation times U50 (U90) 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 m and/or il 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 (m=1.5; il=20; α=1.0) 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 il. U50, 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%.


## 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 Figure4. The sizes of Be, bs, and bw 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 kh0, kv0, ks0, and kw0 are listed in Table 3, in which parameters corresponding to the initial stress fields γ′ and ϕ′ are also presented.Table 3
Parameters and their values.
B(m)
H (m)
d
s (cm)
d
w (cm)
E  (MPa)
μ (/)
γ
′ (kN/m3)
ϕ
′ (°)
k
h
0 (m/s)
k
h
0/kv0 (/)
k
h
0/ks0 (/)
k
w
0/kh0 (/)
1.15
10.0
15.0
7.0
3.0
0.35
10
30
10−9
1.0
4.0
104As shown in Table4, nine calculation conditions will be conducted in which m is set as 1.5 because this value is validated by several studies, whereas il and a 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 q0 is related to the hydraulic gradient. The higher the value of q0, 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.Table 4
Calculation conditions and input parameters.
Conditions
m
i
l
α
=
C
c/Ck
q
0 (kPa)
Condition 1
1.0
0
0.0
100
Condition 2
1.5
5
0.0
100
Condition 3
1.5
5
0.5
100
Condition 4
1.5
5
1.0
100
Condition 5
1.5
5
0.5
50
Condition 6
1.5
5
0.5
10
Condition 7
1.5
10
0.5
100
Condition 8
1.5
20
0.5
100
Condition 9
1.5
20
1.0
100


## 6.2. Parametric Studies

The consolidation data are plotted in terms of average degree of consolidationU (by (38)) with time t 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.Table 5
Time (days) taken to achieve specific degrees of consolidation under different calculation conditions.
Conditions
C1
C2
C3
C4
C5
C6
C7
C8
C9
U
50
22
23
29
39
30
35
31
33
43
U
90
92
102
131
190
133
167
148
176
259
U
50 and U90: the time taken to achieve U=50% and U=90%, respectively.Table 6
The largest difference of consolidation results for different calculation conditions.
Conditions
C1
C2
C3
C4
C5
C6
C7
C8
C9
Δ
Max
⁡
3 (%)
0.0
2.05
6.95
14.58
7.08
13.50
9.64
13.79
22.16
Δ
Max
⁡
3: the largest difference between the FE results for nonlinear flow under different conditions with linear flow model.Figure 7
Result comparisons under different calculation conditions.
(a)
Assess the influence of nonlinear flow parameterα
(b)
Assess the influence of surcharge intensityq0
(c)
Assess the influence of non-Darcy flow parameterilThe influence of the nonlinear flow parameter (α=Cc/Ck), 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 k 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.0, 0.5, and 1.0, consolidation times U50 (U90) 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 α=1.0.Conditions 3, 5, and 6 are used to analyze the effect of loading size by fixingα and varying q0. The comparisons are shown in Figure 7(b) and Table 5. Conclusions can be obtained that (1) with q0 reducing the consolidation rate is retarded and the consolidation time is prolonged; (2) the influence of q0 on consolidation is insignificant when q0 changes from 100 kPa to 50 kPa, whereas sharp retardation happens when q0 further reduces to 10 kPa. It is to be expected that the fluid flow follows non-Darcy flow when q0 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 (q0 = 10 kPa) and condition 1 is 13.5%.Figure7(c) shows the comparison of consolidation results with different non-Darcy flow parameters. For conditions 3, 7, and 8, il varying from 5 to 20, consolidation times U50 (U90) 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 m and/or il 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 (m=1.5; il=20; α=1.0) 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 il. U50, 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α (=Cc/Ck, 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 (il and m) 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.



---

*Source: 3047213-2016-02-10.xml*