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

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 displacement fitting technique.

Keywords:
Displacement Fitting Technique; Dual Boundary Element Method; Three-dimensional fracture problems

Graphical Abstract

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

Figure 1:
Collocation of the boundary integral equations.

2 FUNDAMENTALS

2.1 Dual boundary element method

Figure 1 illustrates an elastic cracked solid with domain $$, external boundary ${}_b$, internal crack boundary ${}_c={}_c^{+}\cup {}_c^{-}$ in which a free traction condition is assumed at ${}_c^{+}$ and ${}_c^{-}$. 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 ${}_c^{-}$ and a displacement boundary integral equation at ${}_c^{+}$. 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 ${}_c^{-}$, whereas the displacement (singular) BIE is applied on collocation points at the opposite crack surface ${}_c^{+}$. The external boundary is discretised by the displacement BIE because of its lower singularity level.

Assuming nil body forces, the traction BIE for collocation points ${\mathbf{y}}^{-}$ at the crack surface ${}_c^{-}$ takes the form

$\frac{1}{2}{t}_i\left({\mathbf{y}}^{-}\right)-\frac{1}{2}{t}_i\left({\mathbf{y}}^{+}\right)+{}_j\left({\mathbf{y}}^{-}\right){\int }_{}^{}{S}_{kij}\left({\mathbf{y}}^{-},\mathbf{x}\right)u_k\left(\mathbf{x}\right)d={}_j\left({\mathbf{y}}^{-}\right){\int }_{}^{}{D}_{kij}\left({\mathbf{y}}^{-},\mathbf{x}\right)t_k\left(\mathbf{x}\right)d$ (1)

where $\mathbf{y}$ and $\mathbf{x}$ are the source and field points, $u_k\left(.\right)$ and $t_k\left(.\right)$ are the displacements and tractions components, ${}_i$ is the direction cosine of the normal outward vector $$, ${\int }_{}^{}\left(\right)d$ denotes the Hadamard Finite Part (hyper singular) integral which is only defined for collocation points at smooth boundaries, ${\int }_{}^{}\left(\right)d$ indicates a Cauchy Principal Value (strong singular) integral,

$S_{kij}\left(\mathbf{y},\mathbf{x}\right)=C_{ikpq}\frac{\partial T_{pj}(\mathbf{y},\mathbf{x})}{\partial y_q}{D}_{kij}\left(\mathbf{y},\mathbf{x}\right)=C_{ikpq}\frac{\partial U_{pj}(\mathbf{y},\mathbf{x})}{\partial y_q}$ (2)

with $C_{ikpq}$ standing for constitutive elastic components and

$T_{ij}\left(\mathbf{y},\mathbf{x}\right)=\frac{1}{8\pi \left(1-\right)}\frac{1}{r^2}\left\{\frac{\partial r}{\partial }\left[\left(1-2\right){\delta }_{ij}+3r_i{r}_j\right]-\left(1-2\right)\left(r_i{}_j+r_j{}_i\right)\right\}$ $U_{ij}\left(\mathbf{y},\mathbf{x}\right)=\frac{1}{16\pi \mu \left(1-\right)}\left[\left(3-4\right){\delta }_{ij}+r_i{r}_j\right]$ (3)

being the traction and displacement Kelvin fundamental solutions (Cruse, 1996). The scalar $r_i$ is a component of the distance vector $\mathbf{r}=\mathbf{y}-\mathbf{x}$, $r=‖\mathbf{r}‖$, and $\partial r/\partial$ is the derivative of $\mathbf{r}$ with respect to the normal $$. The symbol ${\delta }_{ij}$ refers to the Kronecker delta, $$ is the Poisson’s ratio and $\mu$ is the shear modulus.

Figure 2:
(a) Local crack front coordinates; (b) plane crack front coordinates for $\phi =0$.

The displacement BIE discretises the collocation points ${\mathbf{y}}^{+}$ at the crack surface ${}_c^{+}$ integral equation is stated, disregarding the body forces, in the following form

$c_{ij}\left({\mathbf{y}}^{+}\right)u_j\left({\mathbf{y}}^{+}\right)+c_{ij}\left({\mathbf{y}}^{-}\right)u_j\left({\mathbf{y}}^{-}\right)+{\int }_{}^{}{T}_{ij}\left({\mathbf{y}}^{+},\mathbf{x}\right)u_j\left(\mathbf{x}\right)d={\int }_{}^{}{U}_{ij}\left({\mathbf{y}}^{+},\mathbf{x}\right)t_j\left(\mathbf{x}\right)d$ (4)

Its counterpart at ${}_b$ is

$c_{ij}\left(\mathbf{y}\right)u_j\left(\mathbf{y}\right)+{\int }_{}^{}{T}_{ij}\left(\mathbf{y},\mathbf{x}\right)u_j\left(\mathbf{x}\right)d={\int }_{}^{}{U}_{ij}\left(\mathbf{y},\mathbf{x}\right)t_j\left(\mathbf{x}\right)d$ (5)

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 ${}_c^{-}$, ${}_c^{+}$ and ${}_b$ (Chen and Hong, 1999), as usual in BEM. It is worth mentioning that ${\mathbf{y}}^{+}$ and ${\mathbf{y}}^{-}$ 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 $N^{\alpha }$, which are defined over the boundary elements. The functions $N^{\alpha }$ must be $C^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 $C^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

$\mathbf{H}\mathbf{u}=\mathbf{G}\mathbf{t}$ (6)

where $\mathbf{H}$ and $\mathbf{G}$ are influence matrices, $\mathbf{u}$ and $\mathbf{t}$ are vectors containing the displacements and tractions at the collocation points. Notice that there are only half-unknown values of $\mathbf{u}$ and $\mathbf{t}$ in this system because of the knowledge of the boundary conditions. To determine them from (6), one rearranges the system by moving columns of $\mathbf{H}$ and $\mathbf{G}$ from one side to the other, as usual in BEM. Once all unknows are passed to the left-hand side, one writes

$\mathbf{A}\mathbf{\alpha }=\mathbf{b}$ (7)

where $\mathbf{\alpha }$ is the vector of unknown’s displacements and tractions.

2.2 Displacement fitting technique

The displacement fitting technique (DFT) proposed by Gonzalez et al. (2015) is extended herein for three-dimensional problems. Figure 2a illustrates the local Cartesian coordinates $(x\text{'},y\text{'},z\text{'})$ and spherical coordinates $(r,\theta ,\phi )$ systems both with origin at the near crack front point $O\text{'}$. In terms of the system $(r,\theta ,\phi =0)$, illustrated in Figure 2b, the Sih (1991) asymptotic solution for the crack front displacement field is

$\mathbf{u}\left(r,\theta \right)={\mathbf{u}}_{\sqrt{r}}\left(r,\theta \right)+{\mathbf{u}}_r\left(r,\theta \right)+O\left(r^{3/2}\right)={\left[\begin{array}{ccc}{u}_{x^{\mathrm{\text{'}}}}\left(r,\theta \right)& u_{y^{\mathrm{\text{'}}}}\left(r,\theta \right)& u_{z^{\mathrm{\text{'}}}}\left(r,\theta \right)\end{array}\right]}^T$ (8)

where

${\mathbf{u}}_{\sqrt{r}}\left(r,\theta \right)=\frac{1}{\mu \sqrt{8\pi }}\sqrt{r}\left[\begin{array}{ccc}{f}_{1I}\left(\theta \right)& f_{1II}\left(\theta \right)& 0\\ f_{2I}\left(\theta \right)& f_{2II}\left(\theta \right)& 0\\ 0& 0& f_{3III}\left(\theta \right)\end{array}\right]\left[\begin{array}{c}{K}_I\\ K_{II}\\ K_{III}\end{array}\right]=\frac{1}{\mu \sqrt{8\pi }}\sqrt{r}\mathbf{F}\left(\theta \right)\mathbf{K}$ ${\mathbf{u}}_r\left(r,\theta \right)=\frac{1}{\mu \sqrt{8\pi }}r\left[\begin{array}{c}{g}_1\left(\theta \right)\\ g_2\left(\theta \right)\\ g_3\left(\theta \right)\end{array}\right]=\frac{1}{\mu \sqrt{8\pi }}r\mathbf{g}\left(\theta \right)$ (9)

The angular functions $f_{1I}$, $f_{1II}$, $f_{2I}$, $f_{2II}$ and $f_{3III}$ are well-known from the fracture mechanics literature (Sih, 1991). Herein, the function $\mathbf{g}$ that involves terms of order $O\left(r\right)$ is not explicitly handled and takes, by simplicity, a constant value $\mathbf{c}$. Disregarding the terms of order $O\left(r^{3/2}\right)$ in (8), the asymptotic solution at a fixed point $(r_i,{\theta }_i)$ reads

$\mathbf{u}\left(r_i,{\theta }_i\right)=\frac{1}{\mu \sqrt{8\pi }}\left[\begin{array}{cc}\sqrt{r_i}\mathbf{F}\left({\theta }_i\right)& r_i\mathbf{I}\end{array}\right]\left[\begin{array}{c}\mathbf{K}\\ \mathbf{g}\left({\theta }_i\right)\end{array}\right]=\frac{1}{\mu \sqrt{8\pi }}\left[\begin{array}{cc}{\mathbf{M}}_{\sqrt{r}i}& {\mathbf{M}}_{ri}\end{array}\right]\left[\begin{array}{c}\mathbf{K}\\ {\mathbf{c}}_i\end{array}\right]$ (10)

Choosing $i=1\dots N$ arbitrary local points $(r_i,{\theta }_i)$, one writes

$\left[\begin{array}{c}\mathbf{u}\left(r_1,{\theta }_1\right)\\ ⋮\\ \mathbf{u}\left(r_N,{\theta }_N\right)\end{array}\right]=\frac{1}{\mu \sqrt{8\pi }}\left[\begin{array}{cccc}{\mathbf{M}}_{\sqrt{r}1}& {\mathbf{M}}_{r1}& & \\ ⋮& & \ddots & \\ {\mathbf{M}}_{\sqrt{r}N}& & & {\mathbf{M}}_{rN}\end{array}\right]\left[\begin{array}{c}\mathbf{K}\\ {\mathbf{c}}_1\\ ⋮\\ {\mathbf{c}}_N\end{array}\right]=\frac{1}{\mu \sqrt{8\pi }}\mathbf{M}\mathbf{\alpha }\mathbf{}\mathbf{}\mathbf{}$ (11)

in which the vector $\mathbf{\alpha }$ gathers all the unknown quantities. Since the matrix $\mathbf{M}$ is non-square, the best solution is given by (Strang, 2006)

$\mathbf{\alpha }=\frac{1}{\mu \sqrt{8\pi }}{\left({\mathbf{M}}^T\mathbf{M}\right)}^{-1}{\mathbf{M}}^T\left[\begin{array}{c}\mathbf{u}\left(r_1,{\theta }_1\right)\\ ⋮\\ \mathbf{u}\left(r_N,{\theta }_N\right)\end{array}\right]\mathbf{}$ (12)

that provides $K_I$, $K_{II}$ and $K_{III}$ at the crack front point $O\text{'}$. It is worth mentioning that the least-square solution provided by (12) leads always to a well-conditioned system if one sets $r_1=r_2{\dots =r}_N=\stackrel{-}{r}$. The radius $\stackrel{-}{r}$ 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 $\stackrel{-}{r}$ is within the range of the fracture mechanics theory. The DBEM enables to obtain the displacement values at the extraction points. For $N=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 $\phi =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 $\stackrel{-}{r}$ is performed.

Figure 3:
Internal points distribution fixed for the numerical applications.

2.3 Three-dimensional crack growth criteria

Propagation modelling according the LEFM requires the evaluation of an equivalent stress intensity factor $K_{\mathrm{e}\mathrm{q}}=K_{\mathrm{e}\mathrm{q}}\left(K_I,K_{II},K_{III}\right)$ and a propagation angle ${\theta }_p$ that drives the crack growth after $K_{\mathrm{e}\mathrm{q}}$ has reached the material toughness value $K_{I\mathrm{c}}$. Such conditions are mathematically expressed in the form

$F\left({\theta }_p\right)=K_{\mathrm{e}\mathrm{q}}\le K_{I\mathrm{c}}\left({\theta }_p\right)=0$ (13)

The propagation functions $F$ and $$ have been modelled by the maximum energy release rate (MERR) criterion and by the Schöllmann’s criterion (Schöllmann et al., 2002).

2.3.1 The maximum energy release rate criterion

According to this criterion, the energy release rate $G$ governs the crack growth: if the energy release rate associated to the loading pattern reaches the threshold value ${G\le G}_c$ then the crack propagates. The direction drived by the propagation angle ${\theta }_p$ is the one that maximizes the energy release rate $G$. Briefly,

$G\left({\theta }_p\right)=G_{\mathrm{c}}\frac{\partial G\left(\theta ={\theta }_p\right)}{\partial \theta }=0\frac{{\partial }^{2}G\left(\theta ={\theta }_p\right)}{\partial {\theta }^2}<0$ (14)

During the crack extension

$G\left(\theta \right)=G_{\theta }+G_r+G_{z\mathrm{\text{'}}}$ (15)

where

$G_{\theta }\left(\theta \right)=\frac{1}{2}\underset{∆a\to 0}{\lim }{\int }_0^{∆a}{\sigma }_{\theta }\left(a+\right)‖u_{\theta }\left(a+∆a-\right)‖d$ $G_r\left(\theta \right)=\frac{1}{2}\underset{∆a\to 0}{\lim }{\int }_0^{∆a}{\tau }_{\theta z\mathrm{\text{'}}}\left(a+\right)‖u_r\left(a+∆a-\right)‖d$ $G_{z\mathrm{\text{'}}}\left(\theta \right)=\frac{1}{2}\underset{∆a\to 0}{\lim }{\int }_0^{∆a}{\sigma }_{z\mathrm{\text{'}}}\left(a+\right)‖u_{z\mathrm{\text{'}}}\left(a+∆a-\right)‖d$ (16)

in which ${\sigma }_{\theta }$, ${\tau }_{\theta z\text{'}}$ and ${\sigma }_{z\text{'}}$ defines the stress field along the crack surfaces before the propagation and $‖u_i‖$ stands for the displacement discontinuity after the propagation.

Figure 4:
Crack propagation angles: kinking angle ${\theta }_p$ and twisting angle ${}_p$.

A closed-form solution for $G$ is not available in the literature. Chang et al. (2006) present an approximated expression, in which the energy release rate is written for $\theta =0$. Their solution can be readily adapted for mixed-mode propagation $(\theta \ne 0)$ if one assumes that

$K_{I\mathrm{e}\mathrm{q}}\left(\theta \right)=\sqrt{2\pi r}{\sigma }_{\theta }=\frac{\sqrt{2\pi r}}{4}\left[\left(3\cos \frac{\theta }{2}+\cos \frac{3\theta }{2}\right)K_I-3\left(\sin \frac{\theta }{2}+\sin \frac{3\theta }{2}\right)K_{II}\right]$ $K_{II\mathrm{e}\mathrm{q}}\left(\theta \right)=\sqrt{2\pi r}{\tau }_{\theta z\mathrm{\text{'}}}=\frac{\sqrt{2\pi r}}{4}\left[\left(\sin \frac{\theta }{2}+\sin \frac{3\theta }{2}\right)K_I+\left(\cos \frac{\theta }{2}+3\cos \frac{3\theta }{2}\right)K_{II}\right]$ $K_{III\mathrm{e}\mathrm{q}}\left(\theta \right)=\sqrt{2\pi r}{\sigma }_{z\mathrm{\text{'}}}=\sqrt{2\pi r}\cos \frac{\theta }{2}{K}_{III}$ (17)

which brings the following approximated expression:

$G\left(\theta \right)=\frac{1}{2\mu }\left[\left(1-\right)\left(K_{I\mathrm{e}\mathrm{q}}^2+K_{II\mathrm{e}\mathrm{q}}^2\right)+K_{III\mathrm{e}\mathrm{q}}^2\right]$ (18)

Substituting (18) into (14), taking into account (15)-(17), leads to

$F\left({\theta }_p\right)=\frac{1}{\sqrt{2}}\left|\cos \frac{{\theta }_p}{2}\right|\sqrt{\left(1+\cos {\theta }_p\right)K_I^2-4\sin {\theta }_p{K}_I{K}_{II}+\left(5-3\cos {\theta }_p\right)K_{II}^2+\frac{2K_{III}^2}{1-}}$ $\left({\theta }_p\right)=\frac{1+k}{8}\left[\left(\sin \frac{{\theta }_p}{2}+\sin \frac{3{\theta }_p}{2}\right)K_I^2+4\cos \frac{3{\theta }_p}{2}{K}_I{K}_{II}+\left(5\sin \frac{{\theta }_p}{2}-3\sin \frac{3{\theta }_p}{2}\right)K_{II}^2\right]+\sin \frac{{\theta }_p}{2}{K}_{III}^2$ (19)

2.3.2 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 ${\theta }_p$ and the twisting angle ${}_p$. 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

${\sigma }_1^{\mathrm{\text{'}}}=\frac{1}{2}\left({\sigma }_{\theta }+{\sigma }_{z\mathrm{\text{'}}}\right)+\frac{1}{2}\sqrt{{\left({\sigma }_{\theta }+{\sigma }_{z\mathrm{\text{'}}}\right)}^2+4{\tau }_{\theta z\mathrm{\text{'}}}^2}$ (20)

where ${\sigma }_{\theta }$, ${\tau }_{\theta z\text{'}}$ and ${\sigma }_{z\text{'}}$ 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

Figure 5:
Near crack front stress state in polar coordinates $(r,\theta ,z\text{'})$.

factors $\left(K_I,K_{II},K_{III}\right)$. Disregarding the contribution of ${\sigma }_{z\text{'}}$ into the kinking angle ${\theta }_p$, the Schöllmann’s criterion is formulated for ${\sigma }_{z\text{'}}=0$ assuming

$\frac{\partial {\sigma }_1^{\mathrm{\text{'}}}\left(\theta ={\theta }_p\right)}{\partial \theta }=0\frac{{\partial }^2{\sigma }_1^{\mathrm{\text{'}}}\left(\theta ={\theta }_p\right)}{\partial {\theta }^2}<0$ (21)

Finally, substitution of (20) into (21) yields

$F\left({\theta }_p\right)=\frac{1}{2}\left|\cos \frac{{\theta }_p}{2}\right|\left[K_I{\cos }^2\frac{{\theta }_p}{2}-\frac{3}{2}{K}_{II}\sin {\theta }_p+\sqrt{{\left(K_I{\cos }^2\frac{{\theta }_p}{2}-\frac{3}{2}{K}_{II}\sin {\theta }_p\right)}^2+4K_{III}^2}\right]$ $\left({\theta }_p\right)=-\frac{1}{8\sqrt{2\pi r}}\left({}_1+\frac{{}_1{}_2+16\sin {\theta }_p{K}_{III}^2}{\sqrt{{}_3}}\right)$ (22)

where

${}_1\left({\theta }_p\right)=\frac{3}{2}\left({}_4{K}_I+{}_5{K}_{II}\right){}_2\left({\theta }_p\right)={}_6{K}_I-3{}_4{K}_{II}{}_3\left({\theta }_p\right)={}_2^2+{\left(8\cos \frac{{\theta }_p}{2}{K}_{III}\right)}^2$ ${}_4\left({\theta }_p\right)=\sin \frac{{\theta }_p}{2}+\sin \frac{3{\theta }_p}{2}{}_5\left({\theta }_p\right)=\cos \frac{{\theta }_p}{2}+3\cos \frac{3{\theta }_p}{2}{}_6\left({\theta }_p\right)=3\cos \frac{{\theta }_p}{2}+\cos \frac{3{\theta }_p}{2}$ (23)

The ${\sigma }_2^{\text{'}}$ principal stress orientation enables the twisting angle determination

${}_p\left({\theta }_p\right)=\frac{1}{2}{\tan }^{-1}\left(\frac{2{\tau }_{\theta z\mathrm{\text{'}}}}{{\sigma }_{\theta }-{\sigma }_{z\mathrm{\text{'}}}}\right)$ (24)

which can be achieved after the kinking angle evaluation, as illustrated in Figure 6.

Figure 6:
Definition of the twisting angle ${}_p$.

3 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 $\left(K_I,K_{II},K_{III}\right)$. The accuracy of the proposed formulation has been demonstrated by a parametric study with respect to the extraction distance $r$. Moreover, the Maximum Energy Release Rate (MERR) and the Maximum Principal Stress (MPS) criteria determine the propagation angles and the equivalent SIF.

3.1 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, $\sigma =1$ (stress unit), opens the crack. The structure dimensions are functions of the crack length $a=2.5$, as illustrated by the ratios depicted in the figure.

Figure 7:
Mode I fracture problem.

Figure 8:
Meshes utilised for the problem discretisation.

The elastic constants adopted for the material are $E=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.

Figure 9:
Norm of the displacement vector obtained with the model for the three meshes.

Figure 10:
Mode I SIF results.

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 $r/a$, in which $a$ 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 ${\mathbf{u}}_{\sqrt{r}}\left(r,\theta \right)$. In this case, the matrix $\mathbf{M}$ reduces to $M_{\sqrt{r}}$. The DFT-2 considers both the displacement terms ${\mathbf{u}}_{\sqrt{r}}\left(r,\theta \right)$ and ${\mathbf{u}}_r\left(r,\theta \right)$. As shown in Figure 10, the inclusion of the term ${\mathbf{u}}_r\left(r,\theta \right)$ 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 $r$, when compared with the DFT-1. These evaluations have been performed for points near the crack front from$r/a=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 $r/a$ 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.

3.2 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: $L_t=260$, $L=240$, $H=60$ and $t=10$. Moreover, the loading is equal to $F=2$. The initial notch length is $a_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 $E=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.

Figure 11:
Three-point bending beam with an initial kinked notch.

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 $r=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 $‖K_i‖=K_i\left(2H^{2}t/3L\sqrt{\pi a_0}\right)$.

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 $r=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 ${\theta }_p$, ${}_p$ and the equivalent stress intensity factor $K_{\mathrm{e}\mathrm{q}}$ obtained for the MERR and the MPS criteria for $r=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.

Figure 12:
Mesh utilised for the problem discretisation.

Figure 13:
Norm of the displacement vector obateined with the model.

3.3 Penny-shaped crack problem

The last application handles a penny-shaped crack surface, with radius $r$, embedded in an infinity elastic medium and subjected to a remote stress $\sigma$, as shown in Figure 16. The analytical SIF solution for this problem was presented by Murakami (1987), which is described by

$K_I={\frac{2}{\pi }\sin }^2\gamma \sqrt{\pi a}\sigma K_{II}=\frac{2}{\pi \left(2-\right)}\sin 2\gamma \sin \theta \sqrt{\pi a}\sigma K_{III}=\frac{2}{\pi \left(1-\right)\left(2-\right)}\sin 2\gamma \cos \theta \sqrt{\pi a}\sigma$ (25)

where $\gamma$ 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 ($\gamma =90°$) and mixed mode case ($\gamma =45°$). The elastic constants adopted for the analysis are $E=10$ and $=0$. The remote stress was assumed as unitary $\sigma =1$. The penny-shaped crack surface radius is $a=0.1$ (data in consistent but unspecified units).

Figure 14:
SIF results along the thickness of the beam.

Figure 15:
Propagation angles and equivalent SIF along the solid thickness.

Figure 16:
Penny-shaped crack surface in an infinity medium.

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 $K_I$ and the extraction distance $r$. Further, the mixed mode case was studied and it was possible to compare the results of all the SIFs $\left(K_I,K_{II},K_{III}\right)$. Figure 18 presents the deformed crack surfaces for both mode cases with displacements 100 times magnified.

3.3.1 Mode I case

For the mode I ($\gamma =90°$) scenario, it was analysed the influence of the extraction distance $r$ 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 $K_I$ versus $r/a$ 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 $K_I$ versus the $\theta$ 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 ${\theta }_p$, ${}_p$ and the equivalent stress intensity factor $K_{\mathrm{e}\mathrm{q}}$ obtained for the MERR and the MPS criteria for the model quadratic M2.

Figure 17:
Crack surface meshes of DBEM models.

Figure 18:
Deformed crack surfaces for $\gamma =90°$ and $\gamma =45°$.

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 $K_I$.

3.3.2 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 ($\gamma =45°$). The results of $\left(K_I,K_{II},K_{III}\right)$ versus $\theta$ 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 $K_{II}$ and $K_{III}$ with linear and with quadratic elements are attributed to the better approximation of quadratic elements in curved geometries. The lack of accuracy for $K_{II}$ and $K_{III}$ with linear elements is, therefore, due to the lack of accuracy in the definition of the local crack tip coordinate systems.

Figure 19:
Results of $K_I$ versus $r/a$ for the DBEM models.

Figure 20:
Results of $K_I$ versus $\theta$ for the DBEM models.

Figure 23 illustrates the propagation angles ${\theta }_p$, ${}_p$ and the equivalent stress intensity factor $K_{\mathrm{e}\mathrm{q}}$ 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.

Figure 21:
Propagation angles and equivalent stress intensity factor versus $\theta$ for the quadratic M2 model.

Figure 22:
Results of $\left(K_I,K_{II},K_{III}\right)$ versus $\theta$ for the DBEM models.

Figure 23:
Propagation angles and equivalent stress intensity factor versus $\theta$ for the quadratic M2 model.

4 CONCLUSION

This study presented the extension of the Displacement Fitting Technique (DFT) for SIFs assessment in three-dimensional 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 $r/a=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 $r$, when the evaluation points are equally distributed around a circular path of radius $r$ 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.

Acknowledgements

The authors are thankful to the funding provided by the Coordination of Improvement of Higher Education Personnel (CAPES) and to the Sponsorship of this research project by the National Council for Scientific and Technological Development (CNPq) grant number 141565/2015-2.

References

• Balderrama, R., Cisilino, A.P., Martinez, M. (2006). Boundary element method analysis of three-dimensional thermoelastic fracture problems using the energy domain integral, Journal of Applied Mechanics 73(6):959-969.
• Banks-Sills, L., Hershkovitz, I., Wawrzynek, P.A, Eliasi, R., Infraffea, A.R. (2005). Methods for calculating stress intensity factors in anisotropic materials: Part I - z=0 is a symmetric plane, Engineering Fracture Mechanics 72:2328-2358.
• Banks-Sills, L., Sherman, D. (1986). Comparison of methods for calculating stress intensity factors with quarter-point elements, International Journal of Fracture 32:127-140.
• Bezerra, L.M., Medeiros, J.M.S. (2002). Using boundary elements and J-integral for the determination of KI in fracture mechanics. Proceedings of IABEM 2002 Symposium, Austin, USA.
• Chan, S.K., Tuba, I.S., Wilson, W.K. (1970). On the finite element method in linear fracture mechanics, Engineering Fracture Mechanics 2(1):1-17.
• Chang, J., Xu, J.Q., Mutoh, Y. (2006). A general mixed-mode brittle fracture criterion for cracked materials, Engineering Fracture Mechanics 73(9):1249-1263.
• Chen, H.K., Hong, J.T. (1999). Review of dual boundary element methods with emphasis of hypersingular integrals and divergent series, Applied Mechanical Reviews 52(1):17-33.
• Citarella, R., Buchholz, F.G. (2008). Comparison of crack growth simulation by DBEM and FEM for SEN-specimens undergoing torsion or bending loading, Engineering Fracture Mechanics 75(3-4):489-509.
• Cordeiro, S.G.F., Leonel, E.D. (2018). Mechanical modelling of three-dimensional cracked structural components using the isogeometric dual boundary element method, Applied Mathematical Modelling 63:415-444.
• Cordeiro, S.G.F., Leonel, E.D. (2019). An improved computational framework based on the dual Boundary Element Method for three-dimensional mixed-mode crack propagation analyses, Advances in Engineering Software 135:102689.
• Cruse, T.A. (1996). BIE fracture mechanics analysis: 25 years of developments, Computational Mechanics 18:1-11.
• Cruse, T.A., Wilson, R.B. (1978). Advanced applications of the boundary integral equation method, Nuclear Engineering and Design 46(1):223-234.
• Dell’Erba, D.N., Aliabadi, M.H. (2001). BEM analysis of fracture problems in three-dimensional thermoelasticity using J-integral, International Journal of Solids and Structures 38(26-27):4609-4630.
• Dirgantara, T., Aliabadi, M.H. (2002). Stress intensity factors for cracks in thin plates, Engineering Fracture Mechanics 69(13):1465-1486.
• Gonzalez-Albuixech, V.F., Giner, E., Tarancon, J.E., Fuenmayor, F.J., Gravouil, A. (2013). Domain integral formulation for 3-d curved and non-planar cracks with the extended finite element method, Computer Methods in Applied Mechanics and Engineering 264:129-144.
• Gonzalez, M., Teixeira, P., Wrobel, L.C., Martinez, M. (2015). A new displacement-based approach to calculate stress intensity factors with the boundary element method, Latin American Journal of Solids and Structures 12(9):1677-1697.
• Guiggiani, M., Krishnasamy, G., Rudolphi, T.J., Rizzo, F.J. (1992). A general algorithm for the numerical solution of hypersingular boundary integral equations, Journal of Applied Mechanics 59(3):604-614.
• Gupta, P., Duarte, C.A. (2018). Coupled hydro-mechanical-fracture simulations of non-planar three-dimensional hydraulic fracture propagation, International Journal for Numerical and Analytical Methods in Geomechanics 42(1):143-180.
• Gupta, P., Duarte, C.A., Dhankhar, A. (2017). Accuracy and robustness of stress intensity factor extraction methods for the generalized/extended finite element method, Engineering Fracture Mechanics 179:120-153.
• Ingraffea, A.R., Manu, C. (1980). Stress-intensity factor computation in three dimensions with quarter-point elements, International Journal for Numerical Methods in Engineering 15(10):1427-1445.
• Mi, Y., Aliabadi, M.H. (1992). Dual boundary element method for three-dimensional fracture mechanics analysis, Engineering Analysis with Boundary Elements 10(2):161-171.
• Minnebo, H., Maiérus, J., Noels, L. (2010). Displacement extrapolation method: An alternative to J-integral for stress intensity factors computation using X-FEM, Proceedings of IV European Congress on Computational Mechanics: Solids, Structures and Coupled Problems in Engineering, Paris, France.
• Murakami, Y. (1987). Stress Intensity Factors Handbook, Pergamon Press (Oxford).
• Nejati, M., Paluszny, A., Zimmerman, R.Q. (2015). On the use of quarter-point tetrahedral finite elements in linear elastic fracture mechanics, Engineering Fracture Mechanics 144: 194-221.
• Ortiz, J.E., Cisilino, A.P. (2005). Boundary element method for J-integral and stress intensity factor computations in three-dimensional interface cracks, International Journal of Fracture 133:197-222.
• Portela, A., Aliabadi, M.H., Rooke, D.P. (1992). The dual boundary element method: Effective implementation for crack problems, International Journal for Numerical Methods in Engineering 33(6):1269-1287.
• Purbolaksono, J., Dirgantara, T., Aliabadi, M.H. (2012). Fracture mechanics analysis of geometrically nonlinear shear deformable plates, Engineering Analysis with Boundary Elements 36(2):87-92.
• Rigby, R.H., Aliabadi, M.H. (1998). Decomposition of the mixed-mode J-integral - Revisited, International Journal of Solids and Structures 35(17):2073-2099.
• Richard, H.A., Fulland, M., Sander, M. (2005). Theoretical crack path prediction, Fatigue & Fracture of Engineering Materials & Structures 28(1-2):3-12.
• Schätzer, M., Fries, T.P. (2016). Stress Intensity Factors through Crack Opening Displacements in the XFEM, Advances in Discretization Methods: Discontinuities, Virtual Elements, Fictitious Domain Methods, 12, pp. 143-164. Berlin: Springer Vieweg.
• Schöllmann, M., Richard, H.A., Kullmer, G., Fulland, M. (2002). A new criterion for the prediction of crack development in multiaxially-loaded structures, International Journal of Fracture 117(2):129-141.
• Sih, G.C. (1991). Mechanics of Fracture Initiation and Propagation: Surface and volume energy density applied as failure criterion, Kluwer Academic (Massachusetts).
• Strang, G. (2006). Linear Algebra and its applications, 4th edn., Cengage Learning.
• Walters, M.C., Paulino, G.H., Dodds Jr., R.H. (2005). Interaction-integral procedures for 3-D curved cracks including surface tractions, Engineering Fracture Mechanics 72:1635-1663.

Edited by

Guest Editors:

Volnei Tita and Nicholas Fantuzzi.

Publication Dates

• Publication in this collection
  09 Nov 2020
• Date of issue
  2020

History

• Received
  09 Mar 2020
• Reviewed
  10 July 2020
• Accepted
  28 Aug 2020
• Published
  01 Sept 2020