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

0is 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.0e 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.


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. e 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][2][3][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]. e 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][11][12][13][14][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. ey 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. en, 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. e 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. e remainder of this paper is organized as follows: In Section 2, 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. e computational method of SIFs is detailed in Section 5 which is followed by some benchmark testing numerical examples in Section 6. is paper ends with conclusions in Section 7.

The Displacement and Traction Boundary Integral Equations
In Figure 1, a finite region containing a crack is shown. e displacement boundary integral equation is in which u * ij and t * ij are Kelvin solutions of displacement type and traction type. P S (x, y, z) and Q(x, y, z) represent the source and field points. c ij � δ ij /2, if P s locates at a smooth boundary, where δ ij is the Kronecker delta. e expression of u * ij and t * ij can be found in [11]. The traction boundary integral equation is (2) e expression of U * ijk and T * ijk can also be found in [11]. In the DBEM, the following relations are considered: in which Q + ∈ C + , Q − ∈ C − , and t k (Q + ) and t k (Q − ) represent tractions on the two crack surfaces, respectively. Considering equations (3)- (7), (1) and (2) can be rewritten into where, in the crack surface, the CODs are Δu k (Q) which is displacements on the two crack surfaces. Equations (8) and (9) are collocated on noncrack surfaces and crack surfaces, respectively.

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. us, 6-node discontinuous triangular elements (6DTEs) are employed in the crack surfaces, while nearby the crack front, we applied some specially constructed elements. e 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.
In 6DTEs, the geometric shape functions ϕ i geo are different from the physical shape functions ϕ i coll , 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 that From the work of Yates et al. [19], equation (10) obtained is an approximation formula, which omits high-order terms (o(r)), 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 r 3/2 or o(r 3/2 ) are omitted to obtain the relations between stress intensity factors and displacements in the near-front field. us, from those pioneering works in [17], it can be found that the crack opening displacements should have the following form: 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. is type of shape function is designed to approximate the physical variables nearby the crack front.
in which i � 0, . . ., 5. In order to determine the value of the coefficients in equation (11), the following condition is imposed: 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 a i j 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 |η|. e shape functions for physical quantities should be where i � 0, . . ., 5. Inversing the 6 × 6 matrix obtained by equation (13), the coefficients a i j 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(ξ, η) − Q(0, η)| is in proportion to |ξ|. e shape functions for physical quantities should be where i � 0, . . ., 5. Inversing the 6 × 6 matrix obtained by equation (14), the coefficients a i j 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: Figure 2: 6-node discontinuous triangular element.

Mathematical Problems in Engineering
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.
Collocating P s through all field points, the boundary conditions are imposed and the linear equations are rearranged into where x contains the unknowns u i , t i , and Δu i and f is obtained by the given boundary quantities.

Regularization of Singular Integrals
To achieve high computational accuracy, the integrals with singularity in the DBEM must be accurately computed [9,12,15,[19][20][21][22][23][24][25][26][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][20][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].

Triangular Element Subdivision Technique.
To cope with the singular integrals, the element subdivision is necessary in the 3D BEM [9,[19][20][21][22][23][24][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. us, 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 Figures 3(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, π).
us, the accuracy of numerical results will become poor which can be found in Section 6.1. is 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.
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 Figure 3(e). e integrals over these triangles around the source point are singular and hypersingular. en, the polygons are further subdivided into several subtriangles, over which the integrals are treated as nearly singular integrals or regular integrals. e element subdivision scheme can be described as follows. Firstly, we should draw parallel lines of three edges of the triangle through the source point P, 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 is d 1 � 0.5 d, around the source point P. In Figure 3(e), it can be found that 6 triangles can be obtained, i.e., ΔPI 1 I 2 , ΔPI 2 I 3 , ΔPI 3 I 4 , ΔPI 4 I 5 , ΔPI 5 I 6 , and ΔPI 6 I 1 . It is worth noting that if the angles around P are larger than 4 Mathematical Problems in Engineering Mathematical Problems in Engineering 2π/3, the corresponding triangle can be further subdivided into two subtriangles. Without loss of generality, assume that 〈I 2 PI 3 〉2π/3, taking the midpoint of the line segment I 2 I 3 as I 23 . ΔPI 2 I 3 can be divided into two triangles, i.e., ΔPI 2 I 23 and ΔPI 23 I 3 . is 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 subpolygon I 2 I 3 O 3 CO 2 for example; three subtriangles can be obtained, i.e., I 2 O 2 C, I 2 CI 3 , and I 3 CO 3 . e integrals over these subtriangular patches are treated as nearly singular integrals. e calculation accuracy is also very important for the final results. ese integrals can be accurately calculated by the method in [21]. 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:

Mapping a Subtriangle from the Local Intrinsic Coordinates (ξ, η) to an Isosceles Right Triangle in the Integration
Using polar transformation in the integration space (a, b), the singularity subtraction can be directly applied. It should be noted that, for each subtriangle Δ m around the source point, It should be noted that the term S e u * ij (P S , Q(ξ, η))N α (ξ, η)J(ξ, η)dξdη in equation (13) possesses singularity when P S locates in the integration element, and singular integrals arise. e 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 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: where J ab (a, b) is the Jacobian from the (ξ, η) space to the integration space (a, b).

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 highorder singularity should be regularized. e singular term of the integrand is separated in an analytical or a semianalytical form. e remaining part is regular, which can be directly calculated. In equation (16), the expression of the Hadamard finite part integral is C + e T * ijk (P C + , 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 For each subtriangle Δ m around the source point, using the polar transformation in the integration space (a, b), the following equation can be obtained: 6 Mathematical Problems in Engineering Using Taylor series expansions of kernel, shape functions, and the Jacobian of the transformation, the integrand Using equation (23), equation (22) can be transformed into the following form: In the right-hand side of equation (24), the first integral is regular and can be integrated accurately by the standard Gaussian quadrature method. e 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].

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 that O is an arbitrary point at the crack front. We construct a local coordinate system in which the center is O. e 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.
e SIFs can be evaluated as r approaches zero: Δu n , where Δu b � Δu · b, Δu n � Δu · n, Δu t � Δ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 In equation (26), ϕ is the angle of P 1 P ′ (P 2 P ′ or P 3 P ′ ) and n. r is the distance between points P and P'. is is because that P 1 P ′ (P 2 P ′ or P 3 P ′ ) may not be perpendicular to the crack front, where the CODs Δu P are evaluated at point P (such as P 1 , P 2 , and P 3 ), as shown in Figure 6. Δu P b , Δu P n , and Δu P t are the projections of Δu P 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 Figure 6, points P 1 , P 2 , P 3 , and P′ are in a line. Firstly, we compute the distance r P 1 � P 1 P ′ , r P 2 � P 2 P ′ , and r P 3 � P 3 P ′ . en, we calculate the SIFs using equation (27) through the CODs of P 1 , P 2 , and P 3 . We denote the results as K p 1 , K p 2 , and K p 3 . Finally, using a linear extrapolation of K p 2 and K p 3 , we obtain the SIFs of P′.

Singular Integrals on Triangles Considering Vertex Angles and Aspect Ratios of Two Edge Lengths.
In this example, the following integral is computed: Considering singular integrals I 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 . us, we choose θ � 2π/3 as the demarcation point. en, we assume the three vertexes of the triangles are (0, 0), (1, 0), and (b cos θ, b sin θ), 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. us, 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.

A Penny-Shaped Crack in an Infinite Space under Uniform
Loading. Firstly, we consider an embedded penny-shaped crack in Figure 9. e 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. e crack surface is discretized into 58 elements in Figure 10, including 16 crack front elements. e 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. e 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 results with different offset parameters from 0.05 to 0.4 agree well with the exact solutions. e largest error is within 1%. e 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. e largest error arises when the offset parameter k is 0.05, which is about 2%. e 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 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]. e 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.
is numerical example is illustrated clearly. e relative error in our method is much lower than that in the comparison work [30]. e shape of the element may affect the accuracy of the integral. But the influence is not quite significant. Figure 13, 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 c � π/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. e elliptic equation is x 2 /a 2 + y 2 /b 2 � 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. e 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 π. e convergence of the SIFs computed by our method is studied with different numbers of elements. In this example, we assume the offset parameter k � 0.25. We treat the solutions of Irwin et al. as the exact solutions [28].

Mathematical Problems in Engineering
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. e largest error is about 1.5%. And the largest error is about 1% when the element number is 108.      [30], which may have errors.
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.

A Penny-Shaped Crack in an Infinite
Space under Nonuniform Loading. In this example, an embedded penny- Exact solutions K I using 74 elements K II using 74 elements K III using 74 elements K I using 108 elements K II using 108 elements K III using 108 elements shaped crack under nonuniform loading is considered [31]. e geometry can be found in Figure 9, and three mesh models can also be found in Figure 10(a). ree nonuniform loadings are considered. σ zz � − y, σ zz � − 2xy, or σ zz � y 2 − x 2 . e 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.

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. ese 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. e SIFs obtained by the present method were in very good agreement with previously published results.