Discrete crack dynamics: A planar model of crack propagation and crack-inclusion interactions in brittle materials

The Multipole Method (MPM) is used to simulate the many-body self-consistent problem of interacting elliptical micro-cracks and inclusions in single crystals. A criterion is employed to determine the crack propagation path based on the stress distribution; the evolution of individual micro-cracks and their interactions with existing cracks and inclusions is then predicted using what we coin the Discrete Crack Dynamics (DCD) method. DCD is fast (semi-analytical) and particularly suitable for the simulation of evolving low-speed crack networks in brittle or quasi-brittle materials. The method is validated against finite element analysis predictions and previously published experimental data. © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The path followed by a slowly propagating crack depends on the microstructure of the material, comprising features such as inclusions, grain boundaries and other cracks. Several techniques have been developed to predict crack paths. The finite element method (FEM) is one of the most widely used tools for studying stationary cracks and fracture along a pre-defined path ( Strang and Fix, 1973;Zienkiewicz et al., 1977;Reddy, 1993 ), but localised features such as cracks, holes, and inclusions are not efficiently resolved by mesh refinement, and general path crack growth cannot be modelled accurately and efficiently using conventional FEM, which in part led to the development of the extended finite element method (XFEM) by Dolbow and Belytschko (1999) , Belytschko and Black (1999) and, Belytschko et al. (2009) . Another alternative to FEM is the boundary element method (BEM), also referred to as the boundary integral equation method (BIEM) ( Aliabadi and Rooke, 1991;Mogilevskaya and Linkov, 1998 ). The main attraction of BEM is the reduction of the spatial dimensions of the problem since only the boundary needs to be discretised. Although the simulation of crack growth including interactions with inclusions is possible with the XFEM and BEM methods ( Loehnert and Belytschko, 2007;Belytschko et al., 2008;Natarajan, 2011;Natarajan et al., 2014 ), these methods become computationally ex- * Corresponding author. pensive when the number of localised features is increased. Therefore, there is scope to develop more efficient tools to predict the evolution of complex crack/defect networks. In this paper, a semianalytical method is proposed that fulfils this objective.
Discrete Crack Dynamics (DCD) was developed to predict the evolution of crack networks in heterogeneous solids containing elliptical inclusions. Similar to other computational methods, such as Molecular Dynamics or Discrete Dislocation Dynamics ( Alder and Wainwright, 1959;Cleveringa et al., 1997;Van der Giessen and Needleman, 1995 ), the validity of the calculations relies upon resolving the elastic interactions of defects in the system. In DCD, elliptical inclusions and cracks are the basic constituent elements of the system, and play a role analogous to that played by dislocations in discrete dislocation dynamics.
Micro-crack and inclusion interactions have been studied using analytical and numerical methods before: notably the Kachanov method ( Kachanov, 1985;Chudnovsky, 1983;Kachanov, 1993;Kachanov et al., 1990;Kachanov and Montagut, 1986;Kachanov, 1986;Sevostianov and Kachanov, 2010 ), the modified Kachanov method for closely interacting cracks ( Gorbatikh et al.,20 07;Li et al.,20 03 ), the polynomial approximation method ( Horii and Nemat-Nasser, 1985 ), Maxwell's method ( Maxwell, 1881;McCartney and Kelly, 2008;McCartney, 2010;Mogilevskaya and Crouch, 2013;Mogilevskaya et al., 2010 ) and the Multipole method (MPM) ( Kushch, 2013a;2013b;1998a;1998b;Kushch and Sangani, 20 0 0;Kushch et al., 2005 ). Kachanov considered the average traction exerted on a crack by the other cracks in the system; this was later modified ( Gorbatikh et al., 2007;Li et al., 2003 ) https to account for linearly varying traction along the crack. Another method approximated the traction using a polynomial expansion ( Horii and Nemat-Nasser, 1985 ), which further improved the accuracy of the representation of the traction along the crack. This method was used and developed in both Maxwell's method and the Multipole method. All of these methods have their strengths, and weaknesses. Here, the MPM is chosen because it can accurately and semi-analytically capture the elastic interactions of multiple cracks and inclusions in one system; the MPM is used to simulate the many-body self-consistent problem of elliptical microcrack and inclusion interactions in heterogeneous solids ( Kushch, 2013a;2013b;1998a;1998b;Kushch and Sangani, 20 0 0;Kushch et al., 2005 ). It is a series method that reduces the boundary value problem in the homogeneous domain to an ordinary system of linear algebraic equations.
By characterising the localised stress fields and crack/inclusion interactions using the multipole method, we establish the DCD method. This method is based on the assumption that in brittle materials incremental crack growth can be represented by a sequence of micro-crack initiation events ahead of the original crack tip, which are subsequently merged with the original crack ( Carpinteri and Ingraffea, 2012;Hoagland et al., 1973;Li and Shah, 1994;Ortiz, 1988;Pompe et al., 1978;Wu et al., 1978 ). In this paper, the basic notation of the multipole method is introduced, and a general solution for a system comprising a finite array of elliptical inclusions and cracks is proposed. Once the distribution of stresses and strains is computed using the MPM, a crack propagation criterion is employed to determine the crack propagation path. The algorithm of crack propagation is elaborated in Section 3 and the results of numerical test problems are summarised in Section 4 in addition to a validation against an experimental result for a crack interacting with a circular inclusion.

Methodology
The multipole method provides a fast, accurate solution for multi-inclusion systems (where cracks may be thought of as collapsed elliptical inclusions). The method evaluates the local stress field in the vicinity of the inclusions using the linear superposition principle. The local field is the linear sum of the incident far-field and the fields generated by all inclusions. The problem has to be solved for the entire array of inclusions simultaneously. With the exception of the basis functions used in the expansion, which must be chosen according to the geometry of the inclusions, the solution procedure is identical for different shapes of inclusions ( Kushch, 2013a;2013b;1998a;1998b;Kushch and Sangani, 20 0 0;Kushch et al., 2005 ).
A general problem contains a distribution of shapes of cracks and inclusions, and the latter may deviate considerably from perfectly circular. Using elliptical inclusions instead of circular is of particular interest because both straight cracks and inclusions can be treated at the same time and with the same formalism, by regarding the former as collapsed, zero stiffness inclusions. We apply the multipole method here to a brittle material containing cracks, and elliptical inclusions with initially perfectly bonded ( i . e . displacement compatible) interfaces. Basis functions for elliptical inclusions were introduced by Kushch (2013a,b) .

Geometry
Consider a single elliptical inclusion, embedded at the origin of the Cartesian coordinate system ( x 1 Ox 2 ), in an infinite isotropic elastic solid. The Ox 1 and Ox 2 axes are along the major and minor axes of the ellipse, which have lengths 2 l 1 and 2 l 2 respectively. The aspect ratio and the distance between the foci of the ellipse are e = l 2 /l 1 < 1 and d = l 2 1 − l 2 2 respectively. Conformal mapping is a mathematical technique, applied here, that converts one problem into another that is easier to solve. The most convenient conformal mapping for problems containing elliptical objects is Kushch et al. (2005) Setting ζ equal to a constant, ζ 0 , the Cartesian axes become x 1 = d cosh ζ 0 cos η and x 2 = d sinh ζ 0 sin η, which describes an ellipse in the x 1 − x 2 plane with major and minor axes 2 d cosh ζ 0 and 2 d sinh ζ 0 respectively, for which the aspect ratio is e = tanh ζ 0 , and d is again the distance between the foci. Thus the locus of the ellipse in the z -plane corresponds to a straight line ζ = ζ 0 in the ξ -plane. Perfect bonding is assumed at the matrixinclusion interface ( ζ = ζ 0 ), thus displacement u and traction t are continuous there: Upper indices 0 and 1 refer to the matrix and inclusion materials, at the interface, respectively.
When there are N elliptical inclusions in the material, we define N parallel local coordinate systems ( x 1 p , x 2 p ), each with its origin at the centroid of the ellipse; another parallel coordinate system is designated as the global coordinate system ( x 1 , x 2 ). Defining z p = x 1 p + ix 2 p and z = x 1 + ix 2 , the coordinates of a point expressed in the global coordinate system ( x 1 , x 2 ) relates to the coordinates of the point expressed in the local system ( x 1 p , x 2 p ) by z = z p + Z p , where Z p is the coordinate of the centroid of the p th inclusion in the global coordinate system.
To describe the differences in the relative orientations of the ellipses another set of N local rotated coordinate systems ( y 1 p , y 2 p ) is required, aligned along the major and minor axes of each ellipse. If the coordinate axes ( y 1 p , y 2 p ) are related by an anticlockwise rotation of θ p relative to the local coordinate axes ( x 1 p , x 2 p ) then y p = y 1 p + iy 2 p = e −iθp z p . The conformal transformation y p = D p cosh ξ p , where 2 D p is the distance between the foci of the p th ellipse, introduces the local elliptic coordinates ξ p = ζ p + iη p . Then z p = D p e iθp cosh ξ p .

Formulation of the 2D elastic field problem in terms of complex potentials
We introduce the Muskhelishvili-Kolosov complex analytic potential functions ϕ( z ) and ˜ ψ (z) , in terms of which the stress tensor components σ ij and displacement vector u = (u 1 , u 2 ) may be expressed in plane strain as follows: where μ is the shear modulus, ν is the Poisson's ratio, and κ is Kolosov's constant factor equal to 3 − 4 ν for plane strain and 3 −ν 1+ ν for plane stress. The complex conjugate of a function f ( z ) is written as f (z) , and the primes denote differentiation. The displacement vector u in the complex form is u = u 1 + iu 2 .
(2) -( 4 ) become: 2.3.1. Arbitrary array of inclusions and cracks (pseudo-inclusions) For a 2D solid that contains a finite array of arbitrarily oriented elliptical inclusions, the displacement fields in the vicinity of any inclusion comprise disturbance fields generated by all inclusions in the system. Using the linear superposition principle, the displacement u ( z ) for the multiply connected matrix in the global coordinate system is: where u far and u dis are displacement fields caused by an external far-field load, and a total disturbance field induced by all inclusions respectively. Fields generated by inclusions are written in their local coordinate systems. However, it is necessary to write all the fields in a single coordinate system, hence in this section fields generated by the p th inclusion are expanded again in the q th coordinate system.
The displacement field of the p th inclusion in the local parallel coordinate system z p is u p (z p ) = u p (x 1 p , x 2 p ) . When this displacement field is expressed in the local rotated coordinate system y p it has the form u p (y p ) = u p (y 1 p , y 2 p ) . Here, a function is defined to be regular if it and its derivatives are bounded within the region, otherwise, it is singular (irregular). The displacement field induced by the p th inclusion in its coordinate system is singular as the p th inclusion is approached, but is regular in the vicinity of other inclusions. Upper indices s and r refer to singular and regular functions. The complex displacement caused by the p th inclusion in its coordinate system is denoted by u s p ; in the local q th coordinate system, this field is denoted by u r pq .
Eq. (9) becomes u s p (y p ) e iθp = u r pq (y q ) e iθq in the local y coordinate systems. The displacements u s p (y p ) and u r pq (y q ) can be expressed in terms of complex potentials using Eq. (7) : where ϕ r pq (y q ) = ∂ ϕ r pq /∂ y q and ϕ s p (y p ) = ∂ ϕ s p /∂ y p . The scalar potentials that solve the problem in the p and q coordinate systems are: By equating similar terms in u s p (z p ) and u r pq (z q ) and transforming the basis functions, the expansion coefficients a mpq and b mpq in the q coordinate system can be derived in terms of A np and B np . Details of the basis functions v q in different coordinate systems and the derivation of the coefficients a and b can be found in Appendix A and Appendix B . Superimposing the disturbance fields of all inclusions in the vicinity of the q th inclusion and adding it to the far-field gives the total disturbance for a finite array of N non-overlapping elliptical inclusions: Transforming all the disturbance fields into one coordinate system gives: Finally, displacement u ( z ) in the local z q coordinate system becomes: which means that the condition experienced by each inclusion is that of a single equivalent inclusion. In other words, the problem is reduced to the one inclusion problem with a modified farfield traction and a 4 N × n set of linear equations. By applying the boundary conditions of a perfectly bonded inclusion-matrix interface, the unknown coefficients A and B can be obtained. Once all coefficients including a npq and b npq are obtained, the scalar potentials and consequently the stress and displacement fields can be calculated ( Fig. 1 ). For simplicity, the problem of one inclusion is solved in the local y coordinate system. Therefore, the far-field for each inclusion has to be rotated according to its y local coordinate system. For instance, for a uniform far-field tension, the scalar potentials ϕ 0 far and ψ 0 far have to be expanded as follows: The non-zero coefficients of the field in the local y coordinate system are as below; details of the derivation can be found in B.3 : where d q = D q e iθq . After the problem is solved in the local y q coordinate system, a rotational transformation is applied to the calculated fields to transfer them into the local z q coordinate system using σ (z q ) = R σ (y q ) R T where R is the rotational transformation tensor. Additional translation is required to convert the stress fields from the local z q system into the global coordinate system.

Stress fields on the perimeter
To verify the reliability of the multipole method, stress fields of one, two and three mis-oriented elliptical holes are analysed. All simulations are performed under a uniform tensile load in the x 2 direction. In the first set of simulations, stress fields on the perimeter of one elliptical hole oriented at 45 degrees with respect to the global axes and subjected to remote tensile load ( σ ∞ ) in the x 2 direction are calculated using the commercial finite element package ABAQUS (FEM), and the multipole method. Stress fields are illustrated along the perimeter of the inclusion starting from the point (y 1 , y 2 ) = (0 , l 1 ) , following the inclusion in the anticlockwise direction in the local y coordinate system.
Another set of simulations was run to calculate the stress fields on the perimeter of two and three arbitrarily placed elliptical holes with aspect ratio 0.5 ( l 1 = 1 μm and l 2 = 0 . 5 μm ). Fig. 2 illustrates the stress component σ x 2 x 2 /σ ∞ along the perimeter of the q th inclusion starting from the point (y 1 q , y 2 q ) = (0 , l 1 ) , following the inclusion in the anticlockwise direction in the local y q coordinate system; in Fig  (shown in blue) starting from the point: (y 11 , y 21 ) = (x 11 , x 21 ) = (0 , l 1 ) . Fields calculated using the multipole method and FEM are represented by the blue and black lines, respectively.

Stress intensity factors
Analytical expressions for the stress intensity factors of the cracks and stress concentrations of the inclusions can be derived using the complex potential technique ( Kushch, 2013a ). Here, K + and K − refer to the right and left crack tips with respect to the local coordinate system. For inclusions the stress concentrations are calculated at the rounded edge intercepting the major axis of the inclusion. In linear elastic fracture mechanics the SIFs are defined as: where the term ( z q ∓d q ) defines the distance between the branch points and the surface of the inclusion. For the case where e q → 0, the SIFs become: After some algebra, and considering the fact that complex potentials are singular in the local coordinate system, one can obtain the formula for the case where ζ q → 0: This formulation enables calculating the SIFs without any additional numerical effort. The formula is valid for inclusions with nearly zero aspect ratio and arbitrary far load. The SIFs cannot be determined for an ellipse unless we assume that a crack is emanating from the ellipse. Analytical solutions for cracks emanating from elliptical inclusions are very difficult to obtain. Here the evaluation of SIFs for cracks emanating from elliptical inclusions is addressed by averaging over a short distance from where the crack emanation is predicted. Once the distribution of stress fields are calculated, the average SIFs are computed The influence of the length l is investigated and results of calculations are compared with an analytical method given by Panasyuk and Buina (1968) , Berezhnitskii (1967) , and Panasyuk and Berezhnitskii (1966) for an elastic plane weakened by an elliptical hole with a crack terminating at its edge. Fig. 3 shows a schematic representation of the problem; a crack emanating from an elliptical hole, in an unbounded solid subjected to remote Mode I loading. The host material is assumed to behave elastically and isotropically. SIFs are calculated for an ellipse with various aspect ratios. The SIFs are averaged over the lengths l/l 1 = 0 . 05 and 0.005. Results of simulations are presented in Fig. 4 . For comparison purposes, the figure also shows results obtained by Berezhnitskii. It is observed that the present numerical results are in excellent agreement with the analytical results, especially for the voids with larger aspect ratios.

Crack propagation: discrete crack dynamics (DCD)
Discrete crack dynamics, a new tool to study fracture in brittle materials, is developed based on four postulates: (i) each crack can be modelled as a collapsed zero-stiffness inclusion, (ii) crack growth is incremental and quasi-static, (iii) crack growth can be modelled by introducing a small crack (or pseudo-inclusion) ahead of the existing crack tip in the predetermined propagation direction at each step. Crack growth can be thought of comprising sequential micro-crack initiation events ahead of the main crack tip followed by merging of the crack ensemble ( Carpinteri and Ingraffea, 2012;Hoagland et al., 1973;Li and Shah, 1994;Ortiz, 1988;Pompe et al., 1978;Wu et al., 1978 ), (iv) two cracks (psudoinclusions) can be approximated by a single kinked crack if their sizes and the distance between them satisfy specific criteria, which will be studied in detail later.
The first postulate makes it possible to find the interactions between cracks and inclusions of the system simultaneously. The sec- ond and third postulates allow the introduction of a small crack at the existing tip to model incremental crack propagation. The direction of crack propagation can be determined using a suitable mixed mode crack propagation criterion. The fourth postulate is valid only if the SIF at the tip of the second crack can approximate the SIF of the two-crack ensemble, which might not be straight. A consequence of the fourth postulate is that the stress fields between two adjacent cracks becomes critically high and the field at the "junction" between the two cracks (see Fig. 5 ) affects the stress fields and the SIFs at crack tips far away from the junction. Thus, a modification has to be applied to derive the SIFs of kinked cracks from two adjacent cracks; details of the correction are provided in Section 3.1 .
Studying crack propagation in elastic materials requires solving three distinct problems. Firstly, a method to calculate interactions between cracks and inclusions and determine the stress and displacement fields around the crack must be available. Secondly, the method needs to be capable of computing the fracture parameters, such as the stress intensity factors or the energy release rate. In this study, the multipole method is used to address the first two problems. Thirdly, a criteria is needed to characterise the conditions leading to crack propagation, as well as the direction and

Crack propagation criteria
In most models of brittle fracture the quasi-static evolution of a crack is governed by the energy release rate along the crack path and it is closely related to the metastability of the crack ( Griffith, 1921;Negri and Ortner, 2008 ). Another class of mathematical model for fracture in brittle materials, based primarily on global stability, was proposed by Francfort and Marigo (1998) and has subsequently been developed by Larsen (2003) , Dal Maso et al. (2005) , and recently by Sutula et al. (2017c,b,a) for multiple crack propagation in a linear elastic solid under quasistatic conditions. The former formulation of fracture based on the energy release rate relies on local stability while the latter requires global stability, i.e. global minimisation of the free energy. Although the formulations of the two approaches are quite different, it was shown by Francfort and Marigo (1998) and Negri and Ortner (2008) that the resulting solutions should coincide in many cases. In this paper, we use the energy release rate concept and the local energy minimisation to model the propagation of individual microcracks.
A fracture criterion is required to accurately determine the crack propagation path. A number of fracture criteria has been developed for a crack under critical mixed mode loading ( Erdogan and Sih, 1963;Palaniswamy, 1971;Palaniswamy and Knauss, 1972;Hussain et al., 1973;Zhang and Eckert, 2005 ). Here the general energy release rate fracture criterion is used. An approximate expression of the energy release rate ( G ( θ )), as a function of angledependant SIFs, was presented by Chang et al. (2006) .

Crack propagation: length scale
Predictions of crack propagation based on the stress fields at the notch root (and/or the crack tip) are notoriously inaccurate. The calculated maximum stress on the surface is unsuitable for crack propagation prediction in situations where the gradient of the stress near the notch is high ( Nowell and Dini, 2003;Dini and Hills, 2004;Dini et al., 2006 ). Such high stress fields decrease significantly over short distances; it is not the entire stress field that is important in studying crack propagation, only the field over a critical finite distance ahead of the crack tip. This critical distance relates to the microstructure of the material and it can be determined using e.g. the Finite Fracture Method (FFM) ( Cornetti et al., 2006;Taylor et al., 2005 ).
FFM is based on the assumption that crack growth is not smooth and continuous, but instead occurs in a discontinuous manner. The size of the extension ( L ) is determined by the microstructure and deformation behaviour of the material. In FFM the failure condition is obtained using an energy balance similar to that of Griffith, but assuming a finite amount of crack extension. For the case of a sharp crack, the value of L = 1 π ( K C σ 0 ) 2 is expressed as a function of two other material constants, the fracture toughness K C and the inherent strength parameter σ 0 . This formula may require corrections for crack blunting. A similar material constant was proposed by Peterson (1959) , Neuber et al. (1946) , andEl Haddad et al. (1980) using the point method, line method, and imaginary crack method respectively.

Stress intensity factors of kinked cracks
It is important at this point to confirm the validity of the last postulate, namely that two adjacent cracks can be approximated by a single kinked crack. In general, two cracks cannot be approximated as one larger crack ( Rubinstein, 1985 ) unless modifications are applied to the system. In this section, a method to derive the SIFs of a kinked crack comprising two adjacent cracks is discussed.
Consider a pair of cracks, p and q where crack p is short compared to crack q , and their separation δ is small compared to the length of either crack ( Fig. 5 a). Under the external loading, crack   6. The ratio of the optimised distance between two inclusions, δ, to the second crack/inclusion length, L 2 , is plotted versus the ratio of the second to the first crack/inclusion length, L 2 / L 1 . q produces disturbance fields in the vicinity of Z p and visa versa. Kushch et al. (2005) and Sevostianov (2015, 2016) derived a formula for the average stress in the p th inclusion < σ > p , caused by the q th inclusion in terms of expansion coefficients a n , q , p and b n , q , p . Berezhnitskii et al. (1973) showed in a system containing two cracks that under the critical external load, cracks grow initially one towards the other then merge, and consequently a single kinked crack is formed. In DCD, cracks are propagated by inserting a microcrack a small distance ahead of the crack tip, in the direction of crack propagation, as depicted in Fig. 5 a. Therefore, it is necessary to establish that this ensemble can, to a very good approximation, represent the equivalent contiguous kinked crack (in terms of the SIFs and energy release rate) shown in Fig. 5 b, by suitable choice of the microcrack length and separation distance. Hence SIFs and energy release rate are computed for the p th crack at tip B (see Fig. 5 a) also incorporating the disturbance field at tip B caused by the longer q th crack. To simplify the calculation, the average disturbance fields introduced in Kushch et al. (2005) are used.
The average disturbance fields in one crack in the binary crack system caused by the other depends on the distance δ ( Fig. 5 a).
If cracks are positioned very close to each other, the large stress fields of the adjacent tips will reduce the accuracy of the approximation. Also, if they are positioned too far from each other, the ensemble in Fig. 5 a is incapable of representing the equivalent kinked crack in 5 b. Therefore, there is an optimum δ that makes the approximate SIFs and energy release rate of Fig. 5 a maximally accurate. The case of two aligned cracks in Fig. 5 c is more sensitive to δ than any other configuration ( Mehdi-Soozani et al., 1987;Peng and Sung, 2003;Gdoutos, 2012 ), hence δ is optimised for that case and the value is used for all configurations.
Since the optimisation of δ also depends upon the crack lengths, a number of simulations were run using different values of the crack lengths and junction size δ, to find the optimum δ. The optimum δ is the distance for which the two crack ensemble of Fig. 5 a (of length L 1 + L 2 + δ) differs by no more than 1% in its SIFSs at tip B compared to those of the equivalent contiguous crack for the aligned cracks case (see Fig. 5 c). Fig. 6 shows the optimum value of δ for different crack lengths L 1 and L 2 .  Table 1 present results for a successful validation test performed using kinked crack systems. Normalised SIFs of the kinked cracks calculated by analytical methods using the complex potential formulation and Greens function technique ( Lo, 1978;Chen, 1999 ), and the discrete dislocation density method Hills and Nowell (1990) are compared with results obtained using the multipole method for two adjacent cracks separated by the optimum distance δ; the value of δ was obtained from the results of Fig. 6 for the corresponding b / a ratios used in the comparison, not-ing that a = L 1 and b = L 2 + δ. Analysing the results shows that normalised SIFs calculated using the MPM never differ by more than 5% for the φ = 15 • , 30 °, 45 °and 60 °degree cases, compared to the results of Lo (1978) , Hills and Nowell (1990) and Chen (1999) , respectively. Accuracy of the MPM depends on the number of harmonic terms, n , used in the calculations. According to the study performed by Kushch (2013a) , n = 25 gives accurate results, which was also confirmed in this study by comparing predictions to FEM results; hence that value was used here. Data in Table 1 shows that discrepancies are larger for larger angles φ.

Crack propagation: simulation steps
In the DCD method, crack propagation is a discretised process (based on the concepts of FFM Cornetti et al. (2006) ; Taylor et al. (2005) ), in which the crack advances incrementally at each simulation step. The general algorithm to study the evolution of a crack is as follows. If the crack propagation condition is satisfied, a micro-crack is added ahead of the main crack-tip at a distance equal to the optimal spacing determined in 3.1 . The length of the micro-crack (the crack extension) is determined from the theory of FFM; the angle of orientation of the micro-crack is prescribed according to the adopted fracture criterion. The arrangement of the original crack and the new crack is an approximation of a kinked crack, using the theory developed in Section 3.1 . Accordingly, the SIFs and the energy of the approximated kinked crack are calculated.
A problem arises in the next step of the simulation when the kinked crack propagates further. The discretised process of crack propagation requires introducing another micro-crack at the tip of the crack/micro-crack ensemble. At this point, calculating the SIFs of the propagated kinked crack comprising one main crack and a sequence of two micro-cracks, as shown in Fig. 8 a, is not straightforward.
In fact, the crack propagation path cannot be determined simply by adding a succession of small micro-cracks to the main crack-tip, because of the shielding effect associated with the micro-cracks in this configuration. Crack-tip shielding caused by micro-cracking was studied by Chudnovsky et al. (1984) , Hutchinson (1987) and Shum and Hutchinson (1990) . An analytical study was done by Rubinstein (1985) for an array of aligned micro-cracks ahead of a main crack-tip. It was proved that by increasing the number of micro-cracks used to represent the crack extension, the effect of the main crack on the local SIFs and stress field distributions of the leading micro-crack progressively diminishes. Hence the estimated crack propagation direction and energy release rate based on the leading micro-crack becomes unreliable. Precise prediction of the SIFs and energy release rate are necessary to determine the crack propagation path. Therefore, another method to determine the crack propagation path is required. Here a solution is given for a system of three cracks, which is then generalised for longer crack propagation in the DCD method.
The problem of a propagating kinked crack may be addressed by approximating a set of three cracks, comprising a main crack and two micro-cracks, using just two cracks: either by replacing the main crack and the first micro-crack with an equivalent crack ( Fig. 8 b), or by merging the two micro-cracks ( Fig. 8 c) representing the crack extension. Before using either of these approximations, it is necessary to assess their accuracy and limitations. The point of this approach is to always represent a propagating kinked crack using only a main crack and a single micro-crack ahead of it, since it was proved in 3.1 that the SIFs and energy release rate for this configuration are an accurate representation of the actual propagating kinked crack provided the optimal spacing between them is used, and the micro-crack length is based upon FFM theory.
(a) (b) Fig. 7. a) Kinked crack benchmark problem. Appropriate δ is selected according to Fig. 6 ; b) Normalised SIF of a kinked crack with b/a = 1 using the MPM method.

Table 1
Variation of the normalised SIF ( K/σ √ πa ) with angle for a kinked crack of the same size as the original crack under a uniform tensile stress. Suffix B refers to the endpoint (tip B) of the kinked crack (see Fig. 7 ). Methods 1, 2, 3 and 4 are the MPM method, the complex potential formulation, discrete dislocation density and Greens function technique, respectively ( Lo, 1978;Chen, 1999;Hills and Nowell, 1990 ). It was shown in Bilby and Cardew (1975) , Cotterell and Rice (1980) , and Suresh (1983) that the SIFs of a kinked crack depend on the inclination angle of the crack extension and the SIFs of the main crack tip prior to its extension. Kitagawa et al. (1975) and Bechtle et al. (2010) showed that, from the standpoint of the SIFs, a kinked crack can be approximated by a straight line crack in which the orientation is the same as that of the crack extension; and Chen (1999) showed that the angle between the crack-tip and load directions determines the proportion of Mode I and II, while the distance between the two extents of the crack ensemble is a dominant factor on the SIFs.
Considering the finding of these studies and the fact that the proportion of Mode I and II at the lead crack tip determines the orientation of the next increment of crack extension ( He and Hutchinson, 1989;Hutchinson, 1990 ), it is clear that the orientation of the last micro-crack in the problem of three cracks shown in Fig. 8 a must remain the same in any approximate representation. If the main crack and the first micro-crack are replaced with an equivalent crack as shown in Fig. 8 b, the length of the equivalent crack is roughly the same as the distance between the extents of the crack ensemble, thus, according to Chen (1999) the SIFs in Fig. 8 b should be a good approximation of those at the actual propagated crack-tip. Fig. 8 d shows the orientations of the last two micro-cracks in their local coordinate systems. It might be more accurate to merge the two micro-cracks ( Fig. 8 c) rather than replace the main crack and first micro-crack with an equivalent crack. The criterion for whether to merge the two consecutive micro-cracks is a function of the relative angle between them ( θ = θ 2 − θ 1 ): if the mismatch in orientation ( θ ) of the two micro-cracks is less than a threshold amount, then the cracks will be merged. The threshold angular mismatch is that which produces a better estimate of the SIFs at the lead crack tip by merging the two micro-cracks, than what is obtained from the alternative combination of equivalent crack and single micro-crack; it depends on the material and geometry of the problem. The same methodology can be used in all of the iterative steps in predicting the crack trajectory. The above methodology is adequate to deal with individual micro-cracks propagating within a loaded system in the presence of macro-cracks and inclusions. This can be extended to deal with a network of micro-cracks for example by performing minimisation of the total energy of the mechanical system with respect to the crack extension directions and crack extension lengths to solve for the evolution of the mechanical system over time, as proposed in Sutula et al. (2017c,b,a) .

Finite-size problems
The aim of this section is to extend DCD to finite-sized problems using a linear superposition scheme to incorporate boundary effects; the approach is the same as that commonly used in discrete dislocation dynamics modelling. The superposition principle is exact at each simulation step and it requires no modification to the multipole method described here. Let S be the domain of the boundary-value problem with boundary Γ delimiting a solid body characterised by stress and displacement fields σ( x , y ) and u ( x , y ).
By virtue of the linear superposition principle, these can be written as the sum of two fields: σ = ˜ σ + ˆ σ and u = ˜ u + ˆ u where ˜ σ and ˜ u represent fields of the infinite domain containing the inclusions, and ˆ σ and ˆ u represent the fields of the finite-size inclusion-free medium ˆ S coinciding with S , in response to the corrective boundary conditions. As a result of inclusion fields, surface ˜ Γ , which coincides with Γ , experiences a traction ˜ t and displacement ˜ u . In order for the superposition of the fields in ˜ S and ˆ S to equate to the fields in S , the surface ˆ Γ must have, for every simulation step, -˜ t and -˜ u applied over it in addition to the boundary conditions applied on Γ .
To summarise, the problem under consideration is a 2D linear elastic body of area S comprising elastic inclusions and holes ( Fig. 9 ) subjected to the remote far-field displacement and traction boundary conditions. The elastic properties of the matrix material are given by the shear modulus μ 0 and Poisson's ratio ν 0 , and the elastic properties of the i th inclusion are given by the shear modulus μ i and Poisson's ratio ν i . The stress and displacement fields of the system are determined by decomposing the problem of the finite body with inclusions into two problems: the problem of interacting inclusions in the homogeneous infinite solid and the complementary problem for the homogeneous body without inclusions subject to the corrective boundary conditions. The former is solved using the multipole method (MPM) formulation, and the latter can be solved using a numerical scheme such as the finite element or boundary element methods. The complex variables boundary element method (CVBEM) Aliabadi and Rooke (1991) , Linkov (1998) , andZografos (2011) uses the same mathematical framework as the MPM and is therefore chosen to solve the finite boundary value problem in the absence of inclusions. In this study, the CVBEM code developed by Zografos (2011) and Zografos and Dini (2009) was integrated with the DCD method to determine the effect of finite geometry on crack propagation.

Numerical results
Simulations in this section were conducted to establish the validity of the DCD method, which could become a standard micromechanics model for dealing with crack propagation in networks of inclusions or populations of other objects that can be represented as inclusions or collapsed inclusions.
To start, the DCD method was used to predict the propagation of an inclined crack under uniaxial tensile loading in the x 2 direction. Fig. 10 shows the predicted propagation path of an inclined crack for various inclination angles. Fig. 10 shows that the crack propagation path tends toward the direction of maximum Mode I loading. This agrees with previous calculations ( Patrıcio and Mattheij, 2007;Cherepanov, 1979;Meggiolaro et al., 2005 ), as well as accepted theoretical and experimental understanding of crack propagation in brittle materials.
The same problem was solved under uniform compression in the x 2 direction. Fig. 11 graphically compares the numerical results obtained using DCD and Boundary Element methods given by Haeri et al. (2014) , and demonstrates the accuracy and validity of the DCD method. Results are in agreement with other predictions given by Horii and Nemat-Nasser (1986) and Nemat-Nasser and Horii (1982) .
Sets of simulations were run to establish the effect of impurity and micro-crack interactions on the deflection of the main crackpath. The crack propagation trajectory in the vicinity of circular soft and hard inclusions was studied on a plate with side length 10 μm. The plate was loaded remotely in tension in the x 2 direction. The initial crack length a and the inclusion radius r were set to 1 μm, and the inclusion was centred at (5, 1.5) μm. Fig. 12 a shows the accuracy of the predicted crack propagation path in an infinite system for a hard inclusion case with inclusion stiffness 10 times that of the matrix. The green ( n = 1 ), blue  ( n = 3 ) and red ( n = 9 ) curves represent crack paths for a different number of harmonic terms in the expansion. The crack paths are indistinguishable for n = 9 and higher values of n ; this value of n was used in subsequent calculations. Fig. 12 b and c show crack path deflection in the presence of circular hard and soft inclusions in a finite plate. The plate was fixed at the bottom edge with load distributed uniformly on the top edge. The relative stiffness of the inclusion to the matrix was taken to be 2, 10, 100, and 1/2, 1/10 for the hard and soft inclusion cases, respectively. Fig. 12 d shows that DCD simulations are comparable to modelling results in Nielsen et al. (2012) . The results also agree with other simulations and experimental predictions given by Bouchard (2005) , Patton (1991) , Patton and Santare (1993) , and Misseroni et al. (2015) .
Finally, crack paths predicted by the DCD method are compared to the asymptotic model and experimental result presented in Misseroni et al. (2015) . A 120 × 95 mm plate containing elliptical holes was subjected to Mode-I loading; elliptical holes with major and minor axes of a = 4 and b = 0 . 6 mm were placed at (x 101 , x 201 ) = (75 , −9) and (x 102 , x 202 ) = (105 , 8) mm locations. The inclinations of the elliptical voids are θ 1 = 0 and θ 2 = π / 2 , respectively. The experimental set up shown in Fig. 13 was reproduced in the simulations. In Fig. 13 , the blue and black  curves show the crack path predicted by the asymptotic model ( Misseroni et al., 2015 ) and the DCD method, respectively, whereas the red curve shows the actual crack path obtained in the experiment. Excellent agreement between the analytical results of the asymptotic solution, the DCD method and the experiments is observed, further validating the DCD method.

Conclusion and further work
In this paper, the so-called Discrete Crack Dynamics (DCD) method was developed, based upon the multipole method (MPM), to simulate the interactions between a finite array of arbitrarily oriented cracks and inclusions with varying aspect ratios and material properties in a finite geometry. The method provides an accurate semi-analytical tool to derive full stress and displacement fields, as well as stress intensity factors and energy release rates of cracks, and to predict low-speed crack propagation through complicated defect networks in brittle or quasi-brittle materials in a very efficient way that is much faster than computational methods such as the extended finite element method. A set of test problems with known analytical, computational and experimental results was used to assess the method, and it was found to be very accurate for these cases. Forthcoming studies will examine in detail evolving low-speed crack networks in brittle or quasi-brittle materials, cracks developing very near inclusion interfaces and grain boundaries, and predict the paths of cracks propagating through grain networks with a background inclusion distribution.
The second term can be written as:  On the other hand, complex potential ψ r pq (y q ) is defined in Eq. (12) as : Using Lemma B.18 , ψ r pq (y q ) becomes:  b −npq = b npq − 2 a npq | n | sinh 2 ξ 0 q (B.34)

B.3. Expansion coefficient of uniform far-field load
Complex potentials for uniform far-field loading have to be selected as a linear function of z due to constraints imposed by the displacement and traction in the complex potential form. Potentials are chosen as: Here, u 0 q , and u 00 q are displacement fields caused by the farfield loading in the global ( z ) coordinate, and in the local y q coordinate systems respectively. By substituting ϕ 0 , and ψ 0 into u 0 q and using formulas z = Z q + z q , and z q = y q e iθq , u 0 q is derived in terms of u 00 q . u 0 q (z) = κϕ 0 (z) − (z −z ) ϕ 0 (z) − ψ 0 (z) = κ 1 z − (z −z ) 1 − 2 z = κ 1 Z q − (Z q − Z q ) 1 − 2 Z q + κ 1 z q − (z q − z q ) 1 − 2 z q = u 00 q (Z q ) + u 01 q (z q ) (B.36) As far-field load is a regular field at infinity, complex potentials should obey where 0 denotes the far-field load. Using formula a mpq = a −mpq , and b mpq = b −mpq + 2 ma mpq sinh 2 ζ 0 q , complex potentials are derived as: ϕ r 0 q (y q ) = a 00 q + a 10 q (v −1 q + v q ) ψ r 0 q (y q ) = b 00 q + 2 z q d q b −10 q + a 10 q 1 − 1 v 2 0 q (B.38) Using the u 0 q (z) = e iθq u q (y q ) relation between displacements in coordinate systems z , and y q , complex potentials in the y q coordinate are predicted. Displacements are calculated by: Displacement u q contains two parts, a constant term and a linear function of coordinate. u q (y q ) = κϕ 00 − (y q −ȳ q ) ϕ 00 − ψ 00 + κϕ 01 (y q ) − (y q −ȳ q ) ϕ 01 (y q ) − ψ 01 (y q )