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

Authors: Mina Torabi; Manuel Pastor; Miguel Martín Stickle
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

An accurate, stable, and efficient three-step predictor-corrector time integration method is considered, for the first time, to obtain numerical solution for the one-dimensional consolidation equation within a finite 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.


---

## Body

## 1. 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 time-dependent 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 four-step 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. Therefore, Stickle and Pastor method is not restricted to a particular type of loading. The 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–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.The aim of the present study was to provide a three-step predictor-corrector method for time discretization of one-dimensional consolidation equation. The 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 functionFt, the steps are defined as follows:(1)Ft+Δt3=Ft+Δt3∂F∂tt,(2)Ft+Δt2=Ft+Δt2∂F∂tt+Δt3,(3)Ft+Δt=Ft+Δt∂F∂tt+Δt2.As can be seen, equations (1)–(3) are equivalent to the third-order Taylor expansion [21, 22]. Therefore, 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 time-stepping 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. This 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. The 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.The organization of this paper is as follows. In Section2, the fundamental properties of the consolidation model are summarized. The 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.


## 2. Basic Properties of the Consolidation Equation for Saturated Soils

Section2 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 ofL. The load is applied at the upper boundary of the column which is perfectly drained. The 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. The bottom of the column is considered impermeable and fixed. Therefore, a standard PTIB (permeable top impermeable bottom) column is considered. Zero vertical displacement is assumed through the column at the initial time.Figure 1
Geometry and boundary conditions of the investigated problem.Under the abovementioned assumptions, one-dimensional consolidation equations can be expressed as(4)λ+2μ∂2u∂z2z,t−∂p∂zz,t=0,(5)∂2u∂t∂zz,t−kρg∂2p∂z2z,t=0,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, ∂z stands for partial differentiation with respect to the spatial coordinate z and ∂t 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:(6)p0,t=0,∂p∂zL,t=0,σ0,t=−ft,uL,t=0.where σ is the total Cauchy stress and ft is the external load. Also the initial conditions are considered as follows:(7)pz,0=f0,uz,0=0.By integrating equation (4) with respect to z and using boundary conditions the following relation between u and p can be obtained:(8)λ+2μ∂u∂zz,t−pz,t=−ft.Taking the partial derivative with respect tot in (8) and bearing (5) in mind, the well known diffusion like equation governing the pore pressure distribution is obtained as follows:(9)∂p∂tz,t=cν∂2p∂z2z,t+f′t,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:(10)boundary  conditions,p0,t=0∂p∂zL,t=0.(11)initial condition,pz,0=f0.


## 3. Analysis of the Proposed Method

Section3 explains the main structure of the three-step predictor-corrector method, along with the use of two finite element method.


### 3.1. Time Discretization

Consider equation (9) with the boundary and initial conditions (10) and (11). Assume that n≥0 and Δt denotes the time step such that tn=nΔt,n=0,1,…Nt. By using Taylor expansion, the value of function px,t at the time tn+1 can be expressed as follows:(12)pn+1=pn+Δt∂p∂tn+Δt22∂2p∂t2n+Δt36∂3p∂t3n+OΔt4,where symbols pn and ∂p/∂tn represent px,tn and ∂p/∂tx,tn, respectively.The 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:(13)∂p∂tn=px,tn+1−px,tnΔt+OΔt.Therefore a semidiscrete equation was obtained:(14)pn+1=pn+Δt∂p∂tn.The 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:(15)ℐ+Δt∂∂tℐ+Δt2∂∂tℐ+Δt3∂∂tpn=pn+Δt∂∂tpn+Δt2∂∂tpn+Δt3∂p∂tn,where the symbol ℐ is the identity operator. Accordingly, the proposed predictor-corrector method was obtained using equation (15), employing a new notation as follows:(16)pn+1/3=pn+Δt3∂p∂tn,(17)pn+1/2=pn+Δt2∂p∂tn+1/3,(18)pn+1=pn+Δt∂p∂tn+1/2.It should be noted thatpn+1/3, pn+1/2, and pn+1 represent the computed solution at time level tn+Δt/3, tn+Δt/2, and (tn+Δ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].


### 3.2. Spatial Discretization

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


#### 3.2.1. 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. Therefore, the polynomial interpolation functions are first order. We divide the one-dimensional solution domain,0≤z≤L, into NE subdomains. The element size is constant and equal to h=L/NE. We define the global shape interpolation functions ϕjz for j=1,…,NE+1. The jth function is supported by the interval zj−1,zj+1, takes the value of unity at the jth node, drops linearly to zero at the adjacent nodes, zj−1 and zj+1, and remains zero outside the supporting interval. By introducing the nodal variables p^jt, the pore pressure distribution in equation (9) can be expanded as(19)pz,t≃∑j=1NE+1p^jtϕjz=P^t⋅Φx,where P^t and Φx are determined as follows:(20)P^t=p^1t,p^2t,…,p^NE+1t,Φx=ϕ1z,ϕ2z,…,ϕNE+1z.The initial condition (11) gives P^0, and the first Dirichlet boundary condition, i.e., p0,t=0 leads to(21)p0,t=p^1t=0.Following the Galerkin method and considering equation (21), both sides of equation (9) were multiplied by each one of the global interpolation functions ϕiz for i=2,3,…,NE+1, and then integrated to make a weak form of the equation as follows:(22)∫0Lϕiz∂p∂tz,tdz=cν∫0Lϕiz∂2p∂z2dz+f′t∫0Lϕizdz,i=2,3,…,NE+1.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:(23)∫0Lϕiz∂p∂tz,tdz=−cν∫0Ldϕidzdpdzdz+f′t∫0Lϕizdz,i=2,3,…,NE+1.Then, equation (19) is substituted into equation (23), leading to a new equation:(24)∑j=1NE+1∫0Lϕjzϕizdp^jdttdz=−cν∑j=1NE+1∫0Lp^jtdϕjdzdϕidzdz+f′t∫0Lϕizdz,i=2,3,…,NE+1.Equation (24) can be written in matrix form using equation (21) as presented as follows:(25)M⋅dP^dtt=−cνD⋅P^t+Ft,where M and D are the global mass and diffusion matrix, respectively. These matrices are defined as follows:(26)Mi,j=∫0Lϕiϕjdz,(27)Di,j=∫0Ldϕidzdϕjdzdz.For example, the structure of the global mass matrix is as follows:(28)M=13h16h0…00016h23h16h…000016h23h…000⋮⋮⋮⋱⋮⋮⋮000…23h16h0000…16h23h16h000…016h13h.Also, theith-entry of the vector Ft is defined as(29)Fit=f′t∫0Lϕizdz,i=2,3,…,NE+1.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. This simplification is known as mass lumping [31].The differential equation (25) can be evaluated at any chosen time (t). The 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(30)P^n+1/3=P^n+Δt3M−1−cνD⋅P^n+Fn.(31)P^n+1/2=P^n+Δt2M−1−cνD⋅P^n+1/3+Fn+1/3.(32)P^n+1=P^n+ΔtM−1−cνD⋅P^n+1/2+Fn+1/2.Therefore, after applying equations (30)–(32), and imposing the initial condition, vector P^n+1 can be specified in each time (tn) for n=0,1,2,…,Nt.


#### 3.2.2. Discretization with Spectral Elements

High-order finite element methods use high-order polynomial expansions over the individual elements. The polynomial expansions themselves are defined by an appropriate number of element interpolation nodes. In this section, we use the notationNE and NG 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. These constraints ensure the C0 continuity of the finite element expansion [31]. It should be noted that m may vary across the NE elements.Therefore, 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(33)zξ=12z2l+z1l+12z2l−z1l,where z1l is the first element end-node and z2l 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.Then, we approximate the pore pressure distribution over thelth-element with an mth-degree polynomial, Pml, expressed in the modal form(34)pξ,t≃Pmlξ,t=∑i=1m+1p^iltζiξ,where ζiξi=1m+1 is a set of m+1 bases functions, which are defined as follows:(35)ζ1ξ=121−ξ,ζm+1ξ=121+ξ,ζiξ=1−ξ2Loi−2ξ,i=2,…,m,where Loi−2 in the equation (35) is the i−2-degree Lobatto polynomial. The Lobatto polynomials, Loi, are the derivatives of the Legendre polynomials, Li, according to the relation(36)Loiξ≡Li+1′ξ.Hence, they have similar properties. To ensureC0 continuity of the solution, we require that the end-node values P^1l and P^m+1l 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.The above bases are in local variable. It is possible to extend them into the global variable by(37)ϕilz=ζiξ,ifz∈z1l,z2l,0,ifz∉z1l,z2l.Note that for eachl=1,…,NE−1, ϕm+1lz and ϕ1l+1z are the two parts of a finite element function defined according to the three nodes z1l,z2l,z2l+1. By the continuous condition, the local degree of freedom corresponds to ϕm+1lz, and ϕ1l+1z refers to the same global degree of freedom. So the global basis set is(38)V=ϕjzj=1mNE+1=ϕm+1lz+ϕ1l+1zl=1NE−1∪ϕilz,i=2,…,ml=1NE∪ϕ11z,ϕm+1NEz.Now, we can approximate the pore pressure distribution over the spatial domain,0,L, and implement the Galerkin projection by using the global basis, ϕjzj=1mNE+1. 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(39)∑j=1m+1∫z1lz2lϕjlzϕilzdp^jldttdz=−cν∑j=1m+1∫z1lz2lp^jltdϕjldzdϕildzdz+f′t∫z1lz2lϕilzdz,i=1,2,…,m+1,l=1,2,…,NE.We have to consider the domain integrals in (26) and (27) are assembled in terms of corresponding element integrals. The lth-element diffusion and mass matrices are given by(40)Di,jl=2hl∫−11dζidξdζjdξdξ,Mi,jl=hl2∫−11ζiζjdξ,where hl=z2l−z1l is the element size and i,j=1,…,m+1. In other words, global diffusion and mass matrices are NG×NG matrices, that are obtained by assembly of matrices Dll=1NE and Mll=1NE, respectively.To implement the Dirichlet boundary condition, after an extended preliminaryNG×NG system has been obtained for all nodes, the first constituent equation is discarded and replaced with the Dirichlet boundary condition. Therefore, 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 NG−1 ordinary differential equations for the time at all but the first node.Finally, the three-step equations, explained in equations (16)–(18), were used for time discretization, and equations (30)–(32) were similarly achieved.


## 3.1. Time Discretization

Consider equation (9) with the boundary and initial conditions (10) and (11). Assume that n≥0 and Δt denotes the time step such that tn=nΔt,n=0,1,…Nt. By using Taylor expansion, the value of function px,t at the time tn+1 can be expressed as follows:(12)pn+1=pn+Δt∂p∂tn+Δt22∂2p∂t2n+Δt36∂3p∂t3n+OΔt4,where symbols pn and ∂p/∂tn represent px,tn and ∂p/∂tx,tn, respectively.The 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:(13)∂p∂tn=px,tn+1−px,tnΔt+OΔt.Therefore a semidiscrete equation was obtained:(14)pn+1=pn+Δt∂p∂tn.The 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:(15)ℐ+Δt∂∂tℐ+Δt2∂∂tℐ+Δt3∂∂tpn=pn+Δt∂∂tpn+Δt2∂∂tpn+Δt3∂p∂tn,where the symbol ℐ is the identity operator. Accordingly, the proposed predictor-corrector method was obtained using equation (15), employing a new notation as follows:(16)pn+1/3=pn+Δt3∂p∂tn,(17)pn+1/2=pn+Δt2∂p∂tn+1/3,(18)pn+1=pn+Δt∂p∂tn+1/2.It should be noted thatpn+1/3, pn+1/2, and pn+1 represent the computed solution at time level tn+Δt/3, tn+Δt/2, and (tn+Δ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].


## 3.2. Spatial Discretization

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


### 3.2.1. 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. Therefore, the polynomial interpolation functions are first order. We divide the one-dimensional solution domain,0≤z≤L, into NE subdomains. The element size is constant and equal to h=L/NE. We define the global shape interpolation functions ϕjz for j=1,…,NE+1. The jth function is supported by the interval zj−1,zj+1, takes the value of unity at the jth node, drops linearly to zero at the adjacent nodes, zj−1 and zj+1, and remains zero outside the supporting interval. By introducing the nodal variables p^jt, the pore pressure distribution in equation (9) can be expanded as(19)pz,t≃∑j=1NE+1p^jtϕjz=P^t⋅Φx,where P^t and Φx are determined as follows:(20)P^t=p^1t,p^2t,…,p^NE+1t,Φx=ϕ1z,ϕ2z,…,ϕNE+1z.The initial condition (11) gives P^0, and the first Dirichlet boundary condition, i.e., p0,t=0 leads to(21)p0,t=p^1t=0.Following the Galerkin method and considering equation (21), both sides of equation (9) were multiplied by each one of the global interpolation functions ϕiz for i=2,3,…,NE+1, and then integrated to make a weak form of the equation as follows:(22)∫0Lϕiz∂p∂tz,tdz=cν∫0Lϕiz∂2p∂z2dz+f′t∫0Lϕizdz,i=2,3,…,NE+1.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:(23)∫0Lϕiz∂p∂tz,tdz=−cν∫0Ldϕidzdpdzdz+f′t∫0Lϕizdz,i=2,3,…,NE+1.Then, equation (19) is substituted into equation (23), leading to a new equation:(24)∑j=1NE+1∫0Lϕjzϕizdp^jdttdz=−cν∑j=1NE+1∫0Lp^jtdϕjdzdϕidzdz+f′t∫0Lϕizdz,i=2,3,…,NE+1.Equation (24) can be written in matrix form using equation (21) as presented as follows:(25)M⋅dP^dtt=−cνD⋅P^t+Ft,where M and D are the global mass and diffusion matrix, respectively. These matrices are defined as follows:(26)Mi,j=∫0Lϕiϕjdz,(27)Di,j=∫0Ldϕidzdϕjdzdz.For example, the structure of the global mass matrix is as follows:(28)M=13h16h0…00016h23h16h…000016h23h…000⋮⋮⋮⋱⋮⋮⋮000…23h16h0000…16h23h16h000…016h13h.Also, theith-entry of the vector Ft is defined as(29)Fit=f′t∫0Lϕizdz,i=2,3,…,NE+1.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. This simplification is known as mass lumping [31].The differential equation (25) can be evaluated at any chosen time (t). The 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(30)P^n+1/3=P^n+Δt3M−1−cνD⋅P^n+Fn.(31)P^n+1/2=P^n+Δt2M−1−cνD⋅P^n+1/3+Fn+1/3.(32)P^n+1=P^n+ΔtM−1−cνD⋅P^n+1/2+Fn+1/2.Therefore, after applying equations (30)–(32), and imposing the initial condition, vector P^n+1 can be specified in each time (tn) for n=0,1,2,…,Nt.


### 3.2.2. Discretization with Spectral Elements

High-order finite element methods use high-order polynomial expansions over the individual elements. The polynomial expansions themselves are defined by an appropriate number of element interpolation nodes. In this section, we use the notationNE and NG 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. These constraints ensure the C0 continuity of the finite element expansion [31]. It should be noted that m may vary across the NE elements.Therefore, 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(33)zξ=12z2l+z1l+12z2l−z1l,where z1l is the first element end-node and z2l 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.Then, we approximate the pore pressure distribution over thelth-element with an mth-degree polynomial, Pml, expressed in the modal form(34)pξ,t≃Pmlξ,t=∑i=1m+1p^iltζiξ,where ζiξi=1m+1 is a set of m+1 bases functions, which are defined as follows:(35)ζ1ξ=121−ξ,ζm+1ξ=121+ξ,ζiξ=1−ξ2Loi−2ξ,i=2,…,m,where Loi−2 in the equation (35) is the i−2-degree Lobatto polynomial. The Lobatto polynomials, Loi, are the derivatives of the Legendre polynomials, Li, according to the relation(36)Loiξ≡Li+1′ξ.Hence, they have similar properties. To ensureC0 continuity of the solution, we require that the end-node values P^1l and P^m+1l 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.The above bases are in local variable. It is possible to extend them into the global variable by(37)ϕilz=ζiξ,ifz∈z1l,z2l,0,ifz∉z1l,z2l.Note that for eachl=1,…,NE−1, ϕm+1lz and ϕ1l+1z are the two parts of a finite element function defined according to the three nodes z1l,z2l,z2l+1. By the continuous condition, the local degree of freedom corresponds to ϕm+1lz, and ϕ1l+1z refers to the same global degree of freedom. So the global basis set is(38)V=ϕjzj=1mNE+1=ϕm+1lz+ϕ1l+1zl=1NE−1∪ϕilz,i=2,…,ml=1NE∪ϕ11z,ϕm+1NEz.Now, we can approximate the pore pressure distribution over the spatial domain,0,L, and implement the Galerkin projection by using the global basis, ϕjzj=1mNE+1. 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(39)∑j=1m+1∫z1lz2lϕjlzϕilzdp^jldttdz=−cν∑j=1m+1∫z1lz2lp^jltdϕjldzdϕildzdz+f′t∫z1lz2lϕilzdz,i=1,2,…,m+1,l=1,2,…,NE.We have to consider the domain integrals in (26) and (27) are assembled in terms of corresponding element integrals. The lth-element diffusion and mass matrices are given by(40)Di,jl=2hl∫−11dζidξdζjdξdξ,Mi,jl=hl2∫−11ζiζjdξ,where hl=z2l−z1l is the element size and i,j=1,…,m+1. In other words, global diffusion and mass matrices are NG×NG matrices, that are obtained by assembly of matrices Dll=1NE and Mll=1NE, respectively.To implement the Dirichlet boundary condition, after an extended preliminaryNG×NG system has been obtained for all nodes, the first constituent equation is discarded and replaced with the Dirichlet boundary condition. Therefore, 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 NG−1 ordinary differential equations for the time at all but the first node.Finally, the three-step equations, explained in equations (16)–(18), were used for time discretization, and equations (30)–(32) were similarly achieved.


## 3.2.1. 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. Therefore, the polynomial interpolation functions are first order. We divide the one-dimensional solution domain,0≤z≤L, into NE subdomains. The element size is constant and equal to h=L/NE. We define the global shape interpolation functions ϕjz for j=1,…,NE+1. The jth function is supported by the interval zj−1,zj+1, takes the value of unity at the jth node, drops linearly to zero at the adjacent nodes, zj−1 and zj+1, and remains zero outside the supporting interval. By introducing the nodal variables p^jt, the pore pressure distribution in equation (9) can be expanded as(19)pz,t≃∑j=1NE+1p^jtϕjz=P^t⋅Φx,where P^t and Φx are determined as follows:(20)P^t=p^1t,p^2t,…,p^NE+1t,Φx=ϕ1z,ϕ2z,…,ϕNE+1z.The initial condition (11) gives P^0, and the first Dirichlet boundary condition, i.e., p0,t=0 leads to(21)p0,t=p^1t=0.Following the Galerkin method and considering equation (21), both sides of equation (9) were multiplied by each one of the global interpolation functions ϕiz for i=2,3,…,NE+1, and then integrated to make a weak form of the equation as follows:(22)∫0Lϕiz∂p∂tz,tdz=cν∫0Lϕiz∂2p∂z2dz+f′t∫0Lϕizdz,i=2,3,…,NE+1.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:(23)∫0Lϕiz∂p∂tz,tdz=−cν∫0Ldϕidzdpdzdz+f′t∫0Lϕizdz,i=2,3,…,NE+1.Then, equation (19) is substituted into equation (23), leading to a new equation:(24)∑j=1NE+1∫0Lϕjzϕizdp^jdttdz=−cν∑j=1NE+1∫0Lp^jtdϕjdzdϕidzdz+f′t∫0Lϕizdz,i=2,3,…,NE+1.Equation (24) can be written in matrix form using equation (21) as presented as follows:(25)M⋅dP^dtt=−cνD⋅P^t+Ft,where M and D are the global mass and diffusion matrix, respectively. These matrices are defined as follows:(26)Mi,j=∫0Lϕiϕjdz,(27)Di,j=∫0Ldϕidzdϕjdzdz.For example, the structure of the global mass matrix is as follows:(28)M=13h16h0…00016h23h16h…000016h23h…000⋮⋮⋮⋱⋮⋮⋮000…23h16h0000…16h23h16h000…016h13h.Also, theith-entry of the vector Ft is defined as(29)Fit=f′t∫0Lϕizdz,i=2,3,…,NE+1.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. This simplification is known as mass lumping [31].The differential equation (25) can be evaluated at any chosen time (t). The 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(30)P^n+1/3=P^n+Δt3M−1−cνD⋅P^n+Fn.(31)P^n+1/2=P^n+Δt2M−1−cνD⋅P^n+1/3+Fn+1/3.(32)P^n+1=P^n+ΔtM−1−cνD⋅P^n+1/2+Fn+1/2.Therefore, after applying equations (30)–(32), and imposing the initial condition, vector P^n+1 can be specified in each time (tn) for n=0,1,2,…,Nt.


## 3.2.2. Discretization with Spectral Elements

High-order finite element methods use high-order polynomial expansions over the individual elements. The polynomial expansions themselves are defined by an appropriate number of element interpolation nodes. In this section, we use the notationNE and NG 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. These constraints ensure the C0 continuity of the finite element expansion [31]. It should be noted that m may vary across the NE elements.Therefore, 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(33)zξ=12z2l+z1l+12z2l−z1l,where z1l is the first element end-node and z2l 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.Then, we approximate the pore pressure distribution over thelth-element with an mth-degree polynomial, Pml, expressed in the modal form(34)pξ,t≃Pmlξ,t=∑i=1m+1p^iltζiξ,where ζiξi=1m+1 is a set of m+1 bases functions, which are defined as follows:(35)ζ1ξ=121−ξ,ζm+1ξ=121+ξ,ζiξ=1−ξ2Loi−2ξ,i=2,…,m,where Loi−2 in the equation (35) is the i−2-degree Lobatto polynomial. The Lobatto polynomials, Loi, are the derivatives of the Legendre polynomials, Li, according to the relation(36)Loiξ≡Li+1′ξ.Hence, they have similar properties. To ensureC0 continuity of the solution, we require that the end-node values P^1l and P^m+1l 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.The above bases are in local variable. It is possible to extend them into the global variable by(37)ϕilz=ζiξ,ifz∈z1l,z2l,0,ifz∉z1l,z2l.Note that for eachl=1,…,NE−1, ϕm+1lz and ϕ1l+1z are the two parts of a finite element function defined according to the three nodes z1l,z2l,z2l+1. By the continuous condition, the local degree of freedom corresponds to ϕm+1lz, and ϕ1l+1z refers to the same global degree of freedom. So the global basis set is(38)V=ϕjzj=1mNE+1=ϕm+1lz+ϕ1l+1zl=1NE−1∪ϕilz,i=2,…,ml=1NE∪ϕ11z,ϕm+1NEz.Now, we can approximate the pore pressure distribution over the spatial domain,0,L, and implement the Galerkin projection by using the global basis, ϕjzj=1mNE+1. 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(39)∑j=1m+1∫z1lz2lϕjlzϕilzdp^jldttdz=−cν∑j=1m+1∫z1lz2lp^jltdϕjldzdϕildzdz+f′t∫z1lz2lϕilzdz,i=1,2,…,m+1,l=1,2,…,NE.We have to consider the domain integrals in (26) and (27) are assembled in terms of corresponding element integrals. The lth-element diffusion and mass matrices are given by(40)Di,jl=2hl∫−11dζidξdζjdξdξ,Mi,jl=hl2∫−11ζiζjdξ,where hl=z2l−z1l is the element size and i,j=1,…,m+1. In other words, global diffusion and mass matrices are NG×NG matrices, that are obtained by assembly of matrices Dll=1NE and Mll=1NE, respectively.To implement the Dirichlet boundary condition, after an extended preliminaryNG×NG system has been obtained for all nodes, the first constituent equation is discarded and replaced with the Dirichlet boundary condition. Therefore, 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 NG−1 ordinary differential equations for the time at all but the first node.Finally, the three-step equations, explained in equations (16)–(18), were used for time discretization, and equations (30)–(32) were similarly achieved.


## 4. 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. This solution remains bounded when we fix the final time,tF=nΔt, and let n⟶∞. Therefore, a stability analysis applies a restriction to time step.The 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:(41)M⋅P^n+1/3−P^n=−cvΔt3D⋅P^n,(42)M⋅P^n+1/2−P^n=−cvΔt2D⋅P^n+1/3,(43)M.P^n+1−P^n=−cvΔtD⋅P^n+1/2,where P^n=p^1n,p^2n,…,p^NE+1n 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.The 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:(44)p^jn+1/3−p^jn=−cvΔt3h2−p^jn+2p^j+1n−p^j+2n,j=2,…,NE.Equation (44) can be written as follows:(45)p^jn+1/3−p^jn=Cr3δ2p^jn,j=2,…,NE,where Cr=cνΔt/h2 is the diffusion number, and δ2 is the standard notation for the difference operator such that(46)δ2p^jn=p^j−1n−2p^jn+p^j+1n.The round-off error can be defined as follows:(47)εjn=Njn−p^jn,j=2,…,NE,where p^jn and Njn are the numerical and exact nodal values, respectively, which satisfy equation (45). Hence, substituting εjn in equation (45) results in(48)εjn+1/3−εjn=Cr3δ2εjn,j=2,…,NE.In addition, based on the Von Neumann stability analysis, error variation can be displayed by using Fourier series as follows:(49)εjn=∑m=−∞∞Ameimπjhgmπh,Crn,where g≡gmπh,Cr 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(50)εjn=gneikjh,where k=mπ is the wave number. Also, we can write(51)εjn+1/3=gn+1/3eikjhafter substituting εjn and εjn+1/3 in equation (48), we have(52)gn+1/3eikjh−gneikjh=Cr3δ2gneikjh,j=2,…,NE,where the difference operator δ2 can determined through equation (46). After dividing each term in equation (52) by gneikjh, the following equation is obtained:(53)g1/3−1=Cr3e−ikh−2+eikh.Sinceeikh=coskh+isinkh, we have(54)e−ikh−2+eikh=−2+2coskh=−21−coskh.Substituting equation (54) in equation (53) results in(55)g1/3kh,Cr≡g1/3=1−23Cr1−coskh.Equation (55) holds true for j=2,…,NE, therefore standing for the vector of nodal values P^n+1/3. The 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)εjn+1/2−εjn=Cr2δ2εjn+1/3,εjn+1−εjn=Crδ2εjn+1/2.Thus, the amplification factors for the second and third steps are achieved as follows:(57)g1/2kh,Cr≡g1/2=1−Cr1−coskh.g1/3,(58)gkh,Cr≡g=1−2Cr1−coskh.g1/2.Ifξ=kh, by combining equations (55)–(58), the following equation is obtained for the amplification factor g:(59)g=1−2Cr1−cosξ+2Cr21−cosξ2−43Cr31−cosξ3.For a present model problem, we shall say that a method is stable if there exists a constantK, independent of k, such that:(60)gkh,Crn≤K,fornΔt≤tF,∀k.Evidently, the following condition is achieved for stability based on Von Neumann analysis:(61)gkh,Cr≤1+K′Δt,∀k.Accordingly, the necessary and sufficient condition for the error to remain bounded is [32]:(62)gkh,Cr≤1.Finally, the numerical stability condition for the three-step method using equation (59) is achieved as follows:(63)Cr=cνΔth2≤0.6282.Therefore, we have(64)Δt≤0.6282h2cv.Based on [31], the stability condition for the Euler explicit methods is Cr≤0.5. Therefore, 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. Therefore, it has solutions of the following form:(65)ϕz,t=eikz−ωt,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:(66)ω=−ik2.Upon substituting this relation back in to the plane wave, the following solution for the consolidation equation is obtained:(67)ϕz,t=e−cνk2teikz.On the other hand, by makingξ⟶0 in equation (59), the following equation is achieved:(68)g=1−2Crξ22+2Cr2ξ222−43Cr3ξ223.Substituting the corresponding values ofCr and ξ in equation (68) results in(69)g=1−cvk2Δt+12cv2Δt2k4−16cv3Δt3k6.This expression agrees with the first four terms of the power series ofe−cvk2Δt; therefore, the propose three-step method is third-order accurate.


## 5. Numerical Examples

This 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:(70)L∞=max1≤zi≤NGpexactzi,tfinal−p¯zi,tfinal,where p¯ is the approximate solution, obtained from the numerical method, and zii=1NG 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]. Their method can obtain the solution for the one-dimensional consolidation equation under any loading profile (ft), whether the external load is smooth or piecewise smooth. In addition, pore pressure, i.e., pz,t, along with the time derivative of the external loading (f′t) are expressed in Stickle and Pastor’s method as(71)pz,t=∑n=0∞pntsin1+2nπz2L,f′t=∑n=0∞fntsin1+2nπz2L,where the coefficients fnt are the Fourier coefficients of f′t, and pnt can be computed by(72)pnt=4π2n+1e−Nntf0+∫0te−Nnt−τf′τdτ,where Nn=cν1+2nπ/2L2 is the coefficient of consolidation in the vertical direction.The abovementioned equation works very well for external loads. When the loading profile was piecewise smooth, it defined an adequate partition of the interval0,tmax, where tmax is the maximum time that the evolution of pore pressure should be computed. Then, ft was expanded by continuous and Heaviside functions on each interval of the partition. The distribution of ft 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. The 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. Therefore, in this example, we have f′t=0 for the whole consolidation process. By using equation (72), the coefficients pnt are given by(73)pnt=4e−Nntf0π2n+1.
The value off0 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, pz,t.
Table1 indicates a comparison between the three-step solutions with the Euler explicit solutions using linear elements. The 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(74)tfinal=L2Tνcν.
A comparison between the three-step solutions and Euler explicit solutions is presented in Table2 using spectral elements. In the proposed method, 15 elements were considered in total. The pore pressure distribution with second- and fourth-order polynomials was approximated over the even and odd elements, respectively. In this case, 47 total interpolation nodes, i.e., NG=47 existed. Numerical results demonstrated that the three-step method leads to higher stability criteria compared to Euler explicit method when using spectral elements.
Table3 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 (Figures2–3). The three-step 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.
Figures4 and 5 illustrate a comparison between three-step 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.Table 1
A comparison of different time integration methods using linear elements for Example1 (NE=20).
T
ν
Δ
t
Euler explicit method
Three-step method
T
ν
=
0.1
0.0001
4.3241
×
10
−
2
4.3273
×
10
−
2
T
ν
=
0.1
0.004
Unstable
4.2507
×
10
−
2
T
ν
=
0.2
0.0001
3.0363
×
10
−
2
3.0357
×
10
−
2
T
ν
=
0.2
0.004
Unstable
2.8188
×
10
−
2
T
ν
=
0.4
0.0001
1.8127
×
10
−
2
1.8107
×
10
−
2
T
ν
=
0.4
0.004
Unstable
1.6850
×
10
−
2
T
ν
=
0.7
0.0001
8.5773
×
10
−
3
8.5606
×
10
−
3
T
ν
=
0.7
0.004
Unstable
8.4832
×
10
−
3
T
ν
=
0.9
0.0001
5.2176
×
10
−
3
5.2044
×
10
−
3
T
ν
=
0.9
0.004
Unstable
5.2044
×
10
−
3Table 2
A comparison of different time integration methods using spectral elements for Example1 (NE=15).
T
ν
Δ
t
Euler explicit method
Three-step method
T
ν
=
0.1
0.00006
2.8494
×
10
−
4
2.3736
×
10
−
4
T
ν
=
0.1
0.0001
Unstable
2.6388
×
10
−
4
T
ν
=
0.2
0.00006
2.0273
×
10
−
4
1.6847
×
10
−
4
T
ν
=
0.2
0.0001
Unstable
1.7672
×
10
−
4
T
ν
=
0.4
0.00006
1.0152
×
10
−
4
1.1452
×
10
−
4
T
ν
=
0.4
0.0001
Unstable
1.2294
×
10
−
4
T
ν
=
0.7
0.00006
5.4264
×
10
−
5
5.5951
×
10
−
5
T
ν
=
0.7
0.0001
Unstable
5.1930
×
10
−
5
T
ν
=
0.9
0.00006
3.3041
×
10
−
5
3.2254
×
10
−
5
T
ν
=
0.9
0.0001
Unstable
3.7163
×
10
−
5Table 3
A comparison between linear and spectral elements, for Example1 (Δt=0.0001).
T
ν
Three-step method using linear elements
Three-step method using spectral elements
T
ν
=
0.1
4.2507
×
10
−
2
2.6388
×
10
−
4
T
ν
=
0.2
2.8188
×
10
−
2
1.7672
×
10
−
4
T
ν
=
0.4
1.6850
×
10
−
2
1.2294
×
10
−
4
T
ν
=
0.7
8.4832
×
10
−
3
5.1930
×
10
−
5
T
ν
=
0.9
5.2044
×
10
−
3
3.7163
×
10
−
5Figure 2
Solution of Example1 generated by the three-step method with Tν=0.9 and Δt=0.004.Figure 3
Solution of Example1 generated by the Euler explicit method with Tν=0.9 and Δt=0.004.Figure 4
Three-step solutions of Example1 with linear elements for Δt=0.0001. The solid lines correspond to Stickle and Pastor’s analytical solution [10].Figure 5
Three-step solutions of Example1 with spectral elements for Δt=0.0001. The solid lines correspond to Stickle and Pastor’s analytical solution [10].Example 2.
In this example, we consider equation (9) with a cyclic haversine loading. This 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 ft for a haversine loading might be given as(75)ft=qsin2πttc,where tc=0.1L2/cν. This type of loading profile is smooth. Taking time derivative of (75), f′t has the following expression for a haversine loading:(76)f′t=qπtcsin2πttc.
Bearing in mind thatf0=0 in this case of study, substituting (76) in (72) and computing the integral, the following expression is obtained for the coefficients pnt:(77)pnt=4qNntc2n+1Nn2tc2+4π2sin2πttc−2πNntccos2πttc+2πNntce−Nnt.
Substituting (77) in (71) provides the analytical solution for the pore fluid pressure distribution pz,t.
Table4 presents a comparison between the three-step solutions with Euler explicit solutions using linear elements. The first two rows of the table indicates the results obtained for the final time tc/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 Table5, spectral elements are used to compare the three-step and Euler explicit solutions. Similar to the first example, NE=15 and NG=47 were selected.
Table6 presents a comparison between linear and spectral elements with the three-step method for fix Δt=0.0001. The 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.
Figures6 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.
Figures8 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.Table 4
A comparison of different time integration methods using linear elements for Example2 (NE=20).
Final time
Δ
t
Euler explicit method
Three-step method
t
=
t
c
/
2
0.0001
1.9220
×
10
−
3
1.7717
×
10
−
3
t
=
t
c
/
2
0.004
Unstable
5.1297
×
10
−
2
t
=
5
t
c
/
2
0.0001
1.1194
×
10
−
3
1.0416
×
10
−
3
t
=
5
t
c
/
2
0.004
Unstable
1.8287
×
10
−
2
t
=
39
t
c
/
2
0.0001
1.0733
×
10
−
3
9.9699
×
10
−
4
t
=
39
t
c
/
2
0.004
Unstable
1.5856
×
10
−
2
t
=
t
c
0.0001
1.0342
×
10
−
3
9.5968
×
10
−
4
t
=
t
c
0.004
Unstable
1.8101
×
10
−
2
t
=
3
t
c
0.0001
1.0550
×
10
−
3
9.7973
×
10
−
4
t
=
3
t
c
0.004
Unstable
1.6715
×
10
−
2
t
=
20
t
c
0.0001
1.0737
×
10
−
3
9.9719
×
10
−
4
t
=
20
t
c
0.004
Unstable
1.5972
×
10
−
2Table 5
A comparison of different time integration methods using spectral elements for Example2 (NE=15).
Final time
Δ
t
Euler explicit method
Three-step method
t
=
t
c
/
2
0.00001
1.3637
×
10
−
4
1.3578
×
10
−
4
t
=
t
c
/
2
0.0001
Unstable
1.3604
×
10
−
3
t
=
5
t
c
/
2
0.00001
4.7078
×
10
−
5
3.2489
×
10
−
5
t
=
5
t
c
/
2
0.0001
Unstable
4.2143
×
10
−
4
t
=
39
t
c
/
2
0.00001
4.3553
×
10
−
5
2.9253
×
10
−
5
t
=
39
t
c
/
2
0.0001
Unstable
3.6589
×
10
−
4
t
=
t
c
0.00001
4.8212
×
10
−
5
3.4286
×
10
−
5
t
=
t
c
0.0001
Unstable
4.4562
×
10
−
4
t
=
3
t
c
0.00001
4.4852
×
10
−
5
3.0785
×
10
−
5
t
=
3
t
c
0.0001
Unstable
3.8548
×
10
−
4
t
=
20
t
c
0.00001
4.3515
×
10
−
5
2.9196
×
10
−
5
t
=
20
t
c
0.0001
Unstable
3.6534
×
10
−
4Table 6
A comparison between linear and spectral elements, for Example2 (Δt=0.0001).
Final time
Three-step method using linear elements
Three-step method using spectral elements
t
=
t
c
/
2
1.7717
×
10
−
3
1.3604
×
10
−
3
t
=
5
t
c
/
2
1.0416
×
10
−
3
4.2143
×
10
−
4
t
=
39
t
c
/
2
9.9699
×
10
−
4
3.6589
×
10
−
4
t
=
t
c
9.5968
×
10
−
4
4.4562
×
10
−
4
t
=
3
t
c
9.7973
×
10
−
4
3.8548
×
10
−
4
t
=
20
t
c
9.9719
×
10
−
4
3.6534
×
10
−
4Figure 6
Solution of Example2 generated by the three-step method with t=tc/2 and Δt=0.004.Figure 7
Solution of Example2 generated by the Euler explicit method with t=tc/2 and Δt=0.004.Figure 8
Three-step solutions of Example2 with linear elements for Δt=0.0001. The solid lines correspond to Stickle and Pastor’s analytical solution [10].Figure 9
Three-step solutions of Example2 with spectral elements for Δt=0.0001. The solid lines correspond to Stickle and Pastor’s analytical solution [10].


## 6. Conclusions

In the present study, a new numerical method was proposed in order to solve one-dimensional consolidation equation. This 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. The proposed method was compared with the Euler explicit and semianalytical methods, presented by Stickle and Pastor [10]. The 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 high-dimensional consolidation equations.



---

*Source: 2873869-2020-04-28.xml*