Three-step Predictor-Corrector Finite Element Schemes for Consolidation Equation

An accurate, stable, and eﬃcient three-step predictor-corrector time integration method is considered, for the ﬁrst time, to obtain numerical solution for the one-dimensional consolidation equation within a ﬁnite and spectral element framework. Theoretical order of accuracy and stability conditions are provided. The three-step predictor-corrector time integration method is third-order accurate and shows a larger stability region than the forward Euler method when applied to the one-dimensional consolidation equation. Furthermore, numerical results are in agreement with analytical solutions previously derived by the authors.


Introduction
Consolidation is a gradual process which involves, simultaneously, drainage, compression, and stress transfer [1]. In other words, consolidation is a quasistatic transient process in which there is coupling between the fluid flow and soil skeleton while effective stresses are updated [2]. A theory to explain the one-dimensional consolidation process, where the pore fluid flow and the skeleton deformation are constrained to take place in the vertical direction, was first proposed by Terzaghi [3,4]. Based on Terzaghi work, Biot [5] proposed a close analytical solution for the one-dimensional consolidation process in a finite length column of soil under constant loading.
In practice, external surface loads may be time-dependent. To analyze time-dependent loading, a graphical construction method was suggested by Terzaghi [4]. External loads can be smooth (harmonic, haversine, etc.) and piecewise smooth (step, ramp, triangular loading, etc.).
Several analytical solutions for the under different timedependent external loads have been presented in the specialized literature. For instance, Olson [6] presented charts for one-dimensional consolidation when the external load is the simple ramp loading. Razouki et al. [7] proposed an analytical solution for when the external load is the haversine cyclic loading without a rest period while Müthing et al. [8] achieved an analytical solution using Fourier harmonic analysis for the haversine cyclic loading with the rest period. In another study, Walker and Indraratna [9] suggested an analytical solution to three-multilayer consolidation problems based on the spectral-Galerkin method for when loads may vary in depth and time. Recently, a fourstep strategy was presented by Stickle and Pastor [10] based on the eigenfunction expansion in order to provide an analytical solution to the one-dimensional consolidation equation, which can be applied under general loading profiles. erefore, Stickle and Pastor method is not restricted to a particular type of loading. e present study compared the numerical results of the proposed method with the analytical solution provided by Stickle and Pastor. A large number of researchers have investigated the numerical solutions of the linear and nonlinear consolidation equations [11][12][13][14][15][16]. Several of the studies have adopted a finite element formulation for spatial discretization and a finite difference approach for time discretization. Huang and Griffiths [17] used the linear finite element method to solve coupled, uncoupled, and classical Terzaghi's [3] one-dimensional consolidation equation. In addition, Desai and Johnson [18] solve numerically the one-dimensional consolidation equation over a short time interval considering linear and cubic finite element formulations for spatial discretization. Feldkamp [19] discussed a semidiscrete formulation with regard to the finite element algorithm, utilizing quadratic Lagrange basis functions to handle the highly nonlinear PDE of one-dimensional large-strain consolidation. In a recent study by Chen et al. [20], backward Euler weak Galerkin finite element scheme was used to solve Biot's [5] consolidation. e aim of the present study was to provide a three-step predictor-corrector method for time discretization of onedimensional consolidation equation. e first step works as a predictor for the second step and the second step works as a predictor for the third step, which is considered as the final corrector step. For the suitable differentiable function F(t), the steps are defined as follows: As can be seen, equations (1)-(3) are equivalent to the third-order Taylor expansion [21,22]. erefore, the proposed method is third-order accurate for t. Equations (1)- (3) are usually accompanied by the Galerkin finite element method, known as the three-step Taylor-Galerkin method [23,24].
Done [25] was the first to introduce high-order timestepping schemes based on the second-and third-order Taylor expansion, namely, Taylor-Galerkin methods. In comparison with ordinary time-stepping methods, Taylor-Galerkin methods generate high phase accuracy with improved stability properties. However, the one-step Taylor-Galerkin schemes use too many terms in the third-order time derivative, especially for the nonlinear multidimensional equations.
Multistep Taylor-Galerkin methods can offer high accuracy for nonlinear multidimensional problems without computational complexity. For instance, Quartapelle [26] presented a two-step Taylor-Galerkin method based on fourth-order Taylor expansion, which can effectively handle nonlinear problems.
Jiang and Kawahara [21] are the pioneers of using the three-step Taylor-Galerkin method. In a study conducted by Kashiyama et al. [27], the abovementioned three steps, along with the finite element methods, were followed to handle computing shallow water flows. In addition, supercomputers were used to implement the related massively parallel finite element computations. Further, Kumar and Mehra [28] used Daubechies wavelets with this type of numerical integration to solve partial differential equations. Torabi and Hosseini [22] used the three-step approach and Legendre wavelets to solve the time-dependent equation.
In the present study, performance of two different spatial discretization techniques for the finite element method is compared. One of the methods employs linear interpolation functions and the other employs Lobatto interpolation polynomials as the shape functions. Both methods enjoy advantages and disadvantages.
Explicit schemes for time integration are very efficient as they might not require solving algebraic system of equations to advance in time. However, explicit integration schemes are conditionally stable. is is because the time step in an explicit integration scheme requires an upper bound in order to avoid round-off error accumulation that might lead to poor numerical results. On the contrary, implicit integration schemes might be endowed with larger time step than explicit ones. However, implicit integration schemes require solving an algebraic system of equations to advance in time, making them less efficient than the explicit ones. Classical explicit integration methods, like forward Euler, present too restrictive upper bounds for the time step in order to achieve stability. e use of efficient explicit integration schemes with improved accuracy and stability properties is a central task regarding numerical solutions of practical engineering problems. In this respect, a comparison was drawn between the proposed three-step method and explicit forward Euler solutions. e organization of this paper is as follows. In Section 2, the fundamental properties of the consolidation model are summarized. e proposed method is presented in Section 3. Stability and accuracy analysis is explained in Section 4.
Some numerical examples are presented in Section 5. Finally, Section 6 provides the conclusions of the study.

Basic Properties of the Consolidation Equation for Saturated Soils
Section 2 provides a summary of the consolidation equation.
More details can be found in [10]. We assume that a porous medium is fully saturated. Solid and fluid phases are intrinsically incompressible. Neither mass nor heat exchanges between solid and fluid phases are considered. Internal viscosity force within the fluid is neglected as compared with the viscous resistance opposed by the internal walls of the porous medium. Permeability is considered isotropic. Porous solid is considered isotropic and linear elastic under a small strain regime. All dynamic terms including convective terms are neglected. Only loading by external forces is considered and body forces are excluded.
Consider the one-dimensional column with a length of L. e load is applied at the upper boundary of the column which is perfectly drained. e motion of the solid skeleton within the column is constrained to take place along the vertical direction, defined by the z axis, which is considered normal to the top boundary and pointing downward as illustrated in Figure 1. e bottom of the column is considered impermeable and fixed. erefore, a standard PTIB (permeable top impermeable bottom) column is considered. Zero vertical displacement is assumed through the column at the initial time.
Under the abovementioned assumptions, one-dimensional consolidation equations can be expressed as where u is the solid vertical displacement, p is the pore fluid pressure, k is Darcy's permeability coefficient, and ρ is the intrinsic density of the fluid phase. Also, zz stands for partial differentiation with respect to the spatial coordinate z and zt means partial differentiation with respect to time. Appropriate initial and boundary conditions should be specified along with equations (4) and (5). Due to the above description of the problem, the boundary conditions for p and u are achieved as follows: where σ is the total Cauchy stress and f(t) is the external load. Also the initial conditions are considered as follows: By integrating equation (4) with respect to z and using boundary conditions the following relation between u and p can be obtained: Taking the partial derivative with respect to t in (8) and bearing (5) in mind, the well known diffusion like equation governing the pore pressure distribution is obtained as follows: where c ] � (k(λ + 2μ)/ρg) is the coefficient of consolidation in vertical direction while λ + 2μ is none other than the oedometer modulus. In order to solve (9), the following initial and boundary conditions are considered:

Analysis of the Proposed Method
Section 3 explains the main structure of the three-step predictor-corrector method, along with the use of two finite element method. (9) with the boundary and initial conditions (10) and (11). Assume that n ≥ 0 and Δt denotes the time step such that t n � nΔt, n � 0, 1, . . . N t . By using Taylor expansion, the value of function p(x, t) at the time t n+1 can be expressed as follows:

Time Discretization. Consider equation
where symbols p n and (zp/zt) n represent p(x, t n ) and ((zp/zt)(x, t n )), respectively. e forward Euler explicit method was obtained using Taylor's first-order expansion. In addition, time derivative in the given differential equations was approximated by Euler's formula as follows: erefore a semidiscrete equation was obtained: e proposed predictor-corrector method included three steps for time discretization, derived from applying a factorization process to the right side of equation (12) as follows: where the symbol I is the identity operator. Accordingly, the proposed predictor-corrector method was obtained using equation (15), employing a new notation as follows: It should be noted that p n+(1/3) , p n+ (1/2) , and p n+1 represent the computed solution at time level (t n + (Δt/3)), (t n + (Δt/2)), and (t n + Δt), respectively. As can be seen from equations (16)-(18), the first equation plays the predictor role for the second equation and the second one is used as a predictor for the third equation. Accordingly, the third-order accuracy in time without handling higher-order derivatives was achieved through equations (16)- (18) in time discretization [21].

Spatial Discretization.
e spatial derivatives of p(x, t) are approximated by finite element methods. is section describes the use of linear and spectral elements in the spatial discretization of equation (9). e coefficient of consolidation is assumed constant along the column. For more information on finite element methods, readers can refer to [29][30][31].

Discretization with Linear Elements.
In finite element methods, the solution domain is discretized into elementary units. For linear elements, the elements are first-order elements. erefore, the polynomial interpolation functions are first order. We divide the one-dimensional solution domain, 0 ≤ z ≤ L, into N E subdomains. e element size is constant and equal to h � (L/N E ). We define the global shape interpolation functions ϕ j (z) for j � 1, . . . , N E + 1. e jth function is supported by the interval [z j−1 , z j+1 ], takes the value of unity at the jth node, drops linearly to zero at the adjacent nodes, z j−1 and z j+1 , and remains zero outside the supporting interval. By introducing the nodal variables p j (t), the pore pressure distribution in equation (9) can be expanded as where P(t) and Φ(x) are determined as follows: e initial condition (11) gives P(0), and the first Dirichlet boundary condition, i.e., p(0, t) � 0 leads to Following the Galerkin method and considering equation (21), both sides of equation (9) were multiplied by each one of the global interpolation functions ϕ i (z) for i � 2, 3, . . . , N E + 1, and then integrated to make a weak form of the equation as follows: By applying the integration by parts on the right side of equation (22) and imposing the boundary condition (10), we can achieve the following equation: en, equation (19) is substituted into equation (23), leading to a new equation: Equation (24) can be written in matrix form using equation (21) as presented as follows: where M and D are the global mass and diffusion matrix, respectively. ese matrices are defined as follows: For example, the structure of the global mass matrix is as follows: 4 Mathematical Problems in Engineering Also, the ith-entry of the vector F(t) is defined as It should be noted that, by transferring the off-diagonal elements of the mass matrix (M) in each row to the corresponding diagonal element, numerical computations can be reduced. is simplification is known as mass lumping [31]. e differential equation (25) can be evaluated at any chosen time (t). e time derivative on the left-hand side is approximated by the three-step method, discussed in Section 3.1. Substituting equation (25) into equations (16)-(18) results in erefore, after applying equations (30)- (32), and imposing the initial condition, vector P n+1 can be specified in each time (t n ) for n � 0, 1, 2, . . . , N t .

Discretization with Spectral Elements.
High-order finite element methods use high-order polynomial expansions over the individual elements. e polynomial expansions themselves are defined by an appropriate number of element interpolation nodes. In this section, we use the notation N E and N G to show the total number of elements and unique global interpolation nodes, respectively.
In the spectral element method, the element interpolation nodes are distributed at the zeros of an appropriate family of orthogonal polynomials over the standard interval of the ξ axis, subject to the mandatory constraints that the first node is placed at ξ 1 � 1 and the last node is placed at ξ m+1 � −1, where m is the order of the polynomial expansion defined by m + 1 nodes. ese constraints ensure the C 0 continuity of the finite element expansion [31]. It should be noted that m may vary across the N E elements. erefore, to implement the spectral element method for the consolidation equation (9), we map the lth-element to the standard interval of the dimensionless ξ axis, [−1, 1], by the function where z (l) 1 is the first element end-node and z (l) 2 is the second element end-node. As ξ increases from −1 to 1, the physical point z is shifted from the first element end-node to the second element end-node. en, we approximate the pore pressure distribution over the lth-element with an mth-degree polynomial, P (l) m , expressed in the modal form where i�1 is a set of m + 1 bases functions, which are defined as follows: where Lo i−2 in the equation (35) is the (i − 2)-degree Lobatto polynomial. e Lobatto polynomials, Lo i , are the derivatives of the Legendre polynomials, L i , according to the relation Hence, they have similar properties. To ensure C 0 continuity of the solution, we require that the end-node values P (l) 1 and P (l) m+1 are shared by neighboring elements at the corresponding positions [31]. Also, we choose the element interior interpolation nodes, ξ j , for j � 2, . . . , m as the zeros of (m − 1)-degree Lobatto polynomial. e above bases are in local variable. It is possible to extend them into the global variable by Note that for each l � 1, . . . , N E − 1, ϕ (l) m+1 (z) and ϕ (l+1) 1 (z) are the two parts of a finite element function defined according to the three nodes [z (l) ]. By the continuous condition, the local degree of freedom Mathematical Problems in Engineering 5 corresponds to ϕ (l) m+1 (z), and ϕ (l+1) 1 (z) refers to the same global degree of freedom. So the global basis set is Now, we can approximate the pore pressure distribution over the spatial domain, [0, L], and implement the Galerkin projection by using the global basis, ϕ j (z) So, we follow the similar process discussed in the linear elements discretization, i.e, equations (19)- (25). For example, equation (24) can be expanded by We have to consider the domain integrals in (26) and (27) are assembled in terms of corresponding element integrals.
e lth-element diffusion and mass matrices are given by where h l � z (l) 2 − z (l) 1 is the element size and i, j � 1, . . . , m + 1. In other words, global diffusion and mass matrices are N G × N G matrices, that are obtained by assembly of matrices D (l) N E l�1 and M (l) N E l�1 , respectively. To implement the Dirichlet boundary condition, after an extended preliminary N G × N G system has been obtained for all nodes, the first constituent equation is discarded and replaced with the Dirichlet boundary condition. erefore, the unknown vector, P, contains the values of the solution at all nodes, except at the very first node where the Dirichlet boundary condition is specified. So, the spectral finite element method provides us with a system of N G − 1 ordinary differential equations for the time at all but the first node.

Stability and Accuracy Analysis
A numerical method is said to be stable when the numerical errors, e.g., due to round-off, are not amplied and the approximate solution remains bounded. is solution remains bounded when we fix the final time, t F � nΔt, and let n ⟶ ∞. erefore, a stability analysis applies a restriction to time step. e most popular technique for stability analysis is the Von Neumann method [32]. Jiang and Kawahara [21] employed this method to evaluate the stability condition of the three-step and two-step methods for unsteady incompressible flows. Done [25] used Von Neumann method for analyzing stability for different Taylor-Galerkin methods.
WWe present the stability analysis for the one-dimensional purely consolidation equation. After spatial discretization with linear shape functions, the following equations are achieved: where P n � [p n 1 , p n 2 , . . . , p n N E +1 ] is the vector of nodal values of unknowns and the superscript n represents the number of time steps. Also M is the lumped global mass matrix. e present study only explained the stability analysis of equation (41) and the stability of the other equations can be analyzed in the same way. By considering the structure of the lumped global mass (M) and diffusion (D) matrices and rearranging indices, equation (44) is obtained as follows: Equation (44) can be written as follows: where C r � (c ] Δt/h 2 ) is the diffusion number, and δ 2 is the standard notation for the difference operator such that e round-off error can be defined as follows: where p n j and N n j are the numerical and exact nodal values, respectively, which satisfy equation (45). Hence, substituting ε n j in equation (45) results in In addition, based on the Von Neumann stability analysis, error variation can be displayed by using Fourier series as follows: where g ≡ g(mπh, C r ) is the amplification factor and i � �� � −1 √ . Since the difference equation (48) is linear, it is sufficient to consider the growth of error in a typical term, which is where k � mπ is the wave number. Also, we can write after substituting ε n j and ε n+(1/3) j in equation (48), we have where the difference operator δ 2 can determined through equation (46). After dividing each term in equation (52) by g n e ik(jh) , the following equation is obtained: Since e ikh � cos(kh) + i sin(kh), we have e −ikh − 2 + e ikh � (−2 + 2 cos(kh)) � −2(1 − cos(kh)).

(54)
Substituting equation (54) in equation (53) results in Equation (55) holds true for j � 2, . . . , N E , therefore standing for the vector of nodal values (P n+ (1/3) ). e amplification factors for equations (42) and (43) are identified in a similar process. For example, for the second and third steps, the following equations are achieved: (56) us, the amplification factors for the second and third steps are achieved as follows: If ξ � kh, by combining equations (55)-(58), the following equation is obtained for the amplification factor g: For a present model problem, we shall say that a method is stable if there exists a constant K, independent of k, such that: Evidently, the following condition is achieved for stability based on Von Neumann analysis: g kh, C r ≤ 1 + K ′ Δt, ∀k.
(61) Accordingly, the necessary and sufficient condition for the error to remain bounded is [32]: g kh, C r ≤ 1. (62) Finally, the numerical stability condition for the threestep method using equation (59) is achieved as follows: erefore, we have Based on [31], the stability condition for the Euler explicit methods is C r ≤ 0.5. erefore, the proposed method is stable for a bigger time step (Δt).
To prove the third-order accuracy, we noticed the consolidation equation (9) is linear, homogeneous with constant coefficients. erefore, it has solutions of the following form: where ω is the temporal frequency. Such solutions are called plane waves. When the plane wave form is substituted into the consolidation equation, we obtain the following dispersion relation: Upon substituting this relation back in to the plane wave, the following solution for the consolidation equation is obtained: On the other hand, by making ξ ⟶ 0 in equation (59), the following equation is achieved: Substituting the corresponding values of C r and ξ in equation (68) results in is expression agrees with the first four terms of the power series of e − c v k 2 Δt ; therefore, the propose three-step method is third-order accurate.

Numerical Examples
is section discusses some numerical examples in the form of equation (9) with initial and boundary conditions (10) and (11). Error function is defined as maximum error (L ∞ ) as follows:

Mathematical Problems in Engineering
where p is the approximate solution, obtained from the numerical method, and z i N G i�1 denotes the total interpolation element nodes.
In order to assess the accuracy and robustness of the proposed method, the results were compared with the new semianalytical solutions derived by Stickle and Pastor [10]. eir method can obtain the solution for the one-dimensional consolidation equation under any loading profile (f(t)), whether the external load is smooth or piecewise smooth. In addition, pore pressure, i.e., p(z, t), along with the time derivative of the external loading (f ′ (t)) are expressed in Stickle and Pastor's method as where the coefficients f n (t) are the Fourier coefficients of f ′ (t), and p n (t) can be computed by where N n � c ] ((1 + 2n)π/2L) 2 is the coefficient of consolidation in the vertical direction. e abovementioned equation works very well for external loads. When the loading profile was piecewise smooth, it defined an adequate partition of the interval [0, t max ], where t max is the maximum time that the evolution of pore pressure should be computed. en, f(t) was expanded by continuous and Heaviside functions on each interval of the partition. e distribution of f(t) was effectively used to compute the integral term in equation (72). Finally, a finite number of terms were added in expression (71) in order to obtain the pore pressure distribution. e present study only considered the smooth loading profile, such as constant and cyclic haversine. Example 1. First example refers to the classical theory developed by Terzhagi [3,4] for the one-dimensional consolidation under constant loading. erefore, in this example, we have f ′ (t) � 0 for the whole consolidation process. By using equation (72), the coefficients p n (t) are given by e value of f(0) in (73) corresponds to the initial value of the load profile, which is held constant in the present case of analysis. Substituting (73) in (71) provides the analytical solution for the pore fluid pressure distribution, p(z, t). Table 1 indicates a comparison between the three-step solutions with the Euler explicit solutions using linear elements. e first two rows of this table present the results obtained for T ] � 0.1, where T ] is the dimensionless time factor and the final time is obtained through A comparison between the three-step solutions and Euler explicit solutions is presented in Table 2 using spectral elements. In the proposed method, 15 elements were considered in total. e pore pressure distribution with secondand fourth-order polynomials was approximated over the even and odd elements, respectively. In this case, 47 total interpolation nodes, i.e., N G � 47 existed. Numerical results demonstrated that the three-step method leads to higher stability criteria compared to Euler explicit method when using spectral elements. Table 3 compares the linear and spectral elements with the three-step method for the fix Δt � 0.0001, which indicates that accuracy improves after using spectral elements compared to linear elements. However, three-step method with spectral elements is stable for the smaller time step.
Further, the three-step and Euler explicit solution were compared using linear elements (Figures 2-3). e threestep solution using linear elements, generated every 100 time steps, are shown in Figure 2 while Figure 3 presents the Euler explicit solution using linear elements, generated every 15 time steps. By comparing these two figures, we can conclude that the three-step solution remains bounded as time increases whereas Euler explicit solution tend to infinity. Figures 4 and 5 illustrate a comparison between threestep solutions and Stickle-Pastor solutions [10]. Linear and spectral elements are used in Figures 4 and 5, respectively. Accordingly, a visual comparison between linear and spectral elements is obtained, showing that the accuracy of the spectral elements is greater than that of linear elements.
Example 2. In this example, we consider equation (9) with a cyclic haversine loading.
is case has been studied by Razouki and coworkers [7] in 2013 proposing an analytical expression to describe the response under this type of cyclic loading. Following Huang [33], the profile f(t) for a haversine loading might be given as where t c � 0.1(L 2 /c ] ). is type of loading profile is smooth. Taking time derivative of (75), f ′ (t) has the following expression for a haversine loading: Bearing in mind that f(0) � 0 in this case of study, substituting (76) in (72) and computing the integral, the following expression is obtained for the coefficients p n (t): Substituting (77) in (71) provides the analytical solution for the pore fluid pressure distribution p(z, t). Table 4 presents a comparison between the three-step solutions with Euler explicit solutions using linear elements. e first two rows of the table indicates the results obtained for the final time (t c /2), revealing that the Euler explicit method is stable for the smaller time step (Δt � 0.0001) while the three-step method is stable for (Δt � 0.004). Other parameters can be similarly analyzed.
In Table 5, spectral elements are used to compare the three-step and Euler explicit solutions. Similar to the first example, N E � 15 and N G � 47 were selected. Table 6 presents a comparison between linear and spectral elements with the three-step method for fix Δt � 0.0001. e table shows that spectral elements provide better accuracy than linear elements. However, three-step method with spectral elements is stable for the smaller time step. Figures 6 and 7 display three-step and Euler explicit solutions using linear elements, generated every 5 time steps, indicating that the three-step solutions remain bounded without oscillations whereas Euler explicit solutions have oscillations. Figures 8 and 9 draw a comparison between three-step solutions and Stickle-Pastor solutions [10], where linear elements are used in Figure 8 and spectral elements in Figure 9. Depth/column length

Conclusions
In the present study, a new numerical method was proposed in order to solve one-dimensional consolidation equation. is method included three steps for time discretization and presented third-order accuracy. Linear and spectral finite element methods were used for spatial discretization. Further, the theoretical aspect of the stability condition was discussed.
e proposed method was compared with the Euler explicit and semianalytical methods, presented by Stickle and Pastor [10]. e numerical results suggested that the three-step solutions with spectral elements have higher accuracy than the three-step solutions with linear elements. In addition, the proposed method shows more stability than Euler explicit solutions and can be easily used for highdimensional consolidation equations.

Data Availability
e calculated data used to support the findings of this study are included within the article. Some basic MATLAB finite element software library can be downloaded from the website http://dehesa.freeshell.org/FSELIB. Other pieces of information are available from the corresponding author upon request. Mathematical Problems in Engineering 13