Large-Strain Nonlinear Consolidation of Sand-Drained Foundations Considering Vacuum Preloading and the Variation in Radial Permeability Coefficient

Abstract

The vacuum preloading method effectively strengthens soft soil foundations with vertical drainage, which produces a smear effect when laying sand drains. Meanwhile, the seepage of pore water and soil deformation during consolidation exhibit nonlinear characteristics. Therefore, based on Gibson’s 1D large-strain consolidation theory, this paper developed a more generalized large-strain radical consolidation model of sand-drained soft foundations under free-strain assumptions. In this system, the double logarithmic compression permeability relationships for soft soils with large-strain properties, the variation in the radical permeability coefficient in the smear zone, and the effect of the non-Darcy flow were all included. Then, the partial differential control equations were numerically solved by the finite difference method and validated with existing radical consolidation test results and derived analytical solutions. Finally, the influences of relevant model parameters on consolidation are discussed. The analysis shows that the greater the maximum dimensionless vacuum negative pressure P₀, the faster the consolidation rate of sand-drained foundations. Meanwhile, the decrease in the negative pressure transfer coefficient k₁ will result in a decreasing final settlement amount. Moreover, the consolidation rate of sand-drained foundations is slower considering the non-Darcy flow, but the final settlement is unaffected.

1. Introduction

The vacuum preloading method is one of the most effective methods to treat soft soil ground with the drainage consolidation method. Its principle aim is to produce negative pressure in the soil as a load through vacuum pumping, and then the water in the soil pores is drained so that the ground soil is consolidated and the strength is improved. Meanwhile, combined with vacuum pressure, the load can be quickly applied at once, with fast reinforcement speed and a short construction period. Previous studies [1,2,3,4] have shown that the vacuum preloading method will produce large deformations when dealing with soft soils with large void ratios, and the permeability coefficient decreases with the decrease in the soil–void ratio. Meanwhile, the construction disturbance of the installed sand drain produces a smear effect, so the variation in the radical permeability coefficient within the smear zone should also be considered. Meanwhile, it is also indisputable that water seepage in soft clay deviates from Darcy’s flow [5,6,7]. Therefore, a more generalized large-strain nonlinear consolidation model of soft soil foundations with sand drains is established in this paper, considering the coupled effect of the nonlinear compression relationship and permeability relationships, the variable permeability coefficient in the smear zone, and the non-Darcy flow effect. It is expected to simulate the radial drainage consolidation of deep clayey foundations more truly.

For the study of vacuum preloading on sand-drained foundations, relevant scholars have conducted much research and achieved rich results [8,9,10,11,12]. Among them, studies [8,9] developed Terzaghi’s 1D consolidation system, incorporating vacuum preloading. Chai et al. [10] established a method for estimating the final vacuum pressure distribution and proposed the corresponding model for the vacuum-drain consolidation-induced settlement curves. Indraratna et al. [11] and Geng et al. [12], respectively, deduced an analytical solution for vertical drain consolidation foundations with linear attenuation of vacuum preloading: conventional and point array vacuum preloading. However, these studies have neglected the large-strain consolidation characteristics of soft clayey soils, and research reveals that the small-strain theory would not be suitable for highly compressible soil, such as soft clayey soils with large-stain properties [13,14]. In this regard, McVay et al. [15] investigated large-strain nonlinear consolidation of sand-drained foundations in clay by assuming that the void ratio varies exponentially with effective stress. Perera et al. [16] derived a large-strain consolidation analytical solution for vacuum preloading incorporating the semi-logarithmic nonlinear compression relationship and soil disturbance. Geng et al. [17] and Nguyen et al. [18] respectively studied the nonlinear large-strain consolidation of sand-drained foundations considering self-weight and well resistance decay over time by using a semi-logarithmic relationship (i.e., $e-\lg {\sigma }^{\prime }$ and $e-\lg k$). Walker et al. [19] presented an approximate solution for the consolidation of sand-drained foundations considering smear zone overlap. In more in-depth research, R. Butterfield [20] proposed to apply double logarithmic compression and permeability relationships (i.e., Equations (2) and (3)) to analyze the consolidation problem and pointed out that the double-logarithmic form of relationships can be well-consistent with the test results for soft clayey soils under a large-strain scenario. Tavenas et al. [21] also demonstrated that the semi-logarithmic form of relationships only satisfies linear relationships in low-strain categories below 20%. Based on double logarithmic compression and permeability relationships, Li et al. [22] initially developed the 1D large-strain consolidation for soft ground and derived the corresponding analytical solutions. However, the above studies on sand-drained foundation consolidation often assume that pore water seepage in the soil is consistent with Darcy’s flow.

In terms of numerical methods in large-strain problems, compared with common, familiar numerical methods such as the finite element method and finite volume method, the finite difference method is more popular and widely used. In addition, in recent years, the use of meshless methods for complex nonlinear partial differential problems has gradually gained the favor of scholars [23,24,25,26]. This type of method has significant advantages in dealing with some special problems, such as simulating large deformation problems, and can be one of the key research directions for future numerical methods.

According to the research results of indoor tests [5,6,7], it was found that pore water seepage exhibited nonlinear characteristics under low hydraulic gradients. Among them, considering the non-Newtonian fluid characteristics of soil, Swartzendurber [6] summarized different regions’ minor data and proposed a continuous mathematical expression flow rule with only two flow parameters. This flow relationship could fit well with the test data and easily be applied in practical engineering. Based on this, scholars [27,28] have conducted a 1D consolidation considering this flow relationship and obtained a more generalized understanding of 1D consolidation characteristics.

In summary, although the research on sand-drain consolidation under vacuum preloading has achieved certain results, the influence of the nonlinearity compression and permeability relationships, the radical permeability coefficient variation, and non-Darcy flow still have not been considered simultaneously. Therefore, based on Gibson’s [13] 1D large-strain consolidation theory and Barron’s [29] free-strain assumption, an improved large-strain vacuum preloading consolidation equation of soft ground with sand drain is established and solved by the finite difference method. The double logarithmic compression and permeability relationships for large voids ratio to soft clays, the three different permeability variation modes in the smear zone, and the non-Darcy flow with continuous function form were also included in this system. Based on the proposed vertical drain model solutions, relevant parameters affecting the consolidation of soft foundations with sand drains are further analyzed.

2. Large Radical Strain Consolidation Model with Sand Drains

2.1. Problem Description

As shown in Figure 1, assuming the thickness of sand-drained foundations is H, a sand drain is installed in the foundation for drainage consolidation. Considering the smear effect caused by construction disturbance, the whole influence area of radial seepage is divided into a smear zone and an undisturbed zone. The radius of the sand drain, smear area, and influence area are r_w, r_s, and r_e, respectively. r and z denote the radial and vertical coordinates, respectively, and −p_v(z) is the vacuum negative pressure acting on the sand-drain zone. The edge of the sand drain is a fully pervious boundary, and the outermost periphery of the weak smear area is an impermeable boundary. q(t) is the overload load on the foundation surface (see Figure 2), where t₁ is the completion of the single-loading time, q₀ is the initial load value, and q_u is the final value.

2.2. Basic Assumptions

To facilitate the solution, the following assumptions are made on the basis of Gibson’s 1D large deformation consolidation theory and Barron’s free-strain consolidation theory:

(1)

The soil particles only move vertically, and the radical geometric deformation is ignored.

(2)

Carrillo [30] concluded that the consolidation of sand-drained foundations can be decomposed into radial and vertical seepage for calculation, making the radial consolidation theory more practical. Therefore, in parallel with the existing research [29], only radical seepage is considered in the process of soft foundations with sand-drain consolidation.

(3)

As shown in Figure 3, it is assumed that the seepage of pore water corresponds to the non-Newtonian index flow [6], which is:

$$v=-k_r\left[i-i_0(1-e^{-i/i_0})\right]$$

(1)

where v is the seepage velocity; i is the hydraulic gradient; k_r is the radial permeability coefficient; i₀ is the non-Newtonian index.

The nonlinear relationship of soft soils can be simulated by the double logarithmic compression and permeability relationship proposed by R. Butterfield (1979) [20], i.e.,

$$\lg \left(1+e\right)-\lg (1+e_0)=I_c\left(\lg {\sigma }_0^{\prime }-\lg {\sigma }^{\prime }\right)$$

(2)

$$\lg \left(1+e\right)-\lg \left(1+e_0\right)=\frac{1}{\alpha }\left(\lg k_{\mathrm{r}}-\lg k_{\mathrm{r}0}\right)$$

(3)

where e₀ and e are the initial void ratio and the void ratio at any time, respectively. ${\sigma }_0^{\prime }$ and σ′ are the initial effective stress and the initial effective stress at any time. k_r is the radical permeability coefficient at the radius of the undisturbed area at any time, and k_r0 is the initial radical permeability coefficient. I_c is the slope of the lg(1 + e) − lg σ′ linear relationship, α is the reciprocal of the slope of the linear relationship lg(1 + e) − lg k_r. For convenience, I_c and α are expressed by the compression index and permeability index.

(4)

Suppose the vacuum negative pressure at the sand well boundary changes from −p₀ to −k₁p₀ along the depth, i.e., Equation (4):

$$-p_v(z)=-p_0\left[1-\left(1-k_1\right)\frac{z}{H}\right]$$

(4)

where p₀ is the vacuum preloading value at the foundation’s top surface, k₁ is the negative pressure transfer coefficient, and its value is not greater than 1.

Considering the variation in the initial radial permeability coefficient along the radial direction within the smear zone, i.e., the smear effect due to the construction disturbance. That is, letting k_r(r) = k_rf(r) and assuming three variation modes of radial permeability coefficient, as shown in Figure 4. Mode 1 is the radical permeability coefficient that keeps constant at the smear zone; modes 2 and 3 assume that the radical permeability coefficient increases linearly and parabolically with the increase in radial radius r of the zone, respectively.

Equations (5)–(7) are mathematical functional relations of modes 1, 2, and 3, respectively.

Mode 1:

$$f(r)=\{\begin{array}{cc}\delta ,& r_w\le r\le r_s\\ 1,& r_s<r\le r_e\end{array}$$

(5)

Mode 2:

$$f(r)=\frac{r-r_w}{r_e-r_w}(1-\delta )+\delta ,r_w\le r\le r_e$$

(6)

Mode 3:

$$f(r)=\{\begin{array}{cc}(1-\delta )\left(a_1-b_1+c_1\frac{r}{r_w}\right)\left(a_1+b_1-c_1\frac{r}{r_w}\right),& r_w\le r\le r_s\\ 1,& r_s<r\le r_e\end{array}$$

(7)

where $a_1=1/\sqrt{1-\delta }$, $b_1=s/(s-1)$, $c_1=1/(s-1)$, $s=r_s/r_w$, δ = k_s/k_r.

2.3. Control Equations

According to the above basic assumptions and the continuous flow conditions of radial sand-drain drainage consolidation, there are:

$$\frac{\partial v}{\partial r}+\frac{v}{r}=\frac{1}{1+e_{}}\frac{\partial e}{\partial t}$$

(8)

Meanwhile, Equation (2) can be transformed to:

$$e=\left(1+e_0\right){\left(\frac{{{\sigma }^{\prime }}_0}{{\sigma }^{\prime }}\right)}^{I_c}-1$$

(9)

By taking the partial derivative of Equations (9) and (1) into Equation (8), we can obtain:

$$-\frac{1}{r}\frac{\partial }{\partial r}\left(r{k}_r\left(\mathrm{r}\right)\frac{\partial u}{{\gamma }_w\partial r}\right)=\frac{I_c}{{{\sigma }_0}^{\prime }}\frac{{{\sigma }_0}^{\prime }}{\overline{{\sigma }^{\prime }}}\left(\frac{\partial q(t)}{\partial t}-\frac{\partial u}{\partial t}\right)$$

(10)

Combining Equations (2) and (3) gives:

$$\frac{k_r}{k_{r0}}={(\frac{{{\sigma }^{\prime }}_0}{{\sigma }^{\prime }})}^{\alpha I_c}$$

(11)

Bringing Equation (11) into Equation (10), we can obtain:

$$\begin{array}{l}\frac{{{\sigma }^{\prime }}_0{k}_{r0}}{I_c}(\frac{{\sigma }^{\prime }}{{\sigma }_0^{\prime }})\{\frac{f(r)}{r}{\left(\frac{{\sigma }^{\prime }}{{\sigma }_0^{\prime }}\right)}^{-\alpha I_c}\left[\frac{\partial u}{{\gamma }_w\partial r}-i_0\left(1-e^{-\frac{1}{i_0}\frac{\partial u}{{\gamma }_w\partial r}}\right)\right]\\ +\frac{\partial }{\partial r}\left[f(r){\left(\frac{{\sigma }^{\prime }}{{\sigma }_0^{\prime }}\right)}^{-\alpha I_c}\left[\frac{\partial u}{{\gamma }_w\partial r}-i_0\left(1-e^{-\frac{1}{i_0}\frac{\partial u}{{\gamma }_w\partial r}}\right)\right]\right]\}=\frac{\partial u}{\partial t}-\frac{\partial q(t)}{\partial t}\end{array}$$

(12)

The corresponding boundary conditions are:

$$\{\begin{array}{l}u(r_w,z,t)=-p_v(z)\\ \frac{\partial u_r}{\partial r}{|}_{r=r_e}=0\\ u(r,z,0)=q(0)\end{array}$$

(13)

3. Solution of Equations

3.1. Dimensionless Parameters

To facilitate the solution, dimensionless processing is performed:

$$\begin{array}{c}S=\frac{{{\sigma }^{\prime }}_0}{{{\sigma }^{\prime }}_0},Q=\frac{q(t)}{{{\sigma }^{\prime }}_0},U=\frac{u}{{{\sigma }^{\prime }}_0},P_0=\frac{p_0}{{{\sigma }^{\prime }}_0},T_v=\frac{C_{v0}t}{4r_e^2},T_{v1}=\frac{C_{v0}{t}_1}{4r_e^2},C_{v0}=\frac{k_{r0}{{\sigma }_0}^{\prime }}{{\gamma }_w{I}_c},\hfill \\ I_0=\frac{i_0{\gamma }_w{r}_e}{{{\sigma }_0}^{\prime }},R=\frac{r}{r_e},R_s=\frac{r_s}{r_e},R_w=\frac{r_w}{r_e},Z=\frac{z}{H},P(Z)=\frac{p_v(z)}{{{\sigma }^{\prime }}_0}\hfill \end{array}$$

(14)

By introducing the above dimensionless equation into control Equation (13), we can obtain:

$$\begin{array}{l}\frac{1}{R}{(S+Q-U)}^{-\alpha I_c+1}\left[\frac{\partial U}{\partial R}-I_0\left(1-e^{-\frac{1}{I_0}\frac{\partial U}{\partial R}}\right)\right]\\ +\left(S+Q-U\right)\frac{\partial }{\partial R}\left[{(S+Q-U)}^{-\alpha I_c}\left(\frac{\partial U}{\partial R}-I_0\left(1-e^{-\frac{1}{I_0}\frac{\partial U}{\partial R}}\right)\right)\right]=\frac{\partial U}{\partial T}-\frac{\partial Q}{\partial T}\end{array}$$

(15)

The dimensionless form of the boundary conditions is as follows:

$$\{\begin{array}{l}U(R_w,Z,T_v)=-P(Z)\\ \frac{\partial u_r}{\partial r}(1,Z,T_v)=0\\ U(R,Z,0)=Q^0\end{array}$$

(16)

3.2. Differential Control Equations

The soil layer is divided into n equal parts along the radial direction, each equal part is divided into ΔR, the vertical direction is divided into m equal parts, each equal part is divided into ΔZ, and then the time is discretized by the ΔT_v. The Equation (15) can be discrete as:

$$\begin{array}{l}{U}_{i,m}^{j+1}=U_{i,b}^j+Q^{j+1}-Q^j+\frac{\mathsf{\Delta }{T}_v}{R_i}{\alpha }_{i,b}^j{\beta }_{i,b}^j\left[A_{i,b}^j-I_0\left(1-e^{-\frac{A_{i,b}^j}{I_0}}\right)\right]+\\ \frac{\mathsf{\Delta }{T}_v}{\mathsf{\Delta }R}\left[{\alpha }_{i,b}^j{\beta }_{i+1/2,b}^j\left(A_{i+1/2,b}^j-I_0\left(1-e^{-\frac{A_{i+1/2,b}^j}{I_0}}\right)\right)-{\alpha }_{i,b}^j{\beta }_{i-1/2,b}^j\left(A_{i-1/2,b}^j-I_0\left(1-e^{-\frac{A_{i-1/2,b}^j}{I_0}}\right)\right)\right]\end{array}$$

(17)

in which i, j, and b are the number of radial space nodes, time nodes, and vertical space nodes: $R_i=R_w+i\mathsf{\Delta }R$, ${\alpha }_{i,b}^j=(S+Q^j-U_{i,b}^j)$, ${\beta }_{i,b}^j={(S+Q^j-U_{i,b})}^{-\alpha I_c}$, $A_{i,b}^j=\frac{U_{i+1,b}^j-U_{i-1,b}^j}{2\mathsf{\Delta }R}$, $A_{i+1/2,b}^j=\frac{U_{i+1,b}^j-U_{i,b}^j}{\mathsf{\Delta }R}$, $A_{i-1/2,b}^j=\frac{U_{i,b}^j-U_{i-1,b}^j}{\mathsf{\Delta }R}$.

Boundary condition Equation (14) can be discretized as:

$$\{\begin{array}{c}{U}_{0,b}^j=-P_0[1-(1-k_1)(b\mathsf{\Delta }Z)]\\ U_{n+1,b}^j=U_{n-1,b}^j\\ U_{i,b}^0=Q^0\end{array}$$

(18)

Equations (17) and (18) constitute a closed system of equations, which can obtain the excess pore water pressure at any time and depth by iteration calculation. Meanwhile, for soil layers with a thickness of H, the settlement S_t at the moment of T_v can be expressed as:

$$S_t={\displaystyle {\int }_0^H{\displaystyle {\int }_{r_w}^{r_e}\frac{2r}{r_e^2-r_w^2}}\frac{(e_0-e)}{(1+e_0)}\mathrm{d}r}\mathrm{dz}$$

(19)

Substituting Equation (9) into Equation (19) discretely gave:

$$S_t={\displaystyle \sum _{b=0}^{m-1}\frac{G(b,T_v)+G(b+1,T_v)}{2}}$$

(20)

where $G(b,T_v)=\frac{{\displaystyle \sum _{i=0}^{n-1}{R}_i\left(1-{\left(S+Q-U_{i,b}\right)}^{-I_c}\right)+R_{i+1}\left(1-{\left(S+Q-U_{i+1,b}\right)}^{-I_c}\right)}}{(1+R_w^{})n}$.

The degree of consolidation U_s defined by deformation is:

$$U_s=\frac{S_t}{S_{\infty }}\times 100\%$$

(21)

In addition, the average consolidation degree Up defined by the pore pressure of the soft foundation with sand drains is:

$$U_p=1-\frac{{\displaystyle {\int }_0^H{\displaystyle {\int }_{r_w}^{r_e}2\pi rudr}}zdz}{{\displaystyle {\int }_0^H{\displaystyle {\int }_{r_w}^{r_e}2\pi r{q}_{u}drdz}}}=1-\frac{{\displaystyle {\int }_0^1{\displaystyle {\int }_{R_w}^{1}RUdRdZ}}}{{\displaystyle {\int }_0^1{\displaystyle {\int }_{R_w}^{1}R{Q}_{u}dRdZ}}}=1-{\displaystyle \sum _{b=0}^{m-1}\frac{{\displaystyle \sum _{i=0}^{n-1}{Y}_{i,b}+{\displaystyle \sum _{i=0}^{n-1}{Y}_{i,b+1}}}}{(1-R_w^2)Q_u}}$$

(22)

where $Y_{i,b}=\frac{U_{i,b}{R}_{i+1}-U_{i+1,b}{R}_i}{2\mathsf{\Delta }R}\left(R_{i+1}^2-R_i^2\right)+\frac{U_{i+1,b}-U_{i,b}}{3\mathsf{\Delta }R}\left(R_{i+1}^3-R_i^3\right)$.

4. Solution Verification

4.1. Verification with an Indoor Radial Penetration Test

Herein, the average consolidation degree defined by the settlement is compared with the indoor radial consolidation test results of Indraratna et al. [11] to verify the reliability of the proposed finite difference solutions.

The proposed system parameters need to be degraded to ensure the proposed model is consistent with the indoor test conditions. That is, reducing the proposed model with Darcy’s flow and 0 vacuum preloading scenario, i.e., I₀ = 0.00001, P₀ = 0 kPa. Meanwhile, the loading mode is set as constant loading, and the variation mode of the radial permeability coefficient takes mode 1 and $\delta =2/3$. Moreover, as shown in Figure 5 and Figure 6, the four sets of original data given by laboratory tests are re-fitted in the double logarithmic coordinate systems, i.e., determining I_c = 0.061 and α = 10.42. Other parameters are consistent with the test results, as summarized in

Table 1. The geotechnical parameters from indoor tests.

rw(m)	re(m)	rs(m)	kr0(m/s)	δ	σ0′(kPa)	e0	qu(kPa)
0.033	0.225	0.1	4.6 × 10−10	2/3	20	1	30

. Then, the degree of consolidation defined by the settlement is calculated according to Equation (21) and compared with the indoor monitored data, as shown in Figure 7. As observed, the consolidation degree U_s(%) calculated theoretically is in good agreement with the settlement results monitored by laboratory tests, preliminarily verifying the reliability of the proposed solutions.

4.2. Verification with Analytical Solutions under Equal Strain Assumption

Since there is currently no analytical solution for the large-strain sand-drained foundations under a double-logarithmic coordinate system, the analytical solution of equal strain is derived in this section by referring to the consolidation solution of the sand-drained foundations under equal strain conditions derived by Xie et al. [31]. The derived analytical solution assumption is consistent with the above-mentioned assumptions, except that it satisfies the equal strain assumption. For the sake of simplification, it is reduced to not considering the variation in radical permeability coefficient; that is, selecting the radical change in the permeability coefficient of the smear zone as mode 1, so that δ = 1, α = 0. Meanwhile, it degenerates into Darcy’s flow with I₀ = 0, and the loading mode is instantaneous loading. The boundary conditions, model parameters, and corresponding dimensionless variables are consistent with those mentioned above. Then, the radial consolidation equation is:

$$-\frac{1}{r}\frac{\partial }{\partial r}\left(r{k}_{r0}\frac{\partial u}{{\gamma }_w\partial r}\right)=-m_v\frac{\partial \overline{u}}{\partial t}=\frac{1}{1+e}\frac{de}{d\overline{{\sigma }^{\prime }}}\frac{\partial \overline{u}}{\partial t}$$

(23)

where m_v is the volume compression coefficient, $\overline{\mathrm{u}}$ is the average excess pore pressure of the sand-drained foundations, and the expression $\overline{\mathrm{u}}$ is as follows:

$$\overline{u}=\frac{1}{(r_e^2-r_w^2)}{\displaystyle {\int }_0^H{\displaystyle {\int }_{r_w}^{r_e}2ru\mathrm{d}r}}\mathrm{d}z$$

(24)

Bringing Equation (9) into Equation (24) gives:

$$-\frac{1}{r}\frac{\partial }{\partial r}\left(r{k}_{r0}\frac{\partial u_s}{{\gamma }_w\partial r}\right)=-\frac{I_c}{{{\sigma }_0}^{\prime }}\left(\frac{{{\sigma }_0}^{\prime }}{\overline{{\sigma }^{\prime }}}\right)\frac{\partial \overline{u_s}}{\partial t}$$

(25)

Combining Equation (14) and integrating both sides of Equation (25) yields:

$$\frac{\partial u_s}{\partial r}=-\frac{{{\sigma }_0}^{\prime }}{C_{v0}\overline{{\sigma }^{\prime }}}\frac{\partial \overline{u_s}}{\partial t}\frac{r_e^2-r^2}{2r}$$

(26)

Combining Equation (14) and integrating both sides of Equation (26) then gives:

$$u=\left(-\frac{{{\sigma }_0}^{\prime }}{2C_{v0}\overline{{\sigma }^{\prime }}}\frac{\partial \overline{u_s}}{\partial t}\right)\left(r_e^2\ln r-r_e^2\ln r_w-\frac{r^2-r_w^2}{2}\right)-p_v(z)$$

(27)

Bringing Equation (27) into Equation (24):

$$\overline{u}=\left(-\frac{{{\sigma }_0}^{\prime }}{2C_{v0}\overline{{\sigma }^{\prime }}}\frac{\partial \overline{u_s}}{\partial t}\right)G-\frac{1}{2}{p}_0$$

(28)

where $G=\left(\frac{r_e^4\ln r_e-r_e^4\ln r_w-\frac{3r_e^4}{4}-\frac{r_w^4}{4}+r_e^2{r}_w^2}{r_e^2-r_w^2}\right)$.

Separating the variables in Equation (28) yields:

$$\frac{\partial \overline{u_s}}{\partial t}=-\frac{2C_{v0}\left(\overline{{\sigma }_0^{\prime }}+q_u-\overline{u}\right)}{\overline{{\sigma }^{\prime }}}\frac{\left(\overline{u}+p_0\right)}{G}$$

(29)

Integrating t on both sides of Equation (31) gives:

$$\overline{u}=\overline{{\sigma }_0^{\prime }}+q_u-\frac{\overline{{\sigma }_0^{\prime }}(\overline{{\sigma }_0^{\prime }}+q_u+p_0/2)}{\overline{{\sigma }_0^{\prime }}+(q_u+p_0/2)e^{-\frac{2\overline{({\sigma }_0^{\prime }}+q_u+p_0/2)C_{v0}}{\overline{{\sigma }_0^{\prime }}G}t}}$$

(30)

The degree of consolidation defined by excess pore water pressure is:

$$U_p=1-\frac{\overline{u}}{q_u}$$

(31)

Substituting Equation (30) into Equation (31) and making it dimensionless then gives:

$$U_p=1-\frac{\overline{u}+p_0/2}{q_u+p_0/2}=-\frac{1}{Q+P_0/2}+\frac{1+Q+P_0/2}{\left(Q+P_0/2\right)\left(1+\left(Q+P_0/2\right)e^M\right)}$$

(32)

where $M=-8\left(1-R_w^2\right)\left(1+Q+\frac{P_0}{2}\right)T_v/\left(\mathrm{n}\frac{1}{R_w}-\frac{3}{4}-\frac{R_w^4}{4}+R_w^2\right)$.

The dimensionless parameters S = 1, Q =10, R_w =0.2, ΔR = 0.01, ΔT_v = 10⁻⁶, and P₀ = 5 are selected to verify the results obtained by finite difference solutions with the analytical solution under equal strain assumption. According to the comparison of solutions under the assumptions of equal strain and free strain given by Li [32], the impact of the two extreme assumptions on the consolidation of sand drains is mainly reflected in the early stage of consolidation. That is, the degree of consolidation obtained under the assumption of free strain is greater in the early consolidation stage but gradually reduces over time. As seen, the results shown in Figure 8 are consistent with this, further verifying the reliability of the proposed numerical solutions.

5. Parametric Analysis

In this part, the corresponding influencing factors of the proposed soft soil consolidation system with sand drain were analyzed, including the change modes of the radial permeability coefficient, the I_c and α obtained by linear fitting in the double-logarithmic compression and permeability logarithmic coordinates, and seepage parameter i0 in non-Darcy flow rule. To facilitate the calculation, and if there is no special statement, the basic parameters for the calculation and analysis of the sand-drained foundations are as follows: R_w = 0.2, Rs = 0.6, σ₀′ = 10 kPa, q₀ = 50 kPa, q_u = 100 kPa, ΔR = 0.01, P₀ = 5, I₀ = 1, Ic = 0.1, α = 10, δ = 2/3, ΔT_v = 10⁻⁶, ΔZ = 0.05, k₁ = 0, H = 1 m, the radial permeability coefficient variation mode is taken as pattern 1, and the loading mode is constant loading.

5.1. Excess Pore Water Pressure Distribution Analysis

Figure 9 and Figure 10 present the distribution curves of excess pore water pressure at Z = 0.5 under different loading modes when the dimensionless vacuum preloading value P₀ = 5 and 20, respectively. As observed, the excess pore pressure curves under instantaneous loading are smaller than the linear loading values. The slower the linear loading rate, the greater the corresponding excess pore water pressure values. In addition, by comparing Figure 9 and Figure 10, it can be observed that an increase in vacuum pressure simultaneously reduces the excess pore pressure values at the same location. The pore pressure at dimensionless radius R = 0.2 is the vacuum negative pressure value, while the pore pressure in other radical places gradually approaches the vacuum negative pressure over time. That is to say, the vacuum negative pressure is the final excess pore water pressure value, essentially at the completion of sand-drained foundation consolidation under vacuum preloading.

Figure 11 shows the excess pore pressure distribution at the intermediate soil layer (i.e., Z = 0.5) under three radial permeability coefficient variation modes. Meanwhile, Figure 12, Figure 13 and Figure 14 depict the excess pore pressure distribution curves along the radial radius r at the intermediate soil layer under a different compression index (I_c), permeability index ($\alpha$), and seepage parameter (I₀), respectively. It can be seen from Figure 11 that the pore pressure value in mode 1 is the highest, followed by mode 2, and that mode 3 is slightly lower than mode 2. That is to say, the consolidation rate is the fastest in the parabolic variation mode of the radial permeability coefficient within the smear zone, followed by the linear variation mode, and the slowest in the constant variation mode. Further, it can be seen that, with the increase in the compression index, permeability index, and non-Newtonian index, the excess pore pressure curve values all increase under the same radical locations. That is to say, the larger the compression index, permeability index, and non-Newtonian index, the slower the excess pore pressure dissipation proceeding the sand-drained foundations.

5.2. Analysis of Average Consolidation Degree and Settlement

Figure 15 shows the influence curves of three radial permeability coefficient variation modes on average consolidation and soil settlement. As seen, the permeability coefficient of the smear area shows a parabolic change, and the consolidation rate is the fastest (mode 3). And in mode 1, where the permeability coefficient of the smear zone is constant, the consolidation rate is the slowest. In essence, when the permeability coefficients within the undisturbed zone are the same, the larger the average permeability coefficient of the smear zone is, the larger the whole consolidation rate of the sand-drained foundations will be. For the sand-drained foundations’ settlement, although the final settlement amount achieved by the three modes is consistent, mode 3 takes the shortest time to reach the same settlement amount under settlement proceeding, and mode 1 takes the longest time. As a further comparison with mode 1, mode 2 has a maximum relative deviation of 32.1% for the average consolidation degree and 18.9% for soil settlement compared to mode 1, and the maximum relative deviation of the two consolidation indexes in mode 3, respectively, reached 47.8% and 29.6% compared to mode 1. That is to say, different radical permeability coefficient modes within the smear zone still have a significant impact on the consolidation rate of the sand-drained foundations, which also indicates that the permeability modes should be reasonably adopted in practical engineering.

Figure 16 shows the influence curves of the top surface dimensionless vacuum negative pressure P₀ on the radial consolidation of sand-drained foundations under different loading modes. As seen, the maximum vacuum negative pressure has a relatively minor impact on the consolidation rate of sand-drained foundations under instantaneous loading. For sand-drained foundations under single-stage loading, the larger the top surface vacuum negative pressure is, the faster the consolidation rate will be. Moreover, with the increase in the dimensionless single-stage loading completion time, T_v1, the influence of the top surface vacuum negative pressure on the consolidation rate becomes larger. Therefore, it can be concluded that, for actual soft foundation treatment projects with longer loading duration, more attention should be paid to the impact of vacuum preloading on the consolidation rate.

Li and Qiu’s [22] work has shown that the compression index (I_c) and permeability index (α) have the most concentrated range of variation, ranging from 0.08 to 0.12 and concentrated between 6 and 14, respectively, so this section adopts I_c = 6, 10, 14 and α = 0.08, 0.1, 0.12 for analysis. Figure 17 and Figure 18 show the influence of the average consolidation degree and settlement with time factor under different compression index (I_c) and permeability index (α) values. As observed, the whole consolidation rate slows with the increase in the compression index or permeability index. As for the settlement amount, the compressive index has little influence on the settlement amount at the early consolidation stage. However, as time goes on, the larger the compression index, the more significant the settlement at the same time, and the final settlement also becomes more significant. Meanwhile, with the increase in the permeability index (α), the time to reach the same settlement amount decreases, but the final settlement amount is unaffected. Nevertheless, it also should be noted that, compared to α = 6, the maximum relative deviation of consolidation degree reached 54.9% at α = 10, and it reached as high as 76.6% when α = 14. This indicates that although the permeability index, α, does not affect the final settlement, it has a significant impact on the whole consolidation process.

Figure 19 shows the average consolidation degree U_p corresponding to a different drain spacing ratio n, n = r_e/r_w. As can be seen, the radial consolidation rate decreases with an increase in the drain spacing ratio. Especially when n is greater than 40, the ratio has little influence on the consolidation rate. This also means that foundation soil beyond a certain range will not be affected by the spacing ratio n, which reflects that the range of the influence zone is limited. Therefore, in practical engineering, for rapid drainage and soft foundation consolidation with sand drains, the spacing of the drains should be reasonably arranged; n less than 40 is recommended.

Figure 20 illustrates the influence curves of the seepage parameter I₀ on the average degree of consolidation, U_p, and foundation settlement, St, at I₀ = 0, 1, 3, and 5. As seen, the consolidation and settlement rate with Darcy’s flow (I₀ = 0) is faster than that of the non-Newtonian index flow, and the whole consolidation and settlement rate becomes slower as I₀ increases. For example, the corresponding U_p values for increasing I₀ are 58.61%, 69.33%, 79.74%, and 86.29%, respectively, and the maximum relative deviation reached 38.2% when I₀ = 5 compared with Darcy’s flow. Therefore, without considering the influence of the non-Newtonian index flow, the consolidation and settlement rates of the sand-drained foundation will both be overestimated. Moreover, the settlement curves tend to be the same with time at different seepage parameters, I₀, indicating that the non-Darcy flow does not influence the final settlement.

Figure 21 shows the average consolidation and settlement curves under different negative pressure transfer coefficients k₁. It can be found that the change in k₁ has almost no effect on the whole consolidation degree of the sand-drained foundations. In contrast, the k₁ undeniably impacts foundation settlement at the middle and late consolidation stages. As seen, the soil settlement rate and final settlement decrease with the decrease of k₁, indicating that the long settlement deformation of actual engineering should also be paid more attention under vacuum preloading scenarios.

6. Limitations

In this study, a developed large-strain consolidation system has been presented by considering vacuum preloading, double-logarithmic compression permeability relationships, the non-Darcy flow effect, and the radial permeability coefficient variation. Nevertheless, both the vertical seepage of pore water and soil rheological properties during the consolidation process cannot be reflected in the system. Therefore, the simplified assumptions can be further improved by incorporating the effects of vertical seepage and rheological model relationships to reveal more realistic soil creep consolidation behavior [33,34,35]. Additionally, the complex boundary and initial conditions of stone column-improved soft clay [36,37,38] combined with the vacuum preloading method can also be an essential research branch for larger-strain consolidation of sand-drained foundations.

7. Conclusions

In this paper, an improved large-strain consolidation system of soft foundations with sand drain is established considering vacuum preloading, double-logarithmic compression permeability relationships, and the non-Darcy flow effect. The finite difference method was used to solve it, and the validity of the numerical solution of the model was verified. Based on this, the consolidation behavior of sand-drain-improved soft foundations was deeply explored, and the following main conclusions were obtained:

(1)

The larger the average permeability coefficient of the smear area, the greater the consolidation rate of the sand-drained foundations. Compared with the constant permeability coefficient, when the permeability coefficient of the smear area changes linearly, the maximum relative deviations of consolidation degree and settlement amount reach 32.1% and 18.9%, respectively. When the permeability coefficient of the smear area changes in a parabolic manner, the maximum relative deviations between the degree of consolidation and the settlement amount reach 47.8% and 29.6%, respectively.

(2)

The greater the maximum dimensionless vacuum negative pressure, P₀, the faster the consolidation rate of sand-drained foundations. Meanwhile, the more extended the loading duration, the greater the influence of vacuum negative pressure on the consolidation rate. This indicates that, in practical engineering projects that need to ensure foundation stability and fast consolidation, the vacuum preloading value can be adjusted promptly based on the proceedings of the construction period.

(3)

With an increase in the drain spacing ratio, n, the consolidation rate of soft foundations with sand drains decreases. Especially when the diameter ratio n > 40, increasing the healthy diameter ratio has little effect on the consolidation rate. Therefore, for practical engineering with a diameter ratio greater than 40, 40 can be used for calculations.

(4)

With the increase in the compression index and permeability index, the consolidation rate of sand-drained foundations slows down; the final settlement increases with the increase in the compression index, while the permeability index does not affect the final settlement.

(5)

Considering the seepage parameter, I₀, the consolidation rate of sand-drained foundations becomes slower, and the maximum deviation of the consolidation degree can reach 38.2%. Therefore, ignoring the non-Darcy flow effect will significantly overestimate the consolidation rate of the sand-drained foundations. Nevertheless, the seepage parameter I₀ does not affect the final settlement, indicating that the prediction of the final settlement result can be captured using a concise Darcy’s law.

(6)

The negative pressure transfer coefficient k₁ has little effect on the consolidation rate of sand-drain foundations, but a decrease in the negative pressure transfer coefficient, k₁, will result in a decrease in both the settlement rate and the final settlement amount of sand-drained foundations.