Early prediction of macrocrack location in concrete, rocks and other granular composite materials

Heterogeneous quasibrittle composites like concrete, ceramics and rocks comprise grains held together by bonds. The question on whether or not the path of the crack that leads to failure can be predicted from known microstructural features, viz. bond connectivity, size, fracture surface energy and strength, remains open. Many fracture criteria exist. The most widely used are based on a postulated stress and/or energy extremal. Since force and energy share common transmission paths, their flow bottleneck may be the precursory failure mechanism to reconcile these optimality criteria in one unified framework. We explore this in the framework of network flow theory, using microstructural data from 3D discrete element models of concrete under uniaxial tension. We find the force and energy bottlenecks emerge in the same path and provide an early and accurate prediction of the ultimate macrocrack path \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {C}}$$\end{document}C. Relative to all feasible crack paths, the Griffith’s fracture surface energy and the Francfort–Marigo energy functional are minimum in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {C}}$$\end{document}C; likewise for the critical strain energy density if bonds are uniformly sized. Redundancies in transmission paths govern prefailure dynamics, and predispose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {C}}$$\end{document}C to cascading failure during which the concomitant energy release rate and normal (Rankine) stress become maximum along \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {C}}$$\end{document}C.


List of symbols α(e)
Bond capacity of link e α E (e) Bond capacity for energy of link e α F (e) Bond capacity for force of link ē α Residual bond capacitȳ α E (ǫ) Residual bond energy Ŵ An arbitrary, feasible crack path (cut) in N Ŵ min Minimum cut in N Ŵ max Maximum cut in N γ (e) Griffith's fracture surface energy of link e γ (Ŵ) Griffith's fracture surface energy for path (cut) Ŵ γ (Ŵ) Griffith's fracture surface energy density for path Ŵ ν Poisson's ratio of contact ν ITZ Poisson's ratio of cement-ITZ or ITZ-ITZ contacts ν CM Poisson's ratio of cement-cement contacts µ Inter-particle friction angle µ ITZ Inter-particle friction angle for cement-ITZ or ITZ-ITZ contacts µ CM Inter-particle friction angle for cement-cement contacts δ + (v) Set of arcs emanating from v δ − (v) Set of arcs terminating at v ǫ Time stage �(Ŵ) Energy release rate at Ŵ  Radius of the smaller of the two grains sharing the bond r A , r B (r = r A ≤ r B ) Radii of the two grains sharing the bond R + 0 Non-negative real numbers S(Ŵ) Edges in N that are saturated s Source node t Sink node T n (e) Critical normal tensile stress of bond or link e T ITZ n Critical normal tensile stress for cement-ITZ or ITZ-ITZ contacts T CM n Critical normal tensile stress for cement-cement contacts U(e, ǫ) Potential energy of bond or link e at loading stage ǫ u F Speed of applied loading force V Set of nodes v Node V B u Set of grains in the upper part of the bottleneck V C u Set of grains in the upper part of the macrocrack W, W ′ Two distinct node sets of V V T All the grains in the sample Force and strain energy share common pathways for flow. For many everyday quasibrittle materials like concrete, ceramics, asphalt mixtures, and rocks, these pathways are embodied in a disordered network of intergranular bonds with different geometric and constitutive properties [1][2][3][4] . While there is broad recognition that the bulk strength and fracture of granular composites are governed by subscale heterogeneities, the details of these dependencies remain mired in controversy 1,5,6 . Many engineering challenges depend on a detailed understanding of these microstructural-to-bulk relationships. In the built environment, the growth in demand for more resilient and sustainable materials continues to outpace development 7 . Although only 2-3% of the Earth's landmasses, cities consume 60-80% of its energy, produce more than 75% of its greenhouse gases, and consume materials at Scientific Reports | (2020) 10:20268 | https://doi.org/10.1038/s41598-020-76616-y www.nature.com/scientificreports/ a level that is projected to reach 90 billion tonnes by 2050 8 . These trends have seen a push for technologies that can harness disorder and heterogeneities in rational design and engineering of high performance and green materials 3,4,6,[9][10][11] . Great impetus for improvements in early detection of damage in materials and structures has also arisen, with the confluence of reduced cost and advances in nondestructive microstructural sensing technology, the rapid growth in construction and manufacturing activities, and the emergence of continuous structural health monitoring to meet government regulations on new and aging infrastructure 12,13 . These challenges have prompted calls for a more nuanced complex systems view of materials, which couple experimental observations with data-driven modeling and simulation, especially in the archetypal composite granular material like concrete 7 . Here, on the centennial anniversary of the founding Griffith's theory for fracture 14 , we do so in a manner complementary to traditional fracture mechanics. Specifically, this study seeks to answer: Can the path of the crack that leads to failure be predicted from the connectivity of transmission paths and their microstructural fracture properties and, if so, how far in advance? To achieve this, we break with the tradition of modeling fracture propagation and attendant energy and force transmission as flows in a continuum 1,6,15 . Instead, we model these as flows in a network [15][16][17][18][19][20] . To the best of our knowledge, this is the first study that: (i) delivers a method for early and accurate prediction of the ultimate macrocrack location from known microstructural bond properties (i.e., fracture strength and surface energy, size), disorder in the network connectivity (an effect of packing and grain size distribution), and stage of loading; and (ii) reconciles and consolidates, in a single framework, the precursory failure mechanism of energy and force flow bottlenecks with widely used fracture criteria. With this in mind, this study may help generalize fracture mechanics theory to complex materials where discreteness, disorder, path redundancy and heterogeneity are not only the norm rather than the exception-but are the salient features that govern bulk strength and failure.
The approach of modeling energy and force transmission as flows in a network holds several advantages. First it lends well to proven tools for modeling and characterization of transmission dynamics germane to capacitated complex networks (e.g., power grid, telecommunication and road networks etc.) 21,22 . In such networks, each component (e.g., link) has a finite capacity for flow, and the most vulnerable part of the network to cascading failure is the so-called bottleneck where congestion occurs. Flow bottlenecks are a fundamental concept of network resilience and have been considered as precursors for cascading failure in various cyber, transport and infrastructure systems 23 . Prior evidence also demonstrate that these emergent structures can reliably predict, early in the prefailure regime, the ultimate locale of cascading failure for various natural and synthetic granular media 16,[17][18][19][20] . Second, the problem of finding where the flow bottleneck emerges and how this path is influenced by heterogeneities in transmission pathways can be expressed in terms of an optimization problem. Consequently, this opens up the opportunity to consolidate network flow approaches with widely used fracture mechanics criteria, which are similarly couched in terms of a postulated energy or stress extrema. Third, as shown recently in the work of Patel et al. 24 on complex biostructures and van der Linden et al. 17 on porous granular media flow, synergies between continuum finite element methods and discrete network flow analysis can produce insights on subscale heterogeneities not otherwise accessible from a purely continuum mechanics approach.
Heterogeneity in transmission pathways has major implications for the spread of failure in a broad range of networked systems 21,22 . In quasibrittle granular composites, this was recently investigated with respect to the coupled evolution of force and damage propagation in the framework of network flow theory 16 . It showed that the majority of force chains develop in the most direct (shortest possible) routes for force transmission, while the ultimate macrocrack emerged in the recurrent force bottleneck that persisted from the nascent stages of the prefailure regime. Different from the work in 16 , here we study the general case of energy, force and their residual transmission processes; the endmost deals with the proximity to breakage of bonds, viz. the amount of energy or force flow that each bond could still transmit given its capacity and current flow.
We proceed in four steps: measure, summarize, predict and characterize ( Fig. 1). In Step 1, the input data to the analysis is obtained by measuring the microstructural properties of various samples at distinct stages of a uniaxial tension test. This is achieved through a combined DEM simulation and experimental test campaign which builds on prior investigations on fracture in concrete [24][25][26][27][28][29][30] (see also Supplementary file). 3D DEM simulations of real physical tests on concrete were conducted, incorporating the real aggregate size and shape distribution, positions and volume (Fig. 2); these included uniaxial compression and tension, three-point bending and four-point bending tests to ensure the simulations can reproduce realistic quasibrittle fracture patterns in random heterogeneous 3-phase materials composed of aggregate particles, cement matrix and interfacial transitional zones (ITZs), as well as in 2-phase cement composites with no ITZs. In this study, for simplicity, we use data from uniaxial tension test simulations where the particles (aggregate and mortar) are spheres in planar (one grain thick samples) and fully 3D DEM samples. In Step 2, we summarize the properties of the sample: network connectivity and individual bond properties (i.e., fracture strength, fracture surface energy, size, realized tensile force and strain energy). This summary is compiled into an input data to Step 3, where a data-driven and multiscale network flow analysis is undertaken to predict the ultimate macrocrack path C based on a prescribed fracture criteria across the different loading stages from pre-peak to the post-peak softening regime. Eight fracture criteria are examined with respect to their ability to deliver an early prediction of C in the prefailure regime. Criteria 1 and 2 correspond to the bottleneck of the strain energy flow network and its residual, respectively; while Criteria 3 and 4 focus on the bottleneck of the tensile force flow network and its residual. We recast into a network flow formulation four widely used fracture criteria: minimum energy functional of Francfort and Marigo 31 , minimum critical strain energy density 15 , maximum energy release rate 32 , maximum normal stress (Rankine stress) 33 . In Step 4, assuming specific conditions hold, we present closed-form relationships between the eight criteria which, in conjunction with the numerical predictions from Step 3, characterize explicitly the influences of disorder in the network connectivity, bond heterogeneities in fracture strength and surface energy, and stage of loading on C . Our hypothesis is that since force and energy are correlated and share common transmission paths, their flow  www.nature.com/scientificreports/ bottleneck may be the precursory failure mechanism that can reconcile and consolidate these optimality criteria in one unified framework.

Theory
Construction of the flow networks. At each stage of loading ǫ , we construct three flow networks F , F 1 and F from the bond contact network N . The flow network F = (G, s, t, c) is given in terms of: which is a directed graph that consists of a set of nodes v ∈ V and a set of arcs e ∈ A ; a source node s; a sink node t; and a capacity function c : A → R + 0 . A flow is a function that assigns non-negative real numbers to the arcs of G, φ : A → R + 0 , subject to the following conditions.
Capacity constraint: for every arc, e ∈ A , the flow cannot exceed the capacity Conservation of flow: for every node except the source and sink, v ∈ V − {s, t} , the amount of flow entering a node is the same as that leaving the node viz.
where δ − (v) is the set of arcs terminating at v, and δ + (v) is the set of arcs emanating from v.
Next we relate F to the bond contact network N by assigning the nodes v to represent the grains while the arcs e represent the bonds (Fig. 3). Each bond in N has two properties: a bond flow f which changes with the stage of loading and a fixed capacity α . Both f and α are known from the DEM data and form part of the input to the network flow analysis. We consider two types of flow: (i) energy flow, (ii) tensile force flow. These flows can be reasonably assumed to obey the above two conditions for flow 16 .
For energy transmission, the bond flow corresponds to the potential energy in the bond U: f (e, ǫ) = f E (e, ǫ) = U(e, ǫ) . The bond capacity is given by the Griffith's fracture surface energy γ 14 : where K n = E 2r A r B r A +r B is the normal contact stiffness and F n min is the bond strength, the force needed to break the bond, given by www.nature.com/scientificreports/ E is the Young's elastic modulus of the contact, T n is the critical normal tensile stress for the contact, b is the bond contact area b = πr 2 A , and r A and r B are the radii of the grains sharing the bond such that r A ≤ r B (for convenience, we refer to r A simply as r in the rest of the paper). For force transmission, the bond flow is the tensile force magnitude f (e, ǫ) = f F (e, ǫ) = ||F n (e)|| , while the capacity is given by the bond strength α(e) = α F (e) = F n min (e) . All the material properties have been calibrated against the physical experiments of van Vliet and van Mier 34 (see Nitka and Tejchman 35 for complete details). In all the samples, we assume all bonds are in tension and no new contacts are created, following Tordesillas et al. 16 A feasible crack path (cut) of F can be defined as two distinct partitions {W, W ′ } of V such that s ∈ W and t ∈ W ′ . Every cut of F determines a set of links Ŵ of G, where Ŵ contains all arcs emanating from a node in W and terminating on a node in W ′ . The capacity of such a cut Ŵ is defined as c(Ŵ) = e∈Ŵ α(e). Here we are interested in finding cuts Ŵ of F which satisfy specific optimality conditions on their capacity c(Ŵ) . We focus on the following two types of cuts of F = (G, s, t, c): This cut corresponds to the bottleneck of the flow network, B(F(G, s, t, c)) = Ŵ min . This is obtained by solving the Max-flow Min-cut Problem using the Ford-Fulkerson algorithm 36 .
The problem of finding an optimal solution is NP-hard, so we use the probabilistic heuristic optimization algorithm of simulated annealing to solve this approximately for F 37 .
In addition to F , we define two additional flow networks F 1 and F with the same topology as F but which differ in their capacities as follows: unit capacity flow network F 1 = (G, s, t, 1) and the residual flow network F = (G, s, t,ᾱ) where ᾱ = α − f . Criterion 1: The energy flow bottleneck as the path of minimum Griffith's fracture surface energy 14 . Griffith postulated that a crack will propagate when there is sufficient stored potential energy to overcome the fracture surface energy of the material. The outstanding question is: Once this condition is met, along which path will the ultimate macrocrack propagate? Here we propose that this path is given by B(F E ) , the bottleneck of the flow network F E = (G, s, t, γ ) : the cut where the Griffith's fracture surface energy γ (Ŵ) is minimum over all feasible cuts Ŵ of F E . Criterion 2: The residual energy flow bottleneck. Criterion 1 can be extended to account for the state of flow in the system at ǫ , by incorporating the realized strain energy for the bond U(ǫ) in the capacity function: ᾱ E (ǫ) = γ − U(ǫ) . Thus an alternative criterion for the crack path is that path which is closest to breaking point relative to the fracture surface energy of the system: the bottleneck B(F E ) of the flow network  proposed that the macrocrack develops along the path where the total energy functional E tot is minimum, where E strain (N \ Ŵ, ǫ) is the sum of elastic strain energies stored in all the bonds in N excluding those in the virtual cut path Ŵ , and γ (Ŵ) is Griffith's fracture surface energy (or energy dissipated due to crack formation) of the cut Ŵ . Equation 3 can be expressed as (3)  1 can instead be expressed in terms of the critical strain energy density in the spirit of Sih 15 : α E (e) =γ (e) = γ (e) π(r(e)) 2 . This yields a new potential crack path which is given by the minimum cut of the flow network F = (G, s, t,γ ) . Thus, the path where the Griffith's fracture surface energy density γ (Ŵ) is minimum is given by the minimum cut of F (G, s, t, γ πr 2 ) , B F (G, s, t, γ πr 2 ) . Here we present Step 4 (recall Fig. 1), where we assume that conservation of flow is exactly satisfied which implies that: (i) the total flow across any two feasible cuts of a flow network is the same; (ii) the bottleneck of any two flow networks with capacities that differ only by a multiplicative constant is the same. Below, subject to these two conditions, we establish closed-form relationships between the energy flow bottleneck B(F E ) from Criterion 1, the force flow bottleneck B(F F ) from Criterion 3, and the paths corresponding to all the other fracture criteria. Note that corollaries A-E apply to all stages of loading history, whereas corollary F is confined to the post-peak softening regime. C. If the bonds are uniformly sized, then the path with the minimum critical strain energy density is the path with the minimum fracture surface energy.. If the bonds are uniformly sized, then the two flow networks F (G, s, t, γ π r 2 ) and F (G, s, t, γ ) will have the same bottleneck as their capacity functions differ only by a constant multiple. That is, the critical strain energy density is minimized along B F E from Criterion 1.
D. When the bonds have the same critical tensile bond stress, then the path that satisfies the Rankine stress criterion is the path with the minimum fracture strength.. Observe that, regardless of stage of loading, the path which maximizes the tensile normal stress in F F (G, s, t, F n min ) is that path with the least total bond area since all cut flows in the flow network are the same. This path is given by B F F (G, s, t, b) , the minimum cut of the flow network F F (G, s, t, b) with edge capacity equal to the area of the corresponding bond b. Since F n min (e) = T n (e) π b(e) , in general, B F F (G, s, t, F n min ) = B F F (G, s, t, b) . But if the critical tensile normal contact stress T n is the same for every bond in N , then the capacity functions will be different only by a multiplicative constant. Thus  www.nature.com/scientificreports/ s, t, b) . That is, when the critical tensile stress T n is uniform in N , the path that satisfies the Rankine stress criterion is the force bottleneck B F F from Criterion 3.

E. When bond properties are uniform, and the only source of heterogeneity is the bond connectivity-viz. disordered but otherwise homogeneous network-then all the fracture criteria except for the maximum release rate (Criterion 7) are satisfied along that path with the least number of bonds..
In the special case when bonds in N are homogeneous, the capacity functions α E and α F are both constants and B(F E ) = B(F F ) . Moreover, the paths from Criteria 1-6 and 8 converge to one path, which is given by the minimum cut of F 1 (G, s, t, 1) , the flow network with unit capacity. This path, the so-called minimum edge cut, is the cut with the least number of bonds 38 . The number of bonds in this cut, p min , gives a measure of path redundancy in N : the minimum number of all available flow pathways (joint and disjoint) between the top and bottom walls of the specimen.

Results and discussion
We sought to address an open question: Can the path of the crack that leads to failure be predicted from the connectivity of transmission paths and their microstructural fracture properties and, if so, how far in advance? That the energy and force bottleneck distinguish themselves from other feasible partitions of the specimen-even before the onset of damage-suggests this is possible (Figs. 5, 6 and 7). For each criterion, we quantified the error between the predicted location and the actual macrocrack path C as follows: where V B u , V C u are the set of grains in the upper part of the bottleneck and the macrocrack respectively, V T is the set of all the grains in the sample, and | | denotes the cardinality (number of grains). Hence P quantifies the number of offset grains in the prediction relative to the macrocrack, normalized by the total number of grains. www.nature.com/scientificreports/ 1. We showed that the bottlenecks of the energy and force flow network and those from their respective residual networks are essentially the same. This implies that at any given stage of loading, the bottlenecks distinguish themselves from all other feasible crack paths in two respects: by having the least total capacity (Griffith's fracture surface energy in F E , or fracture strength in F F ), and by being closest to this capacity. That is, the fracture criterion that delivers an early prediction of C must be related to the flow bottleneck, the path of least resistance to fracture, where the force and energy flows are also closest to their critical fracture values. On the other hand, the realized (actual) flows are an outcome of an optimized transmission process that maximizes the global throughput subject to the redundancies of, and the capacities of member segments or bonds along, the percolating paths through the specimen. Results here are consistent with recent experiments on ceramics 3, 4 , metallic glasses 6 and graphite 39 .    Fig. 1). 3. The paths of minimum critical strain energy density, maximum energy release rate and maximum normal stress (Criteria 6-8) do not provide a reliable early prediction of the failure location. This is due to path redundancy in N in the prefailure regime, which is not captured by any of these fracture criteria. Redundancies in transmission paths are what enables energy and force to be rerouted to higher capacity bonds to prevent bond breakage. As shown recently 16 , damage in the early stages of of the prefailure regime mainly occurs away from the macrocrack path C . Similar stress redistributions around strain concentration sites have been observed to divert damage away from the crack path (viz., the crack tip) into the bulk of the sample 9,41 . Given that C essentially coincides with the energy bottleneck B(F E ) , this may initially seem counterintuitive given the close proximity to fracture of the bonds in B(F E ) . It turns out there is an interplay between force chain stability in the energy bottleneck and path redundancy, as we now demonstrate. www.nature.com/scientificreports/ • We first quantify path redundancy through p min 36 , the minimum number of available transmission pathways between the top and bottom walls of the specimen (Fig. 8). p min is purely a topological property of N and is equal to the number of percolating link-disjoint paths (paths that have no common links) through N between the source and the sink. Broadly, p min is one of the fundamental concepts in measuring resilience in transmission networks, viz. the ability of a network to withstand 'attacks' (loss of links or nodes). For instance, in road networks, p min may be used to quantify the extent to which cars can be rerouted to alternative paths when a road is closed to traffic (e.g., roadworks or accident). The monotonic decrease in p min reflects the degradation of flow paths caused by the spread of damage (Fig. 3). Damage disrupts the transmission of energy and force through the specimen since the breakage of a bond disconnects percolating paths for flow. In the presence of redundant paths, this disruption triggers a reroute or redistribution of flow. We highlight this process in Fig. 8-inset in terms of the optimized force route P , the most direct or shortest possible percolating routes for force transmission, where most force chains develop 16 . As seen in Fig. 8-inset, the system recovers in the presence of damage. New links in P replace old links therein which either break or can no longer be accessed. When p min is relatively high, just after the onset of damage in N (stages 6-7 for D1; stages 31-32 for D2), the system compensates by replacing contacts that can no longer be accessed in P (numbers of red and blue contacts are roughly balanced). By contrast, in the transition to the post-peak softening regime (stages 9-10 for D1; stages 45-46 for D2), P rapidly degrades without a matching recovery: observe the surge in population of blue contacts as that of feasible replacement contacts in red dwindle. • Remarkably, while the bottleneck force chains are the most obvious suspect locales for incipient failure, we find these endure during the prefailure regime, contributing mainly to the accumulation of stored energy and consequent energy release in the terminal phase of rapid fracture (Fig. 9). Observe in Fig. 9inset that prefailure damage is confined to low capacity links, suggesting that forces are spread out across member contacts (Fig. 10a-d). As previously demonstrated 16 , there is a process of cooperative evolution between the preferred paths for damage (i.e., bottleneck) and the preferred paths for force transmission (i.e., force chains) that underlies material robustness in the prefailure regime. Damage being confined to low capacity links in the bottleneck means that not only is the reduction in the global transmission capacity effectively minimized but also that high capacity contacts are left to support the remaining tensile force chains in the bottleneck. This coevolution among preferred paths for damage and force is an example of the so-called compromise-in competition between dominant mechanisms in the mesoregime of complex systems 42 . Evidence here suggests that this may be a contributing factor to the so-called 'crack shielding' or crack toughening that has been reported in the literature on granular composites 9,41,43 . • The presence or lack of ITZs has no effect on the outcome of the analysis. When ITZs are present, the case with D2, we found that 60% of the contacts in the bottleneck are ITZ. This is consistent with prior studies which showed that ITZs tend to form attractors for macrocracks [27][28][29][30] . That is, they present paths www.nature.com/scientificreports/ of least resistance for crack growth, which evolve to ultimately span the sample through coalescence of ITZ microcracks (interconnected broken links along ITZ contacts) 24-30 . 4. While the antecedent dynamics described above protects the bottleneck from damage in the prefailure regime, it inevitably elicits the opposite effect by predisposing the bottleneck to rapid fracture in the softening regime. As path redundancy progressively drops with the spread of damage-the remaining bottleneck bonds are collectively brought closer and closer to their respective breaking points (Fig. 10b). A cascade of bond breakages then ensues, precipitating the abrupt transition to post-peak softening regime, with two attendant mechanisms in the bottleneck B(F E ): • energy release rate becomes maximum (Fig. 9). This burst to a peak in the energy released in the transition to, and during, post-peak softening can be explained by the remaining high capacity bonds of tensile force chains and their lateral supports, since these store the highest levels of potential strain energy while being closest to fracture (Fig. 10). In contrast, concurrent energy release rates in all other feasible crack paths are much less, due to member bonds being further away from their respective fracture surface energies and hence are able to support an increase in strain energy without breaking. Note the higher energy released during failure cascade in D2 compared to that in D1, evident in Fig. 9, due to the higher pathway redundancy in D2 (Fig. 8). Tensile force chains in D2 are significantly more supported (higher connectivity means higher redundancy) than those in D1 16 , in turn enabling the D2 bottleneck to store higher strain energies in the stages preceding rapid fracture. This influence of connectivity on the stability of force chains was recently observed by Tang et al. 43 in their study of packing and grain size distribution for optimal performance of reactive materials like aluminumpolytetrafluoroethylene (Al-PTFE) granular composites. • normal stress becomes maximum as the tensile force flow through the bottleneck becomes close to the maximum force flow while the total bond area becomes minimum.

5.
In the previous section, we derived closed-form relationships based on the result that force and energy bottlenecks provide an early and accurate prediction of the ultimate crack path C (step 4 in Fig. 1). These relationships naturally shed light on the specific influences of heterogeneities in bond geometry, strength and surface energy, topological disorder in the bond network, as well as the stage of loading on C , as outlined in the previous section. Additionally, it can help explain why all eight fracture criteria approximately agree in their prediction of C in the post-peak softening regime. That is, the collective breakage of bonds along C in  www.nature.com/scientificreports/ this regime results in C being identified: by criterion 8 due to corollary F, by criterion 7 due to corollary D combined with C becoming the path with the least total bond area, and by criterion 6 since C also becomes the path with the least the number of bonds.
To conclude, we have shown that the energy or force bottleneck, the path of least resistance to fracture, provides an early prediction of the path of the mode I crack that leads to failure in heterogeneous quasibrittle granular materials like concrete. The relationships between the flow bottlenecks and the various paths determined from key fracture criteria (viz. Minimum strain energy density, Maximum energy release rate and Maximum normal (Rankine) stress) highlight the salient influences of disorder in the bond network, heterogeneities in bond properties and stage of loading on the ultimate crack path. Building on recent work 24-30 , a direct comparison between flow bottlenecks and crack paths from experiments with more complex loading conditions (e.g., 3-point and 4-point bending tests) as well as cracking mechanisms in other materials like clay is now the subject of an ongoing investigation. Upscaling the method from laboratory to field is also being explored using a wide range of data sets, including proxies for force and energy 19,20 , with a view toward developing tools for practical decision-making, especially in material design and in Early Warning Systems (EWS) for failure hazard monitoring in natural and man-made structures 44 .

Input data and DEM simulations
Data for the virtual samples D1-D3, all submitted to uniaxial tension, came from a family of discrete element (DEM) models of fracture in concrete (D2) and other granular composites (D1, D3) [24][25][26][27][28][29][30]35 . These models comprise planar (1 grain thick) and fully 3D samples for 2-phase cement composites (aggregate, cement matrix), 3-phase concrete (aggregate, cement matrix, interfacial transitional zones (ITZs)), using the explicit 3D spherical,   27 . Good agreement between numerical and experimental results on real concrete was achieved. D1 is from a "dog-bone" shaped specimen (2-phase): 150 mm high and 100 mm wide (60 mm at the midheight). It consists of 4942 spherical grains: 704 aggregate grains (diameter range of 2-10 mm) and 4238 cement matrix grains (diameter range of 0.5-2 mm). D2 is from a rectangular concrete specimen: 150 mm high and 100 mm wide with two diagonally opposite U-shaped notches. Each notch is of size 15 mm × 5 mm: the notch on the left (right) boundary is 50 mm (100 mm) from the bottom boundary. The specimen is modeled as a 3-phase material composed of aggregate, cement matrix and ITZs. Aggregate grains possess ITZs which are simulated as contacts between aggregate and cement matrix grains. The cement matrix grains have no ITZs. The ITZs are weaker by 30% than the cement matrix. There are 200 aggregate spherical grains (diameter range of 2-16 mm) and 8000 cement matrix spherical grains (diameter range of 0.25-2 mm). The 3D specimen D3, here 2-phase for simplicity, is 100 mm high, 100 mm wide (60 mm at the mid-height) and 100 mm deep. There are 44,386 aggregate spherical grains (diameter range 2-10 mm) and 199,148 cement matrix spherical grains (diameter range of 0.75-2 mm). Overall, the grains comprise 95% of the specimen in D1, D2 and D3. All the material parameters were calibrated against the experiment of van Mier 47 (Table 1).
We investigated the influence of the different parameters and ITZ properties on the macrocrack geometry and location in prior studies 26,30 (see also the Supplementary file). For completeness, a brief summary of the key findings is given in Fig. 11 using D1 as the reference sample. As shown, the DEM parameters, viz., loading speed u F , Poisson's ratio ν , inter-particle friction angle µ and cohesive contact stress C, have little to no influence Figure 11. (Color online) Influence of various fracture properties and loading conditions on the crack geometry and location with D1 as the reference specimen: (a-f) without ITZ, (g) with ITZ. A list of parameter values for samples D1-D8 is given in Table 1.  The parameters are the same as that for D6 except T ITZ n = 18.75 MPa D8 The parameters are the same as that for D6 except T ITZ n = 12.5 MPa Scientific Reports | (2020) 10:20268 | https://doi.org/10.1038/s41598-020-76616-y www.nature.com/scientificreports/ on the macrocrack trajectory under uniaxial tension test conditions (Fig. 11a-f). By contrast, the presence of ITZs has a significant influence on the fracture pattern (Fig.11g). Note that in our past experimental studies of ordinary concrete, we observed ITZs on the concrete surface by means of scanning electron microscope (SEM) and nano-indentation tests 30,48,49 . They exist adjacent to aggregates and have a width of about 20-100 µ m; when compared to the cement matrix, they reveal pronounced compositional differences which are strongest next to the aggregate surface and gradually diminish with distance away from the aggregate interface, though negligible beyond 15-100µ m. Although debate continues about the properties and mechanical effects of ITZs 11,50,51 , we found these to comprise more and larger pores, smaller particles, and a reduced stiffness and strength relative to the bulk phase 27,29 . Thus, they cannot be ignored in the DEM simulations.

Data availability
The data in this paper are available upon reasonable request to Michał Nitka and Jacek Tejchman.