An auxiliary crack approach for efficient approximative crack tip loading analyses

An efficient approach based on the path-independent interaction-integral (I-integral) is suggested for assessing the crack tip loading in elastic structures without having to geometrically model a physical crack. Exploiting just the elastic solution of the uncracked structure, the I-integral is adapted incorporating the closed formulation of crack tip stress and displacement fields of an auxiliary crack, which in this approach is interpreted as hypothetical physical crack. Different specimens and crack configurations are investigated, sparing the expensive numerical modeling of discontinuities, and stress intensity factors are assessed according to the auxiliary crack approach. Various results are verified based on classic crack tip loading analyses.


Introduction
The fracture assessment plays a key role in the development and operation of light-weight structures, particularly in safety-relevant systems. This incorporates transportation systems for passengers and goods on P. O. Judt · A. Ricoeur (B) Institute of Mechanics, University of Kassel, 34125 Kassel, Germany email: ricoeur@uni-kassel.de ground, water or in the air, as well as hazardous facilities in the chemical and power industry. In this context important issues are where, into which direction and at which load cracks will evolve and if stable or unstable crack growth is to be expected.
The numerical modeling of structures with cracks in terms of a finite element (FE) framework is timeconsuming, particularly if several hypothetical crack configurations are investigated, and in general goes along with special mesh requirements to account for the stress and strain singularities at the crack tip. Special crack tip elements have been introduced to account for the singular stress fields (Barsoum 1974(Barsoum , 1976. In the extended FE method (XFEM), adapted shape functions are included in the element formulation to represent the discontinuity of crack faces and the typical stress and displacement fields at crack tips (Belytschko and Black 1999;Moës et al. 1999). A different approach is taken in the phase field method (Francfort and Marigo 1998). Here, the crack is represented by a diffuse phase field, continuously changing from solid to crack phase.
In classical numerical fracture mechanics, the crack tip loading analysis is a challenging task and different approaches have been presented. Here, pathindependent integrals are beneficial and can be evaluated along arbitrary integration contours far from the crack tip. First, the J -integral was introduced by Cherepanov (1967) and Rice (1968), representing a crack tip loading quantity in elastic and elastic-plastic materials. Later, J = J 1 was extended by a second coordinate J 2 in plane mixed-mode problems (Budian-sky and Rice 1973). In linear elastic fracture mechanics (LEFM) J k is related to the energy release rate (ERR) (Griffith 1921;Irwin and Kies 1952) and the stress intensity factors (SIF) (Bergez 1974). Later, the interaction integral, in the following denoted as I -integral, was introduced for plane elastic problems (Stern et al. 1976;Chen and Shield 1977), employing Betti's reciprocity theorem (Betti 1872) and the superposition of a physical and an auxiliary crack tip loading.
In the case of curved cracks, crack face integrals must not be neglected for the accurate calculation of Jand I -integrals Ricoeur 2013, 2015a). Especially when applying the I -integral in connection with straight auxiliary cracks and large integration contours, this leads to the necessity of crack face integrals along both the physical curved crack and the auxiliary straight crack, tangentially approaching each other and finally coinciding at the crack tip (Judt and Ricoeur 2015a;Gosz and Moran 2002); see Fig. 1. The auxiliary stress field at the position of the physical crack has to be incorporated as tractions on the physical crack faces and vice versa. The approach has successfully been applied to the loading analysis at single and multiple crack systems (Judt and Ricoeur 2015b).
The submodeling technique (Kuna 2013;Gloger et al. 2012) has been developed to avoid the explicit modeling of cracks in complex boundary value problems (BVP), providing an approximate crack tip loading. In a global analysis, the BVP without a crack is solved and resulting displacements are applied to the boundary nodes of a simple submodel containing the crack. The crack's loading state is subsequently determined applying a classical loading analysis approach. Applying the technique of crack weight functions (Bueckner 1970;Rice 1972) in connection with the principle of linear superposition, a further increase of efficiency can be achieved, exploiting stresses just in the crack plane of the global analysis and inserting them into SIF solutions of an arbitrary test load. Although the computation of complex crack problems thus is mapped onto much simpler BVP, the modeling of cracks is basically unavoidable.
In the following, an efficient approach is suggested for numerically determining approximate crack tip loadings. This approach is beneficial compared to others, as neither a crack is explicitly modeled in the BVP of the investigated structure, nor a submodel is required. In contrast to our approach in (Judt and Ricoeur 2015a), here the I -integral is applied to the numerical solu-tion of the defect-free BVP, where the physical crack is introduced in terms of discontinuous stress and displacement near tip fields of an auxiliary crack. With the freedom in the choice of location, orientation and length of the crack, very low effort is required for the assessment of crack tip loadings. The so-called auxiliary crack approach (ACA), on the other hand, lacks a rigorous theoretical background, and thus still requires fundamental research. Various plane structures are investigated applying the ACA. The promising results of SIF are compared to those of classical fracture mechanics calculations.

Interaction integral technique and the auxiliary crack approach
Crack tip loading quantities such as the J -integral, ERR or SIF provide information about crack growth and deflection behavior. The path-independent J k -integral (Rice 1968) is defined as where represents an infinitesimal and a finite integration contour enclosing the crack tip. In the latter case, includes the integration along crack faces and volume or inertia forces have to be excluded. The integrand Q k j is Eshelby's energy-momentum tensor (Eshelby 1975) and n j the unit normal vector on . The energy-momentum tensor is composed of the elastic strain energy density with stress and strain tensors σ i j and ε i j , respectively, and of the displacement gradient u i,k . The identity tensor is represented by δ k j . The superpositions of field variables of two loading conditions (a) and (b) in a linear elastic structure are obtained as Substituting the superimposed field variables into the energy-momentum tensor of Eq. (1) and considering Eq.
(2) provides The J k -integral of two loading conditions results from integrating Eq. (4) according to Eq. (1): The first coordinate (k = 1) of the last term in Eq. (5) represents the I -integral (Stern et al. 1976;Chen and Shield 1977). According to Irwin (1958) the SIF are related to the energy release rate G of a straight crack extension and thus to the J -integral according to with the generalized Young's modulus E = E for plane stress and E = E/(1 − ν 2 ) for plain strain conditions. In LEFM the superposition of two crack tip loading conditions (a) and (b) leads to where the index i represents the opening modes I or II for plane problems. Substituting the superimposed SIF Fig. 1 Integration paths containing the external contour 0 , auxiliary straight crack faces aux C and actual curved crack faces phys C , corresponding crack tip coordinate system x k of auxiliary and actual crack, global coordinate system x k into Eq. (6) provides From Eqs. (5) and (8) a relation of I -integral and SIF is obtained: The SIF K i and associated stress and displacement fields are related to possible solutions of the crack problem, no matter if going back to the actual loading or just a fictitious one, or even a fictitious auxiliary crack, as long as the crack faces (a) and (b) meet at a joint crack tip, see Fig. 1. Typically, solution (a) is obtained from numerical calculations representing the physical solution of the BVP with a crack exposed to the actual external loading. The fields (b) result from a closed-form analytical solution, e.g. the near tip stress and displacement fields (Irwin 1958;Williams 1957), representing an auxiliary solution in the I -integral. The unknown mode-I SIF K (a) I of the physical problem is then determined from the numerically evaluated Iintegral according to Eq. (9) with the auxiliary field variables related to a unit mode-I crack tip loading, i.e. K (b) Similar to the mode-I SIF, K II is determined choosing the auxiliary SIF according to a unit mode-II loading with K (b) Typically, to apply the I -integral to curved cracks, small integration contours are employed in order to avoid the implementation of additional line integrals, however, resulting in a lower accuracy. Judt and Ricoeur (2015a) showed, that in the case of curved cracks and large integration contours, additional crack face integrals along both physical and auxiliary cracks are necessary to maintain the path-independence of the I -integral see Fig. 1. The stresses along the physical and auxiliary crack surfaces, respectively, being attributed to both cracks in each case, are interpreted as surface tractions (Judt and Ricoeur 2015a;Gosz and Moran 2002). A superposition of SIF according to Eq. (7), however, is only valid if both physical and auxiliary cracks are approaching each other tangentially at the joint crack tip, which is intrinsically satisfied for straight cracks. Only in this case, near tip stresses on the ligament σ 0 i j (r ) = σ i j (r, ϕ = 0) may be superimposed, being controlled by loadsσ It is clear that at BVP without cracks, the crack length is a (a) = 0 and therefore K (a) i = 0. As a result, the interaction integrals according to Eqs. (9) or (12) must vanish and on the right hand side of Eq. (12) the contribution of the integration along the nonexisting physical crack face phys C vanishes independently. Accordingly, for a pure auxiliary crack problem the integrals along 0 and aux C are related as: For the numerical assessment of the crack tip loading following the ACA, the integration contour 0 completely encloses the auxiliary crack faces, the latter representing an equivalent physical crack, both intersecting at the external boundary of the structure, see Fig. 2. Meeting these requirements, the integral I aux 1 is independent of the chosen path 0 . The procedure now is as follows. Numerically, just the integral along 0 in Eq. (14) is calculated including continuous stress, strain and displacement fields σ (a) mn , ε (a) mn and u (a) i , see Eq. (4), from the FE solution of the crack-free BVP. The discontinuity of an equivalent crack problem is incorporated just by the auxiliary fields σ (b) mn , ε (b) mn and u (b) i taken from closed-form near tip solutions (Irwin 1958;Williams 1957). A single unit mode-I or mode-II auxiliary crack tip loading is chosen and substituted into Eq. (10) or (11) with I 1 = I aux 1 , finally obtaining SIF K I or K II of the presumed crack. For a vanishing contour aux C , Eq. (14) illustrates that I aux 1 = 0, thus the limiting condition K I/II = 0 for a → 0 is satisfied. While Eqs. (9) and (14) are exact for themselves, their equating is an approximation. The results in Sect. 3 will show that the approach denoted as ACA definitely provides suitable results for a variety of crack problems.

Examples of numerical crack tip loading analyses
Different crack problems have been analyzed with the ACA. Reference values of SIF have been generated computing the J -integral according to Eqs. (1) and (6) for the physical crack.
3.1 Crack in CT-specimen with U-notch Figure 2 shows a CT-specimen (h/w = 5/8) with a U-notch ( /w = 0.525) subjected to a symmetric displacement loadingū 2 . The normalized mode-I SIF is furthermore depicted vs. the normalized crack length. It is noted that the latter is calculated starting from the root of the notch, since the more common measurement starting at the points of load incidence obscures the results of shorter cracks. The solutions of the ACA are labeled as auxiliary crack and the reference values are denoted as physical crack. Both, the conventional analysis of the physical crack and the ACA provide similar results. The error depends on the crack length, reaching from 3% at a/w = 0.075 to 7% at a/w = 0.25, see Fig. 2, disregarding larger errors for very short cracks, where K I asymptotically approaches zero.

Crack in plate with hole
A plate with a hole and applied tensile stressσ 22 , as shown in Fig. 3a, was analytically investigated by Kirsch (1898). The stress field exhibits a stress concentration factor of S F = σ/σ 22 = 3 at the points of vertical tangents and S F = −1 at those of horizontal tangents. Thus, forσ 22 > 0 cracks may initiate along the horizontal and forσ 22 < 0 along the vertical axis. In Fig. 3b-d, normalized SIF of cracks emanating from the hole are depicted vs. the normalized crack lengths. The ratio of hole diameter to plate width in all specimens is d/w = 0.2. Figure 3b is obtained from a crack along the horizontal symmetry axis (ϕ = 0) emanating from the right side of the hole with tensile external loadingσ 22 . Just as with the CT-specimen in Fig. 2, the error of the ACA exhibits a pronounced dependence on the crack length with appropriate results for all except for very short cracks. The problem in the latter case is due to the large gradient of K I (a) for a → 0, which is not apparent in Fig. 3b and d due to the normalized graphs.
The normalized SIF of a crack emanating at ϕ = π/4 from the hole with again tensile loadingσ 22 are depicted vs. the normalized crack length in Fig. 3c. K II is obtained following Eq. (11), applying a unit mode-II auxiliary loading. Alternatively, K II is provided identically from the second coordinate I I 2 applying a unit mode-I auxiliary loading (Judt and Ricoeur 2015a). As a result of both ACA and physical crack approaches, K II (a) exhibits a horizontal tangent at a = 0, while K I (a) starts with a vertical one, albeit not being obvious in the normalized depiction of Fig. 3. As a consequence thereof, the analysis of K II based on the ACA is sufficiently accurate for all investigated crack lengths, whereas for K I errors decrease with increasing crack length, taking more than 10% just for very short cracks 2a/d < 0.2 and being almost exact for 2a/d > 0.7.
Finally, in Fig. 3d the mode-I SIF due to a compressive loadingσ 22 of a crack emanating vertically from the hole at ϕ = π/2 is plotted vs. the normalized crack length. Again, the error of K I from the ACA decreases quickly with increasing crack length, asymptotically approaching zero.
3.3 Crack due to Hertzian pressure Indentation fracture was first experimentally and analytically investigated by Hertz (1882a, b). The observed cone cracks in brittle materials such as glass or ceramics have further been analyzed by several researchers (Frank and Lawn 1967;Mouginot and Maugis 1985;Fischer-Cripps 1997;Strobl et al. 2017). In the state of crack initiation, however, the circumferential crack is cylindrical.
The cut-away view of a cylindrical specimen exposed to a cylindrical indenter punch loaded withσ 22 is depicted in Fig. 4. The indenter radius is r = 0.025R, the crack formation is supposed to occur at d = 0.06R and the cylinder height is h = 2R. The problem is simplified investigating just a planar model instead of a cylindrical one. The induced stress distribution is such that tensile stresses are present close to the upper specimen's surface, allowing for crack initiation at a distance d > r . With increasing distance to the upper surface of the specimen, the stresses field changes from tensile to compressive and thus crack growth is stopped.
This specific characteristic is also observed in the mode-I SIF of a crack growing vertically into the specimen. After a short range of positive crack driving force (K I > 0), the SIF becomes negative at a/ h ≈ 0.001, see Fig. 4. Contact of crack faces is not considered in the auxiliary fields and therefore unphysical overlapping occurs providing negative values. According to Fig. 4, the results from the ACA and the conventional analysis of the physical crack qualitatively coincide very well. Relative errors are only unacceptable for crack lengths going along with vanishing K I .

Conclusions
The auxiliary crack approach (ACA) is based on the interaction integral of fracture mechanics, commonly applied to calculate stress intensity factors (SIF) from the numerical solution of a boundary value problem with a crack in connection with auxiliary stress and displacement fields emanating from arbitrary auxiliary crack solutions. The idea is to exclude the crack physically from the structure and mathematically from the interaction integral, leaving the auxiliary crack and its discontinuous fields, adopted from available closedform near-tip solutions. The auxiliary crack then taking the role of the physical one, the SIF of a cracked structure can be determined reasonably, which has been demonstrated with examples. Qualitatively, the behavior of SIF vs. crack length is always predicted correctly by the ACA. Quantitatively, errors depend on the crack length, mostly being below 10%. Being an approximation on the one hand, the ACA gives the opportunity of very efficient fracture mechanical assessment of a structure on the other. Based on the data of e.g. a FE solution of the uncracked structure, arbitrary configurations of straight cracks are easily investigated without any further numerical effort.
From these promising results, the ACA might be established in fracture mechanical loading analyses in the future. Presently, the research in this field, however, is a work in progress, leaving open questions both from the theoretical and application-oriented points of view. A reliable application of the ACA in engineering fracture mechanics requires deeper insight into the nature of errors and their influencing factors.