# Analysis of Interaction of Multiple Cracks Based on Tip Stress Field Using Extended Finite Element Method

Authors: Yuxiao Wang; Akbar A. Javadi; Corrado Fidelibus
Journal: Journal of Applied Mathematics
 (2022)

---

## 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.


---

## Body

## 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 [1–8]. 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) [13–15], 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) [18–21], 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 [24–26]. 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 J 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) [37–40]. 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 Section2, 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Ωw enclosing all the crack segments (each segment assumed linear), by checking the coordinates of the crack extremities (see Figure 1). The area Ωw can be represented by the set S:
(1)S=x,y∈Ωwxs1≤x≤xs2,fxs1≤y≤fxs2.Figure 1
Sketch of a simulation domain withΩw, Ωsubi, and segments i and i+1; the equation of the line of segment i is Aix+Biy+Ci=0.Each crack segmenti, enclosed by the subdomain Ωsubi, can be represented by a linear equation, as follows:
(2)Aix+Biy+Ci=0,whereAi, Bi, and Ci are the coefficients of the straight line of segment i of the crack.SubdomainΩsubi can be represented by the set Ssubi:
(3)Ssubi=x,y∈Ωsubixsubm≤x≤xsubn,fxsubm≤y≤fxsubn.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 e has a given area and can be represented by the following set:
(4)Γe=x,y∈Ωexe1≤x≤xe2,fxe1≤fx≤fxe2.Substitutingxe1 and xe2 into Equation (2) yields ye1 and ye2:
(5)ye1=−Ci−Aixe1Bi,ye2=−Ci−Aixe2Bi.Ifye1, or ye2, is in the range of fxe1,fxe2, the crack certainly intersects the element e; therefore, e is an enriched element. Alternatively, by substituting fxe1 and fxe2 in the same equation, one obtains
(6)xe1′=−Ci−Aifxe1Bi,xe2′=−Ci−Aifxe1Bi.Ifxe1′, or xe2′, is in the range of xe1,xe2, element e is enriched.An example follows about the use of the technique for rectangular elements. With reference to Figure2, where crack segment 1 and enclosing subdomain Ωsub1 of Figure 1 are shown, subdomain Ωsub1 can be represented by the set S1:
(7)S1=x,y∈Ωsub1x1≤x≤x3,fx1≤y≤fx3.Figure 2
Searching for elements intersecting crack segment 1.The equation of the straight line of crack segment 1, enclosed by subdomainΩsub1, is
(8)A1x+B1y+C1=0.The setΓeA of element A is
(9)ΓeA=x,y∈ΩeAx11≤x≤x22,fx11≤fx≤fx22,whereΩeA is the domain of element A. By substituting the x coordinates of this element in Equation (2), the coordinate y (P) of the intersection node P can be obtained; it is
(10)yP=−C1−A1x11B1.Since the coordinatesx11,yP of node P are included by the set ΓeA, it is proven that the crack segment 1 certainly goes through element A. Similarly, by substituting the y coordinate of the other intersection node Q (still in Figure 2) in Equation (2), the relative x coordinate is
(11)xQ=−C1−B1fx22A1.It can be shown that bothxQ and yQ 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 (A1, B1, and C1, respectively) and 0.65, -1, and -1.122 (A2, B2, and C2, respectively) for segment 1 and segment 2, respectively.Figure 3
Sketch for the comparison of the proposed tracking technique and the narrow band method.The distanced1 between node P2,3 and crack segment 1 is
(12)d1=0.05×2−1×3+1.70.052+1=1.199.The intersection node N between the linex=6 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 distanced2 between node Q7,2 and crack segment 2 is
(13)d2=0.65×7−1×2−1.1220.652+1=1.199.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 Figure3 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 Figure3, the intersection node M of line x=2 and crack segment 1 has the coordinates (2, 1.9). The y 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 y coordinate of the left edge of element F is from 1 to 2, including the y 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.


## 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 (A1, B1, and C1, respectively) and 0.65, -1, and -1.122 (A2, B2, and C2, respectively) for segment 1 and segment 2, respectively.Figure 3
Sketch for the comparison of the proposed tracking technique and the narrow band method.The distanced1 between node P2,3 and crack segment 1 is
(12)d1=0.05×2−1×3+1.70.052+1=1.199.The intersection node N between the linex=6 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 distanced2 between node Q7,2 and crack segment 2 is
(13)d2=0.65×7−1×2−1.1220.652+1=1.199.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 Figure3 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 Figure3, the intersection node M of line x=2 and crack segment 1 has the coordinates (2, 1.9). The y 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 y coordinate of the left edge of element F is from 1 to 2, including the y 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 Figure2, 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φj of the shape function Nj can be adopted to represent the jump of displacements across the crack surfaces [42]:
(14)φjx,y=Hx,y,where Hx,y is the step function that, for example, in the case of a horizontal crack, is equal to 1 when y is greater than 0, and equal to -1 when y is less than 0.For tip-enriched nodes, the enrichmentφl can be expressed as a tip branch function [2, 43]:
(15)φlr,θ=Blr,θ=rsinθ2,rsinθ2sinθ,rcosθ2,rcosθ2sinθ,where Bl is the tip branch function (l=1,2,3,4) and r and θ are the local polar coordinates at the crack tip.The approximation for the displacementu can be expressed by using the enrichments of Equations (14) and (15):
(16)ux,y=∑i=1nsNix,yui+∑j=1ncutNjx,yHx,y−Hxj,yjaj+∑k=1ntip1Nkx,y∑l=14Bl1x,y−Bl1xk,ykbkl+∑k=1ntip2Nkx,y∑l=14Bl2x,y−Bl2xk,ykbkl,where ui are the standard degrees of freedom (dofs); aj, bk1, and bk2 are additional dofs; Hx,y are the enrichment values of Heaviside-enriched nodes; Bl1 are the enrichment values of nodes around crack tip 1; Bl2 are the enrichment values of nodes around crack tip 2; Ni are the standard shape functions; ns is the number of standard nodes; ncut is the number of Heaviside-enriched nodes; and ntip1 and ntip2 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 ana-long single-edge crack (Figure 4(a)). It is horizontal and runs from the left vertical boundary of a rectangular plate of width (w) 100 m and height (h) 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 50×50, 100×100, and 200×200 elements. The elastic modulus (E) 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 σvm for the case with a 7.5 m long crack and 100×100 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 KI is [44]
(17)KI=Cσaπ,where σ is the dominant vertical tensile stress (equal to the vertical tractions applied in the example), a is the crack length, and C is a correction coefficient, equal to
(18)C=1.12−0.231aw+10.55aw2−21.71aw3+30.382aW4.Figure 4
Single edge crack: (a) schematic of aw×L plate with an a-long horizontal edge crack and loaded by vertical tractions σ; (b) shaded plot of the von Mises stress σvm for a 100×100 m2 plate mesh and a=7.5 m.
(a)(b)Table 1
Comparison between the SIF values (unit: 107 Pa/m2) of the numerical solution with those of the analytical solution for a plate with an edge crack.
Length (m)SIF numerical solutionSIF analytical solutionMesh50×50Mesh100×100Mesh200×2003.51.35691.35251.37011.36305.51.72961.71851.70961.70867.52.06992.00542.00101.99509.52.39762.37122.35952.245011.52.72472.69332.68242.4700


### 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 100×100-element mesh is shown. The corresponding analytical solution for the Mode-I SIF KI is [45]
(19)KI=σaπ,with σ here the dominant isotropic stress (equal to the normal tractions applied in the example).Figure 5
Central crack: (a) schematic of aw×L plate with an a-long horizontal central crack and loaded by normal tractions σ on all the edges; (b) shaded plot of the von Mises stresses σvm for a 100×100 m2 plate mesh and a=15 m.
(a)(b)In Table2, the comparison between the analytical solution and the numerical solution for different meshes is reported.Table 2
Comparison between the SIF values (unit: 107 Pa/m2) of the numerical solution with those of the analytical solution for a plate with a central crack.
Length (m)SIF numerical solutionSIF analytical solutionMesh50×50Mesh100×100Mesh200×2003.01.07890.75980.77260.79657.01.26191.18651.19711.216611.01.57601.50131.50371.525115.01.84041.77781.78321.781019.02.07522.03712.01892.0049Observing the results in Tables1 and 2, one may notice that when the mesh size is 100×100, 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.


## 3.1. Edge Crack

The first example refers to ana-long single-edge crack (Figure 4(a)). It is horizontal and runs from the left vertical boundary of a rectangular plate of width (w) 100 m and height (h) 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 50×50, 100×100, and 200×200 elements. The elastic modulus (E) 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 σvm for the case with a 7.5 m long crack and 100×100 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 KI is [44]
(17)KI=Cσaπ,where σ is the dominant vertical tensile stress (equal to the vertical tractions applied in the example), a is the crack length, and C is a correction coefficient, equal to
(18)C=1.12−0.231aw+10.55aw2−21.71aw3+30.382aW4.Figure 4
Single edge crack: (a) schematic of aw×L plate with an a-long horizontal edge crack and loaded by vertical tractions σ; (b) shaded plot of the von Mises stress σvm for a 100×100 m2 plate mesh and a=7.5 m.
(a)(b)Table 1
Comparison between the SIF values (unit: 107 Pa/m2) of the numerical solution with those of the analytical solution for a plate with an edge crack.
Length (m)SIF numerical solutionSIF analytical solutionMesh50×50Mesh100×100Mesh200×2003.51.35691.35251.37011.36305.51.72961.71851.70961.70867.52.06992.00542.00101.99509.52.39762.37122.35952.245011.52.72472.69332.68242.4700


## 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 100×100-element mesh is shown. The corresponding analytical solution for the Mode-I SIF KI is [45]
(19)KI=σaπ,with σ here the dominant isotropic stress (equal to the normal tractions applied in the example).Figure 5
Central crack: (a) schematic of aw×L plate with an a-long horizontal central crack and loaded by normal tractions σ on all the edges; (b) shaded plot of the von Mises stresses σvm for a 100×100 m2 plate mesh and a=15 m.
(a)(b)In Table2, the comparison between the analytical solution and the numerical solution for different meshes is reported.Table 2
Comparison between the SIF values (unit: 107 Pa/m2) of the numerical solution with those of the analytical solution for a plate with a central crack.
Length (m)SIF numerical solutionSIF analytical solutionMesh50×50Mesh100×100Mesh200×2003.01.07890.75980.77260.79657.01.26191.18651.19711.216611.01.57601.50131.50371.525115.01.84041.77781.78321.781019.02.07522.03712.01892.0049Observing the results in Tables1 and 2, one may notice that when the mesh size is 100×100, 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 Table3). The nodes of a crack k are reported in a column set from 2+10k−1 to 11+10k−1; column 1 contains the numbers of the standard nodes; a is the number of total standard nodes; m represents the number of the last enriched node in the last crack; m+1 and m+2 are the numbers of Heaviside-enriched nodes; and m+3, m+4, m+5, and m+6 are the numbers of the tip-enriched nodes.Table 3
Numbering scheme for the node-numbering matrix for thekth crack.
Columns1…2+10k−13+10k−14+10k−15+10k−1…10+10k−111+10k−1…1…2…3…m+1e.v.m+3e.v.…m+6e.v.…4…m+2e.v.…a−1…a…e.v.: enriched value.The governing equation for the whole domain can be written as(20)Ku=P,where P is the global load matrix, u is the global displacement matrix of Equation (16), and K is the global stiffness matrix, consisting of the standard stiffness matrix and the enriched stiffness matrices of the cracks; in formula,
(21)K=K0K0,1K0,2⋯K0,n−2K0,n−1K0,nK1,0K10⋯000K2,00K2⋯000⋯⋯⋯⋯⋯⋯⋯Kn−2,000⋯Kn−200Kn−1,000⋯0Kn−10Kn,000⋯00Kn,in which K0 is the standard stiffness matrix; K1, K2,…, Kn are the enriched stiffness matrices of cracks 1 to n, with n being the number of cracks; and K0,n is the cross-term between the standard elements and the enriched elements. The matrices are as follows:
(22)K0=∫ΩBiTDBjdΩ,i,j=u,Kn=∫ΩBiTDBjdΩ,i,j=a,b,K0,n=∫ΩBiTDBjdΩ,i,j=u,a,b,where D is the elasticity matrix, Bu, Ba, and Bb are strain differential operator matrices, respectively defined as:
(23)Bu=∂Ni∂x00∂Ni∂y∂Ni∂y∂Ni∂x,Ba=∂NiH−Hxi,yi∂x00∂NiH−Hxi,yi∂y∂NiH−Hxi,yi∂y∂NiH−Hxi,yi∂x,Bb=∂NiB1−B1xi,yi∂x0∂NiB2−B2xi,yi∂x00∂NiB1−B1xi,yi∂y0∂NiB2−B2xi,yi∂y∂NiB1−B1xi,yi∂y∂NiB1−B1xi,yi∂x∂NiB2−B2xi,yi∂y∂NiB2−B2xi,yi∂x∂NiB3−B3xi,yi∂x0∂NiB4−B4xi,yi∂x00∂NiB3−B3xi,yi∂y0∂NiB4−B4xi,yi∂y∂NiB3−B3xi,yi∂y∂NiB3−B3xi,yi∂x∂NiB4−B4xi,yi∂y∂NiB4−B4xi,yi∂x,with B1, B2, B3, and B4 the tip branch functions around the crack tips and B1xi,yi, B2xi,yi, B3xi,yi, and B4xi,yi 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 [48–51]. IIM is often used to calculate SIF values since it has a high degree of accuracy compared with other methods [52–55], although the involved process is quite complex; for the application, it is required to solve an energy integral, based on the J 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 x1,x2 of the crack, with origin at the tip and x1 aligned with it, the interaction integral I can be written as [2]
(24)I=∫ΓW1,2δ1j−σij1∂ui2∂x1−σij2∂ui1∂x1njdΓ,where δ is the Dirac delta function, n is the vector normal to Γ, and W is the interaction strain energy, equal to [1, 3]
(25)W1,2=σij1εij2=σij2εij1.The integralI can also be expressed by using the values of SIF KI and KII, respectively, for Mode-I and Mode-II, for both the present state and the auxiliary state, as follows:
(26)I=2EKI1KI2+KII1KII2,where E is Young’s modulus, and
(27)KI=E2I1,Mode‐I,(28)KII=E2I1,Mode‐II.Equation (24) can be transformed as follows:
(29)I=∫Ωσij1∂ui2∂x1+σij2∂ui1∂x1−W1,2δij∂q∂xjdΩ,where qx1,x2 is a weighting function. The solution of the integral in Equation (29) is relatively straightforward.With reference to Figure6(a), the radius R of the zone is assumed equal to the square root of three times the element area [56]. The weighting function q is as follows:
(30)qr=0,r≥R,1r<R,where r originates from the crack tip. When a portion of the zone is outside the domain, as in Figure 6(b), the weighting function (qx) on the edge of the domain is assumed equal to zero [57].Figure 6
Zone of integration around a tip for IMM: (a) zone inside the domain; (b) zone outside the domain.
(a)(b)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,h and l, 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.Figure 7
Schematic of the two-crack example: (a) case 1; (b) case 2.
(a)(b)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 (Figure7(b)) and Equation (29) for crack h becomes
(31)Ih=∫Ωhσh,ij1∂uh,i2∂xh,1+σh,ij2∂uh,i1∂xh,1−Wh1,2δh,1j∂qh∂xh,jdΩh,and, for crack l,
(32)Il=∫Qσl,ij1∂ul,i2∂xl,1+σl,ij2∂ul,i1∂xl,1−Wl1,2δl,1j∂ql∂xl,jdQl.The overlapped integral domainIlh is then
(33)Ilh=∫Ωlhσlh,ij1∂ulh,i2∂xlh,1+σlh,ij2∂ulh,i1∂xlh,1−Wlh1,2δlh,1j∂qlh∂xlh,jdΩlh,where Ωlh is the intersection zone. For Equations (31), (32), and (33), subscripts h and l indicate crack h and crack l, respectively, based on Equation (24). The remaining equations follow the same notation rule, so the comprehensive interaction integral for crack h is equal to the superposition of Ih and Ilh:
(34)Ihh=Ih+Ilh,and the comprehensive interaction integral for crack l is
(35)Ill=Il+Ilh.Thus, by employing Equations (27) and (28), SIFs in Mode-I and Mode-II for the two cracks are
(36)KI,l=E2Ill,Mode‐I,KII,l=E2Ill,Mode‐II,KI,h=E2Ihh,Mode‐I,KII,h=E2Ihh,Mode‐II.The results of the XFEM calculations in terms of the von Mises stress for two cases are shown in Figures8 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].Figure 8
Results of the two-crack example, case 1: (a) mesh; (b) shaded plot of the von Mises stressσvm (lengths in m).
(a)(b)Figure 9
Results of the two-crack example, case 2: (a) mesh; (b) shaded plot of the von Mises stressσvm (lengths in m).
(a)(b)In Figures10(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, h) and a 10 m long slanted crack (inclined 45° with respect to the x-axis) (the “major” crack, l). In case 3, the major crack is rather far from the minor crack (50 m from the tip of h to the center of l). In case 4, the two cracks are closer. It is assumed that the propagation of crack h is affected by the other crack when this one lies within the domain of the interaction integral of h. The major crack l can be seen as a border “arresting” the effect on the stress field of the minor crack h, so the actual interaction integral domain for crack h is equivalent to the entire circular area minus the shaded area (Figure 10(b)). The interaction integral for crack h is as Equation (31).Figure 10
Schematic of the two-crack example: (a) case 3; (b) case 4.
(a)(b)The interaction integralIΔh of the shaded area Δh of case 4 is (domain ΩΔh)
(37)IΔh=∫ΩΔhσΔh,ij1∂uΔh,i2∂xΔh,1+σΔh,ij2∂uΔh,i1∂xΔh,1−WΔh1,2δΔh,1j∂qΔh∂xΔh,jdΩΔh.Thus, the final interaction integral of crackh is equal to
(38)Ih′=Ih−IΔh.The results for these two cases are shown in Figure11. 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].Figure 11
Shaded plots of the von Mises stressσvm for the two-crack example: (a) case 3; (b) case 4 (lengths in m).
(a)(b)


### 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 Figure12(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.Figure 12
Shaded plots of the von Mises stressσvm for the multiple-crack example: (a) 5 parallel edge cracks (case 1); (b) 5 randomly located cracks (case 2) (lengths in m).
(a)(b)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 (Figure12(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.Figure 13
Shaded plots of the von Mises stressσvm for the multiple-crack example, case 2, for an increased propagation (lengths in m).


## 5.1. Two-Crack Example

In this example, two cracks,h and l, 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.Figure 7
Schematic of the two-crack example: (a) case 1; (b) case 2.
(a)(b)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 (Figure7(b)) and Equation (29) for crack h becomes
(31)Ih=∫Ωhσh,ij1∂uh,i2∂xh,1+σh,ij2∂uh,i1∂xh,1−Wh1,2δh,1j∂qh∂xh,jdΩh,and, for crack l,
(32)Il=∫Qσl,ij1∂ul,i2∂xl,1+σl,ij2∂ul,i1∂xl,1−Wl1,2δl,1j∂ql∂xl,jdQl.The overlapped integral domainIlh is then
(33)Ilh=∫Ωlhσlh,ij1∂ulh,i2∂xlh,1+σlh,ij2∂ulh,i1∂xlh,1−Wlh1,2δlh,1j∂qlh∂xlh,jdΩlh,where Ωlh is the intersection zone. For Equations (31), (32), and (33), subscripts h and l indicate crack h and crack l, respectively, based on Equation (24). The remaining equations follow the same notation rule, so the comprehensive interaction integral for crack h is equal to the superposition of Ih and Ilh:
(34)Ihh=Ih+Ilh,and the comprehensive interaction integral for crack l is
(35)Ill=Il+Ilh.Thus, by employing Equations (27) and (28), SIFs in Mode-I and Mode-II for the two cracks are
(36)KI,l=E2Ill,Mode‐I,KII,l=E2Ill,Mode‐II,KI,h=E2Ihh,Mode‐I,KII,h=E2Ihh,Mode‐II.The results of the XFEM calculations in terms of the von Mises stress for two cases are shown in Figures8 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].Figure 8
Results of the two-crack example, case 1: (a) mesh; (b) shaded plot of the von Mises stressσvm (lengths in m).
(a)(b)Figure 9
Results of the two-crack example, case 2: (a) mesh; (b) shaded plot of the von Mises stressσvm (lengths in m).
(a)(b)In Figures10(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, h) and a 10 m long slanted crack (inclined 45° with respect to the x-axis) (the “major” crack, l). In case 3, the major crack is rather far from the minor crack (50 m from the tip of h to the center of l). In case 4, the two cracks are closer. It is assumed that the propagation of crack h is affected by the other crack when this one lies within the domain of the interaction integral of h. The major crack l can be seen as a border “arresting” the effect on the stress field of the minor crack h, so the actual interaction integral domain for crack h is equivalent to the entire circular area minus the shaded area (Figure 10(b)). The interaction integral for crack h is as Equation (31).Figure 10
Schematic of the two-crack example: (a) case 3; (b) case 4.
(a)(b)The interaction integralIΔh of the shaded area Δh of case 4 is (domain ΩΔh)
(37)IΔh=∫ΩΔhσΔh,ij1∂uΔh,i2∂xΔh,1+σΔh,ij2∂uΔh,i1∂xΔh,1−WΔh1,2δΔh,1j∂qΔh∂xΔh,jdΩΔh.Thus, the final interaction integral of crackh is equal to
(38)Ih′=Ih−IΔh.The results for these two cases are shown in Figure11. 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].Figure 11
Shaded plots of the von Mises stressσvm for the two-crack example: (a) case 3; (b) case 4 (lengths in m).
(a)(b)


## 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 Figure12(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.Figure 12
Shaded plots of the von Mises stressσvm for the multiple-crack example: (a) 5 parallel edge cracks (case 1); (b) 5 randomly located cracks (case 2) (lengths in m).
(a)(b)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 (Figure12(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.Figure 13
Shaded plots of the von Mises stressσvm for the multiple-crack example, case 2, for an increased propagation (lengths in m).


## 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.



---

*Source: 1010174-2022-12-22.xml*