Frequency-dependent viscoelastic models for passive vibration isolation systems

Even though there is a growing interest in active vibration isolation systems, passive approaches are still the best choice in many cases because they are inherently the simplest and of lowest cost. Moreover, better comprehension of the dynamics and specially of the damping behavior in passive systems is required for successful implementation of active schemes. In the vast literature of passive isolation systems, there are not many works that consider damping models more elaborated than the widely used complex modulus. In this work a passive isolation system composed of a base and two isolators, modelled as Timoshenko beams, and a vibration source, modelled as a rigid body, is considered. For the isolators, two different viscoelastic models are considered: the Anelastic Displacement Fields (ADF) and Fractional Calculus (FC), which will be compared with the complex modulus model. The results show that both ADF and FC models lead to better approximation of dissipated energy, since they account for frequency-dependence of the viscoelastic isolators. Analysis of the curve-fitting of material parameters, using ADF and FC models has shown that generally less parameters are needed by FC model, for the same fitting quality, although optimization results depends strongly on the initial guess for the solution.


Introduction
Vibration isolation systems are used in a wide variety of applications in order to reduce the transmission of mechanical vibrations from the equipments (source) to the foundation (base) or vice versa.Isolation is obtained by inserting a mechanical component (isolator) that links the source and the base.Hence, it is very important to accurately model the system source/base/isolator.However, it is not easy to find a model that describes suitably the response at all frequencies.Single-degree-of-freedom models are not satisfactory and even models considering multiple mounts, each one treated as a lumped system, are incomplete.Additionally, the material stiffness E and loss factor η normally vary with the frequency, this is specially noticeable in isolators made of viscoelastic materials.In this work, the simple frequencyindependent complex modulus model is compared with Anelastic Displacement Fields (ADF) [6] and Fractional Calculus (FC) [1] models, which are able to predict the response even in highly damped structures.
On the other hand, there are many methods used to model flexible bodies, the Finite Element Method (FEM) is largely the most widespread, however alternative methods like the Spectral Element Method (SEM) [4] may be most suitable for analysis in medium and high frequency ranges, mainly because it has the advantage of considering exact solutions in the frequency domain.These two modelling methods can be easily implemented in conjunction with the viscoelastic models aforementioned.
The object of the work is to present an analysis of an isolation system with two isolators made of viscoelastic material.We pay special attention to the model formulation and the frequency-dependent viscoelastic material properties.

Isolation system model
Considering the Fig. 1, expressions based in a mobility-impedance formulation will be introduced.
All the forces ( ) have 3 components, though, in general, we show just one for simplicity.Note that the displacements, velocities, accelerations and forces will be represented by x = Xe iωt , where x is a harmonically time-varying quantity and X is an amplitude that in the more general case will be a complex number.Each component is a quantity that varies harmonically, so it can be written in the following way: Expressions that relate the velocities ( v) and the forces (f ) in the source and in the base can be expressed as [5] Where we have the relations The following notation will be used: the mobilities (M) with subscript 1 correspond to the internal forces (f s1 , f s2 , f b1 , f b2 ) and those with subscript 2 to the external forces (q s , q b ).Meanwhile, the velocities ( v) or the forces (f ) with subscript 1 correspond to the left isolator and the ones with subscript 2 to the right isolator.
For the isolators, the forces and the velocities in their ends will be related by where Z i1 is the impedance matrix.Moreover Then, the first two equations can be grouped as Using the action-reaction principle for the forces and continuity for the velocities, we obtain two additional equations T is a transformation matrix that relates displacements, velocities or forces at the contact points between the source (or the base) and the isolators.
Using the expressions above, we can find the velocities ( vsb ) and the forces (f sb ) in the source and in the base as a function of the external forces (q sb ) vsb = M1 v q sb (5) where The source-base mobility expressions due to internal forces (M sb1 ) and due to external forces (M sb2 ) will be calculated using a matrix approach [3].Additionally, the two isolators impedance matrix (Z i1 ) will be calculated using the Spectral Element Method (SEM).It is worth noting that using the SEM we will obtain a dynamic stiffness matrix (K i1 ) that is easily related with the impedance matrix, Z i1 = iωK i1 .
In this work, the vibration transmission from the source to the base will be quantified through the evaluation of the input power due to the source and the power transmitted to the base.The input power to the isolation system is given by the product of f s times the hermitian transpose of the velocity vs In addition, we can evaluate the transmitted power to the base which is given by the product of f b times the hermitian transpose of vb The difference between the input power P I and the transmitted power P T is the power dissipated in the isolator.This quantity allows us to quantify the performance of the isolation system.Notice also that both P I and P T are obtained from the velocities and loads at the source and the base.

The spectral element method (SEM)
The major advantage of the SEM [4] is that the elementary dynamic stiffness matrix is calculated from the analytical solution, this can be done by transforming the equation of motion to the frequency domain.Since the isolators are considered as a combination of a rod and a beam, in this section we will derive some basic relations of these structural elements.

Simple rod
The transformed homogeneous equation for a simple rod is where û is the transformed displacement, E is the Young modulus, A is the cross-section area, ρ is the material density and ω is the angular frequency.In this paper we will compare some models to describe the frequency-dependence of the Young modulus.These models will be briefly addressed in the following section.
In order to obtain the dynamic stiffness matrix, we will use the dynamic shape functions approach.These shape functions are essentially interpolation functions, but instead of being simple polynomials, they are the exact displacement distributions.The general longitudinal displacement for a rod is where k 1 = ω ρ E is the wavenumber and L is the element length.
The displacement end conditions are After solving for A and B, the longitudinal displacement can be written as where The member loads at each end of the rod are related to the displacements by where primes indicate differentiation with respect to x.These two equations can be written in the matrix form

Timoshenko beam
The transformed equations for a Timoshenko beam are [4] v and φ are the transformed displacement and rotation.K 1 is the shear correction factor and K 2 is a coefficient of the rotation inertia.The Mindlin coefficient 12 for shear correction is considered here, while setting K 2 to zero allows us to neglect the rotation inertia for analysis.However, the results presented in this work were obtained with K 2 = 1.G is the frequencydependent shear modulus and I is the second area moment of the cross section.The frequency-dependence of the shear modulus G will be treated in detail in the next section.
The general transverse displacement is where A, B, C and D are coefficients determined from the boundary conditions of the element and R i are the amplitude ratios The end conditions of the element are After solving a system of equations we get where Ĝ is a [4 × 4] matrix.The displacements and rotations can be written as where The shear force and the moment expressions are the primes indicating differentiation with respect to x.
Evaluation at x = 0 and We will consider V1 = − V (0), M1 = − M (0), V2 = V (L) and M2 = M (L).The relation between the loads and element degrees of freedom can be written as where We use the two stiffness matrices derived, coming from the rod k R and beam k B spectral elements, to obtain a global stiffness matrix for each isolator.After, the global stiffness matrices of the two isolators need to be assembled, leading to the dynamic stiffness matrix K i1 required for the evaluation of the impedance matrix Z i1 in Eqs ( 7) and ( 9).This process is detailed in Coronado [3] and is not repeated here for the sake of brevity.

Dissipation models
In this section, three viscoelastic models are presented, namely the Anelastic Displacement Fields (ADF), proposed by Lesieutre and Bianchini [6], the Fractional Calculus (FC) put forward by Bagley and Torvik [1] and the well-known complex modulus approach.These are needed to model the frequencydependence of the two viscoelastic isolators considered in this work.

Anelastic displacement fields (ADF)
The linear viscoelasticity theory [2] states the following law where σ is the stress tensor and ε is the strain tensor.In the special case of linear elastic materials, c(t) is constant, i.e. c(t) = c o .For viscoelastic materials, c(t) is not constant, so the integration must be calculated in the range [−∞, t].Considering that the material is initially at rest, i.e. ε j (t) = 0 for −∞ < t < 0, the j th component of the stress tensor is where G is the shear modulus.If we define h(t) = G(t) − G 0 , then The modulus after the material relaxation will be G 0 = lim t→∞ G(t).If the initial conditions are equal to zero, we observe that the last equation is the expression of the Laplace transform with h(s) = sh(s).Additionally, G 0 εj (s) accounts for the material elasticity and h(s)ε j (s) is related to dissipation.Because h(s) functions give the dissipative behavior, many authors have developed mathematical representation for these functions, such as [6] and [7].Lesieutre and Bianchini [6] represent h(s) as a series of functions in the Laplace-domain which may be written in the frequency-domain as G * (ω) is the complex frequency-dependent shear modulus, Ω i is the inverse of the characteristic time of relaxation at constant deformation, and i is the relaxation resistance.These series of two parameters and the static modulus G 0 will be calculated from experimental data making use of a non-linear least-squares optimization.Additional details and comparison between this model and some others can be found in [11].

Fractional calculus (FC)
In the standard linear viscoelastic model, a series of time derivatives applied to the stress fields is related to a series of time derivatives on the strain fields which is cumbersome to use.For viscoelastic materials having mechanical properties that are strongly frequency dependent, the number of time derivatives, M and N , in the series can be large.Consequently the number of empirical parameters in the model is large.Observations of the mechanical properties of viscoelastic materials indicated that the stress relaxation phenomenon appeared to be proportional to time raised to fractional powers [1].The general form of the present fractional derivative viscoelastic model is Experimental observations indicate that many viscoelastic materials can be modelled by retaining only the first fractional derivative term in each series [1].The result is a viscoelastic model with five parameters, b, G 0 , G 1 , α and β.
The fractional derivatives are defined by (32) where Γ is the gamma function.The fractional derivatives have the following property in the Fourier domain That is, the Fourier transform of the fractional derivative of order α of x(t) is (iω) α times the Fourier transform of x(t).Hence, taking the Fourier transform of the five parameter viscoelastic model, yields Factoring and dividing terms, one may obtain a relation between stress and strain, through a complex  frequency-dependent modulus G * (ω)

Complex modulus
This approach is motivated by observing the relationship between sinusoidal stress and sinusoidal strain in viscoelastic materials.The strain lags the stress and the imaginary part of the complex constant adequately describes this phenomenon.
The expression for the complex shear modulus can be written in terms of the frequency ω as where G (ω) is the elastic shear modulus and η G (ω) is  the loss factor which is defined by Finally, the Young modulus E (ω) can be easily calculated from the well known relation where ν is the Poisson ratio.Here, for simplicity, ν(ω) is considered frequency-independent so that the stiffness of the material is proportional to G * (ω).

Numerical examples
The properties of the isolation system (Fig. 1) were taken from the work of Gardonio et al. [5] and are  The frequency responses of vibration isolation systems have some interesting characteristics [8].For the above defined system we can recognize three distinct  areas (Fig. 7): below 10 Hz there are the 'suspension modes', where the source acts like a rigid body; between 10 and 50 Hz the 'isolation area' is located, where there are no resonances, and above 50 Hz there are the 'flexible modes', where we can find the base and isolator resonances.
For demonstration purposes, we will consider that the 'operational area' of a given machine coincides with the 'isolation area' defined above (between 10 and 50 Hz).However, the fitting of the damping material properties and the response of the isolation system are shown on a broader frequency range, for comparison purposes.Additionally, from Fig. 4 it is seen that the maximum value of the loss factor η associated with the transition frequency [9] occurs near 10 4 Hz, which is out of the 'isolation area' frequencies, so that the material can be considered as Type I [9].The fitting curves of the ADF and FC models were carried out using the Optimization Toolbox of the MATLAB .Specifically, we used the lsqnonlin function which uses a non-linear least squares algorithm.
In the Figs 2 and 3, we show the measured [10] and fitted curves of the shear modulus G and the loss factor η, in a 1-1000 Hz frequency range, for the ADF and FC models.The two models show a great fit with the measured values, although the ADF presents a slight difference in the 1000 Hz neighborhood.The ADF model owns 7 parameters (3 series) and the FC model owns 5 parameters, additionally the later one has in general a fastest convergence, but it depends hardly on the initial values of the optimization.Whereas, in this work we will carry simulations on the 1-1000 Hz range, we will show the results obtained for fitting process that have a wider range.
The Figs 4 and 5 show a broader frequency range, 0.1 − 10 5 Hz.In Fig. 4, we observe that the ADF model fits the measured curve in a acceptable manner in the low frequencies, but becomes poorer at the high extreme.On the other hand, in Fig. 5, we note that the FC model lacks some fitness in the low frequencies, but over 10 Hz it becomes highly suitable.As an additional comparison, it is shown in Fig. 6 an ADF fitting with 9 parameters (4 series).There is improvement mainly in the high frequencies, but in the medium range some error remains.
In the following figures we compare the input and transmitted power in the isolation system.The Figs 7 and 8 show a very similar behavior, with small differences in the high frequencies, near 1000 Hz.
The three latest Figs 9, 10 and 11 show very dissimilar responses obtained with the frequency-independent complex modulus model in comparison with the two former models.These differences can be noted already at the rigid source modes (near 10 Hz) and are more evident at the other frequencies.The use of a small constant loss factor (η = 0.1, Fig. 9) is intended to capture the system behavior in the lower frequencies.On the other hand, the use of a large constant loss factor (η = 0.7, Fig. 11) is intended to capture the system behavior in the higher frequencies.However, not only Figs 9 and 11 furnish wrong results for the higher and lower frequencies, respectively, but also both provide wrong results for the power dissipated in the isolation area.Moreover, Fig. 10 shows that even an interme-diate value of η = 0.3 does not capture the system behavior adequately.

Conclusions
The influence of the isolator stiffness and damping frequency dependence was analyzed, this was carried out on an isolation system with two isolators and a flexible base.In order to excite many vibration modes as an attempt to consider more realistic situations, two excitation forces, in y and z directions, and one excitation moment, in x direction, were considered.This fact allow us to see easily the different responses obtained using the approaches early discussed.
The isolation systems using Anelastic Displacement Fields (ADF) and Fractional Calculus (FC) models have shown very similar responses in the frequencies analyzed here.On the other hand the responses obtained using the frequency-independent complex modulus model were clearly different, which indicates that this model is not suitable when we work in a broad frequency range and with highly dissipative isolators.
Additionally, we can conclude that for a frequency analysis like the one done in the present work, the FC model has some advantage with respect to ADF, since the FC requires less parameters and its fitting converges quickly.