Analytical Solutions for the Radial Consolidation of Unsaturated Foundation with Prefabricated Vertical Drain Based on Fourier Series Expansion Theory

Abstract

This paper presents the analytical solution of the radial consolidation of a prefabricated vertical drain (PVD) foundation under the unsaturated condition. In the proposed modeling, air and water phases in the foundation are thought to dissipate horizontally toward to the drain, and the smear effect, drain resistance and external time-dependent loading are fully considered. The analytical mathematical tools, namely the general integration method, Fourier series expansion method, decoupling method and the constant variation method, are utilized to solve the partial differential equations. Moreover, the current solutions are verified with existing solutions in the literature. Finally, a case study considering the ramp loading and exponential loading is conducted to investigate the consolidation patterns under various loading parameters. The results show that smear effect and drain resistance can significantly hinder the dissipation process of excess pore pressures, and different external loading types will lead to various dissipation characteristics (i.e., peak values).

1. Introduction

The PVD technique is a traditional foundation treatment method for its high efficiency to accelerate the drainage of air and water in soil [1,2,3,4,5], and thus shorten the consolidation process of a foundation. In fact, the shorter drainage paths along the horizontal direction will be produced when a PVD is installed in a soil layer, resulting in the efficient radial dissipation of excess pore pressures. For the PVD, however, the construction process itself is an external disturbance, densifying the soil around the drain (i.e., smear effect). Besides, drain resistance will occur even though the drain is filled with highly permeable materials, which have been studied comprehensively in saturated soils [6,7,8,9,10,11,12,13,14]. Nevertheless, the soils in semi-arid regions around the world are always under a partially saturated state, whose consolidation results calculated with saturated consolidation theory could directly produce unaccepted errors.

In drought regions, the soils of foundation are always under the unsaturated state because of the low groundwater levels [15]. Under this circumstance, the designs of the foundations based on the consolidation theory in saturated soils may be uneconomical, because the dissipation time of excess pore pressures is misestimated. Therefore, Fredlund and his colleagues [16,17,18] researched the consolidation theory of unsaturated soils by the two-stress variables constitutive model, based on which the consolidation characteristics were analyzed by later scholars. For example, the variations of the excess pore-air and pore-water pressures were first researched by Qin et al. [19], in which the Laplace transform method is adopted to solve the partial differential equations of the unsaturated soil layer. Li et al. [20] proposed the semi-analytical solutions considering that the soil layer is multi-layered, and the boundary is time-dependent. Later, Shan et al. [21,22,23] analyzed the one-dimensional consolidation patterns of the soil layer under the unsaturated condition by considering that the foundation is single- and multi-layered, and the mixed boundary condition was investigated. With the consideration of various initial and boundary conditions, Zhou et al. [24,25,26] proposed more generalized investigations of the one-dimensional consolidation. In real foundation engineering, the upper surfaces are not ideally permeable or impermeable, and Wang et al. [27,28] comprehensively investigated the third type of boundary condition in unsaturated soils with semi-analytical methods. Afterwards, Ho and Fatahi [29] and Wang et al. [30,31] extended the solutions of consolidation from one-dimensional to two-dimensional situations. By taking the unsaturated and saturated conditions into account, Li et al. [32] presented semi-analytical solutions to illustrate the distributions of excess pore pressures in depth.

In 2010, a vertical drain foundation modeling under the partially saturated condition with the radial dissipation of air and water phases was first put forward by Qin et al. [33], in which the Laplace transform and the Crump method were applied to obtain the semi-analytical solution. Zhou presented the numerical solutions [34] (with the utilization of the differential quadrature method), and the analytical solutions [13] (by the decoupling method) of the unsaturated vertical drain foundation considering the radial dissipation. According to the difference of the initial values of the excess pore pressures along the depth, Ho et al. [35] derived the analytical solution to explore the relatively more practical consolidation characteristics with the linear-distributed initial excess pore pressures. Moreover, Ho and Fatahi [36]; Ho and Fatahi [37] researched the consolidation of an unsaturated vertical drain foundation considering the vertical dissipation, free-strain and equal-strain assumptions with analytical methods. More practically, Huang et al. [15,38] and Li et al. [39] analyzed the consolidation of a PVD foundation in unsaturated soil with the smear effect and drain resistance. Specifically, the analytical procedure in literature [38], based on the Fourier expansion method, the decoupling method and the constant variation method, which can efficiently solve the consolidation problems in the PVD foundation, is adopted in this paper.

In this paper, the analytical solution is proposed to investigate the consolidation of an unsaturated PVD foundation considering the smear effect, drain resistance and time-dependent loading. Firstly, the governing equations are rearranged as a new set of equations under the complete equal-strain assumption, and the boundary condition related to drain resistance is then adopted to obtain the equations concerning the excess pore pressures in the drain. Afterwards, the above-mentioned Fourier expansion method, the decoupling method and the constant variation method are conducted upon the equations to obtain the solutions in the time domain. It should be emphasized that the current solutions are not only neat and clear, but also can be computed with high efficiency because of the absence of the complicated numerical inversions.

2. Mathematical Modeling

The computational modeling is shown in Figure 1 [15,38], where a PVD is installed in the unsaturated soil layer, with a time-dependent loading acting on the upper surface. The vertical coordinate is downward, and the radial coordinate is horizontal. r_w, r_s and r_e are the largest radii of the drain, smeared zone and the influenced area, respectively. The thickness of the layered foundation is H. k_ad, k_wd k_as, k_ws, k_a and k_w are the permeability coefficients of air and water phases in the drain, smeared zone and influenced areas, respectively.

2.1. Basic Assumptions

In this paper, the solutions are derived based on the following assumptions:

• The soil layer is homogeneous.

• Solid particles and the water phase are incompressible.

• Air and water phases can only flow along the radial direction.

• The air and water phases are continuous and independent.

• Permeability coefficients related to air and water phases are constant.

2.2. Governing Equations

The governing equations of excess pore pressures in the unsaturated PVD foundation under the equal-strain assumption and time-dependent loading are concluded as the following matrix form:

(1)$C{\overline{u}}_t=-C_v\left(u_{rr}+\frac{1}{r}{u}_r\right)+C_q{q}_t$

where, $C=\left[\begin{array}{cc}1& C_{\mathrm{a}}\\ C_{\mathrm{w}}& 1\end{array}\right]$, $C_v=\left[\begin{array}{cc}{C}_{\mathrm{v}}^{\mathrm{a}}& 0\\ 0& C_{\mathrm{v}}^{\mathrm{w}}\end{array}\right]$, $C_q=\left[\begin{array}{cc}{C}_{\mathrm{q}}^{\mathrm{a}}& 0\\ 0& C_{\mathrm{q}}^{\mathrm{w}}\end{array}\right]$, $u=\left[\begin{array}{c}{u}_{\mathrm{a}}\\ u_{\mathrm{w}}\end{array}\right]$ and $q=\left[\begin{array}{c}q\\ q\end{array}\right]$. ${\left(•\right)}_r$ and ${\left(•\right)}_{rr}$ are the first and second derivatives of a function upon r, and ${\left(•\right)}_t$ denotes the first derivative of a function upon t. u_a and u_w are excess pore-air and water pressures respectively, q is the external loading, C_a and C_w, respectively, are the interaction constants of air and water phases, $C_{\mathrm{v}}^{\mathrm{a}}$ and $C_{\mathrm{v}}^{\mathrm{w}}$ are the consolidation coefficients of air and water phases and ${\overline{u}}_{\mathrm{a}}$ and ${\overline{u}}_{\mathrm{w}}$, respectively, are average excess pore-air and water pressures, defined as:

(2)$\overline{u}=\frac{{\displaystyle {\int }_{r_{\mathrm{w}}}^{r_{\mathrm{s}}}{u}_{s}2\pi rdr}+{\displaystyle {\int }_{r_{\mathrm{s}}}^{r_{\mathrm{e}}}u2\pi rdr}}{\pi \left(r_{\mathrm{e}}^2-r_{\mathrm{w}}^2\right)}$

where, $u_{\mathrm{s}}={\left[\begin{array}{cc}{u}_{{\mathrm{a}}_{\mathrm{s}}}& u_{{\mathrm{w}}_{\mathrm{s}}}\end{array}\right]}^{\mathrm{T}}$, u_as denotes excess pore-air pressure in the smeared area and u_ws is the excess pore-water pressure in the smeared area. The consolidation parameters are expressed as:

(3a)$C_{\mathrm{a}}=\frac{m_2^{\mathrm{a}}\left(u_{\mathrm{a}}^0+u_{\mathrm{atm}}\right)}{\left(m_1^{\mathrm{a}}-m_2^{\mathrm{a}}\right)\left(u_{\mathrm{a}}^0+u_{\mathrm{atm}}\right)-n(1-S_{\mathrm{r}})}$

(3b)$C_{\mathrm{v}}^{\mathrm{a}}=\frac{k_{\mathrm{a}}RT}{gM\left[\left(m_1^{\mathrm{a}}-m_2^{\mathrm{a}}\right)\left(u_{\mathrm{a}}^0+u_{\mathrm{atm}}\right)-n\left(1-S_{\mathrm{r}}\right)\right]}$

(3c)$C_{\mathrm{q}}^{\mathrm{a}}=\frac{m_1^{\mathrm{a}}\left(u_{\mathrm{a}}^0+u_{\mathrm{atm}}\right)}{\left(m_1^{\mathrm{a}}-m_2^{\mathrm{a}}\right)\left(u_{\mathrm{a}}^0+u_{\mathrm{atm}}\right)-n(1-S_{\mathrm{r}})}$

(3d)$C_{\mathrm{w}}=\frac{m_1^{\mathrm{w}}-m_2^{\mathrm{w}}}{m_2^{\mathrm{w}}}$

(3e)$C_{\mathrm{v}}^{\mathrm{w}}=\frac{k_{\mathrm{w}}}{{\gamma }_{\mathrm{w}}{m}_2^{\mathrm{w}}}$

(3f)$C_{\mathrm{q}}^{\mathrm{w}}=\frac{m_1^{\mathrm{w}}}{m_2^{\mathrm{w}}}$

where, $m_1^{\mathrm{a}}$, $m_2^{\mathrm{a}}$, $m_1^{\mathrm{w}}$ and $m_2^{\mathrm{w}}$, respectively, are the variations of net stress (q₀-u_a) and suction (u_a-u_w), u_atm is the atmospheric pressure, and u_atm = 101.3 kPa, M is the air mass molecular, and M = 0.029 kg/mol, and n₀ and S_r are the initial porosity and saturation, respectively. R is the universal air constant, R = 8.314 J/(mol·K), T is the absolute temperature, and T = 293 K, g is the gravitational acceleration, that is g = 9.8 m/s², and γ_w is the unit weight of water, and γ_w = 9.8 kN/m³.

2.3. Initial and Boundary Conditions

2.3.1. Initial Conditions

(4)${\left.u\right|}_{t=0}={\left.u_d\right|}_{t=0}=u_0$

where, $u_0={\left[\begin{array}{cc}{u}_{\mathrm{a}}^0& u_{\mathrm{w}}^0\end{array}\right]}^{\mathrm{T}}$, $u_d={\left[\begin{array}{cc}{u}_{{\mathrm{a}}_{\mathrm{d}}}& u_{{\mathrm{w}}_{\mathrm{d}}}\end{array}\right]}^{\mathrm{T}}$. $u_{\mathrm{a}}^0$ and $u_{\mathrm{w}}^0$, respectively, are the initial excess pore-air and pore-water pressures under the applied loading at t = 0, and u_ad and u_wd are excess pore-air and pore-water pressures in the drain, respectively.

2.3.2. Boundary Conditions

Continuation conditions at r = r_w:

(5)${\left.u_s\right|}_{r=r_w}=u_d$

(6)${\left.u_{s,r}\right|}_{r=r_w}=-\frac{r_w}{2}K{u}_{d,zz}$

where, $K=\left[\begin{array}{cc}\frac{k_{{\mathrm{a}}_{\mathrm{d}}}}{k_{{\mathrm{a}}_{\mathrm{s}}}}& 0\\ 0& \frac{k_{{\mathrm{w}}_{\mathrm{d}}}}{k_{{\mathrm{w}}_{\mathrm{s}}}}\end{array}\right]$.

It is noted that the excess pore pressures u_ad and u_wd, which are defined as the functions of variables z and t, are induced by the vertical dissipation in the drain instead of the foundation.

Continuation conditions at r = r_s:

(7)${\left.u_s\right|}_{r=r_{\mathrm{s}}}={\left.u\right|}_{r=r_{\mathrm{s}}}$

Boundary conditions at r = r_e:

(8)${\left.u_r\right|}_{r=r_e}=0$

Boundary conditions at z = 0:

(9)${\left.u\right|}_{z=0}={\left.u_d\right|}_{z=0}=0$

Boundary conditions z = H:

(10)${\left.u_z\right|}_{z=H}={\left.u_{d,z}\right|}_{z=H}=0$

3. Analytical Solutions

Firstly, Equation (1) needs to be adopted to transform u_a and u_w into ${\overline{u}}_{\mathrm{a}}$ and ${\overline{u}}_{\mathrm{w}}$ respectively, to introduce the equal-strain theorem into the solution. Rearranging Equation (1) produces:

(11)$\frac{\partial }{\partial r}\left(r{u}_r\right)=-r{C}_v{}^{-1}\left(C{\overline{u}}_t-C_q{q}_t\right)$

Incorporated with Equations (5), (7) and (8), the following expressions of excess pore pressures in the smeared and influenced areas can be obtained by integrating Equation (11) upon r in intervals [r_w, r] and [r_s, r]:

(12)$u_s=\frac{1}{2}\left[r_e^2\ln \left(\frac{r}{r_w}\right)-\frac{r^2-r_w^2}{2}\right]C_{v_s}{}^{-1}\left(C{\overline{u}}_{,t}-C_q{q}_{,t}\right)+u_d$

(13)$u=\frac{1}{2}\left[r_e^2\ln \left(\frac{r}{r_s}\right)-\frac{r^2-r_s^2}{2}\right]C_{v_r}{}^{-1}\left(C{\overline{u}}_{,t}-C_q{q}_{,t}\right)+{\left.u_s\right|}_{r=r_s}$

By substituting Equations (12) and (13) into Equations (1) and (2), one can produce the expression of excess pore pressure under the equal-strain assumption [15]:

(14)$\overline{u}=\frac{r_e^2}{2}F{C}_v{}^{-1}\left(C{\overline{u}}_t-C_q{q}_t\right)+u_{\mathrm{d}}$

where, $F=\left[\begin{array}{cc}\begin{array}{l}{F}_{\mathrm{a}}\\ 0\end{array}& \begin{array}{l}0\\ F_{\mathrm{w}}\end{array}\end{array}\right]$. F_a and F_w are the intermediate variables concerning smear effects:

(15)$F_{\mathrm{a}}=\frac{N^2}{N^2-1}\left[\left(\ln \frac{N}{S}-\frac{3}{4}-\frac{S^4}{4N^4}+\frac{S^2}{N^2}\right)+{\alpha }_{\mathrm{a}}\left(\frac{S^4-1}{4N^4}-\frac{S^2-1}{N^2}+\ln S\right)\right]$

(16)$F_{\mathrm{w}}=\frac{N^2}{N^2-1}\left[\left(\ln \frac{N}{S}-\frac{3}{4}-\frac{S^4}{4N^4}+\frac{S^2}{N^2}\right)+{\alpha }_{\mathrm{w}}\left(\frac{S^4-1}{4N^4}-\frac{S^2-1}{N^2}+\ln S\right)\right]$

α_a and α_w are the smear coefficients regarding air and water phases respectively, that is α_a = k_a/k_as and α_w = k_w/k_ws, and S = r_s/r_w and N = r_e/r_w.

It is obvious that the derivation of item $u_{\mathrm{d}}$ in Equation (14) is the key to obtaining the final solutions of $\overline{u}$. Therefore, substituting Equation (12) into Equation (6) affords:

(17)$u_{d,zz}=\left(1-N^2\right)K^{-1}{C}_{v_s}{}^{-1}\left(C{\overline{u}}_t-C_q{q}_t\right)$

Combining Equations (14) and (17) produces the following equation:

(18)$u_{d,zz}=\mathsf{\rho }\left(u_d-\overline{u}\right)$

where, $\mathsf{\rho }=\left[\begin{array}{cc}{\rho }_{\mathrm{a}}& 0\\ 0& {\rho }_{\mathrm{w}}\end{array}\right]$, ${\rho }_{\mathrm{a}}=\frac{8\left(N^2-1\right)}{N^2{H}^2}\frac{G_{\mathrm{a}}}{F_{\mathrm{a}}}$, ${\rho }_{\mathrm{w}}=\frac{8\left(N^2-1\right)}{N^2{H}^2}\frac{G_{\mathrm{w}}}{F_{\mathrm{w}}}$. G_a and G_w are the factors of drain resistance of air and water phases respectively, and $G_{\mathrm{a}}={\left(\frac{H}{d_{\mathrm{w}}}\right)}^2\frac{k_{\mathrm{a}}}{k_{{\mathrm{a}}_{\mathrm{d}}}}$, $G_{\mathrm{w}}={\left(\frac{H}{d_{\mathrm{w}}}\right)}^2\frac{k_{\mathrm{w}}}{k_{{\mathrm{w}}_{\mathrm{d}}}}$; $d_{\mathrm{w}}=2r_{\mathrm{w}}$.

Furthermore, combining Equations (14) and (18) provides the exact expression of u_d:

(19)$u_{d,zzt}=\mathsf{\rho }\left(u_{d,t}+\mathsf{\lambda }{\mathsf{\rho }}^{-1}{u}_{d,zz}-{\mathsf{\lambda }}_q{q}_t\right)$

where, $\mathsf{\lambda }=\frac{2}{r_{\mathrm{e}}^2}{C}^{-1}{C}_{v_r}{F}^{-1}$, and ${\mathsf{\lambda }}_q=C^{-1}{C}_q$.

It is noted that the only function $u_d$ to be solved appears in Equation (19), thus the solution of $u_d$ can be obtained explicitly hereinafter. Based on the Fourier series expansion theory and boundary conditions (i.e., Equations (9) and (10)), $u_d$ is expressed as:

(20)$u_d={\displaystyle \sum _{k=1}^{\infty }{U}_d\left(t\right)\sin \left(\frac{K}{H}z\right)}$

where, $K=\frac{2k+1}{2}\pi , k=1,2,3\dots$

By substituting Equations (19) into (20), one can obtain:

(21)${\displaystyle \sum _{k=1}^{\infty }{U}_{d,t}\left(t\right){\left(\frac{K}{H}\right)}^2\sin \left(\frac{K}{H}z\right)}=-\mathsf{\rho }\left[{\displaystyle \sum _{k=1}^{\infty }{U}_{d,t}\left(t\right)\sin \left(\frac{K}{H}z\right)}-\mathsf{\lambda }{\mathsf{\rho }}^{-1}{\displaystyle \sum _{k=1}^{\infty }{U}_d\left(t\right){\left(\frac{K}{H}\right)}^2\sin \left(\frac{K}{H}z\right)}-{\mathsf{\lambda }}_q{q}_t\right]$

By applying the orthogonality of trigonometric functions, Equation (21) can be converted into:

(22)$U_{d,t}\left(t\right)=\left[\begin{array}{cc}{A}_{\mathrm{a}}& -A_{\mathrm{w}}\\ -W_{\mathrm{a}}& W_{\mathrm{w}}\end{array}\right]U_d\left(t\right)+\left[\begin{array}{cc}{A}_{\mathrm{q}}& 0\\ 0& W_{\mathrm{q}}\end{array}\right]q_t$

where,

$A_{\mathrm{a}}=\frac{{\lambda }_{\mathrm{a}}{K}^2}{K^2+{\rho }_{\mathrm{a}}{H}^2},A_{\mathrm{w}}=\frac{C_{\mathrm{a}}{\lambda }_{\mathrm{w}}{\rho }_{\mathrm{a}}{K}^2}{{\rho }_{\mathrm{w}}\left(K^2+{\rho }_{\mathrm{a}}{H}^2\right)},W_{\mathrm{a}}=\frac{C_{\mathrm{w}}{\lambda }_{\mathrm{a}}{\rho }_{\mathrm{w}}{K}^2}{{\rho }_{\mathrm{a}}\left(K^2+{\rho }_{\mathrm{w}}{H}^2\right)};$

$W_{\mathrm{w}}=\frac{{\lambda }_{\mathrm{w}}{K}^2}{K^2+{\rho }_{\mathrm{w}}{H}^2},A_{\mathrm{q}}=-\frac{2}{K}\frac{{\rho }_{\mathrm{a}}{\lambda }_{\mathrm{a}}^{\mathrm{q}}{H}^2}{K^2+{\rho }_{\mathrm{a}}{H}^2},W_{\mathrm{q}}=-\frac{2}{K}\frac{{\rho }_{\mathrm{w}}{\lambda }_{\mathrm{w}}^{\mathrm{q}}{H}^2}{K^2+{\rho }_{\mathrm{w}}{H}^2}.$

It is appointed that the expression of $U_d\left(t\right)$, as shown is Equation (22), is arranged as the form in accordance with that of [15] due to the similarity of the partial differential equations, and the differences occur in the intermediate parameters (i.e., A_a, A_w, W_a, W_w, A_q and W_q).

Furthermore, the solution for $U_d\left(t\right)$ according to the decoupling method proposed by Qin et al. (2010), and the constant variation method by Huang et al. [38], is:

(23)$U_d\left(t\right)=\frac{1}{q_{12}{q}_{21}-1}\left\{\left[\begin{array}{cc}-C_1& q_{21}{C}_2\\ q_{12}{C}_1& -C_2\end{array}\right]\left[\begin{array}{c}{e}^{Q_{1}t}\\ e^{Q_{2}t}\end{array}\right]+{\displaystyle {\int }_0^t\left[\begin{array}{cc}-{\beta }_1& q_{21}{\beta }_2\\ q_{12}{\beta }_1& -{\beta }_2\end{array}\right]\left[\begin{array}{c}{e}^{Q_1\left(t-\tau \right)}\\ e^{Q_2\left(t-\tau \right)}\end{array}\right]\frac{\partial q}{\partial \tau }}d\tau \right\}$

where, $q_{12}=\frac{W_{\mathrm{a}}}{A_{\mathrm{a}}-Q_2}=\frac{W_{\mathrm{w}}-Q_2}{A_{\mathrm{w}}}$, $q_{21}=\frac{A_{\mathrm{a}}-Q_1}{W_{\mathrm{a}}}=\frac{A_{\mathrm{w}}}{W_{\mathrm{w}}-Q_1}$, ${\beta }_1=A_{\mathrm{q}}+q_{21}{W}_{\mathrm{q}}$, ${\beta }_2=q_{12}{A}_{\mathrm{q}}+W_{\mathrm{q}}$; $Q_{1,2}=\frac{1}{2}\left[\left(A_{\mathrm{a}}+W_{\mathrm{w}}\right)\pm \sqrt{{\left(A_{\mathrm{a}}-W_{\mathrm{w}}\right)}^2+4W_{\mathrm{a}}{A}_{\mathrm{w}}}\right]$. C₁ and C₂ are the unknown constants determined by initial conditions.

Finally, the analytical solutions for the variations of excess pore pressures can be easily obtained based on Equation (18) and the initial conditions (i.e., Equation (4)):

(24)$\overline{u}=\frac{1}{q_{12}{q}_{21}-1}{\displaystyle \sum _{m=1}^{\infty }\left(\Lambda \left[\begin{array}{c}{\mathrm{e}}^{Q_{1}t}\\ {\mathrm{e}}^{Q_{2}t}\end{array}\right]+\Theta \right)\sin \left(\frac{M}{H}z\right)}$

where,

$\Lambda =\left[\begin{array}{cc}-{\mathsf{\Psi }}_{\mathrm{a}}& {\mathsf{\Omega }}_{\mathrm{a}}\\ {\mathsf{\Omega }}_{\mathrm{w}}& -{\mathsf{\Psi }}_{\mathrm{w}}\end{array}\right],\Theta =\left[\begin{array}{c}{\mathsf{\Theta }}_{\mathrm{a}}\\ {\mathsf{\Theta }}_{\mathrm{w}}\end{array}\right],{\mathsf{\Omega }}_{\mathrm{a}}=\frac{2}{K}{q}_{21}\left(q_{12}{u}_{\mathrm{a}}^0+\frac{{\chi }_{\mathrm{a}}}{{\chi }_{\mathrm{w}}}{u}_{\mathrm{w}}^0\right),{\mathsf{\Psi }}_{\mathrm{a}}=\frac{2}{K}\left(u_{\mathrm{a}}^0+q_{21}\frac{{\chi }_{\mathrm{a}}}{{\chi }_{\mathrm{w}}}{u}_{\mathrm{w}}^0\right);$

${\mathsf{\Omega }}_{\mathrm{w}}=\frac{2}{K}{q}_{12}\left(\frac{{\chi }_{\mathrm{w}}}{{\chi }_{\mathrm{a}}}{u}_{\mathrm{a}}^0+q_{21}{u}_{\mathrm{w}}^0\right),{\mathsf{\Psi }}_{\mathrm{w}}=\frac{2}{K}\left(q_{12}\frac{{\chi }_{\mathrm{w}}}{{\chi }_{\mathrm{a}}}{u}_{\mathrm{a}}^0+u_{\mathrm{w}}^0\right),{\chi }_{\mathrm{a}}=1+\frac{K^2}{{\rho }_{\mathrm{a}}{H}^2},{\chi }_{\mathrm{w}}=1+\frac{K^2}{{\rho }_{\mathrm{w}}{H}^2};$

${\mathsf{\Theta }}_{\mathrm{a}}={\chi }_{\mathrm{a}}{\displaystyle {\int }_0^t\left[q_{21}{\beta }_2{\mathrm{e}}^{Q_2\left(t-\tau \right)}-{\beta }_1{\mathrm{e}}^{Q_1\left(t-\tau \right)}\right]\frac{\partial q}{\partial \tau }}d\tau ,{\mathsf{\Theta }}_{\mathrm{w}}={\chi }_{\mathrm{w}}{\displaystyle {\int }_0^t\left[q_{12}{\beta }_1{\mathrm{e}}^{Q_1\left(t-\tau \right)}-{\beta }_2{\mathrm{e}}^{Q_2\left(t-\tau \right)}\right]\frac{\partial q}{\partial \tau }}d\tau .$

Based on the consolidation theory by Fredlund and Hasan [16], the average consolidation degree can be expressed as:

(25)$\overline{U}=\frac{\left(m_2^{\mathrm{s}}-m_1^{\mathrm{s}}\right)\left(\frac{1}{H}{\displaystyle {\int }_0^H{\overline{u}}_{\mathrm{a}}dz}-u_{\mathrm{a}}^0\right)-m_2^{\mathrm{s}}\left(\frac{1}{H}{\displaystyle {\int }_0^H{\overline{u}}_{\mathrm{w}}dz}-u_{\mathrm{w}}^0\right)+m_1^{\mathrm{s}}\left(q-q_0\right)}{\left(m_1^{\mathrm{s}}-m_2^{\mathrm{s}}\right)u_{\mathrm{a}}^0+m_2^{\mathrm{s}}{u}_{\mathrm{w}}^0++m_1^{\mathrm{s}}\left(q_{\mathrm{f}}-q_0\right)}$

where, $m_1^{\mathrm{s}}=m_1^{\mathrm{a}}+m_1^{\mathrm{w}}$, $m_2^{\mathrm{s}}=m_2^{\mathrm{a}}+m_2^{\mathrm{w}}$; $q_{\mathrm{f}}$ is the final value of loading when time approaches to infinity.

4. Verification Work

In this section, the validity of the proposed solutions is verified by existing results in the literature [15]. It should be noted that the solutions were derived based on the instantaneous loading by Huang et al. [15], thus q = q₀ is adopted in current solutions to perform the verification work. Practical values of variables are as follows:

$$\begin{gathered} q_0=100\mathrm{kPa},u_{\mathrm{a}}^0=20\mathrm{kPa},u_{\mathrm{w}}^0=40\mathrm{kPa},k_{\mathrm{w}}={10}^{-10}\mathrm{m}/\mathrm{s},{\alpha }_{\mathrm{a}}={\alpha }_{\mathrm{w}}=1; \\ {r}_{\mathrm{w}}=r_{\mathrm{s}}=0.2\mathrm{m},r_{\mathrm{e}}=1.8\mathrm{m},H=10\mathrm{m},u_{\mathrm{atm}}=101.3\mathrm{kPa},m_1^{\mathrm{a}}=-2\times {10}^{-4}\mathrm{kPa}; \\ {m}_1^{\mathrm{w}}=-5\times {10}^{-5}\mathrm{kPa},S_{\mathrm{r}}=80\%,n_0=50\%,m_2^{\mathrm{w}}=-2\times {10}^{-4}\mathrm{kPa},m_2^{\mathrm{a}}=-1\times {10}^{-4}\mathrm{kPa}. \end{gathered}$$

The validity of the current solution is conducted by the existing solutions of reference [15], in which the consolidation modeling is similar to the one proposed in this paper. The relevant parameters are adopted as: S = 1, N = 9, α_a = α_w = 1 and G_a = G_w = 0.2. It can be obviously observed from Figure 2 that the results in this paper and those in [15] are perfectly consistent, which can enhance the reliability of the solutions proposed in this paper. Moreover, the dissipation speed of excess pore-air (whole period) and pore-water pressures (at the earlier period) varies with the k_a/k_w ranging from 0.1 to 1000. In particular, the increasing values of the ratio k_a/k_w imply the faster dissipation of excess pore pressures, and vice versa. On the other hand, the proposed solutions can be smoothly converted into the reduced solutions under certain boundary conditions, such as those of [15], since the two research frames (i.e., the proposed solutions and those in the literature) are similar.

5. Case Study

The dissipation performance of the excess pore pressures is investigated in the subsequent section, in which the common ramp and exponential loading in practical engineering are adopted. For the influencing factors, modeling sizes (i.e., S and N) and hydraulic-related parameters (i.e., α_a, α_w, G_a and G_w) upon the consolidation of the unsaturated PVD foundation have been carried out comprehensively in the existing literature [15], and these parameters will not be discussed in this paper. The two time-dependent loadings are shown in Equations (26) and (27).

Ramp loading:

(26)$q\left(t\right)=\left\{\begin{array}{l}{q}_0+at,0<t<\frac{q_0}{a}\\ 2q_0,\frac{q_0}{a}\le t\end{array}\right.$

Exponential loading:

(27)$q\left(t\right)=q_0+q_0\left(1-e^{-bt}\right)$

where a and b are the loading parameters of ramp and exponential loadings, respectively. Other parameter values are consistent with those noted above in the verification section if they are not involved in the special discussion. Substituting Equations (26) and (27) into Equations (23)–(25) can produce the final solutions under certain circumstances.

The variations of the excess pore pressures in the unsaturated PVD foundation under different ramp loading parameters are illustrated in Figure 3, and the loading parameters $a=5\times {10}^{-2}$, $5\times {10}^{-3}$, $5\times {10}^{-4}$ and $5\times {10}^{-5}$ kPa/s are adopted. It can be concluded from Figure 3 that the loading parameters can lead to different dissipation rates of excess pore pressures during the dissipation process. Nevertheless, the final closure time seems to not be influenced, and this phenomenon is more obvious in Figure 3b. The first dissipation period of excess pore-water pressure (i.e., in the time interval 10¹ < t < 10⁶) is influenced by the change of excess pore-air pressure, and the variations of the dissipation rates are attributed to the changing disparities between the applied loading and the excess pore-air (or pore-water) pressure.

In Figure 4, the dissipation process of excess pore pressures under the exponential loading is researched, with loading parameters $b=5\times {10}^{-2}$, $5\times {10}^{-3}$, $5\times {10}^{-4}$, $5\times {10}^{-5}$ and $5\times {10}^{-6}$ s⁻¹. Similarly, the loading parameter b of exponential loading poses significant influences upon the excess pore pressures, and the bigger values of b symbolize faster dissipation of excess pressures during a certain period.

Specifically, the upward trends in the dissipation curves appear when the increase of the ramp or exponential loadings exceeds the reduction of excess pore pressures, and the excess pore pressures will dissipate (or reduce) gradually when the loading curves rise to the constant level. More importantly, the higher peak values occur when slower dissipation of the excess pore pressures is caused by the ramp or exponential loading (or the increase of the loading parameters). Above all, the increase of the loading parameters means the faster rise or decline of excess pore pressures at the loading influenced period, which is due to the same terminal dissipation time with different peak values. This phenomenon is harmful for real construction, because the sudden variations of excess pore pressures will lead to unpredicted accidents (i.e., the sudden rise of the excess pore pressures will lead to liquefaction of soil, and the sudden decline will cause a drop of the foundation at short notice).

The consolidation characteristics of the unsaturated PVD foundation are shown in Figure 5, in which the loading parameters, $a=5\times {10}^{-2}$, $5\times {10}^{-3}$, $5\times {10}^{-4}$ and $5\times {10}^{-5}$ kPa/s, $b=5\times {10}^{-2}$, $5\times {10}^{-3}$, $5\times {10}^{-4}$, $5\times {10}^{-5}$ and $5\times {10}^{-6}$ and s⁻¹ are applied. It can be observed from Figure 5 that the bigger values of loading parameters lead to faster consolidation speed during the first period. In particular, the consolidation curves can be divided into three periods when $a>5\times {10}^{-3}$ kPa/s and $b>5\times {10}^{-4}$ s⁻¹, and only two parts are involved in other cases. We think that the risk of the sudden drop of the foundation should take the excess pore pressure variations, settlement magnitude and the drop speed into account. Although the consolidation process under loading parameter $a=5\times {10}^{-2}$ kPa/s might be more reliable in practical engineering, the above-mentioned conclusions in Figure 3 and Figure 4 should also be considered. In conclusion, the loading parameters should be chosen appropriately to meet the safety and reliability in practical PVD foundation engineering.

6. Concluding Remarks

By taking the smear effect, drain resistance and external time-dependent loading into consideration, this paper presented analytical solutions for the PVD foundation under the unsaturated condition. The Fourier series expansion method, decoupling method and constant variation method were utilized in the derivative process, which can present the final solutions in neat and explicit forms, and verification work was conducted to prove the correctness of the solutions via existing solutions. To some extent, the proposed solutions filled the research gap of the unsaturated PVD foundation in analytical procedures by considering different external boundary conditions and time-dependent loading simultaneously.

Smear effect and drain resistance can significantly influence the dissipation process of excess pore pressures, and the dissipation curves resemble certain external loading types. Compared with those semi-analytical solutions with complex numerical transforms and inversion processes (i.e., the Laplace transform and the inverse Laplace transform), the current solutions can be computed with high efficiency and accuracy. Additionally, the proposed analytical procedure may provide powerful guidelines to solve relevant coupled partial differential equations.

Author Contributions

Conceptualization, T.L.; writing—original draft preparation, Q.M.; writing—review and editing, Y.C., Q.X. and X.L. All authors have read and agreed to the published version of the manuscript.

Funding

We sincerely acknowledge the National Natural Science Foundation of China (Grant No. 41807232) for financial support.

Institutional Review Board Statement

This study did not require ethical approval.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study.

Conflicts of Interest

The authors declare that there are no conflict of interest.

Nomenclature

$$q_0$$	the applied external loading at t = 0
$$a$$	the ramp loading parameter
$$b$$	the exponential loading parameter
$$r_{\mathrm{w}}$$	the radius of the drain
$$r_{\mathrm{s}}$$	the largest radius of the smeared zone
$$r_{\mathrm{e}}$$	the largest radius of the influenced region
$$H$$	the thickness of the foundation
$$u_{\mathrm{a}}$$	the excess pore-air pressure
$$u_{\mathrm{w}}$$	the excess pore-water pressure
$$u_{{\mathrm{a}}_{\mathrm{s}}}$$	the excess pore-air pressure in the smeared zone
$$u_{{\mathrm{w}}_{\mathrm{s}}}$$	the excess pore-water pressure in the smeared zone
$$u_{\mathrm{a}}^0$$	the initial excess pore-air pressure
$$u_{\mathrm{w}}^0$$	the initial excess pore-water pressure
$${\overline{u}}_{\mathrm{a}}$$	the average excess pore-air pressure
$${\overline{u}}_{\mathrm{w}}$$	the average excess pore-water pressure
$$k_{\mathrm{a}}$$	the air-related permeability coefficient of influenced area
$$k_{\mathrm{w}}$$	the water-related permeability coefficient of influenced area
$$k_{{\mathrm{a}}_{\mathrm{s}}}$$	the air-related permeability coefficient of smeared zone
$$k_{{\mathrm{w}}_{\mathrm{s}}}$$	the water-related permeability coefficient of smeared zone
$$k_{{\mathrm{a}}_{\mathrm{d}}}$$	the air-related permeability coefficient in the drain
$$k_{{\mathrm{w}}_{\mathrm{d}}}$$	the water-related permeability coefficient in the drain
$${\alpha }_{\mathrm{a}}$$	the air-related smear coefficient
$${\alpha }_{\mathrm{w}}$$	the water-related smear coefficient
$$G_{\mathrm{a}}$$	the air-related drain resistance factor
$$G_{\mathrm{w}}$$	the water-related drain resistance factor
$$C_{\mathrm{a}}$$	the water-related interaction constant related to air phase
$$C_{\mathrm{w}}$$	the water-related interaction constant related to water phase
$$C_{\mathrm{v}}^{\mathrm{a}}$$	the air-related consolidation coefficient
$$C_{\mathrm{v}}^{\mathrm{w}}$$	the water-related consolidation coefficient
$$m_1^{\mathrm{a}}$$	the air-related parameter of net normal stress, q0-ua
$$m_2^{\mathrm{a}}$$	the air-related parameter of the net normal stress, q0-ua
$$m_1^{\mathrm{w}}$$	the water-related parameter of the net normal stress, q0-ua, and suction, ua-uw
$$m_2^{\mathrm{w}}$$	the water-related parameter of the net normal stress, q0-ua, and suction, ua-uw
$$u_{\mathrm{atm}}$$	the atmospheric pressure
$$M$$	the air mass molecular
$$n$$	the porosity
$$S_{\mathrm{r}}$$	the saturation
$$R$$	the universal air constant
$$T$$	the absolute temperature
$$g$$	the gravitational acceleration
$${\gamma }_{\mathrm{w}}$$	the water unit weight

Footnotes

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Consolidation modeling of unsaturated PVD foundation under time-dependent loading.

Figure 2. The comparison results of the current solutions versus existing solutions upon: (a) excess pore-air and (b) pore-water pressures.

Figure 3. Dissipation of (a) excess pore-air and (b) pore-water pressures against various ramp loading parameters, a.

Figure 4. Dissipation of (a) excess pore-air and (b) pore-water pressures against various exponential loading parameters, b.

Figure 5. Consolidation patterns of the PVD foundation under (a) ramp and (b) exponential loadings against various loading parameters of a and b.

References

1. Leo, C.J. Consolidation of fine-grained compressible soils by vertical drains based on the free strain approach. Appl. Mech.; 2002; pp. 207-212. [DOI: https://dx.doi.org/10.1142/9789812777973_0033]

2. Nogami, T.; Li, M.X. Consolidation of clay with a system of vertical and horizontal drains. J. Geotech. Geoenviron. Eng.; 2003; 129, pp. 838-848. [DOI: https://dx.doi.org/10.1061/(ASCE)1090-0241(2003)129:9(838)]

3. Leo, C.J. Equal strain consolidation by vertical drains. J. Geotech. Geoenviron. Eng.; 2004; 130, pp. 316-327. [DOI: https://dx.doi.org/10.1061/(ASCE)1090-0241(2004)130:3(316)]

4. García-Ros, G.; Alhama, I.; Cánovas, M. Powerful software to simulate soil consolidation problems with prefabricated vertical drains. Water; 2018; 10, 242. [DOI: https://dx.doi.org/10.3390/w10030242]

5. García-Ros, G.; Alhama, I.; Cánovas, M. Use of discriminated nondimensionalization in the search of universal solutions for 2-D rectangular and cylindrical consolidation problems. Open Geosci.; 2018; 10, pp. 209-221. [DOI: https://dx.doi.org/10.1515/geo-2018-0016]

6. Chen, Z.; Ni, P.P.; Mei, G.X.; Chen, Y.F. Semi-Analytical Solution for Consolidation of Ground with Partially Penetrating PVDs under the Free-Strain Condition. J. Eng. Mech.; 2021; 147, 04020148. [DOI: https://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0001884]

7. Deng, Y.B.; Xie, K.H.; Lu, M.M.; Tao, H.B.; Liu, G.B. Consolidation by prefabricated vertical drains considering the time dependent well resistance. Geotext. Geomembr.; 2013; 36, pp. 20-26. [DOI: https://dx.doi.org/10.1016/j.geotexmem.2012.10.003]

8. Fu, C.W.; Lei, G.H. A general solution for vertical-drain consolidation with impeded drainage boundaries. J. Cent. South Univ.; 2016; 23, pp. 934-943. [DOI: https://dx.doi.org/10.1007/s11771-016-3141-x]

9. Geng, X.Y. Non-linear consolidation of soil with vertical and horizontal drainage under time-dependent loading. Proceedings of the 2008 International Conference on Advanced Computer Theory and Engineering; Phuket, Thailand, 20–22 December 2008; pp. 800-804.

10. Indraratna, B.; Sathananthan, I.; Rujikiatkamjorn, C.; Balasubramaniam, A.S. Analytical and numerical modeling of soft soil stabilized by prefabricated vertical drains incorporating vacuum preloading. Int. J. Geomech.; 2005; 5, pp. 114-124. [DOI: https://dx.doi.org/10.1061/(ASCE)1532-3641(2005)5:2(114)]

11. Lei, G.H.; Zheng, Q.; Ng, C.W.W.; Chiu, A.C.F.; Xu, B. An analytical solution for consolidation with vertical drains under multi-ramp loading. Géotechnique; 2015; 65, pp. 531-547. [DOI: https://dx.doi.org/10.1680/geot.13.P.196]

12. Zhou, W.H.; Lok, T.M.H.; Zhao, L.S.; Mei, G.X.; Li, X.B. Analytical solutions to the axisymmetric consolidation of a multi-layer soil system under surcharge combined with vacuum preloading. Geotext. Geomembr.; 2017; 45, pp. 487-498. [DOI: https://dx.doi.org/10.1016/j.geotexmem.2017.06.003]

13. Zhou, W.H.; Zhao, L.S.; Lok, T.M.H.; Mei, G.X.; Li, X.B. Analytical Solutions to the Axisymmetrical Consolidation of Unsaturated Soils. J. Eng. Mech.; 2018; 144, 04017152. [DOI: https://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0001382]

14. Zhu, G.F.; Yin, J.H. Consolidation analysis of soil with vertical and horizontal drainage under ramp loading considering smear effects. Geotext. Geomembr.; 2004; 22, pp. 63-74. [DOI: https://dx.doi.org/10.1016/S0266-1144(03)00052-9]

15. Huang, Y.C.; Li, T.Y.; Fu, X.L.; Chen, Z. Analytical solutions for the consolidation of unsaturated foundation with prefabricated vertical drain. Int. J. Numer. Anal. Methods Geomech.; 2020; 44, pp. 2263-2282. [DOI: https://dx.doi.org/10.1002/nag.3127]

16. Fredlund, D.G.; Hasan, J.U. One-dimensional consolidation theory: Unsaturated soils. Can. Geotech. J.; 1979; 16, pp. 521-531. [DOI: https://dx.doi.org/10.1139/t79-058]

17. Fredlund, D.G.; Morgenstern, N.R. Constitutive relations for volume change in unsaturated soils. Can. Geotech. J.; 1976; 13, pp. 261-276. [DOI: https://dx.doi.org/10.1139/t76-029]

18. Fredlund, D.G.; Rahardjo, H. Soil Mechanics for Unsaturated Soils; Wiley: Hoboken, NJ, USA, 1993.

19. Qin, A.F.; Sun, D.A.; Tan, Y.W. Semi-analytical solution to one-dimensional consolidation in unsaturated soils. Appl. Math. Mech.; 2010; 31, pp. 215-226. [DOI: https://dx.doi.org/10.1007/s10483-010-0209-9]

20. Li, L.Z.; Qin, A.F.; Jiang, L.H.; Mei, G.X. Semianalytical Solution for One-Dimensional Consolidation in a Multilayered Unsaturated Soil System with Exponentially Time-Growing Permeable Boundary. J. Eng. Mech.; 2021; 147, 04021025. [DOI: https://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0001934]

21. Shan, Z.D.; Ling, D.S.; Ding, H.J. Exact solutions for one-dimensional consolidation of single-layer unsaturated soil. Int. J. Numer. Anal. Methods Geomech.; 2012; 36, pp. 708-722. [DOI: https://dx.doi.org/10.1002/nag.1026]

22. Shan, Z.D.; Ling, D.S.; Ding, H.J. Analytical solution for 1D consolidation of unsaturated soil with mixed boundary condition. J. Zhejiang Univ. -Sci. A; 2013; 14, pp. 61-70. [DOI: https://dx.doi.org/10.1631/jzus.A1200167]

23. Shan, Z.D.; Ling, D.S.; Ding, H.J. Analytical solution for the 1D consolidation of unsaturated multi-layered soil. Comput. Geotech.; 2014; 57, pp. 17-23. [DOI: https://dx.doi.org/10.1016/j.compgeo.2013.11.009]

24. Zhou, W.H.; Zhao, L.S. One-dimensional consolidation of unsaturated soil subjected to time-dependent loading with various initial and boundary conditions. Int. J. Geomech.; 2014; 14, pp. 291-301. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000314]

25. Zhou, W.H.; Zhao, L.S.; Li, X.B. A simple analytical solution to one-dimensional consolidation for unsaturated soils. Int. J. Numer. Anal. Methods Geomech.; 2014; 38, pp. 794-810. [DOI: https://dx.doi.org/10.1002/nag.2231]

26. Zhou, W.H.; Zhao, L.S.; Garg, A.; Yuen, K.V. Generalized analytical solution for the consolidation of unsaturated soil under partially permeable boundary conditions. Int. J. Geomech.; 2017; 17, 04017048. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000942]

27. Wang, L.; Sun, D.A.; Qin, A.F. General semi-analytical solutions to one-dimensional consolidation for unsaturated soils. Appl. Math. Mech. -Engl.; 2017; 38, pp. 831-850. [DOI: https://dx.doi.org/10.1007/s10483-017-2209-8]

28. Wang, L.; Sun, D.A.; Qin, A.F.; Xu, Y.F. Semi-analytical solution to one-dimensional consolidation for unsaturated soils with semi-permeable drainage boundary under time-dependent loading. Int. J. Numer. Anal. Methods Geomech.; 2017; 41, pp. 1636-1655. [DOI: https://dx.doi.org/10.1002/nag.2694]

29. Ho, L.; Fatahi, B. Analytical solution for the two-dimensional plane strain consolidation of an unsaturated soil stratum subjected to time-dependent loading. Comput. Geotech.; 2015; 67, pp. 1-16. [DOI: https://dx.doi.org/10.1016/j.compgeo.2015.02.011]

30. Wang, L.; Xu, Y.F.; Xia, X.H.; He, Y.L. Semi-analytical solutions of two-dimensional plane strain consolidation in unsaturated soils subjected to the lateral semipermeable drainage boundary. Int. J. Numer. Anal. Methods Geomech.; 2019; 43, pp. 2628-2651. [DOI: https://dx.doi.org/10.1002/nag.2986]

31. Wang, L.; Xu, Y.F.; Xia, X.H.; He, Y.L.; Li, T.Y. Semi-analytical solutions to two-dimensional plane strain consolidation for unsaturated soils under time-dependent loading. Comput. Geotech.; 2019; 109, pp. 144-165. [DOI: https://dx.doi.org/10.1016/j.compgeo.2019.01.002]

32. Li, L.Z.; Qin, A.F.; Jiang, L.H. Semi-analytical solution for one-dimensional consolidation of a two-layered soil system with unsaturated and saturated conditions. Int. J. Numer. Anal. Methods Geomech.; 2021; 45, pp. 2284-2300. [DOI: https://dx.doi.org/10.1002/nag.3266]

33. Qin, A.F.; Sun, D.A.; Yang, L.P.; Weng, Y.F. A semi-analytical solution to consolidation of unsaturated soils with the free drainage well. Comput. Geotech.; 2010; 37, pp. 867-875. [DOI: https://dx.doi.org/10.1016/j.compgeo.2010.07.006]

34. Zhou, W.H.; Tu, S. Unsaturated consolidation in a sand drain foundation by differential quadrature method. Procedia Earth Planet. Sci.; 2012; 5, pp. 52-57.

35. Ho, L.; Fatahi, B.; Khabbaz, H. Analytical solution to axisymmetric consolidation in unsaturated soils with linearly depth-dependent initial conditions. Comput. Geotech.; 2016; 74, pp. 102-121. [DOI: https://dx.doi.org/10.1016/j.compgeo.2015.12.019]

36. Ho, L.; Fatahi, B. Axisymmetric consolidation in unsaturated soil deposit subjected to time-dependent loadings. Int. J. Geomech.; 2017; 17, 04016046. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000686]

37. Ho, L.; Fatahi, B. Analytical solution to axisymmetric consolidation of unsaturated soil stratum under equal strain condition incorporating smear effects. Int. J. Numer. Anal. Methods Geomech.; 2018; 42, pp. 1890-1913. [DOI: https://dx.doi.org/10.1002/nag.2838]

38. Huang, Y.C.; Li, T.Y.; Fu, X.L. Consolidation of Unsaturated Drainage Well Foundation with Smear Effect under Time-Dependent Loading. KSCE J. Civ. Eng.; 2021; 25, pp. 768-781. [DOI: https://dx.doi.org/10.1007/s12205-021-0086-x]

39. Li, T.Y.; Qin, A.F.; Pei, Y.C.Q.; Sun, D.A.; Fu, X.L. Semianalytical solutions to the consolidation of drainage well foundations in unsaturated soils with radial semipermeable drainage boundary under time-dependent loading. Int. J. Geomech.; 2020; 20, 04020150. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0001767]