1.INTRODUCTION

Numerical methods are necessary to solve many fracture mechanics problems. The Finite Element Method (FEM) and the Boundary Element Method (BEM) have become very popular for the analysis of fracture mechanics in solids. The BEM has been widely used in recent years because it allows a very accurate stress analysis along crack faces and modelling of crack propagation without re-meshing (^{Aliabadi, 2002}). Under a linear elastic approach, fracture mechanics allows the determination of the stress field at the crack tip using stress intensity factors (SIFs) as a function of the crack geometry and loading.

The application of the BEM to fracture mechanics problems was initiated by Cruse through two works presented in 1970 and 1971 (^{Cruse, 1996}). These early works reported inaccurate SIF results (^{Aliabadi, 2002}). Later, (^{Cruse and Wilson, 1977}) implemented quarter-point elements to improve the accuracy of the BEM calculations, but the method had other difficulties for its application to crack problems. In general, these initial applications of the BEM to crack problems were limited by the fact that the two surfaces that form the crack were coplanar, generating a mathematical degeneration (^{Cruse, 1996}).

In the early nineties, (^{Portela et al., 1992}) for two dimensions, and (^{Mi and Aliabadi, 1992}) for 3D solids, proposed the Dual Boundary Element Method (DBEM) in which a displacement boundary integral equation (BIE) is applied on a surface of the crack and a traction BIE is applied on the other crack surface, thus avoiding the degeneracy in the Kelvin formulation found by (^{Cruse, 1996}).

From there, many works have been developed in the area, such as (^{dell'Erba and Aliabadi, 2001}) who developed a DBEM methodology to solve 3D thermo-elasticity problems using the J-integral for evaluating the SIF. (^{Dirgantara and Aliabadi, 2002}) used the DBEM to obtain mixed mode SIF values for cracked thin plates using crack surface displacement extrapolation and the J-integral technique. (^{Purbolaksono et al., 2012}) calculated the SIF in deformable plates using the DBEM and displacement extrapolation techniques. (^{Wen and Aliabadi, 2012}) developed an algorithm to model smooth curve cracks using the DBEM.

Several alternative methods have been proposed for calculating SIF values at the crack tip using the FEM or the BEM as primary methods to solve the linear elasticity problem, such as:

•.
Displacement Extrapolation (DE), which consists in extrapolating the numerical displacement field with the analytical solution to obtain the SIF. (^{Cruse and Wilson, 1978}) used the DE with quarter-point elements, obtaining reasonable results.

•.
Strain Energy Release, which calculates the strain energy of the deformed body or the external work done by loads for small crack advances in order to differentiate it and extract the SIF. This method was used by (^{Cruse, 1988}) but was computationally expensive due to the small crack advance needed to achieve a reasonable accuracy.

•.
J-integral, a path-independent integral proposed by (^{Rice, 1968}), which is a contour integral that measures the strain energy flux across its boundary. This technique has been used to compute the SIF in many FEM and BEM works, including (^{Rigby and Aliabadi, 1998}) who proposed a decomposition technique to extract the mixed mode SIF, (^{Bezerra and Medeiros, 2002}) who proposed an alternative numerical scheme to implement the J-integral calculation, (^{Ortiz and Cisilino, 2006}) who developed a J-integral based methodology for 3D cracks.

•.
M-integral, also called interaction integral, a variant of the J-integral used by (^{Walters et al., 2005}) with the Galerkin BEM to obtain the SIF for 3D curved loaded cracks.

•.
Energy Domain Integral, an approach used by (^{Balderrama et al., 2006}) that measures the total change of potential energy (including thermal strains) when the crack advances.

•.
Crack Closure Integral, originally developed for the FEM, which is a stress-based approach to extract the SIF by measuring the force needed to close the crack. (^{Singh et al., 1998}) developed a formulation to be applied with the BEM, obtaining good results.

•.
The Least Squares method was used by (^{Ju, 1998}) to extract the KIII SIF through a least square fit of the stress solution obtained from the FEM.

•.
The Generalized Displacement Correlation method, recently developed by (^{Fu et al., 2012}) for the FEM, uses the displacement solution at crack surfaces for an explicit calculation of the mixed mode SIF.

In this work, a new displacement-based technique is proposed to calculate SIFs for different geometric configurations using the DBEM. Different schemes are proposed to be used with this new technique, and numerical results are compared with those obtained using the J-integral technique in order to compare their accuracy and computing performance.

2.THE BOUNDARY ELEMENT METHOD

The formulation of the BEM is based on Betti's reciprocal theorem where the following integral equation relating displacements u with tractions t for the boundary S is used (^{Becker, 1992}),

where C*ij =* δ*ij/*2 for smooth surfaces, with δ*ij* the Kronecker delta, *T _{ij} * and

*U*are traction and displacement kernels for the displacement integral equation,

_{ij}*p*is the collocation point and

*Q*is a generic boundary point.

The boundary geometry is discretized using quadratic elements (adopted in this work) and then Eq. (^{1}) is written for each node, generating a square system of equations after known boundary conditions are applied.

However, for a cracked body, the crack geometry schematized inFigure 1, defined byand, has the same nodal coordinates if the crack faces are coplanar. This generates an ill-posed problem since Eq. (^{1}) written for thenodes is linearly dependent on theequations. To overcome this issue, (^{Portela et al., 1992}) developed the DBEM for two-dimensional problems.

The DBEM consists on applying Eq. (^{1}) to the non-crack boundary S and one crack face, while the traction integral equation below is applied on the other crack face,

This will result in a well-posed system of equations that can be solved to obtain displacement and tractions fields over the boundary. In Eq. (^{2}), Skij and Dkij are traction and displacement kernels for the traction integral equation.

The kernels for Eq. (^{1}) and Eq. (^{2}) are shown below (^{Aliabadi, 2002}). The different singular behaviors of the integrands require special treatment to obtain meaningful and accurate results.

In this work, the DBEM implementation is done using isoparametric quadratic elements; regular integration is performed using Gauss quadrature, and the Cauchy principal value and the Hadamard finite-part regularization are used to evaluate the kernels with singular integrals.

3 STRESS INTENSITY FACTORS

The validity of the linear elastic fracture mechanics (LEFM) assumption resides in the small scale yielding hypothesis, meaning that plastic strains are only developed, at the crack tip, in a small region compared to the whole geometry, thus they can be neglected since their contribution to the global response is negligible. The stress and displacement fields are given by,

where G is the shear modulus, k = 3 - 4v for plane strain and k = (3 - v)/(1 + v) for plane stress, in accordance with the crack geometry coordinate system shown inFigure 2(^{Gross and Seelig, 2011}).

Now, the fracture mechanics problem with the LEFM approach is reduced to the determination of stress intensity factors. The easiest way to calculate the SIF is by obtaining stress values directly at the crack tip, but this task is unsuitable since numerical results near the crack tip are generally imprecise.

4 COMPUTATION TECHNIQUES FOR SIFS

In order to obtain stress intensity factors for cracked bodies, the classic approach begins by calculating the stress field near the crack tip and extrapolating the results to the crack tip with Eq. (^{7}). This approach leads to numerical inaccuracies due to the singular behavior of the stress field near the crack tip, which is usually underestimated in the elastic solution of the problem by numerical methods including the BEM (^{Cruse, 1996}).

4.1 The J-Integral

The J-integral is a contour integral that measures the strain energy flux across its boundary. Setting the integration contour far from the crack tip, the strain and stress fields can be accurately computed to evaluate the contour integral and obtain the SIF values. The J-Integral for plane problems is calculated as follows:

The relationship between J and the SIF for the LEFM approach is given by

A circumferential integration contour centered at the crack tip is defined to implement the J- integral technique, as shown inFigure 3. The contour is discretized into quadratic elements (^{Bezerra and Medeiros, 2002}) parameterized with an intrinsic variable (z). After some manipulations of Eq. (^{9}), the integrand is rewritten in quadratic form as a function of the displacement gradient as,

The first array in Eq. (^{11}) corresponds to the constitutive matrix that relates stresses to deformations, where λ and μ are Lame´s parameters, while the second array contains the arrangement for the contour outward normal. The contour integration is parameterized by means of the Jacobian jc and the J-integral can be evaluated by the summation of the values at the Gauss points (ng) times their weights wi for each element of the contour ne.

The gradient of the displacement field u (∇u), defined below, is written in vector form and can be decomposed into its symmetric and anti-symmetric parts using the method proposed by (^{Rigby and Aliabadi, 1998}) to obtain KI,

4.2 Displacement-based methods

Displacement extrapolation

From the displacement field solution, the easiest way to obtain the SIF is by using the displacement extrapolation method (DE). This could be accomplished by taking the nearest crack tip opening displacement value obtained immediately after solving the problem and using Eq. (^{13}) to retrieve the SIF directly (^{Aliabadi, 2002}),

This technique is very efficient because the numerical solution at crack nodes is immediately available from the BEM and no internal nodes need to be evaluated.

The displacement field near the crack tip is well behaved and can be retrieved with good accuracy using the BEM even at the crack tip, but the DE results are very sensitive to displacement numerical errors.

Displacement fitting

The displacement field given in Eq. (^{8}) is valid at any point near the crack tip because it only includes the leading termof the power series. However, from the complete crack tip solution, it is known that the next term in the power series is *r*, which must be included in the approximation in order to consider its contribution to the displacement field accompanied with a φ function, as shown in Eq. (^{14}) below

An internal or surface mesh with an arbitrary set of nodes, as shown inFigure 4, can be used to apply this methodology; the displacement field can be retrieved using Eq. (^{1}) at each of the nodes from the BEM solution.

The displacement field can be decomposed to decouple the SIF in Eq. (^{14}), in order to take advantage of its symmetry properties (^{Aliabadi, 2002}). The decomposition can be carried out using Eq. (^{15}). For the given set of twelve internal nodes inFigure 4, the displacement field (modes I and II) is known from the BEM numerical solution. The position of these nodes in the crack tip coordinate system is also known.

Writing Eq. (^{14}) in matrix form, separating unknowns from geometric parameters, the following equation is obtained

The *r* coefficients in Eq. (^{14}) are considered as unknown. The numerical displacement field for each node is equaled to Eq. (14) considering the idea of the DE method. This leads to a system of equations with four unknowns (*K _{I}, K_{II}, c*

_{1},

*c*

_{2}),

The system of equations (^{17}) is solved through the least square method. Renaming the terms in Eq. (^{17}) gives,

The fitting solution is retrieved using the pseudo-inverse approach as the equation system is linear, as follows,

Thus, *K _{I} * and

*K*can be retrieved from a completely arbitrary node distribution, by fitting the numerical solution to the analytic field. This procedure is suitable to the BEM where the internal solution for the displacement field is easily calculated in a post-processing routine using Eq. (

_{II}^{1}), which is more efficient than evaluating the stress integral equation needed in the J-integral calculations.

5.METHODOLOGY

5.1 General

In order to compare the performance and accuracy of the different schemes for calculating the SIF, six examples are solved and compared with their respective reference solutions, which can be found in (^{Tada et al., 1985}) for specimen cases, (^{Shahani and Tabatabaei, 2008}) for the FPB specimen and in (^{API 579-1/ASME FFS-1, 2007}) for cylinder cases. The material properties are set to *E* = 200*GPa* and *v* = 0.3 to solve the linear elasticity problem.

The geometries are generated and meshed using an automatic crack growth algorithm developed in MatLab(. All geometries were modelled from a relation varying from *a/t=0.1* until *a/t=0.6*. *K _{I } *values are calculated at each crack growth step using the following schemes:

•. J-integral: Evaluated using four symmetric elements and a contour radius of one element length to guarantee a straight crack inside the integrating contour.

•. Displacement Fitting Technique (DFT). The following schemes were used to fit the solution considering that this technique could be applied to any arbitrary set of nodes:

○. Surface nodes (S nodes): In the crack tip element, there are three boundary nodes because the element is quadratic as shown inFigure 5. Displacements are directly available from the BEM solution and the DFT can be used to calculate the SIF.

○.
Contour nodes (J nodes): The nodes resulting from the discretization of the circular contour to evaluate the J-integral are taken as follows: a node belongs to the surface of the crevice and the remaining nodes are internal (seeFigure 6). These nodes are evaluated using Eq.(^{1}), obtaining another particular approach to apply the DFT.

○. Internal nodes (M nodes): A symmetric mesh with twelve internal nodes is used for calculating the SIF (Figure 7). Several internal nodes at different angles and radius are used.

Six different geometries are evaluated to determine the performance of the numerical schemes used in this work. The values obtained are compared to the corresponding reference solutions. The comparison is carried out using the following expression:

5.2 Single Edge Notched Tension (SENT)

A SENT geometry is solved with dimensions *t* = 1*m* and *L* = 3*m* as shown inFigure 8with an initial crack length *a* _{o} = 0,1*m* and a crack growth advance of Δ*a*/*t* =0.05. The boundary conditions and loads are also schematized with the BEM contour mesh.

5.3 Three and Four Point Bend (TPB and FPB)

As a second case, a TPB specimen is modelled using the dimensions *t* = 1*m* and *L* = 4*m* (keeping the relation *L*/*t* = 4. On the other hand, the FPB specimen is solved using the dimensions* t* = 1*m, L* = 4*m* and *d* = 1.5*m* An initial crack length *a* _{0} = 0.1*m* is also used and a crack growth advance Δ*a*/*t* = 0.05 is established. The boundary conditions and loading are also schematized in the BEM contour meshes shown inFigure 9.

5.4 Compact Specimen (CS)

The CS geometry is modelled with the required dimensions for the experimental testing given by ASTM E647. However, some simplifications are carried out to model the CS, e.g. the loading pin holes are neglected and the tensile load is directly applied as shown inFigure 10. The geometry thickness is set to *t* = 1*m, a _{n} = 0.*2

*t, L =*1.2

*t,*and the initial crack length

*a*

_{0}=.0.2

*t*

5.5 Thick and Thin Walled Cylinders (CYL1 and CYL50)

Finally, a thick walled cylinder (*R _{i} */

*t*= 1) and a thin walled cylinder (

*R*/

_{i}*t*= 50) are modelled, both with an infinite long radial crack subjected to internal pressure (Figure 11). The boundary conditions are defined at symmetry planes and a pressure loaded crack is located at an angle of 45º measured from the ends. The thickness is set to

*t*= 1

*m*and the initial crack length is

*a*

_{0}=.01

*m*.

6.RESULTS

6.1 Single Edge Notched Tension (SENT)

Results obtained for the SENT geometry are shown inFigure 12. The *K _{I} * values calculated with the J-integral method and the DFT applied with the different schemes used (S-nodes, J-nodes and M-nodes) are compared with the respective reference solution.

First, the general response of the *K _{I } *behavior is well retrieved by both approaches (the J-integral and DFT). On the other hand, comparing the numerical results versus the reference solution by means of Eq. (

^{20}), a numerical error lesser than 2% is obtained for all the proposed schemes except for the S-nodes approach for small cracks. This difference could be attributed to numerical errors associated to the closeness of the crack tip. Nevertheless, the other DFT results are in very good agreement with the reference solution.

Figure 13shows the deformation pattern of the SENT geometry.

6.2 Three Point Bend (TPB)

Figure 14shows the results for the TPB specimen. There, *K _{I} * has a sharper behavior compared with the SENT results, but all the schemes successfully capture this response. In this case, the J-integral method gives more accurate results compared with displacement based ones, although a difference lesser than 4% is found. The S-nodes scheme again produced the highest differences for small crack sizes (

*a/t*≤0.4). These differences are attributed to the model and the boundary closeness, leading to the conclusion that the S-nodes scheme is very sensitive to this error.

The geometry deformation pattern is shown inFigure 15, where the symmetry of the solution and the crack opening mode could be appreciated.

Figure 16shows the results for the FPB specimen. In this example, the crack opening mode II is the most important, and *K _{II} * values are also calculated and compared to the reference solution (

^{Shahani and Tabatabaei, 2008}). In this example, for small cracks, the J-integral method gives more accurate results compared with displacement based ones. For longer cracks, the displacement based methods show comparable precision to the J-Integral.

InFig. 17, a non-symmetric displacement behavior leading to a crack opening mode II is observed.

6.3 Compact Specimen (CS)

Figure 18shows the results obtained for the CS geometry. The J-integral method has a higher difference (almost 7%) to the reference solution; this difference is due to the effect of the initial notch, which is not considered in the reference solution and is mitigated for longer cracks. In this example, the DFT results are more accurate than the J-integral results, with the exception of the S-node scheme. These results justify the simplification made to the original CS geometry. The compact specimen displacement solution is plotted inFigure 19.

6.4.Thick Walled Cylinder (CYL1)

The thick walled cylinder is a more complicated geometry because of its configuration and loading, and a finer mesh is needed to achieve convergence. Results for this geometry are shown inFigure 20, using the reference solution obtained from the API 579 standard for comparison.

The pressure load acting at the inner cylinder face lead to a symmetric load with respect to the crack surface as can be seen inFigure 21. The crack growth mode is K_{I} since the crack grows straight for this load case.

6.5 Thin Walled Cylinder (CYL50)

Finally, a thin walled cylinder is solved. Higher values for *K _{I} * are found in comparison with the previous geometries, and both the J-integral and displacement-based methods give good results with differences of around 1% for long cracks, as can be seen inFigure 22. The deformation pattern of the thin walled cylinder is quite different from the thick one.Figure 23shows a deformed

*R*/

*t*= 50 cylinder, the symmetry is well retrieved and the displacement variation through the thickness is negligible.

6.6 Computing time

The performance of the different schemes in terms of computing time and accuracy is analyzed in the context of the BEM. After solving the elasticity problem, boundary displacements and tractions are known and the solution at any internal point can be retrieved in a post-process routine by means of Eq.(^{1}).

The specific SIF calculations for each scheme are small compared to the calculation time of the internal points, and proportional to the number of evaluating points.Figure 24shows a comparison between the schemes where the number of required evaluating points and internal points for each method are plotted.

The faster scheme corresponds to S nodes because this scheme only uses information from boundary nodes, which are directly available from the BEM solution, followed by J nodes and M nodes, which require less internal points than the J integral. The time spent in the evaluation of internal points depends of many factors, but is independent of the SIF calculation method.

The accuracy of the schemes, however, is inversely proportional to the computation time. The most accurate scheme is the J-integral, closely followed by fitted M nodes and J nodes. S nodes are less accurate than the other schemes as can be seen inFigure 25, where the mean relative error for the different test specimens is shown.

It is noteworthy that the accuracy is proportional to the number of evaluating points. All schemes studied give good accuracy and can be used to estimate K_{I}. Each method has different advantages and is suitable for implementation in the BEM as has been demonstrated. The proposed methodology is also applicable to FEM models.

7 CONCLUSIONS

In this work, different techniques for evaluating stress intensity factors (SIF) in cracked bodies were compared by solving six different geometries using the BEM. The J-integral and three different new schemes based on the Displacement Fitting Technique (DFT) have been developed for the BEM: Surface nodes (S nodes), Contour nodes (J nodes) and Internal nodes (M nodes).

The comparison carried out showed that the J-integral and the M-node method are the most accurate techniques. In terms of efficiency, the computing time of the J-integral method is the highest, due to the calculation of displacement gradients, compared with the displacement-based schemes whose computing time in the BEM context is proportional to the number of internal points. The faster scheme corresponds to the DFT with surface nodes, but it is also the least accurate, concluding that a compromise between accuracy and speed need to be assessed to select the appropriate scheme.

The BEM showed to be an efficient tool for solving fracture problems independently of the scheme used for the calculation of the SIF.