Fracture of disordered and stochastic lattice materials

TL;DR

This paper develops a probabilistic framework to analyze how disorder and stochastic properties influence failure and fracture behavior in lattice mechanical metamaterials, validated through finite element analysis.

Contribution

It introduces a novel probabilistic approach that couples stochastic failure and geometric disorder to predict failure in disordered lattice materials.

Findings

Disorder and stochasticity reduce average failure loads.

Variance in failure loads increases with disorder.

Fracture length and volume size constraints are identified.

Abstract

The failure of mechanical metamaterials is a function of the interplay between the properties of the base material and the microstructural geometry. Stochastic failure properties of the base material and disordered microstructural geometries can contribute to variations in the global failure mechanics that are not captured in traditional analyses of ordered, deterministic architected materials. We present a probabilistic framework that couples stochastic material failure and geometric disorder to predict failure in lattice mechanical metamaterials. These predictions are verified through finite element analysis, which confirm that disorder and stochasticity affect both the mean and variance of the damage initiation load in a lattice, with average failure loads being generally reduced and variance increasing with higher levels of disorder and stochasticity. The fracto-cohesive length and representative volume element size are also predicted and constrain the minimum defect and lattice sizes, respectively, for failure to be considered a fracture process. The framework is extended to consider the fracture behavior of the lattice, the development of damage zones, and their impact on the steady-state fracture toughness.

Equations

$${\sigma }_a(r,\theta )=R_{1i}{R}_{1j}{\tilde{\sigma }}_{ij}(r,\theta ),\text{ with }{R}_{ij}=\left\{\begin{array}{ccc}\cos (\varphi )& & -\sin (\varphi )\\ \sin (\varphi )& & \cos (\varphi )\end{array}\right\},$$

$$\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\end{array}\right\}=\frac{K_I}{\sqrt{2\pi r}}\cos \left(\frac{\theta }{2}\right)\left\{\begin{array}{c}1-\sin \left(\frac{\theta }{2}\right)\sin \left(\frac{3\theta }{2}\right)\\ 1+\sin \left(\frac{\theta }{2}\right)\sin \left(\frac{3\theta }{2}\right)\\ \sin \left(\frac{\theta }{2}\right)\cos \left(\frac{3\theta }{2}\right)\end{array}\right\},$$

$${\sigma }_a=\frac{{\tilde{K}}_I}{2\sqrt{2\pi r}}\cos \left(\frac{\theta }{2}\right)\left[\cos (2\theta -2\varphi )-\cos (\theta -2\varphi )+2\right],$$

$$P_{\sigma }(\sigma |m)=1-e^{-\frac{1}{V_o}{\int }_V{\left(\frac{\sigma }{{\sigma }_o}\right)}^{m}dV},$$

$$P_{\sigma }(r,\theta |\varphi ,m)=1-e^{-{\left(\frac{{\sigma }_a(r,\theta |\varphi )}{{\sigma }_o}\right)}^m},$$

$$f(\varphi |b,\mu {)}_{{\varphi }_i}=\frac{1}{2\pi I_0(b)}{e}^{b\cos (\varphi -\mu )},$$

$$P_f(r,\theta ,\bar{\sigma }|m,{\sigma }_o,b)=1-\prod _{i=1}^6\left[1-{\int }_{-\pi }^{\pi }{P}_{\sigma }(r,\theta ,\varphi ,\bar{\sigma }||m,{\sigma }_o)p_{\varphi }(\varphi |b,{\mu }_i)d\varphi \right].$$

$$F_f(\bar{\sigma }|m,{\sigma }_o,b)\approx 1-\prod _{j=1}^s\left[1-P_f(r_j,{\theta }_j,\bar{\sigma }|m,{\sigma }_o,b)\right].$$

$$\lambda =\frac{G_c}{W_f}$$

$$\lambda ={\left(\frac{K_{IC}}{{\sigma }_f^{\ast }}\right)}^2.$$

$$\xi =\frac{1}{2\pi }{\left(-\ln (1-T)\right)}^{-\frac{2}{m}}{\left(\frac{K_{IC}}{{\sigma }_0}\right)}^2.$$

$$W_c=\frac{K_{IC}^2}{E^{\ast }}\left(\frac{1-{\nu }^{\ast }}{2\pi \xi }\right),$$

$$W_c=\frac{{\sigma }_f^2}{E^{\ast }}(1-{\nu }^{\ast }).$$

$$h=\left(\frac{2\pi }{1-{\nu }^{\ast }}\right)\xi .$$

$$\frac{G_c}{G_o}=\frac{{\mathbb{E}}_{{\sigma }^2}}{{\sigma }_o^2}\left(\frac{\Omega }{{\Omega }_o}\right),$$

Full text

††journal: TBD

\affiliation

[label1] organization=Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, city=Philadelphia, state=PA, country=USA \affiliation[label2] organization=Department of Mechanical Engineering and Production, Aarhus University, city=Aarhus, country=Denmark

1 Introduction

The fracture properties of many materials emerge from their structural features at varying length scales. Fundamentally, the localized nature of fracture leads to the nucleation and spread of damage throughout fine-scale microstructures before ultimately resulting in global failure of the material. These microstructures, whether naturally occurring in high toughness materials like bone [1, 2] or in low toughness materials such as blood clots [3, 4], play a critical role in how damage develops near a crack tip, which in turn impacts the material’s toughness. In traditional engineering materials, knowledge of material microstructure can be used [5, 6] to inform predictions of fracture toughness; however, these structures, whether regular or random, are often inherent to the material or controlled rudimentarily through manufacturing process controls. Conversely, in structural mechanical metamaterials, precisely engineered fine-scale structures are designed and leveraged to control and enhance failure properties [7, 8, 9, 10, 11, 12, 13, 14] by altering the stresses around cracks or defects. In particular, lattice structures with fine-scale features are often engineered in such a way that damage at the microscale governs the global failure behavior of the material [15, 16, 17, 18, 19, 20]. While the failure of simple structures, such as uniform lattices, can be precisely forecast and linked to broader failure characteristics [13, 14, 20, 21], predicting failure in more complex or disordered structures remains a significant challenge. Disordered materials have shown substantial enhancements in toughness, relative to equivalent ordered materials [22], which they achieve through the generation of distributed damage [23]. In these cases global failure analysis is typically performed by measuring the local damage in simulations or experiments [22, 24] or inferred using data-driven methods such as graph neural networks [25]. However, both of these approaches fail to elucidate the underlying mechanisms of microscale damage and its correlation to global fracture properties.

Additionally, the interplay between architecture and non-deterministic material properties must be considered to accurately predict the failure of structured materials. For brittle materials, stochasticity [26] often plays a significant role in the failure process. This has been shown to result in significant impacts for structured lattice materials, with highly stochastic lattices exhibiting large losses in average strength [27]. For fracture, it has also been shown [28] that stochastic failure has a significant effect on the fracture length scale [29], $\lambda$, also called the fracto-cohesive length [30] of the lattice material. Sufficiently stochastic materials may achieve critical length scales such that cracks that would otherwise cause a fracture-like failure no longer effect the failure loads of the material [28]. Thus, in order to fully leverage material architecture to control and enhance the fracture toughness of brittle structured materials, a model is required that integrates the effects of the architected material geometry and stochastic failure.

In this work, we present a framework for predicting failure of elastic-brittle materials featuring disordered and irregular lattice-like microstructures. The framework includes considerations for stochastic material failure and demonstrates the connection between geometric variability in the material’s structure, and stochasticity in its failure properties. Loads to cause the initiation of damage are predicted rigorously, along with estimates of the fracto-cohesive length and representative volume element size that control whether a fracture-like process occurs. The model is extended to predict crack-tip damage zones that may also be used to estimate the steady-state fracture toughness. Results are verified through finite element simulations of a triangular lattice, which reveal how toughness is enhanced through geometric disorder, while stochastic failure can enhance toughness of ordered lattices but generally reduces the toughness of disordered lattices.

2 Theoretical Framework

2.1 Stresses in a Stretch-Dominated Lattice

A lattice material can be assumed to have a pseudo-continuous global stress state, $\sigma_{ij}$, analogous to a homogeneous material, provided that the length scale of the lattice’s geometric features are much smaller than the specimen size and the characteristic length scale of the stress distribution. Under these conditions, for a 2-D structured material with a lattice microstructure that is stretch-dominated [31], the axial stress of any isolated ligament oriented at an angle $\phi$ can be calculated as the component of the stress field in that direction. Under the assumption of a stretch-dominated material, this is the only pertinent stress component. Thus, the ligament stress, $\sigma_{a}$, at a position $(r,\theta)$ is given by

[TABLE]

where $\phi(r,\theta)$ is the angle made by the ligament to the $x$-axis, and $\tilde{\sigma}_{ij}=\bar{\rho}\sigma_{ij}$ is the local stress state resulting from the relative density, $\bar{\rho}$. Given the linear scaling of the local stress with the relative density, it is not difficult to infer behaviors of lattices at varying densities so the model is presented for the general case of $\tilde{\sigma}_{ij}$.

As an example, we consider the case of mode I fracture of a linear elastic-brittle material as shown in Fig. 1(a). In this case, the global stress distribution is singular around a sharp crack tip [32]. Under plane stress conditions using polar coordinates centered at the crack tip, these stresses are given [29] by

[TABLE]

where $K_{I}=\bar{\sigma}\sqrt{\pi a}$ is the mode I stress intensity factor, with $\bar{\sigma}$ being the far-field applied stress and $a$ the crack length.

The characteristic length scale of the distribution was shown by Irwin [29] to be related to the size of the process zone, $\lambda_{0}$. If the length scale of the lattice is much smaller than the process zone size, $\lambda_{0}$, it can be assumed that the crack-tip stress field of a homogeneous specimen will be largely maintained at the length scale of the crack and the specimen. Predicting $\lambda_{0}$ is considered later in this work. Thus, if $\sigma_{ij}$ is given by eq. (2), the axial stress in a ligament becomes

[TABLE]

where $\tilde{K}_{I}=\tilde{\sigma}\sqrt{\pi a}$. Failure of any individual ligament occurs when its axial stress $\sigma_{a}=\sigma_{f}$, where $\sigma_{f}$ is the failure stress for a brittle material. This is in contrast to an elastic-brittle homogeneous material, where failure via fracture is determined when the stress intensity factor reaches a critical value, $K_{IC}$, the fracture toughness, which is a material property. This is related to the material’s toughness [29] (i.e., the critical energy release rate) by $G_{c}=K_{IC}^{2}/\bar{E}$, where $\bar{E}=E$ in plane stress and $E/(1-\nu^{2})$ in plane strain, where $E$ is the Young’s modulus of the material, and $\nu$ is the Poisson’s ratio.

In this work, we consider the example case of a triangular 2-D lattice, as shown in Fig. 1(b), although extensions to other stretch-dominated lattices are also discussed. In a regular triangular lattice with a microscale unit cell size, we can assume that at any location $(r,\theta)$, there will be six ligaments oriented at fixed angles, $\phi=-2\pi/3,-\pi/3,0,\pi/3,2\pi/3$, or $\pi$ (assuming no geometric disorder). Thus, ligament failure is deterministically given by the condition $\sigma_{a}(r,\theta|\phi)\geq\sigma_{f}$, for any of the six possible values of $\phi$.

2.2 Effect of Material Stochasticity

The above analysis is rigorous for an idealized elastic-brittle material; however, many real brittle materials exhibit stochastic failure due to their sensitivity to cracks and defects. This sensitivity is the result of the critical crack size to cause failure by fracture being on the order of the typical defect size in the material. Thus, these materials have three distinct length scales; the specimen and stress distribution scale, the lattice scale, and the microstructural defect scale. For a lattice composed of a brittle material, the resulting variability in failure can significantly affect the distribution of damage during failure, which changes the analysis from a deterministic problem to a probabilistic one. In the case of tensile failure [27], it has been shown that repeated tests of the same lattice results in damage propagation along varying paths, leading to a range of macroscopic failure behaviors. Stochastic material failure is also shown to significantly reduce the average tensile strength of a lattice [27].

Stochastic failure of materials is typically described by a cumulative distribution function (CDF), such as the Weibull distribution [26], which predicts a failure probability, $P_{\sigma}$, of a material as

[TABLE]

where $V$ is the volume of the material, $\sigma_{o}$ is the reference stress measured on a sample of volume $V_{o}$, and $m$ is the Weibull modulus, which describes the variance in failure behavior. Materials with larger Weibull moduli are less stochastic, as shown in Fig. 2(a). Many engineering ceramics, for example, have a Weibull modulus near $10$ [33]. For simplicity, in this work we will assume that the ligaments are under uniaxial tension and that their volume is comparable to the reference volume of the material, $V\approx V_{0}$. This simplifies the integral in eq. (4), although the framework presented can be easily repurposed without these simplifications. A Weibull distribution is used in this work due to its wide usage for brittle materials, but the analysis is directly adaptable to other distributions suitable for a particular system. The CDF is considered here, rather than the probability density function (PDF), as most real loading conditions will result in a stress increasing from zero to a finite non-zero value, rather than instantaneously achieving the final magnitude.

The CDF provides the probability that failure will occur as the material is loaded to some final stress state. For a lattice material, the probability that a ligament a point $(r,\theta)$ and oriented at an angle $\phi$ will fail at a given load is given by

[TABLE]

where the integral in eq. (4) was simplified using the assumption of uniaxial stress, and the ligament volume is assumed to be approximately equal to $V_{o}$. As the failure of the lattice is no longer deterministic, ligament failures cannot be predicted at any location, only the probability that a ligament will have failed can be determined, given the stress in the region and the orientation of the ligament.

2.3 Effect of Geometric Variability

The most difficult systems to predict where damage will form are materials with random or disordered structures. Unlike in a regular or periodic system, where a small number of possible features are available and the crack front is self-similar during propagation, a disordered system must consider a variety of possible structures and the probability of their appearance in the material as the crack advances. These details can be encoded in the analysis of the lattice microstructure by considering the ligament angles, $\phi_{i}$, as a set of random variables described by probability density functions, $f_{\phi_{i}}$, rather than as a set of deterministic discrete quantities. These probability density functions provide the likelihood of finding a ligament at any particular angle and can be tuned to match the microstructure of the material.

For the case of a triangular lattice, disorder can be introduced to the structure by perturbing the node locations. This will result in ligaments oriented at a set of random angles centered around the original discrete angles of the regular lattice but does not change the connectivity, preserving the assumption of stretch-dominated mechanics. The distribution of ligament angles can be described by a von Mises distribution [34], which is analogous to a normal distribution but is $\pi$-periodic, as shown in Fig. 2(b). The probability of finding a ligament at an angle $\phi$ is then given by

[TABLE]

where $I_{0}$ is the modified Bessel function of the first kind of order zero, $\mu$ is the mean around which the angle is centered, and $b$ controls the variance, with $b^{-1}$ being qualitatively analogous to the variance in a normal distribution. A geometrically deterministic lattice is therefore given by $b\rightarrow\infty$ (Fig. 2(b)). This distribution is provided as an example for a randomly disordered triangular lattice, but the approach is adaptable to other distributions that are tailored to specific material systems.

If a node of a lattice has a connected set of ligaments that are likely to be found distributed around a set of $N$ average angles $\mu={\mu_{1},\mu_{2},...,\mu_{N}}$, then the probability of finding a ligament at angle $\phi$ is given by the set of probability density functions, $f_{\phi}={f_{\phi_{1}},f_{\phi_{2}},...,f_{\phi_{N}}}$. If ligaments are described by a single distribution, $N=1$, but for the random triangular lattice $N=6$. As the possibility of each of these ligaments failing can be taken as independent of each other, the total probability of ligaments failing can be found by the product of their probabilities of failure. The total probability of failure can then be found by considering the complement of the probabilities of survival. Thus, for a node in the lattice at location $(r,\theta)$, the total probability of failure of any of the adjoining ligaments is then given by

[TABLE]

If the lattice is geometrically deterministic (i.e., $b\rightarrow\infty$), equation (7) simplifies to $P_{f}(r,\theta|m,\sigma_{o})=\prod\limits_{i=1}^{N}P_{\sigma}(r,\theta|\mu_{i},m)$, and only depends on the material stochasticity. If the ordered lattice is deterministic in its failure (i.e., $m\rightarrow\infty$), equation (7) further simplifies to a binary, where $P_{f}(r,\theta)=0$ if $\sigma_{a}\leq\sigma_{0}$, or $P_{f}(r,\theta)=1$ otherwise, which is the process outlined at the beginning of this work for a regular, deterministic lattice.

To quantify the probability that failure has occurred anywhere in the lattice, we treat the failure of the ligaments around any node as distinct from one another and consider the product of their probabilities. For a large lattice, the number of nodes tends towards infinity, suggesting an infinite product. To remain general to varying probability density functions that may not converge for an infinite product, we instead propose a Monte-Carlo approximation of the problem. Considering $s$ nodes located as a set of points $S=\{(r_{j},\theta_{j})|j=1,...,s,0<r_{j}\leq r_{max},-\pi\leq\theta\leq\pi\}$, where $r_{max}$ is the furthest radial extent of the specimen, the total probability of failure occurring anywhere in the specimen is given by

[TABLE]

The above approximation becomes an equality when $s$ accounts for all nodes in the lattice, with $s\rightarrow\infty$ as the lattice becomes large. This relationship provides the CDF that a ligament has failed anywhere in a random, stochastic lattice. However, the expression does not account for changes to the global stress distribution that arises from these failures, and so it can only be rigorously applied to the case of damage initiation, or the failure of one ligament. We demonstrate its predictive capabilities in this capacity through numerical modeling below. When considering a large number of ligament failures, such as if there is significant distributed damage during a fracture process, the model is only approximate, and becomes less accurate for larger number of ligament failures as the stress state will deviate in increasing amounts from the initial assumed state. Instead, we propose an alternative extension of the model to considering distributed damage later in the text.

3 Numerical Modeling

Lattice geometries were generated in the open-source computer-aided design software FreeCAD, implemented through a Python script (v. 3.10). A 30x30 triangular lattice with unit cells of length $L=5$ mm and ligament thickness $t=0.5$ mm were used, resulting in a 150 mm by 130 mm lattice. Two solid 10 mm beams were mounted on the top and bottom to apply a uniform displacement to the lattice, as shown in Fig. 3(a). To generate disordered lattices, nodes were perturbed at random angles with a controlled distribution of radial magnitudes to match the case of a von Mises distribution with $b=10$, as shown in Fig. 3(a). Models were then exported as 2-D DXF drawing files and imported into the finite element software.

Finite element simulations were performed in ABAQUS (2020, Providence, RI) using 2-D plane stress conditions (specimen thickness was taken as 6 mm), with the lower mounting point in the solid beams being fixed in $x$ and $y$, and one under displacement boundary conditions in $y$ and fixed in $x$. The models were meshed with bilinear plane-stress quadrilateral elements (CPS4), with a free mesh in the lattice and a uniform mesh in the beams, with a mesh size of 0.05 mm. Approximately 160,000 elements were used per simulation, depending on the exact geometry.

The material was modeled with an elastic modulus of $E=1$ GPa, and the reference stress was taken as $\sigma_{o}=20$ MPa. A 1 mm displacement was applied to the upper solid beam. The average von Mises stress over an area with a radius of $2t$ was calculated at the center of every ligament. The probability of survival was calculated for every ligament using eq. (4) and compared to a set of randomly selected real numbers between 0 and 1 using SciPy’s inbuilt random number generator. If any of the ligaments had a probability of survival lower than their associated random number, they are identified as a failed ligament. If multiple ligaments fail, then the ligament with the lowest probability of survival was taken as the broken ligament. If no ligaments fail, then the displacement is incremented by 0.5 mm, and a new probability of survival for each ligament is calculated. This process is iterated until a ligament fails. One hundred simulations were performed for the cases of an ordered lattice with $m=20$, an ordered lattice with $m=5$, and a disordered lattice ($b=10$) with $m=5$, for a total of three hundred simulations.

Validation of the finite element simulations and the analytic model was performed by comparing the prevalence of ligament failures in the simulations as a function of radial distance to that predicted by eq. (7). To model the failure throughout the entire lattice structure, a piecewise global stress distribution appropriate for an edge-cracked specimen was used. In this stress distribution, the K-field stresses given by eq. (2) are used near the crack tip which then decays to the the far-field applied stress, which is then taken as a fixed value determined by the boundary conditions.

The probability of ligament failure of the analytic model and the finite element simulations are shown in Fig. 3, as a function of radial distance from the crack tip. Results are given for two example cases of ordered lattices with two different Weibull moduli, $m=20$ and $m=5$. The simulations and analytic model show good agreement and demonstrate how damage is localized near the crack tip for the more deterministic case of a high Weibull modulus, and becomes more distributed for a low Weibull modulus. This is expected since the high stresses near the crack tip should result in a large probability of failure, but the lower stresses farther from the crack tip only result in a meaningful number of failures for the more stochastic material. There is some small disagreement between the model and simulations for the more stochastic case at short distances from the crack tip. This is expected due to the assumption of a singular stress state that results in a probability of failure approaching 1 for any Weibull modulus. As the finite element model with a crack tip of finite radius does not achieve a true singular stress, it is expected that the model should under-predict the probability of failure near the crack tip, relative to the singular analytic model, as is observed.

4 Results: Stochasticity, Disorder, and Length Scales

4.1 Failure Loads

The model in eq. (8) is rigorous for the initiation of damage (i.e., the first ligament failure of the lattice). Figure 4(a) shows the prevalence of failure loads measured from the finite element simulations for ordered lattices with Weibull moduli $m=20$ (pseudo-deterministic), and $m=5$ (stochastic), as well as a disordered lattice also with Weibull modulus $m=5$. As shown in Fig. 4(a), stochasticity has a significant impact on the average failure load and its variation, with the $m=20$ ordered lattice having a failure load of $750\pm 43$ N (mean $\pm$ std. dev.), while the ordered lattice with $m=5$ has a failure load of $420\pm 95$ N. Note that the relative magnitude of the standard deviation to the mean failure load increases by nearly a factor of four between the two cases. The introduction of disorder to the $m=5$ case results in no significant change in failure load or variation, with an average failure load of $425\pm 80$ N.

By taking a discrete derivative of eq. (8) with respect to the far-field stress, we can calculate the PDF of the failure loads. These predictions are given in Figure 4(a) and show excellent agreement with the simulation results in all cases, accurately capturing both the mean and variance of the failure loads. Also shown are best-fit models to the data, considering either a Weibull distribution or normal distribution, as was done in [27]. A maximum log-likelihood criterion was used to determine the quality of the fit. For the ordered lattices the best-fit model is the Weibull distributions with fitted moduli values of 20.0 and 5.0, exactly. This is expected from eq. (8), as it is simple to see that in the limit of an ordered lattice the probability density functions for $\phi$ collapse to delta functions, resulting in a simple product of Weibull distributions. For the disordered lattice, the normal distribution achieved a better quality of fit, as was also observed in [27]. This is now explained by eq. (8), which shows that, while the data is not expected to truly follow a normal distribution, it will also not follow a pure Weibull distribution as the stochastic behavior is modulated by the geometric disorder.

4.2 Fracto-cohesive Length

The results in Fig. 4(a) demonstrate that the analytic model accurately predicts the required load to initiate damage in a lattice material. This result can be utilized to calculate additional parameters that characterize the failure of the lattice, such as the fracto-cohesive length, $\lambda$. If any defect in the lattice is larger than $\lambda$, the material will undergo a fracture-based failure, rather than a strength-based failure [30]. Knowing which regime a system is in can be useful for making predictions about failure. For example, if a lattice will undergo a fracture-based failure, the failure loads will strongly depend on the defect size, unlike in the strength-based failure regime.

The fracto-cohesive length is defined [30] as

[TABLE]

where $W_{f}=\int_{0}^{\epsilon_{c}}\sigma(\epsilon)d\epsilon$ is the work of fracture for an uncracked material, where $\epsilon_{f}$ is the failure strain. For a linear elastic material that fails in a brittle manner under uniaxial tension, this expression simplifies to $W_{f}=\frac{\sigma_{f}^{\ast^{2}}}{2E}$, where $\sigma_{f}^{\ast}$ is the failure stress of the lattice material. Note that this is distinct from the failure stress of the base material of the lattice, $\sigma_{f}$. This can be calculated directly from eq. (8), where the local stress field is taken as a constant ($\sigma(r,\theta)=\bar{\sigma}$).

The toughness, $G_{c}$, can be estimated from eq. (8) as well using the Irwin-Kies relationship, with $G_{c}=K_{IC}^{2}/E$, in plane stress, where $K_{IC}=Y\sigma_{c}^{\ast}\sqrt{\pi a}$ is the critical stress intensity factor, with $Y\approx 1$ is a geometry factor, and $\sigma_{c}^{\ast}$ is the critical stress to cause fracture for the given crack length, which can be predicted from eq. (8). Therefore, the fracto-cohesive length can be estimated as

[TABLE]

Figure 4(b) shows predictions of the fracto-cohesive length from eq. (9), where $G_{c}$ is calculated from eq. (8) using the K-field stress distribution, and $W_{c}$ is calculated using a uniform tensile stress. Results are normalized by the predicted fracto-cohesive length for an ordered, deterministic lattice of equivalent density, $\lambda_{0}$. The results show that the fracto-cohesive length depends strongly on the Weibull modulus, as shown in [28], but that had not been quantitatively explained before. Further, there is a dependence on the geometric disorder, although it’s effect is far smaller. This makes sense since the stresses can only be affected marginally by the angles of the ligaments if they are randomly oriented, whereas the changes in stochasticity can significantly impact the probability of failure given a local stress.

4.3 Representative Volume Element

Beyond the fracto-cohesive length, the size of a full representative volume element (RVE) for a lattice can be calculated. This has a theoretical utility as a test of whether a lattice is sufficiently large for the model to be accurate. If a lattice is much smaller than its RVE, the assumption of a pseudo-continuous stress field (eq. (3)) will not apply. A practical utility of the calculation is to determine if a lattice is much larger than the RVE, as its mechanics can then be homogenized, which is useful for characterizing the global behavior of the material, as shown below.

Considering a simple rectangular RVE that spans the width of the material, $b$, with an area $h\xi$, where $\xi$ is the length of the RVE along the crack path ($\theta=0$) and $h$ is the height of the RVE perpendicular to the crack path ($\theta=\pi/2$). In a strict sense, the length of the RVE, $\xi$, needs to only span the distance over which the crack advances. In the case of an ordered, deterministic lattice material, the crack will jump in discrete intervals equal to the length of a unit cell, $L$. For a lattice material with disorder, this may not be the case and damage may occur in a unit cell not directly in front of the crack tip. By defining the RVE length using the first-order Irwin fracture stress length scale, $\lambda_{0}=(2\pi)^{-1}(K_{IC}/\sigma_{f})^{2}$, ensures that the entire region of high stress around the crack tip is encompassed where most damage is expected to occur. Inspired by the approach of Irwin, this length scale can be adapted to a stochastic material by considering a length ahead of the crack tip, $\xi$, where the probability of failure is $T$ (taken to be sufficiently close to $1$). From eq. (5) for $\theta=0$, it is easy to show that

[TABLE]

It is clear that as $m\rightarrow\infty$, $\xi/\lambda_{0}\rightarrow 1$, as expected. The choice of $T$ is arbitrary, and so this length scale should not be treated as an inherent material property like $\lambda$; however, for reasonable values of $m$ ($m>1$), choices of $T$ between 90% and 99.99% are of the same order of magnitude and can therefore provide an estimate of the RVE. In this work we will consider the conservative case of $T=0.9$, which will create a larger RVE.

The height of the RVE defines the region above the crack in which the majority of the energy is released during fracture. Following the calculation in [4], the height of the RVE can be found as $h=G_{c}/W_{c}$, where $W_{c}=\frac{1}{2}\sigma_{ij}\epsilon_{ij}$ is the critical strain energy density in the RVE at the point of crack propagation. Using the K-field stresses given in eq. (2), and the corresponding strains for the case of plane strain, the critical strain energy density in front of the crack ($\theta=0$), becomes

[TABLE]

where $E^{\ast}$ and $\nu^{\ast}$ are the effective modulus and Poisson’s ratio of the lattice material. For simplicity, we consider the case of a geometrically ordered, stochastic lattice, and estimate the fracture toughness, $K_{IC}$, using the stress at which a ligament fails at $r=\xi$. (Note, a more thorough estimate of $K_{IC}$ is presented below.) In an ordered lattice, at the point $r=\xi$ and $\theta=0$, the first ligament to fail is the one oriented at $\phi=\pi/3$, resulting in $K_{IC}\approx\sigma_{f}\sqrt{2\pi\xi}$. Substitution of this result into eq. (12) leads to

[TABLE]

Using the Irwin-Kies relationship to define the lattice toughness as $G_{c}=K_{IC}^{2}/E^{\ast}$, the height of the RVE can be found as

[TABLE]

For the case of a triangular lattice, $\nu^{\ast}\approx 1/3$ and $h=3\pi\xi$, or the aspect ratio of the RVE is $h/\xi=3\pi$. For triangular lattices of with unit cells of width $L$, the height of a unit cell is $\sqrt{3}L/2$. Thus, if the RVE width is $\xi=N$ unit cells, then the height is $h=2\sqrt{3}\pi N\approx 11N$. For the lattices considered in our finite element simulations, the vast majority of ligament failures are localized to the first $2.5\pm 0.5$ unit cells ahead of the crack tip, as shown in Fig. 3, which suggests that the RVE height must be $27\pm 5$ unit cells high.

Lattices of varying numbers of unit cells were simulated to determine what the minimum RVE size should be and to verify the prediction of $h$ from the model. Ordered lattices with deterministic failure were considered. Minimum RVE size can be determined independently by measuring the local ligament stresses and comparing to the predictions of a K-field stress. Fig. 5(a) shows that the average difference in ligament axial stresses between the finite element model and the analytic predictions from eq. (3) as a function of the radial distance from the crack tip. Results are given for four lattices of varying size, with $N=12,18,24$ and $30$, as shown in Fig. 5(b-e). For all of the lattices there is disagreement between the analytic model and the finite element model for radial distances smaller than one unit cell in length as $r\rightarrow 0$. This is expected since the idealized K-field stresses become infinite, which is not physical. However, at greater radial distances it is clear that the large lattices agree more accurately with the analytic model. The $N=12$ lattice maintains an approximately 20% deviation from the analytic model across the entire lattice, and the $N=18$ lattice has a nearly 10% error across most of the lattice. It is only for $N=24$ that the error drops below 5% and for $N=30$ the error converges to $<1\%$ for larger radial distances. This is in line with the predictions of $h$ required to obtain the asymptotic crack tip field in an architected material and is an independent verification of the model prediction.

5 Results: Toughness and Damage Zones

For the case of a lattice that is much larger than its RVE, its mechanics can be homogenized, which allows the model to be extended beyond predicting failure loads to estimating fracture toughness. As noted above, eq. (8) is rigorous in predicting failure for the initiation of damage in a lattice structure. However, it is apparent from the results in Fig. 3(b) that significant distributed damage can occur in these lattices, suggesting that the steady-state toughness of the material would likely not be achieved until a large number of ligament failures occur, such that the K-field stresses shift along the crack path. As more ligament failures occur, eq. (8) will become less and less accurate as it does not account for the changing stress state due to the change in local connectivity throughout the lattice. However, in the case of a large lattice where the mechanics are homogenized, the model can instead be extended to predict the toughness by using the distribution of damage, as has been shown previously to be an accurate predictor of the effective toughness of the material [14, 22].

While in a stochastic lattice damage may theoretically occur anywhere, the singular nature of the K-field stresses ensures that the majority of damage occurs near the crack-tip (see Fig. 3(b)). Since the stress fields are monotonically decreasing as the radial distance increases, it enables us to define a region within which there is a threshold level of probability of failure, $T$, outside of which the probability of failure is lower. This is then analogous to a process zone in a homogeneous material.

The probability of failure is calculated using eq. (7) as the radial distance from the crack tip increases and for all $\theta$. A threshold probability of failure, $T$, is chosen to define the damage zone. In this work a probability of failure $T\geq 0.9$ is used. The choice of the threshold is arbitrary, although we note that there was not a significant change in the damage zone sizes as the threshold was increased due to the singular nature of the stress fields. As a control study, the plane stress deterministic damage zones for a homogeneous material, and an ordered lattice are shown in Fig. 6(a). The lattice damage zone follows closely with the homogeneous material, with perturbations due to geometric coupling with the six individual ligament angles.

Figure 6(b) shows the process zones inside of which the probability of failure $F\geq 0.9$ for the cases given in Fig. 4(a) ($m=20$ and $m=5$ ordered lattices and a $m=5$ disordered lattice, $b=10$). It is apparent that as the network becomes more geometrically disordered, the process zone continues to increase in size as a result of the variability in the ligament stresses due to the geometric disorder, and the variability in failure stresses that result from the stochastic failure. It is also apparent how the introduction of geometric disorder and stochastic failure results in a more isotropic process zone, with the angular dependence seen in the regular lattice replaced by a more uniform zone.

5.1 Estimating Material Toughness

As shown in Fig. 6, increasing geometric variability and material stochasticity both tend to increase the size of the process zone in a lattice material. Larger process zones result in an increase in the number of ligament failures, increasing the energy dissipated during crack propagation. As shown in [13, 14], the toughness depends on the size of the process zone and maximum local strain energy density, with $G_{c}\propto\rho_{e}\Omega$, where $\rho_{e}$ is the energy density of the process zone during failure, and $\Omega$ is the volume of the process zone.

For a linear elastic fiber or lattice material, the maximum energy density in the process zone depends on the local stress before failure. If failure is stochastic and determined by a probability, $P_{s}$, then the average ligament energy density will be $\rho_{e}=\mathbb{E}_{\sigma^{2}}$, where $\mathbb{E}_{\sigma^{2}}=\int_{0}^{\infty}P_{s}^{\prime}(\sigma|m)\sigma^{2}d\sigma$ is the expectation value of $\sigma^{2}$, where $P_{s}^{\prime}(\sigma|m)=\frac{d}{d\sigma}P_{s}(\sigma|m)$. Thus, the toughness of a stochastic disordered lattice, relative to its deterministic regular analog is

[TABLE]

where $G_{o}$ and $\Omega_{o}$ are the toughness and process zone size of the deterministic lattice, respectively. If $P_{s}$ is the Weibull distribution, the expectation value will be $\mathbb{E}_{\sigma^{2}}=\sigma_{o}^{2}\Gamma\left(\frac{2+m}{m}\right)$, where $\Gamma$ is the Gamma function, which will be $\leq 1$ for all positive values of $m$. Thus, $\mathbb{E}_{\sigma^{2}}/\sigma_{o}^{2}\leq 1$ and the introduction of material stochasticity will reduce the local strain energy density, which in turn will have the tendency to lower the toughness. However, as shown in Fig. 6(b), the process zone volume increases with stochasticity, tending to increase the toughness. Introducing geometric variability will also increase the process zone but has a much smaller effect on the local energy density. Thus, geometric variability is expected to generally increase the toughness, while stochasticity introduces competing effects. Specifically for the Weibull distribution in a triangular lattice, the loss in local strain energy will quickly become comparable to the increase in process zone size as $m$ decreases, suggesting limited enhancements in toughness as a result of the loss of strain energy, even as the damage zones are enlarged.

6 Conclusions

The failure of lattice metamaterials with geometric disorder and stochastic material failure is predicted from an analytic model that reveals the effects of local geometry and material properties on the mean and variance in the failure loads of the material, as well as the fracto-cohesive length and the minimum RVE size that control whether a fracture-like failure occurs in the material. Increasing stochasticity significantly reduces the failure loads and results in large fracto-cohesive lengths controlling whether failure depends on initial defect size. Geometric disorder has a smaller effect on failure but changes the variation in failure loads from being a pure Weibull-like distribution to a more complicated model better approximated by a normal distribution. An extension to a full fracture model is presented that uses predictions of the damage zones around the crack tip to estimate the toughness of the lattice. Both stochasticity and disorder result in larger damage zones around the crack tip, but stochasticity significantly reduces the local strain energy density, which reduces the toughness and suggests that while disorder does not significantly influence the initiation of failure, it provides an avenue for enhancing the steady state toughness of the lattice.

These results demonstrate how the failure of lattice mechanical metamaterials is influenced by the interplay of the base material failure properties, such as stochastic failure, and the microstructural geometry. Results are presented for a precracked triangular lattice with Weibull stochastic failure, but the framework is easily adaptable to varying lattice geometries, failure criteria, and stress distributions. An analytical model provides predictions of the failure stress and effective toughness of a lattice, as well as details about the nature of the failure (i.e., strength-based or fracture-based). Unlike with data-driven methods, the analytical approach here provides insights into how the global mechanical behavior of metamaterials arises from the combination of material properties and microstructural geometry, along with the quantitative predictions of the failure.

7 Acknowledgments

This research was funded primarily by the National Science Foundation (NSF) MRSEC program [awards DMR-2309043 and DMR-1720530] to K.T.T and P.K.P. S.F. acknowledges support from the US Department of Defense (DoD) through the National Defense Science & Engineering Graduate (NDSEG) Fellowship Program.