Graph theoretical approaches for the characterization of damage in hierarchical materials

TL;DR

This paper explores how graph theory can characterize damage in hierarchical materials, introducing novel indicators like eigenvector localization and topological dimension to predict failure, especially where traditional methods fall short.

Contribution

It presents new graph-theoretic indicators for damage detection in hierarchical materials, addressing failure modes not well captured by existing percolation or crack growth models.

Findings

Eigenvector localization increases near failure

Topological dimension drops as damage progresses

Indicators can predict failure without avalanche activity

Abstract

We discuss the relevance of methods of graph theory for the study of damage in simple model materials described by the random fuse model. While such methods are not commonly used when dealing with regular random lattices, which mimic disordered but statistically homogeneous materials, they become relevant in materials with microstructures that exhibit complex multi-scale patterns. We specifically address the case of hierarchical materials, whose failure, due to an uncommon fracture mode, is not well described in terms of either damage percolation or crack nucleation-and-growth. We show that in these systems, incipient failure is accompanied by an increase in eigenvector localization and a drop in topological dimension. We propose these two novel indicators as possible candidates to monitor a system in the approach to failure. As such, they provide alternatives to monitoring changes in the precursory avalanche activity, which is often invoked as a candidate for failure prediction in materials which exhibit critical-like behavior at failure, but may not work in the context of hierarchical materials which exhibit scale-free avalanche statistics even very far from the critical load.

Figures

Equations

$$\sum _{i=1}^N{\mathrm{A}}_{ij}(V_i-V_j)=0,$$

$$\mathrm{L}\mathbf{V}=0,$$

$$\mathcal{L}\mathbf{v}=\mathbf{b},$$

$$\mathbf{v}=\sum _{m=1}^{N^{\prime }}\frac{1}{{\lambda }_m}({\mathbf{w}}_m\cdot \mathbf{b}){\mathbf{w}}_m,$$

$$I_m=\frac{{\sum }_i(w_{mi}{)}^4}{{\sum }_i(w_{mi}{)}^2},$$

Full text

1 Introduction

Materials with hierarchical microstructures are characterized by microstructure patterns that repeat on different length scales in a self-similar fashion. This type of organization is encountered in nature in systems such as collagen [1], bone [2, 3], wood [4], and nacre [5, 6], and is believed to provide enhanced toughness, contain damage and impede crack propagation. Recent numerical results for hierarchical fuse networks confirm that the hierarchical organization is responsible for the systematic confinement of fracture patterns [7]. The resulting crack profiles differ from the self-affine scaling behavior that is found in experimental studies for a wide range of materials from metals and glasses[8, 9] to geomaterials [10] and paper [11] as well as in theoretical investigations [12, 13] which represent fracture using random fuse network (RFN) or random beam network models. In hierarchical materials, by contrast, the fluctuations of the fracture surface profile cannot be described in the sense of stochastic roughening and quantified by standard roughness exponents. Such materials can be modelled in terms of hierarchical fuse networks (HFN) or hierarchical beam networks. For different methods of constructing such networks, see [7]. In HFN, the fracture surfaces exhibit discontinuities in the form of abrupt jumps which are power law distributed in size. At the same time, in such hierarchical systems the fracture precursor activity which characterizes damage accumulation before failure exhibits scale-invariant behavior at every loading stage. In non-hierarchical RFN, microcrack accumulation proceeds in avalanches whose power law distribution has an exponential cut-off which diverges at the peak current $I_{\textrm{p}}$. This allows one to envisage failure as a critical phenomenon and large avalanches as precursors of failure. In HFN, by contrast, scale-invariant avalanche-size distributions are encountered well before $I_{\textrm{p}}$, suggesting that the statistical signatures of damage accumulation do not undergo any qualitative changes as one approaches the peak load. The same is true for the post-peak behavior in fracture simulations that are carried out under voltage control: also in this regime, microcracking of HFN occurs in bursts which continue to exhibit power law size distributions and no apparent cut-off. We may thus speak of generic scale-invariant behavior that, unlike critical behavior, is not contingent upon fine tuning the system to a critical point (in RFN: the peak load $I_{p}$). Figure 1 illustrates the difference between RFN - mimicking statistically homogeneous disordered materials - and HFN mimicking hierarchically structured materials, both in terms of the load dependent avalanche statistics and the characteristic damage and stress pattern during failure.

The observation of generic scale invariant behavior in hierarchically patterned networks is by no means new. Activity propagation in the brain, which also exhibits a hierarchical structure, is in fact believed to benefit from such generic scale invariance in dynamic behavior [14, 15]. In the case of materials, generic scale invariance may however imply that the precursor activity cannot be easily used to assess the proximity to failure. The behavior is in fact analogous to the properties of percolation transitions in hierachical networks. In this case too, hierarchical networks do not display a clear percolation threshold [16, 17]: upon removing a fraction $p$ of links, power law distributions of connected component sizes are encountered for every value of $p$, suggesting that connected component statistics is not a suitable tool for assessing the state of damage, and even less so, the resilience of a hierarchical material.

Here we aim at filling the conceptual gap that results from the inadequacy of standard methods for the study of damage in hierarchical materials. To do so, we resort to methods borrowed from network science and graph theory. The idea of using networks to model materials is of course not new. Networks have been used in the past to describe force chains e.g. in granular materials or contact patterns in colloidal dispersions. Models of hierarchical materials may benefit from recent advances in the application of concepts of spectral theory and topology to the study of hierarchical networks. In the following we will focus in particular on the role of eigenvalue localization and Hausdorff dimension, which we believe may provide better indicators of the state of damage in such systems than conventional indicators such as divergence of avalanche sizes.

2 Hierarchical fiber networks

Hierarchical fiber network models are characterized by two main features: i) a modular organization, where modules are densely connected components which exhibit weaker mutual connections; ii) a hierarchical encapsulation, where smaller modules are embedded into larger less connected modules, in an iterative fashion. It was shown that hierarchical fuse networks (HFN) of this type exhibit generic critical-like behavior, regardless of their deterministic or randomized nature [7]. The distinctive feature of a sheet material with hierarchical load carrying structure is that the material is, parallel to the loading direction, sub-divided by cuts of power law distributed lengths, which partition the sample into load carrying modules with power law distributed sizes. Example structures of this type are shown in Figure 1, details of the construction process are discussed in Ref. [7].

From the point of view of graph theory each elementary segment in Figure 1a is an edge of the graph. Each edge is delimited by two nodes. In the simple case in Figure 1 nodes are laid out as in a 2D square lattice, whereas more complex and higher dimensional arrangements are possible. We call $N$ the total number of nodes, $N_{x}$ ($N_{y}$) the number of nodes in the horizontal (vertical) direction, and $E$ the number of edges. The HFN in Figure 1 is periodically continued in the horizontal direction, and its linear size in both the $x$ and the $y$ direction is $L=N_{x}$ (in adimensional units). As for the vertical direction, the top and bottom rows act as boundaries. According to the specific construction implemented in Figure 1, $N_{y}=N_{x}+1$ and the total number of non-boundary nodes is $N^{\prime}=N-2N_{x}=N_{x}(N_{x}-1)$.

To simulate damage in these structures, we follow the standard quasi-static simulation method for the random fuse model (RFM). We consider each edge as a load carrying element of unit conductivity. Voltages $V_{i}$ can be measured at each node $i$ and result in a current $I_{ij}=V_{i}-V_{j}$ carried by edge $ij$. Boundary conditions are applied as fixed voltages at the top boundary nodes, while each link $ij$ is initially assigned a fixed random current threshold $t_{ij}$. In each iteration of the simulation, the boundary voltage is set to $V=1$, and equilibrium voltages at each non-boundary node are calculated solving numerically the Kirchhoff nodes’ law equations (algebraic sum of incident currents is zero) for each node. Currents are computed accordingly. The edge carrying the highest $r=I_{ij}/t_{ij}$ will be the first one to break as soon as the external $V$ is set to $1/r$. Consequently, the weak edge is removed, the global $I$ and $V$ (correctly rescaled by $1/r$) are stored, and the simulation moves to the next iteration, until formation of a system wide crack sets the global conductivity to zero. While results generally depend on the choice of the probability distribution function $P(t_{ij})$, it was shown that the qualitative features of the damage accumulation process, viz., the emergence of super-rough cracks and the observation of generic scale invariance in the precursor statistics, are not affected by changing $P(t_{ij})$ from a uniform distribution (widely used in the statistical physics literature) to a Weibull distribution (more realistic, according to reasoning in terms of materials strength) [7]. The numerical results shown in the following consider the case of a uniform threshold distribution, for the deterministic HFN introduced above, but the findings carry over to Weibull distributed thresholds.

3 Spectral theory and localization

In this section we introduce the spectral graph analysis for fuse networks and highlight its main results for the HFN models considered. To this end, we will introduce two different matrices: i) the Laplacian matrix $\mathrm{L}$, which appears in several graph theoretical studies; and ii) the Dirichlet Laplacian $\mathcal{L}$, which arises from $\mathrm{L}$ when applying boundary conditions and provides the correct description of our RFM equations. We recall that the Laplacian nature of the RFM equilibrium equations can be seen by writing the generic Kirchhoff nodes’ law for each node $i$ as

[TABLE]

where $\mathrm{A}_{ij}$ is the generic element of the symmetric $N\times N$ adjacency matrix $\mathrm{A}$: by definition $\mathrm{A}_{ij}=1$ if $i$ is neighbour of $j$, or zero otherwise. Eq. 1 can be rewritten as

[TABLE]

where $\mathbf{V}$ is the column vector of generic element $V_{i}$ and $\mathrm{L}$ is the Laplacian matrix of the network (more precisely, the symmetric non-normalized graph Laplacian), with $\mathrm{L}_{ij}=\delta_{ij}\sum_{l}\mathrm{A}_{jl}-\mathrm{A}_{ij}$ . Eq. 2 simply states that the equilibrium equation for our scalar elasticity is a discretized Laplace equation of the type $-\nabla^{2}\mathbf{V}=0$. By construction, if the network consists of a unique connected component, the Perron-Frobenius theorem ensures that $\mathrm{L}$ has eigenvalues $0=\mu_{1}<\mu_{2}\leq\mu_{3}...$. Being $\mu_{1}=0$, Eq. 2 has infinitely many solutions, arising from the (unique) zero eigenvalue $\mu_{1}$. This property simply reflects the fact that infinitely many equilibrium states are possible when no boundary conditions are applied. By applying boundary conditions (i.e. by fixing the values of $V_{i}$ for $1\leq i\leq N_{x}$ and $N-N_{x}+1\leq i\leq N$), Eq. 2 is reduced to a $N^{\prime}-$dimensional system of equations of the form

[TABLE]

where we recall that $N^{\prime}=N-2N_{x}$ is the number of non-boundary nodes. Here $\mathrm{\mathcal{L}}$ is a Dirichlet Laplacian, with lowest eigenvalue $\lambda_{1}>0$, which ensures a unique solution to the electric/elastic equilibrium problem. Here $\mathbf{v}$ is a vector whose $v_{i}$ element is the voltage at non-boundary node $i$, while $\mathbf{b}$ is fixed and depends on the voltages at the boundary nodes.

While Eq. 3 represents the normal form in which the equilibrium problem is fed to a numerical solver in order to compute $\mathbf{v}$, here we calculate its analytical solution, which can be formally expressed in terms of the eigenvalues $\lambda_{1},\lambda_{2},...$ and corresponding normalized eigenvectors $\mathbf{w}_{1},\mathbf{w}_{2},...$ of $\mathrm{\mathcal{L}}$ as

[TABLE]

where $\cdot$ indicates the standard scalar product. It can be seen from Eq. 4 that the equilibrium configuration can be expressed as a linear combination of eigenvectors of $\mathrm{\mathcal{L}}$, in which eigenvectors corresponding to near-zero eigenvalues (the lower spectral edge) have a high weight because of the $1/\lambda_{m}$ factor. While these considerations may be of little relevance for regular lattices, which do not exhibit exotic spectral properties, they become crucial in hierarchical networks like the ones considered here. In particular, hierarchical modular networks: (i) possess a number of exceedingly small eigenvalues; (ii) these eigenvalues are associated with localized eigenvectors [15]. Eigenvector localization [18, 19] is the property by which the components of an eigenvector are all vanishing except for a subset (see Fig. 2 for an eigenvector of a HFN). Property (ii) tells us that the building blocks of the final $\mathbf{v}$ are localized patches in the networks, while property (i) ensures that only a few of those patches contribute to $\mathbf{v}$, and thus that voltages/displacements and currents/loads are localized in specific regions of the system. Localization of loads implies that, rather than growing a critical crack, the system develops localized patterns of damage, which eventually coalesce into a super-rough crack as the one shown in Fig. 1.

In order to quantify eigenvector localization, it is customary to consider the inverse participation ratio (IPR) for the $m-$th eigenvector as

[TABLE]

with higher values of $I_{m}$ pointing to stronger localization [18]. Fig. 3 shows values of the IPR for eigenvectors corresponding to the $20$ lowest eigenvalues, for a HFN and for a reference square lattice, at three different stages of the loading curve. While our results confirm the intuition of an increase of localization near the peak load, such an increase is much more significant in the hierarchical case, where it nears a one-order-of-magnitude jump in IPR at $I_{p}$.

A clearer picture of how effectively localization signals the vicinity to $I_{p}$ in HFN is provided in Figure 4 for much larger system sizes. The inverse participation ratio determined for an eigenvector in the lower spectral edge (in this case, the eigenvector corresponding to the lowest eigenvalue) starts to increase exactly at the peak current, providing a clear-cut indicator of the system’s incipient failure. Figure 4 (left) also provides important information regarding the size dependence of damage induced localization, as we detail in the following. We recall that, for large system sizes, $I_{m}\sim 1/N$ for delocalized states; any slower scaling is considered a signal of localization. Given the rescaling on the vertical axis of Figure 4, we can conclude that damage induced localization is a robust phenomenon in HFN. The robustness of our results is further confirmed in Figure 4 (right), where the above results are averaged over different network realizations, showing little variation with respect to the single realization in Figure 4 (left).

4 Hausdorff dimension

We have seen above how a sharp increase in eigenvector localization signals failure in HFNs, in a more reliable way than standard methods based on avalanche statistics. More importantly, this type of indicator is a structural one, that is one that can be inferred from the morphology of a sample, rather than from its response to load. The network dimension is in our view another candidate damage indicator, which may exhibit significant changes as failure advances. The dimension of a network is a generalization of the concept of lattice dimension. As various definitions exist, we focus here on the intrinsic Hausdorff dimension $d_{t}$, or topological dimension, which can be computed as follows: i) for each node $i$ in the network, compute the number $n_{ir}$ of nodes that are $r$ steps away from $i$, for all possible $r$; ii) take the node average $\bar{n}_{r}$; iii) if $\bar{n}_{r}\sim r^{d_{t}}$, $d_{t}$ is the topological dimension of the network. It is apparent that in regular lattices with nearest neighbor connections $d_{t}$ equals the lattice dimension, up to finite size effects (the rigorous definition of $d_{t}$ implies an infinite size limit). It is also clear that any damage process on any network structure can only decrease its value of $d_{t}$ (up to fluctuations due to the limited accuracy of exponent measurements). Figure 5 shows the evolution of $d_{t}$ for increasing damage in a HFN, throughout the entire loading curve. The peak current is accompanied by a sharp decrease in $d_{t}$. For comparison, we also show the evolution of $s_{1}$, the size of the largest connected component, which is normally monitored in percolation theory. As expected from the above considerations, $s_{1}$ is only mildly related to the softening of the system, and in particular it overestimates the peak current. The topological dimension $d_{t}$, instead, matches the value of $I_{p}$ with greater precision.

5 Discussion and conclusions

While network-related methods are usually of limited relevance in the study of materials, they become more relevant in the case of materials with complex microstructres, which can be represented by complex contact patterns. This appears to be the norm in biological materials, but also may be relevant in architectured and porous materials, as well as in metamaterials. Such systems can be naturally described as complex network structures. The understanding of the mechanical properties of such materials requires the introduction of quantities whose evolution under damage can be monitored, providing meaningful indicators of failure. In this work we focused on the case of materials with hierarchical microstructures and their modeling as hierarchical fuse networks. Damage in these systems has a peculiarity: it is poorly described by theories that rely on phase transitions and critical points, since scale invariant behavior is observed at any loading stage. Here we focused instead on structural features of the networks at hand, and highlighted their relationships with the damage state of a system under load. We find that the peak current in a HFN is accompanied by a sharp increase in eigenvector localization and a decrease in dimensionality. While both tendencies are easily to understand from an intuitive viewpoint, we find that their onsets mark precisely the peak current values in our simulations.

Our results provide general guidelines regarding what spectral and topological properties materials should possess in order to recover the behavior of HFN. A question of more practical relevance would be how these properties can be measured in experiments. The gradual shift of low eigenvalues to lower values for increasing damage is probably the easiest to measure, as each eigenvalue $\lambda_{i}$ of the lattice Laplacian is proportional to the square of the $i-$th lowest oscillation frequency, $\omega_{i}^{2}$. It is also easy to show that in the case of over-damped relaxation, $\lambda_{i}$ is inversely proportional to the relaxation time of the $i-$th longest relaxation time $\lambda_{i}\sim 1/\tau_{i}$. These quantities are amenable to direct measurement e.g. in terms of acoustic spectral properties of a component or structure. As for the possibility to measure eigenvectors, Eq. 4 states that the voltage field is approximately a linear combination of the eigenvectors in the lower spectral edge: the localization properties of individual eigenvectors should carry over to voltages (in terms of mechanical analogue: local strains) and, by extension, to any property whose spatial structure is governed by the Laplacian operator. To visualize this, Figure 6 shows results from a proof-of-concept study, in which we simulate current-controlled load in a HFN and we quantify the localization of the strain profile at each loading stage. Localization is measured as the inverse participation ratio of the vector $\mathbf{v}$. As the system reaches the peak current and softens, strain localization increases substantially, and this localization may be experimentally accessible e.g. from an analysis of digital image correlation data.

Beyond linear properties, localized features also emerge within inherently non-linear processes associated with patterning and damage in multi-scale systems, with examples ranging from activity patterns in brain networks [20, 21] and pattern retrieval in artifical neural networks [22], to shear band formation in plastically deforming amorphous materials [23], or detachment and frictional sliding of bio-mimetic hierarchically patterned materials [24]. A generic approach to identify damage-induced changes in connectivity patterns associated with some observable $A$ might be to discretize measured spatial patterns of $A$ and compute the associated correlation matrix $C_{ij}$ which may then be appropriately sparsified by applying a physically motivated threshold value. The resulting sparse matrix can be interpreted as the adjacency matrix of a correlation network (a functional network, in the neuroscience jargon), whose spectral properties and topological dimension depend on the state of damage and can be measured as presented in this paper. We hope that the main concepts discussed here can stimulate further interest in the usage of this class of methods to measure damage in broader classes of materials and related to a wider range of structural or functional materials properties.

Acknowledgements

The authors acknowledge support from DFG under grants ZA 171/9-1, MO 3049/1-1 and 3049/3-1, as well as by the European commission under H2020-MSCA-RISE project 734485 FRAMED.