A Review on the Mullins Effect in Tough Elastomers and Gels

Abstract

Tough elastomers and gels have garnered broad research interest due to their wide-ranging potential applications. However, during the loading and unloading cycles, a clear stress softening behavior can be observed in many material systems, which is also named as the Mullins effect. In this work, we aim to provide a complete review of the Mullins effect in soft yet tough materials, specifically focusing on nanocomposite gels, double-network hydrogels, and multi-network elastomers. We first revisit the experimental observations for these soft materials. We then discuss the recent developments of constitutive models, emphasizing novel developments in the damage mechanisms or network representations. Some phenomenological models will also be briefly introduced. Particular attention is then placed on the anisotropic and multiaxial modeling aspects. It is demonstrated that most of the existing models fail to accurately predict the multiaxial data, posing a significant challenge for developing future anisotropic models tailored for tough gels and elastomers.

1 Introduction

Synthetic elastomers and gels are typically brittle, which cannot satisfy the demand for stretchable and tough materials to support mechanical loading. Intensive research has been conducted to improve the mechanical properties of elastomers and gels in the past few decades [1, 2]. Although various strategies have been employed, such as the incorporation of crystalline domains [3, 4], the entanglement effect [5], and physical cross-linked networks [6], the two most commonly used strategies are using fillers and interpenetrating networks. For the former, using carbon black nanocomposites to reinforce natural rubbers has been extensively investigated due to the wide usage of these materials for tires. Nanocomposite fillers have also been introduced into gels to achieve extremely higher stretchability and toughness compared to covalently cross-linked gels [7]. For the latter, the pioneering work of Gong et al. [8] developed extremely tough hydrogels by interpenetrating a short-chain cross-linked network with a long-chain network, named as double-network hydrogels. Similar methods have been employed to fabricate double/multi-network elastomers [9,10,11]. Regardless of double-network hydrogels or multi-network elastomers, the small volume ratio of the first network acts as the filler, while the other networks can be treated as the matrix. Thus, these soft materials are also referred to as molecular soft composites in Millereau et al. [10].

The high toughness originates from the existence of energy dissipation mechanisms in these soft composites [1]. For filled rubbers and nanocomposite gels, polymer chains can attach to the particle surfaces. The detachment of these chains can gradually dissipate the stored elastic energy. For the double-network hydrogels and multi-network elastomers, the chemically cross-linked filler network can be progressively damaged with the rupture of polymer chains to allow energy dissipation. Thus, all these soft materials exhibit a reduced stress response in the reloading process compared to the virgin specimens. This stress softening phenomenon was intensively investigated by Mullins and coworkers for filled rubbers [12,13,14,15] and subsequently named as the Mullins effect. The Mullins effect is manifested as a hysteresis loop during the loading–unloading process, which is mainly characterized in uniaxial deformation conditions [16, 17]. However, several studies have also investigated the biaxial deformation behaviors of filled rubbers and double-network hydrogels [18,19,20], which reveal that the damage behaviors in these materials are inherent anisotropic.

Although damage mechanisms vary depending on the material type, the mechanical responses share similar features. This allows for the application of the same theoretical framework to model the Mullins effect in different soft materials. Developing damage models for filled rubbers is of great interest to the mechanics community due to its widespread engineering applications. Some representative models include the continuum damage model [21], the pseudo-elastic model [22], the two-phase model [23], the network alternation theory [24], and the progressive damage model [18]. Diani et al. [25] provided an excellent review on this topic, while some recent developments on modeling the Mullins effect in filled rubbers can be found in Xiang et al. [26]. In recent years, various groups, including the authors’ own group, have developed novel damage models for emerging soft composites [27,28,29,30,31,32], which have distinct features compared to classical damage models. However, there is still a lack of a comprehensive review that systematically introduces these developments. Thus, this work aims to provide a summary of these new developments for the Mullins effect in emerging soft materials, with a focus on the theoretical modeling aspects.

The remainder of this paper is organized as follows. In Sect. 2, we clarify the main features of the Mullins effect in various tough elastomers and gels, including nanocomposite hydrogels, double-network hydrogels, and multi-network elastomers. In Sect. 3, we summarize some of the recently developed models for simulating the Mullins effect, with emphasis on new damage mechanisms and network representatives. The constitutive modeling of multiaxial and anisotropic features is then discussed in sections that follow. Finally, concluding remarks are presented in Sect. 6.

2 The Mullins Effect in Soft Composites

It has long been desirable to improve the mechanical properties of soft polymers. Two distinct reinforcing strategies have been widely used. The first is to incorporate particle fillers into the polymer matrix, such as filled rubbers and nanocomposite hydrogels. The other is to synthesize multiple polymer networks, which can be referred to as double-network hydrogels and multi-network elastomers. As discussed in Millereau et al. [10], the rigid network in multiple networks essentially takes a similar role as fillers. Thus, these materials are also named as elastomeric molecular composites. Compared to pure polymer matrix, these soft composites have astonishingly enhanced strength and fracture toughness. In exchange, however, they may lose the elastic feature, exemplified by the lower unloading/reloading stress than the virgin loading. Although the details may be different, such a stress softening effect is indeed similar to that observed in filled rubbers and hence also referred to as the Mullins effect. In this section, we shall revisit the main features of the Mullins effect in these emerging tough soft composites based on the experimental observations reported.

2.1 Filled Rubbers

The addition of nanocomposite particle fillers, such as carbon black, has long been recognized to significantly enhance the mechanical properties of natural rubbers. This leads to filled rubbers that are commonly used in industrial products [33,34,35]. Filled rubbers typically demonstrate the Mullins effect, which includes stress softening, residual strain, and a certain recovery effect, as shown in Fig. 1a. Extensive research has been conducted on the Mullins effect in filled rubbers [19, 36,37,38,39,40], and for a more detailed review, readers can refer to Diani et al. [25].

Fig. 1

Mullins effect in soft composites: a filled rubbers; b nanocomposite hydrogels; c double-network hydrogels; d multi-network elastomers

2.2 The Nanocomposite Hydrogels

Hydrogels exhibit remarkable properties such as elasticity, low modulus, swelling/deswelling, shock-absorbing ability, and low sliding friction [41]. As a result, they have found considerable applications in drug delivery [42], actuators for optics and fluidics [43], artificial polymeric materials [44], and particularly in biological tissue engineering [45]. However, traditional hydrogels are relatively soft and unable to sustain high stress. Moreover, they are also relatively brittle due to the lack of internal dissipation mechanisms [2]. This drawback in mechanical strength and toughness often limits their further applications.

In analogy to the reinforcing scheme of filled rubbers, the strength and toughness of hydrogels can also be improved with the incorporation of nanofillers, resulting in nanocomposite (NC) hydrogels [46, 47]. These nanofillers can be particles or layers with high surface-to-volume aspect ratios, such as nanoclays [46], graphenes [48], carbon nanotubes [49], and hydroxyapatite [50]. Due to the high aspect area, these nanofillers can adsorb a large number of polymer chains through physical/chemical interactions to form nanostructures, which dissipate the stored mechanical energy as the chains detach or slide on the fillers, endowing the NC hydrogels with high toughness [51, 52].

Consequently, such dissipative processes in NC hydrogels give rise to the inelastic stress softening in the unloading and reloading path compared to the virgin loading, i.e., the Mullins effect. The features of Mullins effect in NC hydrogels are quite similar to those in filled rubbers, exemplified by a stress softening, residual strain, and recovery effect, as shown in Fig. 1a, b. However, compared to the vast literature on filled rubbers, only a few results concerning the Mullins effect in NC hydrogels have been reported. Rose et al. [53] found that PDMA/silica NC gels dissipated more energy as strain rate and filler concentration increased, which is analogous to that reported in filled rubbers [15]. Tang et al. [54] reported that the dissipative behaviors of NC hydrogels can be tuned by changing the interaction between polymer chains and clay platelet surfaces. The idealized NC hydrogels should show large hysteresis and good recovery ability so that high toughness and repeatability can be ensured. Yang et al. [55] and Li et al. [56] reported that PAAm/Laponite and PAAm/GO gels both exhibit the self-recovery property at room temperature. However, the hysteresis loop of PAAm/GO gels is negligible at small strains. Lin et al. [57] found that PDMA/silica NC gels can fully recover their elastic modulus during the successive loading–unloading cycles without any resting, and such a mechanical recovery was found to be rate-dependent [53]. This is attributed to the reversible physical adsorption and desorption of PDMA chains on silica surfaces, which prevents permanent damage to network structures [58].

2.3 Double-Network Hydrogels

Another approach toward developing stiff and tough hydrogels is the double-network (DN) scheme. Gong et al. [8] first reported the exceptional resistance to fracture in a novel hydrogel composed of double networks. They synthesized a second network of neutral polymer in the presence of a first covalently cross-linked polyelectrolyte network. The secondary network, with a much lower cross-linker density and higher volume fraction orders than the first network, serves as the soft matrix while the first network acts as the tough filler. Compared to each single network, the DN hydrogels demonstrated a significant enhancement in both strength and toughness. For further different types of DN hydrogels, readers can refer to the work of Nakajima et al. [59], Zheng et al. [60], Muroi et al. [61], Bakarich et al. [62], and the reviews of Chen et al. [63] and Xu et al. [41].

In DN hydrogels, the rigid first network serves as sacrificial chains through the rupture of covalent bonds, which dissipates a substantial amount of mechanical energy, while the soft secondary network maintains material configuration during large strain deformation [2, 64]. This bond fracture process can lead to a complete change in the virgin elastic property, as demonstrated by the large hysteresis observed in the first loading–unloading cycle [16, 17, 27]. When the reloading strain exceeds the maximum strain previously applied, the stress returns to the same path as the virgin material. This phenomenon is similar to that observed in filled rubbers and also named as the Mullins effect. As shown in Fig. 1c, the hysteresis loop of DN hydrogels in the second and subsequent reloading–unloading cycles remains almost undetectable if deformation does not exceed the maximum value of the first cycle. This is different from filled rubbers and NC hydrogels, in which the reloading stress is often higher than the previous unloading stress, indicating certain recovery from the damage, as shown in Fig. 1a, b. This recovery effect can be more pronounced when the sample is relaxed for a few hours [25, 65]. However, stress softening in DN hydrogels cannot be recovered after the sample is left for even one week [16]. The irrecoverable nature of the Mullins effect in DN hydrogels demonstrates the discrepancy of the underlying damage mechanism compared to filled rubbers. As mentioned above, most studies have attributed the stress softening in DN hydrogels mainly to permanent rupture of the chemical covalent bonds in the rigid network [2, 64], while in filled rubbers, this issue is much more complex and can be explained by various mechanisms [25].

The irreversible rupture of covalent bonds results in limited repeatability of DN hydrogels after initial loading. To address this issue, some groups used non-covalent interactions to replace the covalent cross-link in the rigid network or in both networks, leading to the so-called hybrid or physical DN hydrogels [41, 66]. These non-covalent interactions, including metal ion coordination [67, 68], hydrogen bonding [69, 70], hydrophobic interactions [71, 72], and host–guest interactions [73, 74], act as reversible sacrificial bonds to dissipate energy while imparting self-recovery properties. Although each single non-covalent interaction is weaker than the covalent bond, a large amount of them provide toughness comparable to classical covalent DN hydrogels. However, their mechanical strength is often lower [41].

Typically, the stress response in filled rubbers also depends on the strain rate, with residual strains able to recover after hours of relaxation [15, 19, 34]. In contrast, the stress response of the covalent DN hydrogels shows negligible dependence on the deformation rate, and the residual strain is fairly small [16, 17, 75]. Hence, unloading and reloading are often treated as elastic processes without appreciable viscoelastic effects, as the deformation does not exceed its historical maximum value. However, the hybrid or physical DN gels with self-recovery effect can exhibit larger residual strains and viscoelastic effects [75, 76].

2.4 Multi-network Elastomers

Inspired by DN hydrogels, Ducrot et al. [9] successfully synthesized multi-network (MN) elastomers with high strength and toughness, exceeding traditional particle reinforcing schemes. They synthesized the DN elastomer by swelling the cross-linked single network with a second monomer and a small amount of cross-linker, followed by polymerization. This procedure can be repeated to synthesize triple-network and quadruple-network elastomers. As discussed in Ducrot et al. [9], the toughing mechanism is similar to that of DN hydrogels: the fracture of covalent bonds in the first network controls stress levels (and hence strength), whereas the second and third networks prevent the material from cracking. Consequently, MN elastomers also exhibit a hysteresis loop under cycling loading, as shown in Fig. 1d. However, Ducrot et al. [9] observed that the hysteresis loop of their DN elastomer is much smaller than that of the triple-network elastomer or quadruple-network elastomer. It is shown that the prestretching ratio of the filler network is the most determining parameter for mechanical performance. By replacing the C–C bond in the first network with BADOBA (a chemoluminescent bond that can emit light when it breaks), they found that the bond breaking only occurs above a certain level of the applied load in the loading procedure of the first cycle. Moreover, the damage zone is much larger in the MN elastomer than that in the single-network elastomer, which is consistent with the results in Chen et al. [77]. Compared to filled rubbers, the damage of MN elastomers shows nearly constant yield stress, which is related to the bond break stress of the first network.

3 Recent Developments toward Modeling the Mullins Effect in Soft Composites

As discussed above, fabricating tough soft composites (filled rubbers, NC gels, DN gels, MN elastomers) mostly relies on creating dissipative mechanisms to release the stored elastic energy of the polymer network, which can be achieved through either particle reinforcing schemes or designing multiple networks. Along with the enhanced strength and toughness, such dissipative mechanisms also result in substantial softening upon the first loading. Hence, the constitutive modeling of these soft composites is often beyond the scope of classical hyperelasticity theory, which was well established for rubber-like materials in the literature [78,79,80]. In this section, we present an overview of recent developments in quantitatively describing the Mullins effect in these soft composites. As already discussed in the review of Diani et al. [25], various constitutive models for the Mullins effect have been reported. These sophisticated models were mostly developed for simulating filled rubbers. Some of them, in a practical sense, may also be applied to simulate the stress softening behavior of emerging soft composite regardless of the underlying mechanisms. On the other hand, more specialized models concerned with tough gels and multi-network polymers have also been developed in recent years. In the current study, we will pay particular attention to these specialized studies. We shall also cover several of the newest developments for filled rubbers reported after Diani et al. [25].

The adequate description of the Mullins effect depends on the proper understanding of damage mechanisms and a reliable polymer network model, which usually appears as a hyperelastic model. Recent developments were mainly achieved based on new concepts regarding these two aspects. In the following, we summarize recently developed damage models for soft composites. Rather than chronicling them, we will first discuss the damage mechanism and then come to the models based on new network representations. For simplicity and readability, some original formulas and symbols are rewritten within the unified definitions of the current work. In most cases, constructing damage models for soft materials is to define a specified function of the free energy (also termed as strain energy) density. Stress can be worked out through standard procedure in continuum mechanics. We only discuss the construction of free energy functions here for simplicity. The detailed derivation of related stress can be found in textbooks such as Holzapfel [81] and Gurtin [82].

It should be mentioned that the list of models discussed below is surely not exhaustive. There are other excellent damage models that may not be included due to our focused idiosyncrasy and negligence. Readers can also refer to Lei et al. [83] for a review of recently developed hydrogel network models and Xiang et al. [26] for a review of constitutive models for both elastic and damage behavior of soft materials.

3.1 Models with New Damage Mechanism

Understanding and modeling the damage mechanism have long been pursued for soft composites. For NC gels, the reinforcing scheme is quite similar to filled rubbers. They both incorporate fillers that attach polymer chains and exhibit stress softening, residual strain, and recovery effect. Hence, it might be expected that they should share some similar origins of the Mullins effect. In fact, the Mullins effect in filled rubbers has already been explained by various physical mechanisms in the literature, such as bond rupture [33, 84], molecular slip [85], network rearrangement [24], etc. However, as discussed in Diani et al. [25], none of these explanations has achieved general agreement. Compared to filled rubbers and NC hydrogels, the stress softening in DN gels and MN elastomers seems to receive a more general consensus in underlying mechanisms. It is found that the damage observed in the first loading–unloading cycle is mostly attributed to bond scission in the filler network [9, 64]. In the following, we will give a brief introduction to recently proposed models concerning these distinct damage mechanisms in soft composites.

3.1.1 Model of Wang et al.

Wang et al. [86] developed a constitutive model that accounts for the pullout effect in the matrix network of DN hydrogels. In their model, the first network (i.e., the filler network) with shorter polymer chains (usually a polyelectrolyte) is assumed to be highly stretched even in the free swelling state. Hence, it fails much earlier than the second network with longer chains. Due to the heterogeneity of the polymer network, the partial damage splits the first network into subnetworks. Subnetworks of the divided first network, along with the compliant second network, form a new composite structure with a compliant matrix and stiff inclusions. The damage process is sketched in Fig. 2.

Fig. 2

Schematic diagram of the partial damage process in a DN gel: the first network fractures and forms stiff subnetworks, from which the second network is pulled out. This mechanism was adopted in Wang et al. [86]

The damage model is then developed based on the classical Gent model [87]:

$$ W\left( {J_{1} ,\eta ,J_{{\text{m}}} } \right) = - \frac{\mu }{2}\eta J_{{\text{m}}} \ln \left( {1 - \frac{{J_{1} }}{{J_{{\text{m}}} }}} \right) $$

(1)

where $\mu$ is the initial modulus of the virgin material, and $J_{1} = \lambda_{1}^{2} + \lambda_{2}^{2} + \lambda_{3}^{2} - 3$ with $\lambda_{i}$ being the principal stretch ratio. The variables $\eta$ and $J_{{\text{m}}}$ are two internal variables that characterize the damage process, of which $\eta$ measures the modulus reduction due to the split of the first network, and $J_{{\text{m}}}$ characterizes the stretching limit of the damaged phase.

The variable $J_{{\text{m}}}$ evolves due to the pullout effect of the second network. Wang et al. [86] assumed that $\eta$ depends on the historical maximum value of the principal stretch, namely $\eta = \eta \left( {\lambda_{\max } } \right)$, while $J_{{\text{m}}}$ depends on the historical maximum of the invariant $J_{{\text{m}}} = J_{{\text{m}}} \left( {J_{1\max } } \right)$. The specified formulations are respectively chosen as a log-normal distribution and a linear equation as

$$ \begin{gathered} \eta = 1 - \frac{{\eta_{0} }}{2}\left[ {{\text{erf}} \left( {\frac{1}{\sqrt 2 m}\ln \frac{{\lambda_{\max } - 1}}{{\lambda_{0} - 1}}} \right) + 1} \right] \hfill \\ J_{{\text{m}}} = J_{0} + \alpha J_{1\max } \hfill \\ \end{gathered} $$

(2)

where $\lambda_{0}$ is the most probably breaking stretch of polymer chains, $m$ is a positive parameter related to the width of the probability distribution, $\eta_{0}$ is a dimensionless normalization parameter, $J_{0}$ is the stretching limit of the virgin material, and $\alpha$ is a dimensionless material parameter. In the above, ${\text{erf}} (x) = \frac{2}{\sqrt \pi }\int_{0}^{x} {\text{e}^{{ - y^{2} }} } {\text{d}{y}}$ is the error function.

This model can capture the Mullins effect and the stable necking phenomenon of a DN gel. As discussed in Nakajima et al. [88], it was found that the two networks of DN hydrogels are partially covalently bonded to each other rather than forming separate interpenetrated networks. Shams et al. [89] further reported that DN hydrogels with high mechanical properties are only obtained when the chains of the second network are grafted on the first network. These developments demonstrate stronger interactions between the two networks than the interpenetrated models. However, very few existing models deal with these interaction effects. We note that another damage model in the work of Zhao [27] also accounts for the interplays between two networks based on a newly developed elastic interpenetrating network model, which will be detailed in the next subsection.

3.1.2 Model of Tang et al.

Tang et al. [90] developed a pseudo-elastic model for NC hydrogels. They assumed that during unloading, polymer chains prefer to contract to their initial state, while clay sheets cannot recover to their random configuration. Thus, the chains will probably slip on the surface of clay sheets, leading to desorption and a dissipation of elastic energy during unloading. This plausible mechanism stands up for the debatable assumption of the pseudo-elasticity theory that damage only evolves in the loading stage. As depicted in Fig. 3, the clay sheets are randomly distributed in the reference state (a). When the sample is loaded in state (b), polymer chains stretch, and in the meantime, clay sheets and chains orient along the tensile direction. In the unloaded state (c), polymer chains tend to recover their initial state, while clay sheets cannot recover to their random configuration. Thus, polymer chains slip on the surface of clay sheets, leading to desorption (the red chains), which causes the failure of physical bonds between clay sheets and polymer chains and dissipates mechanical energy. In the relaxed state (d), clay sheets maintain their oriented configuration, and some of the desorbed chains re-attach to the surface of the clay sheets (the green chains). When the sample is reloaded for the second cycle in state (e), the load-bearing of re-adsorbed chains results in a recovery effect between unloading and reloading stresses. As the reloading exceeds the previous maximum stretch, relaxed polymer chains become tight, and the stress follows the primary loading curve again. In the re-unloaded state (f), more polymer chains slip on the surface of clay sheets, inducing a larger hysteresis than in the first cycle.

Fig. 3

A schematic diagram of the damage mechanism of NC hydrogels

The reduction in the free energy is represented by multiplying a softening factor $\eta$ to the Ogden model:

$$ \tilde{W}\left( {\lambda_{1} ,\lambda_{2} ,\lambda_{3} ,\eta } \right) = \eta W\left( {\lambda_{1} ,\lambda_{2} ,\lambda_{3} } \right) = \eta \cdot \sum_{n} {\frac{{\mu_{n} }}{{\alpha_{n} }}} \left( {\lambda_{1}^{{\alpha_{n} }} + \lambda_{2}^{{\alpha_{n} }} + \lambda_{3}^{{\alpha_{n} }} - 3} \right) $$

(3)

Here, $W\left( {\lambda_{1} ,\lambda_{2} ,\lambda_{3} } \right) = \sum\limits_{n} {\frac{{\mu_{n} }}{{\alpha_{n} }}} \left( {\lambda_{1}^{{\alpha_{n} }} + \lambda_{2}^{{\alpha_{n} }} + \lambda_{3}^{{\alpha_{n} }} - 3} \right)$ is the undamaged Ogden model. The softening function is further specified as

$$ \eta = \frac{1}{2}\left[ {1 + {\text{erf}} \left( {\frac{1}{\sqrt 2 m}\ln \left( {\frac{{W_{0} }}{{W\left( {\lambda_{\max } } \right) - W(\lambda )}}} \right)} \right)} \right] $$

(4)

with two parameters $W_{0}$ and $m$, respectively, related to the most probable bond energy density and the distribution width. Here, $\lambda_{\max }$ is the historical maximum uniaxial stretch.

The model is demonstrated to properly describe the experimental hysteresis upon different deformations of NC hydrogels in the uniaxial case. Unlike the model in Ogden et al. [22], the free energy density in Tang et al. [90] is not modified by adding an additional damage function to form a pseudo-elastic energy, thus bypassing zero dissipation in the classical pseudo-elastic theory [91]. A similar approach was also adopted in Qi et al. [92] to describe the fracture toughness and hysteresis of soft materials.

3.1.3 Model of Lavoie et al.

Lavoie et al. [93] developed a damage model for soft elastomers based on the concept of progressive damage, in which polymer chains are assumed to have varying lengths and damage is induced progressively from shorter to longer chains. In the model, the total strain energy density is given by the sum of contributions from chains of different lengths:

$$ W = \sum_{j = 1}^{M} {P_{j} } w_{j} - p(J - 1) $$

(5)

where $P_{j}$ is the number fraction of chains with a length of $N_{j}$ (here the chain length refers to the number of chain segments), $w_{j}$ is the undamaged strain energy of chains with a length of $N_{i}$, and $M$ is the maximum possible chain length. Moreover, $J$ is the determinant of the deformation gradient, and $p$ is a Lagrange multiplier to enforce incompressibility.

The undamaged strain energy is given by the Langevin single-chain model and the eight-chain network representation:

$$ w_{j} = k_{B} TN_{j} \left( {\frac{\lambda }{{\sqrt {N_{j} } }}\beta_{j} + \ln \frac{{\beta_{j} }}{{\sinh \beta_{j} }}} \right),\quad \lambda = \sqrt {\frac{{I_{1} }}{3}} $$

(6)

where $k_{B}$ is the Boltzmann constant, $T$ is the absolute temperature, $\beta_{j} = {\mathcal{L}}^{ - 1} \left( {\lambda /\sqrt N_{j} } \right)$, ${\mathcal{L}}^{ - 1} ()$ is the inverse Langevin function, and $I_{1}$ is the first invariant of Cauchy–Green tensor. It should be noted that Eq. (6) is rewritten from the original formula by replacing $\frac{r}{L}$ in Lavoie et al. [93] with the chain stretch $\lambda$ to ensure a consistent formulation throughout the paper.

The theory of mechanochemistry [94] is further employed to calculate the density of surviving polymer chains due to scission:

$$ - \dot{P}_{j} = \frac{{N_{j} P_{j} }}{\tau }{\text{e}}^{{{L\beta_{j}{{\text{max}}}} /b}} $$

(7)

where $\tau$ and $L$ are fixed parameters, respectively, representing the relaxation time and activation length for bond dissociation, $\beta_{j}^{{{\text{max}}}}$ is the maximum $\beta_{j}$ during deformation, and $b$ is the Kuhn length. The above rate equation was solved with an initial condition of Maxwell–Boltzmann distribution. Lavoie et al. [93] demonstrated the successful application of their model in simulating the Mullins effect for the GR-S tread vulcanizate and triple-network elastomers.

3.1.4 Model of Bacca et al.

Bacca et al. [95] also developed a damage model for MN elastomers based on a similar progressive damage concept. The model assumes that the filler network has a non-uniform distribution for chain extensibility, and the least extensible chains (i.e., the shortest) break first. Different from the model of Lavoie et al. [93], the ruptured shorter chains do not simply disappear. Instead, they transfer to more extensible ones to bear the applied stress. As the total number of chain monomers is constant, such progressive damage leads to longer average extension limits, analogous to the network alternation model proposed by Marckmann et al. [24]. The undamaged network is treated as a hybrid form combining the neo-Hookean model and the Gent model as:

$$ W_{0} = \frac{{\mu_{0} }}{6}\left[ {I_{1} - 2\ln J - 3 - 2I\ln \left( {\frac{{I - I_{1} + 2\ln J}}{I - 3}} \right)} \right] - p\left( {J - 1} \right) $$

(8)

In the above, $\mu_{0}$ is the shear modulus, $I_{1}$ is the first invariant of the Cauchy–Green tensor, $J = {\text{det}}{\varvec{F}}$ is the swelling ratio of the MN elastomer, $p$ is the hydrostatic pressure, and $I$ is the chain extension limit.

To model the damage behavior of the MN elastomers, the filler networks are assumed to possess a range of chain extensibility constraints from $I_{s}$ (shortest chains) to $I_{l}$ (longest chains). Then, the strain energy density is rewritten as:

$$ W_{f} = \frac{{\mu_{0} }}{6}\left[ {I_{1} - 2\ln J - 3 - 2\int_{{I_{s} }}^{{I_{l} }} I \ln \left( {\frac{{I - I_{1} + 2\ln J}}{I - 3}} \right)f(I){\text{d}{I}}} \right] - p\left( {J - 1} \right) $$

(9)

where $f(I)$ is the weight of chains with extensibility of $I$, and $\int_{{I_{s} }}^{{I_{l} }} f (I){\text{d}}I = {1}$.

In traditional progressive damage models, a distribution of chain extensibility is often assumed. When deformation is applied, a progressive rupture of chains is observed from short chains to long chains, which describes damage evolution. The process terminates when the longest chains are broken, and total rupture of the network is obtained. In the model of Bacca et al. [95], however, an alternative idea is adopted to simulate this progressive damage. It is assumed that in each loading step, chains with extensibility shorter than $I_{c}$ are ruptured, while chains with longer extensibility are loose without bearing stress. Here $I_{c}$ is associated with the least extensible chains yet available in the microstructure. Therefore, during the damage process, the weight factor can be approximated by the Dirac delta function: $f(I) = \delta \left( {I - I_{c} } \right)$. Then, the strain energy density in Eq. (9) can be recast as

$$ W_{f} = \frac{{\mu_{0} }}{6}\left[ {I_{1} - 2\ln J - 3 - 2I_{c} \ln \left( {\frac{{I_{c} - I_{1} + 2\ln J}}{{I_{c} - 3}}} \right)} \right] - p\left( {J - 1} \right).\quad I_{s} \le I \le I_{l} $$

(10)

Consequently, the undamaged polymer has $I_{c} = I_{s}$. As damage initiates and progresses, $I_{c}$ will continuously increase up to the value of $I_{l}$ when the total rupture is obtained. A damage parameter $d$ is adopted to describe the evolving $I_{c}$:

$$ I_{c} = I_{s} (1 - d) + I_{l} d,\quad 0 \le d \le 1 $$

(11)

The damage law is specified as

$$ d = \left\{ {\begin{array}{*{20}l} {0,} \hfill & {I_{1} \le I_{1}^{{\text{i}}} } \hfill \\ {\left( {\frac{{I_{1} - I_{1}^{{\text{i}}} }}{{I_{1}^{{\text{r}}} - I_{1}^{{\text{i}}} }}} \right)^{q} ,} \hfill & {I_{1}^{{\text{i}}} \le I_{1} \le I_{1}^{{\text{r}}} } \hfill \\ \end{array} } \right. $$

(12)

where $I_{1}^{{\text{i}}}$ is the invariant when damage is initiated, $I_{1}^{{\text{r}}}$ is the invariant of total rupture, and $q$ is a dimensionless parameter. The above damage model was further applied to simulate MN elastomers by incorporating the swelling and mixing effect. Since we are mainly concerned with the damage modeling, these details are not introduced here. Bacca et al. [95] demonstrated that their model can successfully describe the stress-stretch curves for a triple-network elastomer subject to cyclic loading.

3.1.5 Model of Wang et al.

Wang et al. [28] developed a constitutive model for NC hydrogels with both damage and recovery effects. In their model, the polymer network of NC hydrogel is described by chains of inhomogeneous length between two nanoparticles. These chains can detach from the particles, becoming either inactive or longer active chains in the loading process. During unloading, long chains can be adsorbed by particles, and previously detached short chains can re-attach to the particles. This network evolution process is depicted in Fig. 4.

Fig. 4

Evolution of the polymer network of the nanocomposite hydrogel: a the polymer network in the undeformed state; b under loading, polymer chains can detach from the particles, becoming either inactive (blue and black chains) or longer active chains (red chains); c under unloading, long chains can graft on the middle particles (red and green chains), and previously detached short chains can re-attach to the particles (blue chain). (Color figure online)

Wang et al. [28] used the classical eight-chain model to describe the network structure of the NC hydrogel: in a unit cube, eight nanoparticles are located in the corners and one nanoparticle at the center. Each corner particle and the center particle are connected by a set of polymer chains with non-uniform chain length distribution. The free energy density of chains with a length of $N_{i}$ is given by:

$$\begin{aligned} w_{i} &= k_{B} TN_{i} \left( {\frac{{\beta_{i} }}{{\tanh \beta_{i} }} + \ln \frac{{\beta_{i} }}{{\sinh \beta_{i} }}} \right),\quad\\ \beta_{i} &= {\mathcal{L}}^{ - 1} \left( {\lambda /\sqrt N_{i} } \right),\quad \lambda = \sqrt {I_{1} /3}\end{aligned} $$

(13)

Noting that the Langevin function is defined by ${\mathcal{L}}(x) = \coth x - \frac{1}{x}$, the above formula is essentially identical to Eq. (6) used in Lavoie et al. [93].

During loading, the short chains detach from the particles, leaving a renewed network with only long chains between two particles. Thus, the minimum chain length can be calculated by equating the chain force to the critical bonding strength, leading to

$$ N_{c} = \frac{1}{{{\mathcal{L}}\left( {\frac{{bf_{{\text{c}}} }}{{k_{B} T}}} \right)}}\frac{D}{b}\sqrt {\frac{{I_{1} }}{3}} $$

(14)

where $f_{c}$ is the bonding strength between the chains and the nanoparticles, $b$ is the Kuhn length, and $D$ is the distance between the corner particle and center particle.

After detachment, chains shorter than $N_{c}$ transform into chains longer than $N_{c}$. Therefore, the number of active long chains increases. An exponential function is used to model this effect:

$$ \frac{{n_{ia} }}{{n_{i} }} = \exp \left[ {\alpha_{1} \left( {\sqrt {\frac{{I_{1}^{\max } }}{3}} - 1} \right)} \right] $$

(15)

where $n_{i}$ and $n_{ia}$ are the numbers of chains with a length of $N_{i}$ before and after network evolution, and $I_{1}^{{{\text{max}}}}$ is the historical maximum value of $I_{1}$. Here, the chain length is restricted to an interval of $N_{c} \le N_{i} \le N_{m}$ with $N_{m}$ being the maximum chain length.

During unloading, longer chains may re-attach to the particles, which is simulated by another exponential function:

$$ \frac{{n_{ia} }}{{n_{i} }} = \exp \left[ {\alpha_{1} \left( {\sqrt {\frac{{I_{1}^{\max } }}{3}} - 1} \right)} \right]\exp \left[ { - \alpha_{2} \left( {\sqrt {\frac{{I_{1}^{\max } }}{3}} - \sqrt {\frac{{I_{1} }}{3}} } \right)} \right] $$

(16)

The above equation is restricted to $N_{c} \le N_{i} \le N_{{{\text{mc}}}}$ with $N_{{{\text{mc}}}}$ being the critical length corresponding to the maximum $I_{1}^{{{\text{max}}}}$. On the other hand, the attachment of longer chains contributes to the increase in short chains

$$ \frac{{n_{ia} }}{{n_{i} }} = \exp \left[ { - \alpha_{3} \left( {\sqrt {\frac{{I_{1}^{\max } }}{3}} - \sqrt {\frac{{I_{1} }}{3}} } \right)} \right] $$

(17)

with $N_{{{\text{mc}}}} \le N_{i} \le N_{{\text{c}}}$. In the above, $\alpha_{1}$ and $\alpha_{2}$ are two positive parameters, while $\alpha_{3}$ can be positive for $n_{ia} < n_{i}$ or negative for $n_{ia} > n_{i}$.

In the initial state, the weight of chains with a length of $N_{i}$ is defined as $P_{i} = n_{i} /{\sum\limits_{i = 1}^{m}} {n_{i} }$. Following Govindjee et al. [96] and Dargazany et al. [18], the weight distribution is given by

$$\begin{aligned} P_{i} &= \overline{P}\sqrt {\frac{{3\kappa^{2} }}{{2\pi N_{i} }}} \exp \bigg\{ - G - \frac{\kappa }{\sqrt \pi }\bigg[\sqrt {6N_{i} } {\text{e}}^{ - G} + 3\frac{D}{b}\sqrt {\pi N_{i} }\\ &\quad {\text{erf}} (\sqrt G ) - \sqrt 6 {\text{e}}^{{ - GN_{i} }} - 3\frac{D}{b}\sqrt \pi {\text{erf}} \bigg( {\sqrt {GN_{i} } } \bigg) \bigg] \bigg\} \end{aligned}$$

(18)

with $G = 3D^{2} /\left( {2N_{i} b^{2} } \right)$ and $\kappa$ being a parameter related to the particle surface activity. $\overline{P}$ is a factor that ensures the total density equals 1.

Finally, the total free energy density is given by

$$ W = 4\rho n\sum_{i = c}^{m} {\left( {w_{i} P_{i} } \right)} + W_{m} $$

(19)

where $\rho$ is the particle concentration, $n = \sum_{i = 1}^{m} {n_{i} }$ is the average chain number between two particles, and $W_{m}$ is a constant corresponding to the mixing free energy of polymer chains and water due to solvent interactions. Later, the same group [97] developed another rate-dependent theory for the self-recovery effect in NC gels, in which the free polymer chains diffuse across the interface and reform cross-links to bridge the interface. It is revealed that the healed strength of nanocomposite hydrogels increases with the healing time in an error-function-like form. The equilibrium self-healing time of the full-strength recovery decreases with the temperature and increases with the nanoparticle concentration. Moreover, the healed interfacial strength decreases with increasing delaying time before the healing process.

It is still unclear whether the eight-chain network representation used by Wang et al. [28, 97], as well as other classical network models (e.g., the three-chain model, four-chain model, and full network model) can properly describe the network structure of NC gels, given the large number of polymer chains connected to the nanofillers, resulting in more complex network structures than traditional filled rubbers. Nevertheless, the theoretical results in Wang et al. [28, 97] align quantitatively with several experiments on NC hydrogels. It is noted that Iyer et al. [98, 99] have employed molecular dynamic simulations to study the mechanical deformation and strain recovery of nanoparticle cross-linked polymer networks. However, due to the limited computational size scales, it remains unclear how well the computational results match the macroscopic stress–strain behaviors of NC hydrogels.

3.1.6 Model of Vernerey et al.

Vernerey et al. [100] formulated a statistical model in which the damage behavior of MN polymers is captured by tracking the evolution of chain configurations. In the model, damage manifested as softening and hysteresis under cyclic loading is explained by the scission of chains in the rigid network, which is further described by an accumulative probability function that characterizes the likelihood of failure of a chain with a specified end-to-end distance (Fig. 5). For a set of chains with the end-to-end vector ${\varvec{R}}$ in the initial configuration, the chain density $\Phi_{0}$ is assumed to follow the Maxwell distribution:

$$ \Phi_{0} ({\varvec{R}}) = c_{0} \left( {\frac{3}{{2\pi Nb^{2} }}} \right)^{3/2} \exp \left( { - \frac{{3{\varvec{R}} \cdot {\varvec{R}}}}{{2Nb^{2} }}} \right) $$

(20)

with $c_{0}$ being the total number of chains. Here, the chain length $N$ is a constant. When assuming a non-uniform distribution $\Phi (N)$ for chain length $N$, the initial end-to-end distance will also follow a distribution of $\Phi \left( {\frac{{R^{2} }}{{b^{2} }}} \right)$ if the chain stretch is defined as $\lambda = \frac{R}{\sqrt N b} = 1$ in the standard entropic elasticity theory. In this sense, the model of Vernerey et al. [100] shares a similarity with the progressive damage theory.

Fig. 5

Illustration of the evolution for the chain distribution. Deformation can stretch the distribution without changing its overall area. Rupture of chains tends to decrease the number of chains with larger end-to-end distances

Distinct from classical progressive damage models, Vernerey et al. [100] used a function to describe the damage probability rather than adopting specific damage criteria. The chains with the current end-to-end distance $r$ are assumed to have a probability of failure by $P(r)$, given as

$$ P(r) = \frac{1}{{\sqrt {2\pi } \sigma_{{\text{d}}} }}\int_{0}^{r} {\exp } \left[ { - \frac{1}{{2\sigma_{{\text{d}}}^{2} }}\left( {\xi - r_{{\text{c}}} } \right)^{2} } \right]{\text{d}}\xi $$

(21)

with ${\varvec{r}}$ being specified by the affine assumption ${\varvec{r}} = {\varvec{FR}}$. In the above, $\sigma_{d}$ is the distribution deviation, and $r_{c}$ is the critical end-to-end distance. The function $P(r)$ remains nearly zero when the end-to-end distance is small and increases rapidly to 1 as the end-to-end distance approaches $r_{c}$.

The fraction of ruptured chains with the initial end-to-end vector $\varvec{R}$ at time $t$ is defined as $\Delta ({\varvec{R}},t)$. Then, the effective chain density $\Phi ({\varvec{R}},t)$ is given by eliminating the ruptured chains:

$$ \Phi ({\varvec{R}},t) = \Phi_{0} ({\varvec{R}})\left( {1 - \Delta ({\varvec{R}},t)} \right) $$

(22)

or alternatively in the current configuration:

$$ \phi ({\varvec{r}},t) = \phi_{0} ({\varvec{r}},t)\left( {1 - \delta ({\mathbf{r}},t)} \right) $$

(23)

with $\delta ({\varvec{r}},t) = \Delta \left( {{\varvec{F}}^{ - 1} {\varvec{r}},t} \right)$ and $\phi_{0} ({\varvec{r}},t) = \Phi_{0} \left( {{\varvec{F}}^{ - 1} {\varvec{r}},t} \right)$.

The damage evolution rule is specified by

$$ \begin{gathered} \dot{\Delta }({\varvec{R}},t) = \nabla P({\varvec{r}}) \cdot {{\varvec {r}}}, \, \nabla P({\varvec{r}}) \cdot {{\varvec{r}}} \ge 0 \hfill \\ \dot{\Delta }({\varvec{R}},t) = 0, \, \nabla P({r}) \cdot {{\varvec{r}}} < 0 \hfill \\ \end{gathered} $$

(24)

The free energy is finally given by the increased part between the initial and deformed states:

$$ W = \int {\left[ {\phi - \Phi } \right]} k_{B} TN\left( {\frac{\lambda }{\sqrt N }\beta + \frac{\beta }{\sinh \beta }} \right){\text{d}}\Omega_{r} $$

(25)

The above integral is calculated by the vector ${\mathbf{r}}$ passing through all configurations in the three-dimensional space, which can be done through numerical approaches, for example, Bazant et al. [101] and Miehe et al. [102].

3.1.7 Models with Recovery Effects

Some novel DN hydrogels with hybrid/physical bonds exhibit high strength/toughness, as well as rate-dependent viscoelasticity and recovery effect. These properties go beyond those of classical covalent DN hydrogels and are more similar to the Mullins effect in the filled rubbers and NC hydrogels. Yu et al. [103] developed a general theory that mathematically describes the damage and recovery behavior of dynamic physical bonds that can autonomously reform after fracture or dissociation. It is assumed that self-healing polymers comprise interpenetrating networks cross-linked by dynamic bonds. The network chains follow inhomogeneous chain length distributions, and the dynamic bonds obey force-dependent chemical kinetics. During the self-healing process, the polymer chains diffuse across the interface to reform the dynamic bonds, which is modeled by a diffusion–reaction theory. The theories can predict the stress-stretch behaviors of original and self-healed DPNs, as well as the healing strength as a function of healing time. Lin et al. [6] also developed a constitutive model to capture the complex mechanical behaviors of physical hydrogels. The theoretical framework incorporates two main physical mechanisms: the microscopic kinetic dissociation/association of cross-links and inelastic conformational changes of gel networks. The model predictions agree well with the experimental data of gel responses under different scenarios, which simultaneously capture the time-dependent stress response, relaxation, hysteresis, self-healing behavior, and ultimate failure. In the current study, we will not further discuss self-healing models with rate dependence. Readers can refer to Yu et al. [103] and Lin et al. [6] and the references therein.

3.2 Developments Based on New Network Models

Besides the damage mechanism, the modeling of the Mullins effect also strongly depends on the proper description of the polymer network. Usually, the network model is concerned with two basic issues: the single-chain model and the micro–macro-transition. The single-chain model specifies the free energy function with respect to the microscopic chain stretch, while the micro–macro-transition relates the microscopic chain stretch to the detectable macroscopic deformation. However, it has been found that classical single-chain models, together with the widely used micro–macro-transitions, do not adequately describe polymeric materials under complex multiaxial loading [80, 104]. Therefore, new network models have been proposed, and several related damage models have been developed based thereon.

3.2.1 Zhao’s Model

Zhao [27] developed a theory for the large deformation and damage behaviors of DN hydrogels based on an interpenetrating polymer network (IPN) model. In the model, the network interpenetrating is described by two effects. As one network is stretched, the polymer chains in the other network are also stretched isotropically and homogeneously without externally applied forces, as shown in Fig. 6. Meanwhile, the density of chains in the other network decreases. The main formula of the model is summarized as follows. First, the total free energy of an IPN with $m$ networks is given by

$$ W_{{{\text{total}}}} = \sum_{i = 1}^{m} {W_{i} } + W_{{\text{M}}} $$

(26)

where $W_{i}$ is the free energy due to the deformation of the $i$-th network, and $W_{{\text{M}}}$ is the free energy due to mixing polymers and solvents.

Fig. 6

a Schematics of the interpenetrating eight-chain network model; b schematics of the breaking of chains and cross-links in the shorter-chain network

In Zhao’s model [27], a constant $W_{{\text{M}}}$ is further assumed without involving concentration change during deformation. Then, the total free energy is the sum of each network. For each network, the free energy is specified by the Langevin single-chain model and the eight-chain network representation:

$$ {\begin{aligned} & W_{i} = C_{i} n_{i} k_{B} TN_{i} \left( {\frac{{\beta_{i} }}{{\tanh \beta_{i} }} + \ln \frac{\beta }{{\sinh \beta_{i} }}} \right)\\ & \beta_{i} = {\mathcal{L}}^{ - 1} \left( {\overline{\lambda }_{i} /\sqrt N_{i} } \right),\quad \overline{\lambda }_{i} = C_{i}^{ - 1/3} \sqrt {\frac{{I_{1} }}{3}} \end{aligned}}$$

(27)

In the above, $C_{i}$ is the volume concentration of the $i$-th network in the IPN, and $n_{i}$ is the chain density. Here, the stretch $\overline{\lambda }_{i}$ consists of two parts: $C_{i}^{ - 1/3}$ due to the interpenetration of other networks and solvent, and $\sqrt {I_{1} /3}$ due to the deformation of IPN.

The above network model was adopted to simulate DN hydrogels. The damage in the rigid first network is described by the network alternation theory [24], with increasing chain length and decreasing chain density specified by exponential functions [105]:

$$ \begin{array}{*{20}l} {N_{A} = N_{A0} \exp \left[ { - p\left( {\overline{\lambda }_{A}^{\max } - C_{A}^{ - 1/3} } \right)} \right]} \hfill \\ {n_{A} = n_{A0} \exp \left[ {q\left( {\overline{\lambda }_{A}^{\max } - C_{A}^{ - 1/3} } \right)} \right]} \hfill \\ \end{array} $$

(28)

where $N_{A0}$ and $n_{A0}$ are the average chain length and chain density of the rigid network (denoted as network A) in the undeformed state, respectively; $p$ and $q$ are two material parameters with the relation of $p \ge q \ge 0$ due to thermodynamic restriction. This is different from the original work of Marckmann et al. [24], according to which the parameters $p$ and $q$ should be equal to enforce the condition that $nN = {\text{constant}}$. Through several examples, Zhao [27] demonstrated the capability of the model in quantitatively characterizing the Mullins effect and necking instability experimentally observed in DN hydrogels.

3.2.2 Model of Mao et al.

In an early publication of Lake and Thomas [106], it was asserted that the rupture of polymer chains was primarily due to internal energy rather than entropy. However, traditional single-chain models do not account for the extension of monomer units and, hence, internal energy. To address this issue, Mao et al. [107] proposed a model considering bond deformation based on the Langevin chain model, allowing for both entropic and energetic contributions. In this model, the total chain stretch is given by a multiplicative decomposition with the configuration stretch $\lambda_{{\text{c}}}$, and the monomer stretch $\lambda_{{\text{m}}}$:

$$ \lambda = \lambda_{{\text{c}}} \lambda_{{\text{m}}} $$

(29)

as depicted in Fig. 7.

Fig. 7

Single-chain model with bond stretch

The free energy of a single chain is then given as

$$ {\begin{aligned} & w = w\left( {\lambda_{{\text{m}}} ,\lambda_{{\text{c}}} } \right) = k_{B} TN\left[ {\frac{{\lambda_{c} \beta }}{\sqrt N } + {\text{ln}}\frac{\beta }{{{\text{sinh}}\beta }}} \right] + \varepsilon \left( {\lambda_{{\text{m}}} } \right)\\ & \beta = {\mathcal{L}}^{ - 1} \left( {\frac{{\lambda_{c} }}{\sqrt N }} \right) \end{aligned}}$$

(30)

The first term in the above function is the configuration entropic part given by the Langevin model. The term of $\varepsilon \left( {\lambda_{{\text{m}}} } \right)$ is the energetic part due to bond deformation.

To enforce a minimum total energy, the bond deformation is related to the total stretch through the stationary condition

$$ \frac{{\partial w(\lambda ,\lambda_{\text{m}} )}}{{\partial \lambda_{\text{m}} }} = 0 $$

(31)

The Hencky elastic model with a quadratic form of the logarithmic strain is adopted to model the energetic part:

$$ \varepsilon \left( {\lambda_{{\text{m}}} } \right) = \frac{1}{2}NE_{\text{m}} ({\text{ln}}\lambda_{\text{m}} )^{2} $$

(32)

where $E_{\text{m}}$ is the monomer stiffness.

The above single-chain model is further extended to a continuum model through the eight-chain network representation with

$$ W = nk_{B} TN\left[ {\frac{{\lambda_{c} \beta }}{\sqrt N } + {\text{ln}}\frac{\beta }{{{\text{sinh}}\beta }}} \right] + \varepsilon \left( {\lambda_{{\text{m}}} } \right),\quad \lambda = \sqrt {I_{1} /3} $$

(33)

with $n$ being the chain density.

Mao et al. [107] postulated that the scission of a chain occurs when the internal energy reaches a critical value. Specifically, the polymer chain will break as the bond is stretched to $\lambda_{\text{m}} \approx 1.4$ at room temperature, indicating a total chain stretch of $\lambda \approx 1.4\sqrt N$ that considerably exceeds the locking stretch of the classical model. The internal energy of bond stretching occupies 95% of the total free energy at the point of scission, which agrees with the assumption of Lake and Thomas [106]. In fact, the idea of stretchable bonds dates back to the early work of Smith et al. [108] on the elastic response of DNA molecules under overstretching. We note that other physically based formulations for the internal energy of the chains have also been developed [109, 110], which is not detailed here due to the lack of relevant studies concerning the Mullins effect.

3.2.3 Model of Lu et al.

In another relevant study, Lu et al. [111] introduced an amplification factor $D$ to account for the release of the folded domains in a single polymer chain, resulting in a similar expression of the chain free energy:

$$ {\begin{aligned} & w(\lambda ,D) = kT\left[ {N\left( {\frac{\lambda }{D\sqrt N }\beta + \ln \frac{\beta }{\sinh \beta }} \right) + \Phi (D)} \right],\\& \beta = {\mathcal{L}}^{ - 1} \left( {\frac{\lambda \beta }{{D\sqrt N }}} \right) \end{aligned}} $$

(34)

Comparing Eq. (34) with Eqs. (29)–(32), it is found that the internal variable $D$ plays a similar role to the monomer stretch $\lambda_{\text{m}}$. Lu et al. [111] then regarded the internal variable $D$ as an amplification factor to contain some folded domains, i.e., the “hidden length.” Prior to the stiffening of the polymer chain in the main network, the hidden length is released, leading to stress softening. Similar to Mao et al. [107], the stationary condition

$$ \frac{\partial w(\lambda ,D)}{{\partial D}} = 0 $$

(35)

and the following function

$$ \Phi (D) = \frac{1}{2}\kappa ({\text{ln}}D)^{2} $$

(36)

are adopted. The single-chain model is then mapped to a macroscopic one through the eight-chain network model. Lu et al. [111] showed the capability of the model in simulating the Mullins effect of DN hydrogels under cyclic loading conditions. This specific free energy density function can also be employed to model the mechanical behaviors of MN elastomers.

3.2.4 Model (a) of Xiao et al.

Xiao et al. [29] successfully applied the free energy of a single chain (Eq. (30)) to simulate the Mullins effect for MN elastomers. It is assumed that only the rigid filler network is damaged during deformation, while the matrix network exhibits a hyperelastic response without damage. For the latter, the classical hyperelastic eight-chain model is used. For the filler network, the single-chain model is given by Eqs. (29)–(32) according to Mao et al. [107]. The progressive damage approach is further adopted to formulate the damage behavior, with the free energy density written as:

$$ \begin{aligned} & W_{{\text{f}}} \left( {\lambda ,\lambda_{j}^{b} } \right) \\ & = n\sum_{j} {P_{j} } \left[ {k_{B} TN_{j} \left( {\frac{\lambda }{{\lambda_{j}^{b} \sqrt N_{j} }}\beta_{j} + \ln \frac{{\beta_{j} }}{{\sinh \beta_{j} }}} \right) + \frac{1}{2}N_{j} E_{b} \left( {\ln \lambda_{\text{m}}{j} } \right)^{2} } \right], \\ & \lambda = \sqrt {\frac{{I_{1} }}{3}} \end{aligned} $$

(37)

where $n$ is the chain density of the filler network, and $P_{j}$ is the weight fraction of chains with chain segment of $N_{j}$. The chain length distribution is assumed to be in the same form as Eq. (18).

The damage criterion is given by

$$ P_{j} = 0,{\text{ when }}\frac{{\lambda^{\max } }}{{\sqrt {N_{j} } }} \ge \lambda_{{{\text{cri}}}} $$

(38)

where $\lambda^{{{\text{max}}}} = \sqrt {I_{1}^{{{\text{max}}}} /3}$ is the maximum equivalent stretch in the synthesized and loading history, and $\lambda_{{\text{cri }}}$ is the critical stretch ratio.

The total free energy is then obtained through the sum of filler network energy specified in Eq. (37) and the matrix energy given by the eight-chain model. Xiao et al. [29] showed that the model can be successfully applied to describe the stress response of MN elastomers with different values of prestretching ratios of the filler network and the Mullins effect under cyclic loading conditions.

3.2.5 Model of Talamini et al.

In another development, Talamini et al. [112] introduced the damage effect in the energetic part based on a nonlocal damage approach. They first replaced the ε(λ_m) in the single-chain model of Eq. (30) with the following form:

$$ \varepsilon (d,\lambda_{\text{m}} ) = g(d)\frac{1}{2}NE_{\text{m}} \left( {\lambda_{\text{m}} - 1} \right)^{2} $$

(39)

Here, $d$ is a damage variable. In analogy to the phase field fracture theory (see, for example, Miehe et al. [113] and Wu et al. [114]), $g(d)$ is a monotonically decreasing degradation function with the property of $g(0) = 1,g(1) = 0,g^{\prime}(1) = 0$. The widely used quadratic form is adopted as

$$ g(d) = (1 - d)^{2} $$

(40)

The single-chain model is then extended to the macroscopic model by defining

$$ \lambda = \sqrt {\frac{{{\text{tr}}\overline{{\varvec{C}}} }}{3}} $$

(41)

where $\overline{{\varvec{C}}} = J^{ - 2/3} {\varvec{C}}$ is the distortional right Cauchy–Green tensor. The macroscopic damaged internal energy is further rewritten as:

$$ \overline{\varepsilon }(d,\lambda_{\text{m}} ) = (1 - d)^{2} \left[ {\frac{1}{2}E\left( {\lambda_{\text{m}} - 1} \right)^{2} + \frac{1}{2}K(J - 1)^{2} } \right] + \frac{1}{2}\varepsilon^{f} \ell^{2} |\nabla d|^{2} $$

(42)

where the second term in the square brackets is the volumetric contribution, the last term is a nonlocal part according to the phase field theory, $\varepsilon^{f}$ is a parameter related to the chain broken energy, and $\ell$ is an intrinsic length scale for the damage process.

The damage evolution rule is governed by the balance equation:

$$ \zeta \dot{d} = 2(1 - d){\mathcal{H}} - \varepsilon^{f} \left( {d - \ell^{2} \Delta d} \right) $$

(43)

with $\zeta$ being a kinetic modulus. In the above,${\mathcal{H}}(t)$ is a history field function defined by

$$ {\mathcal{H}}(t)\mathop = \limits^{{\text{ def }}} \max_{s \in [0,t]} \left\langle {\overline{\varepsilon }_{0} \left( {\lambda_{\text{m}} (s),J(s)} \right) - \frac{1}{2}\varepsilon^{f} } \right\rangle $$

(44)

with $\overline{\varepsilon }_{0} \left( {\lambda_{\text{m}} ,J} \right) = \frac{1}{2}E\left( {\lambda_{\text{m}} - 1} \right)^{2} + \frac{1}{2}K(J - 1)^{2}$.

The above boundary value problem can be solved through standard numerical procedures with the phase field method. It is noticed that Eq. (30) is a single-chain model at the microscopic level. At the continuum level, the fracture of materials is well explained by the Griffth theory that the reduction in elastic energy balances with rupture energy (the energy needed to create the crack surface). The Lake-Thomas model [106] explains how the entropic free energy function is used in fracture mechanics: the reduction in free energy, which is entropic for rubber, balances with the rupture of chains, which is energetic. In this sense, it seems that a nonlocal approach may be more suitable for the fracture behavior of materials: the chain is described by the model of Mao et al. in the crack region and by the entropic model in other regions. Hence, the model of Talamini et al. [112] seems to be a proper extension of Mao et al. [107] in the spirit of fracture mechanics. Most damage models for the Mullins effect were developed for the stable stage of materials before softening and necking. We note that Zhao’s model [27] and the model of Vernerey et al. [100] have shown the ability to describe the softening feature before necking. However, these local damage models may suffer from excessive mesh dependence due to the lack of an internal length scale in the numerical implementations, as discussed in Wu et al. [114]. Here, with the help of the nonlocal damage approach, the model of Talamini et al. [112] has the potential advantage for simulating the Mullins effect that accompanies the softening and necking process in the large strain region.

3.2.6 Model of Zhu et al.

Besides the Langevin single-chain model and its variants, Zhu et al. [115] developed a damage model for DN hydrogels based on another single-chain model, which is termed as the worm-like-chain (WLC) model:

$$ w = \frac{nkTNb}{{4A}}\left( {\frac{2\lambda }{N} + \frac{\sqrt N }{{\sqrt N - \sqrt \lambda }} - \sqrt {\frac{\lambda }{N}} } \right) $$

(45)

where $A$ is the persistence length of the chain.

Following the eight-chain network representation, the free energy density of the two networks are given as

$$ W_{j} = \frac{{C_{j} n_{j} kTN_{j} b}}{4A}\left( {\frac{{2\Lambda_{j} }}{{N_{j} }} + \frac{{\sqrt {N_{j} } }}{{\sqrt {N_{j} } - \sqrt {\Lambda_{j} } }} - \sqrt {\frac{{\Lambda_{j} }}{{N_{j} }}} } \right),\quad \Lambda_{j} = C_{j}^{ - 1/3} \sqrt {\frac{{I_{1} }}{3}} $$

(46)

where the quantities with the subscript $j = s$ or $j = l$ represent those of the short-chain or the long-chain network, and $C_{j}$ is the volume concentration of the $j$-th network.

The damage is assumed to be induced in the short-chain network with evolving chain density $n_{s}$ and chain length $N_{s}$. The mass conservation in both networks yields:

$$ n_{s} N_{s} kTb/(4A) = n_{s0} N_{s0} kTb/(4A) = \Phi_{s} $$

and representing as

$$ n_{l} N_{l} kTb/(4A) = \Phi_{l} $$

$$ n_{s} = \frac{{\Phi_{s} }}{{N_{s} }},\quad n_{l} = \frac{{\Phi_{l} }}{{N_{l} }},\quad \dot{n}_{s} = - \frac{{\Phi_{s} }}{{N_{s}^{2} }}\dot{N}_{s} $$

(47)

For uniaxial deformation, the energy dissipation function $\xi (\epsilon)$ is derived as

$$ \xi^{\prime } (\epsilon) = \left\{ {\begin{array}{*{20}l} {\left[ { - \frac{\partial W}{{\partial N_{s} }} + \frac{{\Phi_{s} }}{{N_{s}^{2} }}\frac{\partial W}{{\partial n_{s} }}} \right]N_{s}^{\prime } \left( {\epsilon_{\max } } \right)} \hfill & {{\text{ if }}\epsilon = \epsilon_{\max } \ge \epsilon_{{{\text{cri}}}} } \hfill \\ 0 \hfill & {\text{ otherwise }} \hfill \\ \end{array} } \right. $$

(48)

with $\epsilon = \lambda - 1$ being the uniaxial strain. Here, $\epsilon_{\max }$ is the historical maximum strain, and $\epsilon_{{{\text{cri}}}}$ is the critical strain for damage occurrence.

Given the specific form of energy dissipation $\xi (\epsilon)$, Eqs. (47)–(48) can give the evolution laws of two internal variables $n_{s} = n_{s} \left( {\epsilon_{\max } } \right)$ and $N_{s} = N_{s} \left( {\epsilon_{\max } } \right)$. A three-order polynomial function is used to fit the accumulated dissipation energy:

$$ \xi \left( {\epsilon_{\max } } \right) = a_{0} + a_{1} \epsilon_{\max } + a_{2} \epsilon_{\max }^{2} + a_{3} \epsilon_{\max }^{3} $$

(49)

with $a_{0}$, $a_{1}$, $a_{2}$, and $a_{3}$ being four coefficients to be determined through specific experimental data.

The capability and efficiency of the model has been demonstrated in simulating the Mullins effect of DN gels. In a more recent study, Zhu et al. [116] further proposed another form for the dissipation energy for the multiaxial deformation cases:

$$ \xi (J) = \frac{D}{2}\left[ {{\text{erf}} \left( {\frac{1}{\sqrt 2 m}\ln \left( {\frac{J}{{J_{0} }}} \right)} \right) + 1} \right] $$

(50)

with $J = I_{1} - 3$ being the multiaxial equivalent strain. Here, $D$, $m$, and $J_{0}$ are three material parameters. The multiaxial model has demonstrated good predictive capability for both uniaxial tension and pure shear.

3.2.7 Model of Lavoie et al.

Lavoie et al. [117] developed another progressive damage model for MN elastomers using the non-predetermined chain energy, which is alternatively arrived through integrating uniaxial force-extension relations. The total strain energy density of the model is represented by the sum of contributions from each network weighted by its volume fraction, which is written as

$$ W_{M} ({\varvec{B}}) = \sum_{i = 1}^{M} {\phi_{M}^{i} } W_{M}^{i} \left( {{\varvec{B}}_{M}^{i} } \right) $$

(51)

here, ${\varvec{B}}$ is the left Cauchy–Green tensor, $W_{M}^{i}$ is the free energy of the $i$-th network in the MN elastomers with $M$ networks, ${\varvec{B}}_{M}^{i}$ is the left Cauchy–Green tensor of MN elastomers with respect to the relaxed configuration of the $i$-th network. It takes into account swelling, drying, and mechanical deformation. $\phi_{M}^{i}$ is the volume fraction of the $i$-th network.

To capture the progressive damage that allows for the breaking of shorter chains before longer chains, the strain energy of the filler network, i.e., the first network, is assumed to be

$$ W_{M}^{1} = \mu \int_{1}^{\infty } P \left( {N_{j} } \right)d\left( {\lambda_{\max }^{*} ,N_{j} } \right)N_{j} E^{*} \left( {r^{*} } \right){\text{d}}N_{j} $$

(52)

here, $\mu$ is a modulus defined by $\mu = vk_{B} T$ with $v$ being the volumetric density of chains (before any swelling), $N_{j}$ is the number of chain segments, and $P\left( {N_{K} } \right)$ is a probability density function of chain length distribution. $\lambda^{*} = \frac{\lambda }{{\sqrt {N_{j} } }} = \frac{r}{{N_{j} b}}$ is the fractional chain extension, where $r$ is the end-to-end distance of the polymer chain and $b$ is the Kuhn length of the chain segment. $d\left( {\lambda_{\max }^{*} ,N_{K} } \right)$ is the damage function that depends on the chain length $N_{j}$ and the historical maximum fractional extension $\lambda_{\max }^{*}$, and $E^{*} \left( {\lambda^{*} } \right) = E\left( {\lambda^{*} } \right)/k_{B} T$ is the nondimensional chain energy to be determined.

The above model is related to the continuum level through the eight-chain network representation, with the chain extension specified as

$$ \lambda^{*} = \frac{\lambda }{{\sqrt {N_{j} } }} = \sqrt {\frac{{I_{1} \left( {{\varvec{B}}_{M}^{1} } \right)}}{{3N_{j} }}} $$

(53)

Rather than directly defining a specified chain energy, Lavoie et al. [117] introduced a dimensionless chain force $F^{*} \left( {\lambda^{*} } \right) = bF\left( {\lambda^{*} } \right)/k_{B} T$, and the dimensionless chain energy is then calculated from

$$ E^{*} \left( {\lambda^{*} } \right) = \int_{{\lambda_{0}^{*} }}^{{\lambda^{*} }} {F^{*} } \left( {\lambda^{*} } \right){\text{d}}\lambda^{*} $$

(54)

where $\lambda_{0}^{*} = \sqrt {1/N_{j} }$ is the fractional extension of the chain in the stress-free state.

For the poly(ethyl acrylate) (PEA) chain, the force-extension relationship $F^{*} \left( {\lambda^{*} } \right)$ can be approximated by polynomial functions as follows:

$$ \begin{aligned} F^{*} &= \frac{1}{2}\left( {1 - \lambda^{*} } \right)^{ - 2} - \frac{1}{2} - 2\lambda^{*} \quad \left({\lambda^{*} < 0.9}\right), \\ F^{*} &= \frac{513}{{20}} + 501\left( {\lambda^{*} - 0.9} \right) + 26238\left( {\lambda^{*} - 0.9} \right)^{2}\\ &\quad + 68436\left( {\lambda^{*} - 0.9} \right)^{3} \quad \left(\lambda^{*} \ge 0.9\right) \\ \end{aligned} $$

(55)

Moreover, the damage function is chosen as a Weibull cumulative distribution function

$$ d = 1 - \exp \left( { - \left( {\frac{{F_{cr} - F_{\max } }}{{c_{2} }}} \right)^{{c_{1} }} } \right) $$

(56)

where $F_{\max }$ is the maximum chain force related to $\lambda_{\max }^{*}$, the constants $F_{cr}$, $c_{1}$, and $c_{2}$ depend on the material, the chain length, and the chain extension rate.

Finally, the generalized neo-Hookean model is adopted for individual matrix networks:

$$ W_{M}^{i} \left( {I_{1} \left( {{\varvec{B}}_{M}^{i} } \right)} \right) = \frac{{\mu_{i} }}{{2n_{i} }}\left\{ {\left[ {1 + \left( {I_{1} \left( {{\varvec{B}}_{M}^{i} } \right) - 3} \right)} \right]^{{n_{i} }} - 1} \right\},\quad i \ge 2 $$

(57)

where $\mu_{i}$ is the shear modulus and $n_{i}$ is an additional parameter related to the nonlinearity of $W_{M}^{i}$. Lavoie et al. [117] showed that their model can closely match the cyclic uniaxial loading data for the MN elastomers in Ducrot et al. [9].

3.2.8 Tube-Like Chain Model with Entangled Effect

The single-chain models above do not account for any interactions from other surrounding chains, allowing the considered chain to freely rotate in three-dimensional space and pass through neighboring chains. This apparently unsatisfactory assumption is often deemed to be responsible for the poor biaxial performance of many classical hyperelastic models [80].

To address this issue, Doi et al. [118] simulated the topological constraint from neighboring chains through an imagined tube that encloses the polymer chain, and kinematics of the chain is restricted to be inside the tube, as shown in Fig. 8. Hence, the topological constraint, also termed as the entanglement effect, increases as the tube shrinks. The free energy of a single chain is then given by the sum of the cross-linked entropic part $w_{{\text{c}}}$ and the entangled part $w_{{\text{e}}}$:

$$ w = w_{{\text{c}}} + w_{{\text{e}}} $$

(58)

Fig. 8

Tube model of a single chain

Representative tube-like entanglement models can be referred to, for example, Edwards [119], Heinrich et al. [120], Kaliske et al. [121], Miehe et al. [102], and Xiang et al. [122]. This concept has also been advocated in the constitutive modeling of the Mullins effect in polymer networks, as discussed below.

3.2.9 Model of Khiem et al.

Khiem et al. [123] developed a tube model based on the averaging network for deformation-induced stress softening of filled elastomers. The total free energy density is decomposed into three parts:

$$ W = W_{{\text{E}}} + W_{{\text{M}}} + W_{{\text{H}}} $$

(59)

where $W_{{\text{E}}}$ corresponds to the elastic pure rubber network, the M network is in charge of the Mullins effect, and the H network is responsible for the hysteresis.

For the rubber network $W_{{\text{E}}}$, the single-chain free energy is given by the sum of the cross-linked part and the entangled part. The former is specified by a novel form in Khiem et al. [124]:

$$ w_{{\text{c}}} = k_{B} TN\ln \frac{{\sin \frac{\pi }{\sqrt N } \cdot \lambda }}{{\sin \left( {\frac{\pi }{\sqrt N }\lambda } \right)}} $$

(60)

which is derived through a coarse-grained model based on the solution of finitely extensible dumbbells in quantum mechanics [124, 125]. The entangled part is given as

$$ w_{{\text{e}}} = k_{B} T\omega v $$

(61)

where $\omega$ is a geometric parameter of the tube, and $v$ is the contraction ratio of the tube’s cross section.

Given a set of predefined direction specified by unit vectors ${\varvec{E}}_{i} (i = 1,2,3...,m)$, the chain distribution along each predefined direction is assumed to follow the Mises-Arnold-Fisher distribution function:

$$ \rho_{i} \left( {\theta_{i} ,\varsigma_{i} ,{\varvec{E}}_{i} } \right) = \frac{{\varsigma_{i} }}{{4\pi \sinh \left( {\varsigma_{i} } \right)}}\cosh \left( {\varsigma_{i} \cos \theta_{i} } \right) $$

(62)

with $\theta_{i}$ being the angle between the chain direction and the predefined direction ${\varvec{E}}_{i}$ and $\varsigma_{i}$ being the concentration parameter. Then, the average chain stretch $\overline{\lambda }$ and average tube contraction $\overline{v}$ are given by their root mean square, respectively:

$$ \overline{\lambda } = \Lambda^{q} ,\quad \Lambda = \left[ {\int_{S} \rho ({\varvec{n}})\lambda_{{\varvec{n}}}^{2} {\text{d}}S} \right]^{\frac{1}{2}} = \left[ {\sum_{i = 1}^{m} \frac{1}{m} \left( {\left( {1 - w_{i} } \right)\frac{{I_{1} }}{3} + w_{i} {\varvec{C}}:{\varvec{E}}_{i} \otimes {\varvec{E}}_{i} } \right)} \right]^{\frac{1}{2}} $$

(63)

$$ \overline{v} = \left[ {\int_{S} \rho ({\varvec{n}{\varvec{\Pi}}})v_{{\varvec{n}}}^{2} {\text{d}}S} \right]^{\frac{1}{2}} = \left[ {\sum_{i = 1}^{m} \frac{1}{m} \left( {\left( {1 - w_{i} } \right)\frac{{I_{2} }}{3} + w_{i} ({\text{det}}{\varvec{C}}){\varvec{C}}^{\text{T}} :{\varvec{E}}_{i} \otimes {\varvec{E}}_{i} } \right)} \right]^{\frac{1}{2}} $$

(64)

with $w_{i} = 1 + \frac{{3 - 3\coth \left( {\varsigma_{i} } \right)s_{i} }}{{\varsigma_{i}^{2} }}$.

The free energy density of the elastic pure rubber network is then written as

$$ W_{{\text{E}}} = \mu_{c} \psi (N,\overline{\lambda }) + \mu_{t} \overline{v} $$

(65)

where $\mu_{c} = n_{c} k_{B} T$ is the shear modulus, $n_{c}$ is the chain density, $\psi (N,\overline{\lambda }) = N\ln \frac{{\sin \frac{\pi }{\sqrt N } \cdot \overline{\lambda }}}{{\sin \left( {\frac{\pi }{\sqrt N }\overline{\lambda }} \right)}}$, and $\mu_{t} = n_{c} k_{B} T\omega$ is the topological shear modulus.

A progressive damage approach is further adopted to simulate the damage behaviors (Fig. 9). During deformation, the elastic rubber network E induces anisotropy in three principal directions. The damaged network is then considered as a composition of several one-dimensional subnetworks with length distribution dispersed around three principal directions. Each directional subnetwork of polymer chains is considered a series of M polymer chains connecting two adjacent aggregates. The length distribution function is given by

$$ P\left( {N_{j} } \right) = \frac{{\alpha \sin^{2} \left( {\frac{{\pi N_{j} }}{{nN_{m} }}} \right)}}{{\exp \left( {\frac{\alpha }{4\pi }\left[ {nN_{m} \left( {\sin \frac{{2\pi N_{c} }}{{nN_{m} }} - \sin \frac{2\pi }{{nN_{m} }}} \right) + 2\pi \left( {1 - N_{j} } \right)} \right]} \right)}} $$

(66)

where $N_{j} (j = 1,2,3...,m)$ is the chain length, $\alpha$ is the average functionality of the active adsorption site of the aggregate surface, and $N_{m}$ is the maximum chain length. During deformation, subnetworks oriented in the $i$-th principal direction with the chain length $N_{j}$ shorter than $\Lambda_{i}^{\max } L$ will break and no longer contribute to the network entropy. Here, $\Lambda_{{\text{i}}}^{\max }$ is the maximum value of the $i$-th principal stretch, and $L$ is the average referential end-to-end distance of the subnetworks of polymer chains.

Fig. 9

Evolution of the set of available chains under deformation

Hence, the fraction of available chains in one direction can be calculated as

$$\begin{aligned} \Phi_{i} \left( {\Lambda_{i}^{\max } } \right) &= \int_{{\Lambda_{{\text{i}}}^{\max } L}}^{{nN_{m} }} P \left( {N_{j} } \right){\text{d}}N_{j}\\ &= {A}_{{D}} \exp \left( {\frac{{\overline{\alpha }}}{{{4}\pi }}\left[ {n\sin \frac{{2\pi \Lambda_{i}^{\max } \overline{L}}}{n} - 2\pi \Lambda_{i}^{\max } \overline{L}} \right]} \right)\end{aligned} $$

(67)

where $A_{D}$ is a normalization factor, and $\overline{\alpha } = \alpha M_{{\text{c}}}^{\max }$ and $\overline{L} = \frac{L}{{M_{{\text{c}}}^{\max } }}$ are treated as material constants.

Then, the number of available M network chains and the corresponding shear modulus in the $i$-th direction are given respectively by

$$ n_{{{\text{Mi}}}} = \frac{\delta }{3}N_{{\text{p}}} \Phi_{{\text{i}}} \left( {\Lambda_{{\text{i}}}^{\max } } \right),\quad \mu_{{{\text{Mi}}}} = N_{{\text{M}}} k_{{B}} T\kappa \Phi_{{\text{i}}} \left( {\Lambda_{{\text{i}}}^{\max } } \right) $$

(68)

here, $n_{{\text{p}}} = n_{{\text{c}}} \frac{{C_{0} }}{{1 - C_{0} }}$ and $N_{{\text{M}}} = \frac{\delta }{3}\frac{{C_{0} }}{{\left( {1 - C_{0} } \right)}}N_{{\text{C}}}$, with $C_{0}$ being the virgin volume fraction of filler. $N_{{\text{M}}} k_{{B}} T\kappa$ is the cross-linked shear modulus of the M network in the reference configuration.

The fraction of available chains within the H network in the $i$-th direction is given by

$$ \Omega_{{\text{i}}} \left( {\Lambda_{{\text{i}}}^{{{\text{cmax}}}} } \right) = A_{D} \exp \left( {\frac{{\overline{\alpha }}}{4\pi }\left[ {n\sin \frac{{2\pi \Lambda_{{\text{i}}}^{{{\text{cmax}}}} \overline{L}}}{n} - 2\pi \Lambda_{{\text{i}}}^{{{\text{cmax}}}} \overline{L}} \right]} \right) $$

(69)

here, $\Lambda_{{\text{i}}}^{{{\text{cmax}}}}$ is the maximal cyclic stretch in the $i$-th direction defined by $\Lambda_{{\text{i}}}^{{{\text{cmax}}}} = \max_{{\tau \in \left[ {t_{{\text{i}}}^{,} ,t} \right]}} \Lambda_{{\text{i}}} (\tau )$, where the ending time of the closest directional loading cycle $t_{{\text{i}}}^{c}$ is defined by the local minimum of $\Lambda_{{\text{i}}}$ within the time interval $\tau \in [ - \infty ,t]$.

Then, the number of available chains in the H network and the corresponding shear modulus in the $i$-th direction is given as

$$ n_{{{\text{Hi}}}} = \frac{(1 - \delta )}{3}n_{{\text{p}}} \Omega_{{\text{i}}} \left( {\Lambda_{{\text{i}}}^{{{\text{cmax}}}} } \right),\quad \mu_{{{\text{Hi}}}} = n_{{\text{H}}} k_{{B}} T\kappa \Omega_{{\text{i}}} \left( {\Lambda_{{\text{i}}}^{{{\text{cmax}}}} } \right) $$

(70)

with $N_{{\text{H}}} = \frac{(1 - \delta )}{3}\frac{{C_{0} }}{{\left( {1 - C_{0} } \right)}}N_{{\text{C}}}$.

Finally, the topological shear moduli of M and H networks can be given by

$$ \mu_{{{\text{tM}}}} = N_{{\text{M}}} k_{{B}} T\omega \prod\limits_{i = 1}^{3} {\Phi_{{\text{i}}} } \left( {\Lambda_{{\text{i}}}^{\max } } \right),\quad \mu_{{{\text{tH}}}} = N_{{\text{H}}} k_{{B}} T\omega \prod\limits_{i = 1}^{3} {\Omega_{{\text{i}}} } \left( {\Lambda_{{\text{i}}}^{{{\text{cmax}}}} } \right) $$

(71)

Moreover, the kinematics of subnetworks differs from the stretch ratio of the whole sample, which is further specified by

$$ \lambda_{{{\text{Mi}}}} = \Lambda_{{\text{i}}}^{{q/(1 - C_{{{\text{Mi}}}} )}} $$

(72)

where $q$ is the stretch amplification exponent for pure rubber and $C_{{{\text{Mi}}}} = \left[ {C_{0} \Phi_{{\text{i}}} \left( {\Lambda_{{\text{i}}}^{\max } } \right)} \right]^{X}$ is the effective directional fraction of filler with the exponent $X$ depending on the structure of the filler network.

For the H network, the following expressions are used:

$$ \lambda_{{{\text{Mi}}}} = \Lambda_{{\text{i}}}^{{q/(1 - C_{{{\text{Mi}}}} )}} ,\quad C_{{{\text{Mi}}}} = \left[ {C_{0} \Phi_{{\text{i}}} \left( {\Lambda_{{\text{i}}}^{{{\text{cmax}}}} } \right)} \right]^{X} $$

(73)

The tube contraction is given as

$$ v_{{{\text{Mi}}}} = \overline{v}_{{\text{i}}}^{{q_{{{\text{Mi}}}} }} ,\quad v_{{{\text{Hi}}}} = \overline{v}_{{\text{i}}}^{{q_{{{\text{Hi}}}} }} $$

(74)

Then, the total free energy density of the elastomer is written as

$$ \begin{aligned} W & = W_{{\text{E}}} + W_{{\text{M}}} + W_{{\text{H}}} = \mu_{{\text{c}}} \psi_{{\text{c}}} (N,\overline{\lambda }) + \mu_{{\text{t}}} \overline{v} \\ &\quad + \sum_{i = 1}^{3} {\mu_{{{\text{Mi}}}} } \psi_{{\text{c}}} \left( {N,\lambda_{{{\text{Mi}}}} } \right) + \mu_{{{\text{tM}}}} v_{{{\text{Mi}}}} + \mu_{{{\text{Hi}}}} \psi_{{\text{c}}} \left( {n,\lambda_{{{\text{Hi}}}} } \right) + \mu_{{{\text{tH}}}} v_{{{\text{Hi}}}} \\ \end{aligned} $$

(75)

Compared to other damage models, the parameters in the model of Khiem et al. [123] have clear physical motivations, and the model achieves good agreement with multi-dimensional experimental data of soft composites.

3.2.10 Model of Zhong et al.

Zhong et al. [126] developed a damage model for the Mullins effect for soft elastomers based on the tube-like entangled models. In their model, the damage of the entangled part is governed by adding an exponential reduction factor to the three-chain-based entangled model of Xiang et al. [122]:

$$ W_{{\text{e}}} = \frac{{G_{{{\text{e0}}}} }}{{\sqrt {{\text{exp}}\left[ {k\sqrt {I_{1}^{{{\text{max}}}} /3} } \right]} }}\left( {\frac{1}{{\lambda_{1} }} + \frac{1}{{\lambda_{2} }} + \frac{1}{{\lambda_{3} }}} \right) $$

(76)

where $G_{\text{e}}{0}$ is the initial entangled modulus, $k$ is a parameter related to the damage speed, and $I_{1}^{{{\text{max}}}}$ is the historical maximum $I_{1}$.

Different from Xiang et al. [122], in which the cross-linked part is given by the eight-chain model, the cross-linked part in the model of Zhong et al. [126] is characterized by the affine microsphere model:

$$ W_{{\text{c}}} = \frac{1}{S}\int {n_{{\varvec{n}}} } k_{B} TN_{{\varvec{n}}} \left( {\frac{{\lambda_{{\varvec{n}}} \beta }}{\sqrt N } + {\text{ln}}\frac{\beta }{{{\text{sinh}}\beta }}} \right){\text{d}{S}},\quad \lambda_{{\varvec{n}}} { = }\sqrt {{\varvec{nCn}}} $$

(77)

where $S$ is the sphere area, and $n_{{\varvec{n}}}$, $N_{{\varvec{n}}}$, and $\lambda_{{\varvec{n}}}$ are respectively the chain density, chain length, and chain stretch in direction ${\varvec{n}}$.

The network alternation theory [24] is applied to characterize the damage of the cross-linked part:

$$ n_{{\varvec{n}}} = \frac{{n_{0} }}{{p(\lambda_{{\varvec{n}}}^{{{\text{max}}}} - 1)^{2} + 1}},\quad N_{{\varvec{n}}} = N_{0} \left( {p\left( {\lambda_{{\varvec{n}}}^{{{\text{max}}}} - 1} \right)^{2} + 1} \right) $$

(78)

where $n_{0}$ and $N_{0}$ are the initial chain density and chain segment number, respectively, and $\lambda_{{\varvec{n}}}^{{{\text{max}}}}$ is the historical maximum stretch in direction ${\varvec{n}}$. Here, $n_{{\varvec{n}}} N_{{\varvec{n}}} = n_{0} N_{0}$ is assumed to enforce the conservation of chain segments.

The success of the models of Xiang et al. [122] and Zhong et al. [126] is mostly attributed to the three-chain model adopted for the entangled part, which is capable of describing various elastomers with relatively simple formulations. We noticed that the tube model was also adopted in the early work of Goktepe et al. [127] toward modeling the Mullins effect. They used the affine assumption for both the chain stretch and tube shrinkage to simulate the Mullins effect in filled rubbers. In their model, however, the damage is only induced in the particle-to-particle network, and hence, the entangled free energy in the cross-link-to-cross-link network does not decrease. This is different from the work of Zhong et al. [126].

3.2.11 Model of Shen et al.

In fact, the main novelty of the models of Zhong et al. [126] and Xiang et al. [123] lies in the new micro–macro-transition used to relate the single-chain models to the continuum ones. Regarding this issue, Shen et al. [128] also developed a network model with the entanglement effect based on the average deformation of a unit microsphere. The cross-linked and entangled parts of a single chain are given as

$$ \begin{gathered} w_{{\text{c}}} = k_{B} TN\left( {\frac{\lambda }{\sqrt N }\beta + \ln \frac{\beta }{\sinh \beta }} \right) \hfill \\ w_{{\text{e}}} = \alpha k_{B} TN\left( {\frac{b}{{d_{0} }}} \right)^{2} \hfill \\ \end{gathered} $$

(79)

The single-chain model is then extended to the macroscopic level as

$$ W = G_{{\text{c}}} N\left( {\frac{{\overline{\lambda }}}{\sqrt N }\beta + \ln \frac{\beta }{\sinh \beta }} \right) + 3G_{{\text{e}}} \overline{v} $$

(80)

In the above, the average chain stretch $\overline{\lambda }$ and average tube shrink $\overline{v}$ are specified by:

$$ \begin{gathered} \overline{\lambda } = \sqrt {\frac{1}{S}\int {\lambda_{{{\text{affine}}}}^{2} } {\text{d}}S} = \sqrt {\frac{{I_{1} }}{3}} \hfill \\ \overline{v} = \sqrt {\frac{1}{S}\int {v_{{{\text{affine}}}}^{2} } {\text{d}}S} = \sqrt {\frac{{I_{2} }}{3}} \hfill \\ \end{gathered} $$

(81)

where $\lambda_{{{\text{affine}}}}$ and $v_{{{\text{affine}}}}$ are the affine chain stretch and affine area shrink in a specific direction, respectively.

The above network model, with simple forms directly related to the first and second invariants of the Cauchy–Green tensor, is further applied to characterize the Mullins effect in filled rubbers using the network alternation theory [24] with evolving chain length and modulus for the cross-linked part:

$$ \begin{gathered} G_{\text{c}} = G_{{{\text{c0}}}} /\exp \left[ {a_{1} \left( {\sqrt {\frac{{I_{1}^{\max } }}{3}} - 1} \right)} \right] \hfill \\ N = N_{0} *\exp \left[ {a_{1} \left( {\sqrt {\frac{{I_{1}^{\max } }}{3}} - 1} \right)} \right] \hfill \\ \end{gathered} $$

(82)

The damage of entangled part is specified by the following relation:

$$ G_{\text{e}} = \frac{{G_{{{\text{e0}}}} }}{{1 + a_{2} \left( {\sqrt {\frac{{I_{2}^{\max } }}{3}} - 1} \right)}} $$

(83)

The above model can properly describe the stress response of filled rubbers under uniaxial and biaxial loading conditions. It captures the tremendous stress softening behaviors under biaxial loading conditions for rubbers with a large fraction of fillers, which is beyond the capability of traditional damage models. Moreover, this model is easy for applications in commercial finite element software since it only depends on the strain invariants.

3.2.12 Model of Diani et al.

Other new concepts concerning the micro–macro-transition are also adopted for modeling the Mullins effect. We start with the model of Diani et al. [129], which employs an equal-force model and the network alternation to construct a damage model that successfully captures the main features of the Mullins effect in styrene-butadiene rubber (SBR) gum with varying concentrations of carbon black fillers. The equal-force model, originally proposed in Tkachuk et al. [130], is an implicit network model with minimum free energy subject to the kinematic constraint that the macro-deformation gradient is the average of the micro-directional stretch in all directions within a unit microsphere:

$$ \frac{1}{3}{\varvec{F}} = \frac{1}{S}\int {{{\varvec{\lambda}}}_{{\varvec{n}}} } \otimes {\varvec{n}}{\text{d}}S $$

(84)

where ${\varvec{F}}$ is the deformation gradient, ${{\varvec{\uplambda}}}_{{\varvec{n}}}$ is the micro-stretch vector in the direction denoted by unit vector ${\varvec{n}}$, and $S$ is the microsphere surface.

The constitutive relations can be worked out by solving the minimum energy problem:

$$ L({{\varvec{\lambda}}}_{{\varvec{n}}} ,{{\varvec{\Pi}}}) = \frac{1}{S}\int {G_{{\varvec{n}}} } N_{{\varvec{n}}} \left( {\frac{{\lambda_{{\varvec{n}}} \beta }}{\sqrt N } + {\text{ln}}\frac{\beta }{{{\text{sinh}}\beta }}} \right){\text{d}S} - {{\varvec{\Pi}}}:\left( {\frac{{1}}{{\text{S}}}\int {{{\varvec{\lambda}}}_{{\varvec{n}}} } \otimes {\varvec{n}}{\text{ds}} - \frac{{1}}{{3}}{\varvec{F}}} \right) $$

(85)

The Lagrange multiplier is proved to be the first Piola–Kirchhoff stress, namely ${{\varvec{\Pi}}} = {\varvec{P}}$. The above model then indirectly yields an assumption of

$$ f_{{\varvec{n}}} = \frac{\partial w}{{\partial \lambda_{{\varvec{n}}} }} = \sqrt {{\varvec{Pn}} \cdot {\varvec{Pn}}} $$

(86)

where $f_{{\varvec{n}}}$ is the chain force in direction ${\varvec{n}}$. This is in line with the equal-force model proposed by Verron et al. [131], and we hence refer to both models in Tkachuk et al. [130] and Verron et al. [131] as equal-force models.

Toward modeling the stress softening, it is assumed that the number of chain segments $N_{{\varvec{n}}}$ increases with the maximum force $f_{{\varvec{n}}}^{{{\text{max}}}}$ undergone over the loading history in each direction ${\varvec{n}}$:

$$ \left\{ {\begin{array}{*{20}l} {N_{{\varvec{n}}} (t) = \left[ {1 + \alpha \left( {f_{{\varvec{n}}}^{{{\text{max}}}} - f^{{{\text{ini}}}} } \right)} \right]N_{{\varvec{n}}}^{0} } \hfill \\ {\dot{N}_{{\varvec{n}}} = 0{\text{ when }}f_{{\varvec{n}}} < f_{{\varvec{n}}}^{{{\text{max}}}} } \hfill \\ \end{array} } \right. $$

(87)

where $f^{{{\text{ini}}}}$ is the chain force in the undeformed state. This model can reproduce the behavior of SBR gum filled with four levels of carbon black fillers. The equal-force network model has also been applied in Rastak et al. [132] for the strain-crystallizing rubber-like materials, as well as in Li et al. [133] for the brittle fracture of polydisperse elastomers.

3.2.13 Model (b) of Xiao et al.

The equal-force network representation, together with the eight-chain model and the affine model, was applied in Xiao et al. [30] toward modeling the Mullins effect of DN hydrogels under complex multiaxial loading conditions. The damage is assumed to be induced in the first network. Moreover, the stress contribution in the secondary network is ignored since it is significantly lower than that of the first network.

For the eight-chain network model, the free energy and the damage rule are given as

$$ \left\{ {\begin{array}{*{20}l} {W = GN\left[ {\frac{{\overline{\lambda }\beta }}{\sqrt N } + {\text{ln}}\frac{\beta }{{{\text{sinh}}\beta }}} \right],\quad \overline{\lambda } = \sqrt {\frac{{I_{1} }}{3}} } \hfill \\ {G = \frac{{G_{{0}} }}{{\exp \left[ {a_{1} \left( {I_{1}^{\max } - 3} \right)} \right]}},\quad N = N_{0} *\exp \left[ {a_{1} \left( {I_{1}^{\max } - 3} \right)} \right]} \hfill \\ \end{array} } \right. $$

(88)

For the affine microsphere model, the free energy and the damage rule are represented as

$$ \left\{ {\begin{array}{*{20}l} {\frac{1}{S}\int {G_{{\varvec{n}}} } N_{{\varvec{n}}} \left( {\frac{{\lambda_{{\varvec{n}}} \beta }}{{\sqrt N_{{\varvec{n}}} }} + {\text{ln}}\frac{\beta }{{{\text{sinh}}\beta }}} \right){\text{d}S},\quad \lambda_{{\varvec{n}}} { = }\sqrt {{\varvec{nCn}}} } \hfill \\ {G_{{\varvec{n}}} = \frac{{G_{{{\varvec{n}}0}} }}{{1 + a_{2} \left( {\lambda_{{\varvec{n}}}^{\max } - 1} \right)^{2} }},\quad N_{{\varvec{n}}} = N_{{{\varvec{n}}0}} \left( {1 + a_{2} \left( {\lambda_{{\varvec{n}}}^{\max } - 1} \right)^{2} } \right)} \hfill \\ \end{array} } \right. $$

(89)

For the equal-force model, the free energy is given through the following minimum problem:

$$ L(\lambda_{{\varvec{n}}} ,{{\varvec{\Pi}}}) = \frac{1}{S}\int {G_{{\varvec{n}}} } N_{{\varvec{n}}} \left( {\frac{{\lambda_{{\varvec{n}}} \beta }}{{\sqrt N_{{\varvec{n}}} }} + {\text{ln}}\frac{\beta }{{{\text{sinh}}\beta }}} \right){\text{d}}S - {{\varvec{\Pi}}}:\left( {\frac{{1}}{{\text{S}}}\int {{{\varvec{\uplambda}}}_{{\varvec{n}}} } \otimes {\varvec{n}}{\text{d}S - }\frac{{1}}{{3}}{\varvec{F}}} \right) $$

(90)

and the damage rule is specified by

$$ \begin{gathered} G_{{\varvec{n}}} = \frac{{G_{{{\varvec{n}}0}} }}{{1 + a_{3} \left( {\beta_{{\varvec{n}}}^{{{\text{max}}}} - \beta_{{\varvec{n}}}^{{{\text{ini}}}} } \right)^{2} }} \hfill \\ N_{{\varvec{n}}} = N_{{{\varvec{n}}0}} \left( {1 + a_{3} \left( {\beta_{{\varvec{n}}}^{{{\text{max}}}} - \beta_{{\varvec{n}}}^{{{\text{ini}}}} } \right)^{2} } \right) \hfill \\ \end{gathered} $$

(91)

In the above, $a_{1}$, $a_{2}$, and $a_{3}$ are three pre-factors of the damage function, $\beta_{{\varvec{n}}}^{{{\text{ini}}}} = {\mathcal{L}}^{ - 1} \left( {1/\sqrt {N_{0} } } \right)$, and $\beta_{{\varvec{n}}}^{{{\text{max}}}}$ is the maximum $\beta_{{\varvec{n}}}$ over all the deformation history in the ${\varvec{n}}$ direction. Since the chain force of the primary cross-linked network is proportional to $\beta_{j}^{{\text{p}}}$, the last model, i.e., the equal-force model, is essentially a stress-based damage model, which is different from the other two deformation-based ones.

As shown in Xiao et al. [30], the equal-force microsphere damage model can provide more satisfactory results for the multiaxial behaviors of DN hydrogels than the eight-chain model and affine microsphere model. We attribute this success mostly to the stress-based damage criterion used therein. However, the stress in a deformed body is a prior unknown and not detectable by experiments, particularly for incompressible cases. This leads to extra numerical complexities, which may be difficult for further practical implementation.

3.2.14 Model of Zhan et al.

As shown above, the equal-force model is essentially developed through a force-based micro–macro-mapping, which is completely different from traditional network models. Most recently, Zhan et al. [104] proposed a new micro–macro-transition rule that relates the chain stretch to the macro-deformation for the microsphere network representation: $\lambda_{{\varvec{n}}} = {\varvec{nUn}},$ where $\lambda_{{\mathbf{n}}}$ is the micro-chain stretch in direction ${\mathbf{n}}$ and ${\mathbf{U}}$ is the right stretch tensor. Zhan et al. [104] developed the micro–macro-mapping through an average of affine stretch projection along a specified direction. The same micro–macro-transition has also been proposed by Amores et al. [134] using the free fluctuating network assumption. With this new micro–macro-mapping, the full chain model can well predict the multiaxial responses for various hyperelastic materials, along with the parameters pre-calibrated through the uniaxial data.

This new micro–macro-mapping was further applied to develop a damage model in Zhan et al. [135] for the Mullins effect of soft composites. It starts with a one-dimensional damage model:

$$ w = \text{e}^{ - l} w_{0} (\lambda ) $$

(92)

where $l$ is the logarithmic damage variable, and $w_{0} (\lambda )$ is the elastic one-dimensional free energy.

Analogous to the plasticity theory, the one-dimensional damage potential is defined as

$$ F = w - r(l) $$

(93)

and the damage law is derived as

$$ \dot{l} = \frac{m}{{w_{0} + \text{e}^{l} \frac{\partial r}{{\partial l}}}}\frac{{\partial w_{0} }}{\partial \lambda }\dot{\lambda } $$

(94)

through the consistency condition $\dot{F} = 0$. The damage criterion is given by the Kuhn–Tucker condition:

$$ \begin{gathered} {\text{i}}. \, m = 1\qquad {\text{for}}\quad \dot{\lambda } > {0}\quad {\text{and}}\quad F = {0} \hfill \\ {\text{ii}}. \, m = 0\qquad {\text{for}}\quad {\text{othercases}} \hfill \\ \end{gathered} $$

(95)

The above one-dimensional model is extended to yield a three-dimensional multiaxial model through the combination of microsphere representation and the newly developed micro–macro-mapping:

$$ W = \frac{1}{S}\int {\text{e}^{{ - l_{{\varvec{n}}} }} } w_{0} (\lambda_{{\varvec{n}}} ){\text{d}}S,\quad \lambda_{{\varvec{n}}} = {\varvec{nUn}} $$

(96)

where $l_{{\varvec{n}}}$ and $\lambda_{{\varvec{n}}}$ are the one-dimensional logarithmic damage variable and stretch in direction ${\varvec{n}}$, respectively.

The following exponential form is then used for the elastic one-dimensional free energy:

$$ w_{0} (\lambda ) = \frac{G}{u}\text{e}^{u(\lambda - 1)} $$

(97)

with G the modulus and $u$ a parameter related to the nonlinearity. The stored energy limit $r(l)$ is further specified by

$$ r(l) = w_{0} (\lambda_{c} ) + a[1 - {\text{sech}} (kl)] $$

(98)

where $\lambda_{c}$ is the damage threshold stretch, and $a$ is a hardening parameter.

In the above model, with the help of the utilization of a logarithmic damage variable, the thermodynamic conjugate to damage can be defined with clear physical meanings, and the damage evolution rule can be constructed analogously to plasticity theory, which has been a challenge for traditional damage variables. The model demonstrated general capability to describe the Mullins effect in various soft composites. It also exhibits enhanced multiaxial performance compared to models based on affine or eight-chain networks.

A summary of the damage mechanisms and network representations of different models is shown in Table 1.

Table 1 Damage mechanisms and network representations of models for Mullins effect in soft composites

3.3 Phenomenological Approaches

In the current study, we do not intend to provide a detailed introduction to phenomenological approaches since most do not relate to specific damage mechanisms. They can be applied to various soft materials, including traditional filled rubbers and emerging soft composites. This leads to little demand for new specialized phenomenological models beyond the traditional approaches toward modeling emerging soft composites. In the following, we only briefly discuss two classical approaches widely used for phenomenologically simulating the Mullins effect: the continuum damage mechanics approach and the pseudo-elastic approach. Readers can also refer to Diani et al. [25] for more detailed information about these two approaches. Moreover, we note that the basic concepts of these phenomenological approaches have already been adopted in several models discussed above. However, these models are also related to physically based damage mechanisms or network representations. It should also be noted that various physically based models are based on certain assumptions that may not be justifiable. Thus, it can be difficult to distinguish between phenomenological models and physically based models in some cases.

3.3.1 The CDM Approach

The continuum damage mechanics (CDM) approach, pioneered by Kachanov [136], is widely used for hard materials such as metals, concretes, and rocks. Gurtin et al. [137] first introduced it to model the Mullins effect in a one-dimensional case. Simo [21] later proposed a three-dimensional CDM model with a reduction parameter in the free energy function, specified by the form of $\left( {1 - d} \right)W$. This form has been widely adopted by many elastic damage models for filled rubbers [138,139,140]. In general, these well-established continuum frameworks for filled rubbers can also be applied to simulate the Mullins effect in new tough soft composites, as they are flexible to contain various possible damage mechanisms and need not be restricted to any specified physical mechanisms. As mentioned above, many physically based models have also adopted the CDM concept. For example, the one-dimensional damage model of Zhan et al. [135] was essentially constructed through the CDM approach, then combined with the microstructural network to form a “mixed” continuum model. Moreover, the microstructural model developed by Vernerey et al. [100] adopted a damage variable derived through the integration of the chain fracture fraction, in line with the CDM models. The damage model proposed by Talamini et al. [112] is also based on nonlocal CDM approaches, although it starts with physically based chain energy and the eight-chain network model. These models were constructed based on the physical network representations, while the damage feature is simulated through phenomenological CDM approaches.

3.3.2 The Pseudo-Elasticity Theory

Besides the CDM approach, the pseudo-elasticity theory is also widely used for simulating the Mullins effect in soft materials. The concept of pseudo-elasticity was first proposed by Fung [141] toward modeling biological tissues. Ogden et al. [22] extended this concept to simulate the Mullins effect in filled rubbers, which was then followed by many recent studies [91, 142,143,144,145]. In contrast to the CDM models, the pseudo-elasticity theory assumes that damage only evolves in the unloading stage. This key assumption has often been criticized by researchers [25, 139]. For emerging soft composites, Ducrot et al. [9] found that the bond rupture in MN elastomers only occurs in the loading process, with the unloading process being purely elastic. In this sense, the pseudo-elasticity models may not reflect reality, at least for the MN polymers. However, as is discussed in Diani et al. [25], the CDM model intrinsically assumes that the ratio of two unloading stresses corresponding to two specified maximum loads should be constant. This theoretical property does not coincide with experimental observations. In fact, most CDM models cannot well reproduce the unloading stress–strain curve, while the pseudo-elastic models often perform much better. Recently, Tang et al. [90] and Lu et al. [111] have adopted the pseudo-elastic theory together with physically based damage mechanisms to model emerging soft composites, as discussed above.

We shall also mention the ancient idea of the “two-phase model” first proposed by Mullins et al. [13], which was then extended in a series of studies [23, 38, 146, 147] for filled rubbers. In the two-phase model, rubber-like materials are assumed to be composited of hard and soft regions. During loading, the hard regions irreversibly transform into soft regions. In the unloading and reloading prior to the previous maximum loading, the phase fractions remain unchanged. We find that this idea is coincident with MN polymers composed of hard and soft networks. However, the two-phase models have not been applied in this area so far.

4 Modeling Anisotropic Features of the Mullins Effect

Most theoretical models reported in the literature are constructed to describe the isotropic Mullins effect. However, the pioneering work of Mullins [12] has revealed the anisotropic features of stress softening phenomena. The softening effect is much weaker in an orthogonal direction than that in the previous loading direction. Many experiments [18, 20, 34, 36, 148,149,150,151,152,153] have now observed that the Mullins effects in various soft composites are strongly dependent on the loading directions. Such anisotropic features, along with large deformation and nonlinear inelastic behaviors, further complicate theoretical modeling aspects. In this subsection, we will discuss several representative anisotropic models for the Mullins effect reported in the literature.

4.1 Anisotropic Model within the Microsphere Framework

We shall first pay particular attention to the microsphere model developed by Miehe et al. [102] for the hyperelastic behavior of rubbers, in which the polymer chains are assumed to be uniformly distributed in all directions of the unit microsphere element. The macroscopic free energy is then obtained by integrating the average of the directional micro-chain energies. This idea is essentially in line with the full network concept [154, 155]. Throughout this paper, we refer to both as the microsphere framework. For isotropic rubber elasticity, Miehe et al. [102] proposed a new micro-to-macro-transition between the macro-deformation field and the average network stretch based on the p-root averaging operator. In a companion publication [127], the affine microsphere model was further developed for the inelastic Mullins effect, incorporating a damage mechanism of bond breakage between particles and chains, allowing for formulating distinct chain kinematics and damages in each direction. Goktepe et al. [127] demonstrated that the directional damages of the polymer network automatically result in the permanent set and induce anisotropy as macroscopic phenomena. Since then, the microsphere framework has been widely adopted to construct anisotropic models in conjunction with different damage mechanisms [18, 30, 100, 129, 156, 157]. As discussed in Merckel et al. [158], the Mullins effect can be induced as long as any one of the three-dimensional directions is stretched beyond its historical maximum stretch. Such a direction-dependent damage criterion can be easily implemented in microsphere models. In fact, almost all existing microstructural anisotropic models are developed within the microsphere framework. These models address two main issues: determining the appropriate damage mechanism and related damage evolution rules, and relating directional kinematics in the microscale to detectable deformation in the macroscale.

In the microsphere model of Dargazany et al. [18], damage to the polymer-filler network is considered as a consequence of chain sliding on or debonding from aggregates. Damage in each direction is governed by a network evolution concept that describes changes in the inter-aggregate distribution of polymer chains, and micro-directional stretch is related to macro-stretch through a phenomenological amplification function. The model includes several physically motivated material constants and demonstrates good agreement with its own experimental data on subsequent uniaxial tensions in two orthogonal directions. This model was further extended in another study [156] to incorporate the cross-link-particle network, accounting for the damage of filler aggregates in cyclic deformation as the source of hysteresis energy loss. The accuracy of the resulting model is evaluated in comparison to a new set of experimental data.

In the previous section, we have discussed the model of Zhong et al. [126], which was developed based on the entangled network model of Xiang et al. [122]. In fact, Zhong et al. [126] adopted the microsphere model for the cross-linked part to replace the isotropic eight-chain network in Xiang et al. [122], and the directional chain stretch is related to macro-deformation through the affine assumption. Therefore, this model is also capable of simulating the anisotropic feature of the Mullins effect. The model of Xiao et al. [30] is also an anisotropic damage model within the microsphere framework. The damage process is described by the network alteration theory, while the micro–macro-transition is achieved through the equal-force model. This model successfully describes the stress-stretch responses of the DN hydrogels in all loading conditions. Moreover, it provides a detailed discussion of anisotropic damage by analyzing the modulus reduction in different directions. In the studies of Diani et al. [129] and Zhan et al. [157], the microsphere framework is also applied to simulate the Mullins effect with other different damage mechanisms. However, the anisotropic features are not discussed in detail.

In addition to physically based models, the microsphere framework can be combined with other non-physical damage rules to yield phenomenological anisotropic models. For example, Merckel et al. [159] developed a phenomenological damage model within the microsphere framework to simulate the anisotropic Mullins effect of filled rubbers. Distinct from traditional models, the directional elastic free energy and damage rule in [159] are obtained through a computational procedure based on the given experimental data without treating any previously defined constitutive functions and undetermined parameters. The resulting model is shown to accurately reproduce the stress-stretch responses for uniaxial and biaxial loading conditions. Moreover, the model has been tested for several filled rubbers and provides a reasonable estimation of the anisotropic features for the Mullins effects. This damage model is reminiscent of some other parameter-free constitutive approaches, such as the machine-learning-based constitutive modeling [160] and the polynomial interpolation-based hyperelasticity modeling [161, 162]. Other phenomenological models within the microsphere framework can be referred to, for example, Diani et al. [163], Machado et al. [164], and Zhan et al. [135].

It should be noted that the anisotropic microsphere framework is computationally expensive. It needs to discretize the sphere surface with a set of Gaussian points with specified weight factors and orientations. The integration is then approximated by a sum of weighted values of integrands on these Gaussian points. Consequently, such numerical integration schemes induce discontinuous anisotropy in the loaded materials. The reloading along a direction beyond the previous discrete Gaussian points may need further proper approximate treatment. Moreover, large errors in the stress–strain relation of microsphere models may be observed due to the numerical integration over the unit sphere [165]. Hence, a tensor-based anisotropic model may be expected to bypass these numerical difficulties, as will be discussed below.

4.2 Tensor/Principal Stretch-Based Anisotropic Models

In addition to microsphere models, only a few models based on tensor/principal stretch for the anisotropic Mullins effect have been reported. To the best of our knowledge, these models include Shariff [166, 169], Itskov et al. [167], Dorfmann et al. [168], Khiem et al. [123], and Anssari et al. [170].

Shariff [166] developed a direction-dependent model that describes the anisotropic behaviors of the Mullins effect. The model is built based on a symmetric orthotropic damage tensor, in which the three principal values correspond to the damage degree of the principal directions. The damage rules are formulated in each direction, similar to the CDM approach. This model was further extended in another development [169] to incorporate residual strains. The new model is directly constructed on the orthotropic principal axis, while a direction-dependent free energy function, written explicitly in terms of principal stretches, is postulated, with a direction-dependent damage parameter to simulate anisotropic stress softening in rubbers. The models show good consistency with several experiments reported in the literature. Dorfmann et al. [168] extended the pseudo-elasticity model of Ogden et al. [22] to direction-dependent stress softening by using an orthotropic constitutive model. The pseudo-elastic energy is formulated by an isotropic damaged variable and three additive anisotropic ones in the mutually orthogonal preferred directions, specified by the corresponding directional strain invariants, respectively. When active, these variables modify the form of the energy function during the deformation process and, therefore, change the material response. However, this model has not been tested by any experimental results. Itskov et al. [167] also developed an anisotropic model that considers the permanent set. The formulation is based on an anisotropic three-dimensional softening criterion and a scalar damage function, both formulated in terms of the principal stretches. The damage function is used to approximate the stress differences in three principal directions, which are formulated in terms of material parameters depending on the maximal principal stretch. The model is shown to be in accordance with the experimental data concerning the anisotropic features and residual strains. The model developed by Khiem et al. [123], which we have discussed in the previous section, is essentially an anisotropic one with mathematical expressions in the form of three principal stretches. It is derived through an average approach from the microsphere model, thus bypassing the complex numerical integration procedure. The obtained mean values in the three principal directions are considered in the anisotropic permanent and recoverable networks. The Mullins effect is accounted for by debonding of polymer chains from the filler surface in the permanent damage network, while the hysteresis is assumed to result from the deformation-induced recovery of the elastic network.

In fact, all the above principal stretch-based models are orthotropic damage models that treat the damaged sample as an ortho-anisotropic material. The damage laws are formulated in each principal direction, which is somewhat in line with the damage models based on the three-chain network models. However, the three-chain network model does not seem to perform satisfactorily in multiaxial cases [78, 79].

It should be mentioned that although classical continuum damage approaches have enjoyed a great deal of success in modeling the isotropic Mullins effect, they can hardly be extended to anisotropic cases. Most of the existing anisotropic CDM models need to introduce high-order damage tensors as well as the related additional constitutive assumptions (such as the strain equivalence and energy equivalence assumptions [171,172,173,174,175]). However, as has been discussed in Carol et al. [176], such tensor-based anisotropic approaches cannot be easily accommodated in the thermodynamic framework due to the lack of physical meanings of damage tensors’ work conjugate stresses, which constitute the generalized stress space wherein the damage criteria and damage rules are defined. Moreover, the existing anisotropic CDM models are commonly developed for geomaterials or metals with linear elastic relations in the small-strain region. The generalization to soft elastomers with nonlinearity and large deformation may involve exceeding complexities and has not been reported yet, to our knowledge. Recently, Reese et al. [177] developed a general framework for modeling anisotropic damage behavior through structural tensors. However, the theory is mathematically complicated and lacks concrete testing. It appears that the proper description of large deformation anisotropic behavior still remains a challenge.

5 Multiaxial Performance of the Models

Many well-established hyperelastic models for soft polymers, despite their extensive popularity, have very limited predictive abilities for multiaxial behaviors [79, 80]. As most damage models are constructed on the base of these hyperelastic models, they may inherit this potential defect. The complexities in characterizing damage features can exacerbate the situation. However, compared to the vast studies concerning the multiaxial performance of hyperelastic models, few damage models have been tested in multiaxial experiments. This issue cannot be overlooked, as most soft composite are in a complex multiaxial loading state in service.

5.1 Multiaxial Performance of Elastic/Damage Model

The multiaxial performance of a hyperelastic model is known to rely on the proper description of the network structure. Regarding this issue, various network models have been proposed, such as the three-chain model [178], the four-chain model [179], the eight-chain model [180], and the microsphere or full network models [102, 154]. Distinct from other discrete chain models, the microsphere model assumes uniform distribution of chains in three-dimensional directions, which intuitively aligns with the realistic microstructure of isotropic polymer networks. However, the microsphere model with affine assumption [154, 155] is blamed for its poor multiaxial performance in describing rubbers [78, 180]. As discussed in Boyce and Arruda [78], the apparently more reasonable microsphere model is not as successful as the eight-chain model. It is argued that the chains in the direction of maximum principal stretch should no longer stretch as long as they reach their limiting extensibility, and the chains in other directions may, therefore, stretch more than that suggested by the affine assumption, thus compensating for the total macro-deformation. To overcome this remedy, some non-affine assumptions have been proposed for the microsphere model [102, 124, 130]. Nevertheless, many damage models within the microsphere framework are still developed based on the affine assumption [100, 126, 127, 181]. Therefore, the multiaxial performance of these models needs further examination.

So far, most relevant experimental studies are presented based on uniaxial tests [16, 17, 182, 183]. Only a few results concerning the multiaxial features of the Mullins effect have been reported in the literature [18,19,20, 153]. Zhan et al. [135] demonstrated that the damage model within the affine full network framework cannot properly predict the multiaxial data of the DN hydrogels reported in Mai et al. [153], as shown in Fig. 10. The multiaxial performance of the damage model based on the eight-chain network model is even worse (Fig. 11). Compared to the affine model, the stress and damage level are both overestimated by the eight-chain model. In fact, the anisotropic features of the Mullins effect clearly demonstrate the limitation of the eight-chain network model with an isotropic average stretch in all directions. Similar results were also reported in Xiao et al. [30] and Zhan et al. [157]. As shown in Table 1, most existing models are constructed based on the network models with affine assumption or eight-chain representation. The results in Zhan et al. [135, 157] and Xiao et al. [30] clearly demonstrate the possible limitations of these models in describing multiaxial deformations.

Fig. 10

Comparing the measured data and predicted result using the affine microsphere model for multiaxial damage behaviors of DN hydrogels: a uniaxial tension; b equibiaxial tension; c unequibiaxial tension in the e₁ direction; and d unequibiaxial tension in the e₂ direction (reproduced from Zhan et al. [104])

Fig. 11

Comparing the measured stress response and the predicted stress response of the damage model based on the eight-chain network representation for DN hydrogels: a uniaxial tension; b equibiaxial tension; c unequibiaxial tension in the e₁ direction; d unequibiaxial tension in the e₂ direction (reproduced from Zhan et al. [104])

In the well-established elastoplasticity or continuum damage theory, multiaxial models are constructed in the first place with three-dimensional tensors and tensor invariants. Some of these models have been validated for multiaxial loading in small-strain cases [174, 184]. However, the situation becomes more complex for the large deformation nonlinear behavior of soft polymers. Usually, classical tensors cannot be easily related to the macromolecular chain features, which intrinsically dominate macroscopic behaviors. Phenomenological models can be constructed without treating this issue. However, most phenomenological models use empirical constitutive functions to describe nonlinear behaviors in experiments, and then relate the simple experimental results to more general multiaxial loading cases directly using the three-dimensional tensor invariants. Such a procedure, of course, does not seem to demonstrate convincing rationality, and the multiaxial validity cannot be ensured. In fact, the microsphere framework is essentially a rational approach to construct multiaxial models from one-dimensional constitutive laws, and the micro–macro-mapping, i.e., the relation between directional micro-stretch and macro-deformation, plays a key role in the multiaxial performance of the microsphere models [104, 135, 181, 185]. Regarding this issue, Amores et.al. [134] and Zhan et al. [104] formulated a novel micro–macro-transition rule: $\lambda_{{\varvec{n}}} = {\varvec{nUn}}$, where $\lambda_{{\varvec{n}}}$ is the micro-chain stretch in direction n, and U is the right stretch tensor. With this new micro–macro-mapping, the microsphere model can well predict the multiaxial responses for various hyperelastic materials, along with the parameters pre-calibrated through uniaxial data. This new micro–macro-mapping was further applied to develop an anisotropic damage model [135], which also exhibits enhanced multiaxial performance compared to the affine model and eight-chain models, as shown in Fig. 12. In classical models, the micro–macro-mapping is endowed no more significant than defining the micro-chain stretch through macro-deformation. However, within the microsphere framework, it plays a more conspicuous role in extending one-dimensional models to three-dimensional cases.

Fig. 12

Comparing the measured stress response and the predicted stress response of the model of Zhan et al. for DN hydrogels: a uniaxial tension; b equibiaxial tension; c unequibiaxial tension in the e₁ direction; d unequibiaxial tension in the e₂ direction (reproduced from Zhan et al. [104])

5.2 Stress- and Strain-Based Formulations

To further study the multiaxial anisotropic features of the Mullins effect, we consider the results of pure shear tests, in which the sample is stretched in one direction and fixed in a perpendicular direction. According to Mai et al. [153], the stress softening effect of the DN hydrogels in the fixed direction is comparable with the primarily stretched direction, as shown in Fig. 13. Similar results were also reported in Mai et al. [19] for filled rubbers. Since the stretch in a fixed direction is completely suppressed, it appears that stress softening is due to the stress applied in this direction and irrelevant to the stretch. We shall mention that in Harwood et al. [186], a similar stress softening was already observed in both the filled and natural rubbers when they were loaded to the same stress, while previously the stress softening was not found in the natural rubber when applying the same stretch [12]. These results, as discussed in Diani et al. [25], suggest relating the Mullins effect to a stress level rather than a stretch level. This poses a challenge for the widely used maximum-deformation-based damage criterion in the literature. Although constitutive relations can be constructed in either strain or stress spaces for solid materials, a strain-based model is often preferred for soft materials since deformation is a fundamental measurable quantity. In fact, most of the existing hyperelastic models are constructed based on the deformation quantities. One may notice that most successful plasticity and damage theories in the literature are constructed with stress-based criteria, as discussed in Simo and Ju [187]. This fact, as well as the stress sensitivity reported in Harwood et al. [186] and Diani et al. [25], delivers strong motivation to develop stress-based damage models for the Mullins effect of soft composites.

Fig. 13

Mullins effect of the DN hydrogel under pure shear deformation: stress-stretch responses a in the primary stretched direction e₁; b in the fixed direction e₂ (reproduced from Mai et al. [153])

Within the microsphere framework, the widely used affine stretch assumption written as [127, 154, 155]

$$ {{\varvec{\lambda}}}_{{\varvec{n}}} = {\varvec{Fn}} $$

(99)

and the equal-force assumption written as [30, 129, 130]

$$ {\varvec{f}}_{{\varvec{n}}} = {\varvec{Pn}} $$

(100)

jointly form a dual formulation of the micro–macro-mapping in the stretch and stress space. Therefore, it is possible to develop strain- and stress-based damage models within the microsphere framework along with the single-chain behavior related to the macro-responses through either kinematic relations (the stretch mapping) or static relations (the stress mapping). We have tested several stretch-based damage models within the microsphere framework, which are the model of Zhan et al. [135], the affine model, and the eight-chain model. Although these models successfully describe stress-stretch responses in the primary stretched direction, they all fall short in predicting responses in the fixed direction, as depicted in Fig. 14. Both the affine model and the model of Zhan et al. [135] formulate the free energy density and the damage law based on the directional stretch, which actually vanishes in the fixed direction, leading to an underestimation of stress and damage level. On the contrary, the eight-chain model overestimates stress and damage due to the average stretch assumption.

Fig. 14

Comparing the measured stress response and predicted stress response of damage models with different micro–macro-mappings for DN hydrogels under pure shear deformation: the model of Zhan et al. in a the stretched direction e₁ and b fixed direction e₂; the affine model in c the stretched direction e₁ and d fixed direction e₂; the eight-chain model in e the stretched direction e₁ and f fixed direction e₂ (the experimental data are extracted from Mai et al. [153])

Xiao et al. [30] demonstrated that a stress-based microsphere damage model can provide a more satisfactory result for the pure shear test, as depicted in Fig. 15. Despite the higher predicted loading stress, the damage level aligns with the experimental data. We attribute this success mostly to the stress-based damage criterion used therein. However, the stress in a deformed body is a prior unknown and not detectable through experiments. This leads to extra numerical complexities, which may pose challenges for further practical implementation. It should also be noted that the equal-force network model is less accurate in predicting the multiaxial elastic behavior of soft polymers, as demonstrated in Zhan et al. [104]. This can be responsible for overestimating the loading stress of the model of Xiao et al. [30]. Therefore, it may be inferred that a proper description of the multiaxial Mullins effect needs both adequate network models and damage criteria.

Fig. 15

Comparing the measured stress response and predicted stress response of DN hydrogels under pure shear deformation a in the stretched direction e₁; b in the fixed direction e₂ (reproduced from Xiao et al. [30])

6 Conclusion

In this paper, we provide a comprehensive overview of the Mullins effect in emerging soft composites, including nanocomposite gels, double-network hydrogels, and multi-network elastomers. We first revisit the experimental observations for these soft composites. Although all these soft materials have exhibited clear stress softening behavior, the detailed performance still varies with material types, primarily due to inherent differences in damage mechanisms. We then summarize the formulations of various models developed in recent years for soft composites, which are mainly constructed based on new concepts concerning the damage mechanisms and network representations, as summarized in Table 1. We also discuss the performances of several anisotropic models, most of which are developed based on the microsphere framework, although several models have also been constructed based on principle stretches. Finally, we compare the predictions of models to multiaxial deformation data of soft composites. Results reveal that the widely used classical network representation and maximum-deformation-based damage criteria do not provide a satisfactory performance in capturing the directional-damage behaviors. We hope that this review can shed light on the development of novel damage models with improved prediction abilities, particularly for complex loading states. The current study focuses on the rate-independent damage behaviors of tough elastomers and gels. The Mullins effect is also related to the viscous effect in certain material systems. In addition, the damage can also continuously evolve with fatigue loading. For these features, interested readers can refer to, for example, Lin et al. [6], Lu et al. [188], Lin et al. [189], Wang et al. [190], and the references therein.