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.
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 . 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 . 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 . Comparisons among some of these methods are made in Refs. [24–26].
Introduced by Rice , 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 . 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. , 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. . 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. ; in Ref. , 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. , 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 , rather than, e.g., the domain truncation method , 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 with radius ɛ, direction d, and orientation α, we obtain the cracked domain . The boundary of is given by where is the crack surface with outward normal vector m and crack front .1
3 Asymptotic Analysis
4 Energy Release Rate Algorithm
In summary, we perform the following steps to compute the energy release rate approximation :
Notably, we know the value of for any crack radius ɛ, location , direction d, orientation α, and position s from a single analysis on the non-cracked domain Ω.
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 , direction d, orientation α, and position s. To verify the proposed approximation, reference energy release rates GBE and stress intensity factors 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.
The stress intensity factors and energy release rate are first calculated for when cracks of radius ɛ/T = 0.01, direction , and orientation are introduced. Figures 7(a)–7(c) compare the stress intensity factor approximations , , and with the reference boundary element values , , and . The difference between and does not exceed . When a crack is introduced at , 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 , is within of and is within of . This excludes positions near s = 0.5 where and 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 . 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 determined using the topological derivative and boundary element method are the same.
Figure 8(a) compares and , where it is seen that similar critical crack directions are determined using the topological derivative and boundary element methods. One can also see that is in good agreement with . is within of outside of the range . Within this range, the small values of and result in larger discrepancies.
The energy release rate is then calculated for when cracks of radius ɛ/T = 0.01 and direction are introduced. and are compared in Fig. 8(b), where it is again seen that the critical crack orientations determined using the topological derivative and boundary element methods are similar. Once more one can see that is in good agreement with . When , is within of . The difference between and grows as large as when . One conjecture is that when , 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 when cracks of various radii ɛ are introduced with direction and orientation . The energy release rate calculations and are compared in Fig. 9. When the crack radius is less than of the thickness, i.e., ɛ/T ≤ 0.05, the difference between and is small. As expected, for larger cracks, the difference between and increases. However, remains a reasonable approximation of even for cracks of considerable radius, e.g., ɛ/T = 0.25.
As indicated in Fig. 9, the energy release rate approximation associated with the boundary location is more accurate for longer crack radii than the same approximation associated with the boundary location . 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 [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. . These higher-order approximations take into account the boundary curvature, boundary location , 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 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 and crack radius ɛ. In contrast, the energy release rate approximation 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 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.
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 when cracks of various radii ɛ are introduced at boundary location with direction d = 1.00e3 and orientation . The stress intensity factors versus the position s are shown in Figs. 12(a)–12(c) for a crack of radius ɛ = (R2 − R1)/100. As in the previous example, one can see the stress intensity factors , , and computed using the topological derivative are in good agreement with the reference boundary element values , , and . The energy release rates appear in Fig. 12(d), where it is seen that the proposed approximation GTD is within of the reference boundary element value GBE. The critical crack positions determined using the topological derivative and boundary element methods are again the same. Figure 13 shows the energy release rates versus the crack radius ɛ/(R2 − R1) for the positions s = 0.3, s = 0.5, and s = 0.7. When the crack radius is less than of the annular radius, i.e., ɛ/(R2 − R1) ≤ 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.
An energy release rate approximation for half-penny cracks at any boundary location , 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 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.
In the equations that follow, the ɛ, , d, and α dependencies of the domain Ωɛ and boundaries ∂Ωɛ, γɛ, and ∂γɛ are dropped for conciseness.
In the equations that follow, the j = 1 indices are dropped for conciseness.
In the equations that follow, the s dependencies of the directions en, et, and eb are dropped for conciseness.
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.
- d =
- e =
- m =
outward normal vector on γɛ
- n =
outward normal vector on ∂Ω
- u =
- x =
- y =
scaled position vector
- f =
- s =
position on γɛ or γ1
- D =
- E =
- G =
energy release rate
- J =
- T =
outer stress (Cauchy stress tensor in Ω)
scaled half-space domain with a crack
prescribed boundary traction on γ1
inner stress (Cauchy stress tensor in
- hij =
- DT =
- GC =
- Kλ =
mode λ stress intensity factor
normalized mode λ stress intensity factor
- Tɛ =
Cauchy stress tensor in Ωɛ
- tp =
prescribed boundary traction on ∂Ω
- tR =
prescribed remainder traction on ∂Ω
prescribed remainder traction on γɛ
infinitesimal strain tensor
Subscripts or Superscripts
- b =
binormal direction on ∂γɛ
- BE =
reference (“exact”) numerical values
- (j) =
- n =
outward normal direction on ∂γɛ
- r =
outward normal direction at s = 0.5
- t =
tangential direction on ∂γɛ
- TD =
- z =
tangential direction at s = 0.5
- ɛ =
response quantity evaluated in Ωɛ
- θ =
binormal direction at s = 0.5