Adaptive CT XIGA Using LR B-Splines for Efficient Fracture Modeling

TL;DR

This paper introduces a new method for simulating fractures in materials using advanced mathematical techniques to improve accuracy and efficiency.

Contribution

The novel adaptive CT XIGA framework using LR B-splines enables efficient fracture modeling without remeshing.

Findings

Adaptive CT XIGA achieves superior error convergence rates in fracture simulations.

The framework significantly reduces computational costs compared to existing methods.

Crack propagation simulations benefit from high accuracy and robustness using LR B-splines.

Abstract

This paper presents a novel adaptive crack-tip extended isogeometric analysis (adaptive CT XIGA) framework based on locally refined B-splines (LR B-splines) for efficient and accurate fracture modeling in two-dimensional solids. The XIGA method facilitates crack modeling without requiring the specific locations of crack faces and enables crack propagation simulation without remeshing by employing localized enrichment functions. LR B-splines, as an advanced extension of B-splines and NURBS, offer high-order continuity, precise geometric representation, and local refinement capabilities, thereby enhancing computational accuracy and efficiency. Various local mesh refinement strategies, designed based on crack and crack-tip locations, are investigated. Among these strategies, the crack-tip topological refinement strategy is adopted for local refinement in the adaptive CT XIGA framework. Stress intensity factors (SIFs) are evaluated using the contour interaction integral technique, while the maximum circumferential stress criterion is adopted to predict the crack growth direction. Numerical examples demonstrate the accuracy, efficiency, and robustness of adaptive CT XIGA. The results confirm that the proposed framework achieves superior error convergence rates and significantly reduces computational costs compared to a-posteriori-error-based adaptive XIGA methods, particularly in crack propagation simulations. These advantages establish adaptive CT XIGA as a powerful and efficient tool for addressing complex fracture problems in solid mechanics.

Figures

Full text

1. Introduction

Cracks in engineering structures have a significant impact on their load-bearing capacity, potentially compromising structural integrity and safety. To ensure the safe operation of such structures and to prevent catastrophic failures, gaining a comprehensive understanding of the mechanical behavior associated with fractures is essential. Numerical simulation has become an indispensable tool for investigating the mechanical behavior associated with fracture problems. Several numerical methods, such as the boundary element method (BEM) [1,2], meshfree methods (MMs) [3,4], phase-field method (PFM) [5,6], and extended finite element method (XFEM) [7,8], have been developed for crack problems. The XFEM is widely regarded as one of the most effective numerical approaches for fracture analysis. In the XFEM, standard finite element approximation is locally enriched by incorporating discontinuous displacement fields and the leading terms of the asymptotic crack-tip displacement fields, all within the framework of the partition of unity. These enrichment functions enable the XFEM to accurately capture the discontinuity across crack faces and the singularity near crack tips. As a result, the XFEM enables the entire crack to be represented independently of the computational mesh, thereby eliminating the need for remeshing during crack propagation simulations.

However, in the XFEM, geometric modeling relies on piecewise polynomial approximations, which introduce discretization errors that may reduce computational accuracy. The XFEM provides only $[eqn]$ continuity when analyzing fracture problems. To address the limitations of the XFEM, its concept has been integrated into isogeometric analysis (IGA) [9], leading to the development of extended isogeometric analysis (XIGA) [10]. The core idea of IGA is to employ the spline basis functions used in CAD for geometric representation as the shape functions in finite element analysis. This approach enables seamless integration between CAD and CAE. XIGA retains the advantages of the XFEM, including the ability to model cracks without needing to specify crack positions and to simulate crack propagation without remeshing. Additionally, it incorporates the benefits of IGA, such as higher-order continuity and precise geometric modeling. De Luycker et al. [11] employed XIGA for Mode-I crack analysis and achieved the optimal convergence rate using blending correction. Ghorashi et al. [12] were the first to use XIGA to simulate quasi-static crack propagation in 2D linear elasticity. Their results demonstrated excellent agreement between the computed SIFs and crack propagation paths in comparison to the analytical solutions and results obtained using the XFEM. Singh et al. [13] simulated three-dimensional crack problems in linear elastic materials with XIGA based on Bézier extraction. Zhong et al. [14] developed a three-dimensional rotating crack model using the XIGA method to perform vibration analysis of blades containing inclined single and multiple cracks. Khatir et al. [15] studied single- and multiple-crack identification through POD-RBF XIGA and the Jaya algorithm. Extended isogeometric analysis has also found applications in various fields, including thermo-mechanical fracture [16,17], fatigue crack growth [18,19,20], and orthotropic fracture materials [21,22,23]. Other notable applications include cohesive fracture [24,25], cracks in piezoelectric materials [26,27,28], T-stress evaluation [29,30], and optimization [31,32].

XIGA has found extensive applications in the study of fracture problems in plate and shell structures, attributed to its inherent high-order continuity. Nguyen-Thanh et al. [33] proposed an extended isogeometric analysis formulation to evaluate through-the-thickness cracks in thin shell structures based on Kirchhoff–Love theory. Nguyen et al. [34] investigated the vibration behavior of functionally graded microplates containing cracks. The study employed strain gradient theory with a single material length scale parameter and an additional micro-inertia term in conjunction with XIGA. Singh et al. [35] utilized XIGA to investigate the behavior of cracked functionally graded material (FGM) plates based on a generalized higher-order shear deformation theory. Yang et al. [36] integrated XIGA with the finite cell method to study the vibration and buckling behavior of FGM plates containing multiple internal cracks and cutouts under combined thermal and mechanical loads. Singh et al. [37] presented an extended isogeometric analysis for cracked functionally graded magneto-electro-elastic materials. Yadav et al. [38] conducted a pioneering study on the fracture behavior of carbon-nanotube-reinforced functionally graded structures with discontinuities under thermo-mechanical loading utilizing XIGA. Kumar et al. [39] conducted a comprehensive fracture analysis of edge-cracked porous functionally graded structures, considering various porosity distribution patterns under mechanical loading.

To accurately simulate fracture problems, fine elements are required in the vicinity of the crack. While using fine elements throughout the entire domain would ensure high computational accuracy, it would also result in excessively high computational costs. To improve computational efficiency without compromising accuracy, fine elements are employed near the crack region, while coarser elements are utilized in regions farther from the crack. Adaptive XIGA or multi-scale XIGA has been successfully applied to address fracture problems under such scenarios. Ghorashi et al. [40] were the first to apply multi-scale XIGA for the fracture analysis of orthotropic cracked media. Their approach utilized T-splines, which enable efficient local refinements to enhance computational accuracy and efficiency. Nguyen-Thanh et al. [41] utilized adaptive XIGA based on PHT-splines to investigate crack propagation in the vicinity of inclusions. Gu et al. [42,43,44] employed adaptive XIGA utilizing LR B-splines to model fracture problems and simulate crack propagation in both 2D isotropic and orthotropic materials. Yuan et al. [45] studied crack propagation in complex Mindlin–Reissner plates by developing an adaptive multi-patch XIGA according to LR NURBS and Nitsche’s method. The adaptive XIGA method, based on LR B-splines or LR NURBS, has also been applied to simulate arbitrary holes in orthotropic media [46], perform upper-bound limit analysis [47,48], and model cracked composite FG Mindlin–Reissner plates [49]. Schmidt et al. [50] developed the adaptive XIGA framework to predict the behavior of multi-material multiphysics problems with complex geometries, employing locally refined discretizations based on hierarchical B-splines. Jiang et al. [51] proposed adaptive XIGA with the strong imposition of essential boundary conditions for linear elastic fracture problems based on B++ splines. Qiu et al. [52] presented an adaptive XIGA method based on PHT-splines for solving phase-field fracture mechanics problems.

In our previous works [42,43,44], we developed an adaptive XIGA method based on LR B-splines for fracture analysis, focusing on the computation of SIFs and crack propagation. Adaptive local mesh refinement was achieved using a posteriori error estimation, in which the error was computed through recovery techniques. We observed that adaptive XIGA based on a posteriori error estimation facilitates fine refinement near the crack. In contrast, coarse-scale mesh refinement is applied in regions that are farther from the crack. This approach achieves high accuracy with reduced degrees of freedom, yielding an error convergence rate that is superior to global refinement. However, the application of a-posteriori-error-based adaptive XIGA to fracture problems, particularly crack propagation, entails significant computational costs. The primary reason lies in the computational demands of a-posteriori-error-based adaptive XIGA with local mesh refinement for calculating SIFs. This process requires multiple displacement solutions, along with stress or strain recovery, where each recovery step incurs a computational cost equivalent to that of solving a displacement solution. Furthermore, our previous studies revealed that adaptive local refinement predominantly occurs near the crack tip during crack propagation.

In this paper, we propose an adaptive CT XIGA framework based on LR B-splines for solving fracture problems in 2D elastic solids, including the calculation of stress intensity factors and the simulation of crack propagation. LR B-splines, as an extension of B-splines and NURBS, inherit precise geometric modeling capabilities and high-order continuity while enabling local mesh refinement. The adaptive CT XIGA framework incorporates the generalized Heaviside function and crack-tip enrichment functions to effectively capture the discontinuities and singularities induced by cracks. To identify regions for local refinement, we propose multiple mesh refinement strategies based on the positions of cracks and crack tips. We compare the error convergence rates and computational times of each local refinement strategy with those of both a-posteriori-error-based mesh refinement and uniform global mesh refinement. Based on these comparisons, we identify the most effective local refinement strategy, termed the crack-tip topological refinement strategy. The adaptive XIGA framework employing this strategy provides the following advantages: (i) crack modeling is independent of the crack position, and crack propagation can be simulated without remeshing. (ii) Compared to the XFEM, the proposed method achieves precise geometric modeling with higher accuracy and convergence rates. (iii) With local mesh refinement near the crack tip, the proposed method achieves higher accuracy at the same degrees of freedom compared to traditional XIGA. (iv) The crack-tip topological refinement strategy simplifies the implementation of local mesh refinement compared to a-posteriori-error-based methods. (v) Compared to a-posteriori-error-based adaptive XIGA, the proposed approach achieves higher error convergence rates while significantly reducing computational time, particularly in crack propagation simulations. The advantages of the proposed method are validated through numerical examples.

The structure of the paper is as follows: Section 2 begins with a brief overview of LR B-splines, B-splines, and NURBS. This is followed by the formulation of the XIGA approach applied to crack problems. Then, the computational formulas for SIFs and crack propagation are outlined. In Section 3, several numerical examples are provided to demonstrate the effectiveness of the proposed method. Finally, the conclusions of this study are presented in Section 4.

2. Adaptive XIGA Based on LR B-Splines for Cracked Media

2.1. LR B-Spline-Based Local Refinement

B-splines and non-uniform rational B-splines (NURBS), owing to their properties, including higher-order continuity and precise geometric representation, are widely used as basis functions in IGA and XIGA. However, their tensor-product structure poses challenges for performing local mesh refinement. LR B-splines were introduced by Dokken et al. [53] to overcome certain limitations associated with NURBS and B-splines. LR B-splines are a generalization of NURBS and B-splines, retaining their advantageous properties while enabling local refinement. The B-spline basis functions are defined using a knot vector $[eqn]$ , where $[eqn]$ are non-decreasing coordinates in the parameter space. Here, n denotes the number of basis functions, and p represents their degree. The ith B-spline basis function $[eqn]$ of degree p is recursively defined as follows:

$$[eqn]$$ $$[eqn]$$

NURBS are derived from B-splines and associated weights, enabling the exact representation of a wide variety of complex geometries. Given a B-spline basis function $[eqn]$ and weight $[eqn]$ , a univariate NURBS of degree p is constituted by

$$[eqn]$$

A bivariate NURBS is expressed by

$$[eqn]$$

where $[eqn]$ represents the weight; $[eqn]$ is the B-spline basis function of degree p in the $[eqn]$ direction, which is defined on the knot vector $[eqn]$ ; $[eqn]$ is the B-spline basis function of degree q in the $[eqn]$ direction, which is defined on the knot vector $[eqn]$ .

An LR B-spline basis function of degree p in $[eqn]$ direction is a weighted B-spline basis function with minimum support on an LR mesh and is expressed as follows:

$$[eqn]$$

where $[eqn]$ is a scalar weight factor ensuring that LR B-splines maintain the partition of unity property; $[eqn]$ denotes a B-spline basis function defined on a local knot vector $[eqn]$ . The set of local knot vectors $[eqn]$ with $[eqn]$ knots, extracted from the global knot vector $[eqn]$ , is defined as

$$[eqn]$$

A bivariate LR B-spline basis function of degree p in $[eqn]$ direction and degree q in $[eqn]$ direction is defined as follows:

$$[eqn]$$

Similar to the construction of NURBS using B-splines, LR NURBS can be constructed using LR B-splines, expressed as follows:

$$[eqn]$$

where $[eqn]$ denotes the number of control points.

Since LR NURBS can be formulated from LR B-splines and LR B-splines are a special case of LR NURBS where all weights $[eqn]$ are equal to 1, LR B-splines will be used to represent both LR B-splines and LR NURBS throughout the remainder of this text. For further details on LR B-splines and LR NURBS, the reader is referred to Ref. [53].

2.2. XIGA Formulation for Cracked Media

The XIGA is formed by incorporating the concepts of the XFEM into conventional IGA, enabling the simulation of cracks without considering the locations of cracks, and modeling crack propagation without remeshing. In XIGA, the isogeometric approximation of the displacement field is locally enriched to capture the discontinuity across crack faces and the singularity near crack tips. The XIGA approximation for a 2D crack problem is mathematically expressed as follows [42]:

$$[eqn]$$

with

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the generalized Heaviside function employed to capture the discontinuity across crack faces, which takes the value $[eqn]$ on one side of the crack and $[eqn]$ on the opposite side; $[eqn]$ denotes crack-tip enrichment functions to simulate singularity in the vicinity of crack tips; $[eqn]$ , $[eqn]$ , and $[eqn]$ represent the LR B-spline basis functions; $[eqn]$ , $[eqn]$ , and $[eqn]$ are the displacement control variables related to the standard IGA, the additional degrees of freedom (DOFs) associated with modeling crack faces, and the additional DOFs associated with crack tips, respectively; $[eqn]$ , $[eqn]$ , and $[eqn]$ denote the sets of LR B-splines in the whole discretization domain, supporting the elements containing crack faces (excluding crack tips), and whose support domain contains crack tips, respectively. The elements containing crack tips are referred to as crack-tip elements, while the elements containing crack faces (excluding crack tips) are referred to as crack-face elements. The control points whose corresponding basis functions are enriched by the Heaviside function and crack-tip enrichment functions are referred to as crack-face-enriched control points and crack-tip-enriched control points, respectively. These are illustrated in Figure 1b.

The displacement, strain, and stress at any point are written as follows:

$$[eqn]$$

where $[eqn]$ stands for the basis function matrix and $[eqn]$ is the strain–displacement matrix. For the details of $[eqn]$ and $[eqn]$ , the reader is referred to Ref. [42].

2.3. SIF Evaluation and Direction of Crack Growth

Fracture criteria and crack growth direction are critical topics in crack propagation research. Various criteria are employed to determine the direction of crack propagation. Considering efficiency and ease of implementation, the maximum circumferential stress criterion [54] is commonly used to predict the crack propagation direction. According to the maximum circumferential stress criterion, the crack growth direction $[eqn]$ is written as

$$[eqn]$$

In Equation (14), $[eqn]$ and $[eqn]$ denote the Mode-I and Mode-II SIFs, respectively. The interaction integral technique is applied to compute the mixed-mode SIFs. The formulation for mixed-mode SIFs is provided by [51]

$$[eqn]$$

where

$$[eqn]$$

and

$$[eqn]$$

In Equation (16), $[eqn]$ and $[eqn]$ represent the actual and auxiliary states, respectively. A denotes the area enclosed by the contour. The function q is defined as 1 within the integration area A and 0 outside. $[eqn]$ represents the interaction strain energy. In Equation (17), E and $[eqn]$ are the elastic modulus and Poisson’s ratio, respectively.

3. Numerical Examples

To validate the accuracy and efficiency of the methodology proposed herein, six numerical examples are presented in this section. The first three examples focus on the calculation of SIFs, while the subsequent three examples simulate crack propagation. The validity and effectiveness of the adaptive extended isogeometric analysis, based on crack-tip locations for addressing fracture issues, are confirmed by comparing the computed results with exact or reference solutions. This study investigates quasi-static crack propagation under the assumption that all cracks are assumed to propagate freely. Unless otherwise specified, all examples in this section adopt LR B-spline basis functions with third-order polynomial directions ( $[eqn]$ ), consistent with Refs. [42,43]. The red line represents the crack in the figures below. Plane strain conditions are assumed. Convergence curves of error are plotted on a natural logarithmic scale $[eqn]$ , and seconds are used as the time unit for calculations. The structured mesh refinement strategy [55] is employed to locally refine the LR mesh. The numerical integration is conducted in a manner consistent with previous works [42,43]. The adaptive extended isogeometric analysis program in this study is implemented using a hybrid MATLAB R2023b–C++ 11.4.0 programming approach. All simulations are performed on a single workstation equipped with a 14th-generation Intel Core i9 processor (3.2 GHz) and 64 GB of memory. To ensure optimal performance and accuracy, the workstation is exclusively allocated for running the program, with no other tasks handled during execution.

3.1. Infinite Plate with a Center Crack

An infinite plate with a central straight crack, subjected to a unidirectional uniform tensile force of $[eqn]$ units, is illustrated in Figure 1a. The crack length is defined as $[eqn]$ units. To effectively deal with the infinite domain, a square region $[eqn]$ enclosing the crack is selected for analysis. The square has a side length of 10 units, with the crack having a length of $[eqn]$ units. A coordinate system is established with vertex A as the origin, where the direction $[eqn]$ serves as the x-axis and $[eqn]$ serves as the y-axis. The coordinates of the crack tip are specified as $[eqn]$ . The material properties include elastic modulus of $[eqn]$ units and Poisson’s ratio $[eqn]$ . Using a local polar coordinate system centered at the crack tip, the analytical solutions for the displacement and stress fields are provided as follows [42]:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ denotes the polar coordinate system centered at the crack tip; $[eqn]$ denotes the Model-I SIF.

In the numerical analysis, exact displacement boundary conditions are imposed on the top, bottom, and right edges of the square $[eqn]$ in accordance with the analytical solution for displacements outlined in Equation (18). Exact stress boundary conditions are applied to the left edge of the square $[eqn]$ according to the analytical stress solution presented in Equation (19). The initial computational mesh consists of $[eqn]$ elements, as illustrated in Figure 1b.

Firstly, the localized refinement regions around the crack are investigated. As the structured mesh refinement strategy is employed, the focus is on determining which LR B-spline basis functions require refinement. Two main categories are considered as follows:

1. Refinement is applied only to the localized region near the crack tip, ignoring the crack-face region.
2. Refinement is applied both to the localized region near the crack tip and the area surrounding the crack face.

For mesh refinement in the localized region around the crack tip, two scenarios are considered as follows:

• Refinement is applied to the LR B-spline basis functions supporting the crack-tip element, referred to as the crack-tip topological refinement strategy.
• A refinement radius $[eqn]$ , where $[eqn]$ and $[eqn]$ represent the area of the crack-tip element. Refinement is applied to the LR B-spline basis functions supporting elements within the circular region of radius $[eqn]$ centered at the crack tip. This approach is referred to as the crack-tip geometric refinement strategy. The crack-tip topological refinement strategy can be considered a special case of the crack-tip geometric refinement strategy where $[eqn]$ . This study specifically examines cases where $[eqn]$ or 3. Regions with $[eqn]$ , being farther from the crack tip, are excluded from refinement.

For the crack-face region, refinement is utilized for all non-zero LR B-spline basis functions supporting the crack-face elements. The various refinement strategies are classified and analyzed. Using these refinement strategies in accordance with the structured mesh refinement strategy, the first three steps of adaptive local refinement meshes are obtained, as shown in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7. These strategies are classified into six main cases:

Case 1: The crack-tip topological refinement strategy. Adaptive XIGA with this strategy is referred to as adaptive crack-tip XIGA (adaptive CT XIGA), as shown in Figure 2.

Case 2: The crack-tip geometric refinement strategy with $[eqn]$ . Adaptive XIGA with this strategy is referred to as adaptive crack-tip 2 XIGA (adaptive CT2 XIGA), as shown in Figure 3.

Case 3: The crack-tip geometric refinement strategy with $[eqn]$ . Adaptive XIGA with this strategy is referred to as adaptive crack-tip 3 XIGA (adaptive CT3 XIGA), as shown in Figure 4.

Case 4: A combined strategy involving crack-tip topological refinement and crack-face refinement. Adaptive XIGA with this strategy is referred to as adaptive crack-tip and crack-face XIGA (adaptive CTCF XIGA), as shown in Figure 5.

Case 5: A combined strategy involving crack-tip geometric refinement with $[eqn]$ and crack-face refinement. Adaptive XIGA with this strategy is referred to as adaptive crack-tip 2 and crack-face XIGA (adaptive CT2CF XIGA), as shown in Figure 6.

Case 6: A combined strategy involving crack-tip geometric refinement with $[eqn]$ and crack-face refinement. Adaptive XIGA with this strategy is referred to as adaptive crack-tip 3 and crack-face XIGA (adaptive CT3CF XIGA), as shown in Figure 7.

Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 show that, for all refinement strategies, both the mesh sizes around the crack tip and the refinement zones decrease as the number of refinement steps increases. For the crack-tip geometric refinement strategy, the refinement zone expands as $[eqn]$ increases. Among all the strategies, the crack-tip topological refinement strategy yields the smallest refinement zone.

Then, the convergence rate of relative errors in the H1 and energy norms, as a function of the number of DOFs, is analyzed for adaptive XIGA under various refinement strategies, as displayed in Figure 8 and Figure 9. The formulations for the relative errors in the H1 and energy norms are provided in Ref. [42]. For comparison, the results of the XIGA with adaptive local refinement based on a posteriori error estimation and uniform global refinement [42] are also presented. The adaptive XIGA approach using local refinement based on a posteriori error estimation with the recovery technique proposed by Zienkiewicz and Zhu [56] is called adaptive ZZ XIGA, while the XIGA with uniform global refinement is called uniform XIGA. All mesh refinement strategies exhibit a decrease in relative errors as the number of DOFs increases, converging to smaller values after three mesh refinement steps. Notably, the natural logarithm values for the relative errors in H1 norm and energy norm stabilize around $[eqn]$ and $[eqn]$ , respectively, while their actual values approach 0.111% and 0.691%. Such minimal errors underscore the accuracy and efficiency of adaptive XIGA based on LR B-splines under various refinement strategies for addressing fracture problems.

Among the various strategies examined, the error convergence rate of the crack-tip mesh refinement strategy is generally superior to that of the combined strategy involving both crack-tip and crack surface mesh refinement. It is also faster than that of uniform global refinement and a-posteriori-error-based adaptive mesh refinement strategies, as shown in Figure 8 and Figure 9. Furthermore, the crack-tip topological mesh refinement strategy exhibits the highest error convergence rate among all the refinement strategies. The adaptive XIGA method adopts the topological mesh refinement strategy, performing two and three local mesh refinement steps, respectively. This approach achieves relative errors of $[eqn]$ and $[eqn]$ for Mode-I SIF, demonstrating its high accuracy and efficiency.

Finally, the computation time required to determine the Mode-I SIF at different refinement steps derived from adaptive XIGA with various refinement strategies is discussed, as illustrated in Figure 10. For all the refinement strategies, the computation time increases as the number of refinement steps increases. However, the computation time growth rate of adaptive XIGA with six refinement strategies proposed in this paper is significantly lower than that of uniform XIGA and adaptive ZZ XIGA. In particular, adaptive CT XIGA shows the slowest growth rate in computation time and has the shortest computation time at the 1st, 2nd, and 3rd refinement steps among all the strategies. The computation times for adaptive CT XIGA to calculate the Mode-I SIF at the 1st, 2nd, and 3rd refinement steps are 0.75 s, 0.95 s, and 1.13 s, respectively. It is noteworthy that the longer computation time for uniform XIGA results from the increased number of mesh elements and degrees of freedom caused by global uniform refinement. Although adaptive ZZ XIGA can reduce the number of elements and degrees of freedom to some extent, the computational cost remains high due to multiple recovery steps and displacement field calculations [42]. The computation times for uniform XIGA and adaptive ZZ XIGA to calculate the SIF at the 3rd refinement step are 15.36 s and 11.33 s, respectively, which are 11.33 times and 10.04 times longer than that of adaptive CT XIGA. Due to its high convergence rate and low computational cost, adaptive CT XIGA is adopted for the remaining numerical examples.

3.2. Circular Plate with a Central Crack

To illustrate the capability of the adaptive CT XIGA method based on LR B-splines in accurately representing geometric structures with curved boundaries, this example considers a circular plate with a central crack. The radius of the plate is $[eqn]$ unit, and the crack length is $[eqn]$ units, passing through the center, as depicted in Figure 11. A constant traction force of $[eqn]$ units is applied along the circumferential normal direction. The elastic modulus is $[eqn]$ units, and the Poisson’s ratio is $[eqn]$ . This example focuses on the analysis of Mode-I SIF, with the theoretical solution referenced in [57]

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ represents the ratio of the crack length to the diameter of the circular plate and is expressed as $[eqn]$ .

The adaptive XIGA approach adopts the crack-tip topological refinement strategy and structured mesh refinement strategy to locally refine the mesh for solving the Mode-I SIF in this example. The initial computational mesh, which consists of $[eqn]$ elements, and the first three steps of refined meshes are presented in Figure 12. The geometry of the circular plate is accurately modeled by adaptive CT XIGA based on LR B-splines. As expected, local refinement primarily occurs near the two crack tips. As the number of refinement steps increases, the local refinement regions near the two crack tips become progressively smaller and the element sizes in the vicinity of the two crack tips continue to decrease.

The convergence curve of the relative error of the Mode-I SIF versus the number of DOFs obtained using adaptive CT XIGA is presented in Figure 13, while the convergence curves corresponding to the uniform XIGA and adaptive ZZ XIGA methods [42] are also included for comparison. These results clearly demonstrate the superior performance of adaptive CT XIGA, achieving a significantly lower relative error with fewer DOFs compared to the other approaches. Figure 14 illustrates the computational time required to solve the SIFs for the two crack tips at different refinement steps using three methods: uniform XIGA, adaptive ZZ XIGA, and adaptive CT XIGA. The results reveal that adaptive CT XIGA exhibits a much lower computational cost compared to uniform XIGA and adaptive ZZ XIGA, particularly with increasing refinement steps. This underscores the superior efficiency of adaptive CT XIGA in terms of computational performance.

To assess the performance of adaptive CT XIGA in addressing the Mode-I SIF, Mode-I SIFs derived from adaptive CT XIGA are compared with exact solutions for various crack lengths, specifically $[eqn]$ = 0.04 m, 0.06 m, 0.08 m, 0.10 m, 0.12 m, and 0.14 m, as illustrated in Figure 15. There is remarkable agreement between the computed values and their exact solutions, thereby demonstrating that adaptive CT XIGA exhibits efficiency in resolving fracture problems.

3.3. Square Plate with a Center Curved Crack

To demonstrate the capability of the proposed method for simulating both curved cracks and small cracks, a square plate with a side length of $[eqn]$ units containing a central curved crack is analyzed, as depicted in Figure 16. The plate has a large edge-to-crack length ratio (>10), indicating a relatively small crack size compared to the overall model. The central curved crack is modeled as a circular arc with endpoints at $[eqn]$ and $[eqn]$ . The crack geometry is further characterized by a radius $[eqn]$ units and a subtended angle $[eqn]$ . The elastic modulus is specified as $[eqn]$ units, while the Poisson’s ratio is set to $[eqn]$ . The top and bottom edges of the plate are subjected to uniform tensile stress along the y-direction with a magnitude of $[eqn]$ unit. For the central curved crack in an infinite plate, the analytical solutions for the mixed-mode SIFs are provided by [42]

$$[eqn]$$

where $[eqn]$ .

The convergence behavior for the mixed-mode SIFs is studied by setting $[eqn]$ . The initial computational mesh, consisting of $[eqn]$ elements, is presented in Figure 17, along with the locally refined meshes of adaptive CT XIGA at the first, third, and fifth refinement steps. Due to the small size of the curved crack, the local refinements during the first three steps are focused around the crack face, while the fourth and fifth refinement steps predominantly occur near the crack tips.

Figure 18 illustrates the convergence behavior of the mixed-mode SIFs, $[eqn]$ and $[eqn]$ , with respect to the degrees of freedom. The results achieved with adaptive CT XIGA are compared with those from uniform XIGA and adaptive ZZ XIGA [42]. All the methods demonstrate convergence in computing the mixed-mode SIFs. Specifically, the Mode-I SIF $[eqn]$ converges after the fourth mesh refinement, whereas the Mode-II SIF $[eqn]$ requires five mesh refinements to converge. Among the three methods, adaptive CT XIGA achieves the fastest convergence rate. When calculating the mixed-mode SIFs using the interaction integral method, sufficiently small mesh sizes near the crack tips are crucial, particularly for small cracks. While all three methods—adaptive ZZ XIGA, uniform XIGA, and adaptive CT XIGA—can adequately refine the mesh near the crack tips, their computational time requirements differ significantly. Figure 19 illustrates the computational time required by each method at different refinement steps to calculate the mixed-mode SIFs. For four refinement steps, uniform XIGA requires substantially more computational time than adaptive ZZ XIGA and adaptive CT XIGA. Moreover, adaptive CT XIGA shows much lower computational time and a slower increase in computational time per refinement step compared to the other two methods. This highlights the advantage of adaptive CT XIGA in terms of computational efficiency when analyzing small cracks.

To evaluate the performance of adaptive CT XIGA in modeling curved cracks, the mixed-mode SIFs are investigated for various values of $[eqn]$ . As $[eqn]$ increases, the Mode-I SIF ( $[eqn]$ ) exhibits a decreasing trend, whereas the Mode-II SIF ( $[eqn]$ ) shows a corresponding increase. The results computed by adaptive CT XIGA with five local refinements for the mixed-mode SIFs are in very good agreement with the exact solutions, as shown in Figure 20. This agreement highlights the robustness and accuracy of adaptive CT XIGA in resolving the mixed-mode SIFs associated with a curved crack.

3.4. Cantilever Beam with an Edge Crack

To demonstrate the applicability of adaptive CT XIGA to crack propagation simulation, this example investigates the propagation of an edge crack in a cantilever beam. The geometric configuration of the cantilever beam, with dimensions $[eqn]$ units and $[eqn]$ units, is illustrated in Figure 21. A planar Cartesian coordinate system is defined, with its origin at the bottom-left corner of the beam. An edge crack of length $[eqn]$ units is positioned slightly above the mid-plane of the left boundary, offset by $[eqn]$ units. The beam is subjected to concentrated vertical loads of $[eqn]$ unit, applied vertically at its top-left and bottom-left corners, while the beam’s right boundary is fixed. The material properties are defined under plane stress conditions, with the elastic modulus of $[eqn]$ units and the Poisson’s ratio of $[eqn]$ .

The adaptive CT XIGA approach based on LR B-splines is employed to simulate edge crack propagation in the cantilever beam. The initial mesh consists of $[eqn]$ uniform elements, and the crack increment length is set to $[eqn]$ units, consistent with Ref. [43]. The total number of crack propagation steps is 10. According to Ref. [43], to ensure accuracy in the crack propagation path, the interaction integral radius $[eqn]$ should be smaller than the crack increment length $[eqn]$ . Therefore, adaptive CT XIGA requires at least three local refinement steps.

Figure 22 illustrates the locally refined meshes during crack propagation at step 0, step 4, step 7, and step 9. The figure shows that, as the crack propagates, local mesh refinement primarily occurs in the vicinity of the crack tip, ensuring that the element size near the crack tip remains sufficiently small. Starting from step 4 (Figure 22b), the crack propagation direction begins to deviate from its initial trajectory, extending upward from the horizontal axis. This behavior is consistent with the actual propagation scenario, where the initial crack is slightly above the centerline of the cantilever beam and is subjected to concentrated loading. By step 7 (Figure 22c), the deviation of the crack propagation direction from its initial trajectory becomes more pronounced. Additionally, at step 9 (Figure 22d), the crack has propagated near the model boundary, where further propagation and refinement are significantly affected by boundary effects. This is one reason that the number of propagation steps is limited to 10, with no further extension considered.

Figure 23 presents the crack propagation paths achieved via uniform XIGA, adaptive ZZ XIGA, and adaptive CT XIGA with three refinement steps, alongside the crack propagation path computed by the XFEM with 4800 elements [58]. The close agreement among the crack propagation paths generated by these methods highlights the effectiveness of adaptive CT XIGA in accurately simulating edge crack propagation. However, as illustrated in Figure 24, the computational time required to simulate crack propagation differs significantly among these methods.

Figure 24 shows the computational time for each crack propagation step using uniform XIGA, adaptive ZZ XIGA, and adaptive CT XIGA with three refinement steps. The computational times for uniform XIGA and adaptive ZZ XIGA are comparable. In contrast, adaptive CT XIGA requires substantially less computational time per propagation step than both uniform XIGA and adaptive ZZ XIGA. The primary reason for this difference is attributed to the computational time required for SIF evaluation. Adaptive CT XIGA significantly reduces the computational cost of SIF evaluation compared to uniform XIGA and adaptive ZZ XIGA, as shown in earlier examples. The smaller the crack increment length $[eqn]$ , the smoother the crack propagation path becomes. However, since the crack increment length $[eqn]$ must be greater than the interaction integral radius $[eqn]$ , the element size near the crack tip must be sufficiently small. If fine elements are used throughout the entire domain, as is the case with uniform XIGA, the computational cost increases substantially, especially when the crack size is relatively small compared to the overall model or when the crack increment length is set to a small value. A more computationally efficient approach involves using fine elements near the crack tip and coarser elements in regions farther from the crack tip, as implemented in adaptive ZZ XIGA and adaptive CT XIGA.

However, adaptive ZZ XIGA requires multiple recovery steps and displacement variable computations to calculate the SIFs, leading to substantially higher computational time compared to adaptive CT XIGA. Consequently, adaptive CT XIGA is highly suitable for simulating crack propagation, delivering both superior accuracy and computational efficiency.

3.5. Square Plate with a Center-Inclined Crack

To further demonstrate the effectiveness of the adaptive CT XIGA method in simulating crack propagation, this example investigates the propagation of a centrally inclined crack in a square plate under uniaxial tension, as shown in Figure 25. The length of the square plate is $[eqn]$ units, and the crack length is $[eqn]$ unit. The angle of inclination of the crack relative to the horizontal direction is represented by $[eqn]$ . The material parameters are $[eqn]$ units and $[eqn]$ . A uniaxial uniform tensile stress of $[eqn]$ unit is applied to the top and bottom edges of the square plate in the y-direction.

First, the crack inclination angle $[eqn]$ is set to 45, with a crack growth increment $[eqn]$ of $[eqn]$ units. The crack propagates over a total of three steps for both crack tips, in accordance with Ref. [58]. The initial computational mesh consists of $[eqn]$ elements. To ensure that the interaction integral radius, $[eqn]$ , remains smaller than the crack growth increment $[eqn]$ , at least three local refinement steps are required for adaptive CT XIGA. The meshes generated by adaptive CT XIGA during the first three crack growth steps are presented in Figure 26. As the crack propagates, local mesh refinement consistently occurs in the vicinity of the two crack tips.

Figure 27 presents the crack growth paths obtained using uniform XIGA, adaptive ZZ XIGA, adaptive CT XIGA with three refinement steps, and the XFEM method [58]. The results demonstrate strong agreement between the crack growth paths produced by these methods, highlighting the effectiveness and accuracy of adaptive CT XIGA in simulating inclined crack propagation. A comparison of the computation times for each crack propagation step using uniform XIGA, adaptive ZZ XIGA, and adaptive CT XIGA with three refinement steps is illustrated in Figure 28. It is evident that both the time required for each propagation step and the time increase rate indicate that adaptive CT XIGA outperforms the other methods.

Then, the crack growth increment remains unchanged, and the total number of crack propagation steps for both tips is set to 6. The crack propagation paths for different crack inclination angles, obtained using adaptive CT XIGA with three refinement steps, are shown in Figure 29. For each crack inclination angle, the crack tip propagation directions are symmetric, and, in all cases, the cracks ultimately propagate in the horizontal direction.

3.6. Square Plate with Two Edge Cracks

To evaluate the proposed method’s effectiveness in simulating multi-crack propagation, the final numerical example analyzes a representative example of double-crack propagation. Figure 30 illustrates a square plate with a side length of $[eqn]$ units. The plate features two edge cracks of length $[eqn]$ units, positioned symmetrically on the left and right edges. The material properties of the plate are defined by the elastic modulus $[eqn]$ unit and the Poisson’s ratio $[eqn]$ . Vertical displacements, $[eqn]$ unit, are applied along both the top and bottom boundaries of the plate.

The initial computational mesh is set to $[eqn]$ elements, and the crack growth increment is selected as $[eqn]$ , with a total of 18 crack propagation steps, consistent with Ref. [43]. To ensure that the interaction integral radius $[eqn]$ does not exceed the crack growth increment $[eqn]$ , adaptive CT XIGA requires at least four local refinement steps. Figure 31 presents meshes of adaptive CT XIGA for the crack growth at the initial stage, the 8th step, and the 17th step. As expected, with crack propagation, local mesh refinement predominantly occurs near the two crack tips, ensuring sufficiently small element sizes.

The crack growth paths simulated using uniform XIGA, adaptive ZZ XIGA, adaptive CT XIGA with four refinement steps, and the SGBEM method [59] are presented in Figure 32. The results indicate that the crack growth paths obtained by these methods show strong agreement, highlighting the effectiveness and accuracy of adaptive CT XIGA in simulating multiple-crack propagation. Similar to the single-crack propagation simulation, the computational time for uniform XIGA, adaptive ZZ XIGA, and adaptive CT XIGA differs significantly, as illustrated in Figure 33. Because of the use of four refinement steps, uniform XIGA involves a large number of elements (up to 65536), resulting in significantly higher computational time per propagation step than adaptive ZZ XIGA and adaptive CT XIGA. Furthermore, compared to adaptive ZZ XIGA, adaptive CT XIGA requires significantly less time per propagation step and for the entire 18 steps. This highlights the significant computational efficiency advantage of adaptive CT XIGA in simulating multiple-crack propagation.

4. Conclusions

In this study, an adaptive CT XIGA framework based on LR B-splines has been presented for efficient and accurate fracture modeling in two-dimensional solids. By utilizing localized enrichment functions, the proposed method enables the representation of cracks without requiring the computational mesh to conform to the crack geometry, thereby eliminating the need for remeshing during crack propagation. The integration of LR B-splines, which offer high-order continuity and flexible local refinement capabilities, greatly enhances both the computational accuracy and efficiency of fracture simulations. Several local mesh refinement strategies have been designed and developed based on the positions of cracks and crack tips. The error convergence rates and computational times of different refinement strategies have been systematically compared. Among these, the crack-tip topological refinement strategy has been identified as the most effective strategy due to its ability to achieve the highest convergence rate while demonstrating the lowest computational cost. Numerical examples have been presented to validate the accuracy and effectiveness of the adaptive CT XIGA method in computing SIFs and simulating crack propagation. Furthermore, compared to a-posteriori-error-based adaptive XIGA, the adaptive CT XIGA approach exhibits superior error convergence rates and achieves substantial reductions in computational time, making it especially advantageous for crack propagation simulations. In addition, the crack-tip topological refinement strategy offers a simpler and more streamlined implementation of local mesh refinement relative to a-posteriori-error-based methods. By integrating the geometric precision of isogeometric analysis with concepts inspired by the XFEM and leveraging the efficiency of local refinement, the proposed method offers a robust and computationally efficient framework for addressing complex fracture problems in solid mechanics. It has potential applications in crack growth modeling for plates, shells, and three-dimensional structures.