Most tough hydrogels suffer accumulated damages under cyclic loads. The damages may stem from breakage of covalent bonds, unzipping of ionic crosslinks, or desorption of polymer chains from nanoparticle surfaces. Recent experiments report that when a tough hydrogel is subject to cyclic loads, the stress–stretch curves of tough hydrogels change cycle by cycle and approach a steady-state after thousands of cycles, denoted as the shakedown phenomenon. In this paper, we develop a phenomenological model to describe the shakedown of tough hydrogels under prolonged cyclic loads for the first time. We specify a new evolution law of damage variable in multiple cycles, motivated by the experimental observations. We synthesize nanocomposite hydrogels and conduct the cyclic tests. Our model fits the experimental data remarkably well, including the features of Mullins effect, residual stretch and shakedown. Our model is capable of predicting the stress–stretch behavior of subsequent thousands of cycles by using the fitting parameters from the first and second cycle. We further apply the model to polyacrylamide (PAAM)/poly(2-acrylanmido-2-methyl-1-propanesulfonic acid) (PAMPS) and PAAM/alginate double-network hydrogels. Good agreement between theoretical prediction and experimental data is also achieved. The model is hoped to serve as a tool to probe the complex nature of tough hydrogels, through cyclic loads.

## Introduction

Hydrogels, aggregating a three-dimensional (3D) polymer network and water molecules, are abundant in plants and animal tissues. Traditional hydrogels like jello and tofu are fragile with very low fracture toughness. Over recent decades, various strategies have been proposed to design hydrogels with extraordinary mechanical properties. Examples include double-network hydrogels [1–3] and nanocomposite hydrogels [4–6]. The fracture toughness has been increased from ∼10 J/m^{2} to ∼10,000 J/m^{2}. The discovery of tough hydrogels greatly expanded the application scope of hydrogels in the field of stretchable ionotronics [7–12], tissue engineering [13–15], soft robotics [16–21], and tough adhesives [22,23].

Many hydrogel devices in use are prone to cyclic loads. Recent papers reported the mechanical behavior, especially fatigue, of tough hydrogels under cyclic loads [24–31]. Three important features of the stress–stretch curves of tough hydrogels are demonstrated in Fig. 1. First, Mullins effect is clearly observed. In particular, the stress level in the second cycle is much lower than that in the first cycle. The stress–stretch hysteresis in the first cycle is the largest, indicating that most energy dissipation takes place during the initial loading and unloading. Second, a finite residual stretch is usually observed after each loading cycle. The residual stretch depends on the loading history, such as the strain rate, maximum stretch and number of cycles. Third, during the subsequent cycles, the sacrificial bonds in the hydrogel keep breaking, and lead to a gradual decrease of the maximum stress over cycles. After thousands of cycles, the change of stress–stretch curve becomes negligible and the hydrogel approaches a steady-state.

The mechanical models to describe Mullins effect and residual stretch of tough hydrogels have been developed recently, referring to models of Mullins effect of elastomers [32–34]. For example, Wang and Hong analyzed Mullins effect of tough hydrogels by introducing a damage evolution law [35]. Zhang et al. developed a model of Mullins effect to predict fracture energies and crack-tip fields of tough hydrogels [36]. Lu et al. developed a model for tough hydrogels by coupling viscoelasticity and Mullins effect [37]. Tang et al. developed a model to describe the residual stretch of nanocomposite hydrogels [38]. However, all the existing models only characterize the deformation behaviors in the first few loading cycles, but do not account for the gradual decrease of maximum stress induced by thousands of loading cycles.

In continuum mechanics, shakedown represents the steady-state that the material is either pure elastic in a closed elastic-plastic loop or of no further increase of hysteresis [39–42]. The shakedown of metals results from the plastic flow. When the magnitude of cyclic stress is beyond the yielding stress of a metal, three different responses can be observed: elastic shakedown, plastic shakedown and ratcheting. The understanding of shakedown of metals has greatly aided their applications, such as the design of nuclear reactor fuel elements and thin film structures [43–45]. The stress–stretch curves of tough hydrogels change cycle by cycle and approach a steady-state after thousands of cycles, which was also named as shakedown [24,26]. For example, the polyacrylamide (PAAM)/poly(2-acrylanmido-2-methyl-1-propanesulfonic acid) (PAMPS) double-network hydrogel consists a tightly crosslinked PAMPS network as sacrificial bonds, and a loosely crosslinked PAAM network that provides elasticity [1]. Upon stretching, the damage of the PAMPS network leads to large energy dissipation and significant stress–stretch hysteresis [46]. The hysteresis and stress level appear high in the first cycle of loading, but continuously decrease with thousands of cycles and reach a steady-state [26]. The steady-state of shakedown is potentially important for the concept of *intrinsic fracture toughness* [47], which refers to the minimum energy taken to propagate a crack per unit area. This enables us to remove the energy dissipation from the bulk material and study the fracture of tough hydrogels with respect to their elastic networks when the materials reach a steady-state under cyclic loads [24,26]. While the shakedown of metals has been extensively studied, this shakedown of tough hydrogels still lacks theoretical understanding, especially with respect to mechanical modeling.

In this paper, we develop a model to describe shakedown of tough hydrogels for the first time. Motivated from the experimental observations, we make a distinct treatment between the first loading cycle and the subsequent cycles in the model formulation, and incorporate the cycle number into the damage variable to describe the gradual decrease of the maximum stress. We compare our theory to the experimental data of three types of tough hydrogels, including nanocomposite hydrogels, PAAM/PAMPS hydrogels and PAAM/alginate hydrogels. The model fits all the experimental data well.

## Model Formulation

### Mullins Effect in Multiple Cycles.

Two types of theoretical models have been developed over the decades to describe Mullins effect. One is the phenomenological model and the other is the physical model. The physical models provide various possible mechanisms of Mullins effect, such as bond rupture, molecules slipping, filler rupture, and disentanglement [33]. The physical models can capture the main feature of Mullins effect, but usually fail to fit the experimental data well due to the much more complex nature in real materials. The phenomenological models usually introduce an additional parameter, the softening variable *η* or the damage variable (1 − *η*), to the function of the free energy density. Mullins effect is characterized by specifying an evolution law of *η*. Existing phenomenological models for Mullins effect mainly focus on the first loading cycle. Dorfmann and Ogden developed a model describing the Mullins effect in the first and the second loading cycle [48]. In this work, we develop a phenomenological model to fit Mullins effect in tough hydrogels particularly under prolonged cyclic loads.

Because of incompressibility, the third principle stretch is calculated as $\lambda 3=(\lambda 1\lambda 2)\u22121$. The nominal stresses are $s1\u2212(s3/\lambda 12\lambda 2)=(\u2202W(\lambda 1,\lambda 2,\eta )/\u2202\lambda 1)$ and $s2\u2212(s3/\lambda 1\lambda 22)=(\u2202W(\lambda 1,\lambda 2,\eta )/\u2202\lambda 2)$.

*η*does not change in the loading process, but evolves in the unloading process as a function of the principle stretches $(\lambda 1,\lambda 2)$, as well as the maximum stretch of the cycle $(\lambda 1m,\lambda 2m)$. Following the models developed by Ogden and Roxburgh [32] as well as Zhang et al. [36],

*η*evolves with the deformation in the unloading process as

where $erf(\u22c5)$ is the error function, $Wm$ is the maximum strain energy density before unloading, *r*, *c*_{0}, and *β* are constants that characterize the damage of the material.

*η*. The experimental data of the subsequent cycles indicate that the maximum stress evolves with the number of cycles

*N*following a power law approximately

where $smax0$ is the maximum nominal stress in the first loading cycle, *κ* and *γ* are dimensionless parameters that represent the stress decay.

*N*th cycle of unloading, denoted as $\eta uN$, follows the similar form of Eq. (3) with a different parameter

*c*:

_{u}We prescribe that in each subsequent cycle the softening variable follows Eqs. (6) and (5) in the reloading and unloading process, respectively. Under these assumptions, the gradual decrease of maximum stress follows a power law. The fitting parameters are $(r,c0,\beta )$ for the first cycle, $(cr,cu)$ for the subsequent cycles, and $(\kappa ,\gamma )$ for the gradual decrease of maximum stress.

### Residual Deformation.

Under cyclic loads, a tough hydrogel exhibits a residual stretch. The similar residual stretch observed in particle-reinforced rubbers is attributed to the damage-induced anisotropy [49]. For tough hydrogels with reversible weak interactions, usually denoted as *self-recovery hydrogels*, the residual stretch is likely due to the limited time of reformation of the broken bonds during unloading. If the unloading is too fast, and the reformation of bonds is not complete, a residual stretch is generated. This rate-dependent residual stretch has been reported in self-recovery tough hydrogels, and has been characterized theoretically [27,50].

The proportional factor for reformed chains is lumped into the parameter $W2(\lambda 1,\lambda 2)$. In this work, we do not consider the rate-dependent unloading, so that Eq. (7) holds at a fixed unloading rate.

*μ*and

*α*are material parameters. To describe the residual stretch, we choose a suitable form for $W2(\lambda 1,\lambda 2)$ with

*λ*

_{1}and

*λ*

_{2}asymmetric to represent anisotropy, following [49]:

where *ν*_{1}, *ν*_{2}, and *ν*_{3} are material parameters.

When the stress returns to zero, the residual stretch is calculated by setting *s* = 0. When the stretch returns to 1, the first term in Eq. (10) becomes zero, and the second term gives a negative stress.

## Experiments

We use the developed model to fit the experimental data of three kinds of tough hydrogels, including a PAAM/PAMPS hydrogel, a PAAM/alginate hydrogel, and a nanocomposite hydrogel. The data of the PAAM/PAMPS hydrogel and PAAM/alginate hydrogel are taken from Refs. [24,26]. In addition, we synthesize the nanocomposite hydrogel and measure its cyclic loading stress–stretch curves. The synthesizing method is as follows. We purchased from Aladdin (China) the following substances: *N*,*N*-dimethylacrylamide (DMAA), *N*,*N*′-methylenebisacrylamide (MBAA, crosslinker), ammonium persulfate and *N*,*N*,*N*′,*N*′-tetramethylethylenediamine (TEMED). Clay Laponite^{®} XLG was purchased from BYK Additives and Instruments, Germany. Clay Laponite XLG of 1.14 g was dispersed in 30 mL distilled water and stirred for 1 h to form a homogeneous solution. Afterwards, DMAA of 3 ml was added. The solution was stirred for another 0.5 h. TEMED of 30 *μ*L and ammonium persulfate of 0.33 g were added to trigger the polymerization. The precursor was poured into acrylic molds and sealed by two glass plates for 12 h to complete the reactions.

The synthesized hydrogels were cut into samples of rectangular sheet with length of 50 mm and width of 10 mm by a laser cutter (Universal Laser Systems, Scottsdale, AZ). We used Shimadzu tester (Autograph AGS-X, Japan) to conduct the pure shear test on the prepared samples. The cyclic loading, unloading and reloading tests were performed using a constant loading rate of 960 mm/min. The total number of cycles was 10,000. To minimize dehydration of the samples during the test, we used an acrylic chamber to cover the entire sample and the grips. The force and displacement in each cycle were recorded and converted to the stress–stretch curves (Fig. 2). The nominal stress is calculated as the force divided by the cross-sectional area of the sample in the undeformed state. The stretch is calculated as the length of the sample in the deformed state divided by that in the undeformed state.

Figure 2(a)–2(d) show the stress–stretch curves of nanocomposite hydrogels at different maximum stretches. The polymer chains of nanocomposite hydrogels are attached on the clay sheets by physical interactions [4,52]. Upon stretching during the first cycle, the polymer chains detach, leading to the softening of the stress. As the sample undergoes more cycles, more polymer chains detach, causing the continuous softening. Figure 2(e) shows the decay of the maximum stress for each cycle. The hysteresis is significant in the initial cycle, and decays in subsequent cycles. After thousands of cycles, the hysteresis becomes negligible. Figure 2(f) shows the stress–stretch curve of the 2000th cycle for different maximum stretches. Residual stretch is observed in each cycle when the nominal stress returns to zero.

## Results and Discussion

### Comparison Between Theory and Experiment.

We apply the model to fit the stress–stretch curves of nanocomposite hydrogels under cyclic loads (Fig. 3(a)). For the loading curve in the first cycle, the softening variable is *η* = 1. We fit the curve by the two-parameter Ogden model with *μ* = 3.97 kPa and *α* = 2.158. For the unloading curve in the first cycle, the softening variable evolves with (*r*, *c*_{0}, *β*). *r* reflects the damage extent of the material and *c*_{0} represents an energy threshold for activating dissipation. *r* and *c*_{0} determine the shape of the unloading curve. *β* is a positive number to avoid over stiff response at the initiation of unloading from relatively large stretch levels [36,37]. We obtain the fitting parameters *r* = 1.75, *c*_{0} = 10 kJ/m^{3} and *β* = 0.1. To fit the residual stretch, we select a suitable $W2(\lambda 1,\lambda 2)$ with the parameter *ν*_{3} larger than *ν*_{1} to represent anisotropy [49]. We fit the experiment and obtain *ν*_{1} = 0.04 and *ν*_{3} = 0.72. With the model, Mullins effect and the residual stretch are well captured in the first cycle. Based on the evolution of maximum stress under prolonged cyclic loads in the experiment, we obtain *κ* = 0.85 and *γ* = 0.16 by using Eq. (4). The parameters *c _{r}* and

*c*have the similar function with

_{u}*c*

_{0}in fitting the subsequent cycles. We fit the second cycle and obtain parameters

*c*= 45 kJ/m

_{r}^{3}and

*c*= 25 kJ/m

_{u}^{3}. The other parameters in the second cycle are the same as those used in the first cycle. Figure 3(c) shows the fitting result for the second loading–unloading cycle. We obtain the fitting parameters (

*r*,

*c*

_{0},

*β*,

*ν*

_{1},

*ν*

_{3},

*c*,

_{r}*c*) from the first two cycles. Our model is capable of predicting the stress–stretch behavior of subsequent thousands of cycles by using fitting parameters from the first and second cycle. For the subsequent loading cycles, the parameter

_{u}*η*(

*N*) evolves according to Eqs. (5) and (6). Figure 3(e) and 3(g) show the results for the third cycle and the

*N*th cycle. Good agreement is achieved between the model and experimental data. Even for the 8000th cycle, the model still well captures the stress–stretch curve and the residual stretch. Since the hysteresis becomes negligible after thousands of cycles, we only plot the unloading curves to reflect the stress decay and the residual stretch in Fig. 3(e). The gradual decrease of stress is shown in Fig. 3(i). With the same set of parameters, we compare the prediction of the model with the experimental data under a different maximum stretch ratio

*λ*

_{max}= 4.5 in Figs. 3(b), 3(d), 3(f), 3(h), and 3(j).

We extend the current model to two other materials, the PAAM/PAMPS hydrogel and PAAM/alginate hydrogel. For the PAAM/PAMPS hydrogel, the PAMPS chain is tight and short, and the PAAM chain is loose and long. During stretching, the short PAMPS network becomes highly stretched and ruptures much earlier than the long PAAM network. The damage of the short network induces hysteresis. Following the same strategy for modeling the nanocomposite hydrogel, we fit the experimental data of PAAM/PAMPS hydrogels under cyclic loads in Figs. 4(a), 4(c), and 4(e). For the PAAM/alginate hydrogel, the alginate network is ionically crosslinked by calcium ions. During stretching, the ionic crosslinks unzip, resulting in energy dissipation and hysteresis. The comparison between our model and the experimental data is shown in Figs. 4(b), 4(d), and 4(f). The fitting parameters for all the three hydrogels are listed in Table 1. Although the residual stretches for the three hydrogels are different, the parameters $(\nu 1,\nu 3)$ can be the same, indicating that the residual stretch is mainly influenced by the parameters $(r,c0)$ in the evolution law.

### Parametric Study.

The damage of the material is mainly controlled by the parameters (*r*, *c*_{0}, *β*, *c _{r}*,

*c*), which has been extensively studied in literatures [32,36,37]. In this work we focus on the influences of parameters

_{u}*γ*that controls the decrease of maximum stress and (

*ν*

_{1},

*ν*

_{3}) that control the residual stretch. Figure 5(a) shows the stress–stretch curves of the first, second, and 1000th cycle using the fitting parameters for nanocomposite hydrogels. We change the parameter

*γ*that represents stress decay to

*γ*= 0.02 and keep all other parameters fixed. The stress–stretch curves are shown in Fig. 5(b). We can see that a larger

*γ*leads to a faster stress decay over cycles. We change the parameters that determine the residual stretch to

*ν*

_{1}= 0.04 and

*ν*

_{3}= 4.0 in Fig. 5(c). We can see that a larger ratio of

*ν*

_{3}to

*ν*

_{1}causes a larger residual stretch. The extreme condition

*ν*

_{1}=

*ν*

_{3}will eliminate the anisotropy of Eq. (8) and lead to zero residual stretch.

Figure 6 shows the stress–stretch curves of the 2000th cycle under different maximum stretches using the fitting parameters for nanocomposite hydrogels. When the maximum stretch is *λ*_{max} = 2.0, the hysteresis of the material becomes negligible after 2000 cycles. The softening variables *η*(*N*) approaches to a constant and the material shakes down to a steady-state. The material behaves like a purely elastic material. When the maximum stretch is *λ*_{max} = 5.0, the hysteresis is still apparent after 2000 cycles and more cycles are needed to approach to a steady-state.

## Conclusions

We develop a phenomenological model to describe the shakedown of tough hydrogels under prolonged cyclic loads for the first time. We separately treat the evolution law of the softening variable in the first cycle and subsequent cycles. Motivated by the experimental observation of the stress decay, we specify the evolution law in the *N*th cycle by a power-law function of the number of cycles. We synthesize nanocomposite hydrogels and conduct cyclic tests. Our model fits the experimental data well, and captures the three major features of the stress–stretch curves including Mullins effect, residual stretch and shakedown. Besides, our model is capable of predicting the stress–stretch behavior of subsequent thousands of cycles by using fitting parameters from the first and second cycle. We further apply the model to PAAM/PAMPS hydrogels and PAAM/alginate hydrogels. Good agreements are again achieved between theory and experiment. It is hoped that this phenomenological model can motivate further development of physical models to capture the shakedown of tough hydrogels.

## Acknowledgment

T. Lu acknowledges the valuable discussions with Professor Z Suo at Harvard University.

## Funding Data

NSFC (Grant Nos. 11772249, 11702208).

Program for Postdoctoral Innovative Talents (No. BX201700192).

Young Elite Scientist Sponsorship Program by CAST (2016QNRC001).

Key subject “Computational Solid Mechanics” of the China Academy of Engineering Physics.

NSF MRSEC (DMR-14-20570) at Harvard University.