A New Displacement-based Approach to Calculate Stress Intensity Factors With the Boundary Element Method

The analysis of cracked brittle mechanical components considering linear elastic fracture mechanics is usually reduced to the evaluation of stress intensity factors (SIFs). The SIF calculation can be carried out experimentally, theoretically or numerically. Each methodology has its own advantages but the use of numerical methods has become very popular. Several schemes for numerical SIF calculations have been developed, the J-integral method being one of the most widely used because of its energy-like formulation. Additionally, some variations of the J-integral method, such as displacement-based methods, are also becoming popular due to their simplicity. In this work, a simple displacement-based scheme is proposed to calculate SIFs, and its performance is compared with contour integrals. These schemes are all implemented with the Boundary Element Method (BEM) in order to exploit its advantages in crack growth modelling. Some simple examples are solved with the BEM and the calculated SIF values are compared against available solutions, showing good agreement between the different schemes.


Latin American Journal of Solids and
The application of the BEM to fracture mechanics problems was initiated by Cruse through two works presented in 1970 and1971 (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.

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 ‫ܥ‬ = ߜ 2 ⁄ for smooth surfaces, with ߜ the Kronecker delta, ܶ and ܷ are traction and displacement kernels for the displacement integral equation, 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 in Figure 1, defined by S େ ା and S େ ି , has the same nodal coordinates if the crack faces are coplanar.This generates an ill-posed problem since Eq.(1) written for the S େ ା nodes is linearly dependent on the S େ ି equations.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), ܵ and ‫ܦ‬ 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.

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 ‫ܩ‬ is the shear modulus, ݇ = 3 − 4ߥ for plane strain and ݇ = ሺ3 − ߥሻ ሺ1 + ߥሻ ⁄ for plane stress, in accordance with the crack geometry coordinate system shown in Figure 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.

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).

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: Solids and Structures 12 (2015) 1677-1697 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 in Figure 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 ‫ܿܬ‬ and the J-integral can be evaluated by the summation of the values at the Gauss points ሺ݊݃ሻ times their weights ‫ݓ‬ for each element of the contour ݊݁.
The gradient of the displacement field ‫ݑ‬ ‫,)ݑ∇(‬ 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)

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 term ‫ݎ√‬ of the power series.However, from the complete crack tip solution, it is known that the next term in the power series is ‫,ݎ‬ 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 in Figure 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 in Figure 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 ‫ݎ‬ 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 ‫ܭ(‬ ூ , ‫ܭ‬ ூூ , ܿ ଵ , ܿ ଶ ), 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, Latin American Journal of Solids and Structures 12 (2015) 1677-1697 Thus, ‫ܭ‬ ூ and ‫ܭ‬ ூூ 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. ( 1), which is more efficient than evaluating the stress integral equation needed in the J-integral calculations.

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 ‫ܧ‬ = ‫ܽܲܩ002‬ and ߥ = 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.‫ܭ‬ ூ 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: o Surface nodes (S nodes): In the crack tip element, there are three boundary nodes because the element is quadratic as shown in Figure 5. Displacements are directly available from the BEM solution and the DFT can be used to calculate the SIF.o 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 (see Figure 6).These nodes are evaluated using Eq. ( 1), obtaining another particular approach to apply the DFT.o 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:

Single Edge Notched Tension (SENT)
A SENT geometry is solved with dimensions ‫ݐ‬ = 1݉ and ‫ܮ‬ = 3݉ as shown in Figure 8 with an initial crack length ܽ = 0.1݉ and a crack growth advance of ∆ܽ ‫ݐ‬ ⁄ = 0.05.The boundary conditions and loads are also schematized with the BEM contour mesh.

Three and Four Point Bend (TPB and FPB)
As a second case, a TPB specimen is modelled using the dimensions ‫ݐ‬ = 1݉ and ‫ܮ‬ = 4݉ (keeping the relation ‫ܮ‬ ‫ݐ‬ ⁄ = 4).On the other hand, the FPB specimen is solved using the dimensions ‫ݐ‬ = 1݉, ‫ܮ‬ = 6݉ and ݀ = 1.5݉.An initial crack length ܽ = 0.1݉ is also used and a crack growth vance ∆ܽ ‫ݐ‬ ⁄ = 0.05 is established.The boundary conditions and loading are also schematized in the BEM contour meshes shown in Figure 9.

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 in Figure 10.The geometry thickness is set to ‫ݐ‬ = 1݉, ܽ = ‫,ݐ2.0‬‫ܮ‬ = ‫,ݐ2.1‬and the initial crack length ܽ = ‫.ݐ2.0‬

Thick and Thin Walled Cylinders (CYL1 and CYL50)
Finally, a thick walled cylinder (ܴ ‫ݐ‬ ⁄ = 1) and a thin walled cylinder ሺܴ ‫ݐ‬ ⁄ = 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 ‫ݐ‬ = 1݉ and the initial crack length is ܽ = 0.1݉.

Single Edge Notched Tension (SENT)
Results obtained for the SENT geometry are shown in Figure 12.The ‫ܭ‬ ூ values calculated with the J-integral method and the DFT applied with the different schemes used (S-nodes, J-nodes and Mnodes) are compared with the respective reference solution.First, the general response of the ‫ܭ‬ ூ 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 13 shows the deformation pattern of the SENT geometry.

Three Point Bend (TPB)
Figure 14 shows the results for the TPB specimen.There, ‫ܭ‬ ூ 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 in Figure 15, where the symmetry of the solution and the crack opening mode could be appreciated.Figure 16 shows the results for the FPB specimen.In this example, the crack opening mode II is the most important, and ‫ܭ‬ ூூ 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.In Fig. 17, a non-symmetric displacement behavior leading to a crack opening mode II is observed.

Compact Specimen (CS)
Figure 18 shows 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 in Figure 19.

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 in Figure 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 in Figure 21.The crack growth mode is ‫ܭ‬ ூ since the crack grows straight for this load case.

Thin Walled Cylinder (CYL50)
Finally, a thin walled cylinder is solved.Higher values for ‫ܭ‬ ூ 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 in Figure 22.The deformation pattern of the thin walled cylinder is quite different from the thick one.Figure 23 shows a deformed ‫ݐ/ܴ‬ = 50 cylinder, the symmetry is well retrieved and the displacement variation through the thickness is negligible.

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 24 shows 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 less accurate than the other schemes as can be seen in 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 vantages and is suitable for implementation in the BEM as has been demonstrated.The proposed methodology is also applicable to FEM models.Computing load for the different schemes.
ause 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 ends 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 integral, closely followed by fitted M nodes and J nodes.S nodes are less accurate than the other schemes as can be seen in Figure 25, where the mean relative error for It is noteworthy that the accuracy is proportional to the number of evaluating points.All studied give good accuracy and can be used to estimate ‫ܭ‬ ூ .Each method has different advantages and is suitable for implementation in the BEM as has been demonstrated.The proposed e error for the different schemes.

Evaluating points
Internal Nodes M nodes

% Error J nodes
Internal Mesh

Figure 24 :
Figure 24: Computing load for the different schemes

Figure 25 :
Figure 25: Mean relative error for the different schemes