domain: spi phys chim

text:
A Vattré
email: aurelien.vattre@onera.fr

Kinked and forked crack arrays in anisotropic elastic bimaterials

Keywords: Continuous distributions of dislocations, anisotropic elasticity theory, heterogeneity, crack arrays, kinked and branched cracks, mixed-mode configurations

The fracture problem of multiple branched crack arrays in anisotropic bimaterials is formulated by use of the Stroh formalism to the linear elasticity theory of dislocations. The general full-field solutions are obtained from the standard technique of distributed dislocations along finite-sized cracks of arbitrary shapes, which are embedded in dissimilar anisotropic half-spaces under far-field stress loading conditions. The bimaterial boundary-value problem leads to a set of coupled integral equations of Cauchy-type that is numerically solved by using the Gauss-Chebyshev quadrature scheme with appropriate boundary conditions for kinked and forked crack arrays. The path-independent J k -integrals as crack propagation criterion are therefore evaluated for equally-spaced cracks, while the limiting configuration of individual cracks is theoretically described by means of explicit expressions of the local stress intensity factors K for validation and comparison purposes on several crack geometries. The short-range interactions resulting from the idealized configurations of infinitely periodic cracks are investigated as well as various size-and heterogeneity-effects on the mixed-mode cracks in complex stress-state environments. The influences of anisotropic elasticity, elastic mismatch, applied stress direction, inter-crack spacings and crack length ratios on the predictions from the J k -and Kbased fracture criteria are examined in the light of different configurations from the single kinked crack case in a homogeneous medium to the network of closely-spaced forked cracks in presence of bimaterial interfaces.

Introduction

The problem of cracks and fracture mechanics is an important wide-ranging research field of engineering science and mathematical physics since the pioneering theory developed by [START_REF] Griffith | The phenomenon of rupture and flow in solid[END_REF], which forms the foundation of the modern continuum fracture mechanics. The cracks of various types and scales are found in many manufactured components, while their effects on the mechanical strength of solids are of continued concern. In many problems, failures and progressive damages are often related to interactions of multiple cracks in brittle or quasi-brittle materials such as concrete, rocks, ceramics and interface-dominated composite materials in bonded structures. In particular, the interaction effect among cracks has received considerable attention in the past decades, mainly because the disastrous failure events are preceded by the interaction and coalescence of those cracks. The theoretical analysis of the energetically favorable conditions for branched cracks to merge is therefore of importance in assessing structural integrity in many situations, including stress corrosion cracking, corrosion fatigue and thermal fatigue. The present work is part of the long-standing problems in fracture mechanics, namely 1) the determination of the stress fields produced by an infinite array of finite-length kinked and forked cracks of arbitrary geometry in anisotropic bimaterials, as well as 2) the corresponding configurational driving forces acting at the crack singularities. The bimaterials are subjected to arbitrary far-field stress state, and the arrays of branched cracks are represented as suitable continuous distributions of regularly-spaced dislocations, which are determined by a Cauchy-type singular integral equation and solved using series of Chebyshev polynomials. The isolated straight crack in an isotropic and homogeneous elastic medium is therefore a special limiting case, and is considered for validation purposes.

In many fracture mechanics problems in infinite isotropic solids, a large class of analytical solutions under plane stress and strain conditions in two dimensions is based on the fundamental work of complex potentials of [START_REF] Muskhelishvili | Some Basic Problems of Mathematical Theory of Elasticity[END_REF]. This mathematical theory of plane elasticity has been used to treat the finite-length branched cracks [START_REF] Chatterjee | The stress field in the neighborhood of a branched crack in an infinite elastic sheet[END_REF][START_REF] Kitagawa | Crack-Morphological Aspects in Fracture Mechanics[END_REF][START_REF] Lo | Analysis of branched cracks[END_REF] as well as the doubly-periodic array of straight cracks in an infinite isotropic medium [START_REF] Fil'shtinskii | Interaction of a doubly periodic system of rectilinear cracks in an isotropic medium[END_REF][START_REF] Ioakimidis | Doubly-periodic array of cracks in an infinite isotropic medium[END_REF][START_REF] Panasyuk | A general method of solution of two-dimensional problems in the theory of cracks[END_REF], within which the complex potential functions are taken as the unknown functions, and the Cauchy-type singular integral equations are established by obtaining these functions subjected to appropriate boundary conditions. For most of these solutions, the analysis is such that the limit for infinitesimally small kinks cannot be obtained readily from the analysis for finite kinks, leading to asymptotic expansions for stress intensity factors at the tips of infinitesimally small kinks [START_REF] Amestoy | Crack paths in plane situations-II. Detailed form of the expansion of the stress intensity factors[END_REF][START_REF] Bilby | The crack with a kinked tip[END_REF][START_REF] Cotterell | Slightly curved or kinked cracks[END_REF][START_REF] Lo | Analysis of branched cracks[END_REF][START_REF] Wu | Maximum-energy-release-rate criterion applied to a tension-compression specimen with crack[END_REF]. Due to mathematical difficulties in many complicated branched-crack situations under non-uniform loading conditions, however, the Muskhelishvili complex potentials can not be readily applicable, while the alternative dislocation-based formalism [START_REF] Bilby | Dislocations and the Theory of Fracture[END_REF][START_REF] Delameter | Weakening of an elastic solid by a rectangular array of cracks[END_REF][START_REF] Stroh | Dislocations and Cracks in Anisotropic Elasticity[END_REF] can be used for more practical boundary value problems of structural materials containing branched cracks.

Dislocations and the corresponding elasticity of extrinsic and intrinsic dislocations are introduced into the theory of fracture mechanics in two directions. First, the nucleation and plastic flow of dislocations from the stressed crack tips are of great importance to determine the growth resistance of ductile metals (Kelly et al., 1967, Rice and[START_REF] Rice | Ductile versus brittle behaviour of crystals[END_REF], exhibiting the physical competition between Griffith cleavage and dislocation-mediated plastic shear at crack tips. Secondly, planar dislocations are viewed as mathematical elements to describe the displacement discontinuities along faults within the Earth's crust [START_REF] Chinnery | The dislocation fault model with a variable discontinuity[END_REF]Petrak, 1967, Weertman, 1964), semicoherent interfaces [START_REF] Bilby | Continuous distributions of dislocations: New application of the methods of non-Riemannian geometry[END_REF][START_REF] Bullough | Continuous distributions of dislocations: Surface dislocations and the crystallography of martensitic transformations[END_REF][START_REF] Vattré | Determining the Burgers vectors and elastic strain energies of interface dislocation arrays using anisotropic elasticity theory[END_REF] as well as crack surfaces [START_REF] Bilby | Dislocations and the Theory of Fracture[END_REF][START_REF] Bilby | The crack with a kinked tip[END_REF][START_REF] Weertman | Continuum distribution of dislocations on faults with finite friction[END_REF] in classical elasticity theory. From the crystallographical picture, the crack plane dislocations of a mode I crack are climb edge dislocations, while the glide edge dislocations and screw dislocations are related to the mode II and III cracks, respectively. In practice, a complex crack problem leads also to a combinaison of continuum distribution of infinitesimal edge-and screw-type dislocations, such that the proper stress field produced by the cracks are calculated as a convolution of the multiple stress fields of dislocations with unknown distribution functions of dislocation densities. These densities are determined by satisfying the crack surface traction-free conditions, while the mixed-mode stress intensity factors are then written with respect to the local dislocation densities [START_REF] Hills | Solution of Crack Problems: The Distributed Dislocation Technique[END_REF].

Since the seminal work of [START_REF] Volterra | Sur l'équilibre des corps élastiques multiplement connexes[END_REF] that lays out the theory of dislocation fields generated by a discontinuity of displacement in linear elastic bodies, much attention has been paid to the determination of the elastic interaction between single arbitrary-oriented dislocations and traction-free surfaces or internal interfaces, which is one the most important part of the elasticity theory of dislocations [START_REF] Hirth | Theory of Dislocations[END_REF]. Extensions of the single dislocation problems to periodic arrays of equally-spaced extrinsic and intrinsic dislocations have been addressed in both isotropic [START_REF] Chou | Instability of edge dislocation walls in a two-phase isotropic medium[END_REF][START_REF] Grekov | Interaction of edge dislocation array with bimaterial interface incorporating interface elasticity[END_REF][START_REF] Gutkin | Equilibrium position of misfit dislocations[END_REF][START_REF] Lubarda | Energy analysis of dislocation arrays near bimaterial interfaces[END_REF] and anisotropic [START_REF] Chou | Interaction of parallel dislocations in a hexagonal crystal[END_REF][START_REF] Chu | Elastic fields due to dislocation arrays in anisotropic bimaterials[END_REF][START_REF] Gosling | The energy of arrays of dislocation in an anisotropic half-space[END_REF][START_REF] Vattré | Three-dimensional interaction and movements of various dislocations in anisotropic bicrystals with semicoherent interfaces[END_REF] bimaterial solids, and this is the latter that is explored in the present two-dimensional analysis. Although the numerous discrete dislocation dynamics codes in the current literature are capable of exploring the collective evolution of random-oriented forest dislocations, lattice-and nodal-based simulations are essentially resorted to using the isotropic elastic solutions. In comparison with results from the full anisotropic elasticity, the elastic field solutions for dislocations [START_REF] Asaro | A further synthesis of sextic and integral theories for dislocations and line forces in anisotropic media[END_REF], the corresponding elastic energy and driving forces (Aubry et al., 2011, Barnett and[START_REF] Barnett | The elastic energy of a straight dislocation in an infinite anisotropic elastic medium[END_REF], the reactions between glide dislocations and forest dislocations [START_REF] Püschl | Reactions between glide dislocations and forest dislocations in anisotropic bcc metals[END_REF] as well as the local relaxations between coplanar dislocations (Vattré, 2017a,b) can be altered by the isotropic elastic approximation. These significant consequences on the fundamental behavior of such elementary defects in solids suggest unequivocally that similar important effects due to elastic anisotropy are also prevalent and transferable in fracture mechanics.

In dislocation mechanics, the analysis of periodic dislocation arrays is of importance in modeling grain boundaries and heterophase interfaces as well as in understanding persistent slip bands that have dislocation cell structures. On the other hand, examination of periodic arrays of cracks is important in investigating the fracture processes of crack accumulation and coalescence [START_REF] Kachanov | Elastic solids with many cracks: A simple method of analysis[END_REF][START_REF] Kachanov | On the problems of crack interactions and crack coalescence[END_REF], specially with specific collinear and straight cracks, for which the solutions based on the Muskhelishvili theory have been obtained by [START_REF] Koiter | An infinite row of collinear cracks in an infinite elastic sheet[END_REF], and more recently, in predicting the propagation of segmented crack fronts in mode I+III fracture [START_REF] Lazarus | Mode I+III multiscale cohesive zone model with facet coarsening and overlap: Solutions and applications to facet orientation and toughening[END_REF][START_REF] Leblond | Development of fracture facets from a crack loaded in mode I+III: Solution and application of a model 2D problem[END_REF][START_REF] Vattré | Mechanical interactions between semicoherent heterophase interfaces and free surfaces in crystalline bilayers[END_REF]. The crack front segmentation in mixed mode I+III and I+II+III is therefore an illustrative example, for which the crack kinking and forking processes cannot be ignored when understanding the experimental crack front profiles and trajectories. Changes in crack path are generally induced by multi-axial stresses and/or the presence of microstructural inhomogeneities close to the crack tips, which can have a beneficial effect in the sense that crack kinking and forking improve the fracture toughness as well as the resistance to crack growth [START_REF] Suresh | Plastic near-tip fields for branched cracks[END_REF]. Furthermore, non-straight cracks with zig-zag growth of crack and the crack branching are also often be found in materials subjected to stress corrosion cracking and corrosion fatigue under mixed-mode loading conditions [START_REF] Kitagawa | Crack-Morphological Aspects in Fracture Mechanics[END_REF]. Such features have motivated and still motivate the development of various failure stress-and energy-based criteria for predicting the direction of crack propagation [START_REF] Amestoy | Crack paths in plane situations-II. Detailed form of the expansion of the stress intensity factors[END_REF][START_REF] Azhdari | Energy-release rate and crack kinking in anisotropic brittle solids[END_REF][START_REF] Cotterell | Slightly curved or kinked cracks[END_REF][START_REF] Erdogan | On the crack extension in plates under plane loading and transverse shear[END_REF][START_REF] Gao | Slightly curved or kinked cracks in anisotropic elastic solids[END_REF][START_REF] Goldstein | Brittle fracture of solids with arbitrary crack[END_REF][START_REF] Ichikawa | A critical analysis of the relationship between the energy release rate and the stress intensity factors for non-coplanar crack extension under combined mode loading[END_REF][START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF][START_REF] Leguillon | Crack deflection in a biaxial stress state[END_REF][START_REF] Lin | Criterion for initiation of cracks under mixed-mode I+III loading[END_REF], Sih, 1974), importantly designed to investigate the configurational bifurcation and stability of cracks loaded under mixed-mode conditions.

By use of the linear elasticity theory of dislocations, a general singular integral equation formalism for two-dimensional fracture analysis of equally-spaced configurations of multiple branched cracks in bimaterials is proposed. The bimaterials are composed of two dissimilar orthotropic half-spaces under arbitrary far-field stress loading conditions. The content of the paper is organized as follows. Section 2 describes the boundary-value problem, within which the continuously distributed dislocations along the crack segments lead to a set of coupled integral equations of Cauchy-type that is numerically solved by the Gauss-Chebyshev quadrature scheme with appropriate boundary conditions for kinked and forked cracks. The full-field solutions are therefore used to evaluate the path-independent J k -integrals as crack propagation criterion. The limiting case of individual cracks in homogeneous and isotropic materials is formulated in section 3, while the corresponding stress intensity factors are provided for validation purposes and further analyzed to investigate the interaction effects resulting from the idealized arrangements of infinitely periodic cracks. More advanced application examples are investigated to illustrate various size-and heterogeneityeffects on the elastic properties of mixed-mode cracks with arbitrarily branched geometries in anisotropic bimaterials.

Problem formulation

In the two-dimensional Cartesian coordinate system, the elastic field solutions due a network of periodically spaced dislocations with a inter-dislocation distance L are used to solve the boundary value problems of finite kinked and branched cracks in anisotropic bimaterials. For a single dislocation in bimaterials, which is not strictly located at the interfaces, the derivation 2 of the Green's functions has been presented by [START_REF] Barnett | An image force theorem for dislocations in bicrystals[END_REF], [START_REF] Belov | Elastic fields of dislocations piercing the interface of an anisotropic bicrystal[END_REF] and [START_REF] Ting | Image singularities for anisotropic elastic half-spaces and bimaterials[END_REF][START_REF] Ting | Anisotropic elasticity: theory and applications[END_REF], while the dislocation arrays have been formulated by Chou et al. (1975), [START_REF] Hirth | Stress fields of dislocation arrays at interfaces in bicrystals[END_REF] and more recently by [START_REF] Chu | Elastic fields due to dislocation arrays in anisotropic bimaterials[END_REF].

Using the Green's functions for a dislocation array with arbitrary characters, the crack arrays are therefore considered as the continuous distributions of dislocations, for which the density of defects must fulfill the traction-free boundary conditions at the crack faces. The singular integral equations of Cauchy-type are obtained, and finally solved by the Gauss-Chebyshev quadrature scheme developed by [START_REF] Erdogan | Numerical solution of singular integral equations[END_REF]. From the field solutions, the contour J k -integrals as crack propagation criterion are consistently evaluated along the crack tips of branched cracks. Although the present general formalism is dedicated to multiple straight crack segments with arbitrary geometries, such that the curved cracks can be suitably represented through a series of straight cracks, a particular attention is paid to kinked and forked cracks.

Stroh formalism and Green's functions for dislocation arrays in anisotropic bimaterials

As pictured in Fig. (1a), the bicrystals are formed by joining two linear anisotropic elastic materials A (upper material) and B (lower), which are defined by means of specific orientation relations along fixed crystal directions. The global and fixed coordinate system is conveniently defined by (O, x 1 , x 2 , x 3 ), the perfectly matched interface is located at x 2 = 0, with x 2 > 0 for material A, and x 2 < 0 for material B. Without loss of generality, the individual dislocations are embedded in the upper material A, and are aligned with the x 3 -axis. The prescribed Burgers vector b A is the same in magnitude and direction for all dislocations, and is defined in two dimensions without screw dislocation components, such that b A = b x 1 , b x 2 , 0 t , in accordance with the present mixed-mode loading at long distances. The field stresses σ(x 1 , x 2 ) = σ ij (x 1 , x 2 ) and the displacements u(x 1 , x 2 ) = u i (x 1 , x 2 ) that are produced by the internal displacement jump in both half-spaces A and B are related by the full anisotropic Hooke law for constitutive relations, as follows

σ ij (x 1 , x 2 ) = c ijkl (x 2 ) u k,l (x 1 , x 2 ) , (1)
where a comma stands for differentiation, with repeated indices denoting summation convention ranging from 1 to 3, unless stipulated otherwise. The anisotropic elastic constants of the fourth-order stiffness tensor are fully symmetric, i.e., c ijkl = c jikl = c ijlk = c klij , and the partial differential equation of mechanical equilibrium without body forces that is fulfilled in both crystals A and B in terms of the elastic displacement fields is given by

A c ijkl A u k,jl (x 1 , x 2 ) = 0 B c ijkl B u k,jl (x 1 , x 2 ) = 0 , (2)
where the elastic constants A c ijkl and B c ijkl are both defined by the local coordinates relative to the local crystal orientation and crystallography of each misoriented material. In the following, the pre-subscripts A and B in the elastic properties and also the field expressions will be omitted for clarity if no distinction between the dissimilar materials is required. The displacement field in the upper material A produced by n extrinsic dislocations from -∞ to ∞ is obtained by using the standard anisotropic elastic solution for a single dislocation and changing the dislocation location at (x D 1 , x D 2 ) by (x D 1 + nL, x D 2 ), as already formulated (Chou et al., 1975[START_REF] Chu | Elastic fields due to dislocation arrays in anisotropic bimaterials[END_REF][START_REF] Hirth | Stress fields of dislocation arrays at interfaces in bicrystals[END_REF][START_REF] Vattré | Three-dimensional interaction and movements of various dislocations in anisotropic bicrystals with semicoherent interfaces[END_REF]. Thus, the summation of the displacement fields of the individual dislocations in A leads to the displacement fields due to an infinite dislocation network parallel to the interface, A u array dis k (x 1 , x 2 ), with k = {1, 2, 3}, as follows A u array dis (x 1 , x 2 ;

x D 1 , x D 2 ) = 1 π ∞ ∑ n = -∞ Im A A ln(x 1 -(x D 1 + nL) + A p † (x 2 -x D 2 )) A q ∞ + 1 π ∞ ∑ n = -∞ Im 3 ∑ α = 1 A A ln(x 1 -(x D 1 + nL) + A p † x 2 -A p α * x D 2 ) A q α , (3)
where * indicates complex conjugation. The complex matrices f (z † ) with a logarithmic form are introduced, i.e.

f (z † ) = diag f (z 1 ), f (z 2 ), f (z 3 ) = f (z 1 ) f (z 2 ) f (z 3 ) , (4)
for which the symbol † is used to identify the element in the diagonal matrix. The first term in eq. ( 3) corresponds to the full-plane dislocation Green's functions in A, and both unknown constant vectors A q ∞ and A q α are determined by the perfectly matched interface conditions [START_REF] Ting | Anisotropic elasticity: theory and applications[END_REF], as defined by

A q ∞ = A H t b A A q α = A A -1 ( A M + B M * ) -1 ( B M * -A M * ) A A * = AA N I α A q ∞ * = AA N I α A H t * b A , (5)
where the diagonal matrices I α are defined by

I 1 = diag 1, 0, 0 , I 2 = diag 0, 1, 0 , and, I 3 = diag 0, 0, 1 , (6)
so that α = {1, 2, 3}. In eqs.

(3) and ( 5), the complex non-singular eigenmatrices A A = A a 1 , A a 2 , A a 3 are associated with the Stroh eigenvectors A a α as columns and the corresponding Stroh eigenvalues A p α are obtained by solving the following eigensystem of equations, i.e.

Q + p R + R t + p 2 T a = Π a = 0 , (7)
by virtue of the standard eigenvalue problem in anisotropic elasticity theory [START_REF] Stroh | Dislocations and Cracks in Anisotropic Elasticity[END_REF][START_REF] Stroh | Steady state problems in anisotropic elasticity[END_REF]. In eq. ( 7), the superscript t denotes the matrix transpose, while Q ik = c i1k1 , R ik = c i1k2 , and, T ik = c i2k2 , which are properly rotated with respect to the coordinate systems spanned by the x 1 -and x 2 -axes. A non-trivial solution can be found only if the determinant of Π is zero, i.e.

det Π = 0 , (8)
leading to a polynomial equation in p of degree six, with real coefficients. Due to the positive definiteness of elastic strain energy, the solutions of eq. ( 8) have six imaginary roots, which occur in complex conjugates [START_REF] Eshelby | Anisotropic elasticity with applications to dislocation theory[END_REF]. In particular, the roots are conveniently arranged such that the three first eigenvalue solutions p α have positive imaginary parts, indexed by superscripts α = {1, 2, 3}. The remaining three solutions have negative imaginary parts, so that p α+3 = p α * . The eigenvectors a α = a α k are also complex conjugates with a α+3 = a α * = a α k * , so that A A * in eq. ( 5) is defined by A A * = A a 1 * , A a 2 * , A a 3 * . Furthermore, the complex vectors h α of the subsidiary complex matrices A H = A h 1 , A h 2 , A h 3 in eq. ( 5) are related to the vectors a α by

p α h α = p α R t + p α T a α = -(Q + p α R) a α , (9)
while the eigenvectors a α and h β satisfy the central normalization relations in the Stroh formalism, as follows

(h β ) t a α + (a β ) t h α = δ αβ , (10)
with δ αβ the 3 × 3 Kronecker delta. Futhermore, the displacement field in material B is given by

B u array dis (x 1 , x 2 ; x D 1 , x D 2 ) = 1 π ∞ ∑ n = -∞ Im 3 ∑ α = 1 B A ln(x 1 -(x D 1 + nL) + B p † x 2 -A p α x D 2 ) B q α , (11)
within which the eigenmatrix B A is accordingly associated with the Stroh eigenvalues B p α , while the unknown vectors B q α for perfect interfaces are given by

B q α = B A -1 ( A M * + B M) -1 ( A M + A M * ) A A = BA N I α A q ∞ = BA N I α A H t b A , (12)
where the positive-definite Hermitian impedance tensors χ M in eqs. ( 5) and ( 12) are defined by

χ M = -i χ H χ A -1 , (13)
with χ = {A, B}, while χ h α = i χ M χ a α ,which offers a second relation between the displacement with χ a α and the traction with χ h α . Using index notation, the derivative of the displacement field components in A, given by eq. ( 3), with respect to the space coordinates is written as

A u array dis k,l (x 1 , x 2 ; x D 1 , x D 2 ) = 1 π ∞ ∑ n = -∞ Im 3 ∑ m = 1 A A km A q ∞ m (x 1 -(x D 1 + nL) + A p m (x 2 -x D 2 )) -1 (δ l1 + A p m δ l2 ) = + 1 π ∞ ∑ n = -∞ Im 3 ∑ m = 1 3 ∑ α = 1 A A km A q α m (x 1 -(x D 1 + nL) + A p m x 2 -A p α * x D 2 ) -1 (δ l1 + A p m δ l2 ) , (14)
which can be summed to give strains. The sum in eq. ( 14) over n from -∞ to ∞ has an explicit solution (Hirth et al., 1979, Morse and[START_REF] Morse | Methods of Theoretical Physics[END_REF]

, i.e., ∞ ∑ n=-∞ (z + nL) -1 = π L ctg π L z , (15)
so that eq. ( 14) can be expressed in a closed form as

A u array dis k,l (x 1 , x 2 ; x D 1 , x D 2 ) = 1 L Im A A km A q ∞ m ctg π L x 1 -x D 1 + A p m (x 2 -x D 2 ) (δ l1 + A p m δ l2 ) = + 1 L Im 3 ∑ α = 1 A A km A q α m ctg π L x 1 -x D 1 + A p m x 2 -A p α * x D 2 (δ l1 + A p m δ l2 ) , (16)
where summation over the repeated index m is applied from 1 to 3, while the displacement gradients in B are analogously given by

B u array dis k,l (x 1 , x 2 ; x D 1 , x D 2 ) = 1 L Im 3 ∑ α = 1 B A km B q α m ctg π L x 1 -x D 1 + B p m x 2 -A p α x D 2 (δ l1 + B p m δ l2 ) , (17)
by differentiating eq. ( 11) with respect to the space coordinates. By virtue of the constitutive Hooke law in eq. ( 1), the stress state in material A can also determined, as follows

A σ array dis ij (x 1 , x 2 ; x D 1 , x D 2 ) = 1 L Im A c ijk1 + A p m A c ijk2 A A km A q ∞ m ctg π L x 1 -x D 1 + A p m (x 2 -x D 2 ) = + 1 L Im 3 ∑ α = 1 A c ijk1 + A p m A c ijk2 A A km A q α m ctg π L x 1 -x D 1 + A p m x 2 -A p α * x D 2 , (18)
using eq. ( 1), while the stress fields in B are given by

B σ array dis ij (x 1 , x 2 ; x D 1 , x D 2 ) = 1 L Im 3 ∑ α = 1 B c ijk1 + B p m B c ijk2 B A km B q α m ctg π L x 1 -x D 1 + B p m x 2 -A p α x D 2 , (19)
which result from the presence of an infinite periodically spaced arrays of lattice dislocations with inter-dislocation spacings L embedded in the upper material A.

The singular integral equations for a network of kinked and forked cracks

The two-dimensional boundary-value problem for kinked and forked cracks in anisotropic bimaterials is treated by use of continuous distribution of dislocations with specific Burgers vectors along each crack segment to be determined by the tractionfree surface conditions [START_REF] Bilby | Continuous distributions of dislocations: New application of the methods of non-Riemannian geometry[END_REF][START_REF] Bilby | Dislocations and the Theory of Fracture[END_REF][START_REF] Bilby | The crack with a kinked tip[END_REF][START_REF] Hills | Solution of Crack Problems: The Distributed Dislocation Technique[END_REF]. As illustrated in Fig. (1a), the elementary branched crack, which is composed of N Max crack segments, with N Max = 2, and N Max > 2 for the kinked, and forked crack problem, respectively, and is also part of the infinite array of cracks is arbitrary oriented in the upper material A. Each rectilinear segment N , with N = 1, . . . , N Max , of the elementary branched crack is arbitrarily defined by θ N the oriented angle between the horizontal x 1 -axis of the global coordinate system (x 1 , x 2 ) and the local xN 1 -axis, which is collinear with the crack segment. A local Cartesian coordinate system ( xN 1 , xN 2 ) is also introduced and attached to each crack segment, centered at the point of coordinate c N 1 , c N 2 with half crack-length a N , subjected to an externally and uniformly applied stress field. These displacements are represented by unknown distributions of dislocation densities along each crack segment, which are determined by requiring that the traction-free surface conditions along the crack segments are fulfilled under plane strain conditions.

On the crack surfaces, the applied stress fields are therefore removed by the equivalent and opposite stresses produced by the sliding and opening displacement components from the crack arrays. The uniform load components σ ∞ ij applied at infinity are given by

σ ∞ ij =   σ ∞ 11 σ ∞ 12 0 σ ∞ 12 σ ∞ 22 0 0 0 0   , (20)
without loss of generality. The non-zero stress field components are conveniently expressed in the local axis of each crack segment N to

σ∞N 11 = σ ∞ 11 cos 2 θ N + σ ∞ 22 sin 2 θ N + 2 σ ∞ 12 sin θ N cos θ N , σ∞N 22 = σ ∞ 11 sin 2 θ N + σ ∞ 22 cos 2 θ N -2 σ ∞ 12 sin θ N cos θ N , σ∞N 12 = (σ ∞ 22 -σ ∞ 11 ) sin θ N cos θ N + σ ∞ 12 cos 2θ N , (21)
for which the superimposed symbol ˆrepresents any quantities transformed into the local coordinate system of the N th crack of interest. On the other hand, the corresponding stress field components produced by an array of extrinsic dislocations are given by introducing the vector A τ array dis (x 1 , x 2 ;

x D 1 , x D 2 ), as follows A τ array dis (x 1 , x 2 ; x D 1 , x D 2 ) = A σ array dis 11 (x 1 , x 2 ; x D 1 , x D 2 ), A σ array dis 22 (x 1 , x 2 ; x D 1 , x D 2 ), A σ array dis 12 (x 1 , x 2 ; x D 1 , x D 2 ) t , (22)
which, according to eqs. ( 18) and ( 5), is written in matrix notation as

A τ array dis (x 1 , x 2 ; x D 1 , x D 2 ) = 1 L Im A Λ t A A ctg π L x 1 -x D 1 + A p † (x 2 -x D 2 ) A H t b A = + 1 L Im 3 ∑ α = 1 A Λ t A A ctg π L x 1 -x D 1 + A p † x 2 -A p α * x D 2 AA N I α A H t * b A , (23)
where the elasticity-based matrix A Λ 3×3 is given by

A Λ = A Q 1k + A p k A R 1k , A R k2 + A p k A T 2k , A R k1 + A p k A T 1k , (24)
with k = {1, 2, 3}. Similarly to eq. ( 21), the shear and normal stress components produced by the dislocation networks are transformed from the global coordinate system to the local coordinate system dedicated to the specific N th crack segment, by use of the Mohr transformation matrix P N , i.e.

A σ array dis N

12

( xN 1 , xN

2 ; x D 1 , x D 2 ) A σ array dis N 22 ( xN 1 , xN 2 ; x D 1 , x D 2 ) = A ˜τ array dis N ( xN 1 , xN 2 ; x D 1 ,x D 2 ) = -sin θ N cos θ N sin θ N cos θ N cos 2θ N sin 2 θ N cos 2 θ N -2 sin θ N cos θ N = P N    A σ array dis 11 (x 1 , x 2 ; x D 1 , x D 2 ) A σ array dis 22 (x 1 , x 2 ; x D 1 , x D 2 ) A σ array dis 12 (x 1 , x 2 ; x D 1 , x D 2 )    = A τ array dis (x 1 ,x 2 ; x D 1 ,x D 2 ) , (25)
which are the only two components required to solve the branched crack problem. Here and in the following, the tilde symbol in subscripts represents the non-zero shear 12 and normal 22 components of the local traction stresses only, thus without the 11 component. Combining eq. ( 24) with eq. ( 25), the local traction stresses at

( xN 1 , xN 2 ) become A ˜τ array dis N ( xN 1 , xN 2 ; x D 1 , x D 2 ) = 1 L Im P N A Λ t A A ctg π L x 1 -x D 1 + A p † (x 2 -x D 2 ) A H t b A = + 1 L Im 3 ∑ α = 1 P N A Λ t A A ctg π L x 1 -x D 1 + A p † x 2 -A p α * x D 2 AA N I α A H t * b A , (26)
for which the discrete Burgers vectors b A are located at (x D 1 + nL, x D 2 ) in the upper material A. Following the concept of continuously distributed infinitesimal dislocations to analyze crack problems, the non-uniform displacements across surfaces are described in terms of continuous infinitesimal dislocations. The continuous distribution of dislocations is specified by a single-valued density function B( η) situated at point η, such that B( η)d η represents the total length of the Burgers vectors of the infinitesimal dislocations lying between the points η and η + d η on the crack plane. These artificial and mathematical dislocations with infinitesimal Burgers vectors db A = B A ( η)d η are also conveniently used in eq. ( 23) to picture the present branched cracks. For any particular point located along the N th crack segment among an arbitrary total number of segments N Max , the shear and normal traction across the crack segment N in the local system is also obtained by superposing the self-interaction of all infinitesimal dislocations along the segment as well as the contribution from other segments M N . In other words, these stress components at any point ( xN 1 , xN 2 = 0) on the N th crack segment are therefore due to the integrated effect from all distributed dislocation densities and can be expressed in the form of Fredholm integral equations, as follows

A ˜τ crack N ( xN 1 , xN 2 = 0) = 1 L ˆaN -a N Im P N A Λ t A A ctg π L x N 1 -x D N 1 + A p † (x N 2 -x D N 2 ) A H t B A ( ηN ) = + Im 3 ∑ α = 1 P N A Λ t A A ctg π L x N 1 -x D N 1 + A p † x N 2 -A p α * x D N 2 AA N I α A H t * B A ( ηN ) d ηN = + 1 L N Max ∑ M =1 M N ˆaM -a M 1 L Im P N A Λ t A A ctg π L x N 1 -x D M 1 + A p † (x N 2 -x D M 2 ) A H t B A ( ηM ) = + Im 3 ∑ α = 1 P N A Λ t A A ctg π L x N 1 -x D M 1 + A p † x N 2 -A p α * x D M 2 AA N I α A H t * B A ( ηM ) d ηM , 6 (27)
with respect to the local coordinates, where the parameters needed to compute the components of the diagonal matrices with cotangent terms are written with respect to coordinates of the crack segments, i.e.

x

N 1 = c N 1 + xN 1 cos θ N x N 2 = c N 2 + xN 1 sin θ N
, and,

         x D N 1 = c N 1 + ηN cos θ N x D N 2 = c N 2 + ηN sin θ N x D M 1 = c M 1 + ηM cos θ M x D M 2 = c M 2 + ηM sin θ M , (28)
merely derived from geometrical characteristics. By virtue of the superposition principle, the applied remote stresses are combined with the traction solutions produced by the continuously distributed dislocation densities to determine the resultant stress state, and the essential corresponding traction-free condition along all crack segment N , with N = 1, . . . , N Max , reads

∀ | xN 1 | < a N :                  0 = σ∞N=1 12 + A σ array N = 1 12 ( xN=1 1 , 0) 0 = σ∞N=1 22 + A σ array N = 1 22 ( xN=1 1 , 0) = . . . 0 = σ∞N=N Max 12 + A σ array N = N Max 12 ( xN=N Max 1 , 0) 0 = σ∞N=N Max 22 + A σ array N = N Max 22 ( xN=N Max 1 , 0) , (29)
where the components 12 and 22 induced by the dislocation densities on the right-hand side are given by the integral relations in eq. ( 27). The absence of traction vector components along all cracks leads to 2N Max equations, for 2N Max unknowns, for which the system of integral equations that cannot be analytically solved by inversion. Before applying the numerical integration treatment in the next section 2.3 to solve complex crack configurations, the following scaling rules for the continuous variables are conveniently introduced as follows

       t = xN 1 a N ŝ = ηN a N = ηM a M , (30)
where t ∈ [-1, 1] and ŝ ∈ [-1, 1] are dimensionless curvilinear coordinates along the crack segments. Furthermore, anticipating the inverse square root singularities at each end of each crack N , with N = 1, . . . , N Max , the unknowns Burgers vectors B A ( ηN ) in eq. ( 27) are also rewritten as B N A ( ŝ), and subsequently replaced by choosing the following interpolation functions, i.e.

B A ( ηN ) = B A ( ηN ( ŝ)) = B N A ( ŝ) = ω( ŝ) φ N ( ŝ) = 1 -ŝ2 -1/2 φ N ( ŝ) , (31)
where ω( ŝ) is the fundamental function that is continuously differentiable with one-sided derivative at the end points, while φ N ( ŝ) are nonsingular unknown functions that reflect the strength of singularities. In accordance with the Gauss-Chebyshev quadrature method, the interpolation functions are therefore related to the Chebyshev polynomials, such that eqs. (27)(28)(29)(30)(31) are also discretized by a series of algebraic relations to determine the suitable strength functions φ N ( ŝ) of the dislocation distributions for every crack segment N , as formulated in the following section.

Numerical solutions based on the Gauss-Chebyshev collocation method

The local displacements that are represented by the unknown distributions of dislocation densities B N A ( ŝ) along each crack segment, the unknown functions φ N ( ŝ) are determined by requiring that the traction-free state along the crack segments is fulfilled at a given and appropriate set of collocation points. These collocation points are are placed in between the integration points where the values of the sliding and opening displacement components are numerically computed. The total number of N Int integration points imposed on the N th crack branch and of N Int -1 collocation points are denoted by ŝi and tk , with i = 1, . . . , N Int , and k = 1, . . . , N Int -1, as defined by ŝi = cos

((i -1/2) π/N Int ) tk = cos (k π/N Int ) , (32)
respectively. The discrete collocation relations of the Fredholm-type integral equation with Cauchy kernel for the shear and normal traction across the crack segment N in eq. ( 27) are given in series form, ∀ k = 1, . . . , N Int -1, by

A ˜τ crack N ( tk , 0)

a N L NInt ∑ i = 1 W i ( ŝi ) Im P N A Λ t A A ctg π L a N ( tk -ŝi ) cos θ N + A p † a N ( tk -ŝi ) sin θ N A H t φ N ( ŝi ) = + Im 3 ∑ α = 1 P N A Λ t A A ctg π L a N ( tk -ŝi ) cos θ N + A p † (c N 2 + a N tk sin θ N ) -A p α * (c N 2 + a N ŝi sin θ N ) AA N I α A H t * φ N ( ŝi ) = + NMax ∑ M =1 M N a M L NInt ∑ i = 1 W i ( ŝi ) Im P N A Λ t A A ctg π L (c N 1 -c M 1 + a N tk cos θ N -a M ŝi cos θ M ) + A p † (c N 2 -c M 2 + a N tk sin θ N -a M ŝi sin θ M ) A H t φ M ( ŝi ) = + Im 3 ∑ α = 1 P N A Λ t A A ctg π L (c N 1 -c M 1 + a N tk cos θ N -a M ŝi cos θ M ) + A p † (c N 2 + a N tk sin θ N ) -A p α * (c M 2 + a M ŝi sin θ M ) AA N I α A H t * φ M ( ŝi ) , (33)
by use of eqs. ( 28) and (30)(31)(32), where W i ( ŝi ) = π/N Int is the appropriate weight function with respect to the fundamental function ω( ŝi ), as described by [START_REF] Erdogan | Numerical solution of singular integral equations[END_REF]. Furthermore, the corresponding net stress state on each segment N in eq. ( 33) is also reduced to a standard matrix equation, i.e.

∀ k = 1, . . . , N Int -1 :                  0 = σ∞N=1 12 + A σ crack N = 1 12 ( tk , 0) 0 = σ∞N=1 22 + A σ crack N = 1 22 ( tk , 0) = . . . 0 = σ∞N=N Max 12 + A σ crack N = N Max 12 ( tk , 0) 0 = σ∞N=N Max 22 + A σ crack N = N Max 22 ( tk , 0) , (34)
where the relevant external and crack-induced stress contributions are given by eqs. ( 21) and ( 33), respectively. The latter includes the dislocation density strengths φ N ( ŝi ) that are also evaluated at the N Int zeros of the Chebyshev polynomials of the first kind, while the crack tractions are specified at the N Int zeros of the Chebyshev polynomials of the second kind. The tractionfree boundary conditions provided by eq. ( 34) leads to an incomplete system of 2N Max (N Int -1) equations, for 2N Max N Int unknowns, so that 2N Max additional relations are needed to solve the system of equations. For a crack configuration with N segments, the specific unknowns are defined in the index notation by

∀ i = 1, . . . , N Int : φ 1 x1 ( ŝi ), φ 1 x2 ( ŝi ), . . . , φ N Max x1 ( ŝi ), φ N Max x2 ( ŝi ) , ( 35
) such that φ N ( ŝi ) = φ N x1 ( ŝi ), φ N x2 ( ŝi ), 0 t
, where φ N x1 ( ŝi ) and φ N x2 ( ŝi ) are equivalent to the sliding and opening of the material across the local coordinate axis x1 and x2 , respectively.

A consistent closure condition is used to prevent any perturbation of the external applied stress state at long range by the completely embedded cracks, formulated in the integral form as follows

N Max ∑ N = 1 ˆaN -a N B A ( ηN ) d ηN = 0 , (36)
so that the total Burgers vector content along the branched cracks is zero, i.e. no net displacement occurs. In the discretized scheme, eq. ( 36) reads

         N Max ∑ N =1 a N N Int ∑ i=1 W i ( ŝi ) φ N x1 ( ŝi ) = 0 N Max ∑ N =1 a N N Int ∑ i=1 W i ( ŝi ) φ N x2 ( ŝi ) = 0 , (37)
leading to two additional equations, after projection on the local coordinate axis x1 and x2 , respectively. For problems dedicated to periodically spaced arrays of single cracks, thus without branched segments, the system that combines together eqs. ( 34) with (37) contains enough linearly dependent equations to solve the unknown dislocation density distributions defined by eq. ( 31) with eq. ( 35). For N Int ≥ 2, however, extra equations are needed, and two types of continuity constraints that have to be fulfilled at crack kinking and branching are defined, respectively, such that the collocation integral equations and the additional constraints close the system of equations.

For the specific configuration of kinked cracks, with N Int = 2, the sliding and opening displacements are continuous at the kinking nodes between two adjacent crack segments, e.g., the principal segment indexed by N = 1 and the finite kinked segment by N = 2, for which the continuity condition is expressed as

lim η1 → a 1 B A ( η1 ) = lim η2 →-a 2 B A ( η2 ) ⇒ lim ŝ→ 1 B 1 A ( ŝ) = lim ŝ→ -1 B 2 A ( ŝ) , (38)
where the sign in ηN → ±a N and ŝ → ±1 depends on the left-and right-hand sides of the specific crack segment that is connected to the kink node. By virtue of the interpolation functions given by in eq. ( 31) in the normalized and discretized form, the equality of the dislocation densities at the kink locations from eq. ( 38) reads

     lim ŝ1 → 1 φ 1 x1 ( ŝ1 ) = lim ŝN Int →-1 φ 2 x1 ( ŝN Int ) lim ŝ1 → 1 φ 1 x2 ( ŝ1 ) = lim ŝN Int →-1 φ 2 x2 ( ŝN Int ) , (39)
which yields two extra boundary conditions after projection. From the square-root asymptotic field analysis near the crack tips, the unknown functions {φ N x1 , φ N x2 } remain finite at the end-points ŝ → ±1, as follows

lim ŝ1 → 1 φ 1 xj ( ŝ1 ) = 1 N Int N Int ∑ i=1 sin (π (2N Int -1) (2i -1) / (4N Int )) sin (π (2i -1) / (4N Int )) φ 1 xj (s i ) , lim ŝN Int →-1 φ 2 xj ( ŝN Int ) = 1 N Int N Int ∑ i=1 sin (π (2N Int -1) (2i -1) / (4N Int )) sin (π (2i -1) / (4N Int )) φ 2 xj (s N Int +1-i ) , (40)
concisely obtained by the Krenk interpolation formulae [START_REF] Hills | Solution of Crack Problems: The Distributed Dislocation Technique[END_REF][START_REF] Krenk | On the use of interpolation polynomials for solutions of singular integral equations[END_REF].

On the other hand, a stronger continuity constraint than the relations in eq. ( 39) is imposed at the branched nodes with multiple intersecting crack segments, when N > 2, by ensuring the amplitude of the singular dislocation density of the branch crack equal to zero at the point where all main and branched segments intersect, as suggested by Cleary and co-authors [START_REF] Barr | Thermoelastic fracture solutions using distributions of singular influence functions -Determining crack stress fields from dislocations distributions[END_REF]Cleary, 1981, Narendran andCleary, 1984) for application problems to rock fracture mechanics. Thus, eq. ( 39), for which the net opening and sliding displacements cancel out at crack barchning, becomes

     lim ŝ1 → 1 φ 1 x1 ( ŝ1 ) = lim ŝN Int →-1 φ 2 x1 ( ŝN Int ) = . . . = lim ŝN Int →-1 φ N Max -1 x1 ( ŝN Int ) = 0 lim ŝ1 → 1 φ 1 x2 ( ŝ1 ) = lim ŝN Int →-1 φ 2 x2 ( ŝN Int ) = . . . = lim ŝN Int →-1 φ N Max -1 x2 ( ŝN Int ) = 0 , (41)
leading to (2N Max -2) extra equations for the branched crack configurations. Combining the latter equations with the two relations in eq. ( 37) that guarantee the zero net Burgers vector condition as well as the 2N Max (N Int -1) relations from eqs. (34), the complete system that is composed of 2N Max linear equations can be solved and the corresponding unknowns in eq. ( 35) consistently be determined. Once the dislocation density distributions are known, the resultant stress fields σ array crack ij (x 1 , x 2 ) at any point in the upper material A is obtained as follows

A σ array crack (x 1 , x 2 ) σ ∞ + π LN Int NMax ∑ N = 1 a N NInt ∑ i = 1 Im A Λ t A A ctg π L x 1 -(c N 1 + a N ŝi cos θ N ) + A p † (x 2 -(c N 2 + a N ŝi sin θ N )) A H t φ N ( ŝi ) + Im 3 ∑ α = 1 A Λ t A A ctg π L x 1 -(c N 1 + a N ŝi cos θ N ) + A p † x 2 -A p α * (c N 2 + a N ŝi sin θ N ) AA N I α A H t * φ N ( ŝi ) , (42)
including the A σ array crack 33 (x 1 , x 2 ) component, where AA N and A Λ are defined in eqs. ( 5) and ( 24), respectively. Furthermore, the corresponding stress state in material B is given by

B σ array crack (x 1 , x 2 ) σ ∞ + π LN Int NMax ∑ N = 1 a N NInt ∑ i = 1 Im 3 ∑ α = 1 B Λ t B A ctg π L x 1 -(c N 1 + a N ŝi cos θ N ) + B p † x 2 -A p α (c N 2 + a N ŝi sin θ N ) BA N I α A H t φ N ( ŝi ) , (43)
where BA N is defined in eq. ( 12), and the elasticity-based matrix B Λ by

B Λ = B Q 1k + B p k B R 1k , B R k2 + B p k B T 2k , B R k1 + B p k B T 1k , (44) while
B Q ik = B c i1k1 , B R ik = B c i1k2
, and, B T ik = B c i2k2 . In particular, eq. ( 42) is of vital importance for computing the J k -integrals [START_REF] Budiansky | Conservation laws and energy-release rates[END_REF]Rice, 1973, Rice, 1968) and also predicting the crack growth of the different crack branches in anisotropic bimaterials.

Calculation of J k -integrals for branched crack arrays in anisotropic bimaterials

Path-independent relations are widely used to determine the two-dimensional J k -integral [START_REF] Budiansky | Conservation laws and energy-release rates[END_REF]Rice, 1973, Rice, 1968) as a criteria to predict crack growth direction [START_REF] Hellen | The calculation of stress intensity factors for combined tensile and shear loading[END_REF], where J 1 is the standard path-independent J-integral introduced by Rice (1968) that is related to the energy release rate per unit crack advance [START_REF] Gurtin | Configurational forces and the basic laws for crack propagation[END_REF]. On the other hand, J 2 is path-independent only by adding an extra nontrivial contribution on the crack faces [START_REF] Herrmann | On energy release rates for a plane crack[END_REF], and is also required to the mechanical fields near the crack tips. Basically, for a two-dimensional elastic body containing a crack that lies in the direction of x 1 -axis, the expression of the J k -integral components in the global coordinate system are

J Γ → 0 k = lim Γ → 0 ˛Γ A W (x 1 , x 2 )δ kj -A σ array crack ij (x 1 , x 2 ) A u array crack i,k (x 1 , x 2 ) = A L n j dΓ = lim Γ → 0 ˛Γ 1 2 A σ array crack mn (x 1 , x 2 ) A u array crack m,n (x 1 , x 2 ) n k -A σ array crack ij (x 1 , x 2 ) n j A u array crack i,k (x 1 , x 2 ) dΓ , ( 45
)
where A L is the Eshelby energy-momentum tensor [START_REF] Eshelby | The force on an elastic singularity[END_REF], A W (x 1 , x 2 ) is the strain energy density that is determined with respect of the Cartesian components of the stress and displacement gradient tensors from the previous section 2.3, and Γ is a contour in the (x 1 , x 2 ) plane that encloses the crack tips, while n j denotes the corresponding unit outward normal to Γ at (x 1 , x 2 ). The equality in eq. ( 45) holds only when the contour Γ shrinks onto the crack tips Γ → 0. However, if the contour Γ N for the N th crack segment does not tend to the tips, the integral relation in eq. ( 45) becomes

J Γ N k = ˆΓN ( 1 2 A σ array crack mn (x 1 , x 2 ) A u array crack m,n (x 1 , x 2 ) nk -A σ array crack ij (x 1 , x 2 ) n j ∂ x k A u array crack i (x 1 , x 2 ) = A u array crack i,k (x 1 ,x 2 ) dΓ - ˆΓ N A W (x 1 , x 2 ) + -n k dΓ (46)
where A W (x 1 , x 2 ) + -stands for the jump in the strain energy density across the crack faces ΓN , while the superscripts {-, +} refer to the two half-space subdomains connected to the negative and positive crack faces with respect to the x 2 -axis. The additional integral energy-based term on the right-hand side must be taken into account on the crack faces, which vanishes for J 1 , and also influences J 2 only. For kinked and forked cracks, eq. ( 46) is consistently transformed into the local coordinate system, as follows

ĴΓ N k = ˆΓN 1 2 A σ array crack mn (x 1 , x 2 ) A u array crack m,n (x 1 , x 2 ) nk -A σ array crack ij (x 1 , x 2 ) n j ∂ xk A u array crack i (x 1 , x 2 ) dΓ - ˆΓ N A W (x 1 , x 2 ) + -nk dΓ = -Ĵ2ex , ( 47
)
where nk = cos n j , xk is the projection of n j on the xk -axes. Furthermore, the gradient terms in eq. ( 47) are given by

∂ x1 A u array crack i (x 1 , x 2 ) ∂ x2 A u array crack i (x 1 , x 2 ) = cos θ N sin θ N -sin θ N cos θ N ∂ x 1 A u array crack i (x 1 , x 2 ) ∂ x 2 A u array crack i (x 1 , x 2 ) , (48)
with positive θ N , so that both J Γ N 1 -and J Γ N 2 -integrals are defined by

                   ĴΓ N 1 = ˆΓN 1 2 cos n j , x1 A σ array crack mn (x 1 , x 2 ) A u array crack m,n (x 1 , x 2 ) -A σ array crack ij (x 1 , x 2 ) n j cos θ N ∂ x 1 A u array crack i (x 1 , x 2 ) + sin θ N ∂ x 2 A u array crack i (x 1 , x 2 ) dΓ ĴΓ N 2 = ˆΓN 1 2 cos n j , x2 A σ array crack mn (x 1 , x 2 ) A u array crack m,n (x 1 , x 2 ) -A σ array crack ij (x 1 , x 2 ) n j -sin θ N ∂ x 1 A u array crack i (x 1 , x 2 ) + cos θ N ∂ x 2 A u array crack i (x 1 , x 2 ) dΓ + Ĵ2ex , ( 49
)
where Γ N is a circular path with suitable polar coordinates centered at crack tips of the N th segment, while ΓN in the integral term Ĵ2ex from eq. ( 47) is a convenient rectangular contour along the crack lips of width r 0 , as schematically shown in Fig. (1b). The latter is introduced in the calculation of ĴΓ N 2 only, and corresponds to the vicinity of dislocation cores for which the present linear elasticity theory fails.

Numerical examples

Various application examples from the single kinked crack in isotropic homogeneous materials to the closely-spaced network of forked cracks in anisotropic bimaterials as well as some effects from the inter-crack spacings to the elastic heterogeneity are investigated. The former limiting case of single branched cracks is theoretically formulated and the corresponding asymptotic solutions are compared to existing results reported in the literature. The present solutions are numerically evaluated in terms of stress intensity factors K and the path-independent Ĵk -integrals using N Int = 256 integration points, which give excellent results compared to the mathematical crack solutions from [START_REF] Kitagawa | Crack-Morphological Aspects in Fracture Mechanics[END_REF]. The numerical accuracy of the the shear and normal tractions with the Gauss-Chebyshev quadrature scheme ranges from ten to twelve significant digits for both pure tensile and shear conditions at long range. Application examples are performed on Copper, which is moderately anisotropic, where the elastic constants are c 11 = 168.4 GPa, c 12 = 121.4 GPa, and c 44 = 75.4 GPa, indexed in Voigt notation, while the dislocation core width is arbitrary given by the lattice parameter of Copper as r 0 = 0.3615 nm. Unless mentioned otherwise, the Kand Ĵk -based quantities are normalized by K I = σ ∞ 0 √ πa 1 and J I = K 2 I (1 -ν 2 iso )/E iso , where σ ∞ 0 is the constant tensile or shear stress applied at infinity, ν iso is the Poisson ratio and E iso is Young modulus, obtained using the Voigt averaging procedure.

The limiting problem of single kinked and forked cracks

Several theoretical studies have investigated the implications of crack kinking and forking to the fracture behavior of engineering materials, specially examined on the basis of elastic stress intensity factors. However, considerable discrepancies of the stress intensity factor solutions have been recognized in the literature, as stipulated by as mentioned by [START_REF] Bilby | Stress intensity factors at the tips of kinked and forked cracks[END_REF], [START_REF] Suresh | Plastic near-tip fields for branched cracks[END_REF], and [START_REF] He | Kinking of a crack out of an interface[END_REF], for instance. The results for a single crack in an isotropic and infinite medium from [START_REF] Bilby | Stress intensity factors at the tips of kinked and forked cracks[END_REF], using a plane theory of elasticity derived by [START_REF] Khrapkov | The first basic problem for a notch at the apex of an infinite wedge[END_REF], are widely referenced by the fracture community.

Without completely repeating the entire derivation of the equations from the previous section 2, the unknown functions in eq. ( 35) are reduced to 4N Int (to 6N Int ) for the kinked * (forked * * ) crack problems, which can also be determined by solving 4N Int -4 (6N Int -6) linear equations from eq. ( 34), respectively, as follows

∀ k = 1, . . . , N Int -1 : * *                      *          0 = σ∞N=1 12 + A σ crack N = 1 12 ( tk , 0) 0 = σ∞N=1 22 + A σ crack N = 1 22 ( tk , 0) 0 = σ∞N=2 12 + A σ crack N = 2 12 ( tk , 0) 0 = σ∞N=2 22 + A σ crack N = 2 22 ( tk , 0) , 0 = σ∞N=3 12 + A σ crack N = 3 12 ( tk , 0) 0 = σ∞N=3 22 + A σ crack N = 3 22 ( tk , 0) , (50)
combined with 2 additional equations given by eq. ( 37), plus 2 (4) more relations by eqs. ( 39) and ( 41) with the aid of eq. ( 40) to complete the system * ( * * ). The limiting problem of single kinked and forked cracks is also analyzed by use of the standard limit when L → ∞, i.e.

lim L → ∞ 1 L ctg π L z = (πz) -1 , (51)
such that the crack-induced stress fields in eq. ( 50) are obtained by reducing eq. ( 33) to

A ˜τ single crack N ( tk , 0)

a N N Int NInt ∑ i = 1 Im P N A Λ t A A a N ( tk -ŝi ) cos θ N + A p † a N ( tk -ŝi ) sin θ N -1 A H t φ N ( ŝi ) = + Im 3 ∑ α = 1 P N A Λ t A A a N ( tk -ŝi ) cos θ N + A p † (c N 2 + a N tk sin θ N ) -A p α * (c N 2 + a N ŝi sin θ N ) -1 AA N I α A H t * φ N ( ŝi ) = + 1 N Int NMax ∑ M =1 M N a M NInt ∑ i = 1 Im P N A Λ t A A (c N 1 -c M 1 + a N tk cos θ N -a M ŝi cos θ M ) + A p † (c N 2 -c M 2 + a N tk sin θ N -a M ŝi sin θ M ) -1 A H t φ M ( ŝi ) = + Im 3 ∑ α = 1 P N A Λ t A A (c N 1 -c M 1 + a N tk cos θ N -a M ŝi cos θ M ) + A p † (c N 2 + a N tk sin θ N ) -A p α * (c M 2 + a M ŝi sin θ M ) -1 AA N I α A H t * φ M ( ŝi ) , (52)
with N Max = 2 ( * ) or = 3 ( * * ), respectively, for which the associated stress components 12 and 22 are arranged as the term on the left-hand side of eq. ( 25). Once the system is solved and the corresponding vectorial unknowns {φ N ( ŝi )}, with i = 1, . . . , N Int , are determined, eq. ( 42) is used to calculate the mixed-mode stress intensity factors A KN in the suitable local coordinate system affiliated to the tips of the N th crack, as follows

A KN =    A KN 2 A KN 1 A KN 3    = lim x1 →±a N 2π (±a N -x1 )    cos θ N sin θ N 0 -sin θ N cos θ N 0 0 0 1    = Ω N 3×3    A σ single crack N 12 ( xN 1 , 0) A σ single crack N 22 ( xN 1 , 0) A σ single crack N 32 ( xN 1 , 0) = 0    = A τ single array N ( t,0) = ± 2πa N Ω N 3×3 lim t→±1 1 ∓ t A τ single array N ( t, 0) , (53)
thus A KN 3 = 0, in the present work. By applying eq. ( 51) to eq. ( 18), the shear and normal stress components produced by a single dislocation are also given by

A σ array dis i2 (x 1 , x 2 ; x D 1 , x D 2 ) = 1 π Im x 1 -x D 1 + A p m (x 2 -x D 2 ) -1 ( A R ki + A p m A T ik ) A A km A q ∞ m = + 1 π Im 3 ∑ α = 1 x 1 -x D 1 + A p m x 2 -A p α * x D 2 -1 ( A Q ik + A p m A T ik ) A A km A q α m , (54)
so that the curly bracketed part in eq. ( 53) is obtained by considering the integrated effect from all distributed dislocation densities, and also including eqs. ( 5), ( 9) and (31) into eq. ( 54). For the crack tip on the right-hand side of the N th segment, i.e. for positive t, the bracketed term in eq. ( 53) reads

lim t→1 1 -t A τ single crack N ( t) = 1 π Im A H A H t lim t→1 1 - t ˆ1 -1 B N A ( ŝ) t - ŝ d ŝ = + 1 π Im 3 ∑ α = 1 A H AA N I α A H t * lim t→1 1 - t ˆ1 -1 B N A ( ŝ) t - ŝ d ŝ , (55)
within which the limits of the integral fonctions with Cauchy kernels can be explicitly evaluated [START_REF] Hills | Solution of Crack Problems: The Distributed Dislocation Technique[END_REF][START_REF] Huang | Mixed-mode stress intensity factors for cracks located at or parallel to the interface in bimaterial half planes[END_REF][START_REF] Yang | Kinked crack in anisotropic bodies[END_REF], as follows

lim t→1 1 - t ˆ1 -1 B N A ( ŝ) t - ŝ d ŝ = π √ 2 φ N (+1) , (56)
with φ N (+1) = φ N x1 (+1), φ N x2 (
+1), 0 t . Similar expression for the left-hand side of the tips are similarly derived for negative t, so that eq. ( 53) reads

A KN =    A KN 2 A KN 1 A KN 3 = 0    = ± √ πa N Ω N 3×3 Im A H A H t φ N (±1) + 3 ∑ α = 1 A H AA N I α A H t * φ N (±1) , (57)
where the finite values at the end-points, i.e.

φ N (±1) = φ N x1 (±1), φ N x2 (±1), 0 t , (58)
are obtained by the Krenk interpolation formulas in eq. ( 40). [START_REF] Kitagawa | Crack-Morphological Aspects in Fracture Mechanics[END_REF] present work [START_REF] Kitagawa | Crack-Morphological Aspects in Fracture Mechanics[END_REF] present work with the induced (a) kinked and (b) forked angle θ 2 under uniform tension, and comparison with the results from [START_REF] Kitagawa | Crack-Morphological Aspects in Fracture Mechanics[END_REF] for a 2 /a 1 = 0.1, while K I is the nominal stress intensity factor.

A K 2 / K I ^2θ2 2a1 2a2 θ2 2a1 2a2 a) b) N = 2 A K 2 / K I ^N = 2 A K 1 / K I ^N = 2 A K 1 / K I ^N = 2

Validation with existing results reported in the literature

The limiting case for single kinked and branched cracks in two-dimensional elastically homogeneous and isotropic materials is investigated for validation purposes, for which the isotropic elastic constants are obtained using the Voigt averaging procedure. Under uniaxial uniform traction, the mixed-mode stress intensity factors are normalized by the corresponding nominal stress intensity factor

K I = σ ∞ 22 √
π for a straight crack, where the crack length is projected in the direction perpendicular to the tensile x 2 -axis, so that = a 1 + a 2 cos θ 2 , with θ 1 = 0 • and a 2 /a 1 = 0.1. The variations of the normalized stress intensity factors 2a) and (b), respectively. The present solutions are in excellent agreement with the results introduced by [START_REF] Kitagawa | Crack-Morphological Aspects in Fracture Mechanics[END_REF], also presented by [START_REF] Suresh | Plastic near-tip fields for branched cracks[END_REF]. In particular, the maximum value of A KN=2 1 for the kinked cracks is found using isotropic elasticity at θ 2 = 60.1 • , while the in-plane shear mode vanishes at 2θ 2 = 28.4 • , as similarly obtained by [START_REF] Bilby | Stress intensity factors at the tips of kinked and forked cracks[END_REF].

Effect of the anisotropic elasticity under pure tensile and shear conditions

The previous configuration of kinked cracks, with θ 1 = 0 • and a 2 /a 1 = 0.1, is used to analyze the effect of anisotropic elasticity in the field solutions, as focused in Fig. (3). Here and in the following, the superscripts N and subscripts A in the Kand Ĵk -based quantities will be omitted for clarity in the figure legends since the schematics explicitly illustrate where these tip quantities are evaluated. The solid (dotted) lines represent the Kand Ĵk -based results obtained by the anisotropic (isotropic) elasticity theory under pure traction and shear loadings at long range.

The differences between isotropic and anisotropic elasticity under traction occur for large values of θ 2 , i.e., θ 2 > 40 • and θ 2 > 25 • for K1 and K2 , respectively, as illustrated in Fig. (3a). The maximum value of K2 using anisotropic elasticity corresponds to θ 2 = 69.1 • , which also differs from the isotropic calculation with θ 2 = 60.1 • . On the other hand, Fig. (3b) shows that the stress intensity factor K2 for pure shear possesses two maximum values in magnitude with opposite sign, which reveal the crackopening stress state under tension for negative values of θ 2 . The anisotropic elasticity tends to decrease the magnitudes of K1 and K2 for the kinked cracks under shear loads, especially for the large kinked angles, when |θ 2 | > 30 • .

The Ĵk -integrals in Fig. (3c) and (d) exhibit more significant quantitative and qualitative differences between isotropic and anisotropic results than the stress intensity factors K. For instance, the maximum magnitude in the J 1 -integrals for pure shear occurs at Ĵ1max = ±59.4 • (= ±58.9 • ) using anisotropic (isotropic) elasticity, which differs from the commonly used zero-K2 criterion for finite kinked cracks since K2 = 0 occurs at θ 2 = ±66.9 • (±70.3 • ). Interestingly, Ĵ 1 (θ 2 = 0 • ) = 0 and Ĵ 1 (θ 2 = 0 • ) < 0, for the anisotropic case, where denotes the derivative with respect to θ 2 , while Ĵ 1 (θ 2 = 0 • ) = 0 and Ĵ 1 (θ 2 = 0 • ) > 0, for the isotropic case under shear loads.

Overall, the extra Ĵ2ex terms in the Ĵk -integrals, as developped in section 2.4, are non-zero for all presented crack configurations. The Ĵ2 -integrals also play a critical role in the Ĵ-integral evaluations [START_REF] Lee | Contradiction of J k integral[END_REF], as defined by Ĵ = ( Ĵ2 1 + Ĵ2 2 ) -1/2 , while the corresponding maximum magnitude of Ĵmax is found to occur when θ 2 = 21.6 • (= 16.9 • ) using the anisotropic (isotropic) elasticity calculation under traction, as pictured by * in Fig. (3c), while Ĵ1max occurs at θ 2 = 0 • . The criterion based on the Ĵ-integrals also means that the actual path is not a pure mode-I path, which may also have significant consequences in predicting of the crack stability as well as the directions of crack propagation as the compared with the well-know Ĵ1 -integrals in relation with the corresponding maximum energy-release rates for cracks in anisotropic solids. Under shear loads, Ĵmax = ±42.6 • and ±37.4 • using anisotropic and isotropic solutions, respectively, exhibiting further discrepancies performed by use of the isotropic approximation in the elastic material properties.

According to the aforementioned fracture criteria, the most significant dissimilarity among the predictions occurs when the loading is related to pure shears, for which the zero-K2 criterion as a fracture criterion for finite kinked cracks in anisotropic solids is very questionable.

Effect of the crack length ratios

Figure ( 4) illustrates the influence of the a 2 /a 1 ratio on the stress intensity factors K as well as the Ĵk -integral calculations. The previous results with a 2 /a 1 = 0.1 are pictured with dotted lines for comparison, while the field solutions for a 2 /a 1 = 1 are J 1 iso J 2 iso J 2 ex iso J iso aniso aniso aniso aniso (4b) shows that the zero-K2 criterion is satisfied for θ 2 = ±50.4 • , which can also be compared to ±66.9 • for a 2 /a 1 = 0.1. Such significant difference demonstrates the great importance in considering the finite size of kinked cracks in appropriate crack problems as well as the possible hazards that can emerge when comparing the present results with solutions from asymptotic and perturbation analysis. According to Fig. (4c) and (d), the maximum Ĵ1 -integral magnitudes are reached when θ 2 = 0 • for both external loads using anisotropic elasticity theory, which also dramatically differ from the previous prediction based on the zero-K2 evaluation for the pure shear case. Clearly, such further discrepancy indicates that the K-based criteria are debatable as the fracture criteria for finite-length branched cracks. On the other hand, the previous comment on the discrepancies resulting from the Ĵversus Ĵ1integral criteria is even more evident for the case with equal lengths in both traction or shear loads, mainly due to the effect of Ĵ2 -terms that exhibit more oscillating characteristics for a 2 /a 1 = 1 with respect to θ 2 than the solutions with a 2 /a 1 = 0.1.

Stress intensity factors K/ K I K 1 iso K 2 iso K 1 aniso K 2 aniso 0 θ 2 (

Pure traction

Effect of the inter-crack spacings

The second size effect in Fig. ( 5a) and (b) is related to the influence on the Ĵk -integrals of the inter-crack spacings L of the kinked crack arrays with a 2 /a 1 = 0.1, under pure tensile and shear conditions, respectively. The dotted lines represent the calculations for single kinked cracks in anisotropic materials, thus for L → ∞, while the corresponding arrays of closely-spaced cracks with L = 5a 1 /2 and θ 1 = 0 • are depicted by solid lines. For information, the particular crack configuration with L = 2a 1 and θ 2 0 • , corresponds to a crack network that cuts horizontally the entire solids. The amplitude of both Ĵ1 -and Ĵ2 -integrals increases with decreasing L for both loading conditions. While the trends of the variations are qualitatively unchanged for the pure traction case in Fig. (5a), the influence of L is more noticeable for the shear loads. For instance, Ĵ1max and Ĵmax occur at θ 2 ± 48.0 • and ±36.9 • for the network of cracks, compared to ±59.4 • and ±42.6 • for the single kinked cracks in anisotropic material, respectively. Without of need for extensive statistical aspects, this significant difference shows that the present ideal elastic interactions between the cracks at short range are also relevant to take into account for applications to crack coalescence and more generally, to the problem of effective elastic properties of anisotropic solids containing specific crack arrangements and densities with mixed-mode interactions. Furthermore, Ĵ 1 (θ 2 = 0 • ) = 0 and Ĵ 1 (θ 2 = 0 • ) > 0 for the network of cracks under shear loads, which exhibits a change of sign of Ĵ 1 compared with the corresponding single crack case.

Effect of heterogeneity on the elastic properties

Figure ( 6) illustrates the influence of elastic interaction that arises from the mismatch of the elastic properties on the Ĵkintegrals for both the single kinked crack configuration as well as the infinite arrays of kinked cracks under pure shear loads. As 13 = a 1 /5, for which the tip of the kinked segment is therefore close to the interface when θ 2 = -90 • . The elastic properties in the upper material A are fixed, while the elastic constants of the lower material B are fictitiously multiplied or divided by three. In comparison with the previous solutions using homogeneous elasticity, as pictured by dotted lines, the Ĵk profiles are asymmetrical with respect to θ 2 = 0 • , as shown by solid lines, especially for the Ĵ1 -integrals. In accordance with Fig. (4b), for which the valid tensile opening mode with K1 > 0 is related to θ 2 < 0 • , thus with a kinked segment taken downward towards material B, the driving force by evaluating Ĵ1 at the kinked tips for the single crack is larger for the case where material B is softer, thus in turn, when material A is stiffer, by comparing Fig.

0 1 K 1 a 2 /a 1 = 0.1 K 2 a 2 /a 1 = 0.1 K 1 a 2 /a 1 = 1 K 2 a 2 /a 1 = 1 J 1 a 2 /a 1 = 0.1 J 2 a 2 /a 1 = 0.1 J 2 ex a 2 /a 1 = 0.1 J a 2 /a 1 = 0.1 a 2 /a 1 = 1 a 2 /a 1 = 1 a 2 /a 1 = 1 a 2 /a 1 = 1 ^* Fi

The linear elastic problem of single and closely-spaced forked cracks

As illustrated by the schematics in Fig. ( 7), the single forked crack in the homogeneous anisotropic material is composed of three finite-length segments, for which two segments are fixed, i.e., crack 1 and crack 3, arbitrarily oriented with θ 1 = 0 • and θ 3 = 60 • , while θ 2 varies from -120 • to 20 • . Furthermore, a 1 and a 2 are fixed as well, with a 2 /a 1 = 0.1, and different configurations with three lengths for a 3 are investigated, i.e., a 3 = a 2 /2 (case 1), a 3 = a 2 (case 2), and a 3 = 2a 2 (case 3). Figure (7a) shows that K1 cancels at θ 2 = -92.4 • for crack 2 and becomes negative below for the three cases. When θ 2 > -60 • , K1 becomes increasingly sensitive as crack 2 approaches crack 3, while the magnitude is larger as the length a 3 is smaller. Furthermore, K2 varies quadratically with θ 2 for crack 2, until a change of sign for larger angles. From all results in terms of the stress intensity factors on crack 3, K1 becomes negative for case 3 of the smallest a 3 when θ 2 > -8.4 • , while the other K-based quantities vary more monotonically, as illustrated Fig. (7b). Figure (7c) shows that the closer crack 2 is to crack 3, the larger the magnitude of Ĵ1 , while Ĵ2 reaches a maximum magnitude at θ 2 = -30.9 • , -35.5 • , and -39.6 • , for case 1, 2, and 3, respectively, before becoming negative for θ 2 = 1.7 • , -1.5 • , -9.3 • , respectively. Interestingly, both Ĵ1 -and Ĵ2integrals for crack 3 behave oppositely in sign until the magnitudes converge to almost zero for the largest of θ 2 , for instance when θ 2 > -14.1 • for case 1, as depicted by Fig. (7d). The most advanced application example of the present formalism is introduced in Fig. ( 8), illustrating the effects of both intercrack spacings and elastic mismatch on the Ĵk -integrals at two extreme tips of forked crack 2 and crack 3 in homogeneous and heterogeneous materials under traction. The specific case 2 from Fig. ( 7) is treated, thus with a 2 /a 1 = a 2 /a 3 and θ 2 = 60 • , while θ 2 varies from -120 • to 20 • . For comparison, the dotted lines are associated with simplified crack configurations. According 14

J 1 L = ∞ J 2 L = ∞ J 2 ex L = ∞ J L = ∞ L = 50 L = 50 L = 50 L = 50
^F igure 5: Effect of the inter-crack spacings on the Ĵk -integral evaluation for kinked cracks under pure tensile and shear conditions. The single kinked crack is depicted with dotted lines, while the closely-spaced network of kinked cracks is illustrated using solid lines.

J case 1 J 1 case 2

J 2 case 1 J 1 case 1 J case 2 J 2 case
K 2 case 1 K 1 case 2 K 1 case 1 K 2 case 2 K 1 case 3 K 2 case 3 J 1 case 3 J case 3 J 2 case 3 σ 22 ∞ crack 1 σ 22 ∞ θ 3 = 60° (fixed) 2a 1 2a 3 θ 2 2a 2 crack 3 θ 3 = 60° (fixed) 2a 1 2a 3 θ 2 2a 2 K 1 K 2 K 2 K 1 crack 2
crack 1 ^F igure 7: Variation of the stress intensity factors K and Ĵk -integrals for various configurations of single forked cracks. Three cases with different crack lengths are investigated, while θ 2 is the varying angle of crack 2, as pictured in the schematics. and discussed with respect to the predictions indicated by the crack-tip K-based criteria for anisotropic solids. The influence played by the short-range interactions, anisotropic elasticity, elastic mismatch, applied stress direction, inter-crack spacings and crack length ratios of various kinked and forked crack geometries are investigated in the view to revealing significant differences between both classical J k -and Kbased predictions in facture mechanics.

Extensions of the present work to the fundamental case of interfacial cracks interacting with free surfaces using appropriate field solutions of intrinsic dislocation networks [START_REF] Vattré | Mechanical interactions between semicoherent heterophase interfaces and free surfaces in crystalline bilayers[END_REF][START_REF] Vattré | Elastic interactions between interface dislocations and internal stresses in finite-thickness nanolayered materials[END_REF] is left to follow-on analyses.

Figure 1 :

 1 Figure1: (a) An infinite periodically spaced network of branched cracks with three segments is embedded in the upper half space of an anisotropic bimaterial in presence of a perfectly bounded internal interface. The mismatch between both anisotropic elastic constants is taken into account, such that discontinuous stress components can be obtained by the present field solutions. (b) Definition of integration contours around a crack tip for evaluation of the path-independent J k -integrals.

Figure

  Figure 2: Variation of normalized A KN=2 1

  the kinked and forked angle θ 2 are shown in Fig. (

  Effect of the anisotropic elasticity on the stress intensity factors K and Ĵk -integrals for kinked cracks under pure tensile and shear conditions. depicted with solid lines. When θ 2 < 79.5 • , K1 becomes negative for the pure traction case with a 2 /a 1 = 1, as pictured by * in Fig.(4a), meaning that the kinked crack tip closes and values of K1 for greater values of θ 2 > 79.5 • are invalid. The shear loading case in Fig.

gure 4 :

 4 Effect of of the crack length ratio on the stress intensity factors K and Ĵk -integrals for kinked cracks under pure tensile and shear conditions. illustrated in Fig.(1a), the elastic mismatch leads to a discontinuity in the in-plane stresses, for which the magnitude increases by increasing material mismatch of the anisotropic bimaterials. The results are presented for θ 1 = 0 • and a 2 /a 1 = 0.1, as previously, while c N = 2 2

  (6a) and (b). The kinked crack part is also attracted to the soft materials. For elastically soft material B, the Ĵ1 and Ĵ components reach maximum values in Fig. (6b) when the kinked segment is located at θ 2 = -73.3 • and -56.8 • , respectively, while Ĵ1max and Ĵmax occur at θ 2 = -58.6 • and -46.1 • for an infinity network with closely spaced cracks, as shown in Fig. (6c).

Figure 9 :

 9 Figure 9: Dimensionless stress field solutions produced by a network of equally-and closely-spaced forked cracks in an anisotropic bimaterial under traction that corresponds to the particular configuration * in Fig. (8c). (a) Contours of non-zero and dimensionless stress components. (b) Stress profiles along the vertical x 2 -axis, as depicted by the dotted lines in (a). The traction-free conditions on the crack planes are also satisfied.

Fi gure 6: Effect of the heterogeneous elasticity on the Ĵk -integral evaluation for a single kinked crack as well an infinite network of kinked cracks under pure shear loads. The elastic mismatch is introduced by fictitiously multiplied or divided by three the elastic constants of the lower material B.

to Fig. (8a), the closer crack 2 is to crack 3, the higher the Ĵ1 amplitude of the former crack 2, so that crack 3 enhances the propagation of crack 2 for the largest values of θ 2 . This situation is more pronounced in the case of crack network. In particular, the Ĵk -based quantities for crack 2 are small in magnitude when crack 2 is diametrically opposed to the former crack 3. On the contrary, Fig. (8b) shows that the Ĵk -based quantities for crack 3 tend towards zero when the crack 2 is closest to the fixed crack 3, which corresponds to a stress state of crack 3 that is partly shielded by the elastic stress field produced by the crack 2. Figures (8c) and (d) demonstrate that the elastic heterogeneity only plays a role on crack 2, while the elastic mismatch has also no effect on the behavior of the crack 3. Interestingly, the elastic mismatch leads to the presence of a zero derivative for Ĵ1 , resulting in the maximum value of Ĵ1max at θ 2 = -18.9 • for crack 2, as pictured by * in Fig. (8c). This specific state does not exist in the homogeneous elasticity context, and the corresponding non-zero and dimensionless stress components, defined by 9a) and (b) depict the large discontinuities of the in-plane stress component 11 across both crack and interface planes as well as the traction-free conditions for 22 and 12 along the main crack plane that are therefore fully satisfied, as required.

Concluding remarks

Using the Stroh formalism to the anisotropic elasticity theory of extrinsic dislocations, the two-dimensional fracture problem of multiple branched crack arrays in anisotropic bimaterials is formulated by means of coupled integral equations. These equations are related to the arbitrarily-oriented configurations of infinitely periodic cracks in dissimilar orthotropic half-spaces under arbitrary far-field stress loading conditions, combined with appropriate boundary conditions in terms of dislocation density distributions along the crack segments. The full-field solutions for kinked and forked crack arrays are solved by employing the Gauss-Chebyshev quadrature and collocation formulae on Chebyshev nodes, while the limiting case of individual cracks in homogeneous and isotropic materials is theoretically derived for validation and comparison purposes. Explicit expressions of the local stress intensity factors for the single branched cracks in anisotropic bimaterials are obtained by interpolation, for which the results of the asymptotic limiting case are in excellent agreement with existing solutions reported in the literature. The path-independent J k -integrals as crack propagation criterion are subsequently evaluated for mixed-mode crack configurations