**Micromechanics of Stress-Softening and Hysteresis of Filler Reinforced Elastomers with Applications to Thermo-Oxidative Aging**

#### **Jan Plagge and Manfred Klüppel \***

Deutsches Institut für Kautschuktechnologie e. V., Eupener Str. 33, D-30519 Hannover, Germany; jan.plagge@dikautschuk.de

**\*** Correspondence: manfred.klueppel@dikautschuk.de

Received: 25 May 2020; Accepted: 8 June 2020; Published: 15 June 2020

**Abstract:** A micromechanical concept of filler-induced stress-softening and hysteresis is established that describes the complex quasi-static deformation behavior of filler reinforced rubbers upon repeated stretching with increasing amplitude. It is based on a non-affine tube model of rubber elasticity and a distinct deformation and fracture mechanics of filler clusters in the stress field of the rubber matrix. For the description of the clusters we refer to a three-dimensional generalization of the Kantor–Webman model of flexible chain aggregates with distinct bending–twisting and tension deformation of bonds. The bending–twisting deformation dominates the elasticity of filler clusters in elastomers while the tension deformation is assumed to be mainly responsible for fracture. The cluster mechanics is described in detail in the theoretical section, whereby two different fracture criteria of filler–filler bonds are considered, denoted "monodisperse" and "hierarchical" bond fracture mechanism. Both concepts are compared in the experimental section, where stress–strain cycles of a series of ethylene–propylene–diene rubber (EPDM) composites with various thermo-oxidative aging histories are evaluated. It is found that the "hierarchical" bond fracture mechanism delivers better fits and more stable fitting parameters, though the evolution of fitting parameters with aging time is similar for both models. From the adaptations it is concluded that the crosslinking density remains almost constant, indicating that the sulfur bridges in EPDM networks are mono-sulfidic, and hence, quite stable—even at 130 ◦C aging temperature. The hardening of the composites with increasing aging time is mainly attributed to the relaxation of filler–filler bonds, which results in an increased stiffness and strength of the bonds. Finally, a frame-independent simplified version of the stress-softening model is proposed that allows for an easy implementation into numerical codes for fast FEM simulations

**Keywords:** filled elastomers; stress softening; filler-induced hysteresis; cluster mechanics; FEM simulation

#### **1. Introduction**

Nanoscopic fillers like carbon black or silica play an important role in the mechanical reinforcement of elastomers [1,2]. They make the elastomer stiffer and tougher, leading to a pronounced reduction of crack propagation rates and wear, which is accompanied by an increased life time of rubber goods [3]. But the incorporation of fillers results in a nonlinear dynamic-mechanical response, which is reflected e.g., by the amplitude-dependence of the dynamic moduli. This so called Payne effect was investigated by several authors like Payne [4] and Medalia [5]. A related phenomenon is the stress softening effect under quasi-static cyclic deformation, which is also called Mullins effect due to the intensive studies by Mullins [6]. Accordingly, a drop in stress appears if the loading goes beyond the previous maximum. Most of the stress drop at a certain strain occurs in the first cycle, and in the following cycles the specimen approaches a steady state stress–strain curve. A second characteristic effect caused

by fillers is the pronounced hysteresis which is related to the dissipation of mechanical energy during every cycle.

Based on the investigations of filler network morphology by different experimental techniques, a model of rubber reinforcement by flexible filler clusters has been proposed that allows for a quantitative understanding of the complex mechanical response under quasi-static and dynamic deformations [7–10]. This model of rubber reinforcement refers to the kinetics of cluster–cluster aggregation (CCA) of filler networking in elastomers, which represents a reasonable theoretical basis for analyzing the linear viscoelastic properties of reinforced rubbers. According to this approach, filler networks consist of a space-filling configuration of fractal CCA-clusters with characteristic mass fractal dimension *df* ≈ 1.8. The mechanical response of such filler networks at small strain depends purely on the fractal connectivity of the CCA-clusters. It can be evaluated by referring to the Kantor–Webman model [11] of flexible chains of filler particles that allows for a micromechanical description of the elastic properties of tender CCA-clusters in elastomers. The main contribution of the elastically stored energy in the strained filler clusters results from the bending–twisting deformation of filler–filler bonds, considered by an elastic constant *G*. Based on this approach a power law behavior of the small strain modulus vs. filler concentration can be derived. The predicted exponent 3.5 is in good agreement with experimental data of Payne [4] for carbon black filled butyl rubber. This power law behavior has been confirmed by further investigations of carbon black and silica filled rubbers as well as composites with polymeric model fillers (microgels) [7–10].

Based on this approach, a micromechanical material model of filler reinforced elastomers has been put forward, denoted dynamic flocculation model (DFM) [7,12–16]. Similar models have been proposed by other authors [17–19]. The DFM describes the complex quasi-static stress–strain response of filler reinforced elastomers up to large strains in fair agreement with experimental data. It is based on a non-affine tube model of rubber elasticity and considers hydrodynamic amplification of the rubber matrix by a fraction of hard rigid filler clusters with filler–filler bonds in the unbroken, virgin state. The filler-induced hysteresis is described by the cyclic breakdown and re-aggregation of the residual fraction of more soft filler clusters with already broken, but reformed filler–filler bonds. The difference between hard and soft filler–filler bonds arises due to the softening of glassy-like polymer bridges between adjacent filler particles [7,20–22]. If they break under stress the recovery to a virgin bond will take time and/or high temperatures, as applied during vulcanization [23]. A finite element implementation of the DFM was established [24] by referring to the concept of representative directions introduced by Ihlemann [25]. Thereby also temperature and time dependent effects have been considered [26].

It is the aim of this study to consider the underlying fracture mechanics of filler clusters entering the DFM more closely, whereby we will investigate two different fracture criteria for filler–filler bonds. The differences between both approaches will be evaluated by comparing fitting results obtained for stress–strain cycles of ethylene–propylene–diene monomer rubber (EPDM) composites with various thermo-oxidative aging histories, which will be analyzed by the variation of fitting parameters. This is in close relation to recent investigations of the structure–property relationships of silica/silane formulations in different rubber composites [27] delivering microscopic information about the material properties of the systems. Finally, we will introduce a frame-invariant version of the stress-softening part of the DFM that allows for an easy implementation into a finite element code for fast finite element (FEM) applications of the isotropic discontinuous damage effects in engineering rubber science.

#### **2. Modeling Approach**

#### *2.1. Basic Assumptions of the Dynamic Flocculation Model*

The (microscopic) free energy density of the dynamic flocculation model (DFM) consists of two contributions, which are weighted by the effective filler volume fraction Φeff:

$$\mathcal{W}\{\varepsilon\_{\mu}\} = (1 - \Phi\_{\text{eff}})\mathcal{W}\_{\text{R}}\{\varepsilon\_{\mu}\} + \Phi\_{\text{eff}}\,\mathcal{W}\_{\text{A}}\{\varepsilon\_{\mu}\} \tag{1}$$

The first addend is the equilibrium energy density stored in the extensively strained rubber matrix, which includes hydrodynamic amplification by a fraction of rigid filler clusters. The second addend considers the energy stored in the strained soft filler clusters and is responsible for the filler-induced hysteresis. The symbol εμ is defined in this work as the macroscopic strain in direction μ. The rubber elastic part is modeled by the free energy density of the extended non-affine tube model [12,28]:

$$\mathcal{W}\_{\mathbb{R}}(\varepsilon\_{\mu}) = \frac{G\_{\mathbb{C}}}{2} \left\{ \frac{\left(\sum\_{\mu=1}^{3} \lambda\_{\mu}^{2} - 3\right) \left(1 - \frac{T\_{\text{e}}}{n\_{\text{e}}}\right)}{1 - \frac{T\_{\text{e}}}{n\_{\text{e}}} \left(\sum\_{\mu=1}^{3} \lambda\_{\mu}^{2} - 3\right)} + \ln\left[1 - \frac{T\_{\text{e}}}{n\_{\text{e}}} \left(\sum\_{\mu=1}^{3} \lambda\_{\mu}^{2} - 3\right)\right] \right\} + 2G\_{\mathbb{e}} \left(\sum\_{\mu=1}^{3} \lambda\_{\mu}^{-1} - 3\right) \tag{2}$$

with *n*<sup>e</sup> being the number of statistical chain segments between neighboring entanglements and *T*<sup>e</sup> is the trapping factor (0 < *T*<sup>e</sup> < 1) characterizing the fraction of elastically active entanglements. We define λμ to be the microscopic strain ratio (or stretch) on the nanoscale in direction μ. Thus, for unfilled rubbers the usual relation λμ = 1 + εμ holds. The first addend in Equation (2) considers the constraints due to interchain junctions, with an elastic modulus *G*c proportional to the density of network junctions. The second addend considers topological constraints in densely packed polymer networks, whereby *G*<sup>e</sup> is proportional to the entanglement density of the rubber. The parenthetical expression in the first addend considers the finite chain extensibility of polymer networks by referring to an approach of Edwards and Vilgis [29]. For the limiting case *n*e/*T*<sup>e</sup> = λ<sup>2</sup> <sup>μ</sup> − 3 a singularity is obtained for the free energy density *W*R, indicating the maximum extensibility of the network. This is reached when the chains between successive trapped entanglements are fully stretched out. In the limit *n*<sup>e</sup> → ∞ the original Gaussian formulation of the non-affine tube model, derived by Heinrich et al. [30] for infinite long chains, is recovered.

The presence of tightly bonded (virgin bonds) rigid filler clusters gives rise to hydrodynamic reinforcement of the rubber matrix. This is specified by the strain amplification factor *X* as proposed by Mullins and Tobin [6], which relates the external, macroscopic strain εμ of the sample to the internal, microscopic strain ratio λμ of the rubber matrix, λμ = 1 + *X* εμ,max εμ. For strain amplified rubbers this strain has to be used in the free energy density Equation (2). The microscopic stress of the rubber matrix is then obtained by differentiation with respect to the internal strain λμ:

$$\sigma\_{\mathcal{R},\,\nu} \equiv \frac{\partial \mathcal{W}\_{\mathcal{R}}}{\partial \lambda\_{\mathcal{V}}} = G\_{\mathcal{C}} \sum\_{\mu=1}^{3} \frac{\partial \lambda\_{\mu}}{\partial \lambda\_{\mathcal{V}}} \lambda\_{\mu} \left\{ \frac{\left(1 - \frac{T\_{\mathcal{C}}}{\eta\_{\mathcal{C}}}\right)}{\left(1 - \frac{T\_{\mathcal{C}}}{\eta\_{\mathcal{C}}} \left(\sum\_{\mu=1}^{3} \lambda\_{\mu}^{2} - 3\right)\right)^{2}} - \frac{\frac{T\_{\mathcal{C}}}{\eta\_{\mathcal{C}}}}{1 - \frac{T\_{\mathcal{C}}}{\eta\_{\mathcal{C}}} \left(\sum\_{\mu=1}^{3} \lambda\_{\mu}^{2} - 3\right)} \right\} - \mathcal{L}G\_{\mathcal{C}} \sum\_{\mu=1}^{3} \frac{\partial \lambda\_{\mu}}{\partial \lambda\_{\mathcal{V}}} \lambda\_{\mu}^{-2} \tag{3}$$

This is the microscopic stress between the filler clusters that can be identified with the macroscopically measured engineering stress (1. PK stress) in equilibrium. For uniaxial deformations (λ<sup>2</sup> = λ<sup>3</sup> = λ−1/2 <sup>1</sup> and ∂λ<sup>2</sup> <sup>=</sup> ∂λ<sup>3</sup> <sup>=</sup> <sup>−</sup>1/2 <sup>λ</sup>−3/2 <sup>1</sup> ∂λ1) we obtain for the engineering stress in stretching direction:

$$\sigma\_{\rm R, \,1} = G\_{\rm c} \{ \lambda\_1 - \lambda\_1^{-2} \} \left| \frac{\left(1 - \frac{T\_{\rm e}}{n\_{\rm e}}\right)}{\left(1 - \frac{T\_{\rm e}}{n\_{\rm e}} (\lambda\_1^2 + 2/\lambda\_1 - 3)\right)^2} - \frac{\frac{T\_{\rm e}}{n\_{\rm e}}}{1 - \frac{T\_{\rm e}}{n\_{\rm e}} (\lambda\_1^2 + 2/\lambda\_1 - 3)} \right| + 2G\_{\rm e} \{ \lambda\_1^{-\frac{1}{2}} - \lambda\_1^{-2} \} \tag{4}$$

In the case of a preconditioned sample and for strains smaller than the previous straining (εμ < εμ,max), the materials microscopic structure is already adjusted to the maximum load and the strain amplification factor *X* is independent of strain. In that case it is determined by εμ,max (*X* = *X*(εμ,max)). We relate this to the irreversible fracture of filler clusters (see below). A relation for the strain amplification factor of overlapping fractal clusters of size ξ was derived by Huber and Vilgis [31]. By using path integral methods they found *X* = 1 + *c*Φ2/(3−*d*f) ξ*d*w−*d*<sup>f</sup> where Φ is the filler volume fraction and *c* is a constant of order one. With this, *X*(εμ,max) can be evaluated by averaging over the size distribution of hard clusters in all space directions. In the case of preconditioned samples this yields:

$$X\{\varepsilon\_{\mu,\max}\} = 1 + c\Phi\_{\text{eff}}^{\frac{2}{3-d\_{\text{f}}}} \sum\_{\mu=1}^{3} \frac{1}{d} \left[ \int\_{0}^{\xi\_{\mu,\min}} \left(\frac{\xi\_{\mu}'}{d}\right)^{d\_{\text{W}} - d\_{\text{f}}} \varphi(\xi\_{\mu}') d\xi\_{\mu}' + \int\_{\xi\_{\mu,\min}}^{\infty} \varphi(\xi\_{\mu}') d\xi\_{\mu}' \right] \tag{5}$$

Here, *d* is the particle size, ξμ is the cluster size in spatial direction μ and ξμ,min is the minimum cluster size which will be calculated later on. The fractal exponents are determined as *d*<sup>f</sup> ≈ 1.8 for the mass fractal dimension and *d*<sup>w</sup> = 3.1 for the anomalous diffusion exponent of CCA-clusters [2]. Note that the effective filler volume fraction Φeff > Φ is used in Equations (1) and (5), which considers the effective volume of the rigid phase of structured filler particles, e.g., carbon black or silica, according to the "occluded rubber concept" of Medalia [32]. Occluded rubber is defined as the rubber part of the rubber matrix that penetrates into the voids of the particles, which partially shields it from deformation. The second addend of Equation (5) takes into account that also fully broken clusters contribute to the strain amplification factor by the remaining particles. ϕ(ξμ) is the normalized cluster size distribution:

$$\varphi(\xi\_{\mu}) = 4d \frac{\xi\_{\mu}}{\langle \xi\_{\mu} \rangle^2} \exp \left( -2 \frac{\xi\_{\mu}}{\langle \xi\_{\mu} \rangle} \right) \tag{6}$$

This is a peaked cluster size distribution with ξμ being the ensemble average in spatial direction μ. It is motivated by analytical results referring to Smoluchowski's equation for the kinetics of cluster–cluster aggregation of colloids [33–35] (comp. also [7]). In the undeformed state it is assumed to be isotropic, i.e., ϕ(ξ1) = ϕ(ξ2) = ϕ(ξ3) ≡ ϕ(ξ).

The model of stress softening and hysteresis assumes that the breakdown of filler clusters during the first deformation of the virgin samples is reversible, though the initial virgin state of filler–filler bonds is not recovered. This implies that, on the one side, the fraction of hard (virgin) filler clusters decreases with increasing pre-strain, leading to pronounced stress softening after the first deformation cycle. On the other side, the fraction of soft (reaggregated) filler clusters increases with rising pre-strain, which affects the filler-induced hysteresis. A schematic view of the decomposition of filler clusters in hard and soft units for preconditioned samples is shown in Figure 1.

**Figure 1.** (**a**) Schematic view of the decomposition of filler clusters in hard and soft units for preconditioned samples and (**b**) cluster size distribution with the pre-strain dependent boundary size *x*min between hard and soft clusters.

The second addend of Equation (1) describes the filler-induced hysteresis. It considers the energy stored in the substantially strained soft filler clusters, which break under stress and reaggregate on retraction

$$\mathcal{W}\_{\mathcal{A}}(\varepsilon\_{\mu}) = \sum\_{\mu}^{\partial \varepsilon / \partial t > 0} \frac{1}{2d} \int\_{\xi\_{\mu} \text{min}}^{\xi\_{\mu}(\varepsilon\_{\mu})} \mathcal{G}\_{\mathcal{A}}(\xi\_{\mu}') \varepsilon\_{\mathcal{A},\mu}^{2}(\xi\_{\mu}', \varepsilon\_{\mu}) \phi(\xi\_{\mu}') \, \text{d}\xi\_{\mu}' \tag{7}$$

*G*<sup>A</sup> is the elastic modulus and εA,<sup>μ</sup> is the strain of the fragile filler clusters in spatial direction μ. These quantities and their dependence on cluster size *x*<sup>μ</sup> and external strain εμ will be specified in the next sections. In addition. The integral boundaries of Equations (5) and (7) have to be described more closely, which requires the consideration of elasticity and fracture of filler clusters in stretched elastomers.

#### *2.2. Elasticity and Fracture of Filler Clusters in Stretched Elastomers*

For consideration of filler network breakdown in stretched rubbers, the elasticity and failure properties of tender filler clusters have to be evaluated in dependence of cluster size. This will be obtained by referring to the two-dimensional Kantor–Webman model of flexible chains of arbitrary connected filler particles [11] as represented in Figure 2a. We apply here a simplified generalization of this model to three dimensions, where on-plane bending, and off-plain twisting deformations of bonds are considered by a single bending–twisting term [20]. By identifying the three-dimensional flexible chain with the backbone of a CCA-cluster, the model can be applied for modeling the small-strain modulus of fractal filler networks, consisting of a space-filling configuration of CCA-clusters [7–10]. Here, we use it for a micromechanical description of CCA-clusters that are deformed in the stress field of a strained rubber matrix. Note that this is possible because the CCA-cluster backbone is not branched on large length scales, which is a typical result of cluster–cluster aggregation.

**Figure 2.** Schematic view of the Kantor–Webman model of flexible chains of arbitrary connected filler particles (**a**) and mechanical equivalence between a filler cluster and a series of soft and stiff molecular springs (**b**). The two springs with force constants *k*<sup>s</sup> and *k*<sup>h</sup> are related to bending–twisting- and tension deformations of filler–filler bonds referring to elastic constants *G* and *Q*, respectively.

In our model two kinds of deformations of filler–filler bonds are considered, bending–twisting- and tension deformations. This corresponds to a mechanical equivalence between a filler cluster and a series of two molecular springs depicted schematically in Figure 2b. We will see that the bending–twisting deformation governs the elasticity while the tension deformation is sensitive for fracture. The total force constant of a cluster of size ξ<sup>0</sup> with *NB* particles in the backbone reads:

$$k\_{\xi} = \left(\frac{1}{k\_{\mathrm{h}}} + \frac{1}{k\_{\mathrm{s}}}\right)^{-1} \tag{8}$$

with the tension part given by:

$$k\_{\rm h} = \frac{Q}{\sum\_{i=1}^{N\_{\rm B}} \left(\overrightarrow{\vec{F}} \cdot \overrightarrow{\vec{b}}\_{i}\right)^{2}} = \frac{Q}{N\_{\rm B}d^{2} \oint\_{\left(\overrightarrow{\vec{F}} \cdot \overrightarrow{\vec{b}}\right)^{2}} \frac{Q}{d} \bigg| \frac{d}{\xi\_{0}}\bigg|^{d/B}\tag{9}$$

The bending–twisting part reads:

$$k\_{\rm s} = \frac{\overline{G}}{\sum\_{i=1}^{N\_{\rm B}} \left[ \left( \frac{\overrightarrow{r}}{\overline{r}} \times \overrightarrow{z} \right) \overrightarrow{(R\_{i-1} - \overrightarrow{R}\_{N\_{\rm B}}) \right]^2} = \frac{\overline{G}}{N\_{\rm B} S\_{\perp}^2} = \frac{\overline{G}}{d^2 g'} \left( \frac{d}{\xi\_0} \right)^{2 + d\_{f, \rm B}} \tag{10}$$

Here, *<sup>d</sup>* is the bond length (particle size), <sup>→</sup> *<sup>F</sup>* is the force and *<sup>g</sup>* <sup>≡</sup> %  <sup>→</sup> *F F* · → *b d* 2 d*S* is the average projection

of bond vectors <sup>→</sup> *b <sup>i</sup>* on the direction of the force (0 < *g* < 1). *S*<sup>2</sup> <sup>⊥</sup> <sup>≡</sup> <sup>1</sup> *N*<sup>B</sup> *N*<sup>B</sup> *i*=1 & <sup>→</sup> *F <sup>F</sup>* <sup>×</sup> <sup>→</sup> *z* → *Ri*−<sup>1</sup> − → *RN*<sup>B</sup> '2 is the average squared radius of gyration in direction perpendicular to the force and includes a unit vector → *<sup>z</sup>* pointing perpendicular to the connecting vector <sup>→</sup> *Ri*−<sup>1</sup> − → *RN*<sup>B</sup> . It scales with the squared cluster size, *S*2 <sup>⊥</sup> <sup>=</sup> *<sup>g</sup>* ξ2 <sup>0</sup>, with a scaling factor 0 < *g* < 1. *G* and *Q* are elastic constants due to bending–twisting—and tension deformations, respectively. For the particle number the scaling relation *N*<sup>B</sup> = (ξ0/*d*) *<sup>d</sup>*f,B was used with *d*f,B = 1.3 being the backbone fractal dimension of CCA-clusters.

By comparing the exponents of Equations (9) and (10) one finds that the force constant *k*s decreases much more rapidly with cluster size ξ<sup>0</sup> than the force constants *k*h. Accordingly, Equation (8) implies *k*<sup>ξ</sup> ≈ *k*<sup>s</sup> for sufficient large clusters, i.e., the stiffness of the cluster is determined by the bending–twisting deformations of bonds. This determines the following scaling law for the elastic modulus entering Equation (7):

$$G\_{\mathcal{A}}(\xi\_0) \approx \xi\_0^{-1} k\_{\mathcal{S}} = \frac{\overline{G}}{d^3 g'} \left(\frac{d}{\xi\_0}\right)^{3+d\_{\mathcal{S}B}} \simeq \frac{\overline{G}}{d^3} \left(\frac{d}{\xi\_0}\right)^{3+d\_{\mathcal{S}B}} \tag{11}$$

This approximation without the tension term can be applied for sufficient large clusters with (ξ0/*d*) <sup>2</sup>+*d*f,B (ξ0/*d*) *<sup>d</sup>*f,B . For the evaluation of the scaling factor *g* in Equation (11) we have to consider the ensemble average of clusters. However, this will not be considered here, because we are mainly interested in the scaling exponents.

The stretching of the clusters can be evaluated in the same approximation [7]:

$$
\Delta l\_{\xi} = \Delta l\_{\text{s}} + \Delta l\_{\text{h}} = \Delta l\_{\text{b}} \left( \frac{k\_{\text{b}}}{k\_{\text{s}} \left( \xi\_{\text{0}} \right)} + gN\_{\text{B}} \right) \approx \Delta l\_{\text{b}} \frac{g'Q}{\overline{G}} \left( \frac{\xi\_{\text{0}}}{d} \right)^{2 + d\_{\text{fB}}} \tag{12}
$$

Here, we have introduced the stretching <sup>Δ</sup>*l*<sup>b</sup> and force constant *<sup>k</sup>*<sup>b</sup> <sup>≡</sup> *<sup>Q</sup>*/*d*<sup>2</sup> related to stretching of the single bonds. In addition. we used the equilibrium conditions for the force Δ*l*s/Δ*l*<sup>b</sup> = *k*b/*k*s(ξ0) and Δ*l*h/Δ*l*<sup>b</sup> = *k*b/*k*h(ξ0) = *gN*B. In the next section we will use Equation (12) for describing the fracture of filler clusters by relating it to the fracture of bonds under tension.

Examples:

The cluster mechanics described by Equations (8)–(10) shall be illustrated by two simple examples which are depicted in Figure 3. In Figure 3a a linear chain is considered with the force <sup>→</sup> *F* pointing perpendicular or alternatively in direction of the chain. In the first case we have *<sup>g</sup>* <sup>=</sup> 0 and <sup>→</sup> *<sup>F</sup>* <sup>×</sup> <sup>→</sup> *z* points into the direction of the chain. This implies for the force constant and the deformation:

$$k\_{\xi} = k\_{\mathfrak{s}} = 12 \frac{\overline{G}}{N\_{\mathrm{B}} \left( N\_{\mathrm{B}} d \right)^{2}} \rightarrow \ \Delta \overline{l}\_{s} = \frac{\overrightarrow{F}}{k\_{\mathfrak{s}}} = \frac{\overrightarrow{F}}{12 d \overline{G}} \left( N\_{\mathrm{B}} d \right)^{3} \tag{13}$$

**Figure 3.** Illustration of the cluster mechanics by simple examples: (**a**) Linear chain and (**b**) one dimensional random walk.

With *g* = 1/12 being the ratio between the squared radius of gyration and the squared length *L*<sup>2</sup> = (*N*B*d*) <sup>2</sup> of a linear chain. Accordingly, the force constants *k*<sup>s</sup> drops with the 3rd power of the length.

In the second case, where <sup>→</sup> *<sup>F</sup>* points parallel to the chain, we have *<sup>g</sup>* <sup>=</sup> 1 and <sup>→</sup> *<sup>F</sup>* <sup>×</sup> <sup>→</sup> *z* points perpendicular to the chain. This implies *S*<sup>⊥</sup> = 0 yielding:

$$k\_{\xi} = k\_{\mathrm{h}} = \frac{\mathcal{Q}}{N\_{\mathrm{B}}d^{2}} \rightarrow \Delta \stackrel{\rightarrow}{l}\_{\mathrm{h}} = \frac{\stackrel{\rightarrow}{F}}{k\_{\mathrm{h}}} = N\_{\mathrm{B}} \frac{d^{2}\stackrel{\rightarrow}{F}}{Q} \tag{14}$$

As expected, the deformation increases linear with the number of bonds.

In Figure 3b we consider the case where a random walk structure of the chain is realized, corresponding to three 1-dimensional random walks with *N*B/3 particles. The average projection *g* is then given by *<sup>g</sup>* <sup>=</sup> %  <sup>→</sup> *F F* · → *b d* 2 *dS* = 1/3. The ratio between the ensemble average of the squared radius of gyration *<sup>R</sup>*<sup>2</sup> <sup>g</sup> and the squared end-to-end distance *<sup>R</sup>*<sup>2</sup> = *<sup>N</sup>*B*d*<sup>2</sup> is evaluated as *<sup>g</sup>* = 1/6 (see e.g., Chapter 2 in [36]). This implies for the force constant:

$$k\_{\xi} = \left(\frac{N\_B^2 d^2}{6\overline{G}} + \frac{N\_B d^2}{3Q}\right)^{-1} \tag{15}$$

For the deformation under tension of the 1-dimensional random walk shown in Figure 3b one obtains:

$$
\Delta \stackrel{\rightarrow}{l}\_{\text{h}} = \frac{\stackrel{\rightarrow}{F}}{k\_{\text{h}}} = \frac{N\_{\text{B}}}{3} \frac{d^2 \stackrel{\rightarrow}{F}}{Q} \tag{16}
$$

For the more general case that the force points in arbitrary direction also the bending–twisting deformations of bonds must be taken into account by referring to the full force constant *k*<sup>ξ</sup> of Equation (15).

#### *2.3. Evaluation of Boundary Cluster Size and Cluster Stress*

In view of introducing a fracture criterion for strained clusters, we assume that the tension of bonds is a much more critical deformation compared to bending and twisting, since it separates the filler particles from each other. Equation (12) relates the total stretching of a cluster to the stretching of the bonds and can therefore be used for evaluation of the failure strain of the cluster by defining a fracture criterion for the bonds. We will introduce here two different fracture criteria, which will be denoted "monodisperse" and "hierarchical".

In the first approach all bonds are considered to be equal (monodisperse) having the same strength. Then, the failure strain εf,bof the bonds is given by the critical stretch of the bonds Δ*l*f,b in relation to the bond length: ε*<sup>m</sup>* f,b ≡ Δ*l*f,b/*d*. This implies for the failure strain of the cluster:

$$
\varepsilon\_{\mathbf{f},\xi} \equiv \frac{\Delta l\_{\mathbf{f},\xi}}{\xi\_0} \approx \varepsilon\_{\mathbf{f},\mathbf{b}}^m \frac{\mathbf{g}' \mathbf{Q}}{\overline{\mathbf{G}}} \left(\frac{\xi\_0}{d}\right)^{1+d\_{\mathbf{f},\mathbf{b}}} \tag{17a}
$$

In contrast, the hierarchical model takes into account that a hierarchy of bond strengths develops during cluster–cluster aggregation, because the mobility of the clusters decreases with cluster size. Accordingly, the first bond formed between two particles is the strongest while successive bonds formed between the growing sub-clusters become weaker and weaker. The last bond formed in the cluster is the weakest and will break first under tension. This effect is taken into account by the hierarchical fracture criterion, where the failure strain εf,b of the bonds is defined in relation to the cluster size, which is the only relevant length scale in our model: ε<sup>h</sup> f,b ≡ Δ*l*f,b/ξ0. This implies that the failure strain of the cluster increases more rapidly with cluster size compared to the monodisperse case:

$$
\varepsilon\_{\mathbf{f},\xi} \equiv \frac{\Delta l\_{\mathbf{f},\xi}}{\xi\_0} \approx \varepsilon\_{\mathbf{f},\mathbf{b}}^h \frac{\mathbf{g}' \mathbf{Q}}{\overline{\mathbf{G}}} \left(\frac{\xi\_0}{d}\right)^{2+d\_{\mathbf{f},\mathbf{b}}} \tag{17b}
$$

For the evaluation of the boundary cluster size between broken and unbroken clusters in stretched rubbers, we assume that a stress equilibrium is realized between the strain amplified rubber matrix and the clusters σR,<sup>μ</sup> εμ = *G*AεA,μ(εμ). With the scaling relation Equation (11) for the elastic modulus of the clusters this delivers for the cluster strain:

$$\varepsilon\_{\rm A,\mu}(\varepsilon\_{\mu}) = G\_{\rm A}^{-1} \hat{\sigma}\_{\rm R,\mu}(\varepsilon\_{\mu}) \approx \frac{g'd^3}{\overline{G}} \left(\frac{\xi\_0}{d}\right)^{\mathfrak{Z}+d\_{\rm fB}} \hat{\sigma}\_{\rm R,\mu}(\varepsilon\_{\mu}) \tag{18}$$

Here we have replaced the rubber stress by a relative stress with respect to the minimum strain:

$$\left(\sigma\_{\mathcal{R},\mu}\left(\varepsilon\_{\mu}\right) \equiv \left(\sigma\_{\mathcal{R},\mu}\left(\varepsilon\_{\mu}\right) - \sigma\_{\mathcal{R},\mu}\left(\varepsilon\_{\mu,\min}\right)\right) \tag{19}$$

This ensures that the stretching of clusters in spatial direction μ starts at the minimum strain εμ,min for each cycle. Here, we assume that clusters reaggregate into a stress-free state at minimum strain.

A comparison of the exponents in Equations (18) and (17) makes clear that the strain of the clusters under external strain increases faster with cluster size than the failure strain, in both cases. This implies that large clusters break first followed by smaller ones, i.e., the boundary cluster size between broken and unbroken clusters ξμ εμ moves from larger to smaller values with increasing strain. It is obtained by equating the cluster strain to the failure strain. This yields for the two fracture criteria:

$$\left(\frac{\mathcal{E}\_{\mu}(\varepsilon\_{\mu})}{d}\right)^{2} = \frac{Q\varepsilon\_{\mathrm{f,b}}^{\mathrm{m}}}{d^{3}\hat{\sigma}\_{\mathrm{R,\,\mu}}(\varepsilon\_{\mu})} \equiv \frac{s\_{\mathrm{d}}}{\hat{\sigma}\_{\mathrm{R,\,\mu}}(\varepsilon\_{\mu})}\tag{20a}$$

And

$$\frac{\mathcal{E}\_{\mu}(\varepsilon\_{\mu})}{d} = \frac{Q \varepsilon\_{\mathrm{f,b}}^{\mathrm{h}}}{d^{3} \hat{\sigma}\_{\mathrm{R,\,\,\mu}}(\varepsilon\_{\mu})} \equiv \frac{s\_{\mathrm{d}}}{\hat{\sigma}\_{\mathrm{R,\,\,\mu}}(\varepsilon\_{\mu})} \tag{20b}$$

Here, *s*<sup>d</sup> is defines as the fracture stress under tension of bonds, i.e., the tensile strength of damaged filler–filler bonds. The boundary cluster size ξμ εμ applies for the integral boundaries of Equation (7), describing the filler-induced hysteresis due to the successive breakdown of soft filler clusters with damaged filler–filler bonds.

Similar expressions are found for the upper boundary of Equation (5), but now the tensile strength *s*v of virgin filler–filler bonds is entering:

$$\left(\frac{\xi\_{\mu,\text{min}}}{d}\right)^2 = \frac{\vec{Q}\,\hat{\varepsilon}\_{\text{f,b}}^{\text{m}}}{d^3 \hat{\sigma}\_{\text{R, }\mu}\left(\varepsilon\_{\mu,\text{max}}\right)} \equiv \frac{s\_{\text{V}}}{\hat{\sigma}\_{\text{R, }\mu}\left(\varepsilon\_{\mu,\text{max}}\right)}\tag{21a}$$

And

$$\frac{\mathcal{E}\_{\mu,\text{min}}}{d} = \frac{\vec{Q}\varepsilon\_{\text{f,b}}^{\text{th}}}{d^3 \mathfrak{d}\_{\text{R, }\mu}(\varepsilon\_{\mu,\text{max}})} \equiv \frac{s\_{\text{V}}}{\hat{\mathfrak{o}}\_{\text{R, }\mu}(\varepsilon\_{\mu,\text{max}})} \tag{21b}$$

*Polymers* **2020** , *12* , 1350

The elastic constant and failure strains are denoted by *Q*˜ and ε˜f,b, respectively. Note that the tensile strength of virgin filler–filler bonds must be larger than the tensile strength of damaged bonds, i.e., *s*<sup>v</sup> > *s*d. Equation (21a) or (21b) together with Equation (5) define the amplification of the rubber matrix, and thus stress—softening effects of the model. Solving this set of equations for σ requires iterative methods, e.g., Newton iteration.

By referring to the stress equilibrium between the strain amplified rubber matrix and the clusters, σˆR,<sup>μ</sup> εμ = *GA*εA,<sup>μ</sup> εμ , the cluster stress σA,<sup>μ</sup> responsible for filler-induced hysteresis is obtained by differentiation of Equation (7) with respect to cluster strain:

$$\begin{split} \boldsymbol{\sigma}\_{\mathbf{A},\boldsymbol{\nu}} &= \frac{\frac{\partial \boldsymbol{\varepsilon}}{\partial \boldsymbol{\varepsilon}} > 0}{\frac{1}{\mu} \boldsymbol{\varepsilon}\_{\boldsymbol{\Lambda},\boldsymbol{\mu}}} \frac{1}{d} \int\_{d(\frac{\boldsymbol{\varepsilon}\_{d}}{\boldsymbol{\delta}\_{\mathbf{R},\boldsymbol{\mu}}(\boldsymbol{\varepsilon}\_{\boldsymbol{\mu}})})^{\boldsymbol{\alpha}}} \, \boldsymbol{G}\_{\mathbf{A}} \Big{(}\boldsymbol{\xi}\_{\boldsymbol{\mu}}\big{)} \boldsymbol{\varepsilon}\_{\boldsymbol{\Lambda},\boldsymbol{\mu}} \Big{(}\boldsymbol{\xi}\_{\boldsymbol{\mu}}\big{)} \boldsymbol{\varepsilon}\_{\boldsymbol{\Lambda},\boldsymbol{\mu}} \Big{(}\boldsymbol{\xi}\_{\boldsymbol{\mu}}\big{)} \boldsymbol{\varepsilon}\_{\boldsymbol{\mu}}\big{)} \boldsymbol{\varepsilon}\_{\boldsymbol{\mu}} \Big{)} \boldsymbol{\varepsilon}\_{\boldsymbol{\mu}} \Big{(}\boldsymbol{\xi}\_{\boldsymbol{\mu}}\big{)} \, \mathrm{d}\boldsymbol{\xi}\_{\boldsymbol{\mu}} = \\ & \overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{\partial}}{\overset{\scriptstyle{$$

The exponent α takes the two fracture criteria into account, i.e., α = 1/2 for the "monodisperse" model and α = 1 for the "hierarchical" model. In addition. we assume that the clusters, on average, deform like the sample:

$$
\langle \frac{\partial \varepsilon\_{\rm A, \mu}}{\partial \varepsilon\_{\rm A, \nu}} \rangle = \frac{\partial \varepsilon\_{\mu}}{\partial \varepsilon\_{\nu}} \tag{23}
$$

The sum in Equation (22) runs over stretching directions, only, implying that the up- and down cycles are different. The cluster stress of the upcycle is positive while the downcycle gives a negative contribution, producing the filler-induced hysteresis.

For uniaxial deformations, realized on microscales λ<sup>2</sup> = λ<sup>3</sup> = λ−1/2 <sup>1</sup> ; ∂λ<sup>2</sup> <sup>=</sup> ∂λ<sup>3</sup> <sup>=</sup> <sup>−</sup>1/2 <sup>λ</sup>−3/2 <sup>1</sup> ∂λ1) and macroscales 1 + ε<sup>2</sup> = 1 + ε<sup>3</sup> = (1 + ε1) <sup>−</sup>1/2; ∂ε<sup>2</sup> <sup>=</sup> ∂ε<sup>3</sup> <sup>=</sup> <sup>−</sup>1/2 (<sup>1</sup> <sup>+</sup> <sup>ε</sup>1) <sup>−</sup>3/2∂ε1, the cluster stress in stretching direction for the upcycle (∂ε1/∂*t* > 0) is obtained as:

$$
\sigma\_{\rm A,1}^{\rm up} = \hat{\sigma}\_{\rm R,1}(\varepsilon\_1) \frac{1}{d} \int\_{d(\frac{s\_d}{\beta\_{\rm R,1}(\varepsilon\_1, \max)})}^{d(\frac{s\_d}{\beta\_{\rm R,1}(\varepsilon\_1)})} \, \_{\rm a}\rho(\xi\_1) \, \mathrm{d}\xi\_1 \tag{24}
$$

For the down cycle, the lateral directions contribute to the cluster stress (∂ε2/∂*t* > 0; ∂ε3/∂*t* > 0):

$$\sigma\_{\rm A,1}^{\rm down} = 2\hat{\sigma}\_{\rm R,2}(\varepsilon\_1) \Big( -\frac{1}{2} (1 + \varepsilon\_1)^{-\frac{3}{2}} \Big) \frac{1}{d} \int\_{d(\frac{s\_d}{\delta\_{\rm R,2}(\varepsilon\_1)})}^{d(\frac{s\_d}{\delta\_{\rm R,2}(\varepsilon\_1)})^a} \rho(\xi\_2) \,\mathrm{d}\xi\_2 \tag{25}$$

With ϕ(ξ1) = ϕ(ξ2) = ϕ(ξ3) ≡ ϕ(ξ). This gives a negative stress contribution, which must be subtracted from the rubber stress. It can also be expressed by the rubber stress σR,1 in stretching direction. By assuming that the same energy is needed for stretching in 1-direction and compressing in 2- and 3-direction to obtain a final deformed state, the following relation is derived:

$$
\sigma\_{\rm R,2}(\varepsilon\_2) \equiv \sigma\_{\rm R,2}(\varepsilon\_2) - \sigma\_{\rm R,2}(\varepsilon\_{2,\min}) = -\lambda\_1^{3/2}\sigma\_{\rm R,1}(\varepsilon\_1) - \lambda\_{1,\max}^{3/2}\sigma\_{\rm R,1}(\varepsilon\_{1,\max})\tag{26}
$$

Finally, for the evaluation of the (measured) total stress we have to consider an additional set stress σset that appears as a remaining stress in the undeformed state after stretching and retraction. Note that this is also found for unfilled rubbers and depends on temperature and stretching rate. It probably results from long time relaxation effects of the polymer network. We introduce it in a purely empirical manner for the case of uniaxial deformations:

$$
\sigma\_{\text{set},1} = s\_{\text{set},0} \Big( \sqrt{\varepsilon\_{1,\text{max}}} - \sqrt{\varepsilon\_{1,\text{min}}} \Big) \tag{27}
$$

Then for uniaxial deformations the total stress reads:

$$
\sigma\_{\text{bot},1}(\varepsilon\_1) = \sigma\_{\text{R,1}}(\varepsilon\_1) + \sigma\_{\text{A,1}}^{\text{up/down}}(\varepsilon\_1) + \sigma\_{\text{set},1}(\varepsilon\_{1,\text{max}}, \varepsilon\_{1,\text{min}}) \tag{28}
$$

The stress of the rubber matrix σR,1 is given by Equation (4) with λμ = 1 + *X* εμ,max εμ and strain amplification factor *X* εμ,max specified by Equation (5). The cluster stress σA,1 depends on the direction of straining and is determined by Equations (24) and (25) for up and down, respectively. The theory presented here describes the complex quasi-static deformation behavior of filler reinforced elastomers for repeated stretching with increasing amplitude. A more general formulation of the DFM that applies for arbitrary deformation histories requires an additional term for the free energy density of soft filler clusters Equation (7), which considers the relaxation of cluster stress upon retraction [15,16]. For a test of the theory the present formulation for repeated stretching with increasing amplitude is sufficient because all open parameters are entering already and the extension to arbitrary deformation histories requires no additional fitting parameters.

#### *2.4. Frame-Independent Formulation of Stress-Softening for Fast FEM Simulations*

The implementation of the DFM into a FEM algorithm faces several problems. First, the DFM is a microscopic theory, where stresses are calculated by differentiation of the free energy density with respect to the internal strain variables λμ and λ*A*,μ, respectively. This is in discrepancy to continuum mechanical considerations, where corresponding differentiations have to be performed with respect to external strain variables. Second, the DFM is formulated in the main axis system and requires stress contributions of different directions, especially for description of filler-induced hysteresis, which all sum up to produce close cycles. This can hardly be transferred to a pure tensorial formulation. Nevertheless, a FEM implementation of the DFM was obtained by referring to the concept of representative directions, which considers uniaxial deformations along fibers in different spatial directions [24,25]. However, the computational cost of this workaround is very large, and the efficiency is low. Therefore, we want to focus here on a frame-independent tensor formulation of the stress-softening part of the DFM for fast and efficient FEM simulations.

In a first step we put the strain amplification factor *X*max ≡ *X*(εμ,max) in front of the deformation invariants, appearing in the free energy density of the extended non-affine tube model Equation (2) and replace the internal strain λμ = 1 + *X*(εμ,max)εμ by the external strain λμ = 1 + εμ. This follows the ideas of Einstein [37] and Domurath et al. [38] and is done in close correlation to the evaluations in [22]. The free energy density reads:

$$\mathcal{W}\_{\rm R} = \frac{\mathcal{G}\_{\rm c}}{2} \left\{ \frac{X\_{\rm max} \tilde{I}\_1 \left( 1 - \frac{T\_{\rm e}}{n\_{\rm e}} \right)}{1 - \frac{T\_{\rm e}}{n\_{\rm e}} X\_{\rm max} \tilde{I}\_1} + \ln \left[ 1 - \frac{T\_{\rm e}}{n\_{\rm e}} X\_{\rm max} \tilde{I}\_1 \right] \right\} + 2 \mathcal{G}\_{\rm e} X\_{\rm max} \tilde{I} \tag{29}$$

with the (frame-independent) first invariant of the left Cauchy–Green tensor:

$$
\overline{I}\_1 \equiv \lambda\_1^2 + \lambda\_2^2 + \lambda\_3^2 - 3 \tag{30}
$$

and the (frame-independent) generalized invariant:

$$
\overrightarrow{I} \equiv \lambda\_1^{-1} + \lambda\_2^{-1} + \lambda\_3^{-3} - 3 \tag{31}
$$

Here, λμ = 1 + εμ is the external strain of the sample. A frame-independent formulation of the strain amplification factor *X*max is obtained similar to Equation (5) by replacing the relative stress σˆR,μ(εμ,max) used for the calculation of the boundary cluster size *x*μ, min by the Frobenius norm σR(εmax) of the engineering stress.

$$X(\varepsilon\_{\text{max}}) = 1 + c\Phi\_{\text{eff}}^{\frac{2}{3-d\_{\text{f}}}} \frac{1}{d} \left[ \int\_{0}^{d(s\_{\text{v}}/\sigma\_{\text{R}}(\varepsilon\_{\text{max}}))^{a}} \left( \frac{\xi}{d} \right)^{d\_{\text{w}} - d\_{\text{f}}} \varphi(\xi) d\xi + \int\_{d(s\_{\text{v}}/\sigma\_{\text{R}}(\varepsilon\_{\text{max}}))^{a}}^{\infty} \varphi(\xi) d\xi \right] \tag{32}$$

The iteration procedure for the evaluation of stresses (compare Equation (21a) or (21b) together with Equation (5)) is then replaced by its tensorial analog:

$$\sigma\_{\mathbb{R},\mu}(\varepsilon\_{\mu,\max}) = f(\mathbf{X}\_{\max}(\sigma\_{\mathbb{R},\mu}(\varepsilon\_{\mu,\max}))) \to \sigma\_{\mathbb{R}}(\varepsilon\_{\max}) = f(\mathbf{X}\_{\max}(\sigma\_{\mathbb{R}}(\varepsilon\_{\max}))) \tag{33}$$

The free energy density Equation (29) can be used in a standard continuum mechanical sense for the evaluation of stresses and tangent vectors. It can be further simplified by omitting the logarithmic term, which gives a minor contribution to the stress upturn. In addition. The generalized invariant can be approximated by the square root of the second invariant, which avoids the calculation of eigenvalues [39].

#### **3. Experimental and Fitting Procedure**

#### *3.1. Materials*

EPDM/A ("aging") samples were prepared by using an amorphous-type ethylene–propylene-diene rubber (EPDM, Keltan 4450) filled with 50-phr (parts per hundred mass parts rubber) carbon black (N339). The rubber has a Mooney viscosity of 46 MU (at 125 ◦C) and consists of 52% ethylene and 4.3% ethylidene norbornene. Moreover, 3 phr zinc oxide, 1 phr stearic acid and 1.5 phr *N*-isopropyl-*N*'-phenyl-1,4-phenylenediamine (IPPD, aging protection) were added. The curing system consists of 1.8-PHR sulfur, 1.5 phr 1,3-diphenylguanidine (DPG) and 2.4 phr *N*-cyclohexyl-2-benzothiazolylsulfenamide (CBS, accelerator). EPDM/CB ("carbon black") consists of the same polymer, additives and curing system, but has varying amounts of carbon black (N339): 20, 40 and 60 phr.

#### *3.2. Mixing and Sample Preparation*

All ingredients except of the vulcanization system were mixed in an internal mixer of type GK 1,5E (Werner & Pfleiderer Gummitechnik GmbH, Freudenberg, Germany) at a loading of 75% and a temperature of less than 140 ◦C. The vulcanization system was added in a second step on a 150\*350RR roller mill (KraussMaffei Berstorff GmbH, Hannover, Germany).

All samples were cured in a heated press of type WLP 63/3,5/3 (Wickert Maschinenbau GmbH, Landau, Germany) at 160 ◦C up to the t90% time, where 90% of the torque obtained from a vulcameter measurement is reached (17:23 min). One minute curing time was added per millimeter sample thickness to account for heat diffusion.

From EPDM/A dumbbell test specimen were prepared. These were stored under thermo-oxidative aging conditions at 130 ◦C in an air-ventilated oven for 0, 1, 3, 7 and 14 days.

#### *3.3. Test Methods and Fitting*

Quasi-static multi-hysteresis experiments (multiple deformation cycles up to different strain levels) were performed in a Zwick 1445 (Zwick Roell, Ulm, Germany) universal testing machine using a crosshead speed of 20 mm/min. Every strain level was repeated five times. From EPDM/CB tensile test specimen of S2 type were prepared to achieve strains, which are not accessible using dumbbell samples. Multi-hysteresis experiments at 100 mm/min crosshead speed were performed in the same stretching machine. In both cases, the 5th cycles were separated for fittings with the DFM, which can be considered as equilibrium cycles.

The adaptation of stress-strain cycles to Equation (28) was performed by minimization of the error functional χ<sup>2</sup> = cycles *n y*mod,*<sup>n</sup>* − *y*exp, *<sup>n</sup>* 2 , where *y*mod, *<sup>n</sup>* and *y*exp, *<sup>n</sup>* represent the *n*th model and experimentally obtained 1. Piola–Kirchhoff ("engineering") stress data point. The set stress parameter sset,0 was included in the fitting procedure as defined by referring to Equation (27). The remaining two stress contributions of Equation (28), which involve seven fitting parameters, were obtained iteratively by using Equations (4)–(6) for the intrinsic stress of the rubber matrix σR,1 and Equations (24) and (25) for the evaluation of cluster stress σA,1. Note that both stress contributions involve the same cluster size distribution, which stabilizes the fitting procedure, significantly. The front factor of Equation (5) was fixed as *c* = 2.5 according to the Einstein equation [37] and the exponent was approximated as *d*<sup>w</sup> − *d*<sup>f</sup> ≈ 1 to allow for an analytical solution of the integral. The three fitting parameters describing the rubber elastic network are the crosslink modulus *G*c, the topological constraint modulus *G*<sup>e</sup> and the effective chain length for finite extensibility *n*eff ≡ *n*e/*T*e. The filler clusters are described by four fitting parameters, i.e., the effective filler volume fraction Φeff, the average cluster size *x*<sup>0</sup> ≡ *x*<sup>μ</sup>/*d*, the tensile strength of virgin filler–filler bonds *s*<sup>v</sup> and the tensile strength of damaged filler–filler bonds *s*d.

#### **4. Results and Discussion**

#### *4.1. Micromechanical Investigations of Thermo-Oxidative Aging*

The thermo-oxidative aging of elastomers is of high technological interest for the rubber industry, because it mostly increases the hardness of the samples and has a negative effect on the fracture toughness or crack resistance. This limits the life time of rubber goods, significantly. The reason for this property losses are mainly seen in a change of the rubber–elastic network due to post-curing, but the often-accompanied increase in electrical conductivity indicates that also the carbon black network is altering during thermo-oxidative aging of elastomers. This is hardly to distinguish by standard measurement techniques since changes of the polymer- or filler network structure are difficult to detect. We therefore refer here to the evaluation of micromechanical material parameters, which are obtained by fitting the stress–strain response of the aged samples to the DFM.

Figure 4 shows a series of fittings of multi-hysteresis stress–strain cycles of EPDM/A samples stored under thermo-oxidative aging conditions at 130 ◦C for 0, 1, 3, 7 and 14 days. In Figure 4a, the "monodisperse" bond fracture model (α = 1/2) is used while Figure 4b refers to the "hierarchical" bond fracture model (α = 1). Before we discuss the effect of thermomechanical aging on material parameters, we will first consider the quality of the fits for the two different bond fracture criteria. The "monodisperse" model depicted in Figure 4a delivers fair agreement between fits and experimental data with correlation coefficients between *R*<sup>2</sup> = 0.994 and 0.995, but systematic deviations are seen, e.g., for the peak stresses in the medium strain regime. For the "hierarchical" model shown in Figure 4b, the fits are significantly better with correlation coefficients between *R*<sup>2</sup> = 0.995 and *R*<sup>2</sup> = 0.998. This indicates that the "hierarchical" model is more suited for describing the stress–strain cycles of filler reinforced elastomers. Even in the case of aged samples we get excellent adaptations in the small and medium strain regime up to 100% strain. Therefore, we will mainly focus on the "hierarchical" model with bond fracture exponent α = 1 for the discussion of thermo-oxidative aging effects on a microscopic level. Nevertheless, we will see that the "monodisperse" model delivers similar values and trends of the fitting parameters.

**Figure 4.** Fit of stress–strain cycles of ethylene–propylene-diene rubber EPDM/A samples for various aging times (**a**) with the "monodisperse" bond fracture model (α = 1/2) and (**b**) with the "hierarchical" bond fracture model (α = 1).

*Polymers* **2020** , *12* , 1350

(**a**)

The stress–strain data in Figure 4 show that the average stress level increases with increasing aging time. The reason for this hardening of the samples shall be analyzed by referring to the evolution of fitting parameters that are depicted in Figure 5. The fitting parameters obtained with the "monodisperse" bond fracture model (α = 1/2) are shown in Figure 5a and those from the "hierarchical" bond fracture model (α = 1) in Figure 5b. Obviously, the "hierarchical" bond fracture model in Figure 5b delivers a smoother evolution of fitting parameters, which correlates with the higher correlation coefficients of the fits in Figure 4b. This confirms our view that the "hierarchical" bond fracture model is more suited for the discussion of aging effects. Looking first at the crosslink and topological constraint moduli of the polymer network, *G*<sup>c</sup> and *G*e, we see that the former remains almost constant while the later increases successively with aging time. This indicates that the post-curing effect is not pronounced for the EPDM samples used in this study, but the topological constraints of the chains increase with aging time possibly due to an increasing number of entanglements close to the carbon black particles (surface-induced entanglements). This correlates with the observed decrease of the effective chain length, *n*eff ≡ *n*e/*T*e, since the segment number *n*<sup>e</sup> between successive entanglements decreases with aging time if *G*<sup>e</sup> increases (*G*<sup>e</sup> ∼ 1/*n*e). However, it must be noted that it is difficult to distinguish between the two parameters *G*c and *G*e in the frame of the DFM, since both act in a similar way.

**Figure 5.** Evolution of fitting parameters obtained with (**a**) the "monodisperse" bond fracture model (α = 1/2) and (**b**) the "hierarchical" bond fracture model (α = 1) of the EPDM/A samples from Figure 4 for various aging times.

*n*eff ≡ *n*e/*T*e. The filler clusters are described by four fitting parameters, i.e., the effective filler volume fraction Φeff, the average cluster size *x*<sup>0</sup> ≡ *x*<sup>μ</sup>/*d*, the tensile strength of virgin filler–filler bonds *s*<sup>v</sup> and the tensile strength of damaged filler–filler bonds *s*d.

An additional significant effect of thermo-oxidative aging is observed for the strength of virgin and damaged filler–filler bonds, *s*<sup>v</sup> and *s*d, which both increase systematically with aging time. This is clearly seen in the case of the "hierarchical" bond fracture model (α = 1). It indicates that a relaxation of bonds takes place during heating of the samples at 130 ◦C, leading to more stable filler–filler joints. Note that this is related to confined polymer between the filler particles, which is assumed to be in a glassy-like state at room temperature [7] but can relax at elevated temperature. This relaxation process has been investigated recently by online dielectric spectroscopy during heat treatment of carbon black filled EPDM [40]. The increase of *s*<sup>v</sup> and *s*<sup>d</sup> with aging time observed in Figure 5b results in stronger stress-softening and hysteresis effects of the aged samples. Two further parameters describing the aging of the filler network are the effective filler volume fraction Φeff and the mean cluster size *x*0. The former decreases slightly with aging time, but remains larger than the real filler volume fraction, as expected (Φeff > Φ ≈ 0.2). This indicates a slight decrease of the occluded rubber during aging, which is hidden in the voids of the filler particles and acts like additional filler. The mean cluster size lies in a reasonable range of about 10 to 15 particle diameters and goes through a weak maximum with increasing aging time. This can be related to flocculation effects and restructuring of the filler network

#### *4.2. Frame-Independent Model of Stress Softening*

For the discussion of the frame-independent model of stress softening we will focus on the "hierarchical" bond fracture model with α = 1, only. For the analysis of the stress-softening effect, described here, we refer to the EPDM/CB samples with varying amount of carbon black. Note that the hysteresis is not included in this model.

Figure 6a,b shows fits of stress–strain cycles of the EPDM/CB samples with the frame-independent model (α = 1) up to 100% and 200%, respectively. In both cases the stress-softening effect is reproduced fairly well, though for the 200% fit systematic deviations are seen for the sample with 60-PHR N339 in the small strain regime. The effect of filler concentration on the fitting parameters is depicted for both cases in Figure 7a,b, respectively. All values of the parameters are found in a reasonable range. However, they strongly depend on the range of the fitted stress–strain data up to 100% and 200%, respectively. This indicates that the DFM cannot simply be extended to strains up to 200% because additional mechanisms of stress softening may appear at larger strains, e.g., detachment of the polymer from the filler surface. Nevertheless, the quite reasonable fits in Figure 6b show that the simplified DFM can be used as an empirical model also for larger strains up to rupture of the samples. Due to the simplifications of the frame independent model compared to the original DFM it makes no sense to discuss the dependence of fitting parameters on filler concentration in detail. The more interesting point is the high stability of the fitting procedure with parameters that all are positive and can easily be reproduced. Most important is the ability to implement the simplified model into a finite element code for fast FEM simulations by using standard methods. This is of major interest for the rubber industry and will be a task of future work.

(**a**)

**Figure 6.** Fit of stress–strain cycles of EPDM/CB samples filled with various amounts of carbon black (N339) with the frame-independent model of stress softening (α = 1) data up to 100% strain and (**b**) up to 200% strain.

**Figure 7.** Fitting parameters obtained with the frame independent model of stress softening of the EPDM/2.4-PHR *N*-cyclohexyl-2-benzothiazolylsulfenamide (CBS) samples as shown in Figure 6, using (**a**) data up to 100% strain and (**b**) up to 200% strain.

(**b**)

#### **5. Conclusions**

A micromechanical model of stress-softening and hysteresis of filler reinforced elastomers was presented, which is based on a non-affine tube model of rubber elasticity and a generalized three-dimensional Kantor–Webman model of flexible chain aggregates, describing the deformation and fracture of filler clusters in the stress field of the rubber matrix. This dynamic flocculation model (DFM) has been shown to reproduce the complex quasi-static deformation behavior of filler reinforced elastomers upon repeated stretching with increasing amplitude fairly well. It is described in some detail in the theoretical section, whereby two different fracture mechanisms of filler–filler bonds, denoted "monodisperse" and "hierarchical" bond fracture mechanism, are considered. In the first approach all bonds are considered to be equal (monodisperse) having the same strength. In the second approach a hierarchy of bond strengths is realized, because during cluster–cluster aggregation the mobility of the clusters decreases with cluster size.

In the experimental section the DFM is adapted to a series of aged EPDM samples which were treated in an oven at 130 ◦C for different thermo-oxidative aging times. The fitting parameters indicate that the crosslinking density remains almost constant while the entanglement density increases slightly. This indicates that the sulfur bridges in EPDM networks are quite stable even at 130 ◦C aging temperature, which can be related to the mono-sulfidic nature of the crosslinks. The observed hardening of the composites with increasing aging time is mainly attributed to the relaxation of filler–filler bonds. This results in an increased strength of the bonds, which produces a larger stiffness and filler-induced hysteresis of the composites.

The two different bond fracture mechanisms are investigated by separate adaptations to the full series of aged EPDM composites. They show that the "hierarchical" bond fracture mechanism delivers better fits and more stable fitting parameters, though the evolution of fitting parameters with aging time is similar for both models. Therefore, it is concluded that the "hierarchical" bond fracture mechanism, which takes into account that the mobility of clusters decreases with cluster size, appears to be realized in filler reinforced elastomers.

In the last section a frame-independent simplified version of the DFM is proposed that focuses on an easy implementation of the stress-softening effect of filled rubbers into a finite element codes by using standard methods. The model is shown to reproduce the stress-softening effect of EPDM samples with varying amount of carbon black fairly well. Therefore, it appears to be well suited for performing fast FEM simulations of highly filled rubber goods, where stress-softening cannot be neglected.

**Author Contributions:** J.P. developed the curve fitting algorithm, made the graphs and provided the ideas for the frame-independent DFM. M.K. wrote the major part of the introduction, theory and results section. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work received no special funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
