Dynamic stiffness of chemically and physically ageing rubber vibration isolators in the audible frequency range

The constitutive equations of chemically and physically ageing rubber in the audible frequency range are modelled as a function of ageing temperature, ageing time, actual temperature, time and frequency. The constitutive equations are derived by assuming nearly incompressible material with elastic spherical response and viscoelastic deviatoric response, using Mittag-Leffler relaxation function of fractional derivative type, the main advantage being the minimum material parameters needed to successfully fit experimental data over a broad frequency range. The material is furthermore assumed essentially entropic and thermo-mechanically simple while using a modified William–Landel–Ferry shift function to take into account temperature dependence and physical ageing, with fractional free volume evolution modelled by a nonlinear, fractional differential equation with relaxation time identical to that of the stress response and related to the fractional free volume by Doolittle equation. Physical ageing is a reversible ageing process, including trapping and freeing of polymer chain ends, polymer chain reorganizations and free volume changes. In contrast, chemical ageing is an irreversible process, mainly attributed to oxygen reaction with polymer network either damaging the network by scission or reformation of new polymer links. The chemical ageing is modelled by inner variables that are determined by inner fractional evolution equations. Finally, the model parameters are fitted to measurements results of natural rubber over a broad audible frequency range, and various parameter studies are performed including comparison with results obtained by ordinary, non-fractional ageing evolution differential equations.

In general, physical ageing of rubber is a process where the molecular structure, initially brought out of thermodynamical equilibrium, spontaneously advances to regenerate its thermodynamical equilibrium state [6]. Physically, the molecular structure processes include volume alterations, trapping and freeing of polymer chain ends and other polymer chain reorganizations. Suitable thermodynamic variables to describe the processes include free volume and enthalpy [6,33]. The advantage of enthalpy is its ability to mutually include volumepreserving and non-preserving configurational molecular alterations [33]. There are various definitions of free volume of molecular structure of which the excess volume compared to that found in the crystalline state of the same rubber is the most suitable since it does not contain the concept of atom volume being physically difficult to define [33]. Free volume brought out of thermodynamical equilibrium by a sudden temperature alteration spontaneously advances to regenerate its thermodynamical equilibrium free volume state. This process is present throughout the temperature range including glassy, transition and rubbery regions. However, the process is fast in the rubbery region due to the short relaxation times of the configurational molecular processes involved, resulting in a fast free volume equilibrium state completion. On the contrary, the process is significantly slower in the transition and glassy regions due to longer relaxation times of the configurational molecular processes involved, resulting in a slow free volume equilibrium state completion. Moreover, the complexity of the molecular structure including the covalent cross-bonds and van der Waals bonds further delays the completion of the free volume equilibrium state [33]. Furthermore, the free volume deviation from its equilibrium state is frequently significantly larger than the spontaneous fluctuations of the molecular structure resulting in a further increase of the time to reach free volume equilibrium state in the transition and glassy regions [6]. The physical ageing time is the time to reach free volume equilibrium and depends on the temperature, its history, including other thermodynamical histories, and the molecular structure. The physical ageing process is the process to attain thermodynamical equilibrium, and since several mechanical, electrical and calorimetric properties depend on the free volume, those will also alter during the course of physical ageing. The memory of physical ageing is possible to erase by thermal rejuvenation whereby the rubber sample is typically heated well above the glass transition temperature for a sufficiently extended time [33]. There are several models developed for physical ageing of which the simplest relies on a linear differential equation describing the evolution of the free volume resulting in an exponential free volume decay and rises with time, depending on the history, while using a simple one-parameter relaxation time. A more elaborated model, known as Kohlrausch-Williams-Watts (KWW) function [25,38], applies an additional stretching parameter to account for the non-symmetric evolution of free volume observed between a jump-down and jump-up of the temperature history [3]. Furthermore, the relaxation time is possible to extend from a simple one-parameter into an exponential function depending on temperature, on a parameter specifying the departure of the system from its equilibrium state, on thermal expansion coefficient and on a nonlinear parameter specifying the relative contribution from temperature and structure on the relaxation time resulting in a model known as Kovacs-Aklonis-Hutchinson-Ramos (KAHR) model [26]. An analogue model, known as Tool-Narayanaswamy-Moynihan (TNM) model [29,32,36], applies a further parameter to the relaxation time dependence named as fictive temperature, being the minimum temperature for an aged rubber where a full free volume recovery is possible, instead of a parameter describing the departure of the system from the equilibrium as in the KAHR-model. Greiner and Schwarzl [13] apply Doolittle equation [8,10] showing an reversed exponential free volume dependence of the relaxation time of the free volume differential evolution equation, thus resulting in a nonlinear free volume differential evolution equation. However, the free volume evolution differential equations mentioned are all based on integer differential operators, while constitutive equations relating stress and strain that successfully fit measurement results of rubber frequently contain fractional differential operators; see for example Bagley and Torvik [1,2] and Kari et al. [22]. In this paper, a novel free volume evolution fractional differential equation is developed that is connected to the previously developed constitutive equations while using an identical equilibrium relaxation time since the physical origin for the free volume evolution and that of stress and strain response of the molecular structure is similar.
Chemical ageing by oxygen is generally an irreversible diffusion and reaction-driven process where the molecular network is permanently altered by link degeneration and reformation in contrast to physical ageing being a reversible process. Generally, chemical ageing is a slower process than the corresponding physical ageing since the operating temperature range for rubber is typically well above the glass transition temperature showing a fast physical ageing. The most common chemical ageing process is attributed to exposure of oxygen from outside resulting in diffusion and eventually oxygen reaction either damaging the network by scission or reforming new links. The chemical ageing properties are then reflected by altered stress-strain properties of the rubber. The chemical ageing is frequently modelled by molecular dynamics [5] and internal inner variables that are determined by inner differential evolution equations [18][19][20][21]27,30], recently extended into their nonlinear counterpart [16]. In this paper, novel inner differential evolution equations are developed for the inner variables while containing fractional differential operators in order to allow for a wider spectrum of relaxation times found in real molecular networks, thus giving a richer structure to the model.
The developed physical and chemical ageing models are implemented into constitutive stress and strain relations for rubber in the audible frequency range, and various parameter studies are performed including comparison with results obtained by ordinary, non-fractional ageing evolution differential equations. The developed models can be used to describe dynamic stiffness of ageing rubber vibration isolators in the audible frequency range. This is performed in part 2 of this paper [23].

Constitutive preliminaries
The rubber is assumed isotropic and nearly incompressible while obeying the principle of fading memory [9,37]. The analysis is confined to infinitesimal strain and pre-strain, usually being the condition for structureborne sound applications, where Fletcher-Gent [11] effect is negligible. A convolution integral, expressed as a constitutive relaxation relation at natural time t and temperature T while being physically aged a time t ph at temperature T and chemically aged a time t ch at temperature T ch , 1 is additively decomposed into a spherical part and a deviatoric part where σ is stress, u is displacement, the nearly incompressibility material constant a 1, μ ∞ is equilibrium elastic modulus, the relaxation function μ 1 obeys lim t→∞ μ 1 (t) = 0, ∇ is covariant derivative, the operators tr, div and dev are trace, divergence and deviator, respectively, defined as tr [ · ] = i : : : and where : : : denotes double contraction, i is second-order unit tensor, 1 = (i i + i i)/2 is fourth-order unit tensor, and the tensor products ⊗, and are defined through (a ⊗ b) : : : c = (b : : : c)a, (a b) : : : c = acb T and (a b) : : : c = ac T b T , where a, b and c are arbitrary second-order tensors and T denotes transpose. The rubber is essentially considered entropically elastic where ρ T = ρ T (t ph ) is actual density at temperature T and ∞ ρ T 0 = lim t ph →∞ ρ T 0 is thermodynamical equilibrium density at reference temperature T 0 . Furthermore, the rubber is assumed essentially thermo-mechanically simple with a kernel of fractional standard linear solid [24] showing a vanishing limes as t → ∞, where the non-dimensional relaxation intensity 1, μ ν and 0 < α ≤ 1 are material constants, h is Heaviside step function, the Mittag-Leffler function reads is Gamma function, and the William-Landel-Ferry (WLF) shift function reads using the method of reduced variables, where A 1 and A 2 are material constants, T = T − T 0 and the physical ageing shift function X ph > 0 is to be defined later while obeying X ph = 1 at the thermodynamical equilibrium state. The derived constitutive equations follow standard derivation, however, extended with physical and chemical ageing compared to those previously developed models, such as Kari et al. [22], the extension being the study focus of this paper. Temporal Fourier transformations where i is imaginary unit, ω is angular frequency, of the constitutive Eqs. (1) and (2) with (6), (7) and (8) and where (iX ph X T ω) α def = (iX ph X T ω + 0 + ) α = |X ph X T ω| α exp[sign(ω)iπα/2]. Consequently, the rubber is displaying an admissible behaviour with a finite instantaneous shear modulus where shear modulusμ(T ch , The constitutive relaxation relation for ordinary standard linear solid is obtained by letting α → 1 resulting in for the deviatoric part, while the corresponding spherical part remains the same and where the relaxation time is In the frequency domain with identical instantaneous shear modulus (12) as for its fractional counterpart.

Relaxation spectrum
To further investigate the relaxation function (7), its relaxation spectrum γ T is now considered, being defined through where τ is relaxation time. Spectrum is studied in numerous fields including visco-elasticity, crack growth, molecular dynamics and their relation to macroscopic properties. It may be interpretable as the characteristic molecular time response range, being extremely broad while making a logarithmic time scale log e (τ ) suitable. The relaxation spectrum (16) is possible to rewrite as where the relaxation spectrum γ T is obtained from the shear modulus μ T [14] as where denotes imaginary part, withμ while using (11) with a normalizing relaxation time τ T ≡ X ph X T μ ν / 1/α . Clearly, the spectrum is continuous while displaying a maximum of γ T max = tan(πα/2) /2π at τ = τ T . Obviously, the maximum value of relaxation spectrum for γ T is independent of temperature, including the shift factors X ph and X T , while its relaxation time locus is temperature dependent on those shift factors. Consequently, the normalizing relaxation time above is equivalent to the locus of maximum relaxation time and may therefore serve as the generalized relaxation time of the fractional derivative constitutive model presented in this paper for the deviatoric stress response. Finally, the relaxation spectrum for an ordinary standard linear solid α → 1 is where δ is Dirac delta function and relaxation time τ T is given by (14). Clearly, the relaxation spectrum is singular and non-continuous, thus resulting in a single relaxation time for an ordinary standard linear solid.

Physical ageing
The novel free volume evolution fractional differential equation for t ph ≥ 0 developed in this paper reads where f T is actual fractional free volume at temperature T , v T is actual specific volume at temperature T , cry v T is specific volume at crystalline state and temperature T , is specific volume at the thermodynamical equilibrium state and temperature T , ∞ f T = lim t ph →∞ f T is fractional free volume at the thermodynamical equilibrium state and temperature T , is the Caputo fractional derivative of f T and order 0 < α < 1 [7] with respect to physical ageing time t ph , and τ T is the actual relaxation time at temperature T , conveniently reading since the physical origin for the free volume evolution and that of stress and strain response of the molecular structure is similar and the generalized relaxation time for the latter may be interpreted as (25); see Sect. 2.1. Apparently, the fractional free volume evolution differential equation is a fractional derivative extension of the evolution equation given by Greiner and Schwarzl [13] using ordinary derivative. The reason for using Caputo fractional derivative instead of other fractional time derivative definition such as that of Riemann-Liouville is its ability to use ordinary, integer order instead of fractional order initial value condition, a useful property while studying fractional order differential equations. The physical ageing shift factor X ph is still unspecified and is considered next. To this end, the general derivation by Greiner and Schwarzl [13] to relate actual relaxation time τ T at temperature T to any relaxation time, including its thermodynamic equilibrium value ∞ τ T = lim t ph →∞ τ T at the same temperature, is followed. First, the relaxation time is assumed to follow the Doolittle equation [8,10] where A 3 and A 4 are material constants. In passing, it is noted that the relaxation time is depending on the fractional free volume, thus making the free volume evolution fractional differential Eq. (21) nonlinear. Next, the assumption made by Gibbs and DiMarzio [12] that the fractional free volume at thermodynamical equilibrium state is linearly dependent on absolute temperature within a broad temperature range including the rubbery region and a broad part of the transition region, is stated where the fractional free volume thermal expansion coefficient β rubber and β crystal are the thermal volume expansion coefficient of rubber at the rubbery and crystalline state, respectively. Furthermore, it is assumed that the relaxation time for the previously derived constitutive equations relating stress and strain follows the Doolittle Eq. (26) and that the temperature shift factor is derived at the thermodynamical equilibrium state giving and using (27) log 10 where the material constants in (9) are readily identified as and Hence, the temperature shift factor (WLF) is established relating the thermodynamic equilibrium state at temperature T to T 0 , which was already shown by Greiner and Schwarzl [13], Kovacs et al. [26] among others and is shown here for completeness. Next, the physical ageing factor is derived relating the actual state at temperature T to its corresponding thermodynamic equilibrium state at the same temperature This is rewritten, using Doolittle equation (26), (27), identities (32) and (33), into The thermodynamic shift from thermodynamic equilibrium state at reference temperature T 0 into actual thermodynamic state at reference temperature T 0 is achieved by multiplying the shift factors into reading using (35) and (9), where the evolution of fractional free volume is given by (21). Clearly, the physical ageing contribution to the relaxation function of the constitutive stress and strain relations is now determined. However, there is also a physical ageing affected factor in the constitutive stress and strain relations ρ T / ∞ ρ T 0 that relates the actual density to its thermodynamical equilibrium value at the reference temperature T 0 . This factor is considered next.
This factor is possible to multiplicatively split into two factors within the considered temperature range including rubbery and transition regions, and using (22) and (23), giving using (27) and (33), where the evolution of fractional free volume is given by (21).

Chemical ageing
Chemical ageing is most important for rubber since the operating temperature range is typically well above the glass transition temperature showing a fast, almost non-perceptible, physical ageing.

Network link scission
The novel network link scission inner variable fractional differential evolution equation for t ch ≥ 0 developed in this paper reads with initial condition q sci (0) = 0 while obeying 0 ≤ q sci ≤ 1 and d q sci / d t ch ≥ 0 and used in where D α sci is Caputo fractional derivative of order 0 < α sci < 1 (24) with respect to chemical ageing time t ch and μ ∞ 0 is equilibrium elastic modulus at reference temperature T 0 and at no chemical ageing; t ch = 0. The scission process is displaying an Arrhenius temperature dependence resulting in a network link scission relaxation time where τ sci 0 is network link scission relaxation time at reference temperature T 0 , E sci is network link scission activation energy per mole, and R = 8.314 J/mol K is the universal molar gas constant. The solution to (42) is as shown by, for example, Ishteva [17]. The developed network link scission inner variable fractional differential evolution equation is an extension of standard ordinary differential evolution equations [18][19][20][21]27,30] in order to allow for a wider spectrum of relaxation times found in real molecular network link scission processes, thus giving a richer structure to the model. In fact, the ordinary solution is redeemed while α sci → 1 resulting in

Network link reformation
The novel network link reformation inner variable fractional differential evolution equation for t ch ≥ 0 developed in this paper reads where τ ref 0 is network link reformation relaxation time at reference temperature T 0 and E ref is network link reformation activation energy per mole. The solution to (47) is The developed network link reformation inner variable fractional differential evolution equation is, like for the previous scission modelling, an extension of standard ordinary differential evolution equations [18][19][20][21]27,30] in order to allow for a wider spectrum of relaxation times found in real molecular network link reformation processes. Likewise, the ordinary solution is recovered while α ref → 1 resulting in Finally, the total evolution of equilibrium elastic modulus including both network link scission and reformation is reducing into its ordinary counterpart while α sci → 1 and

Resulting constitutive relations including physical and chemical ageing
The resulting spherical part of the constitutive relaxation relation, physically and chemically aged a time t ph and t ch at temperature T and T ch , respectively, is while its deviatoric part is where f T (t ph ), τ sci , τ ref and X ph ( f T (t ph ))X T are given by (21), (44), (49) and (37), respectively. Temporal Fourier transformations of (54) and (55) result in constitutive equations in frequency domain where bulk modulusκ and shear modulusμ treating physical and chemical time independent from time involved in the temporal Fourier transformation.

Results and discussion
The results are calculated and presented by means of MATLAB .

Material
The material applied in this study is unfilled sulphur cured standard Malaysian natural rubber, identical to that of Kari et al. [22], where equilibrium density ∞ ρ T 0 ≈ 984 kg/m 3 at reference temperature T 0 = 298 K and WLF material constants of (9) are determined to A 1 = 5.940 and A 2 = 151.6 K for T > 243 K, while the shift factor is manually set for T ≤ 243 K to, for example, X 213 K = 10 9.80 , using a dynamic mechanical thermal analyser and frequency-temperature shifts. Furthermore, the material constants of shear modulus (57) at no chemical ageing and at thermodynamical equilibrium state (X ph = 1) are optimized, minimizing the relative least square errors between measurement and modelling results over a broad frequency range, evaluated at reference temperature T 0 ; reading μ ∞0 = 8.25 × 10 5 N/m 2 , μ ν = 1.53 × 10 −5 s, α = 0.657 and = 276, in addition to the nearly incompressible material constant here set to a = 2.222 × 10 3 . The reader is referred to Kari et al. [22] regarding details over measurements and optimization performed including rubber material used.

Relaxation time, fractional free volume evolution and physical ageing
The volume expansion coefficient is set to β rubber = 6.60 × 10 −4 K −1 , being identical to that of pure vulcanized natural rubber, tabulated in Brandrup et al. [4] and used in Kari et al. [22], while the fractional free volume expansion coefficient is properly approximated [13,18] as β ≈ [2/3] β rubber = 4.40 × 10 −4 K −1 . The resulting equilibrium density, relaxation time and fractional free volume are given in Table 1. Apparently, the equilibrium density and fractional free volume differ about 8 and 280 %, respectively, while the equilibrium relaxation time is displaying a variation of approximately 11 decades over the studied temperature range of −60 to +60 • C. The WLF material constants A 1 and A 2 , subsequently applied in this paper, are shown in same Table 1. Clearly, the constants are temperature independent down to −25 • C and adjusted at temperature −60 • C, where A 1 is attuned to meet the manual fit result of shift function X 213 K to follow that of measurement. This is a plausible adjustment as shift function undergoes a general transformation from a WLF type of function above transition temperature T g to an Arrhenius type of function below transition temperature [33], being approximately T g ≈ −63 • C according to Brandrup et al. [4] for pure vulcanized natural rubber.
Next, fractional free volume evolution from Eq. (21) is determined for rubber material exposed by a temperature alteration while initially being at thermodynamical equilibrium state, thus constituting a typical temperature history involving physical ageing. Although being practically impossible to achieve, the temperature alteration is here considered immediate in order to enable a full exploration of physical ageing in all temperature alterations studied. In particular, the studied temperature histories are T (t ph ) = T 0 + (T k − T 0 )h(t ph ), where k = 1, 2, 3, 4 and T k = 213, 248, 273, 333 K. To this end, the numerical procedure specified in "Appendix" is applied.
The resulting evolution of fractional free volume f T for the studied temperature histories is shown in Fig. 1 versus a very broad physical ageing time t ph range, starting from 10 −14 s which is the same order as the period time of fastest atom-to-atom oscillations in the molecules to 10 6 s which is of order of weeks (more precisely, 11 and half days). Clearly, fractional free volume brought out of thermodynamical equilibrium by a sudden temperature alteration spontaneously advances to regenerate its thermodynamical equilibrium fractional free volume state, the latter represented by the horizontal lines in Fig. 1. The fractional free volume at physical Horizontal lines are equilibrium fractional free volume at those temperatures ageing time 0 + s is still identical to that of equilibrium fractional free volume at 25 • C although the temperature is altered from T 0 to T k , where k = 1, 2, 3, 4. The larger the difference between initial and final temperature (for T k < T 0 ), the longer is the ageing time to reach thermodynamical equilibrium fractional free volume state; more precisely, the physical ageing time to reach within ±1 % from its equilibrium free volume is 6.26 × 10 −9 , 8.72 × 10 −7 , 1.82 × 10 −4 and 4.71 × 10 3 s for final temperature +60, ±0, −25 and −60 • C. On the contrary, the larger the difference between initial and final temperature (again for T k < T 0 ), the higher is the absolute rate value of fractional free volume change at initial physical ageing. The main reason is that the larger the temperature difference, the larger is the initial out of thermodynamical equilibrium of the system and, thus, the larger is the molecular system driving force to make it into an equilibrium. From mathematical point of view, this property is reflected by Eq. (21) where fractional time derivative of fractional free volume is proportional to the difference between final and initial equilibrium fractional free volume. However, at longer physical ageing and lower temperatures, the molecular reconfigurations involved are significantly slower due to longer relaxation times, less free space available and further delays due to the complexity of the molecular structure including the covalent cross-bonds and van der Waals bonds. This is quite different from the temperature history with step to T 4 = 333 K, where the short relaxation time and amply free space available result in a very fast fractional free volume completion. In passing, it is observed from Fig. 1 that fractional free volume evolution involves a wide spectrum of time scales due to not only the nonlinear relaxation time dependence on fractional free volume via Doolittle equation (26), but also as a result of applying Caputo fractional time derivative in the differential evolution equation causing itself a broad range of time scales. In conclusion, the physical ageing time to reach sufficiently close to free volume equilibrium is very short, almost non-perceptible for all temperatures studied expect for the temperature history with step to T 1 = 213 K, where physical ageing time is significantly longer and needs in general to be accounted for, while for higher temperatures the fractional free volume evolution using as a simple step function directly to its equilibrium value is in general sufficient. The resulting evolution of relaxation time τ T for the studied temperature histories is shown in Fig. 2 versus the same broad physical ageing time t ph range as for the corresponding evolution of fractional free volume. Clearly, the relaxation time spans over a very broad scale, thus displaying the strong nonlinearity of Eq. (21) where a small perturbation in fractional free volume results in a large alteration of relaxation time via Doolittle equation (26) and where the horizontal lines correspond to thermodynamical equilibrium relaxation time at each temperature studied. In particular, the covered short relaxation time displayed for temperature steps studied up to T 4 = 333 K numerically explains the short, almost non-perceptible physical ageing time for equilibrium state completion of the fractional free volume evolution previously studied in Fig. 1. Likewise,  Fig. 3 Evolution of relaxation spectrum γ 213 K versus relaxation time τ at physical ageing time 0 + , 10 −9 , 10 −4 , 10 1 and 10 6 s, associated with the fractional free volume f 213 K evolution starting at equilibrium fractional free volume corresponding to temperature +25 • C and ending at equilibrium value at −60 • C the long relaxation time displayed for temperature steps studied down to T 1 = 213 K numerically explains the very long physical ageing time for equilibrium state completion of the fractional free volume evolution also previously studied in Fig. 1.
The resulting evolution of relaxation spectrum γ 213 K for temperature history with step to T 1 = 213 K is shown in Fig. 3 versus relaxation time τ at physical ageing time 0 + , 10 −9 , 10 −4 , 10 1 and 10 6 s. Clearly, the peak relaxation spectrum of relaxation function μ 1 is independent of physical ageing time due to the normalization performed for μ 1 in Eq. (2), with special case of Eq. (16), while its relaxation time locus is strongly dependent on physical ageing time, the relaxation time loci being identical to the relaxation time for τ 213 K curve in Fig. 2 evaluated at same physical ageing times as above. Furthermore, relaxation spectrum at each physical ageing time shows a broad, continuous spectrum enabling a wider time scale found in real molecular networks to be taken into account as compared to classical models embodying discrete spectra, thus providing a plausible explanation of the successful fit of Mittag-Leffler function to measurement results while also explaining the broader times scales involved in the fractional free volume and relaxation time evolution results using Caputo fractional time derivatives in the differential evolution equations.

Numerical procedure and shear modulus evolution
The numerical procedure specified in "Appendix" applies a minimum time step of t = 10 −14 s, while the maximum relative tolerance error allowed is set to = 10 −8 . Furthermore, the exponential coefficient in the auxiliary function y T is selected τ T k . The minimum time step above is applied over several decades while using the numerical procedure in "Appendix". Then the time step is increased by a factor of 10, and the same procedure is repeated over a broader time span (more decades) while using the already determined values in the preceding evaluations. By this technique several decades (here from 10 −14 to 10 6 s) are effectively determined. The resulting evolution of fractional free volume f 213 K for temperature history with step to T 1 = 213 K is shown in Fig. 4 versus physical ageing time t ph together with exponential auxiliary function y 213 K and numerically determined function g applied in the numerical procedure. Clearly, the exponential auxiliary function satisfies the same initial value condition as f 213 K while lim t ph →∞ y 213 K = 0, whereas numerically determined function satisfies zero initial value condition g(0) = 0 and lim t ph →∞ g = f 213 K , that is, approaching sought function f 213 K for increasing physical ageing. In passing, it may be noted that the time span in Fig. 4 where auxiliary function y 213 K and function g interchange is relatively short, forcing the time step to be sufficiently short in order to get numerically stable solution within that time span. A possible improvement of the numerical procedure is to apply a low α-order Mittag-Leffler function instead of an exponential function to allow for a smoother and broader time span for interchange of function g and y 213 K , thus resulting in a more stable numerical procedure allowing for a longer time step.
The resulting evolution of shear modulus magnitude |μ| and loss factor μ/ μ, where denotes real part, for temperature history with step to T 1 = 213 K is shown in Fig. 5 versus the whole audible frequency range from 20 to 20 000 Hz at physical ageing time 0 + , 10 −9 , 10 −4 , 10 1 and 10 6 s. Clearly, the shear modulus for two shortest physical ageing times 0 + and 10 −9 s is mainly within the rubber region, displaying a low-level magnitude plateau in the low frequency range with a gradual increase above approximately 10 3 Hz, the rise being slightly more pronounced for the 10 −9 s curve. The loss factor increases with frequency in the whole audible frequency range, the curves mainly being vertically separated; the 10 −9 s curve shows the highest loss factor, while the 0 + s curve the lowest. The high-frequency fragment of the 10 −9 and 0 + s curves also enters  Fig. 4 Evolution of fractional free volume f 213 K versus physical ageing time t ph starting at equilibrium fractional free volume corresponding to temperature +25 • C and ending at equilibrium value corresponding to −60 • C together with exponential auxiliary function y 213 K and function g, the latter satisfying zero initial value condition

Shear modulus factor [-]
ρ-factor 1 (no factor) f-factor total-factor T-factor Fig. 6 Evolution of shear modulus ρ-factor ( and total factor (ρ × f × T -factor) versus physical ageing time t ph , associated with the fractional free volume f 213 K evolution starting at equilibrium fractional free volume corresponding to temperature +25 • C and ending at equilibrium value at −60 • C. The reference value 1 (no factor) is also shown the transition zone, whereas 10 −4 s curve lies almost entirely at the high-frequency end of that zone, the latter displaying a steep magnitude increase of more than 1 decade, while loss factor is showing a maximum of more than 100% at 20 Hz. Finally, the shear modulus for two longest physical ageing times 10 1 and 10 6 s is entirely within the glassy region, displaying a high-level magnitude plateau, almost overlapping each other, with only a minor increase with frequency, while the loss factor curves mainly being vertically separated show a decrease with frequency, down to surprisingly low value of about 0.005% at 20 000 Hz. In conclusion, the shear modulus undergoes a large alteration during the course of physical ageing, displaying a maximum magnitude and loss factor change of more than 2 and 4 decades, respectively, evaluated at a specific frequency.
There is a factor Tρ T /T 0 ∞ ρ T 0 in the shear modulus; see for example Eq. (6). This factor is multiplicatively split into ρ-factor (1−β rubber T ), f -factor (1+ f 213 K − β[A 2 + T ]) −1 and T -factor (1+ T /T 0 ), see for example Eq. (57), and is shown in Fig. 6 for temperature history with step to T 1 = 213 K versus physical ageing time t ph , together with the total factor (ρ × f × T -factor) and reference value 1 (no factor). Clearly, the entropy associated T -factor has the largest and f -factor the smallest influence on shear modulus magnitude. The evolution of f -factor starts approximately from (1 at t ph = 0 + , and ends approximately as unity (≈ 1) at t ph = 10 6 s, since lim t ph →∞ f 213 K = ∞ f 213 K . However, the total factor (ρ × f × T -factor) in the shear modulus has a smaller influence than the physical ageing shift factor X ph , as is clearly observed while comparing the results in Fig. 6 with those in Fig. 5.

Chemical ageing
The material constants applied for chemical ageing are E sci = 1.05 × 10 5 J mol −1 , E ref = 7.37 × 10 4 J mol −1 , τ sci0 = 8.08 × 10 8 s and τ ref0 = 1.00 × 10 8 s, experimentally determined by Johlitz et al. [21] for filled natural rubber using classical, differential evolution equations with ordinary derivatives, where the two latter material constants are here recalculated from other constants used in Johlitz et al. [21]. Moreover, applied fractional derivative material constants for the novel differential evolution equations developed in this paper read α sci = 0.5 and α ref = 0.5 together with maximum network link reformation possible q ref = 0.80 (80%).
Since the physical origin of chemical ageing-here in specific covering oxygen reaction kinetics with polymer network-is not exactly related to that of stress and strain response of the molecular network, at least not as much as the latter is related to corresponding physical ageing, the fractional order material constants for chemical ageing is therefore set to 1/2-a common fractional derivative order for mechanical rubber modelling. In passing, it is noted that the total evolution of equilibrium elastic modulus including both network link scission and reformation (52) is possible to reformulate [15] as for the selected fractional derivative order, where scaled complementary error function and error function The scaled complementary error function is a numerically stable function reducing possible under-and overflow during its numerical estimation. The resulting evolution of stress ratio dev[σ ] ageing / dev[σ ] no ageing = tr[σ ] ageing / tr[σ ] no ageing = μ ∞ (T ch , t ch )/μ ∞ 0 , see Eq. (58), at applied constant harmonic infinitesimal strain is shown in Fig. 7 versus a very broad chemical ageing time t ch range, starting from 10 2 s, which is of order of minutes (more precisely; 1 min 40 s), to 10 14 s, which is of order of tens of thousands years (more precisely, 31 700 years), at −60, −25, ±0, +25 and +60 • C for molecular network link scission only, reformation only and both simultaneously. Clearly, stress ratio is 1 at scission and no chemical ageing. Subsequently, it diminishes gradually with increased chemical ageing time; the higher the temperature, the faster is the decrease, where the −60 • C curve is displaying a vanishing chemical ageing even after 31,700 years, while the corresponding +60 • C curve shows already a reduction of stress ratio of about 4% after 3 h and drop of 28% after 11.5 days. All stress ratio curves will eventually diminish to zero although the decrease rate is strongly dependent on temperature. On the contrary, stress ratio is 0 at reformation and no chemical ageing. Subsequently, it grows gradually with increased chemical ageing time; the higher the temperature, the faster is the increase, where the −60 • C curve is displaying only a minor increase of 19% after 31,700 years, while the corresponding +60 • C curve shows already a growth of stress ratio of about 4% after 3 h and a rise to 30% after 11.5 days. All stress ratio curves will eventually increase to 80%, being the selected maximum network link reformation possible, although the growth rate is strongly dependent on temperature. Finally, the stress ratio at simultaneous scission and reformation is 1 at no chemical ageing. However, the stress ratio curves display a strong temperature dependence for increased chemical ageing; the −60 and −25 • C curves are laying above 1, +60 • C curve is almost entirely below 1, while ±0 and +25 • C curves are fluctuating around 1 within the considered chemical time range. In particular, the +25 • C curve displays a maximum stress ration at about 10 8 s (3.17 years) chemical ageing time due to previously increased shear modulus magnitude resulting from a prevailing molecular network link reformation before the link scission process starts to dominate. All stress ratio curves will eventually approach 0.80 although the chemical ageing time required is strongly dependent on temperature. The same stress ratio as above versus chemical ageing time t ch , at applied harmonic infinitesimal strain, is shown in Fig. 8 at temperature +25 • C, while using fractional time derivatives order of 0.25, 0.50, 0.75 and 1.00 in the inner differential evolution equation for molecular network link scission only, reformation only and both simultaneously. Apparently, the chemical ageing behaviour is largely similar for the studied derivative orders. However, the time band displaying the largest chemical ageing variation is strongly dependent on derivative order; the ordinary, integer derivative order (1.00) shows the shortest, while fractional derivative order of 0.25 the largest time band, as is clearly visible in the separate scission and reformation curves. In conclusion, the developed fractional differential equations for chemical ageing allow for a wider range of time scales most likely found in real oxygen reaction kinetics with polymer network.
Finally, the evolution of shear modulus magnitude |μ| and loss factor μ/ μ versus the whole audible frequency range from 20 to 20 000 Hz is shown in Fig. 9 at temperature +25 • C and chemical ageing time 0, 10 6 , 10 8 , 10 10 and 10 12 s, for simultaneous molecular network link scission and reformation. Clearly, shear modulus is mainly within the rubber region displaying a low-level magnitude plateau and low-level loss factor in the low frequency range with a gradual increase above approximately 10 3 Hz. The magnitude curve is vertically shifted at increased chemical ageing; the 10 6 and 10 8 s curves are shifted up and 10 10 and 10 12 s curves are shifted down with respect to no ageing curve. In particular, the 10 8 s curve displays the maximum shear modulus magnitude due to increased shear modulus magnitude resulting from a prevailing molecular network link reformation before the link scission process starts to dominate, as is clearly seen in Fig. 7. Finally, the loss factor curve is independent on chemical ageing due to the chemical ageing model applied, operating equally on elastic μ and loss part μ of the shear modulusμ.

Conclusion
The developed constitutive models for rubber including physical and chemical ageing while using fractional differential evolution equations provide a richer model structure, apply a small material parameter number and Evolution of shear modulus magnitude |μ| and loss factor μ/ μ versus the whole audible frequency range from 20 to 20 000 Hz, at chemical ageing time 0, 10 6 , 10 8 , 10 10 and 10 12 s, for simultaneous molecular network link scission and reformation at temperature +25 • C display a wide spectrum of time scales found in real polymer network. The constitutive models are possible to use in various applications, for example, to investigate physical ageing under specific temperature histories, relative contribution from scission and reformation at chemical ageing, the temperature influence on chemical ageing and to investigate how physical and chemical ageing influence dynamic stiffness of rubber vibration isolators within the audible frequency range. The latter is considered in part 2 of this paper [23].

Appendix: Numerical procedure to compute f T
The relaxation time at temperature T reads using (25) and (37), and (21) becomes which is a nonlinear fractional differential evolution equation for the fractional free volume f T , with a generally nonzero initial condition, f T (0) = 0 f T ≥ 0, while requiring a numerical procedure for its solution. To this end, a special case of a general procedure to solve nonlinear fractional differential equations [40] is applied: 1 Initiating (a) Transform the differential equation (62) into a differential equation with zero initial value condition: i Select an exponential auxiliary function y T = 0 f T exp (λ t) that satisfies the same initial condition as the original function f T and where λ < 0 is an arbitrary constant to be à priori selected ii Let f T def = g + y T , where function g satisfies zero initial value condition; g(0) = 0 iii The differential equation (62) transforms into where Y (t ph ) = D α [y T ] [17] and the generalized, two-parameter Mittag-Leffler function , that is, the right-hand side of equation (63) (c) Set time step-size t = t ph /n resulting in a discrete time set 0 ≤ k t ≤ t ph , where n is total number of steps and 0 ≤ k ≤ n (d) Let  g using Grünwald-Letnikov operator [39] for fractional Caputo derivative: where fractional binomial coefficient recursively reads Ω 0 = 1 and Ω m = 1 − 1+α