# Modeling the Payne effect in filled rubbers: accounting for time effect and partial reversibility

Rebeca Cedeno Madera, Stéphane Lejeunes, Eric Euchler

► To cite this version:

Rebeca Cedeno Madera, Stéphane Lejeunes, Eric Euchler. Modeling the Payne effect in filled rubbers: accounting for time effect and partial reversibility. Polymer Bulletin, In press, 10.1007/s00289-025-05680-z. hal-04964987

HAL Id: hal-04964987

https://hal.science/hal-04964987v1

Submitted on 25 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Modeling the Payne effect in filled rubbers: accounting for time effect and partial reversibility

Rebeca Cedeno Madera<fnref n="1" /> · Stéphane Lejeunes<fnref n="1" /> · Eric Euchler<fnref n="2" />

Received: 15 November 2024 / Revised: 4 February 2025 / Accepted: 10 February 2025
© The Author(s) 2025

## Abstract

Carbon black-filled rubbers have been used in the industry since the early 20th century, and their mechanical advantages are well known. However, describing and predicting their mechanical behavior can become a non-trivial task. In this work, we present an approach to a new mechanical model that is capable of capturing the principal nonlinearities in rubber behavior when subjected to dynamic loadings, focusing on the Payne effect. Inspired by phase-change mechanisms, we introduce internal scalar variables to account for the carbon filler network dynamics in the presence of dynamic stretch. Reversibility of the Payne effect is taken into account within a clear and standard thermodynamical framework. This approach aims to enhance our understanding of the mechanical responses of carbon-filled rubbers and improve predictive capabilities.

**Keywords** Payne effect · Filled elastomers · Damage · Reversibility · Preload

## Introduction

Technical elastomers used in the industry are mostly reinforced vulcanized polymers with rigid particles (fillers) to enhance their mechanical properties. These fillers are added in the polymer matrix during the compounding and before the vulcanization. The final mechanical properties depend on many parameters that can be related to

<figure><img src="image_3.png" /></figure>

Stéphane Lejeunes
lejeunes@lma.cnrs-mrs.fr

Rebeca Cedeno Madera
cedeno@lma.cnrs-mrs.fr

Eric Euchler
euchler-eric@ipfdd.de

<footnote n="1">¹ LMA, Aix-Marseille Univ, CNRS, UPR 7051, Centrale Méditerranée, 13402 Marseille Cedex 20, France</footnote>

<footnote n="2">² Leibniz-Institut für Polymerforschung Dresden e.V., 01069 Dresden, Germany</footnote>
the process (mixing conditions, temperature and duration of curing, pressure, etc.) to the recipe (volume ratio of each component, etc.) and the ingredients (type of fillers and matrix, vulcanization agents, etc.). Furthermore, the number of scales, in time and space, involved in the mechanical response makes it difficult to precisely relate experimental observations to a simple mechanism within the material even with advanced experimental techniques (see [1, 2] for a more complete overview). Nevertheless, it is well admitted from the literature that the strong nonlinear viscoelastic behavior of filled rubber can be related to filler-filler or filler-matrix interactions.

Among the specific effects that these materials can exhibit, Mullins effect [3], thermal softening, etc., we are particularly interested in the Payne effect [4]. For CB-reinforced elastomers, many (indirect) observations in the literature converge to a scenario of cyclic breakdown and re-aggregation of filler-filler interactions at least at small strains (e.g., [5–7]). The Payne effect is a dynamical softening effect that is partially reversible, due to the possible re-aggregation of filler-filler bonds, sensitive to the temperature and to the preload in non-relaxed configuration (see for instance [8]). With Fig. 1, we propose a naive schematization of the probable mechanisms. In the undeformed configuration, the fillers (carbon blacks) are organized into networks or agglomerates. A thin layer of glassy polymer encapsulates the fillers and constitutes rigid links between fillers that progressively breaks or becomes rubbery while increasing the dynamic strain (leading to a softening). The opposite phenomena takes place when the dynamics strain decreases (but not necessarily on the same time scale as observed by [8]). A partial recovery can also take place under constant stretch when stress relaxation allows for more chain mobility (also observed by [8] for instance).

From the modeling point of view, this scenario of cyclic breakdown and re-aggregation of filler-filler bonds led to the development of the Krauss model (e.g., [9]). This model furnishes an expression for the storage and loss modulus in a power-law form of the applied amplitude by considering an equilibrium between the rate of breakage and reformation. Further Krauss-like models have been proposed, and interested readers can have a look at [10] for a more detailed review on these models. All these models provide interesting ideas but they are empirical and cannot be used in a general manner (for FE implementation for instance). In [11], the authors have proposed a specific strategy to incorporate an empirical dependency on the dynamic amplitude within a consistent standard viscoelastic model in the case of harmonic loadings. They explicitly defined a dynamic amplitude within a finite element simulation by the help of a static simulation. In [12, 13], a fully thermodynamically consistent approach was proposed to incorporate

<figure><img src="image_4.png" /><figcaption>Fig. 1 Schematic representation of the filler-filler mechanisms involved in the Payne effect</figcaption></figure>
a strain amplitude dependency. The idea proposed consist in introducing two
relaxation mechanisms in each Maxwell elements (called process-dependent vis-
cosity by the authors). This model was latter improved in [14] to account for the
reversibility and have the advantage of being transposable in the frequency space.
In [15, 16], we have proposed a model inspired by continuum damage mechan-
ics. As proposed by Krauss, we account for a progressive breakage of filler-filler
bonds within a generalized Maxwell model as the dynamic amplitude increases.
However, in this model, we do not account for reversibility neither reorganiza-
tion of filler-filler bonds under relaxation. Therefore this approach is not fully
consistent with the experiments and in particular do not capture correctly the
Payne effect in tension with preload. Nevertheless, in [17], we have also proposed
a strategy to use the model in the frequency domain. Other approaches have been
also proposed in [18–20]. In these approaches, the authors have proceeded differ-
ently, but the main idea is to introduce a nonlinear viscosity effect (sometimes in
the Maxwell damper other times in the spring element) but without adding addi-
tional internal variables or mechanisms from which the Payne effect could origi-
nate. Nevertheless, none of the previous mentioned approaches is able to accu-
rately capture all the effects and dependencies related to the dynamic softening
observed from an experimental point of view.

It is proposed here the development of a new constitutive model. To our point of view, compared to what exist in the literature, the main originality comes from two aspects. First, we account for the Payne effect and its partial reversibility by introducing appropriate internal variables with a clear aggregation/disaggregation scenario of filler-filler bonds, linked to observations from the literature. Second, we also account for the influence of the pre-strain or the preload which is commonly encountered in many industrial applications of rubber. This is done in a thermodynamically consistent way, and the proposed model can be implemented in a finite element code in a straightforward way. The evolution laws obtained have a structure similar to that of standard viscoelastic models, and their numerical integration poses no difficulty. The modeling is based on a decomposition of the non-equilibrium response into two contributions, one attached to the filler-filler interactions and the other one to the filler-matrix interaction. A phase-changing-like model allows to switch between the two in a reversible manner (such as schematized in Fig. 1). A generalized Maxwell model is used as the basis for the modeling to account for the frequency dependence, and the progressive breakage and re-arrangement of filler-filler bonds is integrated by a superposition of Maxwell elements that incorporate the previously mentioned phase-changing mechanism. An appropriate strategy for the identification of the material parameters is proposed.

The manuscript is organized as follows: A first section is devoted to the
description of the material considered in this study (natural rubber filled with car-
bon blacks) and the results of experimental tensile tests. The proposed modeling
approach is presented in a second section followed by the description of a specific
identification strategy in a third section. In the last section, we discuss about the
reliability of the proposed model compared to observed behavior.
Experimental characterization

For this study, we have selected a natural rubber filled with carbon blacks because this material is widely used by industry in many applications concerned with harmonic loadings and damping by viscoelasticity while maintaining good resistance to crack propagation. Tires are a typical example for such applications. The exact material formulation is given in Table 1. The material processing and the production of tensile samples have been done at the Leibniz-Institut für Polymerforschung e.V. Dresden, Germany.

Experimental tensile tests were performed to determine and characterize the mechanical behavior. Mechanical tests include quasi-static test and DMA (dynamic mechanical analysis). All the tests were performed at room temperature. In the following, for the sake of simplicity, we use the term “stress” in place of first Piola–Kirchhoff stress and “strain” in place of stretch minus one ($\lambda - 1$).

Quasi-static tests

The tests were carried out on normalized specimens (DIN 53504/S2 or ISO37 type 2), and the loading is defined in two sections: (i) Six loading–unloading cycles at 120% strain amplitude with a strain rate $\dot{\lambda} = 0.01 \text{ s}^{-1}$, corresponding to the *de-mul-<u>linization</u>* stage (Fig. 2), and (ii) multi-relaxation steps, where an increasing incremental strain of 10% with a strain rate $\dot{\lambda} = 0.125 \text{ s}^{-1}$ followed by a hold of the stretch during 300 s for each step. After attaining 10 steps, the strain amplitude decreased gradually in the same manner (see Fig. 3). From the stress versus time results (see Fig. 3), it can be seen that the viscosity is still present after 300 s. This relaxation time of 300 s results from a compromise between the test duration and the remaining part of the viscosity.

In Fig. 4, we have reported the apparent stiffness (current Piola-Kirchhoff stress divided by the current strain) on a logarithmic time scale. These curves illustrates that the stress relaxation is nonlinear (i.e., it depends on the strain level). This effect is widely observed in filled elastomers.

<table><caption>**Table 1** Rubber blend composition (value in phr, in weight)</caption><thead><tr><th></th><th>Component</th><th>phr</th></tr></thead><tbody><tr><td>Elastomer</td><td>NR CV60-3 L</td><td>100</td></tr><tr><td>Carbon black</td><td>N234</td><td>50</td></tr><tr><td>Antiozonant</td><td>IPPD</td><td>3</td></tr><tr><td>Antioxidant</td><td>TMQ</td><td>2</td></tr><tr><td>Wax</td><td>wax2124</td><td>2</td></tr><tr><td>Accelerator</td><td>TBBS</td><td>1.75</td></tr><tr><td>Vulcanizing agent</td><td>Sulfur</td><td>2.5</td></tr><tr><td>Activator</td><td>ZnO</td><td>5</td></tr><tr><td></td><td>Stearin</td><td>3</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Loading-unloading cycles performed on the rubber sample prior to the multi-stress relaxation steps</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Multi-stress relaxation steps involving both incremental and decremental strain tension</figcaption></figure>

A very standard result consist in plotting the so-called equilibrium hysteresis. This curve is obtained by keeping only the experimental points obtained at the end of each relaxation steps. If the material behavior is perfectly viscoelastic and not sensitive to any thermomechanical couplings, the response should be perfectly hyperelastic without any hysteresis effect. This result is reported at Fig. 5. It can be seen that an hysteresis effect is obtained highlighting that we have either long-time viscoelastic effects, non-neglectible thermomechanical couplings or other irreversible mechanisms (like viscoplasticity or damage for instance). This result is also widely observed in the literature for filled elastomers (even for very long relaxations), see for instance [21, 22].
<figure><img src="image_7.png" /><figcaption>Fig. 4 Apparent Young modulus (stress relaxation steps normalized by the corresponding strain step level) versus relaxation time in log scale</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5 Equilibrium hysteresis at room temperature</figcaption></figure>

## DMA experiments

Dynamical mechanical analysis experiments under tension were conducted, on rec-
tangular rubber strips, to investigate the frequency and strain amplitude dependen-
cies of the rubber blend. As we have used strip-like specimens, we also have to take
into account for a pre-stretch to avoid buckling of the specimen during the harmonic
load. To achieve this, the DMA tests were carried out at 1% and 10% static pre-strain
for the frequency sweep tests and 60% static pre-strain for the strain amplitude sweep tests.

Frequency sweep tests: In Fig. 6a, b, we report the results of the storage modulus ($E'$) and the loss factor ($\tan(\delta)$) with different applied pre-strain and different dynamic strain amplitudes. It can be seen that the material appears stiffer for the lowest applied pre-strain, and if we consider the same pre-strain, the behavior is softer for the larger dynamic amplitude. Furthermore, it can be seen that the loss factor is sensitive to the frequency and to the dynamic strain amplitude. The pre-strain effect on the loss factor is illustrated in Fig. 8b.

<figure><img src="image_9.png" /><figcaption>(a) Storage Modulus $E'$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Loss factor $\tan(\delta)$</figcaption></figure>

**Fig. 6** DMA frequency sweep tests for NR at different static ($\lambda_{st}$) and dynamic stretches ($\lambda_{dyn}$) at room temperature
In line with previous observations in the literature, we can see that the apparent modulus is sensitive to the pre-strain. (The curves in Fig. 6a are not superimposable.)

Amplitude sweep tests: Strain amplitude sweep tests were conducted at constant frequencies of 0.5, 5, and 10 Hz, with strain amplitude ranging from 0.1 to 30% (see Fig. 7a for the storage modulus and Fig. 7b for the loss factor). The

<figure><img src="image_11.png" /><figcaption>(a) Storage modulus $E'$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) tan $\delta$ results</figcaption></figure>

**Fig. 7** DMA strain amplitude sweep tests at frequencies 0.5, 5 and 10 Hz at room temperature and 60% pre-strain (strain amplitudes are in %)
Payne effect is observed (softening effect while increasing dynamic amplitudes), and it can be seen that this effect is dependent on the applied frequency.

The static pre-strain has also an important effect in the dynamic moduli values (here we have reported only the storage modulus $E'$), particularly evident in the strain amplitude tests results (Fig. 8a). This effect is also observable in the frequency sweep tests (Fig. 6a). In both DMA tests, smaller static pre-strains correlate with larger observed dynamic moduli. Figure 8b shows that increasing the pre-strain leads to lower loss factor.

<figure><img src="image_13.png" /><figcaption>Fig. 8 DMA strain amplitude sweep tests at 5 Hz and room temperature: pre-strain effect</figcaption></figure>
Modeling

To model the Payne effect, we propose to consider a scenario of breakage and re-
organization of the filler-filler network which is accompanied by a phase-change
mechanism of the bound rubber (which change from glassy to rubbery state).
Considering a generalized Maxwell approach with *n* viscous branches, we pro-
pose to introduce *n* scalar internal variables, *ω*<sup>*i*</sup>(**X**, *t*) ∈ [0, 1], that can be viewed
as the percentage of active filler-filler bounds associated with the *i*th viscous
mechanisms.

This leads to a progressive breakdown and re-organization of each Maxwell ele-
ments depending on the applied strain (or more precisely the strain energy) and on
some intrinsic characteristic times. The mechanical properties of the Maxwell ele-
ments are affected (viscosity and elastic stiffness) by this breakdown/re-organization
mechanisms which implies a dependence of these properties upon the variable ω<sup>*i*</sup>.

Figure 9 shows a schematic representation of the considered model. Each
Maxwell branch is put in parallel with a purely elastic branch. The branches are

<figure><img src="image_14.png" /><figcaption>**Fig. 9** Generalized Maxwell model with phase-changing elements for the Payne effect</figcaption></figure>

<figure><img src="image_15.png" /></figure>
connected with a progressive switch which represent the progressive breakdown and re-organization mechanisms.

## Thermodynamics

We use here the so-called local state approach, based on the definition of an appropriate set of state variables, free energy potentials and Onsager-like relations such as to fulfill the Clausius–Duhem inequality. The behavior is assumed to be isotropic and nearly incompressible, and viscosity is introduced with the help of intermediate states as proposed by [23, 24]. Furthermore, we neglect here all thermal effects and thermal dependencies (i.e., isothermal and adiabatic hypotheses). We briefly recall in this section the main steps to formulate such models.

Denoting by $\mathbf{X}$ the position of a material point in the reference frame and $\mathbf{x}(\mathbf{X}, t)$ its position in the current frame, it can be defined the deformation gradient as: $\mathbf{F}(\mathbf{X}, t) = \partial\mathbf{x}(\mathbf{X}, t)/\partial\mathbf{X}$. To take into account for the nearly incompressible behavior, we follow [25] and we adopt the so-called isochoric/volumic decomposition of the deformation gradient, and we split the isochoric part into $n$ elastic and $n$ viscous parts, such that (e.g., [26]):

$$ \mathbf{F} = \bar{\mathbf{F}}J^{1/3} = \bar{\mathbf{F}}_e^i \bar{\mathbf{F}}_v^i J^{1/3}, \quad i = 1 \dots n \quad (1) $$

where $J$ stands for the volume variation: $J = \det(\mathbf{F})$.

From Eq. (1), we can define the left isochoric Cauchy–Green tensors: $\bar{\mathbf{B}} = \bar{\mathbf{F}}\bar{\mathbf{F}}^T$ and $\bar{\mathbf{B}}_e^i = \bar{\mathbf{F}}_e^i \bar{\mathbf{F}}_e^i J^i$.

We assume that the Helmholtz free energy, $\psi$, can be additively decomposed into a purely hyperelastic part and viscous parts as done for a standard generalized Maxwell model:

$$ \psi(\bar{\mathbf{B}}, \bar{\mathbf{B}}_e^i, \omega^i, J) = \psi_{\text{hyp}}(\bar{\mathbf{B}}) + \psi_{\text{vol}}(J) + \sum_{i=1}^{n} \psi_{\text{vis}}^i(\bar{\mathbf{B}}, \bar{\mathbf{B}}_e^i, \omega^i) \quad (2) $$

where $\psi_{\text{hyp}}$ and $\psi_{\text{vol}}$ are the initial hyperelastic and volumic energies. The $i$th viscous free energy is then assumed to be composed of two contributions. The first one is related to the filler-filler network that is assumed to progressively transform into an hyperelastic contribution as the chains locally change of state from glassy to rubbery. The opposite phenomena takes place when the strain amplitude is decreased, the chain mobility is reduced, and the filler-filler network is strengthened. This leads to the following supplementary additive decomposition:

$$ \psi_{\text{vis}}^i(\bar{\mathbf{B}}, \bar{\mathbf{B}}_e^i, \omega^i) = \underbrace{\omega^i \psi_{j-f}^i(\bar{\mathbf{B}}_e^i)}_{\text{filler-filler}} + (1 - \omega^i) \underbrace{\psi_{c-f}^i(\bar{\mathbf{B}})}_{\text{chain-filler}} \quad (3) $$

Combining the first and second fundamental thermodynamics principles, the following form of the Clausius–Duhem (or Coleman–Noll) inequality can be obtained:
$$ \phi = \sigma : \mathbf{D} - \rho \dot{\psi} \geq 0 \qquad (4) $$

where $\phi$ is the intrinsic dissipation, $\sigma$ is the Cauchy stress tensor, $\mathbf{D}$ the Eulerian rate of deformation, and $\rho$ the actual density. The time variation of the Helmholtz free energy is given by:

$$ \dot{\psi} = \left( \frac{\partial \psi_{\text{hyp}}}{\partial \bar{\mathbf{B}}} + \sum_{i=1}^{n} (1 - \omega^i) \frac{\partial \psi_{c-f}^i}{\partial \bar{\mathbf{B}}} \right) : \dot{\bar{\mathbf{B}}} + \frac{\partial \psi_{\text{vol}}}{\partial J} \dot{J} + \sum_{i=1}^{n} \omega^i \frac{\partial \psi_{f_f}}{\partial \bar{\mathbf{B}}_e^i} : \dot{\bar{\mathbf{B}}}_e^i \\ + \sum_{i=1}^{n} (\psi_{f-f}^i - \psi_{c-f}^i) \dot{\omega}^i \qquad (5) $$

and the time variation of $\dot{\bar{\mathbf{B}}}, J$ and $\bar{\mathbf{B}}_e^i$ can be computed as follows:

$$ \dot{\bar{\mathbf{B}}} = \mathbf{L}\bar{\mathbf{B}} + \bar{\mathbf{B}}\mathbf{L}^T - \frac{2}{3}(\mathbf{1} : \mathbf{D})\bar{\mathbf{B}} \qquad (6) $$

$$ \dot{J} = J(\mathbf{1} : \mathbf{D}) \qquad (7) $$

$$ \dot{\bar{\mathbf{B}}}_e^i = \mathbf{L}\bar{\mathbf{B}}_e^i + \bar{\mathbf{B}}_e^i\mathbf{L}^T - \frac{2}{3}(\mathbf{1} : \mathbf{D})\bar{\mathbf{B}}_e^i - 2\bar{\mathbf{V}}_v^i\bar{\mathbf{D}}_v^{\prime o}\bar{\mathbf{V}}_e^i \qquad (8) $$

where $\mathbf{1}$ is the second-order unit tensor, $\mathbf{L} = \dot{\mathbf{F}}\mathbf{F}^{-1}$ is the Eulerian velocity gradient, $\bar{\mathbf{V}}_e^i$ comes from the polar decomposition, $\bar{\mathbf{F}}_e^i = \bar{\mathbf{V}}_e^i\bar{\mathbf{R}}_e^i$, and $\bar{\mathbf{D}}_v^{\prime o}$ is an objective (frame independent) inelastic rate defined from: $\bar{\mathbf{D}}_v^{\prime o} = \bar{\mathbf{R}}_e^i\bar{\mathbf{D}}_v^{\prime o}\bar{\mathbf{R}}_e^{\prime T}$.

Introducing Eqs. (6)–(8) in Eq. (4) and using symmetry properties, one can rewrite the dissipation as:

$$ \phi = \left( \sigma - 2\rho \left( \bar{\mathbf{B}} \frac{\partial \psi_{\text{hyp}}}{\partial \bar{\mathbf{B}}} + \sum_{i=1}^{n} (1 - \omega^i) \bar{\mathbf{B}} \frac{\partial \psi_{c-f}^i}{\partial \bar{\mathbf{B}}} + \sum_{i=1}^{n} (\omega^i) \bar{\mathbf{B}}_e^i \frac{\partial \psi_{f-f}^i}{\partial \bar{\mathbf{B}}_e^i} \right)^{\mathrm{D}} - \rho J \frac{\partial \psi_{\text{vol}}}{\partial J} \mathbf{1} \right) : \mathbf{D} \\ + 2\rho \sum_{i=1}^{n} \omega^i \left( \bar{\mathbf{B}}_e^i \frac{\partial \psi_{f-f}^i}{\partial \bar{\mathbf{B}}_e^i} \right) : \bar{\mathbf{D}}_v^{\prime o} + \sum_{i=1}^{n} -\rho (\psi_{f-f}^i - \psi_{c-f}^i) \dot{\omega}^i \geq 0 \qquad (9) $$

the exponent, $^{D}$, stands for the deviatoric operator.

Assuming that the only sources of dissipation are related to $\bar{\mathbf{B}}_e^i$ and $\omega^i$, it can be obtained from the previous equation a constitutive expression for the stress tensor (e.g., [27] for a similar argument):

$$ \sigma = 2\rho \left( \bar{\mathbf{B}} \frac{\partial \psi_{\text{hyp}}}{\partial \bar{\mathbf{B}}} + \sum_{i=1}^{n} (1 - \omega^i) \bar{\mathbf{B}} \frac{\partial \psi_{c-f}^i}{\partial \bar{\mathbf{B}}} + \sum_{i=1}^{n} \omega^i \bar{\mathbf{B}}_e^i \frac{\partial \psi_{f-f}^i}{\partial \bar{\mathbf{B}}_e^i} \right)^{\mathrm{D}} - \rho J \frac{\partial \psi_{\text{vol}}}{\partial J} \mathbf{1} \qquad (10) $$
Evolution laws

For the viscoelastic contribution, we consider an Onsager-like relation (see [28]) such that the objective inelastic deformation rate is directly proportional to its associated thermodynamics force, so that:

$$ \bar{\mathbf{D}}_v^{\omega^i} = \frac{\omega^i}{\eta_i(\omega^i)} \left( 2\rho \bar{\mathbf{B}}_e^i \frac{\partial \psi_{f-f}^i}{\partial \bar{\mathbf{B}}_e^i} \right)^D \quad (11) $$

where $\eta_i$ is a viscosity coefficient that can depend on the percentage of active filler-filler network. As the filler-filler bonds become lower, i.e., as $\omega^i$ is decreasing, we consider that the viscosity should also decrease such as to become null if all the bonds are broken. Inserting Eq. (11) into Eq. (8), the viscoelastic flow can be rewritten as:

$$ \dot{\bar{\mathbf{B}}}_e^i = \mathbf{L}\bar{\mathbf{B}}_e^i + \bar{\mathbf{B}}_e^i\mathbf{L}^T - \frac{2}{3}(\mathbf{1}:\mathbf{D})\bar{\mathbf{B}}_e^i - \rho \frac{4\omega^i}{\eta_i(\omega^i)} \bar{\mathbf{V}}_e^i \left( \bar{\mathbf{B}}_e^i \frac{\partial \psi_{f-f}^i}{\partial \bar{\mathbf{B}}_e^i} \right)^D \bar{\mathbf{V}}_e^i \quad (12) $$

For the percentage of active filler-filler bonds, $\omega^i$, it can be seen from Eq. (9) that an admissible evolution law could be such that the sign of $\dot{\omega}^i$ is given by the sign of $\psi_{c-f}^i - \psi_{f-f}^i$. Therefore, the softening effect, due to the phase changing of glassy bonds to rubbery one or to the breakage of this bonds, will take place if the strain energy of the filler-filler network is higher than the one of the chain-filler and the opposite for the stiffening or recovery effect. To account both of the time dependency for the softening or recovery (see for instance [14]) and the effect of the maximum applied amplitude, we propose to use a model inspired by viscoplasticity as done in our previous approaches (see [15, 16]). Therefore, we propose the following law:

$$ \dot{\omega}^i = \begin{cases} -\frac{1}{h_i} < \omega^i - 1 + \tanh\left(\frac{\rho_0(\psi_{f-f}^i - \psi_{c-f}^i)}{b_i}\right) > & \text{if } \psi_{f-f}^i - \psi_{c-f}^i \ge 0 \\ \frac{S}{h_i} < 1 - \omega^i + \tanh\left(\frac{\rho_0(\psi_{f-f}^i - \psi_{c-f}^i)}{b_i}\right) > & \text{else} \end{cases} \quad (13) $$

where $< \cdot >$ denotes the positive part of the argument (Macauley brackets), $h_i$ is a characteristic time that corresponds to the stabilizing time for the Payne effect, $b_i$ is an energetic parameter that control the rate of damage/healing depending on the difference of strain energy in the filler-filler and chain-filler contribution, and $S$ is a dimensionless factor that control the difference in term of characteristic time for the breakage or recovery phenomena. We will assume in the following that longer time is needed for the recovery compared to the damage phase (i.e., $S \le 1$). Experimentally, it is nearly impossible to separate the stabilization time due to viscosity and the one due to Payne (if there exist one) but one can think about this as a sort of
homogenization procedure: the characteristic time, $h_i$, could be interpreted as a specific time for the phase-change phenomena inside the filler network.

For the sake of simplicity, we assume that the factor $S$ is the same for all viscous mechanisms. As $h_i$ and $S$ are positive parameters, it is straightforward to see that the Clausius–Duhem inequality is always satisfied. It can also be remarked that contrarily to the model we have previously developed in [17], evolution laws defined by Eqs. (12) and (13) are fully coupled and form a nonlinear differential system to be solved for each Maxwell element.

## Free energies

As already mentioned, we assume an isotropic behavior. It is chosen to adopt an invariant representation of the Helmholtz free energies previously defined. For the hyperelastic contribution, without loss of generality, we have selected a polynomial form (that could be attributed to [29]):

$$ \rho_0 \psi_{\text{hyp}} = c_{10}(I_1(\bar{\mathbf{B}}) - 3) + c_{20}(I_1(\bar{\mathbf{B}}) - 3)^2 + c_{01}(I_2(\bar{\mathbf{B}}) - 3) \quad (14) $$

$$ \rho_0 \psi_{\text{vol}} = \frac{k}{2}(\log J)^2 \quad (15) $$

where $c_{10}$, $c_{20}$ and $c_{01}$ control the deviatoric part of the stress and $k$ is the bulk modulus.

For the filler–filler and chain–filler energies, for the sake of simplicity, it adopted neo-Hookean forms:

$$ \rho_0 \psi_{f-f}^i = \mu_i (I_1(\bar{\mathbf{B}}_e^i) - 3) \quad (16) $$

$$ \rho_0 \psi_{c-f}^i = A \mu_i (I_1(\bar{\mathbf{B}}) - 3) \quad (17) $$

where $0 \le A \le 1$ is a infinitesimal shear modulus ratio of the chain–filler and filler–filler contribution.

Using the free energy of Eq. (16) in Eq. (12) and choosing a viscosity function such that $\eta_i(\omega^i) = \omega^i \tilde{\eta}_i$, one can obtain the following Maxwell-like equation for the viscous left Cauchy–Green tensor

$$ \dot{\bar{\mathbf{B}}}_e^i = L \bar{\mathbf{B}}_e^i + \bar{\mathbf{B}}_e^i L^T - \frac{2}{3} (1 : D) \bar{\mathbf{B}}_e^i - \frac{1}{\omega^i \tau_i J} \bar{\mathbf{B}}_e^i (\bar{\mathbf{B}}_e^i)^D \quad (18) $$

where $\tau_i = \tilde{\eta}_i/4\mu_i$ are the so-called characteristic time of viscosity. The evolution equation for the variable $\omega^i$ takes the following form with the previous choices for the free energies
$$ \dot{\omega}^i = \begin{cases} -\frac{1}{h_i} < \omega^i - 1 + \tanh \left( \frac{\mu_i}{b_i} (I_1(\bar{\mathbf{B}}_e^i) - 3 - A(I_1(\bar{\mathbf{B}}) - 3)) \right) > & \text{if } I_1(\bar{\mathbf{B}}_e^i) - 3 - A(I_1(\bar{\mathbf{B}}) - 3) \ge 0 \\ \frac{S}{h_i} < 1 - \omega^i + \tanh \left( \frac{\mu_i}{b_i} (I_1(\bar{\mathbf{B}}_e^i) - 3 - A(I_1(\bar{\mathbf{B}}) - 3)) \right) > & \text{else} \end{cases} \quad (19) $$

To simplify the definition of material parameters, we can adopt the following notation: $\tilde{b}_i = \mu_i/b_i$ such as $\tilde{b}_i$ is determined by the identification procedure rather than $b_i$. The previous equation can be rewritten as:

$$ \dot{\omega}^i = \begin{cases} -\frac{1}{h_i} < \omega^i - 1 + \tanh \left( \frac{I_1(\bar{\mathbf{B}}_e^i) - 3 - A(I_1(\bar{\mathbf{B}}) - 3)}{\tilde{b}_i} \right) > & \text{if } I_1(\bar{\mathbf{B}}_e^i) - 3 - A(I_1(\bar{\mathbf{B}}) - 3) \ge 0 \\ \frac{S}{h_i} < 1 - \omega^i + \tanh \left( \frac{I_1(\bar{\mathbf{B}}_e^i) - 3 - A(I_1(\bar{\mathbf{B}}) - 3)}{\tilde{b}_i} \right) > & \text{else} \end{cases} \quad (20) $$

It can be remarked from the previous equation that if $A = 0$, then only the damage can take place, and Eq. (20) becomes:

$$ \dot{\omega}^i = \frac{1}{h_i} < \omega^i - 1 + \tanh \left( \frac{I_1(\bar{\mathbf{B}}_e^i) - 3}{\tilde{b}_i} \right) > \quad (21) $$

To illustrate the influence of the material parameter $\tilde{b}_i$, we can consider the case of harmonic loadings in simple shear. In this case, the deformation gradient can be written as $\mathbf{F} = \mathbf{1} + \gamma(t)\mathbf{e}_1 \otimes \mathbf{e}_2$, where $\mathbf{e}_1, \mathbf{e}_2$ are the two first elements of the considered Cartesian basis (in the initial configuration) and $\gamma(t) = \gamma_0 \sin(2\Pi ft)$. We also restrict the model to one Maxwell branch with an infinite relaxation time such that $\bar{\mathbf{B}}_e = \bar{\mathbf{B}}$. The first invariant is in this case $I_1(\bar{\mathbf{B}}_e) = I_1(\bar{\mathbf{B}}) = 3 + \gamma(t)^2$. Neglecting time effect and considering only the equilibrium value of $\omega$ (i.e., value reached at stabilization, noted $\omega_{eq}$), an analytical solution can be obtained from Eq. (20), such that (we consider the case $A = 0$):

$$ \omega_{eq} = 1 - \tanh \left( \frac{\gamma 0^2}{\tilde{b}} \right) \quad (22) $$

Figure 10 illustrates the role of $\tilde{b}$ using the previous expression, and this parameter allows to control the strain level at which the softening will take place.

In the case $A = 1$ and if the relaxation time is infinite (or very large compared to the loading rate), there are no damage or healing and $\dot{\omega} = 0$.
<figure><img src="image_16.png" /><figcaption>Fig. 10 Influence of b for a cyclic shear test with only damage and infinite relaxation time. The curves correspond to the stabilized value of $\omega$ for different shear amplitudes</figcaption></figure>

# Identification of the material parameters

## An identification strategy

We must develop an optimal and practical strategy to identify the model parameters for each branch, as the increase in Maxwell branches also increases the number of parameters, potentially making the calculations lengthy. For each branch, we fixed a combination of $\tau_i$, $h_i$, and $b_i$. The values for S and A were selected in a preliminary stage of the identification process, in which we tested different combinations of these parameters for a simplified version of the model with a single Maxwell element (and considering only a DMA amplitude sweep test).

$$ \tau_i = \{10^{-3}, 10^{-2}, 10^{-1}, 10^{0}, 10^{1}, 10^{2}, 10^{3}, 10^{4}\}(s) \quad (23) $$

$$ h_i = \{10^{-2}, 10^{-1}\}(s) \quad (24) $$

$$ b_i = \{10^{-3}, 10^{-2}, 10^{-1}, 10^{0}, 10^{1}, 10^{2}, 10^{3}\}(MPa) \quad (25) $$

$$ S = 10^{-3} \quad (26) $$

$$ A = 10^{-2} \quad (27) $$

Then, the identification problems is defined such that to minimize a least square distance between the experimental data and the one obtained by the model for the viscous modulus $\mu_i$ of each branch and the hyperelastic coefficients $c_{ij}$. With the
previous chosen values, there are 80 possible Maxwell branches but it has to be remarked that the corresponding minimization problem is linear.

The least square error is calculated using multi-step relaxation and dynamic mechanical analysis (DMA) experimental data.

$$f_{\text{Total-err}} = f_{\text{MSR-err}} + f_{\text{DMA-err}} \qquad (28)$$

where $f_{\text{Total-err}}$ is the objective function, $f_{\text{MSR-err}}$ is the difference between the total stress from the multi-step relaxation tests and the total stress results from the theoretical model at all time points ($N_{\text{temp,time}}$), and $f_{\text{DMA-err}}$ is the difference between experimental and theoretical real modulus $E'$ and hysteresis ($A_{\text{hyst}}$) from DMA tests in both frequency and strain amplitude sweep.

$$f_{\text{MSR-err}} = \sum_{i}^{N_{\text{temp,time}}} \frac{(\Pi_{\text{exper}}^{i} - \Pi_{\text{theor}}^{i})^{2}}{(\Pi_{\text{exper}}^{i})^{2}} \qquad (29)$$

$$f_{\text{DMA-err}} = \sum_{i}^{N_{\text{temp,time}}} \left( \frac{(E_{\text{exper}}^{i} - E_{\text{theor}}^{i})^{2}}{(E_{\text{exper}}^{i})^{2}} + \frac{(A_{\text{hyst,exper}}^{i} - A_{\text{hyst,theor}}^{i})^{2}}{(A_{\text{hyst,exper}}^{i})^{2}} \right) \qquad (30)$$

To solve this problem, we have used Mathematica. It allows for symbolic and numeric computations (which is useful to express the model response for quasi-static and DMA test). The minimization of the objective function given at Eq. (28) is realized under constraints. This constraints are essential to obtain an acceptable behavior (in the thermodynamics and physical sense), and we simply impose [$c_{ij} \ge 0, \mu_i \ge 0$]. At the end, we can eliminate inactive branches (the one that do not contribute significantly to the mechanical response, for instance if $\mu_i$ is very small). This approach allows us to determine the optimal parameters corresponding to the active branches while limiting the time required for the identification due to the fact that we limit ourselves to a linear minimization problem under constraints.

## Identification results

The identified values of the material parameters are given in Tables 2 and 3. It can be seen that the final number of Maxwell branches has been reduced to 15 thanks to our identification strategy.

Figures 11 and 12 show the results obtained from the multi-step experiment. In the identification process, we have used only one half of the experimental tests to compute the objective function (corresponding to only increasing step amplitude). The idea was to find the minimal set of data needed to obtain a correct agreement between model and experiments (especially for these step relaxation tests that are time expansive).

For the DMA tests, it can be seen from Figs. 13 and 14 that the model can capture with a relative good agreement the frequency sweep behavior for the
<table><caption>**Table 2** Identification of viscoelastic parameters</caption><thead><tr><th>μ<sub>i</sub> (MPa)</th><th>τ<sub>i</sub> (s)</th><th>b<sub>i</sub> (MPa)</th><th>h<sub>i</sub> (s)</th></tr></thead><tbody><tr><td>0.105</td><td>10</td><td>0.001</td><td>0.01</td></tr><tr><td>1.020</td><td>1000</td><td>0.001</td><td>0.01</td></tr><tr><td>0.082</td><td>0.1</td><td>0.01</td><td>0.01</td></tr><tr><td>0.840</td><td>1</td><td>0.01</td><td>0.01</td></tr><tr><td>0.020</td><td>1000</td><td>0.1</td><td>0.01</td></tr><tr><td>0.179</td><td>10,000</td><td>0.1</td><td>0.01</td></tr><tr><td>0.306</td><td>0.01</td><td>1</td><td>0.01</td></tr><tr><td>0.006</td><td>100</td><td>10</td><td>0.01</td></tr><tr><td>0.016</td><td>1000</td><td>10</td><td>0.01</td></tr><tr><td>0.601</td><td>0.001</td><td>100</td><td>0.01</td></tr><tr><td>0.540</td><td>0.01</td><td>100</td><td>0.01</td></tr><tr><td>0.022</td><td>0.1</td><td>100</td><td>0.01</td></tr><tr><td>0.011</td><td>10</td><td>1000</td><td>0.01</td></tr><tr><td>0.059</td><td>1000</td><td>1000</td><td>0.01</td></tr><tr><td>0.111</td><td>10,000</td><td>1000</td><td>0.01</td></tr></tbody></table>

<table><caption>**Table 3** Identification of hyperelastic coefficients (in MPa)</caption><tr><td>c<sub>10</sub></td><td>0.022</td></tr><tr><td>c<sub>01</sub></td><td>0.572</td></tr><tr><td>c<sub>20</sub></td><td>0.052</td></tr></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 11** Experimental data fit results for multi-step stress relaxation at room temperature</figcaption></figure>

storage modulus. However, the loss factor is not correctly described for higher frequencies (c.f. Fig. 14). Considering the loss factor, this result may have different reasons, the two main hypothesis being model definition or identification process and should be investigated in future studies. For the identification
<figure><img src="image_18.png" /><figcaption>Fig. 12 Experimental data fit results for multi-step stress relaxation (time) at room temperature</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 13 Experimental data fit results for storage modulus in frequency sweep at room temperature and $\lambda_{st} = 10\%$, $\lambda_{dyn} = 0.2\%$</figcaption></figure>

process, we feel that weights could be added to Eqs. (28)–(30) to improve the results for the loss factor, but a compromise must be found so as not to deteriorate the results too much for the storage modulus and multi-step stress relaxation.

Figures 15 and 16 show the results for the amplitude sweep. Similar results can be observed: The storage modulus is accurately fitted but the loss factor is not perfectly described.

Overall, as expected, the model effectively replicates the nonlinear behavior observed in the experimental tests and used for the identification process. The discrepancies observed when focusing on the loss factor may be due to an inappropriate
<figure><img src="image_20.png" /><figcaption>**Fig. 14** Experimental data fit results for tan δ in frequency sweep at room temperature and λst = 10%, λdyn = 0.2%</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 15** Experimental data fit results for storage modulus in strain amplitude sweep, with a static strain of 60% and at a frequency of 0.5 Hz</figcaption></figure>

or incomplete set of initial parameters. In particular, the choice of decades for the characteristic times of viscosity may have a strong impact to fully describe the relaxation spectrum.
<figure><img src="image_22.png" /><figcaption>**Fig. 16** Experimental data fit results for tan δ in strain amplitude sweep, with a static strain of 60% and at a frequency of 0.5 Hz</figcaption></figure>

# Discussions about the model performance

## Payne effect recovery

As mentioned in the introduction and in the modeling part, one important new aspect (compared to our previous works) of the proposed model is to account for Payne reversibility.

To assess the reversibility of the Payne effect, we can set a case example, from the modeling point of view, where a dynamic load ($\lambda(t)$) with variable amplitude is applied (Fig. 17) at a low frequency (0.5 Hz) with a static strain of 60%, incrementally increasing and then subsequently decreasing the dynamic amplitude. In this

<figure><img src="image_23.png" /><figcaption>**Fig. 17** Amplitude variable dynamic load: $\lambda_{dyn1} = 0.05$, $\lambda_{dyn2} = 0.1$, $\lambda_{dyn3} = 0.2$, $t_{pre-stretch} = 500$ s</figcaption></figure>
manner, it is possible to assess the stiffness recovery, indicating that the Payne effect
is reversible to a certain extent (Fig. 18).

Our findings indicate that the Payne effect exhibits a degree of reversibility, though this reversibility is influenced by several factors, including frequency. Specifically, as the frequency increases, the recovery percentage of E' decreases, suggesting that the material's ability to recover stiffness diminishes under higher-frequency conditions. This highlights the complex interplay between dynamic loading conditions and the inherent properties of the material. Here, we do not have a direct comparison between model and experiments; however, this effect has been observed from an experimental point view for other materials in the literature, e.g., [12, 16].

Multi-step stress relaxation in loading and unloading

As detailed in the experimental tests section, we conducted multi-step stress relaxa-
tion during both the loading and unloading phases. These tests aimed to systemati-
cally evaluate the material’s stress response under different conditions. Our results
confirm that the model accurately simulates the stress response not only during the
loading steps, but also during the unloading phase. This demonstrates the model’s
robustness and reliability in capturing the material’s complex behavior under both
increasing and decreasing stress conditions (see Figs. 19 and 20). Nevertheless, it
has to be noted that the model underestimates the viscosity at higher stretches with-
out any clear explanation for this.

Pre-strain effect

One of the primary objectives of the proposed model is to accurately capture the
effect of static strain applied prior to dynamic loading. This is crucial for under-
standing how preexisting deformation influences material behavior under sub-
sequent dynamic conditions. As illustrated in Fig. 8a, b, the magnitude of the

<figure><img src="image_24.png" /><figcaption>Fig. 18 E' recovery during the amplitude variable dynamic load (obtained from the current model)</figcaption></figure>

<figure><img src="image_25.png" /></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 19** Multi-step stress relaxation in loading-unloading model simulation (line) and experimental data (dashed) comparison, stress vs time at room temperature</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 20** Multi-step stress relaxation in loading-unloading model simulation (line) and experimental data (dashed) comparison, stress vs stretch at room temperature</figcaption></figure>

pre-strain significantly impacts the dynamic response, with higher pre-strains leading to a noticeable decrease in dynamic moduli. The model effectively replicates this observed trend, demonstrating its capability to account for the interplay between static pre-strain and dynamic load responses (Figs. 21 and 22).
<figure><img src="image_28.png" /><figcaption>**Fig. 21** Real modulus *E'* simulation results at 10% (red), 20% (green) and 60% (blue) pre-strain and 0.5 Hz (color figure online)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 22** Tan δ simulation results at 10% (red), 20% (green) and 60% (blue) pre-strain and 0.5 Hz (color figure online)</figcaption></figure>

**Conclusions**

In this study, we propose a new thermodynamically consistent mechanical model for filler-reinforced vulcanized rubbers. By incorporating phase-changing mechanisms, the model provides deeper insights into the cyclic breakdown and re-organization of filler-filler bonds, offering a better understanding of the underlying nonlinear viscoelastic behavior, particularly the Payne effect, associated with the filler network and matrix interactions.

Through an extensive experimental campaign, where CB-reinforced NR was tested under various loading conditions, we have obtained comprehensive experimental observations. These observations confirm the results already known from the literature and have guided us in the development a robust model that accurately
replicates the observed behavior of NR filled with carbon blacks. We have demonstrated the model's capability to replicate the impact of load histories and the partial reversibility of the Payne effect.

While focused on carbon-filled rubbers, the framework of the model can potentially be adapted for other types of filled elastomers (obviously strain energies and identification process may have to be adapted) thus expanding its utility in the field. Furthermore, the model's thermodynamically consistent formulation allows for straightforward implementation in finite element simulations, facilitating its use in a broad range of applications in engineering design and analysis.

Overall, the proposed model represents a significant advancement in the field of rubber mechanics. This model and the identification procedure still need to be improved to take proper account for the evolution of the loss factor for amplitude sweep tests. Future work could involve further refinement of the model (in particular to account for a thermomechanical coupling and temperature dependence of the behavior) and its application to a broader range of loading conditions and material compositions.

**Acknowledgements** This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant Agreement No. 958174. This project is co-financed by tax revenues on the basis of the budget adopted by the Saxon State Parliament and by the French National Research Agency (ANR) under the Grant Agreement ANR-22-MER3-0007-01.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Heinrich G, Klüppel M, Vilgis TA (2002) Reinforcement of elastomers. Curr Opin Solid State Mater Sci 6(3):195-203. https://doi.org/10.1016/S1359-0286(02)00030-X

2. Klüppel M (2003) In: Capella B, Geuss M, Klüppel M, Munz M, Schulz E, Sturm H (eds) The role of disorder in filler reinforcement of elastomers on various length scales. Springer, Berlin, pp 1-86. https://doi.org/10.1007/b11054

3. Mullins L (1948) Effect of stretching on the properties of rubber. Rubber Chem Technol 21:281-300. https://doi.org/10.5254/1.3546914

4. Payne AR (1960) A note on the existence of a yield point in the dynamic modulus of loaded vulcanizates. J Appl Polym Sci 3(7):127-127. https://doi.org/10.1002/app.1960.070030721

5. Klüppel M, Schuster RH, Heinrich G (1997) Structure and properties of reinforcing fractal filler networks in elastomers. Rubber Chem Technol 70(2):243-255

6. Fröhlich J, Niedermeier W, Luginsland H-D (2005) The effect of filler-filler and filler-elastomer interaction on rubber reinforcement. Compos Part A: Appl Sci Manuf 36(4):449-460. https://doi.org/10.1016/j.compositesa.2004.10.004

7. Warasithinon N, Genix A-C, Sztucki M, Oberdisse J, Robertson CG (2019) The payne effect: primarily polymer-related or filler-related phenomenon? Rubber Chem Technol 92(4):599-611
8. Chazeau L, Brown JD, Yanyo LC, Sternstein SS (2000) Modulus recovery kinetics and other insights into the payne effect for filled elastomers. Polym Compos 21(2):202–222

9. Kraus G, Childers CW, Rollmann KW (1966) Stress softening in carbon black reinforced vulcanizates: strain rate and temperature effects. Rubber Chem Technol 39(5):1530–1543. https://doi.org/10.5254/1.3547068

10. Heinrich G, Klüppel M (2002) Recent advances in the theory of filler networking in elastomers. Springer, Berlin, pp 1–44. https://doi.org/10.1007/3-540-45362-8_1

11. Gil-Negrete N, Viñolas J, Kari L (2006) A simplified methodology to predict the dynamic stiffness of carbon-black filled rubber isolators using a finite element code. J Sound Vib 296(4):757–776. https://doi.org/10.1016/j.jsv.2006.03.038

12. Rendek M, Lion A (2010) Amplitude dependence of filler-reinforced rubber: experiments, constitutive modelling and FEM—implementation. Int J Solids Struct 47(21):2918–2936. https://doi.org/10.1016/j.ijsolstr.2010.06.021

13. Rendek M, Lion A (2010) Strain induced transient effects of filler reinforced elastomers with respect to the payne-effect: experiments and constitutive modelling. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 90(5):436–458

14. Höfer P, Lion A (2009) Modelling of frequency- and amplitude-dependent material properties of filler-reinforced rubber. J Mech Phys Solids 57(3):500–520. https://doi.org/10.1016/j.jmps.2008.11.004

15. Delattre A, Lejeunes S, Méo S, Lacroix F, Richard C (2014) On the multiaxial amplitude and frequency dependent behavior of rubber: experiments and constitutive modeling. Rubber Chem Technol 87(3):557–578. https://doi.org/10.5254/rct.14.85995

16. Delattre A, Lejeunes S, Lacroix F, Méo S (2016) On the dynamical behavior of filled rubbers at different temperatures: experimental characterization and constitutive modeling. Int J Solids Struct 90:178–193. https://doi.org/10.1016/j.ijsolstr.2016.03.010

17. El Archi Y, Lahellec N, Lejeunes S, Jouan A (2024) Modeling the nonlinear dynamic behavior of filled elastomers: application to the vibration of structures with large preloads. Eur J Mech A/Solids. https://doi.org/10.1016/j.euromechsol.2024.105367

18. Narayan SPA, Palade LI (2020) Modeling payne effect with a framework of multiple natural configurations. Int J Eng Sci 157:103396. https://doi.org/10.1016/j.ijengsci.2020.103396

19. Kim S, Park S, Lee W, Kaliske M (2023) Numerical investigation of radial passenger car tire standing waves. Finite Elem Anal Des 220:103945. https://doi.org/10.1016/j.finel.2023.103945

20. Jalocha D (2020) Payne effect: a constitutive model based on a dynamic strain amplitude dependent spectrum of relaxation time. Mech Mater 148:103526. https://doi.org/10.1016/j.mechmat.2020.103526

21. Lion A (1997) A physically based method to represent the thermo-mechanical behaviour of elastomers. Acta Mecahnica 123:1–25. https://doi.org/10.1007/BF01178397

22. Haupt P, Sedlan K (2001) Viscoplasticity of elastomeric materials: experimental facts and constitutive modelling. Arch Appl Mech 71(2):89–109. https://doi.org/10.1007/s004190000102

23. Sidoroff F (1974) Un modéle viscoélastique non linéaire avec configuration intermédiaire. J Méc 13(4):679–713

24. Sidoroff F (1975) Variables internes en viscoélasticité, 2. milieux avec configuration intermédiaire. J Méc 14(4):571–595

25. Flory RJ (1961) Thermodynamic relations for highly elastic materials. Trans Faraday Soc 57:829–838. https://doi.org/10.1039/TF9615700829

26. Sidoroff F (1973) The geometrical concept of intermediate configuration and elastic–plastic finite strain. Arch Mech 25(2):299–308

27. Coleman BD, Noll W (1963) The thermodynamics of elastic materials with heat conduction and viscosity. Arch Ration Mech Anal 13:167–178. https://doi.org/10.1007/BF01262690

28. Onsager L (1931) Reciprocal relations in irreversible processes I. Phys Rev 37:405–426. https://doi.org/10.1103/PhysRev.37.405

29. Biderman VL (1958) Calculation of rubber parts. Raschetti Na Prochnost, Moscow ((in Russian))