Abstract

A new method is presented to study the interaction of multiple cracks, especially for the areas near crack tips by using the extended finite element method. In order to track the cracks, a new geometric tracking technique is proposed to track enriched elements and nodes along the crack instead of using the narrow band level set method. This allows to accurately determine enriched elements and nodes and calculate enrichment values. A method is proposed for constructing a multicrack matrix, which involves numbering enriched nodes of multiple cracks and solving the global stiffness matrix. In this approach, the stress fields around multiple cracks can be studied. The interaction integral method is employed to study the crack propagation and its direction by calculating the stress intensify factor. The developed model has been coded in MATLAB environment and validated against analytical solutions. The application of the model in the crack interaction study is demonstrated through a number of examples. The results illustrate the influence of the interaction of multiple cracks as they approach each other.

1. Introduction

Crack propagation can be efficiently simulated by resorting to the eXtended Finite Element Method (XFEM), saving a great deal of computational time and costs, given that no remeshing and refinement are performed [18]. Great advantages are gained over other methods such as the standard Finite Element Method (FEM) [9, 10], Element-Free Galerkin (EFG) method [11, 12], Boundary Element Method (BEM) [1315], and Discrete Element Method (DEM) [16, 17]. With XFEM, by enriching the nodes along the crack, additional degrees of freedom are given to reproduce the jump of the variables across the discontinuities. XFEM is associated with the Partition of Unity Method (PUM) [1821], and the shape functions for the enriched elements consist of standard shape functions, Heaviside step functions, and near-tip asymptotic functions [2].

In the context of XFEM, as a crack grows, the elements at the crack tips are sectioned and new elements containing a crack tip are formed. In order to deal with extending cracks, the narrow band level set method [22, 23] is used to track the crack path and to search for the advancing crack tips. This method can be sometimes inaccurate, especially when a crack forms a kink angle (see Section 2.1), so there is the need to find alternatives in order to overcome the shortcomings of the narrow band level set method.

Many researchers in the past decades studied the growth and mutual interaction of multiple cracks, mainly focusing on method developments and practical field applications [2426]. With reference to the mutual interaction of two cracks, Lawler [27] and Tanaka et al. [28] experimentally defined a relation between the crack propagation rate and the integral [29]. Carpinteri and Monetto [30] employed BEM to model the propagation of multiple cracks and implemented a global nonlinear stress-strain relationship associated with the geometry changes for the extension, intersection, and coalescence of the cracks. Budyn et al. [31] studied the growth of multiple cracks in brittle materials and presented a displacement equation for intersecting cracks. Daux et al. [32] proposed a theory for branched and intersecting cracks based on XFEM. Fageehi and Alshoaibi [33] studied nonplanar crack growth of multiple cracks by developing a code using FEM. Broumand [34] presented two methods, based on XFEM, to accurately detect multiple cracks of diverse size, shape, and orientation in 2D elastic bodies. Pham and Weijermars [35] used the linear superposition method to calculate stress tensor fields with multiple pressure-loaded fractures and further developed analytical solutions for the domains with large numbers of internally pressurized fractures and boreholes in a homogeneous elastic medium. Santoro et al. [36] discussed the response from a beam under static loads with multiple cracks and proposed an approach to evaluate the response in the presence of multiple cracks.

Even though an abundant literature is available concerning the numerical analysis of multiple cracks, the application of XFEM in this context is very limited. In this paper, a new method is proposed, based on XFEM, for the analysis of the propagation and interaction of multiple cracks in two-dimensional domains. The enriched nodes of the cracks are tracked by a proposed Geometric Tracking Technique (GTT). The enrichments for shape functions are ascribed to the enriched nodes of all fractures. Different rules for the Gauss points apply for the elements crossed by cracks and for those containing tips. An approach is presented to solve multiple cracks. In this approach, each crack in the domain is set as a unit prone to coalesce into the whole framework of the multiple cracks by combining all enriched stiffness matrices into one matrix. A node numbering rule is presented for multiple cracks. The interaction between cracks is based on domain forms of the Interaction Integral Method (IIM) [3740]. By studying the distancing and intersection of integral areas around tips, the crack behaviors can be predicted. The method is implemented into a MATLAB code. In what follows, the theoretical basis, steps of the analysis, and numerical implementation of the method are reported. In addition, numerical simulations for validation are illustrated to show the effectiveness and robustness of the method.

The paper is structured as follows: in Section 2, GTT is presented and an example is set up to verify the effectiveness with respect to the narrow band level set method; in Section 3, combining GTT and XFEM, the numerical results of two standard models are compared to corresponding analytical solutions; in Section 4, the simulation scheme for multiple cracks and the features of the XFEM solution are also presented; finally, in Section 5, the application of the method to multiple crack interaction and crack propagation prediction is illustrated. Concluding remarks are reported in Section 6.

2. Geometric Tracking Technique

In order to track the cracks in a domain, GTT is introduced for the accurate search and location of crack nodes and elements. The proposed technique can be used with different types of elements according to the given mesh and geometry. It is implemented in several steps. The first step involves of the selection of a rectangular area enclosing all the crack segments (each segment assumed linear), by checking the coordinates of the crack extremities (see Figure 1). The area can be represented by the set :

Each crack segment , enclosed by the subdomain , can be represented by a linear equation, as follows:

where , , and are the coefficients of the straight line of segment of the crack.

Subdomain can be represented by the set :

The enriched elements crossed by a crack can be located by finding the intersection between the crack linear segments (Equation (2)) and the boundaries of the elements. Each element has a given area and can be represented by the following set:

Substituting and into Equation (2) yields and :

If , or , is in the range of , the crack certainly intersects the element ; therefore, is an enriched element. Alternatively, by substituting and in the same equation, one obtains

If , or , is in the range of , element is enriched.

An example follows about the use of the technique for rectangular elements. With reference to Figure 2, where crack segment 1 and enclosing subdomain of Figure 1 are shown, subdomain can be represented by the set :

The equation of the straight line of crack segment 1, enclosed by subdomain , is

The set of element A is

where is the domain of element A. By substituting the coordinates of this element in Equation (2), the coordinate (P) of the intersection node P can be obtained; it is

Since the coordinates of node P are included by the set , it is proven that the crack segment 1 certainly goes through element A. Similarly, by substituting the coordinate of the other intersection node Q (still in Figure 2) in Equation (2), the relative coordinate is

It can be shown that both and are inside the domain of element B, so element B is also crossed by the crack segment 1; therefore, element B is an enriched element.

Note that GTT, in association with XFEM, could be potentially used also in Nonlinear Fracture Mechanics (NLFM) problems. In fact, the technique is apt to track crack paths following any geometric shape, irrespective of the cracking process (linear or nonlinear).

2.1. Comparison with the Narrow Band Level Set Method

In what follows, an example is given to explain the advantages in using GTT when dealing with cracks having large kink angles if compared to the narrow band level set method.

In the narrow band level set method, circular ranges around crack tips, supposedly containing all the enriched elements, are defined first [41]. However, within those circles, some of the elements inside may not intersect the crack, so a narrow band is introduced for further sifting the enriched elements. The distances between nodes and cracks can be calculated by using the nodal coordinates and the linear equations of the crack segments. Normally, the nodes within the range of the narrow band are regarded as target nodes and are saved in the set of the enriched nodes. The level set functions of the enriched nodes are used to ascertain if these nodes are located within the narrow band. With reference to Figure 3, where domain A and domain B encompass all the enriched elements, the linear equations of the two segments of the crack (Equation (2)) have coefficients 0.05, -1, and 1.7 (, , and , respectively) and 0.65, -1, and -1.122 (, , and , respectively) for segment 1 and segment 2, respectively.

The distance between node and crack segment 1 is

The intersection node N between the line and crack segment 2 has coordinates (6, 2.778). As 2.778 is between 2 and 3 and node N is on the edge of element D, the element D is an enriched element and node Q is an enriched node.

The distance between node and crack segment 2 is

According to the above calculations, node P and node Q have the same distance to crack segments 1 and 2, respectively. Since node Q is on the enriched element D, the distance from node Q to the crack should be equal to or less than the size of the narrow band. As node P has the same distance to the crack as node Q, the distance from the crack to node P is also equal to or less than the size of the narrow band. However, it is shown in Figure 3 that node P is not an enriched node. So, in this case, the application of the narrow band level set method produces an error.

Different from the narrow band level set method, the proposed GTT determines all the enriched elements that crack segments pass through by judging whether the element edges intersect with the crack segment. All the elements are examined for this quest with respect to each crack segment, so it is ensured that all the enriched elements are picked as target elements, no extra elements are included, and no enriched elements are ignored.

With reference to Figure 3, the intersection node M of line and crack segment 1 has the coordinates (2, 1.9). The coordinate of element C ranges from 2 to 3, so element C is not regarded as an enriched element; therefore, node P is not be counted as an enriched node. However, the range of the coordinate of the left edge of element F is from 1 to 2, including the coordinate of node M. So element F can be regarded as an enriched element, and the relative nodes are enriched nodes.

Based on the above comparison, one may state that the proposed GTT is accurate and effective and constitutes a robust procedure for tracking the crack growth trajectories.

3. Numerical Simulation of Single Cracks

With reference to the single crack in a solid domain in Figure 2, according to GTT, the elements A, B, C, and D can be part of a set of Heaviside-enriched elements (crack passing through) and tip-enriched elements (crack tips stay). This set contains element numbers and corresponding node numbers. In the next step, those elements are classified into two subsets. The elements including a crack tip are removed from the initial set (except the ones that contain crack tips bordering the edge of the domain) and then inserted in another set. The remaining elements of the initial set are then considered Heaviside-enriched elements. If an element includes both a Heaviside-enriched node and a tip-enriched node, it is denoted as a tip-enriched element or as the last Heaviside-enriched element before a tip-enriched element, so these two types of elements are differentiated. The enriched nodes normally include three types: Heaviside-enriched nodes, tip-enriched nodes, and mix-enriched nodes. Different types of enriched nodes are separately dealt with, according to the location of the elements.

After rearranging the Heaviside-enriched elements and the tip-enriched elements into different sets, an enrichment of the shape function can be adopted to represent the jump of displacements across the crack surfaces [42]: where is the step function that, for example, in the case of a horizontal crack, is equal to 1 when is greater than 0, and equal to -1 when is less than 0.

For tip-enriched nodes, the enrichment can be expressed as a tip branch function [2, 43]: where is the tip branch function () and and are the local polar coordinates at the crack tip.

The approximation for the displacement can be expressed by using the enrichments of Equations (14) and (15): where are the standard degrees of freedom (dofs); , , and are additional dofs; are the enrichment values of Heaviside-enriched nodes; are the enrichment values of nodes around crack tip 1; are the enrichment values of nodes around crack tip 2; are the standard shape functions; is the number of standard nodes; is the number of Heaviside-enriched nodes; and and are the sets of tip-enriched nodes for the first and second crack tips, respectively.

In order to demonstrate the validity of GTT combined with XFEM, in what follows, the solutions of two standard examples are reported. A comparison is made among the obtained values of the Stress Intensify Factor (SIF) with those descending from analytical solutions. IIM is used to calculate the SIF values.

3.1. Edge Crack

The first example refers to an -long single-edge crack (Figure 4(a)). It is horizontal and runs from the left vertical boundary of a rectangular plate of width () 100 m and height () 100 m, loaded by vertical tractions of 3.67 MPa on the top and bottom edges. Four-node quadrilateral elements are used. In order to examine the sensitivity of the model results to the element size, different structured meshes, all having equal-size square elements, are tested of , , and elements. The elastic modulus () and Poisson’s ratio () are 15 GPa and 0.25, respectively. Cracks having lengths 3.5, 5.5, 7.5, 9.5, and 11.5 m are considered. The shaded plot of the von Mises stress for the case with a 7.5 m long crack and elements is shown in Figure 4(b). The comparison between the numerical solution and analytical solution in terms of SIF values is reported in Table 1. An analytical solution for the Mode-I SIF is [44] where is the dominant vertical tensile stress (equal to the vertical tractions applied in the example), is the crack length, and is a correction coefficient, equal to

3.2. Central Crack

A plate with a central crack is loaded by an isotropic tensile stress of 3.67 MPa (Figure 5(a)). The discrete length values selected for the crack are 3, 7, 11, 15, and 19 m. In Figure 5(b), the shaded plot of the von Mises stress predicted by using XFEM for the 15 m long crack and the -element mesh is shown. The corresponding analytical solution for the Mode-I SIF is [45] with here the dominant isotropic stress (equal to the normal tractions applied in the example).

In Table 2, the comparison between the analytical solution and the numerical solution for different meshes is reported.

Observing the results in Tables 1 and 2, one may notice that when the mesh size is , the numerical solution is very close to the analytical solution, so there is no point in further refining the mesh, given that enough accuracy is gained. One may also state that with the proposed method, the crack growth propagation and direction for the two examples is effectively predicted.

4. Numerical Simulation of Multiple Cracks

For a multiple-crack problem, a specific numbering rule is applied (see Table 3). The nodes of a crack are reported in a column set from to ; column 1 contains the numbers of the standard nodes; is the number of total standard nodes; represents the number of the last enriched node in the last crack; and are the numbers of Heaviside-enriched nodes; and , , , and are the numbers of the tip-enriched nodes.

The governing equation for the whole domain can be written as where is the global load matrix, is the global displacement matrix of Equation (16), and is the global stiffness matrix, consisting of the standard stiffness matrix and the enriched stiffness matrices of the cracks; in formula, in which is the standard stiffness matrix; , ,…, are the enriched stiffness matrices of cracks 1 to , with being the number of cracks; and is the cross-term between the standard elements and the enriched elements. The matrices are as follows: where is the elasticity matrix, , , and are strain differential operator matrices, respectively defined as: with , , , and the tip branch functions around the crack tips and , , , and the enrichment values of the tip branch functions at the tip-enriched nodes.

A Delaunay triangulation is adopted herein to partition the enriched elements [46]. For the solutions of the integrals, the number of Gauss points is 9 for the standard elements, 3 for the triangular Heaviside-enriched elements, and 7 for the triangular tip-enriched elements.

5. Interaction of Multiple Cracks

The proposed method is used to analyze the interaction among multiple cracks through the following illustrative examples.

The stress field around a crack tip is influenced by other cracks in proximity, thus leading to stress-shadowing effects and redirection of the crack propagation [47]. SIF is useful to check if a crack propagates and to predict the propagation trajectory [4851]. IIM is often used to calculate SIF values since it has a high degree of accuracy compared with other methods [5255], although the involved process is quite complex; for the application, it is required to solve an energy integral, based on the integral method [29], covering a zone around the crack tip. Two states are considered: present state (1) and auxiliary state (2). If is the domain of integration (the zone around the tip), bounded by , for a local reference system of the crack, with origin at the tip and aligned with it, the interaction integral can be written as [2] where is the Dirac delta function, is the vector normal to , and is the interaction strain energy, equal to [1, 3]

The integral can also be expressed by using the values of SIF and , respectively, for Mode-I and Mode-II, for both the present state and the auxiliary state, as follows: where is Young’s modulus, and

Equation (24) can be transformed as follows: where is a weighting function. The solution of the integral in Equation (29) is relatively straightforward.

With reference to Figure 6(a), the radius of the zone is assumed equal to the square root of three times the element area [56]. The weighting function is as follows: where originates from the crack tip. When a portion of the zone is outside the domain, as in Figure 6(b), the weighting function () on the edge of the domain is assumed equal to zero [57].

In the following subsections, two examples of interaction are described, with two and five interacting cracks, respectively.

5.1. Two-Crack Example

In this example, two cracks, and , are considered. Domain, loading, material and mesh are the same used in Section 3.1. In a first case (case 1, Figure 7(a)), the cracks are parallel, both 5.5 m long and 25 m mutually distant. In a second case (case 2, Figure 7(b)), they are still parallel only 5 m distant. A third and a fourth case (case 3, case 4) with a slanted crack follow. For all the cases, four-node quadrilateral elements are employed. The vertical traction at the top and at the bottom boundaries is 3.67 MPa.

In case 1, the stress fields around the crack tips do not significantly interact. In case 2, the two stress fields overlap with mutual interference, producing a shadowing effect and causing an alteration of the SIFs resulting in the redirection of the crack propagation; in fact, given the mesh size, the zones for IIM overlap (Figure 7(b)) and Equation (29) for crack becomes and, for crack ,

The overlapped integral domain is then where is the intersection zone. For Equations (31), (32), and (33), subscripts and indicate crack and crack , respectively, based on Equation (24). The remaining equations follow the same notation rule, so the comprehensive interaction integral for crack is equal to the superposition of and : and the comprehensive interaction integral for crack is

Thus, by employing Equations (27) and (28), SIFs in Mode-I and Mode-II for the two cracks are

The results of the XFEM calculations in terms of the von Mises stress for two cases are shown in Figures 8 and 9. It can be seen that the cracks propagate horizontally (along the alignment of the cracks) only in case 1, whereas in case 2, the propagation directions deviate due to the mutual interaction. The results are in good agreement with those reported in the literature [58].

In Figures 10(a) and 10(b), the schemes of cases 3 and 4 are reported, respectively. There is a 5.5 m long horizontal crack (the “minor” crack, ) and a 10 m long slanted crack (inclined 45° with respect to the -axis) (the “major” crack, ). In case 3, the major crack is rather far from the minor crack (50 m from the tip of to the center of ). In case 4, the two cracks are closer. It is assumed that the propagation of crack is affected by the other crack when this one lies within the domain of the interaction integral of . The major crack can be seen as a border “arresting” the effect on the stress field of the minor crack , so the actual interaction integral domain for crack is equivalent to the entire circular area minus the shaded area (Figure 10(b)). The interaction integral for crack is as Equation (31).

The interaction integral of the shaded area of case 4 is (domain )

Thus, the final interaction integral of crack is equal to

The results for these two cases are shown in Figure 11. It can be seen that in case 3, the minor crack propagates along the initial horizontal direction, while in case 4, when the interaction integral domain includes the other crack, the propagation trajectory deviates from the horizontal along the shortest pathway towards the major crack, according to the results reported in the literature [3].

5.2. Multiple-Crack Example

In order to study the interaction of more complicated patterns with more than two cracks, two cases are considered in what follows. In case 1, there are five parallel edge cracks. Again, the domain, loading, material and mesh are the same as the previous example. The results in terms of von Mises stress are shown in Figure 12(a). It can be noticed that the cracks far from the horizontal middle line have a higher stress concentration and also show a larger deflection angle, while the cracks near to it have a smaller stress concentration and a mild deflection. Obviously, given the symmetry, the crack right along the middle line has no deflection. In front of the crack tips, the stress fields join to form a relatively large symmetric shadow area.

In case 2, five slanted cracks are randomly distributed. Similar to case 1, the farther the crack is from the middle line, the stronger the stress concentration is; some of the stress fields intersect each other to form a large shadow area (Figure 12(b)). However, the stress concentration area is smaller than that in case 1, for the weaker stress-shadowing effect, given the larger distances among the cracks. In fact, as the cracks propagate, they get closer to each other, and the shadowing effect starts to increase and become visible as shown in Figure 13. The consequential effect is that the shadowing areas extend and intersect with the other ones, and also, the values of the stress concentration increase.

6. Conclusion

This paper presents the development of a method for the analysis of propagation and interaction of multiple cracks. The method consists of the combination of a Geometric Tracking Technique (GTT) with an eXtended Finite Element Method (XFEM) formulation. The application of the method is illustrated by simulating three different examples, single crack, two cracks, and multiple cracks. The results of the first two examples are in good agreement with the corresponding solutions available in the literature. The third example is proposed to simulate the interaction among multiple cracks. The Interaction Integral Method (IIM) is employed to derive the Stress Intensify Factors (SIFs), useful for the prediction of the magnitude and direction of the crack propagation. From the results, one may derive that when the stress field around a crack tip is influenced by another crack in proximity, the propagation trajectory deviates off the crack alignment, thus producing a stress-shadowing effect. The closer the cracks are, the stronger this effect is. The patterns of the computed stress distributions well reflect the interaction among the cracks and as such are proofs of the good quality of the proposed method.

Data Availability

The data supporting this paper are from previously reported studies and datasets, which have been cited. The processed data are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no potential conflict of interests.

Acknowledgments

This study is financially supported by the Chinese Scholarship Council and University of Exeter.