Applications of Symmetric and Nonsymmetric MSSOR Preconditioners to Large-Scale Biot's Consolidation Problems with Nonassociated Plasticity

Abstract

Two solution schemes are proposed and compared for large 3D soil consolidation problems with nonassociated plasticity. One solution scheme results in the nonsymmetric linear equations due to the Newton iteration, while the other leads to the symmetric linear systems due to the symmetrized stiffness strategies. To solve the resulting linear systems, the QMR and SQMR solver are employed in conjunction with nonsymmetric and symmetric MSSOR preconditioner, respectively. A simple footing example and a pile-group example are used to assess the performance of the two solution schemes. Numerical results disclose that compared to the Newton iterative scheme, the symmetric stiffness schemes combined with adequate acceleration strategy may lead to a significant reduction in total computer runtime as well as in memory requirement, indicating that the accelerated symmetric stiffness method has considerable potential to be exploited to solve very large problems.

1. Introduction

The wide applications of large-scale finite element analyses for soil consolidation settlements demand efficient iterative solutions. Among various iterative solvers the Krylov subspace iterative methods have enjoyed great popularity so that it was ranked as one of top ten algorithms invented in the 20th century [1]. In large-scale finite element computations, the success of Krylov subspace iterative methods may be partially attributed to the matrix-vector multiplications which can be implemented sparsely in the iterative process, but it is the preconditioning technique that makes Krylov subspace iterative methods become practically useful.

For an indefinite linear system some preconditioning methods (such as incomplete LU or ILU) originally proposed for general problems may encounter slow convergence or even breakdown during the iterative process (e.g., [2]). In the past decade, many researchers have focused on constructing preconditioners for indefinite problems either by revising a standard preconditioner or by designing a block preconditioner according to the block matrix structure. For example, for the indefinite linear equation stemming from Biot’s consolidation a generalized Jacobi (GJ) preconditioner was proposed by Phoon et al. [3] to overcome the unbalanced diagonal scaling of standard Jacobi, while to obviate unstable triangular solves a modified symmetric successive overrelaxation (MSSOR) preconditioner [4] was designed to enhance the pivots embedded within the standard symmetric successive overrelaxation (SSOR) factorization. The successes of GJ and MSSOR undoubtedly are attributed to the fact that both of them were developed based on the block matrix structure. On a separate track, many other researchers have been working on block preconditioners. Perugia and Simoncini [5] proposed symmetric indefinite block diagonal preconditioners for mixed finite element formulations, and, furthermore, Simoncini [6] proposed block triangular preconditioners to couple with symmetric quasiminimal residual (SQMR) [7] for symmetric saddle-point problems. Keller et al. [8] and Bai et al. [9] developed their own block constraint preconditioners for indefinite linear systems. For three types of block preconditioners, an in-depth comparison and eigenvalue study has been carried out by Toh et al. [10], and it was concluded that the block constraint preconditioner and the diagonal GJ preconditioner could be more promising for large-scale Biot’s linear system. Recently, more attentions have been paid to the constraint block structure. For instance, Bergamaschi et al. [11], Ferronato et al. [12], and Haga et al. [13] defined their own block constraint preconditioners, respectively, and validated the good performance of their block constraint preconditioners. Furthermore, Janna et al. [14] developed a parallel block preconditioner for symmetric positive definite (SPD) matrices by coupling the generalized factored sparse approximate inverse (FSAI) with ILU factorization, and it was concluded that in many realistic cases the FSAI-ILU preconditioner is more effective than both FSAI and ILU preconditioners. To resolve the soil-structure problems in which material stiffness contrast is significant, the partitioned block diagonal preconditioners were proposed by Chaudhary [15], and numerical results indicate that the proposed preconditioners are not sensitive to the material stiffness contrast ratio. Based on the above review, it seems that constructing a preconditioner based on the block structure of the coefficient matrix has become a popular strategy for large-scale Biot’s indefinite linear systems.

Soil consolidation accompanied with soil dilatancy may be frequently encountered in practice, and commonly the soil dilatancy may be modeled by nonassociated plasticity. The main contribution of the article is that for coupled consolidation problems involving nonassociated plasticity, two solution strategies are proposed and evaluated, respectively, by employing a strategy of combining a nonlinear iterative scheme with a corresponding linear iterative method (e.g., [16, 17]). The paper is organized as follows. In Section 2, the nonlinear iterative methods are formulated for elastoplastic soil consolidation problems, and the coupled 2×2 nonsymmetric indefinite Biot’s finite element linear system of equation due to nonassociated plastic flow is derived. In Section 3, the recently proposed block preconditioners are reviewed and commented. In Section 4, both symmetric version and nonsymmetric version of MSSOR are implemented within the proposed two schemes, respectively, to solve the soil consolidation examples with nonassociated plasticity, and their performances are investigated and compared. Finally, some useful observations and conclusions are summarized in Section 5.

2. Coupled Linear System of Equation Arising from Elastoplastic Biot’s Consolidation Problems

Soil consolidation is a soil settlement process coupled with dissipation of excess pore water pressure, and this process may be physically modeled by the widely accepted Biot’s consolidation theory [18]. Recall that for a fully saturated porous media, the volumetric fluid content variation within the soil skeleton is solely related to the deformation of solid skeleton, that is, ∇⋅𝝈(𝐱,𝑡)−𝛾𝐛=0,(𝐱,𝑡)∈Ω×(0,𝑇], (2.1) in which the Terzaghi’s effective stress principle 𝝈=𝝈^′+𝑝 is applied; 𝐛=[0,0,1]^𝑇 is the unit body force vector; γ is the total unit weight; Ω is the solution domain; 𝑇 is the total analyzing time. For the fully saturated porous media, the pore fluid should comply with the fluid continuity equation or the mass-conservation equation, that is,∇⋅𝐮̇(𝐱,𝑡)+∇⋅𝐯_𝑓(𝐱,𝑡)=0, (2.2) where ∇⋅𝐮̇=𝜀_𝑣̇ with 𝜀_𝑣 as the volumetric strain, the dot symbol over any symbols means time differentiation, and u is the displacement vector; the fluid velocity 𝐯_𝑓 is described by Darcy’s law𝐯_𝑓[𝐤]=−𝛾_𝑓∇𝑝−𝛾_𝑓𝐛. (2.3) Here, 𝛾_𝑓 is the unit weight of pore fluid; [k] = [k_s] is the hydraulic conductivity tensor (it is a diagonal matrix [𝐤_𝑠]=diag(𝑘_𝑠,𝑥,𝑘_𝑠,𝑦,𝑘_𝑠,𝑧) when orthogonal hydraulic conductivity properties are assumed) for saturated flow. The solution domain Ω has the complementary boundary condition𝑠 may be given as 𝜕Ω=𝜕Ω_𝐷∪𝜕Ω_𝑁 and 𝜕Ω_𝐷∩𝜕Ω_𝑁=∅ with 𝜕Ω_𝐷 for the Dirichlet (or prescribed displacement and pore pressure) boundary and 𝜕Ω_𝑁 for the Neumann (or prescribed force and flu𝑥) boundary. After discretizing (2.1) and (2.2) in space and time domain, respectively, the incremental form of Biot’s finite element equation is derived as (e.g., [19]):𝐋𝐊𝐋_𝑇−𝜃Δ𝑡𝐆Δ𝐮Δ𝐩^𝑒𝑥=Δ𝐟Δ𝑡𝐆𝐩_𝑛^𝑒𝑥 (2.4) in which the subblocks within the coupled matrix are defined, respectively, as𝐊=_𝑒𝑉𝑒𝐁^𝑇_𝑢𝐃_𝑒𝑝𝐁_𝑢𝑑𝑉,𝐋=_𝑒𝑉𝑒𝐁^𝑇_𝑢𝟏𝐍_𝑝𝑑𝑉,𝐆=_𝑒𝑉𝑒𝐁^𝑇_𝑝[𝐤]𝛾_𝑓𝐁_𝑝𝑑𝑉, (2.5) where 𝐃_𝑒𝑝 is the constitutive matrix; 𝐁_𝑢, 𝐁_𝑝 is the gradient of shape function 𝐍_𝑢 and 𝐍_𝑝 used for interpolating the displacement u and the excess pore pressure 𝐩^𝑒𝑥, respectively; 𝟏=[111000]^𝑇; Δ𝑡 is the time increment; θ is the time integration parameter ranging from 0 to 1; θ =1/2 corresponds to the second-order accuracy Crank-Nicolson method, while θ =1 leads to the fully implicit method possessing the first-order accuracy.

Equation (2.4) is the discretized finite element equation for each time increment, and to be convenient the following weak form of the residual equation is used to derive the nonlinear iterative process,𝐑𝐮;𝐩^𝑒𝑥=𝐑_𝑢𝐑_𝑝=𝐅^ext_𝑢𝐅^ext_𝑝−𝐅_𝑢^int𝐅_𝑝^int=𝟎𝟎, (2.6) where R is the residual or out-of-balance force vector derived from the total potential energy Φ(𝐮;𝐩^𝑒𝑥); 𝐅^ext is the applied external force at current time increment, and its two parts corresponding to u and 𝐩^𝑒𝑥, respectively, have been presented in the RHS of (2.4). Noting that the deformation of soil skeleton is solely caused by the effective stress, the internal forces 𝐅_𝑢^int and 𝐅_𝑝^int may be expressed, respectively, as below,𝐅_𝑢^int=_𝑉𝐁^𝑇_𝑢𝝈^′𝑑𝑉+_𝑉𝐁^𝑇_𝑢𝐩^𝑒𝑥𝐅𝑑𝑉,_𝑝int=_𝑉𝐍^𝑇_𝑝𝜺𝑑𝑉+𝜃Δ𝑡_𝑉𝐁^𝑇_𝑝𝐯_𝑓𝑑𝑉. (2.7) Applying the first-order Taylor expansion to (2.6) leads to the following Newton-Raphson (NR) iteration:𝐑𝐮_𝑘−1;𝐩^𝑒𝑥_𝑘−1+𝐀_𝑘−1𝛿𝐮_𝑘𝛿𝐩_𝑘^𝑒𝑥𝐮=0,_𝑘𝐩_𝑘^𝑒𝑥=𝐮_𝑘−1𝐩^𝑒𝑥_𝑘−1+𝛿𝐮_𝑘𝛿𝐩_𝑘^𝑒𝑥, (2.8)(2.9) where𝐀_𝑘−1=𝜕𝐑_𝑢𝜕𝐮_𝑘−1𝜕𝐑_𝑢𝜕𝐩^𝑒𝑥𝜕𝐑_𝑝𝜕𝐮𝜕𝐑_𝑝𝜕𝐩^𝑒𝑥,𝜕𝐑_𝑢𝜕𝐩^𝑒𝑥=𝜕𝐑_𝑝𝜕𝐮^𝑇=−𝐋,𝜕𝐑_𝑝𝜕𝐩^𝑒𝑥=𝜃Δ𝑡𝐆 (2.10) in which it is apparent that L is independent of time increment and iterative process due to small-strain finite element computation, and G is independent of time increment and iterative process due to the saturated soil assumption. In 𝐀_𝑘−1, the gradient of 𝐑_𝑢 about u is the tangential solid stiffness matrix,𝐊_𝑇=−𝜕𝐑_𝑢=𝜕𝐮𝜕𝐅_𝑢^int(𝐮),𝜕𝐮 (2.11) and 𝐊_𝑇 is assembled by the element stiffness 𝐤_𝑒𝑝, that is,𝐊_𝑇=_elements𝜕∫_𝑉𝑒𝐁^𝑇𝐃_𝑒𝑝𝐁𝑑𝑉𝐮_𝑒𝜕𝐮=_elements𝐤_𝑒𝑝with𝐤_𝑒𝑝=_𝑉𝑒𝐁^𝑇𝐃_𝑒𝑝𝐁𝑑𝑉. (2.12) Evidently, the symmetry of 𝐊_𝑇 or the 2×2 coupled matrix 𝐀_𝑘−1 in (2.10) is governed by 𝐃_𝑒𝑝, which may be expressed as𝐃_𝑒𝑝=𝐃_𝑒−𝐃_𝑝=𝐃_𝑒−𝐃_𝑒𝐛𝐚^𝑇𝐃_𝑒𝐚^𝑇𝐃_𝑒𝐛−𝐜^𝑇𝐡, (2.13) where 𝐚=𝜕_𝜎𝐹 represents the yield surface normal vector, 𝐛=𝜕_𝜎𝐺 denotes the plastic flow direction vector, while 𝐜=𝜕_𝐪𝐹 is the gradient of yield function (i.e., 𝐹) about the internal variables of the soil model. In conventional elastoplastic soil models, the physically associated plastic flow leads to symmetric 𝐃_𝑒𝑝 due to 𝐛=𝐚, while the nonassociated plastic flow results in nonsymmetric 𝐃_𝑒𝑝 due to 𝐛≠𝐚. For convenience, (2.8) may be further expressed as𝐊_𝑘−1𝐋𝐋^𝑇−𝜃Δ𝑡𝐆𝛿𝐮𝛿𝐩^𝑒𝑥_𝑘=𝐑_𝑘−1=𝐅^ext−𝐅^int_𝑘−1, (2.14) where the superscript 𝑘 signifies the nonlinear iteration count within each time increment. To distinguish the nonlinear iterative process from the Krylov subspace iterative process, the iteration count for the first process as described by (2.14) is called nonlinear iteration count, while the second Krylov subspace process is a linear iterative process, and the iteration count for this process is called linear iteration count.

In the nonlinear iterative process, the search direction 𝐝_𝑘 may be computed inexactly (i.e., 𝐝_𝑘≈[𝛿𝐮;𝛿𝐩^𝑒𝑥]^𝑇_𝑘), leading to the so-called inexact Newton method whose convergence is governed by (e.g., [16, 17]),𝐑_𝑘−1+𝐀_𝑘−1𝐝_𝑘≤𝜂_𝑘−1𝐑_𝑘−1, (2.15) where 𝜂_𝑘−1∈ [0, 1) is the forcing term. With the computed search direction 𝐝_𝑘, the displacement is updated by𝐮_𝑘𝐩_𝑘^𝑒𝑥=𝐮_𝑘−1𝐩^𝑒𝑥_𝑘−1+𝜒_𝑘𝛿𝐮_𝑘𝛿𝐩_𝑘^𝑒𝑥 (2.16) in which 𝜒_𝑘 is the step-length parameter, which may be determined by some optimization strategies. The inexact computation of 𝐝_𝑘 may arise from the approximate solve of 𝐝_𝑘 or the approximation to 𝐀_𝑘−1.

For ease of presentation, the linear system of equation at each nonlinear iteration as shown by (2.14) is expressed concisely as𝐵𝐾𝐵_𝑇𝑥𝑦=𝑓−𝐶𝑔𝐵with𝐴=𝐾𝐵_𝑇−𝐶_𝑁×𝑁. (2.17)

3. Block-Structured Preconditioners

Here, the block-structured preconditioners are defined as those developed according to the coefficient matrix block structure, and hence they include the block preconditioners and those modified preconditioners inspired by the matrix block structure.

3.1. Block Preconditioners for Biot’s Linear System of Equation

Block preconditioners can be classified into three types as mentioned previously. According to Toh et al. [10], it could be convenient to derive the preconditioners based on the matrix inverse. Note that the inverse of the 2×2 matrix in (2.17) may be expressed as𝐾𝐵₂𝐵^𝑇₁−𝐶⁻¹=𝐾⁻¹−𝐾⁻¹𝐵₂𝑆⁻¹𝐵^𝑇₁𝐾⁻¹𝐾⁻¹𝐵₂𝑆⁻¹𝑆⁻¹𝐵^𝑇₁𝐾⁻¹−𝑆⁻¹ (3.1) in which 𝑆=𝐶+𝐵^𝑇₁𝐾⁻¹𝐵₂ (with 𝐵₁=𝐵₂) is the Schur complement matrix. Given a vector [𝑢; 𝑣] and a block preconditioner 𝑀≈𝐴 with the approximation 𝐾 and 𝑆 to 𝐾 and 𝑆, respectively. Hence, the preconditioning step may be written as the following.

Preconditioning Step 𝑀⁻¹[𝑢;𝑣]:
solve𝑝=𝐾⁻¹𝑢,solve𝑞=𝑆⁻¹(𝐵^𝑇₁𝑝−𝑣),compute𝑀⁻¹[𝑢,𝑣]^𝑇=[𝐾⁻¹(𝑢−𝐵₂𝑞);𝑞].

The above preconditioning step has already been introduced in [4, 10] for the block constraint preconditioners; however, it is highlighted here as a unified computational framework for all block preconditioners. It can be observed that within the above preconditioning step, canceling the terms associated with off-diagonal subblocks 𝐵₁ and 𝐵₂ corresponds to the block diagonal preconditioning, while canceling either 𝐵₁ or 𝐵₂ corresponds to the block triangular types. Depending on how 𝐾 and 𝑆 are approximated, various block preconditioners may be developed. Regardless of the type of preconditioners adopted, solving the linear systems associated with 𝐾 and 𝑆 can not be obviated. Consequently, how 𝐾 and 𝑆 are approximated and how these linear systems are solved efficiently are the key features distinguishing recently proposed block preconditioners. For example, Phoon et al. [3] proposed the diagonal approximations to both 𝐾 and 𝑆 so that the inverses of 𝐾 and 𝑆 can be directly attained. While in the block constraint preconditioner proposed by Toh et al. [10], because it is observed that the size of 𝐾 is significantly larger than that of 𝑆, a diagonal approximation is employed for 𝐾, and the incomplete or sparse Cholesky factorization is recommended for 𝑆. In another version of block constraint preconditioner proposed by Bergamaschi et al. [11, 12], the incomplete Cholesky factorization and factorized approximate inverse are adopted for 𝐾 and 𝑆, respectively. Different from the previous studies, Haga et al. [13] recommended solving the linear systems 𝐾 and 𝑆 using the algebraic multigrid (AMG) method. To summarize the key differences, these block preconditioners with the proposed 𝐾 and 𝑆 as well as the corresponding solution schemes are tabulated in Table 1. In the table, all block preconditioners are categorized into three types including the block diagonal preconditioners, block triangular preconditioners, and block constraint preconditioners which are denoted as “diagonal,” “constraint,” and “triangular,” respectively.

Table 1 

Block preconditioners recently proposed for Biot’s linear system.

3.2. Modified Preconditioners Based on Matrix Block Structure

There are some advantages to developing preconditioners based on the matrix block structure, particularly in the case of block indefinite linear systems. It is natural to hope that modifying some standard or general preconditioning techniques may lead to improved versions. For instance, the generalized Jacobi preconditioner [3] may be viewed as an improved version of standard Jacobi preconditioner by referring to the block diagonal. The MSSOR preconditioner [4] may also be regarded as a good example that demonstrates how to develop a preconditioner based on a standard preconditioner to suit indefinite problems. The general nonsymmetric form of MSSOR may be expressed as below𝑀_MSSOR=𝐿_𝐴+𝐷𝜔𝐷𝜔⁻¹𝑈_𝐴+𝐷𝜔=𝐿_𝐴+𝐷𝐷⁻¹𝑈_𝐴+𝐷, (3.2) where 𝜔∈[1,2]is the relaxation parameter;𝐿_𝐴, 𝑈_𝐴 is the strictly lower and upper part of matrix 𝐴, respectively, that is, 𝐴=𝐿_𝐴+𝑈_𝐴+𝐷_𝐴. For symmetric linear equation, 𝐿_𝐴=𝑈^𝑇_𝐴 leads to the symmetric MSSOR. 𝐷 is the modified diagonal and it is recommended to take GJ, that is,𝐷=𝑀_GJ=(𝐾)0𝑆diag0𝛼diag (3.3) in which α is a scaling factor which is recommended to be −4 according to the eigenvalue study. It is clear that MSSOR preconditioner can be implemented efficiently when combined with the Eisenstat trick [20]. In the following part, both symmetric version and nonsymmetric version of MSSOR preconditioner will be employed for large-scale nonlinear soil consolidation involving nonassociated plasticity.

4. Numerical Experiments

To simulate soil dilatancy during the soil consolidation process, an elastoplastic soil model with nonassociated plasticity should be used [21]. The nonassociated plasticity theory is a generalization of the classical elastoplastic theory by introducing a new plastic potential function [22], and in the past decades the nonassociated plasticity theory is widely applied in practical finite element analyses.

It is also well known that when the nonassociated plastic soil models are employed in finite element analysis, the resulting linear systems of equations are nonsymmetric. Solving the linear equations separately from the nonlinear iterative procedure may not be wise, because an appropriate combination between the outer nonlinear iterative scheme and an inner linear iterative solution may lead to a significant reduction in computer runtime without sacrificing accuracy [17]. In this work, two solution schemes are proposed and compared for the target problem.

Scheme 1. Applying the Newton-Krylov iterative method. As shown by (2.8), the resultant linear system is nonsymmetric at each Newton iteration, a nonsymmetric Krylov subspace iterative method, such as quasiminimal residual (QMR) [23], is adopted, and hence the nonsymmetric MSSOR preconditioner is used to accelerate its convergence.

Scheme 2. Compared to the nonsymmetric linear systems, solving the symmetric linear systems could lead to a saving in required memory and computer runtime. Therefore, accelerated nonlinear solution schemes with symmetrized stiffness are attempted here, that is, at each nonlinear iteration one attempts to solve 𝐑𝐮_𝑘−1;𝐩^𝑒𝑥_𝑘−1+𝐀_𝑘−1𝛿𝐮_𝑘𝛿𝐩_𝑘^𝑒𝑥=0, (4.1) and the solution vector may be updated according to (2.16). The difficulty associated with the scheme is how to construct the symmetric linear systems. In the accelerated symmetric stiffness method proposed by Chen and Cheng [24], the idea is to construct the symmetric constitutive matrix as sym𝐃_𝑒𝑝=𝐃_𝑒−sym𝐃_𝑝, (4.2) where sym(·) is a symmetrizing symbol. When sym(𝐃_𝑒𝑝)=𝐃_𝑒, it corresponds to the initial stiffness (IS) method proposed by Zienkiewicz et al. [25]. When it is defined that sym𝐃_𝑒𝑝=𝐃^𝐺_𝑒𝑝=𝐃_𝑒−𝐃_𝑒𝐛𝐛^𝑇𝐃_𝑒𝐛^𝑇𝐃_𝑒𝐛−𝐜^𝑇_𝐺𝐡_𝐺 (4.3) which corresponds to the so-called accelerated 𝐾_𝐺 approach. In (2.16), the step-length parameter 𝜒_𝑘 can be determined by some optimization strategies such as the Chen’s method [26] and Thomas’ method [27]. Chen’s method is derived by minimizing the out-of-balance load at next iteration in the least-squares sense, while the Thomas method is proposed by minimizing the “symmetric” displacement at the next iteration in the least-squares sense. In this study, the solution vector consists of coupled displacement and excess pore water pressure degrees of freedom (DOFs). The Chen method may be more straightforward to apply in this case, which will be demonstrated by the following numerical experiments. As a result, by using the present scheme, the symmetric MSSOR preconditioner may be adopted to accelerate the convergence of the SQMR solver.

4.1. Convergence Criteria

In the examples to be studied, the nonlinear iteration is terminated if the relative residual force criterion is satisfied,𝐑_𝑘2𝐅^ext₂≤Tol_NLor𝑘≥Maxit_NL,(𝑘=0,1,…) (4.4) in which Tol_NL and Maxit_NL are the prescribed tolerance and the maximum nonlinear iterative number, respectively, and in this study Tol_NL and Maxit_NL are set as 0.01 and 50. For the employed Krylov subspace iterative method, the relative residual convergence criterion is adopted,𝐑_𝑘−1−𝐀_𝑘−1𝛿𝐮;𝛿𝐩^𝑒𝑥𝑇_𝑘,𝑗2𝐑_𝑘−12≤Tol_Lin,or𝑗≥Maxit_Lin,𝑗=0,1,…, (4.5) in which zero initial guess is assumed; 𝑗 is the linear iterative index; Tol_Lin and Maxit_Lin are the prescribed tolerance and the maximum iterative number, respectively, and in this study Maxit_Lin is set as 50000. While to achieve a better performance for the adopted nonlinear scheme, a combined tolerance (in which Tol_Lin = 1.E−5 is used if the residual load is the external applied force, or Tol_Lin = 1.E−3 is used if the residual load is the out-of-balance force) is employed, as recommended by Chen and Cheng [24]. Based on our numerical experiences, this recommendation is reasonable because the deformation of a soil body induced by an external load is usually large and should be solved more accurately, while corrections of the deformation to resolve out-of-balance loads are relatively smaller in magnitude and hence may be solved with a less strict tolerance.

In addition, it should be mentioned that the uniform substepping method with 𝑛substep=500 (i.e., the number of substeps) is adopted for stress-strain integration. For more advanced automatic substepping algorithms, see Abbo [28]. In the present study, an ordinary personal computer platform equipped with a 2.4 GHz Intel Core(TM)2 Duo CPU and 3 GB physical memory is used.

4.2. Homogeneous Flexible Footing

To investigate and compare the two schemes proposed, a simple flexible footing problem with homogeneous soil material is simulated. For the homogenous soil property, the hydraulic conductivity is assumed to be 𝑘_𝑠,𝑥 = 𝑘_𝑠,𝑦 = 𝑘_𝑠,𝑧 = 10⁻⁹ m/s, the effective Young’s modulus 𝐸^′ as 10 MPa, and the Poisson’s ratio as 𝜈^′=0.3. The nonassociated Mohr-Coulomb soil model is used to simulate soil plasticity with the cohesion 𝑐^′=10 kPa, the internal friction angle 𝜙^′ = 30° and the dilatancy angle φ = 5°. The 𝐾₀ approach is used to generate the initial field stress with 𝐾₀=1−sin𝜙^′ and soil unit weight is γ = 18 kN/m³. Due to the geometric symmetry, only a quadrant of the footing is modeled with the solution domain discretized by 8000 (20×20×20) 20-node hexahedral consolidation elements as shown in Figure 1, and the resultant total number of DOFs is 107180. For the boundary conditions, standard displacement fixities are assumed, and only the ground surface is drained. The uniform pressure load is applied on a 1×1 m² area. It is increased incrementally using a total of 5 load steps. Each load increment is 0.1 MPa per second.

Figure 1 

The 20×20×20 finite element mesh of the flexible footing.

Table 2 provides the numerical performance of two solution schemes proposed above, that is, the Newton-Krylov (i.e., Newton-QMR) solution scheme and the Chen accelerated symmetric stiffness scheme. In the two schemes, the nonsymmetric MSSOR and symmetric MSSOR are used in conjunction with QMR and SQMR solver, respectively. Because the combined tolerance is used, the required linear iterative count for the first nonlinear iteration is higher than the following nonlinear iterations. When comparing the two schemes, it is found that the average iterative count required by MSSOR-preconditioned QMR is about two times of that required by MSSOR preconditioned SQMR solver, explaining the final results. In addition, the fact that two matrix-vector products are required in each iteration of QMR contributes to the longer computer runtime. Note that the average iterative count required by SQMR in Chen accelerated IS method is slightly smaller than that required by SQMR in Chen accelerated 𝐾_𝐺 method, indicating that the initial stiffness (i.e., the elastic stiffness) could be better conditioned. When observing the required nonlinear iterations by each scheme, it is interesting to note that three solution strategies (i.e., the Newton scheme, the Chen accelerated IS method, and Chen accelerated 𝐾_𝐺 method) exhibit similar nonlinear iterative behavior, although the Chen accelerated 𝐾_𝐺 method leads to slightly less nonlinear iterations (i.e., 18) than the other two counterparts. The exhibited similarity in nonlinear convergence indicates that the nonlinearity in the present problem is not strong. Obviously, for the homogeneous flexible footing problem the Chen accelerated 𝐾_𝐺 method coupled with MSSOR preconditioned SQMR solver may lead to a 66% reduction in total computer runtime than the Newton-QMR scheme preconditioned by nonsymmetric MSSOR. Even for the Chen accelerated IS method, a 57% reduction of total computer runtime can be achieved. It is noteworthy that the Thomas’ acceleration strategy appears to be more effective in a past study [24]. However, for the coupled consolidation problem discussed in the present study, the Thomas accelerated IS method requires 1, 4, 7, 8, 8 nonlinear iterations, respectively, for the five load steps, while the Thomas accelerated 𝐾_𝐺 method takes 1, 5, 8, 6 nonlinear iterations, respectively, for the first four load steps, but it fails to converge at the last load step.

Table 2 

Performance of different solution schemes for the flexible footing example.

4.3. Heterogeneous Soil-Structure Interaction Example

A pile-group example is used to demonstrate the interaction between soil and structure. Because of the significant contrast in material stiffness between soils and concrete, the performance of the MSSOR preconditioner may deteriorate with the increasing contrast ratio, as noticed and investigated by Chaudhary [15]. As shown in Figure 2, the pile-group finite element mesh has total 4944 elements, in which 430 are concrete elements. The material properties for soils and concrete are tabulated in Table 3. From the table, it is seen that the concrete material is modeled by linear elastic consolidation elements with extremely small permeability, but the soils are still simulated by the Mohr-Coulomb model.

Table 3 

Material properties for the pile-group foundation.

Figure 2 

The finite element mesh of heterogeneous pile-group example.

Table 4 provides numerical results of these solution schemes for the pile-group example. To be efficient for the accelerated symmetric stiffness methods, the combined tolerance introduced in Section 4.1 is still adopted. In the pile-group example, three uniform load steps are simulated and the pressure increment of 50 kPa is applied on the pile cap for a 2-day time increment. Similar to the flexible footing example, the iteration number spent by nonsymmetric MSSOR-preconditioned QMR is about two times of that of MSSOR-preconditioned SQMR solver. From the nonlinear iteration counts, it is seen that the Chen accelerated 𝐾_𝐺 method achieves a similar convergence rate as that of the Newton method, but at the third load step QMR solver dose not converge at one nonlinear iteration, which denoted by the bracketed number. Hence, at that load step the average iterative number for each linear system is remarkably increased because of Maxit_Lin = 50000. Compared to the Newton method, the Chen accelerated symmetric stiffness methods show better convergence behaviors, while the Chen accelerated IS method may be markedly slower than Newton method, indicating that the soil-structure interaction involving significant contrast in material stiffness may produce a stronger nonlinearity than the simple homogeneous footing problem. Upon close examination of the computer runtime, it is noticed that even though it is slower in convergence, the Chen accelerated IS method may lead to a reduction of 63% in computer runtime compared with the Newton method. When using 𝐾_𝐺 symmetric stiffness, the resultant convergence rate is similar to the Newton method, the saving in computer runtime is more impressive and a reduction of 82% may be attained. Furthermore, the Thomas acceleration scheme is examined for the initial stiffness and the 𝐾_𝐺 stiffness matrix, respectively. Both nonlinear solution strategies associated with the Thomas acceleration scheme fail, indicating that the Thomas acceleration scheme may not be suitable for problems involving coupled pore water pressure and the displacement degrees of freedom.

Table 4 

Performance of different solution schemes for the pile-group example.

5. Conclusions

Soil consolidation accompanied with soil dilatancy may be frequently encountered in practice, and usually it can be modeled by nonassociated soil models. Due to the nonassociated soil model, the resultant finite element linear equation is nonsymmetric. To solve the nonsymmetric coupled Biot’s linear systems of equations, two schemes in conjunction with MSSOR preconditioner are proposed and examined. Some useful observations are summarized as follows.(1)Two schemes are proposed for the Biot’s consolidation problem involving nonassociated plasticity. Depending on the discretized linear equations, both nonsymmetric and symmetric MSSOR preconditioners are adopted for such problems,(2)In the accelerated symmetric stiffness methods for coupled consolidation problems, the Thomas’ acceleration strategy does not exhibit better convergence behaviors as observed in the previous studies. On the other hand, the Chen’s acceleration strategy is more effective, and hence it is the recommended approach for such coupled consolidation problems.(3)Compared to the Newton solution scheme which adopts the QMR solver preconditioned by nonsymmetric MSSOR preconditioner, the Chen accelerated symmetric stiffness approaches, which use the SQMR solver preconditioned by the symmetric MSSOR, may lead to a significant reduction in computer runtime. Based on the above numerical experiments, it may be concluded that the Chen accelerated symmetric stiffness methods have considerable potential to be exploited for solution of large-scale Biot’s consolidation problems.

Acknowledgment

The research is supported in part by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, no. (2010) 1561, the Fundamental Research Funds for the Central Universities, no. 2011JBM070, and the Research Fund Program of State Key Laboratory of Hydroscience and Engineering, no. SKL-2010-TC-2.