Abstract

An efficient boundary element-free method is established for 2-D crack problems by combining a pair of boundary integral equations and the moving-least square approximation. The displacement boundary integral equation is collated on the on-crack boundary, and a new traction boundary integral equation is applied on the crack surface without the separate consideration of the upper and lower sides. In virtue of integration by parts, only singularity in order is involved in the integral kernels of new traction boundary integral equation, which brings convenience to the numerical implementation. Meanwhile, the integration by parts produces the new variables, the displacement density, and displacement dislocation density, and they are the coexisting unknowns along with the displacement and displacement dislocation. With the high-order continuity of the moving-least square approximation, these new variables are directly approximated with the nodal displacement or displacement dislocation, and the final system of equations contains the unknowns of nodal displacements and displacement dislocations only. The boundary element-free computational scheme is established, and several examples show the efficiency and flexibility of the proposed method.

1. Introduction

In comparison with the domain-type methods such as finite element method (FEM), the boundary integral equation (BIE) method show some special advantages for the fracture mechanics problems because it only requires the discretization of the general boundary and crack surface rather than the domain. But the traditional displacement BIE cannot directly be applied to crack bodies because the geometrical overlapping of the upper and lower crack surfaces leads to an indeterminacy of equations [1, 2]. To overcome this difficulty, a direct way is to derive the traction BIE of cracked bodies [36], but the kernels of the traction BIE involve hypersingular terms and their evaluation is very arduous. A different approach is the dual BIE method, in which the displacement BIE is applied to the outside boundary and to one side of crack surface, and the traction BIE is applied to another side of crack surface [79]. To reduce the number of equations in the final system matrix, Pan and Amadei proposed a new pair of BIEs so that the displacement BIE is collocated only on the no-crack boundary, and the traction BIE is applied only on one side of crack surface [10, 11]. But the evaluation of hypersingular integrals is still a vital task in the above methods. Actually, the hypersingular integral can be avoided by using integration by parts, and this skill has been used to derive a regular new traction BIE by Chau and Wang [12]. This skill has also been applied to the anisotropic medium and bimaterials by Sun et al. [1315]. In this paper, we propose that an alternative computational scheme by combining the new traction BIE and displacement BIE, that is, the displacement BIE, is collocated on the no-crack boundary, and the traction BIE only with the singular integral kernel in order of is collocated along the entire crack surface (no need for the separate discretization of the upper and lower crack surfaces).

In the derivation of the new traction BIE, two new variables, the displacement density and displacement dislocation density, are introduced, and they are the basic coexisting unknowns with the displacement and displacement dislocation in the proposed computational scheme. So the new technique is needed to efficiently bridge these unknowns in the numerical simulation. Mesh-free/element-free methods emerge from the nonlocal interpolation/approximation techniques such as the moving-least square (MLS) approximation [16] that has no dependence on the elements. The domain-type element-free methods, such as element-free Galerkin method [1719], can conveniently model the crack discontinuity with the enriched nodal degrees of freedom [2022]. With the enlightenment of element-free methods, the extended FEM method has been proposed for the crack problems through the embedding of displacement jumps within elements [2325]. BIE can also be combined with the nonlocal approximations to generate the boundary-type element-free methods, for example, the called boundary node method [26, 27] and the called boundary element-free method [2830], in which the unknowns are approximated with the MLS approximation [26], improving MLS approximation [2830] or Shepard and Taylor interpolation [27]. Another advantage of element-free methods is that the shape function stratifies the higher-order continuity automatically, and the derivative of shape function can be constructed easily, which brings convenience to the higher-order continuity problems [3133]. For example, for the strain gradient structures, the strain gradient can be approximated directly with nodal displacement, and nodal displacements are the only unknowns in the final system matrix [33]. This advantage encourages us to establish a boundary element-free method to implement the numerical simulation, in which the displacement density and displacement dislocation density are directly approximated with the nodal displacement and displacement dislocation.

This paper will firstly present the new traction BIE only with singularity in order by using integration by parts [12] and then combine it to the displacement BIE to construct a pair of BIEs for 2-D crack problems. The displacement BIE is collocated on the on-crack boundary, and the unknowns are the nodal displacements and displacement dislocations. The traction BIE is collocated on crack surface, and the unknown displacement density and displacement dislocation density are directly approximated with the nodal displacements and displacement dislocations, so that the final system of equations contains the nodal displacements and displacement dislocations as the only unknowns. Since the shape functions are digitally known in the present method and the singular integrals must be evaluated in a pure numerical scheme, the method of Torino [34] is employed to efficiently evaluate the Cauchy integrals.

2. New Boundary Integral Equations

Consider small displacement in two-dimensional, homogeneous, isotropic, and linear-elastic solids, the displacement, strain, and stress are, respectively, denoted as , , and (). The governing equations are given aswhere is the body force; ; G is the shear modulus; equals for plane strain and for plane stress with being Poisson’s ratio. Equation (2) can also be written as

The fundamental solutions and are the th displacement and traction at the field point caused by a unit point force along the th direction at a source point , and they can be expressed as follows [1, 2, 36]:where is the Kronecker delta; is the distance between points and ; is the unit normal along the boundary.

Consider a crack embedding in the elastic body as Figure 1, if its upper and lower surfaces are viewed as being geometrically overlapping, that is, , the Somigliana formula expresses the displacement at an interior point in terms of the traction and displacement on the boundary point [1, 12]:where is the arc length along the boundary or crack surface ; is the sum of the tractions acting on the upper and lower crack surfaces; is the difference of the displacements between the upper and lower crack surfaces.

Generally, in (6) is equal to zero, but it is not directly omitted at present, and this will lead to a different equation in comparison with the general traction boundary integral equation [35]. Let trend to the no-crack boundary point; we obtain the following boundary integral equation:where is the coefficient that depends only upon the local boundary geometry, and it is equal to for the smooth boundary.

It is well known, however, that for a cracked body, (7) does not have a unique solution because of the geometry singularity of the crack surface [2, 10, 11]. To overcome this difficulty, the traction integral equation is generally needed to be supplemented. In the dual BIE method [79], the displacement BIE is collocated on the no-crack boundary and on one side of crack surface while the traction integral equation is collocated on the other side of crack surface. To reduce the number of equations in the final system matrix, Pan and Amadei proposed a new pair of BIEs, so that the displacement integral equation is collocated only on the no-crack boundary, and a different traction integral equation is applied only on one side of the crack surface [10, 11]. In the above methods, the integral kernels of traction BIE involve the hypersingular terms and their evaluation is very arduous. In the following, we will present a new traction integral equation [12] that involves the singular integral kernel only in order and then combine it with (7) to solve the crack problems.

Since (6) is valid at all internal points of the linear-elastic body, we can differentiate with respect to to give Using , we haveSubstituting (9) into (4), we havehere, is used.

The fundamental solution should also satisfy the equilibrium equation, so we have the relationshipFurthermore, we can obtainThe last step of (13) resulted from the identity (see Figure 1)

Similar procedure givesSubstituting (8) into (2) or (4), using (11), (13), and (14), we havein which and is 2-D permutation tensor; that is, and .

Applying integration by parts to the second and fourth integrals of (15) and then letting tend to the crack upper or lower surface, the limit of giveswhere is the coefficient similar to the one in (7) and is the difference of the tractions between the upper and lower crack surfaces.

Equation (16) does not involve the hypersingular integral kernels, and it is very helpful for the numerical computation. On the other hand, the integration by parts produces the new variables and with the physical meaning: the displacement densities and dislocation densities. The crack boundary conditions on the upper and lower surface have been incorporated into (16), so there is no need to discretize the upper and lower crack surfaces separately in the numerical implementation. Although, in (6) is generally equal to zero, the traction difference may not be zero and may need to be considered, so it is vital to not directly omit in (6). In the following boundary element-free method, the displacement BIE (7) is applied to on-crack boundary, and (16) is applied to the crack surface.

3. The MLS Approximation for the Unknowns

The MLS approximation is a flexible nonlocal interpolation technique. Particularly, the MLS approximation can produce the higher-order continuity shape function automatically, and the derivative of unknowns can directly be approximated with the nodal parameters. While (7) is used, the following approximations are used:where and are the nodal parameters; is the MLS shape function; is the number of nodes which are covered by the effect domain of the evaluated point.

While (16) is used, the displacement density on the on-crack boundary is approximated asin which is the derivative of the MLS shape function.

Referring to previous researches [13, 15], the crack displacement dislocation can be approximated usingwhere and are the arc coordinates of two crack tips and is the nodal parameter. The displacement dislocation density can thus be approximated as

The MLS approximation is summarily described here, and the related details can be learned from the related literatures [1618]. In MLS approximation, the trail function is written aswhere , are monomial basis functions; is the number of terms in the basis; and are coefficients of the basis functions. Examples of commonly used bases are the linear basisand the quadratic basisUnknown coefficients in (21) can be determined by minimizing the weighted discrete normwhere is the weight function with compact support; is the number of nodes with ; and is the nodal parameter. The minimum of in (24) with respect to leads to a set of linear equationswhereThus, unknown coefficients can be obtained from (25) asSubstituting (27) into (21), the MLS approximation can be written in a standard form aswhere the MLS shape function is defined aswith

The derivatives of shape function can be obtained asin which

In the present work, the quadratic basis is used, and the weight function is chosen as [16]The compact support domain of the evaluation point is determined by choosing the suitable quantity in the formula

4. Boundary Element-Free Method

To carry out the numerical integration, the outside boundary and crack surface are firstly separated into a series of subdomains that are also called “integral cells” [1618]. Some nodes are then selected on each cell. The choice of cells is arbitrary and independent on nodes, but each cell must contain enough nodes to ensure that all compact supports domains cover the total boundary.

In this paper, a crack is uniformly collocated nodes and is simultaneously separated into integral cells uniformly. For the outside boundary, the integral cells are arranged in accordance with the nodal collocation; that is, two adjacent nodes denote an integral cell. The source point is chosen as the midpoint of integral cells in turn after the substitution of (17)–(20) into (7) and (16), and the numerical integrals are performed on every cell, which gives a series of linear equations. The solution of linear equations gives the nodal parameters, and the nodal displacements, displacement dislocations, displacement densities, and displacement dislocation densities can be evaluated with (17)–(20).

The integral cells are only used so that the integral can be evaluated numerically. The cell is not a boundary element, and the shape function is independent of it. The present method is called the boundary element-free method.

Consider the expressions of , , , and , when the source point is located inside an integral cell that does not contain any crack tip; we need to evaluate the singular integrals in the formswhere is the smooth function on the integral cells.

If the integral cell includes a crack tip at its left or right end, we need to evaluate the singular integrals in the formswhere is the length of the integral cells that include the crack tip and is the smooth function on the integral cells. Two integrals in (38) show the weak singularity at the crack tip. By using a simple transformation , this weak singularity can be eliminated and these two integrals can be transformed as the integrals in forms of (37); that is, we need to find the appropriate approaches to evaluate the logarithmic singularity integral and Cauchy singularity integral.

Since the shape functions are digitally known in the present method, the singular integrals must be evaluated in a pure numerical scheme. The logarithmic singularity integral can be numerically evaluated using Gauss quadrature rule of logarithmic functions [37]. For the Cauchy singular integral, we find the method of Torino [34] to be very effective, and it has been used in our previous researches [13, 38]. This method is simply described as follows.

As shown in Figure 2, the singular pole is located on the integral interval from to (assuming that point is nearer to than to ). A suitable point is selected so that point is equidistant from and b. point Gauss ruler points are used to draw the segments from to and from to b, and the total integral is then evaluated with the degree of exactness aswherewith and being the standard Gauss quadrature point and weight, respectively.

Around a crack tip, the relation between the crack opening displacements and the extended stress intensity factors is [11, 13] where is the distance ahead of the crack tip. So and can be evaluated with and corresponding to the crack tip nodes based on (19) and (20).

5. Numerical Examples

Four numerical examples are presented in this section. The first one is to verify the convergence and stability, and the other examples are to check the efficiency and advantages of the present method. During the numerical simulations, the scale is fixed as two times the distance between the adjacent nodes to ensure that four nodes are involved for the ordinary evaluation points.

5.1. Straight Crack in a Square Plate

The first example is a square plate that contains a horizontal straight crack along the -axis (Figure 3). The length of the crack is , and the origin of the Cartesian coordinate system is located at the crack center. The plate is under a uniform tensile stress along . The plate dimension is fixed as . The computation indicates that good results can be obtained while 80 nodes are used on the outside boundary. Table 1 gives normalized stress intensity factors for various numbers of nodes collocated on crack surface. It can be seen that a stable convergence can be achieved, and the error is just 0.3% in comparison with [35] while 20 nodes are used for crack surface.

5.2. A Kinked Crack in a Rectangular Plate

Now, consider a rectangular plate with an internal kinked crack (Figure 4). The horizontal segment has the length a, and the other segment has the length and makes an angle of with the horizontal one. The horizontal projection of the total crack is . The kink of crack is at the center of plate, which has a height equal to twice the width; that is, . The plate is under a uniform tensile stress along . is fixed as 0.1, and is taken as 0.2, 0.4, and 0.6. The normalized stress intensity factors are defined as and . 120 nodes are used on the outside boundary, and 30 nodes are used for the horizontal segment and 6, 12, or 18 nodes are used, respectively, for the inclined segment in the cases of = 0.2, 0.4, or 0.6. The stress intensity factors for crack tips and are listed in Tables 2 and 3. The results show a good agreement with those in [8].

5.3. Two Collinear Cracks in an Infinite Plate

This example studies the interaction between two collinear cracks with the same length in an infinite elastic plate (Figure 5), and the plate is subjected to the uniform remote tension along the direction. The crack length is , and the distance between two crack centers is . When the distance between the centers of two cracks is , the obtained normalized stress intensity factors are as follows: and (20 integral cells are uniformly collocated on each crack surface), and these results are very close to the reported results [35]. Figure 6 plots and versus , which indicates that the interaction between two collinear cracks is considerably large only when they are very close to each other.

5.4. Circular-Arc Crack in an Unbounded Domain

To demonstrate the validity of the present method for the curved crack, a circular-arc crack in an infinite elastic medium is considered (Figure 7). The crack region is under uniform tension in two perpendicular directions. The normalized stress intensity factors are denoted as , and they are plotted as a function of the semiangle in Figure 7. The present results are very close to the report results [35], and it is shown that the method is very effective for the curved crack since the equation is built on the arc coordinate system.

6. Conclusions

The integration by parts is used to obtain a new traction BIE in which the hypersingularity integral is avoided, and it is combined with the displacement BIE to construct a pair of BIEs for the crack problems. The proposed method brings convenience to the numerical implementation. The negative consequence of integration by parts is the introduction of two new variables, the displacement density and displacement dislocation density, and they are the coexisting unknowns along with the displacement and displacement dislocation. In virtue of the high-order continuity of the MLS approximation, we propose to approximate the new variables directly with the nodal displacement and displacement dislocation, so that nodal displacement density and displacement dislocation density are not contained in the system of equations of the established element-free method. The method is proved to be very efficient for the general crack problems, and it is also shown that the MLS approximation can be well used in the problems that have the higher-order continuity requirement.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by National Natural Science Foundation of China [Grant no. 11472316], Plan for Scientific Innovation Talent of Henan Province, China [Grant no. 164200510020], and Natural Science Foundations from Department of Education of Henan Province, China [Grant nos. 13B560304 and 15A560045].