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

: 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 veriﬁed 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 signiﬁcantly hinder the dissipation process of excess pore pressures, and different external loading types will lead to various dissipation characteristics (i.e., peak values).


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 2 of 12 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 semianalytical 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.

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.
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. rw, rs and re are the largest radii of the drain, smeared zone and the influenced area, respectively. The thickness of the layered foundation is H. ka d , kw d ka s , kw s , ka and kw are the permeability coefficients of air and water phases in the drain, smeared zone and influenced areas, respectively.

Basic Assumptions
In this paper, the solutions are derived based on the following assumptions: 1. The soil layer is homogeneous. 2. Solid particles and the water phase are incompressible. 3. Air and water phases can only flow along the radial direction. 4. The air and water phases are continuous and independent. 5. Permeability coefficients related to air and water phases are constant.

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:

Basic Assumptions
In this paper, the solutions are derived based on the following assumptions: 1.
The soil layer is homogeneous.

2.
Solid particles and the water phase are incompressible.

3.
Air and water phases can only flow along the radial direction. 4.
The air and water phases are continuous and independent. 5.
Permeability coefficients related to air and water phases are constant.

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: (•) r and (•) rr are the first and second derivatives of a function upon r, and (•) 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 a v and C w v are the consolidation coefficients of air and water phases and u a and u w , respectively, are average excess pore-air and water pressures, defined as: where, u s = u a s u w s T , u a s denotes excess pore-air pressure in the smeared area and u w s is the excess pore-water pressure in the smeared area. The consolidation parameters are expressed as: where, m a 1 , m a 2 , m w 1 and m w 2 , respectively, are the variations of net stress (q 0 -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 0 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 2 , and γ w is the unit weight of water, and γ w = 9.8 kN/m 3 .

Initial Conditions
u 0 a and u 0 w , respectively, are the initial excess pore-air and pore-water pressures under the applied loading at t = 0, and u a d and u w d are excess pore-air and pore-water pressures in the drain, respectively.

Boundary Conditions
Continuation conditions at r = r w : where, It is noted that the excess pore pressures u a d and u w d , 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 : Boundary conditions at r = r e : Boundary conditions at z = 0: Boundary conditions z = H: Appl. Sci. 2021, 11, 9285 5 of 12

Analytical Solutions
Firstly, Equation (1) needs to be adopted to transform u a and u w into u a and u w respectively, to introduce the equal-strain theorem into the solution. Rearranging Equation (1) produces: 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]: 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]: where, F = F a 0 0 F w . F a and F w are the intermediate variables concerning smear effects: α 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 d in Equation (14) is the key to obtaining the final solutions of u. Therefore, substituting Equation (12) into Equation (6) affords: Combining Equations (14) and (17) produces the following equation:  (14) and (18) provides the exact expression of u d : where, λ = 2 r 2 e C −1 C v r F −1 , and λ 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: where, K = 2k+1 2 π, k = 1, 2, 3 . . . By substituting Equations (19) into (20), one can obtain: By applying the orthogonality of trigonometric functions, Equation (21) can be converted into: where, It is appointed that the expression of U d (t), 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 (t) according to the decoupling method proposed by , and the constant variation method by Huang et al. [38], is: C 1 and C 2 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)): where, Based on the consolidation theory by Fredlund and Hasan [16], the average consolidation degree can be expressed as: where, m s 1 = m a 1 + m w 1 , m s 2 = m a 2 + m w 2 ; q f is the final value of loading when time approaches to infinity.

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 0 is adopted in current solutions to perform the verification work. Practical values of variables are as follows: q 0 = 100 kPa, u 0 a = 20 kPa, u 0 w = 40 kPa, k w = 10 −10 m/s, α a = α w = 1; r w = r s = 0.2 m, r e = 1.8 m, H = 10 m, u atm = 101.3 kPa, m a 1 = −2 × 10 −4 kPa; m w 1 = −5 × 10 −5 kPa, S r = 80%, n 0 = 50%, m w 2 = −2 × 10 −4 kPa, m a 2 = −1 × 10 −4 kPa. 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.

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 hydraulicrelated 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: Exponential loading: 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 × 10 −2 , 5 × 10 −3 , 5 × 10 −4 and 5 × 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 1 < t < 10 6 ) 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 × 10 −2 , 5 × 10 −3 , 5 × 10 −4 , 5 × 10 −5 and 5 × 10 −6 s −1 . 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 × 10 −2 , 5 × 10 −3 , 5 × 10 −4 and 5 × 10 −5 kPa/s, b = 5 × 10 −2 , 5 × 10 −3 , 5 × 10 −4 , 5 × 10 −5 and 5 × 10 −6 and s −1 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 × 10 −3 kPa/s and b > 5 × 10 −4 s −1 , 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 × 10 −2 kPa/s might be more reliable in practical engineering, the above-mentioned conclusions in Figures 3 and 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.

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.; writingreview and editing, Y.C., Q.X. and X.L. All authors have read and agreed to the published version of the manuscript.