The stress intensity factor assessment in three-dimensional problems by the displacement fitting technique and the dual Boundary Element Method

Corresponding author https://doi.org/10.1590/1679-78256002 Abstract This work presents an extension of the displacement fitting technique for the assessment of stress intensity factors (SIFs) of three-dimensional linear elastic fracture problems using the dual Boundary Element Method. The developed framework accounts for higher-order terms of the asymptotic displacement solution near crack front. The number and location of points surrounding the crack front are properly defined in order to accurately evaluate the SIFs. Three-dimensional benchmarks demonstrate the efficiency of the proposed framework. Moreover, two different fracture criteria illustrate the influence of SIFs values with respect to the crack propagation angle and equivalent factors calculations. The proposed higher-order technique has demonstrated superior performance in comparison with the conventional Three-dimensional fracture


INTRODUCTION
The structural modelling of cracked materials has major importance for the mechanical integrity assessment. The physical discontinuities named as cracks lead to the stress concentrations, which trigger the mechanical collapses. In this context, fatigue and the failure of brittle materials have been handled accurately by the linear elastic fracture mechanics (LEFM). In such case, the stress intensity factors (SIFs) govern the crack propagation and crack growth rates because of the stress singularity at the crack tip. These parameters can be obtained by experimental, theoretical and numerical approaches. However, real-world problems, which involve complex geometry and boundary conditions, are properly handled solely in the context of numerical methods.
Various numerical strategies have been developed for the solution of fracture mechanics problems via boundary element method (BEM) (Portela et al., 1992;Mi and Aliabadi, 1992;Cruse, 1996;Dell´Erba and Aliabadi, 2001;Dirgantara and Aliabadi, 2002;Purbolaksono et al., 2012). Among them, the most robust approach is known as the dual BEM (DBEM), which requires two independent integral equations. Then, the external boundary and one of the crack surfaces are discretised by the collocation of the displacement integral equation whereas the opposite crack surface is discretised by the collocation of the traction integral equation. Portela et al. (1992) and Mi and Aliabadi (1992) proposed this solution technique for two and three-dimensions problems, respectively. In addition to the accuracy and robustness, the DBEM avoids the singularities in the final algebraic system of equations, as previously highlighted by Cruse (1996), in spite of the requirement of a hyper-singular integral kernel. Dell´Erba and Aliabadi (2001) presented a DBEM methodology for the solution of 3D thermo-elasticity problems, in which the J-integral evaluates the SIFs. Dirgantara and Aliabadi (2002) modelled thin crack plates using DBEM. The SIFs were extracted by J-integral and crack surface displacement extrapolation. Rigby and Aliabadi (1998) have decomposed the J-integral into components for mixed modes SIFs assessment. Bezerra and Medeiros (2002) proposed an alternative and efficient numerical scheme for SIF assessment in mode I using J-integral. Ortiz and Cisilino (2005) have developed the J-integral like methodology for 3D interface cracks. It leads to the so-called M-integral, which is a variant of the J-integral and has been used to approximate SIFs for 3D curved loaded cracks (Walters et al., 2005).
It is worth emphasizing that the J-integral requires internal stresses values, which can be achieved straightforwardly at internal points. The BEM handles this accurately because of the non-requirement of the domain mesh. However, this operation usually takes more computational time consuming then the solution of the boundary value problem itself. Thus, alternative J-integral like schemes are largely necessary, which is one contribution of the present study.
The displacement correlation/extrapolation techniques consist of correlating or extrapolating the numerical displacement field at the crack tip with the asymptotic local displacement analytical solution. These solution techniques have been used in its classical form to extract SIFs from the finite element method (FEM) (Chan et al., 1970;Ingraffea and Manu, 1980;Banks-Sills and Sherman, 1986;Banks-Sills et al., 2005;Nejati et al., 2015) and BEM (Cruse and Wilson, 1978;Balderrama et al., 2006;Purbolaksono et al., 2012;Gonzalez et al., 2015;Cordeiro and Leonel, 2018;Cordeiro and Leonel, 2019). Recently, SIFs values have also been obtained from the extended/generalized finite element method (X/GFEM) using such scheme (Minnebo et al., 2010;Gonzalez-Albuixech et al., 2013;Schätzer and Fries, 2016;Gupta et al., 2017;Gupta and Duarte, 2018). Gonzalez et al. (2015) proposed a variation of the displacement correlation technique for two-dimensional problems, which has been named as displacement fitting technique (DFT). This technique incorporates higher-order terms of the displacement asymptotic solutions near the crack front. These high order terms improve the DFT accuracy in comparison with the classical approach, which disregards them. As previously demonstrated by Gonzalez et al. (2015), the DFT is as accurate as the J-integral approaches for SIFs assessment. However, its computational requirements are considerably lower than the J-integral ones.
In the context of 3D fracture problems, the coupling of DFT and DBEM has major potential. From one side, BEM provides accurate modelling of the internal mechanical fields because of the non-requirement of a domain mesh. On the other side, the high order terms enable the proper assessment of SIFs by DFT. However, some studies in the literature report responses instabilities of displacement correlation/fitting techniques (Cordeiro and Leonel, 2019), which can be largely dependent on the extraction points. This subject is particularly handled in the present study and a parametric analysis suggests the better positions for extracting SIFs.
The present study extends the DFT proposed by Gonzalez et al. (2015) for three-dimensional fracture problems. The three-dimensional DBEM formulation provides the solution for internal fields and for the boundary value problem. Three applications demonstrate the accuracy of the proposed solution technique, in which the DBEM results are compared with reference results available in the literature. Moreover, the implementation aspects, as well as numerical stability features, are properly discussed.
Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e310 3/20  Figure 1 illustrates an elastic cracked solid with domain Ω, external boundary Γ , internal crack boundary Γ = Γ + ∪ Γ − in which a free traction condition is assumed at Γ + and Γ − . The discretisation of both crack surfaces leads to the collocation points at a same geometrical position. Thus, singularities appear in the final algebraic system of equations. The DBEM formulation avoids such singularities by setting the equilibrium in terms of a traction boundary integral equation at Γ − and a displacement boundary integral equation at Γ + . The collocation strategy utilised by the DBEM is shown in the figure. The traction (hyper singular) boundary integral equation (BIE) is applied solely for collocation points over Γ − , whereas the displacement (singular) BIE is applied on collocation points at the opposite crack surface Γ + . The external boundary is discretised by the displacement BIE because of its lower singularity level.

Dual boundary element method
Assuming nil body forces, the traction BIE for collocation points − at the crack surface Γ − takes the form where and are the source and field points, with standing for constitutive elastic components and being the traction and displacement Kelvin fundamental solutions (Cruse, 1996). The scalar is a component of the distance vector = − , = ‖ ‖, and η ⁄ is the derivative of with respect to the normal η. The symbol refers to the Kronecker delta, ν is the Poisson's ratio and is the shear modulus. The displacement BIE discretises the collocation points + at the crack surface Γ + integral equation is stated, disregarding the body forces, in the following form Because of the nature singularity of the fundamental solutions, the hyper and strong singular integrals are regularized by the Guiggiani method (Guiggiani et al., 1992). The free terms in (1), (4) and (5) appear during the limit evaluation of the boundary integral equations at Γ − , Γ + and Γ (Chen and Hong, 1999), as usual in BEM. It is worth mentioning that + and − share the same physical coordinates. However, their normal outward vectors are opposite.
The continuous BIE's are the exact representation of the linear elastostatic problem. In the DBEM formulation, boundary geometry and mechanical fields are approximated by basis functions , which are defined over the boundary elements. The functions must be 1 continuous at the collocation points to ensure the existence of the hyper singular integrals. According to Mi and Aliabadi (1992), discontinuous elements at the crack surfaces enable 1 continuity at internal collocation points.
The numerical integration of the BIE's over the discretised boundaries and crack surfaces leads to the DBEM linear system of equations, which is as follows where and are influence matrices, and are vectors containing the displacements and tractions at the collocation points. Notice that there are only half-unknown values of and in this system because of the knowledge of the boundary conditions. To determine them from (6), one rearranges the system by moving columns of and from one side to the other, as usual in BEM. Once all unknows are passed to the left-hand side, one writes where is the vector of unknown's displacements and tractions.
The stress intensity factor assessment in three-dimensional problems by the displacement fitting technique and the dual Boundary Element Method

Displacement fitting technique
The displacement fitting technique (DFT) proposed by Gonzalez et al. (2015) is extended herein for threedimensional problems. Figure 2a illustrates the local Cartesian coordinates ( ′, ′, ′) and spherical coordinates ( , , ) systems both with origin at the near crack front point ′. In terms of the system ( , , = 0), illustrated in Figure 2b, the Sih (1991) asymptotic solution for the crack front displacement field is The angular functions 1 , 1 , 2 , 2 and 3 are well-known from the fracture mechanics literature (Sih, 1991). Herein, the function that involves terms of order ( ) is not explicitly handled and takes, by simplicity, a constant value . Disregarding the terms of order � 3/2 � in (8), the asymptotic solution at a fixed point ( , ) reads Choosing = 1 … arbitrary local points ( , ), one writes in which the vector gathers all the unknown quantities. Since the matrix is non-square, the best solution is given by (Strang, 2006) = 1 that provides , and at the crack front point ′. It is worth mentioning that the least-square solution provided by (12) leads always to a well-conditioned system if one sets 1 = 2 … = = ̅ . The radius ̅ should be of the same magnitude of the characteristic crack dimension (for instance, the largest crack dimension). This procedure is accurate because the extraction distance ̅ is within the range of the fracture mechanics theory. The DBEM enables to obtain the displacement values at the extraction points. For = 7, for instance, all points can be arranged into a uniform circular pattern as presented in Figure 3 (the points are equally spaced around the crack front in the plane = 0).
It is worth emphasizing that DFT is largely suitable for coupling with BEM because the displacements at internal points provided by BEM are straightforwardly utilised into SIFs assessment, without requirements of any additional displacements derivatives, strains or stress evaluations. This feature makes the DFT a computationally efficient approach for SIFs evaluation whenever a parametric analysis with respect to ̅ is performed.

Three-dimensional crack growth criteria
Propagation modelling according the LEFM requires the evaluation of an equivalent stress intensity factor eq = eq ( , , ) and a propagation angle that drives the crack growth after eq has reached the material toughness value c . Such conditions are mathematically expressed in the form The propagation functions and Φ have been modelled by the maximum energy release rate (MERR) criterion and by the Schöllmann's criterion (Schöllmann et al., 2002).

The maximum energy release rate criterion
According to this criterion, the energy release rate governs the crack growth: if the energy release rate associated to the loading pattern reaches the threshold value ≤ then the crack propagates. The direction drived by the propagation angle is the one that maximizes the energy release rate . Briefly, During the crack extension where The stress intensity factor assessment in three-dimensional problems by the displacement fitting technique and the dual Boundary Element Method in which , ′ and ′ defines the stress field along the crack surfaces before the propagation and ‖ ‖ stands for the displacement discontinuity after the propagation.

A closed-form solution for
is not available in the literature. Chang et al. (2006) present an approximated expression, in which the energy release rate is written for = 0. Their solution can be readily adapted for mixed-mode propagation ( ≠ 0) if one assumes that eq ( ) = √2 ′ = √2 cos 2 which brings the following approximated expression: Substituting (18) into (14), taking into account (15)-(17), leads to

The maximum principal stress criterion
The criterion proposed by Schöllmann et al. (2002) is known for its accuracy in mixed mode fracture problems with pronounced mode III. As illustrated in Figure 4, two propagation angles are involved: the kinking angle and the twisting angle ψ . The stress state surrounding the crack tip is sought as sketched in Figure 5. According to this criterion, the propagation direction maximizes the principal stress The stress intensity factor assessment in three-dimensional problems by the displacement fitting technique and the dual Boundary Element Method where , ′ and ′ are the near crack front stress components solutions (Richard et al., 2005). These solutions are largely presented in fracture mechanics theory and have been written as functions of the fundamental stress intensity

NUMERICAL TESTS
The previously mentioned approaches have been implemented in a Fortran code by the authors. Three applications demonstrate the accuracy and performance of the high order DFT in three-dimensional fracture problems. The first application addresses a mode I problem with analytical SIF solution (Murakami, 1987). The second application deals with a mixed mode fracture problem, in which Citarella and Buchholz (2008) provided a numerical SIF reference solution. The last application refers to a penny-shaped crack in an infinity medium, for which analytical solutions are available for the fundamental SIFs ( , , ). The accuracy of the proposed formulation has been demonstrated by a parametric study with respect to the extraction distance . Moreover, the Maximum Energy Release Rate (MERR) and the Maximum Principal Stress (MPS) criteria determine the propagation angles and the equivalent SIF.

Mode I crack problem
The first application handles the pure mode I problem presented in Figure 7. The structure geometry is parallelepiped with an edge notch. Moreover, a unitary tensile load, = 1 (stress unit), opens the crack. The structure dimensions are functions of the crack length = 2.5, as illustrated by the ratios depicted in the figure.  The elastic constants adopted for the material are = 1000 (stress unit) and ν = 0.3. Figure 8 illustrates the boundary element meshes adopted for discretising the structure illustrated in Figure 7.
The meshes are composed of: Mesh 1 (695 collocation points), Mesh 2 (1591 collocation points) and Mesh 3 (3895 collocation points). It is worth mentioning that the crack surfaces were discretised by discontinuous boundary elements, which assure the existence of the hyper singular integrals in (1). The norm of the displacement vector obtained with the three meshes are illustrated in Figure 9. Because of the problem symmetry, the crack surfaces displacements are symmetric and the displacement discontinuities develop in the Mode I direction.  The performance of DFT is illustrated in Figure 10. The SIFs responses were evaluated at the centre crack front point. Moreover, the numerical results have been compared against the analytical predictions described in Murakami (1987). Figure 10 illustrates the SIFs evolution as a function of the dimensionless extraction distance / , in which is the characteristic crack length.
Two versions of the DFT are accounted. The DFT-1 is the conventional approach, in which the near crack front asymptotic displacement expansion (8) involves only the term √ ( , ). In this case, the matrix reduces to r M . The DFT-2 considers both the displacement terms √ ( , ) and ( , ). As shown in Figure 10, the inclusion of the term ( , ) into the DFT system improves considerably the evaluation of the SIF. Moreover, the DFT-2 version becomes less sensitive with respect to the adopted extraction distance , when compared with the DFT-1. These evaluations have been performed for points near the crack front from / = 0.01 to 0.4, where the asymptotic displacement behaviour should be observed. Notice that for this range of the extraction distance, the DFT-1 did not converged for the analytical solution while the higher-order DFT-2 converged for the solution.
From the DFT-2 results, it is possible to observe that the SIFs results converged for the analytical solution at extraction distances closer to the crack front as the mesh refinement was performed. For meshes 1, 2 and 3, the SIF converged at extraction distances / equal to 0.4, 0.35 and 0.20, respectively. This result is expected because as the mesh is refined, the displacement fields singularities near the crack front are better reproduced numerically. Actually, for the ideal analytical situation, the limit of the extraction distance towards zero should reproduces the exact value of the SIF result.

Mixed-mode crack problem
The second application concerns the three-point bending beam with an initial central notch, as illustrated in Figure 11. The following geometrical parameters were adopted: = 260, = 240, = 60 and = 10.
Moreover, the loading is equal to = 2. The initial notch length is 0 = 20 and it is kinked by β = 45° with respect to the beam axis. It leads to mixed mode condition. The beam material is governed by = 70.656 and ν = 0.34 (data in consistent but unspecified units). Numerical reference solutions are given by Citarella and Buchholz (2008), which applied the commercial boundary element software BEASY. The reference solutions account for DBEM formulation and J-integral. The numerical model is composed of 7453 linear boundary elements (10012 collocation points) and the adopted mesh is illustrated in Figure 12. Moreover, the obtained norm of the displacement vector is illustrated in Figure 13, in which the mixed mode condition is highlighted. In this mixed mode problem, it is possible to evaluate non-nil SFIs values. Because of the accuracy demonstrated in the previous application, solely the higher-order DFT-2 is adopted herein. Moreover, extraction distances equal to = 0.1, 0.3, 0.5 are adopted for the displacement assessment, in which these results are presented in terms of normalized stress intensity factor ‖ ‖ = �2 2 /3 � 0 �. Figure 14 presents the SIF evolution along the crack front throughout the solid thickness, in which the responses of the proposed formulation and reference are accounted. One observes that the DFT-2 reproduced properly the SIF responses for the extraction distance approximately equal to = 0.3. The responses were sensible with respect to the extraction distance. Notice that the main differences occur near the external boundaries, as expected. Figure 15 illustrates the propagation angles , ψ and the equivalent stress intensity factor eq obtained for the MERR and the MPS criteria for = 0.3. The results obtained by the proposed formulation are compared against the predictions of Citarella and Buchholz (2008). One observes the good agreement among the responses. As expected, the main differences occur near the external boundaries.

Penny-shaped crack problem
The last application handles a penny-shaped crack surface, with radius , embedded in an infinity elastic medium and subjected to a remote stress , as shown in Figure 16. The analytical SIF solution for this problem was presented by Murakami (1987), which is described by where is the angle between the direction of the remote stress and the plane perpendicular to the crack surface. Two distinct cases were analysed: mode I case ( = 90°) and mixed mode case ( = 45°). The elastic constants adopted for the analysis are = 10 and ν = 0. The remote stress was assumed as unitary = 1. The penny-shaped crack surface radius is = 0.1 (data in consistent but unspecified units).   Four different meshes with linear and quadratic Lagrange elements were adopted to discretise the crack surfaces: Linear M1 (480 collocation points), Linear M2 (992 collocation points), Quadratic M1 (1056 collocation points) and Quadratic M2 (2208 collocation points). Figure 17 illustrates the crack surface discontinuous elements of the DBEM models. The SIF results for all models were analysed with the DFT-2. First the mode I case was studied to obtain the relations between and the extraction distance . Further, the mixed mode case was studied and it was possible to compare the results of all the SIFs ( , , ). Figure 18 presents the deformed crack surfaces for both mode cases with displacements 100 times magnified.

Mode I case
For the mode I ( = 90°) scenario, it was analysed the influence of the extraction distance in the SIFs assessment. Because of the excellent performance in the previous applications, solely the DFT-2 was herein utilised. Figure 19 presents the results of versus / for the different meshes. It is possible to observe the optimum extraction distances for each model. Those values were fixed to evaluate the SIF along the crack front. Figure 20 presents the results of versus the angle (illustrated in Figure 16) obtained by fixing the referred extraction distances obtained for each numerical model. On the other hand, Figure 21 illustrates the propagation angles , ψ and the equivalent stress intensity factor eq obtained for the MERR and the MPS criteria for the model quadratic M2.  These results are compared with the responses obtained from the analytical solution. Because the crack is subjected to pure mode I, both propagation criteria lead to nil propagation angles and the equivalent SIF resulted equal to .

Mixed mode case
Fixing the extraction distances defined in the mode I scenario, it was possible to evaluate the SIFs results for the mixed mode case ( = 45°). The results of ( , , ) versus are presented in Figure 22 for the four DBEM models. Among the different meshes, the most accurate results are obtained by quadratic elements, which provided better geometrical approximation for the circular crack front. The differences observed among the responses of and with linear and with quadratic elements are attributed to the better approximation of quadratic elements in curved geometries. The lack of accuracy for and with linear elements is, therefore, due to the lack of accuracy in the definition of the local crack tip coordinate systems.   Figure 23 illustrates the propagation angles , ψ and the equivalent stress intensity factor eq obtained for the MERR and the MPS criteria for the model quadratic M2. These results are compared with the responses achieved from the analytical solution. Because of the different assumptions assumed by each of the criteria, they lead to different values of propagation angles and equivalent SIF, as expected. Therefore, the accurate crack propagation modelling obviously requires accurate assessment of SIFs. However, the adequate choice of the propagation criteria has major importance on the determination of the crack growth directions and equivalent SIF. Thus, inaccurate crack propagation approaches cannot be attributed solely to fundamental SIFs assessment, but mainly, to the fracture criteria adopted.

CONCLUSION
This study presented the extension of the Displacement Fitting Technique (DFT) for SIFs assessment in threedimensional LEFM problems. The higher-order terms of the asymptotic near crack front solutions were accounted into the fitting scheme and lead to accurate SIFs results in comparison with the classical first order correlation techniques, which is largely applied in the literature. This improvement is mainly emphasized in the first application, where both the conventional DFT (DFT-1) and the higher-order DFT (DFT-2) were adopted. For extraction distances ranging from / = 0.01 to 0.4 the DFT-2 converged in all analysis for the analytical SIF solution. On the other hand, the DFT-1 did not converged. Thus, the improvement that the inclusion of the higher-order terms brings to the SIF solution is evident. It is worth mentioning that the DFT is computationally efficient because displacement derivatives and stresses are not required for the SIF evaluation. Moreover, its coupling with DBEM is adequate because DBEM provides accurate responses for the internal mechanical fields.
As expected, the numerical SIFs responses obtained with the DFT are still sensible with respect to the extraction distance , when the evaluation points are equally distributed around a circular path of radius around the crack front. Nevertheless, the present study suggested distance ranges, in which accurate SIFs values are obtained by the technique. Two propagation criteria were considered to evaluate the propagation angles and the equivalent SIF. Then, the accuracy of the fundamental SIFs was analysed in terms of global variables. Moreover, the last application demonstrated that the SIFs accuracy has major importance on the crack propagation modelling. However, the proper choice of the propagation criterion should be of major concern by the analysts, because huge behaviour changes are observed according to the chosen criterion.
Finally, the amount and position of the internal near crack front points for displacement evaluation can be properly performed by DBEM, because of its non-requirement of a domain mesh. Moreover, such modifications do not lead to further mesh modifications, which is a huge advantage over domain mesh methods such as FEM. Such improvements could be explored in future studies.