The J-integral for mixed-mode loaded cracks with cohesive zones

The J-integral quantifies the loading of a crack tip, just as the crack tip opening displacement (CTOD) emanating from the cohesive zone model. Both quantities, being based on fundamentally different interpretations of cracks in fracture mechanics of brittle or ductile materials, have been proven to be equivalent in the late 60s of the previous century, however, just for the simple mode-I loading case. The relation of J and CTOD turned out to be uniquely determined by the constitutive law of the cohesive zone in front of the physical crack tip. In this paper, a J-integral vector is derived for a mixed-mode loaded crack based on the cohesive zone approach, accounting for the most general case of a mode-coupled cohesive law. While the J1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_1$$\end{document}-coordinate, as energy release rate of a straight crack extension, is uniquely related to the cohesive potential at the physical crack tip and thus to the CTOD, the J2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_2$$\end{document}-coordinate depends on the solution of the specific boundary value problem in terms of stresses and displacement gradients at the cohesive zone faces. The generalized relation is verified for the Griffith crack, employing solutions of the Dugdale crack based on improved holomorphic functions.


Introduction
In linear elastic fracture mechanics (LEFM) two different crack models are essentially used to describe fracture processes. In the first model all inelastic processes are assumed to take place in a sufficiently small zone in front of a sharp crack tip, the fracture process zone, and the crack tip exhibits singularities in all stresses. For this model, which will be denoted as "classical" model from now on, the stress intensity factors were introduced by Irwin (1957). Other crack tip loading quantities are the energy release rate (Griffith 1921) and the J-integral, which was independently introduced by Cherepanov (1967) and Rice (1968b). The latter loading quantity is motivated by the work of Eshelby (1956) which characterizes generalized forces on dislocations and point defects in elastic fields. The scalar J-integral was first limited to straight cracks, however, was extended to the J k -integral vector by Budiansky and Rice (1973), resulting in a path independent formulation for arbitrary crack configurations. This extension also gave rise to the J-integral vector criterion of crack deflection (Strifors 1974). The second model is the cohesive zone model (CZM) originally introduced by Barenblatt (1959) and more prominently known from (Barenblatt 1962), in which the fracture process zone in front of the crack tip is assumed be an extended cohesive zone, in which restraining stresses arise due to atomic separation, leading to the absence of singular crack driving stresses. The loading quantity in the CZM is the crack tip open-ing displacement (CTOD). As the classical model and the CZM both intend to describe the same problem, a relation between the loading quantities should exist, in case of an equivalency of the approaches. Respective proof for a straight mode I loaded crack was achieved by Rice (1968b), relating the scalar J-integral to the CTOD, by enclosing the cohesive zone by an integration contour. Just a few years before, Burdekin and Stone (1966) showed that for a Dugdale crack, the energy release rate and the CTOD are related. Further proof of this relation was again given by Rice (1968a), comparing the J-integral of a classical Griffith crack with crack tip singularity to a Griffith crack with a cohesive zone in front of the crack tip, adopting a Dugdale crack model (Dugdale 1960) with constant restraining stress. He showed that the length of the cohesive zone approaches zero for a vanishing ratio of applied stress to yield stress. The small scale yielding approximation being satisfied in this limiting case, the J-integral calculated for the CZM approaches the reference value of the equivalent classical Griffith crack, thus verifying his relation derived shortly before.
For the case of mixed-mode loading, where the scalar J-integral is replaced by a vector, a relation to the CZM or CTOD, respectively, has not been achieved before 2019 by Scheel and Ricoeur (2019), who investigated a matrix crack in the vicinity of a cohesive interface crack, introducing a J-integral vector for the cohesive crack tip. For a plane problem it holds two non-zero coordinates in general and the derived equation is considered to be the mixed-mode generalization of Rice's mode I relation (Rice 1968b). Mode coupling, in terms of shear stresses depending on normal separation and vice versa, is disregarded in (Scheel and Ricoeur 2019). While the J 1 -coordinate, as energy release rate of a straight crack extension, is straightforwardly derived, the second coordinate J 2 depends on the solution of the specific boundary value problem. For the J 1 -coordinate Nicholson (1993) also introduced a mixed-mode formulation disregarding mode coupling, however, fundamentally limiting his model to constant restraining stresses.
The J-integral for a classical bi-material interface crack was first introduced by Smelser and Gurtin (1977), focusing on the J 1 -coordinate of the vector. They stated that the J-integral for bi-materials with a straight bonding line is equal to the formulation for a homogeneous body. Khandelwal and Kishen (2006) extended the formulation to a J-integral vector, con-cluding that the J 2 -coordinate is path dependent and even non-existent if the integration path approaches the crack tip. Consequently, a geometric parameter has to be introduced requiring physical interpretation.
In this work, the J-integral vector of cohesive zones is derived, accounting for mode coupling, bi-material crack faces and arbitrary cohesive laws, represented by cohesive potentials. The fundamental goal is to give the opportunity of an alternative calculation of the Jintegral vector based on data provided by CZM and to investigate the equivalence with respect to the classical crack tip approach with path independent contour integrals. Applying the Leibniz integral rule, J 1 is complemented by a coupling integral term, still being uniquely related to the CTOD of normal and tangential separations. It is further shown that J 1 equals the cohesive potential at the physical crack tip. J 2 , on the other hand, turns out not to be uniquely related to the CTOD but is calculated from an integration along the cohesive zone. This might give the opportunity of a path independent calculation of J 2 for bi-material interface cracks.
A verification of the provided relations requires stresses and displacement gradients at the cohesive zone faces. These are taken from a Griffith crack with cohesive zones in front of the physical crack tips, thus constituting a Dugdale crack. Kolosov's equations (Kolosov 1909) and holomorphic functions (Muskhelishvili 1963) or alternatively Westergaard stress functions (Westergaard 1939) are the basis of the closedform solution. For the Dugdale crack, a variety of functions are provided by literature (e.g Burdekin and Stone (1966); Hayes and Williams (1972); Becker and Gross (1988); Tada et al. (2000)). These, however, require case-by-case analyses to cover the entire domain. The holomorphic functions introduced in Sect. 3.2 cope with this drawback, yielding a monolithic solution.
2 The J-integral vector of a mixed-mode crack with a cohesive zone Rice (1968b) derived a relation between the J-integral and the CTOD by shrinking an integration contour to the cohesive zone in front of a single mode I loaded crack tip, yielding the well-known relation where δ t is the opening separation at the physical crack tip and σ (δ) is the restraining stress. In order to derive a generalized mixed-mode formulation, an analogous procedure is employed, starting with the J-integral vector, neglecting body forces, which is given by Budiansky and Rice (1973): Here, S represents the arbitrarily chosen integration contour, w is the potential energy density, n k is the normal vector, t i is the traction vector and u i,k is the displacement gradient. Tensors and differential operations are depicted in index notation, implying summation over repeated indices, holding values one and two for plane mode I/II problems and one to three for the spatial mode I/II/III case outlined in the Appendix B.
A comma further denotes a partial spatial derivative. It has to be noted that the J k equals a configurational force F k (Gurtin and Podio-Guidugli 1996) with the relation J k = −F k , thus holding the unit [N]. A surface integration was replaced by a line integration in Eq.
(2), thus J k has to be interpreted as crack driving force per unit thickness. The integration contour of Eq.
(2) is now shrunk to the cohesive crack faces. In order to keep it general, an interface crack in an arbitrarily loaded bi-material body is considered, see Fig. 1, yielding where the superscripts +/-refer to the two subdomains connected to the positive or negative crack face. With an idealized cohesive zone of vanishing lateral extension, the x 2 -coordinate of S ± is zero, i.e. x − 2 = x + 2 = 0, and opposing orientations of the integral paths along the interface yield dS + = −dS − = −dx 1 . Eq. (3) is thus rewritten as where Δa is the length of the cohesive zone and position of the fictitious crack tip, respectively. Since the traction vectors on the positive and the negative cohesive crack faces are continuous, they are replaced by a Fig. 1 Interface crack in an arbitrarily loaded bi-material body with cohesive zone of length Δa, J-integral vector J and integration contours S, S + , S − , coordinate system (x 1 , x 2 ) at the physical crack tip and crack tip opening displacements δ t i restraining traction vector, i.e.
yielding the following formulation of the J-integral vector of the crack with cohesive zone: In order to complement the interface conditions, tractions and displacements on the ligament are continuous and in the cohesive zone, i.e. between physical and fictitious crack tips, displacements are discontinuous.

The J-integral vector coordinate J 1
The first coordinate of the normal vector being zero in the coordinate system of Fig. 1, i.e. n ± 1 = 0, the potential energy density in Eq. (6) vanishes for k = 1, leaving the scalar product of the traction vector and the displacement gradient for the coordinate J 1 : where δ i (x 1 ) is the local separation. For the special case of restraining stresses only depending on the mode associated separation, i.e.
is rewritten, applying the substitution rulē The integration limits of Eq. (8) then turn into the separations at the beginning (x 1 = 0) and the end (x 1 = Δa) of the cohesive zone, i.e.
finally leading to the mixed-mode generalization of J 1 being This relation between J 1 and the CTOD for a decoupled mixed-mode loading case has already been given by Scheel and Ricoeur (2019) and is a rather straightforward generalization of Eq. (1), just adding a shear term. The integrals describe the areas underneath the traction-separation curves in the interval [0, δ t i ] and the dissipated surface energy density in case of crack growth, i.e. δ t i = δ iC , with a critical separation δ iC . In the more general case of restraining stresses depending on all separations, i.e. t R i = f (δ 1 , δ 2 ), a simple substitution according to Eq. (9) is not applicable. J 1 is rather formulated employing a generalized formulation of the Leibniz integral rule, see Appendix B, from which follows Inserting Eq. (12) into Eq.(7) finally yields The decoupled case of Eq. (11) is obtained from Eq. (13), as for t R i (δ i ) the derivatives ∂t R 1 /∂δ 2 and ∂t R 2 /∂δ 1 vanish.
The fundamental extension in this generalized formulation are the two double integral terms. At first glance, due to the gradients of the separations δ i,1 , Eq. (13) requires the solution of the specific boundary value problem, thus J 1 apparently does not uniquely depend on the CTOD as in the decoupled case. This issue is illuminated introducing the cohesive potential Π coh , from which the restraining tractions are extracted as: Inserting Eq. (14) into Eq. (7) yields Introducing the total derivative of the cohesive potential, Eq. (15) is rewritten according to With the separations at the beginning and the end of the cohesive zone from Eq. (10), the J 1 -coordinate finally reads J 1 obviously does not depend on the specific boundary value problem, and Eq. (18) is valid for any kind of traction separation law and fracture process, as long as it is uniquely described by a cohesive potential Π coh (δ 1 , δ 2 ). Eq. (18) also provides the crack tip loading quantity for a fracture criterion of mixed-mode loaded cohesive cracks, in terms of where G c is the critical energy release rate and crack growth resistance, respectively. Reducing the two loading quantities δ t 1 , δ t 2 to a single energy density Π coh (δ t 1 , δ t 2 ) introduces an equivalent mixed-mode crack tip loading quantity. Hence, the cohesive potential function needs to be formulated carefully, based on empiric data of crack growth initiation.
Eq. (18) is equivalent to Eq. (11) for the decoupled case t R i (δ i ), where the cohesive potential is decomposed into the mode-related contributions, according to For the general case of t R i = f (δ 1 , δ 2 ), Eqs. (18) and (13) are equivalent, revealing that the double integrals in Eq. (13), apparently depending on the boundary value problem via δ i,1 , are not problem-specific. A general proof for this equivalency is given in the Appendix A.

The J-integral vector coordinate J 2
The second coordinate of the normal vector on the cohesive zone faces in Eq. (6) is non-zero in the coordinate system of Fig. 1, i.e. n + 2 = −n − 2 = 1; therefore, the potential energy density in J 2 does not vanish either. Furthermore, the displacement separation cannot be introduced in the same way as δ i is for J 1 in Eq. (7), since the derivative δ i,2 requires the coordinate x 2 at x ± 2 = 0. Therefore, a separationδ i (x 1 , x ± 2 ) between two arbitrary points below and above the cohesive crack faces needs to be introduced, subsequently taking the limit x 2 → 0, i.e.
The partial derivatives in Eq. (21), being executed with respect to x 2 , and the integration, on the other hand, being performed with respect to x 1 , doesn't allow for an application of the Leibniz integral rule. The same holds for an alternative formulation based on the cohesive potential, which was appropriate just with J 1 , whereupon akin to Eq. (46) in Appendix A is inserted into Eq. (21), yielding where Eq. (14) again introduces the restraining stresses t R i . Thus, it is not constructive to introduce a displacement separation in J 2 . According to Eq. (23) the CTOD of a cohesive zone approach are not uniquely related to the coordinate J 2 of the J-integral. Starting from Eq. (6) with n + 2 = −n − 2 = 1, the potential energy density with Ω ± 12 = (u ± 1,2 − u ± 2,1 )/2, ε ± i j and σ ± i j being the rigid body rotations, strains and stresses on the cohesive crack faces, respectively. J 2 , in contrast to J 1 , depends on the solution of the boundary value problem, requiring integration along the cohesive zone. The J-integral vector then finally is given by its coordinates according to Eqs. (18) and (24), reading for the general mode-coupled case, whereas for the decoupled case, t R i = f (δ i ), Eqs. (11) and (24) yield with e k being the unit vectors corresponding to the coordinate system of Fig. 1. A generalization including the anti-plane shear mode III is given in the Appendix B.
Albeit not being the primary focus of the paper, it shall finally be noted, that the calculation of J 2 of a bimaterial interface crack by integrating along Δa does not allow for any path dependence as it is known from the classical contour integral formulation (Khandelwal and Kishen 2006). Due to the absence of a sharp crack tip, the typical oscillating stress and strain singularities appearing in the integrand and giving rise to discussions on the physical interpretation of J 2 in the classical approach, are not an issue in Eq. (24).

Verification for a Griffith crack based on a Dugdale crack solution
In order to verify the generalized J-integral vector of a mixed-mode loaded cohesive zone, a Griffith crack is the appropriate example for a closed-form solution. In Fig. 2 a Griffith crack is depicted with the J-integral vector indicated at one of the physical crack tips. Cohesive zones with lengths Δa I /I I are introduced constituting a Dugdale crack. The restraining stresses are thus acting as surface tractions on the Dugdale crack in the intervals −a − Δa ≤ x ≤ −a and a ≤ x ≤ a + Δa. Unlike depicted in Fig. 2, the mode-associated cohesive zone lengths Δa I and Δa I I are different in general, just as the restraining stresses are functions of the position. While assuming a cohesive law akin to perfect plasticity, i.e. t R 1 = τ 0 and t R 2 = σ 0 , yields a closed-form solution, it excludes the coupled case of t R i = f (δ 1 , δ 2 ). Simplifying assumptions are furthermore a homogeneous mono-material body, equal mode I and II loading σ ∞ = τ ∞ = Σ ∞ and equal cohesive tractions σ 0 = τ 0 = Σ 0 . The latter assumption is dispensable; however, dissimilar σ 0 and τ 0 do not provide deeper insight within the context of verification.
The procedure now is to calculate the J-integral vector with the generalized relation according to Eq. (26) on the one hand and by stress intensity factors of the classical Griffith crack on the other, the latter acting as a reference. A similar procedure based on the same assumptions was applied by Rice (1968a) to verify the mode I relation according to Eq. (1). While the verification of the mixed-mode J 1 is straightforward and surprising results may not be expected, J 2 requires the displacement gradient and stress solutions of the Dugdale crack and a verification of the derived relation is pending at this point.

Verification of J 1
In order to calculate J 1 in Eqs. (11), (13) and (18), respectively, the crack tip opening displacements are required, depending on loading and crack length. Their calculation necessitates the determination of the cohesive zone lengths first. They are obtained by the requirement that there is no singularity of the crack driving stresses at the fictitious crack tips. Applying the superposition principle, separating loads at infinity Σ ∞ and the cohesive zone Σ 0 , see Fig. 3, the cohesive zone lengths Δa I /I I are calculated from the above requirement, introducing the stress intensity factors of the two subproblems by means of crack weight functions (Kuna 2013) or complex potentials (Tada et al. 2000), leading to the following equation: Calculating the integrals yields the cohesive zone lengths With the simplifying assumptions made (τ 0 = σ 0 = Σ 0 and τ ∞ = σ ∞ = Σ ∞ ), it follows that the cohesive zone lengths are equal, i.e. Δa I = Δa I I = Δa. The CTOD as the displacements at the physical crack tip are obtained as with E = E for plane stress and E = E/(1 − ν 2 ) for plane strain, where E is Young's modulus and ν is Poisson's ratio. In the decoupled case the J 1 -coordinate is calculated according to Eq. (11), where inserting Eq. (29) yields: On the other hand, the J 1 -coordinate of the Griffith crack is classically calculated as The ratio of the two J 1 -coordinates of Eqs. (30) and (31) and the normalized cohesive zone length Δa/a are plotted in Fig. 4 versus the ratio of applied to restraining stress. With the cohesive zone vanishing for a stress ratio tending to zero, the small-scale yielding approximation is satisfied. J L E 1 being exact in that limiting case, the two functions of Eqs. (30) and (31) asymptotically approach each other, see Fig. 4. The coordinate J 1 of the generalized relation of Eqs. (26) and (11), respectively, can be considered verified.  The calculation of these fields is based on Kolosov's equations (Kolosov 1909), where z = x + iy is the complex coordinate of the coordinate system in Figs. 2 and 3, respectively, μ is the shear modulus and Φ, Ψ are holomorphic functions with their derivatives Φ , Φ , Ψ . Bars on quantities denote the complex conjugates, e.g.z = x − iy. From Eq. (32) the stresses are obtained as: For holomorphic functions, the Cauchy-Riemann equations are satisfied, e.g. for Φ reading with the Wirtinger derivatives (Wirtinger 1927) being defined as From Eq. (34) it follows that which, inserted into the Wirtinger derivatives of Eq.
(32) and (37) the displacement gradient is calculated yielding the coordinates The gradients u 1,1 = ε 11 and u 2,2 = ε 22 , of course, could alternatively be calculated with Hooke's law. The holomorphic functions provided in literature often only cover parts of the whole solution or require case-by-case analysis. Improved functions have thus been set up here, reading providing the required fields in the entire plate, with c = a + Δa. The stresses and displacement gradients are finally obtained inserting Eq. (39) into Eqs. (33) and (38). Figures 5 and 6 illustrate the relevant quantities in the plane of the crack, i.e. at y = 0. In Fig. 6 only the positive crack face is considered and f + (−x) = f − (x) holds for the negative one in the mixed-mode loading case. For this example the physical crack length was chosen a = 2 mm, the cohesive zone length Δa = 1 mm and κ = 1.6 under plane strain conditions. It is self-evident that σ ± 12 and σ ± 22 are zero for pure mode I and II, respectively, loading. The crack driving stresses σ ± i2 are nonsingular at both the physical and the fictitious crack tips in all loading cases; however, the stresses σ ± 11 are singular at the physical crack tip at mode II and mixed-mode loading, see Fig. 5 b, c, being attributed to the discontinuity of crack face stresses at |x| = a. This aspect was also found and shortly discussed by Becker and Gross (1988), pointing out a logarithmic singularity instead of the well-known r −1/2 -behavior of stresses at sharp crack tips. While the stresses σ ± 11 are axisymmetric for mode I, see Fig.  5 a, and point symmetric for mode II, see Fig. 5 b, the mixed-mode loading exhibits an asymmetry of these stresses, see Fig. 5 c. Singularities in the displacement gradients at the physical crack tip are the consequence of the singularity in the σ 11 stresses, easily reproduced applying Hooke's law. To verify the graphs of Figs. 5 and 6 the method of distributed dislocations (Bilby et al. 1963;Bilby and Eshelby 1968;Weertman 1996), as a common approach in LEFM, has been employed in this work, confirming the results derived from Eqs. (33), (38) and (39).
The superposition principle gives the opportunity to further simplify the J 2 -coordinate, separating the stresses, strains and rotations in Eq. (24) into moderelated parts, leading to: Taking into account the following conditions of symmetry and antisymmetry of the mode-related contributions, which can be derived for the stresses from ; Ω I I + 12 = Ω I I − 12 , u I + 1,2 = −u I + 2,1 → Ω I + 12 = u I + 1,2 , Analogous to the verification of J 1 , the ratio of the J 2 -coordinates is examined and plotted in Fig. 7 versus the ratio of applied stress to restraining stress. In contrast to J 1 in Fig. 4, the J 2 -coordinates are on the one hand not asymptotically approaching each other and on the other hand J 2 /J L E 2 ≤ 1. After all, the coordinates become equal for the depicted vanishing cohesive zone length, thus the formulation of the J 2 -coordinate of Eqs. (24) and (42), respectively, is considered verified. Finally, it shall be noted that both J 1 /J L E 1 and J 2 /J L E 2 in Figs. 4 and 7 tend to infinity for Σ ∞ /Σ 0 → 1, where an interpretation of failure within the context of fracture mechanics gives way to a global collapse.

Conclusion
The vector of the J-integral has been derived from a cohesive zone approach of a crack subject to mixedmode loading, accounting for arbitrary mode-coupled cohesive laws. Just as for the pure mode I case, the coordinate J 1 is uniquely connected to the normal and now also tangential CTOD emanating from the CZM. The generalized relation is either given by means of the cohesive potential or includes line integrals along the cohesive zone. The conjunction of J 1 and the CTOD constitutes the equivalence of the CZM and the classi-cal singular crack tip approach of LEFM with respect to the energy release rate of a straight crack extension. Restricting to J 1 , the equivalence of the approaches holds for arbitrary traction separation laws, including elasto-plasticity, not being limited by the cohesive zone size. The coordinate J 2 , on the other hand, is not uniquely related to the CTOD, but is calculated by integration of stress and displacement gradients along the cohesive zone. Consequently, two loading scenarios being equivalent with respect to the CTOD, might yield different J-integral vectors and thus stress intensity factors, if applicable. In the mixed-mode case, CZM and classical singular crack tip approaches of LEFM thus would have to be considered as non-equivalent.
In particular concerning J 2 of bi-material interface cracks, a calculation based on CZM might yield path independent results providing a unique physical interpretation, which will have to be investigated in the future.
Acknowledgements The authors would like to thank the German National Science Foundation (DFG) for financial support, Dr. R. Boukellif for providing the calculations with the distributed dislocations technique and Dr. D. Wallenta for mathematical support with the Leibniz integral rule.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.
Inserting Eqs. (44) and (47) into Eq. (13) finally gives Appendix B: J k -integral vector for mixed-mode I, II, III loading The J-integral vector along the front of a 3D crack requires integration along the surfaces of the cohesive zone, see Fig. 8. With dS + = −dS − = −ΔBdx 1 and for indices holding values one to three, Eq. (7) is complemented according to whereupon J 1 /ΔB is the average within a crack front segment ΔB and δ i as well as t R i are assumed constant in ΔB, representing averages on their part.