Machine-Precise Evaluation of Stress Intensity Factors with the Consistent Boundary Element Method

Fracture mechanics phenomena have been systematically investigated since the seminal studies by Inglis [1] in the year 1913, who showed the importance of analyzing the stress concentration around straight cracks in plane plates. Griffith [2] introduced in 1920 thermodynamics con- cepts to formulate a fracture mechanics theory based on energy balance. Later on, Irwin extended Griffith’s ideas to metals and proposed the concept of energy release rate [3] and, on the basis of the developments by Westergaard [4], also the stress intensity factor (SIF) as a measure to be definitely considered in the fracture mechanics theory [5]. The concept of a crack opening displacement (COD) was proposed by Dugdale in 1960 to account for yield- ing. Shortly thereafter, Rice [7] showed that the energy release rate may be obtained by means of a line integral, for two-dimensional (2D) problems, which became known as the J-integral. These concepts have been since then generalized for mixed-mode phenomena at a crack’s tip, whether or not in the elastic range (e.g. Jiangbo et al. [8]).

The application of these concepts to a practical problem depends on the adequate evaluation of the stress state around a crack tip. The first boundary element development in the field was made by Cruse and Van Buren in the year 1971 [9]. However, owing to space restrictions, we cannot carry out a literature review on the subject.

The second author and previous collaborators have already, in part, successfully arrived at an attempt on the basis of the hybrid boundary element method to adequately represent the stress field around a crack tip in terms of Williams’ series as well as of generalizations of Westergaard’s stress functions [10–12].

More recently, a consistent formulation – actually both conceptually and numerically long due to correction – of the conventional collocation boundary element method (CBEM) was proposed [13,14], according to which a convergence theorem could be envisaged and the evaluation of results turns out to be achievable within machine precision for internal points arbitrarily close to an elastic body’s boundary, regardless of the topology issues. The numer- ical implementation considers arbitrarily high-order boundary elements for 2D bodies of any shape and topology as well as flat boundary elements for 3D bodies [15]. The application of this consistent formulation to 2D fracture mechanics problems was the subject of the second author’s M.Sc. thesis, which showed that a generically curved crack may be numerically simulated exactly as if in a real-life experimentation, with the faces distinctly modeled for arbitrarily small openings as long as the continuum mechanics is deemed applicable.

The following outline extends the authors’first attempt to report some academic applications of the CBEM to linear fracture mechanics [17], keeping for conciseness much of this first paper’s introductory part as well as the same numerical examples – only better elaborated. Some more examples adapted from [16] are also presented.

As applied to an elastostatic problem, boundary nodal displacements d and traction attributes

t are interrelated by means of the single- and double-layer potential matrices G and H as

In this equation, it is also considered that body forces may be at least approximately expressed in terms of equivalent nodal displacements $\boldsymbol{d}^p$ and boundary traction parameters $\boldsymbol{t}^p$, which is not always a simple task $[13,14,16]$. This equation turns out to be an application of Somigliana's identity that converts boundary data into domain displacements. For the reallife case of inaccurate boundary data, in general, as the result of approximated values of $\boldsymbol{d} \equiv d_f$ and $\boldsymbol{t} \equiv d_{\ell}$ for a given problem as well as of the piecewise boundary interpolation, an error term $\varepsilon$ has to be considered in this equation [13,14], from which the concept of an admissible matrix $\boldsymbol{G}_{\text {ad }}$ is derived - a matrix that does not transform traction forces out of balance. (Owing to space restrictions, this concept is not further explored herein.) The matrices $\boldsymbol{G}$ and $\boldsymbol{H}$ are expressed in terms of the boundary integrals

where $u_{i s}^*\left(\boldsymbol{x}-\boldsymbol{x}_s\right)$ and $\sigma_{j i s}^*\left(\boldsymbol{x}-\boldsymbol{x}_s\right)$ are the displacement and stress (Kelvin's) fundamental solutions of the elastic problem, respectively - which have global support - and $\Gamma(\boldsymbol{x})$ is the integration boundary. In the above equation and in the following, repeated indices mean summation. The vector $\boldsymbol{x} \equiv(x, y, z)$ stands for the Cartesian coordinates of a given point, in the case a field point, and $n_j(\boldsymbol{x}(\xi, \eta))$ are the Cartesian components of the unity outward vector from $\vec{n}(\boldsymbol{x}(\xi, \eta))$ to $\Gamma(\boldsymbol{x}(\xi, \eta))$, in terms of parametric variables $(\xi, \eta)$, for a general 3D problem. The subscript $s$ refers to a given source node (at which the unit point force of the singular fundamental solution is applied), and the subscripts $f$ (which stands for field) and $\ell$ (also a field reference) indicate respectively to which node or surface point the displacement-interpolation function $u_{i f}(\boldsymbol{x}(\xi, \eta))$ or the traction-interpolation function $t_{i \ell}(\boldsymbol{x}(\xi, \eta))$ - both with local support - is referred. $u_{i f}(\boldsymbol{x}(\xi, \eta))$ comes from the piecewise interpolation of displacements $u_i(\boldsymbol{x}(\xi, \eta))$ along the boundary, $u_i(\boldsymbol{x})=u_{i f}(\boldsymbol{x}(\xi, \eta)) d_f$, where $d_f$ is the nodal displacement. In a practical finite element or boundary element implementation, $u_{\text {if }}(\boldsymbol{x}(\xi, \eta))$ is actually represented by polynomial shape functions $N_f(\xi, \eta)$ :

In the expression of $\boldsymbol{H}$, the Jacobian used in the definition of $n_j(\boldsymbol{x}(\xi, \eta))$ cancels out with the Jacobian of $\mathrm{d} \Gamma(\boldsymbol{x})=|J(\xi, \eta)| \mathrm{d} \xi \mathrm{d} \eta$. For the single-layer potential matrix $\boldsymbol{G}$, it is proposed that the usual (as found in the literature) interpolation polynomials $t_{i \ell}$ of traction forces in eqn (2) be replaced with

where $|J|_{\text {(at } \ell)}$ is the value of the Jacobian at the point characterized by the subscript $\ell[13,14]$. Nothing changes formally in the development of BEM for curved boundary segments (and, of course, nothing changes numerically for the trivial cases of straight or flat boundary segments), except that the evaluation of $\boldsymbol{G}$ becomes much easier and actually more consistent as compared to the proposed implementations given in the technical literature. In fact, $|J|$ cancels out in the product $t_{i \ell} \mathrm{~d} \Gamma$ in eqn (2) for $t_{i \ell}$ defined as suggested, and the integrand of $\boldsymbol{G}$ becomes a polynomial that multiplies the assumed kernel $u_{i s}^*$.

As proposed above, the integrands in eqn (2) are given as products of the kernels with polynomial terms, regardless of the element order for either 2D or 3D problems. For 2D problems, it has been shown that boundary elements of any shape may be dealt with in a unified procedure that enables machine precision results for the evaluation not only of $\boldsymbol{G}$ and $\boldsymbol{H}$ but also for displacement and stress results at internal points - in the latter case, allowing for eventual quasi hypersingularities - no matter how close source and field points may be.

The outline of the numerical evaluation procedures is out of scope in the present developments, but it may be worth illustrating in Fig. 1 the kinds of singularity that must be taken into account. For the curved element of the figure, given in terms of the parametric variable $\xi \in[0,1]$, an actual singularity occurs for the source point at A , thus for $0 \leq \xi \leq 1$. For the source point at B , represented by $\xi>1$ but still close to the integration interval, the integrals of eqn (2) must be dealt with for the special case of a real quasi singularity. On the other hand, if the source point is at C , its complex location $\xi=a \pm i b$ must be found for the correct consideration of the corresponding procedure [14, 18].

specified, we consider a material with transversal elasticity modulus G = 80,000 units of stress and Poisson’s ratio v = 0.25 for plane stress state. All singularity and quasi-singularity cases are evaluated within machine precision, but the use of ng = 4 Gaussian points for the regular integrals leads to final precision of results that cannot be bet- ter than 10-7. In fact, several numerical assessments conducted using four, six and eight Gaussian points combined with 20, 30 and 40 precision digits implemented in a Maple™ code (Maple-15, Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario) showed that all results of interest could be obtained within a precision of 10-5 for ng = 4points and 30 digits [16], which will be used here for all regular integrals along generally curved quadratic ele- ments. Mode I SIFs are evaluated in terms of J-integrals along circles centered at the crack tip, subdivided into 10 sectors and also using ng = 4 Gaussian points per sector, as well as in terms of crack tip opening displacement (CTOD).

The simplest and best representative case of fracture mechanics is the one of a straight crack in the open domain for a constant stress state applied transversally to the crack at infinity, with stress disturbances to be assessed at the crack’s tips and mode I SIF to be evaluated. We represent this crack with length 2a = 2 units and already as an elliptic opening with 2b = 2 × 10−3 units, which, only in the limit case of b → 0, would exactly correspond to the problem originally proposed and solved by Westergaard [4].

In eqn (1), this case corresponds to a far-field stress state σij(p) evaluated as traction parame-

ters tp and nodal displacements dp = 0 at the crack borders. Moreover, the occurrence of the crack implies that t = 0, so that eqn (1) directly leads to the solution of the crack borders’ nodal displacements d as

since H is non-singular for an open domain problem, although we might expect some ill-con- ditioning due to the proximity of the crack faces (this is assessed in the next example). After evaluation of d, results at internal points are obtained, with final stress values given as the superposition with the applied far-field particular solution σij(p) .

A convergence study is shown on the leftin Fig.2forsimulationswith4–60quadraticele-ments along the crack faces (ellipse semiaxisb= 10-3), with the SIF evaluated in terms ofthe J-integral and compared with the analytical result KI No advantage istaken ofthe problem’s double symmetry.) The node spacing along the crack varies geometri-cally fromthe cracktiptothe center accordingtothe indicatedratios.We seethat increasingly agglutinating the nodes close to the crack tip (ratio = 1.5) does not necessarily lead to better results, which are for ratio 1.35, with errors tending to 10-4. This is mainly related to the geometric approximation of the assumed elliptic initial crack shape in terms of successive quadratic segments, which provokes an undesired angulation between adjacent elements instead of a smooth surface.

In order to assess how the size of the ellipse’s minor axis affects results, the right graphic of Fig. 2 shows relative errors in the evaluation of SIFs for b = 10-3, 10-4 and 10-5 units for different ratios of the geometric increase of node spacing from the tip to the middle of the crack. There is a tendency to arrive at the SIF for the fracture mechanics theory, although the angulation between adjacent elements causes too much disturbance: the best result is the error of about 10-6 for the opening b = 10-5 with ratio 1.1 for the internode space increasing.

Theangulationbetweenadjacentquadraticelementsalongasmoothboundarymaybe measured in termsof how much thediscontinuouspartof thematrixHforanodebetween elements differs from 0.5. Sinceallsingular evaluationsarecarriedoutwithinmachinepre- cision, this deviation from 0.5isareliableaccountof theintroducedgeometryerrors.Suchangulationoccursfornodes3and119of theinsertinFig.3andisgivenaserrorsof themeasuredvalues (1- θ/2π) relativeto0.5fortwodifferentincreasingratiosofthedis-tancebetweenadjacentnodesinthecaseofthehalfopening b = 10-5 and atotalof60 quadratic elements in the crack discretization, which corresponds to 31 nodes from the tip toface’s middle. Node1 is at the crack tip and the results shown are absolute values thatwouldtend to zero as b→ 0.The geometry errors are not large, in general, and smaller for nodes nottoo differently spaced, but larger when closer to the crack tip, where spurious stress gradients should better be avoided. The conclusion is that results in the values of the SIF cannot be expected to improve unless this angulation issue is adequately solved.

Figure 4 presents on the left normal stress results at 40 points spaced at geometrically increasing distances r / a along a straight line extending from the crack tip, as obtained for the discretization with b = 10-3 and 60 elements with node spacing ratio 1.35, and are com- pared with the reference values according to Westergaard’s developments. The closest point is at a relative distance 10-3. The 1/ r singularity tendency is well followed. The comparison with the target results is better displayed in terms of the error norm on the right of the figure, where both vertical and horizontal axes are in logarithmic scale.

As assessed by Amaral Neto [16] for this case of an immersed straight crack as well as for some other crack configurations, the SIF results are actually dependent on the assumed Pois- son’s ratio, a fact generally ignored in the literature, maybe because the numerical differences are smaller than the built-in errors of the hitherto implemented simulations.

The normalized modes I and II SIFs, KI and KII , of the linear fracture mechanics problem may be evaluated in terms of the CTODs Δu丄 and Δu=, that is, the relative opening and slip displacements of the two opposite face points that come as close to the crack tip as possible, such as for nodes 1 and 120 in the insert of Fig. 3,

for plane strain problems, where G is the shear modulus and v is the Poisson’s ratio. In the present CBEM, we might evaluate displacements in the domain (just the Somigliana’s iden- tity is accurately evaluated) as close to node 1 in the insert of Fig. 3 and as close to the crack border as possible. This has not been done yet, but we have rather considered for the evalua- tion of KI in eqn (7), the opposite nodal displacements along the crack faces – directly obtained in eqn (6) – such as in the sequence for nodes 2 through 31 in the two graphics of Fig. 5 for crack openings 2 × 10-3 and 2 × 10-5 using either four or eight quadrature points per quadratic element, both cases with the ratio 1.35 of the geometric increase of node spac- ing (which justifies using the logarithmic scale for the indicated node distances in the horizontal axes). This example shows that all evaluations are sufficiently accurate with just four quadrature points per element. Moreover, the procedure of using eqn (7) is by far simpler than the tedious evaluation of the J-integral around the crack tip. On the other hand, the spu- rious angulosities already reported and assessed by using Fig. 3 may affect the results’ accuracy more than when compared with the ones using the J-integral, as in Fig. 2. For a more elaborated crack simulation with a total of 120 quadratic elements, the results shown in Fig. 6, which are to be directly compared with the ones of Fig. 5, follow quite smoothly the tendency of eqn (7), but also present some disturbances close to crack tip owing to the reported lack of smoothness. Using eqn (7) is by far more expedite than using the J-integral, but leads to results that seem to be at least one order of magnitude less accurate. Observe that since the original resource to the CTOD applies to a crack with initial zero opening, eqn (7) must be accrued of the inherent error O (∆θ)2 of considering that our proposed degenerated crack ellipse already has some opening.

The degenerated elliptic crack between two holes in the open domain shown on the left in Fig. 7 with the opening $2 b=2 \times 10^{-3}$ is analyzed as in the previous section. For all simulations, each circle is discretized with a fixed number of 40 quadratic elements, with internode distance increasing from the horizontal symmetry axes with the ratio 1.2 , while varying from 4 to 60 quadratic elements the discretization of the elliptic crack and using the internode increasing ratios shown on the right of the figure. The mode I SIF reference value for the proposed geometric configuration (but straight crack) and loading is $K_I=2.0402701$ [19]. The SIF relative error obtained in our evaluations by means of the J -integral tends to be about $10^{-3}$, as shown in Fig. 7, a result to be interpreted with all the caveats of the previous section.

The rectangular plate with an inclined edge crack of Fig. 8 (drawing not in scale) was exam- ined by Amaral Neto [16] for a longitudinal constant load and edge openings 2 × 10-3, as indicated, as well as 2 × 10-4 and 2 × 10-5. The discretization of the four external edges con- sisted of a fixed number of 80 uniformly distributed quadratic elements. In the numerical analyses, the numbers of quadratic elements along each crack face varied from 1 to 15, and results for several internode space increasing ratios from the crack tip to the edge were inves- tigated, as described for the previous example. Overall results tended to be better when such ratios were not too large, but just about 1.2–1.3. The obtained results were quite insensitive to discretizations with more than four elements along each crack face.

Since this is the case of a finite elastic body, a very objective assessment of the imple- mented numerical solution is worth being undertaken. In fact, as already used by the first author in several numerical investigations, Amaral Neto [16] submitted the plate of Fig. 8 to a family of displacement fields given by polynomials of degrees one, two, three, four and five, a total of 20 fields (four for each degree) that are non-singular fundamental solutions of the proposed elastic problem. For each polynomial, boundary displacements and traction forces can be evaluated and inserted into eqn (1) in order to check its accuracy as stated in this equa- tion as well as in terms of the system solution indicated in eqn (6) – after duly compensating for rigid body displacements [16]. In fact, while sheer discretization accuracy can be assessed by means of eqn (1), the matrix inversion in eqn (6) makes it possible to measure the eventual onset of ill-conditioning.

According to Dumont’s convergence theorem for the CBEM[14], the numerical results areexpected to beaccurate within machine precision forlineardisplacementfields.(Itwouldbethecaseforquadraticdisplacementfieldsaswell,if the bordersweregiven bystraightseg-ments,whichisstrictlyspeakingnotthecasehere,asthecrackfacenodesarenotevenlyspaced inside a quadraticboundary element and the Jacobian ofthe coordinates transformation is, therefore, not constant.) Since the regular parts of the integrals are evaluated numerically, the expected global results cannot be more accurate than the quadratures allow for. Figure 9 shows the error results for eqn (1), Gtpol - Hdpol , as well as for eqn (6), dpol - dnum , in the cases of crack openings 2 × 10-3, on the left, and 2 × 10-5, using four, six and eight quadrature points. The horizontal axes show the average results for the four polynomials of each degree used in the analyses. The results for the first-degree polynomials – for which analytical results are expected – are multiplied by some constants in order to keep the graphics within a reasonable interval. It is expected that errors increase as the polynomial degrees increase. In fact, as the quintic polynomials, for instance, give rise to overall high stress gradients, the displayed results are a threshold of the accuracy we should expect in the fracture mechanics assessments. It is worth observing some ill-conditioning related to eqn (6) – which increases as the crack opening gets smaller – but way below the approximation errors, which allows us to attest the robustness and reliability of the proposed numerical simulations.

A configuration similar to the one of Fig. 9 was analyzed by Namakian et al. [20] using two quite similar meshless implementations (results here are referred to as (a) and (b)), with nor- malized modes I and II SIFs KI and KII compared with another meshless model by Zhuang et al. [21] as well as with the key reference results by Murakami [22]. The plate considered by these authors has dimensions 444.5 × 177.8mm2, with the 45。slanted crack of (very large) length 177.8/ 2 mm starting 177.8mm from one border, so that two pieces of preserved 177.8 × 177.8mm2 squares are kept before and after the crack formation. The plate is submit- ted to uniform normal stress along its short edges for plane strain state and Poisson’s ratio v = 0.25. The (KI, KII ) results found by these authors were apud [20] respectively (1.830, 0.814) 20](a), (1.833, 0.819) [20](a), (1.778, 0.799) [21] and (1.856, 0.834) [22].

We ran two analyses for edge openings 2 × 10-3 and 2 × 10-5 mm. The plate was simulated with 8 + 20 + 8 + 20 = 56 equally spaced quadratic elements along the borders and 22 quad- ratic elements along each of the crack faces, for internode distance increasing from crack tip by a factor 1.15: the distance from the crack tip to the first node was 0.04034 mm (≈ 0.00032 relative to crack length).

Figure 10 shows mode I (on the left) and mode II SIF results evaluated according to eqn (7) for the two simulations, as compared with the values given by Murakami [22]. Results seem to be quite stable for nodes around a relative distance r / a ≈ 0.011 from the crack tip, node 141, which corresponds to the opposite nodes 128 and 154, 2.204 × 10-5 and 2.204 × 10-7 mm apart for the edge openings 2 × 10-3 and 2 × 10-5 mm, respectively. The corresponding (KI, KII ) values are (1.845, 0.824) and (1.837, 0.824), which are confronted in Table 1 with the previously reported results. Our values lie between the results obtained by Namakian et al. [20] and Murakami [22] and are most probably the more reliable ones.

It has been shown that a code implementation of the consistent CBEM can be used to analyze crack configurations of the 2D fracture mechanics as they are expected to exist in real life [16,17]. In this code, a unified integration scheme deals with singularities and quasi singular- ities of an arbitrary degree of severity, as source and field points can be infinitely close to each other, the only limitation being the machine’s capability of representing numbers so as to avoid uncontrollable round-off errors [14]. One might think of working with a code that com- bines Kelvin’s singular and hypersingular fundamental solutions to model a crack’s face and its opposite, respectively, a classical usage to represent cracks of zero opening. The present code actually already allows for such simulation, as the evaluation of stress results at internal points deals with quasi hypersingularities within machine precision. However, modeling real crack configurations opens up the possibility of implementing a procedure to investigate the crack propagation phenomenon from its onset and during successive steps of the crack becoming wider and larger – and eventually kinking and bifurcating. The combination with cohesive elements is also quite natural in the present code implementation.

An issue investigated in the example given in Section 3.1 is the spurious angulations between adjacent elements in the attempt to represent an elliptic geometry as a succession of quadratic elements – aggravated at the proximity of the crack tip as some stress disturbance is unwillingly introduced. (An isogeometric implementation is not a solution, as the consist- ent BEM must rely on an isoparametric formulation [14].) The main concern is, however, how to assume the shape of an unloaded crack in a real sample other than making local meas- urements and just trying to have representative models. Westergaard [4] developed an analytical solution for a straight crack with initial zero opening, for which a far-field constant stress state is shown to cause an elliptic opening. This is the reason to suppose in Section 3.1 that the unloaded straight crack is elliptic with the minor axis tending to zero. In the config- uration for a = 10-3 , the distance between opposite nodes closest to the crack tip is as small as 1.8 × 10-5 for the model with 60 boundary elements. The assumption of a straight edge crack of Fig. 8 seems quite reasonable, but the linear narrowing from the edge to the tip is also in need of justification. In such a case, there are no spurious angulations and the crack tip is very narrow indeed: for the configuration with the indicated opening 2 × 10-3 and a total of 30 boundary elements, the distance between the opposite nodes closest to the crack tip is 8.6 × 10-8, which is already going beyond the continuous mechanics assumptions, but is still doable mathematically.

A possibility that is already under investigation is the modeling of the crack as arbitrarily close, parallel (or not) interfaces interconnected by curved transition elements that preserve surface smoothness wherever required (this is done by just adjusting the loci of the internal element nodes). The simplest case of the straight crack of Section 3.1 would be modeled as in Fig. 11, for quadratic, cubic or quartic boundary elements.

The model proposed in this figure seems to enable a friendlier manipulation of the numer- ical simulations as a crack propagates and eventually kinks or bifurcates. This is the subject of a research work in progress.