Using Waveguides to Model the Dynamic Stiffness of Pre-Compressed Natural Rubber Vibration Isolators

Abstract

A waveguide model for a pre-compressed cylindrical natural rubber vibration isolator is developed within a wide frequency range—20 to 2000 Hz—and for a wide pre-compression domain—from vanishing to the maximum in service, that is 20%. The problems of simultaneously modeling the pre-compression and frequency dependence are solved by applying a transformation of the pre-compressed isolator into a globally equivalent linearized, homogeneous, and isotropic form, thereby reducing the original, mathematically arduous, and complex problem into a vastly simpler assignment while using a straightforward waveguide approach to satisfy the boundary conditions by mode-matching. A fractional standard linear solid is applied as the visco-elastic natural rubber model while using a Mittag–Leffler function as the stress relaxation function. The dynamic stiffness is found to depend strongly on the frequency and pre-compression. The former is resulting in resonance phenomena such as peaks and troughs, while the latter exhibits a low-frequency magnitude stiffness increase in addition to peak and trough shifts with increased pre-compressions. Good agreement with nonlinear finite element results is obtained for the considered frequency and pre-compression range in contrast to the results of standard waveguide approaches.

1. Introduction

Natural rubber is an interesting polymeric material. It exhibits frequency- and pre-strain-dependent mechanical properties, including its shear storage and loss modula, and an immense dissimilarity between the transversal and longitudinal wave velocities [1]. The latter is due to its nearly incompressibility and displays similarities to that of water: comparatively soft in shear but hard in compression—in addition to a density close to 1000 kg/m${}^3$. Natural rubber has further excellent elastic and mechanical properties including abrasion resistance, tensile strength, and elongation resistance. Natural rubber vibration isolators are frequently used to insulate receiving structures from a source thanks to their efficiency in attenuating structural vibrations. To ensure an optimized utilization, the dynamic stiffness of natural rubber vibration isolators needs to be accurately determined. While steel springs have a well-established mechanical behavior, natural rubber vibration isolators still include features that require more modeling work.

Simplified models such as the classical long rod example are extensively applied to vibration isolators [2,3,4,5], despite some major limitations. The long rod model in which plane cross sections perpendicular to the axial co-ordinate axis are assumed to remain plane and perpendicular to the axial co-ordinate axis and where the stress is assumed to be uniaxial and uniform while radial expansions and contractions are neglected fails to give satisfactory predictions [6]. In particular, the lack of success is frequent for a vastly common vibration isolator—namely, the cylindrical rubber isolator with a large diameter-to-length ratio. The increased influence of higher-order modes, which are necessary to satisfy the boundary conditions, and the increased influence of dispersion explain the lack of correctness of the approximate solution. The dynamic stiffness results of a refined waveguide model [7,8] accounting for those effects are shown to diverge substantially from all the investigated, traditional simple models.

However, the above approaches do not take into account the influence of pre-compression, which is from a practical point of view, an inescapable parameter since most vibration isolators carry the static load of the source or isolated system they are attached to. These effects involve geometrical and material non-linearities, rendering the modeling task substantially more complex. Indeed, the original, free outer surface of the isolator is transformed into a deformed shape, clearly non-identical to the original surface, by the application of the static load—a geometrical non-linearity. The geometrical non-linearity results in a nonlinearity in the axial stiffness due to geometrical changes. In addition, the response to a superimposed dynamic excitation displays a seemingly anisotropic and non-homogeneous feature in the transformed state—a material non-linearity. The material non-linearity results in a nonlinearity in the axial stiffness due to nonlinearity in the stress–strain relation for the material itself. A way to account for the pre-compression dependence is by adopting a nonlinear approach based on single integral viscoelastic constitutive theory. The mechanical response of two natural rubber compounds were examined [9,10] under steady-state torsion oscillations combined with a large static extension while additionally assuming a time–strain separability of the relaxation tensor and displayed good agreement with experiments. Using a similar approach, the dynamic stiffness of a cylindrical vibration isolator was derived [11] for different pre-compressions within an audible frequency range obtaining excellent agreements with a series of experiments. To describe the audible short term response of the rubber, a fractional derivative type relaxation function was successfully adopted [7,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74], significantly reducing the number of required material parameters. Visco-elastic models with fractional derivatives have sound bases in polymer molecular theory [75]. The dynamic shear modulus was derived by using the principle of virtual work and by linearizing the constitutive relation at a given prestrain [76]. The model was applied to rubber bushings under simple shear and comparisons to experiments for different pre-compressions displayed good agreement although the frequency range was narrow: 0 to $100$Hz [76].

Although good results were obtained using the aforementioned models, the intrinsic non-linearity of the theory of pre-strained isolators, involving a series of complex mathematical difficulties, renders its application arduous. Furthermore, the use of non-linear finite element methods makes the results difficult to interpret as no closed-form solutions are obtained.

The present paper proposes a waveguide approach for the dynamic stiffness modeling of natural rubber vibration isolators accounting for both the pre-compression and frequency dependencies. This original, mathematically cumbersome task is simplified by transforming the complex shape of the pre-compressed natural rubber isolator into a globally equivalent, homogeneous, and isotropic form, combined with a straightforward waveguide approach to satisfy the boundary conditions by mode matching. A fractional standard linear solid is applied while using a Mittag–Leffler function as the stress relaxation function [14]. The proposed method is applied to a cylindrical rubber isolator, where the influence of pre-compression, cylinder length, and diameter are investigated. Comparisons with the results of non-linear finite element method are made.

2. Materials and Methods

2.1. Vibration Isolator

The studied natural rubber isolator in Figure 1 consists of a vulcanized natural rubber cylinder of height H and diameter $D=2R$ in its natural state, stress free and undeformed, firmly bonded to two parallel, identical metal plates of mass $m_{\mathrm{mp}}$.

2.2. Kinematics

Consider the vibration isolator, a simple body consisting of continuously distributed natural rubber material occupying a reference configuration $\mathcal{B}$ in its natural state, and refer to $\mathcal{S}$ as an ambient space in which the evolution of the body takes place, where $\mathcal{B}$ and $\mathcal{S}$$\subset {\mathbb{R}}^3$ are sufficiently smooth, oriented open manifolds. Consider an excitation source applied to the studied natural rubber isolator; the isolator undergoes a regular motion $\phi :\mathcal{B}\to \mathcal{S}$ from the reference configuration $\mathcal{B}$ to the current $\mathcal{S}$, which is decomposed into $\phi =\Delta \phi \circ {\phi }_0:\mathcal{B}\to \mathcal{S}\to \mathcal{S}$, where the excitation source’s mass load and other applied external static forces result in a finite motion ${\phi }_0:\mathcal{B}\to \mathcal{S}$, being the quasi-static natural rubber isolator pre-compression, while a superimposed infinitesimal excitation results in a motion $\Delta \phi :\mathcal{S}\to \mathcal{S}$. Consequently, the compression displacement boundary condition $\mathbf{u}=u{\mathbf{e}}_z$ at the upper bonded natural rubber surface in $\mathcal{S}$ is decomposable into $\mathbf{u}=(u_0-\Delta u\left(t\right)){\mathbf{e}}_z$, where $\Delta u\left(t\right)$ is a superimposed axial displacement, being the infinitesimal dynamic source excitation, and t is time.

2.3. Equivalent Transformation

Accounting for the nearly incompressibility of natural rubber, an approximate globally equivalent configuration is derived from the pre-compressed natural rubber isolator configuration $\mathcal{S}$, as depicted in Figure 2. The unknown deformed shape obtained after the pre-compression is transformed into a straight cylinder of identical volume. Consider the pre-compression ${\phi }_0$ in Figure 1, resulting in an axial displacement of $\mathbf{u}=u_0{\mathbf{e}}_z$ at the upper bonded natural rubber surface, the equivalent cylinder height and radius in terms of those of the original, natural state cylinder $\mathcal{B}$, read

$$h_0\left(u_0\right)=H-u_0$$

(1)

and

$$r_0\left(u_0\right)=\frac{d_0\left(u_0\right)}{2}=R\sqrt{\frac{H}{H-u_0}},$$

(2)

respectively, where $d_0$ is diameter.

2.4. Quasi-Static Natural Rubber Isolator Pre-Compression

The non-linear relation between the quasi-static axial force $F_0$ needed to pre-compress the natural rubber vibration isolator a quasi-static axial displacement of $u_0$ [77] reads

$$F_0=\frac{3\pi }{4}\frac{{\mu }_{\infty }{D}^2{u}_0}{H-u_0}\left[1+\frac{D^2(2H-u_0)}{16H^2(H-u_0)}\right]$$

(3)

and the corresponding quasi-static stiffness is

$$k_0=\frac{\mathrm{d}{F}_0}{\mathrm{d}{u}_0}=\frac{3\pi }{4}\frac{{\mu }_{\infty }{D}^{2}H}{{(H-u_0)}^2}\left[1+\frac{D^2}{8H(H-u_0)}\right],$$

(4)

where ${\mu }_{\infty }$ is the equilibrium (static) shear modulus. An interested reader is referred to Kari [77] regarding the derivations and validations of the previous quasi-static formulas. The relation (3) is advantageously applied while transforming a known axial pre-load into a quasi-static axial displacement, and vice versa, whilst the relation (4) is successfully applied while finding the quasi-static stiffness at a known axial pre-compression, and vice versa.

2.5. General Solution for the Superimposed Motion

The equation of infinitesimal motion in $\mathcal{S}$ for an isotropic viscoelastic material in absence of body forces [78,79] reads

$$\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}={\int }_{-\infty }^t\left[\kappa (t-\tau )+\frac{\mu (t-\tau )}{3}\right]\frac{\partial \nabla \left[\nabla ·\mathbf{u}\left(\tau \right)\right]}{\partial \tau }\mathrm{d}\tau +{\int }_{-\infty }^t\mu (t-\tau )\frac{\partial {\nabla }^2\mathbf{u}\left(\tau \right)}{\partial \tau }\mathrm{d}\tau ,$$

(5)

where $\mathbf{u}$ is the vector displacement field in $\mathcal{S}$; ${\nabla }^2$ is the Laplacian; $\rho$ is the mass per unit volume; and $\kappa$ and $\mu$ are the bulk and shear relaxation function, respectively. A convenient bulk relaxation for natural rubber reads

$$\kappa \left(t\right)=\alpha {\mu }_{\infty }h\left(t\right),$$

(6)

due to the near incompressibility of rubber, where the Heaviside step function is h, the nearly incompressibility material parameter is $\alpha \gg 1$, and the long-term relaxation function (equilibrium shear modulus) is ${\mu }_{\infty }={lim}_{t\to \infty }\mu \left(t\right)$. Furthermore, a suitable shear relaxation function reads

$$\mu \left(t\right)=\left[1+▵E_{\beta }\left(-{\left\{\frac{t}{\tau }\right\}}^{\beta }\right)\right]{\mu }_{\infty }h\left(t\right),$$

(7)

while modelling a fractional standard linear solid [14] in shear, where the Mittag–Leffler function is

$$E_{\beta }\left(x\right)=\sum _{n=0}^{\infty }\frac{x^n}{\Gamma (1+n\beta )},$$

(8)

where the Gamma function $\Gamma \left(z\right)={\int }_0^{\infty }{s}^{z-1}{e}^{-s}\mathrm{d}s$, $z>0$, and where the non-dimensional relaxation density $▵\gg 1:{lim}_{t\to 0^{+}}\mu \left(t\right)={\mu }_{\infty }(1+▵)$, the fractional derivative order $\beta :0<\beta ⩽1$, and the generalized relaxation time $\tau ⩾0$ are material parameters.

The frequency domain solution to the equation can be expressed in terms of scalar and vector potential functions $\tilde{\varphi }$ and $\tilde{\mathit{\psi }}$, respectively, where $\left(\widetilde{·}\right)={\int }_{-\infty }^{\infty }(·)exp(-\mathrm{i}\omega t)\mathrm{d}t$ denotes temporal Fourier transformation, $\mathrm{i}$ is an imaginary unit, and $\omega$ is an angular frequency, according to the Helmholtz theorem [79]

$$\tilde{\mathbf{u}}=grad\tilde{\varphi }+curl\tilde{\mathit{\psi }},$$

(9)

where $grad\tilde{\varphi }=\nabla \tilde{\varphi }$, in which $\tilde{\varphi }$ and $\tilde{\mathit{\psi }}$ satisfy the Helmholtz equations:

$${\nabla }^2\tilde{\varphi }+k_{\mathrm{L}}^2\tilde{\varphi }=0$$

(10)

and

$${\nabla }^2\tilde{\mathit{\psi }}+k_{\mathrm{T}}^2\tilde{\mathit{\psi }}=\mathbf{0},$$

(11)

and where $k_{\mathrm{L}}$ and $k_{\mathrm{T}}$ are the longitudinal and transversal wave number, respectively. They read

$$k_{\mathrm{L}}=\omega \sqrt{\frac{3\rho }{3\widehat{\kappa }+4\widehat{\mu }}}$$

(12)

and

$$k_{\mathrm{T}}=\omega \sqrt{\frac{\rho }{\widehat{\mu }}},$$

(13)

where the bulk modulus is

$$\widehat{\kappa }=\alpha {\mu }_{\infty }$$

(14)

and the shear modulus is

$$\widehat{\mu }=\left(1+\frac{▵{\left(\tau \mathrm{i}\omega \right)}^{\beta }}{1+{\left(\tau \mathrm{i}\omega \right)}^{\beta }}\right){\mu }_{\infty }.$$

(15)

In passing, alternative and equivalent forms of the equilibrium shear modulus and the non-dimensional relaxation density are ${\mu }_{\infty }={lim}_{\omega \to 0}\widehat{\mu }\left(\omega \right)$ and $▵\gg 1:{lim}_{\omega \to \infty }\widehat{\mu }\left(\omega \right)={\mu }_{\infty }(1+▵)$, respectively. The shear modulus (15) is a fractional extension of the shear modulus to the classical standard linear solid, where the viscous element with the stress proportional to the time derivative of the strain is substituted with an element with the stress proportional to the fractional time derivative of order $\beta$ of the strain. Likewise, the exponential stress relaxation function of the classical standard linear solid is substituted with the Mittag–Leffler stress relaxation function for the fractional standard linear solid. The models coincide as $\beta \to 1$. The interested reader is referred to Koeller [14] for more details.

The displacement field is provided by the Helmholtz theorem (9), reading, in a cylindrical co-ordinate system,

$${\tilde{u}}_r=\frac{\partial \tilde{\varphi }}{\partial r}+\mathrm{i}{k}_z{\tilde{\psi }}_{\phi },$$

(16)

$${\tilde{u}}_{\phi }=0$$

(17)

and

$${\tilde{u}}_z=-\mathrm{i}{k}_z\tilde{\varphi }+\frac{{\tilde{\psi }}_{\phi }}{r}+\frac{\partial {\tilde{\psi }}_{\phi }}{\partial r},$$

(18)

for the axially symmetric non-torsional case, where the axial dependence is assumed to be $\propto exp(-\mathrm{i}{k}_{z}z)$ and $k_z$ is the axial wave number. The corresponding stress field is

$${\tilde{\sigma }}_{rr}=\left(\alpha {\mu }_{\infty }-\frac{2\widehat{\mu }}{3}\right)\left(\frac{3\alpha {\mu }_{\infty }+4\widehat{\mu }}{3\alpha {\mu }_{\infty }-2\widehat{\mu }}\frac{\partial {\tilde{u}}_r}{\partial r}+\frac{{\tilde{u}}_r}{r}-\mathrm{i}{k}_z{\tilde{u}}_z\right),$$

(19)

$${\tilde{\sigma }}_{r\phi }={\tilde{\sigma }}_{\phi r}=0,$$

(20)

$${\tilde{\sigma }}_{rz}={\tilde{\sigma }}_{zr}=\widehat{\mu }\left(-\mathrm{i}{k}_z{\tilde{u}}_r+\frac{\partial {\tilde{u}}_z}{\partial r}\right),$$

(21)

$${\tilde{\sigma }}_{\phi \phi }=\left(\alpha {\mu }_{\infty }-\frac{2\widehat{\mu }}{3}\right)\left(\frac{\partial {\tilde{u}}_r}{\partial r}+\frac{3\alpha {\mu }_{\infty }+4\widehat{\mu }}{3\alpha {\mu }_{\infty }-2\widehat{\mu }}\frac{{\tilde{u}}_r}{r}-\mathrm{i}{k}_z{\tilde{u}}_z\right),$$

(22)

$${\tilde{\sigma }}_{\phi z}={\tilde{\sigma }}_{z\phi }=0$$

(23)

and

$${\tilde{\sigma }}_{zz}=\left(\alpha {\mu }_{\infty }-\frac{2\widehat{\mu }}{3}\right)\left(\frac{\partial {\tilde{u}}_r}{\partial r}+\frac{{\tilde{u}}_r}{r}-\mathrm{i}{k}_z\frac{3\alpha {\mu }_{\infty }+4\widehat{\mu }}{3\alpha {\mu }_{\infty }-2\widehat{\mu }}{\tilde{u}}_z\right).$$

(24)

Using the traction-free boundary condition along the radial surface ${\partial }_{\mathrm{t}}\mathcal{S}$ in Figure 2, the dispersion relation for axially symmetric and non-torsional waves in an infinite cylinder [80,81] is given by

$$0=4k_z^2\left[k_{\mathrm{L}}^2-k_z^2\right]\chi \left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_z^2)}{H-u_0}}\right)+{\left[k_{\mathrm{T}}^2-2k_z^2\right]}^2\chi \left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_z^2)}{H-u_0}}\right)-2k_{\mathrm{T}}^2\left[k_{\mathrm{L}}^2-k_z^2\right],$$

(25)

where the Onoe function $\chi \left(x\right)=x{\mathrm{J}}_0\left(x\right)/{\mathrm{J}}_1\left(x\right)$ and ${\mathrm{J}}_s$ is the Bessel function of first kind and order s.

The dispersion relation solution (25) is calculated by a Newton–Raphson method, with initial values given by a winding integral method, where the search domain is split into branch-cut absent sub-domains [8]. The axial wave number reads $k_{z,n}$, where $n\in {\mathbb{Z}}_{+}$ and labels different solutions, arranged in increasing spatial order, as $|\Im k_{z,n}|\le |\Im k_{z,m}|;n<m$ and ${\mathbb{Z}}_{+}$ denotes the set of positive integers.

The potential fields solutions to Helmholtz equations (Equations (10) and (11)) can be formulated as follows [7]:

$$\tilde{\varphi }=\sum _{n=1}^{\infty }\left[A_n^{+}{e}^{-\mathrm{i}{k}_{z,n}z}+A_n^{-}{e}^{\mathrm{i}{k}_{z,n}z}\right]{\mathrm{J}}_0\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right),$$

(26)

$${\tilde{\psi }}_{\phi }=\sum _{n=1}^{\infty }\left[B_n^{+}{e}^{-\mathrm{i}{k}_{z,n}z}+B_n^{-}{e}^{\mathrm{i}{k}_{z,n}z}\right]{\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right)$$

(27)

and

$${\tilde{\psi }}_r={\tilde{\psi }}_z=0,$$

(28)

where $r\in [0,r_0\left(u_0\right)]$ and $z\in [-h_0\left(u_0\right)/2,h_0\left(u_0\right)/2]$, with the origin of the local system being taken at the center of the natural rubber isolator. The coefficients are interrelated by the traction free boundary condition as $A_n^{+}={\Lambda }_n\left(u_0\right)B_n^{+}$ and $A_n^{-}=-{\Lambda }_n\left(u_0\right)B_n^{-}$, where

$${\Lambda }_n\left(u_0\right)=\frac{k_{\mathrm{T}}^2-2k_{z,n}^2}{2\mathrm{i}{k}_{z,n}\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}}\frac{{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right)}{{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right)}.$$

(29)

The boundary conditions at the cylinder ends are satisfied by mode matching the above fields.

2.6. Boundary Conditions

The boundary conditions for the pre-compressed, transformed natural rubber isolator undergoing a superimposed infinitesimal motion $\Delta \phi :\mathcal{S}\to \mathcal{S}$ are vanishing shear and axial stress on ${\partial }_{\mathrm{t}}^1\mathcal{S}$ and ${\partial }_{\mathrm{t}}^2\mathcal{S}$ while restraining the radial and axial displacement on ${\partial }_{\mathrm{d}}^2\mathcal{S}$ and imposing an axial infinitesimal displacement of $-\Delta \tilde{u}{\mathbf{e}}_z$ as well as a vanishing radial displacement on ${\partial }_{\mathrm{d}}^1\mathcal{S}$. The boundary conditions hereby applied to $\overline{{\partial }_{\mathrm{d}}^1\mathcal{S}\cup {\partial }_{\mathrm{t}}^1\mathcal{S}\cup {\partial }_{\mathrm{t}}\mathcal{S}\cup {\partial }_{\mathrm{d}}^2\mathcal{S}\cup {\partial }_{\mathrm{t}}^2\mathcal{S}}=\partial \mathcal{S}$ are locally non-mixed and read

$$\tilde{\mathbf{u}}·{\mathbf{e}}_r=0$$

(30)

and

$$\tilde{\mathbf{u}}·{\mathbf{e}}_z=-\Delta \tilde{u},$$

(31)

given on ${\partial }_{\mathrm{d}}^1\mathcal{S}$;

$$\tilde{\mathbf{u}}·{\mathbf{e}}_r=0$$

(32)

and

$$\tilde{\mathbf{u}}·{\mathbf{e}}_z=0,$$

(33)

given on ${\partial }_{\mathrm{d}}^2\mathcal{S}$;

$${\mathbf{e}}_r·\left(\tilde{\mathit{\sigma }}{\mathbf{e}}_z\right)=0,$$

(34)

and

$${\mathbf{e}}_z·\left(\tilde{\mathit{\sigma }}{\mathbf{e}}_z\right)=0,$$

(35)

given on ${\partial }_{\mathrm{t}}^1\mathcal{S}$ and ${\partial }_{\mathrm{t}}^2\mathcal{S}$ and, finally,

$$\tilde{\mathit{\sigma }}{\mathbf{e}}_r=\mathbf{0},$$

(36)

given on ${\partial }_{\mathrm{t}}\mathcal{S}$, with the latter satisfied by the dispersion relation (25). Using (16), (18), (26), (27), and (29), the requirements (30) and (31) become

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}+C_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)V_n^r\left(r;u_0\right)=0$$

(37)

and

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}-C_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)V_n^z\left(r;u_0\right)=-\Delta \tilde{u},$$

(38)

at $z=-h_0\left(u_0\right)/2$, (32) and (33) become

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}+C_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)V_n^r\left(r;u_0\right)=0$$

(39)

and

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}-C_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)V_n^z\left(r;u_0\right)=0,$$

(40)

at $z=h_0\left(u_0\right)/2$, where the coefficients $X_n\left(u_0\right)C_n^{+}=B_n^{+}$ and $-X_n\left(u_0\right)C_n^{-}=B_n^{-}$ with

$$X_n\left(u_0\right)=2\mathrm{i}\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}{k}_{z,n}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right)$$

(41)

and

$$\begin{array}{c}{V}_n^r\left(r;u_0\right)=\left(2k_{z,n}^2-k_{\mathrm{T}}^2\right)\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-2k_{z,n}^2\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right)\end{array}$$

(42)

and

$$\begin{array}{c}{V}_n^z\left(r;u_0\right)=\mathrm{i}{k}_{z,n}\left[2k_{z,n}^2-k_{\mathrm{T}}^2\right]{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_0\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+2\mathrm{i}{k}_{z,n}\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_0\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right),\end{array}$$

(43)

$r\in [0,R]$. Similarly, using (16), (18), (21), (26), (27), and (29), the requirement (34) becomes

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}-C_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)S_n^{zr}\left(r;u_0\right)=0,$$

(44)

at $z=-h_0\left(u_0\right)/2$, and

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}-C_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)S_n^{zr}\left(r;u_0\right)=0,$$

(45)

at $z=h_0\left(u_0\right)/2$, where

$$\begin{array}{c}{S}_n^{zr}\left(r;u_0\right)=2\mathrm{i}\widehat{\mu }{k}_{z,n}\left(k_{\mathrm{T}}^2-2k_{z,n}^2\right)\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}[{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill -{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right)],\end{array}$$

(46)

$r\in ]R,r_0\left(u_0\right)]$. Finally, using (16), (18), (24), (26), (27), and (29), the requirement (35) becomes

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}+C_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)S_n^{zz}\left(r;u_0\right)=0,$$

(47)

at $z=-h_0\left(u_0\right)/2$, and

$$\sum _{n=1}^{\infty }\left(C_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}+C_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right)S_n^{zz}\left(r;u_0\right)=0,$$

(48)

at $z=h_0\left(u_0\right)/2$, where

$$\begin{array}{c}{S}_n^{zz}\left(r;u_0\right)=\widehat{\mu }\left(k_{\mathrm{T}}^2-2k_{\mathrm{L}}^2+2k_{z,n}^2\right)\left(2k_{z,n}^2-k_{\mathrm{T}}^2\right){\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_0\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill +4\widehat{\mu }{k}_{z,n}^2\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_0\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right),\end{array}$$

(49)

$r\in ]R,r_0\left(u_0\right)]$.

The simplest way is to determine the unknown coefficients in (37) to (40), (44), (45), (47), and (48) is to apply a point-matching technique, but in most cases, more accurate results are obtained using the subregion method, dividing the domain into sub-domains. The cylinder ends may be divided into $P_{\mathrm{d}}+P_{\mathrm{t}}\in {\mathbb{Z}}_{+}$ axially symmetric circular bands between the radii $0=a_0<a_1<a_2<\cdots <a_{P_{\mathrm{d}}}=R<a_{P_{\mathrm{d}}+1}<a_{P_{\mathrm{d}}+2}<\cdots <a_{P_{\mathrm{d}}+P_{\mathrm{t}}}=r_0\left(u_0\right)$. By performing integrations between subsequent radii and using the relations $\int x{\mathrm{J}}_0\left(x\right)\mathrm{d}x=x{\mathrm{J}}_1\left(x\right)$ and $\int x{\mathrm{J}}_1\left(x\right)\mathrm{d}x=\pi x({\Xi }_0\left(x\right){\mathrm{J}}_1\left(x\right)-{\Xi }_1\left(x\right){\mathrm{J}}_0\left(x\right))/2$, where ${\Xi }_j$ is a Struve function of order j [82], the conditions (37)–(40) read

$$\sum _{n=1}^N\left[D_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}+D_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[U_n^r(a_p;u_0)-U_n^r(a_{p-1};u_0)\right]=0,$$

(50)

$$\sum _{n=1}^N\left[D_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}-D_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[U_n^z(a_p;u_0)-U_n^z(a_{p-1};u_0)\right]=\pi \left[a_{p-1}^2-a_p^2\right],$$

(51)

$$\sum _{n=1}^N\left[D_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}+D_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[U_n^r(a_p;u_0)-U_n^r(a_{p-1};u_0)\right]=0$$

(52)

and

$$\sum _{n=1}^N\left[D_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}-D_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[U_n^z(a_p;u_0)-U_n^z(a_{p-1};u_0)\right]=0,$$

(53)

where

$$\begin{array}{c}{U}_n^r\left(r;u_0\right)={\pi }^{2}r[\left(k_{\mathrm{T}}^2-2k_{z,n}^2\right){\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right)\mathrm{W}\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill +2k_{z,n}^2\sqrt{\frac{k_{\mathrm{L}}^2-k_{z,n}^2}{k_{\mathrm{T}}^2-k_{z,n}^2}}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right)\mathrm{W}\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right)],\end{array}$$

(54)

$$\begin{array}{c}{U}_n^z\left(r;u_0\right)=2\mathrm{i}\pi k_{z,n}r[\frac{\left(2k_{z,n}^2-k_{\mathrm{T}}^2\right)}{\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill +2\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right)],\end{array}$$

(55)

$\mathrm{W}={\mathrm{J}}_0{\Xi }_1-{\mathrm{J}}_1{\Xi }_0$, $D_n^{+}\Delta \tilde{u}=C_n^{+},D_n^{-}\Delta \tilde{u}=C_n^{-}$ and $p=1,2,\dots ,P_{\mathrm{d}}$ with $\Delta \tilde{u}\ne 0$. Similarly, the conditions (44), (45), (47), and (48) read

$$\sum _{n=1}^N\left[D_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}-D_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[T_n^{zr}(a_q;u_0)-T_n^{zr}(a_{q-1};u_0)\right]=0,$$

(56)

$$\sum _{n=1}^N\left[D_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}-D_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[T_n^{zr}(a_q;u_0)-T_n^{zr}(a_{q-1};u_0)\right]=0,$$

(57)

$$\sum _{n=1}^N\left[D_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}+D_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[T_n^{zz}(a_q;u_0)-T_n^{zz}(a_{q-1};u_0)\right]=0,$$

(58)

and

$$\sum _{n=1}^N\left[D_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}+D_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]\left[T_n^{zz}(a_q;u_0)-T_n^{zz}(a_{q-1};u_0)\right]=0,$$

(59)

where

$$\begin{array}{c}{T}_n^{zr}\left(r;u_0\right)=2{\pi }^2\mathrm{i}{k}_{z,n}r\widehat{\mu }\left(k_{\mathrm{T}}^2-2k_{z,n}^2\right)\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}[{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right)\frac{\mathrm{W}\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right)}{\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill -{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right)\frac{\mathrm{W}\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)}{\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}}],\end{array}$$

(60)

$$\begin{array}{c}{T}_n^{zz}\left(r;u_0\right)=2\pi r\widehat{\mu }[\left(k_{\mathrm{T}}^2-2k_{\mathrm{L}}^2+2k_{z,n}^2\right)\left(2k_{z,n}^2-k_{\mathrm{T}}^2\right){\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{T}}^2-k_{z,n}^2)}{H-u_0}}\right)\frac{{\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}\right)}{\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill +4k_{z,n}^2\sqrt{k_{\mathrm{L}}^2-k_{z,n}^2}{\mathrm{J}}_1\left(R\sqrt{\frac{H(k_{\mathrm{L}}^2-k_{z,n}^2)}{H-u_0}}\right){\mathrm{J}}_1\left(r\sqrt{k_{\mathrm{T}}^2-k_{z,n}^2}\right)]\end{array}$$

(61)

and $q=P_{\mathrm{d}}+1,P_{\mathrm{d}}+2,\dots ,P_{\mathrm{d}}+P_{\mathrm{t}}$. To obtain a finite number of equations, the series (50) to (53) and (56) to (59) are truncated after $N\in {\mathbb{Z}}_{+}$ terms, where $P_{\mathrm{d}}+P_{\mathrm{t}}\ge N/2$.

The relations (50) to (53) and (56) to (59) may be written in a matrix form as $\mathbf{A}\mathbf{x}=\mathbf{b}$, where $\mathbf{A}$ is a $4(P_{\mathrm{d}}+P_{\mathrm{t}})\times 2N$ system matrix, $\mathbf{x}$ is an unknown coefficient vector ${(D_1^{+},D_1^{-},D_2^{+},D_2^{-},\dots ,D_N^{+},D_N^{-})}^{\mathrm{T}}$, and $\mathbf{b}$ is a known right-hand vector

$$\mathbf{b}=\pi {(\underset{P_{\mathrm{d}}}{\underbrace{0,0,\dots ,0}},\underset{P_{\mathrm{d}}}{\underbrace{a_0^2-a_1^2,a_1^2-a_2^2,\dots ,a_{P_{\mathrm{d}}-1}^2-a_{P_{\mathrm{d}}}^2}},\underset{2P_{\mathrm{d}}+4P_{\mathrm{t}}}{\underbrace{0,0,\dots ,0}})}^{\mathrm{T}},$$

(62)

where T is a transpose. The equation system is solved by singular-value decomposition after rescaling and equilibration [7].

Finally, the transfer and driving point stiffness read

$${\tilde{k}}_{\mathrm{t}}\left(\omega ;u_0\right)=\frac{1}{\Delta \tilde{u}}{\int }_{{\partial }_{\mathrm{d}}^2\mathcal{S}}{\tilde{\sigma }}_{zz}\left(\omega ;u_0\right)\mathrm{d}\left(\partial \mathcal{S}\right)$$

(63)

and

$${\tilde{k}}_{\mathrm{dp}}\left(\omega ;u_0\right)=\frac{1}{\Delta \tilde{u}}{\int }_{{\partial }_{\mathrm{d}}^1\mathcal{S}}{\tilde{\sigma }}_{zz}\left(\omega ;u_0\right)\mathrm{d}\left(\partial \mathcal{S}\right)-m_{\mathrm{mp}}{\omega }^2,$$

(64)

respectively, where frequency and pre-compression dependencies are explicitly stated with an assumed rigid body motion of the metal mass $m_{\mathrm{mp}}$. Consequently, they read

$${\tilde{k}}_{\mathrm{t}}\left(\omega ;u_0\right)=\sum _{n=1}^N{T}_n(\omega ;u_0)\left[D_n^{+}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}+D_n^{-}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]$$

(65)

and

$${\tilde{k}}_{\mathrm{dp}}\left(\omega ;u_0\right)=\sum _{n=1}^N{T}_n(\omega ;u_0)\left[D_n^{+}{e}^{\mathrm{i}{k}_{z,n}(H-u_0)/2}+D_n^{-}{e}^{-\mathrm{i}{k}_{z,n}(H-u_0)/2}\right]-m_{\mathrm{mp}}{\omega }^2,$$

(66)

where $k_{z,n}=k_{z,n}(\omega ;u_0)$, $D_n^{+}=D_n^{+}(\omega ;u_0)$, $D_n^{-}=D_n^{-}(\omega ;u_0)$ and $T_n(\omega ;u_0)\equiv T_n^{zz}(R;u_0)$, with the latter given by expression (61) with $r=R$. The series (65) and (66) model the dynamic stiffness to any desired accuracy for a sufficient number of modes N while satisfying the traction free boundary condition at the radial surface exactly and the displacement and traction boundary conditions at the cylinder ends in mean. They globally represent the dynamic stiffness for a pre-compressed natural rubber vibration isolator as given by a waveguide model, thereby reducing an original, mathematically arduous, and complex non-linear problem into a vastly simpler assignment.

3. Results and Discussion

The formulas above are implemented in Lahey Fortran 90^® (Lahey Computer Systems Inc., Incline Village, NV, USA) with all calculations performed in double precision and the results presented using Matlab^® (MathWorks Inc., Natick, MA, USA). Two test objects are studied with a diameter of $D=100$ mm and heights of $H=25$ and $50$ mm, respectively. The two test objects selected represent a normal range of shape factors $\mathcal{S}=D/4L$ [1] for this vastly common vibration isolator—with $\mathcal{S}$ spanning from 0.5 to 1.

3.1. Natural Rubber Material

The natural rubber material applied in this numerical investigation is supposed to be identical to that of Kari et al. [20]—namely, unfilled sulphur cured natural rubber General Purpose Standard Malaysian Rubber (GP SMR) with the ingredients given in

Table 1. Ingredients of the natural rubber material and their concentrations [20].

Ingredient	Kind	Concentration [phr]
Natural rubber	GP SMR	100
Stabilizers	Wax	1
	Antioxidant	1
	Antiozonant	1
Activators	Zinc oxide	5
	Stearic acid	1
Vulcanizing agent	Sulphur	3
Accelerator	CBS ${}^1$	2
Processing oils	Aromatic	5
	Paraffinic	1

¹n-Cyclohexyl-2-Benzothiazole Sulfenamide.

. The fitted material parameters are given in

Table 2. Applied natural rubber material parameters in the model [20].

Quantity	Variable	Value
Density	$\rho$	984 kg/m${}^3$
Nearly incompressible parameter	$\alpha$	2222
Equilibrium shear modulus	${\mu }_{\infty }$	825 kN/m${}^2$
Relaxation density	Δ	276
Generalized relaxation time	$\tau$	2.94 ns
Fractional derivative order	$\beta$	0.657

. The interested reader is referred to Kari [20] for more details regarding the material manufacturing process, the measurement methods and equipments, and the material parameter fitting methods.

3.2. Model Results

To obtain accurate stiffness results, the number of modes is required to be large; $N=128$ is sufficient, while the equation system is taken to be threefold over-determined; $P=P_{\mathrm{d}}+P_{\mathrm{t}}=256$, with subregion applied radii

$$a_{\xi }\left(u_0\right)=\frac{\xi R}{P}\sqrt{\frac{H}{H-u_0}},$$

(67)

for all integers $\xi \in [0,P]$ obeying $\xi \ne P_{\mathrm{d}}$, together with radius $a_{P_{\mathrm{d}}}=R$, where the specific $P_{\mathrm{d}}\left(u_0\right)\in [0,P]$ satisfies

$$\underset{P_{\mathrm{d}}}{min}\left|\frac{P_{\mathrm{d}}}{P}\sqrt{\frac{H}{H-u_0}}-1\right|.$$

(68)

The expression (67) slices the radius from 0 to $r_0$ (see Equation (2)) into $(P-2)$ equally distant segments. The remaining two segments between the points $(P_{\mathrm{d}}-1)$ to $P_{\mathrm{d}}$ and $P_{\mathrm{d}}$ to $(P_{\mathrm{d}}+1)$ are forced to coincide with R at the point $P_{\mathrm{d}}$. The latter is performed by setting $a_{P_{\mathrm{d}}}=R$, with $P_{\mathrm{d}}$ given by the expression (68).

For all frequencies considered, 20 to $2000$Hz with a step of $10$Hz, the effective rank of the conditionalized system matrix is full, that is, 256. The frequency range selected covers two important frequency decades within the audible frequency range.

The stiffness magnitude $|\tilde{k}|$ and phase $180\angle \tilde{k}/\pi$ for $H=25$ mm natural rubber isolator at pre-compressions $u_0=0,1,2,3,4$ and $5$ mm—that is, from vanishing to 20% compression—are presented in Figure 3 and Figure 4, where $\tilde{k}=\left|\tilde{k}\right|exp(\mathrm{i}\angle \tilde{k})$. Figure 3 displays the transfer stiffness, while the driving point stiffness is in Figure 4. The direction for increased pre-compression is marked by arrows. The corresponding transfer and driving point stiffness results for $H=50$ mm natural rubber isolator are presented in Figure 5 and Figure 6 at pre-compressions $u_0=0,2,4,6,8$, and $10$ mm. The driving point stiffness is presented here and in the following without the plate mass contribution. To include this contribution, ${\tilde{k}}_{\mathrm{dp}}{|}_{\mathrm{no}\mathrm{mass}}$ is simply replaced by ${\tilde{k}}_{\mathrm{dp}}{{|}_{\mathrm{mass}}={\tilde{k}}_{\mathrm{dp}}|}_{\mathrm{no}\mathrm{mass}}-m_{\mathrm{mp}}{\omega }^2$, as in Equation (66). The transfer stiffness is however independent of that mass.

Clearly, the low-frequency stiffness magnitude in Figure 3 is plateau-like at the vanishing pre-compression. Then, the curve drops steeply into a trough at $860$Hz and rises sharply while displaying minor oscillations. With respect to the magnitude curve and to the behavior of the phase curve, the trough most likely corresponds to a resonance. At a resonance for elastic materials (that is, with no material damping), the stiffness shows a magnitude trough and a phase jump of $+{180}^{\circ }$, while displaying a magnitude peak and a phase jump of $-{180}^{\circ }$ at an anti-resonance. Introducing damping, the resonances and anti-resonances are blunted; in general, the magnitude troughs and peaks are rounded, while the sudden phase jumps disappear, showing a ‘slower’ phase shift. Therefore, it may become difficult to distinguish individual resonances and anti-resonances, particularly for high damping material (such as rubber) at closely spaced resonances and anti-resonances. In passing, it may be observed that it is not possible to predict a transfer stiffness resonance by traditional, simplified models—such as the classical long rod example [2,3,4,5]—in contrast to the waveguide approach applied in this paper.

Regarding pre-compression dependence in Figure 3, the low-frequency magnitude plateau shows an increase in compression, as expected. This magnitude increase is also evident in the high-frequency region. However, the trough frequency is almost pre-compression independent, only showing an initial, slight frequency decrease with compression.

The driving point stiffness result in Figure 4 displays a similar behavior as its transfer counterpart. However, the stiffness trough frequency is lower, $570$Hz at the vanishing pre-compression, while displaying a clear frequency increase with compression.

The stiffness results for $H=50$ mm natural rubber isolator in Figure 5 and Figure 6 are similar to those of the shorter natural rubber isolator: displaying a low-frequency stiffness magnitude plateau that increases with pre-compression. However, the transfer stiffness magnitude in Figure 5 rises steeply to peak at the right end of the plateau region while displaying a strong frequency dependence—resulting in peaks and troughs—after this initial anti-resonance. On the other hand, the driving point stiffness magnitude in Figure 6 drops into a deep trough after the low-frequency plateau region—in line with the shorter natural rubber isolator—while rising steeply after this initial resonance, though eventually showing a stronger high-frequency oscillation compared to that of the shorter natural rubber isolator within the considered frequency range. The driving point stiffness resonance and anti-resonance succeed each other repeatedly in Figure 4 and Figure 6, as $0\le \angle {\tilde{k}}_{\mathrm{dp}}\le \pi$ from thermodynamically point of view. However, this alternating order is not necessary for the transfer stiffness in Figure 3 and Figure 5 as no such phase restriction exists for that stiffness.

3.3. Finite Element Comparison

The developed waveguide model for a pre-compressed natural rubber isolator is verified by calculating the dynamic stiffness of the original, pre-compressed natural rubber isolator in Figure 1 using a non-linear finite element method [11]. To this end, a hyperelastic neo-Hookean material model with a fractional standard linear solid is adopted for the studied, unfilled rubber, where parameters are fitted to match those of the waveguide model under infinitesimal strains, that is, resulting in identical shear modulus, bulk modulus, and density. The mesh employed consists of $3000(50\times 60)$ axially symmetric hybrid elements, gradually refined towards the rubber cylinder corners. The force needed for the stiffness calculation is evaluated by spatial integration of extrapolated Gauss point stresses at massless, almost rigid, elastic endplates bonded to the rubber cylinder, thereby avoiding problems at the rubber corners while pre-compressed, resulting in enhanced force estimations. Furthermore, the Lagrangian 9-node (biquadratic displacement–bilinear pressure) and the 4-node (bilinear displacement–constant pressure) hybrid elements, result in infinitesimal stiffness differences, thereby in practice verifying the non-linear finite element solution convergence. The interested reader is referred to Kari [11] for more details concerning the implementation and validation of the non-linear finite element method.

The transfer and driving point stiffness magnitude and phase for $H=25$ mm natural rubber isolator at vanishing and maximum pre-compressions; that is, $u_0=0$ and $5$ mm, are presented in Figure 7 and Figure 8. The corresponding stiffness results for $H=50$ mm natural rubber isolator are presented in Figure 9 and Figure 10 at pre-compressions $u_0=0$ and $10$ mm.

Obviously, the waveguide model results closely follow the non-linear finite element solutions in Figure 7, Figure 8, Figure 9 and Figure 10. The stiffness magnitude and phase curves are in practice overlapping at the vanishing pre-compression throughout the whole frequency range. Furthermore, the waveguide stiffness magnitude curves display a slight low-frequency over-estimation while both waveguide stiffness magnitude and phase curves display comparatively close high-frequency agreements with the corresponding non-linear finite element results at the maximum pre-compression. The developed waveguide model presented in this paper is thereby in practice verified. Furthermore, it may be noted that the close agreement between the waveguide model results and those of the non-linear finite element model at the maximum pre-compressions implies that the non-linear geometrical effects prevails over the non-linear material effects in this study since the waveguide model presumes linear material as opposed to the finite element model. This is not surprisingly as the non-linear material model applied is hyperelastic neo-Hookean, being the straightforward non-linear extension of the linear elastic Hookean model and a suitable model for the unfilled natural rubber material virtually employed.

3.4. Convergence Properties

The stiffness convergence of the presented waveguide model is investigated by extending the number of radii $P=P_{\mathrm{d}}+P_{\mathrm{t}}$ and modes N.

The P-extension is processed by calculating the stiffness for $P=2^{\gamma }$, where $\gamma =6,7,8$, while using $N=128$; that is, the equation system is exactly determined, single, and threefold over-determined, respectively. The transfer and driving point stiffness magnitude for $H=25$ mm natural rubber isolator at maximum pre-compression $u_0=5$ mm are presented in Figure 11. The corresponding stiffness results for $H=50$ mm natural rubber isolator are presented in Figure 12 at maximum pre-compression $u_0=10$ mm.

Clearly, the $P=64$ stiffness curves in Figure 11 and Figure 12 oscillate slightly around the well-behaved, almost coinciding $P=128$ and 256 curves. In practice, it is therefore sufficient to apply a single-fold or threefold over-determined equation system for the stiffness calculation of the object under study focus, while an exactly determined equation system results in an insufficient stiffness solution.

The N-extension is processed by calculating the stiffness for $N=2^{\epsilon }$, where $\epsilon =4,5,6,7$, while using a threefold over-determined equation system; that is, the number of radii is concurrently extended according to $P=2N$. The transfer and driving point stiffness magnitude for $H=25$ mm natural rubber isolator at maximum pre-compression $u_0=5$ mm are presented in Figure 13, while the corresponding results for $H=50$ mm natural rubber isolator are presented in Figure 14 at maximum pre-compression $u_0=10$ mm.

Clearly, the stiffness converges for larger N, and the difference between the $N=64$ and 128 solutions is negligible. In passing, it should be noted that variations in manufacturing process parameters such as curing time, temperature, and pressure; the unbridled compound ingredients; and mixing result in stiffness variations of the individual natural rubber vibration isolators, thus, in practice, limiting any exact agreements between the modelling and measurement results of the individual samples [83]. Furthermore, although the stiffness solution is shown to converge for larger N, the necessary number of modes is rather high: $N\ge 64$. Physically, these high-order modes contribute to the fulfillment of the displacement and stress boundary conditions at the cylinder ends. The stress field at these ends in turn determines the force needed for the dynamic stiffness calculation, with the latter being the global variable for which the convergence properties are investigated. In practice, the stiffness solution convergence of the developed waveguide model is thereby verified.

3.5. Edge Boundary Condition Influence

The applied boundary condition on the edges ${\partial }_{\mathrm{t}}^1\mathcal{S}$ and ${\partial }_{\mathrm{t}}^2\mathcal{S}$ for the pre-compressed natural rubber isolator in Figure 2 are free, that is, vanishing shear and axial stress. However, other boundary conditions may be imposed at the edges: fixed, that is, imposing an axial infinitesimal displacement of $-\Delta \tilde{u}{\mathbf{e}}_z$ and vanishing radial displacement on ${\partial }_{\mathrm{t}}^1\mathcal{S}$ and vanishing axial and radial displacement on ${\partial }_{\mathrm{t}}^2\mathcal{S}$—physically representing an extension of the constraining metal plates over the whole end surfaces of the rubber cylinder in Figure 2; slip, that is, imposing an axial infinitesimal displacement of $-\Delta \tilde{u}{\mathbf{e}}_z$ and vanishing shear stress on ${\partial }_{\mathrm{t}}^1\mathcal{S}$ and vanishing axial displacement and shear stress on ${\partial }_{\mathrm{t}}^2\mathcal{S}$; and no-edge, that is, assuming a preserved radius of the natural rubber isolator during pre-compression $u_0$. Although the latter condition is simple to apply, it is not physically attractive as the incompressibility condition is not fulfilled. The transfer and driving point stiffness magnitude for $H=25$ mm natural rubber isolator at maximum pre-compression $u_0=5$ mm are presented in Figure 15 while using the different boundary conditions—fixed, slip, free, FEM, and no-edge, where FEM represents the non-linear finite element solution from Section 3.3 of the original, pre-compressed natural rubber isolator—from Figure 1. The corresponding stiffness results for $H=50$ mm natural rubber isolator are presented in Figure 16 at maximum pre-compression $u_0=10$ mm. The number of modes applied is $N=128$ while using a threefold over-determined equation system with $P=256$.

Clearly, the waveguide model applying the free edge boundary condition predicts a stiffness close to that of the non-linear finite element model in Figure 15 and Figure 16. Indeed, its solution is very close to that of the finite element model almost throughout the whole frequency region—deviating only slightly within a small frequency range—whereas the fixed and slip boundary conditions overestimate the low-frequency stiffness magnitude while the no-edge condition underestimates the same. Not surprising, as a fixed and slip boundary condition should—while constraining the motion more than the free boundary condition—result in a statically stiffer natural rubber isolator whereas the thinner, no-edge natural rubber isolator results in a softer natural rubber isolator. Furthermore, the no-edge solution underestimates the stiffness magnitude in the high-frequency end, as is clearly shown in Figure 15 and Figure 16.

In summary, the developed waveguide model presented in this paper, applying a free(-free) edge boundary condition on ${\partial }_{\mathrm{t}}^1\mathcal{S}$ and ${\partial }_{\mathrm{t}}^2\mathcal{S}$, for the pre-compressed natural rubber isolator in Figure 2, is hereby in practice shown to result in a superior stiffness estimate.

3.6. Displacement and Stress Fields

The displacement and stress fields at the cylinder ends and $2000$Hz are calculated by using $N=128$ modes and a threefold over-determined equation system with $P=256$. The normalized axial and radial displacement ${\tilde{u}}_z/\Delta \tilde{u}$ and ${\tilde{u}}_r/\Delta \tilde{u}$, respectively, at $z=\mp h_0\left(u_0\right)/2$ for $H=25$ mm natural rubber isolator versus normalized radius $r/R$, at maximum pre-compression $u_0=5$ mm, are presented in Figure 17. The corresponding fields for $H=50$ mm natural rubber isolator, at maximum pre-compression $u_0=10$ mm, are presented in Figure 18.

Clearly, the displacement boundary conditions are well satisfied for $0\le r/R\le 1$ in Figure 17 and Figure 18, displaying only a slight oscillation close to $r/R\approx 1$.

Likewise, the normalized axial and shear stress ${\tilde{\sigma }}_{zz}/{\mu }_{\infty }/\Delta \tilde{u}$ and ${\tilde{\sigma }}_{zr}/{\mu }_{\infty }/\Delta \tilde{u}$, respectively, at $z=\mp h_0\left(u_0\right)/2$ for $H=25$ mm natural rubber isolator versus normalized radius $r/R$, at maximum pre-compression $u_0=5$ mm, are presented in Figure 19. The corresponding fields for $H=50$ mm natural rubber isolator, at maximum pre-compression $u_0=10$ mm, are presented in Figure 20.

Clearly, the stress boundary conditions are well satisfied for $1<r/R\le r_0\left(u_0\right)/R$ in Figure 19 and Figure 20. In addition, the shear stress ${\tilde{\sigma }}_{zr}={\tilde{\sigma }}_{rz}$ vanishes identically at $r/R=r_0\left(u_0\right)/R$ due to the stress boundary condition applied at the radial surface while deriving the dispersion relation (25). However, the axial stress is singular $|{\tilde{\sigma }}_{zz}|\to \infty$ at $r/R=1$. This point value is impossible to satisfy using the proposed method while using non-singular terms unless $N\to \infty$. However, the axial force needed for the dynamic stiffness evaluation is determined by a spatial integration of the axial stress (see Equations (63) and (64)) which is not that sensitive to singular points, as is clearly shown in Section 3.3 and Section 3.4.

In addition to the unfilled natural rubber material studied in this paper, the developed model is most likely suitable for other incompressible and nearly incompressible polymer materials where, moreover, the Mittag-Leffler stress relaxation function is suitable, such as for styrene-butadiene rubber (SBR) [45] and tough, doubly cross-linked, single network hydrogels with both chemical and physical cross-links [64] and where, furthermore, the geometrical non-linearity prevails over the material non-linearity, such as for not overly filled polymers.

4. Conclusions

A waveguide model is developed for the dynamic stiffness of pre-compressed natural rubber vibration isolators by transforming the isolator into a globally equivalent homogeneous and isotropic form, where the boundary conditions are satisfied by mode matching. The conclusions are as follows:

• The dynamic stiffness solutions are shown to converge for a moderately high number of modes—such as $N\ge 64$ for the studied vibration isolators within the considered frequency range 20 to 2000 Hz and pre-compression range up to 20%;
• The dynamic stiffness is found to depend strongly on the frequency displaying resonance phenomena such as peaks and troughs;
• The dynamic stiffness is found to depend strongly on the pre-compression, displaying low-frequency stiffness magnitude stiffness increase in addition to peak and trough shifts with increased pre-compressions;
• The waveguide model is shown to yield dynamic stiffness results close to those of the non-linear finite element method;
• The non-linear finite element method was previously shown to give results close to those of the experiments for similar natural rubber vibration isolators [11];
• The applied Mittag–Leffler shear relaxation function was previously shown to give shear modulus results close to those of the experiments for the same natural rubber material as studied in this paper [20]; and
• The free-free cylinder edge boundary condition is shown to result in superior pre-compressed dynamic stiffness results.

An interesting continuation of the work presented is to apply the proposed method to other excitation directions, such as torsional and shear, and to other natural rubber vibration isolators, such as strips and rectangular shaped natural rubber vibration isolators. Another possible extension is to additionally include physically and chemically aging of the natural rubber, for example, by applying shear modulus models that incorporate physical and chemical aging [52,53].