---
abstract: 'The toughness of a polymer material can increase significantly if two networks are combined into one material. This toughening effect is a consequence of a transition from a brittle to a ductile failure response. Although this transition and the accompanying toughening effect have been demonstrated in hydrogels first, the concept has been proven effective in elastomers and in macroscopic composites as well. This suggests that the transition is not caused by a specific molecular architecture, but rather by a general physical principle related to the mechanical interplay between two interpenetrating networks. Here we employ theory and computer simulations, inspired by this general principle, to investigate how disorder controls the brittle-to-ductile transition both at the macroscopic and the microscopic level. A random spring network model featuring two different spring types, enables us to study the joined effect of initial disorder and network-induced stress heterogeneity on this transition. We reveal that a mechanical force balance gives a good description of the brittle-to-ductile transition. In addition, the inclusion of disorder in the spring model predicts four different failure regimes along the brittle-to-ductile response in agreement with experimental findings. Finally, we show that the network structure can result in stress concentration, diffuse damage and loss of percolation depending on the failure regime. This work thus provides a framework for the design and optimization of double network materials and underlines the importance of network structure in the toughness of polymer materials.'
author:
- Justin Tauber
- Simone Dussi
- Jasper van der Gucht
bibliography:
- 'ms.bib'
title: Microscopic insights into the failure of elastic double networks
---

Introduction
============

Polymer networks, such as rubbers and gels, can undergo large deformation without losing elasticity. This property makes them ideal for numerous applications, for example in the biomedical field or in soft robotics. One of the factors that limits their applicability, however, is the brittleness of many polymer networks, which leads to inferior mechanical performance, such as low fracture toughness and strength [@Creton2017]. In recent years, several strategies have been developed to toughen polymer networks by introducing dissipation mechanisms that delay crack propagation. One of the most successful strategies relies on the combination of two different polymer networks into one material, to create a so-called double network. In particular, extremely tough double network hydrogels have been produced by interpenetrating a stiff and weak first network with a soft and extensible second network [@Gong2010; @Nakajima2013; @Gong2003; @Weng2008; @Sun2012; @Nakajima2012; @Xin2013; @Nakajima2013a; @Chen2014; @Yan2017] (Fig. \[fig:Mechanism\](a)). Later, this toughening strategy was also shown to be effective for multi-network elastomers [@Ducrot2014; @Millereau2018] or macroscopic composites [@Takahashi2018; @King2019](Fig. \[fig:Mechanism\](b)). This suggests that the underlying principle that leads to toughening is not strongly dependent on molecular details, but governed more generally by an interplay between two mechanically different networks.

![ Double networks across the scale. (a) By combining a sacrificial network that is stiff and weak (blue), with a matrix network that is soft and strong (red), a (molecular) composite both stiff and strong (tough) can be created. (b) Typical double networks are made of elastomers [@Ducrot2014] and hydrogels [@Gong2010], but they can also be made in the form of macroscopic composites [@Takahashi2018]. (c) Schematic of a brittle (blue) and a ductile (yellow) stress-strain response. []{data-label="fig:Mechanism"}](figure1.eps)

The increase in toughness of double networks is related to a transition from brittle to ductile failure. This brittle-to-ductile transition (BDT) becomes apparent when looking at the stress-strain response of double networks subjected to tensile deformation, as schematically shown in Fig. \[fig:Mechanism\](c). In contrast to single polymer networks that typically fail in a brittle fashion (blue curve), tough double networks fail in a ductile manner with significant softening and plastic deformation (leading to ‘necking’) prior to failure (yellow curve) [@Matsuda2016; @Millereau2018]. This macroscopic softening is believed to be due to the progressive breaking of bonds in the stiff network (which is also called the sacrificial network), while the soft matrix formed by the second network stays intact. To obtain tough materials, it is thus crucial that the rupture of the first network occurs without leading to macroscopic cracks in the material.

One way to go forward is to include molecular details, for example using molecular dynamics simulations [@Jang2007; @Wang2017; @Higuchi2018]. However, given the computational costs it is difficult to do this for large system sizes and to explore a large range of parameters. Alternatively, statistical models, such as the fiber bundle model, can give insight into the role of disorder on network failure [@Kovacs2013; @Roy2017]. In such models, however, a priori assumptions must be made about the (re)distribution of stress among the different elements, so that they cannot explain the spatial evolution of damage in the materials and its relation to the structure of the networks.

Here, we develop a random double spring network model, which explicitly takes into account disorder and the resulting heterogeneous stress distributions. Rather than making assumptions about the distribution of strain and stress in the material, this distribution emerges naturally from the condition of mechanical equilibrium between the different networks. Because a network is included explicitly, the model provides information on the failure process both at the macroscopic and the microscopic level. We show how this model predicts a simple criterion for the BDT and we show how the nature of this transition is influenced by disorder and the resulting stress localization. We also compare our network model to a simple 1D multi-spring model, which can be solved analytically.

\[sec:Models\] Models and Methods
=================================

The BDT is the universal mechanical feature of double network materials. In a brittle material, a small microcrack in the sacrifial network directly leads to macroscopic failure. On the contrary, a ductile material remains intact when microcracks develop. To describe this transition, we first present a simple 1D multi-spring model, which is an extension of a previously used two-spring model [@Ahmed2014; @Okumura2004] and which can be solved analytically (Fig. \[fig:ModelOverview\](a)). While heterogeneity can be included in this 1D model, it does not take into account the network structure of the material and therefore cannot give an accurate description of the stress distribution in the material. We therefore consider a more realistic double random spring model that does allow for a heterogeneous stress distribution in the material (Fig. \[fig:ModelOverview\](b)).

 Multi-spring model
-------------------

![Minimal models for elastic double networks and their reaction to failure. The springs represent elements of the sacrificial network (blue) and matrix network (red). (a) Multi-springs (MS) model. (b) Random Spring Network (RSN) model. Simulations are performed for networks of $50\times50$ nodes.[]{data-label="fig:ModelOverview"}](figure2.eps)

In the multi-spring model (MS), the double network is modeled as a series of elements, each consisting of two parallel springs, one representing the sacrificial network (with spring constant $\mu_S$) and one representing the matrix (with spring constant $\mu_M$) (see Fig. \[fig:ModelOverview\](a)). Each spring is assumed to be a linear Hookean spring, so that the force acting on it is given by $F_i=\mu_i\Delta l_i$, with $\Delta l_i$ the extension of spring $i$. Removal of a sacrificial spring in an element leads to the formation of a microcrack, which is bridged by the remaining matrix spring. As the elements are connected in series, the load on each element is equal. This force balance between the elements implies

$$\label{forcebalance}
\sigma=(\mu_M+\mu_S)\Delta l_D = \mu_M\Delta l_M \; ,$$

where $\Delta l_D$ and $\Delta l_M$ represent the extension of the network in the intact double network and in the matrix spanning the microcrack, respectively, and where $\sigma$ denotes the stress. The overall strain of the system $\epsilon$ can be written as

$$\label{strain}
\epsilon=\Delta l_D(1-\phi)+\Delta l_M \phi  \; ,$$

where $\phi$ denotes the fraction of broken sacrificial bonds (i.e. the fraction of microcracks).

Rupture of a bond is assumed to occur instantaneously when the extension of a spring $\Delta l_i$ exceeds its threshold $\lambda_i$. In the absence of disorder, all sacrificial bonds have the same threshold $\lambda_S$ and all matrix bonds have the same threshold $\lambda_M$. Brittle failure after formation of a microcrack occurs when the bridging matrix bond reaches its threshold before the other sacrificial bonds, while ductile failure occurs when the sacrificial bonds reach their threshold first. From equation \[forcebalance\] it immediately follows that the brittle-to-ductile transition occurs when $(\mu_S+\mu_M)\lambda_S=\mu_M\lambda_M$. We can therefore predict which system parameters are required for a system at the BDT. For example, the sacrificial bond threshold $\lambda_S^{*}$ for which the rupture force of a (matrix-reinforced) sacrificial bond is equal to the rupture force of a matrix bond is given by

$$\lambda_S^{*} = \frac{\mu_M\lambda_M}{\mu_S + \mu_M}  \;.$$

More generally, we can define a parameter $\Delta\alpha$ to quantify the distance from the BDT:

$$\label{Da}
\Delta\alpha=\frac{\mu_M\lambda_M}{(\mu_M+\mu_S)\lambda_S}-1  \; .$$

If $\Delta\alpha=0$ the system is at the BDT. Brittle failure occurs for $\Delta\alpha<0$; in this case only one microcrack is enough to cause global failure so the fraction of broken sacrificial bonds at failure is $\phi_f=1/N$ with $N$ the system size. Ductile failure occurs for $\Delta\alpha>0$, where all sacrificial bonds break before the system fails globally, so $\phi_f=1$. It follows from this analysis that for ductile failure we need $\lambda_S<\lambda_M$, which explains why the creation of tough, ductile networks requires the sacrificial network to be much weaker than the matrix.

For a perfectly homogeneous system, the BDT in this model is an abrupt transition. Experiments, however, show that the descriptors of failure vary in a continuous manner through the BDT [@Matsuda2016]. We therefore introduce disorder by assuming that the thresholds of the sacrificial network vary according to a certain distribution $P(\lambda_S)$. All other parameters are taken the same. In this case, the sacrificial bonds will fail progressively, with the weakest bonds breaking first and the stronger bonds breaking later. At a given strain $\epsilon$, all sacrificial bonds for which $\lambda_S<\Delta l_D$ are broken. In the continuous limit (for $N\gg 1$), the fraction of broken sacrificial bonds can thus be written as

$$\label{dist}
\phi(\epsilon)=\int_0^{\Delta l_D} P(\lambda_S)\textrm{d}\lambda_D  \; .$$

Here, we will consider a Gaussian threshold distribution with mean $\langle\lambda_S\rangle$ and standard deviation $\delta\lambda$. For this case, it follows from equations \[forcebalance\] and \[dist\] that at the moment of macroscopic failure (where $\Delta l_M=\lambda_M$), the fraction of broken sacrificial bonds is given by

$$\label{dist2}
\phi_f=\frac{1}{2}\left[1+\textrm{erf}\left(\Delta\tilde{\alpha}\right)\right]  \; ,$$

with

$$\label{Datilde}
\Delta\tilde{\alpha}=\frac{\Delta\alpha\langle\lambda_S\rangle}{\delta\lambda\sqrt{2}}  \; ,$$

a normalized parameter to quantify the distance to the BDT, which takes into account the disorder in the thresholds. The fraction of broken sacrificial bonds thus increases gradually from 0 to 1 along the BDT and the parameter $\Delta\tilde{\alpha}$ sets the steepness of this transition. With increasing disorder, the transition becomes more gradual. The strain at break follows from equations \[forcebalance\] and \[strain\]:

$$\label{eq:failureStrain}
\epsilon_f=\lambda_M\left(\frac{\mu_M+\phi_f\mu_S}{\mu_M+\mu_S}\right)  \;.$$

In the brittle regime ($\Delta\tilde{\alpha}\ll-1$, $\phi_f\approx0$) we find $\epsilon_f=\mu_M\lambda_M/(\mu_M+\mu_S)$, while in the ductile regime ($\Delta\tilde{\alpha}\gg1$, $\phi_f\approx1$) we have $\epsilon_f=\lambda_M$.

Random Spring Network Model
---------------------------

Being a 1D model, the multi-spring model cannot account for stress heterogeneity and the resulting localization of stresses. We therefore consider a random spring network model (RSN) as shown in Fig. \[fig:ModelOverview\](b) composed of $L \times L$ nodes arranged on a triangular lattice whose nearest neighbours are connected by the same element we introduced for the multi-spring model ($L=50$ for all RSN simulations). Now, upon applying an external load, stress concentration is possible due to the topological restrictions imposed by the network structure.

Following the experimentally found guidelines for making tough double networks [@Gong2010] we consider materials in which the sacrificial network is stiff and weak, while the matrix is soft and strong. This means that the elastic constant of sacrificial springs $\mu_S$ is always higher than the one of the matrix springs $\mu_M$. We therefore vary the ratio of the stiffnesses such that $0<\mu_M/\mu_S<1$. To reflect the asymmetry in network strengths, we fix $\lambda_M=4.0$ and vary $\lambda_S$ such that $\lambda_S<\lambda_M$. In all cases, we keep the parameters of all matrix springs the same. To implement disorder at the bond-level, the failure thresholds $\lambda_S$ are picked from a Gaussian distribution as in the MS model. Nodes at the bottom and top are fixed in the $y$-direction, but they can slide along the $x$-direction. Along the $x$-axis periodic boundary conditions are implemented.

We characterize the mechanical response of the network by applying an extensional strain $\epsilon$ along the y-axis in small steps of 0.1% strain. We consider quasi-static loading and assume that the networks remain in mechanical equilibrium at each step, i.e. it settles in its minimum energy state. Therefore, after every step, the total energy of the network is minimized by displacing the nodes using the FIRE algorithm [@Bitzek2006] with a tolerance of $F_{rms}\leq1\cdot10^{-5}$. Here, $F_{rms}$ is the maximum root mean squared force allowed in the system [@Bitzek2006]. To simulate failure, we break all bonds which exceed their failure threshold $\lambda_i$ consecutively. After the failure of every single bond the energy of the system is minimized. Once all overstretched bonds are broken, the next strain step is determined according to the bond that is closest to its failure threshold with a minimum step size of 0.001% strain.

For every strain in the simulations the virial stress is calculated from the forces exerted by the springs on the nodes. As the system is elongated along the $y$-direction, we consider the $yy$-component of the virial tensor as a measure for the stress $\sigma$. We define the softening strain $\epsilon_{soft}$ as the strain where the stiffness drops below the initial stiffness, and the failure strain $\epsilon_f$ and failure stress $\sigma_f$ as the strain and stress for which the system becomes disconnected along the y-direction. For every data point 50 simulations are performed, errorbars represent standard deviation.

\[sec:ResultsDiscussion\] Results & discussion
==============================================

![Failure strain, $\epsilon_f$, versus the average failure threshold $\langle \lambda_S \rangle$ in systems with $\delta\lambda=0.250$. Results are shown for the RSN model (solid lines) and the MS model (dashed lines) for stiffness ratios $\mu_M/\mu_S =$ 0.10 (purple), 0.50 (green) and 1.00 (yellow)).[]{data-label="fig:FailureStrain"}](figure3.eps)

![Typical failure response of double networks obtained via the RSN model. (a) The four distinct stress-strain responses: brittle (B), quasibrittle (QB), the first ductile regime (D1), and the second ductile regime (D2). (b)-(e) Corresponding failure patterns. Every line represents a broken sacrificial bond, color-coded according to the strain at which it failed (see color bar above).[]{data-label="fig:StressStrain"}](figure4.eps)

\[sec:FailureRegimes\] Failure regimes of double networks
---------------------------------------------------------

The aim of this paper is to systematically explore the effect of disorder and stress heterogeneity on the failure of double networks, focusing especially on the role of the network. To do this in a systematic fashion, we will extensively compare the results of the random spring network (RSN) model with the 1D multi-spring (MS) model lacking any network structure.

We start our comparison by checking a prediction from the MS model for the dependence of the failure strain $\epsilon_f$, described in equation \[eq:failureStrain\], on the sacrificial bond properties $\lambda_S$ and $\mu_M/\mu_S$ in the presence of disorder (Fig. \[fig:FailureStrain\]). Although it is clear that both the MS model and the RSN model follow the same trend, the overall trend is not immediately obvious. Intuitively, it is expected that the failure strain should increase with an increase in $\lambda_S$. However, only at the lowest stiffness ratio, we see that the $\epsilon_f$ increases with $\lambda_S$. At high stiffness ratios, this trend is reversed. This transition is a consequence of the brittle-to-ductile transition (BDT). From equation \[Da\] it can be deduced that an increase in $\lambda_S$ brings the system from the ductile into the brittle regime. The transition is not observed at the lowest stiffness ratio, because all those systems are already in the brittle regime. Similar trends can be found for the material strength and the work of extension (see Supplemental Material [@supp]).

These results show that we can tune the failure response based on the properties of the individual sacrificial bonds. Taking a broader view at the entire stress-strain response of the RSN, we can differentiate between four distinct responses, shown in Fig. \[fig:StressStrain\](a). There are two brittle responses, brittle (B) and quasibrittle (QB), that can be distinguished by the (slight) softening that occurs before macroscopic failure in the QB case. There are also two ductile responses (D1 and D2) that both show a clear plateau in the stress-strain response after softening. In case of D1, macroscopic failure occurs within this plateau. For D2, a second increase in the stress precedes final failure. These same four responses are also found in the MS model (see Supplemental Material [@supp]).

By contrast, at the microscopic level the failure behaviour of the RSN is distinct from the MS model. In Fig. \[fig:StressStrain\](b-e), we plot the broken bonds of the sacrificial network color-coded according to the strain at which they break. We clearly observe that the number of broken sacrificial bonds increases from B to D2 type of fracture. In B, all broken bonds are part of a single crack at the location where macroscopic failure occurs. A similar pattern is observed for QB, where however more bonds fail homogeneously throughout the sample before the final crack appears. It is the failure of these bonds that causes the softening observed in the QB stress-strain response. The ductile regimes are characterized by many bonds breaking simultaneously. Such avalanches occur over a large strain range, from the moment of softening until macroscopic failure. In case of D2 fracture, all the non-horizontal bonds of the sacrificial network, which are the elastically active bonds, are broken. In the MS model we can not differentiate between failure patterns, because the load is homogeneous throughout the system, so that stress can not concentrate and damage will always be homogeneous.

We note that both the RSN model and the MS model capture the four different mechanical responses observed in experiments on elastomers [@Millereau2018], hydrogels [@Gong2010] and macroscopic composites [@Takahashi2018; @King2019]. The clear difference in fraction of broken sacrificial bonds between (quasi)brittle and ductile responses is also well established experimentally [@Matsuda2016; @Millereau2018]. Furthermore, the softening response, which is irreversible both in our model and in the experiments, is attributed to early damage of bonds in the sacrificial network [@Vernerey2018] which has been shown to occur homogeneously throughout the network well before macroscopic failure [@Millereau2018]. Finally, also a plateau in the stress response is typically observed in experiments after the yielding of the material, often accompanied by a necking region where only sacrificial bonds break [@Millereau2018]. We conclude that both the RSN and MS model represent a minimal yet insightful model to study fracture in elastic double networks both at the macroscopic and microscopic level and therefore we can proceed to a more detailed analysis.

\[sec:MacroResponse\] Macroscopic characterization of the failure regimes
-------------------------------------------------------------------------

![image](figure5.eps)

The results in Fig. \[fig:FailureStrain\] clearly demonstrate the dependence of the failure response on the properties of the sacrificial bonds: $\langle \lambda_S \rangle$, and $\mu_M/\mu_S$. To study the joined effect of these properties together with the disorder in the failure threshold ($\delta\lambda$), we construct phase diagrams of the failure response (Fig. \[fig:PhaseDiagram\]) based on the four different regimes described above (see Supplemental Material [@supp]). For each point in the diagram we use a pie-chart symbol indicating the frequency of the four different regimes observed over 50 independent simulations. We first focus on the case with no disorder $\delta\lambda=0$, shown in Fig. \[fig:PhaseDiagram\](a), and observe that both $\mu_M/\mu_S$ and $\lambda_S$ control the sharp transition from the BDT. Also, the BDT of the RSN model coincides with the BDT of the MS model (Eq. \[Da\]) as indicated by the black line.

Next, we consider the cases with disorder, plotted in Fig. \[fig:PhaseDiagram\] (b) and (c), and we observe the appearance of the intermediate regimes QB and D1. These two regimes appear around the BDT and upon increasing disorder $\delta \lambda$ progressively span the entire diagram, which is a clear indication that disorder controls the position of these transitions. In addition, the transitions between the regimes are less sharp when the disorder is larger, as testified by the less homogeneously colored pie-chart symbols. Indeed, the mechanical responses of networks with large $\delta \lambda$ depend on the exact realization of the network (exact distribution of threshold values and their spatial organization), even when having the same $\mu_M/\mu_S$ and $\left<\lambda_S\right>$, especially around boundaries between different regimes. A similar dependence of the location of the B-to-QB and D1-to-D2 transitions on disorder is also predicted by the MS model (dashed lines, see Supplemental Material for further information [@supp]). At low disorder, the location of the boundaries corresponds well between the models, but at high disorder the boundaries move away from each other. This comparison reveals that network structure is not essential for the occurrence of any of the four regimes, but does have an effect on its location. Furthermore, disorder clearly plays an important role also in the experimental systems since QB and D1 responses are typically observed experimentally and we show that these occur only in presence of disorder. Our results confirm that knowing the average value of the bond strengths might not be enough to predict the failure regime, in agreement with recent phenomenological models aiming to capture softening and the Mullins effect in double networks [@Vernerey2018].

\[sec:MicroResponse\] Failure at the microscopic level
------------------------------------------------------

In the previous section, we found that the macroscopic force balance reasonably captures the transition from (quasi)brittle to ductile fracture in the RSN model and that disorder influences this transition. Furthermore, we have seen in Fig. \[fig:StressStrain\](b)-(e) that as long as some disorder is present distinct failure patterns develop in the RSN that suggest a tight link between the macroscopic and microscopic failure process. In this section, we explore several indicators of microscopic failure in disordered RSN systems, with a special focus on stress heterogeneity, and we consider to what extent the force balance influences the microscopic failure processes in double networks.

### The fraction of broken bonds

![Rescaling of the fraction of broken bonds $\phi_f$. (a) The fraction of total broken bonds $\phi_{f}$ obtained in simulations of RSN is plotted against $\Delta\alpha$, which quantifies the distance from the BDT predicted by the equal load sharing models. (b) $\phi_{f}$ from RSN versus $\Delta\tilde{\alpha}$. The dashed black line indicates the theoretical prediction for $\phi_f$ according to the MS model (equation \[Da\]) multiplied by a factor $2/3$ to correct for the fraction of vertically aligned bonds. The symbols correspond to the $\lambda_{S}$ value (see legend panel (b)) color-coded according to the amount of disorder $\delta\lambda$ (see legend panel (a)). Error bars indicate standard deviation based on 50 configurations.[]{data-label="fig:BrokenBondRescaling"}](figure6.eps)

The first parameter we explore is the fraction of broken sacrificial bonds at final failure $\phi_f$, as we have a direct prediction from the MS model. In Fig. \[fig:BrokenBondRescaling\](a), we plot $\phi_{f}$ as a function of $\Delta\alpha$, the distance from the BDT as introduced in equation \[Da\] for various amount of disorder $\delta\lambda$. Clearly, we see a strong increase in the number of broken bonds upon going from the brittle to the ductile regime. For $\delta\lambda=0$ we observe that all data points collapse on a step function exactly at the BDT, $\Delta\alpha=0$. When disorder is introduced, the transition becomes more gradual, corresponding to the shift of the B-to-QB and D1-to-D2 transition in Fig. \[fig:PhaseDiagram\]. We can understand these effects by pointing out that both weak and strong bonds are introduced due to the spread in sacrificial bond strength. For $\Delta\alpha<0$, $\phi_{f}$ is increased by the disorder due to the presence of weak springs that break before brittle failure occurs. Similarly, for $\Delta\alpha>0$, $\phi_{f}$ decreases due to disorder that introduces bonds that are too strong to break before the matrix fails.

According to the prediction for the MS model in equation \[dist2\] we can collapse all this data on one master curve for all values of $\lambda_S$, $\mu_M/\mu_S$ and $\delta\lambda$, if we consider the rescaled distance to the BDT: $\Delta\tilde{\alpha}$. Indeed we obtain a good collapse for $\phi_{f}$, as shown in Fig. \[fig:BrokenBondRescaling\](b). However, in the brittle regime $\phi_f$ is significantly higher than the prediction of the MS model (dashed line). This could be an indication that the load sharing via the network enhances the failure of sacrificial bonds, however we can not exclude that the discrepancy is due to a finite size effect.

### \[sec:stressConc\] Stress concentration in the spring network model

![image](figure7.eps)

In the MS model it is assumed that stress is distributed homogeneously throughout the entire system, even after failure of sacrificial bonds. Hence, stress can never concentrate within the MS model. In the RSN model stress concentration is possible and therefore it can be used to study the nucleation and propagation of failure in double networks. It is fascinating that this fundamental difference between the MS model and the RSN model does not result in big differences in macroscopic failure or the fraction of broken sacrificial bonds.

Interestingly, a second look at the type of sacrificial bonds that break does reveal a striking difference between the MS model and the RSN model caused by the concentration of stress. By construction, in the MS model sacrificial bonds break in the order of their strength, from weak to strong, until a certain threshold $\lambda_{S}^{*}$ is reached and macroscopic failure occurs (see Fig. \[fig:StressConc\](a)). At this threshold, $\lambda_{S}^{*}$, the force required to break a (matrix-reinforced) sacrificial bond is equal to the force required to break a matrix bond. However, in the RSN model we see that sacrificial bonds stronger than $\lambda_{S}^{*}$ also break (see Fig. \[fig:StressConc\](b) and (c)). In particular, in Fig. \[fig:StressConc\](d) we plot the fraction of broken strong sacrificial bonds $\phi_{strong,f}$ (i.e. bonds with $\lambda_S > \lambda^{*}_{S}$ that break) as a function of the distance from the BDT $\Delta\tilde{\alpha}$ for various amounts of disorder $\delta \lambda$. At low $\delta\lambda$ and far below the BDT, failure of strong sacrificial bonds is required in order to achieve macroscopic failure, since there are just not enough weak sacrificial bonds (Fig. \[fig:StressConc\](b)). However, upon approaching the transition at $\Delta\tilde{\alpha}=0$, $\phi_{strong,f}$ grows, even though a decrease could be expected based on the total concentration of strong sacrificial bonds (Fig. \[fig:StressConc\](c)). Interestingly, we observe a maximum in $\phi_{strong,f}$ before fully entering the ductile regime $\Delta\tilde{\alpha}>0$ where mostly weak sacrificial bonds control the fracture process. Strong bonds can only break if the force they carry exceeds the maximum force that is expected based on the macroscopic force balance. Thus, we conclude that the network structure plays a crucial role in (re)distributing the load during deformation, giving rise to stress concentration. The development of stress heterogeneity is tightly bound to disorder that is present in the initial structure. We expect that by including more disorder in the network (e.g. with inhomogeneities in the connectivity of the nodes or in the stiffness of the bonds) the stress heterogeneity would be even larger [@Dussi2020].

Because in the RSN the initial structure is homogeneous, stress concentration does not occur from the onset of deformation, but develops during deformation as bonds are ruptured. We can quantify this development via the difference between the failure threshold of a broken sacrificial bond $\lambda_{S,fail}$, a measure for the actual local extension, and the extension estimated from the global strain assuming all deformations are affine, $dl_{aff}$. We therefore define $\Lambda=\lambda_{S,fail}-dl_{aff}$. To compare with the MS model, it must be noted that in this model the stress is homogeneous, but the strain field is not. In fact, the microcracks take up additional strain to compensate for the absence of the sacrificial network. This means that the deformation of the remaining sacrificial bonds is less than the affine strain, so that $\Lambda$ is always smaller that 1 in the MS model (See dashed line in Fig. \[fig:StressConc\](e)). If at any time $\Lambda>0$, stress concentration must therefore be present.

This analysis reveals different modes of stress concentration during the failure of a RSN system as shown in in Fig. \[fig:StressConc\](e). At first nothing breaks, this regime coincides with the reversibly elastic regime found at the macroscopic level (regime A). Then, sacrificial bonds start to break exactly as predicted by the affine deformation. Here both the average stress concentration, $\langle\Lambda\rangle$, and its variance are zero (regime B). Upon further increasing the deformation, the variance of $\Lambda$ starts to grow indicating stress heterogeneity in the system, quickly followed by an increase of $\langle\Lambda\rangle$ indicating stress concentration (regime C). As a result strong sacrificial bonds also start to break and failure of weak sacrificial bonds is postponed or even prevented (see Supplemental Material [@supp]). Finally we arrive at a peak in $\langle\Lambda\rangle$, after which the $\langle\Lambda\rangle$ decreases until macroscopic failure is reached (regime D). We attribute this decrease to structural relaxation of the sacrificial network via the formation of microcracks, leading to loss of rigidity and eventually loss of percolation in the sacrificial network. It must be noted that the variance in $\langle\Lambda\rangle$ remains constant, indicating that heterogeneity in stress is present until the final failure event.

### Crack development and stress delocalization

![image](figure8.eps)

The development of stress heterogeneity and stress concentration during deformation, is the result of microscopic damage evolution in the sacrificial network: the nucleation and propagation of cracks. We also suggested that the development of stress concentration is tightly linked to the disorder in the initial, undeformed, system. How crack nucleation and propagation are influenced by both initial disorder and stress heterogeneity during deformation is not trivial. We have already seen that in the RSN model not all sacrificial bonds that break are part of the largest crack (Fig. \[fig:StressStrain\]). Thus, it seems that stress can be delocalized away from the crack tip of the largest crack throughout the system via diffuse damage. In this section, we quantify the evolution of (micro)cracks in the sacrificial network, focusing on diffuse damage and its coupling with stress delocalization.

In Fig. \[fig:Diffuse\](a), we plot the diffuse damage, quantified by the fraction of broken sacrificial bonds (with respect to the total number of initial sacrificial bonds) that are not part of the largest crack $\phi_{diff}$ as a function of strain $\epsilon$ for three representative systems below, around and above the BDT (see Supplemental Material [@supp]). In the brittle case ($\Delta\tilde{\alpha}=-2.9$), $\phi_{diff}$ barely increases and fracture occurs at smaller strain. Closer to the BDT ($\Delta\tilde{\alpha}=-0.7$), the diffuse damage significantly increases, reaches a maximum, and just before final failure slightly decreases. This decrease in the curve indicates that microcracks merge into the largest crack that eventually breaks the system (Fig. \[fig:Diffuse\](a)). Beyond the BDT ($\Delta\tilde{\alpha}=1.1$), the diffuse damage reaches a higher maximum of more than 15% of the total amount of sacrificial bonds before a significant decrease that extends for a large strain interval. Finally, $\phi_{diff}$ reaches zero, indicating that all microcracks are now merged in a single large crack or damage zone, well before the final failure and all subsequent failing sacrificial bonds join this large crack.

With respect to the global stress-strain response, diffuse failure is expected to be most influential before the BDT. In the ductile regime, broken bonds will eventually join to form one large damage zone. However, in the brittle regime, many bonds can break in a diffuse fashion without being part of the largest crack, potentially diminishing stress concentration. We therefore focus on the final ratio of the broken sacrificial bonds belonging to the diffuse damage $\phi_{diff,f}$ as an indicator of stress delocalization and plot it versus $\Delta\tilde{\alpha}$ for various amounts of disorder in Fig. \[fig:Diffuse\](b). We observe a non-monotonic behavior with a peak slightly before the BDT. At the BDT $\phi_{diff,f}$ drops drastically as all bonds eventually join the largest crack. We conclude that just before the BDT, the diffuse damage is maximum and therefore hypothesize that the crack nucleation is delayed the most in this region. A consequence of this delayed nucleation is the widening and branching of the largest crack as shown in Fig. \[fig:Diffuse\] (a), snapshot in the middle. The widening of the damage zone in the sacrificial network around a defect is often used as an explanation for the enhanced fracture toughness in double network materials[@Brown2007; @Gong2010]. Although our simulations do not measure the fracture toughness directly, our results suggest that widening of the damage zone and therefore fracture toughness is maximal just before the BDT.

\[sec:MatProp\] Relating microscopic events to the macroscopic failure regimes
------------------------------------------------------------------------------

![Characteristic strains corresponding to microscopic and macroscopic events as a function of $\Delta\tilde{\alpha}$ ($\langle<\lambda_S\rangle=1.4$ and $\delta\lambda=0.250$) . Strain at maximum stress concentration $\epsilon_{\Lambda,max}$ (purple triangle), strain $\epsilon_{diff,max}$ at which the diffuse damage is maximal (pink circle) and strain $\epsilon_{perc}$ at which percolation in the sacrificial network is lost (gray square). For the macroscopic properties we show the softening strain $\epsilon_{soft}$ (light gray downward triangle) and the failure strain $\epsilon_{f}$ (dark gray star). The softening strain is defined as the onset of softening in the stress-strain curve. (see Supplemental Material for details [@supp]). The vertical lines indicate the transitions between the failure regimes estimated from the phase diagram in Fig. \[fig:PhaseDiagram\](b).[]{data-label="fig:MacroToMicro"}](figure9.eps)

The network and the heterogeneous stress distribution emerging upon deformation, clearly have a huge impact on the microscopic failure mechanism. In particular, we have seen that which microscopic failure mechanisms occur is controlled by the distance from the BDT. Here we will unify our microscopic insights with the macroscopic failure regimes described earlier. These results are summarized in Fig. \[fig:MacroToMicro\].

*B regime -* In the brittle regime, far below the BDT, the system deforms homogeneously and every microcrack immediately develops into a macroscopic crack, leading to global failure. Therefore all failure events, both macroscopically and microscopically, occur at the failure strain $\epsilon_f$.

*QB regime -* Approaching the BDT from below, we arrive at the QB regime, where global failure is preceded by softening in the stress-strain curve. The onset of this softening behaviour $\epsilon_{soft}$ becomes distinguishable from the failure strain $\epsilon_f$ around the transition from B-to-QB. Interestingly, the softening strain does not depend strongly on the distance from the BDT, $\Delta\tilde{\alpha}$ (Fig. \[fig:MacroToMicro\]). If we direct our attention to the microscopic level, we see that the onset of softening coincides with the peak in stress concentration $\Lambda$ (Fig. \[fig:StressConc\](e)). Also for $\epsilon_{\Lambda,max}$ we only see a weak dependence on the distance from the BDT (Fig. \[fig:MacroToMicro\]). Here we would like to stress that even though the softening strain is the first clear macroscopic sign of failure, it is not the strain at which the first bonds break. Actually, a significant amount of bonds break before the softening strain and a large heterogeneity develops in the stress distribution, as we can see from the $\Lambda$ curves. The question that remains is: what happens at the softening strain? We hypothesize that at the softening strain a shift occurs in how additional deformation is handled by the network. Before $\epsilon_{soft}$ bonds break homogeneously throughout the sample. These are weak bonds introduced via the spread in $\lambda_S$. With an increase in deformation, stress concentrates around these microcracks, leading to stress heterogeneity. Because the network is largely intact around these microcracks, the microcracks can not expand significantly to compensate for the applied deformation. We think that only after the softening point enough bonds are broken, such that a large part of the added deformation can be absorbed completely by structural relaxation. After $\epsilon_{soft}$ enough stress is concentrated around the microcracks to also cause the failure of strong sacrificial bonds. Still, the failure of these bonds mostly induces new microcracks or causes a minor expansion of existing microcracks as indicated by the strong rise of diffuse failure (Fig. \[fig:Diffuse\]). However, just before final failure, most of these microcracks merge in one large macroscopic crack, leading to global failure as is visible in the snapshot in (Fig. \[fig:Diffuse\]).

*D1 regime -* It is only after the BDT, that the merging of microcracks starts to occur well before system failure, as can be identified from a peak at $\epsilon_{diff,max}$ in the diffuse failure as a function of strain (see Fig. \[fig:Diffuse\](a)). Once we enter the ductile regime, the merging of cracks always occurs at the same strain and therefore seems to be independent of the distance from the BDT. Even though microscopic cracks in the sacrificial network start to merge, this does not mean that the sacrificial network looses percolation. In fact, in the D1 regime for a system with sufficient disorder in $\lambda_S$ loss of percolation in the sacrificial network $\epsilon_{perc}$ only occurs around the failure strain (Fig. \[fig:MacroToMicro\]). Thus, in a system with a network structure such as the RSN, the criterion for entering a plateau region does not have to be loss of percolation in the sacrificial network, but rather the opportunity for structural relaxation within the sacrificial network. This observation provides an alternative to failure mechanisms proposed in literature, where loss of percolation was thought to be necessary [@Sumiyoshi2006; @Nakajima2013a].

*D2 regime -* In the D2 regime, the stress-strain curve enters a matrix dominated regime after the plateau. Microscopically this transition is marked by the loss of percolation in the sacrificial network. From the loss of percolation onward, the stress in the material is carried by the matrix as is also apparent from the constant slope in the stress-strain curve (Fig. \[fig:StressStrain\](a)). However, we do also observe that $\epsilon_f$ continues to grow with increasing $\Delta\tilde{\alpha}$, until a plateau is reached as we have observed earlier in Fig. \[fig:FailureStrain\]. The initial increase of $\epsilon_f$ in the D2 regime, shows that although percolation is lost in the sacrificial network and the matrix carries most of the stress, the sacrificial network can still contribute to the mechanical properties. Finally, we emphasize that during all these processes the matrix network remains largely intact, only very close to to macroscopic failure event the matrix bonds start to break (see Supplemental Material [@supp]).

\[sec:SummConclusions\]Summary and Conclusions
==============================================

In this paper we studied both the macroscopic and microscopic failure behavior of double network materials using a spring network model (RSN) with emergent load sharing. By comparing the results of the RSN with a model based on equal load sharing, we reveal that: (i) The location of the BDT, defined on the basis of macroscopic stress-strain behavior, is captured by a simple force balance. (ii) Disorder introduces intermediate failure regimes but it can be incorporated in the parameter $\Delta\tilde{\alpha}$ to correctly describe the distance from the BDT, allowing for rescaling of the number of total broken sacrificial bonds. (iii) At the microscopic level stress concentration and delocalization reveal a markedly different picture compared to global load sharing.

The overall picture that emerges from the RSN is that the force balance, a central feature of double networks, has significant control over both the macroscopic and microscopic failure behaviour, irrespective of how stress is (re)distributed. By contrast, the nucleation and propagation of (micro)cracks is also highly dependent on the mode of stress (re)distribution. In particular we have identified how stress concentration, diffuse damage and loss of percolation are related to the transitions from B-to-QB, QB-to-D1, and D1-to-D2, respectively.

We highlight that, because many microcracks can form before global failure due to the stabilization by the matrix, the load sharing in double networks becomes highly non-linear due to the interaction between these microcracks. Therefore double networks provide a unique opportunity to exploit these non-linearities in microcrack interaction. For example, this knowledge about the microscopic failure process could aid the development of robust self-healing double networks, as the healing of diffuse damage is easier than the healing of macroscopic cracks [@Jia2016; @Liu2017; @Yan2017].

By extending the RSN model, additional features of double network gel failure could be studied, such as the influence of pre-stress or structural disorder in the sacrificial network. Furthermore the introduction of disorder in the matrix would allow to also study the enhancement of strength and toughness, including the role of microscopic failure. In conclusion, this paper demonstrates that a random spring network model provides the opportunity to systematically study the microscopic failure process within double networks.

The work of J.T., S.D. and J.v.d.G. is part of the SOFTBREAK project funded by the European Research Council (ERC Consolidator Grant).
