Finite Element Modeling and Experiments of Systems with Viscoelastic Materials for Vibration Attenuation

Abstract

This chapter addresses the finite element modeling methodologies intended for performance evaluation, analysis, and design of viscoelastic systems. The mathematical models widely used to represent the frequency and temperature dependent behavior of viscoelastic materials are also considered, namely the complex modulus approach, the fractional derivative model, the Golla-Hughes-McTavish (GHM) model, and the anelastic displacement fields (ADFs) model. The straightforward strategies to integrate the viscoelastic effects into finite element matrices of structural systems such as three-layer sandwich plates, intended for the modeling of medium and large-scale engineering structures, are presented. In the same context, emphasis is placed on the condensation methods for the reduction of the order of the finite element matrices to perform frequency-response functions, complex eigenvalue problem, and time domain analyses. Based on the fact that for viscoelastic materials subjected to dynamic loadings superimposed on static preloads, the classical modeling assuming isothermal conditions can lead to poor designs, since the energy dissipated within the volume of the material leads to temperature rises, an experimental investigation of the self-heating phenomenon is also addressed.

1. Introduction

It is well known that viscoelastic materials can be used with advantage to mitigate undesirable vibrations [1, 2] and, consequently, to increase the fatigue life of engineering structures to avoid catastrophes [3, 4]. As a result, they have been applied in a number of engineering applications such as robots, automobiles, airplanes, communication satellites, tall buildings, and space structures [5, 6].

However, one of the major difficulties regarding the use of viscoelastic materials for vibration mitigation is the fact that their mechanical properties vary strongly with environment and operational factors for which, the excitation frequency and temperature are usually considered to be the most relevant parameters [1]. This is a reason for which in the last decades, a number of viscoelastic models exist in the open literature to be used in conjunction with the finite element (FE) discretization procedure. Among those viscoelastic models, the fractional derivative model (FDM) is a more complex relationship between the stress and strain than the standard linear viscoelastic model [1], which is based on the use of fractional derivatives [7–9] in order to reduce the number of terms required by the generalized standard viscoelastic model. The so-called Golla-Hughes-McTavish (GHM) model initially proposed by Golla and Hughes [10] is based on the use of internal variables to represent the dissipation mechanism of viscoelastic materials. In the Laplace domain, the resulting GHM complex modulus function is a series of mini-oscillators terms similar to that of a damped single degree of freedom system. McTavish and Hughes [11] have demonstrated later the FE modeling strategy of a truss structure incorporating a structural viscoelastic damper by using the GHM model. Another viscoelastic model normally used is the so-called anelastic displacement field (ADF) model suggested by Lesieutre and his co-authors [12–14]. The strategy is the use of anelastic fields to represent the dissipative effects of viscoelastic materials similar to the GHM model. However, as opposed to the GHM, the ADF model can be formulated directly in the time domain.

For large FE models of real-world engineering structures incorporating viscoelastic materials typically composed by many thousands of degrees of freedom (DOFs), the inclusion of internal variables by using ADF or GHM models leads to global systems of equations of motion whose numbers of DOFs largely exceeds the order of the associated undamped system. As a result, if such evaluations are made based on response computations performed on the full finite element matrices, the computational time required to obtain the solutions is high [15–18]. It must be also reminded that if the interest is to perform active control techniques by using the resulting viscoelastic models, the dimension is a typical problem and must be not disregarded [19–23]. Thus, it is interesting to perform model condensation techniques especially adapted to viscoelastic systems with the aim of reducing the computational burden [24, 25].

The simplest model reduction method very useful to reduce static systems is the well-known Guyan static method [19], where the condensation is performed by separating the physical coordinates of the static system in master and slaves coordinates. However, as discussed in [20], this strategy is not adapted to viscoelastic systems [23], even when the reduced coordinates are a subset of the initial physical coordinates. In this case, the internal balancing method proposed by Moore [24] is more interesting to be used, but leads to reduced coordinates that are not necessarily a subset of the original physical coordinates. In order to circumvent this problem, Yae and Inman [19] have proposed a modified internal balancing method adapted to viscoelastic systems.

The other very simple reduction strategy, is the so-called modal reduction method [20], where in the construction of the reduction basis it is considered only the most relevant eigenvectors characterizing the dynamic behavior of the system in the frequency band of interest. However, it is shown here that since the viscoelastic stiffness matrix is frequency- and temperature-dependent, this strategy must be modified to generate a constant reduction basis formed by the static residues associated with the external loads and the viscoelastic damping forces. Clearly, the dimension of the reduced viscoelastic system of matrices corresponds to the number of eigenvectors and the static residues kept in the truncated series to form the reduction basis [18].

A natural extension of the modeling capability of engineering systems incorporating a viscoelastic damping device is the optimization aiming the reduction of cost and/or maximization of its performance. In the quest for optimization, the engineers are frequently faced with large FE models of real-world systems that require a great number of evaluations of the cost functions involved [19]. Consequently, if such computations are performed based on the full FE matrices, the time required to obtain the dynamic responses may be high. Here, a general strategy to construct a constant reduction basis for viscoelastic systems is suggested, composed by the eigenmodes of the associated conservative viscoelastic behavior enriched by the static residues associated with the external loadings and the viscoelastic damping forces. Also, the reduction basis can be easily extended to the case of robust condensation of viscoelastic systems subjected to parametric uncertainties if the static residues due to the small modifications are included into the basis. This robust condensation strategy in the frequency domain is a very interesting tool to be integrated into numerical optimization and/or model updating processes [25].

Another problem regarding the practical application of a viscoelastic damper to mitigate unwanted vibrations is the fact that the assumption of assuming a constant and uniform temperature distribution within the viscoelastic material can lead to a poor design and to unexpected damping performance of it due to the self-heating phenomenon [26–29]. As a result, it is expected that the temperature changes induced by the self-heating when the viscoelastic damper is subjected to cyclic excitations have a strong influence on the stiffness degradation and damping capacity of it. Moreover, in applications such as engine mounts, it must be considered the effects of the cyclic loadings superimposed on the static strains in the self-heating modeling procedure in which the prediction of the temperature evolution inside the viscoelastic material volume is an interesting thermoviscoelastic problem to be solved. Here, the influence of the cyclic loading superimposed on static preloads on the self-heating phenomenon is investigate numerically and experimentally for a three-dimensional translational viscoelastic mount used for vibration attenuation.

In the remainder, after the presentation of the theoretical foundations of the methodology, the description of some numerical studies of engineering systems incorporating passive constrained viscoelastic layers and discrete viscoelastic damping devices is addressed. The main interest is to illustrate those viscoelastic models and topics described in the methodology, intended to design and performance analyses of viscoelastically damped systems. In addition, the results of some experimental investigations with a freely suspended plate partially treated by passive constraining damping layer are carried out to validate the viscoelastic models and to confirm the effectiveness of the viscoelastic materials applied as a passive control strategy. Finally, an experimental investigation of the self-heating phenomenon on a three-dimensional translational viscoelastic mount subjected to dynamic loadings superimposed on static preloads is also addressed.

2. Review of viscoelastic models

2.1. The concept of complex modulus approach

According to the linear theory of viscoelasticity [30], the complex modulus of viscoelastic materials in the frequency domain is expressed as follows:

$G\left(\omega \right)=G^{\prime }\left(\omega \right)+i{G}^{″}\left(\omega \right)=G^{\prime }\left(\omega \right)\left[1+i\eta \left(\omega \right)\right]$E1

where G′(ω), G″(ω) and η(ω) = G″(ω)/G′(ω) designates the storage modulus, loss modulus and loss factor, respectively. As reported by [31], these parameters can be used to characterize the dynamic behavior of linear viscoelastic materials.

This model is adopted in the study reported herein since it enables the direct use of the data of viscoelastic materials commonly provided by the manufacturers [1]. Within this context, in the open literature, various mathematical models have been developed in order to represent the material modulus function, G(ω), as summarized later in this chapter.

However, it is widely recognized that the temperature also exerts a strong influence upon the properties of viscoelastic materials. Hence, it becomes important to account its variations in the modeling procedures of engineering systems incorporating viscoelastic materials. According to Nashif et al. [1] and Christensen [30], this can be done by making use of the frequency–temperature superposition principle (FTSP), where the damping properties of linear viscoelastic materials as functions of frequency at various temperatures can be collapsed on a single master curve, as illustrated in Figure 1, if appropriate horizontal shifts along the frequency axes are used. This establishes the well-known shift factor and reduced frequency concepts, symbolically expressed by the following relations [18]:

Figure 1.

Illustration of the frequency–temperature superposition principle.

$G\left({\omega }_r,T_0\right)=G\left({\alpha }_T\omega ,T_0\right),\eta \left({\omega }_r,T_0\right)=\eta \left({\alpha }_T\omega ,T_0\right)$E2

where T is the temperature, ω_r = α_T(T)ω is the reduced frequency, ω is the excitation frequency, α_T(T) is the shift factor, and T₀ is a reference value of temperature.

Functions G(ω_r) and α_T(T) can be obtained from experimental tests for specific viscoelastic materials [1]. Drake and Soovere [32] suggest analytical expressions for the complex modulus and shift factor for various commercially available viscoelastic materials, as functions of reduced frequency and temperature. The following equation represents the complex modulus for the 3M ISD112™ viscoelastic material [33], as provided by the authors:

$G\left({\omega }_r\right)=0.4307+\frac{1200}{1+3.24\times {\left(\frac{i{\omega }_r}{1543000}\right)}^{-0.18}+{\left(\frac{i{\omega }_r}{1543000}\right)}^{-0.6847}}\left[\mathrm{M}\mathrm{P}\mathrm{a}\right]$E3

where $\alpha \left(T\right)={10}^{\left(-3758.4\times \left(\frac{1}{T}-0.00345\right)-225.06\times log\left(0.00345\times T\right)+0.23273\times \left(T-290\right)\right)}$ and ω_r = α(T)ω.

Figure 2 shows the curves representing the variations of the storage modulus and loss factor as functions of the reduced frequency, as obtained from Eq. (3).

Figure 2.

Master curves and shift factor for the 3M™ ISD112 material.

2.2. Golla-Hughes-McTavish model (GHM)

According to the developments made by Golla and Hughes [10], and McTavish and Hughes [11], the viscoelastic material modulus function is expressed under the form:

$G\left(s\right)=G_0\left(1+{\displaystyle \sum _{i=1}^{N_G}{\alpha }_i\frac{s^2+2{\zeta }_i{\omega }_{i}s}{s^2+2{\zeta }_i{\omega }_{i}s+{\omega }_i^2}}\right)$E4

where G₀ is the static modulus, s designates the Laplace variable and $G_{\infty }=G_0\left(1+{\displaystyle {\sum }_{i=1}^{N_G}{\alpha }_i}\right)$ is the higher frequency modulus.

It can be clearly seen the similarity of each N_G term of the series appearing in Eq. (4) with the transfer function of a damped single DOF system known as mini-oscillator formed by the design variables (α_i, ω_i, ζ_i), identified for a specific viscoelastic material [20].

2.3. Anelastic displacement field model (ADF)

The ADF model, developed by Lesieutre and co-workers [12–14], is similar, in some aspects, to the GHM model. The modulus function is represented in Laplace domain by:

$G\left(s\right)=G_0\left(1+{\displaystyle \sum _{i=1}^{N_A}{\Delta }_i\frac{s}{{\omega }^2+{\Omega }_i^2}}\right)$E5

where G₀ is the low frequency modulus, N_A is the number of anelastic displacement fields, each of which is represented by parameters Ω_i, the inverse of the characteristic relaxation time at constant deformation, Δ_i is the relaxation magnitude, and $G_{\infty }=G_0\left(1+{\displaystyle {\sum }_{i=1}^{N_A}{\Delta }_i}\right)$ is the higher frequency modulus.

2.4. Fractional derivative model (FDM)

The FDM model proposed by Bagley and Torvik [7–9] can be viewed as the generalization of the standard viscoelastic model, by the introduction of non-integer time derivatives in the differential constitutive equation relating the stresses to strains as follows:

$\sigma \left(t\right)+{\displaystyle \sum _{m=1}^M{b}_m{D}^{{\beta }_m}}\sigma \left(t\right)=G_0\epsilon \left(t\right)+{\displaystyle \sum _{n=1}^P{G}_n{D}^{{\alpha }_n}\epsilon \left(t\right)}$E6

where σ_n and β_m are non-integer numbers. As stated by Bagley and Torvik [7], experimental observations have indicated that the behavior of many commercially viscoelastic materials can be adequately modeled by retaining only one term in each of the series appearing in Eq. (6). Through this simplified form, the modulus requires only five parameters, resulting in the following expression:

$G\left(s\right)=\left(G_0+G_{\infty }{s}^{\alpha }\right)/\left(1+b{s}^{\beta }\right)$E7

3. Curve fitting of material parameters

One important aspect regarding the use of GHM, ADF, and FDM models presented in previously section is the identification of the model parameters from experimental data. In most cases, manufacturers provide data sheets containing the material storage modulus G′(ω, T) and loss factor η(ω, T), defined in Eq. (1), as functions of excitation frequency ω and temperature T. Thus, it becomes necessary to express, for each model, the material modulus expressions as complex functions by making, s = iω :

$\left(\mathrm{G}\mathrm{H}\mathrm{M}\right):\mathrm{R}\mathrm{e}\left[G\left(\omega \right)\right]=G_0\left[1+{\displaystyle \sum _{i=1}^{N_G}{\alpha }_i\frac{\left[{\omega }^2\left({\omega }^2-{\omega }_i^2\right)+4{\zeta }_i^2{\omega }_i^2{\omega }^2\right]}{{\left({\omega }^2-{\omega }_i^2\right)}^2+4{\zeta }_i^2{\omega }_i^2{\omega }^2}}\right]$E8.a

$\left(\mathrm{G}\mathrm{H}\mathrm{M}\right):\mathrm{I}\mathrm{m}\left[G\left(\omega \right)\right]=G_0\left[{\displaystyle \sum _{i=1}^{N_G}{\alpha }_i\frac{2{\zeta }_i{\omega }_i^3\omega }{{\left({\omega }^2-{\omega }_i^2\right)}^2+4{\zeta }_i^2{\omega }_i^2{\omega }^2}}\right]$E8.b

$\left(\mathrm{A}\mathrm{D}\mathrm{F}\right):\mathrm{R}\mathrm{e}\left[G\left(\omega \right)\right]=G_0\left[1+{\displaystyle \sum _{i=1}^{N_A}{\Delta }_i\frac{{\left(\omega /{\Omega }_i\right)}^2}{1+{\left(\omega /{\Omega }_i\right)}^2}}\right]$E9.a

$\left(\mathrm{A}\mathrm{D}\mathrm{F}\right):\mathrm{I}\mathrm{m}\left[G\left(\omega \right)\right]=G_0\left[{\displaystyle \sum _{i=1}^{N_A}{\Delta }_i\frac{{\left(\omega /{\Omega }_i\right)}^2}{1+{\left(\omega /{\Omega }_i\right)}^2}}\right]$E9.b

$\left(\mathrm{F}\mathrm{D}\mathrm{M}\right):\mathrm{R}\mathrm{e}\left[G\left(\omega \right)\right]=\frac{G_0\left[1+b{\omega }^{\beta }cos\left(\pi \beta /2\right)\right]+G_{\infty }\left[{\omega }^{\alpha }cos\left(\pi \alpha /2\right)+b{\omega }^{\alpha +\beta }cos\left(\pi \left(\alpha -\beta \right)/2\right)\right]}{1+2b{\omega }^{\beta }cos\left(\pi \beta /2\right)+{\left(b{\omega }^{\beta }\right)}^2}$E10.a

$\left(\mathrm{F}\mathrm{D}\mathrm{M}\right):\mathrm{I}\mathrm{m}\left[G\left(\omega \right)\right]=\frac{G_{\infty }\left[{\omega }^{\alpha }sin\left(\pi \alpha /2\right)+b{\omega }^{\alpha +\beta }cos\left(\pi \left(\alpha -\beta \right)/2\right)\right]-G_{0}b{\omega }^{\beta }sin\left(\pi \beta /2\right)}{1+2b{\omega }^{\beta }cos\left(\pi \beta /2\right)+{\left(b{\omega }^{\beta }\right)}^2}$E10.b

From the equations above, for each viscoelastic model, the determination of the material parameters can be carried out by formulating a deterministic optimization problem in which the objective function represents the difference between the experimental data points and the corresponding model predictions. Clearly, the number of design variables depends on the previous choice of a model order, which is assumed to be sufficient to represent the frequency-dependent behavior in the frequency band of interest: $N_{par}^G=1+3N_G$ for the GHM; $N_{par}^A=1+2N_A$ for the ADF; and $N_{par}^F=5$ for the FDM.

4. Incorporation of the viscoelastic behavior into FE models

Based on the formulation presented in the previous sections, and including external viscous damping, the FE equations of motion in the frequency domain of a viscoelastic structure containing N DOFs can be expressed as [18]:

$\mathbf{M}\ddot{\mathbf{q}}\left(t\right)+\mathit{D}\stackrel{.}{\mathbf{q}}\left(t\right)+{\mathbf{K}}^{\ast }\mathbf{q}\left(t\right)=\mathbf{f}\left(t\right),\mathbf{f}\left(t\right)=\mathbf{b}\mathbf{u}\left(t\right)\mathrm{and}\mathbf{y}\left(t\right)=\mathbf{c}\mathbf{q}\left(t\right)$E11

where M, D, K^∗ ∈ R^N × N are, respectively, the mass, viscous damping, and complex stiffness matrices. q(t) ∈ R^N and f(t) ∈ R^Nare the vectors of displacements and external loads. y(t) ∈ R^c and u(t) ∈ R^f are the vectors of responses external loads. b ∈ R^N × f and c ∈ R^c × N are matrices, which enable to select, among the FE DOFs, those in which the external excitations are applied and responses are computed, respectively.

It should be emphasized that the dependency of the stiffness matrix on frequency and temperature is a consequence of the dependency of the material moduli with respect to these parameters as expressed by Eq. (1). Thus, by assuming harmonic excitation conditions, f(t) = Fe^iωt, u(t) = Ue^iωt, and steady-state harmonic responses, q(t) = Qe^iωt, y(t) = Ye^iωt, the following equations in the frequency domain are obtained:

$\left[\mathit{K}\left(\omega ,T\right)+i\omega \mathit{D}-{\omega }^2\mathbf{M}\right]\mathit{Q}\left(\omega ,T\right)=\mathit{F}\mathit{,}\mathit{F}=\mathit{b}\mathit{U}\mathrm{and}\mathit{Y}\left(\omega ,T\right)=\mathit{c}\mathit{Q}\left(\omega ,T\right)$E12

It is assumed that the model contains both elastic and viscoelastic elements. Thus, the elastic–viscoelastic correspondence principle [1] is applied leading to:

$\mathit{K}\left(\omega ,T\right)={\mathit{K}}_e+{\mathit{K}}_v\left(\omega ,T\right)$E13

where K_e and K_v(ω, T) are the stiffness matrices associated with the elastic and viscoelastic substructures, respectively. Also, taking into account the stress state for the specific viscoelastic element assumed in the analysis, one of the moduli (through the relation, G(ω, T) = E(ω, T)/2(1 + ν)) can be factored out of, ${\mathit{K}}_v\left(\omega ,T\right)=G\left(\omega ,T\right){\overline{\mathit{K}}}_v$, where ${\overline{\mathit{K}}}_v$ is the frequency- and temperature-independent matrix. Then, Eqs. (12) and (13) can be combined to give the following complex receptance matrix for the viscoelastic system:

$\mathit{H}\left(\omega ,T\right)=\mathit{c}\mathit{Z}{\left(\omega ,T\right)}^{\mathit{‐}1}\mathit{b}$E14

where $\mathit{Z}\left(\omega ,T\right)={\mathit{K}}_e+G\left(\omega ,T\right){\overline{\mathit{K}}}_v+i\omega \mathit{D}-{\omega }^2\mathit{M}$ is the complex dynamic stiffness matrix.

Clearly, the difficulty in predicting the dynamic responses for viscoelastic systems comes from the fact that the stiffness matrix depends on frequency and temperature. As a result, one has an eigenvalue problem that must be solved iteratively [34]. Some other procedures for dealing with this problem have been suggested, such as the contributions by Palmeri and Ricciardelli [35] and by Palmeri et al. [36], where the eigenvalue problem of viscoelastic systems has been derived in the time domain with the help of the novel concept of modal relaxation functions. In the papers proposed by Yuan and Agrawal [37] and Wagner and Adhikari [38], an alternative state-space approach has been proposed for the time-domain analysis of viscoelastic structures. Others alternatives have been suggested based on the adoption of particular representations for the frequency-dependent behavior of the viscoelastic materials [39]. Such an approach is used in the FDM, GHM, and ADF models, which enable to transform the equations of motion of a viscoelastic system in the time-domain into state-space forms, with frequency-independent state matrices, at the expense of a typically high increase in the order of the system matrices.

4.1. GHM model

Applying the Laplace Transformation to Eq. (11) and replacing G(s) by the modulus function given by Eq. (4), the following governing equation of motion is obtained:

$\left[s^2\mathit{M}+s\mathit{D}+{\mathit{K}}_e+G_0{\overline{\mathit{K}}}_v\left(1+{\displaystyle \sum _{i=1}^{N_G}{\alpha }_i\frac{s^2+2{\zeta }_i{\omega }_{i}s}{s^2+2{\zeta }_i{\omega }_{i}s+{\omega }_i^2}}\right)\right]\mathit{Q}\left(\omega ,T\right)=\mathit{F}$E15

A series of dissipation coordinates ${\mathit{Q}}_i^G\left(i=1,\dots ,N_G\right)$ are defined as:

${\mathit{Q}}_i^G=\frac{{\omega }_i^2}{s^2+2{\zeta }_i{\omega }_{i}s+{\omega }_i^2}\mathit{Q}\left(\omega ,T\right)$E16

Upon introduction of Eq. (16) into (15), after some mathematical manipulations and back transformation to time domain, the following coupled system of equations is obtained:

${\mathit{M}}_G{\ddot{\mathit{q}}}_G\left(t\right)+{\mathit{D}}_G{\stackrel{.}{\mathit{q}}}_G\left(t\right)+{\mathit{K}}_G{\mathit{q}}_G\left(t\right)={\mathit{f}}_G\left(t\right)$E17

where

${\mathit{M}}_G=\left[\begin{array}{cccc}\hfill \mathit{M}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \frac{{\alpha }_1}{{\omega }_1^2}{\mathit{K}}_v^0\hfill & \hfill \mathit{0}\hfill & \hfill ⋮\hfill \\ \hfill ⋮\hfill & \hfill \mathit{0}\hfill & \hfill \ddots \hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \frac{{\alpha }_{N_G}}{{\omega }_{N_G}^2}{\mathit{K}}_v^0\hfill \end{array}\right],{\mathit{D}}_G=\left[\begin{array}{cccc}\hfill \mathit{D}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \frac{2{\zeta }_1{\alpha }_1}{{\omega }_1}{\mathit{K}}_v^0\hfill & \hfill \mathit{0}\hfill & \hfill ⋮\hfill \\ \hfill ⋮\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \frac{2{\zeta }_{N_G}{\alpha }_{N_G}}{{\omega }_{N_G}}{\mathit{K}}_v^0\hfill \end{array}\right]\text{,}$E18

${\mathit{K}}_G=\left[\begin{array}{cccc}\hfill {\mathit{K}}_e+{\mathit{K}}_v^{\infty }\hfill & \hfill -{\alpha }_1{\mathit{K}}_v^0\hfill & \hfill \cdots \hfill & \hfill -{\alpha }_{N_G}{\mathit{K}}_v^0\hfill \\ \hfill -{\alpha }_{N_G}{{\mathit{K}}_v^0}^T\hfill & \hfill {\alpha }_1{\mathit{K}}_v^0\hfill & \hfill \mathit{0}\hfill & \hfill ⋮\hfill \\ \hfill ⋮\hfill & \hfill \mathit{0}\hfill & \hfill \ddots \hfill & \hfill \mathit{0}\hfill \\ \hfill -{\alpha }_{N_G}{{\mathit{K}}_v^0}^T\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill {\alpha }_{N_G}{\mathit{K}}_v^0\hfill \end{array}\right],{\mathit{q}}_G={\left[\mathit{q},{\mathit{q}}_1^G,\dots ,{\mathit{q}}_{N_G}^G\right]}^T\mathrm{and}{\mathit{f}}_G={\left[\mathit{f},\mathit{0},\dots ,\mathit{0}\right]}^T\text{.}$E19

where ${\mathit{M}}_G,{\mathit{D}}_G,{\mathit{K}}_G\in R^{T_{G}x{T}_G}$ with T_G = N(1 + N_G), ${\mathit{K}}_v^0=G_0{\overline{\mathit{K}}}_v$ and ${\mathit{K}}_v^{\infty }=G_{\infty }^G{\overline{\mathit{K}}}_v\text{.}$.

4.2. ADF model

Following the procedure outlined above for the GHM model, the equations of motion obtained by considering the ADF model is expressed by Eq. (5) [12]:

$\left[s^2\mathit{M}+s\mathit{D}+{\mathit{K}}_e+G_0{\overline{\mathit{K}}}_v\left(1+{\displaystyle \sum _{i=1}^{N_A}{\Delta }_i}\frac{s}{{\omega }^2+{\Omega }_i^2}\right)\right]\mathit{Q}\left(\omega ,T\right)=\mathit{F}$E18

In the ADF model the coordinates are separated into an elastic part, Q^e, which is instantaneously proportional to the stress, and an anelastic part, ${\mathit{Q}}_i^A\left(i=1,\dots ,N_A\right)\text{,}$ that represents the dissipation effects of the viscoelastic materials [13]:

$\mathit{Q}\left(\omega ,T\right)={\mathit{Q}}^e+{\displaystyle \sum _{i=1}^{N_A}{\mathit{Q}}_i^A}$E19

The system of equations is adopted for describing the time-evolution of the anelastic fields:

$\frac{C_i}{{\Omega }_i}{\mathit{K}}_v^{\infty }{\stackrel{.}{\mathit{Q}}}_i^A+C_i{\mathit{K}}_v^{\infty }{\mathit{Q}}_i^A-{\mathit{K}}_v^{\infty }\mathit{Q}=\mathit{0}$E20

Introducing Eq. (19) into (18), transforming the resulting equation to time domain and combining it with Eq. (20), the following coupled system of equations is obtained:

${\mathit{M}}_A{\ddot{\mathit{q}}}_A\left(t\right)+{\mathit{D}}_A{\stackrel{.}{\mathit{q}}}_A\left(t\right)+{\mathit{K}}_A{\mathit{q}}_A\left(t\right)={\mathit{f}}_A\left(t\right)$E21

where:

$\begin{array}{c}{\mathit{M}}_A=\left[\begin{array}{ccc}\hfill \mathit{M}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \end{array}\right],{\mathit{K}}_A=\left[\begin{array}{cccc}\hfill {\mathit{K}}_e+{\mathit{K}}_v^{\infty }\hfill & \hfill -{\mathit{K}}_v^{\infty }\hfill & \hfill \cdots \hfill & \hfill -{\mathit{K}}_v^{\infty }\hfill \\ \hfill -{{\mathit{K}}_v^{\infty }}^T\hfill & \hfill C_1{\mathit{K}}_v^{\infty }\hfill & \hfill \mathit{0}\hfill & \hfill ⋮\hfill \\ \hfill ⋮\hfill & \hfill \mathit{0}\hfill & \hfill \ddots \hfill & \hfill \mathit{0}\hfill \\ \hfill -{{\mathit{K}}_v^{\infty }}^T\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill C_{N_A}{\mathit{K}}_v^{\infty }\hfill \end{array}\right],\\ {\mathit{D}}_A=\left[\begin{array}{cccc}\hfill \mathit{D}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \frac{C_1}{{\Omega }_1}{\mathit{K}}_v^{\infty }\hfill & \hfill \mathit{0}\hfill & \hfill ⋮\hfill \\ \hfill ⋮\hfill & \hfill \mathit{0}\hfill & \hfill \ddots \hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \frac{C_{N_A}}{{\Omega }_{N_A}}{\mathit{K}}_v^{\infty }\hfill \end{array}\right],{\mathit{q}}_A={\left[\mathit{q},{\mathit{q}}_1^A,\dots ,{\mathit{q}}_{N_G}^A\right]}^T,{\mathit{f}}_A={\left[\mathit{f},\mathit{0},\dots ,\mathit{0}\right]}^T\text{.}\end{array}$E27

where, ${\mathit{M}}_A,{\mathit{D}}_A,{\mathit{K}}_A\in R^{T_{A}x{T}_A}$ with $T_A=N\left(1+N_A\right),{\mathit{K}}_v^0=G_0{\overline{\mathit{K}}}_v$ and ${\mathit{K}}_v^{\infty }=G_{\infty }^A{\overline{\mathit{K}}}_v$.

4.3. FDM model

Introducing Eq. (7) into the Laplace transform of Eq. (11) and multiplying the resulting equation by, (1 + bs^β), the following system of equations of motion is obtained:

$\left(s^{2+\beta }b\mathit{M}+s^{1+\beta }b\mathit{D}+s^2\mathit{M}+s\mathit{D}+s^{\beta }b{\mathit{K}}_e+s^{\alpha }{G}_{\infty }{\overline{\mathit{K}}}_v+{\mathit{K}}_e+G_0{\overline{\mathit{K}}}_v\right)\mathit{Q}=\left(1+b{s}^{\beta }\right)\mathit{F}$E22

According to Bagley and Torvik [8], the system of Eq. (22) can be written under the compact form, ${\displaystyle {\sum }_{j=0}^J{s}^{\raisebox{1ex}{$j$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}{\mathit{A}}_j\mathit{Q}=\left(1+b{s}^{\beta }\right)\mathit{F}}$, where J = m(2 + β), m is the smallest common multiple of the denominators of the fractional terms α and β, matrices A_j are computed by direct comparison between Eq. (22) and the compact representation:

$\left(s^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}{\mathit{M}}_F+{\mathit{K}}_F\right){\mathit{Q}}_F={\mathit{F}}_F$E23

where

$\begin{array}{c}{\mathit{M}}_F=\left[\begin{array}{ccccc}\hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill {\mathit{A}}_j\hfill \\ \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill {\mathit{A}}_j\hfill & \hfill {\mathit{A}}_{j-1}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋰\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill \mathit{0}\hfill & \hfill {\mathit{A}}_j\hfill & \hfill \cdots \hfill & \hfill {\mathit{A}}_3\hfill & \hfill {\mathit{A}}_2\hfill \\ \hfill {\mathit{A}}_j\hfill & \hfill {\mathit{A}}_{j-1}\hfill & \hfill \cdots \hfill & \hfill {\mathit{A}}_2\hfill & \hfill {\mathit{A}}_1\hfill \end{array}\right],{\mathit{K}}_F=\left[\begin{array}{ccccc}\hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill -{\mathit{A}}_j\hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill -{\mathit{A}}_{j-1}\hfill & \hfill \mathit{0}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋰\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill -{\mathit{A}}_j\hfill & \hfill -{\mathit{A}}_{j-1}\hfill & \hfill \cdots \hfill & \hfill -{\mathit{A}}_2\hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill A_0\hfill \end{array}\right],\\ {\mathit{Q}}_F={\left[s^{\left(j-1\right)/m}\mathit{Q},s^{\left(j-2\right)/m}\mathit{Q},\dots ,s^{1/m}\mathit{Q},\mathit{Q}\right]}^T,{\mathit{F}}_F={\left[\mathit{0},\dots ,\left(1+b{s}^{\beta }\right)\mathit{F}\right]}^T,{}_F,{\mathit{K}}_F\in R^{T_{F}x{T}_F},T_F=JN\text{.}\end{array}$E30

5. Modal reduction methods for viscoelastic systems

It can be seen that the internal non-physical coordinates used by the GHM and ADF models to represent the viscoelastic dissipation mechanism lead to large system of equations of motion. Thus, it requires a high computational cost in the computation of dynamic responses of the viscoelastic system [18]. To develop the formulation pertaining the modal reduction method, it is convenient to transform the Eqs. (17) and (21) into an equivalent first-order form (space-state model) with an output equation added as follows:

$\stackrel{.}{\mathit{x}}=\mathit{Ax}+\mathit{B}\mathit{u}\mathit{,}\mathit{y}=\mathit{C}\mathit{x}$E24

where A assumes the following forms for the GHM and ADF models, respectively:

$\begin{array}{c}\mathit{A}=\left[\begin{array}{cccccccc}\hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \mathit{I}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \mathit{I}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \cdots \hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \cdots \hfill & \hfill ⋮\hfill \\ \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{I}\hfill \\ \hfill -{\mathit{M}}^{-1}{\mathit{K}}_T\hfill & \hfill {\alpha }_1{\mathit{M}}^{-1}\mathit{R}\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{N_G}{\mathit{M}}^{-1}\mathit{R}\hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill {\omega }_1^2{\mathit{\Lambda }}^{-1}{\mathit{R}}^T\hfill & \hfill {\omega }_1^2\mathit{I}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill -2{\zeta }_1{\omega }_1\mathit{I}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill {\omega }_{N_G}^2{\mathit{\Lambda }}^{-1}{\mathit{R}}^T\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill {\omega }_{N_G}^2\mathit{I}\hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill -2{\zeta }_{N_G}{\omega }_{N_G}\mathit{I}\hfill \end{array}\right],\mathit{x}=\left\{\begin{array}{c}\hfill {\mathit{q}}^T\hfill \\ \hfill {{\mathit{q}}_1^{GG}}^T\hfill \\ \hfill ⋮\hfill \\ \hfill {{\mathit{q}}_{N_G}^{GG}}^T\hfill \\ \hfill {\stackrel{.}{\mathit{q}}}^T\hfill \\ \hfill ⋮\hfill \\ \hfill {{\stackrel{.}{\mathit{q}}}_1^{GG}}^T\hfill \\ \hfill {{\stackrel{.}{\mathit{q}}}_{N_G}^{GG}}^T\hfill \end{array}\right\}\\ \mathit{A}=\left[\begin{array}{ccccc}\hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill \mathit{0}\hfill & \hfill \mathit{I}\hfill \\ \hfill \frac{{\Omega }_1}{C_1}{\mathit{\Lambda }}^{-1}{\mathit{R}}^T\hfill & \hfill -{\Omega }_1\mathit{I}\hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill & \hfill \mathit{0}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill \frac{{\Omega }_{N_A}}{C_{N_A}}{\mathit{\Lambda }}^{-1}{\mathit{R}}^T\hfill & \hfill \mathit{0}\hfill & \hfill \cdots \hfill & \hfill -{\Omega }_{N_A}\mathit{I}\hfill & \hfill \mathit{0}\hfill \\ \hfill -{\mathit{M}}^{-1}{\mathit{K}}_T\hfill & \hfill {\mathit{M}}^{-1}\mathit{R}\hfill & \hfill \cdots \hfill & \hfill {\mathit{M}}^{-1}\mathit{R}\hfill & \hfill \mathit{0}\hfill \end{array}\right],\mathit{x}=\left\{\begin{array}{c}\hfill {\mathit{q}}^T\hfill \\ \hfill {{\mathit{q}}_1^{\mathit{AA}}}^T\hfill \\ \hfill {{\mathit{q}}_{N_A}^{\mathit{AA}}}^T\hfill \\ \hfill ⋮\hfill \\ \hfill {\stackrel{.}{\mathit{q}}}^T\hfill \end{array}\right\}\mathrm{and}{\mathit{K}}_T={\mathit{K}}_e+{\mathit{K}}_v^{\infty }\text{.}\end{array}$E32

Since A is a non-symmetric matrix, the following eigenvalue problems are formulated:

$\mathit{A}{\mathit{R}}_r=\mathit{\Lambda }{\mathit{R}}_r,{\mathit{A}}^T{\mathit{R}}_l=\mathit{\Lambda }{\mathit{R}}_l$E25

where Λ is the spectral matrix composed by the complex eigenvalues, and R_r and R_l are the modal matrices whose columns contain the right- and left-hand-side eigenvectors, respectively, normalized by the relation, ${\mathit{R}}_l^T{\mathit{R}}_r=\mathit{{\rm I}}$. Thus, the left- and right-hand-side eigenvectors can be separated into a structural and dissipative eigenvectors represented, respectively, by the subscripts e and d according to the following expressions:

${\mathit{R}}_r=\left[{\mathit{R}}_{re},{\mathit{R}}_{\mathit{rd}}\right],{\mathit{R}}_l=\left[{\mathit{R}}_{le},{\mathit{R}}_{ld}\right]$E26

The eigenvectors related to the internal variables are overdamped with a small contribution to the dynamic behavior. As result, the state of the system can be approached by the contributions of the elastic eigenvectors, x ≈ R_rex_e, so that the Eq. (24) can be expressed under the following form [24],

and y = C_rx_e, where ${\mathit{A}}_r={\mathit{R}}_{le}^T\mathit{A}{\mathit{R}}_{re}$, ${\mathit{B}}_r={\mathit{R}}_{le}^T\mathit{B}$ and C_r = CR_re are input and output state reduced matrices of the system.

This modal reduction technique is very simple to implement since the choice of eigenmodes to be retained is based only on the frequency band of interest. However, for some viscoelastic systems, where the elastic modes may be overdamped, care must be taken, and an eigenfrequency analysis a priori must be performed in order to identify these eigenmodes. Also, the main disadvantage of the reduced state-space system, intended to apply control techniques, is that the matrices are complex. However, since all overdamped (relaxation/dissipative) modes were neglected, all elements of

are composed of complex conjugates pairs, such that they can be rewritten as follows:

${\mathit{A}}_r=\left[\begin{array}{cccc}\hfill \ddots \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill {\lambda }_i\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill {\overline{\lambda }}_i\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{array}\right],{\mathit{B}}_r=\left[\begin{array}{c}\hfill ⋮\hfill \\ \hfill {\varphi }_i\hfill \\ \hfill {\overline{\varphi }}_i\hfill \\ \hfill ⋮\hfill \end{array}\right],{\mathit{C}}_r=\left[\begin{array}{cccc}\hfill \cdots \hfill & \hfill {\psi }_i\hfill & \hfill {\overline{\psi }}_i\hfill & \hfill \cdots \hfill \end{array}\right]$E27

where λ_i and ${\overline{\lambda }}_i$ are the retained elastic eigenvalues and their complex conjugates.

According to Friot and Bouc [40], to construct a real representation of the state-space system represented by Eq. (24), one can use a state transformation matrix such as ${\mathit{x}}_e=\mathit{T}{\overline{\mathit{x}}}_e$, where T is defined as:

$\mathit{T}=\left[\begin{array}{cccc}\hfill \ddots \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill -\frac{1}{2}\frac{j}{Im\left({\lambda }_i\right)}\hfill & \hfill \frac{1}{2}\frac{j}{Im\left({\lambda }_i\right)}\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \frac{1}{2}-\frac{1}{2}\frac{jRe\left({\lambda }_i\right)}{Im\left({\lambda }_i\right)}\hfill & \hfill \frac{1}{2}+\frac{1}{2}\frac{jRe\left({\lambda }_i\right)}{Im\left({\lambda }_i\right)}\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{array}\right]$E28

Thus, the real state-space system can be constructed,

and $\mathit{y}={\overline{\mathit{C}}}_r{\overline{\mathit{x}}}_e$, where Ā_r = TA_rT^− 1, ${\overline{\mathit{B}}}_r=\mathit{T}{\mathit{B}}_r$, and ${\overline{\mathit{C}}}_r={\mathit{C}}_r{\mathit{T}}^{-1}$.

6. Internal balancing method

The internal balancing method is another interesting reduction method of viscoelastic systems. However, it does not guarantee that the reduced coordinates are a subset of the original coordinates of the system. This method is based on the controllability and observability of each balanced state.

Based on the work by Moore [24], the balanced internal system is defined such that the grammians are diagonal and equal. For a viscoelastic system expressed by Eq. (24), the controllability and observability grammians, denoted by W_c and W₀, are defined in order to satisfy the Lyapunov stability equations, AW_c + W_cA^T = − BB^T and A^TW₀ + W₀A = − C^TC, where the Cholesky decomposition of matrix W_c is performed by the relation, ${\mathit{W}}_c={\mathit{L}}_c{\mathit{L}}_c^T$. Also, a transformation matrix P is introduced, P = L_cUΛ^− 1/2, where Λ and U are the eigenvalues and eigenvectors of the eigenvalue problem, ${\mathit{L}}_c{\mathit{W}}_0{\mathit{L}}_c^T$, in such a way that the internal balanced model is given as:

${\stackrel{.}{\mathit{x}}}_r={\mathit{A}}_r{\mathit{x}}_r+{\mathit{B}}_r\mathit{u}\mathit{,}\mathit{y}={\mathit{C}}_r{\mathit{x}}_r$E29

where A_r = P^− 1AP, B_r = P^− 1B, C_r = CP, and x_r = P^− 1x.

According to the definition of an internally balanced system, W_c(P) = W₀(P) = diag(σ_i), where σ_i(1, 2, …, N) is the controllability of a state i, and N is the number of DOFs. Thus, according to the high and small values of states, σ_i, the internally balanced system (29) can be partitioned into retained states, x_rr, and reduced states, x_rd, as follows:

$\left[\begin{array}{c}\hfill {\stackrel{.}{\mathit{x}}}_{rr}\hfill \\ \hfill {\stackrel{.}{\mathit{x}}}_{\mathit{rd}}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill {\mathit{A}}_{rr}\hfill & \hfill {\mathit{A}}_{\mathit{rd}}\hfill \\ \hfill {\mathit{A}}_{dr}\hfill & \hfill {\mathit{A}}_{dd}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\mathit{x}}_{rr}\hfill \\ \hfill {\mathit{x}}_{\mathit{rd}}\hfill \end{array}\right]+\left[\begin{array}{c}\hfill {\mathit{B}}_{rr}\hfill \\ \hfill {\mathit{B}}_{dd}\hfill \end{array}\right]\mathit{u}$E30

Hence, the undesirable states, x_rd, must be removed by performing a simple static reduction and retaining only the contributions of states, x_rr.

7. Robust condensation strategy of viscoelastic systems

For real-world systems incorporating viscoelastic elements, it is practically impossible to perform time and frequency analyses using the GHM, ADF or FDM models, owing to the prohibitive computation times and storage memory required to evaluate the augmented equations of motion. This fact motivates the use of alternatives strategies with the aim of diminishing the model dimensions (and the computational burden, as a result), while keeping a reasonable predictive capacity of the numerical models. This can be done based on the assumption that the exact responses, given by the resolution of Eq. (12), can be approached by solutions in a reduced space as follows [18]:

$\mathit{Q}=\mathit{T}{\mathit{Q}}_r$E31

where T ∈ C^N × NR is the vector basis (or Ritz basis) and Q_r ∈ C^NR where NR ≪ N is the number of vectors in the basis. Hence, the transfer function can be rewritten as follows:

${\mathit{H}}_r\left(\omega ,T\right)=\mathit{c}{\mathit{Z}}_r{\left(\omega ,T\right)}^{\mathit{‐}\mathit{1}}\mathit{b}$E32

where ${\mathit{Z}}_r\left(\omega ,T\right)={\mathit{T}}^T{\mathit{K}}_e\mathit{T}+G\left(\omega ,T\right){\mathit{T}}^T{\overline{\mathit{K}}}_v\mathit{T}-{\omega }^2{\mathit{T}}^T\mathit{M}\mathit{T}$ is reduced dynamic stiffness matrix.

It can be seen that for viscoelastic systems, the construction of the basis is an issue, since the stiffness matrix is frequency- and temperature-dependent. Thus, three solutions are possible: (a) one can neglect this dependence by considering the stiffness matrix as independent from frequency and temperature [15, 19]. In this case, the basis is also constant; (b) one can use a reduction basis obtained by the resolution of the nonlinear eigenvalue problem [22, 23]; (c) it is possible to use an iterative method, which allows the re-actualization of the basis according to frequency [15, 16].

The strategy adopted here consists in using a constant reduction basis computed by using the tangent stiffness matrix representing the static behavior of the viscoelastic material. As can be seen in Figure 2, on the low frequency range, by prolonging the modulus and loss factor curves by asymptotes, the extrapolation leads to the real values G₀ and η₀ = 0. The tangent stiffness matrix can thus be calculated by the relation, ${\mathit{K}}_{\mathit{0}}={\mathit{K}}_e+G_0{\overline{\mathit{K}}}_v$ [15].

The nominal basis containing the first retained modes of the viscoelastic system can thus be obtained by the resolution of the following eigenvalue problem [18]:

$\begin{array}{ll}\left({\mathit{K}}_0-{\lambda }_i\mathit{M}\right){\varphi }_i=\mathit{0}\mathit{,}\hfill & i=1,\dots ,N\hfill \\ {\varphi }_0=\left[\begin{array}{cccc}\hfill {\varphi }_1\hfill & \hfill {\varphi }_2\hfill & \hfill \dots \hfill & \hfill {\varphi }_{NR}\hfill \end{array}\right],\hfill & {\mathit{\Lambda }}_0=\mathit{diag}\left({\lambda }_1,\dots ,{\lambda }_{NR}\right)\hfill \end{array}$E33

where ϕ₀ contains the mode shapes of the conservative associated system, which is further enriched by introducing the following residues associated with the external loads and the viscoelastic damping forces, respectively:

$\mathit{R}={\mathit{K}}_0^{\mathit{‐}1}\mathit{b}\mathit{,}{\mathit{R}}_v^{\mathit{0}}={\mathit{K}}_0^{\mathit{‐}1}{\overline{\mathit{K}}}_v{\varphi }_0$E34

Thus, the enriched basis of reduction is represented as follows:

${\mathit{T}}_0=\left[\begin{array}{ccc}\hfill {\varphi }_0\hfill & \hfill \mathit{R}\hfill & \hfill {\mathit{R}}_v^{\mathit{0}}\hfill \end{array}\right]$E35

The basis (35) can be used to reduce viscoelastic systems with reasonable accuracy, but it is not “robust” enough to taking into account parametric uncertainties and to be used in conjunction with optimization and/or model updating processes. This means that the basis given by Eq. (35) should in principle be updated successively to guarantee a satisfactory accuracy of model reduction during iterations. However, according to the authors in reference [18], the strategy consists in using a fixed reduction basis evaluated from the initial model [as given by Eq. (35)], to be further enriched by a set of residual vectors depending on the parametric modifications introduced by the viscoelastic treatment. In most practical applications, the viscoelastic surface treatments are not applied to the entire structure, but only in specific zones. Thus, based on Eq. (12), one can write the following modified system:

$\left[\mathit{\Delta }\mathit{Z}\left(\omega ,T\right)+\mathit{Z}\left(\omega ,T\right)\right]\mathit{Q}=\mathit{F}$E36

where ΔZ(ω, T) = ΔΚ_v(ω, T) − ω²ΔΜ_v is the dynamic stiffness matrix due the perturbations.

Hence, Eq. (36) can be interpreted as the equilibrium equation of the nominal model, subjected to forces of modifications, Z(ω, T)Q = F + f_Δ(ω, T), where f_Δ(ω, T) = − ΔZ(ω, T)Q, $\mathit{\Delta }{\mathit{{\rm M}}}_v={\displaystyle \sum _{i=1}^{n\_mp}\Delta m_i{\mathit{{\rm M}}}_v}$ and $\mathit{\Delta }{\overline{\mathit{{\rm K}}}}_v={\displaystyle \sum _{i=1}^{n\_mp}\Delta k_i{\overline{\mathit{{\rm K}}}}_v}$. Δm_i and Δk_i are the mass and stiffness variations, ${\mathit{{\rm M}}}_v={\displaystyle \underset{i}{\cup }{\mathit{{\rm M}}}_{v_i}}$ and ${\overline{\mathit{{\rm K}}}}_v={\displaystyle \underset{i}{\cup }{\overline{\mathit{{\rm K}}}}_{v_i}}$ are the mass and stiffness matrices of the zones, ${\mathit{{\rm M}}}_{v_i}$ and ${\overline{\mathit{{\rm K}}}}_{v_i}$, are the elementary viscoelastic mass and stiffness matrices, and ΔΜ_v and $\mathit{\Delta }{\overline{\mathit{{\rm K}}}}_v$ are the matrices to be reduced, which are in general sparse and nonlinear functions of the design parameters.

7.1. Basis of displacements associated with the structural modifications

The vector of forces of modifications, f_Δ(ω, T), depends on the response of the modified system, Q. Since this response is unknown a priori, the forces associated with the modifications cannot be computed exactly. The strategy is to generate the forces due to the small modifications by f_Δ(ω, T) using a truncated basis of normal modes of the mean model according to Eq. (40); next, the static responses of the mean model can be obtained by using the estimated forces of modifications. These two steps must be performed for each design variable subjected to small modifications.

${\mathit{f}}_{\Delta }\left(\omega ,T\right)\approx -\mathit{\Delta }\mathit{Z}\left(\omega ,T\right){\varphi }_0{\mathit{Q}}_r$E37

Hence, after obtaining the basis of forces, one can calculate a series of static responses of the system based on the tangent stiffness matrix according to the following form:

${\mathit{R}}_{\Delta }={\mathit{K}}_0^{\mathit{‐}1}{\mathit{f}}_{\Delta }$E38

and the final robust basis of reduction taking into account the small modifications is as follows:

$\mathit{T}=\left[\begin{array}{cc}\hfill {\mathit{T}}_0\hfill & \hfill {\mathit{R}}_{\Delta }\hfill \end{array}\right]$E39

Figure 3 illustrates a cycle of optimization process by using the robust basis, where the standard Ritz basis is increased by the static residues associated with the external loadings and the forces associated with the viscoelastic modifications. This procedure is used to approximate the behavior of the modified viscoelastic system without the re-actualization of the nominal basis, leading to a drastic reduction of the time required for computing the response of the large-scale viscoelastic systems.

Figure 3.

Block diagram showing the standard and the robust optimization procedures.

8. Review of FE modeling of passive constraining layer damping

One type of structure of particular interest in terms of practical viscoelastic applications is the three-layer sandwich plate illustrated in Figure 4. In the present work, the FE modeling procedure is summarized based on the original contribution made by Khatua and Cheung [41] and implemented by de Lima et al. [42]. The sandwich element is composed by four nodes and seven DOFs per node, as depicted in the figure, where u and v are the in-plane displacements in directions x and y, respectively, w is the transverse displacement, and θ_x and θ_y are the rotations.

Figure 4.

Illustration of the three-layer sandwich plate FE.

In the development of the theory, it has been assumed that the three layers are perfectly bounded and the materials involved are considered to be isotropic with linear behavior. These assumptions are reasonable [18], since, in practice, the most commercially available viscoelastic materials for vibration attenuation are self-adhesive. The analysis also assumes the Kirchhoff’s theory for the base plate and constraining layer with the same rotations, and only for the viscoelastic core, the transverse shear is included (Mindlin’s theory). The transverse displacement is assumed to be same for all the layers.

A number of approaches have considered in the open literature to describe with reasonable accuracy the shear behavior of constrained-layer damping treatments. However, the assumptions adopted herein are often used to model moderately thin sandwich beam and plate structures with reasonable accuracy [43].

The displacements are discretized by using linear shape functions for the in-plane displacements of the base plate and constraining layer, and a cubic shape function for the transverse displacement, by the expression, u(x, y, t) = N(x, y)u_(e)(t), where N(x, y) represents the matrix of the interpolation functions, and ${\mathit{u}}_{\left(e\right)}\left(t\right)={\left[u_1^i{v}_1^i{u}_3^i{v}_3^i{w}^i{\theta }_x^i{\theta }_y^i\right]}^T$ with i = 1 to 4 is vector formed by the mechanical nodal DOFs. According to the theory of elasticity, the strain–displacement relations are formulated, ε(x, y, z, t) = B(x, y, z) u_(e)(t), where the strains for elastic layers, ${\mathit{\epsilon }}_k={\left[{\epsilon }_x^{\left(k\right)}{\epsilon }_y^{\left(k\right)}{\gamma }_{xy}^{\left(k\right)}\right]}^T$ with (k = 1, 3), and for the viscoelastic core, ${\mathit{\epsilon }}_2={\left[{\epsilon }_x^{\left(2\right)}{\epsilon }_y^{\left(2\right)}{\gamma }_{xy}^{\left(2\right)}{\gamma }_{xz}^{\left(2\right)}{\gamma }_{yz}^{\left(2\right)}\right]}^T$, are obtained. Thus, based on the stress states assumed for each layer and the stress–strain relations, the stress responses of the system can be obtained as follows:

${\mathit{{\rm M}}}^{\left(e\right)}={\displaystyle \sum _{k=1}^3{\rho }_k{h}_k{\displaystyle \underset{x=0}{\overset{b}{\int }}{\displaystyle \underset{y=0}{\overset{a}{\int }}{\mathit{N}}^T\left(x,y\right)\mathit{N}\left(x,y\right)\mathit{dydx}}}}$E40.a

${\mathit{{\rm K}}}_e^{\left(e\right)}={\displaystyle \sum _{k=1,3}{\displaystyle \underset{z=0}{\overset{h_k}{\int }}{\displaystyle \underset{x=0}{\overset{b}{\int }}{\displaystyle \underset{y=0}{\overset{a}{\int }}{\mathit{B}}_k^T\left(x,y,z\right){\mathit{C}}_k{\mathit{B}}_k\left(x,y,z\right)\mathit{dydxdz}}}}}$E40.b

${\mathit{{\rm K}}}_2^{\left(e\right)}\left(\omega ,T\right)={\displaystyle \underset{z=0}{\overset{h_2}{\int }}{\displaystyle \underset{x=0}{\overset{b}{\int }}{\displaystyle \underset{y=0}{\overset{a}{\int }}{\mathit{B}}_2^T\left(x,y,z\right){\mathit{C}}_2\left(\omega ,T\right){\mathit{B}}_2\left(x,y,z\right)\mathit{dydxdz}}}}$E40.c

where b and a designate, respectively, the dimensions of the rectangular plate element in directions x and y, and h_k and ρ_k represent the thickness and the mass density of the kth layer, respectively. The stiffness matrices, ${\mathit{K}}_e^{\left(e\right)}={\mathit{K}}_1^{\left(e\right)}+{\mathit{K}}_3^{\left(e\right)}$ and ${\mathit{K}}_v^{\left(e\right)}\left(\omega ,T\right)={\mathit{K}}_2^{\left(e\right)}\left(\omega ,T\right)$, give, respectively, the contributions of the purely elastic and viscoelastic parts of the sandwich structure. Hence, the elementary equations of motion are given as follows:

${\mathit{M}}^{\left(e\right)}{\ddot{\mathit{q}}}_{\left(e\right)}\left(t\right)+\left[{\mathit{K}}_e^{\left(e\right)}+{\mathit{K}}_v^{\left(e\ast \right)}\right]{\mathit{q}}_{\left(e\right)}\left(t\right)={\mathit{f}}^{\left(e\right)}\left(t\right)$E41

where ${\mathit{M}}^{\left(e\right)}\in R^{N_e\times N_e}$ is the mass (symmetric, positive definite) matrix, and ${\mathit{K}}_e^{\left(e\right)}\in R^{N_e\times N_e}$ and ${\mathit{K}}_v^{\left(e\ast \right)}={\mathit{K}}_v^{\left(e\right)}\left(\omega ,T\right)\in R^{N_e\times N_e}$ are the stiffness matrices (symmetric, non-negative definite) corresponding to the purely elastic and viscoelastic substructures, respectively. ${\mathit{q}}_{\left(e\right)}\left(t\right)\in R^{N_e}$ and ${\mathit{f}}^{\left(e\right)}\left(t\right)\in R^{N_e}$ are displacements and load vectors, respectively.

9. FE modeling of discrete viscoelastic dampers

From the practical standpoint, the use of viscoelastic materials in mounts and joints is an interesting alternative [1, 31]. Figure 5a illustrates the two mostly used configurations of viscoelastic mounts with the corresponding geometrical parameters. The placement of those mounts in structures is illustrated in Figure 5b. The mounts can be conveniently represented by springs, meaning that a translation mount produces damping forces while a rotational mount generates damping moments. In the same figure, the translational and rotational stiffness coefficients, K^t(s) and K^r(s), are given.

Figure 5.

Sketches (a) and springs representation (b) of discrete viscoelastic devices.

Designating by p the order of the coordinate following the which the viscoelastic mount works, the inclusion of the viscoelastic effect into the equations of motion can be easily done by using the concept of dyadic structural modifications [44]. Thus, the equations of motion of the structural system with viscoelastic mount can be written as follows:

$\mathbf{M}\ddot{\mathbf{q}}\left(t\right)+{\mathbf{K}}_e\mathbf{q}\left(t\right)+G\left(s\right){\overline{\mathbf{K}}}_v\mathbf{q}\left(t\right)=\mathbf{f}\left(s\right)$E42

where ${\mathbf{K}}_v\left(s\right)={\mathbf{K}}^t\left(s\right){\mathit{I}}_p^T{\mathit{I}}_p$ for a translation mount, and ${\mathbf{K}}_v\left(s\right)={\mathbf{K}}^r\left(s\right){\mathit{I}}_p^T{\mathit{I}}_p$ for a rotation mount, and I_p designates the pth column of the identity matrix of order N.

Hence, the global system of equations of motion can be expressed under the form:

10. Numerical examples

The purpose of this section is to perform numerical examples in order to illustrate the main features and capabilities of the viscoelastic modeling procedures intended to design and performance analysis of the viscoelastic damping treatments presented herein. In addition, experimental investigations with a freely suspended rectangular plate were performed, where frequency-response functions (FRFs) and modal analysis have been performed to demonstrate the accuracy of the viscoelastic models and to confirm the effectiveness of the viscoelastic materials applied in the context of vibration attenuation.

10.1. Curve fitting of the viscoelastic model parameters

In the simulations that follow the viscoelastic characteristics of commercially available ISD112 manufactured by 3M [33] have been used. The material data provided by the manufacturer, in terms of storage and loss moduli, at 25°C in the frequency band [8–8000 Hz], have been used to identify the parameters for each viscoelastic model. Eq. (3) was used to form the objective function, which was minimized with respect to the unknown set of model parameters. Such objective function is symbolically defined as follows:

$F={\displaystyle \sum _{i=1}^{N_f}\frac{\left|Re\left[G_{\mathit{model}}\left({\omega }_i\right)\right]-Re\left[G_{\mathit{data}}\left({\omega }_i\right)\right]\right|}{Re\left[G_{\mathit{data}}\left({\omega }_i\right)\right]}}+\frac{\left|Im\left[G_{\mathit{model}}\left({\omega }_i\right)\right]-Im\left[G_{\mathit{data}}\left({\omega }_i\right)\right]\right|}{Im\left[G_{\mathit{data}}\left({\omega }_i\right)\right]}$E43

Optimization was carried out by using genetic algorithms [45], with populations of 800 individuals, allowing for 200 generations and using side constraints. For illustration, Figure 6 shows the storage modulus, loss modulus, and loss factor functions reconstructed from the identified parameters only for the GHM model with five mini-oscillators, superimposed to the experimental counterparts. As can be seen, good quality of the curve fitting could be achieved. The same quality could be obtained for the FDM model and the ADF model with five anelastic fields. Negligible improvement was obtained by increasing the order of those models. The values of the parameters obtained for the models are defined in Table 1.

Figure 6.

Curve fitting of modulus functions for 3M ISD112 according to the GHM model.

10.2. Model validation

To verify the model summarized in Section 8, experimental tests were performed on a freely suspended plate made of aluminum with a constraining damping layer made of a thin ISD112 viscoelastic material and an outer thin aluminum sheet. The experiments consisted in obtaining a set of 20 FRFs corresponding to point I, indicated in Figure 7. Only the average FRF is shown here. The number of elements used to generate the model is shown in the same figure, formed by 378 elastic DOFs, and the anelastic displacements are computed according to the mini-oscillator terms defined in Table 1.

Table 2 provides the physical and geometrical properties used to generate the FE model.

GHM				ADF				FDM
Gr (MPa)	αi	ωi (rad/s)	ζi	Gr (MPa)	Δi	Ωi (rad/s)	Ci	$\begin{array}{l}{G}_0=0.428\left(\mathrm{M}\mathrm{P}\mathrm{a}\right)\\ G_1=0.0088\left(\mathrm{M}\mathrm{P}\mathrm{a}\right)\\ \alpha =0.67\\ \beta =0.41\end{array}$
0.4623	0.26	991.33	4.575	0.4680	0.205	103.48	460.7
	0.95	6986.15	4.27		0.682	103.48	138.2
	2.04	103,437.5	2.702		1.942	638.33	48.5
Gu(MPa)	3.69	22,950.1	1.923	Gu(MPa)	7.062	3054.43	13.4
28.49	53.7	266,466.6	1.299	44.08	83.37	17,583.2	1.13

Table 1.

Identified parameters for the GHM, FDM, and ADF models.

Figure 7.

Illustration of the FE model for the plate with partial viscoelastic treatment.

Base plate	Viscoelastic core	Constrained layer
A = 20 × 10^−2 m; B = 25 × 10^−2 m	C = 2 × 10^−2 m	C = 2 × 10^−2 m
hp = 5 × 10^−4 m	hv = 20 × 10^−5 m	hc = 5 × 10^−4 m
E = 70.3 × 10^9 N/m	ρ = 1099.5 kg/m^3	E = 70.3 × 10^9 N/m
ρ = 2750 kg/m^3	ν = 0.5	ρ = 2750 kg/m^3

Table 2.

Physical and geometrical characteristics of the plate FE model.

Figure 8 shows the amplitudes of the average FRFs calculated from the experiments, compared to the numerically acquired counterparts. It can be seen the efficiency of the surface damping treatment in mitigating the amplitudes of vibrations in the frequency band of interest. Also, it can be noted the accuracy of the model in predicting the dynamic response of the viscoelastic system.

Figure 8.

FE and experimental FRFs of the system with and without treatment.

Mode	1		2		3
Modal parameters	ω1 (Hz)	ζ1	ω2 (Hz)	ζ2	ω3 (Hz)	ζ3
Experimental	184.37	8.1 × 10^−3	242.47	1.9 × 10^−3	429.38	18.2 × 10^−3
FE prediction	184.38	7.5 × 10^−3	238.75	1.8 × 10^−3	434.4	14.4 × 10^−3
Deviations (%)	0.00	8.00	1.60	5.50	1.30	26.40

Table 3.

Experimental natural frequencies and modal damping factors of the plate.

Table 3 compares numerical and experimental frequencies and damping factors obtained by applying the half-power bandwidth method [1]. It can be noted that the two sets are reasonably close to each other. However, the differences observed are mostly due to the identification procedure of the mini-oscillator parameters from the experimental data for the ISD112 material; the theory adopted in the FE model such as the perfectly bounded conditions; the variations on the temperature during the tests; and the boundary conditions.

10.3. Two-dimensional truss with a translational viscoelastic mount

Figure 9 shows the two-dimensional truss FE model in which a translational mount is applied on node 7 with the direction indicated on the same figure.

Figure 9.

Two-dimensional truss with a translational viscoelastic mount.

Figure 10.

Natural frequencies and damping ratios (a) and FRFs (b) for the truss.

Using the theory presented in Section 9 combined with the FDM model, the complex eigenvalue problem was performed to obtain the natural frequencies and damping ratios. The results corresponding to the five vibration modes in the frequency band [80–800 Hz] are presented in Figure 10a. In Figure 10b, the amplitudes of the FRFs of the systems with and without viscoelastic damper are compared. Again, it is possible to evaluate the influence of the damping on the response amplitudes and the influence of the frequency on damping and stiffness of the structure. The FRFs are related to the vertical displacement of node 7 indicated in Figure 9.

10.4. Internally balanced method

Figure 12 shows the results for a beam-like structure partially treated with constrained-layer damping, as illustrated in Figure 11, in terms of the controllability and observability grammians in the balanced realization. It can be noted that W_c and W₀ must be equal and diagonal, as predicted by the theory of internally balanced reduction method detailed in Section 6. Also, since the retained states must be composed by the states with major controllability and observability indices, the first five modes will be considered in the model reduction system. The time response to an impulse excitation applied at point P (indicated in Figure 11), and the amplitudes of the FRFs of the beam before and after reduction are depicted in Figure 13. It can be clearly perceived the efficiency of the internally balanced method in predicting both time and frequency responses of the viscoelastic system by the appropriate choice of the major controllability and observability indices of the states.

Figure 11.

Illustration of the beam partially treated with constraining viscoelastic layer.

Figure 12.

W_c and W₀ versus internally balanced states for the viscoelastic beam.

Figure 13.

Time and FRFs for the full and reduced systems—internally balanced method.

Figure 14.

Time and FRFs for the full and reduced systems—constant enriched basis.

10.5. Modal reduction method

Figure 14 shows the time and frequency domain responses of the reduced beam system by using the enriched modal reduction method compared with the full system. In this case, it has been considered only the first five modes of the associated conservative viscoelastic system, ϕ₀, enriched with the static residues associated with the external loads, R, and the viscoelastic damping forces, ${\mathit{R}}_v^{\mathit{0}}$. As can be seen, both impulse responses and FRFs appear as expected when compared with the time and frequency responses obtained by the internally balanced method, leading to conclude that the reduction method by applying the constant enriched basis (35) is also a viable method to reduce viscoelastic systems.

11. The self-heating phenomenon

The good damping performance and inherent stability of viscoelastic materials in relatively broad frequency bands, besides cost-effectiveness, offers many possibilities for practical engineering applications. However, some drawbacks must be dealt with, such as ageing and chemical instability in the presence of some substances, the mass added and the fact that in most traditional design procedures of viscoelastic dampers subjected to cyclic loadings, uniform and constant temperature is generally assumed and does not take into account the self-heating phenomenon. Also, for viscoelastic dampers subjected to dynamic loadings superimposed on static preloads, especially when good isolation characteristics are required at high frequencies, traditional design guidelines can lead to poor designs or even severe failures, since it is observed a rapidly increasing rate of temperature change and an accompanying stiffness reduction.

The self-heating can cause temperature increases in viscoelastic materials, affecting significantly their damping capacity [26–28]. Thus, in applications in which the viscoelastic materials are subjected to cyclic loadings superimposed on static preloads, such as engine mounts and tall buildings, the interest to obtain high isolation characteristics becomes essential, since the vibration amplitudes are directly related to fatigue and, consequently, to structural integrity [4, 29, 35]. Moreover, depending on the magnitude of the applied loadings, the vibration energy of the viscoelastic material is converted to heat at a rate faster than the heat is conducted away, leading to a rapidly increasing rate of local temperature change known as thermal runaway phenomenon [26]. Thus, it is expected that it can have a strong influence on the stiffness and damping properties of viscoelastic materials, leading to unexpected damping performance or even severe failures of viscoelastic damping devices.

Figure 15 shows the experimental results obtained for a viscoelastic damper subjected to a vertical cyclic loadings during 3396 s, superimposed on different values of static displacement applied to the specimen by the screws shown in the same figure.

Figure 15.

Time evolution of the temperature inside the viscoelastic material and the experimental setup.

Figure 16.

Temperature contours for one half of the damper at t = 100 s for δ = 250 N.

One can conclude that as the static preload increases, the self-heating becomes more pronounced. As a result, an increasing in the temperature values of the viscoelastic material is observed, leading to a significantly reduction of its damping capacity or even its complete failure in practical engineering applications. Moreover, it is not possible to identify a progressive stabilization of the temperatures in the loading phase, indicating the occurrence of the so-called thermal runaway phase [29].

Figure 16 enables to conclude that the assumption of assuming a constant and uniform temperature distribution for viscoelastic materials subjected to cyclic loadings is not correct, since the temperatures are not constant and vary from one point to another.

12. Concluding remarks

A comprehensive review of the modeling strategies of engineering structures incorporating viscoelastic materials has been showed. The FE modeling procedure of two-dimensional sandwich plates treated with viscoelastic materials as a passive constrained-layer damping and a modeling strategy of discrete viscoelastic damping devices including translational and rotational mounts have been also implemented. As can be noted, the modeling of viscoelastic materials was conceived so as to encompass different designs, regarding the type of treatment applied as surface or discrete viscoelastic vibration dampers. The GHM, ADF, and FDM models were used to include the frequency- and temperature-dependent viscoelastic behavior into FE matrices, in spite of the significant increase in the order of the system’s augmented matrices, entailed by the inclusion of internal variables especially for the GHM and ADF models. Moreover, the separation of the material modulus function of each viscoelastic model into real and imaginary parts to enable the identification of the material modulus parameters from experimental data has also been addressed and illustrated for the ISD112 viscoelastic material as detailed in the examples.

The ongoing work aims at developing a user-friendly computer code incorporating various modeling tools available to date to be used for the design, performance analysis, and optimization of different types of viscoelastic vibration dampers taking into account the self-heating phenomenon, as can be available in numerical examples. Also, the implementation of efficient numerical procedures as model reduction methods for the resolution of the equations of motion for modal and frequency-domain analyses of more complex engineering systems incorporating viscoelastic materials was addressed.

In general, the numerical simulations presented enabled to illustrate the application of the modeling procedure as a tool to evaluate the damping effectiveness in terms of eigenvalue and frequency response analysis. Based on the obtained results, one can conclude about the convenience of using more elaborate viscoelastic models in combination with FE models of complex medium- to large-scale structural systems.

Currently, the modeling procedure is being extended to include other types of structural elements, such as three-dimensional beams, plates, and shells. Also, the implementation of efficient numerical and experimental procedures of the self-heating phenomenon and the thermal runaway phase in viscoelastic materials is a topic under investigation.

Acknowledgments

The authors are grateful to the state agencies FAPEMIG and FAPEG for research funding of their research projects. A.M.G. de Lima is thankful to CNPq and CAPES for the financial support to their research activities and for the grant of doctorate and postdoctorate scholarships.