Domain decomposition methods for 3D crack propagation problems using XFEM

Abstract

The extended finite element method (XFEM) has been successfully implemented in solving crack propagation problems by enriching the polynomial basis functions of standard finite elements with specialized non-smooth functions. The resulting approximation space can be used to solve problems with moving discontinuities, such as cracks, without the need of remeshing in the vicinity of the crack. The enrichment of the displacement field in XFEM inflicts a substantial increase in the ellipticity of the discretized problem. As a consequence, the resulting algebraic systems become strongly ill-conditioned, rendering the convergence of iterative solvers very slow. On the other hand, direct solvers may become inefficient in 3D problems, due to the increased bandwidth of the system matrix. In this paper, two of the most efficient domain decomposition solvers, namely the FETI-DP and P-FETI-DP, are proposed for solving the linear systems resulting from XFEM crack propagation analysis in large-scale 3D problems. Both solvers are amenable to parallelization and can be implemented in modern parallel computing environments, with multicore processors and distributed memory systems, following appropriate modifications, to achieve a drastic reduction of memory and computing time in computationally intensive problems.

Introduction

The extended finite element method (XFEM) is one of the most popular methods to simulate fracture phenomena and was originally developed by Belytschko and Black [1] and Moës et al. [2] for analyzing brittle crack propagation in 2D problems. In XFEM, the polynomial approximation space is enriched with problem specific functions, in order to accurately model the discontinuous displacement field and the singular stress field around a crack. These enriched basis functions correspond to additional degrees of freedom (dof), which are introduced around a crack and are called enriched dof to differentiate from the standard dof that express nodal displacements. Thus, simple meshes that do not conform to the crack geometry can be used without the need of remeshing, whenever the crack grows, and mapping the displacement field between the old and new meshes. Ever since, XFEM has been adopted for a wide range of applications, such as fluid–structure interaction [3], contact problems [4], topology optimization [5], probabilistic shape optimization [6] and heat transfer in composite materials [7].

Despite its general success, XFEM has certain shortcomings, especially with regard to the formulation and structure of the resulting linear systems of equations and the ill-conditioning of the corresponding stiffness matrices. The enriched basis functions inflict a substantial increase on the ellipticity of the discretized systems of equations, which become ill-conditioned and inhibit the convergence of iterative-based solvers. On the other hand, although the performance of direct solvers is not generally affected by the ill-conditioning of the system, they may become inefficient due to the increased topological dispersity of stiffness matrix entries, especially in 3D crack propagation problems. Consequently, various specialized solvers have been developed for the solution of XFEM crack propagation problems.

A preconditioning scheme to improve the convergence rate of standard iterative solvers was developed by Bechet et al. [8], based on the Cholesky decomposition of certain node-level submatrices of the global stiffness matrix. Similarly, in [9], a geometric preconditioner constructed from the nodal basis functions is proposed to eliminate the ill-conditioning caused by the Heaviside enrichment in problems with material interfaces. An exact reanalysis direct solver which updates the factorized matrix at each crack propagation step, instead of rebuilding and refactorizing at the global level, is implemented by Pais et al. [10]. Another reanalysis-type algorithm is featured in [11], where the transfer operations of a geometric multigrid solver are established at the beginning of the analysis and then reused at each crack propagation step. The method presented in [12] was also based on geometric multigrid, but added finer mesh patches around small cracks. Moreover, the algebraic multigrid solver of [13] modifies the sparsity pattern of the prolongator operator to prevent interpolation across cracks. In contrast, Gerstenberger and Tuminaro [14] introduced a simple modification of algebraic multigrid, in order to use black-box AMG software.

Domain decomposition methods (DDM) are widely considered as the most computationally efficient solvers for large-scale problems, particularly in parallel computing architectures. The first DDM for XFEM crack propagation was proposed in [15] with the aim of reusing existing FEM software. It is based on the FETI method [16], using only two subdomains, a large uncracked subdomain, assigned to a general purpose FEM software and a smaller subdomain, located around a crack which is modeled with XFEM. A more performance-oriented approach was proposed by Menk and Bordas [17], where the domain was separated into one subdomain containing all standard dof and multiple subdomains containing the enriched ones. The Cholesky factorization was applied to stiffness matrices of enriched subdomains and the QR factorization to matrices connecting them with the large monolithic subdomain with standard dof. The resulting DD matrix was used as a preconditioner, which was effective at reducing the number of iterations in 2D problems, but did not scale well, since the convergence rate decreased when increasing the number of subdomains. Furthermore, Waisman and Berger-Vergiat [18] implemented a multiplicative Schwarz domain decomposition preconditioner to accelerate the convergence of a generalized minimum residual solver in 2D problems. The domain was partitioned into one uncracked subdomain, which was treated with an algebraic multigrid approach, and many smaller subdomains defined around cracks, which were concurrently solved with direct methods. A similar approach for 2D crack propagation problems was presented by Chen and Cai [19], where the preconditioner was based on an additive Schwarz DD solver. This method used a LU factorization for subdomains with enriched dof and an incomplete LU for subdomains with standard dof. Even though the cracked subdomains were still dependent on the locations of cracks, multiple uncracked subdomains could be used. Neither [18] nor [19] scaled well with the number of subdomains, since by increasing the subdomains, an increase on the required iterations was observed. Recently, Agathos et al. [20] developed a deflated conjugate gradient solver, where the deflation space is composed of the rigid body modes of subdomains, as well as additional deflation vectors due to the XFEM enrichments near a crack. The deflation is also combined with a block-Jacobi preconditioner, which uses the subdomains as blocks.

In the present paper we propose the use of two well-established domain decomposition solvers, namely the FETI-DP (Farhat et al. [21]) and P-FETI-DP (Fragakis and Papadrakakis [22]), in order to solve the linear systems resulting from 3D XFEM crack propagation analysis. Instead of decomposing the domain into cracked and uncracked subdomains, the domain is partitioned into an arbitrary number of load-balanced subdomains, independently from the location of cracks, while treating both standard and enriched dof consistently. The customizable coarse problem of FETI-DP and P-FETI-DP allows us to introduce XFEM-specific modifications that eliminate the ill-conditioning caused by the enrichment functions, as well as any singularities that the cracks may induce as they propagate throughout subdomains. The reduced bandwidth of subdomain-level stiffness matrices, afforded by the flexible partitioning, results in a drastic decrease of computing time and memory requirements for solving large-scale 3D crack propagation problems. In contrast to the domain decomposition methods for XFEM crack propagation in [15], [17], [18], [19], the solvers proposed in this work exhibit excellent numerical scalability, when increasing the number of subdomains. This scalability is essential for efficiently implementing them in modern parallel computing environments, where the available memory and processing power can be arbitrarily increased by including additional multicore CPUs and GPUs in distributed memory systems.