Dynamic stiffness of chemically and physically ageing rubber vibration isolators in the audible frequency range: Part 2—waveguide solution

The dynamic stiffness of a chemically and physically ageing rubber vibration isolator in the audible frequency range is modelled as a function of ageing temperature, ageing time, actual temperature, time, frequency and isolator dimension. In particular, the dynamic stiffness for an axially symmetric, homogeneously aged rubber vibration isolator is derived by waveguides where the eigenmodes given by the dispersion relation for an infinite cylinder satisfying traction free radial surface boundary condition are matched to satisfy the displacement boundary conditions at the lateral surface ends of the finite rubber cylinder. The constitutive equations are derived in a companion paper (Part 1). The dynamic stiffness is calculated over the whole audible frequency range 20–20,000 Hz at several physical ageing times for a temperature history starting at thermodynamic equilibrium at $$+25\,^{\circ }\hbox {C}$$+25∘C and exposed by a sudden temperature step down to $$-60\,^{\circ }\hbox {C}$$-60∘C and at several chemical ageing times at temperature $$+25\,^{\circ }\hbox {C}$$+25∘C with simultaneous molecular network scission and reformation. The dynamic stiffness results are displaying a strong frequency dependence at a short physical ageing time, showing stiffness magnitude peaks and troughs, and a strong physical ageing time dependence, showing a large stiffness magnitude increase with the increased physical ageing time, while the peaks and troughs are smoothed out. Likewise, stiffness magnitude peaks and troughs are frequency-shifted with increased chemical ageing time. The developed model is possible to apply for dynamic stiffness prediction of rubber vibration isolator over a broad audible frequency range under realistic environmental condition of chemical ageing, mainly attributed to oxygen exposure from outside and of physical ageing, primarily perceived at low-temperature steps.


Introduction
Vibrations from sources such as engines, machines and fans frequently cause unwanted emission of noise while transmitting into receiving structure. Noise is generated by structure-borne sound, that is, vibrations in the audible frequency range that connects to surrounding air at surfaces receiving the transmitted vibrations from the source through the structure. The ability to reduce noise emission is a key market discriminator for presumptive consumers of devices containing any vibration source, and a simple measure is to decouple the source from receiving structure by mounting it upon vibration isolators. Vibration isolators are frequently made of rubber due to its high elasticity and moderate material damping and the possibility of designing rubber vibration Communicated by Andreas Öchsner.
L. Kari (B) The Marcus Wallenberg Laboratory for Sound and Vibration Research (MWL), KTH Royal Institute of Technology, 100 44 Stockholm, Sweden E-mail: leifkari@kth.se isolators with desirable properties along with their low need of maintenance. The basic principle of vibration isolation is to provide mechanical mismatch between source and receiving structure, the mismatch acting as a mechanical mirror, reflecting vibrations backwards to the source while minimizing transmitted vibrations and thus the noise emitted. As a general rule of thumb, sources and receiving structures are hard, thus requiring a soft vibration isolator for appropriate vibration transmission reduction; the softer the vibration isolator is, the better is the attenuation. However, vibration isolators also offer stability and convey the static load of source resulting in a strongly constrained multi-dimensional design parameter space in the attenuation optimization procedure, thus clearly excluding a vanishing vibration isolator stiffness solution. Furthermore, source and receiving structure undergo resonances and anti-resonances within the audible frequency range resulting in a break down of the rule of thumb while showing both soft and stiff behaviour depending on the studied frequency. Finally, vibration isolators display internal resonances and anti-resonances within the same audible frequency range also resulting in alternating soft and stiff behaviour depending on frequency [11,27]. As a result, the mechanical properties of source, vibration isolator and receiving structure need to be properly determined for appropriate vibration transmission attenuation within the audible frequency range where generally the dynamic stiffness modelling of vibration isolators provides the greatest challenge.
There are various techniques to model dynamic stiffness of vibration isolators; a simple, straightforward method is to apply lumped parameter models with Maxwell-Voigt elements [32,33] and fractional damping elements [1,2,8], the latter shown to fit rubber material experiments with a minimum number of material parameters. In order to include internal resonances and anti-resonances, a multi-dimensional mass-spring system [28] and classical, essentially one-dimensional waveguides [20][21][22] may be applied. However, the wavelength in rubber is generally short within the audible frequency range resulting in a break down of the assumption on which the classical, essentially one-dimensional waveguides methods rely on, thus requiring more elaborated methods such as waveguide methods derived from the full-wave equations that are also valid for short wavelengths [7,15,16,19,34,35] and are essentially based on the possible eigenmodes found in the Pochhammer-Chree solution for waves in infinite solid cylinders [6,26]. For more elaborated geometries, boundary elements method [9] and nonlinear finite element methods [17,29,30] are applicable, the latter also taking into account nonlinear viscoelasticity. However, none of the above methods includes physical and chemical ageing, the former being a reversible process including such as free volume alterations, trapping and freeing of polymer chain ends and other polymer chain reorganizations [4,12,13,24], while the latter being an irreversible process mainly attributed to oxygen reaction with polymer network either damaging the network by scission or reformation of new polymer links [13,14]. Rubber vibration isolators undergo ageing while frequently being exposed for harsh environmental conditions under a long time including oxygen, water, ozone, oil, contaminants and radioactive radiation exposures, usually at various temperatures, resulting in a dynamic stiffness alteration from that predicted by methods assuming a virgin state, disregarding any ageing.
In this paper, dynamic stiffness within the whole audible frequency range 20 to 20,000 Hz for an axially symmetric, homogeneously aged rubber vibration isolator is derived by waveguides while using constitutive equations, derived in a companion paper [18], including physical and chemical ageing. The resulting model is suitable for stiffness prediction of rubber vibration isolators over a broad audible frequency range under realistic environmental condition of chemical ageing, mainly attributed to oxygen exposure from outside, and of physical ageing, primarily perceived at low-temperature steps.

Vibration isolator
The cylindrical vibration isolator in Fig. 1 of radius r 0 and length l 0 , consisting of isotropic, homogeneously aged rubber of density 0 at reference temperature T 0 , is axially excited a displacement u excitation with a force f excitation at one end while being blocked u blocked = 0 at the other end generating a blocking force f blocked . The radial displacement is enforced to be zero at both ends, in practical achieved by gluing a circular, thin metal plate on the excitation end, while the other end of rubber cylinder is glued to the hard, brown blocking end plate shown in Fig. 1. The dynamic driving point stiffness k D = f excitation / u excitation and dynamic transfer stiffness k T =f blocked /ũ excitation , where temporal Fourier transformation [ · ] = ∞ −∞ [ · ] exp(−iωt) d t, t is time, i is imaginary unit and ω is angular frequency.

Waveguides
The constitutive equations for rubber are derived in a companion paper [18], including physical and chemical ageing, giving bulk modulus and shear moduluŝ where actual temperature is T , temperature difference T = T − T 0 , physical ageing time is t ph , and chemical ageing time is t ch . Moreover, network link scission relaxation time network link scission relaxation time is τ sci 0 at reference temperature T 0 , chemical ageing temperature is T ch , network link scission activation energy per mole is E sci , universal molar gas constant R = 8.314 J/mol K, network link reformation relaxation time Gamma function is Γ (z) = ∞ 0 s z−1 e −s d s, network link scission fractional order is α sci , network link reformation fractional order is α ref , maximum network link reformation possible is q ref , thermal volume expansion coefficient of rubber at the rubbery and crystalline state is β rubber and β crystal , respectively, fractional free volume thermal expansion coefficient β ≈ β rubber − β crystal , nearly incompressibility material constant is a 1, material constants in shear are μ ν and α, respectively, non-dimensional relaxation intensity is 1 and equilibrium elastic modulus is μ ∞ 0 at reference temperature T 0 and at no chemical ageing. William-Landel-Ferry (WLF) shift function using the method of reduced variables, A 1 and A 2 are material constants, physical ageing shift factor evolution of actual fractional free volume f T at temperature T is given by fractional differential equation fractional free volume at the thermodynamical equilibrium state is and finally, Caputo fractional derivative of f T and order α [5] with respect to physical ageing time t ph is The density relation where actual density is T and equilibrium density is ∞ T 0 at temperature T , results in approximate expressions for the actual vibration isolator length and radius while assuming an homogeneous thermal expansion, where the relative equilibrium volume expansion at T is β rubber T , equilibrium fractional free volume at T 0 and T are β A 2 and β A 2 + β T , respectively. The waveguide method is based on mode-matching technique [7,15,16,19,34,35], here extended to include physical and chemical ageing. To this end, the axial dependence is separated, the eigenvalues and eigenmodes of the cylinder cross section are calculated, where the total field is obtained by eigenmode superposition, and finally they are matched with the cylinder end boundary conditions given in Sect. 2 and Fig. 1. The reader is referred to Kari [15,16] for details regarding the complete derivation of the waveguide guide solution for axial waves in solid rubber cylinders and their resulting dynamic stiffness solution, however restricted to nonageing material at reference temperature T 0 . The transcendental, dispersion relation for axially symmetric and non-torsional waves in an infinite rubber cylinder, extended to include physical and chemical ageing, is while satisfies the traction free boundary condition on the radial surface in Fig. 1, where k z = k S , k z = k L , Onoe function of first kind and first order is ϑ(z) = zJ 0 (z)/J 1 (z) [23,25], Bessel function of first kind and order n is J n , longitudinal wave number and shear wave number are and respectively, and actual density The solution to Eq. (14) is eigenvalue k z , also known as axial wave number. The solution is calculated by a modified Newton-Raphson method with initial values given by a winding integral method based on argument principle where search domain is spit into branch-cut absent sub-domains, resulting in eigenvalues where n = 1, 2, . . . , M, total number of eigenmodes is M, eigenvalues are ordered in increasing imaginary part value as | k z,1 | ≤ | k z,2 | ≤ · · · ≤ | k z,M | and imaginary part is denoted . The search domain is restricted to k z,n ≥ 0 as −k z,n is also a solution to Eq. (14), where real part is denoted . The total field is obtained by eigenmode superposition using a sufficient number of eigenvalues M, where rubber cylinder end displacement boundary conditions given in Sect. 2 and Fig. 1 are satisfied by modematching, and where the resulting dynamic stiffness solution is readily calculated from the resulting total eigenmode field [15,16].

Results and discussion
The computer code is written in Lahey Fortran 95/90/77 with support for Fortran 2003 and 2008 standards and with all calculations performed in quadruple precision using a Windows 7 Professional operative system implemented in a HP EliteBook 840 G3 laptop with Intel Core TM i7-6500 CPU @ 2.50 GHz and 16.0 GB of RAM. Graphically, the results are presented by means of Matlab .

Eigenvalue and dynamic stiffness evolution at physical ageing of vibration isolator
The applied temperature history is in order to study eigenvalue and dynamic stiffness evolution at physical ageing, where h is Heaviside step function, while using the numerical procedure specified in Appendix of companion paper [18] for evaluation of Eq. (8) with minimum step of t = 10 −14 s, maximum relative tolerance error allowed set to = 10 −8 and exponential coefficient in the auxiliary function The resulting evolution of the complete set of the 64 first axial wave numbers (eigenvalues) k z,n applied in the dynamic stiffness calculation, where n = 1, 2, . . . , 64, is shown in Fig. 2 displaying real and imaginary part of axial wave number k z,n and k z,n , respectively, versus the whole audible frequency range from 20 to 20,000 Hz at physical ageing time 0 + , 10 −4 and 10 6 s. The complete 3D spectrum is shown at top of Fig. 2, while the corresponding slice at frequency 10,000 Hz is shown to left in Fig. 2, together with shear and longitudinal wave number k T and k L , respectively. Note that also −k z,n , n = 1, 2, . . . , 64, are used in the stiffness calculation, thus making the total eigenmode number to be 128. In passing, it was noted during the simulations that the difference between this stiffness solution and a solution with more eigenmodes, for example totally 256 eigenmodes, is negligible, thus verifying the solution accuracy and stiffness convergence regarding sufficient number of eigenmodes. The reader is referred to Kari [16] regarding details over convergence properties of the waveguide method applied to a similar vibration isolator although disregarding any ageing. Apparently, the axial wave number spectrum is seemingly complex, displaying eigenvalue solutions to Eq. (14) in all four quadrants of the complex plane while also counting −k z axial wave numbers. The number of propagating eigenmodes at physical ageing time 0 + s grows rapidly, being 1 + 1 at 20 Hz, 4 + 4 at 10,000 Hz and 6 + 6 at 20,000 Hz, where x and y in x + y refer to number of propagating eigenmodes with real part of eigenvalue being positive and negative, respectively. Formally, propagating eigenmodes do not exist for rubber and other viscoelastic materials. However, propagating eigenmodes do formally exist for purely elastic materials with no material damping. Whether an eigenmode is propagating or not is therefore determined from the cut-on frequency for its elastic counterpart, derived by means of the formal substitutionμ ← μ in the dispersion relation (14). On the contrary, the number of propagating eigenmodes at physical ageing time 10 −4 and 10 6 s is constant throughout the considered audible frequency range, that is, being 1 + 1 at 20 Hz, 1 + 1 at 10,000 Hz and 1 + 1 at 20,000 Hz. Not surprisingly, the magnitude of shear modulus grows rapidly from physical ageing time 0 + s to 10 −4 and 10 6 s, as is readily shown in Fig. 5 of Kari [18], thus resulting in a rapid increase in the cut-on frequency of the higher-order eigenmodes, being above the audible frequency range at physical ageing time 10 −4 and 10 6 s. However, there are always 1 + 1 continuously propagating eigenmodes with their two eigenvalues starting from origin k z = 0 at vanishing frequency, showing continuously increasing real part with increased frequency for the first eigenvalue (shown in Fig. 2) and increasing negative real part with increased frequency for the second eigenvalue (not shown in Fig. 2), respectively, both being close to the plane k z = 0 throughout the considered audible frequency range. The shear wave number follows closely the eigenvalue of the first continuously propagating eigenmode for lower frequencies and then close to the eigenvalue of the second propagating eigenmode for higher frequencies. In passing, it is noted that the first continuously propagating eigenmode at physical ageing time 0 + s approaches the Rayleigh wave, while the second propagating eigenmode (lowest higher-order eigenmode) approaches the shear wave as their short wavelength limits when frequency goes to infinity [10]. However, the longitudinal wave number is close to the origin; not surprisingly, the magnitude of bulk modulus is significantly larger than the corresponding magnitude of shear modulus due to the nearly incompressible material model, resulting in a significantly smaller magnitude of the longitudinal wave number than the corresponding magnitude of the shear wave number. This is readily visible for ageing time 0 + s, while the shear and longitudinal wave numbers are more clustered together close to the origin k z = 0 at physical ageing time 10 −4 and 10 6 s due to the rapid increase in shear modulus magnitude.   Evolution of real and imaginary part of axial wave number k z,n and k z,n , respectively, where n = 1, 2, . . . , 64, versus the whole audible frequency range from 20 to 20,000 Hz, at an applied temperature history undergoing a sudden temperature drop from equilibrium at +25 to −60 • C and evaluated at physical ageing time 0 + , 10 −4 and 10 6 s. The complete 3D spectrum is shown above, while the corresponding slice at frequency 10,000 Hz is shown below, together with shear and longitudinal wave number k T and k L , respectively Moreover, a propagating eigenmode for purely elastic material has a purely real axial wave number, while the non-propagating eigenmodes have non-vanishing imaginary parts. The latter wave numbers are possible to split into complex and purely imaginary axial wave numbers; both eigenmodes represent near-fields, but those for complex axial wave numbers also involve spatially oscillating factors. Furthermore, complex axial wave numbers occur in fours, one in each of the four quadrants, obeying the additional symmetry k * z ↔ k z , in addition to the usual point symmetry with respect to origin (k z = 0), that is, −k z ↔ k z , where * denotes complex conjugate, whereas the real and imaginary axial wave numbers occur only in pairs (−k z ↔ k z ). Since energy dissipation in purely elastic material is impossible, a complex wave number must be augmented with its dual from other quadrant. Extended to the material under study focus, that is, viscoelastic rubber, the purely imaginary and complex axial wave numbers are clustered together into smoothly banana-shaped axial wave number pattern in the complex plane; the larger the imaginary part is, the more frequency and physical ageing time independent is the eigenvalue. Not surprisingly, the complex eigenvalue solution to the dispersion relation (14) depends approximately only on radius of cylinder r and large eigenvalue order p 1 : p ∈ N as k z ≈ [± log e (4 pπ) ± 2 pπi]/2 r according to Zemanek [31]. The complex eigenvalues are approximately overlapping at large values of p for physical ageing time 10 −4 and 10 6 s; the small discrepancies are partly due to different radius at those physical ageing times alongside with other approximations involved. However, the required value of p to satisfyingly approximate the complex eigenvalues is very large for physical ageing time 0 + s; merely valid only for complex eigenvalues well outside the important eigenvalues ±k z,n needed for sufficiently accurate stiffness determination of the considered vibration isolator, where n = 1, 2, . . . , 64. This is clearly seen in Fig. 2 where the complex eigenvalue pattern for physical ageing time 0 + s is positioned well from the corresponding complex eigenvalue pattern for physical ageing time 10 −4 and 10 6 s and for frequencies above approximately 5000 Hz. Not surprisingly, the large number of propagating and close to be propagating eigenmodes, for physical ageing time 0 + s, results in a somewhat jumbled pattern for the complex eigenvalues in comparison to physical ageing time 10 −4 and 10 6 s, where only the continuously propagating eigenmodes are propagating and the cut-on frequency for lowest higher-order eigenmode is very high about 50,000 Hz, thus far away from the audible frequency range. The resulting evolution of driving point stiffness magnitude |k D | and angle k D , where k D = |k D | exp(i k D ), is shown in Fig. 3 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, driving point stiffness magnitude is a plateau at the low-frequency range for physical ageing time 0 + and 10 −9 s and at an extended frequency range for physical ageing time 10 1 and 10 6 s; the magnitude ratio between the latter and former curves being of the same order as relaxation intensity . Not surprisingly, the former curve portions are within the rubber region, while the latter curve portions are within the glassy region. Furthermore, the wavelength of the propagating eigenmode within the corresponding frequency ranges is significantly longer than the length of the vibration isolator, thus making a lumped massless stiffness model suitable, essentially displaying a shear modulus proportional stiffness showing neither resonances nor anti-resonances. However, driving point stiffness magnitude is gradually decreasing from about 500 Hz for physical ageing time 0 + and 10 −9 s, the decrease being slightly more pronounced for 0 + s curve, while approaching a damped stiffness trough at about 1200 Hz, interpretable as a resonance. For higher frequencies, driving point stiffness magnitude is displaying damped stiffness peaks (anti-resonances) and troughs (resonances) in alternating order for physical ageing time 0 + and 10 −9 s. Likewise, the driving point stiffness magnitude for physical ageing time 10 1 and 10 6 s is displaying a deep trough at about 18,300 Hz, the drop being more than four decades. Not surprisingly, shear modulus loss factor is less than 0.01 % around that resonance frequency for physical ageing time 10 6 s, see Fig. 5 in Kari [18]. The driving point stiffness curves for physical ageing time 10 1 and 10 6 s are almost overlapping; the only practical difference is the depth of the stiffness trough and the low-frequency driving point stiffness angle where 10 6 s curve is displaying about a decade smaller angle as compared to the 10 1 s curve. However, the latter difference is not perceptible in Fig. 3 since both low-frequency driving point stiffness angles are very close to 0 • . The driving point stiffness magnitude for physical ageing time 10 −4 s is displaying a gradual increase in the low-frequency range while displaying a very damped stiffness trough at about 17,400 Hz. The latter 10 −4 s curve is mainly within the transition region as is clearly seen in the corresponding driving point stiffness angle curve for Hz, corresponding to a shear modulus loss factor of more than 100% while using a suitable lumped stiffness model. In passing, it is noted that the maximum shear modulus loss factor indeed is larger than 100% for the rubber material under study focus, see Fig. 5 in Kari [18]. Finally, the driving point stiffness angle is restricted to be within 0 and 180 • due to physical reasons. Obviously, this restriction is clearly satisfied for the studied vibration isolator. The resulting evolution of transfer stiffness magnitude |k T | and angle k T , where k T = |k T | exp(i k T ), is shown in Fig. 4 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. Essentially, transfer stiffness is displaying similar behaviour as the corresponding driving point stiffness. However, stiffness peaks and troughs are generally not located at the same frequencies nor is peak and trough order necessarily alternating. The latter statement is a result from dropping the stiffness angle restriction to be within 0 and 180 • , a restriction not physically necessary for transfer stiffness angle and obviously also not for the studied vibration isolator. In passing, it is noted that the general transfer stiffness frequency dependence is coherent with previous simulations [19] and measurements [15].
In conclusion, driving point and transfer stiffness are undergoing a fast physical ageing for the studied temperature history, displaying a large stiffness magnitude increase in order in the low-frequency range, while the high-frequency stiffness, initially showing stiffness magnitude peaks and troughs, is transformed into a high stiffness magnitude while displaying a sharp and deep driving point stiffness magnitude trough at about 18,300 Hz. The physical stiffness ageing is merely completed within 10 s, although the physical ageing process continues for a longer time, showing almost no stiffness change after 10 6 s. Similarly, eigenvalue spectrum is undergoing a fast physical ageing for the studied temperature history, in particular displaying a fast alteration of higher-order propagating eigenmodes into non-propagating eigenmodes with complex eigenvalues in addition to moving the complex eigenvalue pattern into a smooth configuration which is less frequency dependent. Likewise, the physical eigenvalue ageing is merely completed within 10 s, although the physical ageing process continues for a longer time, showing almost no eigenvalue change after 10 6 s.

Eigenvalue and dynamic stiffness evolution at chemical ageing of vibration isolator
In order to study eigenvalue and dynamic stiffness evolution at chemical ageing, the vibration isolator is exposed to a constant temperature of 25 • C while undergoing simultaneous and homogeneous molecular network scission and reformation.
The resulting evolution of all 64 first axial wave numbers (eigenvalues) k z,n applied in the dynamic stiffness calculation, where n = 1, 2, . . . , 64, is shown in Fig. 5 displaying real and imaginary part of axial wave number k z,n and k z,n , respectively, versus the whole audible frequency range from 20 to 20,000 Hz at chemical ageing time 0, 10 8  Evolution of real and imaginary part of axial wave number k z,n and k z,n , respectively, where n = 1, 2, . . . , 64, versus the whole audible frequency range from 20 to 20,000 Hz, at temperature 25 • C undergoing simultaneous molecular network scission and reformation and evaluated at chemical ageing time 0, 10 8 and 10 12 s. The complete 3D spectrum is shown above, while the corresponding slice at frequency 10,000 Hz is shown below, together with shear and longitudinal wave number k T and k L , respectively frequency 10,000 Hz is shown to left in Fig. 5, together with shear and longitudinal wave number k T and k L , respectively. Note again that also −k z,n , n = 1, 2, . . . , 64, are used in the stiffness calculation, thus making the total eigenmode number to be 128. In passing, it was again noted during the simulations that the difference between this stiffness solution and a solution with more eigenmodes, for example totally 256 eigenmodes, is negligible, thus verifying the solution accuracy and stiffness convergence regarding sufficient number of eigenmodes. Apparently, the number of propagating eigenmodes grows rapidly, being 1 + 1 at 20Hz and 3 + 3 at 10,000 Hz at chemical ageing time 0, 10 8 and 10 12 s, 5 + 5 at 20,000 Hz at chemical ageing time 0 and 10 8 s and, finally, 6 + 6 at 20,000 Hz at chemical ageing time 10 12 s, where again x and y in x + y refer to number of propagating eigenmodes with real part of eigenvalue being positive and negative, respectively. Real parts of eigenvalue for the first continuously propagating eigenmode and of shear wave number are larger for chemical ageing time 10 12 s than those corresponding to chemical ageing time 0 s, which in turn are larger than those corresponding to chemical ageing time 10 8 s. Not surprisingly, shear modulus magnitude for chemical ageing time 10 12 s is smaller than that corresponding to chemical ageing time 0 s, which in turn is smaller than that corresponding to chemical ageing time 10 8 s at a specific frequency, see Fig. 9 in Kari [18]. Furthermore, the complex eigenvalue pattern is in general more jumbled and non-symmetric with respect to plane k z = 0 as compared to the complex wave pattern at physical ageing just studied. Not surprisingly, as the number of propagating eigenmodes are large throughout the studied chemical ageing time range, thus chiefly influencing the complex eigenvalue pattern. Moreover, longitudinal wave number for all chemical ageing time 0, 10 8 and 10 12 s is clustered together and close to origin (k z = 0), mainly due to the nearly incompressible material model applied in this study, resulting in a significantly larger bulk modulus magnitude compared to shear modulus magnitude and thus significantly smaller longitudinal wave number magnitude than shear wave number magnitude. Finally, the resulting evolution of driving point stiffness magnitude |k D | and angle k D is shown in Fig. 6 and transfer stiffness magnitude |k T | and angle k T in Fig. 7, 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. The low-frequency stiffness magnitude is a plateau while displaying troughs (resonances) and peaks (anti-resonances) at higher frequencies. The lowfrequency stiffness magnitude for chemical ageing time 10 8 s is larger than that corresponding to chemical ageing time 10 6 s, which in turn is larger than that corresponding to chemical ageing time 0 s, while the lowfrequency stiffness magnitude for chemical ageing time 10 12 s is smaller than that corresponding to chemical ageing time 10 10 s, which in turn is smaller than that corresponding to chemical ageing time 0 s. Likewise, the stiffness magnitude trough and peak frequencies for chemical ageing time 10 8 s are larger than those corresponding to chemical ageing time 10 6 s, which in turn are larger than those corresponding to chemical ageing time 0 s, while the trough and peak frequencies for chemical ageing time 10 12 s are smaller than those corresponding to chemical ageing time 10 10 s, which in turn are smaller than those corresponding to chemical ageing time 0 s. This is not surprising, as shear modulus magnitude for chemical ageing time 10 8 s is larger than that corresponding to chemical ageing time 10 6 s, which in turn is larger than that corresponding to chemical ageing time 0 s at a specific frequency, and as shear modulus magnitude for chemical ageing time 10 12 s is smaller than that corresponding to chemical ageing time 10 10 s, which in turn is smaller than that corresponding to chemical ageing time 0 s, at a specific frequency, see Fig. 9 in Kari [18]. The stiffness magnitude trough and peak frequency shifts are readily seen also in the stiffness angle subplots where the stiffness magnitude troughs and peaks are manifested by +180 • and −180 • angle shifts, respectively, although smeared out due to material damping and strong resonance and anti-resonance overlapping. It is finally noted that the driving point stiffness angle is within the physically allowed range 0 to 180 • regardless of frequency and chemical ageing time, while the corresponding angle for transfer stiffness is frequently outside that range; not surprisingly, no physical restriction range exists for the latter angle.

Conclusion
The developed waveguide model for dynamic stiffness of physically and chemically ageing rubber vibration isolator within the audible frequency range is possible to apply in various applications. An example is to investigate the stiffness alteration due to physical ageing under specific temperature histories. A second example is to investigate the relative contribution to stiffness change from scission and reformation at chemical ageing. A third example is to investigate the temperature influence on stiffness due to chemical ageing and, finally, to investigate stiffness change due to simultaneous physical and chemically ageing including simultaneous scission and reformation, the latter case being a more realistic environmental condition for vibration isolators in practice compared to stiffness modelling assuming an eternally lasting virgin vibration isolator. An interesting continuation of the work performed is to extend the model with oxygen diffusion resulting in a chemically heterogeneous ageing of the vibration isolator, being an even more realistic case for vibration isolators in practice.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.