# Numerical Solution of Nonlinear Large Strain Consolidation Based on Non-Darcian Flow

Authors: Xu-dong Zhao; Wen-hui Gong
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this paper, based on non-Darcian flow, the governing equation of 1D nonlinear large strain consolidation is established, which comprehensively accounts for vertical strain, soil self-weight, geometrical nonlinearity, continuity of pore water flow, the relative velocity of the fluid and solid phases, and changing compressibility and hydraulic conductivity during consolidation. Then the numerical solution is obtained with the finite difference method (FDM). Verification of the FDM solution shows excellent accuracy. On this basis, we investigate the influence of the non-Darcian flow on consolidation behavior. The results show that the increase of the non-Darcian exponent will accelerate the consolidation rate in the beginning, while slowing down the consolidation rate in the end. However, it has no effect on the final settlement of the soil layer. In addition, boundary drainage conditions have a huge impact on the consolidation rate, whether it is Darcian or non-Darcian flow.


---

## Body

## 1. Introduction

The deformation of the soil layer has always been one of the most important fields in geotechnical engineering and civil engineering, which is closely related to the theory of large strain consolidation. If the settlement of the soil layer is not predicted accurately, a large amount of manpower and material resources will be wasted, especially in terms of reclamation projects. Therefore, it is important to establish the theory of 1D large strain consolidation that can reflect the settlement of the soil layer more precisely.So far, in the field of large strain consolidation, plenty of great progress has been made in terms of analytical theory, numerical simulations, experimental verifications, and other aspects [1–8]. Most of these methods, however, were based on the assumption that the water flow in the soils obeyed Darcy’s law. In fact, Darcy’s law was obtained from the indoor sand experiment in the 1850s, which was initially established for coarse-grained soils [9]. The inaccuracy of Darcy’s law in fine-grained soils (e.g., soft clay and dredged fill) has been confirmed by many studies [10–15]. As fine-grained soils are widespread in the civil engineering, it is of theoretical and practical significance to study one-dimensional large strain consolidation based on non-Darcian flow.Slepicka [16] presented a non-Darcian flow law as follows:(1)v=kimin which v signifies the seepage velocity,k represents the vertical hydraulic conductivity,i represents the hydraulic gradient, and m is the non-Darcian flow exponent. Based on (1), many studies have incorporated this flow law into the theories of consolidation [17–21]. All studies above, however, accounted for small strain consolidation rather than large strain consolidation. Zhao et al. [22] first proposed a piecewise linear model, referred to as CNDF1, which has been proven to be both reliable and robust in the analysis of large strain consolidation with the non-Darcian flow.In this paper, we aim to propose a numerical solution for 1D large strain consolidation with nonlinear flow law. The governing equation comprehensively accounts for the large strain, geometrical nonlinearity, self-weight of the soil, the relative velocity of the fluid and solid phases, nonlinear compression, and permeability relationship, and the numerical solution is verified with CNDF1 model. Although the method presented in this paper is totally different from that in the previous study [22], the results of the comparison are extremely close. Furthermore, the effects of the exponential flow on consolidation behavior are further analyzed.


## 2. Mathematical Model

### 2.1. Coordinate Systems

Figure1 illustrates a saturated soil stratum of initial height a0 with the bottom fixed and referenced to both the Lagrangian and the convective coordinate systems. The Lagrangian coordinate a and the convective coordinate ξ are both defined as positive downward.Figure 1
Lagrangian and convective coordinates: (a) initial configuration att =0; (b) configuration at timet.The convective coordinateξ has its origin at the initial top surface of the stratum, which is referred to as the datum plane.The Lagrangian coordinatea and time t are independent variables while the convective coordinate ξ(a,t) is a variable depending on a and t. The relationship between a and ξ is provided by Gibson [1]:(2)∂ξ∂a=1+e1+eoin which e=e(a,t) represents the void ratio and eo represents the initial void ratio.


### 2.2. Basic Assumptions

We adopted the same basic assumptions in the theory of one-dimensional large strain consolidation [1–5, 7, 8], except for water flow law:( 1 ) The soil layer is totally saturated and homogeneous.( 2 ) The pore water and soil particles are incompressible.( 3 ) The settlement of the soil layer is considered to be vertical only.( 4 ) The classical nonlinear constitutive relationships of effective stress, σ′, and vertical hydraulic conductivity, k, are adopted here:(3)e-e1=Cclg⁡σ′1σ′(4)e-e1=Cklg⁡kk1where Cc is the compressibility index and Ck is the permeability index, which are both constant. σ1′ is a known effective stress and k1 is a known hydraulic conductivity, which are user-specified and both constant.( 5 ) The external stress is applied on the layer instantaneously whent=0, which is assumed to be constant during the process of consolidation.


### 2.3. Governing Equation of Large Strain Consolidation with Non-Darcian Flow

Several basic equations for 1D large strain consolidation in the convective coordinate system are present here [1, 2]:( 1 ) Equation for vertical equilibrium of the soil mass:(5)∂σ∂ξ=Gs+e1+eγwin whichσ is the total stress, Gs is the specific gravity of the solid particles (constant), and γw is the unit weight of water (constant).( 2 ) Equation for equilibrium of the pore water:(6)∂u∂ξ=∂uex∂ξ+γwin which u is the total pore pressure and uex is the excess pore pressure.( 3 ) Principle of effective stress:(7)σ=σ′+uin which σ′ is the vertical effective stress.( 4 ) Expression of the hydraulic gradient:(8)i=1γw∂uex∂ξ( 5 ) Equation for continuity of pore water flow:(9)∂∂ξe1+evw-vs=-11+e∂e∂tin which vw and vs are the velocities of the pore water and soil particle relative to datum plane.When considering the relative motion of solid and liquid phases, the flow law as illustrated in (1) can be further described as follows:(10)e1+evw-vs=-kimTherefore, the governing equation can be developed with the combination of all above relations:(11)∂e∂t=1+e0∂∂akGs-11+e-dσ′de∂e∂a1+e01+e1γwmwhere k and σ′ both have relation to void ratio as illustrated in (3) and (4).


### 2.4. Boundary Conditions

The saturated and homogeneous soil layer is taken to have an initially uniform density, and the initial void ratio can be considered as constant throughout the whole layer [23]:(12)ea,0=eoDrainage conditions in the vertical direction are either pervious top and pervious bottom (PTPB) or pervious top and impervious bottom (PTIB), as illustrated in Figure2.Figure 2
The consolidation problem studied.The top boundary is generally pervious, where the excess pore pressure is zero. Therefore, the boundary condition should be as follows:(13)e0,t=e1+Cclg⁡σ′1σ′0,t(14)σ′0,t=q+σ0′where q is the instantaneous external stress applied on the surface of the soil stratum at the beginning of consolidation, which is assumed to be constant. σ0′ is the initial effective stress:(15)σ0′=σ1′exp⁡e1-eoCcSimilarly, for the pervious bottom boundary, the boundary condition should be(16)ea0,t=e1+Cclg⁡σ′1σ′a0,tSince the settlement of the fix bottom boundary is always zero, it can be assumed that the bottom boundary is always in the final state of consolidation [23]. Therefore,(17)σ′a0,t=σ′a0,∞=γs-γwa01+e0+q+σ0′where γs=Gs•γw is the unit weight of solid particles (constant).The bottom boundary can also be impervious, where the hydraulic gradient is zero. Therefore, for the impervious bottom boundary, the boundary condition should be(18)∂e∂aa=a0=-γs-γw1+e0avwhere av=-de/dσ′ is the coefficient of compressibility.


## 2.1. Coordinate Systems

Figure1 illustrates a saturated soil stratum of initial height a0 with the bottom fixed and referenced to both the Lagrangian and the convective coordinate systems. The Lagrangian coordinate a and the convective coordinate ξ are both defined as positive downward.Figure 1
Lagrangian and convective coordinates: (a) initial configuration att =0; (b) configuration at timet.The convective coordinateξ has its origin at the initial top surface of the stratum, which is referred to as the datum plane.The Lagrangian coordinatea and time t are independent variables while the convective coordinate ξ(a,t) is a variable depending on a and t. The relationship between a and ξ is provided by Gibson [1]:(2)∂ξ∂a=1+e1+eoin which e=e(a,t) represents the void ratio and eo represents the initial void ratio.


## 2.2. Basic Assumptions

We adopted the same basic assumptions in the theory of one-dimensional large strain consolidation [1–5, 7, 8], except for water flow law:( 1 ) The soil layer is totally saturated and homogeneous.( 2 ) The pore water and soil particles are incompressible.( 3 ) The settlement of the soil layer is considered to be vertical only.( 4 ) The classical nonlinear constitutive relationships of effective stress, σ′, and vertical hydraulic conductivity, k, are adopted here:(3)e-e1=Cclg⁡σ′1σ′(4)e-e1=Cklg⁡kk1where Cc is the compressibility index and Ck is the permeability index, which are both constant. σ1′ is a known effective stress and k1 is a known hydraulic conductivity, which are user-specified and both constant.( 5 ) The external stress is applied on the layer instantaneously whent=0, which is assumed to be constant during the process of consolidation.


## 2.3. Governing Equation of Large Strain Consolidation with Non-Darcian Flow

Several basic equations for 1D large strain consolidation in the convective coordinate system are present here [1, 2]:( 1 ) Equation for vertical equilibrium of the soil mass:(5)∂σ∂ξ=Gs+e1+eγwin whichσ is the total stress, Gs is the specific gravity of the solid particles (constant), and γw is the unit weight of water (constant).( 2 ) Equation for equilibrium of the pore water:(6)∂u∂ξ=∂uex∂ξ+γwin which u is the total pore pressure and uex is the excess pore pressure.( 3 ) Principle of effective stress:(7)σ=σ′+uin which σ′ is the vertical effective stress.( 4 ) Expression of the hydraulic gradient:(8)i=1γw∂uex∂ξ( 5 ) Equation for continuity of pore water flow:(9)∂∂ξe1+evw-vs=-11+e∂e∂tin which vw and vs are the velocities of the pore water and soil particle relative to datum plane.When considering the relative motion of solid and liquid phases, the flow law as illustrated in (1) can be further described as follows:(10)e1+evw-vs=-kimTherefore, the governing equation can be developed with the combination of all above relations:(11)∂e∂t=1+e0∂∂akGs-11+e-dσ′de∂e∂a1+e01+e1γwmwhere k and σ′ both have relation to void ratio as illustrated in (3) and (4).


## 2.4. Boundary Conditions

The saturated and homogeneous soil layer is taken to have an initially uniform density, and the initial void ratio can be considered as constant throughout the whole layer [23]:(12)ea,0=eoDrainage conditions in the vertical direction are either pervious top and pervious bottom (PTPB) or pervious top and impervious bottom (PTIB), as illustrated in Figure2.Figure 2
The consolidation problem studied.The top boundary is generally pervious, where the excess pore pressure is zero. Therefore, the boundary condition should be as follows:(13)e0,t=e1+Cclg⁡σ′1σ′0,t(14)σ′0,t=q+σ0′where q is the instantaneous external stress applied on the surface of the soil stratum at the beginning of consolidation, which is assumed to be constant. σ0′ is the initial effective stress:(15)σ0′=σ1′exp⁡e1-eoCcSimilarly, for the pervious bottom boundary, the boundary condition should be(16)ea0,t=e1+Cclg⁡σ′1σ′a0,tSince the settlement of the fix bottom boundary is always zero, it can be assumed that the bottom boundary is always in the final state of consolidation [23]. Therefore,(17)σ′a0,t=σ′a0,∞=γs-γwa01+e0+q+σ0′where γs=Gs•γw is the unit weight of solid particles (constant).The bottom boundary can also be impervious, where the hydraulic gradient is zero. Therefore, for the impervious bottom boundary, the boundary condition should be(18)∂e∂aa=a0=-γs-γw1+e0avwhere av=-de/dσ′ is the coefficient of compressibility.


## 3. Derivation of Numerical Solution

### 3.1. Alternative Form of the Governing Equation

For the convenience of following studies, an equivalent Darcy hydraulic conductivity,Kd=kim-1, is adopted in this paper. In this way, (1) can be expressed in a form similar to Darcy’s law, namely, v=Kdi. The expression of Kd is as follows:(19)Kd=kGs-11+e-dσ′de∂e∂a1+e01+e1γwm-1Therefore, (11) can be further expressed as(20)∂e∂t=Gs-1ddeKd1+e01+e∂e∂a-1+e02γw∂∂aKd1+edσ′de∂e∂aor in the alternative form:(21)∂e∂t=A∂e∂a+B∂e∂a2+C∂2e∂a2where(22)A=GS-1ddeKd1+e01+e(23)B=-1+e02γwdσ′de•ddeKd1+e+Kd1+ed2σ′de2(24)C=-1+e02Kdγw1+ed2σ′de2The derivation process of (21) is specifically illustrated in the appendix.


### 3.2. Discretization of the Governing Equation

As stated above, the finite difference method (FDM) is introduced to develop the numerical solution for (21). Firstly, the spatial domain 0≤a≤a0 is subdivided into M equal elements from the top down by the proportional spacing Δa. If the nodal point of spatial domain is denoted as ai(25)ai=iΔaMeanwhile, if the nodal point of time is denoted astj,(26)tj=jΔtThe Crank-Nicolson difference scheme is used in this paper for its relative stability. Therefore, the difference scheme of (21) can be obtained:(27)eij-eij-1Δt=Aij-1ei+1j-1-ei-1j-12Δa+Bij-1ej-1i+1-ei-1j-12Δa2+12Cij-1×ei+1j+ei-1j-2eijΔa2+ei+1j-1+ei-1j-1-2eij-1Δa2,i=1,2,…,Mwhere eij denotes the void ratio at ai=iΔa and tj=jΔt. Aij-1, Bij-1, and Cij-1 are the corresponding value at ai=iΔa and tj-1=(j-1)t. For the sake of calculation, (27) can be rewritten in the alternative form as follows:(28)ei+1j-21+Δa2ΔtCij-1eij+ei-1j=-ei+1j-1+21-Δa2ΔtCij-1eij-1-ei-1j-1-Aij-1Cij-1Δaei+1j-1-ei-1j-1-Bij-12Cij-1ei+1j-1-ei-1j-12,i=1,2,…,MThe initial condition can be expressed in terms of discrete nodal points:(29)ei0=eo,i=1,2,…,Mand the boundary conditions can also be expressed as(30)PTIB:e0j=e1+Cclg⁡σ′1q+σ1′exp⁡e1-e0/CceM+1j=eM-1j-2Δaγs-γw1+e0avj(31)PTPB:e0j=e1+Cclg⁡σ′1q+σ1′exp⁡e1-e0/CceM+1j=e1+Cclg⁡σ′1γs-γwa0/1+e0+q+σ1′exp⁡e1-e0/Ccwhere e0j denotes the void ratio at the top boundary and eM+1j denotes the void ratio at the bottom boundary.For PTIB, by substituting (30) into (28), equation (28) can be expressed in terms of matrix, namely, [D]{e}=[E]:(32)D111D211D31⋱⋱⋱1DM-112DMj-1e1e2e3⋮eM-1eMj=E1E2E3⋮EM-1EMj-1where(33)Di=-21+Δa2ΔtCij-1(34)E1=-2e0j-1+21-Δa2ΔtC1j-1e1j-1-e2j-1-ΔaA1j-1C1j-1e2j-1-e0j-1-B1j-12C1j-1e2j-1-e0j-12(35)Ei=-ei+1j-1+21-Δa2ΔtCij-1eij-1-ei-1j-1-ΔaAij-1Cij-1ei+1j-1-ei-1j-1-Bij-12Cij-1ei+1j-1-ei-1j-12(36)EM=21-Δa2ΔtCMj-1eMj-1-2eM-1j-1-2ΔaAMj-1CMj-1ΓMj-1-2Δa2BMj-1CMj-1ΓMj-12+4ΔaΓMj-1where(37)ΓM=∂e∂aa=a0=-γs-γw1+e0avSimilarly, for PTPB, (28) can be expressed as follows:(38)D111D211D31⋱⋱⋱1DM-111DMj-1e1e2e3⋮eM-1eMj=E1E2E3⋮EM-1Fj-1where(39)F=-eM-1j-1+21-Δa2ΔtCMj-1eMj-1-2eM+1j-1-ΔaAMj-1CMj-1eM+1j-1-eM-1j-1-BMj-12CMj-1eM+1j-1-eM-1j-12and the other symbols are the same as those in (32).It can be seen that coefficient matrices [D] are known when t=(j-1)Δt, which are both symmetrical three-diagonal. Therefore, the value of eij at any moment tj=jΔt can be obtained once we know the value of the void ratio at the previous moment tj-1=(j-1)Δt. With Matlab programming, it is easy to obtain the distribution of the void ratio during the process of consolidation, as illustrated in Figure 3.Figure 3
Flow chart for computer procedure.


### 3.3. Excess Pore Pressure, Discharge Velocity, Settlement, and Average Degree of Consolidation

Once the value ofeij is obtained, the value of excess pore stress can also be determined according to (6) and (7):(40)uex,ij=γs-γwaij1+e0+q+σ0-σ′ijwhere the effective stress σ′ij can be calculated from (3):(41)σ′ij=σ1′expe1-eijCcSimilarly,kij can be calculated from (4):(42)kij=k1exp⁡eij-e1CkTherefore, the velocity of water flow in the soil layer can be obtained according to (10):(43)vij=kijGs-11+eij-dσ′deej-1i+1-ei-1j-12Δa1+e01+e1γwmi=1,2,…,MAccording to (2), the settlement of the layer can be obtained:(44)St=a0-H=a0-11+e0a0+∑i=1M+1eij+ei-1j2Δaand the average degree of consolidation, Uavg, is as follows:(45)Uavg=StS∞where S∞ is the ultimate settlement:(46)S∞=a0-11+e0a0+∑i=1M+1ei∞+ei-1∞2Δa


## 3.1. Alternative Form of the Governing Equation

For the convenience of following studies, an equivalent Darcy hydraulic conductivity,Kd=kim-1, is adopted in this paper. In this way, (1) can be expressed in a form similar to Darcy’s law, namely, v=Kdi. The expression of Kd is as follows:(19)Kd=kGs-11+e-dσ′de∂e∂a1+e01+e1γwm-1Therefore, (11) can be further expressed as(20)∂e∂t=Gs-1ddeKd1+e01+e∂e∂a-1+e02γw∂∂aKd1+edσ′de∂e∂aor in the alternative form:(21)∂e∂t=A∂e∂a+B∂e∂a2+C∂2e∂a2where(22)A=GS-1ddeKd1+e01+e(23)B=-1+e02γwdσ′de•ddeKd1+e+Kd1+ed2σ′de2(24)C=-1+e02Kdγw1+ed2σ′de2The derivation process of (21) is specifically illustrated in the appendix.


## 3.2. Discretization of the Governing Equation

As stated above, the finite difference method (FDM) is introduced to develop the numerical solution for (21). Firstly, the spatial domain 0≤a≤a0 is subdivided into M equal elements from the top down by the proportional spacing Δa. If the nodal point of spatial domain is denoted as ai(25)ai=iΔaMeanwhile, if the nodal point of time is denoted astj,(26)tj=jΔtThe Crank-Nicolson difference scheme is used in this paper for its relative stability. Therefore, the difference scheme of (21) can be obtained:(27)eij-eij-1Δt=Aij-1ei+1j-1-ei-1j-12Δa+Bij-1ej-1i+1-ei-1j-12Δa2+12Cij-1×ei+1j+ei-1j-2eijΔa2+ei+1j-1+ei-1j-1-2eij-1Δa2,i=1,2,…,Mwhere eij denotes the void ratio at ai=iΔa and tj=jΔt. Aij-1, Bij-1, and Cij-1 are the corresponding value at ai=iΔa and tj-1=(j-1)t. For the sake of calculation, (27) can be rewritten in the alternative form as follows:(28)ei+1j-21+Δa2ΔtCij-1eij+ei-1j=-ei+1j-1+21-Δa2ΔtCij-1eij-1-ei-1j-1-Aij-1Cij-1Δaei+1j-1-ei-1j-1-Bij-12Cij-1ei+1j-1-ei-1j-12,i=1,2,…,MThe initial condition can be expressed in terms of discrete nodal points:(29)ei0=eo,i=1,2,…,Mand the boundary conditions can also be expressed as(30)PTIB:e0j=e1+Cclg⁡σ′1q+σ1′exp⁡e1-e0/CceM+1j=eM-1j-2Δaγs-γw1+e0avj(31)PTPB:e0j=e1+Cclg⁡σ′1q+σ1′exp⁡e1-e0/CceM+1j=e1+Cclg⁡σ′1γs-γwa0/1+e0+q+σ1′exp⁡e1-e0/Ccwhere e0j denotes the void ratio at the top boundary and eM+1j denotes the void ratio at the bottom boundary.For PTIB, by substituting (30) into (28), equation (28) can be expressed in terms of matrix, namely, [D]{e}=[E]:(32)D111D211D31⋱⋱⋱1DM-112DMj-1e1e2e3⋮eM-1eMj=E1E2E3⋮EM-1EMj-1where(33)Di=-21+Δa2ΔtCij-1(34)E1=-2e0j-1+21-Δa2ΔtC1j-1e1j-1-e2j-1-ΔaA1j-1C1j-1e2j-1-e0j-1-B1j-12C1j-1e2j-1-e0j-12(35)Ei=-ei+1j-1+21-Δa2ΔtCij-1eij-1-ei-1j-1-ΔaAij-1Cij-1ei+1j-1-ei-1j-1-Bij-12Cij-1ei+1j-1-ei-1j-12(36)EM=21-Δa2ΔtCMj-1eMj-1-2eM-1j-1-2ΔaAMj-1CMj-1ΓMj-1-2Δa2BMj-1CMj-1ΓMj-12+4ΔaΓMj-1where(37)ΓM=∂e∂aa=a0=-γs-γw1+e0avSimilarly, for PTPB, (28) can be expressed as follows:(38)D111D211D31⋱⋱⋱1DM-111DMj-1e1e2e3⋮eM-1eMj=E1E2E3⋮EM-1Fj-1where(39)F=-eM-1j-1+21-Δa2ΔtCMj-1eMj-1-2eM+1j-1-ΔaAMj-1CMj-1eM+1j-1-eM-1j-1-BMj-12CMj-1eM+1j-1-eM-1j-12and the other symbols are the same as those in (32).It can be seen that coefficient matrices [D] are known when t=(j-1)Δt, which are both symmetrical three-diagonal. Therefore, the value of eij at any moment tj=jΔt can be obtained once we know the value of the void ratio at the previous moment tj-1=(j-1)Δt. With Matlab programming, it is easy to obtain the distribution of the void ratio during the process of consolidation, as illustrated in Figure 3.Figure 3
Flow chart for computer procedure.


## 3.3. Excess Pore Pressure, Discharge Velocity, Settlement, and Average Degree of Consolidation

Once the value ofeij is obtained, the value of excess pore stress can also be determined according to (6) and (7):(40)uex,ij=γs-γwaij1+e0+q+σ0-σ′ijwhere the effective stress σ′ij can be calculated from (3):(41)σ′ij=σ1′expe1-eijCcSimilarly,kij can be calculated from (4):(42)kij=k1exp⁡eij-e1CkTherefore, the velocity of water flow in the soil layer can be obtained according to (10):(43)vij=kijGs-11+eij-dσ′deej-1i+1-ei-1j-12Δa1+e01+e1γwmi=1,2,…,MAccording to (2), the settlement of the layer can be obtained:(44)St=a0-H=a0-11+e0a0+∑i=1M+1eij+ei-1j2Δaand the average degree of consolidation, Uavg, is as follows:(45)Uavg=StS∞where S∞ is the ultimate settlement:(46)S∞=a0-11+e0a0+∑i=1M+1ei∞+ei-1∞2Δa


## 4. Verification of Numerical Solutions

The following example is to verify the accuracy of numerical method with CNDF1 model. A soil layer in [24] is investigated with a0=10 m and Gs=2.78. The initial void ratio, e0=1.571, is constant throughout the whole stratum, and an instantaneous loadingq=100kPa is applied on the top surface of the stratum whent=0. Constitutive relationships are taken from the experimental data [24], e1= 1.57, σ′1=50 kPa, k1=1 × 10−8, Cc=0.6378, and Ck=0.8. Taking the non-Darcian exponentm=1.2 for an example, the comparison results of FDM solutions and CNDF1 simulations are as follows.Figure4 illustrates the distribution of void ratio at several different moments, with Figure 4(a) for PTIB and Figure 4(b) for PTPB. In the figure, the dimensionless elevation, (a0-ξ)/H, is adopted to describe the relative position of each element in the layer. Figure 5 shows the uex profiles at several different moments obtained by two methods, with Figure 5(a) for PTIB and Figure 5(b) for PTPB. Figures 6 and 7 show the curves of settlement and the average degree of consolidation, respectively. Figure 8 shows the curves of discharge velocity at pervious boundaries; the signal of negative means that the direction of velocity is against with gravity.Figure 4
Verification for the void ratio: (a) PTIB; (b) PTPB.
(a) (b)Figure 5
Verification for the excess pore pressure: (a) PTIB; (b) PTPB.
(a) (b)Figure 6
Verification for the settlement of layer: (a) PTIB; (b) PTPB.
(a) (b)Figure 7
Verification for the average degree of consolidation: (a) PTIB; (b) PTPB.
(a) (b)Figure 8
Verification for the discharge velocity at the boundaries: (a) PTIB; (b) PTPB.
(a) (b)The results of the comparison show excellent agreement between two methods. It means that FDM solution has good accuracy in the analysis of consolidation behavior with non-Darcian flow.


## 5. Influence of Boundary Conditions

According to Figures4 and 5, the results indicate that the pervious bottom boundary (PTPB) results in smaller values of the void ratio and residual excess pore pressure at the same depth. The reason for this phenomenon is that a greater deformation has occurred. This means that pervious boundary can speed up the process of consolidation. Therefore, the consolidation rate under the double-drained condition is obviously faster than that under the single-drained condition, just as illustrated in Figures 6, 7, and 8.In addition, it is interesting that the value of void ratio at the pervious boundary is constant, while the value at the impervious bottom boundary is changed with time as shown in Figure4(a). The reason can be easily illustrated by (18).


## 6. Influence of Non-Darcian Flow on the Consolidation Behavior

A soil layer with initial height of 8 m is affected by the external loadingq=100kPa whent=0. The initial void ratio is 1.445 throughout the layer and the specific gravity is 2.61. The constitutive relationships are taken from [25], where e1= 1.57, σ′1=50 kPa, k1=1 × 10−8, Cc=0.6378, and Ck=0.8. The effects of non-Darcian flow are investigated with m increasing from 0.75 to 1.5.Figure9 shows the curves of settlement, with Figure 9(a) for PTIB and Figure 9(b) for PTPB. There is an interesting phenomenon that the increase of the non-Darcian exponent will accelerate the consolidation rate in the beginning, while slowing down the consolidation rate in the end. The reason for this phenomenon is that the value of hydraulic gradient will decrease during the process of consolidation, and the discharge velocity will decrease with increasing the value of m when the hydraulic gradient is less than 1. It is noted that the consolidation rate under the double-drained condition is obviously faster than that under the single-drained condition, but the ultimate settlement is the same in both cases, indicating that the drainage conditions at the bottom boundary have no effect on the final settlement.Figure 9
The settlement of layer with different value ofm: (a) PTIB; (b) PTPB.
(a) (b)Figures10 and 11 show the distribution of void ratio and excess pore pressure for PTIB, respectively. In Figures 10(a) and 11(a), a larger value of m results in smaller values of the void ratio residual excess pore pressure at the same depth. The reason for this phenomenon is that a greater deformation has occurred at this time (t=1 year). This means that the increase of the value of m can speed up the process of consolidation in the beginning. However, Figures 10(b) and 11(b) show an opposite trend. The reason for this phenomenon is that at the later stage, the increase of the value of m will slow down the consolidation rate in the end, as illustrated in the previous study [22].Figure 10
The distribution of void ratio for PTIB with different value ofm: (a)t=1 year; (b)t=5 years.
(a) (b)Figure 11
The distribution of excess pore pressure for PTIB with different value ofm: (a)t=1 year; (b)t=5 years.
(a) (b)Figure12 illustrates the distribution of void ratio for PTPB. It can be observed that the non-Darcian flow m also has a similar effect for PTPB. This means that for both PTIB and PTPB, a larger value of m results in smaller values of the void ratio and residual excess pore pressure at the same depth in the beginning of consolidation and larger values of these parameters in the end. Besides, it is noted that the void ratio at the boundaries of stratum is constant during the process of consolidation. The reason is that the void ratio is only related to the external load and the soil self-weight, as illustrated in (17). Figure 13 illustrates the uex profiles for PTPB. It can be seen that the value of uex at the permeable boundaries is always zero, just as mentioned before.Figure 12
The distribution of void ratio for PTPB with different value ofm: (a)t=1 year; (b)t=2 years.
(a) (b)Figure 13
The distribution of excess pore pressure for PTPB with different value ofm:(a)t=1 year; (b)t=2 years.
(a) (b)


## 7. Influence of Soil Self-Weight on the Non-Darcian Consolidation

The third example is to study the influence of soil self-weight on the non-Darcian consolidation. Figures14 and 15 illustrate the settlement and Uavg of the layer, respectively. It can be seen that the results of the final settlement are smaller than actual if soil self-weight is negligible, and there is also a significant difference in the results of the consolidation rate. When the soil self-weight is not considered, the consolidation rate is slower than the actual. However, this phenomenon is different from that in [26]. The reason is that in this paper, we adopt the same initial condition as [23], which is not mentioned in [26]. This means that the initial void ratio is constant throughout the whole layer, but the excess pore pressure is not zero under the initial condition, which results in a totally different conclusion of the effect of soil self-weight on the consolidation behavior.Figure 14
Effects of soil self-weight on the settlement: (a) PTIB;(b) PTPB.
(a) (b)Figure 15
Effects of soil self-weight on the average degree of consolidation: (a) PTIB;(b) PTPB.
(a) (b)


## 8. Conclusions

In this paper, the governing equation of 1D large strain consolidation is first established, with the consideration of non-Darcian flow and geometrical nonlinearity. With finite difference method, numerical solution is obtained. Comparisons of FDM solution and CNDF1 simulation show excellent agreement. The influence of the nonlinear flow law on consolidation behavior is investigated.( 1 ) The FDM solution is in high accuracy, even though the complicated coefficient matrices [D] take a large amount of time to calculate in the computer procedure. Besides, it is difficult to obtain the analytical solution of the governing equation. Therefore, we can adopt the FDM solution instead of the analytical solution to some extent.( 2 ) The increase of the non-Darcian exponent will accelerate the consolidation rate in the beginning, while slowing down the consolidation rate in the end.( 3 ) Although non-Darcian seepage will affect the consolidation rate, it has no effect on the final settlement of the soil layer.( 4 ) Boundary drainage conditions have a huge impact on the consolidation rate, whether it is Darcian or non-Darcian flow. The consolidation rate under the double-drained condition is obviously faster than that under the single-drained condition.( 5 ) The influence of soil self-weight is significant. However, the different initial condition may result in a totally different conclusion of effect of soil self-condition on the consolidation behavior.



---

*Source: 5745068-2019-05-27.xml*