Network evolution controlling strain-induced damage and self-healing of elastomers with dynamic bonds

The equilibrated DPN model is loaded in uniaxial tension along the $x$-direction at a strain rate of $8.7\times 10^{5}\,{\rm s}^{-1}$ until $400\%$ strain, followed by unloading (i.e. reversing the sign of strain rate) until zero stress is reached. The volume is conserved during the loading-unloading process. The time period of loading from $0$ to $400\%$ strain is $~{}\sim 4.6\,\mu{\rm s}$. The residual strain at the end of unloading is about $100\%$. As shown in Fig. 1(b), the stress-strain curves during loading and unloading exhibit a hysteresis which is commonly referred to as the Mullins effect [20]. Although the Mullins effect is commonly applied to the behavior exhibited by filled polymers (polymer networks with filer particles, the phenomenon is common to all polymers, including polymers without fillers. The stress-strain curve during unloading is lower than that during loading, indicating that certain bond-breaking events have happened during loading that weakened the polymer network.

We made three copies of the configuration obtained at the end of the loading-unloading cycle and subjected them to different loading paths. The first configuration is immediately loaded in tension along the $x$-direction again. The resulting stress-strain curve is shown in Fig. 1(b) (red line). The second configuration is immediately loaded in tension along the $y$-direction. The resulting stress-strain curve is shown in Fig. 1(c) (red line). The differences in these two stress-strain curves show that the material becomes anisotropic after the first loading-unloading cycle. It is also interesting to note that the stress-strain curve along the second loading in the $y$-direction is lower than the initial stress-strain curve. This means that the first loading-unloading cycle along $x$ also weakens the material in the $y$-direction.

The third configuration is allowed to relax under zero stress for a long period of $107.4\,\mu$s. Remarkably, during this period the strain in the $x$-direction is seen to go back to close to zero ($\sim 2\%$ strain), and the simulation cell eventually returns close to its initial cubic shape. Furthermore, when the resulting configuration is subjected to uniaxial loading again in $x$ and $y$ directions, the stress-strain curves recover the original stress-strain curve, as shown in Fig. 1(b) and (c), respectively (blue lines). This means that our CGMD model successfully captures the self-healing behavior of DPNs. After the self-healing period, the elastomer not only returns close to its initial shape but also regains its original mechanical properties. In this work by self-healing we mean the process of DPN recovering its original shape and stress-strain properties after the strain-induced damage by the first loading-unloading cycle, not the process of rejoining two halves of the elastomer after being cut by a blade.

To make a qualitative comparison of the kinetics of the self-healing process with experiments, we plot the fraction residual strain $\bar{\epsilon}=\epsilon/\epsilon_{\rm s}$ as a function of dimensionless time $\bar{\tau}=t/\tau$, where $\epsilon_{\rm s}$ is the residual strain at the beginning of the self-healing period, and $\tau$ is a characteristic time obtained by fitting the strain $\epsilon(t)$ to an exponential function ($\epsilon_{\rm s}\exp(-t/\tau)$). Although due to computational limitations, the absolute time scale of self-healing in the CGMD simulations is much faster than that in the experiments, the non-dimensionalized relaxation curve $\bar{\epsilon}(\bar{\tau})$ predicted by CGMD shows qualitative agreement with some of the experiments, as shown in Fig. 2(a). In particular, the relaxation is characterized by a rapid period where a large fraction (e.g. 60%) of the strain is recovered over a short period of time (e.g. $\bar{\tau}\in[0,0.1]$), followed by a much longer period where the remaining strain is recovered. We also extract the configurations from the self-healing process after it has proceeded for different amounts of time, and subject these configurations to uniaxial tensile loading in $x$-direction again. Fig. 2(b) plots the corresponding stress-strain curves; they are qualitatively similar to the experimental results [21, 22, 23] shown in the inset. In particular, the initial slope and the general shape of this family of stress-strain curves are similar to each other. The longer the healing time, the higher is the peak stress reached at re-loading (at the final strain of $400\%$).

We have seen that the CGMD model with reversible (A-B-A type) cross-links is able to capture both the stress-strain hysteresis and the self-healing behavior of DPNs. This CGMD model thus provides us with an opportunity to understand the physical origin of these behaviors. We start with three general remarks on how we may approach this problem. First, we can establish that these behaviors are caused by the bond-breaking and re-forming events that alter the network topology, because these behaviors disappear when we change the cross-links to be permanent (i.e. unbreakable) resulting in the absence of bond-breaking events in the CGMD model [24].

Second, many existing theories that link the stress-strain responses of polymer networks to bond-breaking events are based on the length of local chains between two neighboring cross-link sites on the same backbone. For example, the network alteration theory [25, 26] attributes the changes in the stress-strain curve to the changes in the distribution of local chain lengths, with more recent works extending this analysis to polymers with dynamic bonds [10, 27]. However, it can be easily shown that the local chain length cannot possibly explain the behaviors reported here. In our CGMD simulations, most of the reversible bonds immediately form (with other partners) after breaking, so that at any instant, $>99\%$ of all Type-A (and Type-B) beads remain bonded. Since the Type-A beads are regularly distributed on the polymer backbone (with a fixed period of 15 beads), the local chain length remains very close to 15 during the entire process of loading, unloading, and self-healing. Therefore, it can be concluded that the rich behaviors exhibited by DPNs are not related to the local chain length (which remains at 15) but are the result of the changing network topology. In other words, although every Type-A bead stays connected (through the A-B-A cross-link) to another Type-A bead nearly all the time, what is important is which ones are connected to each other. This connectivity information can be specified by a connectivity matrix $\mathbf{A}$, which is an $n\times n$ matrix, with $n$ being the number of Type-A beads; $A_{ij}=1$ if bead $i$ and bead $j$ are cross-linked and $A_{ij}=0$ if they are not. We can thus envision the scenario that during deformation (or healing), the $\mathbf{A}$ matrix evolves with strain (or time), and hence changes the stress-strain characteristics of the DPN. While we believe that the connectivity matrix formally captures the network characteristics responsible for the stress-strain hysteresis and self-healing behaviors, it is a microscopic description of the network that still contains many irrelevant details. For example, another realization of the CGMD model of the DPN with the same general procedure but different initial configurations (e.g. using different random seeds) would result in a different connectivity matrix but a DPN model with nearly identical deformation and self-healing behaviors. Therefore, we still need to pinpoint the essential microstructural features contained in the connectivity matrix that are responsible for the deformation and self-healing behaviors of DPNs.

Third, a recent study has shown that the distribution of shortest-path lengths between far-away beads in the polymer network is the controlling microstructural parameter responsible for the stress-strain hysteresis of elastomer with breakable and irreversible cross-links [24]. We note that the shortest-path (SP) lengths between far-away beads are a non-local property of the polymer network; hence an SP-based theory will not suffer from the same limitation as the theories based on local chain lengths mentioned above. Therefore, it is of interest to see if SP lengths also control the deformation and self-healing behaviors of DPNs reported here. However, the reversible bonds in DPNs make the situation much more complex than traditional elastomers with irreversible bonds. In the language of the connectivity matrix, all that can happen for an elastomer with irreversible bonds is that certain matrix element $A_{ij}$ can change from $1$ to $0$ (and stays $0$ afterwards). However, in DPN a matrix element $A_{ij}$ can also change from $0$ to $1$. In fact, for the choice of the interaction potential in our DPN model, the total number of non-zero elements stays nearly the same all the time.

To reduce the computational cost without the loss of information about the network connectivity, we define a graph that contains only Type-A beads as vertices and weighted edges between them. Two vertices that are connected by a cross-link result in an edge with weight $1$. Two vertices that are connected by a polymer backbone result in an edge with a weight equal to the length of the polymer chain between them (which is always $15$ in our model). We then choose an offset vector $\vec{q}$. For each vertex (e.g. A), we identify the vertex (e.g. A’) that is nearest to the location offset from A by $\vec{q}$ and then find the shortest path (SP), which is the contour length traveled over all the bonds between them. This shortest path search is performed for all vertices for the same offset vector $\vec{q}$ [24]. When the polymer network is deformed, the offset vector $\vec{q}$ is subjected to the affine transformation consistent with the macroscopic strain. In this work, we focus on the cases where $\vec{q}$ is a repeat vector of the simulation box, $\vec{L}_{x}$ in the $x$-direction or $\vec{L}_{y}$ in the $y$-direction, as illustrated in Fig. 3. For this specific choice of $\vec{q}$, vertex A’ will always be the periodic image of vertex A in $x$-direction or $y$-direction, respectively. According to this definition, the shortest paths would not change with deformation if there were no bond-breaking or bond-forming events that would cause a change in the network connectivity.

Fig. 4(a) shows the evolution of the SP length distribution of beads offset by $\vec{q}=\vec{L}_{x}$ during the deformation shown in Fig. 1(b). Before initial loading, the SP length follows a Gaussian-like distribution, reflecting the randomness of the initial polymer structure. During the initial loading, this histogram noticeably shifts to the right (towards longer SPs) due to a large number of bond-breaking events, which destroy the original SPs, and replace them with longer SPs. The SP distribution changes negligibly during unloading. When the polymer is re-loaded in $x$ immediately after unloading, the SP length distribution also stays nearly unchanged. That the SP length distribution (for $\vec{q}$ along the loading direction) increases during loading and remains nearly unchanged during unloading and immediate reloading is similar to the observation in polymer networks with irreversible bonds [24].

Figure 4(b) shows that loading in $x$ causes the SP length distribution in $y$ to shift to the left (for $\vec{q}=\vec{L}_{y}$), i.e., the SPs in $y$ become shorter. This is different from polymer networks with irreversible bonds where the SP lengths in the transverse direction remain nearly unchanged during loading [24]. We attribute this unique behavior of the anisotropy in the SP evolution observed in the DPNs to the reversible bonds, which are able to create shorter SPs in $y$ while the cell size in $y$ decreases with strain in $x$.

Remarkably, the SP length distributions in both $x$ and $y$ directions return to the initial pattern (i.e., a pronounced shift to the left) during the healing process. During subsequent reloading, the SP length distributions change with strain in the same way as during the initial loading. This indicates that the damage and anisotropy caused by the initial loading are eliminated by the healing process. The strong correlation between SP lengths and stress-strain hysteresis during loading, unloading, and self-healing leads us to conclude that the SP is a key microstructural feature controlling the mechanical properties of DPNs.

The stress-strain behavior and bond dynamics of DPN are more complicated than of polymer networks with irreversible bonds. For comparison purposes, let us first give a quick summary of the situation for polymer networks with irreversible bonds. There, the stress-strain curve remains reversible and SP distribution does not change during unloading and reloading as long as the previous maximum strain is not exceeded, because no new bonds break and no new bonds can form [28]. Once the previous maximum strain is exceeded, bond-breaking events occur resulting in an irreversible softening of the material as is the characteristic of the Mullins effect [24]. Hence the entire stress-strain behavior for loading/unloading along a given direction can be described by a family of curves parameterized by the maximum strain. This relatively simple picture no longer applies to DPNs due to their ability to form new bonds. To probe the behavior of DPNs, we performed CGMD simulations and the SP analysis along different loading/unloading paths.

Firstly, we analyze the shortest path evolution during the loading trajectory shown in Fig. 1(b-c). Fig. 5(a) shows the average SP (in both $x$ and $y$ directions) as a function of strain during loading to $400\%$, unloading, healing and reloading. The SP-strain map shown here provides a more concise summary of the SP evolution shown in Fig. 4. It is interesting to note that on the SP-strain map the healing process can be clearly seen to occur in two steps. The first step corresponds to a (rapid) reduction of strain with nearly no change of SP in $x$, and the second step corresponds to a (slow) reduction of both SP in $x$ and strain. The first, rapid, step of healing can also be seen in Fig. 2. Fig. 5(b) shows the average SP in $x$-direction during loading after the DPN has been allowed to heal for different amounts of time $\tau_{H}$. Interestingly, it appears that they can be well approximated by a family of linear curves with different initial points but a common final point (at $400\%$ strain).

Secondly, we consider a loading path where the sample is repeatedly unloaded to zero stress and then re-loaded to a higher strain. Fig. 6(a) shows the stress-strain curves where significant hysteresis can be observed, similar to the experimental observation (inset). Fig. 6(b) shows the corresponding map of SP (in $x$) versus strain. It can be seen that during immediate reloading the SP appears more or less constant until the previous maximum strain is reached, beyond which the SP increases with strain at almost the same rate as during the initial loading. Also shown in Fig. 6(b) are SP-strain curves during healing from different unloaded states.

Finally, we consider loading paths where the sample is stretched to different extents of maximum strain, followed by unloading and healing. Fig. 7(a) shows the stress-strain curves and Fig. 7(b) shows the corresponding map of SP (in $x$) versus strain.

Although the evolution of SP depends on the loading path in a complex way, there are some general rules that can be stated based on the observations from the CGMD simulations. These empirical rules can serve as a basis for the future development of a comprehensive theory on the microstructural evolution of DPNs. For simplicity, we limit the scope of our discussion to uniaxial tensile loading along $x$-direction while keeping the volume constant.