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 GC. 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 [57], 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 [1820] 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. [2426].

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 [3234] 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 [3840], and structural and topology optimization [4143]. 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 tp, 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 Ωε(ε,x^,d,α). The boundary of Ωε(ε,x^,d,α) is given by Ωε(ε,x^,d,α)=Ωγε(ε,x^,d,α) where γε(ε,x^,d,α) is the crack surface with outward normal vector m and crack front γε(ε,x^,d,α).1

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
ψ(Ωε)=12ΩεsuεTεdVΩtpuεdS
(1)
where u is the displacement vector, su:=12(u+uT) is the infinitesimal strain tensor, T=C[su] 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.,
G(Ωε)=ddεψ(Ωε)
(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 Ωɛ:
divTε=0inΩεTεm=0onγεTεn=tponΩ
(3)
The prescribed traction tp is defined such that global equilibrium is satisfied:
Ωtp(x)dS=0Ω(xz)tp(x)dS=0
(4)
where is the cross product and z is an arbitrary location, not necessarily in the domain.
Our objective is to evaluate 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=1nf(j)(ε)DT(j)ψ(x^,d,α,s)+R(f(n)(ε))
(5)
where R is the remainder function, DT(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=Lπε
(6)
where 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
limε0f(k)(ε)f(j)(ε)=0,k>jandlimε0R(f(n)(ε))f(n)(ε)=0
(7)
In this work, we evaluate the topological derivatives using the so-called topological shape sensitivity method where a crack of boundary location 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
DT(j)ψ(x^,d,α,s)=limε0{1f(j)(ε)(ddεψ(Ωε)i=1j1f(i)(ε)DT(i)ψ(x^,d,α,s))}
(8)
We see from Eq. (8) that the gauge functions f(j) must be defined such that DT(j)ψ(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. (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)ψ(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. (5) with respect to ɛ:
ddεψ(Ωε)=j=1nf(j)(ε)DT(j)ψ(x^,d,α,s)+R(f(n)(ε))f(n)(ε)
(9)
Equations (2) and (9) are then combined to express G as a function of the topological derivatives:
G(ε,x^,d,α,s)=j=1nf(j)(ε)DT(j)ψ(x^,d,α,s)R(f(n)(ε))f(n)(ε)
(10)
Equation (10) is significant because as seen shortly, we are able to easily evaluate it for any crack radius ɛ, location 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(ε,x^,d,α,s)=ψ(Ω)+f(1)(ε)DT(1)ψ(x^,d,α,s)
(11)
DT(1)ψ(x^,d,α,s)=limε0{1f(1)(ε)ddεψ(Ωε)}
(12)
GTD(ε,x^,d,α,s)=f(1)(ε)DT(1)ψ(x^,d,α,s)
(13)
which result from Eqs. (5), (8), and (10) for the case of j = 1.2
In order to proceed, we define the local directions en(s), et(s), and eb(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 en such that the shape derivative is defined as [51]
ddεψ(Ωε)=γεΣεnηendS
(14)
where η = η(s) defines the size of the perturbation and Σ is the energy-momentum tensor
Σε=12(suεTε)IsuεTTε
(15)
Equations (14) and (15) define the path-independent J-integral [27]. Alternatively, the J-integral may be evaluated as [53]
J(Σε)=KI2+KII2E¯+KIII22μ
(16)
where Kλ are the mode λ stress intensity factors which describe the displacement and stress fields in the neighborhood of the crack [2], 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. (16), we introduce the following stress composite expansion [54]:
Tε(x)=T(x)+T~(y)+Tˇε(x)
(17)
where T(x) denotes the outer stress, T~(y) denotes the inner stress, and Tˇε(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.,
Tεn=Tn+T~|y=x/εn+Tˇεn=tp,xϵΩ
(18)
and, since by definition eb = m,
Tεeb=Teb+T~|y=x/εeb+Tˇεeb=0,xϵγε
(19)
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,
divT=0inΩTn=tponΩ
(20)
In general, T(x) does not satisfy the traction-free boundary condition on the crack surface γɛ (Eq. (19)).
Since the crack radius ɛ is small, we expand the outer stress T(x) about the crack initiation site x^ in the en direction, i.e., xx^=δen, to obtain
T(x)=T(x^)+δDT(x^)[en]+δ22!D2T(x^)[en,en]+forxϵγε
(21)
where δ=|xx^|. D represents the Fréchet derivative, which is defined as the linear operators
DT(x)[u]=T(x+u)T(x)+o(|u|)D2T(x)[u,v]=DT(x+v)[u]DT(x)[u]+o(|v|)
(22)
We subsequently use the scaled coordinate ξ = δ/ɛ (Fig. 3) to express Eq. (21) as
T(x)=T(x^)+εξDT(x^)[en]+ε2ξ22!D2T(x^)[en,en]+forxϵγε
(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 γ1). There is a non-zero crack face traction on γ1 and zero far-field traction, i.e.,
divT~=0inHT~eb=t~onγ1T~0as|y|
(24)
The inner stress T~(y) annihilates the leading order term of the traction T(x)eb on the crack surface γɛ due to the outer stress. Therefore, using Eq. (23), the prescribed traction t~ is defined as
t~=T(x^)eb,yγ1
(25)
The remainder stress Tˇε(x) is defined on the cracked domain Ωɛ and satisfies the boundary value problem,
divTˇε=0inΩεTˇεeb=tεRonγεTˇεn=tRonΩ
(26)
The remainder stress ensures that Eqs. (18) and (19) are satisfied. To this end, the prescribed tractions tεR and tR are defined as
tεR={T(x)T(x^)}eb,xγεtR=T~(x/ε)n,xΩ
(27)
In the (en, et, eb) coordinate system, the inner stress can be expressed as
T~(ρ,ϕ)=ερKF(ϕ)+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 eb, and F(ϕ) is a dimensionless tensor function of ϕ [55]. Notably, T~(x/ε)=O(ε) far from the crack front, in particular for x ∈ ∂Ω. Furthermore, we observe from Eq. (23) that
T(x)T(x^)=εξDT(x^)[en]+=O(ε)forxϵγε
(29)
Thus, Tˇε=O(ε) and we rewrite Eq. (17) as
Tε(x)=T(x)+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 K¯I, K¯II, and K¯III such that
(K¯I(x^,d,α,s)K¯II(x^,d,α,s)K¯III(x^,d,α,s))1πε(KI(ε,x^,d,α,s)KII(ε,x^,d,α,s)KIII(ε,x^,d,α,s))
(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 er, ez, and eθ which are the outward normal, tangent, and binormal directions at the position s = 0.5 (Fig. 2), i.e.,
{er,ez,eθ}={en,et,eb}|s=0.5
(32)
Using Eq. (25) to relate the prescribed traction t~ to the outer stress T(x), the normalized factors are expressed as
(K¯I(x^,d,α,s)K¯II(x^,d,α,s)K¯III(x^,d,α,s))=(h11(α,s)h12(α,s)h13(α,s)h21(α,s)h22(α,s)h23(α,s)h31(α,s)h32(α,s)h33(α,s))×(Tθθ(x^,d,α)Trθ(x^,d,α)Tzθ(x^,d,α))
(33)
where hij(α, s) are the weight functions, dimensionless functions of the orientation α and location s and, e.g., the Trθ=erTeθ are the crack face traction components t~=T(x^)eθ, cf. Eq. (25).
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 t~=eθ, t~=er, and t~=ez were then solved. These problems were used to calculate hi1(α, s), hi2(α, s), and hi3(α, s) for i = {1, 2, 3}, respectively. For example, when t~=eθ then Trθ(x^,d,α)=0 and Tzθ(x^,d,α)=0, and Eq. (33) reduces to
(K¯I(x^,α,s)K¯II(x^,α,s)K¯III(x^,α,s))=(h11(α,s)h21(α,s)h31(α,s))Tθθ(x^,α)
(34)
The hij(α, s) were calculated for 75degα75deg 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 ±1.13%.
Having the normalized factors, the topological derivative is expressed by combining Eqs. (8), (16), and (31) as
DTψ(x^,d,α,s)=limε0{πεf(j)(ε)(K¯I2(x^,d,α,s)+K¯II2(x^,d,α,s)E¯+K¯III2(x^,d,α,s)2μ)}
(35)
where the gauge function f(ɛ) is defined so that Eq. (35) evaluates to a non-zero finite number cf. Eq. (8)
f(ε)=ε22
(36)
From Eq. (35) and the above definition of f(ɛ), the topological derivative becomes
DTψ(x^,d,α,s)=π{K¯I2(x^,d,α,s)+K¯II2(x^,d,α,s)E¯+K¯III2(x^,d,α,s)2μ}
(37)
Finally, we obtain the energy release rate approximation GTD(ε,x^,d,α,s) by combining Eqs. (13) and (37) to give
GTD(ε,x^,d,α,s)=πε{K¯I2(x^,d,α,s)+K¯II2(x^,d,α,s)E¯+K¯III2(x^,d,α,s)2μ}
(38)

In summary, we perform the following steps to compute the energy release rate approximation GTD(ε,x^,d,α,s):

  1. Perform an analysis on the non-cracked domain Ω with prescribed traction tp to evaluate the outer stress T(x).

  2. Evaluate K¯I(x^,d,α,s), K¯II(x^,d,α,s), and K¯III(x^,d,α,s) via Eq. (33) and the tabulated hij(α, s) values.

  3. Evaluate GTD(ε,x^,d,α,s) via Eq. (38).

Notably, we know the value of GTD(ε,x^,d,α,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 GTD for any crack radius ɛ, location x^, direction d, orientation α, and position s. To verify the proposed approximation, reference energy release rates GBE 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 GTD is in good agreement with GBE.

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=2m, Ri = 5 m, Ro = 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 denotes the (d, α, s) that maximizes GTD for a given ɛ, hereafter referred to as the critical crack parameters, i.e.,
{d(x^),α(x^),s(x^)}=argmax(d,α,s)C{DTψ(x^,d,α,s)}D={d|dR3anddn(x^)=0andd2=1}A={α|αRand65degα65deg}S={s|sRand0.1s0.9}C=D×A×S
(39)
where 2 denotes the Euclidean norm and × denotes the Cartesian product. The ranges of the orientation α 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 sS when cracks of radius ɛ/T = 0.01, direction d(x^), and orientation α(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 ±1.44%. When a crack is introduced at x^1, the stress intensity factors KII and KIII are two orders of magnitude less than KI, i.e., the energy release rate is controlled by KI. When a crack is introduced at x^2, KIITD is within ±1.32% of KIIBE and KIIITD is within ±1.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 GTD with the reference boundary element value GBE. The difference between GTD and GBE does not exceed ±1.67%. Notably, the energy release rates GTD and GBE 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(x^) determined using the topological derivative and boundary element method are the same.

The energy release rate is next calculated for (d,s)(D×S) when cracks of radius ɛ/T = 0.01 and orientation α(x^) are introduced. To visualize these results, we define the angle β(x^,d) which defines d such that
β(x^,d)={arccos(de3)ifx^=x^1arccos(de2)ifx^=x^2
(40)
and the functional Gmax(ε,x^,d,α) which extracts the positional dependence from the energy release rate, i.e.,
Gmax(ε,x^,d,α)=maxsS{G(ε,x^,d,α,s)}
(41)

Figure 8(a) compares GmaxTD and GmaxBE, where it is seen that similar critical crack directions d(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 ±3.84% of GmaxBE outside of the range 100degβ(x^,d)120deg. Within this range, the small values of GmaxTD and GmaxBE result in larger discrepancies.

The energy release rate is then calculated for (α,s)(A×S) when cracks of radius ɛ/T = 0.01 and direction d(x^) are introduced. GmaxTD and GmaxBE are compared in Fig. 8(b), where it is again seen that the critical crack orientations α(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 α40deg, GmaxTD is within ±2.77% of GmaxBE. The difference between GmaxTD and GmaxBE grows as large as ±7.97% when α50deg. One conjecture is that when α50deg, 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 sS when cracks of various radii ɛ are introduced with direction d(x^) and orientation α(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.

The critical topological derivative DTψ(x^), which is the maximum of DTψ for a given crack radius ɛ, i.e.,
DTψ(x^)=DTψ(x^,d(x^),α(x^),s(x^))
(42)
is computed on the boundary ∂Ω by solving an optimization problem of Eq. (39). To visualize these results, we define the three surfaces
Ω1={x^|x^R3andx^Ωand(x^e1)2+(x^e2)2=Ri}Ω2={x^|x^R3andx^Ωandx^e3=T}Ω3={x^|x^R3andx^Ωandx^e3=T}
(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 GTD is linearly proportional to the topological derivative DTψ, cf. Eqs. (13) and (36). Thus, cracks of radius ɛ will propagate at locations x^ if εGC/DTψ(x^), where GC is the fracture toughness of the material. As seen in Fig. 10, there are two regions where DTψ(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 Ω.
We can also evaluate the critical direction d(x^) on ∂Ω. To do this, the definition of β(x^,d) in Eq. (40) is redefined on these boundaries, i.e.,
β(x^,d)={arccos(de3)ifx^Ω1arccos(de2)ifx^Ω2arccos(de2)ifx^Ω3
(44)
On the boundary ∂Ω1, d(x^) is generally between 0deg and 40deg. On the boundaries ∂Ω2 and ∂Ω3, d(x^) varies between the lower bound of 0deg and upper bound of 180deg. For most locations x^, d(x^) is the direction which maximizes the local opening stress component Tθθ(x^,d,α) and mode I stress intensity factor approximation KITD. Contrary to the results for d(x^), there is limited variation in the critical orientation α(x^) on ∂Ω1, ∂Ω2, and ∂Ω3. Excluding the edges of the boundaries, α(x^) is generally between 6deg and 6deg, i.e., nearly perpendicular to the boundaries. For all locations 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 hij(α, 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, R1 = 1 m, R2 = 1.5 m, R3 = 2 m, R4 = 0.075 m, h1 = 1 m, h2 = 0.5 m, and σ = 100 MPa. The stress intensity factors and energy release rate are calculated for sS when cracks of various radii ɛ are introduced at boundary location x^3=1.00e1+0.250e3 with direction d = 1.00e3 and orientation α=45deg. The stress intensity factors versus the position s are shown in Figs. 12(a)12(c) for a crack of radius ɛ = (R2R1)/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 GTD is within ±5.62% of the reference boundary element value GBE. The critical crack positions s(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 ɛ/(R2R1) 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., ɛ/(R2R1) ≤ 0.05, the difference between GTD and GBE is small. Once again, the difference between GTD and GBE increases for larger cracks. We note that the difference between GTD and GBE increases at different rates for the three positions shown in Fig. 13. For large cracks, the difference between GTD and GBE 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 GTD with reference values GBE that were calculated using the boundary element method. In all cases, GTD and GBE are in good agreement, particularly when the crack radius ɛ is small. The computational time required to compute GTD is orders of magnitude smaller than the computational time required to compute GBE.

Footnotes

1

In the equations that follow, the ɛ, 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 en, et, and eb 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)

  •  
  • hij =

    weight functions

  •  
  • DT =

    topological derivative

  •  
  • GC =

    fracture toughness

  •  
  • Kλ =

    mode λ stress intensity factor

  •  
  • K¯λ =

    normalized mode λ stress intensity factor

  •  
  • Tˇε =

    remainder stress

  •  
  • Tɛ =

    Cauchy stress tensor in Ωɛ

  •  
  • tp =

    prescribed boundary traction on ∂Ω

  •  
  • tR =

    prescribed remainder traction on ∂Ω

  •  
  • tεR =

    prescribed remainder traction on γɛ

  •  
  • su =

    infinitesimal strain tensor

Greek Letters

     
  • α =

    crack orientation

  •  
  • γɛ =

    crack surface in Ωɛ

  •  
  • γɛ =

    crack front in Ωɛ

  •  
  • γ1 =

    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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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)
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
53.
Tada
,
H.
,
Paris
,
P. C.
, and
Irwin
,
G. R.
,
2000
,
The Stress Analysis of Cracks Handbook
, 3rd ed.,
American Society of Mechanical Engineers
.
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