Creep failure of hierarchical materials

TL;DR

This paper studies how hierarchical structures in materials affect their long-term strength and failure under constant stress.

Contribution

The study introduces a new failure mode in hierarchical materials that differs from traditional crack-based failure.

Findings

Hierarchical patterning increases creep strain and lifetime while reducing variation.

The new failure mode aligns with experimental results on patterned paper sheets.

Hierarchical structures induce different failure mechanisms compared to random cross-linking.

Abstract

Creep failure of hierarchical materials is investigated by simulation of beam network models. Such models are idealizations of hierarchical fibrous materials where bundles of load-carrying fibers are held together by multi-level (hierarchical) cross-links. Failure of individual beams is assumed to be governed by stress-assisted thermal activation over local barriers, and beam stresses are computed by solving the global balance equations of linear and angular momentum across the network. Disorder is mimicked by a statistical distribution of barrier heights. Both initially intact samples and samples containing side notches of various length are considered. Samples with hierarchical cross-link patterns are simulated alongside reference samples where cross-links are placed randomly without hierarchical organization. The results demonstrate that hierarchical patterning may strongly increase creep strain and creep lifetime while reducing the lifetime variation. This is due to the fact that hierarchical patterning induces a failure mode that differs significantly from the standard scenario of failure by nucleation and growth of a critical crack. Characterization of this failure mode demonstrates good agreement between the present simulations and experimental findings on hierarchically patterned paper sheets.

Figures

Top: iterative “top-down” construction of a deterministic hierarchical beam network: a beam network is divided by load-parallel cuts (removal of vertically adjacent cross links) into four highest-level modules. These form two groups of two modules loaded in series, connected by a system spanning connector (the row of cross links shown in red); next, each of the four modules is again divided into four lower-level modules plus a module-spanning lateral connector, etc.; Bottom: Network variants as discussed in the text: Deterministic Hierarchical Beam Network (DHBN), Stochastic Hierarchical Beam

Creep strain versus normalized time curves for structures of size $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=512$$\end{document}$ ; left: deterministic hierarchical beam network (DHBN): (a) un-notched, (c) notch length = 0.2L, (e) notch length = 0.4L; right: random beam network (RBN): (b) un-notched, (d) notch length = 0.2L, (f) notch length = 0.4L; each plot shows results of 20 simulated samples, the blue

Cumulative distribution of creep failure time for different structures of size $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=512$$\end{document}$ , deterministic hierarchical beam network (DHBN) and random beam network (RBN): (a) un-notched, (b) notch length = 0.2L, (c) notch length = 0.4L; The constant applied stress on both DHBN and RBN is $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasys

Axial strain $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _{yy}$$\end{document}$ map, left: deterministic hierarchical beam network (DHBN), right: random beam network (RBN) structures of size $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt}

Avalanche size distributions (probability vs. number of broken beams) during the different stages of the creep curve.

Full text

Introduction

Hierarchical materials have microstructures whose geometrical features repeat on different length scales in a self-similar fashion. Biological materials provide compelling examples: Materials such as nacre, bivalve shell, bone, and tendon collagen exhibit multi-level hierarchical microstructures which extend from the nanoscale to the macroscale. Such microstructures serve a dual purpose: By efficiently bridging from the nanoscale, on which the smallest hierarchical features are found, to the macroscale, they directly harness nanoscale size effects, which endow the smallest features with increased strength as compared to macroscopic structures consisting of the same material^1^. At the same time, it has been suggested that hierarchical organization may delay or prevent the spreading of flaws driven by local stress concentrations e.g. at crack tips, which often controls failure of non-hierarchical materials. This may happen because, as shown by^2^, hierarchical architecture may extend the fracture process zone in such a manner that it ultimately reaches the specimen dimensions and stress/strain concentrations are avoided. As a consequence, hierarchical materials may exhibit remarkable stiffness, good toughness, and high strength, even when they are made of brittle constituents^3–8^.

Numerical modeling has been envisaged as a tool for designing hierarchical structures and predicting their properties. Simple modelling approaches such as fiber bundle models have been extended to hierarchical materials^9^ but cannot adequately represent crack-tip stress concentrations which are an essential aspect in the failure of materials by crack propagation. Spatial stress concentrations and crack propagation can be described by fuse and beam network models. As an example, Moretti et al.^10^ used a random fuse model (RFM) to comparatively study failure behaviour of materials with and without hierarchical microstructure. The failure behaviour of hierarchical materials, as well as the resulting super-rough crack morphology, was found to differ significantly from the behavior of non-hierarchical ones. Hosseini et al.^11^ introduced a hierarchical version of a beam network model (BNM) to simulate hierarchically patterned materials. The findings of these authors show that the failure of hierarchically structured materials is caused by local damage nucleation followed by damage percolation rather than critical crack propagation, a finding which holds for both two-dimensional^11^ and three-dimensional^12^ structures. In both cases, crack profiles show large deflections as failure arises from the coalescence of widely separated micro-cracks. Zaiser et al.^13^ conducted both simulations and experiments on two-dimensional quasi-brittle sheets loaded along a single axis. It was demonstrated that hierarchical patterning considerably increased crack propagation resistance both in terms of the peak stress and the work of failure of pre-cracked samples. The differences in failure mode between hierarchical and non-hierarchically patterned 2D material were further investigated by Pournajar et al.^14^, who studied failure precursors in patterned paper sheets under both displacement-controlled and creep loading. Failure precursors were characterized in terms of acoustic emission and digital image correlation to map the evolution of strain patterns in the run-up to failure. The results confirm the idea that hierarchical patterning efficiently mitigates against stress/strain concentrations and promotes a diffuse mode of failure that is insensitive to the existence of even large flaws.

The hierarchical structures considered in the above mentioned works can be alternatively envisaged as assemblies of laterally connected load carrying fibers which are bundled into hierarchically nested “modules”, or as solid sheets (in 3D: blocks) of a matrix material that are divided into modules by load-parallel gaps and that exhibit a power-law distribution of gap lengths (in 3D: gap areas). From the perspective of composite mechanics, such structures can be considered as limit cases of “calcite-like” hierarchical composites as studied by Buehler and co-workers^8,15^ where weak and compliant, lamellar inclusions are embedded into a hard and brittle matrix in a hierarchical manner: In our work the inclusions have zero elastic modulus and near-zero volume fraction as they correspond to the gaps, whereas the matrix corresponds to the matrix material.

Most simulation studies of failure of hierarchical materials have focused on load-driven failure, notably under uni-axial loading conditions, where a load or displacement is imposed on the system and increased until the system loses connectivity. Typical of such scenarios, the process of damage accumulation and ultimate failure does not depend on loading rate, i.e. the material behavior is time independent. In real materials, however, an important failure pathway is that of creep, or sub-critical failure, where the system accumulates damage and approaches failure over an extended period of time. Here we thus investigate creep failure of hierarchical materials, i.e. the time evolution of damage and fracture under conditions of constant sub-critical load.

Methods

Construction and geometrical properties of the studied beam networks

The construction method which we use has been described in detail elsewhere^10,11^. Here, we only briefly summarize the main features. In all cases, we start out from a simple square lattice of beams where the loading axis is aligned with one of the cubic lattice directions. The size of the lattice is L which we take to be a power of 2, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=2^n$$\end{document} , where n will define the number of hierarchical levels. There are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} load-carrying beams in loading (henceforth “vertical”) direction forming L load-carrying columns, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} beams that are initially unloaded (“cross links”) in lateral (“horizontal”) direction forming L rows. The points at which the beams are connected are denoted as nodes; the endpoints of the load carrying columns are denoted as upper and lower boundary nodes. In all cases, loading is carried out by keeping the bottom boundary nodes fixed and displacing the top boundary nodes rigidly upwards such as to create a prescribed total load. Periodic boundary conditions are applied in lateral direction. From this basic set-up, different variants are constructed by removal of beams while the node locations are always kept the same.

To construct deterministic hierarchical beam networks (DHBN), a hierarchical structure is obtained by removing cross links such as to create gaps which recursively sub-divide the structure into load-carrying modules of decreasing order as illustrated in Fig. 1, top.Figure 1. Top: iterative “top-down” construction of a deterministic hierarchical beam network: a beam network is divided by load-parallel cuts (removal of vertically adjacent cross links) into four highest-level modules. These form two groups of two modules loaded in series, connected by a system spanning connector (the row of cross links shown in red); next, each of the four modules is again divided into four lower-level modules plus a module-spanning lateral connector, etc.; Bottom: Network variants as discussed in the text: Deterministic Hierarchical Beam Network (DHBN), Stochastic Hierarchical Beam Network (SHBN) and Random Beam Network (RBN).

This construction leads to a hierarchical pattern of modules separated by gaps, and of connecting cross links. In addition to DHBN, we could consider a randomized variant, so called stochastic hierarchical beam networks, SHBN, constructed by first randomly reshuffling whole columns of cross-links and gaps, and then the resulting rows. Previous studies have addressed the comparison of DHBN and SHBN structures, and have found that both variants usually exhibit comparable fracture properties^10,11^.For this reason, here we limit our analysis to the deterministic hierarchical variant DHBN.

Furthermore, as non hierarchical reference structures, random beam networks (RBN) are constructed by removing randomly the same number of cross links as in a DHBN or SHBN structure. The different network variants are illustrated in Fig. 1, bottom.

In all cases, we may introduce cracks of integer length a into the structure by removing, at a randomly chosen location, a row of a adjacent load carrying beams.

Stress re-distribution

Displacement of the boundary nodes induces forces and moments in the beam network. To evaluate these forces we envisage straight, identical beams of unit modulus of elasticity, unit length and square unit cross section which are capable of transmitting axial and shear forces as well as bending moments. The beams are rigidly clamped to the nodes, and each node has two translational degrees of freedom (node displacements u and v along the global x and y axes) and one rotational degree of freedom (rotation angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} about the global z axis). The beams are assumed to deform in a linearly elastic manner, their deformation is described using Timoshenko beam theory. Evaluating the resulting balance equations for linear and angular momentum at all nodes in the quasi-static limit and using a small-strain approximation results in a linear system of equations for the nodal degrees of freedom, whose solution gives us the local forces \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_i$$\end{document} and bending moments \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_i$$\end{document} acting on all beams (for details see Hosseini et al.^11^).

Failure criterion

Beams deform in a linearly elastic manner and undergo time-dependent failure where the failure time is controlled by the maximum deviatoric stress magnitude acting on the beam. We define the maximum deviatioric stress acting on the beam connecting nodes i and j as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sigma _{ij} = \sqrt{ \Big ( \dfrac{F_{i}n_i}{A}+\textrm{max}(|M_{i}|,|M_{j}|)\dfrac{y_{\textrm{max}}}{I} \Big ) ^ {2}+3 \Big ( \dfrac{Q_{i}}{A} \Big ) ^ {2} } \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i$$\end{document} indicates the outward normal direction of the beam end surface connecting to node i, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_i$$\end{document} is the normal force which can be tensile ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_i n_i > 0$$\end{document} ) or compressive ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_i n_i < 0$$\end{document} ), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q_i$$\end{document} is the shear force, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_i$$\end{document} is the moment acting on this surface. A derivation of Eq. (1) based on the Maximum Distortion Energy of Failure (von Mises) criterion was given by Hosseini et al.^11^. From the deviatoric stress, the beam failure rate is calculated under the assumption that it is controlled by a thermally activated process:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nu _{ij} = \nu _0 \exp \left[ -\frac{H_{ij} - \sigma _{ij}V_{\textrm{a}}}{k_{\textrm{B}} T}\right] = \nu _0 \exp \left[ -\frac{t_{ij} - \sigma _{ij}}{\Theta }\right] \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _0$$\end{document} is a characteristic attempt frequency, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_{ij}$$\end{document} is the activation barrier for beam failure, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{\textrm{ij}}=H_{ij}/V_{\textrm{a}}$$\end{document} is the corresponding beam failure stress, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{\textrm{a}}$$\end{document} is an activation volume, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta = k_{\textrm{B}} T/V_{\textrm{a}}$$\end{document} is the characteristic thermal energy density therein. Equation (2) holds as long as the beam stress is below the beam failure stress, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{ij} < t_{ij}$$\end{document} (subcritical failure), if on the other hand \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{ij} \ge t_{ij}$$\end{document} , then failure occurs instantaneously.

Mimicking material heterogeneity, beam failure stresses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{ij}$$\end{document} are randomly assigned based on a Weibull probability distribution function with probability density

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(t_{ij}) = \dfrac{\beta }{\eta } \Big ( \dfrac{t_{ij}}{\eta } \Big ) ^ {\beta - 1} \exp \left( -\Big ( \dfrac{t_{ij}}{\eta } \Big )^\beta \right) \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta > 0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta > 0$$\end{document} are the shape and scale parameters of the distribution, respectively. In this paper we normalize all stresses by the mean failure threshold stress and adjust the shape parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} such as to have a fixed mean value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle t_{ij} \rangle = 1.0$$\end{document} . The shape parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} may be varied to implement different degrees of disorder.

Simulation protocol

Simulations are performed to mimic creep conditions, i.e., deformation under a stationary temporally constant global load \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\textrm{c}}$$\end{document} . However, as in experiment, imposing this load requires an initial, non stationary loading stage.

Initial loading

The simulation starts with an initial loading stage where the external axial displacement is increased until the global stress level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} (evaluated as the sum of axial beam forces on the top nodes, divided by L) reaches the prescribed level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma = \sigma _{\textrm{c}}$$\end{document} or a beam breaks because \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{ij} \le \sigma _{ij}$$\end{document} . If beam failure occurs during initial loading, the displacement is then kept fixed while load is re-distributed which may trigger further instantaneous failures. Load re-calculation is repeated iteratively in such a way that one beam is removed in every iteration, until all beams are in a subcritical state ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{ij} > t_{ij})$$\end{document} ; this sequence of failures in general leads to a stress relaxation. The axial displacement is then again increased until either the prescribed creep stress is reached or beams fail, and this is repeated until the system reaches a sub-critical state at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma = \sigma _{\textrm{c}}$$\end{document} or the system becomes disconnected; in the latter case, we terminate the simulation and note instantaneous system failure during initial loading.

Creep behavior

Once the system has settled into a sub-critical state, we use a Kinetic Monte Carlo Algorithm to determine the location of the next beam failure. We count the number of thermally activated failures by the discrete index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_T$$\end{document} . To identify the failing beam, we evaluate the failure rates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _{ij}$$\end{document} for all remaining intact beams and choose one of them with probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{ij} = \nu _{ij}/\nu _{\textrm{tot}}$$\end{document} where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _{\textrm{tot}} = \sum _{ij} \nu _{ij}$$\end{document} . This beam is removed from the system, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_T$$\end{document} is increased by on, and time is increased by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t(n_T)$$\end{document} which we take to be an exponentially distributed random variable of average \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _{\textrm{tot}}^{-1}$$\end{document} .

Subsequent to a thermally activated beam failure, the stresses on all beams are re-computed, and beams where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{ij}> t_{ij}$$\end{document} are removed instantaneously. The external displacement is then adjusted to maintain the global stress at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma = \sigma _{\textrm{c}}$$\end{document} , and the process is repeated until either all remaining beams are in a sub-critical state or the system disconnects. The number of beams that fail instantaneously after the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_T$$\end{document} th thermally activated failure is denoted the avalanche size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s(n_T)$$\end{document} , the displacement after the avalanche is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d(n_T) = L \epsilon (N_T)$$\end{document} where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon (N_T)$$\end{document} is the total creep strain, and the avalanche time is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sum _{n_T}\Delta t(n_T)$$\end{document} .

After an avalanche has terminated, we again use a Kinetic Monte Carlo step to determine the next thermally activated beam failure, and this sequence is repeated until the system globally disconnects. The corresponding time is called the creep lifetime \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{\textrm{c}}$$\end{document} .

Results

Creep curves and time to failure

Figure 2 shows typical creep curves (creep strain vs time) obtained for DHBN and RBN samples with different initial crack lengths. In all simulations, temperature was set so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta = 0.01$$\end{document} and the Weibull parameter of the barrier height distribution was \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta = 4$$\end{document} . In the plots, the time axis is normalized by the failure time, which allows to average curves pertaining to different samples to obtain average creep curves (blue lines in Fig. 2).Figure 2. Creep strain versus normalized time curves for structures of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=512$$\end{document} ; left: deterministic hierarchical beam network (DHBN): (a) un-notched, (c) notch length = 0.2L, (e) notch length = 0.4L; right: random beam network (RBN): (b) un-notched, (d) notch length = 0.2L, (f) notch length = 0.4L; each plot shows results of 20 simulated samples, the blue lines represent the average creep strain curves; for each crack length, the applied stress amounts to 80% of the mean RBN failure stress at the respective crack length: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\sigma /\langle \sigma _{\textrm{f}}\rangle )$$\end{document} = 0.04, 0.06, 0.16 for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 0.4L, 0.2L, 0$$\end{document} respectively; for other parameters, see text.

With exception of the uncracked RBN, all samples show typical three-stage creep curves. During initial loading, an instantaneous strain builds up which is almost identical for hierarchical DHBN and non-hierarchical RBN samples. This is followed by a regime of decreasing creep rate (creep stage I), then an extended regime of almost constant slope (creep stage II) and an accelerating regime in the run-up to failure (creep stage III).

Stage II shows clear differences between DHBN and RBN samples. First of all, the overall strain that accumulates in stage II is significantly higher in the hierarchical samples. Second, the deformation process of hierarchical samples is punctuated by large avalanches, where the rate of deformation suddenly jumps upwards and then relaxes back to a low level. In single samples, one might speak of multi-stage creep, however, the random timing of the strain bursts leads on average to a smooth average response (blue lines in Fig. 2, left) which is still well described as stage-II creep with a well defined average creep rate. The observation of large strain bursts in the stage-II creep curves of hierarchical samples in our simulations matches experimental observations of Pournajar et al.^14^, who reported such bursts in hierarchically patterned paper samples but not in their non-hierarchical counterparts.

Stage III is less pronounced in RBN samples where the onset of failure is more sudden than in DHBN samples. In uncracked RBN samples this stage is almost absent. These observations match the findings for the same sample types in displacement-controlled uni-axial loading: In RBN samples, the pattern of failure is essentially brittle – loading is almost elastic up to the critical failure stress, and after the peak load one observes an almost vertical load drop with little or no post-peak activity. These findings match the observation in the creep simulations of a small stage-II creep strain before the onset of softening and almost no strain accumulation after the onset of softening. In DHBN samples, on the other hand, displacement controlled tests are characterized by significantly higher strains-to-failure which accumulate both during loading, where damage accumulates in a series of avalanche-like bursts, and in the form of significant post-peak activity, leading to a much more gradual mode of failure as compared to the RBN samples.

Regarding time to failure, a complex picture emerges. For notched samples, the lifetime of DHBN samples significantly exceeds that of notched RBN samples subjected to the same creep load (Fig. 3b,c). At the same time, the lifetime scatter is very significantly reduced. Thus, in samples containing macroscopic flaws, hierarchical patterning not only provides enhanced flaw tolerance but also ensures a more predictable material response. The picture is different in case of samples without initial macroscopic notch (Fig. 3a). In this case, the logarithmic scatter of the lifetimes is comparable for DHBN and RBN samples, while the absolute creep lifetimes of the RBN samples are higher than those of the DHBN samples.Figure 3. Cumulative distribution of creep failure time for different structures of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=512$$\end{document} , deterministic hierarchical beam network (DHBN) and random beam network (RBN): (a) un-notched, (b) notch length = 0.2L, (c) notch length = 0.4L; The constant applied stress on both DHBN and RBN is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$80\%$$\end{document} of the peak stresses of RBN for each crack length, which corresponds to creep stresses in units of the mean beam failure stress \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \sigma _{B}\rangle $$\end{document} of 0.04, 0.06 and 0.16 for systems with crack lengths \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a= 0.4L, a= 0.2L$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=0$$\end{document} , respectively. Each plot represents the results of 20 simulations. For other parameters, see text.

Spatial and temporal patterns of damage accumulation

In order to analyze the reasons for the characteristic differences between the behavior of hierarchical and non-hierarchical samples, we look at the spatial and temporal patterns of damage accumulation that emerge in the run-up to creep failure. Figure 4 shows typical strain patterns in creep stage II in pre-cracked samples (crack length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=0.2 L$$\end{document} ), at about half the sample lifetime. For interpretation of the patterns it is important to note that, in a homogeneous sample and in absence of creep damage, the imposed creep stress of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\textrm{c}} = 0.06$$\end{document} would induce a corresponding homogeneous elastic strain of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _{yy} = 0.06$$\end{document} . Because of the crack in the sample center, the stress and strain above and below the crack are significantly reduced (shielding) while the strain ahead of the crack tips is enhanced. At the same time, we observe very significant differences between DHBN and RBN samples. In the RBN samples, damage accumulation is confined to two process zones right ahead of the crack tips, whereas in the rest of the sample, deformation remains almost elastic. Damage accumulation at the crack tip leads to a gradual crack advance, and the transition to failure is related to the pre-existing crack becoming supercritical. In DHBN samples, on the other hand, we observe a much higher degree of global damage accumulation. This damage is not localized at the tips of the pre-existing crack but distributed over the entire sample with exception of the shielded regions above and below the crack. The diffuse nature of damage accumulation without supercritical crack propagation allows the hierarchical material to sustain a much higher creep strain. Again these findings are in excellent agreement with the experimental observations reported by Pournajar et al.^14^ who used digital image correlation to map the strain distribution in paper samples with both hierarchical and non hierarchical cut patterns.Figure 4. Axial strain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _{yy}$$\end{document} map, left: deterministic hierarchical beam network (DHBN), right: random beam network (RBN) structures of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=512$$\end{document} -during creep steady state step at 0.5 failure time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{c}$$\end{document} in a creep simulation; The constant applied stress on both DHBN and RBN is 80% of the peak stress of RBN which is 0.06 in unit of mean beam failure stress \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \sigma _{B}\rangle $$\end{document} for systems having 0.2L precrack length. The precrack is marked in gray, for other parameters, see text.

We next look at the statistics of damage bursts or avalanches occurring during creep deformation. Failure of one beam leads to load being re-distributed onto other beams, and this can give rise to secondary failures, and ultimately to a cascade of beam failures (“avalanches”). While most thermally activated beam failures trigger only very few sequels, the largest avalanches may encompass hundreds of beam failures and are directly visible in form of sudden large strain increments on the creep curves (see Fig. 2). In the following, we define the avalanche size S as the number of instantaneous, “athermal” beam failures after any given thermally activated beam failure.

Avalanches emerging in failure of disordered materials often exhibit power-law statistics (for an overview, see^16^). This is also true for hierarchical materials subject to load-controlled deformation, as found in simulation studies^10,11^ as well as in experiments^14^.

We determine probability distributions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(\Delta )$$\end{document} by logarithmic binning, approximating

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P(\Delta ) \approx \frac{N(S,\Delta _S)}{N_{\infty } \Delta _S} \quad ,\quad N(S,\Delta _S) = N(S_i|(S-\Delta _S/2)< S_i < (S+ \Delta _S/2)) \end{aligned}$$\end{document}

where the bin sizes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta _S$$\end{document} increase exponentially with S to ensure adequate statistics in the large-S tails of the distributions. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{\infty }$$\end{document} is the total number of avalanches in the sample. Avalanche statistics was determined from 1000 creep simulations conducted for the same parameters as used in Fig. 4. Distributions were determined separately for the first \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10\%$$\end{document} of the failure time (roughly corresponding to creep stage I, Fig. 5a), the second 80% (creep stage II, Fig. 5b), and the last 10% (creep stage III, Fig. 5c), as well as for the total of all avalanches (Fig. 5d).Figure 5. Avalanche size distributions (probability vs. number of broken beams) during the different stages of the creep curve.

A first conspicuous observation is that the avalanche statistics shows, within the margins of statistical error, no differences between hierarchical and non-hierarchical samples. This is also observed in experiment^14^, but apparently contradicts the direct observations showing large jumps on the DHBN but not the RBN creep strain vs time curves, which moreover match the observations on experimental creep samples. The reason is simply that the largest avalanches, while conspicuous on the creep curves where they account for a large fraction of the creep strain, carry insufficient statistical weight to shape the statistical distributions.

Contrary to the findings in load-controlled simulations^10–12^, the distributions are not well described as truncated power laws. Rather, they are composed of a rapid decay at small avalanche sizes below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S=5$$\end{document} , and a much slower decaying long tail in the form of a stretched exponential which extends, in creep stages I and II, to maximum avalanche sizes of about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S=100$$\end{document} . In creep stage III, the maximum avalanche sizes are larger and range up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S\approx 600$$\end{document} . The total distribution is practically identical to the distribution of stage-III avalanches (Fig. 5d), for the simple reason that most avalanches occur in creep stage III.

If more substantial differences in avalanche size distributions exist between hierarchical and non-hierarchical systems, they might become visible at larger sizes ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=1024$$\end{document} and above), where one expects to find a larger avalanche size range before the onset of the stretched exponential regime. Such system sizes have been explored in the past in simplified hierarchical fuse models^10^, but achieving them with the current beam elasticity implementation, coupled with a Kinetic Monte Carlo algorithm, would require a dedicated large-scale numerical study. Given that the present results are in good agreement with the available experimental data (which inevitably also suffer from size restrictions), one may ask whether the added insight obtained from such a study is worth the effort.

Discussion and conclusions

Our simulation study shows characteristic differences in the creep behavior of hierarchical beam network samples as compared to non-hierarchical reference samples. These differences corroborate, for the case of creep deformation, the idea that hierarchical patterning is a powerful means to enhance the flaw tolerance of materials, by enhancing the creep lifetime in presence of flaws while reducing the lifetime scatter. This beneficial effect arises from the strong strain de-localization induced by the hierarchical patterning, which reduces crack-tip stress concentrations and reduces statistical scatter by distributing damage accumulation evenly across the sample, thus mitigating against the effect of strength fluctuations near the crack tip which otherwise have a strong impact on sample lifetime. These findings are in good agreement with the experimental data reported by Pournajar et al.^14^. At the same time, it must be noted that hierarchical patterning is not always beneficial; in samples without initial notches the creep lifetime is reduced, rather than enhanced, in hierarchical as compared to non hierarchical reference samples.

In line with experimental studies, our simulations indicate that the statistics of avalanche precursors to creep failure is not strongly influenced by hierarchical patterning. In the context of experimental data reported by Pournajar et al.^14^, this could be understood by noting that the avalanche size distributions took the shape of truncated power laws \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\Delta ) \propto \Delta ^{-\delta }$$\end{document} with an avalanche exponent close to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta = 1.5$$\end{document} , typical of mean-field behavior. However, the present simulation results are inconsistent with mean-field behavior and further work is needed to understand the nature and statistics of avalanches in the present model.

The results of the present study were obtained under the assumptions of perfectly elastic-brittle behavior of the constituent beams. While such an assumption is motivated by brittle collagen-based bioceramics such as nacre and bone, other protein-based hierarchical materials exist, such as spider silk, where a much broader variety in constituent behavior is observed^17,18^. In these cases we expect the purely elastic response to be well described by our model, in particular in regard to the mitigation of stress concentrations at crack tips. However the large-scale fracture behavior would result from the interplay of elasticity and more complex constitutive laws, a topic that we are considering for future work.

Our work has been based on the application of a simple stochastic model to the predication and statistical analysis of creep failure in hierarchical materials. An alternative approach consists in the use of machine learning techniques. Machine learning can help to predict failure times of individual metamaterial samples based on their creep records^19^, but also for the purpose of further optimizing the design of ordered or disordered hierarchical metamaterials in view of enhancing their mechanical performance^20–23^. Further developing such approaches beyond the realm of mechanical properties towards the formulation of automated design strategies for multi-functional hierarchical materials remains an important task for future investigations.