A model based on discrete unit events cou- pled with a graph search algorithm is developed to pre- dict intergranular fracture. The model is based on two hypotheses: (i) the key unit event associated with inter- granular crack propagation is the interaction of a grain boundary crack with a grain boundary segment located at an angle with the initial crack plane; and (ii) for a given crack path, the overall crack growth resistance can be calculated using the crack growth resistance of a collection of unit events. Next, using a directed graph containing the connectivity of grain boundary junctions and the distances between them, and crack de ection versus crack growth resistance data, a directed graph in theJ-resistance space is created. This graph contains information on the crack growth resistance for all possi- ble crack paths in a given grain microstructure. Various crack growth resistance curves are then calculated in- cluding those corresponding to: (i) a local resistance minimum; (ii) a global minimum; and (iii) for veri ca- tion, a path speci ed by microstructure-based nite el- ement calculations. The results show that the proposed method based on discrete unit events and graph search can predict the crack path and the crack growth resis- tance for cracks that propagate from one grain bound- ary junction to another. The proposed computationally inexpensive model can be used to design material mi- S. Osovski Faculty of Mechanical Engineering, Technion - Israel Insti- tute of Technology, Haifa, Israel A. Needleman Department of Materials Science and Engineering, Texas A&M University, College Station, TX, USA A. Srivastava Department of Materials Science and Engineering, Texas A&M University, College Station, TX, USA E-mail: ankit.sri@tamu.educrostructures with improved intergranular fracture re- sistance, and/or to assess the overall crack growth re- sistance of materials with a known distribution of grain morphology. Keywords Grain Boundaries Crack PathFracture ToughnessGraph SearchMicrostructure Design 1 Introduction The increasingly demanding requirements for lower fuel consumption and for lower CO 2emissions drives design toward lightweight structures in weight-critical applica- tions such as for aircraft and for long-haul trucks. High levels of reliability are also required in these applica- tions. Hence, there is an increasing desire to use struc- tural materials with high speci c strength (strength to weight ratio). Examples of materials with high spe- ci c strength include Al-Li alloys and high strength metastable titanium alloys. Unfortunately, these ma- terials are prone to grain boundary ductile fracture that limits their use in fracture-critical components or struc- tures. Here, we propose a novel methodology for de- signing material microstructures with improved crack growth resistance for materials that are prone to grain boundary ductile fracture. The proposed methodology can also be used to assess the overall crack growth re- sistance of materials with a known distribution of key microstructural features. The past few decades have seen rapid advances in inverse design methodologies that facilitate material de- sign to achieve a targeted performance metric (Olson, 1997; McDowell, 2007; McDowell and Olson, 2009; Full- wood et al., 2010; Liu et al., 2015). However, relatively few studies have been aimed at designing high tough- ness ductile materials (Hao et al., 2003, 2004; Kulka- 2 Shmuel Osovski et al. rni et al., 2004). In Hao et al. (2003, 2004) the focus was on fracture initiation in steels and the e ects of microstructural heterogeneities, such as void nucleat- ing particle distributions and grain morphology. Also, in these works the material microstructure entered the design space in a homogenized way, i.e. the material mi- crostructure was smeared out and/or de ned in terms of correlation functions. In Srivastava et al. (2017) it was shown that the duc- tile crack growth resistance of a material is strongly af- fected by the crack path and to the extent it is possible to engineer the crack path, the crack growth resistance can be increased by designing the material microstruc- ture. Thus, material design involving crack path engi- neering to increase the crack growth resistance involves optimizing the topological features of the material's microstructure. Various topology optimization schemes have been applied to continuum and discrete structures with the aim of designing optimal structural layouts for a variety of design objectives, e.g. Bends e and Kikuchi (1988); Kirsch (1990); Cheng and Jiang (1992); Chap- man et al. (1994); Bends e and Sigmund (2004); van Dijk et al. (2013); James and Waisman (2014); Xia and Breitkopf (2015). Topology optimization techniques have been used to design composite microstructures with desirable prop- erties, e.g. Sigmund and Torquato (1997); Larsen et al. (1997); Sigmund et al. (1998); Hyun and Torquato (2001); Torquato and Hyun (2001); Torquato (2015). Topological optimization of composite materials has usu- ally been focused on optimizing properties that can be found by homogenizing the local elds in the mi- crostructure. This, to an extent, allows using approxi- mate (but ecient) nite element models of the com- posite microstructure coupled with an iterative opti- mization scheme within a single computational frame- work. However, full eld microstructure-based nite el- ement calculations of crack growth in a speci ed mi- crostructure as in Osovski et al. (2015a); Srivastava et al. (2017) are very computationally expensive so that direct coupling of such calculations with an iterative optimization scheme is not practical. The microstructure-based nite element modeling in Osovski et al. (2015a) focused on understanding the e ects of material properties and grain microstructural features on the crack growth resistance of metastable titanium alloys. The microstructure of metastable ti- tanium alloys consists of hard grains of phase embed- ded within phase laths, and with a soft phase layer along the grain boundaries (L utjering and Williams, 2007). This microstructure leads to grain boundary duc- tile fracture at room temperature (Williams et al., 1978; L utjering and Williams, 2007; Mantri et al., 2015). Theresults of the microstructure-based nite element mod- eling (Osovski et al., 2015a) show that changes in grain morphology a ect both the crack growth initiation tough- ness and the tearing modulus (slope of the crack growth resistance curve). However, all the calculated crack growth resistance curves exhibited one common feature, namely a piecewise linear relation between crack growth re- sistance and crack extension that is occasionally per- turbed by jumps in crack growth resistance. The slope of the crack growth resistance curve was correlated with the morphology of the grain boundary network, whereas the jumps in crack growth resistance were found to be due to the interaction of the propagating crack with grain boundary junctions. The analyses in Osovski et al. (2015a) suggest that the crack growth resistance of a crack propagating through a grain boundary network can be described by discrete unit events. Based on these observations we formulate two hypothesis: (i) the key unit event associated with intergranular crack propagation is the interaction of a grain boundary crack with a grain boundary segment located at an angle with the initial crack plane; and (ii) for a given crack path, the crack growth resistance curve can be calculated using the crack growth resis- tance of a collection of unit events given the position, orientation and length of each grain boundary segment along the crack path is known. Here, the hypotheses formulated based on the anal- yses of Osovski et al. (2015a) are used to develop a sim- ple model based on discrete unit events coupled with a graph search algorithm to predict intergranular frac- ture in various grain microstructures. The relationship between the angle of crack de ection (due to the inter- action of a grain boundary crack with a grain boundary junction) and the associated crack growth resistance for a unit length of crack is carried out via a series of unit event microstructure-based nite element calculations. Next, the network of grain boundaries is constructed via a centroidal Voronoi tessellation scheme and is mapped onto a graph, were each grain boundary junction is a node on the graph. The graph is constructed in Eu- clidean space and all trajectories are extracted using a breadth- rst search algorithm (Lee, 1961). To calcu- late the crack growth resistance curve of a given crack path along the grain boundary network, a similar graph is built in the Jresistance space using the crack de- ection versus crack growth resistance data. Once con- structed, the graph contains information as to all pos- sible crack growth paths in a given microstructure and their corresponding crack growth resistance curves. Fi- nally, the path of least resistance is found by utilizing Dijkstra's algorithm (Dijkstra, 1959). Dijkstra's algo- rithm is a graph search algorithm that produces the Intergranular Fracture Prediction and Microstructure Design 3 Fig. 1 Sketch of the initially cracked slice of material ana- lyzed, showing the grain microstructure in the region in front of the initial crack tip. minimal distance between a given node and any other node on the graph. The distances in this case are the increments in crack growth resistance resulting from a crack traveling between two nodes. Calculations are carried out for grain microstruc- tures representative of metastable titanium alloys. The predictions are found to be in close agreement with the predictions of the full eld microstructure-based - nite element crack growth calculations. The comparison of the predictions of our simple model and the nite el- ement calculations also shows that the nal crack path in a speci ed microstructure is not necessarily the path of least resistance. Because of its computational e- ciency, the model proposed here has the potential for being employed in an iterative optimization scheme to design material microstructures with improved inter- granular fracture resistance, and/or to assess the over- all crack growth resistance of materials with a known distribution of grain boundaries. 2 Microstructure-based Material Model and Numerical Implementation The microstructure-based numerical modeling is car- ried out to quantify the e ect of grain morphology on crack growth resistance, as well as to calibrate and validate the simpli ed model formulated based on the graph search method presented in Section 3. The ma- terial model and the numerical implementation are thesame as in Needleman et al. (2012); Srivastava et al. (2014); Osovski et al. (2015a,b); Srivastava et al. (2017). Here, for completeness, we brie y describe the formu- lation. More complete descriptions are given in the ref- erences cited. As in Needleman et al. (2012); Srivastava et al. (2014); Osovski et al. (2015a,b); Srivastava et al. (2017), we analyze a mode I small scale yielding boundary value problem for a slice of material with an initial crack, Fig. 1. The region analyzed has dimensions hxhyhz with boundary conditions on z0= 0;hzgiving over- all plane strain. The initial and boundary displacement and velocity elds corresponding to the isotropic elastic mode I plane strain singular eld are applied. For small scale yielding the value of Jis related to the imposed stress intensity factor KIby Rice (1968) J=K2 I 12 E(1) whereis Poisson's ratio and Eis Young's modulus. A monotonically increasing remote stress intensity factor, KI(t) is prescribed with _KI= 1:2107MPapms1. A nite deformation, nite element formulation based on the dynamic principle of virtual work is used. The nite element formulation is three dimensional, using twenty-node brick elements and the internal force con- tributions are integrated using eight point Gaussian in- tegration. The explicit Newmark -method, with = 0, (Belytschko et al., 1976) is used for time integra- tion. Here, the microstructure is taken to be constant through the slice thickness with only one element through the thickness hzso that the deformation eld is essen- tially two dimensional plane strain. A uniform in-plane (x0y0plane) mesh is used in a region in front of the initial crack tip in which the material microstruc- ture is modeled as shown schematically in Fig. 1. The material modeled consists of a soft phase layer sur- rounding hard grains. The hard grains, referred to as grains, mainly consist of Ti but have embedded phase laths. Here, the grains are modeled as homoge- neous and isotropic. The grain microstructure in the region in front of the initial crack tip is generated using a centroidal Voronoi tessellation scheme as described in Osovski et al. (2015a). First,Ngrandom points are assigned in the region where the microstructure is modeled, Fig. 1. Each ran- dom point acts as a generator for a Voronoi cell (grain). For each cell, the distance between its generator and its center of mass is calculated and a new generator is created to minimize that distance. The asymptotic limit of this iterative scheme in 2D results in a perfectly equiaxed microstructure. To generate the grain bound- aries each Voronoi cell is shrunk so that a distance equal 4 Shmuel Osovski et al. to half the desired layer width (grain boundary) is generated between the original and new cell faces. For all the calculations reported here the width of grain boundary layer varies from 6 m to 10m and the av- erage grain size is 330m. The obtained list of poly- hedrons is then projected onto the ne meshed region in front of the initial crack tip and each Gauss integration point in the nite element mesh is assigned the mate- rial/constitutive properties associated with the phase in which it lies. This procedure gives a microstructure withNggrains. The material model is a modi ed Gurson elastic- viscoplastic constitutive relation for a progressively cav- itating solid as used in Osovski et al. (2015a). The ow potential given by Gurson (1975) =2 e 2+ 2q1fcosh3q2h 2 1(q1f)2= 0 (2) whereq1= 1:25,q2= 1:0 are parameters introduced in Tvergaard (1981, 1982a), fis the void volume fraction, is the matrix ow strength and e2=3 20:0; h=1 3:I;0=hI (3) withfgiven by Eq. (4) as introduced in Tvergaard and Needleman (1984) f=f; f <f c fc+ (1=q1fc)(ffc)=(fffc); ffc (4) where the values fc= 0:12 andff= 0:25 are used. The rate of deformation tensor is written as the sum of an elastic part, de=L1:^, a viscoplastic part, dp, and a thermal straining contribution, d= _I d=L1:^+ _I+dp(5) Here, ^is the Jaumann rate of Cauchy stress, is the temperature, = 1105/K is the thermal expan- sion coecient and Lis the tensor of isotropic elastic moduli. The plastic part of the strain rate, dp, as given in Pan et al. (1983) is dp=" (1f)_ :@ @# @ @(6) The matrix plastic strain rate, _, is given by _= _0 g(;)1=m ; g(;) =0G() [1 + =0]N(7) with =R_dtand0=0=E. The elastic constants E= 116GPa and = 0:3 were taken to be the samefor both phases as well as the hardening and rate sen- sitivity parameters given by N= 0:1,m= 0:01 and _0= 103s1. The ow strength of the phase grains,  0is taken to be 1200MPa and for the phase layer  0= 800MPa. Adiabatic conditions are assumed so that cp@ @t= :dp(8) with= 4500 kg =m3= 4:5103MPa=(m=sec)2, cp= 520 J/(kgK), = 0:9, and the temperature- dependence of the ow strength is given by G() = 1+bGexp(c[0273]) [exp(c[0])1] (9) withbG= 0:1406 andc= 0:00793=K and0= 293K. The initial void volume fraction is taken to be zero and the evolution of the void volume fraction is gov- erned by _f= (1f)dp:I+_fnucl (10) where the rst term on the right hand side of Eq. (10) accounts for void growth and the second term for void nucleation. Since fracture in metastable Ti alloys mainly oc- curs along grain boundaries (in the phase layer) the phase layer is taken to nucleate voids at early stages of the deformation process following the stress controlled nucleation criterion _fstress nucl =A(_+ _h); A=fstress N sstress Np 2exp 1 2 +hN sstress N2(11) if (+h) is at its maximum over the deformation history. Otherwise A= 0. Here, fstress N = 0:06,N= 1:50andsN= 0:3. Plastic strain controlled void nucleation is taken to occur both in the grains and in the grain boundary phase layer via _fstrain nucl =D_; D=fstrain N sstrain Np 2exp 1 2 N sstrain N2(12) withfstrain N = 0:04,N= 0:2 andsN= 0:2 for both phases. The constitutive updating is based on the rate tan- gent modulus method in Peirce et al. (1984) and failure is implemented using the element vanish technique of Tvergaard (1982b). When f= 0:9ff, the value of f is kept xed so that the material deforms with a very small ow strength. The entire element is taken to have failed (zero stress carrying capacity) once three Gauss Intergranular Fracture Prediction and Microstructure Design 5 integration points in the element have met this condi- tion. The extent of crack growth, a, beyond the initial crack is de ned as the maximum projected length on thex-axis of the void volume fraction contour f= 0:1 as in Srivastava et al. (2014); Osovski et al. (2015a,b); Srivastava et al. (2017). 3 Modeling Intergranular Cracking as a Graph Search The motivation for the graph search model comes from the numerical simulations in Osovski et al. (2015a). Those simulations of ductile crack growth in metastable titanium alloys gave results for the e ect of changes in grain morphology on: (i) the crack growth initiation toughness, JIc; and (ii) the tearing modulus (the slope of theJRcurve),TR= E=2 0 dJ=d (a), (Paris et al., 1979). All the crack growth resistance curves in Osovski et al. (2015a) exhibited a piecewise linear rela- tion between Jandathat are occasionally perturbed by jumps in the value of J. The variation in the tear- ing modulus TRwas correlated with the morphology of the grain boundaries that acted as a network through which the crack propagates. The jumps in the value of Jwere found to be due to the interaction of the propa- gating crack with the grain boundary junctions. When a propagating crack reached a grain boundary junction, it remained there until sucient weakening occurred on one of the grain boundary facets comprising the (triple) junction. This suggested that the crack growth resis- tance for a crack propagating through a grain boundary network could be described as a collection discrete unit events. Based on the results in Osovski et al. (2015a), our rst hypothesis is that the unit event to be modeled is the interaction of a grain boundary crack with grain boundary layers located symmetrically at an angle about the crack plane (the inset of Fig. 2a shows this con guration for the speci c case = 45). To obtain a relationship between the angle of crack growth, , and the associated change in the value of Jrequired for a unit length of propagation, we carried out a se- ries of nite element calculations using the formulation presented in Section 2. Calculations were carried out for the con guration shown in the inset of Fig. 2a with values of ranging from 0to 85. The variation of tearing modulus, TR, and increment in crack growth resistance, J, with the crack de ection angle, , are shown in Fig. 2a. The results shown in Fig. 2a provide the unit event data for our proposed model to predict intergranular cracking. Our second hypothesis is that knowing the position, orientation and length of each grain boundary seg-ment for a given crack path, the crack growth resis- tance curve can be calculated using the unit event data in Fig. 2a. Thus, for a given grain morphology, the crack growth resistance curve, Jversusa, can be calculated for any crack path emanating from the initial crack tip. To this end, the network of grain boundary is con- structed the previously described Voronoi tessellation scheme. The obtained network of grain boundary is then mapped onto a graph, were each grain boundary junction is represented by a node on the graph. Follow- ing this, the connectivity map is constructed such that ow on the graph is restricted to the direction of crack propagation (i.e. an increasing value of the crack length projected on the xaxis). The graph is constructed in Euclidean space so that all trajectories which traverse the graph up to a given distance can be extracted using a breadth- rst search algorithm (Lee, 1961). Next, to calculate the crack growth resistance curve for a given crack path along the grain boundary network, the ener- getic cost of crack propagating between the two nodes on the graph has to be known, i.e. a similar graph has to be built in Jresistance space. Using the unit event data in Fig. 2a in conjunction with the directed graph containing the connectivity of grain boundary junctions and the distances between them, a directed graph in the Jresistance space is created. The distances between the nodes in the Euclidean space and the propagation angle between them are input into the unit event model and the edge value of the new graph (i.e. the change in the value of Jrequired to transverse the given edge) is obtained. A schematic representation of the directed graph inJresistance space is shown in Fig. 2b. The directed graph in the Jresistance space, Fig. 2b, contains information of the crack growth resistance curve for all possible crack paths in the microstructure. How- ever, one needs to know the exact crack path to esti- mate the resultant crack growth resistance. The crack path for a speci ed amount of crack extension in the crack propagation direction can be estimated by nding the path with the least crack growth resistance, i.e. the path of global minimum resistance. The crack path for a speci ed amount of crack extension in the crack prop- agation direction can also be estimated by propagating the crack by a few grain boundary facets at a time based on some local minimum crack growth resistance crite- rion. The path of least resistance for the entire amount of the speci ed crack extension or for the cracks prop- agating a few grain boundary facets at a time can be estimated using Dijkstra's algorithm (Dijkstra, 1959). Dijkstra's algorithm is a graph search algorithm that produces the minimal distance between two nodes. The distances in this case are the increments in Jfor a crack traversing between two nodes. 6 Shmuel Osovski et al. (a) (b) Fig. 2 (a) The tearing modulus, TR, of a crack propagating along a grain boundary oriented at an angle with respect to the plane of the propagating crack (closed squares) and the increment in crack growth initiation resistance, J, resulting from a crack reaching a triple junction as a function of (closed circles). The inset in (a) shows a symmetric triple junction placed directly in front of the initial crack with an opening angle = 45. (b) A schematic representation of the graph describing the change in Jrequired for the crack to propagate through the grain boundary network. The edges of the graph representing grain boundary network are colored by the value of the required Jto propagate through the grain boundary facet. (a) (b) Fig. 3 The (a) crack growth resistance curves and corresponding (b) crack paths predicted using three local minimum crack growth resistance criteria and the global minimum crack growth resistance criterion. In Osovski et al. (2015a) it was shown that, for given material properties and a given grain microstructure, the value of Jstrongly depends on the o set between the initial crack tip and the nearest grain boundary that is not orthogonal to the initial crack. Thus, the value ofJat the initiation of crack growth depends on a single event and can be highly variable. In Fig. 3a and subsequently, we focus on the e ect of microstruc- ture on crack growth following crack initiation so that the values of Jplotted have the value of Jat the ini- tiation of crack growth subtracted. The values of Jassociated with crack growth following initiation de- pend on a collection of events associated with the entire grain microstructure and not just a single event. Fig- ure 3a compares the crack growth resistance curves, J versusa, for crack paths estimated using four crack path selection criteria: (i) the global minimum resis- tance crack path for a speci ed amount of crack ex- tension,a= 11mm, (ii) a local minimum resistance crack path obtained by propagating the crack one grain boundary facet at a time and choosing the propagation path following the value of TRat the adjacent edges, Intergranular Fracture Prediction and Microstructure Design 7 (iii) a local minimum resistance crack path obtained by propagating one grain boundary facet at a time and choosing the propagation path by nding the nearest connected neighbor, and (iv) a local minimum resis- tance crack path obtained by propagating two grain boundary facets at a time and choosing the propaga- tion path based on the shortest path in terms of the change inJ. The criteria (ii) can be thought of as the \most local" crack path selection criteria. The crack paths estimated using the four crack path selection cri- teria are shown in Fig. 3b. The Jacurves obtained using all the four criteria in Fig. 3a are nearly identi- cal fora4:5mm but for a > 4:5mm the global minimum criteria under predicts the crack growth re- sistance. One possibility is because the choice of the incremental crack path at a speci c junction is not bi- ased from the knowledge of what is to come next. An- other possibility is that the unit event model does not account for the e ects of events, such as micro-cracking away from the main crack, that can occur ahead of the main crack tip. 4 Model Prediction and Validation Two sets of crack growth resistance curves, Jversusa, for two realizations of the same material microstructure are compared in Fig. 4. One set is obtained from full microstructure-based nite element calculations while the others were obtained using the unit event model and graph search. The two realizations of the grain mi- crostructure are identical in terms of material proper- ties and microstructural length scales (average grain boundary thickness and grain size), and only di er in the initial array of random generators for the Voronoi tessellation process. In Figs. 4a and b, Jversusa curves are predicted using the proposed unit event model for three crack paths: (i) the crack path obtained from the full eld microstructure-based nite element calcu- lations (model prediction), (ii) the crack path of mini- mum crack growth resistance for a speci ed crack exten- sion,a= 8mm (global minimum), and (iii) the crack path dictated by the local minimum crack growth resis- tance at each grain boundary junction (local minimum). The local minimum crack growth resistance crack path is determined by propagating the crack at each junction along the grain boundary giving the smallest change in the value of Jfor a givenauntil the crack has reached the speci ed a= 8mm. As noted in Section 3, the value of Jat initia- tion i.e. the value at a= 0, depends strongly on the o set between the initial crack tip and the nearest grain boundary facet and is highly variable. Because our focus is on the e ect of microstructure on the crackpath, and, in turn, on the e ect of crack path on crack growth resistance, the values of Jplotted in Figs. 4a and 4b have the value of Jata= 0 in the unit event model and graph search predictions scaled to match the nite element results. The crack paths obtained from the microstructure-based nite element calcula- tions and predicted using the proposed model based on global minimum and local minimum crack growth re- sistance criteria for the two realizations of the speci ed grain microstructure are shown in Figs. 5a and 5b, re- spectively. As shown in Fig. 4a for realization 1, the Jversus acurve predicted using the proposed model for the three crack paths and microstructure-based nite ele- ment calculations are nearly identical for a < 4mm. The model predictions for all three crack paths di er from the results of the microstructure-based nite el- ement calculations for a > 4mm. This is because for realization 1, the nite element calculation exhibits crack branching, see Fig. 5a. Branching acts as addi- tional toughening mechanism that is not account for in our model. On the other hand, the microstructure- based nite element calculation for realization 2 does not exhibit crack branching, Fig. 5b. In this case, the Jversusacurve predicted using (i) the same crack path as in the microstructure-based nite element cal- culations (model prediction) and (ii) the crack path dic- tated by the local minimum crack growth resistance at each grain boundary junction (local minimum) are es- sentially identical to the Jversusacurve obtained from the microstructure-based nite element calcula- tions, Fig. 4b. Furthermore, for realization 2, the pro- posed model based on the global minimum criteria under- predicts the Jversusacurve. This is because there is no a priori reason for a crack propagating in a given microstructure to choose the path of least global resis- tance. 5 An Implication for Microstructure Design The results presented in Section 4 clearly show that the computationally ecient unit event based model can predict the crack growth resistance of material mi- crostructures that are prone to intergranular fracture. Possible routes to designing microstructures that in- crease the crack growth resistance are to increase the extent of crack de ection from the crack propagation direction and/or to promote crack meandering. Both of these increase the amount of new free surface created for a given amount of crack growth a. To quantify the scaling between the crack growth resistance and the ex- tent of crack de ection from the crack propagation di- rection, we generated three realizations of a speci ed 8 Shmuel Osovski et al. (a) (b) Fig. 4 Comparison of crack growth resistance curves obtained from the microstructure-based nite element calculations (FE results) and predicted using the proposed model when following the same crack path as obtained from the nite element calculations (Model prediction), when following the crack path of global minimum crack growth resistance (Global minimum) and when following the path of local minimum crack growth resistance (Local minimum) for two realizations of the speci ed grain microstructure: (a) realization 1 and (b) realization 2. initial microstructure. For each realization the starting point of crack growth is chosen to be the node closest to the symmetry plane y= 0. All of the reachable nodes atx=aandjyj=yare chosen as targets and the crack path as well as the crack growth resistance are cal- culated using the global minimum resistance criterion. Figure 6a shows the variation in the value of Jwith the overall deviation from the symmetry plane, y, i.e. the height of the target node with respect to y= 0, at a speci edafor the three realizations of the speci ed microstructure. A few representative crack paths corre- sponding to the crack growth resistance for realization 1 of the microstructure are shown in Fig.6b. As seen in Fig.6a, the Jvs.ycurves for all three realizations fol- low the same master curve with J=J0+k1yk2, where J0= 62:2 is the overall crack growth resistance for the crack path with negligible crack de ection, i.e. y0. The tting parameters are k1= 1:68 andk2= 2. The dependence of Jony2is likely due to the fact that J depends on the square of the actual crack path, S2, such thatJ/S2=a2+y2. The results presented in Fig.6 highlight the motiva- tion to control microstructures in order to increase the extent of crack de ection. One way to achieve this is by introducing anisotropic grain morphology in which the weak grain boundaries are preferentially inclined at an angle with respect to the cross section of the crack growth direction. To illustrate this and to as- sess the contribution of grain orientation on enhancing crack growth resistance, a grain structure with an as- pect ratio of 1 : 5 in the crack growth direction wasgenerated. Next, the grain microstructure was rotated such that the crack growth direction is oriented at an anglegwith respect to the long axis of the grains. The elongated grain morphology and the speci ed mi- crostructure rotated at an angle of g= 45are shown in Fig. 7b. A comparison of the crack growth resistance curves,Ja, for the grain microstructures shown in Fig. 7a, obtained using the microstructure-based full eld nite element calculations and the model using unit event calculations and graph search are presented in Fig. 7b. As seen in Fig. 7b the crack de ection arising as a result of highly anisotropic grain morphology has a signi cant toughening e ect. In addition, the results also highlight the predictive capabilities of the model using unit event calculations and graph search. In the previous example the focus was on large crack de ection away from the crack propagation direction. The crack growth resistance can also be increased by increasing the extent of crack meandering. The model proposed here only requires a series of discrete events, namely the de ection angle and grain boundary facet length. To illustrate the e ect of crack meandering, we focus on the de ection angle, , and assume a rather uniform grain boundary facet length. In particular, we assume that during the crack propagation process, the angle of the grain boundaries the crack encounters have a known distribution and take the distribution to be Gaussian. The Gaussian distribution is taken to be cen- tered around mean = 60, corresponding to an equiaxed hexagonal structure. Various values of the standard de- viation are considered, = 5, 15,20and the grain Intergranular Fracture Prediction and Microstructure Design 9 (a) (b) Fig. 5 Crack paths obtained from the microstructure-based nite element calculations and predicted using the proposed model based on global minimum crack growth resistance and local minimum crack growth resistance for two realizations of the speci ed grain microstructure: (a) realization 1 and (b) realization 2. facet size to taken be 160 m. The de ection angle dis- tributions were sampled and the obtained angle along with facet size is used to calculate the increment in Jresulting from a propagation event. The process is repeated until achieving a speci ed amount of crack growth (in this case a= 6mm). For each distribu- tion, 10;000 possible trajectories were extracted andthe mean, along with the 95% con dence level were calculated. The results are presented in Fig. 8. Increas- ing the probability of encountering large de ection an- gles leads to a broadening of the 95% con dence in- terval. In addition, the lower bound increases with in- creasing standard deviation . This implies that the larger the deviation is from an equiaxed structure, the 10 Shmuel Osovski et al. (a) (b) Fig. 6 (a) The dependence of overall crack growth resistance, J, for a speci ed amount of projected crack extension, a, along the crack propagation direction ( x-axis) on the overall crack de ection, y, from the crack propagation direction for three realizations of a speci ed grain microstructure. (b) The crack paths corresponding to the crack growth resistance for realization 1 of the grain microstructure. (a) (b) Fig. 7 (a) The elongated grain microstructure with an aspect ratio of 1:5 and rotated at an angle of g= 0and 45with respect to the crack growth direction. (b) Comparison of the crack growth resistance curves obtained using the microstructure- based nite element calculations and predicted using the proposed model when following the path of local minimum crack growth resistance for the two microstructures. tougher the material is. To validate the model predic- tions microstructure-based full eld nite element cal- culations were carried out for a microstructure consist- ing of grains with a mean facet size 160 m and an an- gle distribution with mean = 60and= 17. The crack growth resistance curve obtained from the nite element calculation are compared with the model pre- dictions in Fig. 8. The crack growth resistance curve obtained from nite element calculation falls within the con dence limits of the = 15and= 20. The re-sults presented in Fig. 8 show that the model based on unit events and graph search can also be used to as- sess the overall crack growth resistance of the materials with a known distribution of grain morphology. 6 Concluding Remarks With the emergence of advanced manufacturing tech- nologies such as additive manufacturing (Herderick, 2011; Intergranular Fracture Prediction and Microstructure Design 11 Fig. 8 The 95% con dence envelope for the crack growth resistance curves of a grain boundary network with a facet length ofda= 160mand a Gaussian distribution of de ec- tion angles with a mean of = 60 and a standard deviation of:= 5 (blue); = 15 (green); = 20 (red). The results of a FE simulation for a grain microstructure with the same average facet length and = 17 is shown for comparison (black). Frazier, 2014) or others (Hosokawa et al., 2013) it is be- coming increasingly possible to control the microstruc- tural features of metallic alloys. In addition, the emer- gence of 3D tomography techniques (Uchic et al., 2007; Ueda et al., 2014; Kahziz et al., 2016) make it possible to characterize the microstructural features of a mate- rial in great detail. The capability to characterize and control microstructural features provides a motivation for developing methodologies for designing material mi- crostructures with improved properties. In this work, for materials that are prone to grain boundary ductile fracture, we have developed a methodology for design- ing grain microstructures with improved crack growth resistance. The proposed model is based on discrete unit events and graph search. The discrete unit event associated with grain boundary fracture is the propa- gation of a crack from one grain boundary junction to the another. The \cost" in terms of the change in J, associated with each unit event is used to construct a graph in both Euclidean and Jresistance space. The directed graph in Jresistance space contains informa- tion on the crack growth resistance curve of all possible crack paths in the microstructure. The crack path for a speci ed amount of crack growth is then estimated as the crack path dictated by either global minimum or local minimum crack growth resistance criteria. The predictive capability of the proposed model is validated against results of full eld microstructure-based nite element calculations. The comparison of the results ob- tained from the nite element calculations and thosepredicted using the proposed model highlights the fact that the crack path is largely dictated by the local min- imum crack growth resistance at each grain boundary junction and not by the global minimum crack growth resistance. In a broad range of areas of microstructural design, grain boundaries have attracted considerable attention. Grain boundaries are an extremely important hetero- geneity in many structural alloys and metals that can be controlled (in a statistical manner) both in terms of spatial arrangement and grain boundary character. The ability of polycrystalline materials to resist failure under extreme environmental and mechanical loads is strongly correlated with the grain boundaries network they poses (Kim et al., 2003; Kobayashi et al., 2014, 2016). Stress corrosion cracking, hydrogen embrittle- ment and low damage tolerance are all dominated by the grain boundaries (Kim et al., 2006; Kobayashi et al., 2012; Ara n and Szpunar, 2009). Over the past decades, our understanding of the re- lation between speci c grain boundaries and their re- sistance to the aforementioned processes has increased along with the ability to choose processing routes lead- ing to speci c grain boundary character distributions, geometrical properties of the grains and chemical pa- rameters (Randle and Owen, 2006; Cantwell et al., 2016). However, designing an optimum grain boundary net- work to improve material properties involves optimiz- ing the topological features of the material microstruc- ture. The (limited) available computational techniques are extremely time consuming so that direct coupling of such calculations with an iterative optimization scheme is not practical. We hope that our computationally in- expensive model can be used in an iterative, non-gradient topological optimization scheme to design material mi- crostructures with improved intergranular fracture re- sistance. Additionally, the model can be used to predict the path and resistance of any crack growth process that discretely propagates from one grain boundary junction to another, or to predict the overall crack growth resis- tance of materials with a known distribution of grain boundary orientations and grain boundary lengths. Acknowledgements The nancial support provided by the Pazy foundation young researchers award grant # 1176 (SO), U.S. National Science Foundation grant CMMI - 1663130 (AS), and European Union's Horizon2020 Programme (Ex- cellent Science, Marie-Sklodowska - Curie Actions, H2020 - MSCA - RISE - 2017) under REA grant agreement 777896 (Project QUANTIFY, SO and AS) are gratefully acknowl- edged. We are also grateful for the high performance research computing resources provided by Texas A&M University. 12 Shmuel Osovski et al.