## Abstract

The topological derivative describes the variation of a response functional with respect to infinitesimal changes in topology, such as the introduction of an infinitesimal crack or hole. In this three-dimensional fracture mechanics work, we propose an approximation of the energy release rate field associated with a small surface crack of any boundary location, direction, and orientation combination using the topological derivative. This work builds on the work of Silva et al. (“Energy Release Rate Approximation for Small Surface-Breaking Cracks Using the Topological Derivative,” J. Mech. Phys. Solids **59**(5), pp. 925–939), in which the authors proposed an approximation of the energy release rate field which was limited to two-dimensional domains. The proposed method is computationally advantageous because it only requires a single analysis. By contrast, current boundary element and finite element-based methods require an analysis for each crack length-location-direction-orientation combination. Furthermore, the proposed method is evaluated on the non-cracked domain, obviating the need for refined meshes in the crack tip region.

## 1 Introduction

The energy release rate *G*, which is defined as the loss of total potential energy as an infinitesimal crack grows, has been the principal quantity used to predict the failure of cracked bodies since its introduction by Griffith [1]. It represents the energy that is available for the crack extension, which takes place when *G* is equal to the material’s fracture toughness *G*_{C}. Due to the importance of the energy release rate in the failure analysis and design of structures, multiple methods have been proposed to compute *G* for stationary cracks in linear elastic domains, which is the subject of this work.

The first set of methods rely on the relation between *G* and the stress intensity factors *K*_{λ}, which describes the displacement and stress fields in the neighborhood of the crack [2]. In general, the stress intensity factor computations require a precise description of the near-tip stress and strain fields. The most prominent method used to that effect is the finite element method. However, it is well known that crack tip stress singularities limit the use of standard finite element methods and necessitate a refined mesh in the crack tip region. Specialty finite element methods with adaptive meshing [3,4], enriched p-refinement, basis functions, and fractal meshes [5–7], least-square enhancements [8,9], quarter-point and crack tip displacement discontinuity elements [10,11], and submodels and subregions [12,13] have been used. Finite element method variants including the element free Galerkin method [14,15] and extended finite method [16,17] have also been used. Other prominent methods to compute *K*_{λ} include the boundary element method [18–20] and weight function method [21,22]. In the weight function method, the stress intensity factors for any loading system applied to a body can be calculated by integration. Universal weight functions have been introduced to simplify these calculations [23]. Comparisons among some of these methods are made in Refs. [24–26].

Introduced by Rice [27], the path-independent *J*-integral can be used to calculate the energy release rate [28,29]. The *J*-integral relates the energy flux through an arbitrary contour enclosing the crack front to the far-field loading. The interaction energy integral method [30,31] can also be used to calculate the energy release rate. This method uses superposition to combine the responses corresponding to several boundary value problems. Evaluation of a conservation integral for this superimposed response produces expressions for the energy release rate. The virtual crack extension method [32–34] models crack extension by simulating the responses on the cracked domain and a perturbed domain in which the crack is extended. The energy release rate is subsequently calculated in a post-processing operation. In this method, the energy release rate and its higher-order derivatives for multiple cracked systems can be computed in one analysis.

The topological derivative describes the variation of a response functional with respect to infinitesimal changes in topology, such as the introduction of an infinitesimal crack or hole [35]. The first application of the topological derivative was in the so-called bubble method. In this method, holes are nucleated to reduce the weight of the structure without reducing its integrity. The holes are then enlarged and reconfigured by traditional shape optimization methods. The topological derivative has since been applied to image processing [36,37], inverse problems [38–40], and structural and topology optimization [41–43]. It has also been applied to numerous fracture mechanics problems; in Ref. [44], the topological derivative is used in an inverse heat conduction problem to detect and locate cracks. The shape sensitivities and topological derivatives are reported for cracks in elastic bodies on boundaries of rigid inclusions in Ref. [45]. Crack nucleation and crack growth criterion are introduced through an asymptotic expansion based on the topological derivative associated with the total potential energy of a linear elastic body in Ref. [46]; in Ref. [47], the authors used higher-order topological derivatives to approximate the energy release rate for edge cracks of arbitrary length in two-dimensional domains.

In the present work, we develop an approximation of the energy release rate field associated with a small half-penny shaped surface crack of any boundary location, direction, and orientation combination using the topological derivative in a linear elastic isotropic three-dimensional domain; in Ref. [48], the authors developed a similar approximation of the energy release rate field which was limited to two-dimensional domains. The proposed method is computationally advantageous because it only requires a single analysis. By contrast, current boundary element and finite element-based methods require an analysis for each crack length-location-direction-orientation combination. Furthermore, the proposed method is evaluated on the non-cracked domain, obviating the need for refined meshes in the crack tip region. There are multiple methods to calculate the topological derivative, including those presented in Refs. [35,41,49]. In this work, we use the shape sensitivity method [42], rather than, e.g., the domain truncation method [49], to compute the topological derivative due to its simplicity and generality with respect to the hole shape and hole boundary conditions.

The topological derivative and its application to fracture mechanics are described in Sec. 2. An asymptotic analysis of the stress field is presented in Sec. 3. The energy release rate approximation is detailed in Sec. 4. In Sec. 5, examples that demonstrate the efficacy of the method are presented. Finally, conclusions are detailed in Sec. 6.

## 2 Topological Derivative

We consider an isotropic, linearly elastic, and three-dimensional domain Ω, cf. Fig. 1. The boundary, outward normal vector, and prescribed traction are given by ∂Ω, **n**, and **t**^{p}, respectively. Without loss of generality, we limit our discussion to a single surface crack geometry, the half-penny. When a small half-penny crack is introduced at boundary location $x^$ with radius ɛ, direction **d**, and orientation *α*, we obtain the cracked domain $\Omega \epsilon (\epsilon ,x^,d,\alpha )$. The boundary of $\Omega \epsilon (\epsilon ,x^,d,\alpha )$ is given by $\u2202\Omega \epsilon (\epsilon ,x^,d,\alpha )=\u2202\Omega \u222a\gamma \epsilon (\epsilon ,x^,d,\alpha )$ where $\gamma \epsilon (\epsilon ,x^,d,\alpha )$ is the crack surface with outward normal vector **m** and crack front $\u2202\gamma \epsilon (\epsilon ,x^,d,\alpha )$.^{1}

*G*is defined as the loss of total potential energy as an infinitesimal crack grows [1]. The total potential energy is defined as

**u**is the displacement vector, $\u2207su:=12(\u2207u+\u2207uT)$ is the infinitesimal strain tensor, $T=C[\u2207su]$ is the symmetric Cauchy stress tensor, and $C$ is the elasticity tensor. Response quantities which are defined on the cracked domain Ω

_{ɛ}, e.g.,

**u**, are denoted by the subscript

_{ɛ}*ɛ*. The energy release rate can be expressed as the shape sensitivity of Eq. (1) with respect to the crack radius

*ɛ*[50,51], i.e.,

_{ɛ}:

**t**

^{p}is defined such that global equilibrium is satisfied:

**z**is an arbitrary location, not necessarily in the domain.

*G*of Eq. (2) for a half-penny crack of any boundary location $x^$, crack radius

*ɛ*, direction

**d**, and orientation

*α*. To this end, a topological asymptotic expansion is constructed to express the total potential energy as:

*j*th-order topological derivative of the total potential energy

*ψ*[52]. The parameter

*s*∈ [0, 1], which describes the position on the crack front ∂

*γ*

_{ɛ}, is defined as

*L*is the arc length (Fig. 2). The gauge functions

*f*

_{(j)}in Eq. (5) are functions of the crack radius

*ɛ*and the crack surface boundary conditions. They monotonically tend to zero as

*ɛ*tends to zero and satisfy

*ɛ*, direction

**d**, and orientation

*α*is assumed to exist [42]. A shape sensitivity analysis is performed on the crack with respect to

*ɛ*and the limit is then taken as

*ɛ*goes to zero to give

*f*

_{(j)}must be defined such that $DT(j)\psi (x^,d,\alpha ,s)$ is finite and nonzero, as we assume the total potential energy

*ψ*(Ω

_{ɛ}) is always bounded. Furthermore, we define a negative valued function

*f*

_{(j)}such that the topological derivative of the total potential energy defined in Eq. (1) is positive. We also see from Eqs. (5) and (8) that if

*f*

_{(j)}is a suitable gauge function, then

*ζf*

_{(j)}is also suitable for any constant

*ζ*> 0 and $DT(j)\psi (x^,d,\alpha ,s)$ is unaffected by the choice of

*ζ*.

*ψ*as a function of the topological derivatives, we differentiate Eq. (5) with respect to

*ɛ*:

*G*as a function of the topological derivatives:

*ɛ*, location $x^$, direction

**d**, and orientation

*α*combination using response quantities computed on the non-cracked domain Ω.

*j*= 1.

^{2}

**e**

_{n}(

*s*),

**e**

_{t}(

*s*), and

**e**

_{b}(

*s*) which are the outward normal, tangent, and binormal directions on the crack front ∂

*γ*

_{ɛ}(Fig. 2).

^{3}The direction of crack growth is presumed to be in the normal direction

**e**

_{n}such that the shape derivative is defined as [51]

*η*=

*η*(

*s*) defines the size of the perturbation and

**Σ**is the energy-momentum tensor

*J*-integral [27]. Alternatively, the

*J*-integral may be evaluated as [53]

*K*

_{λ}are the mode

*λ*stress intensity factors which describe the displacement and stress fields in the neighborhood of the crack [2], $E\xaf=E/(1\u2212\nu 2)$,

*E*is the Young’s modulus,

*ν*is the Poisson’s ratio, and

*μ*is the shear modulus.

## 3 Asymptotic Analysis

**T**(

**x**) denotes the outer stress, $T~(y)$ denotes the inner stress, and $T\u02c7\epsilon (x)$ denotes the remainder stress (Fig. 3). Note that the inner stress $T~(y)$ is computed at the scaled position vector

**y**=

**x**/

*ɛ*. By construction, the boundary conditions of Eq. (3) are satisfied by the stress composite expansion, i.e.,

**e**

_{b}=

**m**,

**T**(

**x**) satisfies the nominal problem, i.e., it is defined on the non-cracked domain Ω and satisfies the boundary value problem,

**T**(

**x**) does not satisfy the traction-free boundary condition on the crack surface

*γ*

_{ɛ}(Eq. (19)).

*ɛ*is small, we expand the outer stress

**T**(

**x**) about the crack initiation site $x^$ in the

**e**

_{n}direction, i.e., $x\u2212x^=\delta en$, to obtain

*D*represents the Fréchet derivative, which is defined as the linear operators

*ξ*=

*δ*/

*ɛ*(Fig. 3) to express Eq. (21) as

**y**, it satisfies the boundary value problem of a half-space $H$ with a unit radius half-penny crack (crack surface

*γ*

_{1}). There is a non-zero crack face traction on

*γ*

_{1}and zero far-field traction, i.e.,

**T**(

**x**)

**e**

_{b}on the crack surface

*γ*

_{ɛ}due to the outer stress. Therefore, using Eq. (23), the prescribed traction $t~$ is defined as

_{ɛ}and satisfies the boundary value problem,

**t**

^{R}are defined as

**e**

_{n},

**e**

_{t},

**e**

_{b}) coordinate system, the inner stress can be expressed as

*ρ*is the radial distance from the crack front,

**K**is a tensor,

*ϕ*is the angle between

**e**

_{ρ}and the crack plane with binormal vector

**e**

_{b}, and

**F**(

*ϕ*) is a dimensionless tensor function of

*ϕ*[55]. Notably, $T~(x/\epsilon )=O(\epsilon )$ far from the crack front, in particular for

**x**∈ ∂Ω. Furthermore, we observe from Eq. (23) that

## 4 Energy Release Rate Algorithm

**e**

_{r},

**e**

_{z}, and $e\theta $ which are the outward normal, tangent, and binormal directions at the position

*s*= 0.5 (Fig. 2), i.e.,

**T**(

**x**), the normalized factors are expressed as

*h*

_{ij}(

*α*,

*s*) are the weight functions, dimensionless functions of the orientation

*α*and location

*s*and, e.g., the $Tr\theta =er\u22c5Te\theta $ are the crack face traction components $t~=T(x^)e\theta $, cf. Eq. (25).

*W*× 2

*W*× 2

*W*cubic domain. It was determined that the dimension

*W*≥ 100

*ɛ*was required for the numerical approximation to sufficiently represent a half-space. Three traction boundary value problems with crack face tractions of $t~=e\theta $, $t~=er$, and $t~=ez$ were then solved. These problems were used to calculate

*h*

_{i1}(

*α*,

*s*),

*h*

_{i2}(

*α*,

*s*), and

*h*

_{i3}(

*α*,

*s*) for

*i*= {1, 2, 3}, respectively. For example, when $t~=e\theta $ then $Tr\theta (x^,d,\alpha )=0$ and $Tz\theta (x^,d,\alpha )=0$, and Eq. (33) reduces to

*h*

_{ij}(

*α*,

*s*) were calculated for $\u221275deg\u2264\alpha \u226475deg$ in increments of $1deg$. In this manner, more than 200 numerical analyses were required to calculate the weight functions presented in Fig. 4. We emphasize that this is an off-line computation that need not be repeated. In Fig. 5, the function values $h11(0deg,s)$ are presented and compared with values computed by Tada et al. using the alternating and finite element methods [53]. The difference between the two computations does not exceed $\xb11.13%$.

*f*(

*ɛ*) is defined so that Eq. (35) evaluates to a non-zero finite number cf. Eq. (8)

*f*(

*ɛ*), the topological derivative becomes

In summary, we perform the following steps to compute the energy release rate approximation $GTD(\epsilon ,x^,d,\alpha ,s)$:

Perform an analysis on the non-cracked domain Ω with prescribed traction

**t**^{p}to evaluate the outer stress**T**(**x**).Evaluate $K\xafI(x^,d,\alpha ,s)$, $K\xafII(x^,d,\alpha ,s)$, and $K\xafIII(x^,d,\alpha ,s)$ via Eq. (33) and the tabulated

*h*_{ij}(*α*,*s*) values.Evaluate $GTD(\epsilon ,x^,d,\alpha ,s)$ via Eq. (38).

Notably, we know the value of $GTD(\epsilon ,x^,d,\alpha ,s)$ for any crack radius *ɛ*, location $x^$, direction **d**, orientation *α*, and position *s* from a single analysis on the non-cracked domain Ω.

## 5 Examples

In the examples presented hereafter, we perform a single boundary element analysis on the non-cracked domain Ω to compute the energy release rate approximation *G*^{TD} for any crack radius *ɛ*, location $x^$, direction **d**, orientation *α*, and position *s*. To verify the proposed approximation, reference energy release rates *G*^{BE} and stress intensity factors $K\lambda BE$ are calculated by using the boundary element method on the cracked domain Ω_{ɛ} with a refined mesh in the crack tip region. We show that *G*^{TD} is in good agreement with *G*^{BE}.

The numerical analyses in this work were performed on a workstation running Microsoft^{®} Windows^{®} 10 Enterprise with a Dual Intel^{®} Xeon^{®} Processor E5-2609 v4 and 128 GB DDR4 RAM using BEASY software’s boundary element solver and fracture simulation tools. Convergence studies were performed on the cracked domain Ω_{ɛ} and non-cracked domain Ω to ensure accurate results.

### 5.1 Curved Beam Subjected to Shear and Tensile Loads.

*E*= 207 GPa,

*ν*= 0.3, $T=2m$,

*R*

_{i}= 5 m,

*R*

_{o}= 10 m,

*L*= 25 m, and

*σ*= 100 MPa. The stress intensity factors and energy release rate are calculated when cracks are introduced at boundary locations $x^1$ and $x^2$. The superscript $\u22c6$ denotes the (

**d**,

*α*,

*s*) that maximizes

*G*

^{TD}for a given

*ɛ*, hereafter referred to as the critical crack parameters, i.e.,

*α*and position

*s*in Eq. (39) are chosen to reduce free surface effects.

The stress intensity factors and energy release rate are first calculated for $s\u2208S$ when cracks of radius *ɛ*/*T* = 0.01, direction $d\u22c6(x^)$, and orientation $\alpha \u22c6(x^)$ are introduced. Figures 7(a)–7(c) compare the stress intensity factor approximations $KITD$, $KIITD$, and $KIIITD$ with the reference boundary element values $KIBE$, $KIIBE$, and $KIIIBE$. The difference between $KITD$ and $KIBE$ does not exceed $\xb11.44%$. When a crack is introduced at $x^1$, the stress intensity factors *K*_{II} and *K*_{III} are two orders of magnitude less than *K*_{I}, i.e., the energy release rate is controlled by *K*_{I}. When a crack is introduced at $x^2$, $KIITD$ is within $\xb11.32%$ of $KIIBE$ and $KIIITD$ is within $\xb11.26%$ of $KIIIBE$. This excludes positions near *s* = 0.5 where $KIITD$ and $KIIBE$ are both 0. Figure 7(d) compares the proposed energy release rate approximation *G*^{TD} with the reference boundary element value *G*^{BE}. The difference between *G*^{TD} and *G*^{BE} does not exceed $\xb11.67%$. Notably, the energy release rates *G*^{TD} and *G*^{BE} are proportional to the sum of squares of the stress intensity factors, cf. Eqs. (2), (16), and (38). As a result, the difference in the energy release rates is generally greater than that in the stress intensity factors. Finally, we note that the critical crack positions $s\u22c6(x^)$ determined using the topological derivative and boundary element method are the same.

*ɛ*/

*T*= 0.01 and orientation $\alpha \u22c6(x^)$ are introduced. To visualize these results, we define the angle $\beta (x^,d)$ which defines

**d**such that

Figure 8(a) compares $GmaxTD$ and $GmaxBE$, where it is seen that similar critical crack directions $d\u22c6(x^)$ are determined using the topological derivative and boundary element methods. One can also see that $GmaxTD$ is in good agreement with $GmaxBE$. $GmaxTD$ is within $\xb13.84%$ of $GmaxBE$ outside of the range $100deg\u2264\beta (x^,d)\u2264120deg$. Within this range, the small values of $GmaxTD$ and $GmaxBE$ result in larger discrepancies.

The energy release rate is then calculated for $(\alpha ,s)\u2208(A\xd7S)$ when cracks of radius *ɛ*/*T* = 0.01 and direction $d\u22c6(x^)$ are introduced. $GmaxTD$ and $GmaxBE$ are compared in Fig. 8(b), where it is again seen that the critical crack orientations $\alpha \u22c6(x^)$ determined using the topological derivative and boundary element methods are similar. Once more one can see that $GmaxTD$ is in good agreement with $GmaxBE$. When $\alpha \u226440deg$, $GmaxTD$ is within $\xb12.77%$ of $GmaxBE$. The difference between $GmaxTD$ and $GmaxBE$ grows as large as $\xb17.97%$ when $\alpha \u226550deg$. One conjecture is that when $\alpha \u226550deg$, the free surface effects become significant and the boundary element analyses on the cracked domain Ω_{ɛ} are less accurate.

To test this conjecture, the energy release rate is finally calculated for $s\u2208S$ when cracks of various radii *ɛ* are introduced with direction $d\u22c6(x^)$ and orientation $\alpha \u22c6(x^)$. The energy release rate calculations $GmaxTD$ and $GmaxBE$ are compared in Fig. 9. When the crack radius is less than $5%$ of the thickness, i.e., *ɛ*/*T* ≤ 0.05, the difference between $GmaxTD$ and $GmaxBE$ is small. As expected, for larger cracks, the difference between $GmaxTD$ and $GmaxBE$ increases. However, $GmaxTD$ remains a reasonable approximation of $GmaxBE$ even for cracks of considerable radius, e.g., *ɛ*/*T* = 0.25.

As indicated in Fig. 9, the energy release rate approximation $GmaxTD$ associated with the boundary location $x^1$ is more accurate for longer crack radii than the same approximation associated with the boundary location $x^2$. This result is associated with the accuracy of the stress composite expansion introduced in Sec. 3. This accuracy is influenced by the boundary curvature and boundary location $x^$ [47,48]. Higher-order approximations of the energy release rate can be developed using higher-order topological derivatives. Such approximations were developed by the authors for edge cracks in two-dimensional domains in Ref. [47]. These higher-order approximations take into account the boundary curvature, boundary location $x^$, and derivatives of the stress field. They allow the analyst to accurately consider cracks of larger radii and determine the crack radii for which the first-order approximation is accurate. On the other hand, the computation of higher-order approximations requires additional analyses.

A boundary element analysis on the cracked domain Ω_{ɛ} was required to compute each reference value of $GmaxBE$ presented in Fig. 9. The average computational time of these analyses was 25.5 min, i.e., the computational time required to produce the reference values in Fig. 9 was 1275 min per crack location. The average number of boundary elements was 17,600, depending on the boundary location $x^$ and crack radius *ɛ*. In contrast, the energy release rate approximation $GmaxTD$ was computed using one boundary element analysis with 1260 boundary elements on the non-cracked domain Ω. The computational times of this analysis and the post-processing required to compute $GmaxTD$ were 2.5 min and 0.25 min per crack, respectively. In other words, the computational speedup is 464 for the first crack location. For each additional crack location, the computational speedup is 5100.

*D*

_{T}

*ψ*for a given crack radius

*ɛ*, i.e.,

*G*

^{TD}is linearly proportional to the topological derivative

*D*

_{T}

*ψ*, cf. Eqs. (13) and (36). Thus, cracks of radius

*ɛ*will propagate at locations $x^$ if $\epsilon \u2265GC/DT\psi \u22c6(x^)$, where

*G*

_{C}is the fracture toughness of the material. As seen in Fig. 10, there are two regions where $DT\psi \u22c6(x^)$ is especially large, located on the ∂Ω

_{1}− ∂Ω

_{2}and ∂Ω

_{1}− ∂Ω

_{3}boundaries. Small cracks appearing in these regions will lead to failure. We note that the results presented in Fig. 10 were calculated from a single boundary element analysis on the non-cracked domain Ω. To achieve the same plots using the boundary element method on the cracked domain Ω

_{ɛ}would require millions of simulations, each more expensive than this analysis on the non-cracked domain Ω.

_{1}, $d\u22c6(x^)$ is generally between $0deg$ and $40deg$. On the boundaries ∂Ω

_{2}and ∂Ω

_{3}, $d\u22c6(x^)$ varies between the lower bound of $0deg$ and upper bound of $180deg$. For most locations $x^$, $d\u22c6(x^)$ is the direction which maximizes the local opening stress component $T\theta \theta (x^,d,\alpha )$ and mode I stress intensity factor approximation $KITD$. Contrary to the results for $d\u22c6(x^)$, there is limited variation in the critical orientation $\alpha \u22c6(x^)$ on ∂Ω

_{1}, ∂Ω

_{2}, and ∂Ω

_{3}. Excluding the edges of the boundaries, $\alpha \u22c6(x^)$ is generally between $\u22126deg$ and $6deg$, i.e., nearly perpendicular to the boundaries. For all locations $x^$, the critical position $s\u22c6$ is equal to 0.1 or 0.9, i.e., it is as close to the free surface as allowed by the bounds of $S$. This is predictable, as the absolute values for seven of the nine weight functions

*h*

_{ij}(

*α*,

*s*) are maximized when

*s*equals 0.1 and 0.9, cf. Fig. 4.

### 5.2 Coupling Subjected to Axial and Radial Loads.

In this example, we study a coupling subjected to axial and radial loads (Fig. 11). The parameters for this numerical analysis are *E* = 207 GPa, *ν* = 0.3, *R*_{1} = 1 m, *R*_{2} = 1.5 m, *R*_{3} = 2 m, *R*_{4} = 0.075 m, *h*_{1} = 1 m, *h*_{2} = 0.5 m, and *σ* = 100 MPa. The stress intensity factors and energy release rate are calculated for $s\u2208S$ when cracks of various radii *ɛ* are introduced at boundary location $x^3=1.00e1+0.250e3$ with direction **d** = 1.00**e**_{3} and orientation $\alpha =45deg$. The stress intensity factors versus the position *s* are shown in Figs. 12(a)–12(c) for a crack of radius *ɛ* = (*R*_{2} − *R*_{1})/100. As in the previous example, one can see the stress intensity factors $KITD$, $KIITD$, and $KIIITD$ computed using the topological derivative are in good agreement with the reference boundary element values $KIBE$, $KIIBE$, and $KIIIBE$. The energy release rates appear in Fig. 12(d), where it is seen that the proposed approximation *G*^{TD} is within $\xb15.62%$ of the reference boundary element value *G*^{BE}. The critical crack positions $s\u22c6(x^)$ determined using the topological derivative and boundary element methods are again the same. Figure 13 shows the energy release rates versus the crack radius *ɛ*/(*R*_{2} − *R*_{1}) for the positions *s* = 0.3, *s* = 0.5, and *s* = 0.7. When the crack radius is less than $5%$ of the annular radius, i.e., *ɛ*/(*R*_{2} − *R*_{1}) ≤ 0.05, the difference between *G*^{TD} and *G*^{BE} is small. Once again, the difference between *G*^{TD} and *G*^{BE} increases for larger cracks. We note that the difference between *G*^{TD} and *G*^{BE} increases at different rates for the three positions shown in Fig. 13. For large cracks, the difference between *G*^{TD} and *G*^{BE} is much greater at *s* = 0.7 than it is at *s* = 0.3 or *s* = 0.5.

## 6 Conclusions

An energy release rate approximation for half-penny cracks at any boundary location $x^$, radius *ɛ*, direction **d**, and orientation *α* combination has been developed using the topological derivative. The proposed method only requires a single numerical analysis, in contrast to current boundary element and finite element-based methods which require a numerical analysis for each combination. In addition, the proposed method is evaluated on the non-cracked domain, obviating the need for refined meshes in the crack tip region. Plots of the so-called critical topological derivative values reveal the crack initiation locations $x^$ most susceptible to failure by fracture.

The stress field **T _{ɛ}** was expanded using a stress composite expansion. The topological derivative has been calculated by using this stress composite expansion along with the so-called topological-shape sensitivity method. To this end, the stress intensity factors were also calculated, using the stress composite expansion together with the weight functions for a half-penny crack in a half-space. These weight functions have been evaluated in an off-line computation and tabulated along the crack front for half-penny cracks of any orientation

*α*. The proposed approximation can be applied to surface cracks of arbitrary geometries by calculating the corresponding weight functions. In the same manner, the proposed approximation could be extended to anisotropic materials or embedded cracks.

Numerical examples compare the proposed approximation of the energy release rate *G*^{TD} with reference values *G*^{BE} that were calculated using the boundary element method. In all cases, *G*^{TD} and *G*^{BE} are in good agreement, particularly when the crack radius *ɛ* is small. The computational time required to compute *G*^{TD} is orders of magnitude smaller than the computational time required to compute *G*^{BE}.

## Footnotes

In the equations that follow, the *ɛ*, $x^$, **d**, and *α* dependencies of the domain Ω_{ɛ} and boundaries ∂Ω_{ɛ}, *γ*_{ɛ}, and ∂*γ*_{ɛ} are dropped for conciseness.

In the equations that follow, the *j* = 1 indices are dropped for conciseness.

In the equations that follow, the *s* dependencies of the directions **e**_{n}, **e**_{t}, and **e**_{b} are dropped for conciseness.

## Acknowledgment

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. The authors gratefully acknowledge the support of the National Science Foundation (Award CMI-1200086). The authors also want to recognize insightful comments from the reviewers and from Dr. Mariana Silva, and the many helpful conversations with Tom Curtin of BEASY USA – Computational Mechanics Inc.

## Nomenclature

**d**=crack direction

**e**=unit vector

**m**=outward normal vector on

*γ*_{ɛ}**n**=outward normal vector on ∂Ω

**u**=displacement vector

**x**=position vector

**y**=scaled position vector

*f*=gauge function

*s*=position on

*γ*_{ɛ}or*γ*_{1}- $C$ =
elasticity tensor

*D*=Fréchet derivative

*E*=Young’s modulus

*G*=energy release rate

*J*=*J*-integral**T**=outer stress (Cauchy stress tensor in Ω)

- $H$ =
scaled half-space domain with a crack

- $R$ =
remainder function

- $x^$ =
crack location

- $t~$ =
prescribed boundary traction on

*γ*_{1} - $T~$ =
inner stress (Cauchy stress tensor in $H)$

*h*=_{ij}weight functions

*D*_{T}=topological derivative

*G*_{C}=fracture toughness

*K*_{λ}=mode

*λ*stress intensity factor- $K\xaf\lambda $ =
normalized mode

*λ*stress intensity factor - $T\u02c7\epsilon $ =
remainder stress

**T**=_{ɛ}Cauchy stress tensor in Ω

_{ɛ}**t**^{p}=prescribed boundary traction on ∂Ω

**t**^{R}=prescribed remainder traction on ∂Ω

- $t\epsilon R$ =
prescribed remainder traction on

*γ*_{ɛ} - $\u2207su$ =
infinitesimal strain tensor

### Greek Letters

### Subscripts or Superscripts

*b*=binormal direction on ∂

*γ*_{ɛ}- BE =
reference (“exact”) numerical values

- (
*j*) =*j*th-order *n*=outward normal direction on ∂

*γ*_{ɛ}*r*=outward normal direction at

*s*= 0.5*t*=tangential direction on ∂

*γ*_{ɛ}- TD =
first-order approximation

*z*=tangential direction at

*s*= 0.5*ɛ*=response quantity evaluated in Ω

_{ɛ}*θ*=binormal direction at

*s*= 0.5- $\u22c6$ =
critical parameter