A New and Efficient Boundary Element-Free Method for 2-D Crack Problems

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 1/r 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 newvariables are directly approximatedwith 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.


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 [3][4][5][6], 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 [7][8][9].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. [13][14][15].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 1/ 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 domaintype element-free methods, such as element-free Galerkin method [17][18][19], can conveniently model the crack discontinuity with the enriched nodal degrees of freedom [20][21][22].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 [23][24][25].BIE can also be combined with the nonlocal approximations to generate the boundary-type elementfree methods, for example, the called boundary node method [26,27] and the called boundary element-free method [28][29][30], in which the unknowns are approximated with the MLS approximation [26], improving MLS approximation [28][29][30] 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 [31][32][33].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 1/ 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.

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   (,  = 1, 2).The governing equations are given as where   is the body force;  = (3 − )/( − 1); G is the shear modulus;  equals 3 − 4] for plane strain and (3 − ])/(1 + ]) for plane stress with ] being Poisson's ratio.Equation (2) can also be written as The fundamental solutions   (x 0 , x) and   (x 0 , x) are the th displacement and traction at the field point x( 1 ,  2 ) caused by a unit point force along the th direction at a source point x 0 ( 01 ,  02 ), and they can be expressed as follows [1,2,36]: where   is the Kronecker delta;  = |x − x 0 | is the distance between points x and x 0 ;   (x) 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, n − = −n + , the Somigliana formula expresses the displacement at an interior point x 0 ( 01 ,  02 ) in terms of the traction   (x) =     and displacement   (x) on the boundary point x( 1 ,  2 ) [1,12]: where (x) is the arc length along the boundary  or crack surface Γ; ∑   (x) =  +  (x) +  −  (x) is the sum of the tractions acting on the upper and lower crack surfaces; Δ  (x) =  +  (x) −  −  (x) is the difference of the displacements between the upper and lower crack surfaces.
Generally, ∑   (x) 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 [3][4][5].Let x 0 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 1/2 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 [7][8][9], 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 1/ 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   (x 0 ) with respect to  0 to give Using (  −  0 )/ 0 = −(  −  0 )/  , we have (x 0 , x) Substituting ( 9) into (4), we have here,   (x 0 , x) =   (x 0 , x) is used.The fundamental solution   (x 0 , x) should also satisfy the equilibrium equation, so we have the relationship Furthermore, we can obtain The last step of (13) resulted from the identity − 2 (/ 1 ) +  1 (/ 2 ) = / (see Figure 1) Similar procedure gives Substituting ( 8) into ( 2) or (4), using ( 11), (13), and ( 14), we have in which Applying integration by parts to the second and fourth integrals of (15) and then letting x 0 tend to the crack upper or lower surface, the limit of   (x 0 )  (x 0 ) gives where   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   ( + 0 ) −   ( − 0 ) 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 oncrack boundary, and ( 16) is applied to the crack surface.

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;   (x) 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 oncrack boundary is approximated as ,   (18) in which  , (x) is the derivative of the MLS shape function.
Referring to previous researches [13,15], the crack displacement dislocation can be approximated using where  1 and  2 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 [16][17][18].In MLS approximation, the trail function is written as where   (),  = 1, 2, . . .,  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 basis and the quadratic basis Unknown coefficients   () in ( 21) can be determined by minimizing the weighted discrete  2 norm where ( −   ) is the weight function with compact support;  is the number of nodes with ( −   ) > 0; and  I is the nodal parameter.The minimum of  in (24) with respect to a() leads to a set of linear equations where Thus, unknown coefficients a() can be obtained from (25) as Substituting ( 27) into ( 21), the MLS approximation can be written in a standard form as where the MLS shape function  I () is defined as with The derivatives of shape function can be obtained as in 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  0 in the formula

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" [16][17][18].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   (x 0 , x),   (x 0 , x),   (x 0 , x), and   (x 0 , x), 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 forms where () 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 forms Mathematical Problems in Engineering where   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 1/√  /2 −  at the crack tip.By using a simple transformation  = √  /2 − , 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 ∫  0 () ln(1/) and Cauchy singularity integral∫   (()/||) ( < 0,  > 0).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. 2 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 4 as [ (1)    (  ) +  (2)    (  )] where with   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  1 and  2 corresponding to the crack tip nodes based on (19) and (20).

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  0 is fixed as two times the distance between the adjacent nodes to ensure that four nodes are involved for the ordinary evaluation points.

Straight
Crack in a Square Plate.The first example is a square plate that contains a horizontal straight crack along the  1 -axis (Figure 3).The length of the crack is 2, and the origin of the Cartesian coordinate system  1  2 is located at the crack center.The plate is under a uniform tensile stress  along  2 .The plate dimension is fixed as / = 2.0.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    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].

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  2 direction.The crack length is 2, and the distance between two crack centers is 2.When the distance between the centers of two cracks is  = 2, the obtained normalized stress intensity factors are as follows:   =   /√ = 1.0259 and   =   /√ = 1.0466 (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.

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  , =  , /√ sin(), 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.

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.

Figure 1 :
Figure 1: The positive direction and normal direction for the outside boundary  and crack surface Γ.

Figure 2 :
Figure 2: The integral interval is divided into two segments for the evaluation of the singular integrals.

Figure 3 :
Figure 3: A horizontal straight crack in a square elastic plate.

Figure 4 :
Figure 4: An internal kinked crack in a rectangular plate.

Figure 5 :
Figure 5: Two collinear cracks in an infinite elastic plate subjected to the remote uniform tension stress.

Table 1 :
The normalized stress intensity factors   =   /√ for a square plate.

Table 2 :
The normalized stress intensity factors for the crack tips .

Table 3 :
The normalized stress intensity factors for the crack tips .
= +/ √ 2. The kink of crack is at the center of plate, which has a height equal to twice the width; that is, ℎ = 2.The plate is under a uniform tensile stress  along  2 ./ 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