# A Generation of Special Triangular Boundary Element Shape Functions for 3D Crack Problems

Authors: Guizhong Xie; Fenglin Zhou
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper focuses on tackling the two drawbacks of the dual boundary element method (DBEM) when solving crack problems with a discontinuous triangular element: low accuracy of the calculation of integrals with singularity and crack front element must be utilized to model the square-root property of displacement. In order to calculate the integrals with higher order singularity, the triangular elements are segmented into several subregions which consist of subtriangles and subpolygons. The singular integrals in those subtriangles are handled by the singularity subtraction technique in the integration space and can be regularized and accurately calculated. For the nearly singular integrals in those subpolygons, the element subdivision technique is employed to improve the calculation accuracy. In addition, considering the location of the crack front in the element, special crack front elements are constructed based on a 6-node discontinuous triangular element, in which the displacement extrapolation method is introduced to obtain the stress intensity factors (SIFs) without consideration of orthogonalization of the crack front mesh. Several numerical results are investigated to fully verify the validation of the presented approach.


---

## Body

## 1. Introduction

Accurate computation of the SIFs is of great importance for analysis of 3D fracture mechanical problems. For such problems, it is necessary to accurately calculate the SIFs, which can characterize the fracture property. The major difficulty in the calculation of SIFs is approximation of the displacement distribution nearby the crack front.In order to tackle this problem, numerical methods including the finite element method (FEM) and the boundary element method (BEM) have been widely applied [1–4]. In the application of the FEM or BEM in crack problems, special crack front elements are usually employed to model the square-root distribution of the displacements nearby the crack front. Many special crack front elements have been defined to capture the asymptotic behavior of a specific node in the FEM [1]. In the 3D BEM, however, the available crack front elements are only of quadrilateral type, such as 8-node crack front elements proposed by Mi and Aliabadi [5], 9-node crack front elements in the pioneering work of Li et al. [6], and also crack front elements in the work of Pan and Yuan [7]. The crack front elements of the triangular type have not been found by authors [8, 9]. Due to the general geometric adaptability of the triangular mesh, a special triangular crack front element of which one edge lies in the crack front for the BEM is also necessary.For solutions of crack problems, the dual boundary element method (DBEM), which is originated from the BEM, is a more general and efficient method than other extension methods [6, 10–15] in the range of BEMs, such as the multidomain BEM [16], the symmetric Galerkin boundary element method [17, 18], and the displacement discontinuity method [11]. In this method, crack front elements are widely employed. They are of great importance to the accuracy and efficiency of the method [14, 19]. In this paper, a type of special discontinuous triangular crack front element is constructed to model the square-root distribution of the displacement nearby the crack front. It should be noted that, in Frangi et al.’s pioneer work [18], they have already constructed a series of mixed elements including both the quadrilateral element and the triangular element and obtained accurate results in a symmetric Galerkin BEM framework. In their work, however, the elements employed in the crack front are also quadrilateral elements. Generally, the quadrilateral element could provide a more accurate interpolation result than the triangular one. In the case of the crack of complicated shape, however, the crack surface can be more easily meshed by the triangular elements rather than by the quadrilateral ones. In other words, the triangular element could be generated more easily with high quality and could be applied to model the crack front with arbitrary shape. Furthermore, it is believed that the computed stress intensity factors (SIFs), which are usually computed through an extrapolation method, in the triangular element are less than those in the quadrilateral one.To improve the computational accuracy of integrals with high-order singularities, the integral patch is divided into several subpatches which consist of subtriangles and subpolygons at first. Then, the subtriangles are mapped into normalized isosceles right triangles. In the normalized space, a polar coordinate transformation is finally performed to cancel the singularity of the integral. With these three steps, the high-order singular integrals can be calculated accurately. The integrals in other types of subpatches can also be computed through an element subdivision method, in which the subpolygons are subdivided into several triangular patches.In the computation of the SIFs, a crack opening displacement (COD) extrapolation method is adopted without consideration of orthogonalization of the mesh nearby the crack front.The remainder of this paper is organized as follows: In Section2, the displacement and traction boundary integral equations which are involved in the DBEM are briefly introduced. A 6-node discontinuous triangular element modeling strategy is introduced in Section 3. Section 4 mainly describes the regularization of singular integrals. The computational method of SIFs is detailed in Section 5 which is followed by some benchmark testing numerical examples in Section 6. This paper ends with conclusions in Section 7.


## 2. The Displacement and Traction Boundary Integral Equations

In Figure1, a finite region containing a crack is shown. The displacement boundary integral equation is(1)cijPSujPS=∫S+C++C−uij∗PS,QtjQdSQ−∫S+C++C−tij∗PS,QujQdSQ,in which uij∗ and tij∗ are Kelvin solutions of displacement type and traction type. PSx,y,z and Qx,y,z represent the source and field points. cij=δij/2, if Ps locates at a smooth boundary, where δij is the Kronecker delta. The expression of uij∗ and tij∗ can be found in [11].Figure 1
A finite region containing a crack.The traction boundary integral equation is(2)cijtjPS=niPS∫S+C++C−Uijk∗PS,QtkQdSQ−niPS∫S+C++C−Tijk∗PS,QukQdSQ.The expression ofUijk∗ and Tijk∗ can also be found in [11].In the DBEM, the following relations are considered:(3)tkQ+=−tkQ−,(4)uij∗PS,Q+=uij∗PS,Q−,(5)tij∗PS,Q+=−tij∗PS,Q−,(6)Uijk∗PS,Q+=Uijk∗PS,Q−,(7)Tijk∗PS,Q+=−Tijk∗PS,Q−,in which Q+∈C+, Q−∈C−, and tkQ+ and tkQ− represent tractions on the two crack surfaces, respectively. Considering equations (3)–(7), (1) and (2) can be rewritten into(8)cijPSujPS=∫Suij∗PS,QtjQdS−∫Stij∗PS,QujQdS−∫C+tij∗PS,Q+ΔujQdS,(9)tjPC+=niPC+∫SUijk∗PC+,QtkQdSQ−niPC+∫STijk∗PC+,QukQdSQ−niPC+∫C+Tijk∗PC+,QΔukQdSQ,where, in the crack surface, the CODs are ΔukQ which is ukQ+−ukQ−, in which ukQ+ and ukQ− represent displacements on the two crack surfaces. Equations (8) and (9) are collocated on noncrack surfaces and crack surfaces, respectively.


## 3. 6-Node Discontinuous Triangular Elements

Due to the existence of Cauchy and Hadamard principle value integrals in the DBEM [20], the shape function derivatives of elements on crack surfaces must be Holder continuous. Thus, 6-node discontinuous triangular elements (6DTEs) are employed in the crack surfaces, while nearby the crack front, we applied some specially constructed elements. The 6-node discontinuous triangular element is illustrated in Figure 2, in which k is the offset parameter, varying from 0.05 to 0.4 in our implementation.Figure 2
6-node discontinuous triangular element.In 6DTEs, the geometric shape functionsϕgeoi are different from the physical shape functions ϕcolli, which are both listed in Appendix A (equations (9) and (A.2)).For the crack front elements, special shape functions with asymptotic properties are constructed to approximate the distribution of the displacement in the vicinity of the crack front. In the first case, the edge, the formulation of which isξ+η=1 in the parametric space, of the special crack front element lies in the crack front.From the geometric meaning of area coordinates about the triangular element, we can find thatr=Pξ,η−Qξ,1−ξ is in proportion to 1−ξ−η. From the work of Yates et al. [19], equation (10) obtained is an approximation formula, which omits high-order terms (or), and two terms in the asymptotic expansion of the near-front relative crack-face displacement can be accurately captured, i.e., r and r. In [19], it can also be found that the terms of order r3/2 or or3/2 are omitted to obtain the relations between stress intensity factors and displacements in the near-front field. Thus, from those pioneering works in [17], it can be found that the crack opening displacements should have the following form:(10)Δuξ,η=L2¯ξr+L1¯ξr+L0¯ξ=L2ξ1−ξ−η+L1ξ1−ξ−η+L0ξ.Based on the idea that is similar to that in the pioneer work of Li et al. [6], we construct the shape function through equation (10) in which 1−ξ−η and 1−ξ−η are included. This type of shape function is designed to approximate the physical variables nearby the crack front.(11)Mi=a1i+a2iξ+a3iξ2+a4i1−ξ−η+a5iξ1−ξ−η+a6i1−ξ−η,in which i = 0, …, 5. In order to determine the value of the coefficients in equation (11), the following condition is imposed:(12)Miξj,ηj=δij,i=0,…,5,in which ξj,ηj are the collocation node coordinates of the element in the ξ,η coordinate system. Inversing the 6×6 matrix obtained by equation (12), the coefficients aji can be determined. Assume that k=0.25; the physical shape functions for the crack front element are given in Appendix B (equation (A.3)).In the second case, the edgeη=0 lies in the crack front. It can be found that r=Pξ,η−Qξ,0 is in proportion to η. The shape functions for physical quantities should be(13)Mi=a1i+a2iξ+a3iξ2+a4iη+a5iξη+a6iη,where i = 0, …, 5. Inversing the 6×6 matrix obtained by equation (13), the coefficients aji can be obtained. We assume k=0.25; the physical shape functions for the crack front element are given in Appendix C (equation (A.4)).In the final case, the edgeξ=0 lies in the crack front. It can be found that r=Pξ,η−Q0,η is in proportion to ξ. The shape functions for physical quantities should be(14)Mi=a1i+a2iη+a3iη2+a4iξ+a5iηξ+a6iξ,where i = 0, …, 5. Inversing the 6×6 matrix obtained by equation (14), the coefficients aji can be obtained. We assume k=0.25; the physical shape functions for the crack front element are given in Appendix D (equation (A.5)).Substituting the special shape functions into equations (8) and (9), the following discretized boundary integral equations can be obtained:(15)cijPSujPS=∑e=1ne1∑α=05tjαQα∫Seuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη−∑e=1ne1∑α=05ujαQα∫Setij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη−∑e=1ne2∑α=05ΔujαQα∫Ce+tij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη,(16)tjPC+=niPC+∑e=1ne1∑α=05tkαQα∫SeUijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη−niPC+∑e=1ne1∑α=05ukαQα∫SeTijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη−niPC+∑e=1ne2∑α=05ΔukαQα∫Ce+Tijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη,where ne1 and ne2 are the total number of elements on S or C+. Qα denotes the collocation nodes, while Q is the inner field points. Jξ,η is the Jacobian. For 6DTEs, Nαξ,η=ϕCollα. For the crack front elements, Nαξ,η=Mα.CollocatingPs through all field points, the boundary conditions are imposed and the linear equations are rearranged into(17)Ax=f,where x contains the unknowns ui, ti, and Δui and f is obtained by the given boundary quantities.


## 4. Regularization of Singular Integrals

To achieve high computational accuracy, the integrals with singularity in the DBEM must be accurately computed [9, 12, 15, 19–27]. In this section, the singular integral is gradually computed through three steps. First, the considered integral patch is divided into several subpatches in which both subtriangles and subpolygons are involved. Second, the subpolygon patches are further subdivided into triangular patches and the subtriangles are normalized into some isosceles right triangles. Finally, in the normalized space, a classical singularity subtraction method [15, 19–21] based on the polar coordinate transformation is implemented to eliminate the singularity of the integral. After these three steps of operation, the Cauchy principle value integral and the Hadamard finite part integral can be accurately calculated. A more detailed discussion about these two types of integrals can be found in [9, 15].


### 4.1. Triangular Element Subdivision Technique

To cope with the singular integrals, the element subdivision is necessary in the 3D BEM [9, 19–25]. Moreover, a suitable integral area is of great importance to improving the accuracy of singular integral calculation, especially for Cauchy or Hadamard finite part integrals. Thus, in this part, considering the location of the source point x in the element, the element can be segmented into several subregions which contain subtriangles and subpolygons. Adaptive integration based on element subdivisions for the integrals on subpolygons is employed to improve the calculation accuracy of the singularity subtraction method.Considering the location of the source point, the source point might locate near the corner node or edge node, as shown in Figures3(a) or 3(b). It can be found that if the triangular element is directly subdivided into three subtriangles, one angle of subtriangles may belong to 2π/3,π. Thus, the accuracy of numerical results will become poor which can be found in Section 6.1. This will cause large errors in the integral calculation over the singular integration elements. In order to improve the calculation accuracy, an element subdivision method, in which the source points, corner points, and edge middle points are all considered the vertexes in the subdivided triangular elements, is applied, as shown in Figures 3(c) or 3(d). As illustrated in numerical examples, however, the computational accuracy is still very sensitive to the value of the offset parameter k.Figure 3
Subdivisions of the triangular element according to the location of the source point.
(a)
(b)
(c)
(d)
(e)Furthermore, to overcome this drawback, an adaptive patch subdivision scheme for an arbitrary triangular element has been developed. In this scheme, the triangular element is firstly divided into several triangular and polygon patches around the source point, as illustrated in Figure3(e). The integrals over these triangles around the source point are singular and hypersingular. Then, the polygons are further subdivided into several subtriangles, over which the integrals are treated as nearly singular integrals or regular integrals. The element subdivision scheme can be described as follows.Firstly, we should draw parallel lines of three edges of the triangle through the source pointP, compute the distances in the Cartesian coordinate system from P to the intersection point between the parallel lines and the triangular element edges, and find the minimum distance d.Secondly, we construct isosceles subtriangles, the edge length of which isd1=0.5d, around the source point P. In Figure 3(e), it can be found that 6 triangles can be obtained, i.e., ΔPI1I2, ΔPI2I3, ΔPI3I4, ΔPI4I5, ΔPI5I6, and ΔPI6I1. It is worth noting that if the angles around P are larger than 2π/3, the corresponding triangle can be further subdivided into two subtriangles. Without loss of generality, assume that I2PI32π/3, taking the midpoint of the line segment I2I3¯ as I23. ΔPI2I3 can be divided into two triangles, i.e., ΔPI2I23 and ΔPI23I3. This step is to ensure fine shape of the subtriangles, which is important for accurate calculation of the singular integrals, especially for hypersingular integrals. In each triangular patch, the singular integrals are calculated by the singularity subtraction method.Finally, we construct additional subtriangles which the source point is not located in. Take subpolygonI2I3O3CO2 for example; three subtriangles can be obtained, i.e., I2O2C, I2CI3, and I3CO3. The integrals over these subtriangular patches are treated as nearly singular integrals. The calculation accuracy is also very important for the final results. These integrals can be accurately calculated by the method in [21].


### 4.2. Mapping a Subtriangle from the Local Intrinsic Coordinatesξ,η to an Isosceles Right Triangle in the Integration Space a,b

To implement the singularity subtraction method, firstly we introduce a coordinate transformation from the intrinsic coordinateξ,η space to the integration space a,b. As shown in Figure 4, any subtriangle for singular or hypersingular integrals in the ξ,η space can be transformed into an isosceles right triangle in the integration space a,b as follows:(18)ξ=1−a−bξ0+aξ1+bξ2,η=1−a−bη0+aη1+bη2.Figure 4
A coordinate transformation from theξ,η space to the integration space a,b.Using polar transformation in the integration spacea,b, the singularity subtraction can be directly applied. It should be noted that, for each subtriangle Δm around the source point, θ1,θ2=0,π/2 and 0,Rθ=0,2/cosθ−π/4. It should be noted that the term ∫Seuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη in equation (13) possesses singularity when PS locates in the integration element, and singular integrals arise. The most common methods to remove the weak singularity are based on variable transformations [9, 19, 20, 24]. In our application, the polar transformation is applied. Combined with element subdivision, the integrals with weak singularity fall into two parts: singular integrals and nearly singular integrals, as in(19)∫Seuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη=∑m6∫SΔmuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη+∫Se−∑m6SΔmuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη.In each subtriangleΔm around the source point, the polar transformation is performed to eliminate the singularity of the integral, and the following equation can be obtained:(20)∫SΔmuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη=∫SΔmuij∗PS,Qα,βNαa,bJa,bJaba,bdadb=∫0π/4∫02/cosθ−π/4uij∗PS,Qρ,θNαρ,θJρ,θJabρ,θρdρdθ,where Jaba,b is the Jacobian from the ξ,η space to the integration space a,b.


### 4.3. Evaluation of the Cauchy Principle Value Integral and the Hadamard Finite Part Integral

Employing the simple solution method and the polar transformation, calculation of singular integrals in the DBEM can be solved when the collocation point is not on the crack surface. However, if the collocation point belongs to the crack surfaces, the integrands with high-order singularity should be regularized. The singular term of the integrand is separated in an analytical or a semianalytical form. The remaining part is regular, which can be directly calculated. In equation (16), the expression of the Hadamard finite part integral is ∫Ce+Tijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη, which is similar to equation (19); the polar transformation is applied in our method. Combined with element subdivision, the integrals with high-order singularity fall into two parts: singular and nearly singular integrals, as found in(21)∫Ce+Tijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη=∑m6∫SΔmTijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη+∫Ce+−∑m6SΔmTijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη.For each subtriangleΔm around the source point, using the polar transformation in the integration space a,b, the following equation can be obtained:(22)∫SΔmTijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη=∫SΔmTijk∗PC+,Qξ,ηNαa,bJa,bJaba,bdadb=∫0π/4∫02/cosθ−π/4Tijk∗PC+,Qρ,θNαρ,θJρ,θJabρ,θρdρdθ.Using Taylor series expansions of kernel, shape functions, and the Jacobian of the transformation, the integrandFρ,θ=Tijk∗PC+,Qρ,θNαρ,θJρ,θJabρ,θρ can be expanded as(23)Fρ,θ=F−2θρ2+F−1θρ+O1.Using equation (23), equation (22) can be transformed into the following form:(24)∫SΔmTijk∗Ps+,Qξ,ηNαξ,ηJξ,ηdξdη=∫0π/4∫02/cosθ−π/4Fρ,θ−F−2θρ2−F−1θρdρdθ+∫0π/4∫02/cosθ−π/4F−2θρ2+F−1θρdρdθ.In the right-hand side of equation (24), the first integral is regular and can be integrated accurately by the standard Gaussian quadrature method. The second integral has the same singularity as the original integrand with a simpler form which can be integrated analytically or semianalytically. More details of the singularity subtraction technique can be found in [15].


## 4.1. Triangular Element Subdivision Technique

To cope with the singular integrals, the element subdivision is necessary in the 3D BEM [9, 19–25]. Moreover, a suitable integral area is of great importance to improving the accuracy of singular integral calculation, especially for Cauchy or Hadamard finite part integrals. Thus, in this part, considering the location of the source point x in the element, the element can be segmented into several subregions which contain subtriangles and subpolygons. Adaptive integration based on element subdivisions for the integrals on subpolygons is employed to improve the calculation accuracy of the singularity subtraction method.Considering the location of the source point, the source point might locate near the corner node or edge node, as shown in Figures3(a) or 3(b). It can be found that if the triangular element is directly subdivided into three subtriangles, one angle of subtriangles may belong to 2π/3,π. Thus, the accuracy of numerical results will become poor which can be found in Section 6.1. This will cause large errors in the integral calculation over the singular integration elements. In order to improve the calculation accuracy, an element subdivision method, in which the source points, corner points, and edge middle points are all considered the vertexes in the subdivided triangular elements, is applied, as shown in Figures 3(c) or 3(d). As illustrated in numerical examples, however, the computational accuracy is still very sensitive to the value of the offset parameter k.Figure 3
Subdivisions of the triangular element according to the location of the source point.
(a)
(b)
(c)
(d)
(e)Furthermore, to overcome this drawback, an adaptive patch subdivision scheme for an arbitrary triangular element has been developed. In this scheme, the triangular element is firstly divided into several triangular and polygon patches around the source point, as illustrated in Figure3(e). The integrals over these triangles around the source point are singular and hypersingular. Then, the polygons are further subdivided into several subtriangles, over which the integrals are treated as nearly singular integrals or regular integrals. The element subdivision scheme can be described as follows.Firstly, we should draw parallel lines of three edges of the triangle through the source pointP, compute the distances in the Cartesian coordinate system from P to the intersection point between the parallel lines and the triangular element edges, and find the minimum distance d.Secondly, we construct isosceles subtriangles, the edge length of which isd1=0.5d, around the source point P. In Figure 3(e), it can be found that 6 triangles can be obtained, i.e., ΔPI1I2, ΔPI2I3, ΔPI3I4, ΔPI4I5, ΔPI5I6, and ΔPI6I1. It is worth noting that if the angles around P are larger than 2π/3, the corresponding triangle can be further subdivided into two subtriangles. Without loss of generality, assume that I2PI32π/3, taking the midpoint of the line segment I2I3¯ as I23. ΔPI2I3 can be divided into two triangles, i.e., ΔPI2I23 and ΔPI23I3. This step is to ensure fine shape of the subtriangles, which is important for accurate calculation of the singular integrals, especially for hypersingular integrals. In each triangular patch, the singular integrals are calculated by the singularity subtraction method.Finally, we construct additional subtriangles which the source point is not located in. Take subpolygonI2I3O3CO2 for example; three subtriangles can be obtained, i.e., I2O2C, I2CI3, and I3CO3. The integrals over these subtriangular patches are treated as nearly singular integrals. The calculation accuracy is also very important for the final results. These integrals can be accurately calculated by the method in [21].


## 4.2. Mapping a Subtriangle from the Local Intrinsic Coordinatesξ,η to an Isosceles Right Triangle in the Integration Space a,b

To implement the singularity subtraction method, firstly we introduce a coordinate transformation from the intrinsic coordinateξ,η space to the integration space a,b. As shown in Figure 4, any subtriangle for singular or hypersingular integrals in the ξ,η space can be transformed into an isosceles right triangle in the integration space a,b as follows:(18)ξ=1−a−bξ0+aξ1+bξ2,η=1−a−bη0+aη1+bη2.Figure 4
A coordinate transformation from theξ,η space to the integration space a,b.Using polar transformation in the integration spacea,b, the singularity subtraction can be directly applied. It should be noted that, for each subtriangle Δm around the source point, θ1,θ2=0,π/2 and 0,Rθ=0,2/cosθ−π/4. It should be noted that the term ∫Seuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη in equation (13) possesses singularity when PS locates in the integration element, and singular integrals arise. The most common methods to remove the weak singularity are based on variable transformations [9, 19, 20, 24]. In our application, the polar transformation is applied. Combined with element subdivision, the integrals with weak singularity fall into two parts: singular integrals and nearly singular integrals, as in(19)∫Seuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη=∑m6∫SΔmuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη+∫Se−∑m6SΔmuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη.In each subtriangleΔm around the source point, the polar transformation is performed to eliminate the singularity of the integral, and the following equation can be obtained:(20)∫SΔmuij∗PS,Qξ,ηNαξ,ηJξ,ηdξdη=∫SΔmuij∗PS,Qα,βNαa,bJa,bJaba,bdadb=∫0π/4∫02/cosθ−π/4uij∗PS,Qρ,θNαρ,θJρ,θJabρ,θρdρdθ,where Jaba,b is the Jacobian from the ξ,η space to the integration space a,b.


## 4.3. Evaluation of the Cauchy Principle Value Integral and the Hadamard Finite Part Integral

Employing the simple solution method and the polar transformation, calculation of singular integrals in the DBEM can be solved when the collocation point is not on the crack surface. However, if the collocation point belongs to the crack surfaces, the integrands with high-order singularity should be regularized. The singular term of the integrand is separated in an analytical or a semianalytical form. The remaining part is regular, which can be directly calculated. In equation (16), the expression of the Hadamard finite part integral is ∫Ce+Tijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη, which is similar to equation (19); the polar transformation is applied in our method. Combined with element subdivision, the integrals with high-order singularity fall into two parts: singular and nearly singular integrals, as found in(21)∫Ce+Tijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη=∑m6∫SΔmTijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη+∫Ce+−∑m6SΔmTijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη.For each subtriangleΔm around the source point, using the polar transformation in the integration space a,b, the following equation can be obtained:(22)∫SΔmTijk∗PC+,Qξ,ηNαξ,ηJξ,ηdξdη=∫SΔmTijk∗PC+,Qξ,ηNαa,bJa,bJaba,bdadb=∫0π/4∫02/cosθ−π/4Tijk∗PC+,Qρ,θNαρ,θJρ,θJabρ,θρdρdθ.Using Taylor series expansions of kernel, shape functions, and the Jacobian of the transformation, the integrandFρ,θ=Tijk∗PC+,Qρ,θNαρ,θJρ,θJabρ,θρ can be expanded as(23)Fρ,θ=F−2θρ2+F−1θρ+O1.Using equation (23), equation (22) can be transformed into the following form:(24)∫SΔmTijk∗Ps+,Qξ,ηNαξ,ηJξ,ηdξdη=∫0π/4∫02/cosθ−π/4Fρ,θ−F−2θρ2−F−1θρdρdθ+∫0π/4∫02/cosθ−π/4F−2θρ2+F−1θρdρdθ.In the right-hand side of equation (24), the first integral is regular and can be integrated accurately by the standard Gaussian quadrature method. The second integral has the same singularity as the original integrand with a simpler form which can be integrated analytically or semianalytically. More details of the singularity subtraction technique can be found in [15].


## 5. Calculation of SIFs

Using the asymptotic behavior of displacement and elasticity solutions obtained by the DBEM, the SIFs can be evaluated by results of any point nearby the crack front. Assume thatO is an arbitrary point at the crack front. We construct a local coordinate system in which the center is O. The unit tangent vector of the crack front at point O is t, the unit normal of the crack surface is b, and n=b×t/b×t, whose direction points into the body, as shown in Figure 5.Figure 5
Local coordinate system at the crack front.The SIFs can be evaluated asr approaches zero:(25)KI=E41−v2limr⟶0π2rΔub,KII=E41−v2limr⟶0π2rΔun,KIII=E41+vlimr⟶0π2rΔut,where Δub=Δu⋅b, Δun=Δu⋅n, Δut=Δu⋅t, and v and E are Poisson’s ratio and Young’s modulus. Using the discontinuous elements, CODs are obtained at the collocation nodes inside the element. As shown in Figure 6, employing the one-point formula, SIFs are expressed as(26)KIP′=E41−v2π2rcosϕΔubP,KIIP′=E41−v2π2rcosϕΔunP,KIIIP′=E41+vπ2rcosϕΔutP.Figure 6
Calculation of the SIFs at the crack front.In equation (26), ϕ is the angle of P1P′ (P2P′ or P3P′) and n. r is the distance between points P and P'. This is because that P1P′ (P2P′ or P3P′) may not be perpendicular to the crack front, where the CODs ΔuP are evaluated at point P (such as P1, P2, and P3), as shown in Figure 6. ΔubP, ΔunP, and ΔutP are the projections of ΔuP in the coordinate directions of the local crack front coordinate system, as shown in Figure 5.In this method, the COD extrapolation applied for calculating SIFs is as follows: As found in Figure6, points P1, P2, P3, and P′ are in a line. Firstly, we compute the distance rP1=P1P′, rP2=P2P′, and rP3=P3P′. Then, we calculate the SIFs using equation (27) through the CODs of P1, P2, and P3. We denote the results as Kp1, Kp2, and Kp3. Finally, using a linear extrapolation of Kp2 and Kp3, we obtain the SIFs of P′.(27)KP′=rP3KP2−rP2KP3rP3−rP2.


## 6. Numerical Examples

### 6.1. Singular Integrals on Triangles Considering Vertex Angles and Aspect Ratios of Two Edge Lengths

In this example, the following integral is computed:(28)I=∫S1rdS.Considering singular integralsI in equation (24), firstly, assume S is an isosceles triangle with 0,0, 1,0, and cosθ,sinθ being its three vertexes and a=b=1. θ varies from 4π/9 to 8π/9. In Figure 7, it can be found that as θ becomes larger, the computational accuracy decreases quickly from 10−5 to 10−1. When θ=2π/3, the computational accuracy is about 1.9×10−3, while when θ=13π/18, the computational accuracy is about 5.0×10−3. Thus, we choose θ=2π/3 as the demarcation point.Figure 7
Errors of singular integrals on isosceles triangles with differentθ values.Then, we assume the three vertexes of the triangles are0,0, 1,0, and bcosθ,bsinθ, respectively. In this example, a=1 and b varies from 1 to 10. We fix θ=iπ/6,i=1,…,5. It can be found that, in Figure 8, as the values of b/a become larger, the computational accuracy decreases quickly; especially, θ>2π/3 and b/a>5. Thus, the distortion and aspect ratio of special crack elements will also affect the accuracy of the computed stress intensity factors along the crack front. In order to obtain stress intensity factors with high accuracy, the fine shape of elements should be ensured.Figure 8
Errors of singular integrals on triangles with different ratios ofb/a.


### 6.2. A Penny-Shaped Crack in an Infinite Space under Uniform Loading

Firstly, we consider an embedded penny-shaped crack in Figure9. The parameters of this problem are as follows: a=1.0. In the z-direction, a far-field unit stress is imposed. Material constants are chosen as E=1.0 and v=0.25. The crack surface is discretized into 58 elements in Figure 10, including 16 crack front elements. The results of our method (element subdivision in Figure 3(e)) with those obtained by element subdivision in Figures 3(c) and 3(d) combined with the singularity subtraction technique are compared first to test the efficiency of our method. The results in Table 1 show that the results of our method are more accurate and less affected by the offset parameter of the collocation nodes. In Table 1, method A represents the singularity subtraction method with element subdivision in Figure 3(e), while method B represents the singularity subtraction method with element subdivision in Figures 3(c) and 3(d). The results with different offset parameters from 0.05 to 0.4 agree well with the exact solutions. The largest error is within 1%. The exact solutions of Δu are given in [28, 29]. As shown in Figure 11, the results of our method are in good agreement with the exact solutions with different offset parameters. The largest error arises when the offset parameter k is 0.05, which is about 2%. The convergence of Δu calculated by our method is shown in Figure 12. Assuming the offset is 0.25, it can be found that the curves of Δu obtained by our method are nearly indistinguishable from the exact solutions when the number of elements varies from 58 to 244. And the largest error of our method is within 0.8%. It should be noted that if we use element subdivision in Figure 3(a) and 3(b), the results are very poor. Thus, we compare the results in the case of Figure 3(a) to those in the cases of Figures 3(c) and 3(d). From Table 2, we note that, in Liu’s method, the SIF results converge within 2% of the analytical solutions very quickly using 864 elements and then improve slowly afterwards [30]. The result with the finest mesh (4704 elements) still has a relative error of 1.66%. However, in our method, with the help of special crack front elements, the results agree well with the analytical solutions. This numerical example is illustrated clearly. The relative error in our method is much lower than that in the comparison work [30]. The shape of the element may affect the accuracy of the integral. But the influence is not quite significant.Figure 9
A penny-shaped crack in an infinite space.Figure 10
Boundary element discretization of a penny-shaped crack in an infinite space: (a) 58 elements; (b) 162 elements; (c) 244 elements.
(a)
(b)
(c)Table 1
Comparisons with other methods with differentk; the normalized SIF (KI/πa) is 0.6366.
k
SIFs by methodA
SIFs by methodB
0.05
0.6372
Wrong
0.1
0.6349
Wrong
0.15
0.6358
Wrong
0.2
0.6367
0.5755
0.25
0.6373
0.6554
0.3
0.6372
0.6401
0.4
0.6366
0.6380Figure 11
Comparisons ofΔu over the crack surface with different offset parameters using 58 elements.Figure 12
Convergence ofΔu calculated by our method when k = 0.25.Table 2
Comparison with Liu’s method for different number of mesh.
Methods
Mesh
Errors (%)
Method of Liu [30]
96
4.63∗
384
2.60∗
864
1.96∗
1536
1.92∗
2400
1.69∗
3456
1.69∗
4704
1.66
Our method
58
0.11
162
0.20
244
0.17
∗Values are obtained from the figures in [30], which may have errors.


### 6.3. Elliptical Crack in an Infinite Space

As shown in Figure13, the elliptical crack in an infinite solid is subjected to uniform inclined traction. Young’s Modulus (E) is 1, and Poisson’s ratio v is 0.25, when γ=π/4,w=0. 6-node triangular elements are used to discretize the crack surface, and special crack front elements proposed in our paper are also employed. The elliptic equation is x2/a2+y2/b2=1,a≥b. In this paper, a/b varies from 2 to 5. We compute SIFs of three modes at different sample points along the crack front. The position of the crack front is denoted by the angle ϕ from the positive direction of the largest radius of the ellipse to the negative direction of the largest radius of the ellipse. ϕ varies from 0 to π.Figure 13
An elliptical crack in an infinite solid under inclined traction [28].The convergence of the SIFs computed by our method is studied with different numbers of elements. In this example, we assume the offset parameterk = 0.25. We treat the solutions of Irwin et al. as the exact solutions [28]. From Figures 14–18, it can be found that the results of our method agree well with the exact solutions. And SIFs of three modes obtained by our method show a high convergence rate as the number of elements increases. In Figure 14, for the case a/b=2, when 74 elements are used, it can be found that the largest error arises when ϕ approaches zero. The largest error is about 1.5%. And the largest error is about 1% when the element number is 108. In Figure 15, for a/b=2.5, the error is 1.7% when 96 elements are employed and the error is about 1% using 206 elements. In Figure 16, for the case a/b=3, the error is about 1.9% employing 72 elements and the error is about 1.6% employing 172 elements. In Figure 17, in the case a/b=4, the error is 2.0% when 128 elements are employed, while the error is reduced to 1.5% when 212 elements are employed. In Figure 18, in the case a/b=5, the relative error is 2.5% if we employ 126 elements, while the error is 2.0% when we use 42 more elements.Figure 14
SIFs along the front of the elliptical crack in the finite space whena = 2b.Figure 15
SIFs along the front of the elliptical crack in the finite space whena = 2.5b.Figure 16
SIFs along the front of the elliptical crack in the finite space whena = 3b.Figure 17
SIFs along the front of the elliptical crack in the finite space whena = 4b.Figure 18
SIFs along the front of the elliptical crack in the finite space whena = 5b.


### 6.4. A Penny-Shaped Crack in an Infinite Space under Nonuniform Loading

In this example, an embedded penny-shaped crack under nonuniform loading is considered [31]. The geometry can be found in Figure 9, and three mesh models can also be found in Figure 10(a). Three nonuniform loadings are considered. σzz=−y, σzz=−2xy, or σzz=y2−x2. The stress intensity factors along the front of the penny-shaped crack can be found in Figure 19. And the largest errors are all within 3%. In order to test the accuracy, we also list the results in Tables 3–5.Figure 19
SIFs along the front of a penny-shaped crack in the finite space under nonuniform loadings.Table 3
Stress intensity factors along the crack front under nonuniform loadings compared with analytical solutions [31] for the case σzz=−y.
ϕ
Our method
Analytical solutions
Relative errors (%)
0.1608
0.06834
0.06797
0.54
0.3566
0.14788
0.14816
0.19
0.5532
0.22185
0.22298
0.51
0.7497
0.28726
0.2892
0.67
0.9455
0.34171
0.34412
0.70
1.1422
0.38287
0.38602
0.82
1.3388
0.40941
0.41304
0.88Table 4
Stress intensity factors along the crack front under nonuniform loadings compared with analytical solutions [31] for the case σzz=−2xy.
ϕ
Our method
Analytical solutions
Relative errors (%)
0.1608
0.10464
0.10735
2.52
0.3566
0.21707
0.22214
2.28
0.5532
0.29645
0.30356
2.34
0.7497
0.33104
0.33866
2.25
0.9455
0.31525
0.32226
2.18
1.1422
0.25091
0.25669
2.25
1.3388
0.14878
0.15196
2.09Table 5
Stress intensity factors along the crack front under nonuniform loadings compared with analytical solutions [31] for the case σzz=y2−x2.
ϕ
Our method
Analytical solutions
Relative errors (%)
0.1608
0.31486
0.32211
2.25
0.3566
0.25127
0.25678
2.15
0.5532
0.1488
0.15209
2.16
0.7497
0.02395
0.02423
1.16
0.9455
−0.10432
−0.1069
2.41
1.1422
−0.21689
−0.22224
2.41
1.3388
−0.29664
−0.30363
2.30


## 6.1. Singular Integrals on Triangles Considering Vertex Angles and Aspect Ratios of Two Edge Lengths

In this example, the following integral is computed:(28)I=∫S1rdS.Considering singular integralsI in equation (24), firstly, assume S is an isosceles triangle with 0,0, 1,0, and cosθ,sinθ being its three vertexes and a=b=1. θ varies from 4π/9 to 8π/9. In Figure 7, it can be found that as θ becomes larger, the computational accuracy decreases quickly from 10−5 to 10−1. When θ=2π/3, the computational accuracy is about 1.9×10−3, while when θ=13π/18, the computational accuracy is about 5.0×10−3. Thus, we choose θ=2π/3 as the demarcation point.Figure 7
Errors of singular integrals on isosceles triangles with differentθ values.Then, we assume the three vertexes of the triangles are0,0, 1,0, and bcosθ,bsinθ, respectively. In this example, a=1 and b varies from 1 to 10. We fix θ=iπ/6,i=1,…,5. It can be found that, in Figure 8, as the values of b/a become larger, the computational accuracy decreases quickly; especially, θ>2π/3 and b/a>5. Thus, the distortion and aspect ratio of special crack elements will also affect the accuracy of the computed stress intensity factors along the crack front. In order to obtain stress intensity factors with high accuracy, the fine shape of elements should be ensured.Figure 8
Errors of singular integrals on triangles with different ratios ofb/a.


## 6.2. A Penny-Shaped Crack in an Infinite Space under Uniform Loading

Firstly, we consider an embedded penny-shaped crack in Figure9. The parameters of this problem are as follows: a=1.0. In the z-direction, a far-field unit stress is imposed. Material constants are chosen as E=1.0 and v=0.25. The crack surface is discretized into 58 elements in Figure 10, including 16 crack front elements. The results of our method (element subdivision in Figure 3(e)) with those obtained by element subdivision in Figures 3(c) and 3(d) combined with the singularity subtraction technique are compared first to test the efficiency of our method. The results in Table 1 show that the results of our method are more accurate and less affected by the offset parameter of the collocation nodes. In Table 1, method A represents the singularity subtraction method with element subdivision in Figure 3(e), while method B represents the singularity subtraction method with element subdivision in Figures 3(c) and 3(d). The results with different offset parameters from 0.05 to 0.4 agree well with the exact solutions. The largest error is within 1%. The exact solutions of Δu are given in [28, 29]. As shown in Figure 11, the results of our method are in good agreement with the exact solutions with different offset parameters. The largest error arises when the offset parameter k is 0.05, which is about 2%. The convergence of Δu calculated by our method is shown in Figure 12. Assuming the offset is 0.25, it can be found that the curves of Δu obtained by our method are nearly indistinguishable from the exact solutions when the number of elements varies from 58 to 244. And the largest error of our method is within 0.8%. It should be noted that if we use element subdivision in Figure 3(a) and 3(b), the results are very poor. Thus, we compare the results in the case of Figure 3(a) to those in the cases of Figures 3(c) and 3(d). From Table 2, we note that, in Liu’s method, the SIF results converge within 2% of the analytical solutions very quickly using 864 elements and then improve slowly afterwards [30]. The result with the finest mesh (4704 elements) still has a relative error of 1.66%. However, in our method, with the help of special crack front elements, the results agree well with the analytical solutions. This numerical example is illustrated clearly. The relative error in our method is much lower than that in the comparison work [30]. The shape of the element may affect the accuracy of the integral. But the influence is not quite significant.Figure 9
A penny-shaped crack in an infinite space.Figure 10
Boundary element discretization of a penny-shaped crack in an infinite space: (a) 58 elements; (b) 162 elements; (c) 244 elements.
(a)
(b)
(c)Table 1
Comparisons with other methods with differentk; the normalized SIF (KI/πa) is 0.6366.
k
SIFs by methodA
SIFs by methodB
0.05
0.6372
Wrong
0.1
0.6349
Wrong
0.15
0.6358
Wrong
0.2
0.6367
0.5755
0.25
0.6373
0.6554
0.3
0.6372
0.6401
0.4
0.6366
0.6380Figure 11
Comparisons ofΔu over the crack surface with different offset parameters using 58 elements.Figure 12
Convergence ofΔu calculated by our method when k = 0.25.Table 2
Comparison with Liu’s method for different number of mesh.
Methods
Mesh
Errors (%)
Method of Liu [30]
96
4.63∗
384
2.60∗
864
1.96∗
1536
1.92∗
2400
1.69∗
3456
1.69∗
4704
1.66
Our method
58
0.11
162
0.20
244
0.17
∗Values are obtained from the figures in [30], which may have errors.


## 6.3. Elliptical Crack in an Infinite Space

As shown in Figure13, the elliptical crack in an infinite solid is subjected to uniform inclined traction. Young’s Modulus (E) is 1, and Poisson’s ratio v is 0.25, when γ=π/4,w=0. 6-node triangular elements are used to discretize the crack surface, and special crack front elements proposed in our paper are also employed. The elliptic equation is x2/a2+y2/b2=1,a≥b. In this paper, a/b varies from 2 to 5. We compute SIFs of three modes at different sample points along the crack front. The position of the crack front is denoted by the angle ϕ from the positive direction of the largest radius of the ellipse to the negative direction of the largest radius of the ellipse. ϕ varies from 0 to π.Figure 13
An elliptical crack in an infinite solid under inclined traction [28].The convergence of the SIFs computed by our method is studied with different numbers of elements. In this example, we assume the offset parameterk = 0.25. We treat the solutions of Irwin et al. as the exact solutions [28]. From Figures 14–18, it can be found that the results of our method agree well with the exact solutions. And SIFs of three modes obtained by our method show a high convergence rate as the number of elements increases. In Figure 14, for the case a/b=2, when 74 elements are used, it can be found that the largest error arises when ϕ approaches zero. The largest error is about 1.5%. And the largest error is about 1% when the element number is 108. In Figure 15, for a/b=2.5, the error is 1.7% when 96 elements are employed and the error is about 1% using 206 elements. In Figure 16, for the case a/b=3, the error is about 1.9% employing 72 elements and the error is about 1.6% employing 172 elements. In Figure 17, in the case a/b=4, the error is 2.0% when 128 elements are employed, while the error is reduced to 1.5% when 212 elements are employed. In Figure 18, in the case a/b=5, the relative error is 2.5% if we employ 126 elements, while the error is 2.0% when we use 42 more elements.Figure 14
SIFs along the front of the elliptical crack in the finite space whena = 2b.Figure 15
SIFs along the front of the elliptical crack in the finite space whena = 2.5b.Figure 16
SIFs along the front of the elliptical crack in the finite space whena = 3b.Figure 17
SIFs along the front of the elliptical crack in the finite space whena = 4b.Figure 18
SIFs along the front of the elliptical crack in the finite space whena = 5b.


## 6.4. A Penny-Shaped Crack in an Infinite Space under Nonuniform Loading

In this example, an embedded penny-shaped crack under nonuniform loading is considered [31]. The geometry can be found in Figure 9, and three mesh models can also be found in Figure 10(a). Three nonuniform loadings are considered. σzz=−y, σzz=−2xy, or σzz=y2−x2. The stress intensity factors along the front of the penny-shaped crack can be found in Figure 19. And the largest errors are all within 3%. In order to test the accuracy, we also list the results in Tables 3–5.Figure 19
SIFs along the front of a penny-shaped crack in the finite space under nonuniform loadings.Table 3
Stress intensity factors along the crack front under nonuniform loadings compared with analytical solutions [31] for the case σzz=−y.
ϕ
Our method
Analytical solutions
Relative errors (%)
0.1608
0.06834
0.06797
0.54
0.3566
0.14788
0.14816
0.19
0.5532
0.22185
0.22298
0.51
0.7497
0.28726
0.2892
0.67
0.9455
0.34171
0.34412
0.70
1.1422
0.38287
0.38602
0.82
1.3388
0.40941
0.41304
0.88Table 4
Stress intensity factors along the crack front under nonuniform loadings compared with analytical solutions [31] for the case σzz=−2xy.
ϕ
Our method
Analytical solutions
Relative errors (%)
0.1608
0.10464
0.10735
2.52
0.3566
0.21707
0.22214
2.28
0.5532
0.29645
0.30356
2.34
0.7497
0.33104
0.33866
2.25
0.9455
0.31525
0.32226
2.18
1.1422
0.25091
0.25669
2.25
1.3388
0.14878
0.15196
2.09Table 5
Stress intensity factors along the crack front under nonuniform loadings compared with analytical solutions [31] for the case σzz=y2−x2.
ϕ
Our method
Analytical solutions
Relative errors (%)
0.1608
0.31486
0.32211
2.25
0.3566
0.25127
0.25678
2.15
0.5532
0.1488
0.15209
2.16
0.7497
0.02395
0.02423
1.16
0.9455
−0.10432
−0.1069
2.41
1.1422
−0.21689
−0.22224
2.41
1.3388
−0.29664
−0.30363
2.30


## 7. Conclusions

A novel discontinuous triangular crack front element for calculating SIFs was presented in this paper. In our numerical implementation, the crack surface was modeled with 6-node discontinuous triangular elements which can satisfy the existence of the finite part integrals. And an element subdivision technique for the singular element was introduced. Combined with the element subdivision technique and singularity subtraction technique in the integration space, the high-order singular integrals can be calculated with high accuracy. Furthermore, special crack front elements were constructed. These specially constructed elements can successfully approximate the distribution of displacements nearby the crack front. And the SIFs can be obtained by a COD extrapolation method without consideration of orthogonalization of the crack front mesh. The SIFs obtained by the present method were in very good agreement with previously published results.



---

*Source: 4629761-2020-01-24.xml*