Exact XFEM Reanalysis
Introduction
While the extended finite element method has reduced many of the complications involved with modeling crack growth in the finite element framework, there are still limitations on the number of crack growth iterations one may consider. In particular when modeling crack growth there are two main approaches. The first assumes a constant increment of crack growth at each cycle1. The other option is to assume a constant number of cycles and apply a fatigue crack growth law such as Paris Law to predict the crack growth increment for the fixed number of cycles2. Depending upon the mesh and incremental crack growth or fatigue cycle the path of crack growth can be different.
The ideal situation would be to have a very refined mesh with a very small increment of crack growth or fatigue cycles, however, this is not computationally feasible. Each of these approaches is in place due to limitations in the amount of computational time available for solving problems. The algorithm presented remedies some of these issues such that a smaller increment of crack growth or fatigue cycles can be used during the analysis.
Algorithm3
Recall that the XFEM displacement approximation for a domain with a crack takes the following form.
As a crack grows in a domain, two main changes occur, the crack tip moves to a new element and elements which were previously enriched with the crack tip enrichment function are then enriched with the Heaviside function. Note from the above equation that the continuous displacement associated with the uI degrees of freedom is independent of the crack tip location. In addition, once an element is enriched with the Heaviside enrichment, it is for all future iterations independent of the crack tip location and non-changing. This leads to a convenient property when one considers the global XFEM stiffness matrix.
Thus, in effect for all iterations, Kuu, Kua and Kaa are constant. New Heaviside elements may be appended onto the end of Kua and Kaa and at each iteration of crack growth all stiffness terms associated with the crack tip enrichment Kub, Kab and Kbb must be recomputed. As the local region of crack growth may be very small this can in effect lead to approximately calculating the stiffness matrix in only ten elements per iteration, nine crack tip elements and a single Heaviside element which previously contained the crack tip.
Here, the Cholesky factorization from the initial iteration is modified directly through the use of a custom interface between CHOLMOD4 and MATLAB where the factorization of K and the row add and row delete operations as well as the solve of the triangular system of equations occur via MEX file.
Example
Consider an edge crack in a finite plate as shown in the figure below. For the case shown the a square domain of length 2 is meshed with square elements with sides of length 1/80. An edge crack of length 0.1 is located mid-height along the left edge of the domain. Fifteen iterations of crack growth are considered, with the incremental crack growth for each increment being equal to 0.1. The geometry and loading is shown on the left and the comparison of the needed computational time for assembly of the stiffness matrix with and without the reanalysis algorithm are provided. The assembly savings are instantaneous, while it takes a while to recover the cost of building and finding the fill-reducing ordering of the XFEM stiffness matrix where all enriched degrees of freedom are active.
Additionally, it was of interest to see how the cost of reinitialization depends upon the mesh density for both the reanalysis and factorization and solve. For the assembly the traditional trend is quadratic with respect to mesh density while the reanalysis algorithm is approximately linear. For the factoring and solving the traditional trend is cubic, while the reanalysis trend is approximately quadratic.
References
1. Moes, N., Dolbow, J., Belytschko, T. (1999) "A finite element method for crack growth without remeshing," International Journal for Numerical Methods in Engineering, 46, 131-150.
2. Gravouil, A., Moes, N., Belytschko, T., (2002) "Non-planar 3D crack growth by the extended finite element and level sets - Part II: level set update," International Journal for Numerical Methods in Engineering, 53, 2569-2586.
3. Pais, M., Kim, N.H., Davis, T. (2010) "Reanalysis of the extended finite element method for crack initiation and propagation," 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Orlando, Florida.
4. Davis, T. "CHOLMOD: Supernodal Sparse Cholesky Factorization and Update/Downdate," available at http://www.cise.ufl.edu/research/sparse/cholmod/.