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.
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 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. 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 Figure 1, a finite region containing a crack is shown. The displacement boundary integral equation isin which and are Kelvin solutions of displacement type and traction type. and represent the source and field points. , if Ps locates at a smooth boundary, where is the Kronecker delta. The expression of and can be found in [11].
The traction boundary integral equation is
The expression of and can also be found in [11].
In the DBEM, the following relations are considered:in which , , and and represent tractions on the two crack surfaces, respectively. Considering equations (3)–(7), (1) and (2) can be rewritten intowhere, in the crack surface, the CODs are which is , in which and 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.
In 6DTEs, the geometric shape functions are different from the physical shape functions , 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 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 is in proportion to . From the work of Yates et al. [19], equation (10) obtained is an approximation formula, which omits high-order terms (), and two terms in the asymptotic expansion of the near-front relative crack-face displacement can be accurately captured, i.e., and r. In [19], it can also be found that the terms of order or 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:
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 and are included. This 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 are the collocation node coordinates of the element in the coordinate system. Inversing the matrix obtained by equation (12), the coefficients can be determined. Assume that ; the physical shape functions for the crack front element are given in Appendix B (equation (A.3)).
In the second case, the edge lies in the crack front. It can be found that is in proportion to . The shape functions for physical quantities should bewhere i = 0, …, 5. Inversing the matrix obtained by equation (13), the coefficients can be obtained. We assume ; the physical shape functions for the crack front element are given in Appendix C (equation (A.4)).
In the final case, the edge lies in the crack front. It can be found that is in proportion to . The shape functions for physical quantities should bewhere i = 0, …, 5. Inversing the matrix obtained by equation (14), the coefficients can be obtained. We assume ; 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:where ne1 and ne2 are the total number of elements on S or C+. denotes the collocation nodes, while Q is the inner field points. is the Jacobian. For 6DTEs, . For the crack front elements, .
Collocating Ps through all field points, the boundary conditions are imposed and the linear equations are rearranged intowhere contains the unknowns , , and 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 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 . 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.
(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 Figure 3(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 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 , around the source point P. In Figure 3(e), it can be found that 6 triangles can be obtained, i.e., , , , , , and . It is worth noting that if the angles around P are larger than , the corresponding triangle can be further subdivided into two subtriangles. Without loss of generality, assume that , taking the midpoint of the line segment as . can be divided into two triangles, i.e., and . 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 subpolygon for example; three subtriangles can be obtained, i.e., , , and . 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
To implement the singularity subtraction method, firstly we introduce a coordinate transformation from the intrinsic coordinate space to the integration space . 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 as follows:
Using polar transformation in the integration space , the singularity subtraction can be directly applied. It should be noted that, for each subtriangle around the source point, and It should be noted that the term in equation (13) possesses singularity when 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
In each subtriangle around the source point, the polar transformation is performed to eliminate the singularity of the integral, and the following equation can be obtained:where is the Jacobian from the space to the integration space .
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 , 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 around the source point, using the polar transformation in the integration space , the following equation can be obtained:
Using Taylor series expansions of kernel, shape functions, and the Jacobian of the transformation, the integrand can be expanded as
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. 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 that O 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 , whose direction points into the body, as shown in Figure 5.
The SIFs can be evaluated as r approaches zero:where , , , and 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 ( or ) and . is the distance between points P and P'. This is because that ( or ) may not be perpendicular to the crack front, where the CODs are evaluated at point P (such as P1, P2, and P3), as shown in Figure 6. , , and are the projections of 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 P1, P2, P3, and P′ are in a line. Firstly, we compute the distance , , and . Then, we calculate the SIFs using equation (27) through the CODs of P1, P2, and P3. We denote the results as , , and . Finally, using a linear extrapolation of and , we obtain the SIFs of P′.
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:
Considering singular integrals I in equation (24), firstly, assume S is an isosceles triangle with , , and being its three vertexes and . varies from to . In Figure 7, it can be found that as becomes larger, the computational accuracy decreases quickly from 10−5 to 10−1. When , the computational accuracy is about , while when , the computational accuracy is about . Thus, we choose as the demarcation point.
Then, we assume the three vertexes of the triangles are , , and , respectively. In this example, and b varies from 1 to 10. We fix . It can be found that, in Figure 8, as the values of b/a become larger, the computational accuracy decreases quickly; especially, and . 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.
6.2. A Penny-Shaped Crack in an Infinite Space under Uniform Loading
Firstly, we consider an embedded penny-shaped crack in Figure 9. The parameters of this problem are as follows: . In the z-direction, a far-field unit stress is imposed. Material constants are chosen as and . 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 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 calculated by our method is shown in Figure 12. Assuming the offset is 0.25, it can be found that the curves of 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.
(a)
(b)
(c)
6.3. Elliptical Crack in an Infinite Space
As shown in 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 is 0.25, when . 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 . In this paper, 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 .
The 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]. 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 , 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 , 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 , 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 , 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 , the relative error is 2.5% if we employ 126 elements, while the error is 2.0% when we use 42 more elements.
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. , , or . 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.
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.
Appendix
A. Shape Functions for Geometric Quantities and Physical Quantities in Normal Elements
The shape functions used for 6DTEs in Figure 2 are as follows:
The functional shape functions for 6DTEs in Figure 2 are as follows:
B. Shape Functions for Physical Quantities in the Crack Front Element for the First Case
C. Shape Functions for Physical Quantities in the Crack Front Element for the Second Case
D. Shape Functions for Physical Quantities in the Crack Front Element for the Third Case
Data Availability
The numerical result data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported partly by the National Natural Science Foundation of China (11602229, 11602082, and U1804141), partly by a project supported by the scientific research fund of Hunan Provincial Education Department (19B145), partly by the Key Scientific and Technological Project of Henan Province (192102210227), partly by the Scientific and Technological Innovation Team of Colleges and Universities in Henan Province in 2020 (20IRTSTHN015), and partly by the Key Discipline Laboratory Foundation of Noise and Vibration Control of Ship Equipment at Shanghai Jiao Tong University (VSN201903).