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

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


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

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 : 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]: with n e being the number of statistical chain segments between neighboring entanglements and T e is the trapping factor (0 < T e < 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 e 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 e = λ 2 µ − 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 e → ∞ 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 λ µ : 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 (λ 2 = λ 3 = λ −1/2 1 and ∂λ 2 = ∂λ 3 = −1/2 λ −3/2 1 ∂λ 1 ) we obtain for the engineering stress in stretching direction: 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 f 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: 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 f ≈ 1.8 for the mass fractal dimension and d w = 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: 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][34][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.
Polymers 2020, 12, x FOR PEER REVIEW 4 of 20 Here, is the particle size, is the cluster size in spatial direction and , is the minimum cluster size which will be calculated later on. The fractal exponents are determined as ≈ 1.8 for the mass fractal dimension and = 3.1 for the anomalous diffusion exponent of CCA-clusters [2]. Note that the effective filler volume fraction Φ > Φ 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: 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][34][35] (comp. also [7]). In the undeformed state it is assumed to be isotropic, i.e., ( ) = ( ) = ( ) ≡ ( ).
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.  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 G A is the elastic modulus and ε A,µ is the strain of the fragile filler clusters in spatial direction µ. These quantities and their dependence on cluster size x µ 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.

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][8][9][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.
Polymers 2020, 12, x FOR PEER REVIEW 5 of 20 more closely, which requires the consideration of elasticity and fracture of filler clusters in stretched elastomers.

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][8][9][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 clustercluster aggregation. In our model two kinds of deformations of filler-filler bonds are considered, bending-twistingand 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 bendingtwisting deformation governs the elasticity while the tension deformation is sensitive for fracture. The total force constant of a cluster of size 0 with NB particles in the backbone reads: with the tension part given by: The bending-twisting part reads: 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 s and k h 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 ξ 0 with N B particles in the backbone reads: with the tension part given by: The bending-twisting part reads: Here, d is the bond length (particle size), is the average squared radius of gyration in direction perpendicular to the force and includes a unit vector → z pointing perpendicular to the connecting vector It scales with the squared cluster size, , 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 By comparing the exponents of Equations (9) and (10) one finds that the force constant k s decreases much more rapidly with cluster size ξ 0 than the force constants k h . Accordingly, Equation (8) implies k ξ ≈ k s 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): This approximation without the tension term can be applied for sufficient large clusters with 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]: Here, we have introduced the stretching ∆l b and force constant k b ≡ Q/d 2 related to stretching of the single bonds. In addition. we used the equilibrium conditions for the force ∆l s /∆l In the next section we will use Equation (12)   length.
In the second case, where ⃗ points parallel to the chain, we have = 1 and ⃗ ⃗ points perpendicular to the chain. This implies ̅ = 0 yielding: As expected, the deformation increases linear with the number of bonds.
(a) (b)  With g = 1/12 being the ratio between the squared radius of gyration and the squared length L 2 = (N B d) 2 of a linear chain. Accordingly, the force constants k s drops with the 3rd power of the length.
In the second case, where → F points parallel to the chain, we have g = 1 and → F × → z points perpendicular to the chain. This implies S ⊥ = 0 yielding: 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 The ratio between the ensemble average of the squared radius of gyration R 2 g and the squared end-to-end distance R 2 = N B d 2 is evaluated as g = 1/6 (see e.g., Chapter 2 in [36]). This implies for the force constant: For the deformation under tension of the 1-dimensional random walk shown in Figure 3b one obtains: 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 ξ of Equation (15).

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,b of the bonds is given by the critical stretch of the bonds ∆l f,b in relation to the bond length: ε m f,b ≡ ∆l f,b /d. This implies for the failure strain of the cluster: 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: ε h 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: 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,µ ε µ = G A ε A,µ (ε µ ). With the scaling relation Equation (11) for the elastic modulus of the clusters this delivers for the cluster strain: Here we have replaced the rubber stress by a relative stress with respect to the minimum strain: 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: Here, s d 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: And The elastic constant and failure strains are denoted byQ 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 v > 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,µ ε µ = G A ε A,µ ε µ , the cluster stress σ A,µ responsible for filler-induced hysteresis is obtained by differentiation of Equation (7) with respect to cluster strain: 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: 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 λ 2 = λ 3 = λ −1/2 and macroscales 1 + ε 2 = 1 + ε 3 = (1 + ε 1 ) −1/2 ; ∂ε 2 = ∂ε 3 = −1/2 (1 + ε 1 ) −3/2 ∂ε 1 , the cluster stress in stretching direction for the upcycle (∂ε 1 /∂t > 0) is obtained as: For the down cycle, the lateral directions contribute to the cluster stress (∂ε 2 /∂t > 0; ∂ε 3 /∂t > 0): 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: 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: Then for uniaxial deformations the total stress reads: 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.

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: with the (frame-independent) first invariant of the left Cauchy-Green tensor: and the (frame-independent) generalized invariant: 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.
The iteration procedure for the evaluation of stresses (compare Equation (21a) or (21b) together with Equation (5)) is then replaced by its tensorial analog: 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].

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 t 90% 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.

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 χ 2 = cycles n y mod,n − y exp, n 2 , where y mod, n and y exp, n represent the nth model and experimentally obtained 1. Piola-Kirchhoff ("engineering") stress data point. The set stress parameter s set,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 w − d f ≈ 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 e 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 0 ≡ x µ /d, the tensile strength of virgin filler-filler bonds s v and the tensile strength of damaged filler-filler bonds s d .

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 2 = 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 2 = 0.995 and R 2 = 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.  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 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 c 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 e between successive entanglements decreases with aging time if G e increases (G e ∼ 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.
Polymers 2020, 12, x FOR PEER REVIEW 14 of 20 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, and , 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, ≡ / , since the segment number between successive entanglements decreases with aging time if increases ( ∼ 1/ ). However, it must be noted that it is difficult to distinguish between the two parameters and in the frame of the DFM, since both act in a similar way. ≡ / . The filler clusters are described by four fitting parameters, i.e., the effective filler volume fraction Φ , the average cluster size ≡ 〈 〉/ , the tensile strength of virgin filler-filler bonds and the tensile strength of damaged filler-filler bonds . An additional significant effect of thermo-oxidative aging is observed for the strength of virgin and damaged filler-filler bonds, and , 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 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 0 ≡ x µ /d, the tensile strength of virgin filler-filler bonds s v 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 v 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 v and s d 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

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. 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.

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

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.