Energy Release Rate Approximation for Small Surface Cracks in Three-Dimensional Domains Using the Topological Derivative

Alidoost, Kazem; Feng, Meng; Geubelle, Philippe H.; Tortorelli, Daniel A.

doi:10.1115/1.4045793

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 $\hat{x}$ with radius ɛ, direction d, and orientation α, we obtain the cracked domain $Ω_{ε} (ε, \hat{x}, d, α)$ ⁠. The boundary of $Ω_{ε} (ε, \hat{x}, d, α)$ is given by $\partial Ω_{ε} (ε, \hat{x}, d, α) = \partial Ω \cup γ_{ε} (ε, \hat{x}, d, α)$ where $γ_{ε} (ε, \hat{x}, d, α)$ is the crack surface with outward normal vector m and crack front $\partial γ_{ε} (ε, \hat{x}, d, α)$ ⁠.^¹

Fig. 1

View large Download slide

(a) Cracked domain Ω_ɛ and (b) crack surface γ_ɛ and crack front ∂γ_ɛ

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

ψ (Ω_{ε}) = \frac{1}{2} \int_{Ω_{ε}} \nabla^{s} u_{ε} \cdot T_{ε} dV - \int_{\partial Ω} t^{p} \cdot u_{ε} dS

(1)

where u is the displacement vector,

\nabla^{s} u : = \frac{1}{2} (\nabla u + \nabla u^{T})

is the infinitesimal strain tensor,

T = C [\nabla^{s} u]

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. with respect to the crack radius ɛ [50,51], i.e.,

G (Ω_{ε}) = - \frac{d}{d ε} ψ (Ω_{ε})

(2)

For this traction boundary value problem, we assume that there are no body forces and the crack faces are traction free.

The governing equations of linear elasticity are satisfied by the response quantities on Ω_ɛ:

\begin{aligned} div T_{ε} & = 0 in Ω_{ε} \\ T_{ε} m & = 0 on γ_{ε} \\ T_{ε} n & = t^{p} on \partial Ω \end{aligned}

(3)

The prescribed traction t^p is defined such that global equilibrium is satisfied:

\begin{aligned} \int_{\partial Ω} t^{p} (x) dS = 0 \\ \int_{\partial Ω} (x - z) \land t^{p} (x) dS = 0 \end{aligned}

(4)

where

\land

is the cross product and z is an arbitrary location, not necessarily in the domain.

Our objective is to evaluate G of Eq. for a half-penny crack of any boundary location

\hat{x}

⁠, crack radius ɛ, direction d, and orientation α. To this end, a topological asymptotic expansion is constructed to express the total potential energy as:

ψ (Ω_{ε}) = ψ (Ω) + \sum_{j = 1}^{n} f_{(j)} (ε) D_{T}^{(j)} ψ (\hat{x}, d, α, s) + R (f_{(n)} (ε))

(5)

where

R

is the remainder function,

D_{T}^{(j)} ψ

represents the jth-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

s = \frac{L}{π ε}

(6)

where L is the arc length (Fig. 2). The gauge functions f_(j) in Eq. are functions of the crack radius ɛ and the crack surface boundary conditions. They monotonically tend to zero as ɛ tends to zero and satisfy

lim_{ε \to 0} \frac{f_{(k)} (ε)}{f_{(j)} (ε)} = 0, k > j and lim_{ε \to 0} \frac{R (f_{(n)} (ε))}{f_{(n)} (ε)} = 0

(7)

Fig. 2

View large Download slide

Local directions and position s on the crack front ∂γ_ɛ

In this work, we evaluate the topological derivatives using the so-called topological shape sensitivity method where a crack of boundary location

\hat{x}

⁠, radius ɛ, 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

D_{T}^{(j)} ψ (\hat{x}, d, α, s) = lim_{ε \to 0} {\frac{1}{f_{(j)}^{'} (ε)} (\frac{d}{d ε} ψ (Ω_{ε}) - \sum_{i = 1}^{j - 1} f_{(i)}^{'} (ε) D_{T}^{(i)} ψ (\hat{x}, d, α, s))}

(8)

We see from Eq. that the gauge functions f_(j) must be defined such that

D_{T}^{(j)} ψ (\hat{x}, d, α, 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. is positive. We also see from Eqs. and that if f_(j) is a suitable gauge function, then ζf_(j) is also suitable for any constant ζ > 0 and

D_{T}^{(j)} ψ (\hat{x}, d, α, s)

is unaffected by the choice of ζ.

In order to express the shape sensitivity of ψ as a function of the topological derivatives, we differentiate Eq. with respect to ɛ:

\frac{d}{d ε} ψ (Ω_{ε}) = \sum_{j = 1}^{n} f_{(j)}^{'} (ε) D_{T}^{(j)} ψ (\hat{x}, d, α, s) + R^{'} (f_{(n)} (ε)) f_{(n)}^{'} (ε)

(9)

Equations and are then combined to express G as a function of the topological derivatives:

G (ε, \hat{x}, d, α, s) = - \sum_{j = 1}^{n} f_{(j)}^{'} (ε) D_{T}^{(j)} ψ (\hat{x}, d, α, s) - R^{'} (f_{(n)} (ε)) f_{(n)}^{'} (ε)

(10)

Equation is significant because as seen shortly, we are able to easily evaluate it for any crack radius ɛ, location

\hat{x}

⁠, direction d, and orientation α combination using response quantities computed on the non-cracked domain Ω.

In this work, we restrict ourselves to first-order analysis and define the first-order approximations

ψ^{TD} (ε, \hat{x}, d, α, s) = ψ (Ω) + f_{(1)} (ε) D_{T}^{(1)} ψ (\hat{x}, d, α, s)

(11)

D_{T}^{(1)} ψ (\hat{x}, d, α, s) = lim_{ε \to 0} {\frac{1}{f_{(1)}^{'} (ε)} \frac{d}{d ε} ψ (Ω_{ε})}

(12)

G^{TD} (ε, \hat{x}, d, α, s) = - f_{(1)}^{'} (ε) D_{T}^{(1)} ψ (\hat{x}, d, α, s)

(13)

which result from Eqs. , , and for the case of j = 1.^²

In order to proceed, we define the local directions 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).^³ The direction of crack growth is presumed to be in the normal direction e_n such that the shape derivative is defined as [51]

\frac{d}{d ε} ψ (Ω_{ε}) = - \int_{\partial γ_{ε}} Σ_{ε} n \cdot η e_{n} dS

(14)

where η = η(s) defines the size of the perturbation and Σ is the energy-momentum tensor

Σ_{ε} = \frac{1}{2} (\nabla^{s} u_{ε} \cdot T_{ε}) I - \nabla^{s} u_{ε}^{T} T_{ε}

(15)

Equations and define the path-independent J-integral [27]. Alternatively, the J-integral may be evaluated as [53]

J (Σ_{ε}) = \frac{K_{I}^{2} + K_{I I}^{2}}{\bar{E}} + \frac{K_{I I I}^{2}}{2 μ}

(16)

where K_λ are the mode λ stress intensity factors which describe the displacement and stress fields in the neighborhood of the crack [2],

\bar{E} = E / (1 - ν^{2})

⁠, E is the Young’s modulus, ν is the Poisson’s ratio, and μ is the shear modulus.

3 Asymptotic Analysis

To compute the stress intensity factors in Eq. , we introduce the following stress composite expansion [54]:

T_{ε} (x) = T (x) + \tilde{T} (y) + {\overset{ˇ}{T}}_{ε} (x)

(17)

where T(x) denotes the outer stress,

\tilde{T} (y)

denotes the inner stress, and

{\overset{ˇ}{T}}_{ε} (x)

denotes the remainder stress (Fig. 3). Note that the inner stress

\tilde{T} (y)

is computed at the scaled position vector y = x/ɛ. By construction, the boundary conditions of Eq. are satisfied by the stress composite expansion, i.e.,

T_{ε} n = T n + \tilde{T} |_{y = x / ε} n + {\overset{ˇ}{T}}_{ε} n = t^{p}, x ϵ \partial Ω

(18)

and, since by definition e_b = m,

T_{ε} e_{b} = T e_{b} + \tilde{T} |_{y = x / ε} e_{b} + {\overset{ˇ}{T}}_{ε} e_{b} = 0, x ϵ γ_{ε}

(19)

Fig. 3

View large Download slide

Composite expansion expressed as the sum of the responses on a domain without the crack Ω, scaled half-space domain with a crack $H$ ⁠, and domain with a crack Ω_ɛ subject to remainder tractions

The outer stress T(x) satisfies the nominal problem, i.e., it is defined on the non-cracked domain Ω and satisfies the boundary value problem,

\begin{aligned} div T & = 0 in Ω \\ T n & = t^{p} on \partial Ω \end{aligned}

(20)

In general, T(x) does not satisfy the traction-free boundary condition on the crack surface γ_ɛ (Eq. ).

Since the crack radius ɛ is small, we expand the outer stress T(x) about the crack initiation site

\hat{x}

in the e_n direction, i.e.,

x - \hat{x} = δ e_{n}

⁠, to obtain

T (x) = T (\hat{x}) + δ D T (\hat{x}) [e_{n}] + \frac{δ^{2}}{2!} D^{2} T (\hat{x}) [e_{n}, e_{n}] + \dots for x ϵ γ_{ε}

(21)

where

δ = | x - \hat{x} |

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

\begin{aligned} D T (x) [u] & = T (x + u) - T (x) + o (| u |) \\ D^{2} T (x) [u, v] & = D T (x + v) [u] - D T (x) [u] + o (| v |) \\ ⋮ \end{aligned}

(22)

We subsequently use the scaled coordinate ξ = δ/ɛ (Fig. 3) to express Eq. as

T (x) = T (\hat{x}) + ε ξ D T (\hat{x}) [e_{n}] + ε^{2} \frac{ξ^{2}}{2!} D^{2} T (\hat{x}) [e_{n}, e_{n}] + \dots for x ϵ γ_{ε}

(23)

As the inner stress is evaluated at the scaled position vector y, it satisfies the boundary value problem of a half-space

H

with a unit radius half-penny crack (crack surface γ₁). There is a non-zero crack face traction on γ₁ and zero far-field traction, i.e.,

\begin{aligned} div \tilde{T} & = 0 in H \\ \tilde{T} e_{b} & = \tilde{t} on γ_{1} \\ \tilde{T} & \to 0 as | y | \to \infty \end{aligned}

(24)

The inner stress

\tilde{T} (y)

annihilates the leading order term of the traction T(x)e_b on the crack surface γ_ɛ due to the outer stress. Therefore, using Eq. , the prescribed traction

\tilde{t}

is defined as

\tilde{t} = - T (\hat{x}) e_{b}, y \in γ_{1}

(25)

The remainder stress

{\overset{ˇ}{T}}_{ε} (x)

is defined on the cracked domain Ω_ɛ and satisfies the boundary value problem,

\begin{aligned} div {\overset{ˇ}{T}}_{ε} & = 0 in Ω_{ε} \\ {\overset{ˇ}{T}}_{ε} e_{b} & = t_{ε}^{R} on γ_{ε} \\ {\overset{ˇ}{T}}_{ε} n & = t^{R} on \partial Ω \end{aligned}

(26)

The remainder stress ensures that Eqs. and are satisfied. To this end, the prescribed tractions

t_{ε}^{R}

and t^R are defined as

\begin{aligned} t_{ε}^{R} & = - {T (x) - T (\hat{x})} e_{b}, x \in γ_{ε} \\ t^{R} & = - \tilde{T} (x / ε) n, x \in \partial Ω \end{aligned}

(27)

In the (e_n, e_t, e_b) coordinate system, the inner stress can be expressed as

\tilde{T} (ρ, ϕ) = \frac{\sqrt{ε}}{\sqrt{ρ}} K F (ϕ) + O (ρ^{1 / 2})

(28)

where ρ 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,

\tilde{T} (x / ε) = O (ε)

far from the crack front, in particular for x ∈ ∂Ω. Furthermore, we observe from Eq. that

T (x) - T (\hat{x}) = ε ξ D T (\hat{x}) [e_{n}] + \dots = O (ε) for x ϵ γ_{ε}

(29)

Thus,

{\overset{ˇ}{T}}_{ε} = O (ε)

and we rewrite Eq. as

T_{ε} (x) = T (x) + \tilde{T} (y) + O (ε)

(30)

4 Energy Release Rate Algorithm

In order to extract the crack radius dependence from the stress intensity factors, we define the normalized factors

{\bar{K}}_{I}

⁠,

{\bar{K}}_{I I}

⁠, and

{\bar{K}}_{I I I}

such that

(\begin{matrix} {\bar{K}}_{I} (\hat{x}, d, α, s) \\ {\bar{K}}_{I I} (\hat{x}, d, α, s) \\ {\bar{K}}_{I I I} (\hat{x}, d, α, s) \end{matrix}) \equiv \frac{1}{\sqrt{π ε}} (\begin{matrix} K_{I} (ε, \hat{x}, d, α, s) \\ K_{I I} (ε, \hat{x}, d, α, s) \\ K_{I I I} (ε, \hat{x}, d, α, s) \end{matrix})

(31)

There are no normalized factors associated with the nominal problem, i.e., the outer problem, as this problem does not contain a crack. The normalized factors associated with the inclined unit radius half-penny crack problem in the half-space, i.e., the inner problem, may be expressed in terms of the weight functions for half-penny cracks in a half-space. To this end, we define the directions e_r, e_z, and

e_{θ}

which are the outward normal, tangent, and binormal directions at the position s = 0.5 (Fig. 2), i.e.,

{e_{r}, e_{z}, e_{θ}} = {e_{n}, e_{t}, e_{b}} |_{s = 0.5}

(32)

Using Eq. to relate the prescribed traction

\tilde{t}

to the outer stress T(x), the normalized factors are expressed as

\begin{aligned} (\begin{matrix} {\bar{K}}_{I} (\hat{x}, d, α, s) \\ {\bar{K}}_{I I} (\hat{x}, d, α, s) \\ {\bar{K}}_{I I I} (\hat{x}, d, α, s) \end{matrix}) & = (\begin{matrix} h_{11} (α, s) & h_{12} (α, s) & h_{13} (α, s) \\ h_{21} (α, s) & h_{22} (α, s) & h_{23} (α, s) \\ h_{31} (α, s) & h_{32} (α, s) & h_{33} (α, s) \end{matrix}) \\ \times (\begin{matrix} T_{θ θ} (\hat{x}, d, α) \\ T_{r θ} (\hat{x}, d, α) \\ T_{z θ} (\hat{x}, d, α) \end{matrix}) \end{aligned}

(33)

where h_ij(α, s) are the weight functions, dimensionless functions of the orientation α and location s and, e.g., the

T_{r θ} = e_{r} \cdot T e_{θ}

are the crack face traction components

\tilde{t} = T (\hat{x}) e_{θ}

⁠, cf. Eq. .

The weight functions for a half-penny crack in a half-space were calculated numerically using the boundary element method (Fig. 4). Half-penny cracks were introduced at the center of a face on a 2W × 2W × 2W 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

\tilde{t} = e_{θ}

⁠,

\tilde{t} = e_{r}

⁠, and

\tilde{t} = e_{z}

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

\tilde{t} = e_{θ}

then

T_{r θ} (\hat{x}, d, α) = 0

and

T_{z θ} (\hat{x}, d, α) = 0

⁠, and Eq. reduces to

(\begin{matrix} {\bar{K}}_{I} (\hat{x}, α, s) \\ {\bar{K}}_{I I} (\hat{x}, α, s) \\ {\bar{K}}_{I I I} (\hat{x}, α, s) \end{matrix}) = (\begin{matrix} h_{11} (α, s) \\ h_{21} (α, s) \\ h_{31} (α, s) \end{matrix}) T_{θ θ} (\hat{x}, α)

(34)

The h_ij(α, s) were calculated for

- 75 \deg \leq α \leq 75 \deg

in increments of

1 \deg

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

h_{11} (0 \deg, 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

\pm 1.13 %

⁠.

Fig. 4

View large Download slide

Weight functions h_ij for a half-penny shaped surface crack in a half-space as a function of the orientation α and position s

Fig. 5

View large Download slide

Weight function $h_{11} (0 \deg, s)$ computed by the authors using the boundary element method (line) and computed by Tada et al. using the alternating and finite element methods (symbols) [53]

Having the normalized factors, the topological derivative is expressed by combining Eqs. , , and as

D_{T} ψ (\hat{x}, d, α, s) = lim_{ε \to 0} {\frac{- π ε}{f_{(j)}^{'} (ε)} (\frac{{\bar{K}}_{I}^{2} (\hat{x}, d, α, s) + {\bar{K}}_{I I}^{2} (\hat{x}, d, α, s)}{\bar{E}} + \frac{{\bar{K}}_{I I I}^{2} (\hat{x}, d, α, s)}{2 μ})}

(35)

where the gauge function f(ɛ) is defined so that Eq. evaluates to a non-zero finite number cf. Eq.

f (ε) = \frac{- ε^{2}}{2}

(36)

From Eq. and the above definition of f(ɛ), the topological derivative becomes

D_{T} ψ (\hat{x}, d, α, s) = π {\frac{{\bar{K}}_{I}^{2} (\hat{x}, d, α, s) + {\bar{K}}_{I I}^{2} (\hat{x}, d, α, s)}{\bar{E}} + \frac{{\bar{K}}_{I I I}^{2} (\hat{x}, d, α, s)}{2 μ}}

(37)

Finally, we obtain the energy release rate approximation

G^{TD} (ε, \hat{x}, d, α, s)

by combining Eqs. and to give

G^{TD} (ε, \hat{x}, d, α, s) = π ε {\frac{{\bar{K}}_{I}^{2} (\hat{x}, d, α, s) + {\bar{K}}_{I I}^{2} (\hat{x}, d, α, s)}{\bar{E}} + \frac{{\bar{K}}_{I I I}^{2} (\hat{x}, d, α, s)}{2 μ}}

(38)

In summary, we perform the following steps to compute the energy release rate approximation $G^{TD} (ε, \hat{x}, d, α, s)$ ⁠:

Perform an analysis on the non-cracked domain Ω with prescribed traction t^p to evaluate the outer stress T(x).
Evaluate ${\bar{K}}_{I} (\hat{x}, d, α, s)$ ⁠, ${\bar{K}}_{I I} (\hat{x}, d, α, s)$ ⁠, and ${\bar{K}}_{I I I} (\hat{x}, d, α, s)$ via Eq. (33) and the tabulated h_ij(α, s) values.
Evaluate $G^{TD} (ε, \hat{x}, d, α, s)$ via Eq. (38).

Notably, we know the value of $G^{TD} (ε, \hat{x}, d, α, s)$ for any crack radius ɛ, location $\hat{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 $\hat{x}$ ⁠, direction d, orientation α, and position s. To verify the proposed approximation, reference energy release rates G^BE and stress intensity factors $K_{λ}^{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.

In this example, we study a curved beam subjected to shear and tensile loads (Fig. 6). For this numerical analysis, the parameters are E = 207 GPa, ν = 0.3,

T = 2 m

⁠, 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

{\hat{x}}_{1}

and

{\hat{x}}_{2}

⁠. The superscript

⋆

denotes the (d, α, s) that maximizes G^TD for a given ɛ, hereafter referred to as the critical crack parameters, i.e.,

\begin{aligned} {d^{⋆} (\hat{x}), α^{⋆} (\hat{x}), s^{⋆} (\hat{x})} & = \underset{(d, α, s) \in C}{argmax} {D_{T} ψ (\hat{x}, d, α, s)} \\ D & = {d | d \in R^{3} and d \land n (\hat{x}) = 0 and ∥ d ∥_{2} = 1} \\ A & = {α | α \in R and - 65 \deg \leq α \leq 65 \deg} \\ S & = {s | s \in R and 0.1 \leq s \leq 0.9} \\ C & = D \times A \times S \end{aligned}

(39)

where

∥ ∥_{2}

denotes the Euclidean norm and × denotes the Cartesian product. The ranges of the orientation α and position s in Eq. are chosen to reduce free surface effects.

Fig. 6

View large Download slide

Geometry and loading for a curved beam with shear and tensile loads

The stress intensity factors and energy release rate are first calculated for $s \in S$ when cracks of radius ɛ/T = 0.01, direction $d^{⋆} (\hat{x})$ ⁠, and orientation $α^{⋆} (\hat{x})$ are introduced. Figures 7(a)–7(c) compare the stress intensity factor approximations $K_{I}^{TD}$ ⁠, $K_{I I}^{TD}$ ⁠, and $K_{I I I}^{TD}$ with the reference boundary element values $K_{I}^{BE}$ ⁠, $K_{I I}^{BE}$ ⁠, and $K_{I I I}^{BE}$ ⁠. The difference between $K_{I}^{TD}$ and $K_{I}^{BE}$ does not exceed $\pm 1.44 %$ ⁠. When a crack is introduced at ${\hat{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 ${\hat{x}}_{2}$ ⁠, $K_{I I}^{TD}$ is within $\pm 1.32 %$ of $K_{I I}^{BE}$ and $K_{I I I}^{TD}$ is within $\pm 1.26 %$ of $K_{I I I}^{BE}$ ⁠. This excludes positions near s = 0.5 where $K_{I I}^{TD}$ and $K_{I I}^{BE}$ 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 $\pm 1.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^{⋆} (\hat{x})$ determined using the topological derivative and boundary element method are the same.

Fig. 7

View large Download slide

(a)–(c) Stress intensity factors and (d) energy release rate versus position s. The star symbols in (d) denote the maxima of G^TD.

The energy release rate is next calculated for

(d, s) \in (D \times S)

when cracks of radius ɛ/T = 0.01 and orientation

α^{⋆} (\hat{x})

are introduced. To visualize these results, we define the angle

β (\hat{x}, d)

which defines d such that

β (\hat{x}, d) = {\begin{matrix} \arccos (d \cdot e_{3}) & if \hat{x} = {\hat{x}}_{1} \\ \arccos (d \cdot e_{2}) & if \hat{x} = {\hat{x}}_{2} \end{matrix}

(40)

and the functional

G_{max} (ε, \hat{x}, d, α)

which extracts the positional dependence from the energy release rate, i.e.,

G_{max} (ε, \hat{x}, d, α) = max_{s \in S} {G (ε, \hat{x}, d, α, s)}

(41)

Figure 8(a) compares $G_{max}^{TD}$ and $G_{max}^{BE}$ ⁠, where it is seen that similar critical crack directions $d^{⋆} (\hat{x})$ are determined using the topological derivative and boundary element methods. One can also see that $G_{max}^{TD}$ is in good agreement with $G_{max}^{BE}$ ⁠. $G_{max}^{TD}$ is within $\pm 3.84 %$ of $G_{max}^{BE}$ outside of the range $100 \deg \leq β (\hat{x}, d) \leq 120 \deg$ ⁠. Within this range, the small values of $G_{max}^{TD}$ and $G_{max}^{BE}$ result in larger discrepancies.

Fig. 8

Energy release rate versus (a) angle β(x^,d) and (b) orientation α. The star symbols denote the maxima of GmaxTD.

View large Download slide

Energy release rate versus (a) angle $β (\hat{x}, d)$ and (b) orientation α. The star symbols denote the maxima of $G_{max}^{TD}$ ⁠.

The energy release rate is then calculated for $(α, s) \in (A \times S)$ when cracks of radius ɛ/T = 0.01 and direction $d^{⋆} (\hat{x})$ are introduced. $G_{max}^{TD}$ and $G_{max}^{BE}$ are compared in Fig. 8(b), where it is again seen that the critical crack orientations $α^{⋆} (\hat{x})$ determined using the topological derivative and boundary element methods are similar. Once more one can see that $G_{max}^{TD}$ is in good agreement with $G_{max}^{BE}$ ⁠. When $α \leq 40 \deg$ ⁠, $G_{max}^{TD}$ is within $\pm 2.77 %$ of $G_{max}^{BE}$ ⁠. The difference between $G_{max}^{TD}$ and $G_{max}^{BE}$ grows as large as $\pm 7.97 %$ when $α \geq 50 \deg$ ⁠. One conjecture is that when $α \geq 50 \deg$ ⁠, 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 \in S$ when cracks of various radii ɛ are introduced with direction $d^{⋆} (\hat{x})$ and orientation $α^{⋆} (\hat{x})$ ⁠. The energy release rate calculations $G_{max}^{TD}$ and $G_{max}^{BE}$ are compared in Fig. 9. When the crack radius is less than $5 %$ of the thickness, i.e., ɛ/T ≤ 0.05, the difference between $G_{max}^{TD}$ and $G_{max}^{BE}$ is small. As expected, for larger cracks, the difference between $G_{max}^{TD}$ and $G_{max}^{BE}$ increases. However, $G_{max}^{TD}$ remains a reasonable approximation of $G_{max}^{BE}$ even for cracks of considerable radius, e.g., ɛ/T = 0.25.

As indicated in Fig. 9, the energy release rate approximation $G_{max}^{TD}$ associated with the boundary location ${\hat{x}}_{1}$ is more accurate for longer crack radii than the same approximation associated with the boundary location ${\hat{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 $\hat{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 $\hat{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 $G_{max}^{BE}$ 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 $\hat{x}$ and crack radius ɛ. In contrast, the energy release rate approximation $G_{max}^{TD}$ 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 $G_{max}^{TD}$ 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.

Fig. 9

View large Download slide

Energy release rate versus crack radius ɛ/T

The critical topological derivative

D_{T} ψ^{⋆} (\hat{x})

⁠, which is the maximum of D_Tψ for a given crack radius ɛ, i.e.,

D_{T}^{⋆} ψ (\hat{x}) = D_{T} ψ (\hat{x}, d^{⋆} (\hat{x}), α^{⋆} (\hat{x}), s^{⋆} (\hat{x}))

(42)

is computed on the boundary ∂Ω by solving an optimization problem of Eq. . To visualize these results, we define the three surfaces

\begin{aligned} \partial Ω_{1} & = {\hat{x} | \hat{x} \in R^{3} and \hat{x} \in \partial Ω and (\hat{x} \cdot e_{1})^{2} + (\hat{x} \cdot e_{2})^{2} = R_{i}} \\ \partial Ω_{2} & = {\hat{x} | \hat{x} \in R^{3} and \hat{x} \in \partial Ω and \hat{x} \cdot e_{3} = T} \\ \partial Ω_{3} & = {\hat{x} | \hat{x} \in R^{3} and \hat{x} \in \partial Ω and \hat{x} \cdot e_{3} = - T} \end{aligned}

(43)

Such plots can be used to detect the failure by fracture just as plots of the von Mises stress are used to predict failure by yielding. Indeed, the energy release rate approximation G^TD is linearly proportional to the topological derivative D_Tψ, cf. Eqs. and . Thus, cracks of radius ɛ will propagate at locations

\hat{x}

if

ε \geq G_{C} / D_{T} ψ^{⋆} (\hat{x})

⁠, where G_C is the fracture toughness of the material. As seen in Fig. 10, there are two regions where

D_{T} ψ^{⋆} (\hat{x})

is especially large, located on the ∂Ω₁ − ∂Ω₂ and ∂Ω₁ − ∂Ω₃ 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 Ω.

Fig. 10

Surface plots of the critical topological derivative DTψ⋆(x^) on the boundaries (a) ∂Ω1, (b) ∂Ω1 (e2 − e3 plane view), (c) ∂Ω2, and (d) ∂Ω3. Note the same color bars are used for each surface plot.

View large Download slide

Surface plots of the critical topological derivative $D_{T} ψ^{⋆} (\hat{x})$ on the boundaries (a) ∂Ω₁, (b) ∂Ω₁ (e₂ − e₃ plane view), (c) ∂Ω₂, and (d) ∂Ω₃. Note the same color bars are used for each surface plot.

We can also evaluate the critical direction

d^{⋆} (\hat{x})

on ∂Ω. To do this, the definition of

β (\hat{x}, d)

in Eq. is redefined on these boundaries, i.e.,

β (\hat{x}, d) = {\begin{matrix} \arccos (d \cdot e_{3}) & if \hat{x} \in \partial Ω_{1} \\ \arccos (d \cdot e_{2}) & if \hat{x} \in \partial Ω_{2} \\ \arccos (d \cdot e_{2}) & if \hat{x} \in \partial Ω_{3} \end{matrix}

(44)

On the boundary ∂Ω₁,

d^{⋆} (\hat{x})

is generally between

0 \deg

and

40 \deg

⁠. On the boundaries ∂Ω₂ and ∂Ω₃,

d^{⋆} (\hat{x})

varies between the lower bound of

0 \deg

and upper bound of

180 \deg

⁠. For most locations

\hat{x}

⁠,

d^{⋆} (\hat{x})

is the direction which maximizes the local opening stress component

T_{θ θ} (\hat{x}, d, α)

and mode I stress intensity factor approximation

K_{I}^{T D}

⁠. Contrary to the results for

d^{⋆} (\hat{x})

⁠, there is limited variation in the critical orientation

α^{⋆} (\hat{x})

on ∂Ω₁, ∂Ω₂, and ∂Ω₃. Excluding the edges of the boundaries,

α^{⋆} (\hat{x})

is generally between

- 6 \deg

and

6 \deg

⁠, i.e., nearly perpendicular to the boundaries. For all locations

\hat{x}

⁠, the critical position

s^{⋆}

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 m, R₂ = 1.5 m, R₃ = 2 m, R₄ = 0.075 m, h₁ = 1 m, h₂ = 0.5 m, and σ = 100 MPa. The stress intensity factors and energy release rate are calculated for $s \in S$ when cracks of various radii ɛ are introduced at boundary location ${\hat{x}}_{3} = 1.00 e_{1} + 0.250 e_{3}$ with direction d = 1.00e₃ and orientation $α = 45 \deg$ ⁠. The stress intensity factors versus the position s are shown in Figs. 12(a)–12(c) for a crack of radius ɛ = (R₂ − R₁)/100. As in the previous example, one can see the stress intensity factors $K_{I}^{T D}$ ⁠, $K_{I I}^{T D}$ ⁠, and $K_{I I I}^{T D}$ computed using the topological derivative are in good agreement with the reference boundary element values $K_{I}^{B E}$ ⁠, $K_{I I}^{B E}$ ⁠, and $K_{I I I}^{B E}$ ⁠. The energy release rates appear in Fig. 12(d), where it is seen that the proposed approximation G^TD is within $\pm 5.62 %$ of the reference boundary element value G^BE. The critical crack positions $s^{⋆} (\hat{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₂ − R₁) 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₂ − R₁) ≤ 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.

Fig. 11

View large Download slide

Geometry and loading for a coupling with axial and radial loads

Fig. 12

View large Download slide

(a)–(c) Stress intensity factors and (d) energy release rate versus position s. The star symbol in (d) denotes the maximum of G^TD.

Fig. 13

Energy release rate versus crack radius ɛ/(R2 − R1)

View large Download slide

Energy release rate versus crack radius ɛ/(R₂ − R₁)

6 Conclusions

An energy release rate approximation for half-penny cracks at any boundary location $\hat{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 $\hat{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

1

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

2

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

3

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 γ₁
$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
$\hat{x}$ =: crack location
$\tilde{t}$ =: prescribed boundary traction on γ₁
$\tilde{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
${\bar{K}}_{λ}$ =: normalized mode λ stress intensity factor
${\overset{ˇ}{T}}_{ε}$ =: remainder stress
T_ɛ =: Cauchy stress tensor in Ω_ɛ
t^p =: prescribed boundary traction on ∂Ω
t^R =: prescribed remainder traction on ∂Ω
$t_{ε}^{R}$ =: prescribed remainder traction on γ_ɛ
$\nabla^{s} u$ =: infinitesimal strain tensor

Greek Letters

α =: crack orientation
γ_ɛ =: crack surface in Ω_ɛ
∂γ_ɛ =: crack front in Ω_ɛ
γ₁ =: crack surface in $H$
ɛ =: crack radius
μ =: shear modulus
ν =: Poisson’s ratio
ψ =: total potential energy
Σ =: energy momentum tensor
Ω =: non-cracked domain
Ω_ɛ =: cracked domain
∂Ω =: boundary of Ω
∂Ω_ɛ =: boundary of Ω_ɛ

Subscripts or Superscripts

b =: binormal direction on ∂γ_ɛ
BE =: reference (“exact”) numerical values
(j) =: jth-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
$⋆$ =: critical parameter

References

1.

Griffith

,

A. A.

,

1921

, “

The Phenomena of Rupture and Flow in Solids

,”

Philos. Trans. R. Soc. London. Ser. A, Containing Papers Math. Phys. Charact.

,

221

, pp.

163

–

198

. 10.1098/rsta.1921.0006

Google Scholar

Crossref

2.

Williams

,

M. L.

,

1957

, “

On the Stress Distribution At the Base of a Stationary Crack

,”

ASME J. Appl. Mech.

,

24

(

1

), pp.

109

–

114

.

3.

Palani

,

G.

,

Iyer

,

N. R.

, and

Dattaguru

,

B.

,

2006

, “

New a Posteriori Error Estimator and Adaptive Mesh Refinement Strategy for 2-D Crack Problems

,”

Eng. Fract. Mech.

,

73

(

6

), pp.

802

–

819

. 10.1016/j.engfracmech.2005.10.003

Google Scholar

Crossref

4.

Souiyah

,

M.

,

Alshoaibi

,

A.

,

Muchtar

,

A.

, and

Ariffin

,

A.

,

2009

, “

Two-Dimensional Finite Element Method for Stress Intensity Factor Using Adaptive Mesh Strategy

,”

Acta Mech.

,

204

(

1–2

), pp.

99

–

108

. 10.1007/s00707-008-0054-2

Google Scholar

Crossref

5.

Ayhan

,

A. O.

,

2007

, “

Stress Intensity Factors for Three-Dimensional Surface Cracks Using Enriched Finite Elements

,”

Int. J. Solids Struct.

,

44

(

25-26

), pp.

8579

–

8599

. 10.1016/j.ijsolstr.2007.06.022

Google Scholar

Crossref

6.

Munaswamy

,

K.

, and

Pullela

,

R.

,

2007

, “

Computation of Stress Intensity Factors for Through Cracks in Plates Using P-Version Finite Element Method

,”

Commun. Numer. Methods Eng.

,

24

(

12

), pp.

1753

–

1780

. 10.1002/cnm.1065

Google Scholar

Crossref

7.

Treifi

,

M.

,

Oyadiji

,

S. O.

, and

Tsang

,

D. K.

,

2008

, “

Computations of Modes I and II Stress Intensity Factors of Sharp Notched Plates Under In-Plane Shear and Bending Loading by the Fractal-Like Finite Element Method

,”

Int. J. Solids Struct.

,

45

(

25–26

), pp.

6468

–

6484

. 10.1016/j.ijsolstr.2008.08.013

Google Scholar

Crossref

8.

Ju

,

S.

,

1998

, “

Simulating Three-Dimensional Stress Intensity Factors by the Least-Squares Method

,”

Int. J. Numer. Methods Eng.

,

43

(

8

), pp.

1437

–

1451

. 10.1002/(SICI)1097-0207(19981230)43:8<1437::AID-NME477>3.0.CO;2-9

Google Scholar

Crossref

9.

Yu

,

T.

, and

Shi

,

L.

,

2012

, “

Determination of Sharp V-Notch Stress Intensity Factors Using the Extended Finite Element Method

,”

J. Strain Anal. Eng. Des.

,

47

(

2

), pp.

95

–

103

. 10.1177/0309324711433981

Google Scholar

Crossref

10.

Fehl

,

B.

, and

Truman

,

K.

,

1999

, “

An Evaluation of Fracture Mechanics Quarter-Point Displacement Techniques Used for Computing Stress Intensity Factors

,”

Eng. Struct.

,

21

(

5

), pp.

406

–

415

. 10.1016/S0141-0296(97)00221-6

Google Scholar

Crossref

11.

Yan

,

X.

,

2005

, “

Stress Intensity Factors for Asymmetric Branched Cracks in Plane Extension by Using Crack-Tip Displacement Discontinuity Elements

,”

Mech. Res. Commun.

,

32

(

4

), pp.

375

–

384

. 10.1016/j.mechrescom.2004.10.005

Google Scholar

Crossref

12.

Diamantoudis

,

A. T.

, and

Labeas

,

G.

,

2005

, “

Stress Intensity Factors of Semi-Elliptical Surface Cracks in Pressure Vessels by Global-Local Finite Element Methodology

,”

Eng. Fract. Mech.

,

72

(

9

), pp.

1299

–

1312

. 10.1016/j.engfracmech.2004.10.004

Google Scholar

Crossref

13.

Si

,

Y.

,

Yongjun

,

X.

, and

Williams

,

F.

,

2007

, “

Computation of Stress Intensity Factors by the Subregion Mixed Finite Element Method of Lines

,”

Acta Mech. Solida Sin.

,

20

(

2

), pp.

149

–

162

. 10.1007/s10338-007-0718-9

Google Scholar

Crossref

14.

Wen

,

P.

, and

Aliabadi

,

M.

,

2011

, “

A Variational Approach for Evaluation of Stress Intensity Factors Using the Element Free Galerkin Method

,”

Int. J. Solids Struct.

,

48

(

7–8

), pp.

1171

–

1179

. 10.1016/j.ijsolstr.2011.01.002

Google Scholar

Crossref

15.

Muthu

,

N.

,

Maiti

,

S.

,

Yan

,

W.

, and

Falzon

,

B.

,

2017

, “

Modelling Interacting Cracks Through a Level Set Using the Element-Free Galerkin Method

,”

Int. J. Mech. Sci.

,

134

, pp.

203

–

215

. 10.1016/j.ijmecsci.2017.10.009

Google Scholar

Crossref

16.

Junior

,

A. M.

,

Sales

,

R.

,

Osterno

,

F.

,

Timbo

,

K.

, and

de Deus

,

E.

,

2018

, “

Failure Analysis of E-Glass Reinforced PMMA Using XFEM

,”

Sci. Technol. Mater.

,

30

(

1

), pp.

34

–

50

. 10.1016/j.stmat.2018.11.002

Google Scholar

Crossref

17.

Tian

,

H.

,

Wen

,

L.

, and

Wang

,

L.

,

2019

, “

Three-Dimensional Improved XFEM (IXFEM) for Static Crack Problems

,”

Comput. Methods Appl. Mech. Eng.

,

343

, pp.

339

–

367

. 10.1016/j.cma.2018.08.029

Google Scholar

Crossref

18.

Mi

,

Y.

, and

Aliabadi

,

M.

,

1992

, “

Dual Boundary Element Method for Three-Dimensional Fracture Mechanics Analysis

,”

Eng. Anal. Boundary Elem.

,

10

(

2

), pp.

161

–

171

. 10.1016/0955-7997(92)90047-B

Google Scholar

Crossref

19.

Portela

,

A.

,

Aliabadi

,

M.

, and

Rooke

,

D.

,

1992

, “

The Dual Boundary Element Method: Effective Implementation for Crack Problems

,”

Int. J. Numer. Methods Eng.

,

33

(

6

), pp.

1269

–

1287

. 10.1002/nme.1620330611

Google Scholar

Crossref

20.

Dong

,

Y.

,

Wang

,

Z.

, and

Wang

,

B.

,

1997

, “

On the Computation of Stress Intensity Factors for Interfacial Cracks Using Quarter-Point Boundary Elements

,”

Eng. Fract. Mech.

,

57

(

4

), pp.

335

–

342

. 10.1016/S0013-7944(97)00057-X

Google Scholar

Crossref

21.

Moftakhar

,

A.

, and

Glinka

,

G.

,

1992

, “

Calculation of Stress Intensity Factors by Efficient Integration of Weight Functions

,”

Eng. Fract. Mech.

,

43

(

5

), pp.

749

–

756

. 10.1016/0013-7944(92)90005-Y

Google Scholar

Crossref

22.

Mikkola

,

T. P.

,

1993

, “

Applications of the Weight Function Method

,”

Eng. Fract. Mech.

,

45

(

2

), pp.

209

–

231

. 10.1016/0013-7944(93)90187-W

Google Scholar

Crossref

23.

Glinka

,

G.

, and

Shen

,

G.

,

1991

, “

Universal Features of Weight Functions for Cracks in Mode I

,”

Eng. Fract. Mech.

,

40

(

6

), pp.

1135

–

1146

. 10.1016/0013-7944(91)90177-3

Google Scholar

Crossref

24.

Kobayashi

,

A. S.

,

1984

, “

Numerical Analysis in Fracture Mechanics

,”

Application of Fracture Mechanics to Materials and Structures

,

G. C.

Sih

,

E.

Sommer

and

W.

Dahl

, eds., pp.

27

–

56

. 10.1007/978-94-009-6146-3_2

25.

Schmitt

,

W.

,

1987

,

Numerical Methods in Fracture Mechanics

,

Springer

,

Netherlands, Dordrecht

, pp.

47

–

74

.

26.

Mohammadnejad

,

M.

,

Liu

,

H.

,

Chan

,

A.

,

Dehkhoda

,

S.

, and

Fukuda

,

D.

,

2018

, “

An Overview on Advances in Computational Fracture Mechanics of Rock

,”

Geosyst. Eng.

, pp.

1

–

24

. 10.1080/12269328.2018.1448006

27.

Rice

,

J. R.

,

1968

, “

A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks

,”

ASME J. Appl. Mech.

,

35

(

2

), pp.

379

–

386

. 10.1115/1.3601206

Google Scholar

Crossref

28.

Shih

,

C.

,

1981

, “

Relationships Between the J-Integral and the Crack Opening Displacement for Stationary and Extending Cracks

,”

J. Mech. Phys. Solids

,

29

(

4

), pp.

305

–

326

. 10.1016/0022-5096(81)90003-X

Google Scholar

Crossref

29.

deLorenzi

,

H.

,

1982

, “

On the Energy Release Rate and the J-Integral for 3-D Crack Configurations

,”

Int. J. Fract.

,

19

(

3

), pp.

183

–

193

. 10.1007/BF00017129

Google Scholar

Crossref

30.

Gosz

,

M.

, and

Moran

,

B.

,

2002

, “

An Interaction Energy Integral Method for Computation of Mixed-Mode Stress Intensity Factors Along Non-Planar Crack Fronts in Three Dimensions

,”

Eng. Fract. Mech.

,

69

(

3

), pp.

299

–

319

. 10.1016/S0013-7944(01)00080-7

Google Scholar

Crossref

31.

Walters

,

M.

,

Paulino

,

G.

, and

Dodds

,

R.

, Jr.,

2006

, “

Computation of Mixed-Mode Stress Intensity Factors for Cracks in Three-Dimensional Functionally Graded Solids

,”

J. Eng. Mech.

,

132

(

1

), pp.

1

–

15

. 10.1061/(ASCE)0733-9399(2006)132:1(1)

Google Scholar

Crossref

32.

deLorenzi

,

H.

,

1985

, “

Energy Release Rate Calculations by the Finite Element Method

,”

Eng. Fract. Mech.

,

21

(

1

), pp.

129

–

143

. 10.1016/0013-7944(85)90060-8

Google Scholar

Crossref

33.

Lin

,

S.-C.

, and

Abel

,

J.

,

1988

, “

Variational Approach for a New Direct-Integration Form of the Virtual Crack Extension Method

,”

Int. J. Fract.

,

38

(

3

), pp.

217

–

235

.

34.

Hwang

,

C.

,

Wawrzynek

,

P.

,

Tayebi

,

A.

, and

Ingraffea

,

A.

,

1998

, “

On the Virtual Crack Extension Method for Calculation of the Rates of Energy Release Rate

,”

Eng. Fract. Mech.

,

59

(

4

), pp.

521

–

542

. 10.1016/S0013-7944(97)00103-3

Google Scholar

Crossref

35.

Sokolowski

,

J.

, and

Zochowski

,

A.

,

1999

, “

On the Topological Derivative in Shape Optimization

,”

SIAM J. Control Optim.

,

37

(

4

), pp.

1251

–

1272

. 10.1137/S0363012997323230

Google Scholar

Crossref

36.

Auroux

,

D.

, and

Masmoudi

,

M.

,

2006

, “

A One-Shot Inpainting Algorithm Based on the Topological Asymptotic Analysis

,”

J. Comput. Appl. Math.

,

25

(

2–3

), pp.

251

–

267

.

37.

Larrabide

,

I.

,

Feijoo

,

R. A.

,

Novotny

,

A. A.

, and

Taroco

,

E. A.

,

2008

, “

Topological Derivative: A Tool for Image Processing

,”

Comput. Struct.

,

86

(

13–14

), pp.

1386

–

1403

. 10.1016/j.compstruc.2007.05.004

Google Scholar

Crossref

38.

Feijoo

,

G. R.

,

2004

, “

A New Method in Inverse Scattering Based on the Topological Derivative

,”

Inverse Prob.

,

20

(

6

), pp.

1819

–

1840

. 10.1088/0266-5611/20/6/008

Google Scholar

Crossref

39.

Bonnet

,

M.

, and

Constantinescu

,

A.

,

2005

, “

Inverse Problems in Elasticity

,”

Inverse Prob.

,

21

(

2

), pp.

R1

–

R50

. 10.1088/0266-5611/21/2/R01

Google Scholar

Crossref

40.

Masmoudi

,

M.

,

Pommier

,

J.

, and

Samet

,

B.

,

2005

, “

The Topological Asymptotic Expansion for the Maxwell Equations and Some Applications

,”

Inverse Prob.

,

21

(

2

), pp.

547

–

564

. 10.1088/0266-5611/21/2/008

Google Scholar

Crossref

41.

Cea

,

J.

,

Garreau

,

S.

,

Guillaume

,

P.

, and

Masmoudi

,

M.

,

2000

, “

The Shape and Topological Optimization Connection

,”

Comput. Methods Appl. Mech. Eng.

,

188

(

4

), pp.

713

–

726

. 10.1016/S0045-7825(99)00357-6

Google Scholar

Crossref

42.

Novotny

,

A. A.

,

Feijoo

,

R. A.

,

Taroco

,

E.

, and

Padra

,

C.

,

2003

, “

Topological Sensitivity Analysis

,”

Comput. Methods Appl. Mech. Eng.

,

192

(

7–8

), pp.

803

–

829

. 10.1016/S0045-7825(02)00599-6

Google Scholar

Crossref

43.

Norato

,

J. A.

,

Bendsoe

,

M. P.

,

Haber

,

R. B.

, and

Tortorelli

,

D. A.

,

2007

, “

A Topological Derivative Method for Topology Optimization

,”

Struct. Multidiscip. Optim.

,

44

(

4-5

), pp.

375

–

386

. 10.1007/s00158-007-0094-6

Google Scholar

Crossref

44.

Amstutz

,

S.

,

Horchani

,

I.

, and

Masmoudi

,

M.

,

2005

, “

Crack Detection by the Topological Gradient Method

,”

Control Cybern.

,

34

(

1

), pp.

81

–

101

.

45.

Khludnev

,

A. M.

,

Novotny

,

A. A.

,

Sokolowski

,

J.

, and

Zochowski

,

A.

,

2009

, “

Shape and Topology Sensitivity Analysis for Cracks in Elastic Bodies on Boundaries of Rigid Inclusions

,”

J. Mech. Phys. Solids

,

57

(

10

), pp.

1718

–

1732

. 10.1016/j.jmps.2009.07.003

Google Scholar

Crossref

46.

Goethem

,

N. V.

, and

Novotny

,

A. A.

,

2010

, “

Crack Nucleation Sensitivity Analysis

,”

Math. Methods Appl. Sci.

,

33

(

16

), pp.

1978

–

1994

.

47.

Alidoost

,

K.

,

Geubelle

,

P. H.

, and

Tortorelli

,

D. A.

,

2018

, “

Energy Release Rate Approximation for Edge Cracks Using Higher-Order Topological Derivatives

,”

Int. J. Fract.

,

210

(

1–2

), pp.

187

–

205

. 10.1007/s10704-018-0271-1

Google Scholar

Crossref

48.

Silva

,

M.

,

Geubelle

,

P. H.

, and

Tortorelli

,

D. A.

,

2011

, “

Energy Release Rate Approximation for Small Surface-Breaking Cracks Using the Topological Derivative

,”

J. Mech. Phys. Solids

,

59

(

5

), pp.

925

–

939

. 10.1016/j.jmps.2011.03.005

Google Scholar

Crossref

49.

Garreau

,

S.

,

Guillaume

,

P.

, and

Masmoudi

,

M.

,

2001

, “

The Topological Asymptotic for PDE Systems: The Elasticity Case

,”

SIAM J. Control Optim.

,

39

(

6

), pp.

1756

–

1778

. 10.1137/S0363012900369538

Google Scholar

Crossref

50.

Feijoo

,

R. A.

,

Padra

,

C.

,

Saliba

,

R.

,

Taroco

,

E.

, and

Venere

,

M. J.

,

2000

, “

Shape Sensitivity Analysis for Energy Release Rate Evaluation and Its Application to the Study of Three-Dimensional Cracked Bodies

,”

Comput. Methods Appl. Mech. Eng.

,

188

(

4

), pp.

649

–

664

. 10.1016/S0045-7825(99)00353-9

Google Scholar

Crossref

51.

Taroco

,

E.

,

2000

, “

Shape Sensitivity Analysis in Linear Elastic Fracture Mechanics

,”

Comput. Methods Appl. Mech. Eng.

,

188

(

4

), pp.

692

–

712

. 10.1016/S0045-7825(99)00356-4

Google Scholar

Crossref

52.

de Faria

,

J. R.

,

Novotny

,

A. A.

,

Feijoo

,

R. A.

,

Taroco

,

E.

, and

Padra

,

C.

,

2007

, “

Second Order Topological Sensitivity Analysis

,”

Int. J. Solids Struct.

,

44

(

14–15

), pp.

4958

–

4977

. 10.1016/j.ijsolstr.2006.12.013

Google Scholar

Crossref

53.

Tada

,

H.

,

Paris

,

P. C.

, and

Irwin

,

G. R.

,

2000

,

The Stress Analysis of Cracks Handbook

, 3rd ed.,

American Society of Mechanical Engineers

.

Google Scholar

Crossref

54.

Kozlov

,

V.

,

Maz’ya

,

V.

, and

Movchan

,

A.

,

1999

,

Asymptotic Analysis of Fields in MultiStructures

,

Oxford University Press

,

Oxford

.

55.

Kassir

,

M. K.

, and

Sih

,

G. C.

,

1966

, “

Three-Dimensional Stress Distribution Around An Elliptical Crack Under Arbitrary Loadings

,”

ASME J. Appl. Mech.

,

33

(

3

), pp.

601

–

611

. 10.1115/1.3625127

Google Scholar

Crossref

Energy Release Rate Approximation for Small Surface Cracks in Three-Dimensional Domains Using the Topological Derivative

Abstract

1 Introduction

2 Topological Derivative

3 Asymptotic Analysis

4 Energy Release Rate Algorithm

5 Examples

5.1 Curved Beam Subjected to Shear and Tensile Loads.

5.2 Coupling Subjected to Axial and Radial Loads.

6 Conclusions

Footnotes

Acknowledgment

Nomenclature

Greek Letters

Subscripts or Superscripts

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Energy Release Rate Approximation for Small Surface Cracks in Three-Dimensional Domains Using the Topological Derivative

Abstract

1 Introduction

2 Topological Derivative

3 Asymptotic Analysis

4 Energy Release Rate Algorithm

5 Examples

5.1 Curved Beam Subjected to Shear and Tensile Loads.

5.2 Coupling Subjected to Axial and Radial Loads.

6 Conclusions

Footnotes

Acknowledgment

Nomenclature

Greek Letters

Subscripts or Superscripts

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only