Modeling Payne effect on basis of linearization of a visco-hyperelastic model

Rubber components are widely used as elastic connections in engineering structures to achieve a combination of stiffness and vibration isolation [1]. The best example is the engine mounts in the automotive and aircraft industries. Also, the tire performance is related to the thermal and viscoelastic properties of rubber of tires. Indeed, Ebbot et al [2] have proposed a coupled thermo-mechanical model, in which the stiffness and the loss are functions of the strain, temperature, and loading frequency.

Accordingly, the design of these rubber parts requires considering various complex effects: frequency, pre-load, thixotropy, etc. Much research has been devoted to study all these dependencies, particularly for carbon black-filled rubber (CBR) which usually present complex behaviors [3, 4]. The effects of Mullins and Payne are intrinsic to the behavior of reinforced rubbers (by carbon black or silica). The stress-softening phenomenon or the Mullins effect is characterized by a loss of stiffness during the first cycle of loading [5, 6]. This loss of rigidity will then be common to the following cycles (loading/unloading), that is provided that the applied strain does not go beyond the maximum strain reached in the history of deformation. Several authors have attributed the Mullins effect to bond rupture [7], molecular disentanglement [8], and multiple networks [9, 10]. Phenomenological descriptions of the Mullins effect was proposed in the framework of the continuum damage mechanics [11–14]. The reader may refer to a recent review article by Diani et al [15]. We note that the Mullins effect remains a hurdle in the formulation of hyperelastic models. Indeed, one should remove it experimentally in order to reach a stable material behavior, i.e. the hyperelastic behavior. Consequently, the rubbers are subjected to preloading conditions [16]. Hyperelastic models are usually used to represent the elastic behavior of rubbers undergoing large deformations [17, 18]. The reader may refer to a recent review article by [19]. Rubber-like materials exhibit also strain stiffening effects arising from limiting molecular chains extensibility; the simplest model of this type is due to Gent [20]. This model was used to represent the hyperelastic response of dielectric elastomers [21, 22]. Any bibliography on the subject of modeling rubber-like materials remains incomplete.

Under time-harmonic cyclic loading of increasing amplitude, filled elastomers show a decrease in storage modulus and a peak in loss modulus. This strain amplitude dependence is referred to as the Payne effect [23]. It was also shown that increasing the magnitude of strain for a given frequency will decrease the storage modulus monotonically; while the loss modulus increases reaching a maximum and then decreases. Both the storage and loss moduli decrease with a decrease in frequency for a fixed value of the amplitude of strain. Changes in dynamic moduli of filled elastomers seem to be related to the aggregation of filler particles into clusters and networks [24, 25]. The Payne effect has been a subject of numerous theories [26–30].

For instance, Lion and Kardelky [31] have developed a constitutive approach of finite viscoelasticity to represent the Payne effect in a continuum framework. The constitutive relationship is of integral form, in which the hyperelasticity is modeled by a neo-Hookean model, and the physical time variables are replaced by process-dependent intrinsic times. Then, the integral constitutive equation is transformed using a complicated computation into a differential equation similar to that of the classical Maxwell model of linear viscoelasticity. The differential constitutive model was linearized in order to compute the harmonic complex Young's modulus. Consequently, the loss modulus does not depend on the pre-stretch. This result seems not to be consistent with experimental results [32–37]. Recently, Narayan and Palade [38] have noticed a similarity between the behaviors of thixotropic fluids and filled elastomers. Therefore, the internal microstructure of a material can be characterized by a scalar parameter or second-order tensorial called 'structure parameter' originated from changes of the thixotropy, i.e. on changes of the viscosity [39, 40]. As a result, two approaches have been developed to model the Payne effect [38]. The first was based on a thermodynamic framework involving multiple natural configurations, while the second used the concept of 'structure parameter'. These authors have shown that, their models fall short of capturing the Payne effect apparently due to inappropriate choice of both the Helmholtz potential and rate of dissipation.

In experiments, the dynamic properties of elastomers can be described by the linearized steady-state response, i.e. small sinusoidal vibrations that are superimposed on large static deformation. The static stress is due to pre-strain and determined only by the last state of static deformation; so that; it does not depend on the deformation history. Consequently, the viscoelastic effects depend only on superimposed small vibratory deformations [41, 42]. Classically, viscoelastic models are derived from two approaches, namely integral and internal variables approaches. So, many models have been proposed for large strain viscoelasticity and their numerical implementation in finite elements (FE) [14, 34, 35, 43–58] to mention a few. Furthermore, the response of soft materials and biological tissues can be strain-rate dependent; the Helmholtz free energy function is written in terms of external thermodynamic state variables. As a result, the nominal stress is expressed as the sum of the hyperelastic nominal stress corresponding to the quasi-static equilibrium state, the viscous overstress arising from high strain rate and stress due to long time memory effect [59]. We point out that a detailed review of the various visco-hyperelastic modeling approaches as well as an external-state variable based models have been proposed by [60]. One drawback of the integral models is probably their high number of material parameters, which are necessary to reproduce experiments. In engineering designs, the constitutive models should be simple with few material parameters, which are easy to determine using experimental data of standard testing. We recall that many models have extended the generalized Maxwell (GM) model of the linear viscoelasticity to finite strains. The GM model can capture almost a full frequency range of the experimental data. However, it produces local 'oscillations' [61, 62], and existing literature on this subject does not provide a convincing interpretation of this phenomenon. To overcome this difficulty, the fractional derivative-type models seem to be more attractive because the relaxation functions of elastomeric materials can be predicted with fewer material parameters.

In the present work, we have developed the integral version of the constitutive equation, i.e. equation (24) of Bouzidi et al [63], in which the thixotropy was not accounted for. To model the Payne effect, the so-called intrinsic time or shift time concept was used, which arises from a change of the material microstructure. Additionally, the Mooney–Rivlin hyperelastic model (instead of the neo-Hookean one) was considered to investigate the coupling between material nonlinearity and relaxation. So, the range of large strains, i.e. the order of 100% could be reached. As a result, the constitutive relationship was linearized in order to highlight the effects of pre-strain, the magnitude of dynamic strain, and frequency. Thereby, special attention will be paid to the small harmonic vibrations that are superimposed on a large static deformation; this loading is relevant in many engineering designs.

This paper is organized as follows: in section 2, we briefly summarize the basics of the visco-hyperelastic model already presented and analyzed in previous works. The third section covers the linearization of the proposed model, whereas section 4 is devoted to the computation of the harmonic complex Young's modulus in the pre-deformed configuration. The identification of the material parameters and checking predictions of the model are presented in section 5.

Let F = ∂x/∂X be the deformation gradient tensor, where X is the position vector of a material particle in the undeformed configuration and x is the corresponding position in the deformed configuration.

The right and left Cauchy–Green strain tensors are denoted by C = F^T F and B = FF^T , respectively. The principal invariants of C (or B) are given by ${I}_{1}={I}_{1}\left(\mathbf{C}\right)={I}_{1}\left(\mathbf{B}\right)=\mathrm{t}\mathrm{r}(\mathbf{C})$, ${I}_{2}={I}_{2}\left(\mathbf{C}\right)={I}_{2}\left(\mathbf{B}\right)=\frac{1}{2}\left({I}_{1}^{2}-\mathrm{t}\mathrm{r}({\mathbf{C}}^{2})\right)$, ${I}_{3}={I}_{3}\left(\mathbf{C}\right)={I}_{3}\left(\mathbf{B}\right)=\mathrm{det}(\mathbf{C})={\left(\mathrm{det}\enspace \mathbf{F}\right)}^{2}={J}^{2}$ [17].$\mathbf{E}=\frac{1}{2}\left(\mathbf{C}-\mathbf{I}\right)$ and $\mathbf{e}=\frac{1}{2}\left({\mathbf{C}}^{-1}-\mathbf{I}\right)$ are the Green–Lagrange strain tensor and Piola strain tensor, respectively, and I is the identity tensor.

Combining both the one-dimensional fractional linear constitutive equation of standard linear solid (SLS) and the concept of dual variables [64], Bouzidi et al [63] obtained the following constitutive equation

$$\begin{equation}\mathbf{S}(t)=-p\left(t\right){\mathbf{C}}^{-1}\left(t\right)+\underset{0{}}{\overset{t}{\int }}K(t-s){\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(s)\right]\mathrm{d}s-2{\mu }_{0}{\xi }^{-1}\underset{0{}}{\overset{t}{\int }}r(t-s)\frac{\mathrm{d}\mathbf{e}(s)}{\mathrm{d}s}\mathrm{d}s,\end{equation} \tag{ 1 }$$

where $\mathbf{S}\left(t\right)$ is the second Piola–Kirchhoff stress tensor; ${\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(t)\right]=2\frac{\partial W}{\partial \mathbf{C}(t)}$ is the second Piola–Kirchhoff stress tensor arising from the theory of hyper elasticity, W = W(I₁, I₂) is the strain-energy density function for an incompressible and isotropic material, and $p\left(t\right)$ is a Lagrange multiplier that was introduced to ensure the incompressibility assumption.${\xi }^{-1}={\left({\tau }_{\mathrm{C}}/{\tau }_{\mathrm{R}}\right)}^{\alpha } > 1$ is a material parameter, τ_R and τ_C are the relaxation time and creep time, respectively. The functions $K(t)={t}^{-1}{\left(t/{\tau }_{\mathrm{R}}\right)}^{\alpha }{E}_{\alpha ,\alpha }\left(-{\left(t/{\tau }_{\mathrm{R}}\right)}^{\alpha }\right)$ and $r(t)={E}_{\alpha ,1}\left(-{\left(t/{\tau }_{\mathrm{R}}\right)}^{\alpha }\right)$ are the relaxation functions. Recall that the generalization of the Mittag–Leffler function is defined by the series expansion ${E}_{\alpha ,\phi }(z)={\sum }_{k=0}^{\infty }\left({z}^{k}/{\Gamma}\left(\alpha k+\phi \right)\right)$ where $\alpha > 0,\enspace \phi > 0,\left\vert z\right\vert < \infty$ and ${\Gamma}\left(z\right)=\underset{0{}}{\overset{\infty }{\int }}{t}^{z-1}\mathrm{exp}(-t)\mathrm{d}t$ is the Gamma function [65]. We note that forz > 0, 0 < α ⩽ 1 and ϕ ⩾ α, the generalized Mittag–Leffler function on the negative axis E_α,ϕ (−z) is completely monotone.

Based on the relationship between the relaxation functions $r\left(t\right)$ and $K\left(t\right)$ [65], we may write

$$\begin{equation}K(t-s)=\frac{\mathrm{d}}{\mathrm{d}s}r(t-s).\end{equation} \tag{ 2 }$$

Notice that the dimension of the function $K\left(t\right)$ is 1/time and $r\left(t\right)$ is a scalar function. Figure 1 shows the plots of $r\left(t\right)$ and $t\enspace K\left(t\right)$. It can be seen that the function $K\left(t\right)$ presents a weak singularity at t = 0, and points out that $\underset{0{}}{\overset{\infty }{\int }}K\left(s\right)\mathrm{d}s=1$.

Zoom In Zoom Out Reset image size

By changing the variable in equation (1), carrying out the integration by part, and substituting equation (2) into equation (1); equation (1) can be recast as follows

$$\begin{align}\hfill \mathbf{S}(t)& =-p\left(t\right){\mathbf{C}}^{-1}\left(t\right)+{\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(t)\right]\hfill \\ \hfill & \quad \hspace{3pt}\underset{{\widetilde{\boldsymbol{S}}}^{\left(\mathrm{v}\right)}(t)}{\underbrace{-{\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(0)\right]r(t)-\underset{0{}}{\overset{t}{\int }}r(t-s)\frac{\mathrm{d}}{\mathrm{d}s}\left({\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(s)\right]\right)\mathrm{d}s-2{\mu }_{0}{\xi }^{-1}\underset{0{}}{\overset{t}{\int }}r(t-s)\frac{\mathrm{d}\mathbf{e}(s)}{\mathrm{d}s}\enspace \mathrm{d}s}},\hfill \end{align} \tag{ 3 }$$

where ${\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(t)\right]$ is the immediate contribution (elastic behavior) as it supports the contribution of the deformation at actual time, and ${\widetilde{\boldsymbol{S}}}^{\left(\mathrm{v}\right)}\left(t\right)$ is the viscous stress.

We note that the integral terms depend on the past history of the material, and the viscous stress is expressed as the sum of three components. The first integral of equation (3) represents a coupling between hyperelasticity (the material nonlinearity) and relaxation. Thereby, in a relaxation test, when time is tending towards 'infinity', the stress response at the equilibrium state is ${\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(t)\right]$. As a consequence, this stress is meaningful for the long term highlighting the equilibrium stress. The second integral of equation (3) arises from the so-called fading memory; it represents the instantaneous elasticity at zero-time. We recall that the principle of fading memory was introduced by Coleman and Noll [66] in the sense that the current or recent configurations have a much greater influence on the present state of stress than configurations that existed in the distant past. Valanis [67] showed that fading memory is a consequence of the positive rate of irreversible entropy production. Haupt and Lion [48] proved that the relaxation function of the Mittag–Leffler is consistent with fading memory concept in agreement with the theory of Truesdell and Noll [68]. Moreover, the exponential function for the relaxation is compatible with the principle of thermodynamics.

In order to describe the frequency dependence properties of rubber-like materials, the evolution of material microstructure can be incorporated using the concept of intrinsic time, so-called, z(t) which is driven by the deformation history but not by temperature. Indeed, for thermo-rheological simple materials, the physical meaning of shift time arises from the TTS (time-temperature superposition) principle. The latter states that the effects of time and temperature on the mechanical properties of a viscoelastic material are equivalent [69]. Notice that McGuirt and Lianis [70] generalized the concept of thermo-rheological simple materials to finite deformations in the framework of the theory of thermodynamics of materials with memory.

As a result, we replace the physical time variables, t and s by the variables z and x, [31, 71, 72], and equation (3) becomes

$$\begin{align}\hfill \mathbf{S}(t)& =-p\left(t\right){\mathbf{C}}^{-1}\left(t\right)+{\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(t)\right]\hfill \\ \hfill & \quad \hspace{3pt}\underset{{\widetilde{\boldsymbol{S}}}^{\left(\mathrm{v}\right)}(t)}{\underbrace{-{\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(0)\right]\enspace r\left(z(t)\right)-\hspace{-2pt}\underset{0{}}{\overset{z(t)}{\int }}\hspace{-2pt}r(z-x)\frac{\mathrm{d}}{\mathrm{d}x}\left({\mathbf{S}}^{\left(\mathrm{e}\right)}\left[\mathbf{C}(x)\right]\right)\mathrm{d}x-2{\mu }_{0}{\xi }^{-1}\underset{0{}}{\overset{z(t)}{\int }}r(z-x)\frac{\mathrm{d}\mathbf{e}(x)}{\mathrm{d}x}\enspace \mathrm{d}x}},\hfill \end{align} \tag{ 4 }$$

where $r(z)={E}_{\alpha ,1}\left(-{\left(z/{\tau }_{\mathrm{R}}\right)}^{\alpha }\right)$.

The intrinsic time variable z(t) depends on the structural variable, so-called q(t) as follows

$$\begin{equation}\frac{\mathrm{d}z}{\mathrm{d}t}=1+bq(t)\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\;\;z(0)=0,\end{equation} \tag{ 5 }$$

where b ⩾ 0 is a material parameter.

The variable q(t) is a phenomenological measure in the current state of the material's microstructure, which is determined by the evolution equation [31]

$$\begin{equation}\frac{\mathrm{d}q(t)}{\mathrm{d}t}+\frac{1}{\varsigma }q(t)=\frac{1}{\varsigma }\left({\tau }^{\beta }\sqrt{\frac{2}{3}}{\Vert}\mathbf{F}(t)\left\{\frac{{\mathrm{d}}^{\beta }}{\mathrm{d}{t}^{\beta }}\mathbf{e}(t)\right\}{\mathbf{F}}^{T}(t){\Vert}\right)\end{equation} \tag{ 6 }$$

with the initial condition q(0) = 0 and ${\Vert}\mathbf{X}{\Vert}={\left(\mathbf{X}\cdot \mathbf{X}\right)}^{1/2}=\sqrt{{\sum }_{i,k=1}^{3}\;{X}_{ik}^{2}}$.

According to Lion [73], the parameter ς > 0 is a relaxation time that describes the dynamic behavior of the material's microstructure, and the constant τ = 1 s is introduced for dimensional reasons.

The evolution law given by equation (5) is assumed to be linear function of a structural variable q(t). A more general differential equation was used in [51, 56] governing the intrinsic time

$$\begin{equation}{\dot{z}}_{k}(t)={g}_{k}\left.\mathbf{E},\dot{\mathbf{E}},\enspace {q}_{1},\dots \dots .{q}_{L}\right.\geqslant 0\end{equation} \tag{ 7 }$$

where k = 1, ...‥, M.

The solving of equation (7) could require long computation times [56]. Also, it is difficult to provide physical mechanisms justifying a relatively large number of model parameters.

Let us linearize the constitutive equation in the neighborhood of a static pre-deformation. Following Lion and Kardelky [31], the deformation gradient can be decomposed into a time-independent and a time-dependent part (see figure 2):

$$\begin{equation}\mathbf{F}(t)=\mathbf{f}(t){\mathbf{F}}_{0},\end{equation} \tag{ 8 }$$

where F₀ and f(t) are respectively the static part and the dynamic part of the deformation gradient; F₀ = const, f(t) = I + h(t) and the magnitude of the incremental displacement gradient h is assumed to be sufficiently small, so $\delta =\underset{-\infty \leqslant t\prec \infty }{\mathrm{max}}\left({\Vert}\mathbf{h}\left(t\right){\Vert}\right)\ll 1$.

Zoom In Zoom Out Reset image size

We note R the difference between a second order tensor Y and his linearization form, we have ${\Vert}\mathbf{R}{\Vert}\leqslant M{\delta }^{2}$ for δ → 0 and we write $\mathbf{R}=\mathbf{O}\left({\delta }^{2}\right)$, where $\mathbf{O}$ is the Landau symbol [31].

Based on these assumptions, the linearization of all kinematical tensors reads as

$$\begin{equation}{\mathbf{F}}^{-1}={\mathbf{F}}_{0}^{-1}\left(\mathbf{I}-\mathbf{h}\right)+\mathbf{O}\left({\delta }^{2}\right)\approx {\mathbf{F}}_{0}^{-1}\left(\mathbf{I}-\mathbf{h}\right)\end{equation} \tag{ 9 }$$

$$\begin{equation}\mathbf{C}={\mathbf{C}}_{0}+2{\mathbf{F}}_{0}^{T}{\mathbf{E}}_{L}{\mathbf{F}}_{0}+\mathbf{O}\left({\delta }^{2}\right)\approx {\mathbf{C}}_{0}+2{\mathbf{F}}_{0}^{T}{\mathbf{E}}_{L}{\mathbf{F}}_{0}\end{equation} \tag{ 10 }$$

$$\begin{equation}\mathbf{B}={\mathbf{B}}_{0}+\mathbf{h}{\mathbf{B}}_{0}+{\mathbf{B}}_{0}{\mathbf{h}}^{T}+\mathbf{O}\left({\delta }^{2}\right)\approx {\mathbf{B}}_{0}+\mathbf{h}{\mathbf{B}}_{0}+{\mathbf{B}}_{0}{\mathbf{h}}^{T}\end{equation} \tag{ 11 }$$

$$\begin{equation}{\mathbf{C}}^{-1}={\mathbf{C}}_{0}^{-1}-2{\mathbf{F}}_{0}^{-1}{\mathbf{E}}_{L}{\mathbf{F}}_{0}^{T-1}+\mathbf{O}\left({\delta }^{2}\right)\approx {\mathbf{C}}_{0}^{-1}-2{\mathbf{F}}_{0}^{-1}{\mathbf{E}}_{L}{\mathbf{F}}_{0}^{T-1}\end{equation} \tag{ 12 }$$

$$\begin{equation}{\mathbf{B}}^{-1}={\mathbf{B}}_{0}^{-1}-{\mathbf{h}}^{T}{\mathbf{B}}_{0}^{-1}-{\mathbf{B}}_{0}^{-1}\mathbf{h}+\mathbf{O}\left({\delta }^{2}\right)\approx {\mathbf{B}}_{0}^{-1}-{\mathbf{h}}^{T}{\mathbf{B}}_{0}^{-1}-{\mathbf{B}}_{0}^{-1}\mathbf{h}\end{equation} \tag{ 13 }$$

$$\begin{equation}\mathbf{e}={\mathbf{e}}_{0}-{\mathbf{F}}_{0}^{-1}{\mathbf{E}}_{L}{\mathbf{F}}_{0}^{T-1}+\mathbf{O}\left({\delta }^{2}\right)\approx {\mathbf{e}}_{0}-{\mathbf{F}}_{0}^{-1}{\mathbf{E}}_{L}{\mathbf{F}}_{0}^{T-1}\end{equation} \tag{ 14 }$$

$$\begin{equation}\dot{\mathbf{e}}=-{\mathbf{F}}_{0}^{-1}{\dot{\mathbf{E}}}_{L}{\mathbf{F}}_{0}^{T-1}+\mathbf{O}\left.{\delta }^{2}\right.\approx -{\mathbf{F}}_{0}^{-1}{\dot{\mathbf{E}}}_{L}{\mathbf{F}}_{0}^{T-1},\end{equation} \tag{ 15 }$$

where ${\mathbf{C}}_{0}={\mathbf{F}}_{0}^{T}{\mathbf{F}}_{0}$ _, ${\mathbf{B}}_{0}={\mathbf{F}}_{0}{\mathbf{F}}_{0}^{T}$ _, ${\mathbf{e}}_{0}=\frac{1}{2}\left({\mathbf{C}}_{0}^{-1}-\mathbf{I}\right)$ and ${\mathbf{E}}_{L}=\frac{1}{2}\left(\mathbf{h}+{\mathbf{h}}^{T}\right)$ is the incremental strain tensor.

The linearization version of the evolution equation (6) is given by

$$\begin{equation}\frac{\mathrm{d}q(t)}{\mathrm{d}t}+\frac{1}{\varsigma }q(t)=\frac{1}{\varsigma }\left({\tau }^{\beta }\sqrt{\frac{2}{3}}{\Vert}\frac{{\mathrm{d}}^{\beta }}{\mathrm{d}{t}^{\beta }}{\mathbf{E}}_{L}(t){\Vert}\right).\end{equation} \tag{ 16 }$$

The elastic potential of Mooney–Rivlin is considered in the constitutive model of equation (1), so

$$\begin{equation}W={c}_{10}\left({I}_{1}-3\right)+{c}_{01}\left({I}_{2}-3\right)=\frac{{\mu }_{0}}{2}\left[\chi \left({I}_{1}-3\right)+\left(1-\chi \right)\left({I}_{2}-3\right)\right],\end{equation} \tag{ 17 }$$

where ${\mu }_{0}=2\left({c}_{10}+{c}_{01}\right)$ is the shear modulus at infinitesimal deformations and 0 < χ ⩽ 1 is a dimensionless constant; if χ = 1 then the neo-Hookean model is recovered.

The Cauchy stress tensor can be inferred from equation (4) as follows

$$\begin{equation}\mathbf{\sigma }(t)=-p(t)\mathbf{I}+{\mathbf{\sigma }}^{(1)}(t)+{\mathbf{\sigma }}^{(2)}(t)+{\mathbf{\sigma }}^{(3)}(t),\end{equation} \tag{ 18 }$$

where

$$\begin{align*}\hfill {\mathbf{\sigma }}^{(1)}(t)& =2\left({c}_{10}\mathbf{B}(t)-{c}_{01}{\mathbf{B}}^{-1}(t)\right)-2\left(\mathbf{F}(t)\enspace \left[{c}_{10}\mathbf{I}+{c}_{01}\left({I}_{1}(\mathbf{C}(0))\mathbf{I}-\mathbf{C}(0)\right)\right]\enspace {\mathbf{F}}^{T}(t)\right)r(z(t))\hfill \end{align*}$$

$$\begin{equation}{\mathbf{\sigma }}^{(2)}(t)=-2{\mu }_{0}{\xi }^{-1}\mathbf{F}(t)\left(\underset{0{}}{\overset{z(t)}{\int }}r(z-x)\frac{\mathrm{d}\mathbf{e}(x)}{\mathrm{d}x}\mathrm{d}x\right){\mathbf{F}}^{T}(t)\end{equation} \tag{ 19 }$$

$$\begin{equation*}{\mathbf{\sigma }}^{(3)}(t)=-2{c}_{01}\mathbf{F}(t)\left(\underset{0{}}{\overset{z(t)}{\int }}r(z-x)\frac{\mathrm{d}}{\mathrm{d}x}\left[{I}_{1}\left(\mathbf{C}(x)\right)\mathbf{I}-\mathbf{C}(x)\right]\mathrm{d}x\right){\mathbf{F}}^{T}(t).\end{equation*}$$

The second term i.e. σ ⁽¹⁾(t) of the constitutive equation embedded both the transient and steady-state responses; the long-term viscous stress (the term for fading memory), i.e. σ ⁽²⁾(t) depends on the strain rate of the superposed state and the overall effect of the past history of the material, whereas r(t) is actually a relaxation function, which could be related to the rearrangement of polymer chains. Indeed, the transient network theory of Tanaka and Edwards [74] states that polymeric chains are divided into permanent networks and transient ones, the latter being rearranged upon deformation. Therefore, the number of chains (active) that contribute to the polymeric network has evolved with time. Accordingly, the viscoelastic response is related to the rearrangement of polymer chains, that is, detachment of active chains or attachment of dangling chains to temporary junctions. Noting that the stress relaxation time is longer compared to the time it takes for the internal microstructure to rearrange. The stress σ ⁽³⁾(t) shows a coupling between the integral term and the material nonlinearity arising from the Mooney Rivlin model $\left({c}_{01}\ne 0\right)$. This stress could be related to the free-fluctuating of the deformed network due to the applied pre-strain. Recalling that the elastic and relaxation responses of a deformed rubber are a priori physically distinct, the former being caused by a change in the equilibrium configuration free energy of the rubber network and the latter being associated with the kinetics of network chain configuration rearrangements.

For sufficiently large times (in stationary state), practically t tends to infinity, so that$r\left(\infty \right)=0$ and equation (18) becomes

$$\begin{align}\hfill \mathbf{\sigma }(t)& =-p(t)\mathbf{I}+\underset{{\mathbf{\sigma }}^{(\mathrm{e}\mathrm{q})}(t)}{\underbrace{2\left({c}_{10}\mathbf{B}(t)-{c}_{01}{\mathbf{B}}^{-1}(t)\right)}}\hfill \\ \hfill & \hspace{3pt}\quad \underset{{\mathbf{\sigma }}^{(\mathrm{v})}(t)}{\underbrace{-2\mathbf{F}(t)\left(\underset{0{}}{\overset{\infty }{\int }}r(z-x)\frac{\mathrm{d}}{\mathrm{d}x}\left({c}_{01}\left[{I}_{1}\left(\mathbf{C}(x)\right)\mathbf{I}-\mathbf{C}(x)\right]+{\mu }_{0}{\xi }^{-1}\mathbf{e}(x)\right)\mathrm{d}x\right){\mathbf{F}}^{T}(t)}}.\hfill \end{align} \tag{ 20 }$$

The linearization of the first invariant scalar of the tensor C is as follows

$$\begin{equation}{I}_{1}\left.\mathbf{C}\right.\approx {I}_{1}\left.{\mathbf{C}}_{0}\right.+2\enspace \mathrm{t}\mathrm{r}\left.{\mathbf{F}}_{0}^{T}{\mathbf{E}}_{L}{\mathbf{F}}_{0}\right..\end{equation} \tag{ 21 }$$

Using these approximations and omitting the higher-order terms in the displacement gradient, we may obtain the linearized version of the constitutive equation with respect to the intrinsic time as follows

$$\begin{align}\hfill \mathbf{\sigma }(t)& =-p(t)\mathbf{I}+\underset{{\mathbf{\sigma }}^{(0)}}{\underbrace{2\left({c}_{10}{\mathbf{B}}_{0}-{c}_{01}{\mathbf{B}}_{0}^{-1}\right)}}+\underset{{\mathbf{\sigma }}^{\left(\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\right)}(t)}{\underbrace{2{c}_{10}\left(\mathbf{h}(t){\mathbf{B}}_{0}+{\mathbf{B}}_{0}{\mathbf{h}}^{T}(t)\right)+2{c}_{01}\left({\mathbf{h}}^{T}(t){\mathbf{B}}_{0}^{-1}+{\mathbf{B}}_{0}^{-1}\mathbf{h}(t)\right)}}\hfill \\ \hfill & \hspace{3pt}\quad \underset{{\mathbf{\sigma }}^{\left(\mathrm{v}\right)}(t)}{\underbrace{-4{c}_{01}{\mathbf{B}}_{0}\underset{0{}}{\overset{\infty }{\int }}r(z-x)\left(\mathrm{t}\mathrm{r}\left({\mathbf{F}}_{0}^{T}\frac{\mathrm{d}\left({\mathbf{E}}_{L}(x)\right)}{\mathrm{d}x}{\mathbf{F}}_{0}\right)\mathbf{I}-\frac{d\left({\mathbf{E}}_{L}(x)\right)}{\mathrm{d}x}{\mathbf{B}}_{0}\right)\mathrm{d}x+2{\mu }_{0}{\xi }^{-1}\underset{0{}}{\overset{\infty }{\int }}r(z-x)\frac{\mathrm{d}\left({\mathbf{E}}_{L}(x)\right)}{\mathrm{d}x}\enspace \mathrm{d}x}}\hfill \end{align} \tag{ 22 }$$

We observe that the 'equilibrium' stress ${\mathbf{\sigma }}^{\left(\mathrm{e}\mathrm{q}\right)}(t)$ is a sum of the static equilibrium stress, i.e. ${\mathbf{\sigma }}^{\left(0\right)}$ that depends only on the pre-deformation, i.e. F₀ and of the incremental part of the equilibrium stress, i.e. ${\mathbf{\sigma }}^{\left(\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\right)}\left(t\right)$, which represents the coupling effects of the basic state and superposed state at the instant t. Pointing out that ${\mathbf{\sigma }}^{\left(0\right)}$ and ${\mathbf{\sigma }}^{\left(\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\right)}\left(t\right)$ have been highlighted based on another theory [51]. The stress arising from the first integral depends on the past history of the material, the strain rate, and the strain-energy density. The second integral term depends on the strain rate and the past history of the material as in linear viscoelasticity.

Let us consider the simple tension test defined by the deformation gradient

$$\begin{equation}\mathbf{F}(t)=\lambda (t)\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right)+\frac{1}{\sqrt{\lambda (t)}}\left({\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right),\end{equation} \tag{ 23 }$$

where λ(t) is the stretch ratio.

The deformation gradient F₀ of the static pre-deformation and the incremental deformation gradient f(t) are respectively given by

$$\begin{equation}{\mathbf{F}}_{0}={\lambda }_{0}\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right)+\frac{1}{\sqrt{{\lambda }_{0}}}\left({\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right)\end{equation} \tag{ 24 }$$

and

$$\begin{equation}\mathbf{f}(t)=\left(1+{\varepsilon }_{L}(t)\right)\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right)+\frac{1}{\sqrt{1+{\varepsilon }_{L}(t)}}\left({\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right)\end{equation} \tag{ 25 }$$

then

$$\begin{equation}{\mathbf{B}}_{0}={\lambda }_{0}^{2}\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right)+\frac{1}{{\lambda }_{0}}\left({\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right)\end{equation} \tag{ 26 }$$

and

$$\begin{equation}{\mathbf{B}}_{0}^{-1}={\lambda }_{0}^{-2}\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right)+{\lambda }_{0}\left({\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right).\end{equation} \tag{ 27 }$$

The linearization of ${\left(1+{\varepsilon }_{L}\right)}^{-1/2}\approx 1-\frac{1}{2}{\varepsilon }_{L}$, which leads to the following set of equations

$$\begin{equation}\mathbf{f}(t)\approx \left(1+{\varepsilon }_{L}(t)\right)\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right)+\left(1-\frac{1}{2}{\varepsilon }_{L}(t)\right)\left({\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right)\end{equation} \tag{ 28 }$$

$$\begin{equation}\mathbf{h}(t)=\mathbf{f}(t)-\mathbf{I}\approx {\varepsilon }_{L}(t)\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}-\frac{1}{2}\left[{\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right]\right)\end{equation} \tag{ 29 }$$

$$\begin{equation}{\mathbf{E}}_{L}(t)\approx {\varepsilon }_{L}(t)\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}-\frac{1}{2}\left[{\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right]\right).\end{equation} \tag{ 30 }$$

Under simple tension (or compression), the Cauchy stress tensor is given by

$$\begin{equation}\mathbf{\sigma }(t)=\sigma (t)\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right).\end{equation} \tag{ 31 }$$

We obtain by substituting the equations (26), (27) and (29)–(31) into equation (22):

$$\begin{align}\hfill {\sigma }_{11}(t)=\sigma (t)& =-p+2\left({c}_{10}{\lambda }_{0}^{2}-{c}_{01}{\lambda }_{0}^{-2}\right)+4\left({c}_{10}{\lambda }_{0}^{2}+{c}_{01}{\lambda }_{0}^{-2}\right){\varepsilon }_{L}\left(t\right)\hfill \\ \hfill & \quad +\left(4{c}_{01}{\lambda }_{0}+2{\mu }_{0}{\xi }^{-1}\right)\underset{0{}}{\overset{\infty }{\int }}r(z-x)\frac{\mathrm{d}{\varepsilon }_{L}(x)}{\mathrm{d}x}\enspace \mathrm{d}x\hfill \end{align} \tag{ 32 }$$

$$\begin{align}\hfill {\sigma }_{22}(t)={\sigma }_{33}(t)=0& =-p+2\left({c}_{10}{\lambda }_{0}^{-1}-{c}_{01}{\lambda }_{0}\right)-2\left({c}_{10}{\lambda }_{0}^{-1}+{c}_{01}{\lambda }_{0}\right){\varepsilon }_{L}\left(t\right)\hfill \\ \hfill & \quad -\left(2{c}_{01}\left(2{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+{\mu }_{0}{\xi }^{-1}\right)\underset{0{}}{\overset{\infty }{\int }}r(z-x)\frac{\mathrm{d}{\varepsilon }_{L}(x)}{\mathrm{d}x}\enspace \mathrm{d}x.\hfill \end{align} \tag{ 33 }$$

By eliminating the Lagrange multiplier, i.e. p, we may obtain the Cauchy stress

$$\begin{align}\hfill \sigma (t)& =2\left({c}_{10}+{c}_{01}{\lambda }_{0}^{-1}\right)\left({\lambda }_{0}^{2}-{\lambda }_{0}^{-1}\right)+2\left({c}_{10}\left(2{\lambda }_{0}^{2}+{\lambda }_{0}^{-1}\right)+{c}_{01}\left(2{\lambda }_{0}^{-2}+{\lambda }_{0}\right)\right){\varepsilon }_{L}\left(t\right)\hfill \\ \hfill & \quad +\left(2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right)\underset{0{}}{\overset{\infty }{\int }}r(z-x)\frac{\mathrm{d}{\varepsilon }_{L}(x)}{\mathrm{d}x}\enspace \mathrm{d}x.\hfill \end{align} \tag{ 34 }$$

We will consider the first Piola–Kirchhoff tensor acting on the statically pre-deformed configuration as follows

$$\begin{equation}{\mathbf{T}}_{\mathrm{R}0}=\mathbf{\sigma }{\mathbf{f}}^{T-1}\approx \mathbf{\sigma }\left(\mathbf{I}-{\mathbf{h}}^{T}\right).\end{equation} \tag{ 35 }$$

In the case of simple tension, we have

$$\begin{equation}{\left({T}_{\mathrm{R}0}\right)}_{11}={T}_{\mathrm{R}0}\approx \sigma \left(1-{\varepsilon }_{L}\right).\end{equation} \tag{ 36 }$$

So

$$\begin{align}\hfill {T}_{\mathrm{R}0}(t)& =2\left({c}_{10}+{c}_{01}{\lambda }_{0}^{-1}\right)\left({\lambda }_{0}^{2}-{\lambda }_{0}^{-1}\right)+2\left({c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right){\varepsilon }_{L}\left(t\right)\hfill \\ \hfill & \quad +\left(2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right)\underset{0{}}{\overset{\infty }{\int }}r(z-x)\frac{\mathrm{d}{\varepsilon }_{L}(x)}{\mathrm{d}x}\enspace \mathrm{d}x.\hfill \end{align} \tag{ 37 }$$

Let us evaluate equation (37) for stationary sinusoidal excitations in the form

$$\begin{equation}{\varepsilon }_{L}(t)={\Delta}\varepsilon \enspace \mathrm{sin}\left(\omega t\right),\end{equation} \tag{ 38 }$$

where ω is the angular frequency and Δɛ the strain magnitude $\left({\Delta}\varepsilon \ll 1\right)$.

Inserting equation (30) into equation (16), we obtain the one dimensional formulation of the evolution equation

$$\begin{equation}\dot{q}(t)=\frac{1}{\varsigma }\left({\tau }^{\beta }\left\vert \frac{{\mathrm{d}}^{\beta }}{\mathrm{d}{t}^{\beta }}{\varepsilon }_{L}(t)\right\vert -q(t)\right).\end{equation} \tag{ 39 }$$

According to Lion and Kardelky [31], the solution of equation (39) is

$$\begin{equation}q(t)=\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }.\end{equation} \tag{ 40 }$$

Inserting equation (40) into equation (5), and after integration, we obtain the intrinsic time variable

$$\begin{equation}z(t)=\left(1+b\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }\right)t=at.\end{equation} \tag{ 41 }$$

Noticing that z(t) is a linear function of time and dynamic amplitude, but nonlinearly depends on the frequency.

Let us develop equation (37) by substituting equation (41) into the Green function $r\left(z-x\right)$ as follows

$$\begin{equation}r\left(z-x\right)={E}_{\alpha ,1}\left(-{\left(\frac{z-x}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }\right)={E}_{\alpha ,1}\left(-{\left(\frac{t-s}{{\tau }_{0}}\right)}^{\alpha }\right)={r}_{0}\left(t-s\right),\end{equation} \tag{ 42 }$$

where τ₀ = τ_R/a.

As a result, we obtain

$$\begin{equation}\begin{gathered}{T}_{\mathrm{R}0}(t)=2\left({c}_{10}+{c}_{01}{\lambda }_{0}^{-1}\right)\left({\lambda }_{0}^{2}-{\lambda }_{0}^{-1}\right)+2\left({c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right){\varepsilon }_{L}\left(t\right)\\ +\left(2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right)\underset{0{}}{\overset{\infty }{\int }}{r}_{0}(t-s)\frac{\mathrm{d}{\varepsilon }_{L}(s)}{\mathrm{d}s}\mathrm{d}s\end{gathered}.\end{equation} \tag{ 43 }$$

After calculations, we get

$$\begin{equation}\begin{aligned}\hfill {T}_{\mathrm{R}0}(t)& =2\left({c}_{10}+{c}_{01}{\lambda }_{0}^{-1}\right)\left({\lambda }_{0}^{2}-{\lambda }_{0}^{-1}\right)+2\left[{c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right]{\Delta}\varepsilon \enspace \mathrm{sin}\left(\omega t\right)\hfill \\ \hfill & \quad +\left[2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right]{\Delta}\varepsilon \left(\frac{{\left(\omega {\tau }_{0}\right)}^{\alpha }}{1+{\left(\omega {\tau }_{0}\right)}^{2\alpha }+2{\left(\omega {\tau }_{0}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}\right)\hfill \\ \hfill & \quad \times \left[\left(\mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\omega {\tau }_{0}\right)}^{\alpha }\right)\mathrm{sin}\left(\omega t\right)+\mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\mathrm{cos}\left(\omega t\right)\right].\hfill \end{aligned}\end{equation} \tag{ 44 }$$

Equation (44) can be written as follows

$$\begin{equation}{T}_{\mathrm{R}0}(t)=2({c}_{10}+{c}_{01}{\lambda }_{0}^{-1})({\lambda }_{0}^{2}-{\lambda }_{0}^{-1})+{E}^{\prime }{\Delta}\varepsilon \enspace \mathrm{sin}(\omega t)+{E}^{\prime \prime }{\Delta}\varepsilon \enspace \mathrm{cos}(\omega t).\end{equation} \tag{ 45 }$$

The first term of equation (45) represents the static stress due to the static pre-deformation λ₀, whereas the second term is the dynamic response depending on λ₀, ω and Δɛ. So, the storage and loss moduli are given by

$$\begin{align}\hfill \begin{aligned}\hfill {E}^{\prime }\left({\lambda }_{0},\omega ,{\Delta}\varepsilon \right)& =2\left[{c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right]+\left[2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right]\hfill \\ \hfill & \quad \times \left(\frac{{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{2\alpha }}{1+{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{2\alpha }+2{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}\right)\hfill \end{aligned}\end{align} \tag{ 46 }$$

and

$$\begin{equation}{E}^{{\prime\prime}}\left({\lambda }_{0},\omega ,{\Delta}\varepsilon \right)=\left[2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right]\left(\frac{{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)}{1+{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{2\alpha }+2{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}\right),\end{equation} \tag{ 47 }$$

where

$$\begin{equation}a=a\left(\omega ,{\Delta}\varepsilon \right)=\left(1+b\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }\right).\end{equation} \tag{ 48 }$$

Consequently, both the storage and the loss moduli depend on the magnitude of pre-strain in the case of Mooney–Rivlin model.

For the neo-Hookean model, i.e. c₀₁ = 0, we get

$$\begin{equation}{E}^{\prime }\left({\lambda }_{0},\omega ,{\Delta}\varepsilon \right)={\mu }_{0}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{\mu }_{0}{\xi }^{-1}\left(\frac{{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{2\alpha }}{1+{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{2\alpha }+2{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}\right)\end{equation} \tag{ 49 }$$

and

$$\begin{equation}{E}^{{\prime\prime}}\left(\omega ,{\Delta}\varepsilon \right)=3{\mu }_{0}{\xi }^{-1}\left(\frac{{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)}{1+{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{2\alpha }+2{\left(\frac{\omega {\tau }_{\mathrm{R}}}{a}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}\right).\end{equation} \tag{ 50 }$$

For this case, the loss modulus, i.e. E'' does not depend on the static pre-deformation, i.e. λ₀. We recall that the complex moduli in [31] arise from solving a fractional differential equation (with respect to the time variable) governing the overstress, which requires many computations. However, the proposed description was proved to be more flexible because it necessitates fewer computations and the Lion and Kardelky model [31] could be recovered.

Following Nashif et al [37], we may introduce the loss factor of the material, i.e.$\eta =\widetilde{\eta }\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)$ as follows

$$\begin{equation}{E}^{\ast }\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)={E}^{\prime }\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)\left(1+i\eta \right).\end{equation} \tag{ 51 }$$

We may observe that, if λ₀ → 1, we then get

$$\begin{equation}\underset{{\lambda }_{0}\to 1}{\mathrm{lim}}\enspace {E}^{\ast }\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)={E}^{\ast }\left({\Delta}\varepsilon ,\omega \right).\end{equation} \tag{ 52 }$$

The real and imaginary parts of the harmonic complex Young's modulus E* can be inferred from equations (46) and (47) by setting λ₀ = 1. As a result, the strain-dependence of harmonic moduli that arises from this approach could be compared, for instance, with moduli in [26]; however, it does not the scope of this paper.

Let us consider the loss factor that is given by

$$\begin{equation}\eta \left({\lambda }_{0},{{\Omega}}_{n}\right)=\frac{{E}^{{\prime\prime}}\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)}{{E}^{\prime }\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)}=\frac{\left[\left(1-\chi \right)\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\xi }^{-1}\right]\left(\mathrm{sin}\left(\alpha \pi /2\right)\right){{\Omega}}_{n}^{\alpha }}{A+B{{\Omega}}_{n}^{\alpha }+C{{\Omega}}_{n}^{2\alpha }},\end{equation} \tag{ 53 }$$

where

$$\begin{equation*}A=\chi \left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}-3{\lambda }_{0}^{-2}\right)+3{\lambda }_{0}^{-2}\end{equation*}$$

$$\begin{equation*}B=\left[\chi \left(2{\lambda }_{0}^{2}-4{\lambda }_{0}+4{\lambda }_{0}^{-1}-5{\lambda }_{0}^{-2}\right)+4{\lambda }_{0}+5{\lambda }_{0}^{-2}+3{\xi }^{-1}\right]\left(\mathrm{cos}\left(\alpha \pi /2\right)\right)\end{equation*}$$

$$\begin{equation*}C=\chi \left({\lambda }_{0}^{2}-4{\lambda }_{0}+2{\lambda }_{0}^{-1}-2{\lambda }_{0}^{-2}\right)+4{\lambda }_{0}+2{\lambda }_{0}^{-2}+3{\xi }^{-1}\end{equation*}$$

and Ω_n = ωτ_R/a.

Figure 3 shows the plots of the loss factor without the Payne effect, i.e. a = 1 and Ω_n = ωτ_R = Ω in the framework of tension and compression mode loadings. In tensile mode, the predictions seem to be qualitatively consistent with the results of [37]. The loss factor decreases with increasing the amplitude of pre-deformation. It can be seen that the plots representing the loss factor are shifted by increasing the pre-strain; to our knowledge, this result was not highlighted in the literature. In compression, the loss factor decreases significantly due to stiffening of the material; the material seems to be stiffening more in compression. Besides, we will consider the influence of the Payne effect a ≠ 1 on the loss factor in the next section.

Zoom In Zoom Out Reset image size

To check the proposed model's capability of simulating the Payne effect, experimental data of Lion and Kardelky [31] and of Peng et al [36] have been used. Lion and Kardelky [31] have performed their experimental tests on CBR, the testing temperature was about 283 K and the pre-strain was ɛ₀ = −0.1. The magnitude of harmonic excitation of small magnitude varied between 0.1% and 5%, and the frequency between 10 Hz and 60 Hz. Peng et al [36] used a dynamic mechanical analyzer to measure the complex Young's modulus of chloroprene rubber (CR). The specimens were pre-stretched sinusoidally with an amplitude λ₀ = 1.05 at various frequencies (1–30 Hz), the dynamic amplitude ranges from 0.1% to 3%, and the temperature is 296 K.

The linearization of the model is based on the assumption that the magnitude of the incremental displacement gradient h is small in comparison to 1, i.e. δ ≪ 1.

The maximum value of the norm of h(t) is given by

$$\begin{equation*}\delta =\underset{0\leqslant t\leqslant T}{\mathrm{max}}\left({\Vert}\mathbf{h}\left(t\right){\Vert}\right)=\underset{0\leqslant t\leqslant T}{\mathrm{max}}\sqrt{{\sum }_{i,k=1}^{3}{h}_{ik}^{2}\left(t\right)},\end{equation*}$$

$$\begin{equation*}\text{where}\;\;\;\mathbf{h}(t)={\varepsilon }_{L}(t)\left({\mathbf{e}}_{1}\otimes {\mathbf{e}}_{1}\right)+\left(\frac{1}{\sqrt{1+{\varepsilon }_{L}(t)}}-1\right)\left({\mathbf{e}}_{2}\otimes {\mathbf{e}}_{2}+{\mathbf{e}}_{3}\otimes {\mathbf{e}}_{3}\right).\end{equation*}$$

The maximum value of the magnitude of harmonic excitation is Δɛ = 0.05 and we have $\mathrm{sin}\left(\omega t\right)\in \left[-1,\enspace 1\right]$, then $\delta =\mathrm{max}\left({\left[{\left({\varepsilon }_{L}(t)\right)}^{2}+2{\left(\frac{1}{\sqrt{1+{\varepsilon }_{L}(t)}}-1\right)}^{2}\right]}^{1/2}\right)\approx 0.0621$. which can be neglected in comparison to 1

In these tests, the equilibrium state is assumed to be reached after a sufficiently long time (when the Mullins effect is removed), and before applying dynamic small vibrations on the specimen. We first observe that the material constants c₁₀ and c₀₁ may, in principle, be determined from simple tension tests or/and equi-biaxial tension tests under the quasi-static conditions [75, 76]. These constants should be valid to characterize the multi-axial mechanical behavior in the steady-state. Thus, it is difficult to determine the material constants c₁₀ and c₀₁ from experimental data of the simple tension tests, which are performed in the pre-deformed configuration. It should be pointed out that both equilibrium and instantaneous responses are elastic and are linked to the viscosity domain. As a result, the present identification would induce imprecision on the values of these material parameters. Notice that an agreement with the uni-axial test data does not guarantee that the proposed model will correctly predict actual material behavior for multi-axial stress fields. For this reason, a biaxial test needs to be performed in order to obtain experimental data, which could be used to validate the model predictions. In fact, for incompressible materials, the biaxial test is the most general since any triaxial stress state may be reduced to a biaxial one (without changing the state of deformation) by the superposition of an appropriate hydrostatic pressure [77]. Besides, the relaxation function r(t) does not depend on the pre-stretch, i.e. λ₀. However, for filled rubbers, experiments show that the relaxation depends on the static deformation [78]. Thereby, a correction factor was introduced in order to account for the non-separability of time and deformation in the relaxation function [34].

The material parameters were estimated by using the standard minimization technique. To this end, one needs to define the objective function ${\Lambda}\left(\kappa \right)$ as the quadratic squared sum of the residuals, i.e. the difference between the experimental moduli data and the corresponding theoretical values:

$$\begin{equation}{\Lambda}\left(\kappa \right)={w}_{1}\sum\limits _{i=1}^{M}{\left({{E}^{\prime }}_{\mathrm{M}\mathrm{o}\mathrm{d}}\left(\kappa ;\enspace i\right)-{{E}^{\prime }}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(\kappa ;\enspace i\right)\right)}^{2}+{w}_{2}\sum\limits _{j=1}^{{M}^{\prime }}{\left({{E}^{{\prime\prime}}}_{\mathrm{M}\mathrm{o}\mathrm{d}}\left(\kappa ;\enspace j\right)-{{E}^{{\prime\prime}}}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(\kappa ;\enspace j\right)\right)}^{2},\end{equation} \tag{ 54 }$$

where M and M' are the numbers of data points, $\left(\kappa \right)={\left({\mu }_{0},\enspace \chi ,\enspace \xi ,\enspace \alpha ,\enspace \beta ,\enspace {\tau }_{\mathrm{R}},\enspace b\right)}^{T}$ the unknown parameters set, and ${E}_{\mathrm{M}\mathrm{o}\mathrm{d}}^{\prime }\left(\kappa ;\enspace i\right)$ and ${E}_{\mathrm{M}\mathrm{o}\mathrm{d}}^{{\prime\prime}}\left(\kappa ;\enspace j\right)$ are the theoretical moduli and computed according to equations (46) and (47). The weighted factors w₁ and w₂ are introduced in order to properly scale the residuals terms, i.e. $\left({{E}^{\prime }}_{\mathrm{M}\mathrm{o}\mathrm{d}}\left(\kappa ;\enspace i\right)-{{E}^{\prime }}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(\kappa ;\enspace i\right)\right)$ and $\left({{E}^{{\prime\prime}}}_{\mathrm{M}\mathrm{o}\mathrm{d}}\left(\kappa ;\enspace j\right)-{{E}^{{\prime\prime}}}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(\kappa ;\enspace j\right)\right)$. When the two residual terms have different orders of magnitude, the values of weight factors should not be equal. It can be seen that the magnitude of the dissipation modulus is approximate to 1/3 for [31] and to 1/4 for [36] of that of the storage modulus. So, the weighting factors are chosen as w₁ = 1, w₂ = 3 and w₁ = 1, w₂ = 4, respectively. The minimization problem consists of

$$\begin{equation}\underset{\left(\kappa \right)}{\mathrm{min}}\enspace {\Lambda}\left(\kappa \right)\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \left(\kappa \right)\in \enspace {\Xi}.\end{equation} \tag{ 55 }$$

With ${\Xi}=\left\{\kappa :\enspace {\kappa }^{\mathrm{min}}\leqslant \kappa \leqslant {\kappa }^{\mathrm{max}}\right\}$ the solution space, $\left({\kappa }^{\mathrm{min}}\right)$ and $\left({\kappa }^{\mathrm{max}}\right)$ are the lower and upper bounds, respectively. The physical values of each material parameter were estimated to define the bounds of the vector $\left(\kappa \right)$.

Noting that global optimization or gradient-based local optimization methods are often used to calibrate the hyperelastic model constants [79, 80]. Apart from gradient-based methods, other gradient-free schemes are applied, for example, neural networks [81] or with a direct search method [82]. To obtain the model parameters of equation (54), we use a nonlinear least-square method (NLS), which requires an initial guess for the solution. We follow the gradient-based local optimization method, the objective function ${\Lambda}\left(\kappa \right)$ is minimized using the Lsqcurvefit implemented in the Optimization Toolbox of Matlab, and it is called the function fmincon and Lsqcurvefit. For instance, the combination of the Lsqcurvefit function and the trust-region-reflective method was used in the framework of visco-hyperelasticity [60] to obtain an optimal model parameter set corresponding to a global optimum. Accordingly, we have used a multi-start approach to running the optimization function with varying initial guesses. To this end, each parameter of the model should be bounded based on physical arguments, so that

$$\begin{equation}{\kappa }_{p=1\dots 7}^{\mathrm{min}}\leqslant {\kappa }_{p=1\dots 7}\leqslant {\kappa }_{p=1\dots 7}^{\mathrm{max}}.\end{equation} \tag{ 56 }$$

The constants of hyperelastic models must be determined carefully to obtain stable constitutive equations for all deformation states, and so to avoid numerical instabilities. Reasonable restrictions on material parameters based on physical arguments seem to be an unsolved question. Several inequalities such as the Backer–Ericksen inequalities, Coleman–Noll inequalities, strong ellipticity, polyconvexity, etc, have been proposed in the literature. The reader may refer to a recent paper for details [83]. Practically, we can assume that the model constants must meet the Drucker-stability criterion, i.e. the fourth-order elasticity tensor must be positive definite. Consequently, the choice for stable Mooney Rivlin model constants would be choosing them to be positive-valued [84]. Recently, Upadhyay et al [83] have proved that this choice is consistent based on the so-called T–C inequality. As a result, in equation (17), the Mooney Rivlin's constants are positive-valued. Moreover, for soft materials, the shear modulus is of the order of a few MPa, so that, we may assume that its upper value is 10 MPa, therefore; $0< {\mu }_{0}\leqslant 10\enspace \left(\mathrm{M}\mathrm{P}\mathrm{a}\right)$. Note that the material parameters μ₀ and χ must be determined from the data of both simple tension and equi-biaxial tension (or compression test).

The fractional derivatives theory imposes the values of the real number $0< \alpha < 1,0< \beta < 1$ and we have $0< \xi < 1$ [63]. We note that the relaxation behavior of filled rubbers is characterized by a very rapid decrease of the stress at the beginning of a hold time, followed by an extremely slow decay during several hours [71]. So, the values of the relaxation time, i.e. τ_R are assumed to be in the range of 10⁻³ s to 60 s. Moreover, the values of the material parameter b must be delimited; Lion and Kardelky [31] have estimated b of the order of 5 × 10⁴. We choose the upper bound for the parameter b equal to 10⁹ which can be considered a very great value (see appendix A).

The 'best' solution $\left({\kappa }^{\ast }\right)$ is that corresponding to the globally minimum. Practically, the number of 3000 iterations has been performed, and it was assumed as a criterion of the convergence to quasi-unique and 'stable' model parameters.

The quality of the solution is validated by computing the root-mean-square (RMS) as follows

$$\begin{equation}\mathrm{R}\mathrm{M}\mathrm{S}=\frac{1}{\sqrt{n}}\sqrt{\sum\limits _{i=1}^{i=n}{\left(\frac{{E}_{\mathrm{M}\mathrm{o}\mathrm{d}}\left(\kappa ;i\right)-{E}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(\kappa ;i\right)}{{E}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(\kappa ;i\right)}\right)}^{2}}.\end{equation} \tag{ 57 }$$

The RMS gives indications on how the model reproduces the experimental data. However, we cannot say anything about the quality of the parameters found [85]. Notice that the previous boundings exclude a priori to find numerical values without physical meanings of the model parameters.

Both the RMS corresponding to the storage and dissipation moduli were computed for each frequency. It is worth noting that, the experimental data used in this study are available in the literature. The determined values of material parameters are listed in tables 1 and 2, whereas the fitting results are shown in figures 4(a)–7(a). Figures 4(b)–7(b) show the relative error versus strain amplitude for different values of frequency. The storage and loss moduli give qualitatively good results concerning the shape of the two curves versus the strain magnitude and the frequency. These figures show that the storage modulus i.e. ${E}^{\prime }\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)$ decreases monotonously with strain, i.e. Δε, while the loss modulus, i.e. ${E}^{{\prime\prime}}\left({\lambda }_{0},{\Delta}\varepsilon ,\omega \right)$ displays a non-monotonic variation. Hence, they reflect the influence of the Payne effect. For the results of Kardelky and Lion, the worst fit for the loss modulus was achieved for low magnitudes of the dynamic strain Δɛ ⩽ 0.3% and low frequencies, i.e. 10 Hz–20 Hz, while the RMS for the storage modulus does not exceed 12%. For the results of Peng et al, the error is less than 20%. These high percentages of error arise from the linearity of the differential equation governing the intrinsic time variable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To improve the predictions of the storage and loss moduli, a nonlinear differential equation can be used for the intrinsic time. Therefore, the number of the model parameters must increase, and the computation of the moduli becomes very complicated (see appendix B).

For the data of Lion and Kardelky (2004), the investigated storage and loss moduli show that the RMS is found between (7.92%–11.55%) and (4.62%–19.61%) respectively. For the data of Peng et al (2020), the RMS of storage and loss moduli are varied respectively between (3.32%–10.1%) and (3.92%–10.74%). It can be seen that a relatively good agreement was obtained between theory and experiment.

Let us consider the influence of the Payne effect a ≠ 1 on the loss factor.

We have Ω_n = ωτ_R/a where $a=1+2\frac{b}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }$ and ω = 2πf is the radial frequency.

The plots of η versus log Ω_n are displayed in figures 8(a) and (b) for the two extreme values of the frequencies used in their experiments, i.e. 10 Hz and 60 Hz for the results of Lion and Kardelky (2004) and 1 Hz and 30 Hz for the results of Peng et al (2000). The curves seem to reproduce the master curve that can be compared to that from arising by varyingtemperatures.

Zoom In Zoom Out Reset image size

We recall that the data of uni-axial tensile alone are used for the calibration of the model parameters. As a result, these parameters should not ensure to predict the multi-axial stress state. Hereafter, the same set of model parameters was used to predict the first Piola–Kirchhoff stress in the pre-deformed configuration, i.e. equation (44). The plots are shown in figures 9 and 10 for the compression loading and in figures 11 and 12 for the tensile loading, the trends of the results reproduce the experimental data and seem to be qualitatively satisfactory and conforming with the results of the literature [56].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work, the Payne effect has been modeled based on the visco-hyperelastic model presented and analyzed in previous work. The proposed model was formulated on the grounds of the concept of dual variables in the framework of extension of the fractional constitutive equation of the SLS to finite strains. So, the proposed constitutive equation seems to include the Lion and Kardelky model as a special case. The present model was used to evaluate the harmonic complex Young's modulus with regards to the Payne effect. To this end, experimental data of literature were employed to check the model's capability. Indeed, the model was calibrated from experimental data of compression and tensile loading tests, which were performed for a pre-deformation, a wide range of frequencies, and small harmonic strains that were superimposed on a steady state. We have shown that the proposed description is more flexible, given that the computation is easier. Thereafter, the same set of parameters identified from the storage and loss moduli was used to predict the stress in the pre-deformed configuration. The present modeling has certain accuracy in predicting the Payne effect. Hence, the proposed constitutive equation could be suitable for modeling the mechanical behavior of unfilled rubbers. To improve the predictions of the storage and loss moduli, a nonlinear differential equation can be used for the intrinsic time. Therefore, the number of the model parameters must increase, and the computation of the moduli becomes very complicated. So, a balance should be found between the flexibility of the calculation of moduli and the type of differential equation. Besides, the predictions of the present model can be improved by considering the effect of pre-stretch on the relaxation function. The present model can also be implemented in FE codes in order to simulate the response of materials and structures.

We would like to thank the anonymous Reviewers for their detailed comments, which have significantly improved the content of this paper.

All data that support the findings of this study are included within the article (and any supplementary files).

The procedure of calibration of optimal model parameters is based on the gradient-based local optimization method. Therefore, the objective function is minimized by using the Levenberg–Marquardt algorithm implemented in the Optimization Toolbox of Matlab, and it is called the function fmincom and lsqcurvefit. Then, multi-start approaches, i.e. different initial guesses were used to obtain the global optimal parameter set.

We have chosen the Latin hypercube sampling procedure to generate N starting guesses. To do this, each parameter of the model should be bounded based on physical arguments. The values of the positive material parameter b must be delimited, so we write

$$\begin{equation*}0< {\mu }_{0}\leqslant 10\enspace \left(\mathrm{M}\mathrm{P}\mathrm{a}\right),\quad 0< \chi \leqslant 1,\quad 0< \alpha < 1,\quad 0< \beta < 1,\quad 0< \xi < 1,\quad 0< {\tau }_{\mathrm{R}}< 60\left(s\right),\quad 0\leqslant b\leqslant 1{0}^{9}.\end{equation*}$$

In order to find a global solution, we have used the lsqcurvefit function of Matlab together with MultiStart by varying initial guesses. The lhsdesign (N, 7) is used to generate a Latin hypercube sampling that contains N values of each of the seven (7) parameters. The N values of each colon of the matrix T are randomly distributed with one from each interval (0, 1/N), (1/N, 2/N), .........((N − 1)/N, 1). We create a matrix A that contains N initial guesses as follows

$$\begin{equation}{A}_{ij}={T}_{ij}\times l{b}_{ij}+\left(1-{T}_{ij}\right)\times u{b}_{ij},\end{equation} \tag{ A.1 }$$

where lb_ij and ub_ij represent respectively the lower and upper bound for a parameter κ_j .

We execute multiStart using the lsqcurvefit function as the local solver with N initial guesses given by the matrix A. The iterations of the NLS algorithm can stop when the number of iteration surpasses 3000.

The tables A.1 and A.2 give, respectively, the identified parameters of the experimental data of Lion and Kardelky (2004) and of Peng et al (2020) for N = 10, N = 100, N = 1000,N = 3000.

We can see that we find a similar value of these parameters $\left(\xi ,\enspace \alpha ,\enspace \beta ,\enspace {\tau }_{\mathrm{R}},\enspace b\right)$ for different values of iterations but different values for $\left({\mu }_{0},\enspace \chi \right)$, this latter may be determined under the quasi-static conditions.

We have identified the model parameters by considering other values for the upper bound of the parameter τ_R and the number of iterations is taken as N = 3000.

Let us consider the error R₀ that is defined by

$$\begin{equation}{R}_{0}=\sum\limits _{l=1}^{m}\left\{\sum\limits _{i=1}^{{M}_{l}}{\left({{E}^{\prime }}_{\mathrm{M}\mathrm{o}\mathrm{d}}\left(l,i\right)-{{E}^{\prime }}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(l,i\right)\right)}^{2}+\sum\limits _{j=1}^{{{M}^{\prime }}_{l}}{\left({{E}^{{\prime\prime}}}_{\mathrm{M}\mathrm{o}\mathrm{d}}\left(l,j\right)-{{E}^{{\prime\prime}}}_{\mathrm{e}\mathrm{x}\mathrm{p}}\left(l,j\right)\right)}^{2}\right\},\end{equation} \tag{ A.2 }$$

where m represents the number of frequencies used in the test.

The plots of the error R₀ are shown in figure A.1 versus the value of the upper bound of the relaxation time, i.e. ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$, in the framework of both the Kardelky and Lion and of Peng et alresults.

Zoom In Zoom Out Reset image size

The graph of figure A.1(a) shows that the function R₀ decreases for the low values of ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$ and it remains quasi-constant for ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}\geqslant 60\enspace \mathrm{s}$. In figure A.1(b), the variation of values of ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$ do not affect the error R₀.

In tables A.3 and A.4, we present a comparison between the numerical values of material parameters that arise from identification and of the values of RMS in the context of choosing two upper bounds for the relaxation time, ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=60\enspace \mathrm{s}$ and ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=400\enspace \mathrm{s}$.

Table A.3. Comparison of the consistency of computed numerical values of material parameters depending on the value of ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$.

Kardelky and Lion results			Peng et al results
Parameters	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=60\enspace \mathrm{s}$	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=400\enspace \mathrm{s}$	Parameters	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=60\enspace \mathrm{s}$	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=400\enspace \mathrm{s}$
μ0 (MPa)	1.956	2.042	μ0 (MPa)	1.897	1.941
χ	0.707	0.903	χ	0.603	0.375
ξ	0.179	0.182	ξ	0.188	0.202
α	0.451	0.451	α	0.349	0.349
β	0.249	0.249	β	0.424	0.424
τR (s)	59.997	399.994	τR (s)	0.0036	0.0036
b	3542 485.04	23 619 428.21	b	236.22	236.23

Table A.4. Comparison of the consistency of computed numerical values of RMS depending on the value of ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$.

RMS (E', E'') % Kardelky and Lion results			RMS (E', E'') % Peng et al results
Frequency (Hz)	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=60\enspace \mathrm{s}$	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=400\enspace \mathrm{s}$	Frequency (Hz)	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=60\enspace \mathrm{s}$	${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}=400\enspace \mathrm{s}$
10	(9.10, 19.61)	(9.10, 19.61)	1	(10.10, 3.92)	(10.10, 3.92)
20	(8.41, 13.48)	(8.41, 13.48)	5	(8.81, 10.74)	(8.81, 10.74)
30	(7.92, 6.74)	(7.92, 6.74)	10	(3.32, 7.47)	(3.32, 7.47)
40	(10.91, 4.62)	(10.91, 4.62)	20	(4.23, 9.59)	(4.23, 9.59)
50	(10.75, 6.1)	(10.75, 6.1)	30	(4.62, 10.63)	(4.62, 10.63)
60	(11.55, 6.95)	(11.55, 6.95)

One can see that the values of the parameters b and τ_R are sensitive on choosing the value of ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$ in the framework of Kardelky and Lion results, but their ratio, i.e. b/τ_R ≈ 5.90 × 10⁴ s⁻¹ is quasi-constant. Noting that the RMS values of both the storage and loss moduli do not depend on the value of ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$. However, no difference was highlighted in analyzing the Peng et alresults.

Let us consider the expressions of the storage and loss moduli equations (46) and (47) in order to explain in depth the difference arising from Kardelky and Lion results and Peng et alresults.

We recall that the storage and loss moduli are functions of the variable ωτ₀ that is defined by $\omega {\tau }_{0}=\frac{\omega }{{\tau }_{\mathrm{R}}^{-1}+\left(\frac{b}{{\tau }_{\mathrm{R}}}\right)\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }}$.

If ${\tau }_{\mathrm{R}}^{-1}\ll \left(\frac{b}{{\tau }_{\mathrm{R}}}\right)\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }$ then we obtain $\omega {\tau }_{0}\approx \frac{\omega }{\left(\frac{b}{{\tau }_{\mathrm{R}}}\right)\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }}$.

As a result, the variable ωτ₀ depends on the ratio $\left(b/{\tau }_{\mathrm{R}}\right)$.

For the ratio corresponding to the Kardelky and Lion results, the values of parameters b and τ_R cannot be unique.

Otherwise, if ${\tau }_{\mathrm{R}}^{-1}$ is not negligible compared to$\left(\frac{b}{{\tau }_{\mathrm{R}}}\right)\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }$, the values of the parameters b and τ_R seem to be stable, that is the case of the results of Peng et al for which the results do not depend on the value of ${\tau }_{\mathrm{R}}^{\mathrm{M}\mathrm{a}\mathrm{x}}$.

In experimental testing, the material behavior of filler-reinforced rubber under dynamic deformation shows a weakly dependence on the frequency of deformation but a pronounced dependence on the dynamic strain amplitude. This well known Payne effect can be described as follows

• Increasing frequency f leads to an increase in the storage and loss moduli,
• The storage modulus decreases with increasing dynamic strain amplitude while the loss modulus runs through a distinct maximum.

Mathematically, we can write

$$\begin{equation}\forall \enspace {\Delta}\varepsilon :\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}f} > 0\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \forall \enspace f:\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}{\Delta}\varepsilon }< 0\end{equation} \tag{ B.1 }$$

and

$$\begin{align}\hfill \forall \enspace {\Delta}\varepsilon :\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}f} > 0\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \forall \enspace f:{\left.\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}{\Delta}\varepsilon }\right\vert }_{{\Delta}\varepsilon ={\Delta}{\varepsilon }_{0,f}}=0,\quad {\left.\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}{\Delta}\varepsilon }\right\vert }_{{\Delta}\varepsilon < {\Delta}{\varepsilon }_{0,f}} > 0,\\ \quad {\left.\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}{\Delta}\varepsilon }\right\vert }_{{\Delta}\varepsilon > {\Delta}{\varepsilon }_{0,f}}< 0\end{align} \tag{ B.2 }$$

where Δɛ_0,f represents the value of strain amplitude in which the loss modulus is maximum.

The storage and loss moduli given by equations (46) and (47) depend on seven parameters and we have

$$\begin{equation}\forall \enspace {\Delta}\varepsilon :\quad \frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}f}=A\times \left[2\pi \alpha {\tau }_{\mathrm{R}}{\left(\omega {\tau }_{0}\right)}^{\alpha -1}\left(1+\left(1-\beta \right)b\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }\right)\right]\end{equation} \tag{ B.3 }$$

and

$$\begin{equation}\forall \enspace f:\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}{\Delta}\varepsilon }=A\times \left(\alpha \omega {\tau }_{\mathrm{R}}{\left(\omega {\tau }_{0}\right)}^{\alpha -1}\right)\times \left[-b\frac{2}{\pi }{\left(\omega \tau \right)}^{\beta }\right],\end{equation} \tag{ B.4 }$$

where

$$\begin{equation}{\tau }_{0}=\frac{{\tau }_{\mathrm{R}}}{a},\quad a=\left(1+b\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }\right)\end{equation} \tag{ B.5a }$$

and

$$\begin{equation}A=\frac{\left(2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right)\times \left[\mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\omega {\tau }_{0}\right)}^{2\alpha }\mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+2{\left(\omega {\tau }_{0}\right)}^{\alpha }\right]}{{\left[1+{\left(\omega {\tau }_{0}\right)}^{2\alpha }+2{\left(\omega {\tau }_{0}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)\right]}^{2}{\left(1+b\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }\right)}^{2}} > 0.\end{equation} \tag{ B.5b }$$

For any value of parameters ${\mu }_{0} > 0,\enspace 0< \chi \leqslant 1,\enspace 0< \xi < 1,\enspace \enspace 0< \alpha < 1,\enspace 0< \beta < 1,\enspace {\tau }_{\mathrm{R}} > 0,\enspace b > 0$, we have

$$\begin{equation}\forall \enspace {\Delta}\varepsilon :\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}f} > 0,\qquad \forall \enspace f:\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}{\Delta}\varepsilon }< 0.\end{equation} \tag{ B.6 }$$

For the loss modulus, we have

$$\begin{equation}\forall \enspace {\Delta}\varepsilon :\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}f}=B\times \left[2\pi \alpha {\tau }_{\mathrm{R}}{\left(\omega {\tau }_{0}\right)}^{\alpha -1}\left(1+\left(1-\beta \right)b\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }\right)\right]\times \left[1-{\left(\omega {\tau }_{0}\right)}^{2\alpha }\right]\end{equation} \tag{ B.7 }$$

and

$$\begin{equation}\forall \enspace f:\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}{\Delta}\varepsilon }=B\times \left(\alpha \omega {\tau }_{\mathrm{R}}{\left(\omega {\tau }_{0}\right)}^{\alpha -1}\right)\times \left[1-{\left(\omega {\tau }_{0}\right)}^{2\alpha }\right]\times \left[-b\frac{2}{\pi }{\left(\omega \tau \right)}^{\beta }\right],\end{equation} \tag{ B.8 }$$

where

$$\begin{equation}B=\frac{\left(2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right)\times \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)}{{\left[1+{\left(\omega {\tau }_{0}\right)}^{2\alpha }+2{\left(\omega {\tau }_{0}\right)}^{\alpha }\mathrm{cos}\left(\frac{\alpha \pi }{2}\right)\right]}^{2}{\left(1+b\frac{2}{\pi }{\Delta}\varepsilon {\left(\omega \tau \right)}^{\beta }\right)}^{2}} > 0.\end{equation} \tag{ B.9 }$$

It can be seen that if $\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}f} > 0$, i.e. $\left[1-{\left(\omega {\tau }_{0}\right)}^{2\alpha }\right] > 0$, implies that$\frac{\mathrm{d}{E}^{{\prime\prime}}}{\mathrm{d}{\Delta}\varepsilon }< 0$.

Then, if we impose the constraint of monotonicity of loss modulus versus frequency, the loss modulus decreases with dynamic strain amplitude which disagrees with the experiment. It can be concluded that for the present model, the two conditions cannot be verified simultaneously for low dynamic strains.

If we consider another type of evolution equation to represent the process dependence of the Payne effect, for example, the evolution equation for the rate of the intrinsic time given in [51] or in [56]:

$$\begin{equation}\dot{z}(t)=\delta \left(\tau \left\vert \dot{\varepsilon }\left(t\right)\right\vert +1\right)+\left(1-\delta \right)\left({\sum }_{i=1}^{l}{b}_{i}{q}_{i}+1\right)\end{equation} \tag{ B.10a }$$

$$\begin{equation}\dot{z}(t)=\delta \left({\tau }^{\gamma }\left\vert \frac{{\mathrm{d}}^{\gamma }\varepsilon (t)}{\mathrm{d}{t}^{\gamma }}\right\vert +1\right)+\left(1-\delta \right)\left(bq+1\right),\end{equation} \tag{ B.10b }$$

where 0 ⩽ δ ⩽ 1, and 0 < γ ⩽ 1.

The computation of the storage and loss moduli becomes very complicated and the number of material parameters increases.

In order to illustrate the effect of incorporation of the first term of equation (B.10) in the evolution equation for the rate of the intrinsic time, we consider a special choice given by equation (B.10a).

In the expression of the stress T_R0 (equation (37)), we can replace the integral term

$$\begin{equation}{{T}_{\mathrm{R}0}}_{(v)}={T}_{v}=\left(2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right)\underset{0{}}{\overset{\infty }{\int }}r(z-x)\frac{\mathrm{d}{\varepsilon }_{L}(x)}{\mathrm{d}x}\mathrm{d}x\end{equation} \tag{ B.11 }$$

by the differential equation

$$\begin{equation}{\tau }_{\mathrm{R}}^{-\alpha }{T}_{v}+\frac{{\mathrm{d}}^{\alpha }{T}_{v}}{\mathrm{d}{z}^{\alpha }}=m\enspace \frac{{\mathrm{d}}^{\alpha }\varepsilon }{\mathrm{d}{z}^{\alpha }}\end{equation} \tag{ B.12 }$$

where

$$\begin{equation}m=\left(2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right).\end{equation} \tag{ B.13 }$$

We note that the equation (B.11) is the solution for the fractional differential equation (B.12) in the stationary state.

This differential equation is formulated with reference to a new time scale z, which can either evolve faster or slower than the real time t, this formulation is also equivalent to a varying relaxation time in equation (B.12) [51].

We assume that

$$\begin{equation}{\tau }_{\mathrm{R},\mathrm{e}\mathrm{ff}}^{\alpha }=\frac{{\tau }_{\mathrm{R}}^{\alpha }}{{\left[\dot{z}(t)-\delta \tau \left\vert \dot{\varepsilon }\left(t\right)\right\vert \right]}^{\alpha }+\delta \tau \left\vert \dot{\varepsilon }\left(t\right)\right\vert }=\frac{{\tau }_{\mathrm{R}}^{\alpha }}{{\left[\left(1-\delta \right){\sum }_{i=1}^{l}\hspace{2pt}{b}_{i}{q}_{i}+1\right]}^{\alpha }+\delta \tau \left\vert \dot{\varepsilon }\left(t\right)\right\vert }\end{equation} \tag{ B.14 }$$

and the differential equation (B.12) becomes

$$\begin{equation}{\tau }_{\mathrm{R},\mathrm{e}\mathrm{ff}}^{-\alpha }{T}_{v}+\frac{{\mathrm{d}}^{\alpha }{T}_{v}}{\mathrm{d}{t}^{\alpha }}=m\enspace \frac{{\mathrm{d}}^{\alpha }\varepsilon }{\mathrm{d}{t}^{\alpha }}.\end{equation} \tag{ B.15 }$$

To calculate the dynamic moduli of the proposed model with the hypothesis of the effective relaxation time equation (B.14), we assume the stationary sinusoidal deformation

$$\begin{equation}{\varepsilon }_{L}(t)={\Delta}\varepsilon \enspace \mathrm{sin}\left(\omega t\right).\end{equation} \tag{ B.16 }$$

We have

$$\begin{equation}{T}_{v}(t)={\Delta}\varepsilon \left[A\enspace \mathrm{sin}\left(\omega t\right)+B\enspace \mathrm{cos}\left(\omega t\right)\right]\end{equation} \tag{ B.17 }$$

and

$$\begin{equation}{T}_{\mathrm{R}0}(t)=2\left({c}_{10}+{c}_{01}{\lambda }_{0}^{-1}\right)\left({\lambda }_{0}^{2}-{\lambda }_{0}^{-1}\right)+2\left({c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right){\varepsilon }_{L}\left(t\right)+{T}_{v}(t).\end{equation} \tag{ B.18 }$$

We use the method of Galerkin to calculate the functions A and B

$$\begin{equation}\underset{0{}}{\overset{2\pi /\omega }{\int }}\left({\tau }_{\mathrm{R},\mathrm{e}\mathrm{ff}}^{-\alpha }{T}_{v}+\frac{{\mathrm{d}}^{\alpha }{T}_{v}}{\mathrm{d}{t}^{\alpha }}\right){N}_{j}(t)\mathrm{d}t=m\underset{0{}}{\overset{2\pi /\omega }{\int }}\left(\frac{{\mathrm{d}}^{\alpha }\varepsilon }{\mathrm{d}{t}^{\alpha }}\right){N}_{j}(t)\mathrm{d}t.\end{equation} \tag{ B.19 }$$

For ${N}_{1}(t)=\mathrm{sin}\left(\omega t\right)$ and ${N}_{2}(t)=\mathrm{cos}\left(\omega t\right)$, we obtain a system of linear equations for the functions A and B

$$\begin{align}& \underset{0{}}{\overset{2\pi /\omega }{\int }}\left[{\tau }_{\mathrm{R}}^{-\alpha }\left({\left(\left(1-\delta \right){\sum }_{i=1}^{l}{b}_{i}{q}_{i}+1\right)}^{\alpha }+\delta \tau \left\vert \dot{\varepsilon }\left(t\right)\right\vert \right)\left(A\enspace \mathrm{sin}\left(\omega t\right)+B\enspace \mathrm{cos}\left(\omega t\right)\right)\right]\mathrm{sin}\left(\omega t\right)\mathrm{d}t\\ & \qquad +{\omega }^{\alpha }\underset{0{}}{\overset{2\pi /\omega }{\int }}\left[A\enspace \mathrm{sin}\left(\omega t+\frac{\alpha \pi }{2}\right)+B\enspace \mathrm{cos}\left(\omega t+\frac{\alpha \pi }{2}\right)\right]\mathrm{sin}\left(\omega t\right)\mathrm{d}t \\ & \qquad =m{\omega }^{\alpha }\underset{0{}}{\overset{2\pi /\omega }{\int }}\left[\mathrm{sin}\left(\omega t+\frac{\alpha \pi }{2}\right)\right]\mathrm{sin}\left(\omega t\right)\mathrm{d}t \end{align} \tag{ B.20 }$$

$$\begin{align} &\underset{0{}}{\overset{2\pi /\omega }{\int }}\left[{\tau }_{\mathrm{R}}^{-\alpha }\left({\left(\left(1-\delta \right){\sum }_{i=1}^{l}{b}_{i}{q}_{i}+1\right)}^{\alpha }+\delta \tau \left\vert \dot{\varepsilon }\left(t\right)\right\vert \right)\left(A\enspace \mathrm{sin}\left(\omega t\right)+B\enspace \mathrm{cos}\left(\omega t\right)\right)\right]\mathrm{cos}\left(\omega t\right)\mathrm{d}t\\ & \qquad +{\omega }^{\alpha }\underset{0{}}{\overset{2\pi /\omega }{\int }}\left[A\enspace \mathrm{sin}\left(\omega t+\frac{\alpha \pi }{2}\right)+B\enspace \mathrm{cos}\left(\omega t+\frac{\alpha \pi }{2}\right)\right]\mathrm{cos}\left(\omega t\right)\mathrm{d}t \\ & \qquad =m{\omega }^{\alpha }\underset{0{}}{\overset{2\pi /\omega }{\int }}\left[\mathrm{sin}\left(\omega t+\frac{\alpha \pi }{2}\right)\right]\mathrm{cos}\left(\omega t\right)\mathrm{d}t. \end{align} \tag{ B.21 }$$

After calculations, we obtain

$$\begin{equation}\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{4\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]A-\left[{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\right]B=m{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)\end{equation} \tag{ B.22 }$$

$$\begin{equation}\left[{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\right]A+\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{8\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]B=m{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right),\end{equation} \tag{ B.23 }$$

where

$$\begin{equation}\left(1-\delta \right){\sum }_{i=1}^{l}{b}_{i}{q}_{i}+1=c.\end{equation} \tag{ B.24 }$$

Solving these equations, we obtain

$$\begin{equation}A=m\frac{{\omega }^{2\alpha }+{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)\left[{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{8\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]}{\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{4\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{8\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]+{\left[{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\right]}^{2}}\end{equation} \tag{ B.25 }$$

$$\begin{equation}B=m\frac{{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\left[{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{4\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]}{\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{4\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{8\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]+{\left[{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\right]}^{2}}.\end{equation} \tag{ B.26 }$$

As a result, the dynamic moduli are given by

$$\begin{equation}\begin{aligned}\hfill & {E}^{\prime }\left({\lambda }_{0},\omega ,{\Delta}\varepsilon \right)\hfill \\ \hfill & \qquad =2\left[{c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right]\hfill \\ \hfill & \qquad \quad +m\frac{{\omega }^{2\alpha }+{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)\left[{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{8\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]}{\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{4\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{8\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]+{\left[{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\right]}^{2}}\hfill \end{aligned}\end{equation} \tag{ B.27 }$$

$$\begin{equation}{E}^{{\prime\prime}}\left({\lambda }_{0},\omega ,{\Delta}\varepsilon \right)=m\frac{{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\left[{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{4\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]}{\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{4\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]\left[{\omega }^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\frac{c}{{\tau }_{\mathrm{R}}}\right)}^{\alpha }+\frac{8\delta \left(\omega \tau \right)}{3\pi {\tau }_{\mathrm{R}}^{\alpha }}{\Delta}\varepsilon \right]+{\left[{\omega }^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)\right]}^{2}}.\end{equation} \tag{ B.28 }$$

If δ = 0 and l = 1 then, we obtain c = bq + 1 = a and ${\tau }_{\mathrm{R},\mathrm{e}\mathrm{ff}}^{\alpha }=\frac{{\tau }_{\mathrm{R}}^{\alpha }}{{\left(bq+1\right)}^{\alpha }}=\frac{{\tau }_{\mathrm{R}}^{\alpha }}{{a}^{\alpha }}={\tau }_{0}^{\alpha }$

In this case, the dynamic moduli reduce to the equations (46) and (47)

$$\begin{equation}\begin{aligned}\hfill {E}^{\prime }\left({\lambda }_{0},\omega ,{\Delta}\varepsilon \right)& =2\left[{c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right]+m\left(\frac{{\omega }^{2\alpha }+{\omega }^{\alpha }{\tau }_{0}^{-\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}{{\omega }^{2\alpha }+2{\omega }^{\alpha }{\tau }_{0}^{-\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\tau }_{0}^{-2\alpha }}\right)\hfill \\ \hfill \quad \quad \quad & =2\left[{c}_{10}\left({\lambda }_{0}^{2}+2{\lambda }_{0}^{-1}\right)+3{c}_{01}{\lambda }_{0}^{-2}\right]\hfill \\ \hfill & \quad +\left[2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right]\times \left(\frac{{\left(\omega {\tau }_{0}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\left(\omega {\tau }_{0}\right)}^{2\alpha }}{1+{\left(\omega {\tau }_{0}\right)}^{2\alpha }+2{\left(\omega {\tau }_{0}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}\right)\hfill \end{aligned}\end{equation} \tag{ B.29 }$$

$$\begin{align}{E}^{{\prime\prime}}\left({\lambda }_{0},\omega ,{\Delta}\varepsilon \right)&=m\left(\frac{{\omega }^{\alpha }{\tau }_{0}^{-\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)}{{\omega }^{2\alpha }+2{\omega }^{\alpha }{\tau }_{0}^{-\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)+{\tau }_{0}^{-2\alpha }}\right)\\ & =\left[2{c}_{01}\left(4{\lambda }_{0}-{\lambda }_{0}^{-2}\right)+3{\mu }_{0}{\xi }^{-1}\right]\times \left(\frac{{\left(\omega {\tau }_{0}\right)}^{\alpha }\enspace \mathrm{sin}\left(\frac{\alpha \pi }{2}\right)}{1+{\left(\omega {\tau }_{0}\right)}^{2\alpha }+2{\left(\omega {\tau }_{0}\right)}^{\alpha }\enspace \mathrm{cos}\left(\frac{\alpha \pi }{2}\right)}\right).\end{align} \tag{ B.30 }$$

An equilibrium between the flexibility of the moduli's computations and the differential equation of the evolution law should, therefore, be found.