Numerical Comparison Research on the Solution of Stress Intensity Factors of Multiple Crack Problems

Guo Zhao

College of Architecture Engineering and Planning, Jiujiang University, Jiujiang, China.

DOI: 10.4236/apm.2020.1012044                

Abstract

A newly developed approach without crack surface discretization for modeling 2D solids with large number of cracks in linear elastic fracture mechanics is proposed with the eigen crack opening displacement (COD) boundary integral equations in this paper. The eigen COD is defined as a crack in an infinite domain under fictitious traction acting on the crack surface. Respect to the computational accuracies and efficiencies, the multiple crack problems in finite and infinite plates are solved and compared numerically using three different kinds of boundary integral equations (BIEs): 1) the dual BIEs require crack surface discretization; 2) the BIEs with numerical Green’s functions (NGF) without crack surface discretization, but have to solve a complementary matrix; 3) the eigen crack opening displacement (COD) BIEs in the present paper. With the concept of eigen COD, the multiple crack problems can be solved by using a conventional displacement discontinuity boundary integral equation in an iterative fashion with a small size of system matrix as that in the NGF approach, but without troubles to determine the complementary matrix. Solution of the stress intensity factors of multiple crack problems is solved and compared in some numerical examples using the above three computational algorithms. Numerical results clearly demonstrate the numerical models of eigen COD BIEs have much higher efficiency, providing a newly numerical technique for multiple crack problems. Not only the accuracy and efficiency of computation can be guaranteed, but also the overall properties and local details can be obtained. In conclusion, the numerical models of eigen COD BIEs realize the simulations for multiple crack problems with large quantity of cracks.

Keywords

Multiple Crack Problems, Boundary Integral Equations, Eigen Crack Opening Displacements, Eshelby Matrix, Stress Intensity Factors

1. Introduction

Solid materials that contain large number of cracks can be the fracture mechanical models for many structural materials, especially for lots of brittle and quasi-brittle materials, such as cementitious materials in ref. [1], concretes structure in ref. [2], rocks in ref. [3], ceramics in refs. [4] [5], biological bones in ref. [6] and metallic materials in ref. [7] etc. The cracks may be formed by natural growth, also may be produced in the process of working. It is easy to be a fracture failure under external load for these solid materials. In other words, the existence of cracks especially multiple cracks is the main cause of fracture failure for these solid materials. Therefore, the prediction of multiple crack behaviors and the strong interaction between cracks need to be researched. The fracture failure always happens due to strong interactions among cracks or crack interconnections each other, where the strong interactions or interconnections depend on the applied external loading the geometrical configuration of cracks, such as crack sizes, crack shapes, crack locations as well as crack orientations. However, it is impossible to know the distribution information of all the cracks beforehand, particularly when the quantity of cracks up to hundreds even thousands. As a result, the analysis of strong interactions among cracks or crack interconnections can be a numerical tool for researchers and engineers in designing general structures.

In 1957, Irwin first postulated the concept of stress intensity factor (SIF). SIF became one of the main parameters which need to be obtained, depending on the stress field around the crack tips. Based on the value of SIFs, crack behaviors such as fatigue propagation can be investigated. The analysis of SIFs plays a greatly important role and real concern in linear elastic fracture mechanisms (LEFM), since solids with multiple cracks become one class of the most important problems in fracture mechanics. However, as mentioned earlier, the exact crack sizes, shapes, locations as well as orientations are unable to know in advance, it is an enormous challenge for researchers to seek an appropriate numerical model for LEFM problems, leading to huge complexities such as the computational efficiency and storage of solids with multiple cracks. As a result, most of them are incredibly difficult to solve by analytical procedures in ref. [8]. It is urgently needed to develop more reliable, accurate and efficient numerical processing techniques.

The finite element method (FEM), now widely available in many fields, not only has achieved a great reputation after decades of development in the applications of fracture mechanics, but also has many mature commercial software tools, such as ANSYS and ABAQUS,etc. Despite of its widely spread popularity, solids with cracks in large number are probably one class of the most difficult problems to simulate since the crack tips need to be discretized. To the author’s best knowledge, the advantage of the FEM is to model the solids as a whole, not enough to understand the local information.

Owing to the advantage of boundary discretization only, the boundary element method (BEM) is subsequently considered as an efficient numerical method to deal with multiple crack problems. In other words, the BEM reduces the dimension of original problems compared with the FEM, enjoyed absolute advantages in improving the speed and accuracy of calculation. Inevitably, the conventional BEM also has its own weaknesses in solving large-scale cracks due to its dense and asymmetric of the final system matrix. It is well known there are generically two difficulties for numerical modeling in ref. [9]. The first difficulty is how to accurately describe the stress field around the crack tips. The second difficulty is how to deal with the degeneration of the boundary integral equation resulting from the geometric coincidence of the upper crack surface and lower crack surface. To overcome the former difficulty, special discretization element as in the FEM can be used as a good refinement of the mesh in ref. [10], while sub-region technique can be usually applied to avoid the degeneration of the integral equations for the latter in refs. [11] [12].

A combination or dual boundary integral equations (BIEs) in refs. [10] [13] [14] [15] are described alternatives and compared in the present paper, where conventional displacement integral equation and traction integral equation are simultaneous utilized correspondingly to the upper and lower crack surface, respectively. Since all unknown information of the crack surfaces and outer boundaries are contained in the final system matrix, the size of dual BIEs will be a sharply rise when the quantity of cracks increase. The fast multipole expansion technology embedded in the conventional BEM, which is known as the fast multipole boundary element method (FMBEM), may be employed to solve multiple crack problems in refs. [16] [17].

The analytical solutions for most of multiple crack problems are restricted to obtain, except a few regular cracks. In order to eliminate the unknowns appearing on the crack surfaces, an analytical Green’s functions, known as the Erdogan function in ref. [18] was used as the fundamental solution in conventional BEM in ref. [19]. Unfortunately, most of the above Erdogan function was always in plural, only restricted to lower dimension or single crack with simple size and shape, not available for solid materials with large number of cracks.

In 1995, Telles et al. put forward the numerical Green’s function (NGF) procedure, which can be modeled as the superposition of a basic plus a complementary solution problem in refs. [20] [21], where the basic part is known as the Kelvin analytical solution problem. The complementary problem can be described as a single crack embedded into an infinite domain under far-field loading with the negative of the former Kelvin traction. However, the size of the final system matrix of the complementary solution problem of the NGF will grow large immediately with the number of crack increases, which have to be determined numerically step by step by using the conventional BEM. The processing technology seems very troublesome. Kachanov in refs. [22] [23] also proposed an algorithm to simulate the interaction between cracks. The original crack problem in an infinite domain under surface loading can be decomposed into a superposition model of uniform and non-uniform terms. It is assumed that the interaction between non-uniform terms can be ignored in ref. [24]. In fact, it is not advisable to ignore this term since the interaction between cracks will be much stranger when the distance between cracks is quite close.

To overcome these weaknesses of the above mentioned algorithms, the concept of the eigen crack opening displacement (COD) are firstly proposed in author's previous paper in ref. [25] and reintroduced in the present paper as shown in Section 2.1. Based on the concept of eigen COD, the basic format of eigen COD BIEs and its iterative fashion are established step by step with a small size of final system matrix in refs. [26] [27].

In order to efficiently and numerically simulate the interaction effects between the cracks, a superposition technique as in NGF is applied and the local Eshelby matrix derived from the traction BIE in discretized form is introduced subsequently. Due to limited space, more details will be discussed in Section 3.2. The eigen strain BIEs for multiple inclusion problems can be also found in refs. [28] [29]. In fact, a generalized inclusion can be degenerated into a crack in physical. Correspondingly, the eigen strain BIEs will degenerate into the eigen COD BIEs indirectly.

In the present paper, two dimensional multiple crack problems in finite/infinite plates are numerically considered and discretized by using the boundary point method (BPM) in refs. [30] [31]. The rest structure of this paper is organized as follows: in Section 2, the basic formulations of the above three kinds of BIEs are briefly given, respectively. In Section 3, the concept of eigen COD is reintroduced, and also the computational model is established. After that, the local Eshelby matrix is presented with the numerical treatments in details. In Section 4, solution procedures of the eigen COD algorithm are given with its iteration fashion and convergence criterion. In Section 5, several numerical examples in large number of cracks (the crack quantity up to hundreds) in finite/infinite plates are calculated and compared to show the accuracy and the efficiency with the proposed three kinds of boundary integral equations.

2. Three Kinds of Boundary Integral Equations

2.1. Eigen COD Boundary Integral Equations

Without loss of generality, one two dimensional elastic domain Ω with outer boundary Гwhich contains multiple N_C cracks is considered, where N_C is the total number of cracks. Suppose there is a source point y, the displacements at y can be expressed with the unknown displacement discontinuities as follows in refs. [22] [23]:

$\begin{array}{l}\gamma u_i\left(y\right)={\displaystyle \underset{\Gamma }{\int }{\tau }_j\left(x\right)u_{ij}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}-{\displaystyle \underset{\Gamma }{\int }{u}_j\left(x\right){\tau }_{ij}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}\\ -{\displaystyle \underset{m=1}{\overset{N_C}{\sum }}{\displaystyle \underset{A_m^{+}}{\int }\Delta u_j\left(x\right){\tau }_{ij}^{*}\left(x,y\right)\text{d}A}\left(x\right)},y\in \left(\Gamma \cup \Omega \right)\backslash A_m\left(m=1,\cdots ,N_C\right)\end{array}$ (1)

where $u_{ij}^{*}$ and ${\tau }_{ij}^{*}$ are respectively the displacement and traction fundamental solutions in ref. [9]. A_m represents the crack number m. x and y represent the field point and the source point respectively. γis the free term coefficient, which depends on the location of the source point and the geometry of the boundary. γ = 1 when the source point y is inside in the domain Ω, γ = 0.5 if it is on the outer boundary Гand only if the outer boundary is smooth and continuous. Δu_i are the displacement discontinuities, or the CODs at the crack surfaces, which can be defined as follows:

$\Delta u_i\left(x\right)=u_i\left(x\right)|{}_{x\in A^{+}}-u_i\left(x\right)|{}_{x\in A^{-}}$ (2)

where A⁺ represents the upper surface of one crack, A^- represents the lower surface, respectively. It is obvious from Equation (1) that if all the unknown displacement discontinuities (CODs) are obtained in advance, the multiple crack problems can be solved by Equation (1) in a discrete form as in conventional BEM in ref. [9]. The stress and traction at source point y can be derived by the stress BIE and traction BIE, respectively, with the unknown CODs as follows:

$\begin{array}{l}\gamma {\sigma }_{ij}\left(y\right)={\displaystyle \underset{\Gamma }{\int }{\tau }_k\left(x\right)u_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}-{\displaystyle \underset{\Gamma }{\int }{u}_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}\\ -{\displaystyle \underset{m=1}{\overset{N_C}{\sum }}{\displaystyle \underset{A_m^{+}}{\int }\Delta u_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}A}\left(x\right)},y\in \left(\Gamma \cup \Omega \right)\backslash \left(\cup A_m\right)\end{array}$ (3)

$\begin{array}{l}\gamma {\tau }_i\left(y\right)=n_j\left(y\right){\displaystyle \underset{\Gamma }{\int }{\tau }_k\left(x\right)u_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}-n_j\left(y\right){\displaystyle \underset{\Gamma }{\int }{u}_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}\\ -n_j\left(y\right){\displaystyle \underset{m=1}{\overset{N_C}{\sum }}{\displaystyle \underset{A_m^{+}}{\int }\Delta u_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}A}\left(x\right)},y\in \left(\Gamma \cup \Omega \right)\backslash \left(\cup A_m\right)\end{array}$ (4)

where n_j represents the unit outer direction cosine at the source point y. Equations (1), (3) as well as (4) are designated as the eigen COD BIEs in the present paper. The CODs in those equations are unknown and should be computed with the concept of eigen COD, which will be shown in Section 3.1. Anyway, as it mentioned before, because there are general two difficulties in conventional BEM, i.e., accurately describe the stress field around the crack tips and the degeneration of the boundary integral equations resulting from crack surfaces coincide, Equation (1) cannot be employed alone to solve multiple crack problems.

2.2. Dual Boundary Integral Equations

For the numerical solution of multiple crack problems with dual formulations, the boundary integral equations used are as follows for the crack surfaces, together with the use of Equation (1) for the outer boundary:

$\begin{array}{l}\frac{1}{2}\left\{u_i\left(y^{+}\right)+u_i\left(y^{-}\right)\right\}\\ ={\displaystyle \underset{\Gamma }{\int }{\tau }_j\left(x\right)u_{ij}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}-{\displaystyle \underset{\Gamma }{\int }{u}_j\left(x\right){\tau }_{ij}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}\\ -{\displaystyle \underset{m=1}{\overset{N_C}{\sum }}{\displaystyle \underset{A_m^{+}}{\int }\Delta u_j\left(x\right){\tau }_{ij}^{*}\left(x,y\right)\text{d}A}}\left(x\right),y\in A_m\left(m=1,\cdots ,N_C\right)\end{array}$ (5)

$\begin{array}{l}{\tau }_i\left(y\right)=n_j\left(y\right){\displaystyle \underset{\Gamma }{\int }{\tau }_k\left(x\right)u_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}-n_j\left(y\right){\displaystyle \underset{\Gamma }{\int }{u}_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}\\ -n_j\left(y\right){\displaystyle \underset{m=1}{\overset{N_C}{\sum }}{\displaystyle \underset{A_m^{+}}{\int }\Delta u_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}A}}\left(x\right),y\in A_m\left(m=1,\cdots ,N_C\right)\end{array}$ (6)

which means that both the displacement and the traction equations should be employed when the source point y is placed on the crack surface, thus avoiding the degeneration of the system matrix and getting rid of interfaces between the sub-regions into the modeling scheme in refs. [11] [12]. The limitation of the use of the dual BIEs for multiple crack problems is that the size of the finial system matrix will grow large shapely when the number of crack increases, since all the unknowns are placed both on the crack surfaces and outer boundaries of the domain.

2.3. Numerical Green’s Function

Since the analytical solutions for most of multiple crack problems are restricted to obtain, Telles et al. put forward the numerical Green’s function in 1995, which can be modeled as the superposition of a basic plus a complementary solution problem in refs. [20] [21], which became an important alternative in solving fracture mechanics problem. The main model of NGF is represented as follows:

$\gamma u_i\left(y\right)={\displaystyle \underset{\Gamma }{\int }{\tau }_j\left(x\right)u_{ij}^G\left(x,y\right)\text{d}\Gamma \left(x\right)}-{\displaystyle \underset{\Gamma }{\int }{u}_j\left(x\right){\tau }_{ij}^G\left(x,y\right)\text{d}\Gamma \left(x\right)}$ (7)

where the superscript G means the Green's function fundamental solution as follows:

$u_{ij}^G\left(x,y\right)=u_{ij}^{*}\left(x,y\right)+u_{ij}^C\left(x,y\right)$ (8)

${\tau }_{ij}^G\left(x,y\right)={\tau }_{ij}^{*}\left(x,y\right)+{\tau }_{ij}^C\left(x,y\right)$ (9)

The complementary parts of the NGF, $u_{ij}^C$ and ${\tau }_{ij}^C$ need to be determined by numerical means by using the traction BIE in refs. [20] [21]. It can be seen from Equation (7) that the discretization of crack surfaces has been elegantly avoided. It is a merit of the usage of the NGF.

3. The Computational Model of Eigen COD BIEs

3.1. Concept of Eigen COD

In Section 2.1, the eigen COD BIEs have been defined. The CODs are the unknowns and need to be computed step by step with the concept of eigen COD, which will be discussed in what follows.

For convenience, we consider an infinite domain containing a traction free crack A under far-field loading σ as shown in Figure 1(a). The crack has a COD response. This cracked model can be decomposed equivalently into two parts: the first part is that an enclosed crack under far-field loading, and both of the tractions with a negative minus sign along the upper and lower surface of the crack as shown in Figure 1(b). It is important that there is no COD response in

Figure 1. Fictitious tractions on crack surface.

this part. The second part is that an opened crack under the fictitious tractions which applied on the upper and lower surface of the crack but without far-field loading as shown in Figure 1(c). Unlike the first part, the opened crack in Figure 1(c) has same COD response as that of Figure 1(a). After this treatment, the fictitious tractions applied on the crack surfaces can be numerically calculated with no difficulty regardless the existence of crack by using Equation (4), i.e., tractions boundary integral equation with the eigen COD unknowns. The concept of eigen COD, Δu_i, redefined as the crack opening displacements of current crack A under the fictitious tractions, τ_i, acting on the surface.

To better explain the following processes, an infinite domain contains only one crack is considered. By making the source point y on the crack surface, the expression of the traction, a hypersingular traction integral equation, can be derived from Equation (4) in the global coordinate after a limiting process as follows:

$n_j\left(y\right)\text{HFP}{\displaystyle \underset{A^{+}}{\int }\Delta u_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}A\left(x\right)}=-{\tau }_i\left(y\right),y\in A^{+}$ (10)

where the capital HFP in the left side means the Hadamard finite part when the source y and the field point x coincide, which has hypersingularity. Without loss of generality, the single crack can be modeled as the case of a straight crack in refs. [20] [21]. It is noticed that there are always two respectively loading modes, i.e., the opening mode and the sliding mode, having similarly mathematical expression in fact. By the way, there is an explicitly analytical solution, ${\tau }_{ijk}^{*}$ , for a single crack in infinite domain in ref. [12]. Suppose the straight crack with length of 2a. Embedding the explicitly analytical solution for a single crack into Equation (10), a simplified form can be reduced as follows:

$\frac{\mu }{2\pi a\left(1-\nu \right)}\text{HFP}{\displaystyle \underset{-1}{\overset{+1}{\int }}\frac{{\delta }_i}{r^2\left(x,y\right)}\text{d}A\left(x\right)}=-{\tau }_i\left(y\right)$ (11)

where μ and ν are the shear modulus and Poisson’s ratio of the solid material, respectively. r represents the distance between field point x and source point y. δ_i stands for the COD for convenience.a in denominator means the half length of the crack. It can be seen from Equation (11) that the COD, δ_i, has to be solved numerically.

To eliminate the subscript i, the explicitly analytical expression of the HFP integral in ref. [30] can be employed into Equation (11) and transform the global coordinate to the local coordinate, the discrete form of Equation (11) can be obtained as follows:

$\begin{array}{l}\frac{\mu }{2\pi a\left(1-\nu \right)}\{{\displaystyle \underset{j=1,j\ne k}{\overset{N_G}{\sum }}\frac{w_j}{{\left({\xi }_j-{\xi }_k\right)}^2}{\delta }_j}-\left[{\displaystyle \underset{j=1,j\ne k}{\overset{N_G}{\sum }}\frac{w_j}{{\left({\xi }_j-{\xi }_k\right)}^2}}+\frac{2}{1-{\xi }_i^2}\right]{\delta }_k\\ \underset{}{\overset{}{\begin{array}{c}\\ \end{array}}}-\left[{\displaystyle \underset{j=1,j\ne k}{\overset{N_G}{\sum }}\frac{w_j}{{\xi }_j-{\xi }_k}}-\ln \left(\frac{1-{\xi }_k}{1+{\xi }_i}\right)\right]\frac{\partial {\delta }_k}{\partial \xi }+\frac{w_k}{2}\frac{{\partial }^2{\delta }_k}{\partial {\xi }^2}\}=-{\tau }_k,k=1,\cdots ,N_G\end{array}$ (12)

where the subscript j and k are respectively the collocation points and the Gauss stations. w_j and ξ_j are respectively the Gauss weight functions and station functions. N_G is the total discrete Gauss number. It is noticed that in Equation (12) there exists the derivatives of the COD, $\partial {\delta }_k$, which can be computed by using a Lagrange polynomial interpolation as follows:

$\delta ={\displaystyle \underset{k=1}{\overset{N_G}{\sum }}{l}_k{\delta }^k}$ (13)

where l_k represent the coefficients of the Lagrange polynomials interpolation with an order number of (N_G + 2), making the Equation (13) satisfy that $\delta \left(+a\right)=\delta \left(-a\right)=0$.

Finally, by embedding Equation (13) into Equation (12), an equation in matrix form can be written for both opening mode and sliding mode of the CODs as follows:

$a^{-1}{S}^0\delta =\tau$, $S^0=\left[\begin{array}{cc}S& 0\\ 0& S\end{array}\right]$ (14)

where the vector form $\delta$ and fictitious tractions $\tau$ have a size of (2N_G × 1). $S$ is discretized from Equation (12) in matrix form with a size of (N_G×N_G), which can be defined as the basic matrix for multiple crack problems in the present paper.

3.2. The Local Eshelby Matrix

In this section, an infinite domain contains multiple cracks is considered as shown in Figure 2. To efficiently and numerically simulate the interaction effects among cracks, the superposition technique as in the NGF is applied and the local Eshelby matrix derived from the traction BIE in discrete form is introduced in what follows.

The main idea is to first select a crack as the current crack A (research object), and then divide all cracks into two groups according to a distance of other cracks to the current crack, defined as the adjacent group and the far-field group. It needs to be explained that the intermediate crack is usually selected as the current crack

Figure 2. The group definition for multiple crack problems.

to facilitate the calculation in the present paper. The adjacent group (interior of the dashed circle) is consisted of the cracks with relatively small distances while the far-field group (exterior of the dashed circle) is made up of those with relatively large distances. Correspondingly, the adjacent group has strong interaction effects while the far-field group has week interaction effects to the current crack, respectively.

In order to better explain the derivation of local Eshelby matrix, only the adjacent group is considered, that is, the dashed circle around the current crack A contains all the cracks in infinite domain while the far-field group just contains crack number of zero. The number of cracks in the adjacent group being denoted as N_L (N_L = N_C). By placing the source point y on the crack surfaces, Equation (10) can be rewritten as follows:

$\begin{array}{l}{n}_j\left(y\right)\text{HFP}{\displaystyle \underset{A_l^{+}}{\int }\Delta u_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}A\left(x\right)}\\ +n_j\left(y\right){\displaystyle \underset{m=1,m\ne l}{\overset{N_L}{\sum }}{\displaystyle \underset{A_m^{+}}{\int }\Delta u_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}A\left(x\right)}}=-{\tau }_i\left(y\right),y\in A_l\end{array}$ (15)

The first integral in the left side of Equation (15) has hypersingularity with the same structure as that of Equation (11), which can be discretized and computed by using Equation (12) and Equation (14), respectively. The second integral in the left side of Equation (15) is regular, which is easy to compute without any difficulty. After discretization and arrangement, Equation (15) can be rewritten in matrix form as follows:

$\left[\begin{array}{cccc}{a}_1^{-1}{S}^0& S_{1,2}& \cdots & S_{1,N_L}\\ S_{2,1}& a_2^{-1}{S}^0& \cdots & S_{2,N_L}\\ ⋮& ⋮& \ddots & ⋮\\ S_{N_L{}_{,1}}& S_{N_L{}_{,2}}& \cdots & a_{N_L}^{-1}{S}^0\end{array}\right]\left\{\begin{array}{c}{\delta }_1\\ {\delta }_2\\ ⋮\\ {\delta }_{N_L}\end{array}\right\}=\left\{\begin{array}{c}{\tau }_1\\ {\tau }_2\\ ⋮\\ {\tau }_{N_L}\end{array}\right\}$ (16)

where a_k in the left side represent the half length of the kth crack. $S_{km}$ represent the sub-matrices which derived from the discrete form corresponding to the second regular integral in the left side of Equation (15).

By inversing the total square matrix to the right side, the eigen COD, δ_k, of the kth crack can be obtained as follows:

${\delta }_{\left(k\right)}=S_k{\tau }_{\left(k\right)},k=\text{1},\text{2},\cdots ,N_C$ (17)

where the right vector ${\tau }_{\left(k\right)}={\left\{{\tau }_1,{\tau }_2,\cdots ,{\tau }_{N_L}\right\}}_k^{\text{T}}$ are the fictitious tractions of all cracks in the adjacent group with a size of ((2N_G × N_L) × 1). The left vector ${\delta }_{\left(k\right)}$ are the eigen COD of the kth crack with a size of (2N_G × 1). $S_k$ are the inverse matrix of Equation (16), which is named as the local Eshelby matrix with a size of (2N_G × (2N_G × N_L)) in the present paper. Since the discrete Gauss numberN_G and the adjacent number N_L are always a limit number, the matrix $S_k$ have a small size.

It should be noted that the above local Eshelby matrix $S_k$ are distinctly different for each crack. Once the radius of the dashed circle is defined in Figure 2, the adjacent group will be given for the current crack, then the local Eshelby matrix $S_k$ which reflect the interaction among cracks will be correspondingly given. In this way, there are many local Eshelby matrices with small sizes for different current crack. However, from the perspective of computational efficiency, the computational efficiency of solving a cracked problem with many small matrices must be much higher than that of a single but huge matrix as in dual BIEs and NGF.

3.3. Stress Intensity Factors

In order to avoid directly compute stresses or displacements around the crack tips due to its singularity, the SIFs are computed by using a polynomial expansion of fictitious tractions in this section.

Under normal circumstance, we can define the normal direction (opening mode) as n and the tangential direction (sliding mode) as t, respectively, with respect to the crack in the local coordinate as shown in Figure 1(c). The fictitious tractions, τ, in either normal or tangential direction, can be fitted approximately by using a polynomial expansion of an order N_P as follows:

$\tau =c_0+c_1\frac{\eta }{a}+c_2{\left(\frac{\eta }{a}\right)}^2+\cdot \cdot \cdot +c_{N_P}{\left(\frac{\eta }{a}\right)}^{N_P},\eta \in \left[-a,+a\right]$ (18)

where c_i represent the coefficients of the polynomial.

According to the analytical solution of fracture mechanics [32] [33], the formula of stress intensity factors for an infinite domain contains one crack under loading on the crack surfaces are as follows ( $\eta \in \left[-a,+a\right]$ ):

$K_{\mathrm{I}}^R=\frac{1}{\sqrt{\pi a}}{\displaystyle {\int }_{-a}^a\frac{p\left(\eta \right)\left(a-\eta \right)}{\sqrt{a^2-{\eta }^2}}\text{d}\eta }$ (19)

$K_{\mathrm{I}}^L=\frac{1}{\sqrt{\pi a}}{\displaystyle {\int }_{-a}^a\frac{p\left(\eta \right)\left(a+\eta \right)}{\sqrt{a^2-{\eta }^2}}\text{d}\eta }$ (20)

where p is the loading on the crack surfaces, and p is the same with the fictitious tractions τ in the present paper.

Then embedded Equation (18) into Equation (19) and (20), respectively, the SIFs at the two crack tips of a crack can be easily obtained as follows:

$\begin{array}{l}\left\{\begin{array}{c}{K}^R\\ K^L\end{array}\right\}=\{c_0+\frac{1}{2}\left(\pm c_1+c_2\right)+\frac{1\cdot 3}{2\cdot 4}\left(\pm c_3+c_4\right)\\ +\cdots +\frac{1\cdot 3\cdot \cdot \cdot \left(N_P-1\right)}{2\cdot 4\cdot \cdot \cdot N_P}\left(\pm c_{N_P-1}+c_{N_P}\right)\}\sqrt{\pi a}\end{array}$ (21)

where K^R and K^L represent the SIFs of right and left tips of the crack, respectively. It is important to note that the troubles of computing stresses or displacements around the crack tips can be avoided by using Equation (21), which only related to the coefficients of the polynomial expansion of the fictitious tractions.

4. Solution Procedures of the Eigen COD BIEs

It is obvious that all the unknown CODs and the fictitious tractions are related to the outer loading mode, the geometries, the quantities as well as the distributions of cracks in finite/infinite domain, which having mutual interactions among them. The solution procedures of the eigen COD BIEs are mainly divided into four stages as follows, i.e., the initiation stage, the iteration stage, the convergence check stage and the post process stage, respectively.

4.1. The Initiation Stage

1) Initializing all the data related to the finite/infinite domain, such as boundary conditions and geometric information of cracks, etc.;

2) Computing the local Eshelby matrices, $S_k$ ;

3) Computing the boundary unknowns by using the displacement eigen COD boundary integral Equation (1);

4) Computing the fictitious tractions of all cracks by using the traction eigen COD boundary integral Equation (4);

5) Computing the initial SIFs by using Equation (21).

4.2. The Iteration Stage

Since all cracks are divided into two groups, adjacent group and far-field group, the eigen CODs of all cracks are correspondingly calculated with respect to the fictitious tractions by two parts as follows:

1) The first part is calculated by using the fictitious tractions of the current crack and other cracks in the adjacent group via the local Eshelby matrix Equation (17);

2) The second part is calculated by using the fictitious tractions of all the cracks in the far-field group. A modified traction boundary integral equation of Equation (4) can be obtained as follows:

$\begin{array}{l}{\tau }_i\left(y\in A_l^{+}\right)=n_j\left(y\right){\displaystyle \underset{\Gamma }{\int }{\tau }_k\left(x\right)u_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}-n_j\left(y\right){\displaystyle \underset{\Gamma }{\int }{u}_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}\Gamma \left(x\right)}\\ -n_j\left(y\right){\displaystyle \underset{m=1,m\notin N_L^l}{\overset{N_C}{\sum }}{\displaystyle \underset{A_m^{+}}{\int }\Delta u_k\left(x\right){\tau }_{ijk}^{*}\left(x,y\right)\text{d}A}}\left(x\right),l=1,\cdots ,N_C\end{array}$ (22)

where $N_L^l$ represent the cracks of adjacent group centered at the current crack l. Then computing the corresponding SIFs of all cracks using Equation (21).

4.3. The Convergence Check Stage

In this block, we defined the maximum iteration error as follows:

$K_{\max }=\max \left|K^{\left(k\right)}-K^{\left(k-1\right)}\right|$ (23)

where K_max represents the maximum difference of the SIFs between the two iterations. The superscript k is the iteration count.

The convergence criterion of the present eigen COD BIEs is chosen as follows:

$K_{\max }/\sigma \sqrt{\pi a}\le {10}^{-3}$ (24)

1) If the convergence criterion is not satisfied Equation (24), it should be solved the boundary unknowns again and then return to the previous iteration stage;

2) If the convergence criterion is satisfied Equation (24), then go to the next post process stage in what follows.

4.4. The Post Process Stage

The post process stage can be carried out according to the computational needs or the interests of the researchers such as the overall properties and the fracture properties of these solid materials, as well as the local stresses or strain field around the crack tips, etc.

In conclusion, the solution procedures of the present eigen COD BIEs can be described in the flow chart as shown in Figure 3.

5. Numerical Examples

In this section, solutions of stress intensity factors of multiple crack problems in finite/infinite plates are presented to demonstrate the accuracy and the efficiency

Figure 3. The flow chart of the solution procedures of eigen COD.

by using the above three kinds of boundary integral equations. The numerical examples are solved by using a desk-top computer Dell with Intel Core Dual CPU, 2.50 GHz, 8 Gb of memory.

5.1. Stress Intensity Factors in Finite Square Plates

5.1.1. Square Plate with Four Cracks

The first example is a square plate (W = H) in tension contains four cracks (N_C = 4) with the same length (a = 0.2 W) as shown in Figure 4, which was previously analyzed by Chen & Chang in ref. [34] by using FEM. The outer boundaries of the square plate and the cracks are respectively discretized by using 200 and 9 nodes in dual BIE and eigen COD algorithm for this example. The number of cracks in adjacent group is set as N_L = 4, the same as that of the total cracks (N_L = N_C) in the eigen COD algorithm.

The stress intensity factors are computed with the variation of angle θ of the left and the right cracks symmetrically while the upper and the bottom cracks are kept stationary. The computed results are presented and compared in Figure 5 for

Figure 4. A square plate with four cracks of equal size.

Figure 5. Normalized stress intensity factors as a function of tilting angle θ.

the stress intensity factors at the crack tips A, B and C. The results by Chen & Chang in ref. [34] at θ = 0˚, 30˚ and 60˚ are also presented in Figure 5 for comparison. It is seen from Figure 5 that the good agreement exists among the dual BIE and the eigen COD algorithms and those by Chen & Chang in ref. [34], showing the effectiveness of the eigen COD BIE approach.

5.1.2. Square Plate with Multiple Cracks

The second example is a square plate in tension with multiple cracks, and two typical examples are as illustrated in Figure 6. The total number of cracks are chosen as $N_C={\left(2n+1\right)}^2,n=1,2,\cdots$, where n is an integer. In this block, n is chosen as 3 and 5 as shown in Figure 6(a) and Figure 6(b), respectively. The size of each crack is set as 2a = 0.4 W/n, where W = H are the width and the height of the plate, respectively. The stress intensity factors of the central crack are computed and compared by using the three algorithms. The number of cracks in adjacent group is set as N_L = 9 in the eigen COD algorithm.

The effects of the number adopted for crack discretization, N_G, are compared in Table 1 with the total crack number N_C = 9 (the tilting angle θ = 55.86˚) and compared in Table 2 with the total crack number N_C = 25 (the tilting angle θ = -55.75˚) for the three algorithms. It is seen from the Table 1 and Table 2 that the computed results are not sensitive to the number of N_G, which become stable when the number for crack discretization is equal to and greater than 13 (N_G ≥ 13) in all these examples.

The CPU time are compared as listed in Table 3 of the three algorithms. It can be obvious seen that the calculation efficiency of the proposed eigen COD approach is at least 20 times and 40 times higher compared with the NGF and the Dual BIEs, respectively, when the total crack number is set N_C = 121 (the tilting angle θ = 60.50˚). In Dual BIEs, the final system matrix contains all unknown information both on boundaries and crack surfaces, the size will be grow large with the total crack number increases, while in the NGF algorithm, the size of final system matrix to determine the complementary solution also grow large with the total crack number increases.

5.2. Stress Intensity Factors in Infinite Plates

In this section, infinite plates contain multiple rows of periodical cracks are

Figure 6. The square plates in tension with multiple cracks of equal size. (a)N_C = 9; (b)N_C = 25.

Table 1. The normalized SIF with different N_G (N_C = 9, θ = 55.86˚).

Table 2. The normalized SIF with different N_G (N_C = 25, θ = −55.75˚).

Table 3. The CPU time (s) of the three algorithms.

considered. In order to make the problem being much convenient, the crack problems are composed of cracks with the same size, orientation and spacing. The cracks in arrays considered are placed periodically in a number of configurations. Since all the cracks in arrays are under the same loading condition, one crack in the array is designated as the representative crack while the influences of all the other cracks, serving as an infinite series, on the representative crack can be summed up.

5.2.1. Periodical Collinear Cracks in Horizontal Line

The first example, as shown in Figure 7, is a row of periodic collinear cracks in horizontal line under far-field uniform tension perpendicular to the crack surfaces. In this computational model, up to 501 cracks (the total number of crack N_C = 501) are taken into consideration instead of using an infinite number. It should be noted that the total number of crack is set to be an odd number so that the current crack can be centrally arranged for convenience.

The calculated results for I-mode SIFs are expressed as shown in Table 4. It is important to point out that the unknowns are corrected step by step in an iterative fashion. Since our concern is infinite plate problems, with the crack number N_C increasing, the absolute error between computed results and theoretical solutions is decreasing. It can be seen that the computing results are matched well by using the three algorithms, even though at a/b = 0.9, the error is still less than 1%. It is so trivial that we could assume N_C = 501 is enough to describe infinite array of periodically spaced collinear cracks problem. If the local number of cracks N_L and the gauss points N_G increase, the accuracy will be even more satisfactory.

5.2.2. Double Periodical Collinear Cracks

The last example in this block is an infinite plate containing double periodical collinear cracks with the same length under far-field uniform tension perpendicular to the crack surfaces as shown in Figure 8. The number of row and column are equally set to be an odd number N so that the geometrical center of the

Figure 7. Periodical collinear cracks with equal size.

Table 4. The normalized I-SIF comparsions (N_L = 9,N_G = 13).

Figure 8. Double periodical collinear cracks with equal length.

problem is located at the center. The total number of crack is (N × N).

The normalized stress intensity factors of the central crack at different number of rows N are calculated and compared as listed in Table 5 with different number of row N. It can be seen that the numerical results agree well with ref. [35], showing the effectiveness of the proposed approaches.

The CPU time are also compared as listed in Table 6, showing high efficiency of the eigen COD algorithm. This is because, as mentioned previously, the final system matrix of the dual BIE algorithm contains all unknowns both on the boundaries and crack surfaces, where the size will be always largely increase with the total crack number. On the other side, the sizes of the final system matrices of either the NGF or the eigen COD algorithm remain unchanged. However, in the NGF algorithm, the size of the matrix to determine the complementary solution also increases. Much computional time is required for solving this matrix. In this way, the eigen COD algorithm has improved greatly in dealing with large number of crack problems, providing a newly numerical technique for multiple crack problems. Not only the accuracy and efficiency of computation can be

Table 5. The normalized I-SIF comparsions (N_L = 9,N_G = 13).

Table 6. The CPU time (s) of the three algorithms.

guaranteed, but also the overall properties and local details can be obtained. In conclusion, the numerical models of eigen COD BIEs realize the simulations for multiple crack problems with large quantity of cracks using ordinary desk-top computers.

6. Conclusion

A new algorithm based on the eigen COD BIEs is proposed to simulate solid materials with cracks in large number, where the eigen COD is defined as the crack opening displacements of a crack under the fictitious tractions loading on the crack surfaces. With the concept of eigen COD, multiple crack problems can be easily solved in an iterative fashion with a relatively small size of final system matrix compared with the Dual BIEs and the NGF, showing the practical significance of the present approach. Through the division of adjacent group and far-field group, the local Eshelby matrix, reflecting the strong interactions among cracks and avoiding numerical iteration divergence in the case of dense cracks, derived from the traction boundary integral equation in discrete form is introduced. Numerical examples verify the feasibility of the eigen COD BIEs for the simulation of multiple crack problems from two aspects of the numerical calculation accuracy and efficiency of computation. The numerical results of stress intensity factors show that the eigen COD algorithm has intrinsic scientific and rationality. Not only ensures the accuracy, also greatly improves the efficiency. The overall properties such as the rigidities and the local details such as the stresses around the crack tips can be achieved in future research, which have important theoretical significance and engineering application value.

Acknowledgments

This work is supported by the Funds Area of the National Natural Science Foundation of China (Grant No. 11662005), and the Science Foundation of Jiangxi Province (Grant No. 20202BABL201015).

The first author would like to respectively gratefully acknowledge Prof. Yijun Liu of the University of Cincinnati and his doctorial advisor Prof. Hang Ma of Shanghai University, for their advices in the research on boundary element method in modeling multiple crack problems.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Dolado, J.S. and Breugel, K.V. (2011) Recent Advances in Modeling for Cementitious Materials. Cement and Concrete Research, 41, 711-726. https://doi.org/10.1016/j.cemconres.2011.03.014
[2]	Wu, J.Y. and Xu, S.L. (2011) An Augmented Multicrack Elastoplastic Damage Model for Tensile Cracking. International Journal of Solids and Structures, 48, 2511-2528. https://doi.org/10.1016/j.ijsolstr.2011.05.001
[3]	Xiao, H.T. and Yue, Z.Q. (2011) A Three-Dimensional Displacement Discontinuity Method for Crack Problems in Layered Rocks. International Journal of Rock Mechanics and Mining Sciences, 48, 412-420. https://doi.org/10.1016/j.ijrmms.2011.02.005
[4]	Chmelik, F., Anton Trnik, A., Stubna, I. and Pesick, J. (2011) Creation of Microcracks in Porcelain during Firing. Journal of the European Ceramic Society, 31, 2205-2209. https://doi.org/10.1016/j.jeurceramsoc.2011.05.045
[5]	Hu, G.L., Ramesh, K.T., Cao, B.Y. and McCauley, J.W. (2011) The Compressive Failure of Aluminum Nitride Considered as a Model Advanced Ceramic. Journal of the Mechanics and Physics of Solids, 59, 1076-1093. https://doi.org/10.1016/j.jmps.2011.02.003
[6]	Sobelman, O.S., Gibeling, J.C., Stover, S.M., Hazelwood, S.J., Yeh, O.C., Shelton, D.R. and Martina, R.B. (2004) Do Microcracks Decrease or Increase Fatigue Resistance in Cortical Bone. Journal of Biomechanics, 37, 1295-1303. https://doi.org/10.1016/j.jbiomech.2003.12.034
[7]	Dhanasekar, M., Han, J.J. and Qin, Q.H. (2007) Interaction of Surface and Cracks in Railhead. Computer Assisted Mechanics and Engineering Sciences, 14, 79-89.
[8]	Chen, Y.Z. (2007) Integral Equation Methods for Multiple Crack Problems and Related Topics. Applied Mechanics Review, 60, 172-194. https://doi.org/10.1115/1.2750671
[9]	Brebbia, C.A., Telles, J.C.F. and Wrobel, L.C. (1984) Boundary Element Techniques—Theory and Applications in Engineering. Springer, Berlin. https://doi.org/10.1007/978-3-642-48860-3
[10]	Saez, A., Gallego, R. and Dominguez, J. (1995) Hypersingular Quarter-Point Boundary Elements for Crack Problems. International Journal for Numerical Methods in Engineering, 38, 1681-1701. https://doi.org/10.1002/nme.1620381006
[11]	Blandford, G.E., Ingraffea, A.R. and Liggett, J.A. (1981) Two-Dimensional Stress Intensity Factor Computation Using the Boundary Element Method. International Journal for Numerical Methods in Engineering, 17, 387-404. https://doi.org/10.1002/nme.1620170308
[12]	Aliabadi, M.H. (1997) Boundary Element Formulations in Fracture Mechanics. Applied Mechanics Reviews, 50, 83-96. https://doi.org/10.1115/1.3101690
[13]	Hong, H.K. and Chen, J.T. (1988) Derivations of Integral Equations of Elasticity. ASCE Journal of Engineering Mechanics, 114, 1028-1044. https://doi.org/10.1061/(ASCE)0733-9399(1988)114:6(1028)
[14]	Chen, J.T. and Hong, H.K. (1999) Review of Dual Boundary Element Methods with Emphasis on Hypersingular Integrals and Divergent Series. Applied Mechanics Reviews, 52, 17-33. https://doi.org/10.1115/1.3098922
[15]	Portela, A., Aliabadi, M.H. and Rooke, D.P. (1992) Dual Boundary Elements Analysis of Cracked Plates: Singularity Subtraction Technique. International Journal of Fracture, 55, 17-28. https://doi.org/10.1007/BF00018030
[16]	Greengard, L.F. and Rokhlin, V. (1987) A Fast Algorithm for Particle Simulations. Journal of Computational Physics, 73, 325-348. https://doi.org/10.1016/0021-9991(87)90140-9
[17]	Liu, Y.J. (2009) Fast Multipole Boundary Element Method—Theory and Applications in Engineering. Cambridge University Press, London. https://doi.org/10.1017/CBO9780511605345
[18]	Erdogan, F. (1983) Stress Intensity Factors. ASME Journal of Applied Mechanics, 50, 992-1002. https://doi.org/10.1115/1.3167212
[19]	Ang, W.T. and Clements, D.L. (1987) A Boundary Integral Equation Method for the Solution of a Class of Crack Problems. Journal of Elasticity, 17, 9-21. https://doi.org/10.1007/BF00042444
[20]	Telles, J.C.F., Castor, G.S. and Guimaraes, S. (1995) A Numerical Green’s Function Approach for Boundary Elements Applied to Fracture Mechanics. International Journal for Numerical Methods in Engineering, 38, 3259-3274. https://doi.org/10.1002/nme.1620381906
[21]	Telles, J.C.F. and Guimaraes, S. (2000) Green’s Function: A Numerical Generation for Fracture Mechanics Problems via Boundary Elements. Computer Methods in Applied Mechanics and Engineering, 188, 847-858. https://doi.org/10.1016/S0045-7825(99)00366-7
[22]	Kachanov, M. (1987) Elastic Solids with Many Cracks: A Simple Method of Analysis. International Journal of Solids and Structures, 23, 23-43. https://doi.org/10.1016/0020-7683(87)90030-8
[23]	Kachanov, M. (1993) Elastic Solids with Many Cracks and Related Problems. Advances in Applied Mechanics, 30, 259-445. https://doi.org/10.1016/S0065-2156(08)70176-5
[24]	Feng, X.Q. and Yu, S.W. (2010) Damage Micromechanics for Constitutive Relations and Failure of Microcracked Quasi-Brittle Materials. International Journal of Damage Mechanics, 19, 911-948. https://doi.org/10.1177/1056789509359662
[25]	Guo, Z. and Ma, H. (2011) Solution of Stress Intensity Factors of Multiple Cracks in Plane Elasticity with Eigen COD Formulation of Boundary Integral Equation. Journal of Shanghai University (English Edition), 15, 173-179. https://doi.org/10.1007/s11741-011-0716-1
[26]	Ma, H., Guo, Z., Dhanasekar, M., Yan, C. and Liu, Y.J. (2013) Efficient Solution of Multiple Cracks in Great Number Using Eigen COD Boundary Integral Equations with Iteration Procedure. Engineering Analysis with Boundary Elements, 37, 487-500. https://doi.org/10.1016/j.enganabound.2012.12.007
[27]	Ma, H., Guo, Z., Yan, C. and Dhanasekar, M. (2013) Numerical Solution of Stress Intensity Factors of Multiple Cracks in Great Number with Eigen COD Boundary Integral Equations. Australian Journal of Mechanical Engineering, 11, 1-10. https://doi.org/10.7158/M12-054.2013.11.1
[28]	Ma, H., Yan, C. and Qin, Q.H. (2009) Eigenstrain Formulation of Boundary Integral Equations for Modeling Particle-Reinforced Composites. Engineering Analysis with Boundary Elements, 33, 410-419. https://doi.org/10.1016/j.enganabound.2008.06.002
[29]	Ma, H., Fang, J.B. and Qin, Q.H. (2011) Simulation of Ellipsoidal Particle-Reinforced Materials with Eigenstrain Formulation of 3D BIE. Advances in Engineering Software, 42, 750-759. https://doi.org/10.1016/j.advengsoft.2011.05.013
[30]	Ma, H. and Qin, Q.H. (2007) Solving Potential Problems by a Boundary-Type Meshless Method—The Boundary Point Method Based on BIE. Engineering Analysis with Boundary Elements, 31, 749-761. https://doi.org/10.1016/j.enganabound.2007.03.001
[31]	Ma, H., Zhou, J. and Qin, Q.H. (2010) Boundary Point Method for Linear Elasticity Using Constant and Quadratic Moving Elements. Advances in Engineering Software, 41, 480-488. https://doi.org/10.1016/j.advengsoft.2009.10.006
[32]	Tada, H., Paris, P. and Irwin, G. (1985) The Stress Analysis of Cracks Handbook. Professional Engineering Publishers, Suffolk.
[33]	Sih, G.C. (1973) Handbook of Stress Intensity Factors. Lehigh University, Bethlehem.
[34]	Chen, W.H. and Chang, C.S. (1989) Analysis of Two-Dimensional Mixed-Crack Problems for Finite Element Alternating Method. Computers & Structures, 6, 1451-1458. https://doi.org/10.1016/0045-7949(89)90485-9
[35]	Wang, G.S. and Feng, X.T. (2001) The Interaction of Multiple Rows of Periodical Cracks. International Journal of Fracture, 110, 73-100. https://doi.org/10.1023/A:1010824317724